text:
struct FFTData <: JAXTAMData
    mission        :: Mission
    instrument     :: Symbol
    obsid          :: String
    e_range        :: Tuple{Float64,Float64}
    bin_time       :: Real
    bin_size       :: Int
    bin_count      :: Int
    gti_index      :: Int
    gti_start_time :: Real
    group          :: Int
    src_ctrate     :: Array
    bkg_ctrate     :: Array
    power          :: Array
    avg_power      :: Array
    freq           :: Array
end

function _FFTData(gti, freq, power, fspec_bin_size, bin_count, src_ctrate)

    return FFTData(gti.mission, gti.instrument, gti.obsid, gti.e_range, gti.bin_time, fspec_bin_size, bin_count,
        gti.gti_index, gti.gti_start_time, gti.group, src_ctrate, zeros(size(src_ctrate)), power, mean(power, dims=2)[:], freq)
end

function _fft(counts::Array, times::StepRangeLen, bin_time::Real, fspec_bin_size; leahy=true)
    FFTW.set_num_threads(4)

    spec_no = Int(floor(length(counts)/fspec_bin_size))

    counts = counts[1:Int(spec_no*fspec_bin_size)]
    counts = reshape(counts, fspec_bin_size, spec_no)

    src_ctrate = mean(counts, dims=1)./bin_time

    power = abs.(rfft(counts))
    freq = Array(rfftfreq(fspec_bin_size, 1/bin_time))

    if leahy
        power = (2 .*(power.^2)) ./ sum(counts)
    end

    return freq, power, spec_no, src_ctrate
end

function _best_gti_pow2(gtis::Dict{Int64,JAXTAM.GTIData})
    gti_pow2s = [prevpow(2, length(gti[2].counts)) for gti in gtis]
    min_pow2  = Int(median(gti_pow2s))
    exclude   = gti_pow2s .< min_pow2

    instrument = gtis[collect(keys(gtis))[1]].instrument

    @info "`$instrument` median `pow2` length is $min_pow2, excluded $(count(exclude)) gtis,
        $(length(gti_pow2s) - count(exclude)) remain"

    good_gtis = [gti for gti in gtis][.!exclude]

    return Dict(good_gtis), min_pow2
end

function gtis_pow2(mission_name::Symbol, gtis::Dict{Symbol,Dict{Int64,JAXTAM.GTIData}})
    instruments = config(mission_name).instruments

    gtis_pow2 = Dict{Symbol,Dict{Int64,JAXTAM.GTIData}}()

    for instrument in instruments
        gtis_pow2[Symbol(instrument)] = _best_gti_pow2(gtis[Symbol(instrument)])
    end

    return gtis_pow2
end

function _fspec_save(fspec_data::Dict{Int64,JAXTAM.FFTData}, fspec_dir::String)
    mkpath(fspec_dir)
    log_entries = Dict{Union{Symbol,Int},String}()

    fspec_indecies = collect(keys(fspec_data))
    fspec_starts   = [t.gti_start_time for t in values(fspec_data)]
    fspec_b_counts = [t.bin_count for t in values(fspec_data)]
    fspec_groups   = [o.group for o in values(fspec_data)]
    fspec_example  = fspec_data[fspec_indecies[1]]

    fspec_basename = string("$(fspec_example.instrument)_lc_$(fspec_example.bin_time)_fspec")
    fspec_meta     = DataFrame(mission=_mission_name(fspec_example.mission),
    instrument=String(fspec_example.instrument), obsid=fspec_example.obsid,
    bin_time=fspec_example.bin_time, bin_size=fspec_example.bin_size, groups=fspec_groups,
    indecies=fspec_indecies, starts=fspec_starts, b_counts=fspec_b_counts)

    Feather.write(joinpath(fspec_dir, "$(fspec_basename)_meta.feather"), fspec_meta)

    for index in fspec_indecies
        fspec_savepath  = joinpath(fspec_dir, "$(fspec_basename)_$(index).feather")

        fspec_data_df = DataFrame(fspec_data[index].power)
        fspec_data_df[:freq]   = fspec_data[index].freq
        fspec_data_df[:avg_power] = fspec_data[index].avg_power

        Feather.write(fspec_savepath, fspec_data_df)
    end
    @warn "Logging not implemented yet as saving is off by default" "Contact dev to request fix"
    # TODO: Fix this... eventually
end

function _fspec_load(fspec_dir::String, instrument::Symbol, e_range::Tuple{Float64,Float64}, bin_time::Number, fspec_bin::Number)
    bin_time = float(bin_time)

    fspec_basename  = string("$(instrument)_lc_$(bin_time)_fspec")
    fspec_meta_path = joinpath(fspec_dir, "$(fspec_basename)_meta.feather")

    fspec_meta = Feather.read(fspec_meta_path)

    fspec_data = Dict{Int64,JAXTAM.FFTData}()

    for fspec_row in eachrow(fspec_meta)
        fspec_mission = _mission_symbol_to_type(Symbol(fspec_row[:mission]))
        fspec_inst    = Symbol(fspec_row[1, :instrument])
        fspec_obsid   = fspec_row[1, :obsid]
        fspec_bin_t   = fspec_row[1, :bin_time]
        fspec_bin_sze = fspec_row[1, :bin_size]
        fspec_idx     = fspec_row[1, :indecies]
        fspec_starts  = fspec_row[1, :starts]
        fspec_b_count = fspec_row[1, :b_counts]
        fspec_group   = fspec_row[1, :groups]
        fspec_file    = Feather.read(joinpath(fspec_dir, "$(fspec_basename)_$(fspec_idx).feather"))
        fspec_freq   = Array(fspec_file[:freq]); delete!(fspec_file, :freq)
        fspec_ampavg  = Array(fspec_file[:avg_power]); delete!(fspec_file, :avg_power)
        if fspec_idx != -1
            fspec_power = convert(Array, fspec_file[:, names(fspec_file)])
        else # If avg_power is loaded, fspec_power = [], throws Argument Error on convert
            fspec_power = []
        end

        fspec_data[fspec_idx] = FFTData(fspec_mission, fspec_inst, fspec_obsid, e_range, fspec_bin_t,
            fspec_bin_sze, fspec_b_count, fspec_idx, fspec_starts, fspec_group, fspec_power, fspec_ampavg,
            fspec_freq)
    end

    return fspec_data
end

function _fspec(gtis::Dict{Int64,JAXTAM.GTIData}, fspec_bin::Real; pow2=true, fspec_bin_type=:time)
    powspec = Dict{Int64,FFTData}()

    first_gti = gtis[collect(keys(gtis))[1]]

    if fspec_bin_type == :time
        fspec_bin_size = Int(fspec_bin/first_gti.bin_time)
    elseif fspec_bin_type == :length
        fspec_bin_size = Int(fspec_bin)
    end

    @info "                       -> fspec_bin_size = $fspec_bin_size"

    if !ispow2(fspec_bin_size)
        @warn "`fspec_bin_size` not `pow2`: $fspec_bin_size, `fspec_bin` = $fspec_bin ($(1/fspec_bin))"
    end

    for gti in values(gtis)
        if length(gti.counts) < fspec_bin_size
            @info "                       -> skipped gti $(lpad(gti.gti_index, 3)) (under by: $(length(gti.counts) - fspec_bin_size))"
            continue
        end

        freq, power, bin_count, src_ctrate = _fft(gti.counts, gti.times, gti.bin_time, fspec_bin_size)

        powspec[gti.gti_index] = _FFTData(gti, freq, power, fspec_bin_size, bin_count, src_ctrate)
    end

    return powspec
end

function _scrunch_sections(instrument_data::Dict{Int64,JAXTAM.FFTData}; append_mean=true)
    gtis = keys(instrument_data)
    gti_example = instrument_data[collect(gtis)[1]]

    joined_powspecs = Dict{Symbol,Array}()
    joined_together = Array{Float64,2}(undef, size(gti_example.power, 1), 0)

    bin_count_sum = 0

    mean_src_ctrate = Array{Float64,1}(undef, 1)
    mean_bkc_ctrate = Array{Float64,1}(undef, 1)
    for gti in values(instrument_data)
        if gti.gti_index <= 0 # Don't try and scrunch special, non-gti values
            continue
        end

        bin_count_sum += gti.bin_count
        joined_together = hcat(joined_together, gti.power)
        mean_src_ctrate = hcat(mean_src_ctrate, gti.src_ctrate)
        mean_bkc_ctrate = hcat(mean_bkc_ctrate, gti.bkg_ctrate)
    end

    mean_src_ctrate = [mean(mean_src_ctrate[2:end])]
    mean_bkc_ctrate = [mean(mean_bkc_ctrate[2:end])]

    if append_mean
        instrument_data[-1] = FFTData(
            gti_example.mission, gti_example.instrument, gti_example.obsid, gti_example.e_range,
            gti_example.bin_time, gti_example.bin_size, bin_count_sum, -1, -1, -1,
            mean_src_ctrate, mean_bkc_ctrate,
            [], mean(joined_together, dims=2)[:], gti_example.freq
        )
    end

    return instrument_data
end

function fspec(mission_name::Symbol, gtis::Dict{Symbol,Dict{Int64,JAXTAM.GTIData}}, fspec_bin::Real; pow2=true, fspec_bin_type=:time, scrunched=true)
    instruments = config(mission_name).instruments

    instrument_fspecs = Dict{Symbol,Dict{Int,JAXTAM.FFTData}}()

    for instrument in instruments
        instrument_fspecs[Symbol(instrument)] = _fspec(gtis[Symbol(instrument)], fspec_bin, pow2=pow2, fspec_bin_type=fspec_bin_type)

        if scrunched
            instrument_fspecs = _scrunch_sections(instrument_fspecs[Symbol(instrument)])
        end
    end

    return instrument_fspecs
end

function fspec(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame},
        bin_time::Real, fspec_bin::Real; overwrite=true, overwrite_gtis=false, save_fspec=false,
        pow2=true, fspec_bin_type=:time, scrunched=true, e_range::Tuple{Float64,Float64}=_mission_good_e_range(mission),
        gtis_data::Dict{Symbol,Dict{Int64,JAXTAM.GTIData}}=Dict{Symbol,Dict{Int64,JAXTAM.GTIData}}()
    )
    obsid       = obs_row[:obsid]
    instruments = _mission_instruments(mission)
    fspec_files = _log_query(mission, obs_row, "data", :fspec, e_range)

    fspec_path  = abspath(
        _obs_path_local(mission, obs_row; kind=:jaxtam),
        "JAXTAM",
        _log_query_path(; category=:data, kind=:fspec, e_range=e_range, bin_time=bin_time)
    )

    fspecs = Dict{Symbol,Dict{Int64,JAXTAM.FFTData}}()
    for instrument in instruments
        if ismissing(fspec_files) || !haskey(fspec_files, instrument) || overwrite
            @info "Missing fspec files for $instrument"

            if !haskey(gtis_data, instrument)
                gtis_data = JAXTAM.gtis(mission, obs_row, bin_time, overwrite=overwrite_gtis, e_range=e_range)
            end

            @assert _recursive_first(gtis_data).e_range == e_range "gti energy range does not match input: $(_recursive_first(gtis_data).e_range) != $e_range"

            @info "Generating fspec for $instrument"
            fspec_data = _fspec(gtis_data[instrument], fspec_bin; pow2=pow2, fspec_bin_type=fspec_bin_type)

            if length(fspec_data) == 0
                throw(JAXTAMError("No GITs left in fspec", :fspec))
            end

            if scrunched
                @info "                       -> scrunching gtis"
                fspec_data = _scrunch_sections(fspec_data)
            end

            if save_fspec
                @info "                       -> saving $instrument -> $fspec_path"
                _fspec_save(fspec_data, fspec_path)
            else
                @info "                       -> not saving fspec file, `save_fspec` is false"
            end

            fspecs[instrument] = fspec_data
        else
            @info "Loading $instrument fspecs"
            fspecs[instrument] = _fspec_load(fspec_path, instrument, bin_time, fspec_bin)
        end
    end

    return fspecs
end

function fspec(mission::Mission, obsid::String, bin_time::Number, fspec_bin::Real;
    overwrite_fs=false, overwrite_gtis=false, pow2=true, fspec_bin_type=:time, scrunched=true, e_range=_mission_good_e_range(mission))

    obs_row = master_query(mission, :obsid, obsid)

    fs = fspec(mission, obs_row, bin_time, fspec_bin; overwrite=overwrite_fs, overwrite_gtis=overwrite_gtis,
        pow2=pow2, fspec_bin_type=fspec_bin_type, scrunched=scrunched, e_range=e_range)

    return fs
end

function _fspec_rebin(power, freq, bin_count::Int, bin_size, bin_time, rebin=(:log10, 0.01))
    rebin_type     = rebin[1]
    rebin_factor   = rebin[2]

    freq_intervals = missing
    if ismissing(bin_time)
        freq_intervals = mean(diff(freq))
    else
        freq_intervals = 1/(bin_size*bin_time) #freq[3] - freq[2]
    end

    # Frequency rebinning (e.g. binning to 1Hz, 2Hz, etc... intervals) is
    # just linear rebinning, with the factor determined by the desired frequency
    binary_thresh  = missing
    if rebin_type == :freq
        rebin_type   = :linear
        rebin_factor = round(Int, bin_size*bin_time*rebin_factor)
    elseif rebin_type == :freq_binary
        rebin_factor  = round(Int, bin_size*bin_time*rebin_factor)
        binary_thresh = rebin[3]
    end

    if rebin_type == :log10
        freq_max = freq[end]

        scale_start  = log10(rebin_factor)

        scale = exp10.(range(scale_start, step=rebin_factor, stop=log10(freq_max)))
        scale = [0; scale[1:end-1]]
        scale = scale./freq_intervals
        scale = ceil.(Int, scale[:, 1])
        scale[1, 1] = 1
        scale = unique(scale)

        final_scale_value = ceil(Int, exp10(log10(scale[end])+rebin_factor))
        final_scale_value > length(power) ? final_scale_value=length(power) : ""
        scale = [scale [scale[2:end]; final_scale_value].+1]

        errors = power
        errors = [mean(errors[scale[i, 1]:scale[i, 2]]) for i = 1:size(scale, 1)]
        errors = errors ./ sqrt.(diff(scale, dims=2).*bin_count)

        rebinned_fspec = [mean(power[scale[i, 1]:scale[i, 2]]) for i = 1:size(scale, 1)]

        freq_scale = freq[round.(Int, (scale[:, 2]+scale[:, 1])./2)]

        rebinned_fspec[rebinned_fspec.==0] .= NaN
    elseif rebin_type == :linear
        if rebin_factor < 1 || typeof(rebin_factor) != Int
            error("Rebin factor must be a >1 integer for linear rebinning")
        end

        in_rebin = floor(Int, size(power, 1)/rebin_factor).*rebin_factor

        errors = power ./ sqrt(bin_count*rebin_factor)
        errors = mean(reshape(errors[1:in_rebin], rebin_factor, :), dims=1)'

        rebinned_fspec = mean(reshape(power[1:in_rebin], rebin_factor, :), dims=1)'
        rebinned_fspec = rebinned_fspec[:]

        freq_scale = mean(reshape(freq[1:in_rebin], rebin_factor, :), dims=1)'
        freq_scale = freq_scale[:]

        rebinned_fspec[freq_scale .== 0] .= NaN
    elseif rebin_type == :freq_binary
        if rebin_factor < 1 || typeof(rebin_factor) != Int
            error("Rebin factor must be a >1 integer for linear rebinning")
        end

        in_rebin = floor(Int, size(power, 1)/rebin_factor).*rebin_factor

        errors = power ./ sqrt(bin_count*rebin_factor)
        errors = zeros(in_rebin) # Errors don't really have a meaning for this

        binned_power    = reshape(power[1:in_rebin], rebin_factor, :)
        rebinned_fspec = zeros(1, size(binned_power, 2))
        rebinned_fspec[any(binned_power .>= binary_thresh, dims=1)] .= binary_thresh
        rebinned_fspec = rebinned_fspec'

        freq_scale = mean(reshape(freq[1:in_rebin], rebin_factor, :), dims=1)'
    else
        error("Invalid rebin type $rebin_type, must be :log10, :linear, :freq, or :freq_binary")
    end

    return freq_scale[:], rebinned_fspec, errors
end

function fspec_rebin(fs::JAXTAM.FFTData; rebin=(:log10, 0.01))
    rebinned_data = _fspec_rebin(fs.avg_power, fs.freq, fs.bin_count, fs.bin_size, fs.bin_time, rebin)

    return rebinned_data
end

function fspec_rebin_sgram(fs::Dict{Int64,JAXTAM.FFTData}; rebin=(:log10, 0.01))
    fs = sort(fs) # Ensure po=wer spectra is in GTI order
    fs_groups       = unique([f[2].group for f in fs if f[1]>0])
    fs_group_bounds = cumsum([sum([f[2].bin_count for f in fs if f[2].group==g]) for g in fs_groups])
    #fs_group_bounds = cumsum([f[2].bin_count for f in fs if f[1]>0])

    fs_freq   = fs[-1].freq
    fs_power  = hcat([f[2].power for f in fs if f[1] > 0]...) # Skip special < 0 keys
    fs_power[1, :] .= NaN # NaN 0 Hz

    bin_size = fs[-1].bin_size; bin_time = fs[-1].bin_time
    fs_rebinned_power  = hcat([_fspec_rebin(fs_power[:, i], fs_freq, 1, bin_size, bin_time, rebin)[2] for i in 1:size(fs_power, 2)]...)
    fs_rebinned_freq = fspec_rebin(fs[-1], rebin=rebin)[1]

    return fs_rebinned_freq, fs_rebinned_power, fs_group_bounds, fs_groups
end

function fspec_pulses(fs::Dict{Int64,JAXTAM.FFTData};
        freq_bin=1, power_limits=[10, 20, 100])

    rebinned_data = Dict{Int,Tuple}()
    for p in power_limits
        rebinned_data[p] = fspec_rebin_sgram(fs;
            rebin=(:freq_binary, freq_bin, p))
    end

    return rebinned_data
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 61078404-224f-44f1-a75a-9b3537ecc95c
# we write a function that decomposes a square matrix A as A = LU
# where L is lower triangular and U is upper triangular
function LUdecomp(A)

    # we get the size of A
    n, m = size(A)

    # we check if A is indeed a square matrix
    if n != m
        error("Input must be a square matrix")
    end

    # we initialize the output matrices
    L = zeros(n, n)
    U = zeros(n, n)

    # we make a copy of A in order not to modify the original
    B = copy(A)

    for k = 1:n

        # we set the entries of the kth row of U
        for j = k:n
            U[k,j] = B[k,j]
        end

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("** A^(k-1)_{k,k} == 0 in LU decomp")
        end

        # we set the entries of the kth column of L
        for i = k:n
           L[i,k] = B[i,k] / U[k,k]
        end

        # we modify A for the next iteration
        for i = k:n
            for j = k:n
                B[i,j] -= L[i,k] * U[k,j]
            end
        end
    end

    # we return L and U
    return (L, U)
end

# ╔═╡ 4996c2fe-068b-11ec-221d-5fdd71b0f050
md"Simple matrix for which there is no LUdecomp:"

# ╔═╡ 41b2c106-e317-4d65-bb61-94cecb4683a2
A = [0.0 2.0; 3.0 1.0];

# ╔═╡ 75d691e3-34e6-4506-b9b6-559c2412d791
md"This will cause an error:"

# ╔═╡ ffa8307d-751a-42d7-8a5d-ce8a3e92ca0d
LUdecomp(A);

# ╔═╡ a2c516ea-c23f-4f25-8e39-b5213e6fb94e
# we write a function that returns the row index of the (absolutely) largest element in column k of the matrix A
function findLargestInCol(A, k)

	# we get column k
	col = A[:,k]

	# we take elementwise absolute value
	absCol = broadcast(abs, col)

	# we get the index of the largest element
	index = argmax(absCol)

	# we return index
	return index
end

# ╔═╡ ca343c99-8fa8-4489-8c57-3e8afc1fbe71
# we write a more efficient function that swaps rows u and v of some matrix A
function swapRowsV2!(A, u, v)

	# we directly swap rows u and v of A without copying
	temp = A[u, :]
	A[u, :] = A[v, :]
	A[v, :] = temp
end

# ╔═╡ e10c36a3-d0c0-4f63-ac43-f1f04eff1e19
# we write another function that composes A into A = P^{-1}LU where P is stored as a vector and swapRows is more efficient
function PALUdecompV2(A)

    # we get the size of A
    m, n = size(A)

    # we check that A is indeed a square matrix
    if m != n
        error("Input must be a square matrix")
    end

    # we initialize L and U to zero, P to {1,2, ..., n}
    L = zeros(n, n)
    U = zeros(n, n)
    p = collect(1:n)

    # we make a copy of A to work with
    B = copy(A)

    for k = 1:n

        # we find the pivot
        maxk = findLargestInCol(B, k)
        if maxk != k
            swapRowsV2!(B, maxk, k)
            swapRowsV2!(L, maxk, k)

            # we swap two elements of P instead of two rows
            temp = p[maxk]
            p[maxk] = p[k]
            p[k] = temp
        end

        # we proceed with the normal LU decomposition algorithm
        for j = k:n
            U[k,j] = B[k,j]
        end

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("** A^(k-1)_{k,k} == 0 in PALU decomp")
        end

        for i = k:n
           L[i,k] = B[i,k] / U[k,k]
        end

        for i = k:n
            for j = k:n
                B[i,j] -= L[i,k] * U[k,j];
            end
        end
    end

	P = zeros(n,n)
	for i = 1:n
		P[i,p[i]] = 1
	end

    # we return P, L and U
    return P, L, U
end

# ╔═╡ a0f23d95-18a5-427f-b52c-ebdb44df8fa3
PA, LA, UA = PALUdecompV2(A);

# ╔═╡ 23babfbd-d571-454e-8598-3fc9a9044baa
LA*UA == PA*A

# ╔═╡ a5570224-5282-48a3-8f5a-7f5f7a9b741c
md"""
Slightly less trivial example where the error appears in step $k=2$. The problem arises because $B_{22} = B_{21} B_{12} / B_{11}$ so in the usual LU decomposition $A^{1}_{22}$ is zero:
"""

# ╔═╡ 1dd12bff-9da3-4f26-ab3c-23a0e2671aa1
B = [1.0 2.0 3.0; 2.0 4.0 2.0; 1.0 3.0 2.0];

# ╔═╡ 5a34190b-5c77-4732-9134-c6bd432c053c
LUdecomp(B)

# ╔═╡ 751479d0-269c-4be3-893c-1317f9707b5f
PB, LB, UB = PALUdecompV2(B);

# ╔═╡ f5229f99-e78b-48e0-853e-ea1b50aa2633
LB*UB == PB*B

# ╔═╡ 3304fcdb-923f-4b21-932d-cdbfa1c4b243
md"Check that LUdecomp of $PB$ gives the same answer:"

# ╔═╡ b76c8925-323f-423d-98b0-a42bb7ef4f60
LUdecomp(PB*B) == (LB, UB)

# ╔═╡ 344175ec-c622-4967-8f4d-4c1507886f7a
md"Bigger matrix..."

# ╔═╡ 7bfae705-78ff-4f4b-8f85-4f623488d632
C = rand(800, 800);

# ╔═╡ 5c497b8f-d2a3-4609-bd5d-891fe21a4714
PC, LC, UC = PALUdecompV2(C);

# ╔═╡ 5c8e5b61-137e-401b-96f5-df3c9ffaa9dd
md"Estimate 'error' as sum of squares of elements of $LU-PA$:"

# ╔═╡ f76c12c0-3a8e-4ef1-b8dc-a5b2eaf8a761
sum((LC*UC-PC*C).^2)

# ╔═╡ Cell order:
# ╠═61078404-224f-44f1-a75a-9b3537ecc95c
# ╟─4996c2fe-068b-11ec-221d-5fdd71b0f050
# ╠═41b2c106-e317-4d65-bb61-94cecb4683a2
# ╟─75d691e3-34e6-4506-b9b6-559c2412d791
# ╠═ffa8307d-751a-42d7-8a5d-ce8a3e92ca0d
# ╠═a2c516ea-c23f-4f25-8e39-b5213e6fb94e
# ╠═ca343c99-8fa8-4489-8c57-3e8afc1fbe71
# ╠═e10c36a3-d0c0-4f63-ac43-f1f04eff1e19
# ╠═a0f23d95-18a5-427f-b52c-ebdb44df8fa3
# ╠═23babfbd-d571-454e-8598-3fc9a9044baa
# ╟─a5570224-5282-48a3-8f5a-7f5f7a9b741c
# ╠═1dd12bff-9da3-4f26-ab3c-23a0e2671aa1
# ╠═5a34190b-5c77-4732-9134-c6bd432c053c
# ╠═f5229f99-e78b-48e0-853e-ea1b50aa2633
# ╠═751479d0-269c-4be3-893c-1317f9707b5f
# ╟─3304fcdb-923f-4b21-932d-cdbfa1c4b243
# ╠═b76c8925-323f-423d-98b0-a42bb7ef4f60
# ╟─344175ec-c622-4967-8f4d-4c1507886f7a
# ╠═7bfae705-78ff-4f4b-8f85-4f623488d632
# ╠═5c497b8f-d2a3-4609-bd5d-891fe21a4714
# ╟─5c8e5b61-137e-401b-96f5-df3c9ffaa9dd
# ╠═f76c12c0-3a8e-4ef1-b8dc-a5b2eaf8a761

ext: jl

max_stars_repo_path: 1a/LU with pivot/pivotTest.jl

max_issues_repo_path: 1a/LU with pivot/pivotTest.jl

max_forks_repo_path: 1a/LU with pivot/pivotTest.jl


text:
#=

  Futoshiki problem in Julia ConstraintSolver.jl

  http://en.wikipedia.org/wiki/Futoshiki
  """
  The puzzle is played on a square grid, such as 5 x 5. The objective
  is to place the numbers 1 to 5 (or whatever the dimensions are) such
  that each row, and column contains each of the digits 1 to 5. Some
  digits may be given at the start. In addition, inequality
  constraints are also initially specifed between some of the squares,
  such that one must be higher or lower than its neighbour. These
  constraints must be honoured as the grid is filled out.
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function futoshiki(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    p = problem[:problem]
    less_than = problem[:less_than]
    n,_ = size(p)

    @variable(model, 1 <= x[1:n,1:n] <= n, Int)

    latin_square(model, x)

    # fill the data
    for i in 1:n, j in 1:n
        if p[i,j] > 0
            @constraint(model, x[i,j] == p[i,j])
        end
    end

    # less than hints
    for (i1,j1,i2,j2) in eachrow(less_than)
        @constraint(model, x[i1,j1] <= x[i2,j2])
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                # println("x:$x_val")
                display(x_val)

            end
        end
    else
        println("status:$status")
    end

    return status
end

futoshiki_problems = Dict(
  :1 => Dict(
    # Example from Tailor model futoshiki.param/futoshiki.param
    #
    # Solution:
    # 5 1 3 2 4
    # 1 4 2 5 3
    # 2 3 1 4 5
    # 3 5 4 1 2
    # 4 2 5 3 1
    #
    # Futoshiki instance, by Andras Salamon
    #
    :problem =>
       resize_matrix([[0,0,3,2,0],  # problem grid
        [0,0,0,0,0],
        [0,0,0,0,0],
        [0,0,0,0,0],
        [0,0,0,0,0]]),

    :less_than =>
       resize_matrix([[1,2,1,1], # [i1,j1, i2,j2] requires that values[i1,j1] < values[i2,j2]
        [1,4,1,5],
        [2,3,1,3],
        [3,3,2,3],
        [3,4,2,4],
        [2,5,3,5],
        [3,2,4,2],
        [4,4,4,3],
        [5,2,5,1],
        [5,4,5,3],
        [5,5,4,5]])
   ),

    # Example from http://en.wikipedia.org/wiki/Futoshiki
    # Solution:
    # 5 4 3 2 1
    # 4 3 1 5 2
    # 2 1 4 3 5
    # 3 5 2 1 4
    # 1 2 5 4 3
    #
:2 => Dict(
    :problem =>
       resize_matrix([[0,0,0,0,0],
        [4,0,0,0,2],
        [0,0,4,0,0],
        [0,0,0,0,4],
        [0,0,0,0,0]]),
    :less_than =>
       resize_matrix([[1,2, 1,1],
        [1,4, 1,3],
        [1,5, 1,4],
        [4,4, 4,5],
        [5,1, 5,2],
        [5,2, 5,3]])
),

    # From http://www.sudoku-puzzles.net/futoshiki06x6.html
   :3 => Dict(
    :problem =>
       resize_matrix([[5,0,2,0,1,0],
        [6,0,1,0,0,0],
        [0,0,3,0,0,0],
        [0,0,4,0,0,0],
        [0,0,0,4,2,5],
        [0,0,0,1,3,0]]),
    :less_than =>
       resize_matrix([[1,6, 2,6],
        [4,1, 3,1]])

    )

)

for p in sort(collect(keys(futoshiki_problems)))
    println("\nProblem $p")
    @time futoshiki(futoshiki_problems[p])
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# maxSizeSlices(slice::Array{Int64,1}) = maxInner(slice, 1, length(slice), div(length(slice), 3), 1)

# # use closure
# function maxInner(slice::Array{Int64}, i, j, left, cycle=0)
#     cache = Dict{Tuple{Int64,Int64,Int64}, Int64}()
#     function inside(slice::Array{Int64}, i, j, left, cycle=0)
#         try
#             return cache[(i,j,left)]
#         catch end
#         left == 1 && return maximum(slice[i:j])
#         left == 0 && return 0
#         (j - i + 1) < (2 * left - 1) && return -2147483648

#         result = max(
#             (slice[j] + inside(slice, i + cycle, j - 2, left - 1)),
#             inside(slice, i, j - 1, left)
#         )
#         cache[(i,j,left)] = result
#         result
#     end
#     inside(slice, i, j, left, cycle)
# end

## another solution of closure
## as same as before solution
function maxSizeSlices(slice::Array{Int64,1})
    a = maxInner()
    a(slice, 1, length(slice), div(length(slice), 3), 1)
end

function maxInner()
    cache = Dict{Tuple{Int64,Int64,Int64}, Int64}()
    function inside(slice::Array{Int64}, i, j, left, cycle=0)
        try
            return cache[(i,j,left)]
        catch end
        left == 1 && return maximum(slice[i:j])
        left == 0 && return 0
        (j - i + 1) < (2 * left - 1) && return -2147483648

        result = max(
            (slice[j] + inside(slice, i + cycle, j - 2, left - 1)),
            inside(slice, i, j - 1, left)
        )
        cache[(i,j,left)] = result
        result
    end
    inside
end

using Test
@test maxSizeSlices([3, 1, 2]) == 3
@test maxSizeSlices([4, 1, 2, 5, 8, 3, 1, 9, 7]) == 21

@test maxSizeSlices([8, 9, 8, 6, 1, 1]) == 16
@test maxSizeSlices([1, 2, 3, 4, 5, 6]) == 10

@test maxSizeSlices([10, 9, 1, 10, 8, 5, 10, 2, 2]) == 30
@test maxSizeSlices([9, 8, 1, 7, 7, 9, 5, 10, 7, 9, 3, 8, 3, 4, 8]) == 45

@test maxSizeSlices([9, 5, 1, 7, 8, 4, 4, 5, 5, 8, 7, 7]) == 30

@test maxSizeSlices([
    2, 4, 3, 1, 10, 6, 1, 2, 10, 2, 8, 9, 4, 8, 8, 8, 10, 3, 9, 10, 7, 9, 4, 5, 4, 3, 1,
    10, 6, 1
]) == 85

@test maxSizeSlices([
    6, 3, 1, 2, 6, 2, 4, 3, 10, 4, 1, 4, 6, 5, 5, 3, 4, 7, 6, 5, 8, 7, 3, 8, 8, 1, 7, 1, 7,
    8
]) == 70

@test maxSizeSlices([
    7, 8, 5, 6, 9, 10, 1, 6, 5, 10, 8, 10, 5, 4, 7, 2, 8, 5, 9, 7, 5, 9, 3, 7, 7, 2, 2, 10,
    7, 6, 4, 6, 5, 7, 7, 9, 6, 8, 10, 7, 5, 7, 2, 5, 4, 9, 6, 10, 10, 2, 10
]) == 150

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
#=
Created on Thurday 20 Feb 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Process the results of the host-parasite data.
=#

using StatsPlots, HypothesisTests, BSON, NamedArrays

BSON.@load "results/hp.bson" Pstrain Pstest M kl_ot_ab  kl_ot_a kl_ot_b kl_ot_softmax kl_ot_neutral γ

# HYPOTHESIS TESTS
# ----------------

Δkl = kl_ot_neutral - kl_ot_ab

mean(Δkl) / (std(Δkl) / √(length(Δkl)))

@show wsr = SignedRankTest(kl_ot_neutral, kl_ot_ab)
@show pvalue(wsr, tail=:right)

# FIGURES
# -------

p_hp = boxplot(kl_ot_neutral, label="neutral", legend=:none)
boxplot!(kl_ot_ab, label="OT (both)")
boxplot!(kl_ot_a, label="OT (parasites fixed)")
boxplot!(kl_ot_b, label="OT (hosts fixed)")
boxplot!(kl_ot_softmax, label="OT (none)")
ylabel!("Kullback-Leibler (nats)")
xticks!(1:5, ["neutral", "OT (both)", "OT (parasites fixed)", "OT (hosts fixed)", "OT (none)"])
title!("Predicted versus observed coupling\nhost-parasite networks")

savefig(p_hp, "figures/hostparasite.svg")

# on the training data

kl_ot_ab = Float64[]
kl_ot_a = Float64[]
kl_ot_b = Float64[]
kl_ot_softmax = Float64[]
kl_ot_neutral = Float64[]

cutoff = 1e-6  # remove sp. with no interactions

for P in Pstrain
    P = Matrix(P)
    P ./=sum(P)
    a, b = marginals(P)
    # remove empty rows/columns
    P = P[a.>cutoff, b.>cutoff]
    Mss = M[a.>cutoff, b.>cutoff]
    a, b = marginals(P)
    # optimal transport both
    Q = optimaltransport(Mss, a, b)
    push!(kl_ot_ab, KL(P, Q))
    Q .= optimaltransport(Mss, a, nothing)
    push!(kl_ot_a, KL(P, Q))
    Q .= optimaltransport(Mss, nothing, b)
    push!(kl_ot_b, KL(P, Q))
    Q .= optimaltransport(Mss, nothing, nothing)
    push!(kl_ot_softmax, KL(P, Q))
    Q .= optimaltransport(a, b)
    push!(kl_ot_neutral, KL(P, Q))
end

p_hp = boxplot(kl_ot_neutral, label="neutral", legend=:none)
boxplot!(kl_ot_ab, label="OT (both)")
boxplot!(kl_ot_a, label="OT (parasites fixed)")
boxplot!(kl_ot_b, label="OT (hosts fixed)")
boxplot!(kl_ot_softmax, label="OT (none)")
ylabel!("Kullback-Leibler (nats)")
xticks!(1:5, ["neutral", "OT (both)", "OT (parasites fixed)", "OT (hosts fixed)", "OT (none)"])
title!("Predicted versus observed coupling\nhost-parasite networks (training)")

savefig(p_hp, "figures/hostparasite_train.svg")

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
include("roblib.jl")

find_mrk(r)


text:
# note: grammar was manually brought to Chomsky normal form, as there were only 2 rules to be modified

function parseRule(str)
    a = split(str, ": ")
    parent = parse(Int64, a[1])
    children = a[2]
    options = []
    for rhs in split(children, " | ")
        r = []
        for child in split(rhs, " ")
            try
                push!(r, parse(Int64, child))
            catch
                push!(r, child[2:end-1])
            end
        end
        println(r)
        @assert(length(r) == 2 || r[1][1] == 'a' || r[1][1] == 'b')
        push!(options, r)
    end
    (parent, options)
end

function findParent(rules, rhs)
    try
        return rules[rhs]
    catch
        return []
    end
    parents = []
    for rule in rules
        for r in rule[2]
            if r == rhs
                push!(parents, rule[1])
            end
        end
    end
    parents
end

# prints solution of any part, if the appropriate modification is made in the input
function x()
    f = open("19.txt")
    rules = Dict{Array{Any,1}, Array{Int64,1}}()
    for i=1:139
        rule = parseRule(readline(f))
        for rhs in rule[2]
            rules[rhs] = push!(get(rules, rhs, []), rule[1])
        end
    end
    println(rules)
    readline(f)
    words = []
    for i=1:490
        push!(words, readline(f))
    end
    println(words)

    c = 0
    for word in words
        print(word)
        n = length(word)
        pyramid = []
        lvl1 = Array{Array{Int64,1}}(undef, n)
        fill!(lvl1, [])
        for i=1:n
            lvl1[i] = findParent(rules, [string(word[i])])
        end
        push!(pyramid, lvl1)
        for i=2:n
            lvl = Array{Array{Int64,1}}(undef, n+1-i)
            for j = 1:n+1-i
                lvl[j] = []
                for k = 1:i-1
                    rhs = [pyramid[k][j], pyramid[i-k][j+k]]
                    for a=1:length(rhs[1]), b = 1:length(rhs[2])
                        rh = [rhs[1][a], rhs[2][b]]
                        append!(lvl[j], findParent(rules, rh))
                    end
                end
            end
            push!(pyramid, lvl)
        end
        if 0 in pyramid[end][1]
            println(": OK")
            c += 1
        else
            println()
        end
    end
    c
end


text:
"""
    gamma_rng(α, β, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Gamma(`α`, `β`) distribution. Optionally you can set a specific seed.

# Notes

The Gamma distribution is given:

``f(x,α,β) = \\frac{β^α x^{α-1} e^{-βx}}{Γ(α)} \\quad x ≥ 0``

# Examples

```julia-repl
julia> gamma_rng(1,1)
1-element Vector{Float64}:
 0.5190236735858542

julia> gamma_rng(1,1,4)
4-element Vector{Float64}:
 0.3035517926878862
 0.5765419737109622
 0.44121996206333797
 0.7325887616559309

julia> gamma_rng(1,1,(2,2))
2×2 Matrix{Float64}:
 0.228818  0.88849
 0.665729  1.01668

```

# References

Law, A. Simulation modeling and analysis, 5th Ed. McGraw Hill Education, Tuscon, 2013.
"""
function gamma_rng(α::Real, β::Real=1, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    seed_setter(seed)
    X = β .* get_gamma_rv(α, shape)
    return X
end

function get_gamma_rv(α::Real, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    seed_setter(seed)
    if α < 0
        throw(ArgumentError("Not implemented for α < 0."))
    elseif α == 1
        return expon_rng(1, shape)
    elseif α < 1
        X = zeros(shape)
        X .= get_gamma_ad.(α)
        return X
    else
        X = zeros(shape)
        X .= get_gamma_gb.(α)
        return X
    end
end

function get_gamma_ad(α::Real; max_iter=300)
    b = (ℯ + α)/ℯ
    i = 0
    while i < max_iter
        U₁ = get_std_uniform()[1]
        P = b*U₁
        if P > 1
            Y = -log((b-P)/α)
            U₂ = get_std_uniform()[1]
            if U₂ ≤ Y^(α-1)
                return Y
            end
        else
            Y = P^(1/α)
            U₂ = get_std_uniform()[1]
            if U₂ ≤ ℯ^(-Y)
                return Y
            end
        end
        i +=1
    end
end

function get_gamma_gb(α::Real; max_iter=300)
    a = 1 / sqrt(2*α - 1)
    b = α - log(4)
    q = α + 1/a
    θ = 4.5
    d = 1 + log(θ)
    i = 0
    while i < max_iter
        U₁, U₂ = get_std_uniform(2)
        V = a * log(U₁/(1-U₁))
        Y = α * ℯ^V
        Z = U₁^2 * U₂
        W = b + q * V - Y
        if W + d - θ * Z ≥ 0
            return Y
        elseif W ≥ log(Z)
            return Y
        end
        i += 1
    end
end


text:
using Plots
using Printf
using SpectralMethodFFT

function xy_kl_init(c::ConfigFFT)

    xy_X = xy_Xgen(c)
    xy_Y = xy_Ygen(c)

    xy_r = sqrt(xy_X^2 + xy_Y^2)
    xy_u = sinc(xy_r / 2)
    kl_u = kl_xy(xy_u)

    return xy_u, kl_u

end

function kl_develop_kl(kl_u::KFunc, dt::Real)

    c = kl_u.config
    xlen, ylen = xlens_xranges(c.xranges)

    kl_K = kl_Kgen(c)
    kl_L = kl_Lgen(c)

    return exp( - 4 * π^2 * dt * (
        (kl_K/xlen)^2 + (kl_L/ylen)^2
    )) * kl_u

end

function main()

    # -------------
    # Parameters
    # -------------
    # *** space ***
    xranges = ((-10., 10.), (-10., 10.))
    ngrids = (128, 128)
    c = ConfigFFT(ngrids, xranges)
    # *** time ***
    t_st = 0.
    t_ed = 2.0
    nt = 200
    dt = (t_ed - t_st) / nt

    # -------------
    # Time Integration
    # -------------
    xy_u, kl_u = xy_kl_init(c)
    xy_X = xy_Xgen(c)
    xy_Y = xy_Ygen(c)

    anim = Animation()
    for it = 0:nt
        t = @sprintf("%3.2f", t_st + it*dt)
        xy_u = xy_kl(kl_u)
        plt = plot(
                xy_X[:,1], xy_Y[1,:], xy_u,
                title="t=$t",
                zlims=(-0.25, 1.1),
                xlabel="x", ylabel="y", zlabel="u",
                st=:wireframe)
        frame(anim, plt)
        if it < nt
            kl_u = kl_develop_kl(kl_u, dt)
        end
    end
    gif(anim, "heat_2d.gif", fps=20)

end

main()

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
using Plots,ApproxFun

#The following samples eigenvalues of an n = 2 Gaussian Unitary Ensemble:

ff=(x,y)->(x-y)^2*exp(-x^2/2.-y^2/2)
f=Fun(ff,[-4.,4.],[-4.,4.])
r=ApproxFun.sample(f,5000)

#We can compare the histogram to the 1-point correlation
plot(sum(f,1)/sum(f))
histogram!(vec(r);normed=true)

#We can compare the histograms of x with the GUE:
using RandomMatrices

histogram(vcat([sqrt(2)eigvalrand(GaussianHermite(2),2) for k=1:10000]...);normed=true)
plot!(sum(f,1)/sum(f))

max_stars_repo_path: examples/GUE Sampling.jl

max_issues_repo_path: examples/GUE Sampling.jl

max_forks_repo_path: examples/GUE Sampling.jl


text:
# Tex
# ---

abstract type LaTeXFormat <: WeaveFormat end

function set_format_options!(docformat::LaTeXFormat; keep_unicode = false, template = nothing, _kwargs...)
    docformat.keep_unicode |= keep_unicode
    docformat.template =
        get_mustache_template(isnothing(template) ? normpath(TEMPLATE_DIR, "md2pdf.tpl") : template)
end

# very similar to export to html
function render_chunk(docformat::LaTeXFormat, chunk::DocChunk)
    out = IOBuffer()
    io = IOBuffer()
    for inline in chunk.content
        if isa(inline, InlineText)
            write(io, inline.content)
        elseif !isempty(inline.rich_output)
            clear_buffer_and_format!(io, out, WeaveMarkdown.latex)
            write(out, addlines(inline.rich_output, inline))
        elseif !isempty(inline.figures)
            write(io, inline.figures[end], inline)
        elseif !isempty(inline.output)
            write(io, addlines(inline.output, inline))
        end
    end
    clear_buffer_and_format!(io, out, WeaveMarkdown.latex)
    out = take2string!(out)
    return unicode2latex(docformat, out)
end

render_output(docformat::LaTeXFormat, output) = unicode2latex(docformat, output, true)

render_code(docformat::LaTeXFormat, code) = unicode2latex(docformat, code, true)

render_termchunk(docformat::LaTeXFormat, chunk) =
    string(docformat.termstart,
            unicode2latex(docformat, chunk.output, true),
            docformat.termend, "\n")

# from julia symbols (e.g. "\bfhoge") to valid latex
const UNICODE2LATEX = let
    function texify(s)
        return if occursin(r"^\\bf[A-Z]$", s)
            replace(s, "\\bf" => "\\bm{\\mathrm{") * "}}"
        elseif startswith(s, "\\bfrak")
            replace(s, "\\bfrak" => "\\bm{\\mathfrak{") * "}}"
        elseif startswith(s, "\\bf")
            replace(s, "\\bf" => "\\bm{\\") * "}"
        elseif startswith(s, "\\frak")
            replace(s, "\\frak" => "\\mathfrak{") * "}"
        else
            s
        end
    end
    Dict(unicode => texify(sym) for (sym, unicode) in REPL.REPLCompletions.latex_symbols)
end

function unicode2latex(docformat::LaTeXFormat, s, escape = false)
    # Check whether to convert at all and return input if not
    docformat.keep_unicode && return s
    for (unicode, latex) in UNICODE2LATEX
        body = "\\ensuremath{$(latex)}"
        target = escape ? string(docformat.escape_starter, body, docformat.escape_closer) : body
        s = replace(s, unicode => target)
    end
    return s
end

function render_figures(docformat::LaTeXFormat, chunk)
    fignames = chunk.figures
    caption = chunk.options[:fig_cap]
    width = chunk.options[:out_width]
    height = chunk.options[:out_height]
    f_pos = chunk.options[:fig_pos]
    f_env = chunk.options[:fig_env]
    result = ""
    figstring = ""

    if isnothing(f_env) && !isnothing(caption)
        f_env = "figure"
    end

    (isnothing(f_pos)) && (f_pos = "!h")
    # Set size
    attribs = ""
    isnothing(width) || (attribs = "width=$(md_length_to_latex(width,"\\linewidth"))")
    (!isempty(attribs) && !isnothing(height)) && (attribs *= ",")
    isnothing(height) || (attribs *= "height=$(md_length_to_latex(height,"\\paperheight"))")

    if !isnothing(f_env)
        result *= "\\begin{$f_env}"
        (!isempty(f_pos)) && (result *= "[$f_pos]")
        result *= "\n"
    end

    for fig in fignames
        if splitext(fig)[2] == ".tex" # Tikz figures
            figstring *= "\\resizebox{$width}{!}{\\input{$fig}}\n"
        else
            if isempty(attribs)
                figstring *= "\\includegraphics{$fig}\n"
            else
                figstring *= "\\includegraphics[$attribs]{$fig}\n"
            end
        end
    end

    # Figure environment
    if !isnothing(caption)
        result *= string("\\center\n", "$figstring", "\\caption{$caption}\n")
    else
        result *= figstring
    end

    if !isnothing(chunk.options[:label]) && !isnothing(f_env)
        label = chunk.options[:label]
        result *= "\\label{fig:$label}\n"
    end

    if !isnothing(f_env)
        result *= "\\end{$f_env}\n"
    end

    return result
end

function md_length_to_latex(def, reference)
    if occursin("%", def)
        _def = tryparse(Float64, replace(def, "%" => ""))
        isnothing(_def) && return def
        perc = round(_def / 100, digits = 2)
        return "$perc$reference"
    end
    return def
end

function render_doc(docformat::LaTeXFormat, body, doc)
    return Mustache.render(
        docformat.template;
        body = body,
        highlight = "",
        tex_deps = docformat.tex_deps,
        [Pair(Symbol(k), v) for (k, v) in doc.header]...,
    )
end

# minted Tex
# ----------

Base.@kwdef mutable struct LaTeXMinted <: LaTeXFormat
    description = "LaTeX using minted package for code highlighting"
    extension = "tex"
    codestart = "\\begin{minted}[texcomments = true, mathescape, fontsize=\\small, xleftmargin=0.5em]{julia}"
    codeend = "\\end{minted}"
    termstart = "\\begin{minted}[texcomments = true, mathescape, fontsize=\\footnotesize, xleftmargin=0.5em]{jlcon}"
    termend = "\\end{minted}"
    outputstart = "\\begin{minted}[texcomments = true, mathescape, fontsize=\\small, xleftmargin=0.5em, frame = leftline]{text}"
    outputend = "\\end{minted}"
    mimetypes = ["application/pdf", "image/png", "text/latex", "text/plain"]
    fig_ext = ".pdf"
    out_width = "\\linewidth"
    out_height = nothing
    fig_pos = "htpb"
    fig_env = "figure"
    # specials
    keep_unicode = false
    template = nothing
    tex_deps = "\\usepackage{minted}"
    # how to escape latex in verbatim/code environment
    escape_starter = "|\$"
    escape_closer = reverse(escape_starter)
end
register_format!("texminted", LaTeXMinted())

# Tex (directly to PDF)
# ---------------------

abstract type WeaveLaTeXFormat <: LaTeXFormat end

function set_format_options!(docformat::WeaveLaTeXFormat; template = nothing, highlight_theme = nothing, keep_unicode = false, _kwargs...)
    docformat.template =
        get_mustache_template(isnothing(template) ? normpath(TEMPLATE_DIR, "md2pdf.tpl") : template)
    docformat.highlight_theme = get_highlight_theme(highlight_theme)
    docformat.keep_unicode |= keep_unicode
end

function render_output(docformat::WeaveLaTeXFormat, output)
    # Highligts has some extra escaping defined, eg of $, ", ...
    output_escaped = sprint(
        (io, x) ->
            Highlights.Format.escape(io, MIME("text/latex"), x, charescape = true),
        output,
    )
    return unicode2latex(docformat, output_escaped, true)
end

function render_code(docformat::WeaveLaTeXFormat, code)
    ret = highlight_code(MIME("text/latex"), code, docformat.highlight_theme)
    unicode2latex(docformat, ret, false)
end

function render_termchunk(docformat::WeaveLaTeXFormat, chunk)
    if should_render(chunk)
        ret = highlight_term(MIME("text/latex"), chunk.output, docformat.highlight_theme)
        unicode2latex(docformat, ret, true)
    else
        ""
    end
end

function render_doc(docformat::WeaveLaTeXFormat, body, doc)
    return Mustache.render(
        docformat.template;
        body = body,
        highlight = get_highlight_stylesheet(MIME("text/latex"), docformat.highlight_theme),
        tex_deps = docformat.tex_deps,
        [Pair(Symbol(k), v) for (k, v) in doc.header]...,
    )
end

Base.@kwdef mutable struct WeaveLaTeX <: WeaveLaTeXFormat
    description = "Weave-styled LaTeX"
    extension = "tex"
    codestart = ""
    codeend = ""
    termstart = codestart
    termend = codeend
    outputstart = "\\begin{lstlisting}"
    outputend = "\\end{lstlisting}\n"
    mimetypes = ["application/pdf", "image/png", "image/jpg", "text/latex", "text/markdown", "text/plain"]
    fig_ext = ".pdf"
    out_width = "\\linewidth"
    out_height = nothing
    fig_pos = nothing
    fig_env = nothing
    # specials
    highlight_theme = nothing
    template = nothing
    keep_unicode = false
    tex_deps = ""
    # how to escape latex in verbatim/code environment
    escape_starter = "(*@"
    escape_closer = "@*)"
end
register_format!("md2tex", WeaveLaTeX())

# will be used by `write_doc`
const DEFAULT_LATEX_CMD = ["xelatex", "-shell-escape", "-halt-on-error"]

Base.@kwdef mutable struct LaTeX2PDF <: ExportFormat
    primaryformat = WeaveLaTeX()
    description = "PDF via LaTeX"
    latex_cmd = DEFAULT_LATEX_CMD
end
register_format!("md2pdf", LaTeX2PDF())
register_format!("minted2pdf", LaTeX2PDF(primaryformat=LaTeXMinted()))

function set_format_options!(docformat::LaTeX2PDF; latex_cmd = DEFAULT_LATEX_CMD, _kwargs...)
    docformat.latex_cmd = latex_cmd
    set_format_options!(docformat.primaryformat; _kwargs...)
end


text:
#=
Created on Thurday 16 Jan 2020
Last update:

@author: Michiel Stock
michielfmstock@gmail.com

Solution of the Monge problem
=#

using STMO
using STMO.OptimalTransport
using Plots, LaTeXStrings

p = scatter(X1[:,1], X1[:,2], color=myorange, label="location cells at t1")
xlabel!("\$x\$")
ylabel!("\$y\$")
scatter!(X2[:,1], X2[:,2], color=mygreen, label="location cells at t2")

C = dist(X1, X2)

perm, cost = monge_brute_force(C)

for (i, j) in enumerate(perm)
    plot!(p, [X1[i,1], X2[j,1]], [X1[i,2], X2[j,2]], label="", color=myred,
                alpha=0.8, lw=2)
end

title!(p, "Monge solution with cost=$cost")

savefig(p, "figures/mongesol.png")


text:
using GeometricIntegrators.Config
using GeometricIntegrators.Solvers
using Test

n = 1
T = Float64

x = [T(π),]
j = reshape(2x, 1, 1)

function F(x::Vector, b::Vector)
    b[:] = x.^2
end

function J(x::Vector, A::Matrix)
    A[:,:] = 2x
end

set_config(:jacobian_autodiff, true)
JPAD = getJacobianParameters(nothing, F, T, n)

set_config(:jacobian_autodiff, false)
JPFD = getJacobianParameters(nothing, F, T, n)

set_config(:jacobian_autodiff, true)
JPUS = getJacobianParameters(J, F, T, n)

@test typeof(JPAD) <: JacobianParametersAD
@test typeof(JPFD) <: JacobianParametersFD
@test typeof(JPUS) <: JacobianParametersUser

function test_jac(j1, j2, atol)
    for i in eachindex(j1,j2)
        @test j1[i] ≈ j2[i] atol=atol
    end
end

jad = zero(j)
jfd = zero(j)
jus = zero(j)

computeJacobian(x, jad, JPAD)
computeJacobian(x, jfd, JPFD)
computeJacobian(x, jus, JPUS)

test_jac(jad, j, eps())
test_jac(jfd, j, 1E-7)
test_jac(jus, j, 0)

jad2 = zero(j)
jfd2 = zero(j)

computeJacobianAD(x, jad2, F)
computeJacobianFD(x, jfd2, F, get_config(:jacobian_fd_ϵ))

test_jac(jad, jad2, 0)
test_jac(jfd, jfd2, 0)


text:
using IntervalOptimisation
using Test

const HeapedVector = IntervalOptimisation.HeapedVector

@testset "heap" begin

	@testset "with standard heaping " begin

    h=HeapedVector([9, 7, 5, 2, 3])

    @test h.data == [2, 3, 7, 9, 5]

    push!(h, 1)
    @test h.data == [1, 3, 2, 9, 5, 7]

    popfirst!(h)
    @test h.data == [2, 3, 7, 9, 5]

  end

  @testset "with ordering function" begin

    h=HeapedVector( [(9,"f"), (7,"c"), (5,"e"), (2,"d"), (3,"b")] , x->x[2])

    @test h.data == [(3, "b"), (7, "c"), (5, "e"), (9, "f"), (2, "d")]

    push!(h, (1,"a"))
    @test h.data == [(1, "a"), (7, "c"), (3, "b"), (9, "f"), (2, "d"), (5, "e")]

    popfirst!(h)
    @test h.data == [(3, "b"), (7, "c"), (5, "e"), (9, "f"), (2, "d")]

  end

end


text:
#=

  Euler #4 in Julia.

  Problem 4
  """
  A palindromic number reads the same both ways. The largest palindrome made
  from the product of two 2-digit numbers is 9009 = 91 × 99.

  Find the largest palindrome made from the product of two 3-digit numbers.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
include("Euler.jl")

# 0.00958293s
function euler4a()
    max = 0
    from = 100
    to = 999
    for i in from:to, j in i:to
        ij = i*j
        if ij > max && palindromic_number(ij)
            max = ij
        end
    end
    return max;
end

# 0.313482985s
function euler4b()
    from = 100
    to = 999
    a = []

    # Note: It seems that one cannot assign a variable (ij=i*j)
    # in a function pipe
    (from:to).|> i->
        (i:to).|> j->
            # ij = i * j; # this don't work! (Why?)
            if palindromic_number(i*j)
                push!(a,i*j)
            end
    return maximum(a)
end

# 0.085608391s
function euler4c()
    from = 100
    to = 999

    a = (from:to).|>i->(i:to).|>j->i*j
    # vcat(a...) flattens the array
    return maximum(filter(ij->palindromic_number(ij),vcat(a...)))
end

# "Full" chaining, much slower
# 0.12105172s
function euler4d()
    from = 100
    to = 999
    a = (from:to).|>i->(i:to).|>j->i*j
    return vcat(a...)|>x->filter(ij->palindromic_number(ij),x)|>maximum
end

run_euler(euler4a)
# run_euler(euler4b)
# run_euler(euler4c)
# run_euler(euler4d)


text:
include("run_sim.jl")

function waiting_diff(λ_H, λ_E, N, M, t_1, time_function, num, T)
    table = DataFrame(λ_H = Array{Float64}(0), λ_E = Array{Float64}(0), w_H = Array{Float64}(0),
                        wait = Array{Float64}(0), wait1 = Array{Float64}(0))
    for (e, λ_E) in enumerate(linspace(1/num,N,num))
        e%(num/10) == 0 && println("$e out of $num")
        for (h, λ_H) in enumerate(linspace(1/num,M,num))
            t_sum = time_function(λ_E + λ_E_1, λ_H + λ_H_1, p_E, p_H, T, tr_prioH)
            t_2 = time_function(λ_E, λ_H, p_E, p_H, T, tr_prioH)
            #t_2 = mean(h_2[Int(T*3/4):T, :w_H])
            #t_sum = mean(h_sum[Int(T*3/4):T, :w_H])
            push!(table, (λ_H, λ_E, t_sum, t_sum - t_1, t_sum - t_2))
        end
    end
    return table
end

function find_boundary(table)
    table[:improvement] = table[:wait] .< 0;
    ve = zeros(Float64, num);
    u = linspace(1/num, N, num)
    for (i, k) in enumerate(linspace(1/num, M, num))
        x = find(! table[table[:λ_H] .== k, :improvement])
        if x == []
            ve[i] = 0
        else
            ve[i] = u[maximum(x)]
        end
    end
    v = []
    for i in 1:num
        append!(v, ve)
    end
    table[:v] = v;
end

# Case lambda_H < lambda_E
λ_E_1 = 1.3
λ_H_1 = 1
p_H = 0.02
p_E = 0.5
N = 3
M = 5
d = 0
T = 10000
num = 50

# Find waiting time for market 1 alone
tic()
t_1 = run_sim(λ_E_1, λ_H_1, p_E, p_H, T, tr_prioH)
toc()

tic()
table1 = waiting_diff(λ_H_1, λ_E_1, N, M, t_1, run_sim, num, T);
toc()

find_boundary(table1)

writetable("csv/fig_2.1.csv", table1)

# Case lambda_H > lambda_E
λ_E_1 = 1
λ_H_1 = 1.3

t_1 = run_sim(λ_E_1, λ_H_1, p_E, p_H, T, tr_prioH);

tic()
table2 = waiting_diff(λ_H_1, λ_E_1, N, M, t_1, run_sim, num, T);
toc()

find_boundary(table2)

writetable("csv/fig_2.2.csv", table2)


text:
export AbstractDataSet, MemoryDataSet, LazyDataSet, rowiter, coliter, entryiter, set!, DataEntry, TrainAndTestSplit, randombatch

import ..@Fun

struct DataEntry{InputType, OutputType}
    inputs::Array{InputType}
    outputs::Array{OutputType}
    DataEntry(inputs::Array{I}, outputs::Array{O}) where {I, O} = new{I, O}(inputs, outputs)
end

abstract type AbstractDataSet{InputType, OutputType} end

Base.size(set::AbstractDataSet) = length(set)
Base.eltype(::Type{A}) where {I, O, A <: AbstractDataSet{I, O}} = DataEntry{I, O}
Base.iterate(set::AbstractDataSet, state) = state <= length(set) ? (set[state], state + 1) : nothing
Base.iterate(set::AbstractDataSet) = Base.iterate(set, 1)
Base.reset(set::AbstractDataSet) = ()
Base.eachindex(set::AbstractDataSet) = 1:length(set)

function Base.print(io::IO, set::AbstractDataSet)
    for item in set
        println(io, "Inputs:", item.inputs, "\tOutputs:", item.outputs)
    end
end

function Base.string(set::AbstractDataSet)
    str = ""
    for item in set
        str *= "Inputs:" * item.inputs * "\tOutputs:" * item.outputs * "\n"
    end
    str
end

function TrainAndTestSplit(x::AbstractDataSet; trainFrac=.8)
    trainLen = Int(round(length(x) * trainFrac))
    return (view(x, 1:trainLen) , view(x, trainLen:length(x)))
end


text:
using RecipesBase

RecipesBase.@recipe function f(var::Symbol, asm::SEAssemblage)
    replace(asm.site.sitestats[var], missing => NaN), getcoords(places(asm))
end


text:
using Test

A    = sprandn(100,10,.1)
pFor = LSparam(A,[])
m0   = randn(10)

# checkDerivative for ForwardProbTypes
passedDerivative, = checkDerivative(m0,pFor)
@test passedDerivative

# adjoint test
passedAdjoint, = adjointTest(m0,pFor)
@test passedAdjoint


text:
module ProximalCore

using LinearAlgebra

is_convex(::Type) = false
is_convex(::T) where T = is_convex(T)

is_generalized_quadratic(::Type) = false
is_generalized_quadratic(::T) where T = is_generalized_quadratic(T)

"""
    gradient!(y, f, x)

In-place gradient (and value) of `f` at `x`.

The gradient is written to the (pre-allocated) array `y`, which should have the same shape/size as `x`.

Returns the value `f` at `x`.

See also: [`gradient`](@ref).
"""
gradient!

"""
    gradient(f, x)

Gradient (and value) of `f` at `x`.

Return a tuple `(y, fx)` consisting of
- `y`: the gradient of `f` at `x`
- `fx`: the value of `f` at `x`

See also: [`gradient!`](@ref).
"""
function gradient(f, x)
    y = similar(x)
    fx = gradient!(y, f, x)
    return y, fx
end

"""
    prox!(y, f, x, gamma=1)

In-place proximal mapping for `f`, evaluated at `x`, with stepsize `gamma`.

The proximal mapping is defined as
```math
\\mathrm{prox}_{\\gamma f}(x) = \\arg\\min_z \\left\\{ f(z) + \\tfrac{1}{2\\gamma}\\|z-x\\|^2 \\right\\}.
```
The result is written to the (pre-allocated) array `y`, which should have the same shape/size as `x`.

Returns the value of `f` at `y`.

See also: [`prox`](@ref).
"""
prox!

prox!(y, f, x) = prox!(y, f, x, 1)

"""
    prox(f, x, gamma=1)

Proximal mapping for `f`, evaluated at `x`, with stepsize `gamma`.

The proximal mapping is defined as
```math
\\mathrm{prox}_{\\gamma f}(x) = \\arg\\min_z \\left\\{ f(z) + \\tfrac{1}{2\\gamma}\\|z-x\\|^2 \\right\\}.
```

Returns a tuple `(y, fy)` consisting of
- `y`: the output of the proximal mapping of `f` at `x` with stepsize `gamma`
- `fy`: the value of `f` at `y`

See also: [`prox!`](@ref).
"""
function prox(f, x, gamma=1)
    y = similar(x)
    fy = prox!(y, f, x, gamma)
    return y, fy
end

struct Zero end

(::Zero)(x) = real(eltype(x))(0)

function gradient!(y, f::Zero, x)
    y .= eltype(x)(0)
    return f(x)
end

function prox!(y, ::Zero, x, gamma)
    y .= x
    return real(eltype(y))(0)
end

is_convex(::Type{Zero}) = true
is_generalized_quadratic(::Type{Zero}) = true

struct IndZero end

function (::IndZero)(x)
    R = real(eltype(x))
    if iszero(x)
        return R(Inf)
    end
    return R(0)
end

is_convex(::Type{IndZero}) = true
is_generalized_quadratic(::Type{IndZero}) = true

function prox!(y, ::IndZero, x, gamma)
    R = real(eltype(x))
    y .= R(0)
    return R(0)
end

struct ConvexConjugate{T}
    f::T
end

is_convex(::Type{<:ConvexConjugate}) = true
is_generalized_quadratic(::Type{ConvexConjugate{T}}) where T = is_generalized_quadratic(T)

function prox_conjugate!(y, u, f, x, gamma)
    u .= x ./ gamma
    v = prox!(y, f, u, 1 / gamma)
    v = real(dot(x, y)) - gamma * real(dot(y, y)) - v
    y .= x .- gamma .* y
    return v
end

prox_conjugate!(y, f, x, gamma) = prox_conjugate!(y, similar(x), f, x, gamma)

prox!(y, g::ConvexConjugate, x, gamma) = prox_conjugate!(y, g.f, x, gamma)

convex_conjugate(f) = ConvexConjugate(f)
convex_conjugate(::Zero) = IndZero()
convex_conjugate(::IndZero) = Zero()

end # module


text:
d = 3    # number of dimensions
K = 50   # number of points along each dimension
N = 100000  # number of points at which to interpolate

A = rand([K for i = 1:d]...)    # filtered coefficients
B = rand(N,d)                    # points at which to evaluate the splines
max(B, minimum(A)+0.01)
min(B, maximum(A)-0.01)

println("Comparison (d=",d,", K=", K, ", N=", N,")")

# interpolations.jl
using Interpolations
itp = interpolate(A, BSpline(Quadratic(Reflect())), OnGrid())

# what is the fastest way to operate on a list of points ?
LB = Vector{Float64}[copy(slice(B,i,:)) for i=1:size(B,1)]

function evaluate(itp, LB)
    s = 0.0
    for i=1:size(LB,1)
        I = LB[i]
        s += itp[I[1],I[2],I[3]]
    end
    return s
end

total_2 = evaluate(itp, LB)

println("interpolations.jl (quadratic)")
@time total_2 = evaluate(itp, LB)
exit()

# splines.jl
using splines

# dummy boundaries
a = [0.0 for i=1:d]
b = [1.0 for i=1:d]
orders = Int64[K for i=1:d]
coeffs = filter_coeffs(a,b,orders,A)

# jit
println("splines.jl: ")
values_1 = eval_UC_spline(a,b,orders,coeffs,B);
tot = sum(values_1)
@time  values_1 = eval_UC_spline(a,b,orders,coeffs,B); tot = sum(values_1)
# not sure whether the summation is taken into accoun
println("splines.jl (with derivatives): ")
values_1, dvalues_1 = eval_UC_spline_G(a,b,orders,coeffs,B); tot = sum(values_1)

@time  values_1, dvalues_1 = eval_UC_spline_G(a,b,orders,coeffs,B);total_1 = sum(values_1)


text:
@testset "smoothing" begin

# Need to test non-box object, such as a sphere and see if subpixel smoothing generates
# the expected smoothed material parameters (instead of simply using kottke_input_simple or
# amean_param or hmean_param.)

@testset "smoothing, box with odd number of voxels" begin
    # Create a grid.
    isbloch = [true, true, true]
    # isbloch = [true, false, false]
    lprim = ([-1.5, -0.5, 0.5, 1.5], [-1.5, -0.5, 0.5, 1.5], [-1.5, -0.5, 0.5, 1.5])
    g3 = Grid(lprim, isbloch)
    N = g3.N

    # Create materials.
    εvac = 1.0
    vac = Material("Vacuum", ε=εvac)

    εdiel = 2.0
    diel = Material("Dielectric", ε=εdiel)

    # Create objects.
    dom_vac = Object(Cuboid(g3.bounds), vac)
    obj_diel = Object(Cuboid([0,0,0], [1,1,1]), diel)
    # obj_diel = Object(Ball([0,0,0], 1), diel)

    # Add objects.
    oind2shp = Shape3[]
    oind2εind = ParamInd[]
    oind2μind = ParamInd[]
    εind2ε = S²ComplexF3[]
    μind2μ = S²ComplexF3[]

    add_obj!(oind2shp, (oind2εind,oind2μind), (εind2ε,μind2μ), dom_vac, obj_diel)

    # Construct arguments and call assign_param!.
    N = g3.N
    ε3d = create_param_array(N)
    εxx_oind3d = create_oind_array(N)
    εyy_oind3d = create_oind_array(N)
    εzz_oind3d = create_oind_array(N)
    εoo_oind3d = create_oind_array(N)

    μ3d = create_param_array(N)
    μxx_oind3d = create_oind_array(N)
    μyy_oind3d = create_oind_array(N)
    μzz_oind3d = create_oind_array(N)
    μoo_oind3d = create_oind_array(N)

    boundft = SVec(EE,EE,EE)
    assign_param!(ε3d, (μxx_oind3d,μyy_oind3d,μzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(ε3d, tuple(μoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, tuple(εoo_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)

    # To-do: test the sanity of the assigned param3d here.  It is relatively easy, and it was very helpful.

    # Perform smoothing.
    smooth_param!(ε3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)
    smooth_param!(ε3d, tuple(εoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)

    ε3dred = view(ε3d, 1:N[1], 1:N[2], 1:N[3], 1:3, 1:3)

    # Construct an expected ε3d.
    ε3dexp = Array{ComplexF64}(undef,3,3,3,3,3)
    rvol = 0.5  # all nonzero rvol used in this test is 0.5
    εh = 1 / (rvol/εdiel + (1-rvol)/εvac)  # harmonic average
    εa = rvol*εdiel + (1-rvol)*εvac  # arithmetic average

    # Initialize ε3dexp.
    for k = 1:N[3], j = 1:N[2], i = 1:N[1]
        ε3dexp[i,j,k,:,:] = εvac * Matrix(I,3,3)
    end

    # Yee's cell at (2,2,2)
    nout = normalize([-1,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,:,:] = εsm  # corner of (2,2,2) cell
    nout = normalize([0,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,1,1] = εsm[1,1]  # x-edge of (2,2,2) cell
    nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,2,2] = εsm[2,2]  # y-edge of (2,2,2) cell
    nout = normalize([-1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,3,3] = εsm[3,3]  # z-edge of (2,2,2) cell

    # Yee's cell at (3,2,2)
    nout = normalize([1,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,:,:] = εsm
    ε3dexp[3,2,2,1,1] = εvac
    nout = normalize([1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,2,2] = εsm[2,2]
    nout = normalize([1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,3,3] = εsm[3,3]

    # Yee's cell at (2,3,2)
    nout = normalize([-1,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,:,:] = εsm
    nout = normalize([0,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,1,1] = εsm[1,1]
    ε3dexp[2,3,2,2,2] = εvac
    nout = normalize([-1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,3,3] = εsm[3,3]

    # Yee's cell at (3,3,2)
    nout = normalize([1,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,:,:] = εsm
    ε3dexp[3,3,2,1,1] = εvac
    ε3dexp[3,3,2,2,2] = εvac
    nout = normalize([1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,3,3] = εsm[3,3]

    # Yee's cell at (2,2,3)
    nout = normalize([-1,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,:,:] = εsm
    nout = normalize([0,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,1,1] = εsm[1,1]
    nout = normalize([-1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,2,2] = εsm[2,2]
    ε3dexp[2,2,3,3,3] = εvac

    # Yee's cell at (3,2,3)
    nout = normalize([1,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,:,:] = εsm
    ε3dexp[3,2,3,1,1] = εvac
    nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,2,2] = εsm[2,2]
    ε3dexp[3,2,3,3,3] = εvac

    # Yee's cell at (2,3,3)
    nout = normalize([-1,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,:,:] = εsm
    nout = normalize([0,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,1,1] = εsm[1,1]
    ε3dexp[2,3,3,2,2] = εvac
    ε3dexp[2,3,3,3,3] = εvac

    # Yee's cell at (3,3,3)
    nout = normalize([1,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,3,:,:] = εsm
    ε3dexp[3,3,3,1,1] = εvac
    ε3dexp[3,3,3,2,2] = εvac
    ε3dexp[3,3,3,3,3] = εvac

    for k = 1:N[3], j = 1:N[2], i = 1:N[1]
        # @info "(i,j,k) = $((i,j,k))"  # uncomment this to know where test fails
        @test @view(ε3dred[i,j,k,:,:]) ≈ @view(ε3dexp[i,j,k,:,:])
        # @info "$(@view(ε3dred[i,j,k,:,:]))"
        # @test issymmetric(@view(ε3dred[i,j,k,:,:]))
        A = @view(ε3dred[i,j,k,:,:])
        @test A ≈ transpose(A)
    end
end  # @testset "smoothing, box with odd number of voxels"

@testset "smoothing, box with even number of voxels" begin
    # Create a grid.
    lprim = ([-2,-1,0,1,2], [-2,-1,0,1,2], [-2,-1,0,1,2])
    isbloch = [true, true, true]
    g3 = Grid(lprim, isbloch)
    N = g3.N

    # Create materials.
    εvac = 1.0
    vac = Material("Vacuum", ε=εvac)

    εdiel = 2.0
    diel = Material("Dielectric", ε=εdiel)

    # Create objects.
    dom_vac = Object(Cuboid(g3.bounds), vac)
    obj_diel = Object(Cuboid([0,0,0], [2,2,2]), diel)
    # obj_diel = Object(Ball([0,0,0], 1), diel)

    # Add objects.
    oind2shp = Shape3[]
    oind2εind = ParamInd[]
    oind2μind = ParamInd[]
    εind2ε = S²ComplexF3[]
    μind2μ = S²ComplexF3[]

    add_obj!(oind2shp, (oind2εind,oind2μind), (εind2ε,μind2μ), dom_vac, obj_diel)

    # Construct arguments and call assign_param!.
    N = g3.N
    ε3d = create_param_array(N)
    εxx_oind3d = create_oind_array(N)
    εyy_oind3d = create_oind_array(N)
    εzz_oind3d = create_oind_array(N)
    εoo_oind3d = create_oind_array(N)

    μ3d = create_param_array(N)
    μxx_oind3d = create_oind_array(N)
    μyy_oind3d = create_oind_array(N)
    μzz_oind3d = create_oind_array(N)
    μoo_oind3d = create_oind_array(N)

    boundft = SVec(EE,EE,EE)
    assign_param!(ε3d, (μxx_oind3d,μyy_oind3d,μzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(ε3d, tuple(μoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, tuple(εoo_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)

    # To-do: test the sanity of the assigned param3d here.  It is relatively easy, and it was very helpful.

    # Perform smoothing.
    smooth_param!(ε3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)
    smooth_param!(ε3d, tuple(εoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)

    ε3dred = view(ε3d, 1:N[1], 1:N[2], 1:N[3], 1:3, 1:3)

    # Construct an expected ε3d.
    ε3dexp = Array{ComplexF64}(undef,4,4,4,3,3)
    rvol = 0.5  # all nonzero rvol used in this test is 0.5
    εh = 1 / (rvol/εdiel + (1-rvol)/εvac)  # harmonic average
    εa = rvol*εdiel + (1-rvol)*εvac  # arithmetic average

    # Initialize ε3dexp.
    for k = 1:N[3], j = 1:N[2], i = 1:N[1]
        ε3dexp[i,j,k,:,:] = εvac * Matrix(I,3,3)
    end

    ## k = 2
    # Yee's cell at (2,2,2)
    nout = normalize([-1,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,:,:] = εsm  # corner of (2,2,2) cell
    nout = normalize([0,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,1,1] = εsm[1,1]  # x-edge of (2,2,2) cell
    nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,2,2] = εsm[2,2]  # y-edge of (2,2,2) cell
    nout = normalize([-1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,2,3,3] = εsm[3,3]  # z-edge of (2,2,2) cell

    # Yee's cell at (3,2,2)
    nout = normalize([0,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,:,:] = εsm
    nout = normalize([0,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,1,1] = εsm[1,1]
    nout = normalize([0,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,2,2] = εsm[2,2]
    nout = normalize([0,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,3,3] = εsm[3,3]

    # Yee's cell at (4,2,2)
    nout = normalize([1,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,2,:,:] = εsm
    ε3dexp[4,2,2,1,1] = εvac
    nout = normalize([1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,2,2,2] = εsm[2,2]
    nout = normalize([1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,2,3,3] = εsm[3,3]

    # Yee's cell at (2,3,2)
    nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,:,:] = εsm
    nout = normalize([0,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,1,1] = εsm[1,1]
    nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,2,2] = εsm[2,2]
    nout = normalize([-1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,3,3] = εsm[3,3]

    # Yee's cell at (3,3,2)
    nout = normalize([0,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,:,:] = εsm
    nout = normalize([0,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,1,1] = εsm[1,1]
    nout = normalize([0,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,2,2] = εsm[2,2]
    ε3dexp[3,3,2,3,3] = εdiel

    # Yee's cell at (4,3,2)
    nout = normalize([1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,2,:,:] = εsm
    ε3dexp[4,3,2,1,1] = εvac
    nout = normalize([1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,2,2,2] = εsm[2,2]
    nout = normalize([1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,2,3,3] = εsm[3,3]

    # Yee's cell at (2,4,2)
    nout = normalize([-1,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,2,:,:] = εsm
    nout = normalize([0,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,2,1,1] = εsm[1,1]
    ε3dexp[2,4,2,2,2] = εvac
    nout = normalize([-1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,2,3,3] = εsm[3,3]

    # Yee's cell at (3,4,2)
    nout = normalize([0,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,2,:,:] = εsm
    nout = normalize([0,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,2,1,1] = εsm[1,1]
    ε3dexp[3,4,2,2,2] = εvac
    nout = normalize([0,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,2,3,3] = εsm[3,3]

    # Yee's cell at (4,4,2)
    nout = normalize([1,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,4,2,:,:] = εsm
    ε3dexp[4,4,2,1,1] = εvac
    ε3dexp[4,4,2,2,2] = εvac
    nout = normalize([1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,4,2,3,3] = εsm[3,3]

    ## k = 3
    # Yee's cell at (2,2,3)
    nout = normalize([-1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,:,:] = εsm
    nout = normalize([0,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,1,1] = εsm[1,1]
    nout = normalize([-1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,2,2] = εsm[2,2]
    nout = normalize([-1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,3,3] = εsm[3,3]

    # Yee's cell at (3,2,3)
    nout = normalize([0,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,:,:] = εsm
    nout = normalize([0,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,1,1] = εsm[1,1]
    ε3dexp[3,2,3,2,2] = εdiel
    nout = normalize([0,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,3,3] = εsm[3,3]

    # Yee's cell at (4,2,3)
    nout = normalize([1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,3,:,:] = εsm
    ε3dexp[4,2,3,1,1] = εvac
    nout = normalize([1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,3,2,2] = εsm[2,2]
    nout = normalize([1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,3,3,3] = εsm[3,3]

    # Yee's cell at (2,3,3)
    nout = normalize([-1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,:,:] = εsm
    ε3dexp[2,3,3,1,1] = εdiel
    nout = normalize([-1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,2,2] = εsm[2,2]
    nout = normalize([-1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,3,3] = εsm[3,3]

    # Yee's cell at (3,3,3)
    ε3dexp[3,3,3,1,1] = εdiel
    ε3dexp[3,3,3,2,2] = εdiel
    ε3dexp[3,3,3,3,3] = εdiel

    # Yee's cell at (4,3,3)
    nout = normalize([1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,3,:,:] = εsm
    ε3dexp[4,3,3,1,1] = εvac
    nout = normalize([1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,3,2,2] = εsm[2,2]
    nout = normalize([1,0,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,3,3,3] = εsm[3,3]

    # Yee's cell at (2,4,3)
    nout = normalize([-1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,3,:,:] = εsm
    nout = normalize([0,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,3,1,1] = εsm[1,1]
    ε3dexp[2,4,3,2,2] = εvac
    nout = normalize([-1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,3,3,3] = εsm[3,3]

    # Yee's cell at (3,4,3)
    nout = normalize([0,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,3,:,:] = εsm
    nout = normalize([0,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,3,1,1] = εsm[1,1]
    ε3dexp[3,4,3,2,2] = εvac
    nout = normalize([0,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,3,3,3] = εsm[3,3]

    # Yee's cell at (4,4,3)
    nout = normalize([1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,4,3,:,:] = εsm
    ε3dexp[4,4,3,1,1] = εvac
    ε3dexp[4,4,3,2,2] = εvac
    nout = normalize([1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,4,3,3,3] = εsm[3,3]

    ## k = 4
    # Yee's cell at (2,2,4)
    nout = normalize([-1,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,4,:,:] = εsm
    nout = normalize([0,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,4,1,1] = εsm[1,1]
    nout = normalize([-1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,4,2,2] = εsm[2,2]
    ε3dexp[2,2,4,3,3] = εvac

    # Yee's cell at (3,2,4)
    nout = normalize([0,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,4,:,:] = εsm
    nout = normalize([0,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,4,1,1] = εsm[1,1]
    nout = normalize([0,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,4,2,2] = εsm[2,2]
    ε3dexp[3,2,4,3,3] = εvac

    # Yee's cell at (4,2,4)
    nout = normalize([1,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,4,:,:] = εsm
    ε3dexp[4,2,4,1,1] = εvac
    nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,2,4,2,2] = εsm[2,2]
    ε3dexp[4,2,4,3,3] = εvac

    # Yee's cell at (2,3,4)
    nout = normalize([-1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,4,:,:] = εsm
    nout = normalize([0,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,4,1,1] = εsm[1,1]
    nout = normalize([-1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,4,2,2] = εsm[2,2]
    ε3dexp[2,3,4,3,3] = εvac

    # Yee's cell at (3,3,4)
    nout = normalize([0,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,4,:,:] = εsm
    nout = normalize([0,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,4,1,1] = εsm[1,1]
    nout = normalize([0,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,4,2,2] = εsm[2,2]
    ε3dexp[3,3,4,3,3] = εvac

    # Yee's cell at (4,3,4)
    nout = normalize([1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,4,:,:] = εsm
    ε3dexp[4,3,4,1,1] = εvac
    nout = normalize([1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,3,4,2,2] = εsm[2,2]
    ε3dexp[4,3,4,3,3] = εvac

    # Yee's cell at (2,4,4)
    nout = normalize([-1,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,4,:,:] = εsm
    nout = normalize([0,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,4,4,1,1] = εsm[1,1]
    ε3dexp[2,4,4,2,2] = εvac
    ε3dexp[2,4,4,3,3] = εvac

    # Yee's cell at (3,4,4)
    nout = normalize([0,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,4,:,:] = εsm
    nout = normalize([0,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,4,4,1,1] = εsm[1,1]
    ε3dexp[3,4,4,2,2] = εvac
    ε3dexp[3,4,4,3,3] = εvac

    # Yee's cell at (4,4,4)
    nout = normalize([1,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[4,4,4,:,:] = εsm
    ε3dexp[4,4,4,1,1] = εvac
    ε3dexp[4,4,4,2,2] = εvac
    ε3dexp[4,4,4,3,3] = εvac

    for k = 1:N[3], j = 1:N[2], i = 1:N[1]
        # @info "(i,j,k) = $((i,j,k))"  # uncomment this to know where test fails
        @test @view(ε3dred[i,j,k,:,:]) ≈ @view(ε3dexp[i,j,k,:,:])
        # @test issymmetric(@view(ε3dred[i,j,k,:,:]))
        A = @view(ε3dred[i,j,k,:,:])
        @test A ≈ transpose(A)
    end
end  # @testset "smoothing, box with even number of voxels"

# Not sure if this test was completed.  It was commented out.  When uncommented, the test failed.
# @testset "smoothing, inspired from SALT3D/test/usage2d_phcwg_test_sym" begin
#     # Create a grid.
#     lprim = ([-0.5,0.5], [-2,-1,0,1,2], [-0.5,0.5])
#     isbloch = [true, true, true]
#     g3 = Grid(lprim, isbloch)
#     N = g3.N
#
#     # Create materials.
#     εvac = 1.0
#     vac = EncodedMaterial(PRIM, Material("Vacuum", ε=εvac))
#
#     εdiel = 2.0
#     diel = EncodedMaterial(PRIM, Material("Dielectric", ε=εdiel))
#
#     # Create objects.
#     dom_vac = Object(Cuboid(g3.bounds), vac)
#     obj_diel = Object(Cuboid([0,0,0], [1,2,1]), diel)
#     # obj_diel = Object(Ball([0,0,0], 1), diel)
#
#     # Add objects.
#     ovec = Object{3}[]
#     paramset = (S²ComplexF3[], S²ComplexF3[])
#     add_obj!(ovec, paramset, dom_vac, obj_diel)
#
#     # Construct arguments and call assign_param!.
#     param3d = create_param_array(N)
#     obj3d = create_p_storage(Object{3}, N)
#     pind3d = create_p_storage(ParamInd, N)
#     oind3d = create_p_storage(ObjInd, N)
#
#     assign_param!(param3d, obj3d, pind3d, oind3d, ovec, g3.ghosted.τl, g3.isbloch)
#     smooth_param!(param3d, obj3d, pind3d, oind3d, g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)
#
#     ε3d = view(param3d[nPR], 1:N[nX], 1:N[nY], 1:N[nZ], 1:3, 1:3)
#
#     # Construct an expected ε3d.
#     ε3dexp = Array{ComplexF64}(undef,4,4,4,3,3)
#     rvol = 0.5  # all nonzero rvol used in this test is 0.5
#     εh = 1 / (rvol/εdiel + (1-rvol)/εvac)  # harmonic average
#     εa = rvol*εdiel + (1-rvol)*εvac  # arithmetic average
#
#     # Initialize ε3dexp.
#     for k = 1:4, j = 1:4, i = 1:4
#         ε3dexp[i,j,k,:,:] = εvac * Matrix(I,3,3)
#     end
#
#     # Yee's cell at (1,2,1)
#     nout = normalize([-1,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[1,2,1,:,:] = εsm  # corner of (2,2,2) cell
#     nout = normalize([0,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[1,2,1,1,1] = εsm[1,1]  # x-edge of (2,2,2) cell
#     nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[1,2,1,2,2] = εsm[2,2]  # y-edge of (2,2,2) cell
#     nout = normalize([-1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[1,2,1,3,3] = εsm[3,3]  # z-edge of (2,2,2) cell
#
#     # Yee's cell at (3,2,2)
#     nout = normalize([1,-1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,:,:] = εsm
#     ε3dexp[3,2,2,1,1] = εvac
#     nout = normalize([1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,2,2] = εsm[2,2]
#     nout = normalize([1,-1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,2,3,3] = εsm[3,3]
#
#     # Yee's cell at (2,3,2)
#     nout = normalize([-1,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,:,:] = εsm
#     nout = normalize([0,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,1,1] = εsm[1,1]
#     ε3dexp[2,3,2,2,2] = εvac
#     nout = normalize([-1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,2,3,3] = εsm[3,3]
#
#     # Yee's cell at (3,3,2)
#     nout = normalize([1,1,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,:,:] = εsm
#     ε3dexp[3,3,2,1,1] = εvac
#     ε3dexp[3,3,2,2,2] = εvac
#     nout = normalize([1,1,0]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,2,3,3] = εsm[3,3]
#
#     # Yee's cell at (2,2,3)
#     nout = normalize([-1,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,:,:] = εsm
#     nout = normalize([0,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,1,1] = εsm[1,1]
#     nout = normalize([-1,0,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,2,3,2,2] = εsm[2,2]
#     ε3dexp[2,2,3,3,3] = εvac
#
#     # Yee's cell at (3,2,3)
#     nout = normalize([1,-1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,:,:] = εsm
#     ε3dexp[3,2,3,1,1] = εvac
#     nout = normalize([-1,0,-1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,2,3,2,2] = εsm[2,2]
#     ε3dexp[3,2,3,3,3] = εvac
#
#     # Yee's cell at (2,3,3)
#     nout = normalize([-1,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,:,:] = εsm
#     nout = normalize([0,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[2,3,3,1,1] = εsm[1,1]
#     ε3dexp[2,3,3,2,2] = εvac
#     ε3dexp[2,3,3,3,3] = εvac
#
#     # Yee's cell at (3,3,3)
#     nout = normalize([1,1,1]); P = nout*nout'; εsm = εh * P + εa * (I-P); ε3dexp[3,3,3,:,:] = εsm
#     ε3dexp[3,3,3,1,1] = εvac
#     ε3dexp[3,3,3,2,2] = εvac
#     ε3dexp[3,3,3,3,3] = εvac
#
#     for k = 1:3, j = 1:3, i = 1:3
#         # @info "(i,j,k) = $((i,j,k))"  # uncomment this to know where test fails
#         @test ε3d[i,j,k,:,:] ≈ ε3dexp[i,j,k,:,:]
#     end
# end

# Need to include a test with three objects with mat1, mat2, mat1.  In the wrong code I
# assumed that I could find the foreground object by sorting voxel corners in terms of pind,
# but this test would show that that is not the case.

end  # @testset "smoothing"


text:
using WaltonForaging
using Test

@testset "WaltonForaging.jl" begin
    # Write your tests here.
end


text:
"""
Calculate ethanol critical properties
T --- temperature in Kelvin
p --- pressure in bar
ro --- density in mol/l
"""
function ethanolCritical()
	sub=newCritical()
	sub.T=514.71
	sub.p=6.268*10	#MPa to bar
	sub.ro=5.93		#mol/l #322 kg/m^3
	return sub
end

"""
Calculate ethanol properties at triple point
T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function ethanolTriple()
	sub=newTriple()
	sub.T=159.0
	sub.p=7.185*10^(-9)	#MPa to bar
	sub.lRo=19.731	#mol/l
	sub.vRo=5.43*10*(-10)	#mol/l
	return sub
end

"""
calculate ethanol properties at saturation line
A Fundamental Equation of State for Ethanol
J. A. Schroeder, S. G. Penoncello, and J. S. Schroeder
Citation: Journal of Physical and Chemical Reference Data 43, 043102 (2014); doi: 10.1063/1.4895394
T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function ethanolSaturate(T)
N=[9.00921, -23.1668, 30.9092, -16.5459, 3.64294]
k=[0.5, 0.8, 1.1, 1.5, 3.3]
Nv=[-1.75362, -10.5323, -37.6407, -129.762]
kv=[0.21, 1.1, 3.4, 10.0]
Np=[-8.94161, 1.61761, -51.1428, 53.1360]
kp=[1.0, 1.5, 3.4, 3.7]

sub=newSaturate()
crit=ethanolCritical()
if (T<crit.T)
	v=1.0-T/crit.T
	sub.p=crit.p*exp(crit.T/T*(Np[1]*v^kp[1]+Np[2]*v^kp[2]+Np[3]*v^kp[3]+Np[4]*v^kp[4]))
	sub.lRo=crit.ro*(1.0+N[1]*v^k[1]+N[2]*v^k[2]+N[3]*v^k[3]+N[4]*v^k[4]+N[5]*v^k[5])
	sub.vRo=crit.ro*exp(Nv[1]*v^kv[1]+Nv[2]*v^kv[2]+Nv[3]*v^kv[3]+Nv[4]*v^kv[4])
	sub.T=T
else
	sub.p=-1
	sub.lRo=-1
	sub.vRo=-1
	sub.T=-1
end
return sub
end

"""
calculate ethanol properties at single phase
A Fundamental Equation of State for Ethanol
J. A. Schroeder, S. G. Penoncello, and J. S. Schroeder
Citation: Journal of Physical and Chemical Reference Data 43, 043102 (2014); doi: 10.1063/1.4895394
Valid for up to 2800 bar and fromn 160 to 200 K

T --- temperature, Kelvin
p --- pressure, bar
ro --- density, mol/l
z ---
s --- entropy,
s0 --- ideal part of entropy,
ds --- nonideal part derived k_B, dimensionless
g --- Gibbs energy
g0 --- ideal part of Gibbs energy
dg --- nonideal part of Gibbs energy
f ---  Helmholtz free energy
f0 --- ideal part
df --- nonideal part, dimensionless
u --- internal energy
du --- nonideal part
h --- entalpy
dh --- nonideal part
w --- speed of sound,
cp ---
cv ---
"""
function ethanolSinglephase(T,ro)

sub=newSingle()
crit=ethanolCritical()
delta=ro/crit.ro
tau=crit.T/T
R=8.3144598
mm=0.04606844

nOne=6
nTwo=10
nThree=9
nFour=0

c=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.000, 1.000, 2.000, 1.000,
2.000, 1.000, 2.000, 1.000, 1.000,
1.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000]

d=[4.0, 1.0, 1.0, 2.0, 2.0,
3.0, 1.0, 1.0, 1.0, 3.0,
3.0, 2.0, 2.0, 6.0, 6.0,
8.0, 1.0, 1.0, 2.0, 3.0,
3.0, 2.0, 2.0, 2.0, 1.0]

t=[1.0, 1.04, 2.72, 1.174,
1.329, 0.195, 2.43, 1.274,
4.16, 3.3, 4.177, 2.5, 0.81,
2.02, 1.606, 0.86, 2.5, 3.72,
1.19, 3.25, 3.0, 2.0, 2.0, 1.0, 1.0]

n=[0.058200796, 0.94391227, -0.80941908,
0.55359038, -1.4269032, 0.13448717,
0.42671978, -1.1700261, -0.92405872,
0.34891808, -0.91327720, 0.022629481,
-0.15513423, 0.21055146, -0.21997690,
-0.0065857238, 0.75564749, 0.1069411,
-0.069533844, -0.24947395, 0.027177891,
-0.00090539530, -0.12310953, -0.089779710,
-0.39512601]

al=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.075, 0.463, 0.876, 1.108,
0.741, 4.032, 2.453, 2.300, 3.143]

bt=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.207, 0.090, 0.581, 0.947,
2.356, 27.010, 4.542, 1.287, 3.090]

gm=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.194, 1.986, 1.583, 0.756,
0.495, 1.002, 1.077, 1.493, 1.542]

ep=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.779, 0.805, 1.869, 0.694,
1.312, 2.054, 0.441, 0.793, 0.313]

phidelta=0.0
phir=0.0
phitau=0.0
phidelta2=0.0
phitau2=0.0
dphideltadtau=0.0

i=0
for ni=1:nOne
	i=i+1
	phidelta=phidelta+n[i]*d[i]*delta^(d[i]-1.0)*tau^t[i]
	phir=phir+n[i]*delta^d[i]*tau^t[i]
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1)
	phidelta2=phidelta2+n[i]*d[i]*(d[i]-1.0)*delta^(d[i]-2)*tau^t[i]
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2.0)
	dphideltadtau=dphideltadtau+n[i]*d[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)
end
for ni=1:nTwo
	i+=1
	phidelta=phidelta+n[i]*exp(-delta^c[i])*(delta^(d[i]-1.0)*tau^t[i]*(d[i]-c[i]*delta^c[i]))
	phir=phir+n[i]*delta^(d[i])*tau^t[i]*exp(-delta^c[i])
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1.0)*exp(-delta^c[i]);
	phidelta2=phidelta2+n[i]*exp(-delta^c[i])*(delta^(d[i]-2)*tau^(t[i])*((d[i]-c[i]*delta^c[i])*(d[i]-1.0-c[i]*delta^c[i])-c[i]^2*delta^c[i]))
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2)*exp(-delta^c[i])
	dphideltadtau=dphideltadtau+n[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)*(d[i]-c[i]*delta^c[i])*exp(-delta^c[i])
end
for pn=1:nThree
	i+=1
	phidelta=phidelta+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta-2.0*al[i]*(delta-ep[i]))
	phir=phir+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)
	phitau=phitau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(t[i]/tau-2*bt[i]*(tau-gm[i]))
	phidelta2=phidelta2+n[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(-2*al[i]*delta^d[i]+4*al[i]^2*delta^d[i]*(delta-ep[i])^2-4*d[i]*al[i]*delta^(d[i]-1)*(delta-ep[i])+d[i]*(d[i]-1)*delta^(d[i]-2))
	phitau2=phitau2+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*((t[i]/tau-2*bt[i]*(tau-gm[i]))^2-t[i]/tau^2-2*bt[i])
	dphideltadtau=dphideltadtau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta -2*al[i]*(delta-ep[i]))*(t[i]/tau-2*bt[i]*(tau-gm[i]))
end

sub.T=T
sub.ro=ro
sub.p=(1.0+phidelta*delta)*ro*R*T/100.0	#bar
sub.z=1.0+phidelta*delta	#
#ideal part
n0=[-12.7531, 9.39094, 3.43069, 2.14326, 5.09206, 6.60138, 5.70777]
gamma0=[0.0, 0.0, 0.0, 0.816771, 2.59175, 3.80408, 8.58736]

#entropy
phi0=log(delta)+n0[1]+n0[2]*tau+n0[3]*log(tau)+n0[4]*log(1.0-exp(-gamma0[4]*tau))+n0[5]*log(1.0-exp(-gamma0[5]*tau))+n0[6]*log(1.0-exp(-gamma0[6]*tau))+n0[7]*log(1.0-exp(-gamma0[7]*tau))
phi0tau=n0[2]+n0[3]/tau+n0[4]*gamma0[4]*(1.0/(1.0-exp(-gamma0[4]*tau))-1.0)+n0[5]*gamma0[5]*(1.0/(1.0-exp(-gamma0[5]*tau))-1.0)+n0[6]*gamma0[6]*(1.0/(1.0-exp(-gamma0[6]*tau))-1.0)+n0[7]*gamma0[7]*(1.0/(1.0-exp(-gamma0[7]*tau))-1.0)
phi0delta=1/delta
phi0delta2=-1/delta^2
phi0tau2=-n0[3]/tau^2-n0[4]*gamma0[4]^2*exp(-gamma0[4]*tau)*(1.0-exp(-gamma0[4]*tau))^(-2)-n0[5]*gamma0[5]^2*exp(-gamma0[5]*tau)*(1.0-exp(-gamma0[5]*tau))^(-2)-n0[6]*gamma0[6]^2*exp(-gamma0[6]*tau)*(1.0-exp(-gamma0[6]*tau))^(-2)-n0[7]*gamma0[7]^2*exp(-gamma0[7]*tau)*(1.0-exp(-gamma0[7]*tau))^(-2)
phi0taudelta=0.0

#Entropy
sub.s=(tau*(phi0tau+phitau)-phi0-phir)*R
sub.s0=(tau*(phi0tau)-phi0)*R
sub.ds=(tau*(phitau)-phir)	#devidet to R

#out.entrfrac=(tau*(phi0tau+phitau)-phi0-phir)/(tau*(phi0tau)-phi0)

dtaudtau=-tau/T
dphidt=dtaudtau*dphideltadtau

#sub.dpdt=(delta*cur_temp*dphidt+delta*phidelta+1)*R*T/100;
#phi0
#phir
#phidelta
#phidelta*delta
#(1+phi0+phir+phidelta*delta)
sub.g=(1+phi0+phir+phidelta*delta)*R*T	#R*cur_temp*
sub.g0=(1+phi0)*R*T
sub.dg=(phir+phidelta*delta)	#deriver RT
#
sub.f=(phi0+phir)*R*T
sub.f0=phi0*R*T
sub.df=phir

#Energy
sub.u=tau*(phi0tau+phitau)*R*T
sub.du=tau*(phitau)
#
sub.h=(1+tau*(phi0tau+phitau)+delta*phidelta)*R*T
sub.dh=1+tau*(phitau)+delta*phidelta;
#speed of sound
sub.w=sqrt(R/mm*T*(1+2*delta*phidelta+delta^2*phidelta2-(1+delta*phidelta-delta*tau*dphideltadtau)^2/(tau^2*(phi0tau2+phitau2))))	#0.018 перевод в килограмы
#
sub.cv=-tau^2*(phi0tau2+phitau2)*R
sub.cp=(-tau^2*(phi0tau2+phitau2)+(1+delta*phidelta-delta*tau*dphideltadtau)^2/(1+2*delta*phidelta+delta^2*phidelta2))*R
sub.dpdt=R/100.0*T*(1.0+2*delta*phidelta+delta^2*phidelta2)
#debug
# println("phi0 $(phi0)")
# println("phi0delta $(phi0delta)")
# println("phi0delta2 $(phi0delta2)")
# println("phi0tau $(phi0tau)")
# println("phi0tau2 $(phi0tau2)")
# println("phi0taudelta $(phi0taudelta)")
#
# println("phir $(phir)")
# println("phidelta $(phidelta)")
# println("phidelta2 $(phidelta2)")
# println("phitau $(phitau)")
# println("phitau2 $(phitau2)")
# println("dphideltadtau $(dphideltadtau)")
#enddebug
return sub
end

"""
calculates maximum of eos works
"""
function ethanolMax(T)
	sub = newViscMax()
	if(T < ethanolCritical().T)
		sub.p = 2800
		init = ethanolSaturate(T).lRo*1.2
	else
		sub.p = 2800
		init = ethanolCritical().ro * 1.2
	end
	sub.ro = find_zero(x -> ethanolSinglephase(T,x).p-sub.p, init)
	return sub
end

"""
uPa*s
Generalized SAFT-DFT/DMT Model for the Thermodynamic,
Interfacial, and Transport Properties of Associating Fluids:
Application for n-Alkanols
S. B. Kiselev and J. F. Ely
I. M. Abdulagatov and M. L. Huber
10.1021/ie050010e

"""
function ethanolVisc(T,ro)
	nu0=-1.03116 + 3.48379e-2*T -6.50264e-6*T^2
	b=[-19.572881, 219.73999, -1015.3226, 2471.0125, -3375.1717, 2491.6597, -787.26086, 14.085455, -0.34664158]
	t=[0.0, -0.25, -0.5, -0.75, -1.0, -1.25, -1.5, -2.5, -5.5]
	bb=0.0
	Tb=T/362.6
	for i=1:9
		bb+=b[i]*Tb^t[i]
	end
	br=bb*0.602*0.453^3

	a=[[0.0, 0.0, 0.0], [0.131194057, -0.382240694, 0.0], [-0.0805700894, 0.153811778, -0.110578307]]
	c=[23.7222995, -3.38264465, 12.7568864]
	Tr=513.9
	ror=5.991
	del=ro/ror
	tau=T/Tr
	del0=c[2]+c[3]*sqrt(tau)

#	println(" i $(3) j $(2) aij $(a[3][3])")
	dh=0
	for j=2:3
		for k=1:3
			dh+=a[j][k]*del^j/tau^(k-1)
		end
	end
	dh=1000*(dh+c[1]*del*(1/(del0-del)-1/del0))
#	println("nu0 $(nu0) br $(br) dh $(dh)")
	return nu0*(1.0+br*ro)+dh
end

"""
ethanol viscosity maximum
"""
function ethanolViscMax(T)
	sub = newViscMax()
	if(T < 480)
		sub.p = 1000
		init = ethanolSaturate(T).lRo*1.2
	else
		sub.p = 1000
		if(T < ethanolCritical().T)
			init = ethanolSaturate(T).lRo*1.2
		else
			init = ethanolCritical().ro * 1.2
		end
	end
	sub.ro = find_zero(x -> ethanolSinglephase(T,x).p-sub.p, init)
	return sub
end

"""
ethanol thermal conductivity, [m W m^-1 K^-1]
Reference Correlation of the Thermal Conductivity of Ethanol from the Triple Point to 600 K and up to 245 MPa
M. J. Assael, E. A. Sykioti, M. L. Huber and R. A. Perkins
10.1063/1.4797368
(Empirical critical enhancement)

"""
function ethanolTherm(T,ro)
	crit=ethanolCritical()
	Tr=T/crit.T
	ror=ro/crit.ro

	la0=(-2.09575 + 19.9045*Tr - 53.964*Tr^2 + 82.1223*Tr^3 - 1.98864*Tr^4 - 0.495513*Tr^5)/(0.17223 - 0.078273*Tr + Tr^2)
	B1=[2.67222e-2, 1.48279e-1, -1.30429e-1, 3.46232e-2, -2.44293e-3]
	B2=[1.77166e-2, -8.93088e-2, 6.84664e-2, -1.45702e-2, 8.09189e-4]
#	println("la0 $(la0)")
	dla=0
	for i=1:5
		dla+=(B1[i]+B2[i]*Tr)*ror^i
#		println("i $(B1[i]) dla $(B2[i])")
	end
	dlac=1.7e-3/(7.0e-2+Tr-1.0)*exp(-(1.7*(ror-1.0))^2)
#	println("dlac $(dlac)")
#more accuracy
#Rd=1.02
#v=0.63
#gam=1.239
#Gam=0.05885
#dz0=1.64296*10.0^-10
#qd=5.3*10.0^-10
#Tref=772.06
#kB=1.38064852*10.0^-23

#eos_crit=ethanolCritical()
#eos_cur=ethanolSinglephase(T,ro)

#println("$(eos_cur.dpdt) --- $(ethanolSinglephase(Tref,ro).dpdt)")
#println("$(eos_cur.dpdt-Tref/T*ethanolSinglephase(Tref,ro).dpdt)")
#println("$((v/gam))")

#if(eos_cur.dpdt-Tref/T*ethanolSinglephase(Tref,ro).dpdt>0)
#	dz=dz0*(eos_crit.p*ro/Gam/eos_crit.ro^2)^(v/gam)*(eos_cur.dpdt-Tref/T*ethanolSinglephase(Tref,ro).dpdt)^(v/gam)
#else
#	dz=0.0
#end
#if(qd*dz>1.7e-7)
#	Om0=2.0/pi*(1-exp(-1/(qd*dz)^-1+(qd*dz*eos_crit.ro/ro)^2/3))
#	Om=2/pi*((eos_cur.cp-eos_cur.cv)/eos_crit.cp*atan(qd*dz)+eos_cur.cv/eos_cur.cp*qd*dz)
#else
#	Om0=0.0
#	Om=0.0
#end

#dlac2=ro*eos_cur.cp*Rd*kB*T/6/pi/ethanolVisc(T,ro)/dz*(Om-Om0)
#println("dlac $(dlac) dlac2 $(dlac2)")

#
	return la0+dla*1000.0+dlac
end

"""

"""
function ethanolThermMax(T)
	sub = newViscMax()
	if(T < 480)
		sub.p = 2500
		init = ethanolSaturate(T).lRo*1.2
	else
		sub.p = 2500
		if(T < ethanolCritical().T)
			init = ethanolSaturate(T).lRo*1.2
		else
			init = ethanolCritical().ro * 1.2
		end
	end
	sub.ro = find_zero(x -> ethanolSinglephase(T,x).p-sub.p, init)
	return sub

end

"""
ethanol meltion line

"""
function ethanolIceBound(T)
	sub = newCritical()	#one point
	sub.p = 4369 * ((T/158.37)^2.6432 -1)
	init = ethanolTriple().lRo*1.02
	sub.ro = find_zero(x -> ethanolSinglephase(T,x).p-sub.p, init)

	return sub
end


text:
module ApproxOperator

import Base: +, -, *, getindex, setindex!, getproperty, setproperty!, length, push!, fill!, issubset, intersect
import InteractiveUtils: subtypes

abstract type AbstractNode end
abstract type AbstractElement{T} end
abstract type SpatialPartition end

include("node.jl")
include("snode.jl")
include("approximation.jl")
include("approximation_mf.jl")
include("integration.jl")
include("operation.jl")
include("operation_thin_plate.jl")
include("approximation_rk.jl")
include("import.jl")

export Node, Element, SNode, ReproducingKernel, getnₚ
export importmsh
export RegularGrid
export Operator, prescribe!, issubset, intersect
export set𝓖!
export set𝝭!, set∇𝝭!, set∇²𝝭!, set∇³𝝭!, set∇̃𝝭!, set∇̃²𝝭!, set∇∇̃²𝝭!, set∇̄𝝭!, set𝒏!, set∇𝑢!, get∇𝑢, get𝝐, set_memory_𝝭!

#debug
include("littletools.jl")
export get𝐴,cal𝗠!,cal𝗚!,get𝒙,get∇𝝭,get𝝭,checkIC, checkCC, checkConsistency

end