text:
#This file will be turned in more flexible algebra routines

##header
import SpecialFunctions

## Polynomial type
"""
Polynomial type
"""
struct Polynomial #TODO: make type parametric use eltype(f)
 coeffs
end

#make Polynomial callable with Horner-scheme evaluation
function (p::Polynomial)(z,k::Int=0)
    p=derive(p,k)
    f=0.0+0.0im #typing
    for i = length(p.coeffs):-1:1
       f*=z
       f+=p.coeffs[i]
    end
    return f
end

import Base.show
import Base.string

function string(p::Polynomial)
    N=length(p.coeffs)
    if N==0
        txt="0"
    else
        txt="$(p.coeffs[1])"
        for n=2:N
            coeff=string(p.coeffs[n])
            if occursin("+",coeff) || occursin("-",coeff)
                txt*="+($coeff)*z^$(n-1)"
            else
                txt*="+$coeff*z^$(n-1)"
            end
        end
    end
    return txt
end

function show(io::IO,p::Polynomial)
    print(io,string(p))
end

import Base.+
import Base.-
import Base.*
import Base.^

function +(a::Polynomial,b::Polynomial)
    if length(a.coeffs)>length(b.coeffs)
        a,b=b,a
    end
    c=copy(b.coeffs)
    for (idx,val) in enumerate(a.coeffs)
        c[idx]+=val
    end
    return Polynomial(c)
end

function -(a::Polynomial,b::Polynomial)
    return +(a,-b)
end

function *(a::Polynomial,b::Polynomial)
    I=length(a.coeffs)
    J=length(b.coeffs)
    K=I+J-1
    c=zeros(ComplexF64,K) #TODO: sanity check for type of b
    idx=1
    for k =1:K
        for i =1:k
            #TODO: figure out when to break the loop
            j=k-i+1
            if i>I ||j>J
                continue
            end
            c[k]+=a.coeffs[i]*b.coeffs[j]
        end
    end
    return Polynomial(c)
end

#scalar multiplication
function *(a::Polynomial,b::Number)
    return a*Polynomial([b])
end
function *(a::Number,b::Polynomial)
    return Polynomial([a])*b
end

function ^(p::Polynomial, k::Int)
    b=Polynomial(copy(p.coeffs))
    for i=1:k-1
        b*=p
    end
    return b
end

function prodrange(start,stop)
    return SpecialFunctions.factorial(stop)/SpecialFunctions.factorial(start)
end

"""
    b=derive(a::Polynomial,k::Int=1)

Compute `k`th derivative of Polynomial `a`.
"""
function derive(a::Polynomial,k::Int=1)
    N=length(a.coeffs)
    if k>=N
        return Polynomial([]) #TODO: type
    end
    d=zeros(typeof(a.coeffs[1]),N-k)
    for i=1:N-k
        d[i]=a.coeffs[i+k]*prodrange(i-1,i+k-1)
    end
    return Polynomial(d)
end

"""
    g=shift(f,Δz)

Shift Polynomial with respect to its argument. (Taylor shift)
"""
function shift(f::Polynomial,Δz)
    g=Array{eltype(f)}(undef,length(f.coeffs))
    for n=0:length(f.coeffs)-1
        g[n+1]=f(Δz)/SpecialFunctions.factorial(n*1.0)
        f=Polynomial(f.coeffs[2:end])
        f.coeffs.*=1:length(f.coeffs)
    end
    return Polynomial(g)
end

function scale(p::Polynomial,a::Number)
        p=p.coeffs
        s=1
        for idx =1: length(p)
            p[idx]*=s
            s*=a
        end
        return Polynomial(p)
end

#TODO: write macro nto create !-versions

## tests
a=Polynomial([0.0,1.0])
b=a*a
c=derive(a)
d=a+b+c
e=d^3
f=shift(e,2)
e(2)-f(0)
## rational Polynomial
struct rational_Polynomial
    num::Polynomial
    den::Polynomial
end

function derive(a::rational_Polynomial)
    return rational_Polynomial(derive(a.num)*a.den-a.num*derive(a.den),a.den^2)
end

"""
    derive(a::Any, k:Int)

Take the `k`th derivative of `a` by calling the derive function multiple times.

# Notes
This is a fallback implementation.
"""
function derive(a::Any,k::Int)
    for i=1:k
        a=derive(a)
    end
    return a
end
##
function ndd(vararg...)
    #parse input TODO: sanity check that z is unique
    #get number of points
    #(confluent points are counted multiple times)
    N=0
    for i = 1:length(vararg)÷2
        N+=length(vararg[2*i])
    end
    #initialize lists
    z=Array{ComplexF64}(undef,N) #TODO: proper type rules
    coeffs=Array{ComplexF64}(undef,N)
    cnfl=Array{Int64}(undef,N) # confluence counter array
    n=0 #total counter
    for i = 1:length(vararg)÷2
        z[1+n:n+length(vararg[2*i])] .= vararg[2*i-1]
        coeffs[1+n:n+length(vararg[2*i])]  = vararg[2*i]
        cnfl[1+n:n+length(vararg[2*i])] = 0:length(vararg[2*i])-1
        n+=length(vararg[2*i])
    end

    #pyramid scheme for computation of divided-difference scheme
    println("###########")
    println(cnfl)
    Z=copy(z)
    for i=1:N-1
        println(coeffs)
        for j=N:-1:(i+1)
            Δz=z[j]-z[j-i]
            if Δz!=0 # this is a non-confluent point
               coeffs[j]=(coeffs[j]-coeffs[j-1-cnfl[j-1]])/Δz
            end
            if cnfl[j-1]!=0
                cnfl[j-1]-=1
            end
        end
    end
    println(coeffs)

    # construct newton polynomial from divided difference coefficients
    basis=Polynomial([1])
    p=Polynomial([0])
    for i=1:N
        p+=basis*coeffs[i]
        basis*=Polynomial([-z[i],1])
    end

    return p
end

##
m0=Polynomial([1])
m1=Polynomial([-1,1])
m2=Polynomial([-2,1])
m3=Polynomial([-3,1])

p=3*m0+(-1)*m1+2.5*m1*m2
##
p=ndd(1,[3,],2,[2,],3,[6,])
p.([1,2,3])
##
p=ndd(2,[6,3,5])
p(2,3)
##
p=ndd(1,[1,2],2,[4])
p.([1,2],0)
##
vals=(1,[1,2,3],2,[4,5,0,7,8,9],3,[60,1])
p=ndd(vals...)
##
println("#########Here##########")
for i=1:length(vals)÷2
    z=vals[2*i-1]
    for (k,coeff) in enumerate(vals[2*i])
        k=k-1
        test=p(z,k)/factorial(k)-coeff==0
        if !test
            println("z:$z k:$k coeff:$coeff")
        end
    end
end

## TODO next: Find reference for this algorithm.
function newton_divided_difffunction(vararg...)
    #hässliches steigungsschema...
    coeffs=Dict()
    z=Array{ComplexF64}(undef,length(vararg)÷2)
    comb=[]
    #initialize divided difference scheme
    #use taylor coefficient at confluent points
    for (idx,n) in enumerate(1:2:length(vararg))
        z[idx]=vararg[n]
        taylor_coeffs=vararg[n+1]
        loc_tab=repeat([idx],length(taylor_coeffs))
        push!(comb,loc_tab...)
        for (i,j) in enumerate(1:length(loc_tab)) #TODO: enumeration and i are not required
            coeffs[loc_tab[1:j]]=taylor_coeffs[j]
        end
    end

    #fill missing values in the newton divided difference table
    for n =1:length(comb)
        for j =1:length(comb)-n
            idx=comb[j:j+n]
            if idx ∉ keys(coeffs)
                a=coeffs[idx[1:end-1]]
                b=coeffs[idx[2:end]]
                z_a=z[idx[1]]
                z_b=z[idx[end]]
                coeffs[idx]=(b-a)/(z_b-z_a)
            end
        end
    end

    #construct newton_Polynomial
    basis_function=Polynomial([1])
    idx=comb[1:1]

    c=Polynomial(coeffs[idx])
    p=basis_function*c
    for n=2:length(comb)
        idx=comb[1:n]
        c=Polynomial(coeffs[idx])
        basis_function=basis_function*Polynomial(-z[idx[end-1]],1]
        p=p+ᴾc*ᴾbasis_function
    end
    basis_function=basis_function*ᴾ[-z[comb[end]],1]

    #Kronecker's Algorithm to find the interpolating Pade approximant
    #initialize Kronecker's algorithm
    # (see Chapter 7.1 (page 341) in G. A. Baker and P. Graves-Morris.
    # Padé Approximants, 2nd Edition)
    p0=basis_function
    q0=[]
    q=[1]
    pade_table=[]
    pade_table=push!(pade_table,(p,q))

    while length(p)>1
        K=2 #system size

        A=zeros(ComplexF64,K,K)
        y=zeros(ComplexF64,K)
        for k=1:K
            y[k]=p0[end-k+1]
            for i=1:k
                    A[k,i]=p[end-k+i]
            end
        end
        x=A\y
        #TODO: sanity check!
        reverse!(x)
        p,p0=x*ᴾp-ᴾp0,p
        q,q0=x*ᴾq-ᴾq0,q
        p=p[1:end-K]
        push!(pade_table,(p,q))
    end

    return pade_table

end

## helper functions to directly use arrays
function *ᴾ(a,b)
    I=length(a)
    J=length(b)
    K=I+J-1
    c=zeros(ComplexF64,K)
    idx=1
    for k =1:K
        for i =1:k
            #TODO: figure out when to break the loop
            j=k-i+1
            if i>I ||j>J
                continue
            end
            c[k]+=a[i]*b[j]
        end
    end
    return c
end

function +ᴾ(a,b)
    if length(a)>length(b)
        a,b=b,a
    end
    for (idx,val) in enumerate(a)
        b[idx]+=val
    end
    return b
end

function -ᴾ(a,b)
    return +ᴾ(a,-b)
end

function derive(a)
    N=length(a)
    d=zeros(ComplexF64,N-1)
    for i=1:N-1
        d[i]=a[i+1]*i
    end
    return d
end
##
function compute_newton_polynomial(vararg...)
    #hässliches steigungsschema...
    coeffs=Dict()
    z=Array{ComplexF64}(undef,length(vararg)÷2)
    comb=[]
    #initialize divided difference scheme
    #use taylor coefficient at confluent points
    for (idx,n) in enumerate(1:2:length(vararg))
        z[idx]=vararg[n]
        taylor_coeffs=vararg[n+1]
        loc_tab=repeat([idx],length(taylor_coeffs))
        push!(comb,loc_tab...)
        for (i,j) in enumerate(1:length(loc_tab))
            coeffs[loc_tab[1:j]]=taylor_coeffs[j]
        end
    end

    #fill missing values in the newton divided difference table
    for n =1:length(comb)
        for j =1:length(comb)-n
            idx=comb[j:j+n]
            if idx ∉ keys(coeffs)
                a=coeffs[idx[1:end-1]]
                b=coeffs[idx[2:end]]
                z_a=z[idx[1]]
                z_b=z[idx[end]]
                coeffs[idx]=(b-a)/(z_b-z_a)
            end
        end
    end

    #construct newton_polynomial
    basis_function=[1]
    idx=comb[1:1]

    c=[coeffs[idx]]
    p=basis_function*ᴾc
    for n=2:length(comb)
        idx=comb[1:n]
        c=[coeffs[idx]]
        basis_function=basis_function*ᴾ[-z[idx[end-1]],1]
        p=p+ᴾc*ᴾbasis_function
    end
    basis_function=basis_function*ᴾ[-z[comb[end]],1]

    #Kronecker's Algorithm to find the interpolating pade approximant
    #initialize Kronecker's algorithm
    p0=basis_function
    q0=[]
    q=[1]
    pade_table=[]
    pade_table=push!(pade_table,(p,q))

    while length(p)>1
        K=2 #system size

        A=zeros(ComplexF64,K,K)
        y=zeros(ComplexF64,K)
        for k=1:K
            y[k]=p0[end-k+1]
            for i=1:k
                    A[k,i]=p[end-k+i]
            end
        end
        x=A\y
        #TODO: sanity check!
        reverse!(x)
        p,p0=x*ᴾp-ᴾp0,p
        q,q0=x*ᴾq-ᴾq0,q
        println("#####")
        println(p[end-K+1:end])
        p=p[1:end-K]
        push!(pade_table,(p,q))
    end

    return pade_table

end

##
polyval(p,z)= sum(p.*(z.^(0:length(p)-1)))

function taylor_shift(f,Δz)
    g=Array{eltype(f)}(undef,length(f))
    for n=0:length(f)-1
        g[n+1]=polyval(f,Δz)/SpecialFunctions.factorial(n*1.0)
        f=f[2:end]
        f.*=1:length(f)
    end
    return g
end

#Test
f=[-42.168, 2im, π, 1, 0.05]
z= π*1im-5.923779
g=taylor_shift(f,z)
polyval(f,z)-polyval(g,0)

#TODO implement Kronecker's algorithm
function multi_point_pade(L,M,vararg...;Z0=0)
    #TODO sanity checks length vararg should be divisibl #and the sum of the entries should amound to L+M

    #Build linear system

end

function newton_polynomial(vararg...)
    N=0
    for idx =1:2:length(vararg)
        N+=length(vararg[idx+1])
    end
    A=zeros(ComplexF64,N,N)
    y=zeros(ComplexF64,N)
    eqidx=1
    for idx =1:2:length(vararg)
        coeffs=ones(ComplexF64,N)
        z0=vararg[idx]
        z=[z0^i for i in 0:N-1]
        k=1
        fac=1
        for val in vararg[idx+1]
            y[eqidx]=val*fac
            A[eqidx,:]=coeffs.*z
            coeffs[k]=0
            z[k]=0
            for (i,j) in enumerate(k+1:N)
                coeffs[j]*=i
                z[j]=z0^(i-1)
            end
            fac*=k
            k+=1
            eqidx+=1
        end

    end

    x=A\y
    return x

end

function newton_divided_difffunction(vararg...)
    #hässliches steigungsschema...
    coeffs=Dict()
    z=Array{ComplexF64}(undef,length(vararg)÷2)
    comb=[]
    #initialize divided difference scheme
    #use taylor coefficient at confluent points
    for (idx,n) in enumerate(1:2:length(vararg))
        z[idx]=vararg[n]
        taylor_coeffs=vararg[n+1]
        loc_tab=repeat([idx],length(taylor_coeffs))
        push!(comb,loc_tab...)
        for (i,j) in enumerate(1:length(loc_tab))
            coeffs[loc_tab[1:j]]=taylor_coeffs[j]
        end
    end

    #fill missing values in the newton divided diffrence table
    for n =1:length(comb)
        for j =1:length(comb)-n
            idx=comb[j:j+n]
            if idx ∉ keys(coeffs)
                a=coeffs[idx[1:end-1]]
                b=coeffs[idx[2:end]]
                z_a=z[idx[1]]
                z_b=z[idx[end]]
                coeffs[idx]=(b-a)/(z_b-z_a)
            end
        end
    end

    #construct newton_polynomial
    basis_function=[1]
    idx=comb[1:1]

    c=[coeffs[idx]]
    p=basis_function*ᴾc
    for n=2:length(comb)
        idx=comb[1:n]
        c=[coeffs[idx]]
        basis_function=basis_function*ᴾ[-z[idx[end-1]],1]
        p=p+ᴾc*ᴾbasis_function
    end
    basis_function=basis_function*ᴾ[-z[comb[end]],1]

    #Kronecker's Algorithm to find the interpolating pade approximant
    #initialize Kronecker's algorithm
    p0=basis_function
    q0=[]
    q=[1]
    pade_table=[]
    pade_table=push!(pade_table,(p,q))

    while length(p)>1
        K=2 #system size

        A=zeros(ComplexF64,K,K)
        y=zeros(ComplexF64,K)
        for k=1:K
            y[k]=p0[end-k+1]
            for i=1:k
                    A[k,i]=p[end-k+i]
            end
        end
        x=A\y
        #TODO: sanity check!
        reverse!(x)
        p,p0=x*ᴾp-ᴾp0,p
        q,q0=x*ᴾq-ᴾq0,q
        println("#####")
        println(p[end-K+1:end])
        p=p[1:end-K]
        push!(pade_table,(p,q))
    end

    return pade_table

end

pade_tab=newton_divided_difffunction(1,[1],2,[pi,3,-1],0,[1,2,5,1])
p,q=pade_tab[3]
test=polyval(p,2)/polyval(q,2)
##
q=newton_polynomial(1,[1],2,[0,3,-1],0,[1,2,5,1])

#TODO: write testing scheme


text:
#Automatically name series, so that downstream functions like legend which use names have one
function seriesnames!(ec::EChart)

    for (i, x) in enumerate(ec.series)
        x.name == nothing ? x.name = "Series $i" : nothing
    end

    return ec

end

function seriesnames!(ec::EChart, names::AbstractVector{String})

    length(ec.series) != length(names) ? error("Names not equal to number of Series") : nothing

    for i in 1:length(ec.series)
        ec.series[i].name = names[i]
    end

    #Modify legend attribute to use new seriesnames
    ec.legend.data = names

    return ec

end

seriesnames!(ec::EChart, names::AbstractVector) = seriesnames!(ec, String[string(x) for x in names])

ext: jl


text:
@inline function calculate_Φ_lab(K::Int)
    Φ_lab = K > 1 ? Matrix{Int}(undef, binomial(K, 2), 2) : [1 1]
    if K > 1
        i = 1
        for k1 in 1:(K - 1)
            for k2 in (k1 + 1):K
                Φ_lab[i, :] = [k1, k2]
                i += 1
            end
        end
    end
    return Φ_lab
end

@inline function calc_ESS(logweight)
    ESSnum = 0.0
    ESSdenom = 0.0
    max_l = maximum(logweight)
    for l in logweight
        num = exp(l - max_l)
        ESSnum += num
        ESSdenom += num ^ 2
    end
    return (ESSnum ^ 2) / ESSdenom
end

function draw_partstar(logweight, particles)
    u = rand() / particles
    pprob = cumsum(exp.(logweight .- maximum(logweight)))
    partstar = zeros(Int, particles)
    i = 0
    for p = 1:particles
        while pprob[p] / last(pprob) >= u
            u += 1 / particles
            i += 1
            partstar[i] = p
        end
    end
    # Systematic resampling not ideal for CSMC
    # Particle indices are sorted, so can't just set
    # the first as reference trajectory
    # Instead: shuffle, replace first, sort
    shuffle!(partstar)
    partstar[1] = 1
    sort!(partstar)
    return partstar
end

function Φ_upweight!(logweight, sstar, K::Int, Φ, particles)
    Φ_lab = calculate_Φ_lab(K)
    for i in 1:binomial(K, 2)
        Φ_log = log(1 + Φ[i])
        for p in 1:particles
            logweight[p] += (sstar[p, Φ_lab[i, 1]] == sstar[p, Φ_lab[i, 2]]) * Φ_log
        end
    end
    return
end

function align_labels!(s::Array, Φ::Array, γ::Array, N::Int, K::Int)
    K == 1 && return
    Φ_lab = calculate_Φ_lab(K)
    Φ_log = log.(Φ .+ 1)
    # unique_s = unique(s)
    # No need to permute labels for dataset 1
    @inbounds for k = 1:K
        occupied = unique(s[:, k])
        relevant_Φs = Φ_log[(Φ_lab[:, 1] .== k) .| (Φ_lab[:, 2] .== k)]
        for label in occupied
            label_ind  = s[:, k] .== label
            all(label_ind .== false) && continue
            label_rows = s[label_ind, setdiff2(K, k)]
            # Only consider most frequent label in each other dataset
            # as this will always dominate others
            # and saves time
            # However, detailed balance?
            for new_label in 1:N # mapslices(mode, label_rows, dims = 1)
                new_label == label && continue
                new_label_ind   = s[:, k] .== new_label
                new_label_rows  = s[new_label_ind, setdiff2(K, k)]
                log_phi_sum     = sum(count_equals(label_rows, label) .* relevant_Φs + count_equals(new_label_rows, new_label) .* relevant_Φs)
                log_phi_sum_swap = sum(count_equals(label_rows, new_label) .* relevant_Φs + count_equals(new_label_rows, label) .* relevant_Φs)
                accept = exp(log_phi_sum_swap - log_phi_sum)
                if rand() < accept
                    s[label_ind, k]        .= new_label
                    s[new_label_ind, k]    .= label
                    γ[new_label, k], γ[label, k] = γ[label, k], γ[new_label, k]
                    label = new_label
                    label_ind  = s[:, k] .== label
                    label_rows = s[label_ind, setdiff2(K, k)]
                end
            end
        end
    end
end

@inline function count_equals(A::Array, b::Int)
    out = zeros(Float64, size(A, 2))
    for i in 1:size(A, 1)
        for j in 1:size(A, 2)
            if A[i, j] == b
                out[j] += 1.
            end
        end
    end
    return out
end

@inline function setdiff2(K::Int, b::Int)
    # Specifically find diff between set 1:K and
    # integer b
    out = ones(Bool, K)
    out[b] = 0
    return out
end

@inline function find_logical(A::Vector, b::Int)
    # Return a bit-array indicating where A has value b
    out = zeros(Bool, size(A, 1))
    for i in eachindex(A)
        if A[i] == b
            out[i] = true
        end
    end
    return out
end

@inline function findindex(A, b)
    # Find first occurrence of b in A
    for (i,a) in enumerate(A)
        if a == b
            return i
        end
    end
end

@inline function findindices(A, b)
    # Find all occurrences of b in A
    # Specifically for finding occurrences of gamma
    out = Int[]
    for (i, a) in enumerate(A)
        if a == b
            push!(out, i)
        end
    end
    return out
end

function findZindices(k, K, n, N)
    # Instead of searching where the γc_combn = a
    # by searching the vector
    # We can know this in advance
    # This way is >10x faster for larger problems
    out = zeros(Int, N ^ (K - 1))
    start = (n - 1) * N ^ (k - 1) + 1
    ind = 1
    for i in 1:(N ^ (K - k))
        for j in start:(start - 1 + N ^ (k - 1))
            out[ind] = j
            ind += 1
        end
        start += N ^ k
    end
    return out
end

function countn(A, b)
    out = 0
    for a in A
        if a == b
            out += 1
        end
    end
    return out
end

function wipedout0(v1, v2, x)
    # Is the number of occurrences of x greater in
    # v2 than v1?
    if x == 1
        return false
    end
    count1 = 0
    @inbounds for v in v2
        if v == x
            count1 += 1
        end
    end
    for j in 1:size(v1, 2)
    @inbounds for i in 1:size(v1, 1)
            if v1[i, j] == x
                count1 -= 1
                if count1 < 0
                    return false
                else
                    break
                end
            end
        end
    end
    return true
end

function wipedout(v1, v2, x)
    return count(y -> y == x, v2) >= count(y -> y == x, v1)
end

@inline function canonicalise_IDs!(IDs)
    tmp = view(IDs, sortperm(IDs))
    current = - 1
    count = 0
    # tmp[1] = count
    for (i, u) in enumerate(tmp)
        if u != current
            count += 1
            current = u
        end
        tmp[i] = count
    end
end

function draw_uniform_sorted(n)
    @assert n != 0 "WTF"
    u = rand(n)
    u[n] = exp(log(u[n]) / n)
    for i in (n - 1):-1:1
        u[i] = exp(log(u[i]) / i) * u[i + 1]
    end
    return u
end


text:
# """
#     logΣexp, Σexp = logsumexp!(p::WeightedParticles)
# Return log(∑exp(w)). Modifies the weight vector to `w = exp(w-offset)`
# Uses a numerically stable algorithm with offset to control for overflow and `log1p` to control for underflow. `Σexp` is the sum of the weifhts in the state they are left, i.e., `sum(exp.(w).-offset)`.
#
# References:
# https://arxiv.org/pdf/1412.8695.pdf eq 3.8 for p(y)
# https://discourse.julialang.org/t/fast-logsumexp/22827/7?u=baggepinnen for stable logsumexp
# """
# function logsumexp!(p::WeightedParticles)
#     N = length(p)
#     w = p.logweights
#     offset, maxind = findmax(w)
#     w .= exp.(w .- offset)
#     Σ = sum_all_but(w,maxind) # Σ = ∑wₑ-1
#     log1p(Σ) + offset, Σ+1
# end
#
# """
#     sum_all_but(w, i)
#
# Add all elements of vector `w` except for index `i`. The element at index `i` is assumed to have value 1
# """
# function sum_all_but(w,i)
#     w[i] -= 1
#     s = sum(w)
#     w[i] += 1
#     s
# end
#
# """
#     loglik = resample!(p::WeightedParticles)
# Resample the particles based on the `p.logweights`. After a call to this function, weights will be reset to sum to one. Returns log-likelihood.
# """
# function resample!(p::WeightedParticles)
#     N = length(p)
#     w = p.logweights
#     logΣexp,Σ = logsumexp!(p)
#     _resample!(p,Σ)
#     # fill!(p.weights, 1/N)
#     fill!(w, -log(N))
#     logΣexp - log(N)
# end

# """
# In-place systematic resampling of `p`, returns the sum of weights.
# `p.logweights` should be exponentiated before calling this function.
# """
# function _resample!(p::WeightedParticles,Σ)
#     x,w = p.particles, p.logweights
#     N = length(w)
#     bin = w[1]
#     s = rand()*Σ/N
#     bo = 1
#     for i = 1:N
#         @inbounds for b = bo:N
#             if s < bin
#                 x[i] = x[b]
#                 bo = b
#                 break
#             end
#             bin += w[b+1] # should never reach here when b==N
#         end
#         s += Σ/N
#     end
#     Σ
# end

for PT in ParticleSymbols
    @eval begin
        """
            bootstrap([rng::AbstractRNG,] p::Particles, n = nparticles(p))

        Return Particles resampled with replacement. `n` specifies the number of samples to draw. Also works for arrays of Particles, in which case a single set of indices are drawn and used to extract samples from all elements in the array.
        """
        function bootstrap(rng::AbstractRNG, p::$PT, n::Integer = nparticles(p))
            $PT(p.particles[[rand(rng, 1:nparticles(p)) for _ in 1:n]])
        end

        function bootstrap(rng::AbstractRNG, p::AbstractArray{<:$PT}, n::Integer = nparticles(p))
            inds = [rand(rng, 1:nparticles(p)) for _ in 1:n]
            newpart = [p.particles[inds] for p in p]
            $PT.(newpart)
        end
    end
end
bootstrap(p::T, n::Integer = nparticles(p)) where T <: AbstractParticles = bootstrap(Random.GLOBAL_RNG, p, n)
bootstrap(p::MvParticles, n::Integer = nparticles(p)) = bootstrap(Random.GLOBAL_RNG, p, n)


text:
using Base.Test
import WaterData

# preliminary type definitions

"type definitions for testing piecewise EOS"
module test_eos_combined_resources
import WaterData

"1D test EOS"
struct TestEOS1D <: WaterData.FunctionalEOS
    power::Int
end
(e::TestEOS1D)(P::Real) = P^(e.power)

"2D test EOS"
struct TestEOS2D <: WaterData.FunctionalEOS
    power::Int
    add::Int
end
(e::TestEOS2D)(P::Real, T::Real) = (P+T)^(e.power) + (e.add)

end # module resources

# test

@testset "Combined EOS" begin
    res = test_eos_combined_resources

    @testset "1D piecewise EOS" begin
        # three different versions
        identity = res.TestEOS1D(1)
        squared = res.TestEOS1D(2)
        cubed = res.TestEOS1D(3)

        # make piecewise EOS
        eos = WaterData.PressurePiecewiseEOS(
            [identity, squared, cubed], [4, 7, 9, 11])

        @testset "Get correct EOS from the piecewise EOS" begin
            OutOfDomain = WaterData.OutOfDomainEOS

            @test typeof(WaterData.extracteos(eos, 2)) == OutOfDomain
            @test typeof(WaterData.extracteos(eos, 5)) == res.TestEOS1D
            @test typeof(WaterData.extracteos(eos, 12)) == OutOfDomain
        end

        @testset "Evaluate the piecewise EOS to the correct values" begin
            @test_throws DomainError eos(2)
            @test_throws DomainError eos(14)
            @test eos(5) == 5
            @test eos(8) == 8^2
            @test eos(10) == 10^3
        end
    end

    @testset "2D stitched EOS" begin
        # three diferent versions
        identity = res.TestEOS2D(1, 0)
        sqr1 = res.TestEOS2D(2, 1)
        cube2 = res.TestEOS2D(3, 2)

        # bounding boxes that overlap a bit
        BB1 = WaterData.BoundingBox(0, 5, 0, 5)
        BB2 = WaterData.BoundingBox(0, 7, 3, 10)
        BB3 = WaterData.BoundingBox(3, 10, 3, 10)

        # apply the bounds
        e1 = WaterData.BoundedEOS(identity, BB1)
        e2 = WaterData.BoundedEOS(sqr1, BB2)
        e3 = WaterData.BoundedEOS(cube2, BB3)

        # stitch together
        eos = WaterData.StitchedEOS([e1, e2, e3])

        @testset "Values outside the domain are errors" begin
            @test_throws DomainError eos(-3, -7)
            @test_throws DomainError eos(11, 12)
        end

        @testset "Correct values within the domain" begin
            @test eos(2, 2) == (2+2)^1 + 0  # first EOS
            @test eos(2, 8) == (2+8)^2 + 1  # second EOS
            @test eos(8, 8) == (8+8)^3 + 2  # third EOS
        end

        @testset "Correct priority for overlapping EOS" begin
            @test eos(4, 4) == (4+4)^1 + 0  # first EOS only
            @test eos(4, 8) == (4+8)^2 + 1  # second EOS only
            @test eos(2, 8) == (2+8)^2 + 1  # ditto
        end
    end
end


text:
ts = "2019-03-23T23:10:00"
te = "2019-03-23T23:17:00"
chans = [ "CC.JRO..BHZ",
          "CC.LON..BHZ",
          "CC.VALT..BHZ",
          "UW.HOOD..ENZ",
          "UW.HOOD..ENN",
          "UW.HOOD..ENE" ]

printstyled("  IRISWS continuous data requests\n", color=:light_green)
printstyled("    SAC and GeoCSV requests\n", color=:light_green)
for i = 1:length(chans)
  cha = chans[i]
  # println("cha = ", cha)
  S = get_data("IRIS", cha, src="IRIS", s=ts, t=te, fmt="sacbl", v=0, w=true)
  basic_checks(S)
  sleep(1)
  T = get_data("IRIS", cha, src="IRIS", s=ts, t=te, fmt="mseed", v=0, w=true)
  basic_checks(T)
  U = get_data("IRIS", cha, src="IRIS", s=ts, t=te, fmt="geocsv", v=0)
  basic_checks(U)

  if S.n == 0 || T.n == 0 || U.n == 0
    @warn(string(cha, " appears to be offline; trying next."))
    if i == lastindex(chans)
      error("No data for any station; failing test due to connection errors.")
    end
  else
    printstyled("     SAC == MSEED == GeoCSV\n", color=:light_green)
    sync!(S, s=ts, t=te)
    sync!(T, s=ts, t=te)
    sync!(U, s=ts, t=te)
    @test S.x[1] ≈ T.x[1] ≈ U.x[1]

    # :notes and :src will be different; as of 2020-05-28, so will :loc
    # :units is only set in GeoCSV
    for f in Symbol[:id, :name, :fs, :gain, :resp, :t]
      @test getfield(S,f) == getfield(T,f) == getfield(U,f)
    end

    #= Change 2020-05-28
      after IRISWS timeseries issues (first reported to IRIS 2020-05-16),
      :loc and :gain are being set in SAC requests.
    =#

    # printstyled("    GeoCSV\n", color=:light_green)
    break
  end
end

# Test bad data formats
printstyled("    bad request logging\n", color=:light_green)
sta_matrix = vcat(["IU" "ANMO" "00" "BHZ"],["IU" "ANMO" "00" "BHE"])
test_sta = deepcopy(sta_matrix)
ts = "2005-01-01T00:00:00"
te = "2005-01-02T00:00:00"
redirect_stdout(out) do
  try
    # this is an unparseable format
    S = get_data("IRIS", sta_matrix, s=ts, t=te, v=2, fmt="audio")

    # this is a bad request due to the wild card
    get_data!(S, "IRIS", "IU.ANMO.*.*", s=ts, t=te, v=2, fmt="ascii")

    @test S.n == 3

    @test S.id[2] == "XX.FAIL..001"
    @test any([occursin("request failed", n) for n in S.notes[2]])
    @test haskey(S.misc[2], "msg")

    @test S.id[1] == "XX.FMT..001"
    @test any([occursin("unparseable format", n) for n in S.notes[1]])
    @test haskey(S.misc[1], "raw")
  catch err
    @warn(string("Bad request logging test failed; caught error ", err))
  end
end

# check that these aren't modified in-place by the request (very old bug)
@test sta_matrix == test_sta

printstyled("    complicated IRISWS request\n", color=:light_green)
chans = ["UW.TDH..EHZ", "UW.VLL..EHZ", "CC.JRO..BHZ"] # HOOD is either offline or not on IRISws right now
st = -86400.0
en = -86100.0
(d0,d1) = parsetimewin(st,en)
S = get_data("IRIS", chans, s=d0, t=d1, y=true)
if isempty(S)
  @warn(string("No data for channels ", join(chans, ", "), "; test skipped."))
else
  L = [length(S.x[i])/S.fs[i] for i = 1:S.n]
  t = [S.t[i][1,2] for i = 1:S.n]
  L_min = minimum(L)
  L_max = maximum(L)
  t_min = minimum(t)
  t_max = maximum(t)
  try
    @test(L_max - L_min <= maximum(2 ./ S.fs))
    @test(t_max - t_min <= round(Int64, sμ * 2.0/maximum(S.fs)))
  catch
    @warn(string("Unexpected request length; check for partial outage at IRIS for ", join(chans, ", "), "!"))
  end
end


text:
### 3.1 - Loops

n = 0
m = 1

#=
# Major problems. n is not recognized inside the while loop. This makes the loop unusable
while n < 10
    n += 1
    println(m)
end
=#

for i ∈ 0:10
    println(i)
    println(n)
end

# Q1

for n in 1:100
    print(n^2, ", ")
end

# Q2

sq_dict = Dict{Int64, Int64}()

for n in 1:100
    print(n^2, ", ")
    sq_dict[n] = n^2
end

# Q3

squares_arr = [i^2 for i in 1:100]

@assert length(squares_arr) == 100
@assert sum(squares_arr) == 338350

### 3.2 - Conditionals

# Q1

x = 4

if (x % 2) == 0
    println(x)
else
    println("odd")
end

# Q2

(x % 2) == 1 ? println(x) : println("odd")

### 3.3 - Functions

# Q1

add_one = x -> x+1

# Q2 & Q3

M = [i+j for i in 1:3, j in 1:3]
add_one.(M)

max_stars_repo_path: Exercise Solutions/3 - Control Flow.jl

max_issues_repo_path: Exercise Solutions/3 - Control Flow.jl

max_forks_repo_path: Exercise Solutions/3 - Control Flow.jl


text:
using LBFGSB
using Test

@testset "driver2" begin
# translated from driver2.f
# "This driver shows how to replace the default stopping test by other termination
#  criteria. It also illustrates how to print the values of several parameters
#  during the course of the iteration. The sample problem used here is the same
#  as in DRIVER1 (the extended Rosenbrock function with bounds on the variables)."

nmax = 1024   # the dimension of the largest problem to be solved
mmax = 17    # the maximum number of limited memory corrections

task = fill(Cuchar(' '), 60)    # fortran's blank padding
csave = fill(Cuchar(' '), 60)    # fortran's blank padding
lsave = zeros(Cint, 4)
isave = zeros(Cint, 44)
dsave = zeros(Cdouble, 29)

iprint = -1    # suppress the default output

# "f is a DOUBLE PRECISION variable.  If the routine setulb returns with task(1:2)= 'FG',
#  then f must be set by the user to contain the value of the function at the point x."
f = 0.0

# "g is a DOUBLE PRECISION array of length n.  If the routine setulb returns with taskb(1:2)= 'FG',
#  then g must be set by the user to contain the components of the gradient at the point x."
g = zeros(Cdouble, nmax)

wa = zeros(Cdouble, 2mmax*nmax + 5nmax + 11mmax*mmax + 8mmax)
iwa = zeros(Cint, 3*nmax)

# "We suppress both code-supplied stopping tests because the user is providing his own stopping criteria."
factr = 0.0
pgtol = 0.0

n = 25    # the dimension n of the sample problem
m = 5     # the number of m of limited memory correction stored

# provide nbd which defines the bounds on the variables:
nbd = zeros(Cint, nmax)
l = zeros(Cdouble, nmax)    # the lower bounds
u = zeros(Cdouble, nmax)    # the upper bounds
# "First set bounds on the odd-numbered variables."
for i = 1:2:n
    nbd[i] = 2
    l[i] = 1e0
    u[i] = 1e2
end
# "Next set bounds on the even-numbered variables."
for i = 2:2:n
    nbd[i] = 2
    l[i] = -1e2
    u[i] = 1e2
end

# "We now define the starting point."
x = zeros(Cdouble, nmax)
for i = 1:n
    x[i] = 3e0
end

println("Solving sample problem.")
println(" (f = 0.0 at the optimal solution.)")

# "We start the iteration by initializing task."
task[1:5] = b"START"

## ------- the beginning of the loop ----------
let
    @label CALLLBFGSB

    # This is the call to the L-BFGS-B code.
    setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)

    # "the minimization routine has returned to request the function f and gradient g values at the current x."
    if task[1:2] == b"FG"
        # "Compute function value f for the sample problem."
        f = 0.25 * (x[1] - 1)^2
        for i = 2:n
            f = f + (x[i] - x[i-1]^2)^2
        end

        f = 4 * f

        # "Compute gradient g for the sample problem."
        t₁ = x[2] - x[1]^2
        g[1] = 2 * (x[1] - 1) - 1.6e1 * x[1] * t₁
        for i = 2:n-1
            t₂ = t₁
            t₁ = x[i+1] - x[i]^2
            g[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
        end
        g[n] = 8 * t₁

        # "go back to the minimization routine."
        @goto CALLLBFGSB
    end

    # "the minimization routine has returned with a new iterate, At this point have the
    #  opportunity of stopping the iteration or observing the values of certain parameters"
    if task[1:5] == b"NEW_X"
        # "First are two examples of stopping tests."

        # "Note: task(1:4) must be assigned the value 'STOP' to terminate the iteration
        #  and ensure that the final results are printed in the default format. The
        #  rest of the character string TASK may be used to store other information.

        # "1) Terminate if the total number of f and g evaluations exceeds 99."
        isave[34] ≥ 99 && (task[1:4] = b"STOP"; terminate_info="TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT")

        # "2) Terminate if  |proj g|/(1+|f|) < 1.0d-10, where "proj g" denoted the projected gradient"
        dsave[13] ≤ 1e-10 * (1e0 + abs(f)) && (task[1:4] = b"STOP"; terminate_info="THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL")

        # "We now wish to print the following information at each iteration:"
        # "1) the current iteration number, isave(30),"
        # "2) the total number of f and g evaluations, isave(34),"
        # "3) the value of the objective function f,"
        # "4) the norm of the projected gradient,  dsave(13)"

        println("Iterate ", isave[30], "  nfg = ", isave[34], "  f = ", f, "  |proj g| = ", dsave[13])

        # "If the run is to be terminated, we print also the information contained
        #  in task as well as the final value of x."
        if task[1:4] == b"STOP"
            println(terminate_info)
            for i = 1:n
                println("x($i) = ", x[i])
            end
        end

        # "go back to the minimization routine."
        @goto CALLLBFGSB
    end

    # ---------- the end of the loop -------------
    # "If task is neither FG nor NEW_X we terminate execution."
end

tol = 1e-4

@test all(x[1:10] .- 1.0000 .< tol)
@test x[11] - 1.0001 < tol
@test x[12] - 1.0002 < tol
@test x[13] - 1.0003 < tol
@test x[14] - 1.0006 < tol
@test x[15] - 1.0013 < tol
@test x[16] - 1.0026 < tol
@test x[17] - 1.0052 < tol
@test x[18] - 1.0105 < tol
@test x[19] - 1.0210 < tol
@test x[20] - 1.0425 < tol
@test x[21] - 1.0867 < tol
@test x[22] - 1.1810 < tol
@test x[23] - 1.3947 < tol
@test x[24] - 1.9452 < tol
@test x[25] - 3.7837 < tol

end # EOB


text:
using Interpolations, DifferentialEquations, Turing, FillArrays, ReverseDiff

# Define ODE
function mesocosm_metabolism!(du,u,p,t)
    du[1] = p["PAR"](t)*(p["GPP_daily"])[p["day"](t)] +
            (p["ER_daily"])[p["day"](t)] +
            (p["kGas"] + (p["Q"]/27.0))*(p["OSat"](t) - u[1])
end

# Generate some fake data
times = collect(4.0:0.1:75.0)
nDays = 3
dayStart = 4
days = [convert(Int64, floor(times[i]/24.0 - (dayStart/24.0)) + 1) for i in 1:length(times)]
PARData= [if (7<= (times[i] - 24*(days[i]-1)) < 19) 1 else 0 end for i in 1:length(times)]
OSatData = repeat([8.2], length(times))

p = Dict(
    "PAR" => Interpolations.ConstantInterpolation((times,), PARData),
    "OSat" => Interpolations.interpolate((times,), OSatData, Interpolations.Gridded(Interpolations.Linear())),
    "day" => Interpolations.ConstantInterpolation((times,), days),
    "GPP_daily"=> fill(1.8, nDays),
    "ER_daily" => fill(-1.8, nDays),
    "kGas" => 0.6,
    "Q" => 1.794)

tspan = (minimum(times), maximum(times))
u0 = [5.5]

problem = ODEProblem(mesocosm_metabolism!, u0, tspan, p)
sol = solve(problem, alg = AutoVern7(Rodas4()),reltol = 1e-5, maxiters = 1e6, tstops = [7,19], saveat = times)
fakeData = vec(sol) + 0.2*rand(length(times))

# Specify statistical model
@model function f(::Type{T} = Float64; DOdata, ODE) where {T}
    # Priors on daily GPP and ER, k600, observation error, and process error
    GPP_daily ~ MvNormal(FillArrays.Fill(1.80, $nDays), 2.0)
    ER_daily ~ MvNormal(FillArrays.Fill(-1.80, $nDays), 2.0)
    k600 ~ LogNormal(-0.69, 0.53)
    err_obs_iid_sigma ~ truncated(Cauchy(0, 0.03), 0, Inf)
    err_proc_iid_sigma ~ truncated(Cauchy(0, 5.0), 0, Inf)

    # ODE parameters
    p = Dict(
        "PAR" => ODE.p["PAR"],
        "OSat" => ODE.p["OSat"],
        "day" => ODE.p["day"],
        "GPP_daily"=> GPP_daily,
        "ER_daily" => ER_daily,
        "kGas" => k600,
        "Q" => ODE.p["Q"])

    DO_true = Vector{T}(undef, (length(times)))
    DO_true[1] ~ Normal(DOdata[1], err_obs_iid_sigma)

    for i = 2:length(times)
        problem = remake(ODE, p = p, u0 = [DO_true[i-1]], tspan = (times[i-1], times[i]))
        predicted = solve(problem, alg = Euler(), dt = 0.1)
        DO_true[i] ~ Normal(predicted(times[i])[1], err_proc_iid_sigma)
        DOdata[i] ~ Normal(DO_true[i], err_obs_iid_sigma)
    end
end

model = f(DOdata = fakeData, ODE = problem)
setadbackend(:reversediff)
samples = sample(model, NUTS(0.65), 10)

@model function g(::Type{T} = Float64; DOdata, ODE, times, PAR, OSat, day, Q, V) where {T}

   # Priors on daily GPP and ER, k600, observation error, and process error
   GPP_daily ~ MvNormal(FillArrays.Fill(1.80, $nDays), 2.0)
   ER_daily ~ MvNormal(FillArrays.Fill(-1.80, $nDays), 2.0)
   kGas ~ LogNormal(-0.69, 0.53)
   err_obs_iid_sigma ~ truncated(Cauchy(0, 0.03), 0, Inf)
   err_proc_iid_sigma ~ truncated(Cauchy(0, 5.0), 0, Inf)

   DO_mod = Vector{T}(undef, (length(times)))
   DO_true = Vector{T}(undef, (length(times)))
   DO_true[1] ~ Normal(DOdata[1], err_obs_iid_sigma)

   for i = 2:length(times)
       DO_mod[i] = DO_true[i-1] + (times[i] - times[i-1])*
        (PAR(times[i-1])*GPP_daily[day(times[i-1])] +
        ER_daily[day(times[i-1])] +
        (kGas + Q/V)*(OSat(times[i-1]) - DO_true[i-1]))
       DO_true[i] ~ Normal(DO_mod[i], err_proc_iid_sigma)
       DOdata[i] ~ Normal(DO_true[i], err_obs_iid_sigma)
   end

end

V = 27.0
model = g(DOdata = fakeData, ODE = problem, times = times, PAR = PAR, OSat = OSat, day = day, Q = Q, V = V)
samples = sample(model, NUTS(0.65), 10)

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
# pointer.jl
#
# Julia type and memory management for CUDA pointers.  Combined
# from CUDA.jl, CUDArt.jl, and Nick Henderson's fork of CUDArt.jl.

# A raw CUDA pointer
type CudaPtr{T}
    ptr::Ptr{T}
end

# Type alias for previous name
typealias CudaDevicePtr CudaPtr

#############################
# Low-level memory handling #
#############################

CudaPtr() = CudaPtr(C_NULL)
CudaPtr(T::Type) = CudaPtr(unsafe_convert(Ptr{T},C_NULL))
unsafe_convert{P<:Ptr}(::Type{P}, p::CudaPtr) = unsafe_convert(P, p.ptr)
convert{P<:Ptr}(::Type{P}, p::CudaPtr) = unsafe_convert(P, p.ptr)
copy(p::CudaPtr) = CudaPtr(p.ptr)

rawpointer(p::CudaPtr) = p

# Enable both manual and garbage-collected memory management.
# If you need to free resources, you can call free manually.
# cuda_ptrs keeps track of all memory that needs to be freed,
# and prevents double-free (which otherwise causes serious problems).
# key = ptr, val = device id
const cuda_ptrs = Dict{Any,Int}()

function malloc(T::Type, n::Integer)
    p = Ptr{Void}[C_NULL]
    nbytes = sizeof(T)*n
    rt.cudaMalloc(p, nbytes)
    cptr = CudaPtr(unsafe_convert(Ptr{T},p[1]))
    finalizer(cptr, free)
    cuda_ptrs[WeakRef(cptr)] = device()
    cptr
end
malloc(nbytes::Integer) = malloc(UInt8, nbytes)

function free{T}(p::CudaPtr{T})
    cnull = unsafe_convert(Ptr{T}, C_NULL)
    if p.ptr != cnull && haskey(cuda_ptrs, p)
        delete!(cuda_ptrs, p)
        rt.cudaFree(p)
        p.ptr = cnull
    end
end

typealias Ptrs Union{Ptr, CudaPtr, rt.cudaPitchedPtr}
typealias CudaPtrs Union{CudaPtr, rt.cudaPitchedPtr}

cudamemcpykind(dstp::Ptr, srcp::Ptr) = rt.cudaMemcpyHostToHost
cudamemcpykind(dstp::CudaPtrs, srcp::Ptr) = rt.cudaMemcpyHostToDevice
cudamemcpykind(dstp::Ptr, srcp::CudaPtrs) = rt.cudaMemcpyDeviceToHost
cudamemcpykind(dstp::CudaPtrs, srcp::CudaPtrs) = rt.cudaMemcpyDeviceToDevice
cudamemcpykind(dst::Ptrs, src::Ptrs) = error("This should never happen") # prevent a useless ambiguity warning
cudamemcpykind(dst, src::Ptrs) = cudamemcpykind(pointer(dst), src)
cudamemcpykind(dst::Ptrs, src) = cudamemcpykind(dst, pointer(src))
cudamemcpykind(dst, src) = cudamemcpykind(pointer(dst), pointer(src))

## converting pointers to an appropriate unsigned ##

const CUDA_NULL = CudaPtr()

# pointer to integer
convert(::Type{UInt}, x::CudaPtr) = convert(UInt,x.ptr)
convert{T<:Integer}(::Type{T}, x::CudaPtr) = convert(T,unsigned(x))

# integer to pointer
convert{T}(::Type{CudaPtr{T}}, x::Integer) = CudaPtr(convert(Ptr{T},x))

# pointer to pointer
convert{T}(::Type{CudaPtr{T}}, p::CudaPtr{T}) = p
convert{T}(::Type{CudaPtr{T}}, p::CudaPtr) = CudaPtr(convert(Ptr{T},p.ptr))

integer(x::CudaPtr) = convert(UInt, x.ptr)
unsigned(x::CudaPtr) = convert(UInt, x.ptr)

eltype{T}(::CudaPtr{T}) = T

## limited pointer arithmetic & comparison ##

==(x::CudaPtr, y::CudaPtr) = UInt(x) == UInt(y)
-(x::CudaPtr, y::CudaPtr) = UInt(x) - UInt(y)

+(x::CudaPtr, y::Integer) = oftype(x, UInt(UInt(x) + y))
-(x::CudaPtr, y::Integer) = oftype(x, UInt(UInt(x) - y))
+(x::Integer, y::CudaPtr) = y + x

zero{T}(::Type{CudaPtr{T}}) = unsafe_convert(CudaPtr{T}, C_NULL)
zero{T}(x::CudaPtr{T}) = unsafe_convert(CudaPtr{T}, C_NULL)
one{T}(::Type{CudaPtr{T}}) = unsafe_convert(CudaPtr{T}, 1)
one{T}(x::CudaPtr{T}) = unsafe_convert(CudaPtr{T}, 1)


text:
stats = SolverStats(parent=:ALTRO)
N = 10
Altro.reset!(stats, N)
@test length(stats.cost) == N
@test length(stats.iteration) == N
@test sum(stats.cost) == 0
@test sum(stats.iteration) == 0

Altro.record_iteration!(stats, cost=10, c_max=100)
@test stats.cost[1] == 10
@test stats.c_max[1] == 100
@test stats.gradient[1] == 0.0
@test stats.penalty_max[1] == 0.0
@test stats.iteration[1] == 1
@test stats.iteration_outer[1] == 0
@test stats.iteration_pn[1] == false
@test stats.cost[2] == 0
@test stats.c_max[2] == 0

stats.iterations_outer += 1
Altro.record_iteration!(stats, cost=5, penalty_max=1e4)
@test stats.cost[2] == 5
@test stats.c_max[2] == 100
@test stats.gradient[2] == 0.0
@test stats.penalty_max[2] ≈ 1e4
@test stats.iteration[2] == 2
@test stats.iteration_outer[2] == 1
@test stats.iteration_pn[2] == false

Altro.record_iteration!(stats, c_max=2, gradient=1e-1, is_pn=true)
@test stats.cost[3] == 5
@test stats.c_max[3] == 2
@test stats.gradient[3] == 1e-1
@test stats.penalty_max[3] == 1e4
@test stats.iteration[3] == 3
@test stats.iteration_outer[3] == 1
@test stats.iteration_pn[3] == true
@test stats.iterations_pn == 1
@test stats.dJ_zero_counter == 0

Altro.trim!(stats)
@test length(stats.cost) == 3
@test stats.c_max[end] == 2
@test stats.gradient[end] == 1e-1

@test (@ballocated Altro.reset!(stats, 10_000) samples=10 evals=10) == 0
@test (@ballocated Altro.record_iteration!($stats,
    c_max=2, gradient=1e-1, is_pn=true) samples=10 evals=10) == 0

solver = ALTROSolver(Problems.DoubleIntegrator()...)
set_options!(solver, verbose=0, projected_newton=true, show_summary=false)
# solver.opts.verbose = 2
Z0 = copy(get_trajectory(solver))
solve!(solver)

# Check if final stats match output
iters = iterations(solver)
c_max = solver.stats.c_max[end]
J = solver.stats.cost[end]
@test c_max ≈ max_violation(solver.solver_pn)
@test J ≈ cost(solver)

# Test stats to dictionary
@test Dict(solver.stats) isa Dict{Symbol,Vector{<:Any}}

# Make sure ALTRO resets properly
initial_trajectory!(solver, Z0)

solve!(solver)
@test iterations(solver) == iters
@test cost(solver) ≈ J
@test max_violation(solver.solver_pn) ≈ c_max

# Make sure AL solver resets properly
solver = Altro.ALSolver(Problems.DoubleIntegrator()..., show_summary=false)
Z0 = copy(get_trajectory(solver))
solve!(solver)
J = cost(solver)
iters = iterations(solver)
initial_trajectory!(solver, Z0)
Altro.reset!(solver)
solve!(solver)
@test J ≈ cost(solver)
@test iters == iterations(solver)


text:
stress_centrality(g::AbstractGraph, vs=vertices(g); parallel=:distributed) =
parallel == :distributed ? distr_stress_centrality(g, vs) : threaded_stress_centrality(g, vs)

stress_centrality(g::AbstractGraph, k::Integer; parallel=:distributed) =
parallel == :distributed ? distr_stress_centrality(g, sample(vertices(g), k)) :
threaded_stress_centrality(g, sample(vertices(g), k))

function distr_stress_centrality(g::AbstractGraph,
    vs=vertices(g))::Vector{Int64}

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    # Parallel reduction
    stress = @distributed (+) for s in vs
        temp_stress = zeros(Int64, n_v)
        if degree(g, s) > 0  # this might be 1?
            state = LightGraphs.dijkstra_shortest_paths(g, s; allpaths=true, trackvertices=true)
            LightGraphs._stress_accumulate_basic!(temp_stress, state, g, s)
        end
        temp_stress
    end
    return stress
end

function threaded_stress_centrality(g::AbstractGraph,
    vs=vertices(g))::Vector{Int64}

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    # Parallel reduction
    local_stress = [zeros(Int, n_v) for _ in 1:nthreads()]

    Base.Threads.@threads for s in vs
        if degree(g, s) > 0  # this might be 1?
            state = LightGraphs.dijkstra_shortest_paths(g, s; allpaths=true, trackvertices=true)
            LightGraphs._stress_accumulate_basic!(local_stress[Base.Threads.threadid()], state, g, s)
        end
    end
    return reduce(+, local_stress)
end


text: level(r::Reservoir, volume) = level(bathymetry(r), volume)


text:
module HolyTraitPattern

export Residence, Stock, TreasuryBill, Money, Book
export tradable, marketprice

include("asset_types.jl")
include("traits.jl")
include("literature.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
using Pkg, DrWatson

# ╔═╡ e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random
	Random.seed!(1)
end

# ╔═╡ ac3ceabe-1ee4-11eb-3935-313e96aafd62
md"## Listing 2.13"

# ╔═╡ e9fa266e-1eeb-11eb-2630-59930c474ed3
function montyHall(switchPolicy)
    prize, choice = rand(1:3), rand(1:3)
    if prize == choice
        revealed = rand(setdiff(1:3,choice))
    else
        revealed = rand(setdiff(1:3,[prize,choice]))
    end

    if switchPolicy
        choice = setdiff(1:3,[revealed,choice])[1]
    end
    return choice == prize
end

# ╔═╡ 9046a9c8-1f64-11eb-067f-6552cbd8919e
N = 10^6

# ╔═╡ 9046fdce-1f64-11eb-3888-651b7d2d17ac
Text("Success probability with policy I (stay): $(sum([montyHall(false) for _ in 1:N])/N)")

# ╔═╡ 90479df6-1f64-11eb-3e06-af4adc28a9d5
Text("Success probability with policy II (switch): $(sum([montyHall(true) for _ in 1:N])/N)")

# ╔═╡ ea086874-1eeb-11eb-2cb0-4186697d4294
md"## End of listing 2.13"

# ╔═╡ Cell order:
# ╟─ac3ceabe-1ee4-11eb-3935-313e96aafd62
# ╠═de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
# ╠═e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
# ╠═e9fa266e-1eeb-11eb-2630-59930c474ed3
# ╠═9046a9c8-1f64-11eb-067f-6552cbd8919e
# ╠═9046fdce-1f64-11eb-3888-651b7d2d17ac
# ╠═90479df6-1f64-11eb-3e06-af4adc28a9d5
# ╟─ea086874-1eeb-11eb-2cb0-4186697d4294


text:
using Oceananigans.BoundaryConditions: PBC, ZFBC, NFBC, ContinuousBoundaryFunction, DiscreteBoundaryFunction

using Oceananigans.Fields: Face, Center

simple_bc(ξ, η, t) = exp(ξ) * cos(η) * sin(t)

function can_instantiate_boundary_condition(bc, C, FT=Float64, ArrayType=Array)
    success = try
        bc(C, FT, ArrayType)
        true
    catch
        false
    end
    return success
end

@testset "Boundary conditions" begin
    @info "Testing boundary conditions..."

    @testset "Boundary condition instantiation" begin
        @info "  Testing boundary condition instantiation..."

        for C in (Value, Gradient, Flux)
            @test can_instantiate_boundary_condition(integer_bc, C)
            @test can_instantiate_boundary_condition(irrational_bc, C)
            @test can_instantiate_boundary_condition(simple_function_bc, C)
            @test can_instantiate_boundary_condition(parameterized_function_bc, C)
            @test can_instantiate_boundary_condition(field_dependent_function_bc, C)
            @test can_instantiate_boundary_condition(discrete_function_bc, C)
            @test can_instantiate_boundary_condition(parameterized_discrete_function_bc, C)

            for FT in float_types
                @test can_instantiate_boundary_condition(float_bc, C, FT)
                @test can_instantiate_boundary_condition(parameterized_field_dependent_function_bc, C, FT)

                for arch in archs
                    ArrayType = array_type(arch)
                    @test can_instantiate_boundary_condition(array_bc, C, FT, ArrayType)
                end
            end
        end
    end

    @testset "Field and coordinate boundary conditions" begin
        @info "  Testing field and coordinate boundary conditions..."

        # Triply periodic
        ppp_topology = (Periodic, Periodic, Periodic)
        ppp_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=ppp_topology)

        u_bcs = UVelocityBoundaryConditions(ppp_grid)
        v_bcs = VVelocityBoundaryConditions(ppp_grid)
        w_bcs = WVelocityBoundaryConditions(ppp_grid)
        T_bcs = TracerBoundaryConditions(ppp_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa PBC
        @test u_bcs.y.right isa PBC
        @test u_bcs.z.left  isa PBC
        @test u_bcs.z.right isa PBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa PBC
        @test v_bcs.y.right isa PBC
        @test v_bcs.z.left  isa PBC
        @test v_bcs.z.right isa PBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa PBC
        @test w_bcs.y.right isa PBC
        @test w_bcs.z.left  isa PBC
        @test w_bcs.z.right isa PBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa PBC
        @test T_bcs.y.right isa PBC
        @test T_bcs.z.left  isa PBC
        @test T_bcs.z.right isa PBC

        # Doubly periodic. Engineers call this a "Channel geometry".
        ppb_topology = (Periodic, Periodic, Bounded)
        ppb_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=ppb_topology)

        u_bcs = UVelocityBoundaryConditions(ppb_grid)
        v_bcs = VVelocityBoundaryConditions(ppb_grid)
        w_bcs = WVelocityBoundaryConditions(ppb_grid)
        T_bcs = TracerBoundaryConditions(ppb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa PBC
        @test u_bcs.y.right isa PBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa PBC
        @test v_bcs.y.right isa PBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa PBC
        @test w_bcs.y.right isa PBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa PBC
        @test T_bcs.y.right isa PBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        # Singly periodic. Oceanographers call this a "Channel", engineers call it a "Pipe"
        pbb_topology = (Periodic, Bounded, Bounded)
        pbb_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=pbb_topology)

        u_bcs = UVelocityBoundaryConditions(pbb_grid)
        v_bcs = VVelocityBoundaryConditions(pbb_grid)
        w_bcs = WVelocityBoundaryConditions(pbb_grid)
        T_bcs = TracerBoundaryConditions(pbb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa ZFBC
        @test u_bcs.y.right isa ZFBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa NFBC
        @test v_bcs.y.right isa NFBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa ZFBC
        @test w_bcs.y.right isa ZFBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa ZFBC
        @test T_bcs.y.right isa ZFBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        # Triply bounded. Oceanographers call this a "Basin", engineers call it a "Box"
        bbb_topology = (Bounded, Bounded, Bounded)
        bbb_grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1), topology=bbb_topology)

        u_bcs = UVelocityBoundaryConditions(bbb_grid)
        v_bcs = VVelocityBoundaryConditions(bbb_grid)
        w_bcs = WVelocityBoundaryConditions(bbb_grid)
        T_bcs = TracerBoundaryConditions(bbb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa NFBC
        @test u_bcs.x.right isa NFBC
        @test u_bcs.y.left  isa ZFBC
        @test u_bcs.y.right isa ZFBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa ZFBC
        @test v_bcs.x.right isa ZFBC
        @test v_bcs.y.left  isa NFBC
        @test v_bcs.y.right isa NFBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa ZFBC
        @test w_bcs.x.right isa ZFBC
        @test w_bcs.y.left  isa ZFBC
        @test w_bcs.y.right isa ZFBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa ZFBC
        @test T_bcs.x.right isa ZFBC
        @test T_bcs.y.left  isa ZFBC
        @test T_bcs.y.right isa ZFBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        grid = bbb_grid

        T_bcs = TracerBoundaryConditions(grid;
                                           east = BoundaryCondition(Value, simple_bc),
                                           west = BoundaryCondition(Value, simple_bc),
                                         bottom = BoundaryCondition(Value, simple_bc),
                                            top = BoundaryCondition(Value, simple_bc),
                                          north = BoundaryCondition(Value, simple_bc),
                                          south = BoundaryCondition(Value, simple_bc)
                                         )

        # Note that boundary condition setting is valid only in cases where
        # boundary conditionsare same type

        one_bc = BoundaryCondition(Value, 1.0)
        two_bc = BoundaryCondition(Value, 2.0)

        T_bcs = TracerBoundaryConditions(grid;
                                           east = one_bc,
                                           west = one_bc,
                                         bottom = one_bc,
                                            top = one_bc,
                                          north = one_bc,
                                          south = one_bc,
                                         )

        @test T_bcs.east   === one_bc
        @test T_bcs.west   === one_bc
        @test T_bcs.north  === one_bc
        @test T_bcs.south  === one_bc
        @test T_bcs.top    === one_bc
        @test T_bcs.bottom === one_bc

        T_bcs.east   = two_bc
        T_bcs.west   = two_bc
        T_bcs.north  = two_bc
        T_bcs.south  = two_bc
        T_bcs.top    = two_bc
        T_bcs.bottom = two_bc

        @test T_bcs.east   === two_bc
        @test T_bcs.west   === two_bc
        @test T_bcs.north  === two_bc
        @test T_bcs.south  === two_bc
        @test T_bcs.top    === two_bc
        @test T_bcs.bottom === two_bc
    end
end


text:
@testset "AlgAssOrd" begin
  include("AlgAssOrd/CSAMaxOrd.jl")
  include("AlgAssOrd/PicardGroup.jl")
  include("AlgAssOrd/LocallyFreeClassGroup.jl")
end


text:
sciplotlib_style = Dict(
    "axes.spines.top"      => false,
    "axes.spines.right"    => false,
    "axes.grid"            => true,
    "axes.axisbelow"       => true,  # Grid _below_ patches (such as histogram bars), not on top.
    "axes.grid.which"      => "both",
    "grid.linewidth"       => 0.5,        # These are for major grid. Minor grid styling
    "grid.color"           => "#E7E7E7",  # is set in `set!`.

    "xtick.direction"      => "in",
    "ytick.direction"      => "in",
    "xtick.labelsize"      => "small", # Default is "medium"
    "ytick.labelsize"      => "small", # idem
    "legend.fontsize"      => "small", # Default is "medium"
    "axes.titlesize"       => "medium",
    "axes.labelsize"       => 9,
    "xaxis.labellocation"  => "right",
    "axes.titlelocation"   => "right",

    "legend.borderpad"     => 0.6,
    "legend.borderaxespad" => 0.2,

    "lines.solid_capstyle" => "round",

    "figure.facecolor"     => "white",
    "figure.figsize"       => (4, 2.4),
    "figure.dpi"           => 200,
    "savefig.dpi"          => "figure",
    "savefig.bbox"         => "tight",

    "axes.autolimit_mode"  => "data" ,
    "axes.xmargin"         => 0,
    "axes.ymargin"         => 0,
)

"""
Reset Matplotlib's style to `rcParams_original`, then apply the supplied dictionary of
`rcParams` settings. Call without arguments to reset to Matplotlib's defaults. To reset to
Sciplotlib's defaults, pass `sciplotlib_style`.

The initial reset is so that you can experiment with parameters; namely add and then remove
entries. Without the reset, once a parameter was set it would stay set.
"""
function set_mpl_style!(updatedRcParams = nothing)
    merge!(rcParams, rcParams_original)
    isnothing(updatedRcParams) || merge!(rcParams, updatedRcParams)
    return rcParams
end


text:
target = "cpu"
using ExaPF
using ExaPF.Network
using DelimitedFiles

include("pf.jl")
datafile = "test/case14.raw"
gmresblocks =    [2,4,8]
bicgstabblocks = [2,4,8]
outfile = "./plots/case14.txt"
# datafile = "GO-Data/datasets/Trial_3_Real-Time/Network_13R-015/scenario_11/case.raw"
# gmresblocks =    [8, 16, 32, 64, 128]
# bicgstabblocks = [8, 16, 32, 64, 128]
# outfile = "./plots/Network_13R-015.txt"
# datafile = "GO-Data/datasets/Trial_3_Real-Time/Network_30R-025/scenario_1/case.raw"
# gmresblocks =    [16, 32, 64, 128, 256]
# bicgstabblocks = [16, 32, 64, 128, 256]
# outfile = "./plots/Network_30R-025.txt"
linsol_iters = []
first_iters = []
py = []
push!(py, gmresblocks)
for blocks in gmresblocks
  @time sol, conv, res, first_iter, linsol_iter = pf(datafile, blocks, "gmres")
  push!(linsol_iters, linsol_iter)
  push!(first_iters, first_iter)
  println("first gmres iteration: $first_iter, total gmres iterations: $linsol_iter")
end
push!(py, linsol_iters)
push!(py, first_iters)

linsol_iters = []
first_iters = []
for blocks in bicgstabblocks
  @time sol, conv, res, first_iter, linsol_iter = pf(datafile, blocks, "bicgstab")
  push!(linsol_iters, linsol_iter)
  push!(first_iters, first_iter)
  println("first bicgstab iteration: $first_iter, total bicgstab iterations: $linsol_iter")
end
push!(py, linsol_iters)
push!(py, first_iters)
f = open(outfile, "w+")
writedlm(f, py,',')
close(f)


text:
function get_nodes_from_raw(V_keys::Vector, T::Vector)

end

function get_arcs_from_raw(E,N_set)

end

function get_commodities_from_raw()

end


text:
using Plots
using Reproduce
using JLD2
using FileIO
using Statistics

include("plot_utils.jl")

function plot_ic_goal_visitation!(p, ic; step_size=2)
    data = load_results(ic, :visit_order, return_type="array")

    # optimal filtering of the median. This is NOT YET COMPLETE. Still working on it.
    # median_length = median([length(x) for x in data])
    # min_length = min([length(x) for x in data]...)
    # println("min: $(min_length) median: $(median([length(x) for x in data]))")

    # Trimming extra episodes for plotting sake
    min_length = min([length(x) for x in data]...)
    data = [x[1:min_length] for x in data]

    data = cat(data..., dims = 2)
    runs = size(data)[2]

    percentages = []
    for r in 1:runs
        percentage = get_single_goal_percentage(data[:,r]; step_size=step_size)
        push!(percentages, percentage)
    end

    percentages = cat(percentages..., dims = 3)
    percentages = dropdims(mean(percentages, dims=3), dims=3)

    for i in 1:4
        # println(size(percentages[i, :]))
        plot!(p, 1:size(percentages)[2], percentages[i, :])
    end
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module Gaugefield
    #include("../../Wilsonloops/src/Wilsonloops.jl")
    using Requires

    include("./output/verboseprint.jl")
    include("./SUN_generator.jl")
    include("./autostaples/wilsonloops.jl")
    include("./AbstractGaugefields.jl")
    include("./output/io.jl")
    include("./output/ildg_format.jl")
    include("./autostaples/Loops.jl")
    include("./smearing/Abstractsmearing.jl")

    function __init__()
        @require MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" begin
            import .AbstractGaugefields_module:identityGaugefields_4D_wing_mpi,
                    Gaugefields_4D_wing_mpi,calc_rank_and_indices,barrier,comm,
                    setvalue!
        end
    end

    import Wilsonloop:loops_staple_prime

    import .AbstractGaugefields_module:AbstractGaugefields,identitymatrix,Abstractfields,
                                        shift_U,construct_staple!,set_wing_U!,
                                        calculate_Plaquette,substitute_U!,calculate_Polyakov_loop,construct_gauges,
                                        Gaugefields_4D_wing,
                                        identityGaugefields_4D_wing,
                                        add_force!,exp_aF_U!,clear_U!,add_U!,exptU!,
                                        Traceless_antihermitian!,Traceless_antihermitian,Generator,
                                        Staggered_Gaugefields,staggered_U,
                                        Traceless_antihermitian_add!,
                                        IdentityGauges,RandomGauges,Oneinstanton,
                                        construct_Λmatrix_forSTOUT!,
                                        evaluate_gaugelinks_evenodd!,
                                        map_U!

    import .Loops_module:Loops,evaluate_loops,calc_large_wilson_loop!,evaluate_loops!
    import .Wilsonloops_module:Wilson_loop_set,make_staples,Wilson_loop_set,
                make_cloverloops,Tensor_derivative_set, make_loops,
                make_plaq_staple,make_links,make_plaq,
                make_loopforactions,make_plaqloops,make_rectloops,make_polyakovloops,
                make_plaq_staple_prime,
                calc_coordinate,make_plaq_staple_prime,calc_shift,
                Tensor_wilson_lines_set,Tensor_wilson_lines,Tensor_derivative_set,
                get_leftstartposition,get_rightstartposition,Wilson_loop,calc_loopset_μν_name,
                make_originalactions_fromloops
    import .AbstractGaugefields_module:TA_Gaugefields,initialize_TA_Gaugefields
    import .Abstractsmearing_module:Abstractsmearing,Nosmearing,Stoutsmearing,calc_smearedU,
            construct_smearing,Gradientflow,get_tempG,flow!,get_eps,back_prop,CovNeuralnet#gradientflow!
    #import .Verbose_print:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose3
    import .Verbose_print:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose3,println_verbose2,println_verbose1,
    print_verbose1,print_verbose2,print_verbose3
    import .ILDG_format:ILDG,load_gaugefield,load_gaugefield!,save_binarydata
    import .IOmodule:saveU,loadU,loadU!
    import .SUN_generator:Generator

end


text:
export RemotePtr, RemoteCodePtr
import Base: +, -

"""
Represents a debugging session with potentially one or more inferiors (e.g.
  multiple threads/processes in a progess group).
"""
abstract RemoteSession
immutable LocalSession; end

const SupportedWordSizes = Union{UInt32, UInt64}

immutable RemotePtr{T,intptrT<:SupportedWordSizes}
    ptr::intptrT
end
(::Type{RemotePtr{T,S} where S <: SupportedWordSizes})(arg) where {T} = RemotePtr{T,UInt64}(arg)
(::Type{RemotePtr{T,S} where S <: SupportedWordSizes})(arg::RemotePtr{<:Any,S′}) where {T, S′ <: SupportedWordSizes} = RemotePtr{T,S′}(arg)
(::Type{RemotePtr{<:Any,<:SupportedWordSizes}})(x::RemotePtr) = T(x.ptr)
Base.convert(::Type{UInt64}, x::RemotePtr{UInt64}) = UInt64(x.ptr)
Base.convert(::Type{UInt64}, x::RemotePtr) = UInt64(x.ptr)
Base.convert{T<:RemotePtr}(::Type{T}, x::T) = x
Base.convert{T<:RemotePtr}(::Type{T}, x::RemotePtr) =
    ((T === RemotePtr) ? x : T(x.ptr))::T
for f in (:+, :-)
    @eval ($f){T<:RemotePtr}(ptr::T,i::Integer) = T(($f)(ptr.ptr,i))
    @eval ($f){T<:RemotePtr}(i::Integer,ptr::T) = T(($f)(i,ptr.ptr))
end
Base.:(==)(x::Integer,y::RemotePtr) = x == y.ptr
Base.:(==)(x::RemotePtr,y::Integer) = x.ptr == y
Base.:(==)(x::RemotePtr,y::RemotePtr) = x.ptr == y.ptr
Base.isless(x::RemotePtr,y::RemotePtr) = isless(UInt64(x),UInt64(y))

# TODO: It's not entirely clear that this distinction is needed or useful,
# but it's present in RR, so is included in this abstract interface for
# compatibility
immutable RemoteCodePtr
    ptr::UInt64
end
Base.convert(::Type{UInt64},x::RemoteCodePtr) = x.ptr
Base.convert{T}(::Type{RemotePtr{Void,T}},x::RemoteCodePtr) = RemotePtr{Void,T}(x.ptr)
Base.convert{T}(::Type{RemotePtr{T}},x::RemoteCodePtr) = RemotePtr{T}(x.ptr)
Base.convert(::Type{RemotePtr},x::RemoteCodePtr) = RemotePtr{Void}(x.ptr)

function fallible_load
end

function load{T}(session::LocalSession, ptr::RemotePtr{T})
    unsafe_load(Ptr{T}(ptr.ptr))
end

function store!{T}(session::LocalSession, ptr::RemotePtr{T}, val::T)
    unsafe_store!(Ptr{T}(ptr.ptr), val)
end

function write_mem
end

function mapped_file
end

function segment_base
end

# Simple session to put some data at specific addresses
immutable FakeMemorySession{T}
    memory_maps::Vector{Tuple{UInt64,Vector{UInt8}}}
    arch
    asid::T
end
Base.show(io::IO, sess::FakeMemorySession) = print(io,"FakeMemorySession(arch=$(sess.arch),asid=$(sess.asid))")
FakeMemorySession(maps,arch) = FakeMemorySession{Void}(maps, arch, nothing)
current_asid(sess::FakeMemorySession) = sess.asid
getarch(sess::FakeMemorySession) = sess.arch

function load{T}(mem::FakeMemorySession, ptr::RemotePtr{T})
    for (addr, data) in mem.memory_maps
        if addr <= UInt64(ptr) <= addr + sizeof(data) - sizeof(T)
            start = (UInt64(ptr)-addr)
            return reinterpret(T,data[1+start:start+sizeof(T)])[]
        end
    end
    error("Address $ptr not found")
end

function load{T}(mem::FakeMemorySession, ptr::RemotePtr{T}, count)
    size = sizeof(T) * count
    for (addr, data) in mem.memory_maps
        if addr <= UInt64(ptr) <= addr + sizeof(data) - size
            start = (UInt64(ptr)-addr)
            return reinterpret(T,data[1+start:start+size])
        end
    end
    error("Address $ptr not found")
end

function store!{T}(session::FakeMemorySession, ptr::RemotePtr{T}, val::T)
    for (addr, data) in mem.memory_maps
        if addr <= UInt64(ptr) <= addr + sizeof(data)
            start = (UInt64(ptr)-addr)
            data[1+start:start+sizeof(T)] = reinterpret(UInt8,[val])
        end
    end
end

function getarch(s::LocalSession)
    sizeof(Ptr{Void}) == 8 ? X86_64.X86_64Arch() : X86_64.X86_32Arch()
end

function getregs
end

function get_thread_area_base
end

function continue!
end

function single_step!
end

function step_until_bkpt!
end

function read_exe
end