text:
# All functions related to dispersal:
# - dispersal
# - establishment
# - viability checking

"""
    disperse!(world, static)

Dispersal of individuals within the world.
"""
function disperse!(world::Array{Patch,1})
    setting("static") && updatesetting("borders", "mainland") #backward compatibility
    for patch in world
        idx = 1
        while idx <= length(patch.seedbank) # disperse each juvenile
            dispmean = patch.seedbank[idx].traits["dispmean"]
            dispshape = patch.seedbank[idx].traits["dispshape"]
            # scaling so that geometric mean follows original distribution (?)
            xdir = rand([-1,1]) * rand(Logistic(dispmean,dispshape))/sqrt(2)
            ydir = rand([-1,1]) * rand(Logistic(dispmean,dispshape))/sqrt(2)
            xdest = patch.location[1] + Int(round(xdir))
            ydest = patch.location[2] + Int(round(ydir))
            target = (xdest, ydest)
            (setting("borders") != "absorbing") && (target = checkborderconditions(world, xdest, ydest))
            possdest = findfirst(x -> x.location==target, world)
            # Note: the following section allows for a mainland/island scenario in which the
            # mainland is static (i.e. mainland individuals are treated as archetypes from
            # an infinitely large source population).
            if setting("static")
                if !world[possdest].isisland
                    idx += 1 # disperse only to islands
                    continue
                end
                if !patch.isisland
                    indleft = deepcopy(patch.seedbank[idx])
                end
            end
            if !setting("static") || patch.isisland
                # only remove individuals from islands / non-static mainland
                indleft = splice!(patch.seedbank,idx)
                idx -= 1
            end
            # Add the dispersed individual to the new community. If no there is no
            # viable target patch, the individual dies.
            (!isnothing(possdest)) && push!(world[possdest].community, indleft)
            idx += 1
        end
    end
end

"""
    establish!(patch)

Establishment of individuals in patch `p`: Sets the adaptation parameters (~fitness)
according to an individual's adaptation to the niches of the surrounding environment.
Also tests whether the individual is viable.
"""
function establish!(patch::Patch)
    temp = patch.temp
    idx = 1
    while idx <= size(patch.community,1)
        if patch.community[idx].marked
            #TODO make this section configurable?
            # (At the moment, tempadaptation and precadaptation aren't actually needed
            # for the Zosterops experiments, so they are needlessly calculated.)
            opt = patch.community[idx].traits["tempopt"]
            tol = patch.community[idx].traits["temptol"]
            fitness = gausscurve(opt, tol, temp, setting("scaleadaptation"))
            fitness > 1 && (fitness = 1) # should be obsolete
            fitness < 0 && (fitness = 0) # should be obsolete
            patch.community[idx].tempadaptation = fitness
            if setting("nniches") >= 2
                opt = patch.community[idx].traits["precopt"]
                tol = patch.community[idx].traits["prectol"]
                fitness = gausscurve(opt, tol, patch.prec, setting("scaleadaptation"))
                fitness > 1 && (fitness = 1) # should be obsolete
                fitness < 0 && (fitness = 0) # should be obsolete
                patch.community[idx].precadaptation = fitness
            end
            if !isviable(patch.community[idx])
                splice!(patch.community,idx)
                continue
            end
            patch.community[idx].marked = false
        end
        idx += 1
    end
end

"""
    establish!(world)

Carry out establishment for each patch in the world.
"""
function establish!(world::Array{Patch,1})
    for patch in world
        (patch.isisland || !setting("static")) && establish!(patch) # pmap(!,patch) ???
    end
end

"""
    isviable(individual)

Check whether this individual conforms to a basic set of constraints
(i.e. all traits are present and certain properties are >= 0).
"""
function isviable(ind::Individual)
    reason = ""
    ind.size <= 0 && ((reason = "size "); @goto die)
    any(collect(values(ind.traits)) .< 0) && ((reason = "traitvalues"); @goto die)
    ind.traits["repsize"] < ind.traits["seedsize"] && ((reason = "seed/rep"); @goto die)
    ind.tempadaptation < 0 && ((reason = "fitness"); @goto die)
    ind.precadaptation < 0 && ((reason = "fitness"); @goto die)
    in("selfing", keys(ind.traits)) && ind.traits["selfing"] > 1 && ((reason = "selfing"); @goto die)
    ind.traits["seqsimilarity"] > 1 && ((reason = "seqsimilarity"); @goto die)
    !traitsexist(ind.traits) && ((reason = "missingtrait"); @goto die)
    return true
    @label die
    simlog("Individual not viable: $reason. Being killed.", 'w')
    return false
end

"""
    checkviability(world)

Check the viability of all individuals. (Note: this does not have to be called
every update, because viability checking is already integrated into `establish!()`.
However, it may be useful to call once after initialisation.)
"""
function checkviability!(world::Array{Patch,1})
    for patch in world
        idx = 1
        while idx <= length(patch.community)
            if !isviable(patch.community[idx])
                splice!(patch.community,idx)
                continue
            end
            idx += 1
        end
    end
end

"""
    traitsexist(traits)

Check a trait dict to make sure it contains the full set of traitnames required
by the model (as defined in the settings).
"""
function traitsexist(traits::Dict{String, Float64})
    missingtraits = setdiff(setting("traitnames"), keys(traits))
    if length(missingtraits) > 0
        simlog("Missing trait $missingtraits. Individual might be killed.", 'w')
        return false
    end
    true
end

"""
    gausscurve(b, c, x, a=1.0)

Calculate the value of the Gauss function ("bell curve") at point x; with
a being the maximum height of the curve, b the position of the curve center and
c the standard deviation ("width").
"""
function gausscurve(b, c, x, a = 1.0)
    if c != 0 && a != 1.0
        a = 1 / (c * sqrt(2 * pi))
        y = a * exp(-(x-b)^2/(2*c^2))
    elseif c != 0
        y = a * exp(-(x-b)^2/(2*c^2))
    else
        y = 0.0
    end
end

"""
    findisland(w)

Within world `w`, find out in which direction from the continent the island lies.
"""
function findisland(world::Array{Patch,1})
    # Context: this function is one of the oldest pieces of code in the model.
    # It was created when GeMM was still targeted at island/mainland studies,
    # and assumes a lot of things that we do differently now. It has been kept
    # for backward compatibility, but is only used when setting("borders") == "mainland".
    xmin, xmax = extrema(map(x->x.location[1],world))
    ymin, ymax = extrema(map(x->x.location[2],world))
    westernborder = filter(x->x.location[1]==xmin,world)
    northernborder = filter(x->x.location[2]==ymax,world)
    easternborder = filter(x->x.location[1]==xmax,world)
    southernborder = filter(x->x.location[2]==ymin,world)
    if all(map(x->x.isisland,westernborder))
        return "west"
    elseif all(map(x->x.isisland,northernborder))
        return "north"
    elseif all(map(x->x.isisland,easternborder))
        return "east"
    elseif all(map(x->x.isisland,southernborder))
        return "south"
    else
        return "none"
    end
end

"""
    checkborderconditions!(w, x, y)

Check if the coordinates `x` and `y` lie within world `w` and correct if not,
considering the user-defined border conditions.
"""
function checkborderconditions(world::Array{Patch,1}, xdest::Int, ydest::Int)
    # First, figure out whether and by how much a destination coordinate overshoots
    # the world borders
    xmin, xmax = extrema(map(x->x.location[1],world))
    ymin, ymax = extrema(map(x->x.location[2],world))
    xrange = xmax  - xmin + 1 # we're counting cells!
    yrange = ymax  - ymin + 1 # we're counting cells!
    xshift = xdest - xmin + 1 # 1-based count of cells
    yshift = ydest - ymin + 1 # 1-based count of cells
    xshift > 0 ? outofx = abs(xshift) : outofx = abs(xshift) + 1
    while outofx > xrange
        outofx -= xrange
    end
    outofx -= 1
    yshift > 0 ? outofy = abs(yshift) : outofy = abs(yshift) + 1
    while outofy > yrange
        outofy -= yrange
    end
    outofy -= 1
    # If a coordinate lies outside the world, correct it depending on the border conditions
    if setting("borders") == "mainland"
        # Ye who come after: see the comment at `findisland()` to
        # understand the history of this dark place...
        islanddirection = findisland(world::Array{Patch,1})
        if islanddirection == "west"
            xdest > xmax && (xdest = xmax - outofx) # east: reflective
            ydest < ymin && (ydest = ymax - outofy) # south: periodic
            ydest > ymax && (ydest = ymin + outofy) # north: periodic
        elseif islanddirection == "north"
            ydest < ymin && (ydest = ymin + outofy) # south: reflective
            xdest < xmin && (xdest = xmax - outofx) # west: periodic
            xdest > xmax && (xdest = xmin + outofx) # east: periodic
        elseif islanddirection == "east"
            xdest < xmin && (xdest = xmin + outofx) # west: reflective
            ydest < ymin && (ydest = ymax - outofy) # south: periodic
            ydest > ymax && (ydest = ymin + outofy) # north: periodic
        elseif islanddirection == "south"
            ydest > ymax && (ydest = ymax - outofy) # north: reflective
            xdest < xmin && (xdest = xmax - outofx) # west: periodic
            xdest > xmax && (xdest = xmin + outofx) # east: periodic
        else
            ydest > ymax && (ydest = ymin + outofy) # north: periodic
            xdest > xmax && (xdest = xmin + outofx) # east: periodic
            ydest < ymin && (ydest = ymax - outofy) # south: periodic
            xdest < xmin && (xdest = xmax - outofx) # west: periodic
        end
    elseif setting("borders") == "reflective"
        ydest > ymax && (ydest = ymax - outofy) # north: reflective
        xdest > xmax && (xdest = xmax - outofx) # east: reflective
        ydest < ymin && (ydest = ymin + outofy) # south: reflective
        xdest < xmin && (xdest = xmin + outofx) # west: reflective
    end
    xdest, ydest
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
function print_id()
    a::Array{Int64} = zeros(10)

    Threads.@threads for i = 1:10
        a[i] = Threads.threadid()
    end
    println(a)
end

function init(n::Int64, x::Vector{Float64}, y::Vector{Float64})
    for i in 1:n
        x[i] = i-1
        y[i] = 1.0
    end
end

function taxpy(n::Int64, a::Float64, x::Vector{Float64}, y::Vector{Float64})
    @inbounds Threads.@threads for i in 1:n
        if x[i] < n/2.
            y[i] = a*x[i] + y[i]
        else
            y[i] = a*x[i] + y[i] + 10.
        end
    end
end

function saxpy(n::Int64, a::Float64, x::Vector{Float64}, y::Vector{Float64})
    @inbounds for i in 1:n
        if x[i] < n/2.
            y[i] = a*x[i] + y[i]
        else
            y[i] = a*x[i] + y[i] + 10.
        end
    end
end

function driver()
    println("Thread count = $(Threads.nthreads())")
    count = 1
    n::Int64 = 1000000
    a::Float64 = 0.5
    y::Array{Float64} = Array{Float64}(n)
    x::Array{Float64} = Array{Float64}(n)

    init(n, x, y)

    # - thread code is not vectorized
    # - as seen in @code_llvm
    # - and that might be the reason it is
    # - slower than the single thread version
    @time for j in 1:count
        #@code_llvm
        taxpy(n, a, x, y)
    end

    println("y = $(y[n])")

    init(n, x, y)

    @time for j in 1:count
        #@code_llvm
        saxpy(n, a, x, y)
    end

    println("y = $(y[n])")
end

#print_id()
driver()

ext: jl


text:
using BioInfoTheory
using Test

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


text:
#using TensorToolbox

export TTtensor, randTTtensor, innerprod, ewprod, full, minus, mtimes, mtimes!, ndims, plus
export reorth, reorth!, size, TTrank, TTsvd, TTtv, norm

mutable struct TTtensor
  cores::CoreCell
  lorth::Bool
  rorth::Bool
  function TTtensor(cores::CoreCell,lorth::Bool,rorth::Bool)
      N=length(cores)
      #[@assert(ndims(cores[n])==3,"Core tensors must be tensors of order 3!") for n=1:N]
      @assert(size(cores[1],1)==size(cores[N],3)==1,"Core tensors are of incorrect sizes.")
      [@assert(size(cores[n],3)==size(cores[n+1],1),"Dimension mismatch.") for n=1:N-1]
      for G in cores[1:end-1]
          Gn=tenmat(G,row=1:2)
          if norm(eye(size(Gn,2))- transpose(Gn)*Gn)>(size(Gn,1)^2)*eps()
              lorth=false
              break
          end
      end
      for G in cores[2:end]
          Gn=tenmat(G,1)
          if norm(eye(size(Gn,1))- Gn*transpose(Gn))>(size(Gn,2)^2)*eps()
              rorth=false
              break
          end
      end
      new(cores,lorth,rorth)
  end
end
TTtensor(cores::CoreCell)=TTtensor(cores,true,true)
TTtensor(cores::Array{Array{T,3},1}) where {T<:Number}=TTtensor(CoreCell(cores),true,true)
TTtensor(cores::Array{AbstractArray{T,N},1}) where {T<:Number,N}=TTtensor(CoreCell(cores),true,true)

"""
    randTTtensor(I::Vector,R::Vector)
    randTTtensor(I::Integer,R::Integer,N::Integer)

Create random TTtensor of size I and TT-rank rank R, or of order N and size I × ⋯ × I and TT-rank (R,...,R).
"""
function randTTtensor(I::Vector{<:Integer},R::Vector{<:Integer})
  N=length(I)
  @assert(length(R)==N-1,"Dimension mismatch.")
  cores=CoreCell(undef,N) #create radnom cores
  Rfull=[1;R;1]
  [cores[n]=randn(Rfull[n],I[n],Rfull[n+1]) for n=1:N]
  TTtensor(cores)
end
randTTtensor(I::Integer,R::Integer,N::Integer)=randTTtensor(repmat([I],N),repmat([R],N-1));
#For input defined as tuples or nx1 matrices - ranttensor(([I,I,I],[R,R,R]))
function randTTtensor(arg...)
  randTTtensor([arg[1]...],[arg[2]...])
end

"""
---------------------------------------------------
    contract(T::TTtensor,X::Array,start_mode,modes)

Contract modes 1:ndims(X-1) of full tensor X to TTtensor T cores[start_mode:start_mode+modes-1].
Output: TTtensor.
"""
function contract(T::TTtensor,X::Array,start_mode=1,modes=0)
    if modes==0
        modes=ndims(X)-1
    end
    @assert(ndims(X)==modes+1,"Dimension mismatch.")
    for n=start_mode:start_mode+modes-1
        if n==start_mode
            X=contract(T.cores[n],2,X,1)
            X=permutedims(X,[collect(2:ndims(X))...;1])
        else
            X=contract(T.cores[n],collect(1:2),X,collect(1:2))
        end
    end
    X=permutedims(X,[3,2,1])
    N=ndims(T)-modes+1
    G=CoreCell(undef,N)
    for n=1:start_mode-1
        G[n]=T.cores[n]
    end
    G[start_mode]=X
    for n=start_mode+1:N
         G[n]=T.cores[modes+n-1]
    end
    TTtensor(G)
end
contract(X::Array,T::TTtensor,start_mode=1,modes=0)=contract(T,X,start_mode,modes)
function contract(X::CoreCell,squeeze=true)
    N=length(X)
    if N==1
        return X[1]
    end
    @assert(any([size(X[n])[end]==size(X[n+1])[1] for n=1:N-1]),"Dimensions mismatch.")
    Xcontr=contract(X[1],X[2])
    for n=3:N
        Xcontr=contract(Xcontr,X[n])
    end
    squeeze == true ? dropdims(Xcontr) : Xcontr
end
#
function contract(X1::CoreCell,X2::CoreCell)
    N=length(X2)
    @assert(length(X1)==N,"Dimensions mismatch.")
    @assert(size(X1)[1:end-1]==size(X2)[1:end-1],"Dimensions mismatch.")
    Xcontr=transpose(dropdims(contract(X1[1],2,X2[1],2)))
    for n=2:N
        Xcontr=contract(X1[n],1,Xcontr,2,[3,1,2])
        Xcontr=contract(X2[n],[1,2],Xcontr,[1,2])
    end
    Xcontr
end

function innerprod(X1::TTtensor,X2::TTtensor)
    Isz=size(X1)
    @assert(Isz==size(X2),"Dimension mismatch.")
    v=kron(X1.cores[1][:,1,:],X2.cores[1][:,1,:])
    for i=2:Isz[1]
        v+=kron(X1.cores[1][:,i,:],X2.cores[1][:,i,:])
    end
    for n=2:ndims(X1)
        w=VectorCell(undef,Isz[n])
        for i=1:Isz[n]
            w[i]=krontv(transpose(X1.cores[n][:,i,:]),transpose(X2.cores[n][:,i,:]),vec(v))
        end
        v=vec(sum(w))
    end
    v[1]
end

function innerprod(X1::CoreCell,X2::CoreCell)
    Isz=[size(X,2) for X in X1]
    @assert(Isz==[size(X,2) for X in X2],"Dimension mismatch.")
    v=kron(X1[1][:,1,:],X2[1][:,1,:])
    for i=2:Isz[1]
        v+=kron(X1[1][:,i,:],X2[1][:,i,:])
    end
    for n=2:ndims(X1)
        w=VectorCell(undef,Isz[n])
        for i=1:Isz[n]
            w[i]=krontv(transpose(X1[n][:,i,:]),transpose(X2[n][:,i,:]),vec(v))
        end
        v=vec(sum(w))
    end
    v[1]
end

function ewprod(X1::TTtensor,X2::TTtensor)
	Isz=size(X1)
	@assert(Isz == size(X2),"Dimension mismatch.")
	N=ndims(X1)
  	cores=CoreCell(undef,N) #initilize cores
	R=[1,TTrank(X1).*TTrank(X2)...,1]
  	for n=1:N
    	cores[n]=zeros(R[n],Isz[n],R[n+1])
    	for i=1:Isz[n]
        	cores[n][:,i,:]=kron(X1.cores[n][:,i,:],X2.cores[n][:,i,:]);
    	end
	end
  	TTtensor(cores)
end

function full(X::TTtensor)
    contract(X.cores,true)
end

"""
    minus(X::TTtensor,Y::TTtensor)

Subtraction of two TTtensors. Same as: X-Y.
"""
function minus(X1::TTtensor,X2::TTtensor)
	1*X1+(-1)*X2
end
-(X1::TTtensor,X2::TTtensor)=minus(X1,X2)

"""
    mtimes(a,X)

Scalar times TTtensor. Same as: a*X.
"""
function mtimes(α::Number,X::TTtensor)
    # G=deepcopy(X.cores)
    # G[1]=α*G[1]
	# TTtensor(G)
    TTtensor(mtimes(α,X.cores))
end
*(α::Number,X::TTtensor)=mtimes(α,X)
*(X::TTtensor,α::Number)=*(α,X)
function mtimes!(α::Number,X::TTtensor)
    X.cores[1]=α*X.cores[1]
end
function mtimes(α::Number,G::CoreCell)
    Gcopy=deepcopy(G)
    Gcopy[1]=α*G[1]
	Gcopy
end
*(α::T,G::CoreCell) where {T<:Number} = mtimes(α,G)
*(G::CoreCell,α::T) where {T<:Number} = *(α,G)
function mtimes!(α::Number,G::CoreCell)
    G[1]=α*G[1]
end

function ndims(X::TTtensor)
	length(X.cores)
end

#Frobenius norm of a TTtensor. **Documentation in Base.
function norm(X::TTtensor,orth="")
    if X.lorth==true
        sqrt(contract(X.cores[end],[1,2,3],X.cores[end],[1,2,3])[1])
    elseif X.rorth==true
        sqrt(contract(X.cores[1],[1,2,3],X.cores[1],[1,2,3])[1])
    else
        if orth=="reorth"
            Xorth=reorth(X)
            sqrt(abs.(innerprod(Xorth,Xorth)))
        else
            sqrt(abs.(innerprod(X,X)))
        end
    end
end

"""
    plus(X::TTtensor,Y::TTtensor)

Addition of two TTtensors. Same as: X+Y.
"""
function plus(X1::TTtensor,X2::TTtensor)
	@assert(size(X1) == size(X2),"Dimension mismatch.")
    TTtensor(plus(X1.cores,X2.cores))
end
+(X1::TTtensor,X2::TTtensor)=plus(X1,X2)

function plus(X1::CoreCell,X2::CoreCell)
    N=length(X1)
    @assert(N==length(X2),"Dimension mismatch.")
    Xres=CoreCell(undef,N)
    sz=size(X1[1]).+(0,0,size(X2[1])[3])
    Xres[1]=reshape([tenmat(X1[1],1) tenmat(X2[1],1)],sz)
    sz=size(X1[N]).+(size(X2[N])[1],0,0)
    Xres[N]=reshape([X1[N][:,:,1]; X2[N][:,:,1]],sz)
    for n=2:N-1
        sz1=size(X1[n])
        sz2=size(X2[n])
        sz=sz1.+(sz2[1],0,sz2[3])
        C=zeros(sz[1],sz[2]*sz[3])
        C[1:sz1[1],1:sz1[2]*sz1[3]]=tenmat(X1[n],1)
        C[sz1[1]+1:end,sz1[2]*sz1[3]+1:end]=tenmat(X2[n],1);
        Xres[n]=reshape(C,sz)
    end
	Xres
end
+(X1::CoreCell,X2::CoreCell)=plus(X1,X2)

function reorth(X::TTtensor,direction="left",full="no")
    TTtensor(reorth(X.cores,direction,full))
end
function reorth(X::CoreCell,direction="left",full="no")
    G=deepcopy(X)
    N=length(X)
    if direction=="left"
        for n=1:N-1
            Q,Rt=qr(tenmat(G[n],row=1:2))
            Q=Matrix(Q)
            sz=(size(G[n])[1:2]...,size(Q)[2])
            G[n]=reshape(Q,sz)
            G[n+1]=ttm(G[n+1],Rt,1)
        end
        if full=="full"
            # Q,Rt=qr(tenmat(G[N],row=1:2))
            # Q=Matrix(Q)
            # sz=(size(G[N])[1:2]...,size(Q)[2])
            # G[N]=reshape(Q,sz)
            G[N]=G[N]/norm(G[N])
        end
    elseif direction == "right"
        for n=N:-1:2
            Q,Rt=qr(transpose(tenmat(G[n],1)))
            Q=Matrix(Q)
			Rnorm=norm(Rt)
            sz=(size(Q)[2],size(G[n])[2:3]...)
            G[n]=reshape(Matrix(transpose(Q)),sz)
            G[n-1]=ttm(G[n-1],Rt,3)
        end
        if full=="full"
            # Q,Rt=qr(transpose(tenmat(G[1],1)))
            # Q=Matrix(Q)
            # sz=(size(Q)[2],size(G[1])[2:3]...)
            # G[1]=reshape(Matrix(transpose(Q)),sz)
            G[1]=G[1]/norm(G[1])
        end
    else
        error("Invalid direction. Should be either \"left\" or \"right\".")
    end
    G
end

function reorth!(X::TTtensor,direction="left",full="no")
    N=ndims(X)
    if direction=="left"
        for n=1:N-1
            Q,Rt=qr(tenmat(X.cores[n],row=1:2))
            Q=Matrix(Q)
            sz=(size(X.cores[n])[1:2]...,size(Q)[2])
            X.cores[n]=reshape(Q,sz)
            X.cores[n+1]=ttm(X.cores[n+1],Rt,1)
        end
        if full=="full"
            # Q,Rt=Matrix(qr(tenmat(X.cores[N],row=1:2)).Q)
            # sz=(size(X.cores[N])[1:2]...,size(Q)[2])
            # X.cores[N]=reshape(Q,sz)
            X.cores[N]=X.cores[N]/norm(X.cores[N])
        end
        X.lorth=true
    elseif direction == "right"
        for n=N:-1:2
            Q,Rt=qr(transpose(tenmat(X.cores[n],1)))
            Q=Matrix(Q)
            sz=(size(Q)[2],size(X.cores[n])[2:3]...)
            X.cores[n]=reshape(Matrix(transpose(Q)),sz)
            X.cores[n-1]=ttm(X.cores[n-1],Rt,3)
        end
        if full=="full"
            # Q=Matrix(qr(transpose(tenmat(X.cores[1],1))).Q)
            # sz=(size(Q)[2],size(X.cores[1])[2:3]...)
            # X.cores[1]=reshape(Matrix(transpose(Q)),sz)
            X.cores[1]=X.cores[1]/norm(X.cores[1])
        end
        X.rorth=true
    else
        error("Invalid direction. Should be either \"left\" or \"right\".")
    end
end
function reorth!(X::CoreCell,direction="left",full="no")
    N=length(X)
    if direction=="left"
        for n=1:N-1
            Q,Rt=qr(tenmat(X[n],row=1:2))
            Q=Matrix(Q)
            sz=(size(X[n])[1:2]...,size(Q)[2])
            X[n]=reshape(Q,sz)
            X[n+1]=ttm(X[n+1],Rt,1)
        end
        if full=="full"
            # Q=Matrix(qr(tenmat(X[N],row=1:2)).Q)
            # sz=(size(X[N])[1:2]...,size(Q)[2])
            # X[N]=reshape(Q,sz)
            X[N]=X[N]/norm(X[N])
        end
    elseif direction == "right"
        for n=N:-1:2
            Q,Rt=qr(tenmat(X[n],1)')
            Q=Matrix(Q)
            sz=(size(Q)[2],size(X[n])[2:3]...)
            X[n]=reshape(Q',sz)
            X[n-1]=ttm(X[n-1],Rt,3)
        end
        if full=="full"
            # Q=Matrix(qr(tenmat(X[1],1)').Q)
            # sz=(size(Q)[2],size(X[1])[2:3]...)
            # X[1]=reshape(Q',sz)
            X[1]=X[1]/norm(X[1])
        end
    else
        error("Invalid direction. Should be either \"left\" or \"right\".")
    end
end

#Size of a TTtensor. **Documentation in Base.
function size(X::TTtensor)
    tuple([size(X.cores[n],2) for n=1:ndims(X)]...)
end
function size(X::TTtensor,mode::Int)
    size(X.cores[mode],2)
end

"""
    TTrank(X::TTtensor,full=false)
    TTrank(X::CoreCell,full=false)

Representational TT-rank of X. If cores of X are of size R{n-1} x In x Rn, returns (R1,...,R{n-1}) for full=false or (1,R1,....,Rn) for full=true."""
function TTrank(X::TTtensor,full=false)
    TTrank(X.cores,full)
end
function TTrank(X::CoreCell,full=false)
    rnk=[size(X[n])[3] for n=1:length(X)-1]
    if full==true
        rnk=[1,rnk...,1]
    end
    tuple(rnk...)
end

"""
    TTsvd(X;tol=1e-8,reqrank=[])
    TTsvd(X::TTtensor;tol=1e-8,reqrank=[])

Decomposes full tensor X into its TT format. If input is already a TTtensor, recompresses it."""
function TTsvd(X::AbstractArray{<:Number,N};tol=1e-8,reqrank=[],tol_rel=0) where N
    Isz=size(X)
    if length(reqrank)>0
        @assert(length(reqrank)==N-1,"Dimension mismatch. For a N-dim tensor, reqrank has to be a vector of length N-1.")
        rnk=1
    else
        reqrank=[Isz...]
        rnk=0
		tol=tol/sqrt(N-1)
    end
    G=CoreCell(undef,N)
    Xtmp=copy(X)
    for n=1:N-1
        n==1 ?  R=1 : R=reqrank[n-1]
        m=R*Isz[n]
        Xtmp=reshape(Xtmp,m,:)
        U,S,V=svd(Xtmp)
        if rnk==1
            U=U[:,1:reqrank[n]]
            S=S[1:reqrank[n]]
            V=V[:,1:reqrank[n]]
        else
            if tol_rel != 0
                tol=tol_rel*S[1]
            end
            #K=findall(x-> x>tol ? true : false,S)
			Sum=0
			k=0
			while sqrt(Sum)<=tol && length(S)-k>0
				Sum=Sum+S[end-k]^2
				k+=1
			end
			K=length(S)-k+1
            U=U[:,1:K]
            S=S[1:K]
            V=V[:,1:K]
            reqrank[n]=length(S)
        end
        G[n]=reshape(U,(R,Isz[n],reqrank[n]))
        Xtmp=Diagonal(S)*transpose(V)
    end
    G[N]=reshape(Xtmp,(size(Xtmp,1),size(Xtmp,2),1))
    TTtensor(G)
end

function TTsvd(X::TTtensor;reqrank=[],tol=1e-8,tol_rel=0)
    TTtensor(TTsvd(X.cores,reqrank=reqrank,tol=tol,tol_rel=tol_rel))
end

function TTsvd(X::CoreCell;reqrank=[],tol=1e-8,tol_rel=0)
    N=length(X)
	tol=tol/sqrt(N-1)
	Isz=[size(X[n],2) for n=1:N]
 	G=reorth(X,"right")
    if length(reqrank)>0
        @assert(length(reqrank)==N-1,"Requested rank of inapropriate size.")
        if any(reqrank.<TTrank(G))==false
            @warn "Requested rank larger than the current rank."
            return X
        end
        for n=1:N-1
            U,S,V=svd(tenmat(G[n],row=1:2))
            if reqrank[n]<=size(U,2)
                U=U[:,1:reqrank[n]]
                V=V[:,1:reqrank[n]]*Diagonal(S[1:reqrank[n]])
			else
				@warn "Requested rank larger than the actual rank."
				reqrank[n]=size(U,2)
            end
            n==1 ? sz=(1,Isz[1],reqrank[1]) : sz=(reqrank[n-1],Isz[n],reqrank[n])
            G[n]=reshape(U,sz)
            G[n+1]=ttm(G[n+1],transpose(V),1)
        end
    else
        for n=1:N-1
            U,S,V=svd(tenmat(G[n],row=1:2))
            if tol_rel != 0
                tol=tol_rel*S[1]/sqrt(N-1)
            end
            Sum=0
			k=0
			while sqrt(Sum)<tol && length(S)-k>0
				Sum=Sum+S[end-k]^2
				k+=1
			end
			#K=findall(x-> x>δ ? true : false,S)
			if k>0
				K=length(S)-k+1
        		U=U[:,1:K]
        		V=V[:,1:K]
        		S=S[1:K]
			end
			#@show tol, norm(tenmat(G[n],row=1:2)-U*Diagonal(S)*transpose(V))
            G[n]=reshape(U,(size(G[n],1),Isz[n],:))
            G[n+1]=ttm(G[n+1],V*Diagonal(S),1,'t')
        end
    end
    G
end

"""
    TTtv(X::TTtensor,u::VectorCell,mode::Int)
    TTtv(X::TTtensor,U::AbstractMatrix{<:Number},mode::Int)
    TTtv(X::CoreCell,u::VectorCell)
    TTtv(X::CoreCell,U::AbstractMatrix{<:Number})

TT-tensor times vector. For N=ndims(X) and n=[mode,mode+1,...,N] contracts nth core
of the TT-tensor with vector u[N-n+1]. If input is a matrix, then vectors u are its columns.
Called on a subset of cores, it contracts all given cores to vectors.
"""
function TTtv(X::TTtensor,u::VectorCell,mode::Int)
    if mode==1
        TTtv(X.cores,u)
    else
        G=CoreCell(undef,mode-1)
        G[1:mode-1]=deepcopy(X.cores[1:mode-1])
        v=TTtv(X.cores[mode:end],u)
        G[mode-1]=contract(G[mode-1],v)
        mode==2 ? dropdims(G[1]) : TTtensor(G)
    end
end
function TTtv(X::TTtensor,U::AbstractMatrix{<:Number},mode::Int)
    L=size(U,2)
    u=VectorCell(undef,L)
    for l=1:L
        u[l]=U[:,l]
    end
    TTtrv(X,U,mode)
end
function TTtv(X::CoreCell,u::VectorCell)
    N=length(X)
    @assert(length(u)==N,"Dimension mismatch.")
    G=copy(X[N])
    for n=N:-1:1
        G=contract(G,2,u[n],1)
        if n!=1
            G=contract(X[n-1],G)
        end
    end
    G
end
function TTtv(X::CoreCell,U::AbstractMatrix{<:Number})
    Ucell=MatrixCell(undef,1)
    Ucell[1]=U
    TTtrv(X,Ucell)
end


text:
# Rose Teague 06/07/2018
# UpdateEnergyTest.jl
#
# Tests the function 'SiteEnergyFromDipoles' using eigenstate of the system.
push!(LOAD_PATH,"../src/")

using TheDancer
using Base.Test

N, Edisorder, Jdisorder, modelJ, B, dipolestrength, r = init!(3, 0.0, 0.0, 0.0, 298, 8)

"""
Test the on-site energy is correctly calculated as the potential energy felt by
a test charge on site i, in the presence of the electric fields established by
all other sites.

Testing in a symmetric 3-site environment, energy at the central position should
be 0 as effects from dipoles on either side cancel each other out.

External Field should be set to zero.
"""

function main()
    dipoles = [1.0,1.0,1.0]
    E = [0.0,0.0]
    S = [0.0,0.0,0.0]
    S,H = SiteEnergyFromDipoles(dipoles,S,E)

    S_check = 5/(4*r^2)

    ε = 1e-2
    @test S_check ≈ S[1] atol = ε
    @test 0.0 ≈ S[2] atol = ε
end

main()


text:
"""
```julia
	latp_algorithm(G:SimpleWeightedGraph,s::Int)
```
	Perform the *Least-Action trip planning* algorithm on G, starting from s with alpha parameter
"""
latp_algorithm(G::SimpleWeightedGraph,s::Int,alpha::Int)

function latp_algorithm(G::SimpleWeightedGraph,s::Int,alpha::Int,T::DateTime)
	V = Set(vertices(G))
	Random.seed!(time_ns())
	Visited = [s] #It's a vector because we need to maintain visiting order
	Visited_time = [T]
	c = s

	next = 0
	while length(setdiff(V,Visited)) != 0
		#Firstly, get all unvisited nodes next to current node and order them in closeness
		neighbors_ordered = sort([ (G.weights[c,i],i) for i in setdiff(neighbors(G,c),Visited)])
		#Then the probability vector is calculated
		prob_vector = map( d->(1/d[1])^alpha, neighbors_ordered)
		prob_sum = sum( prob_vector )
		prob_vector = map( d -> d / prob_sum, prob_vector)
		#Find node: the probability of selecting a node i is the sum of the probability of all preceiding probs
		p = rand()
		for i in 1:length(neighbors_ordered)
			if p <= sum(prob_vector[1:i])
				next = neighbors_ordered[i][2]
				break
			end
		end
		#TODO: Add pause time
		push!(Visited,next)
		push!(Visited_time,Visited_time[length(Visited_time)] + Dates.Millisecond(10)) #10 Millisecond, just a test
		c = next
	end
	return [ (Visited[i],Visited_time[i]) for i in 1:length(Visited) ] #Combine visited nodes with timestamp
end


text:
"""
    Rule

A `Rule` object contains the information required to apply an
`applyrule` method to every cell of every timestep of a simulation.

Rules are applied to the grid using the [`applyrule`](@ref) method:

```julia
@inline applyrule(data::SimData, rule::YourRule, state, index) = ...
```

Where `index` is a `Tuple` of `Int`, and `state` is a single value, or a `NamedTuple`
if multiple grids are requested. the [`SimData`](@ref) object can be used to access current
timestep and other simulation data and metadata.

Rules can be updated from the original rule before each timestep, in [`modifyrule`](@ref):

```julia
modifyrule(rule::YourRule, data::SimData) = ...
```

Rules can also be run in sequence, often wrapped in a `Tuple` or [`Ruleset`](@ref)s.

DynamicGrids guarantees that:

- `modifyrule` is run once for every rule for every timestep.
    The result is passed to `applyrule`, but not retained after that.
- `applyrule` is run once for every rule, for every cell, for every timestep, unless an
    optimisation like `SparseOpt` is enable to skips empty cells.
- the output of running a rule for any cell does not affect the input of the
    same rule running anywhere else in the grid.
- rules later in the sequence are passed grid state updated by the earlier rules.
- masked areas and wrapped or removed boundary regions are updated between all rules and timesteps

## Multiple grids

The `NamedTuple` keys will match the keys in `R`, which is a type like `Tuple{:key1,:key1}`.
Note the names are user-specified, and should never be fixed by a `Rule`.

They can be retrieved from the type here as `A` and `B` :

```julia
applyrule(data::SimData, rule::YourCellRule{Tuple{R1,R2},Tuple{W1,W2}}, state, index) where {R1,R2,W1,W2}
```

By default the output is written to the current cell in the specified `W` write grid/s.
`Rule`s writing to multiple grids, simply return a `Tuple` in the order specified by
the `W` type params.

## Rule Performance

Rules may run many millions of times during a simulation. They need to be fast.
Some basic guidlines for writing rules are:

- Never allocate memory in a `Rule` if you can help it.
- Type stability is essential. [`isinferred`](@ref) is useful to check
    if your rule is type-stable.
- Using the `@inline` macro on `applyrule` can help force inlining your
    code into the simulation.
- Reading and writing from multiple grids is expensive due to additional load
    on fast cahce memory. Try to limit the number of grids you use.
- Use a graphical profiler, like ProfileView.jl, to check your rules overall
    performance when run with `sim!`.

"""
abstract type Rule{R,W} end

# Rules are also traits - because we use wrapper rules
ruletype(::Rule) = Rule

#= Default constructors for all Rules.
Sets both the read and write grids to `:_default`.

This strategy relies on a one-to-one relationship between fields
and type parameters, besides the initial `R` and `W` params.  =#

# No {R,W} with args or kw
function (::Type{T})(args...; kw...) where T<:Rule
    T{:_default_,:_default_}(args...; kw...)
end
# {R,W} with args
function (::Type{T})(args...) where T<:Rule{R,W} where {R,W}
    # _checkfields(T, args)
    T{map(typeof, args)...}(args...)
end
# Only R specified
function (::Type{T})(args...; kw...) where T<:Rule{R} where R
    T{R}(args...; kw...)
end

@generated function Base.keys(rule::Rule{R,W}) where {R,W}
    Expr(:tuple, QuoteNode.(union(_asiterable(W), _asiterable(R)))...)
end

@inline _writekeys(::Rule{R,W}) where {R,W} = W
@generated function _writekeys(::Rule{R,W}) where {R,W<:Tuple}
    Expr(:tuple, QuoteNode.(W.parameters)...)
end

@inline _readkeys(::Rule{R,W}) where {R,W} = R
@generated function _readkeys(::Rule{R,W}) where {R<:Tuple,W}
    Expr(:tuple, QuoteNode.(R.parameters)...)
end

# Define the constructor for generic rule reconstruction in Flatten.jl and Setfield.jl
function ConstructionBase.constructorof(::Type{T}) where T<:Rule{R,W} where {R,W}
    T.name.wrapper{R,W}
end

# Find the largest radius present in the passed in rules.
function radius(rules::Tuple{Vararg{<:Rule}})
    allkeys = Tuple(union(map(keys, rules)...))
    maxradii = Tuple(radius(rules, key) for key in allkeys)
    return NamedTuple{allkeys}(maxradii)
end
radius(rules::Tuple{}) = NamedTuple{(),Tuple{}}(())
# Get radius of specific key from all rules
radius(rules::Tuple{Vararg{<:Rule}}, key::Symbol) =
    reduce(max, radius(ru) for ru in rules if key in keys(ru); init=0)

radius(rule::Rule, args...) = 0

"""
    RuleWrapper <: Rule

A `Rule` that wraps other rules, altering their behaviour or how they are run.
"""
abstract type RuleWrapper{R,W} <: Rule{R,W} end

"""
    Cellrule <: Rule

A `Rule` that only writes and uses a state from single cell of the read grids,
and has its return value written back to the same cell(s).

This limitation can be useful for performance optimisation,
such as wrapping rules in [`Chain`](@ref) so that no writes occur between rules.

`CellRule` is defined with :

```julia
struct YourCellRule{R,W} <: CellRule{R,W} end
```

And applied as:

```julia
function applyrule(data::SimData, rule::YourCellRule{R,W}, state, index) where {R,W}
    state * 2
end
```

As the `index` is provided in `applyrule`, you can look up an [`aux`](@ref) array
using `aux(data, Val{:auxname}())[index...]` to access cell-specific variables for
your rule.

It's good to add a struct field to hold the `Val{:auxname}()` object instead of
using names directly, so that users can set the aux name themselves to suit the
scripting context.

"""
abstract type CellRule{R,W} <: Rule{R,W} end

ruletype(::CellRule) = CellRule

"""
    SetRule <: Rule

Abstract supertype for rules that manually write to the grid in some way.

These must define methods of [`applyrule!`](@ref).
"""
abstract type SetRule{R,W} <: Rule{R,W} end

ruletype(::SetRule) = SetRule

"""
    SetCellRule <: Rule

Abstract supertype for rules that can manually write to any cells of the
grid that they need to.

`SetCellRule` is applied with a method like:

```julia
function applyrule!(data::SimData, rule::YourSetCellRule, state, index)
     inc = 1
     add!(data, inc, index...)
     return nothing
end
```

Note the `!` bang - this method alters the state of `data`.

To update the grid, you can use atomic operators [`add!`](@ref), [`sub!`](@ref),
[`min!`](@ref), [`max!`](@ref), and [`and!`](@ref), [`or!`](@ref) for `Bool`.
These methods safely combined writes from all grid cells - directly using `setindex!`
would cause bugs.

It there are multiple write grids, you will need to get the grid keys from
type parameters, here `W1` and `W2`:

```julia
function applyrule(data, rule::YourSetCellRule{R,Tuple{W1,W2}}, state, index) where {R,W1,W2}
     inc = 1
     add!(data[W1], inc, index...)
     add!(data[W2], 2inc, index...)
     return nothing
end
```

DynamicGrids guarantees that:

- values written to anywhere on the grid do not affect other cells in
    the same rule at the same timestep.
- values written to anywhere on the grid are available to the next rule in the
    sequence, or the next timestep.
- if atomic operators are always used, race conditions will not occur on any hardware.

"""
abstract type SetCellRule{R,W} <: SetRule{R,W} end

ruletype(::SetCellRule) = SetCellRule

"""
    NeighborhoodRule <: Rule

A Rule that only accesses a neighborhood centered around the current cell.
`NeighborhoodRule` is applied with the method:

```julia
applyrule(data::SimData, rule::YourNeighborhoodRule, state, I)
```

`NeighborhoodRule` must have a `neighborhood` field, that holds
a [`Neighborhood`](@ref) object. `neighbors(rule)` returns an iterator
over the surrounding cell pattern defined by the `Neighborhood`.

For each cell in the grids the neighborhood buffer will be updated
for use in the `applyrule` method, managed to minimise array reads.

This allows memory optimisations and the use of high perforance routines on the
neighborhood buffer. It also means that and no bounds checking is required in
neighborhood code.

For neighborhood rules with multiple read grids, the first is always
the one used for the neighborhood, the others are passed in as additional
state for the cell. Any grids can be written to, but only for the current cell.
"""
abstract type NeighborhoodRule{R,W} <: Rule{R,W} end

ruletype(::NeighborhoodRule) = NeighborhoodRule

neighbors(rule::NeighborhoodRule) = neighbors(neighborhood(rule))
neighborhood(rule::NeighborhoodRule) = rule.neighborhood
offsets(rule::NeighborhoodRule) = offsets(neighborhood(rule))
kernel(rule::NeighborhoodRule) = kernel(neighborhood(rule))
positions(rule::NeighborhoodRule, args...) = positions(neighborhood(rule), args...)
neighborhoodkey(rule::NeighborhoodRule{R,W}) where {R,W} = R
# The first argument is for the neighborhood grid
neighborhoodkey(rule::NeighborhoodRule{<:Tuple{R1,Vararg},W}) where {R1,W} = R1
_buffer(rule::NeighborhoodRule) = _buffer(neighborhood(rule))
@inline _setbuffer(rule::NeighborhoodRule, _buffer) =
    @set rule.neighborhood = _setbuffer(rule.neighborhood, _buffer)
radius(rule::NeighborhoodRule, args...) = radius(neighborhood(rule))

"""
    SetNeighborhoodRule <: SetRule

A Rule that only writes to its neighborhood.

[`positions`](@ref) and [`offsets`](@ref) are useful iterators for modifying
neighborhood values. Atomic

`SetNeighborhood` rules must return a `Neighborhood` object from `neighborhood(rule)`.
By default this is `rule.neighborhood`. If this property exists, no interface methods
are required.
"""
abstract type SetNeighborhoodRule{R,W} <: SetRule{R,W} end

ruletype(::SetNeighborhoodRule) = SetNeighborhoodRule

neighborhood(rule::SetNeighborhoodRule) = rule.neighborhood
offsets(rule::SetNeighborhoodRule) = offsets(neighborhood(rule))
kernel(rule::SetNeighborhoodRule) = kernel(neighborhood(rule))
positions(rule::SetNeighborhoodRule, args...) = positions(neighborhood(rule), args...)
radius(rule::SetNeighborhoodRule, args...) = radius(neighborhood(rule))
neighborhoodkey(rule::SetNeighborhoodRule{R,W}) where {R,W} = R
neighborhoodkey(rule::SetNeighborhoodRule{<:Tuple{R1,Vararg},W}) where {R1,W} = R1

"""
    SetGridRule <: Rule

A `Rule` applies to whole grids. This is used for operations that don't benefit from
having neighborhood buffering or looping over the grid handled for them, or any specific
optimisations. Best suited to simple functions like `rand`(write)` or using convolutions
from other packages like DSP.jl. They may also be useful for doing other custom things that
don't fit into the DynamicGrids.jl framework during the simulation.

Grid rules specify the grids they want and are sequenced just like any other grid.

```julia
struct YourSetGridRule{R,W} <: SetGridRule{R,W} end
```

And applied as:

```julia
function applyrule!(data::SimData, rule::YourSetGridRule{R,W}) where {R,W}
    rand!(data[W])
end
```
"""
abstract type SetGridRule{R,W} <: Rule{R,W} end

ruletype(::SetGridRule) = SetGridRule

"""
    SetGrid{R,W}(f)

Apply a function `f` to fill whole grid/s.

```julia
rule = SetGrid{:a,:b}() do a, b
    b .= a
end
```
"""
struct SetGrid{R,W,F} <: SetGridRule{R,W}
    "Function to apply to the read values"
    f::F
end
SetGrid{R,W}(; kw...) where {R,W} = _nofunctionerror(SetGrid)

@inline function applyrule!(data, rule::SetGrid{R,W}) where {R,W}
    rule.f(
        map(r -> source(getindex(data, r)), _asiterable(R))...,
        map(w -> source(getindex(data, w)), _asiterable(W))...,
    )
end

"""
    Call <: CellRule

    Cell(f)
    Cell{R,W}(f)

A [`CellRule`](@ref) that applies a function `f` to the
`read` grid cells and returns the `write` cells.

Especially convenient with `do` notation.

## Example

Set the cells of grid `:c` to the sum of `:a` and `:b`:

```julia
simplerule = Cell() do a, b
    a + b
end
```

If you need to use multiple grids (a and b), use the `read`
and `write` arguments. If you want to use external variables,
wrap the whole thing in a `let` block, for performance.

```julia
rule = let y = y
    rule = Cell{Tuple{:a,:b},:b}() do a, b
        a + b * y
    end
end
```
"""
struct Cell{R,W,F} <: CellRule{R,W}
    "Function to apply to the read values"
    f::F
end
Cell{R,W}(; kw...) where {R,W} = _nofunctionerror(Cell)

@inline function applyrule(data, rule::Cell, read, I)
    let rule=rule, read=read
        rule.f(read)
    end
end

"""
    Neighbors <: NeighborhoodRule

    Neighbors(f, neighborhood=Moor(1))
    Neighbors{R,W}(f, neighborhood=Moore())

A [`NeighborhoodRule`](@ref) that receives a neighbors object for the first
`read` grid and the passed in neighborhood, followed by the cell values for
the required grids, as with [`Cell`](@ref).

Returned value(s) are written to the `write`/`W` grid.

As with all [`NeighborhoodRule`](@ref), you do not have to check bounds at
grid edges, that is handled for you internally.

Using [`SparseOpt`](@ref) may improve neighborhood performance
when zero values are common and can be safely ignored.

## Example

```julia
rule = let x = 10
    Neighbors{Tuple{:a,:b},:b}() do hood, a, b
        data[:b][I...] = a + b^x
    end
end
```

The `let` block may improve performance.
"""
struct Neighbors{R,W,F,N} <: NeighborhoodRule{R,W}
    "Function to apply to the neighborhood and read values"
    f::F
    "Defines the neighborhood of cells around the central cell"
    neighborhood::N
end
Neighbors{R,W}(; kw...) where {R,W} = _nofunctionerror(Neighbors)
Neighbors{R,W}(f; neighborhood=Moore(1)) where {R,W} =
    Neighbors{R,W}(f, neighborhood)

@inline function applyrule(data, rule::Neighbors, read, I)
    let rule=rule, hood=neighborhood(rule), read=read
        rule.f(hood, read)
    end
end

"""
    SetCell <: SetCellRule

    SetCell(f)
    SetCell{R,W}(f)

A [`SetCellRule`](@ref) to manually write to the array where you need to.
`f` is passed an indexable `data` object, and the index of the current cell,
followed by the required grid values for the index.

To update the grid, you can use: [`add!`](@ref), [`sub!`](@ref) for `Number`,
and [`and!`](@ref), [`or!`](@ref) for `Bool`. These methods safely combined
writes from all grid cells - directly using `setindex!` would cause bugs.

## Example

```julia
rule = let x = 10
    Set{Tuple{:a,:b},:b}() do data, I, a, b
        add!(data[:b], a^x, I...)
    end
end
```
The `let` block greatly improves performance.
"""
struct SetCell{R,W,F} <: SetCellRule{R,W}
    "Function to apply to data, index and read state arguments"
    f::F
end
SetCell{R,W}(; kw...) where {R,W} = _nofunctionerror(Set)

@inline function applyrule!(data, rule::SetCell, read, I)
    let data=data, rule=rule, read=read, I=I
        rule.f(data, read, I)
    end
end

"""
    SetNeighbors <: SetNeighborhoodRule

    SetNeighbors(f, neighborhood=Moor(1))
    SetNeighbors{R,W}(f, neighborhood=Moor(1))

A [`SetCellRule`](@ref) to manually write to the array with the specified
neighborhood. Indexing outside the neighborhood is undefined behaviour.

Function `f` is passed an [`SimData`](@ref) object `data`, the specified
neighborhood object and the index of the current cell, followed by the required
grid values for the index.

To update the grid, you can use: [`add!`](@ref), [`sub!`](@ref) for `Number`,
and [`and!`](@ref), [`or!`](@ref) for `Bool`. These methods can be safely combined
writes from all grid cells.

Directly using `setindex!` is possible, but may cause bugs as multiple cells
may write to the same location in an unpredicatble order. As a rule, directly
setting a neighborhood index should only be done if it always sets the samevalue -
then it can be guaranteed that any writes from othe grid cells reach the same result.

[`neighbors`](@ref), [`offsets`](@ref) and [`positions`](@ref) are useful methods for
`SetNeighbors` rules.

## Example

```julia
SetNeighbors{Tuple{:a,:b},:b}() do data, hood, I, a, b
    for pos in positions(hood)
        add!(data[:b], a^2, pos...)
    end
end
```
"""
struct SetNeighbors{R,W,F,N} <: SetNeighborhoodRule{R,W}
    "Function to apply to the data, index and read values"
    f::F
    "The neighborhood of cells around the central cell"
    neighborhood::N
end
SetNeighbors{R,W}(; kw...) where {R,W} = _nofunctionerror(SetNeighbors)
SetNeighbors{R,W}(f; neighborhood=Moore(1)) where {R,W} = SetNeighbors{R,W}(f, neighborhood)

@inline function applyrule!(data, rule::SetNeighbors, read, I)
    let data=data, hood=neighborhood(rule), rule=rule, read=read, I=I
        rule.f(data, hood, read, I)
    end
end

"""
    Convolution <: NeighborhoodRule

    Convolution(f, neighborhood=Moore(1))
    Convolution{R,W}(f, neighborhood=Moore(1))

A `Rule` that runs a basic convolution kernel over the grid.

# Performance

_Always_ use StaticArrays.jl to define the kernel matrix.

Small radius convolutions in DynamicGrids.jl will be faster or comparable to using
DSP.jl or ImageConvolutions.jl. As the radius increases or grid size gets very large
these packages will be a lot faster.

But `Convolution` is convenient to chain into a simulation, and combined with some other
rules. It should perform reasonably well in all but very large simulations or very large
kernels.

## Example

```julia
rule = Convolution(Kernel(SA[0.05 0.1 0.05; 0.1 0.4 0.1; 0.05 0.1 0.05]))
```
"""
struct Convolution{R,W,N} <: NeighborhoodRule{R,W}
    "The neighborhood of cells around the central cell"
    neighborhood::N
end
Convolution{R,W}(A::AbstractArray) where {R,W} = Convolution{R,W}(Kernel(SMatrix{size(A)...}(A)))
Convolution{R,W}(; neighborhood) where {R,W} = Convolution{R,W}(neighborhood)

@inline function applyrule(data, rule::Convolution, read, I)
    @inbounds neighbors(rule) ⋅ kernel(neighborhood(rule))
end

# Utils

@noinline _nofunctionerror(t) =
    throw(ArgumentError("No function passed to $t. Did you mean to use a do block?"))

# Check number of args passed in as we don't get a normal method
# error because of the splatted args in the default constructor.
# @generated function _checkfields(::Type{T}, args::A) where {T,A<:Tuple}
    # length(fieldnames(T)) == length(fieldnames(A)) ? :(nothing) : :(_fielderror(T, args))
# end

# @noinline function _fielderror(T, args)
    # throw(ArgumentError("$T has $(length(fieldnames(T))) fields: $(fieldnames(T)), you have used $(length(args))"))
# end


text:
using Documenter, RRIFT

makedocs(;
    modules=[RRIFT],
    format=Documenter.HTML(),
    pages=[
        "Introduction" => "index.md",
        "Guide" => Any[
            "Downloading" => "guide/downloading.md",
            "Pre-processing" => "guide/preprocessing.md",
            "Fitting" => "guide/fitting.md"
         ]
    ],
    repo="https://github.com/notZaki/RRIFT.jl/blob/{commit}{path}#L{line}",
    sitename="RRIFT.jl",
    authors="Zaki A",
)

rm("./docs/build/guide/data", recursive = true)

deploydocs(;
    repo="github.com/notZaki/RRIFT.jl",
)


text:
"""
    MPIError

Error thrown when an MPI function returns an error code. The `code` field contains the MPI error code.
"""
struct MPIError <: Exception
    code::Cint
end
function Base.show(io::IO, err::MPIError)
    print(io, "MPIError(", err.code, "): ", error_string(err))
end

"""
    FeatureLevelError

Error thrown if a feature is not implemented in the current MPI backend.
"""
struct FeatureLevelError <: Exception
    function_name::Symbol
    min_version::VersionNumber # minimal MPI version required for this feature to be available
end
function Base.show(io::IO, err::FeatureLevelError)
    print(io, "FeatureLevelError($(err.function_name)): Minimum MPI version is $(err.min_version)")
end

macro mpichk(expr, min_version=nothing)
    if !isnothing(min_version) && expr.args[2].head == :tuple
        fn = expr.args[2].args[1].value
        if isnothing(dlsym(libmpi_handle, fn; throw_error=false))
            return quote
                throw(FeatureLevelError($(QuoteNode(fn)), $min_version))
            end
        end
    end

    expr = macroexpand(@__MODULE__, :(@mpicall($expr)))
    # MPI_SUCCESS is defined to be 0
    :((errcode = $(esc(expr))) == 0 || throw(MPIError(errcode)))
end

function error_string(err::MPIError)
    len_ref = Ref{Cint}()
    str_buf = Vector{UInt8}(undef, Consts.MPI_MAX_ERROR_STRING)
    # int MPI_Error_string(int errorcode, char *string, int *resultlen)
    @mpichk ccall((:MPI_Error_string, libmpi), Cint,
                  (Cint, Ptr{UInt8}, Ref{Cint}), err.code, str_buf, len_ref)
    resize!(str_buf, len_ref[])
    return String(str_buf)
end


text:
# Utilities to rewrite UIDs in .squashfs files.

"""
    getuid()

Wrapper around libc's `getuid()` function
"""
function getuid()
    return ccall(:getuid, Cint, ())
end

"""
    getgid()

Wrapper around libc's `getgid()` function
"""
function getgid()
    return ccall(:getgid, Cint, ())
end

# Note: produce these values by #including squashfs_fs.h from linux in Cxx.jl
# and running the indicated command
const offsetof_id_table_start = 0x30    # offsetof(struct suqashfs_super_block, id_table_start)
const offsetof_no_ids = 0x1a            # offsetof(struct suqashfs_super_block, no_ids)

# From squashfs_fs.h
const SQUASHFS_COMPRESSED_BIT = UInt16(1) << 15
const SQUASHFS_MAGIC = 0x73717368

"""
    generate_per_uid_squashfs(cs, new_uid = getuid())

In order for the sandbox to work well, we need to have the uids of the squashfs
images match the uid of the current unprivileged user.  Unfortunately there is
no mount-time option to do this for us.  Fortunately, squashfs is simple enough
that if the ID table is uncompressed, we can just manually patch the uids to be
what we need.  This function performs this operation, by rewriting all UIDs and
GIDs to the given `new_uid` (which defaults to the current user's UID).
"""
function generate_per_uid_squashfs(cs, new_uid = getuid(); verbose::Bool = false)
    # Because the source .squashfs file comes as an immutable artifact, we will create
    # a copy that is based off of the .squashfs artifact, but with rewritten UIDs
    cache_path = shard_path(cs; uid=new_uid)

    # Regenerate the cache if it doesn't exist.
    if !isfile(cache_path)
        filename = artifact_name(cs)
        # Grab the progenitor path, downloading it if necessary
        progenitor_path = joinpath(artifact_path(shard_source_artifact_hash(cs)), filename)
        if !isfile(progenitor_path)
            error("Compiler shard $(filename) missing from disk at $(progenitor_path)")
        end

        # Copy .squashfs file over to our local directory, make it writable
        cp(progenitor_path, cache_path)
        chmod(cache_path, 0o644)

        open(cache_path, "r+") do file
            # Check magic
            if read(file, UInt32) != SQUASHFS_MAGIC
                error("`$progenitor_path` is not a squashfs file")
            end
            # Check that the image contains only one id (which we will rewrite)
            seek(file, offsetof_no_ids)
            if read(file, UInt16) != 1
                error("`$progenitor_path` uses more than one uid/gid")
            end
            # Find the index table
            seek(file, offsetof_id_table_start)
            offset = read(file, UInt64)
            seek(file, offset)
            # Find the correct metdata block
            index = read(file, UInt64)
            seek(file, index)
            # Read the metadata block
            size = read(file, UInt16)
            # Make sure it's uncompressed (yes, I know that flag is terribly
            # named - it indicates that the data is uncompressed)
            if ((size & SQUASHFS_COMPRESSED_BIT) == 0)
                error("Metadata block is compressed")
            end
            p = position(file)
            uid = read(file, UInt32)
            if uid == new_uid
                return
            end
            if verbose
                @info("Rewriting $(basename(progenitor_path)) from UID $(uid) -> $(new_uid)")
            end
            seek(file, p)
            write(file, UInt32(new_uid))
        end
    end

    # Finally, return the path to that cached squashfs file
    return cache_path
end


text:
include("header.jl")

using Statistics
meanabs(x)=mean(abs,x)
meanabs2(x)=mean(abs2,x)

@testset "statistics" begin
    o = (:delta=>0.0001,:atol=>0.01,:rtol=>0.01)
    for T in (Float32,Float64)
        @test gradcheck(mean, randn(T,2,3); o...)
        @test gradcheck(mean, randn(T,2,3); kw=[:dims=>1], o...)
        @test gradcheck(mean, randn(T,2,3); kw=[:dims=>(1,2)], o...)
        @test gradcheck(meanabs, randn(T,2,3); o...)
        @test gradcheck(meanabs2, randn(T,2,3); o...)
        @test gradcheck(var, randn(T,2,3); o...)
        @test gradcheck(var, randn(T,2,3); kw=[:dims=>1], o...)
        @test gradcheck(var, randn(T,2,3); kw=[:dims=>(1,2)], o...)
        @test gradcheck(std, randn(T,2,3); o...)
        @test gradcheck(std, randn(T,2,3); kw=[:dims=>1], o...)
        @test gradcheck(std, randn(T,2,3); kw=[:dims=>(1,2)], o...)
    end
end


text:
const RL = CommonRLInterface

abstract type AbstractPOMDPsCommonRLEnv <: RL.AbstractEnv end

RL.actions(env::AbstractPOMDPsCommonRLEnv) = actions(env.m)
RL.terminated(env::AbstractPOMDPsCommonRLEnv) = isterminal(env.m, env.s)

mutable struct MDPCommonRLEnv{RLO, M<:MDP, S} <: AbstractPOMDPsCommonRLEnv
    m::M
    s::S
end

"""
    MDPCommonRLEnv(m, [s])
    MDPCommonRLEnv{RLO}(m, [s])

Create a CommonRLInterface environment from MDP m; optionally specify the state 's'.

The `RLO` parameter can be used to specify a type to convert the observation to. By default, this is `AbstractArray`. Use `Any` to disable conversion.
"""
MDPCommonRLEnv{RLO}(m, s=rand(initialstate(m))) where {RLO} = MDPCommonRLEnv{RLO, typeof(m), statetype(m)}(m, s)
MDPCommonRLEnv(m, s=rand(initialstate(m))) = MDPCommonRLEnv{AbstractArray}(m, s)

function RL.reset!(env::MDPCommonRLEnv)
    env.s = rand(initialstate(env.m))
    return nothing
end

function RL.act!(env::MDPCommonRLEnv, a)
    sp, r = @gen(:sp, :r)(env.m, env.s, a)
    env.s = sp
    return r
end

RL.observe(env::MDPCommonRLEnv{RLO}) where {RLO} = convert_s(RLO, env.s, env.m)

RL.@provide RL.clone(env::MDPCommonRLEnv{RLO}) where {RLO} = MDPCommonRLEnv{RLO}(env.m, env.s)
RL.@provide RL.render(env::MDPCommonRLEnv) = render(env.m, (sp=env.s,))
RL.@provide RL.state(env::MDPCommonRLEnv{RLO}) where {RLO} = convert_s(RLO, env.s, env.m)
RL.@provide RL.valid_actions(env::MDPCommonRLEnv) = actions(env.m, env.s)

RL.observations(env::MDPCommonRLEnv{RLO}) where {RLO} = (convert_s(RLO, s, env.m) for s in states(env.m)) # should really be some kind of lazy map that handles uncountably infinite spaces
RL.provided(::typeof(RL.observations), ::Type{<:Tuple{MDPCommonRLEnv{<:Any, M, <:Any}}}) where {M} = static_hasmethod(states, Tuple{<:M})

RL.@provide function RL.setstate!(env::MDPCommonRLEnv{<:Any, <:Any, S}, s) where S
    env.s = convert_s(S, s, env.m)
    return nothing
end

mutable struct POMDPCommonRLEnv{RLO, M<:POMDP, S, O} <: AbstractPOMDPsCommonRLEnv
    m::M
    s::S
    o::O
end

"""
    POMDPCommonRLEnv(m, [s], [o])
    POMDPCommonRLEnv{RLO}(m, [s], [o])

Create a CommonRLInterface environment from POMDP m; optionally specify the state 's' and observation 'o'.

The `RLO` parameter can be used to specify a type to convert the observation to. By default, this is `AbstractArray`. Use `Any` to disable conversion.
"""
POMDPCommonRLEnv{RLO}(m, s=rand(initialstate(m)), o=rand(initialobs(m, s))) where {RLO} = POMDPCommonRLEnv{RLO, typeof(m), statetype(m), obstype(m)}(m, s, o)
POMDPCommonRLEnv(m, s=rand(initialstate(m)), o=rand(initialobs(m, s))) = POMDPCommonRLEnv{AbstractArray}(m, s, o)

function RL.reset!(env::POMDPCommonRLEnv)
    env.s = rand(initialstate(env.m))
    env.o = rand(initialobs(env.m, env.s))
    return nothing
end

function RL.act!(env::POMDPCommonRLEnv, a)
    sp, o, r = @gen(:sp, :o, :r)(env.m, env.s, a)
    env.s = sp
    env.o = o
    return r
end

RL.observe(env::POMDPCommonRLEnv{RLO}) where {RLO} = convert_o(RLO, env.o, env.m)

RL.@provide RL.clone(env::POMDPCommonRLEnv{RLO}) where {RLO} = POMDPCommonRLEnv{RLO}(env.m, env.s, env.o)
RL.@provide RL.render(env::POMDPCommonRLEnv) = render(env.m, (sp=env.s, o=env.o))
RL.@provide RL.state(env::POMDPCommonRLEnv) = (env.s, env.o)
RL.@provide RL.valid_actions(env::POMDPCommonRLEnv) = actions(env.m, env.s)

RL.observations(env::POMDPCommonRLEnv{RLO}) where {RLO} = (convert_o(RLO, o, env.m) for o in observations(env.m)) # should really be some kind of lazy map that handles uncountably infinite spaces
RL.provided(::typeof(RL.observations), ::Type{<:Tuple{POMDPCommonRLEnv{<:Any, M, <:Any, <:Any}}}) where {M} = static_hasmethod(observations, Tuple{<:M})

RL.@provide function RL.setstate!(env::POMDPCommonRLEnv, so)
    env.s = first(so)
    env.o = last(so)
    return nothing
end

Base.convert(::Type{RL.AbstractEnv}, m::POMDP) = POMDPCommonRLEnv(m)
Base.convert(::Type{RL.AbstractEnv}, m::MDP) = MDPCommonRLEnv(m)

Base.convert(::Type{MDP}, env::MDPCommonRLEnv) = env.m
Base.convert(::Type{POMDP}, env::POMDPCommonRLEnv) = env.m


text:
"""
    writeMetric( mesh::Mesh2D, mesher::MesherBAMG,
                      metric::Vector{Float64}, flname::String )

Writes 2D isotropic metric to BAMG format in `flname`.
"""
function writeMetric( mesh::Mesh2D, mesher::MesherBAMG,
                      metric::Vector{Float64}, flname::String )

    np = size(mesh.p,1)
    if np != length(metric)
        error(" Isotropic metric error: size of metric does not correspond with size of mesh")
    end

    fid = open( flname, "w" )

    @printf( fid, "%i 1\n", np )

    # Vertices
    for jj = 1:np
        @printf(fid, "%16.15e\n", metric[jj])
    end

    close(fid)

    # TODO: Corners

end

"""
    writeMetric( mesh::Mesh2D, mesher::MesherBAMG,
                      metric::Matrix{Float64}, flname::String )

Writes 2D anisotropic metric to BAMG format in `flname`.
"""
function writeMetric( mesh::Mesh2D, mesher::MesherBAMG,
                      metric::Matrix{Float64}, flname::String )

    np = size(mesh.p,1)
    if np != size(metric,1)
        error(" Ansotropic metric error: size of metric does not correspond with size of mesh")
    end
    if size(metric,2) != 3
        error(" Ansotropic metric error: wrong number of entries for metric")
    end

    fid = open( flname, "w" )

    @printf( fid, "%i 3\n", np )

    # Vertices
    for jj = 1:np
        @printf(fid, "%16.15e\t%16.15e\t%16.15e\n", metric[jj,1], metric[jj,2], metric[jj,3])
    end

    close(fid)

    # TODO: Corners

end


text:
using WaveletUtils, StaticArrays, Test

@testset "Detail" begin
    include("detail.jl")
end

@testset "Collections" begin
    include("collections.jl")
end

@testset "Trees" begin
    include("trees.jl")
end


text:
d = readline("input.txt")
garbage = false
score, depth, numGarbage = 0, 1, 0
i = 1
while i <= length(d)
    if d[i] == '!'
        i += 1
    elseif garbage && d[i] != '>'
        numGarbage += 1
    elseif d[i] == '<'
        garbage = true
    elseif d[i] == '>'
        garbage = false
    elseif d[i] == '{'
        score += depth
        depth += 1
    elseif d[i] == '}'
        depth -= 1
    end
    i += 1
end

println(score)
println(numGarbage)


text:
module NaturallyUnitful

using Reexport

@reexport using Unitful
export natural, unnatural

using Unitful: c0, ħ, k, ϵ0

const ħc0 = ħ * c0

setDimPow(D::Dict, x::Unitful.Dimension{T}) where {T} = D[T] = x.power

function dimDict(x::Unitful.Dimensions{T}) where {T}
    D = Dict{Symbol,Rational}(
        :Length => 0, :Mass => 0, :Temperature => 0, :Time => 0, :Current => 0
    )
    for t in T
        setDimPow(D, t)
    end
    return D
end

dimDict(x::Unitful.Quantity{T}) where {T} = dimDict(dimension(x))

function dimDict(x)
    return Dict{Symbol,Rational}(
        :Length => 0, :Mass => 0, :Temperature => 0, :Time => 0, :Current => 0
    )
end

gettypeparams(::Unitful.FreeUnits{T,U,V}) where {T,U,V} = T, U, V
const energydimension = gettypeparams(u"GeV")[2]

"""
    natural(q; base=u"eV")

Convert `q` to natural units based on the units specified by `base`. If `base` is
unspecified, `natural` will default to `eV`. Currently, `natural` only supports
`base`s with dimensions of energy. For all other `base` you will need to use `unnatrual`
to convert.

Examples:

    julia> natural(1u"kg")
    5.609588650020686e35 eV

    julia> natural(1u"kg", base=u"GeV")
    5.609588650020685e26 GeV

    julia> natural(1u"m")
    5.067730759202785e6 eV^-1

    julia> natural(1u"m", base=u"GeV")
    5.067730759202785e15 GeV^-1
"""
natural(q; base=u"eV") = _natural(base, q)

function _natural(base::Unitful.FreeUnits{T,energydimension,U}, q) where {T,U}
    D = dimDict(q)
    (α, β, γ, δ, ϕ) = (D[:Length], D[:Mass], D[:Temperature], D[:Time], D[:Current])
    return uconvert(
        base^(-α - δ + β + γ + ϕ),
        q * ħc0^(-α + ϕ / 2) * ħ^(-δ) * c0^(2(β - ϕ / 2)) * k^(γ) * (4π * ϵ0)^(-ϕ / 2),
    )
end

function _natural(base::Unitful.FreeUnits{T,U,V}, q) where {T,U,V}
    return throw(
        """natural(q; base)` where `base` has dimensions `$U` has not yet been implemented. Please use a base with dimensions of `energy` and then use `unnatural` to convert to your desired units.""",
    )
end

"""
    unnatural(targetUnit, q)

Convert a quantity `q` to units specified by `targetUnits` automatically inserting whatever
factors  of `ħ`, `c0`, `ϵ₀` and `k` (Boltzmann constant) are required to make the conversion work.

Examples:

    julia> unnatural(u"m", 5.067730759202785e6u"eV^-1")
    1.0 m

    julia> unnatural(u"m/s", 1)
    2.99792458e8 m s^-1

    julia> unnatural(u"K", 1u"eV")
    11604.522060401006 K
"""
function unnatural(targetUnit, q)
    natTarget = natural(1targetUnit)
    natQ = natural(q)
    ratio = natQ / natTarget
    if typeof(dimension(ratio)) == Unitful.Dimensions{()}
        (ratio)targetUnit
    else
        throw(Unitful.DimensionError)
    end
end

end # module


text:
@testset "Special cases" begin

    T = Float64
    f(l) = n -> gamma(l+1/2)/gamma(2l)*gamma(n+2l)/gamma(n+l+1/2)

    for (P, Q, fn) in ((Gegenbauer{1/2, T}, Legendre{T}, n->1.0),

                       (Gegenbauer{1/4, T}, Jacobi{1/4-1/2, 1/4-1/2, T}, f(1/4)),
                       (Gegenbauer{3/4, T}, Jacobi{3/4-1/2, 3/4-1/2, T}, f(3/4)),
                       (Gegenbauer{1, T}, Jacobi{1-1/2, 1-1/2, T}, f(1))
                       )
        for i in 1:5
            @test convert(Polynomial, Polynomials.basis(P,i)) ≈ fn(i) * convert(Polynomial, Polynomials.basis(Q,i))
        end

    end

end


text:
module NavAbilitySDK

const NVA = NavAbilitySDK
export NVA

# Global imports
using DocStringExtensions
using LinearAlgebra
using JSON
using UUIDs
using Downloads
using HTTP

# for overloading with visualization helpers
import Base: show

# LinearAlgebra pass through exports
export diagm, norm
# UUIDs pass through exports
export uuid4

# Graphql
include("./navability/graphql/Factor.jl")
include("./navability/graphql/Status.jl")
include("./navability/graphql/Variable.jl")
include("./navability/graphql/DataBlobs.jl")

include("./navability/graphql/QueriesDeprecated.jl")
export QUERY_VARIABLE_LABELS
export QUERY_FILES, QUERY_CALIBRATION
export MUTATION_ADDVARIABLE, MUTATION_ADDFACTOR, MUTATION_ADDSESSIONDATA
export MUTATION_SOLVESESSION, MUTATION_SOLVEFEDERATED
export MUTATION_DEMOCANONICALHEXAGONAL
export MUTATION_CREATE_UPLOAD, MUTATION_ABORT_UPLOAD, MUTATION_COMPLETE_UPLOAD
export MUTATION_PROC_CALIBRATION
export SUBSCRIPTION_UPDATES

# Entities
include("./navability/entities/NavAbilityClient.jl")
include("./navability/entities/Client.jl")
include("./navability/entities/Common.jl")
include("./navability/entities/Distributions.jl")
include("./navability/entities/Variable.jl")
include("./navability/entities/InferenceTypes.jl")
include("./navability/entities/Factor.jl")
export NavAbilityClient, NavAbilityWebsocketClient, NavAbilityHttpsClient, QueryOptions, MutationOptions
export Client, Scope
export QueryDetail, LABEL, SKELETON, SUMMARY, FULL
export Distribution, Normal, Rayleigh, FullNormal, Uniform, Categorical
export ManifoldKernelDensity
export VariableType, Variable
export FactorData, PriorData, PriorPose2Data, PriorPoint2Data, LinearRelativeData, Pose2Pose2Data, Pose2AprilTag4CornersData, Pose2Point2BearingRangeData, Point2Point2RangeData, MixtureData
export ScatterAlignPose2Data
export FactorType, Factor

# Services
include("./navability/services/Variable.jl")
include("./navability/services/Factor.jl")
include("./navability/services/Solve.jl")
include("./navability/services/Status.jl")
include("./navability/services/Utils.jl")
include("./navability/services/StandardAPI.jl")
include("./navability/services/DataBlobs.jl")
export getVariable, getVariables, listVariables, ls
export addVariable, addPackedVariable
export getFactor, getFactors, listFactors, lsf
export addFactor, addPackedFactor
export solveSession, solveFederated
export getStatusMessages, getStatusLatest, getStatusesLatest
export waitForCompletion
export GraphVizApp, MapVizApp

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# script to generate tables of common symbols from a system image
# steps to rerun this:
#   1. empty the contents of common_symbols*.inc
#   2. modify dump.c to write and read a '\0' after every symbol name
#   3. build sys.ji
#   4. cd src && ../julia gen_sysimg_symtab.jl ../usr/lib/julia/sys.ji

fname = ARGS[1]

io,_ = open(pipeline(`strings -n 3 $fname`,
                     `tr -d "() \t+-"`,
                     `grep -v '^$'`,  # rm empty lines
                     `grep -v '\\'`,  # avoid backslashes
                     `sort`, `uniq -c`, `sort -g -r`,
                     `head -n 315`))  # 63 + 252

function outputline(io, line)
    row = split(line, " ", keep=false)
    println(io, "jl_symbol(\"", row[2], "\"),")
end

lines = eachline(io)

open(f->foreach(l->outputline(f,l), take(lines,63)), "common_symbols1.inc", "w")
open(f->foreach(l->outputline(f,l), lines), "common_symbols2.inc", "w")


text:
"""
    TerminalState

A type with no fields whose singleton instance `terminalstate` is used to represent a terminal state with no additional information.

This type has the appropriate promotion logic implemented to function like `Missing` when added to arrays, etc.

Note that terminal states NEED NOT be of type `TerminalState`. You can define any state to be terminal by implementing the appropriate `isterminal` method. Solvers and simulators SHOULD NOT check for this type, but should instead check using `isterminal`.
"""
struct TerminalState end

"""
    terminalstate

The singleton instance of type `TerminalState` representing a terminal state.
"""
const terminalstate = TerminalState()

isterminal(m::Union{MDP,POMDP}, ts::TerminalState) = true
Base.promote_rule(::Type{TerminalState}, T::Type) = Union{TerminalState, T}


text:
# Convert a ReactionSystem type into a BondGraph type
function BondGraph(rs::ReactionSystem; chemostats=[])
    bg = BondGraph(rs.name)

    # Create disjoint reaction bondgraphs for each reaction in network
    all_reactions = reactions(rs)
    reaction_num = 1
    for (i, reaction) in enumerate(all_reactions)
        if i > 1 && _is_reverse_off_previous(reaction, all_reactions[i-1])
            # Skip the second reaction
            continue
        end

        Re = Component(:Re, Symbol("R$reaction_num"), numports=2)
        add_node!(bg, Re)
        _half_equation!(bg, reaction.substrates, reaction.substoich, Re, chemostats)
        _half_equation!(bg, reaction.products, reaction.prodstoich, Re, chemostats)

        reaction_num += 1
    end

    # Combine common species across reactions
    species_names = _stringify_species.(species(rs))
    for spcs_name in species_names
        spcs_nodes = getnodes(bg, spcs_name)
        for node in spcs_nodes[2:end]
            merge_nodes!(bg, spcs_nodes[1], node)
        end
    end

    simplify_junctions!(bg)
end

# If this reaction is an exact reverse of the previous reaction,
# then together they form a bi-directional reaction pair.
_is_reverse_off_previous(r1, r2) = Set(r1.substrates) == Set(r2.products) && Set(r2.substrates) == Set(r1.products)

function _half_equation!(bg, species, stoich, Re, chemostats)
    species_names = _stringify_species.(species)

    if length(species) > 1
        one_junction = EqualFlow()
        add_node!(bg, one_junction)

        for (i, spcs) in enumerate(species_names)
            comp = spcs in chemostats ? Component(:Se, Symbol(spcs)) : Component(:Ce, Symbol(spcs))
            add_node!(bg, comp)
            connect!(bg, comp, one_junction)

            if stoich[i] != 1
                tf = Component(:TF, Symbol("$(stoich[i])"), numports=2)
                insert_node!(bg, (comp, one_junction), tf)
            end
        end

        connect!(bg, one_junction, Re)
    else
        spcs = species_names[1]
        comp = spcs in chemostats ? Component(:Se,  Symbol(spcs)) : Component(:Ce,  Symbol(spcs))
        add_node!(bg, comp)
        connect!(bg, comp, Re)

        if stoich[1] != 1
            tf = Component(:TF, Symbol("$(stoich[1])"), numports=2)
            insert_node!(bg, (comp, Re), tf)
        end
    end
end

# removes "(t)" from the end of the species name
_stringify_species(species) = string(species)[1:end-3]

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
abstract type AbstractMetadata{D<:AbstractDict{Symbol,Any}} <: AbstractDict{Symbol,Any} end

Base.haskey(m::AbstractMetadata, k) = haskey(properties(m), k)

Base.empty!(m::AbstractMetadata) = empty!(properties(m))

Base.get(m::AbstractMetadata, k, default) = get(properties(m), k, default)

Base.get!(m::AbstractMetadata, k, default) = get!(properties(m), k, default)

Base.in(k, m::AbstractMetadata) = in(k, propname(m))

Base.getkey(m::AbstractMetadata, k, default) = getkey(properties(m), k, default)

Base.pop!(m::AbstractMetadata, k) = pop!(properties(m), k)

Base.pop!(m::AbstractMetadata, k, default) = pop!(properties(m), k, default)

Base.isempty(m::AbstractMetadata) = isempty(properties(m))

Base.length(m::AbstractMetadata) = length(properties(m))

Base.iterate(m::AbstractMetadata) = iterate(properties(m))

Base.iterate(m::AbstractMetadata, state) = iterate(properties(m), state)

Base.delete!(m::AbstractMetadata, k) = delete!(properties(m), k)

@inline Base.getindex(x::AbstractMetadata, s::Symbol) = getindex(properties(x), s)

@inline function Base.setindex!(x::AbstractMetadata, val, s::Symbol)
    return setindex!(properties(x), val, s)
end

function ImageMetadata.properties(m::AbstractMetadata)
    error("All subtypes of AbstractMetadata must implement a properties method.")
end

"""
    NeuroMetadata

Dictionary structure for storing neuroscience related metadata.
"""
struct NeuroMetadata{D} <: AbstractMetadata{D}
    _properties::D
end

function NeuroMetadata(; kwargs...)
    out = Dict{Symbol,Any}()
    for (k,v) in kwargs
        out[k] = v
    end
    return NeuroMetadata(out)
end

ImageMetadata.properties(nm::NeuroMetadata) = getfield(nm, :_properties)

Base.copy(nm::NeuroMetadata) = NeuroMetadata(copy(properties(nm)))

Base.filter!(pred, m::NeuroMetadata) = filter!(pred, properties(m))

@inline Base.getproperty(x::NeuroMetadata, s::Symbol) = neuroproperty(x, s)

@inline Base.setproperty!(x::NeuroMetadata, s::Symbol, val) = neuroproperty!(x, s, val)

function Base.propertynames(x::NeuroMetadata)
    return isempty(x) ? Tuple(PROPERTIES) : Tuple(merge(PROPERTIES, [keys(x)...]))
end


text:
#
# Factor-DataFrame conversion
#

function Base.convert(::Type{DataFrames.DataFrame}, ϕ::Factor)
    df = DataFrames.DataFrame(pattern_states(ϕ))
    DataFrames.rename!(df, names(df), names(ϕ))
    df[:potential] = ϕ.potential[:]

    DataFrames.pool!(df, names(ϕ))

    return df
end


text:
# residual
function r!(r, z, θ, κ)
    @warn "residual not defined"
    nothing
end

# residual Jacobian wrt z
function rz!(rz, z, θ)
    @warn "residual Jacobian wrt z not defined"
    nothing
end

# residual Jacobian wrt θ
function rθ!(rθ, z, θ)
    @warn "residual Jacobian wrt θ not defined"
    nothing
end

mutable struct ResidualMethods
    r!
    rz!
    rθ!
end

# optimization spaces
abstract type Space end

# Euclidean
struct Euclidean <: Space
    n::Int
end

function candidate_point!(z̄::Vector{T}, ::Euclidean, z::Vector{T}, Δ::Vector{T}, α::T) where T
    z̄ .= z
    z̄ .-= α .* Δ
    return z̄
end

function update_point!(z::Vector{T}, ::Space, z̄::Vector{T}) where T
    z .= z̄
end

function mapping!(δz, s::Euclidean, δzs, z) # TODO: make allocation free
    δz .= δzs
end

# interior-point solver options
@with_kw mutable struct InteriorPointOptions{T}
    r_tol::T = 1.0e-5
    κ_tol::T = 1.0e-5
    ls_scale::T = 0.5
    max_iter::Int = 100
    max_ls::Int = 3
    max_time::T = 1e5
    diff_sol::Bool = false
    reg::Bool = false
    ϵ_min = 0.05 # ∈ [0.005, 0.25]
        # smaller -> faster
        # larger  -> slower, more robust
    κ_reg = 1e-3 # bilinear constraint violation level at which regularization is triggered [1e-3, 1e-4]
    γ_reg = 1e-1 # regularization scaling parameters ∈ [0, 0.1]:
        # 0   -> faster & ill-conditioned
        # 0.1 -> slower & better-conditioned
        # simulation choose γ_reg = 0.1
        # MPC choose γ_reg = 0.0
    solver::Symbol = :lu_solver
    undercut::T = 2.0 # the solver will aim at reaching κ_vio = κ_tol / undercut
        # simulation: Inf
        # MPC: 5.0
    verbose::Bool = false
    warn::Bool = false
end

mutable struct InteriorPoint{T,R,RZ,Rθ}
    s::Space
    idx::IndicesOptimization
    methods::ResidualMethods
    z::Vector{T}               # current point
    z̄::Vector{T}               # candidate point
    r::R                       # residual
    rz::RZ                     # residual Jacobian wrt z
    rθ::Rθ                     # residual Jacobian wrt θ
    Δ::Vector{T}               # search direction
    δz::Matrix{T}              # solution gradients (this is always dense)
    δzs::Matrix{T}             # solution gradients (in optimization space; δz = δzs for Euclidean)
    θ::Vector{T}               # problem data
    solver::LinearSolver
    z_reg::Vector{T}           # current point w/ regularization on cone variables
    reg_val::Vector{T}
    iterations::Int
    opts::InteriorPointOptions{T}
    κ::Vector{T}
    zort::Vector{Vector{T}}
    Δort::Vector{Vector{T}}
    zsoc::Vector{Vector{Vector{T}}}
    Δsoc::Vector{Vector{Vector{T}}}
    ρv::Vector{Vector{Vector{T}}}
    idx_soce::Vector{Vector{Vector{Int}}}
end

function interior_point(z, θ;
        s = Euclidean(length(z)),
        idx = IndicesOptimization(),
        r! = r!, rz! = rz!, rθ! = rθ!,
        r  = zeros(idx.nΔ),
        rz = spzeros(idx.nΔ, idx.nΔ),
        rθ = spzeros(idx.nΔ, length(θ)),
        opts::InteriorPointOptions = InteriorPointOptions()) where T

    rz!(rz, z, θ) # compute Jacobian for pre-factorization
    num_data = length(θ)

    zort = [zeros(length(idx.ortz[1])), zeros(length(idx.ortz[2]))]
    Δort = [zeros(length(idx.ortΔ[1])), zeros(length(idx.ortΔ[2]))]
    zsoc = [[zeros(length(i)), zeros(length(i))] for i in idx.socri]
    Δsoc = [[zeros(length(i)), zeros(length(i))] for i in idx.socri]
    ρv = [[zeros(max(0, length(i) - 1)), zeros(max(0, length(i) - 1))] for i in idx.socri]
    idx_soce = [[collect(2:length(i)), collect(2:length(i))] for i in idx.socri]

    InteriorPoint{typeof.([z[1], r, rz, rθ])...}(
        s,
        idx,
        ResidualMethods(r!, rz!, rθ!),
        z,
        zero(z),
        r,
        rz,
        rθ,
        zeros(idx.nΔ),
        zeros(idx.nz, num_data),
        zeros(idx.nΔ, num_data),
        θ,
        eval(opts.solver)(rz),
        zero(z),
        zeros(1),
        0,
        opts,
        zeros(1),
        zort,
        Δort,
        zsoc,
        Δsoc,
        ρv,
        idx_soce)
end

# interior point solver
function interior_point_solve!(ip::InteriorPoint{T,R,RZ,Rθ}) where {T,R,RZ,Rθ} #TODO: change to solve!
    ip.opts.max_time < 0.1 ? _interior_point_solve_time_budget!(ip) : _interior_point_solve!(ip)
end

function _interior_point_solve!(ip::InteriorPoint{T,R,RZ,Rθ}) where {T,R,RZ,Rθ} #TODO: change to solve!

    # space
    s = ip.s

    # options
    opts = ip.opts
    r_tol = opts.r_tol
    κ_tol = opts.κ_tol
    ls_scale = opts.ls_scale
    max_iter = opts.max_iter
    max_time = opts.max_time
    max_ls = opts.max_ls
    diff_sol = opts.diff_sol
    ϵ_min = opts.ϵ_min
    κ_reg = opts.κ_reg
    γ_reg = opts.γ_reg
    reg = opts.reg
    verbose = opts.verbose
    warn = opts.warn

    # unpack pre-allocated data
    z = ip.z
    z̄ = ip.z̄
    r = ip.r
    rz = ip.rz
    Δ = ip.Δ
    θ = ip.θ

    # regularization
    ip.reg_val[1] = 0.0 # reset

    zort = ip.zort
    Δort = ip.Δort
    zsoc = ip.zsoc
    Δsoc = ip.Δsoc
    ρv = ip.ρv

    # indices
    idx = ip.idx
    ortz = idx.ortz
    ortΔ = idx.ortΔ
    socz = idx.socz
    socΔ = idx.socΔ
    bil = idx.bil
    ortr = idx.ortr
    socr = idx.socr
    sorci = idx.socri
    soce = ip.idx_soce

    # initialization
    solver = ip.solver
    ip.iterations = 0

    # evaluate residual
    ip.methods.r!(r, z, θ, 0.0)

    # evaluate bilinear constraint violation
    κ_vio = bilinear_violation(ip, r)

    # evaluate equality constraint violation
    r_vio = residual_violation(ip, r)

    # begin timer
    elapsed_time = 0.0

    for j = 1:max_iter
        elapsed_time >= max_time && break
        elapsed_time += @elapsed begin
            # check for converged residual
            if (r_vio < r_tol) && (κ_vio < κ_tol)
                break
            end
            ip.iterations += 1

            # Compute regularization level
            regularization!(ip, κ_vio, κ_reg, γ_reg)

            # compute residual Jacobian
            rz!(ip, rz, z, θ, reg=ip.reg_val[1]) # this is not adapted to the second-order cone

            # compute step
            linear_solve!(solver, Δ, rz, r, reg=ip.reg_val[1])

            α_ort = ort_step_length(z, Δ, ortz, ortΔ, τ=1.0)
            α_soc = soc_step_length(z, Δ, socz, socΔ, zsoc, Δsoc, ρv, soce, τ=1.0, verbose=false)

            α = min(α_ort, α_soc)
            μ, σ = centering(z, Δ, ortz, ortΔ, socz, socΔ, zort, Δort, zsoc, Δsoc, α)

            # Compute corrector residual
            ip.κ[1] = max(σ * μ , κ_tol / opts.undercut)
            ip.methods.r!(r, z, θ, ip.κ) # here we set κ = σ*μ, Δ = Δaff
            general_correction_term!(r, Δ, ortr, socr, sorci, ortΔ, socΔ)

            # Compute corrector search direction
            linear_solve!(solver, Δ, rz, r, reg=ip.reg_val[1], fact=false)

            τ = max(0.95, 1 - max(r_vio, κ_vio)^2)

            α_ort = ort_step_length(z, Δ, ortz, ortΔ, τ=τ)
            α_soc = soc_step_length(z, Δ, socz, socΔ, zsoc, Δsoc, ρv, soce, τ=min(τ, 0.99), verbose=false)
            α = min(α_ort, α_soc)

            # reduce norm of residual
            candidate_point!(z, s, z, Δ, α)

            κ_vio_cand = 0.0
            r_vio_cand = 0.0
            for i = 1:max_ls
                ip.methods.r!(r, z, θ, 0.0)
                κ_vio_cand = bilinear_violation(ip, r)
                r_vio_cand = residual_violation(ip, r)
                if (r_vio_cand <= r_vio) || (κ_vio_cand <= κ_vio)
                    break
                end
                verbose && println("linesearch $i")

                # backtracking
                candidate_point!(z, s, z, Δ, -α * ls_scale^i)
            end

            # update
            κ_vio = κ_vio_cand
            r_vio = r_vio_cand

            verbose && iterate_info(j, r, r_vio, κ_vio, Δ, α)
        end
    end

    if (r_vio < r_tol) && (κ_vio < κ_tol)
        # differentiate solution
        regularization!(ip, κ_vio, ip.opts.γ_reg)
        diff_sol && differentiate_solution!(ip, reg=ip.reg_val[1])
        return true
    else
        return false
    end
end

function _interior_point_solve_time_budget!(ip::InteriorPoint{T,R,RZ,Rθ}) where {T,R,RZ,Rθ} #TODO: change to solve!
    # space
    s = ip.s

    # options
    opts = ip.opts
    r_tol = opts.r_tol
    κ_tol = opts.κ_tol
    ls_scale = opts.ls_scale
    max_iter = opts.max_iter
    max_time = opts.max_time
    max_ls = opts.max_ls
    diff_sol = opts.diff_sol
    ϵ_min = opts.ϵ_min
    κ_reg = opts.κ_reg
    γ_reg = opts.γ_reg
    reg = opts.reg
    verbose = opts.verbose
    warn = opts.warn

    # unpack pre-allocated data
    z = ip.z
    z̄ = ip.z̄
    r = ip.r
    rz = ip.rz
    Δ = ip.Δ
    θ = ip.θ

    # regularization
    ip.reg_val[1] = 0.0 # reset

    zort = ip.zort
    Δort = ip.Δort
    zsoc = ip.zsoc
    Δsoc = ip.Δsoc
    ρv = ip.ρv

    # indices
    idx = ip.idx
    ortz = idx.ortz
    ortΔ = idx.ortΔ
    socz = idx.socz
    socΔ = idx.socΔ
    bil = idx.bil
    ortr = idx.ortr
    socr = idx.socr
    sorci = idx.socri
    soce = ip.idx_soce

    # initialization
    solver = ip.solver
    ip.iterations = 0

    # evaluate residual
    ip.methods.r!(r, z, θ, 0.0)

    # evaluate bilinear constraint violation
    κ_vio = bilinear_violation(ip, r)

    # evaluate equality constraint violation
    r_vio = residual_violation(ip, r)

    # begin timer
    elapsed_time = 0.0

    for j = 1:max_iter
        elapsed_time += @elapsed begin
            # check for converged residual
            if (r_vio < r_tol) && (κ_vio < κ_tol)
                break
            end
            ip.iterations += 1
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            # Compute regularization level
            regularization!(ip, κ_vio, κ_reg, γ_reg)
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            # compute residual Jacobian
            rz!(ip, rz, z, θ, reg=ip.reg_val[1]) # this is not adapted to the second-order cone
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            # compute step
            linear_solve!(solver, Δ, rz, r, reg=ip.reg_val[1])
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            α_ort = ort_step_length(z, Δ, ortz, ortΔ, τ=1.0)
            α_soc = soc_step_length(z, Δ, socz, socΔ, zsoc, Δsoc, ρv, soce, τ=1.0, verbose=false)
            α = min(α_ort, α_soc)
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            μ, σ = centering(z, Δ, ortz, ortΔ, socz, socΔ, zort, Δort, zsoc, Δsoc, α)
            # Compute corrector residual
            ip.κ[1] = max(σ * μ , κ_tol / opts.undercut)
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            ip.methods.r!(r, z, θ, ip.κ) # here we set κ = σ*μ, Δ = Δaff
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            general_correction_term!(r, Δ, ortr, socr, sorci, ortΔ, socΔ)
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            # Compute corrector search direction
            linear_solve!(solver, Δ, rz, r, reg=ip.reg_val[1], fact=false)
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            τ = max(0.95, 1 - max(r_vio, κ_vio)^2)

            α_ort = ort_step_length(z, Δ, ortz, ortΔ, τ=τ)
            α_soc = soc_step_length(z, Δ, socz, socΔ, zsoc, Δsoc, ρv, soce, τ=min(τ, 0.99), verbose=false)
            α = min(α_ort, α_soc)
        end
        elapsed_time >= max_time && (return false)

        elapsed_time += @elapsed begin
            # reduce norm of residual
            candidate_point!(z, s, z, Δ, α)
        end
        elapsed_time >= max_time && (return false)

            κ_vio_cand = 0.0
            r_vio_cand = 0.0
            for i = 1:max_ls
                elapsed_time += @elapsed begin
                    ip.methods.r!(r, z, θ, 0.0)
                end
                elapsed_time >= max_time && (return false)

                elapsed_time += @elapsed begin
                    κ_vio_cand = bilinear_violation(ip, r)
                    r_vio_cand = residual_violation(ip, r)
                    if (r_vio_cand <= r_vio) || (κ_vio_cand <= κ_vio)
                        break
                    end
                    verbose && println("linesearch $i")
                end
                elapsed_time >= max_time && (return false)

                elapsed_time += @elapsed begin
                    # backtracking
                    candidate_point!(z, s, z, Δ, -α * ls_scale^i)
                end
                elapsed_time >= max_time && (return false)
            end

            elapsed_time += @elapsed begin
                # update
                κ_vio = κ_vio_cand
                r_vio = r_vio_cand

                verbose && iterate_info(j, r, r_vio, κ_vio, Δ, α)
            end
            elapsed_time >= max_time && (return false)
        end

    if (r_vio < r_tol) && (κ_vio < κ_tol)
        elapsed_time += @elapsed begin
            # differentiate solution
            regularization!(ip, κ_vio, ip.opts.γ_reg)
            diff_sol && differentiate_solution!(ip, reg=ip.reg_val[1])
        end
        elapsed_time >= max_time && (return false)
        return true
    else
        return false
    end
end

function rz!(ip::InteriorPoint, rz::AbstractMatrix{T}, z::AbstractVector{T},
    θ::AbstractVector{T}; reg=0.0) where {T}
    z_reg = ip.z_reg
    ortz = ip.idx.ortz
    socz = ip.idx.socz
    z_reg .= z
    for i in eachindex(ortz) # primal-dual
        zi_reg = @view z_reg[ortz[i]]
        zi_reg .= max.(zi_reg, reg)
    end
    ip.methods.rz!(rz, z_reg, θ)
    return nothing
end

function rθ!(ip::InteriorPoint, rθ::AbstractMatrix{T}, z::AbstractVector{T},
        θ::AbstractVector{T}) where {T}
    ip.methods.rθ!(rθ, z, θ)
    return nothing
end

function general_correction_term!(r::AbstractVector{T}, Δ::AbstractVector{T},
    ortr::Vector{Int}, socr::Vector{Int}, socri::Vector{Vector{Int}},
    ortΔ::Vector{Vector{Int}}, socΔ::Vector{Vector{Vector{Int}}}) where {T}
    Δo1 = @views Δ[ortΔ[1]]
    Δo2 = @views Δ[ortΔ[2]]
    rortr = @views r[ortr]
    rortr .+= Δo1 .* Δo2

    num_cone = length(socri)
    for i = 1:num_cone
        Δso1 = @views Δ[socΔ[i][1]]
        Δso2 = @views Δ[socΔ[i][2]]
        rsocri = @views r[socri[i]]
        cone_product!(rsocri, Δso2, Δso1, reset=false)
    end
    return nothing
end

function regularization!(ip::InteriorPoint, κ_vio::T, κ_reg::T, γ_reg::T) where T
    ip.reg_val[1] = κ_vio < κ_reg ? κ_vio * γ_reg : 0.0
    return nothing
end

function regularization!(ip::InteriorPoint, κ_vio::T, γ_reg::T) where T
    reg = κ_vio * γ_reg
    reg > ip.reg_val[1] && (ip.reg_val[1] = reg)
    return nothing
end

function residual_violation(ip::InteriorPoint, r::AbstractVector{T}) where {T}
    req = @views r[ip.idx.equr]
    return norm(req, Inf)
end

function centering(z::AbstractVector{T}, Δaff::AbstractVector{T},
    ortz::Vector{Vector{Int}}, ortΔ::Vector{Vector{Int}},
    socz::Vector{Vector{Vector{Int}}}, socΔ::Vector{Vector{Vector{Int}}},
    zort::Vector{Vector{T}}, Δort::Vector{Vector{T}},
    zsoc::Vector{Vector{Vector{T}}}, Δsoc::Vector{Vector{Vector{T}}},
    αaff::T) where {T}
    # See Section 5.1.3 in CVXOPT
    # μ only depends on the dot products (no cone product)
    # The CVXOPT linear and quadratic cone program solvers
    n = isempty(ortz[1]) ? 0 : length(ortz[1])
    for sc in socz
        n += isempty(sc[1]) ? 0 : 1 # dimension of each soc "variable" is 1
    end

    # ineq
    zo1 = @views z[ortz[1]]
    zo2 = @views z[ortz[2]]
    Δo1 = @views Δaff[ortΔ[1]]
    Δo2 = @views Δaff[ortΔ[2]]
    μ = zo1' * zo2

    Δort[1] .= Δo1
    Δort[1] .*= -αaff
    Δort[1] .+= zo1

    Δort[2] .= Δo2
    Δort[2] .*= -αaff
    Δort[2] .+= zo2

    μaff = Δort[1]' * Δort[2]

    # soc
    num_cone = length(socz)
    for i = 1:num_cone
        zs1 = @views z[socz[i][1]]
        zs2 = @views z[socz[i][2]]
        Δs1 = @views Δaff[socΔ[i][1]]
        Δs2 = @views Δaff[socΔ[i][2]]

        Δsoc[i][1] .= Δs1
        Δsoc[i][1] .*= -αaff
        Δsoc[i][1] .+= zs1

        Δsoc[i][2] .= Δs2
        Δsoc[i][2] .*= -αaff
        Δsoc[i][2] .+= zs2

        μ += zs1' * zs2
        μaff += Δsoc[i][1]' * Δsoc[i][2]
    end
    μ /= n
    μaff /= n
    σ = clamp(μaff / μ, 0.0, 1.0)^3
    return μ, σ
end

function bilinear_violation(ip::InteriorPoint, r::AbstractVector{T}) where {T}
    rbil = @views r[ip.idx.bil]
    return norm(rbil, Inf)
end

function soc_step_length(λ::AbstractVector{T}, Δ::AbstractVector{T}, ρv::AbstractVector{T}, idx::Vector{Int};
    τ::T=0.99, ϵ::T=1e-14, verbose::Bool=false) where {T}
    # check Section 8.2 CVXOPT

    # Adding to slack ϵ to make sure that we never get out of the cone
    λ0 = λ[1] #- ϵ
    λe = @views λ[idx]
    Δe = @views Δ[idx]
    λ_λ = max(λ0^2 - λe' * λe, 1.0e-25)
    λ_λ += ϵ
    λ_Δ = λ0 * Δ[1] - λe' * Δe + ϵ

    ρs = λ_Δ / λ_λ
    ρv .= Δe
    ρv ./= sqrt(λ_λ)
    β = (λ_Δ / sqrt(λ_λ) .+ Δ[1]) / (λ0 / sqrt(λ_λ) + 1)  / λ_λ
    ρv .-= β .* λe
    # # we make sre that the inverse always exists with ϵ,
    # # if norm(ρv) - ρs) is negative (Δ is pushing towards a more positive cone)
    #     # the computation is ignored and we get the maximum value for α = 1.0
    # # else we have α = τ / norm(ρv) - ρs)
    # # we add ϵ to the denumerator to ensure strict positivity and avoid 1e-16 errors.
    α = 1.0
    if norm(ρv) - ρs > 0.0
        α = min(α, τ / (norm(ρv) - ρs))
    end
    return α
end

function soc_step_length(z::AbstractVector{T}, Δ::AbstractVector{T},
    socz::Vector{Vector{Vector{Int}}}, socΔ::Vector{Vector{Vector{Int}}},
    zsoc::Vector{Vector{Vector{T}}}, Δsoc::Vector{Vector{Vector{T}}},
    ρv::Vector{Vector{Vector{T}}}, idx::Vector{Vector{Vector{Int}}};
    τ::T=0.99, verbose::Bool=false) where {T}
    α = 1.0
    num_cone = length(socz)
    for i = 1:num_cone # primal-dual
        for j in eachindex(socz[i]) # number of cones
            # we need -Δ here because we will taking the step x - α Δ
            zsoc[i][j] .= @views z[socz[i][j]]
            Δsoc[i][j] .= @views Δ[socΔ[i][j]]
            Δsoc[i][j] .*= -1.0
            α = min(α, soc_step_length(zsoc[i][j], Δsoc[i][j], ρv[i][j], idx[i][j], τ=τ, verbose=verbose))
        end
    end
    return α
end

function ort_step_length(z::AbstractVector{T}, Δ::AbstractVector{T},
		ortz::Vector{Vector{Int}}, ortΔ::Vector{Vector{Int}};
        τ::T=0.9995) where {T}
    α = 1.0
    num_cone = length(ortz)
    for i = 1:num_cone # primal-dual
        for j in eachindex(ortz[i])
            k = ortz[i][j] # z
            ks = ortΔ[i][j] # Δz
            if Δ[ks] > 0.0
                α = min(α, τ * z[k] / Δ[ks])
            end
        end
    end
    return α
end

function differentiate_solution!(ip::InteriorPoint; reg=0.0)
    s = ip.s
    z = ip.z
    θ = ip.θ
    rz = ip.rz
    rθ = ip.rθ
    δz = ip.δz
    δzs = ip.δzs

    rz!(ip, rz, z, θ, reg = reg)
    rθ!(ip, rθ, z, θ)

    linear_solve!(ip.solver, δzs, rz, rθ, reg=reg)
    δzs .*= -1.0

    mapping!(δz, s, δzs, z)
    return nothing
end

function iterate_info(j, r, r_vio, κ_vio, Δ, α)
    println("iter:", j,
            "  r: ", scn(norm(r, Inf)),
            "  r_vio: ", scn(r_vio),
            "  κ_vio: ", scn(κ_vio),
            "  Δ: ", scn(norm(Δ)),
            "  α: ", scn(norm(α)))
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#polynomial.jl
#solves an arbitrary 6-degree polynomial, using the arbitrary function solver

include("../unum.jl")
using Unums

const plength = 4
#const parray = [exp(100 * randn()) * (randbool() ? -1 : 1) for idx=1:plength]
const parray = [6, 5, -7, 1]
const tol = 0.0001

function pwr(x::Utype, n::Int)
  if n == 1
    x
  elseif iseven(n)
    sqr(pwr(x, n ÷ 2))
  else
    x * pwr(x, n - 1)
  end
end

function polynomial(v::Utype)
  T = isa(v, Ubound) ? typeof(v.lowbound) : typeof(v)
  sum = convert(T, parray[1])
  for (idx = 2:plength)
    sum += convert(T, parray[idx]) * pwr(v, idx - 1)
  end
  sum
end

a = solve(polynomial, tol, 0, 0)#, verbose = true)
map((x) -> println("result:", describe(x)), a)

println()
println()

print("solved polynomial: ")
for idx=plength:-1:2
  print(parray[idx], "x^", (idx - 1), " + ")
end
println(parray[1], " to relative tolerance ", tol)

println("solutions: ")
map((d) -> println(describe(d)), a)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
function [ qvp,qcp,qrainp,thp,VT ] = kessler(n, qvp,qvtot,qcp,qrainp,dt,qc0,k1,k2,cp,rd,lv,p0,pib,thp,pip,rhou,rhow,tb,qb,VT,nx,nz )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

%% Rain Water Conversions

for k=2: nz-1
    for i=2: nx-1

      autodt=0.0;
      accrdt=0.0;
      revapdt=0.0;

%  AUTOCONVERSION OF CLOUD TO RAIN

        if(qcp[i,k]>0.0)
          auto = k1*max(qcp[i,k]-qc0,0.0);
          autodt = auto*2*dt;
        end

% ACCRETION OF CLOUD BY RAIN

        if qcp[i,k]>0.0 && qrainp[i,k]>0.0
          accr = rhou[k]*k2*qcp[i,k]*(qrainp[i,k]^0.875);
          accrdt = accr*2*dt;
        end

% limit amount of cloud water we can lose

        if autodt+accrdt>qcp[i,k]
          qcsink=autodt+accrdt;
          autodt=qcp[i,k]*autodt/qcsink;
          accrdt=qcp[i,k]*accrdt/qcsink;
        end

% EVAPORATION OF RAIN

        if(qrainp[i,k]> 0.0)
          cvent = 1.6 + 30.39*(rhou[k]*qrainp[i,k])^0.2046;
          pprime = pip[i,k]*cp*rhou[k]*tb[k];
          pbar = p0*pib[k]^(cp/rd);
          pres = pbar + pprime;
          temp = (thp[i,k]+tb[k])*(pip[i,k]+pib[k]);
% Teten's formula
          qvsat = (380.0/pres) *                       ...
                   exp((17.27*(temp-273.0))/(temp-36.0));

          if qvsat>(qb[k]+qvp[i,k])

            revap = cvent *                                           ...
                   (1.0-((qvp[i,k]+qb[k])/qvsat))*                    ...
                   ((rhou[k]*qrainp[i,k])^0.525)/                       ...
                   ( (2.03e4+9.584e6/(qvsat*pres)) * rhou[k] );

%  limit amount of rain we can evaporate

            revapdt = min(revap*2*dt,qrainp[i,k]);
            if (qvp[i,k]+qb[k]+revapdt)>qvsat
              revapdt=qvsat -(qvp[i,k]+qb[k]);
            end

          end
        end

        qvp[i,k] = qvp[i,k] + revapdt;
        qcp[i,k] = qcp[i,k] - (autodt+accrdt);
        qrainp[i,k] = qrainp[i,k] + (autodt+accrdt) - revapdt;
        thp[i,k] = thp[i,k] - revapdt*lv/(cp*pib[k]);

    end
end

%% Compute Rainwater fallsppeds for each grid box
%(to be used in qr next time step)

for k=1: nz
    for i=1: nx

        VT[i,k] = 0.0;
        if qrainp[i,k]>1.0*10^(-12)
            qrr = max(qrainp[i,k]*0.001*rhou[k],0.0);
            vtden = sqrt(rhou(1)/rhou[k]);
            VT[i,k] = 36.34*(qrr^0.1364)*vtden;
        end

    end
end

%% Computation of condensation/evaporation of cloud mass
%

% CONDENSATION/EVAPORATION OF CLOUD MASS
      for k= 2 : nz-1
      for i= 2 : nx-1
        pprime = pip[i,k]*cp*rhou[k]*tb[k];
        pbar = p0*pib[k]^(cp/rd);
        pres = pbar + pprime;
        temp = (thp[i,k]+tb[k])*(pip[i,k]+pib[k]);
% Teten's formula
        qvsat = (380.0/pres) *                       ...
                   exp((17.27*(temp-273.0))/(temp-36.0));

        if qcp[i,k] > 0 || (qvp[i,k]+qb[k] > qvsat)

%  compute the adjustment to total evap/cond owing to slight change in
%    qvsat by the latent heating/chilling
          phi = qvsat*lv*17.27*237.0/(cp * (temp-36.0)^2);
          dqv = (qvsat-qvp[i,k]-qb[k])/(1+phi);

%  if supersaturated, dqv is negative; if subsaturated, dqv is positive;
%     for the subsaturated case, only allow evaporation up to qc, and no more
          dqv = min(dqv,qcp[i,k]);

          qvp[i,k] = qvp[i,k] + dqv;
          qcp[i,k] = qcp[i,k] - dqv;
          thp[i,k] = thp[i,k] - dqv*lv/(cp*pib[k]);

        end

      end
      end

end

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
using GeoStats
using Test, Pkg, Random

# list of tests
testfiles = [
]

@testset "GeoStats.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
## _PARSEDIM
#
# Log:
# 11.05.2021 - Created function [1]
#
# Authors:
# [1] Alexander Smolka (alexander.smolka@tum.de)

"""
    _parsedim(dim::String)
"""
function _parsedim(dim::String)
    if dim == "2Dsym"
        return 2, ("R", "Z")
    end
end


text:
using Mocking

describe_report_definitions_1_patch = @patch function describe_report_definitions(a...; b...)
    return Dict{String, Any}("ReportDefinitions" => Any[Dict{String, Any}("Format" => "textORcsv", "AdditionalArtifacts" => Any["QUICKSIGHT"], "ReportName" => "ExampleReport", "AdditionalSchemaElements" => Any["RESOURCES"], "Compression" => "GZIP", "S3Region" => "us-east-1", "TimeUnit" => "HOURLY", "S3Prefix" => "exampleprefix", "S3Bucket" => "example-s3-bucket"), Dict{String, Any}("Format" => "textORcsv", "AdditionalArtifacts" => Any["QUICKSIGHT"], "ReportName" => "ExampleReport2", "AdditionalSchemaElements" => Any["RESOURCES"], "Compression" => "GZIP", "S3Region" => "us-east-1", "TimeUnit" => "HOURLY", "S3Prefix" => "exampleprefix", "S3Bucket" => "example-s3-bucket")])
end