text:
using Test
using Printf: @sprintf
using StaticArrays: SVector
using LinearAlgebra: dot
using QuantumLattices.Essentials.DegreesOfFreedom
using QuantumLattices.Essentials.QuantumOperators: ID, sequence, id, idtype
using QuantumLattices.Essentials.Spatials: AbstractPID, PID, CPID, Point, Bond, Bonds, Lattice, pidtype, rcoord, icoord
using QuantumLattices.Essentials: kind, update!, reset!
using QuantumLattices.Interfaces: decompose, rank, expand, ⊗
using QuantumLattices.Prerequisites: Float, decimaltostr
using QuantumLattices.Prerequisites.Traits: parameternames, isparameterbound, contentnames, getcontent
using QuantumLattices.Prerequisites.CompositeStructures: NamedContainer

import QuantumLattices.Essentials.DegreesOfFreedom: ishermitian, statistics, latexname, script, couplingcenters, abbr
import QuantumLattices.Prerequisites.VectorSpaces: shape, ndimshape

struct DID{N<:Union{Int, Symbol}} <: SimpleIID
    nambu::N
end
@inline Base.adjoint(sl::DID{Int}) = DID(3-sl.nambu)
@inline statistics(::Type{<:DID}) = :f

struct DFock <: SimpleInternal{DID{Int}}
    nnambu::Int
end
@inline shape(f::DFock) = (1:f.nnambu,)
@inline ndimshape(::Type{DFock}) = 1
@inline DID(i::CartesianIndex, vs::DFock) = DID(i.I...)
@inline CartesianIndex(did::DID{Int}, vs::DFock) = CartesianIndex(did.nambu)
@inline shape(iidspace::IIDSpace{DID{Symbol}, DFock}) = (1:iidspace.internal.nnambu,)
@inline shape(iidspace::IIDSpace{DID{Int}, DFock}) = (iidspace.iid.nambu:iidspace.iid.nambu,)
@inline shape(iidspace::IIDSpace{DID{Symbol}, DFock}) = (1:iidspace.internal.nnambu,)
@inline shape(iidspace::IIDSpace{DID{Int}, DFock}) = (iidspace.iid.nambu:iidspace.iid.nambu,)

function Base.angle(id::OID{<:Index{<:AbstractPID, DID{Int}}}, vectors::AbstractVector{<:AbstractVector{Float}}, values::AbstractVector{Float})
    phase=  (length(vectors) == 1) ? 2pi*dot(decompose(id.icoord, vectors[1]), values) :
            (length(vectors) == 2) ? 2pi*dot(decompose(id.icoord, vectors[1], vectors[2]), values) :
            (length(vectors) == 3) ? 2pi*dot(decompose(id.icoord, vectors[1], vectors[2], vectors[3]), values) :
            error("angle error: not supported number of input basis vectors.")
    (id.index.iid.nambu == 1) ? phase : -phase
end
@inline script(::Val{:site}, index::Index{<:AbstractPID, <:DID}; kwargs...) = index.pid.site
@inline script(::Val{:nambu}, index::Index{<:AbstractPID, <:DID}; kwargs...) = index.iid.nambu==2 ? "\\dagger" : ""

latexname(::Type{<:OID{<:Index{<:AbstractPID, <:DID}}}) = Symbol("OID{Index{AbstractPID, DID}}")
latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:site,)}('d'))

const DCoupling{V, I<:ID{DID}, C<:Subscripts, CI<:SubscriptsID} = Coupling{V, I, C, CI}
@inline Base.repr(tc::(Coupling{V, <:ID{DID}} where V)) = @sprintf "%s ph(%s)" decimaltostr(tc.value) join(tc.cid.nambus, ", ")
@inline couplingcenters(::(Coupling{V, <:ID{DID}} where V), ::Bond, ::Val) = (1, 2)
@inline DCoupling(value, nambus::Tuple{Vararg{Int}}) = Coupling(value, ID(DID, nambus), Subscripts((nambu=Subscript(nambus),)))
@inline DCoupling(value, nambus::Subscript) = Coupling(value, ID(DID, convert(Tuple, nambus)), Subscripts((nambu=nambus,)))

@inline abbr(::Type{<:Term{:Mu}}) = :mu
@inline abbr(::Type{<:Term{:Hp}}) = :hp
@inline ishermitian(::Type{<:Term{:Mu}}) = true
@inline ishermitian(::Type{<:Term{:Hp}}) = false

@testset "IID" begin
    did = DID(1)
    @test statistics(did) == statistics(typeof(did)) == :f

    did₁, did₂ = DID(1), DID(2)
    ciid = CompositeIID(did₁, did₂)
    @test length(ciid) == length(typeof(ciid)) == 2
    @test rank(ciid) == rank(typeof(ciid)) == 2
    @test iidtype(ciid, 1) == iidtype(typeof(ciid), 1) == DID{Int}
    @test iidtype(ciid, 2) == iidtype(typeof(ciid), 2) == DID{Int}
    @test ciid[1]==did₁ && ciid[2]==did₂
    @test ciid.nambus == (1, 2)
    @test string(ciid) == "DID(1) ⊗ DID(2)"
    @test did₁⊗did₂ == ciid
    @test did₁⊗ciid == CompositeIID(did₁, did₁, did₂)
    @test ciid⊗did₁ == CompositeIID(did₁, did₂, did₁)
    @test ciid⊗ciid == CompositeIID(did₁, did₂, did₁, did₂)
end

@testset "SimpleInternal" begin
    it = DFock(2)
    @test it|>eltype == DID{Int}
    @test it|>typeof|>eltype == DID{Int}
    @test it == deepcopy(it)
    @test isequal(it, deepcopy(it))
    @test it|>string == "DFock(nnambu=2)"
    @test it|>collect == [DID(1), DID(2)]
    @test statistics(it) == statistics(typeof(it)) == :f
    @test match(DID(1), it) && match(DID, DFock)
    @test filter(DID(1), it) == filter(DID, it) == it
    @test filter(DID(1), DFock) == filter(DID, DFock) == DFock
end

@testset "CompositeInternal" begin
    it₁, it₂ = DFock(2), DFock(3)
    ci = CompositeInternal(it₁, it₂)
    @test eltype(ci) == eltype(typeof(ci)) == CompositeIID{Tuple{DID{Int}, DID{Int}}}
    @test rank(ci) == rank(typeof(ci)) == 2
    @test string(ci) == "DFock(nnambu=2) ⊗ DFock(nnambu=3)"
    @test shape(ci) == (1:2, 1:3)
    @test ndimshape(ci) == 2
    for i = 1:2, j = 1:3
        @test CartesianIndex(CompositeIID(CartesianIndex(i, j), ci), ci) == CartesianIndex(i, j)
    end
    @test it₁⊗it₂ == ci
    @test it₁⊗ci == CompositeInternal(it₁, it₁, it₂)
    @test ci⊗it₁ == CompositeInternal(it₁, it₂, it₁)
    @test ci⊗ci == CompositeInternal(it₁, it₂, it₁, it₂)
    @test filter(DID(1), ci) == filter(DID, ci) == ci
    @test filter(DID(1), typeof(ci)) == filter(DID, typeof(ci)) == typeof(ci)
end

@testset "Index" begin
    index = Index(CPID('S', 4), DID(1))
    @test index|>pidtype == CPID{Char}
    @test index|>typeof|>pidtype == CPID{Char}
    @test index|>iidtype == DID{Int}
    @test index|>typeof|>iidtype == DID{Int}
    @test index|>adjoint == Index(CPID('S', 4), DID(2))
    @test statistics(index) == statistics(typeof(index)) == :f
    @test ishermitian(ID(index', index)) == true
    @test ishermitian(ID(index, index)) == false
end

@testset "OID" begin
    @test contentnames(CompositeOID) == (:index,)
    @test contentnames(OID) == (:index, :rcoord, :icoord)

    oid = OID(Index(PID(1), DID(1)), [0.0, -0.0], [0.0, 0.0])
    @test indextype(oid) == indextype(typeof(oid)) == Index{PID, DID{Int}}
    @test statistics(oid) == statistics(typeof(oid)) == :f
    @test oid' == OID(Index(PID(1), DID(2)), rcoord=SVector(0.0, 0.0), icoord=SVector(0.0, 0.0))
    @test hash(oid, UInt(1)) == hash(OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 1.0)), UInt(1))
    @test propertynames(ID(oid)) == (:indexes, :rcoords, :icoords)
    @test fieldnames(OID) == (:index, :rcoord, :icoord)
    @test string(oid) == "OID(Index(PID(1), DID(1)), [0.0, 0.0], [0.0, 0.0])"
    @test ID(oid', oid)' == ID(oid', oid)
    @test ishermitian(ID(oid', oid)) == true
    @test ishermitian(ID(oid, oid)) == false
    @test oidtype(DFock, Point{2, PID, Float}, Val(:default)) == OID{Index{PID, DID{Int}}, SVector{2, Float}}
end

@testset "Operator" begin
    opt = Operator(1.0im, ID(
        OID(Index(PID(2), DID(2)), SVector(1.0, 0.0), SVector(2.0, 0.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
        ))
    @test opt' == Operator(-1.0im, ID(
        OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
        OID(Index(PID(2), DID(1)), SVector(1.0, 0.0), SVector(2.0, 0.0))
        ))
    @test ishermitian(opt) == false
    @test string(opt) == "Operator(1.0im, ID(OID(Index(PID(2), DID(2)), [1.0, 0.0], [2.0, 0.0]), OID(Index(PID(1), DID(1)), [0.0, 0.0], [0.0, 0.0])))"

    opt = Operator(1.0, ID(
        OID(Index(PID(1), DID(2)), SVector(0.5, 0.5), SVector(1.0, 1.0)),
        OID(Index(PID(1), DID(1)), SVector(0.5, 0.5), SVector(1.0, 1.0))
        ))
    @test opt' == opt
    @test ishermitian(opt) == true

    opt = Operator(1.0, ID(
        OID(Index(PID(1), DID(2)), SVector(0.5, 0.5), SVector(1.0, 1.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.5), SVector(0.0, 1.0))
        ))
    @test rcoord(opt) == SVector(0.5, 0.0)
    @test icoord(opt) == SVector(1.0, 0.0)

    opt = Operator(1.0, ID(OID(Index(PID(1), DID(2)), SVector(0.5, 0.0), SVector(1.0, 0.0))))
    @test rcoord(opt) == SVector(0.5, 0.0)
    @test icoord(opt) == SVector(1.0, 0.0)
end

@testset "Operators" begin
    opt1 = Operator(1.0im, ID(
        OID(Index(PID(2), DID(2)), SVector(1.0, 0.0), SVector(0.0, 0.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(2.0, 0.0))
        ))
    opt2 = Operator(1.0, ID(
        OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
        ))
    opts = Operators(opt1, opt2)
    @test summary(opts) == "Operators{$(opts|>valtype)}"
    @test opts' == Operators(opt1', opt2')
    @test opts'+opts == Operators(opt1, opt1', opt2*2)
    @test ishermitian(opts) == false
    @test ishermitian(opts'+opts) == true
end

@testset "Hilbert" begin
    map = pid->DFock((pid.site-1)%2+1)
    hilbert = Hilbert(map, [CPID(1, 1), CPID(1, 2)])
    @test hilbert == Hilbert{DFock}(map, [CPID(1, 1), CPID(1, 2)])
    @test convert(Dict, hilbert) == Dict(CPID(1, 1)=>DFock(1), CPID(1, 2)=>DFock(2))
    reset!(hilbert, (CPID(2, 1), CPID(2, 2)))
    @test convert(Dict, hilbert) == Dict(CPID(2, 1)=>DFock(1), CPID(2, 2)=>DFock(2))

    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    @test hilbert[PID(1)] == hilbert[PID(2)] == DFock(2)

    hilbert = Hilbert(PID(1)=>DFock(2), PID(2)=>DFock(3))
    @test hilbert[PID(1)]==DFock(2) && hilbert[PID(2)]==DFock(3)
end

@testset "IIDSpace" begin
    DID₁, DID₂, it = DID(2), DID(:σ), DFock(2)
    iidspace = IIDSpace(DID₁⊗DID₂, it⊗it)
    @test eltype(iidspace) == eltype(typeof(iidspace)) == CompositeIID{Tuple{DID{Int}, DID{Int}}}
    @test kind(iidspace) == kind(typeof(iidspace)) == :info
    @test shape(iidspace) == (2:2, 1:2)
    @test ndimshape(iidspace) == ndimshape(typeof(iidspace)) == 2
    for i = 1:length(iidspace)
        iid = iidspace[i]
        @test iidspace[CartesianIndex(iid, iidspace)] == iid
    end
    @test expand((DID₁, DID₂), (it, it)) == iidspace
    @test collect(iidspace) == [DID(2)⊗DID(1), DID(2)⊗DID(2)]
end

@testset "Metric" begin
    valtype(Metric, AbstractOID) = AbstractOID

    index = Index(CPID('S', 4), DID(1))
    oid = OID(index, SVector(0.5, 0.0), SVector(1.0, 0.0))

    m = OIDToTuple((:scope, :site, :nambu))
    @test m == OIDToTuple(:scope, :site, :nambu)
    @test isequal(m, OIDToTuple(:scope, :site, :nambu))
    @test keys(m) == keys(typeof(m)) == (:scope, :site, :nambu)
    @test filter(≠(:nambu), m) == OIDToTuple(:scope, :site)
    @test OIDToTuple(Index{PID, DID{Int}}) == OIDToTuple(:site, :nambu)
    @test OIDToTuple(OID{Index{CPID{Int}, DID{Int}}}) == OIDToTuple(:scope, :site, :nambu)
    @test OIDToTuple(Hilbert{DFock, CPID{Int}}) == OIDToTuple(:scope, :site, :nambu)
    @test m(index) == ('S', 4, 1) == m(oid)
    @test filter(≠(:scope), m)(index) == (4, 1)
    @test filter(≠(:nambu), m)(index) == ('S', 4)
    @test filter(∉((:site, :nambu)), m)(index) == ('S',)
end

@testset "Table" begin
    @test contentnames(Table) == (:by, :contents)

    by = filter(≠(:nambu), OIDToTuple(Index{PID, DID{Int}}))

    table = Table([Index(PID(1), DID(1)), Index(PID(1), DID(2))], by)
    @test empty(table) == Table{Index{PID, DID{Int}}}(by)
    @test table[Index(PID(1), DID(1))]==1 && table[Index(PID(1), DID(2))]==1

    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    inds1 = (Index(PID(1), iid) for iid in DFock(2))|>collect
    inds2 = (Index(PID(2), iid) for iid in DFock(2))|>collect
    @test Table(hilbert, by) == Table([inds1; inds2], by)
    @test Table(hilbert, by) == union(Table(inds1, by), Table(inds2, by))

    opt = Operator(1.0im, ID(
        OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(1.0, 0.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
        ))
    @test sequence(opt, table) == (1, 1)
    @test haskey(table, opt.id) == (true, true)

    table = Table(hilbert)
    @test table == Table([inds1; inds2])
    @test reset!(empty(table), [inds1; inds2]) == table
    @test reset!(empty(table), hilbert) == table
end

@testset "Subscript" begin
    sub = Subscript(4)
    @test contentnames(typeof(sub)) == (:contents, :rep, :constraint)
    @test getcontent(sub, :contents) == sub.pattern
    @test sub==deepcopy(sub) && isequal(sub, deepcopy(sub))
    @test string(sub) == "[* * * *]"
    @test rank(sub) == rank(typeof(sub)) == 4
    @test match(sub, (2, 2, 2, 2))

    sub = Subscript((1, 2, 3, 4))
    @test string(sub) == "[1 2 3 4]"
    @test sub == subscript"[1 2 3 4]" == subscript"[1 2 3 4]; false"
    @test rank(sub) == 4
    @test match(sub, (1, 3, 3, 4))

    sub = Subscript((1, 2, 2, 1), true)
    @test string(sub) == "[1 2 2 1]"
    @test sub == subscript"[1 2 2 1]; true"
    @test rank(sub) == 4
    @test match(sub, (1, 2, 2, 1)) && !match(sub, (1, 1, 2, 1))

    sub = subscript"[x₁ x₂ x₁ x₂]"
    @test string(sub) == "[x₁ x₂ x₁ x₂]"
    @test rank(sub) == 4
    @test match(sub, (1, 2, 1, 2)) && !match(sub, (1, 2, 2, 1))

    sub = subscript"[x₁ 4 4 x₂](x₁ < x₂)"
    @test string(sub) == "[x₁ 4 4 x₂](x₁ < x₂)"
    @test rank(sub) == rank(typeof(sub)) == 4
    @test match(sub, (1, 4, 4, 2)) && !match(sub, (2, 4, 4, 1))
end

@testset "Subscripts" begin
    subscripts = Subscripts((nambu=subscript"[x y](x > y)",), (nambu=subscript"[x x y y](x < y)",))
    @test string(subscripts) == "nambu[x y](x > y) × nambu[x x y y](x < y)"
    @test repr(subscripts, 1:2, :nambu) == "[x y](x > y)×[x x y y](x < y)"
    @test repr(subscripts, 1, :nambu) == "[x y](x > y)"
    @test repr(subscripts, 2, :nambu) == "[x x y y](x < y)"

    @test rank(subscripts) == rank(typeof(subscripts)) == 6
    @test rank(subscripts, 1) == rank(typeof(subscripts), 1) == 2
    @test rank(subscripts, 2) == rank(typeof(subscripts), 2) == 4
    @test match(subscripts, (DID(2), DID(1), DID(2), DID(2), DID(3), DID(3)))
    @test match(subscripts, DID(2)⊗DID(1)⊗DID(2)⊗DID(2)⊗DID(3)⊗DID(3))
    @test !match(subscripts, DID(1)⊗DID(2)⊗DID(2)⊗DID(2)⊗DID(3)⊗DID(3))
    @test !match(subscripts, DID(2)⊗DID(1)⊗DID(3)⊗DID(3)⊗DID(2)⊗DID(2))
    @test subscripts == Subscripts((nambu=subscript"[x y](x > y)",))*Subscripts((nambu=subscript"[x x y y](x < y)",))

    subscriptsid = SubscriptsID(subscripts)
    @test subscriptsid == SubscriptsID((1:2=>("[x y](x > y)",), 3:6=>("[x x y y](x < y)",)))
    @test idtype(subscripts) == idtype(typeof(subscripts)) == SubscriptsID{NTuple{2, Pair{UnitRange{Int}, Tuple{String}}}}
end

@testset "couplings" begin
    @test parameternames(Coupling) == (:value, :cid, :subscripts, :subscriptsid)
    @test contentnames(Coupling) == (:value, :id, :subscripts)
    @test isparameterbound(Coupling, :cid, Tuple{DID{Int}}) == false
    @test isparameterbound(Coupling, :cid, ID{DID{Int}}) == true
    @test isparameterbound(Coupling, :subscripts, Subscripts{Tuple{NamedTuple{(:nambu,), Tuple{Subscript{Tuple{Int}, typeof(noconstrain)}}}}}) == false
    @test isparameterbound(Coupling, :subscripts, Subscripts) == true
    @test isparameterbound(Coupling, :subscriptsid, SubscriptsID{Tuple{}}) == false
    @test isparameterbound(Coupling, :subscriptsid, SubscriptsID) == true

    tc = DCoupling(2.0, (2,))
    @test id(tc) == ID(CompositeIID(tc.cid), SubscriptsID(tc.subscripts))
    @test rank(tc) == rank(typeof(tc)) == 1
    @test tc == Coupling(2.0, id(tc), tc.subscripts)
    @test ID{SimpleIID}(tc) == tc.cid
    @test Subscripts(tc) == tc.subscripts

    point = Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>DFock(2))
    @test couplingcenters(tc, point, Val(:Mu)) == (1,)
    @test couplingpoints(tc, point, Val(:Mu)) == (point,)
    @test couplinginternals(tc, point, hilbert, Val(:Mu)) == (DFock(2),)

    tc₁ = DCoupling(1.5, (1, 2))
    tc₂ = DCoupling(2.0, subscript"[a b](a < b)")
    @test tc₁*tc₂ == Coupling(3.0, ID(DID(1), DID(2), DID(:a), DID(:b)), Subscripts((nambu=Subscript((1, 2)),), (nambu=subscript"[a b](a < b)",)))

    ex = expand(tc₁, point, hilbert, Val(:info))
    @test eltype(ex) == eltype(typeof(ex)) == Tuple{Float64, NTuple{2, OID{Index{CPID{Int}, DID{Int}}, SVector{2, Float64}}}}
    @test collect(ex) == [
        (1.5, ID(OID(Index(CPID(1, 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
                 OID(Index(CPID(1, 1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0))
                 ))
        ]

    ex = expand(tc₂, point, hilbert, Val(:info))
    @test eltype(ex) == eltype(typeof(ex)) == Tuple{Float64, NTuple{2, OID{Index{CPID{Int}, DID{Int}}, SVector{2, Float64}}}}
    @test collect(ex) == [
        (2.0, ID(OID(Index(CPID(1, 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
                    OID(Index(CPID(1, 1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0))
                    ))
        ]

    @test @couplings(DCoupling(1.0, (1, 1))) == Couplings(DCoupling(1.0, (1, 1)))
    @test @couplings(DCoupling(1.0, (1,))+DCoupling(1.0, (2,))) == DCoupling(1.0, (1,))+DCoupling(1.0, (2,))
end

@testset "TermFunction" begin
    ta = TermAmplitude()
    @test ta() == 1

    ta = TermAmplitude(x->x+3.0)
    @test ta(1.0) == 4.0

    tcs = DCoupling(1.0, (1, 1)) + DCoupling(2.0, (2, 2))
    termcouplings = TermCouplings(tcs)
    @test termcouplings == deepcopy(TermCouplings(tcs))
    @test isequal(termcouplings, deepcopy(TermCouplings(tcs)))
    @test valtype(termcouplings) == valtype(typeof(termcouplings)) == typeof(tcs)
    @test termcouplings() == tcs

    fx = i -> i%2==1 ? (DCoupling(1.0, (1, 1))+DCoupling(1.0, (2, 2))) : (DCoupling(1.0, (2, 1))+DCoupling(1.0, (1, 2)))
    termcouplings = TermCouplings(fx)
    @test termcouplings == TermCouplings{typejoin(typeof(fx(1)), typeof(fx(2)))}(fx)
    @test valtype(termcouplings) == valtype(typeof(termcouplings)) == typejoin(typeof(fx(1)), typeof(fx(2)))
    @test termcouplings(1) == fx(1)
    @test termcouplings(2) == fx(2)

    @test ismodulatable(TermModulate{Val{true}}) == ismodulatable(TermModulate{<:Function}) == true
    @test ismodulatable(TermModulate{Val{false}}) == false
    termmodulate = TermModulate(:t)
    @test termmodulate(t=1) == 1
    @test isnothing(termmodulate(mu=1))
    @test isnothing(termmodulate())
    @test ismodulatable(termmodulate) == true

    termmodulate = TermModulate(:t, t->t*2.0)
    @test ismodulatable(termmodulate) == true
    @test termmodulate(2) == 4
end

@testset "Term" begin
    point = Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>DFock(2))
    term = Term{:Mu}(:mu, 1.5, 0, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=(bond->3), modulate=false)
    @test term|>kind == term|>typeof|>kind == :Mu
    @test term|>id == term|>typeof|>id == :mu
    @test term|>valtype == term|>typeof|>valtype == Float
    @test term|>rank == term|>typeof|>rank == 2
    @test term|>abbr == term|>typeof|>abbr == :mu
    @test term|>ismodulatable == term|>typeof|>ismodulatable == false
    @test term|>ishermitian == term|>typeof|>ishermitian == true
    @test term == deepcopy(term)
    @test isequal(term, deepcopy(term))
    @test string(term) == "Mu{2}(id=mu, value=1.5, bondkind=0, factor=1.0)"
    @test repr(term, point, hilbert) == "mu: 4.5 ph(2, 1)"
    @test +term == term
    @test -term == term*(-1) == replace(term, factor=-term.factor)
    @test 2*term == term*2 == replace(term, factor=term.factor*2)
    @test term/2 == term*(1/2) == replace(term, factor=term.factor/2)

    p1 = Point(PID(1), (0.0, 0.0), (0.0, 0.0))
    p2 = Point(PID(2), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    term = Term{:Mu}(:mu, 1.5, 0,
        couplings=bond->bond.pid.site%2==0 ? Couplings(DCoupling(1.0, (2, 2))) : Couplings(DCoupling(1.0, (1, 1))),
        amplitude=bond->3,
        modulate=true
    )
    @test term|>ismodulatable == term|>typeof|>ismodulatable == true
    @test repr(term, p1, hilbert) == "mu: 4.5 ph(1, 1)"
    @test repr(term, p2, hilbert) == "mu: 4.5 ph(2, 2)"
    @test one(term) == replace(term, value=1.0)
    @test zero(term) == replace(term, value=0.0)
    @test term.modulate(mu=4.0) == 4.0
    @test isnothing(term.modulate(t=1.0))
    @test update!(term, mu=4.25) == replace(term, value=4.25)
    @test term.value == 4.25
end

@testset "expand" begin
    point = Point(PID(1), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>DFock(2))
    term = Term{:Mu}(:mu, 1.5, 0, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(
        Operator(+2.25, ID(
            OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
            OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
            ))
        )
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators*2

    bond = Bond(1, Point(PID(2), (1.5, 1.5), (1.0, 1.0)), Point(PID(1), (0.5, 0.5), (0.0, 0.0)))
    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    term = Term{:Hp}(:t, 1.5, 1, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(
        Operator(4.5, ID(
            OID(Index(PID(1), DID(2)), SVector(0.5, 0.5), SVector(0.0, 0.0)),
            OID(Index(PID(2), DID(1)), SVector(1.5, 1.5), SVector(1.0, 1.0))
            ))
        )
    @test expand(term, bond, hilbert, half=true) == operators
    @test expand(term, bond, hilbert, half=false) == operators+operators'

    lattice = Lattice("Tuanzi", [Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))], vectors=[[1.0, 0.0]], neighbors=1)
    bonds = Bonds(lattice)
    hilbert = Hilbert(pid=>DFock(2) for pid in lattice.pids)
    term = Term{:Mu}(:mu, 1.5, 0, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(
        Operator(+2.25, ID(
            OID(Index(CPID(1, 1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
            OID(Index(CPID(1, 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
            ))
        )
    @test expand(term, bonds, hilbert, half=true) == operators
    @test expand(term, bonds, hilbert, half=false) == operators*2
end

@testset "LaTeX" begin
    @test latexname(OID) == :OID
    @test latexformat(OID{<:Index{<:AbstractPID, <:DID}}) == LaTeX{(:nambu,), (:site,)}('d')

    latex = LaTeX{(:nambu,), (:site,)}('c', vectors=(SVector(1.0, 0.0), SVector(0.0, 1.0)))
    @test superscript(latex|>typeof) == (:nambu,)
    @test subscript(latex|>typeof) == (:site,)

    oid = OID(Index(CPID('d', 1), DID(2)), SVector(0.0, 0.0), SVector(1.0, 0.0))
    @test script(Val(:rcoord), oid) == "[0.0, 0.0]"
    @test script(Val(:icoord), oid) == "[1.0, 0.0]"
    @test script(Val(:integralicoord), oid; vectors=get(latex.options, :vectors, nothing)) == "[1, 0]"
    @test script(Val(:BD), oid, latex) == 'c'
    @test script(Val(:SP), oid, latex) == ("\\dagger",)
    @test script(Val(:SB), oid, latex) == (1,)
    @test repr(oid, latex) == "c^{\\dagger}_{1}"

    opt = Operator(1.0, ID(
        OID(Index(CPID('d', 2), DID(2)), SVector(1.0, 0.0), SVector(2.0, 0.0)),
        OID(Index(CPID('d', 2), DID(2)), SVector(1.0, 0.0), SVector(2.0, 0.0))
        ))
    @test repr(opt) == "(d^{\\dagger}_{2})^2"

    opt = Operator(1.0, ID(
        OID(Index(CPID('d', 2), DID(2)), SVector(1.0, 0.0), SVector(2.0, 0.0)),
        OID(Index(CPID('d', 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
        ))
    @test repr(opt) == "d^{\\dagger}_{2}d^{}_{1}"

    latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:site,)}('c'))
    @test repr(opt) == "c^{\\dagger}_{2}c^{}_{1}"

    latexformat(OID{<:Index{<:AbstractPID, <:DID}},  LaTeX{(:nambu,), (:site,)}('d'))
    @test repr(opt) == "d^{\\dagger}_{2}d^{}_{1}"

    latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:rcoord,)}('d'))
    @test repr(opt) == "d^{\\dagger}_{[1.0, 0.0]}d^{}_{[0.0, 0.0]}"

    latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:integralicoord,)}('d', vectors=(SVector(1.0, 0.0), SVector(0.0, 1.0))))
    @test repr(opt) == "d^{\\dagger}_{[2, 0]}d^{}_{[0, 0]}"

    opts = Operators(
            Operator(1.0-1.0im, ID(
                OID(Index(PID(2), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
                OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
                )),
            Operator(-1.0, ID(
                OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
                OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
                ))
            )
    latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:site,)}('c'))
    @test repr(opts) == "-c^{\\dagger}_{1}c^{}_{1}+(1.0-1.0im)c^{\\dagger}_{2}c^{}_{1}"

    latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:site,)}('d'))
    io = IOBuffer()
    show(io, MIME"text/latex"(), opt)
    @test String(take!(io)) == "\$d^{\\dagger}_{2}d^{}_{1}\$"
    show(io, MIME"text/latex"(), opts)
    @test String(take!(io)) == "\$-d^{\\dagger}_{1}d^{}_{1}+(1.0-1.0im)d^{\\dagger}_{2}d^{}_{1}\$"
end

@testset "Boundary" begin
    opt = Operator(4.5, ID(
        OID(Index(PID(1), DID(2)), SVector(0.5, 0.5), SVector(0.0, 0.0)),
        OID(Index(PID(2), DID(1)), SVector(1.5, 1.5), SVector(1.0, 1.0))
        ))
    bound = Boundary{(:θ₁, :θ₂)}([0.1, 0.2], [[1.0, 0.0], [0.0, 1.0]])
    @test keys(bound) == keys(typeof(bound)) == (:θ₁, :θ₂)
    @test bound==deepcopy(bound) && isequal(bound, deepcopy(bound))
    @test angle(bound, opt) ≈ angle(opt.id, bound.vectors, bound.values) ≈ 0.6pi
    @test bound(opt) ≈ twist(opt, bound.vectors, bound.values) ≈ replace(opt, 4.5*exp(2im*pi*0.3))
    update!(bound, θ₁=0.3)
    @test bound(opt) ≈ replace(opt, 4.5*exp(2im*pi*0.5))

    @test angle(plain, opt) == 0
    @test plain(opt) == opt
    @test update!(plain) == plain
end


text:
using Glutils
using ModernGL
using GLFW
# matrices utils
import Glutils: Mat
# Create the window
window = GLFW.CreateWindow(800, 600, "Basic Example")
# set the current OpenGL context
GLFW.MakeContextCurrent(window)
# create a vertex buffer
vbo = arraybuffer(Cfloat[
    -0.5, -0.5, 0.0,
     0.5, -0.5, 0.0,
     0.0,  0.5, 0.0
], g"static draw")
# create a vertex array
vao = vertexarray()
vertexattrib(Cfloat, 3,
    (size = 3, offset = 0), # 0
)
# creates the shader program
# note: here we use the --- <type> pattern
# to indicate the beginning of each shader
program = program"""
--- vertex
#version 330 core
layout(location = 0) in vec3 pos;

uniform mat4 rotation;

void main() {
    gl_Position = rotation * vec4(pos.x, pos.y, pos.z, 1.0);
}

--- fragment
#version 330 core
out vec4 color;
uniform vec4 our_color;

void main() {
    color = our_color;
}
"""
# get the location of our_color
# and rotation uniforms
@uniforms our_color, rotation = program
t = 0.0
# background color
glClearColor(0.0, 0.0, 0.0, 1.0)
# main loop
while true
    # Break the loop when the window closes
    GLFW.WindowShouldClose(window) && break
    # Wait for events
    GLFW.PollEvents()
    # Clear the screen
    glClear(g"color buffer bit")
    # set the uniform variables
    @uniforms! begin
        # compute the color based on t
        r, g, b = [sin(t), cos(t), sin(t)] ./ 2 .+ 1
        # the 'foo := bar' syntax sets the uniform value of foo
        our_color := Cfloat[r, g, b, 1.0]
        # rotate in the Z-axis
        rotation := Mat.rotationZ(t, Cfloat)
    end
    # Render here
    glBindVertexArray(vao)
    glUseProgram(program)
    glDrawArrays(g"TRIANGLES", 0, 3)
    # Swap the buffers
    GLFW.SwapBuffers(window)
    global t += 0.0001;
end

GLFW.DestroyWindow(window)

ext: jl


text:
using MinimaxCalibratedEBayes
using Test

isa(MarginalDensityTarget, Type{T} where T<: LinearEBayesTarget)


text:
# Tests for invertible neural network module
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020
using InvertibleNetworks, Test

test_suite = "all"   # "all", "layers" or "networks"

layers = ["test_utils/test_objectives.jl",
          "test_utils/test_nnlib_convolution.jl",
          "test_utils/test_activations.jl",
          "test_utils/test_squeeze.jl",
          # Layers
          "test_layers/test_residual_block.jl",
          "test_layers/test_flux_block.jl",
          "test_layers/test_householder_convolution.jl",
          "test_layers/test_coupling_layer_basic.jl",
          "test_layers/test_coupling_layer_irim.jl",
          "test_layers/test_coupling_layer_glow.jl",
          "test_layers/test_coupling_layer_hint.jl",
          "test_layers/test_coupling_layer_slim.jl",
          "test_layers/test_coupling_layer_slim_learned.jl",
          "test_layers/test_conditional_layer_hint.jl",
          "test_layers/test_conditional_layer_slim.jl",
          "test_layers/test_conditional_res_block.jl",
          "test_layers/test_hyperbolic_layer.jl",
          "test_layers/test_actnorm.jl"]

networks = ["test_networks/test_unrolled_loop.jl",
            "test_networks/test_generator.jl",
            "test_networks/test_glow.jl",
            "test_networks/test_hyperbolic_network.jl",
            "test_networks/test_conditional_hint_network.jl"]

if test_suite == "all" || test_suite == "layers"
    for t=layers
        @testset  "Test $t" begin
            include(t)
        end
    end
end

# Networks
if test_suite == "all" || test_suite == "networks"
    for t=networks
        @testset  "Test $t" begin
            include(t)
        end
    end
end


text:
import LinearAlgebra: dot

struct Auto
    μ::Float64
    τ::Int
    M_Δ::Float64
    α_init::Float64
    α::IdDict
    h::IdDict
    n::IdDict
end

Auto(μ, α_init, τ=10000, M_Δ=1.0) = Auto(μ, τ, M_Δ, α_init, IdDict(), IdDict(), IdDict())

# Assumes ϕ and θ are the same size!!!!
function Flux.Optimise.update!(opt::Auto, θ::AbstractArray{F}, ϕ, δ, z, num_demons, num_rows_per_demon) where {F<:AbstractFloat}
    α = get!(()->zero(θ) .+ F(opt.α_init), opt.α, θ)::typeof(θ)
    h = get!(()->zero(θ), opt.h, θ)::typeof(θ)
    n = get!(()->zero(θ) .+ 1, opt.n, θ)::typeof(θ)

    # @show length(findall(!iszero, α .- opt.α_init))
    # @show size(α)

    M_Δ, μ, τ = opt.M_Δ, opt.μ, opt.τ

    δϕ = δ .* ϕ
    # println(size(δϕ))
    abs_ϕ = abs.(ϕ)
    hδϕ = (δϕ).*h
    n .+= (1/τ) * α.*abs_ϕ.*(abs.(hδϕ) - n)

    ϕ_nz_idx = ϕ .!= 0.0
    n_ϕ_nz = @view n[ϕ_nz_idx]
    hδϕ_ϕ_nz = @view hδϕ[ϕ_nz_idx]
    α_ϕ_nz = @view α[ϕ_nz_idx]

    # Δβ = sign.(hδϕ_ϕ_nz) .* min.(abs.(hδϕ_ϕ_nz./n_ϕ_nz), M_Δ)
    # α_ϕ_nz .= min.(α_ϕ_nz .* exp.(μ * Δβ), 1.0 ./(abs_ϕ[ϕ_nz_idx]))
    Δβ = clamp.(hδϕ_ϕ_nz./n_ϕ_nz, -1, 1)
    α_ϕ_nz .= clamp.(α_ϕ_nz .* exp.(μ * Δβ), 1e-6, 1 ./ (abs_ϕ[ϕ_nz_idx]))

    for d ∈ 1:num_demons
        rr = ((d-1)*num_rows_per_demon + 1):(d*num_rows_per_demon)
        α_d = @view α[rr, :]
        z_d = @view z[rr, :]
        if dot(α_d, z_d) > 1
            z_nz_idx = z_d .!= 0.0
            α_z_nz = @view α_d[z_nz_idx]
            α_z_nz .= min.(α_z_nz, 1 ./sum(abs.(z_d)))
        # else
        #     # @show "Not Triggered for $(d)"
        end
    end
    θ .+= α.*δϕ
    # @show "after", sum(θ)
    h .= h.*(1 .- α.*abs_ϕ) + α.*δϕ

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
int_rules_4_7_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.7.5*Inert*trig*functions =#
    @apply_utils Antiderivative((~(a') * (~F)(~(c') + ~(d') * ~x) ^ ~p) ^ ~n, ~x) => With([v = ActivateTrig((~F)(~c + ~d * ~x))], (((~a) ^ IntPart(~n) * (v / NonfreeFactors(v, ~x)) ^ (~p * IntPart(~n)) * (~a * v ^ ~p) ^ FracPart(~n)) / NonfreeFactors(v, ~x) ^ (~p * FracPart(~n))) * Antiderivative(NonfreeFactors(v, ~x) ^ (~n * ~p), ~x)) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (InertTrigQ(~F) && (Not(IntegerQ(~n)) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(a') * (~(b') * (~F)(~(c') + ~(d') * ~x)) ^ ~p) ^ ~(n'), ~x) => With([v = ActivateTrig((~F)(~c + ~d * ~x))], (((~a) ^ IntPart(~n) * (~a * (~b * v) ^ ~p) ^ FracPart(~n)) / (~b * v) ^ (~p * FracPart(~n))) * Antiderivative((~b * v) ^ (~n * ~p), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (InertTrigQ(~F) && (Not(IntegerQ(~n)) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, cos) || EqQ(~F, cos))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, sin) || EqQ(~F, sin))
    @apply_utils Antiderivative(~u * cosh(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * sinh(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, cot) || EqQ(~F, cot))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, tan) || EqQ(~F, tan))
    @apply_utils Antiderivative(~u * coth(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * tanh(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ 2, ~x) => With([d = FreeFactors(tan(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, tan(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, tan(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(tan(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && (NonsumQ(~u) && (EqQ(~F, sec) || EqQ(~F, sec)))
    @apply_utils Antiderivative(~u / cos(~(c') * (~(a') + ~(b') * ~x)) ^ 2, ~x) => With([d = FreeFactors(tan(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, tan(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, tan(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(tan(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && NonsumQ(~u)
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ 2, ~x) => With([d = FreeFactors(cot(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, cot(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cot(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cot(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && (NonsumQ(~u) && (EqQ(~F, csc) || EqQ(~F, csc)))
    @apply_utils Antiderivative(~u / sin(~(c') * (~(a') + ~(b') * ~x)) ^ 2, ~x) => With([d = FreeFactors(cot(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, cot(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cot(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cot(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && NonsumQ(~u)
    @apply_utils Antiderivative(~u * sech(~(c') * (~(a') + ~(b') * ~x)) ^ 2, ~x) => With([d = FreeFactors(tanh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, tanh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, tanh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(tanh(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && NonsumQ(~u)
    @apply_utils Antiderivative(~u * csch(~(c') * (~(a') + ~(b') * ~x)) ^ 2, ~x) => With([d = FreeFactors(coth(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, coth(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, coth(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(coth(~c * (~a + ~b * ~x)) / d, ~u, ~x, True)) <-- FreeQ([~a, ~b, ~c], ~x) && NonsumQ(~u)
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~(n'), ~x) => With([d = FreeFactors(tan(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor(1 / ((~x) ^ ~n * (1 + d ^ 2 * (~x) ^ 2)), tan(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, tan(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(tan(~c * (~a + ~b * ~x)) / d, ~u, ~x, True) && TryPureTanSubst(ActivateTrig(~u) * cot(~c * (~a + ~b * ~x)) ^ ~n, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ(~n) && (EqQ(~F, cot) || EqQ(~F, cot)))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~(n'), ~x) => With([d = FreeFactors(cot(~c * (~a + ~b * ~x)), ~x)], (-1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor(1 / ((~x) ^ ~n * (1 + d ^ 2 * (~x) ^ 2)), cot(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cot(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cot(~c * (~a + ~b * ~x)) / d, ~u, ~x, True) && TryPureTanSubst(ActivateTrig(~u) * tan(~c * (~a + ~b * ~x)) ^ ~n, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ(~n) && (EqQ(~F, tan) || EqQ(~F, tan)))
    @apply_utils Antiderivative(~u * coth(~(c') * (~(a') + ~(b') * ~x)) ^ ~(n'), ~x) => With([d = FreeFactors(tanh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor(1 / ((~x) ^ ~n * (1 - d ^ 2 * (~x) ^ 2)), tanh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, tanh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(tanh(~c * (~a + ~b * ~x)) / d, ~u, ~x, True) && TryPureTanSubst(ActivateTrig(~u) * coth(~c * (~a + ~b * ~x)) ^ ~n, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative(~u * tanh(~(c') * (~(a') + ~(b') * ~x)) ^ ~(n'), ~x) => With([d = FreeFactors(coth(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor(1 / ((~x) ^ ~n * (1 - d ^ 2 * (~x) ^ 2)), coth(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, coth(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(coth(~c * (~a + ~b * ~x)) / d, ~u, ~x, True) && TryPureTanSubst(ActivateTrig(~u) * tanh(~c * (~a + ~b * ~x)) ^ ~n, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && IntegerQ(~n)
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], ShowStep("", "Antiderivative(F(cot(a+b*x)),x)", "-1/b*Subst(Antiderivative(F(x)/(1+x^2),x),x,cot(a+b*x))", Hold( With([d = FreeFactors(cot(v), x)], Dist(-d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), cot(v)/d, u, x), x), x, cot(v)/d), x)))) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(cot(v), x), u, x, True) && TryPureTanSubst(ActivateTrig(u), x)) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], With([d = FreeFactors(cot(v), x)], Dist(-d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), cot(v)/d, u, x), x), x, cot(v)/d), x)) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(cot(v), x), u, x, True) && TryPureTanSubst(ActivateTrig(u), x)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], ShowStep("", "Antiderivative(F(tan(a+b*x)),x)", "1/b*Subst(Antiderivative(F(x)/(1+x^2),x),x,tan(a+b*x))", Hold( With([d = FreeFactors(tan(v), x)], Dist(d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), tan(v)/d, u, x), x), x, tan(v)/d), x)))) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(tan(v), x), u, x, True) && TryPureTanSubst(ActivateTrig(u), x)) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], With([d = FreeFactors(tan(v), x)], Dist(d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), tan(v)/d, u, x), x), x, tan(v)/d), x)) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(tan(v), x), u, x, True) && TryPureTanSubst(ActivateTrig(u), x)))
    @apply_utils Antiderivative((~F)(~(a') + ~(b') * ~x) ^ ~(p') * (~G)(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce(ActivateTrig((~F)(~a + ~b * ~x) ^ ~p * (~G)(~c + ~d * ~x) ^ ~q), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && ((EqQ(~F, sin) || EqQ(~F, cos)) && ((EqQ(~G, sin) || EqQ(~G, cos)) && (IGtQ(~p, 0) && IGtQ(~q, 0))))
    @apply_utils Antiderivative((~F)(~(a') + ~(b') * ~x) ^ ~(p') * (~G)(~(c') + ~(d') * ~x) ^ ~(q') * (~H)(~(e') + ~(f') * ~x) ^ ~(r'), ~x) => Antiderivative(ExpandTrigReduce(ActivateTrig((~F)(~a + ~b * ~x) ^ ~p * (~G)(~c + ~d * ~x) ^ ~q * (~H)(~e + ~f * ~x) ^ ~r), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && ((EqQ(~F, sin) || EqQ(~F, cos)) && ((EqQ(~G, sin) || EqQ(~G, cos)) && ((EqQ(~H, sin) || EqQ(~H, cos)) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && IGtQ(~r, 0))))))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, cos) || EqQ(~F, cos))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, sin) || EqQ(~F, sin))
    @apply_utils Antiderivative(~u * cosh(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * sinh(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor(1, cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, cot) || EqQ(~F, cot))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~F, tan) || EqQ(~F, tan))
    @apply_utils Antiderivative(~u * coth(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * tanh(~(c') * (~(a') + ~(b') * ~x)), ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c)) * Subst(Antiderivative(SubstFor(1 / ~x, cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor((1 - d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2), sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, cos) || EqQ(~F, cos))))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor((1 - d ^ 2 * (~x) ^ 2) ^ ((-(~n) - 1) / 2), sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, sec) || EqQ(~F, sec))))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor((1 - d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2), cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, sin) || EqQ(~F, sin))))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-d / (~b * ~c)) * Subst(Antiderivative(SubstFor((1 - d ^ 2 * (~x) ^ 2) ^ ((-(~n) - 1) / 2), cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, csc) || EqQ(~F, csc))))
    @apply_utils Antiderivative(~u * cosh(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor((1 + d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2), sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && NonsumQ(~u))
    @apply_utils Antiderivative(~u * sech(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor((1 + d ^ 2 * (~x) ^ 2) ^ ((-(~n) - 1) / 2), sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && NonsumQ(~u))
    @apply_utils Antiderivative(~u * sinh(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor((-1 + d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2), cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && NonsumQ(~u))
    @apply_utils Antiderivative(~u * csch(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (d / (~b * ~c)) * Subst(Antiderivative(SubstFor((-1 + d ^ 2 * (~x) ^ 2) ^ ((-(~n) - 1) / 2), cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && NonsumQ(~u))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor((1 - d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2) / (~x) ^ ~n, sin(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sin(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, cot) || EqQ(~F, cot))))
    @apply_utils Antiderivative(~u * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], (-1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor((1 - d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2) / (~x) ^ ~n, cos(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cos(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, tan) || EqQ(~F, tan))))
    @apply_utils Antiderivative(~u * coth(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(sinh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor((1 + d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2) / (~x) ^ ~n, sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, sinh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(sinh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && NonsumQ(~u))
    @apply_utils Antiderivative(~u * tanh(~(c') * (~(a') + ~(b') * ~x)) ^ ~n, ~x) => With([d = FreeFactors(cosh(~c * (~a + ~b * ~x)), ~x)], (1 / (~b * ~c * d ^ (~n - 1))) * Subst(Antiderivative(SubstFor((-1 + d ^ 2 * (~x) ^ 2) ^ ((~n - 1) / 2) / (~x) ^ ~n, cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x), ~x), ~x, cosh(~c * (~a + ~b * ~x)) / d) <-- FunctionOfQ(cosh(~c * (~a + ~b * ~x)) / d, ~u, ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ((~n - 1) / 2) && NonsumQ(~u))
    @apply_utils Antiderivative(~u * (~v + ~(d') * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~(n')), ~x) => With([e = FreeFactors(sin(~c * (~a + ~b * ~x)), ~x)], Antiderivative(ActivateTrig(~u * ~v), ~x) + ~d * Antiderivative(ActivateTrig(~u) * cos(~c * (~a + ~b * ~x)) ^ ~n, ~x) <-- FunctionOfQ(sin(~c * (~a + ~b * ~x)) / e, ~u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (Not(FreeQ(~v, ~x)) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, cos) || EqQ(~F, cos)))))
    @apply_utils Antiderivative(~u * (~v + ~(d') * (~F)(~(c') * (~(a') + ~(b') * ~x)) ^ ~(n')), ~x) => With([e = FreeFactors(cos(~c * (~a + ~b * ~x)), ~x)], Antiderivative(ActivateTrig(~u * ~v), ~x) + ~d * Antiderivative(ActivateTrig(~u) * sin(~c * (~a + ~b * ~x)) ^ ~n, ~x) <-- FunctionOfQ(cos(~c * (~a + ~b * ~x)) / e, ~u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (Not(FreeQ(~v, ~x)) && (IntegerQ((~n - 1) / 2) && (NonsumQ(~u) && (EqQ(~F, sin) || EqQ(~F, sin)))))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], ShowStep("", "Antiderivative(F(sin(a+b*x))*cos(a+b*x),x)", "Subst(Antiderivative(F(x),x),x,sin(a+b*x))/b", Hold( With([d = FreeFactors(sin(v), x)], Dist(d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1, sin(v)/d, u/cos(v), x), x), x, sin(v)/d), x)))) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(sin(v), x), u/cos(v), x)) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], With([d = FreeFactors(sin(v), x)], Dist(d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1, sin(v)/d, u/cos(v), x), x), x, sin(v)/d), x)) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(sin(v), x), u/cos(v), x)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], ShowStep("", "Antiderivative(F(cos(a+b*x))*sin(a+b*x),x)", "-Subst(Antiderivative(F(x),x),x,cos(a+b*x))/b", Hold( With([d = FreeFactors(cos(v), x)], Dist(-d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1, cos(v)/d, u/sin(v), x), x), x, cos(v)/d), x)))) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(cos(v), x), u/sin(v), x)) <-- SimplifyFlag, Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], With([d = FreeFactors(cos(v), x)], Dist(-d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1, cos(v)/d, u/sin(v), x), x), x, cos(v)/d), x)) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(cos(v), x), u/sin(v), x)))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * cos(~(d') + ~(e') * ~x) ^ 2 + ~(c') * sin(~(d') + ~(e') * ~x) ^ 2) ^ ~(p'), ~x) => (~a + ~c) ^ ~p * Antiderivative(ActivateTrig(~u), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && EqQ(~b - ~c, 0)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ 2 + ~(c') * sec(~(d') + ~(e') * ~x) ^ 2) ^ ~(p'), ~x) => (~a + ~c) ^ ~p * Antiderivative(ActivateTrig(~u), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && EqQ(~b + ~c, 0)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ 2 + ~(c') * csc(~(d') + ~(e') * ~x) ^ 2) ^ ~(p'), ~x) => (~a + ~c) ^ ~p * Antiderivative(ActivateTrig(~u), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && EqQ(~b + ~c, 0)
    @apply_utils Antiderivative(~u / ~y, ~x) => With([q = DerivativeDivides(ActivateTrig(~y), ActivateTrig(~u), ~x)], q * log(RemoveContent(ActivateTrig(~y), ~x)) <-- Not(FalseQ(q))) <-- Not(InertTrigFreeQ(~u))
    @apply_utils Antiderivative(~u / (~y * ~w), ~x) => With([q = DerivativeDivides(ActivateTrig(~y * ~w), ActivateTrig(~u), ~x)], q * log(RemoveContent(ActivateTrig(~y * ~w), ~x)) <-- Not(FalseQ(q))) <-- Not(InertTrigFreeQ(~u))
    @apply_utils Antiderivative(~u * (~y) ^ ~(m'), ~x) => With([q = DerivativeDivides(ActivateTrig(~y), ActivateTrig(~u), ~x)], (q * ActivateTrig((~y) ^ (~m + 1))) / (~m + 1) <-- Not(FalseQ(q))) <-- FreeQ(~m, ~x) && (NeQ(~m, -1) && Not(InertTrigFreeQ(~u)))
    @apply_utils Antiderivative(~u * (~y) ^ ~(m') * (~z) ^ ~(n'), ~x) => With([q = DerivativeDivides(ActivateTrig(~y * ~z), ActivateTrig(~u * (~z) ^ (~n - ~m)), ~x)], (q * ActivateTrig((~y) ^ (~m + 1) * (~z) ^ (~m + 1))) / (~m + 1) <-- Not(FalseQ(q))) <-- FreeQ([~m, ~n], ~x) && (NeQ(~m, -1) && Not(InertTrigFreeQ(~u)))
    @apply_utils Antiderivative(~(u') * (~(a') * (~F)(~(c') + ~(d') * ~x) ^ ~p) ^ ~n, ~x) => With([v = ActivateTrig((~F)(~c + ~d * ~x))], (((~a) ^ IntPart(~n) * (v / NonfreeFactors(v, ~x)) ^ (~p * IntPart(~n)) * (~a * v ^ ~p) ^ FracPart(~n)) / NonfreeFactors(v, ~x) ^ (~p * FracPart(~n))) * Antiderivative(ActivateTrig(~u) * NonfreeFactors(v, ~x) ^ (~n * ~p), ~x)) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (InertTrigQ(~F) && (Not(IntegerQ(~n)) && IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~(a') * (~(b') * (~F)(~(c') + ~(d') * ~x)) ^ ~p) ^ ~(n'), ~x) => With([v = ActivateTrig((~F)(~c + ~d * ~x))], (((~a) ^ IntPart(~n) * (~a * (~b * v) ^ ~p) ^ FracPart(~n)) / (~b * v) ^ (~p * FracPart(~n))) * Antiderivative(ActivateTrig(~u) * (~b * v) ^ (~n * ~p), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (InertTrigQ(~F) && (Not(IntegerQ(~n)) && Not(IntegerQ(~p))))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], ShowStep("", "Antiderivative(F(tan(a+b*x)),x)", "1/b*Subst(Antiderivative(F(x)/(1+x^2),x),x,tan(a+b*x))", Hold( With([d = FreeFactors(tan(v), x)], Dist(d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), tan(v)/d, u, x), x), x, tan(v)/d), x)))) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(tan(v), x), u, x)) <-- SimplifyFlag && InverseFunctionFreeQ(u, x) && Not( MatchQ(u, (~v')*((~c')*tan((~w))^(~n')*tan((~z))^(~n'))^(~p') <-- FreeQ([c, p], x) && IntegerQ(n) && LinearQ(w, x) && EqQ(z, 2*w))), Antiderivative((~u), ~x) := With([v = FunctionOfTrig(u, x)], With([d = FreeFactors(tan(v), x)], Dist(d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), tan(v)/d, u, x), x), x, tan(v)/d), x)) <-- Not(FalseQ(v)) && FunctionOfQ(NonfreeFactors(tan(v), x), u, x)) <-- InverseFunctionFreeQ(u, x) && Not( MatchQ(u, (~v')*((~c')*tan((~w))^(~n')*tan((~z))^(~n'))^(~p') <-- FreeQ([c, p], x) && IntegerQ(n) && LinearQ(w, x) && EqQ(z, 2*w))))
    @apply_utils Antiderivative(~u * (~(c') * sin(~v)) ^ ~m, ~x) => With([w = FunctionOfTrig((~u * sin(~v / 2) ^ (2 * ~m)) / (~c * tan(~v / 2)) ^ ~m, ~x)], (((~c * sin(~v)) ^ ~m * (~c * tan(~v / 2)) ^ ~m) / sin(~v / 2) ^ (2 * ~m)) * Antiderivative((~u * sin(~v / 2) ^ (2 * ~m)) / (~c * tan(~v / 2)) ^ ~m, ~x) <-- Not(FalseQ(w)) && FunctionOfQ(NonfreeFactors(tan(w), ~x), (~u * sin(~v / 2) ^ (2 * ~m)) / (~c * tan(~v / 2)) ^ ~m, ~x)) <-- FreeQ(~c, ~x) && (LinearQ(~v, ~x) && (IntegerQ(~m + 1 / 2) && (Not(SumQ(~u)) && InverseFunctionFreeQ(~u, ~x))))
    @apply_utils Antiderivative(~(u') * (~(a') * tan(~(c') + ~(d') * ~x) ^ ~(n') + ~(b') * sec(~(c') + ~(d') * ~x) ^ ~(n')) ^ ~p, ~x) => Antiderivative(ActivateTrig(~u) * sec(~c + ~d * ~x) ^ (~n * ~p) * (~b + ~a * sin(~c + ~d * ~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IntegersQ(~n, ~p)
    @apply_utils Antiderivative(~(u') * (~(a') * cot(~(c') + ~(d') * ~x) ^ ~(n') + ~(b') * csc(~(c') + ~(d') * ~x) ^ ~(n')) ^ ~p, ~x) => Antiderivative(ActivateTrig(~u) * csc(~c + ~d * ~x) ^ (~n * ~p) * (~b + ~a * cos(~c + ~d * ~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IntegersQ(~n, ~p)
    @apply_utils Antiderivative(~u * (~a * (~F)(~(c') + ~(d') * ~x) ^ ~(p') + ~(b') * (~F)(~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(n'), ~x) => Antiderivative(ActivateTrig(~u * (~F)(~c + ~d * ~x) ^ (~n * ~p) * (~a + ~b * (~F)(~c + ~d * ~x) ^ (~q - ~p)) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p, ~q], ~x) && (InertTrigQ(~F) && (IntegerQ(~n) && PosQ(~q - ~p)))
    @apply_utils Antiderivative(~u * (~a * (~F)(~(d') + ~(e') * ~x) ^ ~(p') + ~(b') * (~F)(~(d') + ~(e') * ~x) ^ ~(q') + ~(c') * (~F)(~(d') + ~(e') * ~x) ^ ~(r')) ^ ~(n'), ~x) => Antiderivative(ActivateTrig(~u * (~F)(~d + ~e * ~x) ^ (~n * ~p) * (~a + ~b * (~F)(~d + ~e * ~x) ^ (~q - ~p) + ~c * (~F)(~d + ~e * ~x) ^ (~r - ~p)) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q, ~r], ~x) && (InertTrigQ(~F) && (IntegerQ(~n) && (PosQ(~q - ~p) && PosQ(~r - ~p))))
    @apply_utils Antiderivative(~u * (~a + ~(b') * (~F)(~(d') + ~(e') * ~x) ^ ~(p') + ~(c') * (~F)(~(d') + ~(e') * ~x) ^ ~(q')) ^ ~(n'), ~x) => Antiderivative(ActivateTrig(~u * (~F)(~d + ~e * ~x) ^ (~n * ~p) * (~b + ~a * (~F)(~d + ~e * ~x) ^ -(~p) + ~c * (~F)(~d + ~e * ~x) ^ (~q - ~p)) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (InertTrigQ(~F) && (IntegerQ(~n) && NegQ(~p)))
    @apply_utils Antiderivative(~(u') * (~(a') * cos(~(c') + ~(d') * ~x) + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => Antiderivative(ActivateTrig(~u) * (~a * ℯ ^ ((-(~a) / ~b) * (~c + ~d * ~x))) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative(~u, ~x) => Antiderivative(TrigSimplify(~u), ~x) <-- TrigSimplifyQ(~u)
    @apply_utils Antiderivative(~(u') * (~a * ~v) ^ ~p, ~x) => With([uu = ActivateTrig(~u), vv = ActivateTrig(~v)], (((~a) ^ IntPart(~p) * (~a * vv) ^ FracPart(~p)) / vv ^ FracPart(~p)) * Antiderivative(uu * vv ^ ~p, ~x)) <-- FreeQ([~a, ~p], ~x) && (Not(IntegerQ(~p)) && Not(InertTrigFreeQ(~v)))
    @apply_utils Antiderivative(~(u') * ((~v) ^ ~m) ^ ~p, ~x) => With([uu = ActivateTrig(~u), vv = ActivateTrig(~v)], ((vv ^ ~m) ^ FracPart(~p) / vv ^ (~m * FracPart(~p))) * Antiderivative(uu * vv ^ (~m * ~p), ~x)) <-- FreeQ([~m, ~p], ~x) && (Not(IntegerQ(~p)) && Not(InertTrigFreeQ(~v)))
    @apply_utils Antiderivative(~(u') * ((~v) ^ ~(m') * (~w) ^ ~(n')) ^ ~p, ~x) => With([uu = ActivateTrig(~u), vv = ActivateTrig(~v), ww = ActivateTrig(~w)], ((vv ^ ~m * ww ^ ~n) ^ FracPart(~p) / (vv ^ (~m * FracPart(~p)) * ww ^ (~n * FracPart(~p)))) * Antiderivative(uu * vv ^ (~m * ~p) * ww ^ (~n * ~p), ~x)) <-- FreeQ([~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (Not(InertTrigFreeQ(~v)) || Not(InertTrigFreeQ(~w))))
    @apply_utils Antiderivative(~u, ~x) => With([v = ExpandTrig(~u, ~x)], Antiderivative(v, ~x) <-- SumQ(v)) <-- Not(InertTrigFreeQ(~u))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([w = Block([$ShowSteps = False, $StepCounter = Null], Antiderivative(SubstFor( 1/(1 + FreeFactors(tan(FunctionOfTrig(u, x)/2), x)^2*x^2), tan(FunctionOfTrig(u, x)/2)/ FreeFactors(tan(FunctionOfTrig(u, x)/2), x), u, x), x))], ShowStep("", "Antiderivative(F(sin(a+b*x),cos(a+b*x)),x)", "2/b*Subst(Antiderivative(1/(1+x^2)*F(2*x/(1+x^2),(1-x^2)/(1+x^2)),x),x, tan((a+b*x)/2))", Hold( Module([v = FunctionOfTrig(u, x), d], d = FreeFactors(tan(v/2), x); Dist(2*d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), tan(v/2)/d, u, x), x), x, tan(v/2)/d), x)))) <-- CalculusFreeQ(w, x)) <-- SimplifyFlag && InverseFunctionFreeQ(u, x) && Not(FalseQ(FunctionOfTrig(u, x))), Antiderivative((~u), ~x) := With([w = Block([$ShowSteps = False, $StepCounter = Null], Antiderivative(SubstFor( 1/(1 + FreeFactors(tan(FunctionOfTrig(u, x)/2), x)^2*x^2), tan(FunctionOfTrig(u, x)/2)/ FreeFactors(tan(FunctionOfTrig(u, x)/2), x), u, x), x))], Module([v = FunctionOfTrig(u, x), d], d = FreeFactors(tan(v/2), x); Dist(2*d/Coefficient(v, x, 1), Subst(Antiderivative(SubstFor(1/(1 + d^2*x^2), tan(v/2)/d, u, x), x), x, tan(v/2)/d), x)) <-- CalculusFreeQ(w, x)) <-- InverseFunctionFreeQ(u, x) && Not(FalseQ(FunctionOfTrig(u, x))))
    #= If(TrueQ($LoadShowSteps), Antiderivative((~u),~x) := With([v=FunctionOfTrig(u,x)], ShowStep("","Antiderivative(F(sin(a+b*x),cos(a+b*x)),x)","2/b*Subst(Antiderivative(1/(1+x^ 2)*F(2*x/(1+x^2),(1-x^2)/(1+x^2)),x),x,tan((a+b*x)/2))",Hold( With([d=FreeFactors(tan(v/2),x)], Dist(2*d/Coefficient(v,x,1),Subst(Antiderivative(SubstFor(1/(1+d^2*x^2),tan(v/ 2)/d,u,x),x),x,tan(v/2)/d),x)))) <-- Not(FalseQ(v))) <-- SimplifyFlag &&  InverseFunctionFreeQ(u,x), Antiderivative((~u),~x) := With([v=FunctionOfTrig(u,x)], With([d=FreeFactors(tan(v/2),x)], Dist(2*d/Coefficient(v,x,1),Subst(Antiderivative(SubstFor(1/(1+d^2*x^2),tan(v/ 2)/d,u,x),x),x,tan(v/2)/d),x)) <-- Not(FalseQ(v))) <--  InverseFunctionFreeQ(u,x)) =#
    @apply_utils Antiderivative(~u, ~x) => With([v = ActivateTrig(~u)], CannotIntegrate(v, ~x)) <-- Not(InertTrigFreeQ(~u))
end

max_stars_repo_path: src/rules/4 Trig functions/4.7 Miscellaneous/4.7.5 Inert trig functions.jl


text:
using DynamicMovementPrimitives
using Base.Test
import DynamicMovementPrimitives: get_sched_sig, get_centers_linear, kernel_matrix, _1

# Setup ============================================
Nbasis = 20
αzC = 25.
αzPT = 10.
αx = 1.#αz/3
# bz = 1
# close("all")

optsC = DMPopts(Nbasis,αx,αzC, :canonical)
optsT = DMPopts(Nbasis,αx,αzPT, :time)
optsP = DMPopts(Nbasis,αx,αzPT, :position)
@test optsC.βz == αzC/4

y = [linspace(0,1,100); ones(10)]
T = length(y)
t = linspace(0,T,T)
τ = t[end]/3
h = t[2]-t[1]
y = [y 0.5y]
ẏ = centraldiff(y) / h
ÿ = centraldiff(ẏ) / h
g=y[end,:][:]

# Tests canonical ============================================

x = get_sched_sig(:canonical,αx,τ,t,y,g)
c,σ2    = get_centers_linear(Nbasis,x)
Ψ = kernel_matrix(x,c,σ2)
@test size(x) == (T,)
@test size(c) == (Nbasis,1)
@test size(σ2) == (Nbasis,1)
@test size(Ψ) == (T,Nbasis)

dmp = fit(y,ẏ,ÿ,t,optsC)
f = force(dmp,x)
@test size(f) == (T,2)
@test all(isfinite(f))
f = force(dmp,x[1],1)
@test isa(f,Number)
@test isfinite(f)

a = acceleration(dmp,y[1],ẏ[1],x[1],g[1],1)
@test isa(a,Number)
@test isfinite(a)
a = acceleration(dmp,_1(y),_1(ẏ),x[1],g)
@test size(a) == (2,)
a = acceleration(dmp,y,ẏ,x,g)
@test size(a) == (T,2)
@test all(isfinite(a))

# Tests time     ============================================
x = get_sched_sig(:time,αx,τ,t,y,g)
c,σ2    = get_centers_linear(Nbasis,x)
Ψ = kernel_matrix(x,c,σ2)
@test size(x) == (T,)
@test size(c) == (Nbasis,1)
@test size(σ2) == (Nbasis,1)
@test size(Ψ) == (T,Nbasis)

dmp = fit(y,ẏ,ÿ,t,optsT)
f = force(dmp,x)
@test size(f) == (T,2)
@test all(isfinite(f))
f = force(dmp,x[1],1)
@test isa(f,Number)
@test isfinite(f)

# Tests position ============================================
x = get_sched_sig(:position,αx,τ,t,y,g)
c,σ2    = get_centers_linear(Nbasis,x)
Ψ = kernel_matrix(x[:,1],c[:,1],σ2[:,1])
@test size(x) == (T,2)
@test size(c) == (Nbasis,2)
@test size(σ2) == (Nbasis,2)
@test size(Ψ) == (T,Nbasis)

dmp = fit(y,ẏ,ÿ,t,optsP)
f = force(dmp,x)
@test size(f) == (T,2)
@test all(isfinite(f))
f = force(dmp,x[1],1)
@test isa(f,Number)
@test isfinite(f)

a = acceleration(dmp,y[1],ẏ[1],x[1],g[1],1)
@test isa(a,Number)
@test isfinite(a)
a = acceleration(dmp,_1(y),_1(ẏ),x[1,:][:],g)
@test size(a) == (2,)
a = acceleration(dmp,y,ẏ,x,g)
@test size(a) == (T,2)
@test all(isfinite(a))

dmp = fit(y,ẏ,ÿ,t,optsC)
tout,youtC,ẏout,xout = solve(dmp)
@test tout == t
@test abs(youtC .- y) |> sum < 2
@test abs(ẏout .- ẏ) |> sum < 0.3
# plotdmp(dmp)

dmp = fit(y,ẏ,ÿ,t,optsT)
tout,youtT,ẏout,xout = solve(dmp)
@test abs(youtT .- y) |> sum < 2
@test abs(ẏout .- ẏ) |> sum < 0.3
# plotdmp(dmp)

dmp = fit(y,ẏ,ÿ,t,optsP)
tout,youtP,ẏout,xout = solve(dmp)
@test abs(youtP .- y) |> sum < 6
@test abs(ẏout .- ẏ) |> sum < 0.3
# plotdmp(dmp)
# plotdmpphase(dmp)

# Test twolink
include("../src/two_link.jl")
using TwoLink
p,pd,pdd = traj(0,1,0:100)
@test p[1] == 0
@test p[end] == 1
@test pd[1] == 0
@test pd[end] == 0
@test length(p) == length(pd) == length(pdd) == 101

cpoints = [0.5 -0.5;
        0.5 0;
        1.5 0;
        1.5 -0.5]

ctraj = connect_points(cpoints,20)[1]
jtraj = inverse_kin(ctraj,:up)
ctraj2 = forward_kin(jtraj)

@test ctraj ≈ ctraj2


text:
module stochastic_master

export master, master_dynamic

using ...bases, ...states, ...operators
using ...operators_dense, ...operators_sparse
using ...timeevolution
using LinearAlgebra
import ...timeevolution: integrate_stoch, recast!, QO_CHECKS
import ...timeevolution.timeevolution_master: dmaster_h, dmaster_nh, dmaster_h_dynamic, check_master

const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Nothing}
const DiffArray = Union{Vector{ComplexF64}, Array{ComplexF64, 2}}

"""
    stochastic.master(tspan, rho0, H, J, C; <keyword arguments>)

Time-evolution according to a stochastic master equation.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Deterministic part of the Hamiltonian.
* `J`: Vector containing all deterministic
        jump operators which can be of any arbitrary operator type.
* `C`: Vector containing the stochastic operators for a superoperator
        of the form `C[i]*rho + rho*Cdagger[i]`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
        operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::T, H::AbstractOperator{B,B},
                J::Vector, C::Vector;
                rates::DecayRates=nothing,
                Jdagger::Vector=dagger.(J), Cdagger::Vector=dagger.(C),
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B<:Basis,T<:DenseOperator{B,B}}

    tmp = copy(rho0)

    n = length(C)

    dmaster_stoch(dx::DiffArray, t::Float64, rho::T, drho::T, n::Int) =
            dmaster_stochastic(dx, rho, C, Cdagger, drho, n)

    isreducible = check_master(rho0, H, J, Jdagger, rates) && check_master_stoch(rho0, C, Cdagger)
    if !isreducible
        dmaster_h_determ(t::Float64, rho::T, drho::T) =
            dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
        integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch, rho0, fout, n; kwargs...)
    else
        Hnh = copy(H)
        if typeof(rates) == Matrix{Float64}
            for i=1:length(J), j=1:length(J)
                Hnh -= 0.5im*rates[i,j]*Jdagger[i]*J[j]
            end
        elseif typeof(rates) == Vector{Float64}
            for i=1:length(J)
                Hnh -= 0.5im*rates[i]*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= 0.5im*Jdagger[i]*J[i]
            end
        end
        Hnhdagger = dagger(Hnh)

        dmaster_nh_determ(t::Float64, rho::T, drho::T) =
            dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
        integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch, rho0, fout, n; kwargs...)
    end
end
master(tspan, psi0::Ket, args...; kwargs...) = master(tspan, dm(psi0), args...; kwargs...)

"""
    stochastic.master_dynamic(tspan, rho0, fdeterm, fstoch; <keyword arguments>)

Time-evolution according to a stochastic master equation with a
dynamic Hamiltonian and J.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `fdeterm`: Function `f(t, rho) -> (H, J, Jdagger)` or
        `f(t, rho) -> (H, J, Jdagger, rates)` giving the deterministic
        part of the master equation.
* `fstoch`: Function `f(t, rho) -> (C, Cdagger)` giving the stochastic
        superoperator of the form `C[i]*rho + rho*Cdagger[i]`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch` and all optional functions. If unset, the number
        is calculated automatically from the function outputs. NOTE: Set this
        number if you want to avoid an initial calculation of function outputs!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan::Vector{Float64}, rho0::T, fdeterm::Function, fstoch::Function;
                rates::DecayRates=nothing,
                fout::Union{Function,Nothing}=nothing,
                noise_processes::Int=0,
                kwargs...) where {B<:Basis,T<:DenseOperator{B,B}}

    tmp = copy(rho0)

    if noise_processes == 0
        fs_out = fstoch(0, rho0)
        n = length(fs_out[1])
    else
        n = noise_processes
    end

    dmaster_determ(t::Float64, rho::T, drho::T) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp)
    dmaster_stoch(dx::DiffArray, t::Float64, rho::T, drho::T, n::Int) =
        dmaster_stoch_dynamic(dx, t, rho, fstoch, drho, n)
    integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; kwargs...)
end
master_dynamic(tspan::Vector{Float64}, psi0::Ket, args...; kwargs...) = master_dynamic(tspan, dm(psi0), args...; kwargs...)

# Derivative functions
function dmaster_stochastic(dx::Vector{ComplexF64}, rho::T,
            C::Vector, Cdagger::Vector, drho::T, ::Int) where {B<:Basis,T<:DenseOperator{B,B}}
    recast!(dx, drho)
    operators.gemm!(1, C[1], rho, 0, drho)
    operators.gemm!(1, rho, Cdagger[1], 1, drho)
    drho.data .-= tr(drho)*rho.data
end
function dmaster_stochastic(dx::Array{ComplexF64, 2}, rho::T,
            C::Vector, Cdagger::Vector, drho::T, n::Int) where {B<:Basis,T<:DenseOperator{B,B}}
    for i=1:n
        dx_i = @view dx[:, i]
        recast!(dx_i, drho)
        operators.gemm!(1, C[i], rho, 0, drho)
        operators.gemm!(1, rho, Cdagger[i], 1, drho)
        drho.data .-= tr(drho)*rho.data
        recast!(drho, dx_i)
    end
end

function dmaster_stoch_dynamic(dx::DiffArray, t::Float64, rho::T,
            f::Function, drho::T, n::Int) where {B<:Basis,T<:DenseOperator{B,B}}
    result = f(t, rho)
    QO_CHECKS[] && @assert 2 == length(result)
    C, Cdagger = result
    QO_CHECKS[] && check_master_stoch(rho, C, Cdagger)
    dmaster_stochastic(dx, rho, C, Cdagger, drho, n)
end

function integrate_master_stoch(tspan, df::Function, dg::Function,
                        rho0::DenseOperator, fout::Union{Nothing, Function},
                        n::Int;
                        kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    x0 = reshape(rho0.data, length(rho0))
    state = copy(rho0)
    dstate = copy(rho0)
    integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
end

function check_master_stoch(rho0::DenseOperator{B,B}, C::Vector, Cdagger::Vector) where B<:Basis
    # TODO: replace type checks by dispatch; make types of C known
    @assert length(C) == length(Cdagger)
    isreducible = true
    for c=C
        @assert isa(c, AbstractOperator{B,B})
        if !(isa(c, DenseOperator) || isa(c, SparseOperator))
            isreducible = false
        end
    end
    for c=Cdagger
        @assert isa(c, AbstractOperator{B,B})
        if !(isa(c, DenseOperator) || isa(c, SparseOperator))
            isreducible = false
        end
    end
    isreducible
end

# TODO: Speed up by recasting to n-d arrays, remove vector methods
function recast!(x::Union{Vector{ComplexF64}, SubArray{ComplexF64, 1}}, rho::DenseOperator{B,B,T}) where {B<:Basis,T<:Matrix{ComplexF64}}
    rho.data = reshape(x, size(rho.data))
end
recast!(state::DenseOperator{B,B}, x::SubArray{ComplexF64, 1}) where B<:Basis = (x[:] = state.data)
recast!(state::DenseOperator{B,B}, x::Vector{ComplexF64}) where B<:Basis = nothing

end # module


text:
#!/usr/bin/env julia

using Eigenbroetler

function scaled_sinusoid(x::Number, y::Number)
    r = sqrt(x * x + y * y)
    return sin(100 * r) / (r + 1)
end

eb = fft(Eigenbrot(scaled_sinusoid, 512, 512))
save("scaled_sinusoid.png", eb, ImageSetting(Magn, LogScale),
     colours = joinpath(dirname(@__FILE__), "iron.map"))


text:
@testset "DR types" begin
    grid_smol = Dolo.SmolyakGrid{2}([0.0, 0.0], [1.0, 1.0,], [3,3])
    grid_exo = Dolo.EmptyGrid()
    grid_exo2 = Dolo.UnstructuredGrid{2}(rand(1, 2))
    grid_rand = Dolo.RandomGrid{2}([0.0, 0.0], [1.0, 1.0,], 50)
    pts = rand(200, 2)
    pts_vec = copy(reshape(reinterpret(Dolo.Point{2}, vec(copy(pts'))), (size(pts, 1),)))
    test_sets = [
        (grid_smol, Dolo.SmolyakDR, Dolo.Smolyak),
        (grid_rand, Dolo.CompletePolyDR, Dolo.CompletePolnomial{3}),
        (grid_rand, Dolo.CompletePolyDR, Dolo.CompletePolnomial),
        (grid_rand, Dolo.CompletePolyDR, Dolo.CompletePolnomial{2}),
        (grid_smol, Dolo.CompletePolyDR, Dolo.CompletePolnomial{3}),
        (grid_smol, Dolo.CompletePolyDR, Dolo.CompletePolnomial),
        (grid_smol, Dolo.CompletePolyDR, Dolo.CompletePolnomial{2}),
        (grid_exo2, Dolo.CompletePolyDR, Dolo.CompletePolnomial{3}),
        (grid_exo2, Dolo.CompletePolyDR, Dolo.CompletePolnomial),
        (grid_exo2, Dolo.CompletePolyDR, Dolo.CompletePolnomial{2}),
    ]

    @testset for (grid_endo, dr_name, drT) in test_sets

       g = Dolo.nodes(Matrix, grid_endo)
       sg_vals = [sin.(g[:, 1]) cos.(g[:, 2])]
       sg_vals_vec = copy(reshape(reinterpret(Dolo.Point{2}, vec(copy(sg_vals'))), (size(g,1),)))

        dr1 = dr_name(grid_exo, grid_endo, Val{2})
        dr2 = Dolo.DecisionRule(grid_exo, grid_endo, Val{2}, drT)
        for dr in [dr1, dr2]
            Dolo.set_values!(dr, [sg_vals])
            out = Dolo.evaluate(dr, pts)
            @test maximum(abs, out - Dolo.evaluate(dr, pts_vec)) < 1e-14  ### ?
            Dolo.set_values!(dr, [sg_vals_vec])
            @test maximum(abs, out - Dolo.evaluate(dr, pts_vec)) < 1e-14
            @test maximum(abs, out - Dolo.evaluate(dr, pts)) < 1e-14
            @test_throws ErrorException Dolo.set_values!(dr, [sg_vals_vec, sg_vals_vec])
        end

        dr1 = dr_name(grid_exo2, grid_endo, Val{2})
        dr2 = dr_name(grid_exo2, grid_endo, Val{2})
        for dr in [dr1, dr2]
            Dolo.set_values!(dr, [sg_vals])
            out = Dolo.evaluate(dr, 1, pts)
            @test maximum(abs, out - Dolo.evaluate(dr, 1, pts_vec)) < 1e-14
            Dolo.set_values!(dr, [sg_vals_vec])
            @test maximum(abs, out - Dolo.evaluate(dr, 1, pts_vec)) < 1e-14
            @test maximum(abs, out - Dolo.evaluate(dr, 1, pts)) < 1e-14
            @test_throws MethodError Dolo.evaluate(dr, pts)
            @test_throws MethodError Dolo.evaluate(dr, pts_vec)
            @test_throws MethodError Dolo.evaluate(dr, 1, pts_vec[1])
            @test_throws BoundsError Dolo.evaluate(dr, 2, pts_vec)
            @test_throws ErrorException Dolo.set_values!(dr, [sg_vals_vec, sg_vals_vec])
        end
    end
end


text:
using JuMP
using CPLEX

include("parser.jl")

start = time()

function Mymodel(MyFileName::String)
  P, Nm, Nf, C, G, A, T, init, y = read_instance(MyFileName)
  nbAllele = A * G
  theta = Array{Float64,1}(zeros(T))
  for i in 1:T
    theta[i] = init^((T-i)/(T-1))
  end

  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m, x[1:P] >= 0, Int)
  @variable(m, t[1:nbAllele] >= 0)
  @variable(m, z[1:nbAllele] >= 0)

  ## Constraints
  @constraint(m, sum(x[i] for i in 1:Nm) == P)
  @constraint(m, sum(x[i] for i in Nm+1:P) == P)
  @constraint(m, [j in 1:nbAllele], z[j] >= t[j] - sum(x[i] for i in 1:P if y[i,j] == 2))
  @constraint(m, [j in 1:nbAllele, h in 1:T], log(theta[h]) + ((t[j] - theta[h]) / theta[h]) >= sum(log((1-0.5*y[i,j]))*x[i] for i in 1:P if y[i,j] == 1))
  @constraint(m, [i in 1:P], x[i] <= 3) # limiter nombre d'enfant par individu

  ## Objective
  @objective(m, Min, sum(z[j] for j in 1:nbAllele))

  #resolution
  optimize!(m)

  println(solution_summary(m, verbose=true))
  x_vals = JuMP.getvalue.( m[:x] )
  println(x_vals)

  # penser a calculer proba

end

Mymodel("DivGenetique_ampl.dat")

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
"""
Reset the filter to initial state and covariance/distribution
"""
function reset!(pf::AbstractParticleFilter)
    s = state(pf)
    rn = rng(pf)
    for i = eachindex(s.xprev)
        s.xprev[i] = rand(rn, initial_density(pf))
        s.x[i] = copy(s.xprev[i])
    end
    fill!(s.w, -log(num_particles(pf)))
    fill!(s.we, 1/num_particles(pf))
    s.t[] = 1
end

function predict!(kf::AbstractKalmanFilter, u, t::Integer = index(kf))
    @unpack A,B,x,R,R1 = kf
    if ndims(A) == 3
        At = A[:,:,t]
        Bt = B[:,:,t]
    else
        At = A
        Bt = B
    end
    x .= At*x .+ Bt*u |> vec
    R .= symmetrize(At*R*At') + R1
    kf.t[] += 1
end

@inline function symmetrize(x::SArray)
    x = 0.5 .* (x .+ x')
    Symmetric(x)
end
@inline function symmetrize(x)
    x .+= x'
    x .*= 0.5
    Symmetric(x)
end

function correct!(kf::AbstractKalmanFilter, u, y, t::Integer = index(kf))
    @unpack C,D,x,R,R2 = kf
    if ndims(C) == 3
        Ct = C[:,:,t]
    else
        Ct = C
    end
    # Handle D separately in case it is 0
    if ndims(D) == 3
        Dt = D[:,:,t]
    else
        Dt = D
    end
    e   = y .- Ct*x
    if ! iszero(D)
        e .-= Dt*u
    end
    S   = symmetrize(Ct*R*Ct') + R2
    Sᵪ  = cholesky(S)
    K   = (R*Ct')/Sᵪ
    x .+= K*e
    R  .= symmetrize((I - K*Ct)*R) # WARNING against I .- A
    ll = logpdf(MvNormal(PDMat(S, Sᵪ)), e)# - 1/2*logdet(S) # logdet is included in logpdf
    ll, e
end

"""
    predict!(f, u, t = index(f))
Move filter state forward in time using dynamics equation and input vector `u`.
"""
function predict!(pf, u, t = index(pf))
    s = pf.state
    N = num_particles(s)
    if shouldresample(pf)
        j = resample(pf)
        propagate_particles!(pf, u, j, t)
        reset_weights!(s)
    else # Resample not needed
        s.j .= 1:N
        propagate_particles!(pf, u, t)
    end
    copyto!(s.xprev, s.x)
    pf.state.t[] += 1
end

"""
     ll = correct!(f, u, y, t = index(f))
Update state/covariance/weights based on measurement `y`,  returns loglikelihood.
"""
function correct!(pf, u, y, t = index(pf))
    measurement_equation!(pf, u, y, t)
    ll = logsumexp!(state(pf))
    ll, 0
end

"""
    ll, e = update!(f::AbstractFilter, u, y, t = index(f))
Perform one step of `predict!` and `correct!`, returns loglikelihood and prediction error
"""
function update!(f::AbstractFilter, u, y, t = index(f))
    ll_e = correct!(f, u, y, t)
    predict!(f, u, t)
    ll_e
end

function update!(pf::AuxiliaryParticleFilter, u, y, y1, t = index(pf))
    ll_e = correct!(pf, u, y, t)
    predict!(pf, u, y1, t)
    ll_e
end

function predict!(pf::AuxiliaryParticleFilter, u, y1, t = index(pf))
    s = state(pf)
    propagate_particles!(pf.pf, u, t, nothing)# Propagate without noise
    λ  = s.we
    λ .= 0
    measurement_equation!(pf.pf, u, y1, t, λ)
    s.w .+= λ
    expnormalize!(s.w) # w used as buffer
    j = resample(ResampleSystematic, s.w , s.j, s.bins)
    reset_weights!(s)
    permute_with_buffer!(s.x, s.xprev, j)
    add_noise!(pf.pf)

    s.t[] += 1
    copyto!(s.xprev, s.x)
end

function predict!(pf::AuxiliaryParticleFilter{<:AdvancedParticleFilter},u, y, t = index(pf))
    s = state(pf)
    propagate_particles!(pf.pf, u, t, nothing)# Propagate without noise
    λ  = s.we
    λ .= 0
    measurement_equation!(pf.pf, u, y, t, λ)
    s.w .+= λ
    expnormalize!(s.w) # w used as buffer
    j = resample(ResampleSystematic, s.w , s.j, s.bins)
    reset_weights!(s)
    propagate_particles!(pf.pf, u, j, t)# Propagate with noise and permutation

    s.t[] += 1
    copyto!(s.xprev, s.x)
end

(kf::KalmanFilter)(u, y, t = index(kf)) =  update!(kf, u, y, t)
(kf::UnscentedKalmanFilter)(u, y, t = index(kf)) =  update!(kf, u, y, t)
(pf::ParticleFilter)(u, y, t = index(pf)) =  update!(pf, u, y, t)
(pf::AuxiliaryParticleFilter)(u, y, y1, t = index(pf)) =  update!(pf, u, y, y1, t)
(pf::AdvancedParticleFilter)(u, y, t = index(pf)) =  update!(pf, u, y, t)
(pf::SigmaFilter)(u, y, t = index(pf)) =  update!(pf, u, y, t)

"""
    x,xt,R,Rt,ll = forward_trajectory(kf::AbstractKalmanFilter, u::Vector, y::Vector)

Run a Kalman filter forward

# Returns:
- `x`: predictions
- `xt`: filtered estimates
- `R`: predicted covariance matrices
- `Rt`: filter covariances
- `ll`: loglik
"""
function forward_trajectory(kf::AbstractKalmanFilter, u::AbstractVector, y::AbstractVector)
    reset!(kf)
    T    = length(y)
    x    = Array{particletype(kf)}(undef,T)
    xt   = Array{particletype(kf)}(undef,T)
    R    = Array{covtype(kf)}(undef,T)
    Rt   = Array{covtype(kf)}(undef,T)
    ll   = zero(eltype(particletype(kf)))
    for t = 1:T
        x[t]  = state(kf)      |> copy
        R[t]  = covariance(kf) |> copy
        ll   += correct!(kf, u[t], y[t], t)[1]
        xt[t] = state(kf)      |> copy
        Rt[t] = covariance(kf) |> copy
        predict!(kf, u[t], t)
    end
    x,xt,R,Rt,ll
end

"""
    x,w,we,ll = forward_trajectory(pf, u::AbstractVector, y::AbstractVector)
Run the particle filter for a sequence of inputs and measurements. Return particles, weights, expweights and loglikelihood

If [MonteCarloMeasurements.jl](https://github.com/baggepinnen/MonteCarloMeasurements.jl) is loaded, you may transform the output particles to `Matrix{MonteCarloMeasurements.Particles}` using `Particles(x,we)`. Internally, the particles are then resampled such that they all have unit weight. This is conventient for making use of the [plotting facilities of MonteCarloMeasurements.jl](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/#Plotting-1).
"""
function forward_trajectory(pf, u::AbstractVector, y::AbstractVector)
    reset!(pf)
    T = length(y)
    N = num_particles(pf)
    x = Array{particletype(pf)}(undef,N,T)
    w = Array{Float64}(undef,N,T)
    we = Array{Float64}(undef,N,T)
    ll = 0.
    @inbounds for t = 1:T
        ll += correct!(pf, u[t], y[t], t) |> first
        x[:,t] .= particles(pf)
        w[:,t] .= weights(pf)
        we[:,t] .= expweights(pf)
        predict!(pf, u[t], t)
    end
    x,w,we,ll
end

function forward_trajectory(pf::AuxiliaryParticleFilter, u::AbstractVector, y::AbstractVector)
    reset!(pf)
    T = length(y)
    N = num_particles(pf)
    x = Array{particletype(pf)}(undef,N,T)
    w = Array{Float64}(undef,N,T)
    we = Array{Float64}(undef,N,T)
    ll = 0.
    @inbounds for t = 1:T
        ll += correct!(pf, u[t], y[t], t) |> first
        x[:,t] .= particles(pf)
        w[:,t] .= weights(pf)
        we[:,t] .= expweights(pf)
        t < T && predict!(pf, u[t], y[t+1], t)
    end
    x,w,we,ll
end

"""
    x,ll = mean_trajectory(pf, u::Vector{Vector}, y::Vector{Vector})

This Function resets the particle filter to the initial state distribution upon start
"""
mean_trajectory(pf, u::Vector, y::Vector) = reduce_trajectory(pf, u::Vector, y::Vector, weigthed_mean)
mode_trajectory(pf, u::Vector, y::Vector) = reduce_trajectory(pf, u::Vector, y::Vector, mode)

function reduce_trajectory(pf, u::Vector, y::Vector, f::F) where F
    reset!(pf)
    T = length(y)
    N = num_particles(pf)
    x = Array{particletype(pf)}(undef,T)
    ll = correct!(pf,u[1],y[1],1) |> first
    x[1] = f(state(pf))
    for t = 2:T
        ll += pf(u[t-1], y[t], t) |> first
        x[t] = f(pf)
    end
    x,ll
end

StatsBase.mode(pf::AbstractParticleFilter) = particles(pf)[findmax(expparticles(pf))[2]]

mode_trajectory(x::AbstractMatrix, we::AbstractMatrix) =  reduce(hcat,vec(x[findmax(we, dims=1)[2]]))'

function mean_trajectory(x::AbstractMatrix, we::AbstractMatrix)
    copy(reduce(hcat,vec(sum(x.*we,dims=1)))')
end

"""
    x,u,y = simulate(f::AbstractFilter,T::Int,du::Distribution, [N])
Simulate dynamical system forward in time, returns state sequence, inputs and measurements
`du` is a distribution of random inputs.

If [MonteCarloMeasurements.jl](https://github.com/baggepinnen/MonteCarloMeasurements.jl) is loaded, the argument `N::Int` can be supplied, in which case `N` simulations are done and the result is returned in the form of `Vector{MonteCarloMeasurements.Particles}`.
"""
function simulate(f::AbstractFilter,T::Int,du::Distribution)
    u = [rand(du) for t=1:T]
    simulate(f, u)
end

function simulate(f::AbstractFilter,u)
    y = similar(u)
    x = similar(u)
    x[1] = sample_state(f)
    T = length(u)
    for t = 1:T-1
        y[t] = sample_measurement(f,x[t], u[t], t)
        x[t+1] = sample_state(f, x[t], u[t], t)
    end
    y[T] = sample_measurement(f,x[T], u[T], T)
    x,u,y
end

"""
    x̂ = weigthed_mean(x,we)
Calculated weighted mean of particle trajectories. `we` are expweights.
"""
function weigthed_mean(x,we::AbstractVector)
    @assert sum(we) ≈ 1
    xh = zeros(size(x[1]))
    @inbounds @simd  for i = eachindex(x)
        xh .+= x[i].*we[i]
    end
    return xh
end
function weigthed_mean(x,we::AbstractMatrix)
    N,T = size(x)
    @assert sum(we) ≈ T
    xh = zeros(eltype(x), T)
    for t = 1:T
        @inbounds @simd for i = 1:N
            xh[t] += x[i,t].*we[i,t]
        end
    end
    return xh
end

"""
    x̂ = weigthed_mean(pf)
    x̂ = weigthed_mean(s::PFstate)
"""
weigthed_mean(s) = weigthed_mean(s.x,s.we)
weigthed_mean(pf::AbstractParticleFilter) = weigthed_mean(state(pf))
weigthed_mean(s::SigmaFilter) = weigthed_mean(s.x,s.we)
"""
Similar to [`weigthed_mean`](@ref), but returns covariances
"""
function weigthed_cov(x,we)
    N,T = size(x)
    n = length(x[1])
    [cov(copy(reshape(reinterpret(Float64, x[:,t]),n,N)),ProbabilityWeights(we[:,t]), dims=2, corrected=true) for t = 1:T]
end