text:
### Struct variable
struct OptimalPowerFlowDC
    main::Array{Float64,2}
    flow::Array{Float64,2}
    generation::Array{Float64,2}
end

### DC optimal power flow, created by Ognjen Kundacina
function rundcopf(system, num, settings, info)
    printstyled("Algorithm: DC optimal power flow\n"; bold = true)

    ########## Pre-processing ##########
    busi, type, Pload, Gshunt, Tini, Pmax, Pmin, reactance, transTap, transShift, branchOn, longTermRating, angleMin, angleMax,
    costModel, numOfCoeffs, coeffs = read_opfdcsystem(system, num)

    main = fill(0.0, num.Nbus, 6)
    flow = fill(0.0, num.Nbranch, 5)
    gene = fill(0.0, num.Ngen, 2)
    Va, Pinj, Pbus, Pij, Pji, Pgen = write_dcsystem(system, num, main, flow, gene)

  algtime = @elapsed begin
    ########## Convert in integers ##########
    geni = convert(Array{Int64,1}, system.generator[:, 1])
    genOn = convert(Array{Int64,1}, system.generator[:, 8])
    fromi = convert(Array{Int64,1}, system.branch[:, 2])
    toi = convert(Array{Int64,1}, system.branch[:, 3])

    ########## Numbering and slack bus ##########
    numbering = false
    bus = collect(1:num.Nbus)
    slack = 0
    for i = 1:num.Nbus
        if bus[i] != busi[i]
            numbering = true
        end
        if type[i] == 3
            slack = bus[i]
        end
    end
    if numbering
        println("The new bus numbering is running.")
    end
    if slack == 0
        slack = 1
        println("The slack bus is not found. Slack bus is the first bus.")
    end
    gen = renumber(geni, busi, bus, numbering)
    from = renumber(fromi, busi, bus, numbering)
    to = renumber(toi, busi, bus, numbering)

    ########## Optimization model ##########
    model = Model(Ipopt.Optimizer)
    @variable(model, thetas[bus])
    @constraint(model, thetas[slack] == 0.0)
    @variable(model, Pmin[i] / system.basePower <= pGen[i = 1:num.Ngen] <= Pmax[i] / system.basePower)
    @inbounds for i = 1:num.Ngen
        if genOn[i] == 0
            @constraint(model, pGen[i] == 0.0)
        else
            set_start_value(pGen[i], (Pmax[i] + Pmin[i]) / (2 * system.basePower))
        end
    end

    ########## Ybus matrix ##########
    Ybus, admitance, Pshift = ybusdc(system, num, bus, from, to, branchOn, transTap, reactance, transShift)

    ########## Branch flow constraints ##########
    @inbounds for i = 1:num.Nbranch
        if longTermRating[i] ≉  0 && longTermRating[i] < 10^10
            upf = longTermRating[i] / system.basePower + admitance[i] * transShift[i] * pi / 180
            upt = longTermRating[i] / system.basePower - admitance[i] * transShift[i] * pi / 180
            @constraint(model, -upt <= admitance[i] * thetas[from[i]] - admitance[i] * thetas[to[i]] <= upf)
        end
    end

    ########## Branch angle constraints ##########
    @inbounds for i = 1:num.Nbranch
        if angleMin[i] > -360 && angleMax[i] < 360
            angMinRad = angleMin[i] * (pi / 180)
            angMaxRad = angleMax[i] * (pi / 180)
            @constraint(model, angMinRad <= thetas[from[i]] - thetas[to[i]] <= angMaxRad)
        end
    end

    ########## Power mismatch constraints ##########
    Amis = [Ybus sparse(gen, collect(1:num.Ngen), -1, num.Nbus, num.Ngen)]
    bmis = - (Pload + Gshunt) / system.basePower - Pshift;
    @constraint(model, Amis * [thetas; pGen] .== bmis)

    Qpg, cpg, kpg, nPWL, pwlIds = process_cost_table(system, num, costModel, coeffs, numOfCoeffs)
    if nPWL == 0
        @objective(model, Min, sum(0.5 * Qpg[i] * pGen[i]^2 + cpg[i] * pGen[i] + kpg[i] for i = 1:num.Ngen))
    elseif nPWL > 0
        @variable(model, pwlVars[i = 1:nPWL])
        pgbas = 1
        qgbas = zeros(Int64, 0)
        nq = 0
        ybas = 1 + num.Ngen + nq
        Apwl, bpwl = makeApwl(system, num, nPWL, pwlIds, pgbas, qgbas, ybas, numOfCoeffs)
        for i = 1:(length(bpwl))
            @constraint(model, dot(Apwl[i,:], [pGen; pwlVars]) <= bpwl[i])
        end
        @objective(model, Min, sum(0.5 * Qpg[i] * pGen[i]^2 + cpg[i] * pGen[i] + kpg[i] for i = 1:num.Ngen) + sum(pwlVars[i] for i = 1:nPWL))
    end

    optimize!(model)
    Ti = value.(thetas) .+ (pi / 180) * Tini[slack]

    ########## Post-processing ##########
    for i = 1:num.Nbranch
        if branchOn[i] == 1
            Pij[i] = admitance[i] * (Ti[from[i]] - Ti[to[i]] - (pi / 180) * transShift[i]) * system.basePower
            Pji[i] = -Pij[i]
        end
    end

    Pgen[:] = value.(pGen) .* system.basePower
    for (k, i) in enumerate(gen)
        if genOn[k] == 1
            Pbus[i] += Pgen[k]
        end
    end

    for i = 1:num.Nbus
        Va[i] = (180 / pi) * Ti[i]
        Pinj[i] = Pbus[i] - Pload[i]
        if i == slack
            I = 0.0
            for j in Ybus.colptr[i]:(Ybus.colptr[i + 1] - 1)
                row = Ybus.rowval[j]
                I += Ybus[slack, row] * Ti[row]
            end
            Pinj[slack] = (I + Pshift[slack]) * system.basePower + Gshunt[slack]
        end
    end

 end # algtime

    ########## Results ##########
    results = OptimalPowerFlowDC(main, flow, gene)
    header, group = results_flowdc(system, num, settings, results, slack, algtime)
    if !isempty(settings.save)
        savedata(results, system; info = info, group = group, header = header, path = settings.save)
    end

    return results
end

### Process data needed for the objective function
function process_cost_table(system, num, costModel, coeffs, numOfCoeffs)
    linTerm = fill(0.0, num.Ngen); quadTerm = fill(0.0, num.Ngen); constTerm = fill(0.0, num.Ngen)
    nPWL = 0; pwlIds = zeros(Int64, 0)

    @inbounds for i = 1:num.Ngen
        if costModel[i] == 1
            ########## PWL cost model ##########
            if numOfCoeffs[i] == 2
                ########## Single block pwl cost, convert to polynomial ##########
                x0 = system.generatorcost[i, 5]
                y0 = system.generatorcost[i, 6]
                x1 = system.generatorcost[i, 7]
                y1 = system.generatorcost[i, 8]
                m = (y1 - y0) / (x1 - x0)
                b = y0 - m * x0
                linTerm[i] = m
                constTerm[i] = b
            else
                nPWL = nPWL + 1
                append!(pwlIds, i)
            end
        elseif costModel[i] == 2
            ########## Polynomial ##########
            if numOfCoeffs[i] == 2
                linTerm[i] = coeffs[i, 1]
                constTerm[i] = coeffs[i, 2]
            elseif numOfCoeffs[i] == 3
                ########## Quadratic generator costs ##########
                quadTerm[i] = coeffs[i, 1]
                linTerm[i] = coeffs[i, 2]
                constTerm[i] = coeffs[i, 3]
            end
        end
    end
    Qpg = 2 * quadTerm * system.basePower^2
    cpg = linTerm * system.basePower
    kpg = constTerm

    return Qpg, cpg, kpg, nPWL, pwlIds
end

### Construct parameters for linear basin constraints
function makeApwl(system, num, nPWL, pwlIds, pgbas, qgbas, ybas, numOfCoeffs)
    if nPWL == 0
        Apwl = sparse([], [], [], 0, ybas + nPWL - 1, 0)
        bpwl = zeros(Int64, 0)

        return Apwl, bpwl
    end

    numCostPoints = sum(numOfCoeffs[i] for i in pwlIds)
    Apwl = zeros(Float64, numCostPoints - nPWL, ybas + nPWL - 1)
    bpwl = zeros(Int64, 0)

    k = 1
    @inbounds for i in pwlIds
        ns = numOfCoeffs[i]
        p = @view(system.generatorcost[i, 5:2:(5 + 2 * ns - 1)]) / system.basePower
        c = @view(system.generatorcost[i, 6:2:(5 + 2 * ns)])
        m = diff(c) ./ diff(p)
        if any(diff(p) == 0)
            println("makeApwl: Bad x axis data in the generator cost matrix, row: ", i)
        end
        b = m .* p[1:(ns - 1)] - c[1:(ns - 1)]
        bpwl = [bpwl; b]
        ########## Q and P costs ##########
        if i > num.Ngen
            sidx = qgbas + (i-ng) - 1
        else
            sidx = pgbas + i - 1
        end
        @inbounds for j = 1:(ns - 1)
            Apwl[j + k - 1, sidx] = m[j]
        end
        k = k + ns - 1
    end

    k = 1; j = 1
    @inbounds for i in pwlIds
        ns = numOfCoeffs[i]
        @inbounds for it = k:(k + ns - 2)
            Apwl[it, ybas + j - 1] = -1.0
        end
        k = k + ns - 1
        j = j + 1
    end

    return Apwl, bpwl
end

### Read data
function read_opfdcsystem(system, num)
    busi = @view(system.bus[:, 1])
    type = @view(system.bus[:, 2])
    Pload = @view(system.bus[:, 3])
    Gshunt = @view(system.bus[:, 5])
    Tini = @view(system.bus[:, 9])

    Pmax = @view(system.generator[:, 9])
    Pmin = @view(system.generator[:, 10])

    reactance = @view(system.branch[:, 5])
    transTap = @view(system.branch[:, 10])
    transShift = @view(system.branch[:, 11])
    branchOn = @view(system.branch[:, 12])
    longTermRating = @view(system.branch[:, 7])
    angleMin = @view(system.branch[:, 13])
    angleMax = @view(system.branch[:, 14])

    costModel = @view(system.generatorcost[:, 1])
    numOfCoeffs = convert(Array{Int64,1}, system.generatorcost[:, 4])
    coeffs = @view(system.generatorcost[:, 5:end])

    return busi, type, Pload, Gshunt, Tini,
            Pmax, Pmin,
            reactance, transTap, transShift, branchOn, longTermRating, angleMin, angleMax,
            costModel, numOfCoeffs, coeffs
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
# Generating Set Search as described in Kolda2003:
#  Kolda, Tamara G., Robert Michael Lewis, and Virginia Torczon. "Optimization
#  by direct search: New perspectives on some classical and modern methods."
#  SIAM review 45.3 (2003): 385-482.
#

"""
A supertype for all generating set searcher-like algorithms.
"""
abstract type DirectSearcher <: SteppingOptimizer end

"""
`DirectionGenerator` generates the search directions to use at each step of
a GSS search.
"""
abstract type DirectionGenerator end

struct ConstantDirectionGen <: DirectionGenerator
    directions::Matrix{Float64}

    ConstantDirectionGen(directions) = new(directions)
end

directions_for_k(cg::ConstantDirectionGen, k) =
    cg.directions # For a ConstantDirectionGen it is always the same regardless of k...

# We can easily do a compass search with GSS by generating directions
# individually (+ and -) for each coordinate.
compass_search_directions(n) = ConstantDirectionGen([Matrix{Float64}(I, n,n) -Matrix{Float64}(I, n, n)])

const GSSDefaultParameters = ParamsDict(
    :DeltaTolerance => 1e-10,       # GSS has converged if the StepSize drops below this tolerance level
    :InitialStepSizeFactor => 0.50,        # Factor times the minimum search space diameter to give the initial StepSize
    :RandomDirectionOrder => true,  # Randomly shuffle the order in which the directions are used for each step
    :StepSizeGamma => 2.0,          # Factor by which step size is multiplied if improved point is found. Should be >= 1.0.
    :StepSizePhi => 0.5,            # Factor by which step size is multiplied if NO improved point is found. Should be < 1.0.
    :StepSizeMax => prevfloat(typemax(Float64)) # A limit on the step size can be set but is typically not => Inf.
)

calc_initial_step_size(ss, stepSizeFactor = 0.5) = stepSizeFactor * minimum(dimdelta(ss))

"""
Generating Set Search as described in Kolda2003:
  Kolda, Tamara G., Robert Michael Lewis, and Virginia Torczon. "Optimization
  by direct search: New perspectives on some classical and modern methods."
  SIAM review 45.3 (2003): 385-482.
"""
mutable struct GeneratingSetSearcher{V<:Evaluator, D<:DirectionGenerator, E<:EmbeddingOperator} <: DirectSearcher
    direction_gen::D
    evaluator::V
    embed::E
    search_space::SearchSpace
    n::Int
    k::Int
    random_dir_order::Bool       # shuffle the order of directions?
    step_size_factor::Float64    # initial step size factor
    step_size_gamma::Float64     # step size factor if improved
    step_size_phi::Float64       # step size factor if no improvement
    step_size_max::Float64       # maximal step size
    step_tol::Float64            # step delta tolerance
    step_size::Float64           # current step size
    x::AbstractIndividual
    xfitness::Float64

    function GeneratingSetSearcher(
            evaluator::V, dgen::D, embed::E,
            random_dir_order::Bool, step_size_factor::Float64, step_size_gamma::Float64,
            step_size_phi::Float64, step_size_max::Float64,
            step_tol::Float64) where {V<:Evaluator, D<:DirectionGenerator, E<:EmbeddingOperator}
        # FIXME check parameters ranges
        n = numdims(evaluator)
        ss = search_space(evaluator)
        x = rand_individual(ss)
        new{V, D, E}(dgen, evaluator, embed, ss, n, 0,
            random_dir_order, step_size_factor,
            step_size_gamma, step_size_phi,
            step_size_max, step_tol,
            calc_initial_step_size(ss, step_size_factor),
            x, fitness(x, evaluator))
    end
end

# by default use RandomBound embedder
function GeneratingSetSearcher(problem::OptimizationProblem, parameters::Parameters)
    params = chain(GSSDefaultParameters, parameters)
    GeneratingSetSearcher(ProblemEvaluator(problem),
                        get(params, :DirectionGenerator, compass_search_directions(numdims(problem))),
                        RandomBound(search_space(problem)),
                        params[:RandomDirectionOrder],
                        params[:InitialStepSizeFactor], params[:StepSizeGamma],
                        params[:StepSizePhi], params[:StepSizeMax],
                        params[:DeltaTolerance])
end

# We also include the name of the direction generator.
name(opt::GeneratingSetSearcher) = "GeneratingSetSearcher($(typeof(opt.direction_gen)))"

has_converged(gss::GeneratingSetSearcher) = gss.step_size < gss.step_tol

function step!(gss::GeneratingSetSearcher)
    if has_converged(gss)
        # Restart from a random point
        gss.x = rand_individual(gss.search_space)
        gss.xfitness = fitness(gss.x, gss.evaluator)
        gss.step_size = calc_initial_step_size(gss.search_space, gss.step_size_factor)
    end

    # Get the directions for this iteration
    gss.k += 1
    directions = directions_for_k(gss.direction_gen, gss.k)

    # Set up order vector from which we will take the directions after possibly shuffling it
    order = collect(1:size(directions, 2))
    if gss.random_dir_order
        shuffle!(order)
    end

    # Check all directions to find a better point; default is that no one is found.
    found_better = false
    candidate = zeros(gss.n, 1)

    # Loop over directions until we find an improvement (or there are no more directions to check).
    for direction in order
        candidate = gss.x + gss.step_size .* directions[:, direction]
        apply!(gss.embed, candidate, gss.x)

        if is_better(candidate, gss.xfitness, gss.evaluator)
            found_better = true
            break
        end
    end

    if found_better
        gss.x = candidate
        gss.xfitness = last_fitness(gss.evaluator)
        gss.step_size *= gss.step_size_gamma
    else
        gss.step_size *= gss.step_size_phi
    end
    gss.step_size = min(gss.step_size, gss.step_size_max)

    return gss
end

ext: jl


text:
# init sensor figure
function sensorinit(s::mjSim)
    # set figure to default
    mjv_defaultFigure(s.figsensor)

    # set flags
    s.figsensor[].flg_extend = Cint(1)
    s.figsensor[].flg_barplot = Cint(1)

    s.figsensor[].title = str2vec("Sensor data", length(s.figsensor[].title))

    # y-tick nubmer format
    s.figsensor[].yformat = str2vec("%.0f", length(s.figsensor[].yformat))

    # grid size
    s.figsensor[].gridsize = [2, 3]

    # minimum range
    s.figsensor[].range = [[0 1],[-1 1]]
end

# update sensor figure
function sensorupdate(s::mjSim)
    m = s.m
    d = s.d

    # loop over sensors
    for n=1:m.nsensor
        # get info about this sensor
        cutoff = m.sensor_cutoff[n] > 0 ? m.sensor_cutoff[n] : 1.0
        adr = m.sensor_adr[n]
        dim = m.sensor_dim[n]

        # Print out the sensor name and current values
        print("$(mj_id2name(m, MJCore.mjOBJ_SENSOR, n-1)): ")
        for i=0:(dim-1)
            i > 0 && print(", ")
            print(round(d.sensordata[adr+i+1], digits = 3))
        end
        println()
    end
end

# show sensor figure
function sensorshow(s::mjSim, rect::mjrRect)
    # render figure on the right
    viewport = mjrRect(rect.width - rect.width / 4,
                       rect.bottom,
                       rect.width / 4,
                       rect.height / 3)
    mjr_figure(viewport, s.figsensor, s.con)
end


text:
using Random
include("ED.jl")

@testset "ED checks" begin
    void = State(0)
    up = State(1)
    down = State(2)
    updown = State(3)

    # create up at site 1
    _sign, s = create(void, 1, 1)
    @test _sign == 1.0 &&  s == up
    _sign, s = create(up, 1, 1)
    @test _sign == 0
    _sign, s = create(down, 1, 1)
    @test _sign == 1.0 &&  s == updown
    _sign, s = create(updown, 1, 1)
    @test _sign == 0

    # create down at site 1
    _sign, s = create(void, 1, 2)
    @test _sign == 1.0 &&  s == down
    _sign, s = create(up, 1, 2)
    @test _sign == -1.0 &&  s == updown
    _sign, s = create(down, 1, 2)
    @test _sign == 0
    _sign, s = create(updown, 1, 2)
    @test _sign == 0

    # annihilate up at site 1
    _sign, s = annihilate(void, 1, 1)
    @test _sign == 0
    _sign, s = annihilate(up, 1, 1)
    @test _sign == 1.0 &&  s == void
    _sign, s = annihilate(down, 1, 1)
    @test _sign == 0
    _sign, s = annihilate(updown, 1, 1)
    @test _sign == 1.0 &&  s == down

    # annihilate down at site 1
    _sign, s = annihilate(void, 1, 2)
    @test _sign == 0
    _sign, s = annihilate(up, 1, 2)
    @test _sign == 0
    _sign, s = annihilate(down, 1, 2)
    @test _sign == 1.0 &&  s == void
    _sign, s = annihilate(updown, 1, 2)
    @test _sign == -1.0 &&  s == up

    # check Greens consistency
    model = HubbardModelAttractive(2, 2, U = rand(), mu = rand(), t = rand()
    )
    H = HamiltonMatrix(model)
    for substate1 in 1:2, substate2 in 1:2
        for site1 in 1:model.l.sites, site2 in 1:model.l.sites
            G = expectation_value(
                Greens(site1, site2, substate1, substate2),
                H,
                N_sites = model.l.sites,
            )
            G_perm = expectation_value(
                Greens_permuted(site1, site2, substate1, substate2),
                H,
                N_sites = model.l.sites,
            )
            @test G ≈ G_perm

            # Test G(t, t) = G(0, 0) = G
            UTG = expectation_value(
                s -> annihilate(s, site2, substate2),
                s -> create(s, site1, substate1),
                H, 0.1, 0.1, N_sites=model.l.sites
            )
            @test UTG ≈ real(G) atol=1e-14
            UTG = expectation_value(
                s -> annihilate(s, site2, substate2),
                s -> create(s, site1, substate1),
                H, 0.7, 0.7, N_sites=model.l.sites
            )
            @test UTG ≈ real(G) atol=1e-14
        end
    end

    # Test G(t, 0) + G(0, beta - t) = 0
    UTG1 = calculate_Greens_matrix(H, 0.7, 0.0, model.l, beta=1.0)
    UTG2 = calculate_Greens_matrix(H, 0.0, 1.0 - 0.7, model.l, beta=1.0)
    @test UTG1 ≈ -UTG2 atol=1e-13
end

@testset "Repulsive/Attractive Hubbard Model (ED)" begin
    models = (
        HubbardModelRepulsive(2, 2, U = 1.0, t = 1.0),
        HubbardModelAttractive(2, 2, U = 1.0, mu = 1.0, t = 1.0)
    )

    for model in models
        @testset "$(typeof(model))" begin
            Random.seed!(123)
            dqmc = DQMC(
                model, beta=1.0, delta_tau = 0.1, safe_mult=5, recorder = Discarder,
                thermalization = 10_000, sweeps = 10_000
            )
            @info "Running DQMC ($(typeof(model).name)) β=$(dqmc.p.beta), 10k + 10k sweeps"

            dqmc[:G]    = greens_measurement(dqmc, model)
            dqmc[:Occs] = occupation(dqmc, model)
            dqmc[:CDC]  = charge_density_correlation(dqmc, model)
            dqmc[:Mx]   = magnetization(dqmc, model, :x)
            dqmc[:My]   = magnetization(dqmc, model, :y)
            dqmc[:Mz]   = magnetization(dqmc, model, :z)
            dqmc[:SDCx] = spin_density_correlation(dqmc, model, :x)
            dqmc[:SDCy] = spin_density_correlation(dqmc, model, :y)
            dqmc[:SDCz] = spin_density_correlation(dqmc, model, :z)
            dqmc[:PC]   = pairing_correlation(dqmc, model, K=4)

            # Unequal time
            l1s = [0, 3, 5, 7, 3, 0]
            l2s = [1, 7, 5, 2, 1, MonteCarlo.nslices(dqmc)]
            # l1s = [0, 3, 5, 0]
            # l2s = [1, 7, 5, MonteCarlo.nslices(dqmc)]
            dqmc[:UTG1] = greens_measurement(dqmc, model, GreensAt{l2s[1], l1s[1]})
            dqmc[:UTG2] = greens_measurement(dqmc, model, GreensAt{l2s[2], l1s[2]})
            dqmc[:UTG3] = greens_measurement(dqmc, model, GreensAt{l2s[3], l1s[3]})
            dqmc[:UTG4] = greens_measurement(dqmc, model, GreensAt{l2s[4], l1s[4]})
            dqmc[:UTG5] = greens_measurement(dqmc, model, GreensAt{l2s[5], l1s[5]})
            dqmc[:UTG6] = greens_measurement(dqmc, model, GreensAt{l2s[6], l1s[6]})

            dqmc[:CDS]  = charge_density_susceptibility(dqmc, model)
            dqmc[:SDSx] = spin_density_susceptibility(dqmc, model, :x)
            dqmc[:SDSy] = spin_density_susceptibility(dqmc, model, :y)
            dqmc[:SDSz] = spin_density_susceptibility(dqmc, model, :z)
            dqmc[:PS]   = pairing_susceptibility(dqmc, model, K=4)
            dqmc[:CCS]  = current_current_susceptibility(dqmc, model, K=4)

            # MonteCarlo.enable_benchmarks()

            @time run!(dqmc, verbose=false)

            # Absolute tolerance from Trotter decompositon
            atol = 2dqmc.p.delta_tau^2
            rtol = 2dqmc.p.delta_tau^2
            N = length(lattice(model))

            @info "Running ED"
            @time begin
                H = HamiltonMatrix(model)

                # G_DQMC is smaller because it doesn't differentiate between spin up/down
                @testset "Greens" begin
                    G_DQMC = mean(dqmc.measurements[:G])
                    occs = mean(dqmc.measurements[:Occs])                                   # measuring
                    # occs2 = mean(MonteCarlo.occupations(dqmc.measurements[:Greens]))            # wrapping
                    # occs3 = mean(MonteCarlo.OccupationMeasurement(dqmc.measurements[:Greens]))  # copying
                    G_ED = calculate_Greens_matrix(H, model.l, beta = dqmc.p.beta)
                    for i in 1:size(G_DQMC, 1), j in 1:size(G_DQMC, 2)
                        @test isapprox(G_DQMC[i, j], G_ED[i, j], atol=atol, rtol=rtol)
                    end
                    for i in 1:size(G_DQMC, 1)
                        @test isapprox(occs[i],  1 - G_ED[i, i], atol=atol, rtol=rtol)
                        # @test isapprox(occs2[i], 1 - G_ED[i, i], atol=atol, rtol=rtol)
                        # @test isapprox(occs3[i], 1 - G_ED[i, i], atol=atol, rtol=rtol)
                    end
                end

                @testset "Charge Density Correlation" begin
                    CDC = mean(dqmc.measurements[:CDC])
                    ED_CDC = zeros(size(CDC))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:CDC], dqmc, model)
                        ED_CDC[dir] += expectation_value(
                            charge_density_correlation(trg, src),
                            H, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_CDC/N ≈ CDC atol=atol rtol=rtol
                end

                @testset "Magnetization x" begin
                    Mx = mean(dqmc.measurements[:Mx])
                    for site in 1:length(Mx)
                        ED_Mx = expectation_value(m_x(site), H, beta = dqmc.p.beta, N_sites = N)
                        @test ED_Mx ≈ Mx[site] atol=atol rtol=rtol
                    end
                end
                @testset "Magnetization y" begin
                    My = mean(dqmc.measurements[:My])
                    for site in 1:length(My)
                        ED_My = expectation_value(m_y(site), H, beta = dqmc.p.beta, N_sites = N)
                        @test ED_My ≈ My[site] atol=atol rtol=rtol
                    end
                end
                @testset "Magnetization z" begin
                    Mz = mean(dqmc.measurements[:Mz])
                    ΔMz = std_error(dqmc.measurements[:Mz])
                    for site in 1:length(Mz)
                        ED_Mz = expectation_value(m_z(site), H, beta = dqmc.p.beta, N_sites = N)
                        @test ED_Mz ≈ Mz[site] atol=atol+ΔMz[site] rtol=rtol
                    end
                end

                @testset "Spin density correlation x" begin
                    SDCx = mean(dqmc.measurements[:SDCx])
                    ED_SDCx = zeros(ComplexF64, size(SDCx))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:SDCx], dqmc, model)
                        ED_SDCx[dir] += expectation_value(
                            spin_density_correlation_x(trg, src),
                            H, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_SDCx/N ≈ SDCx atol=atol rtol=rtol
                end
                @testset "Spin density correlation y" begin
                    SDCy = mean(dqmc.measurements[:SDCy])
                    ED_SDCy = zeros(ComplexF64, size(SDCy))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:SDCy], dqmc, model)
                        ED_SDCy[dir] += expectation_value(
                            spin_density_correlation_y(trg, src),
                            H, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_SDCy/N ≈ SDCy atol=atol rtol=rtol
                end
                @testset "Spin density correlation z" begin
                    SDCz = mean(dqmc.measurements[:SDCz])
                    ED_SDCz = zeros(ComplexF64, size(SDCz))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:SDCz], dqmc, model)
                        ED_SDCz[dir] += expectation_value(
                            spin_density_correlation_z(trg, src),
                            H, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_SDCz/N ≈ SDCz atol=atol rtol=rtol
                end

                @testset "Pairing Correlation" begin
                    PC = mean(dqmc.measurements[:PC])
                    ED_PC = zeros(ComplexF64, size(PC))
                    for (dirs, src1, trg1, src2, trg2) in
                            MonteCarlo.EachLocalQuadByDistance{4}(dqmc, model)
                        ED_PC[dirs] += expectation_value(
                            pairing_correlation(src1, trg1, src2, trg2),
                            H, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_PC/N ≈ PC atol=atol rtol=rtol
                end

                ################################################################
                ### Unequal Time
                ################################################################

                @testset "Unequal Time Greens" begin
                    for (i, tau1, tau2) in zip(eachindex(l1s), 0.1l1s, 0.1l2s)
                        UTG = mean(dqmc.measurements[Symbol(:UTG, i)])
                        # ΔUTG = std_error(dqmc.measurements[Symbol(:UTG, i)])
                        M = size(UTG, 1)
                        ED_UTG = calculate_Greens_matrix(H, tau2, tau1, model.l, beta=dqmc.p.beta)

                        @testset "[$i] $tau1 -> $tau2" begin
                            for k in 1:M, l in 1:M
                                @test isapprox(UTG[k, l], ED_UTG[k, l], atol=atol, rtol=rtol)
                            end
                            # println("[$i] $tau1 -> $tau2")
                            # display(UTG)
                            # println()
                            # display(ED_UTG)
                            # println()
                        end
                    end
                end

                # TODO SDS, CDS, PS
                # So all of these are off by more than an order of magnitude...
                # why? G0l needs to be added...
                @testset "Charge Density Susceptibility" begin
                    CDS = mean(dqmc.measurements[:CDS])
                    ED_CDS = zeros(size(CDS))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:CDS], dqmc, model)
                        ED_CDS[dir] += expectation_value_integrated(
                            number_operator(trg), number_operator(src), H,
                            step = dqmc.p.delta_tau, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_CDS/N ≈ CDS atol=atol rtol=rtol
                end

                # SDS
                @testset "Spin density Susceptibility x" begin
                    SDSx = mean(dqmc.measurements[:SDSx])
                    ED_SDSx = zeros(ComplexF64, size(SDSx))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:SDSx], dqmc, model)
                        ED_SDSx[dir] += expectation_value_integrated(
                            m_x(trg), m_x(src), H, step = dqmc.p.delta_tau,
                            beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_SDSx/N ≈ SDSx atol=atol rtol=rtol
                end
                @testset "Spin density Susceptibility y" begin
                    SDSy = mean(dqmc.measurements[:SDSy])
                    ED_SDSy = zeros(ComplexF64, size(SDSy))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:SDSy], dqmc, model)
                        ED_SDSy[dir] += expectation_value_integrated(
                            m_y(trg), m_y(src), H, step = dqmc.p.delta_tau,
                            beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_SDSy/N ≈ SDSy atol=atol rtol=rtol
                end
                @testset "Spin density Susceptibility z" begin
                    SDSz = mean(dqmc.measurements[:SDSz])
                    ED_SDSz = zeros(ComplexF64, size(SDSz))
                    for (dir, src, trg) in MonteCarlo.lattice_iterator(dqmc[:SDSz], dqmc, model)
                        ED_SDSz[dir] += expectation_value_integrated(
                            m_z(trg), m_z(src), H, step = dqmc.p.delta_tau,
                            beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_SDSz/N ≈ SDSz atol=atol rtol=rtol
                end

                @testset "Pairing Susceptibility" begin
                    PS = mean(dqmc.measurements[:PS])
                    ED_PS = zeros(Float64, size(PS))
                    for (dirs, src1, trg1, src2, trg2) in
                            MonteCarlo.EachLocalQuadByDistance{4}(dqmc, model)
                        ED_PS[dirs] += expectation_value_integrated(
                            state -> begin
                                sign1, _state = annihilate(state, trg1, DOWN)
                                sign2, _state = annihilate(_state, src1, UP)
                                p = sign1*sign2
                                p, _state
                            end,
                            state -> begin
                                sign1, _state = create(state, src2, UP)
                                sign2, _state = create(_state, trg2, DOWN)
                                p = sign1*sign2
                                p, _state
                            end,
                            H, step = dqmc.p.delta_tau, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_PS/N ≈ PS atol=atol rtol=rtol
                end

                @testset "Current Current Susceptibility" begin
                    CCS = mean(dqmc.measurements[:CCS])
                    ED_CCS = zeros(Float64, size(CCS))
                    T = dqmc.s.hopping_matrix
                    for (dirs, src1, trg1, src2, trg2) in
                            MonteCarlo.EachLocalQuadBySyncedDistance{4}(dqmc, model)
                        ED_CCS[dirs] += expectation_value_integrated(
                            current_density(src1, trg1, T),
                            current_density(src2, trg2, T),
                            H, step = dqmc.p.delta_tau, beta = dqmc.p.beta, N_sites = N
                        )
                    end
                    @test ED_CCS/N ≈ CCS atol=atol rtol=rtol
                end
            end
        end
    end
end


text:
using LinearAlgebraicRepresentation
using Plasm, DataStructures
using PyCall
Lar = LinearAlgebraicRepresentation
const colorkey = ["orange", "red", "green", "blue", "cyan", "magenta", "yellow",
	"white", "black", "purple", "gray", "brown"]
#=
#1) Lar.depth_first_search:  from Tarjan, '72 (SIAM J. Computing, 1)
#2) Lar.edge_biconnect:  	from Tarjan, '72	(SIAM J. Computing, 1)
#3) Lar.biconnectedComponent:  from Hopcroft, Tarjan, '73 (ACM Communications )
=#

# Example 1 (DFS)

# Graph description as set of edges via EV array
EV = [[1,2],[2,3],[3,4],[4,5],[5,6],[1,6],[6,7],[4,7],[2,7],[5,8],[1,8],[3,8]]
V = hcat([[2.,1],[2,2],[1,2],[0,3],[0,0],[3,0],[3,3],[1,1]]...)

spanningtree, fronds = Lar.depth_first_search(EV) # spanning tree: default to node 1
hpc1 = Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)], spanningtree]))
hpc2 = Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)], fronds]))
Plasm.view([ Plasm.color("cyan")(hpc1) , Plasm.color("magenta")(hpc2) ])

spanningtree, fronds = Lar.depth_first_search(EV,5) # spanning tree:  node 5
hpc1 = Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)], spanningtree]))
hpc2 = Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)], fronds]))
Plasm.view([ Plasm.color("cyan")(hpc1) , Plasm.color("magenta")(hpc2) ])

# Example 2 (edge_biconnect vs biconnectedComponent)

V = hcat([[1.,2],[1,3],[0,3],[1,4],[2,3],[2,2],[1,1],[1,0],[0,1]]...)
EV = [[3,4],[2,4],[2,3],[2,5],[1,2],[5,6],[1,6],[1,7],[7,9],[8,9],[7,8]]

EVs = Lar.edge_biconnect(EV::Lar.Cells) # 2-edge-connected components
hpcs = [ Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)],EVs[i]])) for i=1:length(EVs)]
Plasm.view([ Plasm.color(colorkey[k])(hpcs[k]) for k=1:length(hpcs) ])

V,EVs = Lar.biconnectedComponent((V,EV::Lar.Cells)) # 2-connected components
hpcs = [ Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)],EVs[i]])) for i=1:length(EVs)]
Plasm.view([ Plasm.color(colorkey[k])(hpcs[k]) for k=1:length(hpcs) ])

copEV = Lar.characteristicMatrix(EV)
FE = Lar.Arrangement.biconnected_components(copEV) # 2-connected components (H & T)
EVs = [[EV[e] for e in F] for F in FE]
hpcs = [ Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)],EVs[i]])) for i=1:length(EVs)]
Plasm.view([ Plasm.color(colorkey[k])(hpcs[k]) for k=1:length(hpcs) ])

# Example 3 (edge_biconnect vs biconnectedComponent vs biconnected_components)

V = hcat([[0,4],[4,4],[4,0],[2,0],[0,0],[2,2],[1,2],[3,2],[1,3],[3,3]]...)
EV = [[1,2],[9,10],[1,5],[7,9],[8,10],[2,3],[6,7],[6,8],[4,6],[4,5],[3,4]]

EVs = Lar.edge_biconnect(EV::Lar.Cells) # 2-edge-connected components (T)
hpcs = [ Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)],EVs[i]])) for i=1:length(EVs)]
Plasm.view([ Plasm.color(colorkey[k])(hpcs[k]) for k=1:length(hpcs) ])

V,EVs = Lar.biconnectedComponent((V,EV::Lar.Cells)) # 2-connected components (H & T)
hpcs = [ Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)],EVs[i]])) for i=1:length(EVs)]
Plasm.view([ Plasm.color(colorkey[k])(hpcs[k]) for k=1:length(hpcs) ])

copEV = Lar.characteristicMatrix(EV)
FE = Lar.Arrangement.biconnected_components(copEV) # 2-connected components (H & T)
EVs = [[EV[e] for e in F] for F in FE]
hpcs = [ Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)],EVs[i]])) for i=1:length(EVs)]
Plasm.view([ Plasm.color(colorkey[k])(hpcs[k]) for k=1:length(hpcs) ])

# Example 4

V = hcat([[0.,0],[1,0],[1,1],[0,1],[2,1]]...);
EV = [[1,2],[2,3],[3,4],[4,1],[1,5]];
model = V,EV
W,EW = Lar.fragmentlines(model)
V,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells)) # 2-connected components (H & T)
Plasm.view(Plasm.numbering()((V,[[[k] for k=1:size(V,2)],EVs[1]])))


text:
using ITensors
using ITensorInfiniteMPS
using Random

Random.seed!(1234)

N = 100
s = siteinds("S=1/2", N)

J = 1.0
h = 1.5
model = Model("ising")
H = MPO(model, s; J=J, h=h)
ψ0 = randomMPS(s)

sweeps = Sweeps(10)
maxdim!(sweeps, 10)
cutoff!(sweeps, 1E-10)
energy, ψ = dmrg(H, ψ0, sweeps)
@show energy / N

n1 = N ÷ 2
n2 = n1 + 1

ψ = orthogonalize(ψ, n1)

println("\nFinite MPS energy on bond (n1, n2) = $((n1, n2))")

ϕ12 = ψ[n1] * ψ[n2]
ham12 = ITensor(model, s[n1], s[n2]; J=J, h=h)
@show (ϕ12 * ham12 * prime(dag(ϕ12), commoninds(ϕ12, ham12)))[]

# Get an infinite MPS that approximates
# the finite MPS in the bulk.
# `Nsites` sets the number of sites in the unit
# cell of the infinite MPS.
Nsites = 2
# XXX: not outputing a right orthogonal ψ∞.AR
ψ∞ = infinitemps_approx(ψ; nsites=Nsites, nsweeps=10)

# XXX: this is segfaulting for some reason
#ψ∞ = orthogonalize(ψ∞.AL, :)
C, ψR = ITensorInfiniteMPS.right_orthogonalize(ψ∞.AL)
ψ∞ = InfiniteCanonicalMPS(ψ∞.AL, C, ψR)

println("\nInfinite MPS energy on a few different bonds")
println("Infinite MPS has unit cell of size $Nsites")

# Compute a few energies
for n1 in 1:3
  global n2 = n1 + 1
  @show n1, n2
  global ϕ12 = ψ∞.AL[n1] * ψ∞.AL[n2] * ψ∞.C[n2]
  global s1 = siteind(ψ∞, n1)
  global s2 = siteind(ψ∞, n2)
  global ham12 = ITensor(model, s1, s2; J=J, h=h)
  global ϕ12ᴴ = prime(dag(ϕ12), commoninds(ϕ12, ham12))
  @show (ϕ12 * ham12 * ϕ12ᴴ)[] / norm(ϕ12)
end

s¹ = siteinds(ψ∞, Cell(1))

Σ∞h = InfiniteITensorSum(model, s¹; J=J, h=h)

ψ∞_opt = vumps(Σ∞h, ψ∞)

χ = 6
ψ∞ = InfiniteMPS(only.(s¹); space=χ)
randn!.(ψ∞)

ψ∞ = orthogonalize(ψ∞, :)
@show norm(prod(ψ∞.AL[1:Nsites]) * ψ∞.C[Nsites] - ψ∞.C[0] * prod(ψ∞.AR[1:Nsites]))

vumps(Σ∞h, ψ∞; niter=10)


text:
"""
    niw_hyperparams(κ::Float32, m::AbstractArray{Float32}, ν::Float32, ψ::AbstractArray{Float32})

[Normal Inverse Wishart](https://en.wikipedia.org/wiki/Normal-inverse-Wishart_distribution)
"""
struct niw_hyperparams <: distribution_hyper_params
    κ::Float32
    m::AbstractArray{Float64}
    ν::Float32
    ψ::AbstractArray{Float64}
end

mutable struct niw_sufficient_statistics <: sufficient_statistics
    N::Float32
    points_sum::AbstractArray{Float64,1}
    S::AbstractArray{Float64,2}
end

function calc_posterior(prior::niw_hyperparams, suff_statistics::niw_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    κ = prior.κ + suff_statistics.N
    ν = prior.ν + suff_statistics.N
    m = (prior.m.*prior.κ + suff_statistics.points_sum) / κ
    ψ = (prior.ν * prior.ψ + prior.κ*prior.m*prior.m' -κ*m*m'+ suff_statistics.S) / ν
    ψ = Matrix(Symmetric(ψ))
    ψ = (ψ+ψ')/2
    return niw_hyperparams(κ,m,ν,ψ)
end

function sample_distribution(hyperparams::niw_hyperparams)
    Σ = rand(Distributions.InverseWishart(hyperparams.ν, hyperparams.ν * hyperparams.ψ))
    μ = rand(Distributions.MvNormal(hyperparams.m, Σ / hyperparams.κ))
    invΣ = inv(Σ)
    chol = cholesky(Hermitian(invΣ))
    return mv_gaussian(μ, Σ, invΣ, logdet(Σ), chol.U)
end

function create_sufficient_statistics(
    hyper::niw_hyperparams,
    posterior::niw_hyperparams,
    points::AbstractArray{Float32,2},
    pts_to_group = 0,
)
    if size(points, 2) == 0
        return niw_sufficient_statistics(
            size(points, 2),
            zeros(Float32, length(hyper.m)),
            zeros(Float32, length(hyper.m), length(hyper.m)),
        )
    end
    pts = Array{Float64}(points)
    points_sum = sum(pts, dims = 2)[:]
    S = pts * pts'
    S = 0.5*(S+S')
    return niw_sufficient_statistics(size(points,2),points_sum,S)
end

function log_marginal_likelihood(
    hyper::niw_hyperparams,
    posterior_hyper::niw_hyperparams,
    suff_stats::niw_sufficient_statistics,
)
    D = length(suff_stats.points_sum)
    logpi = log(pi)
    return -suff_stats.N * D * 0.5 * logpi +
           log_multivariate_gamma(posterior_hyper.ν / 2, D) -
           log_multivariate_gamma(hyper.ν / 2, D) +
           (hyper.ν / 2) * (D * log(hyper.ν) + logdet(hyper.ψ)) -
           (posterior_hyper.ν / 2) *
           (D * log(posterior_hyper.ν) + logdet(posterior_hyper.ψ)) +
           (D / 2) * (log(hyper.κ / posterior_hyper.κ))
end

function aggregate_suff_stats(
    suff_l::niw_sufficient_statistics,
    suff_r::niw_sufficient_statistics,
)
    return niw_sufficient_statistics(
        suff_l.N + suff_r.N,
        suff_l.points_sum + suff_r.points_sum,
        suff_l.S + suff_r.S,
    )
end


text:
"""
	SourceParameters

Custom type defining the source parameters of a Fourier spectrum.

Constructed with signature `SourceParameters{S<:Float64, T<:Real}` with fields:
- `Δσ::T` is the stress parameter in bars
- `RΘϕ::S` is the radiation pattern
- `V::S` is the partition factor (for splitting to horizontal components)
- `F::S` is the free surface factor
- `β::S` is the source velocity in units of km/s
- `ρ::S` is the source density in units of t/m³ or g/cm³
- `model::Symbol` identifies the type of source spectrum (`:Brune`, `:Atkinson_Silva_2000`)
"""
struct SourceParameters{S<:Float64,T<:Real}
    # source parameters
    Δσ::T# stressParameter
    RΘϕ::S                            # radiationPattern
    V::S                              # partitionFactor
    F::S                              # freeSurfaceFactor
    β::S                              # sourceVelocity
    ρ::S                              # sourceDensity
    model::Symbol    # source spectrum model
end

SourceParameters(Δσ::T) where {T<:Float64} = SourceParameters(Δσ, 0.55, 1.0 / sqrt(2.0), 2.0, 3.5, 2.75, :Brune)
SourceParameters(Δσ::T) where {T<:Real} = SourceParameters(Δσ, 0.55, 1.0 / sqrt(2.0), 2.0, 3.5, 2.75, :Brune)
SourceParameters(Δσ::T, model::Symbol) where {T<:Real} = SourceParameters(Δσ, 0.55, 1.0 / sqrt(2.0), 2.0, 3.5, 2.75, model)
SourceParameters(Δσ::T, βs::S, ρs::S) where {S<:Float64,T<:Real} = SourceParameters(Δσ, 0.55, 1.0 / sqrt(2.0), 2.0, βs, ρs, :Brune)
SourceParameters(Δσ::T, βs::T, ρs::T) where {T<:Float64} = SourceParameters(Δσ, 0.55, 1.0 / sqrt(2.0), 2.0, βs, ρs, :Brune)

"""
	get_parametric_type(src::SourceParameters{S,T}) where {S,T} = T

Extract type of `T` from parametric `SourceParameters` struct
"""
get_parametric_type(src::SourceParameters{S,T}) where {S,T} = T

"""
	GeometricSpreadingParameters

Struct for geometric spreading parameters.
Holds fields:
- `Rrefi` are reference distances, these are `<:Real` but will generally be `Float64` values
- `γconi` are constant spreading rates, meaning that they will not be free for AD purposes
- `γvari` are variable spreading rates, meaning that they can be represented as `Dual` numbers for AD
- `γfree` is a vector of `Bool` instances, or a `BitVector` that indicates which segments are constant or variable. Variable spreading rates are given `1` or `true`
- `model` is a symbol defining the type of spreading model `:Piecewise`, `:CY14`, `:CY14mod`
"""
struct GeometricSpreadingParameters{S<:Real,T<:Real,U<:AbstractVector{Bool}}
    Rrefi::Vector{S}
    γconi::Vector{S}
    γvari::Vector{T}
    γfree::U
    model::Symbol
end

GeometricSpreadingParameters(Rrefi::Vector{T}, γconi::Vector{T}) where {T} = GeometricSpreadingParameters{T,T,BitVector}(Rrefi, γconi, Vector{T}(), BitVector(zeros(length(γconi))), :Piecewise)
GeometricSpreadingParameters(Rrefi::Vector{T}, γconi::Vector{T}, model::Symbol) where {T} = GeometricSpreadingParameters{T,T,BitVector}(Rrefi, γconi, Vector{T}(), BitVector(zeros(length(γconi))), model)
GeometricSpreadingParameters(Rrefi::Vector{S}, γvari::Vector{T}) where {S<:Float64,T<:Dual} = GeometricSpreadingParameters{S,T,BitVector}(Rrefi, Vector{S}(), γvari, BitVector(ones(length(γvari))), :Piecewise)
GeometricSpreadingParameters(Rrefi::Vector{S}, γvari::Vector{T}, model::Symbol) where {S<:Float64,T<:Dual} = GeometricSpreadingParameters{S,T,BitVector}(Rrefi, Vector{S}(), γvari, BitVector(ones(length(γvari))), model)

"""
	get_parametric_type(geo::GeometricSpreadingParameters{S,T,U}) where {S,T,U} = T

Extract type of `T` from parametric `GeometricSpreadingParameters` struct
"""
get_parametric_type(geo::GeometricSpreadingParameters{S,T,U}) where {S,T,U} = T

"""
	NearSourceSaturationParameters

Struct for near-source saturation parameters. Mimic structure of the `GeometricSpreadingParameters` struct.
Holds fields:
- `mRefi` reference magnitudes
- `hconi` constrained coefficients, not free for AD purposes
- `hvari` variable coefficients, free for AD purposes
- `hfree` is a vector of `Bool` instances, or a `BitVector` indicating which parameters are constant or variable
- `exponent` is the exponent used within equivalent point-source distance calculations: ``r_{ps} = \\left[r_{rup}^n + h(m)^n\\right]^{1/n}``
- `model` is a symbol defining the type of saturation model:
  - `:BT15` is Boore & Thompson (2015)
  - `:YA15` is Yenier & Atkinson (2015)
  - `:CY14` is average Chiou & Youngs (2014)
  - `:None` returns zero saturation length
  - `:ConstantConstrained` is a fixed saturation length not subject to AD operations
  - `:ConstantVariable` is a fixed saturation length that is subject to AD operations (i.e., is a `<:Dual`)

"""
struct NearSourceSaturationParameters{S<:Real,T<:Real,U<:AbstractVector{Bool}}
    mRefi::Vector{S}
    hconi::Vector{S}
    hvari::Vector{T}
    hfree::U
    exponent::Int
    model::Symbol
end

NearSourceSaturationParameters(model::Symbol) = NearSourceSaturationParameters(Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), BitVector(), 2, model)
NearSourceSaturationParameters(exponent::Int, model::Symbol) = NearSourceSaturationParameters(Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), BitVector(), exponent, model)
NearSourceSaturationParameters(mRefi::Vector{T}, hconi::Vector{T}, model::Symbol) where {T} = NearSourceSaturationParameters(mRefi, hconi, Vector{T}(), BitVector(undef, length(hconi)), 2, model)
NearSourceSaturationParameters(mRefi::Vector{S}, hvari::Vector{T}, model::Symbol) where {S,T} = NearSourceSaturationParameters(mRefi, Vector{S}(), hvari, BitVector(ones(length(hvari))), 2, model)

NearSourceSaturationParameters(mRefi::Vector{T}, hconi::Vector{T}) where {T} = NearSourceSaturationParameters(mRefi, hconi, Vector{T}(), BitVector(undef, length(hconi)), 2, :FullyConstrained)
NearSourceSaturationParameters(mRefi::Vector{S}, hvari::Vector{T}) where {S,T} = NearSourceSaturationParameters(mRefi, Vector{S}(), hvari, BitVector(ones(length(hvari))), 2, :FullyVariable)

# specialisation for a constant saturation term
NearSourceSaturationParameters(hcon::Float64) = NearSourceSaturationParameters(Vector{Float64}(), [hcon], Vector{Float64}(), BitVector(undef, 1), 2, :ConstantConstrained)
NearSourceSaturationParameters(hvar::T) where {T<:Dual} = NearSourceSaturationParameters(Vector{Float64}(), Vector{Float64}(), Vector{T}([hvar]), BitVector([1]), 2, :ConstantVariable)
# with exponents
NearSourceSaturationParameters(hcon::Float64, exponent::Int) = NearSourceSaturationParameters(Vector{Float64}(), [hcon], Vector{Float64}(), BitVector(undef, 1), exponent, :ConstantConstrained)
NearSourceSaturationParameters(hvar::T, exponent::Int) where {T<:Dual} = NearSourceSaturationParameters(Vector{Float64}(), Vector{Float64}(), Vector{T}([hvar]), BitVector([1]), exponent, :ConstantVariable)

"""
	get_parametric_type(sat::NearSourceSaturationParameters{S,T,U}) where {S,T,U} = T

Extract type of `T` from parametric `NearSourceSaturationParameters` struct
"""
get_parametric_type(sat::NearSourceSaturationParameters{S,T,U}) where {S,T,U} = T

"""
	AnelasticAttenuationParameters

Struct for anelastic attenuation parameters.
Holds fields:
- `Q0` quality factor at 1 Hz
- `η` quality exponent ∈ [0,1)
- `cQ` velocity (km/s) along propagation path used to determine `Q(f)`
- `rmetric` is a symbol `:Rrup` or `:Rps` to define which distance metric is used for anelastic attenuation
"""
struct AnelasticAttenuationParameters{S<:Real,T<:Real,U<:AbstractVector{Bool}}
    Rrefi::Vector{S}
    Q0coni::Vector{S}
    Q0vari::Vector{T}
    ηconi::Vector{S}
    ηvari::Vector{T}
    cQ::Vector{Float64}
    Qfree::U
    ηfree::U
    rmetric::Symbol
end

AnelasticAttenuationParameters(Q0::T) where {T<:Float64} = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{Float64}(), [0.0], Vector{Float64}(), [3.5], BitVector(zeros(1)), BitVector(zeros(1)), :Rps)
AnelasticAttenuationParameters(Q0::T) where {T<:Real} = AnelasticAttenuationParameters([0.0, Inf], Vector{Float64}(), [Q0], [0.0], Vector{T}(), [3.5], BitVector(ones(1)), BitVector(zeros(1)), :Rps)
AnelasticAttenuationParameters(Q0::T, η::T) where {T<:Float64} = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{T}(), [η], Vector{T}(), [3.5], BitVector(zeros(1)), BitVector(zeros(1)), :Rps)
AnelasticAttenuationParameters(Q0::T, η::T, rmetric::Symbol) where {T<:Float64} = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{T}(), [η], Vector{T}(), [3.5], BitVector(zeros(1)), BitVector(zeros(1)), rmetric)
AnelasticAttenuationParameters(Q0::T, η::T) where {T<:Real} = AnelasticAttenuationParameters([0.0, Inf], Vector{Float64}(), [Q0], Vector{Float64}(), [η], [3.5], BitVector(ones(1)), BitVector(ones(1)), :Rps)
AnelasticAttenuationParameters(Q0::T, η::T, rmetric::Symbol) where {T<:Real} = AnelasticAttenuationParameters([0.0, Inf], Vector{Float64}(), [Q0], Vector{Float64}(), [η], [3.5], BitVector(ones(1)), BitVector(ones(1)), rmetric)
AnelasticAttenuationParameters(Q0::S, η::T) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{T}(), Vector{S}(), [η], [3.5], BitVector(zeros(1)), BitVector(ones(1)), :Rps)
AnelasticAttenuationParameters(Q0::S, η::T, rmetric::Symbol) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{T}(), Vector{S}(), [η], [3.5], BitVector(zeros(1)), BitVector(ones(1)), rmetric)
AnelasticAttenuationParameters(Q0::S, η::T) where {S<:Real,T<:Float64} = AnelasticAttenuationParameters([0.0, Inf], Vector{T}(), [Q0], [η], Vector{S}(), [3.5], BitVector(ones(1)), BitVector(zeros(1)), :Rps)
AnelasticAttenuationParameters(Q0::S, η::T, rmetric::Symbol) where {S<:Real,T<:Float64} = AnelasticAttenuationParameters([0.0, Inf], Vector{T}(), [Q0], [η], Vector{S}(), [3.5], BitVector(ones(1)), BitVector(zeros(1)), rmetric)

AnelasticAttenuationParameters(Q0, η, cQ::Float64) = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{Float64}(), [η], Vector{Float64}(), [cQ], BitVector(zeros(1)), BitVector(zeros(1)), :Rps)
AnelasticAttenuationParameters(Q0, η, cQ::Float64, rmetric::Symbol) = AnelasticAttenuationParameters([0.0, Inf], [Q0], Vector{Float64}(), [η], Vector{Float64}(), [cQ], BitVector(zeros(1)), BitVector(zeros(1)), rmetric)

AnelasticAttenuationParameters(Rrefi::Vector{T}, Q0coni::Vector{T}) where {T<:Float64} = AnelasticAttenuationParameters{T,T,BitVector}(Rrefi, Q0coni, Vector{T}(), zeros(T, length(Q0coni)), Vector{T}(), 3.5 * ones(T, length(Q0coni)), BitVector(zeros(length(Q0coni))), BitVector(zeros(length(Q0coni))), :Rps)
AnelasticAttenuationParameters(Rrefi::Vector{T}, Q0coni::Vector{T}, rmetric::Symbol) where {T<:Float64} = AnelasticAttenuationParameters{T,T,BitVector}(Rrefi, Q0coni, Vector{T}(), zeros(T, length(Q0coni)), Vector{T}(), 3.5 * ones(T, length(Q0coni)), BitVector(zeros(length(Q0coni))), BitVector(zeros(length(Q0coni))), rmetric)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0vari::Vector{T}) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Vector{S}(), Q0vari, zeros(S, length(Q0vari)), Vector{S}(), 3.5 * ones(S, length(Q0vari)), BitVector(ones(length(Q0vari))), BitVector(zeros(length(Q0vari))), :Rps)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0vari::Vector{T}, rmetric::Symbol) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Vector{S}(), Q0vari, zeros(S, length(Q0vari)), Vector{S}(), 3.5 * ones(S, length(Q0vari)), BitVector(ones(length(Q0vari))), BitVector(zeros(length(Q0vari))), rmetric)

AnelasticAttenuationParameters(Rrefi::Vector{T}, Q0coni::Vector{T}, ηconi::Vector{T}) where {T<:Float64} = AnelasticAttenuationParameters{T,T,BitVector}(Rrefi, Q0coni, Vector{T}(), ηconi, Vector{T}(), 3.5 * ones(T, length(Q0coni)), BitVector(zeros(length(Q0coni))), BitVector(zeros(length(ηconi))), :Rps)
AnelasticAttenuationParameters(Rrefi::Vector{T}, Q0coni::Vector{T}, ηconi::Vector{T}, rmetric::Symbol) where {T<:Float64} = AnelasticAttenuationParameters{T,T,BitVector}(Rrefi, Q0coni, Vector{T}(), ηconi, Vector{T}(), 3.5 * ones(T, length(Q0coni)), BitVector(zeros(length(Q0coni))), BitVector(zeros(length(ηconi))), rmetric)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0vari::Vector{T}, ηvari::Vector{T}) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Vector{S}(), Q0vari, zeros(S, length(Q0vari)), ηvari, 3.5 * ones(S, length(Q0vari)), BitVector(ones(length(Q0vari))), BitVector(ones(length(ηvari))), :Rps)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0vari::Vector{T}, ηvari::Vector{T}, rmetric::Symbol) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Vector{S}(), Q0vari, zeros(S, length(Q0vari)), ηvari, 3.5 * ones(S, length(Q0vari)), BitVector(ones(length(Q0vari))), BitVector(ones(length(Q0vari))), rmetric)

AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0coni::Vector{S}, ηvari::Vector{T}) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Q0coni, Vector{T}(), Vector{S}(), ηvari, 3.5 * ones(S, length(Q0coni)), BitVector(zeros(length(Q0coni))), BitVector(ones(length(ηvari))), :Rps)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0coni::Vector{S}, ηvari::Vector{T}, rmetric::Symbol) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Q0coni, Vector{T}(), Vector{S}(), ηvari, 3.5 * ones(S, length(Q0coni)), BitVector(zeros(length(Q0coni))), BitVector(ones(length(ηvari))), rmetric)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0vari::Vector{T}, ηconi::Vector{S}) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Vector{S}(), Q0vari, ηconi, Vector{T}(), 3.5 * ones(S, length(Q0vari)), BitVector(ones(length(Q0vari))), BitVector(zeros(length(ηconi))), :Rps)
AnelasticAttenuationParameters(Rrefi::Vector{S}, Q0vari::Vector{T}, ηconi::Vector{S}, rmetric::Symbol) where {S<:Float64,T<:Real} = AnelasticAttenuationParameters{S,T,BitVector}(Rrefi, Vector{S}(), Q0vari, ηconi, Vector{T}(), 3.5 * ones(S, length(Q0vari)), BitVector(ones(length(Q0vari))), BitVector(zeros(length(ηconi))), rmetric)

"""
	get_parametric_type(anelastic::AnelasticAttenuationParameters{S,T,U}) where {S,T,U}

Extract type most elaborate type `S` or `T` from parametric `AnelasticAttenuationParameters` struct
"""
function get_parametric_type(anelastic::AnelasticAttenuationParameters{S,T,U}) where {S,T,U}
    if S <: Float64
        if T <: Float64
            return S
        else
            return T
        end
    else
        return S
    end
end

"""
	PathParameters

Custom type defining the path parameters of a Fourier spectrum.
Consists of three other custom structs
- `geometric` is a `GeometricSpreadingParameters` type
- `saturation` is a `NearSourceSaturationParameters` type
- `anelastic` is an `AnelasticAttenuationParameters` type

The base constructor is: `PathParameters(geo::G, sat::S, ane::A) where {G<:GeometricSpreadingParameters, S<:NearSourceSaturationParameters, A<:AnelasticAttenuationParameters}`

See also: [`FourierParameters`](@ref)
"""
struct PathParameters{G<:GeometricSpreadingParameters,S<:NearSourceSaturationParameters,A<:AnelasticAttenuationParameters}
    geometric::G
    saturation::S
    anelastic::A
end

PathParameters(geometric::GeometricSpreadingParameters, anelastic::AnelasticAttenuationParameters) = PathParameters(geometric, NearSourceSaturationParameters(:None), anelastic)

"""
	get_parametric_type(path::PathParameters)

Extract type most elaborate type from parametric `PathParameters` struct.
This requires dropping down to lower level structs within `path`.
"""
function get_parametric_type(path::PathParameters)
    T = get_parametric_type(path.geometric)
    U = get_parametric_type(path.saturation)
    V = get_parametric_type(path.anelastic)
    if T <: Float64
        if U <: Float64
            return V
        else
            return U
        end
    else
        return T
    end
end

@doc raw"""
	SiteParameters

Custom type defining the site parameters of a Fourier spectrum

- `κ0::T where T<:Real` is the site kappa in units of s
- `ζ0::T where T<:Real` is the Haendel et al. (2020) ζ parameter (for a reference frequency of ``f_0=1`` Hz)
- `η::T where T<:Real` is the Haendel et al. (2020) η parameter
- `model::Symbol` is a symbol identifying the impedance function

The argument `model` is currently one of:
- `:Unit` for a generic unit amplification
- `:Boore2016` for the Boore (2016) amplification for ``V_{S,30}=760`` m/s
- `:AlAtik2021_ask14_620` for the Al Atik & Abrahamson (2021) inversion of ASK14 for ``V_{S,30}=620`` m/s
- `:AlAtik2021_ask14_760` for the Al Atik & Abrahamson (2021) inversion of ASK14 for ``V_{S,30}=760`` m/s
- `:AlAtik2021_ask14_1100` for the Al Atik & Abrahamson (2021) inversion of ASK14 for ``V_{S,30}=1100`` m/s
- `:AlAtik2021_bssa14_620` for the Al Atik & Abrahamson (2021) inversion of BSSA14 for ``V_{S,30}=620`` m/s
- `:AlAtik2021_bssa14_760` for the Al Atik & Abrahamson (2021) inversion of BSSA14 for ``V_{S,30}=760`` m/s
- `:AlAtik2021_bssa14_1100` for the Al Atik & Abrahamson (2021) inversion of BSSA14 for ``V_{S,30}=1100`` m/s
- `:AlAtik2021_cb14_620` for the Al Atik & Abrahamson (2021) inversion of CB14 for ``V_{S,30}=620`` m/s
- `:AlAtik2021_cb14_760` for the Al Atik & Abrahamson (2021) inversion of CB14 for ``V_{S,30}=760`` m/s
- `:AlAtik2021_cb14_1100` for the Al Atik & Abrahamson (2021) inversion of CB14 for ``V_{S,30}=1100`` m/s
- `:AlAtik2021_cy14_620` for the Al Atik & Abrahamson (2021) inversion of CY14 for ``V_{S,30}=620`` m/s
- `:AlAtik2021_cy14_760` for the Al Atik & Abrahamson (2021) inversion of CY14 for ``V_{S,30}=760`` m/s
- `:AlAtik2021_cy14_1100` for the Al Atik & Abrahamson (2021) inversion of CY14 for ``V_{S,30}=1100`` m/s

See also: [`FourierParameters`](@ref), [`site_amplification`](@ref)
"""
struct SiteParameters{T<:Real,U<:Real,V<:Real}
    # site parameters
    κ0::T           # site kappa
    ζ0::U           # zeta parameter for f0=1 Hz
    η::V            # eta parameter
    model::Symbol   # site amplification model
end

SiteParameters(κ0::T) where {T} = SiteParameters(κ0, NaN, NaN, :AlAtik2021_cy14_760)
SiteParameters(κ0::T, model::Symbol) where {T} = SiteParameters(κ0, NaN, NaN, model)
SiteParameters(ζ0::T, η::U) where {T<:Real,U<:Real} = SiteParameters(NaN, ζ0, η, :AlAtik2021_cy14_760)
SiteParameters(ζ0::T, η::U, model::Symbol) where {T<:Real,U<:Real} = SiteParameters(NaN, ζ0, η, model)

"""
	get_parametric_type(site::SiteParameters{T}) where {T} = T

Extract type of `T` from parametric `SiteParameters` struct
"""
get_parametric_type(site::SiteParameters{T}) where {T} = T

"""
    FourierParameters

Custom type for the parameters the Fourier amplitude spectrum.
This type is comprised of source, path and site types, and so has a base constructor of: `FourierParameters(src::S, path::T, site::U) where {S<:SourceParameters, T<:PathParameters, U<:SiteParameters}`

See also: [`SourceParameters`](@ref), [`PathParameters`](@ref), [`SiteParameters`](@ref)
"""
struct FourierParameters{S<:SourceParameters,T<:PathParameters,U<:SiteParameters}
    source::S   # source parameters
    path::T     # path parameters
    site::U     # site parameters
end

# initialiser focussing upon source and path response only. Uses zero kappa and unit amplification
FourierParameters(src::S, path::T) where {S<:SourceParameters,T<:PathParameters} = FourierParameters(src, path, SiteParameters(0.0, :Unit))

"""
	get_parametric_type(fas::FourierParameters)

Extract type most elaborate type from parametric `FourierParameters` struct.
This requires dropping down to lower level structs within `fas`.
"""
function get_parametric_type(fas::FourierParameters)
    T = get_parametric_type(fas.source)
    U = get_parametric_type(fas.path)
    V = get_parametric_type(fas.site)
    if T <: Float64
        if U <: Float64
            return V
        else
            return U
        end
    else
        return T
    end
end


text:
## Test the Minimum Divergence Problem
## min D
## subject to
## A'p < k
using ProgressMeter
using ModelsGenerators
srand(1)

y, x, z = ModelsGenerators.sim_iv_d01(CP = 50);
y = vec(y)

data = Dict{Symbol, Any}(
  :y => y,
  :x => x,
  :z => z
)

size(x)
G = Array(Float64, 100, 5)

f!(G, theta) = broadcast!(*, G, data[:z], data[:y]-data[:x]*theta)

m = MinimumDivergenceProblem(G, zeros(5))

function loglik(theta::DenseVector, m::MinimumDivergenceProblem, f!::Function)
    ## f! is a mutating function
    f!(m.e.mm.S, theta)
    solve!(m)
    sum(log(m.m.inner.x)) + normlogpdf(0, 10, theta[1])
end

_loglik(theta::DenseVector) = loglik(theta, m, f!)

## MCMC simulation
n = 5000
burnin = 1000
sim = Chains(n, 1, names = ["Θ"])
beta = AMWGVariate([0.0])

p = Progress(n, 1)
for i in 1:n
  amwg!(beta, [.5], _loglik, batchsize = 10, adapt = (i <= burnin))
  sim[i, :, 1] = [beta;]
  next!(p)
end
describe(sim)
p = plot(sim)
draw(p, filename="summaryplot.svg")

#####

    srand(1)

    y, x, z = ModelsGenerators.sim_iv_d01(CP = 50);
    y = vec(y)

    data = Dict{Symbol, Any}(
    :y => y,
    :z => z,
    :x => x
    )

    function f(beta)
        delta = beta[1:5]
        theta = beta[6]
        G = broadcast(*, data[:z], data[:y]-data[:z]*delta.*theta)
        H = broadcast(*, data[:z], data[:x]-data[:z]*delta)
        [G H]
    end

    G = f([ones(5); 0])
    m = MinimumDivergenceProblem(G, zeros(10))

    function loglik(theta::DenseVector, m::MinimumDivergenceProblem, f!::Function)
        ## f! is a mutating function
        m.e.mm.S[:] = f(theta)
        solve!(m)
        if m.m.inner.status == 0
            sum(log(m.m.inner.x)) + norm(theta[1:5])*sqrt((1+theta[6]^2))
        else
            -maxintfloat()
        end
    end

    _loglik_beta(beta::DenseVector) = loglik([delta; beta], m, f!)
    _loglik_delta(delta::DenseVector) = loglik([delta; beta], m, f!)

    V           = vcov(estimate!(MDEstimator(f, [zeros(5); .1])))
    Sigma_delta = cholfact(V[1:5,1:5])
    Sigma_beta = [sqrt(V[6,6])]

    n = 100000
    burnin = 2000
    sim = Chains(n, 6, names = ["γ₁", "γ₂", "γ₃", "γ₄", "γ₅", "Θ"])
    beta  = AMWGVariate([0.01])
    delta = AMMVariate(ones(5)*0.3)
    p = Progress(n, 1)
    for i in 1:n
        amwg!(beta, Sigma_beta, _loglik_beta, batchsize = 50, adapt = (i <= burnin))
        amm!(delta, Sigma_delta, _loglik_delta, adapt = (i <= burnin))
        sim[i, :, 1] = [delta; beta]
        next!(p)
    end
    describe(sim)

sim2 = Chains(80001, 2, names = ["||δ||", "Θ"]);
sim2[:,1,1] = mapslices(u -> norm(u), sim.value[20000:100000,1:5,1], 2)
sim2[:,2,1] = sim.value[20000:100000,6,1]

p = plot(sim)
draw(p, filename="summaryplot.svg")

p = contourplot(sim2, bins = 40);
draw(p, filename="contourplot.svg")

function contourplot(c::AbstractChains; bins::Integer=100, na...)
  nrows, nvars, nchains = size(c.value)
  plots = Plot[]
  offset = 1e4 * eps()
  n = nrows * nchains
  for i in 1:(nvars - 1)
    X = c.value[:, i, :]
    qx = linspace(minimum(X) - offset, maximum(X) + offset, bins + 1)
    mx = map(k -> mean([qx[k], qx[k + 1]]), 1:bins)
    idx = Int[findfirst(k -> qx[k] <= x < qx[k + 1], 1:bins) for x in X]
    for j in (i + 1):nvars
      Y = c.value[:, j, :]
      qy = linspace(minimum(Y) - offset, maximum(Y) + offset, bins + 1)
      my = map(k -> mean([qy[k], qy[k + 1]]), 1:bins)
      idy = Int[findfirst(k -> qy[k] <= y < qy[k + 1], 1:bins) for y in Y]
      density = zeros(bins, bins)
      for k in 1:n
        density[idx[k], idy[k]] += 1.0 / n
      end
      p = plot(x=mx, y=my, z=density, Geom.contour,
               Guide.colorkey("Density"),
               Guide.xlabel(c.names[i], orientation=:horizontal),
               Guide.ylabel(c.names[j], orientation=:vertical))
      push!(plots, p)
    end
  end
  return plots
end


text:
using SplitApplyPlot
using Test

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


text:
# Equations of state

using WaterData
import JLD, Dierckx, VoronoiDelaunay, GeometricalPredicates
import WaterData: extracteos
using IterTools: chain

# The bulk of the EOSes and heat capacities are provided in the WaterData
# package. We define one more EOS which is mass-piecewise (for modelling layered
# planets)
""" Equation of state which is piecewise in the mass coordinate

    * `eoses`: Vector of EOSes for each piece.
    * `edges`: Vector defining the edge of each piece. Evaluating
      past the ends of this vector will just evaluate the nearest EOS. """
struct MassPiecewiseEOS{E<:EOS, M<:Mass} <: WaterData.PiecewiseEOS
    eoses::Vector{E}
    edges::Vector{M}
end
function MassPiecewiseEOS(eoses, M::Mass, mass_fractions)
    layer_edges = vcat([-Inf], cumsum(mass_fractions)) .* M
    layer_edges[end] = Inf
    MassPiecewiseEOS(eoses, layer_edges)
end

# EOS evaluation
# (we loop through a bunch of EOSes to avoid ambiguity clashes)

# Split piecewise EOSes appropriately
function extracteos(eos::MassPiecewiseEOS, vs::ValueSet)
    extracteos(eos, mass(vs))
end
function extracteos(eos::WaterData.PressurePiecewiseEOS, vs::ValueSet)
    extracteos(eos, pressure(vs))
end
for T in (MassRadiusPressure, PhysicalValues)
    (eos::MassPiecewiseEOS)(vs::T) = extracteos(eos, vs)(vs)
end

# call definitions have gone to atmosphere.jl as a workaround


text:
# see: https://www.juliabloggers.com/activating-project-environment-in-julia-repl-automatically/?utm_source=ReviveOldPost&utm_medium=social&utm_campaign=ReviveOldPost
using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end
push!( LOAD_PATH, "src/")

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using JuLIP, BenchmarkTools, Base.Threads

# basic EAM potential
data = joinpath(dirname(@__FILE__), "..", "data") * "/"
eam_W = JuLIP.Potentials.FinnisSinclair(data*"W-pair-Wang-2014.plt", data*"W-e-dens-Wang-2014.plt")

@show nthreads()

for L in [5, 10, 20]
   let at = bulk(:W, cubic=true) * L
      @show length(at)
      print("      energy : ")
      @btime energy($eam_W, $at)
      print("      forces : ")
      @btime forces($eam_W, $at)
      print("  energy_map : ")
      @btime energy_map($eam_W, $at)
      print("  forces_map : ")
      @btime forces_map($eam_W, $at)
   end
end


text:
function collect_node_data!(evolver::AbstractMarketModelEvolver,
                            product::MarketModelMultiProduct,
                            dataProvider::MarketModelBasisSystem,
                            rebate::MarketModelExerciseValue,
                            control::MarketModelExerciseValue,
                            numberOfPaths::Int,
                            collectedData::Vector{Vector{NodeData}})

  numerairesHeld = Vector{Float64}()
  number_of_products(product) == 1 || error("single product required")

  numberCashFlowsThisStep = Vector{Int}(undef, 1)
  cashFlowsGenerated = Vector{Vector{MarketModelCashFlow}}(undef, 1)
  #cashFlowsGenerated[1] = Vector{MarketModelCashFlow}(max_number_of_cashflows_per_product_per_step(product))
  cashFlowsGenerated[1] = MarketModelCashFlow[MarketModelCashFlow() for i = 1:max_number_of_cashflows_per_product_per_step(product)]

  rateTimes = get_evolution(product).rateTimes

  cashFlowTimes = possible_cash_flow_times(product)
  rebateTimes = possible_cash_flow_times(rebate)
  controlTimes = possible_cash_flow_times(control)

  n = length(cashFlowTimes)
  productDiscounters = Vector{MarketModelDiscounter}(undef, n)
  @simd for i in eachindex(productDiscounters)
    @inbounds productDiscounters[i] = MarketModelDiscounter(cashFlowTimes[i], rateTimes)
  end

  n = length(rebateTimes)
  rebateDiscounters = Vector{MarketModelDiscounter}(undef, n)
  @simd for i in eachindex(rebateDiscounters)
    @inbounds rebateDiscounters[i] = MarketModelDiscounter(rebateTimes[i], rateTimes)
  end

  n = length(controlTimes)
  controlDiscounters = Vector{MarketModelDiscounter}(undef, n)
  @simd for i in eachindex(controlDiscounters)
    @inbounds controlDiscounters[i] = MarketModelDiscounter(controlTimes[i], rateTimes)
  end

  evolution = get_evolution(product)
  numeraires = evolver.numeraires

  evolutionTimes = evolution.evolutionTimes

  isProductTime = is_in_subset(evolutionTimes, get_evolution(product).evolutionTimes)
  isRebateTime = is_in_subset(evolutionTimes, rebate.evolution.evolutionTimes)
  isControlTime = is_in_subset(evolutionTimes, control.evolution.evolutionTimes)
  isBasisTime = is_in_subset(evolutionTimes, dataProvider.evolution.evolutionTimes)
  isExerciseTime = falses(length(evolutionTimes))
  v = rebate.isExerciseTime
  exercises = idx = 1

  @inbounds @simd for i in eachindex(evolutionTimes)
    if isRebateTime[i]
      if v[idx]
        isExerciseTime[i] = true
        exercises += 1
      end
      idx += 1
    end
  end

  # collectedData = Vector{Vector{NodeData}}(exercises+1)
  resize!(collectedData, exercises)
  for i in eachindex(collectedData)
    # collectedData[i] = Vector{NodeData}(numberOfPaths)
    collectedData[i] = NodeData[NodeData() for i = 1:numberOfPaths]
  end

  for i = 1:numberOfPaths
    start_new_path!(evolver)
    reset!(product)
    reset!(rebate)
    reset!(control)
    reset!(dataProvider)

    principalInNumerairePortfolio = 1.0
    isDone = false
    nextExercise = 1
    collectedData[1][i].cumulatedCashFlows = 0.0

    while true
      currentStep = evolver.currentStep
      advance_step!(evolver)
      currentState = current_state(evolver)
      numeraire = numeraires[currentStep]

      if isRebateTime[currentStep]
        next_step!(rebate, currentState)
      end

      if isControlTime[currentStep]
        next_step!(control, currentState)
      end

      if isBasisTime[currentStep]
        next_step!(dataProvider, currentState)
      end

      if isExerciseTime[currentStep]
        data = collectedData[nextExercise+1][i]

        exerciseValue = get_value(rebate, currentState)
        data.exerciseValue = exerciseValue.amount * numeraire_bonds(rebateDiscounters[exerciseValue.timeIndex], currentState, numeraire) /
                            principalInNumerairePortfolio
        set_values!(dataProvider, currentState, data.values)

        controlValue = get_value(control, currentState)
        data.controlValue = controlValue.amount * numeraire_bonds(controlDiscounters[controlValue.timeIndex], currentState, numeraire) /
                            principalInNumerairePortfolio

        data.cumulatedCashFlows = 0.0
        data.isValid = true
        nextExercise += 1
      end
      if isProductTime[currentStep]
        isDone = next_time_step!(product, currentState, numberCashFlowsThisStep, cashFlowsGenerated)

        for j = 1:numberCashFlowsThisStep[1]
          cf = cashFlowsGenerated[1][j]
          collectedData[nextExercise][i].cumulatedCashFlows += cf.amount * numeraire_bonds(productDiscounters[cf.timeIndex], currentState, numeraire) / principalInNumerairePortfolio
        end
      end

      if ~isDone
        nextNumeraire = numeraires[currentStep + 1]
        principalInNumerairePortfolio *= discount_ratio(currentState, numeraire, nextNumeraire)
      else
        break
      end
    end

    # fill the remaining (un)collected data with nulls
    for j = nextExercise:exercises-1
      data = collectedData[j+1][i]
      data.exerciseValue = data.controlValue = 0.0
      data.cumulatedCashFlows = 0.0
      data.isValid = false
    end
  end
  evolver, collectedData
end


text:
export arrangeRemoteCallDataIntoLocalData,getSourcesIndicesOfKthWorker

function arrangeRemoteCallDataIntoLocalData(DremoteCall::Union{Array{Array{Float64,2},1},Array{Array{ComplexF64,2},1}})
N = length(DremoteCall);
l = 0;
for k=1:N
	l+=size(DremoteCall[k],2);
end
Dlocal = zeros(eltype(DremoteCall[1]),size(DremoteCall[1],1),l)
for k=1:N
	Dlocal[:,k:N:end] = DremoteCall[k];
end
return Dlocal;
end

function getSourcesIndicesOfKthWorker(nworkers::Int64,k::Int64,nsrc::Int64)
return k:nworkers:nsrc;
end


text:
function git_clone(tmp_dir::AbstractString,
                   previous_directory::AbstractString,
                   url::AbstractString,
                   name::AbstractString)
    cd(tmp_dir) do
        run(`git clone $(url) $(name)`)
    end
    return nothing
end

function download_or_clone(tmp_dir, previous_director, reg_url, registry_path, url, name)
    try
        Pkg.Types.download_verify_unpack(reg_url, nothing, registry_path, ignore_existence = true)
    catch err
        # in case of fail, fallback to clone
        git_clone(tmp_dir, previous_director, url, name)
    end
    return nothing
end

const _pkg_server_registry_url = Base.VERSION >= v"1.7.0-" ? Pkg.Registry.pkg_server_registry_url : Pkg.Types.pkg_server_registry_url

function _get_registry(;
                       use_pkg_server,
                       uuid,
                       registry_urls,
                       tmp_dir,
                       name,
                       previous_directory,
                       url)
    if use_pkg_server
        reg_url, registry_urls = _pkg_server_registry_url(uuid, registry_urls)
        registry_path = joinpath(tmp_dir, name)
        if reg_url !== nothing
            download_or_clone(tmp_dir, previous_directory, reg_url, registry_path, url, name)
        else
            git_clone(tmp_dir, previous_directory, url, name)
        end
    else
        git_clone(tmp_dir, previous_directory, url, name)
    end
    return registry_urls
end

function get_latest_version_from_registries!(dep_to_latest_version::Dict{Package, Union{VersionNumber, Nothing}},
                                             registry_list::Vector{Pkg.RegistrySpec};
                                             use_pkg_server::Bool=false)
    original_directory = pwd()
    num_registries = length(registry_list)
    registry_temp_dirs = Vector{String}(undef, num_registries)
    registry_urls = nothing
    for (i, registry) in enumerate(registry_list)
        tmp_dir = mktempdir(; cleanup = true)
        registry_temp_dirs[i] = tmp_dir
        name = registry.name
        url = registry.url
        uuid = registry.uuid
        previous_directory = pwd()
        registry_urls = _get_registry(;
            use_pkg_server,
            uuid,
            registry_urls,
            tmp_dir,
            name,
            previous_directory,
            url,
        )
    end
    for (registry_temp_dir, registry) in zip(registry_temp_dirs, registry_list)
        previous_directory = pwd()
        name = registry.name
        registry_path = joinpath(registry_temp_dir, name)
        cd(registry_path)
        registry_parsed = TOML.parsefile(joinpath(registry_path, "Registry.toml"))
        packages = registry_parsed["packages"]
        for p in packages
            name = p[2]["name"]
            uuid = UUIDs.UUID(p[1])
            package = Package(name, uuid)
            path = p[2]["path"]
            if package in keys(dep_to_latest_version)
                versions = VersionNumber.(collect(keys(TOML.parsefile(joinpath(registry_path, path, "Versions.toml")))))
                old_value = dep_to_latest_version[package]
                if isnothing(old_value)
                    dep_to_latest_version[package] = maximum(versions)
                else
                    dep_to_latest_version[package] = max(old_value, maximum(versions))
                end
            end
        end
        cd(previous_directory)
    end
    cd(original_directory)
    for tmp_dir in registry_temp_dirs
        rm(tmp_dir; force = true, recursive = true)
    end
    return dep_to_latest_version
end


text:
using Diffusion
import MPI

MPI.Init()
nx = ny   = Cint(8)
timesteps = Cint(100)
T         = zeros(nx, ny)
T_ptr     = pointer(T)
comm      = MPI.COMM_WORLD.val

Diffusion.julia_diffusion(T_ptr, timesteps, nx, ny, comm)

MPI.Finalize()

max_stars_repo_name: JuliaLabs/CallingJuliaFromCWithMPI

max_issues_repo_name: JuliaLabs/CallingJuliaFromCWithMPI

max_forks_repo_name: JuliaLabs/CallingJuliaFromCWithMPI


text:
include("testfunctions.jl")
using DiferentialEvolution
using Test

@testset "DiferentialEvolution.jl" begin

    @testset "Classic Diferential Evolution" begin

        @testset "2-Dimensioanl tests" begin
            AckleyT = DE_classic(Ackley,1000,300, [-32.768,-32.768], [32.768,32.768],0.9,0.5)
            SchwefelT = DE_classic(Schwefel,1000,300, [-500,-500], [500,500],0.9,0.5)
            GriewankT = DE_classic(Griewank,1000,300, [-600,-600], [600,600],0.9,0.5)
            RosenbrockT = DE_classic(Rosenbrock,1000,300, [-5,-5], [10,10],0.9,0.5)
            MichalewiczT = DE_classic(Michalewicz,1000,300, [0,0], [pi,pi],0.9,0.5)

            #Ackley     critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(AckleyT[1],0, atol=1e-8)
            #Schwefel   critic value = -418.9829 d   ///  parameters = [420.9687,420.9687,420.9687,420.9687,420.9687]
            D = 2
            @test isapprox(SchwefelT[1],-418.9829*D, atol=1e-4)
            #Griewank   critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(GriewankT[1],0, atol=1e-8)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(RosenbrockT[1],0, atol=1e-8)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(MichalewiczT[1],-1.8013, atol=1e-4)
        end

        @testset "5-Dimensional tests" begin
            D = 5
            AckleyT = DE_classic(Ackley,1000,300, [-32.768 for i in 1:D], [32.768 for i in 1:D],0.9,0.5)
            SchwefelT = DE_classic(Schwefel,1000,300, [-500 for i in 1:D], [500 for i in 1:D],0.9,0.5)
            GriewankT = DE_classic(Griewank,1000,300, [-600 for i in 1:D], [600 for i in 1:D],0.9,0.5)
            RosenbrockT = DE_classic(Rosenbrock,1500,300, [-5 for i in 1:D], [10 for i in 1:D],0.9,0.5)
            MichalewiczT = DE_classic(Michalewicz,1000,300, [0 for i in 1:D], [pi for i in 1:D],0.9,0.5)

            #Ackley     critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(AckleyT[1],0, atol=1e-8)
            #Schwefel   critic value = -418.9829 d   ///  parameters = [420.9687,420.9687,420.9687,420.9687,420.9687]
            @test isapprox(SchwefelT[1],-418.9829*D, atol=1e-4)
            #Griewank   critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(GriewankT[1],0, atol=1e-8)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(RosenbrockT[1],0, atol=1e-8)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(MichalewiczT[1],-4.687658, atol=1e-6)
        end

        @testset "15 - Dimensional tests" begin
            D = 15
            AckleyT = DE_classic(Ackley,3000,700, [-32.768 for i in 1:D], [32.768 for i in 1:D],0.9,0.5)
            SchwefelT = DE_classic(Schwefel,3000,700, [-500 for i in 1:D], [500 for i in 1:D],0.9,0.5)
            GriewankT = DE_dither(Griewank,3000,700, [-600 for i in 1:D], [600 for i in 1:D],0.9,0.5)
            RosenbrockT = DE_dither(Rosenbrock,3000,700, [-5 for i in 1:D], [10 for i in 1:D],0.9,0.5)
            #MichalewiczT = DE_classic(Michalewicz,1000,300, [0 for i in 1:D], [pi for i in 1:D],0.9,0.5)

            #Ackley     critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(AckleyT[1],0, atol=1e-2)
            #Schwefel   critic value = -418.9829 d   ///  parameters = [420.9687,420.9687,420.9687,420.9687,420.9687]
            @test isapprox(SchwefelT[1],-418.9829*D, atol=1e-2)
            #Griewank   critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(GriewankT[1],0, atol=1e-2)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(RosenbrockT[1],0, atol=1e-2)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            #@test isapprox(MichalewiczT[1],-12.952351, atol=1e-6)
        end
#=
        @testset "30 - Dimensional tests" begin
            D = 30
            AckleyT = DE_classic(Ackley,5000,700, [-32.768 for i in 1:D], [32.768 for i in 1:D],0.9,0.5)
            SchwefelT = DE_classic(Schwefel,5000,700, [-500 for i in 1:D], [500 for i in 1:D],0.9,0.5)
            GriewankT = DE_classic(Griewank,5000,700, [-600 for i in 1:D], [600 for i in 1:D],0.9,0.5)
            RosenbrockT = DE_classic(Rosenbrock,5000,700, [-5 for i in 1:D], [10 for i in 1:D],0.9,0.5)
            #MichalewiczT = DE_classic(Michalewicz,1000,300, [0 for i in 1:D], [pi for i in 1:D],0.9,0.5)

            #Ackley     critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(AckleyT[1],0, atol=1e-2)
            #Schwefel   critic value = -418.9829 d   ///  parameters = [420.9687,420.9687,420.9687,420.9687,420.9687]
            @test isapprox(SchwefelT[1],-418.9829*D, atol=1e-2)
            #Griewank   critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(GriewankT[1],0, atol=1e-2)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(RosenbrockT[1],0, atol=1e-2)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            #@test isapprox(MichalewiczT[1],-4.687658, atol=1e-6)
        end
    =#
    end

    @testset "Dither Diferential Evolution" begin

    end

    @testset "Jither Diferential Evolution" begin

    end

end


text:
# COUNTS LINE NUMBERS IN LIBRARY AND UPDATES THEM IN THE MAIN FILE

# list of all line numbers
line_numbers = Int64[]

# open the main file
main_file_filename = "/home/jattig/PackageDevelopement/LatticePhysics.jl/src/LatticePhysics.jl"

# find the folder
main_file_folder = main_file_filename[1:findlast(main_file_filename, '/')]

# get all files within that folder
all_files = readdir(main_file_folder)

# open the main file
main_file = open(main_file_filename, "r")
# read the complete main file
main_file_lines = String[rstrip(line) for line in readlines(main_file)]
# close the main file
close(main_file)

# include lines
include_lines = Int64[]

# process all lines
for l in 5:length(main_file_lines)
    # look at line lines preceeding line l to find include lines
    if  main_file_lines[l-3] == "################################################################################" &&
        main_file_lines[l-2] == "" &&
        length(main_file_lines[l-1]) > 10 &&
        main_file_lines[l-1][1:10] == "# included"
        # this line contains the include
        push!(include_lines, l)
    end
end

# extract the included files
include_files = String[]
for line in [main_file_lines[i] for i in include_lines]
    # extract the filename
    filename = line[findfirst(line, '"')+1:findlast(line, '"')-1]
    # check if it is in the directory
    if !(filename in all_files)
        println("file not found: $(filename)")
        push!(include_files, "NOTHING.jl")
    else
        # push the filename into the list
        push!(include_files, filename)
    end
end

# extract line number information
for include_filename in include_files
    # get the number of lines in the corresponding file
    f = open("$(main_file_folder)$(include_filename)", "r")
    push!(line_numbers, length(readlines(f)))
    close(f)
    # print the success
    println("$(include_filename)  ($(line_numbers[end]) lines)")
end

# overwrite the lines in the main file
for i in 1:length(include_lines)
    # get the line
    l = include_lines[i]
    count = line_numbers[i]
    # compose the new line
    line = "# included in subfile ($(count) lines)"
    # substitute the line
    main_file_lines[l-1] = line
end

# overwrite the total count string
main_file_lines[end-2] = "# total lines: $(length(main_file_lines))"
for l in line_numbers
    main_file_lines[end-2] = "$(main_file_lines[end-2]) + $(l)"
end
# add the total number
main_file_lines[end-1] = "# = $(sum(line_numbers) + length(main_file_lines)) lines"
# print the total number
println("  --> total: $(sum(line_numbers) + length(main_file_lines)) lines")

# write to file
open(main_file_filename, "w") do main_file
    for line in main_file_lines
        write(main_file, "$(line)\n")
    end
end


text:
"""
 recognize

 """

"""
Functions for recognizing properties of the a LPModel

"""

export recognize

function recognize(::LPModel, ::AbstractRecognizeAttribute) <: Bool end

function recognize(lpmodel::LPModel, ::DifferenceConstraints)
    A, b = get_constraint_matrices(lpmodel)
    variable_count = size(A)[2]
    if (variable_count <= 1 || size(A)[1]==0)
        return false
    end
    compare = zeros(variable_count)
    compare[1] = 1
    compare[2] = -1
    for row in eachrow(A)
        if !issetequal(row, compare)
            return false
        end
    end
    return true
end

function recognize(lpmodel::LPModel, ::AllIntegerVariables)
    integers = get_integer_variables(lpmodel)
    # Filter out the NaNs:
    num_imntegers = length(filter(integers->!isnan(integers),integers))
    return num_imntegers == lpmodel.num_variables_created
end

function recognize(lpmodel::LPModel, ::AllIntegerConstraintBounds)
    A, b = get_constraint_matrices(lpmodel)
    # If after truncation the element is equal to itself, it is an integer:
    for element in b
        if !(trunc(Int, element) == element)
            return false
        end
    end
    return true
end

"""

function recognize(lpmodel::LPModel, ::ConstantObjective)
    c = get_objective_vector(lpmodel)
    return (c == zeros(size(c)[1]))
end"""

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
using Reexport

@reexport module Resampling
    import Interpolations
    import Interpolations: Linear, BoundaryCondition
    import ..bin

    import ..UVAL_COLLECTION_TYPES
    import ..UncertainDatasets:
        AbstractUncertainValueDataset,
        UncertainIndexValueDataset,
        AbstractUncertainDataset
    import ..UncertainValues:
        UncertainValue,
        AbstractUncertainValue

    function resample end

    ###################################
    # Resampling schemes
    ###################################
    include("resampling_schemes/AbstractUncertainDataResampling.jl")
    include("resampling_schemes/BinnedResamplings.jl")
    include("resampling_schemes/ConstrainedResampling.jl")
    include("resampling_schemes/ConstrainedValueResampling.jl")
    include("resampling_schemes/ConstrainedIndexValueResampling.jl")
    include("resampling_schemes/SequentialResampling.jl")
    include("resampling_schemes/SequentialInterpolatedResampling.jl")
    include("resampling_schemes/binned_resamplings.jl")
    include("resampling_schemes/InterpolateAndBin.jl")
    include("resampling_schemes/RandomSequences.jl")

    # Extend some methods to allow easier resampling.
    include("resampling_with_schemes/constrain_with_schemes.jl")

    ###################################
    # Resampling uncertain values
    ###################################
    # Uncertain values based on distributions
    include("uncertain_values/resample_uncertainvalues_distributions.jl")

    # With constraints
    include("uncertain_values/resample_uncertainvalues_theoretical.jl")
    include("uncertain_values/resample_uncertainvalues_theoretical_withconstraints.jl")

    include("uncertain_values/resample_uncertainvalues_kde.jl")
    include("uncertain_values/resample_uncertainvalues_kde_withconstraints.jl")

    include("uncertain_values/resample_certainvalues.jl")
    include("uncertain_values/resample_uncertainvalues_populations.jl")

    include("uncertain_values/resample_measurements.jl")

    #########################################
    # Resampling tuples of uncertain values
    #########################################
    include("uncertain_tuples/uncertain_tuples.jl")

    #########################################
    # Resampling vectors of uncertain values
    #########################################
    include("uncertain_values/resampling_vector_uncertainvalues.jl")

    ###################################
    # Resampling uncertain datasets
    ###################################

    # Element-wise resampling for all subtypes of AbstractUncertainValueDataset
    include("uncertain_dataset/resample_abstractuncertainvaluedataset_elwise.jl")

    # Specialized resampling for each type of dataset.
    include("uncertain_dataset/resample_uncertaindataset_index.jl")
    include("uncertain_dataset/resample_uncertaindataset_value.jl")
    include("uncertain_dataset/resample_uncertaindataset_indexvalue.jl")

    # Resampling vectors of uncertain values
    include("uncertain_vectors/resample_uncertain_vectors.jl")

    #########################################
    # Ordered resampling
    #########################################
    include("ordered_resampling/resample_sequential.jl")
    include("ordered_resampling/strictlyincreasing.jl")
    include("ordered_resampling/strictlydecreasing.jl")

    #########################################
    # Resampling with interpolation
    #########################################
    include("resampling_with_interpolation/resample_linear_interpolation.jl")

    ################################
    # Apply function with resampling
    ################################
    include("apply_func.jl")

    ################################
    # Resampling in-place
    ################################
    include("resampling_inplace.jl")

    ################################
    # Resampling with schemes
    ################################
    include("resampling_with_schemes/resampling_schemes_binned.jl")
    include("resampling_with_schemes/resampling_schemes_interpolated_binned.jl")
    include("resampling_with_schemes/resampling_schemes_constrained.jl")
    include("resampling_with_schemes/resampling_schemes_sequential.jl")

    ################################
    # Interpolation
    ################################
    include("binning/bin_BinnedResampling.jl")
    include("binning/bin_BinnedWeightedResampling.jl")

    export resample, resample!, resample_elwise, bin
end # module

"""
	Resampling

A module defining resampling methods for uncertain values defined in the
`UncertainValues` module and uncertain datasets defined in the
`UncertainDatasets` module.
"""
Resampling


text:
module MPSLabXD

# Imports
using TensorLabXD
const TX = TensorLabXD

using TupleLabXD
const TU = TupleLabXD

using TensorContractionsXD
const TC = TensorContractionsXD

using KrylovKit

using LinearAlgebra

# Export
export AbstractMPS, FiniteMPS
export field, spacetype, sectortype, get_chi, entanglement_entropy

export make_right_canonical, make_left_canonical, norm

export MPO, TransverseFieldIsing

export FiniteEnvironment, initialize_FiniteEnvironment, update_left_env, update_right_env

export finiteDMRG, one_sweep, one_update, build_two_site_H_eff
# source codes
include("MPS/abstractmps.jl")
include("MPS/finitemps.jl")
include("MPO/mpo.jl")
include("Environment/finiteEnvironment.jl")
include("GroundState/finiteDMRG.jl")

end


text:
# test.jl

function twice_array(x::Array{T}) where T <: Number
    return 2x
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer