text:
using Swizzles.StylishArrays
using Swizzles.StylishArrays: square, root, power
using LinearAlgebra

@testset "StylishArrays" begin
    @testset "Squares" begin
        a = sqrt(prevfloat(typemax(Float64)))
        x = square(a)
        @test root(x) ≈ a
        @test root(x + x + x + x) ≈ 2 * a
        @test root(square(0) + x) ≈ a
        @test root(x + square(0)) ≈ a
        @test root(square(0) + square(0)) ≈ 0
        y = square(Inf)
        z = square(-Inf)
        @test root(y) == Inf
        @test root(y + z) == Inf
    end

    @testset "Powers" begin
        a = sqrt(prevfloat(typemax(Float64)))
        x = power(a, 2)
        @test root(x) ≈ a
        @test root(x + x + x + x) ≈ 2 * a
        @test root(power(0, 2) + x) ≈ a
        @test root(x + power(0, 2)) ≈ a
        @test root(power(0, 2) + power(0, 2)) ≈ 0
        y = power(Inf, 2)
        z = power(-Inf, 2)
        @test root(y) == Inf
        @test root(y + z) == Inf
    end

    struct TestStylishArray{T, N, Arr<:AbstractArray{T, N}} <: StylishArray{T, N}
        parent::Arr
    end
    Base.size(arr::TestStylishArray) = size(arr.parent)
    Base.axes(arr::TestStylishArray) = axes(arr.parent)
    Base.getindex(arr::TestStylishArray, i...) = arr.parent[i...]
    Base.getindex(arr::TestStylishArray, i::Int) = arr.parent[i]
    Base.getindex(arr::TestStylishArray, i::CartesianIndex) = arr.parent[i]
    Base.setindex!(arr::TestStylishArray, x, i...) = arr.parent[i...] = x
    Base.setindex!(arr::TestStylishArray, x, i::Int) = arr.parent[i] = x
    Base.setindex!(arr::TestStylishArray, x, i::CartesianIndex) = arr.parent[i] = x

    @testset "Base" begin

        @testset "copyto!" begin
            for A in ([1 2 3; 4 5 6; 7 8 9], TestStylishArray([1 2 3; 4 5 6; 7 8 9]))
                B = TestStylishArray(zeros(3,3))
                @test A == copyto!(B, A)
                @test A == B
                B = TestStylishArray(zeros(9))
                #@test reshape(A, :) == copyto!(B, A) #FIXME
                #@test reshape(A, :) == B #FIXME
            end

            A = TestStylishArray([1 2 3; 4 5 6; 7 8 9])
            B = zeros(3,3)
            @test A == copyto!(B, A)
            @test A == B
            B = zeros(9)
            #@test reshape(A, :) == copyto!(B, A) #FIXME
            #@test reshape(A, :) == B #FIXME
        end

        @testset "getindex and setindex!" begin
            A = [1 2 3; 4 5 6; 7 8 9]
            @test TestStylishArray(A)[1] == A[1]
            @test TestStylishArray(A)[1:3] == A[1:3]
            @test TestStylishArray(A)[2:3, 1:2] == A[2:3, 1:2]
            @test TestStylishArray(A)[2, 1:3] == A[2, 1:3]

            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            TestStylishArray(A)[1] = -1
            B[1] = -1
            @test A == B

            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            TestStylishArray(A)[1:3] = [-1 -1 -1]
            B[1:3] = [-1 -1 -1]
            @test A == B

            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            TestStylishArray(A)[2:3, 1:2] = [-1 -1; -1 -1]
            B[2:3, 1:2] = [-1 -1; -1 -1]
            @test A == B

            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            TestStylishArray(A)[2, 1:3] = [-1 -1 -1]
            B[2, 1:3] = [-1 -1 -1]
            @test A == B
        end

        @testset "map" begin
            A = [1 2 3; 4 5 6; 7 8 9]
            @test map(x -> x == 2, TestStylishArray(A)) == map(x -> x == 2, A)

            A = [1 2 3; 4 5 6; 7 8 9]
            @test map((a, b) -> a + b == 4, TestStylishArray(A), TestStylishArray(A)) == map((a, b) -> a + b == 4, A, A)
        end

        @testset "foreach" begin
            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            C = [1 2 3; 4 5 6; 7 8 9]
            @test foreach(i -> B[i] = -i, TestStylishArray(A)) == foreach(i -> C[i] = -i, A)
            @test B == C

            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            C = [1 2 3; 4 5 6; 7 8 9]
            @test foreach((a, b) -> B[a] = b, TestStylishArray(A), -A) == foreach((a, b) -> C[a] = b, A, -A)
            @test B == C

            A = [1 2 3; 4 5 6; 7 8 9]
            B = [1 2 3; 4 5 6; 7 8 9]
            C = [1 2 3; 4 5 6; 7 8 9]
            @test foreach((a, b) -> B[a] = b, TestStylishArray(A), -A) == foreach((a, b) -> C[a] = b, A, -A)
            @test B == C
        end

        @testset "reductions" begin
            A = [1 2 3; 4 5 6; 7 8 9]

            @test mapreduce(-, +, TestStylishArray(A)) == mapreduce(-, +, A)
            @test mapreduce(-, +, TestStylishArray(A), dims=2) == mapreduce(-, +, A, dims=2)
            @test mapreduce(-, +, TestStylishArray(A), dims=(1, 2)) == mapreduce(-, +, A, dims=(1, 2))
            @test mapreduce(-, +, TestStylishArray(A), dims=:) == mapreduce(-, +, A, dims=:)
            @test mapreduce(-, +, TestStylishArray(A), dims=()) == mapreduce(-, +, A, dims=())
            @test mapreduce(-, +, TestStylishArray(A), init=10) == mapreduce(-, +, A, init=10)

            @test reduce(+, TestStylishArray(A)) == reduce(+, A)
            @test reduce(+, TestStylishArray(A), dims=2) == reduce(+, A, dims=2)
            @test reduce(+, TestStylishArray(A), dims=(1, 2)) == reduce(+, A, dims=(1, 2))
            @test reduce(+, TestStylishArray(A), dims=:) == reduce(+, A, dims=:)
            @test reduce(+, TestStylishArray(A), dims=()) == reduce(+, A, dims=())
            @test reduce(+, TestStylishArray(A), init=10) == reduce(+, A, init=10)

            @test sum(TestStylishArray(A)) == sum(A)
            @test sum(TestStylishArray(A), dims=2) == sum(A, dims=2)
            @test sum(TestStylishArray(A), dims=(1, 2)) == sum(A, dims=(1, 2))
            @test sum(TestStylishArray(A), dims=:) == sum(A, dims=:)
            @test sum(TestStylishArray(A), dims=()) == sum(A, dims=())
            @test sum(TestStylishArray(A), init=10) == reduce(+, A, init=10)
        end
    end

    @testset "LinearAlgebra" begin
        @testset "norm" begin
            for (X, Y) in ((rand(3, 4), rand(3,4)), (rand(3), rand(3)))
                @test dot(TestStylishArray(X), TestStylishArray(Y)) ≈ dot(X, Y)
            end

            big = prevfloat(typemax(Float64))/16
            small = nextfloat(zero(Float64))

            for ps in ((-Inf,), (-3.0,), (-1,), (-0.5,), (0,), (0.5,), (1,), (2,), (3.0,), (Inf), ())
                for X in ([1, 2], [big, big], [small, small], rand(3, 3), rand(3))
                    @test norm(TestStylishArray(X), ps...) ≈ norm(X, ps...)
                end
            end
        end

        @testset "mul!" begin
            A = rand(3,3)
            u = rand(3,1)
            v = rand(1,3)
            w = rand(3)

            for (a, b) in Iterators.product((A, u, v, w), (A, u, v, w))
                valid = false
                y = nothing
                try
                    y = a * b
                    valid = true
                catch MethodError
                end
                if valid
                    @test TestStylishArray(a) * b ≈ a * b
                    @test a * TestStylishArray(b) ≈ a * b
                    @test TestStylishArray(a) * TestStylishArray(b) ≈ a * b

                    @test mul!(identity.(y), TestStylishArray(a), b) ≈ mul!(identity.(y), a, b)
                    @test mul!(identity.(y), a, TestStylishArray(b)) ≈ mul!(identity.(y), a, b)
                    @test mul!(identity.(y), TestStylishArray(a), TestStylishArray(b)) ≈ mul!(identity.(y), a, b)
                    @test mul!(TestStylishArray(identity.(y)), a, b) ≈ mul!(identity.(y), a, b)
                    @test mul!(TestStylishArray(identity.(y)), TestStylishArray(a), b) ≈ mul!(identity.(y), a, b)
                    @test mul!(TestStylishArray(identity.(y)), a, TestStylishArray(b)) ≈ mul!(identity.(y), a, b)
                    @test mul!(TestStylishArray(identity.(y)), TestStylishArray(a), TestStylishArray(b)) ≈ mul!(identity.(y), a, b)
                end
            end
        end

        @testset "transpose and adjoint" begin
            A = rand(Complex{Float64}, 3,3)
            u = rand(Complex{Float64}, 3,1)
            v = rand(Complex{Float64}, 1,3)
            w = rand(Complex{Float64}, 3)

            for (a) in (A, u, v, w)
                @test transpose(TestStylishArray(a)) == transpose(a)
                @test transpose(TestStylishArray(a)) isa Swizzles.SwizzledArray
                @test adjoint(TestStylishArray(a)) == adjoint(a)
                @test adjoint(TestStylishArray(a)) isa Swizzles.SwizzledArray
            end
        end

        @testset "permutedims" begin
            A = rand(3,3)
            u = rand(3,1)
            v = rand(1,3)
            w = rand(3)

            for a in (A, u, v, w)
                @test permutedims(TestStylishArray(a)) == permutedims(a)
                @test permutedims(TestStylishArray(a)) isa Swizzles.SwizzledArray
            end

            A = rand(3,3,3)

            for perm in [(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)]
                @test permutedims(TestStylishArray(A), perm) == permutedims(A, perm)
                @test permutedims(TestStylishArray(A), perm) isa Swizzles.SwizzledArray
            end
        end
    end
end


text:
module Smoothers

using Dierckx
import Base: filter

export filter, hma, loess, sma, stl

# Methods
include("filter.jl")
include("hma.jl")
include("loess.jl")
include("sma.jl")
include("stl.jl")

"""
Collection of smoothing heuristics, models and smoothing related applications. The current available smoothers and applications are:

    `filter`: Linear Time-invariant Difference Equation Filter (Matlab/Octave)
    `hma`:    Henderson Moving Average Filter
    `loess`:  Locally Estimated Scatterplot Smoothing
    `sma`:    Simple Moving Average
    `stl`:    Seasonal and Trend decomposition based on Loess

"""
Smoothers

end

ext: jl


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# 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.
using Documenter, Pkg

push!(LOAD_PATH,"../src/")

using Tesseract

test = true
deploy = true

for arg in ARGS
    if arg == "nodeploy"
        global deploy = false
    elseif arg == "notest"
        global test = false
    else
        println("Argument: $arg")
    end
end

println("Testing: $test")
println("Deploying: $deploy")

makedocs(
    sitename="Tesseract.jl",
    doctest=test,
    pages = [
        "index.md",
        "download.md",
        "load.md",
        "single.md",
        "multiple.md",
        "Reference" => [
            "reference.md",
            hide("Tesseract.download\\_languages" => "ref/download_languages.md"),
            hide("Tesseract.download\\_pdf\\_font" => "ref/download_pdf_font.md"),
            hide("Tesseract.IFF" => "ref/iff.md"),
            hide("Tesseract.lept\\_version" => "ref/lept_version.md"),
            hide("Tesseract.Pix" => "ref/pix.md"),
            hide("Tesseract.PixBox" => "ref/pix_box.md"),
            hide("Tesseract.pix\\_delete" => "ref/pix_delete.md"),
            hide("Tesseract.pix\\_get\\_depth" => "ref/pix_get_depth.md"),
            hide("Tesseract.pix\\_get\\_dimensions" => "ref/pix_get_dimensions.md"),
            hide("Tesseract.pix\\_read" => "ref/pix_read.md"),
            hide("Tesseract.pix\\_read\\_bmp" => "ref/pix_read_bmp.md"),
            hide("Tesseract.pix\\_read\\_gif" => "ref/pix_read_gif.md"),
            hide("Tesseract.pix\\_read\\_jp2k" => "ref/pix_read_jp2k.md"),
            hide("Tesseract.pix\\_read\\_jpeg" => "ref/pix_read_jpeg.md"),
            hide("Tesseract.pix\\_read\\_png" => "ref/pix_read_png.md"),
            hide("Tesseract.pix\\_read\\_pnm" => "ref/pix_read_pnm.md"),
            hide("Tesseract.pix\\_read\\_spix" => "ref/pix_read_spix.md"),
            hide("Tesseract.pix\\_read\\_tiff" => "ref/pix_read_tiff.md"),
            hide("Tesseract.pix\\_read\\_webp" => "ref/pix_read_webp.md"),
            hide("Tesseract.pix\\_write" => "ref/pix_write.md"),
            hide("Tesseract.pix\\_write\\_implied\\_format" => "ref/pix_write_implied_format.md"),
            hide("Tesseract.pix\\_write\\_bmp" => "ref/pix_write_bmp.md"),
            hide("Tesseract.pix\\_write\\_gif" => "ref/pix_write_gif.md"),
            hide("Tesseract.pix\\_write\\_jp2k" => "ref/pix_write_jp2k.md"),
            hide("Tesseract.pix\\_write\\_jpeg" => "ref/pix_write_jpeg.md"),
            hide("Tesseract.pix\\_write\\_pam" => "ref/pix_write_pam.md"),
            hide("Tesseract.pix\\_write\\_pdf" => "ref/pix_write_pdf.md"),
            hide("Tesseract.pix\\_write\\_png" => "ref/pix_write_png.md"),
            hide("Tesseract.pix\\_write\\_pnm" => "ref/pix_write_pnm.md"),
            hide("Tesseract.pix\\_write\\_ps" => "ref/pix_write_ps.md"),
            hide("Tesseract.pix\\_write\\_spix" => "ref/pix_write_spix.md"),
            hide("Tesseract.pix\\_write\\_tiff" => "ref/pix_write_tiff.md"),
            hide("Tesseract.pix\\_write\\_webp" => "ref/pix_write_webp.md"),
            hide("Tesseract.sample\\_pix" => "ref/sample_pix.md"),
            hide("Tesseract.sample\\_tiff" => "ref/sample_tiff.md"),
            hide("Tesseract.TessInst" => "ref/tess_inst.md"),
            hide("Tesseract.TessOutput" => "ref/tess_output.md"),
            hide("Tesseract.TessParam" => "ref/tess_param.md"),
            hide("Tesseract.TessPipeline" => "ref/tess_pipeline.md"),
            hide("Tesseract.tess\\_alto" => "ref/tess_alto.md"),
            hide("Tesseract.tess\\_available\\_languages" => "ref/tess_available_languages.md"),
            hide("Tesseract.tess\\_confidences" => "ref/tess_confidences.md"),
            hide("Tesseract.tess\\_delete" => "ref/tess_delete.md"),
            hide("Tesseract.tess\\_get\\_param" => "ref/tess_get_param.md"),
            hide("Tesseract.tess\\_hocr" => "ref/tess_hocr.md"),
            hide("Tesseract.tess\\_image" => "ref/tess_image.md"),
            hide("Tesseract.tess\\_init" => "ref/tess_init.md"),
            hide("Tesseract.tess\\_initialized\\_languages" => "ref/tess_initialized_languages.md"),
            hide("Tesseract.tess\\_loaded\\_languages" => "ref/tess_loaded_languages.md"),
            hide("Tesseract.tess\\_lstm\\_box" => "ref/tess_lstm_box.md"),
            hide("Tesseract.tess\\_params" => "ref/tess_params.md"),
            hide("Tesseract.tess\\_params\\_parsed" => "ref/tess_params_parsed.md"),
            hide("Tesseract.tess\\_parsed\\_tsv" => "ref/tess_parsed_tsv.md"),
            hide("Tesseract.tess\\_pipeline\\_alto" => "ref/tess_pipeline_alto.md"),
            hide("Tesseract.tess\\_pipeline\\_hocr" => "ref/tess_pipeline_hocr.md"),
            hide("Tesseract.tess\\_pipeline\\_lstm\\_box" => "ref/tess_pipeline_lstm_box.md"),
            hide("Tesseract.tess\\_pipeline\\_pdf" => "ref/tess_pipeline_pdf.md"),
            hide("Tesseract.tess\\_pipeline\\_text" => "ref/tess_pipeline_text.md"),
            hide("Tesseract.tess\\_pipeline\\_tsv" => "ref/tess_pipeline_tsv.md"),
            hide("Tesseract.tess\\_pipeline\\_unlv" => "ref/tess_pipeline_unlv.md"),
            hide("Tesseract.tess\\_pipeline\\_unlv\\_latin1" => "ref/tess_pipeline_unlv_latin1.md"),
            hide("Tesseract.tess\\_pipeline\\_word\\_box" => "ref/tess_pipeline_word_box.md"),
            hide("Tesseract.tess\\_run\\_pipeline" => "ref/tess_run_pipeline.md"),
            hide("Tesseract.tess\\_read\\_config" => "ref/tess_read_config.md"),
            hide("Tesseract.tess\\_read\\_debug\\_config" => "ref/tess_read_debug_config.md"),
            hide("Tesseract.tess\\_recognize" => "ref/tess_recognize.md"),
            hide("Tesseract.tess\\_resolution" => "ref/tess_resolution.md"),
            hide("Tesseract.tess\\_set\\_param" => "ref/tess_set_param.md"),
            hide("Tesseract.tess\\_text" => "ref/tess_text.md"),
            hide("Tesseract.tess\\_text\\_box" => "ref/tess_text_box.md"),
            hide("Tesseract.tess\\_tsv" => "ref/tess_tsv.md"),
            hide("Tesseract.tess\\_unlv" => "ref/tess_unlv.md"),
            hide("Tesseract.tess\\_unlv\\_latin1" => "ref/tess_unlv_latin1.md"),
            hide("Tesseract.tess\\_version" => "ref/tess_version.md"),
            hide("Tesseract.tess\\_word\\_box" => "ref/tess_word_box.md"),
            hide("Tesseract.Tsv" => "ref/tsv.md"),
            hide("Tesseract.update\\_languages" => "ref/update_languages.md"),
            hide("Tesseract.update\\_pdf\\_font" => "ref/update_pdf_font.md"),
            hide("Tesseract.DataFiles.GitHubProject" => "ref/git_hub_project.md")
        ]
    ]
)
if deploy == true
    deploydocs(repo = "github.com/pixel27/Tesseract.jl.git")
end


text:
using Plots
using Distributed
using LinearAlgebra
import YAML
using Statistics
# using Pandas

@everywhere using DPMMSubClusters
@everywhere using Random

include("benchmarks/experiment_utils.jl")

CONF_PATH = "/home/vlad/projects/DPMMSubClusters.jl/benchmarks/configs/basic.yaml"

VERBOSE = true

D = 128

# load config:
cfg = YAML.load_file(CONF_PATH)

SEED = cfg["test_conf"]["seed"]
REPEATS = cfg["test_conf"]["repeats"]
BURNOUT = cfg["test_conf"]["burnout"]
ITERS = cfg["test_conf"]["iterations"]
ALPHA = cfg["test_conf"]["alpha"]

ROOT_DIR = cfg["test_conf"]["root_dir"]
TESTDIR_NAME = cfg["test_conf"]["test_dir"]

real_datasets = cfg["data"]["real_datasets"]

# REPEATS = 5
# BURNOUT = 10
# ITERS = 100
# ALPHA = 100

try
    mkdir(ROOT_DIR*TESTDIR_NAME)
catch e
    println("Folder exists")
end
BASE_PATH = ROOT_DIR*TESTDIR_NAME*"datasets/"

try
    mkdir(BASE_PATH)
catch e
    println("Folder exists")
end

# cp(CONF_PATH, ROOT_DIR*TESTDIR_NAME*"/conf.yaml")

try
    cp(CONF_PATH, BASE_PATH*"/conf.yaml")
catch e
    println("Conf exists")
end

dataset_name = "imagenet"
x, labels = load_data(real_datasets[dataset_name])

dec = 20
gt_data = plot_data(x[:, 1:dec:end],labels[1:dec:end])
path = string(BASE_PATH) * "ds_" * string(dataset_name) * "_gt_data.png"
Plots.savefig(gt_data, path)

res_random = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"none", BURNOUT, ITERS, false, VERBOSE)
res_kmeans = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS, false, VERBOSE)
res_splitnet = repeat_experiment_xy(x, labels, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS, false, VERBOSE)

all_results = Dict(
    "random" => res_random,
    "kmeans" => res_kmeans,
    "splitnet" => res_splitnet,
)

plot_and_save_all_ds(all_results, BASE_PATH, dataset_name)

dataset_name = "mnist"
x, labels = load_data(real_datasets[dataset_name])

dec = 20
gt_data = plot_data(x[:, 1:dec:end],labels[1:dec:end])
path = string(BASE_PATH) * "ds_" * string(dataset_name) * "_gt_data.png"
Plots.savefig(gt_data, path)

data_mean = mean(x, dims=2)[:] .* 1.0
data_cov = Float32.(Statistics.cov(x, dims=2) .* 2.0)
10

res_random = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"none", BURNOUT, ITERS, false, VERBOSE)
res_kmeans = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS, false, VERBOSE)
res_splitnet = repeat_experiment_xy(x, labels, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS, false, VERBOSE)

# # ITERS=10
# VERBOSE=true
# REPEATS=3
# hyper_params = DPMMSubClusters.niw_hyperparams(1.0, data_mean, D + 3, data_cov)
# # labels,clusters,weights = DPMMSubClusters.fit(x, float(ALPHA), iters=ITERS , seed = SEED, gt=labels, init_type="random", burnout=BURNOUT)
# # labels,clusters,weights = DPMMSubClusters.fit(x, float(ALPHA), gt=labels, burnout=BURNOUT, iters=ITERS, seed=SEED, init_type="splitnet_128d", verbose=VERBOSE)
# # labels,clusters,weights = DPMMSubClusters.fit(x, float(ALPHA), gt=labels, burnout=BURNOUT, iters=ITERS, seed=SEED, init_type="kmeans", verbose=VERBOSE)
# # labels,clusters,weights = DPMMSubClusters.fit(x, hyper_params, float(ALPHA), gt=labels, burnout=BURNOUT, iters=ITERS, seed=SEED, init_type="random", verbose=VERBOSE)
# # labels,clusters,weights = DPMMSubClusters.fit(x, hyper_params, float(ALPHA), gt=labels, burnout=BURNOUT, iters=ITERS, seed=SEED, init_type="kmeans", verbose=VERBOSE)
# # labels,clusters,weights = DPMMSubClusters.fit(x, hyper_params, float(ALPHA), gt=labels, burnout=BURNOUT, iters=ITERS, seed=SEED, init_type="splitnet_128d", verbose=VERBOSE)

# # res_random = repeat_experiment_xy(x, labels, ALPHA, 2,"none", 5, 10, false, true, hyperparam=hyper_params)
# # res_random = repeat_experiment_xy(x, labels, ALPHA, 1,"none", 5, 10, false, true, hyper_params)
# res_random = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"none", BURNOUT, ITERS, false, VERBOSE, hyperparam=hyper_params)
# res_kmeans = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS, false, VERBOSE, hyperparam=hyper_params)
# res_splitnet = repeat_experiment_xy(x, labels, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS, false, VERBOSE, hyperparam=hyper_params)

all_results = Dict(
    "random" => res_random,
    "kmeans" => res_kmeans,
    "splitnet" => res_splitnet,
)

# nmi_plot = plot_key_for_iters(all_results_hp, "nmi", "NMI vs. Iterations")
# ari_plot = plot_key_for_iters(all_results_hp, "ari", "ARI vs. Iterations")
# ll_plot = plot_key_for_iters(all_results_hp, "ll", "Log-Likelihood vs. Iterations")
# k_plot = plot_key_for_iters(all_results_hp, "k", "Inferred K vs. Iterations")

# nmi_plot = plot_key_for_iters(all_results, "nmi", "NMI vs. Iterations")
# ari_plot = plot_key_for_iters(all_results, "ari", "ARI vs. Iterations")
# ll_plot = plot_key_for_iters(all_results, "ll", "Log-Likelihood vs. Iterations")
# k_plot = plot_key_for_iters(all_results, "k", "Inferred K vs. Iterations")

plot_and_save_all_ds(all_results, BASE_PATH, dataset_name)

dataset_name = "fmnist"
x, labels = load_data(real_datasets[dataset_name])

dec = 20
gt_data = plot_data(x[:, 1:dec:end],labels[1:dec:end])
path = string(BASE_PATH) * "ds_" * string(dataset_name) * "_gt_data.png"
Plots.savefig(gt_data, path)

res_random = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"none", BURNOUT, ITERS, false, VERBOSE)
res_kmeans = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS, false, VERBOSE)
res_splitnet = repeat_experiment_xy(x, labels, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS, false, VERBOSE)

all_results = Dict(
    "random" => res_random,
    "kmeans" => res_kmeans,
    "splitnet" => res_splitnet,
)

# nmi_plot = plot_key_for_iters(all_results, "nmi", "NMI vs. Iterations")
# ari_plot = plot_key_for_iters(all_results, "ari", "ARI vs. Iterations")
# ll_plot = plot_key_for_iters(all_results, "ll", "Log-Likelihood vs. Iterations")
# k_plot = plot_key_for_iters(all_results, "k", "Inferred K vs. Iterations")

plot_and_save_all_ds(all_results, BASE_PATH, dataset_name)

dataset_name = "cifar10"
x, labels = load_data(real_datasets[dataset_name])

dec = 20
gt_data = plot_data(x[:, 1:dec:end],labels[1:dec:end])
path = string(BASE_PATH) * "ds_" * string(dataset_name) * "_gt_data.png"
Plots.savefig(gt_data, path)

res_random = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"none", BURNOUT, ITERS, false, VERBOSE)
res_kmeans = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS, false, VERBOSE)
res_splitnet = repeat_experiment_xy(x, labels, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS, false, VERBOSE)

all_results = Dict(
    "random" => res_random,
    "kmeans" => res_kmeans,
    "splitnet" => res_splitnet,
)

plot_and_save_all_ds(all_results, BASE_PATH, dataset_name)

dataset_name = "cifar100"
x, labels = load_data(real_datasets[dataset_name])

dec = 20
gt_data = plot_data(x[:, 1:dec:end],labels[1:dec:end])
path = string(BASE_PATH) * "ds_" * string(dataset_name) * "_gt_data.png"
Plots.savefig(gt_data, path)

res_random = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"none", BURNOUT, ITERS, false, VERBOSE)
res_kmeans = repeat_experiment_xy(x, labels, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS, false, VERBOSE)
res_splitnet = repeat_experiment_xy(x, labels, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS, false, VERBOSE)

all_results = Dict(
    "random" => res_random,
    "kmeans" => res_kmeans,
    "splitnet" => res_splitnet,
)

plot_and_save_all_ds(all_results, BASE_PATH, dataset_name)
#


text:
println("... testing multithreading  1 ...")
println("number of available threads: ", nthreads())

clk = PClock()
print(clk)
@test clk.id == 0
m = match(r"Clock thread 1 \(\+ (\d+) ac\)", repr(clk))
@test parse(Int, m.captures[1]) > 0
@test length(clk.ac) >= (nthreads() >>> 1)
@test clk.ac[1].thread == 2

println("... parallel clock identification ...")
c1 = pclock(clk, 1)
@test c1.thread == 2
@test c1.clock.id == 1
m = match(r"Active clock 1 on thrd (\d+)\:", repr(c1))
@test parse(Int, m.captures[1]) > 1

println("... remote error handling ...")
put!(clk.ac[1].forth, Simulate.Clear())
err = diag(clk, 1)
@test err[1] isa ErrorException
@test occursin(r"^step\!\(\:\:ActiveClock.+ at threads\.jl", string(err[2][2]))

println("... testing channel and active clock ...")
a = 0
ev = Simulate.DiscreteEvent(Fun(()->global a+=1),Main,1.0,0.0)
put!(clk.ac[1].forth, Simulate.Register(ev))
sleep(0.1)
@test length(c1.clock.sc.events) == 1
cev = Simulate.DiscreteCond(Fun(≥, Fun(tau, c1), 5), Fun(()->global a+=1), Main)
put!(clk.ac[1].forth, Simulate.Register(cev))
sleep(0.1)
@test length(c1.clock.sc.cevents) == 1
b = 0
sp = Simulate.Sample(Fun(()->global b+=1), Main)
put!(clk.ac[1].forth, Simulate.Register(sp))
sleep(0.1)
@test length(c1.clock.sc.samples) == 1

println("... run parallel clock 1 (thread 2) ...")
put!(clk.ac[1].forth, Simulate.Run(10.0))
sleep(0.5)
@test c1.clock.time ≈ 10
@test c1.clock.scount == 1000
@test a == 2
@test b == 1000

println("... testing register(!) five cases ... ")
reset!(c1.clock)
@test c1.clock.time == 0
ev1 = Simulate.DiscreteEvent(Fun(()->global a+=1),Main,1.0,0.0)
ev2 = Simulate.DiscreteEvent(Fun(()->global a+=1),Main,2.0,0.0)
@test Simulate.register(clk, ev1, 0) == 1.0 # 1. register ev1 to clk
@test Simulate.nextevent(clk) == ev1
Simulate.register(clk, ev1, 1)              # 2. register ev1 to 1st parallel clock
@test Simulate.nextevent(c1.clock) == ev1
Simulate.register(c1, ev2, 1)               # 3. register ev2 directly to 1st parallel clock
@test length(c1.clock.sc.events) == 2
Simulate.register(c1, ev2, 0)               # 4. register ev2 back to master
@test length(clk.sc.events) == 2

if nthreads() > 2                           # This fails on CI (only 2 threads)
    Simulate.register(c1, ev2, 2)           # 5. register ev2 to another parallel clock
    c2 = pclock(clk, 2)
    @test Simulate.nextevent(c2.clock) == ev2
end

println("... collapse! ...")
ac = clk.ac
collapse!(clk)
@test all(x->istaskdone(x.ref[]), ac)
@test all(x->x.ch.state == :closed, ac)


text:
using Optim
# using Zygote

const eps = 1e-7

σ(x) = 1 / (1 + exp(-x))
mse(ŷ, y) = sum((ŷ .- y).^2) / length(y)

function misses(predicted::AbstractVecOrMat, data::AbstractVecOrMat)::Float64
  err = 0.0
  n = length(data)
  for i in 1:n
    err += abs2(predicted[i] - data[i]) > eps ? 1 : 0
  end

  return err / n
end

include("main.jl")
include("enc.jl")
include("pre.jl")

function optim()

  function model(input::Array{Float64,1}, param::Array{Float64,1})::Float64
    ϕ = param[1:end-1]
    b = param[end]

    # 1.0, 0.2471
    # 3.1121371444218893, 0.001393518160188087

    M = I
    M *= ControlledRotation(1, 1, Ry, ϕ[1], pad=2)
    M *= ControlledRotation(2, 2, Ry, ϕ[2])
    M *= ControlledRotation(1, 2, Ry, ϕ[3])
    M *= ControlledRotation(2, 1, Ry, ϕ[4])

    r = M * input

    # return abs2(r[2]) + b > 0.5 ? 1 : 0

    ret = abs2(r[2]) + b
    # ret = σ( (ret - 0.5)*1000 )
    return ret
  end

  tr = encode(pre())
  tr_x = [d[1] for d in tr]
  tr_y = [d[2] for d in tr]

  function target(param::Array{Float64,1})::Float64
    prob_y = map( λ->model(λ, param), tr_x)
    return mse(prob_y, tr_y)
  end

  function target2(param::Array{Float64,1})::Float64
    prob_y = map( λ-> σ( (model(λ, param) - 0.5)*1000 ) , tr_x)
    return mse(prob_y, tr_y)
  end

  # p_init = [3.0, 3.0, 3.0, 3.0, 0.0]
  # p_init = [0.0, 0.0, 0.0, 0.0, 0.0]
  p_init = [1.2524733645306474, -1.141071755919144, -0.08147920673002419, -1.2266097709687462, 0.2010126029766885]
  # p_init = [3.323096449858975, 0.045860433297251896, 2.248264115129544, 2.5039894584527858, 0.12339406935001004]

  @show target2(p_init)
  # prob_y = map( λ->model(λ, p_init), tr_x)
  prob_y = map( λ-> σ( (model(λ, p_init) - 0.5)*1000 ) , tr_x)
  predicted = map(x -> x>0.5 ? 1 : 0, prob_y)
  @show misses(predicted, tr_y)

  # grad = target'
  # grad!(stor,x) = copyto!(stor,grad(x))

  manif = Optim.Flat()
  optimizer = Optim.ConjugateGradient(manifold=manif)
  # optimizer = Optim.GradientDescent(manifold=manif)
  # optimizer = Optim.LBFGS(manifold=manif)
  # res = optimize(target, p_init, optimizer; autodiff = :forward)
  # res = optimize(target, grad!, p_init, optimizer)
  res = optimize(target2, p_init, optimizer, Optim.Options(g_tol = 1e-5, allow_f_increases = false, show_trace=false))

  bestmn = Optim.minimum(res)
  bestmnz = Optim.minimizer(res)

  @show res
  @show bestmn
  @show bestmnz

  @show target2(bestmnz)
  prob_y = map( λ-> σ( (model(λ, bestmnz) - 0.5)*1000 ) , tr_x)
  predicted = map(x -> x>0.5 ? 1 : 0, prob_y)
  @show misses(predicted, tr_y)

end

@time optim()

max_stars_repo_name: dandanua/QML

max_issues_repo_name: dandanua/QML

max_forks_repo_name: dandanua/QML


text:
################################################################################
#
#             EllCrv/QQ.jl : Rational elliptic curves
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
#
################################################################################

export integral_model, istorsion_point, laska_kraus_connell, minimal_model,
       order, tates_algorithm_global, tates_algorithm_local, tidy_model,
       torsion_points_division_poly, torsion_points_lutz_nagell,
       torsion_points, torsion_structure

################################################################################
#
#  Order of a point
#
################################################################################

doc"""
***
    order(P::EllCrvPt{fmpq}) -> fmpz

Returns the order of the point $P$ or $0$ if the order is infinite.
"""
function order(P::EllCrvPt{fmpq})
  Q = P
  for i in 1:12
    if !isfinite(Q)
      return fmpz(i)
    end
    Q = Q + P
  end

  return fmpz(0)
end

################################################################################
#
#  Test if point is torsion point
#
################################################################################

doc"""
    istorsion_point(P::EllCrvPt{fmpq}) -> fmpz

> Returns whether the point $P$ is a torsion point.
"""
function istorsion_point(P::EllCrvPt{fmpq})
  o = order(P)
  return o != 0
end

################################################################################
#
#  Torsion points
#
################################################################################

# via theorem of Lutz-Nagell
doc"""
***
    torsion_points_lutz_nagell(E::EllCrv{fmpq}) -> Array{EllCrvPt, 1}

> Computes the rational torsion points of an elliptic curve using the
>Lutz-Nagell theorem.
"""
function torsion_points_lutz_nagell(F::EllCrv{fmpq})

  if F.short == false
    (G, trafo, ruecktrafo) = short_weierstrass_model(F)
  else
    G = F
  end

  # transform the curve to an equivalent one with integer coefficients, if necessary
  (E, trafo_int, trafo_rat) = integral_model(G)

  #!E.short && error("Elliptic curve must be given by short form")
  #if (denominator(E.coeff[1]) != 1) || (denominator(E.coeff[2]) != 1)
  #  error("Need integer coefficients")
  #end

  res = [infinity(E)]
  d = disc(E)

  # Lutz-Nagell: necessary: y = 0 or y^2 divides d

  ycand = squaredivisors(numerator(d)) # candidates for y-coordinate

  push!(ycand,0)

  pcand = Tuple{fmpz, fmpz}[] # candidates for torsion points

  # Lutz-Nagell: coordinates of torsion points need to be in ZZ
  for i = 1:length(ycand)
    # are there corresponding integer x-values?
    xcand = zeros(_Zx(_x)^3 + (numerator(E.coeff[1]))*_Zx(_x) + (numerator(E.coeff[2])) - ycand[i]^2)
    if length(xcand) != 0
      for j = 1: length(xcand)
        push!(pcand, (xcand[j], ycand[i])) # add to candidates
      end
    end
  end

  # check if candidates are torsion points
  if length(pcand) != 0
    for i = 1:length(pcand)
      P = E([pcand[i][1],pcand[i][2]])

      # Mazur: Order of torsion point is at most 12
      Q = P
      for j = 2:12
        Q = Q + P
        if isinfinite(Q)
          push!(res, P)
          break
        end
      end
    end
  end

  torsionpoints = res

  if F.short == false
    for i = 1:length(torsionpoints)
      torsionpoints[i]= ruecktrafo(trafo_rat(torsionpoints[i]))
    end
  else
    for i = 1:length(torsionpoints)
      torsionpoints[i] = trafo_rat(torsionpoints[i])
    end
  end
  return torsionpoints
end

# via division polynomials
doc"""
***
    torsion_points_division_poly(E::EllCrv{fmpq}) -> Array{EllCrvPt}

> Computes the rational torsion points of an rational elliptic $E$ curve using
> division polynomials.
"""
function torsion_points_division_poly(F::EllCrv{fmpq})

  # if necessary, transform curve into short form
  if F.short == false
    (G, trafo, ruecktrafo) = short_weierstrass_model(F)
  else
    G = F
  end

  # transform the curve to an equivalent one with integer coefficients, if necessary
  (E, trafo_int, trafo_rat) = integral_model(G)

  # check whether we already know the torsion points or not
  if isdefined(E, :torsion_points)
    if F.short == false
      F.torsion_points = EllCrvPt{fmpq}[]
      for i = 1:length(E.torsion_points)
        push!(F.torsion_points, ruecktrafo(trafo_rat(E.torsion_points[i])))
      end
    else
      F.torsion_points = EllCrvPt{fmpq}[]
      for i = 1:length(E.torsion_points)
        push!(F.torsion_points, trafo_rat(E.torsion_points[i]))
      end
    end
    return F.torsion_points::Array{EllCrvPt{fmpq}, 1}
  end

  # curve has integer coefficients
  A = numerator(E.coeff[1])
  B = numerator(E.coeff[2])

  torsionpoints = [infinity(E)]

  # points of order 2 (point has order 2 iff y-coordinate is zero)
  # (note: these points are not detected by the division polynomials)

  s = zeros(_Zx(_x)^3 + A*_Zx(_x) + B) # solutions of x^3 + Ax + B = 0
  if length(s) != 0
    for i = 1:length(s)
      P = E([s[i], 0])
      push!(torsionpoints, P)
    end
  end

  # Mazur: order of torsion point is at most 12
  for n = 7 : 12 # points of smaller order will also be found
    Psi = Psi_polynomial(E,n)
    z = zeros(Psi)

    if length(z) == 0
      continue
    end

    # z contains candidates for x-coordinates of torsion points
    for i = 1:length(z)

      # are there corresponding integer-square y-coordinates?
      ysquarecand = z[i]^3 + A*z[i] + B

      if ysquarecand < 0
        continue
      end

      a = isqrtrem(ysquarecand)

      if a[1] == 0 # then ysquarecand = 0, so have found torsion point
        P = E([z[i], ysquarecand])

        # if not already contained in torsionpoints, add P
        if !(P in torsionpoints)
          push!(torsionpoints, P)
        end

        continue
      end

      # now ysquarecand > 0
      if a[2] == 0 # then y is a square of an integer
        P = E([z[i], a[1]])
        Q = E([z[i], -a[1]])

        # if not already contained in torsionpoints, add P and Q
        if !(P in torsionpoints)
          push!(torsionpoints, P)
          push!(torsionpoints, Q)
        end
      end
    end
  end

  if F.short == false
    for i = 1:length(torsionpoints)
      torsionpoints[i]= ruecktrafo(trafo_rat(torsionpoints[i]))
    end
  else
    for i = 1:length(torsionpoints)
      torsionpoints[i] = trafo_rat(torsionpoints[i])
    end
  end

  return torsionpoints::Array{EllCrvPt{fmpq}, 1}
end

# function for users
doc"""
***
    torsion_points(E::EllCrv{fmpq}) -> Array{EllCrvPt{fmpq}, 1}

> Returns the rational torsion points of $E$.
"""
function torsion_points(E::EllCrv{fmpq})
  if isdefined(E, :torsion_points)
    return E.torsion_points::Array{EllCrvPt{fmpq}, 1}
  end

  t = torsion_points_division_poly(E::EllCrv{fmpq})
  E.torsion_points = t
  return t
end

################################################################################
#
#  Torsion structure
#
################################################################################

doc"""
***
    torsion_structure(E::EllCrv{fmpq}) -> (A::Array{fmpz, 1},
                                           B::Array{EllCrvPt{fmpq}, 1}

> Computes the structure of the rational torsion group of an elliptic curve $E$.
> Then `A` is an array with `A = [n]` resp. `A = [n,m]` such that the
> rational torsion of $E$ is isomorphic to $\mathbf Z/n\mathbf Z$ resp.
> $\mathbf Z/n\mathbf Z \times \mathbf Z/m\mathbf Z$.
> And `B` is an array of points with `B = [P]` and $P$ has order $n$ resp.
> `B = [P, Q]` and $P$ has order $n$, $Q$ has order $m$.
"""
function torsion_structure(E::EllCrv{fmpq})
  T = torsion_points(E)
  grouporder = length(T)
  orders = fmpz[]

  for i in 1:grouporder
    push!(orders, 0)
  end

 # determine orders of group elements
  for i in 1:grouporder
    j = 1
    while (j < 13) && (orders[i] == 0)
      if (j*T[i]).isinfinite == true
        orders[i] = j
      end
    j = j + 1
    end
  end

  # find generators
  if in(grouporder, orders) == true # is the group cyclic?
    k = findfirst(orders, grouporder)
    return (fmpz[grouporder], [T[k]])
  else # group not cyclic
    m = div(grouporder, 2)
    k1 = findfirst(orders, 2)
    k2 = findlast(orders, m) # findlast to get different points if m = 2
    points = [T[k1], T[k2]]
  return (fmpz[2, m], points)
  end
end

################################################################################
#
#  Changing the model
#
################################################################################

doc"""
***
integral_model(E::EllCrv{fmpq}) -> (F::EllCrv{fmpz}, function, function)

> Given an elliptic curve $E$ over $\mathbf Q$ in short form, returns an
> isomorphic curve $F$ with model over $\mathbf Z$. The second and third
> return values are the isomorpisms $E \to F$ and $F \to E$.
"""
function integral_model(E::EllCrv{fmpq})
  A = E.coeff[1]
  B = E.coeff[2]

  if (denominator(A) == 1) && (denominator(B) == 1) # curve already has integer coefficients

    trafo_int = function(P)
     return P
    end

    trafo_rat = function(R)
      return R
    end

    return E, trafo_int, trafo_rat # copy(E)?

  end

  mue = lcm(denominator(A), denominator(B))
  Anew = mue^4 * A
  Bnew = mue^6 * B
  E_int = EllipticCurve([Anew, Bnew])

  trafo_int = function(P) # transformes a point on E into a point on E_int
    if P.isinfinite
      return infinity(E_int)
    end

    xnew = mue^2 * P.coordx
    ynew = mue^3 * P.coordy
    Q = E_int([xnew, ynew])
    return Q
  end

  trafo_rat = function(R) # transformes a point on E_int back into a point on E
    if R.isinfinite
      return infinity(E)
    end

    xnew = divexact(R.coordx, mue^2)
    ynew = divexact(R.coordy, mue^3)
    S = E([xnew, ynew])
    return S
  end

  return E_int::EllCrv{fmpq}, trafo_int, trafo_rat
end

################################################################################
#
#  Algorithm of Laska-Kraus-Connel
#
################################################################################

doc"""
***
    laska_kraus_connell(E::EllCrv{fmpz}) -> Array{Nemo.fmpz}

> Given an elliptic curve over $\mathbf Q$ with integral model, this returns an
> isomorphic elliptic curve over $\mathbf Q$ with minimal discriminant.
"""
# algorithm of Laska-Kraus-Connell
function laska_kraus_connell(E::EllCrv{fmpq})
  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  delta = divexact(c4^3 - c6^2, 1728)

  u = fmpz(1)
  g = gcd(c6^2, delta)

  fac = factor(g)

  for (p, ord) in fac
    d = div(ord, 12)

    if p == 2
      a = divexact(c4, 2^(4*d))
      a = mod(a, 16)
      if a > 8
        a = a - 16
      end

      b = divexact(c6, 2^(6*d))
      b = mod(b, 32)
      if b > 16
        b = b - 32
      end

      if (mod(b, 4) != 3) && !((a == 0) && ((b == 0) || (b == 8)))
        d = d - 1
      end

    elseif p == 3
      ord1 = valuation(c6, 3)

      if (ord1 == 6*d + 2)
        d = d - 1
      end
    end
    u = u * p^d
  end

  c4 = divexact(c4, u^4)
  c6 = divexact(c6, u^6)

  b2 = mod(-c6, 12)
  if b2 > 6
      b2 = b2 - 12
  end

  b4 = divexact(b2^2 - c4, 24)
  b6 = divexact(-b2^3 + 36*b2*b4 - c6, 216)

  na1 = mod(b2, 2)
  na2 = divexact(b2 - na1, 4)
  na3 = mod(b6, 2)
  na4 = divexact(b4 - na1*na3, 2)
  na6 = divexact(b6 - na3, 4)

  return EllipticCurve([na1, na2, na3, na4, na6])::EllCrv{fmpq}
end

################################################################################
#
# T(r, s, t, u) transformation
#
################################################################################

# transformation T(r,s,t,u) as in cremona's book
function transform_rstu(E::EllCrv{fmpq}, T::Array{S, 1}) where S
  r = T[1]
  s = T[2]
  t = T[3]
  u = T[4]

  a1 = E.coeff[1]
  a2 = E.coeff[2]
  a3 = E.coeff[3]
  a4 = E.coeff[4]
  a6 = E.coeff[5]

  a1neu = divexact(a1 + 2*s, u)
  a2neu = divexact(a2 - s*a1 + 3*r - s^2, u^2)
  a3neu = divexact(a3 + r*a1 + 2*t, u^3)
  a4neu = divexact(a4 - s*a3 + 2*r*a2 - (t + r*s)*a1 + 3*r^2 - 2*s*t, u^4)
  a6neu = divexact(a6 + r*a4 + r^2*a2 + r^3 - t*a3 - t^2 - r*t*a1, u^6)

  F = EllipticCurve([a1neu, a2neu, a3neu, a4neu, a6neu])

  trafo1 = function(P) # transformes a point on E into a point on F
    if !isfinite(P)
      return infinity(F)
    end

    xnew = divexact(1, u^2)*(P.coordx - r)
    ynew = divexact(1, u^3)*(P.coordy - s*u^2*xnew - t)
    Q = F([xnew, ynew])
    return Q
  end

  trafo2 = function(P) # transformes a point on F into a point on E
    if !isfinite(P)
      return infinity(E)
    end

    xnew = u^2*P.coordx + r
    ynew = u^3*P.coordy + s*u^2*P.coordx + t
    Q = E([xnew, ynew])
    return Q
  end

  return F, trafo1, trafo2
end

################################################################################
#
#  Tates algorithm
#
################################################################################

doc"""
***
    tates_algorithm_local(E::EllCrv{fmpz}, p::Int) -> EllipticCurve{fmpq}, String, fmpz, fmpz

> Returns a tuple $(\tilde E, K, f, c)$, where $\tilde E$ is a minimal model
> for $E$ at the prime $p$, $K$ is the Kodaira symbol, $f$ is the conduator
> valuation at $p$ and $c$ is the local Tamagawa number at $p$.
"""
# tate's algorithm, see Cremona, p. 66
function tates_algorithm_local(E::EllCrv{fmpq}, p)

  p = FlintZZ(p)

  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  delta = disc(E)
  delta = numerator(delta)

  n = valuation(delta, p)

  # test for type I0
  if n == 0
    return (E, "I0", FlintZZ(0), FlintZZ(1))::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
  end

  # change coordinates so that p | a3, a4, a6
  if p == 2
    if mod(b2, p) == 0
      r = smod(a4, p)
      t = smod(r*(1 + a2 + a4) + a6, p)
    else
      r = smod(a3, p)
      t = smod(r + a4, p)
    end

  elseif p == 3
    if mod(b2, p) == 0
      r = smod(-b6, p)
    else
      r = smod(-b2*b4, p)
    end

    t = smod(a1*r + a3, p)
  else
    if mod(c4, p) == 0
      r = - invmod(FlintZZ(12), p)*b2
    else
      r = - invmod(FlintZZ(12)*c4, p)*(c6 + b2*c4)
    end
      t = - invmod(FlintZZ(2), p)* (a1*r + a3)
      r = smod(r, p)
      t = smod(t, p)
  end

  trans = transform_rstu(E, [r, 0, t, 1])
  E = trans[1]

  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  # test for types In, II, III, IV
  if mod(c4, p) != 0
    if quadroots(1, a1, -a2, p)
      cp = FlintZZ(n)
    elseif mod(n, 2) == 0
      cp = FlintZZ(2)
    else
      cp = FlintZZ(1)
  end

    Kp = "I$(n)"
    fp = FlintZZ(1)

    return (E, Kp, fp, cp)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
  end

  if mod(a6, p^2) != 0
    Kp = "II"
    fp = FlintZZ(n)
    cp = FlintZZ(1)
    return (E, Kp, fp, cp)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
  end

  if mod(b8, p^3) != 0
    Kp = "III"
    fp = FlintZZ(n-1)
    cp = FlintZZ(2)
    return (E, Kp, fp, cp)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
  end

  if mod(b6, p^3) != 0
    if quadroots(1, divexact(a3, p), divexact(-a6, p^2), p)
      cp = FlintZZ(3)
    else
      cp = FlintZZ(1)
    end
    Kp = "IV"
    fp = n - 2
    return (E, Kp, FlintZZ(fp), cp)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
  end

  # change coordinates so that p | a1, a2; p^2 | a3, a4; p^3 | a6
  if p == 2
    s = smod(a2, 2)
    t = 2 * smod(divexact(a6, 4), 2)
  else
    s = -a1 * invmod(FlintZZ(2), p)
    t = -a3 * invmod(FlintZZ(2), p)
  end

  trans = transform_rstu(E, [0, s, t, 1])
  E = trans[1]

  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

  # set up auxililary cubic T^3 + bT^2 + cT + d
  b = divexact(a2, p)
  c = divexact(a4, p^2)
  d = divexact(a6, p^3)
  w = 27*d^2 - b^2*c^2 + 4*b^3*d - 18*b*c*d + 4*c^3
  x = 3*c - b^2

  # test for distinct roots: type I0*
  if mod(w, p) != 0
    Kp = "I0*"
    fp = FlintZZ(n - 4)
    cp = 1 + nrootscubic(b, c, d, p)
    return (E, Kp, fp, FlintZZ(cp))::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}

  # test for double root: type Im*
  elseif mod(x, p) != 0

    # change coordinates so that the double root is T = 0
    if p == 2
      r = c
    elseif p == 3
      r = b*c
    else
      r = (b*c - 9*d) * invmod(FlintZZ(2)*x, p)
    end

    r = p * smod(r, p)

    trans = transform_rstu(E, [r, 0, 0, 1])
    E = trans[1]

    a1 = numerator(E.coeff[1])
    a2 = numerator(E.coeff[2])
    a3 = numerator(E.coeff[3])
    a4 = numerator(E.coeff[4])
    a6 = numerator(E.coeff[5])

    b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

    # make a3, a4, a6 repeatedly more divisible by p
    m = 1
    mx = p^2
    my = p^2
    cp = FlintZZ(0)

    while cp == 0
      xa2 = divexact(a2, p)
      xa3 = divexact(a3, my)
      xa4 = divexact(a4, p*mx)
      xa6 = divexact(a6, mx*my)

      if mod(xa3^2 + 4*xa6, p) !=  0
        if quadroots(1, xa3, -xa6, p)
          cp = FlintZZ(4)
        else
          cp = FlintZZ(2)
        end

      else
        if p == 2
          t = my * xa6
        else
          t = my * smod(-xa3*invmod(FlintZZ(2), p), p)
        end

        trans = transform_rstu(E, [0, 0, t, 1])
        E = trans[1]

        a1 = numerator(E.coeff[1])
        a2 = numerator(E.coeff[2])
        a3 = numerator(E.coeff[3])
        a4 = numerator(E.coeff[4])
        a6 = numerator(E.coeff[5])

        b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

        my = my*p
        m = m + 1
        xa2 = divexact(a2, p)
        xa3 = divexact(a3, my)
        xa4 = divexact(a4, p*mx)
        xa6 = divexact(a6, mx*my)

        if mod(xa4^2 - 4*xa2*xa6, p) != 0
          if quadroots(xa2, xa4, xa6, p)
            cp = FlintZZ(4)
          else
            cp = FlintZZ(2)
          end

        else
          if p == 2
            r = mx * smod(xa6*xa2, 2)
          else
            r = mx * smod(-xa4 * invmod(2*xa2, p), p)
          end

          trans = transform_rstu(E, [r, 0, 0, 1])
          E = trans[1]

          a1 = numerator(E.coeff[1])
          a2 = numerator(E.coeff[2])
          a3 = numerator(E.coeff[3])
          a4 = numerator(E.coeff[4])
          a6 = numerator(E.coeff[5])

          b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

          mx = mx*p
          m = m + 1
        end
      end
    end

    fp = n - m - 4
    Kp = "I$(m)*"

    return (E, Kp, FlintZZ(fp), FlintZZ(cp))::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}

  else
    # Triple root case: types II*, III*, IV* or non-minimal
    # change coordinates so that the triple root is T = 0
    if p == 3
      rp = -d
    else
      rp = -b * invmod(FlintZZ(3), p)
    end

    r = p * smod(rp, p)

    trans = transform_rstu(E, [r, 0, 0, 1])
    E = trans[1]

    a1 = numerator(E.coeff[1])
    a2 = numerator(E.coeff[2])
    a3 = numerator(E.coeff[3])
    a4 = numerator(E.coeff[4])
    a6 = numerator(E.coeff[5])

    b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

    x3 = divexact(a3, p^2)
    x6 = divexact(a6, p^4)

    # Test for type IV*
    if mod(x3^2 + 4* x6, p) != 0
      if quadroots(1, x3, -x6, p)
        cp = FlintZZ(3)
      else
        cp = FlintZZ(1)
      end
      Kp = "IV*"
      fp = FlintZZ(n - 6)

      return (E, Kp, fp, FlintZZ(cp))::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
    else
      if p == 2
        t = x6
      else
        t = x3 * invmod(FlintZZ(2), p)
      end

      t = -p^2 * smod(t, p)

      trans = transform_rstu(E, [0, 0, t, 1])
      E = trans[1]

      a1 = numerator(E.coeff[1])
      a2 = numerator(E.coeff[2])
      a3 = numerator(E.coeff[3])
      a4 = numerator(E.coeff[4])
      a6 = numerator(E.coeff[5])

      b2, b4, b6, b8, c4, c6 = get_b_c_integral(E)

      # Test for types III*, II*
      if mod(a4, p^4) != 0
        Kp = "III*"
        fp = FlintZZ(n - 7)
        cp = FlintZZ(2)

        return (E, Kp, fp, FlintZZ(cp))::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
      elseif mod(a6, p^6) != 0
        Kp = "II*"
        fp = FlintZZ(n - 8)
        cp = FlintZZ(1)

        return (E, Kp, fp, FlintZZ(cp))::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
      else
        E = transform_rstu(E, [0, 0, 0, p])[1]
        return tates_algorithm_local(E, p)::Tuple{EllCrv{fmpq}, String, fmpz, fmpz}
      end
    end
  end
end

doc"""
***
    tates_algorithm_global(E::EllCrv{fmpq}) -> EllCrv{fmpq}

> Returns a global reduced minimal model for $E$ using Tate's algorithm.
"""
function tates_algorithm_global(E)
  delta = abs(numerator(disc(E)))
  fac = factor(delta)

  p_list = [i[1] for i in fac]
  p_list = sort(p_list)

  output = []

  # apply tates algorithm successively for all primes dividing the discriminant
  for p in p_list
    E = tates_algorithm_local(E, p)[1]
  end

  # reduce coefficients (see tidy_model)
  E = tidy_model(E)

  return E::EllCrv{fmpq}
end

################################################################################
#
#  Minimal model(s)
#
################################################################################

doc"""
***
    minimal_model(E::EllCrv{fmpq}) -> EllCrv{fmpq}

> Returns the reduced global minimal model of $E$.
"""
function minimal_model(E::EllCrv{fmpq})
  return tates_algorithm_global(E)
end

doc"""
***
    minimal_model(E::EllCrv{fmpq}) -> EllCrv{fmpq}

> Returns a model of $E$, which is minimal at $p$. It is assumed that $p$
> is prime.
"""
function minimal_model(E::EllCrv{fmpq}, p::Int)
  Ep = tates_algorithm_local(E, p)
  return Ep
end

doc"""
***
    tidy_model(E::EllCrv{fmpz}) -> EllCrv{fmpz}
> Given an elliptic curve with minimal model, this functions returns an
> isomorphic curve with reduced minimal model, that is, $a_1, a_3 \in \{0, 1\}$
> and $a_2 \in \{-1,0,1\}$.
"""
function tidy_model(E::EllCrv{fmpq})

  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  if mod(a1, 2) == 0
    s = -divexact(a1, 2)
  else
    s = -divexact(a1 - 1, 2)
  end

  if mod(-a2 + s*a1 + s^2, 3) == 0
    r = divexact(-a2 + s*a1 + s^2, 3)
  elseif mod(-a2 + s*a1 + s^2, 3) == 1
    r = divexact(-1 - a2 + s*a1 + s^2, 3)
  else
    r = divexact(1 - a2 + s*a1 + s^2, 3)
  end

  if mod(-a3 - r*a1, 2) == 0
    t = divexact(-a3 - r*a1, 2)
  else
    t = divexact(1 - a3 - r*a1, 2)
  end

  E = transform_rstu(E, [r, s, t, 1])[1]

  return E
end

################################################################################
#
#  Integral invariants
#
################################################################################

# this needs to be rewritten
doc"""
    get_b_integral(E::EllCrv{fmpz}) -> Nemo.fmpz
> Computes the invariants b2, b4, b6, b8 of an elliptic curve E with integer coefficients.
"""
function get_b_integral(E)
  a1 = numerator(E.coeff[1])
  a2 = numerator(E.coeff[2])
  a3 = numerator(E.coeff[3])
  a4 = numerator(E.coeff[4])
  a6 = numerator(E.coeff[5])

  b2 = a1^2 + 4*a2
  b4 = a1*a3 + 2*a4
  b6 = a3^2 + 4*a6
  b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

  return b2,b4,b6,b8
end

doc"""
    get_b_c_integral(E::EllCrv{fmpz}) -> Nemo.fmpz
> Computes the invariants b2, b4, b6, b8, c4, c6 of an elliptic curve E with integer coefficients.
"""
function get_b_c_integral(E)
    a1 = numerator(E.coeff[1])
    a2 = numerator(E.coeff[2])
    a3 = numerator(E.coeff[3])
    a4 = numerator(E.coeff[4])
    a6 = numerator(E.coeff[5])

    b2 = a1^2 + 4*a2
    b4 = a1*a3 + 2*a4
    b6 = a3^2 + 4*a6
    b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

    c4 = b2^2 - 24*b4
    c6 = -b2^3 + 36*b2*b4 - 216*b6

    return b2,b4,b6,b8,c4,c6
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e5e0a0da-d45c-11ea-1042-e9b5d0654d4f
md"Fix the value of `c` below to make it `c = a * b`"

# ╔═╡ 4dff4b5e-d461-11ea-29c8-d548fdb5f08b
md"Edit the equation below to calculate the number of pizzas to order using the variables above for **people**, **avg**, and **slices**:"

# ╔═╡ f907e46a-d471-11ea-07e5-f30e2aab3d08
md"""The diameter of a pizza is often stated on a menu so let's define a **formula** to calculate the area of a pizza given the diameter **d**.

We do this by writing a formula like this: `area(d) = pi * (d/2)^2`

Let's write that below:
"""

# ╔═╡ d9575e9c-d472-11ea-1eda-2d335d039f28
md"""Now we have a function called **area** that we can pass any diameter and it will return the area of a pizza (or circle), let's try that with the pizza from before with `area(2*r)` to get the area of the pizza:
"""

# ╔═╡ edb95b14-d473-11ea-3a5a-77382d31f941
md"""## Finding the best pizza deal

Let's see if a larger pizza is a better value by calculating the price per area.  There are 4 sizes: small, medium, large, extra large with the following prices:

Size     | Diameter (inches) | Price ($)
:------- | :---------------: | --------:
small    | 9  | 13.10
medium   | 13 | 20.95
large    | 15 | 24.90
XL       | 17 | 30.95

### 1. How many small pizzas is the same as one XL pizza?

Edit the expression below:
"""

# ╔═╡ 5b07b8fe-d475-11ea-01aa-6b88d6ed8a05
md"""### 2. Calculate the cost per area of each pizza:
"""

# ╔═╡ a42e4eb0-d474-11ea-316a-3d864451bc01
md"Which size of pizza is the best deal?  Write your answer below and assign it to the variable **best_value**."

# ╔═╡ cb419286-d4ff-11ea-1d7f-af5c8574b775
md"""### 3. Is this a good deal?

San Marinos has a special **\"Buy two medium pizzas and save \$5\"**.  Is this a better deal than buying a extra-large pizza?

Calculate the total cost of two medium pizzas deal (saving \$5):"""

# ╔═╡ 0d76d97c-d500-11ea-2433-e96c6fc43b05
md"Calculate the total area of two medium pizzas:"

# ╔═╡ 20a1e9cc-d500-11ea-3d9b-279c71bc20f1
md"Now calculate cost per area by taking the total cost of two medium pizzas and divide by the total area:"

# ╔═╡ 57f024ae-d500-11ea-1cc4-ed28348fdf93
md"""Is it a better deal to get two medium pizzas for \$5 off or to just buy an extra-large?"""

# ╔═╡ 180c8fdc-d503-11ea-04ca-bf2c07fd1c17
md"""### 4. Advanced Problem

A new worker at a pizza shop was getting paid for cutting pizza into pieces.  The pieces of pizza could be any size.  Calculate the maximum number of pieces the worker could make with two cuts of the pizza."""

# ╔═╡ 92b4a012-d503-11ea-15a2-1f3a446d3284
md"Now what about 3 cuts across the pizza?  What is the maximum number of pieces that can be made with **3 cuts**?"

# ╔═╡ 2eb9a560-d507-11ea-3b8b-9d06678fe131
md"Now, how many pieces can be made with **4 cuts**?"

# ╔═╡ d1e3dec0-d507-11ea-1213-d37a9325ee2f
md"Are you starting to see a pattern?  Can you figure out a formula for how many pieces of pizza can be made with \"n\" cuts?  Make a table and fill in the number of pieces for a number of cuts and see if you can find the pattern:

Cuts | Pieces
:--- | ------:
0    |   1
1    |   2
2    |   4
3    |
4    |
"

# ╔═╡ 97bfd13c-dcc2-11ea-0067-ad8c2c6517fc
md"To get an extra hint, figure out how many slices we can get from **5 cuts**:"

# ╔═╡ e0cb2822-dcc2-11ea-2c85-5748bfe526dc
md"Have you found the pattern? Write down the formula below:"

# ╔═╡ 03249876-d508-11ea-16bb-fd5afed37a1f
md"""##### Let's test your formula!"""

# ╔═╡ 14158eb0-d45c-11ea-088f-330e45412320
a = 2

# ╔═╡ 2ed4bb92-d45c-11ea-0b31-2d8e32ce7b44
b = 6

# ╔═╡ 03664f5c-d45c-11ea-21b6-91cd647a07aa
md"# Mathematics in Julia 🍕
This is an introduction to programming.  Let's get started!

## Let's make a calculator!

First let's do some simple math with setting **a = $a**, **b = $b** and **c = a * b**.  What will **c** equal?

Type in the cells (with the coloured background) below and press **`Shift-Enter`** or the click the right-arrow button (▶) to the right to execute the cell after changing the values."

# ╔═╡ 30f0f882-d45c-11ea-2adc-7d84ecf8a7a6
c = 10

# ╔═╡ 262b312a-d460-11ea-26c5-df30459effc5
people = 10

# ╔═╡ 2ea7f162-d460-11ea-0e8e-25340e2e64da
avg = 2.5

# ╔═╡ 3da812c6-d460-11ea-0170-79fbb6a4347c
slices = 8

# ╔═╡ a38cb92e-d45e-11ea-2959-05be909befb2
md"""### Now you have a calculator!

You did multiplication above.  Here's how you do other mathematical operations:

Operation | Type This
:------------ | :-------------:
add | +
subtract | -
multiply | *
divide | /
power | ^

### Pizza Slices

Let's try this out on a problem.  Let's say you want to order pizzas for $people people (**people = $people**) and each person wants $avg slices on average (**avg = $avg**).  A pizza has $slices slices per pizza (**slices = $slices**).  How many pizzas should you order (**pizzas = ?**)?  So we have the following

Meaning | Variable
:------ | :--------:
Number of people | people
Average number of slices each person eats | avg
Number of slices on a piece of pizza | slices

"""

# ╔═╡ 444e2fa4-d460-11ea-12aa-57e0576c2d66
pizzas = 1

# ╔═╡ 3c12f2b4-d471-11ea-2d37-539f061f7cf2
r = 6

# ╔═╡ d9c31dfa-d470-11ea-23b2-838975b71f7c
md"""## Writing your own math functions

The area of a pizza is ``A = \pi r^2``.  Lets try calculating the area of a pizza that has a radius of $r inches (**r = $6**).  Type **pi** to get the value of ``\pi`` and **r^2** to get the radius squared.
"""

# ╔═╡ 50f0f6d6-d471-11ea-304e-8f72e7ef9d7e
A =  r^2

# ╔═╡ cb36a9ee-d472-11ea-1835-bf7963137e18
area(d) = pi * (d / 2)^2

# ╔═╡ 04b010c0-d473-11ea-1767-136c7e26e122
A2 = area(r)

# ╔═╡ 637c26fa-d475-11ea-2c5b-2b0f4775b119
smalls_in_xl = 1

# ╔═╡ 3823d09e-d474-11ea-194e-59b5805f303b
small = 13.10 / area(9)

# ╔═╡ 76c11174-d474-11ea-29c5-81856d47cf74
medium = 20.95 / area(13)

# ╔═╡ 8b12d200-d474-11ea-3035-01eccf39f917
large = 24.90 / area(15)

# ╔═╡ 962e6b86-d474-11ea-11a6-a1d11e33ae42
xl = 30.95 / area(17)

# ╔═╡ 16ec3f32-d4ff-11ea-20e2-5bc6dd5db083
best_value = small

# ╔═╡ f147b6cc-d4ff-11ea-05ad-6f5b441e5d1b
two_medium_cost = 20.95 * 1 - 0

# ╔═╡ 19eb2a82-d500-11ea-3782-596adc689382
two_medium_area = 1 * area(13)

# ╔═╡ 70e85498-d500-11ea-35af-474574f5c011
two_medium_deal = 1

# ╔═╡ 6494e270-d503-11ea-38a7-df96e7f0a241
cuts2 = 1

# ╔═╡ a05aae8e-d506-11ea-190f-57e9ce53b8b9
cuts3 = 1

# ╔═╡ 5a8ede88-d507-11ea-30d9-c99a67243781
cuts4 = 1

# ╔═╡ bae0cb62-dcc2-11ea-0667-512e1c407d40
cuts5 = 1

# ╔═╡ f5f89724-d507-11ea-0a93-6d904f36bbe4
function pieces(n)
	return n
end

# ╔═╡ bd9f3d24-d509-11ea-165d-3d465a0b4542
md"""Move the slider to change the number of cuts:

$(@bind n html"<input type=range max=50>")"""

# ╔═╡ b8644fb0-daa6-11ea-1e94-9bf46e7b0fad
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));

# ╔═╡ 8700d986-d475-11ea-0d0e-790448cf92ba
let
	ans = (pi * (17 / 2)^2) / (pi * (9 / 2)^2)
	if smalls_in_xl == 1
		hint(md"""The diameter of the XL pizza is 17 inches while the diameter of the small pizza is 9 inches.  Use the **area()** function from before to find the area of each and divide them.""")
	elseif smalls_in_xl < ans - 4 * eps(ans)
		md"""Keep trying, your answer is too low."""
	elseif smalls_in_xl > ans + 4 * eps(ans)
		md"""Keep trying, your answer is too high."""
	else
		md"""**Great!** You got it right. Let's move on."""
	end
end

# ╔═╡ 2814a1d4-dcc0-11ea-3d42-f52765e478fe
hint(md"For each extra cut, start out with the solution for the previous number. When you add one extra cut, how many new slices do you get?")

# ╔═╡ 48647ab2-daa5-11ea-0494-ef87be7cbf7c
hint(md"A new cut will create the maximum number of _new slices_ if it intersects all previous cuts.")

# ╔═╡ 8cada086-daa5-11ea-220c-0f660938b604
if cuts5 == 5 + 4 + 3 + 2 + 1 + 1
	hint(md"To get the maximum number of pieces with 5 cuts it will be ``5 + 4 + 3 + 2 + 1``, plus 1 extra for the original pizza with 0 cuts. To find the formula of a sequence of numbers group them like so: ``5 + (4 + 1) + (3 + 2) = 3 * 5``.")
else
	md""
end

# ╔═╡ 4119d19e-dcbc-11ea-3ec8-271e88e1afca
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]));

# ╔═╡ 921bba30-dcbc-11ea-13c3-87554722da8a
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]));

# ╔═╡ 5a6d1a8e-dcbc-11ea-272a-6f769c8d309c
correct(text=md"Great! You got the right answer! Let's move on to the next section.") = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]));

# ╔═╡ 33b1975c-d45c-11ea-035f-ab76e46a31ed
if c == a * b
	correct(md"""**Great!** The value of c = $c.  So you now have a simple computer!

	Now go back above and change the value of **a = $a** to **a = $(a + 3)** and press **`Shift-Enter`**.
	What is the new value of **c**?  Notice how all the values get updated in this notebook!
	""")
else
	keep_working()
end

# ╔═╡ f26d50da-d46b-11ea-0c2d-77ca13532b3d
if pizzas == people * avg / slices
	almost(md"Yes that is right! But we should round $pizzas up to an integer, otherwise the restaurant will be confused.

Try `ceil(...)`!")
elseif pizzas == ceil(people * avg / slices)
	correct(md"Yes that is right, that's a lot of pizza! Excellent, you figured out we need to round up the number of pizzas!")
else
	keep_working()
end

# ╔═╡ 5c4a5f22-d471-11ea-260f-9338d8bfa2d6
if A != pi * r^2
	keep_working(md"Let's fix the above cell before we move on!  Find the formula to calculate the area using **pi** and **r**.")
else
	correct(md"""**Great!**  You figured it out.  Keep going.""")
end

# ╔═╡ a07e5c3e-d476-11ea-308c-718f8f128334
if A2 != pi * (12 / 2)^2
	hint(md"Keep trying to get the right answer.
	**Hint**: you need to multiply the radius by 2 to convert it into the diameter.")
else
	correct()
end

# ╔═╡ 1ba2c208-d4ff-11ea-0a8e-e75bf7e1c3e6
if !isapprox(best_value, xl)
	hint(md"No need to copy these digits yourself - what should we assign to **best_value**?")
else
	correct()
end

# ╔═╡ 6ae748b2-d503-11ea-1c51-6b2df24fd212
if cuts2 != 4
	hint(md"The cuts must go all the way across the pizza!")
else
	correct(md"Awesome!")
end

# ╔═╡ a679bddc-d506-11ea-143a-6d4dcd70e918
if cuts3 == 6
	almost(md"""Close but not quite. Who said that pizza slices need to look like pizza slices?""")
elseif cuts3 == 7
	correct(md"You got it right.  Now for something harder...")
else
	hint(md"Try drawing it out on a piece of paper.")
end

# ╔═╡ 5df7eefc-d507-11ea-0d1f-45b224a04774
if cuts4 == 11
	correct(md"That was a tough question.  How did you figure it out?  You tried hard.")
elseif cuts4 < 10
	hint(md"Draw it out on a piece of paper.  You can make more pieces with 4 cuts.")
elseif cuts4 < 11
	hint(md"Getting close but you can make more pieces with 4 cuts.")
else
	hint(md"That is too high. Only straight lines!")
end

# ╔═╡ e80986c6-d509-11ea-12e3-f79a54b5ab31
if pieces(n) ==  n * (n + 1) / 2 + 1
	md"""_Testing..._

	**For $n cuts, you predict $(pieces(n)) pieces.**

	$(correct(md"Well done!"))"""
else
	md"""_Testing..._

	**For $n cuts, you predict $(pieces(n)) pieces.**

	$(keep_working(md"The answer should be $(Int(n*(n+1)/2+1))."))"""
end

# ╔═╡ Cell order:
# ╟─03664f5c-d45c-11ea-21b6-91cd647a07aa
# ╠═14158eb0-d45c-11ea-088f-330e45412320
# ╠═2ed4bb92-d45c-11ea-0b31-2d8e32ce7b44
# ╟─e5e0a0da-d45c-11ea-1042-e9b5d0654d4f
# ╠═30f0f882-d45c-11ea-2adc-7d84ecf8a7a6
# ╟─33b1975c-d45c-11ea-035f-ab76e46a31ed
# ╟─a38cb92e-d45e-11ea-2959-05be909befb2
# ╠═262b312a-d460-11ea-26c5-df30459effc5
# ╠═2ea7f162-d460-11ea-0e8e-25340e2e64da
# ╠═3da812c6-d460-11ea-0170-79fbb6a4347c
# ╟─4dff4b5e-d461-11ea-29c8-d548fdb5f08b
# ╠═444e2fa4-d460-11ea-12aa-57e0576c2d66
# ╟─f26d50da-d46b-11ea-0c2d-77ca13532b3d
# ╟─d9c31dfa-d470-11ea-23b2-838975b71f7c
# ╠═3c12f2b4-d471-11ea-2d37-539f061f7cf2
# ╠═50f0f6d6-d471-11ea-304e-8f72e7ef9d7e
# ╟─5c4a5f22-d471-11ea-260f-9338d8bfa2d6
# ╟─f907e46a-d471-11ea-07e5-f30e2aab3d08
# ╠═cb36a9ee-d472-11ea-1835-bf7963137e18
# ╟─d9575e9c-d472-11ea-1eda-2d335d039f28
# ╠═04b010c0-d473-11ea-1767-136c7e26e122
# ╟─a07e5c3e-d476-11ea-308c-718f8f128334
# ╟─edb95b14-d473-11ea-3a5a-77382d31f941
# ╠═637c26fa-d475-11ea-2c5b-2b0f4775b119
# ╟─8700d986-d475-11ea-0d0e-790448cf92ba
# ╟─5b07b8fe-d475-11ea-01aa-6b88d6ed8a05
# ╠═3823d09e-d474-11ea-194e-59b5805f303b
# ╠═76c11174-d474-11ea-29c5-81856d47cf74
# ╠═8b12d200-d474-11ea-3035-01eccf39f917
# ╠═962e6b86-d474-11ea-11a6-a1d11e33ae42
# ╟─a42e4eb0-d474-11ea-316a-3d864451bc01
# ╠═16ec3f32-d4ff-11ea-20e2-5bc6dd5db083
# ╟─1ba2c208-d4ff-11ea-0a8e-e75bf7e1c3e6
# ╟─cb419286-d4ff-11ea-1d7f-af5c8574b775
# ╠═f147b6cc-d4ff-11ea-05ad-6f5b441e5d1b
# ╟─0d76d97c-d500-11ea-2433-e96c6fc43b05
# ╠═19eb2a82-d500-11ea-3782-596adc689382
# ╟─20a1e9cc-d500-11ea-3d9b-279c71bc20f1
# ╠═70e85498-d500-11ea-35af-474574f5c011
# ╟─57f024ae-d500-11ea-1cc4-ed28348fdf93
# ╟─180c8fdc-d503-11ea-04ca-bf2c07fd1c17
# ╠═6494e270-d503-11ea-38a7-df96e7f0a241
# ╟─6ae748b2-d503-11ea-1c51-6b2df24fd212
# ╟─92b4a012-d503-11ea-15a2-1f3a446d3284
# ╠═a05aae8e-d506-11ea-190f-57e9ce53b8b9
# ╟─a679bddc-d506-11ea-143a-6d4dcd70e918
# ╟─2eb9a560-d507-11ea-3b8b-9d06678fe131
# ╠═5a8ede88-d507-11ea-30d9-c99a67243781
# ╟─5df7eefc-d507-11ea-0d1f-45b224a04774
# ╟─d1e3dec0-d507-11ea-1213-d37a9325ee2f
# ╟─2814a1d4-dcc0-11ea-3d42-f52765e478fe
# ╟─48647ab2-daa5-11ea-0494-ef87be7cbf7c
# ╟─97bfd13c-dcc2-11ea-0067-ad8c2c6517fc
# ╠═bae0cb62-dcc2-11ea-0667-512e1c407d40
# ╟─8cada086-daa5-11ea-220c-0f660938b604
# ╟─e0cb2822-dcc2-11ea-2c85-5748bfe526dc
# ╠═f5f89724-d507-11ea-0a93-6d904f36bbe4
# ╟─03249876-d508-11ea-16bb-fd5afed37a1f
# ╟─bd9f3d24-d509-11ea-165d-3d465a0b4542
# ╟─e80986c6-d509-11ea-12e3-f79a54b5ab31
# ╟─b8644fb0-daa6-11ea-1e94-9bf46e7b0fad
# ╟─4119d19e-dcbc-11ea-3ec8-271e88e1afca
# ╟─921bba30-dcbc-11ea-13c3-87554722da8a
# ╟─5a6d1a8e-dcbc-11ea-272a-6f769c8d309c

max_stars_repo_path: sample/Basic mathematics.jl

max_issues_repo_path: sample/Basic mathematics.jl

max_forks_repo_path: sample/Basic mathematics.jl


text:
module FEAutodiffTests

using LinearAlgebra
using Gridap.Algebra
using Gridap.FESpaces
using Gridap.Arrays
using Gridap.Fields
using Gridap.Geometry
using Gridap.TensorValues
using Gridap.CellData
using Gridap.ReferenceFEs

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

Ω = Triangulation(model)
dΩ = Measure(Ω,2)

V = FESpace(model,ReferenceFE(lagrangian,Float64,2),conformity=:H1)
U = V

dv = get_cell_shapefuns(V)
du = get_cell_shapefuns_trial(U)
uh = FEFunction(U,rand(num_free_dofs(U)))

ener(uh) = ∫( 0.5*∇(uh)⋅∇(uh) )*dΩ
res(uh) = ∫(∇(uh)⋅∇(dv))*dΩ
jac(uh) = ∫(∇(du)⋅∇(dv))*dΩ

cell_r = get_array(res(uh))
cell_j = get_array(jac(uh))
cell_h = cell_j

cell_r_auto = get_array(gradient(ener,uh))
cell_j_auto = get_array(jacobian(res,uh))
cell_h_auto = get_array(hessian(ener,uh))

test_array(cell_r_auto,cell_r,≈)
test_array(cell_j_auto,cell_j,≈)
test_array(cell_h_auto,cell_h,≈)

Γ = BoundaryTriangulation(model)
dΓ = Measure(Γ,2)

ener(uh) = ∫( 0.5*∇(uh)⋅∇(uh) )*dΓ
res(uh) = ∫( ∇(uh)⋅∇(dv) )*dΓ
jac(uh) = ∫( ∇(du)⋅∇(dv) )*dΓ

cell_r = get_array(res(uh))
cell_j = get_array(jac(uh))
cell_h = cell_j

cell_r_auto = get_array(gradient(ener,uh))
cell_j_auto = get_array(jacobian(res,uh))
cell_h_auto = get_array(hessian(ener,uh))

test_array(cell_r_auto,cell_r,≈)
test_array(cell_j_auto,cell_j,≈)
test_array(cell_h_auto,cell_h,≈)

ener(uh) = ∫( 0.5*∇(uh)⋅∇(uh) )*dΓ + ∫( 0.5*∇(uh)⋅∇(uh) )*dΩ
res(uh) = ∫( ∇(uh)⋅∇(dv) )*dΓ + ∫(∇(uh)⋅∇(dv))*dΩ
jac(uh) = ∫( ∇(du)⋅∇(dv) )*dΓ + ∫(∇(du)⋅∇(dv))*dΩ

cell_r = res(uh)
cell_j = jac(uh)
cell_h = cell_j

cell_r_auto = gradient(ener,uh)
cell_j_auto = jacobian(res,uh)
cell_h_auto = hessian(ener,uh)

test_array(cell_r_auto[Ω],cell_r[Ω],≈)
test_array(cell_j_auto[Ω],cell_j[Ω],≈)
test_array(cell_h_auto[Ω],cell_h[Ω],≈)

test_array(cell_r_auto[Γ],cell_r[Γ],≈)
test_array(cell_j_auto[Γ],cell_j[Γ],≈)
test_array(cell_h_auto[Γ],cell_h[Γ],≈)

const p = 3
j(∇u) = norm(∇u)^(p-2) * ∇u
dj(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2)*∇du
f(x) = 0

res(u,v) = ∫( ∇(v)⋅(j∘∇(u)) - v*f)*dΩ
jac(u,du,v) = ∫( ∇(v)⋅(dj∘(∇(du),∇(u))) )*dΩ

cell_j = get_array(jac(uh,du,dv))
cell_j_auto = get_array(jacobian(u->res(u,dv),uh))

test_array(cell_j_auto,cell_j,≈)

end # module


text:
using .Curve: BN, EP, EP2, EPX, Limb, md_sha256
using .Config: PRIVATE_KEY_SIZE_BITS, ORDER, G1, G2, CAN_STUFF_Y

const Point = EPX

function decode(::Type{T}, bytes::Vector{UInt8}) where {T <: Point}
    tmp = UInt8[
        0x2 | (bytes[1] & 0x80) >> 7,
        bytes[1] & 0x7f,
        bytes[2:end]...
    ]
    return T(tmp)
end

function encode(point::Point)
    bufsize = sizeof(point.x) + 1
    buf = Vector{UInt8}(point)
    @assert CAN_STUFF_Y # other case not implemented
    buf[2] |= (buf[1] & 0x1) << 7
    deleteat!(buf, 1)
    return buf
end

function gensk()
    # Throw away keys that are outside of the allowed key space
    # (mod break the uniform key distribution)
    while true
        sk = rand(BN, bits=PRIVATE_KEY_SIZE_BITS)
        if sk < ORDER
            return sk
        end
    end
end

#hash2curve(msg::Vector{UInt8}) = Curve.curve_map(Config.@EP, msg)
hash2curve(msg::Vector{UInt8}) = BN(md_sha256(msg)) * G1

sign(sk::BN, p::Point) = sk * p
genpk(sk::BN) = sk * G2


text:
module CP
include("dynam_anim.jl")

using Distributions,Luxor

struct Params
  λ::Float64
end

getdefault(xs,i,def) = checkbounds(Bool,xs,i) ? @inbounds(xs[i]) : def

rate(σ::Vector{Bool},p::Params,i::Int)::Float64 = σ[i] ? 1 : p.λ * sum(getdefault(σ,i+j,0) for j in [-1,1])

rates(σ::Vector{Bool},p::Params) = [rate(σ,p,i) for i in eachindex(σ)]

function step(σ::Vector{Bool},p::Params)::Float64
  rs = rates(σ,p)
  r = sum(rs)
  if r == 0
    return Inf
  end
  p = rs/r
  i = rand(Categorical(p))
  σ[i] = !σ[i]
  rand(Exponential(1/r))
end

function setup(σ::Vector{Bool},w,h)
  Tiler(w,h,1,length(σ),margin=0)
end

function draw(σ::Vector{Bool},tiles)
  for (pos,n) in tiles
    setcolor(σ[n] ? "red" : "green")
    box(pos, tiles.tilewidth, tiles.tileheight, :fill)
  end
end

function contact_anim(init,p,frames,steps_per_frame,width,height,path)
  ds = DynamAnim.DynamicalSystem(init,p,step)
  ap = DynamAnim.AnimParams(setup,draw,frames,steps_per_frame,width,height,"Contact Process")
  DynamAnim.anim(ds,ap,path)
end

end


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking
using CSV, DataFrames
using StanSample, MonteCarloMeasurements

ProjDir = @__DIR__

df = DataFrame!(CSV.File(joinpath(ProjDir, "sports.csv"), delim=','));
scale!(df, [:TM, :WUE, :C])
println()

m1 = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] WUE; // Outcome
 vector[N] TM; // Predictor
 vector[N] C; // Predictor
}

parameters {
 real a; // Intercept
 real bC; // Slope (regression coefficients)
 real bTM; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);         //Priors
  bC ~ normal(0, 0.5);
  bTM ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bC * C + bTM * TM;
  WUE ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m1s = SampleModel("m1", m1, [6];
  method=StanSample.Sample(num_samples=1000));

# Input data for cmdstan

m1_data = Dict("N" => size(df, 1), "C" => df[:, :C_s],
    "TM" => df[:, :TM_s], "WUE" => df[:, :WUE_s]);

# Sample using StanSample

rc = stan_sample(m1s, data=m1_data);

if success(rc)

  # Describe the draws

  dfa1 = read_samples(m1s; output_format=:dataframe)

  p1 = Particles(dfa1)
  display(p1)

end


text:
struct SARIMAOrder
    p::Int
    d::Int
    q::Int
    P::Int
    D::Int
    Q::Int
    s::Int
    r::Int # Model order
    n_states_diff::Int
    n_states::Int
    function SARIMAOrder(p::Int, d::Int, q::Int, P::Int, D::Int, Q::Int, s::Int)
        @assert p >= 0
        @assert d >= 0
        @assert q >= 0
        @assert P >= 0
        @assert D >= 0
        @assert Q >= 0
        @assert s >= 0
        r = max(p + s * P, q + s * Q + 1)
        n_states_diff = d + D * s
        n_states = d + D * s + r
        return new(p, d, q, P, D, Q, s, r, n_states_diff, n_states)
    end
end

# This is used to some optimizations such as
# defining the polynomials beforehand
# and avoiding string interpolations such as "ar_L$i"
# for "ar_L1"
mutable struct SARIMAHyperParametersAuxiliary{Fl <: AbstractFloat}
    ar_poly::Vector{Fl}
    ar_names::Vector{String}
    seasonal_ar_poly::Vector{Fl}
    seasonal_ar_names::Vector{String}
    ma_poly::Vector{Fl}
    ma_names::Vector{String}
    seasonal_ma_poly::Vector{Fl}
    seasonal_ma_names::Vector{String}
    function SARIMAHyperParametersAuxiliary{Fl}(order::SARIMAOrder) where Fl
        ar_poly = create_poly(order.p, Fl)
        seasonal_ar_poly = create_seasonal_poly(order.P, order.s, Fl)
        ma_poly = create_poly(order.q, Fl)
        seasonal_ma_poly = create_seasonal_poly(order.Q, order.s, Fl)
        return new{Fl}(
            ar_poly,
            create_names("ar_L", order.p),
            seasonal_ar_poly,
            create_names("s_ar_L", order.P, order.s),
            ma_poly,
            create_names("ma_L", order.q),
            seasonal_ma_poly,
            create_names("s_ma_L", order.Q, order.s)
        )
    end
end

@doc raw"""
    SARIMA(
        y::Vector{Fl};
        order::Tuple{Int,Int,Int} = (1, 0, 0),
        seasonal_order::Tuple{Int, Int, Int, Int} = (0, 0, 0, 0),
        include_mean::Bool = false
    ) where Fl

A SARIMA model (Seasonal AutoRegressive Integrated Moving Average) implemented within the state-space
framework.

The SARIMA model is specified ``(p, d, q) \times (P, D, Q, s)``. We can also consider a polynomial ``A(t)`` to
model a trend, here we only allow to add a constante term with the `include_mean` keyword argument.
```math
\begin{gather*}
    \begin{aligned}
        \phi_p (L) \tilde \phi_P (L^s) \Delta^d \Delta_s^D y_t = A(t) + \theta_q (L) \tilde \theta_Q (L^s) \zeta_t
    \end{aligned}
\end{gather*}
```
In terms of a univariate structural model, this can be represented as
```math
\begin{gather*}
    \begin{aligned}
    y_t & = u_t + \eta_t \\
    \phi_p (L) \tilde \phi_P (L^s) \Delta^d \Delta_s^D u_t & = A(t) + \theta_q (L) \tilde \theta_Q (L^s) \zeta_t
    \end{aligned}
\end{gather*}
```

# Example
```jldoctest
julia> model = SARIMA(rand(100); order=(1,1,1), seasonal_order=(1,2,3,12))
SARIMA(1, 1, 1)x(1, 2, 3, 12) model
```

See more on [Airline passengers](@ref)

# References
 * Durbin, James, & Siem Jan Koopman.
   Time Series Analysis by State Space Methods: Second Edition.
   Oxford University Press, 2012
"""
mutable struct SARIMA <: StateSpaceModel
    order::SARIMAOrder
    hyperparameters_auxiliary::SARIMAHyperParametersAuxiliary
    hyperparameters::HyperParameters
    system::LinearUnivariateTimeInvariant
    results::Results
    include_mean::Bool

    function SARIMA(y::Vector{Fl};
                    order::Tuple{Int,Int,Int} = (1, 0, 0),
                    seasonal_order::Tuple{Int, Int, Int, Int} = (0, 0, 0, 0),
                    include_mean::Bool = false) where Fl
        or = SARIMAOrder(order[1], order[2], order[3],
                        seasonal_order[1], seasonal_order[2], seasonal_order[3], seasonal_order[4])
        hyperparameters_auxiliary = SARIMAHyperParametersAuxiliary{Fl}(or)

        Z = build_Z(or, Fl)
        T = build_T(or, Fl)
        R = build_R(or, Fl)
        d_ss = zero(Fl)
        c = build_c(or, Fl)
        H = zero(Fl)
        Q_ss = zeros(Fl, 1, 1)

        system = LinearUnivariateTimeInvariant{Fl}(y, Z, T, R, d_ss, c, H, Q_ss)

        num_hyperparameters = or.p + or.q + or.P + or.Q + 1
        if include_mean
            num_hyperparameters += 1
        end
        names = Vector{String}(undef, num_hyperparameters)
        for i in 1:(or.p)
            names[i] = "ar_L$i"
        end
        for j in 1:(or.q)
            names[or.p + j] = "ma_L$j"
        end
        for is in 1:(or.P)
            names[or.p + or.q + is] = "s_ar_L$(or.s * is)"
        end
        for js in 1:(or.Q)
            names[or.p + or.q + or.P + js] = "s_ma_L$(or.s * js)"
        end
        if include_mean
            names[end-1] = "mean"
        end
        names[end] = "sigma2_η"

        hyperparameters = HyperParameters{Fl}(names)

        return new(or, hyperparameters_auxiliary, hyperparameters, system, Results{Fl}(), include_mean)
    end
end

# Auxiliary functions
function build_Z(or, Fl)
    seasonal_entries = or.s - 1 >= 0 ? or.s - 1 : 0
    return vcat(
            fill(one(Fl), or.d),
            repeat(vcat(zeros(Fl, seasonal_entries), one(Fl)), or.D),
            one(Fl),
            fill(zero(Fl), or.r - 1)
        )
end

triu_indices(k::Int) = findall(isone, triu(ones(k, k)))
function build_T(or, Fl)
    # zeros matrix in the correct dimension
    T = zeros(Fl, or.n_states, or.n_states)
    if or.r > 0
        T[end - or.r + 1:end, end - or.r + 1:end] = diagm(1 => ones(or.r - 1))
    end
    if or.D > 0
        seasonal_companion = diagm(-1 => ones(or.s - 1))
        seasonal_companion[1, end] = 1
        for l in 1:or.D
            start_idx = or.d + (l - 1) * or.s
            end_idx = or.d + l * or.s
            T[start_idx + 1:end_idx, start_idx + 1:end_idx] .= seasonal_companion
            if l + 1 <= or.D
                T[start_idx + 1, end_idx + or.s] = 1
            end
            T[start_idx + 1, or.n_states_diff + 1] = 1
        end
    end
    if or.d > 0
        idx = triu_indices(or.d)
        T[idx] .= 1
        if or.s > 0
            start_idx = or.d
            end_idx = or.n_states_diff
            T[1:or.d, start_idx+1:end_idx] .= repeat(hcat(zeros(Fl, 1, or.s - 1), one(Fl)), 1, or.D)
        end
        T[1:or.d, or.n_states_diff + 1] .= 1
    end
    return T
end
function build_R(or, Fl)
    return vcat(
        fill(zero(Fl), or.n_states_diff, 1),
        fill(one(Fl), 1, 1),
        fill(zero(Fl), or.r - 1, 1)
    )
end
function build_c(or, Fl)
    return zeros(Fl, or.n_states)
end

function create_names(str::String, dim::Int, s::Int = 1)
    names = Vector{String}(undef, dim)
    for i in 1:dim
        names[i] = str*"$(s * i)"
    end
    return names
end
function create_poly(n::Int, Fl::DataType)
    return fill(one(Fl), n + 1)
end
function create_seasonal_poly(n::Int, s::Int, Fl::DataType)
    seasonal_poly = fill(zero(Fl), (s * n) + 1)
    seasonal_poly[1] = one(Fl)
    for i in 1:n
        seasonal_poly[(i * s) + 1] = one(Fl)
    end
    return seasonal_poly
end
function poly_mul(p1::Vector{Fl}, p2::Vector{Fl}) where Fl
    n, m = length(p1), length(p2)
    if n > 1 && m > 1
        c = zeros(Fl, m + n - 1)
        for i in 1:n, j in 1:m
            c[i + j - 1] += p1[i] * p2[j]
        end
        return c
    elseif n <= 1
        return p2
    else
        return p1
    end
end
get_ar_name(model::SARIMA, i::Int) = model.hyperparameters_auxiliary.ar_names[i]
get_ma_name(model::SARIMA, i::Int) = model.hyperparameters_auxiliary.ma_names[i]
get_seasonal_ar_name(model::SARIMA, i::Int) = model.hyperparameters_auxiliary.seasonal_ar_names[i]
get_seasonal_ma_name(model::SARIMA, i::Int) = model.hyperparameters_auxiliary.seasonal_ma_names[i]

function update_T_terms!(model::SARIMA)
    update_ar_poly!(model)
    update_seasonal_ar_poly!(model)
    # TODO increaase performance
    # We could preallocate this polynomial p3
    p3 = poly_mul(model.hyperparameters_auxiliary.ar_poly, model.hyperparameters_auxiliary.seasonal_ar_poly)
    for i in model.order.n_states_diff + 1:model.order.n_states_diff + model.order.p + model.order.s * model.order.P
        model.system.T[i, model.order.n_states_diff + 1] = p3[i - model.order.n_states_diff + 1]
    end
    return nothing
end

function update_R_terms!(model::SARIMA)
    update_ma_poly!(model)
    update_seasonal_ma_poly!(model)
    # TODO increaase performance
    # We could preallocate this polynomial p3
    p3 = poly_mul(model.hyperparameters_auxiliary.ma_poly, model.hyperparameters_auxiliary.seasonal_ma_poly)
    for i in model.order.n_states_diff + 1:model.order.n_states_diff + model.order.q + model.order.s * model.order.Q
        model.system.R[i + 1] = p3[i - model.order.n_states_diff + 1]
    end
    return nothing
end

function update_ar_poly!(model::SARIMA)
    if model.order.p > 0
        for i in 1:(model.order.p)
            model.hyperparameters_auxiliary.ar_poly[i + 1] = get_constrained_value(
                model, get_ar_name(model, i)
            )
        end
    end
end
function update_seasonal_ar_poly!(model::SARIMA)
    if model.order.P > 0
        for i in 1:(model.order.P)
            model.hyperparameters_auxiliary.seasonal_ar_poly[(i * model.order.s) + 1] = get_constrained_value(
                model, get_seasonal_ar_name(model, i)
            )
        end
    end
end
function update_ma_poly!(model::SARIMA)
    if model.order.q > 0
        for i in 1:(model.order.q)
            model.hyperparameters_auxiliary.ma_poly[i + 1] = get_constrained_value(
                model, get_ma_name(model, i)
            )
        end
    end
end
function update_seasonal_ma_poly!(model::SARIMA)
    if model.order.Q > 0
        for i in 1:(model.order.Q)
            model.hyperparameters_auxiliary.seasonal_ma_poly[(i * model.order.s) + 1] = get_constrained_value(
                model, get_seasonal_ma_name(model, i)
            )
        end
    end
end

function impose_unit_root_constraint(hyperparameter_values::Vector{Fl}) where Fl
    # TODO improve performance possibly preallocating y
    #     Monahan, John F. 1984.
    #    "A Note on Enforcing Stationarity in
    #    Autoregressive-moving Average Models."
    #    Biometrika 71 (2) (August 1): 403-404.
    n = length(hyperparameter_values)
    y = fill(zero(Fl), n, n)
    r = @. hyperparameter_values / sqrt((1 + hyperparameter_values^2))
    if n == 1
        return -r
    end
    @inbounds for k in 1:n
        for i in 1:(k - 1)
            y[k, i] = y[k - 1, i] + r[k] * y[k - 1, k - i]
        end
        y[k, k] = r[k]
    end
    return -@view y[n, :]
end

function relax_unit_root_constraint(hyperparameter_values::Vector{Fl}) where Fl
    #     Monahan, John F. 1984.
    #    "A Note on Enforcing Stationarity in
    #    Autoregressive-moving Average Models."
    #    Biometrika 71 (2) (August 1): 403-404.
    n = length(hyperparameter_values)
    y = fill(zero(Fl), n, n)
    y[n, :] = -hyperparameter_values
    for k in n:-1:2, i in 1:(k - 1)
        y[k - 1, i] = (y[k, i] - y[k, k] * y[k, k - i]) / (1 - y[k, k]^2)
    end
    r = diag(y)
    x = r ./ sqrt.(1 .- r .^ 2)
    return x
end

function SARIMA_exact_initialization!(filter::KalmanFilter, model::SARIMA)
    SARIMA_exact_initialization!(filter.kalman_state, model.system, model.order.r)
    return nothing
end

function SARIMA_exact_initialization!(kalman_state,
                                      system,
                                      num_arma_states::Int)
    # TODO we could preallocate these arrays
    arma_T = system.T[(end - num_arma_states + 1):end, (end - num_arma_states + 1):end]
    arma_R = system.R[(end - num_arma_states + 1):end]
    # Calculate exact initial P1
    arma_P1 = lyapd(arma_T, arma_R * arma_R')
    # Fill P1 in the arma states
    kalman_state.P[(end - num_arma_states + 1):end, (end - num_arma_states + 1):end] .=
        system.Q[1] .* arma_P1
    return nothing
end

function concatenate_on_bottom(X1::Matrix{Fl}, X2::Matrix{Fl}) where Fl
    n = min(size(X1, 1), size(X2, 1))
    return hcat(X1[end-n+1:end, :], X2[end-n+1:end, :])
end

function diff_sarima(y::Vector{Fl}, d::Int, D::Int, s::Int) where Fl
    # Seasonal differencing
    for _ in 1:D
        y = y[s:end] - y[1:end-s+1]
    end
    # Simple differencing
    for _ in 1:d
        y = diff(y)
    end
    return y
end

function constrain_ar!(model::SARIMA, Fl::DataType)
    if model.order.p > 0
        unconstrained_ar = Vector{Fl}(undef, model.order.p)
        for i in 1:(model.order.p)
            unconstrained_ar[i] = get_unconstrained_value(model, get_ar_name(model, i))
        end
        constrained_ar = impose_unit_root_constraint(unconstrained_ar)
        for i in 1:(model.order.p)
            update_constrained_value!(model, get_ar_name(model, i), constrained_ar[i])
        end
    end
    return nothing
end
function constrain_ma!(model::SARIMA, Fl::DataType)
    if model.order.q > 0
        unconstrained_ma = Vector{Fl}(undef, model.order.q)
        for j in 1:(model.order.q)
            unconstrained_ma[j] = get_unconstrained_value(model, get_ma_name(model, j))
        end
        constrained_ma = -impose_unit_root_constraint(unconstrained_ma)
        for j in 1:(model.order.q)
            update_constrained_value!(model, get_ma_name(model, j), constrained_ma[j])
        end
    end
    return nothing
end
function constrain_seasonal_ar!(model::SARIMA, Fl::DataType)
    if model.order.P > 0
        unconstrained_seasonal_ar = Vector{Fl}(undef, model.order.P)
        for i in 1:(model.order.P)
            unconstrained_seasonal_ar[i] = get_unconstrained_value(model, get_seasonal_ar_name(model, i))
        end
        constrained_seasonal_ar = impose_unit_root_constraint(unconstrained_seasonal_ar)
        for i in 1:(model.order.P)
            update_constrained_value!(model, get_seasonal_ar_name(model, i), constrained_seasonal_ar[i])
        end
    end
    return nothing
end
function constrain_seasonal_ma!(model::SARIMA, Fl::DataType)
    if model.order.Q > 0
        unconstrained_seasonal_ma = Vector{Fl}(undef, model.order.Q)
        for j in 1:(model.order.Q)
            unconstrained_seasonal_ma[j] = get_unconstrained_value(model, get_seasonal_ma_name(model, j))
        end
        constrained_seasonal_ma = -impose_unit_root_constraint(unconstrained_seasonal_ma)
        for j in 1:(model.order.Q)
            update_constrained_value!(model, get_seasonal_ma_name(model, j), constrained_seasonal_ma[j])
        end
    end
    return nothing
end
function constrain_mean!(model::SARIMA)
    if model.include_mean
        constrain_identity!(model, "mean")
    end
    return nothing
end

function unconstrain_ar!(model::SARIMA, Fl::DataType)
    if model.order.p > 0
        constrained_ar = Vector{Fl}(undef, model.order.p)
        for i in 1:(model.order.p)
            constrained_ar[i] = get_constrained_value(model, get_ar_name(model, i))
        end
        unconstrained_ar = relax_unit_root_constraint(constrained_ar)
        for i in 1:(model.order.p)
            update_unconstrained_value!(model, get_ar_name(model, i), unconstrained_ar[i])
        end
    end
    return nothing
end
function unconstrain_ma!(model::SARIMA, Fl::DataType)
    if model.order.q > 0
        constrained_ma = Vector{Fl}(undef, model.order.q)
        for j in 1:(model.order.q)
            constrained_ma[j] = get_constrained_value(model, get_ma_name(model, j))
        end
        unconstrained_ma = -relax_unit_root_constraint(constrained_ma)
        for j in 1:(model.order.q)
            update_unconstrained_value!(model, get_ma_name(model, j), unconstrained_ma[j])
        end
    end
    return nothing
end
function unconstrain_seasonal_ar!(model::SARIMA, Fl::DataType)
    if model.order.P > 0
        constrained_seasonal_ar = Vector{Fl}(undef, model.order.P)
        for i in 1:(model.order.P)
            constrained_seasonal_ar[i] = get_constrained_value(model, get_seasonal_ar_name(model, i))
        end
        unconstrained_seasonal_ar = relax_unit_root_constraint(constrained_seasonal_ar)
        for i in 1:(model.order.P)
            update_unconstrained_value!(model, get_seasonal_ar_name(model, i), unconstrained_seasonal_ar[i])
        end
    end
    return nothing
end
function unconstrain_seasonal_ma!(model::SARIMA, Fl::DataType)
    if model.order.Q > 0
        constrained_seasonal_ma = Vector{Fl}(undef, model.order.Q)
        for j in 1:(model.order.Q)
            constrained_seasonal_ma[j] = get_constrained_value(model, get_seasonal_ma_name(model, j))
        end
        unconstrained_seasonal_ma = -relax_unit_root_constraint(constrained_seasonal_ma)
        for j in 1:(model.order.Q)
            update_unconstrained_value!(model, get_seasonal_ma_name(model, j), unconstrained_seasonal_ma[j])
        end
    end
    return nothing
end
function unconstrain_mean!(model::SARIMA)
    if model.include_mean
        unconstrain_identity!(model, "mean")
    end
    return nothing
end

function ar_polinomial(p::Vector{Fl}) where Fl
    return Polynomial([one(Fl); -p])
end

function ma_polinomial(q::Vector{Fl}) where Fl
    return Polynomial([one(Fl); q])
end

function roots_of_inverse_polinomial(poly::Polynomial)
    return roots(poly).^-1
end

function assert_stationarity(p::Vector{Fl}) where Fl
    poly = ar_polinomial(p)
    return all(abs.(roots_of_inverse_polinomial(poly)) .< 1)
end

function assert_invertibility(q::Vector{Fl}) where Fl
    poly = ma_polinomial(q)
    return all(abs.(roots_of_inverse_polinomial(poly)) .< 1)
end

function conditional_sum_of_squares(y_diff::Vector{Fl}, k_ar::Int, k_ma::Int) where Fl
    if (k_ar == 0) && (k_ma == 0)
        return (Fl[], Fl[])
    end
    k = 2 * k_ma
    r = max(k + k_ma, k_ar)
    residuals = nothing
    X = nothing
    if k_ar + k_ma > 0
        # If we have MA terms, get residuals from an AR(k) model to use
        # as data for conditional sum of squares estimates of the MA
        # parameters
        if k_ma > 0
            y = y_diff[k:end]
            X = lagmat(y_diff, k)
            params_ar = X \ y[end - size(X, 1) + 1:end]
            residuals = y[end - size(X, 1) + 1:end] - X * params_ar
        end
        # Run an ARMA(p,q) model using the just computed residuals as
        # data
        y = y_diff[r:end]
        X = Matrix{Fl}(undef, length(y), 0)
        X = concatenate_on_bottom(X, lagmat(y_diff, k_ar))
        if k_ma > 0
            X = concatenate_on_bottom(X, lagmat(residuals, k_ma))
        end
    end

    initial_params = X \ y_diff[end - size(X, 1) + 1:end]
    params_ar = Fl[]
    params_ma = Fl[]
    offset = 1
    if k_ar > 0
        params_ar = initial_params[offset:offset + k_ar - 1]
        offset += k_ar
    end
    if k_ma > 0
        params_ma = initial_params[offset:offset + k_ma - 1]
    end
    return params_ar, params_ma
end

# Obligatory functions
function default_filter(model::SARIMA)
    Fl = typeof_model_elements(model)
    a1 = zeros(Fl, model.order.n_states)
    P1 = Fl(1e6) .* Matrix{Fl}(I, model.order.n_states, model.order.n_states)
    skip_llk_instants = model.order.d + model.order.s * model.order.D
    steadystate_tol = Fl(1e-5)
    return UnivariateKalmanFilter(a1, P1, skip_llk_instants, steadystate_tol)
end

function initial_hyperparameters!(model::SARIMA)
    Fl = typeof_model_elements(model)
    initial_hyperparameters = Dict{String,Fl}()
    # Heuristic inspired in statsmodels from python
    # TODO find a reference to this heuristic
    # conditional sum of squares
    y_diff = diff_sarima(model.system.y, model.order.d, model.order.D, model.order.s)
    y_diff = filter(!isnan, y_diff)
    # Non-seasonal ARMA
    (initial_ar, initial_ma) = conditional_sum_of_squares(y_diff, model.order.p, model.order.q)
    if !assert_stationarity(initial_ar)
        @warn("Conditional sum of squares estimated initial_ar out of the unit circle, using zero as starting params")
        initial_ar .= zero(Fl)
    end
    if !assert_invertibility(initial_ma)
        @warn("Conditional sum of squares estimated initial_ma out of the unit circle, using zero as starting params")
        initial_ma .= zero(Fl)
    end
    # Seasonal ARMA
    (initial_seasonal_ar,
    initial_seasonal_ma) = conditional_sum_of_squares(y_diff, model.order.P, model.order.Q)
    if !assert_stationarity(initial_seasonal_ar)
        @warn("Conditional sum of squares estimated initial_seasonal_ar out of the unit circle, using zero as starting params")
        initial_seasonal_ar .= zero(Fl)
    end
    if !assert_invertibility(initial_seasonal_ma)
        @warn("Conditional sum of squares estimated initial_seasonal_ma out of the unit circle, using zero as starting params")
        initial_seasonal_ma .= zero(Fl)
    end
    for i in 1:(model.order.p)
        initial_hyperparameters[get_ar_name(model, i)] = initial_ar[i]
    end
    for j in 1:(model.order.q)
        initial_hyperparameters[get_ma_name(model, j)] = initial_ma[j]
    end
    for i in 1:(model.order.P)
        initial_hyperparameters[get_seasonal_ar_name(model, i)] = initial_seasonal_ar[i]
    end
    for j in 1:(model.order.Q)
        initial_hyperparameters[get_seasonal_ma_name(model, j)] = initial_seasonal_ma[j]
    end
    if model.include_mean
        initial_hyperparameters["mean"] = mean(y_diff)
    end

    initial_hyperparameters["sigma2_η"] = var(y_diff)
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return nothing
end

function constrain_hyperparameters!(model::SARIMA)
    Fl = typeof_model_elements(model)
    constrain_ar!(model, Fl)
    constrain_ma!(model, Fl)
    constrain_seasonal_ar!(model, Fl)
    constrain_seasonal_ma!(model, Fl)
    constrain_mean!(model)
    constrain_variance!(model, "sigma2_η")
    return nothing
end

function unconstrain_hyperparameters!(model::SARIMA)
    Fl = typeof_model_elements(model)
    unconstrain_ar!(model, Fl)
    unconstrain_ma!(model, Fl)
    unconstrain_seasonal_ar!(model, Fl)
    unconstrain_seasonal_ma!(model, Fl)
    unconstrain_mean!(model)
    unconstrain_variance!(model, "sigma2_η")
    return nothing
end

function fill_model_system!(model::SARIMA)
    update_T_terms!(model)
    update_R_terms!(model)
    if model.include_mean
        model.system.d = get_constrained_value(model, "mean")
    end
    model.system.Q[1] = get_constrained_value(model, "sigma2_η")
    return nothing
end

function fill_model_filter!(filter::KalmanFilter, model::SARIMA)
    SARIMA_exact_initialization!(filter, model)
    return nothing
end

function reinstantiate(model::SARIMA, y::Vector{Fl}) where Fl
    return SARIMA(y;
                order = (model.order.p, model.order.d, model.order.q),
                seasonal_order = (model.order.P, model.order.D, model.order.Q, model.order.s),
                include_mean = model.include_mean
            )
end

has_exogenous(::SARIMA) = false

function model_name(model::SARIMA)
    p, d, q = model.order.p, model.order.d, model.order.q
    P, D, Q, s = model.order.P, model.order.D, model.order.Q, model.order.s
    str = "SARIMA($p, $d, $q)x($P, $D, $Q, $s)"
    if model.include_mean
        str *= " with non-zero mean"
    end
    return str
end


text:
using Distributions

mutable struct RandomWalkAsset <: AbstractAsset
    last_rate::Float64

    # Random Walk Parameters
    α::Float64 # Drift
    ϕ::Float64 # Phi of the random walk (1=stationary)
    σ::Float64 # Standard deviation of a step

    name::String
end

mutable struct FixedAsset <: AbstractAsset
    last_rate::Float64

    rate::Float64

    name::String
    function FixedAsset(rate::Float64; name::String="asset")
        new(NaN, rate, name) # , Array{BidLimitOrder, 1}()
    end
    function FixedAsset(name::String="asset")
        new(NaN, 0.0, name) # , Array{BidLimitOrder, 1}()
    end
end

"Get cashflow of the RandomWalkAsset.
The cashflow would be interest in case of currency
and dividend in case of stock"
function get_rate(asset::RandomWalkAsset)
    return get_random_walk_step(
        asset.α,
        asset.ϕ,
        asset.σ,
        asset.last_interest
    )
end

function get_random_walk_step(α::Float64, ϕ::Float64, σ::Float64, x_prev::Float64)
    μ = 0
    u = rand(Normal(μ, σ))
    return α + ϕ * x_prev + u
end

function get_rate(asset::FixedAsset)
    return asset.rate
end