text:
using Revise # reduce need for recompile
using Plots
using Documenter
using LinearAlgebra
using SparseArrays

# "User defined modules"
push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG
using UnPack

"Approximation parameters"
N   = 3  # order of approximation
K1D = 16 # number of elements per side
CFL = .75
T   = .75

"======= initialize reference element data and mesh ======="

# specify lobatto nodes to automatically get DG-SEM mass lumping
rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))

VX, VY, EToV = uniform_quad_mesh(K1D, K1D)
md = init_mesh((VX,VY),EToV,rd)

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

" ====== set up initial conditions ====== "

@unpack x,y = md
p = @. exp(-100*(x^2+(y-.25)^2))
u = zeros(size(x))
v = zeros(size(x))

"Time integration"
rk4a, rk4b, rk4c = rk45_coeffs()
CN = (N+1)*(N+2)  # estimated trace constant
dt = CFL * 2 / (CN*K1D)
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

"varying wavespeed"
c2 = @. 1 + .5*sin(pi*x)*sin(pi*y)
# c2 = ones(size(x))

function rhs(Q, rd::RefElemData, md::MeshData, params...)

    (p,u,v) = Q
    @unpack Dr,Ds,LIFT,Vf = rd
    @unpack rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ = md
    @unpack mapP,mapB = md

    QM = (x->Vf*x).(Q)
    QP = (x->x[mapP]).(QM)
    (dp,du,dv) = QP.-QM

    "compute central numerical flux"
    # (pM,uM,vM) = QM
    # du[mapB] = -2*uM[mapB]
    # dv[mapB] = -2*vM[mapB]
    pflux = @. du*nxJ + dv*nyJ
    uflux = @. dp*nxJ
    vflux = @. dp*nyJ

    pr,ur,vr = (x->Dr*x).(Q)
    ps,us,vs = (x->Ds*x).(Q)

    px = @. rxJ*pr + sxJ*ps;
    py = @. ryJ*pr + syJ*ps
    ux = @. rxJ*ur + sxJ*us;
    vy = @. ryJ*vr + syJ*vs

    rhsp = (ux+vy) + .5*LIFT*pflux
    rhsu = px + .5*LIFT*uflux
    rhsv = py + .5*LIFT*vflux

    c2 = params[1]
    rhsp = @. c2*rhsp
    return (x->-x./J).((rhsp,rhsu,rhsv))
end

Q = [p,u,v] # make arrays of arrays for mutability
resQ = [zeros(size(x)) for i in eachindex(Q)]
for i = 1:Nsteps
    for INTRK = 1:5
        rhsQ = rhs(Q,rd,md,c2)
        @. resQ = rk4a[INTRK]*resQ + dt*rhsQ
        @. Q    = Q + rk4b[INTRK]*resQ
    end

    if i%10==0 || i==Nsteps
        println("Number of time steps: $i out of $Nsteps")
    end
end

#plotting
gr(aspect_ratio=1,legend=false,
   markerstrokewidth=0,markersize=2)

@unpack Vp = rd # interpolation nodes
vv = Vp*Q[1]
scatter(Vp*x,Vp*y,vv,zcolor=vv,camera=(0,90))

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
module Stochastics

using QuadGK
using Roots

include("DescriptiveStatistics.jl")
include("Inference.jl")
include("Distributions.jl")
include("Regression.jl")

end

#TODO: Documentation and one example.

ext: jl

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
"""
    Base.isequal(c1::Cursor, c2::Cursor)::Bool

Test two `Cursor` objects for equality.

This "loose" form of `Cursor` equality is only determined by the content of the
expression passed to `@with` and reflects the expectation that the same query
(as passed to `withquery`) twice should produce `Cursor` objects that satisfy
`isequal`.
"""
function Base.isequal(c1::Cursor, c2::Cursor)::Bool
    isequal(c1.graph, c2.graph) || return false
    return true
end

"""
    source(c::Cursor)

Return the data source(s) over which `c` is to be collected.
"""
source(c::Cursor) = source(c.graph)

# TODO: Better wording of the following
"""
    graph(q::Cursor)

Return the `QueryNode` graph representation of the query that produced `q`.
"""
graph(c::Cursor) = c.graph

"""
    Base.collect(q::Cursor)

Collect a query against the source wrapped in the base `DataNode` of
`q.graph`.
"""
Base.collect(c::Cursor) = collect(c.graph)


text:
"""
    solve_mf(model, matrices, lap_eig, experiment)

Solve for magnetization using Matrix Formalism.
"""
function solve_mf(model, matrices, lap_eig, experiment)

    # Measure time of function evaluation
    starttime = Base.time()

    # Extract parameters
    @unpack mesh, ρ = model
    @unpack M, M_cmpts = matrices
    @unpack directions, sequences, values, values_type = experiment.gradient
    @unpack ninterval = experiment.mf

    # Laplace eigenmodes
    λ = lap_eig.values
    ϕ = lap_eig.funcs
    Ax = lap_eig.moments
    T = lap_eig.massrelax

    # Sizes
    ndirection = size(directions, 2)
    namplitude = length(values)
    nsequence = length(sequences)
    npoint_cmpts = size.(mesh.points, 2)
    ncompartment = length(mesh.points)

    qvalues, bvalues = get_values(experiment.gradient)

    # Compartment index ranges
    inds_start = cumsum([1; npoint_cmpts[1:end-1]])
    inds_stop = cumsum(npoint_cmpts[1:end])
    get_inds(icmpt) = inds_start[icmpt]:inds_stop[icmpt]

    # Global mass matrix
    M = blockdiag(M_cmpts...)

    # Create initial conditions (enforce complex values)
    ρ = vcat(fill.(complex.(ρ), npoint_cmpts)...)

    # Project initial spin density onto Laplace eigenfunction basis
    ν₀ = ϕ' * (M * ρ)

    # Allocate arrays
    signal = zeros(ComplexF64, ncompartment, namplitude, nsequence, ndirection)
    signal_allcmpts = zeros(ComplexF64, namplitude, nsequence, ndirection)
    magnetization =
        Array{Matrix{ComplexF64}, 4}(undef, ncompartment, namplitude, nsequence, ndirection)
    itertimes = zeros(namplitude, nsequence, ndirection)

    # Laplace operator in Laplace eigenfunction basis
    L = diagm(λ)

    # Iterate over gradient amplitudes, time profiles and directions
    for iamp = 1:namplitude, iseq = 1:nsequence, idir = 1:ndirection

        # Measure time of iteration
        itertime = time()

        # Gradient amplitude
        q = qvalues[iamp, iseq]
        b = bvalues[iamp, iseq]

        # Time profile
        f = sequences[iseq]

        # Gradient direction
        dir = directions[:, idir]

        # Display state of iterations
        @printf "Solving MF with size %d and neig = %d\n" (sum(npoint_cmpts)) length(λ)
        @printf "  Direction %d of %d: g = [%.2f, %.2f, %.2f]\n" idir ndirection dir...
        @printf "  Sequence %d of %d: f = %s\n" iseq nsequence f
        @printf "  Amplitude %d of %d: q = %g, b = %g\n" iamp namplitude q b

        # Gradient direction dependent finite element matrix
        A = sum(dir[i] * Ax[:, :, i] for i = 1:3)

        # Create array for magnetization coefficients
        ν = copy(ν₀)

        if isa(f, PGSE)
            # Constant Bloch-Torrey operator in Laplace eigenfunction basis
            K = L + T + im * q * A
            expmv!(-f.δ, K, ν)
            expmv!(-(f.Δ - f.δ), L + T, ν)
            expmv!(-f.δ, K', ν)
            # edK = exp(-f.δ * K)
            # edL = exp(-(f.Δ - f.δ) * (L + T))
            # ν = edK' * (edL * (edK * ν))
        else
            # Bloch-Torrey operator in Laplace eigenfunction basis for given
            # time profile value
            K(fᵢ) = L + T + im * q * fᵢ * A
            t = LinRange(0, echotime(f), ninterval + 1)
            for i = 1:ninterval
                δᵢ = t[i+1] - t[i]
                fᵢ = (f(t[i+1]) + f(t[i])) / 2
                # ν .= exp(-δᵢ * K(fᵢ)) * ν
                expmv!(-δᵢ, K(fᵢ), ν)
            end
        end

        # Final magnetization coefficients in finite element nodal basis
        mag = ϕ * ν

        # Store results
        for icmpt = 1:ncompartment
            inds = get_inds(icmpt)

            # Store magnetization in compartment
            magnetization[icmpt, iamp, iseq, idir] = mag[inds,:]

            # Integrate magnetization over compartment
            signal[icmpt, iamp, iseq, idir] = sum(M_cmpts[icmpt] * mag[inds])
        end

        itertimes[iamp, iseq, idir] = Base.time() - itertime
    end

    signal_allcmpts[:] = sum(signal, dims = 1)

    totaltime = Base.time() - starttime

    (; magnetization, signal, signal_allcmpts, totaltime, itertimes)
end


text:
# load code
include("flow_su2_self.jl")
include("flow_su2_s.jl")
include("flow_su2_t.jl")
include("flow_su2_u.jl")


text:
abstract type PR78AlphaModel <: AlphaModel end

struct PR78AlphaParam <: EoSParam
    acentricfactor::SingleParam{Float64}
end

@newmodelsimple PR78Alpha PR78AlphaModel PR78AlphaParam

export PR78Alpha
function PR78Alpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    acentricfactor = SingleParam(params["w"],"acentric factor")
    packagedparams = PR78AlphaParam(acentricfactor)
    model = PR78Alpha(packagedparams, verbose=verbose)
    return model
end

function α_function(model::CubicModel,V,T,z,alpha_model::PR78AlphaModel)
    Tc = model.params.Tc.values
    ω  = alpha_model.params.acentricfactor.values
    α = zeros(typeof(T),length(Tc))
    for i in @comps
        ωi = ω[i]
        m = ifelse(ωi<=0.491,
            evalpoly(ωi,(0.37464,1.54226,-0.26992)),
            evalpoly(ωi,(0.379642,1.487503,-0.164423,-0.016666)))
        Tr = T/Tc[i]
        α[i] = (1+m*(1-√(Tr)))^2
    end

    return α

end


text:
function calculate_experiment_results(chain, logger, include_vars, full_results)
    logdir = logger.logdir
    if full_results
        trace_plot = plot(chain[include_vars], seriestype=:traceplot)
        savefig(trace_plot, joinpath(logdir, "traceplots.png"))
        with_logger(logger) do
            @info "Trace Plot" trace_plot
        end

        # density
        density_plot = density(chain[include_vars])
        savefig(density_plot, joinpath(logdir, "densityplots.png"))
        with_logger(logger) do
            @info "Density Plot" density_plot
        end

        # Likelihood
        likelihood_plot = plot(chain[:lp])
        savefig(likelihood_plot, joinpath(logdir, "likelihoodplots.png"))
        with_logger(logger) do
            @info "Likelihood Plot" likelihood_plot
        end

        # Numerical errors
        if :numerical_error in keys(chain)
            num_error = get(chain, :numerical_error)
            prop = 100 * sum(num_error.numerical_error.data) / length(num_error.numerical_error.data)
            numerror_plot = plot(chain[:numerical_error])
            savefig(numerror_plot, joinpath(logdir, "numerrorplots.png"))
            with_logger(logger) do
                @info "Numerical Error Plot" numerror_plot
            end
        else
            prop = missing
        end
    end

    # summary statistics
    sum_stats = describe(chain[include_vars])
    param_names = sum_stats[1][:,1]
    param_mean = sum_stats[1][:,2]
    param_sd = sum_stats[1][:,3]
    param_ess = sum_stats[1][:,6]
    param_rhat = sum_stats[1][:,7]

    CSV.write(
        joinpath(logdir, "sumstats.csv"),
        DataFrame(
            Parameter=param_names,
            Mean=param_mean,
            sd=param_sd,
            ess=param_ess,
            rhat=param_rhat,
            Num_error= ismissing(prop) ? missing : prop * ones(length(include_vars))
        )
    )
    if full_results

        # Save the preconditioner if one is available
        if :samplerstate in keys(chain.info)
            state = chain.info.samplerstate

            # Check if there's a Hamiltonian in the state
            if :hamiltonian in propertynames(state)
                h = state.hamiltonian

                # Get the metric
                mass = h.metric.M⁻¹
                mass_shape = size(mass)

                # If it's diagonal, make it dense
                m = if length(mass_shape) == 1
                    diagm(mass)
                else
                    mass
                end

                # Write it to disk
                writedlm(joinpath(logdir,"preconditioner.csv"), m, ',')
            end
        end

        # cumulative averages
        for var in include_vars
            values = chain[var]
            cum_average = MCMCChains.cummean(values)
            writedlm(
                joinpath(logdir, "cumaverage_$(var).csv"),
                cum_average,
                ','
            )
        end

        # Cumulative mean plots
        cummean_plot = meanplot(chain[include_vars])
        savefig(cummean_plot, joinpath(logdir, "cummean.png"))
        with_logger(logger) do
            @info "Cumulate Means Plots" cummean_plot
        end
    end

    # Log the ESS/sec and rhat.  Nice to show as summary results from tensorboard
    for (i, name) = enumerate(param_names)
        TensorBoardLogger.log_value(
            logger,
            "$(name)_ess_per_sec",
            param_ess[i],
        )
        TensorBoardLogger.log_value(
            logger,
            "$(name)_rhat",
            param_rhat[i],
        )
    end

end


text:
using Oceananigans.Operators: Δzᵃᵃᶜ, Δzᵃᵃᶠ
using Oceananigans.Architectures: device_event
import Oceananigans.Architectures: architecture

struct FourierTridiagonalPoissonSolver{G, B, R, S, β, T}
                          grid :: G
    batched_tridiagonal_solver :: B
                   source_term :: R
                       storage :: S
                        buffer :: β
                    transforms :: T
end

architecture(solver::FourierTridiagonalPoissonSolver) = architecture(solver.grid)

@kernel function compute_main_diagonals!(D, grid, λx, λy)
    i, j = @index(Global, NTuple)
    Nz = grid.Nz

    # Using a homogeneous Neumann (zero Gradient) boundary condition:
    D[i, j, 1] = -1 / Δzᵃᵃᶠ(i, j, 2, grid) - Δzᵃᵃᶜ(i, j, 1, grid) * (λx[i] + λy[j])

    @unroll for k in 2:Nz-1
        D[i, j, k] = - (1 / Δzᵃᵃᶠ(i, j, k+1, grid) + 1 / Δzᵃᵃᶠ(i, j, k, grid)) - Δzᵃᵃᶜ(i, j, k, grid) * (λx[i] + λy[j])
    end

    D[i, j, Nz] = -1 / Δzᵃᵃᶠ(i, j, Nz, grid) - Δzᵃᵃᶜ(i, j, Nz, grid) * (λx[i] + λy[j])
end

function FourierTridiagonalPoissonSolver(grid, planner_flag=FFTW.PATIENT)
    TX, TY, TZ = topology(grid)
    TZ != Bounded && error("FourierTridiagonalPoissonSolver can only be used with a Bounded z topology.")

    Nx, Ny, Nz = size(grid)

    # Compute discrete Poisson eigenvalues
    λx = poisson_eigenvalues(grid.Nx, grid.Lx, 1, TX())
    λy = poisson_eigenvalues(grid.Ny, grid.Ly, 2, TY())

    arch = architecture(grid)
    λx = arch_array(arch, λx)
    λy = arch_array(arch, λy)

    # Plan required transforms for x and y
    sol_storage = arch_array(arch, zeros(complex(eltype(grid)), size(grid)...))
    transforms = plan_transforms(grid, sol_storage, planner_flag)

    # Lower and upper diagonals are the same
    lower_diagonal = CUDA.@allowscalar [1 / Δzᵃᵃᶠ(1, 1, k, grid) for k in 2:Nz]
    lower_diagonal = arch_array(arch, lower_diagonal)
    upper_diagonal = lower_diagonal

    # Compute diagonal coefficients for each grid point
    diagonal = arch_array(arch, zeros(Nx, Ny, Nz))
    event = launch!(arch, grid, :xy, compute_main_diagonals!, diagonal, grid, λx, λy, dependencies=device_event(arch))
    wait(device(arch), event)

    # Set up batched tridiagonal solver
    btsolver = BatchedTridiagonalSolver(grid;
                                        lower_diagonal = lower_diagonal,
                                              diagonal = diagonal,
                                        upper_diagonal = upper_diagonal)

    # Need buffer for index permutations and transposes.
    buffer_needed = arch isa GPU && Bounded in (TX, TY)
    buffer = buffer_needed ? similar(sol_storage) : nothing

    # Storage space for right hand side of Poisson equation
    rhs = arch_array(arch, zeros(complex(eltype(grid)), size(grid)...))

    return FourierTridiagonalPoissonSolver(grid, btsolver, rhs, sol_storage, buffer, transforms)
end

function solve!(x, solver::FourierTridiagonalPoissonSolver, b=nothing)
    !isnothing(b) && set_source_term!(solver, b) # otherwise, assume source term is set correctly

    arch = architecture(solver)
    ϕ = solver.storage

    # Apply forward transforms in order
    [transform!(solver.source_term, solver.buffer) for transform! in solver.transforms.forward]

    # Solve tridiagonal system of linear equations in z at every column.
    solve!(ϕ, solver.batched_tridiagonal_solver, solver.source_term)

    # Apply backward transforms in order
    [transform!(ϕ, solver.buffer) for transform! in solver.transforms.backward]

    ϕ .= real.(ϕ)

    # Set the volume mean of the solution to be zero.
    # Solutions to Poisson's equation are only unique up to a constant (the global mean
    # of the solution), so we need to pick a constant. We choose the constant to be zero
    # so that the solution has zero-mean.
    ϕ .= ϕ .- mean(ϕ)

    copy_event = launch!(arch, solver.grid, :xyz, copy_real_component!, x, ϕ, dependencies=device_event(arch))
    wait(device(arch), copy_event)

    return nothing
end

"""
    set_source_term!(solver, source_term)

Sets the source term in the discrete Poisson equation `solver`
to `source_term` by multiplying it by the vertical grid spacing at z cell centers.
"""
function set_source_term!(solver::FourierTridiagonalPoissonSolver, source_term)
    grid = solver.grid
    arch = architecture(solver)
    solver.source_term .= source_term

    event = launch!(arch, grid, :xyz, multiply_by_Δzᵃᵃᶜ!, solver.source_term, grid, dependencies=Event(device(arch)))
    wait(device(arch), event)

    return nothing
end

@kernel function multiply_by_Δzᵃᵃᶜ!(a, grid)
    i, j, k = @index(Global, NTuple)
    a[i, j, k] *= Δzᵃᵃᶜ(i, j, k, grid)
end


text:
@discord_object struct Role
    id::Snowflake
    name::String
    color::Int
    hoist::Bool
    position::Int
    permissions::Int64
    managed::Bool
    mentionable::Bool
end


text:
using Subsequences
using Base.Test

a = "abc"
b = "abc"
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == "abc"
@test a_range == b_range == 1:3

a = "-abc-"
b = "_abc_"
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == "abc"
@test a_range == b_range == 2:4

a = "abc-abc----"
b = "___ab_c__a_bc"
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == "abcabc"
@test a_range == 1:7
@test b_range == 4:13

# -------

a = "abc"
b = "abc"
seq, a_range, b_range = longest_contiguous_subsequence(a, b)
@test a_range == b_range == 1:3

a = "-abc-"
b = "_abc_"
seq, a_range, b_range = longest_contiguous_subsequence(a, b)
@test seq == "abc"
@test a_range == b_range == 2:4

a = "abc-abc----"
b = "___ab_c__abc"
seq, a_range, b_range = longest_contiguous_subsequence(a, b)
@test seq == "abc"
@test a_range == 1:3
@test b_range == 10:12


text:
using Base.Test
using DataFrames

let
	null_df = DataFrame(0, 0)
	df = DataFrame(Int, 4, 3)

	# Assignment of rows
	df[1, :] = df[1, :]
	df[1:2, :] = df[1:2, :]

	# Broadcasting assignment of rows
	df[1, :] = 1

	# Assignment of columns
	df[1] = datazeros(4)

	# Broadcasting assignment of columns
	df[:, 1] = 1
	df[1] = 3
	df["x3"] = 2

	rbind(null_df)
	rbind(null_df, null_df)
	rbind(null_df, df)
	rbind(df, null_df)
	rbind(df, df)
	rbind(df, df, df)

	alt_df = deepcopy(df)
	rbind(df, alt_df)
	df[1] = datazeros(Int, nrow(df))
	# Fail on non-matching types
	rbind(df, alt_df)

	alt_df = deepcopy(df)
	colnames!(alt_df, ["A", "B", "C"])
	# Fail on non-matching names
	rbind(df, alt_df)

	# df[:, 1] = dvzeros(Int, nrow(df))
end


text:
module CryptoMiniSat

import CryptoMiniSat_jll

export CMSLit, CMS, new_vars, nvars, add_clause, add_xor_clause, solve, get_model, get_conflict, itersolve, Lit, NotLit, XorLit, XNorLit

import Base: show, convert, iterate, IteratorSize

const libcms = CryptoMiniSat_jll.libcryptominisat5

const L_TRUE = 0x00
const L_FALSE = 0x01
const L_UNDEF = 0x02

################################################################
# low-level interface, CMSLit or 0-based

# the main structure
mutable struct CMS
    ptr::Ptr{Cvoid}
end
convert(::Type{Ptr{Cvoid}}, p::CMS) = p.ptr
"""`CMS(variables = 0; num_threads = Sys.CPU_THREADS)`

creates a new CryptoMiniSat problem, with optional number of variables and threas."""
function CMS(nvars::Int = 0; num_threads::Int = Sys.CPU_THREADS)
    p = cmsat_new()
    num_threads > 1 && cmsat_set_num_threads(p, num_threads)
    nvars > 0 && cmsat_new_vars(p, nvars)
    finalizer(cmsat_free, p)
    p
end

function Base.show(io::IO, p::CMS)
    if p.ptr == C_NULL
        print(io, "Unallocated CMS problem")
    else
        print(io, cmsat_nvars(p), "-variable CMS problem")
    end
end

# literals
struct CMSLit
    x::UInt32
end
CMSLit(v::Int) = CMSLit(v > 0 ? UInt32(2v-2) : v < 0 ? UInt32(-2v-1) : error("0 is not a valid literal"))
CMSLit(v::Int, invert::Bool) = CMSLit(UInt32(2v + (invert ? 1 : 0)))
Base.convert(::Type{CMSLit}, x::Int) = CMSLit(x)
Base.convert(::Type{Int}, x::CMSLit) = sign(x) ? -1-var(x) : 1+var(x)
Base.convert(::Type{Int32}, x::CMSLit) = Int32(convert(Int, x))
var(l::CMSLit) = Int(l.x >> 1)
sign(l::CMSLit) = Bool(l.x & 0x01)

function Base.show(io::IO, l::CMSLit)
    print(io, "CMSLit(",var(l),",",sign(l),")")
end

# constructor, destructor
cmsat_new() = ccall((:cmsat_new, libcms), CMS, ())
cmsat_free(p::CMS) = (ccall((:cmsat_free, libcms), Cvoid, (CMS,), p); p.ptr = C_NULL; nothing)

# add variables, clauses
cmsat_nvars(p::CMS) = ccall((:cmsat_nvars, libcms), Cint, (CMS,), p)
cmsat_add_clause(p::CMS, clause::Vector{CMSLit}) = Bool(ccall((:cmsat_add_clause, libcms), UInt8, (CMS,Ptr{CMSLit},Csize_t), p, clause, length(clause)))
cmsat_add_xor_clause(p::CMS, vars::Vector{Cuint}, rhs::Bool) = Bool(ccall((:cmsat_add_xor_clause, libcms), UInt8, (CMS,Ptr{Cuint},Csize_t,Cuint), p, vars, length(vars), rhs))
cmsat_new_vars(p::CMS, n::Int) = ccall((:cmsat_new_vars, libcms), Cvoid, (CMS,Csize_t), p, n)

# solve
cmsat_solve(p::CMS) = ccall((:cmsat_solve, libcms), UInt8, (CMS,), p)
cmsat_solve_with_assumptions(p::CMS, assumptions::Vector{CMSLit}) = ccall((:cmsat_solve_with_assumptions, libcms), UInt8, (CMS,Ptr{CMSLit},Csize_t), p, assumptions, length(assumptions))
cmsat_get_model(p::CMS) = unsafe_wrap(Array, ccall((:cmsat_get_model, libcms), Tuple{Ptr{UInt8},Csize_t}, (CMS,), p)...)
cmsat_get_conflict(p::CMS) = unsafe_wrap(Array, ccall((:cmsat_get_conflict, libcms), Tuple{Ptr{CMSLit},Csize_t}, (CMS,), p)...)

cmsat_print_stats(p::CMS) = ccall((:cmsat_print_stats, libcms), Cvoid, (CMS,), p)
cmsat_set_num_threads(p::CMS, n::Int) = ccall((:cmsat_set_num_threads, libcms), Cvoid, (CMS,Cuint), p, n)
cmsat_set_verbosity(p::CMS, n::Int) = ccall((:cmsat_set_verbosity, libcms), Cvoid, (CMS,Cuint), p, n)
cmsat_set_default_polarity(p::CMS, polarity::Int) = ccall((:cmsat_set_default_polarity, libcms), Cvoid, (CMS,Cint), p, polarity)
cmsat_set_polarity_auto(p::CMS) = ccall((:cmsat_set_polarity_auto, libcms), Cvoid, (CMS,), p)
cmsat_set_no_simplify(p::CMS) = ccall((:cmsat_set_no_simplify, libcms), Cvoid, (CMS,), p)
cmsat_set_no_simplify_at_startup(p::CMS) = ccall((:cmsat_set_no_simplify_at_startup, libcms), Cvoid, (CMS,), p)
cmsat_set_no_equivalent_lit_replacement(p::CMS) = ccall((:cmsat_set_no_equivalent_lit_replacement, libcms), Cvoid, (CMS,), p)
cmsat_set_no_bva(p::CMS) = ccall((:cmsat_set_no_bva, libcms), Cvoid, (CMS,), p)
cmsat_set_no_bve(p::CMS) = ccall((:cmsat_set_no_bve, libcms), Cvoid, (CMS,), p)
cmsat_set_up_for_scalmc(p::CMS) = ccall((:cmsat_set_up_for_scalmc, libcms), Cvoid, (CMS,), p)
cmsat_set_yes_comphandler(p::CMS) = ccall((:cmsat_set_yes_comphandler, libcms), Cvoid, (CMS,), p)
cmsat_simplify(p::CMS, assumptions::Vector{CMSLit}) = Bool(ccall((:cmsat_simplify, libcms), UInt8, (CMS,Ptr{CMSLit},Csize_t), p, clause, length(clause)))
cmsat_set_max_time(p::CMS, max_time::Float64) = ccall((:cmsat_set_max_time, libcms), Cvoid, (CMS,Float64), p, max_time)

################################################################
# higher-level interface, 1-based

"""`new_vars(p::CMS, n::Int)`

Add `n` new variables to the system `p`"""
new_vars(p::CMS, n::Int) = cmsat_new_vars(p, n)

"""`nvars(p::CMS)::Int`

Return the number of variables in the system `p`"""
nvars(p::CMS) = cmsat_nvars(p)

"""`add_clause(p::CMS, clause::Vector{Union{CMSLit,Integer}})::Bool`

Add a clause to the system, namely the disjunction of the entries in `clause`.
Direct terms are integers are numbered from 1, inverted terms are negative numbered from -1. They may also be specified as `CMSLit(variable,invert)`."""
add_clause(p::CMS, clause::Vector{CMSLit}) = cmsat_add_clause(p,clause)
function add_clause(p::CMS, clause::Vector{T}) where T <: Integer
    c = convert(Vector{CMSLit},clause)
    m = -1
    for l in c
        m = max(m,var(l))
    end
    if m ≥ nvars(p)
        new_vars(p, m+1-nvars(p))
    end
    cmsat_add_clause(p,c)
end

"""`add_clauses(p::CMS, clauses::Vector)`

Adds all `clauses` to `p`."""
function add_clauses(p::CMS, clauses)
    for item in clauses
        add_clause(p, item)
    end
end

"""`add_xor_clause(p::CMS, clause::Vector{Union{CMSLit,Int}}, rhs::Bool)::Bool`

Adds an XOR clause to `p`. This is a clause of the form `x_1 ⊻ x_2 ⊻ ... ⊻ x_n = rhs`."""
add_xor_clause(p::CMS, clause::Vector{UInt32}, rhs::Bool = true) = cmsat_add_xor_clause(p, clause, rhs)
function add_xor_clause(p::CMS, clause::Union{Vector{CMSLit},Vector{T}}, rhs = true) where T <: Integer
    cc = UInt32[]
    m = -1
    for c in clause
        if isa(c, CMSLit)
            x = var(c)
            rhs ⊻= sign(c)
        elseif c>0
            x = c-1
        elseif c<0
            x = -c-1
            rhs = !rhs
        else
            error("0 is not a valid literal")
        end
        m = max(m,x)
        push!(cc,x)
    end
    if m ≥ nvars(p)
        new_vars(p, m+1-nvars(p))
    end
    cmsat_add_xor_clause(p, cc, rhs)
end

"""`add_xor_clauses(p::CMS, clauses::Vector)`

Adds all XOR `clauses` to `p`."""
function add_xor_clauses(p::CMS, clauses)
    for item in clauses
        add_xor_clause(p, item)
    end
end

Cbool2Bool(v) = (v == 0x02 ? :undefined : v == 0x01 ? false : true)
"""`solve(p::CMS, assumptions=[])`

returns whether `p` is solvable; possibly returns :undefined."""
solve(p::CMS) = Cbool2Bool(cmsat_solve(p))
function solve(p::CMS, assumptions::Vector{Union{CMSLit,Integer}})
    c = convert(Vector{CMSLit},assumptions)
    m = -1
    for l in c
        m = max(m,var(l))
    end
    if m ≥ nvars(p)
        new_vars(p, m+1-nvars(p))
    end
    Cbool2Bool(cmsat_solve_with_assumptions(p,c))
end

"""`get_model(p::CMS)::Vector{Bool}`

returns the variables in the current solution to `p`. This assumes that `solve(p)` was called and returned `true`."""
function get_model(p::CMS)
    map(Cbool2Bool,cmsat_get_model(p))
end

"""`get_conflict(p::CMS)::Vector{CMSLit}`

returns the conflicting variables in the current solution to `p`. This assumes that `solve(p)` was called and returned `false`"""
function get_conflict(p::CMS)
    cmsat_get_conflict(p)
end

################################################################
# top-level interface, creating / destructing CMS object in one go.
# the syntax is compatible with e.g. PicoSAT

function create_problem(clauses, vars, verbose, num_threads)
    numvars = -1
    for c=clauses, l=c
        numvars = max(numvars,abs(l))
    end
    if vars == -1
        vars = numvars
    end
    vars < numvars && error("There are more variables than specified")

    p = CMS(vars; num_threads = num_threads)
    cmsat_set_verbosity(p, verbose)
    add_clauses(p, clauses)
    p
end

"""
`solve(clauses; vars::Integer=-1, verbose::Integer=0, proplimit::Integer=0, num_threads::Int = Sys.CPU_THREADS)`
   - `vars` - the number of variables (will be inferred if not specified)
   - `verbose` - verbosity
   - `proplimit` - ignored
   - `num_threads` - number of CPU threads to use
Returns a solution if the problem is satisfiable.
Satisfiable solutions are represented as a vector of signed integers.
If the problem is not satisfiable the method returns an `:unsatisfiable` symbol.
```julia
julia> import CryptoMiniSat
julia> cnf = [[1, -5, 4], [-1, 5, 3, 4], [-3, -4]];
julia> CryptoMiniSat.solve(cnf)
5-element Array{Int64,1}:
 -1
 -2
 -3
 -4
 -5
```
"""
function solve(clauses::Vector{Vector{T}};
               vars::Integer=-1,
               verbose::Integer=0,
               proplimit::Integer=0,
               num_threads::Integer=Sys.CPU_THREADS) where T <: Integer
    p = create_problem(clauses, vars, verbose, num_threads)
    if solve(p)
        sol = get_model(p)
        return [sol[i] ? i : -i for i=1:length(sol)]
    else
        return :unsatisfiable
    end
end

struct CMSIterator
    p::CMS
end
"""
`itersolve(clauses; vars::Integer=-1, verbose::Integer=0, proplimit::Integer=0, num_threads::Integer=Sys.CPU_THREADS)`
   - `vars` - the number of variables
   - `verbose` - verbosity
   - `proplimit` - ignored
   - `num_threads` - number of CPU threads to use

Returns an iterable object over all solutions.

```julia
julia> import CryptoMiniSat
julia> cnf = [[1, -5, 4], [-1, 5, 3, 4], [-3, -4]];
julia> collect(CryptoMiniSat.itersolve(cnf))
18-element Array{Any,1}:
 [-1, -2, -3, -4, -5]
 [-1, 2, -3, -4, -5]
 [-1, 2, 3, -4, -5]
 [1, 2, 3, -4, 5]
 [-1, -2, -3, 4, -5]
 [1, -2, 3, -4, 5]
 [-1, -2, -3, 4, 5]
 [-1, 2, -3, 4, -5]
 [1, -2, -3, 4, 5]
 [1, 2, -3, 4, 5]
 [1, 2, -3, 4, -5]
 [-1, 2, -3, 4, 5]
 [1, -2, -3, 4, -5]
 [1, -2, 3, -4, -5]
 [-1, -2, 3, -4, -5]
 [1, 2, -3, -4, 5]
 [1, -2, -3, -4, 5]
 [1, 2, 3, -4, -5]
```
"""
function itersolve(clauses::Vector{Vector{T}};
                   vars::Integer=-1,
                   verbose::Integer=0,
                   proplimit::Integer=0,
                   num_threads::Integer=Sys.CPU_THREADS) where T <: Integer
    p = create_problem(clauses, vars, verbose, num_threads)
    return CMSIterator(p)
end

function Base.iterate(it::CMSIterator, state=nothing)
    if solve(it.p)
        v = get_model(it.p)
        sol = [v[i] ? i : -i for i=1:length(v)]

        # add inverse sol for next iter
        add_clause(it.p, -sol)
        return (sol, nothing)
    else
        return nothing
    end
end

IteratorSize(it::CMSIterator) = Base.SizeUnknown()

################################################################
# super-top level interface, in which the terms are arbitrary, and
# coded into integers by a Dictionary

struct Lit{T}
    x::T
    sign::Int
end

"""
`Lit(x,inverted::Bool = false)`
Creates a literal `x`
"""
Lit(x::T,inverted::Bool = false) where T = Lit{T}(x,inverted)

"""
`NotLit(x)`
Creates a negated literal `x`; this is equivalent to `Lit(x,true)`
"""
NotLit(x::T) where T = Lit{T}(x,true)

"""
`XorLit(x)`
Creates a literal `x`, to be used in an XOR clause
"""
XorLit(x::T) where T = Lit{T}(x,-2)

"""
`XNorLit(x)`
Creates a negated literal `x`, to be used in an XOR clause
"""
XNorLit(x::T) where T = Lit{T}(x,-1)

function Base.show(io::IO,x::Lit{T}) where T
    isodd(x.sign) && print(io,"¬")
    x.sign<0 && print(io,"⊻")
    show(io,x.x)
end

function create_problem(clauses::Vector{Vector{Lit{T}}}, verbose, num_threads) where T
    vardict = Dict{T,Int32}()
    varlist = T[]

    for c=clauses, l=c
        haskey(vardict,l.x) && continue
        push!(varlist,l.x)
        vardict[l.x] = length(varlist)
    end

    p = CMS(length(varlist); num_threads = num_threads)
    cmsat_set_verbosity(p, verbose)

    for c=clauses
        is_xor = any(l->l.sign<0,c)
        is_xor == all(l->l.sign<0,c) || error("Clause $c contains a mix of Xor and usual terms")
        if is_xor
            rhs = true
            newc = UInt32[]
            for l=c
                push!(newc, vardict[l.x]-1)
                rhs ⊻= Bool(l.sign&1)
            end
            add_xor_clause(p, newc, rhs)
        else
            newc = CMSLit[]
            for l=c
                push!(newc, CMSLit(vardict[l.x]-1,Bool(l.sign)))
            end
            add_clause(p, newc)
        end
    end
    (p, vardict, varlist)
end

"""
`solve(clauses; verbose::Integer=0, num_threads::Int = Sys.CPU_THREADS)`
   - `verbose` - verbosity
   - `num_threads` - number of CPU threads to use
Returns a solution if the problem is satisfiable.
Satisfiable solutions are represented as a dictionary.

The first argument is a vector of clauses. Each clause is a vector of
literals of, each encapsulated in `Lit(T)`, `NotLit(T)`, `XorLit(T)`
or `XNorLit(T)` for some type `T`. Each clause containing XorLit and
XNotLit terms is an XOR clause, and the number of satisfied clauses
must have the same parity as the number of NotLit and XNorLit for the
clause to be satified. Clauses containing Lit and NotLit terms are
usual disjunctive clauses.

If the problem is not satisfiable the method returns an `:unsatisfiable` symbol.
```julia
julia> import CryptoMiniSat
julia> # coffee with cream, tea with milk
julia> cnf = [[XorLit(:coffee),XorLit(:tea)],[NotLit(:coffee),Lit(:cream)],[NotLit(:tea),Lit(:milk)]];
julia> CryptoMiniSat.solve(cnf)
Dict{Symbol, Bool} with 3 entries:
  :cream  => 0
  :milk   => 1
  :tea    => 1
  :coffee => 0
```
"""
function solve(clauses::Vector{Vector{Lit{T}}};
               verbose::Integer=0,
               num_threads::Integer=Sys.CPU_THREADS) where T
    (p, vardict, varlist) = create_problem(clauses, verbose, num_threads)
    if solve(p)
        sol = get_model(p)
        return Dict(varlist[i] => sol[i] for i=1:length(sol))
    else
        return :unsatisfiable
    end
end

struct CMSDIterator{T}
    p::CMS
    d::Dict{T,UInt32}
    l::Vector{T}
end

"""
`itersolve(clauses; verbose::Integer=0, num_threads::Integer=Sys.CPU_THREADS)`
   - `verbose` - verbosity
   - `num_threads` - number of CPU threads to use

The first argument is a vector of clauses. Each clause is a vector of
literals of, each encapsulated in `Lit(T)`, `NotLit(T)`, `XorLit(T)`
or `XNorLit(T)` for some type `T`. Each clause containing XorLit and
XNotLit terms is an XOR clause, and the number of satisfied clauses
must have the same parity as the number of NotLit and XNorLit for the
clause to be satified. Clauses containing Lit and NotLit terms are
usual disjunctive clauses.

Returns an iterable object over all solutions.
```julia
julia> import CryptoMiniSat
julia> # coffee with cream, tea with milk
julia> cnf = [[XorLit(:coffee),XorLit(:tea)],[NotLit(:coffee),Lit(:cream)],[NotLit(:tea),Lit(:milk)]];
julia> collect(CryptoMiniSat.itersolve(cnf))
4-element Vector{Any}:
 Dict{Symbol, Bool}(:cream => 0, :milk => 1, :tea => 1, :coffee => 0)
 Dict{Symbol, Bool}(:cream => 1, :milk => 1, :tea => 1, :coffee => 0)
 Dict{Symbol, Bool}(:cream => 1, :milk => 1, :tea => 0, :coffee => 1)
 Dict{Symbol, Bool}(:cream => 1, :milk => 0, :tea => 0, :coffee => 1)
```
"""
function itersolve(clauses::Vector{Vector{Lit{T}}};
                   verbose::Integer=0,
                   num_threads::Integer=Sys.CPU_THREADS) where T
    (p, vardict, varlist) = create_problem(clauses, verbose, num_threads)
    return CMSDIterator{T}(p, vardict, varlist)
end

function Base.iterate(it::CMSDIterator, state=nothing)
    if solve(it.p)
        v = get_model(it.p)

        # add inverse sol for next iter
        add_clause(it.p, [v[i] ? -i : i for i=1:length(v)])
        return (Dict(it.l[i] => v[i] for i=1:length(v)), nothing)
    else
        return nothing
    end
end

IteratorSize(it::CMSDIterator) = Base.SizeUnknown()

end


text:
"""
    ExactReach

ExactReach performs exact reachability analysis to compute the output reachable set for a network.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: HPolytope
3. Output: AbstractPolytope

# Return
`ReachabilityResult`

# Method
Exact reachability analysis.

# Property
Sound and complete.

# Reference
[W. Xiang, H.-D. Tran, and T. T. Johnson,
"Reachable Set Computation and Safety Verification for Neural Networks with ReLU Activations,"
*ArXiv Preprint ArXiv:1712.08163*, 2017.](https://arxiv.org/abs/1712.08163)
"""
struct ExactReach <: Solver end

function solve(solver::ExactReach, problem::Problem)
    reach = forward_network(solver, problem.network, problem.input)
    return check_inclusion(reach, problem.output)
end

forward_layer(solver::ExactReach, layer::Layer, input) = forward_layer(solver, layer, convert(HPolytope, input))

function forward_layer(solver::ExactReach, layer::Layer, input::Vector{<:HPolytope})
    output = Vector{HPolytope}(undef, 0)
    for i in 1:length(input)
        input[i] = affine_map(layer, input[i])
        append!(output, forward_partition(layer.activation, input[i]))
    end
    return output
end

function forward_layer(solver::ExactReach, layer::Layer, input::HPolytope)
    input = affine_map(layer, input)
    return forward_partition(layer.activation, input)
end


text:
using DSP

"""
    dbp(x) = 10 × log₁₀(|x|)

The dB equivalent of the power `x`
"""
function dbp(x)
    return pow2db.(abs.(x))
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
type Arc
    initNode::Int
    termNode::Int
    capacity::Float64
    freeflowtime::Float64
    flow::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64, freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0.)

## Solve an inverse tarffic problem over polynomials of degree at most d

## Optionally use a regularizer from the poly kernel

using JuMP
using Gurobi
using Graphs
using Roots

polyEval(coeffs, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])

polyEval(coeffs::Array{Float64, 1}, pt) = sum([coeffs[i] * pt^(i-1) for i = 1:length(coeffs)])

bpacost(flow::Float64, capacity::Float64, freeflowtime::Float64) = freeflowtime*(1 + .15 * (flow/capacity)^4)
bpacost(flow::Float64, arc) = bpacost(flow, arc.capacity, arc.freeflowtime)
bpacost(arc::Arc) = bpacost(arc.flow, arc)

function setUpFitting(deg::Int, c::Float64)

	m = Model(solver=GurobiSolver(OutputFlag=false))

	@defVar(m, coeffs[1:deg+1])
	@defVar(m, Calphas[1:deg+1])

	#build the graham matrix; cf. Ref. [21] (Regularization Networks and Support Vector Machines), page 47
	samples = linspace(0, 1, deg + 1)
	k(x,y) = (c + x*y)^deg
	K = [ k(x,y) for x = samples, y=samples]
	K = convert(Array{Float64, 2}, K)
	#assert(rank(K) == deg+1)

	C = chol(K + 1e-6* eye(deg+1))
	for i=1:deg + 1
		@addConstraint(m, polyEval(coeffs, samples[i]) == sum{C[j, i] * Calphas[j], j=1:deg+1})
	end

	@defVar(m, reg_term >= 0)
	reg_term_ = QuadExpr(Calphas[:], Calphas[:], ones(deg+1), AffExpr())

	@addConstraint(m, reg_term >= reg_term_)

	return m, coeffs, reg_term

end

function fixCoeffs(m, fcoeffs, coeffs)
	for (fc, c) in zip(fcoeffs, coeffs[:])
		@addConstraint(m, fc == c)
	end
end

function addResid(m, coeffs, ys, demands, arcs, scaling)
    @defVar(m, resid)
	@defVar(m, dual_cost)
	@defVar(m, primal_cost)

	@addConstraint(m, dual_cost == sum{demands[(s,t)] * (ys[(s,t), t] - ys[(s,t), s]), (s,t)=keys(demands)})
	@addConstraint(m, primal_cost == sum{a.flow * a.freeflowtime * polyEval(coeffs, a.flow/a.capacity), a=values(arcs)})

	@addConstraint(m, resid >= (dual_cost - primal_cost) / scaling )
	@addConstraint(m, resid >= (primal_cost - dual_cost) / scaling )
	return resid
end

function addIncreasingCnsts(m, coeffs, arcs; TOL=0.)
	sorted_flows = sort([a.flow / a.capacity for a in values(arcs)])
	@addConstraint(m, polyEval(coeffs, 0) <= polyEval(coeffs, sorted_flows[1]))
	for i = 2:length(sorted_flows)
		@addConstraint(m, polyEval(coeffs, sorted_flows[i-1]) <= polyEval(coeffs, sorted_flows[i]) + TOL)
	end
    @addConstraint(m, coeffs[1] == 1)
end

#equates the total cost of the network to the true total cost
function normalize(m, coeffs, tot_true_cost::Float64, arcs)
	@addConstraint(m,
		sum{a.freeflowtime * a.flow * polyEval(coeffs, a.flow / a.capacity), a=values(arcs)} == tot_true_cost)
end

function normalize(m, coeffs, scaled_flow::Float64, cost::Float64)
	@addConstraint(m, polyEval(coeffs, scaled_flow) == cost)
end

function normalize(m, coeffs, scaled_flows::Array{Float64, 1}, avgCost::Float64)
    @addConstraint(m, sum{polyEval(coeffs, f), f=scaled_flows} == avgCost * length(scaled_flows))
end

function addNetworkCnsts(m, coeffs, demands, arcs, numNodes)
	@defVar(m, ys[keys(demands), 1:numNodes])
	for k = keys(arcs)
		a = arcs[k]
		rhs = a.freeflowtime * polyEval(coeffs, a.flow/a.capacity)
		for od in keys(demands)
			@addConstraint(m, ys[od, k[2]] - ys[od, k[1]] <= rhs)
		end
	end
	return ys
end

############
#Read in the demand file

using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport parameters_julia

out_dir = parameters_julia.out_dir
files_ID = parameters_julia.files_ID
month_w = parameters_julia.month_w
year = parameters_julia.year
instances_ = parameters_julia.instances_ID

instance = readstring(out_dir * "instance_comm.txt")

file = open(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * instance * ".txt")
demands = Dict{(Int64,Int64), Float64}()
s = 0
for line in eachline(file)
    if contains(line, "Origin")
        s = int(split(line)[2])
    else
        pairs = split(line, ";")
        for pair in pairs
            if !contains(pair, "\n")
                pair_vals = split(pair, ":")
                t, demand = int(pair_vals[1]), float(pair_vals[2])
                demands[(s,t)] = demand
            end
        end
    end
end
close(file)

############
#read in the arc files
arcs = Dict{(Int, Int), Arc}()
file = open(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_net_" * month_w * "_" * instance * ".txt")
inHeader=true
for line in eachline(file)
    if inHeader
        inHeader = !contains(line, "Init node")
        continue
    end
    vals = split(line, )
    arcs[(int(vals[1]), int(vals[2]))] = Arc(int(vals[1]), int(vals[2]), float(vals[3]), float(vals[5]))
end
close(file)

#=
###########
#read in the initial flows
file = open("../data_original/SiouxFallsFlow.txt")
ix = 0;
for line in eachline(file)
    ix +=1
    if ix ==1
        continue
    end
    vals = split(line)
    arcs[(int(vals[1]), int(vals[2]))].flow = float(vals[3])
end
close(file)
=#
##########
# Set up demand data and flow data
##########

flow_data = Array(Float64, length(arcs))
flows = Dict{(Int64,Int64), Float64}()
demand_data = Dict{(Int, Int), Array{Float64, 1}}()

numNodes = maximum(map(pair->pair[1], keys(demands)))
g = simple_inclist(numNodes, is_directed=true)
vArcs = Arc[]
for arc in values(arcs)
    add_edge!(g, arc.initNode, arc.termNode)
    push!(vArcs, arc)
end

for odpair in keys(demands)
    if ! haskey(demand_data, odpair)
        demand_data[odpair] = [demands[odpair], ]
    else
        push!(demand_data[odpair], demands[odpair])
    end
end

flow_data = [a.flow::Float64 for a in vArcs]

for a in vArcs
    flows[(a.initNode, a.termNode)] = a.flow
end

#flows

using JSON

#load node-link incidence
nodeLink = readall("node_link_incidence_Sioux.json");
nodeLink = JSON.parse(nodeLink);

link_label_dict = readall("link_label_dict_Sioux.json");
link_label_dict = JSON.parse(link_label_dict);

link_label_dict["1"]

int(split(link_label_dict["1"], ',')[1]), int(split(link_label_dict["1"], ',')[2])

flows[int(split(link_label_dict["1"], ',')[1]), int(split(link_label_dict["1"], ',')[2])]

#string(1)

function addResid_(m, coeffs, ys, demands_, demands, arcs, scaling)
    @defVar(m, resid)
	@defVar(m, dual_cost)
	@defVar(m, primal_cost)

    for (s,t)=keys(demands)
        @addConstraint(m, demands_[(s,t)] >= 0)
    end

	@addConstraint(m, dual_cost == sum{demands_[(s,t)] * (ys[(s,t), t] - ys[(s,t), s]), (s,t)=keys(demands)})
	@addConstraint(m, primal_cost == sum{a.flow * a.freeflowtime * polyEval(coeffs, a.flow/a.capacity), a=values(arcs)})

	@addConstraint(m, resid >= (dual_cost - primal_cost) / scaling )
	@addConstraint(m, resid >= (primal_cost - dual_cost) / scaling )
	return resid
end

##########
#Fitting Funcs
##########

function train_cy(lam::Float64, deg::Int, c::Float64, demands, flow_data, arcs; fcoeffs=nothing)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))
    m, coeffs, reg_term = setUpFitting(deg, c)

    addIncreasingCnsts(m, coeffs, arcs, TOL=1e-8)  #uses the original obs flows

    avgCost = mean( [bpacost(a.flow, a.capacity, 1.0) for a in values(arcs)] )
    normalize(m, coeffs, [a.flow / a.capacity for a in values(arcs)], avgCost)

    resids = Variable[]

    #copy the flow data over to the arcs
    for (ix, a) in enumerate(vArcs)
        a.flow = flow_data[ix]
    end

    #Dual Feasibility
    ys = addNetworkCnsts(m, coeffs, demands, arcs, numNodes)

    #add the residual for this data point
    push!(resids, addResid(m, coeffs, ys, demands, arcs, 1e6))

    if fcoeffs != nothing
        fixCoeffs(m, fcoeffs, coeffs)
    end
    @setObjective(m, Min, sum{resids[i], i = 1:length(resids)} + lam*reg_term)
    solve(m)
    return [getValue(coeffs[i]) for i =1:length(coeffs)], getValue(ys), getValue(resids)
end

#nodeLink["0-75"]

#demands

outfile = open("demands_Sioux.json", "w")

JSON.print(outfile, demands)

close(outfile)

##########
#Fitting Funcs
##########

function train_cd(lam::Float64, deg::Int, c::Float64, ys, flow_data, flows, nodeLink, arcs; fcoeffs=nothing)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))
    m, coeffs, reg_term = setUpFitting(deg, c)

    addIncreasingCnsts(m, coeffs, arcs, TOL=1e-8)  #uses the original obs flows

    avgCost = mean( [bpacost(a.flow, a.capacity, 1.0) for a in values(arcs)] )
    normalize(m, coeffs, [a.flow / a.capacity for a in values(arcs)], avgCost)

    resids = Variable[]

    #copy the flow data over to the arcs
    for (ix, a) in enumerate(vArcs)
        a.flow = flow_data[ix]
    end

    for k = keys(arcs)
		a = arcs[k]
		rhs = a.freeflowtime * polyEval(coeffs, a.flow/a.capacity)
		for od in keys(demands)
			@addConstraint(m, ys[od, k[2]] - ys[od, k[1]] <= rhs)
		end
	end

    @defVar(m, demands_[keys(demands)])

    for (s,t) = keys(demands)
        @addConstraint(m, demands_[(s,t)] >= 0)
        @addConstraint(m, demands_[(s,t)] - demands[(s,t)] <= demands[(s,t)] * 0.05)
        @addConstraint(m, demands[(s,t)] - demands_[(s,t)] <= demands[(s,t)] * 0.05)
    end

    numLinks = length(flows)  # number of links
    @defVar(m, odLinkFlow[keys(demands), 1:numLinks])
    for (s,t) = keys(demands)
        for k = 1:numLinks
            @addConstraint(m, odLinkFlow[(s,t), k] >= 0)
        end
    end

    for k = 1:numLinks
        linkFlowPartial = sum([odLinkFlow[(s,t), k] for (s,t) = keys(demands)])
        key = (int(split(link_label_dict["$(k-1)"], ',')[1]), int(split(link_label_dict["$(k-1)"], ',')[2]))
        @addConstraint(m, linkFlowPartial == flows[key])
    end

    for l = 1:numNodes
        for (s,t) = keys(demands)
            if s == t
                @addConstraint(m, demands_[(s,t)] == 0)
            else
                odLinkFlowPartial = sum([nodeLink["$(l-1)-$(k-1)"] * odLinkFlow[(s,t), k] for k = 1:numLinks])
                if (l == s)
                    @addConstraint(m, odLinkFlowPartial + demands_[(s,t)] == 0)
                elseif (l == t)
                    @addConstraint(m, odLinkFlowPartial - demands_[(s,t)] == 0)
                else
                    @addConstraint(m, odLinkFlowPartial == 0)
                end
            end
        end
    end

    #add the residual for this data point
    push!(resids, addResid_(m, coeffs, ys, demands_, demands, arcs, 1e6))

    if fcoeffs != nothing
        fixCoeffs(m, fcoeffs, coeffs)
    end
    @setObjective(m, Min, sum{resids[i], i = 1:length(resids)} + lam*reg_term)
    solve(m)
    return [getValue(coeffs[i]) for i =1:length(coeffs)], getValue(demands_), getValue(resids)
end

#demands[(1, 1)]

#include("trafficCval.jl")

coeffs_dict = Dict{(Int64,Float64,Float64),Array{Float64,1}}()

deg = 6
c = 3.41
lam = 1.

demands_0 = copy(demands)

fcoeffs, ys, resids = train_cy(lam, deg, c, demands_0, flow_data, arcs)
coeffs_dict[(deg, c, lam)] = fcoeffs

fcoeffs, demands_, resides_ = train_cd(lam, deg, c, ys, flow_data, flows, nodeLink, arcs)
for (s,t) = keys(demands)
    demands_0[(s,t)] = demands_[(s,t)]
end
coeffs_dict[(deg, c, lam)] = fcoeffs

using PyPlot

true_coeffs = [1, 0, 0, 0, .15]

fcoeffs = coeffs_dict[(6, 3.41, 1.)]

xs = linspace(0, 2, 20)
zs_true = map(x->polyEval(true_coeffs, x), xs)

zs = map(x->polyEval(fcoeffs, x), xs)

plot(xs, zs_true, "s-g", label="True")

plot(xs, zs, "^-m", label="deg=6")
legend(loc="upper left",fancybox="true")

grid("on")
xlim(-0.1, 1.6);
ylim(0.9, 2.0);

font1 = ["family"=>"serif","color"=>"darkred","weight"=>"normal","size"=>14]
xlabel("Scaled Flow", fontdict=font1)
ylabel("Scaled Cost", fontdict=font1)

savefig("fitting_Sioux.pdf")

#demands_0

#demands

# based on https://github.com/chkwon/TrafficAssignment.jl

include("load_network_uni-class.jl")

using Graphs

function create_graph(start_node, end_node)
    @assert Base.length(start_node)==Base.length(end_node)

    no_node = max(maximum(start_node), maximum(end_node))
    no_arc = Base.length(start_node)

    graph = simple_inclist(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin
        parent = state.parents[current]

        link_idx = link_dic[parent,current]

        if link_idx != 0
            x[link_idx] = 1
        end

        current = parent
    end

    return x
end

ta_data = load_ta_network("Sioux Falls");

# ta_data.travel_demand;

# unpacking data from ta_data
network_name = ta_data.network_name

number_of_zones = ta_data.number_of_zones
number_of_nodes = ta_data.number_of_nodes
first_thru_node = ta_data.first_thru_node
number_of_links = ta_data.number_of_links

start_node = ta_data.start_node
end_node = ta_data.end_node
capacity = ta_data.capacity
link_length = ta_data.link_length

free_flow_time = ta_data.free_flow_time
B = ta_data.B
power = ta_data.power
speed_limit = ta_data.speed_limit
toll = ta_data.toll
link_type = ta_data.link_type
number_of_zones = ta_data.number_of_zones
total_od_flow = ta_data.total_od_flow
# travel_demand = ta_data.travel_demand
od_pairs = ta_data.od_pairs

toll_factor = ta_data.toll_factor
distance_factor = ta_data.distance_factor

best_objective = ta_data.best_objective

travel_demand = zeros(24, 24)

for (s,t)=keys(demands)
    travel_demand[s,t] = demands_0[(s,t)]
end

# preparing a graph
graph = create_graph(start_node, end_node)
link_dic = sparse(start_node, end_node, 1:number_of_links);

function BPR(x)
    bpr = similar(x)
    for i=1:length(bpr)
        bpr[i] = free_flow_time[i] * polyEval(fcoeffs, x[i]/capacity[i])
#         bpr[i] = free_flow_time[i] * ( 1.0 + B[i] * (x[i]/capacity[i])^power[i] )
    end
    return bpr
end

function all_or_nothing(travel_time)
    state = []
    path = []
    x = zeros(size(start_node))

    for r=1:size(travel_demand)[1]
        # for each origin node r, find shortest paths to all destination nodes
        state = dijkstra_shortest_paths(graph, travel_time, r)

        for s=1:size(travel_demand)[2]
            # for each destination node s, find the shortest-path vector
            # load travel demand
            x = x + travel_demand[r,s] * get_vector(state, r, s, link_dic)
        end
    end

    return x
end

# Finding a starting feasible solution
travel_time = BPR(zeros(number_of_links))
xl = all_or_nothing(travel_time);

max_iter_no = 1e3
l = 1
#average_excess_cost = 1
tol = 1e-6

while l < max_iter_no
    l += 1

    xl_old = xl

    # Finding yl
    travel_time = BPR(xl)

    yl = all_or_nothing(travel_time)

    assert(yl != xl)

    xl = xl + (yl - xl)/l

    xl_new = xl

    relative_gap = norm(xl_new - xl_old, 1) / norm(xl_new, 1)

#     if l % 500 == 0
#         print("l = $l------------------------------------------------\n")
#         print("relative_gap is $relative_gap\n")
#     end

    if relative_gap < tol
        print("l = $l------------------------------------------------\n")
        print("relative_gap is $relative_gap\n")
        break
    end

end

tapFlows = Dict{(Int64,Int64),Float64}()

for i = 1:length(ta_data.start_node)
    key = (ta_data.start_node[i], ta_data.end_node[i])
    tapFlows[key] = xl[i]
end

#tapFlows

#flows

#fcoeffs

function sa(x, a)  # calculate the partial derivatives of c_a w.r.t. x_a
    assert(a <= length(x) && a >= 1)
    n = length(fcoeffs)
    dcdx = 0
    for i=2:n
        dcdx += (i-1) * fcoeffs[i] * (x[a]/capacity[a])^(i-2)
    end
    dcdx *= free_flow_time[a]/capacity[a]
    return dcdx
end

x = zeros(size(start_node))
for k = 1:length(x)
    key = (int(split(link_label_dict["$(k-1)"], ',')[1]), int(split(link_label_dict["$(k-1)"], ',')[2]))
    x[k] = tapFlows[key]
end

x

## Obtain $\left( {\frac{{\partial {c_a}\left( {{g^l}} \right)}}{{\partial {v_a}}};a \in \mathcal{A}} \right)$

saVec = similar(x)
for a = 1:length(x)
    saVec[a] = sa(x, a)
end

# saVec[1:5]

#load OD pair-route incidence
odPairRoute = readall("od_pair_route_incidence_Sioux.json");
odPairRoute = JSON.parse(odPairRoute);

#load link-route incidence
linkRoute = readall("link_route_incidence_Sioux.json");
linkRoute = JSON.parse(linkRoute);

#load OD pair labels
odPairLabel = readall("od_pair_label_dict_Sioux.json");
odPairLabel = JSON.parse(odPairLabel);

odPairLabel_ = readall("od_pair_label_dict__Sioux.json");
odPairLabel_ = JSON.parse(odPairLabel_);

# express the demand data as vector (array)

demandsVec = zeros(length(demands))

for i = 1:length(demandsVec)
    demandsVec[i] = demands_0[(odPairLabel_["$i"][1], odPairLabel_["$i"][2])]
end

demandsVec[1:5]

# convert the demand data into dictionary

demandsDict = similar(demands)

for key = keys(demands)
    demandsDict[key] = demandsVec[odPairLabel["($(key[1]), $(key[2]))"]]
end

# "1-200" in keys(odPairRoute)

numLinks = size(start_node)[1]
numRoutes = length(odPairRoute)
numODpairs = numNodes * (numNodes - 1)


text:
"""
    AbstractLindaParam{T}

Abstract representation of a parameter with valuetype `T`.
"""
abstract type AbstractLindaParam{T} end

"""
    get_param_name(::AbstractLindaParam)

Return name of parameter.
"""
function get_param_name end

"""
    get_param_type(::AbstractLindaParam)

Get parameter's value type.
"""
function get_param_type end

"""
    get_param_value(::AbstractLindaParam)

Return current value of parameter.
"""
function get_param_value end

"""
    set_param_default!(::AbstractLindaParam{T})

Set parameter to its default value.
"""
function set_param_default! end

"""
    set_param_value!(::AbstractLindaParam{T}, v::T)

Check if `v` is an admissible value for parameter and, if so, change parameter's
    value to `v`.
"""
function set_param_value! end

"""
    test_param_value(::AbstractLindaParam{T}, v::T)

Check whether value `v` is admissible for given parameter.
"""
function test_param_value end

"""
    LindaRealParam{T<:Real}

Container for numerical (real-valued) parameters.
"""
mutable struct LindaRealParam{T<:Real} <: AbstractLindaParam{T}
    name::Symbol  # Name of the parameter

    val::T  # Current parameter value
    min_val::T  # Minimum parameter value
    max_val::T  # Maximum parameter value
    def_val::T  # Default value

    LindaRealParam(name::Symbol, vdef::T, vmin::T, vmax::T) where{T<:Real} =
        new{T}(name, vdef, vmin, vmax, vdef)

end

const LindaIntParam = LindaRealParam{Int}
const LindaFloatParam = LindaRealParam{Float64}

get_param_name(par::LindaRealParam) = par.name

get_param_type(par::LindaRealParam{T}) where T = T

get_param_value(par::LindaRealParam{T}) where T = par.val

set_param_default!(par::LindaRealParam) = (par.val = par.def_val)

"""
    set_param_value!(p, v)

Set value of parameter `p` to `v`. Raises an error if `v` is not an admissible value.
"""
function set_param_value!(p::LindaRealParam, v::T) where{T<:Real}

    if test_param_value(p, v)
        p.val = v
    else
        error("$(p.name) must be between $(p.min_val) and $(p.max_val).")
    end

    return nothing
end

"""
    test_param_value(p, v)

Return whether `v` is an admissible value for parameter `p`
"""
test_param_value(p::LindaRealParam, v::T) where{T<:Real} = (p.min_val <= v <= p.max_val)


text:
@testset "construction" begin
    @test basetype(QQ{Int}) == ZZ{Int}
    @test lcunit(QQ(-1, 3)) == 1
    @test lcunit(QQ(0)) == 1
    @test copy(QQ(12//15)) == QQ(4, 5)
    @test QQ{Int32}(QQ(3, 4)) isa QQ{Int32}
    @test QQ{Int}(QQ(3, 4)) != nothing
    @test QQ{BigInt}(ZZ(11)) == 11
    @test ZZ(12) // ZZ(13) == QQ(12, 13)

    @test QQ(1, 2) + QQ{Int8}(1, 3) == 5//6
    @test QQ(1, 3) - ZZ(1) == -2//3
    @test QQ(1, 2) * 2 == 1
    @test QQ(1, 2) - 1//3 == 1//6

    @test QQ(1, 2) / QQ(1, 3) == 3//2
    @test -QQ(1) == -1
    @test divrem(QQ(1,2), QQ(1, 3)) == (3//2, 0//1)
    @test div(QQ(1,2), QQ(1, 3)) == 3//2
    @test iszero(rem(QQ(1,2), QQ(1, 3)))
    @test isunit(QQ(23, 24))
    @test !isunit(QQ(0))
    @test isone(QQ(11, 11))
    @test zero(QQ{Int}) == 0
    @test one(QQ{Int8}) == 1
    @test hash(QQ(13, 14)) == hash(13//14)

    @test sprint(show, QQ(12, 13)) == "12/13"
    @test sprint(show, QQ(12, 1)) === "12"

    @test QQ{Int}(3, 6) == 1//2
    @test QQ(3, 6) == 1//2
    @test QQ(2) == QQ(ZZ(2))
    @test QQ(Int8(2)) == QQ{BigInt}(2//1)
    @test QQ{Int8}(2) == QQ(2)
    @test QQ(1//2) == 1//2

    @test_throws ArgumentError QQ(0, 0)
    @test_throws ArgumentError QQ(1, typemin(Int))
    @test QQ(2, typemin(Int)) == QQ(-1, typemin(Int)÷(-2))
    @test QQ(1, typemin(Int16)) == QQ(-1, -Int(typemin(Int16)))
    @test_throws InexactError QQ(1, typemax(UInt))
    @test_throws InexactError QQ(-1, typemax(UInt))
    #= would be nice to have
    @test_broken QQ(-3, typemax(UInt)) == QQ(-1, Int(typemax(UInt)÷3))
    @test_broken QQ(3, typemax(UInt)) == QQ(1, Int(typemax(UInt)÷3))
    @test_broken QQ(typemax(UInt), -3) == QQ(-1, Int(typemax(UInt)÷3))
    @test_broken QQ(typemax(UInt), 3) == QQ(1, Int(typemax(UInt)÷3))
    =#
end

@testset "arithmetic" begin
    @test QQ(2//3) + QQ(1//2) == QQ(7//6)
end

@testset "gcd calculations" begin
    @test gcd(QQ(1//3), QQ(2//5)) == 1 // 15
    @test gcd(QQ((2*3)//(5*7)), QQ((2*5)//(3*7))) == QQ(2, 3*5*7)
    @test gcdx(QQ((2*3)//(5*7)), QQ((2*5)//(3*7))) == (QQ(2, 3*5*7), -11, 4)
    @test lcm(QQ((2*3)//(5*7)), QQ((2*5)//(3*7))) == QQ(2*3*5, 7)
end


text:
##Home to deprecated functions for future reference
# AKS 2021_02_22

function DHC_compute_old(image::Array{Float64,2}, filter_hash, filter_hash2)
    # Use 2 threads for FFT
    FFTW.set_num_threads(2)

    # array sizes
    (Nx, Ny)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    S20 = zeros(Float64, Nf, Nf)
    S12 = zeros(Float64, Nf, Nf)
    S2  = zeros(Float64, Nf, Nf)  # traditional 2nd order

    # allocate image arrays for internal use
    im_fdf_0_1 = zeros(Float64,           Nx, Ny, Nf)   # this must be zeroed!
    im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    norm_im ./= sqrt(Nx*Ny*S0[2])

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val2   = filter_hash2["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    P   = plan_ifft(im_fd_0)   # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval
                im_fdf_0_1[ind,f] = abs(zval)
            end
            S1[f] = S1tot/(Nx*Ny)  # image power
            im_rd_0_1[:,:,f] .= abs2.(P*zarr)
            zarr[f_i] .= 0
        end
    end
    append!(out_coeff, S1[:])

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    im_rd_0_1 .= sqrt.(im_rd_0_1)

    ## 2nd Order
    #Amat = reshape(im_fdf_0_1, Nx*Ny, Nf)
    #S12  = Amat' * Amat
    #Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
    #S20  = Amat' * Amat

    #append!(out_coeff, S20)
    #append!(out_coeff, S12)

    ## Traditional second order
    for f1 = 1:Nf
        thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
        # println("  f1",f1,"  sum(fft):",sum(abs2.(thisim))/Nx^2, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,f1])))
        # Loop over f2 and do second-order convolution
        for f2 = 1:Nf
            f_i = f_ind2[f2]  # CartesianIndex list for filter
            f_v = f_val2[f2]  # Values for f_i
            # sum im^2 = sum(|fft|^2/npix)
            S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny)
        end
    end
    append!(out_coeff, S2)

    return out_coeff
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
__precompile__()

"""
This package allows you to use a pairwise list as a symmetric matrix.

```julia
using PairwiseListMatrices
```
"""
module PairwiseListMatrices

# standard libraries
using DelimitedFiles
using LinearAlgebra
using Statistics

# external packages
using NamedArrays
using RecipesBase

export  @iterateupper,
        @iteratelist,
        @iteratediag,

        PairwiseListMatrix,
        hasdiagonal,
        getlist, getdiag,
        diagonal,
        lengthlist, ij2k,
        sum_nodiag, mean_nodiag,
        zscore!, zscore,
        getlabels,
        setlabels, setlabels!,
        from_table, to_table,
        to_dict,
        join,
        writedlm

include("macros.jl")
include("pairwiselistmatrix.jl")
include("plotrecipes.jl")

end


text:
module Events

using Dates
using Discord

function on_message_create(c::Client, e::MessageCreate)
    e.message.content == "ping" && reply(c, e.message, "pong")
end

function on_no_handler(c::Client, e::AbstractEvent)
    println("Received unhandled event: $(typeof(e))")
end

function main()
    c = Client(ENV["DISCORD_TOKEN"])

    # This handler will remain valid for one minute.
    add_handler!(c, MessageCreate, on_message_create; tag=:ping, expiry=Minute(1))
    # This handler will run on the first ten events which have no non-default handler.
    add_handler!(c, FallbackEvent, on_no_handler; tag=:fallback, expiry=10)
    # You can also use do syntax.
    add_handler!(c, MessageReactionAdd) do c, e
        println("User $(e.user_id) reacted to message $(e.message_id)")
    end

    open(c)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    c = Events.main()
    wait(c)
end


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

kwf1(ones; tens=0, hundreds=0) = ones + 10*tens + 100*hundreds
@testset "simple keyword args case" begin
    @test kwf1(2) == 2
    @test kwf1(2, hundreds=6) == 602
    @test kwf1(2, tens=6) == 62
    @test kwf1(1, hundreds=2, tens=7) == 271
    @test kwf1(3, tens=7, hundreds=2) == 273

    @test_throws MethodError kwf1()             # no method, too few args
    @test_throws MethodError kwf1(1, z=0)       # unsupported keyword
    @test_throws MethodError kwf1(1, 2)         # no method, too many positional args
end
kwf2(x, rest...; y=1) = (x, y, rest)
@testset "keyword args plus varargs" begin
    @test isequal(kwf2(0), (0, 1, ()))
    @test isequal(kwf2(0,1,2), (0, 1, (1,2)))
    @test isequal(kwf2(0,1,2,y=88), (0, 88, (1,2)))
    @test isequal(kwf2(0,y=88,1,2), (0, 88, (1,2)))
    @test_throws MethodError kwf2(0, z=1)
    @test_throws MethodError kwf2(y=1)
end
@testset "Issue #13919" begin
    test13919(x::Vararg{Int}; key=100) = (x, key)
    @test test13919(1, 1)[1] === (1, 1)
    @test test13919(1, 1)[2] === 100
    @test test13919(1, 1, key=10)[1] === (1, 1)
    @test test13919(1, 1, key=10)[2] === 10
end
@testset "keyword arg with declared type" begin
    kwf3(x; y::Float64 = 1.0) = x + y
    @test kwf3(2) == 3.0
    @test kwf3(2, y=3.0) == 5.0
    @test_throws TypeError kwf3(2, y=3)  # wrong type keyword
end
@testset "function with only keyword args" begin
    kwf4(;a=1,b=2) = (a,b)
    @test isequal(kwf4(), (1,2))
    @test isequal(kwf4(b=10), (1,10))
end
@testset "in-order evaluation of keyword args" begin
    kwf9(;read=true,write=!read) = (read,write)
    @test kwf9() === (true,false)
    @test kwf9(read=false) === (false,true)
    @test kwf9(write=true) === (true,true)
end
# rest keywords
kwdelegator(ones;kw...) = kwf1(ones;kw...)
@test kwdelegator(4,hundreds=8) == 804

@testset "optional positional args" begin
    opaf1(a,b=1,c=2,d=3) = (a,b,c,d)
    @test isequal(opaf1(0), (0,1,2,3))
    @test isequal(opaf1(0,2), (0,2,2,3))
    @test isequal(opaf1(0,2,4), (0,2,4,3))
    @test isequal(opaf1(0,2,4,6), (0,2,4,6))
    @test_throws MethodError opaf1()
    @test_throws MethodError opaf1(0,1,2,3,4)

    @testset "with varargs" begin
        opaf2(a=1,rest...) = (a,rest)
        @test isequal(opaf2(), (1,()))
        @test isequal(opaf2(2), (2,()))
        @test isequal(opaf2(2,3), (2,(3,)))
    end
    @testset "with keyword args" begin
        opkwf1(a=0,b=1;k=2) = (a,b,k)
        @test isequal(opkwf1(), (0,1,2))
        @test isequal(opkwf1(10), (10,1,2))
        @test isequal(opkwf1(10,20), (10,20,2))
        @test_throws MethodError opkwf1(10,20,30)
        @test isequal(opkwf1(10,20,k=8), (10,20,8))
        @test isequal(opkwf1(11;k=8),    (11, 1,8))
        @test isequal(opkwf1(k=8),       ( 0, 1,8))
    end
end
# dictionaries as keywords
@test kwf1(4; Dict(:hundreds=>9, :tens=>5)...) == 954

@testset "with inner function" begin
    function kwf_maker()
        f(;k=0) = 2k+1
    end
    kwf5 = kwf_maker()
    @test kwf5() == 1
    @test kwf5(k=2) == 5
    @test_throws MethodError kwf5(1)
end

extravagant_args(x,y=0,rest...;color="blue",kw...) = (x,y,rest,color,kwf1(6;tens=8,kw...))
@testset "with every feature!" begin
    @test isequal(extravagant_args(1), (1,0,(),"blue",86))
    @test isequal(extravagant_args(1;hundreds=7), (1,0,(),"blue",786))
    @test isequal(extravagant_args(1,2,3;Dict(:color=>"red", :hundreds=>3)...),
                  (1,2,(3,),"red",386))
    # passing junk kw container
    @test_throws BoundsError extravagant_args(1; Any[[]]...)
end
# passing empty kw container to function with no kwargs
@test sin(1.0) == sin(1.0; Dict()...)
f18845() = 2
@testset "issue #18845" begin
    @test (@eval sin(1.0; $([]...))) == sin(1.0)
    @test f18845(;) == 2
    @test f18845(; []...) == 2
    @test (@eval f18845(; $([]...))) == 2
end

@testset "keyword args with static parameters" begin
    kwf6(x; k::T=1) where {T} = T
    @test kwf6(1) === Int
    @test kwf6(1;k=2.5) === Float64

    kwf7(x::T; k::T=1) where {T} = T
    @test kwf7(2) === Int
    @test kwf7(1.5;k=2.5) === Float64
    @test_throws MethodError kwf7(1.5)
    @test_throws TypeError kwf7(1.5; k=2)
end
# try to confuse it with quoted symbol
kwf8(x::MIME{:T};k::T=0) where {T} = 0
@test kwf8(MIME{:T}()) === 0

macro TEST4538()
    quote
        function $(esc(:test4538))(x=1)
            return x
        end
    end
end

macro TEST4538_2()
    quote
        function $(esc(:test4538_2))(;x=1)
            return x
        end
    end
end

module Foo4538
macro TEST()
    quote
        function $(esc(:test4538_foo_2))(;x=1)
            return x
        end
    end
end
end

f4538_3(;x=1) = x
macro TEST4538_3()
    quote
        x = 2
        f4538_3(x=3)
    end
end
@testset "issue #4538" begin
    @TEST4538
    @test test4538() == 1
    @test test4538(2) == 2
    @TEST4538_2
    @test test4538_2() == 1
    @test_throws MethodError test4538_2(2)
    @test test4538_2(x=2) == 2

    # that, but in a module
    @Foo4538.TEST()
    @test test4538_foo_2() == 1
    @test test4538_foo_2(x=2) == 2

    @test (@TEST4538_3) == 3
end
# issue #4801
mutable struct T4801{X}
    T4801{X}(;k=0) where X = new()
end
@test isa(T4801{Any}(k=0), T4801{Any})

# issue #4974
function Foo4974(f; kwargs...)
    (f(), kwargs)
end

function test4974(;kwargs...)
    Foo4974(;kwargs...) do
        2
    end
end

@test test4974(a=1) == (2, pairs((a=1,)))

@testset "issue #7704, computed keywords" begin
    @test kwf1(1; :tens => 2) == 21
    let p = (:hundreds, 3),
        q = (:tens, 1)
        @test kwf1(0; p[1]=>p[2], q[1]=>q[2]) == 310
        @test kwf1(0; q[1]=>q[2], hundreds=4) == 410
    end
end
@testset "with anonymous functions, issue #2773" begin
    f = (x;a=1,b=2)->(x, a, b)
    @test f(0) === (0, 1, 2)
    @test f(1,a=10,b=20) === (1,10,20)
    @test f(0,b=88) === (0, 1, 88)
    @test_throws MethodError f(0,z=1)
end

@test ((a=2)->10a)(3) == 30
@test ((a=2)->10a)() == 20
@test ((a=1,b=2)->(a,b))() == (1,2)
@test ((a=1,b=2)->(a,b))(5) == (5,2)
@test ((a=1,b=2)->(a,b))(5,6) == (5,6)

@testset "issue #9948" begin
    f9948, getx9948 = let x
        x = 3
        h(;x=x) = x
        getx() = x
        h, getx
    end
    @test f9948() == 3
    @test getx9948() == 3
    @test f9948(x=5) == 5
    @test f9948() == 3
    @test getx9948() == 3
end
@testset "issue #17785 - handle all sources of kwargs left-to-right" begin
    g17785(; a=1, b=2) = (a, b)
    let opts = (:a=>3, :b=>4)
        @test g17785(; a=5, opts...) == (3, 4)
        @test g17785(; opts..., a=5) == (5, 4)
        @test g17785(; opts..., a=5, b=6) == (5, 6)
        @test g17785(; b=0, opts..., a=5) == (5, 4)
    end
end
# pr #18396, kwargs before Base is defined
@eval Core.Compiler begin
    f18396(;kwargs...) = g18396(;kwargs...)
    g18396(;x=1,y=2) = x+y
end
@test Core.Compiler.f18396() == 3
@testset "issue #7045, `invoke` with keyword args" begin
    f7045(x::Float64; y=true) = y ? 1 : invoke(f7045,Tuple{Real},x,y=y)
    f7045(x::Real; y=true) = y ? 2 : 3
    @test f7045(1) === 2
    @test f7045(1.0) === 1
    @test f7045(1, y=false) === 3
    @test f7045(1.0, y=false) === 3
end

struct T20804{T}
    y::T
end
(f::T20804)(;x=10) = f.y + x
@testset "issue #20804" begin
    x = T20804(4)
    @test x() == 14
    @test x(x=8) == 12
end
@testset "issue #21147" begin
    function f21147(f::Tuple{A}; kwargs...) where {B,A<:Tuple{B}}
        return B
    end
    @test f21147(((1,),)) === Int
    @test f21147(((1,),), k = 2) === Int
    function g21147(f::Tuple{A}, k = 2) where {B,A<:Tuple{B}}
        return B
    end
    @test g21147(((1,),)) === Int
    @test g21147(((1,),), 2) === Int
end
@testset "issue #21510" begin
    f21510(; Base.@nospecialize a = 2) = a
    @test f21510(a=:b) == :b
    @test f21510() == 2
end
@testset "issue #21518" begin
    a = 0
    f21518(;kw=nothing) = kw
    g21518() = (a+=1; f21518)
    g21518()(kw=1)
    @test a == 1
    g21518()(; :kw=>1)
    @test a == 2
end

@testset "issue #17240 - evaluate default expressions in nested scopes" begin
    a = 10
    f17240(;a=a-1, b=2a) = (a, b)
    @test f17240() == (9, 18)
    @test f17240(a=2) == (2, 4)
    @test f17240(b=3) == (9, 3)
    @test f17240(a=2, b=1) == (2, 1)
end

@testset "issue #9535 - evaluate all arguments left-to-right" begin
    counter = 0
    function get_next()
        counter += 1
        return counter
    end
    f(args...; kws...) = (args, values(kws))
    @test f(get_next(), a=get_next(), get_next(),
            b=get_next(), get_next(),
            [get_next(), get_next()]...; c=get_next(),
            (d = get_next(), f = get_next())...) ==
                ((1, 3, 5, 6, 7),
                 (a = 2, b = 4, c = 8, d = 9, f = 10))
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
###############################################################################
#
#   nmod_rel_series.jl : Power series over flint nmod integers mod n
#
###############################################################################

export nmod_rel_series, NmodRelSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::nmod_rel_series)
   val = pol_length(a) + valuation(a) - 1
   val < 0 && throw(DomainError())
   z = nmod_rel_series(modulus(a), Array{UInt}(0), 0, val, val)
   z.parent = parent(a)
   return z
end

elem_type(::Type{NmodRelSeriesRing}) = nmod_rel_series

parent_type(::Type{nmod_rel_series}) = NmodRelSeriesRing

base_ring(R::NmodRelSeriesRing) = R.base_ring

var(a::NmodRelSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::NmodRelSeriesRing) = R.prec_max

function normalise(a::nmod_rel_series, len::Int)
   if len > 0
      c = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
         (Ref{nmod_rel_series}, Int), a, len - 1)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         c = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
            (Ref{nmod_rel_series}, Int), a, len - 1)
      end
   end
   return len
end

function pol_length(x::nmod_rel_series)
   return ccall((:nmod_poly_length, :libflint), Int, (Ref{nmod_rel_series},), x)
end

precision(x::nmod_rel_series) = x.prec

function polcoeff(x::nmod_rel_series, n::Int)
   R = base_ring(x)
   if n < 0
      return R(0)
   end
   z = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
         (Ref{nmod_rel_series}, Int), x, n)
   return R(z)
end

zero(R::NmodRelSeriesRing) = R(0)

one(R::NmodRelSeriesRing) = R(1)

function gen(R::NmodRelSeriesRing)
   z = nmod_rel_series(modulus(R), [UInt(1)], 1, max_precision(R) + 1, 1)
   z.parent = R
   return z
end

modulus(R::NmodRelSeriesRing) = modulus(base_ring(R))

function deepcopy_internal(a::nmod_rel_series, dict::ObjectIdDict)
   z = nmod_rel_series(a)
   z.prec = a.prec
   z.val = a.val
   z.parent = parent(a)
   return z
end

function renormalize!(z::nmod_rel_series)
   i = 0
   zlen = pol_length(z)
   zval = valuation(z)
   zprec = precision(z)
   while i < zlen && iszero(polcoeff(z, i))
      i += 1
   end
   z.prec = zprec
   if i == zlen
      z.val = zprec
   else
      z.val = zval + i
      ccall((:nmod_poly_shift_right, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int), z, z, i)
   end
   return nothing
end

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

function show(io::IO, a::NmodRelSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{nmod_rel_series}) = true

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::nmod_rel_series)
   z = parent(x)()
   ccall((:nmod_poly_neg, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}),
               z, x)
   z.prec = x.prec
   z.val = x.val
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::nmod_rel_series, b::nmod_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:nmod_poly_set_trunc, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:nmod_poly_shift_left, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:nmod_poly_set_trunc, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:nmod_poly_shift_left, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function -(a::nmod_rel_series, b::nmod_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   lenz = max(lena, lenb)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:nmod_poly_set_trunc, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:nmod_poly_shift_left, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:nmod_poly_neg, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}), z, z)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:nmod_poly_set_trunc, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:nmod_poly_shift_left, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:nmod_poly_sub_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:nmod_poly_sub_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function *(a::nmod_rel_series, b::nmod_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z = parent(a)()
   z.val = a.val + b.val
   z.prec = prec + z.val
   if lena == 0 || lenb == 0
      return z
   end
   lenz = min(lena + lenb - 1, prec)
   ccall((:nmod_poly_mullow, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, a, b, lenz)
   renormalize!(z)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::nmod, y::nmod_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, UInt),
               z, y, data(x))
   renormalize!(z)
   return z
end

*(x::nmod_rel_series, y::nmod) = y * x

function *(x::fmpz, y::nmod_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), x, modulus(y))
   ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, UInt),
               z, y, r)
   renormalize!(z)
   return z
end

function *(x::UInt, y::nmod_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, UInt),
               z, y, mod(x, modulus(y)))
   renormalize!(z)
   return z
end

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

*(x::nmod_rel_series, y::UInt) = y * x

*(x::Integer, y::nmod_rel_series) = fmpz(x)*y

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

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::nmod_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   z = nmod_rel_series(x)
   z.prec = x.prec + len
   z.val = x.val + len
   z.parent = parent(x)
   return z
end

function shift_right(x::nmod_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   xval = valuation(x)
   z = parent(x)()
   if len >= xlen + xval
      z.prec = max(0, x.prec - len)
      z.val = max(0, x.prec - len)
   else
      z.prec = max(0, x.prec - len)
      z.val = max(0, xval - len)
      zlen = min(xlen + xval - len, xlen)
      ccall((:nmod_poly_shift_right, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, x, xlen - zlen)
      renormalize!(z)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::nmod_rel_series, prec::Int)
   prec < 0 && throw(DomainError())
   xlen = pol_length(x)
   xprec = precision(x)
   xval = valuation(x)
   if xprec + xval <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   if prec <= xval
      z = parent(x)()
      z.val = prec
      z.prec = prec
   else
      z.val = xval
      ccall((:nmod_poly_set_trunc, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, x, min(prec - xval, xlen))
   end
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::nmod_rel_series, b::Int)
   b < 0 && throw(DomainError())
   if isgen(a)
      z = parent(a)()
      z = setcoeff!(z, 0, UInt(1))
      z.prec = a.prec + b - 1
      z.val = b
   elseif pol_length(a) == 0
      z = parent(a)()
      z.prec = b*valuation(a)
      z.val = b*valuation(a)
   elseif pol_length(a) == 1
      z = parent(a)([polcoeff(a, 0)^b], 1,
                           (b - 1)*valuation(a) + precision(a), b*valuation(a))
      renormalize!(z)
      return z
   elseif b == 0
      return one(parent(a))
   else
      z = parent(a)()
      z.prec = a.prec + (b - 1)*valuation(a)
      z.val = b*valuation(a)
      ccall((:nmod_poly_pow_trunc, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int, Int),
               z, a, b, z.prec - z.val)
   end
   renormalize!(z)
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::nmod_rel_series, y::nmod_rel_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)
   if prec <= x.val && prec <= y.val
      return true
   end
   if x.val != y.val
      return false
   end
   xlen = normalise(x, min(pol_length(x), prec - x.val))
   ylen = normalise(y, min(pol_length(y), prec - y.val))
   if xlen != ylen
      return false
   end
   return Bool(ccall((:nmod_poly_equal_trunc, :libflint), Cint,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               x, y, xlen))
end

function isequal(x::nmod_rel_series, y::nmod_rel_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || x.val != y.val || pol_length(x) != pol_length(y)
      return false
   end
   return Bool(ccall((:nmod_poly_equal, :libflint), Cint,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               x, y, pol_length(x)))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::nmod_rel_series, y::nmod)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         z = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
                       (Ref{nmod_rel_series}, Int), x, 0)
         return data(y) == z
      else
         return false
      end
   else
      return iszero(data(y))
   end
end

==(x::nmod, y::nmod_rel_series) = y == x

function ==(x::nmod_rel_series, y::fmpz)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), y, modulus(x))
         z = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
                       (Ref{nmod_rel_series}, Int), x, 0)
         return r == z
      else
         return false
      end
   else
      r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), y, modulus(x))
      return r == UInt(0)
   end
end

==(x::fmpz, y::nmod_rel_series) = y == x

function ==(x::nmod_rel_series, y::UInt)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         r = mod(y, modulus(x))
         z = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
                       (Ref{nmod_rel_series}, Int), x, 0)
         return r == z
      else
         return false
      end
   else
      r = mod(y, modulus(x))
      return r == UInt(0)
   end
end

==(x::UInt, y::nmod_rel_series) = y == x

==(x::nmod_rel_series, y::Integer) = x == fmpz(y)

==(x::Integer, y::nmod_rel_series) = y == x

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

function divexact(x::nmod_rel_series, y::nmod_rel_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   yval = valuation(y)
   xval = valuation(x)
   if yval != 0
      if xval >= yval
         x = shift_right(x, yval)
         y = shift_right(y, yval)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec - x.val, y.prec - y.val)
   z = parent(x)()
   z.val = xval - yval
   z.prec = prec + z.val
   if prec != 0
      ccall((:nmod_poly_div_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, x, y, prec)
   end
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::nmod_rel_series, y::nmod)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.val = x.val
   r = inv(y)
   ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, UInt),
               z, x, data(r))
   return z
end

function divexact(x::nmod_rel_series, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), y, modulus(x))
   rinv = inv(base_ring(x)(r))
   ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, UInt),
               z, x, data(rinv))
   return z
end

function divexact(x::nmod_rel_series, y::UInt)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   r = mod(y, modulus(x))
   rinv = inv(base_ring(x)(r))
   ccall((:nmod_poly_scalar_mul_nmod, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, UInt),
               z, x, data(rinv))
   return z
end

divexact(x::nmod_rel_series, y::Integer) = divexact(x, fmpz(y))

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::nmod_rel_series)
   iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ainv.val = 0
   ccall((:nmod_poly_inv_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               ainv, a, a.prec)
   return ainv
end

###############################################################################
#
#   Special functions
#
###############################################################################

function exp(a::nmod_rel_series)
   if iszero(a)
      z = one(parent(a))
      z.prec = precision(a)
      z.val = valuation(a)
      return z
   end
   z = parent(a)()
   R = base_ring(a)
   vala = valuation(a)
   preca = precision(a)
   d = Array{nmod}(preca)
   c = vala == 0 ? polcoeff(a, 0) : R()
   d[1] = exp(c)
   len = pol_length(a) + vala
   z0 = R()
   for k = 1 : preca - 1
      s = R()
      for j = 1 : min(k + 1, len) - 1
         c = j >= vala ? polcoeff(a, j - vala) : z0
         s += j * c * d[k - j + 1]
      end
      !isunit(base_ring(a)(k)) && error("Unable to divide in exp")
      d[k + 1] = divexact(base_ring(a)(s), k)
   end
   z = parent(a)(d, preca, preca, 0)
   ccall((:_nmod_poly_set_length, :libflint), Void,
         (Ref{nmod_rel_series}, Int), z, normalise(z, preca))
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(x::nmod_rel_series)
  ccall((:nmod_poly_zero, :libflint), Void,
                   (Ref{nmod_rel_series},), x)
  x.prec = parent(x).prec_max
  return x
end

function fit!(x::nmod_rel_series, n::Int)
  ccall((:nmod_poly_fit_length, :libflint), Void,
                   (Ref{nmod_rel_series}, Int), x, n)
  return nothing
end

function setcoeff!(z::nmod_rel_series, n::Int, x::fmpz)
   r = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ref{fmpz}, UInt), x, modulus(z))
   ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
                (Ref{nmod_rel_series}, Int, UInt),
               z, n, r)
   return z
end

function setcoeff!(z::nmod_rel_series, n::Int, x::UInt)
   r = mod(x, modulus(z))
   ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
                (Ref{nmod_rel_series}, Int, UInt),
               z, n, r)
   return z
end

function setcoeff!(z::nmod_rel_series, n::Int, x::nmod)
   ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
                (Ref{nmod_rel_series}, Int, UInt),
               z, n, data(x))
   return z
end

function mul!(z::nmod_rel_series, a::nmod_rel_series, b::nmod_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z.val = a.val + b.val
   z.prec = prec + c.val
   lenz = min(lena + lenb - 1, prec)
   if lena <= 0 || lenb <= 0
      lenz = 0
   end
   ccall((:nmod_poly_mullow, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               z, a, b, lenz)
   renormalize!(z)
   return z
end

function addeq!(a::nmod_rel_series, b::nmod_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   modulus = modulus(a)
   if a.val < b.val
      z = nmod_rel_series(modulus)
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:nmod_poly_set_trunc, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:nmod_poly_shift_left, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               a, a, z, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:nmod_poly_truncate, :libflint), Void,
            (Ref{nmod_rel_series}, Int),
            a, max(0, lenz - a.val + b.val))
      ccall((:nmod_poly_shift_left, :libflint), Void,
            (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
            a, a, a.val - b.val)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               a, a, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               a, a, b, lenz)
   end
   a.prec = prec
   a.val = val
   renormalize!(a)
   return a
end

function add!(c::nmod_rel_series, a::nmod_rel_series, b::nmod_rel_series)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenc = max(lena, lenb)
   c.prec = prec
   ccall((:nmod_poly_add_series, :libflint), Void,
                (Ref{nmod_rel_series}, Ref{nmod_rel_series}, Ref{nmod_rel_series}, Int),
               c, a, b, lenc)
   return c
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{nmod_rel_series}, ::Type{T}) where {T <: Integer} = nmod_rel_series

promote_rule(::Type{nmod_rel_series}, ::Type{fmpz}) = nmod_rel_series

promote_rule(::Type{nmod_rel_series}, ::Type{nmod}) = nmod_rel_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::NmodRelSeriesRing)()
   z = nmod_rel_series(modulus(a))
   z.prec = a.prec_max
   z.val = a.prec_max
   z.parent = a
   return z
end

function (a::NmodRelSeriesRing)(b::Integer)
   if b == 0
      z = nmod_rel_series(modulus(a))
      z.prec = a.prec_max
      z.val = a.prec_max
   else
      z = nmod_rel_series(modulus(a), [fmpz(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::NmodRelSeriesRing)(b::fmpz)
   if iszero(b)
      z = nmod_rel_series(modulus(a))
      z.prec = a.prec_max
      z.val = a.prec_max
   else
      z = nmod_rel_series(modulus(a), [b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::NmodRelSeriesRing)(b::nmod)
   if iszero(b)
      z = nmod_rel_series(modulus(a))
      z.prec = a.prec_max
      z.val = a.prec_max
   else
      z = nmod_rel_series(modulus(a), [b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::NmodRelSeriesRing)(b::nmod_rel_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::NmodRelSeriesRing)(b::Array{fmpz, 1}, len::Int, prec::Int, val::Int)
   z = nmod_rel_series(modulus(a), b, len, prec, val)
   z.parent = a
   return z
end

function (a::NmodRelSeriesRing)(b::Array{UInt, 1}, len::Int, prec::Int, val::Int)
   z = nmod_rel_series(modulus(a), b, len, prec, val)
   z.parent = a
   return z
end

function (a::NmodRelSeriesRing)(b::Array{nmod, 1}, len::Int, prec::Int, val::Int)
   z = nmod_rel_series(modulus(a), b, len, prec, val)
   z.parent = a
   return z
end

###############################################################################
#
#   PowerSeriesRing constructor
#
###############################################################################

function PowerSeriesRing(R::NmodRing, prec::Int, s::AbstractString; model=:capped_relative, cached = true)
   S = Symbol(s)

   if model == :capped_relative
      parent_obj = NmodRelSeriesRing(R, prec, S, cached)
   elseif model == :capped_absolute
      error("Not implemented yet")
      # parent_obj = NmodAbsSeriesRing(R, prec, S, cached)
   end
   return parent_obj, gen(parent_obj)
end


text:
"""
	$(SIGNATURES)

For all colleges (to keep track of shared parameters).
"""
abstract type GradRuleSet <: ModelObject end

abstract type GradRuleSwitches end
abstract type GradRule end

@common_fields grad_rule_set_common begin
    objId :: ObjectId
    pvec :: ParamVector
end

include("grad_rule_all.jl");
include("grad_rule_simple.jl");
include("grad_rule_linear.jl");
include("grad_rule_logistic.jl");

## --------  For testing

make_test_grad_rule(switches :: GradRuleSwitches, iCollege :: Integer) =
	make_grad_rule(
		make_grad_rule_set(ObjectId(:test), switches),
		iCollege);

function grad_rule_test_switches()
	return (GradRuleSwitchesSimple(),
		GradRuleSwitchesLogistic(),
		GradRuleSwitchesLinear(useHcLevel = true),
		GradRuleSwitchesLinear(useHcLevel = true, byCollege = false),
		GradRuleSwitchesLinear(useHcLevel = false, hMin = 0.2))
end

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


text:
# This file was generated, do not modify it. # hide
map!(x -> round(x, digits=0), data_nmiss.commissioning_year, data_nmiss.commissioning_year);

# We can now calculate plant age for each plant (worth remembering that the dataset only contains active plants)

current_year = fill!(Array{Float64}(undef, size(data_nmiss)[1]), 2020);
data_nmiss[:, :plant_age] = current_year - data_nmiss[:, :commissioning_year];


text:
include("../include.jl")

using JuMP, Ipopt

function create_1(n_v::Int64)
  m = Model()
  @variable(m, 0.0 <= theta[i=1:n_v] <= pi)
  @variable(m, 0.0 <= r[i=1:n_v] <= 1.0)

  @NLobjective(m, Min, -0.5 * sum(r[i+1] * r[i] * sin(theta[i+1] - theta[i]) for i = 1:(n_v - 1)))
  for i = 1:(n_v-1)
    for j = (i+1):n_v
      @NLconstraint(m, r[i]^2 + r[j]^2 - 2 * r[i] * r[j] * cos(theta[i] - theta[j]) <= 1)
    end
    @NLconstraint(m, theta[i] <= theta[i+1])
  end
  @constraint(m, r[n_v] == 0.0)
  @constraint(m, theta[n_v] == pi)

  return m;
end

m = create_1(20)
nlp_raw = MathProgNLPModel(m);
my_pars = Class_parameters();
it = one_phase_solve(nlp_raw,my_pars);
it.point.x
#=
setsolver(m,IpoptSolver())
solve(m)
=#