text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test
using Distributions
using PDMats

@testset "mvdist_density" begin
    mvt = @inferred MvTDist(1.5, PDMat([2.0 1.0; 1.0 3.0]))
    mvdd = @inferred MvDistDensity(mvt)

    @testset "mvdist_density" begin
        @test typeof(mvdd) <: AbstractDensity
        @test parent(mvdd) == mvt
        @test nparams(mvdd) == 2

        @test typeof(sampler(mvdd)) <: BATMvTDistSampler
    end

    @testset "unsafe_density_logval" begin
        BAT.unsafe_density_logval(mvdd, [0.0, 0.0]) ≈
            -2.64259602
        param = zeros(2)
        ec = @inferred BAT.exec_capabilities(BAT.unsafe_density_logval, mvdd, param)
        @test ec.nthreads == 0
        @test ec.threadsafe == true
        @test ec.nprocs == 0
        @test ec.remotesafe == true

        ret = Array{Float64}(2)
        param = [0.0 0.5; 0.0 -0.5]
        BAT.unsafe_density_logval!(ret, mvdd, param)
        @test ret ≈ [-2.64259602, -3.00960695]
        ec = @inferred BAT.exec_capabilities(BAT.unsafe_density_logval!, ret,
            mvdd, param)
        @test ec.nthreads == 0
        @test ec.threadsafe == true
        @test ec.nprocs == 0
        @test ec.remotesafe == true
    end
end


text:
push!(LOAD_PATH, pwd())
using PyPlot, Colors, ColorSchemes
using Statistics
using Formatting

include("$(pwd())/startup.jl")
include("$(pwd())/includes/DatasetFunctions.jl")

file1 = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP_13-07/results/_result.hdf5"
file2 = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP_14-07/1stFile/results/_result.hdf5"
file3 = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP/results/_result.hdf5"

measResult = ResultFileFunctions.loadResults(file1; useAveragesOnly = false, raw = true)	# to load ApiTOF data!!!
#measResult2 = ResultFileFunctions.loadResults(file2; useAveragesOnly = false, raw = true)	# to load ApiTOF data!!!
#measResult3 = ResultFileFunctions.loadResults(file3; useAveragesOnly = false, raw = true)	# to load ApiTOF data!!!

#measResult = ResultFileFunctions.joinResultsTime(measResult1, measResult2)
#measResult = ResultFileFunctions.joinResultsTime(measResult, measResult3)

steps = [	DateTime(2019,7,13,18,23) -600	1;
		# DateTime(2019,7,13,18,49) -500	2;
		# DateTime(2019,7,13,19,13) -400	3;
		# DateTime(2019,7,13,19,34) -300	4;
		DateTime(2019,7,13,19,52) -250	5;
		# DateTime(2019,7,13,20,20) -650	6;
		# DateTime(2019,7,13,20,38) -750	7;
		# DateTime(2019,7,13,21,10) -700	8;
		# DateTime(2019,7,13,21,30) -650	9;
		# DateTime(2019,7,13,21,42) -600	10;
		# DateTime(2019,7,13,21,54) -500	11;
		# DateTime(2019,7,13,22,08) -400	12;
		# DateTime(2019,7,13,22,22) -300	13;
#=
		DateTime(2019,7,14,14,33) -750	14;
		DateTime(2019,7,14,14,53) -700	15;
		DateTime(2019,7,14,15,37) -800	16;
		DateTime(2019,7,14,15,57) -800	17;
		DateTime(2019,7,14,16,05) -750	18;
		DateTime(2019,7,14,16,17) -700	19;
		DateTime(2019,7,14,20,15) -1000	20;
		DateTime(2019,7,14,20,27) -900	21;
		DateTime(2019,7,14,20,45) -800	22;
		DateTime(2019,7,15,12,40) -800	23;
		DateTime(2019,7,15,13,01) -850	24;
		DateTime(2019,7,15,13,51) -310	25;
		DateTime(2019,7,15,15,56) -750	26;
		DateTime(2019,7,15,16,09) -650	27;
		DateTime(2019,7,15,16,22) -600	28;
		DateTime(2019,7,15,16,36) -500	29;
		DateTime(2019,7,15,16,49) -400	30;
		DateTime(2019,7,15,17,03) -350	31;
		DateTime(2019,7,15,17,21) 0	32;
		DateTime(2019,7,15,17,39) -350	33;
		DateTime(2019,7,15,17,56) -400	34;
		DateTime(2019,7,15,18,10) -500	35;
		DateTime(2019,7,15,18,24) -600	36;
		DateTime(2019,7,15,18,37) -650	37;
		DateTime(2019,7,15,18,54) -750	38;
		DateTime(2019,7,16,10,45) -800	39;
		DateTime(2019,7,16,11,00) -900	40;
		DateTime(2019,7,16,11,18) -700	41;
		DateTime(2019,7,16,21,54) -700	42;
=#
			]

function Dp(m)
	Z = (9.026 .* m .^(-0.318) .-0.315) .*1e-4
	Dp = (Z ./2.2458e-22) .^(1 ./-1.9956)
	Dp_nm = Dp*1e9
	return Dp_nm
end

function tickName(X)
    D = Dp(X)
    ticks = []
    for d in D
	tick = format(d, precision = 2)
        ticks = append!(ticks, [tick])
    end
    return ticks
end

# changing the font type
rcParams = PyDict(matplotlib["rcParams"])
rcParams["font.size"] = 14
rcParams["font.sans-serif"] = "Helvetica"

fig, axs = subplots(2, 1, sharex = true, figsize = (8.27, 11.69))

for i in 1:size(steps)[1]
	ax = axs[i]
	a = (steps[i,1] .< measResult.Times .< steps[i,1]+Minute(15))
	maxa = DatasetFunctions.nanmax(DatasetFunctions.nanmean(measResult.Traces[a,:],1))
	mina = DatasetFunctions.nanmin(DatasetFunctions.nanmean(measResult.Traces[a,:],1))
	println(string("$(i) - max. = ", maxa))
	println(string("$(i) - min. = ", mina))
	if !(isnan(maxa))
		global h = PlotFunctions.massDefectPlot(measResult.MasslistMasses, measResult.MasslistCompositions ,
							mean(measResult.Traces[a,:], dims = 1),
							measResult.MasslistCompositions[6,:]./measResult.MasslistCompositions[1,:],
							"bCP, $(steps[i,1]) - DMA $(steps[i,2]) V", "O:C ratio", ax;
							dotSize = 5, dotbase = 5, maxMass = 1200, maxDefect = 0.6, minConc = 1.0e-2, maxConc = 5,
							sumformulas = false, negativeMD = false, colormap = "terrain", norm = 0)
		ax.annotate("DMA: $(steps[i,2]) V",
			xy=[1200;0.0],
			fontsize=18.0,
			ha="right",
			va="bottom")
	end
end

(axs[2]).set_xlabel("mass [amu]")
#gca() = axs[1]
(axs[1]).set_xlim(0,1200)
ax3 = (axs[1]).twiny()
ax1Ticks = (axs[1]).get_xticks()

ticklabels = tickName(ax1Ticks)
ax3.set_xticks(ax1Ticks)
ax3.set_xticklabels(ticklabels)
ax3.set_xlabel("estimated mobility diameter [nm]")

fig.subplots_adjust(right=0.80)
cbar_ax = fig.add_axes([0.85, 0.1, 0.03, 0.8])
cb = fig.colorbar(h, cax=cbar_ax)
cb["ax"]["set_ylabel"]("O:C ratio")

(axs[1]).grid("on")
(axs[2]).grid("on")
(axs[1]).set_ylabel("mass defect")
(axs[2]).set_ylabel("mass defect")

fig.subplots_adjust(hspace=0.08)
savefig("/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/plots/bCP_massDefect_O-to-C.png")

##########

volatilityColors = matcolors.ListedColormap(["#ddbcdd", "#d6d9d6", "#ffc4c4","#c6ffc6", "#b6b6ff", "#ffffff"])
volatilityBounds = [-2000, -8.3, -4.3, -0.3, 2.7, 6.7, 15]
volatilityNorm = matcolors.BoundaryNorm(volatilityBounds, volatilityColors.N)

fig, axs = subplots(2, 1, sharex = true, figsize = (8.27, 11.69))

for i in 1:size(steps)[1]
	ax = axs[i]
	a = (steps[i,1] .< measResult.Times .< steps[i,1]+Minute(15))
	maxa = DatasetFunctions.nanmax(DatasetFunctions.nanmean(measResult.Traces[a,:],1))
	mina = DatasetFunctions.nanmin(DatasetFunctions.nanmean(measResult.Traces[a,:],1))
	println(string("$(i) - max. = ", maxa))
	println(string("$(i) - min. = ", mina))
	if !(isnan(maxa))
		global h = PlotFunctions.massDefectPlotVolatility(measResult.MasslistMasses, measResult.MasslistCompositions ,
							mean(measResult.Traces[a,:], dims = 1),
							"bCP, $(steps[i,1]) - DMA $(steps[i,2]) V", ax;
							dotSize = 5, dotbase = 5, maxMass = 1200, maxDefect = 0.6, minConc = 1.0e-2, maxConc = 5,
							cmap = volatilityColors, norm = volatilityNorm, sumformulas = false, negativeMD = false)
		ax.annotate("DMA: $(steps[i,2]) V",
			xy=[1200;0.0],
			fontsize=18.0,
			ha="right",
			va="bottom")
	end
end

(axs[2]).set_xlabel("mass [amu]")
#gca() = axs[1]
(axs[1]).set_xlim(0,1200)
ax3 = (axs[1]).twiny()
ax1Ticks = (axs[1]).get_xticks()

ticklabels = tickName(ax1Ticks)
ax3.set_xticks(ax1Ticks)
ax3.set_xticklabels(ticklabels)
ax3.set_xlabel("estimated mobility diameter [nm]")

fig.subplots_adjust(right=0.80)
cbar_ax = fig.add_axes([0.85, 0.1, 0.03, 0.8])

cb = matcolorbar.ColorbarBase(cbar_ax, cmap=volatilityColors,
                                norm=volatilityNorm,
                                boundaries= volatilityBounds,
                                extend="both",
                                extendfrac="auto",
                                ticks=volatilityBounds,
                                spacing="uniform",
                                orientation="vertical")

#cb.set_label("Custom extension lengths, some other units")

#cb = fig.colorbar(h, cax=cbar_ax)
cb["ax"]["set_ylabel"](L"C^*_{300K}")

(axs[1]).grid("on")
(axs[2]).grid("on")
(axs[1]).set_ylabel("mass defect")
(axs[2]).set_ylabel("mass defect")

fig.subplots_adjust(hspace=0.08)
savefig("/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/plots/bCP_massDefect_volatility.png")

ext: jl


text:
mutable struct Sampling <: AbstractSampling
  parameters::Vector{Symbol}

  sensitivity::Real

  study_count::Int

  deck::Union{Void, Symbol}

  defaults::Vector{Real}

  kink_reactors::Vector{AbstractNullableReactor}
  wall_reactors::Vector{AbstractNullableReactor}

  cost_reactors::Vector{AbstractNullableReactor}
  W_M_reactors::Vector{AbstractNullableReactor}
end

function Sampling(parameters::Vector{Symbol}; sensitivity::Real=0.35, study_count::Int=100, deck::Union{Void, Symbol}=nothing, is_parallel::Bool=true, cur_kwargs...)
  cur_dict = merge!(Dict(), Dict(cur_kwargs))
  cur_dict[:constraint] = :beta
  cur_dict[:deck] = deck

  cur_reactor_types = [:kink, :wall, :cost,  :W_M]

  defaults = []
  ranges = OrderedDict()

  reactors_per_study = 4

  work_reactor = Reactor(symbols(:T_bar), cur_dict)

  for cur_parameter in parameters
    cur_med_value = getfield(work_reactor, cur_parameter)
    push!(defaults, cur_med_value)

    cur_min_value = cur_med_value * ( 1 - sensitivity )
    cur_max_value = cur_med_value * ( 1 + sensitivity )
    ranges[cur_parameter] = (cur_min_value, cur_max_value)
  end

  cur_sampling = Sampling(
    parameters, sensitivity, study_count, deck, defaults,
    [ fill!(Vector(study_count), nothing) for cur_index in 1:reactors_per_study]...
  )

  @assert !haskey(ranges, :T_bar)

  cur_constraints = SharedArray{Int}(study_count, reactors_per_study)

  cur_array = SharedArray{Float64}(study_count, reactors_per_study, 3+length(ranges))
  fill!(cur_array, NaN)

  cur_func = function (cur_index::Integer)
    tmp_dict = deepcopy(cur_dict)

    tmp_array = Vector{Float64}(length(ranges))

    for (cur_sub_index, (cur_key, cur_value)) in enumerate(ranges)
      if isa(cur_value, Tuple)
        if length(unique(cur_value)) == 1
          tmp_value = cur_value[1]
        else
          tmp_value = rand(Uniform(cur_value...))
        end
      elseif isa(cur_value, AbstractArray)
        tmp_value = pluck(cur_value)
      elseif isa(cur_value, Number)
        tmp_value = cur_value
      else
        error("Invalid sampling range for: $(cur_key)")
      end

      tmp_array[cur_sub_index] = tmp_value
      tmp_dict[cur_key] = tmp_value
    end

    tmp_reactor = Reactor(symbols(:T_bar), tmp_dict)

    tmp_reactor.is_good || return

    tmp_study = Study(:H, tmp_dict; num_points=1, verbose=false, is_parallel=false)

    for (cur_sub_index, cur_reactor_type) in enumerate(cur_reactor_types)
      tmp_sub_reactors = getfield(tmp_study, Symbol("$(cur_reactor_type)_reactors"))
      isempty(tmp_sub_reactors) && continue

      @assert length(tmp_sub_reactors) == 1
      tmp_sub_reactor = tmp_sub_reactors[1]

      tmp_sub_array = [
        tmp_sub_reactor.T_bar,
        tmp_sub_reactor.I_P,
        tmp_sub_reactor.eta_CD
      ]

      cur_array[cur_index, cur_sub_index, :] = [ tmp_sub_array..., tmp_array... ]
      cur_constraints[cur_index, cur_sub_index] =
        findfirst(tmp_constraint -> tmp_constraint == tmp_sub_reactor.constraint, collect(keys(secondary_limits)))
    end

    return
  end

  if is_parallel
    cur_progress = Progress(study_count)
    pmap(cur_func, cur_progress, shuffle(1:study_count))
  else
    map(cur_func, 1:study_count)
  end

  for cur_index in 1:study_count
    for (cur_sub_index, cur_reactor_type) in enumerate(cur_reactor_types)
      any(isnan, cur_array[cur_index,cur_sub_index,:]) && continue

      cur_T_bar, cur_I_P, cur_eta_CD = cur_array[cur_index,cur_sub_index,1:3]

      tmp_dict = deepcopy(cur_dict)

      for (cur_sub_sub_index, cur_range_name) in enumerate(keys(ranges))
        tmp_dict[cur_range_name] = cur_array[cur_index, cur_sub_index, cur_sub_sub_index+3]
      end

      tmp_dict[:constraint] = collect(keys(secondary_limits))[cur_constraints[cur_index,cur_sub_index]]

      work_reactor = Reactor(cur_T_bar, tmp_dict)
      work_reactor.I_P = cur_I_P
      work_reactor.eta_CD = cur_eta_CD

      work_reactors = getfield(cur_sampling, Symbol("$(cur_reactor_type)_reactors"))
      work_reactors[cur_index] = update!(work_reactor)
    end
  end

  cur_sampling
end


text:
const _ENCODE_AS_UUID_A = (
    Union{Nothing, Arc},
    Union{Nothing, Area},
    Union{Nothing, Bus},
    Union{Nothing, LoadZone},
    Union{Nothing, DynamicInjection},
    Vector{Service},
)

const _ENCODE_AS_UUID_B = (Arc, Area, Bus, LoadZone, DynamicInjection, Vector{Service})
@assert length(_ENCODE_AS_UUID_A) == length(_ENCODE_AS_UUID_B)

should_encode_as_uuid(val) = any(x -> val isa x, _ENCODE_AS_UUID_B)
should_encode_as_uuid(::Type{T}) where {T} = any(x -> T <: x, _ENCODE_AS_UUID_A)

function IS.serialize(component::T) where {T <: Component}
    @debug "serialize" component T
    data = Dict{String, Any}()
    for name in fieldnames(T)
        data[string(name)] = serialize_uuid_handling(getfield(component, name))
    end

    IS.add_serialization_metadata!(data, T)

    # This is a temporary workaround until these types are not parameterized.
    if T <: Reserve || T <: StaticReserveGroup
        data[IS.METADATA_KEY][IS.CONSTRUCT_WITH_PARAMETERS_KEY] = true
    end

    return data
end

"""
Serialize the value, encoding as UUIDs where necessary.
"""
function serialize_uuid_handling(val)
    if should_encode_as_uuid(val)
        if val isa Array
            value = IS.get_uuid.(val)
        elseif val === nothing
            value = nothing
        else
            value = IS.get_uuid(val)
        end
    else
        value = val
    end

    return serialize(value)
end

function IS.deserialize(::Type{T}, data::Dict, component_cache::Dict) where {T <: Component}
    @debug "deserialize Component" T data
    vals = Dict{Symbol, Any}()
    for (name, type) in zip(fieldnames(T), fieldtypes(T))
        val = data[string(name)]
        if val isa Dict && haskey(val, IS.METADATA_KEY)
            vals[name] = deserialize_uuid_handling(
                IS.get_type_from_serialization_metadata(IS.get_serialization_metadata(val)),
                val,
                component_cache,
            )
        else
            vals[name] = deserialize_uuid_handling(type, val, component_cache)
        end
    end

    type = IS.get_type_from_serialization_metadata(data[IS.METADATA_KEY])
    return type(; vals...)
end

function IS.deserialize(::Type{Device}, data::Any)
    error("This form of IS.deserialize is not supported for Devices")
end

"""
Deserialize the value, converting UUIDs to components where necessary.
"""
function deserialize_uuid_handling(field_type, val, component_cache)
    @debug "deserialize_uuid_handling" field_type val
    if val === nothing
        value = val
    elseif should_encode_as_uuid(field_type)
        if field_type <: Vector
            _vals = field_type()
            for _val in val
                uuid = deserialize(Base.UUID, _val)
                component = component_cache[uuid]
                push!(_vals, component)
            end
            value = _vals
        else
            uuid = deserialize(Base.UUID, val)
            component = component_cache[uuid]
            value = component
        end
    elseif field_type <: Component
        value = IS.deserialize(field_type, val, component_cache)
    elseif field_type <: Union{Nothing, Component}
        value = IS.deserialize(field_type.b, val, component_cache)
    elseif field_type <: InfrastructureSystemsType
        value = deserialize(field_type, val)
    elseif field_type <: Union{Nothing, InfrastructureSystemsType}
        value = deserialize(field_type.b, val)
    else
        value = deserialize(field_type, val)
    end

    return value
end


text:
using Pkg
Pkg.activate(".")
Pkg.status()
Pkg.offline()
Pkg.instantiate()

using Pidoh

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
using UtilitiesForMRI, LinearAlgebra, Test

# Cartesian domain
n = (256, 256, 256)
h = (abs(randn()), abs(randn()), abs(randn()))
X = spatial_sampling(Float64, n; h=h)

# Cartesian sampling in k-space
K = kspace_Cartesian_sampling(X; phase_encoding=(1,2))

# Random input
nt, nk = size(K)
d = randn(ComplexF64, nt, nk)

# Down-scaling
d_ = downscale(d, K; fact=1, flat=true)

# Test
n2 = div.(n,2)
n4 = div.(n,4)
cidx = div.(n,2).+1
@test reshape(reshape(d, n)[cidx[1]-n4[1]:cidx[1]+n4[1]-1,cidx[2]-n4[2]:cidx[2]+n4[2]-1,cidx[3]-n4[3]:cidx[3]+n4[3]-1], n2[1]*n2[2],n2[3]).*sqrt(prod(n4)/prod(n2)) ≈ d_ rtol=1e-6

# Consistency check
n = (128,128,128)
h = (abs(randn()), abs(randn()), abs(randn()))
X = spatial_sampling(Float64, n; h=h)
F = nfft_linop(X; phase_encoding=(1,2), tol=1e-6)
u = zeros(ComplexF64, n); u[65-10:65+10,65-10:65+10,65-10:65+10] .= 1
d = F*u
F_h = downscale(F; fact=1)
flat = false
coeff = 10
d_h = downscale(d, F.K; fact=1, flat=flat, coeff=coeff)
u_h = F_h'*d_h
F_2h = downscale(F; fact=2)
d_2h = downscale(d, F.K; fact=2, flat=flat, coeff=coeff)
u_2h = F_2h'*d_2h

# Downscaling images with anti-aliasing
u = zeros(ComplexF64, n); u[65-10:65+10,65-10:65+10,65-10:65+10] .= 1; u[65,65,65] = 10
u_h  = downscale(u, X; fact=1)
u_2h = downscale(u, X; fact=2)
u_h_ = upscale(u, X; fact=1)

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
"""
Sends all log records to the global logger on the master process. An additional
keyword argument, `_pid=myid()` is also passed.
"""
struct WorkerLogger <: AbstractLogger end
function WorkerLogger(logger::AbstractLogger)
    if myid() == 1
        return logger
    end

    return WorkerLogger()
end

function recieve_log_msg(pid, args...; kwargs...)
    handle_message(global_logger(), args...; _pid=pid, kwargs...)
    nothing
end

function handle_message(logger::WorkerLogger, args...; kwargs...)
    remotecall_fetch(EnhancedLogging.recieve_log_msg, 1, myid(), args...; kwargs...)
    nothing
end

shouldlog(::WorkerLogger, args...) = true
min_enabled_level(::WorkerLogger) = Base.CoreLogging.BelowMinLevel
catch_exceptions(::WorkerLogger) = true


text:
using AcuteML: UN, parsehtml, root, nextelement, nodecontent
using DataFrames: DataFrame, groupby
using JLD2: jldsave
using Pkg.Artifacts:
    artifact_hash, artifact_exists, create_artifact, archive_artifact, bind_artifact!
using Pseudopotentials:
    CoreHole, ExchangeCorrelationFunctional, ValenceCoreState, Pseudization

using UnifiedPseudopotentialFormat: UPFFile, analyzename
using UnifiedPseudopotentialFormat.PSlibrary: ELEMENTS, list_elements

const ARTIFACT_TOML = joinpath(dirname(@__DIR__), "Artifacts.toml")
const LIBRARY_URL_BASE = "http://pseudopotentials.quantum-espresso.org/legacy_tables/ps-library/"
const UPF_URL_BASE = "http://pseudopotentials.quantum-espresso.org"
const DATABASE_URL_BASE = "https://github.com/MineralsCloud/PseudopotentialArtifacts/raw/main/pslibrary/"

function getrawdata(element)
    url = LIBRARY_URL_BASE * lowercase(element)
    path = download(url)
    str = read(path, String)
    doc = parsehtml(str)
    primates = root(doc)
    anchors = findall("//table//a", primates)
    return map(anchors) do anchor
        (
            name = strip(nodecontent(anchor)),
            src = UPF_URL_BASE * anchor["href"],
            metadata = nodecontent(nextelement(anchor)),
        )
    end
end

function makedb(element::String)
    database = DataFrame(
        element = String[],
        rel = Bool[],
        corehole = UN{CoreHole}[],
        xc = UN{ExchangeCorrelationFunctional}[],
        cv = UN{Vector{<:ValenceCoreState}}[],
        pseudization = UN{Pseudization}[],
        src = String[],
        name = String[],
    )
    for meta in getrawdata(element)
        info = analyzename(UPFFile(meta.name))
        push!(
            database,
            [
                (
                    getindex(info, i) for i in (
                        :element,
                        :fullrelativistic,
                        :corehole,
                        :xc,
                        :valencecore,
                        :pseudization,
                    )
                )...,
                meta.src,
                meta.name,
            ],
        )
    end
    return database
end
makedb(i::Integer) = makedb(ELEMENTS[i])
function makedb()
    data = mapreduce(makedb, append!, ELEMENTS)
    return groupby(data, :element)
end

function savedb(file, element)
    database = makedb(element)
    return jldsave(file; database)
end
function savedb(file)
    database = makedb()
    return jldsave(file; database)
end

function makeartifact(element::AbstractString)
    dir_hash = artifact_hash(element, ARTIFACT_TOML)
    if dir_hash === nothing || !artifact_exists(dir_hash)
        dir_hash = create_artifact() do artifact_dir
            savedb(joinpath(artifact_dir, "$element.jld2"), element)
        end
        tar_hash = archive_artifact(dir_hash, "$element.tar.gz")
        bind_artifact!(
            ARTIFACT_TOML,
            element,
            dir_hash;
            download_info = [(DATABASE_URL_BASE * "$element.tar.gz", tar_hash)],
            lazy = true,
        )
    end
end
makeartifact(i::Integer) = makeartifact(ELEMENTS[i])
function makeartifact()
    dir_hash = artifact_hash("all", ARTIFACT_TOML)
    if dir_hash === nothing || !artifact_exists(dir_hash)
        dir_hash = create_artifact() do artifact_dir
            savedb(joinpath(artifact_dir, "all.jld2"))
        end
        tar_hash = archive_artifact(dir_hash, "all.tar.gz")
        bind_artifact!(
            ARTIFACT_TOML,
            "all",
            dir_hash;
            download_info = [(DATABASE_URL_BASE * "all.tar.gz", tar_hash)],
            lazy = true,
        )
    end
end


text:
#!/usr/bin/env julia

using SnaggedNecklace

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table! s begin
    "--mass"
        metavar = "M"
        help = "mass (g/mol)"
        arg_type = Float64
        required = true
    "--omega"
        metavar = "W"
        help = "angular frequency (1/ps)"
        arg_type = Float64
        required = true
    "--steps-equil"
        metavar = "N"
        help = "number of equilibration steps (default: 10% of production)"
        arg_type = Int
    "--steps-production"
        metavar = "N"
        help = "number of production steps"
        arg_type = Int
        required = true
    "--temp"
        metavar = "T"
        help = "temperature"
        arg_type = Float64
    "--temp-P"
        metavar = "TP"
        help = "temperature multiplied by P"
        arg_type = Float64
    "-P"
        metavar = "P"
        help = "number of beads"
        arg_type = Int
        required = true
    "--xi"
        metavar = "X"
        help = "value of constrained coordinate"
        arg_type = Float64
        required = true
    "--output-file"
        metavar = "PATH"
        help = "path to main output file"
        arg_type = String
    "--hist-file-template"
        metavar = "PATH"
        help = "template for paths to histogram output files"
        arg_type = String
end
c = parse_args(ARGS, s, as_symbols=true)

const mass = c[:mass]
const omega = c[:omega]

if !isnothing(c[:steps_equil])
    const steps_equil = c[:steps_equil]
else
    const steps_equil = div(c[:steps_production], 10)
end

const steps_production = c[:steps_production]

if isnothing(c[:temp]) == isnothing(c[:temp_P])
    error("One of --temp or --temp-P must be specified.")
elseif !isnothing(c[:temp])
    const temperature = c[:temp]
else
    const temperature = c[:temp_P] / c[:P]
end

const P = c[:P]
const xi_init = c[:xi]

const output_file = c[:output_file]
const hist_file_template = c[:hist_file_template]

pot(q) = 0.5 * mass * omega^2 * q^2

force(q) = -mass * omega^2 * q

const step_size_1 = StepSize1D1()
const step_size_j = StepSize1Dj(0.01)

sys = PimcSystem(pot, force, mass, temperature)
ip = PimcParameters(steps_equil, steps_production, step_size_1, step_size_j)
st = State1D(P, xi_init)

((mean1, err1, acsteps1), (mean2, err2, acsteps2), accept_ratio_1,
    accept_ratio_j, hists) = pimc(sys, ip, st)

if !isnothing(output_file)
    open(output_file, "w") do f
        print(f, "$(mean1) $(err1) $(acsteps1) ")
        print(f, "$(mean2) $(err2) $(acsteps2) ")
        println(f, "$(accept_ratio_1) $(accept_ratio_j)")
    end
end

if !isnothing(hist_file_template)
    for (j, (xis, ws)) in enumerate(hists)
        hist_file = replace(hist_file_template, "{}" => j)
        open(hist_file, "w") do f
            for (xi, w) in zip(xis, ws)
                println(f, "$(xi) $(w)")
            end
        end
    end
end


text:
using Documenter, TableTraitsUtils

makedocs(
	modules=[TableTraitsUtils],
	sitename="TableTraitsUtils.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo="github.com/queryverse/TableTraitsUtils.jl.git"
)


text:
struct SparseKernelData{K<:KernelData} <: KernelData
    Kuu::K
    Kux1::K
    Kux2::K
end

function SparseKernelData(k::Kernel, Xu::AbstractMatrix, X1::AbstractMatrix, X2::AbstractMatrix)
    Kuu = KernelData(k, Xu, Xu)
    Kux1 = KernelData(k, Xu, X1)
    if X1 == X2
        return SparseKernelData(Kuu, Kux1, Kux1)
    else
        Kux2 = KernelData(k, Xu, X2)
        return SparseKernelData(Kuu, Kux1, Kux2)
    end
end

function KernelData(k::Kernel, X1::AbstractMatrix, X2::AbstractMatrix, covstrat::SparseStrategy)
    Xu = covstrat.inducing
    return SparseKernelData(k, Xu, X1, X2)
end


text:
###############################################
#FEniCS tutorial demo program: Convection-diffusion-reaction for a system
#describing the concentration of three species A, B, C undergoing a simple
#first-order reaction A + B --> C with first-order decay of C. The velocity
#is given by the flow field w from the demo navier_stokes_cylinder.py.
#  u_1' + w . nabla(u_1) - div(eps*grad(u_1)) = f_1 - K*u_1*u_2
#  u_2' + w . nabla(u_2) - div(eps*grad(u_2)) = f_2 - K*u_1*u_2
#  u_3' + w . nabla(u_3) - div(eps*grad(u_3)) = f_3 + K*u_1*u_2 - K*u_3
###############################################

module ft09

using FenicsPy

T = 5.0            # final time
num_steps = 500    # number of time steps
Δt = T / num_steps # time step size
eps = 0.01         # diffusion coefficient
K = 10.0           # reaction rate

# Read mesh from file
mesh = Mesh("navier_stokes_cylinder/cylinder.xml.gz")

# Define function space for velocity
W = VectorFunctionSpace(mesh, "P", 2)

# Define function space for system of concentrations
P1 = FiniteElement("P", triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2, v_3 = TestFunctions(V)

# Define functions for velocity and concentrations
w = FeFunction(W)
u = FeFunction(V)
u_n = FeFunction(V)

# Split system functions to access components
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)

# Define source terms
f_1 = Expression("pow(x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0",
                 degree=1)
f_2 = Expression("pow(x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0",
                 degree=1)
f_3 = Constant(0)

# Define expressions used in variational forms
K = Constant(K)
eps = Constant(eps)

# Define variational problem
F = ((u_1 - u_n1) / Δt)*v_1*dx + dot(w, grad(u_1))*v_1*dx +
  eps*dot(grad(u_1), grad(v_1))*dx + K*u_1*u_2*v_1*dx  +
  ((u_2 - u_n2) / Δt)*v_2*dx + dot(w, grad(u_2))*v_2*dx +
  eps*dot(grad(u_2), grad(v_2))*dx + K*u_1*u_2*v_2*dx  +
  ((u_3 - u_n3) / Δt)*v_3*dx + dot(w, grad(u_3))*v_3*dx +
  eps*dot(grad(u_3), grad(v_3))*dx - K*u_1*u_2*v_3*dx + K*u_3*v_3*dx -
  f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx

# Create time series for reading velocity data
timeseries_w = TimeSeries("navier_stokes_cylinder/velocity_series")

# Create VTK files for visualization output
vtkfile_u_1 = File("reaction_system/u_1.pvd")
vtkfile_u_2 = File("reaction_system/u_2.pvd")
vtkfile_u_3 = File("reaction_system/u_3.pvd")

# Create progress bar
#progress = Progress("Time-stepping")
#set_log_level(PROGRESS)

# Time-stepping
t = 0
for i = 1:num_steps

    # Update current time
    global t += Δt

    # Read velocity from file
    timeseries_w.retrieve(w.vector(), t)

    # Solve variational problem for time step
    solve(F == 0, u)

    # Save solution to file (VTK)
    _u_1, _u_2, _u_3 = u.split()
    vtkfile_u_1 << (_u_1, t)
    vtkfile_u_2 << (_u_2, t)
    vtkfile_u_3 << (_u_3, t)

    # Update previous solution
    u_n.assign(u)

    # Update progress bar
    #progress.update(t / T)
end

end # module ft09


text:
# example copied from th0br0 with modifications
using CmdStan

y = collect(1:100);
x = [3:3:300 5:5:500];

model = """
data {
  int<lower = 0> N;
  int<lower = 0> k;
  matrix[N, k] x;
  real<lower=0> y[N];
}
parameters {
  vector[k] beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(x * beta, sigma);
}
""";

# beta parameters are accessed using beta.[1-9+] syntax
data = Dict("N"=>100, "k" => 2, "x" => x, "y" => y);
stanmodel = Stanmodel(monitors = ["beta.1", "beta.2", "sigma"], model=model,
          output_format=:mcmcchains);

rc, sim, cnames = stan(stanmodel, data; diagnostics = false);

cnames |> display
# 3-element Array{String,1}: "beta.1", "beta.2" ,"sigma"
println()

describe(sim)


text:
#############################################################################
# modular.jl
# Handles modular functions of a given vector x.
#############################################################################

export modular
export lovasz
export sign, monotonicity, modularity, evaluate

mutable struct ModularAtom <: SubmodFunc   # TODO: extend to sets other than plain set varables
  head::Symbol
  id_hash::UInt64
  children::Tuple{AbstractExpr, CombiSet}
  size::Tuple{Int, Int}
  param::AbstractExpr
  setvariables::Array{CombiSet}

  function ModularAtom(w::AbstractExpr, S::CombiSet)
    if w.size != (S.cardinality, 1)
      error("Cannot define a modular function when the vector's size is different from the (card(baseset), 1).")
    else
      children = (w, S)
      setvariables = get_sv(S)
      return new(:modular, hash(children), children, (1, 1), w, setvariables)
    end
  end
end

modular(w::Variable, S::CombiSet) = ModularAtom(w, S)

modular(w::Values, S::CombiSet) = ModularAtom(Constant(w), S)

function lovasz(F::ModularAtom, x::Variable)
  if length(F.setvariables) == 1
    if evaluate(F, [])[1] != 0
      error("A combinatorial function should be 0 at the empty set to derive its Lovasz extension.")
    else
      if x.size[1] == F.setvariables[1].cardinality
        if isa(F.children[2], SetVariable)
          return vecdot(F.children[1], x)
        else
          return LovaszExtAtom(F, x)
        end
      else
        error("The size of the continuous variable should be the same as the baseset of the combinatorial variable of the combinatorial function.")
      end
    end
  else
    error("A combinatorial function should be sigle-variant to obtain its Lovasz extension.")
  end
end

function lovasz(F::ModularAtom, x::AbsAtom)
  if length(F.setvariables) == 1
    if length(x.children) == 1 && isa(x.children[1], Variable)
      if evaluate(F, [])[1] != 0
        error("A combinatorial function should be 0 at the empty set to derive its Lovasz extension.")
      else
        var = get_v(x)
        if length(var) != 1
          error("Functions with other than one variable are not supported.")
        else
          if var[1].size[1] == F.setvariables[1].cardinality
            if isa(F.children[2], SetVariable)
              return vecdot(F.children[1], x)
            else
              return LovaszExtAbsAtom(F, x)
            end
          else
            error("The size of the continuous variable should be the same as the baseset of the combinatorial variable of the combinatorial function.")
          end
        end
      end
    else
      error("Only able to define Lovasz extesions on absolute values of a variable.")
    end
  else
    error("A combinatorial function should be sigle-variant to obtain its Lovasz extension.")
  end
end

function sign(F::ModularAtom)
  return sign(F.param)
end

function monotonicity(F::ModularAtom)
  return (sign(F.children[2]) * Nondecreasing(), sign(F.children[1]) * Nondecreasing())
end

function modularity(F::ModularAtom)
  return Modularity()
end

function evaluate(F::ModularAtom) # TODO only allow scalar input
  summ = 0.0
  param = evaluate(F.param)
  elements = []
  elements = vcat(elements, get_elements(F.children[2]))
  elements = Set(elements)
  for i in elements
    summ += param[i]
  end
  # Compatibility with Convex.jl
  return summ
end


text:
# Autogenerated wrapper script for Wigxjpf_jll for x86_64-w64-mingw32
export wigxjpf

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Wigxjpf")
JLLWrappers.@declare_library_product(wigxjpf, "libwigxjpf_shared.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        wigxjpf,
        "bin\\libwigxjpf_shared.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Fidelity in quantum information theory
# This example is inspired from a lecture of John Watrous in the [course on Theory of Quantum Information](https://cs.uwaterloo.ca/~watrous/CS766/LectureNotes/08.pdf).
#
# The Fidelity between two Hermitian semidefinite matrices P and Q is defined as:
#
# $$
# F(P, Q) = \|P^{1/2}Q^{1/2}\|_{\text{tr}} = \max_U \mathrm{tr}(P^{1/2}U Q^{1/2})
# $$
#
# where the trace norm $\|\cdot\|_{\text{tr}}$ is the sum of the singular values, and the maximization goes over the set of all unitary matrices U. This quantity can be expressed as the optimal value of the following complex-valued SDP:
#
# $$
# \begin{array}{ll}
#   \mbox{maximize} &  \frac{1}{2}\text{tr}(Z+Z^\dagger) \\
#   \mbox{subject to} &\\
#   & \left[\begin{array}{cc}P&Z\\{Z}^{\dagger}&Q\end{array}\right] \succeq 0\\
#   & Z \in \mathbf {C}^{n \times n}\\
# \end{array}
# $$
#

using Convex, SCS, LinearAlgebra
if VERSION < v"1.2.0-DEV.0"
     LinearAlgebra.diagm(v::AbstractVector) = diagm(0 => v)
end

n = 20
P = randn(n,n) + im*randn(n,n)
P = P*P'
Q = randn(n,n) + im*randn(n,n)
Q = Q*Q'
Z = ComplexVariable(n,n)
objective = 0.5*real(tr(Z+Z'))
constraint = [P Z;Z' Q] ⪰ 0
problem = maximize(objective,constraint)
solve!(problem, SCSSolver(verbose=0))
computed_fidelity = evaluate(objective)

#-

## Verify that computer fidelity is equal to actual fidelity
P1,P2 = eigen(P)
sqP = P2 * diagm([p1^0.5 for p1 in P1]) * P2'
Q1,Q2 = eigen(Q)
sqQ = Q2 * diagm([q1^0.5 for q1 in Q1]) * Q2'

#-

actual_fidelity = sum(svdvals(sqP * sqQ))

# We can see that the actual fidelity value is very close the computed fidelity value.

max_stars_repo_path: docs/examples_literate/optimization_with_complex_variables/Fidelity in Quantum Information Theory.jl


text:
mutable struct PlateAssessState{U,V}
    score::Float64
    is_empty::Bool
    args::U
    nodes::V
end

function assess_kernel!(gen::Plate{T,U}, key::Int, state::PlateAssessState) where {T,U}
    kernel_args = get_args_for_key(state.args, key)
    node = state.nodes[key]
    subtrace::U = assess(gen.kernel, kernel_args, node)
    state.is_empty = state.is_empty && !has_choices(subtrace)
    call::CallRecord{T} = get_call_record(subtrace)
    state.score += call.score
    (subtrace, call.retval::T)
end

function check_constraint_key!(key, len::Int)
    error("Constraints under key $key not part of trace")
end

function check_constraint_key!(key::Int, len::Int)
    if key < 1 || key > len
        error("Constraints under key $key not part of trace")
    end
end

function assess(gen::Plate{T,U}, args, constraints) where {T,U}
    len = length(args[1])

    # collect constraints, indexed by key
    nodes = Vector{Any}(len)
    for key=1:len
        nodes[key] = get_internal_node(constraints, key)
    end

    # check no other constraints.. (can be skipped if we're trying to be less strict)
    for (key, _) in get_internal_nodes(constraints)
        check_constraint!(key, len)
    end
    for (key, _) in get_leaf_nodes(constraints)
        check_constraint!(key, len)
    end

    # collect initial state
    state = PlateAssessState(0., true, args, nodes)
    subtraces = Vector{U}(len)
    retvals = Vector{T}(len)

    # visit each application and assess it
    for key=1:len
        (subtraces[key], retvals[key]) = assess_kernel!(gen, key, state)
    end

    call = CallRecord{PersistentVector{T}}(state.score, PersistentVector{T}(retvals), args, UnknownChange())
    VectorTrace{T,U}(PersistentVector{U}(subtraces), call, state.is_empty)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
@testset "dNdS" begin
    @testset "NG86" begin
        codonsA = [kmer"ATG",
                   kmer"AAA",
                   kmer"CCC",
                   kmer"GGG",
                   kmer"TTT",
                   kmer"GGG",
                   kmer"TAA"]

        codonsB = [kmer"ATG",
                   kmer"AAA",
                   kmer"CGC",
                   kmer"GGC",
                   kmer"TAC",
                   kmer"GGG",
                   kmer"TAA"]

        @testset "computation" begin
            @testset "expected" begin
                #all_codons = convert(Vector{Kmer{DNA,3}}, collect(UInt64(0):UInt64(63)))
                n_ans = [3.0, 2.666, 2.0, 2.0, 2.666, 2.0, 3.000]
                s_ans = [0.0, 0.333, 1.0, 1.0, 0.333, 1.0, 0.000]
                for i in 1:endof(codonsA)
                    cdn = codonsA[i]
                    @test S_N_NG86(cdn, ncbi_trans_table[1])[1] ≈ s_ans[i] atol=0.001
                    @test S_N_NG86(cdn, ncbi_trans_table[1])[2] ≈ n_ans[i] atol=0.001
                end
            end

            @testset "observed" begin
                function testobserved(a, b, ans)
                    @test DS_DN_NG86(a, b, ncbi_trans_table[1])[1] ≈ ans[1] atol=0.001
                    @test DS_DN_NG86(a, b, ncbi_trans_table[1])[2] ≈ ans[2] atol=0.001
                end
                answers = [(0.0, 0.0),
                           (0.0, 0.0),
                           (0.0, 1.0),
                           (1.0, 0.0),
                           (1.0, 1.0),
                           (0.0, 0.0),
                           (0.0, 0.0)]
                for i in 1:length(answers)
                    testobserved(codonsA[i], codonsB[i], answers[i])
                end
                testobserved(kmer"TTT", kmer"GAC", (1.0, 2.0))
            end
        end

        @testset "lookups" begin
            firstcdn = UInt64(0)
            lastcdn = UInt64(63)

            for i in firstcdn:lastcdn
                cdn = DNACodon(i)
                @test S_N_NG86(cdn, ncbi_trans_table[1])[1] == NaturalSelection.make_S_N_NG86_table(BioSequences.standard_genetic_code)[cdn][1]
                @test S_N_NG86(cdn, ncbi_trans_table[1])[2] == NaturalSelection.make_S_N_NG86_table(BioSequences.standard_genetic_code)[cdn][2]
            end

            for i in firstcdn:lastcdn
                for j in (i + 1):lastcdn
                    a = DNACodon(i)
                    b = DNACodon(j)
                    @test DS_DN_NG86(a, b, ncbi_trans_table[1])[1] == NaturalSelection.make_DS_DN_NG86_table(BioSequences.standard_genetic_code)[a, b][1]
                    @test DS_DN_NG86(a, b, ncbi_trans_table[1])[2] == NaturalSelection.make_DS_DN_NG86_table(BioSequences.standard_genetic_code)[a, b][2]
                end
            end
        end

        @testset "dN/dS" begin
            @test dNdS_NG86(codonsA[1:end-1], codonsB[1:end-1])[1] ≈ 0.154 atol=0.001
            @test dNdS_NG86(dna"ATGAAACCCGGGTTTGGG", dna"ATGAAACGCGGCTACGGG")[1] ≈ 0.154 atol=0.001
            @test dNdS_NG86(codonsA[1:end-1], codonsB[1:end-1])[2] ≈ 0.974 atol=0.001
            @test dNdS_NG86(dna"ATGAAACCCGGGTTTGGG", dna"ATGAAACGCGGCTACGGG")[2] ≈ 0.974 atol=0.001
        end

    end
end


text:
module MultiFieldFESpacesWithLinearConstraintsTests

using Gridap.Algebra
using Gridap.Arrays
using Gridap.Fields
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.TensorValues
using Gridap.MultiField
using Gridap.CellData
using Test

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet_1",[1,2,5])
add_tag_from_tags!(labels,"neumann_1",[6,7,8])
add_tag_from_tags!(labels,"dirichlet_2",[1,3,7])
add_tag_from_tags!(labels,"neumann_2",[5,6,8])

trian = Triangulation(model)
btrian1 = BoundaryTriangulation(model,"neumann_1")
btrian2 = BoundaryTriangulation(model,"neumann_2")
strian = SkeletonTriangulation(model)

quad = CellQuadrature(trian,2)
bquad1 = CellQuadrature(btrian1,2)
bquad2 = CellQuadrature(btrian2,2)
squad = CellQuadrature(strian,2)

bn1 = get_normal_vector(btrian1)
bn2 = get_normal_vector(btrian2)

V1 = FESpace(
  model=model,valuetype=Float64,reffe=:Lagrangian,order=1,
  conformity=:H1,dirichlet_tags="dirichlet_1")

V2 = FESpace(
  model=model,valuetype=Float64,reffe=:Lagrangian,order=1,
  conformity=:H1,dirichlet_tags="dirichlet_2")

sDOF_to_dof = [1,5,-2]
sDOF_to_dofs = Table([[-1,4],[4,6],[-1,-3]])
sDOF_to_coeffs = Table([[0.5,0.5],[0.5,0.5],[0.5,0.5]])
V1c = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V1)

sDOF_to_dof = [4,-2]
sDOF_to_dofs = Table([[2,6],[-1,-3]])
sDOF_to_coeffs = Table([[0.5,0.5],[0.5,0.5]])
V2c = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V2)

u1(x) = x[1] + 2*x[2]
u2(x) = 4*x[1] + x[2]
f1(x) = -Δ(u1)(x) + u2(x)
#f2(x) = -Δ(u2)(x)

U1c = TrialFESpace(V1c,u1)
U2c = TrialFESpace(V2c,u2)

U = MultiFieldFESpace([U1c,U2c])
V = MultiFieldFESpace([V1c,V2c])
test_fe_space(U)
test_fe_space(V)
@test has_constraints(U)
@test has_constraints(V)

a(u,v) = ∇(v)⋅∇(u)
b_Γ(v,u,n_Γ) = v*(n_Γ⋅∇(u))

function A(u,v)
  u1,u2 = u
  v1,v2 = v
  a(v1,u1) + v1*u2 + a(v2,u2)
end

function B(v)
  v1,v2 = v
  v1*f1
end

function B1_Γ(v)
  v1,v2 = v
  b_Γ(v1,u1,bn1)
end

function B2_Γ(v)
  v1,v2 = v
  b_Γ(v2,u2,bn2)
end

# This is only to stress skeleton machinery
# since shape functions are continuous
function As(u,v)
  u1,u2 = u
  v1,v2 = v
  jump(v1)*jump(u2) + jump(v2)*jump(u2)
end

function Bs(v)
  v1,v2 = v
  jump(v2)
end

t_Ω = AffineFETerm(A,B,trian,quad)
t1_Γ = FESource(B1_Γ,btrian1,bquad1)
t2_Γ = FESource(B2_Γ,btrian2,bquad2)
t_s = AffineFETerm(As,Bs,strian,squad)

op = AffineFEOperator(U,V,t_Ω,t1_Γ,t2_Γ,t_s)
uh = solve(op)

uh1, uh2 = uh

e1 = u1 - uh1
e2 = u2 - uh2

#using Gridap.Visualization
#v1 = FEFunction(V1,rand(num_free_dofs(V1)))
#v2 = FEFunction(V2,rand(num_free_dofs(V2)))
#v1c = FEFunction(V1c,rand(num_free_dofs(V1c)))
#v2c = FEFunction(V2c,rand(num_free_dofs(V2c)))
#writevtk(trian,"trian",nsubcells=10,cellfields=["v1"=>v1,"v2"=>v2,"v1c"=>v1c,"v2c"=>v2c])

#using Gridap.Visualization
#writevtk(trian,"trian",nsubcells=10,cellfields=["uh1"=>uh1,"uh2"=>uh2,"e1"=>e1,"e2"=>e2])

l2(u) = u*u
e1_l2 = sqrt(sum(integrate(l2(e1),trian,quad)))
e2_l2 = sqrt(sum(integrate(l2(e2),trian,quad)))

tol = 1.0e-9
@test e1_l2 < tol
@test e2_l2 < tol

end # module


text:
export getSensTMatVec

function getSensTMatVec(v::Vector,m::Vector,pFor::FWIparam)

    # extract pointers
    M  				= pFor.Mesh
    omega 			= pFor.omega
    gamma 			= pFor.gamma
    Q     			= pFor.Sources
    P     			= pFor.Receivers
    Ainv			= pFor.ForwardSolver

	batchSize 		= pFor.forwardSolveBatchSize;
	select  		= pFor.sourceSelection;

	nsrc = size(Q,2); nrec = size(P,2);

	if length(select) > 0
		nsrc = length(select);
	end

	if batchSize > nsrc
		batchSize = nsrc;
	end

	numBatches 	= ceil(Int64,nsrc/batchSize);
	An2cc = getNodalAverageMatrix(M);
	m_nodal = An2cc'*m;
	gamma_nodal = An2cc'*gamma;
	n_nodal = prod(M.n.+1);
	H = spzeros(ComplexF64,n_nodal,n_nodal);
	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		H = GetHelmholtzOperator(M,m_nodal,omega, gamma_nodal, true,useSommerfeldBC);
		Ainv.helmParam = HelmholtzParam(M,gamma_nodal,m_nodal,omega,true,useSommerfeldBC);
		H = H + GetHelmholtzShiftOP(m_nodal, omega,Ainv.shift[1]);
		H = sparse(H');
		# H is actually shifted laplacian now...
	elseif isa(Ainv,JuliaSolver)
		H = GetHelmholtzOperator(M,m_nodal,omega, gamma_nodal, true,useSommerfeldBC);
	end

	JTv = zeros(Float64,prod(M.n))
	Vdatashape = reshape(v,nrec,nsrc);

	for k_batch = 1:numBatches
		batchIdxs = (k_batch-1)*batchSize + 1 : min(k_batch*batchSize,nsrc);
		V = P*Vdatashape[:,batchIdxs];
		V,Ainv = solveLinearSystem(H,V,Ainv,1);
		if pFor.useFilesForFields
			filename = getFieldsFileName(omega);
			file     = matopen(filename);
			U        = read(file,string("Ubatch_",k_batch));
			close(file);
		else
			U = pFor.Fields[:,batchIdxs];
		end
		# Original code:  V = conj((1.0.-1im*vec(gamma./omega)).*U).*V;
		V = An2cc*(conj(U).*(1.0.+1im*vec(gamma_nodal./omega)).*V);
		JTv .+= omega^2*vec(real(sum(V,dims=2)));
	end

	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		clear!(Ainv.MG);
	end

    return JTv
end


text:
"""
Holds neuron index and timestamp for each spike.
"""
struct Spike
    neuron::Int64
    timestamp::Float64
end

"""
Summary information for latent events. See `SeqEvent` for
the full-featured struct that caches things like sufficient
statistics. This summary holds the basic information that
is saved during MCMC sampling and output to the user.

assignment_id :

    Integer specifying the cluster id. Corresponds to
    the `indices` in `SeqEventList` struct.

timestamp :

    Time that the latent event occured. See
    `sampled_timestamp` in `SeqEvent`.

seq_type :

    Integer in {1, 2, ..., R} specifying the sequence
    type. See `sampled_type` in `SeqEvent` struct.

seq_warp:

    Float in [1/τ, ..., τ] specifying the inferred warp
    of the sequence.  See `sampled_warp` in `SeqEvent`;
    here we report the warp value rather than its index.

amplitude :

    Sampled amplitude of the latent event (i.e. the
    expected number of evoked spikes). See
    `sampled_amplitude` in `SeqEvent` struct.
"""
struct EventSummaryInfo
    assignment_id::Int64
    timestamp::Float64
    seq_type::Int64
    seq_warp::Float64
    amplitude::Float64
end

"""
Prior distributions for the SeqModel model.

seq_event_rate :

    Rate parameter, λ, for latent events. The number of
    sequence events, K, follows:

        K ~ Poisson(λ * max_time)

seq_type_proportions :
    Dirichlet prior over the relative proportions
    of R sequence event types.

        π ~ Dirichlet(γ * 1_R)

seq_event_amplitude :

    RateGamma prior on the amplitude (i.e., expected number of
    evoked spikes) of each sequence event.

        A_k ~ RateGamma(α, β)

neuron_response_proportions :

    Dirichlet prior over the assignment probabilities of
    sequence-evoked spikes across neurons.

        a ~ Dirichlet(φ * 1_N)

neuron_response_profile :

    Normal-Inverse-Chi-Squared distribution over the
    mean and variance (i.e. the delay and width) of
    each neuron's respons to a latent event.

        c_nr ~ Inv-χ²(ν, s)
        b_nr ~ Normal(0, c_nr / κ)

bkgd_amplitude :
    RateGamma prior on the rate parameter

        A_∅ ~ RateGamma(α_∅, β_∅)

    The number of background then follows a homogeneous
    Poisson process over the interval [0, T] with
    intensity A_bkgd.

bkgd_proportions :

    Dirichlet prior over the assignment probabilities of
    background spikes to neurons.

        a_∅ ~ Dirichlet(φ_∅ * 1_N)

warp_values :

    Length W vector of warp values, typically log-spaced
    from 1/τ to τ for some τ > 1.

warp_log_proportions :

    Length W log-probability vector over warp values.
"""
struct SeqPriors
    seq_event_rate::Float64

    seq_type_proportions::SymmetricDirichlet
    seq_event_amplitude::RateGamma

    neuron_response_proportions::SymmetricDirichlet
    neuron_response_profile::NormalInvChisq

    bkgd_amplitude::RateGamma
    bkgd_proportions::SymmetricDirichlet

    warp_values::Vector{Float64}
    warp_log_proportions::Vector{Float64}
end

"""
PP-Seq Global Parameters.

seq_type_log_proportions :
    Log-probability vector over sequence types. Element r
    is the probability that a sequence event has type r.

neuron_response_log_proportions :
    N × R matrix, each column is a log-probability vector over
    neurons. Element (n, r) is proportional to the response
    amplitude of neuron n to sequences of type r. Column r is
    denoted `log (a_r)` in the paper.

neuron_response_offsets :
    N × R matrix. Element (n, r) is the mean of the
    univariate Gaussian response profile for neuron n
    to sequences of type r.

neuron_response_widths :
    N × R matrix. Element (n, r) is the variance of the
    univariate Gaussian response profile for neuron n
    to sequences of type r.

bkgd_amplitude :
    The expected number of spikes (pooled across all neurons)
    in a unit time interval. Denoted `A_∅`
    in the paper.

bkgd_neuron_log_proportions :
    Log-probability vector over neurons. Element n is the
    log-probability that a background spike comes from neuron
    n, when chosen randomly. Denoted `a_∅` in the paper.
"""
mutable struct SeqGlobals
    seq_type_log_proportions::Vector{Float64} # R

    neuron_response_log_proportions::Matrix{Float64}  # N × R
    neuron_response_offsets::Matrix{Float64}  # N × R
    neuron_response_widths::Matrix{Float64}  # N × R

    bkgd_amplitude::Float64
    bkgd_log_proportions::Vector{Float64} # N

end

"""
Holds sufficient statistics needed to estimate the time,
amplitude, and type of a sequence occurence / latent event.

Parameters
----------
spike_count :
    Number of spikes assigned to this sequence.

summed_potentials :
   g Size RxW matrix, element (r,w) stores,

        Σᵢ mᵢ

        where mᵢ = (tᵢ - b_{nᵢ,r} * τ_w) / (c_{nᵢ,r} * τ_w)

    where the sum is over spikes (tᵢ, nᵢ) assigned to
    this sequence event. The parameters b_{r,nᵢ} and
    c_{r,nᵢ} correspond to the response offset and width
    of neuron nᵢ to sequences of type r, and they are each
    linearly scaled by the warp value τ_w.

summed_precisions :
    Size RxW matrix, element (r,w) stores,

        Σᵢ vᵢ

        where vᵢ = 1 / (c_{nᵢ,r} * τ_w)

    where the sum is over spikes (tᵢ, nᵢ) assigned to
    this sequence event. The parameter c_{r,nᵢ}
    corresponds to the response width of neuron nᵢ
    to sequences of type r, and it is linearly scaled
    by the warp value τ_w.

summed_logZ:
    Size RxW matrix, element (r,w) stores,

        Σᵢ log[ a_{nᵢ,r} / Z(mᵢ, vᵢ) ]

    where the sum is over spikes (tᵢ, nᵢ) assigned to
    this sequence event. The parameter a_{nᵢ,r} is
    the response amplitude of neuron nᵢ to sequences
    of type r. The term `Z(mᵢ, vᵢ)` denotes the log
    normalization constant for a univariate Gaussian
    with precision `vᵢ` and potential `mᵢ`.  Note that
    the precision and linear potentials depend on τ_w.

    TODO: This naming convention is a little strange
    since we're really storing log a -log Z.

seq_type_posterior:
    Size RxW matrix, element (r,w) is proportional to:

        Σᵢ(log a_{nᵢ,r} - Σᵢ log Z(mᵢ, vᵢ)) + log Z(Σᵢ mᵢ, Σᵢ vᵢ)

    After normalization, the represents the posterior distribution
    over the event sequence type and warp value, given the spikes
    currently assigned to this event.

sampled_type :
    The current type of the event (changes over
    Gibbs samples).

sampled_warp:
    The current warp index of the event (changes over
    Gibbs samples).

sampled_timestamp :
    The current time of the event (changes over
    Gibbs samples).

sampled_amplitude :
    The current amplitude of the event (changes over
    Gibbs samples).
"""
mutable struct SeqEvent
    spike_count::Int64
    summed_potentials::Matrix{Float64} # RxW
    summed_precisions::Matrix{Float64} # RxW
    summed_logZ::Matrix{Float64} # RxW
    seq_type_posterior::Matrix{Float64} # RxW

    sampled_type::Int64
    sampled_warp::Int64
    sampled_timestamp::Float64
    sampled_amplitude::Float64
end

"""
Dynamically re-sized array holding SeqEvent structs.

events :
    Vector of SeqEvent structs, some may be empty.

indices :
    Sorted vector of integers, specifying the indices of
    non-empty SeqEvent structs. Does not contain duplicate
    integers. Note that `length(indices) <= length(events)`,
    with equality only holding if there are no empty events.
"""
struct SeqEventList
    events::Vector{SeqEvent}
    indices::Vector{Int64}    # Indices of occupied sequences.
end

"""
Point Process Sequence Detection (PP-Seq) Model.

max_time :
    All spikes and sequences occur on the time interval
    between zero and max_time.

max_sequence_length :
    Maximum length of a sequence (used to speed up
    parent assignment step of collapsed Gibbs sampling
    -- we don't compute statistics for sequences larger
    than this threshold away from the spike.)

priors :
    Prior distributions.

globals :
    Global variables.

sequence_events :
    List of SeqEvent structs. See `./events.jl` for functionality.

_R_buffer :
    Vector of size `num_sequence_types`.

_RW_buffer :
    Matrix of size `num_sequence_types x num_warps`.

_K_buffer :
    Resized vector, holding probabilities over the number of latent events.
"""
mutable struct SeqModel
    max_time::Float64
    max_sequence_length::Float64

    priors::SeqPriors
    globals::SeqGlobals
    sequence_events::SeqEventList

    # TODO : move these to globals
    new_cluster_log_prob::Float64
    bkgd_log_prob::Float64

    _R_buffer::Vector{Float64}
    _RW_buffer::Matrix{Float64}
    _K_buffer::Vector{Float64}
end

"""
Distributed Point Process Sequence Model
"""
mutable struct DistributedSeqModel
    primary_model::SeqModel
    num_partitions::Int64
    submodels::Vector{SeqModel}
end


text:
using Documenter, DiffEqBase, SciMLBase, DiffEqProblemLibrary, OrdinaryDiffEq

ODEProblemLibrary = DiffEqProblemLibrary.ODEProblemLibrary
ODEProblemLibrary.importodeproblems()

SDEProblemLibrary = DiffEqProblemLibrary.SDEProblemLibrary
SDEProblemLibrary.importsdeproblems()

DDEProblemLibrary = DiffEqProblemLibrary.DDEProblemLibrary
DDEProblemLibrary.importddeproblems()

DAEProblemLibrary = DiffEqProblemLibrary.DAEProblemLibrary
DAEProblemLibrary.importdaeproblems()

include("pages.jl")

makedocs(modules=[DiffEqBase,SciMLBase,DiffEqProblemLibrary,ODEProblemLibrary,SDEProblemLibrary,DDEProblemLibrary,DAEProblemLibrary,OrdinaryDiffEq],
         doctest=false, clean=true,
         format = Documenter.HTML(analytics = "UA-90474609-3",
                                  assets = ["assets/favicon.ico"],
                                  canonical="https://diffeq.sciml.ai/stable/"),
         sitename="DifferentialEquations.jl",
         authors="Chris Rackauckas",
         pages = pages)

#Redirect old links
# cd(joinpath(@__DIR__, "build")) do
#     for (root, dirs, files) in walkdir(".")
#         for file in files
#             path = relpath(joinpath(root, file), ".")
#             m = match(r"(.+)/index\.html$", path)
#             m === nothing && continue
#             redirect = "$(m[1]).html"
#             @info "Adding redirect for $(m[1]) from $(redirect)"
#             isfile(redirect) && (@warn "$redirect exists, skip"; continue)
#             open(redirect, "w") do io
#                 write(io, """
#                 <meta http-equiv="refresh" content="0; url=$(basename(m[1]))/"/>
#                 """)
#             end
#         end
#     end
# end

deploydocs(
   repo = "github.com/SciML/DiffEqDocs.jl.git"
)

max_stars_repo_name: JuliaDiffEq/DiffEqDocs

max_issues_repo_name: JuliaDiffEq/DiffEqDocs

max_forks_repo_name: JuliaDiffEq/DiffEqDocs


text:
using LinearAlgebra, OrdinaryDiffEq, Test, PreallocationTools, ForwardDiff, LabelledArrays, RecursiveArrayTools

#Base Array tests
chunk_size = 5
u0_B = ones(5, 5)
dual_B = zeros(ForwardDiff.Dual{ForwardDiff.Tag{typeof(something), Float64}, Float64, chunk_size}, 2, 2)
cache_B = dualcache(u0_B, chunk_size)
tmp_du_BA = get_tmp(cache_B, u0_B)
tmp_dual_du_BA = get_tmp(cache_B, dual_B)
tmp_du_BN = get_tmp(cache_B, u0_B[1])
tmp_dual_du_BN = get_tmp(cache_B, dual_B[1])
@test size(tmp_du_BA) == size(u0_B)
@test typeof(tmp_du_BA) == typeof(u0_B)
@test eltype(tmp_du_BA) == eltype(u0_B)
@test size(tmp_dual_du_BA) == size(u0_B)
@test typeof(tmp_dual_du_BA) == typeof(dual_B)
@test eltype(tmp_dual_du_BA) == eltype(dual_B)
@test size(tmp_du_BN) == size(u0_B)
@test typeof(tmp_du_BN) == typeof(u0_B)
@test eltype(tmp_du_BN) == eltype(u0_B)
@test size(tmp_dual_du_BN) == size(u0_B)
@test typeof(tmp_dual_du_BN) == typeof(dual_B)
@test eltype(tmp_dual_du_BN) == eltype(dual_B)

#LArray tests
chunk_size = 4
u0_L = LArray((2,2); a=1.0, b=1.0, c=1.0, d=1.0)
zerodual = zero(ForwardDiff.Dual{ForwardDiff.Tag{typeof(something), Float64}, Float64, chunk_size})
dual_L = LArray((2,2); a=zerodual, b=zerodual, c=zerodual, d=zerodual)
cache_L = dualcache(u0_L, chunk_size)
tmp_du_LA = get_tmp(cache_L, u0_L)
tmp_dual_du_LA = get_tmp(cache_L, dual_L)
tmp_du_LN = get_tmp(cache_L, u0_L[1])
tmp_dual_du_LN = get_tmp(cache_L, dual_L[1])
@test size(tmp_du_LA) == size(u0_L)
@test typeof(tmp_du_LA) == typeof(u0_L)
@test eltype(tmp_du_LA) == eltype(u0_L)
@test size(tmp_dual_du_LA) == size(u0_L)
@test typeof(tmp_dual_du_LA) == typeof(dual_L)
@test eltype(tmp_dual_du_LA) == eltype(dual_L)
@test size(tmp_du_LN) == size(u0_L)
@test typeof(tmp_du_LN) == typeof(u0_L)
@test eltype(tmp_du_LN) == eltype(u0_L)
@test size(tmp_dual_du_LN) == size(u0_L)
@test typeof(tmp_dual_du_LN) == typeof(dual_L)
@test eltype(tmp_dual_du_LN) == eltype(dual_L)

#ArrayPartition tests
u0_AP = ArrayPartition(ones(2,2), ones(3,3))
dual_a = zeros(ForwardDiff.Dual{ForwardDiff.Tag{typeof(something), Float64}, Float64, chunk_size}, 2, 2)
dual_b = zeros(ForwardDiff.Dual{ForwardDiff.Tag{typeof(something), Float64}, Float64, chunk_size}, 3, 3)
dual_AP = ArrayPartition(dual_a, dual_b)
cache_AP = dualcache(u0_AP, chunk_size)
tmp_du_APA = get_tmp(cache_AP, u0_AP)
tmp_dual_du_APA = get_tmp(cache_AP, dual_AP)
tmp_du_APN = get_tmp(cache_AP, u0_AP[1])
tmp_dual_du_APN = get_tmp(cache_AP, dual_AP[1])
@test size(tmp_du_APA) == size(u0_AP)
@test typeof(tmp_du_APA) == typeof(u0_AP)
@test eltype(tmp_du_APA) == eltype(u0_AP)
@test size(tmp_dual_du_APA) == size(u0_AP)
@test typeof(tmp_dual_du_APA) == typeof(dual_AP)
@test eltype(tmp_dual_du_APA) == eltype(dual_AP)
@test size(tmp_du_APN) == size(u0_AP)
@test typeof(tmp_du_APN) == typeof(u0_AP)
@test eltype(tmp_du_APN) == eltype(u0_AP)
@test size(tmp_dual_du_APN) == size(u0_AP)
@test typeof(tmp_dual_du_APN) == typeof(dual_AP)
@test eltype(tmp_dual_du_APN) == eltype(dual_AP)


text:
const AnySSAValue = Union{Core.Compiler.SSAValue,JuliaInterpreter.SSAValue}
const AnySlotNumber = Union{Core.Compiler.SlotNumber,JuliaInterpreter.SlotNumber}

# to circumvent https://github.com/JuliaLang/julia/issues/37342, we inline these `isa`
# condition checks at surface AST level
# https://github.com/JuliaLang/julia/pull/38905 will get rid of the need of these hacks
macro isssa(stmt)
    :($(GlobalRef(Core, :isa))($(esc(stmt)), $(GlobalRef(Core.Compiler, :SSAValue))) ||
      $(GlobalRef(Core, :isa))($(esc(stmt)), $(GlobalRef(JuliaInterpreter, :SSAValue))))
end
macro issslotnum(stmt)
    :($(GlobalRef(Core, :isa))($(esc(stmt)), $(GlobalRef(Core.Compiler, :SlotNumber))) ||
      $(GlobalRef(Core, :isa))($(esc(stmt)), $(GlobalRef(JuliaInterpreter, :SlotNumber))))
end

"""
    iscallto(stmt, name)

Returns `true` is `stmt` is a call expression to `name`.
"""
function iscallto(@nospecialize(stmt), name)
    if isa(stmt, Expr)
        if stmt.head === :call
            a = stmt.args[1]
            a === name && return true
            is_global_ref(a, Core, :_apply) && stmt.args[2] === name && return true
            is_global_ref(a, Core, :_apply_iterate) && stmt.args[3] === name && return true
        end
    end
    return false
end

"""
    getcallee(stmt)

Returns the function (or Symbol) being called in a :call expression.
"""
function getcallee(@nospecialize(stmt))
    if isa(stmt, Expr)
        if stmt.head === :call
            a = stmt.args[1]
            is_global_ref(a, Core, :_apply) && return stmt.args[2]
            is_global_ref(a, Core, :_apply_iterate) && return stmt.args[3]
            return a
        end
    end
    error(stmt, " is not a call expression")
end

function callee_matches(f, mod, sym)
    is_global_ref(f, mod, sym) && return true
    if isdefined(mod, sym) && isa(f, QuoteNode)
        f.value === getfield(mod, sym) && return true  # a consequence of JuliaInterpreter.optimize!
    end
    return false
end

function rhs(stmt)
    isexpr(stmt, :(=)) && return (stmt::Expr).args[2]
    return stmt
end

ismethod(frame::Frame)  = ismethod(pc_expr(frame))
ismethod3(frame::Frame) = ismethod3(pc_expr(frame))

ismethod(stmt)  = isexpr(stmt, :method)
ismethod1(stmt) = isexpr(stmt, :method, 1)
ismethod3(stmt) = isexpr(stmt, :method, 3)

# anonymous function types are defined in a :thunk expr with a characteristic CodeInfo
function isanonymous_typedef(stmt)
    if isa(stmt, Expr)
        stmt.head === :thunk || return false
        stmt = stmt.args[1]
    end
    if isa(stmt, CodeInfo)
        src = stmt    # just for naming consistency
        length(src.code) >= 4 || return false
        stmt = src.code[end-1]
        (isexpr(stmt, :call) && is_global_ref(stmt.args[1], Core, :_typebody!)) || return false
        name = (stmt::Expr).args[2]::Symbol
        return startswith(String(name), "#")
    end
    return false
end

function istypedef(stmt)
    isa(stmt, Expr) || return false
    stmt = rhs(stmt)
    isa(stmt, Expr) || return false
    @static if all(s->isdefined(Core,s), structdecls)
        if stmt.head === :call
            f = stmt.args[1]
            if isa(f, GlobalRef)
                f.mod === Core && f.name ∈ structdecls && return true
            end
            if isa(f, QuoteNode)
                (f.value === Core._structtype || f.value === Core._abstracttype ||
                 f.value === Core._primitivetype) && return true
            end
        end
    end
    isanonymous_typedef(stmt) && return true
    return false
end

# Given a typedef at `src.code[idx]`, return the range of statement indices that encompass the typedef.
# The range does not include any constructor methods.
function typedef_range(src::CodeInfo, idx)
    stmt = src.code[idx]
    istypedef(stmt) || error(stmt, " is not a typedef")
    stmt = stmt::Expr
    isanonymous_typedef(stmt) && return idx:idx
    # Search backwards to the previous :global
    istart = idx
    while istart >= 1
        isexpr(src.code[istart], :global) && break
        istart -= 1
    end
    istart >= 1 || error("no initial :global found")
    iend, n = idx, length(src.code)
    have_typebody = have_equivtypedef = false
    while iend <= n
        stmt = src.code[iend]
        if isa(stmt, Expr)
            (stmt.head === :global || stmt.head === :return) && break
            if stmt.head === :call
                if (is_global_ref(stmt.args[1], Core, :_typebody!) ||
                    isdefined(Core, :_typebody!) && is_quotenode_egal(stmt.args[1], Core._typebody!))
                    have_typebody = true
                elseif (is_global_ref(stmt.args[1], Core, :_equiv_typedef) ||
                    isdefined(Core, :_equiv_typedef) && is_quotenode_egal(stmt.args[1], Core._equiv_typedef))
                    have_equivtypedef = true
                    # Advance to the type-assignment
                    while iend <= n
                        stmt = src.code[iend]
                        isexpr(stmt, :(=)) && break
                        iend += 1
                    end
                end
                if have_typebody && have_equivtypedef
                    iend += 1   # compensate for the `iend-1` in the return
                    break
                end
            end
        end
        is_return(stmt) && break
        iend += 1
    end
    iend <= n || (@show src; error("no final :global found"))
    return istart:iend-1
end

"""
    nextpc = next_or_nothing(frame, pc)
    nextpc = next_or_nothing!(frame)

Advance the program counter without executing the corresponding line.
If `frame` is finished, `nextpc` will be `nothing`.
"""
next_or_nothing(frame, pc) = pc < nstatements(frame.framecode) ? pc+1 : nothing
function next_or_nothing!(frame)
    pc = frame.pc
    if pc < nstatements(frame.framecode)
        frame.pc = pc = pc + 1
        return pc
    end
    return nothing
end

"""
    nextpc = skip_until(predicate, frame, pc)
    nextpc = skip_until!(predicate, frame)

Advance the program counter until `predicate(stmt)` return `true`.
"""
function skip_until(predicate, frame, pc)
    stmt = pc_expr(frame, pc)
    while !predicate(stmt)
        pc = next_or_nothing(frame, pc)
        pc === nothing && return nothing
        stmt = pc_expr(frame, pc)
    end
    return pc
end
function skip_until!(predicate, frame)
    pc = frame.pc
    stmt = pc_expr(frame, pc)
    while !predicate(stmt)
        pc = next_or_nothing!(frame)
        pc === nothing && return nothing
        stmt = pc_expr(frame, pc)
    end
    return pc
end

function sparam_ub(meth::Method)
    typs = []
    sig = meth.sig
    while sig isa UnionAll
        push!(typs, Symbol(sig.var.ub))
        sig = sig.body
    end
    return Core.svec(typs...)
end

showempty(list) = isempty(list) ? '∅' : list

# Smooth the transition between Core.Compiler and Base
rng(bb::Core.Compiler.BasicBlock) = (r = bb.stmts; return Core.Compiler.first(r):Core.Compiler.last(r))

function pushall!(dest, src)
    for item in src
        push!(dest, item)
    end
    return dest
end


text:
using Base: axes1, tail
using OffsetArrays
import Statistics
using Statistics: mean, var
import AxisArrays
using ImageMorphology: dilate!, erode!

# Compat.@dep_vectorize_2arg Gray atan2
# Compat.@dep_vectorize_2arg Gray hypot

"""
`M = meanfinite(img, region)` calculates the mean value along the dimensions listed in `region`, ignoring any non-finite values.
"""
meanfinite(A::AbstractArray{T}, region) where {T<:Real} = _meanfinite(A, T, region)
meanfinite(A::AbstractArray{CT}, region) where {CT<:Colorant} = _meanfinite(A, eltype(CT), region)
function _meanfinite(A::AbstractArray, ::Type{T}, region) where T<:AbstractFloat
    inds = Base.reduced_indices(A, region)
    K = similar(Array{Int}, inds)
    S = similar(Array{eltype(A)}, inds)
    fill!(K, 0)
    fill!(S, zero(eltype(A)))
    sumfinite!(S, K, A)
    S./K
end
_meanfinite(A::AbstractArray, ::Type, region) = mean(A, dims=region)  # non floating-point

using Base: check_reducedims, reducedim1, safe_tail
using Base.Broadcast: newindex

"""
    sumfinite!(S, K, A)

Compute the sum `S` and number of contributing pixels `K` for
reductions of the array `A` over dimensions. `S` and `K` must have
identical indices, and either match `A` or have singleton-dimensions for
the dimensions that are being summed over. Only pixels with finite
value are included in the tallies of `S` and `K`.

Note that the pixel mean is just S./K.
"""
function sumfinite!(S, K, A::AbstractArray{T,N}) where {T,N}
    check_reducedims(S, A)
    isempty(A) && return S, K
    axes(S) == axes(K) || throw(DimensionMismatch("S and K must have identical axes"))

    indsAt, indsSt = safe_tail(axes(A)), safe_tail(axes(S))
    keep, Idefault = _newindexer(indsSt)
    if reducedim1(S, A)
        # keep the accumulators as a local variable when reducing along the first dimension
        i1 = first(axes1(S))
        @inbounds for IA in CartesianIndices(indsAt)
            IS = newindex(IA, keep, Idefault)
            s, k = S[i1,IS], K[i1,IS]
            for i in axes(A, 1)
                tmp = A[i, IA]
                if isfinite(tmp)
                    s += tmp
                    k += 1
                end
            end
            S[i1,IS], K[i1,IS] = s, k
        end
    else
        @inbounds for IA in CartesianIndices(indsAt)
            IS = newindex(IA, keep, Idefault)
            for i in axes(A, 1)
                tmp = A[i, IA]
                if isfinite(tmp)
                    S[i, IS] += tmp
                    K[i, IS] += 1
                end
            end
        end
    end
    S, K
end
_newindexer(ax) = Base.Broadcast.shapeindexer(ax)

function Statistics.var(A::AbstractArray{C}; kwargs...) where C<:AbstractGray
    imgc = channelview(A)
    base_colorant_type(C)(var(imgc; kwargs...))
end

function Statistics.var(A::AbstractArray{C,N}; kwargs...) where {C<:Colorant,N}
    imgc = channelview(A)
    colons = ntuple(d->Colon(), Val(N))
    inds1 = axes(imgc, 1)
    val1 = Statistics.var(view(imgc, first(inds1), colons...); kwargs...)
    vals = similar(imgc, typeof(val1), inds1)
    vals[1] = val1
    for i in first(inds1)+1:last(inds1)
        vals[i] = Statistics.var(view(imgc, i, colons...); kwargs...)
    end
    base_colorant_type(C)(vals...)
end

Statistics.std(A::AbstractArray{C}; kwargs...) where {C<:Colorant} = mapc(sqrt, Statistics.var(A; kwargs...))

# Entropy for grayscale (intensity) images
function _log(kind::Symbol)
    if kind == :shannon
        return log2
    elseif kind == :nat
        return log
    elseif kind == :hartley
        return log10
    else
        throw(ArgumentError("Invalid entropy unit. (:shannon, :nat or :hartley)"))
    end
end

"""
    entropy(logᵦ, img)
    entropy(img; [kind=:shannon])

Compute the entropy of a grayscale image defined as `-sum(p.*logᵦ(p))`.
The base β of the logarithm (a.k.a. entropy unit) is one of the following:

- `:shannon ` (log base 2, default), or use logᵦ = log2
- `:nat` (log base e), or use logᵦ = log
- `:hartley` (log base 10), or use logᵦ = log10
"""
entropy(img::AbstractArray; kind=:shannon) = entropy(_log(kind), img)
function entropy(logᵦ::Log, img) where Log<:Function
    hist = StatsBase.fit(Histogram, vec(img), nbins=256, closed=:right)
    counts = hist.weights
    p = counts / length(img)
    logp = logᵦ.(p)

    # take care of empty bins
    logp[Bool[isinf(v) for v in logp]] .= 0

    -sum(p .* logp)
end

function entropy(img::AbstractArray{Bool}; kind=:shannon)
    logᵦ = _log(kind)

    p = sum(img) / length(img)

    (0 < p < 1) ? - p*logᵦ(p) - (1-p)*logᵦ(1-p) : zero(p)
end

entropy(img::AbstractArray{C}; kind=:shannon) where {C<:AbstractGray} = entropy(channelview(img), kind=kind)

"""
`m = minfinite(A)` calculates the minimum value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function minfinite(A::AbstractArray{T}) where T
    ret = sentinel_min(T)
    for a in A
        ret = minfinite_scalar(a, ret)
    end
    ret
end
function minfinite(f, A::AbstractArray)
    ret = sentinel_min(typeof(f(first(A))))
    for a in A
        ret = minfinite_scalar(f(a), ret)
    end
    ret
end

"""
`m = maxfinite(A)` calculates the maximum value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function maxfinite(A::AbstractArray{T}) where T
    ret = sentinel_max(T)
    for a in A
        ret = maxfinite_scalar(a, ret)
    end
    ret
end
function maxfinite(f, A::AbstractArray)
    ret = sentinel_max(typeof(f(first(A))))
    for a in A
        ret = maxfinite_scalar(f(a), ret)
    end
    ret
end

"""
`m = maxabsfinite(A)` calculates the maximum absolute value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function maxabsfinite(A::AbstractArray{T}) where T
    ret = sentinel_min(typeof(abs(A[1])))
    for a in A
        ret = maxfinite_scalar(abs(a), ret)
    end
    ret
end

minfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b < a ? b : a) : b
maxfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b > a ? b : a) : b
minfinite_scalar(a::T, b::T) where {T<:Union{Integer,FixedPoint}} = b < a ? b : a
maxfinite_scalar(a::T, b::T) where {T<:Union{Integer,FixedPoint}} = b > a ? b : a
minfinite_scalar(a, b) = minfinite_scalar(promote(a, b)...)
maxfinite_scalar(a, b) = maxfinite_scalar(promote(a, b)...)

function minfinite_scalar(c1::C, c2::C) where C<:AbstractRGB
    C(minfinite_scalar(c1.r, c2.r),
      minfinite_scalar(c1.g, c2.g),
      minfinite_scalar(c1.b, c2.b))
end
function maxfinite_scalar(c1::C, c2::C) where C<:AbstractRGB
    C(maxfinite_scalar(c1.r, c2.r),
      maxfinite_scalar(c1.g, c2.g),
      maxfinite_scalar(c1.b, c2.b))
end

sentinel_min(::Type{T}) where {T<:Union{Integer,FixedPoint}} = typemax(T)
sentinel_max(::Type{T}) where {T<:Union{Integer,FixedPoint}} = typemin(T)
sentinel_min(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
sentinel_max(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
sentinel_min(::Type{C}) where {C<:AbstractRGB} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C<:AbstractRGB,T} = (s = sentinel_min(T); C(s,s,s))
sentinel_max(::Type{C}) where {C<:AbstractRGB} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C<:AbstractRGB,T} = (s = sentinel_max(T); C(s,s,s))
sentinel_min(::Type{C}) where {C<:AbstractGray} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C<:AbstractGray,T} = C(sentinel_min(T))
sentinel_max(::Type{C}) where {C<:AbstractGray} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C<:AbstractGray,T} = C(sentinel_max(T))

# FIXME: replace with IntegralImage
# average filter
"""
`kern = imaverage(filtersize)` constructs a boxcar-filter of the specified size.
"""
function imaverage(filter_size=[3,3])
    if length(filter_size) != 2
        error("wrong filter size")
    end
    m, n = filter_size
    if mod(m, 2) != 1 || mod(n, 2) != 1
        error("filter dimensions must be odd")
    end
    f = ones(Float64, m, n)/(m*n)
end

# FIXME: do something about this
# more general version
function imlaplacian(alpha::Number)
    lc = alpha/(1 + alpha)
    lb = (1 - alpha)/(1 + alpha)
    lm = -4/(1 + alpha)
    return [lc lb lc; lb lm lb; lc lb lc]
end

accum(::Type{T}) where {T<:Integer} = Int
accum(::Type{Float32})    = Float32
accum(::Type{T}) where {T<:Real} = Float64
accum(::Type{C}) where {C<:Colorant} = base_colorant_type(C){accum(eltype(C))}

graytype(::Type{T}) where {T<:Number} = T
graytype(::Type{C}) where {C<:AbstractGray} = C
graytype(::Type{C}) where {C<:Colorant} = Gray{eltype(C)}

# Simple image difference testing
macro test_approx_eq_sigma_eps(A, B, sigma, eps)
    quote
        if size($(esc(A))) != size($(esc(B)))
            error("Sizes ", size($(esc(A))), " and ",
                  size($(esc(B))), " do not match")
        end
        kern = KernelFactors.IIRGaussian($(esc(sigma)))
        Af = imfilter($(esc(A)), kern, NA())
        Bf = imfilter($(esc(B)), kern, NA())
        diffscale = max(maxabsfinite($(esc(A))), maxabsfinite($(esc(B))))
        d = sad(Af, Bf)
        if d > length(Af)*diffscale*($(esc(eps)))
            error("Arrays A and B differ")
        end
    end
end

# image difference testing (@tbreloff's, based on the macro)
#   A/B: images/arrays to compare
#   sigma: tuple of ints... how many pixels to blur
#   eps: error allowance
# returns: percentage difference on match, error otherwise
function test_approx_eq_sigma_eps(A::AbstractArray, B::AbstractArray,
                         sigma::AbstractVector{T} = ones(ndims(A)),
                         eps::AbstractFloat = 1e-2,
                         expand_arrays::Bool = true) where T<:Real
    if size(A) != size(B)
        if expand_arrays
            newsize = map(max, size(A), size(B))
            if size(A) != newsize
                A = copyto!(zeros(eltype(A), newsize...), A)
            end
            if size(B) != newsize
                B = copyto!(zeros(eltype(B), newsize...), B)
            end
        else
            error("Arrays differ: size(A): $(size(A)) size(B): $(size(B))")
        end
    end
    if length(sigma) != ndims(A)
        error("Invalid sigma in test_approx_eq_sigma_eps. Should be ndims(A)-length vector of the number of pixels to blur.  Got: $sigma")
    end
    kern = KernelFactors.IIRGaussian(sigma)
    Af = imfilter(A, kern, NA())
    Bf = imfilter(B, kern, NA())
    diffscale = max(maxabsfinite(A), maxabsfinite(B))
    d = sad(Af, Bf)
    diffpct = d / (length(Af) * diffscale)
    if diffpct > eps
        error("Arrays differ.  Difference: $diffpct  eps: $eps")
    end
    diffpct
end

"""
BlobLoG stores information about the location of peaks as discovered by `blob_LoG`.
It has fields:

- location: the location of a peak in the filtered image (a CartesianIndex)
- σ: the value of σ which lead to the largest `-LoG`-filtered amplitude at this location
- amplitude: the value of the `-LoG(σ)`-filtered image at the peak

Note that the radius is equal to σ√2.

See also: [`blob_LoG`](@ref).
"""
struct BlobLoG{T,S,N}
    location::CartesianIndex{N}
    σ::S
    amplitude::T
end

"""
    blob_LoG(img, σscales, [edges], [σshape]) -> Vector{BlobLoG}

Find "blobs" in an N-D image using the negative Lapacian of Gaussians
with the specifed vector or tuple of σ values. The algorithm searches for places
where the filtered image (for a particular σ) is at a peak compared to all
spatially- and σ-adjacent voxels, where σ is `σscales[i] * σshape` for some i.
By default, `σshape` is an ntuple of 1s.

The optional `edges` argument controls whether peaks on the edges are
included. `edges` can be `true` or `false`, or a N+1-tuple in which
the first entry controls whether edge-σ values are eligible to serve
as peaks, and the remaining N entries control each of the N dimensions
of `img`.

# Citation:

Lindeberg T (1998), "Feature Detection with Automatic Scale Selection",
International Journal of Computer Vision, 30(2), 79–116.

See also: [`BlobLoG`](@ref).
"""
function blob_LoG(img::AbstractArray{T,N}, σscales::Union{AbstractVector,Tuple},
                  edges::Tuple{Vararg{Bool}}=(true, ntuple(d->false, Val(N))...), σshape=ntuple(d->1, Val(N))) where {T,N}
    sigmas = sort(σscales)
    img_LoG = Array{Float64}(undef, length(sigmas), size(img)...)
    colons = ntuple(d->Colon(), Val(N))
    @inbounds for isigma in eachindex(sigmas)
        img_LoG[isigma,colons...] = (-sigmas[isigma]) * imfilter(img, Kernel.LoG(ntuple(i->sigmas[isigma]*σshape[i],Val(N))))
    end
    maxima = findlocalmaxima(img_LoG, 1:ndims(img_LoG), edges)
    [BlobLoG(CartesianIndex(tail(x.I)), sigmas[x[1]], img_LoG[x]) for x in maxima]
end
blob_LoG(img::AbstractArray{T,N}, σscales, edges::Bool, σshape=ntuple(d->1, Val(N))) where {T,N} =
    blob_LoG(img, σscales, (edges, ntuple(d->edges,Val(N))...), σshape)

blob_LoG(img::AbstractArray{T,N}, σscales, σshape=ntuple(d->1, Val(N))) where {T,N} =
    blob_LoG(img, σscales, (true, ntuple(d->false,Val(N))...), σshape)

@inline function _clippedinds(Router,rstp)
    CartesianIndices(map((f,l)->f:l,
                         (first(Router)+rstp).I,(last(Router)-rstp).I))
end

findlocalextrema(img::AbstractArray{T,N}, region, edges::Bool, order) where {T,N} = findlocalextrema(img, region, ntuple(d->edges,Val(N)), order)

function findlocalextrema(img::AbstractArray{T,N}, region::Union{Tuple{Int,Vararg{Int}},Vector{Int},UnitRange{Int},Int}, edges::NTuple{N,Bool}, order::Base.Order.Ordering) where {T<:Union{Gray,Number},N}
    issubset(region,1:ndims(img)) || throw(ArgumentError("invalid region"))
    extrema = Array{CartesianIndex{N}}(undef, 0)
    edgeoffset = CartesianIndex(map(!, edges))
    R0 = CartesianIndices(axes(img))
    R = _clippedinds(R0,edgeoffset)
    rstp = one(first(R0))
    Rinterior = _clippedinds(R0,rstp)
    iregion = CartesianIndex(ntuple(d->d∈region, Val(N)))
    Rregion = CartesianIndices(map((f,l)->f:l,(-iregion).I, iregion.I))
    z = zero(iregion)
    for i in R
        isextrema = true
        img_i = img[i]
        if i ∈ Rinterior
            # If i is in the interior, we don't have to worry about i+j being out-of-bounds
            for j in Rregion
                j == z && continue
                if !Base.Order.lt(order, img[i+j], img_i)
                    isextrema = false
                    break
                end
            end
        else
            for j in Rregion
                (j == z || i+j ∉ R0) && continue
                if !Base.Order.lt(order, img[i+j], img_i)
                    isextrema = false
                    break
                end
            end
        end
        isextrema && push!(extrema, i)
    end
    extrema
end

"""
`findlocalmaxima(img, [region, edges]) -> Vector{CartesianIndex}`

Returns the coordinates of elements whose value is larger than all of
their immediate neighbors.  `region` is a list of dimensions to
consider.  `edges` is a boolean specifying whether to include the
first and last elements of each dimension, or a tuple-of-Bool
specifying edge behavior for each dimension separately.
"""
findlocalmaxima(img::AbstractArray, region=coords_spatial(img), edges=true) =
        findlocalextrema(img, region, edges, Base.Order.Forward)

"""
Like `findlocalmaxima`, but returns the coordinates of the smallest elements.
"""
findlocalminima(img::AbstractArray, region=coords_spatial(img), edges=true) =
        findlocalextrema(img, region, edges, Base.Order.Reverse)

# restrict for AxisArray and ImageMeta
restrict(img::AxisArray, ::Tuple{}) = img
restrict(img::ImageMeta, ::Tuple{}) = img

function restrict(img::ImageMeta, region::Dims)
    shareproperties(img, restrict(arraydata(img), region))
end

function restrict(img::AxisArray{T,N}, region::Dims) where {T,N}
    inregion = falses(ndims(img))
    inregion[[region...]] .= true
    inregiont = (inregion...,)::NTuple{N,Bool}
    rdata = restrict(img.data, region)
    AxisArray(rdata, map(modax, AxisArrays.axes(img), axes(rdata), inregiont))
end

# FIXME: this doesn't get inferred, but it should be (see issue #628)
function restrict(img::Union{AxisArray,ImageMetaAxis}, ::Type{Ax}) where Ax
    dim = axisdim(img, Ax)
    A = restrict(img.data, dim)
    AxisArray(A, replace_axis(modax(img[Ax], axes(A)[dim]), AxisArrays.axes(img)))
end

replace_axis(newax, axs) = _replace_axis(newax, axnametype(newax), axs...)
@inline _replace_axis(newax, ::Type{Ax}, ax::Ax, axs...) where {Ax} = (newax, _replace_axis(newax, Ax, axs...)...)
@inline _replace_axis(newax, ::Type{Ax}, ax, axs...) where {Ax} = (ax, _replace_axis(newax, Ax, axs...)...)
_replace_axis(newax, ::Type{Ax}) where {Ax} = ()

axnametype(ax::Axis{name}) where {name} = Axis{name}

ofaxes(::Base.OneTo, r) = r
ofaxes(axs::Any, r) = OffsetArray(r, axs)

function modax(ax, Aax)
    v = ax.val
    if iseven(length(v))
        return ax(ofaxes(Aax, range(first(v)-step(v)/2, stop=last(v)+step(v)/2, length=length(v)÷2 + 1)))
    else
        return ax(ofaxes(Aax, range(first(v)/1, stop=last(v)/1, length=length(v)÷2 + 1)))
    end
end

modax(ax, Aax, inregion::Bool) = inregion ? modax(ax, Aax) : ax

function imlineardiffusion(img::Array{T,2}, dt::AbstractFloat, iterations::Integer) where T
    u = img
    f = imlaplacian()
    for i = dt:dt:dt*iterations
        u = u + dt*imfilter(u, f, "replicate")
    end
    u
end

function imgaussiannoise(img::AbstractArray{T}, variance::Number, mean::Number) where T
    return img + sqrt(variance)*randn(size(img)) + mean
end

imgaussiannoise(img::AbstractArray{T}, variance::Number) where {T} = imgaussiannoise(img, variance, 0)
imgaussiannoise(img::AbstractArray{T}) where {T} = imgaussiannoise(img, 0.01, 0)

# image gradients

# forward and backward differences
# can be very helpful for discretized continuous models
forwarddiffy(u::AbstractMatrix) = [u[2:end,:]; u[end:end,:]] - u
forwarddiffx(u::AbstractMatrix) = [u[:,2:end] u[:,end:end]] - u
backdiffy(u::AbstractMatrix) = u - [u[1:1,:]; u[1:end-1,:]]
backdiffx(u::AbstractMatrix) = u - [u[:,1:1] u[:,1:end-1]]
function div(p::AbstractArray{T,3}) where T
    # Definition from the Chambolle citation below, between Eqs. 5 and 6
    # This is the adjoint of -forwarddiff
    inds = axes(p)[1:2]
    out = similar(p, inds)
    Router = CartesianIndices(inds)
    rstp = one(first(Router))
    Rinner = _clippedinds(Router,rstp)
    # Since most of the points are in the interior, compute them more quickly by avoiding branches
    for I in Rinner
        out[I] = p[I,1] - p[I[1]-1, I[2], 1] +
                 p[I,2] - p[I[1], I[2]-1, 2]
    end
    # Handle the edge points
    for I in EdgeIterator(Router, Rinner)
        out[I] = 0
        if I[1] == first(inds[1])
            out[I] += p[I, 1]
        elseif I[1] == last(inds[1])
            out[I] -= p[I[1]-1, I[2], 1]
        else
            out[I] += p[I,1] - p[I[1]-1, I[2], 1]
        end
        if I[2] == first(inds[2])
            out[I] += p[I, 2]
        elseif I[2] == last(inds[2])
            out[I] -= p[I[1], I[2]-1, 2]
        else
            out[I] += p[I,2] - p[I[1], I[2]-1, 2]
        end
    end
    out
end

"""
```
imgr = imROF(img, λ, iterations)
```

Perform Rudin-Osher-Fatemi (ROF) filtering, more commonly known as Total
Variation (TV) denoising or TV regularization. `λ` is the regularization
coefficient for the derivative, and `iterations` is the number of relaxation
iterations taken. 2d only.

See https://en.wikipedia.org/wiki/Total_variation_denoising and
Chambolle, A. (2004). "An algorithm for total variation minimization and applications".
    Journal of Mathematical Imaging and Vision. 20: 89–97
"""
function imROF(img::AbstractMatrix{T}, λ::Number, iterations::Integer) where T<:NumberLike
    # Total Variation regularized image denoising using the primal dual algorithm
    # Also called Rudin Osher Fatemi (ROF) model
    # λ: regularization parameter
    s1, s2 = size(img)
    p = zeros(T, s1, s2, 2)
    # This iterates Eq. (9) of the Chambolle citation
    local u
    τ = 1/4   # see 2nd remark after proof of Theorem 3.1.
    for i = 1:iterations
        div_p = div(p)
        u = img - λ*div_p # multiply term inside ∇ by -λ. Thm. 3.1 relates this to u via Eq. 7.
        grad_u = cat(forwarddiffy(u), forwarddiffx(u), dims=3)
        grad_u_mag = sqrt.(sum(abs2, grad_u, dims=3))
        p .= (p .- (τ/λ).*grad_u)./(1 .+ (τ/λ).*grad_u_mag)
    end
    return u
end

# ROF Model for color images
function imROF(img::AbstractMatrix{<:Color}, λ::Number, iterations::Integer)
    out = similar(img)
    imgc = channelview(img)
    outc = channelview(out)
    for chan = 1:size(imgc, 1)
        outc[chan, :, :] = imROF(imgc[chan, :, :], λ, iterations)
    end
    out
end

# morphological operations for ImageMeta
dilate(img::ImageMeta, region=coords_spatial(img)) = shareproperties(img, dilate!(copy(arraydata(img)), region))
erode(img::ImageMeta, region=coords_spatial(img)) = shareproperties(img, erode!(copy(arraydata(img)), region))

# phantom images

"""
```
phantom = shepp_logan(N,[M]; highContrast=true)
```

output the NxM Shepp-Logan phantom, which is a standard test image usually used
for comparing image reconstruction algorithms in the field of computed
tomography (CT) and magnetic resonance imaging (MRI). If the argument M is
omitted, the phantom is of size NxN. When setting the keyword argument
`highConstrast` to false, the CT version of the phantom is created. Otherwise,
the high contrast MRI version is calculated.
"""
function shepp_logan(M,N; highContrast=true)
  # Initially proposed in Shepp, Larry; B. F. Logan (1974).
  # "The Fourier Reconstruction of a Head Section". IEEE Transactions on Nuclear Science. NS-21.

  P = zeros(M,N)

  x = range(-1, stop=1, length=M)'
  y = range(1, stop=-1, length=N)

  centerX = [0, 0, 0.22, -0.22, 0, 0, 0, -0.08, 0, 0.06]
  centerY = [0, -0.0184, 0, 0, 0.35, 0.1, -0.1, -0.605, -0.605, -0.605]
  majorAxis = [0.69, 0.6624, 0.11, 0.16, 0.21, 0.046, 0.046, 0.046, 0.023, 0.023]
  minorAxis = [0.92, 0.874, 0.31, 0.41, 0.25, 0.046, 0.046, 0.023, 0.023, 0.046]
  theta = [0, 0, -18.0, 18.0, 0, 0, 0, 0, 0, 0]

  # original (CT) version of the phantom
  grayLevel = [2, -0.98, -0.02, -0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]

  if(highContrast)
    # high contrast (MRI) version of the phantom
    grayLevel = [1, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  end

  for l=1:length(theta)
    P += grayLevel[l] * (
           ( (cos(theta[l] / 360*2*pi) * (x .- centerX[l]) .+
              sin(theta[l] / 360*2*pi) * (y .- centerY[l])) / majorAxis[l] ).^2 .+
           ( (sin(theta[l] / 360*2*pi) * (x .- centerX[l]) .-
              cos(theta[l] / 360*2*pi) * (y .- centerY[l])) / minorAxis[l] ).^2 .< 1 )
  end

  return P
end

shepp_logan(N;highContrast=true) = shepp_logan(N,N;highContrast=highContrast)

"""
```
integral_img = integral_image(img)
```

Returns the integral image of an image. The integral image is calculated by assigning
to each pixel the sum of all pixels above it and to its left, i.e. the rectangle from
(1, 1) to the pixel. An integral image is a data structure which helps in efficient
calculation of sum of pixels in a rectangular subset of an image. See `boxdiff` for more
information.
"""
function integral_image(img::AbstractArray)
    integral_img = Array{accum(eltype(img))}(undef, size(img))
    sd = coords_spatial(img)
    cumsum!(integral_img, img, dims=sd[1])
    for i = 2:length(sd)
        cumsum!(integral_img, integral_img, dims=sd[i])
    end
    integral_img
end

"""
```
sum = boxdiff(integral_image, ytop:ybot, xtop:xbot)
sum = boxdiff(integral_image, CartesianIndex(tl_y, tl_x), CartesianIndex(br_y, br_x))
sum = boxdiff(integral_image, tl_y, tl_x, br_y, br_x)
```

An integral image is a data structure which helps in efficient calculation of sum of pixels in
a rectangular subset of an image. It stores at each pixel the sum of all pixels above it and to
its left. The sum of a window in an image can be directly calculated using four array
references of the integral image, irrespective of the size of the window, given the `yrange` and
`xrange` of the window. Given an integral image -

        A - - - - - - B -
        - * * * * * * * -
        - * * * * * * * -
        - * * * * * * * -
        - * * * * * * * -
        - * * * * * * * -
        C * * * * * * D -
        - - - - - - - - -

The sum of pixels in the area denoted by * is given by S = D + A - B - C.
"""
boxdiff(int_img::AbstractArray{T, 2}, y::UnitRange, x::UnitRange) where {T} = boxdiff(int_img, y.start, x.start, y.stop, x.stop)
boxdiff(int_img::AbstractArray{T, 2}, tl::CartesianIndex, br::CartesianIndex) where {T} = boxdiff(int_img, tl[1], tl[2], br[1], br[2])

function boxdiff(int_img::AbstractArray{T, 2}, tl_y::Integer, tl_x::Integer, br_y::Integer, br_x::Integer) where T
    sum = int_img[br_y, br_x]
    sum -= tl_x > 1 ? int_img[br_y, tl_x - 1] : zero(T)
    sum -= tl_y > 1 ? int_img[tl_y - 1, br_x] : zero(T)
    sum += tl_y > 1 && tl_x > 1 ? int_img[tl_y - 1, tl_x - 1] : zero(T)
    sum
end

"""
```
P = bilinear_interpolation(img, r, c)
```

Bilinear Interpolation is used to interpolate functions of two variables
on a rectilinear 2D grid.

The interpolation is done in one direction first and then the values obtained
are used to do the interpolation in the second direction.

"""
function bilinear_interpolation(img::AbstractArray{T, 2}, y::Number, x::Number) where T
    y_min = floor(Int, y)
    x_min = floor(Int, x)
    y_max = ceil(Int, y)
    x_max = ceil(Int, x)

    topleft = chkbounds(Bool, img, y_min, x_min) ? img[y_min, x_min] : zero(T)
    bottomleft = chkbounds(Bool, img, y_max, x_min) ? img[y_max, x_min] : zero(T)
    topright = chkbounds(Bool, img, y_min, x_max) ? img[y_min, x_max] : zero(T)
    bottomright = chkbounds(Bool, img, y_max, x_max) ? img[y_max, x_max] : zero(T)

    if x_max == x_min
        if y_max == y_min
            return T(topleft)
        end
        return T(((y_max - y) * topleft + (y - y_min) * bottomleft) / (y_max - y_min))
    elseif y_max == y_min
        return T(((x_max - x) * topleft + (x - x_min) * topright) / (x_max - x_min))
    end

    r1 = ((x_max - x) * topleft + (x - x_min) * topright) / (x_max - x_min)
    r2 = ((x_max - x) * bottomleft + (x - x_min) * bottomright) / (x_max - x_min)

    T(((y_max - y) * r1 + (y - y_min) * r2) / (y_max - y_min))

end

chkbounds(::Type{Bool}, img, x, y) = checkbounds(Bool, img, x, y)

"""
```
pyramid = gaussian_pyramid(img, n_scales, downsample, sigma)
```

Returns a  gaussian pyramid of scales `n_scales`, each downsampled
by a factor `downsample` > 1 and `sigma` for the gaussian kernel.

"""
function gaussian_pyramid(img::AbstractArray{T,N}, n_scales::Int, downsample::Real, sigma::Real) where {T,N}
    kerng = KernelFactors.IIRGaussian(sigma)
    kern = ntuple(d->kerng, Val(N))
    gaussian_pyramid(img, n_scales, downsample, kern)
end

function gaussian_pyramid(img::AbstractArray{T,N}, n_scales::Int, downsample::Real, kern::NTuple{N,Any}) where {T,N}
    downsample > 1 || @warn("downsample factor should be greater than one")
    # To guarantee inferability, we make sure that we do at least one
    # round of smoothing and resizing
    img_smoothed_main = imfilter(img, kern, NA())
    img_scaled = pyramid_scale(img_smoothed_main, downsample)
    prev = convert(typeof(img_scaled), img)
    pyramid = typeof(img_scaled)[prev]
    if n_scales ≥ 1
        # Take advantage of the work we've already done
        push!(pyramid, img_scaled)
        prev = img_scaled
    end
    for i in 2:n_scales
        img_smoothed = imfilter(prev, kern, NA())
        img_scaled = pyramid_scale(img_smoothed, downsample)
        push!(pyramid, img_scaled)
        prev = img_scaled
    end
    pyramid
end

function pyramid_scale(img, downsample)
    sz_next = map(s->ceil(Int, s/downsample), size(img))
    imresize(img, sz_next)
end

function pyramid_scale(img::OffsetArray, downsample)
    sz_next = map(s->ceil(Int, s/downsample), length.(axes(img)))
#    off = (.-ceil.(Int,(.-iterate.(axes(img).-(1,1))[1])./downsample))
    off = (.-ceil.(Int,(.-iterate.(map(x->UnitRange(x).-1,axes(img)))[1])./downsample))
    OffsetArray(imresize(img, sz_next), off)
end

"""
```
thres = otsu_threshold(img)
thres = otsu_threshold(img, bins)
```

Computes threshold for grayscale image using Otsu's method.

Parameters:
-    img         = Grayscale input image
-    bins        = Number of bins used to compute the histogram. Needed for floating-point images.

"""
function otsu_threshold(img::AbstractArray{T, N}, bins::Int = 256) where {T<:Union{Gray,Real}, N}

    min, max = extrema(img)
    edges, counts = imhist(img, range(gray(min), stop=gray(max), length=bins))
    histogram = counts./sum(counts)

    ω0 = 0
    μ0 = 0
    μt = 0
    μT = sum((1:(bins+1)).*histogram)
    max_σb=0.0
    thres=1

    for t in 1:bins
        ω0 += histogram[t]
        ω1 = 1 - ω0
        μt += t*histogram[t]

        σb = (μT*ω0-μt)^2/(ω0*ω1)

        if(σb > max_σb)
            max_σb = σb
            thres = t
        end
    end

    return T((edges[thres-1]+edges[thres])/2)
end

"""
```
thres = yen_threshold(img)
thres = yen_threshold(img, bins)
```

Computes threshold for grayscale image using Yen's maximum correlation criterion for bilevel thresholding

Parameters:
-    img         = Grayscale input image
-    bins        = Number of bins used to compute the histogram. Needed for floating-point images.

#Citation
Yen J.C., Chang F.J., and Chang S. (1995) “A New Criterion for Automatic Multilevel Thresholding” IEEE Trans. on Image Processing, 4(3): 370-378. DOI:10.1109/83.366472
"""
function yen_threshold(img::AbstractArray{T, N}, bins::Int = 256) where {T<:Union{Gray, Real}, N}

    min, max = extrema(img)
    if(min == max)
        return T(min)
    end

    edges, counts = imhist(img, range(gray(min), stop=gray(max), length=bins))

    prob_mass_function = counts./sum(counts)
    clamp!(prob_mass_function,eps(),Inf)
    prob_mass_function_sq = prob_mass_function.^2
    cum_pmf = cumsum(prob_mass_function)
    cum_pmf_sq_1 = cumsum(prob_mass_function_sq)
    cum_pmf_sq_2 = reverse!(cumsum(reverse!(prob_mass_function_sq)))

    #Equation (4) cited in the paper.
    criterion = log.(((cum_pmf[1:end-1].*(1.0 .- cum_pmf[1:end-1])).^2) ./ (cum_pmf_sq_1[1:end-1].*cum_pmf_sq_2[2:end]))

    thres = edges[findmax(criterion)[2]]
    return T(thres)

end


text:
"""
Our own little home-grown JSON conversion module (Muesli is a cereal...)

Types should overload `to_json` to support serialization
into json, and `from_json` to support deserialization out of json.

`to_json` and `from_json` do not actually convert into the string
format of json, they just convert into a subset of Julia types
(i.e. Dict, Vector, String, Int, Bool) that map cleanly into json.

WebIO handles the conversion between this and actual JSON.

We provide `generic_to_json` and `generic_from_json` that do
what one would expect on structs, mapping between a struct and a
dictionary of fields, converting the fields into/out of json as well.
"""
module Muesli

export to_json, generic_to_json, from_json, generic_from_json

"""
For a lot of types, (i.e. Int, String, Bool, etc.) the correct behavior
is for `to_json` to be the identity.
"""
to_json(x) = x

to_json(x::Symbol) = string(x)

to_json(x::Dict) = Dict(to_json(k) => to_json(v) for (k,v) in x)

to_json(x::Vector) = to_json.(x)

"""
For Dict{Symbol,X}, we want to convert to a dictionary with string keys
"""
to_json(x::Dict{Symbol,T}) where {T} = Dict(string(k)=>to_json(v) for (k,v) in x)

"""
For Dict{Int,X}, we want to convert to a vector of tuples
"""
to_json(x::Dict{Int,T}) where {T} = [(k,to_json(v)) for (k,v) in x]

"""
If x is a struct, this converts x into a dictionary whose keys are
the fieldnames of x, and whose values are the fields of x, converted.
"""
generic_to_json(x) = Dict(string(fn)=>to_json(getfield(x,fn)) for fn ∈ fieldnames(typeof(x)))

"""
`from_json` takes in a JSON value along with a type to convert into
and and attempts to convert to the type passed in.
"""
from_json(x,T) = convert(T,x) # fallback
from_json(x::String, ::Type{Symbol}) = Symbol(x)
from_json(x::String, ::Type{Union{Symbol, Nothing}}) = Symbol(x)
from_json(x::Symbol, ::Type{String}) = string(x)

function from_json(x::Vector{<:Any},::Type{Vector{T}}) where {T}
  map(y -> from_json(y,T), x)
end

function from_json(x::Vector{<:Any},::Type{Tuple{A,B}}) where {A,B}
  (from_json(x[1],A),from_json(x[2],B))
end

function from_json(x::Vector{<:Any},::Type{Pair{A,B}}) where {A,B}
  from_json(x[1],A) => from_json(x[2],B)
end

function from_json(x::Vector{<:Any},::Type{Dict{Int,T}}) where {T}
  Dict{Int,T}(from_json(x,Vector{Pair{Int,T}})...)
end

function from_json(d::Dict{String,<:Any}, ::Type{Dict{Symbol,T}}) where {T}
  Dict(Symbol(k) => from_json(v,T) for (k,v) in d)
end

function from_json(d::Vector{Tuple{Int,T}}, ::Type{Dict{Int,S}}) where {T,S}
  Dict(Int(i) => from_json(x,S) for (i,x) in d)
end

function from_json(s::String, ::Type{T}) where {T <: Real}
  parse(T,s)
end

function generic_from_json(d::Dict{String,<:Any},::Type{T}) where {T}
  augmented = Dict{Symbol,Any}()
  for (i,fn) in enumerate(fieldnames(T))
    sfn = string(fn)
    if sfn ∈ keys(d)
      augmented[fn] = from_json(d[sfn],fieldtypes(T)[i])
    else
      augmented[fn] = missing
    end
  end
  T([augmented[fn] for fn in fieldnames(T)]...)
end

end