text:
struct UnconventionalAlgorithm_Params{sType, lType, kType}
    S_rhs # rhs of S step (core projected rhs)
    L_rhs # rhs of L step (range projected rhs)
    K_rhs # rhs of K step (corange projected rhs)
    S_kwargs
    L_kwargs
    K_kwargs
    S_alg::sType
    L_alg::lType
    K_alg::kType
end

struct UnconventionalAlgorithm{sType, lType, kType} <: AbstractDLRAlgorithm
    alg_params::UnconventionalAlgorithm_Params{sType, lType, kType}
end
function UnconventionalAlgorithm(; S_rhs = nothing, L_rhs = nothing, K_rhs = nothing,
                                   S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(),
                                   S_alg = Tsit5(), L_alg = Tsit5(), K_alg = Tsit5())
    params = UnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, S_kwargs, L_kwargs, K_kwargs, S_alg, L_alg, K_alg)
    return UnconventionalAlgorithm(params)
end

struct UnconventionalAlgorithm_Cache{uType,SIntegratorType,LIntegratorType,KIntegratorType,yType} <: AbstractDLRAlgorithm_Cache
    US::Matrix{uType}
    VS::Matrix{uType}
    M::Matrix{uType}
    N::Matrix{uType}
    QRK::LinearAlgebra.QRCompactWY{uType, Matrix{uType}}
    QRL::LinearAlgebra.QRCompactWY{uType, Matrix{uType}}
    SIntegrator::SIntegratorType
    LIntegrator::LIntegratorType
    KIntegrator::KIntegratorType
    y
    ycurr::yType
    yprev::yType
    Δy::yType
end

function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1])
    # allocate memory for frequently accessed arrays
    tspan = (t0,t0+dt)

    US = u.U*u.S
    QRK = qr(US)
    M = Matrix(QRK.Q)'*u.U

    VS = u.V*u.S'
    QRL = qr(VS)
    N = Matrix(QRL.Q)'*u.V

    if isnothing(alg.alg_params.K_rhs)
        K_rhs = function (US, V, t)
                    return Matrix(prob.f(TwoFactorRepresentation(US,V),t)*V)
                end
    else
        K_rhs = alg.alg_params.K_rhs
    end
    KProblem = ODEProblem(K_rhs, US, tspan, u.V)
    KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...)

    if isnothing(alg.alg_params.L_rhs)
        L_rhs = function (VS, U, t)
                    return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U)
                end
    else
        L_rhs = alg.alg_params.L_rhs
    end
    LProblem = ODEProblem(L_rhs, VS, tspan, u.U)
    LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...)

    if isnothing(alg.alg_params.S_rhs)
        S_rhs = function (S, (U,V), t)
                    return Matrix(U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V)
                end
    else
        S_rhs = alg.alg_params.S_rhs
    end
    SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, (u.U, u.V))
    SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...)

    return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL,
                                         SIntegrator, LIntegrator, KIntegrator,
                                         nothing, nothing, nothing, nothing)
end

function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm,u,dt; t0 = prob.tspan[1])
    n, r = size(u.U)
    m = size(u.V, 1)
    @unpack y = prob

    yprev = y isa AbstractArray ? deepcopy(y[1]) : y(t0)
    ycurr = deepcopy(yprev)
    Δy = similar(yprev)
    US = zeros(eltype(u.U),n,r)
    QRK = qr(US)
    M = zeros(eltype(u.U),r,r)
    VS = zeros(eltype(u.V), m, r)
    QRL = qr(VS)
    N = zeros(eltype(u.U),r,r)
    KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1)
    LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1)
    SIntegrator = MatrixDataIntegrator(Δy, M*u.S*N', u.U, u.V, 1)

    return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL,
                                         SIntegrator, LIntegrator, KIntegrator,
                                         prob.y, ycurr, yprev, Δy)
end

function init(prob::AbstractDLRProblem, alg::UnconventionalAlgorithm, dt)
    t0, tf = prob.tspan
    @assert tf > t0 "Integration in reverse time direction is not supported"
    u = deepcopy(prob.u0)
    # initialize solution
    sol = init_sol(dt, t0, tf, prob.u0)
    # initialize cache
    cache = alg_cache(prob, alg, u, dt, t0 = t0)
    sol.Y[1] = deepcopy(prob.u0)
    return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0)
end

function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem})
    @unpack y, ycurr, yprev, Δy = cache
    update_data!(ycurr, y, t, dt)
    Δy .= ycurr - yprev
    yprev .= ycurr
    unconventional_step!(u, cache, t, dt)
end

function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem})
    unconventional_step!(u, cache, t, dt)
end

function unconventional_step!(u, cache, t, dt)
    @unpack US, VS, M, N, QRK, QRL, KIntegrator, LIntegrator, SIntegrator = cache

    # K step
    mul!(US, u.U, u.S)
    set_u!(KIntegrator, US)
    step!(KIntegrator, dt, true)
    US .= KIntegrator.u
    QRK = qr!(US)
    mul!(M, Matrix(QRK.Q)', u.U)

    # L step
    mul!(VS, u.V, u.S')
    set_u!(LIntegrator, VS)
    step!(LIntegrator, dt, true)
    VS .= LIntegrator.u
    QRL = qr!(VS)
    mul!(N,Matrix(QRL.Q)',u.V)
    u.V .= Matrix(QRL.Q)
    u.U .= Matrix(QRK.Q)

    set_u!(SIntegrator, M*u.S*N')
    step!(SIntegrator, dt, true)
    u.S .= SIntegrator.u
end

function step!(integrator::DLRIntegrator, ::UnconventionalAlgorithm, dt)
    @unpack u, t, iter, cache, probType = integrator
    unconventional_step!(u, cache, t, dt, probType)
    integrator.t += dt
    integrator.iter += 1
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# Maximal Covering Location Problem (MCLP)
# From paper: "The maximal covering location problem"

"""
	solveMclp!(sim::Simulation;
		numAmbs::Int = sim.numAmbs,
		demandWeights::Dict{Priority,Float} = Dict([p => 1.0 for p in priorities]),
		results::Dict = Dict())
Solves the Maximal Covering Location Problem (MCLP) to locate `numAmbs` ambulances at stations for `sim` and returns the number of ambulances to assign to each station, also the converse is returned - a station index for each ambulance.
The problem assumes that all ambulances are equivalent.

# Keyword arguments
- `numAmbs` is the number of ambulances to solve for, should be in `1:sim.numStations`
- `demandWeights` is the weight to apply to each demand priority for the objective function
- `results` will store results of mclp such as the objective value and decision variable values
"""
function solveMclp!(sim::Simulation;
	numAmbs::Int = sim.numAmbs,
	demandWeights::Dict{Priority,Float} = Dict([p => 1.0 for p in priorities]),
	results::Dict = Dict())

	@assert(1 <= numAmbs <= sim.numStations, "mclp requires that numAmbs be in 1:sim.numStations")
	@assert(sim.travel.numSets == 1) # otherwise need to solve mclp for each travel set?

	# initialise demand and demand coverage data if not already initialised
	sim.demand.initialised || initDemand!(sim)
	@assert(sim.demand.numSets == 1) # otherwise need to solve mclp for each demand set?
	sim.demandCoverage.initialised || initDemandCoverage!(sim)

	# get demand point coverage data
	currentTime = sim.startTime # shorthand
	pointData = Dict{Vector{Int},Float}() # pointData[pointStations] = pointDemand
	for demandPriority in priorities
		if !haskey(demandWeights, demandPriority) || demandWeights[demandPriority] == 0
			continue
		end

		pointsCoverageMode = getPointsCoverageMode!(sim, demandPriority, currentTime)
		demandMode = getDemandMode!(sim.demand, demandPriority, currentTime)
		pointSetsDemands = getPointSetsDemands!(sim, demandPriority, currentTime; pointsCoverageMode = pointsCoverageMode) * demandMode.rasterMultiplier * demandWeights[demandPriority]

		for (i,stationSet) in enumerate(pointsCoverageMode.stationSets)
			get!(pointData, stationSet, 0.0)
			pointData[stationSet] += pointSetsDemands[i]
		end
	end
	pointStations = collect(keys(pointData))
	pointDemands = collect(values(pointData))
	# point j has demand pointDemands[j] and is covered by stations pointStations[j]

	stationsNumAmbs = solveMclp(numAmbs, pointDemands, pointStations; results = results)
	deployment = stationsNumAmbsToDeployment(stationsNumAmbs) # convert to a deployment

	return stationsNumAmbs, deployment
end

"""
	solveMclp(n::Int, pointDemands::Vector{Float}, coverMatrix::Array{Bool,2}; results::Dict = Dict())
Solves the Maximal Covering Location Problem (MCLP) to locate `n` facilities, where `pointDemands[j]` is the demand at point `j` and `coverMatrix[i,j] = true` if facility location `i` can cover point `j`.

`results` will store results of mclp such as the objective value and decision variable values (though the decision variables may have changed from the given problem as duplicate columns in `coverMatrix` are removed).
"""
function solveMclp(n::Int, pointDemands::Vector{Float}, coverMatrix::Array{Bool,2}; results::Dict = Dict())
	l = size(coverMatrix,1) # number of potential facility locations
	@assert(all(x -> x >= 0, pointDemands))
	@assert(0 <= n <= l, "mclp requires that n be between 0 and the number of potential facility locations ($l).")
	@assert(length(pointDemands) == size(coverMatrix,2))

	# combine points that are covered by the same facilities
	pointData = Dict{Vector{Bool}, Float}() # pointData[coverMatrixCol] = pointSetDemand
	for (j, pointDemand) in enumerate(pointDemands)
		coverMatrixCol = coverMatrix[:,j]
		get!(pointData, coverMatrixCol, 0.0)
		pointData[coverMatrixCol] += pointDemand
	end

	# Could also combine any facility locations that cover the same points,
	# or remove facility locations with coverage dominated by another location,
	# but I'll leave that as it is unlikely for my use cases.

	# # remove points that can be covered by more facility locations than the number of facilities that will not be opened,
	# # as it is guaranteed that these points will be covered
	# objValDelta = 0.0
	# for (coverMatrixCol, pointDemand) in pointData
		# if sum(coverMatrixCol) > l - n
			# delete!(pointData, coverMatrixCol)
			# objValDelta += pointDemand
		# end
	# end

	pointSetDemands = collect(values(pointData))
	pointSetFacilities = findall.(collect(keys(pointData))) # pointSetFacilities[i] gives indices of facility locations that cover pointSetDemands[i]

	# # remove point set that cannot be covered
	# i = findfirst(isempty, pointSetFacilities)
	# if i != nothing
		# deleteat!(pointSetFacilities, i)
		# deleteat!(pointSetDemands, i)
	# end

	solution = solveMclp(n, pointSetDemands, pointSetFacilities, results = results)
	solution = vcat(solution, zeros(Int, size(coverMatrix,1) - length(solution))) # right-pad with zeros if needed

	# # change results dict to be for values given from function args instead of modified problem (that will be same size as given, or smaller)
	# results[:objVal] += objValDelta
	# # results[:x] should be unchanged, unless bottom row(s) of coverMatrix are all false
	# # results[:y] will not be for given pointDemands unless all columns of coverMatrix are unique

	return solution
end

"""
	solveMclp(n::Int, pointDemands::Vector{Float}, pointFacilities::Vector{Vector{Int}}; results::Dict = Dict())
Solves the Maximal Covering Location Problem (MCLP) to locate `n` facilities, where `pointDemands[i]` is the demand at point `i` and `pointFacilities[i]` are the indices of facility locations that cover point `i`.

`results` will store results of mclp such as the objective value and decision variable values.
"""
function solveMclp(n::Int, pointDemands::Vector{Float}, pointFacilities::Vector{Vector{Int}}; results::Dict = Dict())
	p = length(pointDemands) # number of demand points
	@assert(all(x -> x >= 0, pointDemands))
	l = maximum([isempty(v) ? 0 : maximum(v) for v in pointFacilities]) # number of potential facility locations
	@assert(0 <= n <= l, "mclp requires that n be between 0 and the number of potential facility locations ($l).")
	@assert(length(pointFacilities) == p)

	model = Model()

	# use cbc solver; solve speed not tested
	set_optimizer(model, with_optimizer(Cbc.Optimizer, logLevel=0))

	@variables(model, begin
		(x[i=1:l], Bin) # x[i] = true if facility location i is used, false otherwise
		(0 <= y[j=1:p] <= 1) # y[j] = true if point j is covered
	end)

	@constraints(model, begin
		(useAllFacilities, sum(x) == n)
		(pointCoverage[j=1:p], y[j] <= sum(x[pointFacilities[j]]))
	end)

	@expressions(model, begin
		(coverage, sum(y .* pointDemands))
	end)

	# solve
	@objective(model, Max, coverage)
	optimize!(model)
	@assert(termination_status(model) == MOI.OPTIMAL)
	xValue = JuMP.value.(x); yValue = JuMP.value.(y) # JuMP and LightXML both export value()
	objVal = JuMP.value(coverage)

	results[:x] = xValue
	results[:y] = yValue
	results[:objVal] = objVal

	# extract solution
	solution = convert(Vector{Int}, round.(xValue))

	if checkMode
		# check constraints
		@assert(sum(solution) == n) # useAllFacilities constraint
	end

	return solution
end

ext: jl


text:
export vortex, symmetric_vortex

import TransportBasedInference: Parallel, Serial, Thread

vortex(X, t::Float64, Ny, Nx, cachevels, config; withfreestream::Bool = false) = vortex(X, t, Ny, Nx, cachevels, config, serial, withfreestream = withfreestream)

function vortex(X, t::Float64, Ny, Nx, cachevels, config, P::Serial; withfreestream::Bool=false)
	Nypx, Ne = size(X)
	@assert Nypx == Ny + Nx "Wrong value of Ny or Nx"

	freestream = Freestream(config.U)

	@inbounds for i = 1:Ne
		col = view(X, Ny+1:Nypx, i)
		source = state_to_lagrange(col, config)

		# Compute the induced velocity (by exploiting the symmetry of the problem)
		reset_velocity!(cachevels, source)
		self_induce_velocity!(cachevels, source, t)

		if withfreestream == true
			induce_velocity!(cachevels, source, freestream, t)
		end

		# Advect the system
		advect!(source, source, cachevels, config.Δt)

		X[Ny+1:Nypx, i] .= lagrange_to_state(source, config)
	end

	return X, t + config.Δt
end

symmetric_vortex(X, t::Float64, Ny, Nx, cachevels, config; withfreestream::Bool=false) = symmetric_vortex(X, t, Ny, Nx, cachevels, config, serial; withfreestream = withfreestream)

function symmetric_vortex(X, t::Float64, Ny, Nx, cachevels, config, P::Serial; withfreestream::Bool=false)
	Nypx, Ne = size(X)
	@assert Nypx == Ny + Nx "Wrong value of Ny or Nx"
	freestream = Freestream(config.U)
	@inbounds for i = 1:Ne
		col = view(X, Ny+1:Nypx, i)
		source = state_to_lagrange(col, config)

		# Compute the induced velocity (by exploiting the symmetry of the problem)
		reset_velocity!(cachevels, source)
		self_induce_velocity!(cachevels[1], source[1], t)
		induce_velocity!(cachevels[1], source[1], source[2], t)

		if withfreestream == true
			induce_velocity!(cachevels[1], source[1], freestream, t)
		end
		# @. cachevels[2] = conj(cachevels[1])

		# We only care about the transport of the top vortices
		# Advect the system
		advect!(source[1:1], source[1:1], cachevels[1:1], config.Δt)

		X[Ny+1:Nypx, i] .= lagrange_to_state(source, config)
	end

	return X, t + config.Δt
end


text:
import Base.AbstractSet
type IntListSet  <: AbstractSet{Int}
    list::Vector{Int} # the current list
    indv::Vector{Int} # the current location indicator (for fast deletion)
    n::Int
    IntListSet(indmax::Integer) = new(Vector{Int}(), zeros(Int,indmax), 0)
end

import Base.eltype, Base.similar
eltype(::Type{IntListSet}) = Int
similar(s::IntListSet) = IntListSet(length(s.indv))

""" Change the position of an element in the list of element. """
function swap_positions!(s::IntListSet, i1::Int, i2::Int)
  v1 = s.list[i1]
  v2 = s.list[i2]
  s.list[i1],s.list[i2] = v2,v1 # swap positions in list
  s.indv[v1] = i2
  s.indv[v2] = i1
  return (v1,v2)
end

""" Remove an element by moving it to the end of the list and removing it. """
function pivot_out!(s::IntListSet, listindex::Int)
  v = swap_positions!(s, listindex, s.n)[1]
  s.indv[v] = 0 # remove it
  s.n -= 1 # reduce the size
end

import Base.push!
function push!(s::IntListSet, v::Int)
  0 < v <= length(s.indv) || throw(
    ArgumentError("IntListSet items must be between [1,$(length(s.indv))]"))
  if s.indv[v] == 0 # this is not here
    s.n += 1
    s.indv[v] = s.n
    if length(s.list) < s.n
      push!(s.list, v)
    else
      s.list[s.n] = v
    end
  end
  return s
end

import Base.delete!
function delete!(s::IntListSet, v::Int)
  0 < v <= length(s.indv) || throw(
    KeyError("IntListSet items must be between [1,$(length(s.indv))]"))
  if s.indv[v] > 0 # this is here
    pivot_out!(s, s.indv[v])
  end
  return s
end

import Base.in
function in(v::Int, s::IntListSet)
  if 0 < v <= length(s.indv)
    @inbounds return s.indv[v] > 0
  else
    return false
  end
end

import Base.pop!
pop!(s::IntListSet) = pop!(s::IntListSet, s.list[s.n])
function pop!(s::IntListSet, n::Integer)
  if n in s
    return (delete!(s, n), n)
  else
    throw(KeyError(n))
  end
end
function pop!(s::IntListSet, n::Integer, default)
  if n in s
    return (delete!(s, n), n)
  else
    return default
  end
end

import Base.empty!
function empty!(s::IntListSet)
  for i in 1:s.n
    s.indv[s.list[i]] = 0
  end
  s.n = 0
  return s
end

# Just iterate over the elements in list order
import Base.start, Base.next, Base.done
start(s::IntListSet) = 1
done(s::IntListSet, i) = i > s.n
next(s::IntListSet, i) = (s.list[i], i+1)

import Base.intersect!
function intersect!(s::IntListSet, iterable)
    front = 1
    for i in iterable
        if i in s
          # okay, but pivot to front of list
          swap_positions!(s, front, s.indv[i])
          front += 1
        else
          delete!(s, i)
        end
    end
    # now remove everything from front:s.n
    for i in front:s.n # this captures the value at the start of the loop
      #ideally, just delete!(s, s.list[i]), but that pivots stuff around
      s.indv[s.list[i]] = 0
    end
    s.n = front - 1
    s
end

import Base.length
length(s::IntListSet) = s.n

import Base.Random.rand
rand(r::AbstractRNG,s::IntListSet) = rand(r, @view s.list[1:s.n])
rand(s::IntListSet) = rand(@view s.list[1:s.n])

import Base.show, Base.show_vector
function show(io::IO, s::IntListSet)
  print(io, "IntListSet(")
  show_vector(io,view(s.list,1:s.n),"[","]")
  print(io, ")")
end

export IntListSet


text:
;; --------------------------------------------------------------------
;;	csv_read.jl
;;
;;						Mar/02/2011
;;
;; --------------------------------------------------------------------
(require 'tables)
(load-file "/var/www/data_base/common/rep_common/text_manipulate.jl")
;; --------------------------------------------------------------------
(format t "*** 開始 ***")
(setq file_in (cadr command-line-args))
(format t file_in)
(setq dict_aa (make-table string-hash string=))
(let ((port_in (open-file file_in 'read)))
	(csv_read_proc port_in dict_aa)
)

(table-walk dict_display_proc_single dict_aa)

(format t "*** 終了 ***")

;; --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using Compat, Dates, DSP, HDF5, Logging, Printf, SeisIO, Test
using SeisIO.Quake, SeisIO.RandSeis, SeisIO.SeisHDF
import Dates: DateTime, Hour, now
import DelimitedFiles: readdlm
import Random: rand, randperm, randstring
import SeisIO: BUF, FDSN_sta_xml,
  auto_coords, code2typ, typ2code,
  bad_chars, checkbuf!, checkbuf_8!, datafields, datareq_summ, endtime,
  fillx_i16_le!, fillx_i32_be!, fillx_i32_le!, findhex, formats, get_http_req,
  get_http_post, get_views, int2tstr, mean, minreq!,
  mktaper!, mktime, parse_charr, parse_chstr, parse_sl,
  read_sacpz!, read_sacpz, read_seed_resp!, read_seed_resp, read_station_xml!,
  safe_isdir, safe_isfile, sμ, t_collapse,
  t_expand, t_win, taper_seg!, tnote, trid, tstr2int, w_time, webhdr,
  xtmerge!, μs,
  diff_x!, int_x!,
  code2resptyp, resptyp2code,
  poly, polyval, polyfit
import SeisIO.RandSeis: getyp2codes, pop_rand_dict!
import SeisIO.Quake: unsafe_convert
import SeisIO.SeisHDF:read_asdf, read_asdf!, id_match, id_to_regex
import Statistics: mean

# ===========================================================================
# All constants needed by tests are here
const path = Base.source_dir()
const unicode_chars = String.(readdlm(path*"/SampleFiles/julia-unicode.csv", '\n')[:,1])
const n_unicode = length(unicode_chars)
const breaking_dict = Dict{String,Any}(
  "0" => rand(Char), "1" => randstring(rand(51:100)),
  "16" => rand(UInt8), "17" => rand(UInt16), "18" => rand(UInt32), "19" => rand(UInt64), "20" => rand(UInt128),
  "32" => rand(Int8), "33" => rand(Int16), "34" => rand(Int32), "35" => rand(Int64), "36" => rand(Int128),
  "48" => rand(Float16), "49" => rand(Float32), "50" => rand(Float64),
  "80" => rand(Complex{UInt8}), "81" => rand(Complex{UInt16}), "82" => rand(Complex{UInt32}), "83" => rand(Complex{UInt64}), "84" => rand(Complex{UInt128}),
  "96" => rand(Complex{Int8}), "97" => rand(Complex{Int16}), "98" => rand(Complex{Int32}), "99" => rand(Complex{Int64}), "100" => rand(Complex{Int128}),
  "112" => rand(Complex{Float16}), "113" => rand(Complex{Float32}), "114" => rand(Complex{Float64}),
  "128" => collect(rand(Char, rand(4:24))), "129" => [randstring(rand(4:24)) for i=1:rand(4:24)],
  "144" => collect(rand(UInt8, rand(4:24))), "145" => collect(rand(UInt16, rand(4:24))), "146" => collect(rand(UInt32, rand(4:24))), "147" => collect(rand(UInt64, rand(4:24))), "148" => collect(rand(UInt128, rand(4:24))),
  "160" => collect(rand(Int8, rand(4:24))), "161" => collect(rand(Int16, rand(4:24))), "162" => collect(rand(Int32, rand(4:24))), "163" => collect(rand(Int64, rand(4:24))), "164" => collect(rand(Int128, rand(4:24))),
  "176" => collect(rand(Float16, rand(4:24))), "177" => collect(rand(Float32, rand(4:24))), "178" => collect(rand(Float64, rand(4:24))), "208" => collect(rand(Complex{UInt8}, rand(4:24))),
  "209" => collect(rand(Complex{UInt16}, rand(4:24))), "210" => collect(rand(Complex{UInt32}, rand(4:24))), "211" => collect(rand(Complex{UInt64}, rand(4:24))), "212" => collect(rand(Complex{UInt128}, rand(4:24))),
  "224" => collect(rand(Complex{Int8}, rand(4:24))), "225" => collect(rand(Complex{Int16}, rand(4:24))), "226" => collect(rand(Complex{Int32}, rand(4:24))), "227" => collect(rand(Complex{Int64}, rand(4:24))), "228" => collect(rand(Complex{Int128}, rand(4:24))),
  "240" => collect(rand(Complex{Float16}, rand(4:24))), "241" => collect(rand(Complex{Float32}, rand(4:24))), "242" => collect(rand(Complex{Float64}, rand(4:24)))
  )
const NOOF = "
  HI ITS CLOVER LOL
          ,-'-,  `---..
         /             \\
         =,             .
  ______<3.  ` ,+,     ,\\`
 ( \\  + `-”.` .; `     `.\\
 (_/   \\    | ((         ) \\
  |_ ;  \"    \\   (        ,’ |\\
  \\    ,- '💦 (,\\_____,’   / “\\
   \\__---+ }._)              |\\
   / _\\__💧”)/                  +
  ( /    💧” \\                  ++_
   \\)    ,“  |)                ++  ++
   💧     “💧  (                 *    +***
"
# ===========================================================================
# All functions used by tests are here
Lx(T::GphysData) = [length(T.x[i]) for i=1:T.n]
change_sep(S::Array{String,1}) = [replace(i, "/" => Base.Filesystem.pathsep()) for i in S]
test_fields_preserved(S1::GphysData, S2::GphysData, x::Int, y::Int) =
  @test(minimum([getfield(S1,f)[x]==getfield(S2,f)[y] for f in datafields]))
test_fields_preserved(S1::SeisChannel, S2::GphysData, y::Int) =
  @test(minimum([getfield(S1,f)==getfield(S2,f)[y] for f in datafields]))

function loop_time(ts::Int64, te::Int64; ti::Int64=86400000000)
  t1 = deepcopy(ts)
  j = 0
  while t1 < te
    j += 1
    t1 = min(ts + ti, te)
    s_str = int2tstr(ts + 1)
    t_str = int2tstr(t1)
    ts += ti
  end
  return j
end

function sizetest(S::GphysData, nt::Int)
  @test ≈(S.n, nt)
  @test ≈(maximum([length(getfield(S,i)) for i in datafields]), nt)
  @test ≈(minimum([length(getfield(S,i)) for i in datafields]), nt)
  return nothing
end

function mktestseis()
  L0 = 30
  L1 = 10
  os = 5
  tt = time()
  t1 = round(Int64, tt/μs)
  t2 = round(Int64, (L0+os)/μs) + t1

  S = SeisData(5)
  S.name = ["Channel 1", "Channel 2", "Channel 3", "Longmire", "September Lobe"]
  S.id = ["XX.TMP01.00.BHZ","XX.TMP01.00.BHN","XX.TMP01.00.BHE","CC.LON..BHZ","UW.SEP..EHZ"]
  S.fs = collect(Main.Base.Iterators.repeated(100.0, S.n))
  S.fs[4] = 20.0
  for i = 1:S.n
    os1 = round(Int64, 1/(S.fs[i]*μs))
    S.x[i] = randn(Int(L0*S.fs[i]))
    S.t[i] = [1 t1+os1; length(S.x[i]) 0]
  end

  T = SeisData(4)
  T.name = ["Channel 6", "Channel 7", "Longmire", "September Lobe"]
  T.id = ["XX.TMP02.00.EHZ","XX.TMP03.00.EHN","CC.LON..BHZ","UW.SEP..EHZ"]
  T.fs = collect(Main.Base.Iterators.repeated(100.0, T.n))
  T.fs[3] = 20.0
  for i = 1:T.n
    T.x[i] = randn(Int(L1*T.fs[i]))
    T.t[i] = [1 t2; length(T.x[i]) 0]
  end
  return (S,T)
end

function remove_low_gain!(S::GphysData)
    # Remove low-gain seismic data channels
    i_low = findall([occursin(r".EL?", S.id[i]) for i=1:S.n])
    if !isempty(i_low)
        for k = length(i_low):-1:1
            @warn(join(["Low-gain, low-fs channel removed: ", S.id[i_low[k]]]))
            S -= S.id[i_low[k]]
        end
    end
    return nothing
end

# Test that data are time synched correctly within a SeisData structure
function sync_test!(S::GphysData)
    local L = [length(S.x[i])/S.fs[i] for i = 1:S.n]
    local t = [S.t[i][1,2] for i = 1:S.n]
    @test maximum(L) - minimum(L) ≤ maximum(2.0./S.fs)
    @test maximum(t) - minimum(t) ≤ maximum(2.0./S.fs)
    return nothing
end

function breaking_seis()
  S = SeisData(randSeisData(), randSeisEvent(), randSeisData(2, c=1.0, s=0.0)[2])

  # Test a channel with every possible dict type
  S.misc[1] = breaking_dict

  # Test a channel with no notes
  S.notes[1] = []

  # Need a channel with a very long name to test in show.jl
  S.name[1] = "The quick brown fox jumped over the lazy dog"

  # Need a channel with a non-ASCII filename
  S.name[2] = "Moominpaskanäköinen"
  S.misc[2]["whoo"] = String[]        # ...and an empty String array in :misc
  S.misc[2]["♃♄♅♆♇"] = rand(3,4,5,6)  # ...and a 4d array in :misc

  #= Here we test true, full Unicode support;
    only 0xff can be a separator in S.notes[2] =#
  S.notes[2] = Array{String,1}(undef,6)
  S.notes[2][1] = String(Char.(0x00:0xfe))
  for i = 2:1:6
    uj = randperm(rand(1:n_unicode))
    S.notes[2][i] = join(unicode_chars[uj])
  end

  # Test short data, loc arrays
  S.loc[1] = GenLoc()
  S.loc[2] = GeoLoc()
  S.loc[3] = UTMLoc()
  S.loc[4] = XYLoc()

  # Responses
  S.resp[1] = GenResp()
  S.resp[2] = PZResp()
  S.resp[3] = MultiStageResp(6)
  S.resp[3].stage[1] = CoeffResp()
  S.resp[3].stage[2] = PZResp()
  S.resp[3].gain[1] = 3.5e15
  S.resp[3].fs[1] = 15.0
  S.resp[3].stage[1].b = randn(Float64, 120)
  S.resp[3].i[1] = "{counts}"
  S.resp[3].o[1] = "m/s"

  S.x[4] = rand(Float64,4)
  S.t[4] = vcat(S.t[4][1:1,:], [4 0])

  # Some IDs that I can search for
  S.id[1] = "UW.VLL..EHZ"
  S.id[2] = "UW.VLM..EHZ"
  S.id[3] = "UW.TDH..EHZ"
  return S
end

function basic_checks(T::GphysData)
  # Basic checks
  for i = 1:T.n
    if T.fs[i] == 0.0
      @test size(T.t[i],1) == length(T.x[i])
    else
      @test T.t[i][end,1] == length(T.x[i])
    end
  end
  return nothing
end

function get_edge_times(S::GphysData)
  ts = [S.t[i][1,2] for i=1:S.n]
  te = copy(ts)
  for i=1:S.n
    if S.fs[i] == 0.0
      te[i] = S.t[i][end,2]
    else
      te[i] += (sum(S.t[i][2:end,2]) + dtμ*length(S.x[i]))
    end
  end
  return ts, te
end

function wait_on_data!(S::GphysData; tmax::Real=60.0)
  τ = 0.0
  t = 10.0
  printstyled(string("      (sleep up to ", tmax + t, " s)\n"), color=:green)
  redirect_stdout(out) do

    # Here we actually wait for data to arrive
    sleep(t)
    τ += t
    while isempty(S)
      if any(isopen.(S.c)) == false
        break
      end
      sleep(t)
      τ += t
      if τ > tmax
        show(S)
        break
      end
    end

    # Close the connection cleanly (write & close are redundant, but
    # write should close it instantly)
    for q = 1:length(S.c)
      if isopen(S.c[q])
        if q == 3
          show(S)
        end
        close(S.c[q])
      end
    end
    sleep(t)
  end

  # Synchronize (the reason we used d0,d1 in our test sessions)
  prune!(S)
  if !isempty(S)
    sync!(S, s="first")
  else
    @warn("No data. Is the server down?")
  end
  return nothing
end

function naive_filt!(C::SeisChannel;
  fl::Float64=1.0,
  fh::Float64=15.0,
  np::Int=4,
  rp::Int=10,
  rs::Int=30,
  rt::String="Bandpass",
  dm::String="Butterworth"
  )

  T = eltype(C.x)
  fe = 0.5 * C.fs
  low = T(fl / fe)
  high = T(fh / fe)

  # response type
  if rt == "Highpass"
    ff = Highpass(fh, fs=fs)
  elseif rt == "Lowpass"
    ff = Lowpass(fl, fs=fs)
  else
    ff = getfield(DSP.Filters, Symbol(rt))(fl, fh, fs=fs)
  end

  # design method
  if dm == "Elliptic"
    zp = Elliptic(np, rp, rs)
  elseif dm == "Chebyshev1"
    zp = Chebyshev1(np, rp)
  elseif dm == "Chebyshev2"
    zp = Chebyshev2(np, rs)
  else
    zp = Butterworth(np)
  end

  # polynomial ratio
  pr = convert(PolynomialRatio, digitalfilter(ff, zp))

  # zero-phase filter
  C.x[:] = filtfilt(pr, C.x)
  return nothing
end

function printcol(r::Float64)
  return r ≥ 1.00 ? 1 :
         r ≥ 0.75 ? 202 :
         r ≥ 0.50 ? 190 :
         r ≥ 0.25 ? 148 : 10
end

function safe_rm(file::String)
  try
    rm(file)
  catch err
    @warn(string("Can't remove ", file, ": throws error ", err))
  end
  return nothing
end

function compare_SeisHdr(H1::SeisHdr, H2::SeisHdr)
  for f in fieldnames(EQLoc)
    if typeof(getfield(H1.loc, f)) <: AbstractFloat
      @test isapprox(getfield(H1.loc, f), getfield(H2.loc,f))
    else
      @test getfield(H1.loc, f) == getfield(H2.loc,f)
    end
  end
  @test H1.mag.val ≈ H2.mag.val
  @test H1.mag.gap ≈ H2.mag.gap
  @test H1.mag.src == H2.mag.src
  @test H1.mag.scale == H2.mag.scale
  @test H1.mag.nst == H2.mag.nst
  @test H1.id == H2.id
  @test H1.ot == H2.ot
  @test H1.src == H2.src
  @test H1.typ == H2.typ
  return nothing
end

function compare_SeisSrc(R1::SeisSrc, R2::SeisSrc)
  @test R1.id == R2.id
  @test R1.eid == R2.eid
  @test R1.m0 ≈ R2.m0
  @test R1.mt ≈ R2.mt
  @test R1.dm ≈ R2.dm
  @test R1.gap ≈ R2.gap
  @test R1.pax ≈ R2.pax
  @test R1.planes ≈ R2.planes
  @test R1.src == R2.src
  @test R1.st.desc == R2.st.desc
  @test R1.st.dur ≈ R2.st.dur
  @test R1.st.rise ≈ R2.st.rise
  @test R1.st.decay ≈ R2.st.decay
  return nothing
end

function compare_SeisData(S1::SeisData, S2::SeisData)
  sort!(S1)
  sort!(S2)
  @test S1.id == S2.id
  @test S1.name == S2.name
  @test S1.units == S2.units
  @test isapprox(S1.fs, S2.fs)
  @test isapprox(S1.gain, S2.gain)
  for i in 1:S1.n
    L1 = S1.loc[i]
    L2 = S2.loc[i]
    @test isapprox(L1.lat, L2.lat)
    @test isapprox(L1.lon, L2.lon)
    @test isapprox(L1.el, L2.el)
    @test isapprox(L1.dep, L2.dep)
    @test isapprox(L1.az, L2.az)
    @test isapprox(L1.inc, L2.inc)

    R1 = S1.resp[i]
    R2 = S2.resp[i]
    for f in fieldnames(PZResp)
      @test isapprox(getfield(R1, f), getfield(R2,f))
    end

    @test S1.t[i] == S2.t[i]
    @test isapprox(S1.x[i],S2.x[i])
  end
  return nothing
end

function compare_events(Ev1::SeisEvent, Ev2::SeisEvent)
  compare_SeisHdr(Ev1.hdr, Ev2.hdr)
  compare_SeisSrc(Ev1.source, Ev2.source)
  S1 = convert(SeisData, Ev1.data)
  S2 = convert(SeisData, Ev2.data)
  compare_SeisData(S1, S2)
  return nothing
end

function rse_wb(n::Int64)
  Ev = randSeisEvent(n, s=1.0)

  Ev.source.misc = Dict{String,Any}(
  "pax_desc"    => "azimuth, plunge, length",
  "mt_id"       => "smi:SeisIO/moment_tensor;fmid="*Ev.source.id,
  "planes_desc" => "strike, dip, rake")
  Ev.source.eid = Ev.hdr.id
  Ev.source.npol = 0
  Ev.source.src = Ev.source.id * "," * Ev.source.src
  Ev.source.notes = String[]

  Ev.hdr.int = (0x00, "")
  Ev.hdr.src = "randSeisHdr:" * Ev.hdr.id
  Ev.hdr.loc.src = Ev.hdr.id * "," * Ev.hdr.loc.src
  Ev.hdr.loc.datum = ""
  Ev.hdr.loc.typ = ""
  Ev.hdr.loc.rms = 0.0
  flags = bitstring(Ev.hdr.loc.flags)
  if flags[1] == '1' || flags[2] == '1'
    flags = "11" * flags[3:8]
    Ev.hdr.loc.flags = parse(UInt8, flags, base=2)
  end

  Ev.hdr.mag.src = Ev.hdr.loc.src * ","
  Ev.hdr.notes = String[]
  Ev.hdr.misc = Dict{String,Any}()

  for j in 1:Ev.data.n
    Ev.data.misc[j] = Dict{String,Any}()
    Ev.data.notes[j] = String[]
    Ev.data.loc[j] = GeoLoc(
      lat = (rand(0.0:1.0:89.0) + rand())*-1.0^rand(1:2),
      lon = (rand(0.0:1.0:179.0) + rand())*-1.0^rand(1:2),
      el = rand()*1000.0,
      dep = rand()*1000.0,
      az = (rand()-0.5)*180.0,
      inc = rand()*90.0
    )
    Δ = round(Int64, 1.0e6/Ev.data.fs[j])
    nt = size(Ev.data.t[j],1)
    k = trues(nt)
    for n in 2:nt-1
      if Ev.data.t[j][n,2] ≤ Δ || (Ev.data.t[j][n+1,1]-Ev.data.t[j][n,1] < 2)
        k[n] = false
      end
    end
    Ev.data.t[j] = Ev.data.t[j][k,:]
  end
  return Ev
end

function latlon2xy(xlat::Float64, xlon::Float64)
  s = sign(xlon)
  c = 111194.6976
  y = c*xlat
  d = acosd(cosd(xlon*s)*cosd(xlat))
  x = sqrt(c^2*d^2-y^2)
  return [round(Int32, s*x), round(Int32, y)]
end

# ===========================================================================
# Redirect info, warnings, and errors to the logger
out = open("runtests.log", "a")
logger = SimpleLogger(out)
global_logger(logger)
@info("stdout redirect and logging")

# Set some keyword defaults
SeisIO.KW.comp = 0x00
has_restricted = safe_isdir(path * "/SampleFiles/Restricted/")
keep_log = false
keep_samples = true


text:
#####
##### AbstractShape
#####

abstract type AbstractShape end

is_valid(::AbstractShape) = true

#####
##### functions for identifying where the shape lies with respect to the bounds of the image
#####

get_i_min(image::AbstractMatrix) = firstindex(image, 1)
get_i_max(image::AbstractMatrix) = lastindex(image, 1)
get_i_extrema(x) = (get_i_min(x), get_i_max(x))

get_j_min(image::AbstractMatrix) = firstindex(image, 2)
get_j_max(image::AbstractMatrix) = lastindex(image, 2)
get_j_extrema(x) = (get_j_min(x), get_j_max(x))

function is_outbounds(image::AbstractMatrix, shape::AbstractShape)
    i_min_shape, i_max_shape = get_i_extrema(shape)
    i_min_image, i_max_image = get_i_extrema(image)

    j_min_shape, j_max_shape = get_j_extrema(shape)
    j_min_image, j_max_image = get_j_extrema(image)

    return i_max_shape < i_min_image || i_min_shape > i_max_image || j_max_shape < j_min_image || j_min_shape > j_max_image
end

function is_inbounds(image::AbstractMatrix, shape::AbstractShape)
    i_min_shape, i_max_shape = get_i_extrema(shape)
    i_min_image, i_max_image = get_i_extrema(image)

    j_min_shape, j_max_shape = get_j_extrema(shape)
    j_min_image, j_max_image = get_j_extrema(image)

    return i_min_shape >= i_min_image && i_max_shape <= i_max_image && j_min_shape >= j_min_image && j_max_shape <= j_max_image
end

#####
##### drawing a single pixel on the image
#####

function put_pixel!(image::AbstractMatrix, i, j, color)
    if checkbounds(Bool, image, i, j)
        put_pixel_inbounds!(image, i, j, color)
    end

    return nothing
end

function put_pixel_inbounds!(image::AbstractMatrix, i, j, color)
    @inbounds image[i, j] = color
    return nothing
end

#####
##### DrawingOptimizationStyle to define generic draw! methods with different levels of optimization for drawing shapes
#####

abstract type DrawingOptimizationStyle end

get_drawing_optimization_style(shape) = PUT_PIXEL

draw!(image::AbstractMatrix, shape::AbstractShape, color) = draw!(get_drawing_optimization_style(shape), image, shape, color)

#####
##### PutPixel (least optimized)
#####

struct PutPixel <: DrawingOptimizationStyle end
const PUT_PIXEL = PutPixel()

draw!(::PutPixel, image::AbstractMatrix, shape::AbstractShape, color) = draw!(put_pixel!, image, shape, color)

#####
##### CheckBounds (checks outbounds and inbounds conditions and dispatches accordingly)
#####

struct CheckBounds <: DrawingOptimizationStyle end
const CHECK_BOUNDS = CheckBounds()

function draw!(::CheckBounds, image::AbstractMatrix, shape::AbstractShape, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    if is_outbounds(image, shape)
        return nothing
    end

    if is_inbounds(image, shape)
        draw!(put_pixel_inbounds!, image, shape, color)
    else
        draw!(put_pixel!, image, shape, color)
    end

    return nothing
end

#####
##### Clip (clips the shapes that can be clipped easily and draw them with inbounds optimization)
#####

struct Clip <: DrawingOptimizationStyle end
const CLIP = Clip()

function draw!(::Clip, image::AbstractMatrix, shape::AbstractShape, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    if is_outbounds(image, shape)
        return nothing
    end

    draw!(put_pixel_inbounds!, image, clip(image, shape), color)

    return nothing
end

#####
##### PutPixelInbounds (directly draw assuming shape is inbounds)
#####

struct PutPixelInbounds <: DrawingOptimizationStyle end
const PUT_PIXEL_INBOUNDS = PutPixelInbounds()

draw!(::PutPixelInbounds, image::AbstractMatrix, shape::AbstractShape, color) = draw!(put_pixel_inbounds!, image, shape, color)


text:
@testset "Gaussian_WG_HNE_DL" begin
    θ₀ = [0.7, 0.9, 0.1, 0.3, 0.3]
    ω = 1.1
    @testset "Preallocation" begin
        model = Gaussian_WG_HNE_DL{1,40}(θ₀)
        @test model.norm == model.α/model.κ/sqrt(2π)
        @test model.halfκ⁻² == inv(2model.κ^2)
        @test model.cosϕₘ == cos(model.ϕₘ)
        @test model.sinϕₘ == sin(model.ϕₘ)
        @test model.cos2ϕₘ == cos(2model.ϕₘ)
        @test model.sin2ϕₘ == sin(2model.ϕₘ)
        @test model.invexp2σ² == exp(-2model.σ^2)
        @test model.invexphalfσ² == exp(-0.5model.σ^2)
    end
    @testset "Gradient" begin
        @test approx_gradient(θ -> sdf(Gaussian_WG_HNE_DL{1,40}(θ), ω), θ₀)      ≈ grad_sdf(Gaussian_WG_HNE_DL{1,40}(θ₀), ω)
        @test approx_gradient(θ -> sdf(Gaussian_WG_HNE_DL{1,40}(θ), 0), θ₀)      ≈ grad_sdf(Gaussian_WG_HNE_DL{1,40}(θ₀), 0)
    end
    @testset "Error handling" begin
        @test_throws ArgumentError Gaussian_WG_HNE_DL{1,40}(-0.7, 0.9, 0.1, 0.3, 0.3)
        @test_throws ArgumentError Gaussian_WG_HNE_DL{1,40}(0.7, -0.9, 0.1, 0.3, 0.3)
        @test_throws ArgumentError Gaussian_WG_HNE_DL{1,40}(0.7,  0.9, -0.1, 0.3, 0.3)
        @test_throws ArgumentError Gaussian_WG_HNE_DL{1,40}(0.7,  0.9, 0.1, 0.3, -0.3)
        @test_throws ArgumentError Gaussian_WG_HNE_DL{1,40}(ones(4))
        @test_throws ArgumentError Gaussian_WG_HNE_DL{1,40}(ones(6))
        @test_throws ErrorException Gaussian_WG_HNE_DL(ones(5))
        @test_throws ErrorException Gaussian_WG_HNE_DL(0.7, 0.9, 0.1, 0.3, 0.3)
    end
end


text:
module Haptics

const library = "libHPGE.so"

# lib = Libdl.dlopen("./libHPGE.so") # Open the library explicitly.
# sym = Libdl.dlsym(lib, :my_fcn)   # Get a symbol for the function to call.
# ccall(sym, ...) # Use the pointer `sym` instead of the (symbol, library) tuple (remaining arguments are the same).
# Libdl.dlclose(lib) # Close the library explicitly.

const requiredversion = v"2.3.6"

const string_buffer_length = Cint(100)

function get_version_info()
    major = Ref{Cint}(0); minor = Ref{Cint}(0); patch = Ref{Cint}(0);
    ccall((:get_version_info, "libHPGE.so"), Cvoid,
          (Ref{Cint}, Ref{Cint}, Ref{Cint}),
          major, minor, patch)
    VersionNumber(major[], minor[], patch[])
end

compatible() = get_version_info() == requiredversion

function get_error_msg(error::Integer)
    message = Vector{UInt8}(undef, string_buffer_length)
    res = ccall((:get_error_msg, "libHPGE.so"), Cint,
                (Cint, Cint, Cstring), error, string_buffer_length, pointer(message))
    if res != 0
        return "Could not get error message"
    end
    unsafe_string(pointer(message))
end

function last_error_msg()
    message = Vector{UInt8}(undef, string_buffer_length)
    res = ccall((:last_error_msg, "libHPGE.so"), Cint,
                (Cint, Cstring), string_buffer_length, pointer(message))
    if res != 0
        return "Could not get error message"
    end
    unsafe_string(pointer(message))
end

# int init_logging(
#     int device_coordinates,
#     int position,
#     int velocity, int force,
#     int interaction_forces, int object_number,
#     int[] objects);

#     public static extern int start_logging(int sampling_rate);
#     public static extern int is_logging();
#     public static extern int stop_logging_and_save(byte[] filename);

tick() = ccall((:tick, "libHPGE.so"), Cvoid, ())

count_devices()::Int = ccall((:count_devices, "libHPGE.so"), Cint, ())

function get_device_name(deviceid::Integer)
    devicename = Vector{UInt8}(undef, string_buffer_length)
    res = ccall((:get_device_name, "libHPGE.so"), Cint,
                (Cint, Cint, Cstring),
                deviceid, string_buffer_length, pointer(devicename))
    res != 0 && throw(get_error_msg(res))
    unsafe_string(pointer(devicename))
end

function initialize(deviceid::Integer; ws::Float64 = 10.0, tool::Float64 = 0.05)
    res = ccall((:initialize, "libHPGE.so"), Cint,
                (Cint, Cdouble, Cdouble), deviceid, ws, tool)
    # res != 0 && throw(get_error_msg(res))
    # true
end

function deinitialize()
    res = ccall((:deinitialize, "libHPGE.so"), Cint, ())
    # res != 0 && throw(get_error_msg(res))
    # true
end

function start()
    res = ccall((:start, "libHPGE.so"), Cint, ())
    # res != 0 && throw(get_error_msg(res))
    # true
end

function stop()
    res = ccall((:stop, "libHPGE.so"), Cint, ())
    # res != 0 && throw(get_error_msg(res))
    # true
end

get_loop_frequency() = ccall((:get_loop_frequency, "libHPGE.so"), Cdouble, ())
get_loops() = ccall((:get_loops, "libHPGE.so"), Cint, ())
get_log_frame() = ccall((:get_log_frame, "libHPGE.so"), Clong, ())
is_initialized() = ccall((:is_initialized, "libHPGE.so"), Cint, ()) == 0
is_running() = ccall((:is_running, "libHPGE.so"), Cint, ()) == 0
is_tool_button_pressed(button::Integer) =
    ccall((:is_tool_button_pressed, "libHPGE.so"), Cint, (Cint,), button) == 0

# log_annotate(byte[] note)
# set_hook, remove_hook

# public static extern int enable_dynamic_objects();
# public static extern int disable_dynamic_objects();
# public static extern void enable_wait_for_small_forces();
# public static extern void disable_wait_for_small_forces();
# public static extern void enable_rise_forces();
# public static extern void disable_rise_forces();

end

max_stars_repo_name: HapticPlugin/HPGE-wrappers

max_issues_repo_name: HapticPlugin/HPGE-wrappers

max_forks_repo_name: HapticPlugin/HPGE-wrappers


text:
struct SemForwardDiff <: SemDiff
    algorithm
    options #Optim.Options() call to optimize()
end

# function SemForwardDiff()
#     return SemForwardDiff(Optim.Options())
# end

struct SemReverseDiff{} <: SemDiff
    algorithm
    options
    # For preallocations, see the examples in ReverseDiff
end

struct SemFiniteDiff{} <: SemDiff
    algorithm
    options
end

# function SemFiniteDiff()
#     return SemFiniteDiff(Optim.Options())
# end

function SemAnalyticDiff(algorithm, options, functions)
    return SemAnalyticDiff(algorithm, options, functions, nothing)
end

function (diff::SemAnalyticDiff)(par, grad::AbstractVector, model)
    for i = 1:length(diff.functions)
        diff.functions[i](par, grad, model)
    end
end

function (diff::SemAnalyticDiff)(par, H::AbstractMatrix, model)
    for i = 1:length(diff.functions)
        diff.hessian_functions[i](par, H, model)
    end
end

function (model::Sem{A, B, C, D}
    where {A, B, C, D <: SemAnalyticDiff})(par, grad)
    grad .= zero(eltype(grad))
    model.imply(par, model)
    model.diff(par, grad, model)
end


text:
using GeoData, Test, Statistics, Dates
import ArchGDAL, NCDatasets, HDF5
using GeoData: Time, window, name

testpath = joinpath(dirname(pathof(GeoData)), "../test/")
include(joinpath(testpath, "test_utils.jl"))

# TODO example files without a login requirement
path1 = joinpath(testpath, "data/SMAP_L4_SM_gph_20160101T223000_Vv4011_001.h5")
path2 = joinpath(testpath, "data/SMAP_L4_SM_gph_20160102T223000_Vv4011_001.h5")

if isfile(path1) && isfile(path2)
    @testset "stack" begin
        smapstack = stack(path1)

        @testset "read" begin
            st = read(smapstack)
            @test st isa GeoStack
            @test st.data isa NamedTuple
            @test first(st.data) isa GeoArray
            @test parent(first(st.data)) isa Array
        end

        @testset "conversion to GeoArray" begin
            smaparray = smapstack["soil_temp_layer1"][Y(), X()]
            @test smaparray isa GeoArray{Float32,2}
            @test dims(smaparray) isa Tuple{<:X{<:Array{Float32,1}}, <:Y{<:Array{Float32,1}}}
            @test span(smaparray) isa Tuple{Irregular{Tuple{Float32,Float32}},Irregular{Tuple{Float32,Float32}}}
            @test span(smaparray) == (Irregular((-180.0f0, 180.0f0)), Irregular((-85.04456f0, 85.04456f0)))
            @test bounds(smaparray) == ((-180.0f0, 180.0f0), (-85.04456f0, 85.04456f0))
            @test index(smaparray) isa Tuple{Vector{Float32},Vector{Float32}}
            @test refdims(smaparray) isa Tuple{<:Ti}
            @test missingval(smaparray) == -9999.0
            @test smaparray[1] == -9999.0
            @test name(smaparray) == :soil_temp_layer1
            dt = DateTime(2016, 1, 1, 22, 30)
            step_ = Hour(3)
            @test refdims(smapstack) ==
                (Ti(dt:step_:dt; mode=Sampled(Ordered(), Regular(step_), Intervals(Start()))),)
            # Currently empty
            @test metadata(smaparray) isa Metadata{:SMAP}
        end

        @testset "conversion to GeoStack" begin
            # Stack Constructors
            # This uses too much ram! There is a lingering memory leak in HDF5.
            # geostack = GeoStack(stack)
            # @test Symbol.(Tuple(keys(stack))) == keys(geostack)
            geostack = GeoStack(smapstack; keys=(:baseflow_flux, :snow_mass, :soil_temp_layer1))
            @test keys(geostack) == (:baseflow_flux, :snow_mass, :soil_temp_layer1)
        end

        if VERSION > v"1.1-"
            @testset "copy" begin
                geoA = zero(smapstack[:soil_temp_layer1])
                @test geoA isa GeoArray
                @test geoA != smapstack[:soil_temp_layer1]
                copy!(geoA, smapstack, :soil_temp_layer1)
                @test geoA == smapstack[:soil_temp_layer1]
            end
        end

        @testset "window" begin
            windowedstack = SMAPstack(path1; window=(Y(1:5), X(1:5), Ti(1)))
            @test window(windowedstack) == (Y(1:5), X(1:5), Ti(1))
            windowedarray = windowedstack[:soil_temp_layer1]
            @test size(windowedarray) == (5, 5)
            @test windowedarray[1:3, 2:2] == reshape([-9999.0, -9999.0, -9999.0], 3, 1)
            @test windowedarray[1:3, 2] == [-9999.0, -9999.0, -9999.0]
            @test windowedarray[1, 2] == -9999.0
            windowedstack = SMAPstack(path1; window=(Y(1:5), X(1:5), Ti(1:1)))
            windowedarray = windowedstack[:soil_temp_layer1]
            @test windowedarray[1:3, 2:2, 1] == reshape([-9999.0, -9999.0, -9999.0], 3, 1)
            @test windowedarray[1:3, 2, 1] == [-9999.0, -9999.0, -9999.0]
            @test windowedarray[1, 2, 1] == -9999.0
            windowedstack = SMAPstack(path1; window=(Ti(1),))
            windowedarray = windowedstack[:soil_temp_layer1]
            @test windowedarray[1:3, 2:2] == reshape([-9999.0, -9999.0, -9999.0], 3, 1)
            @test windowedarray[1:3, 2] == [-9999.0, -9999.0, -9999.0]
            @test windowedarray[1, 2] == -9999.0
        end

        @testset "show" begin
            sh1 = sprint(show, MIME("text/plain"), smapstack[:soil_temp_layer1])
            # Test but don't lock this down too much
            @test occursin("GeoArray", sh1)
            @test occursin("Y", sh1)
            @test occursin("X", sh1)
            @test occursin("Ti", sh1)
            sh2 = sprint(show, MIME("text/plain"), smapstack[:soil_temp_layer1][Y(Between(0, 100)), X(Between(1, 100))])
            # Test but don't lock this down too much
            @test occursin("GeoArray", sh2)
            @test occursin("Y", sh2)
            @test occursin("X", sh2)
            @test occursin("Ti", sh2)
        end

    end

    @testset "series" begin
        smapseries = SMAPseries([path1, path2])
        val.(dims(smapseries))
        @test smapseries[1] isa SMAPstack
        @test first(bounds(smapseries, Ti)) == DateTime(2016, 1, 1, 22, 30)
        @test last(bounds(smapseries, Ti)) == DateTime(2016, 1, 3, 1, 30)
        modified_series = modify(Array, smapseries)
        stackkeys = keys(modified_series[1])
        @test typeof(modified_series) <: GeoSeries{<:GeoStack{<:NamedTuple{stackkeys,<:Tuple{<:GeoArray{Float32,2,<:Tuple,<:Tuple,<:Array{Float32,2}},Vararg}}}}

        @testset "read" begin
            geoseries = read(smapseries)
            @test geoseries isa GeoSeries{<:GeoStack}
            @test geoseries.data isa Vector{<:GeoStack}
            @test first(geoseries.data[1].data) isa GeoArray
        end

        @testset "show" begin
            sh = sprint(show, MIME("text/plain"), smapseries)
            # Test but don't lock this down too much
            @test occursin("GeoSeries", sh)
            # @test occursin("SMAPstack", sh)
            # @test occursin("Ti", sh)
        end
    end
end


text:
module Walk1DMDP

using POMDPs, POMDPToolbox, Distributions, Plots, MCTS

export
    Walk1DParams,
    Walk1DState,
    Walk1D,
    RandomGaussian

struct Walk1DParams
    startx::Float64
    threshx::Float64
    t_max::Int
end
Walk1DParams() = Walk1DParams(1.0, 10.0, 20)

struct Walk1DState
    t::Int
    x::Float64
end
mutable struct Walk1D <: MDP{Walk1DState,Float64}
    p::Walk1DParams
    distr::Distribution
end
Walk1D() = Walk1D(Walk1DParams(), Normal(0.0, 1.0))

immutable RandomGaussian <: Policy
    distr::Distribution
    rng::AbstractRNG
end
RandomGaussian(rng::AbstractRNG=MersenneTwister(0), m::Float64=0.0, s::Float64=1.0) = RandomGaussian(Normal(m,s), rng)
function POMDPs.action(pol::RandomGaussian, s::Walk1DState)
    return rand(pol.rng, pol.distr)
end
function POMDPs.initial_state(mdp::Walk1D, rng::AbstractRNG)
    Walk1DState(0, mdp.p.startx)
end
function POMDPs.generate_s(mdp::Walk1D, s::Walk1DState, a::Float64, rng::AbstractRNG)
    Walk1DState(s.t+1, s.x+a)
end
POMDPs.discount(mdp::Walk1D) = 1.0
isevent(mdp::Walk1D, s::Walk1DState) = s.x > mdp.p.threshx
POMDPs.isterminal(mdp::Walk1D, s::Walk1DState) = isevent(mdp, s) || s.t >= mdp.p.t_max
function POMDPs.reward(mdp::Walk1D, s::Walk1DState, a::Float64, sp::Walk1DState)
    r = log(pdf(mdp.distr, a))
    if isterminal(mdp, sp)
        if isevent(mdp, sp)
            r += 0.0
        else
            r += -miss_distance(mdp, sp)
        end
    end
    r
end
function miss_distance(mdp::Walk1D, s::Walk1DState)
    max(mdp.p.threshx-abs(s.x), 0.0)
end
MCTS.next_action(pol::RandomGaussian, sim::Walk1D, s::Walk1DState, snode) = action(pol, s)
Base.string(s::Walk1DState) = "($(s.t),$(round(s.x,2)))"

include("visualization.jl")

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from Euler angle and axis to quaternion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angleaxis_to_quat.jl
# ============================================

# Functions: angleaxis_to_quat
# ----------------------------

@testset "Euler angle and axis => Quaternion" begin
    for T in (Float32, Float64)
        # Sample a random Euler angle and axis.
        v = @SVector randn(T, 3)
        v = v / norm(v)
        a = _rand_ang(T)
        av = EulerAngleAxis(a, v)
        q = angleaxis_to_quat(av)
        s = (cos(a / 2) < 0) ? -1 : 1
        @test eltype(q) === T
        @test q.q0 ≈ s * cos(a / 2)
        @test q.q1 ≈ v[1] * sin(a / 2)
        @test q.q2 ≈ v[2] * sin(a / 2)
        @test q.q3 ≈ v[3] * sin(a / 2)
    end
end


text:
### FlatS02 types
# (FlatS02's are just FieldTuple's of a spin-0 and a spin-2)
const FlatIQUMap{P,T,M}     = FieldTuple{IQUMap,     NamedTuple{(:I,:P),Tuple{FlatMap{P,T,M},    FlatQUMap{P,T,M}}},     T}
const FlatIQUFourier{P,T,M} = FieldTuple{IQUFourier, NamedTuple{(:I,:P),Tuple{FlatFourier{P,T,M},FlatQUFourier{P,T,M}}}, Complex{T}}
const FlatIEBMap{P,T,M}     = FieldTuple{IEBMap,     NamedTuple{(:I,:P),Tuple{FlatMap{P,T,M},    FlatEBMap{P,T,M}}},     T}
const FlatIEBFourier{P,T,M} = FieldTuple{IEBFourier, NamedTuple{(:I,:P),Tuple{FlatFourier{P,T,M},FlatEBFourier{P,T,M}}}, Complex{T}}
# some handy Unions
const FlatS02{P,T,M} = Union{FlatIQUMap{P,T,M},FlatIQUFourier{P,T,M},FlatIEBMap{P,T,M},FlatIEBFourier{P,T,M}}
const FlatIQU{P,T,M} = Union{FlatIQUMap{P,T,M},FlatIQUFourier{P,T,M}}
const FlatIEB{P,T,M} = Union{FlatIEBMap{P,T,M},FlatIEBFourier{P,T,M}}
const FlatS02Map{P,T,M} = Union{FlatIQUMap{P,T,M},FlatIEBMap{P,T,M}}
const FlatS02Fourier{P,T,M} = Union{FlatIQUFourier{P,T,M},FlatIEBFourier{P,T,M}}

spin(::Type{<:FlatS02}) = S02

### convenience constructors
for (F,F2,F0,(X,Y,Z),(X′,Y′,Z′),T) in [
        (:FlatIQUMap,     :FlatQUMap,     :FlatMap,     (:Ix,:Qx,:Ux), (:I,:Q,:U), :T),
        (:FlatIQUFourier, :FlatQUFourier, :FlatFourier, (:Il,:Ql,:Ul), (:I,:Q,:U), :(Complex{T})),
        (:FlatIEBMap,     :FlatEBMap,     :FlatMap,     (:Ix,:Ex,:Bx), (:I,:E,:B), :T),
        (:FlatIEBFourier, :FlatEBFourier, :FlatFourier, (:Il,:El,:Bl), (:I,:E,:B), :(Complex{T}))
    ]

    doc = """
        # main constructors:
        $F(
            $X::AbstractMatrix, $Y::AbstractMatrix, $Z::AbstractMatrix; $(F0==:FlatFourier ? "\n        Nside, # required, size of the map in pixels" : "")
            θpix,  # optional, resolution in arcmin (default: 1)
            ∂mode, # optional, fourier∂ or map∂ (default: fourier∂)
        )
        $F($X′::$F0, $Y′::$F0, $Z′::$F0)

        # more low-level:
        $F{P}($X::AbstractMatrix, $Y::AbstractMatrix, $Z::AbstractMatrix) # specify pixelization P explicilty
        $F{P,T}($X::AbstractMatrix, $Y::AbstractMatrix, $Z::AbstractMatrix) # additionally, convert elements to type $T
        $F{P,T,M<:AbstractMatrix{$T}}($X::M, $Y::M, $Z::M) # specify everything explicilty

    Construct a `$F` object. The top form of the constructors is most convenient
    for interactive work, while the others may be more useful for low-level code.
    """
    @eval begin
        @doc $doc $F
        $F($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array, $Z::AbstractRank2or3Array; kwargs...) =
            $F($F0($X; kwargs...), $F2($Y,$Z; kwargs...))
        $F{P}($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array, $Z::AbstractRank2or3Array) where {P} =
            $F{P}($F0{P}($X), $F2{P}($Y,$Z))
        $F{P,T}($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array, $Z::AbstractRank2or3Array) where {P,T} =
            $F{P,T}($F0{P,T}($X), $F2{P,T}($Y,$Z))
        $F{P,T,M}($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array, $Z::AbstractRank2or3Array) where {P,T,M} =
            $F{P,T,M}($F0{P,T,M}($X), $F2{P,T,M}($Y,$Z))
        $F($X′::$F0, $Y′::$F0, $Z′::$F0) = $F($X′, $F2($Y′,$Z′))
    end
end

### properties
propertynames(f::FlatS02) = (sort([:I, :P, :IP, propertynames(f.I)..., propertynames(f.P)...])...,)
getproperty(f::FlatIQU, s::Val{:IP}) = f
getproperty(f::FlatIQU, s::Union{Val{:Q},Val{:U}}) = getproperty(getfield(f,:fs).P,s)
getproperty(f::FlatIEB, s::Union{Val{:E},Val{:B}}) = getproperty(getfield(f,:fs).P,s)
getproperty(f::FlatS02, s::Union{Val{:Qx},Val{:Ux},Val{:Ex},Val{:Bx},Val{:Ql},Val{:Ul},Val{:El},Val{:Bl}}) = getproperty(getfield(f,:fs).P,s)
getproperty(f::FlatS02, s::Union{Val{:Ix},Val{:Il}}) = getproperty(getfield(f,:fs).I,s)
function getindex(f::FlatS02, k::Symbol; kwargs...)
    @match k begin
        (:IP) => f
        (:I || :P) => getfield(f.fs,k)
        (:Q || :U || :E || :B) => getindex(f.P,k; kwargs...)
        (:Ix || :Il) => getindex(f.I,k; kwargs...)
        (:Qx || :Ux || :Ql || :Ul || :Ex || :Bx || :El || :Bl) => getindex(f.P,k; kwargs...)
        _ => throw(ArgumentError("Invalid FlatS02 index: $k"))
    end
end

# A Flat TEB covariance of the form:
#
#    [ΣTT ΣTE  ⋅
#     ΣTE ΣEE  ⋅
#      ⋅   ⋅  ΣBB]
#
# We store the 2x2 block as a 2x2 SMatrix, ΣTE, so that we can easily call sqrt/inv on
# it, and the ΣBB block separately as ΣB.
struct FlatIEBCov{T,F} <: ImplicitOp{IEBFourier,S02,Pix}
    ΣTE :: SMatrix{2,2,Diagonal{T,F},4}
    ΣB :: Diagonal{T,F}
end

# contructing from Cℓs
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S02}, CℓTT, CℓEE, CℓBB, CℓTE; kwargs...) where {T,P}
    ΣTT, ΣEE, ΣBB, ΣTE = [Cℓ_to_Cov(P,T,S0,Cℓ; kwargs...) for Cℓ in (CℓTT,CℓEE,CℓBB,CℓTE)]
    FlatIEBCov(@SMatrix([ΣTT ΣTE; ΣTE ΣEE]), ΣBB)
end

# applying the operator
*(L::FlatIEBCov, f::FlatS02) =       L * IEBFourier(f)
\(L::FlatIEBCov, f::FlatS02) = pinv(L) * IEBFourier(f)
function *(L::FlatIEBCov, f::FlatIEBFourier)
    (i,e),b = (L.ΣTE * [f.I, f.E]), L.ΣB * f.B
    FlatIEBFourier(i,e,b)
end

# manipulating the operator
adjoint(L::FlatIEBCov) = L
sqrt(L::FlatIEBCov) = FlatIEBCov(sqrt(L.ΣTE), sqrt(L.ΣB))
pinv(L::FlatIEBCov) = FlatIEBCov(pinv(L.ΣTE), pinv(L.ΣB))
simulate(rng::AbstractRNG, L::FlatIEBCov{Complex{T},FlatFourier{P,T,M}}) where {P,T,M} = sqrt(L) * white_noise(rng, FlatIEBFourier{P,T,M})
global_rng_for(::Type{FlatIEBCov{T,F}}) where {T,F} = global_rng_for(F)
diag(L::FlatIEBCov) = FlatIEBFourier(L.ΣTE[1,1].diag, L.ΣTE[2,2].diag, L.ΣB.diag)
similar(L::FlatIEBCov) = FlatIEBCov(similar.(L.ΣTE), similar(L.ΣB))

# getindex
function getindex(L::FlatIEBCov, k::Symbol)
    @match k begin
        :IP => L
        :I => L.ΣTE[1,1]
        :E => L.ΣTE[2,2]
        :B => L.ΣB
        :P => Diagonal(FlatEBFourier(L[:E].diag, L[:B].diag))
        (:QQ || :UU || :QU || :UQ) => getindex(L[:P], k)
        _ => throw(ArgumentError("Invalid FlatIEBCov index: $k"))
    end
end

# FlatIEBCov arithmetic
*(L::FlatIEBCov, D::DiagOp{<:FlatIEBFourier}) = FlatIEBCov(SMatrix{2,2}(L.ΣTE * [[D[:I]] [0]; [0] [D[:E]]]), L.ΣB * D[:B])
+(L::FlatIEBCov, D::DiagOp{<:FlatIEBFourier}) = FlatIEBCov(@SMatrix[L.ΣTE[1,1]+D[:I] L.ΣTE[1,2]; L.ΣTE[2,1] L.ΣTE[2,2]+D[:E]], L.ΣB + D[:B])
*(La::F, Lb::F) where {F<:FlatIEBCov} = F(La.ΣTE * Lb.ΣTE, La.ΣB * Lb.ΣB)
+(La::F, Lb::F) where {F<:FlatIEBCov} = F(La.ΣTE + Lb.ΣTE, La.ΣB + Lb.ΣB)
+(L::FlatIEBCov, U::UniformScaling{<:Scalar}) = FlatIEBCov(@SMatrix[(L.ΣTE[1,1]+U) L.ΣTE[1,2]; L.ΣTE[2,1] (L.ΣTE[2,2]+U)], L.ΣB+U)
*(L::FlatIEBCov, λ::Scalar) = FlatIEBCov(L.ΣTE * λ, L.ΣB * λ)
*(D::DiagOp{<:FlatIEBFourier}, L::FlatIEBCov) = L * D
+(U::UniformScaling{<:Scalar}, L::FlatIEBCov) = L + U
*(λ::Scalar, L::FlatIEBCov) = L * λ
copyto!(dst::Σ, src::Σ) where {Σ<:FlatIEBCov} = (copyto!(dst.ΣB, src.ΣB); copyto!.(dst.ΣTE, src.ΣTE); dst)

function get_Cℓ(f1::FlatS02, f2::FlatS02=f1; which=(:II,:EE,:BB,:IE,:IB,:EB), kwargs...)
    Cℓ = (;[Symbol(x1*x2) => get_Cℓ(getindex(f1,Symbol(x1)),getindex(f2,Symbol(x2)); kwargs...) for (x1,x2) in split.(string.(ensure1d(which)),"")]...)
    which isa Symbol ? Cℓ[1] : Cℓ
end


text:
# The Hausdorff measure generalizes the notion of volume or area to a manifold.
# For manifolds embedded in real space, it can be intuitively understoof as the addition
# of an infinitesimal dimension and the limit of the volume in the embedded space as the
# infinitesimal dimension goes to 0.
# A more useful description comes from the area formula. Given an m-
# dimensional manifold M and a Lipschitz map f: ℝᵐ → M ⊆ ℝⁿ that embeds M in ℝⁿ, then the Hausdorff
# measure dμ(x) on M is is the pushforward of the Lebesgue measure through f, given as
# dμ(x) = det(G(x))^(-1/2) dλ(f⁻¹(x)), where G(x) is the matrix representation of the metric tensor
# in M, given by G(x) = J(x)ᵀ J(x), where J is the Jacobian of f.
# This is the un-normalized Hausdorff measure, which integrates to the area/volume of M.

"""
    Hausdorff{M,A} <: PrimitiveMeasure

The un-normalized Hausdorff measure on a manifold.

The Hausdorff measure generalizes the notion of area or volume to a manifold that is embedded
in a metric space. That is, the mass of the measure over some region of the manifold is the
area/volume of that region in the embedded space.

# Constructors

    Hausdorff(M::AbstractManifold)

Constructs the Hausdorff measure for the manifold `M` using the default embedding of the
manifold.
"""
struct Hausdorff{M,A} <: PrimitiveMeasure
    manifold::M
    atlas::A  # optional, when the Hausdorff measure is constructed from a Lebesgue measure
end
Hausdorff(M) = Hausdorff(M, nothing)

function Base.show(io::IO, mime::MIME"text/plain", μ::Hausdorff)
    print(io, "Hausdorff")
    return show(io, mime, (μ.manifold, μ.atlas))
end

Manifolds.base_manifold(μ::Hausdorff) = μ.manifold

MeasureTheory.logdensity(::Hausdorff, x) = zero(eltype(x))

function Base.rand(rng::AbstractRNG, T::Type, μ::Normalized{<:Hausdorff})
    p = default_point(μ, T)
    return Random.rand!(rng, p, μ)
end

# Stiefel

# In general, given matrix z ∈ 𝔽^(n × k) with IID std normal elements, p=z(z'z)^(-1/2) is drawn
# from the Hausdorff measure on St(n,k,𝔽).
# This implementation uses the unique QR decomposition z=QR where R[i,i] is real and positive
# for all i, so that Q ~ Hausdorff(St(n,k,𝔽)).
# See Theorem 2.3.19 of Gupta AK, Nagar DK. Matrix variate distributions. CRC Press; 2018
function Random.rand!(
    rng::AbstractRNG, p::AbstractMatrix, ::Normalized{<:Hausdorff{<:Stiefel}}
)
    randn!(rng, p)
    Q, _ = qr_unique!(p)
    return Q
end

# Chikuse, 2003 Eq. 1.4.8
function logmass(::Hausdorff{Stiefel{n,k,ℝ}}) where {n,k}
    halfn = n//2
    return k * logtwo + (k * halfn) * logπ - logmvgamma(k, halfn)
end
# vol(St(n,k,𝔽)) = vol(𝔽𝕊ⁿ⁻¹) * vol(St(n-1,k-1,𝔽)), Chikuse, 2003 §1.4.4
# then, vol(St(n,k,𝔽)) = 2ᵏ π^(k(2n-k+1)d/4) / ∏ᵢ₌₁ᵏΓ(d(n-k+i)/2)),
# where d = real_dimension(𝔽)
function logmass(::Hausdorff{Stiefel{n,k,ℂ}}) where {n,k}
    a = n - k
    lfact = logfactorial(a)
    lm = k * logtwo + (k * (2n - k + 1)//2) * logπ - lfact
    for _ in 1:(k - 1)
        a += 1
        lfact += log(a)
        lm -= lfact
    end
    return lm
end
function logmass(::Hausdorff{Stiefel{n,k,ℍ}}) where {n,k}
    a = 2(n - k) + 1
    lfact = logfactorial(a)
    lm = k * logtwo + (k * (a + k)) * logπ - lfact
    for _ in 1:(k - 1)
        a += 2
        lfact += log(a) + log(a - 1)
        lm -= lfact
    end
    return lm
end

# Grassmann

function Random.rand!(
    rng::AbstractRNG, p::AbstractMatrix, ::Normalized{<:Hausdorff{Grassmann{n,k,𝔽}}}
) where {n,k,𝔽}
    return rand!(rng, p, normalize(Hausdorff(Stiefel(n, k, 𝔽))))
end

# because Gr(n,k,𝔽) = St(n,k,𝔽) / U(k,𝔽),
# then vol(Gr(n,k,𝔽)) = vol(St(n,k,𝔽)) / vol(St(k,k,𝔽))
#                     = π^(k((n-k)d/2)) ∏ᵢ₌₁ᵏ Γ(id/2) / Γ((n-k)d/2 + id/2),
# where d = real_dimension(𝔽)
function logmass(::Hausdorff{Grassmann{n,k,ℝ}}) where {n,k}
    return k * (n - k)//2 * logπ + logmvgamma(k, k//2) - logmvgamma(k, n//2)
end
function logmass(::Hausdorff{Grassmann{n,k,ℂ}}) where {n,k}
    a = n - k
    lfact1 = logfactorial(a)
    lfact2 = log(1)
    lm = k * a * logπ + lfact2 - lfact1
    for i in 1:(k - 1)
        a += 1
        lfact1 += log(a)
        lfact2 += log(i)
        lm += lfact2 - lfact1
    end
    return lm
end
function logmass(::Hausdorff{Grassmann{n,k,ℍ}}) where {n,k}
    a = 2(n - k) + 1
    b = 1
    lfact1 = logfactorial(a)
    lfact2 = log(b)
    lm = 2 * k * (n - k) * logπ + lfact2 - lfact1
    for _ in 1:(k - 1)
        a += 2
        b += 2
        lfact1 += log(a) + log(a - 1)
        lfact2 += log(b) + log(b - 1)
        lm += lfact2 - lfact1
    end
    return lm
end

# Circle

default_point(::Normalized{<:Hausdorff{Circle{ℝ}}}, T::Type) = zero(float(real(T)))
default_point(::Normalized{<:Hausdorff{Circle{ℂ}}}, T::Type) = zero(complex(float(real(T))))

function Random.rand!(rng::AbstractRNG, p::Real, ::Normalized{<:Hausdorff{Circle{ℝ}}})
    return rand(rng, typeof(p)) * twoπ - π
end
function Random.rand!(rng::AbstractRNG, p::Complex, ::Normalized{<:Hausdorff{Circle{ℂ}}})
    return sign(randn(rng, typeof(p)))
end

logmass(::Hausdorff{<:Circle}) = log2π

# Sphere

function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, μ::Normalized{<:Hausdorff{<:AbstractSphere}}
)
    return normalize!(randn!(rng, p))
end

function logmass(μ::Hausdorff{<:AbstractSphere{𝔽}}) where {𝔽}
    m = prod(representation_size(base_manifold(μ)))
    ν = real_dimension(𝔽) * m//2
    return logtwo + ν * logπ - loggamma(ν)
end

# ProjectiveSpace

function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, μ::Normalized{<:Hausdorff{<:AbstractProjectiveSpace}}
)
    return normalize!(randn!(rng, p))
end

# because 𝔽ℙⁿ = 𝔽𝕊ⁿ / 𝔽𝕊⁰, then vol(𝔽ℙⁿ) = vol(𝔽𝕊ⁿ) / vol(𝔽𝕊⁰)
function logmass(μ::Hausdorff{<:AbstractProjectiveSpace{ℝ}})
    m = prod(representation_size(base_manifold(μ)))
    ν = m//2
    return ν * logπ - loggamma(ν)
end
function logmass(μ::Hausdorff{<:AbstractProjectiveSpace{ℂ}})
    m = prod(representation_size(base_manifold(μ)))
    n = m - 1
    return n * logπ - logfactorial(n)
end
function logmass(μ::Hausdorff{<:AbstractProjectiveSpace{ℍ}})
    m = prod(representation_size(base_manifold(μ)))
    n = m - 1
    return 2n * logπ - loggamma(2(n + 1))
end

# Rotations/SpecialOrthogonal

function Random.rand!(
    rng::AbstractRNG,
    p::AbstractMatrix,
    ::Normalized{<:Hausdorff{<:Union{Rotations,SpecialOrthogonal}}},
)
    # draw p ~ O(n)
    rand!(rng, p, Hausdorff(Stiefel(n, n)))
    # fix to p ~ SO(n)
    if det(p) < 0
        @views p[:, [1, 2]] = p[:, [2, 1]]
    end
    return p
end

# since O(n) = St(n,n) is essentially two copies of SO(n), then vol(SO(n)) = vol(St(n,n)) / 2
function logmass(::Hausdorff{<:Union{Rotations{n},SpecialOrthogonal{n}}}) where {n}
    return logmass(Hausdorff(Stiefel(n, n))) - logtwo
end


text:
# Abstract type: whether the gating function has random effects
# abstract type RandomGating end
# struct HasRE <: RandomGating end
# struct NonRE <: RandomGating end

# # Abstract type: AnyGating
# abstract type AnyGating{s<:RandomGating} end
# # Logit Gating
# const NonRandomGating = AnyGating{NonRE}

const gating = [
    "logit"
]

for dname in gating
    include(joinpath("gating", "$(dname).jl"))
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
function exclude(a::String,b::String)
    a2 = deepcopy(a)
    for idx in 1:length(b)
        i = b[idx:idx]
        if occursin(i,a2)
            a2 = replace(a2,i=>"")
        end
    end
    a2
end
function getdims(iC,iA,iB,A,B)
    sza = size(A)
    szb = size(B)
    szc = []
    for i in 1:length(iA)
        if occursin(iA[i:i],iC)
            push!(szc,sza[i])
        end
    end
    for i in 1:length(iB)
        if occursin(iB[i:i],iC)
            push!(szc,szb[i])
        end
    end
    szc
end
function contract(A,B,iA,iB)
    # figure out size of C, iC
    contract(A,B,1.0,1.0,iA,iB)
end
function contract(A,B,sA,sB,iA,iB)
    iC = exclude(iA,iB)*exclude(iB,iA)
    C = zeros(eltype(A),getdims(iC,iA,iB,A,B)...)
    contract!(C,A,B,1.0,sA,sB,iC,iA,iB)
    C
end

function contract!(C,A,B,iC,iA,iB)
    contract!(C,A,B,1.0,1.0,1.0,iC,iA,iB)
end

function contract!(C,A,B,sC::T,sA::T,sB::T,iC,iA,iB) where T <: AbstractFloat
    contract!(C,A,B,sC,sA,sB,iC,iA,iB,
              Fermi.ComputeEnvironment.interconnect,Fermi.ComputeEnvironment.communicator,
              Fermi.ComputeEnvironment.accelerator,Fermi.ComputeEnvironment.contractor)
end

function contract!(C,A,B, #tensors
                   sC,sA,sB, #scaling factors
                   iC,iA,iB, #indices
                   ic::I where I <: Fermi.Environments.AbstractInterconnect,
                   comm::Fermi.Environments.NoCommunicator,
                   acc::Fermi.Environments.NoAccelerator,
                   contractor::Fermi.Environments.TTGT)
    # transpose
    # transpose
    # GEMM
    # transpose
end

function contract!(C::T1,
                   A::T2,
                   B::T3,  #tensors
                   sC,sA,sB, #scaling factors
                   iC,iA,iB, #indices
                   ic::I where I <: Fermi.Environments.AbstractInterconnect,
                   comm::Fermi.Environments.NoCommunicator,
                   acc::Fermi.Environments.NoAccelerator,
                   contractor::Fermi.Environments.TBLIS) where { T1 <: Union{Array,Fermi.MemTensor},
                                                                 T2 <: Union{Array,Fermi.MemTensor},
                                                                 T3 <: Union{Array,Fermi.MemTensor}}

    Ap = Fermi.data(A)
    Bp = Fermi.data(B)
    Cp = Fermi.data(C)
    _A = TBLIS.TTensor{eltype(Ap)}(Ap,sA)
    _B = TBLIS.TTensor{eltype(Bp)}(Bp,sB)
    _C = TBLIS.TTensor{eltype(Cp)}(Cp,sC)
    TBLIS.mul!(_C,_A,_B,iA,iB,iC)
end

function contract!(C,A,B, #tensors
                   sC,sA,sB, #scaling factors
                   iC,iA,iB, #indices
                   ic::Fermi.Environments.MPI_IC,
                   comm::Fermi.Environments.MPI_World,
                   acc::Fermi.Environments.AMD_GPU,
                   contractor::Fermi.Environments.TTGT)
    # Super special MPI aware GPU code specialized for AMD products
end


text:
struct E_Gd_H <: System end
Overseer.requested_components(::E_Gd_H) = (Spin, Gd, Mn)
@inline @inbounds function Overseer.update(::E_Gd_H, m::AbstractLedger)
    field = m[H]
    if isempty(field)
        return
    end
    h = field[1]
    spin  = m[Spin]
    p     = singleton(m, Model)
    etot  = singleton(m, Etot)
    if isempty(m[Mn])
        for i = 1:8
            etot.e -= p.χ_Gd * h.v ⋅ spin[i].v
        end
    else
        for e in @entities_in(spin && m[Gd])
            etot.e -= p.χ_Gd * h.v ⋅ spin[i].v
        end
        for e in @entities_in(spin && m[Mn])
            etot.e -= p.χ_L * h.v ⋅ spin[e].v
        end
    end
end

struct E_L_H <: System end

@inline function Overseer.update(::E_L_H, m::AbstractLedger)
    field = m[H]
    if isempty(field)
        return
    end
    etot  = singleton(m, Etot)
    h = field[1]
    p     = singleton(m, Model)
    Lc = m[L]

    for l in Lc
        etot.e += 0.5 * p.χ_L * (l.v ⋅ h.v)^2
    end
end

struct E_LL <: System end

@inline function Overseer.update(::E_LL, m::AbstractLedger)
    etot  = singleton(m, Etot)
    p     = singleton(m, Model)
    Lc = m[L]

    @inbounds for l1 in @entities_in(Lc)
        for l2 in @entities_in(Lc)
            l1 == l2 && continue
            etot.e += 0.5 * p.γ * (Lc[l1].v ⋅ Lc[l2].v)^2
        end
    end
end

struct E_Gd_easy <: System end

@inline function Overseer.update(::E_Gd_easy, m::AbstractLedger)
    etot  = singleton(m, Etot)
    p     = singleton(m, Model)
    spin  = m[Spin]
    easy  = m[EasyAxis]

    @inbounds for e in @entities_in(spin && easy && m[Gd])
        etot.e -= p.K_Gd * (spin[e].v ⋅ easy[e].v)^2
    end
    @inbounds for e in @entities_in(spin && easy && m[Mn])
        etot.e -= p.K_L * (spin[e].v ⋅ easy[e].v)^2
    end
end

struct E_L_easy <: System end

function Overseer.update(::E_L_easy, m::AbstractLedger)
    etot = singleton(m, Etot)
    p    = singleton(m, Model)
    Lc   = m[L]
    easy = m[EasyAxis]

    @inbounds for e in @entities_in(Lc && easy)
        etot.e -= 0.5 * p.K_L * (Lc[e].v ⋅ easy[e].v)^2
    end

end

struct E_Gd_L <: System end

@muladd function Overseer.update(::E_Gd_L, m::AbstractLedger)
    etot     = singleton(m, Etot)
    p        = singleton(m, Model)
    Gd_spins = m[Spin]
    Lc       = m[L]
    L1 = Lc.data[1].v
    L2 = Lc.data[2].v
    J1 = p.J1
    J2 = p.J2
    @inbounds for i=1:8:length(Gd_spins)

        t2 = J1 * L2
        t3 = J2 * L1
        t1 = J1 * L1
        t4 = J2 * L2

        t5 = t2 + t3
        t6 = t1 + t4
        t7 = t1 - t4
        t8 = t3 - t2

        etot.e += (Gd_spins[i].v - Gd_spins[i+4].v) ⋅ t5 +
                  (Gd_spins[i+1].v - Gd_spins[i+5].v) ⋅ t6 +
                  (Gd_spins[i+2].v - Gd_spins[i+6].v) ⋅ t7 +
                  (Gd_spins[i+3].v - Gd_spins[i+7].v) ⋅ t8
    end
end

struct E_Gd_easy4th <: System end

function Overseer.update(::E_Gd_easy4th, m::AbstractLedger)
    etot     = singleton(m, Etot)
    p        = singleton(m, Model)
    Gd_spins = m[Spin]
    easy  = m[EasyAxis]

    @inbounds etot.e -= p.K_Gd_4th * ((Gd_spins[1].v ⋅ easy[1].v)^4 +
                          (Gd_spins[2].v ⋅ easy[2].v)^4 +
                          (Gd_spins[3].v ⋅ easy[3].v)^4 +
                          (Gd_spins[4].v ⋅ easy[4].v)^4 +
                          (Gd_spins[5].v ⋅ easy[5].v)^4 +
                          (Gd_spins[6].v ⋅ easy[6].v)^4 +
                          (Gd_spins[7].v ⋅ easy[7].v)^4 +
                          (Gd_spins[8].v ⋅ easy[8].v)^4)
end

Overseer.requested_components(::E_Gd_L) = (Spin, L, H, Etot, EasyAxis, Model,)

struct E_Dipole <: System end
Overseer.requested_components(::E_Dipole) = (Bond,)

# function Overseer.update(::E_Dipole, m::AbstractLedger)
#     etot = singleton(m, Etot)
#     p = singleton(m, Model)
#     spin = m[Spin]
#     bonds = m[Bond]
#     @inbounds for b in bonds
#         pr = b.r
#         s1 = spin[b.e1].v
#         s2 = spin[b.e2].v
#         etot.e -= 0.5 * p.dipole_strength * (3 * (s1 ⋅ pr) * (s2 ⋅ pr) - b.nr2 * s1⋅s2)*b.inv_nr5 #0.5 for double counting
#     end
# end
@muladd function Overseer.update(::E_Dipole, m::AbstractLedger)
    etot = singleton(m, Etot)
    p = singleton(m, Model)
    spin = m[Spin]
    bonds = m[Bond]
    @inbounds for b in bonds
        t4 = 0.5 * p.dipole_strength * b.inv_nr5 #0.5 for double counting
        pr = b.r
        s1 = spin[b.e1].v
        s2 = spin[b.e2].v
        t1 = 3 * s1 ⋅ pr
        t2 = s2 ⋅ pr
        t3 = b.nr2 * s1 ⋅ s2
        etot.e = etot.e + t4 * (t3 - t1 * t2)
    end
end

struct E_L2H2 <: System end

function Overseer.update(::E_L2H2, m::AbstractLedger)
    etot = singleton(m, Etot)
    p = singleton(m, Model)
    h = singleton(m, H)
    Lc = m[L]
    @inbounds for l in Lc
        etot.e -= p.l2h2_strength * (l.v⋅l.v * h.v ⋅ h.v)
    end
end

struct E_Pb <: System end

function Overseer.update(::E_Pb, m::AbstractLedger)
    etot = singleton(m, Etot)
    p = singleton(m, Model)
    Ls = m[L]
    L1 = Ls[1].v
    L2 = Ls[2].v
    Gd_spins = m[Spin]
    pb = singleton(m, Pb).pb
    for i=1:8:length(Gd_spins)
        etot.e -= ((Gd_spins[i].v   - Gd_spins[i+4].v) ⋅ (p.α*L2 + p.β*L1) +
                  (Gd_spins[i+1].v - Gd_spins[i+5].v) ⋅ (p.α*L1 + p.β*L2) +
                  (Gd_spins[i+2].v - Gd_spins[i+6].v) ⋅ (-p.α*L1 + p.β*L2) +
                  (Gd_spins[i+3].v - Gd_spins[i+7].v) ⋅ (p.α*L2 - p.β*L1)) * pb
    end
    etot.e -= 8p.g * L1 ⋅ L2 * pb
    etot.e += pb^2/2 - pb * singleton(m, E_field).e
end

struct E_Pb_J <: System end

function Overseer.update(::E_Pb_J, m::AbstractLedger)
end

struct E_J <: System end

function Overseer.update(::E_J, m::AbstractLedger)
    s = m[Spin]
    bonds = m[JBond]
    etot = singleton(m, Etot)
    for b in bonds
        etot.e += b.J * s[b.e1].v ⋅ s[b.e2].v
    end
end