text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module BooleanOperations
using NumberTheory

export Not, And, Dif, Cnimp, Xor, Or, Nor, Eqv, Cimp, Imp, Nand
export Bits, BoolOps, BinDigits
export V035327, V003817, V129760, V142151, V080079, V086799, V163617
export V038712, V006257, V048724, V003188, V038554, V048735, V213370
export V080940, V135521, V051933, V135481, V327987
export is327988, L327988, L280172

"""

The 16 boolean operations

* FALSE, AND, DIF, PRJ1, CNIMP, PRJ2, XOR, OR, NOR, EQV, NEG2, CIMP, NEG1, IMP, NAND, TRUE

evaluated bitwise on the binary expansions of integers.

The main function is

* Bits(op::String, n::Int, k::Int, algo=max)

where op are the above acronyms. If the lengths of the binary expansions of ``n`` and ``k`` are different, the parameter 'algo=min' offers the choice to reduce the range of comparison to the smaller range or to the larger range, 'algo=max', which will first pad the representation of the smaller operand with 0's up to the length of the larger operand.

* Not, And, Dif, Cnimp, Xor, Or, Nor, Eqv, Cimp, Imp, Nand, BinDigits, BoolOps, Bits, V035327, V003817, V129760, V142151, V080079, V086799, V038712, V163617, V006257, V048724, V003188, V038554, V048735, V213370, V080940, V135521,  , V280172.
"""
const ModuleBooleanOperations = ""

"""
Return the binary expansions of nonnegative ``n`` and ``k``. If algo=max then pad the resulting int arrays with 0 to make them equally long. For example

    BinDigits``(12, 17, max) = ([0, 0, 1, 1, 0], [1, 0, 0, 0, 1])``.

If algo=min then the arrays are truncated to the length of the smaller operand.

    BinDigits``(12, 17, min) = ([0, 0, 1, 1], [1, 0, 0, 0])``.

These are the lists which are itemwise compared by the logical operators.
"""
function BinDigits(n::Int, k::Int, algo=max)
    N = n == 0 ? [0] : digits(n, base = 2)
    K = k == 0 ? [0] : digits(k, base = 2)
    lenN = length(N)
    lenK = length(K)

    if algo == max
        maxlen = max(lenN, lenK)
        if maxlen > lenN
            N = vcat(N, [0 for i ∈ 1:maxlen-lenN])
        else
            K = vcat(K, [0 for i ∈ 1:maxlen-lenK])
        end
    else
        minlen = min(lenN, lenK)
        if lenN > minlen
            N = N[1:minlen]
        else
            K = K[1:minlen]
        end
    end
    return (N, K)
end

"""

Return the bitwise negation of ``N``.
"""
Not(N) = [1 - n for n ∈ N]

"""

Return the bitwise conjunction of ``N`` and ``K``.
"""
And(N, K) = [N[i] & k for (i, k) ∈ enumerate(K)]

"""

Return the bitwise difference of ``N`` and ``K``.
"""
Dif(N, K) = [N[i] & (1 - k) for (i, k) ∈ enumerate(K)]

"""

Return the bitwise converse nonimplication of ``N`` and ``K``.
"""
Cnimp(N, K) = [(1 - N[i]) & k for (i, k) ∈ enumerate(K)]

"""

Return the bitwise exclusive disjunction of ``N`` and ``K``.
"""
Xor(N, K) = [xor(N[i], k) for (i, k) ∈ enumerate(K)]

"""

Return the bitwise disjunction of ``N`` and ``K``.
"""
Or(N, K) = [N[i] | k for (i, k) ∈ enumerate(K)]

"""

Return the bitwise joint denial of ``N`` and ``K``.
"""
Nor(N, K) = [1 - (N[i] | k) for (i, k) ∈ enumerate(K)]

"""

Return the bitwise equivalence of ``N`` and ``K``.
"""
Eqv(N, K) = [xor(1 - N[i], k) for (i, k) ∈ enumerate(K)]

"""

Return the bitwise converse implication of ``N`` and ``K``.
"""
Cimp(N, K) = [N[i] | (1 - k) for (i, k) ∈ enumerate(K)]

"""

Return the bitwise implication of ``N`` and ``K``.
"""
Imp(N, K) = [(1 - N[i]) | k for (i, k) ∈ enumerate(K)]

"""

Return the bitwise nonconjunction of ``N`` and ``K``.
"""
Nand(N, K) = [1 - (N[i] & k) for (i, k) ∈ enumerate(K)]

"""

Return the bitwise boolean operation represented by ``op`` applied to the binary expansions of the integers ``n`` and ``k``. ``op`` is an integer ``0 <= op < 16`` encoding the result of the operation in terms of the truth table.
"""
function Bits(op::Int, n::Int, k::Int, algo=max)
    if op < 0 || op > 15 || n < 0 || k < 0
        @error("op must be 0 <= op <= 15 and n,k >= 0. op=$op, n=$n, k=$k.")
        throw(DomainError(-1))
    end

    op == 0  && return 0
    op == 15 && return 1

    N, K = BinDigits(n, k, algo)

    if     op == 1   R = And(N, K)
    elseif op == 2   R = Dif(N, K)
    elseif op == 3   R = N
    elseif op == 4   R = Cnimp(N, K)
    elseif op == 5   R = K
    elseif op == 6   R = Xor(N, K)
    elseif op == 7   R = Or(N, K)
    elseif op == 8   R = Nor(N, K)
    elseif op == 9   R = Eqv(N, K)
    elseif op == 10  R = Not(K)
    elseif op == 11  R = Cimp(N, K)
    elseif op == 12  R = Not(N)
    elseif op == 13  R = Imp(N, K)
    elseif op == 14  R = Nand(N, K)
    end

    sum(r << (i-1)  for (i, r) ∈ enumerate(R))
end

"""

The list of acronyms for the 16 binary boolean operations.
"""
const BoolOps = ["FALSE", "AND", "DIF", "PRJ1", "CNIMP", "PRJ2", "XOR", "OR", "NOR", "EQV", "NEG2", "CIMP", "NEG1", "IMP", "NAND", "TRUE"]

"""

Return the bitwise boolean operation represented by ``op`` applied to the binary expansions of the integers ``n`` and ``k``. ``op`` is an acronym as given in the list 'BoolOps'.
"""
function Bits(op::String, n::Int, k::Int, algo=max)
    opnum = findfirst(isequal(op), BoolOps)
    Bits(opnum - 1, n, k, algo)
end

"""

Return ``n`` NAND ``n``.
"""
V035327(n) = Bits("NAND", n, n)
# also : V035327(n) = Bits("NEG1", n, n+1, min)

"""

Return ``n`` EQV ``n``.
"""
V003817(n) = Bits("EQV", n, n)
# sets a(0) = 1

"""

Return ``n`` AND ``n+1``, using max length.
"""
V129760(n) = Bits("AND", n, n+1)
# with offset 0

"""

Return ``n`` CIMP ``n+1``, using max length.
"""
V142151(n) = Bits("CIMP", n, n+1)

#"""
# Return n CNIMP n+1, using max length.
# Offset is 0, shift n left if you want the OEIS offset 1.
# Is exported from Module LandauConstants
#"""
#V006519(n) = Bits("CNIMP", n, n+1)

"""

Return ``n`` NEG1 ``n+1``, using max length.
"""
V080079(n) = Bits("NEG1", n, n+1)
# with offset 0

"""

Return ``n`` OR ``n+1``, using max length.
"""
V086799(n) = Bits("OR", n, n+1)
# with offset 0

"""

Return ``n`` OR ``2n``, using max length.
"""
V163617(n) = Bits("OR", n, n<<1)

"""

Return ``n`` XOR ``n+1``, using max length.
"""
V038712(n) = Bits("XOR", n, n+1)
# with offset 0

"""

Return max``(1, 2n) - (n`` EQV ``n)``, using max length.
"""
V006257(n) = max(1, 2n) - Bits("EQV", n, n)

"""

Return ``n`` XOR ``2n``, using max length.
"""
V048724(n) = Bits("XOR", n, n<<1)

"""

Return ``n`` XOR ``n>>1``, using max length.
"""
V003188(n) = Bits("XOR", n, n>>1)

"""

Return ``n`` XOR ``n>>1``, using min length.
"""
V038554(n) = Bits("XOR", n, n>>1, min)
# with a(1) = 1

"""

Return ``n`` AND ``n>>1``, using min length.
"""
V048735(n) = Bits("AND", n, n>>1, min)

"""

Return ``n`` AND ``n<<1``, using min length.
"""
V213370(n) = Bits("AND", n, n<<1, min)

"""

Return ``n`` CNIMP ``n+1``, using min length.
"""
V080940(n) = Bits("CNIMP", n, n+1, min)
# with a(0) = 1

"""

Return ``n`` XOR ``n+1``, using min length.
"""
V135521(n) = Bits("XOR", n, n+1, min)
# prepend a(0) = 1

"""

Return ``n`` XOR ``k``, using max length.
"""
V051933(n, k) = Bits("XOR", n, k)

"""

Return ``(k-1`` XOR ``n-k) + 1``, using max length.
"""
L280172(n) = [Bits("XOR", k-1, n-k) + 1 for k in 1:n]

"""

Return ``n+1`` CNIMP ``n``, using max length.
"""
V135481(n) = Bits("CNIMP", n+1, n)

"""

Return ∑``_{d|n} d & (n/d)``, where & is the bitwise AND operator.
"""
V327987(n) = sum([Bits("AND", d, div(n, d)) for d ∈ divisors(n)])

"""

Is V327987(n) = ∑``_{d|n} d & (n/d) = 0`` ?
"""
is327988(n) = V327987(n) == 0

"""

Return a list of the zeros of V327987 below (inc.)  ``max``.
"""
L327988(max) = [n for n ∈ 0:max if is327988(n)]

#START-TEST-########################################################

using Test

function test()
    @testset "BoolOps" begin
        A = [1, 1, 0, 0]
        B = [1, 0, 1, 0]
        S(R) = sum(r << (i-1)  for (i, r) ∈ enumerate(R))

        for n ∈ 0:15
            @test n == S([Bits(n, A[k], B[k]) for k ∈ 1:4])
        end
        for n ∈ 0:15
            @test n == S([Bits(BoolOps[n+1], A[k], B[k]) for k ∈ 1:4])
        end

        @test BinDigits(12, 17) == ([0, 0, 1, 1, 0], [1, 0, 0, 0, 1])
        @test [V163617(n) for n ∈ 0:8] == [0, 3, 6, 7, 12, 15, 14, 15, 24]
        @test [V035327(n) for n ∈ 0:8] == [1, 0, 1, 0, 3, 2, 1, 0, 7]
        @test [V129760(n) for n ∈ 0:8] == [0, 0, 2, 0, 4, 4, 6, 0, 8]
        @test [V142151(n) for n ∈ 0:8] == [0, 1, 2, 3, 6, 5, 6, 7, 14]
        @test [V003817(n) for n ∈ 0:8] == [1, 1, 3, 3, 7, 7, 7, 7, 15]
        @test [V080079(n) for n ∈ 0:8] == [1, 2, 1, 4, 3, 2, 1, 8, 7]
        @test [V086799(n) for n ∈ 0:8] == [1, 3, 3, 7, 5, 7, 7, 15, 9]
        @test [V038712(n) for n ∈ 0:8] == [1, 3, 1, 7, 1, 3, 1, 15, 1]
        @test [V006257(n) for n ∈ 0:8] == [0, 1, 1, 3, 1, 3, 5, 7, 1]

        @test [V080940(n) for n ∈ 0:8] == [1, 0, 1, 0, 1, 2, 1, 0, 1]
        @test [V135521(n) for n ∈ 0:8] == [1, 1, 1, 3, 1, 3, 1, 7, 1]

        @test [V048724(n) for n ∈ 0:8] == [0, 3, 6, 5, 12, 15, 10, 9, 24]
        @test [V003188(n) for n ∈ 0:8] == [0, 1, 3, 2, 6, 7, 5, 4, 12]
        @test [V038554(n) for n ∈ 0:8] == [0, 1, 1, 0, 2, 3, 1, 0, 4]
        @test [V048735(n) for n ∈ 0:8] == [0, 0, 0, 1, 0, 0, 2, 3, 0]
        @test [V213370(n) for n ∈ 0:8] == [0, 0, 0, 2, 0, 0, 4, 6, 0]

        @test [V135481(n) for n ∈ 0:8] == [0, 1, 0, 3, 0, 1, 0, 7, 0]
        @test [V327987(n) for n ∈ 0:8] == [0, 1, 0, 2, 2, 2, 4, 2, 0]
        @test L327988(99) == [0, 2, 8, 10, 26, 32, 34, 40, 50, 58, 74, 82]

        @test is327988(121) == false
        @test is327988(122) == true
    end
end

function demo()
    println()
    for n ∈ 0:15
        seq = [Bits(n, k, k+1) for k ∈ 0:12]
        println(BoolOps[n+1], " ", seq)
    end

    println()
    for n ∈ 1:8
        L280172(n) |> println
    end
end

function perf()
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
A327489, A327488, A327490, A001511,
A062383, A265705, A261693, A265716
=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


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

JLLWrappers.@generate_main_file_header("libdrm")
JLLWrappers.@generate_main_file("libdrm", UUID("8e53e030-5e6c-5a89-a30b-be5b7263a166"))
end  # module libdrm_jll

ext: jl


text:
println("\n\n\nPerformance of FinanceLib.Statements functions")

import FinanceLib.Statements as St

@code_warntype St.createCalcSets(St.balanceSheetCalcMap)
@code_warntype St.createCalcSets(St.profitLossCalcMap)

println(St.balanceSheetCalcItems)
println(St.balanceSheetEntries)

println(St.profitLossCalcItems)
println(St.profitLossEntries)

println(St.BsDict(St.Cash => 5.0), St.CfDict(St.ChangeInventories => 2.5))

println("Benchmarks\n")

abra(x::T) where T <: Union{Int64, Float64} =  2x

function abra(x)
  "Not here"

end

@code_warntype abra(5)

println(abra("a"))


text:
export @replay
import Interact, Reactive, DataFrames

const DEBUG = true

function Reactive.foldp(f::Function, inputs::Signal...; typ=Void)
    v0 = f(nothing,  map(value, inputs)...)
    typ == Void && (typ = typeof(v0))
    foldp(f, v0, inputs...; typ=typ)
end

"""
`traces4symvals(symvals::Dict; loc="")`

`symvals` is a Dict of symbols => their value at the desired
trace line hit count (lcount).

Returns:
a vector of TraceItems corresponding to the `@snapall` where all syms had their
respective vals
"""
function traces4symvals(symvals::Dict{String}; loc="", traces=happysnaps)
    query = Dict{Symbol, Any}()
    loc != "" && (query[:location] = loc)
    tis = traceitems(query)
    lcount = 0
    startidx = 0
    endidx = 0
    matched = fill(false, length(symvals))
    poss_match = true
    got_match = false
    for (i, ti) in enumerate(tis)
        if ti.lcount != lcount #XXX a more general way to implement this which would
            #enable `@snap`s in other code locations to be picked up, would be to
            #check for changes in any of the iterator values (symvals) rather than
            #changes in lcount. Implement using a cur_itervals Dict(exprstr => val).
            #empty! it if length(cur_itervals) == length(symvals) &&
            # !all(values(symvals) .== values(cur_itervals)), that way it will always have the seen values for the current iteration.
            if poss_match && all(matched)
                got_match = true
                endidx = i - 1
            elseif got_match
                #previous lcount were matching, now no longer matching, we're done
                break
            else
                startidx = i
                poss_match = true
            end
            fill!(matched, false)
            lcount = ti.lcount
        end
        !poss_match && continue
        for (i, (sym, val)) in enumerate(symvals)
            if ti.exprstr == sym
                #item is for the current `sym`
                if ti.val == val
                    matched[i] = true
                else
                    matched[i] = false
                    poss_match = false
                end
                break
            end
        end
    end
    # handle all matched on the last iteration, slightly more efficient than
    # to move the check to the bottom of the loop and change the i-1 to i
    # since the check currently only happens when the lcount changes (not every outer loop iteration)
    poss_match && all(matched) && (endidx = length(tis))
    tis[startidx:endidx]
end

traces4symvals(symvals::Pair...; loc="") = traces4symvals(Dict(symvals); loc=loc)

"""
`vals4traces(syms::Vector{String}, traces; loc="")`

Given a vector of variables `syms` and a list of traces, returns a vector of the
first value found in `traces` for each sym.
"""
function vals4traces(syms::Vector{String}, traces; loc="")
    vals = Array(Any, length(syms))
    for ti in traces
        i = findfirst(syms, ti.exprstr)
        if i != 0
            vals[i] = ti.val
        end
    end
    vals
end

"""
Values of `itervars[i]` assumed to be `Int`s and for fixed values of itervars[1:i-1]
they should be increasing (as index into `traces` increases). Also assumed is that
the nth trace of itervars[j] will always appear before the nth trace of itervars[j+1]
for all j. I.e. the iter vars are watched from slowest moving (outer-most loop iterator)
to fastest (inner-most loop iterator).
"""
function get_iter_extrema(itervars; traces=RickTracy.happysnaps)
    N = length(itervars)
    ivarstrs = string.(itervars)
    iter2idx = Dict(zip(ivarstrs, 1:N))
    extreme_default = (typemax(Int), typemin(Int)) #min, max worst case scenarios
    iter_extrema = [Dict{NTuple{i-1, Int}, Tuple{Int,Int}}() for i in 1:N]
    tpl_keys = zeros(Int, N)
    minpos, maxpos = (1,2)
    for (tidx, ti) in enumerate(traces)
        !(ti.exprstr in ivarstrs) && continue
        ivar = ti.exprstr
        val = ti.val
        i = iter2idx[ivar]
        tpl_keys[i] = ti.val
        tplkey = (tpl_keys[1:i-1]...)
        imin, imax = get!(iter_extrema[i], tplkey, extreme_default)
        new_min_or_max = false
        if val < imin
            imin = val
            new_min_or_max = true
        end
        if val > imax
            imax = val
            new_min_or_max = true
        end
        new_min_or_max && (iter_extrema[i][tplkey] = (imin, imax))
    end
    iter_extrema
end

function getsliders(itervars)
    iterator_extrema = get_iter_extrema(itervars)
    sliders = Array(Interact.Slider, length(itervars))
    slsigs  = Array(Reactive.Signal, length(itervars))
    minidx, maxidx = 1,2
    # @show itervars
    for (i, ivar) in enumerate(itervars)
        minmax = iterator_extrema[i]
        prev_itervals = (map(s->s.value, sliders[1:i-1])...)
        initial_range = UnitRange(minmax[prev_itervals]...)
        init_val = minmax[prev_itervals][minidx]
        sliders[i] = slider(initial_range;
                            value=init_val, label=string(ivar))
        slsigs[i] = signal(sliders[i])
        if i > 1
            #keep the ranges accurate on sliders, dependent on prev iterator values
            prevsigs = slsigs[1:i-1]
            # @show i prevsigs
            range_sig = Reactive.map(prevsigs...) do ivarvals...
                #unfortunately we can't use ivarvals, because changing the previous
                #slider's range will only push a value to slider's sig asynchronously
                #so ivarvals are old signal vals (whose signals are about to update)
                prev_itervals = (map(s->s.value, sliders[1:i-1])...)
                # @show i ivarvals prev_itervals minmax[prev_itervals]
                # haskey(minmax, (ivarvals...)) && (@show minmax[(ivarvals...)])
                (sliders[i].value < minmax[prev_itervals][minidx]) &&
                    set!(sliders[i], :value, minmax[prev_itervals][minidx])
                (sliders[i].value > minmax[prev_itervals][maxidx]) &&
                    set!(sliders[i], :value, minmax[prev_itervals][maxidx])
                UnitRange(minmax[prev_itervals]...)
            end |> Reactive.preserve
            set!(sliders[i], :range, range_sig)
        end
    end
    sliders, slsigs
end

function kwargs2assignment_block!(kw_ex)
    kw_ex.head = :block #:parameters to block
    kw_symvals = map(kw_ex.args) do ex
        ex.head = Symbol("=")
        ex.args[1], ex.args[2]
    end
    kw_syms, kw_vals = zip(kw_symvals...) |> collect #unzip
    kw_ex, [kw_syms...], [kw_vals...]
end

macro replay(expr)
    @eval begin using Interact, DataFrames end
    fnargs = expr.args[1].args
    numwsyms = 0
    widget_bindings = wsigs_ex = widgs_ex = :()
    if isa(fnargs[1], Expr)
        #kwargs
        widget_bindings, wsyms, widgets = kwargs2assignment_block!(fnargs[1])
        deleteat!(fnargs, 1)
        foreach(wsym->push!(fnargs, wsym), wsyms) #add widgsym vals as arguments
        wsigs_ex = :()
        foreach(wsym->push!(wsigs_ex.args, :(signal($wsym))), wsyms)
        widgs_ex = :()
        foreach(wsym->push!(widgs_ex.args, :($wsym)), wsyms)
        numwsyms = length(wsyms)
    end
    itervars = expr.args[1].args[1:(end-numwsyms)]
    varstrs = String.(unique(RickTracy.tracesdf(), :exprstr)[:exprstr])
    filter!(varstrs) do vstr; !(vstr in string.(itervars)) end
    varsyms = Symbol.(varstrs)
    block = expr.args[2].args[2] #function body
    itervalpairs = :()
    foreach(itervars) do sym
        push!(itervalpairs.args, :($(string(sym)) => $sym))
    end
    unshift!(fnargs, :prev) #for the foldp to return the prev valid value if we run into trouble
    unshift!(block.args,
        # get the traces where the iterators had their particular (current) values
        # and the values of the other syms at those traces
        quote
            _traces = RickTracy.traces4symvals($(itervalpairs)...)
            if isempty(_traces)
                if RickTracy.DEBUG
                    println(STDERR, "traces are empty for ", $itervalpairs)
                end
                return prev
            end
            _symvals = RickTracy.vals4traces($varstrs, _traces)
        end)

    # assign them to local vars with the same name ^_^
    for (i,sym) in enumerate(varsyms)
        symstr = string(sym)
        insert!(block.args, 2,
            quote
                local $sym
                try
                    $sym = _symvals[$i]
                catch e
                    if RickTracy.DEBUG
                        println(STDERR, "Failed to set ", $symstr, ", _symvals: ", _symvals, ", i: ", $i, "_traces: ", _traces)
                    end
                end
            end
        )
    end
    res = quote
        @eval begin using Interact, Reactive, DataFrames end
        sliders, slsigs = RickTracy.getsliders($itervars)
        $widget_bindings
        widglayout = hbox(vbox(sliders...), vbox($widgs_ex...))
        display(widglayout)
        Reactive.foldp($expr,
            slsigs...,
            $(wsigs_ex)...; typ=Any) |> Reactive.preserve
    end |> esc
    # dump(res, maxdepth = 30)
    res
end

max_stars_repo_name: JobJob/RickTracy


text:
using QuasiGeostrophy, FFTW, Test, BenchmarkTools
using LinearAlgebra
using Plots, GraphRecipes
using SymbolicUtils
import SymbolicUtils: Chain, Postwalk
include(pwd() * "/test/test_utils.jl")
const tol = 1e1
gr(size = (300,300))
# Define 1D Test
Ω = S¹(0, 2π)
Nx = Ny = Nz = 2^6;
fourier_grid = FourierGrid(Nx, Ω)
fieldnames = ("u", "σ", "v", "δu")
u, σ, v, δu = create_fields(names = fieldnames, grid = fourier_grid)
∂x = create_operators(fourier_grid)
Δ = ∂x^2
# initialize FourierField data
x = fourier_grid.grid[1]
Δx = x[2] - x[1]
u(sin.(x) .+ 0.5); σ(cos.(x));
v(sin.(x));
starting_energy = (u*u).data[1]

# initialize parameters
const Δt = 0.1 / Nx
const ν = 0.0/(2π)^2 * Δx^2 / Δt

# Forward-Euler
# uⁿ⁺¹ = uⁿ + Δt u̇ⁿ
T = 1.5
steps = round(Int, T / Δt)
plotsteps = round(Int, steps / 10)
for i in 1:steps
    u̇ = - ∂x( 0.5 * u * u) + ν * Δ(u)
    u.data .+= Δt * u̇.data
end
##

# Backward-Euler
# uⁿ⁺¹ = uⁿ + Δt u̇ⁿ⁺¹
# solve uⁿ⁺¹ - Δt *( ∂x( 0.5 * uⁿ⁺¹ * uⁿ⁺¹) - νΔuⁿ⁺¹ ) = uⁿ
# Jacobian
function Jacobian(u, δu, ν, Δt)
    return δu - Δt *( ∂x(u * δu) - ν*Δ(δu) )
end
# Residual
function Residual(uⁿ⁺¹, uⁿ, ν, Δt)
    u̇ⁿ⁺¹ = - ∂x( 0.5 * uⁿ⁺¹ * uⁿ⁺¹) + ν * Δ(uⁿ⁺¹)
    return uⁿ-uⁿ⁺¹ + Δt * u̇ⁿ⁺¹
end

# Create Jacobian Matrix for "Exact" Jacobian
function create_matrix(J, δu)
    n = length(δu.data)
    A = zeros(eltype(δu.data),(n,n))
    δu.data .= 0
    for i in 1:n
        δu.data[i] = 1.0
        jA = J(δu).data
        A[i,:] .+= jA
        δu.data[i] = 0.0
    end
    return A
end

##
u(sin.(x) .+ 0.5); σ(cos.(x));
v(sin.(x));

const Δt = 0.1 / Nx
const ν = 0.01/(2π)^2 * Δx^2 / Δt

# Backward-Euler
# uⁿ⁺¹ = uⁿ + Δt u̇ⁿ⁺¹
T = 1.5
steps = round(Int, T / Δt)
plotsteps = round(Int, steps / 10)

for i in 1:steps
    v.data .= u.data
    # initial guess for δu
    δu.data .= Δt * (- ∂x( 0.5 * u * u) + ν * Δ(u)).data
    #δu.data .= 0.0
    u.data .+= δu.data
    # newton iteration
    newtonN = 1
    for j in 1:newtonN
        residual = Residual(u, v, ν, Δt)
        #println("residual at ", i, " is ", norm(residual))
        J(δu) = Jacobian(u.data, δu, ν, Δt)
        A = create_matrix(J, δu)
        δu.data .= A \ residual.data
        u.data .+= δu.data
    end
    residual = Residual(u, v, ν, Δt)
    println("residual at ", i, " now is ", norm(residual))
end

## Semiplicit
# solve uⁿ⁺¹ - Δt *( -∂x( 0.5 * uⁿ * uⁿ⁺¹) + νΔuⁿ⁺¹ ) = uⁿ
# Semi-Implicit
const Δt = 0.1 / Nx
const ν = 0.0/(2π)^2 * Δx^2 / Δt

T = 1.5
steps = round(Int, T / Δt)
plotsteps = round(Int, steps / 10)

function semiplicit(u, δu, ν, Δt)
    v = u + Δt * (- ∂x( 0.5 * u * u) + ν * Δ(u))
    return δu + Δt *( 0.5*∂x(v * δu) - ν*Δ(δu) )
end

u(sin.(x) .+ 0.5); σ(cos.(x));
v(sin.(x));

for i in 1:steps
    v.data .= u.data
    Snew(δu) = semiplicit(u, δu, ν, Δt)
    A = create_matrix(Snew, δu)
    u.data .= A \ v.data
    println(i)
end

## Backward-Euler with semiplicit guess
# uⁿ⁺¹ = uⁿ + Δt u̇ⁿ⁺¹
u(sin.(x) .+ 0.5); σ(cos.(x));
v(sin.(x));
T = 1.5
const Δt = 0.1 / Nx
const ν = 0.0/(2π)^2 * Δx^2 / Δt
steps = round(Int, T / Δt)
plotsteps = round(Int, steps / 10)

for i in 1:steps
    v.data .= u.data
    # initial guess for δu
    Snew(δu) = semiplicit(u, δu, ν, Δt)
    A = create_matrix(Snew, δu)
    u.data .=  A \ v.data
    # newton iteration
    newtonN = 1
    for j in 1:newtonN
        residual = Residual(u, v, ν, Δt)
        #println("residual at ", i, " is ", norm(residual))
        J(δu) = Jacobian(u.data, δu, ν, Δt)
        A = create_matrix(J, δu)
        δu.data .= A \ residual.data
        u.data .+= δu.data
    end
    residual = Residual(u, v, ν, Δt)
    println("residual at ", i, " now is ", norm(residual))
end

max_stars_repo_name: sandreza/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: sandreza/QuasiGeostrophy


text:
@testset "Basic Functionality" begin

    A = rand(100,100)
    img = DirectImage(A)
    img[1,1]

    img = DirectImage(A, -50:50, -50:50)
    img[0,0]

    img = DirectImage(A, -25:0.5:25, -25:0.5:25)
    img[1.5,1.5]

    img = DirectImage(A, -25m:0.5m:25m, -25m:0.5m:25m)
    img[9.5m,20m]

    img = DirectImage(view(A, 1:10, 1:10), 1:10, 1:10)
    img[5,5]

    img = DirectImage(A, c=1, d=2, f=5)
    img.c == 1

    img = DirectImage(A, f = 1 => "a comment", b = 5 => "another comment")
    img.f == 1
    img[/].f == "a comment"

    # comments must be string, values must be union{number,string,}

    @testset "IO" begin

        img = DirectImage("test.fits")
        img = DirectImage("test.fits", 1, 1:100, 2:5)

        write("test2.fits", img)
    end

    # Need to handle multi-dimensional cubes
    # Want to also be able to stack multiple arrays or views into one thing that looks like
    # a cube.

end

max_stars_repo_name: sefffal/SpaceImages


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Compute the Greenwich Mean Sideral Time (GMST).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] http://www.navipedia.net/index.php/CEP_to_ITRF, accessed 2015-12-01.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export j2000_to_gmst, jd_to_gmst

"""
    j2000_to_gmst(j2000_ut1::Number)

Compute the Greenwich Mean Sideral Time (GMST) \\[rad] given the instant
`j2000_ut1` in J2000.0 reference [UT1].

!!! info
    The algorithm is based in **[1]**.

# References

- **[1]** http://www.navipedia.net/index.php/CEP_to_ITRF, accessed 2015-12-01.
"""
function j2000_to_gmst(j2000_ut1::Number)
	# Julian centuries elapsed from the epoch J2000.0.
	t_ut1 = j2000_ut1 / 36525

	# Greenwich Mean Sideral Time at t_ut1 [s].
    θ_gmst = @evalpoly(
        t_ut1,
        + 67310.54841,
        + (876600.0 * 3600 + 8640184.812866),
        + 0.093104,
        - 6.2e-6
    )

    # Reduce to the interval [0, 86400]s.
    θ_gmst = mod(θ_gmst, 86400)

    # Convert to radian and return.
    return θ_gmst * π / 43200
end

"""
    jd_to_gmst(jd_ut1::Number)

Compute the Greenwich Mean Sideral Time (GMST) \\[rad] for the Julian Day
`jd_ut1` [UT1].

!!! info
    The algorithm is based in **[1]**(p. 188).

# References

- **[1]** Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
jd_to_gmst(jd_ut1::Number) = j2000_to_gmst(jd_ut1 - JD_J2000)


text:
import PSLQ: make_H, make_DH, pslq_step, check_H, check_square, findmaxj, lqrefresh, integertype

function example_1(a, n::Integer)
    x = Vector{eltype(a)}(undef, n)
    b = copy(a)
    for i = 1:n-1
        x[i] = b
        b *= a
    end
    x[n] = sum(x[1:n-1] .* (1:n-1))
    x
end

function example_2(a, n::Integer)
    x = Vector{eltype(a)}(undef, n)
    b = oftype(a, 1)
    for i = 1:n
        x[i] = b
        b *= a
    end
    x
end

function make_Htest(n::Integer; gamma=1, T=Float64)
    H = zeros(T, n, n-1)
    for i = 1:n-1; H[i,i] = gamma^(1-i); end
    for i = 1:n, j = 1:i-1
        H[i,j] = i + j*j
    end
    H
end

@testset "correctness of calculations type $T" for T in (Float64, BigFloat)

    H = make_Htest(1, T=T)
    @test_throws DimensionMismatch check_H(H)
    @test_throws DimensionMismatch check_square(H, 1)

    @test integertype(Float32) == Int
    @test integertype(Float64) == Int
    @test integertype(BigFloat) == BigInt

    H = make_Htest(10, gamma=2.0, T=T)
    @test findmaxj(H, 2.0) == 9
    H[1,1] = nextfloat(H[1,1])
    @test findmaxj(H, 2.0) == 1
    H[7,7] = nextfloat(H[7,7])
    @test findmaxj(H, 2.0) == 7

    H = make_Htest(10, T=T)
    D, DH = make_DH(H)
    @test norm(D * H - Matrix(I, size(H))) == 0

    x = example_1(T(ℯ), 6)
    H = make_H(x)
    P = I - x * x' / norm(x)^2
    Hn = norm(H)
    @test abs(Hn - sqrt(T(size(H,2)))) <= 2eps(Hn)
    @test norm(H'H - I) <= 2eps(Hn)
    @test norm(H * H' - P) <= 2eps(Hn)

    D, DH = make_DH(H)
    @test diag(D) == ones(eltype(D), size(D,2))
    @test istril(D)

    @test_throws DimensionMismatch pslq_step(zeros(4,3), zeros(Int,4,3), zeros(Int,4,3))
    @test_throws DimensionMismatch pslq_step(zeros(4,2), zeros(Int,4,4), zeros(Int,4,4))

    A1 = Matrix(I*1, size(H, 1), size(H,1))
    B1 = copy(A1)
    H1 = H
    H1, A1, B1 = pslq_step(H1, A1, B1)
    H1, A1, B1 = pslq_step(H1, A1, B1)
    H1, A1, B1 = pslq_step(H1, A1, B1)
    H1, A1, B1 = pslq_step(H1, A1, B1)
    @test A1 * B1 == I
    X = H \ (B1 * H1) # verify that H = B1*H1 * X with orthonormal X
    @test norm(X'X - I) <= 10eps(Hn)
end

@testset "iteration algorithm with type $T" for T in (Float64,BigFloat)

    x = example_2(T(2.0)^(T(1)/7), 8)
    r = pslq(x)
    @test r.code == 1
    @test r.iter >= 30
    @test r.col == 1
    @test norm((r.B'x)[r.col]) <= 10*eps(norm(x))

    H = r.H
    H0 = make_H(x)
    H2 = lqrefresh(r.A, H0)
    @test norm(H2 - H) <= 1e8 * eps(norm(x))
end


text:
# Atomic opterations

const atomic_ops = ((:add!, :+), (:sub!, :-), (:min!, :min), (:max!, :max),
                    (:and!, :&), (:or!, :|), (:xor!, :xor))

# Methods for writing to a WritableGridData grid from . These are
# associative and commutative so that write order does not affect the result.
for (f, op) in atomic_ops
    @eval begin
        @propagate_inbounds ($f)(d::AbstractSimData, x, I...) = ($f)(first(d), x, I...)
        @propagate_inbounds ($f)(d::WritableGridData{<:Any,R}, x, I...) where R = ($f)(d, proc(d), x, I...)
        @propagate_inbounds function ($f)(d::WritableGridData{<:Any,R}, ::Processor, x, I...) where R
            @boundscheck checkbounds(dest(d), I...)
            @inbounds _setdeststatus!(d, x, I...)
            @inbounds dest(d)[I...] = ($op)(dest(d)[I...], x)
        end
        @propagate_inbounds function ($f)(
            d::WritableGridData{<:Any,R}, proc::Union{ThreadedCPU,CPUGPU}, x, I...
        ) where R
            lock(proc)
            ($f)(d, SingleCPU(), x, I...)
            unlock(proc)
        end
    end
end


text:
# Linear Isotropic Covariance Function

"""
    LinIso <: Kernel

Isotropic linear kernel (covariance)
```math
k(x, x') = xᵀx'/ℓ²
```
with length scale ``ℓ``.
"""
mutable struct LinIso <: Kernel
    "Squared length scale"
    ℓ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        LinIso(ll::Float64)

    Create `LinIso` with length scale `exp(ll)`.
    """
    LinIso(ll::Float64) = new(exp(2 * ll), [])
end

struct LinIsoData{D} <: KernelData
    XtX::D
end

function KernelData(k::LinIso, X::MatF64)
    XtX=X'*X
    LinearAlgebra.copytri!(XtX, 'U') # make sure it's symmetric
    LinIsoData(XtX)
end
kernel_data_key(k::LinIso, X::MatF64) = "LinIsoData"

_cov(lin::LinIso, xTy) = xTy ./ lin.ℓ2
function Statistics.cov(lin::LinIso, x::VecF64, y::VecF64)
    K = _cov(lin, dot(x,y))
    return K
end

@inline @inbounds function cov_ij(lin::LinIso, X::MatF64, data::LinIsoData, i::Int, j::Int, dim::Int)
    return _cov(lin, data.XtX[i, j])
end
function Statistics.cov(lin::LinIso, X::MatF64, data::LinIsoData)
    K = _cov(lin, data.XtX)
    return K
end
function cov!(cK::MatF64, lin::LinIso, X::MatF64, data::LinIsoData)
    iℓ2 = 1/lin.ℓ2
    @inbounds @simd for I in eachindex(cK,data.XtX)
        cK[I] = data.XtX[I]*iℓ2
    end
    return cK
end

get_params(lin::LinIso) = Float64[log(lin.ℓ2) / 2]
get_param_names(::LinIso) = [:ll]
num_params(lin::LinIso) = 1

function set_params!(lin::LinIso, hyp::VecF64)
    length(hyp) == 1 || throw(ArgumentError("Linear isotropic kernel only has one parameter"))
    lin.ℓ2 = exp(2 * hyp[1])
end

@inline dk_dll(lin::LinIso, xTy::Float64) = -2 * _cov(lin,xTy)
@inline function dKij_dθp(lin::LinIso, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    if p==1
        return dk_dll(lin, dotij(X,i,j,dim))
    else
        return NaN
    end
end
@inline function dKij_dθp(lin::LinIso, X::MatF64, data::LinIsoData, i::Int, j::Int, p::Int, dim::Int)
    if p==1
        return dk_dll(lin, data.XtX[i,j])
    else
        return NaN
    end
end


text:
function sdfcount(root_folder, pattern)
    folders = filter!(isdir, readdir(root_folder, join=true))
    filter!(f-> occursin(pattern, f), folders)
    n = [count(f->endswith(f, ".sdf"), readdir(dir)) for dir in folders]
    idx = findall(!iszero, n)

    Plots.plot(basename.(folders[idx]), n[idx], rotation=45, legend=false)
end

function similar_laser(laser::Type{LaguerreGaussLaser}, t_profile, file)
    fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l
    λ = get_parameter(file, :laser, :lambda)
    m = get_parameter(file, :constant, :m)
    a₀ = float(get_parameter(file, :constant, :a0))
    w₀ = get_parameter(file, :constant, :w_0)u"m" |> unit_l
    τ = get_parameter(file, :constant, :tau_0)u"m" |> unit_t
    t = get_time(file) |> unit_t

    s = setup_laser(laser, :SI_unitful;
        λ,
        m,
        p = 0,
        a₀,
        w₀,
        # ϕ₀=π/2,   # aparently this doesn't do anything
        profile = t_profile(;
            c = c_0,
            τ,
            z₀ = -fx,
            t₀ = fx/(π*c_0) |> unit_t
        )
    )
end

function field_rotation(dir)
    if dir == :x
        RotZ(π)*RotY(-π/2)*RotZ(-π/2)
    elseif dir == :y
        RotY(-π/2)
    else
        1
    end
end

function similar_E(::Val{true}, laser, t, propagation_dir, component::Int, grid; x₀, y₀, z₀)
    map(Iterators.product(grid...)) do (x,y,z)
        r = SVector{3}(x-x₀, y-y₀, z-z₀)
        R = field_rotation(propagation_dir)
        analytic_E = R * (LaserTypes.E(R \ r, laser) * LaserTypes.g((R \ r)[3], t, laser))
        analytic_E[component]
    end
end

function similar_E(::Val{false}, laser, t, propagation_dir, component::Int, grid; x₀, y₀, z₀)
    map(Iterators.product(grid...)) do (x,y,z)
        r = SVector{3}(x-x₀, y-y₀, z-z₀)
        R = field_rotation(propagation_dir)
        analytic_E = R * LaserTypes.E(R \ r, t, laser)
        analytic_E[component]
    end
end

function similar_E(E, t, laser;
    component=1,
    propagation_dir=:x,
    complex=false,
    x₀=zero(recursive_bottom_eltype(getdomain(E))),
    y₀=zero(recursive_bottom_eltype(getdomain(E))),
    z₀=zero(recursive_bottom_eltype(getdomain(E))))

    E_approx = approximate_field(E)

    grid = AxisGrid((getdomain(E_approx).*unit_l...,))
    data = similar_E(Val(complex), laser, t, propagation_dir, component, grid; x₀,y₀,z₀)

    ScalarField(data, grid)
end

function compare_E_slice_with_analytic(laser, profile, file, dir;
    slice_location=0unit_l,
    slice_dir=:x)

    fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l
    t = get_time(file) |> unit_t
    s = similar_laser(laser, profile, file)

    E = uconvert(unit_E, unit_l, file["e$dir"])

    analytic_E = similar_E(E, t, s,
        component=dir_to_idx(dir),
        propagation_dir=:x,
        x₀=fx) |> unit_E

    plt1 = E_slice_plot(analytic_E, dir, slice_location, t;
        legend_title = "Analytic",
        slice_dir,
        issliced=false)
    plt2 = E_slice_plot(E, dir, slice_location, t;
        legend_title = "Numeric",
        slice_dir,
        issliced=false)
    Plots.plot(plt1, plt2, size=(400,600), layout=(2,1))
end


text:
N = 10_000
rngs = Array(MersenneTwister, N)
rngs2 = Array(MersenneTwister, N)

for i in 1:N
    rngs[i] = MersenneTwister(i)
end

for i in 1:N
    srand(rngs[i], i)
end

for i in 1:N
    rngs2[i] = copy(rngs[i])
end

for i in 1:N
    copy!(rngs2[i], rngs[i])
end

println("create new")
@time for i in 1:N
    rngs[i] = MersenneTwister(i)
end

println("srand")
@time for i in 1:N
    srand(rngs[i], i)
end

println("copy")
@time for i in 1:N
    rngs2[i] = copy(rngs[i])
end

println("copy!")
@time for i in 1:N
    copy!(rngs2[i], rngs[i])
end


text:
module AtomicLevels

using UnicodeFun
using Formatting
using Parameters
using BlockBandedMatrices
using WignerSymbols
using HalfIntegers
using Combinatorics

include("common.jl")
include("unicode.jl")
include("parity.jl")
include("orbitals.jl")
include("relativistic_orbitals.jl")
include("spin_orbitals.jl")
include("configurations.jl")
include("excited_configurations.jl")
include("terms.jl")
include("allchoices.jl")
include("jj_terms.jl")
include("intermediate_terms.jl")
include("couple_terms.jl")
include("csfs.jl")
include("jj2lsj.jl")
include("levels.jl")

module Utils
include("utils/print_states.jl")
end

# Deprecations
@deprecate jj2lsj(args...) jj2ℓsj(args...)

end # module


text:
@testset "Objective" begin
    T = 3
    nx = 2
    nu = 1
    nw = 0
    ot = (x, u, w) -> dot(x, x) + 0.1 * dot(u, u)
    oT = (x, u, w) -> 10.0 * dot(x, x)
    ct = Cost(ot, nx, nu, nw, [t for t = 1:T-1])
    cT = Cost(oT, nx, 0, nw, [T])
    obj = [ct, cT]

    J = [0.0]
    grad = zeros((T - 1) * (nx + nu) + nx)
    idx_xu = [collect((t - 1) * (nx + nu) .+ (1:(nx + (t == T ? 0 : nu)))) for t = 1:T]
    x1 = ones(nx)
    u1 = ones(nu)
    w1 = zeros(nw)
    X = [x1 for t = 1:T]
    U = [t < T ? u1 : zeros(0) for t = 1:T]
    W = [w1 for t = 1:T]

    ct.val(ct.val_cache, x1, u1, w1)
    ct.grad(ct.grad_cache, x1, u1, w1)
    @test ct.val_cache[1] ≈ ot(x1, u1, w1)
    @test norm(ct.grad_cache - [2.0 * x1; 0.2 * u1]) < 1.0e-8

    cT.val(cT.val_cache, x1, u1, w1)
    cT.grad(cT.grad_cache, x1, u1, w1)
    @test cT.val_cache[1] ≈ oT(x1, u1, w1)
    @test norm(cT.grad_cache - 20.0 * x1) < 1.0e-8

    @test DirectTrajectoryOptimization.eval_obj(obj, X, U, X) - sum([ot(X[t], U[t], W[t]) for t = 1:T-1]) - oT(X[T], U[T], W[T]) ≈ 0.0
    DirectTrajectoryOptimization.eval_obj_grad!(grad, idx_xu, obj, X, U, W)
    @test norm(grad - vcat([[2.0 * x1; 0.2 * u1] for t = 1:T-1]..., 20.0 * x1)) < 1.0e-8

    # info = @benchmark DirectTrajectoryOptimization.eval_obj($obj, $X, $U, $W)
    # info = @benchmark DirectTrajectoryOptimization.eval_obj_grad!($grad, $idx_xu, $obj, $X, $U, $W)
end


text:
module ContactImplicitMPC

using BenchmarkTools
using Colors
using CoordinateTransformations
using FFMPEG
using FileIO
using ForwardDiff
using GeometryBasics
using IfElse
using InteractiveUtils
using JLD2
using LinearAlgebra
using Logging
using QDLDL
using MeshCat
using MeshCatMechanisms
using Meshing
using MeshIO
using Parameters
using Plots
using Random
using Rotations
using SparseArrays
using StaticArrays
using SuiteSparse
using Symbolics
using Scratch
import Scratch: get_scratch!
using Test
using RoboDojo
import RoboDojo: LinearSolver, LUSolver, Model, ResidualMethods, Space, Disturbances, IndicesZ, InteriorPoint, EmptySolver, Policy, Trajectory, GradientTrajectory, InteriorPointOptions, IndicesOptimization, interior_point, interior_point_solve!, bilinear_violation, residual_violation, general_correction_term!, r!, rz!, rθ!, linear_solve!, lu_solver, empty_policy, empty_disturbances, friction_coefficients, SimulatorStatistics, SimulatorOptions, indices_θ, num_data, initialize_z!, initialize_θ!, indices_z, indices_θ, simulate!, policy, process!, Simulator
using DirectTrajectoryOptimization

# Utilities
include("utils.jl")

# Solver
include("solver/lu.jl") # sparse arrays
include("solver/ldl.jl")

include("solver/qr.jl")
include("solver/schur.jl")

# Environment
include("simulator/environment.jl")

# Dynamics
include("dynamics/model.jl")

# Simulator
include("simulation/index.jl")

# Simulator
include("simulation/contact_methods.jl")
include("simulation/simulation.jl")

include("dynamics/code_gen_dynamics.jl")
include("dynamics/fast_methods_dynamics.jl")

# Models
include("dynamics/quaternions.jl")
include("dynamics/mrp.jl")
include("dynamics/euler.jl")

include("dynamics/particle_2D/model.jl")
include("dynamics/particle/model.jl")
include("dynamics/hopper_2D/model.jl")
include("dynamics/hopper_3D/model.jl")
include("dynamics/quadruped/model.jl")
include("dynamics/flamingo/model.jl")
include("dynamics/pushbot/model.jl")
include("dynamics/walledcartpole/model.jl")
include("dynamics/centroidal_quadruped/model.jl")
include("dynamics/point_foot_quadruped/model.jl")

# Simulator
include("simulator/policy.jl")

include("simulator/disturbances.jl")
include("simulator/simulator.jl")

# Simulation
include("simulation/environments/flat.jl")
include("simulation/environments/piecewise.jl")
include("simulation/environments/quadratic.jl")
include("simulation/environments/slope.jl")
include("simulation/environments/sinusoidal.jl")
include("simulation/environments/stairs.jl")

include("simulation/residual_approx.jl")
include("simulation/code_gen_simulation.jl")

# Controller
include("controller/trajectory.jl")
include("controller/linearized_step.jl")
include("controller/implicit_dynamics.jl")
include("controller/objective.jl")
include("controller/linearized_solver.jl")
include("controller/newton.jl")
include("controller/newton_indices.jl")
include("controller/newton_residual.jl")
include("controller/newton_jacobian.jl")
include("controller/mpc_utils.jl")
include("controller/policy.jl")
include("controller/newton_structure_solver/methods.jl")
include("controller/gains.jl")

# Visuals
include("dynamics/visuals.jl")
include("dynamics/visual_utils.jl")
include("visuals.jl")

include("dynamics/particle_2D/visuals.jl")
include("dynamics/particle/visuals.jl")
include("dynamics/hopper_2D/visuals.jl")
include("dynamics/hopper_3D/visuals.jl")
include("dynamics/quadruped/visuals.jl")
include("dynamics/flamingo/visuals.jl")
include("dynamics/pushbot/visuals.jl")
include("dynamics/walledcartpole/visuals.jl")
include("dynamics/centroidal_quadruped/visuals.jl")
include("dynamics/point_foot_quadruped/visuals.jl")

export
    initial_conditions,
    World,
    LinearizedCone,
    NonlinearCone,
    Model,
    Dimensions,
    BaseMethods,
    DynamicsMethods,
    ResidualMethods,
    Environment,
    environment_2D,
    environment_3D,
    environment_2D_flat,
    environment_3D_flat,
    get_model,
    SparseStructure,
    LinearizedStep,
    bil_addition!,
    r_linearized!,
    rz_linearized!,
    ImplicitTrajectory,
    linearization!,
    implicit_dynamics!,
    TrackingObjective,
    TrackingVelocityObjective,
    second_order_cone_product,
    generate_base_expressions,
    RLin,
    RZLin,
    RθLin,
    ContactTraj,
    Simulation,
    num_var,
    num_data,
    get_simulation,
    get_trajectory,
    interior_point,
    InteriorPointOptions,
    interior_point_solve!,
    r!,
    rz!,
    rθ,
    generate_base_expressions,
    save_expressions,
    instantiate_base!,
    generate_dynamics_expressions,
    save_expressions,
    instantiate_dynamics!,
    environment_3D_flat,
    friction_dim,
    dim,
    sqrt_quat,
    cayley_map,
    L_multiply,
    R_multiply,
    R2,
    R3,
    FrictionCone,
    rotation,
    module_dir,
    open_loop_disturbances,
    disturbances,
    Disturbances,
    NoDisturbances,
    OpenLoopDisturbance,
    impulse_disturbances,
    ImpulseDisturbance,
    RandomDisturbance,
    open_loop_policy,
    policy,
    ci_mpc_policy,
    NewtonOptions,
    CIMPCOptions,
    SimulatorOptions,
    simulator,
    simulate!,
    generate_residual_expressions,
    instantiate_residual!,
    ϕ_func,
    tracking_error,
    repeat_ref_traj,
    Schur,
    schur_factorize!,
    schur_solve!,
    LinearSolver,
    LUSolver,
    lu_solver,
    ldl_solver,
    LDLSolver,
    factorize!,
    linear_solve!,
    IndicesOptimization,
    index_q2,
    index_γ1,
    index_b1,
    index_ψ1,
    index_s1,
    index_η1,
    index_s2,
    index_q0,
    index_q1,
    index_u1,
    index_w1,
    index_μ,
    index_h,
    index_dyn,
    index_imp,
    index_mdp,
    index_fri,
    index_bimp,
    index_bmdp,
    index_bfri,
    linearization_var_index,
    linearization_term_index,
    index_ort,
    index_soc,
    num_var,
    num_data,
    num_bilinear,
    index_equr,
    index_ortr,
    index_socr,
    visualize_meshrobot!,
    visualize_robot!,
    visualize_force!,
    visualize_disturbance!,
    visualize_payload!,
    process!,
    contact_trajectory,
    pack_z,
    pack_θ,
    generate_pusher_traj,
    update_friction_coefficient!
    quadratic_objective,
    particle,
    particle_2D,
    hopper_2D,
    hopper_3D,
    quadruped,
    quadruped_payload,
    flamingo,
    pushbot,
    walledcartpole,
    centroidal_quadruped,
    flat_3D_lc,
    flat_3D_nc,
    quadratic_bowl_3D_lc,
    flat_2D_lc,
    flat_2D_nc,
    slope1_2D_lc,
    sine2_2D_lc,
    stairs3_2D_lc,
    flat_3D_lc,
    sine2_3D_lc,
    sine1_2D_lc,
    piecewise1_2D_lc,
    sine3_2D_lc,
    slope_smooth_2D_lc,
    flat_2D_lc,
    flat_3D_lc,
    plot_surface!,
    stairs!,
    Simulator,
    set_trajectory!,
    initialize_z!

end # module


text:
using CompressedSensing
using Base.Test

srand(0)

#Setup a signal and measurement for testing
Signal=zeros(100)
Signal[int(rand(20)*9+1)]=randn(20)
Measure = randn(40,100)
Noise=randn(40)*.1

#try to rebuild the signal
IRLS_GuessSignal = IRLS(Measure,[Measure*Signal],threshold=1e-10,maxiter=5000)

#see if answer is within 5 digits
@test_approx_eq_eps IRLS_GuessSignal Signal 1e-5
info("IRLS() PASSED")

UIRLS_GuessSignal=UIRLS(Measure,Measure*Signal+Noise,lambda=.1,maxiter=2000)

@test cor(UIRLS_GuessSignal, Signal)>.995
info("UIRLS() PASSED")

A = zeros(100)
A[1]=1

@test GI(A)==.99
info("GI() PASSED")

max_stars_repo_name: dahlend/CompressedSensing


text:
using Stats
using Base.Test

@test_approx_eq autocor([1, 2, 3, 4, 5], 1) 0.4

@test iqr([1, 2, 3, 4, 5]) == [2.0, 4.0]

z = [true, true, false, false, true, false, true, true, true]
values, lengths = rle(z)
@test values == [true, false, true, false, true]
@test lengths == [2, 2, 1, 1, 3]
@test inverse_rle(values, lengths) == z

z = [true, true, false, false, true, false, true, true, true, false]
values, lengths = rle(z)
@test values == [true, false, true, false, true, false]
@test lengths == [2, 2, 1, 1, 3, 1]
@test inverse_rle(values, lengths) == z

X = [1 0; 2 1; 3 0; 4 1; 5 10]
y = [5, 3, 4, 2, 5]
@test_approx_eq cor_spearman(X, y)[1] cor_spearman(X[:,1],y)
@test_approx_eq cor_spearman(X) cor_spearman(X, X)
@test_approx_eq cor_spearman(X, y) [-0.102597835208515, -0.081110710565381]
@test_approx_eq cor_kendall(X,y) [-0.105409255338946, -0.117851130197758]

fnecdf = ecdf(randn(10000000))
@test_approx_eq_eps fnecdf([-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]) [0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975] 1e-3
@test_approx_eq_eps fnecdf(1.96) 0.975 1e-3
@test_approx_eq fnecdf([-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]) map(fnecdf, [-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96])


text:
"""
Implements toy data from paper
"""
function generate_toy_data(N_train=50)

    X = fill(collect(0.005:0.005:1.), 3)
    Y = fill(Float64[], 3)

    Sl = [5, 1, 2]
    Pl = [200, 0.1, 100]

    conv_ker = (Δ, S, P) ->  S * exp(-P * Δ^2)
    u = t -> sum(k -> (1 / k^2) * cos(2 * k * π * t), 1:4)

    for i in 1:3
        y = fill(0., 200)
        for j in 1:200
            y[j] =  quadgk(τ -> conv_ker(X[i][j] - τ, Sl[i], Pl[i]) * u(τ), 0., X[i][j])[1]
        end
        y = sum(x -> y.^x, 1:3)
        Y[i] = y +  randn(200) * sqrt(0.005 * (sum((y .- (sum(y) / 200)).^2) / 200))
    end

    train_X = fill(Float64[], 3)
    test_X = fill(Float64[], 3)
    train_Y = fill(Float64[], 3)
    test_Y = fill(Float64[], 3)
    for i in 1:3
        mix = randperm(200)
        mixtr = sort(mix[1:N_train])
        mixte = sort(mix[N_train:end])

        train_X[i] = X[i][mixtr]
        test_X[i] = X[i][mixte]
        train_Y[i] = Y[i][mixtr]
        test_Y[i] = Y[i][mixte]
    end

    (Data(train_X, train_Y), Data(test_X, test_Y))
end

function load_weather(N_train)

end


text:
module MakiePublication

using CairoMakie
using ColorSchemes

include("defaults.jl")

include("utils.jl")
export
    figsize,
    savefig

include("schemes.jl")
export
    theme_acs,
    theme_acs_1col,
    theme_acs_2col,
    theme_aps,
    theme_aps_1col,
    theme_aps_2col,
    theme_rsc,
    theme_rsc_1col,
    theme_rsc_2col,
    theme_web

end # module


text:
"""
    logit10(x::Real)
Function transforming interval [0,1] to [-Inf, Inf] using logit transformation.
"""
function logit10(x::Real)
    log10(x / (1.0 - x))
end

"""
    logistic10(x::Real)
Function transforming interval [-Inf, Inf] to [0,1] using logistic transformation.
Inverse function for [`logit10`](@ref).
"""
function logistic10(x::Real)
    exp10(x) / (exp10(x) + 1.0)
end

"""
    scaling(x::Real, scale::Symbol = :direct)
Transforms values from specific scale to range [-Inf, Inf] based on option.

## Return
Transformed value.

## Arguments
* `x`: input value.
* `scale`: transformation type: `:direct` (`:lin`), `:log`, `:logit`.
"""
function scaling(x::Real, scale::Symbol = :direct)
    if scale == :direct || scale == :lin
        return x
    elseif scale == :log
        return log10(x)
    elseif scale == :logit
        return logit10(x)
    else
        throw(DomainError(scale, "scale type is not supported"))
    end
end

"""
    unscaling(x::Real, scale::Symbol = :direct)
Transforms values from [-Inf, Inf] to specific scale based on option. Inverse function
for [`scaling`](@ref).

## Return
Transformed value.

## Arguments
* `x`: input value.
* `scale`: transformation type: `:direct` (`:lin`), `:log`, `:logit`.
"""
function unscaling(x::Real, scale::Symbol = :direct)
    if scale == :direct || scale == :lin
        return x
    elseif scale == :log
        return exp10(x)
    elseif scale == :logit
        return logistic10(x)
    else
        throw(DomainError(scale, "scale type is not supported"))
    end
end

scaling(x_tup::Tuple{Float64,Float64}, scale::Symbol = :direct) = scaling.(x_tup, scale)

unscaling(x_tup::Tuple{Float64,Float64}, scale::Symbol = :direct) = unscaling.(x_tup, scale)

scaling(::Nothing, ::Symbol) = nothing

unscaling(::Nothing, ::Symbol) = nothing

"""
    function get_endpoint(
        theta_init::Vector{Float64},
        theta_num::Int,
        loss_func::Function,
        method::Symbol,
        direction::Symbol = :right;

        loss_crit::Float64 = 0.0,
        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_bound::Float64 = unscaling(
            (direction==:left) ? -9.0 : 9.0,
            scale[theta_num]
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 0.,
        local_alg::Symbol = :LN_NELDERMEAD,
        max_iter::Int = 10^5,
        loss_grad::Union{Function, Symbol} = :EMPTY,
        silent::Bool = false
    )

Calculates confidence interval's right or left endpoints for a given parameter `theta_num`.

## Return
[`EndPoint`](@ref) object storing confidence interval's endpoint and intermediate profile points.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values should not necessary be the optimum values of `loss_func` but `loss_func(theta_init)` should be lower than `loss_crit`.
- `theta_num`: index of vector component for identification: `theta_init(theta_num)`.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` for profile likelihood-based (PL) identification. Usually we use log-likelihood for PL analysis: ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to estimate confidence interval's endpoint. Currently the following methods are implemented: `:CICO_ONE_PASS`, `:LIN_EXTRAPOL`, `:QUADR_EXTRAPOL`.
- `direction`: `:right` or `:left` endpoint to estimate.

## Keyword arguments
- `loss_crit` : critical level of loss function. Confidence interval's endpoint value is the intersection point of profile likelihood and `loss_crit` level.
- `scale` : vector of scale transformations for each parameters' component. Possible values: `:direct` (`:lin`), `:log`, `:logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `theta_bounds` : vector of tuple `(lower_bound, upper_bound)` for each parameter. Bounds define the ranges for possible parameter values. Default bounds are `(-Inf,Inf)`.
- `scan_bound` : value which states the area of confidence point analysis.
- `scan_tol` : Absolute tolerance for `theta_num` parameter used as termination criterion.
- `loss_tol` : Absolute tolerance controlling `loss_func` closenes to `loss_crit` (termination criterion). Currently doesn't work for `:CICO_ONE_PASS` method because of limitation in `LN_AUGLAG` interface.
- `local_alg` : algorithm of optimization. Derivative-free and gradient-based algorithms form NLopt package.
- `max_iter` : maximal number of fitter iterations. If reaches the result status will be `:MAX_ITER_STOP`.
- `loss_grad` : For gradient optimization methods it is necessary to set how the gradient of `loss_func` should be calculated.
    There are options:
    - `:EMPTY` (default) no gradient is set. It works only for gradient-free methods.
    - `:AUTODIFF` means autodifferentiation from `ForwardDiff` package is used.
    - `:FINITE` means finite difference method from `Calculus` is used.
    - It is also possible to set gradient function here `function(x::Vector{Float64})` which returns gradient vector.
- `silent` : Boolean argument declaring whether we display the optimization progress. Default is `false`

See also [`get_interval`](@ref)
"""
function get_endpoint(
    theta_init::Vector{Float64},
    theta_num::Int,
    loss_func::Function,
    method::Symbol,
    direction::Symbol = :right;

    loss_crit::Float64 = 0.0,
    scale::Vector{Symbol} = fill(:direct, length(theta_init)), # :direct, :lin, :log, :logit
    theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
        fill((-Inf, Inf), length(theta_init)),
        scale
        ),
    scan_bound::Float64 = unscaling(
        (direction==:left) ? -9.0 : 9.0,
        scale[theta_num]
        ),
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 0.,
    local_alg::Symbol = :LN_NELDERMEAD,
    max_iter::Int = 10^5,
    loss_grad::Union{Function, Symbol} = :EMPTY,
    silent::Bool = false,
    #kwargs... # other options for get_right_endpoint
)
    isLeft = direction == :left
    n_theta = length(theta_init)

    # checking arguments
    # theta_bound[1] < theta_init < theta_bound[2]
    theta_init_outside_theta_bounds = .! [theta_bounds[i][1] < theta_init[i] < theta_bounds[i][2] for i in 1:length(theta_init)]
    if any(theta_init_outside_theta_bounds)
        throw(ArgumentError("theta_init is outside theta_bound: $(findall(theta_init_outside_theta_bounds))"))
    end
    # scan_bound should be within theta_bounds
    !(theta_bounds[theta_num][1] < scan_bound < theta_bounds[theta_num][2]) &&
        throw(ArgumentError("scan_bound are outside of the theta_bounds $(theta_bounds[theta_num])"))
    # theta_init should be within scan_bound
    if (theta_init[theta_num] >= scan_bound && !isLeft) || (theta_init[theta_num] <= scan_bound && isLeft)
        throw(ArgumentError("init values are outside of the scan_bound $scan_bound"))
    end
    # 0 <= theta_bound[1] for :log
    less_than_zero_theta_bounds = (scale .== :log) .& [theta_bounds[i][1] < 0 for i in 1:length(theta_init)]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":log scaled theta_bound min is negative: $(findall(less_than_zero_theta_bounds))"))
    end
    # 0 <= theta_bounds <= 1 for :logit
    less_than_zero_theta_bounds = (scale .== :logit) .& [theta_bounds[i][1] < 0 || theta_bounds[i][2] > 1 for i in 1:length(theta_init)]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":logit scaled theta_bound min is outside range [0,1]: $(findall(less_than_zero_theta_bounds))"))
    end
    # loss_func(theta_init) < loss_crit
    !(loss_func(theta_init) < loss_crit) &&
        throw(ArgumentError("Check theta_init and loss_crit: loss_func(theta_init) should be < loss_crit"))

    # set counter in the scope
    counter::Int = 0
    # set supreme, maximal or minimal value of scanned parameter inside critical
    supreme_gd = nothing

    # progress info
    prog = ProgressUnknown("$direction CP counter:"; spinner=false, enabled=!silent, showspeed=true)

    # transforming loss
    theta_init_gd = scaling.(theta_init, scale)
    if isLeft theta_init_gd[theta_num] *= -1 end # change direction

    function loss_func_gd(theta_gd)
        theta_g = copy(theta_gd)
        if isLeft theta_g[theta_num] *= -1 end # change direction
        theta = unscaling.(theta_g, scale)
        # calculate function
        loss_norm = loss_func(theta) - loss_crit

        # update counter
        counter += 1
        # update supreme ?
        update_supreme = (loss_norm < 0.) &&
            (typeof(supreme_gd)==Nothing || (theta_gd[theta_num] > supreme_gd))
        if update_supreme
            supreme_gd = theta_gd[theta_num]
        end
        # display current
        supreme = if (isLeft && typeof(supreme_gd)!==Nothing)
            unscaling(-supreme_gd, scale[theta_num])
        else
            unscaling(supreme_gd, scale[theta_num])
        end
        ProgressMeter.update!(prog, counter, spinner="⠋⠙⠹⠸⠼⠴⠦⠧⠇⠏"; showvalues = [(:supreme,round(supreme; sigdigits=4))])

        return loss_norm
    end
    theta_bounds_gd = scaling.(theta_bounds, scale)
    if isLeft theta_bounds_gd[theta_num] = (-1*theta_bounds_gd[theta_num][2], -1*theta_bounds_gd[theta_num][1]) end # change direction

    scan_bound_gd = scaling(scan_bound, scale[theta_num])
    if isLeft scan_bound_gd *= -1 end # change direction

    # transform gradient
    loss_grad_gd = loss_grad
    loss_grad_gd = if isa(loss_grad, Function)
        function(theta_gd)
            theta_g = copy(theta_gd)
            if isLeft theta_g[theta_num] *= -1 end # change direction
            theta = unscaling.(theta_g, scale)
            loss_grad_value = loss_grad(theta)
            loss_grad_value_gd = zeros(n_theta)

            for i in 1:n_theta
                if scale[i] == :log
                    loss_grad_value_gd[i] = loss_grad_value[i] * theta[i] * log(10.)
                elseif scale[i] == :logit
                    loss_grad_value_gd[i] = loss_grad_value[i] * theta[i] * (1. - theta[i]) * log(10.)
                else # for :direct and :lin
                    loss_grad_value_gd[i] = loss_grad_value[i]
                end
            end

            if isLeft loss_grad_value_gd[theta_num] *= -1 end # change direction

            return loss_grad_value_gd
        end
    else
        loss_grad
    end

    # calculate endpoint using base method
    (optf_gd, pp_gd, status) = get_right_endpoint(
        theta_init_gd,
        theta_num,
        loss_func_gd,
        Val(method);
        theta_bounds = theta_bounds_gd,
        scan_bound = scan_bound_gd,
        scan_tol,
        loss_tol,
        local_alg,
        max_iter,
        loss_grad = loss_grad_gd
        #kwargs...
    )

    # transforming back
    if (isLeft && typeof(optf_gd)!==Nothing) optf_gd *= -1 end # change direction
    optf = unscaling(optf_gd, scale[theta_num])
    temp_fun = (pp::ProfilePoint) -> begin
        if isLeft pp.params[theta_num] *= -1 end # change direction
        ProfilePoint(
            unscaling(pp.params[theta_num], scale[theta_num]),
            pp.loss + loss_crit,
            unscaling.(pp.params, scale),
            pp.ret,
            pp.counter
        )
    end
    pps = [ temp_fun(pp_gd[i]) for i in 1:length(pp_gd) ]
    # transforming supreme back
    supreme = if (isLeft && typeof(supreme_gd)!==Nothing)
        unscaling(-supreme_gd, scale[theta_num])
    else
        unscaling(supreme_gd, scale[theta_num])
    end
    ProgressMeter.finish!(prog)

    EndPoint(optf, pps, status, direction, counter, supreme)
end

"""
    function get_endpoint(
        theta_init::Vector{Float64},
        scan_func::Function,
        loss_func::Function,
        method::Symbol,
        direction::Symbol = :right;

        loss_crit::Float64 = 0.0,
        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_bound::Float64 = unscaling(
            (direction==:left) ? -9.0 : 9.0,
            scale[theta_num]
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 0.,
        local_alg::Symbol = :LN_NELDERMEAD,
        kwargs...
        )

Calculates confidence interval's right or left endpoints for a function of parameters `scan_func`.

## Return
[`EndPoint`](@ref) object storing confidence interval's endpoint and intermediate profile points.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values should not necessary be the optimum values of `loss_func` but `loss_func(theta_init)` should be lower than `loss_crit`.
- `scan_func`: function of parameters.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` for profile likelihood-based (PL) identification. Usually we use log-likelihood for PL analysis: ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to estimate confidence interval's endpoint. Currently the only supported method is: `:CICO_ONE_PASS`.
- `direction`: `:right` or `:left` endpoint to estimate.

## Keyword arguments
- `loss_crit`: critical level of loss function. Confidence interval's endpoint value is the intersection point of profile likelihood and `loss_crit` level.
- `scale`: vector of scale transformations for each parameters' component. Possible values: `:direct` (`:lin`), `:log`, `:logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `theta_bounds`: vector of tuple `(lower_bound, upper_bound)` for each parameter. Bounds define the ranges for possible parameter values. Default bounds are `(-Inf,Inf)`.
- `scan_bound`: value which states the area of confidence point analysis.
- `scan_tol`: Absolute tolerance for `theta_num` parameter used as termination criterion.
- `loss_tol`: Absolute tolerance controlling `loss_func` closenes to `loss_crit` (termination criterion). Currently doesn't work for `:CICO_ONE_PASS` method because of limitation in `LN_AUGLAG` interface.
- `local_alg`: algorithm of optimization. Derivative-free and gradient-based algorithms form NLopt package.
- `max_iter` : maximal number of fitter iterations. If reaches the result status will be `:MAX_ITER_STOP`.
- `scan_grad` : For gradient optimization methods it is necessary to set how the gradient of `scan_func` should be calculated.
    - `:EMPTY` (default) no gradient is set. It works only for gradient-free methods.
    - `:AUTODIFF` means autodifferentiation from `ForwardDiff` package is used.
    - `:FINITE` means finite difference method from `Calculus` is used.
    - `function(x::Vector{Float64})` which returns gradient vector.
- `loss_grad` : For gradient optimization methods it is necessary to set how the gradient of `loss_func` should be calculated.
    - `:EMPTY` (default) no gradient is set. It works only for gradient-free methods.
    - `:AUTODIFF` means autodifferentiation from `ForwardDiff` package is used.
    - `:FINITE` means finite difference method from `Calculus` is used.
    - `function(x::Vector{Float64})` which returns gradient vector.
- `silent` : Boolean argument declaring whether we display the optimization progress. Default is `false`

See also [`get_interval`](@ref)
"""
function get_endpoint(
    theta_init::Vector{Float64},
    scan_func::Function,
    loss_func::Function,
    method::Symbol,
    direction::Symbol = :right;

    loss_crit::Float64 = 0.0,
    # :direct, :lin, :log, :logit
    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
        fill((-Inf, Inf), length(theta_init)),
        scale
        ),
    scan_bound::Float64 = (direction==:left) ? -1e9 : 1e9, # log scan bound is not implemented
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 0.,
    local_alg::Symbol = :LN_NELDERMEAD,
    max_iter::Int = 10^5,
    scan_grad::Union{Function, Symbol} = :EMPTY,
    loss_grad::Union{Function, Symbol} = :EMPTY,
    silent::Bool = false,
    #kwargs... # other options for get_right_endpoint
)
    isLeft = direction == :left
    n_theta = length(theta_init)
    loss_init = loss_func(theta_init)

    # checking arguments
    # theta_bound[1] < theta_init < theta_bound[2]
    theta_init_outside_theta_bounds = .! [theta_bounds[i][1] < theta_init[i] < theta_bounds[i][2] for i in 1:length(theta_init)]
    if any(theta_init_outside_theta_bounds)
        throw(ArgumentError("theta_init is outside theta_bound: $(findall(theta_init_outside_theta_bounds))"))
    end

    # 0 <= theta_bound[1] for :log
    less_than_zero_theta_bounds = (scale .== :log) .& [theta_bounds[i][1] < 0 for i in 1:length(theta_init)]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":log scaled theta_bound min is negative: $(findall(less_than_zero_theta_bounds))"))
    end
    # 0 <= theta_bounds <= 1 for :logit
    less_than_zero_theta_bounds = (scale .== :logit) .& [theta_bounds[i][1] < 0 || theta_bounds[i][2] > 1 for i in 1:length(theta_init)]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":logit scaled theta_bound min is outside range [0,1]: $(findall(less_than_zero_theta_bounds))"))
    end
    # loss_func(theta_init) < loss_crit
    !(loss_init < loss_crit) &&
        throw(ArgumentError("Check theta_init and loss_crit: loss_func(theta_init) should be < loss_crit"))

    # set counter in the scope
    prog = ProgressUnknown("Fitter counter:"; spinner=false, enabled=!silent, showspeed=true)
    counter::Int = 0
    # set supreme, maximal or minimal value of loss_func inside critical
    supreme_gd = nothing
    scan_val_gd = nothing

    # transforming
    theta_init_gd = scaling.(theta_init, scale)

    function scan_func_gd(theta_gd)
        theta = unscaling.(theta_gd, scale)
        scan_val = scan_func(theta)
        scan_val_gd = isLeft ? (-1)*scan_val : scan_val

        return scan_val_gd
    end

    function loss_func_gd(theta_gd)
        #theta_g = copy(theta_gd) # why copy?
        theta = unscaling.(theta_gd, scale)
        loss_value = loss_func(theta) - loss_crit

        counter += 1
        #println("$loss_value => $supreme_gd")
        should_update = (loss_value < 0.) &&
            !isa(scan_val_gd, Nothing) &&
            (isa(supreme_gd, Nothing) || (scan_val_gd > supreme_gd)) &&
            !isa(scan_val_gd, ForwardDiff.Dual)
        if should_update
            supreme_gd = scan_val_gd
        end
        supreme = if isa(supreme_gd, Nothing)
            "-"
        elseif isLeft
            round(-supreme_gd; sigdigits=4)
        else
            round(supreme_gd; sigdigits=4)
        end
        ProgressMeter.update!(prog, counter; showvalues = [(:supreme,supreme)])

        return loss_value
    end
    theta_bounds_gd = scaling.(theta_bounds, scale)

    # TODO: transformed by scan_scale: scaling(scan_bound, scan_scale)
    scan_bound_gd = scan_bound
    if isLeft scan_bound_gd *= -1 end # change direction

    # transform gradients
    scan_grad_gd = if isa(scan_grad, Function)
        function(theta_gd)
            theta = unscaling.(theta_gd, scale)
            scan_grad_value = scan_grad(theta)
            scan_grad_value_gd = zeros(n_theta)

            for i in 1:n_theta
                if scale[i] == :log
                    scan_grad_value_gd[i] = scan_grad_value[i] * theta[i] * log(10.)
                elseif scale[i] == :logit
                    scan_grad_value_gd[i] = scan_grad_value[i] * theta[i] * (1. - theta[i]) * log(10.)
                else # for :direct and :lin
                    scan_grad_value_gd[i] = scan_grad_value[i]
                end
            end

            return scan_grad_value_gd
        end
    else
        scan_grad
    end

    loss_grad_gd = if isa(loss_grad, Function)
        function(theta_gd)
            theta = unscaling.(theta_gd, scale)
            loss_grad_value = loss_grad(theta)
            loss_grad_value_gd = zeros(n_theta)

            for i in 1:n_theta
                if scale[i] == :log
                    loss_grad_value_gd[i] = loss_grad_value[i] * theta[i] * log(10.)
                elseif scale[i] == :logit
                    loss_grad_value_gd[i] = loss_grad_value[i] * theta[i] * (1. - theta[i]) * log(10.)
                else # for :direct and :lin
                    loss_grad_value_gd[i] = loss_grad_value[i]
                end
            end

            return loss_grad_value_gd
        end
    else
        loss_grad
    end

    # calculate endpoint using base method
    (optf_gd, pp_gd, status) = get_right_endpoint(
        theta_init_gd,
        scan_func_gd,
        loss_func_gd,
        Val(method);

        theta_bounds = theta_bounds_gd,
        scan_bound = scan_bound_gd,
        scan_tol,
        loss_tol,
        local_alg,
        max_iter,
        scan_grad = scan_grad_gd,
        loss_grad = loss_grad_gd,
        #kwargs...
    )

    # transforming back
    temp_fun = (pp::ProfilePoint) -> begin
        value = isLeft ? (-1)*pp.value : pp.value # TODO: scan_scale
        ProfilePoint(
            value,
            pp.loss + loss_crit,
            unscaling.(pp.params, scale),
            pp.ret,
            pp.counter
        )
    end
    pps = [ temp_fun(pp_gd[i]) for i in 1:length(pp_gd) ]
    # transforming supreme back
    if (isLeft && typeof(supreme_gd)!==Nothing) supreme_gd *= -1 end # change direction

    if (isLeft && typeof(optf_gd)!==Nothing) optf_gd *= -1 end # change direction
    # optf = unscaling(optf_gd, scan_scale)
    optf = optf_gd

    EndPoint(optf, pps, status, direction, counter, supreme_gd)
end


text:
#=
diffl.jl
timing tests
=#

using MIRT: diffl!, diffl_adj!, diffl, diffl_adj # new
using MIRT: diff_forw, diff_adj # old
using Test: @test
using BenchmarkTools: @btime

#=
# timing tests show ~40% speedup for in-place version
=#

function diffl_time()

    @info "time vs diff and vs w/o in-place"

    N = (2^10, 2^9)
    x = rand(N...)
    g0 = rand(N...)
    g1 = rand(N...)
    g0 = diff(x, dims=2)
    diffl!(g1, x, 2)
    @test g1[:,2:end] == g0
    @btime diffl!($g1, $x, 2)    # 210 20 672 b
    @btime g0 = diff($x, dims=2) # 303  4 3.99 MiB
    @btime g1 = diffl($x, 2)     # 306 22 4.00 MiB

    @info "time vs diff_forw"
    # in-place version diffl! about 40% faster than diff_forw

    N = (2^5, 2^6, 2^7)
    x = rand(N...)
    d = 2
    g2 = diffl(x, d)
    h2 = diff_forw(x, dims=d)
    @test vec(g2[:,2:end,:]) == h2
    diffl!(g2, x, d)
    @test vec(g2[:,2:end,:]) == h2
    @btime diffl!($g2, $x, $d)     # 92 24 800 b
    @btime diff_forw($x, dims=$d)  # 148 6 1.97 MiB

    @info "time adj" # about 40% speed up

    N = (2^5, 2^6, 2^7)
    x = rand(N...)
    d = 2
    g2 = diffl(x, d)
    h2 = diff_forw(x, dims=d)
    z2 = diff_adj(h2, N ; dims=d)
    z1 = diffl_adj(g2, d)
    @test z1 == z2
    diffl_adj!(z1, g2, d)
    @test z1 == z2

    @btime diffl_adj!($z1, $g2, $d)     # 361 28 928 b
    @btime diff_adj($h2, $N ; dims=$d)  # 607 70 3.94 MiB

    nothing
end

#diffl_time()

#=
Below here are timing tests to compare MIRT.diffl! with StaticKernels
The "in place" versions have nearly identical compute times.
All commented out to avoid the dependency on StaticKernels.jl

using StaticKernels

function diffl_race()

    N = (2^10, 2^9)
    x = rand(N...)
    g1 = similar(x)
    diffl!(g1, x, 2)

    @inline kf(w) = @inbounds w[0,0] - w[0,-1] # for d=2
    k = Kernel{(0:0,-1:0)}(kf)

    g2 = similar(x)
#   g2 .= map(k, extend(x, StaticKernels.ExtensionConstant(0)))
    g2 .= map(k, extend(x, StaticKernels.ExtensionReplicate()))
    @test g1 == g2

#   @btime $g2 .= map($k, extend($x, StaticKernels.ExtensionReplicate()))
#   694.090 μs (3 allocations: 4.00 MiB) [with or without]
    @btime map!($k, $g2, extend($x, StaticKernels.ExtensionReplicate()))
#   283.389 μs (1 allocation: 16 bytes) [without]
#   192.728 μs (1 allocation: 16 bytes) [with inline inbounds]
    @btime diffl!($g1, $x, 2)  # 188 μs 20 672 b

    nothing
end

diffl_race()

=#


text:
int_rules_4_1_1_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.1.1*(a+b*sin)^n =#
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x) ^ ~n, ~x) => (-1 / ~d) * Subst(Antiderivative(expand((1 - (~x) ^ 2) ^ ((~n - 1) / 2), ~x), ~x), ~x, cos(~c + ~d * ~x)) <-- FreeQ([~c, ~d], ~x) && IGtQ((~n - 1) / 2, 0)
    @apply_utils Antiderivative(sin(~(c') + (~(d') * ~x) / 2) ^ 2, ~x) => ~x / 2 - sin(2 * ~c + ~d * ~x) / (2 * ~d) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative((~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~b) * cos(~c + ~d * ~x) * (~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * ~n) + (((~b) ^ 2 * (~n - 1)) / ~n) * Antiderivative((~b * sin(~c + ~d * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && (GtQ(~n, 1) && IntegerQ(2 * ~n))
    @apply_utils Antiderivative((~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (cos(~c + ~d * ~x) * (~b * sin(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~n + 2) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~b * sin(~c + ~d * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && (LtQ(~n, -1) && IntegerQ(2 * ~n))
    @apply_utils Antiderivative(sin(~(c') + pi / 2 + ~(d') * ~x), ~x) => sin(~c + ~d * ~x) / ~d <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * ~x), ~x) => -(cos(~c + ~d * ~x)) / ~d <-- FreeQ([~c, ~d], ~x)
    #= Antiderivative(1/sin((~c')+(~d')*(~x)),~x) := Antiderivative(csc(c+d*x),x) <--  FreeQ([c,d],x) =#
    @apply_utils Antiderivative(sqrt(sin(~(c') + ~(d') * ~x)), ~x) => (2 / ~d) * Elliptic.E((1 / 2) * ((~c - pi / 2) + ~d * ~x), 2) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative(1 / sqrt(sin(~(c') + ~(d') * ~x)), ~x) => (2 / ~d) * Elliptic.F((1 / 2) * ((~c - pi / 2) + ~d * ~x), 2) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative((~b * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((~b * sin(~c + ~d * ~x)) ^ ~n / sin(~c + ~d * ~x) ^ ~n) * Antiderivative(sin(~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~b, ~c, ~d], ~x) && (LtQ(-1, ~n, 1) && IntegerQ(2 * ~n))
    #= Antiderivative(((~b')*sin((~c')+(~d')*(~x)))^(~n),~x) :=  cos(c+d*x)/(b*d*sqrt(cos(c+d*x)^2))*Subst(Antiderivative(x^n/sqrt(1-x^2/b^2),x), x,b*sin(c+d*x)) <-- FreeQ([b,c,d,n],x) && Not(IntegerQ(2*n) ||  IntegerQ(3*n)) =#
    @apply_utils Antiderivative((~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((cos(~c + ~d * ~x) * (~b * sin(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1) * sqrt(cos(~c + ~d * ~x) ^ 2))) * HypergeometricFunctions._₂F₁(1 / 2, (~n + 1) / 2, (~n + 3) / 2, sin(~c + ~d * ~x) ^ 2) <-- FreeQ([~b, ~c, ~d, ~n], ~x) && Not(IntegerQ(2 * ~n))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ 2, ~x) => (((2 * (~a) ^ 2 + (~b) ^ 2) * ~x) / 2 - (2 * ~a * ~b * cos(~c + ~d * ~x)) / ~d) - ((~b) ^ 2 * cos(~c + ~d * ~x) * sin(~c + ~d * ~x)) / (2 * ~d) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandTrig((~a + ~b * sin(~c + ~d * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~n, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (-2 * ~b * cos(~c + ~d * ~x)) / (~d * sqrt(~a + ~b * sin(~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~b) * cos(~c + ~d * ~x) * (~a + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * ~n) + ((~a * (2 * ~n - 1)) / ~n) * Antiderivative((~a + ~b * sin(~c + ~d * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~n - 1 / 2, 0))
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => -(cos(~c + ~d * ~x)) / (~d * (~b + ~a * sin(~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (-2 / ~d) * Subst(Antiderivative(1 / (2 * ~a - (~x) ^ 2), ~x), ~x, (~b * cos(~c + ~d * ~x)) / sqrt(~a + ~b * sin(~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * cos(~c + ~d * ~x) * (~a + ~b * sin(~c + ~d * ~x)) ^ ~n) / (~a * ~d * (2 * ~n + 1)) + ((~n + 1) / (~a * (2 * ~n + 1))) * Antiderivative((~a + ~b * sin(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~n, -1) && IntegerQ(2 * ~n)))
    #= Antiderivative(((~a)+(~b')*sin((~c')+(~d')*(~x)))^(~n),~x) :=  a^2*cos(c+d*x)/(d*sqrt(a+b*sin(c+d*x))*sqrt(a-b*sin(c+d*x)))*Subst( Antiderivative((a+b*x)^(n-1/2)/sqrt(a-b*x),x),x,sin(c+d*x)) <--  FreeQ([a,b,c,d,n],x) && EqQ(a^2-b^2,0) && Not(IntegerQ(2*n)) =#
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => ((-(2 ^ (~n + 1 / 2)) * (~a) ^ (~n - 1 / 2) * ~b * cos(~c + ~d * ~x)) / (~d * sqrt(~a + ~b * sin(~c + ~d * ~x)))) * HypergeometricFunctions._₂F₁(1 / 2, 1 / 2 - ~n, 3 / 2, (1 / 2) * (1 - (~b * sin(~c + ~d * ~x)) / ~a)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(IntegerQ(2 * ~n)) && GtQ(~a, 0)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (((~a) ^ IntPart(~n) * (~a + ~b * sin(~c + ~d * ~x)) ^ FracPart(~n)) / (1 + (~b / ~a) * sin(~c + ~d * ~x)) ^ FracPart(~n)) * Antiderivative((1 + (~b / ~a) * sin(~c + ~d * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(IntegerQ(2 * ~n)) && Not(GtQ(~a, 0))))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => ((2 * sqrt(~a + ~b)) / ~d) * Elliptic.E((1 / 2) * ((~c - pi / 2) + ~d * ~x), (2 * ~b) / (~a + ~b)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && GtQ(~a + ~b, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => ((2 * sqrt(~a - ~b)) / ~d) * Elliptic.E((1 / 2) * (~c + pi / 2 + ~d * ~x), (-2 * ~b) / (~a - ~b)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && GtQ(~a - ~b, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (sqrt(~a + ~b * sin(~c + ~d * ~x)) / sqrt((~a + ~b * sin(~c + ~d * ~x)) / (~a + ~b))) * Antiderivative(sqrt(~a / (~a + ~b) + (~b / (~a + ~b)) * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(GtQ(~a + ~b, 0)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~b) * cos(~c + ~d * ~x) * (~a + ~b * sin(~c + ~d * ~x)) ^ (~n - 1)) / (~d * ~n) + (1 / ~n) * Antiderivative((~a + ~b * sin(~c + ~d * ~x)) ^ (~n - 2) * Simp((~a) ^ 2 * ~n + (~b) ^ 2 * (~n - 1) + ~a * ~b * (2 * ~n - 1) * sin(~c + ~d * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~n, 1) && IntegerQ(2 * ~n)))
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => With([q = Rt((~a) ^ 2 - (~b) ^ 2, 2)], ~x / q + (2 / (~d * q)) * atan((~b * cos(~c + ~d * ~x)) / (~a + q + ~b * sin(~c + ~d * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (GtQ((~a) ^ 2 - (~b) ^ 2, 0) && PosQ(~a))
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => With([q = Rt((~a) ^ 2 - (~b) ^ 2, 2)], -(~x) / q - (2 / (~d * q)) * atan((~b * cos(~c + ~d * ~x)) / ((~a - q) + ~b * sin(~c + ~d * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (GtQ((~a) ^ 2 - (~b) ^ 2, 0) && NegQ(~a))
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(c') + pi / 2 + ~(d') * ~x)), ~x) => With([e = FreeFactors(tan((~c + ~d * ~x) / 2), ~x)], ((2 * e) / ~d) * Subst(Antiderivative(1 / (~a + ~b + (~a - ~b) * e ^ 2 * (~x) ^ 2), ~x), ~x, tan((~c + ~d * ~x) / 2) / e)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / (~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => With([e = FreeFactors(tan((~c + ~d * ~x) / 2), ~x)], ((2 * e) / ~d) * Subst(Antiderivative(1 / (~a + 2 * ~b * e * ~x + ~a * e ^ 2 * (~x) ^ 2), ~x), ~x, tan((~c + ~d * ~x) / 2) / e)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (2 / (~d * sqrt(~a + ~b))) * Elliptic.F((1 / 2) * ((~c - pi / 2) + ~d * ~x), (2 * ~b) / (~a + ~b)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && GtQ(~a + ~b, 0))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (2 / (~d * sqrt(~a - ~b))) * Elliptic.F((1 / 2) * (~c + pi / 2 + ~d * ~x), (-2 * ~b) / (~a - ~b)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && GtQ(~a - ~b, 0))
    @apply_utils Antiderivative(1 / sqrt(~a + ~(b') * sin(~(c') + ~(d') * ~x)), ~x) => (sqrt((~a + ~b * sin(~c + ~d * ~x)) / (~a + ~b)) / sqrt(~a + ~b * sin(~c + ~d * ~x))) * Antiderivative(1 / sqrt(~a / (~a + ~b) + (~b / (~a + ~b)) * sin(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(GtQ(~a + ~b, 0)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~b) * cos(~c + ~d * ~x) * (~a + ~b * sin(~c + ~d * ~x)) ^ (~n + 1)) / (~d * (~n + 1) * ((~a) ^ 2 - (~b) ^ 2)) + (1 / ((~n + 1) * ((~a) ^ 2 - (~b) ^ 2))) * Antiderivative((~a + ~b * sin(~c + ~d * ~x)) ^ (~n + 1) * Simp(~a * (~n + 1) - ~b * (~n + 2) * sin(~c + ~d * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~n, -1) && IntegerQ(2 * ~n)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (cos(~c + ~d * ~x) / (~d * sqrt(1 + sin(~c + ~d * ~x)) * sqrt(1 - sin(~c + ~d * ~x)))) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n / (sqrt(1 + ~x) * sqrt(1 - ~x)), ~x), ~x, sin(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(IntegerQ(2 * ~n)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(c') + ~(d') * ~x) * cos(~(c') + ~(d') * ~x)) ^ ~n, ~x) => Antiderivative((~a + (~b * sin(2 * ~c + 2 * ~d * ~x)) / 2) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.1.1 (a+b sin)^n.jl


text:
"""
    axes3D(n=100)

Draw labelled 3D axes at `(0, 0, 0)` with length `n`.
"""
function axes3D(n=100)
    projection = CURRENTPROJECTION[1]
    @layer begin
        fontsize(16)
        setline(2)
        xaxis1 = project(Point3D(0.1,  0.1,  0.1))
        xaxis2 = project(Point3D(n,    0.1,  0.1))
        sethue("red")
        if (xaxis1 != nothing) && (xaxis2 != nothing) && !isapprox(xaxis1, xaxis2)
            arrow(xaxis1, xaxis2)
            label("X", :N, xaxis2)
        end
        yaxis1 = project(Point3D(0.1,  0.1,  0.1))
        yaxis2 = project(Point3D(0.1,  n,  0.1))
        sethue("green")

        if (yaxis1 != nothing) &&  (yaxis2 != nothing) && !isapprox(yaxis1, yaxis2)
            arrow(yaxis1, yaxis2)
            label("Y", :N, yaxis2)
        end
        zaxis1 = project(Point3D(0.1,  0.1,  0.1))
        zaxis2 = project(Point3D(0.1,  0.1,  n))
        sethue("blue")
        if (zaxis1 != nothing) &&  (zaxis2 != nothing) && !isapprox(zaxis1, zaxis2)
            arrow(zaxis1, zaxis2)
            label("Z", :N, zaxis2)
        end
    end
end

"""
    drawcube(n=10, action=:stroke)

Draw a cube. `drawcube(1)` draws a wireframe unit cube.
"""
function drawcube(n=10, action=:stroke)
    @layer begin
        setline(2)
        fontsize(10)
        setlinecap("butt")
        setlinejoin("bevel")
        p1 = project(Point3D(0,   0,  0))
        p2 = project(Point3D(n,  0,  0))
        p3 = project(Point3D(n, n,  0))
        p4 = project(Point3D(0,  n,  0))
        p5 = project(Point3D(0,   0,  n))
        p6 = project(Point3D(n,  0,  n))
        p7 = project(Point3D(n,  n,  n))
        p8 = project(Point3D(0,  n,  n))

        if all(i -> i != nothing, [p1, p2, p3, p4]) # , p5, p6, p7, p8])
            sethue("red")
            prettypoly([p1, p2, p3, p4], action, close=true, () -> begin
                circle(O, n/100)
            end)
        end

        if all(i -> i != nothing, [p1, p2, p6, p5])
            sethue("purple")
            prettypoly([p1, p2, p6, p5], action, close=true, () -> begin
                circle(O, n/100)
            end)
        end

        if all(i -> i != nothing, [p3, p7, p8, p4])
            sethue("magenta")
            prettypoly([p3, p7, p8, p4], action, close=true, () -> begin
                circle(O, n/100)
            end)
        end

        if all(i -> i != nothing, [p5, p6, p7, p8])
            sethue("green")
            prettypoly([p5, p6, p7, p8], action, close=true, () -> begin
                circle(O, n/100)
            end)
        end
    end
end

"""
    carpet(n; kind=:circular)

Draw a circular carpet centered at the origin, using current Luxor parameters.

If `kind` is not `:circular`, the carpet will be a square.

Points that can't be rendered are not included in the final shape.
"""
function carpet(n=100; kind=:circular)
    projection = CURRENTPROJECTION[1]
    @layer begin
        if kind != :circular
            p1 = project(Point3D(-n/2, -n/2,  0))
            p2 = project(Point3D(n/2,  -n/2,  0))
            p3 = project(Point3D(n/2, n/2,  0))
            p4 = project(Point3D(-n/2,  n/2,  0))
            if all(i -> i != nothing, [p1, p2, p3,  p4])
                poly([p1, p2, p3, p4], :fill, close=true)
            end
        else
            shape3D = [Point3D(n * cos(theta), n * sin(theta), 0) for theta in 0:0.1:2pi]
            shape2D = Point[]
            for i in shape3D
                pt = project(i)
                if pt != nothing
                    push!(shape2D, pt)
                end
            end
            poly(shape2D, :fill, close=true)
        end
    end
end

"""
    crossproduct3D(A::Point3D, B::Point3D)

Find one of these.
"""
function crossproduct3D(A::Point3D, B::Point3D)
    x = (A.y * B.z) - (A.z * B.y)
    y = (A.z * B.x) - (A.x * B.z)
    z = (A.x * B.y) - (A.y * B.x)
    return Point3D(x, y, z)
end

"""
    pointsperpendicular(p1::Point3D, p2::Point3D, radius, angles = [0, π])

Find points perpendicular to a line joining p1 and p2.
Points are `radius` units away from the line.
"""
function pointsperpendicular(p1::Point3D, p2::Point3D, radius, angles = [0, π])
    V3 = p2 - p1 # unit vector
    V3 /= distance(Point3D(0.0, 0.0, 0.0), V3)
    e = Point3D(0.0, 0.0, 1.0)
    V1 = crossproduct3D(e, V3) # unit vector perpendicular to both e and V3
    V1 /= distance(Point3D(0.0, 0.0, 0.0), V1)
    V2 = crossproduct3D(V3, V1) # third unit vector perpendicular to both V3 and V1
    # find points P on the circle θ radians relative to vector V1 (x-axis)
    return [p1 + radius * (cos(θ) * V1 + sin(θ) * V2) for θ in angles]
end


text:
## Reference API ##
tickers_base_url = "https://api.polygon.io/v3/reference/tickers"

ticker_types_base_url = "https://api.polygon.io/v2/reference/types"

ticker_news_base_url = "https://api.polygon.io/v2/reference/news"

markets_base_url = "https://api.polygon.io/v2/reference/markets"

locales_base_url = "https://api.polygon.io/v2/reference/locales"

market_holidays_base_url = "https://api.polygon.io/v1/marketstatus/upcoming"

market_status_base_url = "https://api.polygon.io/v1/marketstatus/now"

stock_exchanges_base_url = "https://api.polygon.io/v1/meta/exchanges"

condition_mappings_base_url = "https://api.polygon.io/v1/meta/conditions"

crypto_exchanges_base_url = "https://api.polygon.io/v1/meta/crypto-exchanges"

ticker_details_base_url = "https://api.polygon.io/v1/meta/symbols"

ticker_details_vX_base_url = "https://api.polygon.io/vX/reference/tickers"

stock_splits_base_url = "https://api.polygon.io/v2/reference/splits"

stock_dividends_base_url = "https://api.polygon.io/v2/reference/dividends"

stock_financials_base_url = "https://api.polygon.io/v2/reference/financials"

## Stock API ##
stocks_trades_base_url = "https://api.polygon.io/v2/ticks/stocks/trades"

stocks_quotes_base_url = "https://api.polygon.io/v2/ticks/stocks/nbbo"

stocks_last_trade_base_url = "https://api.polygon.io/v2/last/trade"

stocks_last_quote_base_url = "https://api.polygon.io/v2/last/nbbo"

stocks_daily_open_close_base_url = "https://api.polygon.io/v1/open-close"

stocks_grouped_daily_bars_base_url = "https://api.polygon.io/v2/aggs/grouped/locale/us/market/stocks"

stocks_previous_close_base_url = "https://api.polygon.io/v2/aggs/ticker"

stocks_aggregates_bars_base_url = "https://api.polygon.io/v2/aggs/ticker"

stocks_snapshot_all_tickers_base_url = "https://api.polygon.io/v2/snapshot/locale/us/markets/stocks/tickers"

stocks_snapshot_ticker_base_url = "https://api.polygon.io/v2/snapshot/locale/us/markets/stocks/tickers"

stocks_snapshot_gainers_losers_base_url = "https://api.polygon.io/v2/snapshot/locale/us/markets/stocks"

## Crypto API ##
crypto_last_trade_pair_base_url = "https://api.polygon.io/v1/last/crypto"

crypto_daily_open_close_base_url = "https://api.polygon.io/v1/open-close/crypto"

crypto_historic_trades_base_url = "https://api.polygon.io/v1/historic/crypto"

crypto_grouped_daily_bars_base_url = "https://api.polygon.io/v2/aggs/grouped/locale/global/market/crypto"

crypto_previous_close_base_url = "https://api.polygon.io/v2/aggs/ticker"

crypto_aggregates_bars_base_url = "https://api.polygon.io/v2/aggs/ticker"

crypto_snapshot_all_tickers_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/crypto/tickers"

crypto_snapshot_ticker_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/crypto/tickers"

crypto_snapshot_ticker_full_book_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/crypto/tickers"

crypto_snapshot_gainers_losers_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/crypto"

## Forex API ##
historic_forex_ticks_base_url = "https://api.polygon.io/v1/historic/forex"

real_time_currency_conversion_base_url = "https://api.polygon.io/v1/conversion"

last_quote_currency_pair_base_url = "https://api.polygon.io/v1/last_quote/currencies"

forex_grouped_daily_bars_base_url = "https://api.polygon.io/v2/aggs/grouped/locale/global/market/fx"

forex_previous_close_base_url = "https://api.polygon.io/v2/aggs/ticker"

forex_aggregates_bars_base_url = "https://api.polygon.io/v2/aggs/ticker"

forex_snapshot_ticker_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/forex/tickers"

forex_snapshot_all_tickers_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/forex/tickers"

forex_snapshot_gainers_losers_base_url = "https://api.polygon.io/v2/snapshot/locale/global/markets/forex"


text:
export FieldStyle, HField, EField

abstract type FieldStyle end
struct HField <: FieldStyle end
struct EField <: FieldStyle end

"""
    _gstcMatrix(::FieldStyle, A, B, C, D, )

Converts gstcs `A(H² - H¹) = B(E² + E¹)` and `C(E² - E¹) = D(H² + H¹)` into a
corresponding scattering matrix for either the E or H field
"""
function _gstcMatrix(::HField, A, B₁, B₂, C₁, C₂, D)
    out = [
        -A + B₁    A - B₂
        C₁ - D     C₂ - D
    ]
    inc = [
        A + B₁     -A - B₂
        C₁ + D     C₂ + D
    ]
    out, inc
end
_gstcMatrix(fs::HField, A, B, C, D) = _gstcMatrix(fs, A, B, B, C, C, D)

function _gstcMatrix(::EField, A₁, A₂, B, C, D₁, D₂)
    out, inc = _gstcMatrix(HField(), C, D₁, D₂, A₁, A₂, B)
    n = Int(size(out, 1)//2)
    prop_sign = Diagonal(vcat(ones(n), -ones(n)))
    out * -prop_sign, inc * prop_sign
end
_gstcMatrix(fs::EField, A, B, C, D) = _gstcMatrix(fs, A, A, B, C, D, D)

function _get_halves(x::AbstractVector)
    n = Int(length(x)//2)
    I₁ = x[1:n]
    I₂ = x[(n+1):2n]
    I₁, I₂
end

function _gstcLinearMap(::HField, A, B₁, B₂, C₁, C₂, D)
    n = 2size(A, 1)
    out = LinearMap(n) do x
        I₁, I₂ = _get_halves(x)
        vcat(
            A*(-I₁ + I₂) + B₁*I₁ - B₂*I₂,
            C₁*I₁ + C₂*I₂ - D*(I₁ + I₂),
        )
    end
    inc = LinearMap(n) do x
        I₁, I₂ = _get_halves(x)
        vcat(
            A*(I₁ - I₂) + B₁*I₁ - B₂*I₂,
            C₁*I₁ + C₂*I₂ + D*(I₁ + I₂),
        )
    end
    out, inc
end

function _gstcLinearMap(::HField, A, B, C, D)
    n = 2size(A, 1)
    out = LinearMap(n) do x
        I₁, I₂ = _get_halves(x)
        sumI₁I₂ = I₁ + I₂
        diffI₁I₂ = I₁ - I₂
        vcat(
            -A*diffI₁I₂ + B*diffI₁I₂,
            C*sumI₁I₂ - D*sumI₁I₂,
        )
    end
    inc = LinearMap(n) do x
        I₁, I₂ = _get_halves(x)
        sumI₁I₂ = I₁ + I₂
        diffI₁I₂ = I₁ - I₂
        vcat(
            A*diffI₁I₂ + B*diffI₁I₂,
            C*sumI₁I₂ + D*sumI₁I₂,
        )
    end
    out, inc
end

function _gstcLinearMap(::EField, A₁, A₂, B, C, D₁, D₂)
    out, inc = _gstcLinearMap(HField(), C, D₁, D₂, A₁, A₂, B)
    n = Int(size(out, 1)//2)
    prop_sign = Diagonal(vcat(ones(n), -ones(n)))
    out * -prop_sign, inc * prop_sign
end

function _gstcLinearMap(::EField, A, B, C, D)
    out, inc = _gstcLinearMap(HField(), C, D, A, B)
    n = Int(size(out, 1)//2)
    prop_sign = Diagonal(vcat(ones(n), -ones(n)))
    out * -prop_sign, inc * prop_sign
end


text:
using MLJTime
using Test

@testset "interval based forest" begin
    X, y = ts_dataset("Chinatown")
    rng = StableRNG(566) # seed to resproduce the results.

    Interval_features, Intervals = MLJTime.IntervalBasedForest.InvFeatureGen(matrix(X[1:5]), 1, 3, rng);

    @test Interval_features  == [[745.4444444444445 505.0103344571245 56.36945304437562 716.3157894736842 506.93973481750487 39.17543859649122 1057.6 311.9797714240104 -88.88484848484839 1197.5 66.31490531295862 13.4;
    643.6111111111111 413.3816266853738 45.493292053663545 619.578947368421 415.16761217468587 31.45964912280703 904.5 265.43057263414266 -75.93333333333331 1046.5 36.27211968808366 9.4;
    944.1666666666666 620.4238826138522 100.38906088751294 939.1052631578947 603.3471356374133 83.81228070175439 1421.2 164.37342310185736 -7.357575757575845 1406.75 169.44689433565904 119.1;
    1023.6666666666666 700.7234916603737 62.30340557275547 982.8421052631579 703.8473842750836 40.710526315789444 1366.3 455.57339937953157 -146.3454545454545 1541.5 139.28507936363224 -102.0;
    659.2222222222222 434.21268244349017 67.04643962848297 660.3157894736842 422.00579941138494 57.31754385964911 982.8 108.05327903914397 -15.757575757575758 922.25 88.08092869628476 7.5]]

    train, test = partition(eachindex(y), 0.7)
    @test unique(y) == [1.0 , 2.0]

    model = TimeSeriesForestClassifier(n_trees=3, random_state=rng)
    mach = machine(model, X[train], y[train])
    fit!(mach)
    y_pred = predict_mode(mach, X[test])
    @test accuracy(y_pred, y[test]) >= 0.75
end

@testset "distances" begin
    rng = StableRNG(566)
    a = rand(rng, Int64, 5)
    b = rand(rng, Int64, 5)
    @test dtw_distance(a, b, -1) == -2.9402683278543684e19
end

@testset "KNN" begin
    X, y = ts_dataset("Chinatown")
    train, test = partition(eachindex(y), 0.7)

    A = rand(StableRNG(566), 5, 5)
    index = MLJTime.select_sort(A, 3)
    @test index == [1.0 2.0 3.0; 1.0 3.0 4.0; 3.0 2.0 5.0; 1.0 5.0 4.0; 5.0 2.0 4.0]

    model = TimeSeriesKNNClassifier()
    mach = machine(model, X[train], y[train])
    fit!(mach)
    y_pred = predict_mode(mach, X[test])
    @test accuracy(y_pred, y[test]) >= 0.70
end