text:
using Flux, CUDA, Test, FluxExtra, FluxExtra.Normalizations

function test_training(model,x,y)
    opt = ADAM()
    loss = Flux.Losses.mse
    losses = Vector{Float32}(undef,2)
    for i = 1:2
        local loss_val
        ps = Flux.Params(Flux.params(model))
        gs = gradient(ps) do
            predicted = model(x)
            loss_val = loss(predicted,y)
        end
        losses[i] = loss_val
        Flux.Optimise.update!(opt,ps,gs)
    end
    if losses[1]==losses[2]
        error("Parameters not updating.")
    end
    return true
end

function test_model(model,x,y)
    test_exp = quote
        @test begin
            @inferred $model($x)
            true
        end
        @test test_training($model,$x,$y)
        @test begin
            @inferred gpu($model)(gpu($x))
            true
        end
        @test test_training(gpu($model),gpu($x),gpu($y))
    end
    return test_exp
end

#---Tests begin-------------------------------------------------------------

@testset verbose=true "Conv as input" begin
    test_layer = Conv((3, 3), 1=>2,pad=SamePad())
    test_layer2 = Conv((3, 3), 1=>2,pad=SamePad())

    @testset "Join" begin
        @testset "Join(1)" begin
            x = ones(Float32,6,6,1,1)
            y = ones(Float32,12,6,2,1)
            model = Chain(Parallel(tuple,(test_layer,test_layer2)),Join(1))
            eval(test_model(model,x,y))
        end

        @testset "Join(2)" begin
            x = ones(Float32,6,6,1,1)
            y = ones(Float32,6,12,2,1)
            model = Chain(Parallel(tuple,(test_layer,test_layer2)),Join(2))
            eval(test_model(model,x,y))
        end

        @testset "Join(dum = 3)" begin
            x = ones(Float32,6,6,1,1)
            y = ones(Float32,6,6,4,1)
            model = Chain(Parallel(tuple,(test_layer,test_layer2)),Join(dim = 3))
            eval(test_model(model,x,y))
        end

        @testset "Join errors" begin
            @test try
                    Join(4)
                catch e
                    if !(e isa DimensionMismatch)
                        @error "Wrong error returned."
                        return e
                    else
                        true
                    end
                end

            @test try
                    Join(dim = 4)
                catch e
                    if !(e isa DimensionMismatch)
                        @error "Wrong error returned."
                        return e
                    else
                        true
                    end
                end
        end

        @testset "Printing" begin
            @test begin
                Base.show(IOBuffer(),Join(1))
                true
            end
        end
    end

    @testset "Split" begin
        @testset "Split(_,1)" begin
            x = ones(Float32,6,6,1,1)
            y = ones(Float32,6,6,2,1)
            model = Chain(Split(2,1),Parallel(tuple,(test_layer,test_layer)),Join(1))
            eval(test_model(model,x,y))
        end

        @testset "Split(_,2)" begin
            x = ones(Float32,6,6,1,1)
            y = ones(Float32,6,6,2,1)
            model = Chain(Split(2,2),Parallel(tuple,(test_layer,test_layer)),Join(2))
            eval(test_model(model,x,y))
        end

        @testset "Split(_,dim = 3)" begin
            x = ones(Float32,6,6,2,1)
            y = ones(Float32,6,6,4,1)
            model = Chain(Split(2,dim = 3),Parallel(tuple,(test_layer,test_layer)),Join(3))
            eval(test_model(model,x,y))
        end

        @testset "Split errors" begin
            @test try
                    Split(2,4)
                catch e
                    if !(e isa DimensionMismatch)
                        @error "Wrong error returned."
                        return e
                    else
                        true
                    end
                end

            @test try
                    Split(2,dim = 4)
                catch e
                    if !(e isa DimensionMismatch)
                        @error "Wrong error returned."
                        return e
                    else
                        true
                    end
                end

            @test try
                    Split(1,1)
                catch e
                    if !(e isa DomainError)
                        @error "Wrong error returned."
                        return e
                    else
                        true
                    end
                end

            @test try
                    Split(1,dim = 1)
                catch e
                    if !(e isa DomainError)
                        @error "Wrong error returned."
                        return e
                    else
                        true
                    end
                end
        end

        @testset "Printing" begin
            @test begin
                Base.show(IOBuffer(),Split(2,3))
                true
            end
        end
    end

    @testset "Addition" begin
        x = ones(Float32,6,6,1,1)
        y = ones(Float32,6,6,1,1)
        model = Chain(Parallel(tuple,(test_layer,test_layer)),Addition())
        eval(test_model(model,x,y))
    end

    @testset "Activation" begin
        x = ones(Float32,4,4,1,1)
        y = ones(Float32,4,4,2,1)
        model = Chain(test_layer,Activation(tanh))
        eval(test_model(model,x,y))

        @testset "Printing" begin
            @test begin
                Base.show(IOBuffer(),Activation(tanh))
                true
            end
        end
    end

    @testset "Flatten" begin
        x = ones(Float32,4,4,1,1)
        y = ones(Float32,32,1)
        model = Chain(test_layer,Flatten())
        eval(test_model(model,x,y))
    end

    @testset "Identity" begin
        x = ones(Float32,4,4,1,1)
        y = ones(Float32,4,4,3,1)
        model = Chain(Parallel(tuple,(test_layer,Identity())),Join(3))
        eval(test_model(model,x,y))
    end
end

@testset verbose=true "Dense as input" begin
    test_layer = Dense(2,3)
    test_layer2 = Dense(2,3)

    @testset "Join" begin
        @testset "Join(1)" begin
            x = ones(Float32,2,1)
            y = ones(Float32,6,1)
            model = Chain(Parallel(tuple,(test_layer,test_layer2)),Join(1))
            eval(test_model(model,x,y))
        end
    end

    @testset "Split" begin
        @testset "Split(_,1)" begin
            x = ones(Float32,4,1)
            y = ones(Float32,6,1)
            model = Chain(Split(2,1),Parallel(tuple,(test_layer,test_layer)),Join(1))
            eval(test_model(model,x,y))
        end

        @testset "Split errors" begin
            x = ones(Float32,2,1)
            layer = Split(2,2)
            @test try
                    layer(x)
                catch e
                    if !(e isa DimensionMismatch)
                        error("Wrong error returned.")
                        return e
                    else
                        true
                    end
                end
        end
    end

    @testset "Addition" begin
        x = ones(Float32,4,1)
        y = ones(Float32,3,1)
        model = Chain(Split(2,1),Parallel(tuple,(test_layer,test_layer)),Addition())
        eval(test_model(model,x,y))
    end

    @testset "Activation" begin
        x = ones(Float32,2,1)
        y = ones(Float32,3,1)
        model = Chain(test_layer,Activation(tanh))
        eval(test_model(model,x,y))
    end

    @testset "Identity" begin
        x = ones(Float32,2,1)
        y = ones(Float32,5,1)
        model = Chain(Parallel(tuple,(test_layer,Identity())),Join(1))
        eval(test_model(model,x,y))
    end
end

@testset "Normalizations" begin
    data = repeat([rand(Float32,5,5,3)],2)
    @test begin
        min_vals,max_vals = norm_01!(data)
        map(x -> norm_01!(x,min_vals,max_vals), data)
        min_vals,max_vals = norm_negpos1!(data)
        map(x -> norm_negpos1!(x,min_vals,max_vals), data)
        mean_vals = norm_zerocenter!(data)
        mean_vals,std_vals = norm_zscore!(data)
        true
    end
end


text:
import Base.length
using LinearAlgebra
using StatsBase
using Plots
using Random
using DelimitedFiles

abstract type KernelSpace end

struct RBFSpace <: KernelSpace
    x::Matrix{Float64}
    sigma::Float64
end

function length(K::RBFSpace)::Int64
    return size(K.x, 2)
end

function(K::RBFSpace)(i::Int64)::Vector{Float64}
    return K.x[:, i]
end

function (K::RBFSpace)(p::Vector{Float64}, q::Vector{Float64})::Float64
    l2 = sum((p-q).^2)
    return exp(-l2/(2*K.sigma^2))
end

function (K::RBFSpace)(i::Int64, j::Int64)::Float64
    return K(K(i), K(j))
end

function(K::RBFSpace)(i::Int64, p::Vector{Float64})::Float64
    return K(K(i), p)
end

function(K::RBFSpace)(p::Vector{Float64}, i::Int64)::Float64
    return K(p, K(i))
end

function GridSpace(xs, sigma::Float64)::RBFSpace
    return RBFSpace(hcat([[p...] for p in Iterators.product(xs...)]...), sigma)
end

struct Sample
    K::KernelSpace
    S::Vector{Int64}
    M::Matrix{Float64} # Inverse of K[S, S]
    U::Matrix{Float64} # K[S, :]
    L::Vector{Float64} # K[., .]-K[., :] M K[:, .]
end

function (s::Sample)(p::Vector{Float64})
    u = s.K.(s.S, [p])
    return s.K(p, p)-u'*s.M*u
end

function value(s::Sample)::Float64
    return det(s.K.(s.S, s.S'))
end

@userplot PointPlot
@recipe function f(inp::PointPlot)
    s = inp.args[1]
    grid := false
    space = hcat(s.K.(1:length(s.K))...)
    showaxis := false
    label := ""
    if size(space, 1)>=1
        xlims := [min(space[1, :]...), max(space[1, :]...)]
    end
    if size(space, 1)>=2
        ylims := [min(space[2, :]...), max(space[2, :]...)]
    end
    if size(space, 1)>=3
        zlims := [min(space[3, :]...), max(space[3, :]...)]
    end
    @series begin
        seriestype := :scatter
        data = hcat(s.K.(s.S)...)
        Tuple(data[i, :] for i=1:size(data, 1))
    end
end

@userplot HeatPlot
@recipe function f(inp::HeatPlot)
    s = inp.args[1]
    space = hcat(s.K.(1:length(s.K))...)
    xrange = range(min(space[1, :]...), max(space[1, :]...), length=100)
    yrange = range(min(space[1, :]...), max(space[1, :]...), length=100)
    @series begin
        seriestype := :heatmap
        xrange, yrange, (x, y)->s([x, y])
    end
end

function length(s::Sample)
    return length(s.S)
end

function initial(K::KernelSpace, S::Vector{Int64}=Int64[])::Sample
    U = K.(S, (1:length(K))')
    M = inv(U[:, S])
    L = K.(1:length(K), 1:length(K))-sum((M*U).*U, dims=1)[1, :]
    return Sample(K, S, M, U, L)
end

function sane(s::Sample)::Sample
    return initial(s.K, s.S)
end

function sanity(s::Sample)::Float64
    t = sane(s)
    eM = maximum(abs.(s.M-t.M))
    eU = maximum(abs.(s.U-t.U))
    eL = maximum(abs.(s.L-t.L))
    return max(eM, eU, eL)
end

function added(s::Sample, ind::Int64)::Sample
    k = length(s)
    S = [s.S; ind]
    U = [s.U; s.K.(ind, 1:length(s.K))']
    M1 = [s.M zeros(k, 1); zeros(1, k) 0]
    v = [-s.M*s.U[:, ind]; 1.]
    M = M1 + (v.*v')/s.L[ind]
    L = s.L-(((v'*U).^2)/s.L[ind])[1, :]
    return Sample(s.K, S, M, U, L)
end

function removed(s::Sample, i::Int64)::Sample
    S = [s.S[1:i-1]; s.S[i+1:end]]
    U = [s.U[1:i-1, :]; s.U[i+1:end, :]]
    v = s.M[:, i]
    M1 = s.M-(v.*v')./v[i]
    M = [M1[1:i-1, 1:i-1] M1[1:i-1, i+1:end]; M1[i+1:end, 1:i-1] M1[i+1:end, i+1:end]]
    L = s.L+(((v'*s.U).^2)./v[i])[1, :]
    return Sample(s.K, S, M, U, L)
end

function greedy_next(s::Sample)::Int64
    return argmax(s.L)
end

function random_next(s::Sample)::Int64
    return sample(1:length(s.K), Weights(max.(s.L, 1e-15)))
end

function output_points(filename::String, s::Sample)
    open(filename, "w") do io
        writedlm(io, ["x" "y"])
        writedlm(io, s.K.(s.S))
    end
end

function output_heatmap(filename::String, s::Sample, samples::Int64=21)
    space = hcat(s.K.(1:length(s.K))...)
    xrange = range(min(space[1, :]...), max(space[1, :]...), length=samples)
    yrange = range(min(space[1, :]...), max(space[1, :]...), length=samples)
    open(filename, "w") do io
        writedlm(io, ["x" "y" "z"])
        for x=xrange
            for y=yrange
                writedlm(io, [x y s([x, y])])
            end
            write(io, "\n")
        end
    end
end

function greedy(K::KernelSpace, k::Int64)::Sample
    s = initial(K)
    while length(s)<k
        s = added(s, greedy_next(s))
    end
    return s
end

function convex(K::KernelSpace, k::Int64)::Sample
    base = initial(K)
    while length(base)<k
        z = Float64[1. for i=1:length(K)]
        for i=base.S
            z[i] = 0.
        end
        s = initial(K, base.S)
        while length(s)<k
            s = added(s, argmax(s.L.*z))
        end
        for iter=1:(length(K)*log(k+1))
            z ./= sum(z)/length(z)
            s = removed(s, rand((length(base)+1):length(s)))
            probs = max.(s.L.*z, 1e-15)
            probs ./= sum(probs)
            ind = sample(1:length(K), Weights(probs))
            s = added(s, ind)
            if iter>2*k*log(k+1)
                z .+= probs
            end
        end
        ind = argmax(z)
        z[ind] = 0.
        base = added(base, ind)
    end
    return base
end

function run3()
  for iters=1:100
    N = 200
    D = 2
    k = 3
    sigma = 1.
    K = RBFSpace(rand(D, N), sigma)

    G = greedy(K, k)
    println(value(G))
    C = convex(K, k)
    println(value(C))
    println()
  end
end

function run20()
  for iters=1:100
    N = 200
    D = 2
    k = 20
    sigma = .2
    K = RBFSpace(rand(D, N), sigma)

    G = greedy(K, k)
    println(value(G))
    C = convex(K, k)
    println(value(C))
    println()
  end
end

run20()

#Random.seed!(0)

#pointplot(C)

#pointplot(convex(K, 10))
# println(greedy(K, 10))

# s = initial(K)
# for i=1:150
#     while true
#         j = rand(1:length(K))
#         if !(j in s.S)
#             global s=added(t, j)
#             break
#         end
#     end
# end
#
# for t=1:1000
#     global s = removed(s, rand(1:length(s)))
#     global s = added(s, random_next(s))
#     output_points("data/points_$(t).dat", s)
#     output_heatmap("data/heat_$(t).dat", s, 51)
# end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 2db004b0-7c45-11eb-2c3a-3f0e111cf63a
using LinearAlgebra, PlutoUI

# ╔═╡ 55f8789f-5361-421f-9f37-98c154fd40cf
begin
	using Dates
	Dates.today()
end

# ╔═╡ 5ea33223-4dad-4869-8724-93b7cb91771f
TableOfContents(title="📚 Table of Contents", aside=true)

# ╔═╡ 8ccdfcf9-d5d1-4b69-a854-15f2715bfafa
md"""
# Julia is Open

`Julia` is an open-source project, source being entirely hosted on [github](https://github.com/JuliaLang/julia).

The code consists of (these are old numbers, actual numbers may differ):

- 29K lines of `C/C++`
- 6K lines of `scheme`
- 68K lines of `julia`

Julia uses [LLVM](http://llvm.org) which itself has 680K lines of code. Therefore, `Julia` is very compact, compared to other languages, like LLVM's `C` compiler `clang` (513K lines of code) or `gcc`  (3,530K lines). This makes it easy to read the actuall code and get full information, in spite the fact that some parts of the documentation are insufficient. `Julia`'s "navigating" system,
consisting of commands `varinfo()`, `methods()` and `@which`, makes this even easier.

Further, the `Base` (core) of `Julia` is kept small, and the rest of the functionality is obtained through packages.
Since packages are written in `Julia`, they are navigated on the same way.

In this notebook, we demonstrate how to get help and navigate the source code.

__Prerequisites__

Basic knowledge of programming in any language.

Read [Methods](https://docs.julialang.org/en/v1/manual/methods/#Methods) section of the `Julia` manual. (5 min)

__Competences__

The reader should be able to read the code and be able to find and understand calling sequences and outputs of any function.

__Credits__

Some examples are taken from [The Julia Manual](https://docs.julialang.org/en/v1/).

"""

# ╔═╡ 4440847c-db90-44a5-be48-08c5794ee5b7
md"""
# Operators `+`, `*` and `⋅`

Consider operators `+`, `*` and `⋅`, the first two of them seem rather basic in any language. The `⋅` symbol is typed as LaTeX command `\cdot` + `Tab`.

`?+` gives some information, which is vary sparse. We would expect more details, and we also suspect that `+` can be used in more ways that just hose two.

`?*` explaind more instances where `*` can be used, but the text itself is vague and not sufficient.

`?⋅` appears to be what we expect fro the dot product off two vectors.
"""

# ╔═╡ 7a55f84e-1439-45ec-b428-5976cac084cd
#?+

# ╔═╡ 7514d88c-d0fa-45e5-947d-10fe0ce47f7d
md"""
# methods()

`Julia` functions have a feature called _multiple dispatch_, which means that the method depends on the name __AND__ the input.
Full range of existing methods for certain function name is given by the `methods()` command.
> Running `methods(+)` sheds a completely differfent light on `+`.

The great Julia feature is that the links to the source code where the respective version of the function is defined, are readily provided.
"""

# ╔═╡ 611770e1-4b5b-47a6-be3a-e0806990796d
md"""
## The `"+"` operator
"""

# ╔═╡ 9c208420-7cdf-11eb-2b99-e32e66b46e85
# methods(+)

# ╔═╡ bfb87460-7cdf-11eb-38ca-0961560c093b
# methodswith(Tridiagonal)

# ╔═╡ cb3e314e-68e0-4fc0-8ece-c5c8bc279f55
md"""
Following the third link, we get the following code snippet:
```
+(x::Bool) = int(x)
-(x::Bool) = -int(x)
+(x::Bool, y::Bool) = int(x) + int(y)
-(x::Bool, y::Bool) = int(x) - int(y)
*(x::Bool, y::Bool) = x & y
```

Therefore:
"""

# ╔═╡ 0f3c49cd-00c7-4654-9dc9-674a89e19a30
+(true), +(false),-(true),-(false)

# ╔═╡ 595878a2-30ab-4a19-a5ad-e621dd37f6c7
xb, yb = BitArray([0,1,0,1,0,1]), BitArray([0,1,1,1,1,0])

# ╔═╡ 2a255a09-3c3f-4fcb-99f2-b31adfb536b5
md"""
The above command is equivalent to
```
xb = bitpack([0,1,0,1,0,1]); yb = bitpack([0,1,1,1,1,0])
```
except that only the last result would be displayed
"""

# ╔═╡ 4f6d3190-45c9-433d-b7f9-c42ed01dc04b
+xb,-(xb)

# ╔═╡ 60a3f371-d9b2-47a8-be8f-fab7cfb3fb48
xb+yb, +(xb,yb)

# ╔═╡ 74d98ae5-0641-4eea-9f8c-8e663a405f32
md"""
## Manipulating dates

We see that one of the `+` methods is adding days to time:

```
 +(x::Date,y::Base.Dates.Day) at dates/arithmetic.jl:77
```
Therefore, the 135-th day from today is:
"""

# ╔═╡ 48783355-31c8-4159-8dd5-8bea047c742c
dd=Dates.today()+Dates.Day(135)

# ╔═╡ c9a0df42-1368-4e5d-a21d-d924f7ffbb75
typeof(dd)

# ╔═╡ 79ca20e7-446e-4264-ae24-c46fbb4b15c1
md"""
More information about the two types can be obtained by `methods(Dates.Date)` and `methods(Dates.Day)`, respectively.
"""

# ╔═╡ 7e6013ea-1e66-4a81-86a8-a69e22cde29a
md"""
## Adding tridiagonal matrices

In the above output of `methods(+)`,  we see that we can add tridiagonal matrices:
```
+(A::Tridiagonal{T}, B::Tridiagonal{T}) at linalg/tridiag.jl:624
```
Following the link, we see that the method separately adds lower, main and upper diagonals, denoted by `dl`, `d` and `du`, respectively:
```
624: +(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl+B.dl, A.d+B.d, A.du+B.du)
```
Let us see how exactly is the type `Tridiagonal` defined:
"""

# ╔═╡ 8591bb73-12e5-474e-b24b-7acfab2e3c5f
# methods(Tridiagonal)

# ╔═╡ 8c0b2316-26fa-4748-b047-3d2c6bf7acc0
md"""
This output seems confusing, but from the second line we conclude that we can define three diagonals, lower, main and upper diagonal, denoted as above. We also know that that the lower and upper diagonals are of size $n-1$. Let us try it out:
"""

# ╔═╡ 4163dd6e-7b5c-42f4-a9f4-c1e1a7faecba
T₁ = Tridiagonal(rand(6),rand(7),rand(6))

# ╔═╡ 89bf9619-7ee7-4ccd-99c0-48583e93bbe1
Matrix(T₁)

# ╔═╡ 684b3f4a-f83b-440d-9b13-d5dae159e77a
T₂ = Tridiagonal(rand(-5.0:5,6),randn(7),rand(-9.0:0,6))

# ╔═╡ 243854d4-495a-43de-88cb-1442275c8a8a
T₃=T₁+T₂

# ╔═╡ d79e7c32-ad00-4373-ae0e-04eeb8f6e333
md"""
This worked as expected, the result is again a `Tridiagonal`. We can access each diagonal by:
"""

# ╔═╡ 22af8627-6713-40a5-b716-81a55a6b7aad
# Output is in the console
println(T₃.dl, T₃.d, T₃.du)

# ╔═╡ 98142e4f-c703-4ea9-b893-5e2c831afa42
md"""
## `@which`

Let us take a closer look at what happens. The `@which` command gives the link to the part of the code which is actually invoked. The argument should be only function, without assignment.
"""

# ╔═╡ dac4c688-cd5f-4ef5-8ba6-3709f945d3f0
@which T₁ + T₂

# ╔═╡ 60469aa5-8d5b-4996-a02d-100a65e0915f
@which Tridiagonal(rand(6),rand(7),rand(6))

# ╔═╡ c4bbf3a7-9497-474a-8174-db9530853486
md"""
In the code, we see that there is a type definition in the `struct` block:
```
## Tridiagonal matrices ##
struct Tridiagonal{T} <: AbstractMatrix{T}
    dl::Vector{T}    # sub-diagonal
    d::Vector{T}     # diagonal
    du::Vector{T}    # sup-diagonal
    du2::Vector{T}   # supsup-diagonal for pivoting
end
```
The `Tridiagonal` structure (or type) consists of __four__ vectors.
In our case, we actually called the function `Tridiagonal()` with __three__ vector arguments. The function creates the type of the same name, setting the fourth reqired vector `du2` to `zeros(T,n-2)`.

The next function with the same name is invoked when the input vectors have different types, in which case the types arer promoted to a most general one, if possible.
"""

# ╔═╡ 88ee9b5a-a1f1-49a1-b6cf-64ebac21634a
md"""
# `size()` and `Matrix()`

For each matrix type we need to define the function which returns the size of a matrix, and the function which converts the matrix of a given type to a full array. These function are listed after the second `Tridiagonal()` function.
"""

# ╔═╡ 5e2c532c-49f3-411c-b110-70c54263f007
T₄ = Tridiagonal([1,2,3],[2,3,4,5],[-1,1,2])

# ╔═╡ 7f026e3e-bf1c-49e5-b815-a9370eb5e3f5
size(T₄)

# ╔═╡ 20f20bfa-9aff-46a0-b358-5fe42640d21a
Matrix(T₄)

# ╔═╡ bd0a3a98-4311-41b0-9209-5e09c4ecd9ad
md"""
## `sizeof()`

Of course, using special types can leasd to much more efficient programs. For example, for `Tridiagonal` type, onlt four diagonals are stored, in comparison to storing full matrix when $n^2$ elements are stored. The storage used is obtained by the `sizeof()` function.
"""

# ╔═╡ e0a4eeaf-c259-453e-82aa-6813538ab536
T₁

# ╔═╡ 41daa15d-9b1f-47df-ad24-b124161ce25a
T₁f=Matrix(T₁)

# ╔═╡ cbb653fc-6ec1-408b-b683-7d6821724051
sizeof(T₁f)  #   392 =  7 * 7 * 8 bytes

# ╔═╡ e6a068a5-cb01-40f8-8f0c-19e84f8adbc0
# This is not yet implemented for Tridiagonal - only the storage
# required for 4 vector variables' names is displayed
sizeof(T₁)

# ╔═╡ 7622b55d-49ba-4f21-8d85-4bfb95f0fe37
md"""
# `struct` is immutable

This means that we can change individual elements of defined parts, but not the parts as a whole (an alternative is to use the `type` construtor). For example:
"""

# ╔═╡ 3e9cf5d9-f97f-468d-ba27-b6dc217418d8
begin
	@show T₄
	T₄.d[2]=123
	@show T₄
	# T₄.dl = [-1, -1 ,1]
end

# ╔═╡ ec7be6ff-07ff-4b6d-bf00-2c48a1535582
md"""
# `methodswith()`

This is the reverse of `methods()` - which methods exist for the given type. For example, what can we do with `Tridiagonal` matrices, or with `Dates.Day`:
"""

# ╔═╡ e54408ea-e984-42b8-926e-74806ffb5aec
# methodswith(Tridiagonal)

# ╔═╡ 017d94a3-c9cb-481f-b3ec-0dc687cb8a46
# methodswith(Dates.Day)

# ╔═╡ 9397e3da-f0b7-43da-9fd5-341a679e46c5
md"""
## The `"*"` operator
"""

# ╔═╡ 32050704-3017-4995-a054-5180f6e6e8d0
# methods(*)

# ╔═╡ 313b771e-95b0-4a6f-add9-e017ae6712a7
md"""
We can multiply various types of numbers and matrices. Notice, however, that there is no multiplication specifically defined for `Tridiagonal` matrices. This would not make much sense, since the product of two tridiagonal matrices is a pentadiagonal matrix, the product of three tridiagonal matrices is septadiagonal matrix, ...

Therefore, two tridiagonal matrices are first converted to full matrices, and then multiplied, as is seen in the source code.
"""

# ╔═╡ e1471fce-4bf3-43d7-8cfc-19e5d4f5b8ec
T₁*T₂

# ╔═╡ 3c8ae987-8cb5-43b4-bda8-bf528e3675af
@which T₁*T₂

# ╔═╡ a04fc04a-c446-4528-abfd-91be7edf0944
T₁*T₂*T₁

# ╔═╡ 8b29ea27-05ef-4b5e-8123-80667c385b70
md"""
## The ⋅ operator
"""

# ╔═╡ 391b705d-0b2d-42b9-86a6-94416d2135cf
# methods(⋅)

# ╔═╡ 0b1557b6-2142-49db-8180-d7d17e1128cf
md"""
By inspecting the source, we see that the `scalar` or the `dot` product of two vectors (1-dimensional arrays) is
computed via `BLAS` function `dot` for real arguments, and the function `dotc` for complex arguments.
"""

# ╔═╡ 76edd9e2-e2c7-4048-b6a4-59f6ec93684a
begin
	x = rand(1:5,5)
	y  = rand(-5:0,5)
	z = rand(5)
	w = rand(5) + im*rand(5)
	@show x, y, z, w
	x⋅y, x⋅z, z⋅x, x⋅w, z⋅w, w⋅z
end

# ╔═╡ b48fd083-d317-4195-8500-c8fa5c5ac4ed
md"""
# `varinfo()`

The command `varinfo()` reveals the content of the specified package or module. It can be invoked either with the package name, or with the package name and a regular expression.
"""

# ╔═╡ 9f93e309-458e-45c8-9a20-dd2179b7b744
# varinfo(Dates)

# ╔═╡ 16ed9927-e764-434a-bbfb-f51fcd339a1f
# varinfo(LinearAlgebra)

# ╔═╡ 7e16a343-d2b6-44c0-9864-46703bb757e5
# Now with a regular expression - we are looking for
# 'eigenvalue' related stuff.
varinfo(LinearAlgebra, Regex("eig"))

# ╔═╡ c084e244-419a-46c1-8bb8-92f58324f7b1
md"""
Funally, let us list all we have in `Julia`'s `Base` module. __It is a long list!__
Notice that `Dates` and `LinearAlgebra` are modules themselves.
"""

# ╔═╡ 198883c4-3f6c-4725-a638-c3be9f7c0c12
# varinfo(Base)

# ╔═╡ a7de4d00-7c4a-11eb-1c76-cd159c6fdf9c
md"
# `code_llvm()` and `code_native()`

It is easy to see the LLVM or assebler code which is executed. The code is displayed in the console.
"

# ╔═╡ d399376e-7c4a-11eb-02e4-936801012749
f(a,b)=2a+b^2

# ╔═╡ e1aa988e-7c4a-11eb-2de6-7b7d5c2fc54c
@code_llvm f(1,2)

# ╔═╡ ef51c630-7c4a-11eb-3331-49602357198f
@code_native f(1,2)

# ╔═╡ 0521934e-7c4b-11eb-37f3-631715a07f3a
@code_native f(1.0,2)

# ╔═╡ 0ba1da00-7c4b-11eb-0abc-c395efc9700e
@code_native f(1,2+2.0im)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.10"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "26b4d16873562469a0a1e6ae41d90dec9e51286d"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.10"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═2db004b0-7c45-11eb-2c3a-3f0e111cf63a
# ╠═5ea33223-4dad-4869-8724-93b7cb91771f
# ╟─8ccdfcf9-d5d1-4b69-a854-15f2715bfafa
# ╟─4440847c-db90-44a5-be48-08c5794ee5b7
# ╠═7a55f84e-1439-45ec-b428-5976cac084cd
# ╟─7514d88c-d0fa-45e5-947d-10fe0ce47f7d
# ╠═611770e1-4b5b-47a6-be3a-e0806990796d
# ╠═9c208420-7cdf-11eb-2b99-e32e66b46e85
# ╠═bfb87460-7cdf-11eb-38ca-0961560c093b
# ╟─cb3e314e-68e0-4fc0-8ece-c5c8bc279f55
# ╠═0f3c49cd-00c7-4654-9dc9-674a89e19a30
# ╠═595878a2-30ab-4a19-a5ad-e621dd37f6c7
# ╟─2a255a09-3c3f-4fcb-99f2-b31adfb536b5
# ╠═4f6d3190-45c9-433d-b7f9-c42ed01dc04b
# ╠═60a3f371-d9b2-47a8-be8f-fab7cfb3fb48
# ╟─74d98ae5-0641-4eea-9f8c-8e663a405f32
# ╠═55f8789f-5361-421f-9f37-98c154fd40cf
# ╠═48783355-31c8-4159-8dd5-8bea047c742c
# ╠═c9a0df42-1368-4e5d-a21d-d924f7ffbb75
# ╟─79ca20e7-446e-4264-ae24-c46fbb4b15c1
# ╟─7e6013ea-1e66-4a81-86a8-a69e22cde29a
# ╠═8591bb73-12e5-474e-b24b-7acfab2e3c5f
# ╟─8c0b2316-26fa-4748-b047-3d2c6bf7acc0
# ╠═4163dd6e-7b5c-42f4-a9f4-c1e1a7faecba
# ╠═89bf9619-7ee7-4ccd-99c0-48583e93bbe1
# ╠═684b3f4a-f83b-440d-9b13-d5dae159e77a
# ╠═243854d4-495a-43de-88cb-1442275c8a8a
# ╟─d79e7c32-ad00-4373-ae0e-04eeb8f6e333
# ╠═22af8627-6713-40a5-b716-81a55a6b7aad
# ╟─98142e4f-c703-4ea9-b893-5e2c831afa42
# ╠═dac4c688-cd5f-4ef5-8ba6-3709f945d3f0
# ╠═60469aa5-8d5b-4996-a02d-100a65e0915f
# ╟─c4bbf3a7-9497-474a-8174-db9530853486
# ╟─88ee9b5a-a1f1-49a1-b6cf-64ebac21634a
# ╠═5e2c532c-49f3-411c-b110-70c54263f007
# ╠═7f026e3e-bf1c-49e5-b815-a9370eb5e3f5
# ╠═20f20bfa-9aff-46a0-b358-5fe42640d21a
# ╟─bd0a3a98-4311-41b0-9209-5e09c4ecd9ad
# ╠═e0a4eeaf-c259-453e-82aa-6813538ab536
# ╠═41daa15d-9b1f-47df-ad24-b124161ce25a
# ╠═cbb653fc-6ec1-408b-b683-7d6821724051
# ╠═e6a068a5-cb01-40f8-8f0c-19e84f8adbc0
# ╟─7622b55d-49ba-4f21-8d85-4bfb95f0fe37
# ╠═3e9cf5d9-f97f-468d-ba27-b6dc217418d8
# ╟─ec7be6ff-07ff-4b6d-bf00-2c48a1535582
# ╠═e54408ea-e984-42b8-926e-74806ffb5aec
# ╠═017d94a3-c9cb-481f-b3ec-0dc687cb8a46
# ╟─9397e3da-f0b7-43da-9fd5-341a679e46c5
# ╠═32050704-3017-4995-a054-5180f6e6e8d0
# ╟─313b771e-95b0-4a6f-add9-e017ae6712a7
# ╠═e1471fce-4bf3-43d7-8cfc-19e5d4f5b8ec
# ╠═3c8ae987-8cb5-43b4-bda8-bf528e3675af
# ╠═a04fc04a-c446-4528-abfd-91be7edf0944
# ╟─8b29ea27-05ef-4b5e-8123-80667c385b70
# ╠═391b705d-0b2d-42b9-86a6-94416d2135cf
# ╟─0b1557b6-2142-49db-8180-d7d17e1128cf
# ╠═76edd9e2-e2c7-4048-b6a4-59f6ec93684a
# ╟─b48fd083-d317-4195-8500-c8fa5c5ac4ed
# ╠═9f93e309-458e-45c8-9a20-dd2179b7b744
# ╠═16ed9927-e764-434a-bbfb-f51fcd339a1f
# ╠═7e16a343-d2b6-44c0-9864-46703bb757e5
# ╟─c084e244-419a-46c1-8bb8-92f58324f7b1
# ╠═198883c4-3f6c-4725-a638-c3be9f7c0c12
# ╟─a7de4d00-7c4a-11eb-1c76-cd159c6fdf9c
# ╠═d399376e-7c4a-11eb-02e4-936801012749
# ╠═e1aa988e-7c4a-11eb-2de6-7b7d5c2fc54c
# ╠═ef51c630-7c4a-11eb-3331-49602357198f
# ╠═0521934e-7c4b-11eb-37f3-631715a07f3a
# ╠═0ba1da00-7c4b-11eb-0abc-c395efc9700e
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: src/Module A - Short Julia Course/L02a Julia is Open.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
using Roots
using Test
import SpecialFunctions.erf

struct SomeInterval
    a
    b
end
Base.extrema(I::SomeInterval) = I.a < I.b ? (I.a, I.b) : (I.b, I.a)

include("./test_find_zero.jl")
include("./test_bracketing.jl")
include("./test_derivative_free.jl")
include("./test_simple.jl")
include("./test_find_zeros.jl")
include("./test_fzero.jl")
include("./test_newton.jl")
include("./test_simple.jl")

#include("./test_composable.jl")
#include("./runbenchmarks.jl")
#include("./test_derivative_free_interactive.jl")


text:
export FrozenLakeEnv

Base.@kwdef mutable struct frozenlakeEnv <: AbstractEnv
	# lake: 3 is the spawn point, 1 is the goal and 2 are obstacles.
	lake::Array{Int8,2} = [1 0 0 0 0
		                   0 2 0 2 0
		                   0 0 0 2 0
						   0 2 0 0 0
		                   0 0 0 0 3]
	reward::Union{Nothing, Float64} = 0
	observation::Int = 24
	is_terminated = false
	num_steps::Int8 = 0
	MAX_STEPS::Int8 = 100
end

RLBase.action_space(env::frozenlakeEnv) = (:up, :right, :down, :left)

begin
    RLBase.state(env::frozenlakeEnv) = env.observation
    RLBase.state_space(env::frozenlakeEnv) = collect(0:24)
    RLBase.is_terminated(env::frozenlakeEnv) = env.is_terminated
end

function RLBase.reset!(env::frozenlakeEnv)
	env.reward = 0
	env.observation = 24
	env.num_steps = 0
	env.is_terminated = false
end

function moveDeterministic(action::String, observation)
	if action == "up"
		observation = clamp((observation÷5)-1, 0, 4)*5 + (observation - 5*(observation÷5))
	elseif action == "down"
		observation = clamp((observation÷5)+1, 0, 4)*5 + (observation - 5*(observation÷5))
	elseif action == "left"
		observation = observation÷5*5 + clamp((observation - 5*(observation÷5))-1, 0, 4)
	elseif action == "right"
		observation = observation÷5*5 + clamp((observation - 5*(observation÷5))+1, 0, 4)
	end
	return observation
end

function (env::frozenlakeEnv)(action)
	r = rand()
	env.num_steps += 1
	if action == :up
		if r < 1/3
			# goes up
			env.observation = moveDeterministic("up", env.observation)
		elseif r < 2/3
			# goes left
			env.observation = moveDeterministic("left", env.observation)
		else
			# goes right
			env.observation = moveDeterministic("left", env.observation)
		end

	elseif action == :right
		if r < 1/3
			# goes right
			env.observation = moveDeterministic("right", env.observation)
		elseif r < 2/3
			# goes up
			env.observation = moveDeterministic("up", env.observation)
		else
			# goes down
			env.observation = moveDeterministic("down", env.observation)
		end

	elseif action == :down
		if r < 1/3
			# goes right
			env.observation = moveDeterministic("right", env.observation)
		elseif r < 2/3
			# goes left
			env.observation = moveDeterministic("left", env.observation)
		else
			# goes down
			env.observation = moveDeterministic("down", env.observation)
		end

	elseif action == :left
		if r < 1/3
			# goes up
			env.observation = moveDeterministic("up", env.observation)
		elseif r < 2/3
			# goes left
			env.observation = moveDeterministic("left", env.observation)
		else
			# goes down
			env.observation = moveDeterministic("down", env.observation)
		end
	else
        @error "unknown action of $action"
	end

	obdiv5 = env.observation÷5
	if env.num_steps>env.MAX_STEPS
		env.is_terminated = true
	end
	if env.lake[obdiv5+1, env.observation - (5*obdiv5) + 1] == 0 || env.lake[obdiv5+1, env.observation - (5*obdiv5) + 1] == 3
		env.reward-=0.01
	elseif env.lake[obdiv5+1, env.observation - (5*obdiv5) + 1] == 2
		env.is_terminated = true
		env.reward-=1
	elseif env.lake[obdiv5+1, env.observation - (5*obdiv5) + 1] == 1
		env.is_terminated = true
		env.reward+=5
    end
	# println(env.reward," ", env.num_steps, " ", env.observation)
end

RLBase.reward(env::frozenlakeEnv) = env.reward

# Env = frozenlakeEnv()
# RLBase.test_runnable!(Env)
# begin
# 	run(RandomPolicy(action_space(Env)), Env, StopAfterEpisode(1))
# end

# begin
#     hook = TotalRewardPerEpisode()
#     run(RandomPolicy(action_space(Env)), Env, StopAfterEpisode(10), hook)
#     plot(hook.rewards)
# end
# p = QBasedPolicy(
#     learner = MonteCarloLearner(;
#             approximator=TabularQApproximator(
#                 n_state = length(state_space(Env)),
#                 n_action = length(action_space(Env)),
#             )
#         ),
#     explorer = EpsilonGreedyExplorer(0.1)
# )


text:
using Documenter, CuArrays

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

function main()
    makedocs(
        sitename = "CuArrays.jl",
        authors = "Tim Besard",
        repo = "$src/blob/{commit}{path}#{line}",
        format = Documenter.HTML(
            # Use clean URLs on CI
            prettyurls = get(ENV, "CI", nothing) == "true",
            canonical = dst,
            assets = ["assets/favicon.ico"],
            analytics = "UA-154489943-3",
        ),
        doctest = false,
        pages = Any[
            "Home" => "index.md",
            "APIs" => Any[
                "memory.md",
            ],
            "Libraries" => Any[
                "lib/blas.md",
                "lib/rand.md",
                "lib/fft.md",
                "lib/solver.md",
                "lib/sparse.md",
                "lib/dnn.md",
                "lib/tensor.md",
            ],
        ]
    )
end

isinteractive() || main()


text:
import Base: Array, convert
import AutoGrad: back
import CUDA: CuArray
using Knet.KnetArrays: DevArray, KnetArray
using AutoGrad: forw, Value, Arg, value

Array(x::Value{K}) where {K<:DevArray}=convert(Array,x)
KnetArray(x::Value{A}) where {A<:AbstractArray}=convert(KnetArray,x)
CuArray(x::Value{A}) where {A<:AbstractArray}=convert(CuArray,x)

convert(::Type{A},x::Value{K}) where {A<:AbstractArray,K<:DevArray}=forw(convert,A,x)
convert(::Type{K},x::Value{A}) where {A<:AbstractArray,K<:DevArray}=forw(convert,K,x)
back(::typeof(convert),::Type{Arg{2}},dy,y,T,x) = convert(typeof(value(x)),dy)

# Array/KnetArray Transfer

# So we will define gradients for convert, KnetArray, Array manually:

# This works but unnecessarily defines new functions:
# cpu2gpu(x::Array)=KnetArray(x)
# @primitive cpu2gpu(x),dy,y (gpu2cpu(dy))
# gpu2cpu(x::KnetArray)=Array(x)
# @primitive gpu2cpu(x),dy,y (cpu2gpu(dy))

# This does not work because !isa(Array,Function)
# @primitive  KnetArray(x::Array),dy  Array(dy)
# @primitive  Array(x::KnetArray),dy  KnetArray(dy)

# This does not work, parametric methods not yet supported, also unnecessary first arg gradient.
# @primitive convert{A<:AbstractArray,K<:KnetArray}(T::Type{K}, x::Value{A}),dy 0 Array(dy)
# @primitive convert{A<:AbstractArray,K<:KnetArray}(T::Type{A}, x::Value{K}),dy 0 KnetArray(dy)

# This gives ambiguity errors:
# @primitive convert(t::Type,x::KnetArray),dy  nothing  convert(KnetArray,dy)
# @primitive convert(t::Type{KnetArray},x::AbstractArray),dy  nothing  convert(Array,dy)


text:
using Assignment04
using Test

@testset "Assignment04" begin

seq = rand(['A','T','G','C'], 20) |> join

@testset "Question 1" begin
    @test normalizeDNA("aatg") == "AATG"
    @test_throws ErrorException normalizeDNA("ZCA")
    @test_throws ErrorException normalizeDNA(42)
    c = normalizeDNA('C')
    @test c == "C"
    @test typeof(c) == String
end

@testset "Question 2" begin
    bc = basecomposition(seq)
    @test bc isa Tuple{Int,Int,Int,Int}
    a,c,g,t = bc

    @test a == count(x-> x == 'A', seq)
    @test c == count(x-> x == 'C', seq)
    @test g == count(x-> x == 'G', seq)
    @test t == count(x-> x == 'T', seq)

    a,c,g,t = basecomposition(lowercase(seq))

    @test a == count(x-> x == 'A', seq)
    @test c == count(x-> x == 'C', seq)
    @test g == count(x-> x == 'G', seq)
    @test t == count(x-> x == 'T', seq)

    @test_throws ErrorException basecomposition(a)
end

@testset "Question 3" begin
    @test gc_content("AATG") == 0.25
    @test gc_content("cccggg") * 100 == 100.0
    @test gc_content("ATta") == 0.0
    @test_throws ErrorException gc_content("ATty")
end

@testset "Question 4" begin
    @test_throws ErrorException kmercount("A", 2)

    a,c,g,t = basecomposition(seq)
    @test kmercount(seq,1) == Dict("A"=>a, "C"=>c, "G"=>g, "T"=>t)
    for k in 3:5
        @test sum(values(kmercount(seq, k))) == 20 - k + 1
    end

    @test kmercount("GCGCGCGCG", 5) == Dict("CGCGC" => 2, "GCGCG" => 3)
end

end # Assignment04


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Mean/variance calculations #
# ========================== #
function outcome_mean(subproblems::Vector{JuMP.Model}, probabilities::AbstractVector, sense::Union{MOI.OptimizationSense, Nothing} = nothing)
    # Sanity check
    length(subproblems) == length(probabilities) || error("Inconsistent number of subproblems and probabilities")
    N = length(subproblems)
    N == 0 && return 0.0
    if sense === nothing || sense == MOI.FEASIBILITY_SENSE
        sense = objective_sense(subproblems[1])
    end
    return mapreduce(+, 1:N) do k
        outcome = subproblems[k]
        val = try
            optimize!(outcome)
            status = termination_status(outcome)
            if !(status ∈ AcceptableTermination)
                if status == MOI.INFEASIBLE
                    objective_sense(outcome) == MOI.MAX_SENSE ? -Inf : Inf
                elseif status == MOI.DUAL_INFEASIBLE
                    objective_sense(outcome) == MOI.MAX_SENSE ? Inf : -Inf
                else
                    error("Outcome model could not be solved, returned status: $status")
                end
            else
                probabilities[k]*objective_value(outcome)
            end
        catch error
            if isa(error, NoOptimizer)
                @warn "No optimizer set, cannot solve outcome model."
                rethrow(NoOptimizer())
            else
                @warn "Outcome model could not be solved."
                rethrow(error)
            end
        end
        val = sense == objective_sense(outcome) ? val : -val
    end
end
function welford(subproblems::Vector{JuMP.Model}, probabilities::AbstractVector, sense::Union{MOI.OptimizationSense, Nothing} = nothing)
    # Sanity check
    length(subproblems) == length(probabilities) || error("Inconsistent number of subproblems and probabilities")
    N = length(subproblems)
    N == 0 && return 0.0, 0.0, 0.0, N
    if sense === nothing || sense == MOI.FEASIBILITY_SENSE
        sense = objective_sense(subproblems[1])
    end
    Q̄ₖ = Sₖ = wₖ = w²ₖ = 0
    for k = 1:N
        π = probabilities[k]
        wₖ = wₖ + π
        Q̄ₖ₋₁ = Q̄ₖ
        problem = subproblems[k]
        try
            optimize!(problem)
            status = termination_status(problem)
            Q = if !(status ∈ AcceptableTermination)
                Q = if status == MOI.INFEASIBLE
                    objective_sense(problem) == MOI.MAX_SENSE ? -Inf : Inf
                elseif status == MOI.DUAL_INFEASIBLE
                    objective_sense(problem) == MOI.MAX_SENSE ? Inf : -Inf
                else
                    error("Subproblem could not be solved during welford algorithm, returned status: $status")
                end
            else
                Q = objective_value(problem)
            end
            Q = sense == objective_sense(problem) ? Q : -Q
            Q̄ₖ = Q̄ₖ + (π / wₖ) * (Q - Q̄ₖ)
            Sₖ = Sₖ + π * (Q - Q̄ₖ) * (Q - Q̄ₖ₋₁)
        catch error
            if isa(error, NoOptimizer)
                @warn "No optimizer set, cannot solve subproblem during welford algorithm."
                rethrow(NoOptimizer())
            else
                @warn "Subproblem could not be solved during welford algorithm."
                rethrow(error)
            end
        end
    end
    correction = N / ((N - 1) * wₖ)
    if isinf(Q̄ₖ)
        correction = 0.0
    end
    return Q̄ₖ, correction * Sₖ, wₖ, N
end
function aggregate_welford(left::Tuple, right::Tuple)
    x̄ₗ, σₗ², wₗ, nₗ = left
    x̄ᵣ, σᵣ², wᵣ, nᵣ = right
    δ = x̄ᵣ - x̄ₗ
    w = wₗ + wᵣ
    n = nₗ + nᵣ
    x̄ = (wₗ * x̄ₗ + wᵣ * x̄ᵣ) / w
    Sₗ = σₗ² * ((nₗ - 1) * wₗ) / nₗ
    Sᵣ = σᵣ² * ((nᵣ - 1) * wᵣ) / nᵣ
    S = Sₗ + Sᵣ + (wₗ * wᵣ / w) * δ * δ
    correction = n / ((n - 1) * w)
    return x̄, correction * S, w, n
end


text:
"""
XDR (machine-independent binary) RData format IO stream wrapper.
"""
struct XDRIO{T<:IO} <: RDAIO
    sub::T             # underlying IO stream
    buf::Vector{UInt8} # buffer for strings

    XDRIO(io::T) where {T <: IO} = new{T}(io, Vector{UInt8}(undef, 1024))
end

readint32(io::XDRIO) = ntoh(read(io.sub, Int32))
readuint32(io::XDRIO) = ntoh(read(io.sub, UInt32))
readfloat64(io::XDRIO) = ntoh(read(io.sub, Float64))

readintorNA(io::XDRIO) = readint32(io)
function readintorNA(io::XDRIO, n::RVecLength)
    v = read!(io.sub, Vector{Int32}(undef, n))
    map!(ntoh, v, v)
end

# this method have Win32 ABI issues, see JuliaData/RData.jl#5
# R's NA is silently converted to NaN when the value is loaded in the register(?)
#readfloatorNA(io::XDRIO) = readfloat64(io)
function readfloatorNA(io::XDRIO, n::RVecLength)
    v = read!(io.sub, Vector{Float64}(undef, n))
    map!(ntoh, v, v)
end

function readfloatorNA!(io::XDRIO, v::AbstractVector{Float64})
    readbytes!(io.sub, reinterpret(UInt8, v))
    map!(ntoh, v, v)
end

readuint8(io::XDRIO, n::RVecLength) = read!(io.sub, Vector{UInt8}(undef, n))

function readnchars(io::XDRIO, n::Int32)  # a single character string
    readbytes!(io.sub, io.buf, n)
    convert(RString, unsafe_string(pointer(io.buf), n))
end


text:
module SteadyState

using QuantumOptics
using IterativeSolvers, LinearMaps, LinearAlgebra, SparseArrays

include("generic_method.jl")
export steadystate_iterative!
include("bicgstab.jl")
export steadystate_bicg, steadystate_bicg!
#include("indexing.jl")

end # module


text:
using Thebes, Luxor

function mygfunction(vertices, faces, labels; action=:fill)
    if !isempty(faces)
        @layer begin
            for (n, p) in enumerate(faces)
                @layer begin
                    setopacity(0.5)
                    randomhue()
                    poly(p, action)
                end

                sethue("black")
                setline(0.5)
                poly(p, :stroke, close=true)

            end
        end
    end
end

@svg begin
    background("black")
    setlinejoin("bevel")
    eyepoint(Point3D(150, 150, 20))
    sethue("pink")
    carpet(200)
    axes3D(20)

    #object = make(Cube, "Cube")
    #object = make(Tetrahedron, "Tetrahedron")
    object = make(Pyramid, "Pyramid")

    scaleby!(object, 100, 100, 100)
    rotateby!(object, Point3D(0, 0, 0), 0, 0, rand())
    sortfaces!(object)
    sethue("orange")
    pin(object, gfunction= mygfunction)
end 500 500 "/tmp/cube.svg"


text:
using HorizonSideRobots
include("roblib.jl")

function inverse(side::HorizonSide)
    HorizonSide(mod(Int(side)+2, 4))
end

function tuda(r)
    for i in 1:12
        obhod(r, Ost)
        putmarker!(r)
    end
    try_move!(r, Nord)
    putmarker!(r)
end

function suda(r)
    for i in 1:12
        obhod(r, West)
        putmarker!(r)
    end
    try_move!(r, Nord)
    putmarker!(r)
end

function mark_all_2(r::Robot, side::HorizonSide)
    while isborder(r, Nord)==false && isborder(r, Ost)==false
        tuda(r, side)
        suda(r, inverse(side))
    end
end

mark_all_2(r, Ost)


text:
# Some discussion on this:
# https://groups.google.com/forum/#!msg/julia-users/YP31LM3Qto0/ET-XjN-vQuAJ
#
# All the model parameters
#
# TODO: improve macro hygiene.

__precompile__()

module Parameters
import Base: @__doc__
import DataStructures: OrderedDict
using Compat

export @with_kw, type2dict, reconstruct, @unpack, @pack

## Parser helpers
#################

# To iterate over code blocks dropping the line-number bits:
immutable Lines
    block::Expr
end
Base.start(lns::Lines) = 1
function Base.next(lns::Lines, nr)
    for i=nr:length(lns.block.args)
        if isa(lns.block.args[i], LineNumberNode)
            continue
        end
        if isa(lns.block.args[i], Symbol) || !(lns.block.args[i].head==:line)
            return lns.block.args[i], i+1
        end
    end
    return -1
end
function Base.done(lns::Lines, nr)
    if next(lns::Lines, nr)==-1
        true
    else
        false
    end
end

# Transforms :(a::b) -> :a
decolon2(a::Expr) = (@assert a.head==:(::);  a.args[1])
decolon2(a::Symbol) = a

# Returns the name of the type as Symbol
function typename(typedef::Expr)
    if isa(typedef.args[2], Symbol)
        return typedef.args[2]
    elseif isa(typedef.args[2].args[1], Symbol)
        return typedef.args[2].args[1]
    elseif isa(typedef.args[2].args[1].args[1], Symbol)
        return typedef.args[2].args[1].args[1]
    else
        error("Could not parse type-head from: $typedef")
    end
end

# Transforms:  Expr(:<:, :A, :B) -> :A
stripsubtypes(s::Symbol) = s
function stripsubtypes(e::Expr)
    e.args[1]
end
stripsubtypes(vec::Vector) = [stripsubtypes(v) for v in vec]

function check_inner_constructor(l)
    if length(l.args[1].args)==1
        error("No inner constructors with zero positional arguments allowed!")
    elseif (length(l.args[1].args)==2 #1<length(l.args[1].args)<=3
            && isa(l.args[1].args[2], Expr)
            && l.args[1].args[2].head==:parameters)
        error("No inner constructors with zero positional arguments plus keyword arguments allowed!")
    end
    nothing
end

## exported functions
#####################
"""
Transforms a type-instance into a dictionary.

```
julia> type T
           a
           b
       end

julia> type2dict(T(4,5))
Dict{Symbol,Any} with 2 entries:
  :a => 4
  :b => 5
```
"""
function type2dict(dt)
    di = Dict{Symbol,Any}()
    for n in @compat fieldnames(dt)
        di[n] = getfield(dt, n)
    end
    di
end

"""
Make a new instance of a type with the same values as
the input type except for the fields given in the associative
second argument or as keywords.

```julia
type A; a; b end
a = A(3,4)
b = reconstruct(a, [(:b, 99)]) # ==A(3,99)
```
"""
function reconstruct{T}(pp::T, di)
    di = !isa(di, Associative) ? Dict(di) : di
    ns = @compat fieldnames(pp)
    args = Array(Any, length(ns))
    for (i,n) in enumerate(ns)
        args[i] = get(di, n, getfield(pp, n))
    end
    T(args...)
end
reconstruct{T}(pp::T; kws...) = reconstruct(pp, kws)

# A type with fields (r,a) in variable aa becomes
# quote
#     r = aa.r
#     a = aa.a
# end
_unpack(binding, fields) = Expr(:block, [:($f = $binding.$f) for f in fields]...)
# Pack fields back into binding using reconstruct:
function _pack(binding, fields)
    kws = [Expr(:kw, f, f) for f in fields]
    :($binding = Main.Parameters.reconstruct($binding, $(kws...)) )
end

"""
This function is called by the `@with_kw` macro and does the AST transformation from:

```julia
@with_kw immutable MM{R}
    r::R = 1000.
    a::R
end
```

into

```julia
immutable MM{R}
    r::R
    a::R
    MM(r,a) = new(r,a)
    MM(;r=1000., a=error("no default for a")) = MM{R}(r,a)
end
MM(m::MM; kws...) = reconstruct(mm,kws)
MM(m::MM, di::Union{Associative, Tuple{Symbol,Any}}) = reconstruct(mm, di)
macro unpack_MM(varname)
    esc(quote
    r = varname.r
    a = varname.a
    end)
end
macro pack_MM(varname)
    esc(quote
    varname = Main.Parameters.reconstruct(varname,r=r,a=a)
    end)
end
```
"""
function with_kw(typedef)
    if typedef.head!=:type
        error("only works on type-defs")
    end
    const err1str = "Field \'"
    const err2str = "\' has no default, supply it with keyword."

    inner_constructors = Any[]

    # parse a few things
    tn = typename(typedef) # the name of the type
    # Returns M{...} (removes any supertypes)
    if isa(typedef.args[2], Symbol)
        typparas = Any[]
    elseif typedef.args[2].head==:<:
        if isa(typedef.args[2].args[1],Symbol)
            typparas = Any[]
        else
            typparas = typedef.args[2].args[1].args[2:end]
        end
    else
        typparas = typedef.args[2].args[2:end]
    end

    # error on types without fields
    lns = Lines(typedef.args[3])
    if done(lns, start(lns))
        error("@with_kw only supported for types which have at least one field.")
    end
    # default type @deftype
    l, i = next(lns, start(lns))
    if isa(l, Expr) && l.head==:macrocall && l.args[1]==Symbol("@deftype")
        has_deftyp = true
        if length(l.args) != 2
            error("Malformed `@deftype` line")
        end
        deftyp = l.args[2]
        if done(lns, i)
            error("@with_kw only supported for types which have at least one field.")
        end
    else
        has_deftyp = false
    end

    # the vars for the unpack macro
    unpack_vars = Any[]
    # the type def
    fielddefs = quote end # holds r::R etc
    fielddefs.args = Any[] # in julia 0.5 this is [:( # /home/mauro/.julia/v0.5/Parameters/src/Parameters.jl, line 228:)]
    kws = OrderedDict{Any, Any}()
    # assertions in the body
    asserts = Any[]
    for (i,l) in enumerate(lns) # loop over body of typedef
        if i==1 && has_deftyp
            continue
        end
        if isa(l, Symbol)  # no default value and no type annotation
            if has_deftyp
                push!(fielddefs.args, :($l::$deftyp))
            else
                push!(fielddefs.args, l)
            end
            sym = l
            syms = string(sym)
            kws[sym] = :(error($err1str * $syms * $err2str))
            # unwrap-macro
            push!(unpack_vars, sym)
        elseif l.head==:(=)  # default value and with or without type annotation
            if isa(l.args[1], Expr) && l.args[1].head==:call # inner constructors
                check_inner_constructor(l)
                push!(inner_constructors, l)
            else
                fld = l.args[1]
                if isa(fld, Symbol) && has_deftyp # no type annotation
                    fld = :($fld::$deftyp)
                end
                push!(fielddefs.args, fld)
                kws[decolon2(fld)] = l.args[2]
                # unwrap-macro
                push!(unpack_vars, decolon2(fld))
            end
        elseif l.head==:macrocall  && l.args[1]==Symbol("@assert")
            # store all asserts
            push!(asserts, l)
        elseif l.head==:function # inner constructor
            check_inner_constructor(l)
            push!(inner_constructors, l)
        else # no default value but with type annotation
            push!(fielddefs.args, l)
            sym = decolon2(l.args[1])
            syms = string(sym)
            kws[sym] = :(error($err1str *$syms * $err2str))
            # unwrap-macro
            push!(unpack_vars, l.args[1])
        end
    end
    # The type definition without inner constructors:
    typ = Expr(:type, deepcopy(typedef.args[1:2])..., deepcopy(fielddefs))

    # Inner keyword constructor.  Note that this calls the positional
    # constructor under the hood and not `new`.  That way a user can
    # provide a special positional constructor (say enforcing
    # invariants) which also gets used with the keywords.
    args = Any[]
    kwargs = Expr(:parameters)
    for (k,w) in kws
        push!(args, k)
        push!(kwargs.args, Expr(:kw,k,w))
    end
    if length(typparas)>0
        innerc = :($tn($kwargs) = $tn{$(stripsubtypes(typparas)...)}($(args...)) )
    else
        innerc = :($tn($kwargs) = $tn($(args...)) )
    end
    push!(typ.args[3].args, innerc)

    # Inner positional constructor: only make it if no inner
    # constructors are user-defined.  If one or several are defined,
    # assume that one has the standard positional signature.
    if length(inner_constructors)==0
        innerc2 = :($tn($(args...)) = new($(args...)))
        prepend!(innerc2.args[2].args, asserts)
        push!(typ.args[3].args, innerc2)
    else
        if length(asserts)>0
            error("Assertions are only allowed in type-definitions with no inner constructors.")
        end
        append!(typ.args[3].args, inner_constructors)
    end

    # Do not provide an outer keyword-constructor as the type
    # parameters need to be given explicitly anyway.

    # Outer positional constructor which does not need explicit
    # type-parameters when called.  Only make this constructor if
    #  (1) type parameters are used at all
    #  (2) all type parameters are used in the fields (otherwise get a
    #      "method is not callable" warning!)
    if typparas!=Any[] # condition (1)
        outer_positional = :(  $tn{$(typparas...)}($(fielddefs.args...))
                             = $tn{$(stripsubtypes(typparas)...)}($(args...)))
        # Check condition (2)
        used_paras = Any[]
        for f in fielddefs.args
            if !isa(f,Symbol) && f.head==:(::)
                push!(used_paras, f.args[2])
            end
        end
        if !issubset(stripsubtypes(typparas), used_paras)
            outer_positional = :()
        end
    else
        outer_positional = :()
    end

    ## outer copy constructor
    ###
    outer_copy = quote
        $tn(pp::$tn; kws... ) = reconstruct(pp, kws)
        # $tn(pp::$tn, di::Union(Associative,Vararg{Tuple{Symbol,Any}}) ) = reconstruct(pp, di) # see issue https://github.com/JuliaLang/julia/issues/11537
        # $tn(pp::$tn, di::Union(Associative, Tuple{Vararg{Tuple{Symbol, Any}}}) ) = reconstruct(pp, di) # see issue https://github.com/JuliaLang/julia/issues/11537
        $tn(pp::$tn, di::Associative ) = reconstruct(pp, di)
        $tn(pp::$tn, di::Vararg{Tuple{Symbol,Any}} ) = reconstruct(pp, di)
    end

    # (un)pack macro from https://groups.google.com/d/msg/julia-users/IQS2mT1ITwU/hDtlV7K1elsJ
    unpack_name = Symbol("unpack_"*string(tn))
    pack_name = Symbol("pack_"*string(tn))
    # Finish up
    quote
        Main.Parameters.@__doc__ $typ # use Main.Parameters.@__doc__ for 0.3 compatibility
        $outer_positional
        $outer_copy
        function Base.show(io::IO, p::$tn)
            println(io, string(typeof(p)))
            for (i, var) in enumerate($unpack_vars)
                print(io, "  " * string(var) * ": $(getfield(p,var))")
                i == length($unpack_vars) || print(io, "\n")
            end
        end
        macro $unpack_name(ex)
            esc(Main.Parameters._unpack(ex, $unpack_vars))
        end
        macro $pack_name(ex)
            esc(Main.Parameters._pack(ex, $unpack_vars))
        end
    end
end

"""
Macro which allows default values for field types and a few other features.

Basic usage:

```julia
@with_kw immutable MM{R}
    r::R = 1000.
    a::Int = 4
end
```

For more details see manual.
"""
macro with_kw(typedef)
    return esc(with_kw(typedef))
end

## @pack and @unpack are independent of the datatype
function parse_pack_unpack(arg)
    if isa(arg, Symbol); error("Need format `t: a, ...`") end
    h = arg.head
    if !(h==:(:) || h==:tuple); error("Need format `t: a, ...`") end
    # var-name of structure
    v = h==:tuple ? arg.args[1].args[1] : arg.args[1]
    # vars to unpack
    up = Any[]
    if h==:tuple
        append!(up, arg.args[2:end])
        push!(up, arg.args[1].args[2])
    else
        push!(up, arg.args[2])
    end
    return v, up # variable holding the structure, variables to (un)-pack
end

"""
Unpacks fields from any datatype (no need to create it with @with_kw):

```julia
type A
    a
    b
end
aa = A(3,4)
@unpack aa: a,b
# is equivalent to
a = aa.a
b = aa.b
```
"""
macro unpack(arg)
    v, up = parse_pack_unpack(arg)
    out = quote end
    for u in up
        push!(out.args, :($u = $v.$u))
    end
    return esc(out)
end

"""
Packs values into a datatype.  The variables need to have the same
name as the fields.  If the datatype is mutable, it will be mutated.
If immutable, a new instance is made with `reconstruct` and assigned
to the original variable.

```julia
type A
    a
    b
end
aa = A(3,4)
b = "ha"
@pack aa: b
# is equivalent to
aa.b = b
```
"""
macro pack(arg)
    v, up = parse_pack_unpack(arg)
    # dict to use with reconstruct:
    di = :([])
    for u in up
        push!(di.args, :($(Base.Meta.quot(u)), $u))
    end
    # assignments for mutables
    ass = quote end
    for u in up
        push!(ass.args, :($v.$u = $u))
    end
    esc(
        quote
        if isimmutable($v)
            $v = Main.Parameters.reconstruct($v, $di)
        else
            $ass
        end
        end
    )
end

end # module


text:
module YouBotWithSphere

using ModiaLang
using Unitful
import Modia3D

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

# some constants

# visual materials
vmat1 = VisualMaterial(color="LightBlue"      , transparency=0.0)  # material of FileMesh
vmat2 = VisualMaterial(color="Grey"           , transparency=0.3)  # material of work piece
vmat3 = VisualMaterial(color="DodgerBlue3"    , transparency=0.0)  # material of table
vmat4 = VisualMaterial(color="LightGoldenRod1", transparency=0.0)  # material of ground
workpiece_Lx = 0.04
workpiece_Ly = 0.04
workpiece_Lz = 0.08
xPosTable = 0.71
tableX    = 0.3
tableY    = 0.3
tableZ    = 0.01
heigthLeg = 0.37
widthLeg  = 0.02
diameter = 0.05
right_finger_offset=0.0

# all needed paths to fileMeshes
base_frame_obj           = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/base_frame.obj")
plate_obj                = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/plate.obj")
front_right_wheel_obj    = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/front-right_wheel.obj")
front_left_wheel_obj     = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/front-left_wheel.obj")
back_right_wheel_obj     = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/back-right_wheel.obj")
back_left_wheel_obj      = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/back-left_wheel.obj")
arm_base_frame_obj       = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_base_frame.obj")
arm_joint_1_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_1.obj")
arm_joint_2_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_2.obj")
arm_joint_3_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_3.obj")
arm_joint_4_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_4.obj")
arm_joint_5_obj          = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_5.obj")
gripper_base_frame_obj   = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/gripper_base_frame.obj")
gripper_left_finger_obj  = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/gripper_left_finger.obj")
gripper_right_finger_obj = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/gripper_right_finger.obj")

# Drive train data: motor inertias and gear ratios
nullRot = nothing #SMatrix{3,3,Float64,9}(Matrix(1.0I, 3, 3))

motorInertia1 = 0.0000135 + 0.000000409
motorInertia2 = 0.0000135 + 0.000000409
motorInertia3 = 0.0000135 + 0.000000071
motorInertia4 = 0.00000925 + 0.000000071
motorInertia5 = 0.0000035 + 0.000000071
gearRatio1    = 156.0
gearRatio2    = 156.0
gearRatio3    = 100.0
gearRatio4    = 71.0
gearRatio5    = 71.0

m1=1.390
translation1 =[0.033,0,0]
rotation1 =Modia3D.rot1(180u"°")

m2=1.318
translation2=[0.155,0,0]
rotation2=Modia3D.rot123(90u"°", 0.0, -90u"°")

m3=0.821
translation3=[0,0.135,0]
rotation3=Modia3D.rot3(-90u"°")

m4=0.769
translation4=[0,0.11316,0]

m5=0.687
translation5=[0,0,0.05716]
rotation5=Modia3D.rot1(-90u"°")

### ----------------- Servo Model -----------------------
# parameters for Link
axisLink = 3
dLink  = 1.0
k1Link = 50.0
k2Link = 0.1
T2Link = 1.0
# parameters for Gripper
axisGripper = 2
dGripper  = 1.0
k1Gripper = 50.0
k2Gripper = 0.1
T2Gripper = 1.0
motorInertiaGripper = 0.0
gearRatioGripper    = 1.0

#### ----------- Path Planning ------------------
referencePath1 = Modia3D.ReferencePath(
    names = ["angle1", "angle2", "angle3", "angle4", "angle5", "gripper"],
    position = [0.0,     0.0,     pi/2,   0.0,    0.0,   0.0],
    v_max =    [2.68512, 2.68512, 4.8879, 5.8997, 5.8997, 1.0],
    a_max =    [1.5, 1.5, 1.5, 1.5, 1.5, 0.5])
#
Modia3D.ptpJointSpace(referencePath = referencePath1, positions =
    [0.0 0.0 pi/2     0.0 0.0 0.0;
    0.0  0.3 pi/2-0.3 0.0 0.0 0.0;
    0.0  0.0 0.0      0.0 0.0 0.0])

getReferencePath() = referencePath1

# Controller Model
Controller = Model(
    # Interface
    refLoadAngle     = input,
    actualMotorAngle = input,
    actualMotorSpeed = input,
    refTorque        = output,
    # Components
    gainOuter     = Gain,
    angleFeedback = Feedback,
    gainInner     = Gain,
    speedFeedback = Feedback,
    PI = PI,

    connect = :[
        (refLoadAngle, gainOuter.u)
        (gainOuter.y, angleFeedback.u1)
        (actualMotorAngle, angleFeedback.u2)
        (angleFeedback.y, gainInner.u)
        (gainInner.y, speedFeedback.u1)
        (actualMotorSpeed, speedFeedback.u2)
        (speedFeedback.y, PI.u)
        (PI.y, refTorque)]
)

Drive = Model(
    # Interface
    refTorque   = input,
    actualSpeed = output,
    actualAngle = output,
    flange      = Flange,
    # Components
    torque       = UnitlessTorque,
    speedSensor  = UnitlessSpeedSensor,
    angleSensor  = UnitlessAngleSensor,
    motorInertia = Inertia,
    idealGear    = IdealGear,

    connect = :[
        (refTorque, torque.tau)
        (torque.flange, speedSensor.flange, angleSensor.flange, motorInertia.flange_a)
        (motorInertia.flange_b, idealGear.flange_a)
        (idealGear.flange_b, flange)
        (speedSensor.w, actualSpeed)
        (angleSensor.phi, actualAngle)
        ]
)

Servo = Model(
    refLoadAngle = input,
    flange       = Flange,
    ppi          = Controller,
    drive        = Drive,

    connect = :[
        # inputs of ppi
        (refLoadAngle, ppi.refLoadAngle)
        (ppi.actualMotorSpeed, drive.actualSpeed)
        (ppi.actualMotorAngle, drive.actualAngle)
        # output of ppi --> input of motor
        (ppi.refTorque, drive.refTorque)
        # Flange of drive
        (drive.flange, flange)
        ]
)

Ground = Model(
    ground = Object3D(parent=:world,
        translation=[0.5, 0.0, -0.005], feature = Solid(shape =Box(lengthX=2.5, lengthY=2.0, lengthZ=0.01), collisionSmoothingRadius=0.001, solidMaterial="DryWood", visualMaterial=vmat4, collision=true))
)

Table = Model(
    plate = Object3D(parent=:world, translation=[xPosTable, 0.0, heigthLeg],
        feature=Solid(shape=Box(lengthX=tableX, lengthY=tableY, lengthZ=tableZ),
        collisionSmoothingRadius=0.001, solidMaterial="DryWood", collision=true, visualMaterial=vmat3)),
    leg1  = Object3D(parent=:plate, translation=[ tableX/2-widthLeg/2,  tableY/2-widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg),
        solidMaterial="DryWood", visualMaterial=vmat3)),
    leg2  = Object3D(parent=:plate, translation=[-tableX/2+widthLeg/2,  tableY/2-widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg), solidMaterial="DryWood", visualMaterial=vmat3)),
    leg3  = Object3D(parent=:plate, translation=[ tableX/2-widthLeg/2, -tableY/2+widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg), solidMaterial="DryWood", visualMaterial=vmat3)),
    leg4  = Object3D(parent=:plate, translation=[-tableX/2+widthLeg/2, -tableY/2+widthLeg/2, -heigthLeg/2],
        feature=Solid(shape=Box(lengthX=widthLeg, lengthY=widthLeg, lengthZ=heigthLeg), solidMaterial="DryWood", visualMaterial=vmat3))
)

Base = Model(
    base_frame = Object3D(parent=:world, translation=[0, 0, 0.084],
        feature=Solid(shape=FileMesh(filename = base_frame_obj), massProperties=MassPropertiesFromShapeAndMass(mass=19.803), visualMaterial=vmat1 )),
    plate = Object3D(parent=:base_frame, translation=[-0.159, 0, 0.046],
        feature=Solid(shape=FileMesh(filename = plate_obj), massProperties=MassPropertiesFromShapeAndMass(mass=2.397), visualMaterial=vmat1)),
    front_right_wheel = Object3D(parent=:base_frame, translation=[0.228, -0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = front_right_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4), visualMaterial=vmat1)),
    front_left_wheel = Object3D(parent=:base_frame, translation=[0.228, 0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = front_left_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4), visualMaterial=vmat1)),
    back_right_wheel = Object3D(parent=:base_frame, translation=[-0.228, -0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = back_right_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4), visualMaterial=vmat1)),
    back_left_wheel   = Object3D(parent=:base_frame, translation=[-0.228, 0.158, -0.034],
        feature=Solid(shape=FileMesh(filename = back_left_wheel_obj), massProperties=MassPropertiesFromShapeAndMass(mass=1.4), visualMaterial=vmat1))
)

servoParameters1 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio1),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia1),
                          idealGear = Map(ratio=gearRatio1)
                      )
                  )

servoParameters2 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio2),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia2),
                          idealGear = Map(ratio=gearRatio2)
                      )
                  )

servoParameters3 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio3),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia3),
                          idealGear = Map(ratio=gearRatio3)
                      )
                  )

servoParameters4 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio4),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia4),
                          idealGear = Map(ratio=gearRatio4)
                      )
                  )

servoParameters5 = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio5),
                          gainInner = Map(k=k1Link),
                          PI        = Map(k=k2Link, T=T2Link)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia5),
                          idealGear = Map(ratio=gearRatio5)
                      )
                  )

featureBody1 = Solid(shape = FileMesh(filename = arm_joint_1_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m1), visualMaterial=vmat1)
featureBody2 = Solid(shape = FileMesh(filename = arm_joint_2_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m2), visualMaterial=vmat1)
featureBody3 = Solid(shape = FileMesh(filename = arm_joint_3_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m3), visualMaterial=vmat1)
featureBody4 = Solid(shape = FileMesh(filename = arm_joint_4_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m4), visualMaterial=vmat1)
featureBody5 = Solid(shape = FileMesh(filename = arm_joint_5_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m5), visualMaterial=vmat1)

linkParameters1 = Map(
                    parent1 = Par(value = :(armBase_b)),
                    featureBody = featureBody1, # Par(value = :featureBody1),
                    initRefPos = referencePath1.position[1],
                    trans = translation1,
                    rota = Par(value = :(rotation1))
)

linkParameters2 = Map(
                    parent1 = Par(value = :(link1.obj2)),
                    featureBody = featureBody2, #Par(value = :featureBody2),
                    m = m2,
                    initRefPos = referencePath1.position[2],
                    trans = translation2,
                    rota = Par(value = :(rotation2))
)

linkParameters3 = Map(
                    parent1 = Par(value = :(link2.obj2)),
                    featureBody = featureBody3, # Par(value = :featureBody3),
                    m = m3,
                    initRefPos = referencePath1.position[3],
                    trans = translation3,
                    rota = Par(value = :(rotation3))
)

linkParameters4 = Map(
                    parent1 = Par(value = :(link3.obj2)),
                    featureBody = featureBody4, # Par(value = :featureBody4),
                    m = m4,
                    initRefPos = referencePath1.position[4],
                    trans = translation4,
                    rota = Par(value = :(nullRot))
)

linkParameters5 = Map(
                    parent1 = Par(value = :(link4.obj2)),
                    featureBody = featureBody5, # Par(value = :featureBody5),
                    m = m5,
                    initRefPos = referencePath1.position[5],
                    trans = translation5,
                    rota = Par(value = :(rotation5))
)

featureVisual1 = Visual(shape=Cylinder(axis=3, diameter=0.01, length=0.12), visualMaterial=vmat1)

Link = Model(
    parent1 = Par(value = :(nothing)),
    featureBody = featureBody1,
    featureVisual = featureVisual1,
    #m = 0.0,
    #initRefPos = Par(value = :(0.0)),
    #index = 1,
    trans = [0,0,0],
    rota = Par(value = :(nullRot)), #rotation1, #Modia3D.rot1(180u"°"),

    obj1 = Object3D(parent=:parent1, rotation=:rota, feature=:featureVisual),    #, visualizeFrame=false),
    body = Object3D(feature = :featureBody),
    obj2 = Object3D(parent =:body, translation = :trans),
)

Gripper = Model(
    obj1 = Object3D(parent=:(link5.obj2), rotation=Modia3D.rot3(-180u"°")),
    gripper_base_frame = Object3D(parent=:obj1, feature=Solid(shape=
        FileMesh(filename = gripper_base_frame_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.199), visualMaterial=vmat1)),
    gripper_left_finger_a = Object3D(parent=:gripper_base_frame),
    gripper_left_finger = Object3D(parent=:gripper_left_finger_a,
        translation=[0, 0.0082, 0], feature=Solid(shape= FileMesh(filename = gripper_left_finger_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.010), visualMaterial=vmat1, contactMaterial="DryWood", collision=true)),
    gripper_right_finger_a = Object3D(parent=:gripper_base_frame,
        translation=[0,-right_finger_offset,0]),
    gripper_right_finger = Object3D(parent=:gripper_right_finger_a, translation=[0, -0.0082, 0],
        feature = Solid(shape = FileMesh(filename = gripper_right_finger_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.010), visualMaterial=vmat1, contactMaterial="DryWood", collision=true))
    ###  needs to be done
    #prism = ControlledPrismatic(gripper_right_finger_a, gripper_left_finger_a, axis=2,
    #                                d=d, k1=k1, k2=k2, T2=T2, gearRatio=gearRatio, referencePath=referencePath,  ref_index=6)
)

YouBot = Model(
    gravField = UniformGravityField(g=9.81, n=[0,0,-1]),
    world = Object3D(feature=Scene(gravityField = :gravField, visualizeFrames=false,
    defaultFrameLength=0.1, enableContactDetection=true, elasticContactReductionFactor=1e-4)),

    base = Base,
    table = Table,
    ground = Ground,

    sphere = Object3D(
        feature=Solid(shape=Sphere(diameter=diameter), visualMaterial=vmat2, solidMaterial="DryWood", collision=true)),

    free = FreeMotion(obj1=:(table.plate), obj2=:sphere, r=Var(init=[-tableX/2+diameter/2, 0.0, diameter/2+tableZ/2]), rot=Var(init=[pi/2, 0.0, 0.0])),

    arm_base_frame = Object3D(parent=:(base.base_frame),
        translation=[0.143, 0.0, 0.046],
        feature = Solid(shape = FileMesh(filename = arm_base_frame_obj), massProperties=MassPropertiesFromShapeAndMass(mass=0.961), visualMaterial=vmat1)),
    armBase_b = Object3D(parent=:arm_base_frame, translation=[0.024, 0, 0.115]),

    link1 = Link | linkParameters1,
    link2 = Link | linkParameters2,
    link3 = Link | linkParameters3,
    link4 = Link | linkParameters4,
    link5 = Link | linkParameters5,
    gripper = Gripper,

    rev1 = RevoluteWithFlange(obj1 = :(link1.obj1), obj2 = :(link1.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 1)), w=Var(init=0.0)),
    rev2 = RevoluteWithFlange(obj1 = :(link2.obj1), obj2 = :(link2.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 1)), w=Var(init=0.0)),
    rev3 = RevoluteWithFlange(obj1 = :(link3.obj1), obj2 = :(link3.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 1)), w=Var(init=0.0)),
    rev4 = RevoluteWithFlange(obj1 = :(link4.obj1), obj2 = :(link4.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 1)), w=Var(init=0.0)),
    rev5 = RevoluteWithFlange(obj1 = :(link5.obj1), obj2 = :(link5.body),
        axis=axisLink, phi = Var(init = getRefPathInitPosition(referencePath1, 1)), w=Var(init=0.0)),

    servo1 = Servo,
    servo2 = Servo,
    servo3 = Servo,
    servo4 = Servo,
    servo5 = Servo,

    equations=:[
        refPath = calculateRobotMovement(getReferencePath(), instantiatedModel),
        servo1.refLoadAngle = getRefPathPosition(refPath, 1),
        servo2.refLoadAngle = getRefPathPosition(refPath, 2),
        servo3.refLoadAngle = getRefPathPosition(refPath, 3),
        servo4.refLoadAngle = getRefPathPosition(refPath, 4),
        servo5.refLoadAngle = getRefPathPosition(refPath, 5)
    ],

    connect = :[
        (servo1.flange, rev1.flange)
        (servo2.flange, rev2.flange)
        (servo3.flange, rev3.flange)
        (servo4.flange, rev4.flange)
        (servo5.flange, rev5.flange)
        ]
)

modelParameters = Map(
    servo1 = servoParameters1,
    servo2 = servoParameters2,
    servo3 = servoParameters3,
    servo4 = servoParameters4,
    servo5 = servoParameters5,
)

youbotModel = buildModia3D(YouBot) | modelParameters

youbot = @instantiateModel(youbotModel, unitless=true, logCode=false, log=false)

stopTime = 5.0
tolerance = 1e-6
requiredFinalStates = [0.5193537637443382, -0.00014620403190263837, -0.3452197601381744, 0.24389499144876015, -2.61707557353202e-5, 9.138175086356761e-8, 1.5711570542895037, -0.00012975080236810933, -26.728944081206805, 0.004096901654603484, -0.00011590306577238007, -9.798531538293224, -9.38328422293869e-8, 9.389543923301749e-8, -1.8339143521277865e-6, 1.8341244885504516e-6, -3.554504768813606e-6, 3.554985652233493e-6, -2.1817977469044572e-6, 2.182160452324344e-6, 8.488376474459234e-9, -8.489788275075988e-9, -0.000716770056313758, 0.22761856966370303, -0.018209520510203175, -0.008528771785063172, 2.952605056036118e-5]
simulate!(youbot, stopTime=stopTime, tolerance=tolerance, log=true, logStates=false, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(youbot, ["free.rot"], figure=1)

end


text:
using Test
using ClimateMachine.Microphysics
using ClimateMachine.MoistThermodynamics

using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, R_v, grav, R_d, molmass_ratio
using CLIMAParameters.Atmos.Microphysics: τ_cond_evap, τ_acnv, q_liq_threshold

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

@testset "RainDropFallSpeed" begin
    # two typical rain drop sizes
    r_small = 0.5 * 1e-3
    r_big = 3.5 * 1e-3

    # example atmospheric conditions
    p_range = [1013.0, 900.0, 800.0, 700.0, 600.0, 500.0] .* 100
    T_range = [20.0, 20.0, 15.0, 10.0, 0.0, -10.0] .+ 273.15
    ρ_range = p_range ./ R_d(param_set) ./ T_range

    # previousely calculated terminal velocity values
    ref_term_vel_small = [4.44, 4.71, 4.96, 5.25, 5.57, 5.99]
    ref_term_vel_big = [11.75, 12.47, 13.11, 13.90, 14.74, 15.85]

    for idx in range(Int(1), stop = Int(6))
        vc = terminal_velocity_single_drop_coeff(param_set, ρ_range[idx])

        term_vel_small = vc .* sqrt(r_small .* grav(param_set))
        term_vel_big = vc .* sqrt(r_big .* grav(param_set))

        @test term_vel_small ≈ ref_term_vel_small[idx] atol = 0.01
        @test term_vel_big ≈ ref_term_vel_big[idx] atol = 0.01
    end
end

@testset "RainFallSpeed" begin

    # eq. 5d in Smolarkiewicz and Grabowski 1996
    # https://doi.org/10.1175/1520-0493(1996)124<0487:TTLSLM>2.0.CO;2
    function terminal_velocity_empir(
        q_rai::FT,
        q_tot::FT,
        ρ::FT,
        ρ_air_ground::FT,
    ) where {FT <: Real}
        rr = q_rai / (1 - q_tot)
        vel = FT(14.34) * ρ_air_ground^FT(0.5) * ρ^-FT(0.3654) * rr^FT(0.1346)
        return vel
    end

    # some example values
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    ρ_air, q_tot, ρ_air_ground = 1.2, 20 * 1e-3, 1.22

    for q_rai in q_rain_range
        @test terminal_velocity(param_set, q_rai, ρ_air) ≈
              terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground) atol =
            0.2 * terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground)
    end
end

@testset "CloudCondEvap" begin

    q_liq_sat = 5e-3
    frac = [0.0, 0.5, 1.0, 1.5]

    for fr in frac
        q_liq = q_liq_sat * fr
        @test conv_q_vap_to_q_liq(
            param_set,
            PhasePartition(0.0, q_liq_sat, 0.0),
            PhasePartition(0.0, q_liq, 0.0),
        ) ≈ (1 - fr) * q_liq_sat / τ_cond_evap(param_set)
    end
end

@testset "RainAutoconversion" begin

    _q_liq_threshold = q_liq_threshold(param_set)
    _τ_acnv = τ_acnv(param_set)
    q_liq_small = 0.5 * _q_liq_threshold
    @test conv_q_liq_to_q_rai_acnv(param_set, q_liq_small) == 0.0

    q_liq_big = 1.5 * _q_liq_threshold
    @test conv_q_liq_to_q_rai_acnv(param_set, q_liq_big) ==
          0.5 * _q_liq_threshold / _τ_acnv
end

@testset "RainAccretion" begin

    # eq. 5b in Smolarkiewicz and Grabowski 1996
    # https://doi.org/10.1175/1520-0493(1996)124<0487:TTLSLM>2.0.CO;2
    function accretion_empir(q_rai::FT, q_liq::FT, q_tot::FT) where {FT <: Real}
        rr = q_rai / (FT(1) - q_tot)
        rl = q_liq / (FT(1) - q_tot)
        return FT(2.2) * rl * rr^FT(7 / 8)
    end

    # some example values
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    ρ_air, q_liq, q_tot = 1.2, 5e-4, 20e-3

    for q_rai in q_rain_range
        @test conv_q_liq_to_q_rai_accr(param_set, q_liq, q_rai, ρ_air) ≈
              accretion_empir(q_rai, q_liq, q_tot) atol =
            0.1 * accretion_empir(q_rai, q_liq, q_tot)
    end
end

@testset "RainEvaporation" begin

    # eq. 5c in Smolarkiewicz and Grabowski 1996
    # https://doi.org/10.1175/1520-0493(1996)124<0487:TTLSLM>2.0.CO;2
    function rain_evap_empir(
        param_set::AbstractParameterSet,
        q_rai::FT,
        q::PhasePartition,
        T::FT,
        p::FT,
        ρ::FT,
    ) where {FT <: Real}

        ts_neq = TemperatureSHumNonEquil(param_set, T, ρ, q)
        q_sat = q_vap_saturation(ts_neq)
        q_vap = q.tot - q.liq
        rr = q_rai / (1 - q.tot)
        rv_sat = q_sat / (1 - q.tot)
        S = q_vap / q_sat - 1

        ag, bg = FT(5.4 * 1e2), FT(2.55 * 1e5)
        G = FT(1) / (ag + bg / p / rv_sat) / ρ

        av, bv = FT(1.6), FT(124.9)
        F =
            av * (ρ / FT(1e3))^FT(0.525) * rr^FT(0.525) +
            bv * (ρ / FT(1e3))^FT(0.7296) * rr^FT(0.7296)

        return 1 / (1 - q.tot) * S * F * G
    end

    # example values
    T, p = 273.15 + 15, 90000.0
    ϵ = 1.0 / molmass_ratio(param_set)
    p_sat = saturation_vapor_pressure(param_set, T, Liquid())
    q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    q_tot = 15e-3
    q_vap = 0.15 * q_sat
    q_ice = 0.0
    q_liq = q_tot - q_vap - q_ice
    q = PhasePartition(q_tot, q_liq, q_ice)
    R = gas_constant_air(param_set, q)
    ρ = p / R / T

    for q_rai in q_rain_range
        @test conv_q_rai_to_q_vap(param_set, q_rai, q, T, p, ρ) ≈
              rain_evap_empir(param_set, q_rai, q, T, p, ρ) atol =
            -0.5 * rain_evap_empir(param_set, q_rai, q, T, p, ρ)
    end
end


text:
# Download the .txt file using curl as follows
# curl -o words.txt  https://raw.githubusercontent.com/BenLauwens/ThinkJulia.jl/master/data/words.txt

# Exercise 9-1
function ex9_1()
    fin = open("words.txt")
    for line in eachline("words.txt")
        if length(line) > 20
            println(line)
        end
    end
end

# ex9_1()

# Exercise 9-2
function hasno_e(word)
    return 'e' ∉ word
end

# println(hasno_e("abcdef"))
# println(hasno_e("abcdf"))

function ex9_2()
    fin = open("words.txt")
    count = 0
    total = 0
    for line in eachline("words.txt")
        if hasno_e(line)
            println(line)
            count += 1
        end
        total += 1
    end
    return count /  total
end

# println(ex9_2())

# Exercise 9-3
function avoids(word, banned)
    for ch in banned
        if ch in word
            return false
        end
    end
    return true
end

# println(avoids("abcdef", "awxy"))
# println(avoids("abcdef", "zwxy"))

function ex9_3()
    println("Enter forbidden characters:")
    banned = readline()

    fin = open("words.txt")
    count = 0
    for line in eachline("words.txt")
        if avoids(line, banned)
            count += 1
        end
    end
    return count
end

# println(ex9_3())

# Exercise 9-4
function useonly(word, letters)
    for ch in word
        if ch == ' '
            continue
        elseif lowercase(ch) ∉ letters
            return false
        end
    end
    return true
end

# println(useonly("Hoe alfalfa?", "acefhlo?"))
# println(useonly("Hoe galfalfa?", "acefhlo?"))

# Exercise 9-5
function useall(word, letters)
    for ch in letters
        if ch ∉ word
            return false
        end
    end
    return true
end

function useall(word, letters)
    return useonly(letters, word)
end

# println(useall("Hoe alfalfa?", "acefhlo?"))
# println(useall("abcdefabcdef", "abcdef"))

function ex9_5()
    check = "aeiou"
    fin = open("words.txt")
    count1, count2 = 0, 0
    for line in eachline("words.txt")
        if useall(line, check)
            count1 += 1
            if 'y' in line
                count2 += 1
            end
        end
    end
    return count1, count2
end

# println(ex9_5())

# Exercise 9-6
function isabecedarian(word)
    index = firstindex(word)
    len = length(word)
    while index < len && index + 1 < len
        if index == 1
            curr = word[index]
        end
        next = word[index + 1]
        if word[index] == word[index + 1] || word[index] + 1 == word[index + 1]
            index += 1
            continue
        else
            return false
        end
    end
    return true
end

# println(isabecedarian("abcdef"))
# println(isabecedarian("abccdef"))
# println(isabecedarian("zabccdef"))

function ex9_6()
    fin = open("words.txt")
    count = 0
    for line in eachline("words.txt")
        if isabecedarian(line)
            count += 1
        end
    end
    return count
end

# println(ex9_6())

# Exercise 9-7
function checkthree_doubles(word)
    count, index = 0, 1
    while index <= length(word) - 1
        if word[index] == word[index + 1]
            count += 1
            if count == 3
                return true
            end
            index += 2
        else
            index += 1 - 2*count
            count = 0
        end
    end
    return false
end

# println(checkthree_doubles("bookkeeper"))
# println(checkthree_doubles("bookkaeeper"))

function ex9_7()
    fin = open("words.txt")
    for line in eachline("words.txt")
        if checkthree_doubles(line)
            println(line)
        end
    end
end

# ex9_7()

# Exercise 9-8
function ispalindrome(number, l, r)
    word = string(number)
    substr = word[l:r]
    return substr == reverse(substr)
end

# println(ispalindrome("abcba", 1, 5))
# println(ispalindrome("zzzzheyayehzzz", 5, 11))
# println(ispalindrome("zzzzheyayehzzz", 3, 12))

function checknumber(n)
    return (
        ispalindrome(n, 3, 6)
        && ispalindrome(n + 1, 2, 6)
        && ispalindrome(n + 2, 2, 5)
        && ispalindrome(n + 3, 1, 6)
    )
end

function ex9_8()
    for i = 100000 : 999996
        if checknumber(i)
            println(i)
        end
    end
end

# ex9_8()


text:
configs = Dict()
for (i, garbage_collection) in enumerate(0:1:30)
    configs[i] = LitteringModel.Config(
        n_nodes=100,
        n0=3,
        k=3,
        agentcount=500,
        initial_trash=100,
        initial_litterers=5,
        garbage_collection=garbage_collection,
        n_iter=100
    )
end


text:
using Base.Meta
using Nabla: boxed_method

function expected_func(sig::Expr, body::Expr)
    return Expr(:function, sig, Expr(:block, LineNumberNode(0), body))
end

@testset "sensitivity" begin
    # # "Test" `Nabla.get_body`. (Not currently unit testing this as it is awkward. Will
    # # change this at some point in the future to be more unit-testable.)
    # let
    #     println()
    #     from_func = Nabla.get_body(:bar, :(Base.bar), :(Tuple{Real}), [:x])
    #     full_expr = Nabla.add_intercept(:bar, :(Base.bar), :(Tuple{Real}))
    #     println(full_expr)
    # end
    # let
    #     println()
    #     from_func = Nabla.get_body(:bar, :(Base.bar), :(Tuple{T} where T<:Real), [:x])
    #     full_expr = Nabla.add_intercept(:bar, :(Base.bar), :(Tuple{T} where T<:Real))
    #     println(full_expr)
    # end
    # let
    #     println()
    #     from_func = Nabla.get_body(:bar, :(Base.bar), :(Tuple{Vararg}), [:x])
    #     full_expr = Nabla.add_intercept(:bar, :(Base.bar), :(Tuple{Vararg}))
    #     println(full_expr)
    # end

    @testset "boxed_method" begin
        let
            from_func = boxed_method(:foo, :(Tuple{Any}), [true], [:x1])
            expected = expected_func(:(foo(x1::Node{<:Any})),
                                     :(Branch(foo, (x1,), getfield(x1, $(quot(:tape))))))
            @test from_func == expected
        end
        let
            from_func = boxed_method(:foo, :(Tuple{T{V}}), [true], [:x1])
            expected = expected_func(:(foo(x1::Node{<:T{V}})),
                                     :(Branch(foo, (x1,), getfield(x1, $(quot(:tape))))))
            @test from_func == expected
        end
        let
            from_func = boxed_method(:foo, :(Tuple{Any, Any}), [true, false], [:x1, :x2])
            expected = expected_func(:(foo(x1::Node{<:Any}, x2::Any)),
                                     :(Branch(foo, (x1, x2), getfield(x1, $(quot(:tape))))))
            @test from_func == expected
        end
        let
            from_func = boxed_method(:foo, :(Tuple{Any, Any}), [true, true], [:x1, :x2])
            expected = expected_func(:(foo(x1::Node{<:Any}, x2::Node{<:Any})),
                                     :(Branch(foo, (x1, x2), getfield(x1, $(quot(:tape))))))
            @test from_func == expected
        end
        let
            from_func = boxed_method(:foo, :(Tuple{Any, Any}), [false, true], [:x1, :x2])
            expected = expected_func(:(foo(x1::Any, x2::Node{<:Any})),
                                     :(Branch(foo, (x1, x2), getfield(x2, $(quot(:tape))))))
            @test from_func == expected
        end
        let
            from_func = boxed_method(:foo, :(Tuple{T} where T), [true], [:x1])
            expected = expected_func(:(foo(x1::Node{<:T}) where T),
                                     :(Branch(foo, (x1,), getfield(x1, $(quot(:tape))))))
            @test from_func == expected
        end
        let
            from_func = boxed_method(:foo, :(Tuple{Any, Any}), [false, true], [:x1, :x2]; a=1, b=2)
            expected = expected_func(:(foo(x1::Any, x2::Node{<:Any}; a=1, b=2)),
                                     :(Branch(foo, (x1, x2), getfield(x2, $(quot(:tape))); a=a, b=b)))
            @test from_func == expected
        end
    end

    # Test `Nabla.branch_expr`.
    import Nabla.Nabla.branch_expr
    let
        from_func = branch_expr(:bar, [true], (Int,), (:x,), :((x,)))
        tape = Expr(:call, :getfield, :x, quot(:tape))
        expected = Expr(:call, :Branch, :bar, :((x,)), tape)
        @test from_func == expected
    end
    let
        from_func = branch_expr(:bar, [true, false], (Int,), (:x, :y), :((x, y)))
        tape = Expr(:call, :getfield, :x, quot(:tape))
        expected = Expr(:call, :Branch, :bar, :((x, y)), tape)
        @test from_func == expected
    end
    let
        from_func = branch_expr(:bar, [false, true], (Int, Int), (:x, :y), :((x, y)))
        tape = Expr(:call, :getfield, :y, quot(:tape))
        expected = Expr(:call, :Branch, :bar, :((x, y)), tape)
        @test from_func == expected
    end
    let
        from_func = branch_expr(:bar, [true, true], (Int, Int), (:x, :y), :((x, y)))
        tape = Expr(:call, :getfield, :x, quot(:tape))
        expected = Expr(:call, :Branch, :bar, :((x, y)), tape)
        @test from_func == expected
    end
    let
        from_func = branch_expr(:bar, [true], ((Leaf{Float64},),), (:x,), :((x...,)))
        tape = Expr(:call, :getfield, Expr(:ref, :x, :1), quot(:tape))
        expected = Expr(:call, :Branch, :bar, :((x...,)), tape)
        @test from_func == expected
    end

    # Test Nabla.invoke_expr
    import Nabla.invoke_expr
    let
        from_func = invoke_expr(:foo, :(Tuple{T, V}), [:x1, :x2])
        expected = Expr(:call, :invoke, :foo, :(Tuple{T, V}), :x1, :x2)
        @test from_func == expected
    end

    # Test `Nabla.get_union_call`.
    import Nabla.Nabla.get_union_call
    let
        from_func = get_union_call(:goo, :(Tuple{Real}))[1]
        expected = :(goo(x1::Union{Real, Node{<:Real}}))
        @test from_func == expected
    end
    let
        from_func = get_union_call(:goo, :(Tuple{Real, Number}))[1]
        expected = :(
        goo(
            x1::Union{Real, Node{<:Real}},
            x2::Union{Number, Node{<:Number}},
        ))
        @test from_func == expected
    end
    let
        arg = :(AbstractArray{V} where V)
        from_func = get_union_call(:goo, :(Tuple{$arg}))[1]
        expected = :(goo(x1::Union{$arg, Node{<:($arg)}}))
        @test from_func == expected
    end
    let
        from_func = get_union_call(:goo, :(Tuple{V, V} where V<:Real))[1]
        expected = :(goo(x1::Union{V, Node{<:V}}, x2::Union{V, Node{<:V}}) where V<:Real)
        @test from_func == expected
    end
    let
        from_func = get_union_call(:goo, :(Tuple{Vararg{Real}}))[1]
        expected = :(goo(x1::Vararg{Union{Real, Node{<:Real}}}))
        @test from_func == expected
    end
    let
        from_func = get_union_call(:goo, :(Tuple{Number, Vararg{Real}}))[1]
        expected = :(
        goo(
            x1::Union{Number, Node{<:Number}},
            x2::Vararg{Union{Real, Node{<:Real}}},
        ))
        @test from_func == expected
    end
    let
        from_func = get_union_call(:goo, :(Tuple{Vararg{V where V<:Real}}))[1]
        expected = :(goo(x1::Vararg{Union{V where V<:Real, Node{<:(V where V<:Real)}}}))
        @test from_func == expected
    end
    let
        arg_type = :(Union{Real, AbstractArray})
        from_func = get_union_call(:goo, :(Tuple{Vararg{$arg_type}}))[1]
        expected = :(goo(x1::Vararg{Union{$arg_type, Node{<:$arg_type}}}))
        @test from_func == expected
    end
    let
        arg_type = :(Union{Real, AbstractArray{V} where V<:Real})
        from_func = get_union_call(:goo, :(Tuple{Vararg{$arg_type}}))[1]
        expected = :(goo(x1::Vararg{Union{$arg_type, Node{<:$arg_type}}}))
        @test from_func == expected
    end
end


text:
# https://github.com/spencerlyon2/PlotlyJS.jl

supported_args(::PlotlyJSBackend) = supported_args(PlotlyBackend())
supported_types(::PlotlyJSBackend) = supported_types(PlotlyBackend())
supported_styles(::PlotlyJSBackend) = supported_styles(PlotlyBackend())
supported_markers(::PlotlyJSBackend) = supported_markers(PlotlyBackend())
supported_scales(::PlotlyJSBackend) = supported_scales(PlotlyBackend())
is_subplot_supported(::PlotlyJSBackend) = true
is_string_supported(::PlotlyJSBackend) = true

# --------------------------------------------------------------------------------------

function _initialize_backend(::PlotlyJSBackend; kw...)
    @eval begin
        import PlotlyJS
        export PlotlyJS
    end

    # # override IJulia inline display
    # if isijulia()
    #     IJulia.display_dict(plt::AbstractPlot{PlotlyJSBackend}) = IJulia.display_dict(plt.o)
    # end
end

# ---------------------------------------------------------------------------

function _create_backend_figure(plt::Plot{PlotlyJSBackend})
    PlotlyJS.plot()
end

function _series_added(plt::Plot{PlotlyJSBackend}, series::Series)
    syncplot = plt.o
    pdict = plotly_series(plt, series)
    typ = pop!(pdict, :type)
    gt = PlotlyJS.GenericTrace(typ; pdict...)
    PlotlyJS.addtraces!(syncplot, gt)
end

function _series_updated(plt::Plot{PlotlyJSBackend}, series::Series)
    xsym, ysym = (ispolar(series) ? (:t,:r) : (:x,:y))
    PlotlyJS.restyle!(
        plt.o,
        findfirst(plt.series_list, series),
        KW(xsym => (series.d[:x],), ysym => (series.d[:y],))
    )
end

# ----------------------------------------------------------------

function _update_plot_object(plt::Plot{PlotlyJSBackend})
    pdict = plotly_layout(plt)
    syncplot = plt.o
    w,h = plt[:size]
    PlotlyJS.relayout!(syncplot, pdict, width = w, height = h)
end

# ----------------------------------------------------------------

function _writemime(io::IO, ::MIME"image/svg+xml", plt::Plot{PlotlyJSBackend})
    writemime(io, MIME("text/html"), plt.o)
end

function _writemime(io::IO, ::MIME"image/png", plt::Plot{PlotlyJSBackend})
    tmpfn = tempname() * ".png"
    PlotlyJS.savefig(plt.o, tmpfn)
    write(io, read(open(tmpfn)))
end

function _display(plt::Plot{PlotlyJSBackend})
    display(plt.o)
end


text:
@symbol_func function K_PC(cur_reactor::AbstractReactor)
  cur_K = K_F(cur_reactor)

  cur_K *= K_n(cur_reactor) ^ 2

  cur_K *= 1 + E_Li / E_F

  cur_K *= cur_reactor.eta_T

  cur_K /= cur_reactor.max_P_E

  cur_K
end