text:
# The code to generate Fig 4
using PyPlot
using SurNoR_2020
using Statistics

PyPlot.svg(true)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["svg.fonttype"] = "none"
rcParams["pdf.fonttype"] = 42

# ------------------------------------------------------------------------------
# Training Algorithms
# ------------------------------------------------------------------------------
Agents_SurNoR = SurNoR_2020.Func_SurNoR_Overall_train_opt_ind_i(3)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 4A: Surprise Time Series
# Only the data for subject 7 is shown in the paper
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
Epi = 1
Block = 2
for Sub = 1:12
    Agent = Agents_SurNoR[Sub,Block,Epi]

    fig = figure(figsize=(7,6)); ax = gca()
    ax.set_xlabel("Time")
    ax.set_ylabel("Surprise")
    ax.plot(Agent.Surprise[2:end])
    ax.set_ylim([0.5,1.5e5])
    ax.set_yscale("log")
    title(string("Surprise over Time, 2nd Block, 1st Epi., Sub. ",string(Agent.Sub)))
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 4B: Surprise HeatMap
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
Epi = 1
Block = 2
Surp_temp = []
# All pairs of state and actions
States = 0:43
for s = States
    push!(Surp_temp,[])
end

# Gathering the maximum surprise values for all pairs of state and actions
for Sub = 1:12
    Agent = Agents_SurNoR[Sub,Block,Epi]
    for s = States
        Ind = 1:length(Agent.State)
        if sum(Agent.State.==s)>0
            Ind = Ind[Agent.State.==s] .+ 1
            append!(Surp_temp[s+1], log10.(findmax(Agent.Surprise[Ind])[1]))
        end
    end
end

for s = States
    Surp_temp[s+1] = Float64.(Surp_temp[s+1]')
end

# Putting the values in a matrix
Surp_Matrix = zeros(10,4)
Map = [1,5,10,8,4,9,3,2,6,7] .+ 1
for s = 1:10
    for a = 1:4
        Surp_Matrix[s,a] = mean(Surp_temp[(Map[s]-1)*4 + a])
    end
end

# Plotting
fig = figure(figsize=(6,6)); ax = gca()
sc = ax.imshow(Surp_Matrix,vmin=0,vmax=5)

ax.set_yticks(0:10)
ax.set_yticklabels(["State 1", "State 2", "State 7 (S)",
                     "State 4", "State 5", "State 6", "State 3 (S)",
                     "State 8", "State 9", "State 10"],rotation=0)

ax.set_xticks(0:3)
ax.set_xticklabels(["A1", "A2", "A3", "A4"])

ax.set_ylim([9.5,-0.5])
ax.set_xlim([-0.5,3.5])

title("Heatmap Average Maximum log Surprise-Value in the 1st Episode of 2nd Block")

cbar= fig.colorbar(sc)
cbar.set_label("Surprise Value")

fig.subplots_adjust(bottom = 0.1, left = -0.3)


text: # import Mux

ext: jl


text:
using Test
using DandelionWebSockets:
    HeaderList, HTTPHandshake, performhandshake,
    issuccessful, HTTPUpgradeResponse
import DandelionWebSockets: dohandshake

mutable struct MockHTTP <: DandelionWebSockets.HTTPAdapter
    sentheaders::AbstractVector{HeaderList}
    senturis::AbstractVector{String}
    accept::String
    status::Int
    io::IO
    excess::AbstractVector{UInt8}
    ex::Union{Exception, Nothing}

    MockHTTP(;
             accept::String = "C/0nmHhBztSRGR1CwL6Tf4ZjwpY=",
             status::Int = 101,
             excess::AbstractVector{UInt8} = b"",
             io::IO = IOBuffer(),
             ex::Union{Exception,Nothing} = nothing) = new([], [], accept, status, io, excess, ex)
end

function dohandshake(m::MockHTTP, uri::String, headers::HeaderList) :: HTTPUpgradeResponse
    push!(m.sentheaders, headers)
    push!(m.senturis, uri)
    if m.ex != nothing
        throw(m.ex)
    end
    HTTPUpgradeResponse(m.io, m.status, [
        "Sec-WebSocket-Accept" => m.accept,
        "Connection" => "Upgrade",
        "Upgrade" => "websocket"],
        m.excess)
end

@testset "Handshake              " begin
    @testset "New handshake" begin
        @testset "Do a handshake; An HTTP request is sent" begin
            # Arrange
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP()
            h = HTTPHandshake(rng, mockhttp)

            # Act
            performhandshake(h, "ws://some/uri")

            # Assert
            @test length(mockhttp.sentheaders) == 1
        end

        @testset "Do a handshake; The URI is the same as provided to performhandshake" begin
            # Arrange
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP()
            h = HTTPHandshake(rng, mockhttp)

            # Act
            performhandshake(h, "ws://some/uri")

            # Assert
            @test mockhttp.senturis[1] == "ws://some/uri"
        end

        @testset "Do a handshake; Standard WebSocket headers are present in the HTTP request" begin
            # Arrange
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP()
            h = HTTPHandshake(rng, mockhttp)

            # Act
            performhandshake(h, "ws://some/uri")

            # Assert
            headers = Dict{String, String}(mockhttp.sentheaders[1])
            @test haskey(headers, "Sec-WebSocket-Key")
            @test haskey(headers, "Sec-WebSocket-Version")
            @test haskey(headers, "Upgrade")
            @test haskey(headers, "Connection")
        end

        @testset "Do a valid handshake; The handshake is successful" begin
            # Arrange
            # These fake random values generates the default accept value for MockHTTP, so this is
            # a valid handshake.
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP()
            h = HTTPHandshake(rng, mockhttp)

            # Act
            handshakeresult = performhandshake(h, "ws://some/uri")

            # Assert
            @test issuccessful(handshakeresult)
        end

        @testset "Do an invalid handshake; The handshake is unsuccessful" begin
            # Arrange
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            invalidstatuscode = 200
            mockhttp = MockHTTP(status = invalidstatuscode)
            h = HTTPHandshake(rng, mockhttp)

            # Act
            handshakeresult = performhandshake(h, "ws://some/uri")

            # Assert
            @test !issuccessful(handshakeresult)
        end

        @testset "The upgraded socket has bytes available; Those bytes can be read in the handshake result" begin
            # Arrange
            # These fake random values generates the default accept value for MockHTTP, so this is
            # a valid handshake.
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP(io = IOBuffer(b"\x01\x02\x03"))
            h = HTTPHandshake(rng, mockhttp)

            # Act
            handshakeresult = performhandshake(h, "ws://some/uri")

            # Assert
            @test readavailable(handshakeresult.io) == b"\x01\x02\x03"
        end

        @testset "Do a valid handshake with excess bytes read; The excess bytes are returned in the handshake result" begin
            # Arrange
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP(; excess = b"\xca\xfe\xba\xbe")
            h = HTTPHandshake(rng, mockhttp)

            # Act
            handshakeresult = performhandshake(h, "ws://some/uri")

            # Assert
            @test handshakeresult.excess == b"\xca\xfe\xba\xbe"
        end

        @testset "The HTTPAdapter throws an exception; The handshake is not successful" begin
            # Arrange
            rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10")
            mockhttp = MockHTTP(ex=EOFError())
            h = HTTPHandshake(rng, mockhttp)

            # Act
            handshakeresult = performhandshake(h, "ws://some/uri")

            # Assert
            @test !issuccessful(handshakeresult)
        end
    end
end


text:
using StaticArrays, DiffEqBase, OrdinaryDiffEq, RecipesBase, DiffEqCallbacks

include("./bodies.jl")
include("./boundary_conditions.jl")
include("./thermostats.jl")
include("./basic_potentials.jl")
include("./nbody_system.jl")

const kb_SI = 1.38e-23 # J/K

# This structure defines conditions under wich we test our system of n-bodies
# With this wrapping we can make such fields as `boundary_conditions` necessary for every simulation
# while allowing one to describe a particular systme of N interacting particles
struct NBodySimulation{sType <: NBodySystem,bcType <: BoundaryConditions,tType <: Real,thermType <: Thermostat}
    system::sType
    tspan::Tuple{tType,tType}
    boundary_conditions::bcType
    thermostat::thermType
    kb::tType
    external_electric_field
    external_magnetic_field
    external_gravitational_field
end

function NBodySimulation(system::BasicPotentialSystem,
    tspan::Tuple{tType,tType},
    boundary_conditions::BoundaryConditions,
    thermostat::Thermostat,
    kb::tType,
    external_electric_field,
    external_magnetic_field,
    external_gravitational_field) where {tType <: Real}

    potential_system = PotentialNBodySystem(system)
    NBodySimulation(potential_system, tspan, boundary_conditions, thermostat, kb, external_electric_field, external_magnetic_field, external_gravitational_field)
end

function NBodySimulation(system::NBodySystem, tspan::Tuple{tType,tType}, boundary_conditions::BoundaryConditions, kb::tType) where {tType <: Real}
    NBodySimulation(system, tspan, boundary_conditions, NullThermostat(), kb, x -> 0, x -> 0, x -> 0)
end

function NBodySimulation(system::NBodySystem, tspan::Tuple{tType,tType}, boundary_conditions::BoundaryConditions, thermostat::Thermostat, kb::tType) where {tType <: Real}
    NBodySimulation(system, tspan, boundary_conditions, thermostat, kb, x -> 0, x -> 0, x -> 0)
end

function NBodySimulation(system::NBodySystem, tspan::Tuple{tType,tType}, boundary_conditions::BoundaryConditions, thermostat::Thermostat) where {tType <: Real}
    NBodySimulation(system, tspan, boundary_conditions, thermostat, kb_SI, x -> 0, x -> 0, x -> 0)
end

function NBodySimulation(system::NBodySystem, tspan::Tuple{tType,tType}, boundary_conditions::BoundaryConditions) where {tType <: Real}
    NBodySimulation(system, tspan, boundary_conditions, NullThermostat(), kb_SI, x -> 0, x -> 0, x -> 0)
end

function NBodySimulation(system::NBodySystem, tspan::Tuple{tType,tType}) where {tType <: Real}
    NBodySimulation(system, tspan, InfiniteBox(), NullThermostat(), kb_SI, x -> 0, x -> 0, x -> 0)
end

function Base.show(stream::IO, s::NBodySimulation)
    print(stream, "Timespan: ")
    show(stream, s.tspan)
    println(stream)
    print(stream, "Boundary conditions: ")
    show(stream, s.boundary_conditions)
    println(stream)
    show(stream, s.system)
end

include("./nbody_to_ode.jl")
include("./nbody_simulation_result.jl")


text:
import DataFrames
import Dates
import PlotlyJS

function _chillPlot(data::DataFrames.DataFrame, temperatureField::Symbol, chillField::Symbol, model::String , units::String)::PlotlyJS.Plot

    function xaxisItem(year::Int = 1, month::Int = 1, day::Int = 1, hour::Int = 0, minute::Int = 0, second::Int = 0)::Dates.DateTime
        return Dates.DateTime(year, month, day, hour, minute, second)
    end
    dataFields::Vector{Symbol} = DataFrames.propertynames(data)
    dateFields::Vector{Symbol} = [:year, :month, :day, :hour, :minute, :second]
    usedFields::Vector{Symbol} = intersect(dateFields, dataFields)

    xaxis::Vector{Dates.DateTime} = xaxisItem.([data[!, field] for field in usedFields]...)

    xMin::Dates.DateTime = Dates.DateTime(
        Dates.year(xaxis[1]),
        Dates.month(xaxis[1]),
        Dates.day(xaxis[1])
    )

    # Start plot on a Sunday.
    while Dates.dayofweek(xMin) != Dates.Sunday
        xMin -= Dates.Day(1)
    end

    xMax::Dates.DateTime = Dates.DateTime(
        Dates.year(xaxis[length(xaxis)]),
        Dates.month(xaxis[length(xaxis)]),
        Dates.day(xaxis[length(xaxis)])
    ) + Dates.Day(1)

    # Major tick every 5.0 degrees.
    tDTick::Float64 = 2.0

    tMin::Float64 = floor(min(data[!, temperatureField]...) / tDTick) * tDTick
    tMax::Float64 = ceil(max(data[!, temperatureField]...) / tDTick) * tDTick

    ticks::Int = ceil(Int, (tMax - tMin) / tDTick)

    tMax = tMin + ticks * tDTick

    cMin::Float64 = 0.0
    cMax::Float64 = max(data[!, chillField]...)

    cDTick::Int = ceil(Int, (cMax - cMin) / ticks)
    cMax = cMin + ticks * cDTick

    trace1 = PlotlyJS.scatter(;
        name = "Temperature",
        x = xaxis,
        y = data[!, temperatureField],
        xaxis = "x",
        yaxis = "y",
        line_color = "lightgreen",
    )

    trace2 = PlotlyJS.scatter(;
        name = "Chill",
        x = xaxis,
        y = data[!, chillField],
        xaxis = "x",
        yaxis = "y2",
        line_color = "blue"
    )

    layout = PlotlyJS.Layout(;
        title = model,
        xaxis_title = "Date",
        xaxis_range = [xMin, xMax],
        xaxis_dtick = 7 * 24 * 60 * 60 * 1000,

        yaxis_title = "Temperature (in °C)",
        yaxis_color = "lightgreen",
        yaxis_range = [tMin, tMax],
        yaxis_dtick = tDTick,

        yaxis2_title = "Chill (in $(units))",
        yaxis2_color = "blue",
        yaxis2_range = [cMin, cMax],
        yaxis2_dtick = cDTick,
        yaxis2_side = "right",
        yaxis2_overlaying = "y",

        showlegend = false
    )

    return PlotlyJS.Plot([trace1, trace2], layout)
end

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
#!/usr/bin/env julia

#Monitors a GitHub repo for push events to master

using Requests
using JSON

owner = "JuliaLang"
repo = "julia"

oldt = Dates.now() - Dates.Hour(24)
#msgs = ["GitHub reporting started $(Dates.now())"]
msgs = []
while true

    t = Dates.now(Dates.UTC)

    #HTTP time format: http://www.w3.org/Protocols/rfc2616/rfc2616-sec3.html
    httptimestamp = Dates.format(oldt, Dates.RFC1123Format)#*" GMT"

    packet = try
        get(URI("https://api.github.com/repos/$owner/$repo/events"),
            headers=Dict("If-Modified-Since"=>httptimestamp))
    catch
        @goto nexttime
    end

    channel = "#bots"
    status = "good"
    if packet.status == 304 #Not modified
        info("HTTP 304: Not modified")
    elseif packet.status != 200 #Anything other than OK
        status = "#6b85dd"
        push!(msgs, "Warning: HTTP $packet.status received from GitHub API")
    else
        response = bytestring(packet.data)

        events = try
            JSON.parse(response)
        catch exc
            status = "error"
            push!(msgs, "Error in parsing GitHub response")
            filename = "events-$(Dates.now()).json"
            open(filename, "w") do f
                print(f, response)
            end
            push!(msgs, "Saved to $filename")
            []
        end

        for event in events
            if event["type"] == "PushEvent"
                user = event["actor"]["login"]
                timestamp = event["created_at"]

                if Dates.DateTime(timestamp[1:end-1], Dates.ISODateTimeFormat) < oldt
                    continue
                end

                branch = event["payload"]["ref"]
                if branch == "refs/heads/master"

                    newhash = event["payload"]["head"]
                    oldhash = event["payload"]["before"]

                    newhashurl = "https://github.com/$owner/$repo/commit/"*newhash
                    oldhashurl = "https://github.com/$owner/$repo/commit/"*newhash

                    newhashs = newhash[1:7]
                    oldhashs = oldhash[1:7]

                    dsize = event["payload"]["size"]
                    ddsize = event["payload"]["distinct_size"]

                    userurl = "https://github.com/"*user
                    branchurl = "https://github.com/$owner/$repo/tree/master"

                    #Use local git repo to compute whether this change is a fast-forward
                    if isdir(repo)
                        run(`git -C $repo fetch`)
                    else #Clone a fresh copy
                        run(`git clone https://github.com/$owner/$repo.git`)
                    end

                    #How many commits different are these?
                    ncommits = parse(Int, readall(pipeline(
                        `git -C $repo rev-list $oldhash ^$newhash`, `wc -l`)))

                    if dsize==0 || ddsize==0 || ncommits > 0
                        #master HEAD changed with no new commits. Force push?
                        status = "#d66661"
                        channel = "#general" #General alert
                    else
                        status = "good"
                    end
                    push!(msgs, "$timestamp: <$userurl|$user> pushed <$oldhashurl|$oldhashs>...<$newhashurl|$newhashs> on <$branchurl|$owner/$repo/master>: $dsize commits ($ddsize distinct; $ncommits non fast-forward)")
                end
            end
        end
    end

    msg = join(msgs, "\n")
    if msg != ""
        info("$t: Posting to Slack:\n"*msg)
        response = post(URI("https://hooks.slack.com/services/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");
            json=Dict("channel"=>channel, "username"=>"Julia-tan", "icon_emoji"=>":juliatan:",
                "attachments" => [
                    Dict("title"=>":julia: GitHub push event",
                        "color"=>status,
                        "fallback"=>msg,
                        "text"=> msg)]))
        info(bytestring(response.data))
    else
        info("$t: Nothing to report")
    end

    msgs = []
    oldt = t
    @label nexttime
    sleep(120) #Wait before polling GitHub again
end


text:
# Bogumił Kamiński, 2019-03-25

function sieve(n)
    isprime = trues(n)
    isprime[1] = false
    for i in 2:isqrt(n)
        if isprime[i]
            for j in i^2:i:n
                isprime[j] = false
            end
        end
    end
    return filter(i -> isprime[i], 1:n)
end

sieve(100)
@time length(sieve(10^6))

max_stars_repo_path: 1. Primes/primes.jl

max_issues_repo_path: 1. Primes/primes.jl

max_forks_repo_path: 1. Primes/primes.jl


text:
using FluxModels
using Test

@testset "FluxModels.jl" begin
    # Write your tests here.
end


text:
"""
    module Datasets

Commonly used datasets and utilities for creating data containers.

ToDos:

- add localization/segmentation datasets
- add labels for classification datasets

## Interface

DC<D> is a data container with observations of type D (i.e. `typeof(getobs(::DC<D>, i)::D)`)

Transformations:

- `mapobs(f::(D -> E), ::DC<D>)::DC<E>`

    Map a function (or a tuple of functions) over a data container.

- `Tuple(DC<D1>, ..., DC<DN>)::DC<(D1,...,DN)>`

    Combine multiple data containers into a single data container that returns tuples
    of the each's observations.

- `filterobs(f, DC<D>)::DC<D>`

    Keep only observations for which `f(obs) === true`.

- `groupobs(f, DC<D>)::(DC1<D>, ..., DCN<D>)` with `N` the unique return values of `f(::D)`

- `joinobs(f, DC<D1>, ..., DC<DN>)::DC<D>`

    Combines N datasets into a single one, "concatenating" them.

Primitive datasets:

- `FileDataset(dir; filterfn)`

    Each file in `dir` is one observation. Currently implemented in DLDatasets.jl with
    FileTrees.jl and observation type `FileTrees.File`.

- `TableDataset(table)`

    Every row in the table is an observation. Could use Tables.jl interface to be compatible
    with tons of packages.

## Examples

Loading and splitting an image classification dataset stored in the
same file structure as ImageNette, i.e.:

- train
    - class1
        - obs1
        - ...
        - obs2
    - class1
        - obs1
        - ...
        - obs2
- valid
    - ...

```julia
# file dataset of images `ds::DC<FileTrees.File>`
ds = FileDataset(DIR; filterfn = file -> extension(file) == "jpg")

# split into train and validation based on grandparent directory
trainds, valds = groupobs(file -> file.parent.parent.name == "train", ds)

# map (file -> input, file -> label) functions over containers to transform to type DC<(image, label)>
trainds = mapobs((FileIO.load, file -> file.parent.name), trainds)
# which is shorthand for
trainds = (
    mapobs(FileIO.load, trainds),
    mapobs(file -> file.parent.name, trainds),
)
```
---

Turning a container of (input, target) into a container of (x, y) and then an iterator
of batches (xs, ys). This is pretty much all `methoddataset` and `methoddataloaders` do.

```julia
# ds of (image, label) for example from above example
ds = ...
method = ImageClassification(...)

xyds = mapobs(ds) do (image, label)
    return encode(method, Training(), (image, label))
end

# data iterator ready to be used in training loop
dl = DataLoader(xyds, 16)

```

---

Loading an image dataset without labels for inference.

```julia
# ds contains original size images
ds = mapobs(FileIO.load(filterobs(file -> extension(file) == "jpg", FileDataset(DIR)))

transform = ProjectiveTransforms((128, 128))
ds = mapobs(ds) do obs
    run(transform, Validation(), obs)
end
# Now each observation is a center-cropped image of size (128, 128)
```
"""
module Datasets

using ..FastAI

using DataDeps
using FilePathsBase
using FilePathsBase: filename
import FileIO
using FileTrees
using MLDataPattern
using MLDataPattern: splitobs
import LearnBase
using Colors
using FixedPointNumbers

include("fastaidatasets.jl")

function __init__()
    initdatadeps()
end

include("containers.jl")
include("transformations.jl")

include("load.jl")

export
    # reexports from MLDataPattern
    splitobs,

    # container transformations
    mapobs,
    filterobs,
    groupobs,
    joinobs,

    # primitive containers
    FileDataset,
    TableDataset,
    Tokenizer,

    # utilities
    isimagefile,
    loadfile,
    filename,
    tokenize_input,

    # datasets
    DATASETS,
    loadtaskdata,
    datasetpath

end  # module


text:
"""
A Julia package to calculate standardized capital requirements
regulated by Brazilian financial system's supervisor (Brazilian Central Bank - BCB or BACEN).
"""
module BaselBrazil

	using Base.Dates
	using BusinessDays

	# Utils
	include("Utils.jl")

	# Market Risk
	## Rates
	include("RwaJur1.jl")
	include("RwaJur2.jl")
	include("RwaJur3.jl")
	include("RwaJur4.jl")
	## Currencies
	include("RwaCam.jl")

end # module


text:
@enum InfFlag::UInt8 FINITE POSINF NEGINF

"""
    Infinite <: Real

A type with two values, `∞` and `-∞` (or `PosInf` and `NegInf`).

Arithmetic with values of other types will be promoted to either the other type if it supports infinity natively or to [`InfExtended`](@ref).
"""
struct Infinite <: Real
  signbit :: Bool
end
const PosInf = Infinite(false)
const NegInf = Infinite(true)
const ∞ = PosInf
Infinite(x::Infinite) = x

"""
    InfMinusInfError()

Infinity was subtracted from infinity.
"""
struct InfMinusInfError <: Exception end

Base.showerror(io::IO, e::InfMinusInfError) = print(io, "∞-∞ is undefined")


text:
include("support.jl")

@testset CustomTestSet "Transports" begin
    @testset "TCP" begin
        @test typeof(tcp()) == Paxos.Transports.TCP.TCPTransport
        @test testRoundtrip(
            tcp(),
            ("localhost", 8000),
            "hello",
            "Ciao! I heard your greeting",
        )
        @test testBadCalls(tcp())
        @test testGoodCalls(tcp())
    end
    @testset "Memory" begin
        @test typeof(memory()) == Paxos.Transports.Memory.MemoryTransport
        @test testRoundtrip(memory(), "memoryTest1", "hello", "Ciao! I heard your greeting")
        @test testBadCalls(memory())
        @test testGoodCalls(memory())
    end
end


text:
__precompile__(true)

module SLAM

export
    # Types
    EKFSlamState,
    Scene,
    SimData,
    SlamState,
    PFSlamState,
    Vehicle,

    # Functions
    add_features,
    associate,
    compute_association,
    compute_landmark_ellipses,
    ellipse,
    local_to_global,
    get_waypoints,
    inbounds,
    initial_pose,
    laser_lines,
    mpi_to_pi,
    predict,
    predict_observation,
    steer!,
    step_vehicle!,
    update

include("common.jl")
include("ekf.jl")
include("data-association.jl")

end # SLAM


text:
#### Test Harness

using Morsel
app = Morsel.app()

route(app, GET, "/liba.tar") do req,sys
    path = dirname(@__FILE__)
    path = joinpath(path,"files")
    readall(`tar c liba`)
end

@async start(app, 4653)

#### Test script
using BinDeps

const testuri = URI("http://localhost:4653/liba.tar")

@BinDeps.setup

deps = [
    liba = library_dependency("liba", aliases = ["liba","liba1","liba.1"])
]

provides(Sources,testuri,liba,SHA="769c43644f239d8825cefc998124060cf9f477f94e8e338f6c3e17839470229d")
provides(BuildProcess,Autotools(libtarget = "liba.$(Libdl.dlext)"),liba)

@BinDeps.install Dict(:liba => :jl_liba)


text:
using SymbolicIntegration, SymbolicUtils
using Symbolics # fopr derivative
SI = SymbolicIntegration

@syms x t y Pi p q r
n =7
a = 2//3
b = -1//5
c = 5//4

#Integration Test Problems from
#https:////rulebasedintegration.org//testProblems.html

#Stewart Problems
# James Stewart - Calculus (1987)

problems = [

# Section 7.1 - Integration by Parts
[x^n,	x,	1,	x^(1+n)/(1+n)],
[exp(x),	x,	1,	exp(x)],
[1/x,	x,	1,	log(x)],
[a^x,	x,	1,	a^x/log(a)],
[sin(x),	x,	1,	-cos(x)],
[cos(x),	x,	1,	sin(x)],
[sec(x)^2,	x,	2,	tan(x)],
[csc(x)^2,	x,	2,	-cot(x)],
[sec(x)*tan(x),	x,	2,	sec(x)],
[cot(x)*csc(x),	x,	2,	-csc(x)],
[sinh(x),	x,	1,	cosh(x)],
[cosh(x),	x,	1,	sinh(x)],
[tan(x),	x,	1,	-log(cos(x))],
[cot(x),	x,	1,	log(sin(x))],
[x*sin(x),	x,	2,	-x*cos(x)+sin(x)],
[log(x),	x,	1,	-x+x*log(x)],
[exp(x)*x^2,	x,	3,	2*exp(x)-2*exp(x)*x+exp(x)*x^2],
[exp(x)*sin(x),	x,	1,	-1//2*exp(x)*cos(x)+1//2*exp(x)*sin(x)],
[atan(x),	x,	2,	x*atan(x)-1//2*log(1+x^2)],
[exp(2*x)*x,	x,	2,	-1//4*exp(2*x)+1//2*exp(2*x)*x],
[x*cos(x),	x,	2,	cos(x)+x*sin(x)],
[x*sin(4*x),	x,	2,	-1//4*x*cos(4*x)+1//16*sin(4*x)],
[x*log(x),	x,	1,	-1//4*x^2+1//2*x^2*log(x)],
[x^2*cos(3*x),	x,	3,	2//9*x*cos(3*x)-2//27*sin(3*x)+1//3*x^2*sin(3*x)],
[x^2*sin(2*x),	x,	3,	1//4*cos(2*x)-1//2*x^2*cos(2*x)+1//2*x*sin(2*x)],
[log(x)^2,	x,	2,	2*x-2*x*log(x)+x*log(x)^2],
[asin(x),	x,	2,	x*asin(x)+sqrt(1-x^2),  SI.NotImplementedError],	# asin not implemented
[t*cos(t)*sin(t),	t,	3,	-1//4*t+1//4*cos(t)*sin(t)+1//2*t*sin(t)^2],
[t*sec(t)^2,	t,	2,	log(cos(t))+t*tan(t)],
[t^2*log(t),	t,	1,	-1//9*t^3+1//3*t^3*log(t)],
[exp(t)*t^3,	t,	4,	-6*exp(t)+6*exp(t)*t-3*exp(t)*t^2+exp(t)*t^3],
[exp(2*t)*sin(3*t),	t,	1,	-3//13*exp(2*t)*cos(3*t)+2//13*exp(2*t)*sin(3*t)],
[cos(3*t)/exp(t),	t,	1,	-1//10*cos(3*t)/exp(t)+3//10*sin(3*t)/exp(t)],
[y*sinh(y),	y,	2,	y*cosh(y)-sinh(y)],
[y*cosh(a*y),	y,	2,	-cosh(a*y)/a^2+y*sinh(a*y)/a],
[t/exp(t),	t,	2,	(-1)/exp(t)-t/exp(t)],
[log(t)*sqrt(t),	t,	1,	-4//9*t^(3//2)+2//3*t^(3//2)*log(t),    SI.NotImplementedError],
[x*cos(2*x),	x,	2,	1//4*cos(2*x)+1//2*x*sin(2*x)],
[x^2/exp(x),	x,	3,	(-2)/exp(x)-2*x/exp(x)-x^2/exp(x)],
[acos(x),	x,	2,	x*acos(x)-sqrt(1-x^2), SI.NotImplementedError], #acos
[x*csc(x)^2,	x,	2,	-x*cot(x)+log(sin(x))],
[cos(5*x)*sin(3*x),	x,	1,	1//4*cos(2*x)-1//16*cos(8*x)],
[sin(2*x)*sin(4*x),	x,	1,	1//4*sin(2*x)-1//12*sin(6*x)],
[cos(x)*log(sin(x)),	x,	2,	-sin(x)+log(sin(x))*sin(x)],
[exp(x^2)*x^3,	x,	2,	-1//2*exp(x^2)+1//2*exp(x^2)*x^2],
[exp(x)*(3+2*x),	x,	2,	-2*exp(x)+exp(x)*(3+2*x)],
[5^x*x,	x,	2,	-5^x/log(5)^2+5^x*x/log(5)],
[cos(log(x)),	x,	1,	1//2*x*cos(log(x))+1//2*x*sin(log(x))],
[exp(sqrt(x)),	x,	3,	-2*exp(sqrt(x))+2*exp(sqrt(x))*sqrt(x), SI.NotImplementedError],
[log(sqrt(x)),	x,	1,	-1//2*x+x*log(sqrt(x)), SI.NotImplementedError],
[sin(log(x)),	x,	1,	-1//2*x*cos(log(x))+1//2*x*sin(log(x))],
[sin(sqrt(x)),	x,	3,	2*sin(sqrt(x))-2*cos(sqrt(x))*sqrt(x),  SI.NotImplementedError],
[x^5*cos(x^3),	x,	3,	1//3*cos(x^3)+1//3*x^3*sin(x^3)],
[exp(x^2)*x^5,	x,	3,	exp(x^2)-exp(x^2)*x^2+1//2*exp(x^2)*x^4],
[x*atan(x),	x,	3,	-1//2*x+1//2*atan(x)+1//2*x^2*atan(x)],
[x*cos(Pi*x),	x,	2,	cos(Pi*x)/Pi^2+x*sin(Pi*x)/Pi, SI.NotImplementedError],	# Pi
[log(x)*sqrt(x),	x,	1,	-4//9*x^(3//2)+2//3*x^(3//2)*log(x), SI.NotImplementedError],

# Section 7.2 - Trigonometric Integrals
[sin(3*x)^2,	x,	2,	1//2*x-1//6*cos(3*x)*sin(3*x)],
[cos(x)^2,	x,	2,	1//2*x+1//2*cos(x)*sin(x)],
[cos(x)^4,	x,	3,	3//8*x+3//8*cos(x)*sin(x)+1//4*cos(x)^3*sin(x)],
[sin(x)^3,	x,	2,	-cos(x)+1//3*cos(x)^3],
[cos(x)^4*sin(x)^3,	x,	3,	-1//5*cos(x)^5+1//7*cos(x)^7],
[cos(x)^3*sin(x)^4,	x,	3,	1//5*sin(x)^5-1//7*sin(x)^7],
[cos(x)^2*sin(x)^4,	x,	4,	1//16*x+1//16*cos(x)*sin(x)-1//8*cos(x)^3*sin(x)-1//6*cos(x)^3*sin(x)^3],
[cos(x)^2*sin(x)^2,	x,	3,	1//8*x+1//8*cos(x)*sin(x)-1//4*cos(x)^3*sin(x)],
[(1-sin(2*x))^2,	x,	1,	3//2*x+cos(2*x)-1//4*cos(2*x)*sin(2*x)],
[cos(x)*sin(1//6*Pi+x),	x,	3,	1//4*x-1//4*cos(1//6*Pi+2*x), SI.NotImplementedError],	# Pi
[cos(x)^5*sin(x)^5,	x,	4,	1//6*sin(x)^6-1//4*sin(x)^8+1//10*sin(x)^10],
[sin(x)^6,	x,	4,	5//16*x-5//16*cos(x)*sin(x)-5//24*cos(x)*sin(x)^3-1//6*cos(x)*sin(x)^5],
[cos(x)^6,	x,	4,	5//16*x+5//16*cos(x)*sin(x)+5//24*cos(x)^3*sin(x)+1//6*cos(x)^5*sin(x)],
[cos(2*x)^4*sin(2*x)^2,	x,	4,	1//16*x+1//32*cos(2*x)*sin(2*x)+1//48*cos(2*x)^3*sin(2*x)-1//12*cos(2*x)^5*sin(2*x)],
[sin(x)^5,	x,	2,	-cos(x)+2//3*cos(x)^3-1//5*cos(x)^5],
[cos(x)^4*sin(x)^4,	x,	5,	3//128*x+3//128*cos(x)*sin(x)+1//64*cos(x)^3*sin(x)-1//16*cos(x)^5*sin(x)-1//8*cos(x)^5*sin(x)^3],
[sin(x)^3*sqrt(cos(x)),	x,	3,	-2//3*cos(x)^(3//2)+2//7*cos(x)^(7//2), SI.NotImplementedError],
[cos(x)^3*sqrt(sin(x)),	x,	3,	2//3*sin(x)^(3//2)-2//7*sin(x)^(7//2), SI.NotImplementedError],
[cos(sqrt(x))^2/sqrt(x),	x,	3,	cos(sqrt(x))*sin(sqrt(x))+sqrt(x), SI.NotImplementedError],
[x*sin(x^2)^3,	x,	3,	-1//2*cos(x^2)+1//6*cos(x^2)^3],
[cos(x)^2*tan(x)^3,	x,	3,	1//2*cos(x)^2-log(cos(x))],
[cot(x)^5*sin(x)^2,	x,	4,	-1//2*csc(x)^2-2*log(sin(x))+1//2*sin(x)^2],
[(1-sin(x))/cos(x),	x,	2,	log(1+sin(x))],
[1/(1-sin(x)),	x,	1,	cos(x)/(1-sin(x))],
[tan(x)^2,	x,	2,	-x+tan(x)],
[tan(x)^4,	x,	3,	x-tan(x)+1//3*tan(x)^3],
[sec(x)^4,	x,	2,	tan(x)+1//3*tan(x)^3],
[sec(x)^6,	x,	2,	tan(x)+2//3*tan(x)^3+1//5*tan(x)^5],
[sec(x)^2*tan(x)^4,	x,	2,	1//5*tan(x)^5],
[sec(x)^4*tan(x)^2,	x,	3,	1//3*tan(x)^3+1//5*tan(x)^5],
[sec(x)^3*tan(x),	x,	2,	1//3*sec(x)^3],
[sec(x)^3*tan(x)^3,	x,	3,	-1//3*sec(x)^3+1//5*sec(x)^5],
[tan(x)^5,	x,	3,	-log(cos(x))-1//2*tan(x)^2+1//4*tan(x)^4],
[tan(x)^6,	x,	4,	-x+tan(x)-1//3*tan(x)^3+1//5*tan(x)^5],
[sec(x)*tan(x)^5,	x,	3,	sec(x)-2//3*sec(x)^3+1//5*sec(x)^5],
[sec(x)^3*tan(x)^5,	x,	3,	1//3*sec(x)^3-2//5*sec(x)^5+1//7*sec(x)^7],
[sec(x)^6*tan(x),	x,	2,	1//6*sec(x)^6],
[sec(x)^6*tan(x)^3,	x,	3,	-1//6*sec(x)^6+1//8*sec(x)^8],
[sec(x)^2/cot(x),	x,	2,	1//2*sec(x)^2],
[sec(x)*tan(x)^2,	x,	2,	-1//2*atanh(sin(x))+1//2*sec(x)*tan(x)],
[cot(x)^2,	x,	2,	-x-cot(x)],
[cot(x)^3,	x,	2,	-1//2*cot(x)^2-log(sin(x))],
[cot(x)^4*csc(x)^4,	x,	3,	-1//5*cot(x)^5-1//7*cot(x)^7],
[cot(x)^3*csc(x)^4,	x,	3,	1//4*csc(x)^4-1//6*csc(x)^6],
[csc(x),	x,	1,	-atanh(cos(x))],
[csc(x)^3,	x,	2,	-1//2*atanh(cos(x))-1//2*cot(x)*csc(x)],
[cos(x)^2/sin(x),	x,	3,	-atanh(cos(x))+cos(x)],
[1/sin(x)^4,	x,	2,	-cot(x)-1//3*cot(x)^3],
[sin(2*x)*sin(5*x),	x,	1,	1//6*sin(3*x)-1//14*sin(7*x)],
[cos(x)*sin(3*x),	x,	1,	-1//4*cos(2*x)-1//8*cos(4*x)],
[cos(3*x)*cos(4*x),	x,	1,	1//2*sin(x)+1//14*sin(7*x)],
[sin(3*x)*sin(6*x),	x,	1,	1//6*sin(3*x)-1//18*sin(9*x)],
[cos(x)^5*sin(x),	x,	2,	-1//6*cos(x)^6],
[cos(x)*cos(2*x)*cos(3*x),	x,	5,	1//4*x+1//8*sin(2*x)+1//16*sin(4*x)+1//24*sin(6*x)],
[(1-tan(x)^2)/sec(x)^2,	x,	2,	cos(x)*sin(x)],
[(cos(x)+sin(x))/sin(2*x),	x,	6,	-1//2*atanh(cos(x))+1//2*atanh(sin(x))],
[sin(x)^2*tan(x),	x,	3,	1//2*cos(x)^2-log(cos(x))],
[cos(x)^2*cot(x)^3,	x,	4,	-1//2*csc(x)^2-2*log(sin(x))+1//2*sin(x)^2],
[sec(x)^3*tan(x),	x,	2,	1//3*sec(x)^3],
[sec(x)^3*tan(x)^3,	x,	3,	-1//3*sec(x)^3+1//5*sec(x)^5],

# Section 7.3 - Trigonometric Substitution
[sqrt(9-x^2)/x^2,	x,	2,	-asin(1//3*x)-sqrt(9-x^2)/x, SI.NotImplementedError],
[1/(x^2*sqrt(4+x^2)),	x,	1,	-1//4*sqrt(4+x^2)/x, SI.NotImplementedError],
[x/sqrt(4+x^2),	x,	1,	sqrt(4+x^2), SI.NotImplementedError],
[1/sqrt(-a^2+x^2),	x,	2,	atanh(x/sqrt(-a^2+x^2)), SI.NotImplementedError],	# sqrt, a
[x^3/(9+4*x^2)^(3//2),	x,	3,	9//16/sqrt(9+4*x^2)+1//16*sqrt(9+4*x^2), SI.NotImplementedError],	# exponent 3//2
[x/sqrt(3-2*x-x^2),	x,	3,	asin(1//2*(-1-x))-sqrt(3-2*x-x^2), SI.NotImplementedError],
[1/(x^2*sqrt(1-x^2)),	x,	1,	-sqrt(1-x^2)/x, SI.NotImplementedError],
[x^3*sqrt(4-x^2),	x,	3,	-4//3*(4-x^2)^(3//2)+1//5*(4-x^2)^(5//2), SI.NotImplementedError],
[x/sqrt(1-x^2),	x,	1,	-sqrt(1-x^2), SI.NotImplementedError],
[x*sqrt(4-x^2),	x,	1,	-1//3*(4-x^2)^(3//2), SI.NotImplementedError],
[sqrt(1-4*x^2),	x,	2,	1//4*asin(2*x)+1//2*x*sqrt(1-4*x^2), SI.NotImplementedError],
[x^3/sqrt(4+x^2),	x,	3,	1//3*(4+x^2)^(3//2)-4*sqrt(4+x^2), SI.NotImplementedError],
[1/sqrt(9+x^2),	x,	1,	asinh(1//3*x), SI.NotImplementedError],
[sqrt(1+x^2),	x,	2,	1//2*asinh(x)+1//2*x*sqrt(1+x^2), SI.NotImplementedError],
[1/(x^3*sqrt(-16+x^2)),	x,	4,	1//128*atan(1//4*sqrt(-16+x^2))+1//32*sqrt(-16+x^2)/x^2, SI.NotImplementedError],
[sqrt(-a^2+x^2)/x^4,	x,	1,	1//3*(-a^2+x^2)^(3//2)/(a^2*x^3), SI.NotImplementedError],	# sqrt, a
[sqrt(-4+9*x^2)/x,	x,	4,	-2*atan(1//2*sqrt(-4+9*x^2))+sqrt(-4+9*x^2), SI.NotImplementedError],
[1/(x^2*sqrt(-9+16*x^2)),	x,	1,	1//9*sqrt(-9+16*x^2)/x, SI.NotImplementedError],
[x^2/(a^2-x^2)^(3//2),	x,	3,	-atan(x/sqrt(a^2-x^2))+x/sqrt(a^2-x^2), SI.NotImplementedError],  # exponent 3//2, a
[x^2/sqrt(5-x^2),	x,	2,	5//2*asin(x/sqrt(5))-1//2*x*sqrt(5-x^2), SI.NotImplementedError],
[1/(x*sqrt(3+x^2)),	x,	3,	-atanh(sqrt(3+x^2)/sqrt(3))/sqrt(3), SI.NotImplementedError],
[x/(4+x^2)^(5//2),	x,	1,	(-1//3)//(4+x^2)^(3//2), SI.NotImplementedError],	# exponent 5//2
[x^3*sqrt(4-9*x^2),	x,	3,	-4//243*(4-9*x^2)^(3//2)+1//405*(4-9*x^2)^(5//2), SI.NotImplementedError],
[x^2*sqrt(9-x^2),	x,	3,	81//8*asin(1//3*x)-9//8*x*sqrt(9-x^2)+1//4*x^3*sqrt(9-x^2), SI.NotImplementedError],
[5*x*sqrt(1+x^2),	x,	2,	5//3*(1+x^2)^(3//2), SI.NotImplementedError],
[1/(-25+4*x^2)^(3//2),	x,	1,	-1//25*x/sqrt(-25+4*x^2), SI.NotImplementedError],
[sqrt(2*x-x^2),	x,	3,	-1//2*asin(1-x)-1//2*(1-x)*sqrt(2*x-x^2), SI.NotImplementedError],
[1/sqrt(8+4*x+x^2),	x,	2,	asinh(1//2*(2+x)), SI.NotImplementedError],
[1/sqrt(-8+6*x+9*x^2),	x,	2,	1//3*atanh((1+3*x)/sqrt(-8+6*x+9*x^2)), SI.NotImplementedError],
[x^2/sqrt(4*x-x^2),	x,	4,	-6*asin(1-1//2*x)-3*sqrt(4*x-x^2)-1//2*x*sqrt(4*x-x^2), SI.NotImplementedError],
[1/(2+2*x+x^2)^2,	x,	3,	1//2*(1+x)/(2+2*x+x^2)+1//2*atan(1+x)],
[1/(5-4*x-x^2)^(5//2),	x,	2,	1//27*(2+x)/(5-4*x-x^2)^(3//2)+2//243*(2+x)/sqrt(5-4*x-x^2), SI.NotImplementedError], #exponent 5//2
[exp(t)*sqrt(9-exp(2*t)),	t,	3,	9//2*asin(1//3*exp(t))+1//2*exp(t)*sqrt(9-exp(2*t)), SI.NotImplementedError],
[sqrt(-9+exp(2*t)),	t,	4,	-3*atan(1//3*sqrt(-9+exp(2*t)))+sqrt(-9+exp(2*t)), SI.NotImplementedError],

# Section 7.4 - Integration of Rational Functions by Partial Fractions
[1//sqrt(a^2+x^2),	x,	2,	atanh(x//sqrt(a^2+x^2)), SI.NotImplementedError],	# sqrt, a
[(5+x)//(-2+x+x^2),	x,	3,	2*log(1-x)-log(2+x)],
[(x+x^3)//(-1+x),	x,	3,	2*x+1//2*x^2+1//3*x^3+2*log(1-x)],
[(-1+2*x+x^2)//(-2*x+3*x^2+2*x^3),	x,	3,	1//10*log(1-2*x)+1//2*log(x)-1//10*log(2+x)],
[(1+4*x-2*x^2+x^4)//(1-x-x^2+x^3),	x,	2,	2//(1-x)+x+1//2*x^2+log(1-x)-log(1+x)],
[(4-x+2*x^2)//(4*x+x^3),	x,	6,	-1//2*atan(1//2*x)+log(x)+1//2*log(4+x^2)],
[(2-3*x+4*x^2)//(3-4*x+4*x^2),	x,	6,	x+1//8*log(3-4*x+4*x^2)+1//4*atan((1-2*x)/sqrt(2))/sqrt(2)],
[(1+x^2+x^3)//((-1+x)*x*(1+x^2)^3*(1+x+x^2)),	x,	14,	1//8*(1+x)//(1+x^2)^2-3//8*(1-x)//(1+x^2)+3//16*x//(1+x^2)+7//16*atan(x)+1//8*log(1-x)-log(x)+15//16*log(1+x^2)-1//2*log(1+x+x^2)-atan((1+2*x)/sqrt(3))/sqrt(3)],
[(1-3*x+2*x^2-x^3)//(x*(1+x^2)^2),	x,	6,	1//2*(-1-2*x)//(1+x^2)-2*atan(x)+log(x)-1//2*log(1+x^2)],
[1//(1+x^2)^2,	x,	2,	1//2*x//(1+x^2)+1//2*atan(x)],
[1//((-1+x)*(2+x)),	x,	3,	1//3*log(1-x)-1//3*log(2+x)],
[7//(-12+5*x+2*x^2),	x,	4,	7//11*log(3-2*x)-7//11*log(4+x)],
[(-4+3*x+x^2)//((-1+2*x)^2*(3+2*x)),	x,	2,	(-9//32)//(1-2*x)+41//128*log(1-2*x)-25//128*log(3+2*x)],
[(-x^2+x^3)//((-6+x)*(3+5*x)^3),	x,	3,	(-12//1375)//(3+5*x)^2+201//15125//(3+5*x)+20//3993*log(6-x)+1493//499125*log(3+5*x)],
[1//(-x^3+x^4),	x,	3,	1//2//x^2+1//x+log(1-x)-log(x)],
[(1-x-x^2+x^3+x^4)//(-x+x^3),	x,	4,	x+1//2*x^2-log(x)+1//2*log(1-x^2)],
[(-2+x^2)//(x*(2+x^2)),	x,	3,	-log(x)+log(2+x^2)],
[(2-4*x^2+x^3)//((1+x^2)*(2+x^2)),	x,	8,	6*atan(x)-1//2*log(1+x^2)+log(2+x^2)-5*atan(x/sqrt(2))*sqrt(2)],
[(1+x^2+x^4)//((1+x^2)*(4+x^2)^2),	x,	6,	-13//24*x//(4+x^2)+25//144*atan(1//2*x)+1//9*atan(x)],
[(1+16*x)//((5+x)^2*(-3+2*x)*(1+x+x^2)),	x,	6,	(-79//273)//(5+x)+200//3211*log(3-2*x)+2731//24843*log(5+x)-481//5586*log(1+x+x^2)+451//2793*atan((1+2*x)/sqrt(3))/sqrt(3)],
[x^4//(9+x^2)^3,	x,	3,	-1//4*x^3//(9+x^2)^2-3//8*x//(9+x^2)+1//8*atan(1//3*x)],
[19*x//((-1+x)^3*(3+5*x+4*x^2)^2),	x,	8,	(-399//736)//(1-x)^2+(-1843//4416)//(1-x)+19//276*(39+44*x)//((1-x)^2*(3+5*x+4*x^2))+209//2304*log(1-x)-209//4608*log(3+5*x+4*x^2)+114437//52992*atan((5+8*x)/sqrt(23))/sqrt(23)],
[(1+x^2+x^3)//(2*x^2+x^3+x^4),	x,	7,	(-1//2)//x-1//4*log(x)+5//8*log(2+x+x^2)+1//4*atan((1+2*x)/sqrt(7))/sqrt(7)],
[1//(-x^3+x^6),	x,	8,	1//2//x^2+1//3*log(1-x)-1//6*log(1+x+x^2)-atan((1+2*x)/sqrt(3))/sqrt(3)],
[x^2//(1+x),	x,	2,	-x+1//2*x^2+log(1+x)],
[x//(-5+x),	x,	2,	x+5*log(5-x)],
[(-1+4*x)//((-1+x)*(2+x)),	x,	2,	log(1-x)+3*log(2+x)],
[1//((1+x)*(2+x)),	x,	3,	log(1+x)-log(2+x)],
[(-5+6*x)//(3+2*x),	x,	2,	3*x-7*log(3+2*x)],
[1//((a+x)*(b+x)),	x,	3,	-log(a+x)/(a-b)+log(b+x)/(a-b)],
[(1+x^2)//(-x+x^2),	x,	3,	x+2*log(1-x)-log(x)],
[(1-12*x+x^2+x^3)//(-12+x+x^2),	x,	5,	1//2*x^2+1//7*log(3-x)-1//7*log(4+x)],
[(3+2*x)//(1+x)^2,	x,	2,	(-1)//(1+x)+2*log(1+x)],
[1//(x*(1+x)*(3+2*x)),	x,	2,	1//3*log(x)-log(1+x)+2//3*log(3+2*x)],
[(-3+5*x+6*x^2)//(-3*x+2*x^2+x^3),	x,	3,	2*log(1-x)+log(x)+3*log(3+x)],
[x//(4+4*x+x^2),	x,	3,	2//(2+x)+log(2+x)],
[1//((-1+x)^2*(4+x)),	x,	2,	1//5//(1-x)-1//25*log(1-x)+1//25*log(4+x)],
[x^2//((-3+x)*(2+x)^2),	x,	2,	4//5//(2+x)+9//25*log(3-x)+16//25*log(2+x)],
[(-2+3*x+5*x^2)//(2*x^2+x^3),	x,	3,	1//x+2*log(x)+3*log(2+x)],
[(18-2*x-4*x^2)//(-6+x+4*x^2+x^3),	x,	2,	log(1-x)-2*log(2+x)-3*log(3+x)],
[(2*x+x^2)//(4+3*x^2+x^3),	x,	1,	1//3*log(4+3*x^2+x^3)],
[1//((-1+x)^2*x^2),	x,	2,	1//(1-x)+(-1)//x-2*log(1-x)+2*log(x)],
[x^2//(1+x)^3,	x,	2,	(-1//2)//(1+x)^2+2//(1+x)+log(1+x)],
[1//(-x^2+x^4),	x,	3,	1//x-atanh(x)],
[(-x+2*x^3)//(1-x^2+x^4),	x,	1,	1//2*log(1-x^2+x^4)],
[x^3//(1+x^2),	x,	3,	1//2*x^2-1//2*log(1+x^2)],
[(-1+x)//(2+2*x+x^2),	x,	4,	-2*atan(1+x)+1//2*log(2+2*x+x^2)],
[x//(1+x+x^2),	x,	4,	1//2*log(1+x+x^2)-atan((1+2*x)/sqrt(3))/sqrt(3)],
[(7+5*x+4*x^2)//(5+4*x+4*x^2),	x,	6,	x+3//8*atan(1//2+x)+1//8*log(5+4*x+4*x^2)],
[(5-4*x+3*x^2)//((-1+x)*(1+x^2)),	x,	5,	-3*atan(x)+2*log(1-x)+1//2*log(1+x^2)],
[(3+2*x)//(3*x+x^3),	x,	6,	log(x)-1//2*log(3+x^2)+2*atan(x/sqrt(3))/sqrt(3)],
[1//(-1+x^3),	x,	6,	1//3*log(1-x)-1//6*log(1+x+x^2)-atan((1+2*x)/sqrt(3))/sqrt(3)],
[x^3//(1+x^3),	x,	7,	x-1//3*log(1+x)+1//6*log(1-x+x^2)+atan((1-2*x)/sqrt(3))/sqrt(3)],
[(-1-2*x+x^2)//((-1+x)^2*(1+x^2)),	x,	5,	1//(-1+x)+atan(x)+log(1-x)-1//2*log(1+x^2)],
[x^4//(-1+x^4),	x,	4,	x-1//2*atan(x)-1//2*atanh(x)],
[(-4+6*x-x^2+3*x^3)//((1+x^2)*(2+x^2)),	x,	6,	-3*atan(x)+3//2*log(1+x^2)+atan(x/sqrt(2))*sqrt(2)],
[(1+x-2*x^2+x^3)//(4+5*x^2+x^4),	x,	7,	-3//2*atan(1//2*x)+atan(x)+1//2*log(4+x^2)],
[(-3+x)//(4+2*x+x^2)^2,	x,	3,	1//6*(-7-4*x)//(4+2*x+x^2)-2//3*atan((1+x)/sqrt(3))/sqrt(3)],
[(1+x^4)//(x*(1+x^2)^2),	x,	3,	1//(1+x^2)+log(x)],
[cos(x)*(-3+2*sin(x))//(2-3*sin(x)+sin(x)^2),	x,	2,	log(2-3*sin(x)+sin(x)^2)],
[cos(x)^2*sin(x)//(5+cos(x)^2),	x,	3,	-cos(x)+atan(cos(x)/sqrt(5))*sqrt(5)],  # result contains algebraic numbers
[1//(-3+2*x+x^2),	x,	3,	1//4*log(1-x)-1//4*log(3+x)],
[1//(-2*x+x^2),	x,	1,	1//2*log(2-x)-1//2*log(x)],
[(1+2*x)//(-7+12*x+4*x^2),	x,	3,	1//8*log(1-2*x)+3//8*log(7+2*x)],
[x//(-1+x+x^2),	x,	3,	1//10*log(1+2*x-sqrt(5))*(5-sqrt(5))+1//10*log(1+2*x+sqrt(5))*(5+sqrt(5))],

# Section 7.5 - Rationalization Substitutions
[(-32+5*x-27*x^2+4*x^3)//(-70-299*x-286*x^2+50*x^3-13*x^4+30*x^5),	x,	6,	-3146//80155*log(7-3*x)-334//323*log(1+2*x)+4822//4879*log(2+5*x)+11049//260015*log(5+x+x^2)+3988//13685*atan((1+2*x)/sqrt(19))/sqrt(19)],
[(8-13*x^2-7*x^3+12*x^5)//(4-20*x+41*x^2-80*x^3+116*x^4-80*x^5+100*x^6),	x,	7,	5828//9075//(2-5*x)+1//1452*(-313-502*x)//(1+2*x^2)-59096//99825*log(2-5*x)+2843//7986*log(1+2*x^2)-251//726*atan(x*sqrt(2))/sqrt(2)+272//1331*atan(x*sqrt(2))*sqrt(2)],
[sqrt(4+x)//x,	x,	3,	-4*atanh(1//2*sqrt(4+x))+2*sqrt(4+x), SI.NotImplementedError],
[1//((-1)//x^(1//3)+sqrt(x)),	x,	9,	6//5*log(1-x^(1//6))-3//10*log(2+x^(1//6)+2*x^(1//3)+x^(1//6)*sqrt(5))*(1-sqrt(5))-3//10*log(2+x^(1//6)+2*x^(1//3)-x^(1//6)*sqrt(5))*(1+sqrt(5))+2*sqrt(x)+3//5*atan((1+4*x^(1//6)-sqrt(5))/sqrt(2*(5+sqrt(5))))*sqrt(2*(5-sqrt(5)))-3//5*atan(1//2*(1+4*x^(1//6)+sqrt(5))*sqrt(1//10*(5+sqrt(5))))*sqrt(2*(5+sqrt(5))), SI.NotImplementedError],
[1//(-4*cos(x)+3*sin(x)),	x,	2,	-1//5*atanh(1//5*(3*cos(x)+4*sin(x)))],
[1//(1+sqrt(x)),	x,	3,	-2*log(1+sqrt(x))+2*sqrt(x), SI.NotImplementedError],
[1//(1+1//x^(1//3)),	x,	3,	3*x^(1//3)-3//2*x^(2//3)+x-3*log(1+1//x^(1//3))-log(x), SI.NotImplementedError], #exponent 1//3
[sqrt(x)//(1+x),	x,	3,	-2*atan(sqrt(x))+2*sqrt(x), SI.NotImplementedError],
[1//(x*sqrt(1+x)),	x,	2,	-2*atanh(sqrt(1+x)), SI.NotImplementedError],
[1//(-x^(1//3)+x),	x,	2,	3//2*log(1-x^(2//3)), SI.NotImplementedError],	# exponent 1//3
[1//(x-sqrt(2+x)),	x,	4,	4//3*log(2-sqrt(2+x))+2//3*log(1+sqrt(2+x)), SI.NotImplementedError],
[x^2//sqrt(-1+x),	x,	2,	4//3*(-1+x)^(3//2)+2//5*(-1+x)^(5//2)+2*sqrt(-1+x), SI.NotImplementedError],
[sqrt(-1+x)//(1+x),	x,	3,	-2*atan(sqrt(-1+x)/sqrt(2))*sqrt(2)+2*sqrt(-1+x), SI.NotImplementedError],
[1//sqrt(1+sqrt(x)),	x,	3,	4//3*(1+sqrt(x))^(3//2)-4*sqrt(1+sqrt(x)), SI.NotImplementedError],
[sqrt(x)//(x+x^2),	x,	3,	2*atan(sqrt(x)), SI.NotImplementedError],
[(1+sqrt(x))//(-1+sqrt(x)),	x,	3,	x+4*log(1-sqrt(x))+4*sqrt(x), SI.NotImplementedError],
[(1+1//x^(1//3))//(-1+1//x^(1//3)),	x,	4,	-6*x^(1//3)-3*x^(2//3)-x-6*log(1-x^(1//3)), SI.NotImplementedError],	# exponent 1//3
[x^3//(1+x^2)^(1//3),	x,	3,	-3//4*(1+x^2)^(2//3)+3//10*(1+x^2)^(5//3), SI.NotImplementedError],	#exponent 1//3
[sqrt(x)//((-1)//x^(1//3)+sqrt(x)),	x,	10,	6*x^(1//6)+x+6//5*log(1-x^(1//6))-3//10*log(2+x^(1//6)+2*x^(1//3)-x^(1//6)*sqrt(5))*(1-sqrt(5))-3//10*log(2+x^(1//6)+2*x^(1//3)+x^(1//6)*sqrt(5))*(1+sqrt(5))-3//5*atan(1//2*(1+4*x^(1//6)+sqrt(5))*sqrt(1//10*(5+sqrt(5))))*sqrt(2*(5-sqrt(5)))-3//5*atan((1+4*x^(1//6)-sqrt(5))/sqrt(2*(5+sqrt(5))))*sqrt(2*(5+sqrt(5))), SI.NotImplementedError],
[1//(1//x^(1//4)+sqrt(x)),	x,	9,	4//3*log(1+x^(1//4))-2//3*log(1-x^(1//4)+sqrt(x))+4*atan((1-2*x^(1//4))/sqrt(3))/sqrt(3)+2*sqrt(x), SI.NotImplementedError],
[1//(1//x^(1//3)+1//x^(1//4)),	x,	4,	12*x^(1//12)-6*x^(1//6)+4*x^(1//4)-3*x^(1//3)+12//5*x^(5//12)+12//7*x^(7//12)-3//2*x^(2//3)+4//3*x^(3//4)-6//5*x^(5//6)+12//11*x^(11//12)-x+12//13*x^(13//12)-6//7*x^(7//6)+4//5*x^(5//4)-12*log(1+x^(1//12))-2*sqrt(x), SI.NotImplementedError],	# exponent 1//3, 1/4
[sqrt((1-x)//x),	x,	5,	-atan(sqrt(-1+1//x))+x*sqrt(-1+1//x), SI.NotImplementedError],
[cos(x)//(sin(x)+sin(x)^2),	x,	2,	log(sin(x))-log(1+sin(x))],
[exp(2*x)//(2+3*exp(x)+exp(2*x)),	x,	4,	-log(1+exp(x))+2*log(2+exp(x))],
[1//sqrt(1+exp(x)),	x,	3,	-2*atanh(sqrt(1+exp(x))), SI.NotImplementedError],
[sqrt(1-exp(x)),	x,	4,	-2*atanh(sqrt(1-exp(x)))+2*sqrt(1-exp(x)), SI.NotImplementedError],
[1//(3-5*sin(x)),	x,	4,	-1//4*log(cos(1//2*x)-3*sin(1//2*x))+1//4*log(3*cos(1//2*x)-sin(1//2*x))],
[1//(cos(x)+sin(x)),	x,	2,	-atanh((cos(x)-sin(x))/sqrt(2))/sqrt(2)],	# result contains algebraic numbers
[1//(1-cos(x)+sin(x)),	x,	2,	-log(1+cot(1//2*x))],
[1//(4*cos(x)+3*sin(x)),	x,	2,	-1//5*atanh(1//5*(3*cos(x)-4*sin(x)))],
[1//(sin(x)+tan(x)),	x,	6,	-1//2*atanh(cos(x))+1//2*cot(x)*csc(x)-1//2*csc(x)^2],
[1//(2*sin(x)+sin(2*x)),	x,	4,	1//4*log(tan(1//2*x))+1//8*tan(1//2*x)^2],
[sec(x)//(1+sin(x)),	x,	4,	1//2*atanh(sin(x))+(-1//2)//(1+sin(x))],
[1//(b*cos(x)+a*sin(x)),	x,	2,	-atanh((a*cos(x)-b*sin(x))/sqrt(a^2+b^2))/sqrt(a^2+b^2)],
[1//(b^2*cos(x)^2+a^2*sin(x)^2),	x,	2,	atan(a*tan(x)/b)/(a*b)],

# Section 7.6 - Strategy for Integration
[x//(-1+x^2),	x,	1,	1//2*log(1-x^2)],
[sqrt(x)*(1+sqrt(x)),	x,	2,	2//3*x^(3//2)+1//2*x^2, SI.NotImplementedError],
[1//(1-cos(x)),	x,	1,	-sin(x)//(1-cos(x))],
[sec(x)*tan(x)^2,	x,	2,	-1//2*atanh(sin(x))+1//2*sec(x)*tan(x)],
[sec(x)^3*tan(x)^3,	x,	3,	-1//3*sec(x)^3+1//5*sec(x)^5],
[exp(sqrt(x)),	x,	3,	-2*exp(sqrt(x))+2*exp(sqrt(x))*sqrt(x), SI.NotImplementedError],
[(1+x^5)//(-10*x-3*x^2+x^3),	x,	3,	19*x+3//2*x^2+1//3*x^3+3126//35*log(5-x)-1//10*log(x)-31//14*log(2+x)],
[1//(x*sqrt(log(x))),	x,	2,	2*sqrt(log(x)), SI.NotImplementedError],
[(5+2*x)//(-3+x),	x,	2,	2*x+11*log(3-x)],
[exp(exp(x)+x),	x,	2,	exp(exp(x))],
[cos(x)^2*sin(x)^2,	x,	3,	1//8*x+1//8*cos(x)*sin(x)-1//4*cos(x)^3*sin(x)],
[(-cos(x)+sin(x))//(cos(x)+sin(x)),	x,	1,	-log(cos(x)+sin(x))],
[x//sqrt(1-x^2),	x,	1,	-sqrt(1-x^2), SI.NotImplementedError],
[x^3*log(x),	x,	1,	-1//16*x^4+1//4*x^4*log(x)],
[sqrt(-2+x)//(2+x),	x,	3,	-4*atan(1//2*sqrt(-2+x))+2*sqrt(-2+x), SI.NotImplementedError],
[x//(2+x)^2,	x,	2,	2//(2+x)+log(2+x)],
[log(1+x^2),	x,	3,	-2*x+2*atan(x)+x*log(1+x^2)],
[sqrt(1+log(x))//(x*log(x)),	x,	4,	-2*atanh(sqrt(1+log(x)))+2*sqrt(1+log(x)), SI.NotImplementedError],
[(1+sqrt(x))^8,	x,	3,	-2//9*(1+sqrt(x))^9+1//5*(1+sqrt(x))^10, SI.NotImplementedError],
[sec(x)^4*tan(x)^3,	x,	3,	-1//4*sec(x)^4+1//6*sec(x)^6],
[x//(2-2*x+x^2),	x,	4,	-atan(1-x)+1//2*log(2-2*x+x^2)],
[x*asin(x),	x,	3,	-1//4*asin(x)+1//2*x^2*asin(x)+1//4*x*sqrt(1-x^2), SI.NotImplementedError], # asin
[sqrt(9-x^2)//x,	x,	4,	-3*atanh(1//3*sqrt(9-x^2))+sqrt(9-x^2), SI.NotImplementedError],
[x//(2+3*x+x^2),	x,	3,	-log(1+x)+2*log(2+x)],
[x^2*cosh(x),	x,	3,	-2*x*cosh(x)+2*sinh(x)+x^2*sinh(x)],
[(1+x+x^3)//(4*x+2*x^2+x^4),	x,	1,	1//4*log(4*x+2*x^2+x^4)],

[cos(x)//(1+sin(x)^2),	x,	2,	atan(sin(x))],
[cos(sqrt(x)),	x,	3,	2*cos(sqrt(x))+2*sin(sqrt(x))*sqrt(x), SI.NotImplementedError],
[sin(Pi*x),	x,	1,	-cos(Pi*x)//Pi, SI.NotImplementedError], #Pi
[exp(2*x)//(1+exp(x)),	x,	3,	exp(x)-log(1+exp(x))],
[exp(3*x)*cos(5*x),	x,	1,	3//34*exp(3*x)*cos(5*x)+5//34*exp(3*x)*sin(5*x)],
[cos(3*x)*cos(5*x),	x,	1,	1//4*sin(2*x)+1//16*sin(8*x)],
[1//(1+x+x^2+x^3),	x,	5,	1//2*atan(x)+1//2*log(1+x)-1//4*log(1+x^2)],
[x^2*log(1+x),	x,	3,	-1//3*x+1//6*x^2-1//9*x^3+1//3*log(1+x)+1//3*x^3*log(1+x)],
[x^5//exp(x^3),	x,	2,	(-1//3)//exp(x^3)-1//3*x^3//exp(x^3)],
[tan(4*x)^2,	x,	2,	-x+1//4*tan(4*x)],
[1//sqrt(-5+12*x+9*x^2),	x,	2,	1//3*atanh((2+3*x)//sqrt(-5+12*x+9*x^2)), SI.NotImplementedError],
[x^2*atan(x),	x,	4,	-1//6*x^2+1//3*x^3*atan(x)+1//6*log(1+x^2)],
[(1-sqrt(x))//x^(1//3),	x,	2,	3//2*x^(2//3)-6//7*x^(7//6), SI.NotImplementedError],
[1//((-1)//exp(x)+exp(x)),	x,	2,	-atanh(exp(x))],
[x//(10+2*x^2+x^4),	x,	3,	1//6*atan(1//3*(1+x^2))],
[1//(1//x^(1//3)+x),	x,	2,	3//4*log(1+x^(4//3)), SI.NotImplementedError],	# exponent 1//3
[cos(x)^4*sin(x)^2,	x,	4,	1//16*x+1//16*cos(x)*sin(x)+1//24*cos(x)^3*sin(x)-1//6*cos(x)^5*sin(x)],
[1//sqrt(5-4*x-x^2),	x,	2,	-asin(1//3*(-2-x)), SI.NotImplementedError],
[x//(1-x^2+sqrt(1-x^2)),	x,	3,	-log(1+sqrt(1-x^2)), SI.NotImplementedError],
[(1+cos(x))*csc(x),	x,	2,	log(1-cos(x))],
[exp(x)//(-1+exp(2*x)),	x,	2,	-atanh(exp(x))],
[1//(-8+x^3),	x,	6,	1//12*log(2-x)-1//24*log(4+2*x+x^2)-1//4*atan((1+x)/sqrt(3))/sqrt(3)],
[x^5*cosh(x),	x,	6,	-120*cosh(x)-60*x^2*cosh(x)-5*x^4*cosh(x)+120*x*sinh(x)+20*x^3*sinh(x)+x^5*sinh(x)],
[log(tan(x))//(cos(x)*sin(x)),	x,	1,	1//2*log(tan(x))^2],
[-2*x+x^2+x^3,	x,	1,	-x^2+1//3*x^3+1//4*x^4],
[(1+exp(x))//(1-exp(x)),	x,	3,	x-2*log(1-exp(x))],
[x//((1+x^2)*(4+x^2)),	x,	4,	1//6*log(1+x^2)-1//6*log(4+x^2)],
[1//(4-5*sin(x)),	x,	4,	-1//3*log(cos(1//2*x)-2*sin(1//2*x))+1//3*log(2*cos(1//2*x)-sin(1//2*x))],
[x*(c+x)^(1//3),	x,	2,	-3//4*c*(c+x)^(4//3)+3//7*(c+x)^(7//3), SI.NotImplementedError],	# exponent 1//3
[exp(x^(1//3)),	x,	4,	6*exp(x^(1//3))-6*exp(x^(1//3))*x^(1//3)+3*exp(x^(1//3))*x^(2//3), SI.NotImplementedError],	#exponent 1//3
[1//(4+x+sqrt(1+x)),	x,	5,	log(4+x+sqrt(1+x))-2*atan((1+2*sqrt(1+x))/sqrt(11))/sqrt(11), SI.NotImplementedError],
[(1+x^3)//(-x^2+x^3),	x,	3,	1//x+x+2*log(1-x)-log(x)],
[(-3+4*x+x^2)*sin(2*x),	x,	8,	7//4*cos(2*x)-2*x*cos(2*x)-1//2*x^2*cos(2*x)+sin(2*x)+1//2*x*sin(2*x)],
[cos(cos(x))*sin(x),	x,	2,	-sin(cos(x))],
[1//sqrt(16-x^2),	x,	1,	asin(1//4*x), SI.NotImplementedError],
[x^3//(1+x)^10,	x,	2,	1//9//(1+x)^9+(-3//8)//(1+x)^8+3//7//(1+x)^7+(-1//6)//(1+x)^6],
[cot(2*x)^3*csc(2*x)^3,	x,	3,	1//6*csc(2*x)^3-1//10*csc(2*x)^5],
[(x+sin(x))^2,	x,	6,	1//2*x+1//3*x^3-2*x*cos(x)+2*sin(x)-1//2*cos(x)*sin(x)],
[exp(atan(x))//(1+x^2),	x,	1,	exp(atan(x))],
[1//(x*(1+x^4)),	x,	4,	log(x)-1//4*log(1+x^4)],
[t^3//exp(2*t),	t,	4,	(-3//8)//exp(2*t)-3//4*t//exp(2*t)-3//4*t^2//exp(2*t)-1//2*t^3//exp(2*t)],
[sqrt(t)//(1+t^(1//3)),	t,	7,	-6*t^(1//6)-6//5*t^(5//6)+6//7*t^(7//6)+6*atan(t^(1//6))+2*sqrt(t), SI.NotImplementedError],
[sin(x)*sin(2*x)*sin(3*x),	x,	5,	-1//8*cos(2*x)-1//16*cos(4*x)+1//24*cos(6*x)],
[log(1//2*x),	x,	1,	-x+x*log(1//2*x)],
[sqrt((1+x)//(1-x)),	x,	3,	2*atan(sqrt((1+x)//(1-x)))-(1-x)*sqrt((1+x)//(1-x)), SI.NotImplementedError],
[x*log(x)//sqrt(-1+x^2),	x,	5,	atan(sqrt(-1+x^2))-sqrt(-1+x^2)+log(x)*sqrt(-1+x^2), SI.NotImplementedError],
[(a+x)//(a^2+x^2),	x,	3,	atan(x//a)+1//2*log(a^2+x^2)],
[sqrt(1+x-x^2),	x,	3,	-5//8*asin((1-2*x)/sqrt(5))-1//4*(1-2*x)*sqrt(1+x-x^2), SI.NotImplementedError],
[x^4//(16+x^10),	x,	2,	1//20*atan(1//4*x^5)],
[(2+x)//(2+x+x^2),	x,	4,	1//2*log(2+x+x^2)+3*atan((1+2*x)/sqrt(7))/sqrt(7)],
[x*sec(x)*tan(x),	x,	2,	-atanh(sin(x))+x*sec(x)],
[x//(-a^4+x^4),	x,	2,	-1//2*atanh(x^2//a^2)//a^2],
[1//(sqrt(x)+sqrt(1+x)),	x,	3,	-2//3*x^(3//2)+2//3*(1+x)^(3//2), SI.NotImplementedError],
[1//(1+(-1)//exp(x)+2*exp(x)),	x,	4,	1//3*log(1-2*exp(x))-1//3*log(1+exp(x))],
[atan(sqrt(x))//sqrt(x),	x,	2,	-log(1+x)+2*atan(sqrt(x))*sqrt(x), SI.NotImplementedError],
[log(1+x)//x^2,	x,	4,	log(x)-log(1+x)-log(1+x)//x],
[1//(-exp(x)+exp(3*x)),	x,	3,	exp(-x)-atanh(exp(x))],
[(1+cos(x)^2)//(1-cos(x)^2),	x,	4,	-x-2*cot(x)],

# Section 7.?
[1//(x*sqrt(-25+2*x)),	x,	2,	2//5*atan(1//5*sqrt(-25+2*x)), SI.NotImplementedError],
[sin(2*x)//sqrt(9-cos(x)^4),	x,	5,	-asin(1//3*cos(x)^2), SI.NotImplementedError],
[x^2//sqrt(5-4*x^2),	x,	2,	5//16*asin(2*x/sqrt(5))-1//8*x*sqrt(5-4*x^2), SI.NotImplementedError],
[x^3*sin(x),	x,	4,	6*x*cos(x)-x^3*cos(x)-6*sin(x)+3*x^2*sin(x)],
[x*sqrt(4+2*x+x^2),	x,	4,	1//3*(4+2*x+x^2)^(3//2)-3//2*asinh((1+x)/sqrt(3))-1//2*(1+x)*sqrt(4+2*x+x^2), SI.NotImplementedError],
[x*(5+x^2)^8,	x,	1,	1//18*(5+x^2)^9],
[cos(x)^2*sin(x)^5,	x,	3,	-1//3*cos(x)^3+2//5*cos(x)^5-1//7*cos(x)^7],
[cos(4*x)//exp(3*x),	x,	1,	-3//25*cos(4*x)//exp(3*x)+4//25*sin(4*x)//exp(3*x)],
[csc(1//2*x)^3,	x,	2,	-atanh(cos(1//2*x))-cot(1//2*x)*csc(1//2*x)],
[sqrt(-1+9*x^2)//x^2,	x,	3,	3*atanh(3*x//sqrt(-1+9*x^2))-sqrt(-1+9*x^2)//x, SI.NotImplementedError],
[sqrt(4-3*x^2)//x,	x,	4,	-2*atanh(1//2*sqrt(4-3*x^2))+sqrt(4-3*x^2), SI.NotImplementedError],
[exp(3*x)*x^2,	x,	3,	2//27*exp(3*x)-2//9*exp(3*x)*x+1//3*exp(3*x)*x^2],
[cos(x)*sin(x)//sqrt(1+sin(x)),	x,	3,	2//3*(1+sin(x))^(3//2)-2*sqrt(1+sin(x)), SI.NotImplementedError],
[x*asin(x^2),	x,	3,	1//2*x^2*asin(x^2)+1//2*sqrt(1-x^4), SI.NotImplementedError],	# asin
[x^3*asin(x^2),	x,	5,	-1//8*asin(x^2)+1//4*x^4*asin(x^2)+1//8*x^2*sqrt(1-x^4), SI.NotImplementedError],	# asin
[exp(x)*sech(exp(x)),	x,	2,	atan(sinh(exp(x)))],
[x^2*cos(3*x),	x,	3,	2//9*x*cos(3*x)-2//27*sin(3*x)+1//3*x^2*sin(3*x)],
[sqrt(5-4*x-x^2),	x,	3,	-9//2*asin(1//3*(-2-x))+1//2*(2+x)*sqrt(5-4*x-x^2), SI.NotImplementedError],
[x^5//(x^2+sqrt(2)),	x,	3,	1//4*x^4+log(x^2+sqrt(2))-x^2/sqrt(2), SI.NotImplementedError],	# integrand contains algebraic number
[sec(x)^5,	x,	3,	3//8*atanh(sin(x))+3//8*sec(x)*tan(x)+1//4*sec(x)^3*tan(x)],
[sin(2*x)^6,	x,	4,	5//16*x-5//32*cos(2*x)*sin(2*x)-5//48*cos(2*x)*sin(2*x)^3-1//12*cos(2*x)*sin(2*x)^5],
[cos(x)*log(sin(x))*sin(x)^2,	x,	4,	-1//9*sin(x)^3+1//3*log(sin(x))*sin(x)^3],
[1//(exp(x)*(1+2*exp(x))),	x,	3,	(-1)//exp(x)-2*x+2*log(1+2*exp(x))],
[sqrt(2+3*cos(x))*tan(x),	x,	4,	2*atanh(sqrt(2+3*cos(x))/sqrt(2))*sqrt(2)-2*sqrt(2+3*cos(x)), SI.NotImplementedError],
[x//sqrt(-4*x+x^2),	x,	3,	4*atanh(x//sqrt(-4*x+x^2))+sqrt(-4*x+x^2), SI.NotImplementedError],
[cos(x)^5,	x,	2,	sin(x)-2//3*sin(x)^3+1//5*sin(x)^5],
[x^4//exp(x),	x,	5,	(-24)//exp(x)-24*x//exp(x)-12*x^2//exp(x)-4*x^3//exp(x)-x^4//exp(x)],
[x^4//sqrt(-2+x^10),	x,	3,	1//5*atanh(x^5//sqrt(-2+x^10)), SI.NotImplementedError],
[exp(x)*cos(4+3*x),	x,	1,	1//10*exp(x)*cos(4+3*x)+3//10*exp(x)*sin(4+3*x)],
[exp(x)*log(1+exp(x)),	x,	4,	-exp(x)+(1+exp(x))*log(1+exp(x))],
[x^2*atan(x),	x,	4,	-1//6*x^2+1//3*x^3*atan(x)+1//6*log(1+x^2)],
[sqrt(-1+exp(2*x)),	x,	4,	-atan(sqrt(-1+exp(2*x)))+sqrt(-1+exp(2*x)), SI.NotImplementedError],
[exp(sin(x))*sin(2*x),	x,	4,	-2*exp(sin(x))+2*exp(sin(x))*sin(x)],
[x^2*sqrt(5-x^2),	x,	3,	25//8*asin(x/sqrt(5))-5//8*x*sqrt(5-x^2)+1//4*x^3*sqrt(5-x^2), SI.NotImplementedError],
[x^2*(1+x^3)^4,	x,	1,	1//15*(1+x^3)^5],
[cos(x)^3*sin(x)^3,	x,	3,	1//4*sin(x)^4-1//6*sin(x)^6],
[sec(x)^4*tan(x)^2,	x,	3,	1//3*tan(x)^3+1//5*tan(x)^5],
[x*sqrt(1+2*x),	x,	2,	-1//6*(1+2*x)^(3//2)+1//10*(1+2*x)^(5//2), SI.NotImplementedError],
[sin(x)^4,	x,	3,	3//8*x-3//8*cos(x)*sin(x)-1//4*cos(x)*sin(x)^3],
[tan(x)^3,	x,	2,	log(cos(x))+1//2*tan(x)^2],
[x^5*sqrt(1+x^2),	x,	3,	1//3*(1+x^2)^(3//2)-2//5*(1+x^2)^(5//2)+1//7*(1+x^2)^(7//2), SI.NotImplementedError],

#=
# David Jeffrey - Rectifying Transformations for Trig Integration (1997)
# Problem (1.2)
[3/(5-4*cos(x)),    x,  2,  x+2*atan(sin(x)/(2-cos(x)))],
# Problem (1.4)
[(1+cos(x)+2*sin(x))/(3+cos(x)^2+2*sin(x)-2*cos(x)*sin(x)), x,  -43,    -atan((2*cos(x)-sin(x))/(2+sin(x)))],
# Problem (1.5)
[(2+cos(x)+5*sin(x))/(4*cos(x)-2*sin(x)+cos(x)*sin(x)-2*sin(x)^2),  x,  -25,    -log(1-3*cos(x)+sin(x))+log(3+cos(x)+sin(x))],
# Problem (3.3)
[(3+7*cos(x)+2*sin(x))/(1+4*cos(x)+3*cos(x)^2-5*sin(x)-cos(x)*sin(x)),  x   ,   -32,    -log(1+cos(x)-2*sin(x))+log(3+cos(x)+sin(x))],
# Problem
[(-1+4*cos(x)+5*cos(x)^2)/(-1-4*cos(x)-3*cos(x)^2+4*cos(x)^3),  x,  -2, x-2*atan(sin(x)/(3+cos(x)))-2*atan((3*sin(x)+7*cos(x)*sin(x))/(1+2*cos(x)+5*cos(x)^2))],
# Problem
[(-5+2*cos(x)+7*cos(x)^2)/(-1+2*cos(x)-9*cos(x)^2+4*cos(x)^3),  x,  -2, x-2*atan(2*cos(x)*sin(x)/(1-cos(x)+2*cos(x)^2))],
# Problem (3.4)
[3/(5+4*sin(x)),    x,  2,  x+2*atan(cos(x)/(2+sin(x)))],
# Problem (3.6)
[2/(1+cos(x)^2),    x,  3,  x*sqrt(2)-atan(cos(x)*sin(x)/(1+cos(x)^2+sqrt(2)))*sqrt(2)],
# Problem (3.8)
[1/(p+q*cos(x)+r*sin(x)),   x,  3,  2*atan((r+(p-q)*tan(1/2*x))/sqrt(p^2-q^2-r^2))/sqrt(p^2-q^2-r^2), SI.NotImplementedError],

# Waldek Hebisch - email May 2013

# Problem #1
[(1-x^3+x^4-x^5+x^6)*exp(x),    x,  25, 871*exp(x)-870*exp(x)*x+435*exp(x)*x^2-145*exp(x)*x^3+36*exp(x)*x^4-7*exp(x)*x^5+exp(x)*x^6],
# Problem #2
#[(2-x^2)*exp(x/(2+x^2))/(2*x+x^3),  x,  -5,Ei(x/(2+x^2))],
#[(2+2*x+3*x^2-x^3+2*x^4)*exp(x/(2+x^2))/(2*x+x^3),  x,  -5, exp(x/(2+x^2))*(2+x^2)+Ei(x/(2+x^2))],
# Problem #3
#[(1+exp(x))*exp(x+exp(x))/(x+exp(x)),  x,   2, Ei(exp(x)+x)],
# Problem #4
[(1-3*x-x^2+x^3)*exp(1/(-1+x^2))/(1-x-x^2+x^3), x,  -6, exp(1/(-1+x^2))*(1+x)],
# Problem #5
[exp(1+1/log(x))*(-1+log(x)^2)/log(x)^2,    x,  1,  exp(1+1/log(x))*x],
[exp(x+1/log(x))*(-1+(1+x)*log(x)^2)/log(x)^2,  x,  -2, exp(x+1/log(x))*x]
=#
]

k = 0
expected_exceptions = 0
unexpected_exceptions = 0
failed = 0
passed = 0
for prob in problems
    global k, passed, failed, expected_exceptions, unexpected_exceptions
    k += 1
    println("--- #$k -----------------------------------------------------------------")
    print("∫", prob[1], "dx = ")
    if length(prob)<=4
        try
            result = integrate(prob[1], prob[2], catchNotImplementedError=false, catchAlgorithmFailedError=false)
            println(result)
            arg = 1.123456789
            err = abs(SymbolicUtils.substitute(prob[1]-Symbolics.derivative(result, prob[2]), prob[2]=>arg))
            if !isa(err, Real) || err>1e-8
                if isa(err, Real)
                    println("FAILED: wrong result, err=$err")
                else
                    println("FAILED: wrong result, complete solution should exist")
                end
                failed += 1
            else
                passed += 1
            end
        catch e
            unexpected_exceptions += 1
            failed += 1
            println("\nFAILED: unexpected exception $e")
        end
        println("expected: ", prob[4])
    else
        try
            result = integrate(prob[1], prob[2], catchNotImplementedError=false, catchAlgorithmFailedError=false)
        catch e
            if e isa prob[5]
                println(e.msg)
                expected_exceptions += 1
                passed += 1
            else
                println("\nFAILED: unexpected exception $e")
                unexpected_exceptions += 1
                failed += 1
            end
        end
        println("expected: ", prob[5])
    end
end

println("----------------------------------------------------------")
total = failed+passed
println("passed: $passed, of which expected exceptions: $expected_exceptions, failed: $failed, of which unexpected exceptions: $unexpected_exceptions, total: $total")


text:
using FileIO
using ArgParse

"""
    parse_commandline(ARGS)

Parse command line arguments and return argument dictionary
"""
function parse_commandline(ARGS)
    s = ArgParseSettings("QXContexts")

    @add_arg_table! s begin
        "input"
            help = "Input file to example"
            required = true
            arg_type = String
    end
    return parse_args(ARGS, s)
end

"""
    main(ARGS)

QXContexts entry point
"""
function main(ARGS)
    args = parse_commandline(ARGS)

    input_file = args["input"]

    @show load(input_file)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end


text:
using Documenter
using Gridap

pages = [
  "Home" => "index.md",
  "Getting Started" => "getting-started.md",
  "Gridap" => "Gridap.md",
  "Gridap.Helpers" => "Helpers.md",
  "Gridap.Io" => "Io.md",
  "Gridap.Algebra" => "Algebra.md",
  "Gridap.Arrays" => "Arrays.md",
  "Gridap.TensorValues" => "TensorValues.md",
  "Gridap.Fields" => "Fields.md",
  "Gridap.Polynomials" => "Polynomials.md",
  "Gridap.ReferenceFEs" => "ReferenceFEs.md",
  "Gridap.Geometry" => "Geometry.md",
  "Gridap.CellData" => "CellData.md",
  "Gridap.Visualization" => "Visualization.md",
  "Gridap.FESpaces" => "FESpaces.md",
  "Gridap.MultiField" => "MultiField.md",
 ]

makedocs(
    sitename = "Gridap.jl",
    format = Documenter.HTML(),
    modules = [Gridap],
    pages = pages
)

deploydocs(
    repo = "github.com/gridap/Gridap.jl.git",
)


text:
@doc raw"""
`VODE`: Variational Ordinary Differential Equation *EXPERIMENTAL*

Defines an implicit initial value problem
```math
\begin{align*}
\dot{q} (t) &= v(t) , &
q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), v(t)) , &
p(t_{0}) &= p_{0} , \\
p(t) &= ϑ(t, q(t), v(t))
\end{align*}
```
with vector field ``f``, the momentum defined by ``p``, initial conditions ``(q_{0}, p_{0})`` and the solution
``(q,p)`` taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``.
This is a special case of a differential algebraic equation with dynamical
variables ``(q,p)`` and algebraic variable ``v``.

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``f`` and ``p``
* `ϑ`: function determining the momentum
* `f`: function computing the vector field
* `g`: function determining the projection, given by ∇ϑ(q)λ
* `v`: function computing an initial guess for the velocity field (optional)
* `t₀`: initial time (optional)
* `q₀`: initial condition for `q`
* `p₀`: initial condition for `p`
* `λ₀`: initial condition for `λ`

The functions `ϑ` and `f` must have the interface
```julia
    function ϑ(t, q, v, p)
        p[1] = ...
        p[2] = ...
        ...
    end
```
and
```julia
    function f(t, q, v, f)
        f[1] = ...
        f[2] = ...
        ...
    end
```
where `t` is the current time, `q` is the current solution vector, `v` is the
current velocity and `f` and `p` are the vectors which hold the result of
evaluating the functions ``f`` and ``ϑ`` on `t`, `q` and `v`.
The funtions `g` and `v` are specified by
```julia
    function g(t, q, λ, g)
        g[1] = ...
        g[2] = ...
        ...
    end
```
and
```julia
    function v(t, q, p, v)
        v[1] = ...
        v[2] = ...
        ...
    end
```
"""
struct VODE{dType <: Number, tType <: Number, ϑType <: Function, fType <: Function, gType <: Function, vType <: Function, ωType <: Function, dHType <: Function, N} <: Equation{dType, tType}
    d::Int
    m::Int
    n::Int
    ϑ::ϑType
    f::fType
    g::gType
    v::vType
    ω::ωType
    dH::dHType
    t₀::tType
    q₀::Array{dType, N}
    p₀::Array{dType, N}
    λ₀::Array{dType, N}
    periodicity::Vector{dType}

    function VODE{dType,tType,ϑType,fType,gType,vType,ωType,dHType,N}(d, n, ϑ, f, g, v, ω, dH, t₀, q₀, p₀; periodicity=[]) where {dType <: Number, tType <: Number, ϑType <: Function, fType <: Function, gType <: Function, vType <: Function, ωType <: Function, dHType <: Function, N}
        @assert d == size(q₀,1) == size(p₀,1)
        @assert n == size(q₀,2) == size(p₀,2)
        @assert dType == eltype(q₀) == eltype(p₀)
        @assert ndims(q₀) == ndims(p₀) == N ∈ (1,2)

        λ₀ = zero(q₀)

        if !(length(periodicity) == d)
            periodicity = zeros(dType, d)
        end

        new(d, d, n, ϑ, f, g, v, ω, dH, t₀, q₀, p₀, λ₀, periodicity)
    end
end

function VODE(ϑ::ϑT, f::FT, g::GT, v::VT, ω::ΩT, dH::HT, t₀::TT, q₀::DenseArray{DT}, p₀::DenseArray{DT}; periodicity=[]) where {DT,TT,ϑT,FT,GT,VT,ΩT,HT}
    @assert size(q₀) == size(p₀)
    VODE{DT, TT, ϑT, FT, GT, VT, ΩT, HT, ndims(q₀)}(size(q₀, 1), size(q₀, 2), ϑ, f, g, v, ω, dH, t₀, q₀, p₀, periodicity=periodicity)
end

function VODE(ϑ::Function, f::Function, g::Function, ω::Function, dH::Function, t₀::Number, q₀::DenseArray, p₀::DenseArray; periodicity=[])
    VODE(ϑ, f, g, function_v_dummy, ω, dH, t₀, q₀, p₀, periodicity=periodicity)
end

function VODE(ϑ::Function, f::Function, g::Function, v::Function, ω::Function, dH::Function, q₀::DenseArray, p₀::DenseArray; periodicity=[])
    VODE(ϑ, f, g, v, ω, dH, zero(Float64), q₀, p₀, periodicity=periodicity)
end

function VODE(ϑ::Function, f::Function, g::Function, ω::Function, dH::Function, q₀::DenseArray, p₀::DenseArray; periodicity=[])
    VODE(ϑ, f, g, function_v_dummy, ω, dH, zero(Float64), q₀, p₀, periodicity=periodicity)
end

Base.hash(ode::VODE, h::UInt) = hash(ode.d, hash(ode.n, hash(ode.ϑ, hash(ode.f, hash(ode.g, hash(ode.v, hash(ode.ω, hash(ode.t₀, hash(ode.q₀, hash(ode.p₀, hash(ode.λ₀, hash(ode.periodicity, h))))))))))))
Base.:(==)(ode1::VODE, ode2::VODE) = (
                                ode1.d == ode2.d
                             && ode1.n == ode2.n
                             && ode1.ϑ == ode2.ϑ
                             && ode1.f == ode2.f
                             && ode1.g == ode2.g
                             && ode1.v == ode2.v
                             && ode1.ω == ode2.ω
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀
                             && ode1.p₀ == ode2.p₀
                             && ode1.λ₀ == ode2.λ₀
                             && ode1.periodicity == ode2.periodicity)

Base.ndims(ode::VODE) = ode.d


text:
# #=
# Copyright 2022 INSIGNEO Institute for in silico Medicine
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# =#

"""
    solveBifurcation(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)

Solve the non-linear system at the bifurcation node between three vessels.
"""
function solveBifurcation(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)
    @fastmath @inbounds U0 = @SArray [v1.u[end],
       v2.u[1],
       v3.u[1],
       sqrt(sqrt(v1.A[end])),
       sqrt(sqrt(v2.A[1])),
       sqrt(sqrt(v3.A[1]))]

    @inbounds k1 = v1.s_15_gamma[end]
    @inbounds k2 = v2.s_15_gamma[1]
    @inbounds k3 = v3.s_15_gamma[1]
    k = @SArray [k1, k2, k3]

    J = calculateJacobianBifurcation(v1, v2, v3, U0, k)
    U = newtonRaphson([v1, v2, v3], J, U0, k, calculateWstarBifurcation, calculateFBifurcation)

    updateBifurcation(U, v1, v2, v3)
end

"""
    calculateJacobianBifurcation(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel, U, k)

Return the Jacobian for bifurcation equations.
"""
function calculateJacobianBifurcation(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel, U, k)
    @fastmath @inbounds U43 = U[4]*U[4]*U[4]
    @fastmath @inbounds U53 = U[5]*U[5]*U[5]
    @fastmath @inbounds U63 = U[6]*U[6]*U[6]

    @fastmath @inbounds J14 = 4.0*k[1]
    @fastmath @inbounds J25 = -4.0*k[2]
    @fastmath @inbounds J36 = -4.0*k[3]

    @fastmath @inbounds J41 = U[4]*U43
    @fastmath @inbounds J42 = -U[5]*U53
    @fastmath @inbounds J43 = -U[6]*U63
    @fastmath @inbounds J44 = 4.0*U[1]*U43
    @fastmath @inbounds J45 = -4.0*U[2]*U53
    @fastmath @inbounds J46 = -4.0*U[3]*U63

    @fastmath @inbounds J54 = 2.0*v1.beta[end]*U[4]*v1.s_inv_A0[end]
    @fastmath @inbounds J55 = -2.0*v2.beta[1]*U[5]*v2.s_inv_A0[1]

    @fastmath @inbounds J64 = 2.0*v1.beta[end]*U[4]*v1.s_inv_A0[end]
    @fastmath @inbounds J66 = -2.0*v3.beta[1]*U[6]*v3.s_inv_A0[1]

    return @SArray [1.0 0.0 0.0 J14 0.0 0.0;
                    0.0 1.0 0.0 0.0 J25 0.0;
                    0.0 0.0 1.0 0.0 0.0 J36;
                    J41 J42 J43 J44 J45 J46;
                    0.0 0.0 0.0 J54 J55 0.0;
                    0.0 0.0 0.0 J64 0.0 J66]
end

"""
    calculateWstarBifurcation(U, k)

Return the Riemann invariants at the bifurcation node.
"""
function calculateWstarBifurcation(U, k)
    @fastmath @inbounds W1 = U[1] + 4.0*k[1]*U[4]
    @fastmath @inbounds W2 = U[2] - 4.0*k[2]*U[5]
    @fastmath @inbounds W3 = U[3] - 4.0*k[3]*U[6]

    return @SArray [W1, W2, W3]
end

"""
    calculateFBifurcation(vessels :: Array{Vessel,1}, U, k, W)
"""
function calculateFBifurcation(vessels :: Array{Vessel,1}, U, k, W)
    v1 = vessels[1]
    v2 = vessels[2]
    v3 = vessels[3]

    @fastmath @inbounds U42 = U[4]*U[4]
    @fastmath @inbounds U52 = U[5]*U[5]
    @fastmath @inbounds U62 = U[6]*U[6]

    @fastmath @inbounds f1 = U[1] + 4*k[1]*U[4] - W[1]

    @fastmath @inbounds f2 = U[2] - 4*k[2]*U[5] - W[2]

    @fastmath @inbounds f3 = U[3] - 4*k[3]*U[6] - W[3]

    @fastmath @inbounds f4 = U[1]*(U42*U42) - U[2]*(U52*U52) - U[3]*(U62*U62)

    f5 = v1.beta[end]*(U42*v1.s_inv_A0[end] - 1.0) -
        (v2.beta[1]*(U52*v2.s_inv_A0[1] - 1.0))

    f6 = v1.beta[end]*(U42*v1.s_inv_A0[end] - 1.0) -
        (v3.beta[1]*(U62*v3.s_inv_A0[1] - 1.0))

    return @SArray [f1, f2, f3, f4, f5, f6]
end

"""
    updateBifurcation(U, v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)

Update the values at the bifurcation node for the three vessels.
"""
function updateBifurcation(U, v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)
    @inbounds v1.u[end] = U[1]
    @inbounds v2.u[1] = U[2]
    @inbounds v3.u[1] = U[3]

    @fastmath @inbounds v1.A[end] = U[4]*U[4]*U[4]*U[4]
    @fastmath @inbounds v2.A[1] = U[5]*U[5]*U[5]*U[5]
    @fastmath @inbounds v3.A[1] = U[6]*U[6]*U[6]*U[6]

    @fastmath @inbounds v1.Q[end] = v1.u[end]*v1.A[end]
    @fastmath @inbounds v2.Q[1] = v2.u[1]*v2.A[1]
    @fastmath @inbounds v3.Q[1] = v3.u[1]*v3.A[1]

    @inbounds v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
    @inbounds v2.P[1] = pressure(v2.A[1], v2.A0[1], v2.beta[1], v2.Pext)
    @inbounds v3.P[1] = pressure(v3.A[1], v3.A0[1], v3.beta[1], v3.Pext)

    @fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
    @fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1])
    @fastmath @inbounds v3.c[1] = waveSpeed(v3.A[1], v3.gamma[1])
end


text:
## polynomial model smoother, borrowed from Friedman's super smoother (Firedman 1984):
## fit a polynomial model to the data in a sliding window of length approximately
## `span_frac * 2 * (maximum(x) - minimum(x))` on x axis,
## set the y value at the middle point of the sliding window
## as the predicted value by the polynomial model.

# Top-level function
function pnmsmu(
    X ::AbstractVector,
    Y ::AbstractVector,
    span_frac ::Real=0.05,
    polynomial_degree ::Integer=1, ## polynomial_degree <= 1 will yield same results as == 1
    x_sortedb ::Bool=true,
    x_evenly_spacedb ::Bool=true
)
    if ndims(X) > 1 || ndims(Y) > 1
        errmsg = "X and Y must be 1-dimensional"
        warn(logger, errmsg)
        throw(ArgumentError(errmsg))
    end

    dlen = length(X)

    if dlen < 5 || length(Y) != dlen
        errmsg = "X and Y must be equal and >= 5 in length"
        warn(logger, errmsg)
        throw(ArgumentError(errmsg))
    end

    if !x_sorted
        indice_sorted = sortperm(X)
        indice_back = sortperm(indice_sorted)
        X = X[indice_sorted]
        Y = Y[indice_sorted]
    end

    # sel_idc = selected indices
    if x_evenly_spaced ## faster
        span_dp = round(Int_T, span_frac * dlen) ## dp = data points
        sel_idc_vec = map(i -> giis_even(dlen, i, span_dp), 1:dlen)
    else
        span_x = span_frac * (maximum(X) - minimum(X))
        sel_idc_vec = map(i -> giis_uneven(X, i, span_x), 1:dlen)
    end

    ## construct `x_pn_df` for to use `pnm_lm`
    # x_pn_df = DataFrame(d1 = X)
    # for d in 2:polynomial_degree ## nothing happends when polynomial_degree == 1
    #     x_pn_df[parse("d$d")] = X .^ d
    # end

    x_ary = hcat(map(d -> X .^ d, 1:polynomial_degree)...)

    Y_smu = map(1:dlen) do i
        sel_idc = sel_idc_vec[i]
        coefs = pnm_lsq(x_ary[sel_idc,:], Y[sel_idc])
        local_x_vec = hcat(1, x_ary[i,:])
        return *(local_x_vec, coefs)[1] ## *(array_1_m, vector_m) == vector_1
    end ## do i

    return x_sorted ? Y_smu : Y_smu[indice_back]
end

## functions called by 'pnmsmu'

## functions to fit a polynomial model on using x to predict y

## using `lm` from "GLM". ~14 sec
function pnm_lm(x_pn_df ::DataFrame, Y ::AbstractVector)
    xy_df = copy(x_pn_df)
    xy_df[:y] = Y
    fml = Formula(:y, Expr(:call, :+, 1, x_pn_df.colindex.names...)) ## note: Formula can also be constructed as eval(parse(computated_str))
    return coef(lm(fml, xy_df))
end

## using native least square method. ~4 sec when using DataFrame, 2-3 sec when using Array
function pnm_lsq(x_ary ::AbstractArray, Y ::AbstractVector)
    x_w1 = hcat(ones(size(x_ary, 1)), x_ary)
    x_w1_T = transpose(x_w1)
    coefs = *(inv(x_w1_T * x_w1), x_w1_T, Y)
end

#

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using MatrixNetworks

type Biconnected_components_output
    bcc_edge_list::Vector{Int64}  #list of bcc edges
    components::Vector{Int64} #Indices of bcc corresponding to the edges in List
    number::Int64 # Total number of Biconnected components
    A::MatrixNetwork #MatrixNetwork
end

I = [1,2,2,5,5,1,5,3,3,6,6,5,4,5]
J = [2,1,5,2,1,5,3,5,6,3,5,6,5,4]
V = [1,1,1,1,1,1,1,1,1,1,1,1,1,1]

A = sparse(I,J,V)

function BiCC_util(A::MatrixNetwork, vertex::Int64, bcc_edges::Vector{Int64}, low::Vector{Int64}, dt::Vector{Int64}, pred::Vector{Int64}, components::Vector{Int64}, cs::Vector{Int64}, cn::Int64, t::Int64, css::Int64, bcc_edges_count::Int64)
    t = t+1
    dt[vertex] = t
    low[vertex] = t
    children = 0
    for i=A.rp[vertex]:A.rp[vertex+1]-1
        adj = A.ci[i];
        if dt[adj] == 0
            children += 1
            pred[adj] = vertex
            css = css+1
            cs[2*css-1] = vertex
            cs[2*css] = adj
            (cn,t,css,bcc_edges_count) = BiCC_util(A, adj, bcc_edges, low, dt, pred, components, cs, cn, t, css, bcc_edges_count)
            if low[vertex] > low[adj]
                low[vertex] = low[adj]
            end
            if ((dt[vertex] == 1 && children > 1) || (dt[vertex] > 1 && low[adj] >= dt[vertex]))
                while css > 0
                    bcc_edges_count += 1
                    bcc_edges[2*bcc_edges_count-1] = cs[2*css-1]
                    bcc_edges[2*bcc_edges_count] = cs[2*css]
                    components[bcc_edges_count] = cn
                    if (cs[2*css-1] == vertex && cs[2*css] == adj)
                        css -= 1
                        break
                    end
                    css -= 1
                end
                cn += 1
            end
        elseif (adj != pred[vertex] && dt[adj] < low[vertex])
            low[vertex] = dt[adj]
            css = css+1
            cs[2*css-1] = vertex
            cs[2*css] = adj
        end
    end
    return (cn,t,css,bcc_edges_count)
end

function BiCC(A::MatrixNetwork)
    bcc_edges=zeros(Int64,2*length(A.ci))
    bcc_edges_count=0
    cn=1
    low=zeros(Int64,A.n)
    dt=zeros(Int64,A.n)
    pred=zeros(Int64,A.n)
    components=-1 * ones(Int64,length(A.ci))
    t=0
    cs=zeros(Int64,2*length(A.ci))
    css=0
    for i=1:A.n
        if dt[i] == 0
            (cn,t,css,bcc_edges_count) = BiCC_util(A, i, bcc_edges, low, dt, pred, components, cs, cn, t, css, bcc_edges_count)
        end
    end
    while css > 0
        bcc_edges_count += 1
        bcc_edges[2*bcc_edges_count-1] = cs[2*css-1]
        bcc_edges[2*bcc_edges_count] = cs[2*css]
	components[bcc_edges_count] = cn
        css -= 1
    end
    print(bcc_edges)
    print(components)
end

function dfs_test()
    print(A);
B = MatrixNetwork(A);
print(B)
   print(BiCC(B))
end

dfs_test();


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

__precompile__(true)

module SolidStateDetectors

using LinearAlgebra
using Random
using Statistics

using ArraysOfArrays
using Interpolations
using IntervalSets
using JSON
using LaTeXStrings
using NamedTupleTools
using ParallelProcessingTools
using ProgressMeter
using RecipesBase
using StaticArrays
using Unitful
using YAML

import Clustering
import Distributions
import Tables
import TypedTables

import Base: size, sizeof, length, getindex, setindex!, axes, getproperty,
             range, ndims, eachindex, enumerate, iterate, IndexStyle, eltype, in
import Base: show, print, println, display, +, -, &
import Base.convert

const SSD = SolidStateDetectors; export SSD
export SolidStateDetector
export SSD_examples

export Grid, CylindricalPoint, CartesianPoint

export ElectricPotential, PointTypes, ChargeDensity, DielectricDistribution, WeightingPotential, ElectricField
export apply_initial_state!
export calculate_electric_potential!, calculate_weighting_potential!, calculate_electric_field!
export update_till_convergence!, refine!
export set_charge_drift_model!, apply_charge_drift_model!
export get_active_volume
export generate_charge_signals, generate_charge_signals!
export VacuumChargeDriftModel, ADLChargeDriftModel
export Simulation, simulate!
export Event, drift_charges!

const SSDFloat = Union{Float16, Float32, Float64}

struct ConfigFileError <: Exception
    msg::AbstractString
end
Base.showerror(io::IO, e::ConfigFileError) = print(io, "ConfigFileError: ", e.msg)

include("Geometries/Geometries.jl")

include("Axes/DiscreteAxis.jl")
include("World/World.jl")
include("Grids/Grids.jl")

include("Types/Types.jl")

include("MaterialProperties/MaterialProperties.jl")
include("Config/Config.jl")
include("ChargeDensityModels/ChargeDensityModels.jl")
include("SolidStateDetector/DetectorGeometries.jl")
include("GeometryRounding.jl")

include("PotentialSimulation/PotentialSimulation.jl")

include("ElectricField/ElectricField.jl")

include("ChargeDriftModels/ChargeDriftModels.jl")
include("ChargeDrift/ChargeDrift.jl")
include("SignalGeneration/SignalGeneration.jl")

include("ChargeStatistics/ChargeStatistics.jl")
include("ChargeClustering/ChargeClustering.jl")

include("Simulation/Simulation.jl")
include("Event/Event.jl")

include("IO/IO.jl")

include("examples.jl")

end # module


text:
using Documenter, Literate
using CUDA

const src = "https://github.com/JuliaGPU/CUDA.jl"
const dst = "https://juliagpu.gitlab.io/CUDA.jl/"

function main()
    @info "Building Literate.jl documentation"
    cd(@__DIR__) do
        Literate.markdown("src/tutorials/introduction.jl", "src/tutorials";
                          repo_root_url="$src/blob/master/docs")
    end

    @info "Generating Documenter.jl site"
    makedocs(
        sitename = "CUDA.jl",
        authors = "Tim Besard",
        repo = "$src/blob/{commit}{path}#{line}",
        format = Documenter.HTML(
            # Use clean URLs on CI
            prettyurls = get(ENV, "CI", nothing) == "true",
            canonical = dst,
            assets = ["assets/favicon.ico"],
            analytics = "UA-154489943-2",
        ),
        doctest = true,
        #strict = true,
        modules = [CUDA],
        pages = Any[
            "Home" => "index.md",
            "Tutorials" => Any[
                "tutorials/introduction.md",
            ],
            "Installation" => Any[
                "installation/overview.md",
                "installation/conditional.md",
                "installation/troubleshooting.md",
            ],
            "Usage" => Any[
                "usage/overview.md",
                "usage/workflow.md",
                "usage/array.md",
                "usage/memory.md",
                "usage/multigpu.md",
            ],
            "Development" => Any[
                "development/profiling.md",
            ],
            "API reference" => Any[
                "api/essentials.md",
                "api/compiler.md",
                "api/kernel.md",
                "api/array.md",
            ],
            "Library reference" => Any[
                "lib/driver.md",
            ],
            "FAQ" => "faq.md",
        ]
    )
end

isinteractive() || main()


text:
"""
    compute_signal(M, ξ) = sum(M * ξ)

Compute signal from magnetization `ξ`, using the mass matrix `M` for integration.

Given a mesh `mesh` and a vector of compartment mass matrices `M_cmpts`, a vector of
compartment signals may be obtained by `compute_signal.(M_cmpts, split_field(mesh, ξ))`.
"""
compute_signal(M, ξ) = sum(M * ξ)


text:
# # Multivariate polynomials implementations

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Extension/typed.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Extension/typed.ipynb)
# **Contributed by**: Benoît Legat

# The SumOfSquares package is built on top of the [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl)
# abstract interface. [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl/)
# is an implementation of this abstract interface so it can be used with
# SumOfSquares. Moreover, any other implementation can be used as well. To
# illustrate, we solve Examples 3.38 of [BPT12] with
# [TypedPolynomials](https://github.com/JuliaAlgebra/TypedPolynomials.jl),
# another implementation of [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl).
#
# [BPT12] Blekherman, G. & Parrilo, P. A. & Thomas, R. R.
# *Semidefinite Optimization and Convex Algebraic Geometry*.
# Society for Industrial and Applied Mathematics, **2012**.

using Test #src
import TypedPolynomials
TypedPolynomials.@polyvar x y
using SumOfSquares
import CSDP
model = SOSModel(CSDP.Optimizer)
con_ref = @constraint(model, 2x^4 + 5y^4 - x^2*y^2 >= -2(x^3*y + x + 1))
optimize!(model)
solution_summary(model)

# We see that the problem is feasible. The Sum-of-Squares decomposition can be
# obtained as follows:

sos_decomposition(con_ref)

# Why is there several implementations ?
# Depending in the use-case, one implementation may be more appropriate than
# another one. [TypedPolynomials](https://github.com/JuliaAlgebra/TypedPolynomials.jl)
# is faster than [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl/)
# but it requires new compilation whenever the list of variables changes.
# This means that [TypedPolynomials](https://github.com/JuliaAlgebra/TypedPolynomials.jl)
# is not appropriate when the number of variables is dynamic or too large.
# However, for a small number of variables, it can be faster.
# When solving Sum-of-Squares programs, the time is mostly taken by the Semidefinite programming solver.
# The time taken by SumOfSquares/JuMP/MathOptInterface are usually negligible
# or it time is taken by manipulation of JuMP or MathOptInterface functions
# therefore using TypedPolynomials over DynamicPolynomials may not make much difference in most cases.

# One case for which using TypedPolynomials might be adequate is when
# using domain defined by equalities (possibly also with inequalities).
# Indeed, in that case, SumOfSquares computes the corresponding Gröbner basis which
# may take a non-negligible amount of time for large systems of equalities.

# To illustrate this, consider the computation of Gröbner basis for the
# following system from [CLO05, p. 17].
# The time taken by TypedPolynomials is below:
#
# [CLO05] Cox, A. David & Little, John & O'Shea, Donal
# *Using Algebraic Geometry*.
# Graduate Texts in Mathematics, **2005**.
# https://doi.org/10.1007/b138611

using BenchmarkTools
@btime let
    TypedPolynomials.@polyvar x y
    S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y
    SemialgebraicSets.computegröbnerbasis!(S.I)
end

# The time taken by DynamicPolynomials is as follows:

import DynamicPolynomials
@btime let
    DynamicPolynomials.@polyvar x y
    S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y
    SemialgebraicSets.computegröbnerbasis!(S.I)
end

# We see that TypedPolynomials is faster.
# The time is still negligible for this small system but for larger systems, choosing TypedPolynomials may be helpful.
# We can use this system in a Sum-of-Squares constraint as follows:

TypedPolynomials.@polyvar x y
S = @set x^3 * y + x == 2x^2 * y^2 && 3x^4 == y
poly = -6x - 4y^3 + 2x*y^2 + 6x^3 - 3y^4 + 13x^2 * y^2
model = Model(CSDP.Optimizer)
con_ref = @constraint(model, poly in SOSCone(), domain = S)
optimize!(model)
solution_summary(model)

# We obtain the following decomposition:

dec = sos_decomposition(con_ref, 1e-6)

# We can verify that it is correct as follows:

# TODO remove `polynomial` #src
@test SumOfSquares.MP.isapproxzero(rem(dec - poly, S.I), ztol = 1e-4) #src
rem(dec - poly, S.I)

# Note that the difference between `dec` and `poly` is larger
# than between the full gram matrix because `dec` is obtained by dropping
# the lowest eigenvalues with the threshold `1e-6`; see [`sos_decomposition`](@ref).

@test SumOfSquares.MP.isapproxzero(rem(gram_matrix(con_ref) - poly, S.I), ztol = 1e-6) #src
rem(gram_matrix(con_ref) - poly, S.I)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares