text:
# Online Outlier Detection
"""
    Euclidean_distance!{tp}(d::Array{tp, 1}, x::AbstractArray{tp, 2}, i::Int, j::Int, dim::Int = 1)

compute the Euclidean distance between x[i,:] and x[j,:] and write the result to d. Memory efficient. dim is the dimension of i and j.
"""
function Euclidean_distance!(d::Array{tp, 1}, x::AbstractArray{tp, 2}, i::Int, j::Int, dim::Int = 1) where {tp}
  d[1] = 0.0
  if dim == 2
    for v = 1:size(x, 1)
      d[1] = d[1] + (x[v, i] - x[v, j]) * (x[v, i] - x[v, j])
    end
  else
    for v = 1:size(x, 2)
      d[1] = d[1] + (x[i, v] - x[j, v]) * (x[i, v] - x[j, v])
    end
  end
  d[1] = sqrt(d[1])
  return d
end

"""
    Mahalanobis_distance!{tp}(d::Array{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, i::Int, j::Int, dim::Int = 1)

compute the Euclidean distance between x[i,:] and x[j,:] and write the result to d. Memory efficient. dim is the dimension of i and j.
"""
function Mahalanobis_distance!(d::Array{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2},i::Int, j::Int, dim::Int = 1) where {tp}
         d[1] = 0.0
         if dim == 2
         for v2 = 1:size(x, 1)
           for v = 1:size(x, 1)
             d[1] = d[1] + (x[v, i] - x[v, j]) * Q[v, v2] *  (x[v2, i] - x[v2, j])
           end
         end
         else
         for v = 1:size(x, 2)
           for v2 = 1:size(x, 2)
             d[1] = d[1] + (x[i, v] - x[j, v]) * Q[v, v2] *  (x[i, v2] - x[j, v2])
           end
         end
         end
         d[1] = sqrt(d[1])
         return d
end

"""
    KDEonline!{tp}(kdescores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2} [, Q::AbstractArray{tp, 2}], σ::tp, dim::Int = 1)

compute (1.0 - Kernel Density Estimates) from x and write it to kdescores with dim being the dimension of the observations.
If teh covariance matrix Q is given, the Mahalanobis distance is used instead of the Euclidean distance.
"""
function KDEonline!(k::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, i::Int, σ::tp, dim::Int = 1) where {tp}
  k[1] = 0.0
  for j = 1:size(x, dim)
    Euclidean_distance!(d, x, i, j, dim)
    innerKDEonline!(k, d, σ)
  end
  finalizeKDEonline!(k, x, dim)
  return k
end

function KDEonline!(k::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, i::Int, σ::tp, dim::Int = 1) where {tp}
  k[1] = 0.0
  for j = 1:size(x, dim)
    Mahalanobis_distance!(d, x, Q, i, j, dim)
    innerKDEonline!(k, d, σ)
  end
  finalizeKDEonline!(k, x, dim)
  return k
end

function innerKDEonline!(k::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, σ::tp) where {tp}
  k[1] = k[1] + exp(-0.5 * d[1] / σ^2)
  return k
end

function finalizeKDEonline!(k::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, dim::Int) where {tp}
  k[1] = 1.0 - (k[1] / size(x, dim))
  return k
end

function KDEonline!(kdescores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, σ::tp, dim::Int = 1) where {tp}
  d = zeros(tp, 1)
  k = zeros(tp, 1)
  for i = 1:size(x, dim)
    KDEonline!(k, d, x, i, σ, dim)
    kdescores[i] = k[1]
  end
  return kdescores
end

function KDEonline!(kdescores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, σ::tp, dim::Int = 1) where {tp}
  d = zeros(tp, 1)
  k = zeros(tp, 1)
  for i = 1:size(x, dim)
    KDEonline!(k, d, x, Q, i, σ, dim)
    kdescores[i] = k[1]
  end
  return kdescores
end

"""
    REConline!{tp}(recscores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2} [, Q::AbstractArray{tp, 2}], ɛ::tp, dim::Int = 1)

compute recurrence scores from x and write it to recscores with dim being the dimension of the observations.
If the covariance matrix Q is given, the mahalanobis distance is used instead of the euclidean distance.
"""
function REConline!(r::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, i::Int, ɛ::tp, dim::Int = 1) where {tp}
    r[1] = 0.0
    for j = 1:size(x, dim)
      Euclidean_distance!(d, x, i, j, dim)
      innerREConline!(r, d, ɛ)
    end
    finalizeREConline!(r, x, dim)
    return r
end

function REConline!(r::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, i::Int, ɛ::tp, dim::Int = 1) where {tp}
    r[1] = 0.0
    for j = 1:size(x, dim)
      Mahalanobis_distance!(d, x, Q, i, j, dim)
      innerREConline!(r, d, ɛ)
    end
    finalizeREConline!(r, x, dim)
    return r
end

function innerREConline!(r::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, ɛ::tp) where {tp}
  if d[1] <= ɛ
    r[1] = r[1] + 1.0
  end
  return r
end

function finalizeREConline!(r::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, dim::Int) where {tp}
  r[1] = 1.0 - ((r[1] - 1.0) / size(x, dim))
  return r
end

function REConline!(recscores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, ɛ::tp, dim::Int = 1) where {tp}
  d = zeros(tp, 1)
  r = zeros(tp, 1)
  for i = 1:size(x, dim)
    REConline!(r, d, x, i, ɛ, dim)
    recscores[i] = r[1]
  end
  return recscores
end

function REConline!(recscores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, ɛ::tp, dim::Int = 1) where {tp}
  d = zeros(tp, 1)
  r = zeros(tp, 1)
  for i = 1:size(x, dim)
    REConline!(r, d, x, Q, i, ɛ, dim)
    recscores[i] = r[1]
  end
  return recscores
end

"""
    KNNonline!{tp}(knnscores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, [Q::AbstractArray{tp, 2},] k::Int, dim::Int = 1)

compute k-nearest neighbor (gamma) scores from x and write it to knnscores with dim being the dimension of the observations.
If the covariance matrix Q is given, the mahalanobis distance is used instead of the euclidean distance.
"""
function KNNonline!(knn::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, i::Int, k::Int, dim::Int = 1) where {tp}
  if size(knn, 1) != k+1 error("size knn != k+1") end
  for l = 1:length(knn)
    knn[l] = 1.0e30#typemax(tp)
  end
  for j = 1:size(x, dim)
    Euclidean_distance!(d, x, i, j, dim)
    innerKNNonline!(knn, d, k)
  end
  finalizeKNNonline!(knn, k)
  return knn
end

function KNNonline!(knn::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, i::Int, k::Int, dim::Int = 1) where {tp}
  if size(knn, 1) != k+1 error("size knn != k+1") end
  for l = 1:length(knn)
    knn[l] = 1.0e30#typemax(tp)
  end
  for j = 1:size(x, dim)
    Mahalanobis_distance!(d, x, Q, i, j, dim)
    innerKNNonline!(knn, d, k)
  end
  finalizeKNNonline!(knn, k)
  return knn
end

function KNNonline!(knnscores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, k::Int, dim::Int = 1) where {tp}
  d = zeros(tp, 1)
  knn = fill(typemax(tp), k+1)
  for i = 1:size(x, dim)
    KNNonline!(knn, d, x, i, k, dim)
    knnscores[i] = knn[1] #mean(knn[2:k+1]) --> in finalize
  end
  return knnscores
end

function KNNonline!(knnscores::AbstractArray{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, k::Int, dim::Int = 1) where {tp}
  d = zeros(tp, 1)
  knn = fill(typemax(tp), k+1)
  for i = 1:size(x, dim)
    KNNonline!(knn, d, x, Q, i, k, dim)
    knnscores[i] = knn[1] #mean(knn[2:k+1]) --> in finalize
  end
  return knnscores
end

function innerKNNonline!(knn::AbstractArray{tp, 1}, d::AbstractArray{tp, 1}, k::Int) where {tp}
    if d[1] <= knn[1]
      for kshift = (k+1):-1:2
        knn[kshift] = knn[kshift-1] # shift old one to k+1
      end
      knn[1] = d[1]   # replace by new distance
    else
      for j = 2:(k+1)
        if d[1] > knn[j-1] && d[1] <= knn[j]
          for kshift = k+1:-1:j
            knn[kshift] = knn[kshift-1] # shift old one to k+1
          end
          knn[j] = d[1]     # replace by new distance
        end
      end
    end
  return knn
end

function finalizeKNNonline!(knn::AbstractArray{tp, 1}, k::Int) where {tp}
  knn[1] = 0.0
  for i = 2:(k+1)
    knn[1] = knn[1] + knn[i]
  end
  knn[1] = knn[1] / k
  return knn
end

"""
    SigmaOnline!{tp}(sigma::Array{tp, 1}, x::AbstractArray{tp, 2}, [Q::AbstractArray{tp, 2}], samplesize::Int = 250, dim::Int = 1)

compute `sigma` parameter as mean of the distances of `samplesize` randomly sampled points along `dim`.
If Q is given the Mahalanobis distance is used instead of Euclidean.
"""
function SigmaOnline!(sigma::Array{tp, 1}, x::AbstractArray{tp, 2}, samplesize::Int = 250, dim::Int = 1) where {tp}
  sigma[1] = 0.0
  d = zeros(tp, 1)
  trainsample1 = zeros(Int, 1)
  trainsample2 = zeros(Int, 1)
  for i = 1:samplesize
    StatsBase.rand!(trainsample1, 1:size(x, dim))[1]
    for j = 1:samplesize
      if j > i
        Euclidean_distance!(d, x, trainsample1[1], StatsBase.rand!(trainsample2, 1:size(x, dim))[1], dim)
        # sum of distances
        sigma[1] = sigma[1] + d[1]
      end
    end
  end
  # get mean
  sigma[1] = sigma[1] / (samplesize^2 / 2 - samplesize / 2)
  return sigma
end

function SigmaOnline!(sigma::Array{tp, 1}, x::AbstractArray{tp, 2}, Q::AbstractArray{tp, 2}, samplesize::Int = 250, dim::Int = 1) where {tp}
  sigma[1] = 0.0
  d = zeros(tp, 1)
  trainsample1 = zeros(Int, 1)
  trainsample2 = zeros(Int, 1)
  for i = 1:samplesize
    StatsBase.rand!(trainsample1, 1:size(x, dim))[1]
    for j = 1:samplesize
      if j > i
        Mahalanobis_distance!(d, x, Q, trainsample1[1], StatsBase.rand!(trainsample2, 1:size(x, dim))[1], dim)
        # sum of distances
        sigma[1] = sigma[1] + d[1]
      end
    end
  end
  # get mean
  sigma[1] = sigma[1] / (samplesize^2 / 2 - samplesize / 2)
  return sigma
end


text:
"""
    validation_input_test(validation::ValidationSettings, benchmark_data::Tuple)

Check whether the settings in `validation` are equal to those in `benchmark_data`.

    validation_input_test(validation_a::ValidationSettings, validation_b::ValidationSettings)

Check whether the settings in `validation_a` are equal to those in `validation_b`.
"""
function validation_input_test(validation::ValidationSettings, benchmark_data::Tuple)

    benchmark_err_type, benchmark_Y, benchmark_n, benchmark_T, benchmark_is_stationary, benchmark_model_struct, benchmark_model_args, benchmark_verb,
        benchmark_verb_estim, benchmark_weights, benchmark_t0, benchmark_subsample, benchmark_max_samples, benchmark_log_folder_path = benchmark_data; # TBD: add test on kwargs

    return ~(false in [validation.err_type == benchmark_err_type;
                      validation.Y === benchmark_Y;
                      validation.n == benchmark_n;
                      validation.T == benchmark_T;
                      validation.is_stationary == benchmark_is_stationary;
                      validation.model_struct == benchmark_model_struct;
                      validation.model_args == benchmark_model_args;
                      validation.verb == benchmark_verb;
                      validation.verb_estim == benchmark_verb_estim;
                      validation.weights == benchmark_weights;
                      validation.t0 == benchmark_t0;
                      validation.subsample == benchmark_subsample;
                      validation.max_samples == benchmark_max_samples;
                      validation.log_folder_path == benchmark_log_folder_path]);
end

function validation_input_test(validation_a::ValidationSettings, validation_b::ValidationSettings)

    return ~(false in [validation_a.err_type == validation_b.err_type;
                      validation_a.Y === validation_b.Y;
                      validation_a.n == validation_b.n;
                      validation_a.T == validation_b.T;
                      validation_a.is_stationary == validation_b.is_stationary;
                      validation_a.model_struct == validation_b.model_struct;
                      validation_a.model_args == validation_b.model_args;
                      validation_a.verb == validation_b.verb;
                      validation_a.verb_estim == validation_b.verb_estim;
                      validation_a.weights == validation_b.weights;
                      validation_a.t0 == validation_b.t0;
                      validation_a.subsample == validation_b.subsample;
                      validation_a.max_samples == validation_b.max_samples;
                      validation_a.log_folder_path == validation_b.log_folder_path]);
end

"""
    grid_input_test(grid::HyperGrid, benchmark_data::Tuple)

Check whether the inputs in `grid` are equal to those in `benchmark_data`.
"""
function grid_input_test(grid::HyperGrid, benchmark_data::Tuple)

    benchmark_p, benchmark_λ, benchmark_α, benchmark_β, benchmark_draws = benchmark_data;

    return ~(false in [grid.p == benchmark_p;
                      grid.λ == benchmark_λ;
                      grid.α == benchmark_α;
                      grid.β == benchmark_β;
                      grid.draws == benchmark_draws]);
end

"""
    validation_settings_tests(model_args::Tuple, weights::FloatVector, t0::Int64, subsample::Float64, max_samples::Int64, log_folder_path::String, benchmark_data::Tuple)

Run basic tests on ValidationSettings.
"""
function validation_settings_tests(model_args::Tuple, weights::FloatVector, t0::Int64, subsample::Float64, max_samples::Int64, log_folder_path::String, benchmark_data::Tuple)

    err_type, Y, n, T, is_stationary, model_struct, _ = benchmark_data;

    # Default optional arguments
    validation_1 = ValidationSettings(err_type, Y, is_stationary, model_struct);
    validation_2 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=());
    validation_3 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true);
    validation_4 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true, verb_estim=false);
    validation_5 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true, verb_estim=false, weights=nothing);
    validation_6 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true, verb_estim=false, weights=nothing, t0=nothing);
    validation_7 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing);
    validation_8 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing);
    validation_9 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=(), verb=true, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_1, benchmark_data);
    @test validation_input_test(validation_1, validation_2);
    @test validation_input_test(validation_1, validation_3);
    @test validation_input_test(validation_1, validation_4);
    @test validation_input_test(validation_1, validation_5);
    @test validation_input_test(validation_1, validation_6);
    @test validation_input_test(validation_1, validation_7);
    @test validation_input_test(validation_1, validation_8);
    @test validation_input_test(validation_1, validation_9);

    # Varying model_args
    validation_10 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args);
    validation_11 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true);
    validation_12 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true, verb_estim=false);
    validation_13 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true, verb_estim=false, weights=nothing);
    validation_14 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true, verb_estim=false, weights=nothing, t0=nothing);
    validation_15 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing);
    validation_16 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing);
    validation_17 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=true, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_10, validation_11);
    @test validation_input_test(validation_10, validation_12);
    @test validation_input_test(validation_10, validation_13);
    @test validation_input_test(validation_10, validation_14);
    @test validation_input_test(validation_10, validation_15);
    @test validation_input_test(validation_10, validation_16);
    @test validation_input_test(validation_10, validation_17);

    # Custom model_args, verb
    validation_18 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false);
    validation_19 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=false);
    validation_20 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=false, weights=nothing);
    validation_21 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=false, weights=nothing, t0=nothing);
    validation_22 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing);
    validation_23 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing);
    validation_24 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=false, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_18, validation_19);
    @test validation_input_test(validation_18, validation_20);
    @test validation_input_test(validation_18, validation_21);
    @test validation_input_test(validation_18, validation_22);
    @test validation_input_test(validation_18, validation_23);
    @test validation_input_test(validation_18, validation_24);

    # Custom model_args, verb, verb_estim
    validation_25 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true);
    validation_26 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=nothing);
    validation_27 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=nothing, t0=nothing);
    validation_28 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=nothing, t0=nothing, subsample=nothing);
    validation_29 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing);
    validation_30 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=nothing, t0=nothing, subsample=nothing, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_25, validation_26);
    @test validation_input_test(validation_25, validation_27);
    @test validation_input_test(validation_25, validation_28);
    @test validation_input_test(validation_25, validation_29);
    @test validation_input_test(validation_25, validation_30);

    # Custom model_args, verb, verb_estim, weights
    validation_31 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights);
    validation_32 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=nothing);
    validation_33 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=nothing, subsample=nothing);
    validation_34 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=nothing, subsample=nothing, max_samples=nothing);
    validation_35 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=nothing, subsample=nothing, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_31, validation_32);
    @test validation_input_test(validation_31, validation_33);
    @test validation_input_test(validation_31, validation_34);
    @test validation_input_test(validation_31, validation_35);

    # Custom model_args, verb, verb_estim, weights, t0
    validation_36 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0);
    validation_37 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=nothing);
    validation_38 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=nothing, max_samples=nothing);
    validation_39 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=nothing, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_36, validation_37);
    @test validation_input_test(validation_36, validation_38);
    @test validation_input_test(validation_36, validation_39);

    # Custom model_args, verb, verb_estim, weights, t0, subsample
    validation_40 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=subsample);
    validation_41 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=subsample, max_samples=nothing);
    validation_42 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=subsample, max_samples=nothing, log_folder_path=nothing);

    @test validation_input_test(validation_40, validation_41);
    @test validation_input_test(validation_40, validation_42);

    # Custom model_args, verb, verb_estim, weights, t0, subsample, max_samples
    validation_43 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=subsample, max_samples=max_samples);
    validation_44 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=subsample, max_samples=max_samples, log_folder_path=nothing);

    @test validation_input_test(validation_43, validation_44);

    # Custom model_args, verb, verb_estim, weights, t0, subsample, max_samples
    validation_45 = ValidationSettings(err_type, Y, is_stationary, model_struct, model_args=model_args, verb=false, verb_estim=true, weights=weights, t0=t0, subsample=subsample, max_samples=max_samples, log_folder_path=log_folder_path);

    @test validation_45.log_folder_path == log_folder_path;
end

@testset "ValidationSettings and HyperGrid: basic functionalities" begin

    # Initialise data and benchmark data
    Y = [0.72 missing 1.86 missing missing 2.52 2.98 3.81 missing 4.36;
         0.95 0.70 missing missing missing missing 2.84 3.88 3.84 4.63];

    n = size(Y,1);
    T = size(Y,2);

    # Tests on ValidationSettings
    benchmark_data_vs = (4, Y, n, T, true, VARSettings, (), true, false, nothing, nothing, nothing, nothing, nothing);
    validation_settings_tests(("something", 123), collect(1.0:0.5:5.5), 100, 0.25, 1000, "some folder", benchmark_data_vs);

    # Tests on HyperGrid
    grid = HyperGrid([1; 10], [0; 0.5], [0; 1], [0; 10], 500);
    @test grid_input_test(grid, ([1; 10], [0; 0.5], [0; 1], [0; 10], 500));
end

@testset "select_hyperparameters: iis and oos" begin

    # Simulate VAR(q) data with mean zero and unit standard deviation
    rng = StableRNG(1);
    q = 2;
    n = 2;
    T = 200;
    true_coefficient = [0.3 0.0 0.0 0.0; 0.0 -0.1 0.0 0.0];

    v = randn(rng, n, T+q);
    Y = zeros(n, T+q);
    for t=1:T+q
        if t <= q
            Y[:,t] .= 0.0;
        else
            Y[:,t] = true_coefficient*vcat([Y[:,t-lag][:] for lag=1:q]...) + v[:,t];
        end
    end
    Y = standardise(Y[:,q+1:end]);

    # Setup in-sample and pseudo out-of-sample validation
    grid = HyperGrid([4, 4], [1, 5], [0, 1], [1, 5], 10);
    validation_settings_iis = ValidationSettings(1, Y, true, VARSettings, verb=false);
    validation_settings_oos = ValidationSettings(2, Y, true, VARSettings, t0=100, verb=false);

    # Run in-sample validation
    candidates_iis, errors_iis = select_hyperparameters(validation_settings_iis, grid, 2);

    # Run pseudo out-of-sample validation
    candidates_oos, errors_oos = select_hyperparameters(validation_settings_oos, grid, 2);

    # Load benchmarks
    benchmark_candidates = read_test_input("./input/validation/candidates");
    benchmark_errors_iis = read_test_input("./input/validation/errors_iis");
    benchmark_errors_oos = read_test_input("./input/validation/errors_oos");

    # Run tests
    @test norm(candidates_iis-benchmark_candidates, Inf) <= 1e-8;
    @test candidates_iis == candidates_oos;
    @test norm(errors_iis-benchmark_errors_iis, Inf) <= 1e-8;
    @test norm(errors_oos-benchmark_errors_oos, Inf) <= 1e-8;
end

ext: jl


text:
"""
  Physical Property Calculator � a software component, which, given the temperature,
  pressure and composition of a Phase of a Material, is able to calculate an additional
  Physical Property or Physical Properties of the Material. Again, a Physical Property
  Calculator will be designed to work with certain kinds of Material. Note that a Physical
  Property Calculator is not called directly by a PME; rather, it is called via a CAPEOPEN Property Package. The purpose of a Physical Property Calculator is to extend or
  to override the list of calculations that a CAPE-OPEN Physical Property Package can
  perform. A CAPE-OPEN Property Calculator can only be used with a Property Package
  which supports the use of Property Calculators.
"""
module PhysicalPropertyCalculator

  using Roots,ECapeExceptions
  export TempPropData, calculate

  type TempPropData
    function TempPropData(prop::String,c::Vector{Float64},compId::UInt16,molweight::Float64)
      new(prop,c,NaN,NaN,compId,molweight,0%UInt8)
    end
    prop::String
    c::Vector{Float64}
    t::Float64
    tc::Float64
    compId::UInt16
    molweight::Float64
    eqno::UInt8
    test::Vector{Float64}
  end

  """
    casRegistryNumber, # Chemical Abstract Service Registry Number
    chemicalFormula, # Chemical formula
    iupacName, # Complete IUPAC Name

    5-  criticalCompressibilityFactor => critical compressibility factor Z
    6-  criticalDensity => critical density in mol/m3
    7-  criticalPressure => critical pressure in Pa
    8-  criticalTemperature => critical temperature in K
    9-  criticalVolume => critical volume in m3/mol
    15- heatOfVaporizationAtNormalBoilingPoint => enthalpy change on vaporization at normal boiling point (101325 Pa) in J/mol
    16- idealGasGibbsFreeEnergyOfFormationAt25C => in J/mol
    17- liquidDensityAt25C => liquid density at 25 ºC in mol/m3
    18- liquidVolumeAt25C => liquid volume at 25 ºC in m3/mol
    20- molecularWeight => relative molar mass
    21- normalBoilingPoint => boiling point temperature at 101325 Pa in K
    24- standardEntropyGas => Standard entropy of gas in J/mol
    28- standardFormationEnthalpyGas => standard enthalpy change on formation of gas in J/mol
    31- standardFormationGibbsEnergyGas => standard Gibbs energy change on formation of gas in J/mol
  """
  function calculate(prop::String, data::Array)
    local floats::Vector{Float64}

    prop=="casRegistryNumber" && (return String(data[1]))
    prop=="chemicalFormula" && (return String(data[2]))
    prop=="iupacName" && (return String(data[3]))

    floats=data[1]
    prop=="criticalCompressibilityFactor" && (return floats[5])
    prop=="criticalDensity" && (return floats[4])
    prop=="criticalPressure" && (return floats[3])
    prop=="criticalTemperature" && (return floats[2])
    prop=="criticalVolume" && (return floats[4])
    prop=="heatOfVaporizationAtNormalBoilingPoint" && (return floats[10])
    prop=="idealGasGibbsFreeEnergyOfFormationAt25C" && (return floats[3])
    prop=="liquidDensityAt25C" && (return floats[10])
    prop=="liquidVolumeAt25C" && (return 1/floats[10])
    prop=="molecularWeight" && (return floats[1])
    prop=="normalBoilingPoint" && (return floats[5])
    prop=="standardEntropyGas" && (return floats[4])
    prop=="standardFormationEnthalpyGas" && (return floats[2])
    prop=="standardFormationGibbsEnergyGas" && (return floats[3])

    throw(ECapeThrmPropertyNotAvailable())
  end

  function calculate(d::TempPropData)
    local tr::Float64
    local ta::Float64
    (d.t<d.test[1] || d.t>d.test[3]) && throw(ECapeOutOfBounds())
    if (d.prop=="heatCapacityOfLiquid")
      if (d.eqno!=2)
        return d.c[1] + d.c[2]*d.t + d.c[3]*d.t^2 + d.c[4]*d.t^3 + d.c[5]*d.t^4  #???? check perry
      else
        tr=1-d.t/d.tc
        return (d.c[1]^2)/tr+d.c[2]-2*d.c[1]*d.c[3]*tr-d.c[1]*d.c[4]*tr^2-(d.c[3]^2*tr^3)/3-(d.c[3]*d.c[4]*tr^4)/2-(d.c[4]^2*tr^5)/5
      end
    elseif (d.prop=="heatOfVaporization")
      # Tr = T/Tc. Heat of vaporization in J/kmol
      # perry 2_150 have not presented c5 although it was a part of formula
      tr=d.t/d.tc
      return d.c[1]*1e7*(1-tr)^(d.c[2]+d.c[3]*tr+d.c[4]*tr^2)
    elseif (d.prop=="idealGasHeatCapacity")
      if (d.eqno==1)
        return d.c[1]+d.c[2]*d.t+d.c[3]*d.t^2+d.c[4]*1e-5*d.t^3+d.c[5]*1e-10*d.t^4
      elseif (d.eqno==2)
        local c1::Float64
        local c2::Float64
        local c3::Float64
        local c4::Float64
        c1=d.c[1]*1e5
        c2=d.c[2]*1e5
        c3=d.c[3]*1e3
        c4=d.c[4]*1e5
        return c1+c2*((c3/d.t)/sinh(c3/d.t))^2+c4*((d.c[5]/d.t)/cosh(d.c[5]/d.t))^2
      end
    elseif (d.prop=="idealGasEnthalpy")
      if(d.eqno==1)
        return d.c[1]*d.t+1/2*d.c[2]*d.t^2+1/3*d.c[3]*d.t^3+1/4*d.c[4]*d.t^4+1/5*d.c[5]*d.t^5 # Integ of Cp Poly
      elseif(d.eqno==2)
        return d.c[1]*d.t+d.c[2]*d.c[3]*coth(d.c[3]/d.t)-d.c[4]*d.c[5]*tanh(d.c[5]/d.t) # Integ of Cp Hyper
      end
    elseif (d.prop=="idealGasEntropy")
      if(d.eqno==1)
        return d.c[2]*d.t+(d.c[3]*d.t^2)/2+(d.c[4]*d.t^3)/3+(d.c[5]*d.t^4)/4+d.c[1]*log(d.t) # Integral of Cp/T Poly
      elseif(d.eqno==2)
        return (d.c[2]*d.c[3]*coth(d.c[3]/d.t)+d.c[1]*d.t*log(d.t)+d.c[4]*d.t*log(cosh(d.c[5]/d.t))-d.c[2]*d.t*log(sinh(d.c[3]/d.t))-d.c[4]*d.c[5]*tanh(d.c[5]/d.t))/d.t # Integral of Cp/T Hyper
      end
    elseif (d.prop=="thermalConductivityOfLiquid")
      # Thermal conductivites are at either 1 atm or the vapor pressure, whichever is higher.
      return d.c[1]+d.c[2]*d.t+d.c[3]*d.t^2+d.c[4]*d.t^3+d.c[5]*d.t^4
    elseif (d.prop=="thermalConductivityOfVapor")
      # Thermal conductivites are at either 1 atm or the vapor pressure, whichever is lower.
      return d.c[1]*(d.t^d.c[2])/(1+d.c[3]/d.t+d.c[4]/(d.t^2))
    elseif (d.prop=="vaporPressure")
      # Vapor pressure in Pa.
      return exp(d.c[1]+d.c[2]/d.t+d.c[3]*log(d.t)+d.c[4]*d.t^d.c[5])
    elseif (d.prop=="viscosityOfLiquid")
      #  Viscosities are at either 1 atm or the vapor pressure, whichever is higher.
      return exp(d.c[1]+d.c[2]/d.t+d.c[3]*log(d.t)+d.c[4]*d.t^d.c[5])
    elseif (d.prop=="viscosityOfVapor")
      # Viscosities are at either 1 atm or the vapor pressure, whichever is lower. in Pa.
      return d.c[1]*(d.t^d.c[2])/(1+d.c[3]/d.t+d.c[4]/(d.t^2))
    elseif (d.prop=="volumeOfLiquid")
      #Liquid dencity
      local ld::Float64
      if (d.eqno==2)
        ld=d.c[1]+d.c[2]*d.t+d.c[3]*d.t^2+d.c[4]*d.t^3 # o-terphenyl and water limited range
      end
      if (d.eqno==3) # For water over the entire temperature range of 273.16 to 647.096 K.
        ta=1-(d.t/647.096)
        ld=17.863+58.606*ta^0.35 - 95.396*ta^(2/3)+213.89*ta- 141.26*ta^(4/3)
      end
      ld=d.c[1]/(d.c[2]^(1+(1-d.t/d.c[3])^d.c[4]))  # The others
      return 1/ld
    end
    throw(ECapeThrmPropertyNotAvailable())
  end

  function calculate(d::TempPropData, p::Float64)
    (p<d.test[2] || p>d.test[4]) && throw(ECapeOutOfBounds())
    if (d.prop=="boilingPointTemperature")
      bp=fzeros(t->vp(d.c,p,t),d.test[1],d.test[3])
      (length(bp)==1) && return bp[1]
    end
    throw(ECapeThrmPropertyNotAvailable())
  end
  function vp(c::Vector{Float64}, p::Float64, t::Float64)
    return exp(c[1] + c[2]/t + c[3]*log(t) + c[4]*t^c[5])-p
  end

end


text:
using Base.Iterators: partition
using Random
using BSON
using Statistics
using MLDatasets
using Flux
using Flux: onehotbatch, onecold, crossentropy
using Flux.Data: DataLoader
using Plots
using ImageInspector
import RDatasets: dataset

# Introduction to Flux

include("utilities.jl")

Random.seed!(666)

iris = dataset("datasets", "iris")

X = Matrix(iris[:, 1:4])
y = iris.Species

X_train, y_train, X_test, y_test, classes = prepare_data(X', y; dims=2)

n_hidden = 5
m = Chain(
    Dense(size(X_train,1), n_hidden, relu),
    Dense(n_hidden, size(y_train,1), identity),
    softmax,
)

m(X_train)

params(m[2])[2] .= [-1;0;1]

L(x,y) = crossentropy(m(x), y)

L(X_train, y_train)

ps = params(m)
grad = gradient(() -> L(X_train, y_train), ps)

grad = gradient(() -> L(X_train, y_train), params(X_train))

size(grad[X_train])

opt = Descent(0.1)
max_iter = 250

acc_test = zeros(max_iter)
for i in 1:max_iter
    gs = gradient(() -> L(X_train, y_train), ps)
    Flux.Optimise.update!(opt, ps, gs)
    acc_test[i] = accuracy(X_test, y_test)
end

plot(acc_test, xlabel="Iteration", ylabel="Test accuracy", label="", ylim=(-0.01,1.01))

# Loading data

T = Float32
X_train, y_train = MLDatasets.MNIST.traindata(T)
X_test, y_test = MLDatasets.MNIST.testdata(T)

# Exercise

# Exercise

# Loading data

function load_data(dataset; T=Float32, onehot=false, classes=0:9)
    X_train, y_train = dataset.traindata(T)
    X_test, y_test = dataset.testdata(T)

    X_train = reshape_data(X_train)
    X_test = reshape_data(X_test)

    if onehot
        y_train = onehotbatch(y_train, classes)
        y_test = onehotbatch(y_test, classes)
    end

    return X_train, y_train, X_test, y_test
end

X_train, y_train, X_test, y_test = load_data(MLDatasets.MNIST; T=T, onehot=true)

# Exercise

# Exercise

# Bonus

batches = map(partition(randperm(size(y, 2)), batchsize)) do inds
    return (X[:, :, :, inds], y[:, inds])
end

[(X[:, :, :, inds], y[:, inds]) for inds in partition(randperm(size(y, 2)), batchsize)]

# Define model

Random.seed!(666)
m = Chain(
    Conv((2,2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, relu),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)),
    softmax,
)

L(X, y) = crossentropy(m(X), y)

# Train model

function train_model!(m, L, X, y;
        opt = Descent(0.1),
        batchsize = 128,
        n_epochs = 10,
        file_name = "")

    batches = DataLoader((X, y); batchsize, shuffle = true)

    for _ in 1:n_epochs
        Flux.train!(L, params(m), batches, opt)
    end

    !isempty(file_name) && BSON.bson(file_name, m=m)

    return
end

# Exercise

# Exercise

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
struct Rocket{T} <: Model{T}
      nq::Int
      nu::Int
      nw::Int
      nc::Int

      mass::T                               # mass
      inertia::Diagonal{T,Vector{T}}        # inertia matrix
      inertia_inv::Diagonal{T,Vector{T}}    # inertia matrix inverse
      gravity::Vector{T}                    # gravity
      length::T                             # length (com to thruster)
end

function f(model::Rocket, z, u, w)
      # states
      x = z[1:3]
      r = z[4:6]
      v = z[7:9]
      ω = z[10:12]

      # force in body frame
      F = u[1:3]

      # torque in body frame
      τ = [model.length * u[2],
           -model.length * u[1],
           0.0]

      [v;
       0.25 * ((1.0 - r' * r) * ω - 2.0 * cross(ω, r) + 2.0*(ω' * r) * r);
       model.gravity + (1.0 / model.mass) * MRP(r[1], r[2], r[3]) * F;
       model.inertia_inv * (τ - cross(ω, model.inertia * ω))]
end

nq, nu, nw = 12, 3, 0

mass = 1.0
len = 1.0

inertia = Diagonal([1.0 / 12.0 * mass * len^2.0, 1.0 / 12.0 * mass * len^2.0, 1.0e-5])
inertia_inv = Diagonal([1.0 / (1.0 / 12.0 * mass * len^2.0), 1.0 / (1.0 / 12.0 * mass * len^2.0), 1.0 / (1.0e-5)])

rocket = Rocket(nq, nu, nw, 0,
                  mass,
                  inertia,
                  inertia_inv,
                  [0.0, 0.0, -9.81],
                  len)