text:
module VerticesTests

using Contour, Test
using Base.MathConstants: π, φ
using LinearAlgebra: Diagonal
using OffsetArrays
using StaticArrays

# Setup test axes that will be shared among the tests

# Shift the axes so that they do not line up with
# integer values
Δ = 0.01
X = collect(0:Δ:4) .+ π
Y = collect(0:Δ:3) .+ φ

# TEST CASE 1
#
# f(x,y) = x^2 + y^2
#
Z = [(x^2 + y^2)::Float64 for x in X, y in Y]
h = rand() * (maximum(Z) - minimum(Z)) + minimum(Z)

contourlevels = Contour.contour(X, Y, Z, h)
for line in contourlevels.lines
    # Contour vertices lie on a circle around the origin
    for v in line.vertices
        @test isapprox(v[1]^2 + v[2]^2, h, atol=0.1Δ)
    end

    # coordinates() returns the correct values
    xs, ys = coordinates(line)
    xs .== [v[1] for v in line.vertices]
    ys .== [v[2] for v in line.vertices]
end

# TEST CASE 1.5
#
# Same as test case 1, but with a shifted center
#
h = 1
x0, y0 = (.5, -1.5)
Z = Float64[(x - x0)^2 + (y - y0)^2 for x in X, y in Y]

contourlevels = Contour.contour(X, Y, Z, h)
for line in contourlevels.lines
    for v in line.vertices
        @test isapprox((v[1] - x0)^2 + (v[2] - y0)^2, h, atol=0.01Δ)
    end
end

# TEST CASE 2
#
# Check that ambigious cells (5, 10) are handled correctly
# Case 5: z_center > h
X = [1.0, 2.0]
Y = [2.0, 1.0]

Z = float([1 0;
           0 1])
h = 0.1

lines = Contour.contour(X, Y, Z, h).lines
@test length(lines) == 2

for line in lines
    @test length(line.vertices) == 2
    d = line.vertices[2] - line.vertices[1]
    @test d[2] / d[1] ≈ -1.0
end

# Case 5: z_center < h
Z = float([1 0;
           0 1])
h = 0.9

lines = Contour.contour(X, Y, Z, h).lines

@test length(lines) == 2

for line in lines
    @test length(line.vertices) == 2
    d = line.vertices[2] - line.vertices[1]
    @test d[2] / d[1] ≈ 1.0
end

# Case 10: z_center > h
Z = float([0 1;
           1 0])
h = 0.1

lines = Contour.contour(X, Y, Z, h).lines
@test length(lines) == 2

for line in lines
    @test length(line.vertices) == 2
    d = line.vertices[2] - line.vertices[1]
    @test d[2] / d[1] ≈ 1.0
end

# Case 10: z_center < h

Z = float([0 1;
           1 0])
h = 0.9

lines = Contour.contour(X, Y, Z, h).lines
@test length(lines) == 2

for line in lines
    @test length(line.vertices) == 2
    d = line.vertices[2] - line.vertices[1]
    @test d[2] / d[1] ≈ -1.0
end

# Test curvilinear coordinates
θ = range(0.0, stop=2π,length=100)
R = range(1.0, stop=2.0, length=100)
ζ = ComplexF64[r*exp(im*ϕ) for ϕ in θ, r in R]
x, y, z = real.(ζ), imag.(ζ), abs.(ζ)

h = 1 + rand()
xs, ys = coordinates(contour(x, y, z, h).lines[1])
@test all(xs.^2 + ys.^2 .≈ h^2)

# Test offset arrays
offset_x, offset_y = -10, 27
z = cumsum(cumsum(randn(20,20); dims=1); dims=2)
zoff = OffsetArray(z, offset_x, offset_y)

x, y = axes(z)
xoff, yoff = axes(zoff)
curves = Contour.contour(x,y,z,0.5)
curves_off = Contour.contour(xoff, yoff, zoff, 0.5)

# sort offset and non-offset curves to the same order
offset = SVector(offset_x, offset_y)
lns = sort(Contour.lines(curves); by=c->sum(sum.(c.vertices.+[offset])))
lns_off = sort(Contour.lines(curves_off); by=c->sum(sum.(c.vertices)))

# verify that each line matches a possibly circularly shifted or reversed
# line from the offset array
opencurve(a) = first(a) == last(a) ? a[1:end-1] : a
for (c1, c2) in zip(lns, lns_off)
    o1 = opencurve(c1.vertices) .+ [offset]
    o2 = opencurve(c2.vertices)
    m = length(o1)
    @test m == length(o2)
    cshifts = [circshift(o2,i) for i=1:m]
    @test any(o1 ≈ c || o1 ≈ reverse(c) for c in cshifts)
end

# Test Known Bugs

# Issue #12
x = float(collect(1:3));
y = copy(x);
z = Diagonal(ones(3));
contours(x, y, z)

# Test handling of saddle points
#
# f(x,y) = x^2 - y^2
#
Δ = 0.01
X = collect(-3:Δ:3)
Y = collect(-3:Δ:3)

Z = [(x^2 - y^2)::Float64 for x in X, y in Y]
h = rand() * (maximum(Z) - minimum(Z)) + minimum(Z)

contourlevels = Contour.contour(X, Y, Z, h)
for line in contourlevels.lines
    # Contour vertices lie on a circle around the origin
    for v in line.vertices
        @test isapprox(v[1]^2 - v[2]^2, h, atol=0.1Δ)
    end

    # coordinates() returns the correct values
    xs, ys = coordinates(line)
    xs .== [v[1] for v in line.vertices]
    ys .== [v[2] for v in line.vertices]
end

# Test range API
#
# f(x,y) = x^2 - y^2
#
Δ = 0.01
X = -3:Δ:3
Y = -3:Δ:3

Z = [(x^2 - y^2)::Float64 for x in X, y in Y]
h = rand() * (maximum(Z) - minimum(Z)) + minimum(Z)

contourlevels = Contour.contour(X, Y, Z, h)
for line in contourlevels.lines
    # Contour vertices lie on a circle around the origin
    for v in line.vertices
        @test isapprox(v[1]^2 - v[2]^2, h, atol=0.1Δ)
    end

    # coordinates() returns the correct values
    xs, ys = coordinates(line)
    xs .== [v[1] for v in line.vertices]
    ys .== [v[2] for v in line.vertices]
end

# Test that closed contours are identified correctly
# when ambiguous cells are involved

Z = float([1 1 1 1 1 1
           1 0 0 1 1 1
           1 0 0 1 1 1
           1 1 1 0 0 1
           1 1 1 0 0 1
           1 1 1 1 1 1])

X = Y = collect(0:0.2:1)
h = 0.75
contourlevels = Contour.contour(X, Y, Z, h)

# There should be only one closed contour
@test length(contourlevels.lines) == 1

# Test contour location on a realistic dataset
include("testdata.jl")
cts = Contour.contours(x, y, z)
@test length(cts.contours) == 10
cts_ct = (8, 8, 8, 8, 126, 7, 5, 5, 5, 4)

# Length need to be sorted as hashing might change see #62
lines_ct = sort!.([[138, 220, 469, 138, 469, 208, 143, 143],
            [220, 475, 140, 210, 146, 475, 140, 146],
            [222, 481, 140, 214, 481, 145, 140, 145],
            [228, 485, 214, 142, 485, 142, 146, 146],
            [9, 7, 515, 501, 9, 9, 9, 7, 7, 9, 7, 15, 7, 18, 39, 9, 7, 7, 9, 9, 7, 7, 7, 7, 9, 7, 18, 9, 7, 7, 9, 9, 35, 9, 9, 9, 7, 7, 9, 35, 9, 7, 9, 7, 9, 9, 7, 7, 7, 9, 7, 7, 7, 7, 9, 7, 9, 7, 9, 7, 7, 9, 9, 9, 7, 9, 7, 7, 9, 7, 7, 9, 7, 9, 7, 7, 7, 7, 7, 9, 7, 7, 7, 9, 7, 7, 7, 7, 9, 7, 7, 7, 7, 7, 7, 7, 9, 7, 9, 9, 7, 7, 7, 7, 7, 7, 9, 7, 7, 7, 7, 7, 7, 7, 9, 7, 9, 7, 7, 7, 7, 9, 7, 7, 7, 7],
            [5, 17, 14, 34, 34, 14, 5],
            [29, 12, 29, 15, 12],
            [26, 10, 13, 26, 10],
            [11, 11, 23, 11, 21],
            [18, 19, 7, 7]])
for i in eachindex(cts_ct)
    @test length(cts.contours[i].lines) == cts_ct[i]
    @test all(lines_ct[i] .== sort!([length(c.vertices) for c in cts.contours[i].lines]))
end

# support non-float z
using StatsBase

N = 10000
x = randn(N)
y = randn(N)
H = fit(Histogram, (x, y), closed = :left)
contours(midpoints.(H.edges)..., H.weights)

end


text:
include("..\\temperature_render.jl")
include("..\\good_plots\\line_and_render.jl")

# choosing the range to animate over
# vars = 75.0:1.0:85.0 # changing angle small
# vars = [-179.0:1.0:-1.0;1.0:1.0:179.0] # changing angle full
vars = 5:2.5:85                 # smaller range of angles
# vars = 1.0:5.0:200 # changing mass
# vars = 0:0.1:0.998 # changing spin

# generating frames
n_frames = length(vars)
anim = @animate for (i, var) in enumerate(vars)
    print("Frame: $(i-1) / $n_frames complete")
    fvar = Float64(var)
    # temperature_render(
    #                     mass=10,
    #                     spin=var,
    #                     tolerance=1e-8,
    #                     resolution=1080,
    #                     size_multiplier=2
    #                     )
    combined_plot(
                    mass=10,
                    spin=0.998,
                    obs_angle=fvar,
                    tolerance=1e-12,
                    size_multiplier=1,
                    fov=10,
                    dtmax=20
    )
end

# saving animation
gif(anim, "combined_gif.gif", fps=5)

ext: jl


text:
# hann.jl
# Hann window function

"""
    hann(x::AbstractFloat)
Hann window function (often called "Hanning" window function).
This version is `0.5 * (1 + cospi(2x))` and used in `[-0.5,+0.5]`.
https://en.wikipedia.org/wiki/Hann_function
"""
hann(x::AbstractFloat) = 0.5 * (1 + cospi(2x))

"""
    hann(n::Int)
Return vector of `n ≥ 2` samples of Hann window function
equally spaced over `[-0.5,+0.5]`.
The first and last samples are zero for `n ≥ 2`.

Caution.
Matlab has both `hann` and `hanning` functions for which
`[0; hanning(n); 0] == hann(n+2)`.

`DSP.hanning(n)` == Matlab.hann(n)`
whereas
`Sound.hann(n) == Matlab.hanning(n)`.
"""
hann(n::Int) = hann.(LinRange(-0.5, 0.5, n+2))[2:end-1]


text:
using Tensorial
using Test, Random

using LinearAlgebra: Symmetric
using StaticArrays: SArray, SVector, SOneTo, SUnitRange

include("Symmetry.jl")
include("Space.jl")
include("AbstractTensor.jl")
include("Tensor.jl")
include("permute.jl")
include("einsum.jl")
include("ops.jl")
include("continuum_mechanics.jl")
include("inv.jl")
include("voigt.jl")
include("ad.jl")
include("broadcast.jl")

include("quaternion.jl")


text:
# See https://web.stanford.edu/~lcambier/fast/demo.php
@testset "Optimize Stock with $solver" for solver in lp_solvers
    isclp(solver) && continue
    numScen = 2
    C = 1
    P = 2
    d = [2, 3]
    @testset "Optimize Stock directly" begin
        # With min ⟨[1, 1], [x, θ]⟩ s.t. no constraint with x >= 0, θ free
        # Clp is unable to return the unbounded ray which should be [0, -1]
        m1 = StructuredModel(num_scenarios=numScen)
        @variable(m1, x >= 0)
        @objective(m1, Min, C * x)

        for ξ in 1:numScen
            m2 = StructuredModel(parent=m1, prob=1/2, id=ξ)
            @variable(m2, s >= 0)
            @constraints m2 begin
                s <= d[ξ]
                s <= x
            end
            @objective(m2, Max, P * s)
        end

        num_stages = 2
        cutmode = StructProg.AvgCutGenerator()
        K = 2
        pereiracoef = 0.1

        # detectlb should change anything here since m2 is min -2s and s is bounded above by x with is unknown so it is unbounded
        for forwardcuts in [false, true]
            for detectlb in [false, true]
                sp = SOI.stochasticprogram(m1, num_stages, solver, AvgCutPruningAlgo(-1), cutmode, detectlb)
                @test SOI.get(sp, SOI.NumberOfPaths(0)) == 1
                @test SOI.get(sp, SOI.NumberOfPaths(1)) == 2
                @test SOI.get(sp, SOI.NumberOfPaths(2)) == 2
                @test sprint(show, StructDualDynProg.StructProg.nodedata(sp, 1)) == "Node of 1 variables\n"
                algo = SDDP.Algorithm(K = K, forwardcuts = forwardcuts, backwardcuts = !forwardcuts)
                sol = SOI.optimize!(sp, algo, SOI.Pereira(0.1) | SOI.IterLimit(10), 0)
                sstats = sprint(show, SOI.last_result(sol))
                @test occursin("Lower Bound: ", sstats)
                @test occursin("Upper Bound: ", sstats)

                # K = 10 is a multiple of 2 so with ProbaPathSampler(true), the sampling is deterministic
                # therefore we can test for sol.attrs[:niter]
                #
                # Iteration | Status    |  LB  |  UB  | #OC |
                #     1     | Unbounded | -Inf |  0.0 |  1  |
                #     2     | Unbounded | -Inf | -5.0 |  1  |
                #     3     | Optimal   | -2.5 | -2.0 |  1  |
                #     4     | Optimal   | -2.0 | -2.0 |  0  |
                @test SOI.niterations(sol) == 4
                @test SOI.last_result(sol).status == :Optimal
                @test SOI.last_result(sol).lowerbound == -2.0
                StructProg.clear(m1)
            end
        end
    end

    # See https://github.com/blegat/StructDualDynProg.jl/issues/10
    if !isgrb(solver) && !iscpx(solver)
        @testset "Add a scenario latter" begin
            m1 = StructuredModel(num_scenarios=numScen)
            @variable(m1, x >= 0)
            @objective(m1, Min, C * x)

            m2 = StructuredModel(parent=m1, prob=1., id=1)
            @variable(m2, s >= 0)
            @constraints m2 begin
                s <= d[2]
                s <= x
            end
            @objective(m2, Max, P * s)

            num_stages = 2
            # Multicut wouldn't work since we are adding a node
            cutmode = StructProg.AvgCutGenerator()
            K = 2
            pereiracoef = 0.1

            sp = SOI.stochasticprogram(m1, num_stages, solver, AvgCutPruningAlgo(-1), cutmode)
            sol = SOI.optimize!(sp, SDDP.Algorithm(K=K), SOI.Pereira(0.1) | SOI.IterLimit(10), 0)
            @test SOI.niterations(sol) == 3
            @test SOI.last_result(sol).status == :Optimal
            @test SOI.last_result(sol).lowerbound == -3.0

            m3 = StructuredModel(parent=m1, id=2)
            @variable(m3, s >= 0)
            @constraints m3 begin
                s <= d[1]
                s <= x
            end
            @objective(m3, Max, P * s)

            root = 1
            newnode = StructProg.createnode(sp, m3, 1, 1, solver, root, AvgCutPruningAlgo(-1), cutmode)
            @test length(sp.out_transitions[1]) == 1
            SOI.set!(sp, SOI.Probability(), sp.out_transitions[1][1], 1/2)
            SOI.add_scenario_transition!(sp, root, newnode, 1/2)
            sol = SOI.optimize!(sp, SDDP.Algorithm(K=K), SOI.Pereira(0.1) | SOI.IterLimit(10), 0)
            # 2 on Mac OS and Windows, 3 otherwise
            @test 2 <= SOI.niterations(sol) <= 3
            @test SOI.last_result(sol).status == :Optimal
            @test SOI.last_result(sol).lowerbound == -2.0
        end
    end
end


text:
#
# grad_tools.jl
#
# Provide tools to load the instrumental parameters and the
# calibrated data, and for the calculus of the RHAPSODIE data fidelity term.
#
# ------------------------------------------------
#
# This file is part of Rhapsodie
#
#
# Copyright (c) 2017-2021 Laurence Denneulin (see LICENCE.md)
#

#------------------------------------------------

const ImageInterpolator{T<:AbstractFloat, K<:Kernel{T}} = TwoDimensionalTransformInterpolator{T,K,K}
const MyKer = LinearInterpolators.CatmullRomSpline(Float64, LinearInterpolators.Flat)
#const MyKer = LinearInterpolators.RectangularSpline(Float64, LinearInterpolators.Flat)
#const MyKer = LinearInterpolators.LinearSpline(Float64, LinearInterpolators.Flat)

struct parameters_table{T<: AbstractFloat}
     cols::NTuple{3,Int64} #Imput size (reconstruction)
     rows::NTuple{2,Int64} #Output size (data)
     dataset_length::Int64
     Nframe::Int64
     Nrot::Int64
     Nangle::Int64
     v::Vector{NTuple{2, NTuple{3, T}}}
     indices::Array{Int64,2}
     center::Array{T,1}
     psf_center::NTuple{2,Array{T,1}}
     epsilon::Vector{NTuple{2,Array{T,1}}}
     derotang::Vector{T}
end

struct FieldTransformOperator{T<:AbstractFloat, L<:Mapping, R<:Mapping} <: LinearMapping
    cols::NTuple{3,Int}
    rows::NTuple{2,Int}
    v_l::NTuple{3,T}
    v_r::NTuple{3,T}
    H_l::L
    H_r::R
end

struct data_table{T<:AbstractFloat,A<:FieldTransformOperator{T}} #autres types
    data::Array{T, 2};
    weights::Array{T, 2};
    H::A
end

const Trans_Table = Vector{NTuple{2,AffineTransform2D}}(); #Contains all the affine transform used
const Parameters = parameters_table[];
get_par()::parameters_table = Parameters[1];
const dataset = data_table[];

const PSF_save = Vector{Array{Float64,2}}();
get_PSF()=PSF_save[1];
const EPSILON_save = Array{Float64,1}(undef,1);
function set_epsilon(epsilon)
    EPSILON_save[1]=epsilon
end
get_epsilon()::Float64=EPSILON_save[1];
const PRECOND_SAVE= Vector{Any}();
U()=PRECOND_SAVE[1];

const MASK_save = Vector{Array{Float64,2}}();
get_MASK()=MASK_save[1];

function load_parameters(size_data::NTuple{3,Int64},
                         Nframe::Int64,
                         Nrot::Int64,
                         Nangle::Int64,
                         center::Array{Float64,1},
                         psf_center::NTuple{2,Array{Float64,1}},
                         epsilon::Vector{NTuple{2,Array{Float64,1}}})
    sets_indices=Indices(size_data[3],Nrot,Nframe)
    sets_v=Set_Vi(sets_indices)
    bbox_output=(0,0)
    Id = AffineTransform2D{Float64}()
    for k=1:size_data[3]
        A_left=translate( epsilon[k][1][1], epsilon[k][1][2], Id)
        A_right=translate( epsilon[k][2][1], epsilon[k][2][2], Id)

        out_left=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]
        out_right=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]

        xmax=max(bbox_output[1],out_left[1],out_right[1]);
        ymax=max(bbox_output[2],out_left[2],out_right[2]);

        bbox_output =(xmax,ymax);
    end
    bbox_output = bbox_output .+ 4;
    push!(Parameters, parameters_table((bbox_output[1],bbox_output[2],3),
                                       (size_data[1], size_data[2]),
                                       size_data[3],
                                       Nframe,
                                       Nrot,
                                       Nangle,
                                       sets_v,
                                       sets_indices,
                                       center,
                                       psf_center,
                                       epsilon,
                                       Float64[]));

        newcenter= (bbox_output .+1)./2;
        centerdiff=newcenter .- (center[1], center[2])
        for k=1:size_data[3]
        A_left=translate( epsilon[k][1][1] + centerdiff[1], epsilon[k][1][2] + centerdiff[2], Id)
        A_right=translate( epsilon[k][2][1] + centerdiff[1], epsilon[k][2][2] + centerdiff[2], Id)
        push!(Trans_Table, (A_left, A_right))
    end
end

function load_parameters(size_data::NTuple{3,Int64},
                         Nframe::Int64,
                         Nrot::Int64,
                         Nangle::Int64,
                         center::Array{Float64,1},
                         psf_center::NTuple{2,Array{Float64,1}},
                         epsilon::Vector{NTuple{2,Array{Float64,1}}},
                         derotang::Vector{Float64})
    sets_indices=Indices(size_data[3],Nrot,Nframe)
    sets_v=Set_Vi(sets_indices)
    bbox_output=(0,0)
    Id = AffineTransform2D{Float64}()
    for k=1:size_data[3]
        A_left=TransRotate(Id, epsilon[k][1], derotang[k], center, center)
        A_right=TransRotate(Id, epsilon[k][2], derotang[k], center, center)

        out_left=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]
        out_right=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]

        xmax=max(bbox_output[1],out_left[1],out_right[1]);
        ymax=max(bbox_output[2],out_left[2],out_right[2]);

        bbox_output =(xmax,ymax);
    end
    bbox_output = bbox_output .+ 4;
    push!(Parameters, parameters_table((bbox_output[1],bbox_output[2],3),
                                       (size_data[1], size_data[2]),
                                       size_data[3],
                                       Nframe,
                                       Nrot,
                                       Nangle,
                                       sets_v,
                                       sets_indices,
                                       center,
                                       psf_center,
                                       epsilon,
                                       derotang));

    newcenter= (bbox_output .+1)./2
    for k=1:size_data[3]
        A_left=inv(TransRotate(Id, epsilon[k][1], derotang[k], center, newcenter))
        A_right=inv(TransRotate(Id, epsilon[k][2], derotang[k], center, newcenter))
        push!(Trans_Table, (A_left, A_right))
    end
end

function load_parameters(size_data::NTuple{3,Int64},
                         Nframe::Int64,
                         Nrot::Int64,
                         Nangle::Int64,
                         center::Array{Float64,1},
                         psf_center::NTuple{2,Array{Float64,1}},
                         epsilon::Vector{NTuple{2,Array{Float64,1}}},
                         v::Array{Float64,2})
    sets_indices=Indices(size_data[3],Nrot,Nframe)
    sets_v=Set_Vi(Nframe, size_data[3], v)
    bbox_output=(0,0)
    Id = AffineTransform2D{Float64}()
    for k=1:size_data[3]
        A_left=translate( epsilon[k][1][1], epsilon[k][1][2], Id)
        A_right=translate( epsilon[k][2][1], epsilon[k][2][2], Id)

        out_left=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]
        out_right=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]

        xmax=max(bbox_output[1],out_left[1],out_right[1]);
        ymax=max(bbox_output[2],out_left[2],out_right[2]);

        bbox_output =(xmax,ymax);
    end
    bbox_output = bbox_output .+ 4;
    push!(Parameters, parameters_table((bbox_output[1],bbox_output[2],3),
                                       (size_data[1], size_data[2]),
                                       size_data[3],
                                       Nframe,
                                       Nrot,
                                       Nangle,
                                       sets_v,
                                       sets_indices,
                                       center,
                                       psf_center,
                                       epsilon,
                                       Float64[]));

        newcenter= (bbox_output .+1)./2;
        centerdiff=newcenter .- (center[1], center[2])
        for k=1:size_data[3]
        A_left=translate( epsilon[k][1][1] + centerdiff[1], epsilon[k][1][2] + centerdiff[2], Id)
        A_right=translate( epsilon[k][2][1] + centerdiff[1], epsilon[k][2][2] + centerdiff[2], Id)
        push!(Trans_Table, (A_left, A_right))
    end
end

function load_parameters(size_data::NTuple{3,Int64},
                         Nframe::Int64,
                         Nrot::Int64,
                         Nangle::Int64,
                         center::Array{Float64,1},
                         psf_center::NTuple{2,Array{Float64,1}},
                         epsilon::Vector{NTuple{2,Array{Float64,1}}},
                         derotang::Vector{Float64},
                         v::Array{Float64,2})
    sets_indices=Indices(size_data[3],Nrot,Nframe)
    sets_v=Set_Vi(Nframe, size_data[3], v)
    bbox_output=(0,0)
    Id = AffineTransform2D{Float64}()
    for k=1:size_data[3]
        A_left=TransRotate(Id, epsilon[k][1], derotang[k], center, center)
        A_right=TransRotate(Id, epsilon[k][2], derotang[k], center, center)

        out_left=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]
        out_right=bbox_size((size_data[1], Int64(size_data[2]/2)), A_left)[1]

        xmax=max(bbox_output[1],out_left[1],out_right[1]);
        ymax=max(bbox_output[2],out_left[2],out_right[2]);

        bbox_output =(xmax,ymax);
    end
    bbox_output = bbox_output .+ 4;
    push!(Parameters, parameters_table((bbox_output[1],bbox_output[2],3),
                                       (size_data[1], size_data[2]),
                                       size_data[3],
                                       Nframe,
                                       Nrot,
                                       Nangle,
                                       sets_v,
                                       sets_indices,
                                       center,
                                       psf_center,
                                       epsilon,
                                       derotang));

    newcenter= (bbox_output .+1)./2
    for k=1:size_data[3]
        A_left=inv(TransRotate(Id, epsilon[k][1], derotang[k], center, newcenter))
        A_right=inv(TransRotate(Id, epsilon[k][2], derotang[k], center, newcenter))
        push!(Trans_Table, (A_left, A_right))
    end
end

function Load_Data(name_data, name_weight, name_psf)
	data=read(FitsArray, name_data);
	weight=read(FitsArray, name_weight);
	psf=read(FitsArray, name_psf);
	A_l=set_fft_op((psf[1:end÷2,:]'[:,:]), get_par().psf_center[1]);
	A_r=set_fft_op((psf[end÷2+1:end,:]'[:,:]), get_par().psf_center[1]);
	ker= MyKer;
	SetCropOperator()
    for k=1:size(data)[3]
        output_size=(get_par().rows[1], Int64(get_par().rows[2]/2));
        input_size= get_par().cols[1:2];
    	T1=TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, Trans_Table[k][1])
    	T2=TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, Trans_Table[k][2])
	    push!(dataset, data_table((data[:,:,k]')[:,:],
	                              (weight[:,:,k]')[:,:],
	                              FieldTransformOperator(get_par().cols,
	                                                     get_par().rows,
	                                                     get_par().v[k][1],
	                                                     get_par().v[k][2],
	                                                     T1,
	                                                     T2,
	                                                     A_l)));#,
	                                                     #A_r)));
	end
end

function Load_Data(name_data, name_weight)
	data=read(FitsArray, name_data);
	weight=read(FitsArray, name_weight);
	ker= MyKer;
	SetCropOperator()
    for k=1:size(data)[3]
        output_size=(get_par().rows[1], Int64(get_par().rows[2]/2));
        input_size= get_par().cols[1:2];
    	T1=TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, Trans_Table[k][1])
    	T2=TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, Trans_Table[k][2])
    	promote
	    push!(dataset, data_table((data[:,:,k]')[:,:],
	                              (weight[:,:,k]')[:,:],
	                              FieldTransformOperator(get_par().cols,
	                                                     get_par().rows,
	                                                     get_par().v[k][1],
	                                                     get_par().v[k][2],
	                                                     T1,
	                                                     T2,
	                                                     LazyAlgebra.Id)));
	end
end

function Indices(S::Int64,Nrot::Int64,Nframe::Int64)
#S = le nombre total de frame gauche ou droite
#Nrot : Nombre de positions du derotateur
#Nframe : Nombre de frame par positions de la lame demi-onde

	Nangle=Int32.(S/(Nrot*Nframe*4)) #Nombre de rotations de la lame demi-onde
	INDICES=zeros(4,Nframe*Nangle*Nrot)
	for i=1:4
		ind=repeat(range(0,stop=4*Nframe*(Nrot*Nangle-1),length=Nrot*Nangle), inner=Nframe)+(Nframe*i .-mod.(range(1,stop=Nframe*Nangle*Nrot,length=Nframe*Nangle*Nrot),Nframe))
		INDICES[i,:]=ind
	end
	INDICES=round.(Int64, INDICES);
	return INDICES
end

function Set_Vi(Indices::AbstractArray{Int64,2}; alpha=[0, pi/4, pi/8, 3*pi/8], psi=[0,pi/2])
	J(a)=[cos(2*a) sin(2*a); sin(2*a) -cos(2*a)]
	P(a)=[cos(a), -sin(a)]
    v=Vector{NTuple{2, NTuple{3, Float64}}}(undef, length(Indices));
	for i=1:4
		v1=J(alpha[i])*P(psi[1])
		vn1=norm(v1)^2;
		v2=J(alpha[i])*P(psi[2])
		vn2=norm(v2)^2;
		for k=1:length(Indices[i,:])
		v[Indices[i,k]] =((round(vn1/2, digits=4), round((abs(v1[1])^2-abs(v1[2])^2)/2, digits=4), round(real(v1[1]*v1[2]), digits=4)),
		                  (round(vn2/2, digits=4), round((abs(v2[1])^2-abs(v2[2])^2)/2, digits=4), round(real(v2[1]*v2[2]), digits=4)));
        end
	end
	return v
end

function Set_Vi(Nframe::Int, dataset_length::Int, mueller_instru::AbstractArray{Float64,2})
    @assert dataset_length == Nframe * size(mueller_instru)[1]

    v=Vector{NTuple{2, NTuple{3, Float64}}}(undef, dataset_length);
    for k=1:size(mueller_instru)[1]
        for l=1:Nframe
            v[Nframe*(k-1) + l] =((mueller_instru[k,1], mueller_instru[k,2], mueller_instru[k,3]),
                                  (mueller_instru[k,4], mueller_instru[k,5], mueller_instru[k,6]));
        end
    end
	return v
end

function reset_instrument(V::Vector{NTuple{2, NTuple{3, T}}})  where {T <: AbstractFloat}
    push!(Parameters, parameters_table(get_par().cols,
                                       get_par().rows,
                                       get_par().dataset_length,
                                       get_par().Nframe,
                                       get_par().Nrot,
                                       get_par().Nangle,
                                       V,
                                       get_par().indices,
                                       get_par().center,
                                       get_par().psf_center,
                                       get_par().epsilon,
                                       get_par().derotang));

    popfirst!(Parameters);
end

function TransRotate(A::AffineTransform2D{Float64}, EPSILON, ANGLE, CENTER, NEWCENTER)
	return translate(rotate(translate(CENTER[1]-EPSILON[1],CENTER[2]-EPSILON[2], A), ANGLE),-NEWCENTER[1], -NEWCENTER[2])
end

function bbox_size(inp_dims::NTuple{2,Integer},
                  A::AffineTransform2D{Float64})

	xmin=typemax(Float64);
	ymin=typemax(Float64);
	xmax=typemin(Float64);
	ymax=typemin(Float64);

	xmin_int=typemin(Int64);
	ymin_int=typemin(Int64);
	xmax_int=typemin(Int64);
	ymax_int=typemin(Int64);

	width=typemin(Float64);
	height=typemin(Float64);

	Ind=[repeat(1:inp_dims[1], inner=inp_dims[2]) repeat(1:inp_dims[2], outer=inp_dims[1])]
	for it=1:(inp_dims[1]*inp_dims[2])
		(x,y)=A(Ind[it,1], Ind[it,2]);
		xmin=min(x, xmin);
		ymin=min(y, ymin);
		xmax=max(x, xmax);
		ymax=max(y, ymax);
	end
	xmin_int=floor(Int64,xmin);
	ymin_int=floor(Int64,ymin);
	xmax_int=ceil(Int64,xmax);
	ymax_int=ceil(Int64,ymax);

	width=xmax_int-xmin_int+1+50;
	height=ymax_int-ymin_int+1+50;
    out_dims=(width, height);

	return(out_dims, (xmin_int, xmax_int, ymin_int, ymax_int))
end

function set_fft_op(PSF::AbstractArray{T,2}, PSFCenter::AbstractArray{T,1}) where {T <: AbstractFloat}
 	MapSize=get_par().cols[1:2];
	MapCenter=floor.(MapSize./2).+1
	MAP=zeros(MapSize);
	ker = LinearInterpolators.CatmullRomSpline(Float64, LinearInterpolators.Flat)

	Id = AffineTransform2D{Float64}()
	Recentering=translate(-(MapCenter[1]-PSFCenter[1]), -(MapCenter[2]-PSFCenter[2]), Id)

	LazyAlgebra.apply!(MAP, ker, Recentering, PSF);
	MAP./=sum(MAP);
	push!(PSF_save,MAP);
	F=FFTOperator(MAP)
	FFT=F\Diag(F*ifftshift(MAP)) .*F;

	return FFT
end

function FieldTransformOperator(cols::NTuple{3,Int64},
                                rows::NTuple{2,Int64},
                                v_l::NTuple{3,T},
                                v_r::NTuple{3,T},
                                T_l::TwoDimensionalTransformInterpolator{T},
                                T_r::TwoDimensionalTransformInterpolator{T},
                                A::M) where {T <: AbstractFloat, M <: Mapping}
    FieldTransformOperator(cols, rows, v_l, v_r, T_l*A, T_r*A)
end

function vcreate(::Type{LazyAlgebra.Direct}, A::FieldTransformOperator{T},
                 x::AbstractArray{T,3}, scratch::Bool = false) where {T <: AbstractFloat}
    @assert !Base.has_offset_axes(x)
    @assert size(x) == A.cols
    Array{T,2}(undef, A.rows)
end

function vcreate(::Type{LazyAlgebra.Adjoint}, A::FieldTransformOperator{T},
                 x::AbstractArray{T,2}, scratch::Bool = false) where {T <: AbstractFloat}
    @assert !Base.has_offset_axes(x)
    @assert size(x) == A.rows
    Array{T,3}(undef, A.cols)
end

function apply!(α::Real,
                ::Type{LazyAlgebra.Direct},
                R::FieldTransformOperator{T},
                src::AbstractArray{T,3},
                scratch::Bool,
                β::Real,
                dst::AbstractArray{T,2}) where {T<:AbstractFloat}
    @assert β==0 && α==1
    @assert size(src) == R.cols
    @assert size(dst) == R.rows
    n = R.rows[2]
    @assert iseven(n)

    z_l = zeros(R.cols[1:2]);
    z_r = zeros(R.cols[1:2]);
    for i=1:length(R.v_l)
        # FIXME: use Generalized Matrix-Vector Multiplication of LazyAlgebra here.
        x_i = view(src,:,:,i)
	z_l .+= R.v_l[i]*x_i;
	z_r .+= R.v_r[i]*x_i;
    end
    dst[:, 1:(n÷2)  ] = R.H_l*z_l;
    dst[:, (n÷2)+1:n] = R.H_r*z_r;
    return dst
end

function apply!(α::Real,
                ::Type{LazyAlgebra.Adjoint},
                R::FieldTransformOperator{T},
                src::AbstractArray{T,2},
                scratch::Bool,
                β::Real,
                dst::AbstractArray{T,3}) where {T<:AbstractFloat}
    @assert β==0 && α==1
    @assert size(src) == R.rows
    @assert size(dst) == R.cols
    n = R.rows[2]
    @assert iseven(n)
    y_l = R.H_l'*view(src, :, 1:(n÷2))
    y_r = R.H_r'*view(src, :, (n÷2)+1:n)
    for i=1:length(R.v_l)
        #dst[:,:,i] = R.v_l[i]*y_l + R.v_r[i]*y_r;
        # FIXME: The following should be equivalent and much faster:
         vcombine!(view(dst,:,:,i), R.v_l[i], y_l, R.v_r[i], y_r)
    end
    return dst;
end

#=
function apply!(α::Real,
                ::Type{LazyAlgebra.Direct},
                R::FieldTransformOperator{T},
                src::AbstractArray{T,3},
                scratch::Bool,
                β::Real,
                dst::AbstractArray{T,2}) where {T<:AbstractFloat}

    @assert size(src) == R.cols
    @assert size(dst) == R.rows
    @assert β==0 && α==1
    vfill!(dst,0);

    dst_left=zeros(R.cols[1:2]);
    dst_right=zeros(R.cols[1:2]);

    for i=1:length(R.v_l)
        dst_left .+=(R.v_l[i]*src[:,:,i]);
        dst_right .+=(R.v_r[i]*src[:,:,i]);
    end
    dst .+= [R.T_l*R.A_l*dst_left R.T_r*R.A_r*dst_right];
    return dst
end

function apply!(α::Real,
                ::Type{LazyAlgebra.Adjoint},
                R::FieldTransformOperator{T},
                src::AbstractArray{T,2},
                scratch::Bool,
                β::Real,
                dst::AbstractArray{T,3}) where {T<:AbstractFloat}
    @assert β==0 && α==1
    @assert size(src) == R.rows
    @assert size(dst) == R.cols
    for i=1:length(R.v_l)
        dst_left=R.A_l'*R.T_l'*(R.v_l[i]*src[:,1:end÷2]);
        dst_right=R.A_r'*R.T_r'*(R.v_r[i]*src[:,end÷2+1:end]);
        dst[:,:,i] = dst_left + dst_right;
    end
    return dst;
end
=#

function mydot(x::AbstractArray{Tx,N},
               y::AbstractArray{Ty,N}) where {Tx<:AbstractFloat,Ty<:AbstractFloat,N}
    @assert axes(x) == axes(y)
    s = zero(promote_type(Tx,Ty))
    @inbounds @simd for i in eachindex(x,y)
       s += x[i]*y[i]
    end
    return s
end

function fg!(x::AbstractArray{T,3},g::AbstractArray{T,3}) where {T<:AbstractFloat}
    local f::Float64 = 0.0;
    vfill!(g,0.0)
    for k=1:length(dataset)
        if sum(dataset[k].weights) !=0
            f += fg!(x, g, dataset[k])
        end
    end
    #crop!(g)
    return f
end

function fg!(x::AbstractArray{T,3},g::AbstractArray{T,3}, d::data_table{T}) where {T<:AbstractFloat}
    r = d.H*x - d.data;
    wr = d.weights .* r;
    g .+= d.H'*wr;
    #f = vdot(Float64, r,wr); #FIXME : mettre une issue !
    f = mydot(r,wr);
    #crop!(g)
    return Float64(f)/2
end

#--------------------------------------------UTILS------------------------------------------------

function SetCropOperator()
    DSIZE=get_par().rows[1]
    MASK=ones(get_par().cols[1:2]);
    for k=1:length(Trans_Table)
    X1=Trans_Table[k][1](1,1);
    X2=Trans_Table[k][1](1,DSIZE);
    X3=Trans_Table[k][1](DSIZE, DSIZE);
    X4=Trans_Table[k][1](DSIZE, 1);
    Mask1=zeros(get_par().cols[1:2]);

    for i=1:get_par().cols[1]
	    for j=1:get_par().cols[2]
		    if ((i-X1[1])*(X2[1]-X1[1]) +(j-X1[2])*(X2[2]-X1[2]) >0)&&((i-X2[1])*(X3[1]-X2[1]) +(j-X2[2])*(X3[2]-X2[2]) >0)&&((i-X3[1])*(X4[1]-X3[1]) +(j-X3[2])*(X4[2]-X3[2]) >0) &&((i-X4[1])*(X1[1]-X4[1]) +(j-X4[2])*(X1[2]-X4[2]) >0)
		    Mask1[i,j,:] .=1;
		    end
	    end
    end
    X1=Trans_Table[k][2](1,1);
    X2=Trans_Table[k][2](1,DSIZE);
    X3=Trans_Table[k][2](DSIZE, DSIZE);
    X4=Trans_Table[k][2](DSIZE, 1);
    Mask2=zeros(get_par().cols[1:2]);

    for i=1:get_par().cols[1]
	    for j=1:get_par().cols[2]
		    if ((i-X1[1])*(X2[1]-X1[1]) +(j-X1[2])*(X2[2]-X1[2]) >0)&&((i-X2[1])*(X3[1]-X2[1]) +(j-X2[2])*(X3[2]-X2[2]) >0)&&((i-X3[1])*(X4[1]-X3[1]) +(j-X3[2])*(X4[2]-X3[2]) >0) &&((i-X4[1])*(X1[1]-X4[1]) +(j-X4[2])*(X1[2]-X4[2]) >0)
		    Mask2[i,j,:] .=1;
		    end
	    end
    end
    MASK .*= Mask1.*Mask2
    end
    #write(FitsFile, "MASK.fits", MASK, overwrite=true)
    push!(MASK_save, MASK)
end

function crop(X::M)  where {T<:AbstractFloat, M<:AbstractArray{T,2}}
    #@assert size(X) .==   get_par().cols
    xymin_1=findfirst(get_MASK() .!= 0.)
    xymin_2=findfirst(get_MASK()' .!= 0.)
    xymax_1=findlast(get_MASK() .!= 0.)
    xymax_2=findlast(get_MASK()' .!= 0.)

    xmin=min(xymin_1[1], xymin_2[2]);
    ymin=min(xymin_1[2], xymin_2[1]);
    xmax=max(xymax_1[1], xymax_2[2]);
    ymax=max(xymax_1[2], xymax_2[1]);

    Y=X[xmin:xmax, ymin:ymax];
    Y[.!isfinite.(Y)].=0
    return Y
end

function crop!(X::M)  where {T<:AbstractFloat, M<:AbstractArray{T,2}}
    #@assert size(X) .==   get_par().cols
    X[.!isfinite.(X)].=0
    X .*= get_MASK()
end

function crop(X::M)  where {T<:AbstractFloat, M<:AbstractArray{T,3}}
    #@assert size(X) .==   get_par().cols
    Y=copy(X);
    return crop!(Y)
end

function crop!(X::M)  where {T<:AbstractFloat, M<:AbstractArray{T,3}}
    #@assert size(X) .==   get_par().cols
    for k in size(X)[3]
        crop!(view(X,:,:,k))
    end
end

function crop(X::PolarimetricMap{T}) where {T<:AbstractFloat}
    #@assert size(X) .==   get_par().cols
    return PolarimetricMap(X.parameter_type,
                           crop(view(X.I,:,:)),
                           crop(view(X.Q,:,:)),
                           crop(view(X.U,:,:)),
                           crop(view(X.Iu,:,:)),
                           crop(view(X.Ip,:,:)),
                           crop(view(X.θ,:,:)))
end

function crop!(X::PolarimetricMap{T})  where {T<:AbstractFloat}
        crop!(view(X.I,:,:));
        crop!(view(X.Q,:,:));
        crop!(view(X.U,:,:));
        crop!(view(X.Iu,:,:));
        crop!(view(X.Ip,:,:));
        crop!(view(X.θ,:,:));
end

max_stars_repo_name: LaurenceDenneulin/Rhapsodie

max_issues_repo_name: LaurenceDenneulin/Rhapsodie

max_forks_repo_name: LaurenceDenneulin/Rhapsodie


text:
module AlgebraInterfacesTests

using Gridap.Algebra
using Test

a = allocate_vector(Vector{Int},10)
@test isa(a,Vector{Int})
@test length(a) == 10

a = allocate_vector(Vector{Int},1:10)
@test isa(a,Vector{Int})
@test length(a) == 10

a = zeros(4,6)

b = allocate_in_range(Vector{Int},a)
@test length(b) == 4

b = allocate_in_domain(Vector{Int},a)
@test length(b) == 6

fill_entries!(b,4)
@test all( b .== 4 )

a = rand(6)
b = rand(6)

copy_entries!(a,b)
@test all( a .== b )

a = rand(6)
c = copy(a)
b = rand(6)

add_entries!(a,b)
@test all( a .== ( c .+ b) )

a = rand(6)
c = copy(a)
b = rand(6)

add_entries!(a,b,-)
@test all( a .== ( c .- b) )

a = rand(6)
c = copy(a)
scale_entries!(a,10)
@test all( a .== 10*c)

a = rand(4,6)
c = rand(4)
d = copy(c)
b = rand(6)

muladd!(c,a,b)

@test all( c .≈ (d .+ a*b ) )

end # module


text:
using Flux
using Flux: @epochs
using Flux: params
using Flux: throttle
using Flux: glorot_uniform
using Flux: @functor
using Plots
using Distributions
using Zygote
using Flux.Optimise: update!
using ForwardDiff
using ReverseDiff
using FiniteDiff
using ZygoteRules

```
Workflow:

# 1. Define network architecture:
    # 1. define 1 LSTM-like layer
    # 2. define 1 dense layer with activation for the start
    # 3. define 1 dense layer without activation for the end
    # 4. connect them all together -> model
# 2. Define the analytical solution function
# 3. Define sampling function
# 4. Define the loss function
# 5. Train the model

```

# define a DGM layer
mutable struct LSTMLayer
    Ug :: Array{Float64}
    Uz :: Array{Float64}
    Ur :: Array{Float64}
    Uh :: Array{Float64}
    Wg :: Array{Float64}
    Wz :: Array{Float64}
    Wr :: Array{Float64}
    Wh :: Array{Float64}
    bg :: Array{Float64}
    bz :: Array{Float64}
    br :: Array{Float64}
    bh :: Array{Float64}
    act1
    act2
end

LSTMLayer(input_dim :: Integer, output_dim :: Integer, activation_1, activation_2) = LSTMLayer(glorot_uniform(output_dim, input_dim), glorot_uniform(output_dim, input_dim),
    glorot_uniform(output_dim, input_dim), glorot_uniform(output_dim, input_dim), glorot_uniform(output_dim, output_dim), glorot_uniform(output_dim, output_dim),
    glorot_uniform(output_dim, output_dim), glorot_uniform(output_dim, output_dim), zeros(output_dim), zeros(output_dim), zeros(output_dim), zeros(output_dim),
    activation_1, activation_2)

function (a :: LSTMLayer)(S,x)

    # assign all opbject variables
    Ug, Uz, Ur, Uh, Wg, Wz, Wr, Wh, bg, bz, br, bh, act1, act2 = a.Ug, a.Uz, a.Ur, a.Uh, a.Wg, a.Wz, a.Wr, a.Wh, a.bg, a.bz, a.br, a.bh, a.act1, a.act2

    # define all operations
    G = act1.(Ug*x + Wg*S .+ bg)
    Z = act1.(Uz*x + Wz*S .+ bz)
    R = act1.(Ur*x + Wr*S .+ br)
    H = act2.(Uh*x + Wh*(S.*R) .+ bh)
    S_new = (1 .- G).*H + Z.*S

    return S_new
end

# define dense layer with activation
mutable struct Dense_start
    W1 :: Array{Float64}
    b1 :: Array{Float64}
    act
end

Dense_start(input_dim::Integer, output_dim::Integer, act) = Dense_start(glorot_uniform(output_dim, input_dim), zeros(output_dim), act)

function (a :: Dense_start)(x)

    # assign all opbject variables
    W1, b1, act = a.W1, a.b1, a.act

    # define all operations
    S_1 = act.(W1*x .+ b1)

    return S_1
end

# define dense layer without activation
mutable struct Dense_end
    W :: Array{Float64}
    b :: Array{Float64}
end

Dense_end(input_dim::Integer, output_dim::Integer) = Dense_end(glorot_uniform(output_dim, input_dim), zeros(output_dim))

function (a :: Dense_end)(x)

    # assign all opbject variables
    W, b = a.W, a.b

    # define all operations
    y = W*x .+ b

    return y
end

# connect them together into one model
mutable struct DGM
    D1 :: Dense_start
    LSTMLayers :: Array{LSTMLayer}
    D2 :: Dense_end
end

DGM(layer_width, n_layers, spatial_dim) = DGM(Dense_start(spatial_dim + 1, layer_width, tanh),
[LSTMLayer(spatial_dim + 1, layer_width, tanh, tanh) for i in 1:n_layers], Dense_end(layer_width, 1))

function (a :: DGM)(t, x)

    # assign object variables
    D1, LSTMLayers, D2 = a.D1, a.LSTMLayers, a.D2

    #concatenate time and space inputs
    X = vcat(t, x)

    # define all operations
    S = D1(X)

    for lstm_layer in LSTMLayers
        S = lstm_layer(S, X)
    end

    y = D2(S)

    return y
end

# convey trainable parameters to Flux
@functor Dense_start
@functor LSTMLayer
@functor Dense_end
@functor DGM

# check if Flux recognizes parameters
m1 = DGM(50, 3, 1)
params(m1)

# define analytical solution function
function analytical_soln(S, K, r, sigma, t)

   d1 = (log.(S./K) .+ (r + sigma^2/2)*(t_term-t))/(sigma*sqrt(t_term-t))
   d2 = d1 .- (sigma*sqrt(t_term-t))
   call_price = S.*cdf.(Normal(0,1), d1) .- K*exp(-r * (t_term-t))*cdf.(Normal(0,1), d2)

   return call_price
end

# define sampling function
function sampling_function(space_internal, space_terminal)

    # within-domain sample
    t_internal = rand(Uniform(t_initial, t_term), space_internal)
    S_internal = rand(Uniform(S_low, S_high*S_multiplier), space_internal)

    # terminal sample
    t_terminal = t_term .*ones(space_terminal)
    S_terminal = rand(Uniform(S_low, S_high*S_multiplier), space_terminal)

    return t_internal', S_internal', t_terminal', S_terminal'
end

# define the loss function
function loss_function()

    t_internal, S_internal, t_terminal, S_terminal = sampling_function(nSim_interior, nSim_terminal)

    # differential operator loss
    model_output_internal = model(t_internal, S_internal)
    ϵ = 0.01

    ∂g∂x = (model(t_internal, S_internal .+ ϵ) - model_output_internal)./ϵ
    ∂g∂t = (model(t_internal .+ ϵ, S_internal) - model_output_internal)./ϵ
    ∂g∂xx = (model(t_internal, S_internal .+ 2*ϵ) - 2*model(t_internal, S_internal .+ ϵ) + model_output_internal)./(ϵ^2)
    operator_loss_vec = ∂g∂t + r.*S_internal.*∂g∂x + (0.5*(sigma^2)).*(S_internal.^2).*∂g∂xx - r.*model_output_internal
    operator_loss = sum(abs2, operator_loss_vec)

    # terminal condition  loss
    target_output_terminal = relu.(S_terminal .- K)
    model_output_terminal = model(t_terminal, S_terminal)
    terminal_loss = sum(abs2, model_output_terminal - target_output_terminal)

    return operator_loss + terminal_loss
end

# initialize all problem parameters
r = 0.05           # Interest rate
sigma = 0.25       # Volatility
K = 50             # Strike
t_term = 1              # Terminal time
S0 = 0.5           # Initial price

# Solution parameters (domain on which to solve PDE)
t_initial = 0 + 1e-10    # time lower bound
S_low = 0.0 + 1e-10  # spot price lower bound
S_high = 2*K         # spot price upper bound

# neural network parameters
num_layers = 3
nodes_per_layer = 50
learning_rate = 0.001

# Training parameters
sampling_stages  = 100   # number of times to resample new time-space domain points
steps_per_sample = 10    # number of SGD steps to take before re-sampling

# Sampling parameters
nSim_interior = 1000
nSim_terminal = 100
S_multiplier  = 1.5   # multiplier for oversampling i.e. draw S from [S_low, S_high * S_multiplier]

# define model
model = DGM(nodes_per_layer, num_layers, 1)
# get params
θ = params(model)

function train_model()
    # set optimizer as ADAM
    evalcb = () -> @show(loss_function())
    opt = ADAM(learning_rate)
    dataset = [(1) for i in 1:20]
    Flux.train!(loss_function, params(model), zip(dataset), opt, cb = throttle(evalcb, 10))
end

@epochs 1000 train_model()

params(model)

# compare with analytical solution
compare_times = [t_initial, t_term/3, 2*t_term/3, t_term-0.1]
n_plots = 41
S_values = [2.5*i for i in 0:40]
S_model = reshape(S_values, 1, length(S_values))

for t in compare_times
    true_values = analytical_soln(S_values, K, r, sigma, t)
    t_model = t.*ones(1, length(S_values))
    fitted_values = model(t_model, S_model)
    fitted_values = reshape(fitted_values, length(fitted_values))
    display(plot(true_values, label = "Truth", linestyle = :dash, lw = 3, title = "t = $(t)",
    xlabel = "Stock prices", ylabel = "Call option prices"))
    display(plot!(fitted_values, label = "DGM", lw = 2))
    savefig("$(t)_euro_call.png")
end


text:
using DataDepsGenerators

const RAW_DATASET_TO_DOI = Dict(
    # {{DATASET}} => https://doi.org/{{DOI}}
    "censusdp1tract" => "10.5281/zenodo.4127049",
    "eia860"  => "10.5281/zenodo.4127027",
    "eia860m" => "10.5281/zenodo.4281337",
    "eia861"  => "10.5281/zenodo.4127029",
    "eia923"  => "10.5281/zenodo.4127040",
    "epacems" => "10.5281/zenodo.4127055",
    "ferc1"   => "10.5281/zenodo.4127044",
    "ferc714" => "10.5281/zenodo.4127101",
)

# TODO: Generate this automatically via Zenodo API?
const RELEASE_VERSION_TO_DOI = Dict(
    v"1.0.0" => "10.5281/zenodo.3653159",
    v"1.1.0" => "10.5281/zenodo.3672068",
)

const LATEST_VERSION = maximum(keys(RELEASE_VERSION_TO_DOI))

const DIR = joinpath(@__DIR__, "datadeps_generated")

function generate_datadeps(dir=DIR)
    ispath(dir) || mkpath(dir)

    for (dataset, doi) in RAW_DATASET_TO_DOI
        filepath = joinpath(dir, "raw_$dataset.jl")
        datadep_str = generate(JSONLD_Web(), "https://doi.org/$doi", "pudl-raw-$dataset")
        open(filepath, "w") do file
            write(file, datadep_str)
        end
    end

    for (version, doi) in RELEASE_VERSION_TO_DOI
        v_str = replace(string(version), '.' => '_')
        filepath = joinpath(dir, "release_v$v_str.jl")
        datadep_str = generate(JSONLD_Web(), "https://doi.org/$doi", "pudl-release-v$version")
        open(filepath, "w") do file
            write(file, datadep_str)
        end

        if version === LATEST_VERSION
            filepath = joinpath(dir, "release_latest.jl")
            open(filepath, "w") do file
                write(file, datadep_str)
            end
        end
    end
    try
        @eval using JuliaFormatter
        format(dir, BlueStyle())
    catch
    end
end


text:
"""
Description

      The uhist() method plots a histogram of a time series in REPL, using the UnicodePlots package.

Usage

      using MarketData
      uhist(cl)
      uhist(ohlc, 4) # the fourth column ("Close") will be plotted

Method Signature(s)

      uhist(t::TimeSeries.TimeArray, n::Int=1)

Details

      UnicodePlots provides simple plots inside of REPL. The uhist() method is a wrapper for the UnicodePlots.histogram()
      method that takes a TimeSeries.TimeArray as the first argument. The second argument defines the column to be
      plotted and defaults to n=1.

References

      https://github.com/Evizero/UnicodePlots.jl

See Also

      uplot() which plots a lineplot of a time series in REPL
"""

function uhist(t::TimeSeries.TimeArray, n::Int=1)
    v = t.values[:,n]
    c = t.colnames[n]
    isa(t.meta, AbstractString) ? m = t.meta :  m = ""
    ttl = string(m, " ", c)
    histogram(v, title=ttl)
end


text:
const rₑ  = 6371220.0 # Radius of the Earth (meters)
const Ωₑ  = 7.292e-5  # Rotational frequency of the Earth (1/second)
const g   = 9.80616   # Gravitational constant (meters/second^2)
function f(x)         # Coriolis term
  2.0*Ωₑ*x[3]/rₑ
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
using Mamba
using Distributions

include("../src/simulation.jl")
include("../gali_bayesian.jl")
include("../hmc/hmc_aux.jl")

yy,shocks = simulate_dsge(GAMMA_0,GAMMA_1,PSI,PI,500)

mm = Model(
    y = Stochastic(1,
    (mu)-> MvNormal(mu,1),
    false
    ),

    mu = Logical(1,
    (par,x) -> dsge_fit(par,x),
    false
    ),

    # alfa = Stochastic(
    # ()-> Beta(3,3),
    # false
    # ),

    bet = Stochastic(
    ()-> Beta(4,2),
    false
    ),

    epsilon = Stochastic(
    () -> Gamma(7,1),
    false
    ),

    sig = Stochastic(
    () -> Gamma(2,1),
    false
    ),

    phi = Stochastic(
    () -> Gamma(2,1),
    false
    ),

    phi_pi = Stochastic(
    () -> Gamma(2,1),
    false
    ),

    phi_y = Stochastic(
    () -> Gamma(2,1),
    false
    ),

    rho_v = Stochastic(
    () -> Beta(1,1),
    false
    ),

    theta = Stochastic(
    () -> Beta(3,3),
    false
    ),

    # s2 = Stochastic(
    # () -> Gamma(2,2),
    # false
    # ),

    par = Logical(1,
    (bet,epsilon,theta,sig,phi,phi_pi,phi_y,rho_v) -> [2/3,bet,epsilon,theta,sig,1,phi,phi_pi,phi_y,rho_v]#() -> [0.3,0.99,6,2/3,2,1,2,1.5,1.5,0.4]#
    )
)

sampling1 = [NUTS([:par])]

sampling2 = [RWM([:par],10)]

dados = Dict{Symbol, Any}(:y => yy[1:499,2], :x => yy[:,2])

inits = [
    Dict(:y => dados[:y],
    #:alfa => rand(Beta(3,3)),
    :bet => rand(Beta(4,2)),
    :epsilon => rand(Gamma(7,1)),
    :sig => rand(Gamma(2,1)),
    :phi => rand(Gamma(2,1)),
    :phi_pi => rand(Gamma(2,1)),
    :phi_y => rand(Gamma(2,1)),
    :rho_v => rand(Beta(1,1)),
    :theta => rand(Beta(3,3)),
    #:s2 => rand(Gamma(2,2))
    )
    for i in 1:2
]

setsamplers!(mm,sampling1)

sim1 = mcmc(mm,dados,inits,5000, burnin=1000, thin=2, chains=2)

describe(sim1)

p = Gadfly.plot(sim1)

Gadfly.draw(p, filename = "plots.svg")

pp = [2/3,inits[1][:bet],inits[1][:epsilon],inits[1][:theta],inits[1][:sig],inits[1][:s2],inits[1][:phi],inits[1][:phi_pi],inits[1][:phi_y],inits[1][:rho_v]]

pp[6] = 0.05

dsge_fit(pp,yy[:,2])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# The content of this file should be upstreamed to Julia proper

function method_instances(@nospecialize(tt::Type), world::UInt=Base.get_world_counter())
    return map(Core.Compiler.specialize_method, method_matches(tt; world))
end

if VERSION >= v"1.7-"
    const hasgenerator = Base.hasgenerator
else
    const hasgenerator = Base.isgenerated
end

function code_lowered_by_type(@nospecialize(tt); generated::Bool=true, debuginfo::Symbol=:default)

    debuginfo = Base.IRShow.debuginfo(debuginfo)
    if debuginfo !== :source && debuginfo !== :none
        throw(ArgumentError("'debuginfo' must be either :source or :none"))
    end
    return map(method_instances(tt)) do m
        if generated && hasgenerator(m)
            if Base.may_invoke_generator(m)
                return ccall(:jl_code_for_staged, Any, (Any,), m)::CodeInfo
            else
                error("Could not expand generator for `@generated` method ", m, ". ",
                      "This can happen if the provided argument types (", t, ") are ",
                      "not leaf types, but the `generated` argument is `true`.")
            end
        end
        code = Base.uncompressed_ir(m.def::Method)
        debuginfo === :none && Base.remove_linenums!(code)
        return code
    end
end

function code_warntype_by_type(io::IO, @nospecialize(tt);
                       debuginfo::Symbol=:default, optimize::Bool=false, kwargs...)
    debuginfo = Base.IRShow.debuginfo(debuginfo)
    lineprinter = Base.IRShow.__debuginfo[debuginfo]
    for (src, rettype) in Base.code_typed_by_type(tt; optimize, kwargs...)
        if !(src isa Core.CodeInfo)
            println(io, src)
            println(io, "  failed to infer")
            continue
        end
        lambda_io::IOContext = io
        p = src.parent
        nargs::Int = 0
        if p isa Core.MethodInstance
            println(io, p)
            print(io, "  from ")
            println(io, p.def)
            p.def isa Method && (nargs = p.def.nargs)
            if !isempty(p.sparam_vals)
                println(io, "Static Parameters")
                sig = p.def.sig
                warn_color = Base.warn_color() # more mild user notification
                for i = 1:length(p.sparam_vals)
                    sig = sig::UnionAll
                    name = sig.var.name
                    val = p.sparam_vals[i]
                    print_highlighted(io::IO, v::String, color::Symbol) =
                        if highlighting[:warntype]
                            Base.printstyled(io, v; color)
                        else
                            Base.print(io, v)
                        end
                    if val isa TypeVar
                        if val.lb === Union{}
                            print(io, "  ", name, " <: ")
                            print_highlighted(io, "$(val.ub)", warn_color)
                        elseif val.ub === Any
                            print(io, "  ", sig.var.name, " >: ")
                            print_highlighted(io, "$(val.lb)", warn_color)
                        else
                            print(io, "  ")
                            print_highlighted(io, "$(val.lb)", warn_color)
                            print(io, " <: ", sig.var.name, " <: ")
                            print_highlighted(io, "$(val.ub)", warn_color)
                        end
                    elseif val isa typeof(Vararg)
                        print(io, "  ", name, "::")
                        print_highlighted(io, "Int", warn_color)
                    else
                        print(io, "  ", sig.var.name, " = ")
                        print_highlighted(io, "$(val)", :cyan) # show the "good" type
                    end
                    println(io)
                    sig = sig.body
                end
            end
        end
        if src.slotnames !== nothing
            slotnames = Base.sourceinfo_slotnames(src)
            lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES => slotnames)
            slottypes = src.slottypes
            nargs > 0 && println(io, "Arguments")
            for i = 1:length(slotnames)
                if i == nargs + 1
                    println(io, "Locals")
                end
                print(io, "  ", slotnames[i])
                if isa(slottypes, Vector{Any})
                    InteractiveUtils.warntype_type_printer(io, slottypes[i], true)
                end
                println(io)
            end
        end
        print(io, "Body")
        InteractiveUtils.warntype_type_printer(io, rettype, true)
        println(io)
@static if VERSION < v"1.7.0"
        Base.IRShow.show_ir(lambda_io, src, lineprinter(src), InteractiveUtils.warntype_type_printer)
else
        irshow_config = Base.IRShow.IRShowConfig(lineprinter(src), InteractiveUtils.warntype_type_printer)
        Base.IRShow.show_ir(lambda_io, src, irshow_config)
end
        println(io)
    end
    nothing
end


text:
using Plots
using Revise
using Wavelets
WT.paul10
supertype(typeof(WT.paul10))
wavelet(WT.morl)
WT.morl
wavelet(WT.morl, s=8)
plotlyjs()
J=11; n = 2^J+395
n = 641
x = testfunction(n,"Bumps")
EmphasizeFrequencyInfo = WT.Morlet(4.7)
y = cwt(x, wavelet(WT.morl))
heatmap(abs.(y)); plot!(x+90,label="")
y = cwt(x, wavelet(WT.dog1))
heatmap(abs.(y)); plot!(x+90,label="")
y = cwt(x, wavelet(WT.paul4))
heatmap(abs.(y)); plot!(x+90,label="")

gr()
pyplot()
waveType = WT.Morlet(6)

Ψ = wavelet(WT.Morlet(5))
daughters, ξ = Wavelets.computeWavelets(n,Ψ)
plot(daughters,legend=false)

Ψ = wavelet(WT.Morlet(2π),decreasing=3.0, averagingLength=2)
daughters, ξ = Wavelets.computeWavelets(n,Ψ)
plot([daughters sum(daughters,dims=2)],legend=false)
size(daughters)

Ψ = wavelet(waveType, s=1,averagingLength=2)
daughters1, ξ = Wavelets.computeWavelets(n,Ψ)
plot([daughters1 sum(daughters1,dims=2)])

Ψ = wavelet(waveType, s=2)
daughters2, ξ = Wavelets.computeWavelets(n,Ψ)
plot([daughters2 sum(daughters2,dims=2)])

Ψ = wavelet(waveType, s=4)
daughters4, ξ = Wavelets.computeWavelets(n,Ψ)
plot([daughters4 sum(daughters4,dims=2)])

waveType = WT.dog6
waveType = WT.paul1

n=321
waveType = WT.Morlet(2π)
Ψ = wavelet(waveType, s=4.0, decreasing=4.0,averagingLength=4)
nOctaves, totalWavelets, sRanges, sWidths =
    WT.getNWavelets(n,Ψ)
totalWavelets
size(cwt(randn(n,5), Ψ, daughters8))
nOctaves
sRanges
nWaveletsInOctave
totalWavelets
daughters8, ω = Wavelets.computeWavelets(n,Ψ)
plt1= plot(abs.([daughters8 sum(daughters8,dims=2)]),legend=false);
#vline!(2 .^ (Ψ.averagingLength:Ψ.averagingLength + nOctaves));
#plt2 = plot(abs.([daughters8 sum(daughters8,dims=2)][1:512,:]),legend=false);
#vline!(min.(2 .^ (Ψ.averagingLength:Ψ.averagingLength + nOctaves), 512));
ψ₁,ψ₂,ϕ₃,σξ = filter_bank(n,computeJ(n),4);
plt3 = plot(abs.([ψ₁ sum(ψ₁,dims=2)])[1:n,:],legend=false)
#plt4 = plot(abs.([ψ₁ sum(ψ₁,dims=2)])[1:512,:],legend=false)
plot(plt1,plt2,plt3, plt4, layout=(4,1))
(totalWavelets, size(ψ₁,2))

Ψ.σ[1] * 2^(6.415+2)
effectiveQualityFactors(sRanges)'
plotSpacing(nOctaves, p, Ψ, totalWavelets)
polySpacing(nOctaves, 4, Ψ.averagingLength, totalWavelets)
nOctaves
effectiveQualityFactors(2 .^ (WT.polySpacing(nOctaves,Ψ.decreasing,Ψ.averagingLength,Ψ)))
WT.polySpacing(nOctaves,Ψ.decreasing,Ψ.averagingLength,Ψ)

p = 3.5; Ψ = wavelet(waveType, s=8.0, decreasing=p, averagingLength=2)
testQSpacing = WT.polySpacing(nOctaves,Ψ);
scatter(testQSpacing, 1:length(testQSpacing), legend=:topleft,label="Q")
testOrigSpacing = polySpacing(nOctaves, p, Ψ.averagingLength, totalWavelets)
scatter!(testOrigSpacing,1:length(testOrigSpacing),label="orig")
println(("Q based", effectiveQualityFactors(2 .^(testQSpacing))',
         length(testQSpacing), minimum(testQSpacing), maximum(testQSpacing)))
println(("Total based", effectiveQualityFactors(2 .^(testOrigSpacing))', length(testOrigSpacing),minimum(testOrigSpacing), maximum(testOrigSpacing)))

scatter!(2 .^ (Ψ.averagingLength:Ψ.averagingLength + nOctaves), fill(.025, Int(Ψ.averagingLength + nOctaves)))

sRanges

function polySpacing(nOct, p, aveLength, totalWavelets)
    samplePoints = 1:totalWavelets
    s =aveLength; T = totalWavelets; O = nOct
    # y= a + b*x^(1/p), solve for a and b with
    # (x₀,y₀)=(1,aveLength)
    # (x₁,y₁)=(totalWavelets, nOct +aveLength)
    a = s - (O) / (T^(1/p) - 1)
    b = O / (T^(1/p) - 1)
    return a .+ b .* (samplePoints).^(1/p)
end

function effectiveQualityFactors(testQSpacing)
    offInt = log2(testQSpacing[end])%1
    effOct = ceil.(Int,log2.(testQSpacing) .- offInt)
    return [count(effOct .==ii) for ii=minimum(effOct):maximum(effOct)]
end
Ψ.σ[1]
Ψ.σ[2]
Ψ.σ[3]
sum((1:size(daughters8,1)) .*daughters8[:,end]./sum(daughters8[:,end]))
sRanges
4π*sRanges[end]/(Ψ.σ[1] + sqrt(2+Ψ.σ[1]^2))

Ψ = wavelet(waveType, s=8.0, decreasing=1.0)
nOctaves, nWaveletsInOctave, totalWavelets, sRanges, sWidths = WT.getNWavelets(n,Ψ)
minimum(sWidths) ./ sWidths
sWidths
totalWavelets
WT.varianceAdjust(Ψ,totalWavelets)
p =8.0
for p=1.0:.25:5.0
    waveType = WT.Morlet(4.5)
    Ψ = wavelet(waveType, s=8.0, decreasing=p)
    nOctaves, nWaveletsInOctave, totalWavelets, sRanges, sWidths =
        WT.getNWavelets(n,Ψ)
    vAdj = WT.varianceAdjust(Ψ,totalWavelets)
    @test vAdj[1] == p
    @test vAdj[end] == 1
    @test length(vAdj) == totalWavelets
end
plot(1 .+a .*(totalWavelets .- (1:21)).^p)
plot(minimum(sWidths) ./ sWidths)
sWidths
sRanges
nWaveletsInOctave
linRange = Array{Array{Float64}}(undef, length(nWaveletsInOctave))
startAt(curOctave,nWavelets) = curOctave > 1 ? 0-1 ./ nWavelets[curOctave-1] : 0
stopAt(curOctave,nWavelets) = curOctave < length(nWavelets) ? 1 +
    1 ./ nWavelets[curOctave+1] : 1
for curOctave = 1:length(nWaveletsInOctave)
    linRange[curOctave] = (range(0,
                                 stop=1, length =
                                nWaveletsInOctave[curOctave]+1) .+ curOctave .+ Ψ.averagingLength .- 1)[2:end]
end
gathered = cat(linRange..., dims=1)

justUse = 1:21#[1,3,5,7, 9:21...]
p = 4; plot(scatter(gathered, 1:length(gathered),xaxis = "log freq", yaxis="index", legend=false), scatter(2 .^ gathered,
                                                                                1:length(gathered), legend=false),
            scatter((5 .+(1.0:21)).^(1/p),
                    1.0:21,legend=false,title="proposed, 'linear'"),scatter(2 .^((5 .+(1.0:21)).^(1/p))[justUse],
                                                                            1:length(justUse),legend=false,title="proposed
                                                                     exp"),
            layout=4)

p = 8; plot(scatter(gathered, 1:length(gathered),xaxis = "log freq", yaxis="index", legend=false), scatter(2 .^ gathered,
                                                                                1:length(gathered), legend=false),
            scatter(WT.polySpacing(nOctaves,p,Ψ.averagingLength, totalWavelets),
                    1.0:totalWavelets,legend=false,xaxis="log freq"),
            scatter(2 .^(WT.polySpacing(nOctaves,p,Ψ.averagingLength, totalWavelets)),
                    1:totalWavelets,legend=false,xaxis="freq"),
            layout=4)
Ψ.averagingLength
function plotSpacing(nOctaves,Ψ,totalWavelets)
    xVals = polySpacing(nOctaves,Ψ.decreasing,Ψ.averagingLength, totalWavelets)
    plt1 = scatter(xVals, 1.0:totalWavelets,legend=false, xaxis="log freq")
    vline!(Array( Ψ.averagingLength: (Ψ.averagingLength + nOctaves)))
    plt2 = scatter(Ψ.σ[1]*2 .^(xVals), 1:totalWavelets, legend=false, xaxis="freq")
    #plot(2.^(1:totalWavelets), 1:totalWavelets)
    #vline!(2 .^ (Ψ.averagingLength:Ψ.averagingLength+totalWavelets))
    plot(plt1, plt2,layout=2)
end

nOctaves + Ψ.averagingLength+1 -1
sp = WT.polySpacing(nOctaves,4,Ψ.averagingLength, totalWavelets)
Δ
Δ = diff(2 .^ sp)./2
[Δ[1]; Δ]
2 .^ sp
scatter(2 .^(sp),1:length(sp), xerror=[Δ[1]; Δ])
gathered
scatter(WT.polySpacing(nOctaves,4,Ψ.averagingLength, totalWavelets),legend=false)
nOctaves

spaced(nPoints,)
sPoints(nO,Q,β,Ψ) = Q .* ((1.0*nO .-
                           (Ψ.averagingLength:(nO-1))).^(-β+1))
sPoints(11,8,.9,Ψ)
length(4:10)
Ψ.averagingLength:(11-1)

function getScales(n1::Integer, c::CFW)
    nOctaves = log2(max(n1, 2)) - c.averagingLength
    nWaveletsInOctave = reverse([max(1, round(Int, c.scalingFactor /
                                              x^(c.decreasing))) for
                                 x=1:round(Int, nOctaves)])
    sj = zeros(sum(nWaveletsInOctave))
    for curOctave = 1:round(Int, nOctaves)
        sRange = (2 .^ (range(0, stop=1, length = nWaveletsInOctave[curOctave]+1) .+
                        curOctave .+ c.averagingLength .- 1))[1:end-1]
        ind = (sum(nWaveletsInOctave[1:curOctave-1])+1):sum(nWaveletsInOctave[1:curOctave])
        @debug "ind = $ind"
        sj[ind] = sRange
    end
    return sj
end


text:
# NOTE: For convenience, these imports are written in a separate file so that
# the main node and worker nodes can source this file in order to import
# libraries. Libraries that are only used in the worker nodes still need to be
# imported in the main node to precompile stale Julia cache files.
# If this doesn't occur, each worker node will try to write these cache
# files simultaneously, causing strage errors.
using Cytof5

import Pkg
Pkg.activate("../../")  # sims

using Random
using Distributions
using DelimitedFiles
import DataFrames
const DF = DataFrames
import CSV
import PyPlot, PyCall, Seaborn
const plt = PyPlot.plt
const sns = Seaborn
using BSON


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

@testset "test_Vec2D" begin
     v1 = Vec2d(5.2 + 1e-11, 6.3)
     v2 = Vec2d(5.2, 6.3 - (2*1e-11))
     @test v1 ≈ v2
end

@testset "test_World" begin
     w = World()
	sph1 = Sphere(translation(VEC_X * 2))
	sph2 = Sphere(translation(VEC_X * 8))
	add_shape!(w, sph1)
	add_shape!(w, sph2)

	intersection1 = ray_intersection(w, Ray(Point(0.0, 0.0, 0.0), VEC_X))
	@test intersection1.world_point ≈ Point(1., 0., 0.)

	intersection2 = ray_intersection(w, Ray(Point(10.0, 0.0, 0.0), -VEC_X))
	@test intersection2.world_point ≈ Point(9., 0., 0.)
end


text:
# ps3000_run_streaming.jl

export ps3000_run_streaming

function ps3000_run_streaming(handle :: Int16, sample_interval_ms, max_samples, windowed)
  ps_status = ccall((:ps3000_run_streaming, ps3000driver), Int16,
    (Int16, Int16, Int32, Int16), handle, sample_interval_ms, max_samples, windowed)
  if ps_status == 0
    error("ps3000 Error: problem occurred or a value out of range")
  end
  return nothing
end

function ps3000_run_streaming(handle :: Int16, sample_interval_ms, max_samples, windowed::Bool)
  ps3000_run_streaming(handle,sample_interval_ms,max_samples,windowed?0x0001:0x0000)
end


text:
using Test, CodeTools
import Base.Docs: Binding

@test CodeTools.getthing("Base.sin") == sin
@test CodeTools.getthing(Base, [:sin]) == sin

@testset "documentation" begin
  let
    function func(x)
      infunc = (y) -> y^x
      return infunc
    end

    @test CodeTools.hasdoc(func(3)) == false
  end

  let
    @doc """
      docfunc(x)
  docs

  more docs which aren't in the description
    """
    function docfunc(x) end

    @test CodeTools.hasdoc(docfunc) == true
    @test chomp(CodeTools.description(Binding(Main, :docfunc))) == "docs"
    @test CodeTools.signature(Binding(Main, :docfunc)) == "docfunc(x)"
    @test CodeTools.completiontype(docfunc) == "λ"
  end

  let
    @doc """
  i'm a constant

  with multiline docs
    """
    foo = :bar

    @test CodeTools.hasdoc(foo) == true
    @test chomp(CodeTools.description(Binding(Main, :foo))) == "i'm a constant"
    @test CodeTools.signature(Binding(Main, :foo)) == nothing
    @test CodeTools.completiontype(foo) == "constant"
  end
end

# module detection tests
@testset "module detection" begin
  code = ["""
  module Mod1
  [x for x=1:2]
  end
  1+1
  """,
  """
  module Mod2
  module Foo
  # for
  end
  1+1
  """,
  """
  module Mod3
  abstract foo
  end
  """,
  """
  module Mod4
  bitstype 8 foo
  (i for i = 1:10)
  end
  """,
  """
  :module Foo
  1+1
  end
  """,
  """
  bla.module
  foo
  end
  """,
  """
  begin
  module Foo
  2+2
  end
  1+1
  """
  ]

  @test CodeTools.codemodule(code[1], 2) == "Mod1"
  @test CodeTools.codemodule(code[1], 4) == ""
  @test CodeTools.codemodule(code[2], 3) == "Mod2.Foo"
  @test CodeTools.codemodule(code[2], 5) == "Mod2"
  @test CodeTools.codemodule(code[3], 2) == "Mod3"
  @test CodeTools.codemodule(code[3], 3) == ""
  @test CodeTools.codemodule(code[4], 2) == "Mod4"
  @test CodeTools.codemodule(code[4], 3) == "Mod4"
  @test CodeTools.codemodule(code[4], 4) == ""
  @test CodeTools.codemodule(code[5], 1) == ""
  @test CodeTools.codemodule(code[5], 2) == ""
  @test CodeTools.codemodule(code[5], 3) == ""
  @test CodeTools.codemodule(code[6], 1) == ""
  @test CodeTools.codemodule(code[6], 2) == ""
  @test CodeTools.codemodule(code[6], 3) == ""
  @test CodeTools.codemodule(code[7], 1) == ""
  @test CodeTools.codemodule(code[7], 3) == "Foo"
  @test CodeTools.codemodule(code[7], 5) == ""

  @test CodeTools.filemodule(normpath(joinpath(@__DIR__, "..", "src", "module.jl"))) == "CodeTools"

  @test CodeTools.includeline(normpath(joinpath(@__DIR__, "..", "src", "CodeTools.jl")), normpath(joinpath(@__DIR__, "..", "src", "utils.jl"))) == 5
end

@testset "completions" begin
  @test CodeTools.prefix("MacroTools.@") == ["MacroTools", "@"]
  @test length(CodeTools.completions("CodeTools.pre")) > 0
  @test "LinearAlgebra" in CodeTools.stdlibs()
end


text:
@testset "ForwardDiff comparison" begin
    Random.seed!(100)
    x0 = srand(3)
    v0 = srand(3)
    q0 = rand(QuatRotation).q
    ϕ0 = srand(3)
    timestep0 = 0.01
    x1 = Dojo.next_position(x0, v0, timestep0)
    q1 = Dojo.next_orientation(q0, ϕ0, timestep0)

    # ∇v
    FD∂i∂v = ForwardDiff.jacobian(v0 -> Dojo.next_position(x0, v0, timestep0), v0)
    @test norm(FD∂i∂v - Dojo.linear_integrator_jacobian_velocity(x0, v0, timestep0), Inf) < 1.0e-8
    # ∇ϕ
    FD∂i∂ϕ = ForwardDiff.jacobian(ϕ0 -> Dojo.vector(Dojo.next_orientation(q0, ϕ0, timestep0)), ϕ0)
    @test norm(FD∂i∂ϕ - Dojo.rotational_integrator_jacobian_velocity(q0, ϕ0, timestep0), Inf) < 1.0e-8

    # ∇x
    FD∂i∂x = ForwardDiff.jacobian(x0 -> Dojo.next_position(x0, v0, timestep0), x0)
    @test norm(FD∂i∂x - Dojo.linear_integrator_jacobian_position(x0, v0, timestep0), Inf) < 1.0e-8
    # ∇q
    FD∂i∂q = ForwardDiff.jacobian(q0 ->
        Dojo.vector(Dojo.next_orientation(Quaternion(q0...), ϕ0, timestep0)),
        Dojo.vector(q0))
    @test norm(FD∂i∂q - Dojo.rotational_integrator_jacobian_orientation(q0, ϕ0, timestep0, attjac=false), Inf) < 1.0e-8

    # ∇vϕ
    vel0 = [v0; ϕ0]
    FDvel = ForwardDiff.jacobian(vel0 ->
        [Dojo.next_position(x0, SVector{3}(vel0[1:3]), timestep0);
        Dojo.vector(Dojo.next_orientation(q0, SVector{3}(vel0[4:6]), timestep0))], vel0)
    @test norm(FDvel - Dojo.integrator_jacobian_velocity(x0, v0, q0, ϕ0, timestep0), Inf) < 1.0e-8

    # ∇xq
    con0 = [x0; Dojo.vector(q0)]
    FDcon = ForwardDiff.jacobian(con0 ->
        [Dojo.next_position(SVector{3}(con0[1:3]), v0, timestep0);
        Dojo.vector(Dojo.next_orientation(Quaternion(con0[4:7]...), ϕ0, timestep0))], con0)
    @test norm(FDcon - Dojo.integrator_jacobian_configuration(x0, v0, q0, ϕ0, timestep0, attjac=false), Inf) < 1.0e-8
    @test norm(FDcon * cat(I(3), Dojo.LVᵀmat(q0), dims=(1,2)) - Dojo.integrator_jacobian_configuration(x0, v0, q0, ϕ0, timestep0, attjac=true), Inf) < 1.0e-8
end


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
using MAT

include("nnutil.jl")
include("NNTrussPull_Domain.jl")

nnname = "Data/$(nntype)/trained_from$(start_id)_ite$(i).mat"
nnname = "Data/$(nntype)/learned_nn_ite$(i).mat"

s = ae_to_code(nnname, nntype)
eval(Meta.parse(s))

# testtype = "PlaneStressPlasticity"
testtype = "NeuralNetwork1D"

prop = Dict("name"=> testtype, "rho"=> 0.1, "E"=> 200, "B"=> 10.0,
            "sigmaY"=>0.300, "K"=>1/9*200, "A0"=> 1.0, "eta"=> 10.0, "nn"=>post_nn)

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)
# need to update state in domain from globdat
updateStates!(domain, globdat)

for i = 1:NT
    @info i, "/" , NT
    solver = NewmarkSolver(Δt, globdat, domain, -1.0, 0.0, 1e-5, 1e-5, 100) # ok
end
@info tid
domain_te = domain
if learn
    @save "Data/learn_domain$(tid)_te.jld2" domain_te
else
    @save "Data/domain$(tid)_te.jld2" domain_te
end


text:
using ZeroInflatedDistributions
using Distributions
using QuadGK
using Random
using Test

@testset "ZeroInflatedDistributions.jl" begin
    @testset "Link functions" begin
        p1 = 0
        p2 = 1
        # Logit-log link tests
        ll = LogitLogLink()
        @test encprob(ll, p1, p2) == 1//2
        @test posrate(ll, p1, p2) == Float64(ℯ)
        # Poisson link tests
        pl = PoissonLink()
        @test isone(pl.offset)
        @test_throws DomainError PoissonLink(-1)
        @test encprob(pl, p1, p2) == 0.6321205588285577
        @test posrate(pl, p1, p2) == 4.300258535328371
        # Identity link tests
        il = IdentityLink()
        @test_throws DomainError encprob(il, -1)
        @test_throws DomainError encprob(il, 1.1)
        @test encprob(il, 0.5) == 0.5
        @test posrate(il, 1.0) == 1
    end

    @testset "Zero-inflated distribution constructors" begin
        p1 = 0
        p2 = 1
        disp = 0.5
        ll = LogitLogLink()
        pl = PoissonLink()
        # Test inner constructor
        zil0 = ZeroInflatedDistribution(Bernoulli(0.5), LogNormal(0.0, 1.0))
        @test logpdf(zil0, 0) == log(0.5)
        @test logpdf(zil0, 1) == logpdf(LogNormal(0.0, 1.0), 1) + log(0.5)

        # Without bias correction, posrate is the median of the log-normal
        zil1 = ZeroInflatedDistribution(ll, LogNormal, p1, p2, disp; biascorrect = false)
        @test median(zil1.posdist) == posrate(ll, p1, p2)

        # With bias correction (default) it is the mean
        zil2 = ZeroInflatedDistribution(ll, LogNormal, p1, p2, disp)
        @test mean(zil2.posdist) == posrate(ll, p1, p2)

        zil3 = ZeroInflatedDistribution(pl, Gamma, p1, p2, disp)
        @test mean(zil3.posdist) == posrate(pl, p1, p2)
        @test std(zil3.posdist) == disp

        zil4 = ZeroInflatedDistribution(pl, InverseGamma, p1, p2, disp)
        @test mean(zil4.posdist) == posrate(pl, p1, p2)
        @test std(zil4.posdist) == disp

        zil5 = ZeroInflatedDistribution(pl, InverseGaussian, p1, p2, disp)
        @test mean(zil5.posdist) == posrate(pl, p1, p2)
        @test shape(zil5.posdist) == disp
    end

    @testset "Zero-inflated log-pdfs" begin
        zil = ZeroInflatedDistribution(Bernoulli(0.5), LogNormal(0.0, 1.0))
        @test logpdf(zil, 0) == log(0.5)
        @test logpdf(zil, 1) == logpdf(LogNormal(0.0, 1.0), 1) + log(0.5)
    end

    @testset "Random zero-inflated data generation" begin
        n = 1_000
        zil = ZeroInflatedDistribution(Bernoulli(0.5), LogNormal(-1/2, 1.0))
        x = rand(zil, n)
        @test any(x .== 0)
        @test all(x .≥ 0)
    end

    @testset "Zero-inflated distribution statistics" begin
        zil = ZeroInflatedDistribution(Bernoulli(0.75), LogNormal(-1/2, 1.0))
        zil0 = ZeroInflatedDistribution(Bernoulli(0), LogNormal(-1/2, 1.0))
        zil1 = ZeroInflatedDistribution(Bernoulli(1), LogNormal(-1/2, 1.0))
        zilgam = ZeroInflatedDistribution(Bernoulli(0.85), Gamma(2, 4))

        @testset "Zero-inflated distribution supports" begin
            @test !insupport(zil, -1)
            @test all(minimum.([zil, zil0, zil1, zilgam]) .== 0)
            @test all(maximum.([zil, zil0, zil1, zilgam]) .== Inf)
        end

        @testset "Test means and variances via simulation" begin
            Random.seed!(13579)
            n = 10_000

            x = rand(zil, n)
            x0 = rand(zil0, n)
            x1 = rand(zil1, n)

            sd = std(zil) / sqrt(n)
            sd1 = std(zil1) / sqrt(n)

            @test isapprox(mean(x), 0.75, atol = 2 * sd)
            @test all(x0 .== 0)
            @test isapprox(mean(x1), 1, atol = 2 * sd1)

            # @test isapprox(var(zil), var(x))
            @test var(zil0) == 0
            @test var(zil1) == var(zil1.posdist)
        end

        @testset "Test means and variances via quadrature" begin
            # 𝔼[X]
            quadex = quadgk(x -> x * pdf(zil, x), 0, Inf)
            quadexgam = quadgk(x -> x * pdf(zilgam, x), 0, Inf)
            # 𝔼[X²]
            quadex2 = quadgk(x -> x^2 * pdf(zil, x), 0, Inf)
            quadexgam2 = quadgk(x -> x^2 * pdf(zilgam, x), 0, Inf)
            # 𝕍[X] = 𝔼[X^2] - 𝔼[X]^2
            quadvar = quadex2[1] - mean(zil)^2
            quadvargam = quadexgam2[1] - mean(zilgam)^2

            @test isapprox(quadex[1], mean(zil), atol = quadex[2])
            @test isapprox(quadvar, var(zil), atol = quadex2[2])
            @test isapprox(quadexgam[1], mean(zilgam), atol = quadexgam[2])
            @test isapprox(quadvargam, var(zilgam), atol = quadexgam2[2])
        end

        @testset "Test CDF function" begin
            @test cdf(zil, 0) == failprob(zil.encdist)
            @test cdf(zil0, 0) == 1
            @test cdf(zil1, 2) == cdf(zil1.posdist, 2)
            @test cdf(zilgam, 10) ==
                failprob(zilgam.encdist) +
                cdf(zilgam.posdist, 10) * succprob(zilgam.encdist)
        end

        @testset "Test quantile function" begin
            @test quantile(zil, 0.2) == 0
            @test quantile(zilgam, 0.1) == 0
            @test quantile(zil, 0.5) ==
                quantile(zil.posdist,
                         (0.5 - failprob(zil.encdist)) / succprob(zil.encdist))
            @test quantile(zil1, 0.9) == quantile(zil1.posdist, 0.9)
            @test quantile(zilgam, 0.75) ==
                quantile(zilgam.posdist,
                         (0.75 - failprob(zilgam.encdist)) /
                         succprob(zilgam.encdist))
        end

        @testset "Test modes function" begin
            @test all(modes(zil) .== [0, mode(LogNormal(-0.5, 1))])
            @test all(modes(zilgam) .== [0, mode(Gamma(2, 4))])
        end
    end
end


text:
"""
Type: NLPStopping

Methods: `start!`, `stop!`, `update_and_start!`, `update_and_stop!`, `fill_in!`, `reinit!`, `status`,
`KKT`, `unconstrained_check`, `unconstrained2nd_check`, `optim_check_bounded`

Specialization of `GenericStopping`. Stopping structure for non-linear optimization models using `NLPModels` ( https://github.com/JuliaSmoothOptimizers/NLPModels.jl ).

Attributes:
- `pb`         : An `AbstractNLPModel`.
- `current_state`      : The information relative to the problem, see `GenericState` or `NLPAtX`.
- (opt) `meta` : Metadata relative to stopping criteria, see `StoppingMeta`.
- (opt) `main_stp` : Stopping of the main loop in case we consider a Stopping
                          of a subproblem.
                          If not a subproblem, then `VoidStopping`.
- (opt) `listofstates` : ListofStates designed to store the history of States.
- (opt) `stopping_user_struct` : Contains any structure designed by the user.

Constructors:
- `NLPStopping(pb::AbstractNLPModel, meta::AbstractStoppingMeta, stop_remote::AbstractStopRemoteControl, state::AbstractState; main_stp::AbstractStopping=VoidStopping(), list::AbstractListofStates = VoidListofStates(), user_struct::AbstractDict = Dict(), kwargs...)`
     The default constructor.
- `NLPStopping(pb::AbstractNLPModel, meta::AbstractStoppingMeta, state::AbstractState; main_stp::AbstractStopping=VoidStopping(), list::AbstractListofStates = VoidListofStates(), user_struct::AbstractDict = Dict(), kwargs...)`
     The one passing the `kwargs` to the `stop_remote`.
- `GenericStopping(pb::AbstractNLPModel, state::AbstractState; stop_remote::AbstractStopRemoteControl = StopRemoteControl(), main_stp::AbstractStopping=VoidStopping(), list::AbstractListofStates = VoidListofStates(), user_struct::AbstractDict = Dict(), kwargs...)`
     The one passing the `kwargs` to the `meta`.
- `GenericStopping(pb::AbstractNLPModel; n_listofstates=, kwargs...)`
     The one setting up a default state `NLPAtX` using `pb.meta.x0`, and initializing the list of states if `n_listofstates>0`. The optimality function is the function `KKT` unless `optimality_check` is in the `kwargs`.

 Notes:
- Designed for `NLPAtX` State. Constructor checks that the State has the required entries.

 """
mutable struct NLPStopping{Pb, M, SRC, T, MStp, LoS} <: AbstractStopping{Pb, M, SRC, T, MStp, LoS}

  # problem
  pb::Pb

  # Common parameters
  meta::M
  stop_remote::SRC

  # current state of the problem
  current_state::T

  # Stopping of the main problem, or nothing
  main_stp::MStp

  # History of states
  listofstates::LoS

  # User-specific structure
  stopping_user_struct::AbstractDict
end

function NLPStopping(
  pb::Pb,
  meta::M,
  stop_remote::SRC,
  current_state::T;
  main_stp::AbstractStopping = VoidStopping(),
  list::AbstractListofStates = VoidListofStates(),
  user_struct::AbstractDict = Dict(),
  kwargs...,
) where {
  Pb <: AbstractNLPModel,
  M <: AbstractStoppingMeta,
  SRC <: AbstractStopRemoteControl,
  T <: AbstractState,
}
  return NLPStopping(pb, meta, stop_remote, current_state, main_stp, list, user_struct)
end

function NLPStopping(
  pb::Pb,
  meta::M,
  current_state::T;
  main_stp::AbstractStopping = VoidStopping(),
  list::AbstractListofStates = VoidListofStates(),
  user_struct::AbstractDict = Dict(),
  kwargs...,
) where {Pb <: AbstractNLPModel, M <: AbstractStoppingMeta, T <: AbstractState}
  stop_remote = StopRemoteControl(; kwargs...) #main_stp == VoidStopping() ? StopRemoteControl() : cheap_stop_remote_control()

  return NLPStopping(pb, meta, stop_remote, current_state, main_stp, list, user_struct)
end

function NLPStopping(
  pb::Pb,
  current_state::T;
  stop_remote::AbstractStopRemoteControl = StopRemoteControl(),
  main_stp::AbstractStopping = VoidStopping(),
  list::AbstractListofStates = VoidListofStates(),
  user_struct::AbstractDict = Dict(),
  kwargs...,
) where {Pb <: AbstractNLPModel, T <: AbstractState}
  if :max_cntrs in keys(kwargs)
    mcntrs = kwargs[:max_cntrs]
  else
    mcntrs = init_max_counters()
  end

  if :optimality_check in keys(kwargs)
    oc = kwargs[:optimality_check]
  else
    oc = KKT
  end

  meta = StoppingMeta(; max_cntrs = mcntrs, optimality_check = oc, kwargs...)

  return NLPStopping(pb, meta, stop_remote, current_state, main_stp, list, user_struct)
end

function NLPStopping(pb::AbstractNLPModel; n_listofstates::Integer = 0, kwargs...)
  #Create a default NLPAtX
  nlp_at_x = NLPAtX(pb.meta.x0)

  if n_listofstates > 0 && :list ∉ keys(kwargs)
    list = ListofStates(n_listofstates, Val{typeof(nlp_at_x)}())
    return NLPStopping(pb, nlp_at_x, list = list, optimality_check = KKT; kwargs...)
  end

  return NLPStopping(pb, nlp_at_x, optimality_check = KKT; kwargs...)
end

"""
init\\_max\\_counters:
initialize the maximum number of evaluations on each of
the functions present in the NLPModels.Counters, e.g.

`init_max_counters(; allevals :: T = typemax(T), obj = allevals, grad = allevals, cons = allevals, jcon = allevals, jgrad = allevals, jac = allevals, jprod = allevals, jtprod = allevals, hess = allevals, hprod = allevals, jhprod = allevals, sum = 11 * allevals, kwargs...)`

`:neval_sum` is by default limited to `|Counters| * allevals`.
"""
function init_max_counters(; allevals::T = typemax(Int), kwargs...) where {T <: Integer}
  entries = [Meta.parse(split("$(f)", '_')[2]) for f in fieldnames(Counters)]
  lim_fields = keys(kwargs)
  cntrs = Dict{Symbol, T}([
    (Meta.parse("neval_$(t)"), t in lim_fields ? kwargs[t] : allevals) for t in entries
  ])
  push!(cntrs, (:neval_sum => :sum in lim_fields ? kwargs[:sum] : typemax(T)))

  return cntrs
end

function max_evals!(stp::NLPStopping, allevals::Integer)
  stp.meta.max_cntrs = init_max_counters(allevals = allevals)
  return stp
end

function max_evals!(stp::NLPStopping; allevals::T = typemax(Int), kwargs...) where {T <: Integer}
  stp.meta.max_cntrs = init_max_counters(allevals = allevals; kwargs...)
  return stp
end

"""
init\\_max\\_counters\\_NLS:
initialize the maximum number of evaluations on each of
the functions present in the `NLPModels.NLSCounters`, e.g.

`init_max_counters_NLS(; allevals = typemax(T), residual = allevals, jac_residual = allevals, jprod_residual = allevals, jtprod_residual = allevals, hess_residual = allevals, jhess_residual = allevals, hprod_residual = allevals, kwargs...)`
"""
function init_max_counters_NLS(; allevals::T = typemax(Int), kwargs...) where {T <: Integer}
  cntrs_nlp = init_max_counters(; allevals = allevals, kwargs...)

  entries =
    [Meta.parse(split("$(f)", '_')[2]) for f in setdiff(fieldnames(NLSCounters), [:counters])]
  lim_fields = keys(kwargs)
  cntrs = Dict{Symbol, T}([
    (Meta.parse("neval_$(t)"), t in lim_fields ? kwargs[t] : allevals) for t in entries
  ])

  return merge(cntrs_nlp, cntrs)
end

"""
fill_in!: (NLPStopping version) a function that fill in the required values in the `NLPAtX`.

`fill_in!( :: NLPStopping, :: Union{AbstractVector, Nothing}; fx :: Union{AbstractVector, Nothing} = nothing, gx :: Union{AbstractVector, Nothing} = nothing, Hx :: Union{MatrixType, Nothing} = nothing, cx :: Union{AbstractVector, Nothing} = nothing, Jx :: Union{MatrixType, Nothing} = nothing, lambda :: Union{AbstractVector, Nothing} = nothing, mu :: Union{AbstractVector, Nothing} = nothing, matrix_info :: Bool = true, kwargs...)`
"""
function fill_in!(
  stp::NLPStopping{Pb, M, SRC, NLPAtX{S, T, HT, JT}, MStp, LoS},
  x::T;
  fx::Union{eltype(T), Nothing} = nothing,
  gx::Union{T, Nothing} = nothing,
  Hx::Union{HT, Nothing} = nothing,
  cx::Union{T, Nothing} = nothing,
  Jx::Union{JT, Nothing} = nothing,
  lambda::Union{T, Nothing} = nothing,
  mu::Union{T, Nothing} = nothing,
  matrix_info::Bool = true,
  kwargs...,
) where {Pb, M, SRC, MStp, LoS, S, T, HT, JT}
  gfx = isnothing(fx) ? obj(stp.pb, x) : fx
  ggx = isnothing(gx) ? grad(stp.pb, x) : gx

  if isnothing(Hx) && matrix_info
    gHx = hess(stp.pb, x).data
  else
    gHx = Hx
  end

  if stp.pb.meta.ncon > 0
    gJx = isnothing(Jx) ? jac(stp.pb, x) : Jx
    gcx = isnothing(cx) ? cons(stp.pb, x) : cx
  else
    gJx = stp.current_state.Jx
    gcx = stp.current_state.cx
  end

  #update the Lagrange multiplier if one of the 2 is asked
  if (stp.pb.meta.ncon > 0 || has_bounds(stp.pb)) && (isnothing(lambda) || isnothing(mu))
    lb, lc = _compute_mutliplier(stp.pb, x, ggx, gcx, gJx; kwargs...)
  elseif stp.pb.meta.ncon == 0 && !has_bounds(stp.pb) && isnothing(lambda)
    lb, lc = mu, stp.current_state.lambda
  else
    lb, lc = mu, lambda
  end

  return update!(stp, x = x, fx = gfx, gx = ggx, Hx = gHx, cx = gcx, Jx = gJx, mu = lb, lambda = lc)
end

function fill_in!(
  stp::NLPStopping{Pb, M, SRC, OneDAtX{S, T}, MStp, LoS},
  x::T;
  fx::Union{T, Nothing} = nothing,
  gx::Union{T, Nothing} = nothing,
  f₀::Union{T, Nothing} = nothing,
  g₀::Union{T, Nothing} = nothing,
  kwargs...,
) where {Pb, M, SRC, S, T, MStp, LoS}
  gfx = isnothing(fx) ? obj(stp.pb, x) : fx
  ggx = isnothing(gx) ? grad(stp.pb, x) : gx
  gf₀ = isnothing(f₀) ? obj(stp.pb, 0.0) : f₀
  gg₀ = isnothing(g₀) ? grad(stp.pb, 0.0) : g₀

  return update!(stp.current_state, x = x, fx = gfx, gx = ggx, f₀ = gf₀, g₀ = gg₀)
end

"""
For NLPStopping, `rcounters` set as true also reinitialize the counters.
"""
function reinit!(
  stp::NLPStopping;
  rstate::Bool = false,
  rlist::Bool = true,
  rcounters::Bool = false,
  kwargs...,
)
  stp.meta.start_time = NaN
  stp.meta.optimality0 = 1.0

  #reinitialize the boolean status
  reinit!(stp.meta)

  #reinitialize the counter of stop
  stp.meta.nb_of_stop = 0

  #reinitialize the list of states
  if rlist && (typeof(stp.listofstates) != VoidListofStates)
    #TODO: Warning we cannot change the type of ListofStates
    stp.listofstates = rstate ? VoidListofStates() : ListofStates(stp.current_state)
  end

  #reinitialize the state
  if rstate
    reinit!(stp.current_state; kwargs...)
  end

  #reinitialize the NLPModel Counters
  if rcounters && typeof(stp.pb) <: AbstractNLPModel
    NLPModels.reset!(stp.pb)
  end

  return stp
end

"""
\\_resources\\_check!: check if the optimization algorithm has exhausted the resources.
                   This is the NLP specialized version that takes into account
                   the evaluation of the functions following the sum_counters
                   structure from NLPModels.

`_resources_check!(:: NLPStopping, :: T)`

Note:
- function uses counters in `stp.pb`, and update the counters in the state.
- function is compatible with `Counters`, `NLSCounters`, and any type whose entries match the entries of `max_cntrs`.
- all the NLPModels have an attribute `counters` and a function `sum_counters(nlp)`.
"""
function _resources_check!(stp::NLPStopping, x::T) where {T <: Union{AbstractVector, Number}}
  max_cntrs = stp.meta.max_cntrs

  if max_cntrs == Dict{Symbol, Int64}()
    return stp.meta.resources
  end

  # check all the entries in the counter
  max_f = false
  if typeof(stp.pb) <: AbstractNLSModel ||
     (:counters in fieldnames(typeof(stp.pb)) && typeof(stp.pb.counters) == NLSCounters)
    for f in intersect(fieldnames(NLSCounters), keys(max_cntrs))
      max_f = f != :counters ? (max_f || (eval(f)(stp.pb) > max_cntrs[f])) : max_f
    end
    for f in intersect(fieldnames(Counters), keys(max_cntrs))
      max_f = max_f || (eval(f)(stp.pb) > max_cntrs[f])
    end
  elseif typeof(stp.pb) <: AbstractNLPModel ||
         (:counters in fieldnames(typeof(stp.pb)) && typeof(stp.pb.counters) == Counters)
    for f in intersect(fieldnames(Counters), keys(max_cntrs))
      max_f = max_f || (eval(f)(stp.pb) > max_cntrs[f])
    end
  end

  # Maximum number of function and derivative(s) computation
  if :neval_sum in keys(stp.meta.max_cntrs)
    max_evals = sum_counters(stp.pb) > max_cntrs[:neval_sum]
  end

  # global user limit diagnostic
  if (max_evals || max_f)
    stp.meta.resources = true
  end

  return stp.meta.resources
end

"""
    `_unbounded_problem_check!`: This is the NLP specialized version that takes into account
                   that the problem might be unbounded if the objective or the
                   constraint function are unbounded.

    `_unbounded_problem_check!(:: NLPStopping, :: AbstractVector)`

Note:
- evaluate the objective function if `state.fx` for NLPAtX or `state.fx` for OneDAtX is `_init_field` and store in `state`.
- if minimize problem (i.e. nlp.meta.minimize is true) check if `state.fx <= - meta.unbounded_threshold`, otherwise check `state.fx ≥ meta.unbounded_threshold`.
"""
function _unbounded_problem_check!(
  stp::NLPStopping{Pb, M, SRC, NLPAtX{S, T, HT, JT}, MStp, LoS},
  x::AbstractVector,
) where {Pb, M, SRC, MStp, LoS, S, T, HT, JT}
  if isnan(stp.current_state.fx)
    stp.current_state.fx = obj(stp.pb, x)
  end

  if stp.pb.meta.minimize
    f_too_large = stp.current_state.fx <= -stp.meta.unbounded_threshold
  else
    f_too_large = stp.current_state.fx >= stp.meta.unbounded_threshold
  end

  if f_too_large
    stp.meta.unbounded_pb = true
  end

  return stp.meta.unbounded_pb
end

function _unbounded_problem_check!(
  stp::NLPStopping{Pb, M, SRC, OneDAtX{S, T}, MStp, LoS},
  x::Union{AbstractVector, Number},
) where {Pb, M, SRC, MStp, LoS, S, T}
  if isnan(stp.current_state.fx)
    stp.current_state.fx = obj(stp.pb, x)
  end

  if stp.pb.meta.minimize
    f_too_large = stp.current_state.fx <= -stp.meta.unbounded_threshold
  else
    f_too_large = stp.current_state.fx >= stp.meta.unbounded_threshold
  end

  return stp.meta.unbounded_pb
end
"""
\\_infeasibility\\_check!: This is the NLP specialized version.

Note:
  - check wether the `current_score` contains Inf.
  - check the feasibility of an optimization problem in the spirit of a convex
  indicator function.
"""
function _infeasibility_check!(stp::NLPStopping, x::T) where {T}
  #=
  #- evaluate the constraint function if `state.cx` is `nothing` and store in `state`.
  #- check the Inf-norm of the violation ≤ stp.meta.atol
   if stp.pb.meta.ncon != 0 #if the problems has constraints, check |c(x)|
    cx = stp.current_state.cx
    if cx == _init_field(typeof(stp.current_state.cx))
     cx = cons(stp.pb, x)
    end
    vio = max.(max.(cx - stp.pb.meta.ucon, 0.), max.(stp.pb.meta.lcon - cx, 0.))
    tol = Inf #stp.meta.atol
    stp.meta.infeasible = _inequality_check(vio, stp.meta.atol, 0.) ? true : stp.meta.infeasible
   end
   =#

  if stp.pb.meta.minimize
    vio = any(z -> z == Inf, stp.current_state.current_score)
    if vio
      stp.meta.infeasible = true
    end
  else
    vio = any(z -> z == -Inf, stp.current_state.current_score)
    if vio
      stp.meta.infeasible = true
    end
  end

  return stp.meta.infeasible
end

################################################################################
# Nonlinear problems admissibility functions
# Available: unconstrained_check(...), optim_check_bounded(...), KKT
################################################################################
include("nlp_admissible_functions.jl")

################################################################################
# line search admissibility functions
#
# TODO: change the ls_admissible_functions and use tol_check et tol_check_neg to
# handle the inequality instead of a max.
################################################################################
include("ls_admissible_functions.jl")

#=
"""
"""
function feasibility_optim_check(pb, state; kwargs...)
     vio = _feasibility(pb, state)
     tol = Inf #stp.meta.atol
     return _inequality_check(vio, tol, 0.)
end
=#

################################################################################
# Functions computing Lagrange multipliers of a nonlinear problem
# Available: _compute_mutliplier(...)
################################################################################
include("nlp_compute_multiplier.jl")


text:
using DbgMacro

function foo()
    m = [1 2; 3 4]
    @dumprt asd(1 + x, :y, $z, k=v)
    42
end

expected_foo_out = """
Expr
  head: Symbol call
  args: Array{Any}((5,))
    1: Symbol asd
    2: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol +
        2: Int64 1
        3: Symbol x
    3: QuoteNode
      value: Symbol y
    4: Expr
      head: Symbol \$
      args: Array{Any}((1,))
        1: Symbol z
    5: Expr
      head: Symbol kw
      args: Array{Any}((2,))
        1: Symbol k
        2: Symbol v
"""

using Suppressor: @capture_out

out = @capture_out global result = foo()
@test result == 42
@test out == expected_foo_out


text:
module EMLreboiler_cost
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...stage_separators.EMLreboiler
	import EMLtypes.length
	include("reboiler_cost/reboiler_cost.jl")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Kissat"
version = v"2.0.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/arminbiere/kissat.git", "00a3a338e3433b54478efb0f7be0a694b01f0eb9"),
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/kissat
./configure -shared
make
mkdir -p ${libdir} ${bindir}
cp build/kissat${exeext} ${bindir}/.
cp build/libkissat.so "$libdir/libkissat.${dlext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=Sys.iswindows, experimental=true)
# The products that we will ensure are always built
products = [
    ExecutableProduct("kissat", :kissat),
    LibraryProduct("libkissat", :libkissat)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is automatically created by the folowing script. Only re-run it if
# AMPL adds new opcodes (which is unlikely).
# ```julia
# using ASL_jll
# open("opcode.jl", "w") do io
#     filename = joinpath(ASL_jll.artifact_dir, "include", "opcode.hd")
#     for line in readlines(filename)
#         items = split(line, c -> c == '\t' || c == ' '; keepempty = false)
#         println(io, "const ", items[2], " = ", items[3])
#     end
# end
# ```

const OPPLUS = 0
const OPMINUS = 1
const OPMULT = 2
const OPDIV = 3
const OPREM = 4
const OPPOW = 5
const OPLESS = 6
const MINLIST = 11
const MAXLIST = 12
const FLOOR = 13
const CEIL = 14
const ABS = 15
const OPUMINUS = 16
const OPOR = 20
const OPAND = 21
const LT = 22
const LE = 23
const EQ = 24
const GE = 28
const GT = 29
const NE = 30
const OPNOT = 34
const OPIFnl = 35
const OP_tanh = 37
const OP_tan = 38
const OP_sqrt = 39
const OP_sinh = 40
const OP_sin = 41
const OP_log10 = 42
const OP_log = 43
const OP_exp = 44
const OP_cosh = 45
const OP_cos = 46
const OP_atanh = 47
const OP_atan2 = 48
const OP_atan = 49
const OP_asinh = 50
const OP_asin = 51
const OP_acosh = 52
const OP_acos = 53
const OPSUMLIST = 54
const OPintDIV = 55
const OPprecision = 56
const OPround = 57
const OPtrunc = 58
const OPCOUNT = 59
const OPNUMBEROF = 60
const OPNUMBEROFs = 61
const OPATLEAST = 62
const OPATMOST = 63
const OPPLTERM = 64
const OPIFSYM = 65
const OPEXACTLY = 66
const OPNOTATLEAST = 67
const OPNOTATMOST = 68
const OPNOTEXACTLY = 69
const ANDLIST = 70
const ORLIST = 71
const OPIMPELSE = 72
const OP_IFF = 73
const OPALLDIFF = 74
const OPSOMESAME = 75
const OP1POW = 76
const OP2POW = 77
const OPCPOW = 78
const OPFUNCALL = 79
const OPNUM = 80
const OPHOL = 81
const OPVARVAL = 82
const N_OPS = 83