text:
@testset "Format Options" begin
    @testset "whitespace in typedefs" begin
        str_ = "Foo{A,B,C}"
        str = "Foo{A, B, C}"
        @test fmt(str_, whitespace_typedefs = true) == str

        str_ = """
        struct Foo{A<:Bar,Union{B<:Fizz,C<:Buzz},<:Any}
            a::A
        end"""
        str = """
        struct Foo{A <: Bar, Union{B <: Fizz, C <: Buzz}, <:Any}
            a::A
        end"""
        @test fmt(str_, whitespace_typedefs = true) == str

        str_ = """
        function foo() where {A,B,C{D,E,F{G,H,I},J,K},L,M<:N,Y>:Z}
            body
        end
        """
        str = """
        function foo() where {A, B, C{D, E, F{G, H, I}, J, K}, L, M <: N, Y >: Z}
            body
        end
        """
        @test fmt(str_, whitespace_typedefs = true) == str

        str_ = "foo() where {A,B,C{D,E,F{G,H,I},J,K},L,M<:N,Y>:Z} = body"
        str = "foo() where {A, B, C{D, E, F{G, H, I}, J, K}, L, M <: N, Y >: Z} = body"
        @test fmt(str_, whitespace_typedefs = true) == str
    end

    @testset "whitespace ops in indices" begin
        str = "arr[1 + 2]"
        @test fmt("arr[1+2]", m = 1, whitespace_ops_in_indices = true) == str

        str = "arr[(1 + 2)]"
        @test fmt("arr[(1+2)]", m = 1, whitespace_ops_in_indices = true) == str

        str_ = "arr[1:2*num_source*num_dump-1]"
        str = "arr[1:(2 * num_source * num_dump - 1)]"
        @test fmt(str_, m = 1, whitespace_ops_in_indices = true) == str

        str_ = "arr[2*num_source*num_dump-1:1]"
        str = "arr[(2 * num_source * num_dump - 1):1]"
        @test fmt(str_, m = 1, whitespace_ops_in_indices = true) == str

        str = "arr[(a + b):c]"
        @test fmt("arr[(a+b):c]", m = 1, whitespace_ops_in_indices = true) == str

        str = "arr[a in b]"
        @test fmt(str, m = 1, whitespace_ops_in_indices = true) == str

        str_ = "a:b+c:d-e"
        str = "a:(b + c):(d - e)"
        @test fmt(str_, m = 1, whitespace_ops_in_indices = true) == str

        # issue 180
        str_ = "s[m+i+1]"
        str = "s[m+i+1]"
        @test fmt(str, m = 1) == str

        str = "s[m + i + 1]"
        @test fmt(str_, m = 1, whitespace_ops_in_indices = true) == str
    end

    @testset "rewrite import to using" begin
        str_ = "import A"
        str = "using A: A"
        @test fmt(str_, import_to_using = true) == str

        str_ = """
        import A,

        B, C"""
        str = """
        using A: A
        using B: B
        using C: C"""
        @test_broken fmt(str_, import_to_using = true) == str

        str_ = """
        import A,
               # comment
        B, C"""
        str = """
        using A: A
        # comment
        using B: B
        using C: C"""
        @test fmt(str_, import_to_using = true) == str

        str_ = """
        import A, # inline
               # comment
        B, C # inline"""
        str = """
        using A: A # inline
        # comment
        using B: B
        using C: C # inline"""
        @test fmt(str_, import_to_using = true) == str

        str_ = """
        import ..A, .B, ...C"""
        str = """
        using ..A: A
        using .B: B
        using ...C: C"""
        @test fmt(str_, import_to_using = true) == str
        t = run_pretty(str_, opts = Options(max_margin = 80, import_to_using = true))
        @test t.len == 13

        # #232
        str = """import A.b"""
        @test fmt(str, import_to_using = true) == str
    end

    @testset "always convert `=` to `in` (for loops)" begin
        str_ = """
        for i = 1:n
            println(i)
        end"""
        str = """
        for i in 1:n
            println(i)
        end"""
        @test fmt(str_, always_for_in = true) == str
        @test fmt(str, always_for_in = true) == str

        str_ = """
        for i = I1, j in I2
            println(i, j)
        end"""
        str = """
        for i in I1, j in I2
            println(i, j)
        end"""
        @test fmt(str_, always_for_in = true) == str
        @test fmt(str, always_for_in = true) == str

        str_ = """
        for i = 1:30, j = 100:-2:1
            println(i, j)
        end"""
        str = """
        for i in 1:30, j in 100:-2:1
            println(i, j)
        end"""
        @test fmt(str_, always_for_in = true) == str
        @test fmt(str, always_for_in = true) == str

        str_ = "[(i,j) for i=I1,j=I2]"
        str = "[(i, j) for i in I1, j in I2]"
        @test fmt(str_, always_for_in = true) == str
        @test fmt(str, always_for_in = true) == str

        str_ = "((i,j) for i=I1,j=I2)"
        str = "((i, j) for i in I1, j in I2)"
        @test fmt(str_, always_for_in = true) == str
        @test fmt(str, always_for_in = true) == str

        str_ = "[(i, j) for i = 1:2:10, j = 100:-1:10]"
        str = "[(i, j) for i in 1:2:10, j in 100:-1:10]"
        @test fmt(str_, always_for_in = true) == str
        @test fmt(str, always_for_in = true) == str
    end

    @testset "rewrite x |> f to f(x)" begin
        @test fmt("x |> f", pipe_to_function_call = true) == "f(x)"

        str_ = "var = func1(arg1) |> func2 |> func3 |> func4 |> func5"
        str = "var = func5(func4(func3(func2(func1(arg1)))))"
        @test fmt(str_, pipe_to_function_call = true) == str
        @test fmt(str_, pipe_to_function_call = true, margin = 1) == fmt(str)
    end

    @testset "function shortdef to longdef" begin
        str_ = "foo(a) = bodybodybody"
        str = """
        function foo(a)
            bodybodybody
        end"""
        @test fmt(str_, 4, length(str_), short_to_long_function_def = true) == str_
        @test fmt(str_, 4, length(str_) - 1, short_to_long_function_def = true) == str

        str_ = "foo(a::T) where {T} = bodybodybodybodybodybodyb"
        str = """
        function foo(a::T) where {T}
            bodybodybodybodybodybodyb
        end"""
        @test fmt(str_, 4, length(str_), short_to_long_function_def = true) == str_
        @test fmt(str_, 4, length(str_) - 1, short_to_long_function_def = true) == str

        str_ = "foo(a::T)::R where {T} = bodybodybodybodybodybodybody"
        str = """
        function foo(a::T)::R where {T}
            bodybodybodybodybodybodybody
        end"""
        @test fmt(str_, 4, length(str_), short_to_long_function_def = true) == str_
        @test fmt(str_, 4, length(str_) - 1, short_to_long_function_def = true) == str
    end

    @testset "always use return" begin
        str_ = "foo(a) = bodybodybody"
        str = """
        function foo(a)
            return bodybodybody
        end"""
        @test fmt(
            str_,
            4,
            length(str_) - 1,
            short_to_long_function_def = true,
            always_use_return = true,
        ) == str

        str_ = """
        function foo()
            expr1
            expr2
        end"""
        str = """
        function foo()
            expr1
            return expr2
        end"""
        @test fmt(str_, 4, length(str_) - 1, always_use_return = true) == str

        str_ = """
        macro foo()
            expr1
            expr2
        end"""
        str = """
        macro foo()
            expr1
            return expr2
        end"""
        @test fmt(str_, 4, length(str_) - 1, always_use_return = true) == str

        str_ = """
        map(arg1, arg2) do x, y
            expr1
            expr2
        end"""
        str = """
        map(arg1, arg2) do x, y
            expr1
            return expr2
        end"""
        @test fmt(str_, 4, length(str_) - 1, always_use_return = true) == str

        str = """
        function foo()
            @macrocall(expr2)
        end"""
        @test fmt(str, 4, 92, always_use_return = true) == str

        str = """
        function foo()
            @macroblock expr2
        end"""
        @test fmt(str, 4, 92, always_use_return = true) == str

        str = """
        function foo()
            for i = 1:10
                println(i)
            end
        end"""
        @test fmt(str, 4, 92, always_use_return = true) == str

        str = """
        function f(a)
            if a > 0
                return -1
            else
                return 1
            end
        end"""
        @test fmt(str, 4, 92, always_use_return = true) == str
    end

    @testset "whitespace in keyword arguments" begin
        str_ = "f(; a = b)"
        str = "f(; a=b)"
        @test fmt(str_, 4, 92, whitespace_in_kwargs = false) == str

        str = "f(; a!) = a!"
        @test fmt(str, 4, 92, whitespace_in_kwargs = false) == str

        # issue 242
        str_ = "f(a, b! = 1; c! = 2, d = 3, e! = 4)"
        str = "f(a, b! = 1; c! = 2, d=3, e! = 4)"
        @test fmt(str_, 4, 92, whitespace_in_kwargs = false) == str

        str_ = "( k1 =v1,  k2! = v2)"
        str = "(k1=v1, k2! = v2)"
        @test fmt(str_, 4, 80, style = YASStyle(), whitespace_in_kwargs = false) == str
        @test fmt(str_, 4, 80, style = DefaultStyle(), whitespace_in_kwargs = false) == str

        str_ = "( k1 =v1,  k2! = v2)"
        str = "(k1 = v1, k2! = v2)"
        @test fmt(str_, 4, 80, style = YASStyle(), whitespace_in_kwargs = true) == str
        @test fmt(str_, 4, 80, style = DefaultStyle(), whitespace_in_kwargs = true) == str
    end

    @testset "annotate untyped fields with `Any`" begin
        str = """
        struct name
            arg::Any
        end"""

        str_ = """
        struct name
            arg
        end"""
        @test fmt(str_) == str

        str_ = """
        struct name
        arg
        end"""
        @test fmt(str_) == str

        str_ = """
        struct name
                arg
            end"""
        @test fmt(str_) == str

        t = run_pretty(str_, 80)
        @test length(t) == 12

        str = """
        mutable struct name
            reallylongfieldname::Any
        end"""

        str_ = """
        mutable struct name
            reallylongfieldname
        end"""
        @test fmt(str_) == str

        str_ = """
        mutable struct name
        reallylongfieldname
        end"""
        @test fmt(str_) == str

        str_ = """
        mutable struct name
                reallylongfieldname
            end"""
        @test fmt(str_) == str

        t = run_pretty(str_, 80)
        @test length(t) == 28

        str = """
        struct name
            arg
        end"""

        str_ = """
        struct name
            arg
        end"""
        @test fmt(str_, annotate_untyped_fields_with_any = false) == str

        str_ = """
        struct name
        arg
        end"""
        @test fmt(str_, annotate_untyped_fields_with_any = false) == str

        str_ = """
        struct name
                arg
            end"""
        @test fmt(str_, annotate_untyped_fields_with_any = false) == str

        t = run_pretty(
            str_,
            opts = Options(max_margin = 80, annotate_untyped_fields_with_any = false),
        )
        @test length(t) == 11

        str = """
        mutable struct name
            reallylongfieldname
        end"""

        str_ = """
        mutable struct name
            reallylongfieldname
        end"""
        @test fmt(str_, annotate_untyped_fields_with_any = false) == str

        str_ = """
        mutable struct name
        reallylongfieldname
        end"""
        @test fmt(str_, annotate_untyped_fields_with_any = false) == str

        str_ = """
        mutable struct name
                reallylongfieldname
            end"""
        @test fmt(str_, annotate_untyped_fields_with_any = false) == str

        t = run_pretty(
            str_,
            opts = Options(max_margin = 80, annotate_untyped_fields_with_any = false),
        )
        @test length(t) == 23
    end

    @testset "format docstrings - basic" begin
        str = """
        \"""
        doc
        \"""
        function f()
            20
        end"""
        t = run_pretty(str, 80)
        @test length(t) == 12

        normalized = """
        \"""
        doc
        \"""
        function f()
            20
        end"""

        str = """
        \"""doc
        \"""
        function f()
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == normalized

        str = """
        \"""
        doc\"""
        function f()
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == normalized

        str = """
        \"""doc\"""
        function f()
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == normalized

        str = """
        "doc
        "
        function f()
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == normalized

        str = """
        "
        doc"
        function f()
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == normalized

        str = """
        "doc"
        function f()
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == normalized

        # test aligning to function identation
        str_ = """
            "doc"
        function f()
            20
        end"""
        str = """
        "doc"
        function f()
            20
        end"""
        @test fmt(str_) == str
        @test fmt(str_, format_docstrings = true) == normalized

        str = """\"""
                 doc for Foo
                 \"""
                 Foo"""
        @test fmt(str) == str
        t = run_pretty(str, 80)
        @test length(t) == 11

        str = """
        \"""
        doc
        \"""
        function f()    #  comment
            20
        end"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        # Issue 157
        str = raw"""
        @doc \"""
           foo()
        \"""
        foo() = bar()"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        str = raw"""
        @doc doc\"""
           foo()
        \"""
        foo() = bar()"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        str = raw"""@doc "doc for foo" foo"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        str = raw"""@doc \"""doc for foo\""" foo"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        str = raw"""@doc doc\"""doc for foo\""" foo()"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        str = raw"""@doc foo"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str

        # issue 160
        str = """
        module MyModule

        import Markdown: @doc_str

        @doc doc\"""
            foo()
        \"""
        foo() = bar()

        end # module"""
        @test fmt(str) == str
        @test fmt(str, format_docstrings = true) == str
    end

    @testset "format docstrings - with code" begin
        unformatted = """
        \"""
        This is a docstring

        ```@example
        a =  1
         b  = 2
         a + b
        ```

        ```jldoctest
        a =  1
        b  = 2
        a + b

        # output

        3
        ```

        ```jldoctest
        julia> a =  1
        1

        julia> b  = 2;

        julia>  a + b
        3

        julia> function test(x)
               x + 1
               x + 2
               end;
        ```
        \"""
        function test(x) x end"""

        formatted = """
        \"""
        This is a docstring

        ```@example
        a = 1
        b = 2
        a + b
        ```

        ```jldoctest
        a = 1
        b = 2
        a + b

        # output

        3
        ```

        ```jldoctest
        julia> a = 1
        1

        julia> b = 2;

        julia> a + b
        3

        julia> function test(x)
                   x + 1
                   x + 2
               end;

        ```
        \"""
        function test(x)
            x
        end"""
        @test fmt(unformatted, format_docstrings = true) == formatted
    end

    @testset "format docstrings - Multi-line indented code-blocks" begin
        unformatted = """
        \"""
            fmt(
            )
        \"""
        function fmt() end"""

        formatted = """
        \"""
            fmt(
            )
        \"""
        function fmt() end"""
        @test fmt(unformatted, format_docstrings = true) == formatted
    end
    @testset "format docstrings - Empty line in docstring" begin
        unformatted = """
        \"""

        \"""
        function test() end"""

        formatted = """
        \"""
        \"""
        function test() end"""
        @test fmt(unformatted, format_docstrings = true) == formatted
    end

    @testset "format docstrings - Indented docstring" begin
        unformatted = """
        begin
            \"""
            Indented docstring

            with multiple paragraphs
            \"""
            indented_item
        end"""
        formatted = """
        begin
            \"""
            Indented docstring

            with multiple paragraphs
            \"""
            indented_item
        end"""
        @test fmt(unformatted, format_docstrings = true) == formatted

        short = """
        begin
            "Short docstring"
            item
        end
        """
        short_formatted = """
        begin
            \"""
            Short docstring
            \"""
            item
        end
        """
        @test fmt(short, format_docstrings = true) == short_formatted
    end
end


text:
using Statistics
using Distributions

# R = mean ranks, n = number of datasets, k = number of models
friedman_test_statistic(R::Vector,n::Int,k::Int) = 12*n/(k*(k+1))*(sum(R.^2) - k*(k+1)^2/4)
# k = number of models
crit_chisq(α::Real, df::Int) = quantile(Chisq(df), 1-α)
friedman_critval(α::Real, k::Int) = crit_chisq(α/2, k-1)

# k = number of groups, df = (N - k) = (total samples - k)
crit_srd(α::Real, k::Real, df::Real) =
	(isnan(k) | isnan(df)) ? NaN : quantile(StudentizedRange(df, k), 1-α)
nemenyi_cd(k::Int, n::Int, α::Real) = sqrt(k*(k+1)/(6*n))*crit_srd(α, k, Inf)/sqrt(2)

# this checks out with our kdd18 paper
nd = 35
nm = 6
α = 0.05
nemenyi_cd(nm, nd, α) # = 1.275

# this is for our new paper
nd = 40
nm = 23
α = 0.05
ncv5 = nemenyi_cd(nm, nd, α) # = 5.48
α = 0.1
ncv10 = nemenyi_cd(nm, nd, α) # = 5.14

include("crit_diag_functions.jl")
nd = 40
#
model_names = ["aae", "avae", "gano", "vae", "wae", "abod", "hbos", "if", "knn",
"loda", "lof", "orbf", "osvm", "pidf", "maf", "rnvp", "sptn", "fmgn", "gan",
"mgal", "dsvd", "vaek", "vaeo"]
model_names_gray = ["aae", "avae", "gano", "vae", "wae", "\\textcolor{gray}{abod}",
"\\textcolor{gray}{hbos}", "\\textcolor{gray}{if}", "\\textcolor{gray}{knn}",
"\\textcolor{gray}{loda}", "\\textcolor{gray}{lof}", "\\textcolor{gray}{orbf}",
"\\textcolor{gray}{osvm}", "\\textcolor{gray}{pidf}", "maf", "rnvp", "sptn", "fmgn",
"gan", "mgal", "dsvd", "vaek", "vaeo"]
model_ranks = [8.8, 13.4, 10.2, 8.0, 6.3, 11.3, 14.4, 14.2, 8.0, 15.6, 11.3, 8.1, 2.8,
	14.1, 8.8, 8.7, 10.3, 11.3, 11.8, 21.4, 16.0, 10.9, 6.7]
plot5 = ranks2tikzcd(model_ranks, model_names_gray,
	nemenyi_cd(length(model_ranks), nd, 0.05);
	scale=0.45, vpad=0.6, levelwidth=0.03, vdist=0.15)
string2file("tabular_auc_auc_cdd5.tex", plot5)

plot10 = ranks2tikzcd(model_ranks, model_names_gray,
	nemenyi_cd(length(model_ranks), nd, 0.1);
	scale=0.45, vpad=0.6, levelwidth=0.03, vdist=0.15)
string2file("tabular_auc_auc_cdd10.tex", plot10)

#
model_ranks_tpr5 = [10.5, 12.5, 9.9, 7.6, 8.2, 12.0, 14.7, 15.9, 9.8, 17.9, 10.9, 10.2,
	3.0, 13.9, 10.0, 9.8, 11.8, 9.6, 10.0, 18.9, 15.9, 11.2, 6.4]
plot5tpr = ranks2tikzcd(model_ranks_tpr5, model_names_gray,
	nemenyi_cd(length(model_ranks), nd, 0.05);
	scale=0.55, vpad=0.5, levelwidth=0.03, vdist=0.13)
string2file("tabular_tpr5_tpr5_cdd5.tex", plot5tpr)

plot10tpr = ranks2tikzcd(model_ranks_tpr5, model_names_gray,
	nemenyi_cd(length(model_ranks), nd, 0.1);
	scale=0.55, vpad=0.5, levelwidth=0.03, vdist=0.13)
string2file("tabular_tpr5_tpr5_cdd10.tex", plot10tpr)

# image data
nd = 40
img_model_names = ["aae", "gano", "skip", "vae", "wae", "\\textcolor{gray}{knn}",
"\\textcolor{gray}{osvm}", "fano", "fmgn", "dsvd", "vaek", "vaeo"]
img_model_ranks = [3.0, 6.2, 11.8, 6.0, 8.2, 8.0, 7.2, 5.8, 3.0, 6.2, 5.0, 2.5]
plot5img = ranks2tikzcd(img_model_ranks,img_model_names,
	nemenyi_cd(length(img_model_ranks), nd, 0.05))
string2file("images_auc_auc_cdd5.tex", plot5img)

plot10img = ranks2tikzcd(img_model_ranks,img_model_names,
	nemenyi_cd(length(img_model_ranks), nd, 0.1))
string2file("images_auc_auc_cdd10.tex", plot10img)

plot10img = ranks2tikzcd(img_model_ranks,img_model_names,
	nemenyi_cd(length(img_model_ranks), 4, 0.1))
string2file("images_auc_auc_cdd10_4_datasets.tex", plot10img)

# img data per class
include("crit_diag_functions.jl")
nd = 40
img_model_names = ["aae", "gano", "skip", "vae", "wae", "\\textcolor{gray}{knn}",
"\\textcolor{gray}{osvm}", "fano", "fmgn", "dsvd", "vaek", "vaeo"]
img_model_ranks_per_class = [5.0, 6.8, 10.6, 5.8, 6.2, 6.3, 7.4, 6.4,
	4.6, 6.5, 3.9, 3.2]
plot10img = ranks2tikzcd(img_model_ranks_per_class,img_model_names,
	nemenyi_cd(length(img_model_ranks_per_class), nd, 0.1))
string2file("images_auc_auc_cdd10_per_class.tex", plot10img)

ext: jl


text:
######### StanSample program example  ###########

using StanSample

ProjDir = dirname(@__FILE__)
cd(ProjDir)

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

tmpdir = ProjDir*"/tmp"

stanmodel = SampleModel("bernoulli", bernoullimodel, tmpdir=tmpdir);

println("\nTest printing of a stanmodel\n")


text:
module DataLoaders

export DataLoader, RandomSampler

include("sampler.jl")
include("loader.jl")

end # module


text:
abstract type Wishart3x3{T} <: AbstractMatrix{T} end

struct InverseWishart{T} <: Wishart3x3{T}
    data::NTuple{8,VE{T}}
end
struct CholInvWishart{T} <: Wishart3x3{T}
    data::NTuple{8,T}
end
struct RevCholWishart{T} <: Wishart3x3{T}
    data::NTuple{8,T}
end

const WishartFactor{T} = Union{CholInvWishart{T},RevCholWishart{T}}

@inline InverseWishart(x::Vararg{T,8}   ) where T = InverseWishart{T}(VE.(x))
@inline InverseWishart{T}(x::Vararg{T,8}) where T = InverseWishart{T}(VE.(x))
# @generated function CholInvWishart(x::Vararg{T,8}) where T
#     quote
#         $(Expr(:meta,:inline))
#         CholInvWishart{$T}($(Expr(:tuple,[:(@inbounds x[$i]) for i ∈ 1:8]...)))
#     end
# end

@inline function LinearAlgebra.logdet(l::Union{CholInvWishart{T},RevCholWishart{T}}) where {T}
    @inbounds SLEEFPirates.log(l[1] * l[4] * l[6])
end
@inline function logdetw(l::Union{CholInvWishart{T},RevCholWishart{T}}) where {T}
    @inbounds SLEEFwrap.log(l[1] * l[4] * l[6])
end

@inline function CholInvWishart{T}(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T) where T
    x_7 = x_8 = one(T)
    CholInvWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function CholInvWishart{T}(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T,x_7::T,x_8::T) where T
    CholInvWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function CholInvWishart(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T) where T
    x_7 = x_8 = one(T)
    CholInvWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function CholInvWishart(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T,x_7::T,x_8::T) where T
    CholInvWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function RevCholWishart{T}(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T) where T
    x_7 = x_8 = one(T)
    RevCholWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function RevCholWishart{T}(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T,x_7::T,x_8::T) where T
    RevCholWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function RevCholWishart{T1}(
            x_1::SVec{W,T},x_2::SVec{W,T},x_3::SVec{W,T},x_4::SVec{W,T},x_5::SVec{W,T},x_6::SVec{W,T},x_7::SVec{W,T},x_8::SVec{W,T}
        ) where {W,T,T1}
    RevCholWishart{SVec{W,T}}((Base.Cartesian.@ntuple 8 x))
end
@inline function RevCholWishart(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T) where T
    x_7 = x_8 = one(T)
    RevCholWishart{T}((Base.Cartesian.@ntuple 8 x))
end
@inline function RevCholWishart(x_1::T,x_2::T,x_3::T,x_4::T,x_5::T,x_6::T,x_7::T,x_8::T) where T
    RevCholWishart{T}((Base.Cartesian.@ntuple 8 x))
end
# @inline RevCholWishart{T}(x::Vararg{T,8}) where T = RevCholWishart{T}(x)
# @inline RevCholWishart(x::Vararg{T,6}   ) where T = RevCholWishart{T}(x,T(1),T(1))
# @inline RevCholWishart{T}(x::Vararg{T,6}) where T = RevCholWishart{T}(x,T(1),T(1))
@inline CholInvWishart{T}(x::Vararg{SVec{W,T},8}) where {W,T} = CholInvWishart{SVec{W,T}}(x)

inv_type(::Type{CholInvWishart{T}}) where T = RevCholWishart{T}
inv_type(::Type{RevCholWishart{T}}) where T = CholInvWishart{T}
# @inline function SArray{Tuple{S},T,N,L}(vs::Vararg{SIMDPirates.SVec{W,T},L}) where {S,T,N,L,W}
#     SArray{S}(vs...)
# end
# @inline Base.size(::Wishart2x2) = (2,2)
# @inline Base.length(::Wishart2x2) = 4
@inline Base.size(::Wishart3x3) = (3,3)
@inline Base.length(::Wishart3x3) = 9
# ScatteredArrays.type_length(::Type{Wishart2x2}) = 3
# ScatteredArrays.type_length(::Type{<:Wishart3x3}) = 6
# ScatteredArrays.type_length(::Type{CholInvWishart{T}}) where T = 6

@inline extractval(x::Core.VecElement{T}) where T = x.value
@inline extractval(x) = x
@inline Base.getindex(w::Wishart3x3, i::Integer) = extractval(@inbounds w.data[i])

@inline Base.getindex(A::Wishart3x3, ::LinearStorage, i::Integer) = @inbounds A.data[i]

@inline function Base.getindex(w::Wishart3x3, i::CartesianIndex{2})
    w[i[1], i[2]]
end
@inline function Base.getindex(iw::InverseWishart, i::Integer, j::Integer)
    i, j = minmax(i, j)
    @boundscheck (j > 3 || i < 1) && throwboundserror()
    im1 = i - 1
    @inbounds iw.data[3im1 - ( (im1*i) >> 1) + j].value
end
@inline function Base.getindex(w::WishartFactor{T}, j::Integer, i::Integer) where T
    j < i && return zero(T)
    @boundscheck (j > 3 || i < 1) && throwboundserror()
    im1 = i - 1
    @inbounds extractval(w.data[3im1 - ( (im1*i) >> 1) + j])
end
function Base.show(io::IO, ::MIME"text/plain", x::InverseWishart{T}) where T
    denom = x[7] > 2 ? x[7]-2 : one(T)
    println(io, "InverseWishart{$T} with ν = $(x[7]+2)")
    Base.print_matrix(io, x / denom)
end
function Base.show(io::IO, ::MIME"text/plain", x::CholInvWishart{T}) where T
    denom = x[7] > 2 ? sqrt(x[7]-2) : one(T)
    println(io, "Cholesky factor of InverseWishart{$T} with ν = $(x[7]+2)")
    Base.print_matrix(io, x / denom)
end

@inline function InverseWishart(x::T, y::T, z::T) where T
    InverseWishart(x * x, x * y, x * z, y * y, y * z, z * z, one(T), one(T))
end

@inline function Base.:+(a::W, b::W) where W <: Wishart3x3
    W(SIMDPirates.vadd(a.data, b.data))
end
@inline function Base.:*(L::WishartFactor{T}, x::SVector{3,T}) where {T}
    @fastmath @inbounds SVector(
        L[1]*x[1],
        L[2]*x[1] + L[4]*x[2],
        L[3]*x[1] + L[5]*x[2] + L[6]*x[3]
    )
end
@inline function Base.:*(L::WishartFactor{SVec{W,T}}, x::SVector{3,T}) where {W,T}
    @fastmath @inbounds SVector(
        L[1]*x[1],
        L[2]*x[1] + L[4]*x[2],
        L[3]*x[1] + L[5]*x[2] + L[6]*x[3]
    )
end
# L2 is treated as a lower triangular matrix
@inline function Base.:*(L1::W, L2::Union{W,NTuple{6,T}}) where {T, W <: WishartFactor{T}}
    @fastmath @inbounds W(
        L1[1]*L2[1],
        L1[2]*L2[1] + L1[4]*L2[2],
        L1[3]*L2[1] + L1[5]*L2[2] + L1[6]*L2[3],
        L1[4]*L2[4],
        L1[5]*L2[4] + L1[6]*L2[5],
        L1[6]*L2[6],
        L1[7],L1[8]
    )
end
@inline function Base.:*(L1::NTuple{6,T}, L2::W) where {T, W <: WishartFactor{T}}
    @fastmath @inbounds W(
        L1[1]*L2[1],
        L1[2]*L2[1] + L1[4]*L2[2],
        L1[3]*L2[1] + L1[5]*L2[2] + L1[6]*L2[3],
        L1[4]*L2[4],
        L1[5]*L2[4] + L1[6]*L2[5],
        L1[6]*L2[6],
        L2[7],L2[8]
    )
end
@inline function Base.:/(a::W, x::T) where {T <: Number, W <: Wishart3x3{T}}
    @fastmath @inbounds begin
        xinv = 1 / x
        W(
            a[1] * xinv, a[2] * xinv, a[3] * xinv, a[4] * xinv, a[5] * xinv, a[6] * xinv, a[7], a[8]
        )
    end
end

@inline extract_ν(iw::Wishart3x3) = iw[7]
@generated function extract_ν(iw::AbstractVector{<:Wishart3x3{T}}, ::Val{NG}) where {NG,T}
    quote
        $(Expr(:meta, :inline))
        ptr_iw = Base.unsafe_convert(Ptr{$T}, pointer(iw))
        $(Expr(:tuple, [:(unsafe_load(ptr_iw, $((ng-1)*8) + 7 ) ) for ng ∈ 1:NG]...))
    end
end
@inline extract_α(iw::Union{InverseWishart,CholInvWishart}) = iw[8]
@generated function extract_α(iw::AbstractVector{<:Union{InverseWishart{T},CholInvWishart{T}}}, ::Val{NG}) where {NG,T}
    quote
        $(Expr(:meta, :inline))
        ptr_iw = Base.unsafe_convert(Ptr{$T}, pointer(iw))
        $(Expr(:tuple, [:(unsafe_load(ptr_iw, $((ng-1)*8) + 8 ) ) for ng ∈ 1:NG]...))
    end
end

@inline function finite_sqrt(x::T, min_value = eps(T)) where {T}
    rx = @fastmath sqrt(x)
    # min_value = eps(T)
    # min_value = T(0.005) # 5
    isfinite(rx) ? rx : min_value
end

# @inline
function inv_and_cholesky!(riwv::AbstractVector{RevCholWishart{T}},
                                    ciwv::AbstractVector{CholInvWishart{T}},
                                    iwv::AbstractVector{InverseWishart{T}}) where T

    @inbounds for i ∈ eachindex(iwv)
        iw = iwv[i]
        @fastmath begin
            min_value = T(0.05) * sqrt(iw[7] - 2)
            L11 = finite_sqrt(iw[1], min_value)
            R11 = 1 / L11
            L21 = R11 * iw[2]
            L31 = R11 * iw[3]
            L22 = finite_sqrt(iw[4] - L21^2, min_value)
            R22 = 1 / L22
            L32 = R22 * (iw[5] - L21 * L31)
            L33 = finite_sqrt(iw[6] - L31^2 - L32^2, min_value)
            R33 = 1 / L33

            R21 = - R22 * L21 * R11
            R31 = - R33 * ( L31*R11 + L32*R21 )
            R32 = - R33 * L32 * R22
        end
        ciwv[i] = CholInvWishart(
            L11, L21, L31, L22, L32, L33, iw[7], iw[8]
        )
        riwv[i] = RevCholWishart(
            R11, R21, R31, R22, R32, R33, iw[7], iw[8]
        )
    end
end
@inline function triangle_inv(L::NTuple{8,T}) where T
    @inbounds @fastmath begin
        L11           = L[1]
        L21, L22      = L[2], L[4]
        L31, L32, L33 = L[3], L[5], L[6]
        R11 = 1 / L11
        R22 = 1 / L22
        R33 = 1 / L33
        R21 = - R22 * L21 * R11
        R31 = - R33 * ( L31*R11 + L32*R21 )
        R32 = - R33 * L32 * R22
    end
    @inbounds (R11, R21, R31, R22, R32, R33, L[7], L[8])
end
@inline function triangle_inv(L::NTuple{6,T}) where T
    @inbounds @fastmath begin
        L11           = L[1]
        L21, L22      = L[2], L[4]
        L31, L32, L33 = L[3], L[5], L[6]
        R11 = 1 / L11
        R22 = 1 / L22
        R33 = 1 / L33
        R21 = - R22 * L21 * R11
        R31 = - R33 * ( L31*R11 + L32*R21 )
        R32 = - R33 * L32 * R22
    end
    R11, R21, R31, R22, R32, R33
end
@inline function Base.inv(w::W) where {T,W<:WishartFactor{T}}
    inv_type(W)(triangle_inv(w.data))
end

@inline function UpperTriangle3(U::WishartFactor{T}) where T
    @inbounds UpperTriangle3(SVector{6,T}(
        U[1], U[2], U[4], U[3], U[5], U[6]
    ))
end

@inline function randinvwishartfactor(rng::AbstractRNG, rcw::RevCholWishart{T}) where T
    rcwU = UpperTriangle3(rcw)
    ν = extract_ν(rcw)
    U11 = randchi(rng, ν - 2)
    U12 = randn(rng)
    U22 = randchi(rng, ν - 1)
    U13 = randn(rng)
    U23 = randn(rng)
    U33 = randchi(rng, ν)
    U =  UpperTriangle3(
        U11, U12, U22, U13, U23, U33
    )
#    @show rcw
#    @show rcwU
#    @show U
    inv(rcwU * U)
end
function randinvwishart(rng::AbstractRNG, rcw::RevCholWishart{T}) where T
    xtx( randinvwishartfactor(rng, rcw) )
end


text:
"""
Continuous SIMP algorithm, see [TarekRay2020](@cite).
"""
@params mutable struct ContinuationSIMP <: AbstractSIMP
    simp
    result
    options
    callback
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::ContinuationSIMP) = println("TopOpt continuation SIMP algorithm")

@params struct CSIMPOptions
    p_cont
    option_cont
    reuse::Bool
    steps::Integer
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::CSIMPOptions) = println("TopOpt continuation SIMP options")

function CSIMPOptions(::Type{T} = Float64;
                        steps = 40,
                        p_gen = nothing,
                        initial_options = Nonconvex.MMAOptions(),
                        pstart = T(1),
                        pfinish = T(5),
                        reuse = false,
                        options_gen = nothing
                    ) where {T}

    if p_gen == nothing
        p_cont = PowerContinuation{T}( ; b = T(1),
                                        start = pstart,
                                        steps = steps + 1,
                                        finish = pfinish
                                    )
    else
        @assert steps == p_gen.length - 1
        p_cont = p_gen
    end

    if options_gen == nothing
        options_cont = MMAOptionsGen(steps = steps, initial_options = initial_options)
    else
        options_cont = options_gen
    end

    return CSIMPOptions(p_cont, options_cont, reuse, steps)
end

function ContinuationSIMP(
    simp::SIMP{T},
    steps::Int = 40,
    options::CSIMPOptions = CSIMPOptions(
        T,
        steps = steps,
        initial_options = deepcopy(simp.optimizer.options),
    ),
    callback = (i) -> (),
) where T
    @assert steps == options.steps
    ncells = getncells(simp.solver.problem)
    result = NewSIMPResult(T, simp.optimizer, ncells)
    return ContinuationSIMP(simp, result, options, callback)
end

# Unused functions
default_ftol(solver) = xmin / 10
function default_grtol(solver)
    @unpack xmin, problem = solver
    return xmin / 1000 * sum(problem.elementinfo.cellvolumes) * YoungsModulus(problem)
end

function MMAOptionsGen(;
    steps::Int = 40,
    initial_options = Nonconvex.MMAOptions(),
    ftol_gen = nothing,
    xtol_gen = nothing,
    grtol_gen = nothing,
    kkttol_gen = nothing,
    maxiter_gen = nothing,
    outer_maxiter_gen = nothing,
    s_init_gen = nothing,
    s_incr_gen = nothing,
    s_decr_gen = nothing,
    store_trace_gen = nothing,
    show_trace_gen = nothing,
    auto_scale_gen = nothing,
    dual_options_gen = nothing
)
    if maxiter_gen == nothing
        maxiter_cont = FixedContinuation(initial_options.maxiter, steps + 1)
    else
        @assert steps == maxiter_gen.length - 1
        maxiter_cont = maxiter_gen
    end

    if outer_maxiter_gen == nothing
        outer_maxiter_cont = FixedContinuation(initial_options.outer_maxiter, steps + 1)
    else
        @assert steps == outer_maxiter_gen.length - 1
        outer_maxiter_cont = outer_maxiter_gen
    end

    if xtol_gen == nothing
        xtol_cont = FixedContinuation(initial_options.tol.x, steps + 1)
    else
        @assert steps == xtol_gen.length - 1
        xtol_cont = xtol_gen
    end

    if ftol_gen == nothing
        ftol_cont = FixedContinuation(initial_options.tol.frel, steps + 1)
    else
        @assert steps == ftol_gen.length - 1
        ftol_cont = ftol_gen
    end

    if kkttol_gen == nothing
        kkttol_cont = FixedContinuation(initial_options.tol.kkt, steps + 1)
    else
        @assert steps == kkttol_gen.length - 1
        kkttol_cont = kkttol_gen
    end

    tol_cont = Nonconvex.Tolerance(x = xtol_cont, f = ftol_cont, kkt = kkttol_cont)

    if s_init_gen == nothing
        s_init_cont = FixedContinuation(initial_options.s_init, steps + 1)
    else
        @assert steps == s_init_gen.length - 1
        s_init_cont = s_init_gen
    end

    if s_incr_gen == nothing
        s_incr_cont = FixedContinuation(initial_options.s_incr, steps + 1)
    else
        @assert steps == s_incr_gen.length - 1
        s_incr_cont = s_incr_gen
    end
    if s_decr_gen == nothing
        s_decr_cont = FixedContinuation(initial_options.s_decr, steps + 1)
    else
        @assert steps == s_decr_gen.length - 1
        s_decr_cont = s_decr_gen
    end

    if store_trace_gen == nothing
        store_trace_cont = FixedContinuation(initial_options.store_trace, steps + 1)
    else
        @assert steps == store_trace_gen.length - 1
        store_trace_cont = store_trace_gen
    end

    if show_trace_gen == nothing
        show_trace_cont = FixedContinuation(initial_options.show_trace, steps + 1)
    else
        @assert steps == show_trace_gen.length - 1
        show_trace_cont = show_trace_gen
    end

    if auto_scale_gen == nothing
        auto_scale_cont = FixedContinuation(initial_options.auto_scale, steps + 1)
    else
        @assert steps == auto_scale_gen.length - 1
        auto_scale_cont = auto_scale_gen
    end

    if dual_options_gen == nothing
        dual_options_cont = FixedContinuation(initial_options.dual_options, steps + 1)
    else
        @assert steps == dual_options_gen.length - 1
        dual_options_cont = dual_options_gen
    end

    return MMAOptionsGen(
        maxiter_cont,
        outer_maxiter_cont,
        tol_cont,
        s_init_cont,
        s_incr_cont,
        s_decr_cont,
        store_trace_cont,
        show_trace_cont,
        auto_scale_cont,
        dual_options_cont
    )
end

function (c_simp::ContinuationSIMP)(
    x0 = copy(c_simp.simp.solver.vars),
)
    @unpack optimizer = c_simp.simp
    @unpack workspace = optimizer
    @unpack options = workspace

    update!(c_simp, 1)
    # Does the first function evaluation
    reset!(workspace, x0)
    c_simp.callback(0)
    r = c_simp.simp(workspace)

    f_x_previous = workspace.solution.prevf
    g_previous = copy(c_simp.simp.optimizer.workspace.solution.g)

    original_maxiter = options.maxiter
    for i in 1:c_simp.options.steps
        c_simp.callback(i)
        update!(c_simp, i+1)
        reset!(workspace)

        maxiter = 1
        workspace.options.maxiter = maxiter
        r = c_simp.simp(workspace)

        if !workspace.solution.convstate.converged
            while !workspace.solution.convstate.converged
                maxiter += 1
                workspace.options.maxiter = maxiter
                c_simp.simp(workspace)
            end
        end
        f_x_previous = workspace.solution.prevf
        g_previous = copy(c_simp.simp.optimizer.workspace.solution.g)

        workspace.options.maxiter = original_maxiter
    end
    c_simp.result = c_simp.simp.result
end

function update!(c_simp::ContinuationSIMP, i)
    p = c_simp.options.p_cont(i)
    setpenalty!(c_simp.simp, p)
    options = c_simp.options.option_cont(c_simp.simp.optimizer.workspace.options, i)
    c_simp.simp.optimizer.workspace.options = options
    return c_simp
end
frac(x) = 2*min(abs(x), abs(x - 1))

function undo_values!(workspace, f_x_previous, g_previous)
    workspace.solution.f = workspace.solution.prevf
    workspace.solution.prevf = f_x_previous
    workspace.solution.g .= g_previous
    workspace.solution.x .= workspace.solution.prevx
    workspace.solution.prevx .= workspace.tempx
    return workspace
end


text:
# type annotations and scope:
x = 1.0 # x is Float64
x = 1 # now x is Int
# y::Float64 = 1.0 # LoadError: "y is not defined"

function scopetest()
	println(x)
    y::Float64 = 1.0 # y must be Float64, this is not possible in global scope
    # y = 1
end

scopetest() #> 1
# println(y) #> ERROR: y not defined

# compound expressions:
x = begin
  a = 5
  2 * a
end # after: x is 10
println(a) #> a is 5
x = (a = 5; 2 * a)

max_stars_repo_path: Module 1/Chapter02/scope.jl

max_issues_repo_path: Module 1/Chapter02/scope.jl

max_forks_repo_path: Module 1/Chapter02/scope.jl


text:
function extract_rates(tree::Tree ,edge_trees::Array{EdgeTreeRates2,1}; id = 1)
    function aux(subtree, rates,ri)
        if subtree.n_nodes == 0
            #return []
        elseif length(subtree.offsprings) == 0
            pushfirst!(rates,subtree.attributes[id]*ri)
        else
            aux(subtree.offsprings[2],rates,ri)
            aux(subtree.offsprings[1],rates,ri)
            pushfirst!(rates,subtree.attributes[id]*ri)
        end
    end

    rates =  Array{Float64,1}(undef,0)
    for i in length(edge_trees):-1:1
        aux(edge_trees[i].tree, rates, edge_trees[i].stem_rate[1])
    end
    pushfirst!(rates, tree.attributes[1])
    return rates
end

function extract_relative_rates(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}, rates::Array{Float64,1}; id = 1)

    function aux(subtree, parent_rate, x)
        if subtree.n_nodes > 0
            if length(subtree.offsprings) == 0
                pushfirst!(x,log(subtree.attributes[id]) - parent_rate)
            else
                log_rate = deepcopy(log(subtree.attributes[id]))
                aux(subtree.offsprings[2], log_rate, x)
                aux(subtree.offsprings[1], log_rate, x)
                pushfirst!(x,log(subtree.attributes[id]) - parent_rate)
            end
        end
    end

    relative_rates = Array{Float64,1}(undef,0)

    i = 0
    for et in edge_trees
        i+= 1
        parent_edge = et.parent_edge

        parent_rate = rates[1]
        if parent_edge > 0
            parent_rate = edge_trees[parent_edge].tip_rate * rates[parent_edge + 1]
        end

        r = log(rates[i + 1]) - log(parent_rate)
        pushfirst!(relative_rates,r)
        if length(et.tree.offsprings) > 0
            aux(et.tree.offsprings[1], 0., relative_rates)
            aux(et.tree.offsprings[2], 0., relative_rates)
        end
    end
    return relative_rates
end

function extract_relative_rates(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}; id = 1)
    function aux(subtree, parent_rate, x)
        if subtree.n_nodes > 0
            if length(subtree.offsprings) == 0
                pushfirst!(x,log(subtree.attributes[id]) - parent_rate)
            else
                log_rate = deepcopy(log(subtree.attributes[id]))
                aux(subtree.offsprings[2], log_rate, x)
                aux(subtree.offsprings[1], log_rate, x)
                pushfirst!(x,log(subtree.attributes[id]) - parent_rate)
            end
        end
    end

    relative_rates = Array{Float64,1}(undef,0)

    for et in edge_trees
        parent_edge = et.parent_edge

        parent_rate = tree.attributes[1]
        if parent_edge > 0
            parent_rate = edge_trees[parent_edge].tip_rate * edge_trees[parent_edge].stem_rate[1]
        end

        r = log(et.stem_rate[1]) - log(parent_rate)
        pushfirst!(relative_rates,r)

        if length(et.tree.offsprings) > 1#0000000
            aux(et.tree.offsprings[1], 0., relative_rates)
            aux(et.tree.offsprings[2], 0., relative_rates)
        end
    end

    return relative_rates
end

function extract_rates(tree::Tree; id=1)
    function aux(subtree, rates)
        if subtree.n_nodes == 0
            #return []
        elseif length(subtree.offsprings) == 0
            pushfirst!(rates,subtree.attributes[id])
        else
            aux(subtree.offsprings[2],rates)
            aux(subtree.offsprings[1],rates)
            pushfirst!(rates,subtree.attributes[id])
        end
    end

    rates = Array{Float64,1}(undef,0)
    aux(tree, rates)
    return(rates)
end

#=
function update_rates!(tree::Tree, rates::Array{T,1} ; id = 1) where {T<:Number}
    function aux(sub_tree, sub_rates)
        if sub_tree.n_nodes < 2
            root_rate = popfirst!(sub_rates)
            sub_tree.attributes[id] = root_rate
        else
            root_rate = popfirst!(sub_rates)
            sub_tree.attributes[id] = root_rate
            t1 = sub_tree.offsprings[1]
            aux(t1, sub_rates)
            t2 = sub_tree.offsprings[2]
            aux(t2, sub_rates)
            #return Tree([t1,t2], sub_tree.branch_length, root_rate, sub_tree.n_nodes)
        end
    end

    new_rates = deepcopy(rates)
    aux(tree, new_rates)
end

function update_rates!(tree::Tree, rate::T ; id = 1) where {T<:Number}
    function aux(sub_tree)
        if sub_tree.n_nodes < 2
            sub_tree.attributes[id] *= rate
        else
            sub_tree.attributes[id] *= rate
            t1 = sub_tree.offsprings[1]
            aux(t1)
            t2 = sub_tree.offsprings[2]
            aux(t2)
        end
    end

    aux(tree)
end
=#

function extract_tip_rates(tree::Tree, edge_trees::Array{EdgeTreeRates2,1},
    tips_id::Array{Bool,1}, rates::Array{Float64,1} ; id = 1, return_extinct = true)

    tip_rates = Array{Float64,1}(undef,0)
    for i in 1:length(edge_trees)
        if tips_id[i]
            r = sample(edge_trees[i].tip_rates) * edge_trees[i].stem_rate[1]
            push!(tip_rates,r)
        end
    end
    return tip_rates
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
function test_generic_operator_interface(op, T)
    ELT = eltype(op)
    @test promote_type(T,ELT) == ELT
    #@test coefficienttype(dest(op)) == ELT
    # This may no longer always be true:
    # @test eltype(src(op)) == ELT

    m = matrix(op)

    # Test elementwise equality
    for i in 1:size(op,1)
        for j in 1:size(op,2)
            @test m[i,j] == op[i,j]
        end
    end

    # Does the operator agree with its matrix?
    r = zeros(src(op))
    for i in eachindex(r)
        r[i] = convert(ELT, rand())
    end
    v1 = zeros(dest(op))
    apply!(op, v1, r)
    if ndims(r)==1
        v2 = m*r
    else
        v2 = m*reshape(r,length(r))
    end
    @test maximum(abs.(v1[:]-v2[:])) < 10*sqrt(eps(T))

    # Verify claim to be in-place
    if isinplace(op)
        v3 = copy(r)
        apply_inplace!(op, v3)
        @test maximum(abs.(v3-v2)) < 10*sqrt(eps(T))
    end

    # Verify that coef_src is not altered when applying out-of-place
    r2 = copy(r)
    v = zeros(dest(op))
    apply!(op, v, r)
    @test maximum(abs.(r-r2)) < eps(T)

    # Test claim to be diagonal
    if isdiag(op)
        for i in 1:size(op,1)
            for j in 1:size(op,2)
                if i != j
                    @test abs(m[i,j]) < 10*sqrt(eps(T))
                end
            end
        end
    end

    # Verify diagonal entries
    d = diag(op)
    for i in 1:min(size(op,1),size(op,2))
        d[i] ≈ m[i,i]
        diag(op, i) ≈ m[i,i]
    end

    ## # Verify eltype promotion
    ## T2 = widen(eltype(op))
    ## if T2 != eltype(op)
    ##     op2 = promote_eltype(op, T2)
    ##     ELT2 = eltype(op2)
    ##     # ELT2 may not equal T, but it must be wider.
    ##     # For example, when T2 is BigFloat, ELT2 could be Complex{BigFloat}
    ##     @test promote_type(T2, ELT2) == ELT2
    ##     @test coefficienttype(dest(op2)) == ELT2
    ## end

    # Verify inverse
    try
        inv_op = inv(op)
        if ~(typeof(inv_op) <: OperatorInverse)
            @test src(inv_op) == dest(op)
            @test dest(inv_op) == src(op)
            m2 = matrix(inv_op)
            I1 = eye(ELT, length(src(op)))
            I2 = eye(ELT, length(dest(op)))
            @test maximum(abs.(m2*m - I1)) < 10*sqrt(eps(T))
            @test maximum(abs.(m*m2 - I2)) < 10*sqrt(eps(T))
        end
    catch MethodError
        # Inverse was not defined
    end

    # Verify transpose
    try
        ct_op = adjoint(opt)
        if ~(typeof(ct_op) <: OperatorTranspose)
            @test src(ct_op) == dest(op)
            @test dest(ct_op) == src(op)
            m2 = matrix(ct_op)
            @test maximum(abs.(m' - m2)) < 10*sqrt(eps(T))
        end
    catch MethodError
        # transpose was not defined
    end
end


text:
macro schema(ex)
    schema_m(ex) |> esc
end

macro schema(typealias::String, ex)
    schema_m(ex, typealias) |> esc
end

macro validate(ex)
    error("@validate must be called inside @schema")
end

macro check(ex)
    error("@check must be called inside @schema")
end

function schema_m(ex::Expr, typealias=nothing)
    def = JLKwStruct(ex, typealias)

    has_plain_constructor(def) && return
    funcs, validate_calls, checker_calls, msgs = split_validate_fn(def)

    ret = Expr(:block)
    # generate name
    for fn in funcs
        push!(ret.args, codegen_ast(fn))
    end

    body = Expr(:block)
    for call in validate_calls
        lhs = Expr(:tuple, call.args[2:end]...)
        push!(body.args, :($lhs = $call))
    end

    for (call, msg) in zip(checker_calls, msgs)
        vars = filter(x->isa(x, Symbol), call.args[2:end])
        got = Expr(:call, :string)
        for v in vars
            push!(got.args, "$v=")
            push!(got.args, v)
        end
        push!(body.args, :($call || throw($ValidationError($msg, $got))))
    end

    fields = [f.name for f in def.fields]
    typevars = name_only.(def.typevars)
    new = isempty(def.typevars) ? :new : Expr(:curly, :new, typevars...)
    push!(body.args, Expr(:call, new, fields...))

    constructor = JLFunction(;
        name = struct_name_plain(def),
        args = fields,
        whereparams = isempty(def.typevars) ? nothing : typevars,
        body = body,
    )

    push!(def.constructors, constructor)
    push!(ret.args, Configurations.codegen_option_type(def))
    return ret
end

function split_validate_or_check(funcs, calls, msg, body)
    if is_function(body)
        fn = JLFunction(body)

        if msg !== nothing
            if fn.name === nothing
                push!(msg, string(prettify(fn.body)))
            else
                push!(msg, fn.name)
            end
        end

        fn.name = gensym(:validate)
        push!(funcs, fn)
        push!(calls, Expr(:call, fn.name, name_only.(fn.args)...))
    else
        body isa Expr && body.head === :call ||
            throw(ArgumentError("expect function definition or function call"))

        push!(calls, body)
        msg === nothing || push!(msg, string(prettify(body)))
    end
    return
end

function split_validate_fn(def::JLKwStruct)
    msg, misc, funcs, validate_calls, checker_calls = [], [], [], [], []

    for each in def.misc
        @switch each begin
            @case Expr(:macrocall, Symbol("@validate"), _, body)
                split_validate_or_check(funcs, validate_calls, nothing, body)
            @case Expr(:macrocall, Symbol("@check"), _, body)
                split_validate_or_check(funcs, checker_calls, msg, body)
            @case _
                push!(misc, each)
        end
    end
    def.misc = misc
    return funcs, validate_calls, checker_calls, msg
end


text:
# Legendre functions

# 4π normalized
# Orthonormalized
# Schmidt semi-normalized
# Unnormalized
# Utilities

"""
    PlmBar!(p::AbstractVector{Cdouble},
            lmax::Integer,
            z::Union{AbstractFloat,Integer};
            csphase::Integer=1,
            cnorm::Integer=0,
            exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the 4π-normalized associated Legendre functions.

See also: [`PlmBar`](@ref), [`PlmBar_d1!`](@ref), [`PlBar!`](@ref),
[`PlmON!`](@ref), [`PlmSchmidt!`](@ref), [`PLegendreA!`](@ref),
[`PlmIndex`](@ref)
"""
function PlmBar! end

"""
    p = PlmBar(lmax::Integer,
               z::Union{AbstractFloat,Integer};
               csphase::Integer=1,
               cnorm::Integer=0,
               exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the 4π-normalized associated Legendre functions.

See also: [`PlmBar!`](@ref), [`PlmBar_d1`](@ref), [`PlBar`](@ref),
[`PlmON`](@ref), [`PlmSchmidt`](@ref), [`PLegendreA`](@ref),
[`PlmIndex`](@ref)
"""
function PlmBar end

"""
    PlmBar_d1!(p::AbstractVector{Cdouble},
               dp::AbstractVector{Cdouble},
               lmax::Integer,
               z::Union{AbstractFloat,Integer};
               csphase::Integer=1,
               cnorm::Integer=0,
               exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the 4π-normalized associated Legendre functions and first
derivatives.

See also: [`PlmBar_d1`](@ref), [`PlmBar_d1`](@ref),
[`PlBar_d1!`](@ref), [`PlmON_d1!`](@ref), [`PlmSchmidt_d1!`](@ref),
[`PLegendreA_d1!`](@ref), [`PlmIndex`](@ref)
"""
function PlmBar_d1! end

"""
    p, dp = PlmBar_d1(lmax::Integer,
                      z::Union{AbstractFloat,Integer};
                      csphase::Integer=1,
                      cnorm::Integer=0,
                      exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the 4π-normalized associated Legendre functions and first
derivatives.

See also: [`PlmBar`](@ref), [`PlmBar_d1!`](@ref), [`PlBar_d1`](@ref),
[`PlmON_d1`](@ref), [`PlmSchmidt_d1`](@ref),
[`PLegendreA_d1`](@ref), [`PlmIndex`](@ref)
"""
function PlmBar_d1 end

"""
    PlmON!(p::AbstractVector{Cdouble},
           lmax::Integer,
           z::Union{AbstractFloat,Integer};
           csphase::Integer=1,
           cnorm::Integer=0,
           exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the orthonormalized associated Legendre functions.

See also: [`PlmON`](@ref), [`PlmON_d1!`](@ref), [`PlON!`](@ref),
[`PlmBar!`](@ref), [`PlmSchmidt!`](@ref), [`PLegendreA!`](@ref),
[`PlmIndex`](@ref)
"""
function PlmON! end

"""
    p = PlmON(lmax::Integer,
              z::Union{AbstractFloat,Integer};
              csphase::Integer=1,
              cnorm::Integer=0,
              exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the orthonormalized associated Legendre functions.

See also: [`PlmON!`](@ref), [`PlmON_d1`](@ref), [`PlON`](@ref),
[`PlmBar`](@ref), [`PlmSchmidt`](@ref), [`PLegendreA`](@ref),
[`PlmIndex`](@ref)
"""
function PlmON end

"""
    PlmON_d1!(p::AbstractVector{Cdouble},
              dp::AbstractVector{Cdouble},
              lmax::Integer,
              z::Union{AbstractFloat,Integer};
              csphase::Integer=1,
              cnorm::Integer=0,
              exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the orthonormalized associated Legendre functions and
first derivatives.

See also: [`PlmON_d1`](@ref), [`PlmON_d1`](@ref), [`PlON_d1!`](@ref),
[`PlmBar_d1!`](@ref), [`PlmSchmidt_d1!`](@ref),
[`PLegendreA_d1!`](@ref), [`PlmIndex`](@ref)
"""
function PlmON_d1! end

"""
    p, dp = PlmON_d1(lmax::Integer,
                     z::Union{AbstractFloat,Integer};
                     csphase::Integer=1,
                     cnorm::Integer=0,
                     exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the orthonormalized associated Legendre functions and first
derivatives.

See also: [`PlmON`](@ref), [`PlmON_d1!`](@ref), [`PlON_d1`](@ref),
[`PlmBar_d1`](@ref), [`PlmSchmidt_d1`](@ref), [`PLegendreA_d1`](@ref),
[`PlmIndex`](@ref)
"""
function PlmON_d1 end

"""
    PlmSchmidt!(p::AbstractVector{Cdouble},
                lmax::Integer,
                z::Union{AbstractFloat,Integer};
                csphase::Integer=1,
                cnorm::Integer=0,
                exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the Schmidt-normalized associated Legendre functions.

See also: [`PlmSchmidt`](@ref), [`PlmSchmidt_d1!`](@ref),
[`PlSchmidt!`](@ref), [`PlmBar!`](@ref), [`PlmON!`](@ref),
[`PLegendreA!`](@ref), [`PlmIndex`](@ref)
"""
function PlmSchmidt! end

"""
    p = PlmSchmidt(lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   csphase::Integer=1,
                   cnorm::Integer=0,
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the Schmidt-normalized associated Legendre functions.

See also: [`PlmSchmidt!`](@ref), [`PlmSchmidt_d1`](@ref),
[`PlSchmidt`](@ref), [`PlmBar`](@ref), [`PlmON`](@ref),
[`PLegendreA`](@ref), [`PlmIndex`](@ref)
"""
function PlmSchmidt end

"""
    PlmSchmidt_d1!(p::AbstractVector{Cdouble},
                   dp::AbstractVector{Cdouble},
                   lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   csphase::Integer=1,
                   cnorm::Integer=0,
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the Schmidt-normalized associated Legendre functions and
first derivatives.

See also: [`PlmSchmidt_d1`](@ref), [`PlmSchmidt_d1`](@ref),
[`PlSchmidt_d1!`](@ref), [`PlmBar_d1!`](@ref), [`PlmON_d1!`](@ref),
[`PLegendreA_d1!`](@ref), [`PlmIndex`](@ref)
"""
function PlmSchmidt_d1! end

"""
    p, dp = PlmSchmidt_d1(lmax::Integer,
                          z::Union{AbstractFloat,Integer};
                          csphase::Integer=1,
                          cnorm::Integer=0,
                          exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the Schmidt-normalized associated Legendre functions and
first derivatives.

See also: [`PlmSchmidt`](@ref), [`PlmSchmidt_d1!`](@ref),
[`PlSchmidt_d1`](@ref), [`PlmBar_d1`](@ref), [`PlmON_d1`](@ref),
[`PLegendreA_d1`](@ref), [`PlmIndex`](@ref)
"""
function PlmSchmidt_d1 end

for op in [:PlmBar, :PlmON, :PlmSchmidt]
    op! = Symbol(op, "!")
    op_d1 = Symbol(op, "_d1")
    op_d1! = Symbol(op, "_d1!")
    @eval begin
        #
        export $op!
        function $op!(p::AbstractVector{Cdouble}, lmax::Integer,
                      z::Union{AbstractFloat,Integer}; csphase::Integer=1,
                      cnorm::Integer=0,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert cnorm ∈ (0, 1)
            @assert length(p) ≥ (lmax + 1) * (lmax + 2) ÷ 2
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Cint, Cdouble, Ref{Cint}, Ref{Cint},
                   Ref{Cint}), p, lmax, z, csphase, cnorm, exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p
        end
        #
        export $op
        function $op(lmax::Integer, z::Union{AbstractFloat,Integer};
                     csphase::Integer=1, cnorm::Integer=0,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            $op!(p, lmax, z; csphase=csphase, cnorm=cnorm,
                 exitstatus=exitstatus)
            return p
        end
        #
        export $op_d1!
        function $op_d1!(p::AbstractVector{Cdouble},
                         dp::AbstractVector{Cdouble}, lmax::Integer,
                         z::Union{AbstractFloat,Integer}; csphase::Integer=1,
                         cnorm::Integer=0,
                         exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert cnorm ∈ (0, 1)
            @assert length(p) ≥ lmax + 1
            @assert length(dp) ≥ lmax + 1
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op_d1)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Ref{Cint},
                   Ref{Cint}, Ref{Cint}), p, dp, lmax, z, csphase, cnorm,
                  exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op_d1!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p, dp
        end
        #
        export $op_d1
        function $op_d1(lmax::Integer, z::Union{AbstractFloat,Integer};
                        csphase::Integer=1, cnorm::Integer=0,
                        exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            dp = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            $op_d1!(p, dp, lmax, z; csphase=csphase, cnorm=cnorm,
                    exitstatus=exitstatus)
            return p, dp
        end
        #
    end
end

"""
    PLegendreA!(p::AbstractVector{Cdouble},
                lmax::Integer,
                z::Union{AbstractFloat,Integer};
                csphase::Integer=1,
                exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the unnormalized associated Legendre functions.

See also: [`PLegendreA`](@ref), [`PLegendreA_d1!`](@ref),
[`PLegendre!`](@ref), [`PlmBar!`](@ref), [`PlmON!`](@ref),
[`PlmSchmidt!`](@ref), [`PlmIndex`](@ref)
"""
function PLegendreA! end

"""
    p = PLegendreA(lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   csphase::Integer=1,
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the unnormalized associated Legendre functions.

See also: [`PLegendreA!`](@ref), [`PLegendreA_d1`](@ref),
[`PLegendre`](@ref), [`PlmBar`](@ref), [`PlmON`](@ref),
[`PlmSchmidt`](@ref), [`PlmIndex`](@ref)
"""
function PLegendreA end

"""
    PLegendreA_d1!(p::AbstractVector{Cdouble},
                   dp::AbstractVector{Cdouble},
                   lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   csphase::Integer=1,
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the unnormalized associated Legendre functions and first
derivatives.

See also: [`PLegendreA_d1`](@ref), [`PLegendreA_d1`](@ref),
[`PLegendre_d1!`](@ref), [`PlmBar_d1!`](@ref), [`PlmON_d1!`](@ref),
[`PlmSchmidt_d1!`](@ref), [`PlmIndex`](@ref)
"""
function PLegendreA_d1! end

"""
    p, dp = PLegendreA_d1(lmax::Integer,
                            z::Union{AbstractFloat,Integer};
                            csphase::Integer=1,
                            exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the unnormalized associated Legendre functions and first
derivatives.

See also: [`PLegendreA`](@ref), [`PLegendreA_d1!`](@ref),
[`PLegendre_d1`](@ref), [`PlmBar_d1`](@ref), [`PlmON_d1`](@ref),
[`PlmSchmidt_d1`](@ref), [`PlmIndex`](@ref)
"""
function PLegendreA_d1 end

for op in [:PLegendreA]
    op! = Symbol(op, "!")
    op_d1 = Symbol(op, "_d1")
    op_d1! = Symbol(op, "_d1!")
    @eval begin
        #
        export $op!
        function $op!(p::AbstractVector{Cdouble}, lmax::Integer,
                      z::Union{AbstractFloat,Integer}; csphase::Integer=1,
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert length(p) ≥ (lmax + 1) * (lmax + 2) ÷ 2
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Cint, Cdouble, Ref{Cint}, Ref{Cint}), p, lmax,
                  z, csphase, exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p
        end
        #
        export $op
        function $op(lmax::Integer, z::Union{AbstractFloat,Integer};
                     csphase::Integer=1,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            $op!(p, lmax, z; csphase=csphase, exitstatus=exitstatus)
            return p
        end
        #
        export $op_d1!
        function $op_d1!(p::AbstractVector{Cdouble},
                         dp::AbstractVector{Cdouble}, lmax::Integer,
                         z::Union{AbstractFloat,Integer}; csphase::Integer=1,
                         exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert csphase ∈ (1, -1)
            @assert length(p) ≥ lmax + 1
            @assert length(dp) ≥ lmax + 1
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op_d1)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Ref{Cint},
                   Ref{Cint}), p, dp, lmax, z, csphase, exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op_d1!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p, dp
        end
        #
        export $op_d1
        function $op_d1(lmax::Integer, z::Union{AbstractFloat,Integer};
                        csphase::Integer=1,
                        exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            dp = Array{Cdouble}(undef, (lmax + 1) * (lmax + 2) ÷ 2)
            $op_d1!(p, dp, lmax, z; csphase=csphase, exitstatus=exitstatus)
            return p, dp
        end
        #
    end
end

"""
    PlBar!(p::AbstractVector{Cdouble},
           lmax::Integer,
           z::Union{AbstractFloat,Integer};
           exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the 4π-normalized Legendre polynomials.

See also: [`PlBar`](@ref), [`PlBar_d1!`](@ref), [`PlmBar!`](@ref),
[`PlON!`](@ref), [`PlSchmidt!`](@ref), [`PLegendre!`](@ref)
"""
function PlBar! end

"""
    p = PlBar(lmax::Integer,
               z::Union{AbstractFloat,Integer};
               exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the 4π-normalized associated Legendre functions.

See also: [`PlBar!`](@ref), [`PlBar_d1`](@ref), [`PlmBar`](@ref),
[`PlON`](@ref), [`PlSchmidt`](@ref), [`PLegendre`](@ref)
"""
function PlBar end

"""
    PlBar_d1!(p::AbstractVector{Cdouble},
               dp::AbstractVector{Cdouble},
               lmax::Integer,
               z::Union{AbstractFloat,Integer};
               exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the 4π-normalized associated Legendre functions and first
derivatives.

See also: [`PlBar_d1`](@ref), [`PlBar_d1`](@ref),
[`PlmBar_d1!`](@ref), [`PlON_d1!`](@ref), [`PlSchmidt_d1!`](@ref),
[`PLegendre_d1!`](@ref)
"""
function PlBar_d1! end

"""
    p, dp = PlBar_d1(lmax::Integer,
                      z::Union{AbstractFloat,Integer};
                      exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the 4π-normalized associated Legendre functions and first
derivatives.

See also: [`PlBar`](@ref), [`PlBar_d1!`](@ref), [`PlmBar_d1`](@ref),
[`PlON_d1`](@ref), [`PlSchmidt_d1`](@ref), [`PLegendre_d1`](@ref)
"""
function PlBar_d1 end

"""
    PlON!(p::AbstractVector{Cdouble},
           lmax::Integer,
           z::Union{AbstractFloat,Integer};
           exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the orthonormalized associated Legendre functions.

See also: [`PlON`](@ref), [`PlON_d1!`](@ref), [`PlmON!`](@ref),
[`PlBar!`](@ref), [`PlSchmidt!`](@ref), [`PLegendre!`](@ref)
"""
function PlON! end

"""
    p = PlON(lmax::Integer,
              z::Union{AbstractFloat,Integer};
              exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the orthonormalized associated Legendre functions.

See also: [`PlON!`](@ref), [`PlON_d1`](@ref), [`PlmON`](@ref),
[`PlBar`](@ref), [`PlSchmidt`](@ref), [`PLegendre`](@ref)
"""
function PlON end

"""
    PlON_d1!(p::AbstractVector{Cdouble},
              dp::AbstractVector{Cdouble},
              lmax::Integer,
              z::Union{AbstractFloat,Integer};
              exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the orthonormalized associated Legendre functions and
first derivatives.

See also: [`PlON_d1`](@ref), [`PlON_d1`](@ref), [`PlmON_d1!`](@ref),
[`PlBar_d1!`](@ref), [`PlSchmidt_d1!`](@ref), [`PLegendre_d1!`](@ref)
"""
function PlON_d1! end

"""
    p, dp = PlON_d1(lmax::Integer,
                     z::Union{AbstractFloat,Integer};
                     exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the orthonormalized associated Legendre functions and first
derivatives.

See also: [`PlON`](@ref), [`PlON_d1!`](@ref), [`PlmON_d1`](@ref),
[`PlBar_d1`](@ref), [`PlSchmidt_d1`](@ref), [`PLegendre_d1`](@ref)
"""
function PlON_d1 end

"""
    PlSchmidt!(p::AbstractVector{Cdouble},
                lmax::Integer,
                z::Union{AbstractFloat,Integer};
                exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the Schmidt-normalized associated Legendre functions.

See also: [`PlSchmidt`](@ref), [`PlSchmidt_d1!`](@ref),
[`PlmSchmidt!`](@ref), [`PlBar!`](@ref), [`PlON!`](@ref),
[`PLegendre!`](@ref)
"""
function PlSchmidt! end

"""
    p = PlSchmidt(lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the Schmidt-normalized associated Legendre functions.

See also: [`PlSchmidt!`](@ref), [`PlSchmidt_d1`](@ref),
[`PlmSchmidt`](@ref), [`PlBar`](@ref), [`PlON`](@ref),
[`PLegendre`](@ref)
"""
function PlSchmidt end

"""
    PlSchmidt_d1!(p::AbstractVector{Cdouble},
                   dp::AbstractVector{Cdouble},
                   lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the Schmidt-normalized associated Legendre functions and
first derivatives.

See also: [`PlSchmidt_d1`](@ref), [`PlSchmidt_d1`](@ref),
[`PlmSchmidt_d1!`](@ref), [`PlBar_d1!`](@ref), [`PlON_d1!`](@ref),
[`PLegendre_d1!`](@ref)
"""
function PlSchmidt_d1! end

"""
    p, dp = PlSchmidt_d1(lmax::Integer,
                          z::Union{AbstractFloat,Integer};
                          exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the Schmidt-normalized associated Legendre functions and
first derivatives.

See also: [`PlSchmidt`](@ref), [`PlSchmidt_d1!`](@ref),
[`PlmSchmidt_d1`](@ref), [`PlBar_d1`](@ref), [`PlON_d1`](@ref),
[`PLegendre_d1`](@ref)
"""
function PlSchmidt_d1 end

"""
    PLegendre!(p::AbstractVector{Cdouble},
                lmax::Integer,
                z::Union{AbstractFloat,Integer};
                exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the unnormalized associated Legendre functions.

See also: [`PLegendre`](@ref), [`PLegendre_d1!`](@ref),
[`PLegendreA!`](@ref), [`PlBar!`](@ref), [`PlON!`](@ref),
[`PlSchmidt!`](@ref)
"""
function PLegendre! end

"""
    p = PLegendre(lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}

Compute all the unnormalized associated Legendre functions.

See also: [`PLegendre!`](@ref), [`PLegendre_d1`](@ref),
[`PLegendreA`](@ref), [`PlBar`](@ref), [`PlON`](@ref),
[`PlSchmidt`](@ref)
"""
function PLegendre end

"""
    PLegendre_d1!(p::AbstractVector{Cdouble},
                   dp::AbstractVector{Cdouble},
                   lmax::Integer,
                   z::Union{AbstractFloat,Integer};
                   exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)

Compute all the unnormalized associated Legendre functions and first
derivatives.

See also: [`PLegendre_d1`](@ref), [`PLegendre_d1`](@ref),
[`PLegendreA_d1!`](@ref), [`PlBar_d1!`](@ref), [`PlON_d1!`](@ref),
[`PlSchmidt_d1!`](@ref)
"""
function PLegendre_d1! end

"""
    p, dp = PLegendre_d1(lmax::Integer,
                            z::Union{AbstractFloat,Integer};
                            exitstatus::Union{Nothing,Ref{<:Integer}}=nothing)
    p::Vector{Cdouble}
    dp::Vector{Cdouble}

Compute all the unnormalized associated Legendre functions and first
derivatives.

See also: [`PLegendre`](@ref), [`PLegendre_d1!`](@ref),
[`PLegendre_d1!`](@ref), [`PlBar_d1`](@ref), [`PlON_d1`](@ref),
[`PlSchmidt_d1`](@ref)
"""
function PLegendre_d1 end

for op in [:PlBar, :PlON, :PlSchmidt, :PLegendre]
    op! = Symbol(op, "!")
    op_d1 = Symbol(op, "_d1")
    op_d1! = Symbol(op, "_d1!")
    @eval begin
        #
        export $op!
        function $op!(p::AbstractVector{Cdouble}, lmax::Integer,
                      z::Union{AbstractFloat,Integer};
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert length(p) ≥ lmax + 1
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Cint, Cdouble, Ref{Cint}), p, lmax, z,
                  exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p
        end
        #
        export $op
        function $op(lmax::Integer, z::Union{AbstractFloat,Integer};
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, lmax + 1)
            $op!(p, lmax, z; exitstatus=exitstatus)
            return p
        end
        #
        export $op_d1!
        function $op_d1!(p::AbstractVector{Cdouble},
                         dp::AbstractVector{Cdouble}, lmax::Integer,
                         z::Union{AbstractFloat,Integer};
                         exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            @assert length(p) ≥ lmax + 1
            @assert length(dp) ≥ lmax + 1
            exitstatus′ = Ref{Cint}()
            ccall(($(QuoteNode(op_d1)), libSHTOOLS), Cvoid,
                  (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Ref{Cint}), p, dp,
                  lmax, z, exitstatus′)
            if exitstatus === nothing
                exitstatus′[] ≠ 0 &&
                    error("$($op_d1!): Error code $(exitstatus′[])")
            else
                exitstatus[] = exitstatus′[]
            end
            return p, dp
        end
        #
        export $op_d1
        function $op_d1(lmax::Integer, z::Union{AbstractFloat,Integer};
                        exitstatus::Optional{Ref{<:Integer}}=nothing)
            @assert lmax ≥ 0
            p = Array{Cdouble}(undef, lmax + 1)
            dp = Array{Cdouble}(undef, lmax + 1)
            $op_d1!(p, dp, lmax, z; exitstatus=exitstatus)
            return p, dp
        end
        #
    end
end

export PlmIndex
"""
    index = PlmIndex(l::Integer,
                     m::Integer)
    index::Int

Compute the index of an array of Legendre functions corresponding to
degree `l` and angular order `m`.
"""
function PlmIndex(l::Integer, m::Integer)
    0 ≤ m ≤ l ||
        error("m must be greater than or equal to zero and less than or equal to l: m=$m, l=$l")
    index = ccall((:PlmIndex, libSHTOOLS), Cint, (Cint, Cint), l, m)
    return Int(index)
end


text:
export TrivialGate, IdentityGate, igate

abstract type TrivialGate{N} <: PrimitiveBlock{N} end

mat(::Type{T}, d::TrivialGate{N}) where {T,N} = IMatrix{1 << N,T}()
Base.adjoint(g::TrivialGate) = g
occupied_locs(g::TrivialGate) = ()

"""
    IdentityGate{N} <: TrivialGate{N}

The identity gate.
"""
struct IdentityGate{N} <: TrivialGate{N} end

"""
    igate(n::Int)

The constructor for identity gate.
"""
igate(n::Int) = IdentityGate{n}()


text:
module TestSharedPtr
module Wrapper
    using CxxCall
    dir = mktempdir()
    lib = joinpath(dir, "libSharedPtr")
    filepath = joinpath(dir, "SharedPtrCxx.cxx")
    eval(cxxsetup())
    eval(cxxnewfile(filepath,
    """
    #include <memory>
    #include <iostream>
    struct S {
        int value;
        int deaths;
        S(int x) : value(x), deaths(0){
        }
        ~S() {
            deaths += 1;
        }
    };

    struct ShareHolder {
        std::shared_ptr<S> s;
        ShareHolder(std::shared_ptr<S> s_) : s(s_) {};
    };
    """))

    struct S end
    CxxCall.tocxx(::Type{S}) = "S"
    @cxx lib function cxxnew(::Type{S}, value::Cint)::Ptr{S}
        """
        S* ret = new S(value);
        return ret;
        """
    end
    @cxx lib free(self::Ptr{S})::Cvoid = "delete self;"

    struct SharedPtr{T} end
    function CxxCall.tocxx(::Type{SharedPtr{T}}) where {T}
        "std::shared_ptr<$(tocxx(T))>"
    end
    for T in [S]
        @cxx lib function cxxnew(::Type{SharedPtr{T}}, obj::Ptr{T})::Ptr{SharedPtr{T}}
            """
            $(tocxx(SharedPtr{T}))* ret = new $(tocxx(SharedPtr{T}))(nullptr);
            ret->reset(obj);
            return ret;
            """
        end
        @cxx lib free(self::Ptr{SharedPtr{T}})::Cvoid = "delete self;"
        @cxx lib use_count(self::Ptr{SharedPtr{T}})::Clong = "return self->use_count();"
    end

    struct ShareHolder end
    CxxCall.tocxx(::Type{ShareHolder}) = "ShareHolder"
    @cxx lib function free(self::Ptr{ShareHolder})::Cvoid
        "delete self;"
    end
    @cxx lib function cxxnew(::Type{ShareHolder}, ps::Ptr{SharedPtr{S}})::Ptr{ShareHolder}
        "return new ShareHolder(*ps);"
    end
    @cxx lib function get_deaths(self::Ptr{ShareHolder})::Cint
        "return self->s.get()->deaths;"
    end
    @cxx lib function get_deaths(self::Ptr{SharedPtr{S}})::Cint
        "return self->get()->deaths;"
    end
    @cxx lib function get_value(self::Ptr{ShareHolder})::Cint
        "return self->s.get()->value;"
    end
    @cxx lib function use_count(self::Ptr{ShareHolder})::Clong
        "return self->s.use_count();"
    end
end#module Wrapper

using Test
import .Wrapper; const W = Wrapper
@testset "std::shared_ptr" begin
    @test !ispath(W.filepath)
    W.cxx_write_code!()
    @test isfile(W.filepath)
    libpath = W.lib * ".so"
    run(`g++ -shared -fPIC $(Wrapper.filepath) -o $libpath`)
    s = W.cxxnew(W.S, Cint(1))
    sp = W.cxxnew(W.SharedPtr{W.S}, s)
    @test W.use_count(sp) == 1
    holder = W.cxxnew(W.ShareHolder, sp)
    @test W.use_count(sp) == 2
    @test W.use_count(holder) == 2
    @test W.get_deaths(holder) == 0
    @test W.get_deaths(sp)     == 0
    holder2 = W.cxxnew(W.ShareHolder, sp)
    @test W.use_count(sp) == 3
    @test W.use_count(holder) == 3
    @test W.use_count(holder2) == 3
    W.free(sp)
    @test W.get_deaths(holder) == 0
    @test W.get_deaths(holder2) == 0
    @test W.use_count(holder) == 2
    @test W.use_count(holder2) == 2
    W.free(holder)
    @test W.use_count(holder2) == 1
    @test W.get_deaths(holder2) == 0
    W.free(holder2)
end
end#module


text:
# Tests for Decision trees/ Random forest
function t_decisiontree()

Ac = j4pr.DataGenerator.fish(20)  # classification dataset
Ar = j4pr.DataGenerator.fishr(20) # regression dataset

Wclass = [j4pr.tree(), j4pr.randomforest(), j4pr.adaboostump()] # only default arguments

Wreg = [j4pr.treer(), j4pr.randomforestr()]			# only default arguments

# Test classification
for w in Wclass
	Test.@test try
		wt1 = Ac |> w
		wt2 = j4pr.strip(Ac) |> w

		result = Ac |> wt1
		result2 = Ac |> wt2

		# Do not check decisions, they may vary based on sampling etc.
		true
	catch
		false
	end
end

# Test regression
for w in Wreg
	Test.@test try
		wt1 = Ar |> w
		wt2 = j4pr.strip(Ar) |> w

		result = Ar |> wt1
		result2 = Ar |> wt2

		# Do not check decisions, they may vary based on sampling etc.
		true
	catch
		false
	end
end

end


text:
# External packages.
import .TikzPictures: TikzPicture

import .TikZ

""" Convert TikZ document to `TikzPicture` object (from `TikzPictures.jl`).
"""
function TikzPicture(doc::TikZ.Document; preamble::String="")::TikzPicture
  data = join(sprint(TikZ.pprint, stmt) for stmt in doc.picture.stmts)
  options = join((sprint(TikZ.pprint, prop) for prop in doc.picture.props), ",")
  preamble = join([
    [ preamble, "\\usepackage{amssymb}" ];
    [ "\\usetikzlibrary{$library}" for library in doc.libraries ];
  ], "\n")
  TikzPicture(data; options=options, preamble=preamble)
end

function Base.show(io::IO, ::MIME"image/svg+xml", doc::TikZ.Document)
  show(io, MIME"image/svg+xml"(), TikzPicture(doc))
end


text:
using Pkg
diffractorspec = PackageSpec(url="https://github.com/JuliaDiff/Diffractor.jl", rev="main")
Pkg.add([diffractorspec])


text:
using MeshCat
using Plots

vis = Visualizer()
open(vis)

mech = Dojo.get_mechanism(:panda,
    timestep=0.05,
    gravity=-0*9.81,
    spring=0.0,
    damper=10.0,
    contact=true,
    limits=true,
    )
# for joint in mech.joints
#     joint.rotational.axis_offset = Quaternion(1,0,0,0.0,true)
# end
function ctrl!(m, k; kp=10.0, kd=5.0)
    nu = input_dimension(m)
    pref = [0,-0.8,0,1.6,0,-2.4,0,0,0]
    y = get_minimal_state(m)
    p = y[1:2:end]
    v = y[2:2:end]
    u = kp * [1,1,1,1,1,1,1,-1,-1] .* (pref - p) - kd * v
    @show norm(p-pref)
    set_input!(m, u*m.timestep)
    return nothing
end

input_dimension(mech.joints[2].rotational)
initialize_panda!(mech,
    joint_angles=[2.0,-0.5,0.3,1.3,0.3,-2.0,0,0,0],
    joint_velocities=[0,0,0,0,0,0,0,0,0.0])

storage = simulate!(mech, 6.0, ctrl!, record=true, opts=SolverOptions(btol=1e-4))
visualize(mech, storage, vis=vis, show_contact=true)

convert_frames_to_video_and_gif("panda_pd_control")

env = get_environment(:panda)
x = srand(minimal_dimension(env.mechanism))
u = srand(input_dimension(env.mechanism))
step(env, x, u)

a = 10
a = 10
a = 10
a = 10
a = 10
a = 10

# mech = get_snake(damper=0.1)
# initialize!(mech, :snake)
mech = get_twister(damper=0.0, joint_type=:Revolute, num_bodies=3)
initialize!(mech, :twister)
# mech = get_snake(damper=0.0, joint_type=:Revolute, num_bodies=3)
# initialize!(mech, :snake)
storage = simulate!(mech, 0.5, record=true, opts=SolverOptions(btol=1e-4))
visualize(mech, storage, vis=vis)

q = rand(QuatRotation).q
Lmat(q) * Tmat() -  Lmat(inv(q))
Tmat() * Lmat(q) -  Lmat(inv(q))
Lmat(q)' -  Lmat(inv(q))

mech = Dojo.get_mechanism(:pendulum,
    timestep=0.01,
    gravity=-0.0,
    spring=0.0,
    damper=0.5)
function ctrl!(m,k)
    set_input!(m, 0.000 * 1 *sones(1))
    return nothing
end
input_dimension(mech)
zero_velocity!(mech)
y = [0.0,10.50]
set_minimal_state!(mech, y)

storage = simulate!(mech, 5.0, Dojo.ctrl!, record=true)
visualize(mech, storage, vis=vis)
get_minimal_state(mech)
storage.ω[1]
plot([-ϕ[1] for ϕ in storage.ω[1]])

xa = szeros(3)
qa = Quaternions.QuaternionF64(1.0, 0.0, 0.0, 0.0, true)
xb = 0*[0.0, 0.4476047997637539, 0.8771773278369399]
qb = mech.bodies[1].state.q2
Q = displacement_jacobian_configuration(:child, mech.joints[1].rotational,
    szeros(3), qa, szeros(3), qb, attjac=true)[2]

VRᵀmat(qa) - VRmat(inv(qa))
LVᵀmat(qa) - LᵀVᵀmat(inv(qa))
VRᵀmat(inv(qb) * qa) * LVᵀmat(inv(qb) * qa) - VRmat(qb) * LᵀVᵀmat(qb) * VRᵀmat(qa) * LVᵀmat(qa)

using Rotations
using FiniteDiff
using Test

mech = Dojo.get_mechanism(:pendulum)
joint0 = mech.joints[1]
rot0 = joint0.rotational
rot0.axis_offset = rand(QuatRotation).q

xa = rand(3)
qa = rand(QuatRotation).q
xb = rand(3)
qb = rand(QuatRotation).q

p0 = rand(3)
J0 = Dojo.impulse_transform_jacobian(:parent, :parent, rot0, xa, qa, xb, qb, p0)
attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
J1 = FiniteDiff.finite_difference_jacobian(
    z -> Dojo.impulse_transform(:parent, rot0, z[1:3], Quaternion(z[4:7]...,true), xb, qb) * p0,
    [xa; Dojo.vector(qa)]
    ) * attjac
@test norm(J0 - J1, Inf) < 1.0e-7

J0 = Dojo.impulse_transform_jacobian(:parent, :child, rot0, xa, qa, xb, qb, p0)
attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
J1 = FiniteDiff.finite_difference_jacobian(
    z -> Dojo.impulse_transform(:parent, rot0, xa, qa, z[1:3], Quaternion(z[4:7]...,true)) * p0,
    [xb; Dojo.vector(qb)]
    ) * attjac
@test norm(J0 - J1, Inf) < 1.0e-7

J0 = Dojo.impulse_transform_jacobian(:child, :parent, rot0, xa, qa, xb, qb, p0)
attjac = cat(I(3),Dojo.LVᵀmat(qa), dims=(1,2))
J1 = FiniteDiff.finite_difference_jacobian(
    z -> Dojo.impulse_transform(:child, rot0, z[1:3], Quaternion(z[4:7]...,true), xb, qb) * p0,
    [xa; Dojo.vector(qa)]
    ) * attjac
@test norm(J0 - J1, Inf) < 1.0e-7

J0 = Dojo.impulse_transform_jacobian(:child, :child, rot0, xa, qa, xb, qb, p0)
attjac = cat(I(3),Dojo.LVᵀmat(qb), dims=(1,2))
J1 = FiniteDiff.finite_difference_jacobian(
    z -> Dojo.impulse_transform(:child, rot0, xa, qa, z[1:3], Quaternion(z[4:7]...,true)) * p0,
    [xb; Dojo.vector(qb)]
    ) * attjac
@test norm(J0 - J1, Inf) < 1.0e-7

rotation_matrix(inv(qb)) * rotation_matrix(qa) - rotation_matrix(inv(qb) * qa)

p0 = rand(3)
q0 = rand(QuatRotation).q
J0 = Dojo.∂rotation_matrix∂q(q0, p0)
J2 = Dojo.∂rotation_matrix∂q(q0, p0, attjac=true)
J1 = FiniteDiff.finite_difference_jacobian(q0 -> Dojo.rotation_matrix(Quaternion(q0...,true))*p0, Dojo.vector(q0))
J3 = J1 * Dojo.LVᵀmat(q0)
norm(J0 - J1, Inf)
norm(J2 - J3, Inf)

J4 = Dojo.∂rotation_matrix_inv∂q(q0, p0)
J6 = Dojo.∂rotation_matrix_inv∂q(q0, p0, attjac=true)
J5 = FiniteDiff.finite_difference_jacobian(q0 -> Dojo.rotation_matrix(inv(Quaternion(q0...,true)))*p0, Dojo.vector(q0))
J7 = J5 * Dojo.LVᵀmat(q0)
norm(J4 - J5, Inf)
norm(J6 - J7, Inf)

qo = Quaternion(1,0,0,0.)

qv = Vmat(inv(qo) * inv(qa) * qb)
Q = Vmat() * Lᵀmat(qo) * Lᵀmat(qa) * LVᵀmat(qb)

qvs = Vmat(inv(qb) * qa * qo)
Qs = Vmat() * Rmat(qo) * Rmat(qa) * Tmat() * LVᵀmat(qb)

x = rand(3)
qra = rand(4)
qrb = rand(4)
qra = Quaternion(qra ./ norm(qra)...)
qrb = Quaternion(qrb ./ norm(qrb)...)
vector_rotate(vector_rotate(x, qra),inv(qrb)) - VRᵀmat(inv(qrb) * qra) * LVᵀmat(inv(qrb) * qra) * x
VRᵀmat(qra) * LVᵀmat(qra) * x - vector_rotate(x, qra)
VLmat(qra) * RᵀVᵀmat(qra) * x - vector_rotate(x, qra)

qra = [1,0,0,1.]
qra = Quaternion(qra ./ norm(qra)..., true)
x = [1,0,0]
vector_rotate(x, qra)
VLmat(qra) * RᵀVᵀmat(qra) * x

VLmat(qb) * RᵀVᵀmat(qb) * x - vector_rotate(x, qb)
VRᵀmat(qb) * LVᵀmat(qb) * x - vector_rotate(x, qb)
VLmat(qb) * RᵀVᵀmat(qb) * x - vector_rotate(x, qb)
LVᵀmat(qrr) * LVᵀmat(qrr)'
VRᵀmat(qb)
Rᵀmat(qb)

x - VRᵀmat(qrr) * LVᵀmat(qrr) * x

VLmat(qb) * RᵀVᵀmat(qb)
qrr = rand(4)
qrr = [1,0,0,0.]
qrr = Quaternion(qrr ./ norm(qrr))
VLmat(qrr)
RᵀVᵀmat(qrr)

vector_rotate([1,0,0], qb)
Q * [1,0,0]

vis = Visualizer()
open(vis)

mech = Dojo.get_mechanism(:panda,
    timestep=0.01,
    gravity=[0,0,-9.0],
    spring=0.0,
    damper=0.1)
# mech.joints[1].rotational.axis_offset = Quaternion(1,0,0,0.0,true)
function ctrl!(m,k)
    nu = minimal_dimension(mech)
    set_input!(m, 0.000 * 1 *sones(nu))
    return nothing
end
mech.joints
zero_velocity!(mech)
# ny = minimal_dimension(mech)
# y = [0.0,10]
# set_minimal_state!(mech, y)

storage = simulate!(mech, 3.0, Dojo.ctrl!, record=true)
visualize(mech, storage, vis=vis)
get_minimal_state(mech)
storage.ω[1]
plot([-ϕ[3] for ϕ in storage.ω[1]])

joint = mech.joints[1]
body = mech.bodies[1]
pbody = get_body(mech, joint.parent_id)
cbody = get_body(mech, joint.child_id)
cbody
damper_impulses(mech, joint, body)
damper_impulses(:child, joint.rotational,
    get_body(mech, joint.parent_id),
    get_body(mech, joint.child_id),
    mech.timestep,
    unitary=false)
damper_force(:child, joint.rotational, current_configuration(pbody.state)[2], pbody.state.ϕsol[2],
    current_configuration(cbody.state)[2], cbody.state.ϕsol[2], 0.01, unitary=false)

function min_vel(qb)
    xa = SVector{3}([0.0, 0.0, 0.0])
    va = SVector{3}([0.0, 0.0, 0.0])
    qa = Quaternions.QuaternionF64(1.0, 0.0, 0.0, 0.0, true)
    ωa = SVector{3}([0.0, 0.0, 0.0])
    xb = SVector{3}([0.0, 0.0, 0.0])
    vb = SVector{3}([0.0, 0.0, 0.0])
    ωb = SVector{3}([0.0, 0.0, 1.0])
    return minimal_velocities(joint.rotational, xa, va, qa, ωa, xb, vb, qb, ωb, 0.05)
end

plot([min_vel(qb) for qb in storage.q[1]])

function get_damping(θ, ω)
    y = [θ, ω]
    set_minimal_state!(mech, y)
    imp = damper_impulses(mech, joint, body)
    # imp = damper_force(:child, joint.rotational, current_configuration(pbody.state)[2], pbody.state.ϕsol[2],
    #     current_configuration(cbody.state)[2], cbody.state.ϕsol[2], 0.01, unitary=false)
    return imp[6]
end
ys = [[2π*x,0] for x in 0:0.005:2]
zs = [minimal_to_maximal(mech, y) for y in ys]
storage = generate_storage(mech, zs)
visualize(mech, storage, vis=vis)

plot([get_damping(0.0, x) for x in 0.0:0.01:1.0])
plot([get_damping(2π*x, 1.0) for x in 0.0:0.01:2])

minimal_dimension(mech)

mech = get_npendulum(num_bodies=1, damper=-0.1)
initialize!(mech, :npendulum)
y = [0.5,4]
set_minimal_state!(mech, y)

storage = simulate!(mech, 3.3, record=true)
visualize(mech, storage, vis=vis)
get_minimal_state(mech)

mech = get_quadruped(damper=0.1, timestep=0.02)
initialize!(mech, :quadruped)

storage = simulate!(mech, 3.3, record=true, verbose=true)
Main.@profiler storage = simulate!(mech, 0.08, record=true, verbose=true)
visualize(mech, storage, vis=vis)


text:
include("../../neat_test.jl")

const DEFAULT_DATASET_PATH = "../../../../src/datasets"

println("--> Dataset 2")
config = "test_dataset2_config.txt"
#config = ARGS[1]
#generations = parse(Int, ARGS[2])
generations = 5000
st2 = run_NEAT_epochs(config, generations, false)
println()
print_best_solution_perfomance(st2)
println()

#print("* Calculating information plane points")
#calculate_information_plane_points(st1)

plot_neat_fitness(st2, imagename="fitness-dt2.png")
plot_output_distribution(st2, imagename="distribution-dt2.png")
plot_complexity(st2, imagename="solution-complexity-dt2.png")
plot_total_accuracy(st2, imagename="total-accuracy-dt2.png")
plot_output_distribution_test(st2, imagename="output-test-distribution-best-NEAT-dt2.png")

print(st2.population.best_fitness[end])

#=
st2 = get_simulation_state("checkpoint-dt2-5000g-relatorio.jld2")

plot_neat_fitness(st2, imagename="fitness-dt2.png")
plot_neat_IP(st2, imagename="IP-dt2.png", nbins=10)
plot_neat_species(st2, imagename="species-dt2.png")
plot_output_distribution(st2, imagename="distribution-dt2.png")
plot_complexity(st2, imagename="solution-complexity-dt2.png")
plot_accuracy(st2, imagename="accuracy-dt2.png")
plot_total_accuracy(st2, imagename="total-accuracy-dt2.png")
plot_output_distribution_test(st2, imagename="output-test-distribution-best-NEAT-dt2.png")
plot_confusion_matrix(st2; imagename="confusion-matrix-testset-dt2.png", set=:test)
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
"""
    Imputor

An imputor stores information about imputing values in `AbstractArray`s and `Tables.table`s.
New imputation methods are expected to subtype `Imputor` and, at minimum,
implement the `_impute!(data::AbstractArrays, imp::<MyImputor>)` method.

While fallback `impute` and `impute!` methods are provided to extend your `_impute!` methods to
n-dimensional arrays and tables, you can always override these methods to change the
behaviour as necessary.
"""
abstract type Imputor end

#=
These default methods are required because @auto_hash_equals doesn't
play nice with Base.@kwdef
=#
function Base.hash(imp::T, h::UInt) where T <: Imputor
    h = hash(Symbol(T), h)

    for f in fieldnames(T)
        h = hash(getfield(imp, f), h)
    end

    return h
end

function Base.:(==)(a::T, b::T) where T <: Imputor
    result = true

    for f in fieldnames(T)
        if !isequal(getfield(a, f), getfield(b, f))
            result = false
            break
        end
    end

    return result
end

"""
    impute(data::T, imp; kwargs...) -> T

Returns a new copy of the `data` with the missing data imputed by the imputor `imp`.
For matrices and tables, data is imputed one variable/column at a time.
If this is not the desired behaviour then you should overload this method or specify a different `dims` value.

# Arguments
* `data`: the data to be impute
* `imp::Imputor`: the Imputor method to use

# Returns
* the input `data` with values imputed

# Example
```jldoctest
julia> using Impute: Interpolate, impute

julia> v = [1.0, 2.0, missing, missing, 5.0]
5-element Vector{Union{Missing, Float64}}:
 1.0
 2.0
  missing
  missing
 5.0

julia> impute(v, Interpolate())
5-element Vector{Union{Missing, Float64}}:
 1.0
 2.0
 3.0
 4.0
 5.0
```
"""
function impute(data, imp::Imputor; kwargs...)
    # NOTE: We don't use a return type declaration here because `trycopy` isn't guaranteed
    # to return the same type passed in. For example, subarrays and subdataframes will
    # return a regular array or dataframe.
    return impute!(trycopy(data), imp; kwargs...)
end

"""
    impute!(data::A, imp; dims=:, kwargs...) -> A

Impute the `missing` values in the array `data` using the imputor `imp`.
Optionally, you can specify the dimension to impute along.

# Arguments
* `data::AbstractArray{Union{T, Missing}}`: the data to be impute along dimensions `dims`
* `imp::Imputor`: the Imputor method to use

# Keyword Arguments
* `dims=:`: The dimension to impute along. `:rows` and `:cols` are also supported for matrices.

# Returns
* `AbstractArray{Union{T, Missing}}`: the input `data` with values imputed

# NOTES
1. Matrices have a deprecated `dims=2` special case as `dims=:` is a breaking change
2. Mutation isn't guaranteed for all array types, hence we return the result
3. `eachslice` is used internally which requires Julia 1.1

# Example
```jldoctest
julia> using Impute: Interpolate, impute!

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute!(M, Interpolate(); dims=1)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  3.0  4.0  5.0
 1.1  2.2  3.3  4.4  5.5

julia> M
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  3.0  4.0  5.0
 1.1  2.2  3.3  4.4  5.5
```
"""
function impute!(
    data::A, imp::Imputor; dims=:, kwargs...
)::A where A <: AbstractArray{Union{T, Missing}} where T
    dims === Colon() && return _impute!(data, imp; kwargs...)

    for x in eachslice(data; dims=dims)
        _impute!(x, imp; kwargs...)
    end

    return data
end

function impute!(
    data::M, imp::Imputor; dims=nothing, kwargs...
)::M where M <: AbstractMatrix{Union{T, Missing}} where T
    dims === Colon() && return _impute!(data, imp; kwargs...)
    # We're calling our `dim` function to throw a depwarn if `dims === nothing`
    d = dim(data, dims)

    for x in eachslice(data; dims=d)
        _impute!(x, imp; kwargs...)
    end

    return data
end

impute!(data::AbstractMatrix{Missing}, imp::Imputor; kwargs...) = data

"""
    impute!(data::T, imp; kwargs...) -> T where T <: AbstractVector{<:NamedTuple}

Special case rowtables which are arrays, but we want to fallback to the tables method.
"""
function impute!(data::T, imp::Imputor)::T where T <: AbstractVector{<:NamedTuple}
    return materializer(data)(impute!(Tables.columns(data), imp))
end

"""
    impute!(data::AbstractArray, imp) -> data

Just returns the `data` when the array doesn't contain `missing`s
"""
impute!(data::AbstractArray, imp::Imputor; kwargs...) = disallowmissing(data)

"""
    impute!(data::AbstractArray{Missing}, imp) -> data

Just return the `data` when the array only contains `missing`s
"""
impute!(data::AbstractArray{Missing}, imp::Imputor; kwargs...) = data

"""
    impute!(table, imp; cols=nothing) -> table

Imputes the data in a table by imputing the values 1 column at a time;
if this is not the desired behaviour custom imputor methods should overload this method.

# Arguments
* `imp::Imputor`: the Imputor method to use
* `table`: the data to impute

# Keyword Arguments
* `cols`: The columns to impute along (default is to impute all columns)

# Returns
* the input `data` with values imputed

# Example
```jldoctest
julia> using DataFrames; using Impute: Interpolate, impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
 Row │ a          b
     │ Float64?   Float64?
─────┼──────────────────────
   1 │       1.0        1.1
   2 │       2.0        2.2
   3 │ missing          3.3
   4 │ missing    missing
   5 │       5.0        5.5

julia> impute(df, Interpolate())
5×2 DataFrame
 Row │ a         b
     │ Float64?  Float64?
─────┼────────────────────
   1 │      1.0       1.1
   2 │      2.0       2.2
   3 │      3.0       3.3
   4 │      4.0       4.4
   5 │      5.0       5.5
```
"""
function impute!(table::T, imp::Imputor; cols=nothing)::T where T
    # TODO: We could probably handle iterators of tables here
    istable(table) || throw(MethodError(impute!, (table, imp)))

    # Extract a columns iterator that we should be able to use to mutate the data.
    # NOTE: Mutation is not guaranteed for all table types, but it avoid copying the data
    columntable = Tables.columns(table)

    cnames = cols === nothing ? propertynames(columntable) : cols
    for cname in cnames
        impute!(getproperty(columntable, cname), imp)
    end

    return table
end

files = [
    "interp.jl",
    "knn.jl",
    "locf.jl",
    "nocb.jl",
    "replace.jl",
    "srs.jl",
    "substitute.jl",
    "svd.jl",
]

for file in files
    include(joinpath("imputors", file))
end


text:
lineprinter = Sys.iswindows() ? "LPT3" : "/dev/lp0"
lp = open(lineprinter, "w")
write(lp, "Hello world")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PrivateMultiplicativeWeights
using CSV
import Plots

weights = CSV.read("our_synthetic_weights.csv");
d = Vector(weights[:weight]);
# Magic number is number of samples (from python notebook)
h = Histogram(d, 38462)

epsilons = reverse([10.0^-n for n in range(1, stop=6)])
mses = []
max_errors = []

for ε in epsilons
    println("Running MWEM with ε=$ε")
    mw = mwem(SeriesRangeQueries(length(d)), h, MWParameters(epsilon=ε, iterations=500));
    mse = mean_squared_error(mw)
    println("Mean squared error = ", mse)
    append!(mses, mse)
end

Plots.plot(epsilons, mses, xlabel = "Epsilon", ylabel="MSE", title="Utility vs Privacy", xaxis=:log)

max_stars_repo_path: Scripts/Round 3/Baysically Measure Zero/mwem.jl