text:
println("Dh) Test of the PhotoExcitationAutoIon module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
setDefaults("print summary: open", "zzz-PhotoExcitationAutoIon.sum")

if  true
    # 1s photoexcitation and subsequent autoionization of Li-like carbon
    asfSettings   = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    excSettings   = PhotoExcitationAutoion.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], true, true, true, true, ExpStokes(),
                                                    [SolidAngle(1.0, 0.0)], 0., 4, PathwaySelection())
    grid          = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 1.0e-2, rbox = 10.0)
    setDefaults("unit: rate", "1/s")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(6.),
                            initialConfigs      = [Configuration("1s^2 2s")],
                            intermediateConfigs = [Configuration("1s 2s 2p")],
                            finalConfigs        = [Configuration("1s^2")],
                            process = PhotoExcAuto(),  processSettings = excSettings  )

    wb = perform(wa)
end
setDefaults("print summary: close", "")


text:
@testset "Strided" begin

    let rng = MersenneTwister(123456), N = 100

        # Test strided matrix-matrix multiplication sensitivities.
        for (TA, TB, tCA, tDA, CA, DA, tCB, tDB, CB, DB) in Nabla.strided_matmul
            A, B, VA, VB = randn.(Ref(rng), [N, N, N, N], [N, N, N, N])
            @test check_errs(*, A * B, (A, B), (VA, VB))
        end
    end
end

ext: jl


text:
function boxfig()
    fig = plot()
        phont = "Times New Roman"
        phontsize = 14
        Gaps0 = collect(Dict_AffProp[d*"_mem0_Gap"][end] for d in dists)
        Gaps1 = collect(Dict_AffProp[d*"_mem1_Gap"][end] for d in dists)
        GapsInf = collect(Dict_AffProp[d*"_memInf_Gap"][end] for d in dists)
        Time0 = collect(Dict_AffProp[d*"_mem0_TimeVec"][end] for d in dists)
        Time1 = collect(Dict_AffProp[d*"_mem1_TimeVec"][end] for d in dists)
        TimeInf = collect(Dict_AffProp[d*"_memInf_TimeVec"][end] for d in dists)

        fig = boxplot(["Mem: 0" "Mem: 1" "Mem: Inf"], [Gaps0 Gaps1 GapsInf], leg = false, title = "Gap", line =1, yaxis = :log)
        display(fig)
        fig = boxplot(["Mem: 0" "Mem: 1" "Mem: Inf"], [Time0 Time1 TimeInf], leg = false, title = "Time", line =1)

        display(fig)
        # savefig(fig, figsDir*"/AffProp_mem$(mem).png")
    end
    # plotAffProp(Dict_AffProp, chosen = DictAffProp, ymin= 0.0005)
    boxfig()

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
### Hastings & Powell 1991
## Reproducing figure 4
using LaTeXStrings

### Figure 4 (Bifurcation diagram)

## Define initial parameters
p = a1, a2, b1, b2, d1, d2 = [5.0, 0.1, 2.0, 2.0, 0.4, 0.01]
u0 = [0.76, 0.16, 9.9]
tspan = (0.0, 10000.0)
# b1 values interval
bs = collect(2.0:0.01:6.2)
# number of values to collect
N = 1000

## Loop to select all values of z for all b1 values
Z = zeros(Float64, (N, length(bs)))
Zmax = zeros(Float64, length(bs))
Zmin = zeros(Float64, length(bs))
errors = zeros(Float64, length(bs))
for j in 1:length(bs)
  # Change b1 value
  local p[3] = local b1 = bs[j]
  # Solve for given b1 value
  local prob = ODEProblem(parameterized_model!, u0, tspan, p)
  local sol = solve(prob)
  # Save values
  if length(sol[3,:]) > N
    Z[:,j] = sol[3, end-(N-1):end]
    Zmax[j] = maximum(Z[:, j])
    Zmin[j] = minimum(Z[:, j])
  else
    # Safeguard for when true solution is unstable
    errors[j] = j
  end
end
Z
Zmax
Zmin
errors[errors .> 0]

## Select only local maxima
Zmaxloc = zeros(Float64, size(Z))
for i in 2:(N-1), j in 1:length(bs)
  # Select maxima values(greater than neighbors)
  if Z[i, j] > Z[i-1, j] && Z[i, j] > Z[i+1, j]
    # Select first local maxima only with arbitrary interval
    Zrange = Zmax[j] - Zmin[j]
    if Z[i, j] > (Zmin[j] + 0.66*(Zrange))
       Zmaxloc[i, j] = Z[i, j]
    end
  end
end
Zmaxloc

## Visualize results
# Vector of b1 values repeated N times
B = repeat(bs, inner = N)
# Match local maxima with respective b1 value
points_maxloc = unique(hcat(B, vec(Zmaxloc)); dims = 1)
# Remove null values
pts_maxloc = points_maxloc[points_maxloc[:,2] .> 0.00, :]
# Fig.4 A)
fig4a = scatter(pts_maxloc[:, 1], pts_maxloc[:, 2],
                markersize = 2, markeralpha = 0.5,
                xlim = (2.2, 3.2), ylim = (9.5, 13.0),
                xlabel = L"b_1", ylabel = L"z_{max}", legend = false,
                grid = :none, c = :black, msw = 0.0, msc = :black,
                title = "a", titleloc = :left
                )
# Fig.4 B)
fig4b = scatter(pts_maxloc[:, 1], pts_maxloc[:, 2],
                markersize = 2, markeralpha = 0.5,
                xlim = (3.0, 6.5), ylim = (3.0, 10.0),
                xticks = 3.0:0.5:6.5, yticks = (3.0:1.0:10.0, string.(3.0:1.0:10.0)),
                xlabel = L"b_1", ylabel = L"z_{max}", legend = false,
                grid = :none, c = :black, msw = 0.0, msc = :black,
                title = "b", titleloc = :left
                )
# Fig.4 C)
fig4c = scatter(pts_maxloc[:, 1], pts_maxloc[:, 2],
                markersize = 2, markeralpha = 0.5,
                xlim = (2.25, 2.6), ylim = (11.4, 12.8),
                xticks = 2.25:0.05:2.60, yticks = 11.1:0.2:12.8,
                xlabel = L"b_1", ylabel = L"z_{max}", legend = false,
                grid = :none, c = :black, msw = 0.0, msc = :black,
                title = "c", titleloc = :left
                )
# Fig.4 - Combine A,B,C
fig4 = plot(fig4a, fig4b, fig4c, layout = (3, 1),
            size = (500, 900), titlefontsize = 11,
            tickfontsize = 8, guidefontsize = 10,
            left_margin = 10Plots.mm, right_margin = 5Plots.mm,
            bottom_margin = 2Plots.mm,
            framestyle = :box
            )

## Export figures
savefig(fig4, joinpath("..", "article", "figures", "fig4.pdf"))


text:
#################################################
# Parser to PowerModels Dict
# Ref: http://www.pserc.cornell.edu/matpower/manual.pdf
#################################################

function gen_dict(; name = nothing,
    ncost = 2, # Number of terms to interpolate in cost funcion
    qc1max= 0, # Maximum reactive power output at PC1 (MVAr)
    pg= 0, # Active Generation (Not Used by PowerModels,but expected)
    model= 2, # Cost model, 1 = piecewise linear, 2 = polynomial
    shutdown= 0, # Shutdown Cost
    startup= 0, # Startup Cost
    qc2max= 0, # Maximum reactive power output at PC2 (MVAr)
    ramp_agc= 0,
    qg= 0,  # Reactive Generation (Not Used by PowerModels,but expected)
    gen_bus= 2, # Bus Of Generator
    pmax= 1, # Maximum Real Power Generation (u.u) (Will be multi by mbase)
    ramp_10= 0,
    vg= 1, # Voltage magnitude setpoint (p.u.)
    mbase= 100, # Total MVA base of machine
    pc2= 0, # Upper real power output of PQ capability curve (MW)
    index= 1, # Index
    cost= [     # Cost terms
        2000,
        0
    ],
    qmax= 999, # Maximun Reacitive
    gen_status= 1, # Status Generator (On 1, OFF 0)
    qmin= -999, # Minimum Reactive
    qc1min= 0, # Minimum reactive power output at PC1 (MVAr)
    qc2min= 0, # Minimum reactive power output at PC2 (MVAr)
    pc1= 0, # Lower real power output of PQ capability curve (MW)
    ramp_q= 0,
    ramp_30= 0,
    pmin= 0, # Minimum Real Power Generation (u.u) (Will be multi by mbase)
    apf= 0 #[not] Area participation factor
    )
    gen = Dict("ncost"=> ncost,
        "qc1max"=> qc1max,
        "pg"=> pg,
        "model"=> model,
        "shutdown"=> shutdown,
        "startup"=> startup,
        "qc2max"=> qc2max,
        "ramp_agc"=> ramp_agc,
        "qg"=> qg,
        "gen_bus"=> gen_bus,
        "pmax"=> pmax,
        "ramp_10"=> ramp_10,
        "vg"=> vg,
        "mbase"=> mbase,
        "pc2"=> pc2,
        "index"=> index,
        "cost"=> cost,
        "qmax"=> qmax,
        "gen_status"=> gen_status,
        "qmin"=> qmin,
        "qc1min"=> qc1min,
        "qc2min"=> qc2min,
        "pc1"=> pc1,
        "ramp_q"=> ramp_q,
        "ramp_30"=> ramp_30,
        "pmin"=> pmin,
        "apf"=> apf
    )

    if name != nothing
        gen["name"] = name
    end
    return gen
end

function branch_dict(; name = nothing,
        br_r= 0.065, # Resistance (p.u.)
        rate_a= 1, #  MVA rating A (long term rating), set to 0 for unlimited
        shift= 0, #  Transformer phase shift angle (degrees), positive ⇒ delay
        br_x= 1, #  Reactance (p.u.)
        g_to= 0,
        g_fr= 0,
        b_fr= 0.225,
        f_bus= 1, # “from” bus number
        br_status= 1, # Initial branch status, 1 = in-service, 0 = out-of-service
        t_bus= 3, # “to” bus number
        b_to= 0.225,
        index= 1, # Index
        angmin= -1.0472, # minimum angle difference, θf − θt (rad)
        angmax= 1.0472, # maximum angle difference, θf − θt (rad)
        transformer= false, # Bool to indicate if it is transformer
        tap= 1) #  transformer off nominal turns ratio

    branch = Dict{String,Any}("br_r"=> br_r,
        "rate_a"=> rate_a,
        "shift"=> shift,
        "br_x"=> br_x,
        "g_to"=> g_to,
        "g_fr"=> g_fr,
        "b_fr"=> b_fr,
        "f_bus"=> f_bus,
        "br_status"=> br_status,
        "t_bus"=> t_bus,
        "b_to"=> b_to,
        "index"=> index,
        "angmin"=> angmin,
        "angmax"=> angmax,
        "transformer"=> transformer,
        "tap"=> tap
    )

    if name != nothing
        branch["name"] = name
    end
    return branch
end

function bus_dict(; name = nothing,
        zone= 1, #  loss zone (positive integer)
        bus_i= 1, # bus number (positive integer)
        bus_type= 1,  # bus type (1 = PQ, 2 = PV, 3 = ref, 4 = isolated)
        vmax= 1.1, # maximum voltage magnitude (p.u.)
        area= 1, # area number (positive integer)
        vmin= 0.9, # minimum voltage magnitude (p.u.)
        index= 1, # index
        va= 0, # voltage angle (degrees)
        vm= 1, # voltage magnitude (p.u.)
        base_kv= 0) # base voltage (kV)

    bus = Dict{String,Any}("zone"=> zone,
        "bus_i"=> bus_i,
        "bus_type"=> bus_type,
        "vmax"=> vmax,
        "area"=> area,
        "vmin"=> vmin,
        "index"=> index,
        "va"=> va,
        "vm"=> vm,
        "base_kv"=> base_kv
    )
    if name != nothing
        bus["name"] = name
    end
    return bus
end

function load_dict(;
        load_bus= 1, # bus where load is present
        status= 1, # status (ON: 1, OFF:0)
        pd= 1, # reactive power demand
        qd= 0.12*pd, # real power demand
        index= 1) # index

    lod = Dict{String,Any}( "load_bus"=> load_bus,
        "status"=> status,
        "qd"=> qd,
        "pd"=> pd,
        "index"=> index)
end

function _powermodels_dict_(; bus::Dict{String,Dict{String,Any}}= Dict{String,Dict{String,Any}}(),
                source_type::String= "matpower",
                name = "CaseX",
                cost_deficit = 1000,
                dcline::Dict{String,Any}= Dict{String,Any}(),
                source_version::Dict{String,Any}= Dict(
                    "major"=> 2,
                    "minor"=> 0,
                    "patch"=> 0,
                    "prerelease"=> [],
                    "build"=> []),
                gen::Dict{String,Dict{String,Any}}= Dict{String,Dict{String,Any}}(),
                branch::Dict{String,Dict{String,Any}}= Dict{String,Dict{String,Any}}(),
                storage::Dict{String,Dict{String,Any}}= Dict{String,Dict{String,Any}}(),
                baseMVA::Float64= 100.0,
                per_unit::Bool= true,
                shunt::Dict{String,Dict{String,Any}}= Dict{String,Dict{String,Any}}(),
                load::Dict{String,Dict{String,Any}}= Dict{String,Dict{String,Any}}())

    pm = Dict(  "name" => name,
        "cost_deficit" =>cost_deficit,
        "bus"=> bus,
        "source_type"=> source_type,
        "dcline"=> dcline,
        "source_version"=> source_version,
        "gen"=> gen,
        "branch"=> branch,
        "storage"=> storage,
        "baseMVA"=> baseMVA,
        "per_unit"=> per_unit,
        "shunt"=> shunt,
        "load"=> load
    )
    return pm

end

function powermodels_dict(;bus::Array{Dict{String,Any}}= Array{Dict{String,Any}}(undef,0),
    source_type::String= "matpower",
    name = "CaseX",
    cost_deficit = 1000,
    dcline::Dict{String,Any}= Dict{String,Any}(),
    source_version::Dict{String,Any}= Dict(
    "major"=> 2,
    "minor"=> 0,
    "patch"=> 0,
    "prerelease"=> [],
    "build"=> []),
    gen::Array{Dict{String,Any}}= Array{Dict{String,Any}}(undef,0),
    branch::Array{Dict{String,Any}}= Array{Dict{String,Any}}(undef,0),
    storage::Array{Dict{String,Any}}= Array{Dict{String,Any}}(undef,0),
    baseMVA::Float64= 100.0,
    per_unit::Bool= true,
    shunt::Array{Dict{String,Any}}= Array{Dict{String,Any}}(undef,0),
    load::Array{Dict{String,Any}}= Array{Dict{String,Any}}(undef,0))

    pm = _powermodels_dict_(bus= Dict(Pair.(string.(collect(1:size(bus,1))) , bus)),
        source_type= source_type,
        name= name,
        cost_deficit = cost_deficit,
        dcline= dcline,
        source_version= source_version,
        gen= Dict(Pair.(string.(collect(1:size(gen,1))),gen)),
        branch= Dict(Pair.(string.(collect(1:size(branch,1))),branch)),
        storage= Dict(Pair.(string.(collect(1:size(storage,1))),storage)),
        baseMVA= baseMVA,
        per_unit= per_unit,
        shunt= Dict(Pair.(string.(collect(1:size(shunt,1))),shunt)),
        load= Dict(Pair.(string.(collect(1:size(load,1))),load)))

end


text:
using YAML

struct FileNotFoundError <: Exception
    path::String
end

function getConfig(filepath::String, env::String)
    if ispath(filepath)
        config = YAML.load_file(filepath)[env]
        user = config["user"]  # default is usually "guest"
        password = config["password"]  # default is usually "guest"
        auth_params = compose_auth(user,password)
        host = config["host"]
        AmqpConnectionDef(host,auth_params)
    else
        throw(FileNotFoundError(filepath))
    end

end

function loadConfig(filepath::String = AMQP_CONN_FILE_PATH, env::String = AMQP_ENV)
    conn_def = getConfig(filepath, env)
    amqp_conn_define!(conn_def)
end


text:
N = 20
Npart = 10
t1 = 1.0
t2 = 0.2
U = 1.0
V1 = 0.5
nsw = 6
sweeps_args = [
  "maxdim" "mindim" "cutoff" "noise"
  50 10 1e-12 1e-7
  100 20 1e-12 1e-8
  200 20 1e-12 1e-10
  400 20 1e-12 0
  800 20 1e-12 1e-11
  800 20 1e-12 0
]


text:
module GizaExample
using Giza

function example1(devname::AbstractString="/xw",
                  prefix::AbstractString="Example 1")
    Giza.open_device(devname, prefix)
    Giza.set_environment(-10, 10, 0, 100, 0, 0)
    x = range(-10, stop=10, length=100);
    y = x.^2;
    Giza.line(x, y)
    Giza.annotate("B", 2.5, 0.5, 0.5, "Furlongs")
    Giza.annotate("L", 2.5, 0.5, 0.5, "Parsecs")
    Giza.close_device()
end

function example2(devname::AbstractString="/xw",
                  prefix::AbstractString="Example 2")
    Giza.open_device(devname, prefix)
    ls = Giza.get_line_style()
    lw = Giza.get_line_width()
    ci = Giza.get_colour_index()
    Giza.set_line_style(GIZA_LS_SOLID)
    Giza.set_line_width(1)
    Giza.set_colour_index(7)
    Giza.set_colour_representation(7, 1.0, 0.5, 0.0)
    Giza.set_environment(-10, 10, 0, 100, 0, 0)
    x = range(-10, stop=10, length=100);
    y = x.^2;
    Giza.line(x, y)
    Giza.annotate("B", 2.5, 0.5, 0.5, "Furlongs")
    Giza.annotate("L", 2.5, 0.5, 0.5, "Parsecs")
    Giza.set_line_style(ls)
    Giza.set_line_width(lw)
    Giza.set_colour_index(ci)
    #Giza.close_device()
end

function example3(devname::AbstractString="/xw",
                  prefix::AbstractString="Example 3")
    imgdim = 400
    windim = 640
    margin = (windim - imgdim)/2
    xmin = ymin = margin/imgdim
    xmax = ymax = (margin + imgdim)/imgdim

    Giza.open_device_size(devname, prefix, windim, windim, GIZA_UNITS_PIXELS)
    Giza.set_viewport(xmin, xmax, ymin, ymax)
    Giza.set_window(0, imgdim, 0, imgdim)
    Giza.set_box("BCTN", 0.0, 0, "BCNST", 0.0, 0)
end

function pli(img::AbstractMatrix{<:Real})
    devwidth, devheight = Giza.get_paper_size(Giza.GIZA_UNITS_PIXELS)
    imgwidth, imgheight = size(img)

    if imgwidth < devwidth && imgheight < devheight
        xmin = (devwidth - imgwidth)/2/devwidth
        xmax = xmin + imgwidth/devwidth
        ymin = (devheight - imgheight)/2/devheight
        ymax = ymin + imgheight/devheight
        Giza.stop_prompting()
        Giza.change_page()
        Giza.set_viewport(xmin, xmax, ymin, ymax)
        Giza.set_window(0, imgwidth, 0, imgheight)
        cimin, cimax = Giza.get_colour_index_range()
        amin, amax = extrema(img)
        if amin < amax && cimin < cimax
            a1 = (convert(Cdouble, cimax) - convert(Cdouble, cimin))
            a2 = (convert(Cdouble, amax) - convert(Cdouble, amin))
            a = a1/a2
            lvl = Array{Cint}(imgwidth, imgheight)
            @inbounds @simd for i in eachindex(lvl, img)
                lvl[i] = cimin + round(Cint, (img[i] - amin)*a)
            end
            Giza.draw_pixels(lvl, 0, imgwidth, 0, imgheight,
                             GIZA_EXTEND_NONE)
        end
    end
end

end # module


text:
include("5.1.1 (a+b arcsin(c x))^n.jl")
include("5.1.2 (d x)^m (a+b arcsin(c x))^n.jl")
include("5.1.3 (d+e x^2)^p (a+b arcsin(c x))^n.jl")
include("5.1.4 (f x)^m (d+e x^2)^p (a+b arcsin(c x))^n.jl")
include("5.1.5 u (a+b arcsin(c x))^n.jl")
include("5.1.6 Miscellaneous inverse sine.jl")

max_stars_repo_path: src/rules/5 Inverse trig functions/5.1 Inverse sine/.jl


text:
module LabJack

# LibLabJackM wrapper module
include("LJM/LJM.jl")
using .LJM
export LJM

# High-level API
include("util.jl")
include("device.jl")
include("analog.jl")
include("digital.jl")

# Common use LJM constants
export dtANY, T7, TSERIES, ctANY, USB, ETHERNET

# High-level types
export LJDevice
export lsdev, open, close, DefaultDev
export write_digital, read_digital
export write_analog, read_analog

global DEFAULT_LJDEVICE

function init_default()
    if lsdev(; quiet = true) !== nothing
        global DEFAULT_LJDEVICE = open(dtANY, ctANY)
    else
        @warn "No LabJack device found."
    end
end

function DefaultDev(dev::LJDevice)
    global DEFAULT_LJDEVICE = dev
end

end # module


text:
const VectorCone = Union{MOI.SecondOrderCone,
                         MOI.RotatedSecondOrderCone,
                         MOI.PowerCone,
                         MOI.DualPowerCone,
                         MOI.ExponentialCone,
                         MOI.DualExponentialCone}
const PositiveSemidefiniteCone = MOI.PositiveSemidefiniteConeTriangle
const LinearFunction = Union{MOI.SingleVariable,
                             MOI.VectorOfVariables,
                             MOI.ScalarAffineFunction,
                             MOI.VectorAffineFunction}
const AffineFunction = Union{MOI.ScalarAffineFunction,
                             MOI.VectorAffineFunction}
const VariableFunction = Union{MOI.ScalarAffineFunction,
                               MOI.VectorAffineFunction}

const ScalarLinearDomain = Union{MOI.LessThan{Float64},
                                 MOI.GreaterThan{Float64},
                                 MOI.EqualTo{Float64},
                                 MOI.Interval{Float64}}
const VectorLinearDomain = Union{MOI.Nonpositives,
                                 MOI.Nonnegatives,
                                 MOI.Reals,
                                 MOI.Zeros}
const LinearDomain = Union{ScalarLinearDomain,VectorLinearDomain}
const ScalarIntegerDomain = Union{MOI.ZeroOne, MOI.Integer}
################################################################################
# ADD CONSTRAINT ###############################################################
################################################################################

MOI.supports_constraint(m::MosekModel, ::Type{<:Union{MOI.SingleVariable, MOI.ScalarAffineFunction}}, ::Type{<:ScalarLinearDomain}) = true
MOI.supports_constraint(m::MosekModel, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction}}, ::Type{<:Union{VectorCone, PositiveSemidefiniteCone, VectorLinearDomain}}) = true
MOI.supports_constraint(m::MosekModel, ::Type{MOI.SingleVariable}, ::Type{<:ScalarIntegerDomain}) = true
#MOI.canaddconstraint(m::MosekModel, ::Type{<:Union{MOI.SingleVariable, MOI.ScalarAffineFunction}}, ::Type{<:ScalarLinearDomain}) = true
#MOI.canaddconstraint(m::MosekModel, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction}}, ::Type{<:Union{VectorCone, PositiveSemidefiniteCone, VectorLinearDomain}}) = true
#MOI.canaddconstraint(m::MosekModel, ::Type{MOI.SingleVariable}, ::Type{<:ScalarIntegerDomain}) = true

function MOI.add_constraint(m   :: MosekModel,
                            axb :: MOI.ScalarAffineFunction{Float64},
                            dom :: D) where {D <: MOI.AbstractScalarSet}

    # Duplicate indices not supported
    axb = MOIU.canonical(axb)

    N = 1
    conid = allocateconstraints(m,N)
    addlhsblock!(m,
                 conid,
                 fill(1, length(axb.terms)),
                 axb.terms)

    if length(m.c_constant) < length(m.c_block)
        append!(m.c_constant,
                zeros(Float64, length(m.c_block) - length(m.c_constant)))
    end
    conidxs = getindexes(m.c_block, conid)
    m.c_constant[conidxs] .= axb.constant

    addbound!(m, conid, conidxs, Float64[axb.constant], dom)
    conref = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, D}(UInt64(conid) << 1)
    select(m.constrmap, MOI.ScalarAffineFunction{Float64}, D)[conref.value] = conid

    return conref
end

function MOI.add_constraint(m   :: MosekModel,
                            axb :: MOI.VectorAffineFunction{Float64},
                            dom :: D) where { D <: MOI.AbstractVectorSet }

    # Duplicate indices not supported
    axb = MOIU.canonical(axb)

    N = MOI.dimension(dom)
    conid = allocateconstraints(m,N)
    addlhsblock!(m,
                 conid,
                 map(t -> t.output_index, axb.terms),
                 map(t -> t.scalar_term, axb.terms))

    conidxs = getindexes(m.c_block, conid)
    m.c_constant[conidxs] .= axb.constants

    addbound!(m, conid, conidxs, axb.constants, dom)
    conref = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, D}(UInt64(conid) << 1)
    select(m.constrmap, MOI.VectorAffineFunction{Float64}, D)[conref.value] = conid

    return conref
end

function trimapL(i, j, n)
    if i < j
        trimapL(j, i, n)
    else
        i + div((2n-j) * (j-1), 2)
    end
end

function trimapU(i::Integer, j::Integer)
    if i < j
        trimapU(j, i)
    else
        div((i-1)*i, 2) + j
    end
end

# Vectorized length for matrix dimension n
sympackedlen(n) = (n*(n+1)) >> 1
# Matrix dimension for vectorized length n
sympackeddim(n) = div(isqrt(1+8n) - 1, 2)

function _sympackedto(x, n, mapfrom, mapto)
    @assert length(x) == sympackedlen(n)
    y = similar(x)
    for i in 1:n, j in 1:i
        y[mapto(i, j)] = x[mapfrom(i, j)]
    end
    y
end
sympackedLtoU(x, n=sympackeddim(length(x))) = _sympackedto(x, n, (i, j) -> trimapL(i, j, n), trimapU)
sympackedUtoL(x, n=sympackeddim(length(x))) = _sympackedto(x, n, trimapU, (i, j) -> trimapL(i, j, n))

function sympackedUtoLidx(x::AbstractVector{<:Integer}, n)
    y = similar(x)
    map = sympackedLtoU(1:sympackedlen(n), n)
    for i in eachindex(y)
        y[i] = map[x[i]]
    end
    y
end

function MOI.add_constraint(m   :: MosekModel,
                            axb :: MOI.VectorAffineFunction{Float64},
                            dom :: PSDCone) where { PSDCone <: PositiveSemidefiniteCone }

    # Duplicate indices not supported
    axb = MOIU.canonical(axb)

    N = dom.side_dimension
    NN = sympackedlen(N)

    conid = allocateconstraints(m,NN)
    addlhsblock!(m,
                 conid,
                 sympackedUtoLidx(map(t -> t.output_index, axb.terms), N),
                 map(t -> t.scalar_term, axb.terms))
    conidxs = getindexes(m.c_block,conid)
    constant = sympackedUtoL(axb.constants, N)
    m.c_constant[conidxs] = constant

    addbound!(m,conid,conidxs,constant,dom)
    conref = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},PSDCone}(UInt64(conid) << 1)
    select(m.constrmap,MOI.VectorAffineFunction{Float64},PSDCone)[conref.value] = conid
    conref
end

################################################################################
# Variable constraints

# We allow following. Each variable can have
# - at most most upper and one lower bound
# - belong to at most one non-semidefinite cone
# - any number of semidefinite cones, which are implemented as ordinary constraints
# This is when things get a bit funky; By default a variable has no
# bounds, i.e. "free". Adding a `GreaterThan` or `Nonnegatives`
# constraint causes it to have a defined lower bound but no upper
# bound, allowing a `LessThan` or ``Nonpositives` constraint to be
# added later. Adding a `Interval` constraint defines both upper and
# lower bounds. Adding a `Reals` constraint will effectively be the
# same as an interval ]-infty;infty[, in that it will define both
# upper and lower bounds, and not allow those to be set afterwards.

domain_type_mask(dom :: MOI.Reals)        = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MOI.Interval)     = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MOI.EqualTo)      = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MOI.Zeros)        = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MOI.GreaterThan)  = boundflag_lower
domain_type_mask(dom :: MOI.Nonnegatives) = boundflag_lower
domain_type_mask(dom :: MOI.LessThan)     = boundflag_upper
domain_type_mask(dom :: MOI.Nonpositives) = boundflag_upper

domain_type_mask(dom :: MOI.SecondOrderCone)        = boundflag_cone
domain_type_mask(dom :: MOI.RotatedSecondOrderCone) = boundflag_cone
domain_type_mask(dom :: MOI.ExponentialCone)        = boundflag_cone
domain_type_mask(dom :: MOI.PowerCone)              = boundflag_cone
domain_type_mask(dom :: MOI.DualExponentialCone)    = boundflag_cone
domain_type_mask(dom :: MOI.DualPowerCone)          = boundflag_cone

domain_type_mask(dom :: MOI.PositiveSemidefiniteConeTriangle) = 0

domain_type_mask(dom :: MOI.Integer) = boundflag_int
domain_type_mask(dom :: MOI.ZeroOne) = (boundflag_int | boundflag_upper | boundflag_lower)

is_positivesemidefinite_set(dom :: PositiveSemidefiniteCone) = true
is_positivesemidefinite_set(dom) = false
is_conic_set(dom :: VectorCone) = true
is_conic_set(dom) = false

function addvarconstr(m::MosekModel, subj::Int, dom::MOI.Interval)
    putvarbound(m.task, subj, MSK_BK_RA, dom.lower, dom.upper)
end
function addvarconstr(m::MosekModel, subj::Int, dom::MOI.EqualTo)
    putvarbound(m.task, subj, MSK_BK_FX, dom.value, dom.value)
end
function addvarconstr(m::MosekModel, subj::Int, dom::MOI.Integer)
    putvartype(m.task, subj, MSK_VAR_TYPE_INT)
end
function addvarconstr(m::MosekModel, subj::Int, dom::MOI.ZeroOne)
    putvartype(m.task, subj, MSK_VAR_TYPE_INT)
    putvarbound(m.task, subj, MSK_BK_RA, 0.0, 1.0)
end

function addvarconstr(m::MosekModel, subj::Int, dom::MOI.LessThan)
    bk, lo, up = getvarbound(m.task, subj)
    bk = if (bk == MSK_BK_FR) MSK_BK_UP else MSK_BK_RA end
    putvarbound(m.task, Int32(subj), bk, lo, dom.upper)
end

function addvarconstr(m::MosekModel, subj::Int, dom::MOI.GreaterThan)
    bk, lo, up = getvarbound(m.task, subj)
    bk = if (bk == MSK_BK_FR) MSK_BK_LO else MSK_BK_RA end
    putvarbound(m.task, Int32(subj), bk, dom.lower, up)
end

addvarconstr(m :: MosekModel, subj::Vector{Int}, dom :: MOI.Reals) = nothing
function addvarconstr(m :: MosekModel, subj::Vector{Int}, dom :: MOI.Zeros)
    bnd = zeros(Float64, length(subj))
    putvarboundlist(m.task, subj, fill(MSK_BK_FX, length(subj)), bnd, bnd)
end
function addvarconstr(m :: MosekModel, subj::Vector{Int}, dom :: MOI.Nonnegatives)
    bkx = Vector{Boundkey}(undef, length(subj))
    blx = zeros(Float64, length(subj))
    bux = zeros(Float64, length(subj))
    for (i,j) in enumerate(subj)
        bk,lo,up = getvarbound(m.task, j)
        bkx[i] = if (bk == MSK_BK_FR) MSK_BK_RA else MSK_BK_LO end
        bux[i] = up
    end
    putvarboundlist(m.task, subj, bkx, blx, bux)
end

function addvarconstr(m :: MosekModel, subj::Vector{Int}, dom :: MOI.Nonpositives)
    bkx = Vector{Boundkey}(undef, length(subj))
    blx = zeros(Float64, length(subj))
    bux = zeros(Float64, length(subj))
    for (i,j) in enumerate(subj)
        bk,lo,up = getvarbound(m.task, j)
        bkx[i] = if (bk == MSK_BK_FR) MSK_BK_RA else MSK_BK_UP end
        blx[i] = lo
    end
    putvarboundlist(m.task, subj, bkx, blx, bux)
end

abstractset2ct(dom::MOI.ExponentialCone)        = MSK_CT_PEXP
abstractset2ct(dom::MOI.DualExponentialCone)    = MSK_CT_DEXP
abstractset2ct(dom::MOI.PowerCone)              = MSK_CT_PPOW
abstractset2ct(dom::MOI.DualPowerCone)          = MSK_CT_DPOW
abstractset2ct(dom::MOI.SecondOrderCone)        = MSK_CT_QUAD
abstractset2ct(dom::MOI.RotatedSecondOrderCone) = MSK_CT_RQUAD

function MOI.add_constraint(
    m   :: MosekModel,
    xs  :: MOI.SingleVariable,
    dom :: D) where {D <: MOI.AbstractScalarSet}

    subj = getindex(m.x_block, ref2id(xs.variable))

    mask = domain_type_mask(dom)
    if mask & m.x_boundflags[subj] != 0
        error("Cannot put multiple bound sets of the same type on a variable")
    end

    xcid = allocatevarconstraints(m, 1)

    xc_sub = getindex(m.xc_block,xcid)

    m.xc_bounds[xcid]  = mask
    m.xc_idxs[xc_sub] = subj

    addvarconstr(m, subj, dom)

    m.x_boundflags[subj] |= mask

    conref = MOI.ConstraintIndex{MOI.SingleVariable,D}(UInt64(xcid) << 1)

    select(m.constrmap,MOI.SingleVariable,D)[conref.value] = xcid
    conref
end

function MOI.add_constraint(m   :: MosekModel,
                            xs  :: MOI.VectorOfVariables,
                            dom :: D) where { D <: PositiveSemidefiniteCone }
    N = dom.side_dimension
    vars = sympackedUtoL(xs.variables, N)
    subj = Vector{Int}(undef, length(vars))
    for i in 1:length(subj)
        getindexes(m.x_block, ref2id(vars[i]), subj, i)
    end

    mask = domain_type_mask(dom)
    if any(mask .& m.x_boundflags[subj[1]] .> 0)
        error("Cannot put multiple bound sets of the same type to a variable")
    end

    if N < 2
        error("Invalid dimension for semidefinite constraint")
    end

    NN = sympackedlen(N)

    if length(subj) != NN
        error("Mismatching variable length for semidefinite constraint")
    end

    id = allocateconstraints(m,NN)

    subi = getindexes(m.c_block,id)

    subii32 = convert(Vector{Int32},subi)
    putaijlist(m.task,
               subii32,
               convert(Vector{Int32},subj),
               ones(Float64,NN))

    addbound!(m,id,subi,zeros(Float64,NN),dom)

    #putconboundlist(m.task,subii32,fill(MSK_BK_FX,NN),zeros(Float64,NN),zeros(Float64,NN))
    #idx = 1
    #for j in 1:N
    #    for i in 1:N
    #        symmatidx = appendsparsesymmat(m.task,N,Int32[i],Int32[j],Float64[-1.0])
    #        putbaraij(m.task,subii32[idx],barvaridx,[symmatidx],Float64[1.0])
    #        idx += 1
    #    end
    #end

    # HACK: We need to return a negative to indicate that this is
    # not, in fact, a real variable constraint, but rather a real
    # constraint, but to resolve return value at compile time we
    # need to disguise it as a variable constraint.
    #id2cref{MOI.VectorOfVariables,D}(-id)

    conref = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.PositiveSemidefiniteConeTriangle}(UInt64((id << 1) | 1))
    select(m.constrmap,MOI.VectorOfVariables,MOI.PositiveSemidefiniteConeTriangle)[conref.value] = id

    conref
end

coneparfromset(dom :: MOI.PowerCone{Float64})     = dom.exponent
coneparfromset(dom :: MOI.DualPowerCone{Float64}) = dom.exponent
coneparfromset(dom :: C) where C <: MOI.AbstractSet = 0.0

function aux_setvardom(m::MosekModel,
                       xcid::Int,
                       subj::Vector{Int},
                       dom::D) where { D <: VectorCone }
    appendcone(m.task,abstractset2ct(dom),  coneparfromset(dom), subj)
    coneidx = getnumcone(m.task)
    m.conecounter += 1
    putconename(m.task,coneidx,"$(m.conecounter)")
    m.xc_coneid[xcid] = m.conecounter
end
function aux_setvardom(m::MosekModel, xcid::Int, subj::Vector{Int},
                       dom :: D) where {D <: MOI.AbstractSet}
    addvarconstr(m, subj, dom)
end

function MOI.add_constraint(m :: MosekModel, xs :: MOI.VectorOfVariables, dom :: D) where {D <: MOI.AbstractSet}
    subj = Vector{Int}(undef,length(xs.variables))
    for i in 1:length(subj)
        getindexes(m.x_block, ref2id(xs.variables[i]),subj,i)
    end

    mask = domain_type_mask(dom)
    if any(mask .& m.x_boundflags[subj] .> 0)
        error("Cannot multiple bound sets of the same type to a variable")
    end

    N = MOI.dimension(dom)
    xcid = allocatevarconstraints(m,N)
    xc_sub = getindexes(m.xc_block,xcid)

    m.xc_bounds[xcid] = mask
    m.xc_idxs[xc_sub] = subj

    aux_setvardom(m, xcid, subj, dom)

    m.x_boundflags[subj] .|= mask

    conref = MOI.ConstraintIndex{MOI.VectorOfVariables,D}(UInt64(xcid) << 1)
    select(m.constrmap,MOI.VectorOfVariables,D)[conref.value] = xcid
    conref
end

################################################################################
################################################################################

# Put the linear left-hand side
function addlhsblock!(m        :: MosekModel,
                      conid    :: Int,
                      conidxs  :: Vector{Int},
                      terms    :: Vector{MOI.ScalarAffineTerm{Float64}})
    consubi = getindexes(m.c_block,conid)
    subj = Vector{Int}(undef,length(terms))
    for i in 1:length(subj)
        getindexes(m.x_block,ref2id(terms[i].variable_index),subj,i)
    end

    N = length(consubi)
    nnz = length(terms)

    msk_subi = convert(Vector{Int32},consubi)

    msk_rowptr = zeros(Int64, N+1)
    for i in conidxs
        msk_rowptr[i+1] += 1
    end
    msk_rowptr[1] = 1
    for i in 2:N+1
        msk_rowptr[i] += msk_rowptr[i-1]
    end

    msk_subj = Array{Int32}(undef,nnz)
    msk_cof  = Array{Float64}(undef,nnz)

    # sort by row
    for i in 1:nnz
        msk_subj[msk_rowptr[conidxs[i]]] = subj[i]
        msk_cof[msk_rowptr[conidxs[i]]]  = terms[i].coefficient
        msk_rowptr[conidxs[i]] += 1
    end

    for i in N:-1:1
        msk_rowptr[i+1] = msk_rowptr[i]
    end
    msk_rowptr[1] = 1

    putarowlist(m.task, msk_subi, msk_rowptr[1:N], msk_rowptr[2:N+1], msk_subj, msk_cof)
end

addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.Reals)                = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FR,MOI.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.Zeros)                = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,MOI.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.Nonnegatives)         = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_LO,MOI.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.Nonpositives)         = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_UP,MOI.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.GreaterThan{Float64}) = putconbound(m.task,Int32(conidxs[1]),MSK_BK_LO,dom.lower-constant[1],dom.lower-constant[1])
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.LessThan{Float64})    = putconbound(m.task,Int32(conidxs[1]),MSK_BK_UP,dom.upper-constant[1],dom.upper-constant[1])
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.EqualTo{Float64})     = putconbound(m.task,Int32(conidxs[1]),MSK_BK_FX,dom.value-constant[1],dom.value-constant[1])
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.Interval{Float64})    = putconbound(m.task,Int32(conidxs[1]),MSK_BK_RA,dom.lower-constant[1],dom.upper-constant[1])

# need a special case for this since MOI's variable order in PEXP cone is the reverse if Mosek's
function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: D) where { D <: Union{MOI.ExponentialCone, MOI.DualExponentialCone} }
    N = MOI.dimension(dom)
    nalloc = ensurefree(m.x_block,N)

    varid = newblock(m.x_block,N)
    numvar = getnumvar(m.task)
    if nalloc > 0
        appendvars(m.task, length(m.x_block) - numvar)
        append!(m.x_boundflags, zeros(Int,length(m.x_block) - numvar))
        append!(m.x_numxc, zeros(Int,length(m.x_block) - numvar))
    end
    subj = getindexes(m.x_block,varid)
    subj = [ subj[3], subj[2], subj[1] ]

    putaijlist(m.task,conidxs,subj,-ones(Float64,N))
    putvarboundlist(m.task,subj,fill(MSK_BK_FR,N),zeros(Float64,N),zeros(Float64,N))
    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,N),-constant,-constant)

    m.c_block_slack[conid] = varid

    appendcone(m.task,abstractset2ct(dom),coneparfromset(dom),subj)
    coneidx = getnumcone(m.task)
    m.conecounter += 1
    #putconename(m.task,coneidx,"$(m.conecounter)")
    m.c_coneid[conid] = m.conecounter
end

function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: D) where { D <: VectorCone }
    N = MOI.dimension(dom)
    nalloc = ensurefree(m.x_block,N)

    varid = newblock(m.x_block,N)
    numvar = getnumvar(m.task)
    if nalloc > 0
        appendvars(m.task, length(m.x_block) - numvar)
        append!(m.x_boundflags, zeros(Int,length(m.x_block) - numvar))
        append!(m.x_numxc, zeros(Int,length(m.x_block) - numvar))
    end
    subj = getindexes(m.x_block,varid)

    putaijlist(m.task,conidxs,subj,-ones(Float64,N))
    putvarboundlist(m.task,subj,fill(MSK_BK_FR,N),zeros(Float64,N),zeros(Float64,N))
    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,N),-constant,-constant)

    m.c_block_slack[conid] = varid

    appendcone(m.task,abstractset2ct(dom),coneparfromset(dom),subj)
    coneidx = getnumcone(m.task)
    m.conecounter += 1
    #putconename(m.task,coneidx,"$(m.conecounter)")
    m.c_coneid[conid] = m.conecounter
end

function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.PositiveSemidefiniteConeTriangle)
    dim = dom.side_dimension
    appendbarvars(m.task,Int32[dim])
    barvaridx = getnumbarvar(m.task)

    idx = 1
    for j in 1:dim
        for i in j:dim
            matrixid = appendsparsesymmat(m.task,Int32(dim), Int32[i], Int32[j], Float64[1.0])
            if i == j
                putbaraij(m.task, conidxs[idx], barvaridx, Int[matrixid], Float64[-1.0])
            else
                putbaraij(m.task, conidxs[idx], barvaridx, Int[matrixid], Float64[-0.5])
            end
            #putconname(m.task,Int32(conidxs[idx]),"bar_slack[$i,$j]")
            idx += 1
        end
    end

    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,length(constant)),-constant,-constant)

    m.c_block_slack[conid] = -barvaridx
end

#function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MOI.PositiveSemidefiniteConeScaled)
#    dim = MOI.dimension(dom)
#    #dim = floor(Int,-.5 + sqrt(.25+2*MOI.dimension(dom)))
#    n = MOI.dimension(dom)
#
#    appendbarvars(m.task,Int32[dim])
#    barvaridx = getnumbarvar(m.task)
#
#    idx = 1
#    for j in 1:dim
#        for i in j:dim
#            matrixid = appendsparsesymmat(m.task,Int32(dim), Int32[i], Int32[j], Float64[1.0])
#            if i == j
#                putbaraij(m.task, conidxs[idx], barvaridx, Int[matrixid], Float64[-1.0])
#            else
#                putbaraij(m.task, conidxs[idx], barvaridx, Int[matrixid], Float64[-sqrt(2.0)/2.0])
#            end
#            idx += 1
#        end
#    end
#
#    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,length(constant)),-constant,-constant)
#
#    m.c_block_slack[conid] = -barvaridx
#end

################################################################################
##  MODIFY #####################################################################
################################################################################

#MOI.canset(m::MosekModel, ::MOI.ConstraintSet, ::Type{MOI.ConstraintIndex{F,D}}) where { F <: LinearFunction, D <: LinearDomain } = true
#MOI.canset(m::MosekModel, ::MOI.ConstraintFunction, ::Type{MOI.ConstraintIndex{F,D}}) where { F <: Union{MOI.SingleVariable,MOI.ScalarAffineFunction}, D <: ScalarLinearDomain } = true
#MOI.canmodify(m::MosekModel, ::Type{MOI.ConstraintIndex{F,D}}, ::Type{<:MOI.AbstractFunctionModification}) where { F <: AffineFunction, D <: MOI.AbstractSet } = true

chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MOI.LessThan{Float64})    = bl,dom.upper-k
chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MOI.GreaterThan{Float64}) = dom.lower-k,bu
chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MOI.EqualTo{Float64})     = dom.value-k,dom.value-k
chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MOI.Interval{Float64})    = dom.lower-k,dom.upper-k

function MOI.set(m::MosekModel,
                 ::MOI.ConstraintSet,
                 xcref::MOI.ConstraintIndex{F,D},
                 dom::D) where { F    <: MOI.SingleVariable,
                                 D    <: ScalarLinearDomain }
    xcid = ref2id(xcref)
    j = m.xc_idxs[getindex(m.xc_block,xcid)]
    bk,bl,bu = getvarbound(m.task,j)
    bl,bu = chgbound(bl,bu,0.0,dom)
    putvarbound(m.task,j,bk,bl,bu)
end

function MOI.set(m::MosekModel,
                 ::MOI.ConstraintSet,
                 cref::MOI.ConstraintIndex{F,D},
                 dom::D) where { F    <: MOI.ScalarAffineFunction,
                                 D    <: ScalarLinearDomain }
    cid = ref2id(cref)
    i = getindex(m.c_block,cid) # since we are in a scalar domain
    bk,bl,bu = getconbound(m.task,i)
    bl,bu = chgbound(bl,bu,0.0,dom)
    putconbound(m.task,i,bk,bl,bu)
end

function set_internal_name(m::MosekModel,c::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},D},name::AbstractString) where {D}
    cid = ref2id(c)
    for i in getindexes(m.c_block, cid)
        putconname(m.task,i,name)
    end
end
function set_internal_name(m::MosekModel,c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D},name::AbstractString) where {D}
    cid = ref2id(c)
    i = getindex(m.c_block, cid)
    putconname(m.task,i,name)
end
function set_internal_name(m::MosekModel, c::MOI.ConstraintIndex{F,D}, name::AbstractString) where {F,D} end

function MOI.set(m    ::MosekModel,
                 ::MOI.ConstraintName,
                 c    ::MOI.ConstraintIndex{F,D},
                 name ::AbstractString) where {F,D}#{F<:MOI.AbstractFunction,D<:AbstractSet}
    if ! haskey(m.constrnames, name)
        m.constrnames[name] = MOI.ConstraintIndex[]
    end
    push!(m.constrnames[name], c)
    set_internal_name(m,c,name)
end

function MOI.modify(m   ::MosekModel,
                    c   ::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D},
                    func::MOI.ScalarConstantChange{Float64}) where {D <: MOI.AbstractSet}

    cid = ref2id(c)

    i = getindex(m.c_block, cid)
    bk,bl,bu = getconbound(m.task,i)
    bl += m.c_constant[i] - func.new_constant
    bu += m.c_constant[i] - func.new_constant
    m.c_constant[i] = func.new_constant
    putconbound(m.task,i,bk,bl,bu)
end

function MOI.modify(m   ::MosekModel,
                    c   ::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D},
                    func::MOI.ScalarCoefficientChange{Float64}) where {D <: MOI.AbstractSet}
    cid = ref2id(c)
    xid = ref2id(func.variable)

    i = getindex(m.c_block,cid)
    j = getindex(m.x_block,xid)

    putaij(m.task,i,j,func.new_coefficient)
end

function MOI.modify(m::MosekModel,
                    c::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},D},
                    func::MOI.VectorConstantChange{Float64}) where {D <: MOI.AbstractSet}
    cid = ref2id(c)
    @assert(cid > 0)

    subi = getindexes(m.c_block, cid)
    bk = Vector{Int32}(undef,length(subi))
    bl = Vector{Float64}(undef,length(subi))
    bu = Vector{Float64}(undef,length(subi))

    bk,bl,bu = getconboundlist(m.task,convert(Vector{Int32},subi))
    bl += m.c_constant[subi] - func.new_constant
    bu += m.c_constant[subi] - func.new_constant
    m.c_constant[subi] = func.new_constant

    putconboundlist(m.task,convert(Vector{Int32},subi),bk,bl,bu)
end

function MOI.modify(m::MosekModel,
                    c::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},D},
                    func::MOI.MultirowChange{Float64}) where {D <: MOI.AbstractSet}
    cid = ref2id(c)
    @assert(cid > 0)

    subi = getindexes(m.c_block, cid)[getindex.(func.new_coefficients, 1)]
    xid = ref2id(func.variable)
    j = getindex(m.x_block,xid)

    putaijlist(m.task,convert(Vector{Int32},subi),fill(j,length(subi)),getindex.(func.new_coefficients,2))
end

#MOI.cantransform(m::MosekModel, c::MOI.ConstraintIndex{F,D1}, newdom::D2) where {F <: VariableFunction, D1, D2 } = false
#MOI.cantransform(m::MosekModel, c::MOI.ConstraintIndex{MOI.VectorAffineFunction,D1}, newdom::D2) where {D1 <: VectorLinearDomain, D2 <: VectorLinearDomain} = false
#function MOI.cantransform(m::MosekModel,
#                          cref::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D1},
#                          newdom::D2) where {D1 <: ScalarLinearDomain,
#                                             D2 <: ScalarLinearDomain}
#    haskey(select(m.constrmap,MOI.ScalarAffineFunction{Float64},D1),cref.value)
#end

function MOI.transform(m::MosekModel,
                       cref::MOI.ConstraintIndex{F,D},
                       newdom::D) where {F <: MOI.AbstractFunction,
                                         D <: MOI.AbstractSet}
    MOI.modify(m,cref,newdom)
    cref
end

function MOI.transform(m::MosekModel,
                       cref::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D1},
                       newdom::D2) where {D1 <: ScalarLinearDomain,
                                          D2 <: ScalarLinearDomain}
    F = MOI.ScalarAffineFunction{Float64}

    cid = ref2id(cref)

    subi = getindexes(m.c_block,cid)

    addbound!(m,cid,subi,m.c_constant[subi], newdom)

    newcref = MOI.ConstraintIndex{F,D2}(UInt64(cid) << 1)
    delete!(select(m.constrmap,F,D1), cref.value)
    select(m.constrmap, F, D2)[newcref.value] = cid
    newcref
end

function MOI.transform(m::MosekModel,
                       cref::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},D1},
                       newdom::D2) where {D1 <: VectorLinearDomain,
                                          D2 <: VectorLinearDomain}
    F = MOI.VectorAffineFunction{Float64}

    cid = ref2id(cref)

    subi = getindexes(m.c_block,cid)

    addbound!(m,cid,subi,m.c_constant[subi], newdom)

    newcref = MOI.ConstraintIndex{F,D2}(UInt64(cid) << 1)
    delete!(select(m.constrmap,F,D1), cref.value)
    select(m.constrmap,F,D2)[newcref.value] = cid
    newcref
end

#MOI.transform(m::MosekModel, c::MOI.ConstraintIndex{F,D1}, newdom::D2) where {F <: MOI.VectorAffineFunction , D1 <: VectorLinearDomain, D2 <: VectorLinearDomain} = false

################################################################################
##  DELETE #####################################################################
################################################################################

#MOI.candelete(
#    m   ::MosekModel,
#    cref::MOI.ConstraintIndex{F,D}) where {F <: Union{MOI.ScalarAffineFunction,
#                                                                       MOI.VectorAffineFunction,
#                                                                       MOI.SingleVariable,
#                                                                       MOI.VectorOfVariables},
#                                                            D <: Union{MOI.LessThan,
#                                                                       MOI.GreaterThan,
#                                                                       MOI.EqualTo,
#                                                                       MOI.Interval,
#                                                                       MOI.Zeros,
#                                                                       MOI.Nonpositives,
#                                                                       MOI.Nonnegatives,
#                                                                       MOI.Reals}} = MOI.isvalid(m,cref)
#
#MOI.candelete(
#    m   ::MosekModel,
#    cref::MOI.ConstraintIndex{F,D}) where {F <: MOI.AbstractFunction,
#                                                            D <: Union{MOI.SecondOrderCone,
#                                                                       MOI.RotatedSecondOrderCone,
#                                                                       MOI.ExponentialCone,
#                                                                       MOI.DualExponentialCone,
#                                                                       MOI.PowerCone,
#                                                                       MOI.DualPowerCone}} = false
#
#MOI.candelete(
#    m   ::MosekModel,
#    cref::MOI.ConstraintIndex{F,D}) where {F <: Union{MOI.SingleVariable,
#                                                                       MOI.VectorOfVariables},
#                                                            D <: Union{MOI.SecondOrderCone,
#                                                                       MOI.RotatedSecondOrderCone,
#                                                                       MOI.ExponentialCone,
#                                                                       MOI.DualExponentialCone,
#                                                                       MOI.PowerCone,
#                                                                       MOI.DualPowerCone}} = true
#

function MOI.delete(
    m::MosekModel,
    cref::MOI.ConstraintIndex{F,D}) where {F <: Union{MOI.ScalarAffineFunction,
                                                                       MOI.VectorAffineFunction},
                                                            D <: Union{MOI.LessThan,
                                                                       MOI.GreaterThan,
                                                                       MOI.EqualTo,
                                                                       MOI.Interval,
                                                                       MOI.Zeros,
                                                                       MOI.Nonpositives,
                                                                       MOI.Nonnegatives,
                                                                       MOI.Reals}}

    delete!(select(m.constrmap, F, D), cref.value)

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)

    n = length(subi)
    subi_i32 = convert(Vector{Int32},subi)
    ptr = fill(Int64(0),n)
    putarowlist(m.task,subi_i32,ptr,ptr,Int32[],Float64[])
    b = fill(0.0,n)
    putconboundlist(m.task,subi_i32,fill(MSK_BK_FX,n),b,b)

    m.c_constant[subi] .= 0.0
    deleteblock(m.c_block,cid)
end

function MOI.delete(
    m::MosekModel,
    cref::MOI.ConstraintIndex{F,D}) where {F <: Union{MOI.VectorAffineFunction},
                                                            D <: Union{MOI.SecondOrderCone,
                                                                       MOI.RotatedSecondOrderCone,
                                                                       MOI.ExponentialCone,
                                                                       MOI.DualExponentialCone,
                                                                       MOI.PowerCone,
                                                                       MOI.DualPowerCone}}

    delete!(select(m.constrmap, F, D), cref.value)
    xcid = ref2id(cref)
    sub = getindexes(m.xc_block,xcid)

    subj = [ getindex(m.x_block,i) for i in sub ]
    N = length(subj)
    m.x_boundflags[subj] .&= ~m.xc_bounds[xcid]
    asgn,coneidx = getconenameindex(m.task,"$(m.xc_coneid[xcid])")
    m.xc_coneid[xcid] = 0
    removecone(m.task,coneidx)

    m.x_numxc[subj]  -= 1
    m.xc_idxs[sub]    = 0
    m.xc_bounds[xcid] = 0

    deleteblock(m.xc_block,xcid)
end

function MOI.delete(
    m::MosekModel,
    cref::MOI.ConstraintIndex{F,D}) where {F <: Union{MOI.SingleVariable,
                                                      MOI.VectorOfVariables},
                                           D <: Union{ScalarLinearDomain,
                                                      VectorLinearDomain,
                                                      ScalarIntegerDomain}}
    delete!(select(m.constrmap, F, D), cref.value)

    xcid = ref2id(cref)
    sub = getindexes(m.xc_block,xcid)

    subj = [ getindex(m.x_block,i) for i in sub ]
    N = length(subj)

    m.x_boundflags[subj] .&= ~m.xc_bounds[xcid]
    if m.xc_bounds[xcid] & boundflag_int != 0
        for i in 1:length(subj)
            putvartype(m.task,subj[i],MSK_VAR_TYPE_CONT)
        end
    end

    if m.xc_bounds[xcid] & boundflag_lower != 0 && m.xc_bounds[xcid] & boundflag_upper != 0
        bnd = fill(0.0, length(N))
        putvarboundlist(m.task,convert(Vector{Int32},subj),fill(MSK_BK_FR,N),bnd,bnd)
    elseif m.xc_bounds[xcid] & boundflag_lower != 0
        bnd = fill(0.0, length(N))
        bk,bl,bu = getvarboundlist(m.task, convert(Vector{Int32},subj))

        for i in 1:N
            if MSK_BK_RA == bk[i] || MSK_BK_FX == bk[i]
                bk[i] = MSK_BK_UP
            else
                bk[i] = MSK_BK_FR
            end
        end

        putvarboundlist(m.task,convert(Vector{Int32},subj),bk,bl,bu)
    elseif m.xc_bounds[xcid] & boundflag_upper != 0
        bnd = fill(0.0, length(N))
        bk,bl,bu = getvarboundlist(m.task, subj)

        for i in 1:N
            if MSK_BK_RA == bk[i] || MSK_BK_FX == bk[i]
                bk[i] = MSK_BK_LO
            else
                bk[i] = MSK_BK_FR
            end
        end

        putvarboundlist(m.task,convert(Vector{Int32},subj),bk,bl,bu)
    else
        @assert(false)
        # should not happen
    end

    m.x_numxc[subj] .-= 1
    m.xc_idxs[sub] .= 0
    m.xc_bounds[xcid] = 0

    deleteblock(m.xc_block, xcid)
end

################################################################################
################################################################################
################################################################################

function allocateconstraints(m :: MosekModel,
                             N :: Int)
    numcon = getnumcon(m.task)
    alloced = ensurefree(m.c_block,N)
    id = newblock(m.c_block,N)

    M = numblocks(m.c_block) - length(m.c_block_slack)
    if alloced > 0
        appendcons(m.task, alloced)
        append!(m.c_constant, zeros(Float64,alloced))
    end
    if M > 0
        append!(m.c_block_slack, zeros(Float64,M))
        append!(m.c_coneid, zeros(Float64,M))
    end
    id
end

function allocatevarconstraints(m :: MosekModel,
                                N :: Int)
    nalloc = ensurefree(m.xc_block,N)
    id = newblock(m.xc_block,N)

    M = numblocks(m.xc_block) - length(m.xc_bounds)
    if M > 0
        append!(m.xc_bounds,zeros(Float64,M))
        append!(m.xc_coneid,zeros(Float64,M))
    end
    if nalloc > 0
        append!(m.xc_idxs, zeros(Float64,nalloc))
    end

    return id
end

function allocatevariable(m :: MosekModel, N :: Int)
    @assert(length(m.x_boundflags) == length(m.x_block))
    numvar = getnumvar(m.task)
    alloced = ensurefree(m.x_block, N)
    if alloced > 0
        appendvars(m.task, length(m.x_block) - numvar)
        append!(m.x_boundflags, zeros(Int,length(m.x_block) - numvar))
        append!(m.x_numxc, zeros(Int,length(m.x_block) - numvar))
    end
    return newblock(m.x_block, N)
end

function MOI.is_valid(model::MosekModel,
                      ref::MOI.ConstraintIndex{F, D}) where {F, D}
    return haskey(select(model.constrmap, F, D), ref.value)
end
function MOI.is_valid(model::MosekModel, ref::MOI.VariableIndex)
    return allocated(model.x_block, ref2id(ref))
end

function getvarboundlist(t::Mosek.Task, subj :: Vector{Int32})
    n = length(subj)
    bk = Vector{Boundkey}(undef,n)
    bl = Vector{Float64}(undef,n)
    bu = Vector{Float64}(undef,n)
    for i in 1:n
        bki,bli,bui = getvarbound(t,subj[i])
        bk[i] = bki
        bl[i] = bli
        bu[i] = bui
    end
    bk,bl,bu
end

function getconboundlist(t::Mosek.Task, subj :: Vector{Int32})
    n = length(subj)
    bk = Vector{Boundkey}(undef,n)
    bl = Vector{Float64}(undef,n)
    bu = Vector{Float64}(undef,n)
    for i in 1:n
        bki,bli,bui = getconbound(t,subj[i])
        bk[i] = bki
        bl[i] = bli
        bu[i] = bui
    end
    bk,bl,bu
end


text:
int_rules_1_2_1_9 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.1.9*P(x)*(d+e*x)^m*(a+b*x+c*x^2)^p =#
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * ~x) ^ (~m + 1) * PolynomialQuotient(~Pq, ~d + ~e * ~x, ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ(PolynomialRemainder(~Pq, ~d + ~e * ~x, ~x), 0))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * ~x) ^ (~m + 1) * PolynomialQuotient(~Pq, ~d + ~e * ~x, ~x) * (~a + ~c * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ(PolynomialRemainder(~Pq, ~d + ~e * ~x, ~x), 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~P2 * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => With([f = Coeff(~P2, ~x, 0), g = Coeff(~P2, ~x, 1), h = Coeff(~P2, ~x, 2)], (h * (~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * ~e * (~m + 2 * ~p + 3)) <-- EqQ((~b * ~e * h * (~m + ~p + 2) + 2 * ~c * ~d * h * (~p + 1)) - ~c * ~e * g * (~m + 2 * ~p + 3), 0) && EqQ((~b * ~d * h * (~p + 1) + ~a * ~e * h * (~m + 1)) - ~c * ~e * f * (~m + 2 * ~p + 3), 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~P2, ~x, 2) && NeQ(~m + 2 * ~p + 3, 0))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~P2 * (~a + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => With([f = Coeff(~P2, ~x, 0), g = Coeff(~P2, ~x, 1), h = Coeff(~P2, ~x, 2)], (h * (~d + ~e * ~x) ^ (~m + 1) * (~a + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * ~e * (~m + 2 * ~p + 3)) <-- EqQ(2 * ~d * h * (~p + 1) - ~e * g * (~m + 2 * ~p + 3), 0) && EqQ(~a * h * (~m + 1) - ~c * f * (~m + 2 * ~p + 3), 0)) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~P2, ~x, 2) && NeQ(~m + 2 * ~p + 3, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * ~Pq * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (PolyQ(~Pq, ~x) && IGtQ(~p, -2))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * ~x) ^ ~m * ~Pq * (~a + ~c * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~m], ~x) && (PolyQ(~Pq, ~x) && IGtQ(~p, -2))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * ~x) ^ (2 * FracPart(~p)))) * Antiderivative((~d + ~e * ~x) ^ ~m * ~Pq * (~b + 2 * ~c * ~x) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * ~Pq * (~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => ~e * Antiderivative((~e * ~x) ^ (~m - 1) * PolynomialQuotient(~Pq, ~b + ~c * ~x, ~x) * (~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1), ~x) <-- FreeQ([~b, ~c, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ(PolynomialRemainder(~Pq, ~b + ~c * ~x, ~x), 0))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => ~d * ~e * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * PolynomialQuotient(~Pq, ~a * ~e + ~c * ~d * ~x, ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && EqQ(PolynomialRemainder(~Pq, ~a * ~e + ~c * ~d * ~x, ~x), 0))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => ~d * ~e * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * PolynomialQuotient(~Pq, ~a * ~e + ~c * ~d * ~x, ~x) * (~a + ~c * (~x) ^ 2) ^ (~p + 1), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && EqQ(PolynomialRemainder(~Pq, ~a * ~e + ~c * ~d * ~x, ~x), 0)))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a * ~e + ~c * ~d * ~x, ~x), f = PolynomialRemainder(~Pq, ~a * ~e + ~c * ~d * ~x, ~x)], (f * (2 * ~c * ~d - ~b * ~e) * (~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / (~e * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum(~d * ~e * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * Q - f * (2 * ~c * ~d - ~b * ~e) * (~m + 2 * ~p + 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (ILtQ(~p + 1 / 2, 0) && GtQ(~m, 0)))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a * ~e + ~c * ~d * ~x, ~x), f = PolynomialRemainder(~Pq, ~a * ~e + ~c * ~d * ~x, ~x)], (-(~d) * f * (~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ (~p + 1)) / (2 * ~a * ~e * (~p + 1)) + (~d / (2 * ~a * (~p + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum(2 * ~a * ~e * (~p + 1) * Q + f * (~m + 2 * ~p + 2), ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (ILtQ(~p + 1 / 2, 0) && GtQ(~m, 0))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, (~d + ~e * ~x) ^ ~m * ~Pq, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (EqQ(~m + Expon(~Pq, ~x) + 2 * ~p + 1, 0) && ILtQ(~m, 0)))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~m * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~c * (~x) ^ 2) ^ ~p, (~d + ~e * ~x) ^ ~m * ~Pq, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (EqQ(~m + Expon(~Pq, ~x) + 2 * ~p + 1, 0) && ILtQ(~m, 0))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x), f = Coeff(~Pq, ~x, Expon(~Pq, ~x))], (f * (~d + ~e * ~x) ^ ((~m + q) - 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * (~e) ^ (q - 1) * (~m + q + 2 * ~p + 1)) + (1 / (~c * (~e) ^ q * (~m + q + 2 * ~p + 1))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((~c * (~e) ^ q * (~m + q + 2 * ~p + 1) * ~Pq - ~c * f * (~m + q + 2 * ~p + 1) * (~d + ~e * ~x) ^ q) + ~e * f * (~m + ~p + q) * (~d + ~e * ~x) ^ (q - 2) * ((~b * ~d - 2 * ~a * ~e) + (2 * ~c * ~d - ~b * ~e) * ~x), ~x), ~x) <-- NeQ(~m + q + 2 * ~p + 1, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0)))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x), f = Coeff(~Pq, ~x, Expon(~Pq, ~x))], (f * (~d + ~e * ~x) ^ ((~m + q) - 1) * (~a + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * (~e) ^ (q - 1) * (~m + q + 2 * ~p + 1)) + (1 / (~c * (~e) ^ q * (~m + q + 2 * ~p + 1))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((~c * (~e) ^ q * (~m + q + 2 * ~p + 1) * ~Pq - ~c * f * (~m + q + 2 * ~p + 1) * (~d + ~e * ~x) ^ q) - 2 * ~e * f * (~m + ~p + q) * (~d + ~e * ~x) ^ (q - 2) * (~a * ~e - ~c * ~d * ~x), ~x), ~x) <-- NeQ(~m + q + 2 * ~p + 1, 0)) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IGtQ(~m, 0))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * ~x) ^ (~m + ~p) * (~a / ~d + (~c / ~e) * ~x) ^ ~p * ~Pq, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && IntegerQ(~p))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * ~x) ^ (~m + ~p) * (~a / ~d + (~c / ~e) * ~x) ^ ~p * ~Pq, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~m], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((~d + ~e * ~x) ^ FracPart(~p) * (~a / ~d + (~c * ~x) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * ~x) ^ (~m + ~p) * (~a / ~d + (~c / ~e) * ~x) ^ ~p * ~Pq, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~p)) && Not(IGtQ(~m, 0))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => ((~a + ~c * (~x) ^ 2) ^ FracPart(~p) / ((~d + ~e * ~x) ^ FracPart(~p) * (~a / ~d + (~c * ~x) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * ~x) ^ (~m + ~p) * (~a / ~d + (~c / ~e) * ~x) ^ ~p * ~Pq, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~p)) && Not(IGtQ(~m, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 1)], ((~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ((f * ~b - 2 * ~a * g) + (2 * ~c * f - ~b * g) * ~x)) / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum((((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * (~d + ~e * ~x) * Q + g * (2 * ~a * ~e * ~m + ~b * ~d * (2 * ~p + 3))) - f * (~b * ~e * ~m + 2 * ~c * ~d * (2 * ~p + 3))) - ~e * (2 * ~c * f - ~b * g) * (~m + 2 * ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && (GtQ(~m, 0) && ((IntegerQ(~p) || (Not(IntegerQ(~m)) || Not(RationalQ(~a, ~b, ~c, ~d, ~e)))) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~b, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0))))))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~c * (~x) ^ 2, ~x), ~x, 1)], ((~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~a * g - ~c * f * ~x)) / (2 * ~a * ~c * (~p + 1)) + (1 / (2 * ~a * ~c * (~p + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum((2 * ~a * ~c * (~p + 1) * (~d + ~e * ~x) * Q - ~a * ~e * g * ~m) + ~c * ~d * f * (2 * ~p + 3) + ~c * ~e * f * (~m + 2 * ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && (GtQ(~m, 0) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0))))))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient((~d + ~e * ~x) ^ ~m * ~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder((~d + ~e * ~x) ^ ~m * ~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder((~d + ~e * ~x) ^ ~m * ~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 1)], (((~b * f - 2 * ~a * g) + (2 * ~c * f - ~b * g) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * (~d + ~e * ~x) ^ -(~m) * Q - (2 * ~p + 3) * (2 * ~c * f - ~b * g) * (~d + ~e * ~x) ^ -(~m), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && ILtQ(~m, 0)))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient((~d + ~e * ~x) ^ ~m * ~Pq, ~a + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder((~d + ~e * ~x) ^ ~m * ~Pq, ~a + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder((~d + ~e * ~x) ^ ~m * ~Pq, ~a + ~c * (~x) ^ 2, ~x), ~x, 1)], ((~a * g - ~c * f * ~x) * (~a + ~c * (~x) ^ 2) ^ (~p + 1)) / (2 * ~a * ~c * (~p + 1)) + (1 / (2 * ~a * ~c * (~p + 1))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum(2 * ~a * ~c * (~p + 1) * (~d + ~e * ~x) ^ -(~m) * Q + ~c * f * (2 * ~p + 3) * (~d + ~e * ~x) ^ -(~m), ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Pq, ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && ILtQ(~m, 0))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~b * ~x + ~c * (~x) ^ 2, ~x), ~x, 1)], ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ((f * ((~b * ~c * ~d - (~b) ^ 2 * ~e) + 2 * ~a * ~c * ~e) - ~a * g * (2 * ~c * ~d - ~b * ~e)) + ~c * (f * (2 * ~c * ~d - ~b * ~e) - g * (~b * ~d - 2 * ~a * ~e)) * ~x)) / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) + (1 / ((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum((((~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2) * Q + f * (((~b * ~c * ~d * ~e * ((2 * ~p - ~m) + 2) + (~b) ^ 2 * (~e) ^ 2 * (~p + ~m + 2)) - 2 * (~c) ^ 2 * (~d) ^ 2 * (2 * ~p + 3)) - 2 * ~a * ~c * (~e) ^ 2 * (~m + 2 * ~p + 3))) - g * (~a * ~e * ((~b * ~e - 2 * ~c * ~d * ~m) + ~b * ~e * ~m) - ~b * ~d * (((3 * ~c * ~d - ~b * ~e) + 2 * ~c * ~d * ~p) - ~b * ~e * ~p))) + ~c * ~e * (g * (~b * ~d - 2 * ~a * ~e) - f * (2 * ~c * ~d - ~b * ~e)) * (~m + 2 * ~p + 4) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~b, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0))))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~c * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~c * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~c * (~x) ^ 2, ~x), ~x, 1)], (-((~d + ~e * ~x) ^ (~m + 1)) * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~a * (~e * f - ~d * g) + (~c * ~d * f + ~a * ~e * g) * ~x)) / (2 * ~a * (~p + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) + (1 / (2 * ~a * (~p + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * ExpandToSum(((2 * ~a * (~p + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2) * Q + ~c * (~d) ^ 2 * f * (2 * ~p + 3)) - ~a * ~e * (~d * g * ~m - ~e * f * (~m + 2 * ~p + 3))) + ~e * (~c * ~d * f + ~a * ~e * g) * (~m + 2 * ~p + 4) * ~x, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~m], ~x) && (PolyQ(~Pq, ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0)))))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~d + ~e * ~x, ~x), R = PolynomialRemainder(~Pq, ~d + ~e * ~x, ~x)], (~e * R * (~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / ((~m + 1) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) + (1 / ((~m + 1) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((((~m + 1) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2) * Q + ~c * ~d * R * (~m + 1)) - ~b * ~e * R * (~m + ~p + 2)) - ~c * ~e * R * (~m + 2 * ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && LtQ(~m, -1))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~m * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~d + ~e * ~x, ~x), R = PolynomialRemainder(~Pq, ~d + ~e * ~x, ~x)], (~e * R * (~d + ~e * ~x) ^ (~m + 1) * (~a + ~c * (~x) ^ 2) ^ (~p + 1)) / ((~m + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) + (1 / ((~m + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) * (~a + ~c * (~x) ^ 2) ^ ~p * ExpandToSum(((~m + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2) * Q + ~c * ~d * R * (~m + 1)) - ~c * ~e * R * (~m + 2 * ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && LtQ(~m, -1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x), k], Antiderivative((~x) ^ ~m * Sum(Coeff(~Pq, ~x, 2k) * (~x) ^ (2k), [k, 0, q / 2]) * (~a + ~b * (~x) ^ 2) ^ ~p, ~x) + Antiderivative((~x) ^ (~m + 1) * Sum(Coeff(~Pq, ~x, 2k + 1) * (~x) ^ (2k), [k, 0, (q - 1) / 2]) * (~a + ~b * (~x) ^ 2) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~p], ~x) && (PolyQ(~Pq, ~x) && (Not(PolyQ(~Pq, (~x) ^ 2)) && (IGtQ(~m, -2) && Not(IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x), f = Coeff(~Pq, ~x, Expon(~Pq, ~x))], (f * (~d + ~e * ~x) ^ ((~m + q) - 1) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * (~e) ^ (q - 1) * (~m + q + 2 * ~p + 1)) + (1 / (~c * (~e) ^ q * (~m + q + 2 * ~p + 1))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((~c * (~e) ^ q * (~m + q + 2 * ~p + 1) * ~Pq - ~c * f * (~m + q + 2 * ~p + 1) * (~d + ~e * ~x) ^ q) - f * (~d + ~e * ~x) ^ (q - 2) * (((~b * ~d * ~e * (~p + 1) + ~a * (~e) ^ 2 * ((~m + q) - 1)) - ~c * (~d) ^ 2 * (~m + q + 2 * ~p + 1)) - ~e * (2 * ~c * ~d - ~b * ~e) * (~m + q + ~p) * ~x), ~x), ~x) <-- GtQ(q, 1) && NeQ(~m + q + 2 * ~p + 1, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~b, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0)))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x), f = Coeff(~Pq, ~x, Expon(~Pq, ~x))], (f * (~d + ~e * ~x) ^ ((~m + q) - 1) * (~a + ~c * (~x) ^ 2) ^ (~p + 1)) / (~c * (~e) ^ (q - 1) * (~m + q + 2 * ~p + 1)) + (1 / (~c * (~e) ^ q * (~m + q + 2 * ~p + 1))) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((~c * (~e) ^ q * (~m + q + 2 * ~p + 1) * ~Pq - ~c * f * (~m + q + 2 * ~p + 1) * (~d + ~e * ~x) ^ q) - f * (~d + ~e * ~x) ^ (q - 2) * ((~a * (~e) ^ 2 * ((~m + q) - 1) - ~c * (~d) ^ 2 * (~m + q + 2 * ~p + 1)) - 2 * ~c * ~d * ~e * (~m + q + ~p) * ~x), ~x), ~x) <-- GtQ(q, 1) && NeQ(~m + q + 2 * ~p + 1, 0)) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (Not(EqQ(~d, 0) && True) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0)))))))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ~Pq * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => With([q = Expon(~Pq, ~x)], (Coeff(~Pq, ~x, q) / (~e) ^ q) * Antiderivative((~d + ~e * ~x) ^ (~m + q) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x) + (1 / (~e) ^ q) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((~e) ^ q * ~Pq - Coeff(~Pq, ~x, q) * (~d + ~e * ~x) ^ q, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~b, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0)))))))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * ~Pq * (~a + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => With([q = Expon(~Pq, ~x)], (Coeff(~Pq, ~x, q) / (~e) ^ q) * Antiderivative((~d + ~e * ~x) ^ (~m + q) * (~a + ~c * (~x) ^ 2) ^ ~p, ~x) + (1 / (~e) ^ q) * Antiderivative((~d + ~e * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ ~p * ExpandToSum((~e) ^ q * ~Pq - Coeff(~Pq, ~x, q) * (~d + ~e * ~x) ^ q, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IGtQ(~m, 0) && (RationalQ(~a, ~c, ~d, ~e) && (IntegerQ(~p) || ILtQ(~p + 1 / 2, 0))))))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/1.2.1.9 P(x) (d+e x)^m (a+b x+c x^2)^p.jl


text:
function meets_compat_for_all_deps(working_directory::AbstractString, pkg, version)
    deps = Pkg.TOML.parsefile(joinpath(working_directory, uppercase(pkg[1:1]), pkg, "Deps.toml"))
    compat = Pkg.TOML.parsefile(joinpath(working_directory, uppercase(pkg[1:1]), pkg, "Compat.toml"))
    # First, we construct a Dict in which the keys are the package's
    # dependencies, and the value is always false.
    dep_has_compat_with_upper_bound = Dict{String, Bool}()
    @debug("We always have julia as a dependency")
    dep_has_compat_with_upper_bound["julia"] = false
    for version_range in keys(deps)
        if version in Pkg.Types.VersionRange(version_range)
            for name in keys(deps[version_range])
                if !is_julia_stdlib(name)
                    @debug("Found a new dependency: $(name)")
                    dep_has_compat_with_upper_bound[name] = false
                end
            end
        end
    end
    # Now, we go through all the compat entries. If a dependency has a compat
    # entry with an upper bound, we change the corresponding value in the Dict
    # to true.
    for version_range in keys(compat)
        if version in Pkg.Types.VersionRange(version_range)
            for compat_entry in compat[version_range]
                name = compat_entry[1]
                value = compat_entry[2]
                if value isa Vector
                    if !isempty(value)
                        value_ranges = Pkg.Types.VersionRange.(value)
                        each_range_has_upper_bound = _has_upper_bound.(value_ranges)
                        if all(each_range_has_upper_bound)
                            @debug("Dependency \"$(name)\" has compat entries that all have upper bounds")
                            dep_has_compat_with_upper_bound[name] = true
                        end
                    end
                else
                    value_range = Pkg.Types.VersionRange(value)
                    if _has_upper_bound(value_range)
                        @debug("Dependency \"$(name)\" has a compat entry with an upper bound")
                        dep_has_compat_with_upper_bound[name] = true
                    end
                end
            end
        end
    end
    meets_this_guideline = all(values(dep_has_compat_with_upper_bound))
    if meets_this_guideline
        return true, ""
    else
        bad_dependencies = Vector{String}()
        for name in keys(dep_has_compat_with_upper_bound)
            if !(dep_has_compat_with_upper_bound[name])
                @error("Dependency \"$(name)\" does not have a compat entry that has an upper bound")
                push!(bad_dependencies, name)
            end
        end
        sort!(bad_dependencies)
        message = string("The following dependencies do not have a compat entry that has an upper bound: ", join(bad_dependencies, ", "))
        return false, message
    end
end

function meets_patch_release_does_not_narrow_julia_compat(pkg::String,
                                                          new_version::VersionNumber;
                                                          registry_head::String,
                                                          registry_master::String)
    old_version = latest_version(pkg, registry_master)
    julia_compats_for_old_version = julia_compat(pkg, old_version, registry_master)
    julia_compats_for_new_version = julia_compat(pkg, new_version, registry_head)
    if Set(julia_compats_for_old_version) == Set(julia_compats_for_new_version)
        return true, ""
    end
    meets_this_guideline = range_did_not_narrow(julia_compats_for_old_version, julia_compats_for_new_version)
    if meets_this_guideline
        return true, ""
    else
        return false, "A patch release is not allowed to narrow the supported range of Julia versions"
    end
end

function meets_name_length(pkg)
    meets_this_guideline = length(pkg) >= 5
    if meets_this_guideline
        return true, ""
    else
        return false, "Name is not at least five characters long"
    end
end

function meets_normal_capitalization(pkg)
    meets_this_guideline = occursin(r"^[A-Z]\w*[a-z][0-9]?$", pkg)
    if meets_this_guideline
        return true, ""
    else
        return false, "Name does not meet all of the following: starts with a capital letter, ASCII alphanumerics only, ends in lowercase"
    end
end

function meets_repo_url_requirement(pkg::String; registry_head::String)
    url = Pkg.TOML.parsefile(joinpath(registry_head, pkg[1:1], pkg, "Package.toml"))["repo"]
    meets_this_guideline = url_has_correct_ending(url, pkg)
    if meets_this_guideline
        return true, ""
    else
        return false, "Repo URL does not end with /name.jl.git, where name is the package name"
    end
end

function meets_sequential_version_number(old_version::VersionNumber, new_version::VersionNumber)
    if new_version > old_version
        diff = difference(old_version, new_version)
        @debug("Difference between versions: ", old_version, new_version, diff)
        if diff == v"1.0.0"
            return true, "", :major
        elseif diff == v"0.1.0"
            return true, "", :minor
        elseif diff == v"0.0.1"
            return true, "", :patch
        else
            return false, "Does not meet sequential version number guideline", :invalid
        end
    else
        return false, "Does not meet sequential version number guideline", :invalid
    end
end

function meets_sequential_version_number(pkg::String,
                                         new_version::VersionNumber;
                                         registry_head::String,
                                         registry_master::String)
    old_version = latest_version(pkg, registry_master)
    return meets_sequential_version_number(old_version, new_version)
end

function meets_standard_initial_version_number(version)
    meets_this_guideline = version == v"0.0.1" || version == v"0.1.0" || version == v"1.0.0"
    if meets_this_guideline
        return true, ""
    else
        return false, "Version number is not 0.0.1, 0.1.0, or 1.0.0"
    end
end

function meets_version_can_be_loaded(working_directory::String,
                                     pkg::String,
                                     version::VersionNumber)
    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    code = """
        import Pkg;
        Pkg.Registry.add(Pkg.RegistrySpec(path=\"$(working_directory)\"));
        @info("Attempting to install package...");
        Pkg.add(Pkg.PackageSpec(name=\"$(pkg)\", version=\"$(string(version))\"));
        @info("Successfully installed package");
        @info("Attempting to import package");
        import $(pkg);
        @info("Successfully imported package");
        """
    # We need to be careful with what environment variables we pass to the child
    # process. For example, we don't want to pass an environment variable containing
    # our GitHub token to the child process. Because if the Julia package that we are
    # testing has malicious code in its __init__() function, it could try to steal
    # our token. So we only pass three environment variables:
    # 1. PATH. If we don't pass PATH, things break. And PATH should not contain any
    #    sensitive information.
    # 2. PYTHON. We set PYTHON to the empty string. This forces any packages that use
    #    PyCall to install their own version of Python instead of using the system
    #    Python.
    # 3. JULIA_DEPOT_PATH. We set JULIA_DEPOT_PATH to the temporary directory that
    #    we created. This is because we don't want the child process using our
    #    real Julia depot. So we set up a fake depot for the child process to use.
    cmd = Cmd(`$(Base.julia_cmd()) -e $(code)`;
              env = Dict("PATH" => ENV["PATH"],
                         "PYTHON" => "",
                         "JULIA_DEPOT_PATH" => tmp_dir))
    @info("Attempting to install the package")
    cmd_ran_successfully = success(pipeline(cmd, stdout=stdout, stderr=stderr))
    rm(tmp_dir; force = true, recursive = true)
    if cmd_ran_successfully
        @info("Successfully installed the package")
        return true, ""
    else
        @error("Was not able to successfully install the package")
        return false, "I was not able to install and import the package. See the Travis logs for details. Please make sure to read all of the Travis logs for all jobs. In particular, make sure to look at the logs for all Julia versions."
    end
end

url_has_correct_ending(url, pkg) = endswith(url, "/$(pkg).jl.git")


text:
@testset "cuRAND" begin

using CuArrays.CURAND

CURAND.seed!()

# in-place
for (f,T) in ((rand!,Float32),
              (randn!,Float32),
              (rand_logn!,Float32),
              (rand_poisson!,Cuint)),
    d in (2, (2,2), (2,2,2))
    A = CuArray{T}(undef, d)
    f(A)
end

# out-of-place, with implicit type
for (f,T) in ((curand,Float32), (curandn,Float32), (curand_logn,Float32),
              (curand_poisson,Cuint)),
    args in ((2,), (2, 2))
    A = f(args...)
    @test eltype(A) == T
end

# out-of-place, with type specified
for (f,T) in ((curand,Float32), (curandn,Float32), (curand_logn,Float32),
              (curand,Float64), (curandn,Float64), (curand_logn,Float64),
              (curand_poisson,Cuint)),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)))
    A = f(args...)
    @test eltype(A) == T
end

# unsupported types that fall back to GPUArrays
for (f,T) in ((curand,Int64),),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)))
    A = f(args...)
    @test eltype(A) == T
end
for (f,T) in ((rand!,Int64),),
    d in (2, (2,2), (2,2,2))
    A = CuArray{T}(undef, d)
    f(A)
end

end


text: (println(1); 1) + (println(2); 2)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function testNHHillY()
close all
nx=2;
ny=60;
lx=2000;
ly=60000;
x0=0;
y0=-30000;
Param.Grav=9.81;
Param.Coriolis=false;
Param.HyperD=0; %1.e15*(nx/30)^3.2;
Param.StrideDamp=6000;
Param.Relax=5.e-1;
Param.Damping=true;
Param.Flat=false;
Param.hC=400;
Param.y0C=0;
Param.aC=1000;
Param.H=15600;
Param.cS=360;
Param.Grav=9.81d0;
Param.Cpd=1004.0d0;
Param.Cvd=717.0d0;
Param.Rd=Param.Cpd-Param.Cvd;
Param.p0=1.0d5;
Param.Cpv=1885.0d0;
Param.Gamma=Param.Cpd/Param.Cvd;
Param.kappa=Param.Rd/Param.Cpd;
Param.Th0=300;
Param.uMax=0;
Param.vMax=10;
Param.NBr=1.e-2;
Param.TopoS='AgnesiCartY';
fig=1;
level=1;
OrdPoly=3;
CG.OrdPoly=OrdPoly;
Boundary.WE='Period';
Boundary.BT='Period';
Param.hS='';
Param.Grid=CartGrid(nx,ny,lx,ly,x0,y0,OrientFaceCart,Boundary,Param);

nz=40;
zP=zeros(nz,1);
z=zeros(nz+1,1);
dz=Param.H/nz;
zP(1)=dz/2;
for i=2:nz
  zP(i)=zP(i-1)+dz;
end
for i=2:nz+1
  z(i)=z(i-1)+dz;
end

Param.Grid.nz=nz;
Param.Grid.z=z;
Param.Grid.zP=zP;
Param.Grid.dz=dz;
[CG.Faces,CG.NumG,CG.NumI,CG.Glob,CG.FaceGlob,CG.Stencil]...
  =NumberingFemCG(Param.Grid,OrdPoly);
Param.JF=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces,nz+1);
Param.XF=zeros(OrdPoly+1,OrdPoly+1,3,Param.Grid.NumFaces,nz+1);
Param.dXdxF=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces,nz+1,3,3);
Param.dXdxIF=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces,nz+1,3,3);
[CG.w,CG.xw]=GaussLobattoQuad(CG.OrdPoly);
[CG.wX,CG.xwX]=GaussLobattoQuad(CG.OrdPoly);
[CG.wY,CG.xwY]=GaussLobattoQuad(CG.OrdPoly);
[CG.DW,CG.DS]=DerivativeMatrixSingle(CG);

for iF=1:Param.Grid.NumFaces
  for iz=1:nz+1
    [X,J,dXdx,dXdxI]=JacobiDG3(CG,Param.Grid.Faces(iF),z(iz),@Topo,Param);
    Param.XF(:,:,:,iF,iz)=X;
    Param.JF(:,:,iF,iz)=J;
    Param.dXdxF(:,:,iF,iz,:,:)=reshape(dXdx,OrdPoly+1,OrdPoly+1,1,3,3);
    Param.dXdxIF(:,:,iF,iz,:,:)=reshape(dXdxI,OrdPoly+1,OrdPoly+1,1,3,3);
  end
end
Param.X=0.5*(Param.XF(:,:,:,:,1:nz)+Param.XF(:,:,:,:,2:nz+1));
Param.J=0.5*(Param.JF(:,:,:,1:nz)+Param.JF(:,:,:,2:nz+1));
Param.dXdx(:,:,:,1:nz,:,:)=0.5*(Param.dXdxF(:,:,:,1:nz,:,:)...
                               +Param.dXdxF(:,:,:,2:nz+1,:,:));
Param.dXdxI(:,:,:,1:nz,:,:)=0.5*(Param.dXdxIF(:,:,:,1:nz,:,:)...
                                +Param.dXdxIF(:,:,:,2:nz+1,:,:));
[CG.M,CG.MW]=MassCG(CG,Param);

U=zeros(CG.NumG,nz,2);
W=zeros(CG.NumG,nz+1,1);
ProfRho='GravityHill';
ProfTheta='GravityHill';
ProfVel='Const';
Param.ProfVel=ProfVel;
Param.ProfRho=ProfRho;
Param.ProfTheta=ProfTheta;
Rho=Project(@fRho,CG,Param);
Slice.Type='XZ';
Slice.y=y0+ly/ny/2;
Slice.iy=1;
Param.SliceXZ=Slice;
Param.SliceYZ.Type='YZ';
Param.SliceYZ.ix=1;
Param.SliceYZ.x=x0+lx/nx/2;
Param.SliceXY.Type='XY';
Param.SliceXY.iz=1;
%PlotCG(Rho,CG,@TransCart,@Topo,Param,fig,Slice);
[U(:,:,1),U(:,:,2)]=ProjectVec(@fVel,CG,Param);
Th=Project(@fTheta,CG,Param).*Rho;
CFL=0.125;
dtau=.4;
time=0;

IntMethod='RungeKutta';
nIter=20000;
%v = VideoWriter ('Galewsky.avi');
%open (v);
switch IntMethod
  case 'RungeKutta'
    for i=1:nIter
      [FRho,FU,FW,FTh]=FcnNHCurlVec(Rho,U,W,Th,CG,Param);
      RhoNeu=Rho+1/3*dtau*FRho;
      UNeu=U+1/3*dtau*FU;
      WNeu=W+1/3*dtau*FW;
      ThNeu=Th+1/3*dtau*FTh;
      [FRho,FU,FW,FTh]=FcnNHCurlVec(RhoNeu,UNeu,WNeu,ThNeu,CG,Param);
      RhoNeu=Rho+1/2*dtau*FRho;
      UNeu=U+1/2*dtau*FU;
      WNeu=W+1/2*dtau*FW;
      ThNeu=Th+1/2*dtau*FTh;
      [FRho,FU,FW,FTh]=FcnNHCurlVec(RhoNeu,UNeu,WNeu,ThNeu,CG,Param);
      Rho=Rho+dtau*FRho;
      U=U+dtau*FU;
      W=W+dtau*FW;
      Th=Th+dtau*FTh;
      time=time+dtau;
      if mod(i,1000)==0

        fig=PlotCG(U(:,:,1),CG,@TransCart,@Topo,Param,fig,Param.SliceYZ);
        fig=PlotCG(U(:,:,2),CG,@TransCart,@Topo,Param,fig,Param.SliceYZ);
        fig=PlotCG(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,@TransCart...
          ,@Topo,Param,fig,Param.SliceYZ);
%         frame=getframe(gcf);
%         writeVideo (v, frame);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)
end


text:
# This file was generated, do not modify it. # hide
ŷ = inverse_transform(box_mach, ẑ)


text:
using Test

include("spiral-matrix.jl")

@testset "Different valid values" begin
    @testset "Empty spiral" begin
        @test spiral_matrix(0) == Matrix{Int}(undef,0,0)
    end
    @testset "Trivial spiral" begin
        @test spiral_matrix(1) == reshape([1],(1,1))
    end
    @testset "Spiral of size 2" begin
        @test spiral_matrix(2) == [1 2; 4 3]
    end
    @testset "Spiral of size 3" begin
        @test spiral_matrix(3) == [1 2 3; 8 9 4; 7 6 5]
    end
    @testset "Spiral of size 4" begin
        @test spiral_matrix(4) == [1 2 3 4; 12 13 14 5; 11 16 15 6; 10 9 8 7]
    end
    @testset "Spiral of size 5" begin
        @test spiral_matrix(5) == [1 2 3 4 5; 16 17 18 19 6; 15 24 25 20 7; 14 23 22 21 8; 13 12 11 10 9]
    end
end