text:
using SparseArrays

function clique_expansion_homophily(N,k_range)
    # Compute homophily index in the clique expansion of the hypergraph
    In1 = 0
    In2 = 0
    Cross = 0
    for k = k_range
        for j = 1:k+1
            i = j-1         # number of women in this type hyperedge
            Hik = N[k,j]    # number of hyperedges of this type
            In1 += binomial(i,2)*Hik    # number of extra edges of this type
            In2 += binomial(k-i,2)*Hik
            Cross += Hik*i*(k-i)
        end
    end

    hw = 2*In1/(2*In1 + Cross)
    hm = 2*In2/(2*In2 + Cross)
    return hw, hm
end

function n_baselines(k,n,alpha)
    n1 = round(Int64,floor(n*alpha))
    n2 = n - n1

    b = zeros(k)
    for t = 1:k
        b[t] = binomial(n1-1,t-1)*binomial(n2,k-t)/binomial(n-1,k-1)
    end
    return b

end

function n1_n2_baselines(k,n1,n2)
    n = n1+n2

    b1 = zeros(k)
    b2 = zeros(k)
    for t = 1:k
        b1[t] = binomial(n1-1,t-1)*binomial(n2,k-t)/binomial(n-1,k-1)
        b2[t] = binomial(n2-1,t-1)*binomial(n1,k-t)/binomial(n-1,k-1)
    end
    return b1, b2
end

function n_alt_baselines(k,n,alpha)
    n1 = round(Int64,floor(n*alpha))
    n2 = n - n1

    b = zeros(k)
    denom = 0
    for t = 1:k
        b[t] = binomial(n1,t)*binomial(n2,k-t)
    end
    return b/sum(b)

end

function arbitrary_baselines(k,alpha)
    """
    Returns asymptotic relative baseline scores for a single class that makes up
        a ratio alpha of the total population. Assumes k-uniform hypergraph.

        k = hyperedge size
        alpha = proportion of nodes in the class of interest
    """
    b = zeros(k)
    for j = 1:k
        b[j] = alpha^(j-1)*(1-alpha)^(k-j)*binomial(k-1,j-1)
    end
    return b
end

function alternative_baselines(k,alpha)
    """
    This is for the alternative definition of affinity score, that doesn't
        degree correct.

    Returns asymptotic relative baseline scores for a single class that makes up
        a ratio alpha of the total population. Assumes k-uniform hypergraph.

        k = hyperedge size
        alpha = proportion of nodes in the class of interest
    """
    b = zeros(k)
    denom = sum(alpha^(i)*(1-alpha)^(k-i)/(factorial(i)*factorial(k-i)) for i = 1:k)
    for j = 1:k
        b[j] = alpha^(j)*(1-alpha)^(k-j)/(denom*(factorial(j)*factorial(k-j)))
    end
    return b
end

function baselines(r,alpha)
    """
    Returns asymptotic relative baseline scores for a single class that makes up
        a ratio alpha of the total population. Does this for all hyperedge sizes
        from 1 to r.
    """
    B = zeros(r,r)
    for k = 1:r
        for t = 1:k
            B[k,t] = alpha^(t-1)*(1-alpha)^(k-t)*binomial(k-1,t-1)
        end
    end
    return B
end

function order_degree(H::SparseMatrixCSC{Float64,Int64})
    """Returns hyperedge sizes and degrees for a hypergraph."""
    order = vec(sum(H,dims=2))
    deg = vec(sum(H,dims=1))
    return order, degree
end

function get_hyperedge_counts(H::SparseMatrixCSC{Float64,Int64},classes::Vector{Int64})
    """
    Compute number of each type of hyperedge. Output is a matrix N, such that
        N[k,j] = number of hyperedge of size k that have j-1 nodes of class 1
    """
    order, degree = order_degree(H)
    r = round(Int64,maximum(order))

    # Ensure classes are 0 and 1
    cs = unique(classes)
    @assert(length(cs) == 2)
    @assert(minimum(cs) == 0 && maximum(cs) == 1)
    Hyp = incidence2elist(H)
    N = zeros(r,r+1)
    for edge in Hyp
        board = classes[edge]
        w = sum(board)
        k = length(board)
        if k > 0
            N[k,w+1] += 1
        end
    end
    return N
end

function relative_affinities(H::SparseMatrixCSC{Float64,Int64},classes::Vector{Int64},r::Int64=0)
    """
    Compute relative affinity scores for both classes in a hypergraph H with
        two classes.

        Aff1[k,t] = type-t relative affinity score for size k hyperedges, for class 1
        Aff0[k,t] = type-t relative affinity score for size k hyperedges, for class 0
    """
    N = get_hyperedge_counts(H,classes)
    if r == 0
        r = size(N,1)   # max hyperedge size
    end
    Aff1 = zeros(r,r)
    Aff0 = zeros(r,r)
    for k = 1:r
        c1_denom = sum(i*N[k,i+1] for i = 1:k)
        c0_denom = sum(i*N[k,k-i+1] for i = 1:k)
        for t = 1:k
            Aff1[k,t] = t*N[k,t+1]/c1_denom
            Aff0[k,t] = t*N[k,k-t+1]/c0_denom
        end
    end
    return Aff1, Aff0
end

function affinities_from_N(N)
    r = size(N,1)
    Aff1 = zeros(r,r)
    Aff0 = zeros(r,r)
    for k = 1:r
        c1_denom = sum(i*N[k,i+1] for i = 1:k)
        c0_denom = sum(i*N[k,k-i+1] for i = 1:k)
        for t = 1:k
            Aff1[k,t] = t*N[k,t+1]/c1_denom
            Aff0[k,t] = t*N[k,k-t+1]/c0_denom
        end
    end
    return Aff1, Aff0
end

function relative_affinities_k(H::SparseMatrixCSC{Float64,Int64},classes::Vector{Int64},k::Int64)
    """
    Compute relative affinity scores for both classes in a hypergraph H with
        two classes, only for hyperedges of size k.

        Aff1[t] = type-t relative affinity score for class 1
        Aff0[t] = type-t relative affinity score for class 0
    """
    N_all = get_hyperedge_counts(H,classes)
    N = N_all[k,:]
    Aff1 = zeros(k)
    Aff0 = zeros(k)
    c1_denom = sum(i*N[i+1] for i = 1:k)
    c0_denom = sum(i*N[k-i+1] for i = 1:k)
    for t = 1:k
        Aff1[t] = t*N[t+1]/c1_denom
        Aff0[t] = t*N[k-t+1]/c0_denom
    end
    return Aff1, Aff0
end

function incidence2elist(Hin::SparseMatrixCSC{Float64,Int64})
    """ Convert a hyperedge list to a hypergraph binary incidence matrix. """
    H = sparse(Hin')
    rp = H.rowval
    ci = H.colptr
    Hyperedges = Vector{Vector{Int64}}()
    n,m = size(H)
    for i = 1:m
        startedge = ci[i]
        endedge = ci[i+1]-1
        edge = rp[startedge:endedge]
        push!(Hyperedges,edge)
    end
    return Hyperedges
end

function elist2incidence(Hyperedges::Vector{Vector{Int64}}, N::Int64)
    """
    Take a list of hyperedges and turn it into a hyperedge incidence matrix
    H. N is the number of nodes in the hypergraph.
    H(e,u) = 1  iff node u is in hyperedge e.
    """
    U = Vector{Int64}()
    E = Vector{Int64}()
    M = length(Hyperedges)
    for enum = 1:length(Hyperedges)
        e = Hyperedges[enum]
        for node in e
            push!(U,node)
            push!(E,enum)
        end
    end

    H = sparse(E,U,ones(length(U)),M,N)
    return H
end

function read_hypergraph_data(dataname::String, maxsize::Int64=25)
    """
    Read data from .txt file to incidence matrix, and extract
    class labels.
    """
    classes = Int64[]
    open("../original-data/$dataname/node-labels-$dataname.txt") do f
        for line in eachline(f)
            push!(classes, parse(Int64, line))
        end
    end
    n = length(classes)

    # hyperedges
    EdgeList = Vector{Vector{Int64}}()
    open("../original-data/$dataname/hyperedges-$dataname.txt") do f
        for line in eachline(f)
            edge = [parse(Int64, v) for v in split(line, ',')]
            sort!(edge)
            push!(EdgeList,edge)
        end
    end

    H = elist2incidence(EdgeList,n)

    return H,classes
end

function Hypergraph_to_Scores(H,classes,r)
    n = size(H,2)
    alpha = sum(classes)/n
    B1 = baselines(r,alpha)
    B2 = baselines(r,1-alpha)
    H1, H2 = relative_affinities(H,classes,r)

    N = get_hyperedge_counts(H,classes)
    N = N[1:r,1:r+1]
    R1 = H1./B1
    R2 = H2./B2

    return alpha, B1, B2, H1, H2, R1, R2, N
end

function GHI(H,classes,r)
    """
    Group homophily index: largest value of j such that the top j affinity
    scores are above baseline.
    """
    alpha, B1, B2, H1, H2, R1, R2, N = Hypergraph_to_Scores(H,classes,r)
    Sets1 = zeros(r)
    Sets2 = zeros(r)
    for k = 1:r
        last1 = k
        last2 = k
        c1searching = true
        c2searching = true

        if R1[k,k] < 1
            c1searching = false
            last1 = k+1
        end

        if R2[k,k] < 1
            c2searching = false
            last2 = k+1
        end

        for j = k:-1:1
            # @assert(R1[k,k] >= 1)
            # @assert(R2[k,k] >= 1)
            if R1[k,j] > 1 && c1searching
                last1 = j
            else
                c1searching = false
            end
            if R2[k,j] > 1 && c2searching
                last2 = j
            else
                c2searching = false
            end
            if R1[k,j] < 1 && R2[k,j] < 1
                break
            end
        end
        Sets1[k] = last1
        Sets2[k] = last2
    end

    ## Set-in plot
    S1 = collect(1:r) - Sets1 .+ 1
    S2 = collect(1:r) - Sets2 .+ 1

    return S1, S2
end

function Bootstrap_Affinities(H,k,B,classes)
    n = size(H,2)
    R1 = zeros(B,k)
    R0 = zeros(B,k)

    A1 = zeros(B,k)
    A0 = zeros(B,k)

    order = vec(sum(H,dims = 2))
    Ek = findall(x->x==k,order)

    alpha = sum(classes)/n
    B1 = arbitrary_baselines(k,alpha)
    B0 = arbitrary_baselines(k,1-alpha)
    for i = 1:B
        sam = sample(Ek,length(Ek))
        Hnew = H[sam,:]
        Elist = incidence2elist(Hnew)
        Nk = zeros(k+1)
        for edge in Elist
            elabels = classes[edge]
            t = sum(elabels)
            Nk[t+1] += 1
        end
        c1_denom = sum(i*Nk[i+1] for i = 1:k)
        c0_denom = sum(i*Nk[k-i+1] for i = 1:k)
        for t = 1:k
            A1[i,t] = (t*Nk[t+1]/c1_denom)
            A0[i,t] = (t*Nk[k-t+1]/c0_denom)
            R1[i,t] = (t*Nk[t+1]/c1_denom)/B1[t]
            R0[i,t] = (t*Nk[k-t+1]/c0_denom)/B0[t]
        end
    end

    # Average of the random trials
    MR1 = vec(mean(R1,dims = 1))
    MR0 = vec(mean(R0,dims = 1))
    MA1 = vec(mean(A1,dims = 1))
    MA0 = vec(mean(A0,dims = 1))

    # Compute standard error
    SR1 = zeros(k)
    SR0 = zeros(k)
    SA1 = zeros(k)
    SA0 = zeros(k)
    for t = 1:k
        SR1[t] = StatsBase.std(R1[:,t])
        SR0[t] = StatsBase.std(R0[:,t])
        SA1[t] = StatsBase.std(A1[:,t])
        SA0[t] = StatsBase.std(A0[:,t])
    end

    return R0, R1, MR1, MR0, SR1, SR0, A0, A1, MA1, MA0, SA1, SA0

end

function Bootstrap_From_N(Nk_old,B,alpha)

    k = length(Nk_old)-1

    R1 = zeros(B,k)
    R0 = zeros(B,k)
    A1 = zeros(B,k)
    A0 = zeros(B,k)
    B1 = arbitrary_baselines(k,alpha)
    B0 = arbitrary_baselines(k,1-alpha)

    for i = 1:B
        Nk = zeros(k+1)
        # @show Nk_old
        for j = 1:round(Int64,sum(Nk_old))
            Nk[sample(1:k+1,Nk_old)] += 1
        end
        c1_denom = sum(i*Nk[i+1] for i = 1:k)
        c0_denom = sum(i*Nk[k-i+1] for i = 1:k)
        for t = 1:k
            A1[i,t] = (t*Nk[t+1]/c1_denom)
            A0[i,t] = (t*Nk[k-t+1]/c0_denom)
            R1[i,t] = (t*Nk[t+1]/c1_denom)/B1[t]
            R0[i,t] = (t*Nk[k-t+1]/c0_denom)/B0[t]
        end
    end

    # Average of the random trials
    MR1 = vec(mean(R1,dims = 1))
    MR0 = vec(mean(R0,dims = 1))
    MA1 = vec(mean(A1,dims = 1))
    MA0 = vec(mean(A0,dims = 1))

    # Compute standard error
    SR1 = zeros(k)
    SR0 = zeros(k)
    SA1 = zeros(k)
    SA0 = zeros(k)
    for t = 1:k
        SR1[t] = StatsBase.std(R1[:,t])
        SR0[t] = StatsBase.std(R0[:,t])
        SA1[t] = StatsBase.std(A1[:,t])
        SA0[t] = StatsBase.std(A0[:,t])
    end

    return R0, R1, MR1, MR0, SR1, SR0, A0, A1, MA1, MA0, SA1, SA0

end


text:
string(d::Date, ::Type{IGCDateFormat}) = Dates.format(d, IGC_DATE_FMT)
string(t::Time, ::Type{IGCTimeFormat}) = Dates.format(t, IGC_TIME_FMT)

string(dt::ZonedDateTime, ::Type{IGCTimeFormat}) = string(Time(DateTime(dt)), IGCTimeFormat)
string(dt::ZonedDateTime, ::Type{IGCDateFormat}) = string(Date(DateTime(dt)), IGCDateFormat)

function string(lat::IGCLatitude)
    ordinal = lat.value < 0 ? 'S' : 'N'
    (d, m) = divrem(abs(lat.value), 1)
    d = lpad(round(Int, d), 2, "0")
    m = lpad(round(Int, m * 60 * 1000), 5, "0")
    return "$d$m$ordinal"
end

function write(stream::IO, lat::IGCLatitude)
    write(stream, string(lat))
end

function string(long::IGCLongitude)
    ordinal = long.value < 0 ? 'W' : 'E'
    (d, m) = divrem(abs(long.value), 1)
    d = lpad(round(Int, d), 3, "0")
    m = lpad(round(Int, m * 60 * 1000), 5, "0")
    return "$d$m$ordinal"
end

function write(stream::IO, long::IGCLongitude)
    return write(stream, string(long))
end

function string(alt::T) where {T <: Union{IGCPressureAltitude, IGCGpsAltitude}}
    digits = 5
    if (alt.value < 10^digits) && (alt.value >=0)
        return lpad(alt.value, digits, "0")
    else
        throw(IGCWriteException("Can't output $alt to $digits digits string"))
    end
end

function write(stream::IO, alt::T) where {T <: Union{IGCPressureAltitude, IGCGpsAltitude}}
    return write(stream, string(alt))
end

function string(rec::A_record)
    s = "A" * rec.manufacturer * rec.id
    if length(rec.id_addition) > 0
        s *= " " * rec.id_addition
    end
    return s
end

function string(rec::B_record)
    s = "B"
    s = s * string(rec.time, IGCTimeFormat)
    s = s * string(rec.latitude)
    s = s * string(rec.longitude)
    s = s * FixValidity.string(rec.validity)
    s = s * string(rec.pressure_alt)
    s = s * string(rec.gps_alt)
    s = s * string(rec.extensions_string)
    return s
end

function string(rec::C_record_task_info)
    s = "C"
    s = s * string(rec.declaration.date, IGCDateFormat)
    s = s * string(rec.declaration.time, IGCTimeFormat)
    s = s * string(rec.flight_date, IGCDateFormat)

    s = s * string(rec.number)

    digits = 2
    if (rec.num_turnpoints < 10 ^ digits) && (rec.num_turnpoints >= 0)
        s = s * lpad(rec.num_turnpoints, digits, "0")
    else
        throw(IGCWriteException("Can't write C record (task info) num_turnpoints $(rec.num_turnpoints)"))
    end

    if length(rec.description) > 0
        s = s * " " * rec.description
    end

    return s
end

function string(rec::C_record_waypoint_info)
    s = "C"

    s = s * string(rec.latitude)
    s = s * string(rec.longitude)

    if length(rec.description) > 0
        s = s * " " * rec.description
    end

    return s
end

function string(rec::D_record)
    s = "D"

    s = s * GpsQualifier.string(rec.qualifier)
    s = s * rec.station_id

    return s
end

function string(rec::E_record)
    s = "E"

    s = s * string(rec.time, IGCTimeFormat)
    s = s * rec.tlc
    s = s * rec.extension_string

    return s
end

function string(rec::F_record)
    s = "F"

    s = s * string(rec.time, IGCTimeFormat)

    for satellite in rec.satellites
        s = s * satellite
    end

    return s
end

function string(rec::G_record)
    s = "G"

    s = s * rec.value

    return s
end

function string(rec::T) where {T <: Union{I_record, J_record}}
    s = flc(rec)

    s = s * lpad(length(rec.extensions), 2, "0")

    for extension in rec.extensions
        s = s * lpad(extension.bytes.start, 2, "0")
        s = s * lpad(extension.bytes.stop, 2, "0")
        s = s * extension.type
    end

    return s
end

function string(rec::K_record)
    s = "K"
    s = s * string(rec.time, IGCTimeFormat)
    s = s * rec.value_string
    return s
end

function string(rec::L_record)
    s = "L"
    s = s * rec.source[1:3]
    s = s * rec.comment
    return s
end

# ToDo... write all H records
function string_fr_header(source::Char, tlc::String, value; subtype_long="")
    s = "H"
    s = s * source
    s = s * tlc
    if length(subtype_long) > 0
        s = s * subtype_long * ":"
    end
    s = s * string(value)
    return s
end

function write_fr_header(stream::IO, source::Char, tlc::String, value; subtype_long="")
    return write(stream, string_fr_header(source, tlc, value, subtype_long=subtype_long))
end

function string(rec::H_record_FiXAccuracy)
    return string_fr_header(rec.source, tlc(rec), lpad(rec.value, 3, "0"))
end

function string(rec::H_record_DaTE)
    return string_fr_header(rec.source, tlc(rec), string(rec.value, IGCDateFormat))
end

function string(rec::T) where {T <: Union{H_record_PiLoT, H_record_Copilot, H_record_GliderType, H_record_GliderRegistration, H_record_GpsDatum, H_record_FirmwareRevision, H_record_HardwareRevision, H_record_CompetitionId, H_record_CompetitionClass, H_record_MeansOfPropulsion, H_record_Site}}
    return string_fr_header(rec.source, tlc(rec), rec.value, subtype_long=subtype_long(rec))
end

function string(rec::H_record_ManufacturerModel)
    return string_fr_header(rec.source, tlc(rec),
        join([rec.manufacturer, rec.model], ","),
        subtype_long=subtype_long(rec)
    )
end

function string(rec::H_record_GpsReceiver)
    return string_fr_header(rec.source, tlc(rec),
        ":" * join([rec.manufacturer,
            rec.model,
            ismissing(rec.channels) ? "" : "$(rec.channels)ch",
            ismissing(rec.max_alt.value) ? "" : "max$(rec.max_alt.value)$(rec.max_alt.unit)"]
        , ",")
    )
end

function string(rec::H_record_PressureAltitudeSensor)
    return string_fr_header(rec.source, tlc(rec),
        join([rec.manufacturer,
            rec.model,
            ismissing(rec.max_alt.value) ? "" : "max$(rec.max_alt.value)$(rec.max_alt.unit)"]
        , ","),
        subtype_long=subtype_long(rec)
    )
end

function string(rec::H_record_UnitsOfMeasure)
    return string_fr_header(rec.source, tlc(rec),
        join(rec.value, ","),
        subtype_long=subtype_long(rec)
    )
end

function string(rec::H_record_TimeZoneOffset)
    (h, m) = divrem(rec.value.offset.std, 60 * Dates.Second(60))
    m = abs(div(m, Dates.Second(60)))
    if m == 0
        s = "$h"
    else
        s = "$h.$(lpad(m, 2, "0"))"
    end
    return string_fr_header(rec.source, tlc(rec),
        s, subtype_long=subtype_long(rec))
end

function write(stream::IO, task::Task)
    write(stream, task.info)
    for wpt in task.waypoints
        write(stream, wpt)
    end
end

function write(stream::IO, rec::T) where {T <: Abstract_IGC_record}
    return write(stream, string(rec))
end

# ========================================

function write(stream::IO, header::Header; eol=EOL_DEFAULT)
    if !ismissing(header.fix_accuracy)
        write(stream, header.fix_accuracy)
        write(stream, eol)
    end

    write(stream, header.date)
    write(stream, eol)
    write(stream, header.pilot_in_charge)
    write(stream, eol)
    write(stream, header.crew2)
    write(stream, eol)
    write(stream, header.glider_type)
    write(stream, eol)
    write(stream, header.glider_id)
    write(stream, eol)
    write(stream, header.gps_datum)
    write(stream, eol)
    write(stream, header.firmware_version)
    write(stream, eol)
    write(stream, header.hardware_version)
    write(stream, eol)
    write(stream, header.flight_recorder_type)
    write(stream, eol)
    write(stream, header.gps_receiver)
    write(stream, eol)
    write(stream, header.pressure_altitude_sensor)

    if !ismissing(header.competition_id)
        write(stream, eol)
        write(stream, header.competition_id)
    end

    if !ismissing(header.competition_class)
        write(stream, eol)
        write(stream, header.competition_class)
    end

end

# ========================================

function write(stream::IO, records::Vector{T}) where {T <: Union{B_record, D_record, E_record, F_record, L_record, G_record}}
    for rec in records
        write(stream, rec)
    end
end

# ========================================

function write(stream::IO, igcdoc::IGCDocument; eol=EOL_DEFAULT)
    # A record
    write(stream, igcdoc.logger_id)
    write(stream, eol)

    # H records
    # @assert iscomplete(igcdoc.header)
    write(stream, igcdoc.header, eol=eol)

    # I records
    if length(igcdoc.fix_record_extensions.extensions) > 0
        write(stream, eol)
        write(stream, igcdoc.fix_record_extensions)
    end

    # J records
    if length(igcdoc.k_record_extensions.extensions) > 0
        write(stream, eol)
        write(stream, igcdoc.k_record_extensions)
    end

    # C records
    if have_task_declared(igcdoc)
        write(stream, igcdoc.task)
    end

    #= flight =#

    # F records
    if length(igcdoc.satellite_records) > 0
        write(stream, igcdoc.satellite_records)
    end

    # B records
    if length(igcdoc.fix_records) > 0
        write(stream, igcdoc.fix_records)
    end

    # D records
    if length(igcdoc.dgps_records) > 0
        write(stream, igcdoc.dgps_records)
    end

    # E records
    if length(igcdoc.event_records) > 0
        write(stream, igcdoc.event_records)
    end

    #= at end =#

    # L records
    if length(igcdoc.comment_records) > 0
        write(stream, igcdoc.comment_records)
    end

    # G records
    if length(igcdoc.security_records) > 0
        write(stream, igcdoc.security_records)
    end
end

#=
    logger_id::A_record
    fix_records::Vector{B_record}
    tasks::Vector{Task}
    dgps_records::Vector{D_record}
    event_records::Vector{E_record}
    satellite_records::Vector{F_record}
    security_records::Vector{G_record}
    header::Header
    fix_record_extensions::Vector{I_record}
    k_record_extensions::Vector{J_record}
    k_records::Vector{K_record}
    comment_records::Vector{L_record}
    errors::Vector{IGCParseException}
=#

ext: jl


text:
# Load Julia packages (libraries) needed for clip

using StatisticalRethinking

ProjDir = @__DIR__

# Define the Stan language model

m_5_3 = "
data {
  int N;
  vector[N] D;
  vector[N] M;
  vector[N] A;
}
parameters {
  real a;
  real bA;
  real bM;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + bA * A + bM * M;
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  D ~ normal( mu , sigma );
}
";

# Define the SampleModel
m5_3s = SampleModel("m5.3", m_5_3);

# Input data

m5_3_data = Dict(
  "N" => size(df, 1),
  "D" => df[:, :D],
  "M" => df[:, :M],
  "A" => df[:, :A]
);

# Sample using cmdstan

rc = stan_sample(m5_3s, data=m5_3_data);

if success(rc)

  # Describe the draws

  dfa3 = read_samples(m5_3s; output_format=:dataframe)
  Particles(dfa3)

end


text:
function BatchedSampler(net, sampler, prob, algo=Gradient();
                        kwargs...)
    if net isa CachedNet
        throw("Only takes standard network.")
    end

    if prob isa KLocalLiouvillian
        return BatchedGradSampler(net, sampler, prob, algo; kwargs...)
    elseif prob isa AbsLinearOperator
        return BatchedValSampler(net, sampler, prob, algo; kwargs...)
    else
        throw("not known")
    end
end


text:
 #-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*---------------------------------------------------------------------
#* Author: Estefane Horn, N�bia do Carmo Ferreira
#*$Id$
#*-------------------------------------------------------------------
#-------------------------------------------------------------------
#* Model of a valve (simplified)
#*--------------------------------------------------------------------
#*
#* Author: Paula B. Staudt
#*--------------------------------------------------------------------
type valve_flow
	valve_flow()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Minimum Flow Allowable in the valve",
				:Default=>0
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Maximum Flow Allowable in the valve",
				:Default=>1000
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.7365,
				:Symbol=>"_{in}",
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.7365,
				:Symbol=>"_{out}",
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Flow Signal",
				:PosX=>0.5,
				:PosY=>0,
				:Protected=>true
			)),
			[
				:(Outlet.F = Inlet.F),
				:(Outlet.T = Inlet.T),
				:(Outlet.P = Inlet.P),
				:(Outlet.h = Inlet.h),
				:(Outlet.v = Inlet.v),
				:(Outlet.z = Inlet.z),
				:(Outlet.F = MinFlow + FlowFraction*(MaxFlow-MinFlow)),
			],
			[
				"Overall Molar Balance","Temperature","Pressure","Energy Balance","Vapour fraction","Componente Molar Balance","Flow computation",
			],
			[:PP,:NComp,:MinFlow,:MaxFlow,],
			[:Inlet,:Outlet,:FlowFraction,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	MinFlow::flow_mol
	MaxFlow::flow_mol
	Inlet::stream
	Outlet::stream
	FlowFraction::control_signal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export valve_flow
function setEquationFlow(in::valve_flow)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
end
function atributes(in::valve_flow,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Valve"
	fields[:Brief]="Model of a very simple valve for setting the flow with a controller."
	fields[:Info]="== ASSUMPTIONS ==
* nothing happens in this valve

== SET ==
* MinFlow: the Minimum Flow Allowable in the valve;
* MaxFlow: the Maximum Flow Allowable in the valve;

== SPECIFY ==
* the Inlet stream
* the FlowFraction (the model requires an inlet signal, also you can use a controller for setting the FlowFraction)
"
	drive!(fields,_)
	return fields
end
valve_flow(_::Dict{Symbol,Any})=begin
	newModel=valve_flow()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#-----------------------------------------------------------------------------# Node
# Like Hyperscript.Node, but with less magic
struct Node{tag}
    attrs::Vector{Pair{String,String}}
    children
    classes::Vector{String}
    style::String
    Node(; tag, attrs, children, classes, style) = new{tag}(attrs, children, classes, style)
end
function Node{T}(children...; kw...) where {T}
    attrs = [string(k) => string(v) for (k,v) in kw]
    Node(; tag=T, attrs, children, classes=String[], style="")
end

Base.getproperty(o::Node{T}, name) where {T} = Base.getproperty(o, Symbol(name))

function Base.getproperty(o::Node{T}, name::Symbol) where {T}
    f = Fields(o)
    classes = vcat(f.classes, split(string(name)))
    Node(; tag=T, attrs=f.attrs, children=f.children, classes, style=f.style)
end

# To switch getproperty back to getfield
struct Fields{T}
    item::T
end
Base.getproperty(f::Fields, prop::Symbol) = getfield(getfield(f, :item), prop)

#-----------------------------------------------------------------------------# htm
function write_opening_tag(io::IO, node::Node{T}) where {T}
    print(io, '<', T)
    for (k,v) in getfield(node, :attrs)
        v == true ? print(io, ' ', k) : print(io, ' ', k, '=', '"', v, '"')
    end
    print(io, " class=\"",  join(getfield(node, :classes), ' '), '"')
    print(io, '>')
end

function write_htm(io::IO, o::Node{T}, i=1) where {T}
    f = Fields(o)
    children = f.children
    for (j, child) in enumerate(children)
        child isa Node && write_htm(io, child, i+j)
    end
    print(io, "C$i = props => {\n    console.log(`C$i with props: \${JSON.stringify(props)}`)\n   return html`")
    write_opening_tag(io, o)
    for (j, child) in enumerate(children)
        if child isa Node
            write(io, "<\${this.C$(i + j)} />")
        elseif child isa String
            write(io, child)
        else
            show(io, MIME"text/html"(), child)
        end
    end
    println(io, "</$T>`\n}")
end

htm(o::Node, i=0) = (io = IOBuffer(); write_htm(io, o, i); String(take!(io)))

#-----------------------------------------------------------------------------# show
Base.show(io::IO, node::Node) = show(io, MIME"text/html"(), node)

function Base.show(io::IO, ::MIME"text/html", node::Node{T}) where {T}
    compact = get(io, :compact, false)
    color = get(io, :tagcolor, 1)
    printstyled(io, "<$T"; color)
    f = Fields(node)
    attrs, classes, children = f.attrs, f.classes, f.children
    for (k,v) in attrs
        if v == "true"
            printstyled(io, ' ', k; color)
        elseif k == "__VERBATIM__"
            printstyled(io, ' ', v; color)
        else
            printstyled(io, ' ', k, '=', '"', v, '"'; color)
        end
    end
    !isempty(classes) && printstyled(io, " class=", "\"$(join(classes, ' '))\""; color)
    printstyled(io, ">"; color)
    compact || println(io)

    for (i, child) in enumerate(children)
        ctx = IOContext(io, :tagcolor => color + i, :compact => true)
        print(io, "  ")
        if child isa Node
            show(ctx, MIME"text/html"(), child)
        else
            try
                show(ctx, MIME"text/html"(), child)
            catch
                print(ctx, child)
            end
            i == length(children) && print(io)
        end
    end
    printstyled(io, "</$T>\n"; color)
end

#-----------------------------------------------------------------------------# H
struct H end
Base.getproperty(::H, x::Symbol) = Node{x}
h = H()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libnabo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libnabo")
JLLWrappers.@generate_main_file("libnabo", UUID("569caa04-3db7-54b1-b5fe-10f3ec93054b"))
end  # module libnabo_jll


text:
using Revise
includet(joinpath("..", "init.jl"))

db_url_in = joinpath(
    "sqlite:///$(spineitemsPath)", "case_study_b2_rolling", "case_study_b2_rolling.sqlite"
)
db_url_out = joinpath(
    "sqlite:///$(spineitemsPath)", "case_study_b2_rolling_out", "case_study_b2_rolling_out.sqlite"
)

m = SpineOpt.run_spinemodel(db_url_in, db_url_out; cleanup=false, log_level=3)


text:
abstract type Combiner end
abstract type CorPol <: Combiner end

struct PCC <: Combiner end
struct PCC2 <: Combiner end
struct SNR_MAX <: Combiner end
struct EGC <: CorPol end

struct MRC_PSD <: CorPol
    fs::Int
end

MRC_PSD() = MRC_PSD(10)

"""
    get_weights(comb, x, ref)

Returns the weights for a section of the PSM by the method comb
"""
function get_weights(comb::PCC, x::AbstractMatrix, ref::Int)
    return cor(x, view(x, :, ref)) |> vec
end

function get_weights(comb::PCC2, x::AbstractMatrix, ref::Int)
    pows = sum(abs2, x, dims=1) .|> sqrt
    cors = @views sum(x .* x[:, ref], dims=1)
    w = cors ./ (pows .* maximum(pows)) |> vec
    return w
end

function get_weights(comb::SNR_MAX, x::AbstractMatrix, ref::Int)
    cors = @views sum(x .* x[:, ref], dims=1)
    w = cors ./ maximum(cors) |> vec
    return w
end

get_weights(comb::Combiner, x::AbstractMatrix) = get_weights(comb, x, choose_ref(x))

function get_weights(comb::MRC_PSD, x::AbstractMatrix{T}) where T
    w = Vector{T}(undef, size(x, 2))
    for (i, s) in enumerate(eachcol(x))
        w[i] = estimate_snr(s, fs=comb.fs)
    end
    return w
end

"""
    combiner(comb, x, ref)

Gets weights from combining method comb with a chosen reference sensor and applies it to x.
"""
function combiner(comb::Combiner, x::AbstractMatrix, ref)
    w = get_weights(comb, x, ref)
    return x * w
end
"""
    combiner(comb, x)

Gets weights from combining method comb and applies them to x
"""
function combiner(comb::Combiner, x::AbstractMatrix)
    w = get_weights(comb, x)
    return x * w
end

function combiner(comb::CorPol, x::AbstractMatrix)
    w = get_weights(comb, x)
    return polarity_flip(x) * w
end

combiner(comb::EGC, x::AbstractMatrix) = sum(polarity_flip(x), dims=2) |> vec

combiner(x::AbstractMatrix) = combiner(SNR_MAX(), x)

"""
    snrmax(x)

Perform SNR-MAX combining on x. Equivalent to combiner(SNR_MAX(), x)
"""
snrmax(x::AbstractMatrix) = combiner(SNR_MAX(), x)

"""
    mrc(x, fs)

Perform MRC-PSD combining on x. Currently x must be less than 512 samples.
Equivalent to combiner(MRC_PSD(fs), x)
"""
mrc(x::AbstractMatrix, fs) = combiner(MRC_PSD(fs), x)

"""
    egc(x)

Perform equal gain combining on x. Equivalent to combiner(EGC(), x)
"""
egc(x::AbstractMatrix) = combiner(EGC(), x)

"""
    pcc(x)

Perform Pearson Correlation Coefficient combining on x.
Equivalent to combiner(PCC2(), x)
"""
pcc(x::AbstractMatrix) = combiner(PCC2(), x)

"""
    selection(x)

Perform selection combining on x.
Equivalent to extract_ref(x)
"""
selection(x::AbstractMatrix) = extract_ref(x)


text:
using LinearRationalExpectations
using LinearAlgebra
using Random
using Test

include("linear_rational_expectations_test.jl")
include("extended_lyapunov_test.jl")
include("variance_test.jl")


text:
module Sessions

using Random
using JLD2: save_object, load_object
using JSON
using Base.Filesystem: basename, splitext, joinpath, rm, abspath, isfile
using Dates: format, now, Time
using NPZ

using ..Agents, ..Environments

export Session, load_session, save_session!, Trial, DMTrial, load_trial, save_trial, save_trial!

mutable struct Session{O,A,E<:AbstractEnvironment{O,A},AG<:AbstractAgent{O,A}}
    name::String
    agent::AG
    frames::Int
    elapsed_frames::Int
    # replications::Int
    # elapsed_reps::Int
    cummulative_rewards::Vector{Float64}
    # steps_per_episode::Vector{Int}
    checkpoint_interval::Time
    keep_all::Bool
    last_checkpoint::String
end

function Session(name, ::Type{E}, agent::AG, frames, checkpoint_interval, keep_all::Bool = false) where {O,A,E<:AbstractEnvironment{O,A},AG<:AbstractAgent{O,A}}
    frames ≤ 0 && throw(DomainError(frames, "number of frames must be greater then 0"))
    # checkpoint_interval > episodes && throw(
    #   DomainError(checkpoint_interval, "interval may not be greater then the number of episodes"))
    Session{O,A,E,AG}(name, agent, frames, 0, Float64[], checkpoint_interval, keep_all, "")
end

function Base.show(io::IO, s::Session{O,A,E,AG}) where {O,A,E,AG}
    print(
        io,
        """ Session $(s.name)
         ├ Environment: $(E)
         ├ Agent : $(nameof(AG))
         ├ Frames: $(s.frames)
         ├ Elapsed frames: $(s.elapsed_frames)
         ├ Checkpoint interval: $(s.checkpoint_interval)
         └ Keep all: $(s.keep_all)
        """
    )
end

function save_session!(s::Session; path::String = pwd())
    path = abspath(path)
    filename = s.name |> basename |> splitext |> first
    filename = filename * "_" * format(now(), "yyyy-mm-ddTHH-MM-SS") * ".jld2"
    absfilename = joinpath(path, filename)

    save_object(absfilename, s)

    !s.keep_all && isfile(s.last_checkpoint) && rm(s.last_checkpoint)

    s.last_checkpoint = absfilename

    nothing
end

function load_session_from_spec(specname::String)
    spec = JSON.parsefile(specname, dicttype = Dict{Symbol,Any})

    name = specname |> basename |> splitext |> first
    env_type = Environments.environment_table[Symbol(spec[:session][:env][:name])]
    agent = make_agent(env_type, spec[:session][:agent])
    frames = spec[:session][:simulation][:frames]
    checkpoint_interval = Time(spec[:session][:simulation][:checkpoint_interval])
    keep_all = spec[:session][:simulation][:keep_all]

    Session(name, env_type, agent, frames, checkpoint_interval, keep_all)
end

load_session_from_jld(name::String) = load_object(name)

function load_session(name::String)
    ext = name |> basename |> splitext |> last

    if ext == ".jld2"
        return load_session_from_jld(name)
    elseif ext == ".json"
        return load_session_from_spec(name)
    else
        throw(ArgumentError("unable to load session from $name"))
    end
end

# ----------------------------------- Trial ---------------------------------- #
mutable struct Trial{O,A,E<:AbstractEnvironment{O,A},AG<:AbstractAgent{O,A}}
    name::String
    agent::AG
    episodes::Int
    elapsed_episodes::Int
    replications::Int
    elapsed_reps::Int
    cummulative_rewards::Matrix{Float64}
    # steps_per_episode::Vector{Int}
    checkpoint_interval::Time
    keep_all::Bool
    last_checkpoint::String
    rng::MersenneTwister
end

function Trial(name, ::Type{E}, agent::AG, episodes, replications, checkpoint_interval, keep_all::Bool = false; seed::Union{Nothing,UInt} = nothing) where {O,A,E<:AbstractEnvironment{O,A},AG<:AbstractAgent{O,A}}
    episodes ≤ 0 && throw(DomainError(episodes, "number of episodes must be greater then 0"))
    # checkpoint_interval > episodes && throw(
    #   DomainError(checkpoint_interval, "interval may not be greater then the number of episodes"))
    Trial{O,A,E,AG}(name, agent, episodes, 0, replications, 0, Matrix{Float64}(undef, episodes, replications), checkpoint_interval, keep_all, "", MersenneTwister(seed))
end

function Base.show(io::IO, t::Trial{O,A,E,AG}) where {O,A,E,AG}
    print(
        io,
        """ Trial $(t.name)
         ├ Environment: $(E)
         ├ Agent : $(nameof(AG))
         ├ Episodes: $(t.episodes)
         ├ Replications: $(t.replications)
         ├ Elapsed replications: $(t.elapsed_reps)
         ├ Elapsed episodes: $(t.elapsed_episodes)
         ├ Checkpoint interval: $(t.checkpoint_interval)
         └ Keep all: $(t.keep_all)
        """
    )
end

function save_trial!(t::Trial; path::String = pwd())
    path = abspath(path)
    filename = t.name |> basename |> splitext |> first
    filename = filename * "_" * format(now(), "yyyy-mm-ddTHH-MM-SS") * ".jld2"
    absfilename = joinpath(path, filename)

    save_object(absfilename, t)

    !t.keep_all && isfile(t.last_checkpoint) && rm(t.last_checkpoint)

    t.last_checkpoint = absfilename

    nothing
end

# ----------------------------- Deepmind-like Trial ----------------------------
mutable struct DMTrial{O,A,E<:AbstractEnvironment{O,A}}
    name::String
    agent_spec::Dict{Symbol,Any}
    total_timesteps::Int
    replications::Int
    eval_freq::Int
    eval_episodes::Int
    timesteps::Vector{Int}
    results::Dict{UInt,Matrix{Float64}}
    rng::MersenneTwister
end

function DMTrial(name, ::Type{E}, agent_spec, total_timesteps, replications, eval_freq, eval_episodes; seed::Union{Nothing,UInt} = nothing) where {O,A,E<:AbstractEnvironment{O,A}}
    total_timesteps ≤ 0 && throw(DomainError(total_timesteps, "number of timesteps must be greater then 0"))
    # checkpoint_interval > episodes && throw(
    #   DomainError(checkpoint_interval, "interval may not be greater then the number of episodes"))
    DMTrial{O,A,E}(name, agent_spec, total_timesteps, replications, eval_freq, eval_episodes, collect(eval_freq*Base.OneTo(fld(total_timesteps, eval_freq))), Dict{Int,Matrix{Float64}}(), MersenneTwister(seed))
end

function Base.show(io::IO, t::DMTrial{O,A,E}) where {O,A,E}
    print(
        io,
        """ Trial $(t.name)
         ├ Environment: $(E)
         ├ Agent spec : $(t.agent_spec)
         ├ Total timesteps: $(t.total_timesteps)
         └ Replications: $(t.replications)
        """
    )
end

# Save it in the format expected by the functions in Python
# A folder named "results"
# Inside, we have folders for each replication, name after their seeds
# Inside those, we have one "evaluations.npz" file for each folder
# The npz file has two fields: timesteps, a array with the timesteps where
# evaluations took place, results, a matrix where each column contains the
# evalutions along the aforementioned timesteps, one column for each seed.
function save_trial(t::DMTrial; path::String = pwd())
    resultsdir = mkpath(joinpath(path, "results"))

    for (seed, results) in t.results
        seeddir = mkpath(joinpath(resultsdir, "$(seed)"))
        npzwrite(
            joinpath(seeddir, "evaluations.npz"),
            Dict(
                "timesteps" => t.timesteps,
                "results" => results
            )
        )
    end

    nothing
end

# ----------------------------- General functions ------------------------------

function load_trial_from_spec(specname::String)
    spec = JSON.parsefile(specname, dicttype = Dict{Symbol,Any})
    dm_like = haskey(spec, :dmtrial)
    spec = spec |> values |> first

    name = specname |> basename |> splitext |> first
    env_type = Environments.environment_table[Symbol(spec[:env][:name])]
    seed = UInt(spec[:simulation][:seed])
    replications = spec[:simulation][:replications]

    if dm_like
        total_timesteps = spec[:simulation][:total_timesteps]
        eval_freq = spec[:simulation][:eval_freq]
        eval_episodes = spec[:simulation][:eval_episodes]

        return DMTrial(name, env_type, spec[:agent], total_timesteps, replications, eval_freq, eval_episodes; seed)
    else
        agent = make_agent(env_type, spec[:agent]; seed)
        checkpoint_interval = Time(spec[:simulation][:checkpoint_interval])
        keep_all = spec[:simulation][:keep_all]
        episodes = spec[:simulation][:episodes]
        return Trial(name, env_type, agent, episodes, replications, checkpoint_interval, keep_all; seed)
    end
end

load_trial_from_jld(name::String) = load_object(name)

function load_trial(name::String)
    ext = name |> basename |> splitext |> last

    if ext == ".jld2"
        return load_trial_from_jld(name)
    elseif ext == ".json"
        return load_trial_from_spec(name)
    else
        throw(ArgumentError("unable to load trial from $name"))
    end
end

end


text:
# The weight_space_approx function needed for ApproximateGPs.pathwise_sample
# Returns both the sampled basis and the distribution over w needed to approximate the GP
function gp_rff_approx(rng, kernel, input_dims, feature_dims)
    ϕ = sample_basis(rng, kernel, input_dims, feature_dims)
    p_w = MvNormal(Diagonal(Fill(1., feature_dims)))
    return ϕ, p_w
end

# Everything necessary to create a RFF approximation to a stationary kernel
# Currently ony supports SqExponentialKernel with variance and lengthscale
struct RFFBasis{Tinner, Touter, Tω, Tτ, Tsample}
    inner_weights::Tinner  # lengthscale
    outer_weights::Touter  # variance (scaled)
    ω::Tω  # Sampled frequencies;               size(ω): (input_dims, num_features)
    τ::Tτ  # Sampled phases;                    size(τ): (num_features,)
    sample_params::Tsample  # Returns a new sample of ω & τ
end

function (ϕ::RFFBasis)(x)
    # ϕ: R^{input_dims} -> R^{num_features}
    x_ = x ./ ϕ.inner_weights
    ωt_x = ϕ.ω'x_
    return ϕ.outer_weights * cos.(ωt_x .+ ϕ.τ)
end

function resample!(ϕ::RFFBasis)
    ω, τ = ϕ.sample_params()
    ϕ.ω .= ω
    ϕ.τ .= τ
end

# Currently need to pass `input_dims` explicitly - will eventually be in KernelFunctions
# https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/issues/16
function sample_basis(rng, kernel, input_dims, num_features=100)
    inner, outer = spectral_weights(kernel)
    outer_scaled = outer * √(2/num_features)
    p_ω = spectral_distribution(kernel, input_dims)

    function sample_params()
        ω = rand(rng, p_ω, num_features)  #
        τ = rand(rng, Uniform(0, 2π), num_features)  #
        return ω, τ
    end

    return RFFBasis(inner, outer_scaled, sample_params()..., sample_params)
end

function spectral_distribution(::SqExponentialKernel, input_dims)
    return MvNormal(Diagonal(Fill(1., input_dims)))
end

spectral_distribution(k::ScaledKernel, args...) = spectral_distribution(k.kernel, args...)
spectral_distribution(k::TransformedKernel, args...) = spectral_distribution(k.kernel, args...)

spectral_distribution(k::Kernel) =  error("Spectral distribution not implemented for kernel:\n$k")

function spectral_weights(::SqExponentialKernel)
    return 1.0, 1.0
end

function spectral_weights(k::ScaledKernel)
    σ² = only(k.σ²)
    inner, outer = spectral_weights(k.kernel)
    return inner, outer * √σ²
end

function spectral_weights(k::TransformedKernel{<:Any, <:ScaleTransform})
    s = only(k.transform.s)
    inner, outer = spectral_weights(k.kernel)
    return inner / s, outer
end

spectral_weights(k::Kernel) = error("Spectral weights not implemented for kernel:\n$k")

# TODO:
# ARDTransform
# ProductKernel
# SumKernel
# MaternKernel


text:
#
# Potential energy function using r2
#

function upair2(r2,data)
  r6 = r2^3
  r12 = r6^2
  return data.ε4*(data.σ12/r12 - data.σ6/r6)
end


text:
using MyExample
using Test
#2x + 3y
#comments
@testset "MyExample.jl" begin
    @test MyExample.my_f(2,1) == 7
    @test MyExample.my_f(2, 3) == 13
    MyExample.my_f(2, 1)
end


text:
int_rules_4_1_3_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.3.1*(a+b*sin)^m*(c+d*sin)^n*(A+B*sin) =#
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(n') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => Antiderivative(ExpandTrig(sin(~e + ~f * ~x) ^ ~n * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~A + ~B * sin(~e + ~f * ~x)), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x) && (EqQ(~A * ~b + ~a * ~B, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m) && IntegerQ(~n))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (~a) ^ ~m * (~c) ^ ~m * Antiderivative(cos(~e + ~f * ~x) ^ (2 * ~m) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n - ~m) * (~A + ~B * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m) && Not(IntegerQ(~n) && (LtQ(~m, 0) && GtQ(~n, 0) || (LtQ(0, ~n, ~m) || LtQ(~m, ~n, 0)))))))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~A * ~c + (~B * ~c + ~A * ~d) * sin(~e + ~f * ~x) + ~B * ~d * sin(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~A * ~b + ~a * ~B) / (2 * ~a * ~b)) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sqrt(~c + ~d * sin(~e + ~f * ~x)), ~x) + ((~B * ~c + ~A * ~d) / (2 * ~c * ~d)) * Antiderivative(sqrt(~c + ~d * sin(~e + ~f * ~x)) / sqrt(~a + ~b * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~B) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~f * (~m + ~n + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (EqQ(~A * ~b * (~m + ~n + 1) + ~a * ~B * (~m - ~n), 0) && NeQ(~m, -1 / 2))))
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (~B / ~d) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1), ~x) - ((~B * ~c - ~A * ~d) / ~d) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~A * ~b - ~a * ~B) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~a * ~f * (2 * ~m + 1)) + ((~a * ~B * (~m - ~n) + ~A * ~b * (~m + ~n + 1)) / (~a * ~b * (2 * ~m + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && ((LtQ(~m, -1 / 2) || ILtQ(~m + ~n, 0) && Not(SumSimplerQ(~n, 1))) && NeQ(2 * ~m + 1, 0))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~B) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~f * (~m + ~n + 1)) - ((~B * ~c * (~m - ~n) - ~A * ~d * (~m + ~n + 1)) / (~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(LtQ(~m, -1 / 2)) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~B * ~c - ~A * ~d) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (EqQ(~m + ~n + 2, 0) && EqQ(~A * (~a * ~d * ~m + ~b * ~c * (~n + 1)) - ~B * (~a * ~c * ~m + ~b * ~d * (~n + 1)), 0)))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-((~b) ^ 2) * (~B * ~c - ~A * ~d) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * (~b * ~c + ~a * ~d)) - (~b / (~d * (~n + 1) * (~b * ~c + ~a * ~d))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1) * Simp((~a * ~A * ~d * ((~m - ~n) - 2) - ~B * (~a * ~c * (~m - 1) + ~b * ~d * (~n + 1))) - (~A * ~b * ~d * (~m + ~n + 1) - ~B * (~b * ~c * ~m - ~a * ~d * (~n + 1))) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (GtQ(~m, 1 / 2) && (LtQ(~n, -1) && (IntegerQ(2 * ~m) && (IntegerQ(2 * ~n) || EqQ(~c, 0))))))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~b) * ~B * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n + 1)) + (1 / (~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n * Simp(~a * ~A * ~d * (~m + ~n + 1) + ~B * (~a * ~c * (~m - 1) + ~b * ~d * (~n + 1)) + (~A * ~b * ~d * (~m + ~n + 1) - ~B * (~b * ~c * ~m - ~a * ~d * (2 * ~m + ~n))) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (GtQ(~m, 1 / 2) && (Not(LtQ(~n, -1)) && (IntegerQ(2 * ~m) && (IntegerQ(2 * ~n) || EqQ(~c, 0))))))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~A * ~b - ~a * ~B) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~a * ~f * (2 * ~m + 1)) - (1 / (~a * ~b * (2 * ~m + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n - 1) * Simp((~A * (~a * ~d * ~n - ~b * ~c * (~m + 1)) - ~B * (~a * ~c * ~m + ~b * ~d * ~n)) - ~d * (~a * ~B * (~m - ~n) + ~A * ~b * (~m + ~n + 1)) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (LtQ(~m, -1 / 2) && (GtQ(~n, 0) && (IntegerQ(2 * ~m) && (IntegerQ(2 * ~n) || EqQ(~c, 0))))))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (~b * (~A * ~b - ~a * ~B) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~a * ~f * (2 * ~m + 1) * (~b * ~c - ~a * ~d)) + (1 / (~a * (2 * ~m + 1) * (~b * ~c - ~a * ~d))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n * Simp(~B * (~a * ~c * ~m + ~b * ~d * (~n + 1)) + ~A * (~b * ~c * (~m + 1) - ~a * ~d * (2 * ~m + ~n + 2)) + ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 2) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (LtQ(~m, -1 / 2) && (Not(GtQ(~n, 0)) && (IntegerQ(2 * ~m) && (IntegerQ(2 * ~n) || EqQ(~c, 0))))))))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-2 * ~b * ~B * cos(~e + ~f * ~x) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (2 * ~n + 3) * sqrt(~a + ~b * sin(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && EqQ(~A * ~b * ~d * (2 * ~n + 3) - ~B * (~b * ~c - 2 * ~a * ~d * (~n + 1)), 0))))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-((~b) ^ 2) * (~B * ~c - ~A * ~d) * cos(~e + ~f * ~x) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * (~b * ~c + ~a * ~d) * sqrt(~a + ~b * sin(~e + ~f * ~x))) + ((~A * ~b * ~d * (2 * ~n + 3) - ~B * (~b * ~c - 2 * ~a * ~d * (~n + 1))) / (2 * ~d * (~n + 1) * (~b * ~c + ~a * ~d))) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && LtQ(~n, -1))))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-2 * ~b * ~B * cos(~e + ~f * ~x) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (2 * ~n + 3) * sqrt(~a + ~b * sin(~e + ~f * ~x))) + ((~A * ~b * ~d * (2 * ~n + 3) - ~B * (~b * ~c - 2 * ~a * ~d * (~n + 1))) / (~b * ~d * (2 * ~n + 3))) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && Not(LtQ(~n, -1)))))
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~(c') + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~A * ~b - ~a * ~B) / ~b) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) + (~B / ~b) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sqrt(~c + ~d * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~B) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~f * (~m + ~n + 1)) + (1 / (~b * (~m + ~n + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n - 1) * Simp(~A * ~b * ~c * (~m + ~n + 1) + ~B * (~a * ~c * ~m + ~b * ~d * ~n) + (~A * ~b * ~d * (~m + ~n + 1) + ~B * (~a * ~d * ~m + ~b * ~c * ~n)) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (GtQ(~n, 0) && (IntegerQ(~n) || EqQ(~m + 1 / 2, 0))))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~B * ~c - ~A * ~d) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2)) + (1 / (~b * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1) * Simp((~A * (~a * ~d * ~m + ~b * ~c * (~n + 1)) - ~B * (~a * ~c * ~m + ~b * ~d * (~n + 1))) + ~b * (~B * ~c - ~A * ~d) * (~m + ~n + 2) * sin(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (LtQ(~n, -1) && (IntegerQ(~n) || EqQ(~m + 1 / 2, 0))))))
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~A * ~b - ~a * ~B) / (~b * ~c - ~a * ~d)) * Antiderivative(1 / sqrt(~a + ~b * sin(~e + ~f * ~x)), ~x) + ((~B * ~c - ~A * ~d) / (~b * ~c - ~a * ~d)) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (~c + ~d * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x))) / (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)), ~x) => (~B / ~d) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) - ((~B * ~c - ~A * ~d) / ~d) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m / (~c + ~d * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && NeQ(~m + 1 / 2, 0))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~A * ~b - ~a * ~B) / ~b) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) + (~B / ~b) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && NeQ(~A * ~b + ~a * ~B, 0))))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ 2 * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~B * ~c - ~A * ~d) * (~b * ~c - ~a * ~d) ^ 2 * cos(~e + ~f * ~x) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~d) ^ 2 * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2)) - (1 / ((~d) ^ 2 * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2))) * Antiderivative((~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1) * Simp((~d * (~n + 1) * (~B * (~b * ~c - ~a * ~d) ^ 2 - ~A * ~d * (((~a) ^ 2 * ~c + (~b) ^ 2 * ~c) - 2 * ~a * ~b * ~d)) - ((~B * ~c - ~A * ~d) * ((~a) ^ 2 * (~d) ^ 2 * (~n + 2) + (~b) ^ 2 * ((~c) ^ 2 + (~d) ^ 2 * (~n + 1))) + 2 * ~a * ~b * ~d * (~A * ~c * ~d * (~n + 2) - ~B * ((~c) ^ 2 + (~d) ^ 2 * (~n + 1)))) * sin(~e + ~f * ~x)) - (~b) ^ 2 * ~B * ~d * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && LtQ(~n, -1))))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-((~b * ~c - ~a * ~d)) * (~B * ~c - ~A * ~d) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2)) + (1 / (~d * (~n + 1) * ((~c) ^ 2 - (~d) ^ 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 2) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1) * Simp(~b * (~b * ~c - ~a * ~d) * (~B * ~c - ~A * ~d) * (~m - 1) + ~a * ~d * ((~a * ~A * ~c + ~b * ~B * ~c) - (~A * ~b + ~a * ~B) * ~d) * (~n + 1) + (~b * (~b * ~d * (~B * ~c - ~A * ~d) + ~a * (~A * ~c * ~d + ~B * ((~c) ^ 2 - 2 * (~d) ^ 2))) * (~n + 1) - ~a * (~b * ~c - ~a * ~d) * (~B * ~c - ~A * ~d) * (~n + 2)) * sin(~e + ~f * ~x) + ~b * (~d * ((~A * ~b * ~c + ~a * ~B * ~c) - ~a * ~A * ~d) * (~m + ~n + 1) - ~b * ~B * ((~c) ^ 2 * ~m + (~d) ^ 2 * (~n + 1))) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (GtQ(~m, 1) && LtQ(~n, -1)))))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~b) * ~B * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n + 1)) + (1 / (~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 2) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n * Simp((~a) ^ 2 * ~A * ~d * (~m + ~n + 1) + ~b * ~B * (~b * ~c * (~m - 1) + ~a * ~d * (~n + 1)) + (~a * ~d * (2 * ~A * ~b + ~a * ~B) * (~m + ~n + 1) - ~b * ~B * (~a * ~c - ~b * ~d * (~m + ~n))) * sin(~e + ~f * ~x) + ~b * (~A * ~b * ~d * (~m + ~n + 1) - ~B * (~b * ~c * ~m - ~a * ~d * (2 * ~m + ~n))) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (GtQ(~m, 1) && Not(IGtQ(~n, 1) && (Not(IntegerQ(~m)) || EqQ(~a, 0) && NeQ(~c, 0)))))))
    @apply_utils Antiderivative((sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x))) / (~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2), ~x) => ((~B * ~d) / (~b) ^ 2) * Antiderivative(sqrt(~b * sin(~e + ~f * ~x)) / sqrt(~c + ~d * sin(~e + ~f * ~x)), ~x) + Antiderivative((~A * ~c + (~B * ~c + ~A * ~d) * sin(~e + ~f * ~x)) / ((~b * sin(~e + ~f * ~x)) ^ (3 / 2) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)
    @apply_utils Antiderivative((sqrt(~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x))) / (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2), ~x) => (~B / ~b) * Antiderivative(sqrt(~c + ~d * sin(~e + ~f * ~x)) / sqrt(~a + ~b * sin(~e + ~f * ~x)), ~x) + ((~A * ~b - ~a * ~B) / ~b) * Antiderivative(sqrt(~c + ~d * sin(~e + ~f * ~x)) / (~a + ~b * sin(~e + ~f * ~x)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (2 * (~A * ~b - ~a * ~B) * cos(~e + ~f * ~x)) / (~f * ((~a) ^ 2 - (~b) ^ 2) * sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~d * sin(~e + ~f * ~x))) + (~d / ((~a) ^ 2 - (~b) ^ 2)) * Antiderivative(((~A * ~b - ~a * ~B) + (~a * ~A - ~b * ~B) * sin(~e + ~f * ~x)) / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~d * sin(~e + ~f * ~x)) ^ (3 / 2)), ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~A, ~B], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~A + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-2 * ~A * (~c - ~d) * tan(~e + ~f * ~x)) / (~f * ~b * (~c) ^ 2)) * Rt((~c + ~d) / ~b, 2) * sqrt((~c * (1 + csc(~e + ~f * ~x))) / (~c - ~d)) * sqrt((~c * (1 - csc(~e + ~f * ~x))) / (~c + ~d)) * Elliptic.E(asin((sqrt(~c + ~d * sin(~e + ~f * ~x)) / sqrt(~b * sin(~e + ~f * ~x))) / Rt((~c + ~d) / ~b, 2)), -((~c + ~d)) / (~c - ~d)) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (EqQ(~A, ~B) && PosQ((~c + ~d) / ~b)))
    @apply_utils Antiderivative((~A + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (-(sqrt(-(~b) * sin(~e + ~f * ~x))) / sqrt(~b * sin(~e + ~f * ~x))) * Antiderivative((~A + ~B * sin(~e + ~f * ~x)) / ((-(~b) * sin(~e + ~f * ~x)) ^ (3 / 2) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (EqQ(~A, ~B) && NegQ((~c + ~d) / ~b)))
    @apply_utils Antiderivative((~A + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-2 * ~A * (~c - ~d) * (~a + ~b * sin(~e + ~f * ~x))) / (~f * (~b * ~c - ~a * ~d) ^ 2 * Rt((~a + ~b) / (~c + ~d), 2) * cos(~e + ~f * ~x))) * sqrt(((~b * ~c - ~a * ~d) * (1 + sin(~e + ~f * ~x))) / ((~c - ~d) * (~a + ~b * sin(~e + ~f * ~x)))) * sqrt((-((~b * ~c - ~a * ~d)) * (1 - sin(~e + ~f * ~x))) / ((~c + ~d) * (~a + ~b * sin(~e + ~f * ~x)))) * Elliptic.E(asin((Rt((~a + ~b) / (~c + ~d), 2) * sqrt(~c + ~d * sin(~e + ~f * ~x))) / sqrt(~a + ~b * sin(~e + ~f * ~x))), ((~a - ~b) * (~c + ~d)) / ((~a + ~b) * (~c - ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (EqQ(~A, ~B) && PosQ((~a + ~b) / (~c + ~d))))))
    @apply_utils Antiderivative((~A + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (sqrt(-(~c) - ~d * sin(~e + ~f * ~x)) / sqrt(~c + ~d * sin(~e + ~f * ~x))) * Antiderivative((~A + ~B * sin(~e + ~f * ~x)) / ((~a + ~b * sin(~e + ~f * ~x)) ^ (3 / 2) * sqrt(-(~c) - ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (EqQ(~A, ~B) && NegQ((~a + ~b) / (~c + ~d))))))
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~A - ~B) / (~a - ~b)) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) - ((~A * ~b - ~a * ~B) / (~a - ~b)) * Antiderivative((1 + sin(~e + ~f * ~x)) / ((~a + ~b * sin(~e + ~f * ~x)) ^ (3 / 2) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && NeQ(~A, ~B))))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => ((~B * ~a - ~A * ~b) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~f * (~m + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / ((~m + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (~n - 1) * Simp((~c * (~a * ~A - ~b * ~B) * (~m + 1) + ~d * ~n * (~A * ~b - ~a * ~B) + (~d * (~a * ~A - ~b * ~B) * (~m + 1) - ~c * (~A * ~b - ~a * ~B) * (~m + 2)) * sin(~e + ~f * ~x)) - ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 2) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (LtQ(~m, -1) && GtQ(~n, 0)))))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-((~A * (~b) ^ 2 - ~a * ~b * ~B)) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ (1 + ~n)) / (~f * (~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n * Simp(((~a * ~A - ~b * ~B) * (~b * ~c - ~a * ~d) * (~m + 1) + ~b * ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 2) + (~A * ~b - ~a * ~B) * (~a * ~d * (~m + 1) - ~b * ~c * (~m + 2)) * sin(~e + ~f * ~x)) - ~b * ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 3) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (RationalQ(~m) && (~m < -1 && (EqQ(~a, 0) && (IntegerQ(~m) && Not(IntegerQ(~n))) || Not(IntegerQ(2 * ~n) && (LtQ(~n, -1) && (IntegerQ(~n) && Not(IntegerQ(~m)) || EqQ(~a, 0))))))))))
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / ((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~A * ~b - ~a * ~B) / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (~a + ~b * sin(~e + ~f * ~x)), ~x) + ((~B * ~c - ~A * ~d) / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (~c + ~d * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x))) / (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)), ~x) => (~B / ~d) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) - ((~B * ~c - ~A * ~d) / ~d) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m / (~c + ~d * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => (-2 * ~B * cos(~e + ~f * ~x) * sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~f * (2 * ~n + 3)) + (1 / (2 * ~n + 3)) * Antiderivative(((~c + ~d * sin(~e + ~f * ~x)) ^ (~n - 1) / sqrt(~a + ~b * sin(~e + ~f * ~x))) * Simp(~a * ~A * ~c * (2 * ~n + 3) + ~B * (~b * ~c + 2 * ~a * ~d * ~n) + (~B * (~a * ~c + ~b * ~d) * (2 * ~n + 1) + ~A * (~b * ~c + ~a * ~d) * (2 * ~n + 3)) * sin(~e + ~f * ~x) + (~A * ~b * ~d * (2 * ~n + 3) + ~B * (~a * ~d + 2 * ~b * ~c * ~n)) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && EqQ((~n) ^ 2, 1 / 4))))
    @apply_utils Antiderivative((~A + ~(B') * sin(~(e') + ~(f') * ~x)) / (sqrt(sin(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x))), ~x) => ((4 * ~A) / (~f * sqrt(~a + ~b))) * Elliptic.Pi(-1, -(asin(cos(~e + ~f * ~x) / (1 + sin(~e + ~f * ~x)))), -((~a - ~b)) / (~a + ~b)) <-- FreeQ([~a, ~b, ~e, ~f, ~A, ~B], ~x) && (GtQ(~b, 0) && (GtQ((~b) ^ 2 - (~a) ^ 2, 0) && EqQ(~A, ~B)))
    @apply_utils Antiderivative((~A + ~(B') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~d * sin(~(e') + ~(f') * ~x))), ~x) => (sqrt(sin(~e + ~f * ~x)) / sqrt(~d * sin(~e + ~f * ~x))) * Antiderivative((~A + ~B * sin(~e + ~f * ~x)) / (sqrt(sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~d, ~A, ~B], ~x) && (GtQ(~b, 0) && (GtQ((~b) ^ 2 - (~a) ^ 2, 0) && EqQ(~A, ~B)))
    @apply_utils Antiderivative((~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~(c') + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (~B / ~d) * Antiderivative(sqrt(~c + ~d * sin(~e + ~f * ~x)) / sqrt(~a + ~b * sin(~e + ~f * ~x)), ~x) - ((~B * ~c - ~A * ~d) / ~d) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)), ~x) => Unintegrable((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n * (~A + ~B * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    #= Antiderivative(((~a)+(~b')*sin((~e')+(~f')*(~x)))^(~m)*((~c)+(~d')*sin((~e')+(~f')*(~x)))^(~n)*((~A')+(~B) .*sin((~e')+(~f')*(~x)))^(~p),~x) :=  a^m*c^m*Antiderivative(cos(e+f*x)^(2*m)*(c+d*sin(e+f*x))^(n-m)*(A+B*sin(e+f*x))^ p,x) <-- FreeQ([a,b,c,d,e,f,A,B,n,p],x) && EqQ(b*c+a*d,0) &&  EqQ(a^2-b^2,0) && IntegerQ(m) && Not(IntegerQ(n) && (LtQ(m,0) && GtQ(n,0) || LtQ(0,n,m) ||  LtQ(m,n,0))) =#
    #= Antiderivative(((~a)+(~b')*cos((~e')+(~f')*(~x)))^(~m)*((~c)+(~d')*cos((~e')+(~f')*(~x)))^(~n)*((~A')+(~B) .*cos((~e')+(~f')*(~x)))^(~p),~x) :=  a^m*c^m*Antiderivative(sin(e+f*x)^(2*m)*(c+d*cos(e+f*x))^(n-m)*(A+B*cos(e+f*x))^ p,x) <-- FreeQ([a,b,c,d,e,f,A,B,n,p],x) && EqQ(b*c+a*d,0) &&  EqQ(a^2-b^2,0) && IntegerQ(m) && Not(IntegerQ(n) && (LtQ(m,0) && GtQ(n,0) || LtQ(0,n,m) ||  LtQ(m,n,0))) =#
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * sin(~(e') + ~(f') * ~x)) ^ ~p, ~x) => ((sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative((~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ (~n - 1 / 2) * (~A + ~B * ~x) ^ ~p, ~x), ~x, sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * cos(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * cos(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * cos(~(e') + ~(f') * ~x)) ^ ~p, ~x) => ((-(sqrt(~a + ~b * cos(~e + ~f * ~x))) * sqrt(~c + ~d * cos(~e + ~f * ~x))) / (~f * sin(~e + ~f * ~x))) * Subst(Antiderivative((~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ (~n - 1 / 2) * (~A + ~B * ~x) ^ ~p, ~x), ~x, cos(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.3.1 (a+b sin)^m (c+d sin)^n (A+B sin).jl


text:
using SimpleContainerGenerator
using Test

import Pkg
import Random: randstring

if get(ENV, "JULIA_VERSION", nothing) == "nightly"
    const JULIA_VERSION = "nightly"
else
    const JULIA_VERSION = SimpleContainerGenerator._default_julia_version()
end

@testset "SimpleContainerGenerator.jl" begin
    @testset "assert.jl" begin
        Test.@test SimpleContainerGenerator.always_assert(true, "") == nothing
        Test.@test_throws SimpleContainerGenerator.AlwaysAssertionError SimpleContainerGenerator.always_assert(false, "")
    end

    @testset "default_values.jl" begin
        Test.@test SimpleContainerGenerator._default_docker_parent_image() isa AbstractString
        Test.@test SimpleContainerGenerator._default_julia_cpu_target() isa AbstractString
        Test.@test SimpleContainerGenerator._default_julia_version() isa Union{AbstractString, VersionNumber}
        Test.@test SimpleContainerGenerator._default_packagecompiler_installation_command() isa AbstractString
        Test.@test SimpleContainerGenerator._default_simplecontainergenerator_installation_command() isa AbstractString
        Test.@test SimpleContainerGenerator._fastest_nonportable_julia_cpu_target() isa AbstractString
        Test.@test SimpleContainerGenerator._predictmd_apt() isa AbstractVector
        Test.@test SimpleContainerGenerator._recommended_julia_cpu_target() isa AbstractString
    end

    @testset "pkg_dir.jl" begin
        Test.@test SimpleContainerGenerator.pkg_dir_simplecontainergenerator() isa AbstractString
    end

    @testset "public.jl" begin
        test_cases = ["Crayons",
                      ["Crayons"],
                      (name = "Crayons",),
                      [(name = "Crayons",)],
                      Dict(:name => "Crayons"),
                      [Dict(:name => "Crayons")]]
        for test_case in test_cases
            SimpleContainerGenerator.with_temp_dir() do tmp_dir
                @test !isfile("Dockerfile")
                tests_must_pass = String["Crayons"]
                SimpleContainerGenerator.create_dockerfile(test_case;
                                                           tests_must_pass = tests_must_pass)
                @test isfile("Dockerfile")
            end
            for make_sysimage in [true, false]
                SimpleContainerGenerator.with_temp_dir() do tmp_dir
                    @test !isfile("Dockerfile")
                    tests_must_pass = String["Crayons"]
                    SimpleContainerGenerator.create_dockerfile(test_case;
                                                               make_sysimage = make_sysimage,
                                                               tests_must_pass = tests_must_pass)
                    @test isfile("Dockerfile")
                end
            end
        end
    end

    @testset "docker.jl" begin
        @test SimpleContainerGenerator._generate_apt_install_command(SimpleContainerGenerator.Config(; override_default_apt = String[])) == ""
    end

    @testset "julia.jl" begin
        @test SimpleContainerGenerator._get_julia_url("nightly") == ("https://julialangnightlies-s3.julialang.org/bin/linux/x64/julia-latest-linux64.tar.gz", "https://julialangnightlies-s3.julialang.org/bin/linux/x64/julia-latest-linux64.tar.gz.asc")
        @test SimpleContainerGenerator._get_julia_url("1.4.1") == ("https://julialang-s3.julialang.org/bin/linux/x64/1.4/julia-1.4.1-linux-x86_64.tar.gz", "https://julialang-s3.julialang.org/bin/linux/x64/1.4/julia-1.4.1-linux-x86_64.tar.gz.asc")
    end

    if get(ENV, "SIMPLECONTAINERGENERATOR_TESTS", "") == "all"
        @testset "docker build" begin
            test_cases = [([Dict(:name => "Example")], ["import Example", "using Example", "import Pkg; Pkg.test(string(:Example))"])]
            for (pkgs, test_eval_exprs) in test_cases
                for make_sysimage in [true, false]
                    SimpleContainerGenerator.with_temp_dir() do tmp_dir
                        @test !isfile("Dockerfile")
                        tests_must_pass = String["Example"]
                        SimpleContainerGenerator.create_dockerfile(pkgs;
                                                                   julia_version = JULIA_VERSION,
                                                                   make_sysimage = make_sysimage,
                                                                   tests_must_pass = tests_must_pass)
                        @test isfile("Dockerfile")
                        image = lowercase("test_simplecontainergenerator_runtests_$(randstring(16))")
                        p_build = run(`docker build -t $(image) .`)
                        wait(p_build)
                        @test success(p_build)
                        for test_eval_expr in test_eval_exprs
                            p_test = run(`docker run $(image) "/usr/bin/julia -e '$(test_eval_expr)'"`)
                            wait(p_test)
                            @test success(p_test)
                        end
                    end
                end
            end
        end
    end
end


text:
function nlp_007_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - infeasible model

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)

    @NLconstraint(model, y == exp(x))
    @constraint(model, x == y^2)

    optimize!(model)

    return check_status(
        model,
        INFEASIBLE_PROBLEM,
        termination_target,
        primal_target,
    )
end


text:
# Reduction functions related to statistics

###################
#
#  Varm & Stdm
#
###################

# varm

function varm(x::ContiguousRealArray, mu::Real)
    !isempty(x) || error("varm: empty array is not allowed.")
    n = length(x)

    @inbounds s2 = abs2(x[1] - mu)
    for i = 2:n
        @inbounds s2 += abs2(x[i] - mu)
    end
    s2 / (n - 1)
end

function _varm_eachcol!{R<:Real}(m::Int, n::Int, dst::ContiguousArray{R}, x, mu)
    c = inv(m - 1)
    o = 0
    for j = 1:n
        s2 = zero(R)
        @inbounds mu_j = mu[j]
        for i = 1:m
            @inbounds s2 += abs2(x[o + i] - mu_j)
        end
        @inbounds dst[j] = s2 * c
        o += m
    end
end

function _varm_eachrow!(m::Int, n::Int, dst::ContiguousRealArray, x, mu)
    o::Int = 0
    for i = 1:m
        @inbounds dst[i] = abs2(x[o + i] - mu[i])
    end
    o += m
    for j = 2:n-1
        for i = 1:m
            @inbounds dst[i] += abs2(x[o + i] - mu[i])
        end
        o += m
    end
    c = inv(n - 1)
    for i = 1:m
        @inbounds v = dst[i] + abs2(x[o + i] - mu[i])
        dst[i] = v * c
    end
end

function varm!(dst::ContiguousRealArray, x::ContiguousRealArray, mu::ContiguousRealArray, dim::Int)
    !isempty(x) || error("varm!: empty array is not allowed.")

    nd = ndims(x)
    1 <= dim <= nd || error("varm: invalid value of dim.")

    shp = size(x)
    rlen = prod(Base.reduced_dims(shp, dim))
    length(dst) == length(mu) == rlen || error("Inconsistent argument dimensions.")

    if dim == 1
        m = shp[1]
        n = succ_length(shp, dim)
        _varm_eachcol!(m, n, dst, x, mu)

    else
        m = prec_length(shp, dim)
        n = shp[dim]
        k = succ_length(shp, dim)

        _varm_eachrow!(m, n, dst, x, mu)
        if k > 1
            mn = m * n
            ro = m
            ao = mn
            for l = 2 : k
                _varm_eachrow!(m, n,
                    offset_view(dst, ro, m), offset_view(x, ao, m, n),
                    offset_view(mu, ro, m))

                ro += m
                ao += mn
            end
        end
    end
    dst
end

function varm(x::ContiguousRealArray, mu::ContiguousRealArray, dim::Int)
    rsiz = Base.reduced_dims(size(x), dim)
    length(mu) == prod(rsiz) || error("Inconsistent argument dimensions.")
    R = fptype(promote_type(eltype(x), eltype(mu)))
    varm!(Array(R, rsiz), x, mu, dim)
end

# stdm

stdm(x::ContiguousRealArray, mu::Real) = sqrt(varm(x, mu))
stdm(x::ContiguousRealArray, mu::ContiguousArray, dim::Int) = sqrt!(varm(x, mu, dim))

function stdm!(dst::ContiguousRealArray, x::ContiguousRealArray, mu::ContiguousRealArray, dim::Int)
    sqrt!(varm!(dst, x, mu, dim))
end

###################
#
#  Var & Std
#
###################

var(x::ContiguousRealArray) = varm(x, mean(x))

function var!(dst::ContiguousRealArray, x::ContiguousRealArray, dim::Int)
    !isempty(dst) || error("var: empty array is not allowed.")
    nd = ndims(x)
    1 <= dim <= nd || error("var: invalid value of dim.")
    shp = size(x)

    if dim == 1
        m = shp[1]
        n = succ_length(shp, 1)
        ao = 0
        for j in 1 : n
            dst[j] = var(offset_view(x, ao, m))
            ao += m
        end
    else
        varm!(dst, x, mean(x, dim), dim)
    end
    dst
end

function var(x::ContiguousRealArray, dim::Int)
    var!(Array(fptype(eltype(x)), Base.reduced_dims(size(x), dim)), x, dim)
end

# std

std(x::ContiguousRealArray) = sqrt(var(x))
std(x::ContiguousRealArray, dim::Int) = sqrt!(var(x, dim))
std!(dst::ContiguousRealArray, x::ContiguousRealArray, dim::Int) = sqrt!(var!(dst, x, dim))

###################
#
#  LogFunsumexp
#
###################

function logsumexp(x::ContiguousRealArray)
    !isempty(x) || error("logsumexp: empty array not allowed.")
    u = maximum(x)
    log(sumfdiff(ExpFun(), x, u)) + u
end

function _logsumexp_eachcol!(m::Int, n::Int, dst::ContiguousRealArray, x::ContiguousRealArray)
    o = 0
    for j in 1 : n
        # compute max
        u = x[o + 1]
        for i in 2 : m
            @inbounds xi = x[o + i]
            if xi > u
                u = xi
            end
        end

        # sum exp
        @inbounds s = exp(x[o + 1] - u)
        for i in 2 : m
            @inbounds s += exp(x[o + i] - u)
        end

        # compute log
        dst[j] = log(s) + u
        o += m
    end
end

function _logsumexp_eachrow!(m::Int, n::Int, dst::ContiguousRealArray,
    u::ContiguousRealArray, x::ContiguousRealArray)

    # compute max
    for i in 1 : m
        @inbounds u[i] = x[i]
    end
    o = m
    for j in 2 : n
        for i in 1 : m
            @inbounds ui = u[i]
            @inbounds xi = x[o+i]
            if xi > ui
                @inbounds u[i] = xi
            end
        end
        o += m
    end

    # sum exp
    for i in 1 : m
        @inbounds dst[i] = exp(x[i] - u[i])
    end
    o = m
    for j in 2 : n
        for i in 1 : m
            @inbounds dst[i] += exp(x[o + i] - u[i])
        end
        o += m
    end

    # compute log
    for i in 1 : m
        @inbounds dst[i] = log(dst[i]) + u[i]
    end
end

function logsumexp!{R<:Real,T<:Real}(dst::ContiguousArray{R}, x::ContiguousRealArray{T}, dim::Int)
    !isempty(x) || error("logsumexp!: empty array not allowed.")
    nd = ndims(x)
    1 <= dim <= nd || error("logsumexp!: invalid value of dim.")
    shp = size(x)

    if dim == 1
        m = shp[1]
        n = succ_length(shp, dim)
        _logsumexp_eachcol!(m, n, dst, x)
    else
        m = prec_length(shp, dim)
        n = shp[dim]
        k = succ_length(shp, dim)

        u = Array(T, m)
        _logsumexp_eachrow!(m, n, dst, u, x)
        if k > 1
            mn = m * n
            ro = m
            ao = mn
            for l = 2 : k
                _logsumexp_eachrow!(m, n, offset_view(dst, ro, m), u, offset_view(x, ao, m, n))
                ro += m
                ao += mn
            end
        end
    end
    dst
end

function logsumexp{T<:Real}(x::ContiguousArray{T}, dim::Int)
    logsumexp!(Array(fptype(T), Base.reduced_dims(size(x), dim)), x, dim)
end

###################
#
#  Softmax
#
###################

function softmax!(dst::ContiguousRealArray, x::ContiguousRealArray)
    !isempty(x) || error("softmax!: empty array is not allowed.")
    n = length(x)
    length(dst) == n || error("Inconsistent argument dimensions.")

    u = maximum(x)
    @inbounds s = dst[1] = exp(x[1] - u)
    for i in 2 : n
        @inbounds s += (dst[i] = exp(x[i] - u))
    end
    c = inv(s)
    for i in 1 : n
        @inbounds dst[i] *= c
    end
    dst
end

softmax(x::ContiguousArray) = softmax!(Array(fptype(eltype(x)), size(x)), x)

function _softmax_eachcol!(m::Int, n::Int, dst::ContiguousRealArray, x::ContiguousRealArray)
    o = 0
    for j in 1 : n
        softmax!(offset_view(dst, o, m), offset_view(x, o, m))
        o += m
    end
end

function _softmax_eachrow!(m::Int, n::Int, dst::ContiguousRealArray, u::ContiguousRealArray, x::ContiguousRealArray)

    # compute max
    for i in 1 : m
        @inbounds u[i] = x[i]
    end
    o = m
    for j in 2 : n
        for i in 1 : m
            @inbounds ui = u[i]
            @inbounds xi = x[o + i]
            if xi > ui
                @inbounds u[i] = xi
            end
        end
        o += m
    end

    # compute sum
    s = view(u, m+1:2*m)

    for i in 1 : m
        @inbounds s[i] = dst[i] = exp(x[i] - u[i])
    end
    o = m

    for j in 2 : n
        for i in 1 : m
            @inbounds s[i] += (dst[o + i] = exp(x[o + i] - u[i]))
        end
        o += m
    end

    rcp!(s)
    o = 0
    for j in 1 : n
        for i in 1 : m
            @inbounds dst[o + i] .*= s[i]
        end
        o += m
    end
end

function softmax!{T<:Real}(dst::ContiguousRealArray, x::ContiguousArray{T}, dim::Int)
    !isempty(x) || error("softmax!: empty array is not allowed.")
    nd = ndims(x)
    1 <= dim <= nd || error("softmax!: invalid value for the dim argument.")
    shp = size(x)

    if dim == 1
        m = shp[1]
        n = succ_length(shp, dim)
        _softmax_eachcol!(m, n, dst, x)

    else
        m = prec_length(shp, dim)
        n = shp[dim]
        k = succ_length(shp, dim)

        u = Array(fptype(T), 2*m)
        _softmax_eachrow!(m, n, dst, u, x)
        if k > 1
            mn = m * n
            o = mn
            for l = 2 : k
                _softmax_eachrow!(m, n, offset_view(dst, o, m, n), u, offset_view(x, o, m, n))
                o += mn
            end
        end
    end
    dst
end

softmax(x::ContiguousRealArray, dim::Int) = softmax!(Array(fptype(eltype(x)), size(x)), x, dim)


text:
## ExoplanetsSysSim/src/limb_darkening.jl
## (c) 2015 Eric B. Ford

abstract LimbDarkeningParamAbstract

immutable LimbDarkeningParamQuadratic <: LimbDarkeningParamAbstract
  u1::Float64
  u2::Float64
  # Demo of how to enforce constraints on parameter when constructing class
  # TODO SCI DETAIL: Repalce with sensible limits on LD params
  LimbDarkeningParamQuadratic(_u1::Real, _u2::Real ) = (!( (-2.0<=_u1<=2.0) && (-2.0<=_u2<=2.0) )) ? error(string("Invalid quadratic limb darkening parameters: ",_u1," & ", _u2)) : new([_u1,_u2])
end


text:
struct ExperimentSummary
    _outputFile::String
    _batchSize::Int16
    _printDuringExperiment::Bool

    _fitnessValues::Array{Float64}
    _bestFitnessValues::Array{Float64}
    _meanFitness::Array{Float64}
    _stdFitness::Array{Float64}
    _bestIndividuals::Array{Any}
end # struct

"""
    getOutputFile(summary::ExperimentSummary)::String

Obtains the output file for the experiment summary.
"""
getOutputFile(summary::ExperimentSummary)::String = summary._outputFile
# function

"""
    getBatchSize(summary::ExperimentSummary)::Int16

Obtains the batch size for results collection.
"""
getBatchSize(summary::ExperimentSummary)::Int16 = summary._batchSize
# function

"""
    printDuringExperiment(summary::ExperimentSummary)::Bool

Checks if printing during the experiment is enabled.
"""
printDuringExperiment(summary::ExperimentSummary)::Bool = summary._printDuringExperiment
# function

"""
    displayFitness(summary::ExperimentSummary)

Checks if the user wants to display the fitness values of the individuals
throughout the generations.
"""
displayFitness(summary::ExperimentSummary) = !(isempty(summary._fitnessValues))
# function

"""
    displayBestFitness(summary::ExperimentSummary)

Checks if the user wants to display the best fitness values of the individuals
throughout the generations.
"""
displayBestFitness(summary::ExperimentSummary) = !(isempty(summary._bestFitnessValues))
# function

"""
    displayMeanFitness(summary::ExperimentSummary)

Checks if the user wants to display the mean of the fitness values of the individuals
throughout the generations.
"""
displayMeanFitness(summary::ExperimentSummary) = !(isempty(summary._meanFitness))
# function

"""
    displaySTDFitness(summary::ExperimentSummary)

Checks if the user wants to display the standard deviation of the fitness values
of the individuals throughout the generations.
"""
displaySTDFitness(summary::ExperimentSummary) = !(isempty(summary._stdFitness))
# function

"""
    somethingToDisplay(summary::ExperimentSummary)

Checks if there is something to display from the summary.
"""
function somethingToDisplay(summary::ExperimentSummary)
    return getBatchSize(summary) >= 0 && (displayFitness(summary) ||
           displayBestFitness(summary) || displayMeanFitness(summary) ||
           displaySTDFitness(summary))
end # function

"""
    anyToDisplay(summary::ExperimentSummary)

Returns the first thing that can be displayed.
"""
function displayAnything(summary::ExperimentSummary)
    if displayFitness(summary)
        return summary._fitnessValues
    elseif displayBestFitness(summary)
        return summary._bestFitnessValues
    elseif displayMeanFitness(summary)
        return summary._meanFitness
    elseif displaySTDFitness(summary)
        return summary._stdFitness
    end
end # function

"""
    saveFitness(summary::ExperimentSummary, fitnesses::Array{Float64},
                currIter::Int16)

Saves the fitness values of the individuals of a generation.
"""
function saveFitness(summary::ExperimentSummary, fitnesses::Array{Float64},
                     currIter::Int, globalFitnesses::Array{Float64} = Array{Float64}(undef, 0))

    index = div(currIter, getBatchSize(summary))

    if !(isempty(globalFitnesses))
        # CAMBIAR ESTO CUANDO SE MODIFIQUE LA ESTRUCTURA DE FILAS Y COLUMNAS DE LOS FITNESS
        # CAMBIARÍA A: summary._fitnessValues[:, i, index] = vcat(fitnesses[:, i], globalFitnesses[i])
        for i = eachindex(globalFitnesses)
            summary._fitnessValues[:, i, index] = vcat(fitnesses[i, :], globalFitnesses[i])
        end
    else
        # CAMBIAR ESTO CUANDO SE MODIFIQUE LA ESTRUCTURA DE FILAS Y COLUMNAS DE LOS FITNESS
        # CAMBIARÍA A: summary._fitnessValues[:, :, index] = fitnesses
        summary._fitnessValues[:, :, index] = transpose(fitnesses)
    end

    return nothing
end # function

"""
    saveBestFitness(summary::ExperimentSummary, fitnesses::Array{Float64},
                    currIter::Int16)

Saves the best fitness values of a generation.
"""
function saveBestFitness(summary::ExperimentSummary, bestFitness::Union{Array{Float64}, Float64},
                         representation::Union{Array, Expr}, currIter::Integer,
                         bestGlobalFitness::Float64 = NaN)

    index = div(currIter, getBatchSize(summary))
    if !(isnan(bestGlobalFitness))
        if length(axes(bestFitness)) == 0
            summary._bestFitnessValues[:, index] = [bestFitness, bestGlobalFitness]
        else
            #println(size(transpose(bestFitness)))
            #println(size([bestGlobalFitness]))
            summary._bestFitnessValues[:, index] = vcat(bestFitness, [bestGlobalFitness])
        end
    else
        if length(axes(bestFitness)) == 0
            summary._bestFitnessValues[:, index] .= bestFitness
        else
            # CAMBIAR ESTO CUANDO SE MODIFIQUE LA ESTRUCTURA DE FILAS Y COLUMNAS DE LOS FITNESS
            # CAMBIARÍA A: summary._fitnessValues[:, :, index] = fitnesses
            summary._bestFitnessValues[:, index] = transpose(bestFitness)
        end
    end

    summary._bestIndividuals[index] = representation

    return nothing
end # function

"""
    saveMeanFitness(summary::ExperimentSummary, fitnesses::Array{Float64},
                    currIter::Int16)

Saves the mean of fitness values of the individuals of a generation.
"""
function saveMeanFitness(summary::ExperimentSummary, fitness::Array{Float64},
                         currIter::Integer, globalFitnesses::Array{Float64} = Array{Float64}(undef, 0))

    index = div(currIter, getBatchSize(summary))
    # CAMBIAR ESTO CUANDO SE MODIFIQUE LA ESTRUCTURA DE FILAS Y COLUMNAS DE LOS FITNESS
    # CAMBIARÍA A: mean(fitness, dims=2)
    if !(isempty(globalFitnesses))
        if length(axes(fitness)) == 1
            summary._meanFitness[:, index] = [Statistics.mean(fitness), Statistics.mean(globalFitnesses)]
        else
            summary._meanFitness[:, index] = vcat(transpose(Statistics.mean(fitness, dims=1)), [Statistics.mean(globalFitnesses)])
        end
    else
        if length(axes(fitness)) == 1
            summary._meanFitness[:, index] .= Statistics.mean(fitness)
        else
            summary._meanFitness[:, index] = transpose(Statistics.mean(fitness, dims=1))
        end
    end

    return nothing
end # function

"""
    saveSTDFitness(summary::ExperimentSummary, fitnesses::Array{Float64},
                   currIter::Int16)

Saves the standard deviation of fitness values of the individuals of a generation.
"""
function saveSTDFitness(summary::ExperimentSummary, fitness::Array{Float64},
                        currIter::Integer, globalFitnesses::Array{Float64} = Array{Float64}(undef, 0))

    index = div(currIter, getBatchSize(summary))
    # CAMBIAR ESTO CUANDO SE MODIFIQUE LA ESTRUCTURA DE FILAS Y COLUMNAS DE LOS FITNESS
    # CAMBIARÍA A: std(fitness, dims=2)
    if !(isempty(globalFitnesses))
        if length(axes(fitness)) == 1
            if displayMeanFitness(summary)
                summary._stdFitness[:, index] = [Statistics.std(fitness,
                                                                mean=summary._meanFitness[1, index]),
                                                Statistics.std(globalFitnesses,
                                                                mean=summary._meanFitness[2, index])]
            else
                summary._stdFitness[:, index] = [Statistics.std(fitness),
                                                Statistics.std(globalFitnesses)]
            end
        else
            if displayMeanFitness(summary)
                summary._stdFitness[:, index] = vcat(transpose(Statistics.std(fitness,
                                                                              mean=summary._meanFitness[1:end-1, index],
                                                                              dims=1)),
                                                     Statistics.std(globalFitnesses,
                                                                    mean=summary._meanFitness[end:end, index],
                                                                    dims=1))
            else
                summary._stdFitness[:, index] = vcat(transpose(Statistics.std(fitness,
                                                                              dims=1)),
                                                     Statistics.std(globalFitnesses,
                                                                    dims=1))
            end
        end
    else
        if length(axes(fitness)) == 1
            if displayMeanFitness(summary)
                summary._stdFitness[:, index] .= Statistics.std(fitness,
                                                                mean=summary._meanFitness[1, index])
            else
                summary._stdFitness[:, index] .= Statistics.std(fitness)
            end
        else
            if displayMeanFitness(summary)
                summary._stdFitness[:, index] = transpose(Statistics.std(fitness,
                                                                         mean=summary._meanFitness[:, index],
                                                                         dims=1))
            else
                summary._stdFitness[:, index] = transpose(Statistics.std(fitness,
                                                                         dims=1))
            end
        end
    end

    return nothing
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia