text:
module TestDataFrame
    using Base.Test
    using DataFrames, Compat
    import Compat.String

    #
    # Equality
    #

    @test isequal(DataFrame(a=[1, 2, 3], b=[4, 5, 6]), DataFrame(a=[1, 2, 3], b=[4, 5, 6]))
    @test !isequal(DataFrame(a=[1, 2], b=[4, 5]), DataFrame(a=[1, 2, 3], b=[4, 5, 6]))
    @test !isequal(DataFrame(a=[1, 2, 3], b=[4, 5, 6]), DataFrame(a=[1, 2, 3]))
    @test !isequal(DataFrame(a=[1, 2, 3], b=[4, 5, 6]), DataFrame(a=[1, 2, 3], c=[4, 5, 6]))
    @test !isequal(DataFrame(a=[1, 2, 3], b=[4, 5, 6]), DataFrame(b=[4, 5, 6], a=[1, 2, 3]))
    @test !isequal(DataFrame(a=[1, 2, 2], b=[4, 5, 6]), DataFrame(a=[1, 2, 3], b=[4, 5, 6]))
    @test isequal(DataFrame(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]),
                  DataFrame(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]))

    # FIXME: equality operators won't work until JuliaStats/NullableArrays#84 is merged
    #@test get(DataFrame(a=[1, 2, 3], b=[4, 5, 6]) == DataFrame(a=[1, 2, 3], b=[4, 5, 6]))
    #@test get(DataFrame(a=[1, 2], b=[4, 5]) != DataFrame(a=[1, 2, 3], b=[4, 5, 6]))
    #@test get(DataFrame(a=[1, 2, 3], b=[4, 5, 6]) != DataFrame(a=[1, 2, 3]))
    #@test get(DataFrame(a=[1, 2, 3], b=[4, 5, 6]) != DataFrame(a=[1, 2, 3], c=[4, 5, 6]))
    #@test get(DataFrame(a=[1, 2, 3], b=[4, 5, 6]) != DataFrame(b=[4, 5, 6], a=[1, 2, 3]))
    #@test get(DataFrame(a=[1, 2, 2], b=[4, 5, 6]) != DataFrame(a=[1, 2, 3], b=[4, 5, 6]))
    #@test get(DataFrame(a=Nullable{Int}[1, 3, Nullable()], b=[4, 5, 6]) !=
    #          DataFrame(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]))
    #@test isnull(DataFrame(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]) ==
    #             DataFrame(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]))
    #@test isnull(DataFrame(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]) ==
    #             DataFrame(a=Nullable{Int}[1, 2, 3], b=[4, 5, 6]))

    #
    # Copying
    #

    df = DataFrame(a = [2, 3], b = Any[DataFrame(c = 1), DataFrame(d = 2)])
    dfc = copy(df)
    dfdc = deepcopy(df)

    df[1, :a] = 4
    get(df[1, :b])[:e] = 5
    names!(df, [:f, :g])

    @test names(dfc) == [:a, :b]
    @test names(dfdc) == [:a, :b]

    @test get(dfc[1, :a]) === 4
    @test get(dfdc[1, :a]) === 2

    @test names(get(dfc[1, :b])) == [:c, :e]
    @test names(get(dfdc[1, :b])) == [:c]

    #

    x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])
    v = DataFrame(a = [5, 6, 7], b = [8, 9, 10])

    z = vcat(v, x)

    z2 = z[:, [1, 1, 2]]
    @test names(z2) == [:a, :a_1, :b]

    #test_group("DataFrame assignment")
    # Insert single column
    x0 = x[Int[], :]
    @test_throws ErrorException x0[:d] = [1]
    @test_throws ErrorException x0[:d] = 1:3

    # Insert single value
    x[:d] = 3
    @test isequal(x[:d], NullableArray([3, 3, 3]))

    x0[:d] = 3
    @test x0[:d] == Int[]

    # similar / nulls
    df = DataFrame(a = 1, b = "b", c = CategoricalArray([3.3]))
    nulldf = DataFrame(a = NullableArray(Int, 2),
                       b = NullableArray(String, 2),
                       c = NullableCategoricalArray(Float64, 2))
    @test isequal(nulldf, similar(df, 2))
    @test isequal(nulldf, DataFrames.similar_nullable(df, 2))

    # Associative methods

    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test haskey(df, :a)
    @test !haskey(df, :c)
    @test get(df, :a, -1) === df.columns[1]
    @test get(df, :c, -1) == -1
    @test !isempty(df)

    @test empty!(df) === df
    @test isempty(df.columns)
    @test isempty(df)

    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test_throws BoundsError insert!(df, 5, ["a", "b"], :newcol)
    @test_throws ErrorException insert!(df, 1, ["a"], :newcol)
    @test isequal(insert!(df, 1, ["a", "b"], :newcol), df)
    @test names(df) == [:newcol, :a, :b]
    @test isequal(df[:a], NullableArray([1, 2]))
    @test isequal(df[:b], NullableArray([3., 4.]))
    @test isequal(df[:newcol], ["a", "b"])

    df = DataFrame(a=[1, 2], b=[3., 4.])
    df2 = DataFrame(b=["a", "b"], c=[:c, :d])
    @test isequal(merge!(df, df2), df)
    @test isequal(df, DataFrame(a=[1, 2], b=["a", "b"], c=[:c, :d]))

    #test_group("Empty DataFrame constructors")
    df = DataFrame(Int, 10, 3)
    @test size(df, 1) == 10
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == NullableVector{Int}
    @test typeof(df[:, 2]) == NullableVector{Int}
    @test typeof(df[:, 3]) == NullableVector{Int}
    @test allnull(df[:, 1])
    @test allnull(df[:, 2])
    @test allnull(df[:, 3])

    df = DataFrame(Any[Int, Float64, String], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == NullableVector{Int}
    @test typeof(df[:, 2]) == NullableVector{Float64}
    @test typeof(df[:, 3]) == NullableVector{String}
    @test allnull(df[:, 1])
    @test allnull(df[:, 2])
    @test allnull(df[:, 3])

    df = DataFrame(Any[Int, Float64, String], [:A, :B, :C], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == NullableVector{Int}
    @test typeof(df[:, 2]) == NullableVector{Float64}
    @test typeof(df[:, 3]) == NullableVector{String}
    @test allnull(df[:, 1])
    @test allnull(df[:, 2])
    @test allnull(df[:, 3])

    df = DataFrame(DataType[Int, Float64, Compat.UTF8String],[:A, :B, :C], [false,false,true],100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == NullableVector{Int}
    @test typeof(df[:, 2]) == NullableVector{Float64}
    @test typeof(df[:, 3]) == NullableCategoricalVector{Compat.UTF8String,UInt32}
    @test allnull(df[:, 1])
    @test allnull(df[:, 2])
    @test allnull(df[:, 3])

    df = convert(DataFrame, zeros(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, ones(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, eye(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    #test_group("Other DataFrame constructors")
    df = DataFrame([@compat(Dict{Any,Any}(:a=>1, :b=>'c')),
                    @compat(Dict{Any,Any}(:a=>3, :b=>'d')),
                    @compat(Dict{Any,Any}(:a=>5))])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == NullableVector{Int}
    @test typeof(df[:,:b]) == NullableVector{Char}

    df = DataFrame([@compat(Dict{Any,Any}(:a=>1, :b=>'c')),
                    @compat(Dict{Any,Any}(:a=>3, :b=>'d')),
                    @compat(Dict{Any,Any}(:a=>5))],
                   [:a, :b])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == NullableVector{Int}
    @test typeof(df[:,:b]) == NullableVector{Char}

    # This assignment was missing before
    df = DataFrame(Column = [:A])
    df[1, :Column] = "Testing"

    # zero-row dataframe and subdataframe test
    df = DataFrame(x=[], y=[])
    @test nrow(df) == 0
    df = DataFrame(x=[1:3;], y=[3:5;])
    sdf = sub(df, df[:x] .== 4)
    @test size(sdf, 1) == 0

    @test hash(convert(DataFrame, [1 2; 3 4])) == hash(convert(DataFrame, [1 2; 3 4]))
    @test hash(convert(DataFrame, [1 2; 3 4])) != hash(convert(DataFrame, [1 3; 2 4]))

    # push!(df, row)
    df=DataFrame( first=[1,2,3], second=["apple","orange","pear"] )

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, Any[3,"pear"])
    @test isequal(df, dfb)

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, (3,"pear"))
    @test isequal(df, dfb)

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, (33.33,"pear"))

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, ("coconut",22))

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, @compat(Dict(:first=>3, :second=>"pear")))
    @test isequal(df, dfb)

    df=DataFrame( first=[1,2,3], second=["apple","orange","banana"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, @compat(Dict("first"=>3, "second"=>"banana")))
    @test isequal(df, dfb)

    df0= DataFrame( first=[1,2], second=["apple","orange"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, @compat(Dict(:first=>true, :second=>false)))
    @test isequal(df0, dfb)

    df0= DataFrame( first=[1,2], second=["apple","orange"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, @compat(Dict("first"=>"chicken", "second"=>"stuff")))
    @test isequal(df0, dfb)

    # delete!
    df = DataFrame(a=1, b=2, c=3, d=4, e=5)
    @test_throws ArgumentError delete!(df, 0)
    @test_throws ArgumentError delete!(df, 6)
    @test_throws KeyError delete!(df, :f)

    d = copy(df)
    delete!(d, [:a, :e, :c])
    @test names(d) == [:b, :d]
    delete!(d, :b)
    @test isequal(d, df[[:d]])

    d = copy(df)
    delete!(d, [2, 5, 3])
    @test names(d) == [:a, :d]
    delete!(d, 2)
    @test isequal(d, df[[:a]])

    # deleterows!
    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test deleterows!(df, 1) === df
    @test isequal(df, DataFrame(a=[2], b=[4.]))

    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test deleterows!(df, 2) === df
    @test isequal(df, DataFrame(a=[1], b=[3.]))

    df = DataFrame(a=[1, 2, 3], b=[3., 4., 5.])
    @test deleterows!(df, 2:3) === df
    @test isequal(df, DataFrame(a=[1], b=[3.]))

    df = DataFrame(a=[1, 2, 3], b=[3., 4., 5.])
    @test deleterows!(df, [2, 3]) === df
    @test isequal(df, DataFrame(a=[1], b=[3.]))

    df = DataFrame(a=NullableArray([1, 2]), b=NullableArray([3., 4.]))
    @test deleterows!(df, 1) === df
    @test isequal(df, DataFrame(a=NullableArray([2]), b=NullableArray([4.])))

    df = DataFrame(a=NullableArray([1, 2]), b=NullableArray([3., 4.]))
    @test deleterows!(df, 2) === df
    @test isequal(df, DataFrame(a=NullableArray([1]), b=NullableArray([3.])))

    df = DataFrame(a=NullableArray([1, 2, 3]), b=NullableArray([3., 4., 5.]))
    @test deleterows!(df, 2:3) === df
    @test isequal(df, DataFrame(a=NullableArray([1]), b=NullableArray([3.])))

    df = DataFrame(a=NullableArray([1, 2, 3]), b=NullableArray([3., 4., 5.]))
    @test deleterows!(df, [2, 3]) === df
    @test isequal(df, DataFrame(a=NullableArray([1]), b=NullableArray([3.])))

    # describe
    #suppress output and test that describe() does not throw
    devnull = is_unix() ? "/dev/null" : "nul"
    open(devnull, "w") do f
        @test nothing == describe(f, DataFrame(a=[1, 2], b=Any["3", Nullable()]))
        @test nothing ==
              describe(f, DataFrame(a=NullableArray([1, 2]),
                                    b=NullableArray(Nullable{String}["3", Nullable()])))
        @test nothing ==
              describe(f, DataFrame(a=CategoricalArray([1, 2]),
                                    b=NullableCategoricalArray(Nullable{String}["3", Nullable()])))
        @test nothing == describe(f, [1, 2, 3])
        @test nothing == describe(f, NullableArray([1, 2, 3]))
        @test nothing == describe(f, CategoricalArray([1, 2, 3]))
        @test nothing == describe(f, Any["1", "2", Nullable()])
        @test nothing == describe(f, NullableArray(Nullable{String}["1", "2", Nullable()]))
        @test nothing == describe(f, NullableCategoricalArray(Nullable{String}["1", "2", Nullable()]))
    end

    #Check the output of unstack
    df = DataFrame(Fish = CategoricalArray(["Bob", "Bob", "Batman", "Batman"]),
                   Key = ["Mass", "Color", "Mass", "Color"],
                   Value = ["12 g", "Red", "18 g", "Grey"])
    # Check that reordering levels does not confuse unstack
    levels!(df[1], ["XXX", "Bob", "Batman"])
    #Unstack specifying a row column
    df2 = unstack(df,:Fish, :Key, :Value)
    #Unstack without specifying a row column
    df3 = unstack(df,:Key, :Value)
    #The expected output
    df4 = DataFrame(Fish = ["XXX", "Bob", "Batman"],
                    Color = Nullable{String}[Nullable(), "Red", "Grey"],
                    Mass = Nullable{String}[Nullable(), "12 g", "18 g"])
    @test isequal(df2, df4)
    @test isequal(df3, df4[2:3, :])
    #Make sure unstack works with NULLs at the start of the value column
    df[1,:Value] = Nullable()
    df2 = unstack(df,:Fish, :Key, :Value)
    #This changes the expected result
    df4[2,:Mass] = Nullable()
    @test isequal(df2, df4)
end


text:
module RecurrentCorticalNetwork
using LightGraphs, SimpleWeightedGraphs, Distributions, Plots
using SpikingNeuralNets, SpikingNeuralNets.LIFs, SpikingNeuralNets.GraphUtils, SpikingNeuralNets.SNNPlot, SpikingNeuralNets.Inputs

# Set up the network
const N = 100
const fE = 0.8
const fI = 1.0 - fE
const nE = round(Int, fE*N)
const nI = round(Int, fI*N)
const f = 0.15  # percentage of excitatory neurons selective to input
const nA = round(Int, f*nE)
const nB = nA
const nNS = nE - nA - nB
const excitatoryP = 1.7   # the weight of strong excitatory connections
const excitatoryD = 1 - f*(excitatoryP-1)/(1-f)    # the weight of weak excitatory connections
const inhibitory = -1.0
const neutral = 1.0

# Vertices:
#   1- A neurons
#   2- B neurons
#   3- Non-selective neurons
#   4- Interneurons
const metagraph = SimpleWeightedDiGraph(4)
add_edge!(metagraph, 1, 1, excitatoryP)  # A connections
add_edge!(metagraph, 1, 2, excitatoryD)
add_edge!(metagraph, 1, 3, neutral)
add_edge!(metagraph, 2, 1, excitatoryD)  # B connections
add_edge!(metagraph, 2, 2, excitatoryP)
add_edge!(metagraph, 2, 3, neutral)
add_edge!(metagraph, 3, 1, excitatoryD)  # non-selective connections
add_edge!(metagraph, 3, 2, excitatoryD)
add_edge!(metagraph, 3, 3, neutral)
add_edge!(metagraph, 1, 4, neutral)      # interneurons have excitatory inputs
add_edge!(metagraph, 2, 4, neutral)
add_edge!(metagraph, 3, 4, neutral)
add_edge!(metagraph, 4, 1, inhibitory)   # and inhibitory outputs
add_edge!(metagraph, 4, 2, inhibitory)
add_edge!(metagraph, 4, 3, inhibitory)
add_edge!(metagraph, 4, 4, inhibitory)
const graph = ComplexGraph(metagraph, [nA, nB, nNS, nI], selfloops=false) # neurons do not excite themselves
add_inputs!(graph, vertices(graph), fill(neutral, nv(graph)))        # add external inputs for all neurons

# Predicates for neuronal populations
const aneurons = map(x -> x <= nA, vertices(graph))
const bneurons = map(x -> x > nA && x <= nA+nB, vertices(graph))
const nsneurons = map(x -> x > nA+nB && x <= nE, vertices(graph))
const pyramidalneurons = @. aneurons | bneurons | nsneurons
const interneurons = map(x -> x > nE && x <= nE+nI, vertices(graph))
const inputneurons = map(x -> x > N, vertices(graph))
const ainputneurons = map(x -> x > N && x <= N+nA+nB, vertices(graph))
const binputneurons = map(x -> x > N && x <= N+nA+nB, vertices(graph))
const nsinputneurons = map(x -> x > N && x <= N+nE, vertices(graph))
const iinputneurons = map(x -> x > N+nE, vertices(graph))

# Define input parameters
const μ0 = 40 # Hz
const ρa = μ0 / 100
const ρb = μ0 / 100
const σ = 4

# Define Leaky Integrate-and-Fire Parameters
const dt = 0.02e-3 # Timestep (0.02ms)
const Vrest = -0.07 # Resting potential (-70 mV)
const Vθ = -0.05 # Firing threshold (-50 mV)
const Vreset = -0.055 # Reset potential (-55 mv)
const Ve = 0.0      # excitatory reversal potential (0 mV)
const Vi = -70e-3   # inhibitory reversal potential (-70 mV)
const Mg = 1e-3    # Mg concentration (1 mM)
const Vext = 2.4e3  # Rate of spontaneous external inputs (2.4 kHz)
const Wext = [1.0]   # weight of spontaneous external inputs

Cm(pyramidal::Bool)::Real = pyramidal ? 0.5e-9 : 0.2e-9  # membrane capacitance (nF)
Rm(pyramidal::Bool)::Real = 1.0 / Cm(pyramidal)
gL(pyramidal::Bool)::Real = pyramidal ? 25e-9 : 20e-9    # membrane leak conductance (nS)
τref(pyramidal::Bool)::Real = pyramidal ? 0.002 : 0.001        # refractory period (mS)
τ(pyramidal::Bool)::Real = Cm(pyramidal) / gL(pyramidal)

# synaptic conductances (nS)
gextAMPA(pyramidal::Bool)::Real = pyramidal ? 1.5e-9 : 1.25e-9 # 2.1e-9 : 1.62e-9
grecAMPA(pyramidal::Bool)::Real = pyramidal ? 0.05e-9 : 0.04e-9 # 0.05e-9 : 0.04e-9
gNMDA(pyramidal::Bool)::Real = pyramidal ? 0.165e-9 : 0.13e-9 # 0.165e-9 : 0.13e-9
gGABA(pyramidal::Bool)::Real = pyramidal ? 4.5e-9 : 4.0e-9 # 1.3e-9 : 1.0e-9
# synaptic decay times
const τAMPA = 2e-3/dt
const α = 500.0*dt
const τNMDAd = 100e-3/dt # decay
const τNMDAr = 2e-3/dt # rise
const τGABA = 5e-3/dt

# Synaptic latency
const synl = 0.5e-3  # delay in seconds
const synt = 0 # round(Int, synl / dt)  # delay in timesteps

# synaptic current and gating functions
function IextAMPA(lif::LIF, n::Integer, s::AbstractVector{<:Real})::Real
    return -gextAMPA(pyramidalneurons[n]) * (voltage(lif, n) - Ve) * s[1]
end
return function sextAMPA(lif::LIF{VT}, n::Integer, s::AbstractVector{VT})::Vector{VT} where {VT <: Real}
    input = spikes(lif, inputs(lif, n); t=synt)
    return @. -(s/τAMPA) + input
end

function IrecAMPA(lif::LIF{VT}, n::Integer, s::AbstractVector{VT}; inputs=false)::Real where {VT <: Real}
    return -grecAMPA(pyramidalneurons[n]) * (voltage(lif, n) - Ve) *
            sum(weights(lif.graph)[excitors(lif, n; inputs=inputs), n] .* s)
end
function srecAMPA(lif::LIF{VT}, n::Integer, s::AbstractVector{VT}; inputs=false)::Vector{VT} where {VT <: Real}
    input = spikes(lif, excitors(lif, n; inputs=inputs); t=synt)
    return @. -(s/τAMPA) + input
end

function INMDA(lif::LIF{VT}, n::Integer, vars::Matrix{VT}; inputs=false)::Real where {VT <: Real}
    v = voltage(lif, n)
    return -gNMDA(pyramidalneurons[n]) * (v - Ve) * sum(weights(lif.graph)[excitors(lif, n; inputs=inputs), n] .* view(vars, :, 1)) /
            (1 + Mg * exp(-0.062*v) / 3.57)
end
function xNMDA(lif::LIF{VT}, n::Integer, x::AbstractVector{VT}; inputs=false)::Vector{VT} where {VT <: Real}
    input = spikes(lif, excitors(lif, n; inputs=inputs); t=synt)
    return @. -(x/τNMDAr) + input
end
function sNMDA(lif::LIF{VT}, s::AbstractVector{VT}, x::AbstractVector{VT})::Vector{VT} where {VT <: Real}
    return @. -(s/τNMDAd) + α*x*(1 - s)
end
@views function varsNMDA(lif::LIF{VT}, n::Integer, vars::Matrix{VT}; inputs=false)::Matrix{VT} where {VT <: Real}
    dvars = similar(vars)
    dvars[:, 1] = sNMDA(lif, vars[:, 1], vars[:, 2])
    dvars[:, 2] = xNMDA(lif, n, vars[:, 2]; inputs=inputs)
    return dvars
end

function IGABA(lif::LIF{VT}, n::Integer, s::AbstractVector{VT}; inputs=false)::Real where {VT <: Real}
    return -gGABA(pyramidalneurons[n]) * (voltage(lif, n) - Vi) * sum(s)
end
function sGABA(lif::LIF{VT}, n::Integer, s::AbstractVector{VT}; inputs=false)::Vector{VT} where {VT <: Real}
    input = spikes(lif, inhibitors(lif, n; inputs=inputs); t=synt)
    return @. -(s/τGABA) + input
end

function singleneurontest(;pyramidal=true, x=x=0:1e-13:10e-11, iterations=25000)
    rates = Vector{Float64}(undef, length(x))
    # Create a network of a single neuron
    g = SimpleWeightedDiGraph(1)
    add_edge!(g, 1, 1, 1.0)
    @Threads.threads for i in eachindex(x)
        lif = LIF{Float64, Int}(g; m=1, dt=dt, Vrest=[Vrest],
                                Vθ=[Vθ], Vreset=[Vreset],
                                τ=[τ(pyramidal)], R=[Rm(pyramidal)],
                                τref=[τref(pyramidal)])
        I(lif, n, s) = x[i]
        s(lif, n, s) = zeros(size(s))
        add_channel!(lif, "input", I, s, [zeros(1)])
        step!(lif, iterations)
        rates[i] = mean(spiketrain(lif, iterations)) / dt
    end
    return plot(x, rates, xlab="Input Current (V)", ylab="Firing Rate (Hz)")
end

function synapsetest(;pyramidal=true, channel="AMPA", x=x=0:500, iterations=10000)
    gating = Matrix{Float64}(undef, length(x), channel=="NMDA" ? 2 : 1)
    rates = Vector{Float64}(undef, length(x))

    # set the current and gating functions, including input neurons
    I = if channel=="AMPA"
         (args...)->IrecAMPA(args...; inputs=true)
     elseif channel=="NMDA"
         (args...)->INMDA(args...; inputs=true)
     else
         (args...)->IGABA(args...; inputs=true)
     end
    s = if channel=="AMPA"
        (args...)->srecAMPA(args...; inputs=true)
    elseif channel=="NMDA"
        (args...)->varsNMDA(args...; inputs=true)
    else
        (args...)->sGABA(args...; inputs=true)
    end
    gate(lif) = vars(lif, channel, 2)

    # Create a network of three neurons (input -> presynaptic neuron -> postsynaptic neuron)
    pyramidal = [channel!="GABA", pyramidal]
    g = SimpleWeightedDiGraph(2)
    add_edge!(g, 1, 2, (pyramidal[1] ? 1.0 : -1.0))
    @Threads.threads for i in eachindex(x)
        input = rand(Poisson(x[i]*dt), 1, iterations)
        lif = LIF{Float64, Int}(g; m=1, dt=dt, Vrest=fill(Vrest, 2),
                                Vθ=fill(Vθ, 2), Vreset=fill(Vreset, 2),
                                τ=τ.(pyramidal), R=Rm.(pyramidal),
                                τref=τref.(pyramidal))
        add_channel!(lif, channel, I, s, fill(channel=="NMDA" ? zeros(1,2) : zeros(1), 2), [2])
        lif, V, G = step!(lif, input, voltages, gate)
        V = reduce(hcat, V)[2,:]
        G = reduce(vcat, G)
        gating[i,:] = mean(G, dims=1)
        rates[i] = mean(spiketrain(lif, iterations)) / dt

        if isone(length(x))
            return plot(plot((1:iterations)*dt, G, xlab="Time (s)", ylab="Gates", ylim=[0.0, 1.0]),
                        plot((1:iterations)*dt, V, xlab="Time (s)", ylab="Voltage (V)"),
                        rasterplot(spiketrain(lif, iterations), timestep=dt), layout=(3,1))
        end
    end
    return plot(plot(x, gating, xlab="Input Rate (Hz)", ylab="Average Gating Variable"),
                plot(x, rates, xlab="Input Rate (Hz)", ylab="Output Firing Rate (Hz)"),
                layout=(2,1))
end

function run(;c::Real=0.0, iterations::Integer=10, window::Integer=100, extStart=1, extStop=iterations)
    GC.gc()
    μa = μ0 + ρa*c
    μb = μ0 - ρb*c
    ext = ConstantPoissonInput(Vext*dt, N, iterations)
    stim = VariablePoissonInput([fill(μa*dt, nA); fill(μb*dt, nB); zeros(nNS+nI)], [fill(σ*dt, nA+nB); zeros(nNS+nI)], iterations; resample=floor(Int, 50e-3/dt), start=extStart, stop=extStop)
    I = SummedInput(ext, stim)

    # create the LIF and add all of the channel types
    lif = LIF{Float64, Int}(graph; m=1, dt=dt, Vrest=fill(Vrest, nv(graph)),
                            Vθ=fill(Vθ, nv(graph)),
                            Vreset=fill(Vreset, nv(graph)),
                            τ=τ.(pyramidalneurons),
                            R=Rm.(pyramidalneurons),
                            τref=τref.(pyramidalneurons))

    # Add excitatory and inhibitory channels
    add_channel!(lif, "extAMPA", IextAMPA, sextAMPA, [zeros(1) for n in neurons(lif)])
    add_channel!(lif, "recAMPA", IrecAMPA, srecAMPA, [zeros(length(excitors(lif, n))) for n in neurons(lif)])
    add_channel!(lif, "NMDA", INMDA, varsNMDA, [zeros(length(excitors(lif, n)), 2) for n in neurons(lif)])
    add_channel!(lif, "GABA", IGABA, sGABA, [zeros(length(inhibitors(lif, n))) for n in neurons(lif)])

    lif, V = step!(lif, I, voltages; threaded=true)
    V = reduce(hcat, V)
    S = spiketrain(lif, iterations)

    #println(lif.S)
    #println(S)

    display(plot(rasterplot(S[aneurons,:]; title="A Neurons", timestep=dt),
                 rateplot(S[aneurons,:]; timestep=dt, window=window),
                 voltageplot(V[aneurons,:]; timestep=dt, ylim=[Vrest, Vθ+0.0025]);
                 xlim=[0.0, iterations*dt], layout=(3,1)))
    display(plot(rasterplot(S[bneurons,:]; title="B Neurons", timestep=dt),
                 rateplot(S[bneurons,:]; timestep=dt, window=window),
                 voltageplot(V[bneurons,:]; timestep=dt, ylim=[Vrest, Vθ+0.0025]);
                 xlim=[0.0, iterations*dt], layout=(3,1)))
    display(plot(rasterplot(S[nsneurons,:]; title="Non-selective Neurons", timestep=dt),
                 rateplot(S[nsneurons,:]; timestep=dt, window=window),
                 voltageplot(V[nsneurons,:]; timestep=dt, ylim=[Vrest, Vθ+0.0025]);
                 xlim=[0.0, iterations*dt], layout=(3,1)))
    display(plot(rasterplot(S[interneurons,:]; title="Interneurons", timestep=dt),
                 rateplot(S[interneurons,:]; timestep=dt, window=window),
                 voltageplot(V[interneurons,:]; timestep=dt, ylim=[Vrest, Vθ+0.0025]);
                 xlim=[0.0, iterations*dt], layout=(3,1)))
    #display(rArBplot(S[aneurons,:], S[bneurons,:], window=window, timestep=dt,
    #                 xlim=[0.0, 50], ylim=[0.0, 50], aspect_ratio=:equal,
    #                 xlab="Firing Rate r_A (Hz)", ylab="Firing Rate r_B (Hz)", leg=false))
    #display(rArBgif(S[aneurons,:], S[bneurons,:], window=window, timestep=dt, tail=100, fps=10, skip=50,
    #                  xlim=[0.0, 1/τref(true)], ylim=[0.0, 1/τref(true)], aspect_ratio=:equal,
    #                  xlab="Firing Rate r_A (Hz)", ylab="Firing Rate r_B (Hz)", leg=false))
    GC.gc()
end

end

ext: jl


text:
# PATH это массив, который хранит путь из левой нижней клетки в начальную.
PATH = []
# path(x,y) это массив, который хранит начальные коррдинаты робота.
path = [0, 0]

# movements_in_angle двигает робота из начальной клетки в левую нижнюю.
function moving_in_angle!(r)::Vector{<:Integer}
    while not(isborder(r,Sud) && isborder(r,West))
        for side in (West, Sud)
            steps = 0
            while not(isborder(r,side))
                move!(r,side)
                steps += 1
            push!(PATH, steps)
        end
    end
    return reverse!(PATH)
end

# pathXY складывает все соответсвтеные координаты пути, записанного в PATH и получает начальные координаты робота.
function pathXY()::Vector{<:Integer}
    for i in lenght(PATH)
        path[2 - (i%2)] += PATH[i]*(i%2==0) + PATH[i]*(i%2!=0)
    end
end

# moving_back_to_start двигает робота из нижней левой клетки в начальную по обратному пути из массива PATH
function moving_back_to_start!(r)::Nothing
    for i in lenght(PATH)
        for a in PATH[i]
            move!(r, HorizonSide((i+2)%4))
        end
    end
end

# tuda_suda делает туда-сюда
function tuda_suda(side)::Integer
    steps = 0
    while not(isborder(side))
        move!(r,side)
        steps += 1
    while not(ismarked(HorizonSide((Int(side)+2)%4))
        move!(r,HorizonSide((Int(side)+2)%4))
    return steps
end

# moving_around двигает робота по внутренним сторанам квадрата
function moving_around!(r)::Nothing
    steps_to_marker = [0,0]

    for steps in path[2]
        move!(r,Nord)
        steps_to_marker[1] += 1
    end
    putmarker!(r)
    while (not(isborder(Nord)))
        move!(r,Nord)
    end

    for steps in path[1]
        move!(r,Ost)
        steps_to_marker[2] += 1
    end
    putmarker!(r)
    while (not(isborder(Ost)))
        move!(r,Ost)
    end

    for steps in steps_to_marker[1]
        move!(r,Sud)
    end
    putmarker!(r)
    while (not(isborder(Sud)))
        move!(r,Sud)
    end

    for steps in steps_to_marker[2]
        move!(r,West)
    end
    putmarker!(r)
    while (not(isborder(West)))
        move!(r,West)
    end
end
# master! главная функция программы
function master(r)::Nothing
    moving_in_angle!(r)
    moving_around!(r)
    moving_back_to_start!(r)
end

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
module Info

using Parameters

@with_kw struct MapInfo
    n_pts::Int64 = 100
    n_transient::Int64 = 10000
    stepsize::Int64 = 1
    noise_frac::Float64 = 0.0
    seed::Array{Int32,1} = reinterpret(Int32, Base.Random.GLOBAL_RNG.seed)
end

export MapInfo
end


text:
@testset "Graph coloring" begin

function normal_49_states(;reverse_constraint = false)
    m = Model(with_optimizer(CS.Optimizer, keep_logs=true))
    num_colors = 20

    @variable(m, 1 <= states[1:50] <= num_colors, Int)
    @variable(m, 1 <= max_color <= num_colors, Int)

    washington,montana,maine,north_dakota,south_dakota,wyoming,wisconsin,idaho,vermont,minnesota,oregon,new_hampshire,
    iowa,massachusetts,nebraska,new_york,pennsylvania,connecticut,rhode_island,new_jersey,indiana,nevada,utah,california,ohio,
    illinois,washington_dc,delaware,west_virginia,maryland,colorado,kentucky,kansas,virginia,missouri,arizona,oklahoma,north_carolina,
    tennessee,texas,new_mexico,alabama,mississippi,georgia,south_carolina,arkansas,louisiana,florida,michigan = states

    @constraint(m, washington != oregon)
    @constraint(m, washington != idaho)
    @constraint(m, oregon != idaho)
    @constraint(m, oregon != nevada)
    @constraint(m, oregon != california)
    @constraint(m, california != nevada)
    @constraint(m, california != arizona)
    @constraint(m, nevada != idaho)
    @constraint(m, nevada != utah)
    @constraint(m, nevada != arizona)
    @constraint(m, idaho != montana)
    @constraint(m, idaho != wyoming)
    @constraint(m, idaho != utah)
    @constraint(m, utah != wyoming)
    @constraint(m, utah != colorado)
    @constraint(m, utah != new_mexico)
    @constraint(m, utah != arizona)
    @constraint(m, arizona != colorado)
    @constraint(m, arizona != new_mexico)
    @constraint(m, montana != north_dakota)
    @constraint(m, montana != south_dakota)
    @constraint(m, montana != wyoming)
    @constraint(m, wyoming != south_dakota)
    @constraint(m, wyoming != nebraska)
    @constraint(m, wyoming != colorado)
    @constraint(m, colorado != nebraska)
    @constraint(m, colorado != kansas)
    @constraint(m, colorado != oklahoma)
    @constraint(m, colorado != new_mexico)
    @constraint(m, new_mexico != oklahoma)
    @constraint(m, new_mexico != texas)
    @constraint(m, north_dakota != minnesota)
    @constraint(m, north_dakota != south_dakota)
    @constraint(m, south_dakota != minnesota)
    @constraint(m, south_dakota != iowa)
    @constraint(m, south_dakota != nebraska)
    @constraint(m, nebraska != iowa)
    @constraint(m, nebraska != missouri)
    @constraint(m, nebraska != kansas)
    @constraint(m, kansas != missouri)
    @constraint(m, kansas != oklahoma)
    @constraint(m, oklahoma != arkansas)
    @constraint(m, oklahoma != texas)
    @constraint(m, texas != arkansas)
    @constraint(m, texas != louisiana)
    @constraint(m, minnesota != wisconsin)
    @constraint(m, minnesota != iowa)
    @constraint(m, iowa != wisconsin)
    @constraint(m, iowa != illinois)
    @constraint(m, iowa != missouri)
    @constraint(m, missouri != illinois)
    @constraint(m, missouri != kentucky)
    @constraint(m, missouri != tennessee)
    @constraint(m, missouri != arkansas)
    @constraint(m, arkansas != tennessee)
    @constraint(m, arkansas != mississippi)
    @constraint(m, arkansas != louisiana)
    @constraint(m, louisiana != mississippi)
    @constraint(m, wisconsin != illinois)
    @constraint(m, wisconsin != michigan)
    @constraint(m, illinois != indiana)
    @constraint(m, illinois != kentucky)
    @constraint(m, kentucky != indiana)
    @constraint(m, kentucky != ohio)
    @constraint(m, kentucky != west_virginia)
    @constraint(m, kentucky != virginia)
    @constraint(m, kentucky != tennessee)
    @constraint(m, tennessee != virginia)
    @constraint(m, tennessee != north_carolina)
    @constraint(m, tennessee != georgia)
    @constraint(m, tennessee != alabama)
    @constraint(m, tennessee != mississippi)
    @constraint(m, mississippi != alabama)
    @constraint(m, michigan != indiana)
    @constraint(m, michigan != ohio)
    @constraint(m, indiana != ohio)
    @constraint(m, alabama != georgia)
    @constraint(m, alabama != florida)
    @constraint(m, ohio != pennsylvania)
    @constraint(m, ohio != west_virginia)
    @constraint(m, maine != new_hampshire)
    @constraint(m, new_hampshire != vermont)
    @constraint(m, new_hampshire != massachusetts)
    @constraint(m, vermont != massachusetts)
    @constraint(m, vermont != new_york)
    @constraint(m, massachusetts != new_york)
    @constraint(m, massachusetts != rhode_island)
    @constraint(m, massachusetts != connecticut)
    @constraint(m, connecticut != rhode_island)
    @constraint(m, new_york != connecticut)
    @constraint(m, new_york != pennsylvania)
    @constraint(m, new_york != new_jersey)
    @constraint(m, pennsylvania != new_jersey)
    @constraint(m, pennsylvania != delaware)
    @constraint(m, pennsylvania != maryland)
    @constraint(m, pennsylvania != west_virginia)
    @constraint(m, new_jersey != delaware)
    @constraint(m, maryland != washington_dc)
    @constraint(m, maryland != west_virginia)
    @constraint(m, maryland != virginia)
    @constraint(m, washington_dc != virginia)
    @constraint(m, west_virginia != virginia)
    @constraint(m, virginia != north_carolina)
    @constraint(m, north_carolina != south_carolina)
    @constraint(m, south_carolina != georgia)
    @constraint(m, georgia != florida)

    if reverse_constraint
        @constraint(m, max_color .>= states)
    else
        @constraint(m, states .<= max_color)
    end

    @objective(m, Min, max_color)

    optimize!(m)

    status = JuMP.termination_status(m)

    com = JuMP.backend(m).optimizer.model.inner

    CS.save_logs(com, "graph_color_optimize.json")
    rm("graph_color_optimize.json")

    @test status == MOI.OPTIMAL
    @test com.best_sol == 4 == JuMP.objective_value(m) == JuMP.objective_bound(m)
    @test maximum([JuMP.value(var) for var in states]) == 4 == JuMP.value(max_color)

    return com
end

@testset "49 US states + DC" begin
    com1 = normal_49_states(;reverse_constraint=true)
    com2 = normal_49_states()
    info_1 = com1.info
    info_2 = com2.info
    @test info_1.pre_backtrack_calls == info_2.pre_backtrack_calls
    @test info_1.backtrack_fixes == info_2.backtrack_fixes
    @test info_1.in_backtrack_calls == info_2.in_backtrack_calls
    @test info_1.backtrack_reverses == info_2.backtrack_reverses
    logs_1 = CS.get_logs(com1)
    logs_2 = CS.get_logs(com2)
    @test CS.same_logs(logs_1[:tree], logs_2[:tree])
end

@testset "49 US states + DC without sorting" begin
    m = Model(with_optimizer(CS.Optimizer, keep_logs=true, backtrack_sorting=false))
    num_colors = 8
    @variable(m, 1 <= max_color <= num_colors, Int)

    @variable(m, 1 <= states[1:50] <= num_colors, Int)
    washington,montana,maine,north_dakota,south_dakota,wyoming,wisconsin,idaho,vermont,minnesota,oregon,new_hampshire,
    iowa,massachusetts,nebraska,new_york,pennsylvania,connecticut,rhode_island,new_jersey,indiana,nevada,utah,california,ohio,
    illinois,washington_dc,delaware,west_virginia,maryland,colorado,kentucky,kansas,virginia,missouri,arizona,oklahoma,north_carolina,
    tennessee,texas,new_mexico,alabama,mississippi,georgia,south_carolina,arkansas,louisiana,florida,michigan = states

    @constraint(m, washington != oregon)
    @constraint(m, washington != idaho)
    @constraint(m, oregon != idaho)
    @constraint(m, oregon != nevada)
    @constraint(m, oregon != california)
    @constraint(m, california != nevada)
    @constraint(m, california != arizona)
    @constraint(m, nevada != idaho)
    @constraint(m, nevada != utah)
    @constraint(m, nevada != arizona)
    @constraint(m, idaho != montana)
    @constraint(m, idaho != wyoming)
    @constraint(m, idaho != utah)
    @constraint(m, utah != wyoming)
    @constraint(m, utah != colorado)
    @constraint(m, utah != new_mexico)
    @constraint(m, utah != arizona)
    @constraint(m, arizona != colorado)
    @constraint(m, arizona != new_mexico)
    @constraint(m, montana != north_dakota)
    @constraint(m, montana != south_dakota)
    @constraint(m, montana != wyoming)
    @constraint(m, wyoming != south_dakota)
    @constraint(m, wyoming != nebraska)
    @constraint(m, wyoming != colorado)
    @constraint(m, colorado != nebraska)
    @constraint(m, colorado != kansas)
    @constraint(m, colorado != oklahoma)
    @constraint(m, colorado != new_mexico)
    @constraint(m, new_mexico != oklahoma)
    @constraint(m, new_mexico != texas)
    @constraint(m, north_dakota != minnesota)
    @constraint(m, north_dakota != south_dakota)
    @constraint(m, south_dakota != minnesota)
    @constraint(m, south_dakota != iowa)
    @constraint(m, south_dakota != nebraska)
    @constraint(m, nebraska != iowa)
    @constraint(m, nebraska != missouri)
    @constraint(m, nebraska != kansas)
    @constraint(m, kansas != missouri)
    @constraint(m, kansas != oklahoma)
    @constraint(m, oklahoma != arkansas)
    @constraint(m, oklahoma != texas)
    @constraint(m, texas != arkansas)
    @constraint(m, texas != louisiana)
    @constraint(m, minnesota != wisconsin)
    @constraint(m, minnesota != iowa)
    @constraint(m, iowa != wisconsin)
    @constraint(m, iowa != illinois)
    @constraint(m, iowa != missouri)
    @constraint(m, missouri != illinois)
    @constraint(m, missouri != kentucky)
    @constraint(m, missouri != tennessee)
    @constraint(m, missouri != arkansas)
    @constraint(m, arkansas != tennessee)
    @constraint(m, arkansas != mississippi)
    @constraint(m, arkansas != louisiana)
    @constraint(m, louisiana != mississippi)
    @constraint(m, wisconsin != illinois)
    @constraint(m, wisconsin != michigan)
    @constraint(m, illinois != indiana)
    @constraint(m, illinois != kentucky)
    @constraint(m, kentucky != indiana)
    @constraint(m, kentucky != ohio)
    @constraint(m, kentucky != west_virginia)
    @constraint(m, kentucky != virginia)
    @constraint(m, kentucky != tennessee)
    @constraint(m, tennessee != virginia)
    @constraint(m, tennessee != north_carolina)
    @constraint(m, tennessee != georgia)
    @constraint(m, tennessee != alabama)
    @constraint(m, tennessee != mississippi)
    @constraint(m, mississippi != alabama)
    @constraint(m, michigan != indiana)
    @constraint(m, michigan != ohio)
    @constraint(m, indiana != ohio)
    @constraint(m, alabama != georgia)
    @constraint(m, alabama != florida)
    @constraint(m, ohio != pennsylvania)
    @constraint(m, ohio != west_virginia)
    @constraint(m, maine != new_hampshire)
    @constraint(m, new_hampshire != vermont)
    @constraint(m, new_hampshire != massachusetts)
    @constraint(m, vermont != massachusetts)
    @constraint(m, vermont != new_york)
    @constraint(m, massachusetts != new_york)
    @constraint(m, massachusetts != rhode_island)
    @constraint(m, massachusetts != connecticut)
    @constraint(m, connecticut != rhode_island)
    @constraint(m, new_york != connecticut)
    @constraint(m, new_york != pennsylvania)
    @constraint(m, new_york != new_jersey)
    @constraint(m, pennsylvania != new_jersey)
    @constraint(m, pennsylvania != delaware)
    @constraint(m, pennsylvania != maryland)
    @constraint(m, pennsylvania != west_virginia)
    @constraint(m, new_jersey != delaware)
    @constraint(m, maryland != washington_dc)
    @constraint(m, maryland != west_virginia)
    @constraint(m, maryland != virginia)
    @constraint(m, washington_dc != virginia)
    @constraint(m, west_virginia != virginia)
    @constraint(m, virginia != north_carolina)
    @constraint(m, north_carolina != south_carolina)
    @constraint(m, south_carolina != georgia)
    @constraint(m, georgia != florida)

    @constraint(m, max_color .>= states)

    @objective(m, Min, max_color)

    optimize!(m)

    status = JuMP.termination_status(m)

    com = JuMP.backend(m).optimizer.model.inner

    CS.save_logs(com, "graph_color_optimize.json")
    rm("graph_color_optimize.json")

    @test status == MOI.OPTIMAL
    @test com.best_sol == 4
    @test all([length(CS.values(m, var)) == 1 for var in states])
    @test maximum([JuMP.value(var) for var in states]) == JuMP.value(max_color) == 4

end

@testset "49 US states + DC only 3 colors" begin
    m = Model(with_optimizer(CS.Optimizer))
    num_colors = 3

    @variable(m, 1 <= states[1:50] <= num_colors, Int)
    washington,montana,maine,north_dakota,south_dakota,wyoming,wisconsin,idaho,vermont,minnesota,oregon,new_hampshire,
    iowa,massachusetts,nebraska,new_york,pennsylvania,connecticut,rhode_island,new_jersey,indiana,nevada,utah,california,ohio,
    illinois,washington_dc,delaware,west_virginia,maryland,colorado,kentucky,kansas,virginia,missouri,arizona,oklahoma,north_carolina,
    tennessee,texas,new_mexico,alabama,mississippi,georgia,south_carolina,arkansas,louisiana,florida,michigan = states

    @constraint(m, washington != oregon)
    @constraint(m, washington != idaho)
    @constraint(m, oregon != idaho)
    @constraint(m, oregon != nevada)
    @constraint(m, oregon != california)
    @constraint(m, california != nevada)
    @constraint(m, california != arizona)
    @constraint(m, nevada != idaho)
    @constraint(m, nevada != utah)
    @constraint(m, nevada != arizona)
    @constraint(m, idaho != montana)
    @constraint(m, idaho != wyoming)
    @constraint(m, idaho != utah)
    @constraint(m, utah != wyoming)
    @constraint(m, utah != colorado)
    @constraint(m, utah != new_mexico)
    @constraint(m, utah != arizona)
    @constraint(m, arizona != colorado)
    @constraint(m, arizona != new_mexico)
    @constraint(m, montana != north_dakota)
    @constraint(m, montana != south_dakota)
    @constraint(m, montana != wyoming)
    @constraint(m, wyoming != south_dakota)
    @constraint(m, wyoming != nebraska)
    @constraint(m, wyoming != colorado)
    @constraint(m, colorado != nebraska)
    @constraint(m, colorado != kansas)
    @constraint(m, colorado != oklahoma)
    @constraint(m, colorado != new_mexico)
    @constraint(m, new_mexico != oklahoma)
    @constraint(m, new_mexico != texas)
    @constraint(m, north_dakota != minnesota)
    @constraint(m, north_dakota != south_dakota)
    @constraint(m, south_dakota != minnesota)
    @constraint(m, south_dakota != iowa)
    @constraint(m, south_dakota != nebraska)
    @constraint(m, nebraska != iowa)
    @constraint(m, nebraska != missouri)
    @constraint(m, nebraska != kansas)
    @constraint(m, kansas != missouri)
    @constraint(m, kansas != oklahoma)
    @constraint(m, oklahoma != arkansas)
    @constraint(m, oklahoma != texas)
    @constraint(m, texas != arkansas)
    @constraint(m, texas != louisiana)
    @constraint(m, minnesota != wisconsin)
    @constraint(m, minnesota != iowa)
    @constraint(m, iowa != wisconsin)
    @constraint(m, iowa != illinois)
    @constraint(m, iowa != missouri)
    @constraint(m, missouri != illinois)
    @constraint(m, missouri != kentucky)
    @constraint(m, missouri != tennessee)
    @constraint(m, missouri != arkansas)
    @constraint(m, arkansas != tennessee)
    @constraint(m, arkansas != mississippi)
    @constraint(m, arkansas != louisiana)
    @constraint(m, louisiana != mississippi)
    @constraint(m, wisconsin != illinois)
    @constraint(m, wisconsin != michigan)
    @constraint(m, illinois != indiana)
    @constraint(m, illinois != kentucky)
    @constraint(m, kentucky != indiana)
    @constraint(m, kentucky != ohio)
    @constraint(m, kentucky != west_virginia)
    @constraint(m, kentucky != virginia)
    @constraint(m, kentucky != tennessee)
    @constraint(m, tennessee != virginia)
    @constraint(m, tennessee != north_carolina)
    @constraint(m, tennessee != georgia)
    @constraint(m, tennessee != alabama)
    @constraint(m, tennessee != mississippi)
    @constraint(m, mississippi != alabama)
    @constraint(m, michigan != indiana)
    @constraint(m, michigan != ohio)
    @constraint(m, indiana != ohio)
    @constraint(m, alabama != georgia)
    @constraint(m, alabama != florida)
    @constraint(m, ohio != pennsylvania)
    @constraint(m, ohio != west_virginia)
    @constraint(m, maine != new_hampshire)
    @constraint(m, new_hampshire != vermont)
    @constraint(m, new_hampshire != massachusetts)
    @constraint(m, vermont != massachusetts)
    @constraint(m, vermont != new_york)
    @constraint(m, massachusetts != new_york)
    @constraint(m, massachusetts != rhode_island)
    @constraint(m, massachusetts != connecticut)
    @constraint(m, connecticut != rhode_island)
    @constraint(m, new_york != connecticut)
    @constraint(m, new_york != pennsylvania)
    @constraint(m, new_york != new_jersey)
    @constraint(m, pennsylvania != new_jersey)
    @constraint(m, pennsylvania != delaware)
    @constraint(m, pennsylvania != maryland)
    @constraint(m, pennsylvania != west_virginia)
    @constraint(m, new_jersey != delaware)
    @constraint(m, maryland != washington_dc)
    @constraint(m, maryland != west_virginia)
    @constraint(m, maryland != virginia)
    @constraint(m, washington_dc != virginia)
    @constraint(m, west_virginia != virginia)
    @constraint(m, virginia != north_carolina)
    @constraint(m, north_carolina != south_carolina)
    @constraint(m, south_carolina != georgia)
    @constraint(m, georgia != florida)

    optimize!(m)

    @test JuMP.termination_status(m) == MOI.INFEASIBLE
end

@testset "Maximization objective" begin
    m = Model(with_optimizer(CS.Optimizer))
    num_colors = 20

    @variable(m, 1 <= states[1:50] <= num_colors, Int)
    @variable(m, 1 <= max_color <= num_colors, Int)

    washington,montana,maine,north_dakota,south_dakota,wyoming,wisconsin,idaho,vermont,minnesota,oregon,new_hampshire,
    iowa,massachusetts,nebraska,new_york,pennsylvania,connecticut,rhode_island,new_jersey,indiana,nevada,utah,california,ohio,
    illinois,washington_dc,delaware,west_virginia,maryland,colorado,kentucky,kansas,virginia,missouri,arizona,oklahoma,north_carolina,
    tennessee,texas,new_mexico,alabama,mississippi,georgia,south_carolina,arkansas,louisiana,florida,michigan = states

    @constraint(m, washington != oregon)
    @constraint(m, washington != idaho)
    @constraint(m, oregon != idaho)
    @constraint(m, oregon != nevada)
    @constraint(m, oregon != california)
    @constraint(m, california != nevada)
    @constraint(m, california != arizona)
    @constraint(m, nevada != idaho)
    @constraint(m, nevada != utah)
    @constraint(m, nevada != arizona)
    @constraint(m, idaho != montana)
    @constraint(m, idaho != wyoming)
    @constraint(m, idaho != utah)
    @constraint(m, utah != wyoming)
    @constraint(m, utah != colorado)
    @constraint(m, utah != new_mexico)
    @constraint(m, utah != arizona)
    @constraint(m, arizona != colorado)
    @constraint(m, arizona != new_mexico)
    @constraint(m, montana != north_dakota)
    @constraint(m, montana != south_dakota)
    @constraint(m, montana != wyoming)
    @constraint(m, wyoming != south_dakota)
    @constraint(m, wyoming != nebraska)
    @constraint(m, wyoming != colorado)
    @constraint(m, colorado != nebraska)
    @constraint(m, colorado != kansas)
    @constraint(m, colorado != oklahoma)
    @constraint(m, colorado != new_mexico)
    @constraint(m, new_mexico != oklahoma)
    @constraint(m, new_mexico != texas)
    @constraint(m, north_dakota != minnesota)
    @constraint(m, north_dakota != south_dakota)
    @constraint(m, south_dakota != minnesota)
    @constraint(m, south_dakota != iowa)
    @constraint(m, south_dakota != nebraska)
    @constraint(m, nebraska != iowa)
    @constraint(m, nebraska != missouri)
    @constraint(m, nebraska != kansas)
    @constraint(m, kansas != missouri)
    @constraint(m, kansas != oklahoma)
    @constraint(m, oklahoma != arkansas)
    @constraint(m, oklahoma != texas)
    @constraint(m, texas != arkansas)
    @constraint(m, texas != louisiana)
    @constraint(m, minnesota != wisconsin)
    @constraint(m, minnesota != iowa)
    @constraint(m, iowa != wisconsin)
    @constraint(m, iowa != illinois)
    @constraint(m, iowa != missouri)
    @constraint(m, missouri != illinois)
    @constraint(m, missouri != kentucky)
    @constraint(m, missouri != tennessee)
    @constraint(m, missouri != arkansas)
    @constraint(m, arkansas != tennessee)
    @constraint(m, arkansas != mississippi)
    @constraint(m, arkansas != louisiana)
    @constraint(m, louisiana != mississippi)
    @constraint(m, wisconsin != illinois)
    @constraint(m, wisconsin != michigan)
    @constraint(m, illinois != indiana)
    @constraint(m, illinois != kentucky)
    @constraint(m, kentucky != indiana)
    @constraint(m, kentucky != ohio)
    @constraint(m, kentucky != west_virginia)
    @constraint(m, kentucky != virginia)
    @constraint(m, kentucky != tennessee)
    @constraint(m, tennessee != virginia)
    @constraint(m, tennessee != north_carolina)
    @constraint(m, tennessee != georgia)
    @constraint(m, tennessee != alabama)
    @constraint(m, tennessee != mississippi)
    @constraint(m, mississippi != alabama)
    @constraint(m, michigan != indiana)
    @constraint(m, michigan != ohio)
    @constraint(m, indiana != ohio)
    @constraint(m, alabama != georgia)
    @constraint(m, alabama != florida)
    @constraint(m, ohio != pennsylvania)
    @constraint(m, ohio != west_virginia)
    @constraint(m, maine != new_hampshire)
    @constraint(m, new_hampshire != vermont)
    @constraint(m, new_hampshire != massachusetts)
    @constraint(m, vermont != massachusetts)
    @constraint(m, vermont != new_york)
    @constraint(m, massachusetts != new_york)
    @constraint(m, massachusetts != rhode_island)
    @constraint(m, massachusetts != connecticut)
    @constraint(m, connecticut != rhode_island)
    @constraint(m, new_york != connecticut)
    @constraint(m, new_york != pennsylvania)
    @constraint(m, new_york != new_jersey)
    @constraint(m, pennsylvania != new_jersey)
    @constraint(m, pennsylvania != delaware)
    @constraint(m, pennsylvania != maryland)
    @constraint(m, pennsylvania != west_virginia)
    @constraint(m, new_jersey != delaware)
    @constraint(m, maryland != washington_dc)
    @constraint(m, maryland != west_virginia)
    @constraint(m, maryland != virginia)
    @constraint(m, washington_dc != virginia)
    @constraint(m, west_virginia != virginia)
    @constraint(m, virginia != north_carolina)
    @constraint(m, north_carolina != south_carolina)
    @constraint(m, south_carolina != georgia)
    @constraint(m, georgia != florida)

    @constraint(m, max_color .<= states)

    @objective(m, Max, max_color)

    optimize!(m)

    status = JuMP.termination_status(m)

    com = JuMP.backend(m).optimizer.model.inner

    @test status == MOI.OPTIMAL
    @test com.best_sol == 17
    @test minimum([JuMP.value(var) for var in states]) == 17 == JuMP.value(max_color)
end

end


text:
using QuasiGeostrophy
using LinearAlgebra, Test, FFTW, BenchmarkTools, Plots

# boiler plate definitions
Ωxy = S¹(0,2π) × S¹(0,2π)
Nx = 2^5; Ny = 2^5;
fourier_grid = FourierGrid((Nx, Ny), Ωxy)
x, y = fourier_grid.grid
kx, ky = fourier_grid.wavenumbers
ψ1 = @. 0 * sin(x) - sin(y) + 0im
ψ2 = @. 0 * sin(x) - sin(y) + 0im
ψ̂1 = copy(ψ1)
ψ̂2 = copy(ψ2)
q̂1 = copy(ψ2)
q̂2 = copy(ψ2)
q1 = copy(ψ2)
q2 = copy(ψ2)
FFTW.set_num_threads(Threads.nthreads())
P = plan_fft(ψ1) # typeof(iP) <: AbstractFFTs.ScaledPlan
iP = plan_ifft(ψ1)
mul!(ψ̂1, P, ψ1)

## Define Parameters
const λ  = 1e-2
const U  = 1e-2
const ν  = 1e-6
const κ  = 1e-0
const Δt = 1e-1
const L  = 1.0
const Q  = 1.0

## Define forcing
forcing = @. 0 * sin(x) + Q * sin(y/L) + 0im
f̂ = copy(forcing)
mul!(f̂, P, forcing)

## Define Operators
∂x = FourierOperator(im .* kx)
∂y = FourierOperator(im .* ky)
Δ = ∂x^2 + ∂y^2
Δ⁵ = Δ^5
G = inv( (Δ + -1 /(2*λ^2) )^2 + -(1/(2*λ^2))^2 )

function get_pv!(q̂1, q̂2, ψ̂1, ψ̂2)
    tmp1 = Δ(ψ̂1) + 1 /(2*λ^2) * (ψ̂2 - ψ̂1)
    q̂1 .= tmp1
    tmp2 = Δ(ψ̂2) + 1 /(2*λ^2) * (ψ̂1 - ψ̂2)
    q̂2 .= tmp2
    return nothing
end

get_pv!(q̂1, q̂2, ψ̂1, ψ̂2)

function get_stream_function!(q1, q2, ψ1, ψ2)
    tmp1 = G(Δ(q1) + 1 /(2*λ^2) * (-q2 - q1))
    ψ1 .= tmp1
    tmp2 = G(Δ(q2) + 1 /(2*λ^2) * (-q1 - q2))
    ψ2 .= tmp2
    return nothing
end

function closure_J(∂x, ∂y, P, iP)
    function J(a,b)
        term_1 = filter_convolve(∂x(a) , ∂y(b), P, iP)
        term_2 = filter_convolve(∂x(b) , ∂y(a), P, iP)
        return  term_1 - term_2
    end
end

J = closure_J(∂x, ∂y, P, iP)

function closure_rhs(∂x, ∂y, Δ, Δ⁵, J, λ, ν, U, κ, forcing)
    function rhs(state)
        q1 = state[1]
        q2 = state[2]
        ψ1 = copy(q1)
        ψ2 = copy(q2)
        get_stream_function!(q1, q2, ψ1, ψ2)
        drag = - 2 .* κ .* (Δ(ψ2))
        q1_rhs = -J(ψ1, q1) - (U/λ^2) .* ∂x(ψ1) - U .* ∂x(q1) + forcing
        q2_rhs = -J(ψ2, q2) + (U/λ^2) .* ∂x(ψ2) + U .* ∂x(q2) - forcing + drag
        return [q1_rhs, q2_rhs]
    end
end

rhs = closure_rhs(∂x, ∂y, Δ, Δ⁵, J, λ, ν, U, κ, f̂)

function closure_lhs(∂x, ∂y, Δ, Δ⁵, J, λ, ν, U, Δt)
    function lhs()
        q1_lhs = 1 + Δt * ν * Δ^4
        return [q1_lhs, q1_lhs]
    end
end

lhs = closure_lhs(∂x, ∂y, Δ, Δ⁵, J, λ, ν, U, Δt)
operators = lhs()
lhs_operator = [inv(operator) for operator in operators]

function imex_step!(state, lhs_operator, rhs, Δt)
    rhs_state = rhs(state)
    current_state = copy(state)
    for i in eachindex(state)
        tmp = lhs_operator[i](current_state[i] .+ Δt .* rhs_state[i])
        @. state[i] =  tmp
    end
end
# @btime imex_step!(state, lhs_operator, rhs, Δt);

##
q1 = @. 1 * sin(x) - cos(y) + 0im
q2 = @. 1 * sin(x) + sin(y) + 0im
fq1 = P * q1
fq2 = P * q2
τ = (q1 - q2) ./ 2
τ̂ = (fq1 - fq2) ./ 2
tmp = P * (400 * sin.(x) .+ 0 * sin.(y))
#fq1 .= good_start1
#fq2 .= good_start2
state = [fq1, fq2]
##
gr(size = (400,400))
for i in 1:2000
    imex_step!(state, lhs_operator, rhs, Δt)
    if (i%100) == 0
        @. τ̂ = (fq1 - fq2) / 2
        mul!(τ, iP, τ̂)
        local p1 = contourf(x[:], y[:], real.(τ)', linewidth = 0, color = :thermometer)
        local p2 = plot(sum(real.(τ), dims = 1)[:] ./ Nx, x[:] )
        display(plot(p1,p2))
        sleep(0.01)
    end
end
#good_start1 = copy(fq1)
#good_start2 = copy(fq2)
##

@. τ̂ = (fq1 - fq2) / 2
mul!(τ, iP, τ̂)

##

# Check PV and Stream Function Calculation
tmp = @. 1 * sin(2*x) + Q * sin(y/L * 3) + 0im
tmp2 = @. 1 * sin(x) - Q * sin(y/L) + 0im
ψ̂1 = P * tmp
ψ̂2 = P * tmp2
q1 = copy(q̂1)
q2 = copy(q̂2)
ψ1 = copy(ψ̂1)
ψ2 = copy(ψ̂2)

get_pv!(q̂1, q̂2, ψ̂1, ψ̂2)
get_stream_function!(q̂1, q̂2, ψ̂1, ψ̂2)
norm(ψ1 - ψ̂1)
norm(ψ2 - ψ̂2)
##
i = 1
@. state[i] *= lhs_operator[i].op

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
"""
    find_best_item(examinee, items)

# Description
Finds the best item for the `examinee` among the vector `items` using maximum information criterion.

# Arguments
  - **`examinee::AbstractExaminee`**: The examinee.
  - **`items:::Vector{<:AbstractItem}`**: Set of items in which searching for the best.
  - **`method`**: Optional. "max-info" | "D-gain".
  max-info takes the item with the maximum latent information,
  D-gain takes the item with the maximum gain in the determinant of the expected information matrix if examinee would answer to the item (D-VC in __Ren, van der Linden, Diao, 2017__).

# Output
It returns an item of generic type.

# References
[^RenDiao2017]

[^RenDiao2017]:Ren H, van der Linden WJ, Diao Q. Continuous online item calibration: Parameter recovery and item calibration. Psychometrika. 2017;82:498–522. doi: 10.1007/s11336-017-9553-1.
"""
function find_best_item(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem};
    method = "max-info",
)
    if method in ["max-info", "D-gain"]
        if method == "max-info"
            infos = latent_information(examinee, items)
        elseif method == "D-gain"
            infos = map(i -> D_gain_method(i, examinee), items)
        end
        return items[findmax(infos)[2]].idx::Int64
    else
        error("only max-info and D-gain optimality are available.")
    end
end

"""
    find_best_examinee(item, examinees; method = "D")

# Description
Finds the best examinee among the `examinees` vector, for the `item` using maximum expected information criterion.

# Arguments
- **`item::AbstractItem`**: Required. Set of items in which searching for the best.
- **`examinees::Vector{<:AbstractExaminee}`**: Required. The examinee.
- **`method`**: Optional. "D" | "A".
  D stands for D-optimality (determinant of the expected information matrix),
  A stands for A-optimality (trace of the expected information matrix).

# Output
It returns the idx of the best examinee.

"""
function find_best_examinee(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee};
    method = "D",
)
    if method in ["D", "A"]
        if method == "D"
            infos = map(e -> D_method(item, e), examinees)
        elseif method == "A"
            infos = map(e -> A_method(item, e), examinees)
        end
        return examinees[findmax(infos)[2]].idx::Int64
    else
        error("only D-optimalit and A-optimality are available.")
    end
end


text:
struct MockField{T,D} <: Field
  v::T
  function MockField{D}(v::Number) where {T,D}
    new{typeof(v),D}(v)
  end
end

function field_cache(f::MockField,x)
  nx = length(x)
  c = zeros(typeof(f.v),nx)
  CachedArray(c)
end

function evaluate_field!(c,f::MockField,x)
  nx = length(x)
  setsize!(c,(nx,))
  for i in eachindex(x)
    @inbounds xi = x[i]
    @inbounds c[i] = f.v*xi[1]
  end
  c
end

@inline function field_gradient(f::MockField{T,D}) where {T,D}
  E = eltype(T)
  P = Point{D,E}
  _p = zero(mutable(P))
  _p[1] = one(E)
  p = Point(_p)
  vg = outer(p,f.v)
  MockField{D}(vg)
end

struct MockBasis{V,D} <: Field
  v::V
  ndofs::Int
  function MockBasis{D}(v::Number,ndofs::Int) where D
    new{typeof(v),D}(v,ndofs)
  end
end

function field_cache(f::MockBasis,x)
  np = length(x)
  s = (np, f.ndofs)
  c = zeros(typeof(f.v),s)
  CachedArray(c)
end

function evaluate_field!(v,f::MockBasis,x)
  np = length(x)
  s = (np, f.ndofs)
  setsize!(v,s)
  for i in 1:np
    @inbounds xi = x[i]
    for j in 1:f.ndofs
      @inbounds v[i,j] = f.v*xi[1]
    end
  end
  v
end

@inline function field_gradient(f::MockBasis{T,D}) where {T,D}
  E = eltype(T)
  P = Point{D,E}
  _p = zero(mutable(P))
  _p[1] = one(E)
  p = Point(_p)
  vg = outer(p,f.v)
  MockBasis{D}(vg,f.ndofs)
end

struct OtherMockBasis{D} <: Field
  ndofs::Int
  function OtherMockBasis{D}(ndofs::Int) where D
    new{D}(ndofs)
  end
end

function field_cache(f::OtherMockBasis,x)
  np = length(x)
  s = (np, f.ndofs)
  c = zeros(eltype(x),s)
  CachedArray(c)
end

function evaluate_field!(v,f::OtherMockBasis,x)
  np = length(x)
  s = (np, f.ndofs)
  setsize!(v,s)
  for i in 1:np
    @inbounds xi = x[i]
    for j in 1:f.ndofs
      @inbounds v[i,j] = 2*xi
    end
  end
  v
end

@inline function field_gradient(f::OtherMockBasis{D}) where D
  E = Float64
  P = Point{D,E}
  p = zero(P)
  vg = 2*one(outer(p,p))
  MockBasis{D}(vg,f.ndofs)
end


text:
module Spectra

export AbstractMatOp, DenseMatProd
export symeigs

include("hessenqr.jl")
include("symeigs.jl")

end # module


text:
@testset "4.3.10 (c+d x)^m (a+b tan)^n" begin
    (a, b, c, d, e, f, m, x, ) = @variables a b c d e f m x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*tan(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(b*tan(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(b*tan(e+f*x))^n=#

    @test_int [x^3*tan(a + b*x), x, 6, (im*x^4)/4 - (x^3*log(1 + ℯ^(2*im*(a + b*x))))/b + (3*im*x^2*PolyLog(2, -ℯ^(2*im*(a + b*x))))/(2*b^2) - (3*x*PolyLog(3, -ℯ^(2*im*(a + b*x))))/(2*b^3) - (3*im*PolyLog(4, -ℯ^(2*im*(a + b*x))))/(4*b^4)]
    @test_int [x^2*tan(a + b*x), x, 5, (im*x^3)/3 - (x^2*log(1 + ℯ^(2*im*(a + b*x))))/b + (im*x*PolyLog(2, -ℯ^(2*im*(a + b*x))))/b^2 - PolyLog(3, -ℯ^(2*im*(a + b*x)))/(2*b^3)]
    @test_int [x^1*tan(a + b*x), x, 4, (im*x^2)/2 - (x*log(1 + ℯ^(2*im*(a + b*x))))/b + (im*PolyLog(2, -ℯ^(2*im*(a + b*x))))/(2*b^2)]
    @test_int [tan(a + b*x)/x^1, x, 0, Unintegrable(tan(a + b*x)/x, x)]
    @test_int [tan(a + b*x)/x^2, x, 0, Unintegrable(tan(a + b*x)/x^2, x)]

    @test_int [x^3*tan(a + b*x)^2, x, 7, -((im*x^3)/b) - x^4/4 + (3*x^2*log(1 + ℯ^(2*im*(a + b*x))))/b^2 - (3*im*x*PolyLog(2, -ℯ^(2*im*(a + b*x))))/b^3 + (3*PolyLog(3, -ℯ^(2*im*(a + b*x))))/(2*b^4) + (x^3*tan(a + b*x))/b]
    @test_int [x^2*tan(a + b*x)^2, x, 6, -((im*x^2)/b) - x^3/3 + (2*x*log(1 + ℯ^(2*im*(a + b*x))))/b^2 - (im*PolyLog(2, -ℯ^(2*im*(a + b*x))))/b^3 + (x^2*tan(a + b*x))/b]
    @test_int [x^1*tan(a + b*x)^2, x, 3, -x^2/2 + log(cos(a + b*x))/b^2 + (x*tan(a + b*x))/b]
    @test_int [tan(a + b*x)^2/x^1, x, 0, Unintegrable(tan(a + b*x)^2/x, x)]
    @test_int [tan(a + b*x)^2/x^2, x, 0, Unintegrable(tan(a + b*x)^2/x^2, x)]

    @test_int [x^3*tan(a + b*x)^3, x, 13, (3*im*x^2)/(2*b^2) + x^3/(2*b) - (im*x^4)/4 - (3*x*log(1 + ℯ^(2*im*(a + b*x))))/b^3 + (x^3*log(1 + ℯ^(2*im*(a + b*x))))/b + (3*im*PolyLog(2, -ℯ^(2*im*(a + b*x))))/(2*b^4) - (3*im*x^2*PolyLog(2, -ℯ^(2*im*(a + b*x))))/(2*b^2) + (3*x*PolyLog(3, -ℯ^(2*im*(a + b*x))))/(2*b^3) + (3*im*PolyLog(4, -ℯ^(2*im*(a + b*x))))/(4*b^4) - (3*x^2*tan(a + b*x))/(2*b^2) + (x^3*tan(a + b*x)^2)/(2*b)]
    @test_int [x^2*tan(a + b*x)^3, x, 9, x^2/(2*b) - (im*x^3)/3 + (x^2*log(1 + ℯ^(2*im*(a + b*x))))/b - log(cos(a + b*x))/b^3 - (im*x*PolyLog(2, -ℯ^(2*im*(a + b*x))))/b^2 + PolyLog(3, -ℯ^(2*im*(a + b*x)))/(2*b^3) - (x*tan(a + b*x))/b^2 + (x^2*tan(a + b*x)^2)/(2*b)]
    @test_int [x^1*tan(a + b*x)^3, x, 7, x/(2*b) - (im*x^2)/2 + (x*log(1 + ℯ^(2*im*(a + b*x))))/b - (im*PolyLog(2, -ℯ^(2*im*(a + b*x))))/(2*b^2) - tan(a + b*x)/(2*b^2) + (x*tan(a + b*x)^2)/(2*b)]
    @test_int [tan(a + b*x)^3/x^1, x, 0, Unintegrable(tan(a + b*x)^3/x, x)]
    @test_int [tan(a + b*x)^3/x^2, x, 0, Unintegrable(tan(a + b*x)^3/x^2, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(b*tan(e+f*x))^(n/2)=#

    @test_int [x^2/tan(a + b*x)^(3/2) - (4*x)/(b*sqrt(tan(a + b*x))) + x^2*sqrt(tan(a + b*x)), x, 76, -((2*x^2)/(b*sqrt(tan(a + b*x))))]
    @test_int [x^2*tan(a + b*x^2)^(3/2) + x^2/sqrt(tan(a + b*x^2)) + sqrt(tan(a + b*x^2))/b, x, -1, (x*sqrt(tan(a + b*x^2)))/b]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*tan(e+f*x))^n*with*a^2+b^2=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+im*a*tan(e+f*x))^n=#

    #= ::Subsubsection:: =#
    #=n>0=#

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [(c + d*x)^3/(a + im*a*tan(e + f*x)), x, 5, (((3*im)/8)*d^3*x)/(a*f^3) - (3*d*(c + d*x)^2)/(8*a*f^2) - ((im/4)*(c + d*x)^3)/(a*f) + (c + d*x)^4/(8*a*d) - (3*d^3)/(8*f^4*(a + im*a*tan(e + f*x))) - (((3*im)/4)*d^2*(c + d*x))/(f^3*(a + im*a*tan(e + f*x))) + (3*d*(c + d*x)^2)/(4*f^2*(a + im*a*tan(e + f*x))) + ((im/2)*(c + d*x)^3)/(f*(a + im*a*tan(e + f*x)))]
    @test_int [(c + d*x)^2/(a + im*a*tan(e + f*x)), x, 4, -(d^2*x)/(4*a*f^2) - ((im/4)*(c + d*x)^2)/(a*f) + (c + d*x)^3/(6*a*d) - ((im/4)*d^2)/(f^3*(a + im*a*tan(e + f*x))) + (d*(c + d*x))/(2*f^2*(a + im*a*tan(e + f*x))) + ((im/2)*(c + d*x)^2)/(f*(a + im*a*tan(e + f*x)))]
    @test_int [(c + d*x)^1/(a + im*a*tan(e + f*x)), x, 3, ((-im/4)*d*x)/(a*f) + (c + d*x)^2/(4*a*d) + d/(4*f^2*(a + im*a*tan(e + f*x))) + ((im/2)*(c + d*x))/(f*(a + im*a*tan(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + im*a*tan(e + f*x))), x, 7, (cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(2*a*d) + log(c + d*x)/(2*a*d) - ((im/2)*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(a*d) - ((im/2)*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a*d) - (sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(2*a*d)]
    @test_int [1/((c + d*x)^2*(a + im*a*tan(e + f*x))), x, 7, ((-im)*f*cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(a*d^2) - (f*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(a*d^2) - (f*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a*d^2) + (im*f*sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a*d^2) - 1/(d*(c + d*x)*(a + im*a*tan(e + f*x)))]
    @test_int [1/((c + d*x)^3*(a + im*a*tan(e + f*x))), x, 8, ((-im/2)*f)/(a*d^2*(c + d*x)) - (f^2*cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(a*d^3) + (im*f^2*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(a*d^3) + (im*f^2*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a*d^3) + (f^2*sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a*d^3) - 1/(2*d*(c + d*x)^2*(a + im*a*tan(e + f*x))) + (im*f)/(d^2*(c + d*x)*(a + im*a*tan(e + f*x)))]

    @test_int [(c + d*x)^3/(a + im*a*tan(e + f*x))^2, x, 10, (-3*d^3*ℯ^((-2*im)*e - (2*im)*f*x))/(16*a^2*f^4) - (3*d^3*ℯ^((-4*im)*e - (4*im)*f*x))/(512*a^2*f^4) - (((3*im)/8)*d^2*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x))/(a^2*f^3) - (((3*im)/128)*d^2*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x))/(a^2*f^3) + (3*d*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x)^2)/(8*a^2*f^2) + (3*d*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x)^2)/(64*a^2*f^2) + ((im/4)*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x)^3)/(a^2*f) + ((im/16)*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x)^3)/(a^2*f) + (c + d*x)^4/(16*a^2*d)]
    @test_int [(c + d*x)^2/(a + im*a*tan(e + f*x))^2, x, 8, ((-im/8)*d^2*ℯ^((-2*im)*e - (2*im)*f*x))/(a^2*f^3) - ((im/128)*d^2*ℯ^((-4*im)*e - (4*im)*f*x))/(a^2*f^3) + (d*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x))/(4*a^2*f^2) + (d*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x))/(32*a^2*f^2) + ((im/4)*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x)^2)/(a^2*f) + ((im/16)*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x)^2)/(a^2*f) + (c + d*x)^3/(12*a^2*d)]
    @test_int [(c + d*x)^1/(a + im*a*tan(e + f*x))^2, x, 7, (((-3*im)/16)*d*x)/(a^2*f) - (d*x^2)/(8*a^2) + (x*(c + d*x))/(4*a^2) + d/(16*f^2*(a + im*a*tan(e + f*x))^2) + ((im/4)*(c + d*x))/(f*(a + im*a*tan(e + f*x))^2) + (3*d)/(16*f^2*(a^2 + im*a^2*tan(e + f*x))) + ((im/4)*(c + d*x))/(f*(a^2 + im*a^2*tan(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + im*a*tan(e + f*x))^2), x, 21, (cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(2*a^2*d) + (cos(4*e - (4*c*f)/d)*CosIntegral((4*c*f)/d + 4*f*x))/(4*a^2*d) + log(c + d*x)/(4*a^2*d) - ((im/4)*CosIntegral((4*c*f)/d + 4*f*x)*sin(4*e - (4*c*f)/d))/(a^2*d) - ((im/2)*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(a^2*d) - ((im/2)*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a^2*d) - (sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(2*a^2*d) - ((im/4)*cos(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(a^2*d) - (sin(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(4*a^2*d)]
    @test_int [1/((c + d*x)^2*(a + im*a*tan(e + f*x))^2), x, 24, -(1/(4*a^2*d*(c + d*x))) - cos(2*e + 2*f*x)/(2*a^2*d*(c + d*x)) - cos(2*e + 2*f*x)^2/(4*a^2*d*(c + d*x)) - (im*f*cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(a^2*d^2) - (im*f*cos(4*e - (4*c*f)/d)*CosIntegral((4*c*f)/d + 4*f*x))/(a^2*d^2) - (f*CosIntegral((4*c*f)/d + 4*f*x)*sin(4*e - (4*c*f)/d))/(a^2*d^2) - (f*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(a^2*d^2) + (im*sin(2*e + 2*f*x))/(2*a^2*d*(c + d*x)) + sin(2*e + 2*f*x)^2/(4*a^2*d*(c + d*x)) + (im*sin(4*e + 4*f*x))/(4*a^2*d*(c + d*x)) - (f*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a^2*d^2) + (im*f*sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a^2*d^2) - (f*cos(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(a^2*d^2) + (im*f*sin(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(a^2*d^2)]

    @test_int [(c + d*x)^3/(a + im*a*tan(e + f*x))^3, x, 14, (-9*d^3*ℯ^((-2*im)*e - (2*im)*f*x))/(64*a^3*f^4) - (9*d^3*ℯ^((-4*im)*e - (4*im)*f*x))/(1024*a^3*f^4) - (d^3*ℯ^((-6*im)*e - (6*im)*f*x))/(1728*a^3*f^4) - (((9*im)/32)*d^2*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x))/(a^3*f^3) - (((9*im)/256)*d^2*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x))/(a^3*f^3) - ((im/288)*d^2*ℯ^((-6*im)*e - (6*im)*f*x)*(c + d*x))/(a^3*f^3) + (9*d*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x)^2)/(32*a^3*f^2) + (9*d*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x)^2)/(128*a^3*f^2) + (d*ℯ^((-6*im)*e - (6*im)*f*x)*(c + d*x)^2)/(96*a^3*f^2) + (((3*im)/16)*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x)^3)/(a^3*f) + (((3*im)/32)*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x)^3)/(a^3*f) + ((im/48)*ℯ^((-6*im)*e - (6*im)*f*x)*(c + d*x)^3)/(a^3*f) + (c + d*x)^4/(32*a^3*d)]
    @test_int [(c + d*x)^2/(a + im*a*tan(e + f*x))^3, x, 11, (((-3*im)/32)*d^2*ℯ^((-2*im)*e - (2*im)*f*x))/(a^3*f^3) - (((3*im)/256)*d^2*ℯ^((-4*im)*e - (4*im)*f*x))/(a^3*f^3) - ((im/864)*d^2*ℯ^((-6*im)*e - (6*im)*f*x))/(a^3*f^3) + (3*d*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x))/(16*a^3*f^2) + (3*d*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x))/(64*a^3*f^2) + (d*ℯ^((-6*im)*e - (6*im)*f*x)*(c + d*x))/(144*a^3*f^2) + (((3*im)/16)*ℯ^((-2*im)*e - (2*im)*f*x)*(c + d*x)^2)/(a^3*f) + (((3*im)/32)*ℯ^((-4*im)*e - (4*im)*f*x)*(c + d*x)^2)/(a^3*f) + ((im/48)*ℯ^((-6*im)*e - (6*im)*f*x)*(c + d*x)^2)/(a^3*f) + (c + d*x)^3/(24*a^3*d)]
    @test_int [(c + d*x)^1/(a + im*a*tan(e + f*x))^3, x, 11, (((-11*im)/96)*d*x)/(a^3*f) - (d*x^2)/(16*a^3) + (x*(c + d*x))/(8*a^3) + d/(36*f^2*(a + im*a*tan(e + f*x))^3) + ((im/6)*(c + d*x))/(f*(a + im*a*tan(e + f*x))^3) + (5*d)/(96*a*f^2*(a + im*a*tan(e + f*x))^2) + ((im/8)*(c + d*x))/(a*f*(a + im*a*tan(e + f*x))^2) + (11*d)/(96*f^2*(a^3 + im*a^3*tan(e + f*x))) + ((im/8)*(c + d*x))/(f*(a^3 + im*a^3*tan(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + im*a*tan(e + f*x))^3), x, 53, (3*cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(8*a^3*d) + (3*cos(4*e - (4*c*f)/d)*CosIntegral((4*c*f)/d + 4*f*x))/(8*a^3*d) + (cos(6*e - (6*c*f)/d)*CosIntegral((6*c*f)/d + 6*f*x))/(8*a^3*d) + log(c + d*x)/(8*a^3*d) - ((im/8)*CosIntegral((6*c*f)/d + 6*f*x)*sin(6*e - (6*c*f)/d))/(a^3*d) - (((3*im)/8)*CosIntegral((4*c*f)/d + 4*f*x)*sin(4*e - (4*c*f)/d))/(a^3*d) - (((3*im)/8)*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(a^3*d) - (((3*im)/8)*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(a^3*d) - (3*sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(8*a^3*d) - (((3*im)/8)*cos(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(a^3*d) - (3*sin(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(8*a^3*d) - ((im/8)*cos(6*e - (6*c*f)/d)*SinIntegral((6*c*f)/d + 6*f*x))/(a^3*d) - (sin(6*e - (6*c*f)/d)*SinIntegral((6*c*f)/d + 6*f*x))/(8*a^3*d)]
    @test_int [1/((c + d*x)^2*(a + im*a*tan(e + f*x))^3), x, 60, -(1/(8*a^3*d*(c + d*x))) - (9*cos(2*e + 2*f*x))/(32*a^3*d*(c + d*x)) - (3*cos(2*e + 2*f*x)^2)/(8*a^3*d*(c + d*x)) - cos(2*e + 2*f*x)^3/(8*a^3*d*(c + d*x)) - (3*cos(6*e + 6*f*x))/(32*a^3*d*(c + d*x)) - (3*im*f*cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(4*a^3*d^2) - (3*im*f*cos(4*e - (4*c*f)/d)*CosIntegral((4*c*f)/d + 4*f*x))/(2*a^3*d^2) - (3*im*f*cos(6*e - (6*c*f)/d)*CosIntegral((6*c*f)/d + 6*f*x))/(4*a^3*d^2) - (3*f*CosIntegral((6*c*f)/d + 6*f*x)*sin(6*e - (6*c*f)/d))/(4*a^3*d^2) - (3*f*CosIntegral((4*c*f)/d + 4*f*x)*sin(4*e - (4*c*f)/d))/(2*a^3*d^2) - (3*f*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/(4*a^3*d^2) + (15*im*sin(2*e + 2*f*x))/(32*a^3*d*(c + d*x)) + (3*sin(2*e + 2*f*x)^2)/(8*a^3*d*(c + d*x)) - (im*sin(2*e + 2*f*x)^3)/(8*a^3*d*(c + d*x)) + (3*im*sin(4*e + 4*f*x))/(8*a^3*d*(c + d*x)) + (3*im*sin(6*e + 6*f*x))/(32*a^3*d*(c + d*x)) - (3*f*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(4*a^3*d^2) + (3*im*f*sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(4*a^3*d^2) - (3*f*cos(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(2*a^3*d^2) + (3*im*f*sin(4*e - (4*c*f)/d)*SinIntegral((4*c*f)/d + 4*f*x))/(2*a^3*d^2) - (3*f*cos(6*e - (6*c*f)/d)*SinIntegral((6*c*f)/d + 6*f*x))/(4*a^3*d^2) + (3*im*f*sin(6*e - (6*c*f)/d)*SinIntegral((6*c*f)/d + 6*f*x))/(4*a^3*d^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*tan(e+f*x))^n*with*m*symbolic=#

    @test_int [(c + d*x)^m*(a + im*a*tan(e + f*x))^2, x, 0, Unintegrable((c + d*x)^m*(a + im*a*tan(e + f*x))^2, x)]
    @test_int [(c + d*x)^m*(a + im*a*tan(e + f*x))^1, x, 0, Unintegrable((c + d*x)^m*(a + im*a*tan(e + f*x)), x)]
    @test_int [(c + d*x)^m/(a + im*a*tan(e + f*x))^1, x, 2, (c + d*x)^(1 + m)/(2*a*d*(1 + m)) + (im*2^(-2 - m)*(c + d*x)^m*Gamma(1 + m, ((2*im)*f*(c + d*x))/d))/(a*ℯ^((2*im)*(e - (c*f)/d))*f*((im*f*(c + d*x))/d)^m)]
    @test_int [(c + d*x)^m/(a + im*a*tan(e + f*x))^2, x, 4, (c + d*x)^(1 + m)/(4*a^2*d*(1 + m)) + (im*2^(-2 - m)*(c + d*x)^m*Gamma(1 + m, ((2*im)*f*(c + d*x))/d))/(a^2*ℯ^((2*im)*(e - (c*f)/d))*f*((im*f*(c + d*x))/d)^m) + (im*4^(-2 - m)*(c + d*x)^m*Gamma(1 + m, ((4*im)*f*(c + d*x))/d))/(a^2*ℯ^((4*im)*(e - (c*f)/d))*f*((im*f*(c + d*x))/d)^m)]
    @test_int [(c + d*x)^m/(a + im*a*tan(e + f*x))^3, x, 5, (c + d*x)^(1 + m)/(8*a^3*d*(1 + m)) + ((3*im)*2^(-4 - m)*(c + d*x)^m*Gamma(1 + m, ((2*im)*f*(c + d*x))/d))/(a^3*ℯ^((2*im)*(e - (c*f)/d))*f*((im*f*(c + d*x))/d)^m) + ((3*im)*2^(-5 - 2*m)*(c + d*x)^m*Gamma(1 + m, ((4*im)*f*(c + d*x))/d))/(a^3*ℯ^((4*im)*(e - (c*f)/d))*f*((im*f*(c + d*x))/d)^m) + (im*2^(-4 - m)*3^(-1 - m)*(c + d*x)^m*Gamma(1 + m, ((6*im)*f*(c + d*x))/d))/(a^3*ℯ^((6*im)*(e - (c*f)/d))*f*((im*f*(c + d*x))/d)^m)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*tan(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*tan(e+f*x))^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(c + d*x)^3*(a + b*tan(e + f*x)), x, 8, (a*(c + d*x)^4)/(4*d) + ((im/4)*b*(c + d*x)^4)/d - (b*(c + d*x)^3*log(1 + ℯ^((2*im)*(e + f*x))))/f + (((3*im)/2)*b*d*(c + d*x)^2*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (3*b*d^2*(c + d*x)*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^3) - (((3*im)/4)*b*d^3*PolyLog(4, -ℯ^((2*im)*(e + f*x))))/f^4]
    @test_int [(c + d*x)^2*(a + b*tan(e + f*x)), x, 7, (a*(c + d*x)^3)/(3*d) + ((im/3)*b*(c + d*x)^3)/d - (b*(c + d*x)^2*log(1 + ℯ^((2*im)*(e + f*x))))/f + (im*b*d*(c + d*x)*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (b*d^2*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^3)]
    @test_int [(c + d*x)^1*(a + b*tan(e + f*x)), x, 6, (a*(c + d*x)^2)/(2*d) + ((im/2)*b*(c + d*x)^2)/d - (b*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f + ((im/2)*b*d*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2]
    @test_int [(a + b*tan(e + f*x))/(c + d*x)^1, x, 0, Unintegrable((a + b*tan(e + f*x))/(c + d*x), x)]
    @test_int [(a + b*tan(e + f*x))/(c + d*x)^2, x, 0, Unintegrable((a + b*tan(e + f*x))/(c + d*x)^2, x)]

    @test_int [(c + d*x)^3*(a + b*tan(e + f*x))^2, x, 15, ((-im)*b^2*(c + d*x)^3)/f + (a^2*(c + d*x)^4)/(4*d) + ((im/2)*a*b*(c + d*x)^4)/d - (b^2*(c + d*x)^4)/(4*d) + (3*b^2*d*(c + d*x)^2*log(1 + ℯ^((2*im)*(e + f*x))))/f^2 - (2*a*b*(c + d*x)^3*log(1 + ℯ^((2*im)*(e + f*x))))/f - ((3*im)*b^2*d^2*(c + d*x)*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^3 + ((3*im)*a*b*d*(c + d*x)^2*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 + (3*b^2*d^3*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^4) - (3*a*b*d^2*(c + d*x)*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/f^3 - (((3*im)/2)*a*b*d^3*PolyLog(4, -ℯ^((2*im)*(e + f*x))))/f^4 + (b^2*(c + d*x)^3*tan(e + f*x))/f]
    @test_int [(c + d*x)^2*(a + b*tan(e + f*x))^2, x, 13, ((-im)*b^2*(c + d*x)^2)/f + (a^2*(c + d*x)^3)/(3*d) + (((2*im)/3)*a*b*(c + d*x)^3)/d - (b^2*(c + d*x)^3)/(3*d) + (2*b^2*d*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f^2 - (2*a*b*(c + d*x)^2*log(1 + ℯ^((2*im)*(e + f*x))))/f - (im*b^2*d^2*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^3 + ((2*im)*a*b*d*(c + d*x)*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (a*b*d^2*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/f^3 + (b^2*(c + d*x)^2*tan(e + f*x))/f]
    @test_int [(c + d*x)^1*(a + b*tan(e + f*x))^2, x, 9, -(b^2*c*x) - (b^2*d*x^2)/2 + (a^2*(c + d*x)^2)/(2*d) + (im*a*b*(c + d*x)^2)/d - (2*a*b*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f + (b^2*d*log(cos(e + f*x)))/f^2 + (im*a*b*d*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 + (b^2*(c + d*x)*tan(e + f*x))/f]
    @test_int [(a + b*tan(e + f*x))^2/(c + d*x)^1, x, 0, Unintegrable((a + b*tan(e + f*x))^2/(c + d*x), x)]
    @test_int [(a + b*tan(e + f*x))^2/(c + d*x)^2, x, 0, Unintegrable((a + b*tan(e + f*x))^2/(c + d*x)^2, x)]

    @test_int [(c + d*x)^3*(a + b*tan(e + f*x))^3, x, 28, (((3*im)/2)*b^3*d*(c + d*x)^2)/f^2 - ((3*im)*a*b^2*(c + d*x)^3)/f + (b^3*(c + d*x)^3)/(2*f) + (a^3*(c + d*x)^4)/(4*d) + (((3*im)/4)*a^2*b*(c + d*x)^4)/d - (3*a*b^2*(c + d*x)^4)/(4*d) - ((im/4)*b^3*(c + d*x)^4)/d - (3*b^3*d^2*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f^3 + (9*a*b^2*d*(c + d*x)^2*log(1 + ℯ^((2*im)*(e + f*x))))/f^2 - (3*a^2*b*(c + d*x)^3*log(1 + ℯ^((2*im)*(e + f*x))))/f + (b^3*(c + d*x)^3*log(1 + ℯ^((2*im)*(e + f*x))))/f + (((3*im)/2)*b^3*d^3*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^4 - ((9*im)*a*b^2*d^2*(c + d*x)*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^3 + (((9*im)/2)*a^2*b*d*(c + d*x)^2*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (((3*im)/2)*b^3*d*(c + d*x)^2*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 + (9*a*b^2*d^3*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^4) - (9*a^2*b*d^2*(c + d*x)*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^3) + (3*b^3*d^2*(c + d*x)*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^3) - (((9*im)/4)*a^2*b*d^3*PolyLog(4, -ℯ^((2*im)*(e + f*x))))/f^4 + (((3*im)/4)*b^3*d^3*PolyLog(4, -ℯ^((2*im)*(e + f*x))))/f^4 - (3*b^3*d*(c + d*x)^2*tan(e + f*x))/(2*f^2) + (3*a*b^2*(c + d*x)^3*tan(e + f*x))/f + (b^3*(c + d*x)^3*tan(e + f*x)^2)/(2*f)]
    @test_int [(c + d*x)^2*(a + b*tan(e + f*x))^3, x, 22, (b^3*c*d*x)/f + (b^3*d^2*x^2)/(2*f) - ((3*im)*a*b^2*(c + d*x)^2)/f + (a^3*(c + d*x)^3)/(3*d) + (im*a^2*b*(c + d*x)^3)/d - (a*b^2*(c + d*x)^3)/d - ((im/3)*b^3*(c + d*x)^3)/d + (6*a*b^2*d*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f^2 - (3*a^2*b*(c + d*x)^2*log(1 + ℯ^((2*im)*(e + f*x))))/f + (b^3*(c + d*x)^2*log(1 + ℯ^((2*im)*(e + f*x))))/f - (b^3*d^2*log(cos(e + f*x)))/f^3 - ((3*im)*a*b^2*d^2*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^3 + ((3*im)*a^2*b*d*(c + d*x)*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (im*b^3*d*(c + d*x)*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (3*a^2*b*d^2*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^3) + (b^3*d^2*PolyLog(3, -ℯ^((2*im)*(e + f*x))))/(2*f^3) - (b^3*d*(c + d*x)*tan(e + f*x))/f^2 + (3*a*b^2*(c + d*x)^2*tan(e + f*x))/f + (b^3*(c + d*x)^2*tan(e + f*x)^2)/(2*f)]
    @test_int [(c + d*x)^1*(a + b*tan(e + f*x))^3, x, 16, -3*a*b^2*c*x + (b^3*d*x)/(2*f) - (3*a*b^2*d*x^2)/2 + (a^3*(c + d*x)^2)/(2*d) + (((3*im)/2)*a^2*b*(c + d*x)^2)/d - ((im/2)*b^3*(c + d*x)^2)/d - (3*a^2*b*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f + (b^3*(c + d*x)*log(1 + ℯ^((2*im)*(e + f*x))))/f + (3*a*b^2*d*log(cos(e + f*x)))/f^2 + (((3*im)/2)*a^2*b*d*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - ((im/2)*b^3*d*PolyLog(2, -ℯ^((2*im)*(e + f*x))))/f^2 - (b^3*d*tan(e + f*x))/(2*f^2) + (3*a*b^2*(c + d*x)*tan(e + f*x))/f + (b^3*(c + d*x)*tan(e + f*x)^2)/(2*f)]
    @test_int [(a + b*tan(e + f*x))^3/(c + d*x)^1, x, 0, Unintegrable((a + b*tan(e + f*x))^3/(c + d*x), x)]
    @test_int [(a + b*tan(e + f*x))^3/(c + d*x)^2, x, 0, Unintegrable((a + b*tan(e + f*x))^3/(c + d*x)^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [(c + d*x)^3/(a + b*tan(e + f*x)), x, 6, (c + d*x)^4/(4*(a + im*b)*d) + (b*(c + d*x)^3*log(1 + ((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2))/((a^2 + b^2)*f) - (3*im*b*d*(c + d*x)^2*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/(2*(a^2 + b^2)*f^2) + (3*b*d^2*(c + d*x)*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/(2*(a^2 + b^2)*f^3) + (3*im*b*d^3*PolyLog(4, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/(4*(a^2 + b^2)*f^4)]
    @test_int [(c + d*x)^2/(a + b*tan(e + f*x)), x, 5, (c + d*x)^3/(3*(a + im*b)*d) + (b*(c + d*x)^2*log(1 + ((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2))/((a^2 + b^2)*f) - (im*b*d*(c + d*x)*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/((a^2 + b^2)*f^2) + (b*d^2*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/(2*(a^2 + b^2)*f^3)]
    @test_int [(c + d*x)^1/(a + b*tan(e + f*x)), x, 4, (c + d*x)^2/(2*(a + im*b)*d) + (b*(c + d*x)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2))/((a^2 + b^2)*f) - (im*b*d*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/(2*(a^2 + b^2)*f^2)]
    @test_int [1/((c + d*x)^1*(a + b*tan(e + f*x))), x, 0, Unintegrable(1/((c + d*x)*(a + b*tan(e + f*x))), x)]
    @test_int [1/((c + d*x)^2*(a + b*tan(e + f*x))), x, 0, Unintegrable(1/((c + d*x)^2*(a + b*tan(e + f*x))), x)]

    @test_int [(c + d*x)^3/(a + b*tan(e + f*x))^2, x, 21, -((2*im*b^2*(c + d*x)^3)/((a^2 + b^2)^2*f)) + (2*b^2*(c + d*x)^3)/((a + im*b)*(im*a + b)^2*(im*a - b + (im*a + b)*ℯ^(2*im*e + 2*im*f*x))*f) + (c + d*x)^4/(4*(a - im*b)^2*d) + (b*(c + d*x)^4)/((im*a - b)*(a - im*b)^2*d) - (b^2*(c + d*x)^4)/((a^2 + b^2)^2*d) + (3*b^2*d*(c + d*x)^2*log(1 + ((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b)))/((a^2 + b^2)^2*f^2) + (2*b*(c + d*x)^3*log(1 + ((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b)))/((a - im*b)^2*(a + im*b)*f) - (2*im*b^2*(c + d*x)^3*log(1 + ((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b)))/((a^2 + b^2)^2*f) - (3*im*b^2*d^2*(c + d*x)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a^2 + b^2)^2*f^3) + (3*b*d*(c + d*x)^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((im*a - b)*(a - im*b)^2*f^2) - (3*b^2*d*(c + d*x)^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a^2 + b^2)^2*f^2) + (3*b^2*d^3*PolyLog(3, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/(2*(a^2 + b^2)^2*f^4) + (3*b*d^2*(c + d*x)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a - im*b)^2*(a + im*b)*f^3) - (3*im*b^2*d^2*(c + d*x)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a^2 + b^2)^2*f^3) - (3*b*d^3*PolyLog(4, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/(2*(im*a - b)*(a - im*b)^2*f^4) + (3*b^2*d^3*PolyLog(4, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/(2*(a^2 + b^2)^2*f^4)]
    @test_int [(c + d*x)^2/(a + b*tan(e + f*x))^2, x, 18, -((2*im*b^2*(c + d*x)^2)/((a^2 + b^2)^2*f)) + (2*b^2*(c + d*x)^2)/((a + im*b)*(im*a + b)^2*(im*a - b + (im*a + b)*ℯ^(2*im*e + 2*im*f*x))*f) + (c + d*x)^3/(3*(a - im*b)^2*d) + (4*b*(c + d*x)^3)/(3*(im*a - b)*(a - im*b)^2*d) - (4*b^2*(c + d*x)^3)/(3*(a^2 + b^2)^2*d) + (2*b^2*d*(c + d*x)*log(1 + ((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b)))/((a^2 + b^2)^2*f^2) + (2*b*(c + d*x)^2*log(1 + ((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b)))/((a - im*b)^2*(a + im*b)*f) - (2*im*b^2*(c + d*x)^2*log(1 + ((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b)))/((a^2 + b^2)^2*f) - (im*b^2*d^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a^2 + b^2)^2*f^3) + (2*b*d*(c + d*x)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((im*a - b)*(a - im*b)^2*f^2) - (2*b^2*d*(c + d*x)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a^2 + b^2)^2*f^2) + (b*d^2*PolyLog(3, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a - im*b)^2*(a + im*b)*f^3) - (im*b^2*d^2*PolyLog(3, -(((a - im*b)*ℯ^(2*im*e + 2*im*f*x))/(a + im*b))))/((a^2 + b^2)^2*f^3)]
    @test_int [(c + d*x)^1/(a + b*tan(e + f*x))^2, x, 5, -((c + d*x)^2/(2*(a^2 + b^2)*d)) + (b*d + 2*a*c*f + 2*a*d*f*x)^2/(4*a*(a + im*b)*(a^2 + b^2)*d*f^2) + (b*(b*d + 2*a*c*f + 2*a*d*f*x)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2))/((a^2 + b^2)^2*f^2) - (im*a*b*d*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(e + f*x)))/(a + im*b)^2)))/((a^2 + b^2)^2*f^2) - (b*(c + d*x))/((a^2 + b^2)*f*(a + b*tan(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + b*tan(e + f*x))^2), x, 0, Unintegrable(1/((c + d*x)*(a + b*tan(e + f*x))^2), x)]
    @test_int [1/((c + d*x)^2*(a + b*tan(e + f*x))^2), x, 0, Unintegrable(1/((c + d*x)^2*(a + b*tan(e + f*x))^2), x)]
    end

max_stars_repo_path: test/4 Trig functions/4.3 Tangent/4.3.10 (c+d x)^m (a+b tan)^n.jl


text:
using Plots, GLPlot; GLPlot.init()
using GLPlot, Reactive, GeometryTypes, Colors
using CUDAnative, GPUArrays
import GPUArrays: GPUArray, GLBackend, CUBackend, cu_map
glctx = GLBackend.init(GLVisualize.current_screen())
cuctx = CUBackend.init()

const cu = CUDAnative

function mandelbulb{T}(x0::T,y0::T,z0::T, n, iter)
    x,y,z = x0,y0,z0
    for i=1:iter
        r = cu.sqrt(x*x + y*y + z*z)
        theta = cu.atan2(cu.sqrt(x*x + y*y) , z)
        phi = cu.atan2(y,x)
        rn = cu.pow(r, n)
        x1 = rn * cu.sin(theta*n) * cu.cos(phi*n) + x0
        y1 = rn * cu.sin(theta*n) * cu.sin(phi*n) + y0
        z1 = rn * cu.cos(theta*n) + z0
        (x1*x1 + y1*y1 + z1*z1) > n && return T(i)
        x,y,z = x1,y1,z1
    end
    T(iter)
end
@target ptx function mandelmap(A, x,y,z, n, iter)
    i = Int((blockIdx().x-1) * blockDim().x + threadIdx().x)
    @inbounds if i <= length(A)
        sz = size(A)
        xi,yi,zi = ind2sub(sz, Int(i))
        A[xi,yi,zi] = mandelbulb(x[xi], y[yi], z[zi], n, iter)
    end
    nothing
end
dims = (200,200,200)
vol_gl = GPUArray(Float32, dims, context=glctx);
xrange, yrange, zrange = ntuple(3) do i
    r = linspace(-1f0, 1f0, dims[i]) # linearly spaced array (not dense) from -1 to 1
end
# create two sliders to interact with the 2 parameters of the mandelbulb function
itslider = GLPlot.play_widget(1:15);
nslider = GLPlot.play_widget(linspace(1f0,30f0, 100));
using GLAbstraction
tex = GPUArray(Texture(Float32, dims), dims, context=glctx);
# register a callback to the sliders with map ("map over updates")
volume = preserve(map(nslider, itslider) do n, it
    # map to cuda memory space
    cu_map(vol_gl) do vol_cu
        CUBackend.call_cuda(mandelmap, vol_cu, xrange, yrange, zrange, n, it)
    end
    unsafe_copy!(tex, vol_gl)
end);

glplot(buffer(tex), color_norm = Vec2f0(0, 50))

function iso_particle(v, isoval)
    particles = Point3f0[]
    sz = 1f0./(Point3f0(size(v))-1f0)
    @inbounds for z=1:size(v, 3), y=1:size(v, 2), x=1:size(v, 1)
        if abs(v[x,y,z] - isoval) <= 0.001f0
            push!(particles, (Point3f0(x-1f0,y-1f0,z-1f0)).*sz)
        end
    end
    particles
end
using GLAbstraction, Meshing
isoval = GLPlot.play_widget(1:50);
particles = async_map2(iso_particle, Point3f0[], volume[2], isoval);
pp = GLBuffer(particles[2])
glplot((Sphere(Point2f0(0), 0.005f0),pp), boundingbox=nothing, color=pp)
mesh = GLNormalMesh(volume[2].value, 4f0)
map!(mesh.vertices) do v
    (v-1f0) ./ 298f0
end
glplot(mesh)

using GLAbstraction
compute_context =