text:
using HTTP, CSV, DataFrames, EzXML, Dates, Plots

file = "data\\iom-land-registry.csv"

#functions to download and analyse the land transactions

"""
    getLatestLandTransactionLink()

get the latest link to the iom government land transaction open data

"""
function getLatestLandTransactionLink()
    url = "https://www.gov.im/about-the-government/government/open-data/economy/land-transactions/"
    urlstem = "https://www.gov.im"

    r = HTTP.request("GET", url)

    # read the body into a String
    status, headers, body = r.status, r.headers , String(r.body)

    doc = parsehtml(body)
    html = root(doc)
    xpath = "//a[text()='Land Transactions']"
    link = findfirst(xpath, html)

    urlstem*link["href"]
end

function updatetransactionsdata()
    fulllink = getLatestLandTransactionLink()

    download(fulllink,file)
end

function readlandtransactiondata(update=false)
    if update
        updatetransactionsdata()
    end

    lr = CSV.read(file,DataFrame)
    lr.Market_Value = parse.(Float64,replace.(coalesce.(lr.Market_Value,"0"),","=>""))
    lr.Consideration = parse.(Float64,replace.(coalesce.(lr.Consideration,"0"),","=>""))
    lr.Acquisition_Date = Date.(lr.Acquisition_Date,"d/m/Y")
    lr.Completion_Date = Date.(lr.Completion_Date,"d/m/Y")
    sort!(lr,:Acquisition_Date)
end

lr = readlandtransactiondata()

byprop = groupby(lr,[:SubUnit_Name,:House_Number,:Postcode])
bybuild = groupby(filter([:SubUnit_Name]=>x->!ismissing(x),lr),[:House_Name,:House_Number,:Postcode])
bystreet = groupby(filter([:SubUnit_Name]=>x->!ismissing(x),lr),[:Street_Name,:Locality])
lr.Postcode4 = SubString.(coalesce.(lr.Postcode,"").*"     ",1,5)
lr.Postcode3 = SubString.(coalesce.(lr.Postcode,"").*"   ",1,3)
bypostcode4 = groupby(filter([:SubUnit_Name]=>x->!ismissing(x),lr),[:Postcode4])

filter([:SubUnit_Name]=>x->occursin("Apartment",coalesce(x,"")),lr)
show(ans,allrows=true)
filter([:SubUnit_Name]=>x->occursin("Garage",coalesce(x,"")),lr)
describe(lr[:,[:SubUnit_Name]])
lr.SubUnit_Name_NoDigits = replace.(coalesce.(lr.SubUnit_Name,""),r"\d"=>"")
unique(lr.SubUnit_Name_NoDigits[lr.Market_Value>20000])
show(ans)
@show names(lr);

filt = lr[(lr.Market_Value.>40000) .& (lr.Market_Value.<750000) .& (lr.Completion_Date .>= Date(2016,1,1)),:]
scatter(filt.Completion_Date, filt.Market_Value, title = "Market_Value vs Date",markersize=1,markerstrokewidth=0)
filtonchan = lr[(lr.Market_Value.>40000) .& (lr.Market_Value.<750000)  .& (lr.Town.=="Onchan").& (lr.Completion_Date .>= Date(2016,1,1)) .&,:]
scatter!(filt.Completion_Date, markersize=1,markercolor=4)

["Apartment|Flat"
"Apartments"
"Roadway"
"land"
"Garage"
"Garages"
"and"
"PENTHOUSE"
]

max_stars_repo_name: dannywinrow/iom-govt-data

max_issues_repo_name: dannywinrow/iom-govt-data

max_forks_repo_name: dannywinrow/iom-govt-data


text:
# using LinearAlgebra
# include("FieldCalc.jl")
# include("LitzACResist.jl")
# include("StrayCap.jl")
# include("InductanceCalc.jl")

"""
This function makes an elliptic solenoid and calculates the inductance, capacitance, and resistance
    The inputs are:
    N - Number of turns
    r₁ - one radius of the ellipse
    r₂ - Second radius of ellipse
    L - length of solenoid (axially)
    kwargs:
    StrandDiam - Conductor diamter(m)
    InsDiam = insulation diamter (m),
    ϵr_ins - relative permittivity of insulation
    NPts_Coil  - number of points that each wire loop will be discretized into,
"""
function MakeEllipSolenoid(
    N, #Number of turns
    r₁, #Radius 1 of ellipse
    r₂,#Radius 2 of ellipse
    L; #Length of solenoid in meters
    StrandDiam = 0.3e-3,
    InsDiam = 0.9*1e-3,
    ϵr_ins = 4,
    NumStrands = 1,
    BundleDiam = StrandDiam)

    FieldCentered = EvalField_Centered(N,r₁,r₂,L)

    MeanRad = (r₁ + r₂) / 2
    # ZCoord = [(wind_num - 1) / N * L for wind_num = 1:N]
    # RCoords = MeanRad .* ones(length(ZCoord))
    # Coords = hcat(RCoords[:], ZCoord[:])
    # L_Coil = CalcInductance(Coords)
    L_Coil = SimpleInduct(N,MeanRad,L)
    C_Coil = CL_DAE_Knight(
        L, #Coil length (axial)
        InsDiam, #Coil diamter
        N, #N turns
        1, #external permittivity
        ϵr_ins) #Internal permittivity
    SRF = 1 / (2 * π * sqrt(L_Coil * C_Coil))

    CoilProps, Rac, SkinDepth = LitzWire_ACRes(
        StrandDiam,
        N,
        L,
        MeanRad,
        NumStrands,
        BundleDiam;
        f = 75e3:1000:125e3,
        PlotOn=false)

    return L_Coil, C_Coil, SRF, FieldCentered,Rac
end

function LfromΦMat(ΦMat)

    N = size(ΦMat)[1]
    LMat = zeros(N,N)
    for i in 1:N
        for j in 1:N
            MagΦᵢ = sqrt(sum(ΦMat[i,i,:].^2))
            LMat[i,j] = sum(ΦMat[i,j,:].*ΦMat[i,i,:]) / MagΦᵢ
        end
    end

    return LMat
end

"""
This function makes an elliptic solenoid and calculates the inductance.
    The function calculates the flux through the center of each loop's cross-section
        Lₛ = Φ/I (self inductance = flux per amp)
        Lₘ₂₁ = (Φ₂₁ ⋅ Φ₂₂/|Φ₂₂|) /I
            -mutual inductance from the 1ˢᵗ wire on 2ⁿᵈ is flux from 1 in the 2ⁿᵈ's cross section
            that is aligned with the flux from the second on itself. All per ampere
        Lₘ₂₁ = Lₘ₁₂
        Total L = ∑ᵢ∑ⱼ(Lᵢⱼ)
    The inputs are:
    N - Number of turns
    r₁ - one radius of the ellipse
    r₂ - Second radius of ellipse
    L - length of solenoid (axially)
    kwargs:
    MeasLayers - Number of concentric layers to calc field
                    The field is calculated over concentric ellipses (linearly spaced)
                    This is the number of concentric test point layers
    PtsPerLayer -  Number of test points per layer
    NPts_Coil  - number of points that each wire loop will be discretized into,
"""
function EvalInduct_Biot(
    N, #Number of turns
    r₁, #Radius 1 of ellipse
    r₂,#Radius 2 of ellipse
    L; #Length of solenoid in meters
    MeasLayers = 55,#Number of concentric layers to calc field
    PtsPerLayer =  6,
    NPts_Coil = 500
    )

    A_All = zeros(MeasLayers*PtsPerLayer,3,N)
    Weights = zeros(MeasLayers*PtsPerLayer)
    NPtsPerSlice = MeasLayers*PtsPerLayer
    for wind_num in 1:N
        A_All[:,1,wind_num],A_All[:,2,wind_num],A_All[:,3,wind_num]  =
         MakeEllipTestPoints(r₁,r₂;Center = [0,0,(wind_num - 1) / N * L],NPts=PtsPerLayer,Layers = MeasLayers)

    end
    SliceArea = π*r₁*r₂
    Weights = sum(A_All[:,1:2,1].^2,dims=2) #A matrix that weights each value based on radius
    Weights = Weights / sum(Weights)*NPtsPerSlice  #normalizing

    ΦAll = zeros(N,N,3)

    for Wᵢ in 1:N
        for Aᵢ in 1:N
            DeconstructA =
                        [([A_All[i,1,Aᵢ],A_All[i,2,Aᵢ],A_All[i,3,Aᵢ]]) for i in 1:length(Weights)]
            Field3DSlice = [
                BiotSav(
                    MakeEllip(
                                r₁, #Radius 1 of ellipse
                                r₂;#Radius 2 of ellipse
                                Center = [0, 0, (Wᵢ - 1) / N * L], #The solenoid's axis is in Z
                                NPts = NPts_Coil, #How discretized the windings are
                                ),
                        TestLoc;
                        MinThreshold = 0.001,
                        )
                        for TestLoc in DeconstructA
                        ]

         Bweighted = vcat( [(Field3DSlice[i].*Weights[i]) for i in 1:length(Weights)]...)
          # Bweighted = vcat( [(Field3DSlice[i]) for i in 1:length(Weights)]...)

         ΦAll[Wᵢ,Aᵢ,:] =sum(Bweighted,dims=1).*SliceArea ./NPtsPerSlice
        # ΦAll[Wᵢ,Aᵢ,:] =sum(Bweighted,dims=1) ./NPtsPerSlice
        end
        println("$Wᵢ of $N windings calculated")
    end

    LMat = LfromΦMat(ΦAll)
    LTotal = sum(LMat)*1e6
end

"""
This function makes an elliptic solenoid and calculates the field at the center of the solenoid
        Lₘ₂₁ = Lₘ₁₂
        Total L = ∑ᵢ∑ⱼ(Lᵢⱼ)
    The inputs are:
    N - Number of turns
    r₁ - one radius of the ellipse
    r₂ - Second radius of ellipse
    L - length of solenoid (axially)
    kwargs:

    NPts_Coil  - number of points that each wire loop will be discretized into,
"""
function EvalField_Centered(
    N, #Number of turns
    r₁, #Radius 1 of ellipse
    r₂,#Radius 2 of ellipse
    L; #Length of solenoid in meters
    NPts_Coil=100
    )

            FieldCentered = [
                BiotSav(
                    MakeEllip(
                                r₁, #Radius 1 of ellipse
                                r₂;#Radius 2 of ellipse
                                Center = [0, 0, (Wᵢ - 1) / N * L], #The solenoid's axis is in Z
                                NPts = NPts_Coil, #How discretized the windings are
                                ),
                        [0, 0, L/2];
                        MinThreshold = 0.001,
                        )[3]
                        for Wᵢ in 1:N
                        ]

    return sum(FieldCentered)
end

ext: jl


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using ClimateMachine.Abstractions
import ClimateMachine.Abstractions: DiscontinuousSpectralElementGrid
using Impero, Printf, MPI, LinearAlgebra, Statistics,Plots
include(pwd() * "/sandbox/test_utils.jl")
ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
Ω = Circle(0,1) × Circle(0,1)

# functional 2D
n = 4
grid = DiscontinuousSpectralElementGrid(Ω, elements = (5,5), polynomialorder = (n,n), array = ArrayType)
x, y, z = coordinates(grid)
M = view(grid.vgeo, :, grid.Mid, :)    # mass matrix
xC, yC, zC = cell_centers(grid)
volumes = sum(M, dims = 1)[:]
check_x = cell_average(x, M = M)
norm(xC - check_x)
power = 10
tmp = cell_average(x.^power, M = M)
sum(tmp .* volumes) - 1/(power+1)

##
# sgeo, the first index is for the ids
# the second index are the entries
# the third index is the face
# the last index is the element
grid.sgeo[grid.sMid, :, :, :]
surface = sum(grid.sgeo[grid.sMid,:,:,:] .* x[grid.vmap⁻], dims = 1)
##
face = 1
element = 1
approx = surface[1, face, element]
# face 1 and 2 are constant x
# face 3 and 4 are constant y
# here we are integrating the function x
if face > 2
    exact = 0.5*(x[grid.vmap⁻[n+1, face, element]]^2 - x[grid.vmap⁻[1, face, element]]^2)
else
    exact = x[grid.vmap⁻[n+1, face, element]] * (y[grid.vmap⁻[n+1, face, element]] - y[grid.vmap⁻[1, face, element]])
end
println(norm(approx - exact))
gr(size=(300,300))
p1 = scatter(x[:], y[:], label = false, color = :black)
scatter!(x[grid.vmap⁻[1:n+1,face,element]], y[grid.vmap⁻[1:n+1,face,element]], color = :red, label = false)

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
#! /usr/bin/julia

# Rosetta Code, Euler's sum of powers conjecture

const lim = 250
const pwr = 5
const p = [i^pwr for i in 1:lim]

x = zeros(Int, pwr-1)
y = 0

for a in combinations(1:lim, pwr-1)
    b = searchsorted(p, sum(p[a]))
    0 < length(b) || continue
    x = a
    y = b[1]
    break
end

if y == 0
    println("No solution found for power = ", pwr, " and limit = ", lim, ".")
else
    s = [@sprintf("%d^%d", i, pwr) for i in x]
    s = join(s, " + ")
    println("A solution is ", s, " = ", @sprintf("%d^%d", y, pwr), ".")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
function wsclient(client, port::Int = 8080, url::String = "ws://localhost")
    ended = Condition()
    listen(client, :connectError) do err
        notify(ended, err)
    end
    listen(client, :connect) do ws
        listen(ws, :close) do reason
            notify(ended, reason)
        end
    end
    @async open(client, "$url:$port"; require_ssl_verification = false)
    wait(ended)
end

function clientconnects(client, port::Int = 8080, url::String = "ws://localhost")
    ended = Condition()
    listen(client, :connectError) do err
        notify(ended, err)
    end
    listen(client, :connect) do ws
        listen(ws, :invalidlistener) do noop
        end
        listen(ws, :close) do reason
            notify(ended, reason)
        end
        close(ws)
    end
    @async open(client, "$url:$port"; require_ssl_verification = false)
    wait(ended)
end

function echoclient(client, port::Int = 8080, url::String = "ws://localhost"; config...)
    config = (; config...)
    ended = Condition()
    listen(client, :connect) do ws
        listen(ws, :close) do reason
            notify(ended, reason)
        end
        listen(ws, :message) do message
            if length(message) < config.maxReceivedMessageSize
                mstring = message
                !(message isa String) && (mstring = String(copy(message)))
                if !startswith(mstring, "This") || !endswith(mstring, "test.")
                    @error "$(ws.id) : $(message isa String ? "string" : "binary") : $(length(message)) : $(Int(config.maxReceivedMessageSize))" message = message
                    throw(error("bad message format"))
                end
                if message isa String
                    message = message*message
                else
                    push!(message, message...)
                end
                if length(message) > config.maxReceivedMessageSize
                    message = message[1:Int(config.maxReceivedMessageSize)]
                end
                send(ws, message)
            elseif length(message) == config.maxReceivedMessageSize
                close(ws)
            else
                throw(error("message size mismatch"))
            end
        end
        #send(ws, "This $(ws.id) test.")
        send(ws, "This is a test.")
    end
    @async open(client, "$url:$port"; require_ssl_verification = false)
    wait(ended)
end

function badclient(client, port::Int = 8080, url::String = "ws://localhost")
    ended = Condition()
    listen(client, :connect) do ws
        listen(ws, :close) do reason
            notify(ended, reason)
        end
        listen(ws, :message) do message
            message = message*message
            send(ws, message)
        end
        send(ws, "This is a test.")
    end
    @async open(client, "$url:$port"; require_ssl_verification = false)
    wait(ended)
end

function timeoutclient(client, port::Int = 8080, url::String = "ws://localhost")
    ended = Condition()
    listen(client, :connect) do ws
        listen(ws, :close) do reason
            notify(ended, reason)
        end
        ws.keepalive[:pingmessage] = "death"
    end
    @async open(client, "$url:$port"; require_ssl_verification = false)
    wait(ended)
end

function pingclient(client, port::Int = 8080, url::String = "ws://localhost")
    ended = Condition()
    listen(client, :connect) do ws
        listen(ws, :pong) do message
            close(ws, 1000, message)
        end
        listen(ws, :close) do reason
            notify(ended, reason.description)
        end
        ping(ws, "testping")
    end
    @async open(client, "$url:$port"; require_ssl_verification = false)
    wait(ended)
end

max_stars_repo_name: citkane/SimpleWebsockets

max_issues_repo_name: citkane/SimpleWebsockets

max_forks_repo_name: citkane/SimpleWebsockets


text:
"""
    skewness(x; method=:moment)

Calculates the skewness using the specified method.

# Methods
- Moment (default)
- Fisher-Pearson
- Sample

# Arguments
- `x`:          Vector of values.
- `method`:     Estimation method: `:moment`, `:fisher_pearson` or `:sample`.
"""
function skewness(x; method::Symbol=:moment)
    n = length(x)
    mean_devs = x .- mean(x)

    if method == :moment
        # Moment
        return mean(mean_devs.^3) / sqrt(mean(mean_devs.^2))^3  # sqrt(x)^3 faster than x^1.5 !
    elseif method == :fisher_pearson
        # Fisher-Pearson
        if n > 2
            return sqrt(n*(n-1))/(n-2) * mean(mean_devs.^3) / sqrt(mean(mean_devs.^2))^3  # sqrt(x)^3 faster than x^1.5 !
        else
            return NaN
        end
    elseif method == :sample
        # Sample
        return n/((n-1)*(n-2)) * sum(mean_devs.^3 / sqrt(mean(mean_devs.^2))^3)
    end

    throw(ArgumentError("Passed method parameter '$(method)' is invalid, must be one of :moment, :fisher_pearson, :sample."))
end

"""
    kurtosis(x; method=:excess)

Calculates the kurtosis using on the specified method.

# Methods
- Excess (default)
- Moment
- Cornish-Fisher

# Arguments
- `x`:          Vector of values.
- `method`:     Estimation method: `:excess`, `:moment` or `:cornish_fisher`.
"""
function kurtosis(x; method::Symbol=:excess)
    n = length(x)
    mean_devs = x .- mean(x)

    if method == :excess
        # Excess
        return sum(mean_devs.^4 / mean(mean_devs.^2)^2 ) / n - 3
    elseif method == :moment
        # Moment
        return sum(mean_devs.^4 / mean(mean_devs.^2)^2 ) / n
    elseif method == :cornish_fisher
        # Cornish-Fisher
        return ((n+1)*(n-1)*((sum(x.^4)/n)/(sum(x.^2)/n)^2 -
            (3*(n-1))/(n+1)))/((n-2)*(n-3))
    end

    throw(ArgumentError("Passed method parameter '$(method)' is invalid, must be one of :excess, :moment, :cornish_fisher."))
end

"""
    lower_partial_moment(returns, threshold, n, method)

This function calculates the Lower Partial Moment (LPM) for a given threshold.

# Arguments
- `returns`:     Vector of asset returns.
- `threshold`:   Scalar value or vector denoting the threshold returns.
- `n`:           `n`-th moment to calculate.
- `method`:      One of `:full` or `:partial`. Indicates whether to use the number of all returns (`:full`), or only the number of returns below the threshold (`:partial`) in the denominator.
"""
function lower_partial_moment(returns, threshold, n, method::Symbol)
    if method == :full
        denominator = length(returns)
    elseif method == :partial
        denominator = count(returns .< threshold)
    else
        throw(ArgumentError("Passed method parameter '$(method)' is invalid, must be one of :full, :partial."))
    end
    excess = threshold .- returns
    sum(map(x -> max(0.0, x)^n, excess)) / denominator
end

"""
    higher_partial_moment(returns, threshold, n, method)

This function calculates the Higher Partial Moment (HPM) for a given threshold.

# Arguments
- `returns`:     Vector of asset returns.
- `threshold`:   Scalar value or vector denoting the threshold returns.
- `n`:           `n`-th moment to calculate.
- `method`:      One of `:full` or `:partial`. Indicates whether to use the number of all returns (`:full`), or only the number of returns above the threshold (`:partial`) in the denominator.
"""
function higher_partial_moment(returns, threshold, n, method::Symbol)
    if method == :full
        denominator = length(returns)
    elseif method == :partial
        denominator = count(returns .> threshold)
    else
        throw(ArgumentError("Passed method parameter '$(method)' is invalid, must be one of :full, :partial."))
    end
    excess = returns .- threshold
    sum(map(x -> max(0.0, x)^n, excess)) / denominator
end


text:
using MultiModalMuSig
using Test

K = [2, 3]
α = [0.1, 0.1]
X = Vector{Matrix{Int}}[
    # document 1
    [
        # modality 1
        [
            1 5;
            2 8
        ],
        # modality 2
        [
            1 2;
            2 5
        ]
    ],
    # document 2
    [
        # modality 1
        [
            3 4;
            4 9
        ],
        # modality 2
        [
            3 4;
            4 6
        ]
    ],
]

function calc_λ_obj(μ, invΣ, λ, ν, ζ, θ, X)
    diff = λ .- μ
    L = (
        -0.5 * (diff' * invΣ * diff)[1] +
        X[1][1][1, 2] * (θ[1][1, 1] * λ[1] + θ[1][2, 1] * λ[2]) +
        X[1][1][2, 2] * (θ[1][1, 2] * λ[1] + θ[1][2, 2] * λ[2]) +
        X[1][2][1, 2] *
            (θ[2][1, 1] * λ[3] + θ[2][2, 1] * λ[4] + θ[2][3, 1] * λ[5]) +
        X[1][2][2, 2] *
            (θ[2][1, 2] * λ[3] + θ[2][2, 2] * λ[4] + θ[2][3, 2] * λ[5]) -
        sum(X[1][1][:, 2]) / ζ[1] *
            (exp(λ[1] + 0.5ν[1]) + exp(λ[2] + 0.5ν[2])) -
        sum(X[1][2][:, 2]) / ζ[2] *
            (exp(λ[3] + 0.5ν[3]) + exp(λ[4] + 0.5ν[4]) + exp(λ[5] + 0.5ν[5]))
    )
    return L
end

function calc_λ_grad(μ, invΣ, λ, ν, ζ, θ, X)
    diff = λ .- μ
    grad = -invΣ * diff
    grad[1] += (
        X[1][1][1, 2] * θ[1][1, 1] + X[1][1][2, 2] * θ[1][1, 2]
        - sum(X[1][1][:, 2]) / ζ[1] * exp(λ[1] + 0.5ν[1])
    )
    grad[2] += (
        X[1][1][1, 2] * θ[1][2, 1] + X[1][1][2, 2] * θ[1][2, 2]
        - sum(X[1][1][:, 2]) / ζ[1] * exp(λ[2] + 0.5ν[2])
    )
    grad[3] += (
        X[1][2][1, 2] * θ[2][1, 1] + X[1][2][2, 2] * θ[2][1, 2]
        - sum(X[1][2][:, 2]) / ζ[2] * exp(λ[3] + 0.5ν[3])
    )
    grad[4] += (
        X[1][2][1, 2] * θ[2][2, 1] + X[1][2][2, 2] * θ[2][2, 2]
        - sum(X[1][2][:, 2]) / ζ[2] * exp(λ[4] + 0.5ν[4])
    )
    grad[5] += (
        X[1][2][1, 2] * θ[2][3, 1] + X[1][2][2, 2] * θ[2][3, 2]
        - sum(X[1][2][:, 2]) / ζ[2] * exp(λ[5] + 0.5ν[5])
    )
    return grad
end

@testset "λ_objective" begin
    model = MultiModalMuSig.MMCTM(K, α, X)

    μ = Float64[1, 1, 2, 2, 1]
    invΣ = model.invΣ
    λ = Float64[1, 2, 3, 4, 1]
    ν = Float64[1, 1, 1, 2, 1]
    ζ = Float64[2, 1]
    θ = Matrix{Float64}[]
    push!(θ, [0.4 0.1; 0.6 0.9])
    push!(θ, [0.3 0.4; 0.3 0.5; 0.4 0.1])

    ∇λ = Array{Float64}(undef, sum(K))
    sumθ = vcat([vec(sum(θ[m] .* X[1][m][:, 2]', dims=2)) for m in 1:model.M]...)
    Ndivζ = vcat([fill(model.N[1][m] / ζ[m], K[m]) for m in 1:model.M]...)
    objective = MultiModalMuSig.λ_objective(λ, ∇λ, ν, Ndivζ, sumθ, μ, invΣ)
    @test objective ≈ calc_λ_obj(μ, invΣ, λ, ν, ζ, θ, X)
    @test ∇λ ≈ calc_λ_grad(μ, invΣ, λ, ν, ζ, θ, X)
end


text:
function g2(byvec, valvec)
    cv = copy(valvec)
    cb = copy(byvec)
    FastGroupBy.grouptwo!(cb, cv)
    cv, cb
end

a = rand([randstring(8) for i=1:100_000], 10_000_000)
x = rand(10_000_000)

gc_enable(false)
@time g2(a,x);
gc_enable(true)

function g3(byvec, valvec)
    bv = collect(zip(Ptr{UInt}.(pointer.(byvec)), valvec))
    sort!(bv, by=x->unsafe_load(x[1]), alg = RadixSort)
end

gc_enable(false)
@time g3(a,x);
gc_enable(true)

import SortLab.load_bits

function ptrload(ptrs, skipbytes = 0)::T where T <: Unsigned
    n = sizeof(s)
    if n < skipbytes
        res = zero(T)
    elseif n - skipbytes >= sizeof(T)
        res = ntoh(unsafe_load(ptrs + skipbytes))
    else
        ptrs  = pointer(s) + skipbytes
        remaining_bytes_to_load = n - skipbytes
        res = load_bits_with_padding(T, ptrs, remaining_bytes_to_load)
    end
    return res
end

# SortingLab.load_bits(UInt, a[1], 2)
function strsort3!(svec)
    ms = maximum(sizeof, svec)
    ms = max(0, ms - sizeof(UInt))

    vs = Ptr{UInt}.(pointer.(svec))
    sort!(vs,by= x->x[1], alg=RadixSort)

    # @time  while ms > 0
    #     sort!(svec, by=x->SortingLab.load_bits(UInt, x, ms) , alg = RadixSort)
    #     ms -= max(0, sizeof(UInt))
    # end
    [vs1[2] for vs1 in vs]
end

strsort3(svec) = strsort3!(copy(svec))

gc_enable(false)
@time bb = strsort3(a)
issorted(bb)
gc_enable(true)

using SortingLab, SortingAlgorithms

srand(1);
svec = rand([randstring(8) for i=1:1_000_000], 100_000_000)
function pointersort(svec)
    # sp =
    @inbounds sv = collect(zip(1:length(svec), unsafe_load.(Ptr{UInt}.(pointer.(svec)))))
    # sv = Vector{Tuple{Ptr{UInt8}, UInt}}(length(svec))
    # for (i,s) in enumerate(svec)
    #     sp = pointer(s)
    #     sv[i] = tuple(sp, ntoh(unsafe_load(Ptr{UInt}(sp))))
    # end
    @inbounds sort!(sv, by=x->x[2], alg=RadixSort)
    @inbounds pos = [sv1[1] for sv1 in sv]
    @inbounds svec[pos]
end

gc_enable(false)
@time pointersort(svec);
@code_warntype pointersort(svec)
gc_enable(true)

gc_enable(false)
@time fsort(svec);
gc_enable(true)

function strsort5(svec)
    ps = Ptr{UInt}.(pointer.(svec))
    vs = collect(zip(ps, ntoh.(unsafe_load.(ps))))
    sort!(vs, by=x->x[2], alg=RadixSort)
    unsafe_string.(Ptr{UInt8}.(getindex.(vs,1)))
end

@time sort(svec)
using BenchmarkTools
@benchmark pointersort($svec) samples=5 seconds=120

gc_enable(false)
gc_enable(true)

gc_enable(false)
@time fsort(svec);
gc_enable(true)

@time fsortperm(svec)

using RCall
R"""
memory.limit(2^31-1)
"""
@rput svec;
rres = R"""
replicate(5, system.time(sort(svec, method="radix"))[3])
"""
mean(rres)

@time a = strsort5(svec);

@code_warntype strsort5(svec);
fsort(svec[1:1])
@time aa = fsort(svec);
all(a == aa)

x = "abc"
pointer(x)
pointer(x) - 8
pointer_from_objref(x) |> unsafe_pointer_to_objref

pointer_from_objref(svec[1]) |> unsafe_pointer_to_objref
pointer(svec)

# SortingLab.load_bits(UInt, a[1], 2)
function strsort!(svec)
    ms = maximum(sizeof, svec)
    ms = max(0, ms - sizeof(UInt))

    @time sort!(svec, by=x->SortingLab.load_bits(UInt, x, ms) , alg = RadixSort)

    @time  while ms > 0
        sort!(svec, by=x->SortingLab.load_bits(UInt, x, ms) , alg = RadixSort)
        ms -= max(0, sizeof(UInt))
    end
    svec
end

strsort(svec) = strsort!(copy(svec))

@time bb = strsort(a)
@code_warntype strsort(a)
issorted(bb)

using SortingLab
@time bb = fsort(a)
issorted(bb)

# SortingLab.load_bits(UInt, a[1], 2)
function strsort2!(svec)
    ms = maximum(sizeof, svec)
    ms = max(0, ms - sizeof(UInt))

    vs = collect(zip(SortingLab.load_bits.(UInt, svec, ms), svec))
    sort!(vs,by= x->x[1], alg=RadixSort)

    # @time  while ms > 0
    #     sort!(svec, by=x->SortingLab.load_bits(UInt, x, ms) , alg = RadixSort)
    #     ms -= max(0, sizeof(UInt))
    # end
    [vs1[2] for vs1 in vs]
end

strsort2(svec) = strsort2!(copy(svec))

gc_enable(false)
@time bb = strsort2(a)
issorted(bb)
gc_enable(true)


text:
# Example 1:
macro macint(ex)
    quote
        println("start")
        $ex
        println("after")
    end
end

@macint println("Where am I?")
#start
#Where am I?
#after

# Example 2:
macro duplicate(ex)
    quote
        $ex
        $ex
    end
end

@duplicate println(sin(42))
# -0.9165215479156338
# -0.9165215479156338

# Example 3:
macro assert(ex)
    :($ex ? nothing : error("Assertion failed: ", $(string(ex))))
end

@assert 1 == 1.0
# @assert 1 == 42 #> ERROR: Assertion failed: 1 == 42

# Example 4:
macro unless(test_cond, branch)
    quote
        if !$test_cond
            $branch
        end
    end
end

arr = [3.14, 42, 'b']
@unless 41 in arr println("arr does not contain 41")
#> array does not contain 41
@unless 42 in arr println("arr does not contain 42") #> nothing
@unless isempty(arr) println("array arr has elements")
#> array arr has elements

macroexpand(Main, :(@unless 41 in arr println("arr does not contain 41")) )
# quote
#     #= REPL[15]:3 =#
#     if !(41 in Main.arr)
#         #= REPL[15]:4 =#
#         (Main.println)("arr does not contain 41")
#     end
# end

# Example 5:
@time factorial(10)
#>   0.023171 seconds (15.36 k allocations: 797.196 KiB)
# 3628800

using Printf
macro timeit(ex)
    quote
        local t0 = time()   # local necessary to not conflict with outside
        local val = $(esc(ex))     # variables; esc makes sure ex is not
        local t1 = time()   # expanded, but is used literally
        print("elapsed time in seconds: ")
        @printf "%.3f" t1 - t0
        val
    end
end

@timeit factorial(10) #> elapsed time in seconds: 0.0003628800
@timeit factorial(15) #> elapsed time in seconds: 0.0001307674368000
a = 42
@timeit a^3 #> elapsed time in seconds: 0.00074088

# Example 6:
using Statistics
macro bench(f)
    quote
        median([@elapsed sleep($f) for x=1:10])
    end
end

@bench 0.1 #> 0.1010291465
@bench 0.2 #> 0.20101891449999998
@bench 0.5 #> 0.501034175


text:
# # Cardiac cell

# ## Model

module cardiac_cell     #jl

using ReachabilityAnalysis, ModelingToolkit

vars = @variables u, v, t

@taylorize function cardiac_cell_on!(dx, x, params, t)
    dx[1] = -0.9 * x[1]^2 - x[1]^3 - 0.9 * x[1] - x[2] + 1
    dx[2] = x[1] - 2 * x[2]
    dx[3] = one(x[3])
    return dx
end

@taylorize function cardiac_cell_off!(dx, x, params, t)
    dx[1] = -0.9 * x[1]^2 - x[1]^3 - 0.9 * x[1] - x[2]
    dx[2] = x[1] - 2 * x[2]
    dx[3] = one(x[3])
    return dx
end

function cardiac_cell_hybrid()

    n = 2 + 1 # variables + time

    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)

    ## mode 1 "approaching"
    invariant = HalfSpace(t <= 5, vars)
    mode1 = @system(x' = cardiac_cell_on!(x), dim:3, x ∈ invariant)
    ## mode 1 "approaching"
    invariant = HalfSpace(t <= 20, vars)
    mode2 = @system(x' = cardiac_cell_off!(x), dim:3, x ∈ invariant)
    modes = [mode1, mode2]

    reset = Dict(n => 0.)

    ## transition on -> off
    guard = HalfSpace(t >= 5)
    trans1 = ConstrainedResetMap(n, guard, reset)
    ## transition off -> on
    guard = HalfSpace(t >= 20)
    trans2 = ConstrainedResetMap(n, guard, reset)
    resetmaps = [trans1, trans2]

    return HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])
end

function model(X0)
    H = cardiac_cell_hybrid()
    return IVP(H, X0)
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
"""
    OSQPSolver
"""
mutable struct OSQPSolver
    options     # Dictionary of options to pass to OSQP

    model::OSQP.Model  # The OSQP solver itself

    function OSQPSolver( opts )
        solver = new()

        solver.options = opts

        return solver
    end
end

function createsolver!( solver::OSQPSolver )

    solver.model = OSQP.Model()

    options = Dict( :verbose => false )

    OSQP.setup!( solver.model; P = sparse( triu( 1.0*I(m) ) ),
                               q = zeros( Float64, (m) ),
                               A = sparse( G ),
                               u = g,
                               solver.options... )

end

function updatesolver!( solver::OSQPSolver )

end

function solve!( solver::OSQPSolver )

    results = OSQP.solve!( solver.model )

end


text:
# ## Datetime plot
# https://github.com/JuliaPlots/Plots.jl/issues/356

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Measures, Random, Dates

Random.seed!(2018)
days = 31
position = cumsum(randn(days))
x = Date(2018,1,1):Day(1):Date(2018,1,31)
ticks = [x[i] for i in 1:5:length(x)]

plot(x, position,
    xlabel="Date",
    ylabel="Position",
    title="Track of random walker",
    xticks=ticks,
    xrotation=45,
    bottom_margin=15mm,
    left_margin=5mm)


text:
#=

# Perlin noise ported from: https://github.com/caseman/noise
# with the license:

Copyright (c) 2008 Casey Duncan

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

const F2 = 0.5 * (sqrt(3.0) - 1.0)
const G2 = (3.0 - sqrt(3.0)) / 6.0

const F4 = (sqrt(5.0) - 1.0) / 4.0
const G4 = (5.0 - sqrt(5.0)) / 20.0

const SIMPLEX = [
    (0, 1, 2, 3), (0, 1, 3, 2), (0, 0, 0, 0), (0, 2, 3, 1), (0, 0, 0, 0),
    (0, 0, 0, 0), (0, 0, 0, 0), (1, 2, 3, 0), (0, 2, 1, 3), (0, 0, 0, 0),
    (0, 3, 1, 2), (0, 3, 2, 1), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0),
    (1, 3, 2, 0), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0),
    (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0), (1, 2, 0, 3),
    (0, 0, 0, 0), (1, 3, 0, 2), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0),
    (2, 3, 0, 1), (2, 3, 1, 0), (1, 0, 2, 3), (1, 0, 3, 2), (0, 0, 0, 0),
    (0, 0, 0, 0), (0, 0, 0, 0), (2, 0, 3, 1), (0, 0, 0, 0), (2, 1, 3, 0),
    (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0),
    (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0), (2, 0, 1, 3), (0, 0, 0, 0),
    (0, 0, 0, 0), (0, 0, 0, 0), (3, 0, 1, 2), (3, 0, 2, 1), (0, 0, 0, 0),
    (3, 1, 2, 0), (2, 1, 0, 3), (0, 0, 0, 0), (0, 0, 0, 0), (0, 0, 0, 0),
    (3, 1, 0, 2), (0, 0, 0, 0), (3, 2, 0, 1), (3, 2, 1, 0)
]

const GRAD4 = NTuple{4, Float32}[
    (0, 1, 1, 1), (0, 1, 1, -1), (0, 1, -1, 1), (0, 1, -1, -1),
    (0, -1, 1, 1), (0, -1, 1, -1), (0, -1, -1, 1), (0, -1, -1, -1),
    (1, 0, 1, 1), (1, 0, 1, -1), (1, 0, -1, 1), (1, 0, -1, -1),
    (-1, 0, 1, 1), (-1, 0, 1, -1), (-1, 0, -1, 1), (-1, 0, -1, -1),
    (1, 1, 0, 1), (1, 1, 0, -1), (1, -1, 0, 1), (1, -1, 0, -1),
    (-1, 1, 0, 1), (-1, 1, 0, -1), (-1, -1, 0, 1), (-1, -1, 0, -1),
    (1, 1, 1, 0), (1, 1, -1, 0), (1, -1, 1, 0), (1, -1, -1, 0),
    (-1, 1, 1, 0), (-1, 1, -1, 0), (-1, -1, 1, 0), (-1, -1, -1, 0)
]

const PERM = [
    151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140,
    36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120,
    234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
    88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71,
    134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133,
    230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161,
    1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130,
    116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250,
    124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227,
    47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44,
    154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98,
    108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34,
    242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14,
    239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121,
    50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243,
    141, 128, 195, 78, 66, 215, 61, 156, 180, 151, 160, 137, 91, 90, 15, 131,
    13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37,
    240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252,
    219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125,
    136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158,
    231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245,
    40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187,
    208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198,
    173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126,
    255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223,
    183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167,
    43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185,
    112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179,
    162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199,
    106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236,
    205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156,
    180
]

dot4(v1, x, y, z, w) = ((v1)[1]*(x) + (v1)[2]*(y) + (v1)[3]*(z) + (v1)[4]*(w))

function noise4(x, y, z, w)
    s = (x + y + z + w) * F4
    i = floor(x + s)
    j = floor(y + s)
    k = floor(z + s)
    l = floor(w + s)
    t = (i + j + k + l) * G4

    x0 = x - (i - t)
    y0 = y - (j - t)
    z0 = z - (k - t)
    w0 = w - (l - t)
    c = (x0 > y0)*32 + (x0 > z0)*16 + (y0 > z0)*8 + (x0 > w0)*4 + (y0 > w0)*2 + (z0 > w0)
    c += 1
    i1 = SIMPLEX[c][1] >= 3
    j1 = SIMPLEX[c][2] >= 3
    k1 = SIMPLEX[c][3] >= 3
    l1 = SIMPLEX[c][4] >= 3
    i2 = SIMPLEX[c][1] >= 2
    j2 = SIMPLEX[c][2] >= 2
    k2 = SIMPLEX[c][3] >= 2
    l2 = SIMPLEX[c][4] >= 2
    i3 = SIMPLEX[c][1] >= 1
    j3 = SIMPLEX[c][2] >= 1
    k3 = SIMPLEX[c][3] >= 1
    l3 = SIMPLEX[c][4] >= 1

    x1 = x0 - i1 + G4
    y1 = y0 - j1 + G4
    z1 = z0 - k1 + G4
    w1 = w0 - l1 + G4
    x2 = x0 - i2 + 2.0*G4
    y2 = y0 - j2 + 2.0*G4
    z2 = z0 - k2 + 2.0*G4
    w2 = w0 - l2 + 2.0*G4
    x3 = x0 - i3 + 3.0*G4
    y3 = y0 - j3 + 3.0*G4
    z3 = z0 - k3 + 3.0*G4
    w3 = w0 - l3 + 3.0*G4
    x4 = x0 - 1.0 + 4.0*G4
    y4 = y0 - 1.0 + 4.0*G4
    z4 = z0 - 1.0 + 4.0*G4
    w4 = w0 - 1.0 + 4.0*G4

    I = Int(i) & 255
    J = Int(j) & 255
    K = Int(k) & 255
    L = Int(l) & 255
    gi0 = PERM[1 + I + PERM[1 + J + PERM[1 + K + PERM[1 + L]]]] & 0x1f
    gi1 = PERM[1 + I + i1 + PERM[1 + J + j1 + PERM[1 + K + k1 + PERM[1 + L + l1]]]] & 0x1f
    gi2 = PERM[1 + I + i2 + PERM[1 + J + j2 + PERM[1 + K + k2 + PERM[1 + L + l2]]]] & 0x1f
    gi3 = PERM[1 + I + i3 + PERM[1 + J + j3 + PERM[1 + K + k3 + PERM[1 + L + l3]]]] & 0x1f
    gi4 = PERM[1 + I + 1 + PERM[1 + J + 1 + PERM[1 + K + 1 + PERM[1 + L + 1]]]] & 0x1f

    t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0
    noise = zeros(5)
    if t0 >= 0.0
        t0 *= t0
        noise[1] = t0 * t0 * dot4(GRAD4[gi0 + 1], x0, y0, z0, w0)
    end
    t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1
    if t1 >= 0.0
        t1 *= t1
        noise[2]= t1 * t1 * dot4(GRAD4[gi1 + 1], x1, y1, z1, w1)
    end
    t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2
    if t2 >= 0.0
        t2 *= t2
        noise[3] = t2 * t2 * dot4(GRAD4[gi2 + 1], x2, y2, z2, w2)
    end
    t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3
    if t3 >= 0.0
        t3 *= t3
        noise[4] = t3 * t3 * dot4(GRAD4[gi3 + 1], x3, y3, z3, w3)
    end
    t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4
    if t4 >= 0.0
        t4 *= t4
        noise[5] = t4 * t4 * dot4(GRAD4[gi4 + 1], x4, y4, z4, w4)
    end

    return 27.0 * sum(noise)
end

function noise4(x, y, z, w, octaves, persistence = 0.5, lacunarity = 2)
    freq = 1.0
    amp = 1.0
    max = 1.0
    total = noise4(x, y, z, w)
    for i = 1:(octaves-1)
        freq *= lacunarity
        amp *= persistence
        max += amp
        total += noise4(x * freq, y * freq, z * freq, w * freq) * amp
    end
    return total / max
end


text:
export Equality, equal

"""
Description:

    An equality constraint factor node

    f([1],[2],[3]) = δ([1] - [2]) δ([1] - [3])

Interfaces:

    1, 2, 3

Construction:

    Equality(id=:some_id)

    The interfaces of an Equality node have to be connected manually.
"""
mutable struct Equality <: DeltaFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function Equality(; id=generateId(Equality))
        self = new(id, Array{Interface}(undef, 3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:3
            self.i[idx] = self.interfaces[idx] = Interface(self)
        end

        return self
    end

    function Equality(out, in1, in2; id=generateId(Equality))
        # Explicit node constructor for creating an equality node where all connected
        # variables are considered separately
        self = new(id, Array{Interface}(undef, 3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[1] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[2] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[3] = self.interfaces[3] = associate!(Interface(self), in2)

        return self
    end
end

slug(::Type{Equality}) = "="

function equal(in1::Variable, in2::Variable)
    out = Variable()
    Equality(out, in1, in2)
    return out
end


text:
function partition(ntriplets::Int64, nthreads::Int64)

    ls = range(1, stop=ntriplets, length=nthreads+1)

    map(1:nthreads) do i
        a = round(Int64, ls[i])
        if i > 1
            a += 1
        end
        b = round(Int64, ls[i+1])
        a:b
    end
end


text:
"""

```
setbin(d,step)
```

Function that sets to which histogram bin a data point pertains
simple, but important to keep consistency over all calls.

"""
setbin(d,step) = trunc(Int, d / step ) + 1


text:
# flag generator
flag_chunksize(size)    = "-c$size"
flag_nworkers(n)        = "-j$n"
flag_latency()          = "-l"
flag_iterations(n)      = "-n$n"
flag_random()           = "-r"
flag_sequential()       = "-s"
flag_throughput()       = "-t"
flag_workingset(n)      = "-w$n"
flag_cpunuma(n)         = "-C$n"
flag_memorynuma(n)      = "-M$n"
flag_read()             = "-R"
flag_write()            = "-W"
flag_stride(n)          = "-S$n"
flag_mmapfile(str)      = "-m$str"
flag_outfile(str)       = "-f$str"

function getflags(;
        chunksize = 256,
        nworkers = 1,
        iterations = 3,
        random = false,
        sequential = false,
        workingset = 256,
        read = false,
        write = false,
        mmapfile = nothing,
        outfile = nothing,
        throughput = false,
        latency = false,
        cpunuma = 0,
        memnuma = 0,
        stride = 1,
    )

    flags = [
        flag_chunksize(chunksize),
        flag_nworkers(nworkers),
        flag_iterations(iterations),
        flag_workingset(workingset),
        flag_cpunuma(cpunuma),
        flag_memorynuma(memnuma),
        flag_stride(stride),
    ]

    read && push!(flags, flag_read())
    write && push!(flags, flag_write())
    mmapfile === nothing || push!(flags, flag_mmapfile(mmapfile))
    outfile === nothing || push!(flags, flag_outfile(outfile))
    random && push!(flags, flag_random())
    sequential && push!(flags, flag_sequential())
    throughput && push!(flags, flag_throughput())
    latency && push!(flags, flag_latency())

    return flags
end

function getoutfile(flags::Vector{String})
    prefix = flag_outfile("")
    # Search for first instance
    for flag in flags
        if startswith(flag, prefix)
            return flag[3:end]
        end
    end

    # Fallback
    outfile = "temp.csv"
    push!(flags, flag_outfile(outfile))
    return outfile
end


text:
import Base:run, |>, show, start, next, done
export AQC, state, run!, |>, clear_process!, copy_process, timer, print, reset!,success_probability

type AQC{A<:AbstractVector,T,N}<:AbstractQC{N}
    state::ComputState{A,T,N}
    status::Int # current processor index
    sys_h::AbstractMatrix # system hamiltonian
    process_list::AbstractArray
end

AQC(n::Int) = AQC{Array{Complex128,1},Complex128,n}(BellState(n), 0, spzeros(2^n,2^n), Array(AQCprocessor,0))
state(aqc::AQC) = aqc.state
hamiltonian(aqc::AQC) = aqc.sys_h

# TODO: aqc processing type for iterations
function start(aqc::AQC)
    return 1
end

function next(aqc::AQC, state)
    return aqc.process_list[state], state+1
end

function done(aqc::AQC, state)
    return state > length(aqc.process_list)
end

function run!(aqc::AQC)
    run_timer = 0.0
    for each_processor in aqc
        aqc.status += 1
        # @show each_processor
        run_timer += run!(each_processor, aqc)
    end
    return run_timer
end

function show(io::IO,aqc::AQC)
    for each_processor in aqc.process_list
        show(io::IO,each_processor)
    end
end

function (|>)(aqc::AQC,aqcprocessor::AQCprocessor)
    push!(aqc.process_list, aqcprocessor)
    return aqc
end

function clear_process!(aqc::AQC)
    aqc.process_list = Array(AQCprocessor,0)
end

function copy_process(aqc::AQC)
    return deepcopy(aqc.process_list)
end

function timer(aqc::AQC)
    ret = 0.0
    for i=1:aqc.status
        ret += timer(aqc.process_list[i])
    end
    return ret
end

function reset!{A,T,N}(aqc::AQC{A,T,N})
    aqc.state = BellState(Int(N))
end

function success_probability(aqc::AQC,tg_ground_state::AbstractVector)
    return abs2(dot(coeff(aqc.state),tg_ground_state))
end

function success_probability(aqc::AQC,pH::AbstractMatrix)
    eigen = eigfact(full(pH))
    val,idx = findmin(eigen[:values])
    return abs2(dot(coeff(aqc.state),eigen[:vectors][:,idx]))
end

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
"""
This module provides concrete implementations of `Number` that represent
1st, 2nd and general order tensors.

## Why

The main feature of this module is that the provided types do not extend from `AbstractArray`, but from `Number`!

This allows one to work with them as if they were scalar values in broadcasted operations on arrays of `VectorValue` objects (also for `TensorValue` or `MultiValue` objects). For instance, one can perform the following manipulations:
```julia
# Assing a VectorValue to all the entries of an Array of VectorValues
A = zeros(VectorValue{2,Int}, (4,5))
v = VectorValue(12,31)
A .= v # This is posible since  VectorValue <: Number

# Broatcasing of tensor operations in arrays of TensorValues
t = TensorValue(13,41,53,17) # creates a 2x2 TensorValue
g = TensorValue(32,41,3,14) # creates another 2x2 TensorValue
B = fill(t,(1,5))
C = inner.(g,B) # inner product of g against all TensorValues in the array B
@show C
# C = [2494 2494 2494 2494 2494]
```

The exported names are:

$(EXPORTS)
"""
module TensorValues

using DocStringExtensions

using StaticArrays
using Base: @propagate_inbounds, @pure

export MultiValue
export TensorValue
export VectorValue

export inner, outer, meas
#export det, inv, tr, dot, norm
export mutable
export symmetic_part
export n_components
export change_eltype

import Base: show
import Base: zero, one
import Base: +, -, *, /, \, ==, ≈
import Base: conj
import Base: sum, maximum, minimum
import Base: getindex, iterate, eachindex
import Base: size, length, eltype
import Base: reinterpret
import Base: convert
import Base: CartesianIndices
import Base: LinearIndices
import Base: adjoint
import Base: transpose

import LinearAlgebra: det, inv, tr, dot, norm

include("Types.jl")

include("Indexing.jl")

include("Operations.jl")

include("Reinterpret.jl")

end # module


text:
int_rules_1_1_1_4 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.1.4*(a+b*x)^m*(c+d*x)^n*(e+f*x)^p*(g+h*x)^q =#
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~e + ~(f') * ~x) * (~(g') + ~(h') * ~x), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) * (~g + ~h * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~m, 0) || IntegersQ(~m, ~n))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n') * (~e + ~(f') * ~x) * (~(g') + ~(h') * ~x), ~x) => (((((~b) ^ 2 * ~d * ~e * ~g - (~a) ^ 2 * ~d * ~f * ~h * ~m) - ~a * ~b * (~d * (~f * ~g + ~e * ~h) - ~c * ~f * ~h * (~m + 1))) + ~b * ~f * ~h * (~b * ~c - ~a * ~d) * (~m + 1) * ~x) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~b) ^ 2 * ~d * (~b * ~c - ~a * ~d) * (~m + 1)) + ((~a * ~d * ~f * ~h * ~m + ~b * (~d * (~f * ~g + ~e * ~h) - ~c * ~f * ~h * (~m + 2))) / ((~b) ^ 2 * ~d)) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ(~m + ~n + 2, 0) && (NeQ(~m, -1) && (SumSimplerQ(~m, 1) || Not(SumSimplerQ(~n, 1)))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~e + ~(f') * ~x) * (~(g') + ~(h') * ~x), ~x) => (((~b) ^ 2 * ~c * ~d * ~e * ~g * (~n + 1) + (~a) ^ 2 * ~c * ~d * ~f * ~h * (~n + 1) + ~a * ~b * (((~d) ^ 2 * ~e * ~g * (~m + 1) + (~c) ^ 2 * ~f * ~h * (~m + 1)) - ~c * ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 2)) + (((~a) ^ 2 * (~d) ^ 2 * ~f * ~h * (~n + 1) - ~a * ~b * (~d) ^ 2 * (~f * ~g + ~e * ~h) * (~n + 1)) + (~b) ^ 2 * (((~c) ^ 2 * ~f * ~h * (~m + 1) - ~c * ~d * (~f * ~g + ~e * ~h) * (~m + 1)) + (~d) ^ 2 * ~e * ~g * (~m + ~n + 2))) * ~x) / (~b * ~d * (~b * ~c - ~a * ~d) ^ 2 * (~m + 1) * (~n + 1))) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) - (((~a) ^ 2 * (~d) ^ 2 * ~f * ~h * (2 + 3 * ~n + (~n) ^ 2) + ~a * ~b * ~d * (~n + 1) * (2 * ~c * ~f * ~h * (~m + 1) - ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 3)) + (~b) ^ 2 * (((~c) ^ 2 * ~f * ~h * (2 + 3 * ~m + (~m) ^ 2) - ~c * ~d * (~f * ~g + ~e * ~h) * (~m + 1) * (~m + ~n + 3)) + (~d) ^ 2 * ~e * ~g * (6 + (~m) ^ 2 + 5 * ~n + (~n) ^ 2 + ~m * (2 * ~n + 5)))) / (~b * ~d * (~b * ~c - ~a * ~d) ^ 2 * (~m + 1) * (~n + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (LtQ(~m, -1) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n') * (~e + ~(f') * ~x) * (~(g') + ~(h') * ~x), ~x) => ((((((~b) ^ 3 * ~c * ~e * ~g * (~m + 2) - (~a) ^ 3 * ~d * ~f * ~h * (~n + 2)) - (~a) ^ 2 * ~b * (~c * ~f * ~h * ~m - ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 3))) - ~a * (~b) ^ 2 * (~c * (~f * ~g + ~e * ~h) + ~d * ~e * ~g * (2 * ~m + ~n + 4))) + ~b * (((~a) ^ 2 * ~d * ~f * ~h * (~m - ~n) - ~a * ~b * (2 * ~c * ~f * ~h * (~m + 1) - ~d * (~f * ~g + ~e * ~h) * (~n + 1))) + (~b) ^ 2 * (~c * (~f * ~g + ~e * ~h) * (~m + 1) - ~d * ~e * ~g * (~m + ~n + 2))) * ~x) / ((~b) ^ 2 * (~b * ~c - ~a * ~d) ^ 2 * (~m + 1) * (~m + 2))) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) + ((~f * ~h) / (~b) ^ 2 - (~d * (~m + ~n + 3) * (((~a) ^ 2 * ~d * ~f * ~h * (~m - ~n) - ~a * ~b * (2 * ~c * ~f * ~h * (~m + 1) - ~d * (~f * ~g + ~e * ~h) * (~n + 1))) + (~b) ^ 2 * (~c * (~f * ~g + ~e * ~h) * (~m + 1) - ~d * ~e * ~g * (~m + ~n + 2)))) / ((~b) ^ 2 * (~b * ~c - ~a * ~d) ^ 2 * (~m + 1) * (~m + 2))) * Antiderivative((~a + ~b * ~x) ^ (~m + 2) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (LtQ(~m, -2) || EqQ(~m + ~n + 3, 0) && Not(LtQ(~n, -2)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n') * (~e + ~(f') * ~x) * (~(g') + ~(h') * ~x), ~x) => (((~a) ^ 2 * ~d * ~f * ~h * (~n + 2) + (~b) ^ 2 * ~d * ~e * ~g * (~m + ~n + 3) + ~a * ~b * (~c * ~f * ~h * (~m + 1) - ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 3)) + ~b * ~f * ~h * (~b * ~c - ~a * ~d) * (~m + 1) * ~x) / ((~b) ^ 2 * ~d * (~b * ~c - ~a * ~d) * (~m + 1) * (~m + ~n + 3))) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) - (((~a) ^ 2 * (~d) ^ 2 * ~f * ~h * (~n + 1) * (~n + 2) + ~a * ~b * ~d * (~n + 1) * (2 * ~c * ~f * ~h * (~m + 1) - ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 3)) + (~b) ^ 2 * (((~c) ^ 2 * ~f * ~h * (~m + 1) * (~m + 2) - ~c * ~d * (~f * ~g + ~e * ~h) * (~m + 1) * (~m + ~n + 3)) + (~d) ^ 2 * ~e * ~g * (~m + ~n + 2) * (~m + ~n + 3))) / ((~b) ^ 2 * ~d * (~b * ~c - ~a * ~d) * (~m + 1) * (~m + ~n + 3))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && ((GeQ(~m, -2) && LtQ(~m, -1) || SumSimplerQ(~m, 1)) && (NeQ(~m, -1) && NeQ(~m + ~n + 3, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~e + ~(f') * ~x) * (~(g') + ~(h') * ~x), ~x) => (-((((~a * ~d * ~f * ~h * (~n + 2) + ~b * ~c * ~f * ~h * (~m + 2)) - ~b * ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 3)) - ~b * ~d * ~f * ~h * (~m + ~n + 2) * ~x)) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~b) ^ 2 * (~d) ^ 2 * (~m + ~n + 2) * (~m + ~n + 3)) + (((~a) ^ 2 * (~d) ^ 2 * ~f * ~h * (~n + 1) * (~n + 2) + ~a * ~b * ~d * (~n + 1) * (2 * ~c * ~f * ~h * (~m + 1) - ~d * (~f * ~g + ~e * ~h) * (~m + ~n + 3)) + (~b) ^ 2 * (((~c) ^ 2 * ~f * ~h * (~m + 1) * (~m + 2) - ~c * ~d * (~f * ~g + ~e * ~h) * (~m + 1) * (~m + ~n + 3)) + (~d) ^ 2 * ~e * ~g * (~m + ~n + 2) * (~m + ~n + 3))) / ((~b) ^ 2 * (~d) ^ 2 * (~m + ~n + 2) * (~m + ~n + 3))) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (NeQ(~m + ~n + 2, 0) && NeQ(~m + ~n + 3, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegersQ(~m, ~n, ~p) || IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => ((~b * ~g - ~a * ~h) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ (~p + 1)) / (~b * (~b * ~e - ~a * ~f) * (~m + 1)) - (1 / (~b * (~b * ~e - ~a * ~f) * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n - 1) * (~e + ~f * ~x) ^ ~p * Simp(~b * ~c * (~f * ~g - ~e * ~h) * (~m + 1) + (~b * ~g - ~a * ~h) * (~d * ~e * ~n + ~c * ~f * (~p + 1)) + ~d * (~b * (~f * ~g - ~e * ~h) * (~m + 1) + ~f * (~b * ~g - ~a * ~h) * (~n + ~p + 1)) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p], ~x) && (ILtQ(~m, -1) && GtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => ((~b * ~g - ~a * ~h) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ (~p + 1)) / (~b * (~b * ~e - ~a * ~f) * (~m + 1)) - (1 / (~b * (~b * ~e - ~a * ~f) * (~m + 1))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n - 1) * (~e + ~f * ~x) ^ ~p * Simp(~b * ~c * (~f * ~g - ~e * ~h) * (~m + 1) + (~b * ~g - ~a * ~h) * (~d * ~e * ~n + ~c * ~f * (~p + 1)) + ~d * (~b * (~f * ~g - ~e * ~h) * (~m + 1) + ~f * (~b * ~g - ~a * ~h) * (~n + ~p + 1)) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p], ~x) && (LtQ(~m, -1) && (GtQ(~n, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => ((~b * ~g - ~a * ~h) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * Simp((((~a * ~d * ~f * ~g - ~b * (~d * ~e + ~c * ~f) * ~g) + ~b * ~c * ~e * ~h) * (~m + 1) - (~b * ~g - ~a * ~h) * (~d * ~e * (~n + 1) + ~c * ~f * (~p + 1))) - ~d * ~f * (~b * ~g - ~a * ~h) * (~m + ~n + ~p + 3) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n, ~p], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => ((~b * ~g - ~a * ~h) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * Simp((((~a * ~d * ~f * ~g - ~b * (~d * ~e + ~c * ~f) * ~g) + ~b * ~c * ~e * ~h) * (~m + 1) - (~b * ~g - ~a * ~h) * (~d * ~e * (~n + 1) + ~c * ~f * (~p + 1))) - ~d * ~f * (~b * ~g - ~a * ~h) * (~m + ~n + ~p + 3) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n, ~p], ~x) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => (~h * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / (~d * ~f * (~m + ~n + ~p + 2)) + (1 / (~d * ~f * (~m + ~n + ~p + 2))) * Antiderivative((~a + ~b * ~x) ^ (~m - 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * Simp((~a * ~d * ~f * ~g * (~m + ~n + ~p + 2) - ~h * (~b * ~c * ~e * ~m + ~a * (~d * ~e * (~n + 1) + ~c * ~f * (~p + 1)))) + (~b * ~d * ~f * ~g * (~m + ~n + ~p + 2) + ~h * (~a * ~d * ~f * ~m - ~b * (~d * ~e * (~m + ~n + 1) + ~c * ~f * (~m + ~p + 1)))) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n, ~p], ~x) && (GtQ(~m, 0) && (NeQ(~m + ~n + ~p + 2, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => (~h * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / (~d * ~f * (~m + ~n + ~p + 2)) + (1 / (~d * ~f * (~m + ~n + ~p + 2))) * Antiderivative((~a + ~b * ~x) ^ (~m - 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * Simp((~a * ~d * ~f * ~g * (~m + ~n + ~p + 2) - ~h * (~b * ~c * ~e * ~m + ~a * (~d * ~e * (~n + 1) + ~c * ~f * (~p + 1)))) + (~b * ~d * ~f * ~g * (~m + ~n + ~p + 2) + ~h * (~a * ~d * ~f * ~m - ~b * (~d * ~e * (~m + ~n + 1) + ~c * ~f * (~m + ~p + 1)))) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n, ~p], ~x) && (GtQ(~m, 0) && (NeQ(~m + ~n + ~p + 2, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => ((~b * ~g - ~a * ~h) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * Simp((((~a * ~d * ~f * ~g - ~b * (~d * ~e + ~c * ~f) * ~g) + ~b * ~c * ~e * ~h) * (~m + 1) - (~b * ~g - ~a * ~h) * (~d * ~e * (~n + 1) + ~c * ~f * (~p + 1))) - ~d * ~f * (~b * ~g - ~a * ~h) * (~m + ~n + ~p + 3) * ~x, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n, ~p], ~x) && (ILtQ(~m + ~n + ~p + 2, 0) && (NeQ(~m, -1) && (SumSimplerQ(~m, 1) || Not(NeQ(~n, -1) && SumSimplerQ(~n, 1)) && Not(NeQ(~p, -1) && SumSimplerQ(~p, 1)))))
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(g') + ~(h') * ~x)) / (~(e') + ~(f') * ~x), ~x) => (((~f * ~g - ~e * ~h) * (~c * ~f - ~d * ~e) ^ (~m + ~n + 1)) / (~f) ^ (~m + ~n + 2)) * Antiderivative((~a + ~b * ~x) ^ ~m / ((~c + ~d * ~x) ^ (~m + 1) * (~e + ~f * ~x)), ~x) + (1 / (~f) ^ (~m + ~n + 2)) * Antiderivative(((~a + ~b * ~x) ^ ~m / (~c + ~d * ~x) ^ (~m + 1)) * ExpandToSum(((~f) ^ (~m + ~n + 2) * (~c + ~d * ~x) ^ (~m + ~n + 1) * (~g + ~h * ~x) - (~f * ~g - ~e * ~h) * (~c * ~f - ~d * ~e) ^ (~m + ~n + 1)) / (~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~m + ~n + 1, 0) && (LtQ(~m, 0) || (SumSimplerQ(~m, 1) || Not(SumSimplerQ(~n, 1)))))
    @apply_utils Antiderivative(((~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x)) / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x)), ~x) => ((~b * ~g - ~a * ~h) / (~b * ~c - ~a * ~d)) * Antiderivative((~e + ~f * ~x) ^ ~p / (~a + ~b * ~x), ~x) - ((~d * ~g - ~c * ~h) / (~b * ~c - ~a * ~d)) * Antiderivative((~e + ~f * ~x) ^ ~p / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative(((~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x)) / (~(a') + ~(b') * ~x), ~x) => (~h / ~b) * Antiderivative((~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p, ~x) + ((~b * ~g - ~a * ~h) / ~b) * Antiderivative(((~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / (sqrt(~(a') + ~(b') * ~x) * sqrt(~c + ~(d') * ~x) * sqrt(~e + ~(f') * ~x)), ~x) => (~h / ~f) * Antiderivative(sqrt(~e + ~f * ~x) / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x)), ~x) + ((~f * ~g - ~e * ~h) / ~f) * Antiderivative(1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (SimplerQ(~a + ~b * ~x, ~e + ~f * ~x) && SimplerQ(~c + ~d * ~x, ~e + ~f * ~x))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x), ~x) => (~h / ~b) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p, ~x) + ((~b * ~g - ~a * ~h) / ~b) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p], ~x) && (SumSimplerQ(~m, 1) || Not(SumSimplerQ(~n, 1)) && Not(SumSimplerQ(~p, 1)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x), ~x) => ((~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~b * (~m + 1)) - (1 / (2 * ~b * (~m + 1))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h + 2 * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h) * ~x + 3 * ~d * ~f * ~h * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegerQ(2 * ~m) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x), ~x) => (2 * (~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~b * (2 * ~m + 5)) + (1 / (~b * (2 * ~m + 5))) * Antiderivative(((~a + ~b * ~x) ^ ~m / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(((3 * ~b * ~c * ~e * ~g - ~a * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h)) + 2 * (~b * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) - ~a * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * ~x) - (3 * ~a * ~d * ~f * ~h - ~b * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegerQ(2 * ~m) && Not(LtQ(~m, -1)))
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)) / sqrt(~(c') + ~(d') * ~x), ~x) => (2 * (~a + ~b * ~x) ^ ~m * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~d * (2 * ~m + 3)) - (1 / (~d * (2 * ~m + 3))) * Antiderivative(((~a + ~b * ~x) ^ (~m - 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(((2 * ~b * ~c * ~e * ~g * ~m + ~a * (~c * (~f * ~g + ~e * ~h) - 2 * ~d * ~e * ~g * (~m + 1))) - (~b * (2 * ~d * ~e * ~g - ~c * (~f * ~g + ~e * ~h) * (2 * ~m + 1)) - ~a * (2 * ~c * ~f * ~h - ~d * (2 * ~m + 1) * (~f * ~g + ~e * ~h))) * ~x) - (2 * ~a * ~d * ~f * ~h * ~m + ~b * (~d * (~f * ~g + ~e * ~h) - 2 * ~c * ~f * ~h * (~m + 1))) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegerQ(2 * ~m) && GtQ(~m, 0))
    @apply_utils Antiderivative((sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)) / ((~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x)), ~x) => (((~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) / (~b) ^ 2) * Antiderivative(1 / ((~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + (1 / (~b) ^ 2) * Antiderivative(Simp(((~b * ~f * ~g + ~b * ~e * ~h) - ~a * ~f * ~h) + ~b * ~f * ~h * ~x, ~x) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)) / sqrt(~(c') + ~(d') * ~x), ~x) => ((~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / ((~m + 1) * (~b * ~c - ~a * ~d)) - (1 / (2 * (~m + 1) * (~b * ~c - ~a * ~d))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(~c * (~f * ~g + ~e * ~h) + ~d * ~e * ~g * (2 * ~m + 3) + 2 * (~c * ~f * ~h + ~d * (~m + 2) * (~f * ~g + ~e * ~h)) * ~x + ~d * ~f * ~h * (2 * ~m + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegerQ(2 * ~m) && LtQ(~m, -1))
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * ~x) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((2 * (~a + ~b * ~x) * sqrt(((~b * ~g - ~a * ~h) * (~c + ~d * ~x)) / ((~d * ~g - ~c * ~h) * (~a + ~b * ~x))) * sqrt(((~b * ~g - ~a * ~h) * (~e + ~f * ~x)) / ((~f * ~g - ~e * ~h) * (~a + ~b * ~x)))) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x))) * Subst(Antiderivative(1 / ((~h - ~b * (~x) ^ 2) * sqrt(1 + ((~b * ~c - ~a * ~d) * (~x) ^ 2) / (~d * ~g - ~c * ~h)) * sqrt(1 + ((~b * ~e - ~a * ~f) * (~x) ^ 2) / (~f * ~g - ~e * ~h))), ~x), ~x, sqrt(~g + ~h * ~x) / sqrt(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ (3 / 2) / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (~b / ~d) * Antiderivative((sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x)) / (sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) - ((~b * ~c - ~a * ~d) / ~d) * Antiderivative(sqrt(~a + ~b * ~x) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (2 * (~b) ^ 2 * (~a + ~b * ~x) ^ (~m - 2) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~d * ~f * ~h * (2 * ~m - 1)) - (1 / (~d * ~f * ~h * (2 * ~m - 1))) * Antiderivative(((~a + ~b * ~x) ^ (~m - 3) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((((~a * (~b) ^ 2 * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h) + 2 * (~b) ^ 3 * ~c * ~e * ~g * (~m - 2)) - (~a) ^ 3 * ~d * ~f * ~h * (2 * ~m - 1)) + ~b * ((2 * ~a * ~b * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h) + (~b) ^ 2 * (2 * ~m - 3) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h)) - 3 * (~a) ^ 2 * ~d * ~f * ~h * (2 * ~m - 1)) * ~x) - 2 * (~b) ^ 2 * (~m - 1) * (3 * ~a * ~d * ~f * ~h - ~b * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IntegerQ(2 * ~m) && GeQ(~m, 2))
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => -2 * Subst(Antiderivative(1 / (Simp((~b * ~c - ~a * ~d) - ~b * (~x) ^ 2, ~x) * sqrt(Simp((~d * ~e - ~c * ~f) / ~d + (~f * (~x) ^ 2) / ~d, ~x)) * sqrt(Simp((~d * ~g - ~c * ~h) / ~d + (~h * (~x) ^ 2) / ~d, ~x))), ~x), ~x, sqrt(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && GtQ((~d * ~e - ~c * ~f) / ~d, 0)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => -2 * Subst(Antiderivative(1 / (Simp((~b * ~c - ~a * ~d) - ~b * (~x) ^ 2, ~x) * sqrt(Simp((~d * ~e - ~c * ~f) / ~d + (~f * (~x) ^ 2) / ~d, ~x)) * sqrt(Simp((~d * ~g - ~c * ~h) / ~d + (~h * (~x) ^ 2) / ~d, ~x))), ~x), ~x, sqrt(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (Not(SimplerQ(~e + ~f * ~x, ~c + ~d * ~x)) && Not(SimplerQ(~g + ~h * ~x, ~c + ~d * ~x)))
    #= Antiderivative(1/(sqrt((~a')+(~b')*(~x))*sqrt((~c')+(~d')*(~x))*sqrt((~e')+(~f')*(~x))*sqrt((~g')+ (~h')*(~x))),~x) :=  -2*(a+b*x)*sqrt((b*g-a*h)*(c+d*x)/((d*g-c*h)*(a+b*x)))*sqrt((b*g-a*h)* (e+f*x)/((f*g-e*h)*(a+b*x)))/ ((b*g-a*h)*sqrt(c+d*x)*sqrt(e+f*x))* Subst(Antiderivative(1/(sqrt(1+(b*c-a*d)*x^2/(d*g-c*h))*sqrt(1+(b*e-a*f)*x^2/ (f*g-e*h))),x),x,sqrt(g+h*x)/sqrt(a+b*x)) <--  FreeQ([a,b,c,d,e,f,g,h],x) =#
    @apply_utils Antiderivative(1 / (sqrt(~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((2 * sqrt(~g + ~h * ~x) * sqrt(((~b * ~e - ~a * ~f) * (~c + ~d * ~x)) / ((~d * ~e - ~c * ~f) * (~a + ~b * ~x)))) / ((~f * ~g - ~e * ~h) * sqrt(~c + ~d * ~x) * sqrt((-((~b * ~e - ~a * ~f)) * (~g + ~h * ~x)) / ((~f * ~g - ~e * ~h) * (~a + ~b * ~x))))) * Subst(Antiderivative(1 / (sqrt(1 + ((~b * ~c - ~a * ~d) * (~x) ^ 2) / (~d * ~e - ~c * ~f)) * sqrt(1 - ((~b * ~g - ~a * ~h) * (~x) ^ 2) / (~f * ~g - ~e * ~h))), ~x), ~x, sqrt(~e + ~f * ~x) / sqrt(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x) ^ (3 / 2) * sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (-(~d) / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + (~b / (~b * ~c - ~a * ~d)) * Antiderivative(sqrt(~c + ~d * ~x) / ((~a + ~b * ~x) ^ (3 / 2) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m / (sqrt(~(c') + ~(d') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((~b) ^ 2 * (~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) - (1 / (2 * (~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((((2 * (~a) ^ 2 * ~d * ~f * ~h * (~m + 1) - 2 * ~a * ~b * (~m + 1) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) + (~b) ^ 2 * (2 * ~m + 3) * (~d * ~e * ~g + ~c * ~f * ~g + ~c * ~e * ~h)) - 2 * ~b * (~a * ~d * ~f * ~h * (~m + 1) - ~b * (~m + 2) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * ~x) + ~d * ~f * ~h * (2 * ~m + 5) * (~b) ^ 2 * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IntegerQ(2 * ~m) && LeQ(~m, -2))
    @apply_utils Antiderivative((sqrt(~(a') + ~(b') * ~x) * sqrt(~(c') + ~(d') * ~x)) / (sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~g + ~h * ~x)) / (~h * sqrt(~e + ~f * ~x)) + (((~d * ~e - ~c * ~f) * ((~b * ~f * ~g + ~b * ~e * ~h) - 2 * ~a * ~f * ~h)) / (2 * (~f) ^ 2 * ~h)) * Antiderivative(1 / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + ((~a * ~d * ~f * ~h - ~b * ((~d * ~f * ~g + ~d * ~e * ~h) - ~c * ~f * ~h)) / (2 * (~f) ^ 2 * ~h)) * Antiderivative(sqrt(~e + ~f * ~x) / (sqrt(~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~g + ~h * ~x)), ~x)) - (((~d * ~e - ~c * ~f) * (~f * ~g - ~e * ~h)) / (2 * ~f * ~h)) * Antiderivative(sqrt(~a + ~b * ~x) / (sqrt(~c + ~d * ~x) * (~e + ~f * ~x) ^ (3 / 2) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * sqrt(~(c') + ~(d') * ~x)) / (sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (2 * ~b * (~a + ~b * ~x) ^ (~m - 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / (~f * ~h * (2 * ~m + 1)) - (1 / (~f * ~h * (2 * ~m + 1))) * Antiderivative(((~a + ~b * ~x) ^ (~m - 2) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp((((~a * ~b * (~d * ~e * ~g + ~c * (~f * ~g + ~e * ~h)) + 2 * (~b) ^ 2 * ~c * ~e * ~g * (~m - 1)) - (~a) ^ 2 * ~c * ~f * ~h * (2 * ~m + 1)) + (((~b) ^ 2 * (2 * ~m - 1) * (~d * ~e * ~g + ~c * (~f * ~g + ~e * ~h)) - (~a) ^ 2 * ~d * ~f * ~h * (2 * ~m + 1)) + 2 * ~a * ~b * ((~d * ~f * ~g + ~d * ~e * ~h) - 2 * ~c * ~f * ~h * ~m)) * ~x) - ~b * (~a * ~d * ~f * ~h * (4 * ~m - 1) + ~b * (~c * ~f * ~h - 2 * ~d * (~f * ~g + ~e * ~h) * ~m)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegerQ(2 * ~m) && GtQ(~m, 1))
    @apply_utils Antiderivative(sqrt(~(c') + ~(d') * ~x) / ((~(a') + ~(b') * ~x) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (~d / ~b) * Antiderivative(1 / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) + ((~b * ~c - ~a * ~d) / ~b) * Antiderivative(1 / ((~a + ~b * ~x) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    #= Antiderivative(sqrt((~c')+(~d')*(~x))/(((~a')+(~b')*(~x))^(3/2)*sqrt((~e')+(~f')*(~x))*sqrt((~g')+ (~h')*(~x))),~x) :=  -2*sqrt(c+d*x)*sqrt((b*g-a*h)*(e+f*x)/((f*g-e*h)*(a+b*x)))/ ((b*g-a*h)*sqrt(e+f*x)*sqrt((b*g-a*h)*(c+d*x)/((d*g-c*h)*(a+b*x))) )* Subst(Antiderivative(sqrt(1+(b*c-a*d)*x^2/(d*g-c*h))/sqrt(1+(b*e-a*f)*x^2/(f* g-e*h)),x),x,sqrt(g+h*x)/sqrt(a+b*x)) <-- FreeQ([a,b,c,d,e,f,g,h],x) =#
    @apply_utils Antiderivative(sqrt(~(c') + ~(d') * ~x) / ((~(a') + ~(b') * ~x) ^ (3 / 2) * sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => ((-2 * sqrt(~c + ~d * ~x) * sqrt((-((~b * ~e - ~a * ~f)) * (~g + ~h * ~x)) / ((~f * ~g - ~e * ~h) * (~a + ~b * ~x)))) / ((~b * ~e - ~a * ~f) * sqrt(~g + ~h * ~x) * sqrt(((~b * ~e - ~a * ~f) * (~c + ~d * ~x)) / ((~d * ~e - ~c * ~f) * (~a + ~b * ~x))))) * Subst(Antiderivative(sqrt(1 + ((~b * ~c - ~a * ~d) * (~x) ^ 2) / (~d * ~e - ~c * ~f)) / sqrt(1 - ((~b * ~g - ~a * ~h) * (~x) ^ 2) / (~f * ~g - ~e * ~h)), ~x), ~x, sqrt(~e + ~f * ~x) / sqrt(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~m * sqrt(~(c') + ~(d') * ~x)) / (sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => (~b * (~a + ~b * ~x) ^ (~m + 1) * sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)) / ((~m + 1) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h)) + (1 / (2 * (~m + 1) * (~b * ~e - ~a * ~f) * (~b * ~g - ~a * ~h))) * Antiderivative(((~a + ~b * ~x) ^ (~m + 1) / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x))) * Simp(((2 * ~a * ~c * ~f * ~h * (~m + 1) - ~b * (~d * ~e * ~g + ~c * (2 * ~m + 3) * (~f * ~g + ~e * ~h))) + 2 * (~a * ~d * ~f * ~h * (~m + 1) - ~b * (~m + 2) * (~d * ~f * ~g + ~d * ~e * ~h + ~c * ~f * ~h)) * ~x) - ~b * ~d * ~f * ~h * (2 * ~m + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IntegerQ(2 * ~m) && LeQ(~m, -2))
    @apply_utils Antiderivative(((~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x) ^ ~q) / ((~(a') + ~(b') * ~x) * (~(c') + ~(d') * ~x)), ~x) => ((~b * ~e - ~a * ~f) / (~b * ~c - ~a * ~d)) * Antiderivative(((~e + ~f * ~x) ^ (~p - 1) * (~g + ~h * ~x) ^ ~q) / (~a + ~b * ~x), ~x) - ((~d * ~e - ~c * ~f) / (~b * ~c - ~a * ~d)) * Antiderivative(((~e + ~f * ~x) ^ (~p - 1) * (~g + ~h * ~x) ^ ~q) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~q], ~x) && LtQ(0, ~p, 1)
    @apply_utils Antiderivative(((~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~n) / (sqrt(~(e') + ~(f') * ~x) * sqrt(~(g') + ~(h') * ~x)), ~x) => Antiderivative(ExpandIntegrand(1 / (sqrt(~c + ~d * ~x) * sqrt(~e + ~f * ~x) * sqrt(~g + ~h * ~x)), (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ (~n + 1 / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IntegerQ(~m) && IntegerQ(~n + 1 / 2))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x) ^ ~q, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && IntegersQ(~p, ~q)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x) ^ ~q, ~x) => (~h / ~b) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ (~q - 1), ~x) + ((~b * ~g - ~a * ~h) / ~b) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ (~q - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p], ~x) && (IGtQ(~q, 0) && (SumSimplerQ(~m, 1) || Not(SumSimplerQ(~n, 1)) && Not(SumSimplerQ(~p, 1))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p') * (~(g') + ~(h') * ~x) ^ ~(q'), ~x) => CannotIntegrate((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p, ~q], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * ~u) ^ ~(m') * (~(c') + ~(d') * ~u) ^ ~(n') * (~(e') + ~(f') * ~u) ^ ~(p') * (~(g') + ~(h') * ~u) ^ ~(q'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p, ~q], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(i') * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~p * (~(g') + ~(h') * ~x) ^ ~q) ^ ~r, ~x) => ((~i * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * (~g + ~h * ~x) ^ ~q) ^ ~r / ((~a + ~b * ~x) ^ (~m * ~r) * (~c + ~d * ~x) ^ (~n * ~r) * (~e + ~f * ~x) ^ (~p * ~r) * (~g + ~h * ~x) ^ (~q * ~r))) * Antiderivative((~a + ~b * ~x) ^ (~m * ~r) * (~c + ~d * ~x) ^ (~n * ~r) * (~e + ~f * ~x) ^ (~p * ~r) * (~g + ~h * ~x) ^ (~q * ~r), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~m, ~n, ~p, ~q, ~r], ~x)
    @apply_utils Antiderivative((~u) ^ ~m, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m, ~x) <-- FreeQ(~m, ~x) && (LinearQ(~u, ~x) && Not(LinearMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n, ~x) <-- FreeQ([~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n') * (~w) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n * ExpandToSum(~w, ~x) ^ ~p, ~x) <-- FreeQ([~m, ~n, ~p], ~x) && (LinearQ([~u, ~v, ~w], ~x) && Not(LinearMatchQ([~u, ~v, ~w], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n') * (~w) ^ ~(p') * (~z) ^ ~(q'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n * ExpandToSum(~w, ~x) ^ ~p * ExpandToSum(~z, ~x) ^ ~q, ~x) <-- FreeQ([~m, ~n, ~p, ~q], ~x) && (LinearQ([~u, ~v, ~w, ~z], ~x) && Not(LinearMatchQ([~u, ~v, ~w, ~z], ~x)))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.4 (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl


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

cell_general = cell()
cell_design_op = cylindrical_cell_designer(cell_general)
cost = cost_default()
cell_general, cost = convert_all(cell_general, cost, mult)

using PyPlot
using ForwardDiff

using3D()

cost.general_costs.no_units_mfg = converter(20.0, mult.units_mfg)
cost.cell_costs.cathode.AM[1]
n = 35

thic = range(125, stop=170, length=n)
poro = range(0.3, stop=0.4, length=n)

normalized_thic = range(0, 1, length = n)
normalized_poro = range(0, 1, length = n)

xgrid = repeat(thic', n, 1)
ygrid = repeat(poro, 1, n)

z = zeros(n,n)

z_grad_x = zeros(n,n)
z_grad_y = zeros(n,n)

legend_data = []

for i in 1:n, j in 1:n

    cell_general.cathode.por  = poro[j]
    cell_general.cathode.th   = converter(thic[i], mult.pos_th)

    dollars_kWh, MWh_per_year = cost_calc(cell_general, cost, system="Cell", cost_verbosity=0)
    z[j, i]                   = dollars_kWh

end

clf()
fig = figure("pyplot_surfaceplot", figsize=(8,8))
plot_surface(xgrid, ygrid, z, rstride=1, edgecolors="b", cstride=1, cmap=ColorMap("gray"), alpha=0.5, linewidth=0.1)
xlabel("Cathode Thickness (um)")
ylabel("Porosity")
zlabel("dollars/kWh")
title("Surface Plot")
figure(2)


text:
using DataStructures
using Formatting
import Statistics

slurp_file(filename) = collect.(readlines(filename))

parens = Dict(
    '{' => '}',
    '<' => '>',
    '(' => ')',
    '[' => ']'
)

opening_chars = Set(keys(parens))
closing_chars = Set(values(parens))

function first_error(line)
    opened = Stack{Char}()
    for char in line
        if isempty(opened)
            if char in closing_chars
                return char
            else
                push!(opened, char)
            end
        else
            opening_char = first(opened)
            if char == parens[opening_char]
                pop!(opened)
            elseif char in closing_chars
                return char
            else
                push!(opened, char)
            end
        end
    end
    return ' '
end

scores = Dict(
    ')' => 3,
    ']' => 57,
    '}' => 1197,
    '>' => 25137,
    ' ' => 0
)

function solve1(filename)
    reduce(+, (c->scores[c]).(first_error.(slurp_file(filename))))
end

solve1("10eg.txt")
solve1("10.txt")

function get_closing(line)
    opened = Stack{Char}()
    for char in line
        if isempty(opened)
            if char in closing_chars
                return char
            else
                push!(opened, char)
            end
        else
            opening_char = first(opened)
            if char == parens[opening_char]
                pop!(opened)
            elseif char in closing_chars
                return char
            else
                push!(opened, char)
            end
        end
    end
    result = []
    for c in opened
        push!(result, parens[c])
    end
    result
end

scores2 = Dict(
    ')' => 1,
    ']' => 2,
    '}' => 3,
    '>' => 4
)

function score2(chars)
    result = 0
    for char in chars
        result *= 5
        result += scores2[char]
    end
    result
end

score2(collect("}}]])})]")) == 288957

function solve2(filename)
    lines = [line for line in slurp_file(filename) if first_error(line) == ' ']
    println(length(lines))
    scores = score2.(get_closing.(lines))
    @show scores
    sorted = sort(scores)
    Statistics.median(sorted)
end
solve2("10eg.txt") == 288957
format(solve2("10.txt"))