text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_isdigit" begin
     @test Raytracing.Interpreter.isdigit("2")
     @test Raytracing.Interpreter.isdigit("1")
     @test !Raytracing.Interpreter.isdigit("10")
     @test !Raytracing.Interpreter.isdigit("a")
     @test !Raytracing.Interpreter.isdigit("")
end

@testset "test_isdecimal" begin
     @test Raytracing.Interpreter.isdecimal("2")
     @test Raytracing.Interpreter.isdecimal("1123129")
     @test !Raytracing.Interpreter.isdecimal("10.9")
     @test !Raytracing.Interpreter.isdecimal("ag")
     @test !Raytracing.Interpreter.isdecimal("10ag")
     @test !Raytracing.Interpreter.isdecimal("")
end

@testset "test_isalpha" begin
     @test Raytracing.Interpreter.isalpha("a")
     @test Raytracing.Interpreter.isalpha("ADcLopnb__o")
     @test !Raytracing.Interpreter.isalpha("1")
     @test !Raytracing.Interpreter.isalpha("aD1sd")
     @test !Raytracing.Interpreter.isalpha("")
end

@testset "test_isalnum" begin
     @test Raytracing.Interpreter.isalnum("adf_g")
     @test Raytracing.Interpreter.isalnum("1239")
     @test Raytracing.Interpreter.isalnum("A1__Def98")
     @test !Raytracing.Interpreter.isalnum("78KO.b")
     @test !Raytracing.Interpreter.isalnum("78KO@b")
     @test !Raytracing.Interpreter.isalnum("")
end

@testset "test_close_bracket" begin
     @test Raytracing.Interpreter.close_bracket("(") == ")"
     @test Raytracing.Interpreter.close_bracket("[") == "]"
     @test Raytracing.Interpreter.close_bracket("{") == "}"
     @test Raytracing.Interpreter.close_bracket("<") == ">"
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket(" (")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket(")")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket("]")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket("}")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket(">")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket("ag ")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket("123")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket("*")
     @test_throws ArgumentError Raytracing.Interpreter.close_bracket("(1")
end

function assert_is_keyword(token::Token, keyword::KeywordEnum)
     @assert isa(token.value, KeywordToken) "Token '$(token.value)' is not a KeywordToken"
     @assert token.value.keyword == keyword "Token '$(token.value)' is not equal to keyword '$(keyword)'"
end

function assert_is_identifier(token::Token, identifier::String)
     @assert isa(token.value, IdentifierToken) "Token '$(token.value)' is not a IdentifierToken"
     @assert token.value.identifier == identifier "expecting identifier '$(identifier)' instead of '$(token.value)'"
end

function assert_is_symbol(token::Token, symbol::String)
     @assert isa(token.value, SymbolToken) "Token '$(token.value)' is not a SymbolToken"
     @assert token.value.symbol == symbol "expecting symbol '$(symbol)' instead of '$(token.value)'"
end

function assert_is_number(token::Token, number::Float64)
     @assert isa(token.value, LiteralNumberToken) "Token '$(token.value)' is not a LiteralNumberToken"
     @assert token.value.number == number "Token '$(token.value)' is not equal to number '$(number)'"
end

function assert_is_string(token::Token, string::String)
     @assert isa(token.value, StringToken) "Token '$(token.value)' is not a StringToken"
     @assert token.value.string == string "Token '$(token.value)' is not equal to string '$(string)'"
end

@testset "test_input_file" begin
     stream = InputStream(IOBuffer("abc   \nd\nef"))

     @test stream.location.line_num == 1
     @test stream.location.col_num == 1

     @test Raytracing.Interpreter.read_char(stream) == "a"
     @test stream.location.line_num == 1
     @test stream.location.col_num == 2

     Raytracing.Interpreter.unread_char(stream, "A")
     @test stream.location.line_num == 1
     @test stream.location.col_num == 1

     @test Raytracing.Interpreter.read_char(stream) == "A"
     @test stream.location.line_num == 1
     @test stream.location.col_num == 2

     @test Raytracing.Interpreter.read_char(stream) == "b"
     @test stream.location.line_num == 1
     @test stream.location.col_num == 3

     @test Raytracing.Interpreter.read_char(stream) == "c"
     @test stream.location.line_num == 1
     @test stream.location.col_num == 4

     skip_whitespaces_and_comments(stream)

     @test Raytracing.Interpreter.read_char(stream) == "d"
     @test stream.location.line_num == 2
     @test stream.location.col_num == 2

     @test Raytracing.Interpreter.read_char(stream) == "\n"
     @test stream.location.line_num == 3
     @test stream.location.col_num == 1

     @test Raytracing.Interpreter.read_char(stream) == "e"
     @test stream.location.line_num == 3
     @test stream.location.col_num == 2

     @test Raytracing.Interpreter.read_char(stream) == "f"
     @test stream.location.line_num == 3
     @test stream.location.col_num == 3

     @test Raytracing.Interpreter.read_char(stream) == ""
end

@testset "test_skip_whitespaces_and_comments" begin
     stream = IOBuffer("""
     # This is a comment
     # This is another comment
     #=
     this is a longer
     comment that i want to avoid
     =#

     #=
     this is anothe long comment where i want to use
     the symbols = and # indifferently, because inside
     this sequence made of # + = and = + # i want to
     ignore anything.
     =#
     FLOAT var(150) # Comment at the end of the line

""")

     scene = parse_scene(InputStream(stream, "test_skip_whitespaces_and_comments"))

     @test length(scene.float_variables) == 1
     @test "var" ∈ keys(scene.float_variables)
     @test scene.float_variables["var"] == 150.0

end

@testset "test_lexer" begin

     stream = IOBuffer("""
     # This is a comment
     # This is another comment
     NEW MATERIAL sky_material(
          DIFFUSE(IMAGE("my file.pfm")),
          <5.0, 500.0, 300.0>
     ) # Comment at the end of the line
""")

     input_file = InputStream(stream, "test_lexer")

     assert_is_keyword(read_token(input_file), Raytracing.Interpreter.NEW)
     assert_is_keyword(read_token(input_file), Raytracing.Interpreter.MATERIAL)
     assert_is_identifier(read_token(input_file), "sky_material")
     assert_is_symbol(read_token(input_file), "(")
     assert_is_keyword(read_token(input_file), Raytracing.Interpreter.DIFFUSE)
     assert_is_symbol(read_token(input_file), "(")
     assert_is_keyword(read_token(input_file), Raytracing.Interpreter.IMAGE)
     assert_is_symbol(read_token(input_file), "(")
     assert_is_string(read_token(input_file), "my file.pfm")
     assert_is_symbol(read_token(input_file), ")")

     @test 1==1
end

@testset "test_parser" begin
     stream = IOBuffer("""
     FLOAT clock(150)

     MATERIAL sky_material(
          DIFFUSE(UNIFORM(<0, 0, 0>)),
          UNIFORM(<0.7, 0.5, 1>)
     )

     # Here is a comment

     MATERIAL ground_material(
          DIFFUSE(CHECKERED(<0.3, 0.5, 0.1>,
                              <0.1, 0.2, 0.5>, 4)),
          UNIFORM(<0, 0, 0>)
     )

     MATERIAL sphere_material(
          SPECULAR(UNIFORM(<0.5, 0.5, 0.5>)),
          UNIFORM(<0, 0, 0>)
     )

     PLANE (sky_material, TRANSLATION([0, 0, 100]) * ROTATION_Y(clock))
     PLANE (ground_material, IDENTITY)

     SPHERE(sphere_material, TRANSLATION([0, 0, 1]))

     CAMERA(PERSPECTIVE, ROTATION_Z(30) * TRANSLATION([-4, 0, 1]), 2.0)
     """)

     scene = parse_scene(InputStream(stream, "test_parser"))

     # Check that the FLOAT variables are ok

     @test length(scene.float_variables) == 1
     @test "clock" ∈ keys(scene.float_variables)
     @test scene.float_variables["clock"] == 150.0

     # Check that the materials are ok

     @test length(scene.materials) == 3
     @test "sphere_material" ∈ keys(scene.materials)
     @test "sky_material" ∈ keys(scene.materials)
     @test "ground_material" ∈ keys(scene.materials)

     sphere_material = scene.materials["sphere_material"]
     sky_material = scene.materials["sky_material"]
     ground_material = scene.materials["ground_material"]

     @test isa(sky_material.brdf, DiffuseBRDF)
     @test isa(sky_material.brdf.pigment, UniformPigment)
     @test sky_material.brdf.pigment.color ≈ RGB(0., 0., 0.)

     @test isa(ground_material.brdf, DiffuseBRDF)
     @test isa(ground_material.brdf.pigment, CheckeredPigment)
     @test ground_material.brdf.pigment.color1 ≈ RGB(0.3, 0.5, 0.1)
     @test ground_material.brdf.pigment.color2 ≈ RGB(0.1, 0.2, 0.5)
     @test ground_material.brdf.pigment.num_steps == 4

     @test isa(sphere_material.brdf, SpecularBRDF)
     @test isa(sphere_material.brdf.pigment, UniformPigment)
     @test sphere_material.brdf.pigment.color ≈ RGB(0.5, 0.5, 0.5)

     @test isa(sky_material.emitted_radiance, UniformPigment)
     @test sky_material.emitted_radiance.color ≈ RGB(0.7, 0.5, 1.0)
     @test isa(ground_material.emitted_radiance, UniformPigment)
     @test ground_material.emitted_radiance.color ≈ RGB(0, 0, 0)
     @test isa(sphere_material.emitted_radiance, UniformPigment)
     @test sphere_material.emitted_radiance.color ≈ RGB(0, 0, 0)

     # Check that the shapes are ok

     @test length(scene.world.shapes) == 3
     @test isa(scene.world.shapes[1], Plane)
     @test scene.world.shapes[1].T ≈ translation(Vec(0, 0, 100)) * rotation_y(150.0)
     @test isa(scene.world.shapes[2], Plane)
     @test scene.world.shapes[2].T ≈ Transformation()
     @test isa(scene.world.shapes[3], Sphere)
     @test scene.world.shapes[3].T ≈ translation(Vec(0, 0, 1))

     # Check that the camera is ok

     @test isa(scene.camera, PerspectiveCamera)
     @test scene.camera.T ≈ rotation_z(30.0) * translation(Vec(-4, 0, 1))
     @test 2.0 ≈ scene.camera.d
end

@testset "test_parser_undefined_material" begin
     # Check that unknown materials raises a GrammarError
     stream = IOBuffer("""
     plane(this_material_does_not_exist, identity)
     """)

     @test_throws GrammarError parse_scene(InputStream(stream))
end

@testset "test_parser_double_camera" begin
     # Check that defining two cameras in the same file raises a GrammarError
     stream = IOBuffer("""
     camera(perspective, rotation_z(30) * translation([-4, 0, 1]), 1.0, 1.0)
     camera(orthogonal, identity, 1.0, 1.0)
     """)

    @test_throws GrammarError parse_scene(InputStream(stream))
end

@testset "test_math_operations_1" begin
     stream = IOBuffer("""
          FLOAT var1(1.0)
          FLOAT var2(2.0)
          FLOAT var3(0.5)
          FLOAT var4(0.5)
          FLOAT var5(-var1 * var2 - 0.5/(var3 + var4))

          ASSERT(var5, -2.5)
     """)

     scene = parse_scene(InputStream(stream, "test_math_operations_1"))
     @test 1==1
end

@testset "test_math_operations_2" begin
     stream = IOBuffer("""
          FLOAT var1(1.0)
          FLOAT var2(2.0)
          FLOAT var3(0.5)
          FLOAT var4(0.5)
          FLOAT var5(-var1 * var2 - 0.5/(var3 + var4))

          ASSERT(var5, 2.5)
     """)

     @test_throws AssertionError parse_scene(InputStream(stream))
end

@testset "test_return_token_value_1" begin
     stream = IOBuffer("""
          VECTOR try( ([1,2,3] + [3,2,1])*2.0 - [1,1,1] )
          ASSERT (try*1.0, -1.0*[-7,-7,-7])
     """)

     scene = parse_scene(InputStream(stream, "test_return_token_value_1"))
     @test 1==1
end

@testset "test_return_token_value_2" begin
     stream = IOBuffer("""
          VECTOR try( ([1,2,3] + [3,2,1])*2.0 - [1,1,1] )
          ASSERT (try, -1.0*[7,-7,-7])
     """)

     @test_throws AssertionError parse_scene(InputStream(stream))
end

@testset "test_parse_VECTOR_1" begin
     stream = IOBuffer("""
          VECTOR v1([1,2,3])
          VECTOR v2([4,5,6] * 1.0 - v1)
          ASSERT(v2-v1*0, 3*[1,1,1] )
     """)

     scene = parse_scene(InputStream(stream, "test_parse_VECTOR_1"))
     @test 1==1
end

@testset "test_parse_VECTOR_2" begin
     stream = IOBuffer("""
          VECTOR v1([1,2,3])
          VECTOR v2([4,5,6] * 1.0 - v1)
          ASSERT(v2-v1*0, 1*[1,1,1] )
     """)

     @test_throws AssertionError parse_scene(InputStream(stream))
end

@testset "test_parse_COLOR_1" begin
     stream = IOBuffer("""
          COLOR c1(<1,2,3>)
          COLOR c2(<4,5,6> * 1.0 - c1)
          ASSERT(c2-c1*0, 3*<1,1,1> )
     """)

     scene = parse_scene(InputStream(stream, "test_parse_COLOR_1"))
     @test 1==1
end

@testset "test_parse_COLOR_2" begin
     stream = IOBuffer("""
          COLOR c1(<1,2,3>)
          COLOR c2(<4,5,6> * 1.0 - c1)
          ASSERT(c2-c1*0, 0*<1,1,1> )
     """)

     @test_throws AssertionError parse_scene(InputStream(stream))
end

@testset "test_ASSERT_1" begin
     stream = IOBuffer("""
          ASSERT(1.0, 1, "=")
          ASSERT(1.0, 1, "==")
          ASSERT(1.0, 2.0, "<")
          ASSERT(1.0, 2.0, "<=")
          ASSERT(3, 2.0, ">")
          ASSERT(3, 2.0, ">=")
     """)

     scene = parse_scene(InputStream(stream, "test_ASSERT_1"))
     @test 1==1
end

@testset "test_ASSERT_2" begin
     @test_throws AssertionError parse_scene(InputStream(IOBuffer("""ASSERT(0.0, 1, "=")""")))
     @test_throws AssertionError parse_scene(InputStream(IOBuffer("""ASSERT(0.0, 1, "==")""")))
     @test_throws AssertionError parse_scene(InputStream(IOBuffer("""ASSERT(1.0, 1, ">")""")))
     @test_throws AssertionError parse_scene(InputStream(IOBuffer("""ASSERT(1.0, 2, ">=")""")))
     @test_throws AssertionError parse_scene(InputStream(IOBuffer("""ASSERT(1.0, 1, "<")""")))
     @test_throws AssertionError parse_scene(InputStream(IOBuffer("""ASSERT(1.0, 0, "<=")""")))
end

@testset "test_tutorial_basic_syntax.txt" begin
     cd("..")
     open("examples/tutorial_basic_syntax.txt") do stream
          scene = parse_scene(InputStream(stream, "tutorial_basic_syntax.txt"))
          @test 1==1
     end

     @test 1==1
     cd("test")
end

@testset "test_demo_world_B.txt" begin
     cd("..")
     open("examples/demo_world_B.txt") do stream
          scene = parse_scene(InputStream(stream, "demo_world_B.txt"))
          @test 1==1
     end

     @test 1==1
     cd("test")
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for pcg.jl

pcg = PCG()

@test pcg.state == UInt64(1753877967969059832)
@test pcg.inc == UInt64(109)

for expected ∈ [2707161783, 2068313097,
                3122475824, 2211639955,
                3215226955, 3421331566]
    @test expected == rand(pcg, UInt32)
end

ext: jl


text:
"""
    gd(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, stepsize::Number=0.1, numepoch::Number=3, scheduling::AbstractString="robbins-monro", g::Number=0.9, epsilon::Number=1.0e-8, lower::Number=0, upper::Number=1.0f+38, expvar::Number=0.1f0, evalfreq::Number=5000, offsetFull::Number=1f-20, offsetStoch::Number=1f-6, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)

Gradient descent method.

Input Arguments
---------
- `input` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `stepsize` : The parameter used in every iteration.
- `numepoch` : The number of epoch.
- `scheduling` : Learning parameter scheduling. `robbins-monro`, `momentum`, `nag`, and `adagrad` are available.
- `g` : The parameter that is used when scheduling is specified as nag.
- `epsilon` : The parameter that is used when scheduling is specified as adagrad.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `evalfreq` : Evaluation Frequency of Reconstruction Error
- `offsetFull` : Off set value for avoding overflow when calculating full gradient
- `offsetStoch` : Off set value for avoding overflow when calculating stochastic gradient
- `initW` : The CSV file saving the initial values of eigenvectors.
- `initV` : The CSV file saving the initial values of loadings.
- `logdir` : The directory where intermediate files are saved, in every 1000 iteration.
- `perm` : Whether the data matrix is shuffled at random.
- `cper` : Count per X (e.g. CPM: Count per million <1e+6>)

Output Arguments
---------
- `W` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `V` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
- stop : Whether the calculation is converged
"""
function gd(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, stepsize::Number=0.1, numepoch::Number=3, scheduling::AbstractString="robbins-monro", g::Number=0.9, epsilon::Number=1.0e-8, lower::Number=0, upper::Number=1.0f+38, expvar::Number=0.1f0, evalfreq::Number=5000, offsetFull::Number=1f-20, offsetStoch::Number=1f-6, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)
    # Initial Setting
    pca = GD()
    if scheduling == "robbins-monro"
        scheduling = ROBBINS_MONRO()
    elseif scheduling == "momentum"
        scheduling = MOMENTUM()
    elseif scheduling == "nag"
        scheduling = NAG()
    elseif scheduling == "adagrad"
        scheduling = ADAGRAD()
    else
        error("Specify the scheduling as robbins-monro, momentum, nag or adagrad")
    end
    pseudocount, stepsize, g, epsilon, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetFull, offsetStoch = init(input, pseudocount, stepsize, g, epsilon, dim, rowmeanlist, rowvarlist, colsumlist, initW, initV, logdir, pca, lower, upper, evalfreq, offsetFull, offsetStoch, cper, scale)
    # Perform PCA
    out = gd(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, stepsize, numepoch, scheduling, g, epsilon, logdir, pca, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetFull, offsetStoch, perm, cper)
    if outdir isa String
        output(outdir, out, expvar)
    end
    return out
end

function gd(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, stepsize, numepoch, scheduling, g, epsilon, logdir, pca, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetFull, offsetStoch, perm, cper)
    # If not 0 the calculation is converged
    stop = 0
    s = 1
    n = 1
    # Each epoch s
    progress = Progress(numepoch)
    while(stop == 0 && s <= numepoch)
        next!(progress)
        # Update Eigen vector
        W, v = gdupdate(scheduling, stepsize, g, epsilon, D, N, M, W, v, s, input, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, offsetFull, offsetStoch, perm, cper)
        # NaN
        checkNaN(W, pca)
        # Retraction
        W .= Array(qr!(W).Q)
        # save log file
        if logdir isa String
            stop = outputlog(s, input, dim, logdir, W, pca, TotalVar, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, lower, upper, stop, cper)
        end
        s += 1
    end

    # Return, W, λ, V
    out = WλV(W, input, dim, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, TotalVar, cper)
    return (out[1], out[2], out[3], out[4], out[5], out[6], stop)
end

# GD × Robbins-Monro
function gdupdate(scheduling::ROBBINS_MONRO, stepsize, g, epsilon, D, N, M, W, v, s, input, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, offsetFull, offsetStoch, perm, cper)
    W .= W .+ ∇f(W, input, D, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, stepsize/s, offsetFull, offsetStoch, perm, cper)
    v = nothing
    return W, v
end

# GD × Momentum
function gdupdate(scheduling::MOMENTUM, stepsize, g, epsilon, D, N, M, W, v, s, input, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, offsetFull, offsetStoch, perm, cper)
    v .= g .* v .+ ∇f(W, input, D, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, stepsize/s, offsetFull, offsetStoch, perm, cper)
    W .= W .+ v
    return W, v
end

# GD × NAG
function gdupdate(scheduling::NAG, stepsize, g, epsilon, D, N, M, W, v, s, input, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, offsetFull, offsetStoch, perm, cper)
    v = g .* v + ∇f(W - g .* v, input, D, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, stepsize/s, offsetFull, offsetStoch, perm, cper)
    W .= W .+ v
    return W, v
end

# GD × Adagrad
function gdupdate(scheduling::ADAGRAD, stepsize, g, epsilon, D, N, M, W, v, s, input, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, offsetFull, offsetStoch, perm, cper)
    grad = ∇f(W, input, D, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, stepsize/s, offsetFull, offsetStoch, perm, cper)
    grad = grad / stepsize
    v .= v .+ grad .* grad
    W .= W .+ stepsize ./ (sqrt.(v) .+ epsilon) .* grad
    return W, v
end


text:
export GaussianCopulaVCObs, GaussianCopulaVCModel
export fitted

"""
GaussianCopulaVCObs
GaussianCopulaVCObs(y, X, V)
A realization of Gaussian copula variance component data.
"""
struct GaussianCopulaVCObs{T <: BlasReal}
    # data
    y::Vector{T}
    X::Matrix{T}
    V::Vector{Matrix{T}}
    n::Int
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # working arrays
    ∇β::Vector{T}   # gradient wrt β
    ∇τ::Vector{T}   # gradient wrt τ
    ∇θ::Vector{T}   # gradient wrt θ
    Hβ::Matrix{T}   # Hessian wrt β
    Hτ::Matrix{T}   # Hessian wrt τ
    res::Vector{T}  # residual vector res_i
    xtx::Matrix{T}  # Xi'Xi
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    storage_n::Vector{T}
    storage_p::Vector{T}
    m1::Vector{T}
    m2::Vector{T}
    Hθ::Matrix{T}   # Hessian wrt variance components θ
end

function GaussianCopulaVCObs(
    y::Vector{T},
    X::Matrix{T},
    V::Vector{Matrix{T}}
    ) where T <: BlasReal
    n, p, m = size(X, 1), size(X, 2), length(V)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇θ  = Vector{T}(undef, m)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    res = Vector{T}(undef, n)
    xtx = transpose(X) * X
    t   = [tr(V[k])/2 for k in 1:m]
    q   = Vector{T}(undef, m)
    storage_n = Vector{T}(undef, n)
    storage_p = Vector{T}(undef, p)
    m1        = Vector{T}(undef, m)
    m2        = Vector{T}(undef, m)
    Hθ  = Matrix{T}(undef, m, m)
    # constructor
    GaussianCopulaVCObs{T}(y, X, V, n, p, m, ∇β, ∇τ, ∇θ, Hβ,
        Hτ, res, xtx, t, q, storage_n, storage_p, m1, m2, Hθ)
end

"""
GaussianCopulaVCModel
GaussianCopulaVCModel(gcs)
Gaussian copula variance component model, which contains a vector of
`GaussianCopulaVCObs` as data, model parameters, and working arrays.
"""
struct GaussianCopulaVCModel{T <: BlasReal} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{GaussianCopulaVCObs{T}}
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    θ::Vector{T}    # m-vector: [θ12, ..., θm2]
    # working arrays
    ∇β::Vector{T}   # gradient from all observations
    ∇τ::Vector{T}
    ∇θ::Vector{T}
    Hβ::Matrix{T}    # Hessian from all observations
    Hτ::Matrix{T}
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    TR::Matrix{T}   # n-by-m matrix with tik = tr(Vi[k]) / 2
    QF::Matrix{T}   # n-by-m matrix with qik = res_i' Vi[k] res_i
    storage_n::Vector{T}
    storage_m::Vector{T}
    storage_θ::Vector{T}
    # hessain with resp to vc
    # asymptotic covariance for inference
    Ainv::Matrix{T}
    Aevec::Matrix{T}
    M::Matrix{T}
    vcov::Matrix{T}
    ψ::Vector{T}
    Hθ::Matrix{T}   # Hessian wrt variance components θ
    penalized::Bool
end

function GaussianCopulaVCModel(gcs::Vector{GaussianCopulaVCObs{T}}; penalized::Bool = false) where T <: BlasReal
    n, p, m = length(gcs), size(gcs[1].X, 2), length(gcs[1].V)
    β   = Vector{T}(undef, p)
    τ   = Vector{T}(undef, 1)
    θ   = Vector{T}(undef, m)
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇θ  = Vector{T}(undef, m)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    XtX = zeros(T, p, p) # sum_i xi'xi
    TR  = Matrix{T}(undef, n, m) # collect trace terms
    ntotal = 0
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        XtX    .+= gcs[i].xtx
        TR[i, :] = gcs[i].t
    end
    QF        = Matrix{T}(undef, n, m)
    storage_n = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, m)
    storage_θ = Vector{T}(undef, m)
    Ainv    = zeros(T, p + m + 1, p + m + 1)
    Aevec   = zeros(T, p + m + 1, p + m + 1)
    M       = zeros(T, p + m + 1, p + m + 1)
    vcov    = zeros(T, p + m + 1, p + m + 1)
    ψ       = Vector{T}(undef, p + m + 1)
    Hθ  = Matrix{T}(undef, m, m)
    GaussianCopulaVCModel{T}(gcs, ntotal, p, m, β, τ, θ,
        ∇β, ∇τ, ∇θ, Hβ, Hτ, XtX, TR, QF,
        storage_n, storage_m, storage_θ, Ainv, Aevec, M, vcov, ψ, Hθ, penalized)
end

function fitted(
    gc::GaussianCopulaVCObs{T},
    β::Vector{T},
    τ::T,
    θ::Vector{T}) where T <: BlasReal
    n, m = length(gc.y), length(gc.V)
    μ̂ = gc.X * β
    Ω = Matrix{T}(undef, n, n)
    @inbounds for k in 1:m
        Ω .+= θ[k] .* gc.V[k]
    end
    θ02 = inv(τ)
    c = inv(1 + dot(θ, gc.t)) # normalizing constant
    V̂ = Matrix{T}(undef, n, n)
    for j in 1:n
        for i in 1:j-1
            V̂[i, j] = c * σ02 * Ω[i, j]
        end
        V̂[j, j] = c * σ02 * (1 + Ω[j, j] + tr(Ω) / 2)
    end
    LinearAlgebra.copytri!(V̂, 'U')
    μ̂, V̂
end

function loglikelihood!(
    gc::GaussianCopulaVCObs{T},
    β::Vector{T},
    τ::T, # inverse of linear regression variance
    θ::Vector{T},
    needgrad::Bool = false,
    needhess::Bool = false;
    penalized::Bool = false
    ) where T <: BlasReal
    # n, p, m = size(gc.X, 1), size(gc.X, 2), length(gc.V)
    needgrad = needgrad || needhess
    if needgrad
        fill!(gc.∇β, 0)
        fill!(gc.∇τ, 0)
        fill!(gc.∇θ, 0)
    end
    needhess && fill!(gc.Hβ, 0)
    # evaluate copula loglikelihood
    sqrtτ = sqrt(abs(τ))
    update_res!(gc, β)
    standardize_res!(gc, sqrtτ)
    rss  = abs2(norm(gc.res)) # RSS of standardized residual
    tsum = dot(abs.(θ), gc.t)
    logl = - log(1 + tsum) - (gc.n * log(2π) -  gc.n * log(abs(τ)) + rss) / 2
    @inbounds for k in 1:gc.m
        mul!(gc.storage_n, gc.V[k], gc.res) # storage_n = V[k] * res
        if needgrad # ∇β stores X'*Γ*res (standardized residual)
            BLAS.gemv!('T', θ[k], gc.X, gc.storage_n, one(T), gc.∇β)
        end
        gc.q[k] = dot(gc.res, gc.storage_n) / 2
    end
    qsum  = dot(θ, gc.q)
    logl += log(1 + qsum)
    # add L2 ridge penalty
    if penalized
        logl -= 0.5 * dot(θ, θ)
    end
    # gradient
    if needgrad
        inv1pq = inv(1 + qsum)
        if needhess
            BLAS.syrk!('L', 'N', -abs2(inv1pq), gc.∇β, one(T), gc.Hβ) # only lower triangular
            gc.Hτ[1, 1] = - abs2(qsum * inv1pq / τ)
            # # hessian of vc vector use with fit_newton_normal.jl
            inv1pt = inv(1 + tsum)
            gc.m1 .= gc.q
            gc.m1 .*= inv1pq
            gc.m2 .= gc.t
            gc.m2 .*= inv1pt
            # hessian for vc
            fill!(gc.Hθ, 0.0)
            BLAS.syr!('U', one(T), gc.m2, gc.Hθ)
            BLAS.syr!('U', -one(T), gc.m1, gc.Hθ)
            copytri!(gc.Hθ, 'U')
        end
        BLAS.gemv!('T', one(T), gc.X, gc.res, -inv1pq, gc.∇β)
        gc.∇β .*= sqrtτ
        gc.∇τ  .= (gc.n - rss + 2qsum * inv1pq) / 2τ
        gc.∇θ  .= inv1pq .* gc.q .- inv(1 + tsum) .* gc.t
        if penalized
            gc.∇θ .-= θ
        end
    end
    # output
    logl
end

function loglikelihood!(
    gcm::GaussianCopulaVCModel{T},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where T <: BlasReal
    logl = zero(T)
    if needgrad
        fill!(gcm.∇β, 0)
        fill!(gcm.∇τ, 0)
        fill!(gcm.∇θ, 0)
    end
    if needhess
        gcm.Hβ .= - gcm.XtX
        gcm.Hτ .= - gcm.ntotal / 2abs2(gcm.τ[1])
    end
    logl = zeros(Threads.nthreads())
    Threads.@threads for i in eachindex(gcm.data)
        @inbounds logl[Threads.threadid()] += loglikelihood!(gcm.data[i], gcm.β,
         gcm.τ[1], gcm.θ, needgrad, needhess; penalized = gcm.penalized)
     end
     @inbounds for i in eachindex(gcm.data)
         if needgrad
             gcm.∇β .+= gcm.data[i].∇β
             gcm.∇τ .+= gcm.data[i].∇τ
             gcm.∇θ .+= gcm.data[i].∇θ
         end
         if needhess
             gcm.Hβ .+= gcm.data[i].Hβ
             gcm.Hτ .+= gcm.data[i].Hτ
             gcm.Hθ .+= gcm.data[i].Hθ
         end
     end
    needhess && (gcm.Hβ .*= gcm.τ[1])
    return sum(logl)
end

# uncomment this and exclude fit_gaussian_vc.jl to fit variance components separately using MM-algorithm instead of Joint Newton.
# function fit!(
#     gcm::GaussianCopulaVCModel,
#     solver=Ipopt.IpoptSolver(print_level=0)
#     )
#     initialize_model!(gcm)
#     optm = MathProgBase.NonlinearModel(solver)
#     lb = fill(-Inf, gcm.p)
#     ub = fill( Inf, gcm.p)
#     MathProgBase.loadproblem!(optm, gcm.p, 0, lb, ub, Float64[], Float64[], :Max, gcm)
#     MathProgBase.setwarmstart!(optm, gcm.β)
#     MathProgBase.optimize!(optm)
#     optstat = MathProgBase.status(optm)
#     optstat == :Optimal || @warn("Optimization unsuccesful; got $optstat")
#     copy_par!(gcm, MathProgBase.getsolution(optm))
#     loglikelihood!(gcm)
#     gcm
# end

# function MathProgBase.initialize(
#     gcm::GaussianCopulaVCModel,
#     requested_features::Vector{Symbol})
#     for feat in requested_features
#         if !(feat in [:Grad, :Hess])
#             error("Unsupported feature $feat")
#         end
#     end
# end

# MathProgBase.features_available(gcm::GaussianCopulaVCModel) = [:Grad, :Hess]

# function MathProgBase.eval_f(
#     gcm::GaussianCopulaVCModel,
#     par::Vector)
#     copy_par!(gcm, par)
#     # maximize σ2 and τ at current β using MM
#     update_Σ!(gcm)
#     # evaluate loglikelihood
#     loglikelihood!(gcm, false, false)
# end

# function MathProgBase.eval_grad_f(
#     gcm::GaussianCopulaVCModel,
#     grad::Vector,
#     par::Vector)
#     copy_par!(gcm, par)
#     # maximize σ2 and τ at current β using MM
#     update_θ!(gcm)
#     # evaluate gradient
#     logl = loglikelihood!(gcm, true, false)
#     copyto!(grad, gcm.∇β)
#     nothing
# end

# function copy_par!(
#     gcm::GaussianCopulaVCModel,
#     par::Vector)
#     copyto!(gcm.β, par)
#     par
# end

# function MathProgBase.hesslag_structure(gcm::GaussianCopulaVCModel)
#     Iidx = Vector{Int}(undef, (gcm.p * (gcm.p + 1)) >> 1)
#     Jidx = similar(Iidx)
#     ct = 1
#     for j in 1:gcm.p
#         for i in j:gcm.p
#             Iidx[ct] = i
#             Jidx[ct] = j
#             ct += 1
#         end
#     end
#     Iidx, Jidx
# end

# function MathProgBase.eval_hesslag(
#     gcm::GaussianCopulaVCModel{T},
#     H::Vector{T},
#     par::Vector{T},
#     σ::T,
#     μ::Vector{T}) where T <: BlasReal
#     copy_par!(gcm, par)
#     # maximize σ2 and τ at current β using MM
#     update_θ!(gcm)
#     # evaluate Hessian
#     loglikelihood!(gcm, true, true)
#     # copy Hessian elements into H
#     ct = 1
#     for j in 1:gcm.p
#         for i in j:gcm.p
#             H[ct] = gcm.Hβ[i, j]
#             ct += 1
#         end
#     end
#     H .*= σ
# end


text:
function cant(gun::Gun)
    h = gun.QE*sind(gun.Φ)
    v =gun.QE*(1-cosd(gun.Φ))
    return h,v
end

function windEffect(target::AbstractTarget, proj::AbstractPenetrator, gun::Gun,wind::Array{Float64,1})
    gun.QE,gun.AZ = QEfinderMPMM!(target, proj, gun)
    u0 = iniCond(gun, proj.calibre)
    proj.position = [u0[1],u0[2],u0[3]]
    proj.velocity = [u0[4],u0[5],u0[6]]
    impactBaseline=trajectoryMPMM(proj, target, gun)[2]
    impactWind=trajectoryMPMM(proj, target, gun, wind)[2]
    return impactBaseline-impactWind

end

"""
    windCrossError(σ,wind,tank, weapon,proj, target, aero;args ...)

Returns the normal distribution for the impact points (horizontal coordinate).
    The function performs a Monte Carlo simulation by sampling the wind (cross field).
        It is assumes that the cross wind variation is a gaussian distribution with
        μ = wind.cross and standard deviation = σ.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function windCrossError(σ::Float64,wind::Wind,tank::Tank, weapon::Gun,proj::AbstractPenetrator, target::AbstractTarget, aero::DataFrame;N=100)
    d = Normal(wind.cross,σ)
#drange = Normal(w.range,σrange)

    x = rand(d, N)
#xrange = rand(drange, 100)

    res = zeros(N)#Array{Float64}#[]
    for i=1:N

        wind.cross = x[i]
#w.range = xrange

        windIn = Cuas.wind(tank, wind)
        res[i]=trajectoryMPMM(proj, target, weapon,aero,w_bar=windIn)[1][3]
#push!(res,tmp)
    end

    σcrossₜ=fit(Normal, res)
    return σcrossₜ
end

"""
    windRangeError(σ,wind,tank, weapon,proj, target, aero;args ...)

Returns the normal distribution for the impact points (vertical coordinate).
    The function performs a Monte Carlo simulation by sampling the wind (range field).
        It is assumes that the range wind variation is a gaussian distribution with
        μ = wind.range and standard deviation = σ.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function windRangeError(σ::Float64,wind::Wind,tank::Tank, weapon::Gun,proj::AbstractPenetrator, target::AbstractTarget, aero::DataFrame;N=100)

    d = Normal(wind.range,σ)
    x = rand(d, N)
    res = zeros(N)#Array{Float64}#[]
    for i=1:N

    #w.cross = xcross[i]
    wind.range = x[i]

    windIn = Cuas.wind(tank, wind)
    res[i]=trajectoryMPMM(proj, target, weapon,aero,w_bar=windIn)[1][2]
    #push!(res,tmp)
    end

    σrangeₜ =fit(Normal, res)

    return σrangeₜ
end

"""
    interTable(range,jumpTable, value)

Returns the 'value' from an interpolation table 'jumpTable' over the 'range'

"""
function interTable(range::Array{Float64,1},jumpTable::Array{Float64,1}, value::Float64)
    inter=LinearInterpolation(range, jumpTable, extrapolation_bc=Flat())
    inter(value)
end

"""
    rangeError(target, σ, tank, weapon, proj, aero; args)

Returns the normal distribution for elevation and azimuth angles.
    The function performs a Monte Carlo simulation by sampling the target range.
        It is assumes that the range variation is a gaussian distribution with
        μ = target.ρ and standard deviation = σ. For each sampled range, a elevation and azimuth angle is computed.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function rangeError(target::AbstractTarget, σ::Float64, tank::Tank, weapon::Gun, proj::AbstractPenetrator, aero::DataFrame; N=100)

    d = Normal(target.ρ,σ)
    x = rand(d, N)
    QE = zeros(N)#Array{Float64}#[]
    AZ = zeros(N)
    for i=1:N
        target.ρ = x[i]
        target.position = targetPos(target, tank)

    QE[i],AZ[i]=QEfinderMPMM!(target, proj, weapon,aero)
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σQE =fit(Normal, QE)
    σAZ = fit(Normal, AZ)

    return σQE,σAZ
end

"""
    elevationError(σ, target,proj,weapon,tank,aero;args )

Returns the normal distribution for the impact points (vertical and horizontal distribution).
    The function performs a Monte Carlo simulation by sampling the elevation angle.
        It is assumes that the elevation angle variation is a gaussian distribution with
        μ = the compute elevation angle and standard deviation = σ.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function elevationError(σ::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame;N=100 )
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero)
    weapon.AZ = rad2deg(AZ)
    d = Normal(QE,σ)
    x = rand(d, N)
    resV = zeros(N)#Arra
    resH = zeros(N)#Arra

    for i=1:N
        weapon.QE=rad2deg(x[i])
        target.position= targetPos(target, tank)
        proj.position=muzzlePos(tank)
        proj.velocity=muzzleVel(tank)

    resV[i], resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][2:3]
    #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σθv =fit(Normal, resV)
    σθh =fit(Normal, resH)

    return σθv ,σθh

end

"""
    azimuthError(σ, target,proj,weapon,tank,aero;args )

Returns the normal distribution for the impact points (vertical and horizontal distribution).
    The function performs a Monte Carlo simulation by sampling the azimuth angle.
        It is assumes that the elevation angle variation is a gaussian distribution with
        μ = the computed azimuth angle and standard deviation = σ.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function azimuthError(σ::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame;N=100 )
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero)
    weapon.QE = rad2deg(QE)
    d = Normal(AZ,σ)
    x = rand(d, N)
    resV = zeros(N)#Arra
    resH = zeros(N)#Arra

    for i=1:N
        weapon.AZ=rad2deg(x[i])
        target.position= targetPos(target, tank)
        proj.position=muzzlePos(tank)
        proj.velocity=muzzleVel(tank)

        resV[i], resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][2:3]
        #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σξv =fit(Normal, resV)
    σξh =fit(Normal, resH)

    return σξv ,σξh

end

"""
    muzzleVelError(σ, target,proj,weapon,tank,aero;args )

Returns the normal distribution for the impact points (vertical and horizontal distribution).
    The function performs a Monte Carlo simulation by sampling the muzzle velocity.
        It is assumes that the muzzle velocity variation is a gaussian distribution with
        μ = weapon.u₀ and standard deviation = σ.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function muzzleVelError(σ::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame;N=100)
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero)
    weapon.QE = rad2deg(QE)
    weapon.AZ = rad2deg(AZ)
    target.position= targetPos(target, tank)
    proj.position=muzzlePos(tank)
    d = Normal(weapon.u₀,σ)
    x = rand(d, N)
    resV = zeros(N)#Arra
    resH = zeros(N)#Arra

    for i=1:N

        weapon.u₀ = x[i]
        proj.velocity=muzzleVel(tank)

        resV[i], resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][2:3]
        #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σv =fit(Normal, resV)
    σh =fit(Normal, resH)

    return σv ,σh
end

"""
    temperatureError(σ, target,proj,weapon,tank,aero, atmosphere;args )

Returns the normal distribution for the impact points (vertical and horizontal distribution).
    The function performs a Monte Carlo simulation by sampling the air temperature.
        It is assumes that the temperature variation is a gaussian distribution with
        μ = atmosphere.t and standard deviation = σ.
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        * atmosphere: the atmospheric characteristics
        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function temperatureError(σ::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame,atmosphere::Air;N=100)
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero,atmosphere=atmosphere)
    weapon.QE = rad2deg(QE)
    weapon.AZ = rad2deg(AZ)
    proj.velocity=muzzleVel(tank)
    target.position= targetPos(target, tank)
    proj.position=muzzlePos(tank)
    d = Normal(atmosphere.t,σ)
    x = rand(d, N)
    resV = zeros(N)#Arra
    resH = zeros(N)#Arra

    for i=1:N

        atmosphere.t = x[i]

        resV[i], resH[i]=trajectoryMPMM(proj, target, weapon,aero,atm=atmosphere)[1][2:3]
        #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σv =fit(Normal, resV)
    σh =fit(Normal, resH)

    return σv ,σh

end

function pressureError(σ::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame,atmosphere::Air;N=100)
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero,atmosphere=atmosphere)
    weapon.QE = rad2deg(QE)
    weapon.AZ = rad2deg(AZ)
    proj.velocity=muzzleVel(tank)
    target.position= targetPos(target, tank)
    proj.position=muzzlePos(tank)
    d = Normal(atmosphere.p,σ)
    x = rand(d, N)
    resV = zeros(N)#Arra
    resH = zeros(N)#Arra

    for i=1:N

        atmosphere.p = x[i]

        resV[i], resH[i]=trajectoryMPMM(proj, target, weapon,aero,atm=atmosphere)[1][2:3]
        #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σv =fit(Normal, resV)
    σh =fit(Normal, resH)

    return σv ,σh

end

function gunPositionError(σ::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame,atmosphere::Air;N=100)
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero,atmosphere=atmosphere)
    weapon.QE = rad2deg(QE)
    weapon.AZ = rad2deg(AZ)
    proj.velocity=muzzleVel(tank)
    target.position= targetPos(target, tank)
    proj.position=muzzlePos(tank)
    dx = Normal(proj.position[1],σ)
    dy = Normal(proj.position[2],σ)
    dz = Normal(proj.position[3],σ)
    x = rand(dx, N)
    y = rand(dy, N)
    z = rand(dz, N)
    resV = zeros(N)#Arra
    resH = zeros(N)#Arra

    for i=1:N

        proj.position = [x[i],y[i],z[i]]

        resV[i], resH[i]=trajectoryMPMM(proj, target, weapon,aero,atm=atmosphere)[1][2:3]
        #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σv =fit(Normal, resV)
    σh =fit(Normal, resH)

    return σv ,σh

end

"""
    cantError(σ, σmeas, target,proj,weapon,tank,aero, atmosphere;args )

Returns the normal distribution for the impact points (vertical and horizontal distribution).
    The function performs a Monte Carlo simulation by sampling the cant angle.
        It is assumes that the cant variation is a gaussian distribution with
        μ = tank.hull.Φ and standard deviation = σ. On this value of cant we perform another sampling to
        take into account the measurement error. We assume also a normal distribution with μ = sampled cant and
        the standard deviation = σmeas
        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile

        One optional argument is N (the number of Monte Carlo simulations). Default = 100

"""
function cantError(σ::Float64, σmeas::Float64, target::AbstractTarget,proj::AbstractPenetrator,weapon::Gun,tank::Tank,aero::DataFrame;N=100)
    QE,AZ=QEfinderMPMM!(target, proj, weapon,aero)
    weapon.QE = rad2deg(QE)
    weapon.AZ = rad2deg(AZ)
    tank.canon.θ = rad2deg(QE)
    tank.turret.ξ = rad2deg(AZ)

    target.position = targetPos(target, tank)
    proj.position=muzzlePos(tank)
    proj.velocity=muzzleVel(tank)

    d1 = Normal(tank.hull.Φ,σ)

    x1 = rand(d1, N)

    resV1 = zeros(N)#Arra
    resH1 = zeros(N)#Arra

    for i=1:N

        #tank.hull.Φ = x[i]

        d2 = Normal(x1[i],σmeas)
        resV2 = zeros(N)#Arra
        resH2 = zeros(N)#Arra

        x2 = rand(d2, N)
        for j=1:N

            tank.hull.Φ = x2[j]

            target.position = targetPos(target, tank)
            proj.position=muzzlePos(tank)
            proj.velocity=muzzleVel(tank)

            resV2[j], resH2[j]=trajectoryMPMM(proj, target, weapon,aero)[1][2:3]
        end
        resV1[i] =fit(Normal, resV2).μ
        resH1[i] =fit(Normal, resH2).μ
        #resH[i]=trajectoryMPMM(proj, target, weapon,aero)[1][3]
    #println("QE", " ", QE[i])
    #println("AZ", " ", AZ[i])
    end
    σv =fit(Normal, resV1)
    σh =fit(Normal, resH1)

    return σv ,σh

end

"""
variableBias(target,tank, weapon, proj, aero, wind, atmosphere;args ...)

Returns the variable bias error as an objet "Error".
optional arguments are all the standard deviation of error sources: σrange, σcross_wind,σrange_wind,
σjump,σfc,σzeroing,σboresight,σMv, σtemp, σpress, σopt, σgun, σcant, σcant_meas
default values for optional parameters are nothing or nul.

        * tank
        * weapon : the launcher of type Gun
        * proj : the projectile
        * target
        * aero : the aerodynamic coefficients of the projectile
        * wind
        * atmosphere

"""
function variableBias(target::AbstractTarget,tank::Tank, weapon::Gun, proj::AbstractPenetrator, aero::DataFrame, wind::Wind, atmosphere::Air;
    σrange=nothing, σcross_wind=nothing,σrange_wind=nothing,σjump=Error(0.0,0.0),σfc=Error(0.0,0.0),σzeroing=Error(0.0,0.0),σboresight=Error(0.0,0.0),
    σMv=nothing, σtemp=nothing, σpress=nothing, σopt=Error(0.0,0.0), σgun=nothing, σcant =nothing, σcant_meas=nothing)

σvariable =Error(0.0,0.0)
σR=Error(0.0,0.0)
σcross_value = 0.0
σrange_value = 0.0
σMVv_value= 0.0
σMVh_value= 0.0
σtempV_value = 0.0
σtempH_value = 0.0
σpressV_value =0.0
σpressH_value =0.0
σgunV_value = 0.0
σgunH_value = 0.0
σcantV_value = 0.0
σcantH_value = 0.0

if isnothing(σrange)==false

dQE,dAZ=rangeError(target,σrange,tank, weapon, proj, aero)

σvEl, σhEl = elevationError(dQE.σ, target,proj,weapon,tank,aero)
σvAz, σhAz = azimuthError(dAZ.σ, target,proj,weapon,tank,aero)

#σR.horizontal = sqrt((rand(σhEl))^2+(rand(σhAz))^2)
#σR.vertical = sqrt((rand(σvEl))^2+(rand(σvAz))^2)

σR.horizontal = sqrt((σhEl.σ)^2+(σhAz.σ)^2)
σR.vertical = sqrt((σvEl.σ)^2+(σvAz.σ)^2)

end

if isnothing(σcross_wind)==false
σcross_dist= windCrossError(σcross_wind,wind,tank,weapon, proj, target, aero) #m
#σcross_value =rand(σcross_dist)
σcross_value =σcross_dist.σ
end

if isnothing(σrange_wind)==false
σrange_dist=windRangeError(σrange_wind,wind,tank,weapon, proj, target, aero) #m
#σrange_value =rand(σrange_dist)
σrange_value =σrange_dist.σ
end

if isnothing(σMv)==false
σMVvₜ, σMVhₜ = muzzleVelError(σMv, target,proj,weapon,tank,aero)

#σMVv_value = rand(σMVvₜ)
#σMVh_value = rand(σMVhₜ)
σMVv_value = σMVvₜ.σ
σMVh_value = σMVhₜ.σ
end

if isnothing(σtemp)==false
σtempV, σtempH = temperatureError(σtemp, target,proj,weapon,tank,aero,atmosphere)

#σtempV_value = rand(σtempV)
#σtempH_value = rand(σtempH)

σtempV_value = σtempV.σ
σtempH_value = σtempH.σ

end

if isnothing(σpress)==false
    σpressV, σpressH = pressureError(σpress, target,proj,weapon,tank,aero,atmosphere)

#    σpressV_value = rand(σpressV)
#    σpressH_value = rand(σpressH)

    σpressV_value = σpressV.σ
    σpressH_value = σpressH.σ
end

if isnothing(σgun)==false
    σgunV, σgunH = gunPositionError(σgun, target,proj,weapon,tank,aero,atmosphere)

#    σgunV_value = rand(σgunV)
#    σgunH_value = rand(σgunH)

    σgunV_value = σgunV.σ
    σgunH_value = σgunH.σ

end

if isnothing(σcant)==false
    σcantV, σcantH = cantError(σcant,σcant_meas, target,proj,weapon,tank,aero)

#    σcantV_value = rand(σcantV)
#    σcantH_value = rand(σcantH)

    σcantV_value = σcantV.σ
    σcantH_value = σcantH.σ
end

σvariable.horizontal = sqrt(σR.horizontal^2+σcross_value^2+σjump.horizontal^2+σfc.horizontal^2+σzeroing.horizontal^2+
σboresight.horizontal^2+σMVh_value^2+σtempH_value^2+σpressH_value^2+σopt.horizontal^2+σgunH_value^2+σcantH_value^2)

σvariable.vertical = sqrt(σR.vertical^2+σrange_value^2+σjump.vertical^2+σfc.vertical^2+σzeroing.vertical^2+σboresight.vertical^2+
σMVv_value^2+σtempV_value^2+σpressV_value^2+σopt.vertical^2+σgunV_value^2+σcantV_value^2)

return σvariable

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Rafael de Pelegrini Soares
#* $Id$
#*--------------------------------------------------------------------
type Supply
	Supply()=begin
		new(
			voltage(),
			voltage(Dict{Symbol,Any}(
				:Default=>0
			)),
			wire(Dict{Symbol,Any}(
				:Brief=>"Inlet",
				:PosX=>0.3923,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			wire(Dict{Symbol,Any}(
				:Brief=>"Outlet",
				:PosX=>0.3984,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			[
				:(outlet.V = V0),
				:(inlet.V - outlet.V = V),
			],
			[
				"","",
			],
			[:V,:V0,],
			[:inlet,:outlet,]
		)
	end
	V::voltage
	V0::voltage
	inlet::wire
	outlet::wire
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Supply
function setEquationFlow(in::Supply)
	addEquation(1)
	addEquation(2)
end
function atributes(in::Supply,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Supply"
	fields[:Brief]="Electrical Supply."
	drive!(fields,_)
	return fields
end
Supply(_::Dict{Symbol,Any})=begin
	newModel=Supply()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
## perofrmance comparison for different TimeArray type definitions
using TimeSeries
using MarketData, MarketTechnicals
using Reactive

# https://github.com/multidis/FinancialSeries.jl/blob/master/src/financialtimeseries.jl
using FinancialSeries

import FinancialSeries.percentchange
function percentchange{T,N, M}(ta::TimeArray{T,N, M}; method="simple")
    logreturn = log(ta.values)[2:end] .- log(lag(ta).values)
#    logreturn = T[ta.values[t] for t in 1:length(ta)] |> log |> diff

    if method == "simple"
      TimeArray(ta.timestamp[2:end], expm1(logreturn), ta.colnames)
    elseif method == "log"
      TimeArray(ta.timestamp[2:end], logreturn, ta.colnames)
    else msg("only simple and log methods supported")
    end
end

aapl = TimeArray(AAPL.timestamp, AAPL.values, AAPL.colnames, FinancialSeries.Stock(FinancialSeries.Ticker("AAPL")))

# log-returns
println("Log returns timing")
for i = 1:10
    @time ret = percentchange(aapl["Adj. Close"], method="log")
end

# reactive signal updates
tsroll = Input(aapl[1:50])
#tsig = lift(ts -> values(sma(ts["Close"], 5))[end] > values(sma(ts["Close"], 40))[end], Bool, tsroll)
tsig = lift(ts -> (sma(ts["Close"], 5) .> sma(ts["Close"], 40)).values[end], Bool, tsroll)

nsteps = 1_000
vsig = Array(Bool, nsteps);
function vsigupd!(vsig::Vector{Bool}, nsteps::Integer)
    for i in 1:nsteps
        push!(tsroll, aapl[i:i+50])
        vsig[i] = tsig.value
    end
end

println(" ")
println("reactive signal timing")
expr_run = parse("@time vsigupd!(vsig, nsteps)")
for i = 1:10
    eval(expr_run)
end


text:
"""
    exponential_solver(rho_u0,rho_v0,LK,LS,LD,t_span,N=0;initial_type = "density")
# Purpose:
- use an exponential integrator to integrate

    ``\\rho_t=-i(L_K+iL_S) \\rho+L_D \rho``

# Argument:
- ``\\rho_{u0},\\rho_{v0}``: ``\\rho0 = ``\\rho_{u0} - i \\rho_{v0}``
- LK: real part of the Hamiltonian operator
- LS: imaginary part of the Hamiltonain operator
- LD: Lindblad operator
- t_span: where the solution will be stored
- initial_type: "density" vectorized density matrix, "states" states vector
# Output:
- rho_u,rho_v: solutions at t_span given by ``\\rho(t_i) = \rho_u(:,i)-i \rho_v(:,i)``
"""
function exponential_solver(rho_u0,rho_v0,
                            LK,LS,LD,
                            t_span,N=0;
                            initial_type = "density")
    if(initial_type == "states")
        rho_u_initial = (rho_u0*transpose(rho_u0)+rho_v0*transpose(rho_v0))[:]
        rho_v_initial = (rho_v0*transpose(rho_u0)-rho_u0*transpose(rho_v0))[:]
        rho_vec0 = [rho_u_initial;rho_v_initial]
        if(N==0)
            N = length(rho_u_initial)
        end
    elseif(initial_type == "density")
        rho_vec0 = [rho_u0;rho_v0];
        if(N==0)
            N = length(rho_u0)
        end
    else
        println("Error! initial_type must be \"density\" or \"states\"")
        return
    end
    #rho_vec = zeros(Float64,length(rho_vec0),length(t_span))
    LL = [LS+LD -LK;
          LK     LS+LD]
    # Zygote does not support functions with setindex!(), so code is organized in this way in order to
    # use auto differentiation
    rho_vec = exp(LL*t_span[1])*rho_vec0
    for time_step = 2:length(t_span)
        rho_vec = [rho_vec exp(LL*t_span[time_step])*rho_vec0]
    end
    return rho_vec[1:N,:],rho_vec[N+1:2*N,:]
end

"""
    exponential_solver_complex(rho_vec0,L,t_span::Array{Float64};initial_type = "density"):

# Argument:
- L: the whole propagation operator
- t_span: where the function will be evaluated stored in
- rho0_vec: initial density/state
- initial_type: decide we aer given an initial density matrix or a state vector

# Output:
- solutions at t_span
"""
function exponential_solver_complex(rho_vec0,L,t_span::Array{Float64};initial_type = "density")
    if(initial_type == "states")
        rho_vec0 = (rho_vec0*adjoint(rho_vec0))[:]
    elseif(initial_type != "density")
        println("Error! initial_type must be \"density\" or \"states\"")
        return
    end
    rho_vec = zeros(ComplexF64,length(rho_vec0),length(t_span))
    for time_step = 1:length(t_span)
        rho_vec[:,time_step] = exp(L*t_span[time_step])*rho_vec0
    end
    return rho_vec
end
"""
    exponential_solver_complex(rho_vec0,L,t_span::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}};
                       initial_type = "density")

# Argument:
- L: the whole propagation operator
- t_span: where the function will be evaluated stored in
- rho0_vec: initial density/state
- initial_type: decide we aer given an initial density matrix or a state vector

# Output:
- Complex-valued solutions at t_span
"""
function exponential_solver_complex(rho_vec0,L,
                            t_span::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}};
                            initial_type = "density")
    if(initial_type == "states")
        rho_vec0 = (rho_vec0*adjoint(rho_vec0))[:]
    elseif(initial_type != "density")
        println("Error! initial_type must be \"density\" or \"states\"")
        return
    end
    rho_vec = zeros(ComplexF64,length(rho_vec0),length(t_span))
    dt = t_span[2]-t_span[1]
    Propagator = exp(L*dt)
    rho_vec[:,1] = exp(L*t_span[1])*rho_vec0
    for time_step = 2:length(t_span)
        rho_vec[:,time_step] = Propagator*rho_vec[:,time_step-1]
    end
    return rho_vec
end

"""
    LindbladODEProblemComplex(rho0,L::Array{ComplexF64,2},time_final::Float64;initial_type = "density")

Function provide interfaces to DifferentialEquations
package to solve: ``\\rho_t = L \\rho``

# Argument:
- L: Lindblad operator
- time: final time or time beining evaluated
- rho0: initial condition
- initial_type: specify initial value is a density matrix/a state vector

# Output:
A complex-valued problem object which we will feed to DifferentialEquations.jl
"""
function LindbladODEProblemComplex(rho0,L::Array{ComplexF64,2},time_final::Float64;initial_type = "density")
    if(initial_type=="states")
        rho0 = (rho0*adjoint(rho0))[:];
    end
    rho0 = convert(Vector{ComplexF64},rho0)
    function ode_problem!(_drho,_rho,_p,_t)
    	_drho .= _p*_rho
    end
    return ODEProblem(ode_problem!,rho0,(0.0,time_final),L)
end

"""
    LindbladODEProblem(rho_u0,rho_v0,LK::Array{Float64,2},LS::Array{Float64},LD::Array{Float64,2},time_final::Float64;initial_type = "density")

Function provides interfaces to DifferentialEquations
package to solve the real-valued Lindblad system:
``(\\rho_u - i \\rho_v)_t = -i (L_K+iL_S)(\\rho_u-i \\rho_v) + L_D(\\rho_u-i \\rho_v)``

# Argument:
- ``\\rho_{u0}, \\rho_{v0}``: initial condition ``\\rho_0 = \\rho_{u0} - i \\rho_{v0}``
- LK: real part of the Hamiltonian operator
- LS: imaginary part of the Hamiltonain operator
- LD: Lindblad operator
- time: final time or time beining evaluated
- initial_type: specify initial value is a density matrix/a state vector

# Output:
A problem object which we will feed to DifferentialEquations.jl
"""
function LindbladODEProblem(rho_u0,rho_v0,LK::Array{Float64,2},LS::Array{Float64},LD::Array{Float64,2},time_final::Float64;initial_type = "density")
    if(initial_type=="states")
        rho_u0,rho_v0 = convert_state_to_density(rho_u0,rho_v0)
    end
    L_real = [LS+LD -LK;
              LK     LS+LD]
    function ode_problem!(_drho,_rho,_p,_t)
    	_drho .= _p*_rho
    end
    return ODEProblem(ode_problem!,[rho_u0;rho_v0],(0.0,time_final),L_real)
end


text:
struct SDFFile{P,B}
    name::String
    header::Header
    blocks::B
    param::Ref{P}
end

include("read.jl")
include("utils.jl")
include("angular_momentum.jl")

Base.keys(sdf::SDFFile) = keys(sdf.blocks)
Base.getindex(sdf::SDFFile, idx::Vararg{AbstractString, N}) where N = sdf[Symbol.(idx)...]


text:
#=
 Methods and types associated with single and double layer potentials
=#

module Layers

using ..Fields
using ..Bodies

export DoubleLayer, SingleLayer, MaskType, Mask, ComplementaryMask

abstract type LayerType{N,NX,NY} end

struct DoubleLayer{N,NX,NY,G,T,DT} <: LayerType{N,NX,NY}
    nds :: VectorData{N,Float64,DT}
    H :: RegularizationMatrix{Edges{G,NX,NY,T},VectorData{N,T,DT}}
end

function DoubleLayer(body::Union{Body,BodyList},H::RegularizationMatrix;weight::Float64 = 1.0)
  ds = ScalarData(dlengthmid(body))*weight
  normals = VectorData(normalmid(body))
  return DoubleLayer(normals∘ds,H)
end

function DoubleLayer(body::Union{Body,BodyList},reg::Regularize{N},w::GridData{NX,NY,T}) where {N,NX,NY,T}

  out = RegularizationMatrix(reg,VectorData(N,dtype=T),Edges(celltype(w),w,dtype=T))
  return DoubleLayer(body,reg._issymmetric ? out[1] : out, weight = sqrt(reg.overdv))
end

(μ::DoubleLayer{N})(p::ScalarData{N}) where {N} = divergence(μ.H*(p∘μ.nds))

function (μ::DoubleLayer{N,NX,NY,G,T,DT})(p::Number) where {N,NX,NY,G,T,DT}
  ϕ = ScalarData(N,dtype=T)
  ϕ .= p
  return μ(ϕ)
end

function Base.show(io::IO, H::DoubleLayer{N,NX,NY,G,T,DT}) where {N,NX,NY,G,T,DT}
    println(io, "Double-layer potential mapping")
    println(io, "  from $N scalar-valued point data of $T type")
    println(io, "  to a $NX x $NY grid of $G nodal data")
end

struct SingleLayer{N,NX,NY,G,T,DT} <: LayerType{N,NX,NY}
    ds :: ScalarData{N,Float64,DT}
    H :: RegularizationMatrix{Nodes{G,NX,NY,T},ScalarData{N,T,DT}}
end

function SingleLayer(body::Union{Body,BodyList},H::RegularizationMatrix)
  ds = ScalarData(dlengthmid(body))
  return SingleLayer(ds,H)
end

function SingleLayer(body::Union{Body,BodyList},reg::Regularize{N},w::GridData{NX,NY,T}) where {N,NX,NY,T}

  out = RegularizationMatrix(reg,ScalarData(N,dtype=T),Nodes(celltype(w),w,dtype=T))
  return SingleLayer(body,reg._issymmetric ? out[1] : out)
end

(μ::SingleLayer{N})(p::ScalarData{N}) where {N} = μ.H*(p∘μ.ds)

function (μ::SingleLayer{N,NX,NY,G,T,DT})(p::Number) where {N,NX,NY,G,T,DT}
  ϕ = ScalarData(N,dtype=T)
  ϕ .= p
  return μ(ϕ)
end

function Base.show(io::IO, H::SingleLayer{N,NX,NY,G,T,DT}) where {N,NX,NY,G,T,DT}
    println(io, "Single-layer potential mapping")
    println(io, "  from $N scalar-valued point data of $T type")
    println(io, "  to a $NX x $NY grid of $G nodal data")
end

abstract type MaskType end

struct Mask{N,NX,NY,G} <: MaskType
  data :: Nodes{G,NX,NY}
end

struct ComplementaryMask{N,NX,NY,G} <: MaskType
  mask :: Mask{N,NX,NY,G}
end

"""
    Mask(b::Body,reg::Regularize,w::GridData)

Returns a `Nodes` mask that sets every grid point equal to 1 if it lies inside the body
`b` and 0 outside. Returns grid data of the same type as `w`. Uses regularization
defined by `reg`.
"""
function Mask(dlayer::DoubleLayer{N,NX,NY,G,T,DT}) where {N,NX,NY,G,T,DT}
  L = plan_laplacian(NX,NY,with_inverse=true,dtype=T)
  return Mask{N,NX,NY,G}(L\dlayer(1))
end

Mask(body::Union{Body,BodyList},regop::Regularize{N},w::GridData{NX,NY,T}) where {N,NX,NY,T} =
    Mask(DoubleLayer(body,regop,w))

(m::Mask{N,NX,NY,G})(w::Nodes{G,NX,NY,T}) where {N,NX,NY,G,T} = m.data ∘ w

(m::ComplementaryMask{N,NX,NY,G})(w::Nodes{G,NX,NY,T}) where {N,NX,NY,G,T} = w - m.mask(w)

end


text:
using Test

@testset "Testing Lab 3" begin
   include("test1.jl")
   include("test2.jl")
   #include("test3.jl")
end;


text:
function rulebased_split_sentences(sentences)
    sentences = replace(sentences, r"([?!.])\s" => Base.SubstitutionString("\\1\n"))

    sentences = postproc_splits(sentences)
    split(sentences, "\n"; keepempty=false)
end

function replace_til_no_change(input, pattern, replacement)
    while(occursin(pattern, input))
        input = replace(input, pattern => replacement)
    end
    input
end

"""
     postproc_splits(sentences)

Applies heuristic rules to repair sentence splitting errors.
Developed for use as postprocessing for the GENIA sentence
splitter on PubMed abstracts, with minor tweaks for
full-text documents.

`sentences` should be a string, with line breaks on sentence boundaries.
Returns a similar string, but more correct.

Based on
https://github.com/ninjin/geniass/blob/master/geniass-postproc.pl
Which is
(c) 2010 Sampo Pyysalo. No rights reserved, i.e. do whatever you like with this.
Which draws in part on heuristics included in Yoshimasa Tsuruoka's
medss.pl script.

"""
function postproc_splits(sentences::AbstractString)
    # Before we do anything remove windows line-ends
    sentences = replace(sentences, "\r" => "")

    # breaks sometimes missing after "?", "safe" cases
    sentences = replace(sentences, r"\b([a-z]+\?) ([A-Z][a-z]+)\b" => Base.SubstitutionString("\\1\n\\2"))
    # breaks sometimes missing after "." separated with extra space, "safe" cases
    sentences = replace(sentences, r"\b([a-z]+ \.) ([A-Z][a-z]+)\b" => Base.SubstitutionString("\\1\n\\2"))

    # no breaks producing lines only containing sentence-ending punctuation
    sentences = replace(sentences, r"\n([.!?]+)\n" => Base.SubstitutionString("\\1\n"))

    # no breaks inside parens/brackets. (To protect against cases where a
    # pair of locally mismatched parentheses in different parts of a large
    # document happens to match, limit size of intervening context. As this
    # is not an issue in cases where there are no intervening brackets,
    # allow an unlimited length match in those cases.)

    # unlimited length for no intevening parens/brackets
    sentences = replace_til_no_change(sentences, r"\[([^\[\]\(\)]*)\n([^\[\]\(\)]*)\]", s"[\1 \2]")
    sentences = replace_til_no_change(sentences, r"\(([^\[\]\(\)]*)\n([^\[\]\(\)]*)\)", s"(\1 \2)")
    # standard mismatched with possible intervening
    sentences = replace_til_no_change(sentences, r"\[([^\[\]]{0,250})\n([^\[\]]{0,250})\]", s"[\1 \2]")
    sentences = replace_til_no_change(sentences, r"\(([^\(\)]{0,250})\n([^\(\)]{0,250})\)", s"(\1 \2)")
    # ... nesting to depth one
    sentences = replace_til_no_change(sentences, r"\[((?:[^\[\]]|\[[^\[\]]*\]){0,250})\n((?:[^\[\]]|\[[^\[\]]*\]){0,250})\]", s"[\1 \2]")
    sentences = replace_til_no_change(sentences, r"\(((?:[^\(\)]|\([^\(\)]*\)){0,250})\n((?:[^\(\)]|\([^\(\)]*\)){0,250})\)", s"(\1 \2)")

    # no break after periods followed by a non-uppercase "normal word"
    # (i.e. token with only lowercase alpha and dashes, with a minimum
    # length of initial lowercase alpha).
    sentences = replace(sentences, r"\.\n([a-z]{3}[a-z-]{0,}[ \.\:\,])" => s". \1")

    # No break after an single letter other than I, which could be an initial in a name
    sentences = replace(sentences, r"(\b[A-HJ-Z]\.)\n" => s"\1 ")

    # no break before CC ...
    sentences = replace(sentences, r"\n(and )" => s" \1")
    sentences = replace(sentences, r"\n(or )" => s" \1")
    sentences = replace(sentences, r"\n(but )" => s" \1")
    sentences = replace(sentences, r"\n(nor )" => s" \1")
    sentences = replace(sentences, r"\n(yet )" => s" \1")
    # or IN. (this is nothing like a "complete" list...)
    sentences = replace(sentences, r"\n(of )" => s" \1")
    sentences = replace(sentences, r"\n(in )" => s" \1")
    sentences = replace(sentences, r"\n(by )" => s" \1")
    sentences = replace(sentences, r"\n(as )" => s" \1")
    sentences = replace(sentences, r"\n(on )" => s" \1")
    sentences = replace(sentences, r"\n(at )" => s" \1")
    sentences = replace(sentences, r"\n(to )" => s" \1")
    sentences = replace(sentences, r"\n(via )" => s" \1")
    sentences = replace(sentences, r"\n(for )" => s" \1")
    sentences = replace(sentences, r"\n(with )" => s" \1")
    sentences = replace(sentences, r"\n(that )" => s" \1")
    sentences = replace(sentences, r"\n(than )" => s" \1")
    sentences = replace(sentences, r"\n(from )" => s" \1")
    sentences = replace(sentences, r"\n(into )" => s" \1")
    sentences = replace(sentences, r"\n(upon )" => s" \1")
    sentences = replace(sentences, r"\n(after )" => s" \1")
    sentences = replace(sentences, r"\n(while )" => s" \1")
    sentences = replace(sentences, r"\n(during )" => s" \1")
    sentences = replace(sentences, r"\n(within )" => s" \1")
    sentences = replace(sentences, r"\n(through )" => s" \1")
    sentences = replace(sentences, r"\n(between )" => s" \1")
    sentences = replace(sentences, r"\n(whereas )" => s" \1")
    sentences = replace(sentences, r"\n(whether )" => s" \1")

    # no sentence breaks in the middle of specific abbreviations
    sentences = replace(sentences, r"(\be\.)\n(g\.)" => s"\1 \2")
    sentences = replace(sentences, r"(\bi\.)\n(e\.)" => s"\1 \2")
    sentences = replace(sentences, r"(\bi\.)\n(v\.)" => s"\1 \2")

    # no sentence break after specific abbreviations
    sentences = replace(sentences, r"(\be\. ?g\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bi\. ?e\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bi\. ?v\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bvs\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bcf\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bDr\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bMr\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bMs\.)\n" => s"\1 ")
    sentences = replace(sentences, r"(\bMrs\.)\n" => s"\1 ")

    # possible TODO: filter excessively long / short sentences
    sentences
end


text:
using LinearAlgebra
using Random
using Distributions
using Cubature
using Roots
using SpecialFunctions

include("kernels/kernel.jl")
include("kernels/mean_basis.jl")
include("transform/nonlinfun.jl")
include("quadrature/quadgrid.jl")
include("param/param_prior.jl")
include("computation/param_grid_comp.jl")
include("computation/post_prob_comp.jl")
include("computation/weight_comp.jl")
include("computation/distribution_comp.jl")

# Inputs:
    # x0: where to predict, testing data
    # training_data: including features/index/x and labels/values/z
    # kernel: a specific RBF kernel
    # nonlinfun: a specific nonlinear transformation, parameterized
    # p: the degree of freedom using in polynomial basis for the prior mean
    #    p = 1 by default (constant prior mean function)

# Outputs:
    # the distribution of function values z0, at a given index x0
    # including mean, median, standard deviation and condidence interval

function model(x0, training_data, kernel, nonlinfun, p = 1)
    # training_data
    x = training_data[1]
    z = training_data[2]

    # define nonlinear transformation
    # both function values and derivatives
    g =  (x, lambda) -> nonlinfun(x, lambda)[1]
    dg = (x, lambda) -> nonlinfun(x, lambda)[2]

    # define polynomial basis
    X = meanbasisMat(x, p)
    X0 = meanbasisMat(x0, p)

    # precompute the constant coefficient in p(z0|z)
    n = size(x, 1) # #training points
    k = size(x0, 1) # #predicting points
    Gamma = gamma((n-p+k)/2)/(gamma((n-p)/2) * pi^(k/2))

    # wrap up basic training info
    traindata = (idx = x, val = z)
    nonlintrans = (fun = g, deriv = dg)
    trainBasicInfo = (data = training_data, kernel = kernel,
                    polydof = p, nonlintrans = nonlintrans)

    ## compute weights
    # compute quadrature grids and weights
    n_quad = [8, 8, 8] # number of quad nodes for each param, could increase if needed
    n_theta = [n_quad[1], n_quad[2]]
    n_lambda = n_quad[3]
    n_param = (theta = n_theta, lambda = n_lambda)
    param_gridInfo = param_grid_comp(n_param, param_priorInfo)

    # compute alpha_ij = w_i w_j h_ij for theta_i and lambda_j
    alpha, sideSetInfo = weight_comp(param_gridInfo, param_priorInfo, trainBasicInfo)

    # compute distribution of z0 at x0
    distribution_z0 = distribution_comp(x0, trainBasicInfo, sideSetInfo, param_gridInfo, alpha, Gamma)

    return distribution_z0, alpha
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    MaybeHotVector{T, U, V} <: AbstractVector{V}

A vector-like structure for representing one-hot encoded variables.
Like `Flux.OneHotVector` but supports `missing` values.

Construct with the [`maybehot`](@ref) function.

See also: [`MaybeHotMatrix`](@ref), [`maybehotbatch`](@ref).
"""
struct MaybeHotVector{T, U, V} <: AbstractVector{V}
    i::T
    l::U

    MaybeHotVector(i::T, l::U) where {T <: Integer, U <: Integer} = new{T, U, Bool}(i, l)
    MaybeHotVector(i::T, l::U) where {T <: Missing, U <: Integer} = new{T, U, Missing}(i, l)
    MaybeHotVector(i::T, l::U) where {T <: Maybe{Integer}, U <: Integer} = new{T, U, Maybe{T}}(i, l)

    function MaybeHotVector{T, U, V}(i, l) where {T <: Maybe{Integer}, U <: Integer, V <: Maybe{Bool}}
        new{T, U, V}(convert(T, i), convert(U, l))
    end
end

Base.size(x::MaybeHotVector) = (x.l,)
Base.length(x::MaybeHotVector) = x.l
Base.getindex(x::MaybeHotVector, i::Integer) = (@boundscheck checkbounds(x, i); x.i == i)
Base.getindex(x::MaybeHotVector, ::Colon) = MaybeHotVector(x.i, x.l)

function Base.convert(::Type{MaybeHotVector{T, U, W}}, x::MaybeHotVector) where {T, U, W}
    MaybeHotVector{T, U, W}(convert(T, x.i), convert(U, x.l))
end

function Base.promote_rule(::Type{MaybeHotVector{T, U, V}},
        ::Type{MaybeHotVector{A, B, C}}) where {T, A, U, B, V, C}
    MaybeHotVector{promote_type(T, A), promote_type(U, B), promote_type(V, C)}
end

function _check_l(xs::AbstractVecOrTuple{MaybeHotVector})
    l = xs[1].l
    if any(!isequal(l), (x.l for x in xs))
        DimensionMismatch(
            "Number of rows of `MaybeHotMatrix`s to hcat must correspond"
        ) |> throw
    end
    l
end

Base.hcat(xs::T...) where T <: MaybeHotVector = _typed_hcat(T, xs)
Base.hcat(xs::MaybeHotVector...) = _typed_hcat(_promote_types(xs...), xs)
function _typed_hcat(::Type{MaybeHotVector{T, U, V}}, xs::Tuple{Vararg{MaybeHotVector}}) where {T, U, V}
    MaybeHotMatrix{T, U, V}([x.i for x in xs], _check_l(xs))
end

Base.reduce(::typeof(hcat), xs::Vector{<:MaybeHotVector}) = _typed_hcat(mapreduce(typeof, promote_type, xs), xs)
function _typed_hcat(::Type{MaybeHotVector{T, U, V}}, xs::AbstractVector{<:MaybeHotVector}) where {T, U, V}
    MaybeHotMatrix{T, U, V}([x.i for x in xs], _check_l(xs))
end

reduce(::typeof(catobs), as::Vector{<:MaybeHotVector}) = reduce(hcat, as)

A::AbstractMatrix * b::MaybeHotVector = (_check_mul(A, b); _mul(A, b))
Zygote.@adjoint A::AbstractMatrix * b::MaybeHotVector = (_check_mul(A, b); Zygote.pullback(_mul, A, b))

_mul(A::AbstractMatrix, b::MaybeHotVector{Missing}) = fill(missing, size(A, 1))
_mul(A::AbstractMatrix, b::MaybeHotVector{<:Integer}) = A[:, b.i]

Flux.onehot(x::MaybeHotVector{<:Integer}) = Flux.onehot(x.i, 1:x.l)

"""
    maybehot(l, labels)

Return a [`MaybeHotVector`](@ref) where the first occurence of `l` in `labels` is set to `1`
and all other elements are set to `0`.

# Examples
```jldoctest
julia> maybehot(:b, [:a, :b, :c])
3-element MaybeHotVector with eltype Bool:
 ⋅
 1
 ⋅

julia> maybehot(missing, 1:3)
3-element MaybeHotVector with eltype Missing:
 missing
 missing
 missing
```

See also: [`maybehotbatch`](@ref), [`MaybeHotVector`](@ref), [`MaybeHotMatrix`](@ref).
"""
maybehot(::Missing, labels) = MaybeHotVector(missing, length(labels))
function maybehot(l, labels)
    i = findfirst(isequal(l), labels)
    isnothing(i) && ArgumentError("Value $l not in labels $labels") |> throw
    MaybeHotVector(UInt32(i), length(labels))
end

Base.hash(x::MaybeHotVector, h::UInt) = hash((x.i, x.l), h)
(x1::MaybeHotVector == x2::MaybeHotVector) = x1.i == x2.i && x1.l == x2.l
Base.isequal(x1::MaybeHotVector, x2::MaybeHotVector) = isequal(x1.i, x2.i) && x1.l == x2.l


text:
#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic
using IntervalContractors

#Preamble
setprecision(53)
setprecision(Interval, Float64)

# Using approximate checks for validaty update later?
import Base.isapprox
isapprox(x::Interval,y::Interval) = isapprox(x.lo,y.lo,atol=1E-4) && isapprox(x.hi,y.hi,atol=1E-4)

@testset "sinh_rev_test" begin
    @test isapprox(sinh_rev(∅, -∞..∞)[2], ∅)
    @test isapprox(sinh_rev(Interval(-10.0, -1.0), -∞..∞)[2], Interval(-2.99823, -0.881373))
    @test isapprox(sinh_rev(Interval(0.0, Inf), -∞..∞)[2], Interval(0.0, ∞))
    @test isapprox(sinh_rev(Interval(0.0, 1.0), -∞..∞)[2], Interval(0, 0.881374))
    @test isapprox(sinh_rev(Interval(-0.5, 1.0), -∞..∞)[2], Interval(-0.481212, 0.881374))
    @test isapprox(sinh_rev(Interval(-1000.0, 1.0), -∞..∞)[2], Interval(-7.60091, 0.881374))
    @test isapprox(sinh_rev(Interval(0.0, 25.0), -∞..∞)[2], Interval(0.0, 3.91243))
    @test isapprox(sinh_rev(Interval(-1.0, 25.0), -∞..∞)[2], Interval(-0.881374, 3.91243))
end

@testset "cosh_rev_test" begin
    @test isapprox(cosh_rev(∅, -∞..∞)[2], ∅)
    @test isapprox(cosh_rev(Interval(-10.0, -1.0), -∞..∞)[2], ∅)
    @test isapprox(cosh_rev(Interval(0.0, Inf), -∞..∞)[2], Interval(0.0, ∞))
    @test isapprox(cosh_rev(Interval(0.0, 1.0), -∞..∞)[2], Interval(0, 0))
    @test isapprox(cosh_rev(Interval(-0.5, 1.0), -∞..∞)[2], Interval(0, 0))
    @test isapprox(cosh_rev(Interval(-1000.0, 1.0), -∞..∞)[2], Interval(0, 0))
    @test isapprox(cosh_rev(Interval(0.0, 25.0), -∞..∞)[2], Interval(0, 3.91163))
    @test isapprox(cosh_rev(Interval(-1.0, 25.0), -∞..∞)[2], Interval(0, 3.91163))
end

@testset "tanh_rev_test" begin
    @test tanh_rev(∅, -∞..∞)[2] == ∅
    @test tanh_rev(Interval(-10.0, -1.0), -∞..∞)[2] == ∅
    @test isapprox(tanh_rev(Interval(0.0, Inf), -∞..∞)[2], Interval(0.0, ∞))
    @test isapprox(tanh_rev(Interval(0.0, 1.0), -∞..∞)[2], Interval(0.0, ∞))
    @test isapprox(tanh_rev(Interval(-0.5, 1.0), -∞..∞)[2], Interval(-0.549307, ∞))
    @test isapprox(tanh_rev(Interval(-1000.0, 1.0), -∞..∞)[2], Interval(-∞, ∞))
    @test isapprox(tanh_rev(Interval(0.0, 25.0), -∞..∞)[2], Interval(0, ∞))
    @test isapprox(tanh_rev(Interval(-1.0, 25.0), -∞..∞)[2], Interval(-∞, ∞))
end


text:
include("./ffestimator.jl")
using .FFEstimator

using FileIO: load, save, loadstreaming, savestreaming
using Unitful
import LibSndFile
import PortAudio
using Plots

freqs, certs = [], []
xlist = []

open("badapple.raw", "r") do inFile
    startFreq = 200u"Hz"
    endFreq = 800u"Hz"

    samplerate = 48000u"Hz"
    nchannels = 2
    eltype = Float32
    ebom = LE_BOM

    windowSize = trunc(Int, convert(Float64, samplerate * 25u"ms"))
    windowStride = 75

    PortAudio.PortAudioStream(1, 1; eltype=Float32, samplerate=ustrip(samplerate |> u"Hz")) do outFile
        for (period, freq, samplePos, timeStamp, certainty) in
                FFreqEstimator{Float64}(inFile, streamFormat{eltype}(nchannels, samplerate, ebom),
                                        startFreq, endFreq, windowSize, windowStride, 750u"ms")
            push!(freqs, freq)
            push!(xlist, timeStamp)
            push!(certs, certainty)
            write(outFile, Float32.(sin.((samplePos .+ (0:windowStride-1)) ./ period .* (2 * pi))))
        end
    end
end

plot(xlist, freqs, ylim=(200, 800))
display(plot!(twinx(), xlist, certs, linecolor="red"))


text:
using LinearAlgebra

function benchmark1()
	n=30
	eigen(reshape(collect(1.0:n^2), n, n))
	return nothing
end

function benchmark2()
	s = collect(1.0:1e4)
	sort(s .* (-1).^s)
	return nothing
end

function benchmark3()
	nthreads = 4
	arraysize = 1e6
	a = Vector{Vector{Float64}}(UndefInitializer(), nthreads)
	thread = Vector{Task}(UndefInitializer(), nthreads)

	s = collect(1.0:arraysize)
	a[1] = s .* (-1).^s
	for i=2:nthreads a[i]=deepcopy(a[1]) end

	for i = 1 : nthreads
		thread[i] = @async sort!(a[i])
	end
	for i = 1 : nthreads
		wait(thread[i])
	end

	return nothing
end

function runbenchmarks()
	for benchmark in [benchmark1, benchmark2, benchmark3]
		# Warm up run
		benchmark()

		# Report the runtime
		@time benchmark()
	end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b0ded90e-6283-11ec-0039-2de048031f30
begin
	using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
	using MLCourse, CSV, MLJ, DataFrames, StatsPlots, MLJMultivariateStatsInterface,
    Plots, LinearAlgebra, Statistics, Random, MLJLinearModels
end

# ╔═╡ 97101f45-c447-4fd1-bbfd-d4fead047118
md"
# Will it rain tomorrow ?
This notebook shows how we dealt with the training dataset given for the project [Will it rain tomorrow ?]
(https://www.kaggle.com/c/bio322-will-it-rain-tomorrow).
We propose some plots of the data and the transformations we used in our models.

First we need to load the necessary packages. We used the MLCourse env in order to access the packages.
"

# ╔═╡ 6bf0c668-d6e4-474e-aef3-79e5fcde5538
md"
## Loading the packages
"

# ╔═╡ 1f8ded80-aee0-47ee-9617-b6cc53b6c51a
md"
## First look at the datasets
"

# ╔═╡ 75c3f94d-8d0a-4f6f-ab50-41e46541c3d8
md"
We load the datasets;
1. The training data set
2. The test data set
"

# ╔═╡ fbd4c00b-9ffe-4656-a98b-09cc80bf4e58
begin
	training_data = CSV.read(joinpath(@__DIR__, "data", "trainingdata.csv"), DataFrame);
	coerce!(training_data, :precipitation_nextday => Multiclass);
end

# ╔═╡ 29b461df-effd-4e0d-8d48-56e63c07913d
begin
	test_data = CSV.read(joinpath(@__DIR__, "data", "testdata.csv"), DataFrame);
end

# ╔═╡ 0d8b1166-b9c8-44ae-abfa-9060e9a7dd33
md"
We can see that our training data consists of 3176 rows and 529 columns. We have 528 predictors in order to predict the 529th column, whether or not it is raining the next day in Pully called :precipitation\_nextday.

The test data is then composed only of the 528 predictors with different data points. We can use this dataset to predict the :precipitation\_nextday output using our models.
"

# ╔═╡ cb8b6542-1ec9-41a2-8508-bb7acc944e32
size(training_data)

# ╔═╡ ef87c1fb-eaa1-459b-b610-737259a34586
size(test_data)

# ╔═╡ ebdf399c-72f4-4aee-92b2-08727423d174
md"
We can use the function describe(::DataFrames) in order to assess the first properties of each predictor. This can help us find the predictors with mean and standard deviation equal to 0, which will be a problem when standardizing.
"

# ╔═╡ 1b27bf59-2f01-4802-bac9-c325ba1681c3
describe(training_data)

# ╔═╡ 81fb7d32-24dc-460b-85df-aecdd8263e91
md"
## Dealing with missing data
"

# ╔═╡ 406831e5-e18a-40ba-a0da-936774d48425
md"
Dealing with missing data points is of the first steps in a Machine Learning project. The models will need a complete dataset to be fitted on.
"

# ╔═╡ c31e9d13-c4fe-4350-9fdf-3ca7450799d9
md"
When dropping the data points containing missing data we are left with 1699 data points.
"

# ╔═╡ dd4ce6ae-2618-4c07-90b6-19c4eea00a37
begin
	data_drop = dropmissing(training_data)
	coerce!(data_drop, :precipitation_nextday => Multiclass);
	size(data_drop)
end

# ╔═╡ d0a77735-69e1-4093-b153-796b395b82b0
md"
Whereas filling the missing data allows us to keep the the 3176 original points.
"

# ╔═╡ 011a56a8-03a5-4580-b79e-8d1feb69e3f7
begin
	data_med = MLJ.transform(fit!(machine(FillImputer(), select(training_data[:,:], Not([:precipitation_nextday])))), select(training_data[:,:], Not([:precipitation_nextday])));
    insertcols!(data_med,:precipitation_nextday=>training_data.precipitation_nextday[:]);
    coerce!(data_med, :precipitation_nextday => Multiclass);
	size(data_med)
end

# ╔═╡ cbfa4870-8db1-4e5d-ae0e-27d1c2b55066
md"
Is there any good reason to put the median in those empty cells ? It is probably better than putting complete random data, but other ways of filling these cells can be studied.
It still allows us to conserve all the data points and thus a lot of information from the original data set.
"

# ╔═╡ 4593a912-1a19-4e65-bdc5-a777681bbcf6
md"
## Visualization
"

# ╔═╡ 76adc00d-081e-4cd0-92a6-07d2d7599fb3
md"
Looking at Pully measurements in order to observe first result correlation in our data.
"

# ╔═╡ 31286f5d-33cb-4399-9744-139655a4e919
@df data_drop corrplot([:PUY_radiation_1 :PUY_delta_pressure_1 :PUY_air_temp_1 :PUY_wind_1], grid = false, fillcolor = cgrad(), size = (1000, 1000))

# ╔═╡ 1f6ce02d-40c5-43e4-a8a6-889fae7fbd77
md"
## Standardization
"

# ╔═╡ 46da38df-4c18-4767-8325-db55454487e6
md"
In order to standardize our data, we need to find the predictors that have a mean and a standard deviation equal to zero. We can use the fit!(..., verbosity=2) function on the Standardizer() machine. This will give us in our terminal which predictors have a mean and standard deviation equal to 0. Having the mean and standard deviation equal to 0 leads to the fraction 0/0 during the standardization, thus creating NaN cells which are a problem when fitting models."

# ╔═╡ e1b4eb7e-f6a3-4510-908d-922b3211f300
md"
In this step we found that :ZER\_sunshine\_1 and :ALT\_sunshine\_4 have a mean and a standard deviation equal to zero. We then drop these two columns in order to standardize the data."

# ╔═╡ cc06a386-7edb-42fa-a4e2-fa43b2d79779
begin
	mach_drop_std = machine(Standardizer(), select(data_drop, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
	fit!(mach_drop_std, verbosity = 2)
	data_drop_std = MLJ.transform(mach_drop_std, select(data_drop, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
end

# ╔═╡ 40430b02-1a70-45fd-9f14-d446c95f77eb
md"
When we look at the size of the data set, we can see that it has 2 predictors less now, so 526 predictors and 1 ouput column."

# ╔═╡ 4224387b-155f-4bff-944d-4254d54e01a7
size(data_drop_std)

# ╔═╡ 419eb119-613e-431c-8126-a408027084e0
md"
Same thing for the filled out data. "

# ╔═╡ bada701a-79eb-40a0-959d-b5ae148c55a9
begin
	mach_med_std = machine(Standardizer(), select(data_med, Not([:ZER_sunshine_1, :ALT_sunshine_4])))
	fit!(mach_med_std)
	data_med_std = MLJ.transform(mach_med_std, select(data_med, Not([:ZER_sunshine_1, :ALT_sunshine_4]))
)
end

# ╔═╡ 24682b9b-505f-4beb-ae32-eafc58449d89
size(data_med_std)

# ╔═╡ 09682503-f59e-4c84-9bcb-bab37806b71f
md"
## Principal Component Analysis (PCA)
"

# ╔═╡ 90133410-7311-45a0-8962-ad293b3f6e14
md"
We perform PCA on our data set to assess the principal features that explain the variance in :precipitation_nextday.
"

# ╔═╡ 1b03e1d8-8e90-4796-913a-3c5864ec093a
pca_drop = fit!(machine(PCA(pratio = 1), data_drop_std));

# ╔═╡ 89a6c6a2-0c2f-4711-af05-7a905cdb57f4
pca_med = fit!(machine(PCA(pratio = 1), data_med_std));

# ╔═╡ b63d208f-915b-45db-963b-f7f4ef23a2eb
report(pca_med)

# ╔═╡ 95589fe4-42e7-4936-9981-78845cccb1ce
report(pca_drop)

# ╔═╡ a1a8a212-b50f-494e-bb63-79e5544c9267
let
	gr()
	biplot(pca_med)
end

# ╔═╡ 14525726-dd52-4ec1-a68b-51d0cd99ef7a
let
	gr()
	biplot(pca_drop)
end

# ╔═╡ 9861d8a5-dce3-4446-a648-357829bb63a1
md"
Below we can see the proportion and cumulative proportion of variance explained using the dataset with median-imputed missing data.
"

# ╔═╡ 886ac080-398e-4060-b419-28b636e77bac
let
    vars = report(pca_med).principalvars ./ report(pca_med).tvar
    p1 = plot(vars, label = nothing, yscale = :log10,
              xlabel = "component", ylabel = "proportion of variance explained")
    p2 = plot(cumsum(vars),
              label = nothing, xlabel = "component",
              ylabel = "cumulative prop. of variance explained")
    plot(p1, p2, layout = (1, 2), size = (700, 400))
end

# ╔═╡ 9cb38d7b-6065-4756-acfb-f3e888539fdd
md"
Below we can see the proportion and cumulative proportion of variance explained using the dataset with the data points containing missing data removed.
"

# ╔═╡ 5af13796-7a27-4cc5-8cc1-90679500f0ac
let
    vars = report(pca_drop).principalvars ./ report(pca_drop).tvar
    p1 = plot(vars, label = nothing, yscale = :log10,
              xlabel = "component", ylabel = "proportion of variance explained")
    p2 = plot(cumsum(vars),
              label = nothing, xlabel = "component",
              ylabel = "cumulative prop. of variance explained")
    plot(p1, p2, layout = (1, 2), size = (700, 400))
end

# ╔═╡ 267b717b-cc3e-4b1a-ac2c-416ffd9379df
md"
This first introduction to the datasets will allow us to tune our models more precisely. An important function called generate(; option, std, valid, test) allows the transformation of the data in the file Data\_Processing.jl that will be used in every ML method."

# ╔═╡ Cell order:
# ╟─97101f45-c447-4fd1-bbfd-d4fead047118
# ╟─6bf0c668-d6e4-474e-aef3-79e5fcde5538
# ╠═b0ded90e-6283-11ec-0039-2de048031f30
# ╟─1f8ded80-aee0-47ee-9617-b6cc53b6c51a
# ╟─75c3f94d-8d0a-4f6f-ab50-41e46541c3d8
# ╠═fbd4c00b-9ffe-4656-a98b-09cc80bf4e58
# ╠═29b461df-effd-4e0d-8d48-56e63c07913d
# ╟─0d8b1166-b9c8-44ae-abfa-9060e9a7dd33
# ╠═cb8b6542-1ec9-41a2-8508-bb7acc944e32
# ╠═ef87c1fb-eaa1-459b-b610-737259a34586
# ╟─ebdf399c-72f4-4aee-92b2-08727423d174
# ╠═1b27bf59-2f01-4802-bac9-c325ba1681c3
# ╟─81fb7d32-24dc-460b-85df-aecdd8263e91
# ╟─406831e5-e18a-40ba-a0da-936774d48425
# ╟─c31e9d13-c4fe-4350-9fdf-3ca7450799d9
# ╠═dd4ce6ae-2618-4c07-90b6-19c4eea00a37
# ╟─d0a77735-69e1-4093-b153-796b395b82b0
# ╠═011a56a8-03a5-4580-b79e-8d1feb69e3f7
# ╟─cbfa4870-8db1-4e5d-ae0e-27d1c2b55066
# ╟─4593a912-1a19-4e65-bdc5-a777681bbcf6
# ╟─76adc00d-081e-4cd0-92a6-07d2d7599fb3
# ╠═31286f5d-33cb-4399-9744-139655a4e919
# ╟─1f6ce02d-40c5-43e4-a8a6-889fae7fbd77
# ╟─46da38df-4c18-4767-8325-db55454487e6
# ╟─e1b4eb7e-f6a3-4510-908d-922b3211f300
# ╠═cc06a386-7edb-42fa-a4e2-fa43b2d79779
# ╟─40430b02-1a70-45fd-9f14-d446c95f77eb
# ╠═4224387b-155f-4bff-944d-4254d54e01a7
# ╟─419eb119-613e-431c-8126-a408027084e0
# ╠═bada701a-79eb-40a0-959d-b5ae148c55a9
# ╠═24682b9b-505f-4beb-ae32-eafc58449d89
# ╟─09682503-f59e-4c84-9bcb-bab37806b71f
# ╟─90133410-7311-45a0-8962-ad293b3f6e14
# ╠═1b03e1d8-8e90-4796-913a-3c5864ec093a
# ╠═89a6c6a2-0c2f-4711-af05-7a905cdb57f4
# ╠═b63d208f-915b-45db-963b-f7f4ef23a2eb
# ╠═95589fe4-42e7-4936-9981-78845cccb1ce
# ╠═a1a8a212-b50f-494e-bb63-79e5544c9267
# ╠═14525726-dd52-4ec1-a68b-51d0cd99ef7a
# ╟─9861d8a5-dce3-4446-a648-357829bb63a1
# ╟─886ac080-398e-4060-b419-28b636e77bac
# ╟─9cb38d7b-6065-4756-acfb-f3e888539fdd
# ╟─5af13796-7a27-4cc5-8cc1-90679500f0ac
# ╟─267b717b-cc3e-4b1a-ac2c-416ffd9379df

max_stars_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_issues_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_forks_repo_name: clara-rossignol/Will-it-rain-tomorrow


text:
module ExistingProcessManagers

import Distributed

export ExistingProcessManager
export addprocs_existing
export hosts_and_ports_to_workerconfigs
export julia_worker_regex
export parse_julia_worker_output_to_hosts_and_ports

include("types.jl")

include("constructors.jl")

include("addprocs_existing.jl")
include("launch.jl")
include("manage.jl")
include("parse_julia_worker_output.jl")

end # end module ExistingProcessManagers


text:
function check_size(x, y)
    x_s, y_s = size(x, 1), size(y, 1)
    if x_s != y_s
        throw(DimensionMismatch("first dimension of the input data is $x_s, and teaching data is $y_s. these dimensions must match."))
    end
end


text:
alpha = 0.85
list = "test/fixtures/dangling.adj"
order = 6

M, absorbing, order, max_ordinal = readadj(Float64, Int64, alpha, list)

expected = sparse([
 0.0  0.85     0.0      0.0      0.0      0.0
 0.0  0.0      0.425    0.425    0.0      0.0
 0.0  0.28333  0.0      0.28333  0.28333  0.0
 0.0  0.28333  0.28333  0.0      0.28333  0.0
 0.0  0.0      0.28333  0.28333  0.0      0.28333
 0.0  0.0      0.0      0.0      0.0      0.0
])

difference = (max(dense(expected) - dense(M)))

@assert isequal(6, order)
@assert isequal(6, max_ordinal)
@assert isequal(0, difference)
@assert isequal([6], absorbing)


text:
import ..Utils: UnionFind, unite!
import ..IncidenceIndex
import ..num_edges

"""
Provide constant time access to children of a node by index operator.

# Example

```jldoctest
julia> cidx = GraphIdx.Tree.ChildrenIndex([1, 1, 2, 1]);

julia> GraphIdx.Tree.root_node(cidx)
1

julia> collect(cidx[1])  # root node 1 has children 2 and 4
2-element Array{Int64,1}:
 2
 4

julia> collect(cidx[3])  # node 3 has no children
Int64[]

```
"""
struct ChildrenIndex
    idx::Vector{Int}
    value::Vector{Int}
end

@inline Base.getindex(c::ChildrenIndex, j::Int) =
    view(c.value, c.idx[j]:c.idx[j+1]-1)

num_nodes(c::ChildrenIndex) = length(c.value)

"""
    root_node(::ChildrenIndex)

Return the root node of the underlying tree.
This is possible because by convention the root node is stored at
`ChildrenIndex.idx[1]`.
"""
root_node(c::ChildrenIndex) = c.value[1]

Base.length(c::ChildrenIndex) = num_nodes(c)

"""
    ChildrenIndex(n::Int)

Allocate enough space for a tree with `n` nodes.
"""
ChildrenIndex(n::Int) =
    ChildrenIndex(Vector{Int}(undef, n+1), Vector{Int}(undef, n))

"""
    ChildrenIndex(parent, [root])

Construct a new index within the tree given by `parent`.
"""
function ChildrenIndex(parent::Vector{Int}, root::Int = 0)
    idx = ChildrenIndex(length(parent))
    reset!(idx, parent, root)
    return idx
end

"""
    reset!(cidx, parent [,root])

Actually compute the index according to `parent`.
"""
function reset!(cidx::ChildrenIndex, parent::Vector{Int}, root::Int = 0)
    if root <= 0
        root = find_root(parent)
    end
    n = length(parent)
    value = cidx.value
    idx = cidx.idx

    @assert 0 < root <= length(parent) "$root: $(length(parent))"
    @assert parent[root] == root "value[$root] == $(value[root])"
    @assert length(value) == n
    @assert length(idx) == n+1

    idx .= 0
    for p in parent     # compute number of children
        idx[p] += 1
    end
    idx[root] -= 1      # root does not have itself as children

    value[1] = root     # root ==> isperm(value) holds (root is not any child)
    acc = 2
    deg_i = 0
    deg_ii = idx[1]
    for i = 1:n
        idx[i] = acc
        acc += deg_i
        deg_i, deg_ii = deg_ii, idx[i+1]
    end
    idx[n+1] = acc
    @assert deg_ii == 0
    @assert(idx[end] + deg_i == n+1,
            "idx[$(length(idx))] + $deg_i = $(idx[end]) != $(n+1)")

    for (v, p) in enumerate(parent)     # collect children values
        if v == p
            continue
        end
        value[idx[p+1]] = v
        idx[p+1] += 1
    end
    @assert(idx[end] == n+1,
            "idx[$(length(idx))] + $deg_i = $(idx[end]) != $(n+1)")
    return cidx
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
using Alea, FactCheck

facts("Multiindex") do

  context("multiindex()") do
    @fact multiindex(0, 4) --> zeros(1, 0)
    @fact multiindex(5, 0) --> zeros(1,5)
    @fact multiindex(1, 3) --> [0 1 2 3]'
    @fact multiindex(2, 2) --> [0 0; 0 1; 1 0; 0 2; 1 1; 2 0]
  end

  context("multiindex_order()") do
    @fact multiindex_order([0 0; 1 2; 3 5]) --> [0, 3, 8]
  end

end