text:
"""
    convert_to_time_series(df::DataFrame)

Converts a `DataFrame` loaded via `load_data_file` from the original format, in which each
row contains a month's worth of data for a particular ELEMENT (variable type e.g. precip,
temperature, etc), to a format in which each row contains a single day's worth of data with
for a single ELEMENT type. Only the element types already present in `df` will be included
in the new `DataFrame`.
"""
function convert_to_time_series(df::DataFrame)

    # A dataframe in which each lines containing a single day _and_ a single element.
    tmp = stack(df, 5:size(df, 2))

    # Compute day-of-month.
    tmp.DAY = map(x -> parse(Int, string(x)[6:end]), tmp.variable)

    # Filter to remove invalid dates.
    filter!(row -> Dates.validargs(Date, row.YEAR, row.MONTH, row.DAY) === nothing, tmp)

    # Merge date information into a single column.
    tmp.DATE = Date.(tmp.YEAR, tmp.MONTH, tmp.DAY)

    # Merge element and variable-type into a single underscore-separated column.
    tmp.col_index = map(zip(tmp.ELEMENT, tmp.variable)) do (element, variable)
        element * "_" * string(variable)[1:5]
    end

    # Drop redundant columns.
    tmp2 = tmp[!, [:DATE, :col_index, :value]]

    # Widen the currently very tall data frame.
    tmp3 = unstack(tmp2, [:DATE], :col_index, :value)

    # Reduce type unions over `Any` columns.
    for c in 2:size(tmp3, 2)
        name = names(tmp3)[c]
        if name[6:end] == "VALUE"
            tmp3[!, name] = convert(Vector{Union{Int, Missing}}, tmp3[!, name])
        else
            tmp3[!, name] = convert(Vector{Union{String, Missing}}, tmp3[!, name])
        end
    end

    # Re-order the data in a more natural order.
    sort!(tmp3, :DATE)

    return tmp3
end

"""
    select_data(
        time_interval::Tuple{Date, Date},
        lat_interval::Tuple{Real, Real},
        lon_interval::Tuple{Real, Real},
        element::String,
    )

Load all of the data from a particular region and time period of type `element` e.g. `TMAX`.
See section III of NOAA's readme [1] for more details on available elements.

Returns the results in a ...

Example:
```julia
using Dates, GHCNData
lat_interval = (50, 60);
lon_interval = (-10, 2);
time_interval = (Date(2010), Date(2020));
element = "TMAX"
data = select_data(time_interval, lat_interval, lon_interval, element)
```

[1] - https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
"""
function select_data(
    time_interval::Tuple{Date, Date},
    lat_interval::Tuple{Real, Real},
    lon_interval::Tuple{Real, Real},
    element::String,
)
    # Load all of the inventories.
    inventories = load_inventories()

    # Filter them to exclude stations whose data don't live inside the specified time period
    # and region.
    df = filter(inventories) do row
        row.ELEMENT == element &&
            lat_interval[1] < row.LATITUDE < lat_interval[2] &&
            lon_interval[1] < row.LONGITUDE < lon_interval[2] &&
            (
                time_interval[1] <= Date(row.FIRSTYEAR) <= time_interval[2] ||
                time_interval[1] <= Date(row.LASTYEAR) <= time_interval[2]
            )
    end

    # Load the data from the individual files.
    big_dfs = map(zip(enumerate(df.ID), df.LATITUDE, df.LONGITUDE)) do ((n, id), lat, lon)

        # Load the data and convert it into a time-series format.
        data = load_data_file(id)
        time_series = convert_to_time_series(data)

        # Extract only the columns requested.
        return DataFrame(
            "date" => time_series.DATE,
            "lat_lon_$n" => fill((lat, lon), size(time_series, 1)),
            "value_$n" => getproperty(time_series, Symbol(element * "_VALUE")),
        )
    end

    # Align all dates.
    joint_data_df = foldl(
        (a, b) -> leftjoin(a, b; on=:date, makeunique=true), big_dfs;
        init=DataFrame(date=time_interval[1]:Day(1):(time_interval[2]-Day(1))),
    )

    # Rearrange into a time x lat-lon matrix, and pull out lat-lon pairs.
    value_vecs = map(eachindex(df.ID)) do n
        return getproperty(joint_data_df, Symbol("value_$n"))
    end
    joint_data_matrix = collect(transpose(hcat(value_vecs...)))

    # Collect lat-lon pairs.
    lat_lon_pairs = map(eachindex(df.ID)) do n

        # These should all be missing or equal.
        lat_lon_vals = getproperty(joint_data_df, Symbol("lat_lon_$n"))

        # Pull out unique lat-lon pair.
        if all(ismissing, lat_lon_vals)
            return missing
        else
            return only(unique(filter(!ismissing, lat_lon_vals)))
        end
    end

    present_idx = findall(!ismissing, lat_lon_pairs)
    joint_data_matrix = joint_data_matrix[present_idx, :]
    lat_lon_pairs = lat_lon_pairs[present_idx]
    ID = df.ID[present_idx]

    return joint_data_matrix, lat_lon_pairs, ID
end


text:
# ---------------------------------------------------------------------------- #
#
#   basisFuncTetLag.jl
#
#   Get Lagrangian basis functions for Tetrahedron
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

function basisFuncTetLag( P::P1, s::Vector{Float64}, t::Vector{Float64}, u::Vector{Float64} )

  dim = 3
  nx  = length( s )

  # Linear

  nphi = 4
  phi  = Array{Float64}(nphi, nx)
  dphi = Array{Float64}(nphi, nx, dim)

  # value
  phi[1,:] = 1 - s - t - u # 1 @ node 0 (s = 0, t = 0, u = 0)
  phi[2,:] =     s         # 1 @ node 1 (s = 1, t = 0, u = 0)
  phi[3,:] =         t     # 1 @ node 2 (s = 0, t = 1, u = 0)
  phi[4,:] =             u # 1 @ node 3 (s = 0, t = 0, u = 1)

  # derivative in x
  dphi[1,:,1] = -1
  dphi[2,:,1] =  1
  dphi[3,:,1] =  0
  dphi[4,:,1] =  0

  # derivative in y
  dphi[1,:,2] = -1
  dphi[2,:,2] =  0
  dphi[3,:,2] =  1
  dphi[4,:,2] =  0

  # derivative in z
  dphi[1,:,3] = -1
  dphi[2,:,3] =  0
  dphi[3,:,3] =  0
  dphi[4,:,3] =  1

  return phi, dphi

end

function basisFuncTetLag( P::P2, s::Vector{Float64}, t::Vector{Float64}, u::Vector{Float64} )

  dim = 3
  nx  = length( s )

  # Quadratric
  nphi = 10
  phi  = Array{Float64}(nphi, nx)
  dphi = Array{Float64}(nphi, nx, dim)

  # precompute
  ss = (s.*s)
  tt = (t.*t)
  uu = (u.*u)

  st = s.*t
  su = s.*u
  tu = t.*u

  # value
  phi[1,:] =  s.*-3.0-t.*3.0-u.*3.0+st.*4.0+su.*4.0+tu.*4.0+ss.*2.0+tt.*2.0+uu.*2.0+1.0
  phi[2,:] =  -s+ss.*2.0
  phi[3,:] =  -t+tt.*2.0
  phi[4,:] =  -u+uu.*2.0
  phi[5,:] =  tu.*4.0
  phi[6,:] =  su.*4.0
  phi[7,:] =  st.*4.0
  phi[8,:] =  t.*4.0-st.*4.0-tu.*4.0-tt.*4.0
  phi[9,:] =  u.*4.0-su.*4.0-tu.*4.0-uu.*4.0
  phi[10,:] =  s.*4.0-st.*4.0-su.*4.0-ss.*4.0

  # derivative in x
  dphi[1,:,1] =  s.*4.0+t.*4.0+u.*4.0-3.0
  dphi[2,:,1] =  s.*4.0-1.0
  dphi[3,:,1] =  0.0
  dphi[4,:,1] =  0.0
  dphi[5,:,1] =  0.0
  dphi[6,:,1] =  u.*4.0
  dphi[7,:,1] =  t.*4.0
  dphi[8,:,1] =  t.*-4.0
  dphi[9,:,1] =  u.*-4.0
  dphi[10,:,1] =  s.*-8.0-t.*4.0-u.*4.0+4.0

  # derivative in y
  dphi[1,:,2] =  s.*4.0+t.*4.0+u.*4.0-3.0
  dphi[2,:,2] =  0.0
  dphi[3,:,2] =  t.*4.0-1.0
  dphi[4,:,2] =  0.0
  dphi[5,:,2] =  u.*4.0
  dphi[6,:,2] =  0.0
  dphi[7,:,2] =  s.*4.0
  dphi[8,:,2] =  s.*-4.0-t.*8.0-u.*4.0+4.0
  dphi[9,:,2] =  u.*-4.0
  dphi[10,:,2] =  s.*-4.0

  # derivative in z
  dphi[1,:,3] =  s.*4.0+t.*4.0+u.*4.0-3.0
  dphi[2,:,3] =  0.0
  dphi[3,:,3] =  0.0
  dphi[4,:,3] =  u.*4.0-1.0
  dphi[5,:,3] =  t.*4.0
  dphi[6,:,3] =  s.*4.0
  dphi[7,:,3] =  0.0
  dphi[8,:,3] =  t.*-4.0
  dphi[9,:,3] =  s.*-4.0-t.*4.0-u.*8.0+4.0
  dphi[10,:,3] =  s.*-4.0

  # transpose

  return phi, dphi

end

function basisFuncTetLag( P::P3, s::Vector{Float64}, t::Vector{Float64}, u::Vector{Float64} )

  dim = 3
  nx  = length( s )

  # Cubic
  nphi = 20
  phi  = Array{Float64}(nphi, nx)
  dphi = Array{Float64}(nphi, nx, dim)

  ss = (s.*s)
  tt = (t.*t)
  uu = (u.*u)

  st = s.*t
  su = s.*u
  tu = t.*u

  ssu = ss.*u
  sst = ss.*t
  tuu =  t.*uu
  ttu = tt.*u
  suu =  s.*uu
  stt =  s.*tt

  stu =  st.*u

  ttt = (t.*t.*t)
  uuu = (u.*u.*u)
  sss = (s.*s.*s)

  # value
  phi[1,:] =  s.*(-1.1E1/2.0)-t.*(1.1E1/2.0)-u.*(1.1E1/2.0)+st.*1.8E1+su.*1.8E1+tu.*1.8E1-stt.*(2.7E1/2.0)-
             sst.*(2.7E1/2.0)-suu.*(2.7E1/2.0)-ssu.*(2.7E1/2.0)-tuu.*(2.7E1/2.0)-ttu.*(2.7E1/2.0)+(s.*
             s).*9.0-sss.*(9.0/2.0)+tt.*9.0-ttt.*(9.0/2.0)+uu.*9.0-uuu.*(9.0/2.0)-stu.*2.7E1+1.0
  phi[2,:] =  s-ss.*(9.0/2.0)+sss.*(9.0/2.0)
  phi[3,:] =  t-tt.*(9.0/2.0)+ttt.*(9.0/2.0)
  phi[4,:] =  u-uu.*(9.0/2.0)+uuu.*(9.0/2.0)
  phi[5,:] =  tu.*(-9.0/2.0)+ttu.*(2.7E1/2.0)
  phi[6,:] =  tu.*(-9.0/2.0)+tuu.*(2.7E1/2.0)
  phi[7,:] =  su.*(-9.0/2.0)+suu.*(2.7E1/2.0)
  phi[8,:] =  su.*(-9.0/2.0)+ssu.*(2.7E1/2.0)
  phi[9,:] =  st.*(-9.0/2.0)+sst.*(2.7E1/2.0)
  phi[10,:] =  st.*(-9.0/2.0)+stt.*(2.7E1/2.0)
  phi[11,:] =  t.*(-9.0/2.0)+st.*(9.0/2.0)+tu.*(9.0/2.0)-stt.*(2.7E1/2.0)-ttu.*(2.7E1/2.0)+tt.*1.8E1-
             ttt.*(2.7E1/2.0)
  phi[12,:] =  t.*9.0-st.*(4.5E1/2.0)-tu.*(4.5E1/2.0)+stt.*2.7E1+sst.*(2.7E1/2.0)+tuu.*(2.7E1/2.0)+
             ttu.*2.7E1-tt.*(4.5E1/2.0)+ttt.*(2.7E1/2.0)+stu.*2.7E1
  phi[13,:] =  u.*9.0-su.*(4.5E1/2.0)-tu.*(4.5E1/2.0)+suu.*2.7E1+ssu.*(2.7E1/2.0)+tuu.*2.7E1+tt.*
             u.*(2.7E1/2.0)-uu.*(4.5E1/2.0)+uuu.*(2.7E1/2.0)+stu.*2.7E1
  phi[14,:] =  u.*(-9.0/2.0)+su.*(9.0/2.0)+tu.*(9.0/2.0)-suu.*(2.7E1/2.0)-tuu.*(2.7E1/2.0)+uu.*1.8E1-
             uuu.*(2.7E1/2.0)
  phi[15,:] =  s.*9.0-st.*(4.5E1/2.0)-su.*(4.5E1/2.0)+stt.*(2.7E1/2.0)+sst.*2.7E1+suu.*(2.7E1/2.0)+
             ssu.*2.7E1-ss.*(4.5E1/2.0)+sss.*(2.7E1/2.0)+stu.*2.7E1
  phi[16,:] =  s.*(-9.0/2.0)+st.*(9.0/2.0)+su.*(9.0/2.0)-sst.*(2.7E1/2.0)-ssu.*(2.7E1/2.0)+ss.*1.8E1-
             sss.*(2.7E1/2.0)
  phi[17,:] =  stu.*2.7E1
  phi[18,:] =  tu.*2.7E1-tuu.*2.7E1-ttu.*2.7E1-stu.*2.7E1
  phi[19,:] =  su.*2.7E1-suu.*2.7E1-ssu.*2.7E1-stu.*2.7E1
  phi[20,:] =  st.*2.7E1-stt.*2.7E1-sst.*2.7E1-stu.*2.7E1

  # derivative in x
  dphi[1,:,1] =  s.*1.8E1+t.*1.8E1+u.*1.8E1-st.*2.7E1-su.*2.7E1-tu.*2.7E1-ss.*(2.7E1/2.0)-tt.*(2.7E1/2.0)-
             uu.*(2.7E1/2.0)-1.1E1/2.0
  dphi[2,:,1] =  s.*-9.0+ss.*(2.7E1/2.0)+1.0
  dphi[3,:,1] =  0.0
  dphi[4,:,1] =  0.0
  dphi[5,:,1] =  0.0
  dphi[6,:,1] =  0.0
  dphi[7,:,1] =  u.*(-9.0/2.0)+uu.*(2.7E1/2.0)
  dphi[8,:,1] =  u.*(-9.0/2.0)+su.*2.7E1
  dphi[9,:,1] =  t.*(-9.0/2.0)+st.*2.7E1
  dphi[10,:,1] =  t.*(-9.0/2.0)+tt.*(2.7E1/2.0)
  dphi[11,:,1] =  t.*(9.0/2.0)-tt.*(2.7E1/2.0)
  dphi[12,:,1] =  t.*(-4.5E1/2.0)+st.*2.7E1+tu.*2.7E1+tt.*2.7E1
  dphi[13,:,1] =  u.*(-4.5E1/2.0)+su.*2.7E1+tu.*2.7E1+uu.*2.7E1
  dphi[14,:,1] =  u.*(9.0/2.0)-uu.*(2.7E1/2.0)
  dphi[15,:,1] =  s.*-4.5E1-t.*(4.5E1/2.0)-u.*(4.5E1/2.0)+st.*5.4E1+su.*5.4E1+tu.*2.7E1+ss.*(8.1E1/2.0)+tt.*
             (2.7E1/2.0)+uu.*(2.7E1/2.0)+9.0
  dphi[16,:,1] =  s.*3.6E1+t.*(9.0/2.0)+u.*(9.0/2.0)-st.*2.7E1-su.*2.7E1-ss.*(8.1E1/2.0)-9.0/2.0
  dphi[17,:,1] =  tu.*2.7E1
  dphi[18,:,1] =  tu.*-2.7E1
  dphi[19,:,1] =  u.*2.7E1-su.*5.4E1-tu.*2.7E1-uu.*2.7E1
  dphi[20,:,1] =  t.*2.7E1-st.*5.4E1-tu.*2.7E1-tt.*2.7E1

  # derivative in y
  dphi[1,:,2] =  s.*1.8E1+t.*1.8E1+u.*1.8E1-st.*2.7E1-su.*2.7E1-tu.*2.7E1-ss.*(2.7E1/2.0)-tt.*(2.7E1/2.0)-
             uu.*(2.7E1/2.0)-1.1E1/2.0
  dphi[2,:,2] =  0.0
  dphi[3,:,2] =  t.*-9.0+tt.*(2.7E1/2.0)+1.0
  dphi[4,:,2] =  0.0
  dphi[5,:,2] =  u.*(-9.0/2.0)+tu.*2.7E1
  dphi[6,:,2] =  u.*(-9.0/2.0)+uu.*(2.7E1/2.0)
  dphi[7,:,2] =  0.0
  dphi[8,:,2] =  0.0
  dphi[9,:,2] =  s.*(-9.0/2.0)+ss.*(2.7E1/2.0)
  dphi[10,:,2] =  s.*(-9.0/2.0)+st.*2.7E1
  dphi[11,:,2] =  s.*(9.0/2.0)+t.*3.6E1+u.*(9.0/2.0)-st.*2.7E1-tu.*2.7E1-tt.*(8.1E1/2.0)-9.0/2.0
  dphi[12,:,2] =  s.*(-4.5E1/2.0)-t.*4.5E1-u.*(4.5E1/2.0)+st.*5.4E1+su.*2.7E1+tu.*5.4E1+ss.*(2.7E1/2.0)+tt.*
             (8.1E1/2.0)+uu.*(2.7E1/2.0)+9.0
  dphi[13,:,2] =  u.*(-4.5E1/2.0)+su.*2.7E1+tu.*2.7E1+uu.*2.7E1
  dphi[14,:,2] =  u.*(9.0/2.0)-uu.*(2.7E1/2.0)
  dphi[15,:,2] =  s.*(-4.5E1/2.0)+st.*2.7E1+su.*2.7E1+ss.*2.7E1
  dphi[16,:,2] =  s.*(9.0/2.0)-ss.*(2.7E1/2.0)
  dphi[17,:,2] =  su.*2.7E1
  dphi[18,:,2] =  u.*2.7E1-su.*2.7E1-tu.*5.4E1-uu.*2.7E1
  dphi[19,:,2] =  su.*-2.7E1
  dphi[20,:,2] =  s.*2.7E1-st.*5.4E1-su.*2.7E1-ss.*2.7E1

  # derivative in z
  dphi[1,:,3] =  s.*1.8E1+t.*1.8E1+u.*1.8E1-st.*2.7E1-su.*2.7E1-tu.*2.7E1-ss.*(2.7E1/2.0)-tt.*(2.7E1/2.0)-
             uu.*(2.7E1/2.0)-1.1E1/2.0
  dphi[2,:,3] =  0.0
  dphi[3,:,3] =  0.0
  dphi[4,:,3] =  u.*-9.0+uu.*(2.7E1/2.0)+1.0
  dphi[5,:,3] =  t.*(-9.0/2.0)+tt.*(2.7E1/2.0)
  dphi[6,:,3] =  t.*(-9.0/2.0)+tu.*2.7E1
  dphi[7,:,3] =  s.*(-9.0/2.0)+su.*2.7E1
  dphi[8,:,3] =  s.*(-9.0/2.0)+ss.*(2.7E1/2.0)
  dphi[9,:,3] =  0.0
  dphi[10,:,3] =  0.0
  dphi[11,:,3] =  t.*(9.0/2.0)-tt.*(2.7E1/2.0)
  dphi[12,:,3] =  t.*(-4.5E1/2.0)+st.*2.7E1+tu.*2.7E1+tt.*2.7E1
  dphi[13,:,3] =  s.*(-4.5E1/2.0)-t.*(4.5E1/2.0)-u.*4.5E1+st.*2.7E1+su.*5.4E1+tu.*5.4E1+ss.*(2.7E1/2.0)+tt.*
             (2.7E1/2.0)+uu.*(8.1E1/2.0)+9.0
  dphi[14,:,3] =  s.*(9.0/2.0)+t.*(9.0/2.0)+u.*3.6E1-su.*2.7E1-tu.*2.7E1-uu.*(8.1E1/2.0)-9.0/2.0
  dphi[15,:,3] =  s.*(-4.5E1/2.0)+st.*2.7E1+su.*2.7E1+ss.*2.7E1
  dphi[16,:,3] =  s.*(9.0/2.0)-ss.*(2.7E1/2.0)
  dphi[17,:,3] =  st.*2.7E1
  dphi[18,:,3] =  t.*2.7E1-st.*2.7E1-tu.*5.4E1-tt.*2.7E1
  dphi[19,:,3] =  s.*2.7E1-st.*2.7E1-su.*5.4E1-ss.*2.7E1
  dphi[20,:,3] =  st.*-2.7E1

  return phi, dphi

end

ext: jl


text:
#=
 PyPlot.jl (Matplotlib) drivers etc. for Pseudospectra.jl

This file is part of Pseudospectra.jl.

Julia implementation
Copyright (c) 2017 Ralph A. Smith

Portions derived from EigTool:
 Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
 of the University of Oxford, and the EigTool Developers. All rights reserved.
 EigTool is maintained on GitHub:  https://github.com/eigtool

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#
module PseudospectraMPL

using ..PyPlot
using Pseudospectra, LinearAlgebra, Printf

export MPLGUIState, psa

# we implement specific methods for these here:
import Pseudospectra: redrawcontour, surfplot, arnoldiplotter!, ewsplotter
import Pseudospectra: plotmode, replzdlg, addmark, fillopts, isheadless
import Pseudospectra: mtxexpsplot, mtxpowersplot
import Pseudospectra: zoomin!, zoomout!, _portrait

# we use these internals here:
import Pseudospectra: vec2ax, expandlevels, oneeigcond, psmode_inv_lanczos
import Pseudospectra: dummyqdlg, replqdlg, transient_bestlb
import Pseudospectra: numrange!

# for convenience (so one could import just this module)
import Pseudospectra: new_matrix, driver!
export new_matrix, driver!

const default_opts = Dict{Symbol,Any}(
    :contourstyle => nothing, :markeig => true,
    :project_fraction => 0.5, :scale_equal => false,
    :no_waitbar => false,
    # good colormaps: "jet" w/ alpha=0.7, "Dark2" w/ alpha=0.8
    :contourkw => Dict{Symbol,Any}(:linewidths => 3,
                                                :extend => "both",
                                                :cmap => "jet"),
    :contourfkw => Dict{Symbol,Any}(:alpha => 0.7,
                                                 :extend => "both",
                                                 :cmap => "jet"),
    :cbarkw => Dict{Symbol,Any}(:extendrect => true),
    :fillcontour => true
)

mutable struct MPLGUIState <: GUIState
    mainph # opaque backend object for main plot object
    # the fig. numbers are redundant, but more convenient
    mainfignum::Int
    secondaryph # opaque backend object for secondary plot object
    secondaryfignum::Int
    drawcmd
    markerlist
    counter
    is_headless::Bool
    do_savefig::Bool
    figfile

    function MPLGUIState(ph=nothing,num=0,specialcmd=nothing;
                           headless=false, savefigs=true)
        if specialcmd == nothing
            dc = headless ? dcheadless : dcinteractive
        else
            dc = specialcmd
        end
        if headless # CHECKME
            ioff()
        else
            ion()
        end
        new(ph,num,nothing,0,dc,[],0,headless,savefigs,"")
    end
end

isheadless(gs::MPLGUIState) = gs.is_headless

function fillopts(gs::MPLGUIState, optsin::Dict{Symbol,Any}=Dict{Symbol,Any}())
    opts = merge(default_opts,optsin)
end

function drawp(gs::MPLGUIState, p, id)
    gs.drawcmd(gs,p,id)
end

# isIJulia() = isdefined(Main, :IJulia) && Main.IJulia.inited

function dcinteractive(gs::MPLGUIState,p,n)
    # IJulia might need special handling here
    nothing
end

function dcheadless(gs::MPLGUIState,p,n)
    if gs.do_savefig
        # TODO: make sure p is current
        if length(gs.figfile) > 0
            fn = gs.figfile
        else
            fn = tempname()
        end
        savefig(fn * ".png")
        println("graph saved in $(fn).png")
    end
    nothing
end

################################################################
# Plotting routines

function redrawcontour(gs::MPLGUIState, ps_data::PSAStruct, opts)
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    ps_dict = ps_data.ps_dict
    levels = expandlevels(zoom.levels)
    Z,x,y = zoom.Z, zoom.x, zoom.y
    if !zoom.computed || isempty(zoom.Z)
        error("should not get here w/ uncomputed zoom")
    end
    eigA = ps_dict[:proj_ews]
        figure(gs.mainfignum) # in case the user looked elsewhere
        clf()
        !isempty(eigA) && plot(real(eigA),imag(eigA),"k.")
        setxylims!(gs.mainph,zoom.ax)
        if isempty(levels)
            contour(x,y,log10.(Z);opts[:contourkw]...)
            contourf(x,y,log10.(Z);opts[:contourfkw]...)
        else
            kwargs = merge(Dict(:levels => levels),opts[:contourkw])
            contour(x,y,log10.(Z);kwargs...)
            if get(opts,:fillcontour,true)
                kwargs = merge(Dict(:levels => levels),opts[:contourfkw])
                contourf(x,y,log10.(Z);kwargs...)
            end
        end
        if get(opts,:colorbar,true)
            colorbar(;opts[:cbarkw]...)
        end
        if get(opts,:showimagax,false)
            plot([0,0],zoom.ax[3:4],"k-")
        end
        if get(opts,:showunitcircle,false)
            plot(cos.((0:0.01:2)*π),sin.((0:0.01:2)*π),"k-")
        end
        if get(opts,:showfov,false)
            if isempty(get(ps_dict,:schur_mtx,[]))
                @warn("showfov set in opts but unavailable")
            else
                if isempty(get(ps_dict,:fov,[]))
                    numrange!(ps_data,get(opts,:fov_npts,20))
                end
                plot(real.(ps_dict[:fov]),imag.(ps_dict[:fov]),"k--")
            end
        end
    drawp(gs,gs.mainph,1)
    nothing
end

function surfplot(gs::MPLGUIState, ps_data::PSAStruct, opts)
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    zoom.computed || return
    nx,ny = size(zoom.Z)
    # TODO: allow for user options
    # TODO: control plotting of mesh, e.g.
    # line_freq = floor(Int,min(nx,ny)/10)
    # TODO: add cylinders for eigenvalues, like upstream
        figure()
        surf(zoom.x,zoom.y,-log10.(zoom.Z))
    drawp(gs,gs.mainph,2)
end

function arnoldiplotter!(gs::MPLGUIState,old_ax,opts,dispvec,infostr,ews,shifts,
                         state)
    if gs.mainph == nothing
        ax = vec2ax(dispvec)
    else
        if get(opts,:ARPACK_auto_ax,true)
            ax = vec2ax(dispvec)
            if isempty(old_ax)
                append!(old_ax,getxylims(gs.mainph))
            else
                copy!(old_ax,getxylims(gs.mainph))
                ax[1] = min(ax[1],old_ax[1])
                ax[3] = min(ax[3],old_ax[3])
                ax[2] = max(ax[2],old_ax[2])
                ax[4] = max(ax[4],old_ax[4])
            end
            # ps_data.zoom_list[ps_data.zoom_pos].ax = ax
        else
            @warn("non-auto axes not implemented")
            # CHECKME: are we sure this is ready?
            # ax = ps_data.zoom_list[ps_data.zoom_pos].ax
        end
    end # ax setting

        # figure(mainfigure)
        print(".") # sometimes apparently needed to force a redraw
        clf()
        # if get(opts,:ProgPSA,false)
        #  TODO: draw pseudospectra
        # end
        plot(real.(shifts),imag.(shifts),"r+")
        plot(real.(ews),imag.(ews),"k.")
        axis(ax)
        #text(?,?,infostr,fontsize=12,fontweight=bold)
        title(infostr)
        drawp(gs,gs.mainph,1)
        draw()
    #   sleep(0.01)
    return nothing
end

"""
    selectfig(gs,main::Bool)

start a new figure, or attach to an existing one.
"""
function selectfig(gs::MPLGUIState,main::Bool)
    if main
        if gs.mainfignum > 0
            figure(gs.mainfignum)
        else
            fh = figure()
            gs.mainph = fh
            gs.mainfignum = fh.number
        end
    else
        if gs.secondaryfignum > 0
            figure(gs.secondaryfignum)
        else
            fh = figure()
            gs.secondaryph = fh
            gs.secondaryfignum = fh.number
        end
    end
end

function ewsplotter(gs::MPLGUIState, ews::Vector, zoom)
    selectfig(gs,true)
    clf()
    plot(real(ews),imag(ews),"k.")
    isempty(zoom.ax) && (zoom.ax = vec2ax(ews))
    setxylims!(gs.mainph,zoom.ax)
    drawp(gs,gs.mainph,1)
end

function plotmode(gs::MPLGUIState,z,A,U,pseudo::Bool,approx::Bool,verbosity)
    m,n = size(A)
    if pseudo
        modestr = "Pseudo"
        niter = 20
        q = normalize!(randn(n) + randn(n)*im)
        Ss,q = psmode_inv_lanczos(A,q,z,1e-15,niter)
        the_col = "magenta"
    else
        modestr = "Eigen"
        Ss,q = oneeigcond(A,z,verbosity,
                          dlg=ifelse(isheadless(gs),dummyqdlg,replqdlg))
        the_cond = Ss
        the_col = "cyan"
    end
    # FIXME: bail out on error

    # transform to account for Schur decomposition
    q = U * q

    # upstream tries to fetch psmode_x_points_ from base workspace
    x = collect(1:length(q))
    showlog = sum(abs.(diff(q))) * (1/length(q)) >= 10*eps()
    if pseudo
        the_str = @sprintf("\$\\epsilon = %15.2e\$",Ss[end])
    else
        the_rel = approx ? "\\approx" : "="
        the_str = "\$\\kappa(\\lambda) " * the_rel *
            @sprintf("%15.2e\$",the_cond)
    end
    λ_str = @sprintf("%12g%+12gi",real(z),imag(z))

    selectfig(gs, false)
    clf()
    showlog && subplot(2,1,1)

    plot(x,abs.(q),color="black",label="abs")
    plot(x,-abs.(q),color="black")
    plot(x,real.(q),color=the_col,label="realpt")
    legend()
    xlim(x[1],x[end])
        # DEVNOTE: upstream puts the_str in the plot area (often obscured)
        #=
        y0,y1 = ylim()
        dx = x[end]-x[1]
        dy = y1-y0
        # upstream also has fontweight="bold", but that's hard w/ math mode
        text(x[1]+dx/40,y0+dy/10,the_str,fontsize=12)
        =#
    title("$(modestr)mode: " *
          "\$\\lambda=" * λ_str * "\$\n" * the_str)
    if showlog
        subplot(2,1,2)
        semilogy(x,abs.(q),"k")
        xlim(x[1],x[end])
    end
    drawp(gs,gs.mainph,2)
end

function mtxpowersplot(gs::MPLGUIState,ps_data::PSAStruct,nmax=50;gradual=false,
                       lbmethod=:none,lbdk=0.25)

    stop_trans = false

    sp_step = 1 # basic stepsize

    A = ps_data.input_matrix^sp_step

    trans = zeros(nmax+1)
    the_pwr = similar(trans)

    trans[1] = 1.0
    the_pwr[1] = 0.0
    ax = [0,20,0,1.5]
    ax2 = copy(ax)
    ax2[3] = 1e-2
    pos = 2
    mtx = Matrix(1.0I,size(A)...)
    max_tp = -Inf
    min_tp = Inf
    ax_factor = 3

    # spectral radius
    alp = maximum(abs.(ps_data.ps_dict[:ews]))
    newfig = true

    function doplot()
        # if plotting gradually, control axis limits to avoid jumpiness
        if newfig
            selectfig(gs,false)
            newfig = false
        end
            subplot(2,1,1)
            gradual && cla()
            plot(the_pwr[1:pos],trans[1:pos])
            plot(the_pwr[1:pos],alp.^the_pwr[1:pos],"g--")
            ylabel("\$\\|A^k\\|\$")
            if gradual
                xlim(ax[1],ax[2])
                ylim(ax[3],ax[4])
            end
            subplot(2,1,2)
            gradual && cla()
            semilogy(the_pwr[1:pos],trans[1:pos])
            plot(the_pwr[1:pos],alp.^the_pwr[1:pos],"g--")
            xlabel("k")
            if gradual
                xlim(ax[1],ax[2])
                ylim(ax2[3],ax2[4])
            end
            drawp(gs,gs.mainph,2)
    end
    while !stop_trans
        mtx = A * mtx
        trans[pos] = norm(mtx)
        the_pwr[pos] = sp_step*(pos-1)
        max_tp = max(max_tp,trans[pos])
        min_tp = min(min_tp,trans[pos])

        if the_pwr[pos] > ax[2]*0.9
            ax[2] *= ax_factor
            ax2[2] *= ax_factor
        end
        if max_tp > ax[4] * 0.9
            ax[4] = max(max_tp * 1.5,ax_factor * ax[4])
        end
        if max_tp > ax2[4] * 0.9
            ax2[4] = ax2[4]^ax_factor
        end
        if (log10(min_tp)-log10(ax2[3])) / (log10(ax2[4]) - log10(ax2[3])) < 0.1
            ax2[3] = ax2[3]^ax_factor
        end

        if gradual
            doplot()
            sleep(0.02)
        end
        if (max_tp > 1e130) || ((min_tp < 1e-130) && (min_tp != 0))
            @warn("stopping: numbers going out of range")
            stop_trans = true
        end
        if min_tp == 0
            @warn("stopping: exactly zero matrix")
            stop_trans = true
        end
        (pos > nmax) && (stop_trans = true)
        pos += 1
    end
    if !gradual
        pos -= 1
        doplot()
    end
    if lbmethod != :none
        pts,bnds,sel_pt = transient_bestlb(ps_data,:powers,0:nmax,
                                           method=lbmethod, dk=lbdk)

        plot(pts,bnds,"m")
        # upstream has something like
        # figure(gs.mainfignum)
        # scatter(sel_pt[:,1],sel_pt[:,2])
    end
end

function mtxexpsplot(gs::MPLGUIState,ps_data::PSAStruct,dt=0.1,nmax=50;
                     gradual=false, lbmethod=:none)

    stop_trans = false
    A = ps_data.input_matrix
    eAdt = exp(dt*A)

    trans = zeros(nmax+1)
    the_time = similar(trans)

    trans[1] = 1.0
    the_time[1] = 0.0
    ax = [0.0,20.0*dt,0,2.0]
    ax2 = copy(ax)
    ax2[3] = 1e-2
    pos = 2
    eAt = Matrix(1.0I,size(A)...)
    max_tp = -Inf
    min_tp = Inf
    ax_factor = 3

    # get spectral abscissa
    alp = maximum(real.(ps_data.ps_dict[:ews]))
    newfig = true

    if lbmethod != :none
        pts,bnds,sel_pt = transient_bestlb(ps_data,:exp,dt*(0:nmax),
                                           method=lbmethod)
    else
        bnds = zeros(0)
    end
    function doplot()
        if newfig
            selectfig(gs,false)
            clf()
            newfig = false
        end
            subplot(2,1,1)
            gradual && cla()
            plot(the_time[1:pos],trans[1:pos])
            plot(the_time[1:pos],exp.(alp*the_time[1:pos]),"g--")
            if gradual
                xlim(ax[1],ax[2])
                ylim(ax[3],ax[4])
            end
            ylabel("\$\\|e^{At}\\|\$")
            subplot(2,1,2)
            gradual && cla()
            semilogy(the_time[1:pos],trans[1:pos])
            plot(the_time[1:pos],exp.(alp*the_time[1:pos]),"g--")
            xlabel("t")
            if !isempty(bnds)
                plot(pts[1:pos],bnds[1:pos],"m")
            end
            if gradual
                xlim(ax[1],ax[2])
                ylim(ax2[3],ax2[4])
            end
        drawp(gs,gs.mainph,2)
    end

    while !stop_trans
        eAt = eAt * eAdt
        trans[pos] = norm(eAt)
        the_time[pos] = (pos-1)*dt
        max_tp = max(max_tp,trans[pos])
        min_tp = min(min_tp,trans[pos])

        if the_time[pos] > ax[2]*0.9
            ax[2] *= ax_factor
            ax2[2] *= ax_factor
        end
        if max_tp > ax[4] * 0.9
            ax[4] = max(max_tp * 1.5,ax_factor * ax[4])
        end
        if max_tp > ax2[4] * 0.9
            ax2[4] = ax2[4]^ax_factor
        end
        if (log10(min_tp)-log10(ax2[3])) / (log10(ax2[4]) - log10(ax2[3])) < 0.1
            ax2[3] = ax2[3]^ax_factor
        end
        if gradual
            doplot()
            sleep(0.02)
        end

        if (max_tp > 1e130) || ((min_tp < 1e-130) && (min_tp != 0))
            @warn("stopping: numbers going out of range")
            stop_trans = true
        end
        (pos > nmax) && (stop_trans = true)
        pos += 1
    end
    if !gradual
        pos -= 1
        doplot()
    end
end

function zoomin!(gs::MPLGUIState, ps_data::PSAStruct,
                 optsin=Dict{Symbol,Any}(); zkw=zeros(0))
    opts = fillopts(gs,optsin)
    if isempty(zkw)
        # User might get z from some other figure, but so what?
        z = replzdlg(gs)
    else
        z = zkw
    end
    figure(gs.mainfignum)
    ax = getxylims(gs.mainph)
    res = zoomin!(ps_data,z,ax)
    if res < 0
        @warn("unable to zoom based on requested point")
    elseif res == 0
        # redraw existing portrait
        redrawcontour(gs, ps_data, opts)
    else
        # compute and draw new portrait
        driver!(ps_data, opts, gs)
    end
    res
end

function zoomout!(gs::MPLGUIState, ps_data::PSAStruct,
                 optsin=Dict{Symbol,Any}(); include_fov=false, zkw=zeros(0))
    opts = fillopts(gs,optsin)
    if isempty(zkw)
        z = replzdlg(gs)
    else
        z = zkw
    end
    figure(gs.mainfignum)
    ax = getxylims(gs.mainph); # current plot axes
    res = zoomout!(ps_data,z,ax,include_fov=include_fov)
    if res < 0
        @warn("unable to zoom based on requested point")
    elseif res == 0
        # redraw existing portrait
        redrawcontour(gs, ps_data, opts)
    else
        # compute and draw new portrait
        driver!(ps_data, opts, gs)
    end
    res
end

################################################################
# Wrappers

"""
    psa(A [,opts]) -> ps_data, graphics_state

Compute and plot pseudospectra of a matrix.
This is a rudimentary driver for the Pseudospectra package.
"""
function psa(A::AbstractMatrix, opts::Dict{Symbol,Any}=Dict{Symbol,Any}())
    allopts = merge(default_opts,opts)
    ps_data = new_matrix(A,allopts)
    fh = figure()
    mainfignum = fh.number
    gs = MPLGUIState(fh,mainfignum)
    driver!(ps_data,allopts,gs=gs)
    ps_data, gs
end

function _portrait(xs,ys,Z,eigA)
    p = contour(xs,ys,log10.(Z),cmap="eigtool")
    plot(real(eigA), imag(eigA), "k.", # label="eigvals",
         markersize=5)
    colorbar()
    return p
end
################################################################
# Utilities

function setxylims!(ph,ax)
    # FIXME: this should insure reference to ph
    axis(ax)
end

function getxylims(ph)
    # FIXME: this should insure reference to ph
    collect(axis())
end

function replzdlg(gs::MPLGUIState; what="a z-value")
    x,y = NaN,NaN
    println("click to select $what in main plot frame")
    figure(gs.mainfignum) # in case the user looked elsewhere
    x,y = ginput()[1]
    z = x + y*im
    return z
end

function addmark(gs::MPLGUIState,z,mykey)
    x,y = real(z),imag(z)
    figure(gs.mainfignum)
    #=
    l = gs.markerlist
    # remove any existing mark of this sort
    # WARNING: logic assumes zero or one mark per class
    for i in eachindex(l)
        if l[i][:class] == mykey
            l[i][:handle][1][:remove]()
            deleteat!(l,i)
            break
        end
    end
    =#
    if mykey == :pseudo
        ltype = "mo"
    elseif mykey == :eigen
        ltype = "co"
    else
        ltype = "k*"
    end
    ph = plot([x],[y],ltype)
    # push!(gs.markerlist,Dict(:class => mykey, :handle => ph))
    nothing
end

function et_cmap()
    cm = [
        0   1.00000000000000   0.10000000000000
        0   0.92500000000000   0.17500000000000
        0   0.85000000000000   0.25000000000000
        0   0.77500000000000   0.32500000000000
        0   0.70000000000000   0.40000000000000
        0   0.62500000000000   0.47500000000000
        0   0.55000000000000   0.55000000000000
        0   0.47500000000000   0.62500000000000
        0   0.40000000000000   0.70000000000000
        0   0.32500000000000   0.77500000000000
        0   0.25000000000000   0.85000000000000
        0   0.17500000000000   0.92500000000000
        0   0.10000000000000   1.00000000000000
        0.06153846153846   0.09230769230769   1.00000000000000
        0.12307692307692   0.08461538461538   1.00000000000000
        0.18461538461538   0.07692307692308   1.00000000000000
        0.24615384615385   0.06923076923077   1.00000000000000
        0.30769230769231   0.06153846153846   1.00000000000000
        0.36923076923077   0.05384615384615   1.00000000000000
        0.43076923076923   0.04615384615385   1.00000000000000
        0.49230769230769   0.03846153846154   1.00000000000000
        0.55384615384615   0.03076923076923   1.00000000000000
        0.61538461538462   0.02307692307692   1.00000000000000
        0.67692307692308   0.01538461538462   1.00000000000000
        0.73846153846154   0.00769230769231   1.00000000000000
        0.80000000000000                  0   1.00000000000000
        0.80769230769231                  0   0.92307692307692
        0.81538461538462                  0   0.84615384615385
        0.82307692307692                  0   0.76923076923077
        0.83076923076923                  0   0.69230769230769
        0.83846153846154                  0   0.61538461538462
        0.84615384615385                  0   0.53846153846154
        0.85384615384615                  0   0.46153846153846
        0.86153846153846                  0   0.38461538461538
        0.86923076923077                  0   0.30769230769231
        0.87692307692308                  0   0.23076923076923
        0.88461538461538                  0   0.15384615384615
        0.89230769230769                  0   0.07692307692308
        0.90000000000000                  0                  0
        0.86923076923077                  0                  0
        0.83846153846154                  0                  0
        0.80769230769231                  0                  0
        0.77692307692308                  0                  0
        0.74615384615385                  0                  0
        0.71538461538462                  0                  0
        0.68461538461538                  0                  0
        0.65384615384615                  0                  0
        0.62307692307692                  0                  0
        0.59230769230769                  0                  0
        0.56153846153846                  0                  0
        0.53076923076923                  0                  0
        0.50000000000000                  0                  0
        0.54166666666667   0.05000000000000                  0
        0.58333333333333   0.10000000000000                  0
        0.62500000000000   0.15000000000000                  0
        0.66666666666667   0.20000000000000                  0
        0.70833333333333   0.25000000000000                  0
        0.75000000000000   0.30000000000000                  0
        0.79166666666667   0.35000000000000                  0
        0.83333333333333   0.40000000000000                  0
        0.87500000000000   0.45000000000000                  0
        0.91666666666667   0.50000000000000                  0
        0.95833333333333   0.55000000000000                  0
        1.00000000000000   0.60000000000000                  0
    ]
    nc = size(cm,1)
    xvals = collect(range(0.0, stop=1.0, length=nc))
    r = [(xvals[i],cm[i,1],cm[i,1]) for i in 1:nc]
    g = [(xvals[i],cm[i,2],cm[i,2]) for i in 1:nc]
    b = [(xvals[i],cm[i,3],cm[i,3]) for i in 1:nc]
    PyPlot.ColorMap("eigtool",r,g,b,
                    Array{Tuple{Float64,Float64,Float64}}(undef, 0),
                    256,1.0)
end

PyPlot.register_cmap("eigtool",et_cmap())

end # module


text:
include("debug.jl")
include("fpump.jl")

function create_root_model!(optimizer::MOI.ModelLike, jp::JuniperProblem; fix_start=false)
    ps = jp.options.log_levels

    jp.model = MOI.instantiate(jp.nl_solver, with_bridge_type=Float64)
    index_map = MOI.copy_to(jp.model, IntegerRelaxation(optimizer); copy_names=false)
    # all continuous we solve relaxation first
    jp.x = [index_map[vi] for vi in MOI.get(optimizer, MOI.ListOfVariableIndices())]
end

function fix_primal_start!(jp::JuniperProblem)
    lb = jp.l_var
    ub = jp.u_var
    x = jp.x
    for i=1:jp.num_var
        if jp.var_type[i] != :Cont && lb[i] <= jp.primal_start[i] <= ub[i]
            set_bounds(jp.model, jp.x[i], jp.primal_start[i], jp.primal_start[i])
        else
            MOI.set(jp.model, MOI.VariablePrimalStart(), jp.x[i], jp.primal_start[i])
        end
    end
end

function unfix_primal_start!(jp::JuniperProblem)
    lb = jp.l_var
    ub = jp.u_var
    x = jp.x
    for i=1:jp.num_var
        set_bounds(jp.model, jp.x[i], lb[i], ub[i])
        MOI.set(jp.model, MOI.VariablePrimalStart(), jp.x[i], jp.primal_start[i])
    end
end

function solve_root_incumbent_model(jp::JuniperProblem)
    MOI.optimize!(jp.model)
    status = MOI.get(jp.model, MOI.TerminationStatus())
    incumbent = nothing
    ps = jp.options.log_levels
    if state_is_optimal(status; allow_almost=jp.options.allow_almost_solved)
        # set incumbent
        objval = MOI.get(jp.model, MOI.ObjectiveValue())
        solution = MOI.get(jp.model, MOI.VariablePrimal(), jp.x)
        if are_type_correct(solution, jp.var_type, jp.disc2var_idx, jp.options.atol)
            if only_almost_solved(status) && jp.options.allow_almost_solved_integral
                @warn "Start value incumbent only almost locally solved. Disable with `allow_almost_solved_integral=false`"
            end

            incumbent = Incumbent(objval, solution, only_almost_solved(status))
            (:All in ps || :Info in ps) && println("Incumbent using start values: $objval")
        end
    else
        (:All in ps || :Info in ps) && println("Start values are not feasible.")
    end
    return incumbent
end

function solve_root_model!(jp::JuniperProblem)
    MOI.optimize!(jp.model)
    status = MOI.get(jp.model, MOI.TerminationStatus())
    jp.relaxation_status = status
    restarts = 0
    max_restarts = jp.options.num_resolve_root_relaxation
    jp.options.debug && debug_init(jp.debugDict)
    while !state_is_optimal(jp.relaxation_status; allow_almost=jp.options.allow_almost_solved) &&
        restarts < max_restarts && time()-jp.start_time < jp.options.time_limit

        # TODO freemodel for Knitro
        restart_values = generate_random_restart(jp)
        jp.options.debug && debug_restart_values(jp.debugDict,restart_values)
        MOI.set(jp.model, MOI.VariablePrimalStart(), jp.x, restart_values)
        MOI.optimize!(jp.model)
        jp.relaxation_status = MOI.get(jp.model, MOI.TerminationStatus())
        restarts += 1
    end
    if only_almost_solved(jp.relaxation_status)
        @warn "The relaxation is only almost solved."
    end
    return restarts
end


text:
"Supertype of identity maps."
abstract type IdentityMap{T} <: Map{T} end

IdentityMap(n::Int) = DynamicIdentityMap(n)
IdentityMap() = StaticIdentityMap()
IdentityMap(::Val{N}) where {N} = StaticIdentityMap(Val(N))

IdentityMap{T}(n::Int) where {T} = DynamicIdentityMap{T}(n)
IdentityMap{T}(n::Int) where {T<:StaticTypes} = StaticIdentityMap{T}()
IdentityMap{T}(::Val{N}) where {N,T} = StaticIdentityMap{T}(Val(N))
IdentityMap{T}() where {T} = StaticIdentityMap{T}()

applymap(map::IdentityMap, x) = x
applymap!(y, map::IdentityMap, x) = y .= x

inverse(m::IdentityMap) = m
inverse(m::IdentityMap, x) = x

islinear(::IdentityMap) = true
isreal(::IdentityMap{T}) where {T} = isreal(T)

isidentity(::IdentityMap) = true
isidentity(m::Map{T}) where {T} = m == StaticIdentityMap{T}()

mapsize(m::IdentityMap{T}) where {T<:Number} = ()
mapsize(m::IdentityMap{T}) where {T} = (euclideandimension(T),euclideandimension(T))

matrix(m::IdentityMap) = identitymatrix(m)
vector(m::IdentityMap) = zerovector(m)

jacobian(m::IdentityMap) = ConstantMap(matrix(m))
jacobian(m::IdentityMap, x) = matrix(m)

jacdet(m::IdentityMap, x) = 1

determinantmap(m::IdentityMap{T}) where {T} = UnityMap{T,prectype(T)}()

mapcompose(m1::IdentityMap) = m1
mapcompose(m1::IdentityMap, maps...) = mapcompose(maps...)
mapcompose2(m1, m2::IdentityMap, maps...) = mapcompose(m1, maps...)

show(io::IO, m::IdentityMap{T}) where {T} = print(io, "x -> x")
Display.object_parentheses(m::IdentityMap) = true

"The identity map for variables of type `T`."
struct StaticIdentityMap{T} <: IdentityMap{T}
end

StaticIdentityMap() = StaticIdentityMap{Float64}()
StaticIdentityMap(::Val{N}) where {N} = StaticIdentityMap{SVector{N,Float64}}()

StaticIdentityMap{T}(n::Int) where {T} =
    (@assert n == euclideandimension(T); StaticIdentityMap{T}())
StaticIdentityMap{T}(::Val{N}) where {N,T} =
    (@assert N == euclideandimension(T); StaticIdentityMap{T}())

similarmap(m::StaticIdentityMap, ::Type{T}) where {T<:StaticTypes} = StaticIdentityMap{T}()
similarmap(m::StaticIdentityMap, ::Type{T}) where {T} =
    DynamicIdentityMap{T}(euclideandimension(T))

convert(::Type{StaticIdentityMap{T}}, ::StaticIdentityMap) where {T} = StaticIdentityMap{T}()

==(m1::StaticIdentityMap, m2::StaticIdentityMap) = true
hash(m::StaticIdentityMap, h::UInt) = hash("StaticIdentityMap", h)

"Identity map with dynamic size determined by a dimension field."
struct DynamicIdentityMap{T} <: IdentityMap{T}
    dimension   ::  Int
end

const EuclideanIdentityMap{N,T} = StaticIdentityMap{SVector{N,T}}
const VectorIdentityMap{T} = DynamicIdentityMap{Vector{T}}

DynamicIdentityMap(dimension::Int) = VectorIdentityMap(dimension)
VectorIdentityMap(dimension::Int) = VectorIdentityMap{Float64}(dimension)

mapsize(m::DynamicIdentityMap) = (m.dimension, m.dimension)

similarmap(m::DynamicIdentityMap, ::Type{T}) where {T} =
    DynamicIdentityMap{T}(m.dimension)
similarmap(m::DynamicIdentityMap, ::Type{T}) where {T<:StaticTypes} =
    StaticIdentityMap{T}()

==(m1::DynamicIdentityMap, m2::DynamicIdentityMap) = m1.dimension == m2.dimension
hash(m::DynamicIdentityMap, h::UInt) = hashrec("DynamicIdentityMap", m.dimension, h)

"The supertype of constant maps from `T` to `U`."
abstract type ConstantMap{T,U} <: TypedMap{T,U} end

applymap(m::ConstantMap, x) = constant(m)

isconstant(m::AbstractMap) = false
isconstant(m::ConstantMap) = true

isreal(m::ConstantMap{T,U}) where {T,U} =
    isreal(T) && isreal(U) && isreal(constant(m))

mapsize(m::ConstantMap) = _constant_mapsize(m, constant(m))
_constant_mapsize(m::ConstantMap{T,U}, c) where {T<:Number,U<:Number} = ()
_constant_mapsize(m::ConstantMap{T,U}, c) where {T<:Number,U} = (length(c),)
_constant_mapsize(m::ConstantMap{T,U}, c) where {T,U<:Number} = (1,euclideandimension(T))
_constant_mapsize(m::ConstantMap{T,U}, c) where {T,U} = (length(c), euclideandimension(T))

matrix(m::ConstantMap) = zeromatrix(m)
vector(m::ConstantMap) = constant(m)

jacobian(m::ConstantMap{T}) where {T} = ConstantMap{T}(matrix(m))
jacobian(m::ConstantMap, x) = matrix(m)

jacdet(::ConstantMap, x) = 0

determinantmap(m::ConstantMap{T}) where {T} = ConstantMap{T}(det(constant(m)))
absmap(m::ConstantMap{T}) where {T} = ConstantMap{T}(abs(constant(m)))

==(m1::ConstantMap, m2::ConstantMap) = constant(m1)==constant(m2)
hash(m::ConstantMap, h::UInt) = hashrec("ConstantMap", constant(m), h)

similarmap(m::ConstantMap, ::Type{T}) where {T} = ConstantMap{T}(constant(m))
similarmap(m::ConstantMap, ::Type{T}, ::Type{U}) where {T,U} = ConstantMap{T,U}(m.c)

ConstantMap() = ConstantMap{Float64}()
ConstantMap(c) = FixedConstantMap(c)
ConstantMap{T}() where {T} = UnityMap{T}()
ConstantMap{T}(c) where {T} = FixedConstantMap{T}(c)
ConstantMap{T,U}() where {T,U} = UnityMap{T,U}()
ConstantMap{T,U}(c) where {T,U} = FixedConstantMap{T,U}(c)

show(io::IO, m::ConstantMap{T}) where {T} = print(io, "x -> $(constant(m))")
Display.object_parentheses(m::ConstantMap) = true

"The zero map `f(x) = 0`."
struct ZeroMap{T,U} <: ConstantMap{T,U}
end
ZeroMap{T}() where {T} = ZeroMap{T,T}()
constant(m::ZeroMap{T,U}) where {T,U} = zero(U)
similarmap(m::ZeroMap{S,U}, ::Type{T}) where {T,S,U} = ZeroMap{T,U}()
similarmap(m::ZeroMap, ::Type{T}, ::Type{U}) where {T,U} = ZeroMap{T,U}()

"The unity map `f(x) = 1`."
struct UnityMap{T,U} <: ConstantMap{T,U}
end
UnityMap{T}() where {T} = UnityMap{T,real(numtype(T))}()
constant(m::UnityMap{T,U}) where {T,U} = one(U)
similarmap(m::UnityMap{S,U}, ::Type{T}) where {T,S,U} = UnityMap{T,U}()
similarmap(m::UnityMap, ::Type{T}, ::Type{U}) where {T,U} = UnityMap{T,U}()

"The constant map `f(x) = c`."
struct FixedConstantMap{T,U} <: ConstantMap{T,U}
    c   ::  U
end
FixedConstantMap{T}(c::U) where {T,U} = FixedConstantMap{T,U}(c)
FixedConstantMap(c::T) where {T} = FixedConstantMap{T}(c)
constant(m::FixedConstantMap) = m.c


text:
export Schedule, step!, clear!, reset!, run!, reshuffle!
export is_complete
export schedule_once!, schedule_once_in!

import Base.Collections: heappush!, heappop!

const EPOCH              = zero(Float64)
const BEFORE_SIMULATION  = EPOCH - 1.0
const AFTER_SIMULATION   = Inf64
const EPOCH_PLUS_EPSILON = eps(EPOCH)

# For values greater than 9.007199254740991E15, the epsilon is > 1.0.
assert(eps(9.007199254740991E15) <= 1.0 && eps(9.007199254740992E15) > 1.0)
const MAXIMUM_TIME = 9.007199254740992E15

type Schedule
    time::Float64
    steps::Int64
    queue::Vector{Event}  # Heap-sorted
end

Schedule() = Schedule(BEFORE_SIMULATION, 0, Vector{Event}())

function push!(schedule::Schedule, actionable, at::ActivationPoint)
    # You cannot schedule an action to occur at precisely the current time.
    # For example, imagine an activated action that reschedules itself for
    # now. Time would never progress. The same action would repeat, ad
    # infinitum. There must be some epsilon between now and activation.
    if schedule.time == at.time && schedule.time != AFTER_SIMULATION
        at = at + eps(at.time)
    end

    at.time >= AFTER_SIMULATION && return false

    at.time < EPOCH && throw(ArgumentError("Scheduled before epoch"))

    isnan(at.time) && throw(ArgumentError("Scheduled at NaN"))

    at.time < schedule.time && throw(ArgumentError("Scheduled in past"))

    heappush!(schedule.queue, Event(at, actionable))

    true
end

function push!(action::Function, schedule::Schedule, at::ActivationPoint)
    push!(schedule, action, at)
end

function step!(schedule::Schedule, env)
    (schedule.time == AFTER_SIMULATION || isempty(schedule)) && return false

    next_time = schedule.queue[1].activate_at.time
    schedule.time = next_time

    while next_time == schedule.time
        event = heappop!(schedule.queue)
        event.action(env, schedule)

        if isempty(schedule)
            break
        else
            next_time = schedule.queue[1].activate_at.time
        end
    end

    schedule.steps += 1

    return true
end

function reshuffle!(schedule::Schedule)
    schedule.queue = [
        Event(ActivationPoint(event.activate_at), event.action)  # new jitter
        for event in schedule.queue
    ]
end

function run!(schedule, env, n::Int64=Int64(-1))
    isempty(schedule) && throw(ArgumentError("Schedule already exhausted"))

    if n > 0
        while step!(schedule, env) && schedule.steps < n
        end
    else
        while step!(schedule, env)
        end
    end

    schedule
end

function run!(f::Function, schedule, env)
    isempty(schedule) && throw(ArgumentError("Schedule already exhausted"))

    while step!(schedule, env)
        if f(env, schedule) == false
            break
        end
    end

    schedule
end

length(schedule::Schedule) = length(schedule.queue)

function show(io::IO, s::Schedule)
    print(io, "Schedule at time=$(s.time) and steps=$(s.steps) ")
    print(io, "with $(length(s)) events.")
end

function clear!(schedule::Schedule)
    n = length(schedule)
    empty!(schedule.queue)
    n
end

function reset!(schedule::Schedule)
    schedule.time = BEFORE_SIMULATION
    schedule.steps = 0
    clear!(schedule)
end

isempty(schedule::Schedule) = isempty(schedule.queue)
is_complete(schedule::Schedule) = isempty(schedule)

function merge!(schedule::Schedule, other::Schedule)
    if other.queue[1].activate_at.time < schedule.time
        error("Merge over items already in the past.")
    end

    for item in other.queue
        heappush!(schedule.queue, item)
    end

    schedule
end

function schedule_once!(schedule::Schedule, actionable,
                        time::AbstractFloat=1.0, order::Integer=1)
    push!(schedule,  actionable, ActivationPoint(time, order))
end

function schedule_once!(f::Function, schedule::Schedule,
                        time::AbstractFloat=1.0, order::Integer=1)
    schedule_once!(schedule, f, time, order)
end

function schedule_once_in!(schedule::Schedule, actionable,
                           Δ::AbstractFloat=1.0, order::Integer=1)
    push!(schedule, actionable, ActivationPoint(schedule.time + Δ, order))
end

function schedule_once_in!(f::Function, schedule::Schedule,
                           Δ::AbstractFloat=1.0, order::Integer=1)
    schedule_once_in!(schedule, f, Δ, order)
end


text:
using DataFrames
using Statistics
using ThreadsX
using Statistics
using LinearAlgebra

"""
The `UpdateRule` abstract type defines mechanisms to update the best composition estimate between
iteration steps.
"""
abstract type UpdateRule end

"""
The `NaiveUpdateRule` implements the 'method of successive approximations' to update
the composition between iteration steps.
"""
struct NaiveUpdateRule <: UpdateRule end

"""
    update(
        ::NaiveUpdateRule,
        prevcomp::Material,
        measured::Vector{KRatio},
        zafs::Dict{Element, Float64}
    )::Dict{Element,Float64}

Determine the next estimate of the composition that brings the estimate k-ratios closer to measured.
"""
update(
    ::NaiveUpdateRule, #
    prevcomp::Material, #
    measured::Vector{KRatio}, #
    zafs::Dict{Element,Float64}, #
)::Dict{Element,Float64} =
    Dict(kr.element => value(nonnegk(kr)) * value(kr.standard[kr.element]) / zafs[kr.element] for kr in measured)

function reset(::NaiveUpdateRule) end

"""
The `WegsteinUpdateRule` implements the very effective method of Reed and Mason (S.J.B. Reed and P.K. Mason,
Transactions ofthe Second National Conference on Electron Microprobe Analysis, Boston, 1967.) for updating
the composition estimate between iteration steps.
"""
mutable struct WegsteinUpdateRule <: UpdateRule
    fallback::UpdateRule
    prevc::Union{Material,Missing}
    prevfn::Union{Dict{Element,Float64},Missing}
    WegsteinUpdateRule() = new(NaiveUpdateRule(), missing, missing)
end

function update( #
    weg::WegsteinUpdateRule,
    prevcomp::Material,
    measured::Vector{KRatio},
    zafs::Dict{Element,Float64},
)::Dict{Element,Float64}
    # Default to naive
    cnp1 = update(weg.fallback, prevcomp, measured, zafs)
    fn = Dict(kr.element => value(kr.standard[kr.element]) / zafs[kr.element] for kr in measured)
    if !ismissing(weg.prevc)
        cn, cnm1, fnm1 = prevcomp, weg.prevc, weg.prevfn
        for mkr in measured
            elm, km = mkr.element, value(nonnegk(mkr))
            if km > 0
                δfδc = (fn[elm] - fnm1[elm]) / (value(cn[elm]) - value(cnm1[elm]))
                den = 1.0 - km * δfδc
                if (abs(δfδc) < 10.0) && (abs(den) > 0.2)
                    Δc = (km * fn[elm] - cn[elm]) / den # Wegstein
                    cnp1[elm] = cn[elm] + Δc
                end
            end
        end
    end
    weg.prevc, weg.prevfn = prevcomp, fn
    return cnp1
end
"""
    reset(weg::WegsteinUpdateRule)

Restart the WegsteinUpdateRule accumulators.
"""
function reset(weg::WegsteinUpdateRule)
    weg.prevc = missing
    weg.prevfn = missing
end

"""
The `RecordingUpdateRule` wraps other `UpdateRule` instances to provide diagnostics which may be tabulated
using asa(DataFrame, rur::RecordingUpdateRule) or plotted using Gadfly's plot(rur::RecordingUpdateRule).
"""
struct RecordingUpdateRule <: UpdateRule
    base::UpdateRule
    comps::Vector{Dict{Element,Float64}}
    """
        RecordingUpdateRule(ur::UpdateRule)

    Wrap an UpdateRule instance with diagnostic recorders.
    """
    RecordingUpdateRule(ur::UpdateRule) = new( ur, Vector{Dict{Element,Float64}}())
end
"""
    NeXLUncertainties.asa(::Type{DataFrame}, rur::RecordingUpdateRule)::DataFrame

Tabulate the iteration steps in a DataFrame.
"""
function NeXLUncertainties.asa(::Type{DataFrame}, rur::RecordingUpdateRule)
    df = DataFrame(Iter = eachindex(rur.comps))
    for elm in keys(rur.comps[1])
        df[!, Symbol("Cest[$(elm.symbol)]")] = [ get(rur.comps[i], elm, 0.0) for i in eachindex(rur.comps) ]
    end
    return df
end

function NeXLMatrixCorrection.update( #
    rur::RecordingUpdateRule,
    prevcomp::Material,
    measured::Vector{KRatio},
    zafs::Dict{Element,Float64},
)::Dict{Element,Float64}
    if length(rur.comps)==0
        push!(rur.comps, Dict(elm=>prevcomp[elm] for elm in keys(prevcomp)))
    end
    res = NeXLMatrixCorrection.update(rur.base, prevcomp, measured, zafs)
    push!(rur.comps, res)
    return res
end

function NeXLMatrixCorrection.reset(rur::RecordingUpdateRule)
    resize!(rur.comps, 0)
    NeXLMatrixCorrection.reset(rur.base)
end

"""
The `ConvergenceTest` abstract type represents mechanisms to decide when the iteration has converged.

    converged(ct::ConvergenceTest, meas::Vector{KRatio}, computed::Dict{Element,Float64})::Bool
"""
abstract type ConvergenceTest end

"""
The `RMSBelowTolerance` `ConvergenceTest` ensures that the root-mean-squared difference between
measured and computed is below a threshold.
"""
struct RMSBelowTolerance <: ConvergenceTest
    tolerance::Float64
end

converged(rbt::RMSBelowTolerance, meas::Vector{KRatio}, computed::Dict{Element,Float64})::Bool =
    sum((value(nonnegk(kr)) - computed[kr.element])^2 for kr in meas) < rbt.tolerance^2

"""
The `AllBelowTolerance` `ConvergenceTest` ensures that the difference between
measured and computed is below a threshold for each k-ratio.
"""
struct AllBelowTolerance <: ConvergenceTest
    tolerance::Float64
end

converged(abt::AllBelowTolerance, meas::Vector{KRatio}, computed::Dict{Element,Float64})::Bool =
    all(abs(value(nonnegk(kr)) - computed[kr.element]) < abt.tolerance for kr in meas)

"""
The `IsApproximate` `ConvergenceTest` checks that the k-ratio differences are either below an absolute threshold
or a relative tolerance.
"""
struct IsApproximate <: ConvergenceTest
    atol::Float64
    rtol::Float64
end

converged(ia::IsApproximate, meas::Vector{KRatio}, computed::Dict{Element,Float64}) = all(
    (abs(1.0 - value(nonnegk(kr)) / computed[kr.element]) < rtol) ||
    (abs(value(nonnegk(kr)) - computed[kr.element]) < atol) for kr in meas
)
"""
Collects the information necessary to define the iteration process including the `MatrixCorrection` and
`FLuorescenceCorrection` algorithms, the iteration `UpdateRule`, the `ConvergenceTest`, and an
`UnmeasuredElementRule`.
"""
struct Iteration
    mctype::Type{<:MatrixCorrection}
    fctype::Type{<:FluorescenceCorrection}
    cctype::Type{<:CoatingCorrection}
    updater::UpdateRule
    converged::ConvergenceTest
    unmeasured::UnmeasuredElementRule

    Iteration(
        mct::Type{<:MatrixCorrection},
        fct::Type{<:FluorescenceCorrection},
        cct::Type{<:CoatingCorrection};
        updater = WegsteinUpdateRule(),
        converged = RMSBelowTolerance(0.00001),
        unmeasured = NullUnmeasuredRule(),
    ) = new(mct, fct, cct, updater, converged, unmeasured)
end

"""
`IterationResult` contains the results of the iteration process including a Label identifying the source of
the k-ratios, the resulting Material, the initial and final k-ratios, whether the iteration converged and the
number of steps.  The results can be output using `asa(DataFrame, ir::IterationResult)`.
"""
struct IterationResult
    label::Label
    comp::Material
    kratios::Vector{KRatio}
    computed::Dict{Element,Float64}
    converged::Bool
    iterations::Int
    iterate::Iteration
end

"""
The source of the k-ratio data as a Label (often a CharXRayLabel).
"""
source(ir::IterationResult)::Label = ir.label

function NeXLUncertainties.asa(::Type{DataFrame}, ir::IterationResult; withZAF::Bool = true)
    elms, mfs, ks, cks, labels = String[], Float64[], Union{Missing,Float64}[], Union{Missing,Float64}[], Label[]
    g, z, a, f = Union{Missing,Float64}[], Union{Missing,Float64}[], Union{Missing,Float64}[], Union{Missing,Float64}[]
    c, gzafc, stds, dmfs, xrays = Union{Missing,Float64}[], Union{Missing,Float64}[], String[], Float64[], String[]
    for elm in keys(ir.comp)
        push!(labels, ir.label)
        push!(elms, elm.symbol)
        rc = round(ir.comp[elm])
        push!(mfs, value(rc))
        push!(dmfs, σ(rc))
        added = false  # computed element
        for kr in ir.kratios
            if elm == kr.element
                push!(ks, value(kr.kratio))
                push!(xrays, repr(kr.xrays))
                if withZAF
                    zafs = _ZAF(ir.iterate, kr.standard, kr.stdProps, kr.xrays)
                    zafu = _ZAF(ir.iterate, ir.comp, kr.unkProps, kr.xrays)
                    push!(cks, σ(kr.kratio))
                    push!(stds, name(kr.standard))
                    push!(c, coating(zafu, zafs, kr.unkProps[:TakeOffAngle], kr.stdProps[:TakeOffAngle]))
                    push!(z, Z(zafu, zafs))
                    push!(a, A(zafu, zafs, kr.unkProps[:TakeOffAngle], kr.stdProps[:TakeOffAngle]))
                    push!(f, F(zafu, zafs, kr.unkProps[:TakeOffAngle], kr.stdProps[:TakeOffAngle]))
                    push!(g, generation(zafu, zafs))
                    push!(gzafc, g[end] * z[end] * a[end] * f[end])
                else
                    push!(cks, ir.computed[elm])
                end
                added = true
                break
            end
        end
        if !added
            push!(ks, missing), push!(cks, missing), push!(stds, missing), push!(c, missing), push!(z, missing),
            push!(a, missing), push!(f, missing), push!(g, missing), push!(gzafc, missing), push!(xrays, missing)
        end
    end
    return withZAF ?
           DataFrame(
        :Label => labels,
        :Element => elms,
        :Standard => stds,
        :Xrays => xrays,
        #:Converged => [ir.converged for elm in keys(ir.comp)],
        #:Iterations => [ir.iterations for elm in keys(ir.comp)],
        Symbol("C") => mfs,
        Symbol("ΔC") => dmfs,
        Symbol("k") => ks,
        Symbol("Δk") => cks,
        :g => g,
        :Z => z,
        :A => a,
        :F => f,
        :c => c,
        :gZAFc => gzafc,
    ) :
           DataFrame(
        :Label => labels,
        :Element => elms,
        :Converged => [ir.converged for elm in keys(ir.comp)],
        :Iterations => [ir.iterations for elm in keys(ir.comp)],
        Symbol("C") => mfs,
        Symbol("ΔC") => dmfs,
        Symbol("k") => ks,
        Symbol("Δk") => cks,
    )
    return res
end

function NeXLUncertainties.asa(::Type{DataFrame}, irs::AbstractVector{IterationResult}; mode = :MassFraction, nominal = nothing)::DataFrame
    if isnothing(nominal)
        asa(DataFrame, map(ir->ir.comp,irs), mode)
    else
        asa(DataFrame, [ ( ir.comp for ir in irs)..., nominal ], mode )
    end
end

DataFrames.describe(irs::AbstractVector{IterationResult}) =
    describe(asa(DataFrame, irs)[:, 2:end], :mean, :std, :min, :q25, :median, :q75, :max)

Base.show(io::IO, itres::IterationResult) = print(
    io,
    itres.converged ? "Converged to $(itres.comp) in $(itres.iterations) steps." :
        "Failed to converge in $(itres.iterations) iterations: Best estimate = $(itres.comp).",
)

NeXLCore.compare(itres::IterationResult, known::Material)::DataFrame = compare(itres.comp, known)

NeXLCore.compare(itress::AbstractVector{IterationResult}, known::Material)::DataFrame =
    mapreduce(itres -> compare(itres, known), append!, itress)
"""
    NeXLCore.material(itres::IterationResult)::Material
"""
NeXLCore.material(itres::IterationResult) = itres.comp
NeXLCore.material(itress::AbstractVector{IterationResult})  = mean(material.(itress))

_ZAF(iter::Iteration, mat::Material, props::Dict{Symbol,Any}, xrays::Vector{CharXRay})::MultiZAF =
    zafcorrection(iter.mctype, iter.fctype, iter.cctype, mat, xrays, props[:BeamEnergy], get(props, :Coating, missing))

"""
    computeZAFs(
        iter::Iteration,
        est::Material,
        stdZafs::Dict{KRatio,MultiZAF}
    )::Dict{Element, Float64}

Given an estimate of the composition compute the corresponding k-ratios.
"""
function computeZAFs(iter::Iteration, est::Material, stdZafs::Dict{KRatio,MultiZAF})::Dict{Element,Float64}
    zaf(kr, zafs) =
        gZAFc(_ZAF(iter, est, kr.unkProps, kr.xrays), zafs, kr.unkProps[:TakeOffAngle], kr.stdProps[:TakeOffAngle])
    return Dict(kr.element => zaf(kr, zafs) for (kr, zafs) in stdZafs)
end

"""
    estimatecoating(substrate::Material, coating::Material, kcoating::KRatio, mc::Type{<:MatrixCorrection}=XPP)::Film

Use the measured k-ratio to estimate the mass-thickness of the coating material on the specified substrate.
Return the result as a Film which may be assigned to the :Coating property for k-ratios associated with the
substrate.

Assumptions:
  * The coating is the same for all measurements of the unknown
  * The coating element k-ratio is assumed to be assigned to the brightest line
"""
function estimatecoating(substrate::Material, coating::Material, kcoating::KRatio, mc::Type{<:MatrixCorrection}=XPP)::Film
    coatingasfilm(mc, substrate, coating, #
            brightest(kcoating.xrays), kcoating.unkProps[:BeamEnergy], #
            kcoating.unkProps[:TakeOffAngle], value(kcoating.kratio))
end

"""
    quantify(
        iter::Iteration,
        label::Label,
        measured::Vector{KRatio};
        maxIter::Int = 100,
        estComp::Union{Nothing,Material}=nothing,
        coating::Union{Nothing, Material}=nothing
    )::IterationResult
    quantify(
        iter::Iteration,
        name::String,
        measured::Vector{KRatio};
        maxIter::Int = 100,
        estComp::Union{Nothing,Material}=nothing
        coating::Union{Nothing, Material}=nothing
    )::IterationResult

Perform the iteration procedurer as described in `iter` using the `measured` k-ratios to produce the best
estimate `Material` in an `IterationResult` object.  The third form makes it easier to quantify the
k-ratios from filter fit spectra.  `estComp` is an optional first estimate of the composition.  This can
be a useful optimization when quantifying many similar k-ratios (like, for example, points on a
hyper-spectrum.)

If `coating` is defined then

"""
quantify(iter::Iteration, name::String, measured::Vector{KRatio}; kwargs...) = quantify(iter, label(name), measured; kwargs...)

function quantify(
    iter::Iteration,
    label::Label,
    measured::Vector{KRatio};
    maxIter::Int = 100,
    estComp::Union{Nothing,Material}=nothing,
    coating::Union{Nothing, Pair{CharXRay, <:Material}}=nothing
)::IterationResult
    # Compute the C = k*C_std estimate
    firstEstimate(meas::Vector{KRatio}) =
        Dict{Element,Float64}(kr.element => value(nonnegk(kr)) * value(kr.standard[kr.element]) for kr in meas)
    # Compute the estimated k-ratios
    computeKs(estComp::Material, zafs::Dict{Element,Float64}, stdComps::Dict{Element,Float64}) =
        Dict{Element,Float64}(elm => estComp[elm] * zafs[elm] / stdComps[elm] for (elm, zaf) in zafs)
    function computefinal(estcomp::Material, meas::Vector{KRatio})
        final = Dict{Element,UncertainValue}(elm=>convert(UncertainValue, estcomp[elm]) for elm in keys(estcomp))
        for kr in meas
            elm = element(kr)
            unc = if value(final[elm])>0.0 && value(kr.kratio) > 0.0 && σ(kr.kratio) > 0.0
                value(final[elm])*fractional(kr.kratio)
            else
                σ(kr.kratio)
            end
            final[elm] = uv(value(final[elm]), unc)
        end
        return material(name(estcomp), final)
    end
    @assert isnothing(coating) || (get(last(coating), :Density, -1.0) > 0.0) "You must provide a positive density for the coating material."
    # Is this k-ratio due to the coating?
    iscoating(k, coatmat) =  (!isnothing(coatmat)) && (first(coatmat) in k.xrays) && (element(k) in keys(last(coatmat)))
    # k-ratios from measured elements in the unknown - Remove k-ratios for unmeasured and coating elements
    kunk = filter(measured) do kr
        !(isunmeasured(iter.unmeasured, element(kr)) || iscoating(kr, coating))
    end
    # k-ratios associated with the coating
    kcoat = filter(kr->iscoating(kr, coating), measured)
    # Compute the k-ratio difference metric
    eval(computed) = sum((value(nonnegk(kr)) - computed[kr.element])^2 for kr in kunk)
    # Compute the standard matrix correction factors
    stdZafs = Dict(kr => _ZAF(iter, kr.standard, kr.stdProps, kr.xrays) for kr in kunk)
    stdComps = Dict{Element,Float64}(kr.element => value(kr.standard[kr.element]) for kr in kunk)
    # First estimate c_unk = k*c_std
    estcomp = something(estComp, material(repr(label), compute(iter.unmeasured, firstEstimate(kunk))))
    # Compute the associated matrix corrections
    zafs = computeZAFs(iter, estcomp, stdZafs)
    bestComp, bestKrs = estcomp, computeKs(estcomp, zafs, stdComps)
    bestEval, bestIter = 1.0e300, 0
    reset(iter.updater)
    for iters in Base.OneTo(maxIter)
        if length(kcoat) >= 1
            coatings = estimatecoating(estcomp, last(coating), first(kcoat), iter.mctype)
            # Previous coatings are replaced on all the unknown's k-ratios
            foreach(k->k.unkProps[:Coating] = coatings, kunk)
        end
        # How close are the calculated k-ratios to the measured version of the k-ratios?
        estkrs = computeKs(estcomp, zafs, stdComps)
        if eval(estkrs) < bestEval
            # If no convergence report it but return closest result...
            bestComp, bestKrs, bestEval, bestIter = estcomp, estkrs, eval(estkrs), iters
            if converged(iter.converged, kunk, bestKrs)
                fc = computefinal(estcomp, kunk)
                return IterationResult(label, fc, measured, bestKrs, true, bestIter, iter)
            end
        end
        # Compute the next estimated mass fractions
        upd = update(iter.updater, estcomp, kunk, zafs)
        # Apply unmeasured element rules
        estcomp = material(repr(label), compute(iter.unmeasured, upd))
        # calculated matrix correction for estcomp
        zafs = computeZAFs(iter, estcomp, stdZafs)
    end
    @warn "$label did not converge in $(maxIter)."
    @warn "   Using best non-converged result from step $(bestIter)."
    return IterationResult(label, bestComp, measured, bestKrs, false, bestIter, iter)
end

function quantify(
    sampleName::String,
    measured::Vector{KRatio};
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
    coating::Union{Nothing, Pair{CharXRay, <:Material}}=nothing
)
    quantify(Iteration(mc, fc, cc), label(sampleName), measured, coating=coating)
end

function quantify(
    lbl::Label,
    measured::AbstractVector{KRatio};
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
    coating::Union{Nothing, Pair{CharXRay, <:Material}}=nothing
)
    quantify(Iteration(mc, fc, cc), lbl, measured, coating=coating)
end

function quantify(
    measured::AbstractVector{<:KRatios};
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = NullFluorescence,
    cc::Type{<:CoatingCorrection} = NullCoating,
    kro::KRatioOptimizer = SimpleKRatioOptimizer(1.5),
    coating::Union{Nothing, Pair{CharXRay, <:Material}}=nothing
)
    @assert all(size(measured[1])==size(krsi) for krsi in measured[2:end]) "All the KRatios need to be the same dimensions."
    # Pick the best sub-selection of `measured` to quantify
    optmeasured = brightest.(optimizeks(kro, measured))
    ThreadsX.map(CartesianIndices(optmeasured[1].kratios)) do ci
        try
            krs = KRatio[ kr[ci] for kr in optmeasured]
            quantify(Iteration(mc, fc, cc), label(ci), krs, coating=coating).comp
         catch ex
             @error ex
             NeXLCore.NULL_MATERIAL
         end
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
using Cairo
using LSystems
using DataStructures

c = CairoRGBSurface(1024, 1024);
cr = CairoContext(c);

save(cr);

set_source_rgb(cr, 1.0, 1.0, 1.0);
rectangle(cr, 0.0, 0.0, 1024.0, 1024.0);
fill(cr);

restore(cr);
save(cr);

x, y, a, d = 50., 512., 0, 15
set_line_width(cr, 0.6);
move_to(cr, x, y);

sierpinski_start = "F+XF+F+XF"
sierpinski_trans = Dict([('X', "XF-F+F-XF+F+XF-F+F-X")])

function determine_new_position(x, y, a, d)
    x -= sind(a) * d
    y -= cosd(a) * d
    return x, y
end

s = Stack(Vector{Float64}) # [x, y, a]

level = 4
drawing_instructions = @task lindenmayer(sierpinski_start, sierpinski_trans, level)

for v in drawing_instructions
    if v == 'F'
        x, y = determine_new_position(x, y, a, d)
        line_to(cr, x, y)
    elseif v == '+'
        a -= 90
    elseif v == '-'
        a += 90
    elseif v == '['
        push!(s, [x, y, a])
    elseif v == ']'
        x, y, a = pop!(s)
        move_to(cr, x, y)
    end
end

stroke(cr);
write_to_png(c, "sierpinski_$level.png");


text:
using Klimakoffer

NT = 48

mesh = Mesh(150)

co2_concentration = 315.0

model = Model(mesh, NT; co2_concentration = co2_concentration, compute_albedo=true)

disc = Discretization(mesh, model, NT)

GlobTemp = compute_equilibrium!(disc)


text:
"""
    AbstractCrystal{N}

Supertype for all types that describe crystal structures.

All subtypes of `AbstractCrystalData{N,T}` consist of a basis, space group with origin, atomic
sites, and information about the choice of cell vectors.
"""
abstract type AbstractCrystal{N}
end

"""
    CrystalStructure{N}

A crystal structure in `N` dimensions. Contains information about the lattice vectors, the space
group (and origin setting), the list of atoms and their positions, and whether the lattice vectors
are primitive and/or conventional.
"""
struct CrystalStructure{N} <: AbstractCrystal{N}
    # Lattice vectors
    basis::SVector{N, SVector{N, Float64}}
    # Space group - set to 0 if there is no information provided
    spgrp::Int
    # Origin of the space group within the cell
    origin::SVector{N, Float64}
    # Names of atoms
    atoms::Vector{String}
    # Generating set of atoms given the space group and lattice vectors
    sites::Vector{SVector{N, Float64}}
    # Whether the lattice is the conventional lattice
    conv::Bool
    # Whether the lattice is primitive
    prim::Bool
    # Both can be true if the space group has primitive centering
    # Both can be false if a non-standard setting is used (e.g. hexagonal cell for cubic crystal)

    # Inner constructor
    function CrystalStructure{N}(;
        basis::AbstractVector{<:AbstractVector{<:Real}} = _basis{N}(),
        spgrp::Integer = 0,
        origin::AbstractVector{<:Real} = zeros(SVector{3, Float64}),
        atoms::AbstractVector{<:AbstractString} = Vector{String}[],
        sites::AbstractVector{<:AbstractVector{<:Real}} = Vector{Float64}[],
        conv::Bool = true,
        prim::Bool = false,
    ) where N
        # Verify the dimensionality of the data
        basis_vl = [length(v) for v in basis]
        sites_vl = [length(v) for v in sites]
        @assert _allsame(basis_vl) "inconsistent lattice vector lengths"
        @assert _allsame(sites_vl) "inconsistent site position vector lengths"
        @assert N == first(basis_vl) "expected $N-dimensional basis vectors, " *
            "got $(first(basis_vl))-dimensional vectors"
        @assert N == first(sites_vl) "expected $N-dimensional site position vectors, " *
            "got $(first(sites_vl))-dimensional vectors"
        @assert N == length(basis) "expected $N basis vectors, got $(length(basis))"
        @assert length(sites) == length(atoms) "number of atoms and sites do not match"
        # Make sure that site and origin vector components are between 0 and 1
        # This may not be needed for CrystalStructure{N}
        # but would be required for CrystalGenerator{N}
        #=
        origin = origin - floor.(origin)
        sites = [v .- floor.(v) for v in sites]
        =#
        # Make sure the space group is valid
        @assert 0 <= spgrp <= NO_SGS[N] "$spgrp is not a valid space group number in $N dimensions"
        return new(basis, spgrp, origin, atoms, sites, conv, prim)
    end
end

# TODO: Equality checking for CrystalStructure{N} is going to be quite complex!
# Two structures may be equal even if their basis vectors and atomic positions are not the same
# because they may have different settings or different choice of basis vectors.
# This is a preliminary attempt but the details are going to take a while.
function Base.:(==)(xtal1::CrystalStructure{N}, xtal2::CrystalStructure{N}) where N
    # Number of atoms should be equal
    size(xtal1.atoms) == size(xtal2.atoms) || return false
    # Types of atoms should be checked next
    if xtal1.basis != xtal2.basis
        # TODO: check that the basis vectors can be converted
        return false
    end
    # If all checks pass, return true
    return true
end

"""
    _enumerate_unique_elements(itr)

Returns a `Vector{Tuple{<:Any, Int}}` that associates every element in `itr` with a numerical index
which counts the number of times that element has occurred counting from the first element.

This is used to count numbers of atoms in a crystal.

# Examples
```jldoctest
julia> _enumerate_unique_elements([2, 2, 3, 2, 3, 3, 1, 1, 2])
9-element Vector{Int}:
 1
 2
 1
 3
 2
 3
 1
 2
 4
```
"""
function _enumerate_unique_elements(v::AbstractVector)
    items = unique(v)
    count = zeros(size(items))
    index = Vector{Int}(undef, length(v))
    for (m, el) in enumerate(v)
        n = findfirst(isequal(el), items)
        index[m] = (count[n] += 1)
    end
    return collect(zip(v, index))
end

printbasis(io::IO, xtal::CrystalStructure{N}) where N = printbasis(io, xtal.basis)

"""
    printsites(io::IO, xtal::CrystalStructure{N}) where Names

Prints the atomic sites for a crystal structure in the following format:

```
Atomic sites:
Ti      0.333333    0.666667    0.250000
Ti      0.666667    0.333333    0.750000
```
"""
function printsites(io::IO, xtal::CrystalStructure{N}) where N
    tostr(x::Number) = lpad(@sprintf("%f", x), 12)
    println(io, "Atomic sites:")
    for n in 1:length(xtal.atoms)
        println(io, rpad(xtal.atoms[n], 4), join(tostr.(xtal.sites[n])))
    end
end

# Display the data in a CrystalStructure{N}
function Base.display(xtal::CrystalStructure{N}) where N
    println(string(typeof(xtal)) * ":")
    # Print the space group
    if xtal.spgrp == 0
        println("No space group specified")
    else
        println("Space group " * string(xtal.spgrp))
    end
    # Print the basis vectors
    println()
    printbasis(stdout, xtal)
    # Print the atomic sites:
    println()
    printsites(stdout, xtal)
end

"""
    formula(xtal::CrystalStructure{N})

Returns the formula string for a crystal.
"""
function formula(xtal::CrystalStructure{N}) where N

end

"""
    AbstractCrystalData{N,T}

Supertype for all types that describe crystal structures with datasets of type `T`.

All subtypes of `AbstractCrystalData{N,T}` must consist of a field containing crystal data and a
field which contains datasets - usually, stored in a dictionary.
"""
abstract type AbstractCrystalData{N,T} <: AbstractCrystal{N}
end

"""
    CrystalStructureWithData{N,T}

A combination of a `CrystalStructure{N}` with `N`-dimensional volumetric datasets. The datasets are
stored in a `Dict{String, DataGrid{N,T}}`.

When a `CrystalStructureWithData{N,T}` is generated with a single dataset, the dataset may be
accessed using the empty string `""`.
"""
struct CrystalStructureWithData{N,T} <: AbstractCrystalData{N,T}
    xtal::CrystalStructure{N}
    data::Dict{String, DataGrid{N,T}}
end

"""
    printdatasets(io, d::AbstractDict{<:Any, DataGrid{N,T}}) where {N,T}

Prints brief info about the datasets contained by a dictionary of datagrids.
"""
function printdatasets(io, d::AbstractDict{<:Any, DataGrid{N,T}}) where {N,T}
    if isempty(d)
        println(io, "No datasets specified")
    else
        str = "Datasets (of type "
        l = lastindex(str) - 1
        println(io, str * string(T) * "):")
        for (k,v) in d
            sz = size(v.data)
            println(io, lpad(k, l) * ":   " * join(string.(sz), '×'))
            # TODO: provide info for dataset vectors and origin
        end
    end
    return nothing
end

function printbasis(io::IO, xtaldata::CrystalStructureWithData{N,T}) where {N,T}
    printbasis(io, xtaldata.xtal.basis)
end

function printsites(io::IO, xtaldata::CrystalStructureWithData{N,T}) where {N,T}
    printsites(io, xtaldata.xtal)
end

function printdatasets(io::IO, xtaldata::CrystalStructureWithData{N,T}) where {N,T}
    printdatasets(io, xtaldata.data)
end

# Display the data in a CrystalStructureWithData{N,T}
function Base.display(xtaldata::CrystalStructureWithData{N,T}) where {N,T}
    println(string(typeof(xtaldata)) * ":")
    # Print the space group
    if xtaldata.xtal.spgrp == 0
        println("No space group specified")
    else
        println("Space group " * string(xtaldata.xtal.spgrp))
    end
    # Print the basis vectors:
    println()
    printbasis(stdout, xtaldata)
    # Print the atomic sites:
    println()
    printsites(stdout, xtaldata)
    # Print the datasets:
    println()
    printdatasets(stdout, xtaldata)
end


text:
using Revise
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

code = @dice begin
    function uniform(b::Int, w::Int) # b is the bits for uniform, w is the bitwidth
        x = Vector(undef, b)
        for i = b:-1:1
            x[i] = flip(0.5)
        end
        return add_bits(DistInt(x), w - b)
    end
    sum = uniform(2, 4)
    for i = 1:4
        sum = (sum + uniform(2, 4))
    end
    println(max_bits(sum[1]))
    sum[1]
    # prob_equals(sum, 3)
end

# BDD analysis
# bdd = compile(code)
# num_flips(bdd)
# num_nodes(bdd)
# infer(code, :bdd)

function execute(code)
    @time infer(code, :bdd)
end

execute(code)

# IR analysis
# to_dice_ir(code)
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)


text:
using Test
using MINLPLib

import JuMP

minlplib_dir = joinpath(dirname(pathof(MINLPLib)), "..")

# Perform Tests
include("loading.jl")

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
################################################################################
#
#  Type
#
################################################################################

embedding_type(::Type{AnticNumberField}) = NumFieldEmbNfAbs

embedding_type(::AnticNumberField) = embedding_type(AnticNumberField)

################################################################################
#
#  Field access
#
################################################################################

number_field(f::NumFieldEmbNfAbs) = f.K

isreal(f::NumFieldEmbNfAbs) = f.isreal

_absolute_index(f::NumFieldEmbNfAbs) = f.i

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

function complex_embeddings(K::AnticNumberField; conjugates::Bool = true)
  res = get_attribute!(K, :complex_embeddings) do
    return __complex_embeddings(K)
  end::Vector{embedding_type(K)}
  if conjugates
    return res
  else
    r, s = signature(K)
    return res[1:(r + s)]
  end
end

function __complex_embeddings(K::AnticNumberField)
  c = conjugate_data_arb(K)
  res = Vector{embedding_type(K)}(undef, degree(K))
  for i in 1:degree(K)
    res[i] = NumFieldEmbNfAbs(K, c, i)
  end
  return res
end

################################################################################
#
#  Conjugate
#
################################################################################

function conj(f::NumFieldEmbNfAbs)
  return complex_embeddings(number_field(f))[f.conjugate]
end

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

function Base.:(==)(f::NumFieldEmbNfAbs, g::NumFieldEmbNfAbs)
  return number_field(f) === number_field(g) &&
      _absolute_index(f) == _absolute_index(g)
end

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

function Base.show(io::IO, ::MIME"text/plain", f::NumFieldEmbNfAbs)
  print(io, "Embedding of\n", number_field(f), "\ncorresponding to ≈ ")
  _print_acb_neatly(io, f.r)
end

function Base.show(io::IO, f::NumFieldEmbNfAbs)
  print(io, "Embedding corresponding to ≈ ")
  _print_acb_neatly(io, f.r)
end

################################################################################
#
#  Evaluation
#
################################################################################

evaluate(x::nf_elem, f::NumFieldEmbNfAbs, p::Int) = f(x, p)

function (f::NumFieldEmbNfAbs)(x::nf_elem, abs_tol::Int = 32)
  K = parent(x)
  d = degree(K)
  r1, r2 = signature(K)
  target_tol = abs_tol
  abs_tol = Int(floor(abs_tol * 1.1))

  i = f.i

  while true
    prec_too_low = false
    c = conjugate_data_arb_roots(K, abs_tol)

    abs_tol < 0 && error("Precision overflow in evaluation of embedding")

    CC = AcbField(abs_tol, cached = false)
    RR = ArbField(abs_tol, cached = false)

    xpoly = arb_poly(parent(K.pol)(x), abs_tol)

    if i <= r1
      o = RR()
      ccall((:arb_poly_evaluate, libarb), Nothing,
            (Ref{arb}, Ref{arb_poly}, Ref{arb}, Int),
             o, xpoly, c.real_roots[i], abs_tol)

      if !isfinite(o) || !radiuslttwopower(o, -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
      else
        cc = CC(o)
        break
      end
    else
      tacb = CC()
      j = i <= r1 + r2 ? i - r1 : i - r1 - r2
      ccall((:arb_poly_evaluate_acb, libarb), Nothing,
            (Ref{acb}, Ref{arb_poly}, Ref{acb}, Int),
             tacb, xpoly, c.complex_roots[j], abs_tol)

      if !isfinite(tacb) || !radiuslttwopower(tacb, -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
      else
        if i <= r1 + r2
          cc = tacb
        else
          cc = conj(tacb)
        end
        break
      end
    end

    if prec_too_low
      continue
    end
  end

  cc = expand!(cc, -target_tol)
  return cc
end

#################################################################################
#
#  Restrict from bigger field
#
################################################################################

function restrict(e::NumFieldEmb, f::NumFieldMor{AnticNumberField, <: Any, <: Any})
  @req number_field(e) === codomain(f) "Number fields do not match"
  L = domain(f)
  emb = complex_embeddings(L)
  a = gen(L)
  # I need to find the embedding of the base_field of L
  b = e(f(a))
  cn = Bool[overlaps(b, e(a)) for e in emb]
  @assert count(cn) == 1
  i = findfirst(cn)
  return emb[i]
end

################################################################################
#
#  Easier creation of embeddings
#
################################################################################

function _find_nearest_embedding(K::AnticNumberField, x::Union{BigFloat, Float64})
  r = real_embeddings(K)
  diffs = [e(gen(K)) - x for e in r]
  t = [abs(z) for z in diffs]
  for i in 1:length(t)
    for j in (i + 1):length(t)
      if overlaps(t[i], t[j])
        possible = [ Float64(real(e.r)) for e in r]
        s = IOBuffer()
        for i in 1:length(possible)
          @printf s "%.2f" possible[i]
          if i < length(possible)
            print(s, ", ")
          end
        end
        ss = String(take!(s))
        error("""Given approximation not close enough to a root. Possible roots are:
                 $ss
              """)
      end
    end
  end
  _, i = findmin(t)
  return r[i]
end

function real_embedding(K::AnticNumberField, x::Union{BigFloat, Float64})
  return _find_nearest_embedding(K, x)
end

_is_contained_in_interval(x::arb, i::Tuple) = i[1] < x && x < i[2]

function _find_nearest_embedding(K::AnticNumberField, x::Tuple)
  r = real_embeddings(K)
  p = 32
  fls = [_is_contained_in_interval(real(i(gen(K), p)), x) for i in r]
  while count(fls) != 1
    p = 2 * p
    fls = [_is_contained_in_interval(real(i(gen(K), p)), x) for i in r]
    if count(fls) > 1
      possible = [ Float64(real(e.r)) for e in r]
      s = IOBuffer()
      for i in 1:length(possible)
        @printf s "%.2f" possible[i]
        if i < length(possible)
          print(s, ", ")
        end
      end
      ss = String(take!(s))
      error("""Interval contains more than one root. Possible roots are:
            $ss
            """)
    end
    p > 2^8 && error("Something wrong!")
  end
  i = findfirst(fls)
  @assert i !== nothing
  return r[i]
end

function real_embedding(K::AnticNumberField, x::Tuple)
  return _find_nearest_embedding(K, x)
end


text:
@testset "738.monotone-increasing-digits.jl" begin
    @test monotone_increasing_digits(10) == 9
    @test monotone_increasing_digits(332) == 299
    @test monotone_increasing_digits(1234) == 1234
    @test monotone_increasing_digits(1234321) == 1233999
    @test monotone_increasing_digits(2333332) == 2299999
end


text:
#####
##### Experimental
#####
using SIMD
import Base.Cartesian: @ntuple
const intlog2 = Core.Intrinsics.cttz_int

const Bucket{NFields, NElems} = NTuple{NFields, Vec{NElems}}

Base.@pure @inline function Base.eltype(T::Type{<:Bucket})
    Tuple{map(eltype, T.parameters)...}
end

struct AoSoA{T, Storage <: Bucket} <: AbstractVector{T}
    storage::Vector{Storage}
    function AoSoA(storage)
        S = eltype(storage)
        new{eltype(S),S}(storage)
    end
end

bucketwidth(::AoSoA{<:Any, <:Bucket{<:Any, NElems}}) where {NElems} = NElems

@generated function makebuckets!(arr, ::Val{bucketwidth}, columns::Vararg{<:Any, N}) where {bucketwidth, N}
    quote
        for k in 0:length(arr)-1
            arr[k+1] = @ntuple $N j ->Vec((@ntuple $bucketwidth i->columns[j][k * $bucketwidth + i]))
        end
        return arr
    end
end

function makeAoSoA(::Val{bucketwidth}, columns...) where {bucketwidth}
    len = length(first(columns))
    exp = intlog2(UInt(bucketwidth))
    @assert 2^exp == bucketwidth
    @assert all(x->length(x) == len, columns)
    @assert len % bucketwidth == 0
    BucketT = Tuple{map(T->Vec{bucketwidth, T}, eltype.(columns))...}
    AoSoA(makebuckets!(Vector{BucketT}(undef, div(len, bucketwidth)), Val{bucketwidth}(), columns...))
end

Base.size(a::AoSoA) = (length(a.storage) * bucketwidth(a),)
@inline function idx2bucketloc(a, i)
    w = UInt(bucketwidth(a))
    nshift = intlog2(w)
    ci = i - 1
    bucketidx = ci >> nshift
    j = ci & (w - 1)
    return bucketidx+1, j+1
end

Base.@propagate_inbounds function Base.getindex(a::AoSoA, i::Integer)
    bucketidx, j = idx2bucketloc(a, i)
    map(arr->arr[j], a.storage[bucketidx])
end

Base.@propagate_inbounds function Base.setindex!(a::AoSoA, vs, i::Integer)
    bucketidx, j = idx2bucketloc(a, i)
    newbucket = map((arr, v)->Base.setindex(arr, v, j), a.storage[bucketidx], vs)
    a.storage[bucketidx] = newbucket
    vs
end

#=

using SIMD

function sumfields(arr)
    T = Vec{4, Float64}
    T2 = Vec{bucketwidth(arr), Float64}
    init = zero(T)
    @inbounds for i in eachindex(arr.storage)
        v1 = T2(arr.storage[i][1])
        v2 = T2(arr.storage[i][2])
        v3 = T2(arr.storage[i][3])
        v4 = T2(arr.storage[i][4])
        init += T((sum(v1), sum(v2), sum(v3), sum(v4)))
    end
    #mapreduce(sumbucket, (x,y)->map(+, x, y), arr.storage, init=map(zero, first(arr)))
    return init
end
=#


text:
module  IOtxt

# types
StringVector = Vector{String}

export
    getFilesInDir,
    readEEG

# read EEG data from a .txt file in LORETA format and put it in a matrix
# of dimension txn, where n=#electrodes and t=#samples.
# if optional keyword argument `msg` is not empty, print `msg` on exit
function readEEG(filename::String; msg::String="")
    S=readlines(filename) # read the lines of the file as a vector of strings
    t=length(S) # number of samples
    n=length(split(S[1])) # get the number of electrodes
    X=Matrix{Float64}(undef, t, n) # declare the X Matrix
    for j=1:t
        x=split(S[j]) # this get the n potentials from a string
        for i=1:n
            X[j, i]=parse(Float64, x[i])
        end
    end
    if !isempty(msg) println(msg) end
    return X
end

# Read several EEG data from .txt files in LORETA format given in `filenames`
# ( a Vector of strings) and put them in a vector of matrices object
# print "read file "*[filenumber]*": "*[filename] after each file has been read.
readEEG(filenames::StringVector, skip::Vector{Int}=[]) =
        [readEEG(filenames[f]; msg="read file $f: "*basename(filenames[f])) for f in eachindex(filenames) if f ∉ skip]

# get the complete path of all files in `dir` as a vector of strings
# `ext` is an optional tuple of file extensions given strings.
# If it is provided, only files with those extensions will be included
# in the returned vector.
# the extensions must be entered in lowercase.
## Examples
# S=getFilesInDir(@__DIR__)
# S=getFilesInDir(@__DIR__; ext=(".txt", ))
# S=getFilesInDir(@__DIR__; ext=(".txt", ".jl"))
function getFilesInDir(dir::String; ext::Tuple=())
    if !isdir(dir) @error "Function `getFilesInDir`: input directory is incorrect!"
    else
        S=[]
        for (root, dirs, files) in walkdir(dir)
            if root==dir
                for file in files
                    if ext==() || ( lowercase(string(splitext(file)[2])) ∈ ext )
                       push!(S, joinpath(root, file)) # complete path and file name
                    end
                end
            end
        end
        isempty(S) && @warn "Function `getFilesInDir`: input directory does not contain any files"
        return StringVector(S)
    end
end

end # module

# Example
# dir="C:\\Users\\congedom\\Documents\\My Data\\EEG data\\NTE 84 Norms"

# Gat all file names with complete path
# S=getFilesInDir(dir)

# Gat all file names with complete path with extension ".txt"
# S=getFilesInDir(@__DIR__; ext=(".txt", ))

# read one file of NTE database and put it in a Matrix object
# X=readEEG(S[1])

# read all files of NTE database and put them in a vector of matrix
# 𝐗=readEEG(S)


text:
##
#
#  File: solvebasis.jl
#
#  Purpose :   To demonstrate the usage of
#              MSK_solvewithbasis on the problem:
#
#              maximize  x0 + x1
#              st.
#                      x0 + 2.0 x1 <= 2
#                      x0  +    x1 <= 6
#                      x0 >= 0, x1>= 0
#
#               The problem has the slack variables
#               xc0, xc1 on the constraints
#               and the variabels x0 and x1.
#
#               maximize  x0 + x1
#               st.
#                  x0 + 2.0 x1 -xc1       = 2
#                  x0  +    x1       -xc2 = 6
#                  x0 >= 0, x1>= 0,
#                  xc1 <=  0 , xc2 <= 0
##

using Mosek

maketask() do task

    putobjname(task,"solvebasis")

    putstreamfunc(task,MSK_STREAM_LOG,msg -> print(msg))

    numcon = 2
    numvar = 2

    c = Float64[1.0, 1.0]
    ptrb = Int64[1,3]
    ptre = Int64[3,4]
    asub = Int32[1,2,
                 1,2]
    aval = Float64[1.0, 1.0,
                   2.0, 1.0]
    bkc = Boundkey[MSK_BK_UP
                   MSK_BK_UP]

    blc = Float64[-Inf
                  -Inf]
    buc = Float64[2.0
                  6.0]

    bkx = Boundkey[MSK_BK_LO
                   MSK_BK_LO]
    blx = Float64[0.0
                  0.0]

    bux = Float64[Inf
                  Inf]
    w1 = Float64[2.0, 6.0]
    w2 = Float64[1.0, 0.0]

    inputdata(task,
              Int32(numcon), Int32(numvar),
              c,
              0.0,
              ptrb,
              ptre,
              asub,
              aval,
              bkc,
              blc,
              buc,
              bkx,
              blx,
              bux)

    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

    r = optimize(task,"mosek://solve.mosek.com:30080")
    if r != MSK_RES_OK
        println("Mosek warning: $r")
    end

    basis = initbasissolve(task)

    #List basis variables corresponding to columns of B
    varsub = Int32[1, 2]

    for i in 1:numcon
        if basis[varsub[i]] < numcon
            println("Basis variable no $i is xc$(basis[i])")
        else
            println("Basis variable no $i is x$(basis[i]-numcon)")

            # solve Bx = w1
            # varsub contains index of non-zeros in b.
            #  On return b contains the solution x and
            # varsub the index of the non-zeros in x.
            nz = 2

            nz = solvewithbasis(task,Int32(0), nz, varsub, w1)
            println("nz = $nz")
            println("Solution to Bx = $w1")

            for i in 1:nz
                if basis[varsub[i]] < numcon
                    println("xc $(basis[varsub[i]]) = $(w1[varsub[i]])")
                else
                    println("x$(basis[varsub[i]] - numcon) = $(w1[varsub[i]])")
                end
            end

            # Solve B^Tx = w2
            nz = 1
            varsub[1] = 1

            nz = solvewithbasis(task,1,nz,varsub,w2)
            println(nz)

            println("Solution to B^Tx = $(w2)")

            for i in 1:nz
                if basis[varsub[i]] < numcon
                    print("xc$(basis[varsub[i]]) = $(w2[varsub[i]])")
                else
                    print("x$(basis[varsub[i]] - numcon) = $(w2[varsub[i]])")
                end
            end
        end
    end

end


text:
using LightGraphs, Distributions, StatsBase, Random, Plots, LaTeXStrings;pyplot()
Random.seed!(0)

function createNetwork(edges)
    network = Graph(maximum(maximum.(edges)))
    for e in edges
        add_edge!(network, e[1], e[2])
    end
    network
end

function uniformRandomEdge(network)
    outDegrees = length.(network.fadjlist)
    randI = sample(1:length(outDegrees),Weights(outDegrees))
    randJ = rand(network.fadjlist[randI])
    randI, randJ
end

function networkLife(network,source,dest,lambda)
    failureNetwork = copy(network)
    t = 0
    while has_path(failureNetwork, source, dest)
        t += rand(Exponential(1/(failureNetwork.ne*lambda)))
        i, j = uniformRandomEdge(failureNetwork)
        rem_edge!(failureNetwork, i, j)
    end
    t
end

lambda1, lambda2 = 0.5, 1.0
roads = [(1,2), (1,3), (2,4), (2,5), (2,3), (3,4), (3,5), (4,5), (4,6), (5,6)]
source, dest = 1, 6
network = createNetwork(roads)
N = 10^6

failTimes1 = [ networkLife(network,source,dest,lambda1) for _ in 1:N ]
failTimes2 = [ networkLife(network,source,dest,lambda2) for _ in 1:N ]

println("Edge Failure Rate = $(lambda1): Mean failure time = ",
	mean(failTimes1), " days.")
println("Edge Failure Rate = $(lambda2): Mean failure time = ",
	mean(failTimes2), " days.")

stephist(failTimes1, bins=200, c=:blue, normed=true, label=L"\lambda=0.5")
stephist!(failTimes2, bins=200, c=:red, normed=true, label=L"\lambda=1.0",
    xlims=(0,5), ylims=(0,1.1), xlabel="Network Life Time", ylabel = "Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# Traversing -------------------------

function traverse(tn::TreeNode, f::Function, args...)

    if !isempty(args)
        f(tn, args...)
    else
        f(tn)
    end

    if !isleaf(tn)
        for (key, chnode) in tn.children
            traverse(chnode, f, args...)
        end
    end
end

function traverse(tn::Nullable{TreeNode}, args...)
    isnull(tn) ? nothing : traverse(get(tn), args...)
end