text:
"""
    ovlp = overlap(distsq, s, w, dirdot)

Calculates the unnormalized overlap between two Gaussian distributions with width `s`,
weight `w', and squared distance `distsq`, and geometric scaling factor `dirdot`.
"""
function overlap(distsq::Real, s::Real, w::Real, dirdot::Real)
    return w * 0.5*(1+dirdot) * exp(-distsq / (2*s)) # / (sqrt2pi * sqrt(s))^ndims
    # Note, the normalization term for the Gaussians is left out, since it is not required that the total "volume" of each Gaussian
    # is equal to 1 (e.g. satisfying the requirements for a probability distribution)
end

"""
    ovlp = overlap(dist, σx, σy, ϕx, ϕy, dirdot)

Calculates the unnormalized overlap between two Gaussian distributions with variances
`σx` and `σy`, weights `ϕx` and `ϕy`, and means separated by distance `dist`, scaled
by the dot product obtained from geometric constraints `dirdot`.
"""
function overlap(dist::Real, σx::Real, σy::Real, ϕx::Real, ϕy::Real, dirdot::Real)
    return overlap(dist^2, σx^2 + σy^2, ϕx*ϕy, dirdot)
end

"""
    ovlp = overlap(x::IsotropicGaussian, y::IsotropicGaussian, xtform=identity)

Calculates the unnormalized overlap between two `IsotropicGaussian` objects.
"""
function overlap(x::AbstractIsotropicGaussian, y::AbstractIsotropicGaussian, s=nothing, w=nothing)
    if isnothing(s)
        s = x.σ^2 + y.σ^2
    end
    if isnothing(w)
        w = x.ϕ*y.ϕ
    end
    if length(x.dirs) == 0 || length(y.dirs) == 0
        cosγ = 1.
    else
        cosγ = -1
        for xdir in x.dirs
            for ydir in y.dirs
                cosγ = max(cosγ, dot(xdir, ydir))
            end
        end
    end
    return overlap(sum(abs2, x.μ.-y.μ), s, w, cosγ)
end

"""
    ovlp = overlap(x::AbstractSingleGMM, y::AbstractSingleGMM, xtform=identity)

Calculates the unnormalized overlap between two `AbstractSingleGMM` objects.
"""
function overlap(x::AbstractSingleGMM, y::AbstractSingleGMM, pσ=nothing, pϕ=nothing)
    # prepare pairwise widths and weights, if not provided
    if isnothing(pσ) || isnothing(pϕ)
        pσ, pϕ = pairwise_consts(x, y)
    end

    # sum overlaps for all pairwise combinations of Gaussians between x and y
    ovlp = zero(promote_type(numbertype(x),numbertype(y)))
    for (i,gx) in enumerate(x.gaussians)
        for (j,gy) in enumerate(y.gaussians)
            ovlp += overlap(gx, gy, pσ[i,j], pϕ[i,j])
        end
    end
    return ovlp
end

"""
    ovlp = overlap(x::AbstractMultiGMM, y::AbstractMultiGMM, xtform=identity)

Calculates the unnormalized overlap between two `AbstractMultiGMM` objects.
"""
function overlap(x::AbstractMultiGMM, y::AbstractMultiGMM, mpσ=nothing, mpϕ=nothing)
    # prepare pairwise widths and weights, if not provided
    if isnothing(mpσ) || isnothing(mpϕ)
        mpσ, mpϕ = pairwise_consts(x, y)
    end

    # sum overlaps from each keyed pairs of GMM
    ovlp = zero(promote_type(numbertype(x),numbertype(y)))
    for k in keys(x.gmms) ∩ keys(y.gmms)
        ovlp += overlap(x.gmms[k], y.gmms[k], mpσ[k], mpϕ[k])
    end
    return ovlp
end

"""
l2dist = distance(x, y)

Calculates the L2 distance between two GMMs made up of spherical Gaussian distributions.
"""
function distance(x::AbstractGMM, y::AbstractGMM)
    return overlap(x,x) + overlap(y,y) - 2*overlap(x,y)
end

"""
tani = tanimoto(x, y)

Calculates the tanimoto distance based on Gaussian overlap between two GMMs.
"""
function tanimoto(x::AbstractGMM, y::AbstractGMM)
    o = overlap(x,y)
    return o / (overlap(x,x) + overlap(y,y) - o)
end


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

using Markdown
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.peek, el) ? Base.peek(el) : missing
        el
    end
end

# ╔═╡ 4be9e1fc-b716-11ea-1a73-afd8c1433865
using CSV, Plots, DataFrames

# ╔═╡ dd5e6e68-b717-11ea-34d7-ff403c5ea379
using Dates

# ╔═╡ 6dc31e06-b716-11ea-17f7-0717e2a64fb1
url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"

# ╔═╡ 7612f77a-b716-11ea-14df-a176322d1945
download(url, "covid_data.csv")

# ╔═╡ 7f1c3796-b716-11ea-3a9a-332f0d30960f
data = CSV.read("covid_data.csv");

# ╔═╡ 5621cc1a-b717-11ea-3753-21cd9218ddd9
rename!(data, 1 => "province", 2 => "country");

# ╔═╡ 68c6bdd0-b717-11ea-39b1-635ec4b86d9a
countries = collect(data[:, 2])

# ╔═╡ c08b5e18-b717-11ea-1ee9-3122ac3a7389
function GetCountryData(country::String, range::Int)::Vector
    countries_row = findfirst(countries .== uppercasefirst(country))

    if !isnothing(countries_row)
        country_data_row = data[countries_row, :]
        country_data = country_data_row[5:range]

        return country_data
    end

    error("Not valid country")
end

# ╔═╡ fb20bbc4-b71a-11ea-2547-93a05833219e
function FormatDates(range::Int)
    format = Dates.DateFormat("m/d/Y")
	date_strings = names(data)[5:range]
    formated_dates = parse.(Date, date_strings, format) .+ Year(2000)

    return formated_dates
end

# ╔═╡ 7faf57f6-b71b-11ea-2391-cf8db547cfe0
paises = ["peru", "mexico", "chile", "brazil"]

# ╔═╡ 1ce1c6ce-b80e-11ea-191e-df7e9fc04238
@bind range html"<input type='range' min=5>"

# ╔═╡ 09df210a-b71b-11ea-0c3f-d3ac3064c3d1
begin
	dates = FormatDates(range)
	paises_data = GetCountryData.(paises, range)
end

# ╔═╡ 1f559534-b71b-11ea-30a5-d3632489c22e
begin
	plot(
		dates,
		paises_data[1],
		xticks=dates[1:5:end],
		xrotation=45,
		leg=:topleft,
		label=uppercasefirst(paises[1]),
		m=2
	)

	for p in 2:length(paises)
		plot!(dates, paises_data[p], m=2, label=uppercasefirst(paises[p]))
	end

	xlabel!("Dates")

	ylabel!("Confirmed")

	title!("COVD-19: Confirmed cases")
end

# ╔═╡ fdb0fe48-b810-11ea-1541-15816e1ca410
begin
	@gif for i in 5:160
		dates = FormatDates(i)
		paises_data = GetCountryData.(paises, i)

		plot(
			dates,
			paises_data[1],
			xticks=dates[1:5:end],
			xrotation=45,
			leg=:topleft,
			label=uppercasefirst(paises[1]),
			m=2
		)

		for p in 2:length(paises)
			plot!(dates, paises_data[p], m=2, label=uppercasefirst(paises[p]))
		end

		xlabel!("Dates")

		ylabel!("Confirmed")

		title!("COVD-19: Confirmed cases")
	end
end

# ╔═╡ Cell order:
# ╠═4be9e1fc-b716-11ea-1a73-afd8c1433865
# ╠═6dc31e06-b716-11ea-17f7-0717e2a64fb1
# ╠═7612f77a-b716-11ea-14df-a176322d1945
# ╠═7f1c3796-b716-11ea-3a9a-332f0d30960f
# ╠═5621cc1a-b717-11ea-3753-21cd9218ddd9
# ╠═68c6bdd0-b717-11ea-39b1-635ec4b86d9a
# ╠═c08b5e18-b717-11ea-1ee9-3122ac3a7389
# ╠═dd5e6e68-b717-11ea-34d7-ff403c5ea379
# ╠═fb20bbc4-b71a-11ea-2547-93a05833219e
# ╠═7faf57f6-b71b-11ea-2391-cf8db547cfe0
# ╠═1ce1c6ce-b80e-11ea-191e-df7e9fc04238
# ╠═09df210a-b71b-11ea-0c3f-d3ac3064c3d1
# ╠═1f559534-b71b-11ea-30a5-d3632489c22e
# ╠═fdb0fe48-b810-11ea-1541-15816e1ca410

ext: jl

max_stars_repo_name: elihugarret/my-notebooks


text:
using DynamicHMC, Turing, Test

@model gdemo(x, y) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
  return s, m
end

mf = gdemo(1.5, 2.0)

chn = sample(mf, DynamicNUTS(2000));

@test mean(chn[:s]) ≈ 49/24 atol=0.2
@test mean(chn[:m]) ≈ 7/6 atol=0.2


text:
module MutableTypes
#Done
if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

export mutable, Mutable, deref

abstract type Mutable{T} end
mutable struct MutableX{T} <: Mutable{T}
    x::T
    MutableX{T}() where {T} = new{T}()
    MutableX{T}(x) where {T} = new{T}(x)
end
struct MutableA{T, N} <: Mutable{T}
    x::Array{T, N}
    i::Int
end
const  MutableV{T} = MutableA{T, 1}

mutable(x::T) where {T} = MutableX{T}(x)
mutable(x::Mutable) = x
mutable(x::Type{T}) where {T} = MutableX{T}()

function mutable(x::Array{T, N}, i = 1) where {T, N}
    if isbitstype(T)
        MutableA{T, N}(x, i)
    else
        mutable(x[i])
    end
end
mutable(x::Array{T, N}, i = 1) where {T <: Ptr, N} = mutable(x[i])
mutable(x::Ptr{T}, i = 1) where {T} = x + (i - 1) * sizeof(T)
mutable(x::T, i) where {T} = (i == 1 ? mutable(x) : error("Object only has one element"))

_addrof(b::T) where {T} = pointer_from_objref(b)
Base.cconvert(::Type{P}, b::MutableX{T}) where {P <: Ptr, T} = isbitstype(T) ? convert(P, _addrof(b)) : convert(P, _addrof(b.x))
Base.cconvert(::Type{P}, b::MutableA{T, N}) where {P <: Ptr, T, N} = convert(P, pointer(b.x, b.i))

deref(b::Ptr) = unsafe_load(b)
deref(b::MutableX) = b.x
deref(b::MutableA) = b.x[b.i]

Base.unsafe_load(b::Mutable) = deref(b)
Base.getindex(b::Mutable) = deref(b)

Base.unsafe_store!(b::MutableX, x) = (b.x = x)
Base.unsafe_store!(b::MutableA, x) = (b.x[b.i] = x)
Base.setindex!(b::MutableX{T}, x::T) where {T} = (b.x = x)
Base.setindex!(b::MutableA{T, N}, x::T) where {T, N} = (b.x[b.i] = x)

end

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
module RouteSequencing

include("init.jl")
include("data_structures.jl")
include("sequencing.jl")
include("io.jl")
include("geometry.jl")

include("Algorithms/Hamiltonian/hamiltonian.jl")
include("Algorithms/MinTime.jl")
include("Algorithms/MinTimeGlobal.jl")
include("Algorithms/MinTimeZ.jl")
include("Algorithms/MinTimeFinal.jl")

include("exports.jl")

end # module


text:
function addWeak(sort_idx, n_types_to_remove, processes_I, processesStd_I, Wall_I, genomeErrors_I, genomesReconstructed_I)

    totalMutTypes = size(Wall_I, 1) + n_types_to_remove
    processes = zeros(totalMutTypes, size(processes_I, 2))
    processesStd = zeros(totalMutTypes, size(processesStd_I, 2))
    Wall = zeros(totalMutTypes, size(Wall_I, 2))
    genomeErrors = zeros(totalMutTypes, size(genomeErrors_I, 2), size(genomeErrors_I, 3))
    genomesReconstructed = zeros(totalMutTypes, size(genomesReconstructed_I, 2), size(genomesReconstructed_I, 3))

    for i = 1 : size(Wall_I, 1)
        idx = sort_idx[i+n_types_to_remove]
        processes[idx, :] = processes_I[i, :]
        processesStd[idx, :] = processesStd_I[i, :]
        Wall[idx, :] = Wall_I[i, :]
        genomeErrors[idx, :, :] = genomeErrors_I[i, :, :]
        genomesReconstructed[idx, :, :] = genomesReconstructed_I[i, :, :]
    end

    return processes, processesStd, Wall, genomeErrors, genomesReconstructed
end

max_stars_repo_name: ypradat/sigprofilerjulia

max_issues_repo_name: ypradat/sigprofilerjulia

max_forks_repo_name: ypradat/sigprofilerjulia


text:
using PMI

include(joinpath(dirname(pathof(PMI)), "..", "examples", "distributed.jl"))

# Everything before this is run on all processes
WireUp.wireup()

# Everything after this is only run on the primary

using Distributed
using Test

@test myid() == 1
@test length(procs()) == PMI.get_size()


text:
module NormalVectors

using Gridap
using Gridap.Helpers
using Gridap.CellFields: NonIterableCellFieldLike

export NormalVector

import Base: length
import Gridap: evaluate

struct NormalVector{Z,D,J,R} <: NonIterableCellFieldLike{Z,VectorValue{D,Float64},1}
  jac_surf::J
  nvec_ref::R
end

function NormalVector(trian::BoundaryTriangulation)
  desc = trian.descriptor
  NormalVector(desc)
end

function NormalVector(trian::SkeletonTriangulation)
  desc = trian.descriptor1
  NormalVector(desc)
end

function NormalVector(desc::BoundaryDescriptor)
  phi = desc.cell_phi
  jac = gradient(phi)
  jac_surf = restrict(jac,desc)
  polytopes = desc.cell_to_polytope
  facet_to_lfacet = desc.facet_to_lfacet
  _normal_vector(jac_surf,polytopes,facet_to_lfacet)
end

length(nv::NormalVector) = length(nv.jac_surf)

function evaluate(nv::NormalVector{Z}, q::CellPoints{Z}) where Z
  jac_surf_q = evaluate(nv.jac_surf,q)
  apply(_map_normal,jac_surf_q,nv.nvec_ref,broadcast=true)
end

function _map_normal(J,n)
  v = inv(J)*n
  m = sqrt(inner(v,v))
  if m < eps()
    return zero(n)
  else
    return v/m
  end
end

function _normal_vector(
  jac_surf,
  polytopes::CellValue{<:Polytope},
  facet_to_lfacet)
  @notimplemented
end

function _normal_vector(
  jac_surf,
  polytopes::ConstantCellValue{<:Polytope{D}},
  facet_to_lfacet) where D

  p = polytopes.value
  n, _ = face_normals(p)

  _n = [ [ni,] for ni in n  ]

  nvec_ref = IndexCompressedCellValue(_n,facet_to_lfacet)

  J = typeof(jac_surf)
  R = typeof(nvec_ref)

  Z = D -1
  NormalVector{Z,D,J,R}(jac_surf,nvec_ref)

end

#@fverdugo this code works does not work here, but works when placed
# in the tests
import Base: show

function show(io::IO,self::NormalVector{Z,D,J,R}) where {Z,D,J,R}
  print(io,"$(nameof(typeof(self))){$Z,$D} object")
end

function show(io::IO,::MIME"text/plain",self::NormalVector{Z,D,J,R}) where {Z,D,J,R}
  show(io,self)
  print(io,":")
  print(io,"\n physdim: $D")
  print(io,"\n refdim: $Z")
  print(io,"\n ncells: $(length(self))")
end

end # module


text:
if !@isdefined(BodePlots)
include("BodePlots.jl")
end

module TLines

using InspectDR
import InspectDR: Plot2D
using NumericIO
using Colors

const VecOrNum = Union{Vector,Number}
const VecOrReal = Union{Vector,Real}

#Colour constants
#-------------------------------------------------------------------------------
const black = RGB24(0, 0, 0); const white = RGB24(1, 1, 1)
const red = RGB24(1, 0, 0); const green = RGB24(0, 1, 0); const blue = RGB24(0, 0, 1)

#Physics constants
#-------------------------------------------------------------------------------
const μ0 = 4pi*1e-7 #F/m
const ϵ0 = 8.854e-12 #H/m

#Convenience aliases
#-------------------------------------------------------------------------------
mag(x) = abs(x); phase(x) = rad2deg(angle(x))
SI(x) = formatted(x, :SI, ndigits=3)
const j = im

#==Equations
===============================================================================#
#=NOTE:
 - Not particularly optimized for speed.
 - Not particularly complete/well named set of functions.
=#

#Compute propagation constant of a transmission line
function γline(f::VecOrReal; α::VecOrReal=0, μ::Real=μ0, ϵ::Real=ϵ0)
	β = f*(2pi*sqrt(μ*ϵ))
	return α+j*β
end

#Compute reflection coefficient of a given impedance
#Zref: Reference impedance for the reflection coefficient
Γ(Z::VecOrNum; Zref::Real=50.0) = (Z .- Zref) ./ (Z .+ Zref)

#Compute impedance of a loaded transmission line
#ZC: Characteristic impedance (not currently modelling f-dependence)
#ZL: Load impendance (termination)
function Zline(ZC::Number, γ::VecOrNum; ZL::VecOrNum=Inf64, ℓ::Real=0.0)
	tanh_γℓ = tanh.(γ*ℓ)
	return ZC .* (ZL .+ ZC .* tanh_γℓ) ./ (ZC .+ ZL .* tanh_γℓ)
end

#==Main functions
===============================================================================#
function newplot()
	w = 500; h = w/1.6 #WANTCONST: Select plot width/height
	smithext = InspectDR.PExtents1D(min=-1.2,max=1.2) #WANTCONST: Padded a bit

	mplot = InspectDR.Multiplot(title="Loaded Transmission Line")
	mplot.layout[:ncolumns] = 2
	#Bode plot looks better with wider aspect ratio:
#	mplot.layout[:halloc_plot] = w; mplot.layout[:valloc_plot] = h/2
	#... But not Smith chart.

	plot_bode = add(mplot, Main.BodePlots.new(InspectDR.Plot))
	strip_mag, strip_phase = plot_bode.strips
		strip_mag.yext_full = InspectDR.PExtents1D(-20, 0)
		strip_phase.yext_full = InspectDR.PExtents1D(-200, 200)
		plot_bode.annotation.title = "Reflection Coefficient (Γ)"

	plot_smith = add(mplot, InspectDR.smithchart(:Z, ref=50))
	strip = plot_smith.strips[1]
		plot_smith.xext_full = smithext; strip.yext_full = smithext
		a = plot_smith.annotation
			a.title = "Z-Smith Chart"
			a.xlabel = "Real(Γ)"
			a.ylabels = ["Imaginary(Γ)"]

	return mplot
end

#Modify previously-generated plot:
#-------------------------------------------------------------------------------
#(Frequencies: in Hz)
function update(mplot::InspectDR.Multiplot, fmin, fmax, ltype::Symbol, R, L, C, ℓ, ZC, αdB)
	npts_bode = 500 #WANTCONST
	ldata = line(color=blue, width=3, style=:solid) #WANTCONST
	plot_bode = mplot.subplots[1] #WANTCONST
	plot_smith = mplot.subplots[2] #WANTCONST
	αDC = (20/log(10)) * αdB #WANTCONST: Nepers/m

	#Add some frequency dependence to α (makes plot more interesting):
	fref = 10e9; αref = 2*αDC #WANTCONST
	fdepα(f, αDC) = αDC .+ (αref/sqrt(fref)) .* sqrt.(f)

	#How to compute load impedance:
	Zload_shunt(jω,R,L,C) = 1 ./ (1 ./ R .+ 1 ./ (jω*L) .+ jω*C)
	Zload_series(jω,R,L,C) = R .+ jω*L .+ 1 ./ (jω*C)
	Zload = :shunt == ltype ? Zload_shunt : Zload_series

	#Generate Bode plot:
	f = 10 .^ range(log10(fmin), stop=log10(fmax), length=npts_bode)
	jω = j*2pi*f
	y = Γ(Zline(ZC, γline(f, α=fdepα(f, αDC)), ZL=Zload(jω,R,L,C), ℓ=ℓ))
	Main.BodePlots.update(plot_bode, f, y)
	strip_mag, strip_phase = plot_bode.strips
		strip_mag.yext = InspectDR.PExtents1D() #Reset
		strip_phase.yext = InspectDR.PExtents1D() #Reset

	#Generate Smith chart:
	plot_smith.data = [] #Clear old data
	β = imag(γline(fmax))
	Θ = β * ℓ; nrev = Θ/pi #revs around the chart
	npts = max(round(Int, 20*nrev), 200) #Ensure sufficient # of points
	f = collect(range(fmin, stop=fmax, length=npts))
	jω = j*2pi*f
	y = Γ(Zline(ZC, γline(f, α=fdepα(f, αDC)), ZL=Zload(jω,R,L,C), ℓ=ℓ))
	wfrm = add(plot_smith, f, y)
	wfrm.line = ldata

	#Clear old annotation/markers:
	plot_bode.userannot = []

	return mplot
end

end #module


text:
module Solutions

    using HDF5
    using SharedArrays

    using Base: TwicePrecision

    using ..CommonFunctions
    using ..Config
    using ..Equations

    export DEFAULT_NSAVE, DEFAULT_NWRITE

    const DEFAULT_NSAVE = 1
    const DEFAULT_NWRITE = 0

    export SolutionVector

    SolutionVector{DT} = Union{Vector{DT}, Vector{TwicePrecision{DT}}}

    export get_data!, set_data!
    export DataSeries, PDataSeries, SDataSeries

    include("solutions/dataseries.jl")

    export TimeSeries, compute_timeseries!

    include("solutions/timeseries.jl")

    export StochasticDataSeries, SStochasticDataSeries

    include("solutions/stochasticdataseries.jl")

    export SemiMartingale
    export WienerProcess, generate_wienerprocess!

    include("solutions/wienerprocess.jl")

    export Solution, StochasticSolution

    include("solutions/solution.jl")
    include("solutions/solutions_common.jl")

    export AtomicSolution, AtomicSolutionODE, AtomicSolutionPODE,
           AtomicSolutionDAE, AtomicSolutionPDAE
    export update!, cut_periodic_solution!

    include("solutions/atomic_solution.jl")
    include("solutions/atomic_solution_ode.jl")
    include("solutions/atomic_solution_pode.jl")
    include("solutions/atomic_solution_dae.jl")
    include("solutions/atomic_solution_pdae.jl")

    export SolutionODE, SolutionPODE, SolutionDAE, SolutionPDAE, SolutionSDE, SolutionPSDE
    export PSolutionPDAE, SSolutionPDAE
    export get_initial_conditions, get_initial_conditions!, set_initial_conditions!,
           create_hdf5

    include("solutions/solution_ode.jl")
    include("solutions/solution_pode.jl")
    include("solutions/solution_dae.jl")
    include("solutions/solution_pdae.jl")
    include("solutions/solution_sde.jl")
    include("solutions/solution_psde.jl")

end


text:
module Nullables

include("nullable.jl")

export Nullable, NullException, isnull, unsafe_get

end # module


text:
function senddata(system::CVSystem,n::Int64,sender::Socket)
    # convert data to scientific notation to maintain precision
    a,b=sc(system.branches.Q[14][11,n+1]);
    c,d=sc(system.branches.P[14][11,n+1]);
    e,f=sc(system.liver.Q[n+1,1]);
    g,h=sc(system.liver.P[n+1,1]);

    # write data to IO buffer
    io = IOBuffer();
    write(io,string(a),string(" "),string(b),string(" "),string(c),string(" "),
        string(d),string(" "),string(e),string(" "),string(f),string(" "),
        string(g),string(" "),string(h),string(" "));

    # send data, wait for confirmation of receipt
    ZMQ.send(sender, Message(io))
    msg = ZMQ.recv(sender);

    # print confirmation every 100 time steps
    if n % 100 == 0
        str = unsafe_string(msg);
        print("Received ")
        print(str)
        print(", length: ")
        println(length(str))
        println("Last character: ")
        print("Converted message to float: ")
        println(parse(Float64, str[1:end-1]))
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
mutable struct Mag
    xy::Complex
    z::Real
end
Base.show(io::IO,M::Mag) = print(io, "Mag(xy = ",round(M.xy,digits=2),", z = ",round(M.z,digits=2),")")
# Contructor
Mag(p::Phantom, dir::Symbol) = dir==:x ? Mag.(p.ρ,0) : Mag.(0,p.ρ)
# Arithmetic operations
+(M1::Mag, M2::Mag) = Mag(M1.xy + M2.xy, M1.z + M2.z) #Vector sum
*(α::Float64, M::Mag) = Mag(α*M.xy, α*M.z) #This was defined to easily define the average M_mean = sum 1/N* M
*(M::Mag, α::Float64) = Mag(α*M.xy, α*M.z)
# Other operations
angle(M::Mag) = angle(M.xy)
abs(M::Mag) = abs(M.xy)
getproperty(x::Vector{Mag}, f::Symbol) = getproperty.(x,f)
# Rotation
@doc raw"""
Spinor (\alpha, \beta) × Magnetization (Mx + i My, Mz)

A vector M = (Mx,My,Mz) can be expressed as a complex 2x2 matrix

M = [Mz Mxy⋆;

------- Mxy  -Mz],	with Mxy = Mx + i My.

Then, to operate with a Spinor M+=RMR⋆, or (α,β)×(Mxy,Mz) = (Mxy+,Mz+) with

Mxy+ = 2α⋆βMz+(α⋆)²Mxy-β²Mxy⋆

and

Mz+ = (|α|² - |β|²)Mz-α⋆ β⋆ Mxy-αβMxy⋆ .
"""
*(s::Spinor, M::Mag) = begin
    # Matricial representation:
    # M_mat = [M.z conj(M.xy); M.xy -M.z]
    # Q = [s.α -conj(s.β); s.β conj(s.α)]
    # M⁺ = Q * M_mat * Q'
    # Mag(M⁺[2,1], real(M⁺[1,1]))
	Mag(
        2*conj(s.α)*s.β*M.z+conj(s.α)^2*M.xy-s.β^2*conj(M.xy),
        (abs(s.α)^2-abs(s.β)^2)*M.z-2real(s.α*s.β*conj(M.xy))#2real(s.α*s.β*conj(M.xy))=conj(s.α)*conj(s.β)*M.xy+s.α*s.β*conj(M.xy)
        )
end
#Operation on vector
# M0 = Mag(0,1)
# Mf = Ry(π/2)*M0
# Mf.xy, Mf.z


text:
include("basic_omnivory_module.jl")
include("pulse.jl")
using PyPlot
pygui(true)

# PULSE SIMULATIONS
# returns dynamics and metrics in different panels, we actually dispatch the
# different pannels between figure 3 and 4 with Inkscape (we did the same
# for figures S3 and S4).

# we pass "true" to this ARGS for supplementary figures
if isempty(ARGS) || ARGS[1] != "true"
    # default for main figures 3 and 4
    par = ModelPar()
    fl_fig = "fig/fig_pulse.svg"
else
    # parameters for SI figures S3 and S4
    par = ModelPar(e_CP = 0.5, e_RC = 0.5, e_RP = 0.5)
    fl_fig = "fig/fig_pulse_si.svg"
end

# run pulse simulation
res = pulse(par, 0.1, 2.0, 2.0)

# temporal ranges of the different phases
rgs = [
    180:0.01:199,
    200:0.01:205,
    205:0.01:325,
    325:0.01:350
    ]
# check degree of omnivory for the different phases
check_doo(res, rgs)

# plot helper
function plot_illustration(sol, eq, ttl_id, leg = true, y_max = 5)
    RCP_cols = ["#1f77b4", "#ff7f0e", "#2ca02c"]
    labs = ["R", "C", "P"]
    ttls = ["Food Chain", "Passive Omnivory", "Active Omnivory"]

    for i in 1:3
        hlines(eq[i], sol.t[1], sol.t[end], color = RCP_cols[i], alpha = 0.5)
    end

    for i in 1:3
        plot(sol.t, sol[i, :], color = RCP_cols[i], label = labs[i])
    end

    if leg
        legend()
    end

    xlim(sol.t[1], sol.t[end],)
    ylim(0, y_max)
    title(ttls[ttl_id])
    ylabel("Density")
end

# Layout
fig = figure(figsize = (8, 9))
st_t = 175
en_t = 350
rg_t = st_t:0.1:en_t
# Dynamics
## FC
subplot(3, 2, 1)
plot_illustration(
    res[1][14](rg_t),
    res[1][12],
    1
)
## PO
subplot(3, 2, 3)
plot_illustration(
    res[2][14](rg_t),
    res[2][12],
    2
)
## RO
subplot(3, 2, 5)
plot_illustration(
    res[3][14](rg_t),
    res[3][12],
    3
)
xlabel("Time")
# Barplots (second Column)
## Stability
subplot(3, 2, 2)
x_loc = [0, 1, 2]
x_labels = ["FC", "PO", "RO"]
plt.bar(x_loc, [res[i][2] for i in 1:3])
plt.xticks(x_loc, x_labels)
plt.tick_params(axis = "x", which = "both", length = 0)
ylabel("Local Return Time")
## Overshoot
subplot(3, 2, 4)
### Set position of bar on X axis
### set width of bar
bar_width = 0.25
r1 = 1:3
r2 = [x + bar_width for x in r1]
r3 = [x + bar_width for x in r2]
### Make the plot
plt.bar(r1, [res[i][7] for i in 1:3], width = bar_width, edgecolor = "white", label = "R")
plt.bar(r2, [res[i][6] for i in 1:3], width = bar_width, edgecolor = "white", label = "C")
plt.bar(r3, [res[i][5] for i in 1:3], width = bar_width, edgecolor = "white", label = "P")
plt.xticks(r2, x_labels)
plt.tick_params(axis = "x", which = "both", length = 0)
ylabel("Degree of Overshoot")
plt.legend()
## Range (max-min)
subplot(3, 2, 6)
plt.bar(r1, [res[i][10] for i in 1:3], width = bar_width, edgecolor = "white", label = "R")
plt.bar(r2,  [res[i][9] for i in 1:3], width = bar_width, edgecolor = "white", label = "C")
plt.bar(r3,  [res[i][8] for i in 1:3], width = bar_width, edgecolor = "white", label = "P")
plt.xticks(r2, x_labels)
plt.tick_params(axis = "x", which = "both", length = 0)
ylabel("Max - Min")

plt.legend()

tight_layout()
savefig(fl_fig)

max_stars_repo_name: McCannLab/Labmnivory

max_issues_repo_name: McCannLab/Labmnivory

max_forks_repo_name: McCannLab/Labmnivory


text:
__precompile__()

"""
Implements the notation `..` for indexing arrays. It's similar to the Python
`...` in that it means 'all of the columns before (or after)'.

`..` slurps dimensions greedily, meaning that the first occurrence
of `..` in an index expression creates as many slices as possible. Other
instances of `..` afterwards are treated simply as slices. Usually, you
should only use one instance of `..` in an indexing expression to avoid
possible confusion.

# Example

```jldoctest
julia> A = Array{Int}(undef,2,4,2);

julia> A[..,1] = [2 1 4 5
                  2 2 3 6];

julia> A[..,2] = [3 2 6 5
                  3 2 6 6];

julia> A[:,:,1] == [2 1 4 5
                    2 2 3 6]
true

julia> A[1,..] = reshape([3 4
                          5 6
                          4 5
                          6 7],1,4,2) # drops singleton dimension
...

julia> B = [3 4
            5 6
            4 5
            6 7];

julia> B == reshape(A[1,..],4,2)
true
```
"""
module EllipsisNotation

using ArrayInterface
using ArrayInterface: indices

import Base: to_indices, tail

struct Ellipsis end
const ..   = Ellipsis()

@inline function to_indices(A, inds::NTuple{M, Any}, I::Tuple{Ellipsis, Vararg{Any, N}}) where {M,N}
    # Align the remaining indices to the tail of the `inds`
    colons = ntuple(n->Colon(), M-N)
    to_indices(A, inds, (colons..., tail(I)...))
end

Base.@propagate_inbounds function ArrayInterface.to_indices(A, inds::Tuple{Vararg{Any,M}}, I::Tuple{Ellipsis,Vararg{Any, N}}) where {M,N}
    return ArrayInterface.to_indices(A, inds, (ntuple(i -> indices(inds[i]), Val(M-N))..., tail(I)...))
end
ArrayInterface.to_indices(A, inds::Tuple{}, I::Tuple{Ellipsis}) = ()
ArrayInterface.is_linear_indexing(A, args::Tuple{Arg}) where {Arg<:Ellipsis} = false

#=
ArrayInterface.can_flatten(::Type{A}, ::Type{T}) where {A,T<:Ellipsis} = true
@inline function ArrayInterface.flatten_args(A, args::Tuple{Arg,Vararg{Any,N}}) where {Arg<:Ellipsis,N}
    return (ntuple(i -> indices(axes(A, i)), Val(ndims(A) - N))..., flatten_args(A, tail(args))...)
end
=#

export ..

end # module


text:
using BenchmarkTools
using WaveSim

function main()
  sim_params = WaveSimParameters(
    focus_depth = 0.03,  # [m]
    steer_angle = 10.0,  # [deg]
    aperture_size = 0.01,  # [m]
    temporal_res = 0.5e-6,  # [s]
    spatial_res = [32, 64]  # [pixels]
  )

  # Compute focusing delays for the elements of the phased array.
  trans_delays = WaveSim.delays_from_focus_and_steer(sim_params)

  # Run the simulation.
  images = WaveSim.wavesim(trans_delays, sim_params)

  return images, sim_params
end

display(@benchmark main())


text:
export TTP, MTT, CBV, CBF, CBF_maxgrad, MTT_fromCBF, lowpasshann, hann, slidingwindowavg

"""
    hann(L::Int)
Hann window of length L, returns array of length L with entries from [0,1]

"""
function hann(L::Int)
  v = zeros(L);
  for i = 1:L
    v[i] = 0.5*(1-cos(2*pi*(i)/(L-1)))
  end
  return v
end

"""
    lowpasshann(tseries[;windowsize=12])
Applies low pass Hann filter pixelwise on the Fourier transformed temporal data. Returns back tranformed array with data
"""
function lowpasshann(tseries::Array;windowsize=12)

  (px,py,pz,pt) = size(tseries)

  hant=zeros(1,1,1,pt)

  hant[1,1,1,round(Int,pt/2)-round(Int,round(Int,windowsize))+1:round(Int,pt/2)]=hann(round(Int,round(Int,windowsize)))
  hant[1,1,1,round(Int,pt/2)+1:round(Int,pt/2)+round(Int,round(Int,windowsize))]=hann(round(Int,round(Int,windowsize)))
  hant[1,1,1,round(Int,pt/2)-round(Int,round(Int,windowsize/2))+1:round(Int,pt/2)+round(Int,round(Int,windowsize/2))]=ones(round(Int,round(Int,windowsize)),1)

  tdata=fft(tseries,4)
  tdata=abs.((ifft(broadcast(*,tdata,ifftshift((hant),4)),4)))

  return tdata

end

function lowpasshann(tseries::ImageMeta;windowsize=12)
  odata = arraydata(tseries)
  return lowpasshann(odata,windowsize=12)
end

"""
    slidingwindowavg(tdata::Array[;filterWidth=20])
Applies a sliding window average on the temporal dimension (fourth dimension)
"""
function slidingwindowavg(tdata::Array;filterWidth=20)
	if (filterWidth > size(tdata,4))
	   filterWidth=floor(Int,size(tdata,4)/2)
           @warn "filterWidth was chosen too large. Reduced to floor(Int,size(tdata,4)/2"
	end
  ldata=tdata;
  for t=1:size(tdata,4)
      if t<size(tdata,4)-filterWidth
          ldata[:,:,:,t]=mean(tdata[:,:,:,t:t+filterWidth],dims=4)
      else
          ldata[:,:,:,t]=mean(tdata[:,:,:,end-filterWidth:end],dims=4)
      end
  end
  return ldata
end
"""
    CBF_maxgrad(tseries::ImageMeta[;positionArtery=[1,1,1], alpha=0.4, alpha2=2,windowsize=6])
Calculates the maximum gradient pixelwise along temporal dimension for pixels inside the mask and divides by maximum concentration at the position of the artery. Alternative: CBF which calculates the cerebral blood flow from the ratio CBV/MTT
"""
function CBF_maxgrad(tseries::ImageMeta;positionArtery=[1,1,1], alpha=0.4, alpha2=2,windowsize=6)
  return CBF_maxgrad(tseries.data;positionArtery=positionArtery, alpha=alpha, alpha2=alpha2,windowsize=windowsize)
end

"""
    CBF_maxgrad(tseries::Array[;positionArtery=[1,1,1], alpha=0.4, alpha2=2,windowsize=6])
Calculates the maximum gradient pixelwise along temporal dimension for pixels inside the mask and divides by maximum concentration at the position of the artery. Alternative: CBF which calculates the cerebral blood flow from the ratio CBV/MTT
"""
function CBF_maxgrad(tseries::Array;positionArtery=[1,1,1], alpha=0.4, alpha2=2,windowsize=6)
  mask = generateMaskFromMIP(tseries[:,:,:,:],alpha, alpha2)
  tmaxgrad, maxgrad = maximumGradient(tseries, mask)

  tseries = lowpasshann(tseries,windowsize=windowsize)
  println(positionArtery)
  CBF = zeros(size(maxgrad))

  for i in CartesianIndices(maxgrad)
    if (mask[i]>0)
      CBF[i] = maxgrad[i]/maximum(tseries[positionArtery[1],positionArtery[2],positionArtery[3],:])
    else
      CBF[i] = 0
    end
  end
  return CBF

end
"""
    CBV(tseries[; positionArtery=[1,1,1], alpha = 0.4, alpha2=0.4,windowsize=12])
Calculates the pixelwise sum along the temporal dimension and normalizes to sum at position of artery
"""
function CBV(tseries; positionArtery=[1,1,1], alpha = 0.4, alpha2=0.4,windowsize=12)
  mask = generateMaskFromMIP(tseries[:,:,:,:],alpha, alpha2)
  tseries = lowpasshann(tseries,windowsize=windowsize)
  dimension = length(size(tseries))
  return sum(tseries,dims=dimension)/sum(vec(tseries[positionArtery[1],positionArtery[2],positionArtery[3],:])).*mask
end

"""
    firstMoment(tseries::Array)
Calculates the first moment (sum(t*c(t))) pixelwise along temporal dimension (4)
"""
function firstMoment(tseries::Array)

  (px,py,pz,pt) = size(tseries)
  time = 1:pt
   firstMom = zeros(px,py,pz,1)
   for ib = 1:px
     for jb = 1:py
       for kb = 1:pz
	       firstMom[ib,jb,kb,1]=sum(tseries[ib,jb,kb,:].*time)
       end
     end
   end

  return firstMom
end
"""
    generateMaskFromMIP(tseries,alpha,alpha2[;windowsize=12])
Not all pixels in the DF FoV experience signal uptake. Mask covers only regions with signal increase during time
"""
function generateMaskFromMIP(tseries,alpha,alpha2;windowsize=12)
  tseries = lowpasshann(tseries,windowsize=windowsize)
  dataMIP = maximum(tseries,dims=4)
  mask = zeros(size(dataMIP))
  mask2 = zeros(size(dataMIP))
  mask2[dataMIP .> maximum(dataMIP)*alpha2] .= 1
  mask[dataMIP .> mean(tseries[:,:,:,1:5],dims=4)*alpha ] .= 1

  return mask.*mask2
end
"""
    MTT(tseries::ImageMeta[; alpha = 0.4, alpha2 =0.4])
Calculates the mean transit time pixelwise from the first moment MTT =(sum(c(t)*t))/sum(c(t)) along the temporal dimension(4)
"""
function MTT(tseries::ImageMeta; alpha = 0.4, alpha2 =0.4,DFPeriodInS=0.02154)
  return MTT(tseries.data; alpha = alpha, alpha2 =alpha2,DFPeriodInS=0.02154)
end
"""
    MTT(tseries::Array[; alpha = 0.4, alpha2 =0.4])
Calculates the mean transit time pixelwise from the first moment MTT =(sum(c(t)*t))/sum(c(t)) along the temporal dimension(4)
"""
function MTT(tseries::Array; alpha = 0.4, alpha2 =0.4,DFPeriodInS=0.02154)

  mask = generateMaskFromMIP(tseries,alpha, alpha2)
  firstMom = firstMoment(tseries)
  IntegratedConcentration = sum(tseries,dims=4)
  mtt = firstMom.*mask./IntegratedConcentration
  for i in CartesianIndices(mtt)
	  if (IntegratedConcentration[i] == 0)
                mtt[i] = 0
          end
  end
  return mtt.*DFPeriodInS
end

"""
    TTP(tseries::ImageMeta[;alpha=0.4,alpha2=0.4,windowsize=12])
Calculates the time in number of repetition times (Multiply with T_R of method to get time in s)
"""
function TTP(tseries::ImageMeta; alpha=0.4,alpha2=0.4,windowsize=12,DFPeriodInS=0.02154)
  dataFG = maximum(tseries, dims=4)
  dataFG.data = TTP(tseries.data,alpha=alpha,alpha2=alpha2,windowsize=windowsize)
  return dataFG
end

function TTP(tseries::Array; alpha=0.4,alpha2=0.4,windowsize=12,DFPeriodInS=0.02154)

  dataFG = maximum(tseries, dims=4)
  dataMIP = deepcopy(dataFG)

  ldata = lowpasshann(tseries,windowsize=windowsize)
  #ldata = slidingwindowavg(ldata)
  (px,py,pz,pt) = size(tseries)

   maxdat = zeros(px,py,pz,1)
   ttp = zeros(px,py,pz,1)
   for ib = 1:px
     for jb = 1:py
       for kb = 1:pz
         buf = vec(ldata[ib,jb,kb,:])
         (maxdat[ib,jb,kb,1],ttp[ib,jb,kb,1])=findmax(buf)
       end
     end
   end
  mask=generateMaskFromMIP(tseries[:,:,:,:],alpha,alpha2)
  dataFG = float(ttp).*mask
  return dataFG.*DFPeriodInS

end

function findFirstEntrySmallerThreshHigherInt(tseries,thresh,minInt)
  comparison = tseries .< thresh
  comparison[1:minInt] .= false
  positionInArray = findmax(comparison)[2]
  if (findmax(comparison)[1] == 0)
    positionInArray = length(tseries)
  end
  return positionInArray
end

function maximumGradient(data,mask)
  maxgrad,tmaxgrad=maximumGradient(data)
  return tmaxgrad.*mask,maxgrad.*mask
end

function maximumGradient(data)
  ldata = lowpasshann(data,windowsize=6)

  (px,py,pz,pt) = size(data)
  normfactors = mean(ldata[:,:,:,:],dims=4)
  normdata = broadcast(-,ldata,normfactors)

  tempdiff = zeros(size(normdata)[1],size(normdata)[2],size(normdata)[3],size(normdata)[4]-1)
  for i = 1:size(tempdiff)[4]
    tempdiff[:,:,:,i] = normdata[:,:,:,i+1]-normdata[:,:,:,i]
  end

  tmaxgrad = zeros(px,py,pz,1)
  maxgrad= zeros(px,py,pz,1)
  for ib = 1:px
    for jb = 1:py
      for kb = 1:pz
        buf = vec(tempdiff[ib,jb,kb,:])
        (maxgrad[ib,jb,kb,1],tmaxgrad[ib,jb,kb,1])=findmax(buf)
      end
    end
  end

  return maxgrad,tmaxgrad
end

"""
    findBeginningEndingBolus(tseries;alpha=0.4, alpha2=0.4)
Determines the beginning of the bolus (bolust1) as the time of maximum derivative and the end (bolust2) as the first time the concentration falls short of concentration(bolust1)
"""
function findBeginningEndingBolus(tseries;alpha=0.4, alpha2=0.4)

  ldata = lowpasshann(tseries,windowsize=6)

  (px,py,pz,pt) = size(tseries)
  normfactors = mean(ldata[:,:,:,:],dims=4)
  normdata = broadcast(-,ldata,normfactors)
  mask = generateMaskFromMIP(tseries[:,:,:,:],alpha,alpha2)
  maxgrad,tmaxgrad=maximumGradient(tseries)

  bolust1=minimum(tmaxgrad.*mask+(1 .-mask).*maximum(tmaxgrad))

  tBolusEnd = zeros(px,py,pz,1)

  for ib = 1:px
    for jb = 1:py
      for kb = 1:pz
        concentrationthresh = normdata[ib,jb,kb,Int(tmaxgrad[ib,jb,kb,1])]
        tBolusEnd[ib,jb,kb,1] = findFirstEntrySmallerThreshHigherInt(normdata[ib,jb,kb,:],concentrationthresh,Int(tmaxgrad[ib,jb,kb,1]))
      end
    end
  end
  bolust2=maximum(tBolusEnd.*mask)

  return bolust1,bolust2

end

"""
    MTT_fromCBF(tseries::ImageMeta[;positionArtery=[1,1,1], alpha=0.4, alpha2=2])
Calculates the mean transit time (MTT) pixelwise from the ratio CBV/CBF_maxgrad. The resulting values seem to be more reasonable from the function MTT.
"""
function MTT_fromCBF(tseries::ImageMeta;positionArtery=[1,1,1], alpha=0.4, alpha2=2,DFPeriodInS=0.02154)
  return  MTT_fromCBF(tseries.data;positionArtery=positionArtery, alpha=alpha, alpha2=alpha2,DFPeriodInS=0.02154)
end

"""
    MTT_fromCBF(tseries::Array[;positionArtery=[1,1,1], alpha=0.4, alpha2=2])
Calculates the mean transit time (MTT) pixelwise from the ratio CBV/CBF_maxgrad. The resulting values seem to be more reasonable from the function MTT.
"""
function MTT_fromCBF(tseries::Array;positionArtery=[1,1,1], alpha=0.4, alpha2=2,DFPeriodInS=0.02154)
  cbf = CBF_maxgrad(tseries,positionArtery=positionArtery,alpha=alpha,alpha2=alpha2)
  cbv = CBV(tseries, positionArtery=positionArtery, alpha = alpha, alpha2=alpha2)
  mask = generateMaskFromMIP(tseries[:,:,:,:],alpha,alpha2)
  return cbv./((cbf.*mask)+(1 .-mask)).*mask.*DFPeriodInS
end

"""
    CBF(tseries, alpha=0.4, alpha2 = 0.4)
Calculates the cerebral blood flow from the ratio CBV/MTT. Alternative: CBF_maxgrad as the normalized maximum gradient of the bolus. (CBF_maxgrad seems to give more realistic values)
"""
function CBF(tseries; positionArtery=[1,1,1], alpha = 0.4, alpha2=0.4,windowsize=12)
  cbv = CBV(tseries,positionArtery=positionArtery,alpha=alpha,alpha2=alpha2,windowsize=windowsize)
  mtt = MTT(tseries,alpha=alpha,alpha2=alpha2)
  cbf =zeros(size(cbv))
  mask = generateMaskFromMIP(tseries[:,:,:,:],alpha,alpha2)
  cbf = cbv.*mask./(mtt.*mask+(1 .-mask))

  return cbf
end


text:
# I'm not too proud of this code duplication -
# but first of all we want to keep dependencies small, and second,
# the bool in Quaternions.Quaternion is annoying for OpenGL + and an array of
# Quaternions.
# TODO replace this file by just `using Quaternions`
using StaticArrays, GeometryTypes, Random
struct Quaternion{T}
    data::NTuple{4, T}
    Quaternion{T}(x::NTuple{4, Any}) where T = new{T}(T.(x))
    Quaternion{T}(x::NTuple{4, T}) where T = new{T}(x)
end

Base.eltype(::Quaternion{T}) where T = T
Base.eltype(::Type{Quaternion{T}}) where T = T
Base.length(::Type{<: Quaternion}) = 4
Base.length(::Quaternion) = 4

const Quaternionf0 = Quaternion{Float32}
function Base.show(io::IO, q::Quaternion)
    pm(x) = x < 0 ? " - $(-x)" : " + $x"
    print(io, q[4], pm(q[1]), "im", pm(q[2]), "jm", pm(q[3]), "km")
end
Random.rand(mt::MersenneTwister, ::Random.SamplerType{Quaternion}) = rand(mt, Quaternion{Float64})
Random.rand(mt::MersenneTwister, ::Random.SamplerType{Quaternion{T}}) where T = Quaternion(rand(mt, T), rand(mt, T), rand(mt, T), 1.0)

@inline Quaternion{T}(x1, x2, x3, s) where T = Quaternion{T}((x1, x2, x3, s))
@inline Base.convert(T::Type{<: Quaternion}, x::NTuple{4, Any}) = T(x)
function Base.convert(T::Type{Quaternion{T1}}, x::Quaternion{T2}) where {T1, T2}
    T(T2.(x.data))
end
@inline Quaternion(x1, x2, x3, s) = Quaternion(promote(x1, x2, x3, s))
@inline Quaternion(x::NTuple{4, T}) where T = Quaternion{T}(x)
@inline Base.getindex(x::Quaternion, i::Integer) = x.data[i]
Base.isapprox(x::Quaternion, y::Quaternion) = all((x.data .≈ y.data))

function qrotation(axis::StaticVector{3}, theta::Real)
    u = normalize(axis)
    s = sin(theta / 2)
    Quaternion(s * u[1], s * u[2], s * u[3], cos(theta / 2))
end
function Base.broadcast(f, arg1::Quaternion, arg2::Quaternion)
    Quaternion(f.(arg1.data, arg2.data))
end

Base.abs(q::Quaternion) = sqrt(sum(q.data .^ 2))

LinearAlgebra.normalize(q::Quaternion) = q / abs(q)

Base.:(/)(q::Quaternion, x::Real) = Quaternion(q[1] / x, q[2] / x, q[3] / x, q[4] / x)

function Base.:(*)(quat::Quaternion, vec::StaticVector{2, T}) where T
    x3 = quat * Vec(vec[1], vec[2], T(0))
    VT = similar_type(vec, StaticArrays.Size(2,))
    VT(x3[1], x3[2])
end

# function (*)(q::Quaternions.Quaternion{T}, v::Vec{3, T}) where T
#     t = T(2) * cross(Vec(q[1], q[2], q[3]), v)
#     v + q[4] * t + cross(Vec(q[1], q[2], q[3]), t)
# end

function Base.:(*)(quat::Quaternion{T}, vec::StaticVector{3}) where T
    num = quat[1] * T(2)
    num2 = quat[2] * T(2)
    num3 = quat[3] * T(2)
    num4 = quat[1] * num
    num5 = quat[2] * num2
    num6 = quat[3] * num3
    num7 = quat[1] * num2
    num8 = quat[1] * num3
    num9 = quat[2] * num3
    num10 = quat[4] * num
    num11 = quat[4] * num2
    num12 = quat[4] * num3
    VT = similar_type(vec, StaticArrays.Size(3,))
    return VT(
        (1f0 - (num5 + num6)) * vec[1] + (num7 - num12) * vec[2] + (num8 + num11) * vec[3],
        (num7 + num12) * vec[1] + (1f0 - (num4 + num6)) * vec[2] + (num9 - num10) * vec[3],
        (num8 - num11) * vec[1] + (num9 + num10) * vec[2] + (1f0 - (num4 + num5)) * vec[3]
    )
end
Base.conj(q::Quaternion) = Quaternion(-q[1], -q[2], -q[3], q[4])

function Base.:(*)(q::Quaternion, w::Quaternion)
    Quaternion(
        q[4] * w[1] + q[1] * w[4] + q[2] * w[3] - q[3] * w[2],
        q[4] * w[2] - q[1] * w[3] + q[2] * w[4] + q[3] * w[1],
        q[4] * w[3] + q[1] * w[2] - q[2] * w[1] + q[3] * w[4],
        q[4] * w[4] - q[1] * w[1] - q[2] * w[2] - q[3] * w[3],
    )
end

const SMat{N, L} = Mat{N, N, T, L} where T

SMat{N, L}(q::Quaternion{T}) where {N, T, L} = Mat{N, N, T, L}(q)

function Mat4{ET}(q::Quaternion{T}) where {T, ET}
    sx, sy, sz = 2q[4]*q[1],  2q[4]*q[2],   2q[4]*q[3]
    xx, xy, xz = 2q[1]^2,    2q[1]*q[2],  2q[1]*q[3]
    yy, yz, zz = 2q[2]^2,    2q[2]*q[3],  2q[3]^2
    T0, T1 = zero(ET), one(ET)
    Mat{4, 4, ET}(
        T1-(yy+zz), xy+sz,      xz-sy,      T0,
        xy-sz,      T1-(xx+zz), yz+sx,      T0,
        xz+sy,      yz-sx,      T1-(xx+yy), T0,
        T0,         T0,         T0,         T1
    )
end

concrete_type(::Type{Any}, ::Type{T}) where T = T
concrete_type(::Type{T}, x) where T = T

function Mat3{ET}(q::Quaternion{T}) where {T, ET}
    sx, sy, sz = 2q[4]*q[1], 2q[4]*q[2],  2q[4]*q[3]
    xx, xy, xz = 2q[1]^2,   2q[1]*q[2], 2q[1]*q[3]
    yy, yz, zz = 2q[2]^2,   2q[2]*q[3], 2q[3]^2
    T0, T1 = zero(ET), one(ET)
    Mat{3, 3, ET}(
        T1-(yy+zz), xy+sz,      xz-sy,
        xy-sz,      T1-(xx+zz), yz+sx,
        xz+sy,      yz-sx,      T1-(xx+yy)
    )
end

function orthogonal(v::T) where T <: StaticVector{3}
    x, y, z = abs.(v)
    other = x < y ? (x < z ? unit(T, 1) : unit(T, 3)) : (y < z ? unit(T, 2) : unit(T, 3))
    return cross(v, other)
end

function rotation_between(u::StaticVector{3, T}, v::StaticVector{3, T}) where T
    k_cos_theta = dot(u, v)
    k = sqrt((norm(u) ^ 2) * (norm(v) ^ 2))
    if (k_cos_theta / k) ≈ T(-1)
        # 180 degree rotation around any orthogonal vector
        Quaternion(normalize(orthogonal(u))..., T(0))
    else
        normalize(Quaternion(cross(u, v)..., k_cos_theta + k))
    end
end


text:
"""
    SequencerResult
        EOSeg::Dict{Tuple,Any} # list of elongation and orderings (BFS,DFS), one per Segment
        EOAlgScale::Dict{Tuple,Any} # elongation and orderings for the cumulated weighted distances
        D::AbstractMatrix # final distance matrix
        mst::LightGraphs.AbstractGraph # final mst
        η::Real # final elongation
        order::AbstractVector #final ordering from bfs

Type containing the results of a Sequencer run.

"""
struct SequencerResult
    EOSeg::Dict{Tuple,Any} # list of elongation and orderings (BFS,DFS), one per Segment
    EOAlgScale::Dict{Tuple,Any} # elongation and orderings for the cumulated weighted distances
    D::AbstractMatrix # final distance matrix
    mst::LightGraphs.AbstractGraph # final mst
    η::Real # final elongation
    order::AbstractVector #final ordering from bfs
end

"""

    D(r::SequencerResult)

Return the final distance matrix from the Sequencer run.
"""
D(r::SequencerResult) = r.D

"""

    mst(r::SequencerResult)

Return the final minimum spanning tree graph from a Sequencer run.
"""
mst(r::SequencerResult) = r.mst

"""

    elong(r::SequencerResult)

Return the elongation coefficient of the final, weighted graph.
"""
elong(r::SequencerResult) = r.η

"""

    order(r::SequencerResult)

Return the result column indices, as determined by the Sequencer algorithm.
"""
order(r::SequencerResult) = r.order

"Sensibly display a SequencerResult object."
show(io::IO, s::SequencerResult) = write(io, "Sequencer Result: η = $(@sprintf("%.4g", elong(s))), order = $(prettyp(order(s),5))")

"""
    sequence(A::VecOrMat{T};
        scales=nothing,
        metrics=ALL_METRICS,
        grid=nothing,
        ) where {T <: Real}

Analyze the provided `m x n` matrix (or m vectors of vectors n) by applying one or more 1-dimensional statistical metrics to
each column of data, possibly after dividing the data into smaller row sets. Columns are compared pairwise for each combination
of metric and scale to create a n x n distance matrix that is analyzed as a graph using a novel algorithm. Details of the algorithm
are provided in the paper cited below, by D. Baron and B. Ménard.

```julia
sequence(A, metrics=ALL_METRICS, scales=nothing)
```

or equivalently:

```julia
sequence(A)
```

as `metrics=ALL_METRICS` and `scales=nothing` are defaults. If you want to specify only one metric, you must wrap it in a 1-tuple.
e.g. to use only KL Divergence, write:

```julia
sequence(A, metrics=(KLD,), scales=nothing
```

Use the `scales` keyword to specify the number of "chunks" into which the data should be divided, as a tuple, e.g.

```julia
sequence(A, scales=(1,3,5))
```
In the `autoscale` mode (enabled by default as `scales=nothing`) SequencerJ will find its own "best scale" based on
running the Sequencer against a sample of columns (10% for now) and picking the scale that results in the greatest elongation.

A 1-D grid may be provided. The grid - whose deltas figure in the distance calculations - must be non-negative
real numbers (Float16, Float32, or Float64). The grid length must equal the number of rows in A.

```julia
julia> sequence(A; metrics=(WASS1D,L2), grid=collect(0.5:0.5:size(A,1))) # grid must equal the size of A along dim 1
```

The paper that describes the Sequencer algorithm and its applications can be found
on [Arxiv](https://arxiv.org/abs/2006.13948).
"""
function sequence(A::VecOrMat{T};
    scales=nothing, #auto-scale by default
    metrics=ALL_METRICS,
    grid=nothing,
    ) where {T <: Real}

    @assert all(sign.(A) .>= 0) && all(!any(isnan.(A)) && !any(isinf.(A))) "Input data cannot contain negative, NaN or infinite values."

    @info "Sequencing data with
    shape: $(size(A))
    metric(s): $(metrics)
    scale(s): $(scales)"

# We need to ensure that A is, or becomes, oriented so that
# observation data run in columns. Each column is a separate data series.
# Each row typically represents a distinct point in time or space where the
# observed data were captured.
    # take a vector of vectors and cat it into a matrix
    A = A isa Vector ? hcat(A...) : A

    # package the arguments to be consumable as iterators
    tuplify(x::PreMetric) = tuple(x)
    tuplify(x) = tuple(x...)

    # replace zeros and Infs with values that work for the math
    clamp!(A, eps(), typemax(eltype(A)))
    @debug "before ensure" grid
    G = ensuregrid!(A, grid)
    @debug "after ensure" G

    metrics = (metrics !== nothing) ? tuplify(metrics) : ALL_METRICS
    if scales === nothing
        with_logger(NullLogger()) do
            scales = autoscale(A, metrics=metrics, grid=G)
        end
        @info "After autoscaling, best scale = $(scales)"
    end
    scales = tuplify(scales)

    return _sequence(A, scales, metrics, G)
end

"Scales used in the autoscale option"
const FIB = [1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,5168]

"""
Sample columns of A and run the Sequencer algorithm against this subspace to identify the
scale at which elongation is greatest. Scales are members of the Fibonacci series, e.g. 1,2,3,5,8,13...

`metrics` and `grid` take the same options as in `sequence`. See [`sequence`](@ref).

p sets the percentage of columns sampled (default = 0.1).
"""
function autoscale(A::AbstractMatrix; scales=FIB, metrics=ALL_METRICS, grid=nothing, p = 0.1)
    M,N = size(A)
    Ŋ = round(Int, N * p, RoundDown)
    @assert Ŋ > 1 "Matrix A contains too few columns ($(N)) to sample with p = $(p). Use a larger p."
    sampidx = sample(collect(axes(A,2)), (Ŋ,), ordered=true)
    Asamp = view(A, :, sampidx)
    minscale = 1
    maxscale = M ÷ 10 # no fewer than 10 rows = observations
    @assert maxscale > 0 "Matrix A contains too few rows ($(M)) for autoscale. The minimum is 10 rows. It is recommended to set scale=(1,)."
    testscales = filter(i->i < maxscale, FIB)
    bestscale = 1
    max_η = 0
    for s in testscales
        r = _sequence(Asamp, s, metrics, ensuregrid!(A, grid))
        η = elong(r)
        if η > max_η
            max_η = η
            bestscale = s
        end
    end
    return bestscale
end

"""
The full algorithm. For internal use only.
"""
function _sequence(A, scales, metrics, grid)

    # rows M and columns N in A
    M, N = size(A)

    MST_all = []
    η_all = []
    EOSeg = Dict{Tuple,Any}() # list of elongation and orderings (BFS,DFS), one per Segment
    EOAlgScale = Dict{Tuple,Any}() # elongation and orderings for the cumulated weighted distances

    @inbounds for k in metrics
        alg = nothing
        for l in scales

            # summary distance matrix for segments
            Dklms = zeros(N, N)
            # Elongation per chunk
            ηs = []
            # BFS sequence per chunk
            orderings = []
            # split the data A and grid into chunks
            # @show size(grid) l
            S, G = _splitnorm(A, grid, l)

            # Each m row in S contains n segments of data,
            # one for each of n data series
            tt = @elapsed for i in eachindex(S, G)
                # m is a chunk of the input data matrix
                # local grid with the same dim 1 dimension as m
                m, lgrid = S[i], G[i]
                @debug "after splitting" size(m), size(lgrid)
                # Some algorithms can be performed on arbitrary real-valued grids
                @debug "k" k
                if k in (EMD, Energy)
                    @debug "adding local grid to alg" k lgrid
                    alg = k((lgrid, lgrid))
                else
                    alg = k
                end
                @debug "alg" alg

                # All the heavy lifting happens here, in the distance calculations
                Dklm = abs.(pairwise(alg, m; dims=2)) .+ ϵ

                # Measure our per-metric, per-scale, per-segment distance matrix
                MSTklm, _, η, bfso = _measure_dm(Dklm)
                # weight the distance matrix by its elongation factor
                Dklm_e = η .* Dklm
                # ensure we stay on the playing field...
                clamp!(Dklm_e, eps(), typemax(eltype(Dklm_e)))
                # update the summary matrix
                Dklms .+= Dklm_e
                # save it for later...don't run away and let me down
                push!(ηs, η)
                push!(orderings, bfso)
            end

            if length(Dklms) < 1
                @warn("Unable to create distance matrices from the given data.")
                return nothing
            end

            # Perform the same MST-based analysis on the weighted results over all segments
            Dkl = Dklms / sum(ηs)

            # measure the per-metric, per-scale average distance matrix
            # minimum spanning tree, elongation, and the breadth-first search tree
            MSTkl, _, ηkl, BFSkl = _measure_dm(Dkl)

            push!(MST_all, MSTkl) # All MSTs (1 per alg,s combo)
            push!(η_all, ηkl) # All elongations (1 per alg,s combo)

            # Store these as intermediate results
            EOSeg[(alg, l)] = (ηs, orderings)
            EOAlgScale[(alg, l)] = (ηkl, BFSkl)

            @info "$(k) at scale $(l): η = $(@sprintf("%.4g", ηkl)) ($(@sprintf("%.2g", tt))s)"
        end
    end

    # Create a sparse N x N proximity matrix to be filled with MST elongation-weighted edge
    # distances
    Pw = _weighted_p_matrix(N, MST_all, η_all)
    # Invert the proximity matrix to get a final distance matrix
    # w/ zeros on the diagonal
    D = inv.(Pw)
    # final minimum spanning tree for analysis
    mstD, stidx, ηD, bfstD = _measure_dm(D)
    # convert the final sequence from tree form into an ordered vector of vertices.
    order = unroll(bfstD, stidx)
    @info "Final: η = $(@sprintf("%.4g", ηD)) sequence = $(prettyp(order,5))"
    return SequencerResult(EOSeg, EOAlgScale, D, mstD, ηD, order)
end

"""
Evaluate a given distance matrix and return its MST graph, elongation, and optimal ordering of vertices. Each vertex
in the result corresponds to a column of data in D.

Internal use only.
"""
function _measure_dm(D::AbstractMatrix)
    g = _mst(D) # minimum spanning tree for D
    stidx = _leastcentralpt(g) # Least central point of MST
    η = elongation(g, stidx)
    bfst = bfs_tree(g, stidx)
    return g, stidx, η, bfst
end

"""
Ensure that the size of the 1-D grid, if provided, is compatible with the data in `A`.
Create a grid if one was not provided.
"""
function ensuregrid!(A, grid=nothing)::AbstractVector
    if isnothing(grid)
        grid = float(collect(axes(A, 1)))
    elseif !(eltype(grid) isa AbstractFloat)
        grid = float(grid)
    end
    @assert length(grid) == size(A, 1) "Grid length must match dims 1 (row count) of A. Got grid:$(length(grid)) and A:$(size(A, 1))"
    grid
end

"Create a N x N proximity matrix from weighted edges of the given graph(s).
Edge weights are provided as η coefficient(s) >= 1.
For internal use only.
"
function _weighted_p_matrix(N, graphs, ηs)
    P = spzeros(N, N)
# fill the proximity matrix with elongation-weighted reciprocal edge weights, treated as distances.
    @inbounds for idx in eachindex(graphs, ηs)
        g = graphs[idx]
        W = LightGraphs.weights(g)
        η = ηs[idx]
        for e in edges(g)
            i, j = src(e), dst(e)
            # proximity is the inverse of weight as distance
            P[i,j] = P[j,i] += η * inv(W[i,j])
        end
    end
    # All P[i] were elongated using an η. Now divide by Ση to get the elongation-weighted average.
    Pw = P / sum(ηs)
    # Put Inf on the 1,1 diagonal to form a true proximity matrix
    Pw .+= sparse(Matrix(Inf * I, size(Pw)...))
    return Pw
end

"Return the index of the least central node in the given graph."
_leastcentralpt(g) = last(findmin(closeness_centrality(g)))

"Return the half-length or mean of the given path lengths (distances)."
_halflen(paths) = mean(paths)

"Return the half-width or mean count of unique values in the given path lengths."
_halfwidth(paths) = mean(count(v -> v == k, paths) for k in unique(paths)) / 2

"""`julia`

`elongation(g, startidx)`

Returns the ratio of the graph g's half-length (mean of path distances) over the
half-width, defined as the mean count of shortest paths from the center node of a
minimum spanning tree over the graph. This function calls the `dijkstra_shortest_paths`
function in `LightGraphs`.
"""
function elongation(G, startidx)
    spaths = dijkstra_shortest_paths(G, startidx, LightGraphs.DefaultDistance())
    # remove those annoying zeros and Infs that ruin it for everyone else
    hlen = clamp(_halflen(spaths.dists), eps(), floatmax(Float32))
    hwid = clamp(_halfwidth(spaths.dists), eps(), floatmax(Float32))
    return hlen / hwid + 1
end

"Return a weighted graph of the minimum spanning tree of the given distance matrix,
calculated using the Kruskal algorithm."
function _mst(dm::AbstractMatrix)
    # Zeros and Infs mess up the calculations...
    clamp!(dm, ϵ, typemax(eltype(dm)))
    # By forcing an unsymmetric distance matrix (from KLD and others) into a Symmetric one
# we are losing 1/2 (the lower half) of the distances
# TODO Should we construct a separate graph for the lower half as well?
# TODO Ask Dalya about this
    g = SimpleWeightedGraph(issymmetric(dm) ? dm : Symmetric(dm))
# Minimum Spanning Tree using the Kruskal algorithm
    mst = kruskal_mst(g; minimize=true)
    Ii = src.(mst)
    Ji = dst.(mst)
    Vi = weight.(mst)
# Populate a result distance matrix with the MST nodes and edge weights and convert it to a graph
    S = sparse(Ii, Ji, Vi, size(dm)...)
    return issymmetric(S) ? SimpleWeightedGraph(S) : SimpleWeightedGraph(Symmetric(S))
end

"""
Split the given `m x n` matrix and grid into `scale` parts along the column dimension such
that the resulting matrices (chunks) are approximately of dimension `k x n` where `k ≊ m / scale`

Each chunk is normalized to sum to 1 (required for balanced tranpsort algorithms.)
"""
function _splitnorm(A::AbstractMatrix{T}, grid, scale) where {T <: Real}

    @assert scale <= length(A[:,1]) "Scale ($(scale)) cannot be larger than the number of data elements ($(length(A[:,1])))."

    M = size(A, 1)
    # break up each column of A into  chunks
    chunklen = cld(M, scale)
    slices = []
    grids = []

    chunkstrt = 1:chunklen:M # the starting indices for each chunk

    for i in chunkstrt
        ii = i + chunklen - one(i) # the end index for the chunk
        ii = clamp(ii, one(ii), M)  # keepin' it real....
        # @show i ii
        S = A[i:ii, :] # a subset of rows, all columns
        G = grid[i:ii]
        # now normalize each column of the slice
        Σslice = sum(S, dims=1)
        push!(slices, S ./ Σslice)
        push!(grids, G)
    end
    return slices, grids
end

EOSeg = Dict{Tuple,Any}() # Dictionary of per-segment intermediate elongation and sequence results
EOAlgScale = Dict{Tuple,Any}() # elongation and orderings for the cumulated weighted distances


text:
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.

using Test
using FAME

using TimeSeriesEcon

if FAME.chli.lib === nothing

    @error "FAME CHLI library not found"
    exit()

end

@testset "workspaces" begin
    w = Workspace(; a = 1, b = TSeries(2020Q1, randn(10)),
        s = MVTSeries(2020M1, (:q, :p), randn(24, 2)),
        c = Workspace(; alpha = 0.1, beta = 0.8,
            n = Workspace(; s = "Hello World")
        ))
    writefame("data.db", w)
    @test length(listdb("data.db")) == 7

    tmp = readfame("data.db")
    @test tmp isa Workspace && length(tmp) == 7
    @test issubset(keys(tmp), (:a, :b, :c_alpha, :c_beta, :c_n_s, :s_p, :s_q))

    tmp = readfame("data.db", "s?")
    @test tmp isa Workspace && length(tmp) == 2
    @test issubset(keys(tmp), (:s_p, :s_q))

    tmp = readfame("data.db", "s?", prefix = "s")
    @test tmp isa Workspace && length(tmp) == 2
    @test issubset(keys(tmp), (:p, :q))

    tmp = readfame("data.db", "c?", collect = "c")
    @test tmp isa Workspace && length(tmp) == 1
    @test issubset(keys(tmp), (:c,))
    @test tmp.c isa Workspace && length(tmp.c) == 3
    @test issubset(keys(tmp.c), (:alpha, :beta, :n_s))

    tmp = readfame("data.db", collect = ["c" => ["n"], "s"])
    @test tmp isa Workspace && length(tmp) == 4
    @test issubset(keys(tmp), (:a, :b, :c, :s))
    @test tmp.c isa Workspace && length(tmp.c) == 3
    @test issubset(keys(tmp.c), (:alpha, :beta, :n))
    @test tmp.c.n isa Workspace && length(tmp.c.n) == 1
    @test issubset(keys(tmp.c.n), (:s,))
    @test tmp.s isa Workspace && length(tmp.s) == 2
    @test issubset(keys(tmp.s), (:p, :q))

    rm("data.db")
end

@testset "missing" begin
    FAME.init_chli()
    pr = TSeries(2020Q1, randn(Float64, 8))
    pr[2020Q4] = NaN
    nu = TSeries(2020Q1, randn(Float32, 8))
    nu[2020Q4] = NaN
    writefame(workdb(), Workspace(; pr, nu))
    for n in ("pr", "nu")
        let b = IOBuffer()
            fame(b, "disp $n")
            seek(b, 0)
            @test sum(Base.Fix1(occursin, r"20:4\s+NC"), readlines(b)) == 1
        end
    end
end


text:
# Parsing tzfiles references:
# - http://man7.org/linux/man-pages/man5/tzfile.5.html
# - ftp://ftp.iana.org/tz/code/tzfile.5.txt

const TZFILE_MAX = unix2datetime(typemax(Int32))

struct TransitionTimeInfo
    gmtoff::Int32     # tt_gmtoff
    isdst::Int8       # tt_isdst
    abbrindex::UInt8  # tt_abbrind
end

function abbreviation(chars::AbstractArray{UInt8,1}, offset::Integer=1)
    unsafe_string(pointer(chars[offset:end]))
end

"""
    read_tzfile(io::IO, name::AbstractString) -> TimeZone

Read the content of an I/O stream and process it as a
[POSIX tzfile](http://man7.org/linux/man-pages/man5/tzfile.5.html). The returned
`TimeZone` will be given the supplied name `name` unless a `FixedTimeZone` is returned.
"""
function read_tzfile(io::IO, name::AbstractString)
    version, tz = read_tzfile_internal(io, name)
    if version != '\0'
        # Another even better transition table after this first one
        version, tz = read_tzfile_internal(io, name, version)
    end
    return tz
end

function read_tzfile_internal(io::IO, name::AbstractString, force_version::Char='\0')
    magic = read(io, 4)  # Read the 4 byte magic identifier
    @assert magic == b"TZif" "Magic file identifier \"TZif\" not found."

    # A byte indicating the version of the file's format: '\0', '2', '3'
    version = Char(read(io, UInt8))
    read(io, 15)  # Fifteen bytes reserved for future use

    tzh_ttisgmtcnt = ntoh(read(io, Int32))  # Number of UTC/local indicators
    tzh_ttisstdcnt = ntoh(read(io, Int32))  # Number of standard/wall indicators
    tzh_leapcnt = ntoh(read(io, Int32))  # Number of leap seconds
    tzh_timecnt = ntoh(read(io, Int32))  # Number of transition dates
    tzh_typecnt = ntoh(read(io, Int32))  # Number of TransitionTimeInfos (must be > 0)
    tzh_charcnt = ntoh(read(io, Int32))  # Number of time zone abbreviation characters

    time_type = force_version == '\0' ? Int32 : Int64

    # Transition time that represents negative infinity
    initial_epoch = time_type == Int64 ? -Int64(2)^59 : typemin(Int32)

    transition_times = Array{time_type}(undef, tzh_timecnt)
    for i in eachindex(transition_times)
        transition_times[i] = ntoh(read(io, time_type))
    end
    lindexes = Array{UInt8}(undef, tzh_timecnt)
    for i in eachindex(lindexes)
        lindexes[i] = ntoh(read(io, UInt8)) + 1 # Julia uses 1 indexing
    end
    ttinfo = Array{TransitionTimeInfo}(undef, tzh_typecnt)
    for i in eachindex(ttinfo)
        ttinfo[i] = TransitionTimeInfo(
            ntoh(read(io, Int32)),
            ntoh(read(io, Int8)),
            ntoh(read(io, UInt8)) + 1 # Julia uses 1 indexing
        )
    end
    abbrs = Array{UInt8}(undef, tzh_charcnt)
    for i in eachindex(abbrs)
        abbrs[i] = ntoh(read(io, UInt8))
    end

    # leap seconds (unused)
    leapseconds_time = Array{time_type}(undef, tzh_leapcnt)
    leapseconds_seconds = Array{Int32}(undef, tzh_leapcnt)
    for i in eachindex(leapseconds_time)
        leapseconds_time[i] = ntoh(read(io, time_type))
        leapseconds_seconds[i] = ntoh(read(io, Int32))
    end

    # standard/wall and UTC/local indicators (unused)
    isstd = Array{Int8}(undef, tzh_ttisstdcnt)
    for i in eachindex(isstd)
        isstd[i] = ntoh(read(io, Int8))
    end
    isgmt = Array{Int8}(undef, tzh_ttisgmtcnt)
    for i in eachindex(isgmt)
        isgmt[i] = ntoh(read(io, Int8))
    end

    # POSIX TZ variable string used for transistions after the last ttinfo (unused)
    if force_version != '\0'
        readline(io)
        posix_tz_str = chomp(readline(io))
    end

    # Now build the time zone transitions
    if tzh_timecnt == 0 || (tzh_timecnt == 1 && transition_times[1] == initial_epoch)
        timezone = FixedTimeZone(name, ttinfo[1].gmtoff)
    else
        # Calculate transition info
        transitions = Transition[]
        utc = dst = 0
        for i in eachindex(transition_times)
            info = ttinfo[lindexes[i]]

            # Since the tzfile does not contain the DST offset we need to
            # attempt to calculate it.
            if info.isdst == 0
                utc = info.gmtoff
                dst = 0
            elseif dst == 0
                # isdst == false and the last DST offset was 0:
                # assume that only the DST offset has changed
                dst = info.gmtoff - utc
            else
                # isdst == false and the last DST offset was not 0:
                # assume that only the GMT offset has changed
                utc = info.gmtoff - dst
            end

            # Sometimes tzfiles save on storage by having multiple names in one for example:
            # "WSST\0" at index 1 turns into "WSST" where as index 2 results in "SST"
            # for "Pacific/Apia".
            abbr = abbreviation(abbrs, info.abbrindex)
            tz = FixedTimeZone(abbr, utc, dst)

            if isempty(transitions) || last(transitions).zone != tz
                if transition_times[i] == initial_epoch
                    utc_datetime = typemin(DateTime)
                else
                    utc_datetime = unix2datetime(Int64(transition_times[i]))
                end

                push!(transitions, Transition(utc_datetime, tz))
            end
        end

        # tzfile's only seem to calculate transitions up to TZFILE_MAX even with version
        # that use 64-bit values. Using a heuristic we can apply a reasonable cutoff for
        # time zones that should probably have one.
        #
        # Note: that without knowing that additional transitions do exist beyond the last
        # stored transition we cannot determine with perfect accuracy what the cutoff should
        # be.
        cutoff = nothing
        if DateTime(2037) <= last(transitions).utc_datetime < TZFILE_MAX
            cutoff = TZFILE_MAX
        end

        timezone = VariableTimeZone(name, transitions, cutoff)
    end

    return version, timezone
end


text:
"""
测试包围费米面
"""

using Test
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Surface

function squ()
    qs = common_square_lattice(0.20)
    Γ4 = TFGamma4(qs, 8.0, 16, 20)
    plt = draw_polygon(Γ4.ltris)
    @savepng plt "refine3_1"
    #
    Γ4 = refine_to_surface(Γ4)
    plt = draw_polygon(Γ4.ltris)
    @savepng plt "refine3_2"
    #
    Γ4 = refine_to_surface(Γ4)
    plt = draw_polygon(Γ4.ltris)
    @savepng plt "refine3_3"
end

squ()


text:
"""
   HermiteWeight()

is a quasi-vector representing `exp(-x^2)` on ℝ.
"""
struct HermiteWeight{T} <: Weight{T} end

HermiteWeight() = HermiteWeight{Float64}()
axes(::HermiteWeight{T}) where T = (Inclusion(ℝ),)
function getindex(w::HermiteWeight, x::Number)
    x ∈ axes(w,1) || throw(BoundsError())
    exp(-x^2)
end

sum(::HermiteWeight{T}) where T = sqrt(convert(T, π))

struct Hermite{T} <: OrthogonalPolynomial{T} end
Hermite() = Hermite{Float64}()
orthogonalityweight(::Hermite{T}) where T = HermiteWeight{T}()

==(::Hermite, ::Hermite) = true
axes(::Hermite{T}) where T = (Inclusion(ℝ), oneto(∞))

"""
     hermiteh(n, z)

computes the `n`-th Hermite polynomial, orthogonal with
respec to `exp(-x^2)`, at `z`.
"""
hermiteh(n::Integer, z::Number) = Base.unsafe_getindex(Hermite{typeof(z)}(), z, n+1)

broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), ::HermiteWeight{T}, ::Hermite{V}) where {T,V} = Weighted(Hermite{promote_type(T,V)}())

# H_{n+1} = 2x H_n - 2n H_{n-1}
# 1/2 * H_{n+1} + n H_{n-1} = x H_n
# x*[H_0 H_1 H_2 …] = [H_0 H_1 H_2 …] * [0    1; 1/2  0     2; 1/2   0  3; …]
jacobimatrix(H::Hermite{T}) where T = Tridiagonal(Fill(one(T)/2,∞), Zeros{T}(∞), one(T):∞)
recurrencecoefficients(H::Hermite{T}) where T = Fill{T}(2,∞), Zeros{T}(∞), zero(T):2:∞

massmatrix(::Hermite{T}) where T = Diagonal(sqrt(convert(T,π)) .* convert(T,2) .^ (0:∞) .* factorial.(convert(T,0):∞))

@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:Hermite}}, B::Hermite) = massmatrix(Hermite{promote_type(eltype(A),eltype(B))}())

##########
# Derivatives
##########

@simplify function *(D::Derivative, H::Hermite)
    T = promote_type(eltype(D),eltype(H))
    D = _BandedMatrix((zero(T):2:∞)', ℵ₀, -1,1)
    H*D
end

@simplify function *(D::Derivative, Q::OrthonormalWeighted{<:Any,<:Hermite})
    X = jacobimatrix(Q.P)
    Q * Tridiagonal(-X.ev, X.dv, X.ev)
end


text:
const AbstractDataset = Union{Dataset, DDataset}

isquotenode(::Any) = false
isquotenode(::QuoteNode) = true

ispair(::Any) = false
ispair(x::Expr) = length(x.args) >=1 && x.args[1] == :(=>)

istuple(::Any) = false
istuple(x::Expr) = x.head == :tuple

parse_function_call(args...) = parse_function_call!([], args...)

parse_function_call!(syms, d, x, func, args...) =
    _parse_function_call!(syms, d, helper_namedtuples_replacement(x), func, args...)

function _parse_function_call!(syms, d, x, func, args...)
    if x == :(_)
        push!(syms, x)
        func(d, args...)
    else
        x
    end
end

function _parse_function_call!(syms, d, x::QuoteNode, func, args...)
    push!(syms, x)
    func(d, x, args...)
end

function _parse_function_call!(syms, d, x::Expr, func, args...)
    if x.head == :. && length(x.args) == 2 && isquotenode(x.args[2])
        Expr(x.head, _parse_function_call!(syms, d, x.args[1], func, args...), x.args[2])
    elseif x.head == :call && length(x.args) == 2 && x.args[1] == :^
            x.args[2]
    elseif x.head == :call && x.args[1] == :cols
        push!(syms, x.args[2])
        func(d, x.args[2], args...)
    else
        Expr(x.head, (_parse_function_call!(syms, d, arg, func, args...) for arg in x.args)...)
    end
end

function extract_anonymous_function(x, func, args...; usekey = false)
    syms = Any[]
    key = gensym()
    data = gensym()
    function_call = parse_function_call!(syms, data, usekey ? replace_key(key, x) : x, func, args...)
    anon_func = usekey ? Expr(:(->), Expr(:tuple, key, data), function_call) :
        Expr(:(->), data, function_call)
    return anon_func, unique(syms)
end

function replace_key(key, ex)
    MacroTools.postwalk(ex) do x
        if x == Expr(:., :(_), QuoteNode(:key))
            return key
        else
            return x
        end
    end
end

# From Query: use curly brackets to simplify writing named tuples
function helper_namedtuples_replacement(ex)
    MacroTools.postwalk(ex) do x
        if x isa Expr && x.head in (:braces, :bracescat)
            new_ex = Expr(:tuple, x.args...)
            for (j,field_in_NT) in enumerate(new_ex.args)
                if isa(field_in_NT, Expr) && field_in_NT.head==:(=)
                    new_ex.args[j] = Expr(:(=), field_in_NT.args...)
                elseif isquotenode(field_in_NT)
                    new_ex.args[j] = Expr(:(=), field_in_NT.value, field_in_NT)
                elseif isa(field_in_NT, Expr)
                    new_ex.args[j] = Expr(:(=), Symbol(filter(t -> t != ':', string(field_in_NT))), field_in_NT)
                elseif isa(field_in_NT, Symbol)
                    new_ex.args[j] = Expr(:(=), field_in_NT, field_in_NT)
                end
            end
            return new_ex
        else
            return x
        end
    end
end

replace_keyword(arg) = (@capture arg x_ = y_) ? Expr(:kw, x, y) : arg

replace_keywords(args) = map(replace_keyword, args)

_table(cols::C) where{C<:Columns} =
        IndexedTable{C}(cols, Int[], IndexedTables.Perm[], fill(missing, length(cols)), nothing)
_table(c) = c

distinct_tuple(args...) = Tuple(IterTools.distinct(args))


text:
using HWfuncapp
using Test

@testset "HWfuncapp.jl" begin
	# test that q1 with 15 nodes has an error smaller than 1e-9
	@test q1(15)[:error] < 1e-9

    # test that the integral of function h [-10,0] ≈ 1.46039789878568
	@test (q3(10)[2] ≈ 1.46039789878568)== false
    #@test false  # by default fail
end


text:
#=
Functions to compute distances from points to primatives
=#

using LinearAlgebra: normalize,norm,⋅
"""
    closest_point_to_tri(p, a, b, c)

This method from Ericson, "Real-time collision detection"
2005, doesn't require the normal or cross products.
"""
@fastmath function closest_point_to_tri(p,a,b,c)
    # is point `a` closest?
    ab,ac,ap = b-a,c-a,p-a
    d1,d2 = ab ⋅ ap, ac ⋅ ap
    d1 <= 0 && d2 <= 0 && return a

    # is point `b` closest?
    bp = p-b
    d3,d4 = ab ⋅ bp, ac ⋅ bp
    d3 >= 0 && d4 <= d3 && return b

    # is point `c` closest?
    cp = p-c
    d5,d6 = ab ⋅ cp, ac ⋅ cp
    d6 >= 0 && d5 <= d6 && return c

    # is segment `ab` closest?
    vc = d1*d4 - d3*d2
    if (vc <= 0 && d1 >= 0 && d3 <= 0)
        d1==d3 && @show d1,d3
        return a + ab * d1 * inv(d1 - d3)
    end

    # is segment `ac` closest?
    vb = d5*d2 - d1*d6
    if (vb <= 0 && d2 >= 0 && d6 <= 0)
        return a + ac * d2 * inv(d2 - d6)
    end

    # is segment `bc` closest?
    va = d3*d6 - d5*d4
    if (va <= 0 && d4 >= d3 && d5 >= d6)
        return b + (c - b) * (d4 - d3) * inv(d4 - d3 + d5 - d6)
    end

    # closest is interior to `abc`
    denom = inv(va + vb + vc)
    v,w = vb * denom, vc * denom
    return a + ab * v + ac * w
end
"""
    closest(p,tri::Triangle)

Determine the closest point on triangle `tri` to point `p`.
"""
closest(p,tri::Triangle) = closest_point_to_tri(p,tri.points.data...)

"""
    absolute_distance(p,tri::Triangle)

Determine the absolute distance from point `p` to the closest point on triangle `tri`.
"""
absolute_distance(p,tri::Triangle) = norm(p-closest(p,tri))

"""
    signed_distance(p,tri::Triangle)

Determine the signed distance from point `p` to the plane defined by triangle `tri`.
Note that the sign depends on triangle point ordering and `d=0` only implies `p`
lies on the plane, not that it is within the triangle.
"""
signed_distance(p,tri::Triangle) = signed_distance_to_tri(p,tri.points.data...)
signed_distance_to_tri(p,a,b,c) = normalize(orthogonal_vector(a,b,c))⋅(p-a)

nonzero(mesh::AbstractMesh) = (i for i in mesh if sum(abs2,orthogonal_vector(i.points.data...))>0)
"""
    absolute_distance(p, mesh::AbstractMesh)

Return the minimum absolute distance from point `p` to any Triangle in `mesh`.
"""
absolute_distance(p,mesh::AbstractMesh) = minimum(t->absolute_distance(p,t), nonzero(mesh))

"""
    signed_distance(p, mesh::AbstractMesh)

Return the signed distance to the geometry defined by the union of planes
passing through each `mesh` face.
"""
signed_distance(p,mesh::AbstractMesh) = maximum(t->signed_distance(p,t), nonzero(mesh))

"""
    signed_distance(p,s::HyperSphere)

Return the signed distance from `p` to `s`.
"""
signed_distance(p,s::HyperSphere) = norm(p-origin(s))-radius(s)

distance_to_corner(q) = norm(max.(q, 0.))+min(maximum(q),0.)
"""
    signed_distance(p,r::Rect)

Return the signed distance from `p` to `r`.
"""
function signed_distance(p,r::Rect)
    # Get vector relative to Rect corner
    q = abs.(p-origin(r).-0.5*width(r)).-0.5*width(r)
    return distance_to_corner(q)
end

"""
    signed_distance(p,c::Cylinder)

Return the signed distance from `p` to `c`.
"""
function signed_distance(p,c::Cylinder)
    e = (extremity(c)-origin(c))/2    # center to edge (defines length direction)
    r = p-(extremity(c)+origin(c))/2  # center to point
    e₁ = norm(e)                      # cylinder half-length
    r₁ = abs(e ⋅ r) / e₁              # projected length to point
    r₂ = √(r ⋅ r - r₁^2)              # height to point
    return distance_to_corner(Point(r₁-e₁,r₂-radius(c)))
end

"""
    signed_distance(p,prim)

Fallback signed distance function.
"""
signed_distance(p,prim::GeometryPrimitive) = signed_distance(p,GeometryBasics.mesh(prim))

"""
    absolute_distance(p,prim) = |signed_distance(p,prim)|

Fallback absolute distance function.
"""
absolute_distance(p,prim) = abs(signed_distance(p,prim))


text:
struct TrackedNeuralDSDE{R,Z,M1,M2,P,RE1,RE2,T,A,K} <: DiffEqFlux.NeuralDELayer
    model1::M1
    model2::M2
    p::P
    re1::RE1
    re2::RE2
    tspan::T
    args::A
    kwargs::K
    len::Int
    nfes::Vector{Int}

    function TrackedNeuralDSDE(model1, model2, tspan, regularize, args...; kwargs...)
        return_multiple = get(kwargs, :save_everystep, false) || haskey(kwargs, :saveat)
        p1, re1 = Flux.destructure(model1)
        p2, re2 = Flux.destructure(model2)
        p = vcat(p1, p2)
        new{
            regularize,
            return_multiple,
            typeof(model1),
            typeof(model2),
            typeof(p),
            typeof(re1),
            typeof(re2),
            typeof(tspan),
            typeof(args),
            typeof(kwargs),
        }(
            model1,
            model2,
            p,
            re1,
            re2,
            tspan,
            args,
            kwargs,
            length(p1),
            [0, 0],
        )
    end
end

@fastmath function (n::TrackedNeuralDSDE{false,true})(x, p = n.p; func = (u, t, int) -> 0)
    function dudt_(u, p, t)
        n.nfes[1] += 1
        n.re1(p[1:n.len])(u)
    end
    function g(u, p, t)
        n.nfes[2] += 1
        n.re2(p[n.len+1:end])(u)
    end
    tspan = _convert_tspan(n.tspan, p)
    ff = SDEFunction{false}(dudt_, g, tgrad = DiffEqFlux.basic_tgrad)
    prob = SDEProblem{false}(ff, g, x, tspan, p)
    sol = solve(prob, n.args...; sensealg = SensitivityADPassThrough(), n.kwargs...)
    arr = diffeqsol_to_3dtrackedarray(sol)::TrackedArray{Float32,3,Array{Float32,3}}
    nfe1, nfe2 = n.nfes
    n.nfes .= 0

    return arr, nfe1, nfe2, nothing
end

@fastmath function (n::TrackedNeuralDSDE{false,false})(x, p = n.p; func = (u, t, int) -> 0)
    function dudt_(u, p, t)
        n.nfes[1] += 1
        n.re1(p[1:n.len])(u)
    end
    function g(u, p, t)
        n.nfes[2] += 1
        n.re2(p[n.len+1:end])(u)
    end
    tspan = _convert_tspan(n.tspan, p)
    ff = SDEFunction{false}(dudt_, g, tgrad = DiffEqFlux.basic_tgrad)
    prob = SDEProblem{false}(ff, g, x, tspan, p)
    sol = solve(prob, n.args...; sensealg = SensitivityADPassThrough(), n.kwargs...)
    arr = diffeqsol_to_trackedarray(sol)::TrackedArray{Float32,2,Array{Float32,2}}
    nfe1, nfe2 = n.nfes
    n.nfes .= 0

    return arr, nfe1, nfe2, nothing
end

@fastmath function (n::TrackedNeuralDSDE{true,true})(
    x,
    p = n.p;
    func = (u, t, integrator) -> integrator.EEst * integrator.dt,
)
    function dudt_(u, p, t)
        n.nfes[1] += 1
        n.re1(p[1:n.len])(u)
    end
    function g(u, p, t)
        n.nfes[2] += 1
        n.re2(p[n.len+1:end])(u)
    end
    tspan = _convert_tspan(n.tspan, p)
    sv = SavedValues(eltype(tspan), eltype(p))
    svcb = SavingCallback(func, sv)
    ff = SDEFunction{false}(dudt_, g, tgrad = DiffEqFlux.basic_tgrad)
    prob = SDEProblem{false}(ff, g, x, tspan, p)
    sol = solve(
        prob,
        n.args...;
        sensealg = SensitivityADPassThrough(),
        callback = svcb,
        n.kwargs...,
    )
    arr = diffeqsol_to_3dtrackedarray(sol)::TrackedArray{Float32,3,Array{Float32,3}}
    nfe1, nfe2 = n.nfes
    n.nfes .= 0

    return arr, nfe1, nfe2, sv
end

@fastmath function (n::TrackedNeuralDSDE{true,false})(
    x,
    p = n.p;
    func = (u, t, integrator) -> integrator.EEst * integrator.dt,
)
    function dudt_(u, p, t)
        n.nfes[1] += 1
        n.re1(p[1:n.len])(u)
    end
    function g(u, p, t)
        n.nfes[2] += 1
        n.re2(p[n.len+1:end])(u)
    end
    tspan = _convert_tspan(n.tspan, p)
    sv = SavedValues(eltype(tspan), eltype(p))
    svcb = SavingCallback(func, sv)
    ff = SDEFunction{false}(dudt_, g, tgrad = DiffEqFlux.basic_tgrad)
    prob = SDEProblem{false}(ff, g, x, tspan, p)
    sol = solve(
        prob,
        n.args...;
        sensealg = SensitivityADPassThrough(),
        callback = svcb,
        n.kwargs...,
    )
    arr = diffeqsol_to_trackedarray(sol)::TrackedArray{Float32,2,Array{Float32,2}}
    nfe1, nfe2 = n.nfes
    n.nfes .= 0

    return arr, nfe1, nfe2, sv
end


text:
struct MyType
    a::Int
    b::Int
    MyType(a) = new(a,1)
end

fun(a::MyType) = a.b

struct MyType2
    a::Int
end

b(a) = 1

fun(a::MyType2) = begin; if a.a==1 || typeof(a)==MyType2; b(a); end; end

a = MyType(2)
a = MyType(2)

a2 = MyType2(2)
a2 = MyType2(2)

@code_llvm MyType(2)
@code_typed MyType(2)

@code_llvm MyType2(2)
@code_typed MyType2(2)

@code_llvm fun(a)
@code_typed fun(a)

@code_llvm fun(a2)
@code_typed fun(a2)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
module ScratchUsage
using Pkg, Scratch

# TODO: Replace these with Pkg.API.get_uuid() and Pkg.API.get_version()
# from https://github.com/JuliaLang/Pkg.jl/pull/1947
const my_uuid = Base.UUID("93485645-17f1-6f3b-45bc-419db53815ea")
function get_version()
    project_toml = Pkg.TOML.parsefile(joinpath(dirname(@__DIR__), "Project.toml"))
    return VersionNumber(project_toml["version"])
end
const my_version = get_version()

# This function will create a bevy of spaces here
function touch_scratch()
    # Create an explicitly version-specific space
    private_space = get_scratch!(
        my_uuid,
        string(my_version.major, ".", my_version.minor, ".", my_version.patch),
    )
    touch(joinpath(private_space, string("ScratchUsage-", my_version)))

    # Create a space shared between all instances of the same major version,
    # using the `@get_scratch!` macro which automatically looks up the UUID
    major_space = @get_scratch!(string(my_version.major))
    touch(joinpath(major_space, string("ScratchUsage-", my_version)))

    # Create a global space that is not locked to this package at all
    # since this code is called from Pkg.test we need to pretend we are
    # at top level and not in the temporary Pkg.test project
    project = Base.ACTIVE_PROJECT[]
    Base.ACTIVE_PROJECT[] = "@v#.#"
    try
        global_space = get_scratch!("GlobalSpace")
        touch(joinpath(global_space, string("ScratchUsage-", my_version)))
    finally
        Base.ACTIVE_PROJECT[] = project
    end
end

end # module ScratchUsage


text:
immutable Conv1d <: SuperRes
    filter :: Function
    filter_grad :: Function
    eval_grid :: Array{Float64}
    approx_grid :: Array{Float64}
    psf_grid :: Matrix{Float64}
    x_max :: Float64
    Conv1d(filter, filter_grad, eval_grid, approx_grid, x_max) =
        new(filter,
            filter_grad,
            eval_grid,
            approx_grid,
            [filter(y - x) for x in approx_grid, y in eval_grid],
            x_max)
end

function psi(model :: Conv1d, theta :: Vector{Float64})
    # This function computes the direct problem for a single point theta
    return Float64[model.filter(x - theta[1]) for x in model.eval_grid]
end

function dpsi(model :: Conv1d, theta :: Vector{Float64})
    # This function computes the gradient of psi
    return Float64[-model.filter_grad(x - theta[1]) for x in model.eval_grid]
end

function getStartingPoint(model :: Conv1d, v :: Vector{Float64})
    # This function gets an initial value for the location
    return [model.approx_grid[indmin(model.psf_grid * v)]]
end

parameterBounds(model :: Conv1d) = [0.0], [model.x_max] # Bounds for the parameters
dim(model :: Conv1d) = 1

immutable DynamicConv1d<: DynamicSuperRes
    static :: Conv1d
    v_max :: Float64
    grid_v :: Vector{Float64}
    times :: Vector{Float64}
    DynamicConv1d(static, v_max, τ, K, num_v = 10) = new(static,
                                                          v_max,
                                                          linspace(-v_max, v_max, num_v),
                                                          -τ*K:τ:τ*K)
end
function getStartingPoint(model :: DynamicConv1d, v :: Vector{Float64})
    # This function gets an initial value for the location
    vals = Float64[dot(v,psi(model, [x; vel])) for x = model.static.approx_grid, vel = model.grid_v]
    ind = indmin(vals);
    i,j = ind2sub(vals, ind)
    return [model.static.approx_grid[i]; model.grid_v[j]]
end

parameterBounds(model :: DynamicConv1d) = [0.0; -model.v_max], [model.static.x_max; model.v_max] # Bounds for the parameters


text:
using .Commands
using .TextObjects
using Match

"""
    Execute the given command, and return whether to refresh the displayed line
"""
function execute(s :: LE.MIState, command :: MotionCommand) :: Bool
    log("executing motion command: $(command.motion)")
    buf = buffer(s)
    for iteration in 1:command.r1
        # call the command's function to generate the motion object
        motion = gen_motion(buf, command.motion)
        if is_stationary(motion)
            @match key(command) begin
                'j' => history_next(s, mode(s).hist)
                'k' => history_prev(s, mode(s).hist)
            end
        end
        # execute the motion object
        motion(s)
    end
    return true
end
function execute(s :: LE.MIState, command :: LineOperatorCommand) :: Bool
    for r in 1:command.r1
        buf = buffer(s)
        line_textobject = line(buf)
        line_motion = Motion(line_textobject...)
        op_fn = operator_fn(command.operator)
        eval(Expr(:call, op_fn, buf, line_motion))
    end
    return true
end

function execute(s :: LE.MIState, command :: InsertCommand) :: Bool
    # cmds = [aAiIoO]
    buf = buffer(s)
    motion = @match command.c begin
        'o' => begin
            endd = line_end(buf).stop
            insert(buf, line_end(buf).stop, '\n')
            if position(buf) == endd
                nothing
            else
                down(buf)
            end
        end
        'O' => begin
            insert(buf, line_zero(buf).stop, '\n')
            up(buf)
        end
        _ => gen_motion(buf, command.c)
    end
    if motion isa Motion
        motion(buf)
    end
    trigger_insert_mode(s)
    return true
end

function execute(s :: LE.MIState, command :: OperatorCommand) :: Bool
    # *5*d2w
    log(command)
    @log op_fn = operator_fn(command.operator)
    for r1 in 1:command.r1
        # 5d*2*w
        for r2 in 1:command.r2
            @log r2
            buf = buffer(s)
            motion = gen_motion(buf, command.action)
            @log result = eval(Expr(:call, op_fn, buf, motion))
        end
    end
    if op_fn === :change
        trigger_insert_mode(s)
    end
    true
end


text:
using MetaGraphs
using LightGraphs
using GraphMLReader
using GraphIO
using TimerOutputs

abstract type MyAbstractCity end

"""
Structure City
Instantiate a structure contening the information about the city (POI, Working Place, etc.)

name_file::String : name of the graph file (in graphml format)
n_node::Int32 : number of nodes of the graph
id_poi::Array{Int32,1} : id of all POI
id_aliment::Array{Int32,1} : id of all restaurants
id_ecole::Array{Int32,1} : id of all school
id_work::Array{Int32,1} : id of all working place
area_work::Array{Float64,1} : area of all working place (expressed in m2)
proba_resto::Array{Float64,2} : probability to go to the restaurant for each working place
"""

export City

struct City <: MyAbstractCity

    name_file::String
    n_node::Int32
    id_poi::Array{Int32,1}
    id_aliment::Array{Int32,1}
    id_ecole::Array{Int32,1}
    id_work::Array{Int32,1}
    area_work::Array{Float64,1}
    proba_resto::Array{Float64,2}

    function City()

        # Read the graph file (from GraphMLReader)
        name_file = ParamFile.GRAPH_FILE ;
        graph = loadgraphml("Data/$(lpad(name_file,1)).graphml") ;
        # Compute the number of nodes (nv: LightGraphs)
        n_node = nv(graph) ;
        # Find the id of POI, food and school
        id_poi, id_aliment, id_ecole = poi(graph,n_node) ;
        # find the id of Working Place and the area expressed in m2
        id_work, area_work = working_place(graph,n_node) ;
        # Compute the probability to go to the restaurant for each working place
        if(Journee.CHOIX_ALEATOIRE_RESTO)
            proba_resto = ones(Float64,(length(id_work),length(id_aliment)))/length(id_aliment) ;
        else
            proba_resto = proba_restaurant() ;
        end

        return new(
        name_file,
        n_node,
        id_poi,
        id_aliment,
        id_ecole,
        id_work,
        area_work,
        proba_resto
        )

    end

end

"""
POI(graph, n_node)

    Function used to find the id node of all POI
    Input:
    - graph: undirected Int64 metagraph with Float64 weights
    - n_node ({Int32}): total number of nodes
    Output:
    - id_poi (Array{Int32,.}): Vector containing the id of all POI

"""
function poi(graph::MetaGraph,n_node::Int)

    # Initialisation of the vector
    id_poi = Array{Int32}([]) ;

    # Find all nodes of POI type
    @inbounds for i=1:n_node
        if(get_prop(graph,i,:type)=="POI")
            append!(id_poi,i) ;
        end
    end

    id_aliment = Array{Int32}([]) ;

    @inbounds for i=1:length(id_poi)
        type_poi = get_prop(graph,id_poi[i],:typeNameICI) ;
        if((type_poi == "loisir.restauration") || (type_poi == "commerce.alimentation"))
            append!(id_aliment,i) ;
        end

    end

    id_aliment = id_poi[id_aliment] ;

    id_ecole = Array{Int32}([]) ;

    @inbounds for i=1:length(id_poi)
        type_poi = get_prop(graph,id_poi[i],:typeNameICI) ;
        if((type_poi == "education.primaire") || (type_poi == "education.secondaire") || (type_poi == "education.superieure") || (type_poi == "education.particuliere") )
        #if((type_poi == "education.primaire") || (type_poi == "education.secondaire"))
            append!(id_ecole,i) ;
        end
    end

    id_ecole = id_poi[id_ecole] ;

    return(id_poi,id_aliment, id_ecole)

end

"""
WorkingPlace(graph, n_node)

    Function used to find the id node of all working place
    Input:
    - graph: undirected Int64 metagraph with Float64 weights
    - n_node ({Int32}): total number of nodes
    Output:
    - id_work (Array{Int32,.}): Vector containing the id of all working place
    - area_work (Array{Int32,.}): Vector containing the area of each working place

"""
function working_place(graph::MetaGraph,n_node::Int)

    id_work = Array{Int32}([]) ;

    @inbounds for i=1:n_node
        if(get_prop(graph,i,:type)=="WorkingPlace")
            append!(id_work,i) ;
        end
    end

    n_work = length(id_work) ;
    area_work = rand(Uniform(50,1500),n_work) ;

    return(id_work,area_work)

end

"""
Function proba_restaurant()
Fonction permettant de définir une probabilité liée à la distance entre le restaurant et l'entreprise

"""
function proba_restaurant()

    distance = readdlm("Data/Distance_resto.txt") ;
    nb_resto = size(distance,2) ;
    nb_travail = size(distance,1) ;
    proba_resto = zeros(Float64,(nb_travail,nb_resto)) ;

    @inbounds for i=1:nb_travail
        Ind = findall(distance[i,:] .!= 0) ;
        nb = sum(1/distance[i,Ind]) ;
        proba_resto[i,Ind] = (1/distance[i,Ind]/nb) ;
    end

    return(proba_resto) ;

end


text:
import GridWorlds
import GridWorlds: GW
import Test
import Random
import ReinforcementLearningBase
import ReinforcementLearningBase: RLBase

ENVS = [GW.EmptyRoomDirected,
        GW.EmptyRoomUndirected,
        GW.GridRoomsDirected,
        GW.GridRoomsUndirected,
        GW.SequentialRoomsDirected,
        GW.SequentialRoomsUndirected,
        GW.MazeDirected,
        GW.MazeUndirected,
        GW.GoToTargetDirected,
        GW.GoToTargetUndirected,
        GW.DoorKeyDirected,
        GW.DoorKeyUndirected,
        GW.CollectGemsDirected,
        GW.CollectGemsUndirected,
        GW.DynamicObstaclesDirected,
        GW.DynamicObstaclesUndirected,
        GW.SokobanDirected,
        GW.SokobanUndirected,
        GW.Snake,
        GW.Catcher,
        GW.TransportDirected,
        GW.TransportUndirected,
       ]

const MAX_STEPS = 3000
const NUM_RESETS = 3

get_terminal_returns(env::GW.EmptyRoomDirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.EmptyRoomUndirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.GridRoomsDirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.GridRoomsUndirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.SequentialRoomsDirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.SequentialRoomsUndirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.MazeDirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.MazeUndirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.GoToTargetDirected) = (env.terminal_reward, env.terminal_penalty)
get_terminal_returns(env::GW.GoToTargetUndirected) = (env.terminal_reward, env.terminal_penalty)
get_terminal_returns(env::GW.DoorKeyDirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.DoorKeyUndirected) = (env.terminal_reward,)
get_terminal_returns(env::GW.CollectGemsDirected) = (env.num_gem_init * env.gem_reward,)
get_terminal_returns(env::GW.CollectGemsUndirected) = (env.num_gem_init * env.gem_reward,)
get_terminal_returns(env::GW.DynamicObstaclesDirected) = (env.terminal_reward, env.terminal_penalty)
get_terminal_returns(env::GW.DynamicObstaclesUndirected) = (env.terminal_reward, env.terminal_penalty)
get_terminal_returns(env::GW.SokobanDirected{T}) where {T} = (T(length(env.box_pos)),)
get_terminal_returns(env::GW.SokobanUndirected{T}) where {T} = (T(length(env.box_pos)),)

get_terminal_returns_win(env::GW.Snake{T}) where {T} = GW.get_terminal_reward(env):GW.get_food_reward(env):convert(T, GW.get_terminal_reward(env) + GW.get_height(env)*GW.get_width(env)*GW.get_food_reward(env))
get_terminal_returns_lose(env::GW.Snake{T}) where {T} = GW.get_terminal_penalty(env):GW.get_food_reward(env):convert(T, GW.get_terminal_penalty(env) + GW.get_height(env)*GW.get_width(env)*GW.get_food_reward(env))

get_terminal_returns(env::GW.Catcher) = env.terminal_reward:env.ball_reward:MAX_STEPS*env.ball_reward
get_terminal_returns(env::GW.TransportDirected) = (GW.get_terminal_reward(env),)
get_terminal_returns(env::GW.TransportUndirected) = (GW.get_terminal_reward(env),)

Test.@testset "GridWorlds.jl" begin
    for Env in ENVS
        Test.@testset "$(Env)" begin
            T = Float32
            env = Env(T = T)
            for _ in 1:NUM_RESETS
                RLBase.reset!(env)
                Test.@test RLBase.reward(env) == zero(T)
                Test.@test RLBase.is_terminated(env) == false

                total_reward = zero(T)
                for i in 1:MAX_STEPS
                    action = rand(RLBase.action_space(env))
                    env(action)
                    total_reward += RLBase.reward(env)

                    Test.@test 1 ≤ GW.get_agent_pos(env)[1] ≤ GW.get_height(env)
                    Test.@test 1 ≤ GW.get_agent_pos(env)[2] ≤ GW.get_width(env)

                    if RLBase.is_terminated(env)
                        if Env == GW.Snake
                            Test.@test (total_reward in get_terminal_returns_win(env) || total_reward in get_terminal_returns_lose(env))
                        else
                            Test.@test total_reward in get_terminal_returns(env)
                        end
                        break
                    end

                    if i == MAX_STEPS
                        @info "$Env not terminated after MAX_STEPS = $MAX_STEPS"
                    end
                end
            end
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "ImageMagick"
version = v"6.9.10-12"

# Collection of sources required to build imagemagick
sources = [
    ArchiveSource("https://github.com/ImageMagick/ImageMagick6/archive/6.9.10-12.tar.gz",
    "efaae51489af9f895762bcb7090636f03194daaa026eda97dae230098d2ccec7"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ImageMagick6*/
./configure --prefix=$prefix --host=$target --without-x --disable-openmp --disable-installed --disable-dependency-tracking --without-frozenpaths --without-perl --disable-docs --disable-static
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct(["libMagickWand", "libMagickWand-6.Q16"], :libwand),
    ExecutableProduct("convert", :imagemagick_convert),
    ExecutableProduct("identify", :identify),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("libpng_jll"),
    Dependency("JpegTurbo_jll"),
    Dependency("Libtiff_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
#!/usr/bin/julia

# Date: 21 August 2016
# https://github.com/trizen

# An efficient algorithm for computing factorial of a large number, modulus a larger number.
# Algorithm from: http://stackoverflow.com/questions/9727962/fast-way-to-calculate-n-mod-m-where-m-is-prime

function facmod(n::Int64, mod::Int64)

    f = 1
    for k in n+1:mod-1
       f *= k
       f %= mod
    end

    (-1 * invmod(f, mod) + mod) % mod
end

println(facmod(100_000_000, 100_000_037))

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts