text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestConstraintMap

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

struct ConstraintDummyBridge <: MOI.Bridges.Constraint.AbstractBridge
    id::Int
end

function test_map()
    map = MOI.Bridges.Constraint.Map()

    @test isempty(map)
    @test !MOI.Bridges.Constraint.has_bridges(map)

    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    b1 = ConstraintDummyBridge(1)
    c1 = MOI.Bridges.Constraint.add_key_for_bridge(map, b1, x, MOI.EqualTo(0.0))

    @test c1.value == x.value
    @test haskey(map, c1)
    @test map[c1] == b1
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c1)) == 1
    @test collect(MOI.Bridges.Constraint.keys_of_type(map, typeof(c1))) == [c1]

    @test length(map) == 1
    @test !isempty(map)
    @test MOI.Bridges.Constraint.has_bridges(map)

    b2 = ConstraintDummyBridge(2)
    vov = MOI.VectorOfVariables([x, y])
    c2 = MOI.Bridges.Constraint.add_key_for_bridge(
        map,
        b2,
        vov,
        MOI.SecondOrderCone(2),
    )
    @test c2.value == -1
    @test haskey(map, c2)
    @test map[c2] == b2
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c2)) == 1
    @test collect(MOI.Bridges.Constraint.keys_of_type(map, typeof(c2))) == [c2]

    @test length(map) == 2
    @test !isempty(map)
    @test MOI.Bridges.Constraint.has_bridges(map)

    b3 = ConstraintDummyBridge(3)
    c3 = MOI.Bridges.Constraint.add_key_for_bridge(
        map,
        b3,
        1.0x + 2.0,
        MOI.EqualTo(0.0),
    )
    @test haskey(map, c3)
    @test map[c3] == b3
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c3)) == 1
    @test collect(MOI.Bridges.Constraint.keys_of_type(map, typeof(c3))) == [c3]

    @test length(map) == 3
    @test !isempty(map)
    @test MOI.Bridges.Constraint.has_bridges(map)

    bridges = collect(values(map))
    @test sort([b.id for b in bridges]) == 1:3
    elements = sort(collect(map), by = el -> el.second.id)
    @test elements[1].first == c1
    @test elements[1].second == b1
    @test elements[2].first == c2
    @test elements[2].second == b2
    @test elements[3].first == c3
    @test elements[3].second == b3

    @test MOI.Bridges.Constraint.list_of_key_types(map) == Set([
        (MOI.VectorOfVariables, MOI.SecondOrderCone),
        (MOI.VariableIndex, MOI.EqualTo{Float64}),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
    ])
    @test MOI.Bridges.Constraint.variable_constraints(map, x) == [c1]
    @test isempty(MOI.Bridges.Constraint.variable_constraints(map, y))

    delete!(map, c1)
    @test length(map) == 2
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 2 constraint bridges"
    @test !isempty(map)
    @test MOI.Bridges.Constraint.has_bridges(map)

    @test MOI.Bridges.Constraint.list_of_key_types(map) == Set([
        (MOI.VectorOfVariables, MOI.SecondOrderCone),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
    ])
    @test isempty(MOI.Bridges.Constraint.variable_constraints(map, x))
    @test isempty(MOI.Bridges.Constraint.variable_constraints(map, y))
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c1)) == 0
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c2)) == 1
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c3)) == 1

    delete!(map, c2)
    @test length(map) == 1
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 1 constraint bridge"
    @test !isempty(map)
    @test MOI.Bridges.Constraint.has_bridges(map)

    @test MOI.Bridges.Constraint.list_of_key_types(map) ==
          Set([(MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})])
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c1)) == 0
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c2)) == 0
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c3)) == 1

    delete!(map, c3)
    @test length(map) == 0
    @test isempty(map)
    @test MOI.Bridges.Constraint.has_bridges(map)

    @test isempty(MOI.Bridges.Constraint.list_of_key_types(map))
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c1)) == 0
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c2)) == 0
    @test MOI.Bridges.Constraint.number_of_type(map, typeof(c3)) == 0

    @test isempty(collect(map))

    map = MOI.Bridges.Constraint.EmptyMap()
    empty!(map)
    @test isempty(map)
    @test isempty(keys(map))
    @test isempty(values(map))
    @test !MOI.Bridges.Constraint.has_bridges(map)
    @test iszero(MOI.Bridges.Constraint.number_of_type(map, typeof(c1)))
    @test iszero(MOI.Bridges.Constraint.number_of_type(map, typeof(c2)))
    @test iszero(MOI.Bridges.Constraint.number_of_type(map, typeof(c3)))
    @test sprint(MOI.Bridges.print_num_bridges, map) == ""
end

end  # module

TestConstraintMap.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
#types
export Term, LinearOperatorFamily, Solution, save
#solvers
export householder, beyn, solve, count_poles_and_zeros
export mslp, inveriter, lancaster, traceiter, rf2s, decode_error_flag
export nicoud, picard
#export mehrmann, juniper, lancaster, count_eigvals
#helper functions for beyn
export pos_test, moments2eigs, compute_moment_matrices, wn, inpoly
#perturbation stuff
export pade!, perturb!, perturb_fast!, perturb_norm!, estimate_pol, conv_radius
#export algebra
export pow0,pow1,pow2,pow_a,pow,exp_delay,z_exp_iaz,z_exp__iaz,exp_pm,exp_az2mzit
export Σnexp_az2mzit, exp_az, generate_stsp_z, generate_z_g_z
export generate_Σy_exp_ikx,generate_gz_hz, generate_exp_az
# export fitting staff
#export fit_ss_wrapper

ext: jl


text:
# # Machine Learning in Julia (continued)

# An introduction to the
# [MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
# toolbox.

# ### Set-up

# Inspect Julia version:

VERSION

# The following instantiates a package environment.

# The package environment has been created using **Julia 1.6** and may not
# instantiate properly for other Julia versions.

using Pkg
Pkg.activate("env")
Pkg.instantiate()

# ## General resources

# - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
# - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
# - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
# - [Data Science Tutorials in Julia](https://juliaai.github.io/DataScienceTutorials.jl/)

# ## Part 5 - Advanced Model Composition

# > **Goals:**
# > 1. Learn how to build a prototypes of a composite model, called a *learning network*
# > 2. Learn how to use the `@from_network` macro to export a learning network as a new stand-alone model type

# Pipelines are great for composing models in an unbranching
# sequence. Another built-in type of model composition is a model
# *stack*; see
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/model_stacking/#Model-Stacking)
# for details. For other more complicated model compositions you'll want to
# use MLJ's generic model composition syntax. There are two main
# steps:

# - **Prototype** the composite model by building a *learning
#   network*, which can be tested on some (dummy) data as you build
#   it.

# - **Export** the learning network as a new stand-alone model type.

# Like pipeline models, instances of the exported model type behave
# like any other model (and are not bound to any data, until you wrap
# them in a machine).

# ### Building a pipeline using the generic composition syntax

# To warm up, we'll do the equivalent of

using MLJ
LogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels
pipe = Standardizer |> LogisticClassifier;

# using the generic syntax.

# Here's some dummy data we'll be using to test our learning network:

X, y = make_blobs(5, 3)
pretty(X)

# **Step 0** - Proceed as if you were combining the models "by hand",
# using all the data available for training, transforming and
# prediction:

stand = Standardizer();
linear = LogisticClassifier();

mach1 = machine(stand, X);
fit!(mach1);
Xstand = transform(mach1, X);

mach2 = machine(linear, Xstand, y);
fit!(mach2);
yhat = predict(mach2, Xstand)

# **Step 1** - Edit your code as follows:

# - pre-wrap the data in `Source` nodes

# - delete the `fit!` calls

X = source(X)  # or X = source() if not testing
y = source(y)  # or y = source()

stand = Standardizer();
linear = LogisticClassifier();

mach1 = machine(stand, X);
Xstand = transform(mach1, X);

mach2 = machine(linear, Xstand, y);
yhat = predict(mach2, Xstand)

# Now `X`, `y`, `Xstand` and `yhat` are *nodes* ("variables" or
# "dynammic data") instead of data. All training, predicting and
# transforming is now executed lazily, whenever we `fit!` one of these
# nodes. We *call* a node to retrieve the data it represents in the
# original manual workflow.

fit!(Xstand)
Xstand() |> pretty

#-

fit!(yhat);
yhat()

# The node `yhat` is the "descendant" (in an associated DAG we have
# defined) of a unique source node:

sources(yhat)

#-

# The data at the source node is replaced by `Xnew` to obtain a
# new prediction when we call `yhat` like this:

Xnew, _ = make_blobs(2, 3);
yhat(Xnew)

# **Step 2** - Export the learning network as a new stand-alone model type

# Now, somewhat paradoxically, we can wrap the whole network in a
# special machine - called a *learning network machine* - before have
# defined the new model type. Indeed doing so is a necessary step in
# the export process, for this machine will tell the export macro:

# - what kind of model the composite will be (`Deterministic`,
#   `Probabilistic` or `Unsupervised`)a

# - which source nodes are input nodes and which are for the target

# - which nodes correspond to each operation (`predict`, `transform`,
#   etc) that we might want to define

surrogate = Probabilistic()     # a model with no fields!
mach = machine(surrogate, X, y; predict=yhat)

# Although we have no real need to use it, this machine behaves like
# you'd expect it to:

Xnew, _ = make_blobs(2, 3)
fit!(mach)
predict(mach, Xnew)

#-

# Now we create a new model type using a Julia `struct` definition
# appropriately decorated:

@from_network mach begin
    mutable struct YourPipe
        standardizer = stand
        classifier = linear::Probabilistic
    end
end

# Instantiating and evaluating on some new data:

pipe = YourPipe()
X, y = @load_iris;   # built-in data set
mach = machine(pipe, X, y)
evaluate!(mach, measure=misclassification_rate, operation=predict_mode)

# ### A composite model to average two regressor predictors

# The following is condensed version of
# [this](https://github.com/alan-turing-institute/MLJ.jl/blob/master/binder/MLJ_demo.ipynb)
# tutorial. We will define a composite model that:

# - standardizes the input data

# - learns and applies a Box-Cox transformation to the target variable

# - blends the predictions of two supervised learning models - a ridge
#  regressor and a random forest regressor; we'll blend using a simple
#  average (for a more sophisticated stacking example, see
#  [here](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/stacking/))

# - applies the *inverse* Box-Cox transformation to this blended prediction

RandomForestRegressor = @load RandomForestRegressor pkg=DecisionTree
RidgeRegressor = @load RidgeRegressor pkg=MLJLinearModels

# **Input layer**

X = source()
y = source()

# **First layer and target transformation**

std_model = Standardizer()
stand = machine(std_model, X)
W = MLJ.transform(stand, X)

box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, y)
z = MLJ.transform(box, y)

# **Second layer**

ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)

forest_model = RandomForestRegressor(n_trees=50)
forest = machine(forest_model, W, z)

ẑ = 0.5*predict(ridge, W) + 0.5*predict(forest, W)

# **Output**

ŷ = inverse_transform(box, ẑ)

# With the learning network defined, we're ready to export:

@from_network machine(Deterministic(), X, y, predict=ŷ) begin
    mutable struct CompositeModel
        rgs1 = ridge_model
        rgs2 = forest_model
    end
end

# Let's instantiate the new model type and try it out on some data:

composite = CompositeModel()

#-

X, y = @load_boston;
mach = machine(composite, X, y);
evaluate!(mach,
          resampling=CV(nfolds=6, shuffle=true),
          measures=[rms, mae])

# ### Resources for Part 5
#
# - From the MLJ manual:
#    - [Learning Networks](https://alan-turing-institute.github.io/MLJ.jl/stable/composing_models/#Learning-Networks-1)
# - From Data Science Tutorials:
#     - [Learning Networks](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/learning-networks/)
#     - [Learning Networks 2](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/learning-networks-2/)
#     - [Stacking](https://juliaai.github.io/DataScienceTutorials.jl/getting-started/stacking/): an advanced example of model composition
#     - [Finer Control](https://alan-turing-institute.github.io/MLJ.jl/dev/composing_models/#Method-II:-Finer-control-(advanced)-1):
#       exporting learning networks without a macro for finer control

#


text:
module FlowPattern

using Blink
using Plots
using Interact

const default_exec_file = joinpath(@__DIR__, "..\\external\\Flowpat.exe")

function calculate_flow_pattern(;
    USᴳ = 0.070, # [m/s]
    USᴸ = 51.70, # [m/s]
    μᴳ = 0.000018, # [Pa.s]
    μᴸ = 0.005, # [Pa.s]
    ρᴳ = 2.8833234073128, # [kg/m³]
    ρᴸ = 997.950268197708, # [kg/m³]
    D = 0.0508, # [m]
    ε = 0, # [m]
    σ = 7.40500030517578E-02, # [N/m]
    θ = 0, # [°]
    exec_file = default_exec_file,
)
    flowpat = `$exec_file`
    input = open("Input.txt", "w")
    flowpat_mode = 2
    println(input, "$flowpat_mode\tMode 1-Transition Lines, 2-FlowPattern Prediction")
    println(input, "$USᴸ\t$USᴳ\t$ρᴸ\t$μᴸ\t$ρᴳ\t$μᴳ\t$σ\t$D\t$θ\t$ε")
    close(input)
    run(flowpat)
    flow_pattern_id, error_code = parse.(Int, split(read("Output.txt", String))[end-1:end])

    if error_code > 0
        return nothing
    end

    flowpattern = ("SS", "SW", "A", "DB", "B", "SL")
    return flowpattern[flow_pattern_id]
end

function calculate_flow_pattern_transitions(
    μᴳ,
    μᴸ,
    ρᴳ,
    ρᴸ,
    D,
    ε,
    σ,
    θ,
    interface;
    exec_file = default_exec_file,
)
    n_transitions = 5
    flowpat = `$exec_file`
    input = open("Input.txt", "w")
    flowpat_mode = 1
    println(input, "$flowpat_mode\tMode 1-Transition Lines, 2-FlowPattern Prediction")
    println(input, "$ρᴸ\tLiquid Density [Kg/m^3]")
    println(input, "$μᴸ\tLiquid Viscosity [Pa s]")
    println(input, "$σ\tSurface Tension [N/m]")
    println(input, "$ρᴳ\tGas Density [Kg/m^3]")
    println(input, "$μᴳ\tGas Viscosity [Pa s]")
    println(input, "$D\tPipe Diameter [m]")
    println(input, "$θ\tPipe inclination  Angle [deg]")
    println(input, "$ε\tPipe Roughness [m]")
    println(input, "$interface\tInterface 1-Smooth, 2-Wavy")
    close(input)
    run(flowpat)

    function init_xy(size)
        x = Vector{Vector{Float64}}()
        y = Vector{Vector{Float64}}()
        for i = 1:size
            push!(x, Vector{Float64}())
            push!(y, Vector{Float64}())
        end
        return x, y
    end

    x, y = init_xy(n_transitions)
    function push_if_available_sync!(arr_x, arr_y, val_x, val_y)
        if val_x > -999 && val_y > -999
            push!(arr_x, val_x)
            push!(arr_y, val_y)
        end
    end
    open("Output.txt") do file
        for (i, ln) in enumerate(eachline(file))
            result = parse.(Float64, split(ln))'
            for j = 1:n_transitions
                push_if_available_sync!(x[j], y[j], result[2j-1], result[2j])
            end
        end
    end

    return x, y
end

function plot_flow_pattern_transitions(θ, D, ρᴳ, ρᴸ, μᴳ, μᴸ, σ, ε, interface)
    x, y = calculate_flow_pattern_transitions(μᴳ, μᴸ, ρᴳ, ρᴸ, D, ε, σ, θ, interface)
    p = plot(
        x,
        y,
        xlabel = "USᴳ [m/s]",
        ylabel = "USᴸ [m/s]",
        xlims = (1e-3, 100),
        ylims = (1e-3, 100),
        xaxis = :log,
        yaxis = :log,
        label = ["Stratified" "Annular" "Wavy" "Dispersed-Bubble" "Bubble"],
        legend = :outertopright,
        grid = :all,
        minorgrid = true,
        formatter = string, # A method which converts a number to a string for tick labeling.
    )
end

function flowpattern_main()
    ##
    # Create sliders and plot
    θ_slider = widget(0:0.5:90, label = "θ [°]"; value = 0.0)
    D_slider = widget(0.01:0.01:1, label = "D [m]"; value = 0.1)
    ρᴳ_slider = widget(0.5:0.1:1100, label = "ρᴳ [kg/m³]"; value = 1.0)
    ρᴸ_slider = widget(0.5:0.1:1100, label = "ρᴸ [kg/m³]"; value = 1000.0)
    μᴳ_slider = widget(1e-5:1e-5:1e-1, label = "μᴳ [Pa.s]"; value = 1e-5)
    μᴸ_slider = widget(1e-5:1e-5:1e-1, label = "μᴸ [Pa.s]"; value = 1e-3)
    σ_slider = widget(1e-3:1e-3:1, label = "σ [N/m]"; value = 1e-2)
    ε_slider = widget(0:1e-5:1e-2, label = "ε [m]"; value = 0)
    interface_selector = widget(Dict("Stratified" => 1, "Wavy" => 2), label = "Interface")

    interactive_plot = map(
        plot_flow_pattern_transitions,
        θ_slider,
        D_slider,
        ρᴳ_slider,
        ρᴸ_slider,
        μᴳ_slider,
        μᴸ_slider,
        σ_slider,
        ε_slider,
        interface_selector,
    )

    ## Create UI layout
    el1 = hbox(
        vbox(θ_slider, ρᴳ_slider, μᴳ_slider, ε_slider),
        vbox(D_slider, ρᴸ_slider, μᴸ_slider, σ_slider),
    )
    el2 = interface_selector
    el3 = Interact.hline()
    el4 = interactive_plot

    ui = vbox(el1, el2, el3, el4)

    ## Create final window
    w = Window()
    size(w, 700, 700)
    title(w, "Flow Pattern Map Calculator")
    body!(w, ui)

    if "--quiet" in ARGS
        close(w)
    else
        # Workaround for keeping Blink window alive when calling `julia <my_script>.jl`
        while active(w)
            sleep(1)
        end
    end
end

function julia_main()
    try
        flowpattern_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

end # module FlowPattern

max_stars_repo_name: arthursoprana/FlowPattern

max_issues_repo_name: arthursoprana/FlowPattern

max_forks_repo_name: arthursoprana/FlowPattern


text:
"""
Convert a FrequencySimulationResult into a TimeSimulationResult by using the inverse fourier transform.
Assumes only positive frequencies and a real time signal
"""
function frequency_to_time(simres::FrequencySimulationResult{T,Dim,FieldDim};
        t_vec::AbstractVector{T} = ω_to_t(simres.ω),
        impulse::ContinuousImpulse{T} = TimeDiracImpulse(zero(T)), #GaussianImpulse(maximum(simres.ω)),
        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse, t_vec, simres.ω),
        method = :dft
    ) where {Dim,FieldDim,T}

    t_vec = discrete_impulse.t

    time_field = frequency_to_time(transpose(field(simres)), simres.ω, t_vec;
        discrete_impulse = discrete_impulse, method = method)

    return TimeSimulationResult(transpose(time_field), simres.x, t_vec)
end

"""
Convert a TimeSimulationResult into a FrequencySimulationResult by using the fourier transform.
Assumes only positive frequencies and a real time signal
"""
function time_to_frequency(timres::TimeSimulationResult{T,Dim,FieldDim};
        ω_vec = t_to_ω(timres.t),
        impulse::ContinuousImpulse{T} = TimeDiracImpulse(zero(T)), #GaussianImpulse(maximum(ω_vec)),
        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse,timres.t, ω_vec),
        method =:dft
    ) where {Dim,FieldDim,T}

    ω_vec = discrete_impulse.ω

    freq_field = time_to_frequency(transpose(field(timres)), timres.t, ω_vec;
        discrete_impulse = discrete_impulse, method = method)

    return FrequencySimulationResult(transpose(freq_field), deepcopy(timres.x), ω_vec)
end

"""
returns an array of time from the frequency array ω_vec.
Uses the same convention for sampling the time as the discrete Fourier transfrom.
Assumes ω_vec is ordered and non-negative.
"""
function ω_to_t(ω_arr::AbstractArray{T}) where T <: AbstractFloat
    N = length(ω_arr)
    if ω_arr[1] == zero(T)
        N -= 1
    elseif minimum(ω_arr) < zero(T)
        error("expected only non-negative values for the frequencies")
    end

    dω = median(abs.((circshift(ω_arr,1) - ω_arr)[2:end]))

    t_arr = LinRange(zero(T),2π/dω,2N+2)[1:(2N+1)]
    return t_arr
end

"The inverse of ω_to_t if ω_vec[1] == 0"
function t_to_ω(t_arr::AbstractArray{T}) where T <: AbstractFloat
    N = Int(round((length(t_arr)-one(T))/T(2)))
    maxt = t_arr[2] - t_arr[1] + t_arr[end] - t_arr[1] # subtract t_arr[1] in case t_arr[1] != zero(T)
    maxω = N*2π/maxt
    ω_vec = LinRange(zero(T),maxω,N+1)
    return ω_vec
end

"""
Returns the first element of array which isn't zero (assumes elements are
increasing and distinct)
"""
function firstnonzero(arr::AbstractArray{T}) where T <: AbstractFloat
    if arr[1] != 0
        return arr[1]
    else
        return arr[2]
    end
end

"""
See also: [`DiscreteImpulse`](@ref), [`ContinuousImpulse`](@ref)

Calculates the time response from the frequency response by approximating an
inverse Fourier transform. The time signal is assumed to be real and the
frequenices ω_vec are assumed to be positive (can include zero) and sorted. The
result is convoluted in time ωith the user specified impulse.

We use the Fourier transform convention:
F(ω) =  ∫ f(t)*exp(im*ω*t) dt
f(t) = (2π)^(-1) * ∫ F(ω)*exp(-im*ω*t) dt

To easily sample any time, the default is not FFT, but a discrete version of the transform above.
"""
function frequency_to_time(field_mat::AbstractArray{Complex{T}}, ω_vec::AbstractVector{T},
        t_vec::AbstractArray{T} = ω_to_t(ω_vec);
        impulse::ContinuousImpulse{T} = TimeDiracImpulse(zero(T)),
        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse, t_vec, ω_vec),
        method=:dft) where T <: AbstractFloat

    # In case the used specifies discrete_impulse but not t_vec
    t_vec = discrete_impulse.t

    if size(field_mat,1) != size(ω_vec,1) error("Vector of frequencies ω_vec expected to be same size as size(field_mat,1)") end

    function f(t::T,j::Int)
        fs = discrete_impulse.in_freq.*field_mat[:,j] .* exp.(-(im*t) .* ω_vec)
        if method == :dft && ω_vec[1] == zero(T)
            fs[1] = fs[1]/T(2) # so as to not add ω=0 tωice in the integral of ω over [-inf,inf]
        end
        fs
    end
    inverse_fourier_integral = (t,j) -> numerical_integral(ω_vec, f(t,j), method)
    u = [inverse_fourier_integral(t,j) for t in discrete_impulse.t, j in axes(field_mat,2)]

    return real.(u)/pi # a constant 1/(2pi) appears due to our Fourier convention, but because we only use positive frequencies, and assume a real time signal, this becomes 1/pi.
end

"""
The inverse of the function frequency_to_time (only an exact inverse when using
:dft integration). We use the Fourier transform convention:
F(ω) =  ∫ f(t)*exp(im*ω*t) dt
"""
function time_to_frequency(field_mat::Union{AbstractArray{T},AbstractArray{Complex{T}}}, t_vec::AbstractVector{T},
        ω_vec::AbstractArray{T} = t_to_ω(t_vec);
        impulse::ContinuousImpulse{T} = TimeDiracImpulse(zero(T)),
        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse, t_vec, ω_vec),
        method=:dft) where T <: AbstractFloat

    # In case the used specifies discrete_impulse but not ω_vec
    ω_vec = discrete_impulse.ω

    # to use an impulse below in time we would need to do a discrete convolution, which we decided against.
    f(ω::T, j::Int) = field_mat[:,j] .* exp.((im*ω) .* t_vec)
    fourier_integral = (ω,j) -> numerical_integral(t_vec, f(ω,j), method)
    uhat = [discrete_impulse.in_freq[i]*fourier_integral(discrete_impulse.ω[i],j) for i in eachindex(discrete_impulse.ω), j in axes(field_mat,2)]

    return uhat
end

function numerical_integral(xs::AbstractArray{T}, fs::Union{AbstractArray{T},AbstractArray{Complex{T}}}, method = :dft) where T <: AbstractFloat

    if method == :trapezoidal
        sum(1:(length(xs)-1)) do xi
            (fs[xi] + fs[xi+1])*(xs[xi+1] - xs[xi])/T(2)
        end
    elseif method == :dft
        fs[1]*(xs[2] - xs[1]) +
        sum(2:length(xs)) do xi
            fs[xi]*(xs[xi] - xs[xi-1])
        end
    else
        error("The method $method for numerical integration is not known.")
    end
end


text:
function horizontal_average_is_correct(arch, FT)
    grid = RegularCartesianGrid(size=(16, 16, 16), length=(100, 100, 100))
    model = Model(grid=grid, architecture=arch, float_type=FT)

    T₀(x, y, z) = 20 + 0.01*z
    set!(model; T=T₀)

    T̅ = HorizontalAverage(model.tracers.T; interval=0.5second)
    computed_profile = T̅(model)
    correct_profile = @. 20 + 0.01 * model.grid.zC

    return all(computed_profile[:][2:end-1] .≈ correct_profile)
end

function nan_checker_aborts_simulation(arch, FT)
    grid=RegularCartesianGrid(size=(16, 16, 2), length=(1, 1, 1))
    model = Model(grid=grid, architecture=arch, float_type=FT)

    # It checks for NaNs in w by default.
    nc = NaNChecker(model; frequency=1, fields=Dict(:w => model.velocities.w.data.parent))
    push!(model.diagnostics, nc)

    model.velocities.w[4, 3, 2] = NaN

    time_step!(model, 1, 1);
end

TestModel(::GPU, FT, ν=1.0, Δx=0.5) =
    Model(grid = RegularCartesianGrid(FT; size=(16, 16, 16), length=(16Δx, 16Δx, 16Δx)),
       closure = ConstantIsotropicDiffusivity(FT; ν=ν, κ=ν),
  architecture = GPU(),
    float_type = FT
)

TestModel(::CPU, FT, ν=1.0, Δx=0.5) =
    Model(grid = RegularCartesianGrid(FT; size=(3, 3, 3), length=(3Δx, 3Δx, 3Δx)),
       closure = ConstantIsotropicDiffusivity(FT; ν=ν, κ=ν),
  architecture = CPU(),
    float_type = FT
)

function max_abs_field_diagnostic_is_correct(arch, FT)
    model = TestModel(arch, FT)
    set!(model.velocities.u, rand(size(model.grid)))
    u_max = FieldMaximum(abs, model.velocities.u)
    return u_max(model) == maximum(abs, model.velocities.u.data.parent)
end

function advective_cfl_diagnostic_is_correct(arch, FT)
    model = TestModel(arch, FT)

    Δt = FT(1.3e-6)
    Δx = FT(model.grid.Δx)
    u₀ = FT(1.2)
    CFL_by_hand = Δt * u₀ / Δx

    model.velocities.u.data.parent .= u₀
    cfl = AdvectiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

function diffusive_cfl_diagnostic_is_correct(arch, FT)
    Δt = FT(1.3e-6)
    Δx = FT(0.5)
    ν = FT(1.2)
    CFL_by_hand = Δt * ν / Δx^2

    model = TestModel(arch, FT, ν, Δx)
    cfl = DiffusiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

get_iteration(model) = model.clock.iteration
get_time(model) = model.clock.time

function timeseries_diagnostic_works(arch, FT)
    model = TestModel(arch, FT)
    iter_diag = Timeseries(get_iteration, model; frequency=1)
    push!(model.diagnostics, iter_diag)
    Δt = FT(1e-16)
    time_step!(model, 1, Δt)

    return iter_diag.time[end] == Δt && iter_diag.data[end] == 1
end

function timeseries_diagnostic_tuples(arch, FT)
    model = TestModel(arch, FT)
    timeseries = Timeseries((iters=get_iteration, itertimes=get_time), model; frequency=2)
    model.diagnostics[:timeseries] = timeseries
    Δt = FT(1e-16)
    time_step!(model, 2, Δt)
    return timeseries.iters[end] == 2 && timeseries.itertimes[end] == 2Δt
end

function diagnostics_getindex(arch, FT)
    model = TestModel(arch, FT)
    iter_timeseries = Timeseries(get_iteration, model)
    time_timeseries = Timeseries(get_time, model)
    model.diagnostics[:iters] = iter_timeseries
    model.diagnostics[:times] = time_timeseries
    return model.diagnostics[2] == time_timeseries
end

function diagnostics_setindex(arch, FT)
    model = TestModel(arch, FT)
    iter_timeseries = Timeseries(get_iteration, model)
    time_timeseries = Timeseries(get_time, model)
    max_abs_u_timeseries = Timeseries(FieldMaximum(abs, model.velocities.u), model; frequency=1)

    push!(model.diagnostics, iter_timeseries, time_timeseries)
    model.diagnostics[2] = max_abs_u_timeseries

    return model.diagnostics[:diag2] == max_abs_u_timeseries
end

@testset "Diagnostics" begin
    @info "Testing diagnostics..."

    for arch in archs
        @testset "Horizontal average [$(typeof(arch))]" begin
            @info "  Testing horizontal average [$(typeof(arch))]"
            for FT in float_types
                @test horizontal_average_is_correct(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "NaN Checker [$(typeof(arch))]" begin
            @info "  Testing NaN Checker [$(typeof(arch))]"
            for FT in float_types
                @test_throws ErrorException nan_checker_aborts_simulation(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "Miscellaneous timeseries diagnostics [$(typeof(arch))]" begin
            @info "  Testing miscellaneous timeseries diagnostics [$(typeof(arch))]"
            for FT in float_types
                @test diffusive_cfl_diagnostic_is_correct(arch, FT)
                @test advective_cfl_diagnostic_is_correct(arch, FT)
                @test max_abs_field_diagnostic_is_correct(arch, FT)
                @test timeseries_diagnostic_works(arch, FT)
                @test timeseries_diagnostic_tuples(arch, FT)
                @test diagnostics_getindex(arch, FT)
                @test diagnostics_setindex(arch, FT)
            end
        end
    end
end


text:
# Macro for checking arguments
macro check_args(K, param, cond, desc=string(cond))
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(K)), ": ", $(string(param)), " = ", $(esc(param)), " does not ",
                "satisfy the constraint ", $(string(desc)), ".")))
        end
    end
end

# Take highest Float among possibilities
# function promote_float(Tₖ::DataType...)
#     if length(Tₖ) == 0
#         return Float64
#     end
#     T = promote_type(Tₖ...)
#     return T <: Real ? T : Float64
# end

check_dims(K,X,Y,featdim,obsdim) = check_dims(X,Y,featdim,obsdim) && (size(K) == (size(X,obsdim),size(Y,obsdim)))

check_dims(X,Y,featdim,obsdim) = size(X,featdim) == size(Y,featdim)

feature_dim(obsdim::Int) = obsdim == 1 ? 2 : 1

base_kernel(k::Kernel) = eval(nameof(typeof(k)))

base_transform(t::Transform) = eval(nameof(typeof(t)))

"""
Will be implemented at some point
```julia
    params(k::Kernel)
    params(t::Transform)
```
For a kernel return a tuple with parameters of the transform followed by the specific parameters of the kernel
For a transform return its parameters, for a `ChainTransform` return a vector of `params(t)`.
"""
#params


text:
@testset "Mutations" begin
    lx = [0.0, 0.0]
    ux = [2.0, 1.0]
    @test PM(lx, ux, Inf)([1.0, 2.0])[1] ∈ [0.0, 2.0]
    @test PM(lx, ux, Inf)([1.0, 2.0])[2] == 1.0
    @test MIPM(lx, ux, Inf, 1.0)(Real[1.0, 2])[1] ∈ [0.0, 2.0]
    @test MIPM(lx, ux, 10.0, Inf)(Real[1.0, 2])[2] == 1.0
end


text:
import Base.getindex

export LinearInterpolation,
	getindex,
	getindex_symbolic

abstract LinearInterpolation <: SimpleInterpolation

@inline function getindex{T<:Real}(interp::Type{LinearInterpolation}, x::T)
	ix = floor(Int, x)
	rx = x-ix
	return ((1-rx, ix,), (rx, ix+1,))
end

function getindex_symbolic(interp::Type{LinearInterpolation}, x::Union{Symbol, Expr})
	@gensym a b
	return :($a = floor($x); $b = floor(Int,$x)), ((:($a+1-$x), :($b)), (:($x-$a), :($b+1)))
end

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
function BackgroundVelocityFields(bg, grid, clock)
    u = :u ∈ keys(bg) ? regularize_background_field(Face, Center, Center, bg[:u], grid, clock) : ZeroField()
    v = :v ∈ keys(bg) ? regularize_background_field(Center, Face, Center, bg[:v], grid, clock) : ZeroField()
    w = :w ∈ keys(bg) ? regularize_background_field(Center, Center, Face, bg[:w], grid, clock) : ZeroField()

    return (u=u, v=v, w=w)
end

function BackgroundTracerFields(bg, tracer_names, grid, clock)
    tracer_fields =
        Tuple(c ∈ keys(bg) ?
              regularize_background_field(Center, Center, Center, getindex(bg, c), grid, clock) :
              ZeroField()
              for c in tracer_names)

    return NamedTuple{tracer_names}(tracer_fields)
end

#####
##### Convenience for model constructor
#####

function BackgroundFields(background_fields, tracer_names, grid, clock)
    velocities = BackgroundVelocityFields(background_fields, grid, clock)
    tracers = BackgroundTracerFields(background_fields, tracer_names, grid, clock)
    return (velocities=velocities, tracers=tracers)
end

"""
    BackgroundField{F, P}

Temporary container for storing information about BackgroundFields.
"""
struct BackgroundField{F, P}
    func:: F
    parameters :: P
end

"""
    BackgroundField(func; parameters=nothing)

Returns a `BackgroundField` to be passed to `IncompressibleModel` for use
as a background velocity or tracer field.

If `parameters` is not provided, `func` must be callable with the signature

```julia
func(x, y, z, t)
```

If `parameters` is provided, `func` must be callable with the signature

```julia
func(x, y, z, t, parameters)
```
"""
BackgroundField(func; parameters=nothing) = BackgroundField(func, parameters)

regularize_background_field(X, Y, Z, f::BackgroundField{<:Function}, grid, clock) =
    FunctionField{X, Y, Z}(f.func, grid; clock=clock, parameters=f.parameters)

regularize_background_field(X, Y, Z, func::Function, grid, clock) =
    FunctionField{X, Y, Z}(func, grid; clock=clock)

function regularize_background_field(X, Y, Z, field::AbstractField, grid, clock)
    if location(field) != (X, Y, Z)
        throw(ArgumentError("Cannot use field at $(location(field)) as a background field at $((X, Y, Z))"))
    end

    return field
end

Base.show(io::IO, field::BackgroundField{F, P}) where {F, P} =
    print(io, "BackgroundField{$F, $P}", '\n',
          "├── func: $(short_show(field.func))", '\n',
          "└── parameters: $(field.parameters)")


text:
test("BWS statistic") do

  test("a small hand-calculated BWS case with 2 elements in each sample") do

    x = [1.0, 2.0]
    n = 2
    y = [3.0, 4.0]
    m = 2

    r = [1, 2] # ranks for values in x
    s = [3, 4] # ranks for values in y

    mn = m+n
    mp1 = m+1
    np1 = n+1

    b_x_part_for_i_1 = (r[1] - 1*mn/n)^2 / ( (1/np1) * (1 - 1/np1) * (m * mn / n) )
    b_x_part_for_i_2 = (r[2] - 2*mn/n)^2 / ( (2/np1) * (1 - 2/np1) * (m * mn / n) )
    b_x = (1 / n) * ( b_x_part_for_i_1 + b_x_part_for_i_2 )

    b_y_part_for_i_1 = (s[1] - 1*mn/m)^2 / ( (1/mp1) * (1 - 1/mp1) * (n * mn / m) )
    b_y_part_for_i_2 = (s[2] - 2*mn/m)^2 / ( (2/mp1) * (1 - 2/mp1) * (n * mn / m) )
    b_y = (1 / m) * ( b_y_part_for_i_1 + b_y_part_for_i_2 )

    b = (b_x + b_y) / 2

    @t FeldtLib.bws_statistic(x, y) == b

    bwst = FeldtLib.BWSTest(x, y)
    @t bwst.b == FeldtLib.bws_statistic(x, y)
    @t bwst.n == n
    @t bwst.m == m

  end

  test("a small hand-calculated BWS case with 3+2 elements") do

    x = [1.0, 2.0, 5.0]
    n = 3
    y = [3.0, 4.0]
    m = 2

    r = [1, 2, 5] # ranks for values in x
    s = [3, 4] # ranks for values in y

    mn = m+n
    mp1 = m+1
    np1 = n+1

    bx1 = (r[1] - 1*mn/n)^2 / ( (1/np1) * (1 - 1/np1) * (m * mn / n) )
    bx2 = (r[2] - 2*mn/n)^2 / ( (2/np1) * (1 - 2/np1) * (m * mn / n) )
    bx3 = (r[3] - 3*mn/n)^2 / ( (3/np1) * (1 - 3/np1) * (m * mn / n) )
    bx = (1 / n) * ( bx1 + bx2 + bx3)

    by1 = (s[1] - 1*mn/m)^2 / ( (1/mp1) * (1 - 1/mp1) * (n * mn / m) )
    by2 = (s[2] - 2*mn/m)^2 / ( (2/mp1) * (1 - 2/mp1) * (n * mn / m) )
    by = (1 / m) * ( by1 + by2 )

    b = (bx + by) / 2

    @t FeldtLib.bws_statistic(x, y) == b

    bwst = FeldtLib.BWSTest(x, y)
    @t bwst.b == FeldtLib.bws_statistic(x, y)
    @t bwst.n == n
    @t bwst.m == m

  end

  test("example from Neuhauser 2004 paper, in table 3") do

    control_group_ranks = [1, 2, 3, 4, 6, 7, 8]
    experimental_group_ranks = [5, 9, 10, 11, 12, 13, 14]
    rs = vcat(control_group_ranks, experimental_group_ranks)
    n = length(control_group_ranks)
    m = length(experimental_group_ranks)
    @t isapprox(FeldtLib.bws_statistic_from_ranks(rs, n, m), 5.132; atol = 0.001)

    b, pvalue, m, sd = r = FeldtLib.bws_test_sampled_from_ranks(rs, n, m, 10000)

    # The exact p-value calculated by Neuhauser was 0.0029 so we should be close
    # to that.
    @t isapprox(pvalue, 0.0029; atol = 0.001)

  end

end


text:
"""
write_xyz(data::Data, name_file::AbstractString)
Do this will write all information in `data` to file `name_file` in `.xyz` formation

# Example
```julia-repl
data_cell = genr_cell([10 10 3])
str = TiO2()
data_atom = genr_atom(data_cell, str)
write_xyz(data_atom, "test.xyz")
```
"""
function write_xyz(data::Data, name_file::AbstractString)
io = open(name_file, "w")
num_atoms = data.data_basic.num_atoms
name_atoms = data.data_str.atom_name
@printf(io, "%.d\n", num_atoms)
@printf(io, ".XYZ file of %s and %s created by lmp_str \n", name_atoms[1], name_atoms[2])
for atom = 1:num_atoms
    @printf(io, "%s", name_atoms[data.vec_atom[atom].typ])
    for dim = 1:3
        @printf(io, "\t%4.4f", data.vec_atom[atom].coord[dim])
    end
    @printf(io, "\n")
end
close(io)
end


text:
module SatelliteData
"""
Module for the Download and processing of descriptors of satellitar data from the Copernicus project databases.
"""

using ArchGDAL
using CombinedParsers
using CombinedParsers.Regexp
using CSV
using Dates
using DataFrames
using Downloads
using HTTP

export getProductsDF, selectProducts, downloadProducts

const agd = ArchGDAL

# Syntax to parse the XMl reppresentation of a product
@syntax product = Repeat(
    Either(
        Sequence(
            "<",
            :type    => re"[^< >]+ ",
            :opening => re"name=\"[^<\">]+\">",
            :content => re"[^<>]+",
            :closing => re"</[^<>]+>"
        ),
        Sequence(
            re"<[^<>]+>",
            re"[^<>]+",
            re"</[^<>]+>"
        ),
        re"<[^<>]+>"
    )
)

# Funzione presa da: https://stackoverflow.com/questions/48104390/julias-most-efficient-way-to-choose-longest-array-in-array-of-arrays
"""
    maxLenIndex( vect::AbstractVector )

Retunr the index of the longest vector inside `vect` ad its length.
"""
function maxLenIndex( vect::AbstractVector )
    len = 0
    index = 0
    @inbounds for i in 1:length(vect)
        l = length( vect[i] )
        l > len && ( index = i; len=l )
    end
    return index
end

"""
    authenticate( username::AbstractString, password::AbstractString[, type::AbstractString ] )

Create an authentication token of type "type" for the user
"""
function authenticate( username::AbstractString, password::AbstractString, type::AbstractString = "Base" )
    if type == "Base"
        return HTTP.Base64.base64encode("$username:$password") #Trasformazione dei dati nel formato valido per la verifica
    end
end

"""
    parseConvert( xmlType::AbstractString, value::AbstractString )

Given the xml type assigned to "value", convert the latter to its correct Julia type
"""
function parseConvert( xmlType::AbstractString, value::AbstractString )
    if xmlType == "str"
        if !isnothing( match( re"^[0-9]+[.,][0-9]+$", value ) )
            return tryparse( Float64, value )
        elseif !isnothing( match( re"^[0-9]+$", value ) )
            return tryparse( Int64, value )
        else
            return value
        end
    end
    if xmlType == "int"
        return tryparse( Int64, value )
    end
    if xmlType == "double"
        return tryparse( Float64, value )
    end
    if xmlType == "date"
        return DateTime( value[1:end-1], "y-m-dTH:M:S.s" )
    end
    throw( DomainError( (xmlType, value), "Undefined type of $value" )  )
end

"""
    getAoi( path::AbstractString )::AbstractString

Return the WKT `String` rappresentation of the geometry contained in the shapefile `path`
"""
function getAoi( path::AbstractString )::AbstractString
    data = agd.read(path)
    layer = agd.getlayer(data, 0)
    geometry = agd.getgeom(collect(layer)[1])
    src_crs = agd.getspatialref(geometry)
    trg_crs = agd.importEPSG(4326)
    agd.createcoordtrans(src_crs, trg_crs) do transform
        agd.transform!(geometry, transform)
    end
    geom = agd.toWKT(geometry)
    return geom[1:7] * replace( geom[9:end], " " => "%20" )
end

#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#                                                            PRODUCTS XML FILES DOWNLOAD AND PROCESSING
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

"""
    getProductsBuffer( authToken::AbstractString[, start::Integer, maxNumber::Union{Integer, Nothing} ] )

Obtain `maxNumber` XML description of products, starting from `start`, through `authToken`, returning the IOBuffer that contains them all
"""
function getProductsBuffer( username::AbstractString, password::AbstractString, aoi::AbstractString; numMonths::Integer=6, start::Integer=0,
                            max::Union{Integer, Nothing}=nothing, last::Bool=false )
 # Definition of the components of The URL
    authToken = authenticate(username, password)
    query2 = "[NOW-$(numMonths)MONTHS%20TO%20NOW]%20AND%20footprint:\"Intersects($aoi)\""
    query = "search?start=0&rows=0&q=ingestiondate:$query2"
    iob = IOBuffer()

 # Get number of total products
    #Download the firs page (This way we can get the number we are interested in and also the first page that we would download anyway)
    Downloads.download( "https://scihub.copernicus.eu/dhus/$query", iob, headers = [ "Authorization" => "Basic $authToken" ] )
    lines = String(take!(iob))
    val = tryparse( Int64, split( split( lines, "totalResults>" )[2], "<" )[1] )
    # "count" has two values: the first is the number of pages (undreds of products) and the second is a number of products smaller than 100 to be downloaded in an additional page
    count = val

    #Check if `start` has a sensible value
    if ( last && start == 0 ) || start > val
        start = val -1
    end

    # Check if `maxNumber` has a sensible value
    if !isnothing(max) && max > 0 && max < val
        count = max
    end

    # Calculate the value of count in terms of pages of 100 products each accounting for the starting point
    if start < count
        count -= start
    end
    count = [ count ÷ 100, count % 100 ]

 # Download the desired number pages and store the in an IOBuffer
    if count[1] > 0
        for i in 0:count[1]-1
            query = "search?start=$(start + ( i * 100 ))&rows=100&q=ingestiondate:$query2"
            Downloads.download( "https://scihub.copernicus.eu/dhus/$query", iob, headers = [ "Authorization" => "Basic $authToken" ] )
        end
    end
    if count[2] > 0
        query = "search?start=$(start + (count[1] * 100))&rows=$(count[2])&q=ingestiondate:$query2"
        Downloads.download( "https://scihub.copernicus.eu/dhus/$query", iob, headers = [ "Authorization" => "Basic $authToken" ] )
    end
    return iob
end
# io = getProductsBuffer( "davidefavaro", "Tirocinio", "POLYGON((9.5000%2047.0000,%2014.0000%2047.0000,%2014.0000%2044.0000,%209.5000%2044.0000,%209.5000%2047.0000))",  max=300 )

"""
    getProductsDicts( fileIO::IO )

Given the IOBuffer obtained through "getProductsBuffer()", return an array of the dictionaries containing the informations on each of the products of the page
"""
function getProductsDicts( fileIO::IO )
    # Get the downloaded XML representations of the products
    original = replace( String( take!(fileIO) ), "\n" => "" )

    # Split the downloaded files into pages
    pages = split(original, r"<\?xml [^<>]+><[^<>]+>", keepempty=false)

    # Split the result in a vector of ready-to-be-parsed strings representing single products
    vector = reduce( vcat, [ split( page, r"</?entry>", keepempty=false )[2:end-1] for page in pages ] )

    # Parse the strings
    products = [ product(x)[8:end] for x in vector ]

    # Generate a vector of dictionaries containing the details for each product of the original page adding to each of them an additional value to account for the
        # original order of the data
    return [ push!(
                Dict( Symbol( join( prod[:opening][7] ) ) => parseConvert( join( prod[:type][1] ), join( prod[:content] ) ) for prod in products[i] ),
                :orderNumber => i,
            )
            for i in 1:length(products) ]
end

#   dict = getProductsDicts(io)

#   io = getProductsBuffer( authenticate("davidefavaro","Tirocinio"), 5115, 50 )
#   res = getProductsDicts(io)

"""
    getProductsDF( authToken::AbstractString[, start::Integer, maxNumber::Union{Integer, Nothing} ] )

Generate the DataFrame containing the data of `maxNumber` products using `username` and `password` for authentication, if `maxNumber` is not specified, the df will include all available products,
if `start` is specified the downloaded data will begin from that index
"""
function getProductsDF( userame::AbstractString, password::AbstractString, aoi::AbstractString; numMonths::Integer=6, start::Integer=0, max::Union{Integer, Nothing}=nothing, last::Bool=false )::DataFrame
 # Download "maxNumber" pages and return the buffer containing them
    io = getProductsBuffer(userame, password, aoi, numMonths=numMonths, start=start, max=max, last=last)
 # Create a vector of dictionaries of the products
    dict_vect = getProductsDicts(io)
 # Obtain the existing subsets of attributes of the products
    keys_groups = unique( keys.(dict_vect) )
 # Divide the dictionaries in groups homogeneus on their attributes
    grouped_vect = [ filter( x -> keys(x) == ks, dict_vect ) for ks in keys_groups ]
 # Turn each group in a DataFrame
    dfs_vect = DataFrame.(grouped_vect)
 # Merge all Dataframes using "DataFrames.append" to create a Dataframe with the union of all the possible columns and the right "missing" values
    data = dfs_vect[1]
    for df in dfs_vect[2:end]
        append!(data, df, cols=:union)
    end
 # Convert `footprint` e `gmlfootprint` columns in geometries
    data[!, :footprint] = agd.fromWKT.( data[:, :footprint] )
    data[!, :gmlfootprint] = agd.fromGML.( replace.( replace.( data[:, :gmlfootprint], "<" => "<" ), ">" => ">" ) )
 # Order the rows based on `orderNumber`, then remove said column
    sort!( data, [:orderNumber] )
    select!( data, Not(:orderNumber) )
 # Create the column indicating wether a product has been downloaded
    insertcols!( data, ( :available => fill( true, nrow(data) ) ) ),( :downloaded => fill( false, nrow(data) ) )
    return data
end

"""
    saveProductsDF( targetDirectory::AbstractString, data::DataFrame; overwrite::Bool=false )

Save "data" in "targetDirectory" if not already existing or if "overwrite" is true, otherwise append its content to "data.csv"
"""
function saveProductsDF( targetDirectory::AbstractString, data::DataFrame; overwrite::Bool=false )
    condition = !overwrite && in("data.csv", readdir(targetDirectory))
    if condition
        # Preexisting data
        old_data = CSV.read( targetDirectory*"\\data.csv", DataFrame )

        # Index of the first entry of the preexisting data
            # It will be used to find the duplicates in "data"
        old_first_id = old_data[1, :uuid]

        # Index of the first duplicate in "data"
        first_dup_idx = findfirst( ==(old_first_id), data[:, :uuid] )

        # Remove all the duplicates from "data"
        filter!( prod -> rownumber(prod) > first_dup_idx, data )
    end
    CSV.write( targetDirectory*"\\data.csv", data, append = condition )
end

"""
    createUnitsDF( columns::AbstractVector{AbstractString}, units::AbstractVector{AbstractString} )

Given the array of column names and the array of their respective units of measure, create a df associating both
"""
function createUnitsDF( columns::AbstractVector{AbstractString}, units::AbstractVector{AbstractString} )
    return Dataframe( Dict( columns .=> units ) )
end

"""
    setDownloaded( fileIDs::Union{ AbstractString, AbstractVector{AbstractString} } )

Mark all products in "fileIDs" as already downloaded
"""
function setDownloaded( fileIDs::Union{ AbstractString, AbstractVector{AbstractString} } )
    updateProductsVal( fileIDs, "products", "downloaded", true )
end

"""
    setUnavailable( fileID::Union{ AbstractString, AbstractVector{AbstractString} } )

Mark all products in "fileIDs" as unavailable
"""
function setUnavailable( fileIDs::Union{ AbstractString, AbstractVector{AbstractString} } )
    updateProductsVal( fileIDs, "products", "available", false )
end

"""
    checkAvailabile( data_path::AbstractString, username::AbstractString, password::AbstractString )

Check and update availability of all the products entries contained in the CSV file `data_path` (the `username` and `password` are necessary for the check and update of the values).
"""
function checkAvailabile( data_path::AbstractString, username::AbstractString, password::AbstractString )
    # Get the last available product as a DataFrame row
    data = getProductsDF( authenticate(username, password), max=1, last=true  )
    # Load the preexisting data
    old_data = CSV.read( path, DataFrame )
    # Obtain the `uuid` of the last available product
    last_available_id = data[end, :uuid]
    # Obtain the index of the last available product of "old_data"
    last_available_idx = findfirst( ==(last_available_id), old_data[:, :uuid] )
    # Set all the products not of level 0, from the last available to the end, as unavailable
    for i in last_available_idx:size(old_data)[1]
        if ismissing(old_data[i, :productlevel]) || old_data[i, :productlevel] != "L0"
            old_data[i, :available] = false
        end
    end
    CSV.write( *( @__DIR__, "\\Dati di prova\\data.csv" ), old_data )
end

#=
test = [ "C:\\Users\\DAVIDE-FAVARO\\Desktop\\XML\\1.xml",
         "C:\\Users\\DAVIDE-FAVARO\\Desktop\\XML\\2.xml",
         "C:\\Users\\Lenovo\\Desktop\\XML\\1.xml",
         "C:\\Users\\Lenovo\\Desktop\\XML\\2.xml",
         "C:\\Users\\Lenovo\\Desktop\\XML\\Prod_Test.xml" ]

out = [ "D:\\Vario\\Stage",
        "C:\\Users\\Lenovo\\Desktop\\XML",
        "C:\\Users\\Lenovo\\Desktop\\XML\\Test",
        "C:\\Users\\DAVIDE-FAVARO\\Desktop",
        "C:\\Users\\DAVIDE-FAVARO\\Desktop\\XML" ]
=#

#   df = getProductsDF( authenticate("davidefavaro","Tirocinio"), 1000 )
#   saveProductsDF( out[2], df )

#   df = getProductsDF( authenticate("davidefavaro","Tirocinio"), max=1, last=true )
#   odf = CSV.read( "D:\\Documents and Settings\\DAVIDE-FAVARO\\My Documents\\GitHub\\Tirocinio\\Dati di prova\\data.csv", DataFrame )
#   odf = CSV.read( "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\Dati di prova\\data.csv", DataFrame )

#   saveProductsDF( "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\Dati di prova", df )
#   saveProductsDF( "D:\\Documents and Settings\\DAVIDE-FAVARO\\My Documents\\GitHub\\Tirocinio\\Dati di prova", df )

#   odf = CSV.read( "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\Dati di prova\\data.csv", DataFrame )
#   checkAvailabile(*( @__DIR__, "\\Dati di prova\\data.csv"), "davidefavaro", "Tirocinio")
#   ndf = CSV.read( "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\Dati di prova\\data.csv", DataFrame )

# :footprint o :gmlfootprint ci dovrebbero dare il poligono
#= La selezione dei prodotti potrebbe avvenire inbase a:
    - :processinglevel / :productlevel ( livello 2A )
    - :instrumentname / :instrumentshortname
    - :size
    - :cloudcoverpercentage + :platformname ( cc minima per Sentinel-2A )
    - :status ?
    - :available
=#

#   sat_file = *( @__DIR__, "\\..\\Mappe\\sat\\sette_sorelle.shp" )
#   aoi = getAoi(sat_file)
#   df = getProductsDF( authenticate("davidefavaro","Tirocinio"), aoi=aoi )
#
# PRENDERE IL PRIMO PRODOTTO PER OGNI MESE POTREBBE NON ESSERE UNA BUON IDEA
#   idxs = [ findfirst( date -> month(date) == m, df[:, :beginposition] ) for m in month( now() - Month(6) ) : month( now() ) ]
#
#   res = df[ idxs, : ]
#
#   aoi_geom = agd.fromWKT( replace(aoi, "%20" => " ") )
#
#   intersections = agd.intersection.()

"""
    selectProducts( userame::AbstractString, password::AbstractString, aoi_path::AbstractString, num_per_month::Integer, from::Integer, to::Integer )

Return products that overlap the geometry specified by `aoi_path` shapefile, using `username` and `password` for authentication return `num_per_month` products for each month
in the `from`-`to` interval.
"""
function selectProducts( userame::AbstractString, password::AbstractString, aoi_path::AbstractString, num_per_month::Integer, from::Integer=0, to::Integer=6 )
    aoi = getAoi(aoi_path)
 # DA AGGIUNGERE LA POSSIBILITA' DI SPECIFICARE IL MESE DA CUI INIZIARE A PRENDERE I PRODOTTI
    #   df::DataFrame = getProductsDF(userame, password, aoi, numMonths=(Months(to) - Months(from)).value)
    df::DataFrame = getProductsDF(userame, password, aoi, numMonths=6 )
    idxs = Vector{Int64}()
    condition( date, platform, clouds, level, m, sat ) = month(date) == m && platform == sat && (platform != "Sentinel-2" || ( !ismissing(clouds) && clouds < 30.0 )) && !ismissing(level) && level == "L2"
    #   for m in month.(collect( Month(from):Month(1):Month(to) ))
    for m in month.(collect( now()-Month(6):Month(1):now() ))
        ind::Int64 = 1
        for i in 1:num_per_month
            for sat in ["Sentinel-1", "Sentinel-2", "Sentinel-3"]
                first::Union{Int64, Nothing} = findfirst( row -> condition(row..., m, sat), eachrow( df[ ind:end, [:beginposition, :platformname, :cloudcoverpercentage, :productlevel] ] ) )
                if !isnothing(first)
                    ind = first
                    push!(idxs, first)
                    ind += 1
                end
            end
        end
    end
    return df[idxs, :]
end

#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#                                                                           PRODUCTS DOWNLOAD
#||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

"""
    dowloadProducts( userame::AbstractString, password::AbstractString, uuid::Int64, output_dir_path::AbstractString )

Download the product identified by `uuid` using `username` and `password` for authentication and save it in `output_dir_path` using `uuid` as filename.
If `show_progress` is set to `true` the progress of the download will be printed, if `verbose` is set to `true` additional informations will be printed.
"""
dowloadProducts( userame::AbstractString, password::AbstractString, uuid::Int64, output_dir_path::AbstractString, show_progress::Bool=false, verbose::Bool=false ) =
    Downloads.download(
        "https://scihub.copernicus.eu/dhus/odata/v1/Products('$uuid')/\$value",
        output_dir_path*"\\$uuid",
        headers = [ "Authorization" => "Basic $(authenticate(username, password))" ],
        progress = show_progress ? ( ( total, now ) -> println("$(now/total*100)% ( $now / $total )") ) : nothing,
        verbose = verbose
    )

"""
    dowloadProducts( userame::AbstractString, password::AbstractString, uuids::Vector{Int64}, output_dir_path::AbstractString )

Download all the products identified by the ids contained in `uuids` using `username` and `password` for authentication and save them in `output_dir_path` using their uuids as filenames.
If `show_progress` is set to `true` the progress of the download will be printed, if `verbose` is set to `true` additional informations will be printed.
"""
function dowloadProducts( userame::AbstractString, password::AbstractString, uuids::Vector{Int64}, output_dir_path::AbstractString, show_progress::Bool=false, verbose::Bool=false )
    authtoken = authenticate(username, password)
    if show_progress
        len = length(uuids)
        for (i, id) in enumerate(uuids)
            println("Product $uuid ($i of $len)")
            Downloads.download(
                "https://scihub.copernicus.eu/dhus/odata/v1/Products('$id')/\$value",
                output_dir_path*"\\$id",
                headers=["Authorization" => "Basic $authtoken"],
                progress = ( total, now ) -> println("\t$(now/total*100)% ( $now / $total )"),
                verbose = verbose
            )
        end
    else
        Downloads.download(
                "https://scihub.copernicus.eu/dhus/odata/v1/Products('$id')/\$value",
                output_dir_path*"\\$id",
                headers=["Authorization" => "Basic $authtoken"],
                verbose = verbose
            )
    end
end

#   downloadProducts( "b57f225e-d288-4e4a-bb35-2a7eb75d60e4", out[3], authenticate("davidefavaro", "Tirocinio")  )

# SI OTTENGONO SOLO Sentinel-3 PERCHE' TUTTI GLI ALTRI PRODOTTI HANNO :productlevel missing
# PER MESE 6 SI HANNO SOLO Sentinel-3
#   n, p = ("davidefavaro", "Tirocinio")
#   dir = "C:\\Users\\DAVIDE-FAVARO\\Desktop"
#   dir = "D:\\Z_Tirocinio_Dati\\Copernicus Data"
#   uuid = "de9494a3-fdd9-46ce-98cb-42a2632e8b87"
#   sat_file = *( @__DIR__, "\\..\\Mappe\\sat\\sette_sorelle.shp" )
#   res = selectProducts(n, p, sat_file, 1, 12, 6)
#=
    downloadProducts(n, p, res[:, :uuid], dir)
    downloadProducts(n, p, uuid, dir)
=#

# @code_warntype res = selectProducts( authToken, sat_file, 1, 12, 6 )

end #module

max_stars_repo_path: src/Data Gathering/SatelliteData.jl

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_path: src/Data Gathering/SatelliteData.jl

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_path: src/Data Gathering/SatelliteData.jl

max_forks_repo_name: DavideFavaro/Tirocinio


text:
mutable struct ConvNN
    reg_coeff
    learning_rate
    transform
    activation

    y    # observation
    x0   # true state

    biases
    kernels
    x    # estimate

    train_step
    sess

    ConvNN() = new()
end

"""
    ConvNN(nLayers::Integer, nAntennas::Integer; transform = (x,_) -> x, learning_rate=1e-3,reg_coeff=1e-6,activation=nn.relu)

Input and output of the NN are complex-valued vectors of dimension `nCoherence, nFilterLength`
where `nFilterLength` is calculated from `nAntennas` and `transform`.
For mini-batch training, these are stacked into three dimensional arrays of dimension
`nBatches, nCoherence, nFilterLength`.
This method performs a random initialization of all kernels and biases.
"""
function ConvNN(nLayers::Integer, nAntennas::Integer; transform = (x,_) -> x, learning_rate=1e-3,reg_coeff=1e-6,activation=nn.relu)

    nFilterLength = size(transform(eye(nAntennas), :notransp),1)

    kernels = [rand(TruncatedNormal(0,0.1,-0.3,0.3), nFilterLength) for i in 1:nLayers]
    biases  = [0.1*ones(nFilterLength) for i in 1:nLayers]

    return ConvNN(kernels, biases, transform = transform, learning_rate=learning_rate, reg_coeff=reg_coeff, activation=activation)
end

function ConvNN(kernels::Array, biases::Array; transform = (x,_) -> x, learning_rate=1e-3,reg_coeff=1e-6,activation=nn.relu)
    nLayers = length(kernels)
    kernels_tf = Array{TensorFlow.Variables.Variable}(nLayers)
    biases_tf  = Array{TensorFlow.Variables.Variable}(nLayers)
    for i in 1:nLayers
        kernels_tf[i] = Variable(Float32.(kernels[i]))
        biases_tf[i]  = Variable(Float32.(biases[i]))
    end
    return ConvNN(kernels_tf, biases_tf; transform = transform, learning_rate=learning_rate,reg_coeff=reg_coeff,activation=activation)

end
# This method takes set values as input for kernels and biases
function ConvNN(kernels::Array{TensorFlow.Variables.Variable}, biases::Array{TensorFlow.Variables.Variable}; transform = (x,_) -> x, learning_rate=1e-3,reg_coeff=1e-6,activation=nn.relu)
    est = ConvNN()
    const nBatches = -1 # special value in TensorFlow (nBatches not known at this point and may also vary during training etc.)

    est.reg_coeff     = reg_coeff
    est.learning_rate = learning_rate
    est.transform     = transform
    est.activation    = activation

    est.y             = placeholder(Complex64) # dims: nBatches, nCoherence, nFilterLength
    est.x0            = placeholder(Complex64) # dims: nBatches, nCoherence, nFilterLength

    est.kernels = kernels
    est.biases  = biases
    nLayers = length(kernels)
    nFilterLength = get(get_shape(kernels[1]).dims[1])

    # Build NN
    intermediates = Array{TensorFlow.Tensor}(nLayers)

    # For each batch b=1,...,nBatches, the input to the network is given as sum_t |y[b,t,:]|^2
    intermediates[1] = reduce_mean(square(real(est.y)) + square(imag(est.y)), axis=[2]) # dims: nBatches, nFilterLength
    for i in 1:nLayers-1
        intermediates[i+1] = activation( circ_conv(intermediates[i], kernels[i]) + biases[i] ) # dims: nBatches, nFilterLength
    end
    # Last layer is without activation function
    filt = circ_conv(intermediates[end], kernels[end]) + biases[end] # dims: nBatches, nFilterLength

    # The output of the NN is an element-wise multiplication of the filter with the input data
    est.x = reshape(filt, [nBatches, 1, nFilterLength]) .* est.y # dims: nBatches, nCoherence, nFilterLength

    # Set learning parameters
    regu = reduce_mean(square(kernels[1])) # regularization term
    for i in 2:nLayers
        regu = regu + reduce_mean(square(kernels[i]))
    end
    cost_function  = reduce_mean(square(real(est.x-est.x0)) + square(imag(est.x-est.x0))) + (est.reg_coeff*regu)
    opt_alg        = tf.train.AdamOptimizer(est.learning_rate)
    est.train_step = tf.train.minimize(opt_alg,cost_function)

    est.sess = Session()
    init = global_variables_initializer()
    run(est.sess,init)

    # close TF-session when network is de-referenced
    # note: All nodes remain in the graph. To reset the graph, restart julia.
    finalizer(est, est -> close(est.sess))
    est
end

# This method takes a (smaller) network as input and interpolates
# the values of the kernels and biases to the (larger) new filter length.
# The TF session of the smaller network is then closed.
function resize( est::ConvNN, nAntennas::Integer; learning_rate=1e-3,reg_coeff=1e-6 )

    nLayers          = length(est.kernels)
    nFilterLength    = size(est.transform(eye(nAntennas), :notransp),1)
    nFilterLengthOld = get(get_shape(est.kernels[1]).dims[1])
    kernels = Array{TensorFlow.Variables.Variable}(nLayers)
    biases  = Array{TensorFlow.Variables.Variable}(nLayers)
    for i in 1:nLayers
        kernel_old       = run(est.sess, est.kernels[i])
        bias_old         = run(est.sess, est.biases[i])

        kernel_old_itp = interpolate(kernel_old, BSpline(Quadratic(Line())), OnGrid())
        bias_old_itp   = interpolate(bias_old,   BSpline(Quadratic(Line())), OnGrid())

        # normalize kernel such that its energy remains the same
        kernel = kernel_old_itp[linspace(1,nFilterLengthOld,nFilterLength)] .* nFilterLengthOld/nFilterLength
        bias   = bias_old_itp[  linspace(1,nFilterLengthOld,nFilterLength)]

        kernels[i] = Variable(Float32.(kernel[:])) # dims: nFilterLength
        biases[i]  = Variable(Float32.(bias[:]))   # dims: nFilterLength
    end

    # create new network and overwrite reference to old network
    ConvNN(kernels, biases, transform = est.transform, learning_rate=learning_rate, reg_coeff=reg_coeff, activation=est.activation)
end

function train!(est::ConvNN, y, x0)
    run(est.sess, est.train_step, Dict([(est.y, cnn_perm(est.transform(y,:notransp))), (est.x0, cnn_perm(est.transform(x0,:notransp)))]))
end
function estimate(est::ConvNN, y)
    x = est.transform( cnn_perm(run(est.sess, est.x, Dict([(est.y, cnn_perm(est.transform(y,:notransp)))]))), :transp)
end

#
# The TensorFlow.Saver functions are somewhat difficult to use with Julia.
# Julia can handle only a single graph, hence all networks are defined in
# the same graph. Thus, loading and saving always concerns all networks
# (if multiple networks are defined)
#
# function save(est::ConvNN, filename)
#     saver = tf.train.Saver()
#     tf.train.save(saver, est.sess, filename)
# end
# function load(est::ConvNN, filename)
#     if isfile(filename)
#         saver = tf.train.Saver()
#         tf.train.restore(saver, est.sess, filename)
#     else
#         warn(@sprintf "Ignoring load... File %s not found" filename)
#     end
# end

function save(est::ConvNN, filename)
    nLayers = length(est.kernels)
    nFilterLength = get(get_shape(est.kernels[1]).dims[1])
    kernels = zeros(nFilterLength,nLayers)
    biases  = zeros(nFilterLength,nLayers)
    for i in 1:nLayers
        kernels[:,i] = run(est.sess, est.kernels[i])
        biases[:,i]  = run(est.sess, est.biases[i])
    end
    data = Dict("kernels" => kernels, "biases" => biases)

    JLD.save(filename, "kernels", kernels, "biases", biases)
end

function load!(est::ConvNN, filename)
    if isfile(filename)
        @printf "Loading %s\n" filename
        kernels_vals = load(filename, "kernels")
        biases_vals  = load(filename, "biases")

        nLayers = length(est.kernels)
        for i in 1:nLayers
            ass = assign(est.kernels[i], Float32.(kernels_vals[:,i]))
            run(est.sess, ass)
            ass = assign(est.biases[i], Float32.(biases_vals[:,i]))
            run(est.sess, ass)
        end
    else
        warn(@sprintf "Ignoring load... File %s not found" filename)
    end
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
function _set_scaling!(con_scale::AbstractVector, jac::SparseMatrixCOO)
    @simd for i in 1:nnz(jac)
        row = @inbounds jac.I[i]
        @inbounds con_scale[row] = max(con_scale[row], abs(jac.V[i]))
    end
end
function _set_scaling!(con_scale::AbstractVector, jac::Matrix)
    for row in 1:size(jac, 1)
        for col in 1:size(jac, 2)
            @inbounds con_scale[row] = max(con_scale[row], abs(jac[row, col]))
        end
    end
end

"""
    scale_constraints!(
        nlp::AbstractNLPModel,
        con_scale::AbstractVector,
        jac::AbstractMatrix;
        max_gradient=1e-8,
    )

Compute the scaling of the constraints associated
to the nonlinear model `nlp`. By default, Ipopt's scaling
is applied. The user can write its own function to scale
appropriately any custom `AbstractNLPModel`.

### Notes

This function assumes that the Jacobian `jac` has been evaluated
before calling this function.

"""
function scale_constraints!(
    nlp::AbstractNLPModel,
    con_scale::AbstractVector,
    jac::AbstractMatrix;
    max_gradient=1e-8,
)
    fill!(con_scale, 0.0)
    _set_scaling!(con_scale, jac)
    con_scale .= min.(1.0, max_gradient ./ con_scale)
end

"""
    scale_objective(
        nlp::AbstractNLPModel,
        grad::AbstractVector;
        max_gradient=1e-8,
    )

Compute the scaling of the objective associated to the
nonlinear model `nlp`. By default, Ipopt's scaling
is applied. The user can write its own function to scale
appropriately the objective of any custom `AbstractNLPModel`.

### Notes

This function assumes that the gradient `gradient` has been evaluated
before calling this function.

"""
function scale_objective(
    nlp::AbstractNLPModel,
    grad::AbstractVector;
    max_gradient=1e-8,
)
    return min(1, max_gradient / norm(grad, Inf))
end

function get_index_constraints(nlp::AbstractNLPModel; fixed_variable_treatment=MAKE_PARAMETER)
    ind_ineq = findall(get_lcon(nlp) .!= get_ucon(nlp))
    xl = [get_lvar(nlp);view(get_lcon(nlp),ind_ineq)]
    xu = [get_uvar(nlp);view(get_ucon(nlp),ind_ineq)]
    if fixed_variable_treatment == MAKE_PARAMETER
        ind_fixed = findall(xl .== xu)
        ind_lb = findall((xl .!= -Inf) .* (xl .!= xu))
        ind_ub = findall((xu .!=  Inf) .* (xl .!= xu))
    else
        ind_fixed = Int[]
        ind_lb = findall(xl .!=-Inf)
        ind_ub = findall(xu .!= Inf)
    end

    ind_llb = findall((get_lvar(nlp) .== -Inf).*(get_uvar(nlp) .!= Inf))
    ind_uub = findall((get_lvar(nlp) .!= -Inf).*(get_uvar(nlp) .== Inf))

    # Return named tuple
    return (
        ind_ineq=ind_ineq,
        ind_fixed=ind_fixed,
        ind_lb=ind_lb,
        ind_ub=ind_ub,
        ind_llb=ind_llb,
        ind_uub=ind_uub,
    )
end


text:
######################################################
#
# Using adjacent frequencies for noise estimate
#
######################################################

"""
    beamformer_lcmv(s::SSR, l::Leadfield)

Linearly constrained minimum variance (LCMV) beamformer for epoched data.

NAI is ratio between frequeny of interest (`foi`) and specified control (`foi + noise_delta`) frequency.

### Input

* s = stimulus condition SSR data with epochs pre calculated
* l = leadfield information

* foi = frequency of interest for cross power spectral density calculations
* fs = sample rate
* n_epochs = number of epochs to average down to with aim of reducing noise

"""
function beamformer_lcmv(s::SSR, l::Leadfield;
                         foi::Real=modulationrate(s), fs::Real=samplingrate(s), n_epochs::Int=0,
                         freq_pm::Real = 0.5, noise_delta::Real=4.0, bilateral::Real = 15, kwargs...)

    Logging.info("Performing LCMV beamforming on signal data with $(foi + noise_delta) Hz as reference")

    if !haskey(s.processing, "epochs")
          Logging.critical("Epochs not calculated")
    end

    if n_epochs > 0
        s.processing["epochs"] = reduce_epochs(s.processing["epochs"], n_epochs)
    end

    l = match_leadfield(l, s)

    C = cross_spectral_density(s.processing["epochs"], foi - freq_pm, foi + freq_pm, fs; kwargs...)
    Q = cross_spectral_density(s.processing["epochs"], foi + noise_delta - freq_pm, foi + noise_delta + freq_pm, fs, ignore = foi; kwargs...)

    @assert size(C) == size(Q)
    @assert C != Q

    Logging.debug("Covariance matrices calculated and of size $(size(Q))")

    V, N, NAI = beamformer_lcmv(C, Q, l.L, l.x, l.y, l.z, bilateral; kwargs...)

    VolumeImage(vec(NAI), "NAI", l.x, l.y, l.z, ones(size(vec(NAI))), "LCMV", Dict(), "Talairach")
end


text:
#=

  Euler #43 in Julia.

  """
  The number, 1406357289, is a 0 to 9 pandigital number because it is made up of
  each of the digits 0 to 9 in some order, but it also has a rather interesting
  sub-string divisibility property.

  Let d1 be the 1st digit, d2 be the 2nd digit, and so on. In this way, we
  note the following:

      * d2d3d4=406 is divisible by 2
      * d3d4d5=063 is divisible by 3
      * d4d5d6=635 is divisible by 5
      * d5d6d7=357 is divisible by 7
      * d6d7d8=572 is divisible by 11
      * d7d8d9=728 is divisible by 13
      * d8d9d10=289 is divisible by 17

  Find the sum of all 0 to 9 pandigital numbers with this property.
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# Using all_permutations:
# 1.95482710s
function euler43a()
    primes = [2,3,5,7,11,13,17];
    perms = all_permutations(Array(0:9))
    sum = 0
    for p in perms
        i = 0;
        found = true
        while i < 7 && found == true
            if (100*p[i+2] + 10*p[i+3] + p[i+4]) % primes[i+1] != 0
                found = false
                break
            end
            i += 1
        end

        if found
            sum += parse(Int,join(p,""))
        end
    end

    return sum;
end

# Using next_permutation
# 0.05306371s
function euler43b()
    primes = [2,3,5,7,11,13,17]
    p = [1,0,2,3,4,5,6,7,8,9]
    sum = 0
    rev = Array(9:-1:0)
     while p != rev
        i = 0
        found = true
        while i < 7 && found == true
            if (100*p[i+2] + 10*p[i+3] + p[i+4]) % primes[i+1] != 0
                found = false
                break
            end
            i += 1
        end

        if found
            sum += parse(Int,join(p,""))
        end
        p = next_permutation(p)
    end

    return sum
end

# run_euler(euler43a)
run_euler(euler43b)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using StaticArrays
struct Dogs{M<:Motion,N,D<:Dog}
    member::SizedVector{N,D}
    motion::M
end # struct

function Base.getindex( d::Dogs, i )
    return d.member[i]
end # function

function Base.getproperty( d::Dogs, s::Symbol )
    if s in fieldnames(Dogs)
        return getfield( d, s )
    else
        return Base.getproperty( d, Val(s) )
    end # if
end # function

function Base.getproperty( d::Dogs{M,N,D}, ::Val{:positions} ) where {M<:Motion,N,D<:Dog}
    T = typeof(d.member[1].position)
    # TODO: gefällt mir nicht (allocations?)
    return SizedVector{N,T}( [dog.position for dog in d.member] )
end # function

function Base.iterate(ds::Dogs, state = nothing)
    if state === nothing
        Base.iterate( ds.member )
    else
        Base.iterate( ds.member, state )
    end
end # function

function Base.length(ds::Dogs)
    Base.length( ds.member )
end # function

function Base.eltype(ds::Dogs)
    Base.eltype( ds.member )
end # function


text:
## IMPORTS ##
using Plots
using FFTW
using NPZ
include(joinpath(@__DIR__,"model.jl"))

## EQUILIBRIUM CORRELATIONS ##

function plot_density(dmrg_results::DMRGResults)
    n = expect(dmrg_results.ground_state, "Ntot")
    plot(1:length(n), n)
    ylims!((0.2,1.2))
    ylabel!("⟨n⟩")
    xlabel!("Site")
    title!("Electron density")
end

function plot_phonon_flux(dmrg_results::DMRGResults)
    n = expect(dmrg_results.ground_state, "Nb")
    plot(1:length(n), n)
    ylims!((0.2,1.2))
    ylabel!("⟨nb⟩")
    xlabel!("Site")
    title!("Phonon density")
end

function plot_equilibrium_correlations(dmrg_results::DMRGResults, corrtype::String,
                                        HH::HubbardHolsteinModel;
                                        start=nothing, stop=nothing)

    N = length(HH.sites)
    if isnothing(start)
        start = floor(Int, N/4)
    end
    if isnothing(stop)
        stop = ceil(Int, N/4*3)
    end

    corrs = equilibrium_correlations(dmrg_results,corrtype,HH,start,stop)

    plot(log10.(collect(0:stop-start)), log10.(abs.(corrs)))
    title!(corrtype*"-"*corrtype*" correlation")
    xlabel!("Distance from centre site (log10)")
    ylabel!("Correlation (log10)")
end

## CHECKING TEBD RESULTS ##

function plot_entropy(tebd_results::TEBDResults)
    ent = tebd_results.entropy
    niters = size(ent)[2]
    ϕ_entropy = ent[1,:]
    ψ_entropy = ent[2,:]
    plot(1:niters, ϕ_entropy, label="ϕ(t)")
    plot!(1:niters, ψ_entropy, label="ψ(t)")
    title!("Von Neumann Entropy")
    xlabel!("Iteration")
end

function plot_overlap(tebd_results::TEBDResults)
    plot(1:length(tebd_results.self_overlap), tebd_results.self_overlap)
end

## TIME-DEPENDENT CORRELATIONS ##

function plot_correlation_function(tebd_results::TEBDResults)
    heatmap(LinearAlgebra.norm.(tebd_results.corrs'),c=:heat)
    title!("Correlation function")
    xlabel!("Site")
    ylabel!("Time")
end

function make_spectral_fcn(corrs, p::Parameters; left_offset::Int=0)
    f1 = fftshift(fft(p.τ * corrs)')
    f2 = fftshift(fft(p.τ * reverse(corrs, dims=1))')

    qs = 2 * π * fftshift(fftfreq(p.N, 1))
    ωs = 2 * π * fftshift(fftfreq(size(corrs)[2], 1/p.τ))

    # Compensate for midpoint
    sqw = zeros(size(f1))
    for i in 1:p.N
        sqw[:,i] = imag.(exp.(1im * qs[i] * (p.mid-left_offset)) * f1[:,i])
                    + imag.(exp.(1im * qs[i] * (p.N-1-p.mid-left_offset)) * f2[:,i])
    end

    # If working with an even number of sites, must average across midline
    if iseven(p.N)
        sqw = (reverse(sqw, dims=2) + sqw) ./2
    end

    return real.(sqw)/π, ωs, qs
end

function decay(u;shape="exponential")
    if shape=="exponential"
        t = collect(1:(length(u)))
        λ = 0.1/sqrt(length(t))
        v = exp.(-λ*t)
    else
        @error "Not implemented"
    end
    u .* v
end

function plot_spectral_function(tebd_results::TEBDResults, p::Parameters;
                                smooth_signal=true, lims=nothing)
    corrs = tebd_results.corrs # num_time_steps x num_sites
    N, U, t, ω, g0, g1 = p.N, p.U, p.t, p.ω, p.g0, p.g1
    # Optionally convolve raw time data with decaying exponential?
    if smooth_signal
        # decaying exponential
        corrs = reverse(hcat(decay.(eachrow(corrs))...)',dims=2)
        # pad with zeros
        corrs = hcat(corrs, zeros(size(corrs)))
    end
    # Calculate the spectral function
    ff, ωs, qs = make_spectral_fcn(corrs, p)

    # Zoom in to the relevant bit
    if isnothing(lims)
        nstep = size(corrs)[2]
        lims = (floor(Int, nstep/2) - floor(Int, 0.07*nstep),floor(Int, nstep/2) + floor(Int, 0.07*nstep))
    end
    ff = ff[lims[1]:lims[2],:]
    ωs = ωs[lims[1]:lims[2]]

    # Plot
    maxval = maximum(abs.(ff))
    heatmap(qs, ωs, abs.(ff), c=:bwr, clims=(-maxval, maxval))
    title!("N=$N, U=$U, t=$t, ω=$ω, g0=$g0, g1=$g1")
    xlabel!("Momentum")
    ylabel!("Frequency")
end

function plot_spectral_function_slice(tebd_results::TEBDResults, p::Parameters; slice=0)
    findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2]

    ff, ωs, qs = make_spectral_fcn(tebd_results.corrs, p)
    ω = findnearest(ωs,slice)
    scatter(qs, abs.(ff[ω,:]),color="blue", label=nothing)
    plot!(qs, abs.(ff[ω,:]),color="blue", label=nothing)
    title!("Slice of S(q,ω) at ω=$slice")
end

function compare_to_ED(tebd_results::TEBDResults, p::Parameters)
    @assert p.τ == 0.01 # Time scale of the ED results
    corrs = tebd_results.corrs
    ED_corrs = npzread("/Users/nicole/Dropbox/Grad school/Devereaux lab/ITensors.jl/examples/dmrg/ed.npy")
    numsteps = size(corrs)[2]
    plot(1:numsteps, real.(corrs[4,:]), label="DMRG real part")
    plot!(1:numsteps, real.(ED_corrs[1:numsteps]), label="ED real part")
    plot!(1:numsteps, imag.(corrs[4,:]), label="DMRG complex part")
    plot!(1:numsteps, imag.(ED_corrs[1:numsteps]), label="ED complex part")
end


text:
module TestCycleGapLas

using Test
import GraphIdx

include("../src/cyclegap.jl")
include("demo3x7.jl")
n1, n2 = size(y)

@testset "GapLas RotateDemo 3x7          " begin
    graph = GraphIdx.Grid.GridGraph(n1, n2)
    lam = GraphIdx.Const(1.0)
    x = CycleGap.gaplas(copy(y), graph, lam, max_iter = 3)
    @test size(x) == size(y)
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
using Ruska
using Documenter

DocMeta.setdocmeta!(Ruska, :DocTestSetup, :(using Ruska); recursive=true)

makedocs(;
    modules=[Ruska],
    authors="Chen Huang <chen1huang2@gmail.com> and contributors",
    repo="https://github.com/chenspc/Ruska.jl/blob/{commit}{path}#{line}",
    sitename="Ruska.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chenspc.github.io/Ruska.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/chenspc/Ruska.jl",
)


text:
using DelimitedFiles
fishdata = readdlm("llist.txt",',',Int)
@show fishdata

function agefish(fishdata,day)
    if day == 81
        @show length(fishdata)
        return
    end

    newfishes = 0
    for fishidx in eachindex(fishdata)
        fish = fishdata[fishidx]
        if fish == 0
            fish = 6
            newfishes += 1
        else
            fish -= 1
        end
        fishdata[fishidx] = fish
    end

    fishdata = [fishdata fill(8,1,newfishes)]
    #@show fishdata

    agefish(fishdata,day+1)
end

agefish(fishdata,1)