text:
#PcPcOnSq(k::Int,i::Int,j::Int) = pc_on_sq[k,int((i+1)*(i)/2)+(j)] # bug
PcPcOnSq(k::Int,i::Int,j::Int) = pc_on_sq[k,int((i-1)*(i)/2)+(j-1)+1]

function make_list( list0::Array{Int,1}, list1::Array{Int,1}, p::Board, gs::GameStatus)
    list2::Array{Int,1} = [0 for x=1:35]::Array{Int,1}
    nlist::Int = 15
    score::Int = 0
    sq_bk0::Int = p.kingposB
    sq_wk0::Int = p.kingposW
    sq_bk1::Int = 82 - p.kingposW
    sq_wk1::Int = 82 - p.kingposB

    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_pawn  + p.BlackPiecesInHands[MJFU]]
    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_lance + p.BlackPiecesInHands[MJKY]]
    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_knight+ p.BlackPiecesInHands[MJKE]]
    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_silver+ p.BlackPiecesInHands[MJGI]]
    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_gold  + p.BlackPiecesInHands[MJKI]]
    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_bishop+ p.BlackPiecesInHands[MJKA]]
    score += kkp[sq_bk0,sq_wk0, 1+ kkp_hand_rook  + p.BlackPiecesInHands[MJHI]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_pawn  + p.WhitePiecesInHands[MJFU]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_lance + p.WhitePiecesInHands[MJKY]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_knight+ p.WhitePiecesInHands[MJKE]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_silver+ p.WhitePiecesInHands[MJGI]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_gold  + p.WhitePiecesInHands[MJKI]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_bishop+ p.WhitePiecesInHands[MJKA]]
    score -= kkp[sq_bk1,sq_wk1, 1+ kkp_hand_rook  + p.WhitePiecesInHands[MJHI]]

    n2::Int = 1
    bb::BitBoard = p.bb[MJGOFU]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_pawn + sq + 1
        list2[n2]    = e_pawn + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_pawn + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJFU]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_pawn + sq + 1
        list2[n2]    = f_pawn + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_pawn + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGOKY]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_lance + sq + 1
        list2[n2]    = e_lance + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_lance + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJKY]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_lance + sq + 1
        list2[n2]    = f_lance + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_lance + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGOKE]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_knight + sq + 1
        list2[n2]    = e_knight + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_knight + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJKE]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_knight + sq + 1
        list2[n2]    = f_knight + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_knight + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGOGI]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_silver + sq + 1
        list2[n2]    = e_silver + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_silver + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJGI]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_silver + sq + 1
        list2[n2]    = f_silver + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_silver + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = (p.bb[MJGOKI]|p.bb[MJGOTO]|p.bb[MJGONY]|p.bb[MJGONK]|p.bb[MJGONG])
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_gold + sq + 1
        list2[n2]    = e_gold + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_gold + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = (p.bb[MJKI]|p.bb[MJTO]|p.bb[MJNY]|p.bb[MJNK]|p.bb[MJNG])
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_gold + sq + 1
        list2[n2]    = f_gold + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_gold + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGOKA]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_bishop + sq + 1
        list2[n2]    = e_bishop + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_bishop + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJKA]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_bishop + sq + 1
        list2[n2]    = f_bishop + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_bishop + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGOUM]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_horse + sq + 1
        list2[n2]    = e_horse + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_horse + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJUM]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_horse + sq + 1
        list2[n2]    = f_horse + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_horse + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGOHI]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_rook + sq + 1
        list2[n2]    = e_rook + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_rook + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJHI]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_rook + sq + 1
        list2[n2]    = f_rook + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_rook + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    n2 = 1
    bb = p.bb[MJGORY]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = f_dragon + sq + 1
        list2[n2]    = e_dragon + ((80 - sq)+1)
        score += kkp[sq_bk0,sq_wk0,kkp_dragon + sq + 1]
        nlist += 1
        n2    += 1
    end

    bb = p.bb[MJRY]
    while bb > uint128(0)
        sq::Int = trailing_zeros(bb)
        bb $= BitSet[sq+1]
        list0[nlist] = e_dragon + sq + 1
        list2[n2]    = f_dragon + ((80 - sq)+1)
        score += kkp[sq_bk1,sq_wk1,kkp_dragon + (80 - sq)+1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    if nlist > 53
        println("nlist is larger than 53.")
    end

    return nlist-1, score
end

function EvalBonanza(nextMove::Int, p::Board, gs::GameStatus)
    list0::Array{Int,1} = [0 for x=1:53]::Array{Int,1}
    list1::Array{Int,1} = [0 for x=1:53]::Array{Int,1}

    list0[ 1] = f_hand_pawn   + p.BlackPiecesInHands[MJFU] + 1
    list0[ 2] = e_hand_pawn   + p.WhitePiecesInHands[MJFU] + 1
    list0[ 3] = f_hand_lance  + p.BlackPiecesInHands[MJKY] + 1
    list0[ 4] = e_hand_lance  + p.WhitePiecesInHands[MJKY] + 1
    list0[ 5] = f_hand_knight + p.BlackPiecesInHands[MJKE] + 1
    list0[ 6] = e_hand_knight + p.WhitePiecesInHands[MJKE] + 1
    list0[ 7] = f_hand_silver + p.BlackPiecesInHands[MJGI] + 1
    list0[ 8] = e_hand_silver + p.WhitePiecesInHands[MJGI] + 1
    list0[ 9] = f_hand_gold   + p.BlackPiecesInHands[MJKI] + 1
    list0[10] = e_hand_gold   + p.WhitePiecesInHands[MJKI] + 1
    list0[11] = f_hand_bishop + p.BlackPiecesInHands[MJKA] + 1
    list0[12] = e_hand_bishop + p.WhitePiecesInHands[MJKA] + 1
    list0[13] = f_hand_rook   + p.BlackPiecesInHands[MJHI] + 1
    list0[14] = e_hand_rook   + p.WhitePiecesInHands[MJHI] + 1

    list1[ 1] = f_hand_pawn   + p.WhitePiecesInHands[MJFU] + 1
    list1[ 2] = e_hand_pawn   + p.BlackPiecesInHands[MJFU] + 1
    list1[ 3] = f_hand_lance  + p.WhitePiecesInHands[MJKY] + 1
    list1[ 4] = e_hand_lance  + p.BlackPiecesInHands[MJKY] + 1
    list1[ 5] = f_hand_knight + p.WhitePiecesInHands[MJKE] + 1
    list1[ 6] = e_hand_knight + p.BlackPiecesInHands[MJKE] + 1
    list1[ 7] = f_hand_silver + p.WhitePiecesInHands[MJGI] + 1
    list1[ 8] = e_hand_silver + p.BlackPiecesInHands[MJGI] + 1
    list1[ 9] = f_hand_gold   + p.WhitePiecesInHands[MJKI] + 1
    list1[10] = e_hand_gold   + p.BlackPiecesInHands[MJKI] + 1
    list1[11] = f_hand_bishop + p.WhitePiecesInHands[MJKA] + 1
    list1[12] = e_hand_bishop + p.BlackPiecesInHands[MJKA] + 1
    list1[13] = f_hand_rook   + p.WhitePiecesInHands[MJHI] + 1
    list1[14] = e_hand_rook   + p.BlackPiecesInHands[MJHI] + 1

    nlist::Int = 0
    score::Int = 0

    nlist, score = make_list( list0, list1, p, gs)
    #for i = 1:nlist
    #    if list0[i] >= fe_end
    #        println("list0[",i,"]=",list0[i])
    #    end
    #    if list1[i] >= fe_end
    #        println("list1[",i,"]=",list1[i])
    #    end
    #end
    sq_bk::Int = p.kingposB
    sq_wk::Int = 82 - p.kingposW

    sum::Int = 0

    for i = 1:nlist
        k0::Int = list0[i]
        k1::Int = list1[i]
        for j = 1:i-1
            l0::Int = list0[j]
            l1::Int = list1[j]
            #try
                sum += PcPcOnSq( sq_bk, k0, l0)
                sum -= PcPcOnSq( sq_wk, k1, l1)
            #catch
            #    println("sq_bk=",sq_bk,",sq_wk=",sq_wk)
            #    println("k0=",k0,",k1=",k1)
            #    println("l0=",l0,",l1=",l1)
            #    println("i,j=",i,",",j)
            #    quit()
            #end
        end
    end

    score += sum
    score += 32 * Eval( SENTE, p, gs)

    score = int(score/32)

    noise = (rand(Uint32) % 6) - 3
    score += noise

    score = (nextMove == SENTE) ? score: -score

    return score
end

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess


text:
using GpABC, DifferentialEquations, Distances, Plots, Distributions
#
# ABC settings
#
true_params =  [2.0, 1.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]
priors = [Uniform(0.2, 5.), Uniform(0.2, 5.), Uniform(10., 20.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.),
            Uniform(75., 125.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.)]
param_indices = [1,3,9]
#param_indices = [2,5,10]
priors = priors[param_indices]

#
# ODE solver settings
#
Tspan = (0.0, 10.0)
x0 = [3.0, 2.0, 1.0]
solver = RK4()
saveat = 0.1
#
# Returns the solution to the toy model as solved by DifferentialEquations
#
GeneReg = function(params::AbstractArray{Float64,1},
    Tspan::Tuple{Float64,Float64}, x0::AbstractArray{Float64,1},
    solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, saveat::Float64)

  if size(params,1) != 10
    throw(ArgumentError("GeneReg needs 10 parameters, $(size(params,1)) were provided"))
  end

  function ODE_3GeneReg(dx, x, par, t)
    dx[1] = par[1]/(1+par[7]*x[3]) - par[4]*x[1]
    dx[2] = par[2]*par[8]*x[1]/(1+par[8]*x[1]) - par[5]*x[2]
    dx[3] = par[3]*par[9]*x[1]*par[10]*x[2]./(1+par[9]*x[1])./(1+par[10]*x[2]) - par[6]*x[3]
  end

  prob = ODEProblem(ODE_3GeneReg, x0 ,Tspan, params)
  Obs = solve(prob, solver, saveat=saveat)

  return Array{Float64, 2}(Obs)
end

function simulator_function(var_params)
    params = copy(true_params)
    params[param_indices] .= var_params
    GeneReg(params, Tspan, x0, solver, saveat)
end

reference_data = simulator_function(true_params[param_indices])
plot(reference_data', xlabel="t", ylabel="C(t)", linewidth=2, labels=["u1(t)", "u2(t)", "u3(t)"])

# Rejection ABC
# Simulation
#
n_particles = 2000
threshold = 1.0
sim_result = SimulatedABCRejection(reference_data, simulator_function, priors, threshold, n_particles;
    write_progress=true)
plot(sim_result)

# Emulation
n_design_points = 500
emu_result = EmulatedABCRejection(reference_data, simulator_function, priors, threshold, n_particles, n_design_points;
    summary_statistic="keep_all", distance_function = Distances.euclidean, write_progress=true)
plot(emu_result)

# ABC-SMC
# Simulation
#
#threshold_schedule = [3.0, 2.0, 1.0, 0.5, 0.2];
threshold_schedule = [1.0, 0.8, 0.6, 0.4, 0.2];
population_colors=["#FF2F4E", "#D0001F", "#A20018", "#990017", "#800013"]

@time sim_abcsmc_res = SimulatedABCSMC(reference_data,
    simulator_function,
    priors,
    threshold_schedule,
    n_particles;
    write_progress = true)

plot(sim_abcsmc_res, population_colors=population_colors)

@time emu_abcsmc_res = EmulatedABCSMC(reference_data,
    simulator_function,
    priors,
    threshold_schedule,
    n_particles,
    n_design_points;
    batch_size=500,
    write_progress=true,
    emulator_retraining = PreviousPopulationThresholdRetraining(n_design_points, 100, 10),
    emulated_particle_selection = MeanVarEmulatedParticleSelection()
    )
plot(emu_abcsmc_res, population_colors=population_colors)

ext: jl


text:
# Addapted from the original https://github.com/JuliaGeo/GADM.jl (MIT Licensed)
# and stripped from all of it's dependencies (ArchGDAL repaced by the GMT GDAL functions).
# Expanded to also return all subregions of a particular administrative entity.
# To (partially) replace the role of the Tables.jl interface an option exists to print the
# names of all subregions of a parent administrative unit.

"""
    gadm(country, subregions...; children=false, names=false, children_raw=false, reportlevels=false)

Returns a GMTdataset for the requested country, or country subregion(s)

1. `country`: ISO 3166 Alpha 3 country code  .
2. subregions: Full official names in hierarchial order (provinces, districts, etc...).
   To know the names of all administrative children of parent, use the option `names`.
3. `children`: When true, function returns all subregions of parent.
4. `children_raw`: When true, function returns two variables -> parent, children, where children is a GDAL object
   E.g. when children is set to true and when querying just the country,
   second return parameter are the states/provinces. If `children` we return a Vector of GMTdataset with
   the polygons. If `children_raw` the second output is a GDAL object much like in GADM.jl (less the Tables.jl)
5. `names`: Return a string vector with all `children` names.
6. `reportlevels`: just report the number of administrative levels (including the country) and exit.

## Examples

```julia
# data of India's borders
data = gadm("IND")

# uttar -> the limits of the Uttar Pradesh state
uttar = gadm("IND", "Uttar Pradesh")

# uttar -> limits of all districts of the  Uttar Pradesh state
uttar = gadm("IND", "Uttar Pradesh", children=true)

# Names of all states of India
gadm("IND", names=true)
```
"""
function gadm(country, subregions...; children::Bool=false, names::Bool=false, children_raw::Bool=false, reportlevels::Bool=false)
	isvalidcode(country) || throw(ArgumentError("please provide standard ISO 3 country codes"))
	ressurectGDAL()			# Some previous GMT modules (or other shits) may have called GDALDestroyDriverManager()
	data_pato = country |> _download		# Downloads and extracts dataset of the given country code
	data = Gdal.read(data_pato)
	!isnothing(data) ? data : throw(ArgumentError("failed to read data from disk"))
	nlayers = Gdal.nlayer(data)
	reportlevels && return nlayers

	function _filterlayer(layer, level, value, all::Bool=false)
		filtered, key = Vector{GMT.Gdal.Feature}(undef,0), "NAME_$(level)"
		for row in layer
			field = Gdal.getfield(row, key)
			if all || occursin(lowercase(value), lowercase(field))
				push!(filtered, row)
			end
		end
		isempty(filtered) && throw(ArgumentError("could not find required region"))
		filtered
	end
	function _getnames(layer, level, parent)
		# Get the names of all direct descendents of 'parent'
		filtered, key, keyP = Vector{String}(undef,0), "NAME_$(level)", "NAME_$(level-1)"
		for row in layer
			field  = Gdal.getfield(row, key)
			if (parent == "")		# Easier, all descendents are what we are looking for
				push!(filtered, field)
			else					# Filter only the descendents of a specific parent
				fieldP = Gdal.getfield(row, keyP)
				contains(lowercase(fieldP), parent) && push!(filtered, field)
			end
		end
		filtered
	end

	function _getlayer(plevel)
		# Get layer of the desired `level` from the `data`.
		nlayers = Gdal.nlayer(data)
		for l = 0:nlayers - 1
			layer = Gdal.getlayer(data, l)
			lname = Gdal.getname(layer)
			llevel = last(split(lname, "_"))
			string(plevel) == llevel && return layer
		end
		error("Asked data for a level ($(plevel+1)) that is lower than lowest data level ($(nlayers))")
	end

	function _get_attrib(feature)
		n = Gdal.nfield(feature)
		attrib = Dict{String, String}()
		[attrib[Gdal.getname(Gdal.getfielddefn(feature, i))] = string(Gdal.getfield(feature, i)) for i = 0:n-1]
		attrib
	end

	function _get_polygs(gdfeature)
		D = gd2gmt(getgeom(gdfeature[1], 0),"")
		att = _get_attrib(gdfeature[1])
		for k = 1:length(D)  D[k].attrib = att;  D[k].colnames = ["Lon", "Lat"]  end
		D[1].attrib = _get_attrib(gdfeature[1])
		((prj = getproj(Gdal.getlayer(data, 0))) != C_NULL) && (D[1].proj4 = toPROJ4(prj))
		for n = 2:length(gdfeature)
			_D = gd2gmt(getgeom(gdfeature[n], 0),"")
			att = _get_attrib(gdfeature[n])
			for k = 1:length(_D)  _D[k].attrib = att;  _D[k].colnames = ["Lon", "Lat"]  end
			append!(D, _D)
		end
		set_dsBB!(D, false)		# Compute only the global BB. The other were computed aready
		D
	end

	# p -> parent, is the requested region
	plevel = length(subregions)
	plevel >= nlayers && throw(ArgumentError("more subregions required than in data"))
	pname = isempty(subregions) ? "" : last(subregions)
	player = _getlayer(plevel)

	if (!children && !children_raw && !names) || (!names && nlayers == plevel + 1)	# Last case is when we have no more levels
		p = _filterlayer(player, plevel, pname, iszero(plevel))
		return !children_raw ? _get_polygs(p) : (_get_polygs(p), nothing)
	end

	# c -> children, is the region 1 level lower than p
	clevel = plevel + 1
	clayer = _getlayer(names || children || children_raw ? clevel : plevel)

	(names) && return _getnames(clayer, clevel, (plevel == 0) ? "" : lowercase(subregions[end]))

	c = _filterlayer(clayer, plevel, pname, iszero(plevel))

	if (children_raw)		# This is close (aside from the GMTdatset vs GDAL features) to the output of the GADM.jl
		p = _filterlayer(player, plevel, pname, iszero(plevel))
		return _get_polygs(p), c
	else
		return _get_polygs(c)
	end
end

# ------------------------------------------------------------------------------------------------------
# Tells whether or not `str` is a valid ISO 3166 Alpha 3 country code. Valid code examples are "IND", "USA", "BRA".
isvalidcode(str) = match(r"\b[A-Z]{3}\b", str) !== nothing

# ------------------------------------------------------------------------------------------------------
function _download(country)
	# Downloads data (for the first call) for `country` and returns its full name.
	ID, name_zip = "gadm36_$(country).gpkg", "gadm36_$(country)_gpkg.zip"
	cache = joinpath(readlines(`gmt --show-userdir`)[1], "cache")
	if !isdir(cache)
		((pato = mkdir(cache)) == "") && error("Failed to create the 'cache' dir where download file would be stored")
	end
	fname = joinpath(cache, ID)
	isfile(fname) && return fname

	println("Downloading geographic data for country $country provided by the https://gadm.org project. It may take a while.")
	println("The file $ID (after uncompressing) will be stored in $cache")
	dlfile = download("https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/$(name_zip)", joinpath(cache, name_zip))
	curr_pato = pwd();	cd(cache)
	@static Sys.iswindows() ? run(`tar -xf $dlfile`) : run(`unzip $dlfile`)
	rm(dlfile)
	cd(curr_pato)
	return fname
end


text:
using DynamicalBilliards, InteractiveChaos, Makie

N = 100
colors = :dense
colors = [Makie.RGBAf0(i/N, 0, 1 - i/N, 0.25) for i in 1:N]

# Uncomment any of the following to get the billiard you want:
# bd = billiard_stadium(1.0f0, 1.0f0)
# bd = Billiard(bd..., Disk(SVector(0.5f0, 0.5f0), 0.2f0))
bd = billiard_mushroom(1.0f0, 0.2f0, 1.0f0, 0.0f0; door = false)
# bd = billiard_hexagonal_sinai(0.5f0, 1.0f0)
# bd = billiard_sinai(0.25f0, 1f0, 1f0)
# bd = Billiard(Antidot(Float32[0, 0], 0.5f0, false))
# bd, = billiard_logo(T = Float32)

ps = [MagneticParticle(1, 0.6f0 + 0.0005f0*i, 0, 1f0) for i in 1:N]
ps = [Particle(1, 0.6f0 + 0.0005f0*i, 0) for i in 1:N]
ps = particlebeam(0.8, 0.6, π/4, N, 0.01, nothing, Float32)

interactive_billiard(bd, 1f0, tail = 1000)

# %% also interact with boundary map
interactive_billiard_bmap(bd)

# %% make it a video
billiard_video("billiard.mp4", bd, 1f0; plot_particles=true)

# %% 3b1b style video:
# Colors of 3b1b
BLUE = "#7BC3DC"
BROWN = "#8D6238"
colors = [BLUE, BROWN]
# Overwrite default color of obstacles to white (to fit with black)
InteractiveChaos.obcolor(::Obstacle) = RGBf0(1,1,1)
bd = billiard_stadium(1.0f0, 1.0f0)
billiard_video("3b1billiard.mp4", bd, 1.0, 0.6, 0;
frames = 600,
speed = 8, backgroundcolor = :black, colors = colors)

# %% static plot of boundary map and billiard (several particles, same color)
ps = [randominside(bd) for i in 1:N]
scene = billiard_bmap_plot(bd, ps; colors = colors, backgroundcolor = RGBf0(1,1,1))
Makie.save("static_billiard_plot.png", scene, dpi = 600)


text:
#!/usr/bin/env julia

# Runs tests of the package

using Pkg

Pkg.activate(".")
Pkg.test("JuliaPkgsList")


text:
"""
    execute_sts_intraday_api_call(requestDictionary::Dict{String,Any};
        logger::Union{Nothing,AbstractLogger} = nothing) -> PSResult

The best method description ever will go here. It will be amazing.
"""
function execute_sts_intraday_api_call(requestDictionary::Dict{String,Any};
    logger::Union{Nothing,AbstractLogger} = nothing)::PSResult

    # initialize -
    api_call_url_string = ""
    required_api_keys = [
        "symbol"        ;   # ticker symbol -
        "function"      ;   # what function are we doing?
        "interval"      ;   # what interval do we need?
        "apikey"        ;   # the API key
        "adjusted"      ;   # whether data is adjsuted or raw (default: true)
        "outputsize"    ;   # Strings compact (default) and full are accepted: compact returns only the latest 100 data points in the intraday time series; full returns the full-length intraday time series.
        "datatype"      ;   # what type of data is returned from API. Default: json
    ]

    try

        # ok, before we get too far, we need to check for the optional data in the requestDictionary -
        if (haskey(requestDictionary,"adjusted") == false)
            requestDictionary["adjusted"] = "true"
        end

        # outputsize -
        if (haskey(requestDictionary,"outputsize") == false)
            requestDictionary["outputsize"] = "compact"
        end

        # datatype -
        if (haskey(requestDictionary,"datatype") == false)
            requestDictionary["datatype"] = "json"
        end

        # get stuff from the requestDictionary and build the url used in the call -
        api_call_url_string = "$(alphavantage_api_url_string)?"
        for (index,key) in enumerate(required_api_keys)
            value = requestDictionary[key]
            if (index == 1)
                api_call_url_string*="$(key)=$(value)"
            else
                api_call_url_string*="&$(key)=$(value)"
            end
        end

        # make the API call -
        api_call_result = http_get_call_with_url(api_call_url_string)
        response_body_string = checkresult(api_call_result)

        # ok, we process the response body differently depending upon json -or- csv type -
        data_type = requestDictionary["datatype"]
        if (Symbol(data_type) == :csv)
            return process_raw_csv_api_data(response_body_string)
        elseif (Symbol(data_type) == :json)
            return process_raw_json_sts_intraday_data(response_body_string)
        else
            throw(PSError("unsupported datatype in request dictionary"))
        end
    catch error
        return PSResult(error)
    end
end


text:
module OpenSecrets

using DataFrames

include("latin1buffer.jl")
include("opensecretsbuffer.jl")
include("data_detection.jl")
include("data_loading.jl")

global _campaign_finance_data_sources = detect_data_sources()

end


text:
################################################################################
### Construct table and write to disk

data_infile    = "my_test_table_in.csv"
data_outfile   = "my_test_table_out.csv"
issues_infile  = "my_test_table_input_issues.tsv"
issues_outfile = "my_test_table_output_issues.tsv"

tbl = DataFrame(
    patientid = [1, 2, 3, 4],
    age       = [11, 22, 33, 444],
    dose      = ["small", "medium", "large", "medium"],
    fever     = [false, true, true, false]
)

CSV.write(data_infile, tbl; delim=',')
tbl = nothing

################################################################################
### Compare DataFrame to Schema

# Schema
patientid = ColumnSchema(:patientid, "Patient ID",  UInt,   !CATEGORICAL, REQUIRED,  UNIQUE, UInt)
age       = ColumnSchema(:age,       "Age (years)", Int,    !CATEGORICAL, REQUIRED, !UNIQUE, Int)
dose      = ColumnSchema(:dose,      "Dose size",   String,  CATEGORICAL, REQUIRED, !UNIQUE, ["small", "medium", "large"])
fever     = ColumnSchema(:fever,     "Had fever",   Bool,    CATEGORICAL, REQUIRED, !UNIQUE, Bool)
ts        = TableSchema(:mytable, "My table", [patientid, age, dose, fever], [:patientid])

# Compare data to schema
compare(ts, data_infile; output_data_file=data_outfile, input_issues_file=issues_infile, output_issues_file=issues_outfile)
issues_in = DataFrame(CSV.File(issues_infile; delim='\t'))
@test size(issues_in, 1) == 0

# Modify schema: Forbid the age column from having values of 120 or above
age.validvalues = 0:120

# Compare again
compare(ts, data_infile; output_data_file=data_outfile, input_issues_file=issues_infile, output_issues_file=issues_outfile)
outdata    = DataFrame(CSV.File(data_outfile))
issues_in  = DataFrame(CSV.File(issues_infile; delim='\t'))
issues_out = DataFrame(CSV.File(issues_outfile; delim='\t'))
@test size(issues_in, 1)  == 1
@test size(issues_out, 1) == 1
@test ismissing(outdata[4, :age]) == true

# Fix input data
indata = DataFrame(CSV.File(data_infile))
indata[4, :age] = 44
CSV.write(data_infile, indata; delim=',')
compare(ts, data_infile; output_data_file=data_outfile, input_issues_file=issues_infile, output_issues_file=issues_outfile)
issues_in = DataFrame(CSV.File(issues_infile; delim='\t'))
@test size(issues_in, 1) == 0

################################################################################
# Test intra-row constraints
function test_row_constraints()
    filename = joinpath(dirname(pathof(Schemata)), "..", "test/schemata/row_constraints.yaml")
    schema   = readschema(filename)
    indata   = DataFrame(
                  patientid = UInt.([1,2,3]),
                  dob=Date.(["1992-10-01", "1988-03-23", "1983-11-18"]),
                  date_of_marriage=[Date("2015-09-13"), missing, Date("1981-11-01")])
    CSV.write(data_infile, indata; delim=',')
    compare(schema.tables[:dates], data_infile; output_data_file=data_outfile, input_issues_file=issues_infile, output_issues_file=issues_outfile)
end
test_row_constraints()
issues_in = DataFrame(CSV.File(issues_infile))
@test size(issues_in, 1) == 1

################################################################################
# Clean up

rm(data_infile)
rm(data_outfile)
rm(issues_infile)
rm(issues_outfile)


text:
"""
    make_folds(n, k, l)

Generate `k` non-overlapping folds.

# Arguments

- n = Total number of observations to split into folds.
- k = Number of folds to create. Defaults to 10. If `k=1`, then all the data
  (along this dimension) will be included in each fold.
- k2 = If `k=1`, then all the data (along this dimension) will be included in
  each fold. `k2` specifies how many folds there are. Defaults to `k`, which
  is kind of silly, but there needs to be a placeholder.

# Value

1d array of length `k` of arrays of indices

"""
function make_folds(n::Int64, k::Int64=10, k2::Int64=k)

    if k > 1
        # When k > 1, run the Kfolds function from MLBase to create folds of length ~ n*(1-1/k)
        return collect(Array{Int64,1}, Kfold(n, k))
    elseif k == 1
        # When k = 1, repeat all indices 1:n k2 times
        return repeat([collect(1:n)], inner=k2)
    else
        error("k must be strictly positive.")
    end
end

"""
    make_folds_conds(conds, k, prop)

Generate `k` folds for a set of conditions, making sure each level of each
condition is represented in each fold.

# Arguments

- conds = 1d array of conditions (strings)
- k = Number of folds to create. Defaults to 10.
- prop = Proportion of each condition level's replicates to include in each
  fold. Defaults to 1/k. Each fold will contain at least one replicate of
  each condition level.

# Value

1d array of length `k` of arrays of indices

"""
function make_folds_conds(conds::AbstractArray{String,1},
                          k::Int64=10, prop::Float64=1/k)

    if k == 1
        print("k is set to 1; may result in unexpected behavior")
    elseif k < 1
        error("k must be strictly positive.")
    end

    # All indices to split into folds
    idx = 1:length(conds)
    # 2d boolean array to store which observations to keep in each fold
    boolFolds = falses(length(conds), k)

    # Iterate through all the conditions
    for cond in unique(conds)
        # For each condition, calculate the number of replicates to sample
        numSamp = convert(Int64, ceil(sum(cond .== conds) * prop))
        # Iterate through all the folds
        for j in 1:k
            # For each fold, randomly keep a sample of each condition
            boolFolds[rand(idx[cond .== conds], numSamp) ,j] .= true
        end
    end

    # Initialze 1d array of arrays of indices
    folds = fill(convert(Array{Int64}, []), k)
    # Convert each column of booleans from boolFolds into indices
    for j in 1:k
        folds[j] = findall(boolFolds[:,j])
    end

    return folds
end

"""
    findnotin(a, b)

Returns elements of `b` that are not present in `a`.

# Arguments

- a = 1d array of integers
- b = 1d array of integers

# Value

1d array of integers

"""
function findnotin(a::AbstractArray{Int64,1}, b::AbstractArray{Int64,1})

    if a != b
        return [x for x=b if !(x in a)]
    else
        return a
    end
end

"""
    calc_mse(MLMNets, data, lambdas, rowFolds, colFolds)

Calculates test MSE for each of the CV folds for each lambda.

# Arguments

- MLMNets = 1d array of Mlmnet objects resulting from running cross validation
- data = RawData object used to generate MLMNets
- lambdasL1 = 1d array of floats consisting of lambda penalties used to
  generate MLMNets
- lambdasL2 = 1d array of floats consisting of lambda penalties used to
  generate MLMNets
- rowFolds = 1d array of arrays containing booleans for the row folds
- colFolds = 1d array of arrays containing booleans for the column folds

# Value

2d array of floats with dimensions equal to the number of lambdas by the
number of folds.

"""
function calc_mse(MLMNets::AbstractArray{Mlmnet,1}, data::RawData,
                  lambdas::AbstractArray{Float64,1},
                  alphas::AbstractArray{Float64,1},
                  rowFolds::Array{Array{Int64,1},1},
                  colFolds::Array{Array{Int64,1},1})

    # Number of folds
    nFolds = length(MLMNets)

    # For each fold, generate a RawData object corresponding to the holdout set
    holdoutFolds = Array{RawData}(undef, nFolds)
    for i in 1:nFolds
        holdoutFolds[i] =
          RawData(Response(get_Y(data)[findnotin(rowFolds[i],
                                                 collect(1:data.n)),
                                       findnotin(colFolds[i],
                                                 collect(1:data.m))]),
                  Predictors(get_X(data)[findnotin(rowFolds[i],
                                                   collect(1:data.n)),:],
                             get_Z(data)[findnotin(colFolds[i],
                                                   collect(1:data.m)),:],
                             data.predictors.hasXIntercept,
                             data.predictors.hasZIntercept))
    end

    # Initialize array to store test MSEs
    mse = Array{Float64}(undef, length(lambdas), length(alphas), nFolds)

    # Iterate through all folds and lambdas to calculate test MSE
    for j in 1:nFolds
        # Residuals for all lambdasL1 and lambdasL2 at fold j
        resids = resid(MLMNets[j], holdoutFolds[j])
        for i in 1:length(lambdas), k in 1:length(alphas)
            # MSE for (lambdaL1 i, lambdaL2 k) and fold j
            mse[i,k,j] = mean(resids[:,:,i,k].^2)
        end
    end

    return mse
end

"""
    calc_prop_zero(MLMNets, lambdas, alphas, dig)

Calculates proportion of zero interaction coefficients for each of the CV
folds for each lambda.

# Arguments

- MLMNets = 1d array of MlmnetDeprecated objects resulting from running cross validation
- lambdas = 1d array of floats consisting of lambda penalties used to
  generate MLMNets

# Keyword arguments

- dig = integer; digits of precision for zero coefficients. Defaults to 12.

# Value

2d array of floats with dimensions equal to the number of lambdas by the
number of folds.

"""
function calc_prop_zero(MLMNets::AbstractArray{Mlmnet,1},
                        lambdas::AbstractArray{Float64,1},
                        alphas::AbstractArray{Float64,1};
                        dig::Int64=12)

    # Number of folds
    nFolds = length(MLMNets)

    # Initialize array to store proportions of zero interactions
    propZero = Array{Float64}(undef, length(lambdas), length(alphas), nFolds)

    # Boolean arrays used to subset coefficients for interactions only
    xIdx = trues(MLMNets[1].data.p)
    xIdx[1] = MLMNets[1].data.predictors.hasXIntercept==false
    zIdx = trues(MLMNets[1].data.q)
    zIdx[1] = MLMNets[1].data.predictors.hasZIntercept==false

    # Iterate through all folds and lambdas to calculate proportion of zero
    # interaction coefficients
    for j in 1:nFolds
        for k in 1:length(alphas),i in 1:length(lambdas)
            # Proportion of zero interaction coefficients for lambda i and
            # fold j
            propZero[i,k,j] = mean(round.(coef(MLMNets[j])[xIdx,zIdx,i,k],
                                        digits=dig) .== 0)
        end
    end

    return propZero
end


text:
using Documenter, LinearForms

makedocs(clean=false)


text:
using Documenter, Potree

Documenter.makedocs(
    format = Documenter.HTML(
		prettyurls = get(ENV, "CI", nothing) == "true"
	),
sitename = "Potree.jl",
pages = [
    "Home" => "index.md",
    "Load" => [
        "bin.jl" => "bin.md",
        "hierarchy.jl" => "hierarchy.md",
        "json.jl" => "json.md",
        "las.jl" => "las.md"
        ],
    "struct.jl" => "struct.md",
    "utilities.jl" => "utilities.md"
    ]
)


text:
"""
A Julia package to send SMS (Short Message Service) using
[BulkSMS API](http://www.bulksms.com/products/sms-api.htm)
(more exactly [HTTP to SMS API](http://www.bulksms.com/products/api/http-to-sms.htm))
"""
module BulkSMS

    export BulkSMSClient, send, shorten, multiple

    using HTTP

    """
        ActionWhenLong

    Enumeration that defines how long message treated

    - `shorten`: if message is too long, it will be shorten.
    - `multiple`: if message is too long, several messages will be sent
    """
    @enum ActionWhenLong begin
        shorten = 1
        multiple = 2
    end

    _DEFAULT_BASE_URL = "http://bulksms.vsms.net:5567"
    _DEFAULT_MAX_MESSAGE_LEN = 160

    """
        BulkSMSClientException(s)

    BulkSMS client exception.
    """
    struct BulkSMSClientException <: Exception
        s::AbstractString
    end

    """
        BulkSMSResponse(response)

    BulkSMS response (from `response::Requests.Response`).
    """
    struct BulkSMSResponse
        status_code::Int
        status_string::AbstractString
        id::Int
        function BulkSMSResponse(response::HTTP.Messages.Response)
            s = String(response.body)
            statusCode, statusString, Id = split(s, "|")
            new(parse(Int64, statusCode), statusString, parse(Int64, Id))
        end
    end

    """
        BulkSMSClient(msisdn, username, password;
                base_url=_DEFAULT_BASE_URL,
                max_message_len=_DEFAULT_MAX_MESSAGE_LEN)

    Create a client that can connect to BulkSMS HTTP to SMS API.
    """
    struct BulkSMSClient
        msisdn::AbstractString
        username::AbstractString
        password::AbstractString

        base_url::AbstractString
        max_message_len::Int

        function BulkSMSClient(msisdn, username, password;
                base_url=_DEFAULT_BASE_URL,
                max_message_len=_DEFAULT_MAX_MESSAGE_LEN)
            new(msisdn, username, password, base_url, max_message_len)
        end
    end

    """
        _crop(msg, max_len) -> cropped_msg
    Crop a `message` to `max_len`.
    """
    function _crop(msg, max_len)
        if max_len > 0 && length(msg) > max_len
            msg[1:max_len-3] * "..."
        else
            msg
        end
    end

    """
        _send(client, message_text, msisdn)

    Send a text message `message_text` using `client`.

    This is a "low level" function. You should probably use `send(...)` insted of `_send(...)`
    """
    function _send(client::BulkSMSClient, message_text::AbstractString, msisdn::AbstractString)
        endpoint = "/eapi/submission/send_sms/2/2.0"

        url = client.base_url * endpoint

        params = Dict(
            "username" => client.username,
            "password" => client.password,
            "message" => message_text,
            "msisdn" => msisdn
        )

        @debug "GET $url with query=$params"
        raw_response = HTTP.request("GET", url; query = params)
        @debug raw_response

        if raw_response.status != 200
            throw(BulkSMSClientException("HTTP status code != 200"))
        end

        response = BulkSMSResponse(raw_response)

        if response.status_code != 0
            throw(BulkSMSClientException("BulkSMS status code != 0"))
        end

        response
    end

    """
        send(client, message_text; msisdn=nothing, action_when_long::ActionWhenLong=shorten)

    Send a text message `message_text` to BulkSMS using HTTP to SMS API with BulkSMSClient `client`.
    """
    function send(client::BulkSMSClient, message_text::AbstractString; msisdn=nothing, action_when_long::ActionWhenLong=shorten)

        if msisdn === nothing
            msisdn = client.msisdn
        end

        response = nothing

        if action_when_long == shorten
            shorten_message_text = _crop(message_text, client.max_message_len)
            response = _send(client, shorten_message_text, msisdn)

        # elseif action_when_long == multiple

        else
            throw(BulkSMSClientException("action_when_long=$action_when_long is not a support action"))
        end

        response
    end

end # module


text:
# Bayesian blocks
# Implementation by Michael P.H. Stumpf and T. Chan
# Based on the Python Code of Jake Vanderplas.

# This version implements Bayesian blocks for histograms, where the
# data are sorted, then treated as event data (see Scargle 2012).
# Defaults as suggested in Scargle 2012 are used.

# References:
# Scargle 2012: http://adsabs.harvard.edu/abs/2012arXiv1207.5578S
# Python implementation: https://github.com/astroML/astroML/blob/master/astroML/density_estimation/bayesian_blocks.py

type DiscretizeBayesianBlocks <: DiscretizationAlgorithm
end

function invVector(input::Array{Float64,})
	len = length(input)
	output = zeros(len)
	for i in 1:len
		output[len+1-i]=input[i]
	end
	return output
end

function binedges{N<:AbstractFloat}(alg::DiscretizeBayesianBlocks, data::AbstractArray{N})

	unique_data = unique(data)
	unique_data = sort(unique_data)

	n = length(unique_data) # Number of observations

	edges = zeros(n+1)
	edges[1] = unique_data[1]
	for i in 1 : (n-1)
		edges[i+1] = 0.5 * (unique_data[i]+unique_data[i+1])
	end
	edges[end] = unique_data[end]
	block_length = unique_data[end] - edges

	if length(unique_data) == length(data)
		nn_vec = ones(length(data))
	else
		nn_vec = convert(Array{Float64}, [length(findin(data, v)) for v in unique_data])
	end
	best = zeros(n)
	last = zeros(Int64,n)

	for K in 1 : n
		width = block_length[1 : K] - block_length[K+1]
		count_vec = invVector(cumsum(invVector(nn_vec[1 : K])))

		# Fitness function (eq. 19 from Scargle 2012)
		fit_vec = count_vec .* (log(count_vec) - log(width))
		# Prior (eq. 21 from Scargle 2012)
		fit_vec -= 4 - log(73.53 * 0.05 * ((K)^-0.478))
		fit_vec[2:end] += best[1 : K-1]

		i_max = indmax(fit_vec)
		last[K] = i_max
		best[K] = fit_vec[i_max]
	end

	change_points = zeros(Int64,n)
	i_cp = n+1
	ind = n+1
	while true
		i_cp -= 1
		change_points[i_cp] = ind
		if ind == 1
			break
		end
		ind = last[ind-1]
	end
	change_points = change_points[i_cp : end]
	edges[change_points]

end
function binedges{N<:Integer}(alg::DiscretizeBayesianBlocks, data::AbstractArray{N})

	data = convert(Array{typeof(0.0)}, data)
	binedges(alg, data)

end


text:
module WaitAndSee

using Compat, Compat.SparseArrays

using StructDualDynProg
using StochOptInterface
const SOI = StochOptInterface

using MathProgBase

struct Algorithm{S<:MathProgBase.AbstractMathProgSolver} <: SOI.AbstractAlgorithm
    solver::S
    K::Int
end
Algorithm(solver::MathProgBase.AbstractMathProgSolver, K=25) = Algorithm(solver, K)

mutable struct WaitAndSeePath{NodeT}
    node::NodeT
    nlds::Vector{StructDualDynProg.StructProg.NLDS}
    z::Float64
    proba::Float64
    K::Int
end

function meanstdpaths(paths::Vector{WaitAndSeePath}, totalK)
    z = Float64[x.z for x in paths]
    proba = Float64[x.proba for x in paths]
    npaths = Int[x.K for x in paths]
    return StructDualDynProg.SDDP.meanstdpaths(z, proba, npaths, totalK)
end

function SOI.optimize!(sp::SOI.AbstractStochasticProgram, algo::Algorithm,
                       stopcrit::SOI.AbstractStoppingCriterion=SOI.IterLimit(0), # Not used
                       verbose::Int=0)
    master = SOI.get(sp, SOI.MasterNode())
    num_stages = SOI.get(sp, SOI.NumberOfStages())
    paths = WaitAndSeePath[WaitAndSeePath(master, StructDualDynProg.StructProg.NLDS[StructDualDynProg.StructProg.nodedata(sp, master).nlds], .0, 1., algo.K)]
    for t in 2:num_stages
        newpaths = WaitAndSeePath[]
        for path in paths
            if isempty(SOI.get(sp, SOI.OutTransitions(), path.node))
                push!(newpaths, path)
            else
                npaths = StructDualDynProg.SDDP.samplepaths(StructDualDynProg.SDDP.ProbaPathSampler(true), sp, path.node, path.K, t, num_stages)
                childs = algo.K == -1 ? (1:length(SOI.get(sp, SOI.OutTransitions(), path.node))) : findall(npaths .> 0)
                for (i, tr) in enumerate(SOI.get(sp, SOI.OutTransitions(), path.node))
                    if algo.K == -1 || npaths[i] > 0
                        push!(newpaths, WaitAndSeePath(SOI.get(sp, SOI.Target(), tr), [path.nlds; StructDualDynProg.StructProg.nodedata(sp, SOI.get(sp, SOI.Target(), tr)).nlds], path.z, path.proba * SOI.get(sp, SOI.Probability(), tr), npaths[i]))
                    end
                end
            end
        end
        paths = newpaths
    end

    sumz = 0
    newpaths = WaitAndSeePath[]
    for path in paths
        nvars = cumsum(Int[nlds.nx for nlds in path.nlds])
        ncons = cumsum(Int[nlds.nπ for nlds in path.nlds])
        # b - Ax in K_1
        # x \in K_2
        A = spzeros(ncons[end], nvars[end])
        bs = Vector{Vector}(undef, length(path.nlds))
        Ks = []
        C = []
        c = similar(StructDualDynProg.StructProg.nodedata(sp, master).nlds.c, nvars[end])
        for i in 1:length(path.nlds)
            offsetnvars = i == 1 ? 0 : nvars[i-1]
            offsetncons = i == 1 ? 0 : ncons[i-1]
            nlds = path.nlds[i]
            c[offsetnvars+1:nvars[i]] = nlds.c
            bs[i] = nlds.h
            A[offsetncons+1:ncons[i], offsetnvars+1:nvars[i]] = nlds.W
            if i > 1
                offset = i == 2 ? 0 : nvars[i-2]
                A[offsetncons+1:ncons[i], offset+1:nvars[i-1]] = nlds.T
            end
            push!(Ks, nlds.K)
            append!(C, [(cone, offsetnvars .+ idx) for (cone, idx) in nlds.C])
        end
        model = MathProgBase.LinearQuadraticModel(algo.solver)
        StructDualDynProg.StructProg._load!(model, c, A, bs, Ks, C)
        MathProgBase.optimize!(model)
        status = MathProgBase.status(model)
        if status == :Error
            error("The solver reported an error")
        elseif status == :UserLimit
            error("The solver reached iteration limit or timed out")
        elseif status == :Infeasible
            error("The problem is infeasible for some realization of the uncertainty")
        elseif status == :Unbounded
            error("The problem is unbounded for some realization of the uncertainty")
        else
            @assert status == :Optimal
            objval = MathProgBase.getobjval(model)
            push!(newpaths, WaitAndSeePath(path.node, path.nlds, objval, path.proba, path.K))
        end
    end
    meanstdpaths(newpaths, algo.K)
end

end


text:
using Query
import Mortchartgen
include("amplt.jl")
if !(isdefined(:frames))
	frames = Mortchartgen.load_frames()
end
deaths = frames[:deaths]
pop = frames[:pop]
pop_csy(country, sex, year) = @from i in pop begin
         @where i.Sex==sex && i.Country==country && i.Year==year
         @select {i.variable,i.value}
         @collect DataFrame
       end
deaths_cscy(country, sex, cause, year) = @from i in deaths begin
         @where i.Sex==sex && i.Country==country && i.Cause==cause && i.Year==year
         @select {i.variable,i.value}
         @collect DataFrame
       end

function caframe(country, sex, cause, year)
	popf = pop_csy(country, sex, year)
	alldeathsf = deaths_cscy(country, sex, "all", year)
	deathsf = deaths_cscy(country, sex, cause, year)
	DataFrame(age=vcat(0:4,5:5:95)
		, pop=popf[:value][2:25]
		, acd=alldeathsf[:value][2:25]
		, cd=deathsf[:value][2:25]
		, cc=fill(0,24)
		)
end

function caframe_sm(caf, smf)
	caf[:aints] = collect(1:24)
	caf[:cd_prop] = zfrac.(caf[:cd], caf[:acd])
	smf[:aints] = [1:5;map(x->fld(x,5)+5, 5:99);fill(24,6)]
	caf_smf = join(caf, smf, on = :aints)
	cd_sm = caf_smf[:dths].*caf_smf[:cd_prop]
	DataFrame(age=smf[:age]
		, pop=smf[:pop]
		, acd=smf[:dths]
		, cd=cd_sm
		, cc=fill(0,106)
		)
end

function chfac(country, sex, cause, year, cf_cd, cf_nocd, smframe = NA)
	caf_0 = caframe(country, sex, cause, year)
	if !(isna(smframe))
		caf = caframe_sm(caf_0, smframe)
	else
		caf = caf_0
	end
	lt = periodlifetable(caf, sex)
	rc = ratechange(caf, cf_cd, 1, cf_nocd)
	lt_rc = periodlifetable(rc, sex)
	caf_amp = amplt(caf, sex, "mort")
	caf_amp_rc = amplt(rc, sex, "mort")
	Dict(:lt => lt, :lt_rc => lt_rc, :changes => changes(caf_amp_rc, caf_amp))
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
export eye

"""
`eye(T,n)` returns an `n`-by-`n` identity matrix whose entries
are numbers of type `T` (`Float64` by default).
"""
function eye(T, n)
    A = zeros(T, n, n)
    for j = 1:n
        @inbounds A[j, j] = one(T)
    end
    return A
end

eye(n) = eye(Float64, n)


text:
# mbe_demo_sitl.jl
#
# This script runs the C implementation of the control law in the simulation.

# Include our custom models.
include("models.jl")

# We work inside a module because it's more flexible than the main workspace.
module MBEDemoSITL

#############
# Resources #
#############

# Use/import the modules we'll need.
using Eratosthenes # The simulation engine (setup, simulate, mc)
using ..MBEModels  # The MBEModels module included above (and "above" this module)
import YAML        # For loading the YAML file directly
import HDF5        # For loading the log file
using Plots        # For plotting our results
using LinearAlgebra

# Create a directory for the log file.
if !isdir("out"); mkdir("out"); end

########################
# Software-in-the-Loop #
########################

# Specify the name of the YAML file that defines our scenario.
yaml_file = joinpath(@__DIR__, "underactuated.yaml")

# Load the YAML file and change the model name to use the SITL version.
# Also, change the log name.
yaml = YAML.load_file(yaml_file)
yaml["vehicles"][1]["computers"][1]["software"][1]["model"] = "ReducedEffortUnderactuatedControllerSITL"
yaml["sim"]["log_file"] = "out/underactuated-sitl.h5"

# The parameters are now stored under constants/parameters, so move them.
yaml["vehicles"][1]["computers"][1]["software"][1]["constants"] = Dict(
    "parameters"=>yaml["vehicles"][1]["computers"][1]["software"][1]["constants"]
);

# Set up the senario and run it.
scenario = setup(Scenario(), yaml, @__MODULE__)
scenario = simulate(scenario)

#########
# Plots #
#########

# Read some data from the logs.
t, q_BI, ω_BI_B = HDF5.h5open(scenario.sim.log_file, "r") do logs
    (dropdims(read(logs, "/smallsat/body/state/time"), dims=1),
     read(logs, "/smallsat/body/state/data/q_BI"),
     read(logs, "/smallsat/body/state/data/ω_BI_B"))
end

# Choose a friendly plotting package.
pyplot()
# plotlyjs()

# Define each plot that we'll need.
display(plot(t, transpose(q_BI),
             label  = ["q1" "q2" "q3" "q4"],
             xlabel = "Time (s)",
             ylabel = "Attitude Quaternion"))
display(plot(t, 180/π * transpose(ω_BI_B),
             label  = ["ω1" "ω2" "ω3"],
             xlabel = "Time (s)",
             ylabel = "Rotation Rate (deg/s)"))

# Compare to the baseline.
(q_BI_baseline, ω_BI_B_baseline) = HDF5.h5open("out/underactuated.h5", "r") do logs
    (read(logs, "/smallsat/body/state/data/q_BI"),
     read(logs, "/smallsat/body/state/data/ω_BI_B"))
end

display(plot(t, transpose(q_BI) - transpose(q_BI_baseline),
    label  = ["q1" "q2" "q3" "q4"],
    xlabel = "Time (s)",
    ylabel = "Quaternion Differences",
    title  = "Differences in Attitude Quaternions"))
display(plot(t, 180/π * (transpose(ω_BI_B) - transpose(ω_BI_B_baseline)),
    label  = ["ω1" "ω2" "ω3"],
    xlabel = "Time (s)",
    ylabel = "Rate (deg/s)",
    title  = "Differences in Rotation Rate"))

end # module


text:
# TODO: Try to make as many types as possible immutable
# TODO: Generalize to arbitrary metric space

abstract AbstractBall

type EmptyBall <: AbstractBall
end

type Ball{T <: Real} <: AbstractBall
	center::Vector{T}
	r::Float64
end

function Base.contains(b1::Ball, b2::Ball)
	d = euclidean(b1.center, b2.center)
	return d + b2.r <= b1.r
end

function intersects(b1::Ball, b2::Ball)
	d = euclidean(b1.center, b2.center)
	return d <= b1.r + b2.r
end


text:
using Chain
using DataFrames
using Dates

using Stonks: AbstractStonksRecord, Symbols, SchemaValidationError, UpdatableSymbol
using Stonks:
  construct_updatable_symbols, create_typed_dataframe, last_workday, get_data, to_dataframe
using Stonks.APIClients: APIClient, APIResource, get_resource
using Stonks.Stores: FileStore

function validate_schema(::Type{T}, data::DataFrame) where {T<:AbstractStonksRecord}
  get_type_param(::Vector{S}) where {S} = S
  ds_types = [
    (name=String(name), type=type) for (name, type) in zip(fieldnames(T), T.types)
  ]
  df_types = [(name=col, type=get_type_param(data[:, col])) for col in names(data)]
  if length(df_types) != length(ds_types)
    return SchemaValidationError("Length of column names of $T != data")
  end
  for i in 1:length(ds_types)
    if ds_types[i].name != df_types[i].name || ds_types[i].type !== df_types[i].type
      return SchemaValidationError(
        "$(ds_types[i].name), $(ds_types[i].type) != $(df_types[i].name), $(df_types[i].type)",
      )
    end
  end
  return true
end

function apply_schema(::Type{T}, data::DataFrame) where {T<:AbstractStonksRecord}
  ds_types = [
    (name=String(name), type=type) for (name, type) in zip(fieldnames(T), T.types)
  ]
  if length(ds_types) != length(names(data))
    return SchemaValidationError("Length of column names of $T != data")
  end
  for (name, CT) in ds_types
    data[!, name] = convert.(CT, data[:, name])
  end
  empty_df = create_typed_dataframe(T)
  return vcat(empty_df, data)
end

function init(ds::FileStore{T}) where {T<:AbstractStonksRecord}
  df = create_typed_dataframe(T)
  if isempty(ds.partitions)
    tr = WriteTransaction([WriteOperation(df, ds.format)], ds.format, ds.path)
    return execute(tr, ds.writer)
  end
  tr = @chain ds.partitions begin
    map(p -> "$p=__init__", _)
    reduce(joinpath, _)
    [WriteOperation(df, ds.format, _)]
    WriteTransaction(_, ds.format, ds.path)
  end
  return execute(tr, ds.writer)
end

function filter_files(
  ds::FileStore{T}, partitions::Dict{String,Vector{String}}
) where {T<:AbstractStonksRecord}
  for p in keys(partitions)
    !(String(p) in map(String, ds.partitions)) &&
      return ErrorException("$p is not a partition key")
  end
  all_files = list_nested_files(ds.path)
  files_to_read = String[]
  for file in all_files
    regex_matches = @chain partitions begin
      [(rmatch=match(Regex("(?<=$p=)[^/]+(?=/)"), file), values=vals) for (p, vals) in _]
      filter(x -> !isnothing(x.rmatch), _)
    end
    if !isempty(regex_matches)
      flags = map(x -> x.rmatch.match in x.values, regex_matches)
      all(flags) && push!(files_to_read, file)
    end
  end
  #@info "After partition filters, got $(length(files_to_read)) / $(length(all_files))."
  return files_to_read
end

# TODO can be done recursively
function list_nested_files(dir::String)::Vector{String}
  nested_files = []
  for (root, dirs, files) in walkdir(dir)
    foreach(file -> push!(nested_files, joinpath(root, file)), files)
  end
  return nested_files
end

function list_partition_nesting(dir::String)::Vector{String}
  partitions = []
  for (root, dirs, files) in walkdir(dir)
    isempty(dirs) && continue
    d = first(dirs)
    !(isdir(joinpath(root, d)) && contains(d, "=")) && continue
    p = first(split(d, "="))
    !in(p, partitions) ? push!(partitions, p) : continue
  end
  return partitions
end

function generate_partition_pairs(df::DataFrame, cols::Vector{AbstractString})
  df_cols = names(df)
  diff = setdiff(cols, df_cols)
  if length(diff) > 0
    cols_excl = join(filter(x -> !in(x, df_cols), cols), ", ")
    return ErrorException(
      "There are partition columns not present in the dataframe: '$cols_excl'"
    )
  end
  gdf = groupby(df, cols)
  partitions = [
    (key=join(map(x -> "$(x[1])=$(x[2])", zip(names(key), values(key))), "/"), data=subdf)
    for (key, subdf) in pairs(gdf)
  ]
  return partitions
end

function symbols_from_df(
  df::AbstractDataFrame,
  ids::AbstractVector{T};
  time_column::AbstractString,
  ref_date::Date=last_workday(),
)::Union{Symbols,Nothing} where {T<:AbstractString}
  groups = @chain df begin
    groupby(_, map(Symbol, ids))
    combine(_, Symbol(time_column) => maximum)
    transform(
      _, Symbol("$(time_column)_maximum") => ByRow(d -> d + Dates.Day(1)) => :lower_limit
    )
    filter(:lower_limit => <(ref_date), _)
  end
  if nrow(groups) == 0
    return nothing
  end
  id_lengths = map(id ->  map(length, df[:, Symbol(id)]), ids) |> maximum
  if length(ids) > 1
    if id_lengths[1] == 3 && id_lengths[2] == 3
      return [
        ("$(row[ids[1]])/$(row[ids[2]])", row[:lower_limit], ref_date) for row in eachrow(groups)
      ]
    else
      return nothing
    end
  end
  return [(row[Symbol(first(ids))], row[:lower_limit], ref_date) for row in eachrow(groups)]
end

function find_update_candidates(
  df::AbstractDataFrame,
  symbols::Union{Symbols,Nothing};
  ids::AbstractVector{T},
  time_column::Union{AbstractString,Missing},
  ref_date=last_workday(),
)::Vector{UpdatableSymbol} where {T<:AbstractString}
  if isnothing(symbols) && nrow(df) == 0
    return []
  elseif isnothing(symbols) && !ismissing(time_column)
    smb = symbols_from_df(df, ids; time_column=time_column, ref_date=ref_date)
    return !isnothing(smb) ? construct_updatable_symbols(smb) : []
  elseif isnothing(symbols) && ismissing(time_column) && length(ids) == 1
    df_vals = unique(df[:, Symbol(first(ids))])
    return construct_updatable_symbols(df_vals)
  elseif isnothing(symbols) && ismissing(time_column) && lenght(ids) == 2
    tickers = @chain begin
      unique(df[:, map(Symbol, ids)])
      ["$(row[1])$(row[2])=X" for row in eachrow(_)]
      construct_updatable_symbols(_)
    end
    return tickers
  else
    ()
  end
  tickers = construct_updatable_symbols(symbols)
  if all(map(t -> isa(t.from, Date) || isa(t.to, Date), tickers))
    return tickers
  end
  df_vals = unique(df[:, map(Symbol, ids)])
  if all(map(t -> !ismissing(t.fx_pair), tickers)) && length(ids) == 2
    g1 = map(t -> first(t.fx_pair), tickers)
    g2 = map(t -> last(t.fx_pair), tickers)
    fx_pairs = map(t -> (t.fx_pair[1], t.fx_pair[2]), tickers)
    g_diff = [
      "$(p[1])/$(p[2])" for p in fx_pairs if
      !(p[1] in df_vals[:, Symbol(ids[1])] && p[2] in df_vals[:, Symbol(ids[2])])
    ]
    data = filter(row -> row[Symbol(ids[1])] in g1 && row[Symbol(ids[2])] in g2, df)
    smb_df = (
      if !isempty(data) && !ismissing(time_column)
        symbols_from_df(data, ids; time_column=time_column, ref_date=ref_date)
      elseif ismissing(time_column) && !isempty(fx_pairs)
        ["$(p[1])/$(p[2])" for p in fx_pairs]
      else
        nothing
      end
    )
    tickers_df = !isnothing(smb_df) ? construct_updatable_symbols(smb_df) : nothing
    tickers_diff = !isempty(g_diff) ? construct_updatable_symbols(g_diff) : nothing
    return vcat(
      !isnothing(tickers_df) ? tickers_df : [], !isnothing(tickers_diff) ? tickers_diff : []
    )
  else
    g1 = map(x -> x.ticker, tickers)
    g_diff = [
      t for t in map(t -> t.ticker, tickers) if !(t in df_vals[:, Symbol(first(ids))])
    ]
    data = filter(row -> row[Symbol(first(ids))] in g1, df)
    smb_df = (
      if !isempty(data) && !ismissing(time_column)
        symbols_from_df(data, ids; time_column=time_column, ref_date=ref_date)
      elseif ismissing(time_column) && !isempty(tickers)
        map(t -> t.ticker, tickers)
      else
        nothing
      end
    )
    tickers_df = !isnothing(smb_df) ? construct_updatable_symbols(smb_df) : nothing
    tickers_diff = !isempty(g_diff) ? construct_updatable_symbols(g_diff) : nothing
    return vcat(
      !isnothing(tickers_df) ? tickers_df : [], !isnothing(tickers_diff) ? tickers_diff : []
    )
  end
end


text:
using DashBootstrapComponents, DashHtmlComponents

card = dbc_card(
    [
        dbc_cardimg(src = "/static/images/placeholder286x180.png", top = true),
        dbc_cardimgoverlay(
            dbc_cardbody([
                html_h4("Card title", className = "card-title"),
                html_p(
                    "An example of using an image in the background of " * "a card.",
                    className = "card-text",
                ),
                dbc_button("Go somewhere", color = "primary"),
            ],),
        ),
    ],
    style = Dict("width" => "18rem"),
)

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
#=Робот - в произвольной клетке ограниченного прямоугольного поля, на котором могут находиться также внутренние прямоугольные перегородки
(все перегородки изолированы друг от друга, прямоугольники могут вырождаться в отрезки)

РЕЗУЛЬТАТ: Робот - в исходном положении и в углах поля стоят маркеры=#
include("roblib.jl")

function mark_angles(r)
    num_steps=[] #создаём массив для num_steps
    while isborder(r,Sud)==false || isborder(r,West)==false #по условию перегородки не могут соприкасаться->соприкосновение-в углах
        push!(num_steps, get_num_movements!(r, West))
        push!(num_steps, get_num_movements!(r, Sud))
    end
#=функция push! определена по умолчанию-она увеличивает кол-во элементов массива
то есть,если элементов было 4,то после get_num_movements! их стало 5(5 элементов где get_num_movements!-даёт 5-ый элемент) =#

    for side in (Nord,Ost,Sud,West)
        movements!(r,side)
        putmarker!(r)
    end

    for (i,n) in enumerate(num_steps)
        side = isodd(i) ? Ost : Nord
        movements!(r,side,n)
    end
#isodd(i::Int)-проверка числа на нечётность(true-если i-нечётный)
#iseven(i::Int)-аналогично,но для чётных
#функция enumerate-генерирует кортеж из 2 элементов(индекс и сам элемент)
end
function moves!(r,side)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end
#уже стандартная функция для подсчёта шагов

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
#!/usr/bin/env julia

using Distributed
if Distributed.nworkers() > 1
    Distributed.rmprocs(workers())
end

max_node_processes = 12

"turns a slurm node list into a list of node names"
function parse_node_names(slurm_node_list::String)
    node_names = String[]
    slurm_node_list_parts = split(slurm_node_list, "[")
    node_name_prefix = slurm_node_list_parts[1]

    #if only one node was found then this length will be 1
    if length(slurm_node_list_parts) == 1
        push!(node_names, node_name_prefix)
    else
        # this assumes all nodes have a common prefix and the nodelist has the form [1-2,4,5-10,12]
        node_ids = split(split(slurm_node_list_parts[2], "]")[1], ",")
        for node_id in node_ids
            node_id_parts = split(node_id, "-")
            start = parse(Int, node_id_parts[1])
            stop = parse(Int, length(node_id_parts) == 1 ? node_id_parts[1] : node_id_parts[2])
            digits = length(node_id_parts[1])
            for i in start:stop
                node_name = "$(node_name_prefix)$(lpad(string(i), digits, "0"))"
                push!(node_names, node_name)
            end
        end
    end
    return node_names
end

function add_local_procs()
    node_processes = min(trunc(Int, Sys.CPU_THREADS*0.75), max_node_processes)
    println("local processes: $(node_processes) of $(Sys.CPU_THREADS)")

    proc_ids = Distributed.addprocs(node_processes)

    return proc_ids
end

function add_remote_procs()
    if haskey(ENV, "SLURM_JOB_NODELIST")
        node_list = ENV["SLURM_JOB_NODELIST"]
    elseif haskey(ENV, "SLURM_NODELIST")
        node_list = ENV["SLURM_NODELIST"]
    else
        @info("unable to find slurm node list environment variable")
        return Int[]
    end

    println(node_list)

    node_names = parse_node_names(node_list)
    println(node_names)

    node_names = [name for name in node_names if name != gethostname()]

    node_processes = min(trunc(Int, Sys.CPU_THREADS*0.75), max_node_processes)
    println("remote processes per node: $(node_processes)/$(Sys.CPU_THREADS)")
    proc_ids = Int[]

    for i in 1:node_processes
        node_proc_ids = Distributed.addprocs(node_names, sshflags="-oStrictHostKeyChecking=no")
        println("process id batch $(i) of $(node_processes): $(node_proc_ids)")
        for npid in node_proc_ids
            push!(proc_ids, npid)
        end
    end
    println("process ids: $(proc_ids)")

    return proc_ids
end

local_proc_ids = add_local_procs()
println("local processes ids: $(local_proc_ids)")

@everywhere function test_it(pid)
    println("pid: $(pid) - $(gethostname()) - $(pwd())")
    return "$(gethostname())"
end

processes = Distributed.nprocs()-1 # save one for the master process
println("nprocs: $(processes)")
result = pmap(test_it, 1:processes)
println(result)

remote_proc_ids = add_remote_procs()
println("remote processes ids: $(remote_proc_ids)")

@everywhere function test_it(pid)
    println("pid: $(pid) - $(gethostname()) - $(pwd())")
    return "$(gethostname())"
end

processes = Distributed.nprocs()-1 # save one for the master process
println("nprocs: $(processes)")
result = pmap(test_it, 1:processes)
println(result)

max_stars_repo_name: ccoffrin/julia-examples

max_issues_repo_name: ccoffrin/julia-examples

max_forks_repo_name: ccoffrin/julia-examples


text:
"""
    forward(hmm, obs, approx, k=3)

Compute forward probabilities for `hmm` with observations `obs` using likelihood approximation `approx` of order `k`. Returns an [`HMMApproximator`](@ref) object.
"""
function forward(hmm::T, obs::Vector, approx::Type{L}, k::Int = 3) where {T<:HMM,L<:LikelihoodApproximation}
  @unpack logP, logπ₀ = hmm

  n = length(obs)
  nk = nstates(hmm)
  r = k - 1

  # initialize likelihood approximation
  likelihood = approx(hmm)

  # iterator over subsequences of length (k-1)
  subseqs = Iterators.product([1:nk for _ in 1:r]...)

  # normalization constants and likelihoods allocation
  zi = [Vector{Float64}(undef, nk^r) for _ in 1:n-r]
  ll = [Vector{Float64}(undef, nk^k) for _ in 1:n-r]

  q   = Vector{Float64}(undef, nk)
  seq = Vector{Int}(undef, k)

  # top sequence
  @inbounds for s in subseqs
    for s1 = 1:nk
      seq = vcat(s1, s...)
      q[s1] = logπ₀[s1] + likelihood(seq[1:1], obs, 1:1) / k
      for i = 2:k
        q[s1] += logP[seq[i-1], seq[i]] + likelihood(seq[1:i], obs, 1:i) / k
      end
      ll[1][encode(seq, nk)] = q[s1]
    end
    zi[1][encode(s, nk)] = logsumexp(q)
  end

  # intermediate subsequences
  @inbounds for t = 2:(n-k)
    range = t:t+r
    for s in subseqs
      for s1 = 1:nk
        seq = vcat(s1, s...)
        q[s1] = likelihood(seq, obs, range) / k + zi[t-1][encode(seq[1:end-1], nk)]
        ll[t][encode(seq, nk)] = q[s1]
      end

      zi[t][encode(s, nk)] = logsumexp(q) + logP[s[end-1], s[end]]
    end
  end

  # last element
  @inbounds for s in subseqs
    for s1 = 1:nk
      seq = vcat(s1, s...)
      q[s1] = likelihood(seq[end:end], obs, n:n) / k

      for i = 1:r
        range = n-k+i:n
        q[s1] += likelihood(seq[i:end], obs, range) / k
      end

      q[s1] += zi[end-1][encode(seq[1:end-1], nk)]
      ll[end][encode(seq, nk)] = q[s1]
    end

    zi[end][encode(s, nk)] = logsumexp(q) + logP[s[end-1], s[end]]
  end

  # normalization constant
  z = zi[end] .- logsumexp(zi[end])

  HMMApproximator(hmm, obs, zi, ll, z, n, nk, k)
end


text:
function get_RealQuantity(xmle::LightXML.XMLElement, str::String)
  dx = 0.0
  x = 0.0
  for z in child_elements(xmle)
    if name(z) == str
      for y in child_elements(z)
        if name(y) == "value"
          x = parse(Float64, content(y))
        elseif name(y) == "uncertainty"
          dx = parse(Float64, content(y))
        end
      end
    end
  end
  return (x, dx)
end

function parse_qml(evt::XMLElement)
  MT = Array{SeisSrc,1}(undef, 0)
  evt_src = attribute(evt, "publicID")
  evt_id = split(evt_src, r"[;:/=]")[end]
  preferredOriginID = ""
  preferredMagnitudeID = ""
  preferredFocalMechanismID = ""

  mags = Array{XMLElement,1}(undef,0)
  origs = Array{XMLElement,1}(undef,0)
  mechs = Array{XMLElement,1}(undef,0)
  eqtype = ""

  for ch in child_elements(evt)
    cn = name(ch)
    if cn == "preferredOriginID"
      preferredOriginID = split(content(ch), r"[;:/=]")[end]
    elseif cn == "preferredMagnitudeID"
      preferredMagnitudeID = split(content(ch), r"[;:/=]")[end]
    elseif cn == "preferredFocalMechanismID"
      preferredFocalMechanismID = split(content(ch), r"[;:/=]")[end]
    elseif cn == "type"
      eqtype = content(ch)
    elseif cn == "magnitude"
      push!(mags, ch)
    elseif cn == "origin"
      push!(origs, ch)
    elseif cn == "focalMechanism"
      push!(mechs, ch)
    end
  end

  # ==========================================================================
  # Focal Mechanism
  n = 0
  while n < length(mechs)
    n         = n+1
    auth      = ""
    mech      = getindex(mechs, n)
    np        = zeros(Float64, 3, 2)
    pax       = zeros(Float64, 3, 3)
    mt        = zeros(Float64, 6)
    dm        = fill!(zeros(Float64, 6), Inf)

    # generate SeisSrc object
    S = SeisSrc()
    ssrc = attribute(mech, "publicID")
    setfield!(S, :id, String(split(ssrc, r"[;:/=]")[end]))

    for child in child_elements(mech)
      cname = name(child)

      if cname == "creationInfo"
        for grandchild in child_elements(child)
          if name(grandchild) ==  "author"
            auth = content(grandchild)
          end
        end

      elseif cname == "azimuthalGap"
        setfield!(S, :gap, parse(Float64, content(child)))

      elseif cname == "stationPolarityCount"
        setfield!(S, :npol, parse(Int64, content(child)))

      elseif cname == "nodalPlanes"
        j = 0
        for np2 in child_elements(child)  # nodalPlane1, nodalPlane2
          if name(np2) != "preferredPlane"
            for grandchild in child_elements(np2) # strike, dip, rake
              for greatgrandchild in child_elements(grandchild)
                if name(greatgrandchild) == "value"
                  j = j+1
                  setindex!(np, parse(Float64, content(greatgrandchild)), j)
                end
              end
            end
          end
        end
        S.misc["planes_desc"] = "strike, dip, rake"

      elseif cname == "principalAxes"
        j = 0
        for grandchild in child_elements(child) # tAxis, pAxis, nAxis
          for greatgrandchild in child_elements(grandchild) # azimuth, plunge, length
            for redheadedstepchild in child_elements(greatgrandchild)
              if name(redheadedstepchild) == "value"
                j = j+1
                setindex!(pax, parse(Float64, content(greatgrandchild)), j)
              end
            end
          end
        end
        S.misc["pax_desc"] = "azimuth, plunge, length"

      elseif cname == "momentTensor"
        S.misc["mt_id"] = attribute(child, "publicID")
        for a in child_elements(child)
          aname = name(a)
          if aname == "scalarMoment"
            setfield!(S, :m0, parse(Float64, content(a)))
          elseif aname == "tensor"
            j = 0
            for k in ("Mrr", "Mtt", "Mpp", "Mrt", "Mrp", "Mtp")
              j = j+1
              (mt[j], dm[j]) = get_RealQuantity(a, k)
            end
          elseif aname == "sourceTimeFunction"
            for b in child_elements(a)
              if name(b) == "type"
                S.st.desc = content(b)
              elseif name(b) == "duration"
                S.st.dur = parse(Float64, content(b))
              elseif name(b) == "riseTime"
                S.st.rise = parse(Float64, content(b))
              elseif name(b) == "decayTime"
                S.st.decay = parse(Float64, content(b))
              end
            end
          elseif aname == "derivedOriginID"
            S.misc["derivedOriginID"] = String(split(content(a), r"[;:/=]")[end])
            S.misc["xmt_derivedOriginID"] = string(a)
          else
            S.misc["xmt_" * aname] = string(a)
          end
        end
      elseif cname == "methodID"
        S.misc["methodID"] = content(child)
      else
        S.misc["xmech_" * cname] = string(child)
      end

    end
    src_str = ssrc * "," * auth
    setfield!(S, :planes, np)
    setfield!(S, :pax, pax)
    setfield!(S, :mt, mt)
    setfield!(S, :dm, dm)
    setfield!(S, :src, src_str)
    note!(S, "+origin ¦ " * src_str)
    push!(MT, S)
  end

  # ==========================================================================
  # Choose a focal mechanism to retain; then determine mag_id and loc_id
  mag_id = ""
  loc_id = ""
  if length(MT) > 0
    # Store moment tensor with lowest M_err
    if isempty(preferredFocalMechanismID)
      sort!(MT, by=x->sum(abs.(x.dm)))
      R = getindex(MT, 1)
    else
      R = getindex(MT, findfirst([occursin(getfield(m, :id), preferredFocalMechanismID) for m in MT]))
    end

    # Set ID
    if isempty(preferredMagnitudeID)
      mag_id = get(getfield(R, :misc), "derivedOriginID", "")
    else
      mag_id = preferredMagnitudeID
    end
    loc_id = isempty(preferredOriginID) ? mag_id : preferredOriginID
  else
    R = SeisSrc()
  end

  # ==========================================================================
  # Magnitude
  H = SeisHdr()
  setfield!(R, :eid, identity(getfield(H, :id)))

  n = 0
  oid = ""
  pid = ""
  auth = ""
  m = -5.0f0
  msc = ""
  nst = 0
  gap = 0.0
  while n < length(mags)
    n = n+1
    mag = getindex(mags, n)
    oid = ""
    pid = ""
    auth = ""
    m = -5.0f0
    msc = ""
    nst = 0
    gap = 0.0

    # src
    originID = get_elements_by_tagname(mag, "originID")
    pid = attribute(mag, "publicID")
    oid = isempty(originID) ? "" : content(first(originID))

    # msc
    mtype = get_elements_by_tagname(mag, "type")
    msc = isempty(mtype) ? "" : content(first(mtype))

    if occursin(oid, mag_id) || n == 1 || startswith(lowercase(msc), "mw")
      for child in child_elements(mag)
        cname = name(child)
        if cname == "mag"
          m = parse(Float32, content(first(get_elements_by_tagname(child, "value"))))
        elseif cname == "stationCount"
          nst = parse(Int64, content(child))
        elseif cname == "azimuthalGap"
          gap = parse(Float64, content(child))
        elseif cname == "creationInfo"
          for grandchild in child_elements(child)
            if name(grandchild) ==  "author"
              auth = content(grandchild)
            end
          end
        end
      end
    end
  end
  MAG = EQMag(m, msc, nst, gap, pid * "," * oid * "," * auth)
  if MAG.src == ",,"
    MAG.src = ""
  end

  # ==========================================================================
  # Location
  n = 0
  nst = 0
  ot = "1970-01-01T00:00:00"
  gap = 0.0
  auth = ""
  ltyp = ""
  loc_src = ""
  locflags = Array{Char, 1}(undef,8)
  fill!(locflags, '0')
  loc = zeros(Float64, 12)
  while n < length(origs)
    n = n+1
    orig = getindex(origs, n)

    # _______________________________________________________
    # Only parse locations corresponding to a desirable ID
    if occursin(attribute(orig, "publicID"), loc_id) || (n == 1)
      loc_src = attribute(orig, "publicID")

      # Reset temp variables
      nst = 0
      gap = 0.0
      auth = ""
      ltyp = ""
      ot = "1970-01-01T00:00:00"
      fill!(locflags, '0')
      fill!(loc, zero(Float64))

      # Try to set location first
      j = 0
      for str in ("latitude", "longitude", "depth")
        j = j + 1
        (loc[j], dloc) = get_RealQuantity(orig, str)
        if j == 1
          loc[5] = dloc
        elseif j == 2
          loc[4] = dloc
        elseif j == 3
          setindex!(loc, getindex(loc, 3)*0.001, 3)
          loc[6] = dloc
        end
      end

      # Now loop
      for child in child_elements(orig)
        cname = name(child)

        #= originUncertainty removed 2019-11-01
          No documentation of originUncertainty is known to exist;
          can't ascertain intended meaning/use of originUncertainty;
          can't tell if observatories use it uniformly
        =#

        if cname == "type"
          ltyp_tmp = content(child)
          if ltyp_tmp != ""
            ltyp = ltyp_tmp
          end

        elseif cname == "time"
          for grandchild in child_elements(child)
            gcname = name(grandchild)
            if gcname == "value"
              ot = content(grandchild)
            elseif gcname == "uncertainty"
              loc[7] = parse(Float64, content(grandchild))
            end
          end

        elseif cname == "quality"
          for grandchild in child_elements(child)
            gcname = name(grandchild)
            if gcname == "standardError"
              loc[8] = parse(Float64, content(grandchild))
            elseif gcname == "azimuthalGap"
              loc[10] = parse(Float64, content(grandchild))
            elseif gcname == "minimumDistance"
              loc[11] = parse(Float64, content(grandchild))
            elseif gcname == "maximumDistance"
              loc[12] = parse(Float64, content(grandchild))
            elseif gcname == "associatedStationCount"
              nst = parse(Int64, content(grandchild))
            end
          end

        elseif cname == "epicenterFixed"
          locflags[1] = content(child)[1]
          locflags[2] = locflags[1]

        elseif cname == "depthType"
          if content(child) == "operator assigned"
            locflags[3] = '1'
          end

        elseif cname == "timeFixed"
          locflags[4] = content(child)[1]

        elseif cname == "creationInfo"
          for grandchild in child_elements(child)
            if name(grandchild) ==  "author"
              auth = content(grandchild)
            end
          end
        elseif (cname in ("latitude", "longitude", "depth")) == false
          H.misc["xloc_" * cname] = string(child)
        end
      end
      # _______________________________________________________

    end
  end
  LOC = EQLoc(loc..., nst, parse(UInt8, join(locflags), base=2), "", ltyp, "", loc_src * "," * auth )

  setfield!(H, :id, String(evt_id))
  setfield!(H, :loc, LOC)
  setfield!(H, :ot, DateTime(replace(ot, r"[A-S,U-Z,a-z]" => "")[1:min(end,23)]))
  setfield!(H, :mag, MAG)
  setfield!(H, :typ, eqtype)
  setfield!(H, :src, evt_src)
  setfield!(R, :eid, String(evt_id))
  note!(H, "+origin ¦ " * evt_src)

  # Post-process: ensure these are empty if not initialized
  if R.mt == zeros(Float64,6)
    R.mt = Float64[]
  end

  if R.dm == Inf64.*ones(Float64,6)
    R.dm = Float64[]
  end

  if R.pax == zeros(Float64, 3, 3)
    R.pax = Array{Float64, 2}(undef, 0, 0)
  end

  if R.planes == zeros(Float64, 3, 2)
    R.planes = Array{Float64, 2}(undef, 0, 0)
  end

  if R.src == ","
    R.src = ""
  end
  return H, R
end

function event_xml!(EvCat::Array{SeisHdr,1}, EvSrc::Array{SeisSrc, 1}, xdoc::XMLDocument)
  qxml  = first(child_elements(root(xdoc)))
  elts  = child_elements(qxml)
  for elt in elts
    if name(elt) == "event"
      H, R = parse_qml(elt)
      push!(EvCat, H)
      push!(EvSrc, R)

    end
  end
  return nothing
end

"""
    EvCat, EvSrc = read_qml(fpat::String)

Read QuakeML files matching string pattern `fpat`. Returns an array of `SeisHdr`
objects as `EvCat` and an array of `SeisSrc` objects as `EvSrc`, such that
`EvCat[i]` and `EvSrc[i]` describe the preferred location (origin) and
preferred event source (focal mechanism or moment tensor) of event `i`.

"""
function read_qml(fpat::String)
  files = safe_isfile(fpat) ? [fpat] : ls(fpat)
  EvCat = Array{SeisHdr,1}()
  EvSrc = Array{SeisSrc,1}()
  for file in files
    xdoc  = parse_file(file)
    k = length(EvCat)
    event_xml!(EvCat, EvSrc, xdoc)
    L = length(EvCat)
    fstr = abspath(file)
    for i in k+1:L
      setfield!(EvCat[i], :src, fstr)
      note!(EvCat[i], "+source ¦ " * fstr)
      setfield!(EvSrc[i], :src, fstr)
      note!(EvSrc[i], "+source ¦ " * fstr)
    end
    free(xdoc)
  end
  return EvCat, EvSrc
end

function is_qml(k::String, s::String)
  (startswith(s, "<") && endswith(s, ">")) || return false
  i1 = first(findnext("<", s,1))
  j1 = last(findnext(">", s,1))
  s1 = s[nextind(s, i1):prevind(s, j1)]
  (split(k, "_")[2] == s1) || return false
  i2 = first(findlast("<", s))
  j2 = last(findlast(">", s))
  s2 = s[nextind(s, i2):prevind(s, j2)]
  ("/" * s1) == s2 || return false
  return true
end

function new_qml!(io::IO)
  write(io, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>
  <quakeml xmlns=\"http://quakeml.org/xmlns/quakeml/1.2\">\n  <eventParameters publicID=\"smi:SeisIO.jl\">\n    <creationInfo>\n      <agencyID>SeisIO</agencyID>\n      <creationTime>")
  print(io, now())
  write(io, "</creationTime>\n    </creationInfo>\n")
  return nothing
end

function write_misc(io::IO, D::Dict{String,Any}, pref::String, p::Int64)
  for k in keys(D)
    if startswith(k, pref)
      xv = D[k]
      if is_qml(k, xv)
        write(io, " "^2p)
        write(io, xv)
        write(io, "\n")
      end
    end
  end
  return nothing
end

function write_real(io::IO, str::String, x::Union{DateTime,AbstractFloat}, p::Int64)
  write(io, " "^2p, "<")
  write(io, str)
  write(io, ">\n", " "^(2p+2), "<value>")
  print(io, x)
  write(io, "</value>\n", " "^2p, "</")
  write(io, str)
  write(io, ">\n")
  return nothing
end

function write_real(io::IO, str::String, x::Union{DateTime,AbstractFloat}, dx::AbstractFloat, p::Int64)
  write(io, " "^2p, "<")
  write(io, str)
  write(io, ">\n", " "^(2p+2), "<value>")
  print(io, x)
  write(io, "</value>\n", " "^(2p+2), "<uncertainty>")
  print(io, dx)
  write(io, "</uncertainty>\n", " "^2p, "</")
  write(io, str)
  write(io, ">\n")
  return nothing
end

function write_pax(io::IO, pax::Array{Float64,2})
  str = ("t", "p", "n")
  nP = size(pax, 2)
  write(io, "        <principalAxes>\n")
  for i in 1:nP
    write(io, "          <", str[i], "Axis>\n")
    write(io, "            <azimuth>\n              <value>")
    print(io, pax[1,i])
    write(io, "</value>\n            </azimuth>\n")
    write(io, "            <plunge>\n              <value>")
    print(io, pax[2,i])
    write(io, "</value>\n            </plunge>\n")
    write(io, "            <length>\n              <value>")
    print(io, pax[3,i])
    write(io, "</value>\n            </length>\n")
    write(io, "          </", str[i], "Axis>\n")
  end
  write(io, "        </principalAxes>\n")
  return nothing
end

function write_qml!(io::IO, HDR::Array{SeisHdr,1}, SRC::Array{SeisSrc,1}, v::Integer)
  Ri = zeros(Int64, length(HDR))
  R_id = Array{String,1}(undef, length(SRC))
  for i in 1:length(SRC)
    R_id[i] = getfield(SRC[i], :eid)
  end

  for i in 1:length(HDR)
    id = getfield(HDR[i], :id)
    for j in 1:length(R_id)
      if id == R_id[j]
        Ri[i] = j
        break
      end
    end
  end

  for i in 1:length(HDR)
    H = HDR[i]
    (v > 0) && println("Writing event ", H.id)

    i_id = findfirst([occursin("origin", n) for n in H.notes])
    id_str = String(split(i_id == nothing ? " ¦  ¦ " : H.notes[i_id], " ¦ ")[end])
    write(io, "    <event publicID=\"")
    print(io, id_str)
    write(io, "\">\n")
    if isempty(H.loc)
      (v > 0) && println("  Skipped location (H.loc empty)")
    else
      loc_orig, loc_auth, xx = split_id(H.loc.src, c=",")
      write(io, "      <preferredOriginID>")
      write(io, H.id)
      write(io, "</preferredOriginID>\n      <type>")
      write(io, H.typ)
      write(io, "</type>\n")

      # ---------------------------------------------------
      # Origin
      write(io, "      <origin publicID=\"", loc_orig, "\">\n")
      if H.loc.dt > 0.0
        write_real(io, "time", H.ot, H.loc.dt, 4)
      else
        write_real(io, "time", H.ot, 4)
      end
      L = H.loc

      # lat, lon, dep
      if L.dy == 0.0
        write_real(io, "latitude", L.lat, 4)
      else
        write_real(io, "latitude", L.lat, L.dy, 4)
      end
      if L.dx == 0.0
        write_real(io, "longitude", L.lon, 4)
      else
        write_real(io, "longitude", L.lon, L.dx, 4)
      end
      if L.dz == 0.0
        write_real(io, "depth", L.dep*1000.0, 4)
      else
        write_real(io, "depth", L.dep*1000.0, L.dz, 4)
      end

      # flags
      flags = falses(4)
      for n = 1:4
        flags[n] = >>(<<(H.loc.flags, n-1),7)
      end
      if flags[3]
        write(io, "        <depthType>operator assigned</depthType>\n")
      end
      if flags[4]
        write(io, "        <timeFixed>1</timeFixed>\n")
      end
      if flags[1] || flags[2]
        write(io, "        <epicenterFixed>1</epicenterFixed>\n")
      end

      # Location quality
      do_qual = false
      for f in loc_qual_fields
        if getfield(L, f) != 0.0
          do_qual = true
          break
        end
      end
      if do_qual || (L.nst > 0)
        write(io, "        <quality>\n")

        if L.nst > 0
          write(io, "          <associatedStationCount>")
          print(io, L.nst)
          write(io, "</associatedStationCount>\n")
        end

        for (i,f) in enumerate(loc_qual_fields)
          j = getfield(L, f)
          if j != 0.0
            write(io, "          <")
            write(io, loc_qual_names[i])
            write(io, ">")
            print(io, j)
            write(io, "</")
            write(io, loc_qual_names[i])
            write(io, ">\n")
          end
        end
        write(io, "        </quality>\n")
      end

      if !isempty(L.typ)
        write(io, "        <type>")
        print(io, L.typ)
        write(io, "</type>\n")
      end

      # Author
      if !isempty(loc_auth)
        write(io, "        <creationInfo>\n          <author>")
        print(io, loc_auth)
        write(io, "</author>\n        </creationInfo>\n")
      end

      # other location properties
      write_misc(io, H.misc, "xloc_", 4)

      # done Origin
      write(io, "      </origin>\n")
    end

    # ---------------------------------------------------
    # Focal Mechanism
    if Ri[i] == 0
      (v > 0) && println("  Skipped focal mechanism (no R.eid matches)\n  Skipped moment tensor (fields empty)")
    else
      j = Ri[i]
      R = SRC[j]

      i_src = findfirst([occursin("origin", n) for n in R.notes])
      src_str = String(split(i_src == nothing ? " ¦  ¦ " : R.notes[i_src], " ¦ ")[end])
      foc_orig, foc_auth, xx = split_id(src_str, c=",")
      write(io, "      <focalMechanism publicID=\"")
      write(io, foc_orig)
      write(io, "\">\n")

      # Nodal planes
      if !isempty(R.planes)
        write(io, "        <nodalPlanes>\n")
        nP = size(R.planes, 2)
        for i in 1:nP
          write(io, "          <nodalPlane")
          print(io, i)
          write(io, ">\n")
          write(io, "            <strike>\n              <value>")
          print(io, R.planes[1,i])
          write(io, "</value>\n            </strike>\n")
          write(io, "            <dip>\n              <value>")
          print(io, R.planes[2,i])
          write(io, "</value>\n            </dip>\n")
          write(io, "            <rake>\n              <value>")
          print(io, R.planes[3,i])
          write(io, "</value>\n            </rake>\n")
          write(io, "          </nodalPlane")
          print(io, i)
          write(io, ">\n")
        end
        write(io, "        </nodalPlanes>\n")
      end

      # Principal axes
      if !isempty(R.pax)
        write_pax(io, R.pax)
      end

      # Azimuthal gap
      if R.gap != 0.0
        write(io, "        <azimuthalGap>")
        print(io, R.gap)
        write(io, "</azimuthalGap>\n")
      end

      # Polarity count
      if R.npol != 0
        write(io, "        <stationPolarityCount>")
        print(io, R.npol)
        write(io, "</stationPolarityCount>\n")
      end

      # methodID
      if haskey(R.misc, "methodID")
        write(io, "        <methodID>")
        write(io, get(R.misc, "methodID", ""))
        write(io, "</methodID>\n")
      end

      # Author
      if !isempty(foc_auth)
        write(io, "        <creationInfo>\n            <author>")
        write(io, foc_auth)
        write(io, "</author>\n          </creationInfo>\n")
      end

      # Moment Tensor
      if (isempty(R.mt) && (R.m0 == 0.0) && isempty(R.st))
        (v > 0) && println("  Skipped moment tensor (fields empty)")
      else
        write(io, "        <momentTensor publicID=\"")
        mt_id = haskey(R.misc, "mt_id") ? R.misc["mt_id"] : "smi:SeisIO/moment_tensor;fmid=" * R.id
        write(io, mt_id)
        write(io, "\">\n")

        if R.m0 != 0.0
          write_real(io, "scalarMoment", R.m0, 5)
        end

        if isempty(R.mt) == false
          write(io, "          <tensor>\n")
          mt_strings = ("Mrr", "Mtt", "Mpp", "Mrt", "Mrp", "Mtp")
          for i = 1:length(R.mt)
            write_real(io, mt_strings[i], R.mt[i], R.dm[i], 6)
          end
          write(io, "          </tensor>\n")
          write_misc(io, R.misc, "xmt_", 5)
        end

        if !isempty(R.st)
          write(io, "          <sourceTimeFunction>\n")
          write(io, "            <type>")
          write(io, R.st.desc)
          write(io, "</type>\n")
          write(io, "            <duration>")
          print(io, R.st.dur)
          write(io, "</duration>\n")
          write(io, "            <riseTime>")
          print(io, R.st.rise)
          write(io, "</riseTime>\n")
          write(io, "            <decayTime>")
          print(io, R.st.decay)
          write(io, "</decayTime>\n")
          write(io, "          </sourceTimeFunction>\n")
        end
        write(io, "        </momentTensor>\n")
      end

      write(io, "      </focalMechanism>\n")
    end

    # ---------------------------------------------------
    # Magnitude
    if isempty(H.mag)
      (v > 0) && println("  Skipped magnitude (H.mag empty)")
    else
      mag_pid, mag_orig, mag_auth, xx = split_id(H.mag.src, c=",")
      write(io, "      <magnitude publicID=\"")
      write(io, mag_pid)
      write(io, "\">\n")
      write_real(io, "mag", H.mag.val, 4)
      write(io, "        <type>", H.mag.scale, "</type>\n")
      isempty(mag_orig) || write(io, "        <originID>", mag_orig, "</originID>\n")

      if !isempty(mag_auth)
        write(io, "        <creationInfo>\n          <author>", mag_auth, "</author>\n        </creationInfo>\n")
      end

      if H.mag.gap != 0.0
        write(io, "        <azimuthalGap>")
        print(io, H.mag.gap)
        write(io, "</azimuthalGap>\n")
      end

      if H.mag.nst > 0
        write(io, "        <stationCount>")
        print(io, H.mag.nst)
        write(io, "</stationCount>\n")
      end

      write(io, "      </magnitude>\n")
    end
    write(io, "    </event>\n")
  end
  write(io, "</eventParameters>\n</quakeml>\n")
  return nothing
end

@doc """
    write_qml(fname, Ev::SeisEvent; v::Integer=0)

Write event metadata from SeisEvent `Ev` to file `fname`.

    write_qml(fname, SHDR::SeisHdr; v::Integer=0)
    write_qml(fname, SHDR::Array{SeisHdr,1}; v::Integer=0)

Write QML to file `fname` from `SHDR`.

If `fname` exists, and is QuakeML, SeisIO appends the existing XML. If the
file exists, but is NOT QuakeML, an error is thrown; the file isn't overwritten.

    write_qml(fname, SHDR::SeisHdr, SSRC::SeisSrc; v::Integer=0)
    write_qml(fname, SHDR::Array{SeisHdr,1}, SSRC::Array{SeisSrc,1}; v::Integer=0)

Write QML to file `fname` from `SHDR` and `SSRC`.

!!! warning

    To write data from `R ∈ SSRC`, it must be true that `R.eid == H.id` for some `H ∈ SHDR`.
""" write_qml
function write_qml(fname::String, HDR::Array{SeisHdr,1}, SRC::Array{SeisSrc,1}; v::Integer=0)
  H0 = SeisHdr[]
  R0 = SeisSrc[]

  if safe_isfile(fname)
    io = open(fname, "a+")

    # test whether file can be appended
    fastseekend(io)
    fastskip(io, -30)
    test_str = String(read(io))
    if test_str == "</eventParameters>\n</quakeml>\n"

      # behavior for files that are QuakeXML as produced by SeisIO
      fastskip(io, -30)
    else
      try

        # file exists and is readable QuakeXML but not produced by SeisIO
        seekstart(io)
        fstart = String(fastread(io, 5))
        if fstart == "<?xml"
          close(io)
          append!(H0, HDR)
          append!(R0, SRC)
          (H1, R1) = read_qml(fname)
          @assert (isempty(H1) == false)
          @assert (isempty(R1) == false)
          append!(H0, H1)
          append!(R0, R1)
        else
          error("incompatible file type!")
        end
      catch err

        # file exists but isn't QuakeXML
        @warn(string(fname, " isn't valid QuakeXML; can't append, exit with error!"))
        rethrow(err)
      end

      io = open(fname, "w")
      new_qml!(io)
    end
  else

    # new file
    io = open(fname, "w")
    new_qml!(io)
  end

  opts = string(", v=", v)
  if isempty(H0) && isempty(R0)
    write_qml!(io, HDR, SRC, v)
    for h in HDR
      fwrite_note_quake!(h, "asdf_wqml", fname, opts)
    end
    for r in SRC
      fwrite_note_quake!(r, "asdf_wqml", fname, opts)
    end
  else
    write_qml!(io, H0, R0, v)
    for h in H0
      fwrite_note_quake!(h, "asdf_wqml", fname, opts)
    end
    for r in R0
      fwrite_note_quake!(r, "asdf_wqml", fname, opts)
    end
  end
  close(io)
  return nothing
end
write_qml(fname::String, H::SeisHdr, R::SeisSrc; v::Integer=0) = write_qml(fname, [H], [R], v=v)
write_qml(fname::String, HDR::Array{SeisHdr,1}; v::Integer=0) = write_qml(fname, HDR, SeisSrc[], v=v)
write_qml(fname::String, H::SeisHdr; v::Integer=0) = write_qml(fname, [H], SeisSrc[], v=v)
write_qml(fname::String, Ev::SeisEvent; v::Integer=0) = write_qml(fname, [Ev.hdr], [Ev.source], v=v)