text:
abstract PulseGenerator

type TwoExponentialPulseGenerator{T<:Real} <: PulseGenerator
	record_length::Int
    pre_rise_points::Int
	rise_points::Float64
	fall_points::Float64
	noise_fwhm::Float64
	event_rate_hz::Float64
    samples_per_second::Float64
    min_points_since_last_pulse::Float64
    amplitude_distribution::Distribution
    numrows::Int
    rownumber::Int
    quiescent_average::Float64
    quiescent_modulation_period_points::Float64
    quiescent_modulation_amplitude::Float64
    last_rowstamp::Int
	last_timestamp::Int
    trigger_points_std::Float64 # std of where it triggers relative to actualy arrival time in points (aka frames)
    d_amp_d_quiescent::Float64 # pretrigger mean/amplitude correlation
    function TwoExponentialPulseGenerator(record_length,pre_rise_points,rise_points,fall_points,noise_fwhm,event_rate_hz,samples_per_second,min_points_since_last_pulse,amplitude_distribution,numrows,rownumber,quiescent_average,quiescent_modulation_period_points,quiescent_modulation_amplitude,last_rowstamp,last_timestamp,trigger_points_std,d_amp_d_quiescent)
    	rise_points<fall_points || error("rise_points should be less than fall points")
    	new(record_length,pre_rise_points,rise_points,fall_points,noise_fwhm,event_rate_hz,samples_per_second,min_points_since_last_pulse,amplitude_distribution,numrows,rownumber,quiescent_average,quiescent_modulation_period_points,quiescent_modulation_amplitude,last_rowstamp,last_timestamp,trigger_points_std,d_amp_d_quiescent)
    end
end
TupacLikeTwoExponentialPulseGenerator{T}(::Type{T},distribution) = TwoExponentialPulseGenerator{T}(520, 100, 50, 200, 5, 13.5,
	104166.6, 520, distribution, 30, 0, 1000, 100_000*3600, 50, 0,0, 0.5, -0.0001)

"Generate samples of a noiseless record made up of one or more two exponential pulses. The number of samples is `points`,
the rising time constant as a number of points is roughly `rise_points`, the falling time constant in number of points is roughly
`fall points`. The next two parameters are tuples of the same length, `arrival_point` and `amplitude` where the first value in each
applies to the first pulse, further values of each (if they exist) add additional pulses to the record."
function two_exponential_pulses(
	points::Int, rise_points, fall_points,quiescent_value, arrival_point::Tuple, amplitude::Tuple)
	pulse = Array(Float64, points)
	fill!(pulse, quiescent_value)
	@assert length(arrival_point) == length(amplitude)
	k=kfactor(1/rise_points,1/fall_points)
	for (ap, amp) in zip(arrival_point, amplitude)
		@inbounds for i = 1:points
			t = i-ap
			pulse[i] += t>0.0?amp*k*(exp(-t/rise_points)-exp(-t/fall_points)):0.0
		end
	end
	pulse
end

# factor for calculating mathematical parameters from physical parameters and maintainging
# constant amplitude
# Otto-von-Guericke Univ. Magdeburg, Magdeburg, Germany
# IEEE Transactions on Electromagnetic Compatibility (Impact Factor: 1.35). 12/2010; DOI: 10.1109/TEMC.2010.2052621
# Estimation of the Mathematical Parameters of Double-Exponential Pulses Using the Nelder–Mead Algorithm - ResearchGate. Available from: http://www.researchgate.net/publication/224153186_Estimation_of_the_Mathematical_Parameters_of_Double-Exponential_Pulses_Using_the_NelderMead_Algorithm [accessed Mar 20, 2015].kfactor(a,b) = 1/(exp(-a*(log(a/b)/(a-b)))-exp(-b*(log(a/b)/(a-b))))
kfactor(a,b) = 1/(exp(-a*(log(a/b)/(a-b)))-exp(-b*(log(a/b)/(a-b)))) # memoize?
# const X = [0.842411, 1.660852, 8.887335, 0.624129]
# const Y = [2.167370, 0.360349, 0.013289, 1.358701, 0.137057, 122.802861, 1.314850]
# alpha(t_rise, t_fwhm) = (X[1]/((t_fwhm/t_rise)^X[2]-X[3])^X[4])/t_rise
# beta(t_rise, t_fwhm) = (Y[1]-Y[2]*exp(-t_fwhm*Y[3]/t_rise)-Y[4]*exp(-t_fwhm*Y[5]/t_rise)-Y[6]*exp(-t_fwhm*Y[7]/t_rise))/t_rise

white_noise(points::Int, noise_fwhm) = rand(points)*noise_fwhm/(2*sqrt(2*log(2))) # fwhm to std
function white_noise!(pulse, noise_fwhm)
	@inbounds for i = 1:length(pulse) pulse[i]+=randn()*noise_fwhm end
	pulse
end
to_type_and_white_noise{T<:Integer}(::Type{T}, noise_fwhm, pulse) = convert(Vector{T}, round(white_noise!(pulse, noise_fwhm)))

function getcleanpulse{T}(pg::TwoExponentialPulseGenerator{T}, amp)
	to_type_and_white_noise(T, pg.noise_fwhm,
		two_exponential_pulses(pg.record_length, pg.rise_points, pg.fall_points, pg.quiescent_average, (pg.pre_rise_points+1,), (amp,)))
end
function gettriggeredpulse!{T}(pg::TwoExponentialPulseGenerator{T})
	points_from_last, points_to_next = rand(Exponential(pg.samples_per_second/pg.event_rate_hz),2)
	arrival_points = (pg.pre_rise_points-pg.min_points_since_last_pulse-points_from_last+1
		, pg.pre_rise_points+rand(Normal(0,pg.trigger_points_std))+1, pg.pre_rise_points+points_to_next+1)

	last_rowstamp = pg.last_rowstamp
	last_timestamp = pg.last_timestamp
	rowstamp = round(last_rowstamp+points_from_last*pg.numrows)
	timestamp = round(last_timestamp + points_from_last*pg.samples_per_second*1e6)
	pg.last_rowstamp = rowstamp
	pg.last_timestamp = timestamp

	quiescent_value_modulation = pg.quiescent_modulation_amplitude*sinpi(2*rowstamp/(pg.quiescent_modulation_period_points*pg.numrows))
	quiescent_value = pg.quiescent_average+quiescent_value_modulation

	amplitudes = rand(pg.amplitude_distribution,3)

	# introduce pretrigger mean/amplitude correlation
	amplitudes ./= 1+pg.d_amp_d_quiescent*quiescent_value_modulation
	amplitudes_tuple = tuple(amplitudes...)

	pulse = two_exponential_pulses(pg.record_length, pg.rise_points, pg.fall_points, quiescent_value, arrival_points, amplitudes_tuple)
	pulse = to_type_and_white_noise(T, pg.noise_fwhm, pulse)
	pulse, rowstamp, timestamp
end
function gettriggeredpulse!{T}(pg::TwoExponentialPulseGenerator{T},n)
	pulses = Array(Vector{T},n)
	rowstamps = Array(Int, n)
	timestamps = Array(Int, n)
	for i = 1:n
		pulse, rowstamp, timestamp = gettriggeredpulse!(pg)
		pulses[i] = pulse
		rowstamps[i] = rowstamp
		timestamps[i] = timestamp
	end
	pulses, rowstamps, timestamps
end
times_s(pg) = range(0,1/pg.samples_per_second, pg.record_length)

# Base.start(pg::TwoExponentialPulseGenerator) = 0.0
# function Base.next(pg::TwoExponentialPulseGenerator, last_arrival::Float64)
# 	timestamp = last_arrival+rand(Exponential(pg.samples_per_second/pg.event_rate_hz),2)
# 	pulse = gettriggeredpulse(pg)
# 	return ((timestamp,pulse), timestamp)
# end
# Base.done(pg::TwoExponentialPulseGenerator, last_arrival::Float64) = false

# function getrowstamp(pg::TwoExponentialPulseGenerator,n::Int, last_rowstamp::Int=0)
# 	rowstamp_diffs = pg.numrows*rand(Exponential(pg.samples_per_second/pg.event_rate_hz),n)
# 	last_rowstamp = pg.numrows*div(last_rowstamp,pg.numrows)+pg.rownumber # make sure rowstamp reflect the correct row
# 	rowstamp_diffs[1]+=last_rowstamp
# 	cumsum(convert(Vector{Int}, round(rowstamp_diffs)))
# end


text:
function solve()
    a, b = [parse(Int, x) for x in split(readline())]
    c, d = [parse(Int, x) for x in split(readline())]

    a*d - b*c
end

println(solve())

ext: jl

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Dates

abstract type Event end

struct Position
    Name::AbstractString
    RequiredQualifications::AbstractArray{AbstractString}
    Position(name::AbstractString,rq::AbstractArray{String}) = new(name,rq)
    Position(name::AbstractString) = new(name,AbstractString[])
end

mutable struct FlightEvent <: Event
    Start::DateTime
    End::DateTime
    Positions::AbstractArray{Position}
    Participants::AbstractArray{Aircrew}
    FlightEvent(
        s::DateTime,
        e::DateTime,
        pos::AbstractArray{Position},
        part::AbstractArray{Aircrew}
    ) = new(s,e,pos,part)
    FlightEvent(
        s::DateTime,
        e::DateTime,
        pos::AbstractArray{Position}
    ) = FlightEvent(s,e,pos,Aircrew[])
end

mutable struct GroundEvent <: Event
    Start::DateTime
    End::DateTime
    Participants::AbstractArray{Aircrew}
    GroundEvent(
        s::Date,
        e::Date,
        part::AbstractArray{Aircrew}
    ) = new(s,e,part)
    GroundEvent(
        s::Date,
        e::Date,
    ) = GroundEvent(s,e,Aircrew[])
end

# isValidPerson(::Position, ::Person)
#
#     Compares if the person has the required qualifications for the position.
#     Returns True if they do
#     Returns False if they do not
function isValidPersonForPosition(event::Event,position::Position,person::Person)
    # Require Aircrew for FlightEvent
    if typeof(event)==FlightEvent && typeof(person)!=Aircrew
        return false
    end
    # Shortcut if no qualifications required, anyone is valid
    if position.RequiredQualifications == []
        return true
    end
    # Cycle through required qualifications
    for qual in position.RequiredQualifications
        # Check if the qualification is in persons qualifications list
        if !(qual ∈ person.Qualifications)
            # If required qualification is not in the persons qualification list, they are not valid
            return false
        else
            # If the required qualification IS in their list, verify its currency
            validQual = false
            # Get the persons qualification object that matches the required qualification
            for pqual in person.Qualifications[qual .== person.Qualifications]
                # If the persons qualification cannot expire OR  its date is valid for the event time, they are valid. No need to keep searching.
                if !pqual.CanExpire || (pqual.Granted < event.Start  && event.End < pqual.Expires)
                    validQual = true
                    break
                end
            end
            if !validQual
                # After looking at the persons qualification, if it is not valid for the event, they are not valid
                return false
            end
        end
    end
    return true
end


text:
"""
    report(res::IS.Results, out_path::String)

This function uses weave to either generate a LaTeX or HTML
file based on the report_design.jmd (julia markdown) file
that it reads. Out_path in the weave function dictates
where the created file gets exported.

# Arguments
- `results::IS.Results`: The results to be plotted
- `out_path::String`: folder path to the location the report should be generated

# Example
```julia
results = solve_op_problem!(OpModel)
out_path = "/Users/downloads"
report(results, out_path)
```

# Accepted Key Words
- `doctype::String = "md2html"`: create an HTML, default is PDF via latex
- `backend::Plots.backend() = plotlyjs()`: sets the plots backend, default is gr()
- `jmd::String = "file_path"`: directs the function to the julia markdown report design, the default
is ".../pwd()/report_design/report_design.jmd". To change the report, use report_design.jmd as a template
and add additional markdowns to the report_design folder.
"""
function report(res::IS.Results, out_path::String; kwargs...)

    doctype = get(kwargs, :doctype, "md2pdf")
    backend = get(kwargs, :backend, Plots.gr())
    default_string = joinpath(pwd(), "examples/report_design/report_design.jmd")
    jmd = get(kwargs, :jmd, default_string)
    args = Dict("res" => res, "variables" => res.variables, "backend" => backend)
    Weave.weave(
        jmd,
        out_path = out_path,
        latex_cmd = "xelatex",
        doctype = doctype,
        args = args,
    )

end

"""
    report(res::IS.Results, generators::Dict, out_path::String)

This function uses weave to either generate a LaTeX or HTML
file based on the report_design.jmd (julia markdown) file
that it reads. Out_path in the weave function dictates
where the created file gets exported.

# Arguments
- `results::IS.Results`: The results to be plotted
- `generators::Dict`: the dictionary of generators and their fuel type made by make_fuel_dictionary
- `out_path::String`: folder path to the location the report should be generated

# Example
```julia
results = solve_op_problem!(OpModel)
generators = make_fuel_dictionary(PSY.system, results)
out_path = "/Users/downloads"
report(results, generators, out_path)
```

# Accepted Key Words
- `doctype::String = "md2html"`: create an HTML, default is PDF via latex
- `backend::Plots.backend() = plotlyjs()`: sets the plots backend, default is gr()
- `jmd::String = "file_path"`: directs the function to the julia markdown report design, the default
is ".../pwd()/report_design/report_design_fuel.jmd". To change the report, use report_design_fuel.jmd as a template
and add additional markdowns to the report_design folder.
"""

function report(res::IS.Results, generators::Dict, out_path::String; kwargs...)

    doctype = get(kwargs, :doctype, "md2pdf")
    backend = get(kwargs, :backend, Plots.gr())
    default_string = joinpath(pwd(), "examples/report_design/report_design_fuel.jmd")
    jmd = get(kwargs, :jmd, default_string)
    args = Dict(
        "res" => res,
        "variables" => res.variables,
        "gen" => generators,
        "backend" => backend,
    )
    Weave.weave(
        jmd,
        out_path = out_path,
        latex_cmd = "xelatex",
        doctype = doctype,
        args = args,
    )

end

"""
    report(res::IS.Results, system::PSY.System, out_path::String)

This function uses weave to either generate a LaTeX or HTML
file based on the report_design.jmd (julia markdown) file
that it reads. Out_path in the weave function dictates
where the created file gets exported.

# Arguments
- `results::IS.Results`: The results to be plotted
- `system::PSY.System`: the system used to create the results for sorting by fuel type
- `out_path::String`: folder path to the location the report should be generated

# Example
```julia
results = solve_op_problem!(OpModel)
out_path = "/Users/downloads"
report(results, PSY.system, out_path)
```

# Accepted Key Words
- `doctype::String = "md2html"`: create an HTML, default is PDF via latex
- `backend::Plots.backend() = plotlyjs()`: sets the plots backend, default is gr()
- `jmd::String = "file_path"`: directs the function to the julia markdown report design, the default
is ".../pwd()/report_design/report_design_fuel.jmd". To change the report, use report_design_fuel.jmd as a template
and add additional markdowns to the report_design folder.
"""

function report(res::IS.Results, system::PSY.System, out_path::String; kwargs...)

    doctype = get(kwargs, :doctype, "md2pdf")
    backend = get(kwargs, :backend, gr())
    default_string = joinpath(pwd(), "examples/report_design/report_design_fuel.jmd")
    jmd = get(kwargs, :jmd, default_string)
    args = Dict(
        "res" => res,
        "variables" => res.variables,
        "gen" => system,
        "backend" => backend,
    )
    Weave.weave(
        jmd,
        out_path = out_path,
        latex_cmd = "xelatex",
        doctype = doctype,
        args = args,
    )

end


text:
###### Weight vector #####
abstract type AbstractWeights{S<:Real, T<:Real, V<:AbstractVector{T}} <: AbstractVector{T} end

"""
    @weights name

Generates a new generic weight type with specified `name`, which subtypes `AbstractWeights`
and stores the `values` (`V<:RealVector`) and `sum` (`S<:Real`).
"""
macro weights(name)
    return quote
        mutable struct $name{S<:Real, T<:Real, V<:AbstractVector{T}} <: AbstractWeights{S, T, V}
            values::V
            sum::S
        end
        $(esc(name))(vs) = $(esc(name))(vs, sum(vs))
    end
end

eltype(wv::AbstractWeights) = eltype(wv.values)
length(wv::AbstractWeights) = length(wv.values)
values(wv::AbstractWeights) = wv.values
sum(wv::AbstractWeights) = wv.sum
isempty(wv::AbstractWeights) = isempty(wv.values)

Base.getindex(wv::AbstractWeights, i) = getindex(wv.values, i)
Base.size(wv::AbstractWeights) = size(wv.values)

@propagate_inbounds function Base.setindex!(wv::AbstractWeights, v::Real, i::Int)
    s = v - wv[i]
    wv.values[i] = v
    wv.sum += s
    v
end

"""
    varcorrection(n::Integer, corrected=false)

Compute a bias correction factor for calculating `var`, `std` and `cov` with
`n` observations. Returns ``\\frac{1}{n - 1}`` when `corrected=true`
(i.e. [Bessel's correction](https://en.wikipedia.org/wiki/Bessel's_correction)),
otherwise returns ``\\frac{1}{n}`` (i.e. no correction).
"""
@inline varcorrection(n::Integer, corrected::Bool=false) = 1 / (n - Int(corrected))

@weights Weights

@doc """
    Weights(vs, wsum=sum(vs))

Construct a `Weights` vector with weight values `vs`.
A precomputed sum may be provided as `wsum`.

The `Weights` type describes a generic weights vector which does not support
all operations possible for [`FrequencyWeights`](@ref), [`AnalyticWeights`](@ref)
and [`ProbabilityWeights`](@ref).
""" Weights

"""
    weights(vs)

Construct a `Weights` vector from array `vs`.
See the documentation for [`Weights`](@ref) for more details.
"""
weights(vs::RealVector) = Weights(vs)
weights(vs::RealArray) = Weights(vec(vs))

"""
    varcorrection(w::Weights, corrected=false)

Returns ``\\frac{1}{\\sum w}`` when `corrected=false` and throws an `ArgumentError`
if `corrected=true`.
"""
@inline function varcorrection(w::Weights, corrected::Bool=false)
    corrected && throw(ArgumentError("Weights type does not support bias correction: " *
                                     "use FrequencyWeights, AnalyticWeights or ProbabilityWeights if applicable."))
    1 / w.sum
end

@weights AnalyticWeights

@doc """
    AnalyticWeights(vs, wsum=sum(vs))

Construct an `AnalyticWeights` vector with weight values `vs`.
A precomputed sum may be provided as `wsum`.

Analytic weights describe a non-random relative importance (usually between 0 and 1)
for each observation. These weights may also be referred to as reliability weights,
precision weights or inverse variance weights. These are typically used when the observations
being weighted are aggregate values (e.g., averages) with differing variances.
""" AnalyticWeights

"""
    aweights(vs)

Construct an `AnalyticWeights` vector from array `vs`.
See the documentation for [`AnalyticWeights`](@ref) for more details.
"""
aweights(vs::RealVector) = AnalyticWeights(vs)
aweights(vs::RealArray) = AnalyticWeights(vec(vs))

"""
    varcorrection(w::AnalyticWeights, corrected=false)

* `corrected=true`: ``\\frac{1}{\\sum w - \\sum {w^2} / \\sum w}``
* `corrected=false`: ``\\frac{1}{\\sum w}``
"""
@inline function varcorrection(w::AnalyticWeights, corrected::Bool=false)
    s = w.sum

    if corrected
        sum_sn = sum(x -> (x / s) ^ 2, w)
        1 / (s * (1 - sum_sn))
    else
        1 / s
    end
end

@weights FrequencyWeights

@doc """
    FrequencyWeights(vs, wsum=sum(vs))

Construct a `FrequencyWeights` vector with weight values `vs`.
A precomputed sum may be provided as `wsum`.

Frequency weights describe the number of times (or frequency) each observation
was observed. These weights may also be referred to as case weights or repeat weights.
""" FrequencyWeights

"""
    fweights(vs)

Construct a `FrequencyWeights` vector from a given array.
See the documentation for [`FrequencyWeights`](@ref) for more details.
"""
fweights(vs::RealVector) = FrequencyWeights(vs)
fweights(vs::RealArray) = FrequencyWeights(vec(vs))

"""
    varcorrection(w::FrequencyWeights, corrected=false)

* `corrected=true`: ``\\frac{1}{\\sum{w} - 1}``
* `corrected=false`: ``\\frac{1}{\\sum w}``
"""
@inline function varcorrection(w::FrequencyWeights, corrected::Bool=false)
    s = w.sum

    if corrected
        1 / (s - 1)
    else
        1 / s
    end
end

@weights ProbabilityWeights

@doc """
    ProbabilityWeights(vs, wsum=sum(vs))

Construct a `ProbabilityWeights` vector with weight values `vs`.
A precomputed sum may be provided as `wsum`.

Probability weights represent the inverse of the sampling probability for each observation,
providing a correction mechanism for under- or over-sampling certain population groups.
These weights may also be referred to as sampling weights.
""" ProbabilityWeights

"""
    pweights(vs)

Construct a `ProbabilityWeights` vector from a given array.
See the documentation for [`ProbabilityWeights`](@ref) for more details.
"""
pweights(vs::RealVector) = ProbabilityWeights(vs)
pweights(vs::RealArray) = ProbabilityWeights(vec(vs))

"""
    varcorrection(w::ProbabilityWeights, corrected=false)

* `corrected=true`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)`
* `corrected=false`: ``\\frac{1}{\\sum w}``
"""
@inline function varcorrection(w::ProbabilityWeights, corrected::Bool=false)
    s = w.sum

    if corrected
        n = count(!iszero, w)
        n / (s * (n - 1))
    else
        1 / s
    end
end

##### Equality tests #####

for w in (AnalyticWeights, FrequencyWeights, ProbabilityWeights, Weights)
    @eval begin
        Base.isequal(x::$w, y::$w) = isequal(x.sum, y.sum) && isequal(x.values, y.values)
        Base.:(==)(x::$w, y::$w)   = (x.sum == y.sum) && (x.values == y.values)
    end
end

Base.isequal(x::AbstractWeights, y::AbstractWeights) = false
Base.:(==)(x::AbstractWeights, y::AbstractWeights)   = false

##### Weighted sum #####

## weighted sum over vectors

"""
    wsum(v, w::AbstractVector, [dim])

Compute the weighted sum of an array `v` with weights `w`, optionally over the dimension `dim`.
"""
wsum(v::AbstractVector, w::AbstractVector) = dot(v, w)
wsum(v::AbstractArray, w::AbstractVector) = dot(vec(v), w)

# Note: the methods for BitArray and SparseMatrixCSC are to avoid ambiguities
Base.sum(v::BitArray, w::AbstractWeights) = wsum(v, values(w))
Base.sum(v::SparseArrays.SparseMatrixCSC, w::AbstractWeights) = wsum(v, values(w))
Base.sum(v::AbstractArray, w::AbstractWeights) = dot(v, values(w))

## wsum along dimension
#
#  Brief explanation of the algorithm:
#  ------------------------------------
#
#  1. _wsum! provides the core implementation, which assumes that
#     the dimensions of all input arguments are consistent, and no
#     dimension checking is performed therein.
#
#     wsum and wsum! perform argument checking and call _wsum!
#     internally.
#
#  2. _wsum! adopt a Cartesian based implementation for general
#     sub types of AbstractArray. Particularly, a faster routine
#     that keeps a local accumulator will be used when dim = 1.
#
#     The internal function that implements this is _wsum_general!
#
#  3. _wsum! is specialized for following cases:
#     (a) A is a vector: we invoke the vector version wsum above.
#         The internal function that implements this is _wsum1!
#
#     (b) A is a dense matrix with eltype <: BlasReal: we call gemv!
#         The internal function that implements this is _wsum2_blas!
#
#     (c) A is a contiguous array with eltype <: BlasReal:
#         dim == 1: treat A like a matrix of size (d1, d2 x ... x dN)
#         dim == N: treat A like a matrix of size (d1 x ... x d(N-1), dN)
#         otherwise: decompose A into multiple pages, and apply _wsum2!
#         for each
#
#     (d) A is a general dense array with eltype <: BlasReal:
#         dim <= 2: delegate to (a) and (b)
#         otherwise, decompose A into multiple pages
#

function _wsum1!(R::AbstractArray, A::AbstractVector, w::AbstractVector, init::Bool)
    r = wsum(A, w)
    if init
        R[1] = r
    else
        R[1] += r
    end
    return R
end

function _wsum2_blas!(R::StridedVector{T}, A::StridedMatrix{T}, w::StridedVector{T}, dim::Int, init::Bool) where T<:BlasReal
    beta = ifelse(init, zero(T), one(T))
    trans = dim == 1 ? 'T' : 'N'
    BLAS.gemv!(trans, one(T), A, w, beta, R)
    return R
end

function _wsumN!(R::StridedArray{T}, A::StridedArray{T,N}, w::StridedVector{T}, dim::Int, init::Bool) where {T<:BlasReal,N}
    if dim == 1
        m = size(A, 1)
        n = div(length(A), m)
        _wsum2_blas!(view(R,:), reshape(A, (m, n)), w, 1, init)
    elseif dim == N
        n = size(A, N)
        m = div(length(A), n)
        _wsum2_blas!(view(R,:), reshape(A, (m, n)), w, 2, init)
    else # 1 < dim < N
        m = 1
        for i = 1:dim-1; m *= size(A, i); end
        n = size(A, dim)
        k = 1
        for i = dim+1:N; k *= size(A, i); end
        Av = reshape(A, (m, n, k))
        Rv = reshape(R, (m, k))
        for i = 1:k
            _wsum2_blas!(view(Rv,:,i), view(Av,:,:,i), w, 2, init)
        end
    end
    return R
end

function _wsumN!(R::StridedArray{T}, A::DenseArray{T,N}, w::StridedVector{T}, dim::Int, init::Bool) where {T<:BlasReal,N}
    @assert N >= 3
    if dim <= 2
        m = size(A, 1)
        n = size(A, 2)
        npages = 1
        for i = 3:N
            npages *= size(A, i)
        end
        rlen = ifelse(dim == 1, n, m)
        Rv = reshape(R, (rlen, npages))
        for i = 1:npages
            _wsum2_blas!(view(Rv,:,i), view(A,:,:,i), w, dim, init)
        end
    else
        _wsum_general!(R, identity, A, w, dim, init)
    end
    return R
end

# General Cartesian-based weighted sum across dimensions
@generated function _wsum_general!(R::AbstractArray{RT}, f::supertype(typeof(abs)),
                                   A::AbstractArray{T,N}, w::AbstractVector{WT}, dim::Int, init::Bool) where {T,RT,WT,N}
    quote
        init && fill!(R, zero(RT))
        wi = zero(WT)
        if dim == 1
            @nextract $N sizeR d->size(R,d)
            sizA1 = size(A, 1)
            @nloops $N i d->(d>1 ? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin
                @inbounds r = (@nref $N R j)
                for i_1 = 1:sizA1
                    @inbounds r += f(@nref $N A i) * w[i_1]
                end
                @inbounds (@nref $N R j) = r
            end
        else
            @nloops $N i A d->(if d == dim
                                   wi = w[i_d]
                                   j_d = 1
                               else
                                   j_d = i_d
                               end) @inbounds (@nref $N R j) += f(@nref $N A i) * wi
        end
        return R
    end
end

@generated function _wsum_centralize!(R::AbstractArray{RT}, f::supertype(typeof(abs)),
                                      A::AbstractArray{T,N}, w::AbstractVector{WT}, means,
                                      dim::Int, init::Bool) where {T,RT,WT,N}
    quote
        init && fill!(R, zero(RT))
        wi = zero(WT)
        if dim == 1
            @nextract $N sizeR d->size(R,d)
            sizA1 = size(A, 1)
            @nloops $N i d->(d>1 ? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin
                @inbounds r = (@nref $N R j)
                @inbounds m = (@nref $N means j)
                for i_1 = 1:sizA1
                    @inbounds r += f((@nref $N A i) - m) * w[i_1]
                end
                @inbounds (@nref $N R j) = r
            end
        else
            @nloops $N i A d->(if d == dim
                                   wi = w[i_d]
                                   j_d = 1
                               else
                                   j_d = i_d
                               end) @inbounds (@nref $N R j) += f((@nref $N A i) - (@nref $N means j)) * wi
        end
        return R
    end
end

# N = 1
_wsum!(R::StridedArray{T}, A::DenseArray{T,1}, w::StridedVector{T}, dim::Int, init::Bool) where {T<:BlasReal} =
    _wsum1!(R, A, w, init)

# N = 2
_wsum!(R::StridedArray{T}, A::DenseArray{T,2}, w::StridedVector{T}, dim::Int, init::Bool) where {T<:BlasReal} =
    (_wsum2_blas!(view(R,:), A, w, dim, init); R)

# N >= 3
_wsum!(R::StridedArray{T}, A::DenseArray{T,N}, w::StridedVector{T}, dim::Int, init::Bool) where {T<:BlasReal,N} =
    _wsumN!(R, A, w, dim, init)

_wsum!(R::AbstractArray, A::AbstractArray, w::AbstractVector, dim::Int, init::Bool) =
    _wsum_general!(R, identity, A, w, dim, init)

## wsum! and wsum

wsumtype(::Type{T}, ::Type{W}) where {T,W} = typeof(zero(T) * zero(W) + zero(T) * zero(W))
wsumtype(::Type{T}, ::Type{T}) where {T<:BlasReal} = T

"""
    wsum!(R, A, w, dim; init=true)

Compute the weighted sum of `A` with weights `w` over the dimension `dim` and store
the result in `R`. If `init=false`, the sum is added to `R` rather than starting
from zero.
"""
function wsum!(R::AbstractArray, A::AbstractArray{T,N}, w::AbstractVector, dim::Int; init::Bool=true) where {T,N}
    1 <= dim <= N || error("dim should be within [1, $N]")
    ndims(R) <= N || error("ndims(R) should not exceed $N")
    length(w) == size(A,dim) || throw(DimensionMismatch("Inconsistent array dimension."))
    # TODO: more careful examination of R's size
    _wsum!(R, A, w, dim, init)
end

function wsum(A::AbstractArray{T}, w::AbstractVector{W}, dim::Int) where {T<:Number,W<:Real}
    length(w) == size(A,dim) || throw(DimensionMismatch("Inconsistent array dimension."))
    _wsum!(similar(A, wsumtype(T,W), Base.reduced_indices(axes(A), dim)), A, w, dim, true)
end

# extended sum! and wsum

Base.sum!(R::AbstractArray, A::AbstractArray, w::AbstractWeights{<:Real}, dim::Int; init::Bool=true) =
    wsum!(R, A, values(w), dim; init=init)

Base.sum(A::AbstractArray{<:Number}, w::AbstractWeights{<:Real}, dim::Int) = wsum(A, values(w), dim)

###### Weighted means #####

"""
    wmean(v, w::AbstractVector)

Compute the weighted mean of an array `v` with weights `w`.
"""
function wmean(v::AbstractArray{<:Number}, w::AbstractVector)
    Base.depwarn("wmean is deprecated, use mean(v, weights(w)) instead.", :wmean)
    mean(v, weights(w))
end

"""
    mean(A::AbstractArray, w::AbstractWeights[, dim::Int])

Compute the weighted mean of array `A` with weight vector `w`
(of type `AbstractWeights`). If `dim` is provided, compute the
weighted mean along dimension `dim`.

# Examples
```julia
w = rand(n)
mean(x, weights(w))
```
"""
mean(A::AbstractArray, w::AbstractWeights) = sum(A, w) / sum(w)

"""
    mean(R::AbstractArray, , A::AbstractArray, w::AbstractWeights[, dim::Int])

Compute the weighted mean of array `A` with weight vector `w`
(of type `AbstractWeights`) along dimension `dim`, and write results to `R`.
"""
mean!(R::AbstractArray, A::AbstractArray, w::AbstractWeights, dim::Int) =
    rmul!(Base.sum!(R, A, w, dim), inv(sum(w)))

wmeantype(::Type{T}, ::Type{W}) where {T,W} = typeof((zero(T)*zero(W) + zero(T)*zero(W)) / one(W))
wmeantype(::Type{T}, ::Type{T}) where {T<:BlasReal} = T

mean(A::AbstractArray{T}, w::AbstractWeights{W}, dim::Int) where {T<:Number,W<:Real} =
    mean!(similar(A, wmeantype(T, W), Base.reduced_indices(axes(A), dim)), A, w, dim)

###### Weighted median #####
function median(v::AbstractArray, w::AbstractWeights)
    throw(MethodError(median, (v, w)))
end

"""
    median(v::RealVector, w::AbstractWeights)

Compute the weighted median of `x`, using weights given by a weight vector `w`
(of type `AbstractWeights`). The weight and data vectors must have the same length.

The weighted median ``x_k`` is the element of `x` that satisfies
``\\sum_{x_i < x_k} w_i \\le \\frac{1}{2} \\sum_{j} w_j`` and
``\\sum_{x_i > x_k} w_i \\le \\frac{1}{2} \\sum_{j} w_j``.

If a weight has value zero, then its associated data point is ignored.
If none of the weights are positive, an error is thrown.
`NaN` is returned if `x` contains any `NaN` values.
An error is raised if `w` contains any `NaN` values.
"""
function median(v::RealVector, w::AbstractWeights{<:Real})
    isempty(v) && error("median of an empty array is undefined")
    if length(v) != length(w)
        error("data and weight vectors must be the same size")
    end
    @inbounds for x in w.values
        isnan(x) && error("weight vector cannot contain NaN entries")
    end
    @inbounds for x in v
        isnan(x) && return x
    end
    mask = w.values .!= 0
    if !any(mask)
        error("all weights are zero")
    end
    if all(w.values .<= 0)
        error("no positive weights found")
    end
    v = v[mask]
    wt = w[mask]
    midpoint = w.sum / 2
    maxval, maxind = findmax(wt)
    if maxval > midpoint
        v[maxind]
    else
        permute = sortperm(v)
        cumulative_weight = zero(eltype(wt))
        i = 0
        for (_i, p) in enumerate(permute)
            i = _i
            if cumulative_weight == midpoint
                i += 1
                break
            elseif cumulative_weight > midpoint
                cumulative_weight -= wt[p]
                break
            end
            cumulative_weight += wt[p]
        end
        if cumulative_weight == midpoint
            middle(v[permute[i-2]], v[permute[i-1]])
        else
            middle(v[permute[i-1]])
        end
    end
end

"""
    wmedian(v, w)

Compute the weighted median of an array `v` with weights `w`, given as either a
vector or an `AbstractWeights` vector.
"""
wmedian(v::RealVector, w::RealVector) = median(v, weights(w))
wmedian(v::RealVector, w::AbstractWeights{<:Real}) = median(v, w)

###### Weighted quantile #####

"""
    quantile(v, w::AbstractWeights, p)

Compute the weighted quantiles of a vector `v` at a specified set of probability
values `p`, using weights given by a weight vector `w` (of type `AbstractWeights`).
Weights must not be negative. The weights and data vectors must have the same length.

With [`FrequencyWeights`](@ref), the function returns the same result as
`quantile` for a vector with repeated values.
With non `FrequencyWeights`,  denote ``N`` the length of the vector, ``w`` the vector of weights,
``h = p (\\sum_{i<= N}w_i - w_1) + w_1`` the cumulative weight corresponding to the
probability ``p`` and ``S_k = \\sum_{i<=k}w_i`` the cumulative weight for each
observation, define ``v_{k+1}`` the smallest element of `v` such that ``S_{k+1}``
is strictly superior to ``h``. The weighted ``p`` quantile is given by ``v_k + \\gamma (v_{k+1} -v_k)``
with  ``\\gamma = (h - S_k)/(S_{k+1}-S_k)``. In particular, when `w` is a vector
of ones, the function returns the same result as `quantile`.
"""
function quantile(v::RealVector{V}, w::AbstractWeights{W}, p::RealVector) where {V,W<:Real}
    # checks
    isempty(v) && error("quantile of an empty array is undefined")
    isempty(p) && throw(ArgumentError("empty quantile array"))

    w.sum == 0 && error("weight vector cannot sum to zero")
    length(v) == length(w) || error("data and weight vectors must be the same size, got $(length(v)) and $(length(w))")
    for x in w.values
        isnan(x) && error("weight vector cannot contain NaN entries")
        x < 0 && error("weight vector cannot contain negative entries")
    end

    # remove zeros weights and sort
    wsum = sum(w)
    nz = .!iszero.(w)
    vw = sort!(collect(zip(view(v, nz), view(w, nz))))
    N = length(vw)

    # prepare percentiles
    ppermute = sortperm(p)
    p = p[ppermute]
    p = bound_quantiles(p)

    # prepare out vector
    out = Vector{typeof(zero(V)/1)}(undef, length(p))
    fill!(out, vw[end][1])

    # start looping on quantiles
    Sk, Skold = zero(W), zero(W)
    vk, vkold = zero(V), zero(V)
    k = 0

    for i in 1:length(p)
        if isa(w, FrequencyWeights)
            h = p[i] * (wsum - 1) + 1
        else
            h = p[i] * (wsum - vw[1][2]) + vw[1][2]
        end
        while Sk <= h
            k += 1
            if k > N
               # out was initialized with maximum v
               return out
            end
            Skold, vkold = Sk, vk
            vk, wk = vw[k]
            Sk += wk
        end
        if isa(w, FrequencyWeights)
            out[ppermute[i]] = vkold + min(h - Skold, 1) * (vk - vkold)
        else
            out[ppermute[i]] = vkold + (h - Skold) / (Sk - Skold) * (vk - vkold)
        end
    end
    return out
end

# similarly to statistics.jl in Base
function bound_quantiles(qs::AbstractVector{T}) where T<:Real
    epsilon = 100 * eps()
    if (any(qs .< -epsilon) || any(qs .> 1+epsilon))
        throw(ArgumentError("quantiles out of [0,1] range"))
    end
    T[min(one(T), max(zero(T), q)) for q = qs]
end

quantile(v::RealVector, w::AbstractWeights{<:Real}, p::Number) = quantile(v, w, [p])[1]

"""
    wquantile(v, w, p)

Compute the `p`th quantile(s) of `v` with weights `w`, given as either a vector
or an `AbstractWeights` vector.
"""
wquantile(v::RealVector, w::AbstractWeights{<:Real}, p::RealVector) = quantile(v, w, p)
wquantile(v::RealVector, w::AbstractWeights{<:Real}, p::Number) = quantile(v, w, [p])[1]
wquantile(v::RealVector, w::RealVector, p::RealVector) = quantile(v, weights(w), p)
wquantile(v::RealVector, w::RealVector, p::Number) = quantile(v, weights(w), [p])[1]


text:
# test_tpch.jl

@testset "Test TPC-H" begin
    sf = "0.01"

    # load TPC-H into DuckDB
    native_con = DBInterface.connect(DuckDB.DB)
    DBInterface.execute(native_con, "CALL dbgen(sf=$sf)")

    # convert all tables to Julia DataFrames
    customer = DataFrame(DBInterface.execute(native_con, "SELECT * FROM customer"))
    lineitem = DataFrame(DBInterface.execute(native_con, "SELECT * FROM lineitem"))
    nation = DataFrame(DBInterface.execute(native_con, "SELECT * FROM nation"))
    orders = DataFrame(DBInterface.execute(native_con, "SELECT * FROM orders"))
    part = DataFrame(DBInterface.execute(native_con, "SELECT * FROM part"))
    partsupp = DataFrame(DBInterface.execute(native_con, "SELECT * FROM partsupp"))
    region = DataFrame(DBInterface.execute(native_con, "SELECT * FROM region"))
    supplier = DataFrame(DBInterface.execute(native_con, "SELECT * FROM supplier"))

    # now open a new in-memory database, and register the dataframes there
    df_con = DBInterface.connect(DuckDB.DB)
    DuckDB.register_data_frame(df_con, customer, "customer")
    DuckDB.register_data_frame(df_con, lineitem, "lineitem")
    DuckDB.register_data_frame(df_con, nation, "nation")
    DuckDB.register_data_frame(df_con, orders, "orders")
    DuckDB.register_data_frame(df_con, part, "part")
    DuckDB.register_data_frame(df_con, partsupp, "partsupp")
    DuckDB.register_data_frame(df_con, region, "region")
    DuckDB.register_data_frame(df_con, supplier, "supplier")
    GC.gc()

    # run all the queries
    for i in 1:22
        #         print("Q$i\n")
        # for each query, compare the results of the query ran on the original tables
        # versus the result when run on the Julia DataFrames
        res = DataFrame(DBInterface.execute(df_con, "PRAGMA tpch($i)"))
        res2 = DataFrame(DBInterface.execute(native_con, "PRAGMA tpch($i)"))
        @test isequal(res, res2)
        #         print("Native DuckDB\n")
        #         @time begin
        #             results = DBInterface.execute(native_con, "PRAGMA tpch($i)")
        #         end
        #         print("DataFrame\n")
        #         @time begin
        #             results = DBInterface.execute(df_con, "PRAGMA tpch($i)")
        #         end
    end

    DBInterface.close!(df_con)
    DBInterface.close!(native_con)
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
"""
    GeoRegion(RegID::AbstractString) -> geo::GeoRegion

Extracts information of the GeoRegion with the ID `RegID`.  If no GeoRegion with this ID exists, an error is thrown.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the GeoRegion.
            If the ID is not valid (i.e. not being used), then an error will be thrown.

Returns
=======

- `geo` : A GeoRegion
"""
function GeoRegion(
    RegID::AbstractString,
    ST = String,
    FT = Float64
)

    @info "$(now()) - GeoRegions.jl - Retrieving information for the GeoRegion defined by the ID $RegID"

    regvec,filevec,typevec = listGeoRegions(); isgeoregion(RegID,regvec)
    ind = findall(RegID.==regvec)[1]
    return getgeoregion(
        RegID,
        joinpath(DEPOT_PATH[1],"files","GeoRegions",filevec[ind]),
        typevec[ind],
        ST,FT
    )

end

function getgeoregion(
    RegID :: AbstractString,
    fname :: AbstractString,
    GType :: AbstractString,
    ST = String,
    FT = Float64
)

    if GType == "PolyRegion"
          return PolyRegion(RegID,fname,ST,FT)
    else; return RectRegion(RegID,fname,ST,FT)
    end

end

function PolyRegion(
    RegID :: AbstractString,
    fname :: AbstractString,
    ST = String,
    FT = Float64
)

    rvec = listpolyregions(fname)
    ind  = findall(rvec.==RegID)[1]
    ind  = (ind) * 4 + 1

    flines = readlines(fname)
    IDinfo = flines[ind];   ParID,RegName = getpolymeta(IDinfo)
    regIDX = flines[ind+1]; X = getpolyx(regIDX)
    regIDY = flines[ind+2]; Y = getpolyy(regIDY)

    is180,is360 = checkbounds(maximum(Y),minimum(Y),maximum(X),minimum(X))

    return PolyRegion{ST,FT}(
        RegID,ParID,RegName,
        maximum(Y),minimum(Y),maximum(X),minimum(X),Point2.(X,Y),
        is180,is360
    )

end

function RectRegion(
    RegID :: AbstractString,
    fname :: AbstractString,
    ST = String,
    FT = Float64
)

    rvec = listrectregions(fname)
    ind  = findall(rvec.==RegID)[1]

    IDinfo = readdlm(fname,',',comments=true,comment_char='#')[ind,:]
    ParID,RegName,regN,regS,regE,regW = IDinfo[[2,7,3,5,6,4]]
    ParID = replace(ParID," "=>"")
    RegName = replace(RegName," "=>"")
    RegName = replace(RegName,"-"=>" ")

    is180,is360 = checkbounds(regN,regS,regE,regW)

    return RectRegion{ST,FT}(
        RegID,ParID,RegName,
        regN,regS,regE,regW,
        is180,is360
    )

end

"""
    templateGeoRegions(; path::AbstractString, overwrite::Bool=false)

Copies the template files for defining GeoRegions in textfiles, that can then be added as a batch using addGeoRegions()

Keyword Arguments
=================

- `path`      : The folder to copy the files to
- `overwrite` : If template files exist in this folder, overwrite?
"""
function templateGeoRegions(;
    path      :: AbstractString=pwd(),
    overwrite :: Bool=false
)

    if !isdir(path); mkpath(path) end
    for fname in ["recttemplate.txt","polytemplate.txt"]

        ftem = joinpath(@__DIR__,"..","extra",fname)
        freg = joinpath(path,fname)

        if !overwrite
            if !isfile(freg)

                @debug "$(now()) - GeoRegions.jl - Unable to find $freg, copying data from $ftem ..."

                open(freg,"w") do io
                    open(ftem) do f
                        for line in readlines(f)
                            write(io,"$line\n")
                        end
                    end
                end

            end
        else

            if isfile(freg)
                @warn "$(now()) - GeoRegions.jl - Overwriting $freg with original file in $ftem ..."
                rm(freg,force=true)
            end

            open(freg,"w") do io
                open(ftem) do f
                    for line in readlines(f)
                        write(io,"$line\n")
                    end
                end
            end

        end

    end

    return nothing

end

"""
    listGeoRegions()

List all GeoRegions, files the data are stored in, and the `Type` of `GeoRegion` they are.

Output
======

- `geovec`  : List of all the GeoRegion IDs
- `filevec` : List of the files that the GeoRegion information is stored in
- `typevec` : List of the `Type` of the `GeoRegion` corresponding to the `geovec` ID list
"""
function listGeoRegions()

    flist   = ["rectlist.txt","polylist.txt","giorgi.txt","srex.txt","ar6.txt"]
    regvec  = []
    filevec = []
    typevec = []

    for fname in flist
        copygeoregions(fname)
        rvec,rtype = listgeoregions(joinpath(DEPOT_PATH[1],"files","GeoRegions",fname))
        regvec = vcat(regvec,rvec)
        nreg = length(rvec)
        fvec = fill(fname,nreg); filevec = vcat(filevec,fvec)
        tvec = fill(rtype,nreg); typevec = vcat(typevec,tvec)
    end

    return regvec,filevec,typevec

end

"""
    resetGeoRegions(; allfiles::Bool = false)

Reset all the files containing GeoRegion information back to the default.

Arguments
=========

- `allfiles` : If `true`, reset the GeoRegions defined in Giorgi & Francisco [2000], AR6 Regions (Iturbide et al., 2020; ESSD) and Seneviratne et al. [2012] as well.  If `false`, only reset the custom GeoRegions.
"""
function resetGeoRegions(;allfiles=false)

    if allfiles

        @info "$(now()) - GeoRegions.jl - Resetting both the master and custom lists of GeoRegions back to the default"
        flist = ["rectlist.txt","polylist.txt","giorgi.txt","srex.txt","ar6.txt"]
    else
        @info "$(now()) - GeoRegions.jl - Resetting the custom lists of GeoRegions back to the default"
        flist = ["rectlist.txt","polylist.txt"]
    end

    for fname in flist
        copygeoregions(fname,overwrite=true)
    end

    return

end

"""
    addGeoRegions(fname::AbstractString)

Extracts information of the GeoRegion with the ID `RegID`.  If no GeoRegion with this ID exists, an error is thrown.

Arguments
=========

- `fname` : name + path of the file containing GeoRegion information
"""
function addGeoRegions(fname::AbstractString)

    @info "$(now()) - GeoRegions.jl - Importing user-defined GeoRegions from the file $fname directly into the custom lists"

    rvec,rtype = listgeoregions(fname)
    for reg in rvec
        if !isGeoRegion(reg,throw=false)
            g = getgeoregion(reg,fname,rtype)
            if rtype == "PolyRegion"
                  _,_,lon,lat = coordGeoRegion(g)
                  PolyRegion(g.regID,g.parID,g.name,lon,lat)
            else; RectRegion(g.regID,g.parID,g.name,[g.N,g.S,g.E,g.W])
            end
        else
            @warn "$(now()) - GeoRegions.jl - The GeoRegion ID $reg is already in use. Please use a different ID, or you can remove the ID using removeGeoRegion()."
        end
    end

    return nothing

end

"""
    isGeoRegion(
        RegID :: AbstractString;
        throw :: Bool = true
    ) -> tf :: Bool

Extracts information of the GeoRegion with the ID `RegID`.  If no GeoRegion with this ID exists, an error is thrown.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the GeoRegion.
        If the ID is not valid (i.e. not being used), then an error will be thrown.
- `throw` : If `true`, then throws an error if `RegID` is not a valid `GeoRegion` identifier instead of returning the Boolean `tf`

Returns
=======

- `tf` : True / False
"""
function isGeoRegion(RegID::AbstractString;throw::Bool=true)

    regvec,filevec,typevec = listGeoRegions()
    return isgeoregion(RegID,regvec;throw=throw)

end

function getpolymeta(metastring::AbstractString)

    metastring = replace(metastring," "=>"")
    metastring = replace(metastring,"RegID"=>"")
    metastring = replace(metastring,"-"=>" ")
    return split(metastring,",")[2:3]

end

function getpolyx(xstring::AbstractString)

    xstring = replace(xstring," "=>"")
    xstring = replace(xstring,"RegX"=>"")
    return parse.(Float64,split(xstring,","))

end

function getpolyy(ystring::AbstractString)

    ystring = replace(ystring," "=>"")
    ystring = replace(ystring,"RegY"=>"")
    return parse.(Float64,split(ystring,","))

end

function checkbounds(
    regN :: Real,
    regS :: Real,
    regE :: Real,
    regW :: Real
)

    if (regN>90) || (regN<-90)
        error("$(now()) - GeoRegions.jl - The latitude of the GeoRegion's northern bound at $regN is not valid.")
    end

    if (regS>90) || (regS<-90)
        error("$(now()) - GeoRegions.jl - The latitude of the GeoRegion's southern bound at $regS is not valid.")
    end

    if (regE>360) || (regE<-180)
        error("$(now()) - GeoRegions.jl - The longitude of the GeoRegion's eastern bound at $regE is not valid.")
    end

    if (regW>360) || (regW<-180)
        error("$(now()) - GeoRegions.jl - The longitude of the GeoRegion's western bound at $regW is not valid.")
    end

    if (regE - regW) > 360
        error("$(now()) - GeoRegions.jl - The GeoRegion cannot be more than 360º in Longitude.")
    end

    if regE < regW
        error("$(now()) - GeoRegions.jl - The eastern bound of the GeoRegion cannot be west of the western bound.")
    end

    if regN < regS
        error("$(now()) - GeoRegions.jl - The northern bound of the GeoRegion cannot be south of the southern bound.")
    end

    if regE > 180; is360 = true; else is360 = false end
    if regW < 0;   is180 = true; else is180 = false end
    if !(is180) && !(is360); is360 = true end

    return is180,is360

end

function copygeoregions(
    fname::AbstractString;
    overwrite::Bool=false
)

    jfol = joinpath(DEPOT_PATH[1],"files","GeoRegions"); mkpath(jfol);
    ftem = joinpath(@__DIR__,"..","extra",fname)
    freg = joinpath(jfol,fname)

    if !overwrite
        if !isfile(freg)

            @debug "$(now()) - GeoRegions.jl - Unable to find $freg, copying data from $ftem ..."

            open(freg,"w") do io
                open(ftem) do f
                    for line in readlines(f)
                        write(io,"$line\n")
                    end
                end
            end

        end
    else

        if isfile(freg)
            @warn "$(now()) - GeoRegions.jl - Overwriting $freg with original file in $ftem ..."
            rm(freg,force=true)
        end

        open(freg,"w") do io
            open(ftem) do f
                for line in readlines(f)
                    write(io,"$line\n")
                end
            end
        end

    end

    return

end

function listgeoregions(fname::AbstractString)

    gtype = readlines(fname)[1]
    if occursin("PolyRegion",gtype)
          return listpolyregions(fname), "PolyRegion"
    else; return listrectregions(fname), "RectRegion"
    end

end

function listpolyregions(fname::AbstractString)

      gcount = -1
      for l in eachline(fname)
          if occursin("RegID",l); gcount += 1 end
      end

      regvec = Vector{String}(undef,gcount)
      gcount = -1
      for l in eachline(fname)
          if occursin("RegID",l)
              gcount += 1
              if gcount > 0
                  l = replace(l," "=>"")
                  l = replace(l,"RegID"=>"")
                  regvec[gcount] = split(l,",")[1]
              end
          end
      end

      return regvec

end

function listrectregions(fname::AbstractString)

      return readdlm(fname,',',comments=true,comment_char='#')[:,1]

end

function isgeoregion(RegID::AbstractString,regvec::AbstractArray;throw::Bool=true)

    @info "$(now()) - GeoRegions.jl - Checking to see if the ID $RegID is in use"

    if sum(regvec.==RegID) == 0
        if throw
            error("$(now()) - GeoRegions.jl - $(RegID) is not a valid GeoRegion identifier, use either RectRegion() or PolyRegion() to add this GeoRegion to the list.")
        else
            @warn "$(now()) - GeoRegions.jl - $(RegID) is not a valid GeoRegion identifier, use either RectRegion() or PolyRegion() to add this GeoRegion to the list."
            return false
        end
    else
        @info "$(now()) - GeoRegions.jl - The ID $RegID is already in use"
        return true
    end

end


text:
"""
测试refine Surface
"""

using Test
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Surface

#@testset "靠近费米面" begin
#    qs = common_square_lattice(0.20)
#    Γ4 = TFGamma4(qs, 8.0, 16, 20)
#    surf = const_energy_triangle(Γ4.ltris, 0.5, Γ4.model.dispersion[1])
#    plt = draw_polygon(surf; pcidx=ones(Int64, length(surf)))
#    @savepng plt "surf_refine2"
#    surf2 = refine_to_surface(surf, 0.5, Γ4.model.dispersion[1])
#    plt = draw_polygon(surf2; pcidx=ones(Int64, length(surf2)))
#    @savepng plt "surf_refine22"
#    surf3 = refine_to_surface(surf2, 0.5, Γ4.model.dispersion[1])
#    plt = draw_polygon(surf3; pcidx=ones(Int64, length(surf3)))
#    @savepng plt "surf_refine23"
#end
#
#
#
#@testset "细化费米面" begin
#    qs = common_square_lattice(0.20)
#    Γ4 = TFGamma4(qs, 8.0, 16, 20)
#    surf = const_energy_triangle(Γ4.ltris, 0.5, Γ4.model.dispersion[1])
#    plt = draw_points([tri.center for tri in surf])
#    @savepng plt "surf_refine24"
#    surf2 = refine_list_triangle(surf, 1)
#    plt = draw_points([tri.center for tri in surf2])
#    @savepng plt "surf_refine25"
#    surf3 = refine_list_triangle(surf2, 3)
#    plt = draw_points([tri.center for tri in surf3])
#    @savepng plt "surf_refine26"
#end
#
#
#@testset "addition" begin
#    qs = common_square_lattice(0.20)
#    Γ4 = TFGamma4(qs, 8.0, 16, 100)
#    posi, nega, away = TFGamma4_addition_ltris_mt(Γ4, 2.0)
#    plt = draw_points([tri.center for tri in away])
#    draw_points!(plt, [tri.center for tri in posi])
#    draw_points!(plt, [tri.center for tri in nega])
#    @savepng plt "surf_refine27"
#    posi, nega, away = TFGamma4_addition_ltris_mt(Γ4, 4.0)
#    plt = draw_points([tri.center for tri in away])
#    draw_points!(plt, [tri.center for tri in posi])
#    draw_points!(plt, [tri.center for tri in nega])
#    @savepng plt "surf_refine28"
#    posi, nega, away = TFGamma4_addition_ltris_mt(Γ4, 8.0)
#    plt = draw_points([tri.center for tri in away])
#    draw_points!(plt, [tri.center for tri in posi])
#    draw_points!(plt, [tri.center for tri in nega])
#    @savepng plt "surf_refine29"
#    pt1 = posi[1].center
#    pt2 = posi[2].center
#    println(Γ4.model.dispersion[1](pt1.x, pt1.y))
#    println(Γ4.model.dispersion[1](pt2.x, pt2.y))
#end

qs = common_square_lattice(0.20)
Γ4 = TFGamma4(qs, 8.0, 16, 100)
nΓ4 = TFGamma4_addition_ltris_mt(Γ4, 4.0)
plt = draw_points([tri.center for tri in nΓ4.ltris])
@savepng plt "surf_refine210"


text:
@testset "FBM" begin
    rng = MersenneTwister(42)
    h = 0.3
    k = FBMKernel(; h=h)
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    @test k(v1, v2) ≈
          (
        sqeuclidean(v1, zero(v1))^h + sqeuclidean(v2, zero(v2))^h -
        sqeuclidean(v1 - v2, zero(v1 - v2))^h
    ) / 2 atol = 1e-5
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"

    test_interface(k)
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"
    test_ADs(FBMKernel; ADs=[:ReverseDiff, :Zygote])
    @test_broken "Tests failing for kernelmatrix(k, x) for ForwardDiff"
    test_params(k, ([h],))
end


text:
# ------------------------------------------------------------------------------------------
# ## Benchmarking tips
#
# In the `Julia is fast` notebook, we saw the package `BenchmarkTools` and used its
# `@benchmark` macro.
#
# In this notebook, we'll explore the importance of "interpolating" global variables when
# benchmarking functions.
#
# We interpolate a global variable by throwing a `$` in front of it. For example, in `Julia
# is fast`, we benchmarked the `sum` function using `Vector` `A` via
#
# ```julia
# @benchmark sum($A)
# ```
#
# not
#
# ```julia
# @benchmark sum(A)
# ```
#
# Let's see if this can make a difference by examining the ratio in execution times of
# `sum($A)` and `sum(A)` for differently sized arrays `A`.
#
# #### Exercise
#
# Call the `sum` function on a pseudo-randomly populated 1D array called `foo` of several
# lengths between 2 and 2^20 (~10^6). For each size of `foo`, determine the ratio of
# execution times for `sum(foo)` and `sum($foo)`. (To determine this ratio, use the minimum
# run times in each case.)
#
# Plot the ratio of execution times for non-interpolated and interpolated `foo` in calls to
# `sum` versus the length of `foo`. Does interpolating `foo` seem to matter? If so, for what
# sizes of `foo`?
#
# #### Solution
#
#
# ------------------------------------------------------------------------------------------

using BenchmarkTools, Plots

gr()

store_ratios = []
foo = rand()
for i in 1:20
    foo = rand(2^i)
    not_interpolated = @benchmark sum(foo)
    interpolated = @benchmark sum($foo)
    min_not_interpolated = minimum(not_interpolated.times)
    min_interpolated = minimum(interpolated.times)
    ratio = min_not_interpolated / min_interpolated
    append!(store_ratios, ratio)
end

plot([i for i in 1:20], store_ratios, legend=false)
scatter!([i for i in 1:20], store_ratios)
xlabel!("The log of the length of `foo`")
ylabel!("Factor of performance increase")
title!("Performance increase from interpolating a global variable")

# ------------------------------------------------------------------------------------------
# In this case, interpolating the global variable is more important when that global
# variable is bound to a relatively small array. For arrays with a few thousand elements and
# more, interpolation doesn't seem to have an impact.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Performance tips -- type stability
#
# One way to optimize code in Julia is to ensure **type stability**. If the type(s) of some
# variables in a function are subject to change or ambiguity, the compiler cannot reason as
# well about those variables, and performance will take a hit. Conversely, we allow the
# compiler to optimize and generate more efficient machine code when we declare variables so
# that their types will be fixed throughout the function body.
#
# For example, let's say we had functions called `baz` and `bar` with the following
# definitions
#
# ```julia
# function baz()
#     s = rand()
#     if s > 2/3
#         return .666667
#     elseif s > 1/3
#         return 1//3
#     else
#         return 0
#     end
# end
# ```
#
# ```julia
# function bar()
#     s = rand()
#     if s > 2/3
#         return .666667
#     elseif s > 1/3
#         return .3333333
#     else
#         return 0.0
#     end
# end
# ```
#
# When I benchmark these via
#
# ```julia
# using BenchmarkTools
# @benchmark baz()
# @benchmark bar()
# ```
#
# I see that `bar` is almost three times as fast as `baz`! The reason is that `bar` is type
# stable while `baz` is not: the compiler can tell that `bar` will always return a `Float`,
# whereas `baz` could return a `Float`, an `Int`, or a `Rational`. When the compiler can
# tell what the types of outputs from a function, or variables declared *within a function*
# are without running the code, it can do much better.
#
# #### Exercise
#
# The following definition for `my_sum` is not type stable.
#
# ```julia
# function my_sum(A)
#     output = 0
#     for x in A
#         output += x
#     end
#     return output
# end
# ```
#
# Copy and execute the above code into a new cell. Benchmark it using `A = rand(10^3)`. Then
# write a new function called `my_sum2` with the same function body as `my_sum`. Update
# `my_sum2` to make it type stable, and benchmark it for a randomly populated array with
# 10^3 entries.
#
# How much does type stability impact performance? If you'd like, try this same exercise for
# multiple sizes of `A` to see if this changes your answer!
#
# #### Solution
# ------------------------------------------------------------------------------------------

function my_sum(A)
    output = 0
    for x in A
        output += x
    end
    return output
end

using BenchmarkTools

A = rand(10^3)
@benchmark my_sum($A)

function my_sum2(A)
    output = zero(eltype(A))
    # or:
    # output = 0.0
    for x in A
        output += x
    end
    return output
end

@benchmark my_sum2($A)

# ------------------------------------------------------------------------------------------
# In my run, `my_sum2` outperformed `my_sum` by a factor of over 20 for an array with 1000
# elements!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise
#
# Make the following code type stable. You'll know your efforts are paying off when you see
# a performance boost!
#
# ```julia
# """
#     my_sqrt(x)
#
# Calculate the square root of `x` with Newton's method.
# """
# function my_sqrt(x)
#     output = 1
#     for i in 1:1000
#         output = .5 * (output + x/output)
#     end
#     output
# end
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Solution
#
# You need to make sure `ouptut`'s type never changes in order to make this function type
# stable!
# ------------------------------------------------------------------------------------------

# The type unstable version
function my_sqrt(x)
    output = 1
    for i in 1:1000
        output = .5 * (output + x/output)
    end
    output
end

@benchmark my_sqrt(100)

# The type stable version, where `output` is always a floating point number.
function my_sqrt(x)
    output = 1.0
    for i in 1:1000
        output = .5 * (output + x/output)
    end
    output
end

@benchmark my_sqrt(100)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using AutoTypeParameters: freeze, thaw

using Base.Test
using AutoHashEquals
using Compat

immutable Type1
    a
    b
end

@auto_hash_equals immutable Type2{A,B}
    a::A
    b::B
end

for (value, param) in [(1, 1),
                       (Int64, Int64),
                       ((1, 2, 3), (1, 2, 3)),
                       (Type1(1, 2), symbol("ATP Type1(1,2)")),
                       ("abc", symbol("ATP \"abc\"")),
                       ([1,2], symbol("ATP [1,2]")),
                       ((1, ("two", 3.0)), symbol("ATP (1,(\"two\",3.0))")),
                       (Type2{Symbol,Int}(:a, 1), symbol("ATP Type2{Symbol,Int64}(:a,1)"))]

    formats = VERSION < v"0.4-" ? (:show,) : (:show, :serialize)

    for format in formats

        println("$(format): $(value) -> $(freeze(value; format=format))")

        # we can invert things correctly
        @test thaw(eval, freeze(value; format=format)) == value

        # what we get is a valid type parameter
        @compat @test Val{freeze(value)} == Val{param}

        if format == :show
            # we get what we expect
            @test freeze(value) == param
        end

    end

end


text:
# ---
# title: Wide data
# cover: assets/wide_data.png
# description: Working with data in the wide format.
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
using AlgebraOfGraphics: density
set_aog_theme!() #src

df = (a=randn(100), b=randn(100), c=randn(100))
labels = ["Trace 1", "Trace 2", "Trace 3"]
plt = data(df) *
    density() *
    mapping([:a, :b, :c] .=> "some label") *
    mapping(color=dims(1) => renamer(labels))
draw(plt)

#

df = (a=rand(100), b=rand(100), c=rand(100), d=rand(100))
labels = ["Trace One", "Trace Two", "Trace Three"]
layers = linear() + visual(Scatter)
plt = data(df) * layers * mapping(1, 2:4 .=> "value", color=dims(1) => renamer(labels))
draw(plt)

# The wide format is combined with broadcast semantics.

# ## Axes are linked when they correspond to the same variable

df = (
    sepal_length = 1 .+ rand(100),
    sepal_width = 2 .+ rand(100),
    petal_length = 3 .+ rand(100),
    petal_width = 4 .+ rand(100)
)
xvars = ["sepal_length", "sepal_width"]
yvars = ["petal_length" "petal_width"]
layers = linear() + visual(Scatter)
plt = data(df) * layers * mapping(xvars, yvars, col=dims(1), row=dims(2))
draw(plt)

# ## Axes can be fully linked or fully unlinked

draw(plt, facet = (; linkxaxes = :all, linkyaxes = :all))

#

fg = draw(plt, facet = (; linkxaxes = :none, linkyaxes = :none))

# save cover image #src
mkpath("assets") #src
save("assets/wide_data.png", fg) #src


text:
# -*- coding: utf-8 -*-
adparam = SNN.ADEXParameter(;a = 8.2,
    b = 7.3,
    cm = 484.2,
    v_rest = -80.0,
    tau_m = 32.5,
    tau_w = 296.8,
    v_thresh = -40.3,
    delta_T = 7.5,
    v_spike = -36.6,
    v_reset = -80.5,
    spike_delta = 44.3)

E = SNN.AD(;N = 1, param=adparam)
E.I = [3.9]
SNN.monitor(E, [:v,:I])
SNN.sim!([E],[], dt=0.1ms, duration=2000ms)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type VesselVolume
	VesselVolume()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"Pi value",
				:Default=>3.141593,
				:Hidden=>true,
				:Symbol=>"\\pi"
			)),
			acceleration(Dict{Symbol,Any}(
				:Brief=>"Gravity Acceleration",
				:Default=>9.81,
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["vertical","horizontal"],
				:Default=>"vertical"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["elliptical","hemispherical","flat"],
				:Default=>"flat"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Vessel diameter",
				:Default=>2,
				:Symbol=>"D_{i}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Side length of the cylinder shell",
				:Default=>6,
				:Symbol=>"L_{vessel}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Elliptical Head Total Volume",
				:Hidden=>true,
				:Symbol=>"V_{head}^{elliptical}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Hemispherical Head Total Volume",
				:Hidden=>true,
				:Symbol=>"V_{head}^{hemispherical}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Cylinder Total Volume",
				:Hidden=>true,
				:Symbol=>"V_{cylinder}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Vessel radius",
				:Hidden=>true,
				:Symbol=>"R_{cylinder}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Vessel total volume",
				:Protected=>true,
				:Symbol=>"V_{total}"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Vessel volume content",
				:Protected=>true,
				:Symbol=>"V_{filled}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"liquid height",
				:Protected=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Vessel cylinder shell Cross section area",
				:Hidden=>true,
				:Symbol=>"A_{cross}"
			)),
			[
				:(Across = 0.25*(Pi*Diameter^2)),
				:(Vtotal = Vhead_elliptical + Vcylinder),
				:(Vfilled = 0.25*Pi*(((Diameter*Level)/(0.25*Diameter))^2)*(0.25*Diameter-Level/3)),
				:(Vfilled = 0.25*Pi*(Diameter^2)*(Level - 0.25*Diameter/3)),
				:(Vtotal = Vhead_hemispherical + Vcylinder),
				:(Vfilled = 0.25*Pi*(Level^2)*(2*Diameter-4*Level/3)),
				:(Vfilled = 0.25*Pi*((2/3)*((0.5*Diameter)^3) - (0.25*(Diameter)^3) + Level*Diameter^2)),
				:(Vtotal = Vcylinder),
				:(Vfilled = Across*Level),
				:(Across = (radius^2)*acos((radius-Level)/radius)/"rad"-(radius-Level)*sqrt((2*radius*Level-Level^2))),
				:(Vtotal = Vhead_elliptical + Vcylinder),
				:(Vfilled = 0.5236*Level^2*(1.5*Diameter-Level) + Across*Lenght),
				:(Vtotal = Vhead_hemispherical + Vcylinder),
				:(Vfilled = 1.0472*Level^2*(1.5*Diameter-Level) + Across*Lenght),
				:(Vtotal = Vcylinder),
				:(Vfilled = Across*Lenght),
			],
			[
				"Vessel Cross Section Area","Vessel Total Volume","Vessel Filled Volume","Vessel Filled Volume","Vessel Total Volume","Vessel Filled Volume","Vessel Filled Volume","Vessel Total Volume","Vessel Filled Volume","Vessel Cross Section Area","Vessel Total Volume","Vessel Filled Volume","Vessel Total Volume","Vessel Filled Volume","Vessel Total Volume","Vessel Filled Volume",
			],
			[:Pi,:Gconst,:Orientation,:Heads,:Diameter,:Lenght,:Vhead_elliptical,:Vhead_hemispherical,:Vcylinder,:radius,],
			[:Vtotal,:Vfilled,:Level,:Across,]
		)
	end
	Pi::positive
	Gconst::acceleration
	Orientation::DanaSwitcher
	Heads::DanaSwitcher
	Diameter::length
	Lenght::length
	Vhead_elliptical::volume
	Vhead_hemispherical::volume
	Vcylinder::volume
	radius::length
	Vtotal::volume
	Vfilled::volume
	Level::length
	Across::area
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export VesselVolume
function set(in::VesselVolume)
	Gconst = 9.81 * "m/(s^2)"
	 Vhead_elliptical = (Pi*Diameter^3)/12
	 Vhead_hemispherical = (Pi*Diameter^3)/6
	 Vcylinder = 0.25*(Pi*Diameter^2)*Lenght
	 radius = 0.5*Diameter

end
function setEquationFlow(in::VesselVolume)
	let switch=Orientation
		if switch=="vertical"
			addEquation(1)
			let switch=Heads
				if switch=="elliptical"
					addEquation(2)
					if Level < 0.25*Diameter
						addEquation(3)
					else
						addEquation(4)
					end
				elseif switch=="hemispherical"
					addEquation(5)
					if Level < 0.5*Diameter
						addEquation(6)
					else
						addEquation(7)
					end
				elseif switch=="flat"
					addEquation(8)
					addEquation(9)
				end
			end
		elseif switch=="horizontal"
			addEquation(10)
			let switch=Heads
				if switch=="elliptical"
					addEquation(11)
					addEquation(12)
				elseif switch=="hemispherical"
					addEquation(13)
					addEquation(14)
				elseif switch=="flat"
					addEquation(15)
					addEquation(16)
				end
			end
		end
	end
end
function atributes(in::VesselVolume,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Model to calculate vessel total volume and vessel filled volume from
different geometries and orientations."
	fields[:Info]="== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;
"
	drive!(fields,_)
	return fields
end
VesselVolume(_::Dict{Symbol,Any})=begin
	newModel=VesselVolume()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module Norms

# A norm

using Base.Test

import Base.+, Base.norm

export Norm
export norm, +
export avg, norm2

immutable Norm
    sum::Float64
    sum2::Float64
    count::Float64
    min::Float64
    max::Float64
    maxabs::Float64
end
Norm() = Norm(0.0, 0.0, 0.0, Inf, -Inf, 0.0)

norm(val::Float64) = Norm(val, val*val, 1.0, val, val, abs(val))
+(x::Norm, y::Norm) = Norm(x.sum + y.sum, x.sum2 + y.sum2, x.count + y.count,
                           min(x.min, y.min), max(x.max, y.max),
                           max(x.maxabs, y.maxabs))

@test norm(+2.0) + norm(+3.0) == Norm( 5.0, 13.0, 2.0,  2.0,  3.0, 3.0)
@test norm(+2.0) + norm(-3.0) == Norm(-1.0, 13.0, 2.0, -3.0,  2.0, 3.0)
@test norm(-2.0) + norm(+3.0) == Norm( 1.0, 13.0, 2.0, -2.0,  3.0, 3.0)
@test norm(-2.0) + norm(-3.0) == Norm(-5.0, 13.0, 2.0, -3.0, -2.0, 3.0)

avg(x::Norm) = x.sum / x.count
norm2(x::Norm) = sqrt(x.sum2 / x.count)

@test avg(norm(2.0)) == 2.0
@test norm2(norm(2.0)) == 2.0

end


text:
using Expokit
using Test
using LinearAlgebra
using SparseArrays

# reference matrix exponential
expm_higham(A) = exp(A)

function ev1(n)

    v = zeros(n)
    v[1] = 1
    return v
end

function test_expmv(n::Int)

    A = sprand(n,n,0.4)
    v = ev1(n)

    t1 = @elapsed w1 = expmv(1.0, A, v)

    t2 = @elapsed w2 = expm_higham(Matrix(A))*v

    return norm(w1-w2)/norm(w2), t1, t2
end

function test_expmv2(n::Int)

    A = sprand(n,n,0.2) + 1im*sprand(n,n,0.2)
    v = ev1(n) + 0im*ev1(n)

    t1 = @elapsed w1 = expmv(1.0, A, v)

    t2 = @elapsed w2 = expm_higham(Matrix(A))*v

    return norm(w1-w2)/norm(w2), t1, t2
end

function test_expmv3()
    e1 = norm(expm_higham(pi/4*[0 1; 1 0])*[1.; 0.] - expmv(pi/4,[0 1; 1 0]|>sparse, [1.; 0.]))
    e2 = norm(expm_higham(-pi/4*1im*[0 1; 1 0])*[1.; 0.] - expmv(pi/4,-1im*[0 1; 1 0]|>sparse, [1.0+0im; 0.]))

    return e1, e2
end

# Testing custom type for A
struct LinearOp
    m
end

LinearAlgebra.mul!(y, lo::LinearOp, x) = mul!(y, lo.m, x)
Base.size(lo::LinearOp, i::Int) = size(lo.m, i)
Base.eltype(lo::LinearOp) = eltype(lo.m)

function test_expmv_linop(n::Int)
    A = LinearOp(sprand(n,n,0.2) + 1im*sprand(n,n,0.2))
    v = ev1(n) + 0im*ev1(n)

    t1 = @elapsed w1 = expmv(1.0, A, v, anorm=norm(A.m, Inf))

    t2 = @elapsed w2 = expmv(1.0, A, v)

    t3 = @elapsed w3 = expm_higham(Matrix(A.m))*v

    return max(norm(w1-w2)/norm(w2), norm(w1-w3)/norm(w3)), t1, t2, t3
end

println("testing real n=100 (first expmv, then expm)")
res, t1, t2 = test_expmv(100)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing complex n=100 (first expmv, then expm)")
res, t1, t2 = test_expmv2(100)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing real n=1000 (first expmv, then expm)")
res, t1, t2 = test_expmv(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing complex n=1000 (first expmv, then expm)")
res, t1, t2 = test_expmv2(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing linear operator n=1000 (first expmv, then expm)")
res, t1, t2, t3 = test_expmv_linop(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing 2x2 cases")
e1, e2 = test_expmv3()
println("residua: $e1, $e2\n")
@test e1 < 1e-10
@test e2 < 1e-10

function test_padm(n::Int)

    A = sprand(n,n,0.001)

    t1 = @elapsed w1 = padm(A)

    t2 = @elapsed w2 = expm_higham(Matrix(A))

    return norm(w1-w2)/norm(w2), t1, t2
end

function test_padm2(n::Int)

    A = sprand(n,n,0.0005) + 1im*sprand(n,n,0.0005)

    t1 = @elapsed w1 = padm(A)

    t2 = @elapsed w2 = expm_higham(Matrix(A))

    return norm(w1-w2)/norm(w2), t1, t2
end

println("testing real n=1000 (first padm, then expm)")
res, t1, t2 = test_padm(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing complex n=100 (first padm, then expm)")
res, t1, t2 = test_padm2(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

function test_chbv(n::Int)

    p = 0.1
    D = Matrix(Diagonal(-rand(n)))
    T = sprandn(n, n, p)
    H = T * D * transpose(T)  # random negative semidefinite symmetric matrix

    vec = randn(n)
    w1 = chbv(H, vec)
    t1 = @elapsed w1 = chbv(H, vec)

    w2 = expm_higham(Matrix(H)) * vec
    t2 = @elapsed  w2 = expm_higham(Matrix(H)) * vec

    return norm(w1-w2)/norm(w2), t1, t2
end

function test_chbv2(n::Int)

    p = 0.1
    D = Matrix(Diagonal(-rand(n)))
    T = sprandn(n, n, p) + sprandn(n, n, p)*im
    H = T * D * adjoint(T)  # random negative semidefinite hermitian matrix

    vec = randn(n) + randn(n)*im
    w1 = chbv(H, vec)
    t1 = @elapsed w1 = chbv(H, vec)

    w2 = expm_higham(Matrix(H)) * vec
    t2 = @elapsed w2 = expm_higham(Matrix(H)) * vec

    return norm(w1-w2)/norm(w2), t1, t2
end

println("testing real n=100 (first chbv, then expm)")
res, t1, t2 = test_chbv(100)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing complex n=100 (first chbv, then expm)")
res, t1, t2 = test_chbv2(100)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing real n=1000 (first chbv, then expm)")
res, t1, t2 = test_chbv(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

#= commented since this test takes sometime difs very long
println("testing complex n=1000 (first chbv, then expm)")
res, t1, t2 = test_chbv2(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6
=#

function test_phimv(n::Int)

    p = 0.1
    found = false
    A = sprand(n, n, p)
    u = rand(n)
    x = similar(u)
    while !found
        try
            copyto!(x, A\u)
            found = true
        catch
            A = sprand(n, n, p)
        end
    end
    vec = ev1(n)

    w1 = phimv(1.0, A, u, vec) # warmup
    t1 = @elapsed w1 = phimv(1.0, A, u, vec)

    w2 = expm_higham(Matrix(A))*(vec+x)-x # warmup
    t2 = @elapsed w2 = expm_higham(Matrix(A))*(vec+x)-x

    return norm(w1-w2)/norm(w2), t1, t2
end

function test_phimv2(n::Int)

    p = 0.1
    found = false
    A = sprand(n, n, p) + 1im*sprand(n, n, p)
    u = rand(n) + 1im*rand(n)
    x = similar(u)
    while !found
        try
            copyto!(x, A\u)
            found = true
        catch
            A = sprand(n, n, p) + 1im*sprand(n, n, p)
        end
    end
    vec = ev1(n) + 0im*ev1(n)

    w1 = phimv(1.0, A, u, vec) # warmup
    t1 = @elapsed w1 = phimv(1.0, A, u, vec)

    w2 = expm_higham(Matrix(A))*(vec+x)-x # warmup
    t2 = @elapsed w2 = expm_higham(Matrix(A))*(vec+x)-x

    return norm(w1-w2)/norm(w2), t1, t2
end

println("testing real n=100 (first phimv, then expm)")
res, t1, t2 = test_phimv(100)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing complex n=100 (first phimv, then expm)")
res, t1, t2 = test_phimv2(100)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing real n=1000 (first phimv, then expm)")
res, t1, t2 = test_phimv(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

println("testing complex n=1000 (first phimv, then expm)")
res, t1, t2 = test_phimv2(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6

function test_phimv_linop(n::Int)
    p = 0.1
    found = false
    A = LinearOp(sprand(n,n,p))
    u = rand(n)
    x = similar(u)
    while !found
        try
            copyto!(x, A.m\u)
            found = true
        catch
            A = LinearOp(sprand(n,n,p))
        end
    end
    vec = ev1(n)

    w1 = phimv(1.0, A, u, vec) # warmup
    t1 = @elapsed w1 = phimv(1.0, A, u, vec)

    w2 = expm_higham(Matrix(A.m))*(vec+x)-x # warmup
    t2 = @elapsed w2 = expm_higham(Matrix(A.m))*(vec+x)-x

    return norm(w1-w2)/norm(w2), t1, t2
end

println("testing real n=1000 (first phimv, then expm), the linear operator version.")
res, t1, t2 = test_phimv_linop(1000)
println("residuum: $res")
println("time dif: $(t1-t2)\n")
@test res < 1e-6


text:
# compiler support for working with run-time libraries

link_library!(mod::LLVM.Module, lib::LLVM.Module) = link_library!(mod, [lib])
function link_library!(mod::LLVM.Module, libs::Vector{LLVM.Module})
    # linking is destructive, so copy the libraries
    libs = [LLVM.Module(lib) for lib in libs]

    for lib in libs
        link!(mod, lib)
    end
end

#
# GPU run-time library
#

# get the path to a directory where we can put cache files (machine-specific, ephemeral)
# NOTE: maybe we should use XDG_CACHE_PATH/%LOCALAPPDATA%, but other Julia cache files
#       are put in .julia anyway so let's just follow suit for now.
function cachedir(depot=DEPOT_PATH[1])
    # this mimicks Base.compilecache. we can't just call the function, or we might actually
    # _generate_ a cache file, e.g., when running with `--compiled-modules=no`.
    entrypath, entryfile = Base.cache_file_entry(Base.PkgId(GPUCompiler))
    abspath(depot, entrypath, entryfile)
end

## higher-level functionality to work with runtime functions

function LLVM.call!(builder, rt::Runtime.RuntimeMethodInstance, args=LLVM.Value[];
                    state::Type=Nothing)
    bb = position(builder)
    f = LLVM.parent(bb)
    mod = LLVM.parent(f)
    ctx = context(mod)

    # get or create a function prototype
    if haskey(functions(mod), rt.llvm_name)
        f = functions(mod)[rt.llvm_name]
        ft = eltype(llvmtype(f))
    else
        ft = convert(LLVM.FunctionType, rt; ctx, state)
        f = LLVM.Function(mod, rt.llvm_name, ft)
    end

    # we may be calling this function after kernel state lowering,
    # in which case we need to manually get and pass the state.
    args = Value[args...]
    if state !== Nothing
        T_state = convert(LLVMType, state; ctx)

        state_intr = kernel_state_intr(mod, T_state)
        state_val = call!(builder, state_intr, Value[], "state")
        pushfirst!(args, state_val)
    end

    # runtime functions are written in Julia, while we're calling from LLVM,
    # this often results in argument type mismatches. try to fix some here.
    for (i,arg) in enumerate(args)
        if llvmtype(arg) != parameters(ft)[i]
            if (llvmtype(arg) isa LLVM.PointerType) &&
               (parameters(ft)[i] isa LLVM.IntegerType)
                # Julia pointers are passed as integers
                args[i] = ptrtoint!(builder, args[i], parameters(ft)[i])
            else
                error("Don't know how to convert ", arg, " argument to ", parameters(ft)[i])
            end
        end
    end

    call!(builder, f, args)
end

## functionality to build the runtime library

function emit_function!(mod, @nospecialize(job::CompilerJob), f, method)
    tt = Base.to_tuple_type(method.types)
    new_mod, meta = codegen(:llvm, similar(job, FunctionSpec(f, tt, #=kernel=# false));
                            optimize=false, libraries=false)
    ft = eltype(llvmtype(meta.entry))
    expected_ft = convert(LLVM.FunctionType, method; ctx=context(new_mod))
    if return_type(ft) != return_type(expected_ft)
        error("Invalid return type for runtime function '$(method.name)': expected $(return_type(expected_ft)), got $(return_type(ft))")
    end

    # recent Julia versions include prototypes for all runtime functions, even if unused
    pm = ModulePassManager()
    strip_dead_prototypes!(pm)
    run!(pm, new_mod)
    dispose(pm)

    temp_name = LLVM.name(meta.entry)
    # FIXME: on 1.6, there's no single global LLVM context anymore,
    #        but there's no API yet to pass a context to codegen.
    # round-trip the module through serialization to get it in the proper context.
    buf = convert(MemoryBuffer, new_mod)
    new_mod = parse(LLVM.Module, buf; ctx=context(mod))
    @assert context(mod) == context(new_mod)
    link!(mod, new_mod)
    entry = functions(mod)[temp_name]

    # if a declaration already existed, replace it with the function to avoid aliasing
    # (and getting function names like gpu_signal_exception1)
    name = method.llvm_name
    if haskey(functions(mod), name)
        decl = functions(mod)[name]
        @assert llvmtype(decl) == llvmtype(entry)
        replace_uses!(decl, entry)
        unsafe_delete!(mod, decl)
    end
    LLVM.name!(entry, name)
end

function build_runtime(@nospecialize(job::CompilerJob); ctx)
    mod = LLVM.Module("GPUCompiler run-time library"; ctx)

    # the compiler job passed into here is identifies the job that requires the runtime.
    # derive a job that represents the runtime itself (notably with kernel=false).
    source = FunctionSpec(identity, Tuple{Nothing}, false, nothing, job.source.world_age)
    job = CompilerJob(job.target, source, job.params)

    for method in values(Runtime.methods)
        def = if isa(method.def, Symbol)
            isdefined(runtime_module(job), method.def) || continue
            getfield(runtime_module(job), method.def)
        else
            method.def
        end
        emit_function!(mod, job, def, method)
    end

    optimize!(job, mod)

    mod
end

const runtime_lock = ReentrantLock()

@locked function load_runtime(@nospecialize(job::CompilerJob); ctx)
    lock(runtime_lock) do
        # find the first existing cache directory (for when dealing with layered depots)
        cachedirs = [cachedir(depot) for depot in DEPOT_PATH]
        filter!(isdir, cachedirs)
        input_dir = if isempty(cachedirs)
            nothing
        else
            first(cachedirs)
        end

        # we are only guaranteed to be able to write in the current depot
        output_dir = cachedir()

        # if both aren't equal, copy pregenerated runtime libraries to our depot
        # NOTE: we don't just lazily read from the one and write to the other, because
        #       once we generate additional runtimes in the output dir we don't know if
        #       it's safe to load from other layers (since those could have been invalidated)
        if input_dir !== nothing && input_dir != output_dir
            mkpath(dirname(output_dir))
            cp(input_dir, output_dir)
        end

        slug = runtime_slug(job)
        name = "runtime_$(slug).bc"
        path = joinpath(output_dir, name)

        lib = try
            if ispath(path)
                open(path) do io
                    parse(LLVM.Module, read(io); ctx)
                end
            end
        catch ex
            @warn "Failed to load GPU runtime library at $path" exception=(ex, catch_backtrace())
            nothing
        end

        if lib === nothing
            @debug "Building the GPU runtime library at $path"
            mkpath(output_dir)
            lib = build_runtime(job; ctx)

            # atomic write to disk
            temp_path, io = mktemp(dirname(path); cleanup=false)
            write(io, lib)
            close(io)
            Base.rename(temp_path, path; force=true)
        end

        return lib
    end
end

# remove the existing cache
# NOTE: call this function from global scope, so any change triggers recompilation.
function reset_runtime()
    lock(runtime_lock) do
        rm(cachedir(); recursive=true, force=true)
        # create an empty cache directory. since we only ever load from the first existing cachedir,
        # this effectively invalidates preexisting caches in lower layers of the depot.
        mkpath(cachedir())
    end

    return
end


text:
using Base.Test
using Inv

a = [2,3,1,5]
#@test Inv.merge_sort(a) == [1,2,3,5]

let left = [1,4,6,8],right = [0,2,5]
    @test Inv.merge(left, right) == [0,1,2,4,5,6,8]
end

let data = [5,8,2,4,1]
    @test Inv.merge_sort(data) == [1,2,4,5,8]
end

let data = [2,1,0,6]
    @test Inv.inv_count(data) == 2+1
end


text:
################################################################################
### Measurements (constructors + kernels)
################################################################################

# Notes on math and stuff:
#=
- Use Wicks theorem to go from one expectation value with 2N operators to
  N expectation values with 2 operators each
- use G_{ij} = ⟨c_i c_j^†⟩
- use ⟨c_i^† c_j⟩ = δ_{ij} - ⟨c_j c_i^†⟩ = (I - G)_{ji}
- we define G[i,j] as spin up, G[i+N,j+N] as spin down (N = number of sites)
- we used (I - G)[i, j+N] = G[i, j+N] etc when possible
- Gl0_{i,j} = ⟨c_i(l) c_j(0)^†⟩
- ⟨c_i(0) c_j(l)^†⟩ = ⟨(c_j(l) c_i(0)^†)^†⟩ = (⟨c_j(l) c_i(0)^†⟩)^†
                    = Gl0_{j,i}^† = Gl0_{j,i}^*
=#
function checkflavors(model, N=2)
    if nflavors(model) != N
        @warn(
            "$N flavors are required, but $(nflavors(model)) have been found"
        )
    end
    nothing
end

# This has lattice_iteratorator = Nothing, because it straight up copies G
function greens_measurement(mc::DQMC, model::Model, greens_iterator=Greens; kwargs...)
    N = length(lattice(model)) * nflavors(model)
    Measurement(
        mc, model, greens_iterator, Nothing, greens_kernel, shape = (N, N); kwargs...
    )
end
greens_kernel(mc, model, G::AbstractArray) = G

function occupation(mc::DQMC, model::Model, kwargs...)
    Measurement(mc, model, Greens, EachSiteAndFlavor, occupation_kernel; kwargs...)
end
occupation_kernel(mc, model, i, G::AbstractArray) = 1 - G[i, i]

function charge_density(
        mc::DQMC, model::Model, greens_iterator;
        lattice_iterator = EachSitePairByDistance, kwargs...
    )
    checkflavors(model)
    Measurement(mc, model, greens_iterator, lattice_iterator, cdc_kernel; kwargs...)
end
charge_density_correlation(mc, m; kwargs...) = charge_density(mc, m, Greens; kwargs...)
function charge_density_susceptibility(mc, m; kwargs...)
    charge_density(mc, m, CombinedGreensIterator; kwargs...)
end

function cdc_kernel(mc, model, i, j, G::AbstractArray)
    N = length(lattice(mc))
    # ⟨n↑n↑⟩
    (1 - G[i, i])       * (1 - G[j, j]) +
    (I[j, i] - G[j, i]) * G[i, j] +
    # ⟨n↑n↓⟩
    (1 - G[i, i]) * (1 - G[j+N, j+N]) -
    G[j+N, i]     * G[i, j + N] +
    # ⟨n↓n↑⟩
    (1 - G[i+N, i+N]) * (1 - G[j, j]) -
    G[j, i+N]         * G[i+N, j] +
    # ⟨n↓n↓⟩
    (1 - G[i+N, i+N])           * (1 - G[j+N, j+N]) +
    (I[j, i] - G[j+N, i+N]) *  G[i+N, j+N]
end
function cdc_kernel(mc, model, i, j, packed_greens::NTuple{4})
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(mc))
    # ⟨n↑(l)n↑⟩
    (1 - Gll[i, i]) * (1 - G00[j, j]) -
    G0l[j, i] * Gl0[i, j] +
    # ⟨n↑n↓⟩
    (1 - Gll[i, i]) * (1 - G00[j+N, j+N]) -
    G0l[j+N, i] * Gl0[i, j+N] +
    # ⟨n↓n↑⟩
    (1 - Gll[i+N, i+N]) * (1 - G00[j, j]) -
    G0l[j, i+N] * Gl0[i+N, j] +
    # ⟨n↓n↓⟩
    (1 - Gll[i+N, i+N]) * (1 - G00[j+N, j+N]) -
    G0l[j+N, i+N] * Gl0[i+N, j+N]
end

"""
    magnetization_measurement(mc, model, dir[; lattice_iterator, kwargs...])

Returns the x, y or z magnetization measurement given `dir = :x`, `:y` or `:z`
respectively.

NOTE:
We're skipping the multiplication by `-1im` during the measurement of the y
magnetization. To get the correct result, multiply the final result by `-1im`.
"""
function magnetization(
        mc::DQMC, model::Model, dir::Symbol;
        lattice_iterator = EachSite, kwargs...
    )
    checkflavors(model)
    if     dir == :x; kernel = mx_kernel
    elseif dir == :y; kernel = my_kernel
    elseif dir == :z; kernel = mz_kernel
    else throw(ArgumentError("`dir` must be :x, :y or :z, but is $dir"))
    end
    Measurement(mc, model, Greens, lattice_iterator, kernel; kwargs...)
end
function mx_kernel(mc, model, i, G::AbstractArray)
    N = length(lattice(model))
    -G[i+N, i] - G[i, i+N]
end
function my_kernel(mc, model, i, G::AbstractArray)
    N = length(lattice(model))
    G[i+N, i] - G[i, i+N]
end
function mz_kernel(mc, model, i, G::AbstractArray)
    N = length(lattice(model))
    G[i+N, i+N] - G[i, i]
end

function spin_density(
        dqmc, model, dir::Symbol, greens_iterator;
        lattice_iterator = EachSitePairByDistance, kwargs...
    )
    checkflavors(model)
    if     dir == :x; kernel = sdc_x_kernel
    elseif dir == :y; kernel = sdc_y_kernel
    elseif dir == :z; kernel = sdc_z_kernel
    else throw(ArgumentError("`dir` must be :x, :y or :z, but is $dir"))
    end
    Measurement(dqmc, model, greens_iterator, lattice_iterator, kernel; kwargs...)
end
spin_density_correlation(args...; kwargs...) = spin_density(args..., Greens; kwargs...)
function spin_density_susceptibility(args...; kwargs...)
    spin_density(args..., CombinedGreensIterator; kwargs...)
end

function sdc_x_kernel(mc, model, i, j, G::AbstractArray)
    N = length(lattice(model))
    G[i+N, i] * G[j+N, j] - G[j+N, i] * G[i+N, j] +
    G[i+N, i] * G[j, j+N] + (I[j, i] - G[j, i]) * G[i+N, j+N] +
    G[i, i+N] * G[j+N, j] + (I[j, i] - G[j+N, i+N]) * G[i, j] +
    G[i, i+N] * G[j, j+N] - G[j, i+N] * G[i, j+N]
end
function sdc_x_kernel(mc, model, i, j, packed_greens::NTuple{4})
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    Gll[i+N, i] * G00[j+N, j] - G0l[j+N, i] * Gl0[i+N, j] +
    Gll[i+N, i] * G00[j, j+N] - G0l[j, i] * Gl0[i+N, j+N] +
    Gll[i, i+N] * G00[j+N, j] - G0l[j+N, i+N] * Gl0[i, j] +
    Gll[i, i+N] * G00[j, j+N] - G0l[j, i+N] * Gl0[i, j+N]
end

function sdc_y_kernel(mc, model, i, j, G::AbstractArray)
    N = length(lattice(model))
    - G[i+N, i] * G[j+N, j] + G[j+N, i] * G[i+N, j] +
      G[i+N, i] * G[j, j+N] + (I[j, i] - G[j, i]) * G[i+N, j+N] +
      G[i, i+N] * G[j+N, j] + (I[j, i] - G[j+N, i+N]) * G[i, j] -
      G[i, i+N] * G[j, j+N] + G[j, i+N] * G[i, j+N]
end
function sdc_y_kernel(mc, model, i, j, packed_greens::NTuple{4})
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    - Gll[i+N, i] * G00[j+N, j] + G0l[j+N, i] * Gl0[i+N, j] +
      Gll[i+N, i] * G00[j, j+N] - G0l[j, i] * Gl0[i+N, j+N] +
      Gll[i, i+N] * G00[j+N, j] - G0l[j+N, i+N] * Gl0[i, j] -
      Gll[i, i+N] * G00[j, j+N] + G0l[j, i+N] * Gl0[i, j+N]
end

function sdc_z_kernel(mc, model, i, j, G::AbstractArray)
    N = length(lattice(model))
    (1 - G[i, i]) * (1 - G[j, j])         + (I[j, i] - G[j, i]) * G[i, j] -
    (1 - G[i, i]) * (1 - G[j+N, j+N])     + G[j+N, i] * G[i, j+N] -
    (1 - G[i+N, i+N]) * (1 - G[j, j])     + G[j, i+N] * G[i+N, j] +
    (1 - G[i+N, i+N]) * (1 - G[j+N, j+N]) + (I[j, i] - G[j+N, i+N]) * G[i+N, j+N]
end
function sdc_z_kernel(mc, model, i, j, packed_greens::NTuple{4})
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    (1 - Gll[i, i])     * (1 - G00[j, j])     - G0l[j, i] * Gl0[i, j] -
    (1 - Gll[i, i])     * (1 - G00[j+N, j+N]) + G0l[j+N, i] * Gl0[i, j+N] -
    (1 - Gll[i+N, i+N]) * (1 - G00[j, j])     + G0l[j, i+N] * Gl0[i+N, j] +
    (1 - Gll[i+N, i+N]) * (1 - G00[j+N, j+N]) - G0l[j+N, i+N] * Gl0[i+N, j+N]
end

function pairing(
        dqmc::DQMC, model::Model, greens_iterator;
        K = 1+length(neighbors(lattice(model), 1)),
        lattice_iterator = EachLocalQuadByDistance{K}, kwargs...
    )
    Measurement(dqmc, model, greens_iterator, lattice_iterator, pc_kernel; kwargs...)
end
pairing_correlation(mc, m; kwargs...) = pairing(mc, m, Greens; kwargs...)
pairing_susceptibility(mc, m; kwargs...) = pairing(mc, m, CombinedGreensIterator; kwargs...)
function pc_kernel(mc, model, src1, trg1, src2, trg2, G::AbstractArray)
    # verified against ED for each (src1, src2, trg1, trg2)
    # Δ_v(src1, trg1) Δ_v^†(src2, trg2)
    # G_{i, j}^{↑, ↑} G_{i+d, j+d}^{↓, ↓} - G_{i, j+d}^{↑, ↓} G_{i+d, j}^{↓, ↑}
    N = length(lattice(model))
    G[src1, src2] * G[trg1+N, trg2+N] - G[src1, trg2+N] * G[trg1+N, src2]
end
function pc_kernel(mc, model, src1, trg1, src2, trg2, packed_greens::NTuple{4})
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    Gl0[src1, src2] * Gl0[trg1+N, trg2+N] - Gl0[src1, trg2+N] * Gl0[trg1+N, src2]
end

#=
Λ^L = Λxx(qx->0, qy=0, iω=0)
    = ∑_r ∫[0, ß] dτ ⟨j_x(r, τ) j_x(0, 0)⟩ exp(iqr) exp(iωτ)
j_x(r, τ) = it ∑_σ [c_{r+e_x, σ}^†(τ) c_{r, σ}(τ) - c_{r, σ}^†(τ) c_{r+e_x, σ}(τ)]

# reframe as:
0 -> r, r -> r + Δr, e_x -> r+e_x, r+e_x -> r+e_x+Δr
then use EachLocalQuadByDistance
Well actually no, because we want e_x to only pick one neighbor.

Λ_{trg-src}(Δr) = \sum_src1 cc_kernel(...)
∑_Δr exp(iq⋅Δr) real(Λ_{trg-src}(Δr)) # Why real? this is silly...
where q = Δq_i = k_1, k_2 (or k_1+k_2?)

                       v- (0, qy)
ρs = 1/8 * (real(Λxxq[1,2]) - real(Λxxq[2,1]) + real(Λyyq[1,2]) - real(Λyyq[2,1]))
                                (qx, 0) -^
This should correspond to
ρ = -K_s - Λ^T = Λxx^L - Λxx^T = Λxx((0, qy)) - Λxx((qx, 0))
  ~ real(Λxxq[1,2]) + real(Λyyq[1,2]) - real(Λxxq[2,1]) - real(Λyyq[2,1])
    ^-        sum over NNs?        -^

So I need to compute reciprocal lattice vectors I guess
So what is qx, qy then? I assume those reciprocal lattice vectors? Or literally
just a dx/dy?

Any real-space runtime computation could be EachLocalQuadByDistance, we can do
summing, direction picking and Fourier afterwards...

To be more efficient:
Should implement new lattice iterator that picks trg1, trg2 to be in the same
direction and skips trg == src
EachSyncedNNQuadByDistance{K}?
=#
function current_current_susceptibility(
        dqmc::DQMC, model::Model;
        K = 1+length(neighbors(lattice(model), 1)),
        greens_iterator = CombinedGreensIterator,
        lattice_iterator = EachLocalQuadBySyncedDistance{K}, kwargs...
    )
    Measurement(dqmc, model, greens_iterator, lattice_iterator, cc_kernel; kwargs...)
end
# current_current_correlation(mc, m; kwargs...) = current_current(mc, m, Greens; kwargs...)
# current_current_susceptibility(mc, m; kwargs...) = current_current(mc, m, CombinedGreensIterator; kwargs...)

function cc_kernel(mc, model, src1, trg1, src2, trg2, packed_greens::NTuple{4})
    # This should compute
    # ⟨j_{trg1-src1}(src1, τ) j_{trg2-src2}(src2, 0)⟩
    # where (trg-src) picks a direction (e.g. NN directions)
    # and (src1-src2) is the distance vector that the Fourier transform applies to
    # From dos Santos: Introduction to Quantum Monte Carlo Simulations
    # j_{trg-src}(src, τ) = it \sum\sigma (c^\dagger(trg,\sigma, \tau) c(src, \sigma, \tau) - c^\dagger(src, \sigma, \tau) c(trg, \sigma \tau))
    # where i -> src, i+x -> trg as a generalization
    # and t is assumed to be hopping matrix element, generalizing to
    # = i \sum\sigma (T[trg, src] c^\dagger(trg,\sigma, \tau) c(src, \sigma, \tau) - T[src, trg] c^\dagger(src, \sigma, \tau) c(trg, \sigma \tau))

	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    T = mc.s.hopping_matrix
    output = zero(eltype(G00))

    # Iterate through (spin up, spin down)
    for σ1 in (0, N), σ2 in (0, N)
        s1 = src1 + σ1; t1 = trg1 + σ1
        s2 = src2 + σ2; t2 = trg2 + σ2
        output +=
            (T[s1, t1] * Gll[t1, s1] - T[t1, s1] * Gll[s1, t1]) *
            (T[s2, t2] * G00[t2, s2] - T[t2, s2] * G00[s2, t2]) +
            T[t1, s1] * T[t2, s2] * (- G0l[s2, t1]) * Gl0[s1, t2] -
            T[s1, t1] * T[t2, s2] * (- G0l[s2, s1]) * Gl0[t1, t2] -
            T[t1, s1] * T[s2, t2] * (- G0l[t2, t1]) * Gl0[s1, s2] +
            T[s1, t1] * T[s2, t2] * (- G0l[t2, s1]) * Gl0[t1, s2]
            # Why no I?
            # T[t1, s1] * T[t2, s2] * (I[s2, t1] - G0l[s2, t1]) * Gl0[s1, t2] -
            # T[s1, t1] * T[t2, s2] * (I[s2, s1] - G0l[s2, s1]) * Gl0[t1, t2] -
            # T[t1, s1] * T[s2, t2] * (I[t2, t1] - G0l[t2, t1]) * Gl0[s1, s2] +
            # T[s1, t1] * T[s2, t2] * (I[t2, s1] - G0l[t2, s1]) * Gl0[t1, s2]

        # Uncompressed Wicks expansion
        # output += T[t1, s1] * T[t2, s2] *
        #     ((I[s1, t1] - Gll[s1, t1]) * (I[s2, t2] - G00[s2, t2]) +
        #     (I[s2, t1] - G0l[s2, t1]) * Gl0[s1, t2])
        # output -= T[s1, t1] * T[t2, s2] *
        #     ((I[t1, s1] - Gll[t1, s1]) * (I[s2, t2] - G00[s2, t2]) +
        #     (I[s2, s1] - G0l[s2, s1]) * Gl0[t1, t2])
        # output -= T[t1, s1] * T[s2, t2] *
        #     ((I[s1, t1] - Gll[s1, t1]) * (I[t2, s2] - G00[t2, s2]) +
        #     (I[t2, t1] - G0l[t2, t1]) * Gl0[s1, s2])
        # output += T[s1, t1] * T[s2, t2] *
        #     ((I[t1, s1] - Gll[t1, s1]) * (I[t2, s2] - G00[t2, s2]) +
        #     (I[t2, s1] - G0l[t2, s1]) * Gl0[t1, s2])
    end

    output
end

function boson_energy_measurement(dqmc, model; kwargs...)
    Measurement(dqmc, model, Nothing, Nothing, energy_boson; kwargs...)
end


text:
# Utils.jl

# [[file:~/Documents/Julia/scrap.org::*Utils.jl][Utils.jl:1]]
walk(x, inner, outer) = outer(x)
walk(x::T, inner, outer) where {T<:AbstractSymExpr} = outer(T(inner(x.op), map(inner, x.args)))
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))
postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

removeiden(ex::AbstractSymExpr) = (ex.op == Sym(:identity)) && (length(ex.args) == 1) ? ex.args[1] : ex
removeiden(x) = x

stripiden(x) = x
stripiden(x::AbstractSymExpr) = postwalk(removeiden, x)

replace_sym(x, p::Pair{AbstractSym, Number}) = x
function replace_sym(ex::Symbolic, p::Pair)
    postwalk(sym -> sym == p.first ? p.second : sym, ex)
end

(ex::SymExpr)(p::Pair) = replace_sym(ex, p)
(ex::Sym)(p::Pair) = replace_sym(ex, p)

Base.promote(::Type{Sym}) = SymExpr
Base.promote(::Type{SymExpr}) = SymExpr

Base.promote(x::T, y::Number) where {T<:Symbolic} = (promote(T)(:identity, [x]), promote(T)(:identity, [y]))
Base.promote(x::Number, y::T) where {T<:Symbolic} = (promote(T)(:identity, [x]), promote(T)(:identity, [y]))
Base.promote(x::Sym, y::SymExpr) = (SymExpr(:identity, [x]), y)
Base.promote(x::SymExpr, y::Sym) = (x, SymExpr(:identity, [y]))

promote_SymForm(x::Number, y::Union{Sym,SymExpr}) = SymExpr
promote_SymForm(x::Union{Sym,SymExpr}, y::Number) = SymExpr
promote_SymForm(x::Union{Sym,SymExpr}, y::Union{Sym,SymExpr}) = SymExpr
# Utils.jl:1 ends here