text:
# Disjoint sets

############################################################
#
#   A forest of disjoint sets of integers
#
#   Since each element is an integer, we can use arrays
#   instead of dictionary (for efficiency)
#
#   Disjoint sets over other key types can be implemented
#   based on an IntDisjointSets through a map from the key
#   to an integer index
#
############################################################

type IntDisjointSets
    parents::Vector{Int}
    ranks::Vector{Int}
    ngroups::Int

    # creates a disjoint set comprised of n singletons
    IntDisjointSets(n::Integer) = new([1:n], zeros(Int, n), n)
end

length(s::IntDisjointSets) = length(s.parents)
num_groups(s::IntDisjointSets) = s.ngroups

# find the root element of the subset that contains x
# path compression is implemented here
#
function find_root(s::IntDisjointSets, x::Integer)
    @inbounds p::Int = s.parents[x]
    @inbounds if s.parents[p] != p
        s.parents[x] = p = find_root(s, p)
    end
    p
end

in_same_set(s::IntDisjointSets, x::Integer, y::Integer) = find_root(s, x) == find_root(s, y)

# merge the subset containing x and that containing y into one
#
function union!(s::IntDisjointSets, x::Integer, y::Integer)
    xroot = find_root(s, x)
    yroot = find_root(s, y)
    if xroot != yroot
        rks::Vector{Int} = s.ranks
        @inbounds xrank::Int = rks[xroot]
        @inbounds yrank::Int = rks[yroot]

        if xrank < yrank
            @inbounds s.parents[xroot] = yroot
        else
            @inbounds s.parents[yroot] = xroot
            if xrank == yrank
                s.ranks[xroot] += 1
            end
        end
        s.ngroups -= 1
    end
end

############################################################
#
#  A forest of disjoint sets of arbitrary value type T
#
#  It is a wrapper of IntDisjointSets, which uses a
#  dictionary to map the input value to an internal index
#
############################################################

type DisjointSets{T}
    intmap::Dict{T,Int}
    internal::IntDisjointSets

    function DisjointSets(xs)    # xs must be iterable
        imap = Dict{T,Int}()
        n = length(xs)
        sizehint(imap, n)
        id = 0
        for x in xs
            imap[x] = (id += 1)
        end
        new(imap, IntDisjointSets(n))
    end
end

length(s::DisjointSets) = length(s.internal)
num_groups(s::DisjointSets) = num_groups(s.internal)

find_root{T}(s::DisjointSets{T}, x::T) = find_root(s.internal, s.intmap[x])

in_same_set{T}(s::DisjointSets{T}, x::T, y::T) = in_same_set(s.internal, s.intmap[x], s.intmap[y])

function union!{T}(s::DisjointSets{T}, x::T, y::T)
    union!(s.internal, s.intmap[x], s.intmap[y])
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _posang{T<:AbstractFloat}(units::Integer, ra1::T, dec1::T, ra2::T, dec2::T)
    # Convert all quantities to radians.
    if units == 0
        # All radians
        ra1_rad  = ra1
        ra2_rad  = ra2
        dec1_rad = dec1
        dec2_rad = dec2
    elseif units == 1
        # Right ascensions are in hours, declinations in degrees.
        ra1_rad  = ra1*pi/12.0
        ra2_rad  = ra2*pi/12.0
        dec1_rad = deg2rad(dec1)
        dec2_rad = deg2rad(dec2)
    elseif units == 2
        # Right ascensions and declinations are in degrees.
        ra1_rad  = deg2rad(ra1)
        ra2_rad  = deg2rad(ra2)
        dec1_rad = deg2rad(dec1)
        dec2_rad = deg2rad(dec2)
    else
        # In any other case throw an error.
        error("units must be 0 (radians), 1 (hours, degrees) or 2 (degrees)")
    end
    radif = ra2_rad - ra1_rad
    angle = atan2(sin(radif), cos(dec1_rad)*tan(dec2_rad) -
                  sin(dec1_rad)*cos(radif))
    if units == 0
        return angle
    else
        return rad2deg(angle)
    end
end

"""
    posang(units, ra1, dec1, ra2, dec2) -> angular_distance

### Purpose ###

Compute rigorous position angle of point 2 relative to point 1.

### Explanation ###

Computes the rigorous position angle of point 2 (with given right ascension and
declination) using point 1 (with given right ascension and declination) as the
center.

### Arguments ###

* `units`: integer, can be either 0, or 1, or 2.  Describes units of inputs and
 output:
    * 0: everything (input right ascensions and declinations, and output
      distance) is radians
    * 1: right ascensions are in decimal hours, declinations in decimal degrees,
      output distance in degrees
    * 2: right ascensions and declinations are in degrees, output distance in
      degrees
* `ra1`:  right ascension or longitude of point 1
* `dec1`: declination or latitude of point 1
* `ra2`: right ascension or longitude of point 2
* `dec2`: declination or latitude of point 2

Both `ra1` and `dec1`, and `ra2` and `dec2` can be given as 2-tuples `(ra1,
dec1)` and `(ra2, dec2)`.

### Output ###

Angle of the great circle containing `[ra2, dec2]` from the meridian containing
`[ra1, dec1]`, in the sense north through east rotating about `[ra1, dec1]`.
See `units` argument above for units.

### Method ###

The "four-parts formula" from spherical trigonometry (p. 12 of Smart's Spherical
Astronomy or p. 12 of Green' Spherical Astronomy).

### Example ###

Mizar has coordinates (ra, dec) = (13h 23m 55.5s, +54° 55' 31'').  Its
companion, Alcor, has coordinates (ra, dec) = (13h 25m 13.5s, +54° 59' 17'').
Find the position angle of Alcor with respect to Mizar.

``` julia
posang(1, ten(13, 25, 13.5), ten(54, 59, 17), ten(13, 23, 55.5), ten(54, 55, 31))
# => -108.46011246802047
```

### Notes ###

* If `ra1`, `dec1` are scalars, and `ra2`, `dec2` are vectors, then the output
 is a vector giving the distance of each element of `ra2`, `dec2` to `ra1`,
 `dec1`.  Similarly, if `ra1`,`de1` are vectors, and `ra2`,` dec2` are scalars,
 then the output is a vector giving the distance of each element of `ra1`,
 `dec1` to `ra2`, `dec2`.  If both `ra1`, `dec1` and `ra2`, `dec2` are vectors
 then the output is a vector giving the distance of each element of `ra1`,
 `dec1` to the corresponding element of `ra2`, `dec2`.
* The function `sphdist` provides an alternate method of computing a spherical
 distance.
* Note that `posang` is not commutative: the position angle between A and B is
  \$\\theta\$, then the position angle between B and A is \$180 + \\theta\$.

Code of this function is based on IDL Astronomy User's Library.
"""
posang(units::Integer, ra1::Real, dec1::Real, ra2::Real, dec2::Real) =
    _posang(units, promote(float(ra1), float(dec1), float(ra2), float(dec2))...)

function posang{R1<:Real, D1<:Real}(units::Integer,
                                   ra1::AbstractArray{R1},
                                   dec1::AbstractArray{D1},
                                   ra2::Real,
                                   dec2::Real)
    @assert length(ra1) == length(dec1)
    dist = similar(ra1, typeof(float(one(R1))))
    for i in eachindex(ra1)
        dist[i] = posang(units, ra1[i], dec1[i], ra2, dec2)
    end
    return dist
end

function posang{R2<:Real, D2<:Real}(units::Integer,
                                   ra1::Real,
                                   dec1::Real,
                                   ra2::AbstractArray{R2},
                                   dec2::AbstractArray{D2})
    @assert length(ra2) == length(dec2)
    dist = similar(ra2, typeof(float(one(R2))))
    for i in eachindex(ra2)
        dist[i] = posang(units, ra1, dec1, ra2[i], dec2[i])
    end
    return dist
end

function posang{R1<:Real, D1<:Real, R2<:Real, D2<:Real}(units::Integer,
                                                       ra1::AbstractArray{R1},
                                                       dec1::AbstractArray{D1},
                                                       ra2::AbstractArray{R2},
                                                       dec2::AbstractArray{D2})
    @assert length(ra1) == length(dec1) == length(ra2) == length(dec2)
    dist = similar(ra1, typeof(float(one(R1))))
    for i in eachindex(ra1)
        dist[i] = posang(units, ra1[i], dec1[i], ra2[i], dec2[i])
    end
    return dist
end

### Tuples input
posang(units::Integer, radec1::Tuple{Real, Real}, ra2::Real, dec2::Real) =
    posang(units, radec1..., ra2, dec2)

posang(units::Integer, ra1::Real, dec1::Real, radec2::Tuple{Real, Real}) =
    posang(units, ra1, dec1, radec2...)

posang(units::Integer, radec1::Tuple{Real, Real}, radec2::Tuple{Real, Real}) =
    posang(units, radec1..., radec2...)

ext: jl


text:
```@meta
DocTestSetup = quote
    using LightPropagation
end
```

##### Diffusion approximation for semi-infinite media (reflectance) #####

#=

@article{contini1997photon,
  title={Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory},
  author={Contini, Daniele and Martelli, Fabrizio and Zaccanti, Giovanni},
  journal={Applied optics},
  volume={36},
  number={19},
  pages={4587--4599},
  year={1997},
  publisher={Optical Society of America}
}

The reflectance for the semi-infinite medium is calculated from Equation 36 from Contini 1997.
The semi-infinite medium is described by retaining the first of the dipole sources (m=0) from Eqn. 36.
=#

"""
    TPSF_DA_semiinf_refl(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)

Compute the time-domain reflectance from a semi-infinite medium.

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction

# Examples
```jldoctest
julia> TPSF_DA_semiinf_refl(0:1:5, [0.1,10.0], 1.0, 1.0, 1.0)
6-element Array{Float64,1}:
 0.0
 0.0001440103022493725
 1.446739954231315e-6
 2.7354735244571076e-8
 6.794070985483474e-10
 1.9657536202689858e-11
```
"""
function TPSF_DA_semiinf_refl(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
  	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	if n > 1.0
		  A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		  5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		  A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		  A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

    Threads.@threads for n in eachindex(t)

        Rt[n] = -(ρ^2/(4D*ν*t[n]))
        Rt[n] = Rt[n] - μa*ν*t[n]
        Rt[n] = -exp(Rt[n])/(2*(4π*D*ν)^(3/2)*sqrt(t[n]*t[n]*t[n]*t[n]*t[n]))

        Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

        if isnan(Rt[n])
            Rt[n] = 0
        end

    end

     return Rt

end

##### Diffusion approximation for Slab Geometry (reflectance) #####

#=
The reflectance for the slab is calculated from Equation 36 from Contini 1997 for the first 11 dipoles (-10:1:10).
The number of sources considered in the summation can be changed with xs. Larger values of m will be needed for large values of SDS and t.
=#

"""
    TPSF_DA_slab_refl(t, β::Array{Float64,1}, ρ::Float64,ndet::Float64, nmed::Float64, s::Float64)

Compute the time-domain reflectance from a slab geometry (x,y->inf, z-> finite).

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `s::Float64`: the thickness (z-depth) of the slab (cm)

# Examples
```jldoctest
julia> TPSF_DA_slab_refl(0:1:5, [0.1,10.0], 1.0, 1.0, 1.0, 2.0)
6-element Array{Float64,1}:
 0.0
 0.00011885510267563147
 3.8254510325745536e-7
 1.5189689097582458e-9
 6.645341571314721e-12
 3.075337855389727e-14

```
"""
function TPSF_DA_slab_refl(t, β::Array{Float64,1}, ρ::Float64,ndet::Float64, nmed::Float64, s::Float64)
	#s is slab thickness

	  n::Float64 = nmed/ndet
      μa::Float64 = β[1]
      μsp::Float64 = β[2]
      D::Float64 = 1/3μsp
      ν::Float64 = 29.9792345/nmed
      xs::UnitRange{Int64} = -10:10

      Rt1 = Array{Float64}(undef, length(t))
	  Rt2 = zeros(Float64, length(t))

		@inbounds begin

			if n == 1.0
				A= 1.0
			elseif n > 1.0
				A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 		5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
			else
				A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
			end

			zs::Float64 = 1/μsp
			ze::Float64 = 2A*D

		    Rt1 = @. -exp(-(ρ^2/(4D*ν*t)) - μa*ν*t)/(2*(4π*D*ν)^(3/2)*t^(5/2))

				Threads.@threads for n in eachindex(t)
					for m in xs

						z3m = -2m.*s .- 4m.*ze .- zs
						z4m = -2m.*s .- (4m .- 2).*ze .+ zs

			 			Rt2[n] += z3m*exp(-(z3m^2 / (4D*ν*t[n]))) - z4m*exp(-(z4m^2 / (4D*ν*t[n])))
					end

				end
			end

	return replace!(Rt1.*Rt2, NaN => 0)
end
##### Diffusion approximation for a turbid parallelepiped (reflectance)  #####

#=
@article{kienle2005light,
  title={Light diffusion through a turbid parallelepiped},
  author={Kienle, Alwin},
  journal={JOSA A},
  volume={22},
  number={9},
  pages={1883--1888},
  year={2005},
  publisher={Optical Society of America}
}

The fluence is given by equation 6 in Kienle 2005. The reflectance as calculated here uses Fick's law where Rt(rho, t) = D∂/∂zϕ(rho, z=0, t)
=#

"""
TPSF_DA_paralpip_refl(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

Compute the time-domain reflectance from a parallelepiped [lx, ly, lz].

# Arguments
- `t`: the time vector (ns).
- `β::Array{Float64,1}`: the optical properties [μa, μs'] (cm⁻¹)
- `ρ::Float64`: the source detector separation (cm⁻¹)
- `ndet::Float64`: the boundary's index of refraction (air or detector)
- `nmed::Float64`: the sample medium's index of refraction
- `rd::Array{Float64,1}`: the location of the detector [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs::Array{Float64,1}`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L::Array{Float64,1}`: the dimenensions [lx, ly, lz] of the parallelepied

# Examples
```jldoctest
julia> TPSF_DA_paralpip_refl(0:1:5, [0.1,10.0], 1.0, 1.0, [2.5,5.0], [5.0,5.0], [10.0,10.0,10.0])
6-element Array{Float64,1}:
 0.0
 3.872453933271761e-5
 7.490933986486167e-7
 1.741880709155373e-8
 4.687789101530684e-10
 1.3825359487332909e-11

```
"""
function TPSF_DA_paralpip_refl(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

	n::Float64 = nmed/ndet
	μa::Float64 = β[1]
	μsp::Float64 = β[2]
	D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	x::Float64 = rd[1]
	y::Float64 = rd[2]
	xu::Float64 = rs[1]
	yu::Float64 = rs[2]
	lx::Float64 = L[1]
	ly::Float64 = L[2]
	lz::Float64 = L[3]

	Rt2 = zeros(Float64, length(t))
	Rt3 = zeros(Float64, length(t))
	Rt4 = zeros(Float64, length(t))

		if n == 1.0
			A= 1.0
		elseif n > 1.0
			A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
			 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
		else
			A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
		end

		zo::Float64 = 1/μsp
		zb::Float64 = 2A*D

	  Rt1 = @. exp(-μa*ν*t)/(2*(4π*D*ν)^(3/2)*t^(5/2))

	  Threads.@threads for n in eachindex(t)
		  for m in xs
			    x1l = 2m.*lx .+ 4m.*zb .+ xu
				x2l = 2m.*lx + (4m-2).*zb .- xu
				y1m = 2m.*ly .+ 4m.*zb .+ yu
				y2m = 2m.*ly .+ (4m.-2).*zb .- yu
				z1n = 2m.*lz .+ 4m.*zb .+ zo
				z2n = 2m.*lz .+ (4m.-2).*zb .- zo

				Rt2[n] += exp(-(x-x1l)^2 / (4D*ν*t[n])) - exp(-(x-x2l)^2 / (4D*ν*t[n]))

				Rt3[n] += exp(-(y-y1m)^2 / (4D*ν*t[n])) - exp(-(y-y2m)^2 / (4D*ν*t[n]))

				Rt4[n] += z1n*exp(-z1n^2 /(4D*ν*t[n])) - z2n'*exp(-z2n^2 /(4D*ν*t[n]))

		  end
	  end
	return replace!(Rt1.*Rt2.*Rt3.*Rt4, NaN => 0)
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
export ADCG

function ADCG(sim :: ForwardModel, lossFn :: Loss, y :: Vector{Float64}, tau :: Float64;
  callback :: Function = (old_thetas,thetas, weights,output,old_obj_val) -> false,
  max_iters :: Int64 = 50,
  min_optimality_gap :: Float64 = 1E-5,
  max_cd_iters :: Int64 = 200,
  fully_corrective :: Bool = false)
  #assert(tau > 0.0)
  if (tau<0.0)
    println("The maximal mass of the total measure (tau) cannot be less than zero")
  end
  bound = -Inf
  thetas = zeros(0,0) #hack
  weights = zeros(0)
  #cache the forward model applied to the current measure.
  output = zeros(length(y))
  for iter = 1:max_iters
    #compute the current residual
    residual = output - y
    #evalute the objective value and gradient of the loss
    objective_value, grad = loss(lossFn, residual)
    #compute the next parameter value to add to the support
    theta,score = lmo(sim,grad)
    #println("The objective value for the current residual is: ", string(objective_value))
    #println("The score for the current residual is: ", string(score))
    #println("The parameter locations are: ", string(theta))
    #score is - |<\psi(theta), gradient>|
    #update the lower bound on the optimal value
    bound = max(bound, objective_value+score*tau-dot(output,grad))
    #println("The lower bound on the optimal value is: ", string(bound))
    #check if the bound is met.
    #if(objective_value > bound + min_optimality_gap || score >= 0.0)
      #print('h', thetas)
    #  return thetas,weights
    #end
    #update the support
    old_thetas = thetas
    thetas = iter == 1 ? reshape(theta, length(theta),1) : [thetas theta]
    #run local optimization over the support.
    old_weights = copy(weights)
    #println("The updated Theta array is: ", string(thetas))
    #println(string(thetas))
    thetas,weights = localUpdate(sim,lossFn,thetas,y,tau,max_cd_iters)
    output = phi(sim, thetas, weights)
    if callback(old_thetas, thetas,weights, output, objective_value)
      return old_thetas, old_weights
    end
  end
#  print(thetas)
  println("Hit max iters in frank-wolfe!")
  return thetas, weights
end

#Default implementation using coordinate descent
#Feel free to override
function localUpdate(sim :: ForwardModel,lossFn :: Loss,
    thetas :: Array{Float64, 2}, y :: Array{Float64,1}, tau :: Float64, max_iters)
  for cd_iter = 1:max_iters
    weights = solveFiniteDimProblem(sim, lossFn, thetas, y, tau)
    #remove points with zero weight
    if any(weights.==0.0)
      println("Removing ",sum(weights.==0.0), " zero-weight points.")
      thetas = thetas[:,weights.!= 0.0]
      weights = weights[weights.!= 0.0]
    end
    #local minimization over the support
    new_thetas = localDescent(sim, lossFn, thetas,weights, y)
    #break if termination condition is met
    if length(thetas) == length(new_thetas) && maximum(abs.(vec(thetas)-vec(new_thetas))) <= 1E-7
        break
    end
    thetas = new_thetas
  end
  #final prune
  weights = solveFiniteDimProblem(sim, lossFn, thetas, y, tau)
  if any(weights.==0.0)
    println("Removing ",sum(weights.==0.0), " zero-weight points.")
    thetas = thetas[:,weights.!= 0.0]
    weights = weights[weights.!= 0.0]
  end
  return thetas, weights
end


text:
module Binning
using Flux

"""
    SmoothStep(x, bᵢ, bⱼ; α=5)

Construct a bin function with edges bᵢ and bⱼ with steepness α.

# Examples
```julia-repl
julia> SmoothStep(5, 0.0, 10.0, α=5.0)
1.0
```
"""
function SmoothStep(x, bᵢ, bⱼ; α=5)
    # assume i < j
    return 0.5(tanh(α*(x-bᵢ))-tanh(α*(x-bⱼ))) # max val 1 min value 0 sharpness controlled by α
end

# binning layer defined by following example from:
# https://fluxml.ai/Flux.jl/stable/models/advanced/

struct DomainBinner{S<:AbstractArray, T<:AbstractArray, N<:Number, M<:Number}
    x::S # input values
    b::T # vector of bin edges
    b₁::N # start bin
    bₑ::M # end bin
end

# define layer constructor
"""
    DomainBinner(x, N)

Construct a DomainBinner layer to be used with a Flux model.
`x` are the current domain values and `N` is the desired number of bins.

# Examples
```julia-repl
julia> DomainBinner(0:1:1000, 100)
DomainBinner{StepRange{Int64, Int64}, Vector{Float64}, Int64, Int64}(0:1:1000, [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0  …  900.0, 910.0, 920.0, 930.0, 940.0, 950.0, 960.0, 970.0, 980.0, 990.0], 0, 1000)
```
"""
function DomainBinner(x::AbstractArray, N::Integer)
    x_min = minimum(x)
    x_max = maximum(x)

    # may wan't to specify different bin initialization according to some distribution
    Δx = (x_max-x_min)/N
    b1 = x_min
    be = x_max
    b = Array(x_min+Δx:Δx:x_max-Δx)

    return DomainBinner(x, b, b1, be)
end

# set the trainable parameters to be bin positions within the boundaries
#Flux.trainable(L::DomainBinner) = (L.b[2:end-1],)
Flux.trainable(L::DomainBinner) = (L.b,)

# define how to call layer as a function
# Note
# function (BL::DomainBinner)(f::AbstractArray)
#     # collect bins and sort in ascending order
#     b = sort(vcat(BL.b₁,BL.b, BL.bₑ))
#     x = BL.x

#     # we need to treat multidimensional data carefully since this isn't a matrix operation
#     # Think of rows as instances of data so that shape is (length data, num features)
#     if ndims(f) == 1
#         fout = [sum([f[j]*SmoothStep(x[j], b[i-1], b[i]) for j ∈ 1:length(x)]) for i ∈ 2:length(b)]
#     else
#         # fout = [sum([f_row[j]*SmoothStep(x[j], b[i-1], b[i]) for j ∈ 1:length(x)]) for   f_row in eachrow(f),  i ∈ 2:length(b)]

#         fout = [sum([f_col[j]*SmoothStep(x[j], b[i-1], b[i]) for j ∈ 1:length(x)]) for i ∈ 2:length(b), f_col in eachcol(f)]
#         # fout = fout'

#     end

#     return fout
# end

function (BL::DomainBinner)(f::AbstractVector)
    # collect bins and sort in ascending order
    b = sort(vcat(BL.b₁,BL.b, BL.bₑ))
    x = BL.x

    fout = [sum([f[j]*SmoothStep(x[j], b[i-1], b[i]) for j ∈ 1:length(x)]) for i ∈ 2:length(b)]
    return fout
end

# function (BL::DomainBinner)(f::AbstractArray)
#     fout = BL.(eachcol(f))
# end

"""
    getBinCenters(BL:DomainBinner)

Return the current bin centers for a `DomainBinner` `BL`

# Examples
```julia-repl
julia> getBinCenters(BL)
100-element Vector{Float64}:
   5.0
  15.0
  25.0
  35.0
  45.0
  55.0
  65.0
  75.0
  85.0
  95.0
   ⋮
 915.0
 925.0
 935.0
 945.0
 955.0
 965.0
 975.0
 985.0
 995.0
```
"""
function getBinCenters(BL::DomainBinner)
    b = sort(vcat(BL.b₁,BL.b, BL.bₑ))
    [(b[i]+b[i-1])/2 for i ∈ 2:length(b)]
end

# not sure exactly what this is doing but it's in the docs
Flux.@functor DomainBinner

export SmoothStep
export DomainBinner
export getBinCenters

end


text:
#Initialize functions

export mymotor
function mymotor()::Array{XIMCMotor,1}
	m = [
           XIMCMotor("gelab-serial03", 2001),
           XIMCMotor("gelab-serial03", 2011)
       ]
	Initialize(m)
	return m
end

function _Initialize(motor::XIMCMotor)
	fbs1 = motor[Feedback, Settings]
	sleep(1)
	fbs1.feedback_flags = 0x01
	motor[Feedback, Settings] = fbs1
	sleep(1)
	# @assert fbs1 == motor[1][Feedback, Settings]

	hms1 = motor[Home, Settings]
	sleep(1)
	hms1.fast_home_spd = 1000
	hms1.home_delta = 0	 	 #was 50 before
	hms1.home_flags = 0x0173 #was 0x0071 before
	motor[Home, Settings] = hms1
	# @assert motor[1][Home, Settings] == hms1

	sleep(1)
	mvs1 = motor[Move, Settings]
	sleep(1)
	mvs1.speed = 500  #was 25 before
	mvs1.accel = 500  #was 25 before
	mvs1.decel = 500  #was 25 before
	mvs1.antiplay_speed = 500 # was 5 before
	motor[Move, Settings] = mvs1
	# @assert mvs1 == motor[1][Move, Settings]

	sleep(1)
	ens1 = motor[Engine, Settings]
	sleep(1)
	ens1.engine_flags = 0x91 #was 0x99 before
	motor[Engine, Settings] = ens1
	# @assert ens1 == motor[1][Engine, Settings]

	pws1 = motor[Power, Settings]
	pws1.power_off_delay = 2
	motor[Power, Settings] = pws1
	# @assert pws1 == motor[1][Power, Settings]
end

export Initialize
function Initialize(motor::Array{XIMCMotor,1})
	_Initialize(motor[1])
	_Initialize(motor[2])
    @info "Motors succefully initialized!"
end

#Calibration functions

function _Calibrate(motor::XIMCMotor)
	motor[Move, Absolute] = :home
	#block until it arrives "home"
	currentpos = _Pos(motor)
	sleep(2)
	while abs(_Pos(motor) - currentpos) > 0.001
		currentpos = _Pos(motor)
		sleep(2)
	end
	pos = motor[Position]
	sleep(1)
	pos.pos = 0
	pos.enc_pos = 0
	motor[Position] = pos
end

export Calibrate
function Calibrate(motor::Array{XIMCMotor,1})
    CalibrateX(motor)
    CalibrateY(motor)
end

export CalibrateX
function CalibrateX(motor::Array{XIMCMotor,1})
    @info "Calibrating X motor stage"
    _Calibrate(motor[1])
end

export CalibrateY
function CalibrateY(motor::Array{XIMCMotor,1})
    @info "Calibrating Y motor stage"
    _Calibrate(motor[2])
end

#Movement functions

function _MoveMM(pos, motor::XIMCMotor; block_till_arrival = true)
	if 0 <= pos <= 100
	    motor[Move, Absolute] = round(pos * -400)
	    if block_till_arrival == true
	        i_cnt = 0
	        while abs(_Pos(motor) - pos) > 0.1
	            if i_cnt%30 == 0
	                motor[Move, Absolute] = round(pos * -400)
	            end
	            sleep(1)
	            i_cnt +=1
	        end
	        motor[Move, Absolute] = round(pos * -400)
	        sleep(1)
	        @info("arrived at $(round(_Pos(motor),digits = 2))mm")
	    end
	else
		@info "Position $(pos)mm is not supported. Please enter a value between 0 and 100"
	end

end

export XMoveMM
function XMoveMM(pos, motor::Array{XIMCMotor,1}; block_till_arrival = true)
    @info("Moving to X-motor to $(pos)mm\n")
	_MoveMM(pos, motor[1], block_till_arrival = block_till_arrival)
end

export YMoveMM
function YMoveMM(pos, motor::Array{XIMCMotor,1}; block_till_arrival = true)
    @info("Moving to Y-motor to $(pos)mm\n")
    _MoveMM(pos, motor[2], block_till_arrival = block_till_arrival)
end

#Position functions

function _Pos(motor::XIMCMotor)::Float64
	if motor[Position].pos == 0 return 0.0 end
    return motor[Position].pos/(-400.)
end

export PosX
function PosX(motor::Array{XIMCMotor,1})
	_Pos(motor[1])
end

export PosY
function PosY(motor::Array{XIMCMotor,1})
	_Pos(motor[2])
end

export Pos
function Pos(motor::Array{XIMCMotor,1})
	println("X: \t$(PosX(motor)) mm")
	println("Y: \t$(PosY(motor)) mm")
end


text:
algebra(a::AbsAlgAssIdl) = a.algebra

iszero(a::AbsAlgAssIdl) = (a.iszero == 1)

###############################################################################
#
#  String I/O
#
###############################################################################

function show(io::IO, a::AbsAlgAssIdl)
  print(io, "Ideal of ")
  print(io, algebra(a))
  println(io, " with basis matrix")
  print(io, basis_mat(a, copy = false))
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AbsAlgAssIdl, dict::IdDict)
  b = typeof(a)(algebra(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :algebra
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  return b
end

################################################################################
#
#  Basis (matrices)
#
################################################################################

function assure_has_basis(a::AbsAlgAssIdl)
  if isdefined(a, :basis)
    return nothing
  end

  A = algebra(a)
  M = basis_mat(a, copy = false)
  a.basis = Vector{elem_type(A)}(undef, nrows(M))
  for i = 1:nrows(M)
    a.basis[i] = elem_from_mat_row(A, M, i)
  end
  return nothing
end

function basis(a::AbsAlgAssIdl; copy::Bool = true)
  assure_has_basis(a)
  if copy
    return deepcopy(a.basis)
  else
    return a.basis
  end
end

function basis_mat(a::AbsAlgAssIdl; copy::Bool = true)
  if copy
    return deepcopy(a.basis_mat)
  else
    return a.basis_mat
  end
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

function in(x::T, a::AbsAlgAssIdl{S, T, U}) where {S, T, U}
  A = algebra(a)
  M = matrix(base_ring(A), 1, dim(A), coeffs(x, copy = false))
  return rank(vcat(basis_mat(a, copy = false), M)) == nrows(basis_mat(a, copy = false)) # so far we assume nrows(basis_mat) == rank(basis_mat)
end

################################################################################
#
#  Test right/left
#
################################################################################

function _test_ideal_sidedness(a::AbsAlgAssIdl, side::Symbol)

  A = algebra(a)
  ba = basis(a, copy = false)
  t = A()
  for i = 1:dim(A)
    for j = 1:length(ba)
      if side == :left
        t = mul!(t, A[i], ba[j])
      elseif side == :right
        t = mul!(t, ba[j], A[i])
      else
        error("side must be either :left or :right")
      end
      if !(t in a)
        return false
      end
    end
  end
  return true
end

function _test_ideal_sidedness(a::AlgAssAbsOrdIdl, side::Symbol)

  O = order(a)
  ba = basis(a, copy = false)
  bO = basis(O, copy = false)
  t = O()
  for i = 1:degree(O)
    for j = 1:length(ba)
      if side == :left
        t = mul!(t, bO[i], ba[j])
      elseif side == :right
        t = mul!(t, ba[j], bO[i])
      else
        error("side must be either :left or :right")
      end
      if !(t in a)
        return false
      end
    end
  end
  return true
end

function isright_ideal(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl })
  if a.isright == 1
    return true
  elseif a.isright == 2
    return false
  end

  if _test_ideal_sidedness(a, :right)
    a.isright = 1
    return true
  end

  a.isright = 2
  return false
end

function isleft_ideal(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl })
  if a.isleft == 1
    return true
  elseif a.isleft == 2
    return false
  end

  if _test_ideal_sidedness(a, :left)
    a.isleft = 1
    return true
  end

  a.isleft = 2
  return false
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
  if iszero(a)
    return deepcopy(b)
  elseif iszero(b)
    return deepcopy(a)
  end

  M = vcat(basis_mat(a), basis_mat(b))
  r = rref!(M)
  if r != nrows(M)
    M = sub(M, 1:r, 1:ncols(M))
  end
  return ideal(algebra(a), M, :nothing, true)
end

function *(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
  if iszero(a)
    return deepcopy(a)
  elseif iszero(b)
    return deepcopy(b)
  end

  A = algebra(a)
  ba = basis(a, copy = false)
  bb = basis(b, copy = false)
  M = zero_matrix(base_ring(A), length(ba)*length(bb), dim(A))
  for i = 1:length(ba)
    ii = (i - 1)*length(bb)
    for j = 1:length(bb)
      elem_to_mat_row!(M, ii + j, ba[i]*bb[j])
    end
  end
  return ideal(algebra(a), M, :nothing)
end

^(A::AbsAlgAssIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AbsAlgAssIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

function one(a::AbsAlgAssIdl)
  return ideal(algebra(a), identity_matrix(base_ring(A), dim(A)), :twosided, true)
end

function Base.copy(a::AbsAlgAssIdl)
  return a
end

################################################################################
#
#  Equality
#
################################################################################

function ==(a::AbsAlgAssIdl, b::AbsAlgAssIdl)
  algebra(a) != algebra(b) && return false
  return basis_mat(a, copy = false) == basis_mat(b, copy = false)
end

################################################################################
#
#  Construction
#
################################################################################

function ideal_from_gens(A::AbsAlgAss, b::Vector{T}, side::Symbol = :nothing) where { T <: AbsAlgAssElem }
  if length(b) == 0
    M = zero_matrix(base_ring(A), 0, dim(A))
    return ideal(A, M, side, true)
  end

  @assert parent(b[1]) == A

  M = zero_matrix(base_ring(A), length(b), dim(A))
  for i = 1:length(b)
    elem_to_mat_row!(M, i, b[i])
  end
  return ideal(A, M, side)
end

function ideal(A::AbsAlgAss, x::AbsAlgAssElem)
  t1 = A()
  t2 = A()
  M = zero_matrix(base_ring(A), dim(A)^2, dim(A))
  for i = 1:dim(A)
    t1 = mul!(t1, A[i], x)
    ii = (i - 1)*dim(A)
    for j = 1:dim(A)
      t2 = mul!(t2, t1, A[j])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end

  return ideal(A, M, :twosided)
end

function ideal(A::AbsAlgAss, x::AbsAlgAssElem, action::Symbol)
  M = representation_matrix(x, action)
  a = ideal(A, M)

  if action == :left
    a.isright = 1
  elseif action == :right
    a.isleft = 1
  end

  return a
end

function ideal(A::AbsAlgAss, M::MatElem, side::Symbol = :nothing, M_in_rref::Bool = false)
  @assert base_ring(M) == base_ring(A)
  @assert ncols(M) == dim(A)
  if !M_in_rref
    r, N = rref(M)
    if r == 0
      a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, zero_matrix(base_ring(A), 0, dim(A)))
      a.iszero = 1
      return a
    end
    if r != nrows(N)
      M = sub(N, 1:r, 1:ncols(N))
    else
      M = N
    end
  end
  if M_in_rref && nrows(M) == 0
    a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, M)
    a.iszero = 1
    return a
  end

  a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, M)
  _set_sidedness(a, side)
  a.iszero = 2
  return a
end

# Helper function to set the side-flags
# side can be :right, :left or :twosided
function _set_sidedness(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl }, side::Symbol)
  if side == :right
    a.isleft = 0
    a.isright = 1
  elseif side == :left
    a.isleft = 1
    a.isright = 0
  elseif side == :twosided
    a.isleft = 1
    a.isright = 1
  else
    a.isleft = 0
    a.isright = 0
  end
  return nothing
end

################################################################################
#
#  Quotient rings
#
################################################################################

function quo(A::S, a::AbsAlgAssIdl{S, T, U}) where { S, T, U }
  @assert A == algebra(a)
  K = base_ring(A)

  # First compute the vector space quotient
  Ma = basis_mat(a, copy = false)
  M = hcat(deepcopy(transpose(Ma)), identity_matrix(K, dim(A)))
  r = rref!(M)
  pivot_cols = Vector{Int}()
  j = 1
  for i = 1:ncols(M)
    if !iszero(M[j, i])
      if i > nrows(Ma)
        push!(pivot_cols, i - nrows(Ma))
      end
      j += 1
      if j > nrows(M)
        break
      end
    end
  end

  # We now have the basis (basis of the quotient, basis of the ideal)
  n = dim(A) - nrows(Ma)
  M = vcat(zero_matrix(K, n, dim(A)), Ma)
  oneK = K(1)
  zeroK = K()
  for i = 1:n
    M[i, pivot_cols[i]] = oneK
  end
  iM = inv(M)

  N = sub(M, 1:n, 1:dim(A))
  NN = sub(iM, 1:dim(A), 1:n)

  # Lift a basis of the quotient to A
  quotient_basis = Vector{elem_type(A)}(undef, n)
  b = zero_matrix(K, 1, n)
  for i = 1:n
    b[1, i] = oneK
    bN = b*N
    quotient_basis[i] = A([ bN[1, i] for i = 1:dim(A) ])
    b[1, i] = zeroK
  end

  # Build the multiplication table
  t = A()
  s = zero_matrix(K, 1, dim(A))
  mult_table = Array{elem_type(K), 3}(undef, n, n, n)
  for i = 1:n
    for j = 1:n
      t = mul!(t, quotient_basis[i], quotient_basis[j])
      elem_to_mat_row!(s, 1, t)
      sNN = s*NN
      mult_table[i, j, :] = [ sNN[1, k] for k = 1:n ]
    end
  end

  B = AlgAss(K, mult_table)
  AtoB = hom(A, B, NN, N)
  return B, AtoB
end

# Assumes b \subseteq a
function quo(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where { S, T, U }
  @assert algebra(a) == algebra(b)
  A = algebra(a)
  K = base_ring(A)

  # First compute the vector space quotient
  Ma = basis_mat(a, copy = false)
  Mb = basis_mat(b, copy = false)
  M = hcat(deepcopy(transpose(Mb)), deepcopy(transpose(Ma)))
  r = rref!(M)
  pivot_cols = Vector{Int}()
  j = 1
  for i = 1:ncols(M)
    if !iszero(M[j, i])
      if i > nrows(Mb)
        push!(pivot_cols, i - nrows(Mb))
      end
      j += 1
      if j > nrows(M)
        break
      end
    end
  end

  # Build the basis matrix for the quotient
  M = zero_matrix(K, dim(A), dim(A))
  n = nrows(Ma) - nrows(Mb)
  for i = 1:n
    for j = 1:dim(A)
      M[i, j] = deepcopy(Ma[pivot_cols[i], j])
    end
  end

  mult_table = _build_subalgebra_mult_table!(A, M)
  # M is now in rref

  B = AlgAss(K, mult_table)
  MM = sub(M, 1:dim(B), 1:dim(A))

  AtoB = AbsAlgAssMor{typeof(A), typeof(B), typeof(MM)}(A, B)

  N = transpose(vcat(MM, Mb)) # Another basis matrix for a
  function _image(x::AbsAlgAssElem)
    t, y = can_solve(N, matrix(K, dim(A), 1, coeffs(x, copy = false)))
    if t
      return B([ y[i, 1] for i = 1:dim(B) ])
    else
      error("Element is not in the domain")
    end
  end

  function _preimage(x::AbsAlgAssElem)
    t = zero_matrix(K, 1, dim(B))
    for i = 1:dim(B)
      t[1, i] = x.coeffs[i]
    end
    tt = t*MM
    return A([ tt[1, i] for i = 1:dim(A) ])
  end

  AtoB.header.image = _image
  AtoB.header.preimage = _preimage
  return B, AtoB
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct WindowsConfiguration <: SwaggerModel
    provisionVMAgent::Any # spec type: Union{ Nothing, Bool } # spec name: provisionVMAgent
    enableAutomaticUpdates::Any # spec type: Union{ Nothing, Bool } # spec name: enableAutomaticUpdates
    timeZone::Any # spec type: Union{ Nothing, String } # spec name: timeZone
    additionalUnattendContent::Any # spec type: Union{ Nothing, Vector{AdditionalUnattendContent} } # spec name: additionalUnattendContent
    patchSettings::Any # spec type: Union{ Nothing, PatchSettings } # spec name: patchSettings
    winRM::Any # spec type: Union{ Nothing, WinRMConfiguration } # spec name: winRM

    function WindowsConfiguration(;provisionVMAgent=nothing, enableAutomaticUpdates=nothing, timeZone=nothing, additionalUnattendContent=nothing, patchSettings=nothing, winRM=nothing)
        o = new()
        validate_property(WindowsConfiguration, Symbol("provisionVMAgent"), provisionVMAgent)
        setfield!(o, Symbol("provisionVMAgent"), provisionVMAgent)
        validate_property(WindowsConfiguration, Symbol("enableAutomaticUpdates"), enableAutomaticUpdates)
        setfield!(o, Symbol("enableAutomaticUpdates"), enableAutomaticUpdates)
        validate_property(WindowsConfiguration, Symbol("timeZone"), timeZone)
        setfield!(o, Symbol("timeZone"), timeZone)
        validate_property(WindowsConfiguration, Symbol("additionalUnattendContent"), additionalUnattendContent)
        setfield!(o, Symbol("additionalUnattendContent"), additionalUnattendContent)
        validate_property(WindowsConfiguration, Symbol("patchSettings"), patchSettings)
        setfield!(o, Symbol("patchSettings"), patchSettings)
        validate_property(WindowsConfiguration, Symbol("winRM"), winRM)
        setfield!(o, Symbol("winRM"), winRM)
        o
    end
end # type WindowsConfiguration

const _property_map_WindowsConfiguration = Dict{Symbol,Symbol}(Symbol("provisionVMAgent")=>Symbol("provisionVMAgent"), Symbol("enableAutomaticUpdates")=>Symbol("enableAutomaticUpdates"), Symbol("timeZone")=>Symbol("timeZone"), Symbol("additionalUnattendContent")=>Symbol("additionalUnattendContent"), Symbol("patchSettings")=>Symbol("patchSettings"), Symbol("winRM")=>Symbol("winRM"))
const _property_types_WindowsConfiguration = Dict{Symbol,String}(Symbol("provisionVMAgent")=>"Bool", Symbol("enableAutomaticUpdates")=>"Bool", Symbol("timeZone")=>"String", Symbol("additionalUnattendContent")=>"Vector{AdditionalUnattendContent}", Symbol("patchSettings")=>"PatchSettings", Symbol("winRM")=>"WinRMConfiguration")
Base.propertynames(::Type{ WindowsConfiguration }) = collect(keys(_property_map_WindowsConfiguration))
Swagger.property_type(::Type{ WindowsConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_WindowsConfiguration[name]))}
Swagger.field_name(::Type{ WindowsConfiguration }, property_name::Symbol) =  _property_map_WindowsConfiguration[property_name]

function check_required(o::WindowsConfiguration)
    true
end

function validate_property(::Type{ WindowsConfiguration }, name::Symbol, val)
end


text:
## ------------------------------------------------------
_parse_gids(syms::Vector{Symbol}, s::Symbol) = push!(syms, s)
_parse_gids(syms::Vector{Symbol}, qn::QuoteNode) = _parse_gids(syms, qn.value)
function _parse_gids(syms::Vector{Symbol}, ex::Expr)
    if Meta.isexpr(ex, [:(.)])
        _parse_gids(syms, ex.args[1])
        _parse_gids(syms, ex.args[2])
    else
        _gids_invalid_syntax_err(ex)
    end
end
_gids_invalid_syntax_err(ex) = error("Invalid syntax, use 'gid.gid...' or 'gid=gid.gid...'. At ex: ", ex)
_parse_gids(::Vector{Symbol}, obj) = _gids_invalid_syntax_err(obj)

## ------------------------------------------------------
function _lglob_macro_ex(Proj, exs...)
    # collect dat
    index = Dict{Symbol, Vector{Symbol}}()
    for ex in exs
        gids = Symbol[]
        if (ex isa Symbol) || Meta.isexpr(ex, [:(.)])
            # @lglob ider
            # @lglob ider1.ider2
            _parse_gids(gids, ex)
            var = last(gids)
        elseif Meta.isexpr(ex, [:(=)])
            # @lglob ider2=...
            lhand = ex.args[1]
            rhand = ex.args[2]

            !(lhand isa Symbol) && _gids_invalid_syntax_err(ex)
            var = lhand
            _parse_gids(gids, rhand)
        else
            _gids_invalid_syntax_err(ex)
        end

        haskey(index, var) && error("Name collision are not allowed. Name: '", var,
            "'. Exs: '", join(string.(index[var]), "."), "' and '",  string(ex), "'."
        )
        index[var] = gids
    end

    # expression
    ex = quote
        !($(esc(Proj)) isa Module) && error($(Proj), "must be a Module identifier")
    end
    for (ider, gids) in index
        ex = quote
            $(ex)
            $(esc(ider)) = $(lglob)($(esc(Proj)), $(gids)...)
        end
    end
    ex = quote
        $(ex)
        nothing
    end

    return ex
end

## ------------------------------------------------------
_sglob_invalid_syntax_err(ex) = error("Invalid syntax, use 'Proj gid' or 'Proj gid.gid...=Expr'. At ex: ", ex)

function _sglob_macro_ex(Proj, exs...)

    # collect dat
    index = Dict{Vector{Symbol}, Any}()
    for ex in exs
        gids = Symbol[]
        if (ex isa Symbol) || Meta.isexpr(ex, [:(.)])
            # @sglob ider2
            _parse_gids(gids, ex)
            index[gids] = last(gids)
        elseif Meta.isexpr(ex, [:(=)])
            # @sglob ider2=Expr
            lhand = ex.args[1]
            rhand = ex.args[2]

            _parse_gids(gids, lhand)
            index[gids] = rhand
        else
            _sglob_invalid_syntax_err(ex)
        end
    end

    # expression
    ex = quote
        !($(esc(Proj)) isa Module) && error($(Proj), "must be a Module identifier")
    end
    for (gids, datexpr) in index
        ex = quote
            $(ex)
            $(sglob)($(esc(Proj)), $(esc(datexpr)), $(gids)...)
        end
    end
    ex = quote
        $(ex)
        nothing
    end

    return ex
end


text:
"""
Chebyshev basis on the Interval [0,1]
"""

using ..BasisDefinition, ..Mod1DynamicDefinition, ..DynamicDefinition
using ValidatedNumerics
import ..BasisDefinition: one_vector, integral_covector, is_integral_preserving

struct Chebyshev{T<:AbstractVector} <:Basis
	p::T
	# TODO: check in constructor that p is sorted, starts with 0 and ends with 1
end

#we suppose ValidatedNumerics computes the Chebyshev points with adequate precision
ChebCouples(n, T) = hcat(
	[Interval{T}(pi); (reverse([Interval{T}(0.0); [j*Interval{T}(pi)/n for j in 1:n-1] ]))],
	[Interval{T}(0.0) ; (reverse([Interval{T}(1.0); [cos(j*Interval{T}(pi)/n) for j in 1:n-1]]).+1)/2])

ChebPoints(n, T) = ChebCouples(n, T)[:, 2]

Chebyshev(n::Integer; T = Float64) = Chebyshev(ChebPoints(n, T))

"""
Return the size of the Chebyshev basis
"""
Base.length(B::Chebyshev) = length(B.p)

"""
Eval the Chebyshev polynomial up to degree n on an array of
points in [-1, 1].

Not satisfactory, the intervals explode

"""
function _eval_T(n, x::Array{T}) where {T}
	k = length(x)
	M = zeros(T, k, n+1)
	M[:, 1] = ones(k)
	M[:, 2] = x
	for i in 3:n+1
		M[:, i] = (2*x).*M[:, i-1]-M[:, i-2]
	end
	return M
end

""" Eval a polynomial in Chebyshev basis, ClenshawBackward, using ball arithmetic
Following Viviane Ledoux, Guillaume Moroz
"Evaluation of Chebyshev polynomials on intervals andapplication to root finding"
"""
function eval_Clenshaw_BackwardFirst(coeff::Vector{Interval{S}}, x::Interval{T}) where {S,T}
	coeff_a = mid.(coeff)
	coeff_r = radius.(coeff)
	a, r = midpoint_radius(x)
	m = length(coeff)
	u = zeros(Interval{T}, m+1)
	ϵ = zeros(Interval{T}, m+1)
	u[m] = coeff_a[m]
	for k in reverse(2:m-1)
		u_temp = 2*a*u[k+1]-u[k+2]+Interval{T}(coeff_a[k])
		u[k], ϵ[k] = midpoint_radius(u_temp)
	end
	u_temp = a*u[2]-u[3]+Interval{T}(coeff_a[1])
	u[1], ϵ[1] = midpoint_radius(u_temp)

	e = zeros(Interval{T}, m+1)
	e[m] = coeff_r[m]
	for k in reverse(2:m-1)
		e[k] = e[k+1]+2*r*abs(u[k+1])+ϵ[k]+coeff_r[k]
	end
	e[1] = e[2]+r*abs(u[1])+ϵ[1]+coeff_r[1]
	γ = e[1].hi
	return u[1]+Interval(-γ, γ)
end
eval_Clenshaw_BackwardFirst(coeff::Vector{Float64}, x::Interval) = eval_Clenshaw_BackwardFirst(Interval.(coeff), x)

function eval_Clenshaw_BackwardSecond(coeff::Vector{Interval{S}}, x::Interval{T}) where {S, T}
	coeff_a = mid.(coeff)
	coeff_r = radius.(coeff)
	a, r = midpoint_radius(x)
	m = length(coeff)
	u = zeros(Interval{T}, m+1)
	ϵ = zeros(Interval{T}, m+1)
	u[m] = coeff_a[m]
	for k in reverse(2:m-1)
		u_temp = 2*a*u[k+1]-u[k+2]+Interval{T}(coeff_a[k])
		u[k], ϵ[k] = midpoint_radius(u_temp)
	end
	u_temp = 2*a*u[2]-u[3]+Interval{T}(coeff_a[1])
	u[1], ϵ[1] = midpoint_radius(u_temp)

	e = zeros(Interval{T}, m+1)
	e[m] = coeff_r[m]
	for k in reverse(2:m-1)
		e[k] = e[k+1]+(k+1)*(2*r*abs(u[k+1])+ϵ[k]+coeff_r[k])
	end
	e[1] = e[2]+2*r*abs(u[1])+ϵ[1]+coeff_r[1]
	γ = e[1].hi
	return u[1]+Interval(-γ, γ)
end

function Clenshaw(coeff, x)
	n = length(coeff)
	u = zeros(typeof(x), n+1)
	u[n] = coeff[n]
	for k in reverse(2:n-1)
		u[k] = coeff[k]+2*x*u[k+1]-u[k+2]
	end
	u[1] = coeff[1]+x*u[2]-u[3]
	return u[1]
end

function ClenshawSecond(coeff, x::T) where {T<:Real}
	n = length(coeff)
	u = zeros(T, n+1)
	u[n] = coeff[n]
	for k in reverse(2:n-1)
		u[k] = coeff[k]+2*x*u[k+1]-u[k+2]
	end
	u[1] = coeff[1]+2*x*u[2]-u[3]
	return u[1]
end

Clenshaw(coeff, x::Interval{T}) where {T} = eval_Clenshaw_BackwardFirst(coeff, x)
function ChebyshevDerivative(coeff, x::Interval{T}) where {T}
	n = length(coeff)
	coeff_der = [(i-1)*Interval{T}(coeff[i]) for i in 2:n]
	#@info coeff
	#@info coeff_der
	return eval_Clenshaw_BackwardSecond(coeff_der, x)
end

evalChebyshev(coeff, x::Interval) = eval_Clenshaw_BackwardFirst(coeff, Interval(mid.(2*x-1)))
evalChebyshevDerivative(coeff, x::Interval) = 2*ChebyshevDerivative(coeff, 2*x-1)
function evalChebyschevCentered(coeff, x::Interval)
	m = Interval(mid.(x))
	return evalChebyshev(coeff, m)+evalChebyshevDerivative(coeff, x)*(x-m)
end

using IntervalOptimisation

function infnormoffunction(B::Chebyshev, v)
	val = 0
	try
		val = maximize(x->abs(evalChebyschevCentered(v, x)), Interval(0, 1))[1]
	catch
		print("Refining grid")
		f(x) = abs(evalChebyshevCentered(v, x))
		ran = range_estimate(f, Interval(0,1), 5)
		Bval = union(val, ran)
	end
	return val
end

function infnormofderivative(B::Chebyshev, v)
	val = Interval(0)
	try
		val = maximize(x->abs(evalChebyshevDerivative(v, x)), Interval(0, 1))[1]
	catch
		print("Refining grid")
		f(x) = abs(evalChebyshevDerivative(v, x))
		ran = range_estimate(f, Interval(0,1), 5)
		val = union(val, ran)
	end
	return val
end

import .C2BasisDefinition: C1Norm

C1Norm(B::Chebyshev, v) = infnormoffunction(B,v)+infnormofderivative(B,v)
rescaling_factor(B::Chebyshev) = log(length(B)+1)

Base.length(S::AverageZero{T}) where T<:Chebyshev = length(S.basis)-1

function Base.iterate(S::AverageZero{T}, state = 1) where T<:Chebyshev
	B = S.basis
	i = state
	if i == length(B)
		return nothing
	end
	v = zeros(length(B))
	v[i+1] = 1
	v[1] = -mid.(integral_covector(B)[i+1])
	return (v, 2*i*i), state+1
end

# makes so that B[j] returns the (i+1)-th Chebyshev polynomial
# on [0, 1]
"""
Returns a function that computes the values of the (i-1)-th Chebyshev polynomial
on [0, 1]
"""
function Base.getindex(B::Chebyshev, i::Int)
 	n = length(B)
 	@boundscheck 1 <= i <= n+1 || throw(BoundsError(B, i))
 	v = zeros(n+1)
	v[i] = 1
	return x-> Clenshaw(v, 2*x-1)
end

struct ChebyshevDual <: Dual
    x::Vector{Interval} #TODO: a more generic type may be needed in future
    xlabel::Vector{Int}
    x′::Vector{Interval}
end

function ChebDualBranch(y, br::Branch, ylabel = 1:length(y), ϵ = 0.0)
	if br.increasing
		endpoint_X = br.X[2]
		der = Contractors.derivative(br.f)(endpoint_X)
		preim_der = preimages_and_derivatives(y, br, ylabel, ϵ)
		return [preim_der[1]; endpoint_X],
			[preim_der[2]; length(preim_der[2])+1],
			[preim_der[3]; der]
	else
		endpoint_X = br.X[2]
		der = Contractors.derivative(br.f)(endpoint_X)
		preim_der = preimages_and_derivatives(B.p, D, 1:length(B.p)-1, ϵ)
		return [preim_der[1]; endpoint_X],
			[preim_der[2]; length(preim_with_der[2])+1],
			[preim_der[3]; der]
	end
end

function Dual(B::Chebyshev, D::PwMap, ϵ)
	@assert is_full_branch(D)
    results = collect(ChebDualBranch(B.p, b, 1:length(B.p)-1, ϵ) for b in branches(D))
    x = vcat((result[1] for result in results)...)
    xlabel = vcat((result[2] for result in results)...)
    x′ = vcat((result[3] for result in results)...)
    return x, xlabel, x′
end

Base.length(dual::ChebyshevDual) = length(dual.x)
Base.eltype(dual::ChebyshevDual) = Tuple{eltype(dual.xlabel), Tuple{eltype(dual.x), eltype(dual.x′)}}
function iterate(dual::ChebyshevDual, state=1)
    if state <= length(dual.x)
        return ((dual.xlabel[state], (dual.x[state], abs(dual.x′[state]))), state+1)
    else
        return nothing
    end
end

using FFTW
function chebtransform(w)
	#@info sum(w)
	n = length(w)-1
	z= fft([reverse(w); w[2: end-1]])/n
	#@info z
	t = real.(z[1:length(w)])
	t[1]/=2
	t[end]/=2
	return Interval.(t)
end

using ProgressMeter
function assemble(B::Chebyshev, D::Dynamic, ϵ=0.0; T = Float64)
	n = length(B.p)
	M = zeros(Interval{T}, (n, n))
	x, labels, x′ = Dual(B, D, ϵ)
	@showprogress for i in 1:n
		ϕ = B[i]
		w = zeros(Interval{Float64}, n)
		for j in 1:length(x)
			w[labels[j]]+=ϕ(x[j])/abs(x′[j])
		end
		#@info w
		M[:, i] = chebtransform(mid.(w))
	end
	return M
end

is_integral_preserving(B::Chebyshev) = false
integral_covector(B::Chebyshev; T= Float64) = [Interval{T}(1); 0; [0.5*Interval{T}((-1)^n+1)/(1-n^2) for n in 2:length(B)-1]]'
one_vector(B::Chebyshev) = [1; zeros(length(B)-1)]

# """
# Return (in an iterator) the pairs (i, (x, |T'(x)|)) where x is a preimage of p[i], which
# describe the "dual" L* evaluation(p[i
# """
# function Base.iterate(S::DualComposedWithDynamic{T, D}, state = (1, 1)) where T<:Hat where D<:Dynamic
# 	@assert is_full_branch(S.dynamic)

# 	i, k = state

# 	if i == length(S.basis)+1
# 			return nothing
# 	end

# 	x = preim(S.dynamic, k, S.basis.p[i], S.ϵ)
# 	absT′ = abs(derivative(S.dynamic, x))

# 	if k == nbranches(S.dynamic)
# 		return ((i, (x, absT′)), (i+1, 1))
# 	else
# 		return ((i, (x, absT′)), (i, k+1))
# 	end
# end

# function BasisDefinition.is_dual_element_empty(::Hat, d)
# 	# TODO: the preim() may indeed be empty, so there could be an additional check here
# 	return false
# end

# BasisDefinition.is_refinement(Bf::Hat, Bc::Hat) = Bc.p ⊆ Bf.p

# function integral_covector(B::Hat)
# 	n = length(B)
# 	return 1/n * ones(Interval{Float64}, n)'
# end

# function one_vector(B::Hat)
# 	return ones(length(B))
# end

# """
# Return the range of indices of the elements of the basis whose support intersects
# with the given dual element (i.e., a pair (y, absT')).
# The range may end with length(B)+1; this must be interpreted "mod length(B)":
# it means that it intersects with the hat function peaked in 0 as well
# (think for instance y = 0.9999).
# """
# function BasisDefinition.nonzero_on(B::Hat, dual_element)
# 	y, absT′ = dual_element
# 	# Note that this cannot rely on arithmetic unless it is verified

# 	y = y ∩ Interval(0.,1.) # we assume it's bona-fide interval in [0,1]
# 	# this should work for preims(), since they are supposed to return
# 	# a number in [0,1]

# 	# finds in which semi-open interval [p[k], p[k+1]) y.lo and y.hi fall
# 	lo = searchsortedlast(B.p, y.lo)
# 	hi = searchsortedlast(B.p, y.hi)
# 	lo = min(lo, length(B)) # lo may be n+1 if y.lo==1
# 	hi = min(hi, length(B)) # hi may be n+1 if y.hi==1
# 	hi = hi + 1 # because the hat centered in p[k] is also nonzero in the interval before

# 	if lo == 1 # 1:N+1 does not make sense and would mean that the first interval is counted twice
# 		hi = min(hi, length(B))
# 	end
# 	return (lo, hi)
# end

# """
# Given a preimage ```y``` of a point ```x```, this iterator returns
# ```\\phi_j(y)/T'(y) ```
# """
# function Base.iterate(S::ProjectDualElement{T,DT}, state = S.j_min) where {T <: Hat,DT}
# 	if state == S.j_max+1
# 		return nothing
# 	end
# 	y, absT′ = S.dual_element
# 	j = state
# 	y_normalized = IntervalOnTorus(y)
# 	n = length(S.basis)

# 	return ((j, S.basis[mod(j, 1:n)](y_normalized) / absT′),
# 		    state+1)
# end

# BasisDefinition.strong_norm(B::Hat) = Lipschitz
# BasisDefinition.weak_norm(B::Hat) = Linf
# BasisDefinition.aux_norm(B::Hat) = L1

# evaluate_integral(B::Hat, i, T) = T(i)/length(B)

# function Base.iterate(S::AverageZero{Hat}, state = 1)
# 	n = length(S.basis)
# 	if state == n
# 		return nothing
# 	end
# 	v = zeros(Float64, n)
# 	v[1] = 1
# 	v[state+1]=-1
# 	return (v, state+1)
# end

# BasisDefinition.weak_projection_error(B::Hat) = 0.5 ⊘₊ Float64(length(B), RoundDown)
# BasisDefinition.aux_normalized_projection_error(B::Hat) = 0.5 ⊘₊ Float64(length(B), RoundDown)
# BasisDefinition.strong_weak_bound(B::Hat) = 2. ⊗₊ Float64(length(B), RoundDown)
# BasisDefinition.aux_weak_bound(B::Hat) = 1.
# BasisDefinition.weak_by_strong_and_aux_bound(B::Hat) = (1., 1.)
# BasisDefinition.bound_weak_norm_from_linalg_norm(B::Hat) = @error "TODO"
# BasisDefinition.bound_linalg_norm_L1_from_weak(B::Hat) = @error "TODO"
# BasisDefinition.bound_linalg_norm_L∞_from_weak(B::Hat) = @error "TODO"

# function BasisDefinition.invariant_measure_strong_norm_bound(B::Hat, D::Dynamic)
# 	A, B = dfly(strong_norm(B), aux_norm(B), D)
# 	@assert A < 1.
# 	return B ⊘₊ (1. ⊖₋ A)
# end

# using RecipesBase

# """
# Plots a function in the Hat basis
# """
# @recipe function f(B::Hat, w::AbstractVector)

# 	legend --> :bottomright

# 	if eltype(w) <: Interval
# 		w = mid.(w)
# 	end

# 	@series begin
# 		seriestype --> :path
# 		label --> L"f_{\delta}"
# 		ylims --> (0, NaN)
# 		B.p, vcat(w, w[end])
# 	end
# end

# """
# Displays error on a function in the Hat basis
# """
# @recipe function f(B::Hat, error::Number, w)

# 	if eltype(w) <: Interval
# 		w = mid.(w)
# 	end

# 	if isfinite(error)
# 		@series begin
# 			seriestype --> :path
# 			seriesalpha --> 0.5
# 			fillrange --> vcat(w, w[end]) .- error
# 			label --> "Error area"
# 			B.p, vcat(w, w[end]) .+ error
# 		end
# 	end
# end


text:
const Prime = UInt64
abstract type PrimesDictAbstract end # used for forward reference
mutable struct PrimesDict <: PrimesDictAbstract
    sieve :: Dict{Prime,Prime}
    baseprimes :: PrimesDictAbstract
    lastbaseprime :: Prime
    q :: Prime
    PrimesDict() = new(Dict())
end
Base.eltype(::Type{PrimesDict}) = Prime
Base.IteratorSize(::Type{PrimesDict}) = Base.SizeUnknown() # "infinite"...
function Base.iterate(PD::PrimesDict, state::Prime = Prime(0) )
    if state < 1
        PD.baseprimes = PrimesDict()
        PD.lastbaseprime = Prime(3)
        PD.q = Prime(9)
        return Prime(2), Prime(1)
    end
    dict = PD.sieve
    while true
        state += 2
        if !haskey(dict, state)
            state < PD.q && return state, state
            p = PD.lastbaseprime # now, state = PD.q in all cases
            adv = p + p # since state is at PD.q, advance to next
            dict[state + adv] = adv # adds base prime composite stream
            # following initializes secondary base strea first time
            p <= 3 && Base.iterate(PD.baseprimes)
            p = Base.iterate(PD.baseprimes, p)[1] # next base prime
            PD.lastbaseprime = p
            PD.q = p * p
        else # advance hit composite in dictionary...
            adv = pop!(dict, state)
            next = state + adv
            while haskey(dict, next) next += adv end
            dict[next] = adv # past other composite hits in dictionary
        end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Basic" begin

    # test mock that accepts no args
    @mockable quack() = "quack"
    apply(quack => () -> "woof") do
        @test quack() == "woof"
    end

    # test mock that accepts some args
    @mockable add(x, y) = x + y
    apply(add => (x, y) -> 0) do
        @test add(1, 2) == 0
    end

    # test mock with fallback (conditional patch)
    apply(add => (x, y) -> x == 1 ? 0 : Fallback()) do
        @test add(1, 2) == 0          # patch returend value
        @test add(2, 2) == 4          # fall back to default
    end

    # test multiple patches
    @mockable mul2(x, y::Int)     = x * y + 2
    @mockable mul2(x, y::Float64) = x * y - 2
    apply(mul2 => (x, y::Int)     -> x == 1 ? 0 : Fallback(),
          mul2 => (x, y::Float64) -> x ^ y) do
        @test mul2(1, 2) == 0          # value returned from patch
        @test mul2(2, 2) == 6          # fall back to default
        @test mul2(2, 3.0) == 8.0      # Float64 patch
    end

    # test functions that uses kwargs
    @mockable div2(x; n = 2) = x / n
    @test div2(10) == 5                # without patch
    apply(div2 => (x; n = 3) -> n) do  # patched
        # This case may be unintutive. Because `n` is not specified,
        # the original value of `n` is passed (n = 2). Hence, the absense
        # of `n` in the call isn't really absent and therefore overrides
        # the different default value in the patch. Hence, it should be
        # advised that
        @test div2(10) == 2
        @test div2(10; n = 4) == 4
    end

    # This is the correct way; patch should not have any default values for kwargs
    apply(div2 => (x; n) -> n) do
        @test div2(10) == 2
        @test div2(10; n = 4) == 4
    end

    # Expect failure when patch does not take kwargs
    apply(div2 => (x) -> 0) do
        if VERSION >= v"1.2"
            @test_throws MethodError div2(10) == 0
        else
            @test_throws ErrorException div2(10) == 0
        end
    end

    # test splatting with Any type
    @mockable splat1(k...) = 1
    apply(splat1 => (k...) -> 2) do
        @test splat1(1,2,3) == 2
    end

    # test splatting with specific type
    @mockable splat2(k::Int...) = 1
    apply(splat2 => (k::Int...) -> 2) do
        @test splat2(1,2,3) == 2
    end

    # Negative test: signature does not match so it cannot be applied
    @mockable splat3(x, y) = x + y
    apply(splat3 => (args...) -> 10) do
        @test splat3(1, 2) == 3
    end

    # Test keyword arg splatting
    @mockable splat_kwarg(; kwargs...) = -1
    apply(splat_kwarg => (; kwargs...) -> length(kwargs)) do
        @test splat_kwarg()               == 0
        @test splat_kwarg(; x = 1)        == 1
        @test splat_kwarg(; x = 1, y = 2) == 2
    end

    # Test unnamed args

    # simple
    @mockable unnamed1(::Int) = 1
    apply(unnamed1 => (::Int) -> 2) do
        @test unnamed1(1) == 2
    end

    # any position
    @mockable unnamed2(x, ::Int, y) = 1
    apply(unnamed2 => (x, ::Int, y) -> 3) do
        @test unnamed2(:a, 1, 2.0) == 3
    end

    # more than one unnamed arg
    @mockable unnamed3(x, ::Int, ::Float64) = 1
    apply(unnamed3 => (x, ::Int, ::Float64) -> 4) do
        @test unnamed3(:a, 1, 2.0) == 4
    end

    # patch has name
    @mockable unnamed4(::Int) = 1
    apply(unnamed4 => (x::Int) -> 5) do
        @test unnamed4(2) == 5
    end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1d4d4c0e-1f81-11ec-1e2a-d9eb1f505571
begin
	using Pkg
	Pkg.activate("..")
	using Revise, Javis, JUtils, Colors, Distributions, Animations, LaTeXStrings
	import Luxor
end

# ╔═╡ 9dd75dc1-dff3-4b4c-945c-cc169ca6e555
mywhite = RGB(254/255, 252/255, 244/255)

# ╔═╡ 82c29842-f1e7-4357-a187-8c75a0714681
function ground(color)
	(args...) -> begin
		background(color)
		sethue("black")
	end
end

# ╔═╡ 4bfea714-7fdd-47cf-8e54-bb30bc2a625c
let
	height = 250
	width = 700
	margin = 50
	n_frames = 200

	prelearn_video = Video(width, height)
	Background(1:n_frames, ground(mywhite))
	mypoints = make_prelearn(n_frames=n_frames, height=height, width=width, margin=margin)
	render(prelearn_video, pathname="../output/prelearn_animation.gif")
end

# ╔═╡ 21e78860-f28b-4475-b38f-01e4ee7d802e
let
	n_frames = 200
	time_steps = range(1, n_frames, length = 6)
	width = 400
	height = 400

	v = Video(300, 300)
	Background(1:n_frames, ground(mywhite))
	make_ner(
		n_frames = n_frames,
		width = width,
		height = height,
		words_in = ["anybody"],
		words_out = ["driver", "golf player", "dog", "astronaut", "rider"],
		margin=50
	)
	render(v, pathname="../output/paper_user.gif")
end

# ╔═╡ 553312d1-74ef-4edc-b293-caf06a65b71b
let
	n_frames = 200
	time_steps = range(1, n_frames, length = 6)
	width = 500
	height = 500

	v = Video(400, 400)
	Background(1:n_frames, ground(mywhite))
	make_ner(
		n_frames=n_frames,
		width=width,
		height=height,
		words_in=["machine"],
		words_out = [
			"coating machine", "coating device", "knife", "computer", "car"
		],
		margin=50
	)
	render(v, pathname="../output/TechNER.gif")
end

# ╔═╡ 67328dd2-4a0b-428c-a4ff-687d76498b69
function postprocess_frame(frame_image, frame, frames)
	T = typeof(frame_image[1,1])
	g(x) = RGB(1 - x.r, 1 - x.g, 1 - x.b)
	convert.(T, g.(frame_image))
end

# ╔═╡ c12aff7b-e936-4e8e-98e4-7f8acde3cf35
function make_neural_net(; n_frames, width, height, margin, n_layers, maxneurons)

	marg_width = width ÷ 2 - margin
	marg_height = height ÷ 2 - margin
	xpos = range(-marg_width, marg_width, length = n_layers)
	ypos = range(-marg_height, marg_height, length = maxneurons)
	neurons_per_layers = [maxneurons - idx for idx in 0:n_layers-1]
	dist = Normal(0, 50)

	layers = [
		[Point(x, y) for y in
				range(
					-(marg_height - 3marg_height/2n_neurons),
					(marg_height - 3marg_height/2n_neurons),
					length=n_neurons
				)
			]
		for (n_neurons, x) in zip(neurons_per_layers, xpos)
	]

	Object(@JShape begin
			for layer in layers
				for neuron in layer
					sethue("red")
					circle(neuron, 15, :fill)
					setcolor("black")
					setline(1)
					circle(neuron, 15, :stroke)
				end
			end
			for i in 1:length(layers) - 1
				for (p1, p2) in Base.product(layers[i], layers[i+1])
					if p1 != p2
						(n1, i1, i2) = intersectionlinecircle(
							p1, p2, p1, 15
						)
						(n2, j1, j2) = intersectionlinecircle(
							p1, p2, p2, 15
						)
						p1 = i1
						p2 = j2
						arrow(p1, p2)
					end
				end
			end
		end
	)
end

# ╔═╡ f429a061-3ed7-41df-b39d-0989bf477da9
begin
	n_frames = 2000
	width = 1600
	height = 900
	margin = 50
	time_steps = floor.(Int, range(1, n_frames, length = 12))

	nn_video = Video(width, height)
	Background(
		1:n_frames,
		ground(mywhite)
	)
	nn = @JLayer 1:n_frames begin
		make_neural_net(
			n_frames = n_frames,
			width=400,
			height=height,
			margin=margin,
			n_layers=3,
			maxneurons=6
		)
	end
	act!(nn, Action(time_steps[2]:time_steps[3], anim_scale(0.5)))
	act!(nn, Action(time_steps[2]:time_steps[3], anim_translate(Point(-width÷3, 0))))
	act!(nn, Action(time_steps[2]:time_steps[3], anim_rotate(-π/2)))
	act!(nn, Action(time_steps[4]:time_steps[5], disappear(:fade)))

	myarrow = Object(
		time_steps[3]:time_steps[5],
		@JShape begin
			arrow(Point(-width ÷ 6, 0), Point(width ÷ 6, 0))
		end
	)

	act!(
		myarrow,
		[
			Action(1:time_steps[2], appear(:fade)),
			Action(time_steps[2]:time_steps[3], disappear(:fade))
		]
	)

	wlstartpoints = [Point(i + width÷3, 0) for i in [-150, -90, -30, 30, 90, 150]]
	wlendpoints = [Point(i, -height÷3) for i in [-180, -105, -35, 35, 105, 180]]
	wlstrings = [
		L"\Large{w_1}",
		L"\Large{w_2}",
		L"\Large{w_i}",
		L"\Large{w_{308}}",
		L"\Large{w_j}",
		L"\Large{w_n}"
		]

	latexvecs = [@JLayer time_steps[4]:time_steps[end] 500 500 wlstartpoint begin
		texvec = Object(
			1:n_frames,
			@JShape begin
				Javis.latex(
						wlstring,
						O,
						:stroke,
						valign=:middle,
						halign=:center
					)
				end
			)
			act!(texvec, Action(1:1, anim_scale(2.0)))
			act!(texvec, Action(time_steps[2]:time_steps[3], anim_scale(2.5)))
		end for (wlstartpoint, wlstring) in zip(wlstartpoints, wlstrings)]

	vecendpoint = Point(0, -height÷3)
	for (idx, latexvec) in enumerate(latexvecs)
		act!(
			latexvec,
			Action(time_steps[2]:time_steps[3], anim_translate(wlstartpoints[idx], wlendpoints[idx] ))
		)
	end
	plots_points = [Point.(
		[rand(-width÷4:width÷4) for i in 1:30],
		[rand(-height÷4:height÷4) for i in 1:30]
	) for i in 1:5]
	plots_rotations = rand(range(0, π/4, length=20), 5)
	plots_colors = [randomhue() for _ in 1:5]
	linear_regressions = @JLayer time_steps[6]:time_steps[end] begin
		for (θ, plot_points, color) in zip(plots_rotations, plots_points, plots_colors)
			Object(
				1:n_frames,
				(args...) -> begin
					bb = box(O, width÷2, height÷2, :stroke)
					rotate(-θ)
					sethue(color)
					# sethue("red")
					line(
						Point(-width÷5, -height÷5),
						Point(width÷5, height÷5),
						:stroke
					)
					# sethue("black")
					for p in plot_points
						rotate(θ)
						if (abs(p.x) > width÷4) || (abs(p.y) > height÷4)
							continue
						else
							circle(p, 3, :fill)
						end
						rotate(-θ)
					end
				end
			)
		end
	end

	act!(linear_regressions, Action(time_steps[2]:time_steps[3], anim_scale(0.5)))
	act!(linear_regressions, Action(time_steps[2]:time_steps[3], anim_translate(Point(width÷3, 0))))

	central_arr = Object(time_steps[7]:time_steps[9],
		@JShape begin
			arrow(vecendpoint + (0, 20), Point(0, height÷3 - 20))
		end
	)
	act!(central_arr, Action(1:time_steps[2], appear(:fade)))

	barpoints = [Point(i, height÷3) for i in [-180, -105, -35, 35, 105, 180]]
	lstrings = [
		L"$\Large{\alpha_{1}}$",
		L"$\Large{\alpha_{2}}$",
		L"$\Large{\alpha_{i}}$",
		L"$\Large{\alpha_{308}}$",
		L"$\Large{\alpha_{j}}$",
		L"$\Large{\alpha_{n}}$"
	]
	texparams = [@JLayer time_steps[6]:time_steps[end] 300 300 p begin
		ob = Object(
		1:n_frames,
		@JShape begin
			Javis.latex(
				Lstr,
				O,
				:stroke,
				valign=:middle,
				halign=:center
			)
		end
	)
		act!(ob, Action(1:1, anim_scale(2.5), keep=true))
		act!(ob, Action(1:time_steps[2], appear(:fade)))
		end
			for (Lstr, p) in zip(lstrings, barpoints)
			]

	[Object(
		time_steps[9]:time_steps[12],
		JLine(x - (0, 10), x - (0, y), linewidth=5, color="red")
	) for (x, y) in zip(barpoints, [80, 40, 12, 200, 60, 20])]

	render(
		nn_video, pathname="../output/bert_lasso.mp4",
		postprocess_frame=postprocess_frame
	)
end

# ╔═╡ Cell order:
# ╠═1d4d4c0e-1f81-11ec-1e2a-d9eb1f505571
# ╟─9dd75dc1-dff3-4b4c-945c-cc169ca6e555
# ╟─82c29842-f1e7-4357-a187-8c75a0714681
# ╟─4bfea714-7fdd-47cf-8e54-bb30bc2a625c
# ╟─21e78860-f28b-4475-b38f-01e4ee7d802e
# ╟─553312d1-74ef-4edc-b293-caf06a65b71b
# ╠═67328dd2-4a0b-428c-a4ff-687d76498b69
# ╟─c12aff7b-e936-4e8e-98e4-7f8acde3cf35
# ╠═f429a061-3ed7-41df-b39d-0989bf477da9


text:
# This script uses the convex penalties with cross-validated lambdas (as obtained in the excel sheet) on a test set
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Compat, Random
include("convexpenalties.jl")

results_template = DataFrame(n=Int[],p=Int[], m=Int[], k=Int[], theSeed=Int[], gamma=Real[], lambda_persp=Real[], lambda_dcl=Real[], lambda_nn=Real[],
gap_persp=Real[], gap_dcl=Real[], t_persp=Real[], t_dcl=Real[], t_nn=Real[], acc_persp=Real[], acc_dcl=Real[], acc_nn=Real[],
mse_persp=Real[], mse_dcl=Real[], mse_nn=Real[], rnk_persp=Int[], rnk_dcl=Int[], rnk_nn=Int[])

results_exact=similar(results_template, 0)

for ARG in ARGS

array_num = parse(Int, ARG)
seeds=1:20

n=50
p=50
ms=5:5:100
k=10
gamma=1e6

lambdas_persp=[0.0001,	0.0001,	0.00784759970351461,	0.00784759970351461,	0.00784759970351461,	0.615848211066026,	0.615848211066026,	0.233572146909012,	0.0885866790410082,	0.00784759970351461,	0.0335981828628378,	0.143844988828766,	0.233572146909012,	0.379269019073224,	0.379269019073224,	0.379269019073224,	0.615848211066026,	0.615848211066026,	0.615848211066026,	0.615848211066026]
lambdas_dcl=[0.0127427498570313,	0.00112883789168468,	0.00784759970351461,	0.0206913808111479,	0.0206913808111479,	0.00483293023857175,	0.00183298071083243,	0.0206913808111479,	0.00297635144163131,	0.0545559478116852,	0.0206913808111479,	0.0206913808111479,	0.00297635144163131,	0.0127427498570313,	0.00483293023857175,	0.00483293023857175,	0.00784759970351461,	0.00784759970351461,	0.00784759970351461,	0.00483293023857175]
lambdas_nn=[0.000001,	0.000001,	0.000001,	0.000001,	0.000001,	0.000885866790410082,	0.000885866790410082,	0.000545559478116852,	0.000335981828628378,	0.000078475997035146,	0.000335981828628378,	0.000545559478116852,	0.000545559478116852,	0.000885866790410082,	0.000885866790410082,	0.000885866790410082,	0.000885866790410082,	0.000885866790410082,	0.00143844988828766,	0.00143844988828766]

n_obs_oos=1000
lambda_persp=lambdas_persp[(array_num-1)%20+1]
lambda_dcl=lambdas_dcl[(array_num-1)%20+1]
lambda_nn=lambdas_nn[(array_num-1)%20+1]

m=ms[(array_num-1)%20+1]

seed_mult=seeds[(array_num-1)÷20+1]

for theSeed in (100+10*(seed_mult-1)+1):(100+10*seed_mult)
        Random.seed!(theSeed)
        X=rand(m,p)
        X_oos=rand(n_obs_oos, p)

        U_true=rand(p,k)
        V_true=rand(n,k)
        beta_true=U_true*V_true'
        E=0.05*rand(m,n)
        E_oos=0.05*rand(n_obs_oos,n)
        Y=X*beta_true.+E
        Y_oos=X_oos*beta_true.+E_oos

        t_persp=@elapsed beta_est,ofv_dual, ofv_primal,=solvepersprelax(X, Y, m, n, p, lambda_persp, gamma)
        @show t_persp
        @show gap_persp=100*abs.(ofv_primal-ofv_dual)/ofv_primal
        @show rel_error_persp=norm(beta_est.-beta_true)/norm(beta_true)
        @show ofv_primal
        @show 1.0/(2.0*gamma)*norm(beta_est)^2
        @show mse_persp=norm(X_oos*beta_est.-Y_oos)^2/norm(X_oos*beta_true.-Y_oos)^2

        u,v,=svd(beta_est)
        @show rnk_persp=sum(v.>1e-4)

        t_dcl=@elapsed beta_est,ofv_dual, ofv_primal,= solvepersprelax_DCL(X, Y, m, n, p, lambda_dcl, gamma)
        @show gap_dcl=100*abs.(ofv_primal-ofv_dual)/ofv_primal
        @show t_dcl
        @show rel_error_dcl=norm(beta_est.-beta_true)/norm(beta_true)

        u,v,=svd(beta_est)
        @show rnk_dcl=sum(v.>1e-4)
        @show mse_dcl=norm(X_oos*beta_est.-Y_oos)^2/norm(X_oos*beta_true.-Y_oos)^2

        t_nn=@elapsed beta_est,=solveNNrelax(X, Y, m, n, p, lambda_nn, gamma)
        @show t_nn
        @show rel_error_nn=norm(beta_est.-beta_true)/norm(beta_true)
        u,v,=svd(beta_est)
        @show rnk_nn=sum(v.>1e-4)
        @show mse_nn=norm(X_oos*beta_est.-Y_oos)^2/norm(X_oos*beta_true.-Y_oos)^2

        push!(results_exact, [n, p, m, k, theSeed, gamma, lambda_persp, lambda_dcl, lambda_nn,
        gap_persp, gap_dcl, t_persp, t_dcl, t_nn, rel_error_persp, rel_error_dcl, rel_error_nn,
        mse_persp, mse_dcl, mse_nn, rnk_persp, rnk_dcl, rnk_nn])

end
CSV.write("createFig1_withfixedlambdas.csv", results_exact, append=true)

end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Statistics, Plots
	Random.seed!(2)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.14"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	n, N, alpha = 10, 10^7, 0.05
	mActual = 0.75
	dist0, dist1 = Uniform(0,1), Uniform(0,mActual)

	ts(sample) = maximum(sample) - minimum(sample)
end

# ╔═╡ d8d55f04-3af8-11eb-2251-75efd46ad5a2
begin
	empiricalDistUnderH0 = [ts(rand(dist0,n)) for _ in 1:N]
	rejectionValue = quantile(empiricalDistUnderH0,alpha)
end;

# ╔═╡ d8d5a0b8-3af8-11eb-1cc3-4db334734491
begin
	sample = rand(dist1,n)
	testStat = ts(sample)
	pValue = sum(empiricalDistUnderH0 .<= testStat)/N

	if testStat > rejectionValue
		print("Didn't reject: ", round(testStat,digits=4))
		print(" > ", round(rejectionValue,digits=4))
	else
		print("Reject: ", round(testStat,digits=4))
		print(" <= ", round(rejectionValue,digits=4))
	end
	Text("p-value = $(round(pValue,digits=4))")
end

# ╔═╡ d8d6283a-3af8-11eb-37cd-5fa24b45e40f
begin
	stephist(empiricalDistUnderH0, bins=100, c=:blue, normed=true, label="")
	plot!([testStat, testStat], [0,4], c=:red, label="Observed test statistic")
	plot!([rejectionValue, rejectionValue], [0,4], c=:black, ls=:dash,
		label="Critical value boundary", legend=:topleft, ylims=(0,4),
			xlabel = "x", ylabel = "Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.14"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═d8d55f04-3af8-11eb-2251-75efd46ad5a2
# ╠═d8d5a0b8-3af8-11eb-1cc3-4db334734491
# ╠═d8d6283a-3af8-11eb-37cd-5fa24b45e40f
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using MDCT, Test

# O(N^2) implementations straight from the definitions

function slow_mdct(x::AbstractVector{T}) where {T<:Number}
    n = length(x)
    m = div(n,2)
    y = Array{T}(undef, m)
    for j = 1:m
        y[j] = sum(x .* cos.(2*pi/n * (j-0.5) * ((1:n) .- 0.5 .+ n/4)))
    end
    return y
end

function slow_imdct(y::AbstractVector{T}) where {T<:Number}
    m = length(y)
    n = m*2
    x = Array{T}(undef, n)
    for k = 1:n
        x[k] = 2/n * sum(y .* cos.(2*pi/n * ((1:m).-0.5) * (k - 0.5 + n/4)))
    end
    return x
end

for N=100:4:108
    X = rand(N)
    Y = mdct(X)
    Z = imdct(Y)
    Yp = plan_mdct(X)*X
    Zp = plan_imdct(Y)*Y
    Ys = slow_mdct(X)
    Zs = slow_imdct(Y)
    @test Y ≈ Ys rtol=1e-13
    @test Z ≈ Zs rtol=1e-13
    @test Yp ≈ Ys rtol=1e-13
    @test Zp ≈ Zs rtol=1e-13
end


text:
@testset "3_of_6" begin
    a = SVector{6,Float64}(1,2,3,4,5,6)
    @test all(first_3_of_6(a) .== SVector{3,Float64}(1,2,3))
    @test all(last_3_of_6(a) .== SVector{3,Float64}(4,5,6))
end

@testset "utility" begin
    p1 = SVector{3,Float64}(1.0, 2.0, 3.0)
    p2 = SVector{3,Float64}(2.0, 3.0, 4.0)
    @test p2 == weightPoly(p1, p2, 1.0, 0.0)
    @test p1 == weightPoly(p1, p2, 0.0, 1.0)
    @test (p1 + p2) * 0.5 == weightPoly(p1, p2, -0.7, 0.7)

    p, d = make_pd_gains(2*pi, 1.0)
    @test p == 1.0
    @test d == 2.0
    p, d = make_pd_gains(1*pi, 1.0)
    @test p == 4.0
    @test d == 4.0

    for k = 1:8
        pd33 = rand_pd(k)
        @test (isa(pd33, Hermitian{Float64,SMatrix{k,k,Float64,k^2}}))
        @test (0.0 < minimum(eigen(pd33).values))
    end
end

@testset "padding" begin
    v3 = SVector{3,Float64}(7, 8, 9)
    v3_1 = SVector{4,Float64}(v3..., 1)
    v3_0 = SVector{4,Float64}(v3..., 0)

    @test unPad(v3_1) == v3
    @test unPad(SMatrix{1,4,Float64,4}(v3_1...)) == v3
    @test onePad(v3) == v3_1
    @test zeroPad(v3) == v3_0

    M = SMatrix{4,4,Float64,16}(randn(16))
    @test mul_then_un_pad(M, v3_1) == unPad(M * v3_1)
    @test one_pad_then_mul(M, v3) == M * onePad(v3)
end


text:
# utils_Taylor1.jl: 1-variable Taylor expansions
#
# Last modification: 2014.06.07
#
# Luis Benet & David P. Sanders
# UNAM
#

## Constructors ##
immutable Taylor{T<:Number} <: AbstractSeries{T,1}
    coeffs :: Array{T,1}
    order :: Int
    ## Inner constructor ##
    function Taylor(coeffs::Array{T,1}, order::Int)
        lencoef = length(coeffs)
        order = max(order, lencoef-1)
        v = zeros(T, order+1)
        #@inbounds v[1:lencoef] = coeffs[1:lencoef]
        for i = 1:lencoef
            @inbounds v[i] = coeffs[i]
        end
        new(v, order)
    end
end
## Outer constructors ##
Taylor{T<:Number}(x::Taylor{T}, order::Int) = Taylor{T}(x.coeffs, order)
Taylor{T<:Number}(x::Taylor{T}) = Taylor{T}(x.coeffs, x.order)
Taylor{T<:Number}(coeffs::Array{T,1}, order::Int) = Taylor{T}(coeffs, order)
Taylor{T<:Number}(coeffs::Array{T,1}) = Taylor{T}(coeffs, length(coeffs)-1)
Taylor{T<:Number}(x::T, order::Int) = Taylor{T}([x], order)
Taylor{T<:Number}(x::T) = Taylor{T}([x], 0)

## Type, length ##
eltype{T<:Number}(::Taylor{T}) = T
length{T<:Number}(a::Taylor{T}) = a.order

## Conversion and promotion rules ##
convert{T<:Number}(::Type{Taylor{T}}, a::Taylor) = Taylor(convert(Array{T,1}, a.coeffs), a.order)
convert{T<:Number, S<:Number}(::Type{Taylor{T}}, b::Array{S,1}) = Taylor(convert(Array{T,1},b))
convert{T<:Number}(::Type{Taylor{T}}, b::Number) = Taylor([convert(T,b)], 0)
promote_rule{T<:Number, S<:Number}(::Type{Taylor{T}}, ::Type{Taylor{S}}) = Taylor{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{Taylor{T}}, ::Type{Array{S,1}}) = Taylor{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{Taylor{T}}, ::Type{S}) = Taylor{promote_type(T, S)}

## Auxiliary function ##
function firstnonzero{T<:Number}(a::Taylor{T})
    order = a.order
    nonzero::Int = order+1
    z = zero(T)
    for i = 1:order+1
        if a.coeffs[i] != z
            nonzero = i-1
            break
        end
    end
    nonzero
end
function fixshape{T<:Number, S<:Number}(a::Taylor{T}, b::Taylor{S})
    order = max(a.order, b.order)
    a1, b1 = promote(a, b)
    return Taylor(a1, order), Taylor(b1, order), order
end

## real, imag, conj and ctranspose ##
for f in (:real, :imag, :conj)
    @eval ($f){T<:Number}(a::Taylor{T}) = Taylor(($f)(a.coeffs), a.order)
end
ctranspose{T<:Number}(a::Taylor{T}) = conj(a)

## zero and one ##
zero{T<:Number}(a::Taylor{T}) = Taylor(zero(T), a.order)
one{T<:Number}(a::Taylor{T}) = Taylor(one(T), a.order)

## Equality ##
function ==(a::Taylor, b::Taylor)
    a1, b1, order = fixshape(a, b)
    return a1.coeffs == b1.coeffs
end

## Addition and substraction ##
for f in (:+, :-)
    @eval begin
        function ($f)(a::Taylor, b::Taylor)
            a1, b1, order = fixshape(a, b)
            v = ($f)(a1.coeffs, b1.coeffs)
            return Taylor(v, order)
        end
        ($f)(a::Taylor) = Taylor(($f)(a.coeffs), a.order)
    end
end

## Multiplication ##
function *(a::Taylor, b::Taylor)
    a1, b1, order = fixshape(a, b)
    T = eltype(a1)
    coeffs = zeros(T,order+1)
    coeffs[1] = a1.coeffs[1] * b1.coeffs[1]
    for k = 1:order
        coeffs[k+1] = mulHomogCoef(k, a1.coeffs, b1.coeffs)
    end
    Taylor(coeffs, order)
end
# Homogeneous coefficient for the multiplication
function mulHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, bc::Array{T,1})
    kcoef == 0 && return ac[1] * bc[1]
    coefhomog = zero(T)
    for i = 0:kcoef
        @inbounds coefhomog += ac[i+1] * bc[kcoef-i+1]
    end
    coefhomog
end

## Division ##
function /(a::Taylor, b::Taylor)
    a1, b1, order = fixshape(a, b)
    orddivfact, cdivfact = divfactorization(a1, b1) # order and coefficient of first factorized term
    T = typeof(cdivfact)
    v1 = convert(Array{T,1}, a1.coeffs)
    v2 = convert(Array{T,1}, b1.coeffs)
    coeffs = zeros(T,order+1)
    coeffs[1] = cdivfact
    for k = orddivfact+1:order
        coeffs[k-orddivfact+1] = divHomogCoef(k, v1, v2, coeffs, orddivfact)
    end
    Taylor(coeffs, order)
end
function divfactorization(a1::Taylor, b1::Taylor)
    # order of first factorized term; a1 and b1 are assumed to be of the same order (length)
    a1nz = firstnonzero(a1)
    b1nz = firstnonzero(b1)
    orddivfact = min(a1nz, b1nz)
    if orddivfact > a1.order
        orddivfact = a1.order
    end
    cdivfact = a1.coeffs[orddivfact+1] / b1.coeffs[orddivfact+1]
    aux = abs2(cdivfact)
    # Is the polynomial factorizable?
    if isinf(aux) || isnan(aux)
        info("Order k=$(orddivfact) => coeff[$(orddivfact+1)]=$(cdivfact)")
        error("Division does not define a Taylor polynomial\n",
            " or its first non-zero coefficient is Inf/NaN.\n")
    ##else orddivfact>0
    ##    warn("Factorizing the polynomial.\n",
    ##        "The last k=$(orddivfact) Taylor coefficients ARE SET to 0.\n")
    end
    return orddivfact, cdivfact
end
# Homogeneous coefficient for the division
function divHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, bc::Array{T,1},
    coeffs::Array{T,1}, ordfact::Int)
    #
    kcoef == ordfact && return ac[ordfact+1] / bc[ordfact+1]
    coefhomog = mulHomogCoef(kcoef, coeffs, bc)
    coefhomog = (ac[kcoef+1]-coefhomog) / bc[ordfact+1]
    coefhomog
end

## Division functions: rem and mod
for op in (:mod, :rem)
    @eval begin
        function ($op){T<:Real}(a::Taylor{T}, x::Real)
            coeffs = a.coeffs
            coeffs[1] = ($op)(a.coeffs[1], x)
            return Taylor( coeffs, a.order)
        end
    end
end
function mod2pi{T<:Real}(a::Taylor{T})
    coeffs = a.coeffs
    coeffs[1] = mod2pi( a.coeffs[1] )
    return Taylor( coeffs, a.order)
end

## Int power ##
function ^(a::Taylor, n::Integer)
    uno = one(a)
    n < 0 && return uno / a^(-n)
    n == 0 && return uno
    if n%2 == 0     # even power
        n == 2 && return square(a)
        pow = div(n, 2)
        return square( a^pow )
    else            # odd power
        n == 1 && return a
        pow = div(n-1, 2)
        return a*square( a^pow )
    end
end
## Rational power ##
^(a::Taylor,x::Rational) = a^(x.num/x.den)
## Real power ##
function ^(a::Taylor, x::Real)
    uno = one(a)
    if x == zero(x)
        return uno
    elseif x == 0.5
        return sqrt(a)
    end
    order = a.order
    # First non-zero coefficient
    l0nz = firstnonzero(a)
    l0nz > order && return zero(a)
    # The first non-zero coefficient of the result; must be integer
    lnull = x*l0nz
    !isinteger(lnull) &&
        error("Integer exponent REQUIRED if the Taylor polynomial is expanded around 0.\n")
    # Reaching this point, it is possible to implement the power of the Taylor polynomial.
    # The last l0nz coefficients are set to zero.
    lnull = itrunc(lnull)
    ##l0nz > 0 && warn("The last k=$(l0nz) Taylor coefficients ARE SET to 0.\n")
    aux = (a.coeffs[l0nz+1])^x
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, order+1)
    coeffs[lnull+1] = aux
    k0 = lnull+l0nz
    for k = k0+1:order
        @inbounds coeffs[k-l0nz+1] = powHomogCoef(k, v, x, coeffs, l0nz)
    end
    Taylor(coeffs,order)
end
# Homogeneous coefficients for real power
function powHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, x::Real,
    coeffs::Array{T,1}, knull::Int)
    #
    kcoef == knull && return (ac[knull+1])^x
    coefhomog = zero(T)
    for i = 0:kcoef-knull-1
        aux = x*(kcoef-i)-i
        @inbounds coefhomog += aux*ac[kcoef-i+1]*coeffs[i+1]
    end
    aux = kcoef - knull*(x+1)
    @inbounds coefhomog = coefhomog / (aux*ac[knull+1])
    coefhomog
end
^{T<:Number,S<:Real}(a::Taylor{T}, x::Complex{S}) = exp( x*log(a) )
^(a::Taylor, b::Taylor) = exp( b*log(a) )

## Square ##
function square{T<:Number}(a::Taylor{T})
    order = a.order
    coeffs = zeros(T,order+1)
    coeffs[1] = a.coeffs[1]^2
    for k = 1:order
        @inbounds coeffs[k+1] = squareHomogCoef(k, a.coeffs)
    end
    Taylor(coeffs,order)
end
# Homogeneous coefficients for square
function squareHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1})
    kcoef == 0 && return ac[1]^2
    coefhomog = zero(T)
    kodd = kcoef%2
    kend = div(kcoef - 2 + kodd, 2)
    for i = 0:kend
        @inbounds coefhomog += ac[i+1]*ac[kcoef-i+1]
    end
    coefhomog = 2coefhomog
    if kodd == 0
        @inbounds coefhomog += ac[div(kcoef,2)+1]^2
    end
    coefhomog
end

## Square root ##
function sqrt(a::Taylor)
    order = a.order
    # First non-zero coefficient
    l0nz = firstnonzero(a)
    if l0nz > order
        return Taylor(zero(T),order)
    elseif l0nz%2 == 1 # l0nz must be pair
        error("First non-vanishing Taylor coefficient must be an EVEN POWER\n",
            "to expand SQRT around 0.\n")
    end
    # Reaching this point, it is possible to implement the sqrt of the Taylor polynomial.
    # The last l0nz coefficients are set to zero.
    ##l0nz > 0 && warn("The last k=$(l0nz) Taylor coefficients ARE SET to 0.\n")
    lnull = div(l0nz, 2)
    aux = sqrt(a.coeffs[l0nz+1])
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, order+1)
    @inbounds coeffs[lnull+1] = aux
    for k = lnull+1:order-l0nz
        @inbounds coeffs[k+1] = sqrtHomogCoef(k, v, coeffs, lnull)
    end
    Taylor(coeffs, order)
end
# Homogeneous coefficients for the square-root
function sqrtHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffs::Array{T,1}, knull::Int)
    kcoef == knull && return sqrt(ac[2*knull+1])
    coefhomog = zero(T)
    kodd = (kcoef - knull)%2
    kend = div(kcoef - knull - 2 + kodd, 2)
    for i = knull+1:knull+kend
        @inbounds coefhomog += coeffs[i+1]*coeffs[kcoef+knull-i+1]
    end
    @inbounds aux = ac[kcoef+knull+1]-2coefhomog
    if kodd == 0
        @inbounds aux = aux - (coeffs[kend+knull+2])^2
    end
    @inbounds coefhomog = aux / (2coeffs[knull+1])
    coefhomog
end

## Exp ##
function exp(a::Taylor)
    order = a.order
    aux = exp( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, order+1)
    coeffs[1] = aux
    for k = 1:order
        @inbounds coeffs[k+1] = expHomogCoef(k, v, coeffs)
    end
    Taylor( coeffs, order )
end
# Homogeneous coefficients for exp
function expHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffs::Array{T,1})
    kcoef == 0 && return exp(ac[1])
    coefhomog = zero(T)
    for i = 0:kcoef-1
        @inbounds coefhomog += (kcoef-i) * ac[kcoef-i+1] * coeffs[i+1]
    end
    coefhomog = coefhomog/kcoef
    coefhomog
end

## Log ##
function log(a::Taylor)
    order = a.order
    l0nz = firstnonzero(a)
    if firstnonzero(a)>0
        error("Not possible to expand LOG around 0.\n")
    end
    aux = log( a.coeffs[1] )
    T = typeof(aux)
    ac = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, order+1)
    coeffs[1] = aux
    for k = 1:order
        @inbounds coeffs[k+1] = logHomogCoef(k, ac, coeffs)
    end
    Taylor( coeffs, order )
end
# Homogeneous coefficients for log
function logHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffs::Array{T,1})
    kcoef == 0 && return log( ac[1] )
    coefhomog = zero(T)
    for i = 1:kcoef-1
        @inbounds coefhomog += (kcoef-i) * ac[i+1] * coeffs[kcoef-i+1]
    end
    @inbounds coefhomog = (ac[kcoef+1] -coefhomog/kcoef) / ac[1]
    coefhomog
end

## Sin and cos ##
sin(a::Taylor) = sincos(a)[1]
cos(a::Taylor) = sincos(a)[2]
function sincos(a::Taylor)
    order = a.order
    aux = sin( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    sincoeffs = zeros(T,order+1)
    coscoeffs = zeros(T,order+1)
    sincoeffs[1] = aux
    coscoeffs[1] = cos( a.coeffs[1] )
    for k = 1:order
        @inbounds sincoeffs[k+1], coscoeffs[k+1] = sincosHomogCoef(k, v, sincoeffs, coscoeffs)
    end
    return Taylor(sincoeffs, order), Taylor(coscoeffs, order)
end
# Homogeneous coefficients for log
function sincosHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1},
    scoeffs::Array{T,1}, ccoeffs::Array{T,1})
    #
    kcoef == 0 && return sin( ac[1] ), cos( ac[1] )
    sincoefhom = zero(T)
    coscoefhom = zero(T)
    for i = 1:kcoef
        @inbounds begin
            x = i * ac[i+1]
            sincoefhom += x * ccoeffs[kcoef-i+1]
            coscoefhom -= x * scoeffs[kcoef-i+1]
        end
    end
    sincoefhom = sincoefhom/kcoef
    coscoefhom = coscoefhom/kcoef
    return sincoefhom, coscoefhom
end

## Tan ##
function tan(a::Taylor)
    order = a.order
    aux = tan( a.coeffs[1] )
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T,order+1)
    coeffst2 = zeros(T,order+1)
    coeffs[1] = aux
    coeffst2[1] = aux^2
    for k = 1:order
        @inbounds begin
            coeffs[k+1] = tanHomogCoef(k, v, coeffst2)
            coeffst2[k+1] = squareHomogCoef(k, coeffs)
        end
    end
    Taylor( coeffs, order )
end
# Homogeneous coefficients for tan
function tanHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, coeffst2::Array{T,1})
    kcoef == 0 && return tan( ac[1] )
    coefhomog = zero(T)
    for i = 0:kcoef-1
        @inbounds coefhomog += (kcoef-i)*ac[kcoef-i+1]*coeffst2[i+1]
    end
    @inbounds coefhomog = ac[kcoef+1] + coefhomog/kcoef
    coefhomog
end

## Differentiating ##
function diffTaylor{T<:Number}(a::Taylor{T})
    order = a.order
    coeffs = zeros(T, order+1)
    coeffs[1] = a.coeffs[2]
    for i = 1:order
        @inbounds coeffs[i] = i*a.coeffs[i+1]
    end
    return Taylor(coeffs, order)
end

## Integrating ##
function integTaylor{T<:Number}(a::Taylor{T}, x::Number)
    order = a.order
    R = promote_type(T, typeof(a.coeffs[1] / 1), typeof(x))
    coeffs = zeros(R, order+1)
    for i = 1:order
        @inbounds coeffs[i+1] = a.coeffs[i] / i
    end
    coeffs[1] = convert(R, x)
    return Taylor(coeffs, order)
end
integTaylor{T<:Number}(a::Taylor{T}) = integTaylor(a, zero(T))
# function definiteIntegralTaylor{T<:Number}(a::Taylor{T}, x1::Number, x2::Number)
#     @assert isfinite(x1) && isfinite(x2)
#     integral = integTaylor(a)
#     v1 = evalTaylor(integral, x1)
#     v2 = evalTaylor(integral, x2)
#     return v2-v1
# end

## Evaluates a Taylor polynomial on a given point using Horner's rule ##
function evalTaylor{T<:Number,S<:Number}(a::Taylor{T}, dx::S)
    R = promote_type(T,S)
    orden = a.order
    @inbounds suma = convert(R, a.coeffs[end])
    for k = orden:-1:1
        @inbounds suma = suma*dx + a.coeffs[k]
    end
    suma
end
evalTaylor{T<:Number}(a::Taylor{T}) = evalTaylor(a, zero(T))

function evalTaylor{T<:Number,S<:Number}(a::Taylor{T}, x::Taylor{S})
    a, x, orden = fixshape(a, x)
    @inbounds suma = a.coeffs[end]
    for k = orden:-1:1
        @inbounds suma = suma*x + a.coeffs[k]
    end
    suma
end

## Returns de n-th derivative of a series expansion
function deriv{T}(a::Taylor{T}, n::Int=1)
    @assert n>= 0
    n > a.order && error(
        "You need to increase the order of the Taylor series (currently ", a.order, ")\n",
        "to calculate its ", n,"-th derivative.")
    res::T = factorial( widen(n) ) * a.coeffs[n+1]
    return res
end


text:
export EQLoc

abstract type ComputedLoc end

function show(io::IO, Loc::T) where {T <: ComputedLoc}
  if get(io, :compact, false) == false
    println(io, T, " with fields:")
    for f in fieldnames(T)
      fn = lpad(String(f), 5, " ")
      if f == :flags
        println(io, fn, ": ", bitstring(getfield(Loc, f)))
      else
        println(io, fn, ": ", getfield(Loc, f))
      end
    end
  else
    c = :compact => true
    print(io, repr(getfield(Loc, :lat), context=c), " N, ",
              repr(getfield(Loc, :lon), context=c), " E, ",
              repr(getfield(Loc, :dep), context=c), " km")
  end
  return nothing
end

isequal(Loc1::T, Loc2::T) where {T <: ComputedLoc} = minimum(
  [isequal(getfield(Loc1, f), getfield(Loc2, f)) for f in fieldnames(T)] )
==(Loc1::T, Loc2::T) where {T <: ComputedLoc} = isequal(Loc1, Loc2)

function hash(Loc::T) where {T<:ComputedLoc}
  h = hash(zero(UInt64))
  for f in fieldnames(T)
    h = hash(getfield(Loc, f), h)
  end
  return h
end

"""
    EQLoc

QuakeML-compliant earthquake location

Field     | Type    | Meaning       | SeisIO conventions/behavior
--------: |:------- |:--------------| :----------
 lat      | Float64 | latitude      | °N = +
 lon      | Float64 | longitude     | °E = +
 dep      | Float64 | depth         | km; down = +
 dx       | Float64 | x error       | uses units of data source (typically km)
 dy       | Float64 | y error       | uses units of data source (typically km)
 dz       | Float64 | z error       | uses units of data source (typically km)
 dt       | Float64 | ot error      | uses units of data source (typically s)
 se       | Float64 | std error     | uses units of data source (typically s)
 rms      | Float64 | rms pick err  | uses units of data source (typically s)
 gap      | Float64 | azimuthal gap | uses units of data source (typically °)
 dmin     | Float64 | min sta dist  | uses units of data source (typically km)
 dmax     | Float64 | max sta dist  | uses units of data source (typically km)
 nst      | Int64   | # of stations |
 flags    | UInt8   | boolean flags | access flag[n] with >>(<<(flags,n-1),7)
 datum    | String  | geog. datum   |
 typ      | String  | location type | freeform (e.g. "centroid", "hypocenter")
 sig      | String  | significance  | freeform (e.g. "95%", "2σ")
          |         | / confidence  |
 src      | String  | source        | freeform (e.g. "HYPOELLIPSE", "HypoDD")

flags (0x01 = true, 0x00 = false)
1. x fixed?
2. y fixed?
3. z fixed?
4. t fixed?
"""
mutable struct EQLoc <: ComputedLoc
  lat   ::Float64
  lon   ::Float64
  dep   ::Float64
  dx    ::Float64
  dy    ::Float64
  dz    ::Float64
  dt    ::Float64
  se    ::Float64
  rms   ::Float64
  gap   ::Float64
  dmin  ::Float64
  dmax  ::Float64
  nst   ::Int64
  flags ::UInt8
  datum ::String
  typ   ::String
  sig   ::String
  src   ::String

  function EQLoc(
                  lat   ::Float64,
                  lon   ::Float64,
                  dep   ::Float64,
                  dx    ::Float64,
                  dy    ::Float64,
                  dz    ::Float64,
                  dt    ::Float64,
                  se    ::Float64,
                  rms   ::Float64,
                  gap   ::Float64,
                  dmin  ::Float64,
                  dmax  ::Float64,
                  nst   ::Int64,
                  flags ::UInt8,
                  datum ::String,
                  typ   ::String,
                  sig   ::String,
                  src   ::String,
                  )
    return new(lat, lon, dep, dx, dy, dz, dt, se, rms, gap, dmin, dmax, nst, flags, datum, typ, sig, src)
  end
end

EQLoc(;
        lat   ::Float64                   = zero(Float64),
        lon   ::Float64                   = zero(Float64),
        dep   ::Float64                   = zero(Float64),
        dx    ::Float64                   = zero(Float64),
        dy    ::Float64                   = zero(Float64),
        dz    ::Float64                   = zero(Float64),
        dt    ::Float64                   = zero(Float64),
        se    ::Float64                   = zero(Float64),
        rms   ::Float64                   = zero(Float64),
        gap   ::Float64                   = zero(Float64),
        dmin  ::Float64                   = zero(Float64),
        dmax  ::Float64                   = zero(Float64),
        nst   ::Int64                     = zero(Int64),
        flags ::UInt8                     = 0x00,
        datum ::String                    = "",
        typ   ::String                    = "",
        sig   ::String                    = "",
        src   ::String                    = "",
        ) = EQLoc(lat, lon, dep, dx, dy, dz, dt, se, rms, gap, dmin, dmax, nst, flags, datum, typ, sig, src)

function write(io::IO, Loc::EQLoc)
  for f in (:lat, :lon, :dep, :dx, :dy, :dz, :dt, :se, :rms, :gap, :dmin, :dmax, :nst, :flags)
    write(io, getfield(Loc, f))
  end
  write(io, sizeof(Loc.datum))
  write(io, Loc.datum)
  write(io, sizeof(Loc.typ))
  write(io, Loc.typ)
  write(io, sizeof(Loc.sig))
  write(io, Loc.sig)
  write(io, sizeof(Loc.src))
  write(io, Loc.src)
  return nothing
end

read(io::IO, ::Type{EQLoc}) = EQLoc(read!(io, Array{Float64,1}(undef, 12))...,
  fastread(io, Int64),
  fastread(io),
  String(fastread(io, fastread(io, Int64))),
  String(fastread(io, fastread(io, Int64))),
  String(fastread(io, fastread(io, Int64))),
  String(fastread(io, fastread(io, Int64)))
  )

function isempty(Loc::EQLoc)
  q::Bool = min(isempty(getfield(Loc, :datum)),
                isempty(getfield(Loc, :typ)),
                isempty(getfield(Loc, :sig)),
                isempty(getfield(Loc, :src)),
                getfield(Loc, :nst) == zero(Int64),
                getfield(Loc, :flags) == 0x00)
  for f in (:lat, :lon, :dep, :dx, :dy, :dz, :dt, :se, :rms, :gap, :dmin, :dmax)
    q = min(q, getfield(Loc, f) == zero(Float64))
  end
  return q
end

sizeof(Loc::EQLoc) = 233 +
                      sizeof(getfield(Loc, :datum)) +
                      sizeof(getfield(Loc, :typ)) +
                      sizeof(getfield(Loc, :sig)) +
                      sizeof(getfield(Loc, :src))