text:
#=

 Unit tests from libieeep1788 for interval boolean operations
 (Original author: Marco Nehmeier)
 converted into portable ITL format by Oliver Heimlich.

 Copyright 2013-2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015-2017 Oliver Heimlich (oheim@posteo.de)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal_is_empty_test" begin
    @test isempty(∅)
    @test isempty(Interval(-Inf, Inf)) == false
    @test isempty(Interval(1.0, 2.0)) == false
    @test isempty(Interval(-1.0, 2.0)) == false
    @test isempty(Interval(-3.0, -2.0)) == false
    @test isempty(Interval(-Inf, 2.0)) == false
    @test isempty(Interval(-Inf, 0.0)) == false
    @test isempty(Interval(-Inf, -0.0)) == false
    @test isempty(Interval(0.0, Inf)) == false
    @test isempty(Interval(-0.0, Inf)) == false
    @test isempty(Interval(-0.0, 0.0)) == false
    @test isempty(Interval(0.0, -0.0)) == false
    @test isempty(Interval(0.0, 0.0)) == false
    @test isempty(Interval(-0.0, -0.0)) == false
end

@testset "minimal_is_empty_dec_test" begin
    @test isempty(nai()) == false
    @test isempty(DecoratedInterval(∅, trv))
    @test isempty(DecoratedInterval(Interval(-Inf, Inf), def)) == false
    @test isempty(DecoratedInterval(Interval(1.0, 2.0), com)) == false
    @test isempty(DecoratedInterval(Interval(-1.0, 2.0), trv)) == false
    @test isempty(DecoratedInterval(Interval(-3.0, -2.0), dac)) == false
    @test isempty(DecoratedInterval(Interval(-Inf, 2.0), trv)) == false
    @test isempty(DecoratedInterval(Interval(-Inf, 0.0), trv)) == false
    @test isempty(DecoratedInterval(Interval(-Inf, -0.0), trv)) == false
    @test isempty(DecoratedInterval(Interval(0.0, Inf), def)) == false
    @test isempty(DecoratedInterval(Interval(-0.0, Inf), trv)) == false
    @test isempty(DecoratedInterval(Interval(-0.0, 0.0), com)) == false
    @test isempty(DecoratedInterval(Interval(0.0, -0.0), trv)) == false
    @test isempty(DecoratedInterval(Interval(0.0, 0.0), trv)) == false
    @test isempty(DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
end

@testset "minimal_is_entire_test" begin
    @test isentire(∅) == false
    @test isentire(Interval(-Inf, Inf))
    @test isentire(Interval(1.0, 2.0)) == false
    @test isentire(Interval(-1.0, 2.0)) == false
    @test isentire(Interval(-3.0, -2.0)) == false
    @test isentire(Interval(-Inf, 2.0)) == false
    @test isentire(Interval(-Inf, 0.0)) == false
    @test isentire(Interval(-Inf, -0.0)) == false
    @test isentire(Interval(0.0, Inf)) == false
    @test isentire(Interval(-0.0, Inf)) == false
    @test isentire(Interval(-0.0, 0.0)) == false
    @test isentire(Interval(0.0, -0.0)) == false
    @test isentire(Interval(0.0, 0.0)) == false
    @test isentire(Interval(-0.0, -0.0)) == false
end

@testset "minimal_is_entire_dec_test" begin
    @test isentire(nai()) == false
    @test isentire(DecoratedInterval(∅, trv)) == false
    @test isentire(DecoratedInterval(Interval(-Inf, Inf), trv))
    @test isentire(DecoratedInterval(Interval(-Inf, Inf), def))
    @test isentire(DecoratedInterval(Interval(-Inf, Inf), dac))
    @test isentire(DecoratedInterval(Interval(1.0, 2.0), com)) == false
    @test isentire(DecoratedInterval(Interval(-1.0, 2.0), trv)) == false
    @test isentire(DecoratedInterval(Interval(-3.0, -2.0), dac)) == false
    @test isentire(DecoratedInterval(Interval(-Inf, 2.0), trv)) == false
    @test isentire(DecoratedInterval(Interval(-Inf, 0.0), trv)) == false
    @test isentire(DecoratedInterval(Interval(-Inf, -0.0), trv)) == false
    @test isentire(DecoratedInterval(Interval(0.0, Inf), def)) == false
    @test isentire(DecoratedInterval(Interval(-0.0, Inf), trv)) == false
    @test isentire(DecoratedInterval(Interval(-0.0, 0.0), com)) == false
    @test isentire(DecoratedInterval(Interval(0.0, -0.0), trv)) == false
    @test isentire(DecoratedInterval(Interval(0.0, 0.0), trv)) == false
    @test isentire(DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
end

@testset "minimal_is_nai_dec_test" begin
    @test isnai(nai())
    @test isnai(DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test isnai(DecoratedInterval(Interval(-Inf, Inf), def)) == false
    @test isnai(DecoratedInterval(Interval(-Inf, Inf), dac)) == false
    @test isnai(DecoratedInterval(Interval(1.0, 2.0), com)) == false
    @test isnai(DecoratedInterval(Interval(-1.0, 2.0), trv)) == false
    @test isnai(DecoratedInterval(Interval(-3.0, -2.0), dac)) == false
    @test isnai(DecoratedInterval(Interval(-Inf, 2.0), trv)) == false
    @test isnai(DecoratedInterval(Interval(-Inf, 0.0), trv)) == false
    @test isnai(DecoratedInterval(Interval(-Inf, -0.0), trv)) == false
    @test isnai(DecoratedInterval(Interval(0.0, Inf), def)) == false
    @test isnai(DecoratedInterval(Interval(-0.0, Inf), trv)) == false
    @test isnai(DecoratedInterval(Interval(-0.0, 0.0), com)) == false
    @test isnai(DecoratedInterval(Interval(0.0, -0.0), trv)) == false
    @test isnai(DecoratedInterval(Interval(0.0, 0.0), trv)) == false
    @test isnai(DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
end

@testset "minimal_equal_test" begin
    @test ==(Interval(1.0, 2.0), Interval(1.0, 2.0))
    @test ==(Interval(1.0, 2.1), Interval(1.0, 2.0)) == false
    @test ==(∅, ∅)
    @test ==(∅, Interval(1.0, 2.0)) == false
    @test ==(Interval(-Inf, Inf), Interval(-Inf, Inf))
    @test ==(Interval(1.0, 2.4), Interval(-Inf, Inf)) == false
    @test ==(Interval(1.0, Inf), Interval(1.0, Inf))
    @test ==(Interval(1.0, 2.4), Interval(1.0, Inf)) == false
    @test ==(Interval(-Inf, 2.0), Interval(-Inf, 2.0))
    @test ==(Interval(-Inf, 2.4), Interval(-Inf, 2.0)) == false
    @test ==(Interval(-2.0, 0.0), Interval(-2.0, 0.0))
    @test ==(Interval(-0.0, 2.0), Interval(0.0, 2.0))
    @test ==(Interval(-0.0, -0.0), Interval(0.0, 0.0))
    @test ==(Interval(-0.0, 0.0), Interval(0.0, 0.0))
    @test ==(Interval(0.0, -0.0), Interval(0.0, 0.0))
end

@testset "minimal_equal_dec_test" begin
    @test ==(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(1.0, 2.0), trv))
    @test ==(DecoratedInterval(Interval(1.0, 2.1), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test ==(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test ==(DecoratedInterval(∅, trv), nai()) == false
    @test ==(nai(), DecoratedInterval(∅, trv)) == false
    @test ==(nai(), nai()) == false
    @test ==(DecoratedInterval(∅, trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test ==(nai(), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test ==(DecoratedInterval(Interval(-Inf, Inf), def), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ==(DecoratedInterval(Interval(1.0, 2.4), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test ==(DecoratedInterval(Interval(1.0, Inf), trv), DecoratedInterval(Interval(1.0, Inf), trv))
    @test ==(DecoratedInterval(Interval(1.0, 2.4), def), DecoratedInterval(Interval(1.0, Inf), trv)) == false
    @test ==(DecoratedInterval(Interval(-Inf, 2.0), trv), DecoratedInterval(Interval(-Inf, 2.0), trv))
    @test ==(DecoratedInterval(Interval(-Inf, 2.4), def), DecoratedInterval(Interval(-Inf, 2.0), trv)) == false
    @test ==(DecoratedInterval(Interval(-2.0, 0.0), trv), DecoratedInterval(Interval(-2.0, 0.0), trv))
    @test ==(DecoratedInterval(Interval(-0.0, 2.0), def), DecoratedInterval(Interval(0.0, 2.0), trv))
    @test ==(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ==(DecoratedInterval(Interval(-0.0, 0.0), def), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ==(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv))
end

@testset "minimal_subset_test" begin
    @test ⊆(∅, ∅)
    @test ⊆(∅, Interval(0.0, 4.0))
    @test ⊆(∅, Interval(-0.0, 4.0))
    @test ⊆(∅, Interval(-0.1, 1.0))
    @test ⊆(∅, Interval(-0.1, 0.0))
    @test ⊆(∅, Interval(-0.1, -0.0))
    @test ⊆(∅, Interval(-Inf, Inf))
    @test ⊆(Interval(0.0, 4.0), ∅) == false
    @test ⊆(Interval(-0.0, 4.0), ∅) == false
    @test ⊆(Interval(-0.1, 1.0), ∅) == false
    @test ⊆(Interval(-Inf, Inf), ∅) == false
    @test ⊆(Interval(0.0, 4.0), Interval(-Inf, Inf))
    @test ⊆(Interval(-0.0, 4.0), Interval(-Inf, Inf))
    @test ⊆(Interval(-0.1, 1.0), Interval(-Inf, Inf))
    @test ⊆(Interval(-Inf, Inf), Interval(-Inf, Inf))
    @test ⊆(Interval(1.0, 2.0), Interval(1.0, 2.0))
    @test ⊆(Interval(1.0, 2.0), Interval(0.0, 4.0))
    @test ⊆(Interval(1.0, 2.0), Interval(-0.0, 4.0))
    @test ⊆(Interval(0.1, 0.2), Interval(0.0, 4.0))
    @test ⊆(Interval(0.1, 0.2), Interval(-0.0, 4.0))
    @test ⊆(Interval(-0.1, -0.1), Interval(-4.0, 3.4))
    @test ⊆(Interval(0.0, 0.0), Interval(-0.0, -0.0))
    @test ⊆(Interval(-0.0, -0.0), Interval(0.0, 0.0))
    @test ⊆(Interval(-0.0, 0.0), Interval(0.0, 0.0))
    @test ⊆(Interval(-0.0, 0.0), Interval(0.0, -0.0))
    @test ⊆(Interval(0.0, -0.0), Interval(0.0, 0.0))
    @test ⊆(Interval(0.0, -0.0), Interval(-0.0, 0.0))
end

@testset "minimal_subset_dec_test" begin
    @test ⊆(nai(), nai()) == false
    @test ⊆(nai(), DecoratedInterval(∅, trv)) == false
    @test ⊆(nai(), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test ⊆(DecoratedInterval(∅, trv), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test ⊆(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-0.0, 4.0), def))
    @test ⊆(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-0.1, 1.0), trv))
    @test ⊆(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-0.1, 0.0), trv))
    @test ⊆(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-0.1, -0.0), trv))
    @test ⊆(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⊆(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(∅, trv)) == false
    @test ⊆(DecoratedInterval(Interval(-0.0, 4.0), def), DecoratedInterval(∅, trv)) == false
    @test ⊆(DecoratedInterval(Interval(-0.1, 1.0), trv), DecoratedInterval(∅, trv)) == false
    @test ⊆(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(∅, trv)) == false
    @test ⊆(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⊆(DecoratedInterval(Interval(-0.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⊆(DecoratedInterval(Interval(-0.1, 1.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⊆(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⊆(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(1.0, 2.0), trv))
    @test ⊆(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test ⊆(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(-0.0, 4.0), def))
    @test ⊆(DecoratedInterval(Interval(0.1, 0.2), trv), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test ⊆(DecoratedInterval(Interval(0.1, 0.2), trv), DecoratedInterval(Interval(-0.0, 4.0), def))
    @test ⊆(DecoratedInterval(Interval(-0.1, -0.1), trv), DecoratedInterval(Interval(-4.0, 3.4), trv))
    @test ⊆(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv))
    @test ⊆(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), def))
    @test ⊆(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ⊆(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, -0.0), trv))
    @test ⊆(DecoratedInterval(Interval(0.0, -0.0), def), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ⊆(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv))
end

@testset "minimal_less_test" begin
    @test ≤(∅, ∅)
    @test ≤(Interval(1.0, 2.0), ∅) == false
    @test ≤(∅, Interval(1.0, 2.0)) == false
    @test ≤(Interval(-Inf, Inf), Interval(-Inf, Inf))
    @test ≤(Interval(1.0, 2.0), Interval(-Inf, Inf)) == false
    @test ≤(Interval(0.0, 2.0), Interval(-Inf, Inf)) == false
    @test ≤(Interval(-0.0, 2.0), Interval(-Inf, Inf)) == false
    @test ≤(Interval(-Inf, Inf), Interval(1.0, 2.0)) == false
    @test ≤(Interval(-Inf, Inf), Interval(0.0, 2.0)) == false
    @test ≤(Interval(-Inf, Inf), Interval(-0.0, 2.0)) == false
    @test ≤(Interval(0.0, 2.0), Interval(0.0, 2.0))
    @test ≤(Interval(0.0, 2.0), Interval(-0.0, 2.0))
    @test ≤(Interval(0.0, 2.0), Interval(1.0, 2.0))
    @test ≤(Interval(-0.0, 2.0), Interval(1.0, 2.0))
    @test ≤(Interval(1.0, 2.0), Interval(1.0, 2.0))
    @test ≤(Interval(1.0, 2.0), Interval(3.0, 4.0))
    @test ≤(Interval(1.0, 3.5), Interval(3.0, 4.0))
    @test ≤(Interval(1.0, 4.0), Interval(3.0, 4.0))
    @test ≤(Interval(-2.0, -1.0), Interval(-2.0, -1.0))
    @test ≤(Interval(-3.0, -1.5), Interval(-2.0, -1.0))
    @test ≤(Interval(0.0, 0.0), Interval(-0.0, -0.0))
    @test ≤(Interval(-0.0, -0.0), Interval(0.0, 0.0))
    @test ≤(Interval(-0.0, 0.0), Interval(0.0, 0.0))
    @test ≤(Interval(-0.0, 0.0), Interval(0.0, -0.0))
    @test ≤(Interval(0.0, -0.0), Interval(0.0, 0.0))
    @test ≤(Interval(0.0, -0.0), Interval(-0.0, 0.0))
end

@testset "minimal_less_dec_test" begin
    @test ≤(nai(), nai()) == false
    @test ≤(DecoratedInterval(∅, trv), nai()) == false
    @test ≤(DecoratedInterval(Interval(1.0, 2.0), trv), nai()) == false
    @test ≤(nai(), DecoratedInterval(Interval(1.0, 2.0), def)) == false
    @test ≤(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test ≤(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(∅, trv)) == false
    @test ≤(DecoratedInterval(∅, trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test ≤(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ≤(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test ≤(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test ≤(DecoratedInterval(Interval(-0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test ≤(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test ≤(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(0.0, 2.0), def)) == false
    @test ≤(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-0.0, 2.0), trv)) == false
    @test ≤(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(0.0, 2.0), trv))
    @test ≤(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(-0.0, 2.0), trv))
    @test ≤(DecoratedInterval(Interval(0.0, 2.0), def), DecoratedInterval(Interval(1.0, 2.0), def))
    @test ≤(DecoratedInterval(Interval(-0.0, 2.0), trv), DecoratedInterval(Interval(1.0, 2.0), trv))
    @test ≤(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(1.0, 2.0), trv))
    @test ≤(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), def))
    @test ≤(DecoratedInterval(Interval(1.0, 3.5), trv), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test ≤(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test ≤(DecoratedInterval(Interval(-2.0, -1.0), trv), DecoratedInterval(Interval(-2.0, -1.0), trv))
    @test ≤(DecoratedInterval(Interval(-3.0, -1.5), trv), DecoratedInterval(Interval(-2.0, -1.0), trv))
    @test ≤(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv))
    @test ≤(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), def))
    @test ≤(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ≤(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, -0.0), trv))
    @test ≤(DecoratedInterval(Interval(0.0, -0.0), def), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test ≤(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv))
end

@testset "minimal_precedes_test" begin
    @test precedes(∅, Interval(3.0, 4.0))
    @test precedes(Interval(3.0, 4.0), ∅)
    @test precedes(∅, ∅)
    @test precedes(Interval(1.0, 2.0), Interval(-Inf, Inf)) == false
    @test precedes(Interval(0.0, 2.0), Interval(-Inf, Inf)) == false
    @test precedes(Interval(-0.0, 2.0), Interval(-Inf, Inf)) == false
    @test precedes(Interval(-Inf, Inf), Interval(1.0, 2.0)) == false
    @test precedes(Interval(-Inf, Inf), Interval(-Inf, Inf)) == false
    @test precedes(Interval(1.0, 2.0), Interval(3.0, 4.0))
    @test precedes(Interval(1.0, 3.0), Interval(3.0, 4.0))
    @test precedes(Interval(-3.0, -1.0), Interval(-1.0, 0.0))
    @test precedes(Interval(-3.0, -1.0), Interval(-1.0, -0.0))
    @test precedes(Interval(1.0, 3.5), Interval(3.0, 4.0)) == false
    @test precedes(Interval(1.0, 4.0), Interval(3.0, 4.0)) == false
    @test precedes(Interval(-3.0, -0.1), Interval(-1.0, 0.0)) == false
    @test precedes(Interval(0.0, 0.0), Interval(-0.0, -0.0))
    @test precedes(Interval(-0.0, -0.0), Interval(0.0, 0.0))
    @test precedes(Interval(-0.0, 0.0), Interval(0.0, 0.0))
    @test precedes(Interval(-0.0, 0.0), Interval(0.0, -0.0))
    @test precedes(Interval(0.0, -0.0), Interval(0.0, 0.0))
    @test precedes(Interval(0.0, -0.0), Interval(-0.0, 0.0))
end

@testset "minimal_precedes_dec_test" begin
    @test precedes(nai(), DecoratedInterval(Interval(3.0, 4.0), def)) == false
    @test precedes(DecoratedInterval(Interval(3.0, 4.0), trv), nai()) == false
    @test precedes(nai(), DecoratedInterval(∅, trv)) == false
    @test precedes(nai(), nai()) == false
    @test precedes(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), def))
    @test precedes(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(∅, trv))
    @test precedes(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test precedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test precedes(DecoratedInterval(Interval(0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test precedes(DecoratedInterval(Interval(-0.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test precedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test precedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test precedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test precedes(DecoratedInterval(Interval(1.0, 3.0), trv), DecoratedInterval(Interval(3.0, 4.0), def))
    @test precedes(DecoratedInterval(Interval(-3.0, -1.0), def), DecoratedInterval(Interval(-1.0, 0.0), trv))
    @test precedes(DecoratedInterval(Interval(-3.0, -1.0), trv), DecoratedInterval(Interval(-1.0, -0.0), trv))
    @test precedes(DecoratedInterval(Interval(1.0, 3.5), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) == false
    @test precedes(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) == false
    @test precedes(DecoratedInterval(Interval(-3.0, -0.1), trv), DecoratedInterval(Interval(-1.0, 0.0), trv)) == false
    @test precedes(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv))
    @test precedes(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), def))
    @test precedes(DecoratedInterval(Interval(-0.0, 0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test precedes(DecoratedInterval(Interval(-0.0, 0.0), def), DecoratedInterval(Interval(0.0, -0.0), trv))
    @test precedes(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(0.0, 0.0), trv))
    @test precedes(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv))
end

@testset "minimal_interior_test" begin
    @test ⪽(∅, ∅)
    @test isinterior(∅, ∅)
    @test ⪽(∅, Interval(0.0, 4.0))
    @test isinterior(∅, Interval(0.0, 4.0))
    @test ⪽(Interval(0.0, 4.0), ∅) == false
    @test isinterior(Interval(0.0, 4.0), ∅) == false
    @test ⪽(Interval(-Inf, Inf), Interval(-Inf, Inf))
    @test isinterior(Interval(-Inf, Inf), Interval(-Inf, Inf))
    @test ⪽(Interval(0.0, 4.0), Interval(-Inf, Inf))
    @test isinterior(Interval(0.0, 4.0), Interval(-Inf, Inf))
    @test ⪽(∅, Interval(-Inf, Inf))
    @test isinterior(∅, Interval(-Inf, Inf))
    @test ⪽(Interval(-Inf, Inf), Interval(0.0, 4.0)) == false
    @test isinterior(Interval(-Inf, Inf), Interval(0.0, 4.0)) == false
    @test ⪽(Interval(0.0, 4.0), Interval(0.0, 4.0)) == false
    @test isinterior(Interval(0.0, 4.0), Interval(0.0, 4.0)) == false
    @test ⪽(Interval(1.0, 2.0), Interval(0.0, 4.0))
    @test isinterior(Interval(1.0, 2.0), Interval(0.0, 4.0))
    @test ⪽(Interval(-2.0, 2.0), Interval(-2.0, 4.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-2.0, 4.0)) == false
    @test ⪽(Interval(-0.0, -0.0), Interval(-2.0, 4.0))
    @test isinterior(Interval(-0.0, -0.0), Interval(-2.0, 4.0))
    @test ⪽(Interval(0.0, 0.0), Interval(-2.0, 4.0))
    @test isinterior(Interval(0.0, 0.0), Interval(-2.0, 4.0))
    @test ⪽(Interval(0.0, 0.0), Interval(-0.0, -0.0)) == false
    @test isinterior(Interval(0.0, 0.0), Interval(-0.0, -0.0)) == false
    @test ⪽(Interval(0.0, 4.4), Interval(0.0, 4.0)) == false
    @test isinterior(Interval(0.0, 4.4), Interval(0.0, 4.0)) == false
    @test ⪽(Interval(-1.0, -1.0), Interval(0.0, 4.0)) == false
    @test isinterior(Interval(-1.0, -1.0), Interval(0.0, 4.0)) == false
    @test ⪽(Interval(2.0, 2.0), Interval(-2.0, -1.0)) == false
    @test isinterior(Interval(2.0, 2.0), Interval(-2.0, -1.0)) == false
end

@testset "minimal_interior_dec_test" begin
    @test ⪽(nai(), nai()) == false
    @test isinterior(nai(), nai()) == false
    @test ⪽(nai(), DecoratedInterval(∅, trv)) == false
    @test isinterior(nai(), DecoratedInterval(∅, trv)) == false
    @test ⪽(nai(), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test isinterior(nai(), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test ⪽(DecoratedInterval(Interval(0.0, 4.0), def), nai()) == false
    @test isinterior(DecoratedInterval(Interval(0.0, 4.0), def), nai()) == false
    @test ⪽(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test isinterior(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test ⪽(DecoratedInterval(∅, trv), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test isinterior(DecoratedInterval(∅, trv), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test ⪽(DecoratedInterval(Interval(0.0, 4.0), def), DecoratedInterval(∅, trv)) == false
    @test isinterior(DecoratedInterval(Interval(0.0, 4.0), def), DecoratedInterval(∅, trv)) == false
    @test ⪽(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test isinterior(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⪽(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test isinterior(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⪽(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test isinterior(DecoratedInterval(∅, trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test ⪽(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test isinterior(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test ⪽(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test isinterior(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test ⪽(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test isinterior(DecoratedInterval(Interval(1.0, 2.0), def), DecoratedInterval(Interval(0.0, 4.0), trv))
    @test ⪽(DecoratedInterval(Interval(-2.0, 2.0), trv), DecoratedInterval(Interval(-2.0, 4.0), def)) == false
    @test isinterior(DecoratedInterval(Interval(-2.0, 2.0), trv), DecoratedInterval(Interval(-2.0, 4.0), def)) == false
    @test ⪽(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(-2.0, 4.0), trv))
    @test isinterior(DecoratedInterval(Interval(-0.0, -0.0), trv), DecoratedInterval(Interval(-2.0, 4.0), trv))
    @test ⪽(DecoratedInterval(Interval(0.0, 0.0), def), DecoratedInterval(Interval(-2.0, 4.0), trv))
    @test isinterior(DecoratedInterval(Interval(0.0, 0.0), def), DecoratedInterval(Interval(-2.0, 4.0), trv))
    @test ⪽(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
    @test isinterior(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
    @test ⪽(DecoratedInterval(Interval(0.0, 4.4), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test isinterior(DecoratedInterval(Interval(0.0, 4.4), trv), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test ⪽(DecoratedInterval(Interval(-1.0, -1.0), trv), DecoratedInterval(Interval(0.0, 4.0), def)) == false
    @test isinterior(DecoratedInterval(Interval(-1.0, -1.0), trv), DecoratedInterval(Interval(0.0, 4.0), def)) == false
    @test ⪽(DecoratedInterval(Interval(2.0, 2.0), def), DecoratedInterval(Interval(-2.0, -1.0), trv)) == false
    @test isinterior(DecoratedInterval(Interval(2.0, 2.0), def), DecoratedInterval(Interval(-2.0, -1.0), trv)) == false
end

@testset "minimal_strictly_less_test" begin
    @test <(∅, ∅)
    @test <(Interval(1.0, 2.0), ∅) == false
    @test <(∅, Interval(1.0, 2.0)) == false
    @test <(Interval(-Inf, Inf), Interval(-Inf, Inf))
    @test <(Interval(1.0, 2.0), Interval(-Inf, Inf)) == false
    @test <(Interval(-Inf, Inf), Interval(1.0, 2.0)) == false
    @test <(Interval(1.0, 2.0), Interval(1.0, 2.0)) == false
    @test <(Interval(1.0, 2.0), Interval(3.0, 4.0))
    @test <(Interval(1.0, 3.5), Interval(3.0, 4.0))
    @test <(Interval(1.0, 4.0), Interval(3.0, 4.0)) == false
    @test <(Interval(0.0, 4.0), Interval(0.0, 4.0)) == false
    @test <(Interval(-0.0, 4.0), Interval(0.0, 4.0)) == false
    @test <(Interval(-2.0, -1.0), Interval(-2.0, -1.0)) == false
    @test <(Interval(-3.0, -1.5), Interval(-2.0, -1.0))
end

@testset "minimal_strictly_less_dec_test" begin
    @test <(nai(), nai()) == false
    @test <(DecoratedInterval(∅, trv), nai()) == false
    @test <(DecoratedInterval(Interval(1.0, 2.0), trv), nai()) == false
    @test <(nai(), DecoratedInterval(Interval(1.0, 2.0), def)) == false
    @test <(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test <(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(∅, trv)) == false
    @test <(DecoratedInterval(∅, trv), DecoratedInterval(Interval(1.0, 2.0), def)) == false
    @test <(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv))
    @test <(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test <(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test <(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test <(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test <(DecoratedInterval(Interval(1.0, 3.5), def), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test <(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), def)) == false
    @test <(DecoratedInterval(Interval(0.0, 4.0), trv), DecoratedInterval(Interval(0.0, 4.0), def)) == false
    @test <(DecoratedInterval(Interval(-0.0, 4.0), def), DecoratedInterval(Interval(0.0, 4.0), trv)) == false
    @test <(DecoratedInterval(Interval(-2.0, -1.0), def), DecoratedInterval(Interval(-2.0, -1.0), def)) == false
    @test <(DecoratedInterval(Interval(-3.0, -1.5), trv), DecoratedInterval(Interval(-2.0, -1.0), trv))
end

@testset "minimal_strictly_precedes_test" begin
    @test strictprecedes(∅, Interval(3.0, 4.0))
    @test strictprecedes(Interval(3.0, 4.0), ∅)
    @test strictprecedes(∅, ∅)
    @test strictprecedes(Interval(1.0, 2.0), Interval(-Inf, Inf)) == false
    @test strictprecedes(Interval(-Inf, Inf), Interval(1.0, 2.0)) == false
    @test strictprecedes(Interval(-Inf, Inf), Interval(-Inf, Inf)) == false
    @test strictprecedes(Interval(1.0, 2.0), Interval(3.0, 4.0))
    @test strictprecedes(Interval(1.0, 3.0), Interval(3.0, 4.0)) == false
    @test strictprecedes(Interval(-3.0, -1.0), Interval(-1.0, 0.0)) == false
    @test strictprecedes(Interval(-3.0, -0.0), Interval(0.0, 1.0)) == false
    @test strictprecedes(Interval(-3.0, 0.0), Interval(-0.0, 1.0)) == false
    @test strictprecedes(Interval(1.0, 3.5), Interval(3.0, 4.0)) == false
    @test strictprecedes(Interval(1.0, 4.0), Interval(3.0, 4.0)) == false
    @test strictprecedes(Interval(-3.0, -0.1), Interval(-1.0, 0.0)) == false
end

@testset "minimal_strictly_precedes_dec_test" begin
    @test strictprecedes(nai(), DecoratedInterval(Interval(3.0, 4.0), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(3.0, 4.0), def), nai()) == false
    @test strictprecedes(nai(), DecoratedInterval(∅, trv)) == false
    @test strictprecedes(nai(), nai()) == false
    @test strictprecedes(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test strictprecedes(DecoratedInterval(Interval(3.0, 4.0), def), DecoratedInterval(∅, trv))
    @test strictprecedes(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test strictprecedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 2.0), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(1.0, 2.0), trv), DecoratedInterval(Interval(3.0, 4.0), trv))
    @test strictprecedes(DecoratedInterval(Interval(1.0, 3.0), def), DecoratedInterval(Interval(3.0, 4.0), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(-3.0, -1.0), trv), DecoratedInterval(Interval(-1.0, 0.0), def)) == false
    @test strictprecedes(DecoratedInterval(Interval(-3.0, -0.0), def), DecoratedInterval(Interval(0.0, 1.0), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(-3.0, 0.0), trv), DecoratedInterval(Interval(-0.0, 1.0), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(1.0, 3.5), trv), DecoratedInterval(Interval(3.0, 4.0), trv)) == false
    @test strictprecedes(DecoratedInterval(Interval(1.0, 4.0), trv), DecoratedInterval(Interval(3.0, 4.0), def)) == false
    @test strictprecedes(DecoratedInterval(Interval(-3.0, -0.1), trv), DecoratedInterval(Interval(-1.0, 0.0), trv)) == false
end

@testset "minimal_disjoint_test" begin
    @test isdisjoint(∅, Interval(3.0, 4.0))
    @test isdisjoint(Interval(3.0, 4.0), ∅)
    @test isdisjoint(∅, ∅)
    @test isdisjoint(Interval(3.0, 4.0), Interval(1.0, 2.0))
    @test isdisjoint(Interval(0.0, 0.0), Interval(-0.0, -0.0)) == false
    @test isdisjoint(Interval(0.0, -0.0), Interval(-0.0, 0.0)) == false
    @test isdisjoint(Interval(3.0, 4.0), Interval(1.0, 7.0)) == false
    @test isdisjoint(Interval(3.0, 4.0), Interval(-Inf, Inf)) == false
    @test isdisjoint(Interval(-Inf, Inf), Interval(1.0, 7.0)) == false
    @test isdisjoint(Interval(-Inf, Inf), Interval(-Inf, Inf)) == false
end

@testset "minimal_disjoint_dec_test" begin
    @test isdisjoint(nai(), DecoratedInterval(Interval(3.0, 4.0), def)) == false
    @test isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), nai()) == false
    @test isdisjoint(DecoratedInterval(∅, trv), nai()) == false
    @test isdisjoint(nai(), nai()) == false
    @test isdisjoint(DecoratedInterval(∅, trv), DecoratedInterval(Interval(3.0, 4.0), def))
    @test isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(∅, trv))
    @test isdisjoint(DecoratedInterval(∅, trv), DecoratedInterval(∅, trv))
    @test isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(Interval(1.0, 2.0), def))
    @test isdisjoint(DecoratedInterval(Interval(0.0, 0.0), trv), DecoratedInterval(Interval(-0.0, -0.0), trv)) == false
    @test isdisjoint(DecoratedInterval(Interval(0.0, -0.0), trv), DecoratedInterval(Interval(-0.0, 0.0), trv)) == false
    @test isdisjoint(DecoratedInterval(Interval(3.0, 4.0), def), DecoratedInterval(Interval(1.0, 7.0), def)) == false
    @test isdisjoint(DecoratedInterval(Interval(3.0, 4.0), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
    @test isdisjoint(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(1.0, 7.0), trv)) == false
    @test isdisjoint(DecoratedInterval(Interval(-Inf, Inf), trv), DecoratedInterval(Interval(-Inf, Inf), trv)) == false
end
# FactCheck.exitstatus()


text:
using Pkg# hideall
Pkg.activate("_literate/EX-airfoil/Project.toml")
Pkg.update()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

# **Main author**: [Ashrya Agrawal](https://github.com/ashryaagr).
#
# ## Getting started
# Here we use the [UCI "Airfoil Self-Noise" dataset](http://archive.ics.uci.edu/ml/datasets/Airfoil+Self-Noise)
# ### Loading and  preparing the data

using MLJ
using PrettyPrinting
import DataFrames
import Statistics
using CSV
using PyPlot
ioff() # hide
using HTTP
using StableRNGs

MLJ.color_off() # hide

req = HTTP.get("https://raw.githubusercontent.com/rupakc/UCI-Data-Analysis/master/Airfoil%20Dataset/airfoil_self_noise.dat");

df = CSV.read(req.body, DataFrames.DataFrame; header=[
                   "Frequency","Attack_Angle","Chord+Length",
                   "Free_Velocity","Suction_Side","Scaled_Sound"
                   ]
              );
df[1:5, :] |> pretty

# inspect the schema:

schema(df)

# unpack into the data and labels:

y, X = unpack(df, ==(:Scaled_Sound));

# Now we Standardize the features using the transformer Standardizer()

X = MLJ.transform(fit!(machine(Standardizer(), X)), X);

# Partition into train and test set

train, test = partition(collect(eachindex(y)), 0.7, shuffle=true, rng=StableRNG(612));

# Let's first see which models are compatible with the scientific type and machine type of our data

for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end

# Note that if we coerce `X.Frequency` to `Continuous`, many more models are available:

coerce!(X, :Frequency=>Continuous)

for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end

# ## DecisionTreeRegressor
#
# We will first try out DecisionTreeRegressor:

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree

dcrm = machine(DecisionTreeRegressor(), X, y)

fit!(dcrm, rows=train)
pred_dcrm = predict(dcrm, rows=test);

# Now you can call a loss function to assess the performance on test set.

rms(pred_dcrm, y[test])

# ## RandomForestRegressor
#
# Now let's try out RandomForestRegressor:

RandomForestRegressor = @load RandomForestRegressor pkg=DecisionTree
rfr = RandomForestRegressor()

rfr_m = machine(rfr, X, y);

# train on the rows corresponding to train

fit!(rfr_m, rows=train);

# predict values on the rows corresponding to test

pred_rfr = predict(rfr_m, rows=test);
rms(pred_rfr, y[test])

# Unsurprisingly, the RandomForestRegressor does a better job.
#
# Can we do even better? Yeah, we can!! We can make use of Model Tuning.
#
# ## Tuning
#
# In case you are new to model tuning using MLJ, refer [lab5](https://alan-turing-institute.github.io/DataScienceTutorials.jl/isl/lab-5/) and [model-tuning](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/model-tuning/)
#
# Range of values for parameters should be specified to do hyperparameter tuning

r_maxD = range(rfr, :n_trees, lower=9, upper=15)
r_samF = range(rfr, :sampling_fraction, lower=0.6, upper=0.8)
r = [r_maxD, r_samF];

# Now we specify how the tuning should be done. Let's just specify a coarse grid tuning with cross validation and instantiate a tuned model:

tuning = Grid(resolution=7)
resampling = CV(nfolds=6)

tm = TunedModel(model=rfr, tuning=tuning,
                resampling=resampling, ranges=r, measure=rms)

rfr_tm = machine(tm, X, y);

# train on the rows corresponding to train

fit!(rfr_tm, rows=train);

# predict values on the rows corresponding to test

pred_rfr_tm = predict(rfr_tm, rows=test);
rms(pred_rfr_tm, y[test])

# That was great! We have further improved the accuracy
#
# Now to retrieve best model, You can use

fitted_params(rfr_tm).best_model

# Now we can investigate the tuning by using report.
# Let's plot a heatmap of the measurements:

r = report(rfr_tm)
res = r.plotting

md = res.parameter_values[:,1]
mcw = res.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(md, mcw, res.measurements)

xlabel("Number of trees", fontsize=14)
ylabel("Sampling fraction", fontsize=14)
xticks(9:1:15, fontsize=12)
yticks(fontsize=12)
plt.savefig(joinpath(@OUTPUT, "airfoil_heatmap.svg")) # hide

# \figalt{Hyperparameter heatmap}{airfoil_heatmap.svg}

PyPlot.close_figs() # hide

ext: jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


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

JLLWrappers.@generate_main_file_header("GR")
JLLWrappers.@generate_main_file("GR", UUID("d2c73de3-f751-5644-a686-071e5b155ba9"))
end  # module GR_jll


text:
# Multivariate Normal distribution

using StatsFuns
using LinearAlgebra
using Random
import Base
using LoopVectorization

import MeasureTheory: logdensity
import MeasureTheory: MvNormal

using StrideArrays
using StaticArrays

const LowerCholesky{T} = Cholesky{T, <:LowerTriangular}
const UpperCholesky{T} = Cholesky{T, <:UpperTriangular}

@inline function logdet_pos(C::Cholesky)
    logdet_pos(getfield(C, :factors))
end

@inline function logdet_pos(A::Union{UpperTriangular{T},LowerTriangular{T}}) where {T}
    # ∑ᵢ log(aᵢ * 2^bᵢ) = log(∏ᵢ aᵢ) + log2 * ∑ᵢ bᵢ

    prod_ai = one(real(T))
    sum_bi = zero(real(T))

    @turbo for i = 1:ArrayInterface.size(A,1)
        diag_i = A.data[i, i]
        ai = significand(diag_i)
        bi = exponent(diag_i)
        prod_ai *= ai
        sum_bi += bi
    end

    return log(prod_ai) + logtwo * sum_bi
end

###############################################################################
# MvNormal(σ)

function logdensity(UnrollQ::MayUnroll, d::MvNormal{(:σ,), <:Tuple{<:Cholesky}}, y)
    logdensity_mvnormal_σ(UnrollQ, KnownSize(getfield(d.σ, :factors)), y)
end

logdensity(d::MvNormal{(:σ,), <:Tuple{<:Cholesky}}, y) = logdensity_mvnormal_σ(getfield(d.σ, :factors), y)

logdensity_mvnormal_σ(UL::Union{UpperTriangular, LowerTriangular}, y::AbstractVector) = logdensity_mvnormal_σ(KnownSize(UL), y)

# TODO: Can we do this statically without a generated function?
@generated function logdensity_mvnormal_σ(U::KnownSize{Tuple{k,k}, <:UpperTriangular}, y::AbstractVector) where {k}
    UnrollQ = k < 20 ? Unroll() : NoUnroll()

    quote
        $(Expr(:meta,:inline))
        logdensity_mvnormal_σ($UnrollQ, U, y)
    end
end

@generated function logdensity_mvnormal_σ(::UnrollQ, U::KnownSize{Tuple{k,k}, <:UpperTriangular}, y::AbstractVector{T}) where {k,T, UnrollQ<:MayUnroll}
    log2π = log(big(2) * π)

    # Solve `y = σz` for `z`. We need this only as a way to calculate `z ⋅ z`

    header = quote
        $(Expr(:meta,:inline))
        U = U.value
        Udata = U.data
        z_dot_z = zero($T)
    end

    body = if UnrollQ <: Unroll
        quote
            z = StrideArray{$T}(undef, ($(StaticInt(k)),))

            @inbounds @fastmath for j ∈ 1:$k
                tmp = zero($T)
                for i = 1:(j-1)
                    tmp += Udata[i, j] * z[i]
                end
                zj = (y[j] - tmp) / Udata[j, j]
                z_dot_z += zj^2
                z[j] = zj
            end
        end
    else
        quote
            @inbounds begin # `@fastmath` needs to be moved inside the `@nexprs`
                Base.Cartesian.@nexprs $k j -> begin
                    tmp = zero($T)
                    z_j = zero($T)
                    Base.Cartesian.@nexprs j - 1 i -> begin
                        @fastmath tmp += Udata[i, j] * z_i
                    end
                    @fastmath z_j = (y[j] - tmp) / Udata[j, j]
                    @fastmath z_dot_z += z_j^2
                end
            end
        end
    end

    footer = quote
        $(T(-k / 2 * log2π)) - logdet_pos(U) - z_dot_z / 2
    end

    return quote
        $(header.args...)
        $(body.args...)
        $(footer.args...)
    end

end

###############################################################################
# MvNormal(ω)

function logdensity(UnrollQ::MayUnroll, d::MvNormal{(:ω,), <:Tuple{<:Cholesky}}, y)
    logdensity_mvnormal_ω(UnrollQ, KnownSize(getfield(d.ω, :factors)), y)
end

logdensity(d::MvNormal{(:ω,), <:Tuple{<:Cholesky}}, y) = logdensity_mvnormal_ω(getfield(d.ω, :factors), y)

logdensity_mvnormal_ω(UL::Union{UpperTriangular, LowerTriangular}, y::AbstractVector) = logdensity_mvnormal_ω(KnownSize(UL), y)

# TODO: Can we do this statically without a generated function?
@generated function logdensity_mvnormal_ω(U::KnownSize{Tuple{k,k}, <:UpperTriangular}, y::AbstractVector) where {k}
    UnrollQ = k < 20 ? Unroll() : NoUnroll()

    quote
        $(Expr(:meta,:inline))
        logdensity_mvnormal_ω($UnrollQ, U, y)
    end
end

@generated function logdensity_mvnormal_ω(::UnrollQ, U::KnownSize{Tuple{k,k}, <:UpperTriangular}, y::AbstractVector{T}) where {k,T, UnrollQ<:MayUnroll}
    log2π = log(big(2) * π)

    header = quote
        $(Expr(:meta,:inline))
        U = U.value
        Udata = U.data
        # if z = Lᵗy, the logdensity depends on `det(U)` and `z ⋅ z`. So we find `z`
        z_dot_z = zero(T)
    end

    body = if UnrollQ <: Unroll
        quote
            @inbounds begin
                Base.Cartesian.@nexprs $k j -> begin
                    zj = zero(T)
                    Base.Cartesian.@nexprs j i -> begin
                        zj += Udata[i, j] * y[i]
                    end
                    z_dot_z += zj^2
                end
            end
        end
    else
        quote
            @fastmath for j ∈ 1:$k
                zj = zero(T)
                for i ∈ 1:j
                    @inbounds zj += Udata[i, j] * y[i]
                end
                z_dot_z += zj^2
            end
        end
    end

    footer = quote
        $(T(-k / 2 * log2π)) + logdet_pos(U) - z_dot_z / 2
    end

    return quote
        $(header.args...)
        $(body.args...)
        $(footer.args...)
    end
end

@inline function logdensity_mvnormal_ω(U::KnownSize{Tuple{nothing, nothing}, <:UpperTriangular}, y::AbstractVector{T}) where {T}
    k = first(size(U))
    @assert length(y) == k

    # if z = Lᵗy, the logdensity depends on `det(U)` and `z ⋅ z`. So we find `z`
    z_dot_z = zero(T)
    for j ∈ 1:k
        zj = zero(T)
        for i ∈ 1:j
            @inbounds zj += U[i, j] * y[i]
        end
        z_dot_z += zj^2
    end

    -k / 2 * log2π + logdet_pos(U) - z_dot_z / 2
end

function logdensity(d::MvNormal{(:μ, :σ)}, y::AbstractArray{T}) where {T}
    x = StrideArray{T}(undef, ArrayInterface.size(y))
    @inbounds for j in eachindex(y)
        x[j] = y[j] - d.μ[j]
    end
    GC.@preserve x logdensity(MvNormal(σ = d.σ), x)
end
using StrideArrays, StaticArrays, LoopVectorization, LinearAlgebra


text:
function smolyak_evaluate(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  unique_multi_index = sort(unique(multi_index))
  unique_orders = m_i(unique_multi_index).-1

  # Below we do the following things:

  #   Generate the polynomial terms for each order
  #   Generate the unique polynomial terms introduced at each higher order
  #   Combine the polynomial terms to construct the Smolyak polynomial

  # Here we construct the base polynomials

  base_polynomials = Array{Array{T,2}}(undef,length(unique_orders))
  for i = 1:length(unique_orders)
    base_polynomials[i] = chebyshev_polynomial(unique_orders[i],node)
  end

  # Compute the unique polynomial terms from the base polynomials

  unique_base_polynomials = Array{Array{T,2}}(undef,length(unique_orders))
  for i = length(unique_orders):-1:2
    unique_base_polynomials[i] = base_polynomials[i][:,size(base_polynomials[i-1],2)+1:end]
  end
  unique_base_polynomials[1] = base_polynomials[1]

  # Construct the first row of the interplation matrix

  polynomials = Array{T,1}(undef,length(weights))

  # Iterate over nodes, doing the above three steps at each iteration

  l = 1
  @inbounds for j = 1:size(multi_index,1)
    new_polynomials = unique_base_polynomials[multi_index[j,1]][1,:]
    for i = 2:size(multi_index,2)
      new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]][i,:])
    end
    m = length(new_polynomials)
    polynomials[l:l+m-1] = new_polynomials
    l += m
  end

  estimate = zero(T)
  for i = 1:length(polynomials)
    estimate += polynomials[i]*weights[i]
  end

  return estimate

end

function smolyak_evaluate(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  node = copy(node)

  if size(domain,2) != length(node)
    error("domain is inconsistent with the number of dimensions")
  end

  d = length(node)
  for i = 1:d
    node[i] = normalize_node(node[i],domain[:,i])
  end

  estimate = smolyak_evaluate(weights,node,multi_index)

  return estimate

end

###########################################################################

function smolyak_evaluate(weights::Array{T,1},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_evaluate(weights,x,multi_index)

  end

  return goo

end

function smolyak_evaluate(weights::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_evaluate(weights,x,multi_index,domain)

  end

  return goo

end


text:
module MCMCBenchmarkSuite
  using Reexport
  @reexport using MCMCBenchmarks

end


text:
"""
    deepcopy_module(m::Module)

Returns a newly constructed Module (`m2`) that is a complete deep-copy of
module `m`, containing deep-copies of all the `name`s in `m`.

All functions contained in `m` are deep-copied, so that the returned module
`m2` contains new generic function objects that have all the same methods as
the functions from `m`, and all code in the new methods will have been modified
so that any references to names in `m` will now be references to names in `m2`.

# Example
```julia-repl
julia> module M;  x = 1;  f() = (global x += 1);  end
Main.M

julia> M2 = deepcopy_module(M)
Main.M

julia> M2.f()  # This only modifies the copy of `x` in M2.
2

julia> M2.x, M.x
(2, 1)

julia> nameof(M2)  # Inside M2, it still believes it's module M
:M
```
"""
function deepcopy_module(m::Module)
    m2 = Module(nameof(m), false)

    @eval m2 include(fname::AbstractString) = Main.Base.include($m2, fname)
    @eval m2 eval(x) = Core.eval($m2, x)

    setparent!(m2, parentmodule(m))
    import_names_into(m2, m)
    return m2
end

usings(m::Module) =
    ccall(:jl_module_usings, Array{Module,1}, (Any,), m)

_full_using_name(modu) = Expr(:., fullname(modu)...)
_import_name_expr(val) = _import_name_expr(parentmodule(val), nameof(val))
_import_name_expr(modu::Module, name::Symbol) = Expr(:(:), _full_using_name(modu), Expr(:., name))

function import_names_into(destmodule, srcmodule)
    for m in usings(srcmodule)
        Core.eval(destmodule, Expr(:using, _full_using_name(m)))
    end
    imports = names(srcmodule, imported=true)
    # don't copy itself into itself; don't copy `include`, defined manually instead
    excluded_names = (nameof(srcmodule), Symbol("#include"), :include)
    imports = setdiff(imports, excluded_names)
    for n in imports
        srcval = Core.eval(srcmodule, n)
        if srcval isa Module
            Core.eval(destmodule, Expr(:import, _full_using_name(srcval)))
        else
            #@show srcmodule, n, typeof(srcval)
            try
                # Try importing from the original source module
                Core.eval(destmodule, Expr(:import, _import_name_expr(srcval)))
            catch
                # Otherwise, import from srcmodule directly
                Core.eval(destmodule, Expr(:import, _import_name_expr(srcmodule, n)))
            end
        end
    end
    ns = setdiff(names(srcmodule, all=true), imports)
    ns = setdiff(ns, excluded_names)
    # Do definitions bottom up to get `Symbol("#foo")` kw defs last
    for n in reverse(ns)
        srcval = Core.eval(srcmodule, n)
        #@show destmodule, srcmodule, n, typeof(srcval)
        # If this name has already been defined by previous steps, skip it
        try Core.eval(destmodule, n) ; continue; catch end
        deepcopy_value(destmodule, srcmodule, n, srcval)
    end
end

# Make new names w/ a copy of the value. For functions, make a new function object.
deepcopy_value(destmodule, srcmodule, name, value) = Core.eval(destmodule, :($name = $(deepcopy(value))))
function deepcopy_value(destmodule, srcmodule, name, value::Function)
    deepcopy_function(destmodule, value)
end
function deepcopy_value(destmodule, srcmodule, name, value::Module)
    m2 = deepcopy_module(value)
    setparent!(m2, destmodule)
    m2
end
# Shallow-copy these as references for now
deepcopy_value(destmodule, srcmodule, name,
               value::Union{Base.Docs.Binding, IdDict}) = Core.eval(destmodule, :($name = $(value)))

# Manually implement "missing" julia Base function to allow setting the parent of a module.
function setparent!(m::Module, p::Module)
    unsafe_store!(Ptr{_Module2}(pointer_from_objref(m)),
                   _Module2(nameof(m), p), 1)
    m
end
# NOTE: This struct must be kept up-to-date with Julia's `_jl_module_t`!
struct _Module2
    name::Symbol
    parent::Module
end


text:
@testset "code generation (relying on a device)" begin

############################################################################################

@testset "SASS" begin

@testset "basic reflection" begin
    @eval sass_valid_kernel() = nothing
    @eval sass_invalid_kernel() = 1

    @test CUDAnative.code_sass(devnull, sass_valid_kernel, Tuple{}) == nothing
    @test_throws CUDAnative.CompilerError CUDAnative.code_sass(devnull, sass_invalid_kernel, Tuple{})
end

end

############################################################################################

end


text:
# data source is: https://pomber.github.io/covid19/timeseries.json
# data is updated thrice daily from Johns Hopkins numbers

using Pkg
pkg"activate .."
using JSON, Statistics, Plots, Dates, LsqFit, TimeSeries
data = JSON.parsefile("../data/covid_data.json");

# Aggregate statistics
function unpack_country(country)
    country_data = data[country]
    confirmed = [datum["confirmed"] for datum in country_data]
    dates = Date.([datum["date"] for datum in country_data])
    diffs = diff(confirmed)
    return country_data, confirmed, dates, diffs
end

function plot_country_aggregates(country)
    country_data, confirmed, dates, diffs = unpack_country(country)
    p = plot(dates[2:end], confirmed[2:end], lw = 3, legend = :topleft, title = "$country Aggregate", label = "Confirmed Cases")
    plot!(p, dates[2:end], diffs, lw = 3, label = "Day-to-Day Difference", xticks = dates[2:10:end])
    savefig(p, "~/Desktop/$(country)_aggregate.jpeg")
end

plot_country_aggregates("Belgium")
# plot_country_aggregates("Canada")
# plot_country_aggregates("France") # because this gives us a longer view

# Regression modeling
function fit_country(country, model, p_guess)
    country_data, confirmed, dates, diffs = unpack_country(country)
    t = 1:length(dates)
    fit = curve_fit(model, t, confirmed, p_guess) # return this
end

# Exponential growth (y(t) = Ae^(kt))
@. model(t, p) = p[1]*exp(p[2]*t)
fr_fit = fit_country("France", model, [0.5, 0.5])
@show coef(fr_fit) # ends up being y(t) = 3.64e^(0.13t)
@show stderror(fr_fit) # stderror(fr_fit) = [0.5807248710582464, 0.002419476241731794]

be_fit = fit_country("Belgium", model, [0.5, 0.5])
@show coef(be_fit) # ends up being y(t) = 0.07*e^(0.17t)
@show stderror(be_fit) # stderror(be_fit) = [0.009881708341035871, 0.0018586282749802605]

prc_fit = fit_country("China", model, [0.5, 0.5])
@show coef(prc_fit) # ends up being y(t) = 30700e^0.017t... bit different than the others...
@show stderror(prc_fit) # stderror(prc_fit) = [3123.3555378369188, 0.0020020903719526564]

# Logistic modeling
# The logistic growth equation is (Carrying Capacity)/(1 + e^(=growth rate*(t - t_peak)))
# Let's try filling in the gaps with stuff we know
L = 0.7*11.4*1e6 # millions
be_rate = 0.17 # from above
@. model(t, p) = L/(1 + exp(-be_rate*(t - p[1])))
be_data, be_conf, be_dates, diffs = unpack_country("Belgium")
t = 1:length(be_dates)
be_logistic_fit = curve_fit(model, t, be_conf, [80.])

@show coef(be_logistic_fit) # 107.08398211493504
@show stderror(be_logistic_fit) # stderror(be_logistic_fit) = [0.03196080987817913]

# Plot regression
p = coef(be_logistic_fit)[1]
restricted_data = model.(1:length(be_dates), Ref(p))
pl = plot(be_dates, restricted_data, lw = 3, label = "Model", title = "Belgium So Far", legend = :topleft)
plot!(be_dates, be_conf, lw = 3, label = "Data", xticks = be_dates[1:10:end])
savefig(pl, "~/Desktop/Belgium_restricted.jpeg")

# Belgium Overall
full_data = model.(1:200, Ref(p))
dates = collect(range(be_dates[1], step = Day(1), length = 200))
pl = plot(dates, full_data, lw = 3, label = "Model", legend = :topleft, title = "Belgium Model Full Timescale", xticks = dates[1:50:end])
peak = dates[1] + Day(floor(p))
vline!(pl, [peak], label = "$peak", lw = 3)
savefig(pl, "~/Desktop/Belgium_model.jpeg")


text:
# compute L2 norm
function sqnorm(x)
  return sum(abs2, x)
end

# compute binary crossentropy loss
function crossentropy(p, l)

  # add epsilon perturbation for floating-point stability
  return -1*(l*log(p+eps(float(eltype(p))))+(1-l)*log(1-p+eps(float(eltype(p)))))
end

# get a batch of indicies for NN training
function get_batches(indices, batchsize, batch)
  if batch*batchsize > length(indices)
    if (batch-1)*batchsize+1 > length(indices)
      error("batchsize too large")
    else
      return indices[1+batchsize*(batch-1):end]
    end
  else
    return indices[1+batchsize*(batch-1):batch*batchsize]
  end
end

# construct data set for NN training based on prototypical probability distributions (cases 1-3)
function get_dataset_proto(p_range, distribution, samples)
  dataset = zeros(eltype(p_range[1]), (3, length(p_range)*length(samples)))
  count = 1
  for i in 1:length(p_range)
    for j in 1:length(samples)
      dataset[:, count] = [j, distribution(samples[j], p_range[i]), i]
      count += 1
    end
  end
  return dataset
end

# construct data set for NN training based on a given probability distributions
function get_dataset(p_range, distribution, samples)
  dataset_ul, n_unique = get_dataset_unlabelled(p_range, distribution, samples)

  dataset = zeros(eltype(p_range[1]), (3, n_unique))
  count = 1
  len_p_train = 0
  for i in 1:length(dataset_ul)
    data = dataset_ul[i]
    p = p_range[i]

    for j in 1:length(data)
      freq = data[j][2]
      if freq != zero(eltype(freq))
        dataset[:, count] = [data[j][1], freq, i*one(eltype(freq))]
        count += 1
      end
    end
  end
  return dataset
end

# construct data set for NN training based on a given probability distributions, where labels are not yet included
function get_dataset_unlabelled(p_range, distribution, samples)
  data, counter = prepare_dataset(p_range[1], distribution, samples)
  datasets = [data]
  for indxp in 2:length(p_range)
    data,count = prepare_dataset(p_range[indxp], distribution, samples)
    push!(datasets, data)
    counter += count
  end
  return datasets, Int(counter)
end

# prepare list of samples and count number of unique samples
function prepare_dataset(p, distribution, samples)
  counter = zero(eltype(length(samples)))
  dist = distribution(samples[1][1],p)
  if dist != zero(eltype(dist))
    counter += one(eltype(counter))
  end

  data = [[samples[1][1], dist]]
  for indx in 2:length(samples)
    dist = distribution(samples[indx][1], p)
    if dist != zero(eltype(dist))
      counter += one(eltype(counter))
    end
    push!(data, [samples[indx][1], dist])
  end
  return data, counter
end

# construct the training data set for applying SL with NNs, where we modify the underlying probability distributions to combat deviations due to finite sample statistics (applicable for data generated from Boltzmann distributions)
function get_modified_dataset_train_SL(distr, samples, p_range, n_samples; gs_index=1)
  new_dataset_train_SL = zeros(eltype(p_range[1]), (3, gs_index+length(samples)))

  for gs_indx in 1:gs_index
    new_dataset_train_SL[:, gs_indx]=[samples[gs_indx], distr(samples[gs_indx], p_range[1]), 1]
  end
  num_samples = n_samples

  for indx in collect(1:gs_index)
    sample = samples[indx]
    p = length(p_range)

    new_dataset_train_SL[:, indx+gs_index]=[sample, distr(sample,p_range[end]), p]
  end

  for indx in collect(gs_index+1:length(samples))
    sample = samples[indx]
    prob = Int(round(distr(sample,p_range[end])*n_samples))

    # add an instance of each sample which has not appeared in the training set to the set of samples associated with the largest sampled temperature value
    if prob == zero(eltype(p_range[1]))
      prob = one(eltype(p_range[1]))
      num_samples += one(eltype(n_samples))
    end
    p = length(p_range)

    new_dataset_train_SL[:, indx+gs_index]=[sample, prob, p]
  end

  new_dataset_train_SL[2, gs_index+1:end] = new_dataset_train_SL[2, gs_index+1:end]./num_samples

  return new_dataset_train_SL
end

# construct the training data set for applying SL with NNs, where we do not modify the underlying probability distributions (applicable for data generated from Boltzmann distributions)
function get_unmodified_dataset_train_SL(distr, samples, p_range, n_samples)
  new_dataset_train_SL = zeros(eltype(p_range[1]), (3, 1+length(samples)))
  new_dataset_train_SL[:, 1] = [samples[1], distr(samples[1], p_range[1]), 1]

  for indx in collect(1:length(samples))
    sample = samples[indx]
    prob = distr(sample, p_range[end])
    p = length(p_range)
    new_dataset_train_SL[:, indx+1] = [sample, prob, p]
  end

  new_dataset_train_SL[2, 2:end] = new_dataset_train_SL[2, 2:end]

  return new_dataset_train_SL
end

# construct the training data set for applying SL with NNs (general probability distributions)
function get_training_data_SL(dataset_SL, p_range, p_max, p_min)
	train = [dataset_SL[:, 1]]
	for i in 2:size(dataset_SL)[2]
		ele = dataset_SL[:, i]
		if p_range[Int(ele[3])] <= p_min || p_range[Int(ele[3])] >= p_max
			push!(train, ele)
		end
	end
	return hcat(train...)
end

# compute mean and standard deviation of input features given a data set (used for standardization)
function get_dataset_stats(dataset, inputs, n_points)
  mu = zeros(eltype(dataset[2, 1]), size(inputs)[1])
  musq = zeros(eltype(dataset[2, 1]), size(inputs)[1])
  for indx in collect(1:size(dataset)[2])
    mu .+= dataset[2, indx]*inputs[Int(dataset[1, indx])]
    musq .+= dataset[2, indx]*inputs[Int(dataset[1, indx])]^2
  end
  mu = mu./n_points
  musq = musq./n_points
  return mu, sqrt.(musq.-mu.^2)
end


text:
using Statistics
using ProgressMeter
using MATLAB
include("atn_module.jl")

N = 10_000          # number of nodes
m = 5               # number of active links per active node
mu = 0.1            # probability I -> S
tmax = 1000         # maximum duration of the simulation
Ni = Int(0.01 * N)  # initial number of infects
gamma = -2.1      # exponent of the power law
NTrials = 1       # number of trials for each value of lambda
TWindow = 500.0      # time window for evaluating steady state
TMinConvergence = 500

LambdaMin = 0.0
LambdaMax = 1.0
LambdaStep = 0.05
Lambda = LambdaMin:LambdaStep:LambdaMax
LambdaSize = length(Lambda)

EtaMin = 0.0
EtaMax = 15
EtaStep = 0.5
Eta = EtaMin:EtaStep:EtaMax
EtaSize = length(Eta)
AvgInfectedLambda = zeros(LambdaSize)
Ni = Int(0.01 * N)  # initial number of infects
gamma = -2.1       # exponent of the power law
TWindow = 500       # time window for evaluating steady state
TMinConvergence = 500
ITrend = zeros(Int, 0)
ActMin = 0.001      # minimum activation probability
ActMax = 1
AgentAct =
    (
        (ActMax .^ (gamma + 1) - ActMin .^ (gamma + 1)) .* rand(N) .+
        ones(N) * ActMin .^ (gamma + 1)
    ) .^ (1.0 / (gamma + 1))
# AvgDegree = 2 * m * mean(AgentAct)
# AgentState = zeros(Bool, N)
# AgentState[1:Ni] .= 1
etaS = 1
eta = 1

# grid=Reshape([(eta, lambda) for eta in Eta, lambda in Lambda],:)
# @showprogress for (ietaI, eta) in enumerate(Eta)
# @showprogress for (ilambda, lambda) in enumerate(Lambda)
lambda = 1
InfectedOverTrials = zeros(0)
InfectedOverATrial = zeros(0)

AgentState = zeros(Int, N)
AgentStateNext = zeros(Int, N)
for k = 1:100
    AgentState[k] = 1
end
# InfectedOverATrial = 0.0
for i_time = 1:tmax
    for l = 1:N
        AgentStateNext[l] = AgentState[l]
    end
    for i = 1:N
        if (AgentState[i] == 1)
            if (rand() < mu)
                AgentStateNext[i] = 0
            end
        end
        # println(i)
        if (AgentAct[i] > rand())
            for j = 1:m
                rInt = rand(1:N)
                # println(rInt)
                if (AgentState[i] == 1)
                    if ((AgentState[rInt] == 0) && (rand() < lambda))
                        AgentStateNext[rInt] = 1
                    end
                else
                    if ((AgentState[rInt] == 1) && (rand() < lambda))
                        AgentStateNext[i] = 1
                    end
                end

            end
        end
    end
    for l = 1:N
        AgentState[l] = AgentStateNext[l]
    end
    NumInf = sum(AgentState)
    push!(InfectedOverTrials, NumInf)
end
# InfectedOverATrial = InfectedOverATrial / TWindow
# push!(InfectedOverTrials, InfectedOverATrial)
# end

mat"figure();set(gcf, 'Position',  [0, 0, 2500, 1500]); hold on;"
# # mat"axis([ -0,1,-0,10 ])"
# # mat"imagesc([0,1],[0,10],$AvgInfectedEtaLambda)"
# # mat"colormap jet"
#
mat"plot(1:3000,$InfectedOverTrials,'.-','MarkerSize',20)"


text:
using Test, MsgPack2, Serialization

function can_round_trip(value, T,
                        expected_typed_output = value,
                        expected_any_output = value)
    bytes = pack(value)
    return isequal(unpack(bytes, T), expected_typed_output) &&
           isequal(unpack(bytes), expected_any_output)
end

# IntegerType

@test can_round_trip(30, UInt8)
@test can_round_trip(-30, Int8)
@test can_round_trip(typemax(UInt8), UInt8)
@test can_round_trip(typemax(UInt16), UInt16)
@test can_round_trip(typemax(UInt32), UInt32)
@test can_round_trip(typemax(UInt64), UInt64)
@test can_round_trip(typemax(Int8), Int8)
@test can_round_trip(typemin(Int8), Int8)
@test can_round_trip(typemax(Int16), Int16)
@test can_round_trip(typemin(Int16), Int16)
@test can_round_trip(typemax(Int32), Int32)
@test can_round_trip(typemin(Int32), Int32)
@test can_round_trip(typemax(Int64), Int64)
@test can_round_trip(typemin(Int64), Int64)

# NilType

@test can_round_trip(nothing, Nothing)
@test can_round_trip(nothing, Missing, missing, nothing)
@test can_round_trip(missing, Missing, missing, nothing)
@test can_round_trip(missing, Nothing, nothing, nothing)

# BooleanType

@test can_round_trip(true, Bool)
@test can_round_trip(false, Bool)

# FloatType

@test can_round_trip(rand(Float32), Float32)
@test can_round_trip(rand(Float64), Float64)

# StringType

@test can_round_trip(join(rand(Char, 9)), String)
@test can_round_trip(join(rand(Char, typemax(UInt8) - 1)), String)
@test can_round_trip(join(rand(Char, typemax(UInt8) + 1)), String)
@test can_round_trip(join(rand(Char, typemax(UInt16) + 1)), String)

# BinaryType

struct ByteVec
    bytes::Vector{UInt8}
end

Base.:(==)(a::ByteVec, b::ByteVec) = a.bytes == b.bytes

MsgPack2.msgpack_type(::Type{ByteVec}) = MsgPack2.BinaryType()
MsgPack2.to_msgpack(::MsgPack2.BinaryType, x::ByteVec) = x.bytes
MsgPack2.from_msgpack(::Type{ByteVec}, bytes::Vector{UInt8}) = ByteVec(bytes)

bytes = rand(UInt8, typemax(UInt8))
@test can_round_trip(ByteVec(bytes), ByteVec, ByteVec(bytes), bytes)
bytes = rand(UInt8, typemax(UInt8) + 1)
@test can_round_trip(ByteVec(bytes), ByteVec, ByteVec(bytes), bytes)
bytes = rand(UInt8, typemax(UInt16) + 1)
@test can_round_trip(ByteVec(bytes), ByteVec, ByteVec(bytes), bytes)

# ArrayType

arr = [30, -30, typemax(UInt8), typemax(UInt16), typemax(UInt32), typemax(UInt64),
       typemax(Int8), typemin(Int8), typemax(Int16), typemin(Int16), typemax(Int32),
       typemin(Int32), typemax(Int64), typemin(Int64), nothing, rand(Float32),
       rand(Float64), true, false, join(rand(Char, typemax(UInt8) - 1))]

push!(arr, deepcopy(arr))

@test can_round_trip(arr, Vector{Any})
for x in arr
    T = Vector{typeof(x)}
    A = MsgPack2.ArrayView{typeof(x)}
    a = fill(x, 9)
    @test can_round_trip(a, T)
    @test can_round_trip(a, A)
    a = fill(x, typemax(UInt8) - 1)
    @test can_round_trip(a, T)
    @test can_round_trip(a, A)
    a = fill(x, typemax(UInt8) + 1)
    @test can_round_trip(a, T)
    @test can_round_trip(a, A)
    a = fill(x, typemax(UInt16) + 1)
    @test can_round_trip(a, T)
    @test can_round_trip(a, A)
end

tup = (arr...,)
@test can_round_trip(tup, Tuple, tup, arr)

set = Set(arr)
@test can_round_trip(set, Set, set, collect(set))

# MapType

dict = Dict(zip(arr, reverse(arr)))

@test can_round_trip(dict, Dict{Any,Any})
for v in values(dict)
    T = Dict{Int,typeof(v)}
    M = MsgPack2.MapView{Int,typeof(v)}
    d = Dict(zip(1:9, fill(v, 9)))
    m = unpack(pack(d), M)
    @test length(m) == length(d)
    for i in 1:length(d)
        @test m[i] == d[i]
        @test get(() -> nothing, m, i + 1) == get(() -> nothing, d, i + 1)
        @test get(m, i + 1, nothing) == get(d, i + 1, nothing)
    end
    @test Set(collect(m)) == Set(collect(d))
    @test Set(collect(keys(m))) == Set(collect(keys(d)))
    @test Set(collect(values(m))) == Set(collect(values(d)))
    @test can_round_trip(d, T)
    @test can_round_trip(d, M)
    d = Dict(zip(1:(typemax(UInt8) - 1), fill(v, typemax(UInt8) - 1)))
    @test can_round_trip(d, T)
    @test can_round_trip(d, M)
    d = Dict(zip(1:(typemax(UInt8) + 1), fill(v, typemax(UInt8) + 1)))
    @test can_round_trip(d, T)
    @test can_round_trip(d, M)
    d = Dict(zip(1:(typemax(UInt16) + 1), fill(v, typemax(UInt16) + 1)))
    @test can_round_trip(d, T)
    @test can_round_trip(d, M)
end

namedtup = (x = arr[1], y = arr[2], z = arr[3], others = arr[4:end])
namedtup_dict = Dict((string(k) => v for (k, v) in pairs(namedtup))...)
@test can_round_trip(namedtup, NamedTuple{keys(namedtup),<:Tuple}, namedtup, namedtup_dict)

arr = [dict]
@test can_round_trip(arr, MsgPack2.ArrayView{MsgPack2.MapView{Any,Any}})

# ImmutableStructType

struct Bar{T}
    a::T
    b::T
end

Base.:(==)(a::Bar, b::Bar) = a.a == b.a && a.b == b.b

MsgPack2.construct(::Type{T}, args...) where {T<:Bar} = T(promote(args...)...)

MsgPack2.msgpack_type(::Type{<:Bar}) = MsgPack2.ImmutableStructType()

struct Foo{T,S}
    x::Union{Nothing,T}
    y::Vector{S}
    z::Bar{T}
end

Base.:(==)(a::Foo, b::Foo) = a.x == b.x && a.y == b.y && a.z == b.z

MsgPack2.msgpack_type(::Type{<:Foo}) = MsgPack2.ImmutableStructType()

foo = Foo{Int,String}(nothing, String["abc", join(rand(Char,typemax(UInt16)))],
                      Bar(rand(Int), rand(Int)))
foo_dict = Dict("x" => foo.x, "y" => foo.y, "z" => Dict("a" => foo.z.a, "b" => foo.z.b))
@test can_round_trip(foo, typeof(foo), foo, foo_dict)

foo = Foo{Float64,Char}(rand(), rand('a':'z', 100), Bar(rand(), rand()))
foo_dict = Dict("x" => foo.x, "y" => map(string, foo.y), "z" => Dict("a" => foo.z.a, "b" => foo.z.b))
@test can_round_trip(foo, typeof(foo), foo, foo_dict)

arr = [foo, foo]
@test can_round_trip(arr, MsgPack2.ArrayView{typeof(foo)}, arr, [foo_dict, foo_dict])

# MutableStructType

mutable struct MBar{T}
    a::T
    b::T
    MBar{T}() where {T} = new{T}()
end

Base.:(==)(a::MBar, b::MBar) = a.a == b.a && a.b == b.b

MsgPack2.msgpack_type(::Type{<:MBar}) = MsgPack2.MutableStructType()

mutable struct MFoo{T,S}
    x::Union{Nothing,T}
    y::Vector{S}
    z::MBar{T}
    MFoo{T,S}() where {T,S} = new{T,S}()
end

Base.:(==)(a::MFoo, b::MFoo) = a.x == b.x && a.y == b.y && a.z == b.z

MsgPack2.msgpack_type(::Type{<:MFoo}) = MsgPack2.MutableStructType()

foo = MFoo{Int,String}()
foo.x = nothing
foo.y = String["abc", join(rand(Char,typemax(UInt16)))]
foo.z = (z = MBar{Int}(); z.a = rand(Int); z.b = rand(Int); z)
foo_dict = Dict("x" => foo.x, "y" => foo.y, "z" => Dict("a" => foo.z.a, "b" => foo.z.b))
@test can_round_trip(foo, typeof(foo), foo, foo_dict)

foo = MFoo{Float64,Char}()
foo.x = rand()
foo.y = rand('a':'z', 100)
foo.z = (z = MBar{Float64}(); z.a = rand(); z.b = rand(); z)
foo_dict = Dict("x" => foo.x, "y" => map(string, foo.y), "z" => Dict("a" => foo.z.a, "b" => foo.z.b))
@test can_round_trip(foo, typeof(foo), foo, foo_dict)

arr = [foo, foo]
@test can_round_trip(arr, MsgPack2.ArrayView{typeof(foo)}, arr, [foo_dict, foo_dict])

# ExtensionType

type = 123
data = (io = IOBuffer(); serialize(io, arr); take!(io))
@test MsgPack2.Ext(type, data) === MsgPack2.Extension(type, data)
ext = MsgPack2.extserialize(type, arr)
@test ext.type == type
@test ext.data == data
ext_type, ext_arr = MsgPack2.extdeserialize(ext)
@test ext_type == ext.type
@test ext_arr == arr

ext = MsgPack2.Extension(4, [0xbb])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(-32, [0x56, 0x8d])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(-123, [0x80, 0x7c, 0x8b, 0xf8])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(111, [0x04, 0x16, 0x94, 0x13, 0x0a, 0x7d, 0x6f, 0x0c])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(79, [0x00, 0x30, 0xd5, 0x64, 0x0f, 0x8d, 0x92, 0x90,
                              0x98, 0x99, 0x14, 0x57, 0x0e, 0x8d, 0xf1, 0x3a])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(-118, [0xc7, 0x00, 0x8a])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(50, [0x62, 0x4c, 0x7c, 0x0f, 0x86, 0x04])
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(78, rand(UInt8, typemax(UInt8) - 1))
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(32, rand(UInt8, typemax(UInt8) + 1))
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(-14, rand(UInt8, typemax(UInt16) + 1))
@test can_round_trip(ext, MsgPack2.Extension)
ext = MsgPack2.Extension(84, UInt8[])
@test can_round_trip(ext, MsgPack2.Extension)


text:
using Base:Float64
using PlantBiophysics
using OrderedCollections
using MutableNamedTuples
using Test
using DataFrames

@testset "File IO" begin
    include("IO.jl")
end

@testset "Structures" begin
    include("structs.jl")
end

@testset "Temperature dependence" begin
    include("temp-dependence.jl")
end

@testset "Stomatal conductance" begin
    include("gs.jl")
end

@testset "Energy balance" begin
    include("energy_balance.jl")
end


text:
using DrWatson
@quickactivate("HyperspectraWithNeXL")
# Load the necessary libraries
using HyperspectraWithNeXL
using NeXLSpectrum
using Gadfly # , Cairo, Fontconfig,
using Unitful
# Load the HyperSpectrum from disk
lt = 0.72*4.0*18.0*3600.0/(1024*1024) # 18.0 hours on 4 detectors
hs = NeXLSpectrum.compress(HyperSpectrum(
    LinearEnergyScale(0.0,10.0),
    Dict{Symbol,Any}(
      :TakeOffAngle => deg2rad(35.0),
      :ProbeCurrent => 1.0,
      :LiveTime => lt,
      :BeamEnergy => 20.0e3,
      :Name => "Mn Nodule"
    ),
    readrplraw(joinpath(datadep"MnNodule","map[15]")),
    fov = [ 4.096u"mm", 4.096u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
))
klms = [
    n"C K-L2", #
    characteristic(n"Ag", ltransitions, 0.1)..., #
    n"Al K-L3", #
    n"Ca K-L3", n"Ca K-M3", #
    n"Ce L3-M5", #
    n"Cl K-L3", #
    n"Cr K-L3", n"Cr K-M3", n"Cr L3-M5", #
    n"Cu K-L3", n"Cu K-M3", n"Cu L3-M5", #
    n"Fe K-L3", n"Fe K-M3", n"Fe L3-M5", #
    n"S K-L3", n"S K-M3", #
    characteristic(n"Ba", ltransitions, 0.03)..., #
    n"P K-L3", #
    n"K K-L3", #
    n"Mg K-L3", #
    n"O K-L3", #
    n"Mn K-L3", n"Mn K-M3", n"Mn L3-M5", #
    n"Na K-L3", #
    n"Ni K-L3", n"Ni K-M3", n"Ni L3-M5", #
    n"Si K-L3", #
    n"Ti K-L3", n"Ti K-M3", n"Ti L3-M4", #
    n"Zn K-L3", n"Zn K-M3", n"Zn L3-M5", #
]
set_default_plot_size(8inch, 3inch)
p1 = vstack(
    plot(maxpixel(hs), klms=klms, xmax=10.0e3),
    plot(sum(hs), klms=klms, xmax=10.0e3)
)
# p1 |> PDF(joinpath(plotsdir(),"MaxPixel.pdf"), 6inch, 4inch)
p1 |> SVG(joinpath(plotsdir(), "MaxPixel.svg"), 6inch, 4inch)

p2=plot(maxpixel(hs), sum(hs), klms=klms, xmax=10.0e3, norm=ScaleSum())
# p2 |> PDF(joinpath(plotsdir(),"MaxPixel2.pdf"), 8inch, 3inch)
p2 |> SVG(joinpath(plotsdir(), "MaxPixel2.svg"), 6inch, 4inch)

mkpath(joinpath(papersdir(),"Figures", "Figure 2"))
p3=plot(maxpixel(hs), klms=klms, xmax=10.0e3)
p3 |> SVG(joinpath(papersdir(),"Figures", "Figure 2","maxpixel.svg"), 8inch, 3inch)


text:
# Simulate a Thomas cluster point process on a rectangle
# Author: H. Paul Keeler, 2019.
# Website: hpaulkeeler.com
# Repository: github.com/hpaulkeeler/posts
# For more details, see the post:
# hpaulkeeler.com/testing-the-julia-language-with-point-process-simulations/

#Note: Need the .+ for adding a scalar to an array
#Also need . for sqrt, exp, cos, sin etc and assinging scalars to arrays
#Best to use vectors instead of 1-D matrices eg x=rand(n),  NOT x=rand(n,1).

using Distributions #for random simulations
using Plots #for plotting

#Simulation window parameters
xMin=-.5;
xMax=.5;
yMin=-.5;
yMax=.5;

#Parameters for the parent and daughter point processes
lambdaParent=10;#density of parent Poisson point process
lambdaDaughter=10;#mean number of points in each cluster
sigma=0.05; #sigma for normal variables (ie random locations) of daughters

#Extended simulation windows parameters
rExt=7*sigma; #extension parameter
#for rExt, use factor of deviation sigma eg 6 or 7
xMinExt=xMin-rExt;
xMaxExt=xMax+rExt;
yMinExt=yMin-rExt;
yMaxExt=yMax+rExt;
#rectangle dimensions
xDeltaExt=xMaxExt-xMinExt;
yDeltaExt=yMaxExt-yMinExt;
areaTotalExt=xDeltaExt*yDeltaExt; #area of extended rectangle

#Simulate Poisson point process
numbPointsParent=rand(Poisson(areaTotalExt*lambdaParent)); #Poisson number of points

#x and y coordinates of Poisson points for the parent
xxParent=xMinExt.+xDeltaExt*rand(numbPointsParent);
yyParent=yMinExt.+yDeltaExt*rand(numbPointsParent);

#Simulate Poisson point process for the daughters (ie final poiint process)
numbPointsDaughter=rand(Poisson(lambdaDaughter),numbPointsParent);
numbPoints=sum(numbPointsDaughter); #total number of points

#Generate the (relative) locations in Cartesian coordinates by
#simulating independent normal variables
xx0=rand(Normal(0,sigma),numbPoints);
yy0=rand(Normal(0,sigma),numbPoints);

#replicate parent points (ie centres of disks/clusters)
xx=vcat(fill.(xxParent, numbPointsDaughter)...);
yy=vcat(fill.(yyParent, numbPointsDaughter)...);

#Shift centre of disk to (xx0,yy0)
xx=xx.+xx0;
yy=yy.+yy0;

#thin points if outside the simulation window
booleInside=((xx.>=xMin).&(xx.<=xMax).&(yy.>=yMin).&(yy.<=yMax));
#retain points inside simulation window
xx=xx[booleInside];
yy=yy[booleInside];

#Plotting
plot1=scatter(xx,yy,xlabel ="x",ylabel ="y", leg=false);
display(plot1);

max_stars_repo_name: hpkeeler/posts

max_issues_repo_name: hpkeeler/posts

max_forks_repo_name: hpkeeler/posts


text:
using Yields
using Test

include("Rate.jl")
include("bootstrap.jl")
include("RateCombination.jl")
include("SmithWilson.jl")

include("ActuaryUtilities.jl")


text:
"""
    unix_timestamp_ms(dt::Union{DateTime, ZonedDateTime}) -> Int64

Get a datetime's representation as a UNIX timestamp in milliseconds.
`DateTime`s with no time zone are assumed to be in UTC.
"""
function unix_timestamp_ms end

unix_timestamp_ms(zdt::ZonedDateTime) = floor(Int64, TimeZones.zdt2unix(zdt) * 1000)
# assume UTC because you have to assume something
unix_timestamp_ms(dt::DateTime) = unix_timestamp_ms(ZonedDateTime(dt, tz"UTC"))

"""
    unix_timestamp_ms() -> Int64

Get the current datetime's representation as a UNIX timestamp in milliseconds.
"""
unix_timestamp_ms() = unix_timestamp_ms(Dates.now(tz"UTC"))

"""
    LogEvent(message::AbstractString, datetime=Dates.now(tz"UTC"))
    LogEvent(message::AbstractString, timestamp)

Log event for submission to CloudWatch Logs.
"""
struct LogEvent
    message::String
    timestamp::Int64

    function LogEvent(message::AbstractString, timestamp::Real)
        message = String(message)

        if isempty(message)
            throw(ArgumentError("Log Event message must be non-empty"))
        end

        if sizeof(message) > MAX_EVENT_SIZE - 26
            # Truncate message and throw a warning
            idx = MAX_EVENT_SIZE - 29
            # Truncated messages include a "..."
            message = string(message[1:idx], "...")
            warn(LOGGER, "CloudWatch Log Event message cannot be more than $MAX_EVENT_SIZE bytes")
        end

        if timestamp < 0
            throw(ArgumentError("Log Event timestamp must be non-negative"))
        end

        new(message, timestamp)
    end
end

function LogEvent(message::AbstractString, dt::Union{DateTime, ZonedDateTime})
    return LogEvent(message, unix_timestamp_ms(dt))
end

LogEvent(message::AbstractString) = LogEvent(message, Dates.now(tz"UTC"))

"""
    aws_size(event::LogEvent) -> Int

Returns the size of a log event as represented by AWS, used to calculate the log batch size.

See the Amazon CloudWatch Logs documentation for [`PutLogEvents`](https://docs.aws.amazon.com/AmazonCloudWatchLogs/latest/APIReference/API_PutLogEvents.html#CWL-PutLogEvents-request-sequenceToken).
"""
aws_size(event::LogEvent) = sizeof(event.message) + 26

message(event::LogEvent) = event.message
timestamp(event::LogEvent) = event.timestamp


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libClp"], :libClp),
    LibraryProduct(prefix, ["libClpSolver"], :libClpSolver),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaOpt/ClpBuilder/releases/download/v1.16.11-1a-static"

download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ClpBuilder.v1.16.11.aarch64-linux-gnu-gcc4.tar.gz", "9ab682ea624549801bf95c1169cd1da055941d9449b53185589e2c818d003e4b"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.aarch64-linux-gnu-gcc7.tar.gz", "92bd67f6f1ab11749b57228cd24134b5e8ca3a4b9ab695291528240ce33e6d68"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.aarch64-linux-gnu-gcc8.tar.gz", "7de4e0d406560ca36d68693371dc5d1aaafc8ccc92922ab2942a3ee17ed5cebb"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ClpBuilder.v1.16.11.arm-linux-gnueabihf-gcc4.tar.gz", "84ecf9257fd7cf35e0e03f70ba91330a120c202ca668bf4ea59ae462e114479d"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.arm-linux-gnueabihf-gcc7.tar.gz", "669c8eb49f1c55ae40f25b3102794771c8b85df4a01282f4de6e03a2d551dc6c"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.arm-linux-gnueabihf-gcc8.tar.gz", "1183075c704d653ba09e3f56971972ff8f8183a7081cca2aa8e244fc98c81be5"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ClpBuilder.v1.16.11.i686-linux-gnu-gcc4.tar.gz", "5bb05730d63755d3c9a79dc8a78d09511568e62ad84ed802bc6f8077120c6b75"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.i686-linux-gnu-gcc7.tar.gz", "a193103da7eee0b9c8ff7f331297d0e71baea40ef2cc25183e426f1f6b0e9e48"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.i686-linux-gnu-gcc8.tar.gz", "5eed843e7db9459d1b47476499b313bfe9a2a020cc2fcb6cca5a6964529e12f1"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc6)) => ("$bin_prefix/ClpBuilder.v1.16.11.i686-w64-mingw32-gcc6.tar.gz", "a5c4036856d399bf06d33e422775a2e031ca906885af092690311a29c907dcee"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.i686-w64-mingw32-gcc7.tar.gz", "e008d435230fa77da0897237df5b594c8f1a8828edbf324a9e1dc54f23db77d4"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.i686-w64-mingw32-gcc8.tar.gz", "67c010fd535e371b44c5ad626107b3c42b24ae00f6f011cb8340cf1a68c8ee9b"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-apple-darwin14-gcc4.tar.gz", "07814295e81f24dd33dc1cfea349a840e086ba420808557161f3c7ada7b3cb03"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-apple-darwin14-gcc7.tar.gz", "a5259539475be446f6ade1c73513ebc8eb9d1bf7b8cd7f53430879ad3f2a8e4a"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-apple-darwin14-gcc8.tar.gz", "2ac3cd8ebc619887f2816b5f7cfdf99fe2d190ffe38eb92252d23d5ebfeb5715"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-linux-gnu-gcc4.tar.gz", "e2895e33c9193e4f547aaf8b75fe602c7506e8f87014dd9a92752c2a71cd6a63"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-linux-gnu-gcc7.tar.gz", "4c4f9c37c790a20d6e7874c684e898a012b5d07d61e125aa3a84fdf80620bcfc"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-linux-gnu-gcc8.tar.gz", "40dc02584678b264dde7be3084245a0b0efa84db4f143959fd1eca1a69c718a8"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc6)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-w64-mingw32-gcc6.tar.gz", "2ee820b182e4f878f973d66e6b8bea0a0cb9b3f9258de930bf7c94cbc80de71c"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-w64-mingw32-gcc7.tar.gz", "b4fe6a96a52ff7db307aab4cc6993f231f77ded6b39aeffa4b4c05cf11a8b770"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ClpBuilder.v1.16.11.x86_64-w64-mingw32-gcc8.tar.gz", "2711d239bcbd3b7da18520694c447923c229c451b7f75c93573518f5424a6f61"),
)

# To fix gcc4 bug in Windows
# https://sourceforge.net/p/mingw-w64/bugs/727/
this_platform = platform_key_abi()
if typeof(this_platform)==Windows && this_platform.compiler_abi.gcc_version == :gcc4
   this_platform = Windows(arch(this_platform), libc=libc(this_platform), compiler_abi=CompilerABI(:gcc6))
end

# no dynamic dependencies until Pkg3 support for binaries
dependencies = [
#    "https://github.com/JuliaOpt/OsiBuilder/releases/download/v0.107.9-1/build_OsiBuilder.v0.107.9.jl",
#     "https://github.com/JuliaOpt/CoinUtilsBuilder/releases/download/v2.10.14-1/build_CoinUtilsBuilder.v2.10.14.jl",
#     "https://github.com/JuliaOpt/COINMumpsBuilder/releases/download/v1.6.0-1/build_COINMumpsBuilder.v1.6.0.jl",
#     "https://github.com/JuliaOpt/COINMetisBuilder/releases/download/v1.3.5-1/build_COINMetisBuilder.v1.3.5.jl",
#     "https://github.com/JuliaOpt/COINLapackBuilder/releases/download/v1.5.6-1/build_COINLapackBuilder.v1.5.6.jl",
#     "https://github.com/JuliaOpt/COINBLASBuilder/releases/download/v1.4.6-1/build_COINBLASBuilder.v1.4.6.jl",
#     "https://github.com/JuliaOpt/ASLBuilder/releases/download/v3.1.0-1/build_ASLBuilder.v3.1.0.jl"
]

custom_library = false
if haskey(ENV,"JULIA_CLP_LIBRARY_PATH")
    custom_products = [LibraryProduct(ENV["JULIA_CLP_LIBRARY_PATH"],product.libnames,product.variable_name) for product in products]
    if all(satisfied(p; verbose=verbose) for p in custom_products)
        products = custom_products
        custom_library = true
    else
        error("Could not install custom libraries from $(ENV["JULIA_CLP_LIBRARY_PATH"]).\nTo fall back to BinaryProvider call delete!(ENV,\"JULIA_CLP_LIBRARY_PATH\") and run build again.")
    end
end

if !custom_library
    # Install unsatisfied or updated dependencies:
    unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)

    dl_info = choose_download(download_info, this_platform)
    if dl_info === nothing && unsatisfied
        # If we don't have a compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something even more ambitious here.
        error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
    end

    # If we have a download, and we are unsatisfied (or the version we're
    # trying to install is not itself installed) then load it up!
    if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
        # Download and install binaries
        # no dynamic dependencies until Pkg3 support for binaries
        # for dependency in reverse(dependencies)          # We do not check for already installed dependencies
        #    download(dependency,basename(dependency))
        #    evalfile(basename(dependency))
        # end
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end
 end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)# using BinaryProvider # requires BinaryProvider 0.3.0 or later


text:
# ```@setup corrections
# # the example command line input for this script
# ARGS = ["example.toml", "143", "143", "TT", "--plot"]
# ```
#
# # [Corrections for Point Sources (corrections.jl)](@id corrections)
# We need to estimate a correction to the covariance matrix due to insufficiently apodized
# point source holes in the likelihood mask, and we had generated spectra in [`signalsim`](@ref signalsim).
# Now we will estimate a smooth correction from these spectra, by calculating the analytic covariance
# and then using a Gaussian process to smoothly interpolate the deviation from 1 between the sample covariance
# and the analytic covariance.
#
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia corrections.jl example.toml 143 143 TT</code></pre>
# ```
#
# # Configuration
# As usual, we set up the various information we'll need later.

configfile, freq1, freq2, spec = ARGS

freqs = [freq1, freq2]
splits = ["1", "2"]
print("$(freqs[1]), $(freqs[2]), $(spec)\n")

## We choose to never do ET, we'll just simulate TE but swap the maps.
transposed = (spec == "ET")
if transposed
    spec = "TE"
    splits = reverse(splits)
    freqs = reverse(freqs)
end

using Plots
using Healpix
using PowerSpectra
using TOML
using UUIDs, JLD2, FileIO
using Statistics
using GaussianProcesses

include("util.jl")

config = TOML.parsefile(configfile)
nside = config["general"]["nside"]
run_name = config["general"]["name"]
simpath = joinpath(config["scratch"], "signalsims",
    "$(freq1)_$(freq2)_$(spec)")
lmax = nside2lmax(nside)
lmax_planck = min(2508, lmax)

# We'll load the existing simulated spectra from disk.

function readsims(path)
    files = [f for f in readdir(simpath; join=false) if f != "summary.jld2"]
    test_spec_ = load(joinpath(path,files[1]), "cl")
    cl_array = zeros((length(test_spec_), length(files)))
    for (i,f) in enumerate(files)
        cl_array[:,i] .= parent(load(joinpath(path,f), "cl"))
    end
    return cl_array
end

@time cl_array = readsims(simpath)
cl_mean = mean(cl_array, dims=2)
cl_var = var(cl_array, dims=2)
@save joinpath(simpath, "summary.jld2") cl_mean cl_var

signal11, th = signal_and_theory(freq1, freq1, config)
signal12, th = signal_and_theory(freq1, freq2, config)
signal21, th = signal_and_theory(freq2, freq1, config)
signal22, th = signal_and_theory(freq2, freq2, config)

# Let's take a look at the results of the simulations compared to the signal.
if "--plot" ∈ ARGS
    plot(cl_mean[1:lmax_planck+1] .* collect(0:lmax_planck).^2,
        yerr=(sqrt.(cl_var[1:lmax_planck+1])  .* collect(0:lmax_planck).^2),
        label="sims", alpha=0.5)
    plot!(signal12[spec] .* collect(0:length(signal11[spec])-1).^2, xlim=(0,2nside),
        label="injected signal")
end

# Next, we'll estimate the analytic signal-only covariance matrix. This requires the masks,
# the signal spectrum, and some zero arrays for the pixel variances.

X, Y = Symbol(spec[1]), Symbol(spec[2])
run_name = config["general"]["name"]
masktype1 = (X == :T) ? "T" : "P"
masktype2 = (Y == :T) ? "T" : "P"
mapid1 = "P$(freq1)hm$(splits[1])"
mapid2 = "P$(freq2)hm$(splits[2])"

maskfileT₁ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_maskT.fits")
maskfileP₁ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_maskP.fits")
maskfileT₂ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_maskT.fits")
maskfileP₂ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_maskP.fits")

maskT₁ = readMapFromFITS(maskfileT₁, 1, Float64)
maskP₁ = readMapFromFITS(maskfileP₁, 1, Float64)
maskT₂ = readMapFromFITS(maskfileT₂, 1, Float64)
maskP₂ = readMapFromFITS(maskfileP₂, 1, Float64);

zero_var_component = HealpixMap{Float64, RingOrder}(zeros(nside2npix(nside)))
zero_var = PolarizedHealpixMap(zero_var_component, zero_var_component, zero_var_component);

m1_signal = CovField(mapid1, maskT₁, maskP₁, zero_var)
m2_signal = CovField(mapid2, maskT₂, maskP₂, zero_var)

## convert a Vector (starting from 0) into a SpectralVector (i.e. a 0-indexed vector)
function format_signal(v::Vector{T}, nside) where T
    result = SpectralVector(zeros(T, nside2lmax(nside)+1))
    last_ind = min(length(v), length(result)) - 1
    for ℓ in 0:last_ind
        result[ℓ] = v[ℓ + 1]
    end
    return result
end

spectra = Dict{SpectrumName, SpectralVector{Float64, Vector{Float64}}}(
    (:TT, mapid1, mapid1) => format_signal(signal11["TT"], nside),
    (:TT, mapid1, mapid2) => format_signal(signal12["TT"], nside),
    (:TT, mapid2, mapid1) => format_signal(signal21["TT"], nside),
    (:TT, mapid2, mapid2) => format_signal(signal22["TT"], nside),

    (:EE, mapid1, mapid1) => format_signal(signal11["EE"], nside),
    (:EE, mapid1, mapid2) => format_signal(signal12["EE"], nside),
    (:EE, mapid2, mapid1) => format_signal(signal21["EE"], nside),
    (:EE, mapid2, mapid2) => format_signal(signal22["EE"], nside),

    (:TE, mapid1, mapid1) => format_signal(signal11["TE"], nside),
    (:TE, mapid1, mapid2) => format_signal(signal12["TE"], nside),
    (:TE, mapid2, mapid1) => format_signal(signal21["TE"], nside),
    (:TE, mapid2, mapid2) => format_signal(signal22["TE"], nside),
);

workspace_signal = CovarianceWorkspace(m1_signal, m2_signal, m1_signal, m2_signal);

@time C = coupledcov(Symbol(spec), Symbol(spec), workspace_signal, spectra)

symspec = (spec == "EE") ? :M⁺⁺ : Symbol(spec)
mask₁ = (spec[1] == 'T') ? maskT₁ : maskP₁
mask₂ = (spec[2] == 'T') ? maskT₂ : maskP₂
@time 𝐌 = mcm(symspec, mask₁, mask₂)
C_decoupled = decouple_covmat(C, 𝐌, 𝐌)

# Once a covariance matrix has been estimated, it needs to be corrected for mode-coupling
# just like the spectra.
# The covariance is possibly transposed, but we only care about the diagonal.

correction = SpectralVector(cl_var[:,1] ./ diag(parent(C_decoupled)))

# Now it's finally time to fit a smooth interpolation over the noisy correction we
# just estimated.

lmin_cut, lmax_cut = get_plic_ellrange(Symbol(spec), freq1, freq2)
correction[lmax_cut:end] .= 1.0

fit_ells = intersect((lmin_cut-20):(lmax_cut), 0:lmax)
u = correction[fit_ells] .- 1
t = 1.0 * collect(fit_ells)

# We use a really generic Gaussian kernel.

se = SE(0.0, 0.0)
m = MeanZero()
gp = GP(t,u,m,se)

using Optim
optimize!(gp)   # Optimise the hyperparameters
correction_gp = predict_y(gp, float.(0:(lmax_cut+20)))[1];
correction_gp[1:(lmin_cut-21)] .= 0.0

correction_path = joinpath(config["scratch"], "point_source_corrections")
mkpath(correction_path)

if transposed
    revspec = reverse(spec)
    correction_file = joinpath(correction_path, "$(freq1)_$(freq2)_$(revspec)_corr.dat")
else
    correction_file = joinpath(correction_path, "$(freq1)_$(freq2)_$(spec)_corr.dat")
end
writedlm(correction_file, correction_gp)

# This smooth curve will be applied to the diagonal of the covariance matrices later.
if "--plot" ∈ ARGS
    plot(correction_gp)
end

max_forks_repo_name: sgiardie/PSpipe


text:
include("tutorial.jl")
include("mathprog.jl")


text:
using Printf

isnumeric(s::AbstractString) = tryparse(Float64, s) isa Number

tests = ["1", "-121", "one", "pi", "1 + 1", "NaN", "1234567890123456789", "1234567890123456789123456789",
        "1234567890123456789123456789.0", "1.3", "1.4e10", "Inf", "1//2", "1.0 + 1.0im"]

for t in tests
    fl = isnumeric(t) ? "is" : "is not"
    @printf("%35s %s a direct numeric literal.\n", t, fl)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export SinusoidalFit
import LsqFit

"""
    SinusoidalFit(s, fs) <: Decomposition
Decompose a timeseries `s` into a **sum** `x + r` (with `x` the periodic component),
by fitting sinuisoidals with given
frequencies `fs` to the signal `s` using the package `LsqFit`. Specifically, fit
```math
s - \\bar{s} \\approx \\sum_i A_i \\sin(2\\pi f_i t + \\phi_i)
```
with ``\\bar{s}`` the mean. The fit happens on the amplitudes
and phases ``A_i, \\phi_i``. After the `decomposition` you can find these in the struct's
fields `amps, phases`. The fit is done on `s` versus `t`, so be sure that you have
transformed `t` appropriately (e.g. if `t` is "days" but your frequencies are multiples of
years, then you should give `t/365.26`).

    SinusoidalFit(fs, φ0s, A0s [, ub, lb])

The quality of the fit depends dramatically on the initial guesses for the phases
and amplitudes, `φ0s, A0s`. This second constructor
gives full control over initial phases and amplitudes, as well as
upper and lower bounds on the amplitudes `ub, lb` (also vectors of `length(fs)`).
In the first constructor `φ0 = 0, A0 = abs(-(extrema(s)...))/2, ub = Inf, lb = -Inf` for
all frequencies. (the bounds for the phases are always ± π)

**Notice**: `LsqFit` performs poorly for fitting sinusoidals and `Fourier`
should be preferred over this method if the signal given is in
multiples of the expected periods (and of course `t` is equally spaced).
"""
mutable struct SinusoidalFit{T<:AbstractFloat} <: Decomposition
    fs::Vector{T}
    φ0s::Vector{T}
    A0s::Vector{T}
    ub::Vector{T}
    lb::Vector{T}
    amps::Vector{T}
    phases::Vector{T}
end

function SinusoidalFit(s::AbstractVector{T}, fs) where {T}
    φ0s = zero(fs)
    A0s = fill(abs(-(extrema(s)...))/2, length(fs))
    ub = fill(T(Inf), length(fs))
    lb = fill(T(-Inf), length(fs))
    SinusoidalFit{T}(fs, φ0s, A0s, ub, lb, zero(fs), zero(fs))
end

SinusoidalFit(fs, φ0s, A0s, ub=fill(Inf, length(fs)), lb=fill(-Inf, length(fs))) =
SinusoidalFit{eltype(fs)}(fs, φ0s, A0s, ub, lb, zero(fs), zero(fs))

function decompose(t, s, method::SinusoidalFit{T}) where {T}
    fs = method.fs
    L = length(fs)
    function sinuisoidal(t, p)
        r = zeros(T, length(t))
        for i in 1:L
            @. r += p[i]*sin(2π * (fs[i]*t + p[i+L]))
        end
        return r
    end

    p0 = vcat(method.A0s, method.φ0s)
    ub = vcat(method.ub, fill( T(π), L))
    lb = vcat(method.lb, fill(-T(π), L))

    fit = LsqFit.curve_fit(sinuisoidal, t, s .- mean(s), p0; lower=lb, upper=ub)
    method.amps .= LsqFit.coef(fit)[1:L]
    method.phases .= LsqFit.coef(fit)[L+1:2L]
    periodic = sinuisoidal(t, LsqFit.coef(fit)) .+ mean(s)
    return periodic, s .- periodic
end


text:
function cub(x::Vector, y::Vector, style...)
    h[i] = x[i] - x[i-1]

    s[i] = 1/h[i] * ( (x[i] - xs)^3 / 6 * M[i-1] + h[i]^3 / 3 * M[i] + y[i-1] - h[i]^2 / 6 * M[i-1] * x[i] - xs + y[i] - h[i]^2 / 6  * M[i] * h[i] )

    return CubicSpline(d, c, b, y, x)
end