text:
# function : current state s, current action a, state space S maps to next state s2
function simplechain_transition(s,a,S)
    if a == 1  # go left
        s2 = s-1 # decrease state value by 1
    elseif a == 2 # go right
        s2 = s+1 # increase state value by 1
    elseif a == -1
        s2 = s
    else
        error("action must be {1,2}, but was : ", a)
    end
    if s2 > maximum(S)  # if s2 is larger than the maximum state
        s2 = minimum(S) # set it to the start state (smallest state)
    elseif s2 < minimum(S) # if s2 is smaller than smallest state
        s2 = minimum(S) # make it the smallest state
    end
    return s2  # return next state
end

# function : current state s, current action a, state space S maps to next state s2
function hardchain_transition(s,a,S)
    if a == 1  # go left
        s2 = mimimum(S) # decrease state value by 1
    elseif a == 2 # go right
        s2 = s+1 # increase state value by 1
    elseif a == -1
        s2 = s
    else
        error("action must be {1,2}, but was : ", a)
    end
    if s2 > maximum(S)  # if s2 is larger than the maximum state
        s2 = minimum(S) # set it to the start state (smallest state)
    elseif s2 < minimum(S) # if s2 is smaller than smallest state
        s2 = minimum(S) # make it the smallest state
    end
    return s2  # return next state
end

function simplechain_perturb_action(a::Int, failchance)
    ϵ = rand()
    if ϵ ≤ 1.0 - failchance
        return a
    elseif ϵ ≤ 1.0 - failchance / 2
        return Int(-1)
    else
        return a==1 ? 2 : 1
    end
end

function create_simple_chain(num_states::Int; stochastic=false, failchance=0.1, hard_transition=false)
    S = 1:num_states
    A = 1:2

    if hard_transition
        tfun = hard_transition
    else
        tfun = simplechain_transition
    end

    if !stochastic
        p = (s,a)->tfun(s,a,S)  # mask out S in above function
    else
        p = (s,a)->tfun(s,simplechain_perturb_action(a,failchance), S)
    end
    d0 = ()->minimum(S)
    r = (s,a,s′)-> -1.0
    γ = (s,a,s′)-> (s==num_states && s′==1) ? 0.0 : 1.0
    meta = Dict{Symbol,Any}()
    meta[:minreward] = -1.0
    meta[:maxreward] = -1.0
    meta[:minreturn] = -Inf
    meta[:maxreturn] = -Float64(num_states)
    meta[:stochastic] = stochastic
    meta[:minhorizon] = num_states
    meta[:maxhorizon] = Inf
    meta[:discounted] = false
    m = MDP(S,A,p,r,γ,d0,meta, ()->nothing)
    return m
end

function simple_chain(num_states::Int; stochastic=false, failchance=0.1, randomize=false)
    if randomize
        failchance = rand() * 0.25
    end
    return create_simple_chain(num_states, stochastic=stochastic, failchance=failchance)
end

function finite_horizon_transition(m0, s, a, maxT)
    t = s[1]
    t += 1
    snum = s[2]
    s′ = m0.p(snum,a)
    if t > maxT || m0.γ(snum, a, s′)==0.0
        t = 1
        s′ = m0.d0()
    end

    return (t, s′)
end

function finitehorizon_chainstep(p, r, s, a, maxT, S)
    t,x = s
    t += 1
    x′ = p(x,a)
    maxS = maximum(S)
    minS = minimum(S)
    γ = 1.0
    if t > maxT || (x==maxS && x′==minS)
        γ = 0.0
    end
    reward = r(x,a,x′)

    return t,x′,reward,γ
end

function chain_constreward(s,a,s′)
    return -1.0
end

function chain_scaledreward(s,a,s′,k)
    return -1.0 * s/k
end

function create_simple_chain_finitetime(num_states::Int; stochastic=false, failchance=0.1, droptime=true, scale_reward=false)
    maxT = num_states * 20
    S = (1:maxT, 1:num_states)
    A = 1:2
    meta = Dict{Symbol,Any}()

    rfun = (s,a,s′)->-1.0#chain_scaledreward(s,a,s′,1.0)#Float64(-1.0)
    meta[:minreward] = -1.0
    meta[:maxreward] = -1.0
    # if scale_reward
    #     # rfun = (s,a,s′)-> Float64(-1.0) / (Float64(s[2]) / Float64(num_states))
    #     rfun = (s,a,s′)->chain_scaledreward(s,a,s′,Float64(num_states))
    #     meta[:minreward] = -1.0
    #     meta[:maxreward] = -1.0 / num_states
    # end

    tfun = (s,a)->simplechain_transition(s,a,S[2])

    if !stochastic
        pfun = (s,a)->finitehorizon_chainstep(tfun, rfun, s, a, maxT, S[2])
    else
        pfun = (s,a)->finitehorizon_chainstep(tfun, rfun, s, simplechain_perturb_action(a,failchance), maxT, S[2])
    end

    if droptime
		X = S[2]
        # function get_outcome1(pfun,s,a)
        function get_outcome1(s,a,maxT,S)
            # @code_warntype (pfun(s,a))
            t, x, r, γ = finitehorizon_chainstep(tfun, rfun, s, a, maxT, S)
            # t, x, r, γ = pfun(s,a)
			s = (t,x)
            return s,x,r,γ
            # return s,s[2],-1.0,1.0
		end
        p = (s,a)->get_outcome1(s, a, maxT, S[2])
        d0 = ()->((1,1),1)
	else
		X = S
		function get_outcome2(pfun,s,a)
			t, x, r, γ = pfun(s,a)
			s = (t,x)
			return s,s,r,γ
		end
        p = (s,a)->get_outcome2(pfun, s, a)
        d0 = ()->((1,1),(1,1))
	end

    meta[:minreward] = -1.0
    meta[:maxreward] = -1.0
    meta[:minreturn] = -Float64(maxT)
    meta[:maxreturn] = -Float64(num_states)
    meta[:stochastic] = stochastic
    meta[:minhorizon] = num_states
    meta[:maxhorizon] = maxT
    meta[:discounted] = false

    render = (state,clearplot=false)->nothing

	m = SequentialProblem(S,X,A,p,d0,meta,render)

    return m
end

function simple_chain_finitetime(num_states::Int; stochastic=false, failchance=0.1, randomize=false, droptime=true, scale_reward=false)
    if randomize
        failchance = rand() * 0.25
    end
    return create_simple_chain_finitetime(num_states, stochastic=stochastic, failchance=failchance, droptime=droptime, scale_reward=scale_reward)
end

function create_hard_chain_finitetime(num_states::Int; stochastic=false, failchance=0.1, droptime=true, scale_reward=false)
    maxT = num_states * 20
    S = (1:maxT, 1:num_states)
    A = 1:2
    m0 = create_simple_chain(num_states; stochastic=stochastic, failchance=failchance, hardchain_transition)
    p = (s,a)->finite_horizon_transition(m0, s, a, maxT)
    γ = (s,a,s′)-> (s[2]==num_states && s′[2]==1 && a==2) || (s[1] ≥ s′[1]) ? 0.0 : 1.0

    meta = Dict{Symbol,Any}()

    if scale_reward
        r = (s,a,s′)->m0.r(s[2],a,s′[2]) / (s[2] / num_states)
        meta[:minreward] = -1.0
        meta[:maxreward] = -1.0 / num_states
    else
        r = (s,a,s′)->m0.r(s[2],a,s′[2])
        meta[:minreward] = -1.0
        meta[:maxreward] = -1.0
    end
    d0 = ()->(1,m0.d0())

    meta[:minreturn] = -Float64(maxT)
    meta[:maxreturn] = -Float64(num_states)
    meta[:stochastic] = stochastic
    meta[:minhorizon] = num_states
    meta[:maxhorizon] = maxT
    meta[:discounted] = false

    if droptime
        X = S[2]
        obs = s->s[2]
        m = POMDP(S,A,X,p,obs,r,γ,d0,meta, ()->nothing)
    else
        m = MDP(S,A,p,r,γ,d0,meta, ()->nothing)
    end

    return m
end

function hard_chain_finitetime(num_states::Int; stochastic=false, failchance=0.1, randomize=false, droptime=true, scale_reward=false)
    if randomize
        failchance = rand() * 0.25
    end
    return create_hard_chain_finitetime(num_states, stochastic=stochastic, failchance=failchance, droptime=droptime, scale_reward=scale_reward)
end


text:
@testset "Comparison" begin
    f0,f1,f2,f3 = -0.9,-1.0,0.0,1.0

    @testset for T in (Posit8, Posit16, Posit32, Posit24_2)
        @test (f0 > f1) == (T(f0) > T(f1))
        @test (f2 >= f2) == (T(f2) >= T(f2))
        @test (f3 < f1) == (T(f3) < T(f3))
        @test T(f0) === T(f0)
    end
end

ext: jl


text:
# TODO: Consider other concept for the names here instead of mapping (ex. partitioning)
module Partitioners

export AbstractPartitioner, partition

abstract type AbstractPartitioner end

struct PartitionerIterator{P <: AbstractPartitioner,T <: AbstractVecOrMat}
    partitioner::P
    X::T
end

# In these generic functions I'm implicitly enforcing that all
# AbstractPartitioner subtypes have `width` and `n` fields. Maybe it would be
# better to enforce in a more explicit manner through traits

function partition(partitioner::P, X::T) where {P <: AbstractPartitioner,T <: AbstractVector}
    length(X) != partitioner.width && throw(
      DimensionMismatch("Expected length of X to be $(partitioner.width), got $(length(X))"))

    return PartitionerIterator{P,T}(partitioner, X)
end

function partition(partitioner::P, X::T) where {P <: AbstractPartitioner,T <: AbstractMatrix}
    size(X, 2) != partitioner.width && throw(
    DimensionMismatch("Expected the number of columns of X to be $(partitioner.width), got $(size(X, 2))"))

    return PartitionerIterator{P,T}(partitioner, X)
end

Base.length(partitioner::P) where {P <: AbstractPartitioner} = ceil(Int, partitioner.width / partitioner.n)

function Base.length(itr::PartitionerIterator{P,T}) where {P <: AbstractPartitioner,T <: AbstractVecOrMat}
    # ceil(Int, size(itr.X)[end] / itr.partitioner.n)
    length(itr.partitioner)
end

Base.eltype(::Type{PartitionerIterator{P,T}}) where {P <: AbstractPartitioner,T <: AbstractVecOrMat} = T

include("RandomPartitioner.jl")
export RandomPartitioner, random_partitioning

include("LinearPartitioner.jl")
export LinearPartitioner
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
using WebDataset
using Test

@testset "WebDataset.jl" begin
    sampletar = "sample.tar"

    # generic_open
    @test read(generic_open("pipe:echo hello"), String) == "hello\n"
    @test read(generic_open("/dev/null"), String) == ""

    # braceexpand
    @test braceexpand("") == [""]
    @test braceexpand("abc") == ["abc"]
    @test braceexpand("a{1..3}bc") == ["a1bc", "a2bc", "a3bc"]
    @test braceexpand("a{1..2}{1..2}bc") == ["a11bc", "a12bc", "a21bc", "a22bc"]

    # tar file iterator
    @test tariterator(open(sampletar)) |> collect |> length == 180

    # sample iteration
    @test sampleiterator(open(sampletar)) |> collect |> length == 90
    @test sampleiterator(open(sampletar)) |> sampleshuffle |> collect |> length == 90

    # transforms
    l = sampleiterator(open(sampletar)) |> it->sampletransforms(it, default_decoders) |> collect
    @test length(l) == 90
    first = l[1][".img"]
    @test size(first) == (28, 28)

    # batching
    l = sampleiterator(open(sampletar)) |> it->samplebatching(it, 10) |> collect
    @test length(l) == 9
    first = l[1][".png"]
    @test length(first) == 10

    # data loader
    ds = DatasetDescriptor(
        sources=[sampletar, sampletar],
        shuffle=10,
        batchsize=10,
        decoding=default_decoders,
        ntasks=2,
    )
    l = dataloader(ds) |> collect
    @test length(l) == 18
    @test length(l[1][".img"]) == 10
    @test size(l[1][".img"][1]) == (28, 28)
end


text:
struct BridgeRotation{T}
    buff::Vector{T}
end

export BridgeRotation

"""
    rotate!(bridge::BridgeRotation, A::AbstractVector, tailpos::Integer)

Rotates `A` on `tailpos` using a small extra memory. ≤ a third of A and ≤ than n+n/3 copy-write ops.
"""
function rotate!(bridge::BridgeRotation, A::AbstractVector, tailpos::Integer)
    m = tailpos - 1
    n = length(A) - m
    if m <= n
        _interchange_right_large!(A, bridge.buff, tailpos)
    else
        _interchange_left_large!(A, bridge.buff, tailpos)
    end

    A
end

"""
    _interchange_left_large!(A::AbstractVector, buff::AbstractVector, tail::Integer)

Interchange and shift two subarrays

```

      buffer
      -----
4 5 6 7 8 9 1 2 3
-----       |
head        tail
```

must produce
```
      prev head
      -----
1 2 3 4 5 6 1 2 3
-----       -----
prev tail   untouched
```

and finally, the buffer should be copied-back
```
      prev head
      -----
1 2 3 4 5 6 7 8 9
-----       -----
prev tail   from buffer
```
"""
@inline function _interchange_left_large!(A::AbstractVector, buff::AbstractVector, tail::Integer)
    m = length(A)
    n = m - tail + 1  # number of elements
    p = 2n
    resize!(buff, m - p)
    _copy!(buff, 1, A, n + 1, length(buff))
    @inbounds for i in 1:n
        n += 1
        A[n] = A[i]
        A[i] = A[tail]
        tail += 1
    end
    _copy!(A, p + 1, buff, 1, length(buff))
    A
end

"""
    _interchange_right_large!(A::AbstractVector, buff::AbstractVector, tail::Integer)

Interchange and shift two subarrays

```

head        buffer
-----       -----
7 8 9 1 2 3 4 5 6
      |
      tail
```

must produce
```
prev tail   prev head
-----       -----
1 2 3 1 2 3 7 8 9
      -----
     untouched
```

The centering block is then filled with the buffer.
"""
@inline function _interchange_right_large!(A::AbstractVector, buff::AbstractVector, tail::Integer)
    m = length(A)
    n = tail - 1  # number of elements
    p = 2n
    resize!(buff, m - p)
    _copy!(buff, 1, A, p + 1, length(buff))

    @inbounds for i in n:-1:1
        A[m] = A[i]
        A[i] = A[p]
        p -= 1
        m -= 1
    end
    _copy!(A, n + 1, buff, 1, length(buff))
    A
end


text:
numbers = [0,14,1,3,7,9]
# numbers = [1,3,2]
# numbers = [2,1,3]

function get_next_number!(spoken)
    idx = findlast(==(spoken[end]), spoken[1:end-1])
    if idx === nothing
        push!(spoken, 0)
    else
        push!(spoken, length(spoken) - idx)
    end
end

spoken_history = [numbers...]
for i in length(numbers)+1:2020
    global spoken_history
    get_next_number!(spoken_history)
end

println(spoken_history[end])


text:
using Statistics, LinearAlgebra, PyPlot

mutable struct Lattice
    name::String
    state::Array{Array{Float64,3},1}
    probability::Array{Array{Float64,3},1}
    Lattice(name::String, state::Array{Array{Float64,3},1}, probability::Array{Array{Float64,3},1}) = new(name, state, probability)
end

"""
	lattice_approximation(bstructure::Array{Int64,1}, path::Function, nIterations::Int64, r::Int64 = 2, dim::Int64 = 1)
Returns a valuated approximated lattice for the stochastic process provided in any dimension. The default dimension is dim = 1 for a lattice in 1 dimension.

Args:
- bstructure - Branching structure of the scenario lattice e.g., bstructure = [1,2,3,4,5] represents a 5-staged lattice
- path - Function generating samples from a known distribution with length equal to the length of bstructure of the lattice.
- nIterations - Number of iterations to be performed.
- r - Parameter for the transportation distance.
- dim - dimension of the lattice to be generated. This depends entirely on the dimension of the samples generated.
"""
function lattice_approximation(bstructure::Array{Int64,1}, path::Function, nIterations::Int64, r::Int64 = 2, dim::Int64 = 1)
    tdist = zeros(dim)         # multistage distance
    T = length(bstructure)     # number of stages in the lattice
    states = [zeros(bstructure[j], 1, dim) for j = 1 : T] # States of the lattice at each time t
    probabilities = vcat([zeros(bstructure[1], 1, dim)],[zeros(bstructure[j-1], bstructure[j], dim) for j = 2 : T]) # Probabilities of the lattice at each time t
    init_path = path()
    #Check the dimension of the sample path if it is equal to the dimension of the lattice.
    if dim != size(init_path)[2]
        @error("Dimension of lattice ($dim) is not equal to the dimension of the input array ($(size(init_path)[2]))")
    end
    # Initialize the states of the nodes of the lattice
    for t = 1 : T
        states[t] .= init_path[t]
    end
    Z = Array{Float64}(undef, T, dim) # Array to hold the new samples generated
    #Stochastic approximation step starts here
    for n = 1 : nIterations
        Z .= path() # Draw a new sample Gaussian path
        last_index = Int64.(ones(dim))
        dist = zeros(dim)
        for t = 1 : T
            for i = 1 : dim
                # corrective action to include lost nodes
                states[t][:, :, i][findall(sum(probabilities[t][:, :, i], dims = 2) .< 1.3 * sqrt(n) / bstructure[t])] .= Z[t,i]
                min_dist, new_index = findmin(abs.(vec(states[t][:, :, i] .- Z[t,i]))) # find the closest lattice entry
                dist[i] += min_dist^2  # Euclidean distance for the paths
                probabilities[t][last_index[i], new_index, i] += 1.0  # increase the counter for the nodes
                # use stochastic approximation algorithm to calculate for the new states of the nodes in the lattice.
                states[t][new_index, :, i] = states[t][new_index, :, i] - r * (min_dist)^(r-1) * (states[t][new_index, :, i] .- Z[t,i]) / ((Z[1] + 30 + n))
                last_index[i] = new_index
            end
        end
        # calculate the multistage distance
        dist = dist .^ (1/2); tdist = (tdist .* (n - 1) + dist .^r )/ n
        #tdist = (tdist*(n-1) + dist)/n
    end
    # calculate the probabilities
    probabilities = probabilities ./ nIterations
    # functions to print the results of the lattices
    # these functions depends on the dimension of the lattice
    if dim == 1
        # Round the results to make sure that they are neat.
        # Remember it is a 3D array, so a loop will do.
        return Lattice("Approximated Lattice with $bstructure branching structure and distance=$(tdist .^ (1 / r)) at $(nIterations) iterations",
        [round.(states[i], digits = 6) for i = 1 : length(states)], [round.(probabilities[j], digits = 6) for j = 1 : length(probabilities)])
    else
        lattices = Lattice[] # create an empty array of type Lattice to hold the resulting lattices
        for i = 1 : length(states[1])
            st = [zeros(bstructure[j], 1, 1) for j = 1 : length(states)]
            pp = vcat([zeros(bstructure[1], 1, 1)],[zeros(bstructure[j-1] , bstructure[j], 1) for j = 2 : length(states)])

            for j = 1 : length(states)
                st[j] .= states[j][:, :, i]
                pp[j] .= probabilities[j][:, :, i]
            end
            sublattice = Lattice("Lattice $i with distance=$((tdist.^(1/r))[i])",
            [round.(st[i], digits = 6) for i = 1 : length(st)], [round.(pp[j], digits = 6) for j = 1 : length(pp)])
            push!(lattices,sublattice)
        end
        return lattices
    end
end

"""
	PlotLattice(lt::Lattice,fig=1)

Returns a plot of a lattice.
"""
function plot_lattice(lt::Lattice,fig = 1)
    if !isempty(fig)
        figure(figsize=(6,4))
    end
    lts = subplot2grid((1,4),(0,0),colspan = 3)
    title("states")
    xlabel("stage,time",fontsize=11)
    #ylabel("states")
    xticks(1:length(lt.state))
    lts.spines["top"].set_visible(false)                                                         # remove the box at the top
    lts.spines["right"].set_visible(false)                                                       # remove the box at the right
    for t = 2:length(lt.state)
        for i=1:length(lt.state[t-1])
            for j=1:length(lt.state[t])
                lts.plot([t-1,t],[lt.state[t-1][i],lt.state[t][j]])
            end
        end
    end

    prs = subplot2grid((1,4), (0,3))
    title("probabilities")
    prs.spines["top"].set_visible(false)
    prs.spines["left"].set_visible(false)
    prs.spines["right"].set_visible(false)
# Use the states and probabilites at the last stage to plot the marginal distribution
    stts = lt.state[end]
    n = length(stts)                                    # length of terminal nodes of the lattice.
    h = 1.05 * std(stts) / (n^0.2) + eps()                  #Silverman rule of thumb
    #lts.set_ylim(minimum(stts)-1.5*h, maximum(stts)+1.5*h)
    #prs.set_ylim(minimum(stts)-3.0*h, maximum(stts)+3.0*h)
    proba = sum(lt.probability[end], dims=1)
    yticks(())                                          #remove the ticks on probability plot
    t = LinRange(minimum(stts) - h, maximum(stts) + h, 100) #100 points on probability plot
    density = zero(t)
    for (i, ti) in enumerate(t)
        for (j, xj) in enumerate(stts)
            tmp = (xj - ti) / h
            density[i] += proba[j]* 35/32 * max(1.0 -tmp^2, 0.0)^3 / h #triweight kernel
        end
    end
    plot(density, t)
    prs.fill_betweenx(t, 0 , density)
end


text:
###############################################################################
#
#   Random Functions
#
###############################################################################

# mainly for testing
function rand(L::Loc{T}, num_scale = (1:1000), den_scale=(1:1000)) where {T <: fmpz}
   num = rand(num_scale)
   den = rand(den_scale)
   while gcd(den,prime(L)) != 1
      den = rand(den_scale)
   end
   return L(num//den)
end

Nemo.promote_rule(::Type{LocElem{T}}, ::Type{T}) where {T} = LocElem{T}


text:
#using LinearAlgebra
#using DelimitedFiles
#using Statistics
#using Plots
#using DelimitedFiles
#using Random
#using PyCall
#using SparseArrays

#plt = pyimport("matplotlib.pyplot")

#include("Geometry_twisted.jl")

### for triangular lattice, use pi-flux hopping!!!

struct get_hamiltonian
    intra
    tx
    ty
    txy
    geometry::geometry_twisted
end

#for twisted lattice, number of sites very big, need sparse matrix, to avoid reallocation of data, create it from [row,col,val], and preassign large enough length to row/col/val.
function get_H_twisted_bilayer(R,r,t_xy,t_z,d,lambda,mode,theta,r0,interlayer_bias=0)  # d: inter-layer distance, lambda: parameter, mode="uniform_hop" or "pi-flux"
    t0=now()
    N=size(R,1)

    #H=spzeros(Complex,N,N)                        # H_QSL is spin-degen
    I,J,value=ones(Int,N*100),ones(Int,N*100),zeros(Complex,N*100)
    n=1
    for i=1:N
        if R[i][3]!=0
           I[n]=i
           J[n]=i
           value[n]=interlayer_bias/2
           n+=1
        else
           I[n]=i
           J[n]=i
           value[n]=-interlayer_bias/2
           n+=1
        end

        for j=1:N
            dR=R[i]-(R[j]+r)
            z=dR[3]
            L=norm(dR)

            if abs(z)>0.01
                if abs(t_z*(z^2/L^2)*exp(-lambda*(L-d)))>0.0001
                    I[n]=i
                    J[n]=j
                    value[n]=t_z*(z^2/L^2)*exp(-lambda*(L-d))
                    n+=1
                    #H[i,j]=H[i,j]+t_z*(z^2/L^2)*exp(-lambda*(L-d))         # from j to i
                end
            elseif 0.9<L<1.1
                t=get_hopping(R[i],R[j]+r,t_xy,mode,theta,r0)
                I[n]=i
                J[n]=j
                value[n]=t
                n+=1
            end
        end
    end
    H=sparse(I,J,value,N,N)
    H=dropzeros(H)
    t1=now()
    println("time to construct H: ",t1-t0)
    return H
end

# for H_inter, H.tx+ty is always 0,only H.tx H.ty and H.tx-ty is nonzero, because theta(a1,a2)=60degree
function get_H_inter_twisted_bilayer(g::geometry_twisted, t_xy,t_z,lambda,mode,interlayer_bias=0,theta2=0)    #theta2: additional rotation of second layer w.r.t. first layer, can be n*2pi/3
    R=g.sites
    d=g.inter_distance
    theta=g.twist_angle
    theta+=theta2
    r0=g.inter_alignment
    H_intra=get_H_twisted_bilayer(R,[0,0,0],t_xy,t_z,d,lambda,mode,theta,r0,interlayer_bias)
    H_tx=get_H_twisted_bilayer(R,g.inter_vector.x,t_xy,t_z,d,lambda,mode,theta,r0)
    H_ty=get_H_twisted_bilayer(R,g.inter_vector.y,t_xy,t_z,d,lambda,mode,theta,r0)
    H_txy=get_H_twisted_bilayer(R,g.inter_vector.xy,t_xy,t_z,d,lambda,mode,theta,r0)
    H=get_hamiltonian(H_intra,H_tx,H_ty,H_txy,g)
    return H
end

function get_Hk(H::get_hamiltonian,k::Array{Float64,1})
    H0=exp(im*dot(k,H.geometry.inter_vector.x))*H.tx+exp(im*dot(k,H.geometry.inter_vector.y))*H.ty+exp(im*dot(k,H.geometry.inter_vector.xy))*H.txy
    Hk=H.intra+H0+adjoint(H0)
    return Hk
end

function get_Hk_ribbon(H::get_hamiltonian,k::Array{Float64,1})   # here k represents the phase difference
    @assert k[1]*k[2]==0 "please input 1D momenta"
    k[1]==0 ? H0=exp(im*k[2])*H.ty : H0=exp(im*k[1])*H.tx

    Hk=H.intra+H0+adjoint(H0)
    return Hk
end

function get_hopping(R,R2,t_xy,mode,theta,r0)
    factor=1
    if mode=="uniform_hop"
        if R[3]==1
            factor=1
        end
        return t_xy*factor
    elseif mode=="pi_flux"
        U=rotation_along_z_3D(-(theta))
        if R[3]==1
            R=U*(R-r0)
            R2=U*(R2-r0)
            factor=1
        end
        x=round(R[1]-R2[1],digits=2)
        y=R[2]-R2[2]
        if x*y<-0.01                             # with this method, it is better that the original lattice passes (0,0)
            hop=t_xy
        elseif x>0.1&&y>0.1
            if mod(round(-R[2]*2/sqrt(3)),2)==0  # from B to A, R==A
                hop=-t_xy
            else
                hop=t_xy
            end
        elseif x<-0.1&&y<-0.1
            if mod(round(-R[2]*2/sqrt(3)),2)==0
                hop=t_xy
            else
                hop=-t_xy
            end
        else
            if mod(round(-R[2]*2/sqrt(3)),2)==0
                hop=t_xy
            else
                hop=-t_xy
            end
        end
        return hop*factor
    else
        println("invalid mode, default to uniform_hop")
        return t_xy
    end
end

function plot_H(H,R)
    plt.figure(figsize=(6,6),dpi=80)
    N=size(H,1)
    for i=1:N
        for j=1:N
            if real(H[i,j])>0.1
                plt.plot((R[i][1],R[j][1]), (R[i][2],R[j][2]),color="red",linewidth=norm(H[i,j])*1)
            elseif real(H[i,j])<-0.1
                plt.plot((R[i][1],R[j][1]), (R[i][2],R[j][2]),color="blue",linewidth=norm(H[i,j])*1)
            end
        end
    end
    plt.axis([-5,5,-5,5])
    plt.xlabel("")
    plt.ylabel("")
    plt.show()
end

function plot_H_2(H_inter,R,inter_vector)
    plt.figure(figsize=(6,6),dpi=80)
    n=size(H_inter,1)

    for i=1:n
        H=H_inter[i]
        r=inter_vector[i]
        R2=[x+r for x in R]

        N=size(H,1)
        for i=1:N
            for j=1:N
                if real(H[i,j])>0.1
                    plt.plot((R[i][1],R2[j][1]), (R[i][2],R2[j][2]),color="red",linewidth=norm(H[i,j])*1)
                elseif real(H[i,j])<-0.1
                    plt.plot((R[i][1],R2[j][1]), (R[i][2],R2[j][2]),color="blue",linewidth=norm(H[i,j])*1)
                end
            end
        end
    end
    plt.axis([-5,5,-5,5])
    plt.xlabel("")
    plt.ylabel("")
    plt.show()
end

function test_tight_binding_twisted()
    n_a1=50
    n_a2=50
    m=1
    r=1
    r0=[0,0,0]
    d=1
    lattice,mode="triangular","pi_flux"
    BC="PBC"

    g=get_twisted_lattice_2(lattice,n_a1,n_a2,m,r,d,r0)
    println(g.twist_angle*180/pi)
    println(size(g.sites,1))

    t_xy=1
    t_z=0
    lambda=10
    H=get_H_inter_twisted_bilayer(g,t_xy,t_z,lambda,mode,0,pi)
    Hk=get_Hk(H,[0.0,0.0,0.0])
    plot_H(Hk,g.sites)
end

#test_tight_binding_twisted()


text:
function create_param_dual_cum_sum(model::ParametricOptimizer)
    return zeros(model.number_of_parameters_in_model)
end

function calculate_dual_of_parameters(model::ParametricOptimizer)
    param_dual_cum_sum = create_param_dual_cum_sum(model)
    update_duals_with_affine_constraint_cache!(param_dual_cum_sum, model)
    update_duals_with_quadratic_constraint_cache!(param_dual_cum_sum, model)
    update_duals_in_affine_objective!(
        param_dual_cum_sum,
        model.affine_objective_cache,
    )
    update_duals_in_quadratic_objective!(
        param_dual_cum_sum,
        model.quadratic_objective_cache_pc,
    )
    empty_and_feed_duals_to_model(model, param_dual_cum_sum)
    return model
end

function update_duals_with_affine_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    model::POI.ParametricOptimizer,
)
    for S in SUPPORTED_SETS
        affine_constraint_cache_inner =
            model.affine_constraint_cache[MOI.ScalarAffineFunction{Float64}, S]
        if !isempty(affine_constraint_cache_inner)
            update_duals_with_affine_constraint_cache!(
                param_dual_cum_sum,
                model.optimizer,
                affine_constraint_cache_inner,
            )
        end
    end
    return
end

function update_duals_with_affine_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    optimizer::OT,
    affine_constraint_cache_inner::DD.WithType{F,S},
) where {OT,F,S}
    for (ci, param_array) in affine_constraint_cache_inner
        calculate_parameters_in_ci!(
            param_dual_cum_sum,
            optimizer,
            param_array,
            ci,
        )
    end
    return
end

function update_duals_with_quadratic_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    model::POI.ParametricOptimizer,
)
    for S in SUPPORTED_SETS
        quadratic_constraint_cache_pc_inner =
            model.quadratic_constraint_cache_pc[
                MOI.ScalarQuadraticFunction{Float64},
                S,
            ]
        if !isempty(quadratic_constraint_cache_pc_inner)
            update_duals_with_quadratic_constraint_cache!(
                param_dual_cum_sum,
                model,
                quadratic_constraint_cache_pc_inner,
            )
        end
    end
    return
end

function update_duals_with_quadratic_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    model::ParametricOptimizer,
    quadratic_constraint_cache_pc_inner::DD.WithType{F,S},
) where {F,S}
    for (poi_ci, param_array) in quadratic_constraint_cache_pc_inner
        moi_ci = model.quadratic_added_cache[poi_ci]
        calculate_parameters_in_ci!(
            param_dual_cum_sum,
            model.optimizer,
            param_array,
            moi_ci,
        )
    end
    return
end

function calculate_parameters_in_ci!(
    param_dual_cum_sum::Vector{Float64},
    optimizer::OT,
    param_array::Vector{MOI.ScalarAffineTerm{T}},
    ci::CI,
) where {OT,CI,T}
    cons_dual = MOI.get(optimizer, MOI.ConstraintDual(), ci)

    for param in param_array
        param_dual_cum_sum[param.variable_index.value-PARAMETER_INDEX_THRESHOLD] +=
            cons_dual * param.coefficient
    end
    return
end

function update_duals_in_affine_objective!(
    param_dual_cum_sum::Vector{Float64},
    affine_objective_cache::Vector{MOI.ScalarAffineTerm{T}},
) where {T}
    for param in affine_objective_cache
        param_dual_cum_sum[param.variable_index.value-PARAMETER_INDEX_THRESHOLD] +=
            param.coefficient
    end
    return
end

function update_duals_in_quadratic_objective!(
    param_dual_cum_sum::Vector{Float64},
    quadratic_objective_cache_pc::Vector{MOI.ScalarAffineTerm{T}},
) where {T}
    for param in quadratic_objective_cache_pc
        param_dual_cum_sum[param.variable_index.value-PARAMETER_INDEX_THRESHOLD] +=
            param.coefficient
    end
    return
end

function empty_and_feed_duals_to_model(
    model::ParametricOptimizer,
    param_dual_cum_sum::Vector{Float64},
)
    empty!(model.dual_value_of_parameters)
    model.dual_value_of_parameters = zeros(model.number_of_parameters_in_model)
    for (vi_val_minus_threshold, param_dual) in enumerate(param_dual_cum_sum)
        model.dual_value_of_parameters[vi_val_minus_threshold] = param_dual
    end
    return
end

struct ParameterDual <: MOI.AbstractVariableAttribute end

function MOI.get(model::ParametricOptimizer, ::ParameterDual, vi_val::Int64)
    if !is_additive(
        model,
        MOI.ConstraintIndex{MOI.SingleVariable,Parameter}(vi_val),
    )
        error("Cannot calculate the dual of a multiplicative parameter")
    end
    return model.dual_value_of_parameters[vi_val-PARAMETER_INDEX_THRESHOLD]
end

function MOI.get(
    model::ParametricOptimizer,
    ::MOI.ConstraintDual,
    cp::MOI.ConstraintIndex{MOI.SingleVariable,POI.Parameter},
)
    if !is_additive(model, cp)
        error("Cannot calculate the dual of a multiplicative parameter")
    end
    return model.dual_value_of_parameters[cp.value-PARAMETER_INDEX_THRESHOLD]
end

function MOI.get(
    model::ParametricOptimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex,
)
    return MOI.get(model.optimizer, attr, ci)
end

function is_additive(model::ParametricOptimizer, cp::MOI.ConstraintIndex)
    if cp.value in model.multiplicative_parameters
        return false
    end
    return true
end


text:
# p13.jl - solve linear BVP u_xx = exp(4x), u(-1)=u(1)=0

N = 16;
(D,x) = cheb(N);
D2 = D^2;
D2 = D2[2:N,2:N];                   # boundary conditions
f = @. exp(4*x[2:N]);
u = D2\f;                           # Poisson eq. solved here
u = [0;u;0];
clf(); axes([.1,.4,.8,.5]);
plot(x,u,".",markersize=6);
xx = -1:.01:1;
uu = polyval(polyfit(x,u),xx);      # interpolate grid data
plot(xx,uu);
grid(true);
exact = @. ( exp(4*xx) - sinh(4)*xx - cosh(4) )/16;
title("max err = $(signif(norm(uu-exact,Inf),4))",fontsize=12);


text:
# FeedForward Networks used after self-attention
"""
    FeedForwardNetwork
"""
struct FFN{L,F,D}
    ffn_layer1::L
    ffn_layer2::L
    activation::F
    dropout::D
end

function FFN(in_dim::Integer,
            hidden_dim::Integer;
            dropout=0.0,
            activation=Flux.relu,
            bias=true,
            init=glorot_uniform)
    ffn_layer1 = Flux.Dense(in_dim, hidden_dim, bias=bias, init=init)
    ffn_layer2 = Flux.Dense(hidden_dim, in_dim, bias=bias, init=init)
    dropout_layer = Flux.Dropout(dropout)
    return FFN(ffn_layer1, ffn_layer2, activation, dropout_layer)
end

@functor FFN

function (ffn::FFN)(x)
    output = ffn.activation.(ffn.ffn_layer1(x))
    output = ffn.dropout(output)
    return ffn.ffn_layer2(output)
end

"""
    GraphTransformerLayer
"""
struct GraphTransformerLayer{M,F,N}
    attention_layer::M
    node_ffn_layer::F
    edge_ffn_layer::F
    node_normalization1::N
    edge_normalization1::N
    node_normalization2::N
    edge_normalization2::N
end

function GraphTransformerLayer(hidden_dim::Integer,
                                num_heads::Integer;
                                dropout=0.0,
                                bias=true,
                                init=glorot_uniform)
    attention_layer = MultiHeadAttentionLayer(hidden_dim, hidden_dim ÷ num_heads, num_heads, bias=bias, init=init)

    # FFN layer for nodes and edges
    node_ffn_layer = FFN(hidden_dim, hidden_dim * 2, dropout=dropout, bias=bias, init=init)
    edge_ffn_layer = FFN(hidden_dim, hidden_dim * 2, dropout=dropout, bias=bias, init=init)

    # Normalization layers
    node_normalization1 = Flux.LayerNorm(hidden_dim)
    edge_normalization1 = Flux.LayerNorm(hidden_dim)
    node_normalization2 = Flux.LayerNorm(hidden_dim)
    edge_normalization2 = Flux.LayerNorm(hidden_dim)

    return GraphTransformerLayer(
                attention_layer,
                node_ffn_layer, edge_ffn_layer,
                node_normalization1, edge_normalization1,
                node_normalization2, edge_normalization2)
end

@functor GraphTransformerLayer

function (GTL::GraphTransformerLayer)(graph::FeaturedGraph)
    init_node_features, init_edge_features = graph.nf, graph.ef
    # Propagate Multi-Head Attention
    node_features, edge_features = GTL.attention_layer(graph)

    # Skip connections
    node_features += init_node_features
    edge_features += init_edge_features

    # Normalization
    normalized_node_features = GTL.node_normalization1(node_features)
    normalized_edge_features = GTL.edge_normalization1(edge_features)

    # FFN transformation
    node_features = GTL.node_ffn_layer(normalized_node_features)
    edge_features = GTL.edge_ffn_layer(normalized_edge_features)

    # Final skip connections
    node_features += normalized_node_features
    edge_features += normalized_edge_features

    # Final normalization
    node_features = GTL.node_normalization2(node_features)
    edge_features = GTL.edge_normalization2(edge_features)

    # Final graph
    return FeaturedGraph(graph, nf=node_features, ef=edge_features)
end

max_stars_repo_name: DwaraknathT/Flux-Graph-Transformer

max_issues_repo_name: DwaraknathT/Flux-Graph-Transformer

max_forks_repo_name: DwaraknathT/Flux-Graph-Transformer


text:
using Base.Threads

"""
    euler_matrix(α::Real, β::Real, γ::Real)

Returns the rotation matrix A based on rotation along the euler angles α, β and γ corresponding to rotations around the x, y, and z-axis respectively.
"""
function euler_matrix(α::Real, β::Real, γ::Real)
    [
        (cos(β)*cos(γ))                         (-cos(β)*sin(γ))                        sin(β)
        (cos(α)*sin(γ) + cos(γ)*sin(α)*sin(β))  (cos(α)*cos(γ) - sin(α)*sin(β)*sin(γ))  (-cos(β)*sin(α))
        (sin(α)*sin(γ) - cos(α)*cos(γ)*sin(β))  (cos(γ)*sin(β) + cos(α)*sin(β)*sin(γ))  (cos(α)*cos(β))

    ]
end

"""
    rotate_3D_quantity(x, α, β, γ)

Rotates a 3D vector along the x-axis with angle α, y-axis with angle β and z-axis with angle γ.
"""
function rotate_3D_quantity(x::Array{<:Real}, α::Real, β::Real, γ::Real)

    A = euler_matrix(α, β, γ)

    return A * x
end

"""
    rotate_3D(x::Array{<:Real}, alpha::Real, beta::Real, gamma::Real)

Rotates and array of 3D positions around the euler angles α, β and γ corresponding to rotations around the x, y, and z-axis respectively.
α, β and γ need to be given in degrees.
"""
function rotate_3D(x::Array{<:Real}, alpha::Real, beta::Real, gamma::Real)

    α = deg2rad(alpha)
    β = deg2rad(beta)
    γ = deg2rad(gamma)

    N = size(x,2)
    ret = Array{eltype(x[1]),2}(undef, 3, N)
    @threads for i = 1:N
        @inbounds ret[:,i] = rotate_3D_quantity(x[:,i], α, β, γ)
    end

    return ret
end

"""
    rotate_3D!(x::Array{<:Real}, alpha::Real, beta::Real, gamma::Real)

Rotates and array of 3D positions around the euler angles α, β and γ corresponding to rotations around the x, y, and z-axis respectively.
α, β and γ need to be given in degrees.
"""
function rotate_3D!(x::Array{<:Real}, alpha::Real, beta::Real, gamma::Real)

    α = deg2rad(alpha)
    β = deg2rad(beta)
    γ = deg2rad(gamma)

    @threads for i = 1:size(x,2)
        @inbounds x[:,i] = rotate_3D_quantity(x[:,i], α, β, γ)
    end

    return x
end

"""
    rotate_to_xz_plane!(x::Array{<:Real})

Rotates an array of 3D positions into the xz-plane.
"""
function rotate_to_xz_plane!(x::Array{<:Real})

    @inbounds for i = 1:size(x,2)
        pos3   = copy(x[2,i])
        x[2,i] = x[3,i]
        x[3,i] = pos3
    end
    x
end

"""
    rotate_to_xz_plane!(x::Array{<:Real})

Rotates an array of 3D positions into the xz-plane.
"""
function rotate_to_xz_plane!(x::Array{<:Real}, x_in::Array{<:Real})

    @inbounds for i = 1:size(x,2)
        x[1,i] = x_in[1,i]
        x[2,i] = x_in[3,i]
        x[3,i] = x_in[2,i]
    end
    x
end

"""
    rotate_to_yz_plane(x::Array{<:Real})

Rotates an array of 3D positions into the yz-plane.
"""
function rotate_to_yz_plane!(x::Array{<:Real})

    @inbounds for i = 1:size(x,2)
        pos3   = copy(x[1,i])
        x[1,i] = x[2,i]
        x[2,i] = x[3,i]
        x[3,i] = pos3
    end
    x
end

"""
    rotate_to_yz_plane(x::Array{<:Real}, x_in::Array{<:Real})

Rotates an array of 3D positions into the yz-plane.
"""
function rotate_to_yz_plane!(x::Array{<:Real}, x_in::Array{<:Real})

    @inbounds for i = 1:size(x,1)
        x[1,i] = x_in[2,i]
        x[2,i] = x_in[3,i]
        x[3,i] = x_in[1,i]
    end
    x
end

"""
    project_along_axis!(x::Array{<:Real}, projection_axis::Integer=3)

Projects and array of 3D along one of the principle axes.
projection_axis ∈ {1, 2, 3} => x, y, z axis.
"""
function project_along_axis!(x::Array{<:Real}, projection_axis::Integer=3)

    # rotation to xy-plane -> nothing is done
    if projection_axis == 3
        return x
    end

    # rotation to xz-plane
    if projection_axis == 2
        return rotate_to_xz_plane!(x)
    end

    # rotation to yz-plane
    if projection_axis == 1
        return rotate_to_yz_plane!(x)
    end
end

"""
    project_along_axis!(x::Array{<:Real}, projection_axis::Integer=3)

Projects and array of 3D along one of the principle axes.
projection_axis ∈ {1, 2, 3} => x, y, z axis.
"""
function project_along_axis!(x::Array{<:Real}, x_in::Array{<:Real}, projection_axis::Integer=3)

    # rotation to xy-plane -> nothing is done
    if projection_axis == 3
        return x_in
    end

    # rotation to xz-plane
    if projection_axis == 2
        return rotate_to_xz_plane!(x, x_in)
    end

    # rotation to yz-plane
    if projection_axis == 1
        return rotate_to_yz_plane!(x, x_in)
    end
end

"""
    project_along_axis!(x::Array{<:Real}, projection_axis::Integer=3)

Projects and array of 3D along one of the principle axes.
projection_axis ∈ {1, 2, 3} => x, y, z axis.
"""
function project_along_axis(x::Array{<:Real}, projection_axis::Integer=3)

    # allocate new array
    ret = Array{eltype(x[1]),2}(undef, 3, size(x,2))

    # rotation to xy-plane -> nothing is done
    if projection_axis == 3
        return x
    end

    # rotation to xz-plane
    if projection_axis == 2
        return rotate_to_xz_plane!(ret, x)
    end

    # rotation to yz-plane
    if projection_axis == 1
        return rotate_to_yz_plane!(ret, x)
    end
end


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(SplitOctonion(a))
end

@test begin
    a = 1
    !isreal(SplitOctonion(a, a, a, a))
end

@test begin
    z = random(SplitOctonion{Int})
    z == +(z)
end

@test begin
    a = CayleyDickson.randomBigFloat()
    real(SplitOctonion(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, SplitOctonion(1.0, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, s: 5.0, is: 6.0, js: 7.0, (ij)s: 8.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, SplitOctonion(1.0, 2, 3, 4, 5, 6, 7))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, s: 5.0, is: 6.0, js: 7.0, (ij)s: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, SplitOctonion(1.0, 2, 3, 4, 5, 6))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, s: 5.0, is: 6.0, js: 0.0, (ij)s: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, SplitOctonion(1.0, 2, 3, 4, 5))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, s: 5.0, is: 0.0, js: 0.0, (ij)s: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, SplitOctonion(1.0, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, j: 3.0, ij: 4.0, s: 0.0, is: 0.0, js: 0.0, (ij)s: 0.0]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    d = CayleyDickson.randomBigFloat()
    l = SplitOctonion(Quaternion(a, b, c, d))
    r = SplitOctonion(a, b, c, d)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = SplitOctonion(Binion(a, b))
    r = SplitOctonion(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = SplitOctonion(a, b, c)
    r = SplitOctonion(a, b, c, 0)
    l == r
end

@test begin
    x = random(SplitOctonion{BigFloat})
    y = random(SplitOctonion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(SplitOctonion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Quaternion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(SplitOctonion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Quaternion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l != r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = alternatorL(x, y)
    r = zero(SplitOctonion{BigInt})
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = alternatorR(x, y)
    r = zero(SplitOctonion{BigInt})
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(SplitOctonion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(SplitOctonion{Rational{BigInt}})
    y = random(SplitOctonion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(SplitOctonion{Rational{BigInt}})
    y = random(SplitOctonion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(SplitOctonion{Rational{BigInt}})
    y = random(SplitOctonion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(SplitOctonion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Quaternion{BigInt})
    x = random(SplitOctonion{BigInt})
    l = commutator(a, x)
    r = zero(SplitOctonion{BigInt})
    l != r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    z = random(SplitOctonion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(SplitOctonion{BigInt})
    quadrance(x) > big(0) || quadrance(x) < big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(SplitOctonion{BigInt})
    y = random(SplitOctonion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(SplitOctonion(0, 0.0))
end

@test_throws ErrorException begin
    random(SplitOctonion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(SplitOctonion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(SplitOctonion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(SplitOctonion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(SplitOctonion{BigInt})
    z == z + zero(SplitOctonion{BigInt})
end

@test begin
    z = random(SplitOctonion{BigInt})
    z == z * one(SplitOctonion{BigInt})
end

@test begin
    z = random(SplitOctonion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(SplitOctonion{BigInt})
    z == z * one(z)
end


text:
export For
using Random
import Base

"""
    For(f, base...)

`For` provides a convenient way to construct a `ProductMeasure`. There are
several options for the `base`. With Julia's `do` notation, this can look very
similar to a standard `for` loop, while maintaining semantics structure that's
easier to work with.

------------

# `For(f, base::Int...)`

When one or several `Int` values are passed for `base`, the result is treated as
depending on `CartesianIndices(base)`.

```
julia> For(3) do λ Exponential(λ) end |> marginals
3-element mappedarray(MeasureTheory.var"#17#18"{var"#15#16"}(var"#15#16"()), ::CartesianIndices{1, Tuple{Base.OneTo{Int64}}}) with eltype Exponential{(:λ,), Tuple{Int64}}:
 Exponential(λ = 1,)
 Exponential(λ = 2,)
 Exponential(λ = 3,)
```

```
julia> For(4,3) do μ,σ Normal(μ,σ) end |> marginals
4×3 mappedarray(MeasureTheory.var"#17#18"{var"#11#12"}(var"#11#12"()), ::CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}) with eltype Normal{(:μ, :σ), Tuple{Int64, Int64}}:
 Normal(μ = 1, σ = 1)  Normal(μ = 1, σ = 2)  Normal(μ = 1, σ = 3)
 Normal(μ = 2, σ = 1)  Normal(μ = 2, σ = 2)  Normal(μ = 2, σ = 3)
 Normal(μ = 3, σ = 1)  Normal(μ = 3, σ = 2)  Normal(μ = 3, σ = 3)
 Normal(μ = 4, σ = 1)  Normal(μ = 4, σ = 2)  Normal(μ = 4, σ = 3)
```

-------

# `For(f, base::AbstractArray...)``

In this case, `base` behaves as if the arrays are `zip`ped together before
applying the map.

```
julia> For(randn(3)) do x Exponential(x) end |> marginals
3-element mappedarray(x->Main.Exponential(x), ::Vector{Float64}) with eltype Exponential{(:λ,), Tuple{Float64}}:
 Exponential(λ = -0.268256,)
 Exponential(λ = 1.53044,)
 Exponential(λ = -1.08839,)
```

```
julia> For(1:3, 1:3) do μ,σ Normal(μ,σ) end |> marginals
3-element mappedarray((:μ, :σ)->Main.Normal(μ, σ), ::UnitRange{Int64}, ::UnitRange{Int64}) with eltype Normal{(:μ, :σ), Tuple{Int64, Int64}}:
 Normal(μ = 1, σ = 1)
 Normal(μ = 2, σ = 2)
 Normal(μ = 3, σ = 3)
```

----

# `For(f, base::Base.Generator)`

For `Generator`s, the function maps over the values of the generator:

```
julia> For(eachrow(rand(4,2))) do x Normal(x[1], x[2]) end |> marginals |> collect
4-element Vector{Normal{(:μ, :σ), Tuple{Float64, Float64}}}:
 Normal(μ = 0.255024, σ = 0.570142)
 Normal(μ = 0.970706, σ = 0.0776745)
 Normal(μ = 0.731491, σ = 0.505837)
 Normal(μ = 0.563112, σ = 0.98307)
```

"""
For(f, dims...) = ProductMeasure(i -> f(i...), zip(dims...))

For(f, inds::AbstractArray) = ProductMeasure(f, inds)

For(f, n::Int) = ProductMeasure(f, 1:n)
For(f, dims::Int...) = ProductMeasure(i -> f(Tuple(i)...), CartesianIndices(dims))

function Base.eltype(d::ProductMeasure{F,I}) where {F,I<:AbstractArray}
    return eltype(d.f(first(d.pars)))
end

# """
#     indexstyle(a::AbstractArray, b::AbstractArray)

# Find the best IndexStyle that works for both `a` and `b`. This will return
# `IndexLinear` if both `a` and `b` support it; otherwise it will fall back on `IndexCartesian`.
# """
# function indexstyle(::A,::B)
#     if IndexStyle(A) == IndexStyle(B) == IndexLinear()
#         return IndexLinear()
#     end

#     return IndexCartesian()
# end

# function Base.rand(rng::AbstractRNG, μ::ForArray{D,N,T,F}) where {F,T<:AbstractArray,D,X}
#     s = size(μ.θ)
#     x = Array{X,length(s)}(undef, s...)
#     rand!(rng, x, μ)
# end

# function logdensity(μ::ForArray{D,N,T,F}, x)
#     getℓ(θⱼ, xⱼ) = logdensity(μ.f(θⱼ), xⱼ)
#     ℓ = mappedarray(getℓ, μ.θ, x)
#     _logdensity(μ, x, indexstyle(μ.θ, x), result_type)
# end

# function _logdensity(μ::ForArray{D,N,T,F}, x, ::IndexLinear, ::Type{R}) where {R<:AbstractFloat}
#     ℓ = zero(R)
#     μ.f(μ.θ)
# end

# function basemeasure(μ::ForArray{D,N,T,F}) where {F,T<:AbstractArray,D,X}

# ForGenerator


text:
module AEEIncidentes

using HTTP: request
using JSON3: JSON3, StructType, Struct, Mutable
using LibPQ: Connection, load!, execute
using Dates: DateTime, now, DateFormat, Hour
using TimeZones: ZonedDateTime, utc_tz, UTC
import Base: getproperty
import DataFrames: DataFrame

const AEEDTF = DateFormat("yyyy-mm-ddTHH:MM:SSz")

conn = Connection("dbname = sdad")
latest = execute(conn, "SELECT MAX(statustime) latest FROM aeepr.incidentes;", not_null = true) |>
    (obj -> round(only(getproperty.(obj, :latest)), Hour))
mutable struct Artifact
    id :: UInt64
    created_at :: String
    Artifact() = new()
end
StructType(::Type{Artifact}) = Mutable()
struct Artifacts
    total_count :: UInt16
    artifacts :: Vector{Artifact}
end
StructType(::Type{Artifacts}) = Struct()
function getproperty(obj::Artifact, sym::Symbol)
    if sym == :created_at
        created_at = floor(DateTime(ZonedDateTime(getfield(obj, sym), AEEDTF), UTC), Hour)
    else
        getfield(obj, sym)
    end
end
function list_artifacts(latest::DateTime = now() - Day(92))
    response = request("GET",
                       "https://api.github.com/repos/Nosferican/AEEIncidentes/actions/artifacts?per_page=100",
                       ["Accept" => "application/vnd.github.v3+json",
                        "Time-Zone" => "GMT",
                       ])
    @assert response.status == 200 "GH API is not working"
    json = JSON3.read(response.body, Artifacts)
    artifacts = json.artifacts
    @assert latest ≥ artifacts[end].created_at "Need manual check: older releases need to be added"
    filter(artifact -> artifact.created_at > latest, artifacts)
end
function DataFrame(obj::Artifact)
    # obj = artifacts[51]
    statustime = floor(obj.created_at, Hour)
    filename = string(string(statustime)[1:13], ".jsonl")
    if !isfile(joinpath("data", filename))
        response = request("GET",
                           "https://api.github.com/repos/Nosferican/AEEIncidentes/actions/artifacts/$(obj.id)/zip",
                           ["Accept" => "application/vnd.github.v3+json",
                            "Time-Zone" => "GMT",
                            "User-Agent" => ENV["GH_USR"],
                            "Authorization" => "token $(ENV["GH_PAT_AEE"])",
                           ])
        write(joinpath("data", "$filename.zip"), response.body)
        run(`unzip -q data/$filename.zip -d data`)
        rm(joinpath("data", "$filename.zip"))
    end
    lns = filter!(!isempty, readlines(joinpath("data", filename)))
    if !isempty(lns)
        data = DataFrame(reduce(vcat, JSON3.read(ln) for ln in lns))
        data[!,:statustime] .= statustime
    else
        data = DataFrame([DateTime, String, String], [:statustime, :zone, :area], 0)
    end
    data
end
artifacts = list_artifacts(latest)
if !isempty(artifacts)
    data = sort!(reduce(vcat, DataFrame(artifact) for artifact in artifacts)[!,[:statustime, :area, :zone]])
    execute(conn, "BEGIN;")
    load!(data,
          conn,
          string("INSERT INTO aeepr.incidentes VALUES(",
                 join(("\$$i" for i in 1:size(data, 2)), ','),
                 ") ON CONFLICT ON CONSTRAINT incidentes_pkey DO NOTHING;"))
    execute(conn, "COMMIT;")
    close(conn)
    foreach(file -> rm(joinpath("data", file)), filter!(x -> occursin(r"\d{4}-\d{2}-\d{2}T\d{2}\.json", x), readdir("data")))
end

end # module

max_stars_repo_name: Nosferican/AEEIncidentes

max_issues_repo_name: Nosferican/AEEIncidentes

max_forks_repo_name: Nosferican/AEEIncidentes


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libtiff"
version = v"4.0.9"

# Collection of sources required to build libtiff
sources = [
    "https://gitlab.com/libtiff/libtiff/-/archive/Release-v4-0-9/libtiff-Release-v4-0-9.tar.gz" =>
    "6656637a1183c32227334fc0d0d08e011e3da74804bcdcaec496723a9612f97a",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libtiff*/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtiff", :libtiff),
    LibraryProduct(prefix, "libtiffxx", :libtiffxx)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/JuliaIO/LibpngBuilder/releases/download/v1.0.3/build_libpng.v1.6.37.jl",
    "https://github.com/SimonDanisch/LibJPEGBuilder/releases/download/v10/build_libjpeg.v9.0.0-b.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
include("test_exports.jl")
include("test_div_by_zero.jl")
include("test_basic_dh.jl")
include("test_geometry_kernel.jl")
include("test_poly_approx.jl")
include("test_utility.jl")
include("test_vector_projections.jl")
include("test_rotations.jl")


text:
"""
    getlast(sol, vars...)
    getlast(sol, vars_vec)

Apply getindex for each var to sol and return a NamedArray.

```jldocstest; output=false
using ModelingToolkit, DifferentialEquations
using MTKHelpers
using NamedArrays
function samplesystem(;name)
    @variables t
    D = Differential(t)
    sts = @variables x(t) RHS(t)  # RHS is observed
    ps = @parameters τ       # parameters
    ODESystem([ RHS  ~ (1 - x)/τ, D(x) ~ RHS ], t, sts, ps; name)
end
@named m = samplesystem()
@named sys = embed_system(m)
prob = ODEProblem(sys, [m.x => 0.0], (0.0,10.0), [m.τ => 3.0])
sol = solve(prob);
res = getlast(sol, m.x, m.RHS)
res == NamedArray([sol[m.x,end], sol[m.RHS,end]], ([m.x, m.RHS],))
# output
true
```
"""
function getlast(sol::SciMLBase.AbstractODESolution, vars...; kwargs...)
    getlast(sol, collect(vars))
end
function getlast(sol::SciMLBase.AbstractODESolution, vars::AbstractVector)
    iend = lastindex(sol,2)
    names_vars_vec = _names_vars(vars)
    #names_vars_vec = names_vars #convert(Array, names_vars)::Vector{eltype(names_vars)}
    a = getindex.(Ref(sol), names_vars_vec, iend)::Array{eltype(sol)}
    #names = ntuple(i -> symbol(vars[i]), length(vars))
    #A = @LArray a names
    NamedArray(a, (names_vars_vec,))
end
_names_vars(vars_vec) = vars_vec
_names_vars(vars_vec::NamedArray) = first(names(vars_vec))


text:
using Base: Float64

#* Flow functions
abstract type RelPerms end

mutable struct SingleRelPerms <: RelPerms
    swr::Float64
    sor::Float64
    krw0::Float64
    kro0::Float64
    nw::Float64
    no::Float64
end

mutable struct MixedRelPerms <: RelPerms
    ww::SingleRelPerms
    ow::SingleRelPerms
end

function SingleRelPerms(;swr::T, sor::T, krw0::T, kro0::T, nw::T, no::T) where T <: Float64
    SingleRelPerms(swr, sor, krw0, kro0, nw, no)
end

#* Reservoir
abstract type Reservoir end

struct Cylinder <: Reservoir
    L   ::Float64
    D   ::Float64
    A   ::Float64
    V   ::Float64
    k   ::Float64
    phi ::Float64
end

struct Block <: Reservoir
    L   ::Float64
    W   ::Float64
    h   ::Float64
    A   ::Float64
    V   ::Float64
    k   ::Float64
    phi ::Float64
end

function Reservoir(L::Float64, D::Float64, k::Float64, phi::Float64)
    A = π * D ^2 /4
    V = A * L
    Cylinder(L, D, A, V, k, phi)
end

function Reservoir(L::Float64, W::Float64, h::Float64, k::Float64, phi::Float64)
    V = L * W * h
    A = W * h
    Block(L, W, h, A, V, k, phi)
end

#* Fluid
struct Fluids{T <: Real}
    μw::T
    μo::T
end
function Fluids(;μw::T, μo::T) where T <: Float64
    Fluids{Float64}(μw, μo)
end

#* Grid
abstract type Grid end

mutable struct Grid1D <: Grid
    Nx  ::Int64
    N   ::Int64
    dx  ::Float64
    x   ::Vector{Float64}
end

mutable struct Grid2D <: Grid
    Nx  ::Int64
    Ny  ::Int64
    N   ::Int64
    dx  ::Float64
    dy  ::Float64
    x   ::Vector{Float64}
    y   ::Vector{Float64}
end

function Grid(res::Reservoir, Nx::Int64)
    dx = res.L / Nx
    x = collect(dx/2:dx:res.L-dx/2)

    Grid1D(Nx, Nx, dx, x)
end

function Grid(res::Reservoir, Nx::Int64, Ny::Int64)
    dx = res.L / Nx
    dy = res.W / Ny
    x = collect(dx/2:dx:res.L-dx/2)
    y = collect(dy/2:dy:res.L-dy/2)
    N = Nx * Ny

    Grid2D(Nx, Ny, N, dx, dy, x, y)
end

# function Grid(;res::Reservoir, Nx::T, Ny) where T <: Int64
#     dx = res.L / Nx
#     x = collect(dx/2:dx:res.L-dx/2)

#     Grid(Nx, dx, x)
# end

struct BoundaryConditions
    type:: Vector{Symbol}
    value:: Vector{Real}
end

abstract type Simulation end

mutable struct Simulation1D <: Simulation
    p   ::Vector{Float64}
    s   ::Vector{Float64}
    θ   ::Vector{Float64}
    v   ::Vector{Float64}
    i   ::Int64
    t   ::Float64
    tmax::Float64
    dt  ::Float64
end

mutable struct Simulation2D <: Simulation
    p   ::Matrix{Float64}
    s   ::Matrix{Float64}
    θ   ::Matrix{Float64}
    v   ::Matrix{Float64}
    i   ::Int64
    t   ::Float64
    tmax::Float64
    dt  ::Float64
end

function Simulation(p::T, s::T, θ::T, tmax::Float64, dt::Float64) where T <: Vector{Float64}
    n = length(p)
    Simulation1D(p, s, θ, zeros(n+1), 1, 0., tmax, dt)
end

function Simulation(p::T, s::T, θ::T, tmax::Float64, dt::Float64) where T <: Matrix{Float64}
    Nx, Ny = size(p)
    Simulation2D(p, s, θ, zeros(Nx, Ny), 1, 0., tmax, dt)
end

mutable struct FlowResults
    t::Vector{Float64}
    p::Matrix{Float64}
    s::Matrix{Float64}
    θ::Matrix{Float64}

   function FlowResults(sim::Simulation, grid::Grid, tmax::Float64, dt::Float64)
    ntimes =  Int(tmax ÷ dt)
    p = zeros(grid.Nx, ntimes+1)
    s = zeros(grid.Nx, ntimes+1)
    θ = zeros(grid.Nx, ntimes+1)

    t = collect(0:dt:ntimes*dt)
    p[:,1] = sim.p
    s[:,1] = sim.s
    θ[:,1] = sim.θ

    new(t, p, s, θ)
    end
end

mutable struct CompTransport{N}
    c:: Array{Float64, N}
    ncomps::Int64
    cinj::Vector{Float64}
end


text:
#####
##### Formulas assuming dry (unsaturated) adiabatic (i.e. isentropic) processes
#####

export DryAdiabaticProcess

export air_pressure_given_θ, air_pressure, air_temperature

"""
    DryAdiabaticProcess

For dispatching to isentropic formulas
"""
struct DryAdiabaticProcess end

"""
    air_pressure_given_θ(θ::FT, Φ::FT, ::DryAdiabaticProcess)

The air pressure for an isentropic process, where

 - `θ` potential temperature
 - `Φ` gravitational potential
"""
air_pressure_given_θ(θ::FT, Φ::FT, ::DryAdiabaticProcess) where {FT} =
  FT(MSLP) * (1 - Φ / (θ * FT(cp_d))) ^ (FT(cp_d) / FT(R_d))

"""
    air_pressure(T::FT, T∞::FT, p∞::FT, ::DryAdiabaticProcess)

The air pressure for an isentropic process, where

 - `T` temperature
 - `T∞` ambient temperature
 - `p∞` ambient pressure
"""
air_pressure(T::FT, T∞::FT, p∞::FT, ::DryAdiabaticProcess) where {FT} =
  p∞ * (T / T∞) ^ (FT(1) / FT(kappa_d))

"""
    air_temperature(p::FT, θ::FT, Φ::FT, ::DryAdiabaticProcess)

The air temperature for an isentropic process, where

 - `p` pressure
 - `θ` potential temperature
"""
air_temperature(p::FT, θ::FT, ::DryAdiabaticProcess) where {FT} =
  (p / FT(MSLP)) ^ (FT(R_d) / FT(cp_d)) * θ

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 1000,)
insts = OrderedDict()
insts["minimal"] = [
    ((:matpoly2, true, true),),
    ((:matpoly5, true, true),),
    ((:matpoly5, true, false),),
    ((:matpoly5, false, false), :SOCExpPSD),
    ]
insts["fast"] = [
    ((:matpoly1, true, true),),
    ((:matpoly1, true, false),),
    ((:matpoly1, false, false), :SOCExpPSD),
    ((:matpoly2, true, true),),
    ((:matpoly2, true, false),),
    ((:matpoly2, false, false), :SOCExpPSD, relaxed_tols),
    ((:matpoly3, true, true),),
    ((:matpoly3, true, false),),
    ((:matpoly3, false, false), :SOCExpPSD),
    ((:matpoly4, true, true),),
    ((:matpoly4, true, false),),
    ((:matpoly4, false, false), :SOCExpPSD),
    ((:matpoly6, true, true),),
    ((:matpoly6, true, false),),
    ((:matpoly6, false, false), :SOCExpPSD),
    ((:matpoly7, true, true), nothing, relaxed_tols),
    ((:matpoly7, true, false),),
    ((:matpoly7, false, false), :SOCExpPSD),
    ]
insts["various"] = insts["fast"]
return (SemidefinitePolyJuMP, insts)


text:
struct RelationType <: AbstractRelationType
    label::Label
    is_root::Bool
end

# Remote functions

function create(x::RemoteConcept{RelationType})
    result = execute(x.transaction, RelationTypeRequestBuilder.create_req(x.concept.label))
    return Relation(result.type_res.relation_type_create_res.relation)
end

function set_relates(transaction::AbstractCoreTransaction,
                    label::Label,
                    role_label::AbstractString,
                    overridden_label::Optional{<:AbstractString} = nothing)

    relates_req = RelationTypeRequestBuilder.set_relates_req(label,
                                    role_label,
                                    overridden_label)

    execute(transaction, relates_req)
    return nothing
end

function get_relates(r::RemoteConcept{RelationType})
    relates_req = RelationTypeRequestBuilder.get_relates_req(r.concept.label)
    res = stream(r.transaction, relates_req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.relation_type_get_relates_res_part.roles for r in res)))
end

function relation_type_get_relates_for_role_label(r::RemoteConcept{RelationType}, for_role::AbstractString)
    role_req = RelationTypeRequestBuilder.relation_type_get_relates_for_role_label_req(r.concept.label,for_role)
    res = execute(r.transaction, role_req)

    if !hasproperty(res.type_res.relation_type_get_relates_for_role_label_res, :role_type)
        return nothing
    end
    return instantiate(res.type_res.relation_type_get_relates_for_role_label_res.role_type)
end

function unset_relates(r::RemoteConcept{RelationType}, role_label::Optional{String})
    unset_req = RelationTypeRequestBuilder.unset_relates_req(r.concept.label, role_label)
    execute(r.transaction, unset_req)
    return nothing
end


text:
module FastRationals2017

export FastRational

import Base: convert, promote_rule, eltype,
    string, show, read, write,
    checked_add, checked_sub, checked_mul, power_by_squaring,
    numerator, denominator, widen, rationalize,
    isinteger, iszero, #isone,
    sign, signbit, copysign, flipsign, typemin, typemax, abs,
    ==, !=, <=, <, >=, >, isless, isequal, cmp,
    -, +, *, inv, /, //, rem, mod, fma, div, fld, cld,
    trunc, floor, ceil, round, ^

import Base.Checked: add_with_overflow, sub_with_overflow, mul_with_overflow

const SignedInt = Union{Int64, Int32, Int16, Int8, Int128, BigInt}

if !isdefined(:NamedTuple)
    using NamedTuples
end

"""
    canonical(numerator::T, denominator::T) where T<:Signed
Rational numbers that are finite and normatively bounded by
the extremal magnitude of an underlying signed type have a
canonical representation.
- numerator and denominator have no common factors
- numerator may be negative, zero or positive
- denominator is strictly positive (d > 0)
""" canonical

function canonical(num::T, den::T) where {T<:SignedInt}
    num, den = canonical_signed(num, den)
    num, den = canonical_valued(num, den)
    return num, den
end

@inline function canonical_signed(num::T, den::T) where {T<:SignedInt}
    return flipsign(num, den), abs(den)
end

@inline function canonical_valued(num::T, den::T) where {T<:SignedInt}
    gcdval = gcd(num, den)
    gcdval === one(T) && return num, den
    num = div(num, gcdval)
    den = div(den, gcdval)
    return num, den
end

abstract type AbstractRational{Q}  <: Real end  # Q is underlying type (Int64 ?Poly)
abstract type SignedRational{Q}    <: AbstractRational{Q} end  # may be reduceable
abstract type CanonicalRational{Q} <: SignedRational{Q}   end  # already reduced

@inline numerator(q::T)   where {Q, T<:AbstractRational{Q}} = q.num
@inline denominator(q::T) where {Q, T<:AbstractRational{Q}} = q.den
@inline value(q::T) where {Q, T<:AbstractRational{Q}} =
    numerator(a), denominator(q)
@inline value(q::T) where {Q, T<:Rational{Q}} =
    numerator(a), denominator(q)

mutable struct SignedRatio{Q} <: SignedRational{Q}
    num::Q    # numerator
    den::Q    # denominator > 0
    can::Bool # in canonical form?
end

struct RationalNum{Q} <: CanonicalRational{Q}
    num::Q
    den::Q
end

# tuples pass the given values to the constructor
@inline SignedRatio(numden::Tuple{I,I}) where I<:SignedInt =
    SignedRatio(numden[1], numden[2], false)
SignedRatio(numden::Tuple{I,I,Bool}) where {I<:SignedInt} =
    SignedRatio(numden...,)
SignedRatio(numden::Tuple{I,I}, can::Bool) where {I<:SignedInt} =
    SignedRatio(numden[1], numden[2], can)

# explicitly type qualified dispatch *may* modify given values
@inline SignedRatio(::Type{I}, num::I, den::I) where I<:SignedInt =
    SignedRatio(num, den)

# applicative constructors work properly
@inline SignedRatio(q::RationalNum{I}) where I<:SignedInt =
    SignedRatio(value(q))

# applicative constructors work properly
@inline SignedRatio(q::Rational{I}) where I<:SignedInt =
    SignedRatio(value(q))

# tuples pass the given values to the constructor
@inline RationalNum(numden::Tuple{I,I}) where I<:SignedInt =
    RationalNum(numden[1], numden[2])

# explicitly type qualified dispatch *may* modify given values
@inline RationalNum(::Type{I}, num::I, den::I) where I<:SignedInt =
    RationalNum(canonical(num, den))

# applicative constructors work properly
@inline RationalNum(q::SignedRatio{I}) where I<:SignedInt =
    q.can ? RationalNum(numerator(q), denominator(q)) :
            RationalNum(canonical(numerator(q), denominator(q)))

# applicative constructors work properly
@inline RationalNum(q::Rational{I}) where I<:SignedInt =
    RationalNum(value(q))

# Rational as an applicative constructor

@inline Rational(q::SignedRatio{I}) where I<:SignedInt =
    Rational(RationalNum(q))

@inline Rational(q::RationalNum{I}) where I<:SignedInt =
    Rational(numerator(q), denominator(q))

const PlainRational = NamedTuple{(:num, :den, :void)}

@inline numerator(x::PlainRational) = x.num
@inline denominator(x::PlainRational) = x.den
@inline value(x::PlainRational) = (x.num, x.den)
eltype(x::PlainRational) = typeof(x.num)

@inline canonical(x::PlainRational) = canonical(numerator(x), denominator(x))

@inline PlainRational(num::T, den::T) where T = PlainRational((num, den, nothing))
@inline PlainRational(numden::Tuple{T,T}) where T = PlainRational((numden..., nothing))

const FastRational = NamedTuple{(:num, :den)}

@inline numerator(x::FastRational) = x.num
@inline denominator(x::FastRational) = x.den
@inline value(x::FastRational) = (x.num, x.den)
eltype(x::FastRational) = typeof(x.num)

@inline canonical(x::FastRational) = values(x)

@inline FastRational(num::T, den::T) where T = FastRational(canonical(num, den))
@inline FastRational(q::Rational{T}) where T = FastRational((q.num, q.den))
@inline FastRational(q::PlainRational) = FastRational(canonical(q.num, q.den))

@inline canonize(q::PlainRational) = FastRational(canonical(q))
@inline canonize(q::FastRational) = q

convert(::Type{FastRational}, q::PlainRational) = FastRational( canonical(q) )
promote_rule(::Type{FastRational}, ::Type{PlainRational}) = FastRational
convert(::Type{FastRational}, q::Rational{T}) where T<:Signed = FastRational(numerator(q), denominator(q))
convert(::Type{Rational{T}}, q::FastRational) where T<:Signed = Rational(numerator(q), denominator(q))
promote_rule(::Type{FastRational}, ::Type{Rational{T}}) where T<:Signed = FastRational

PlainRational(q::Rational{T}) where T = throw(ErrorException("PlainRational(Rational) is disallowed."))
PlainRational(q::FastRational) = throw(ErrorException("PlainRational(FastRational) is disallowed."))
convert(::Type{PlainRational}, x::FastRational) = throw(ErrorException("conversion from FastRational to PlainRational is disallowed."))

@inline sign(x::FastRational)  = oftype(x, sign(numerator(x)))
@inline signbit(x::FastRational)  = signbit(numerator(x))
@inline copysign(x::FastRational, y::Real) = FastRational(copysign(numerator(x), y), denominator(x))
@inline copysign(x::FastRational, y::FastRational)  = FastRational(copysign(numerator(x), numerator(y)), denominator(x))
@inline flipsign(x::FastRational, y::Real)  = FastRational(flipsign(numerator(x), y), denominator(x))
@inline flipsign(x::FastRational, y::FastRational)  = FastRational(flipsign(numerator(x), numerator(y)), denominator(x))

@inline isinteger(x::FastRational) = abs(denominator(x)) === one(T)
@inline iszero(x::FastRational)  = abs(numerator(x)) === zero(T)

@inline abs(x::FastRational)   = FastRational((abs(numerator(x)), denominator(x)))

@inline inv(x::FastRational)  = FastRational((denominator(x), numerator(x)))

@inline function (-)(x::FastRational)
    numerator(x) === typemin(T) && throw(OverflowError())
    return FastRational((-numerator(x), denominator(x)))
end

(/)(x::FastRational, y::FastRational) = x // y

# add

function (+)(x::FastRational, y::FastRational)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x + $y overflowed"))

    return PlainRational(numer, denom)
end

function (+)(x::PlainRational, y::PlainRational)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    x = canonize(x)
    y = canonize(y)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x + $y overflowed"))

    return PlainRational(numer, denom)
end

function (+)(x::FastRational, y::PlainRational)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x + $y overflowed"))

    return PlainRational(numer, denom)
end

function (+)(x::PlainRational, y::FastRational)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = add_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x + $y overflowed"))

    return PlainRational(numer, denom)
end

(+)(x::FastRational, y::Rational{T})  where T = (+)(promote(x, y)...)
(+)(x::Rational{T}, y::FastRational)  where T = (+)(promote(x, y)...)
(+)(x::PlainRational, y::Rational{T}) where T = (+)(promote(x, y)...)
(+)(x::Rational{T}, y::PlainRational) where T = (+)(promote(x, y)...)

# subtract

function (-)(x::FastRational, y::FastRational)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x - $y overflowed"))

    return PlainRational(numer, denom)
end

function (-)(x::PlainRational, y::PlainRational)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    x = canonize(x)
    y = canonize(y)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x - $y overflowed"))

    return PlainRational(numer, denom)
end

function (-)(x::FastRational, y::PlainRational)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x - $y overflowed"))

    return PlainRational(numer, denom)
end

function (-)(x::PlainRational, y::FastRational)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = sub_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x - $y overflowed"))

    return PlainRational(numer, denom)
end

(-)(x::FastRational, y::Rational{T})  where T = (-)(promote(x, y)...)
(-)(x::Rational{T}, y::FastRational)  where T = (-)(promote(x, y)...)
(-)(x::PlainRational, y::Rational{T}) where T = (-)(promote(x, y)...)
(-)(x::Rational{T}, y::PlainRational) where T = (-)(promote(x, y)...)

# multiply

function (*)(x::FastRational, y::FastRational)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x * $y overflowed"))

    return PlainRational(numer, denom)
end

function (*)(x::PlainRational, y::PlainRational)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    x = canonize(x)
    y = canonize(y)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x * $y overflowed"))

    return PlainRational(numer, denom)
end

function (*)(x::FastRational, y::PlainRational)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x * $y overflowed"))

    return PlainRational(numer, denom)
end

function (*)(x::PlainRational, y::FastRational)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = mul_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x * $y overflowed"))

    return PlainRational(numer, denom)
end

(*)(x::FastRational, y::Rational{T})  where T = (*)(promote(x, y)...)
(*)(x::Rational{T}, y::FastRational)  where T = (*)(promote(x, y)...)
(*)(x::PlainRational, y::Rational{T}) where T = (*)(promote(x, y)...)
(*)(x::Rational{T}, y::PlainRational) where T = (*)(promote(x, y)...)

# divide

function (//)(x::FastRational, y::FastRational)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x // $y overflowed"))

    return PlainRational(numer, denom)
end

function (//)(x::PlainRational, y::PlainRational)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    x = canonize(x)
    y = canonize(y)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x // $y overflowed"))

    return PlainRational(numer, denom)
end

function (//)(x::FastRational, y::PlainRational)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x // $y overflowed"))

    return PlainRational(numer, denom)
end

function (//)(x::PlainRational, y::FastRational)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    !ovf && return PlainRational(numer, denom)

    y = canonize(y)
    numer, denom, ovf = div_with_overflow_for_rational(x, y)
    ovf && throw(OverflowError("$x // $y overflowed"))

    return PlainRational(numer, denom)
end

(//)(x::FastRational, y::Rational{T})  where T = (//)(promote(x, y)...)
(//)(x::Rational{T}, y::FastRational)  where T = (//)(promote(x, y)...)
(//)(x::PlainRational, y::Rational{T}) where T = (//)(promote(x, y)...)
(//)(x::Rational{T}, y::PlainRational) where T = (//)(promote(x, y)...)

# core parts of add, sub

@inline function add_with_overflow_for_rational(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), denominator(y)) # here, numer is a temp
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), numerator(y)) # here, denom is a temp
    ovf |= ovfl
    numer, ovfl = add_with_overflow(numer, denom) # numerator of sum
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), denominator(y)) # denominator of sum
    ovf |= ovfl

    return numer, denom, ovf
end

@inline function sub_with_overflow_for_rational(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), denominator(y)) # here, numer is a temp
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), numerator(y)) # here, denom is a temp
    ovf |= ovfl
    numer, ovfl = sub_with_overflow(numer, denom) # numerator of difference
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), denominator(y)) # denominator of difference
    ovf |= ovfl

    return numer, denom, ovf
end

@inline function mul_with_overflow_for_rational(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), numerator(y))
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), denominator(y))
    ovf |= ovfl

    return numer, denom, ovf
end

@inline function div_with_overflow_for_rational(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(numerator(x), denominator(y))
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(denominator(x), numerator(y))
    ovf |= ovfl

    return numer, denom, ovf
end

(==)(x::Rational{T}, y::FastRational) where T =
   numerator(x) == numerator(y) && denominator(x) == denominator(y)
(!=)(x::Rational{T}, y::FastRational) where T = !(x == y)
(==)(x::FastRational, y::Rational{T}) where T =
   numerator(x) == numerator(y) && denominator(x) == denominator(y)
(!=)(x::FastRational, y::Rational{T}) where T = !(x == y)

(==)(x::Rational{T}, y::PlainRational) where T = (numerator(x), denominator(x)) == canonical(y)
(!=)(x::Rational{T}, y::PlainRational) where T = (numerator(x), denominator(x)) != canonical(y)
(==)(x::PlainRational, y::Rational{T}) where T = canonical(x) == (numerator(t), denominator(y))
(!=)(x::PlainRational, y::Rational{T}) where T = canonical(x) != (numerator(y), denominator(y))

(<)(x::Rational{T}, y::FastRational) where T = x < Rational{T}(y)
(<=)(x::Rational{T}, y::FastRational) where T = x <= Rational{T}(y)
(<)(x::FastRational, y::Rational{T}) where T = Rational{T}(x) < y
(<=)(x::FastRational, y::Rational{T}) where T = Rational{T}(x) <= y

(==)(x::FastRational, y::FastRational) =
   numerator(x) === numerator(y) && denominator(x) === denominator(y)
(!=)(x::FastRational, y::FastRational) =
   numerator(x) !== numerator(y) || denominator(x) !== denominator(y)

for F in (:(>), :(>=), :(<=), :(<))
    @eval $F(x::FastRational, y::FastRational) =
   $F(numerator(x)*denominator(y), numerator(y)*denominator(x))
end

for F in (:(==), :(!=), :(>), :(>=), :(<=), :(<))
    @eval $F(x::FastRational, y::PlainRational) = $F(x, canonical(y))
end
for F in (:(==), :(!=), :(>), :(>=), :(<=), :(<))
    @eval $F(x::PlainRational, y::FastRational) = $F(canonical(x), y)
end
for F in (:(==), :(!=), :(>), :(>=), :(<=), :(<))
    @eval $F(x::PlainRational, y::PlainRational) = $F(canonical(x), canonical(y))
end

isequal(x::FastRational, y::FastRational) = x == y
isless(x::FastRational, y::FastRational) = x <= y
isequal(x::FastRational, y::PlainRational) = x == y
isless(x::FastRational, y::PlainRational) = x <= y
isequal(x::PlainRational, y::FastRational) = x == y
isless(x::PlainRational, y::FastRational) = x <= y
isequal(x::PlainRational, y::PlainRational) = x == y
isless(x::PlainRational, y::PlainRational) = x <= y

string(x::FastRational) = string(numerator(x), "//", denominator(x))
@inline string(x::PlainRational) = string(FastRational(x))

function show(io::IO, x::FastRational)
    print(io, string(x))
end
@inline show(io::IO, x::PlainRational) = show(io, FastRational(x))
show(x::FastRational) = show(Base.STDOUT, "text/plain", string(x))
@inline show(x::PlainRational) = show(FastRational(x))

function string(x::Vector{Q}) where Q<:FastRational
    str = string((x...,))
    str = string("FastRational[", str[2:end-1], "]")
    return str
end
string(x::Vector{Q}) where Q<:PlainRational = string(map(FastRational, x))

function show(io::IO, x::Vector{Q}) where Q<:Union{FastRational, PlainRational}
    print(io, string(x))
end
function show(x::Vector{Q}) where Q<:Union{FastRational, PlainRational}
    show(Base.STDOUT, "text/plain", string(x))
end

function read(io::IO, ::Type{T}) where T<:FastRational
    num = read(io,T)
    den = read(io,T)
    return FastRational(num, den)
end
@inline read(io::IO, x::PlainRational) = read(io, FastRational(x))

function write(io::IO, x::FastRational)
    return write(io, numerator(z), denominator(z))
end
@inline write(io::IO, x::PlainRational) = write(io, FastRational(x))

end # module