text:
# # Meridional diffusion

using Oceananigans
using Oceananigans.TurbulenceClosures: Horizontal
using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel, PrescribedVelocityFields

using Statistics
using JLD2
using Printf
using GLMakie

Nx = 360
Ny = 120
latitude = (-60, 60)
longitude = (-180, 180)

# A spherical domain
grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                             radius = 1,
                             latitude = latitude,
                             longitude = longitude,
                             z = (-1, 0))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = :c,
                                    velocities = PrescribedVelocityFields(), # quiescent
                                    closure = HorizontalScalarDiffusivity(κ=1),
                                    buoyancy = nothing)

# Tracer patch for visualization
Gaussian(λ, ϕ, L) = exp(-(λ^2 + ϕ^2) / 2L^2)

# Tracer patch parameters
L = 12 # degree
ϕ₀ = 0 # degrees

cᵢ(λ, ϕ, z) = Gaussian(λ, ϕ - ϕ₀, L)

set!(model, c=cᵢ)

c = model.tracers.c

function progress(s)
    c = s.model.tracers.c
    @info "Maximum(c) = $(maximum(c)), time = $(s.model.clock.time) / $(s.stop_time)"
    return nothing
end

ϕᵃᶜᵃ_max = maximum(abs, ynodes(Center, grid))
Δx_min = grid.radius * cosd(ϕᵃᶜᵃ_max) * deg2rad(grid.Δλ)
Δy_min = grid.radius * deg2rad(grid.Δϕ)
Δ_min = min(Δx_min, Δy_min)

# Time-scale for gravity wave propagation across the smallest grid cell
cell_diffusion_time_scale = Δ_min^2

simulation = Simulation(model,
                        Δt = 0.1cell_diffusion_time_scale,
                        stop_time = 1000cell_diffusion_time_scale,
                        iteration_interval = 100,
                        progress = progress)

output_fields = model.tracers

output_prefix = "spot_tracer_diffusion_Nx$(grid.Nx)_Ny$(grid.Ny)"

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = TimeInterval(10cell_diffusion_time_scale),
                                                      prefix = output_prefix,
                                                      force = true)

run!(simulation)

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

λ = xnodes(Center, grid)
ϕ = ynodes(Center, grid)

λ = repeat(reshape(λ, Nx, 1), 1, Ny)
ϕ = repeat(reshape(ϕ, 1, Ny), Nx, 1)

iter = Node(0)

plot_title = "hi"

c = @lift file["timeseries/c/" * string($iter)][:, :, 1]

set_theme!(Theme(fontsize = 30))

fig = Figure(resolution = (1920, 1080))

title = @lift "Tracer spot on a sphere, t = $(file["timeseries/t/" * string($iter)])"

ax = fig[1, 1] = Axis(fig,
                      xlabel = "λ",
                      ylabel = "ϕ",
                      title = title)

heatmap!(ax, c)

record(fig, "spot_tracer_diffusion_Nx$(grid.Nx)_Ny$(grid.Ny).mp4", iterations, framerate=30) do i
    @info "Animating iteration $i/$(iterations[end])..."
    iter[] = i
end


text:
######## Stan diagnose example  ###########

using Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

bernoulli = "
data {
  int<lower=0> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
    y ~ bernoulli(theta);
}
"

bernoullidata = [
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
  Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
]

stanmodel = Stanmodel(Diagnose(Gradient(epsilon=1e-6)), name="bernoulli", model=bernoulli);

diags = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME);
diags[1]["diagnose"] |> display

end # cd

ext: jl


text:
#  Copyright 2016-2019 Eric S. Tellez <eric.tellez@infotec.mx>
#
#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.

export IHCSearch

# Iterated Hill Climbing Search
struct IHCSearch <: LocalSearchAlgorithm
    restarts::Int
    use_local_improvement::Bool
    IHCSearch() = new(1, false)
    IHCSearch(r; use_local_improvement=false) = new(r, use_local_improvement)
end

function hill_climbing(index::SearchGraph{T}, dist::Function, q::T, res::KnnResult, navigation_state, nodeID::Int64, use_local_improvement::Bool) where T
    omin::Int = -1
    dmin::Float32 = typemax(Float32)
    # @info "STARTING HILL CLIMBING"
    while true
        dmin = typemax(Float32)
        omin = -1
        navigation_state[nodeID] = EXPLORED
        @inbounds for childID in index.links[nodeID]
            S = get(navigation_state, childID, UNKNOWN)
            S != UNKNOWN && continue
            navigation_state[childID] = VISITED
            d = convert(Float32, dist(index.db[childID], q))
            if use_local_improvement  ## this yields to better quality but can't be tuned for early stopping
                push!(res, childID, d)
                if d < dmin
                    dmin = d
                    omin = childID
                end
            else
                if push!(res, childID, d) && d < dmin
                    dmin = d
                    omin = childID
                end
            end
        end

        # @info dmin
        if omin < 0
            break
        else
            nodeID = omin
        end
    end
end

function search(isearch::IHCSearch, index::SearchGraph{T}, dist::Function, q::T, res::KnnResult, navigation_state, hints=EMPTY_INT_VECTOR) where T
    n = length(index.db)
    restarts = min(isearch.restarts, n)
    if length(hints) == 0
        hints = rand(1:n, isearch.restarts)
    end

    @inbounds for start_point in hints
        # start_point = rand(range)
        S = get(navigation_state, start_point, UNKNOWN)
        if S == UNKNOWN
            navigation_state[start_point] = VISITED
            d = convert(Float32, dist(q, index.db[start_point]))
            push!(res, start_point, d)
            hill_climbing(index, dist, q, res, navigation_state, start_point, isearch.use_local_improvement)
        end
    end

    res
end

function opt_expand_neighborhood(fun, algo::IHCSearch, n::Int, iter::Int)
    f_(w) = ceil(Int, w * (rand() - 0.5))
    f(x, w) = max(1, x + f_(w))
    logn = ceil(Int, log(2, n+1))
    w = 1

    while w <= logn  ## log log n
        IHCSearch(f(algo.restarts, w)) |> fun
        w += w
    end
end


text:
# Author: Lukasz Bratos

# Typy zmiennopozycyjne dla ktorych bedziemy wykonywac testy
types = [Float16, Float32, Float64]

# Funkcja przyjmuje typ zmiennopozycyjny i oblicza wartosc wyrazenia 3(4/3 - 1) - 1
function kahanEps(type)
    uno = one(type)
    return 3 * uno * (((4 * uno) / (3 * uno)) - uno) - uno
end

for type in types
    println("Calculated value of Khan's epsilon for $type: ", kahanEps(type), " vs eps($type): ", eps(type))
end


text:
## make the surface struct
using transD_GP, Random
surface = transD_GP.Surface(imagegrid;useML=false, σ=noisegrid)
## lengthscale GP options (stationary)
nminlog10λ, nmaxlog10λ = 2, 150
pnorm = 2.
Klog10λ = transD_GP.GP.Mat32()
log10bounds = [log10(20_000) log10(300_000); log10(20_000) log10(300_000)]
δlog10λ = 0.2
# making sure it is not anisotropic ...
λlog10λ = 0.1*abs.([diff([extrema(xall[2,:])...])[1], diff([extrema(xall[2,:])...])[1]])
demean = false # essential for length scale change approximations
sdev_proplog10λ = 0.05*diff(log10bounds, dims=2)[:]
sdev_poslog10λ = 0.05*abs.([diff([extrema(xall[1,:])...])[1], diff([extrema(xall[2,:])...])[1]])
## Initialize a lengthscale model using these options
Random.seed!(12)
optlog10λ = transD_GP.OptionsStat(nmin = nminlog10λ,
                        nmax = nmaxlog10λ,
                        xbounds = [minimum(xall[1,:]) maximum(xall[1,:]);
                                   minimum(xall[2,:]) maximum(xall[2,:])],
                        fbounds = log10bounds,
                        xall = xall,
                        λ = λlog10λ,
                        δ = δlog10λ,
                        needλ²fromlog = true,
                        updatenonstat = true,
                        demean = demean,
                        sdev_prop = sdev_proplog10λ,
                        sdev_pos = sdev_poslog10λ,
                        pnorm = pnorm,
                        quasimultid = false,
                        K = Klog10λ,
                        timesλ = 4.
                        )
## make options for the nonstationary actual properties GP
nmin, nmax = 2, 150
ftrain = surface.d[surface.select]
fbounds = permutedims([extrema(ftrain)...])
sdev_prop = 0.05*diff(fbounds, dims=2)[:]
sdev_pos = 0.05*abs.([diff([extrema(xall[1,:])...])[1], diff([extrema(xall[2,:])...])[1]])
demean_ns = false
sampledc = true
δ = 1.
K = transD_GP.GP.Mat32()
## Initialize model for the nonstationary properties GP
Random.seed!(13)
opt = transD_GP.OptionsNonstat(optlog10λ,
                        nmin = nmin,
                        nmax = nmax,
                        fbounds = fbounds,
                        δ = δ,
                        demean = demean_ns,
                        sampledc = sampledc,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        K = K
                        )


text:
export ft_resetport

function ft_resetport(ft_handle::Culong)
  ft_status = ccall((:FT_ResetPort, d2xx),Cuint,(Culong,),ft_handle)
  checkstatus(ft_status)
  return nothing
end


text:
# Part of submodule Utils of BetaML _ the Beta Machine Learning Toolkit
# Vatious utils to help in logging/debugging

@enum Verbosity NONE=0 LOW=10 STD=20 HIGH=30 FULL=40

"""
    @codeLocation()

Helper macro to print during runtime an info message concerning the code being executed position

"""
macro codeLocation()
    return quote
        st = stacktrace(backtrace())
        myf = ""
        for frm in st
            funcname = frm.func
            if frm.func != :backtrace && frm.func!= Symbol("macro expansion")
                myf = frm.func
                break
            end
        end
        println("Running function ", $("$(__module__)"),".$(myf) at ",$("$(__source__.file)"),":",$("$(__source__.line)"))
        println("Type `]dev BetaML` to modify the source code (this would change its location on disk)")
    end
end


text:
# TODO: I declare lots of types here because of the lacking ability to do forward declarations in current Julia
# I should move these to the relevant files when something like "extern" is implemented

const AVec = AbstractVector
const AMat = AbstractMatrix
const KW = Dict{Symbol,Any}
const AKW = AbstractDict{Symbol,Any}
const TicksArgs = Union{AVec{T}, Tuple{AVec{T}, AVec{S}}, Symbol} where {T<:Real, S<:AbstractString}

struct PlotsDisplay <: AbstractDisplay end

# -----------------------------------------------------------

struct InputWrapper{T}
    obj::T
end
wrap(obj::T) where {T} = InputWrapper{T}(obj)
Base.isempty(wrapper::InputWrapper) = false

# -----------------------------------------------------------

mutable struct Series
    plotattributes::DefaultsDict
end

attr(series::Series, k::Symbol) = series.plotattributes[k]
attr!(series::Series, v, k::Symbol) = (series.plotattributes[k] = v)

# -----------------------------------------------------------

# a single subplot
mutable struct Subplot{T<:AbstractBackend} <: AbstractLayout
    parent::AbstractLayout
    series_list::Vector{Series}  # arguments for each series
    minpad::Tuple # leftpad, toppad, rightpad, bottompad
    bbox::BoundingBox  # the canvas area which is available to this subplot
    plotarea::BoundingBox  # the part where the data goes
    attr::DefaultsDict  # args specific to this subplot
    o  # can store backend-specific data... like a pyplot ax
    plt  # the enclosing Plot object (can't give it a type because of no forward declarations)
end

Base.show(io::IO, sp::Subplot) = print(io, "Subplot{$(sp[:subplot_index])}")

# -----------------------------------------------------------

# simple wrapper around a KW so we can hold all attributes pertaining to the axis in one place
mutable struct Axis
    sps::Vector{Subplot}
    plotattributes::DefaultsDict
end

mutable struct Extrema
    emin::Float64
    emax::Float64
end
Extrema() = Extrema(Inf, -Inf)

# -----------------------------------------------------------

const SubplotMap = Dict{Any, Subplot}

# -----------------------------------------------------------

mutable struct Plot{T<:AbstractBackend} <: AbstractPlot{T}
    backend::T                   # the backend type
    n::Int                       # number of series
    attr::DefaultsDict            # arguments for the whole plot
    series_list::Vector{Series}  # arguments for each series
    o                            # the backend's plot object
    subplots::Vector{Subplot}
    spmap::SubplotMap            # provide any label as a map to a subplot
    layout::AbstractLayout
    inset_subplots::Vector{Subplot}  # list of inset subplots
    init::Bool
end

function Plot()
    Plot(backend(), 0, DefaultsDict(KW(), _plot_defaults), Series[], nothing,
         Subplot[], SubplotMap(), EmptyLayout(),
         Subplot[], false)
end

# -----------------------------------------------------------------------

Base.getindex(plt::Plot, i::Integer) = plt.subplots[i]
Base.length(plt::Plot) = length(plt.subplots)
Base.lastindex(plt::Plot) = length(plt)

Base.getindex(plt::Plot, r::Integer, c::Integer) = plt.layout[r,c]
Base.size(plt::Plot) = size(plt.layout)
Base.size(plt::Plot, i::Integer) = size(plt.layout)[i]
Base.ndims(plt::Plot) = 2

# attr(plt::Plot, k::Symbol) = plt.attr[k]
# attr!(plt::Plot, v, k::Symbol) = (plt.attr[k] = v)

Base.getindex(sp::Subplot, i::Integer) = series_list(sp)[i]
Base.lastindex(sp::Subplot) = length(series_list(sp))

# -----------------------------------------------------------------------


text:
export InfiniteSumNode, slicesample!

"""
    InfiniteSumNode

An infinite mixture over Bayesian SPNs realized using a stick-breaking representation of a Dirichlet process.

Parameters:

* id::Symbol                        Id
* α::Float64                        Concentration parameter
* w::Vector{<:Real}                 Infinite collection of weights
* z::Vector{Int}                    Assignments of observations
* children::Vector                  Infinite collection of Bayesian SPNs
* ustar::Float64                    Minimum slice variable
* istar::Int                        Observation inducing minimum slice variable
* kstar::Int                        Cluster inducing minimum slice variable
* logpdfbase::Function              Log pdf of base-case: f(x::Matrix) -> Vector{Float64}
* spnprior::Function                Function to draw a region-graph SPN.

"""
mutable struct InfiniteSumNode

    id::Symbol
    α::Float64
    w::Vector{<:Real}
    z::Vector{Int}
    children::Vector{<:AbstractRegionGraphNode}
    ustar::Float64
    istar::Int
    kstar::Int
    logpbase::Function
    spnprior::Function

end

function logpdf(n::InfiniteSumNode, x::AbstractMatrix{T}) where {T<:Real}
    @assert !isempty(n.w)
    lp_ = logpdf.(n.children, Ref(x))
    lp_ = reduce(hcat, lp_)
    return lse(lp_ .+ log.(n.w[1:end-1])', dims=2)
end

## ############## ##
## Slice sampling ##
## ############## ##

"""
    slicesample!(spn::InfiniteSumNode, x::AbstractMatrix{<:Real})

Distributed slice sampler for infinite mixtures of Bayesian SPNs based on Hong Ge et al. Distributed inference for Dirichlet process mixture models, ICML, 2015.
"""
function slicesample!(spn::InfiniteSumNode, x::AbstractMatrix{<:Real})

    (N, D) = size(x)

    if isempty(spn.children)
        # New SPN.
        spn.z[:] = ones(Int, N)
    else
        _slicesample!(spn, x)
    end

    # Update parameters
    K = max(length(spn.children), maximum(spn.z))
    n = map(k -> sum(spn.z .== k), 1:K)
    active = findall(n .!= 0)
    n = map(Float64, n)

    M = sum(n .== 0) + 1
    n[n .== 0] .= spn.α/M
    spn.w = rand(Dirichlet(vcat(n, spn.α/M)))

    # Add new children?
    if K > length(spn.children)
        push!(spn.children, spn.spnprior())
    end

    for (k, child) in enumerate(spn.children)
        ancestralsampling!(child, x; obs = findall(spn.z .== k))
    end

    b = map(k -> rand(Beta(1, n[k])), active)
    u = spn.w[active] .* b

    spn.ustar = minimum(u)
    spn.kstar = active[argmin(u)]
    spn.istar = rand(findall(spn.z .== spn.kstar))

    return spn
end

function _slicesample!(spn::InfiniteSumNode, x::AbstractMatrix{<:Real})
    (N,D) = size(x)

    lpbase = spn.logpbase(x)
    lp = map(child -> logpdf(child, x), spn.children)

    Threads.@threads for i in 1:N
        local w = spn.w[spn.z[i]]
        local u = i == spn.istar ? spn.ustar : rand(Uniform(spn.ustar, w))

        # Resample z_i.
        local ids = findall(spn.w .>= u)
        local K = length(spn.children)
        local p = zeros(length(ids))
        for (j, k) in enumerate(ids)
            if k > K
                p[j] = exp(lpbase[i])
            else
                p[j] = exp(lp[k][i])
            end
        end

        spn.z[i] = randomindex(p)
    end
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
using MolecularGraph: atomsymbol!, atomprop!, atom!

@testset "smarts.atom" begin

@testset "atomsymbol" begin
    # outside of []
    state = SmartsParser("Cl", false)
    Chloride = atomsymbol!(state)
    @test Chloride == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :Cl),
        QueryFormula(:isaromatic, false)
    ]))
    @test state.pos == 3

    state = SmartsParser("Cr", false)
    Chromium = atomsymbol!(state)
    @test Chromium == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :C),
        QueryFormula(:isaromatic, false)
    ]))
    @test state.pos == 2

    state = SmartsParser("p", false)
    aromp = atomsymbol!(state)

    @test aromp == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :P),
        QueryFormula(:isaromatic, true)
    ]))

    state = SmartsParser("*", false)
    anyatom = atomsymbol!(state)
    @test anyatom == QueryFormula(:any, true)
end

@testset "atomprop" begin
    # inside of []
    state = SmartsParser("s", false)
    aroms = atomprop!(state)
    @test aroms == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :S),
        QueryFormula(:isaromatic, true)
    ]))

    state = SmartsParser("123", false)
    iso = atomprop!(state)
    @test iso == QueryFormula(:mass, 123)

    state = SmartsParser("H41", false)
    H4 = atomprop!(state)
    @test H4 == QueryFormula(:hydrogenconnected, 4)
    @test state.pos == 3

    state = SmartsParser("X2", false)
    X2 = atomprop!(state)
    @test X2 == QueryFormula(:connectivity, 2)

    state = SmartsParser("Xe", false)
    Xe = atomprop!(state)
    @test Xe == QueryFormula(:atomsymbol, :Xe)
    @test state.pos == 3

    state = SmartsParser("Na", false)
    Na = atomprop!(state)
    @test Na == QueryFormula(:atomsymbol, :Na)
    @test state.pos == 3

    state = SmartsParser("Yv2", false)
    Yval = atomprop!(state)
    @test Yval == QueryFormula(:atomsymbol, :Y)

    state = SmartsParser("+23", false)
    chg1 = atomprop!(state)
    @test chg1 == QueryFormula(:charge, 2)
    @test state.pos == 3

    state = SmartsParser("----+", false)
    chg2 = atomprop!(state)
    @test chg2 == QueryFormula(:charge, -4)
    @test state.pos == 5

    state = SmartsParser("@+", false)
    stereo1 = atomprop!(state)
    @test stereo1 == QueryFormula(:stereo, :anticlockwise)
    @test state.pos == 2

    state = SmartsParser("@@?", false)
    stereo4 = atomprop!(state)
    @test stereo4 == QueryFormula(:not, QueryFormula(:stereo, :anticlockwise))
    @test state.pos == 4

    state = SmartsParser("\$([CH2]=*)", false)
    rec = atomprop!(state)
    @test rec == QueryFormula(:recursive, "[CH2]=*")
    @test state.pos == 11
end

@testset "atom" begin
    state = SmilesParser("Br", false)
    br = atom!(state)[1]
    @test br.symbol == :Br

    state = SmilesParser("[2H]", false)
    hyd = atom!(state)
    @test length(hyd) == 1
    @test hyd[1].mass == 2

    state = SmilesParser("[H2]", false)
    hmol = atom!(state)
    @test length(hmol) == 2

    state = SmilesParser("[H+]", false)
    proton = atom!(state)[1]
    @test proton.charge == 1

    state = SmilesParser("[14c@@H]", false)
    iso = atom!(state)
    @test iso[1].symbol == :C
    @test iso[1].mass == 14
    @test iso[1].isaromatic
    @test iso[1].stereo == :clockwise
    @test iso[2].symbol == :H
    @test length(iso) == 2
    @test state.pos == 9

    state = SmilesParser("[NH4+]", false)
    ammonium = atom!(state)
    @test length(ammonium) == 5
    @test ammonium[1].charge == 1
    @test ammonium[5].symbol == :H

    state = SmilesParser("[bH+]", false)
    aromB = atom!(state)
    @test aromB[1].symbol == :B
    @test aromB[1].isaromatic
    @test aromB[1].charge == 1
    @test aromB[2].symbol == :H

    state = SmilesParser("[as]", false)
    aromAs = atom!(state)[1]
    @test aromAs.symbol == :As
    @test aromAs.isaromatic

    state = SmilesParser("[Zn++]", false)
    zn = atom!(state)[1]
    @test zn.symbol == :Zn
    @test zn.charge == 2

    state = SmilesParser("[O-2]", false)
    ox = atom!(state)[1]
    @test ox.symbol == :O
    @test ox.charge == -2

end

@testset "smartsatom" begin
    state = SmartsParser("c", false)
    aromc = atom!(state)[1]
    @test aromc.query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :C),
        QueryFormula(:isaromatic, true)
    ]))

    state = SmartsParser("a", false)
    anyarom = atom!(state)[1]
    @test anyarom.query == QueryFormula(:isaromatic, true)

    state = SmartsParser("[]", false)
    @test_throws AssertionError atom!(state)

    state = SmartsParser("[*]", false)
    anyatom = atom!(state)[1]
    @test anyatom.query == QueryFormula(:any, true)

    state = SmartsParser("[#16]", false)
    no16 = atom!(state)[1]
    @test no16.query == QueryFormula(:atomsymbol, :S)

    state = SmartsParser("[CH2]", false)
    methylene = atom!(state)[1]
    @test methylene.query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :C),
        QueryFormula(:isaromatic, false),
        QueryFormula(:hydrogenconnected, 2)
    ]))

    state = SmartsParser("[!C;R]", false)
    ringnotalp = atom!(state)[1]
    @test ringnotalp.query == QueryFormula(:and, Set([
        QueryFormula(:or, Set([
            QueryFormula(:not, QueryFormula(:atomsymbol, :C)),
            QueryFormula(:isaromatic, true)
        ])),
        QueryFormula(:not, QueryFormula(:sssrcount, 0))
    ]))

    state = SmartsParser("[n&H1]", false)
    nh1 = atom!(state)[1]
    @test nh1.query == QueryFormula(:and, Set([
        QueryFormula(:atomsymbol, :N),
        QueryFormula(:isaromatic, true),
        QueryFormula(:hydrogenconnected, 1)
    ]))

    state = SmartsParser("[*r6]", false)
    sixmem = atom!(state)[1]
    @test sixmem.query == QueryFormula(:smallestsssr, 6)

    state = SmartsParser("[35*]", false)
    any35 = atom!(state)[1]
    @test any35.query == QueryFormula(:mass, 35)

    state = SmartsParser("[F,Cl,Br,I]", false)
    fourhalo = atom!(state)[1]
    @test fourhalo.query == QueryFormula(:or, Set([
        QueryFormula(:atomsymbol, :F),
        QueryFormula(:atomsymbol, :Cl),
        QueryFormula(:atomsymbol, :Br),
        QueryFormula(:atomsymbol, :I)
    ]))
end

end # smiles.atom


text:
using FileIO
using JLD2

include("usefulFunctions.jl");
include("DoME.jl");

import Base.OneTo

datasetType=Float64;

Strategies = [StrategyExhaustive, StrategyExhaustiveWithConstantOptimization, StrategySelectiveWithConstantOptimization, StrategySelective];
strategy = 1;

configurations = Array{Int,3}(undef,length(MinimumReductionsMSE),length(MaxNumNodes),numFolds);
n=1; for i1=1:length(MinimumReductionsMSE), i2=1:length(MaxNumNodes), i3=1:numFolds global n; configurations[i1,i2,i3] = n; n += 1; end;
configuration(number::Int) = ( index=findfirst(configurations.==number); (index==nothing)&&(exit()); (MinimumReductionsMSE[index[1]],MaxNumNodes[index[2]],index[3]) );
configurationIndices(minimumReductionMSE, maxNumNodes, numFold) = (findfirst(MinimumReductionsMSE.==minimumReductionMSE), findfirst(MaxNumNodes.==maxNumNodes), numFold);

SourceCodeCompiled = false;

println("Number of executions: ", n-1);

arguments = parseARGS();
if isempty(arguments)
    datasetName = datasetNames()[1];
    configurationsRun = 1:(length(Strategies)*length(MinimumReductionsMSE)*length(MaxNumNodes)*numFolds);
else
    @assert(length(arguments)>=2);
    datasetName = arguments[1];
    strategy = arguments[2];
    if length(arguments)==3
        configurationsRun = arguments[3];
    else
        configurationsRun = 1:(length(Strategies)*length(MinimumReductionsMSE)*length(MaxNumNodes)*numFolds);
    end;
end;

# If the folder "results" does no exist, create it
(!isdir("results")) && mkpath("results");

function loadResults(datasetName, strategy, minimumReductionMSE, maxNumNodes, numFold)
    fileName = resultsFileName(datasetName, strategy);
    results = load(fileName, "results");
    (i1, i2, i3) = configurationIndices(minimumReductionMSE, maxNumNodes, numFold);

    if i1>size(results,1) || i2>size(results,2) || i3>size(results,3)
        return [];
    else
        executionResults = results[i1, i2, i3, :];
        return (any(isnan.(executionResults[[1,3,4,5]])) || any(executionResults[3:end].==0)) ? [] : executionResults;
    end;
end;

function saveResults(datasetName, strategy, minimumReductionMSE, maxNumNodes, numFold, executionResults)
    fileName = resultsFileName(datasetName, strategy);
    results = load(fileName, "results");
    (i1, i2, i3) = configurationIndices(minimumReductionMSE, maxNumNodes, numFold);
    if i1>size(results,1) || i2>size(results,2) || i3>size(results,3)
        newResults = Array{Float64,4}(undef, max(i1,size(results,1)), max(i2,size(results,2)), max(i3,size(results,3)), 5);
        newResults[:] .= NaN;
        if !isempty(results)
            newResults[1:size(results,1), 1:size(results,2), 1:size(results,3), :] .= results[:,:,:,:];
        end;
        results = newResults;
    end;
    results[i1, i2, i3, :] .= executionResults;
    save(fileName, "results", results);
end;

# The function to test
function functionTest(inputs, targets, indicesTest, minimumReductionMSE, maxNumNodes, strategy)
    GC.gc();
    @timed dome(inputs, targets;
        testIndices = indicesTest,
        minimumReductionMSE = minimumReductionMSE,
        maximumHeight = Inf ,
        maximumNodes = maxNumNodes ,
        strategy = Strategies[strategy] ,
        showText = false
    );
end;

if (!isfile(resultsFileName(datasetName, strategy)))
    save(resultsFileName(datasetName, strategy), "results", Array{Float64,4}(undef,0,0,0,0));
end;

(inputs, targets) = loadDataset(datasetName; datasetType=datasetType)

# Create the same cross-validation indices used in the rest of experiments
indicesKFold = crossvalidationIndices(length(targets), numFolds);

println("I'm going to run the following configurations:")
for currentConfiguration in configurationsRun

    (minimumReductionMSE, maxNumNodes, numFold) = configuration(currentConfiguration)
    println("   ", currentConfiguration, ": k-fold with strategy ", strategy, ", minimum reduction in MSE: ", minimumReductionMSE, ", maximum number of nodes: ", maxNumNodes, ", fold ", numFold);

    if isempty(loadResults(datasetName, strategy, minimumReductionMSE, maxNumNodes, numFold))

            global SourceCodeCompiled;
            if (!SourceCodeCompiled)
                # Compile with a "fast" execution
                functionTest(inputs, targets, collect(1:Int(trunc(length(targets)/10))), 1e-4, 50, strategy);
                SourceCodeCompiled = true;
            end;

            ((trainingMSE, validationMSE, testMSE, bestTree), elapsedTime) = functionTest(inputs, targets, findall(indicesKFold.==numFold), minimumReductionMSE, maxNumNodes, strategy);

            executionResults = [trainingMSE, testMSE, elapsedTime, Float64(height(bestTree)), Float64(numNodes(bestTree))];
            saveResults(datasetName, strategy, minimumReductionMSE, maxNumNodes, numFold, executionResults);

            println("      Finished training fold ", numFold, "/", numFolds, " MSE training/test: ", trainingMSE, "/", testMSE, " in ", elapsedTime, " seconds");

    end;

end;


text:
using GeometricIntegrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Utils
using SimpleSolvers
using Test

import GeometricProblems.LotkaVolterra4d
import GeometricProblems.LotkaVolterra4dLagrangian
import GeometricProblems.LotkaVolterra4dLagrangian: Δt, nt, reference_solution

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

SimpleSolvers.set_config(:nls_atol_break, 1E3)
SimpleSolvers.set_config(:nls_rtol_break, 1E3)
SimpleSolvers.set_config(:nls_stol_break, 1E3)

@testset "$(rpad("Lotka-Volterra 4D",80))" begin
    ode  = LotkaVolterra4dLagrangian.lotka_volterra_4d_ode()
    iode = LotkaVolterra4dLagrangian.lotka_volterra_4d_iode()
    idae = LotkaVolterra4dLagrangian.lotka_volterra_4d_idae()

    ref_ode  = LotkaVolterra4d.lotka_volterra_4d_ode()
    ref_iode = LotkaVolterra4d.lotka_volterra_4d_iode()
    ref_idae = LotkaVolterra4d.lotka_volterra_4d_idae()

    ref_equs = get_functions(ref_ode)
    ref_invs = get_invariants(ref_ode)

    v1 = zero(ode.q₀[begin])
    v2 = zero(ref_ode.q₀[begin])
    ode.v(ode.t₀, ode.q₀[begin], v1)
    ref_equs[:v](ref_ode.t₀, ref_ode.q₀[begin], v2)
    @test v1 ≈ v2  atol=1E-14

    h1 = ode.invariants[:h](ode.t₀, ode.q₀[begin])
    h2 = ref_invs[:h](ref_ode.t₀, ref_ode.q₀[begin])
    @test h1 ≈ h2  atol=1E-14

    int = Integrator(ode, TableauGauss(2), Δt)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 8E-4

    ref_int = Integrator(ref_ode, TableauGauss(2), Δt)
    ref_sol = integrate(ref_ode, ref_int, nt)
    @test relative_maximum_error(sol.q, ref_sol.q[end]) < 2E-14

    ref_equs = get_functions(ref_iode)
    ref_invs = get_invariants(ref_iode)

    @test iode.p₀ == ref_iode.p₀

    ϑ1 = zero(iode.q₀[begin])
    ϑ2 = zero(ref_iode.q₀[begin])
    iode.ϑ(iode.t₀, iode.q₀[begin], v1, ϑ1)
    ref_equs[:ϑ](ref_iode.t₀, ref_iode.q₀[begin], v2, ϑ2)
    @test ϑ1 ≈ ϑ2  atol=1E-14

    f1 = zero(iode.q₀[begin])
    f2 = zero(ref_iode.q₀[begin])
    iode.f(iode.t₀, iode.q₀[begin], v1, f1)
    ref_equs[:f](ref_iode.t₀, ref_iode.q₀[begin], v2, f2)
    @test f1 ≈ f2  atol=1E-14

    g1 = zero(iode.q₀[begin])
    g2 = zero(ref_iode.q₀[begin])
    iode.g(iode.t₀, iode.q₀[begin], v1, g1)
    ref_equs[:g](ref_iode.t₀, ref_iode.q₀[begin], v2, g2)
    @test g1 ≈ g2  atol=1E-14

    v1 = zero(iode.q₀[begin])
    v2 = zero(ref_iode.q₀[begin])
    iode.v̄(iode.t₀, iode.q₀[begin], v1)
    ref_equs[:v̄](ref_iode.t₀, ref_iode.q₀[begin], v2)
    @test v1 ≈ v2  atol=1E-14

    f1 = zero(iode.q₀[begin])
    f2 = zero(ref_iode.q₀[begin])
    iode.f̄(iode.t₀, iode.q₀[begin], v1, f1)
    ref_equs[:f̄](ref_iode.t₀, ref_iode.q₀[begin], v2, f2)
    @test f1 ≈ f2  atol=1E-14

    h1 = iode.invariants[:h](iode.t₀, iode.q₀[begin], zero(iode.q₀[begin]))
    h2 = ref_invs[:h](ref_iode.t₀, ref_iode.q₀[begin], zero(ref_iode.q₀[begin]))
    @test h1 ≈ h2  atol=1E-14

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 2E-3

    ref_int = IntegratorVPRKpMidpoint(ref_iode, TableauVPGLRK(2), Δt)
    ref_sol = integrate(ref_iode, ref_int, nt)
    @test relative_maximum_error(sol.q, ref_sol.q[end]) < 8E-14

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 8E-4

    ref_int = IntegratorVPRKpSymmetric(ref_iode, TableauVPGLRK(2), Δt)
    ref_sol = integrate(ref_iode, ref_int, nt)
    @test relative_maximum_error(sol.q, ref_sol.q[end]) < 4E-14

    ref_equs = get_functions(ref_idae)
    ref_invs = get_invariants(ref_idae)

    @test idae.p₀ == ref_idae.p₀

    v1 = zero(idae.q₀[begin])
    v2 = zero(ref_idae.q₀[begin])
    idae.v̄(idae.t₀, idae.q₀[begin], v1)
    ref_equs[:v̄](ref_idae.t₀, ref_idae.q₀[begin], v2)
    @test v1 ≈ v2  atol=1E-14

    f1 = zero(idae.q₀[begin])
    f2 = zero(ref_idae.q₀[begin])
    idae.f̄(idae.t₀, idae.q₀[begin], v1, f1)
    ref_equs[:f̄](ref_idae.t₀, ref_idae.q₀[begin], v2, f2)
    @test f1 ≈ f2  atol=1E-14

    h1 = idae.invariants[:h](idae.t₀, idae.q₀[begin], zero(idae.q₀[begin]))
    h2 = ref_invs[:h](ref_idae.t₀, ref_idae.q₀[begin], zero(ref_idae.q₀[begin]))
    @test h1 ≈ h2  atol=1E-14

    int = Integrator(idae, TableauVSPARKGLRKpMidpoint(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 2E-3

    ref_int = Integrator(ref_idae, TableauVSPARKGLRKpMidpoint(2), Δt)
    ref_sol = integrate(ref_idae, ref_int, nt)
    @test relative_maximum_error(sol.q, ref_sol.q[end]) < 4E-14

    int = Integrator(idae, TableauVSPARKGLRKpSymmetric(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 8E-4

    ref_int = Integrator(ref_idae, TableauVSPARKGLRKpSymmetric(2), Δt)
    ref_sol = integrate(ref_idae, ref_int, nt)
    @test relative_maximum_error(sol.q, ref_sol.q[end]) < 8E-14

end


text:
include("PETSc.jl")

import MPI
MPI.Init()

PetscInitialize(["-ksp_monitor","-malloc","-malloc_debug","-malloc_dump"]);
#
#   Create a vector and put values in it

comm = MPI.COMM_WORLD
comm_size = MPI.Comm_size(MPI.COMM_WORLD)
comm_rank = MPI.Comm_rank(MPI.COMM_WORLD)

println("my comm rank = ", comm_rank)

b = PetscVec(comm);
PetscVecSetType(b,"mpi");
PetscVecSetSizes(b,10, comm_size*10);

for i=1:10
  idxm = [(comm_rank)*10 + i]  # index
  PetscVecSetValues(b,idxm,[10.0 + i],PETSC_ADD_VALUES);
end

#PetscVecSetValues(b,Array(1.0:10.));
println("assigned first set of values")
#PetscVecSetValues(b,[1,2],[11.5,12.5],PETSC_ADD_VALUES);
println("assigned second set of values")
PetscVecAssemblyBegin(b);
println("began vector assembly")
PetscVecAssemblyEnd(b);
println("finished vector assembly")
PetscView(b);

println("finished assembling vector")
#
#    Create a matrix
A = PetscMat(comm);
PetscMatSetType(A,"mpiaij");
PetscMatSetSizes(A,10,10,comm_size*10,comm_size*10);
PetscSetUp(A);
for i=1:10
  idxm = [(comm_rank)*10 + i]  # row index
  idxn = [(comm_rank)*10 + i]  # column index
  PetscMatSetValues(A,idxm, idxn,[10.0 + i],PETSC_ADD_VALUES);
end

println("finished setting matrix values")
PetscMatAssemblyBegin(A,PETSC_MAT_FINAL_ASSEMBLY);
PetscMatAssemblyEnd(A,PETSC_MAT_FINAL_ASSEMBLY);
println("finished assembling matrix")
PetscView(A);

println("finished viewing matrix")


text:
using GPCRAnalysis
using MIToS
using Test

# skip the network-hitting components by setting `skip_download = true` in the global namespace

@testset "GPCRAnalysis.jl" begin
    @test isempty(detect_ambiguities(GPCRAnalysis))
    @testset "utils" begin
        @test @inferred(GPCRAnalysis.residue_range("Q8VGW6_MOUSE/31-308")) == 31:308
        @test @inferred(GPCRAnalysis.strip_residue_range("Q8VGW6_MOUSE/31-308")) == "Q8VGW6_MOUSE"
        @test @inferred(species("Q8VGW6_MOUSE/31-308")) == "MOUSE"
        @test @inferred(uniprotX("Q8VGW6_MOUSE/31-308")) == "Q8VGW6"
        # Examples from https://www.uniprot.org/help/accession_numbers
        @test @inferred(uniprotX("A2BC19")) == "A2BC19"
        @test @inferred(uniprotX("P12345")) == "P12345"
        @test @inferred(uniprotX("A0A023GPI8")) == "A0A023GPI8"
        # X identifiers have strict naming conventions (among other things, only lengths 1-5, 6, or 11 are allowed)
        @test_throws ErrorException uniprotX("Q8VGW67_MOUSE/31-308")
        @test_throws ErrorException uniprotX("q8vgw6_MOUSE/31-308")
        @test_throws ErrorException uniprotX("QV8GW6_MOUSE/31-308")
    end
    if !isdefined(@__MODULE__, :skip_download) || !skip_download
        @testset "AlphaFold" begin
            @test try_download_alphafold("garbage") === nothing
            fn = tempname()
            @test try_download_alphafold("K7N608", fn) == fn
            @test isfile(fn)
            @test getchain(fn) isa AbstractVector{MIToS.PDB.PDBResidue}
            rm(fn)
        end
    end
    @testset "MSA" begin
        # The test file is copied from MIToS/test/data, with gratitude
        pf09645_sto = "PF09645_full.stockholm"
        msa = read(pf09645_sto, MIToS.Pfam.Stockholm)
        @test MIToS.MSA.nsequences(filter_species!(deepcopy(msa), "ATV")) == 1
        @test MIToS.MSA.nsequences(filter_long!(deepcopy(msa), 70)) == 3

        idx = SequenceMapping([0, 4, 5, 0])
        seqvals = fill(NaN, 9)
        seqvals[idx] = [0.1, 0.2, 0.3, 0.4]
        @test seqvals[4] == 0.2
        @test seqvals[5] == 0.3
        @test all(isnan, seqvals[1:3])
        @test all(isnan, seqvals[6:end])
    end
    @testset "analyze" begin
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -42.99754755332367 <= q <= 40.38858433100043)

                     add_NL_constraint(m, :(gelu(-0.8925450194703983 + 0.6110890427800215*gelu(-0.16626426792027038 + 0.44845618413236066*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.1740395626916258*gelu(0.8063063231762757 + 0.4905142959916642*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8017072131372296*gelu(0.9210510999614945 + -0.4935968211971775*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7966067754659663*gelu(-0.7328670881809516 + 0.6102429954796835*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + gelu(-0.7008872510227278 + -0.258173144464545*gelu(-0.16626426792027038 + 0.44845618413236066*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.49550787111793193*gelu(0.8063063231762757 + 0.4905142959916642*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8462943069272191*gelu(0.9210510999614945 + -0.4935968211971775*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.05869131996649157*gelu(-0.7328670881809516 + 0.6102429954796835*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + gelu(-0.2643911691309153 + -0.31660511055293883*gelu(-0.16626426792027038 + 0.44845618413236066*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.09136517532701127*gelu(0.8063063231762757 + 0.4905142959916642*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.36400896662330684*gelu(0.9210510999614945 + -0.4935968211971775*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.4876638787119827*gelu(-0.7328670881809516 + 0.6102429954796835*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + gelu(0.6679332747384503 + -0.6061375406574494*gelu(-0.16626426792027038 + 0.44845618413236066*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7532640510772808*gelu(0.8063063231762757 + 0.4905142959916642*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.9959397069462161*gelu(0.9210510999614945 + -0.4935968211971775*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8544693505755525*gelu(-0.7328670881809516 + 0.6102429954796835*gelu(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*gelu(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*gelu(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*gelu(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@reset
@testset "Coupled Cluster" begin
    wf = @energy CCSD;
    @test isapprox(wf.CorrelationEnergy, -75.0187095714, rtol=tol)
    wf = @energy CCSD(T);
    @test isapprox(wf.correction, -7.38086348153572e-5, rtol=tol)
    wf = @energy bccd;
    @test isapprox(wf.CorrelationEnergy, -75.01901343317104, rtol=tol)
    wf = @energy ecCCSD;
    @test isapprox(wf.CorrelationEnergy, -75.0187251416, rtol=tol)
    wf = @energy ecCCSD(T);
    @test isapprox(wf.correction, -7.382789470598815e-5, rtol=tol)
    @set ci_alg sparse
    wf = @energy ecCCSD;
    @test isapprox(wf.CorrelationEnergy, -75.0188434735, rtol=tol)
end


text:
#These types must be the same as in the QDLDL header
c_int   = Int64;
c_float = Float64;
c_bool  = Uint8;

function qdMatExample()

    i = [1, 7, 2, 3, 6, 10, 2, 3, 10, 4, 8, 5, 2, 6, 1, 7, 9, 4, 8, 7, 9, 2, 3, 10]
    j = [1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5, 6, 6, 7, 7, 7, 8, 8, 9, 9, 10, 10, 10]
    s = [1.000000e+00, 3.503210e-01, 4.606408e-01, -1.211894e-01, -2.900576e-02, 1.824522e-01, -1.211894e-01, 4.179283e-01, -1.565056e+00, 1.778279e-01, -8.453948e-02, 1.000000e-01, -2.900576e-02, -1.000000e+00, 3.503210e-01, -4.410924e-01, 1.786632e-01, -8.453948e-02, -3.162278e-01, 1.786632e-01, -2.990768e-01, 1.824522e-01, -1.565056e+00, -1.000000e-01]

    A = sparse(i,j,s)
end

ldllibH = Libdl.dlopen("../../build/out/libqdldl." * Base.Libdl.dlext)
QDLDL_etree    = Libdl.dlsym(ldllibH, :QDLDL_etree)
QDLDL_factor   = Libdl.dlsym(ldllibH, :QDLDL_factor)
QDLDL_Lsolve   = Libdl.dlsym(ldllibH, :QDLDL_Lsolve)
QDLDL_Ltsolve  = Libdl.dlsym(ldllibH, :QDLDL_Ltsolve)
QDLDL_solve    = Libdl.dlsym(ldllibH, :QDLDL_solve)

A = qdMatExample()::SparseMatrixCSC{c_float,c_int}
A = sparse([1.0 1.0;1.0 -1.0]);
b = cumsum(ones(c_float(A.n)));

A = A::SparseMatrixCSC{c_float,c_int}

println("\n-------------------")
println("Testing: QDLDL_etree")
println("-------------------")
work   = zeros(c_int,A.n)
Lnz    = zeros(c_int,A.n)
etree  = zeros(c_int,A.n)

#Julia matrices are 1-index.  C code is 0 indexed.
#use only the upper triangle
Atriu = triu(A);
Ap = Atriu.colptr-1;
Ai = Atriu.rowval-1;
Ax = Atriu.nzval;

#Construct the elimination tree and column counts for A
sumLnz = ccall(QDLDL_etree,c_int,(c_int,Ptr{c_int},Ptr{c_int},Ptr{c_int},Ptr{c_int},Ptr{c_int}),A.n,Ap,Ai,work,Lnz,etree)

println("\n-------------------")
println("Testing: QDLDL_factor")
println("-------------------  ")

#allocate target matrix data for L
Ln = A.n;
Lp = zeros(c_int,A.n+1);
Li = zeros(c_int,sumLnz);
Lx = zeros(c_float,sumLnz);

#allocate D and its inverse
D      = zeros(c_float,A.n)
Dinv   = zeros(c_float,A.n)

#scratch workspace values
fwork            = zeros(c_float,Ln)
iwork            = zeros(c_int,Ln*3)
bwork            = zeros(c_bool,Ln)

#call the C factorization function
ccall(QDLDL_factor,c_int,
     (     c_int,              #n
           Ptr{c_int},         #Ap
           Ptr{c_int},         #Ai
           Ptr{c_float},       #Ax
           Ptr{c_int},         #Lp
           Ptr{c_int},         #Li
           Ptr{c_float},       #Lx
           Ptr{c_float},       #D
           Ptr{c_float},       #Dinv
           Ptr{c_int},         #Lnz
           Ptr{c_int},         #etree
           Ptr{c_bool},        #bwork
           Ptr{c_int},         #iwork
           Ptr{c_float}),      #fwork
      A.n,Ap,Ai,Ax,Lp,Li,Lx,D,Dinv,Lnz,etree,bwork,iwork,fwork)

#Julia indexing is 1 based
L = SparseMatrixCSC(Ln,Ln,Lp+1,Li+1,Lx)
println("L = ", L)
println("D = ", D)
println("Inf-norm factorisation tolerance : ", norm( (L+I)*diagm(D)*(L'+I) - A, Inf)  );

println("\n-------------------")
println("Testing: QDLDL_solve")
println("-------------------  ")

#call the C solve function
x = copy(b)
ccall(QDLDL_solve,Void,
     (     c_int,              #n
           Ptr{c_int},         #Lp
           Ptr{c_int},         #Li
           Ptr{c_float},       #Lx
           Ptr{c_float},       #Dinv
           Ptr{c_float}),      #x (starts as b)
        Ln,Lp,Li,Lx,Dinv,x)

println("A\b solve tolerance: ", norm(x-A\b, Inf)  );

Libdl.dlclose(ldllibH)

max_issues_repo_name: WatsonZhouAnda/Cheetah-Software

max_forks_repo_name: WatsonZhouAnda/Cheetah-Software


text:
# using Revise, FLPDevelopment, BrowseTables

# Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Test"
# Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/Dev_raw_data"
pathMac = "/Volumes/GoogleDrive/My Drive/Flipping/Datasets/Pharmacology/DarioDevelopment"
pathLin = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Pharmacology/DarioDevelopment"
if ispath(pathMac)
    path = pathMac
elseif ispath(pathLin)
    path = pathLin
else
    error("unknown path")
end
##
Age_d, Age_p, Age_b, Age_s = process_dataset(path)
######### Filtering Datas #############
for df in (Age_p, Age_b, Age_s)
    df[!,:Age] = [x in dario_youngs ? "Juveniles" : "Adults" for x in df.MouseID]
    df[!,:Sex] = [x in females ? "F" : "M" for x in df.MouseID]
    # df[!,:PreInterpoke] = [ismissing(x) ? 0.0 : x for x in df.PreInterpoke]
    transform!(df, :Age => categorical, renamecols=false)
    transform!(df, :Sex => categorical, renamecols=false)
    gd = groupby(df,:Session)
    transform!(gd, :Streak => maximum => :Performance)
end

levels!(Age_p.Age,["Adults", "Juveniles"])
levels!(Age_b.Age,["Adults", "Juveniles"])
levels!(Age_s.Age,["Adults", "Juveniles"])

Afreq = countmap(Age_s.AfterLast)
Aprob = Dict()
for (a,f) in Afreq
    Aprob[a] = round(f/nrow(Age_s),digits = 5)
end
Age_s[!,:IncorrectStart] = [!x for x in Age_s.CorrectStart]
Age_s[!,:IncorrectLeave] = [!x for x in Age_s.CorrectLeave]
Age_s[!,:P_AfterLast] = [Aprob[a] for a in Age_s.AfterLast]
gd = groupby(Age_s,:Session)
transform!(gd, :AfterLast => frequency)
transform!(gd, :AfterLast => (x -> x .<= quantile(x,0.95)) => :Limit)
transform!(gd, :Num_Rewards => cumsum => :Cum_Rewards)
Age_s[!,:RewRate] = Age_s.Cum_Rewards ./ Age_s.Stop
Age_s[!,:TrialRewRate] = Age_s.Cum_Rewards ./ Age_s.Streak
gdp = groupby(Age_p,[:Session, :Streak])
transform!(gdp,:PokeIn => (x -> x .- x[1]) => :In)
transform!(gdp,[:PokeIn, :PokeOut] => ((x,y) -> y .- x[1]) => :Out)
# for (df,name) in zip([Age_p, Age_b, Age_s],["Pokes", "Bouts", "Streaks"])
#     filename = joinpath(path, "Results_" * string(today()),"FullInfo" * name * ".csv")
#     CSV.write(filename,df)
# end


text:
################################################################################
#
#  NfMaxOrd/ResidueRing.jl : Quotients of maximal orders of number fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# (C) 2016 Tommy Hofmann
#
################################################################################

export NfMaxOrdQuoRing, NfMaxOrdQuoRingElem, quo, *, -, ==, deepcopy, divrem,
       gcd, inv, parent, show, divexact, isone, iszero, howell_form,
       strong_echelon_form, triangularize, det, euclid, xxgcd, annihilator

################################################################################
#
#  Assert
#
################################################################################

add_assert_scope(:NfMaxOrdQuoRing)

set_assert_level(:NfMaxOrdQuoRing, 0)

################################################################################
#
#  Field access
#
################################################################################

elem_type(::Type{NfMaxOrdQuoRing}) = NfMaxOrdQuoRingElem

elem_type(::NfMaxOrdQuoRing) = NfMaxOrdQuoRingElem

base_ring(Q::NfMaxOrdQuoRing) = Q.base_ring

ideal(Q::NfMaxOrdQuoRing) = Q.ideal

basis_mat(Q::NfMaxOrdQuoRing) = Q.basis_mat

parent(x::NfMaxOrdQuoRingElem) = x.parent

parent_type(::Type{NfMaxOrdQuoRingElem}) = NfMaxOrdQuoRing

################################################################################
#
#  Functions to allow polynomial and polynomial ring constructions
#
################################################################################

needs_parentheses(::NfMaxOrdQuoRingElem) = true

is_negative(::NfMaxOrdQuoRingElem) = false

Base.promote_rule{S <: Integer}(::Type{NfMaxOrdQuoRingElem},
                                ::Type{S}) = NfMaxOrdQuoRingElem

################################################################################
#
#  Copying
#
################################################################################

deepcopy(x::NfMaxOrdQuoRingElem) =
        NfMaxOrdQuoRingElem(parent(x), deepcopy(x.elem))

################################################################################
#
#  I/O
#
################################################################################

function show(io::IO, Q::NfMaxOrdQuoRing)
  print(io, "Quotient of $(Q.base_ring)")
end

function show(io::IO, x::NfMaxOrdQuoRingElem)
  print(io, "$(x.elem)")
end

################################################################################
#
#  Parent object overloading
#
################################################################################

function call(O::NfMaxOrdQuoRing, x::NfOrdElem)
  parent(x) != base_ring(O) && error("Cannot coerce element into the quotient ring")
  return NfMaxOrdQuoRingElem(O, x)
end

function call(Q::NfMaxOrdQuoRing, x::Integer)
  return Q(base_ring(Q)(x))
end

function call(Q::NfMaxOrdQuoRing, x::fmpz)
  return Q(base_ring(Q)(x))
end

################################################################################
#
#  Quotient function
#
################################################################################

function quo(O::NfMaxOrd, I::NfMaxOrdIdl)
  # We should check that I is not zero
  Q = NfMaxOrdQuoRing(O, I)
  f = NfMaxOrdQuoMap(O, Q)
  return Q, f
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")
  return parent(x)(x.elem + y.elem)
end

function -(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")
  return parent(x)(x.elem - y.elem)
end

function -(x::NfMaxOrdQuoRingElem)
  return parent(x)(-x.elem)
end

function *(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")
  return parent(x)(x.elem * y.elem)
end

function *(x::Integer, y::NfMaxOrdQuoRingElem)
  return parent(y)(x * y.elem)
end

*(x::NfMaxOrdQuoRingElem, y::Integer) = y*x

function *(x::fmpz, y::NfMaxOrdQuoRingElem)
  return parent(x)(x * y.elem)
end

*(x::NfMaxOrdQuoRingElem, y::fmpz) = y*x

################################################################################
#
#  Special elements
#
################################################################################

iszero(x::NfMaxOrdQuoRingElem) = iszero(x.elem)

is_zero(x::NfMaxOrdQuoRingElem) = iszero(x)

isone(x::NfMaxOrdQuoRingElem) = isone(x.elem)

function one(Q::NfMaxOrdQuoRing)
  return Q(one(Q.base_ring))
end

function zero(Q::NfMaxOrdQuoRing)
  return Q(zero(Q.base_ring))
end

################################################################################
#
#  Equality
#
################################################################################

==(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem) = x.elem == y.elem

################################################################################
#
#  Exact division
#
################################################################################

function divexact(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  b, z = isdivisible(x, y)
  @assert b
  return z
end

function isdivisible(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")

  if iszero(x)
    return true, zero(parent(x))
  end

  R = parent(x)
  d = degree(base_ring(R))
  # We cannot solve with non-square matrices.
  # Thus build the matrix
  # ( 1   x    0  )
  # ( 0  M_y   I  )
  # ( 0  M_I   0  ).
  # Compute the UPPER RIGHT HNF ->
  # ( 1   0   u )
  # ( *   *   * )
  # u will be the coefficient vector of the quotient

  V = R.tmp_div
  A = representation_mat(y.elem)
  B = basis_mat(parent(x))

  V[1, 1] = 1

  a = elem_in_basis(x.elem)

  for i in 1:d
    V[1, 1 + i] = a[i]
  end

  _copy_matrix_into_matrix(V, 2, 2, A)   # this really is a copy
  _copy_matrix_into_matrix(V, 2+d, 2, B) # this really is a copy

  for i in 1:d
    V[1 + i, d + 1 + i] = 1
  end

  hnf_modular_eldiv!(V, minimum(parent(x).ideal))

  if !iszero(submat(V, 1:1, 2:(d + 1)))
    ccall((:fmpz_mat_zero, :libflint), Void, (Ptr{fmpz_mat}, ), &V)
    return false, zero(parent(x))
  end

  z = R(-base_ring(R)(fmpz[ deepcopy(V[1, i]) for i in (d + 2):(2*d + 1)]))

  ccall((:fmpz_mat_zero, :libflint), Void, (Ptr{fmpz_mat}, ), &V)

  @hassert :NfMaxOrdQuoRing 1 z*y == x
  return true, z
end

################################################################################
#
#  Strong exact division
#
################################################################################

function _divexact_strong(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  n = euclid(x)
  m = euclid(y)
  @hassert :NfMaxOrdQuoRing 1 mod(n, m) == 0
  target = divexact(n, m)

  #println("target valuation: $target")
  #println("doing first an ordinary divexact with\n $x \n and \n $y")
  q0 = divexact(x, y)
  #println("valuation of first quotient: $(euclid(q0))")

  if euclid(q0) == target
    return q0
  else
    i = annihilator(y)
    #println("generator of annihilator: $i")

    q = q0 + rand(parent(x))*i

    k = 0
    while euclid(q) != target
      k += 1
      q = q0 + rand(parent(x))*i
      #println("current valuation $(euclid(q))")
      if k > 500
        error("Could not find proper quotion for strong division")
      end
    end
  end

  @hassert :NfMaxOrdQuoRing 1 q*y == x
  @hassert :NfMaxOrdQuoRing 1 euclid(q) *euclid(y) == euclid(x)

  return q
end

################################################################################
#
#  Inverse element
#
################################################################################

function inv(x::NfMaxOrdQuoRingElem)
  return divexact(one(parent(x)), x)
end

################################################################################
#
#  Euclidean function
#
################################################################################

function euclid(x::NfMaxOrdQuoRingElem)
  if is_zero(x)
    return fmpz(-1)
  end

  U = parent(x).tmp_euc

  d = degree(base_ring(parent(x)))

  _copy_matrix_into_matrix(U, 1, 1, representation_mat(x.elem))
  _copy_matrix_into_matrix(U, d + 1, 1, parent(x).basis_mat)

  hnf_modular_eldiv!(U, parent(x).ideal.minimum)

  z = fmpz(1)

  for i in 1:degree(base_ring(parent(x)))
    mul!(z, z, U[i, i])
  end

  @hassert :NfMaxOrdQuoRing 1 z == norm(ideal(parent(x.elem), x.elem) + parent(x).ideal)

  return z
end

################################################################################
#
#  Division with remainder
#
################################################################################

function divrem(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  q = rand(parent(x))
  r = x - q*y
  e = euclid(y)

  # This should be only one case and don't do the try/catch crap
  # Write a proper _is_divisible function

  if e == 1
    q = x*inv(y)
    r = x - q*y
    @hassert :NfMaxOrdQuoRing 1 iszero(x - q*y)
    @hassert :NfMaxOrdQuoRing 1 euclid(r) < e
    return q, r
  end

  try q = divexact(x, y)
    r = x - q*y
    @hassert :NfMaxOrdQuoRing 1 iszero(x - q*y)
    @hassert :NfMaxOrdQuoRing 1 euclid(r) < e
    return q, r
  catch
  end

  cnt = 0
  while euclid(r) >= e
    cnt += 1
    q = rand(parent(x))
    r = x - q*y
    if cnt > 1000
      println("Target valuation $e")
      error("Something odd in divrem for $x $y $(parent(x))")
    end
  end

  @hassert :NfMaxOrdQuoRing 1 euclid(r) < e

  return q, r
end

################################################################################
#
#  Random elements
#
################################################################################

function rand(Q::NfMaxOrdQuoRing)
  A = basis_mat(Q)
  B = basis(base_ring(Q))
  z = rand(fmpz(1):A[1,1]) * B[1]

  for i in 2:rows(A)
    z = z + rand(fmpz(1):A[i, i]) * B[i]
  end

  return Q(z)
end

################################################################################
#
#  Annihilator
#
################################################################################

function annihilator(x::NfMaxOrdQuoRingElem)
  I = parent(x).ideal
  O = base_ring(parent(x))
  d = degree(O)
  f = NfMaxOrdQuoMap(O, parent(x))
  U = parent(x).tmp_ann
  # We first compute a Z-basis of the annihilator ideal
  # The basis is the kernel of the following matrix
  # ( M_I )
  # ( M_x )
   _copy_matrix_into_matrix(U, 1, 1, representation_mat(x.elem))
   _copy_matrix_into_matrix(U, d + 1, 1, I.basis_mat)

  m = _kernel(U)
  I = ideal(O, _hnf_modular_eldiv(submat(m, 1:degree(O), 1:degree(O)),
                                  minimum(I), :lowerleft))
  z = f(I)
  @assert iszero(z*x)
  return z
end

################################################################################
#
#  GCD
#
################################################################################

function gcd(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  Q = parent(x)
  O = base_ring(Q)

  I = ideal(O, x.elem) + ideal(O, y.elem)

  f = NfMaxOrdQuoMap(O, parent(x))

  return f(I)
end

################################################################################
#
#  Extended extended GCD
#
################################################################################

function xxgcd(x::NfMaxOrdQuoRingElem, y::NfMaxOrdQuoRingElem)
  Q = parent(x)
  O = base_ring(Q)

  d = degree(O)

  if is_zero(x)
    return deepcopy(y), Q(O(0)), Q(O(1)), Q(O(-1)), Q(O(0))
  elseif is_zero(y)
    return deepcopy(x), Q(O(1)), Q(O(0)), Q(O(0)), Q(O(1))
  end

  g = gcd(x, y)

  e = _divexact_strong(x, g)
  f = _divexact_strong(y, g)

  M_e = representation_mat(e.elem)
  M_f = representation_mat(f.elem)

  M_I = basis_mat(Q)

  # let us build
  # ( 1  Q(1) 0  0 )
  # ( 0  M_e  I  0 )
  # ( 0  M_f  0  I )
  # ( 0  M_I  0  0 )

  a = elem_in_basis(Q(O(1)).elem)

  V = parent(x).tmp_xxgcd

  V[1, 1] = 1

  for i in 1:d
    V[1, 1 + i] = a[i]
  end

  Hecke._copy_matrix_into_matrix(V, 2, 2, M_e)
  Hecke._copy_matrix_into_matrix(V, d + 2, 2, M_f)
  Hecke._copy_matrix_into_matrix(V, 2*d + 2, 2, M_I)

  for i in 1:2*d
    V[1+i, 1 + d + i] = 1
  end

  U = V

  U = hnf_modular_eldiv(U, minimum(Q.ideal))::fmpz_mat

  u = Q(-O([ U[1,i] for i in (d + 2):(2*d + 1)]))
  v = Q(-O([ U[1,i] for i in (2*d + 2):(3*d + 1)]))

  @hassert :NfMaxOrdQuoRing 1 Q(O(1)) == u*e - (v*(-f))

  ccall((:fmpz_mat_zero, :libflint), Void, (Ptr{fmpz_mat}, ), &V)

  return g, u, v, -f, e
end

################################################################################
#
#  Triangularization
#
################################################################################

function _pivot(A, start_row, col)
  if !is_zero(A[start_row, col])
    return 1;
  end

  for j in start_row + 1:rows(A)
    if !is_zero(A[j, col])
      swap_rows!(A, j, start_row)
      return -1
    end
  end

  return 0
end

function _strong_echelon_form(A::GenMat{NfMaxOrdQuoRingElem})
  B = deepcopy(A)

  if rows(B) < cols(B)
    B = vcat(B, MatrixSpace(base_ring(B), cols(B) - rows(B), cols(B))())
  end

  strong_echelon_form!(B)
  return B
end

function strong_echelon_form(A::GenMat{NfMaxOrdQuoRingElem}, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = _strong_echelon_form(_swapcols(A))
    _swapcols!(h)
    _swaprows!(h)
    return h
  elseif shape == :upperright
    return _strong_echelon_form(A)
  else
    error("Not yet implemented")
  end
end

function triangularize!(A::GenMat{NfMaxOrdQuoRingElem})
  n = rows(A)
  m = cols(A)
  d = one(base_ring(A))

  t_isdiv = 0.0
  t_xxgcd = 0.0
  t_arith = 0.0

  row = 1
  col = 1
  tic()
  while row <= rows(A) && col <= cols(A)
    #println("doing row $row")
    t = _pivot(A, row, col)
    if t == 0
      col = col + 1
      continue
    end
    d = d*t
    for i in (row + 1):rows(A)
      if iszero(A[i, col])
        continue
      end

      t_isdiv += @elapsed b, q = isdivisible(A[i, col], A[row, col])

      if b
        for k in col:m
          t_arith += @elapsed A[i, k] = A[i, k] - q*A[row, k]
        end
        @hassert :NfMaxOrdQuoRing 1 A[i, col] == zero(base_ring(A))
      else
        t_xxgcd += @elapsed g,s,t,u,v = xxgcd(A[row, col], A[i, col])
        @hassert :NfMaxOrdQuoRing 1 isone(s*v - t*u)

        for k in col:m
          t_arith += @elapsed t1 = s*A[row, k] + t*A[i, k]
          t_arith += @elapsed t2 = u*A[row, k] + v*A[i, k]
          A[row, k] = t1
          A[i, k] = t2
        end
      end
    end
    row = row + 1;
    col = col + 1;
  end
  #println("  === Time triangularization")
  #println("    isdivisbible: $t_isdiv")
  #println("    xxgcd       : $t_xxgcd")
  #println("    arith       : $t_arith")
  #println("    total time  : $(toc())")
  return d
end

function triangularize(A::GenMat{NfMaxOrdQuoRingElem})
  #println("copying ...")
  B = deepcopy(A)
  #println("done")
  triangularize!(B)
  return B
end

################################################################################
#
#  Strong echelon form
#
################################################################################

function strong_echelon_form!(A::GenMat{NfMaxOrdQuoRingElem})
  #A = deepcopy(B)
  n = rows(A)
  m = cols(A)

  @assert n >= m

  #print("triangularizing ... ")
  triangularize!(A)

  T = MatrixSpace(base_ring(A), 1, cols(A))()

  # We do not normalize!
  for j in 1:m
    if !is_zero(A[j,j]) != 0
      # This is the reduction
      for i in 1:j-1
        if iszero(A[i, j])
          continue
        else
          q, r = divrem(A[i, j], A[j, j])
          for l in i:m
            A[i, l] = A[i, l] - q*A[j, l]
          end
        end
      end

      a = annihilator(A[j, j])

      for k in 1:m
        T[1, k] = a*A[j, k]
      end
    else
      for k in 1:m
        T[1, k] = A[j, k]
      end
    end

    for i in j+1:m

      if iszero(T[1, i])
        continue
      end

      b, q = isdivisible(T[1, i], A[i, i])

      if b
        for k in i:m
          T[1, k] = T[1, k] - q*A[i, k]
        end
        @hassert :NfMaxOrdQuoRing 1 T[1, i] == zero(base_ring(A))
      else
        g,s,t,u,v = xxgcd(A[i, i], T[1, i])

        for k in i:m
          t1 = s*A[i, k] + t*T[1, k]
          t2 = u*A[i, k] + v*T[1, k]
          A[i, k] = t1
          T[1, k] = t2
        end
      end
    end
  end
  return A
end

################################################################################
#
#  Howell form
#
################################################################################

function howell_form!(A::GenMat{NfMaxOrdQuoRingElem})
  @assert rows(A) >= cols(A)

  k = rows(A)

  strong_echelon_form!(A)

  for i in 1:rows(A)
    if is_zero_row(A, i)
      k = k - 1

      for j in (i + 1):rows(A)
        if !is_zero_row(A, j)
          swap_rows!(A, i, j)
          j = rows(A)
          k = k + 1
        end
      end
    end
  end
  return k
end

function howell_form(A::GenMat{NfMaxOrdQuoRingElem})
  B = deepcopy(A)

  if rows(B) < cols(B)
    B = vcat(B, MatrixSpace(base_ring(B), cols(B) - rows(B), cols(B))())
  end

  howell_form!(B)

  return B
end

################################################################################
#
#  Determinant
#
################################################################################

function det(M::GenMat{NfMaxOrdQuoRingElem})
  rows(M) != cols(M) && error("Matrix must be square matrix")
  N = deepcopy(M)
  d = triangularize!(N)
  z = one(base_ring(M))
  for i in 1:rows(N)
    z = z * N[i, i]
  end
  return z*d
  q, r = divrem(z, d)
  @hassert :NfMaxOrdQuoRing 1 iszero(r)
  return divexact(z, d)
end

################################################################################
#
#  Functions for matrix spaces
#
################################################################################

#function call(M::GenMatSpace{NfMaxOrdQuoRingElem}, x::GenMat{NfOrdElem{NfMaxOrd}})
#  base_ring(base_ring(M)) != base_ring(parent(x)) &&
#      error("Base rings do not coincide")
#  z = M()
#  R = base_ring(M)
#  for i in 1:rows(x)
#    for j in 1:cols(x)
#      z[i, j] = R(x[i, j])
#    end
#  end
#  return z
#end

function Base.call(M::GenMatSpace{NfMaxOrdQuoRingElem}, x::GenMat{NfOrdElem{NfMaxOrd}})
  z = map(base_ring(M), x.entries)::Array{NfMaxOrdQuoRingElem, 2}
  return M(z)::GenMat{NfMaxOrdQuoRingElem}
end
################################################################################
#
#  Hensel Lifting
#
################################################################################

## Hensel lifting of roots
# This will fail for too large input
# Need to incorporate the explicit lifting bounds
function _roots_hensel{T}(f::GenPoly{NfOrdElem{T}})
  O = base_ring(f)
  n = degree(O)
  deg = degree(f)

  # First we find a prime ideal such that f is squarefree modulo P
  # (The discriminant of f has only finitely many divisors).

  p = 3

  fder = derivative(f)

  found_prime = false

  # Dummy variable
  Q = NfMaxOrdIdl(O)
  pi_F = NfMaxOrdToFqNmodMor()
  lin_factor = Array{fq_nmod_poly, 1}()

  while !found_prime
    p = next_prime(p)

    if is_index_divisor(O, p) || isramified(O, p)
      continue
    end

    lP = prime_decomposition(O, p)

    for P in lP
      F, pi_F = ResidueField(O, P[1])

      fmodP = pi_F(f)

      if !issquarefree(fmodP)
        continue
      end

      fac = factor(fmodP)

      for i in keys(fac)
        if degree(i) == 1
          push!(lin_factor, i)
        end
      end

      Q = P[1]
      found_prime = true
      break
    end
  end

  fmodQ = pi_F(f)

  # compute the lifting exponent a la Friedrich-Fieker

  R = ArbField(64)
  z = R(0)
  (c1, c2) = norm_change_const(O)
  (r1, r2) = signature(O)

  cc2 = (exp(3*log(R(3))//2) * R(3)^deg)//(const_pi(R) * R(deg))

  bound_root = [ R(0) for i in 1:(r1 + r2) ]

  for j in 0:deg-1
    co = coeff(f, j)
    co_conj = conjugates_arb(co)
    for i in 1:r1+r2
      bound_root[i] += inv(binom(R(deg), UInt(j))) * abs(co_conj[i])^2
    end
  end

  for i in 1:r1
    bound_root[i] = bound_root[i] * cc2
  end

  for i in r1+1:r2
    bound_root[i] = 2 * bound_root[i] * cc2
  end

  boundt2 = sum(bound_root)

  boundk = R(n)*log(R(c1)*R(c2)*boundt2*exp((R(n*(n-1))//4 + 2)*log(R(2))))//(2*Q.splitting_type[2]*log(R(p)))

  t = arf_struct(0, 0, 0, 0)
  ccall((:arf_init, :libarb), Void, (Ptr{arf_struct}, ), &t)

  ss = fmpz()
  ccall((:arb_get_abs_ubound_arf, :libarb), Void, (Ptr{arf_struct}, Ptr{arb}, Clong), &t, &boundk, 64)
  ccall((:arf_get_fmpz, :libarb), Void, (Ptr{fmpz}, Ptr{arf_struct}, Cint), &ss, &t, 1) # 1 is round up
  ccall((:arf_clear, :libarb), Void, (Ptr{arf_struct}, ), &t)

  roots = NfOrdElem[]

  s = Int(ss)

  for j in 1:length(lin_factor)

    zero_mod_Q = - coeff(lin_factor[j], 0)

    # The following should be a uniformizing element
    Q_pi = Q.gen_two

    @hassert :NfMaxOrdQuoRing 1 fmodQ(zero_mod_Q) == 0

    # This is the first step

    R, pi_R = quo(O, Q^2)

    t1 = divexact(pi_R(subst(f, pi_F\zero_mod_Q)), pi_R(Q_pi))
    t2 = -inv(pi_R(subst(fder,pi_F\zero_mod_Q)))
    new_a = pi_R\(pi_R(pi_F\zero_mod_Q) + t1*t2*pi_R(Q_pi))
    #return pi_R(f)

    old_a = new_a

    RR, pi_RR = R, pi_R

    I = Q^2

    reconstructed_new = old_a
    reconstructed_old = reconstructed_new

    stabilized = -1

    for i in 2:s
      if reconstructed_new == reconstructed_old
        stabilized = stabilized + 1
      else
        stabilized = 0
      end

      if stabilized >= 3
        if subst(f, reconstructed_new) == 0
          push!(roots, reconstructed_new)
          break
        else
          stabilized = 0
        end
      end

      if i == s
        if iszero(subst(f, reconstructed_new))
          push!(roots, reconstructed_new)
        end
        break
      end

      reconstructed_old = reconstructed_new
      old_a = new_a
      R, pi_R = RR, pi_RR
      I = I * Q

      # From Q^i -> Q^(i+1)

      RR, pi_RR = quo(O, I)
      t1 = divexact(pi_RR(subst(f, old_a)), pi_RR(Q_pi)^(i))
      t2 = pi_RR(pi_F\(-inv(pi_F(subst(fder, old_a)))))
      new_a = pi_RR\(pi_RR(old_a) + t1*t2*pi_RR((Q_pi))^(i))

      # Try to reconstruct:

      B = basis_mat(I)
      L = lll(B)

      rhs = MatrixSpace(ZZ, degree(O), 1)(elem_in_basis(new_a)'')
      lhs = transpose(L)

      X, d = solve(lhs, rhs)

      zz = [ fmpq(BigInt(X[i, 1])//BigInt(d) - round(BigInt(X[i, 1])//BigInt(d))) for i in 1:degree(O)]

      cden = den(zz[1])

      for i in 2:degree(O)
        cden = lcm(cden, den(zz[i]))
      end

      zz_num = [ num(cden*zz[i]) for i in 1:degree(O) ]

      v = MatrixSpace(FlintZZ, 1, degree(O))(zz_num')

      w = v*L

      # There is no slower function

      reconstructed_new = O(fmpz[ divexact(w[1, i], cden) for i in 1:degree(O) ])
    end
  end
  return roots
end

function probablity(O::NfMaxOrdQuoRing)
  p = 1.0
  I = O.ideal
  f = factor(norm(I))
  for (P, e) in f
    if valuation(I, P) > 0
      p = p * (1 - 1/Float64(norm(P)))
    end
  end
  return p
end


text:
function Base.show(io::IO, cfg::CFG)
    foreach(pairs(cfg.blocks)) do (idx, block)
        println("$idx\t=>\t", join(block.succs, ", "))
    end
end

print_ssa(io::IO, val) = isa(val, SSAValue) ? print(io, "%$(val.id)") : print(io, val)
function print_node(io::IO, idx, stmt, used, maxsize; color = true, print_typ=true)
    if idx in used
        pad = " "^(maxsize-length(string(idx)))
        print(io, "%$idx $pad= ")
    else
        print(io, " "^(maxsize+4))
    end
    if isa(stmt, PhiNode)
        args = map(1:length(stmt.edges)) do i
            e = stmt.edges[i]
            v = !isassigned(stmt.values, i) ? "#undef" :
                sprint() do io′
                    print_ssa(io′, stmt.values[i])
                end
            "$e => $v"
        end
        print(io, "φ ", '(', join(args, ", "), ')')
    elseif isa(stmt, PiNode)
        print(io, "π (")
        print_ssa(io, stmt.val)
        print(io, ", ")
        if color
            printstyled(io, stmt.typ, color=:red)
        else
            print(io, stmt.typ)
        end
        print(io, ")")
    elseif isa(stmt, ReturnNode)
        if !isdefined(stmt, :val)
            print(io, "unreachable")
        else
            print(io, "return ")
            print_ssa(io, stmt.val)
        end
    elseif isa(stmt, GotoIfNot)
        print(io, "goto ", stmt.dest, " if not ")
        print_ssa(io, stmt.cond)
    elseif isexpr(stmt, :call)
        print_ssa(io, stmt.args[1])
        print(io, "(")
        print(io, join(map(arg->sprint(io->print_ssa(io, arg)), stmt.args[2:end]), ", "))
        print(io, ")")
        if print_typ && stmt.typ !== Any
            print(io, "::$(stmt.typ)")
        end
    elseif isexpr(stmt, :new)
        print(io, "new(")
        print(io, join(map(arg->sprint(io->print_ssa(io, arg)), stmt.args), ", "))
        print(io, ")")
    else
        print(io, stmt)
    end
end

function Base.show(io::IO, code::IRCode)
    io = IOContext(io, :color=>true)
    used = Set{Int}()
    println(io, "Code")
    foreach(stmt->scan_ssa_use!(used, stmt), code.stmts)
    foreach(((_a, _b, node, _d),) -> scan_ssa_use!(used, node), code.new_nodes)
    if isempty(used)
        maxsize = 0
    else
        maxused = maximum(used)
        maxsize = length(string(maxused))
    end
    cfg = code.cfg
    max_bb_idx_size = length(string(length(cfg.blocks)))
    bb_idx = 1
    perm = sortperm(code.new_nodes, by = x->x[1])
    new_nodes_perm = Iterators.Stateful(perm)
    for (idx, stmt) in Iterators.enumerate(code.stmts)
        bbrange = cfg.blocks[bb_idx].stmts
        bbrange = bbrange.first:bbrange.last
        bb_pad = max_bb_idx_size - length(string(bb_idx))
        if idx != last(bbrange)
            if idx == first(bbrange)
                print(io, "$(bb_idx) ","─"^(1+bb_pad)," ")
            else
                print(io, "│  "," "^max_bb_idx_size)
            end
        end
        print_sep = false
        if idx == last(bbrange)
            print_sep = true
        end
        floop = true
        while !isempty(new_nodes_perm) && code.new_nodes[Base.peek(new_nodes_perm)][1] == idx
            node_idx = popfirst!(new_nodes_perm)
            _, typ, node, line = code.new_nodes[node_idx]
            node_idx += length(code.stmts)
            if print_sep
                if floop
                    print(io, "$(bb_idx) ","─"^(1+bb_pad)," ")
                else
                    print(io, "│  "," "^max_bb_idx_size)
                end
            end
            print_sep = true
            floop = false
            print_ssa_typ = !isa(node, PiNode) && node_idx in used
            Base.with_output_color(:yellow, io) do io′
                print_node(io′, node_idx, node, used, maxsize; color = false,
                    print_typ=!print_ssa_typ || (isa(node, Expr) && typ != node.typ))
            end
            if print_ssa_typ
                printstyled(io, "::$(typ)", color=:red)
            end
            println(io)
        end
        if print_sep
            if idx == first(bbrange) && floop
                print(io, "$(bb_idx) ","─"^(1+bb_pad)," ")
            else
                print(io, idx == last(bbrange) ? string("└", "─"^(1+max_bb_idx_size), " ") :
                    string("│  ", " "^max_bb_idx_size))
            end
        end
        if idx == last(bbrange)
            bb_idx += 1
        end
        typ = code.types[idx]
        print_ssa_typ = !isa(stmt, PiNode) && idx in used
        print_node(io, idx, stmt, used, maxsize,
            print_typ=!print_ssa_typ || (isa(stmt, Expr) && typ != stmt.typ))
        if print_ssa_typ
            printstyled(io, "::$(typ)", color=:red)
        end
        println(io)
    end
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
using Plots, DelimitedFiles, CSV, DataFrames

mutable struct netAbs
    nodes
    adj
    ep
end

function loadData(;dir="networks")
    files = readdir(pwd() * "/$(dir)/")
    nodesFiles = files[findall(x->endswith(x,"nodes.csv"), files)]
    adjFiles = files[findall(x->endswith(x,"adj.csv"), files)]
    epFiles = files[findall(x->endswith(x,"endpoints.csv"), files)]
    net = Array{netAbs}(undef, size(nodesFiles,1))
    if size(nodesFiles,1) == size(adjFiles,1) == size(epFiles,1)
        for (i,node) in enumerate(nodesFiles)
            # @info i, node, adj, ep
            loadedNodes = Array(CSV.read(pwd() * "/networks/$(node)", DataFrame, header=0 ))[:,1:2]
            loadedAdj = Array(CSV.read(pwd() * "/networks/$(split(node, "_nodes")[1] * "_adj.csv")", DataFrame, header=0 ))
            loadedEndPoints = Array(CSV.read(pwd() * "/networks/$(split(node, "_nodes")[1] * "_endpoints.csv")", DataFrame, header=0 ))[:,1]
            net[i] = netAbs(loadedNodes, loadedAdj, loadedEndPoints)
        end
    else
        @warn "Not same number of nodes, adj and endpoint files."
    end
    return net
end

function plotNetwork(idx)
    fig = Plots.hline(0)
    for i=1:size(net[idx].nodes,1), j=1:size(net[idx].nodes,1)
        if net[idx].adj[j,i]==1
            fig = plot!([net[idx].nodes[j,2],net[idx].nodes[i,2]],[net[idx].nodes[j,1],net[idx].nodes[i,1]],
              markershapes = :circle, lw=2, ms=5, legend=:none, c=:black, markerstrokewidth=0)
        end
    end
    fig = Plots.hline!(0)
    return fig
end

cd("/Users/javier/Desktop/SyntheticVasculature/Data")
net = loadData()

using Interact

@manipulate for i in 1:size(net,1)
    f1 = plot(reshape(sum(net[1].adj, dims=1),:), seriestype=:barhist, width=0, normalize=true, tickfont = font(10, "Helvetica"), frame=:box, size=(700,700), legend=false, xlabel="Node degree", ylabel="PDF")
    f2 = plotNetwork(i)
    plot(f1,f2, size=(2500,2000), layout = (2, 1))
end

# plot(reshape(sum(net[1].adj, dims=1),:), seriestype=:barhist, width=0, normalize=true, tickfont = font(10, "Helvetica"), frame=:box, size=(700,700), legend=false, xlabel="Node degree", ylabel="PDF")

max_stars_repo_name: jquetzalcoatl/SyntheticVasculature

max_issues_repo_name: jquetzalcoatl/SyntheticVasculature

max_forks_repo_name: jquetzalcoatl/SyntheticVasculature


text:
using Documenter, PlaneGeometry, CodeTracking

DocMeta.setdocmeta!(PlaneGeometry, :DocTestSetup, :(using PlaneGeometry); recursive=true)

makedocs(sitename="Plane Geometry in Julia",
         modules = [PlaneGeometry],
         pages = [
                  "index.md",
                  "Geometric Definitions" =>
                  [
                   "Summary" => "definition/index.md",
                   "Elementary" => "definition/elementary.md",
                   "Distance" => "definition/distance.md",
                   "Triangles" => "definition/triangle.md",
                   "Circles" => "definition/circle.md",
                  ],
                  "Theorems" =>
                  [
                   "Summary" => "theorem/index.md",
                   "Napoleon" => "theorem/Napoleon.md",
                   "Centroid" => "theorem/Centroid.md",
                   "Ceva" => "theorem/Ceva.md",
                  ],
                  "Utitlies" => "util.md",
                  "Ackowledgement" => "thanks.md",
                 ],
         format = Documenter.HTML(assets = ["assets/custom.css"]),
         doctest = true
        )


text:
using DimensionalData, Test, Unitful
using DimensionalData: X, Y, Z, Time, Forward, Reverse, Ordered, arrayorder, indexorder, relationorder

a = [1 2  3  4
     5 6  7  8
     9 10 11 12]

@testset "Selectors on IndependentGrid" begin
    da = DimensionalArray(a, (Y((10, 30)), Time((1:4)u"s")))
    dims(da)

    @test At(10.0) == At(10.0, 0.0, Base.rtoldefault(eltype(10.0)))
    x = [10.0, 20.0]
    @test At(x) === At(x, 0.0, Base.rtoldefault(eltype(10.0)))
    @test At((10.0, 20.0)) === At((10.0, 20.0), 0.0, Base.rtoldefault(eltype(10.0)))

    Near([10, 20])

    @test Between(10, 20) == Between((10, 20))

    @testset "selectors with dim wrappers" begin
        @test da[Y<|At([10, 30]), Time<|At([1u"s", 4u"s"])] == [1 4; 9 12]
        @test_throws ArgumentError da[Y<|At([9, 30]), Time<|At([1u"s", 4u"s"])]
        @test view(da, Y<|At(20), Time<|At((3:4)u"s")) == [7, 8]
        @test view(da, Y<|Near(17), Time<|Near([1.5u"s", 3.1u"s"])) == [5, 7]
        @test view(da, Y<|Between(9, 21), Time<|At((3:4)u"s")) == [3 4; 7 8]
    end

    @testset "selectors without dim wrappers" begin
        @test da[At(20:10:30), At(1u"s")] == [5, 9]
        @test view(da, Between(9, 31), Near((3:4)u"s")) == [3 4; 7 8; 11 12]
        @test view(da, Near(22), At([3.0u"s", 4.0u"s"])) == [7, 8]
        @test view(da, At(20), At((2:3)u"s")) == [6, 7]
        @test view(da, Near<|13, Near<|[1.3u"s", 3.3u"s"]) == [1, 3]
        # Near works with a tuple input
        @test view(da, Near<|(13,), Near<|[1.3u"s", 3.3u"s"]) == [1 3]
        @test view(da, Between(11, 20), At((2:3)u"s")) == [6 7]
        # Between also accepts a tuple input
        @test view(da, Between((11, 20)), Between((2u"s", 3u"s"))) == [6 7]
    end

    @testset "more Unitful dims" begin
        dimz = Time<|1.0u"s":1.0u"s":3.0u"s", Y<|(1u"km", 4u"km")
        db = DimensionalArray(a, dimz)
        @test db[Y<|Between(2u"km", 3.9u"km"), Time<|At<|3.0u"s"] == [10, 11]
    end

    @testset "selectors work in reverse orders" begin
        @testset "forward index with reverse relation" begin
            da_ffr = DimensionalArray(a, (Y(10:10:30; grid=RegularGrid(order=Ordered(Forward(), Forward(), Reverse()))),
                                         Time((1:1:4)u"s"; grid=RegularGrid(order=Ordered(Forward(), Forward(), Reverse())))))
            @test indexorder(dims(da_ffr, Time)) == Forward()
            @test arrayorder(dims(da_ffr, Time)) == Forward()
            @test relationorder(dims(da_ffr, Time)) == Reverse()
            @test da_ffr[Y<|At([10, 30]), Time<|At([1u"s", 4u"s"])] == [12 9; 4 1]
            @test da_ffr[Y<|At(20), Time<|At((3.0:4.0)u"s")] == [6, 5]
            @test da_ffr[Y<|Near(7), Time<|Near([1.3u"s", 3.3u"s"])] == [12, 10]
            @test da_ffr[Y<|Between(9, 21), Time<|At((3.0:4.0)u"s")] == [10 9; 6 5]
        end

        @testset "is the same as reverse index with forward realation" begin
            da_rff = DimensionalArray(a, (Y(30:-10:10; grid=RegularGrid(order=Ordered(Reverse(), Forward(), Forward()))),
                                         Time((4:-1:1)u"s"; grid=RegularGrid(order=Ordered(Reverse(), Forward(), Forward())))))
            @test da_rff[Y<|At([10, 30]), Time<|At([1u"s", 4u"s"])] == [12 9; 4 1]
            @test da_rff[Y<|At(20), Time<|At((3.0:4.0)u"s")] == [6, 5]
            @test da_rff[Y<|Near(7), Time<|Near([1.3u"s", 3.3u"s"])] == [12, 10]
            @test da_rff[Y<|Between(9, 21), Time<|At((3.0:4.0)u"s")] == [10 9; 6 5]
        end

    end

    @testset "setindex! with selectors" begin
        c = deepcopy(a)
        dc = DimensionalArray(c, (Y((10, 30)), Time((1:4)u"s")))
        dc[Near(11), At(3u"s")] = 100
        @test c[1, 3] == 100
        dc[Time<|Near(2.2u"s"), Y<|Between(10, 30)] = [200, 201, 202]
        @test c[1:3, 2] == [200, 201, 202]
    end

end

@testset "CategoricalGrid" begin
    dimz = Time([:one, :two, :three]; grid=CategoricalGrid()),
           Y([:a, :b, :c, :d]; grid=CategoricalGrid())
    da = DimensionalArray(a, dimz)
    @test da[Time<|At([:one, :two]), Y<|At(:b)] == [2, 6]
    @test da[At([:one, :three]), At([:b, :c, :d])] == [2 3 4; 10 11 12]
    @test da[At(:two), Between(:b, :d)] == [6, 7, 8]
    # Near doesn't make sense for categories
    @test_throws ArgumentError da[Near(:two), At([:b, :c, :d])]
end

@testset "TranformedGrid " begin
    using CoordinateTransformations

    m = LinearMap([0.5 0.0; 0.0 0.5])

    dimz = Dim{:trans1}(m; grid=TransformedGrid(X())),
           Dim{:trans2}(m, grid=TransformedGrid(Y()))

    @testset "permutedims works on grid dimensions" begin
        @test permutedims((Y(), X()), dimz) == (X(), Y())
    end

    da = DimensionalArray(a, dimz)

    @testset "Indexing with array dims indexes the array as usual" begin
        @test da[Dim{:trans1}(3), Dim{:trans2}(1)] == 9
        # Using selectors works the same as indexing with grid
        # dims - it applies the transform function.
        # It's not clear this should be allowed or makes sense,
        # but it works anyway because the permutation is correct either way.
        @test da[Dim{:trans1}(At(6)), Dim{:trans2}(At(2))] == 9
    end

    @testset "Indexing with grid dims uses the transformation" begin
        @test da[X(Near(6.1)), Y(Near(8.5))] == 12
        @test da[X(At(4.0)), Y(At(2.0))] == 5
        @test_throws InexactError da[X(At(6.1)), Y(At(8))]
        # Indexing directly with grid dims also just works, but maybe shouldn't?
        @test da[X(2), Y(2)] == 6
    end
end


text:
using LinearAlgebra
using Test

#include("test_solve_eye_minus_a_kron_b.jl")
include("test_solve_eye_plus_at_kron_b.jl")