text:
using Revise
using Plots
using Pkg
# Pkg.activate(".")
using RadiativeTransfer
using RadiativeTransfer.Architectures
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using RadiativeTransfer.SolarModel
using InstrumentOperator
using Interpolations
using Polynomials
using ForwardDiff
using Distributions
using NCDatasets

## Atmospheric Radiative Transfer

# Load parameters from file
parameters = vSmartMOM.parameters_from_yaml("RadiativeTransfer/test/helper/O2Parameters.yaml")
FT = Float64

# oco2_L1bScND_18688a_180105_B8100r_180206190633

oco_file = "/net/fluo/data1/group/oco2/L1bSc/oco2_L1bScGL_15258a_170515_B10003r_200214061601.h5"
oco_met_file = "/net/fluo/data1/group/oco2/L2Met/oco2_L2MetGL_26777a_190715_B10003r_200429213029.h5"
ils_file = "/home/rjeyaram/RadiativeTransfer/src/vSmartMOM/ils_oco2.json"

fp = 5
iOrbit = 4400
band = 1

sza_ = 32.4436
vza_ = [0.072]

function conv_spectra_local(m::VariableKernelInstrument, ν, spectrum; stride=1)
    # FT = eltype(m.ν_out)
    # Define grid where to perform convolution:

    # Padding at both sides required:
    off = ceil(Int, size(m.kernel, 1) / 2)
    ind = off:stride:(length(ν) - off)

    # knots where convolution will be applied to
    knots = view(ν, ind)
    te = LinearInterpolation(m.ν_out, Float32.(m.ind_out))
    spec_out = zeros(Real, length(knots));
    for i in eachindex(knots)
        # Simple first, nearest neighbor ILS
        ind_fraction = round(Int, te(knots[i]));
        kernel = view(m.kernel, :, ind_fraction)
        for j in eachindex(kernel)
            spec_out[i] += kernel[j] * spectrum[ind[i] + j]
        end
    end
    # Change this later to only perform conv around output grid!
    fin = LinearInterpolation(ν[ind], spec_out; extrapolation_bc=Interpolations.Flat())
    return fin(m.ν_out)
end;

# Runner is used to set AD fields as duals
function runner!(y, x, parameters=parameters, oco_file=oco_file,
                                              oco_met_file=oco_met_file,
                                              ils_file=ils_file, fp=fp, iOrbit=iOrbit,
                                              sza_=sza_, vza_=vza_)

    # Set parameters fields as the dual numbers
    parameters.brdf = [LambertianSurfaceScalar(x[1])]

    parameters.scattering_params.rt_aerosols[1].τ_ref = x[2];

    # parameters.scattering_params.rt_aerosols[1].aerosol.size_distribution = LogNormal(log(x[3]), log(x[4]), check_args=false)

    # parameters.scattering_params.rt_aerosols[1].aerosol.nᵣ = x[5];
    # parameters.scattering_params.rt_aerosols[1].aerosol.nᵢ = x[6];

    # parameters.scattering_params.rt_aerosols[1].p₀ = x[7];
    # parameters.scattering_params.rt_aerosols[1].σp = x[8];

    # Set profiles properly
    met = Dataset(oco_met_file);
    T_met = met.group["Meteorology"]["temperature_profile_met"][:,fp,iOrbit];
    ak = met.group["MeteorologyDiagnostics"]["ak"][:,fp,iOrbit];
    bk = met.group["MeteorologyDiagnostics"]["bk"][:,fp,iOrbit];
    p_surf = met.group["Meteorology"]["surface_pressure_met"][fp,iOrbit];
    p_half = (ak + bk * p_surf);
    p_half = vcat(p_half, p_surf);
    q = met.group["Meteorology"]["specific_humidity_profile_met"][:,fp,iOrbit];
    parameters.p = p_half / 100
    parameters.q = q # zeros(size(q))
    parameters.T = T_met
    parameters.sza = sza_
    parameters.vza = vza_

    @show parameters.sza, parameters.vza

    # @show parameters.sza
    # plot!(parameters.p)
    # return

    model = model_from_parameters(parameters);

    # Run the model to obtain reflectance matrix
    R = vSmartMOM.rt_run(model, i_band=1);

    # Produce black-body in wavenumber range
    T = 5777
    λ_grid = reverse(1e4 ./ parameters.spec_bands[1]) #collect(757:0.01:777)
    black_body = planck_spectrum_wl(T, λ_grid) * 2.1629e-05 * pi
    black_body = SolarModel.watts_to_photons(λ_grid, black_body)

    # Get solar transmittance spectrum
    solar_transmission = solar_transmission_from_file("/home/rjeyaram/RadiativeTransfer/src/solar_merged_20160127_600_26316_100.out", parameters.spec_bands[1])

    # Get outgoing solar radiation
    sun_out = reverse(solar_transmission) .* black_body

    # Apply Earth reflectance matrix
    earth_out = sun_out .* reverse(R[:])

    # Set up InstrumentOperator
    ils_Δ, ils_in, dispersion = InstrumentOperator.read_ils_table(oco_file, ils_file);

    # Define model grid:
    res = 0.001

    # Just consider the ILS within ± 0.35nm
    grid_x = -0.35e-3:res*1e-3:0.35e-3

    extended_dims = [fp,band] # Footprint, band

    # Re-interpolate I from ν_grid to new grid/resolution
    interp_I = LinearInterpolation(λ_grid, earth_out);
    wl = 757.5:res:771.0
    I_wl = interp_I(wl/1000)

    # Pixels to be used
    ind_out = collect(0:1015);

    # Eventual grid of OCO-2 for Band 1, FP 5:
    dispPoly = Polynomial(view(dispersion, :, extended_dims...))
    ν = Float32.(dispPoly.(1:1016))

    # Prepare ILS table:
    ils_pixel   = InstrumentOperator.prepare_ils_table(grid_x, ils_in, ils_Δ,extended_dims)
    oco2_kernel = VariableKernelInstrument(ils_pixel, ν, ind_out)

    # Convolve input spectrum with variable kernel
    I_conv = conv_spectra_local(oco2_kernel, wl./1e3, I_wl)

    y[:] = I_conv[:]

end

directory = "/net/fluo/data1/group/oco2/L1bSc/"
files_list = filter(x->endswith(x, ".h5") && startswith(x, "oco2_L1bScGL_26"), readdir(directory))

met_files_list = filter(x->endswith(x, ".h5"), readdir("/net/fluo/data1/group/oco2/L2Met/"))

I_convs_all = zeros(1016, length(files_list))
oco2_Abands_all = zeros(1016, length(files_list))

albedos = zeros(length(files_list))
cloud_levels = zeros(length(files_list))

optimal_x = zeros(8, length(files_list))

for i in 1:length(files_list)

    @show i

    oco_file = files_list[i]

    # try
        ocoData = Dataset(directory * oco_file)

        oco_met_file = "/net/fluo/data1/group/oco2/L2Met/" * filter(x->contains(x, oco_file[14:19]), met_files_list)[1]

        oco_file = directory * oco_file
        # @show oco_file
        # @show oco_met_file

        met = Dataset(oco_met_file)
        cloud_levels[i] = met.group["Meteorology"]["cloud_liquid_water_path_met"][fp,iOrbit]

        # println()

        # continue

        ## Getting an OCO spectrum for fun to fit:
        o2AbandSpectra = ocoData.group["SoundingMeasurements"]["radiance_o2"]
        SoundingGeometry = ocoData.group["SoundingGeometry"]
        vza = SoundingGeometry["sounding_zenith"]
        sza = SoundingGeometry["sounding_solar_zenith"]
        lat = SoundingGeometry["sounding_latitude"]
        lon = SoundingGeometry["sounding_longitude"]
        o2AbandSpectra = ocoData.group["SoundingMeasurements"]["radiance_o2"]
        oco2_Aband = o2AbandSpectra[:,fp,iOrbit]

        println("$(round(Float64(lat[fp,iOrbit]), digits=2)) $(round(Float64(lon[fp,iOrbit]), digits=4))")
        # continue

        sza_ = sza[fp,iOrbit]
        vza_ = [vza[fp,iOrbit]]

        x = FT[0.05,
        0.05]#,
        # 1.3,
        # 2.0,
        # 1.3,
        # 0.00000001,
        # 90000,
        # 5000.0]

        # Run FW model:
        # @time runner(x);
        I_conv = zeros(1016)
        dfdx = ForwardDiff.jacobian(runner!, I_conv, x);

        dfdx[:,2] .= 0;

        @show size(dfdx)

        y = oco2_Aband;
        Fx = I_conv;
        K = dfdx;

        # Regular least squares:
        dx = K \ (y-Fx)

        I_convs_all[:,i] = I_conv
        oco2_Abands_all[:,i] = oco2_Aband

        # optimal_x[:,i] = x+dx
        albedos[i] = x[1] + dx[1]

        # @show dx, x+dx

    # catch
    #     println("Error on " * oco_file)
    #     continue
    # end

end


text:
export Table

mutable struct Table{T<:AbstractFloat} <: Environment
  worldAABBmin::Vector{T}
  worldAABBmax::Vector{T}
  keepin_zones::Vector
  keepout_zones::Vector
  obstacle_set::Vector
end

function Table{T}(room::Symbol=:ames) where T
  # 2D granite table

  if room == :ames
    # NASA Ames table (from Andrew Symington)
    # CoM coordinates + robot radius
    radius = 0.15*sqrt(2)
    worldAABBmin, worldAABBmax = [-0.5,-0.75,0.] .- radius, [0.75,0.75,0.001] .+ radius
  else room == :stanford
    # Durand 010 table
    ft2m = 0.3048
    worldAABBmin, worldAABBmax = [0.,0.,0.], [12.,9.,0.001]*ft2m
  end

  return Table(worldAABBmin, worldAABBmax)
end
Table(room::Symbol=:ames,::Type{T} = Float64; kwargs...) where {T} = Table{T}(room; kwargs...)

function Table(worldAABBmin::AbstractArray{T}, worldAABBmax::AbstractArray{T}) where T
  keepin_zones = Vector{HyperRectangle}(undef, 0)
  (length(worldAABBmin) == 2) ? push!(worldAABBmin, T(-1)) : nothing
  (length(worldAABBmax) == 2) ? push!(worldAABBmax, T(1)) : nothing
  push!(keepin_zones, HyperRectangle(worldAABBmin..., (worldAABBmax-worldAABBmin)...))

  # Add obstacles surrounding table
  keepout_zones = Vector{GeometryTypes.GeometryPrimitive}(undef, 0)
  koz_min = []
  koz_max = []
  a = 10.  # keep-out zone size (should be larger than 1.05)
  push!(koz_min, [worldAABBmax[1],   -a, -a])  # x
  push!(koz_max, [worldAABBmax[1]+a,  a,  a])
  push!(koz_min, [worldAABBmin[1]-a, -a, -a])  # -x
  push!(koz_max, [worldAABBmin[1],    a,  a])
  push!(koz_min, [-a, worldAABBmax[2],   -a])  # y
  push!(koz_max, [ a, worldAABBmax[2]+a,  a])
  push!(koz_min, [-a, worldAABBmin[2]-a, -a])  # -y
  push!(koz_max, [ a, worldAABBmin[2],    a])

  for i in 1:length(koz_min)
    push!(keepout_zones, HyperRectangle(koz_min[i]..., (koz_max[i]-koz_min[i])...))
  end
  obstacle_set = Vector{GeometryTypes.GeometryPrimitive}(undef, 0)

  return Table{T}(worldAABBmin, worldAABBmax, keepin_zones, keepout_zones, obstacle_set)
end

ext: jl


text:
module Test_covers

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

println("Test_covers")
@testset "covers " begin
Random.seed!(2123)

@testset "grmcover1" begin

m1 = 0
sys = rdss(0,0,0);
@time sysx, sysy, info = grmcover1(sys, m1)
isys1 = 1:m1; isys2 = m1+1:size(sys,2);

@test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy) == 0   &&
      info.stdim == []

sys1 = rdss(0,0,0); sys2 = rdss(0,0,0);
@time sysx, sysy, info = grmcover1(sys1, sys2)

@test gnrank(sysx-sys1-sys2*sysy) == 0   &&
      info.stdim == []

m1 = 0
sys = rdss(1,1,1);
@time sysx, sysy, info = grmcover1(sys, m1)
isys1 = 1:m1; isys2 = m1+1:size(sys,2);

@test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy) == 0   &&
      info.stdim == []

sys1 = rdss(1,1,0); sys2 = rdss(1,1,1);
@time sysx, sysy, info = grmcover1(sys1, sys2)

@test gnrank(sysx-sys1-sys2*sysy) == 0   &&
      info.stdim == []

A = [0 0 0 0 -1;
1 0 0 0 -1;
0 0 -1 -1 0;
0 0 -1 1 0;
1 1 1 0 0.];
B = [1 0 0 0;
0 0 1 0;
0 1 0 0;
0 1 0 1;
0 -1 0 0.];
C = Matrix{Float64}(I,5,5);
E = Matrix{Float64}(I,5,5);
D = zeros(5,4);
sys = dss(A,E,B,C,D);

for m1 = 1:4
    @time sysx, sysy, info = grmcover1(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0   &&
          info.stdim == [m1, 1]
end

for m1 = 1:4
    sys1 = dss(A,E,B[:,1:m1],C,D[:,1:m1]);
    sys2 = dss(A,E,B[:,m1+1:4],C,D[:,m1+1:4]);
    @time sysx, sysy, info = grmcover1(sys1, sys2)
    @test gnrank(sysx-sys1-sys2*sysy, atol1 = 1.e-7) == 0    &&
          info.stdim == [m1, 1]
end

sys = dss(rand(3,4))
m1 = 2
@time sysx, sysy, info = grmcover1(sys, m1)
isys1 = 1:m1; isys2 = m1+1:size(sys,2);

@test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0 &&
info.stdim == []

fast = true; Ty = Complex{Float64}; Ty = Float64
m = 6; n = 10; p = 3;
m = 6; n = 5; p = 3;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol = 1.e-7,atol2=1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, infinite poles
for m1 = 0:m
    sys = [ rdss(n,p,m1,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)])  rdss(n,p,m-m1,T = Ty,disc=false,id = ones(Int,3)) ] ;
    @time sysx, sysy, info = grmcover1(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2 = 1.e-7) == 0
end

# discrete, descriptor, uncontrollable poles
sys = rdss(n,p,m,T = Ty, disc=true,nfuc = 4);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduc=[3*ones(Int,1);2*ones(Int,1)]);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# discrete, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduc=[3*ones(Int,1);2*ones(Int,1)]);
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

end
end

end # grmcover1

@testset "glmcover1" begin

p1 = 0
sys = rdss(0,0,0);
@time sysx, sysy, info = glmcover1(sys, p1)
osys1 = 1:p1; osys2 = p1+1:size(sys,1);

@test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:]) == 0   &&
      info.stdim == []

sys1 = rdss(0,0,0); sys2 = rdss(0,0,0);
@time sysx, sysy, info = glmcover1(sys1, sys2)

@test gnrank(sysx-sys1-sysy*sys2) == 0   &&
      info.stdim == []

p1 = 0
sys = rdss(1,1,1);
@time sysx, sysy, info = glmcover1(sys, p1)
osys1 = 1:p1; osys2 = p1+1:size(sys,1);

@test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:]) == 0    &&
      info.stdim == []

sys1 = rdss(1,0,1); sys2 = rdss(1,1,1);
@time sysx, sysy, info = glmcover1(sys1, sys2)
@test gnrank(sysx-sys1-sysy*sys2) == 0   &&
      info.stdim == []

A = [0 0 0 0 -1;
1 0 0 0 -1;
0 0 -1 -1 0;
0 0 -1 1 0;
1 1 1 0 0.];
B = [1 0 0 0;
0 0 1 0;
0 1 0 0;
0 1 0 1;
0 -1 0 0.];
C = Matrix{Float64}(I,5,5);
E = Matrix{Float64}(I,5,5);
D = zeros(5,4);
sys = gdual(dss(A,E,B,C,D));

for p1 = 1:4
    @time sysx, sysy, info = glmcover1(sys, p1)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0   &&
          info.stdim == [1, p1]
end

for p1 = 1:4
    sys1 = sys[1:p1,:];
    sys2 = sys[p1+1:end,:];
    @time sysx, sysy, info = glmcover1(sys1, sys2)
    @test gnrank(sysx-sys1-sysy*sys2, atol1 = 1.e-7) == 0    &&
          info.stdim == [1, p1]
end

sys = dss(rand(4,3))
p1 = 2
@time sysx, sysy, info = glmcover1(sys, p1)
osys1 = 1:p1; osys2 = p1+1:size(sys,1);
@test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0   &&
     info.stdim == []

fast = true; Ty = Complex{Float64}; Ty = Float64
p = 6; n = 10; m = 3;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
for p1 = 0:p
    @time sysx, sysy, info = glmcover1(sys, p1)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
for p1 = 0:p
    @time sysx, sysy, info = glmcover1(sys, p1, atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol1 = 1.e-7,atol2=1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
for p1 = 0:p
    @time sysx, sysy, info = glmcover1(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# continuous, descriptor, infinite poles
for p1 = 0:p
    sys = [ rdss(n,p1,m,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);  rdss(n,p-p1,m,T = Ty,disc=false,id = ones(Int,3)) ] ;
    @time sysx, sysy, info = glmcover1(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, unobservable poles
sys = rdss(n,p,m,T = Ty, disc=true,nfuo = 4);
for p1 = 0:p
    @time sysx, sysy, info = glmcover1(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# continuous, descriptor, proper, unobservable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduo=[3*ones(Int,1);2*ones(Int,1)]);
for p1 = 0:p
    @time sysx, sysy, info = glmcover1(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

end
end

end # glmcover1

@testset "grmcover2" begin

m1 = 0
sys = rdss(0,0,0);
@time sysx, sysy, info = grmcover2(sys, m1)
isys1 = 1:m1; isys2 = m1+1:size(sys,2);

@test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy) == 0   &&
      info.stdim == []

sys1 = rdss(0,0,0); sys2 = rdss(0,0,0);
@time sysx, sysy, info = grmcover2(sys1, sys2)

@test gnrank(sysx-sys1-sys2*sysy) == 0   &&
      info.stdim == []

m1 = 0
sys = rdss(1,1,1);
@time sysx, sysy, info = grmcover2(sys, m1)
isys1 = 1:m1; isys2 = m1+1:size(sys,2);

@test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy) == 0   &&
      info.stdim == []

sys1 = rdss(1,1,0); sys2 = rdss(1,1,1);
@time sysx, sysy, info = grmcover2(sys1, sys2)

@test gnrank(sysx-sys1-sys2*sysy) == 0   &&
      info.stdim == []

A = [0 0 0 0 -1;
1 0 0 0 -1;
0 0 -1 -1 0;
0 0 -1 1 0;
1 1 1 0 0.];
B = [1 0 0 0;
0 0 1 0;
0 1 0 0;
0 1 0 1;
0 -1 0 0.];
C = Matrix{Float64}(I,5,5);
E = Matrix{Float64}(I,5,5);
D = zeros(5,4);
sys = dss(A,E,B,C,D);

for m1 = 1:4
    @time sysx, sysy, info = grmcover2(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0   &&
          (info.stdim == [m1] || info.stdim == [m1, 1])
end

for m1 = 1:4
    sys1 = dss(A,E,B[:,1:m1],C,D[:,1:m1]);
    sys2 = dss(A,E,B[:,m1+1:4],C,D[:,m1+1:4]);
    @time sysx, sysy, info = grmcover2(sys1, sys2)
    @test gnrank(sysx-sys1-sys2*sysy, atol1 = 1.e-7) == 0    &&
        (info.stdim == [m1] || info.stdim == [m1, 1])
end

sys = dss(rand(3,4))
m1 = 2
@time sysx, sysy, info = grmcover2(sys, m1)
isys1 = 1:m1; isys2 = m1+1:size(sys,2);

@test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0 &&
info.stdim == []

fast = true; Ty = Complex{Float64}; Ty = Float64
m = 6; n = 10; p = 3;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

#error
# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1, atol = 1.e-7,atol2=1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, infinite poles
for m1 = 0:m
    sys = [ rdss(n,p,m1,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)])  rdss(n,p,m-m1,T = Ty,disc=false,id = ones(Int,3)) ] ;
    @time sysx, sysy, info = grmcover2(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# discrete, descriptor, uncontrollable poles
sys = rdss(n,p,m,T = Ty, disc=true,nfuc = 4);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# continuous, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduc=[3*ones(Int,1);2*ones(Int,1)]);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

# discrete, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduc=[3*ones(Int,1);2*ones(Int,1)]);
for m1 = 0:m
    @time sysx, sysy, info = grmcover2(sys, m1, atol = 1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7) == 0
end

end
end

end # grmcover2

@testset "glmcover2" begin

p1 = 0
sys = rdss(0,0,0);
@time sysx, sysy, info = glmcover2(sys, p1)
osys1 = 1:p1; osys2 = p1+1:size(sys,1);

@test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:]) == 0   &&
      info.stdim == []

sys1 = rdss(0,0,0); sys2 = rdss(0,0,0);
@time sysx, sysy, info = glmcover2(sys1, sys2)

@test gnrank(sysx-sys1-sysy*sys2) == 0   &&
      info.stdim == []

p1 = 0
sys = rdss(1,1,1);
@time sysx, sysy, info = glmcover2(sys, p1)
osys1 = 1:p1; osys2 = p1+1:size(sys,1);

@test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:]) == 0    &&
      info.stdim == []

sys1 = rdss(1,0,1); sys2 = rdss(1,1,1);
@time sysx, sysy, info = glmcover2(sys1, sys2)
@test gnrank(sysx-sys1-sysy*sys2) == 0   &&
      info.stdim == []

A = [0 0 0 0 -1;
1 0 0 0 -1;
0 0 -1 -1 0;
0 0 -1 1 0;
1 1 1 0 0.];
B = [1 0 0 0;
0 0 1 0;
0 1 0 0;
0 1 0 1;
0 -1 0 0.];
C = Matrix{Float64}(I,5,5);
E = Matrix{Float64}(I,5,5);
D = zeros(5,4);
sys = gdual(dss(A,E,B,C,D));

for p1 = 1:4
    @time sysx, sysy, info = glmcover2(sys, p1, atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0   &&
    (info.stdim == [p1] || info.stdim == [1, p1])
end

for p1 = 1:4
    sys1 = sys[1:p1,:];
    sys2 = sys[p1+1:end,:];
    @time sysx, sysy, info = glmcover2(sys1, sys2, atol1=1.e-7,atol2=1.e-7)
    @test gnrank(sysx-sys1-sysy*sys2, atol1 = 1.e-7) == 0    &&
    (info.stdim == [p1] || info.stdim == [1, p1])
end

sys = dss(rand(4,3))
p1 = 2
@time sysx, sysy, info = glmcover2(sys, p1, atol1=1.e-7,atol2=1.e-7)
osys1 = 1:p1; osys2 = p1+1:size(sys,1);
@test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0   &&
     info.stdim == []

fast = true; Ty = Complex{Float64}; Ty = Float64
p = 6; n = 10; m = 3;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
for p1 = 0:p
    @time sysx, sysy, info = glmcover2(sys, p1)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
for p1 = 0:p
    @time sysx, sysy, info = glmcover2(sys, p1, atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
for m1 = 0:m
    @time sysx, sysy, info = grmcover1(sys, m1, atol1 = 1.e-7,atol2=1.e-7)
    isys1 = 1:m1; isys2 = m1+1:size(sys,2);
    @test gnrank(sysx-sys[:,isys1]-sys[:,isys2]*sysy,atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
for p1 = 0:p
    @time sysx, sysy, info = glmcover2(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# continuous, descriptor, infinite poles
for p1 = 0:p
    sys = [ rdss(n,p1,m,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);  rdss(n,p-p1,m,T = Ty,disc=false,id = ones(Int,3)) ] ;
    @time sysx, sysy, info = glmcover2(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# discrete, descriptor, unobservable poles
sys = rdss(n,p,m,T = Ty, disc=true,nfuo = 4);
for p1 = 0:p
    @time sysx, sysy, info = glmcover2(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

# continuous, descriptor, proper, unobservable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduo=[3*ones(Int,1);2*ones(Int,1)]);
for p1 = 0:p
    @time sysx, sysy, info = glmcover2(sys, p1,atol1=1.e-7,atol2=1.e-7)
    osys1 = 1:p1; osys2 = p1+1:size(sys,1);
    @test gnrank(sysx-sys[osys1,:]-sysy*sys[osys2,:],atol1=1.e-7,atol2=1.e-7) == 0
end

end
end

end # glmcover2

end #test covers

end #module


text:
module RocketOperatorTest

using Test
using Rocket

@testset "Operator" begin

    struct DummyType <: AbstractOperator end

    struct NotImplementedOperator <: TypedOperator{Int, String} end

    struct ExplicitlyDefinedOperator <: AbstractOperator end
    Rocket.as_operator(::Type{<:ExplicitlyDefinedOperator}) = TypedOperatorTrait{String, Int}()

    struct IdentityIntOperator <: TypedOperator{Int, Int} end
    Rocket.on_call!(::Type{Int}, ::Type{Int}, operator::IdentityIntOperator, source) = source

    struct LeftTypedStringIdentityOperator <: LeftTypedOperator{String} end
    Rocket.on_call!(::Type{String}, ::Type{String}, operator::LeftTypedStringIdentityOperator, source) = source
    Rocket.operator_right(::LeftTypedStringIdentityOperator, ::Type{String}) = String

    struct LeftTypedIntIdentityNotImplementedOperator <: LeftTypedOperator{Int} end
    Rocket.on_call!(::Type{Int}, ::Type{Int}, operator::LeftTypedIntIdentityNotImplementedOperator, source) = source
    # Should be commented as we are testing this case
    # Rocket.operator_right(::LeftTypedIntIdentityNotImplementedOperator, ::Type{Int}) = Int

    struct RightTypedFloatZeroOperator <: RightTypedOperator{Float64} end
    Rocket.on_call!(::Type{L}, ::Type{Float64}, operator::RightTypedFloatZeroOperator, source) where L = Rocket.from([ 0.0 ])

    struct InferableIntoZeroTupleOperator <: InferableOperator end
    Rocket.on_call!(::Type{L}, ::Type{Tuple{L, L}}, operator::InferableIntoZeroTupleOperator, source) where L = Rocket.from([ (zero(L), zero(L)) ])
    Rocket.operator_right(::InferableIntoZeroTupleOperator, ::Type{L}) where L = Tuple{L, L}

    struct InferableNotImplementedOperator <: InferableOperator end
    Rocket.on_call!(::Type{L}, ::Type{L}, operator::InferableIntoZeroTupleOperator, source) where L = source
    # Should be commented as we are testing this case
    # Rocket.operator_right(::InferableNotImplementedOperator, ::Type{L}) where L = L

    struct SomeSubscribable{D} <: Rocket.Subscribable{D} end

    @testset "as_operator" begin
        # Check if arbitrary dummy type has not valid operator type
        @test as_operator(DummyType) === InvalidOperatorTrait()

        # Check if as_operator returns valid operator type for subtypes of Operator abstract type
        @test as_operator(NotImplementedOperator)          === TypedOperatorTrait{Int, String}()
        @test as_operator(LeftTypedStringIdentityOperator) === LeftTypedOperatorTrait{String}()
        @test as_operator(RightTypedFloatZeroOperator)     === RightTypedOperatorTrait{Float64}()
        @test as_operator(InferableIntoZeroTupleOperator)  === InferableOperatorTrait()

        # Check if as_operator returns valid operator type for explicitly defined types
        @test as_operator(ExplicitlyDefinedOperator) === TypedOperatorTrait{String, Int}()
    end

    @testset "|>" begin
        int_source    = SomeSubscribable{Int}()
        string_source = SomeSubscribable{String}()
        float_source  = SomeSubscribable{Float64}()

        # Check if pipe operator throws an error for invalid operator type
        @test_throws InvalidOperatorTraitUsageError int_source |> DummyType()

        # Check if pipe operator throws an error for invalid source and operator data types
        @test_throws InconsistentSourceOperatorDataTypesError string_source |> IdentityIntOperator()
        @test_throws InconsistentSourceOperatorDataTypesError int_source    |> LeftTypedStringIdentityOperator()

        # Check if pipe operator throws an error for not implemented operator
        @test_throws MissingOnCallImplementationError string_source |> ExplicitlyDefinedOperator()
        @test_throws MissingOnCallImplementationError int_source    |> NotImplementedOperator()

        @test_throws MissingOperatorRightImplementationError int_source  |> LeftTypedIntIdentityNotImplementedOperator()
        @test_throws MissingOperatorRightImplementationError int_source  |> InferableNotImplementedOperator()

        # Check if pipe operator calls on_call! for valid operator
        @test int_source |> IdentityIntOperator() === int_source

        # Check if right typed operator may operate on different input source types
        @test int_source    |> RightTypedFloatZeroOperator() == Rocket.from([ 0.0 ])
        @test string_source |> RightTypedFloatZeroOperator() == Rocket.from([ 0.0 ])

        # Check if inferrable operator may operate on different input source types
        @test int_source    |> InferableIntoZeroTupleOperator() == Rocket.from([ (0, 0) ])
        @test float_source  |> InferableIntoZeroTupleOperator() == Rocket.from([ (0.0, 0.0) ])
    end

    @testset "Operators composition with +" begin
        int_source    = SomeSubscribable{Int}()

        operators = IdentityIntOperator() + IdentityIntOperator()

        @test operators isa OperatorsComposition
        @test operators + IdentityIntOperator() isa OperatorsComposition
        @test IdentityIntOperator() + operators isa OperatorsComposition

        @test int_source |> operators === int_source
    end

    @testset "Operators composition with |>" begin
        int_source    = SomeSubscribable{Int}()

        operators = IdentityIntOperator() |> IdentityIntOperator()

        @test operators isa OperatorsComposition
        @test operators |> IdentityIntOperator() isa OperatorsComposition
        @test IdentityIntOperator() |> operators isa OperatorsComposition

        @test int_source |> operators === int_source
    end

end

end


text:
#=- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 Unit Conversion
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -=#

export @unit_convert!

macro unit_convert!(src, from_unit::Symbol, to_unit::Symbol)
    return Expr(:call, Symbol(from_unit, "2", to_unit, "!"), esc(src))
end

#=- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 Converters
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -=#

thz2nj!(f::Real) = 6.626070150e-13 * f
nj2thz!(e::Real) = e / 6.626070150e-13

thz2meV!(f::Real) = 4.135667696 * f
meV2thz!(e::Real) = e / 4.135667696

nm2nj!(λ::Real) = 1.98644586e-7 / λ
nj2nm!(e::Real) = 1.98644586e-7 / e

nm2eV!(λ::Real) = 1239.84198 / λ
eV2nm!(e::Real) = 1239.84198 / e

cm⁻¹2nm!(λ::Real) = 1e7 / λ
nm2cm⁻¹!(λ::Real) = 1e7 / λ

ps2μm!(t::Real) = round(149.896225 * t)
μm2ps!(d::Real) = d / 149.896225

#=- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 Vectorization
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -=#

for fname! in (
        :thz2nj!, :nj2thz!, :thz2meV!, :meV2thz!,
        :nm2nj!, :nj2nm!, :nm2eV!, :eV2nm!,
        :cm⁻¹2nm!, :nm2cm⁻¹!, :ps2μm!, :μm2ps!
    )
    @eval begin
        function $(fname!)(x::AbstractArray{T}) where T<:Real
            @inbounds @simd for i in eachindex(x)
                x[i] = $(fname!)(x[i])
            end
        end
    end
end


text:
function erzeugeFehlerbehafteteDaten(m,prozent)
if prozent > 1 || prozent < 0
    println("Gebt die Prozentzahl als eine Dezimalzahl ein.")
else
    Messwerte = zeros(size(m,1),1);
    B = randn(size(m,1),1);
    max = maximum(B);
    C = abs.(B/max).*prozent;
            #print(C)
    Fehler(Zeile) = m[Zeile]*C[Zeile];
    for i = 1:size(m,1)
        Messwerte[i] = m[i]+Fehler(i);
    end
end
return Messwerte;
end


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

export submit

include("../data.jl")
include("../submit.jl")

include("nnCostFunction.jl")
include("sigmoidGradient.jl")

function submit()
  parts = [
    Part(1, "Feedforward and Cost Function"),
    Part(2, "Regularized Cost Function"),
    Part(3, "Sigmoid Gradient"),
    Part(4, "Neural Network Gradient (Backpropagation)"),
    Part(5, "Regularized Gradient")
  ]
  conf = Conf("neural-network-learning",
              "Neural Networks Learning", parts, solver)

  submitWithConf(conf)
end

function solver(partId)
  # Random Test Cases
  X = reshape(3 * sin(1:1:30), (3, 10))
  Xm = reshape(sin(1:32), (16, 2)) / 5
  ym = (1 + mod(1:16, 4)')'
  t1 = sin(reshape(1:2:24, (4, 3)))
  t2 = cos(reshape(1:2:40, (4, 5)))
  t  = [t1[:] ; t2[:]]
  if partId == 1
    J, _ = nnCostFunction(t, 2, 4, 4, Xm, ym, 0)
    return @sprintf("%0.5f", J)
  elseif partId == 2
    J, _ = nnCostFunction(t, 2, 4, 4, Xm, ym, 1.5)
    return @sprintf("%0.5f", J)
  elseif partId == 3
    return join(map(x -> @sprintf("%0.5f", x), sigmoidGradient(X)), " ")
  elseif partId == 4
    J, grad = nnCostFunction(t, 2, 4, 4, Xm, ym, 0)
    return @sprintf("%0.5f ", J) * join(map(x -> @sprintf("%0.5f", x), grad), " ")
  elseif partId == 5
    J, grad = nnCostFunction(t, 2, 4, 4, Xm, ym, 1.5)
    return @sprintf("%0.5f ", J) * join(map(x -> @sprintf("%0.5f", x), grad), " ")
  end
end


text:
using HDF5
using CSV
using DataFrames

function GetHitInformation(filename::String)::Union{Nothing, GroupedDataFrame}
    extension = filename[findlast(isequal('.'), filename):end]

    if extension == ".csv"
        return GetSRIMcsvHitInformation(filename)
    elseif extension in [".hdf", ".hdf5", ".h5"]
        return GetSRIMhdfHitInformation(filename)
    end
    return nothing
end

function GetSRIMhdfHitInformation(filename::String)::GroupedDataFrame
    @info "Reading SRIM hits info from $filename"
    fid = h5open(filename, "r")
    ions::Array{Int64, 1} = read(fid, "hits/ion")
    # All distances in angstroms
    # energy in keV
    x::Array{Float64, 1} = read(fid, "hits/X")
    y::Array{Float64, 1} = read(fid, "hits/Y")
    z::Array{Float64, 1} = read(fid, "hits/Z")
    eDep::Array{Float64, 1} = read(fid, "hits/E")

    df = DataFrame(ion = ions, X = x, Y = y, Z = z, E = eDep)
    # convert to mm
    df.X *= 1e-7
    df.Y *= 1e-7
    df.Z *= 1e-7

    return groupby(df, :ion)
end

function GetSRIMcsvHitInformation(filename::String)::GroupedDataFrame
    @info "Reading SRIM hits info from $filename"
    df = CSV.read(filename,DataFrame)
    # convert coordinates from angstroms to mm
    df.X *= 1e-7
    df.Y *= 1e-7
    df.Z *= 1e-7
    # energy in keV
    outdf = df[:, ["ion", "X", "Y", "Z", "E"]]
    return groupby(outdf, :ion)
end

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

abstract type AbstractLinearSolver end
abstract type AbstractIterator end

# dummy solver type
mutable struct EmptyLinearSolver <: AbstractLinearSolver end

# dummy functions
introduce(::EmptyLinearSolver) = ""
factorize!(M::EmptyLinearSolver) = M
solve!(::EmptyLinearSolver,x) = x
is_inertia(::EmptyLinearSolver) = false
inertia(::EmptyLinearSolver) = (0,0,0)
improve!(::EmptyLinearSolver) = false
solve_refine!(y,::AbstractIterator,x) = nothing

# LinearSolverExceptions
struct SymbolicException <: Exception end
struct FactorizationException <: Exception end
struct SolveException <: Exception end
struct InertiaException <: Exception end
LinearSolverException=Union{SymbolicException,FactorizationException,SolveException,InertiaException}

# iterative solvers
include("richardson.jl")

# dense solvers
include("lapack.jl")

# direct solvers
include("umfpack.jl")
BLAS.vendor() == :mkl && include("pardisomkl.jl")


text:
using HCubature
using Distributions
using LinearAlgebra

struct Fun
    f::Function
    n::Int
    y::Vector{Float64}
end

Fun( f::Function, n::Int ) = Fun( f, n, zeros(n) )

function (f::Fun)( x::AbstractVector{Float64} )
    factor = 1.0
    for i = 1:f.n
        y = x[i]
        denom = y*(1-y)
        f.y[i] = (2*y-1)/denom
        factor *= (2*y^2 - 2*y + 1)/denom^2
    end
    return f.f( f.y ) * factor
end

function invtransform( x::Float64 )
    if x == 0.0
        return 0.5
    elseif isinf(x)
        return x == -Inf ? 0.0 : 1.0
    else
        return (x-2+sqrt(x^2+4))/(2*x)
    end
end

function nintegrate( f::Fun, a::Vector{Float64}, b::Vector{Float64}; kwargs... )
    for i = 1:f.n
        a[i] = invtransform(a[i])
        b[i] = invtransform(b[i])
    end
    return hcubature( f, a, b; kwargs... )
end

function lenormalcubature( mean, cov, A, b; kwargs... )
    n = length(mean)
    mvn = MvNormal( mean, cov )
    f = Fun( x -> all(A*x .<= b) ? pdf( mvn, x ) : 0.0, n )
    return nintegrate( f, fill(-Inf,n), fill(Inf,n); kwargs... )
end


text:
@testset "ANN" begin
    function ann(w1, w2, w3, x1)
        xx2 = w1 * x1
        x2 = log.(1.0 .+ exp.(xx2))   # soft RELU unit
        xx3 = w2 * x2
        x3 = log.(1.0 .+ exp.(xx3))   # soft RELU unit
        x4 = w3 * x3
        sum(1.0 ./ (1.0 .+ exp.(-x4)))  # sigmoid output
    end

    test_compare(ann; w1=rand(10,10), w2=rand(10,10), w3=rand(1,10), x1=rand(10))

    f = ann
    inputs = [:w1=>rand(10,10), :w2=>rand(10,10), :w3=>rand(1,10), :x1=>rand(10)]
end


text:
import Base.show
export show

# A Constant is simply a wrapper around a native Julia constant
# Hence, we simply display its value
function show(io::IO, x::Constant)
    print(io, "$(x.value)")
end

# A variable, for example, Variable(3, 4), will be displayed as:
# Variable of
# size: (3, 4)
# sign: NoSign()
# vexity: AffineVexity()
function show(io::IO, x::Variable)
    print(io, """Variable of
          size: ($(x.size[1]), $(x.size[2]))
          sign: $(x.sign)
          vexity: $(x.vexity)""")
    if x.value !== nothing
        print(io, "\nvalue: $(x.value)")
    end
end

# A constraint, for example, square(x) <= 4 will be displayed as:
# Constraint:
# <= constraint
# lhs: ...
# rhs: ...
function show(io::IO, c::Constraint)
    print(io, """Constraint:
          $(c.head) constraint
          lhs: $(c.lhs)
          rhs: $(c.rhs)
          vexity: $(vexity(c))""")
end

# SDP constraints are displayed as:
# Constraint:
# sdp constraint
# expression: ...
function show(io::IO, c::SDPConstraint)
    print(io, """Constraint:
          $(c.head) constraint
          expression: $(c.child)
          """)
end

# An expression, for example, 2 * x, will be displayed as:
# AbstractExpr with
# head: *
# size: (1, 1)
# sign: NoSign()
# vexity: AffineVexity()
function show(io::IO, e::AbstractExpr)
    print(io, """AbstractExpr with
          head: $(e.head)
          size: ($(e.size[1]), $(e.size[2]))
          sign: $(sign(e))
          vexity: $(vexity(e))
          """)
end

# A problem, for example, p = minimize(sum(x) + 3, [x >= 0, x >= 1, x <= 10])
# will be displayed as follows:
#
# Problem: minimize `Expression`
#        subject to
#            Constraint: ...
#            Constraint: ...
#            Constraint: ...
#        current status: not yet solved
#
# Once it is solved, the current status would look like:
#        current status: solved with optimal value of 9.0
function show(io::IO, p::Problem)
    print(io, """Problem:
          $(p.head) $(p.objective)
          subject to
          """)
    join(io, p.constraints, "\n\t\t")
    print(io, "\ncurrent status: $(p.status)")
    if p.status == "solved"
        print(io, " with optimal value of $(round(p.optval, digits=4))")
    end
end


text:
import FileIO
using FileIO: File, @format_str

# Save all the key-value pairs in the dict as top-level variables of the JLD
function fileio_save(f::File{format"JLD3"}, dict::AbstractDict; kwargs...)
    jldopen(FileIO.filename(f), "w"; kwargs...) do file
        wsession = JLDWriteSession()
        for (k,v) in dict
            if isa(k, Symbol)
                @warn "you passed a key as a symbol instead of a string, it has been automatically converted"
            elseif !isa(k, AbstractString)
                throw(ArgumentError("keys must be strings (the names of variables), got $k"))
            end
            write(file, String(k), v, wsession)
        end
    end
end

# Or the names and values may be specified as alternating pairs
function fileio_save(f::File{format"JLD3"}, name::AbstractString, value, pairs...; kwargs...)
    if isodd(length(pairs)) || !isa(pairs[1:2:end], Tuple{Vararg{AbstractString}})
        throw(ArgumentError("arguments must be in name-value pairs"))
    end
    jldopen(FileIO.filename(f), "w"; kwargs...) do file
        wsession = JLDWriteSession()
        write(file, String(name), value, wsession)
        for i = 1:2:length(pairs)
            write(file, String(pairs[i]), pairs[i+1], wsession)
        end
    end
end

fileio_save(f::File{format"JLD3"}, value...; kwargs...) = error("must supply a name for each variable")

# load with just a filename returns a dictionary containing all the variables
function fileio_load(f::File{format"JLD3"}; kwargs...)
    jldopen(FileIO.filename(f), "r"; kwargs...) do file
        loadtodict!(Dict{String,Any}(), file)
    end
end

# When called with explicitly requested variable names, return each one
function fileio_load(f::File{format"JLD3"}, varname::AbstractString; kwargs...)
    jldopen(FileIO.filename(f), "r"; kwargs...) do file
        read(file, varname)
    end
end

fileio_load(f::File{format"JLD3"}, varnames::AbstractString...; kwargs...) =
    fileio_load(f, varnames; kwargs...)

function fileio_load(f::File{format"JLD3"}, varnames::Tuple{Vararg{AbstractString}}; kwargs...)
    jldopen(FileIO.filename(f), "r"; kwargs...) do file
        map((var)->jlread(file, var), varnames)
    end
end


text:
using BenchmarkTools
using Measurements
using SpecialFunctions
using QuadGK

const SUITE = BenchmarkGroup()

const a = 12.3 ± 4.5
const b = 67.8 ± 9.0
const vector = [1 ± 0.1 for _ in 1:10000]

int1(x) = exp(x / a)

SUITE["base"] = BenchmarkGroup()
SUITE["uncertain"] = BenchmarkGroup()
SUITE["vector"] = BenchmarkGroup()
SUITE["specfun"] = BenchmarkGroup()
SUITE["quadgk"] = BenchmarkGroup()

SUITE["base"]["create"] = @benchmarkable 4.7 ± 0.3
SUITE["base"]["sum"] = @benchmarkable $(a) + $(b)
SUITE["base"]["multiplication"] = @benchmarkable $(a) * $(b)
SUITE["base"]["sqrt"] = @benchmarkable sqrt($(b))
SUITE["base"]["sin"] = @benchmarkable sqrt($(a))

SUITE["uncertain"]["zeta"] = @benchmarkable @uncertain zeta($(2 ± 0.13))
SUITE["uncertain"]["log"] = @benchmarkable @uncertain log($(9.4 ± 1.3), $(58.8 ± 3.7))
SUITE["uncertain"]["atan"] = @benchmarkable @uncertain atan($(10), $(13.5 ± 0.8))

SUITE["vector"]["sqrt"] = @benchmarkable sqrt.($(vector))
SUITE["vector"]["sin"] = @benchmarkable sin.($(vector))
SUITE["vector"]["gamma"] = @benchmarkable gamma.($(vector))
SUITE["vector"]["sin²+cos²"] = @benchmarkable cos.($(vector)) .^ 2 .+ sin.($(vector)) .^ 2
SUITE["vector"]["sum - all different"] = @benchmarkable sum($(vector))
SUITE["vector"]["mean - all different"] = @benchmarkable mean($(vector))
SUITE["vector"]["prod - all different"] = @benchmarkable prod($(vector))
SUITE["vector"]["sum - all equal"] = @benchmarkable sum(v) setup = (x = 1 ± 0.1; v = [x for _ in eachindex(vector)])
SUITE["vector"]["mean - all equal"] = @benchmarkable mean(v) setup = (x = 1 ± 0.1; v = [x for _ in eachindex(vector)])
SUITE["vector"]["prod - all equal"] = @benchmarkable prod(v) setup = (x = 1 ± 0.1; v = [x for _ in eachindex(vector)])

SUITE["specfun"]["gamma"] = @benchmarkable gamma($(a))

SUITE["quadgk"]["int1"] = @benchmarkable quadgk(int1, 1, 7)
SUITE["quadgk"]["cos"] = @benchmarkable quadgk(cos, $(a), $(b))
SUITE["quadgk"]["sin"] = @benchmarkable quadgk(sin, $(-a), $(a))


text:
using GeometricBase.Config
using GeometricBase.Utils
using GeometricIntegrators.Integrators
using GeometricIntegrators.Solutions
using GeometricIntegrators.Tableaus
using GeometricProblems.LotkaVolterra2d
using SimpleSolvers
using Test

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

const Δt = 0.01
const nt = 10
const q₀ = [1.0, 1.0]
const parameters = (a₁=1.0, a₂=1.0, b₁=-1.0, b₂=-2.0)

ode  = lotka_volterra_2d_ode(q₀; params=parameters)
iode = lotka_volterra_2d_iode(q₀; params=parameters)
lode = lotka_volterra_2d_lode(q₀; params=parameters)

int  = IntegratorFIRK(ode, TableauGauss(8), Δt)
sol  = integrate(ode, int, nt)
refx = sol.q[end]

@testset "$(rpad("1st Order DVIs",80))" begin

    sol = Solution(iode, Δt, nt)
    int = IntegratorDVIA(iode, Δt)
    integrate!(int, sol)
    @test relative_maximum_error(sol.q, refx) < 1E-1

    sol = Solution(iode, Δt, nt)
    int = IntegratorDVIB(iode, Δt)
    integrate!(int, sol)
    @test relative_maximum_error(sol.q, refx) < 1E-1

end

@testset "$(rpad("2nd Order Centred DVIs",80))" begin

    sol = Solution(iode, Δt, nt)
    int = IntegratorCMDVI(iode, Δt)
    integrate!(int, sol)
    @test relative_maximum_error(sol.q, refx) < 4E-3

    sol = Solution(iode, Δt, nt)
    int = IntegratorCTDVI(iode, Δt)
    integrate!(int, sol)
    @test relative_maximum_error(sol.q, refx) < 4E-3

end


text:
export
    FourierNeuralOperator

function FourierNeuralOperator()
    modes = (16, )
    ch = 64 => 64
    σ = gelu

    return Chain(
        Dense(2, 64),
        FourierOperator(ch, modes, σ),
        FourierOperator(ch, modes, σ),
        FourierOperator(ch, modes, σ),
        FourierOperator(ch, modes),
        Dense(64, 128, σ),
        Dense(128, 1),
        flatten
    )
end


text:
module scoring

export score_population

using CUDA
using DataFrames: DataFrame
using Random: rand
using Base.Threads: @threads, nthreads, threadid

include("constants.jl")
include("biclusterseval.jl")
include("evolution.jl")

using .biclusterseval: evaluate_fitness, compress_chromes

function score_population(
    d_input_data::Vector{CuArray{Float32,2}},
    population::Population;
    gpus_num::Int = 1,
    return_score = true,
)
    compressed_chromes, chromes_ids = compress_chromes(population)

    partial_fitnesses = Vector(undef, nthreads())
    @threads for (dev, d_data_subset) in collect(zip(devices(), d_input_data))
        device!(dev)

        d_fitness = CUDA.zeros(Int32, length(population))

        rows_number::UInt32 = size(d_data_subset, 1)

        blocks_per_chromo = ceil(Int, rows_number / BLOCK_SIZE)

        d_compressed_chromes = CuArray{Int32}(undef, length(compressed_chromes))
        copyto!(d_compressed_chromes, compressed_chromes)
        d_chromes_ids = CuArray{Int32}(undef, length(chromes_ids))
        copyto!(d_chromes_ids, chromes_ids)

        @cuda blocks = (length(population), blocks_per_chromo) threads = (1, BLOCK_SIZE) evaluate_fitness(
            d_fitness,
            d_data_subset,
            rows_number,
            d_compressed_chromes,
            d_chromes_ids,
        )

        synchronize()

        partial_fitnesses[threadid()] = Array(d_fitness)

        CUDA.unsafe_free!(d_fitness)
        CUDA.unsafe_free!(d_compressed_chromes)
        CUDA.unsafe_free!(d_chromes_ids)
    end

    fitness = reduce(+, partial_fitnesses)

    return if return_score
        [chromo => score_chromo(chromo, c_fitness) for
        (chromo, c_fitness) in zip(population, fitness)]
    else
        collect(zip(population, fitness))
    end
end

function score_chromo(chromo::Chromo, fitness)::Float64
    rows = fitness
    rows <= 1 && return 0
    cols = length(chromo)
    return 2.0^min(rows - MIN_NO_ROWS, 0) * log2(max(rows - 1, 0)) * cols
end

end


text:
function reldiff(x, y)
    p = (x - y) ./ abs.(x)
    return norm(p, Inf)
end


text:
using Test

@testset "Periodic BC DO" begin
    N = 100
    d_order = 2
    approx_order = 2
    x = collect(1:1.0:N).^2
    A = DerivativeOperator{Float64}(2,2,1.0,N,:periodic,:periodic)
    boundary_points = A.boundary_point_count

    res = A*x
    @test res[boundary_points[1] + 1: N - boundary_points[2]] == 2.0*ones(N - sum(boundary_points));

    N = 1000
    d_order = 4
    approx_order = 10
    y = collect(1:1.0:N).^4 - 2*collect(1:1.0:N).^3 + collect(1:1.0:N).^2;
    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:periodic,:periodic)
    boundary_points = A.boundary_point_count

    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^-1; # Float64 is less stable

    A = DerivativeOperator{BigFloat}(d_order,approx_order,one(BigFloat),N,:periodic,:periodic)
    y = convert(Array{BigFloat, 1}, y)
    res = A*y;
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^-approx_order;
end

@testset "Periodic BC FD" begin
    N = 100
    d_order = 2
    approx_order = 2
    x = collect(1:1.0:N).^2
    A = FiniteDifference{Float64}(2,2,ones(N-1),N,:periodic,:periodic)
    boundary_points = A.boundary_point_count

    res = A*x
    @test res[boundary_points[1] + 1: N - boundary_points[2]] == 2.0*ones(N - sum(boundary_points));

    N = 1000
    d_order = 4
    approx_order = 10
    y = collect(1:1.0:N).^4 - 2*collect(1:1.0:N).^3 + collect(1:1.0:N).^2;
    A = FiniteDifference{Float64}(d_order,approx_order,ones(N-1),N,:periodic,:periodic)
    boundary_points = A.boundary_point_count

    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^-1; # Float64 is less stable

    A = FiniteDifference{BigFloat}(d_order,approx_order,ones(BigFloat,N-1),N,:periodic,:periodic)
    y = convert(Array{BigFloat, 1}, y)
    res = A*y;
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^-approx_order;
end


text:
module TestNullModels
  using Test
  using EcologicalNetworks

  # Generate some data

  A = [true true false; true true true; false true false]
  N = UnipartiteNetwork(A)

  n1 = null1(N)
  n2 = null2(N)
  n4 = null4(N)
  n3o = null3(N; dims=1)
  n3i = null3(N; dims=2)

  @test typeof(n1) <: UnipartiteProbabilisticNetwork
  @test typeof(n2) <: UnipartiteProbabilisticNetwork
  @test typeof(n4) <: UnipartiteProbabilisticNetwork
  @test typeof(n3o) <: UnipartiteProbabilisticNetwork
  @test typeof(n3i) <: UnipartiteProbabilisticNetwork

  B = BipartiteNetwork(A)
  @test typeof(null1(B)) <: BipartiteProbabilisticNetwork
  @test typeof(null2(B)) <: BipartiteProbabilisticNetwork
  @test typeof(null3(B; dims=2)) <: BipartiteProbabilisticNetwork
  @test typeof(null3(B; dims=1)) <: BipartiteProbabilisticNetwork
  @test typeof(null4(B)) <: BipartiteProbabilisticNetwork

  # Null model 1
  @test n1[1,1] ≈ 6.0 / 9.0
  @test n1[2,2] ≈ 6.0 / 9.0
  @test n1[3,3] ≈ 6.0 / 9.0

  # Null model 2
  @test n2[1,3] ≈ 0.5
  @test n2[2,2] ≈ 1.0
  @test n2[1,1] ≈ 2.0 / 3.0

  # Null model 2 (mult)
  @test n4[1,3] ≈ 2.0 * 1.0 / 36.0
  @test n4[2,2] ≈ 3.0 * 3.0 / 36.0
  @test n4[1,1] ≈ 2.0 * 2.0 / 36.0

  # Null model 3
  @test n3o[1,1] ≈ 2.0 / 3.0
  @test n3o[2,1] ≈ 1.0
  @test n3o[3,3] ≈ 1.0 / 3.0

  @test n3i[1,1] ≈ 2.0 / 3.0
  @test n3i[2,2] ≈ 1.0
  @test n3i[3,3] ≈ 1.0 / 3.0

  # filter

  C = copy(B)
  C[1,1] = false  # remove one interaction
  # make version of the network with the first interaction removed

  # test if the zoo of B is the same as the result of the filter of C at (1,1)
  @test linearfilter(C)[1,1] ≈ linearfilterzoo(B)[1,1]
  @test linearfilter(C)[1,1] != linearfilter(B)[1,1]

end


text:
using Mimi
using Distributions
using Test

@defcomp test1 begin
    p = Parameter(default = 5)
    function run_timestep(p, v, d, t) end
end

@defcomp test2 begin
    p = Parameter(default = 5)
    function run_timestep(p, v, d, t) end
end

@defcomp test3 begin
    a = Parameter(default = 5)
    function run_timestep(p, v, d, t) end
end

##
## Tests for set_translist_modelparams
##

sd = @defsim begin
    sampling(LHSData)
    p = Normal(0, 1)
end

#------------------------------------------------------------------------------
# Test a failure to find the unshared parameter in any components
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test3)

fail_expr1 = :(
    run(sd, m, 100)
)

err1 = try eval(fail_expr1) catch err err end
@test occursin("Cannot resolve because p not found in any of the components of Model1.", sprint(showerror, err1))

#------------------------------------------------------------------------------
# Test a failure due to finding the unshared parameter in more than one component
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test1)
add_comp!(m, test2)

fail_expr2 = :(
    run(sd, m, 100)
)

err2 = try eval(fail_expr2) catch err err end
@test occursin("Cannot resolve because parameter name p found in more than one component of Model1", sprint(showerror, err2))

#------------------------------------------------------------------
# Test a failure due to finding an unshared parameter in one model, but not
# the other

m1 = Model()
set_dimension!(m1, :time, 2000:10:2050)
add_comp!(m1, test1)

m2 = Model()
set_dimension!(m2, :time, 2000:10:2050)
add_comp!(m2, test3)

fail_expr3 = :(
    run(sd, [m1, m2], 100)
)

err3 = try eval(fail_expr3) catch err err end
@test occursin("Cannot resolve because p not found in any of the components of Model2", sprint(showerror, err3))

#------------------------------------------------------------------------------
# Test success cases

# unshared parameter set by default
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test1)
run(sd, m, 100)

# shared parameter
m1 = Model()
set_dimension!(m1, :time, 2000:10:2050)
add_comp!(m1, test1)
add_comp!(m1, test2)
add_shared_param!(m1, :model_p, 5)
connect_param!(m1, :test1, :p, :model_p)
connect_param!(m1, :test2, :p, :model_p)
run(sd, m, 100)

# unshared parameter in both models with different names
m1 = Model()
set_dimension!(m1, :time, 2000:10:2050)
add_comp!(m1, test1)

m2 = Model()
set_dimension!(m2, :time, 2000:10:2050)
add_comp!(m2, test2)

run(sd, [m1, m2], 100)

#------------------------------------------------------------------------------
# Test warning for compname.paramname syntax in a shared parameter case

# shared parameter
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test1) # test1.p is unshared with default of 5
add_comp!(m, test2) # test2.p is unshared with default of 5
add_shared_param!(m, :model_p, 6)
connect_param!(m, :test1, :p, :model_p)
# now test1.p is shared, connected to model_p with a value of 6, test2.p is still unshared

sd = @defsim begin
    sampling(LHSData)
    test2.p = Normal(0, 1)
end
run(sd, m, 2) # no warnings

sd = @defsim begin
    sampling(LHSData)
    test1.p = Normal(0, 1)
end
msg = "Parameter indicated in `defsim` as test1.p is connected to a SHARED parameter model_p. Thus the value will be varied in all component parameters connected to that shared model parameter.  We suggest using model_p = distribution syntax to be transparent about this."
@test_logs (:warn, msg) run(sd, m, 2) #should work but warn

sd = @defsim begin
    sampling(LHSData)
    p = Normal(0, 1)
end
@test_throws ErrorException run(sd, m, 2) #should error, can't resolve

##
## Tests for set_translist_modelparams with a default (not shared)
##

sd = @defsim begin
    sampling(LHSData)
    test1.p = Normal(0, 1)
end

# simple case
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test1)
run(sd, m, 100)

# component not found
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test2)
fail_expr = :(run(sd, m, 100))
err4 = try eval(fail_expr) catch err err end
@test occursin("Component test1 does not exist in Model1.", sprint(showerror, err4))

m1 = Model()
set_dimension!(m1, :time, 2000:10:2050)
add_comp!(m1, test2)
m2 = Model()
set_dimension!(m2, :time, 2000:10:2050)
add_comp!(m2, test1)
fail_expr = :(run(sd, [m1, m2], 100))
err5 = try eval(fail_expr) catch err err end
@test occursin("Component test1 does not exist in Model1.", sprint(showerror, err5))

# transform used only for one of the component's parameters p
m = Model()
set_dimension!(m, :time, 2000:10:2050)
add_comp!(m, test1)
add_comp!(m, test2) # no transform used
run(sd, m, 100)

# two models, both with component parameter pair
m1 = Model()
set_dimension!(m1, :time, 2000:10:2050)
add_comp!(m1, test1)

m2 = Model()
set_dimension!(m2, :time, 2000:10:2050)
add_comp!(m2, test1)

run(sd, [m1, m2], 100)


text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "make cue matrix for latin" begin
    latin = DataFrame(CSV.File(joinpath("data", "latin_mini.csv")))
    latin_cue_obj_train = JudiLing.make_cue_matrix(
        latin,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    latin_val = latin[101:150, :]
    latin_cue_obj_val = JudiLing.make_cue_matrix(
        latin_val,
        latin_cue_obj_train,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )
end

@testset "make cue matrix for french" begin
    french = DataFrame(CSV.File(joinpath("data", "french_mini.csv")))
    french_cue_obj_train = JudiLing.make_cue_matrix(
        french,
        grams = 2,
        target_col = :Syllables,
        tokenized = true,
        sep_token = "-",
        keep_sep = true,
    )

    french_val = french[101:150, :]
    french_cue_obj_val = JudiLing.make_cue_matrix(
        french_val,
        french_cue_obj_train,
        grams = 2,
        target_col = :Syllables,
        tokenized = true,
        sep_token = "-",
        keep_sep = true,
    )
end

@testset "make cue matrix for utterance" begin
    utterance =
        DataFrame(CSV.File(joinpath("data", "utterance_mini.csv")))
    utterance_cue_obj_train = JudiLing.make_cue_matrix(
        utterance,
        grams = 3,
        target_col = :Words,
        tokenized = true,
        sep_token = ".",
        keep_sep = true,
    )

    utterance_val = utterance[101:150, :]
    utterance_cue_obj_val = JudiLing.make_cue_matrix(
        utterance_val,
        utterance_cue_obj_train,
        grams = 3,
        target_col = :Words,
        tokenized = true,
        sep_token = ".",
        keep_sep = true,
    )
end

@testset "make combined cue matrix" begin
    latin_full = DataFrame(CSV.File(joinpath(
        @__DIR__,
        "data",
        "latin_mini.csv",
    )))

    latin_train = latin_full[1:3, :]
    latin_val = latin_full[10:15, :]

    cue_obj_train, cue_obj_val = JudiLing.make_combined_cue_matrix(
        latin_train,
        latin_val,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )

    @test cue_obj_train.C[1, 3] == 1
    @test cue_obj_val.C[1, 3] == 0
    @test cue_obj_train.i2f[3] == "oco"
    @test cue_obj_val.i2f[3] == "oco"
    @test cue_obj_train.f2i["oco"] == 3
    @test cue_obj_val.f2i["oco"] == 3
    @test cue_obj_train.A[1, 3] == 0
    @test cue_obj_val.A[1, 3] == 0
end


text:
"""
Abstract definition of a lattice.
"""
abstract type AbstractLattice end

# AbstractLattice interface: mandatory
# TODO: This needs to be updated. There must be a general way to access sites and bonds.
"""
    length(l::AbstractLattice)

Number of lattice sites.
"""
@inline Base.length(l::AbstractLattice) = MethodError(length, (l))

"""
    size(l::AbstractLattice)

Returns the linear system size of a lattice.
"""
@inline Base.size(l::AbstractLattice) = error("Lattice $(typeof(l)) doesn't implement `length`.")

"""
    neighbors(l::AbstractLattice[, directed=Val(false)])

Returns an iterator over bonds, given as tuples (source index, target index). If
`directed = Val(true)` bonds are assumed to be directed, i.e. both
`(1, 2)` and `(2, 1)` are included. If `directed = Val(false)` bonds are
assumed to be undirected, i.e. `(1, 2)` and `(2, 1)` are assumed to be
equivalent and only one of them will be included.
"""
@inline neighbors(l::AbstractLattice) = neighbors(l::AbstractLattice, Val(false))
@inline function neighbors(l::AbstractLattice, directed::Val{true})
    _neighbors(has_neighbors_table(l), l, directed)
end
@inline function neighbors(l::AbstractLattice, directed::Val{false})
    _neighbors(has_bonds_table(l), l, directed)
end

# AbstractLattice interface: optional
"""
    neighbors(l::AbstractLattice, site_index::Integer)

Returns a list of site indices neighboring the specified `site_index`.
"""
@inline function neighbors(l::AbstractLattice, site_index::Integer)
    _neighbors(has_neighbors_table(l), l, site_index)
end

# Abstract trait
abstract type LatticeProperties end

"""
    HasNeighborsTable trait

Indicates that a lattice `l` has a neighbors lookup table as a field `l.neighs`.

A neighbors lookup table is a `Matrix{Int}` of size `(nneighbors, nsites)` such that
`l.neighs[i,j]` gives the site index of the `i`th neighbor of site `j`.

A lattice should implement `has_neighbors_table(::MyLattice) = HasNeighborsTable()`
to indicate the trait.
"""
struct HasNeighborsTable <: LatticeProperties end

"""
    HasBondsTable trait

Indicates that a lattice `l` has a bonds table as a field `l.bonds`.

A bonds table is a `Matrix{Int}` of size `(nbonds, 3)` where the columns mean
`src`, `trg`, `type`, respectively.

A lattice should implement `has_bonds_table(::MyLattice) = HasBondsTable()` to
indicate the trait.
"""
struct HasBondsTable <: LatticeProperties end

has_neighbors_table(::AbstractLattice) = nothing
has_bonds_table(::AbstractLattice) = nothing

# neighbors(lattice, directed::Val)
@inline function _neighbors(::Nothing, l::LT, directed::V) where {LT<:AbstractLattice, V<:Val}
    error("Lattice $(LT) doesn't implemen `neighbors(::$(LT), ::$V`")
end
@inline function _neighbors(::HasNeighborsTable, l::AbstractLattice, directed::Val{true})
    ([src, trg] for src in 1:length(l) for trg in l.neighs[:, src])
end
@inline function _neighbors(::HasBondsTable, l::AbstractLattice, directed::Val{false})
    (l.bonds[i, 1:2] for i in 1:size(l.bonds, 1))
end

# neighbors(lattice, site_index)
@inline function _neighbors(::Nothing, l::LT, site_index::Integer) where {LT<:AbstractLattice}
    error("Lattice $(LT) doesn't implement `neighbors(::$(LT), site_index`.")
end
@inline function _neighbors(::HasNeighborsTable, l::AbstractLattice, site_index::Integer)
    l.neighs[:, site_index]
end

################################################################################
### Generic stuff
################################################################################

function reciprocal_vectors(lattice::AbstractLattice, L = lattice.L)
    v1, v2 = latticeVectors(getReciprocalUnitcell(unitcell(lattice))) ./ L
    r1 = cross([v1..., 0.0], [0,0,1])
    r2 = cross([v2..., 0.0], [0,0,1])
    return r1, r2
end


text:
# Julia compiler integration

## cache

using Core.Compiler: CodeInstance, MethodInstance

struct CodeCache
    dict::Dict{MethodInstance,Vector{CodeInstance}}
    CodeCache() = new(Dict{MethodInstance,Vector{CodeInstance}}())
end

function Base.show(io::IO, ::MIME"text/plain", cc::CodeCache)
    print(io, "CodeCache with $(mapreduce(length, +, values(cc.dict); init=0)) entries: ")
    for (mi, cis) in cc.dict
        println(io)
        print(io, "  ")
        show(io, mi)

        function worldstr(min_world, max_world)
            if min_world == typemax(UInt)
                "empty world range"
            elseif max_world == typemax(UInt)
                "worlds $(Int(min_world))+"
            else
                "worlds $(Int(min_world)) to $(Int(max_world))"
            end
        end

        for (i,ci) in enumerate(cis)
            println(io)
            print(io, "    CodeInstance for ", worldstr(ci.min_world, ci.max_world))
        end
    end
end

function Core.Compiler.setindex!(cache::CodeCache, ci::CodeInstance, mi::MethodInstance)
    # make sure the invalidation callback is attached to the method instance
    callback(mi, max_world) = invalidate(cache, mi, max_world)
    if !isdefined(mi, :callbacks)
        mi.callbacks = Any[callback]
    else
        if all(cb -> cb !== callback, mi.callbacks)
            push!(mi.callbacks, callback)
        end
    end

    cis = get!(cache.dict, mi, CodeInstance[])
    push!(cis, ci)
end

# invalidation (like invalidate_method_instance, but for our cache)
function invalidate(cache::CodeCache, replaced::MethodInstance, max_world, depth=0)
    cis = get(cache.dict, replaced, nothing)
    if cis === nothing
        return
    end
    for ci in cis
        if ci.max_world == ~0 % Csize_t
            @assert ci.min_world - 1 <= max_world "attempting to set illogical constraints"
            ci.max_world = max_world
        end
        @assert ci.max_world <= max_world
    end

    # recurse to all backedges to update their valid range also
    if isdefined(replaced, :backedges)
        backedges = replaced.backedges
        # Don't touch/empty backedges `invalidate_method_instance` in C will do that later
        # replaced.backedges = Any[]

        for mi in backedges
            invalidate(cache, mi, max_world, depth + 1)
        end
    end
end

const CI_CACHE = CodeCache()

## interpreter

using Core.Compiler: AbstractInterpreter, InferenceResult, InferenceParams, InferenceState, OptimizationParams

struct GPUInterpreter <: AbstractInterpreter
    # Cache of inference results for this particular interpreter
    cache::Vector{InferenceResult}
    # The world age we're working inside of
    world::UInt

    # Parameters for inference and optimization
    inf_params::InferenceParams
    opt_params::OptimizationParams

    function GPUInterpreter(world::UInt)
        @assert world <= Base.get_world_counter()

        return new(
            # Initially empty cache
            Vector{InferenceResult}(),

            # world age counter
            world,

            # parameters for inference and optimization
            InferenceParams(unoptimize_throw_blocks=false),
            OptimizationParams(unoptimize_throw_blocks=false),
        )
    end
end

# Quickly and easily satisfy the AbstractInterpreter API contract
Core.Compiler.get_world_counter(ni::GPUInterpreter) = ni.world
Core.Compiler.get_inference_cache(ni::GPUInterpreter) = ni.cache
Core.Compiler.InferenceParams(ni::GPUInterpreter) = ni.inf_params
Core.Compiler.OptimizationParams(ni::GPUInterpreter) = ni.opt_params
Core.Compiler.may_optimize(ni::GPUInterpreter) = true
Core.Compiler.may_compress(ni::GPUInterpreter) = true
Core.Compiler.may_discard_trees(ni::GPUInterpreter) = true
Core.Compiler.add_remark!(ni::GPUInterpreter, sv::InferenceState, msg) = nothing # TODO

## world view of the cache

using Core.Compiler: WorldView

function Core.Compiler.haskey(wvc::WorldView{CodeCache}, mi::MethodInstance)
    Core.Compiler.get(wvc, mi, nothing) !== nothing
end

function Core.Compiler.get(wvc::WorldView{CodeCache}, mi::MethodInstance, default)
    cache = wvc.cache
    for ci in get!(cache.dict, mi, CodeInstance[])
        if ci.min_world <= wvc.worlds.min_world && wvc.worlds.max_world <= ci.max_world
            # TODO: if (code && (code == jl_nothing || jl_ir_flag_inferred((jl_array_t*)code)))
            return ci
        end
    end

    return default
end

function Core.Compiler.getindex(wvc::WorldView{CodeCache}, mi::MethodInstance)
    r = Core.Compiler.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function Core.Compiler.setindex!(wvc::WorldView{CodeCache}, ci::CodeInstance, mi::MethodInstance)
    Core.Compiler.setindex!(wvc.cache, ci, mi)
end

## codegen/inference integration

Core.Compiler.code_cache(ni::GPUInterpreter) = WorldView(CI_CACHE, ni.world)

# No need to do any locking since we're not putting our results into the runtime cache
Core.Compiler.lock_mi_inference(ni::GPUInterpreter, mi::MethodInstance) = nothing
Core.Compiler.unlock_mi_inference(ni::GPUInterpreter, mi::MethodInstance) = nothing

function ci_cache_populate(mi, min_world, max_world)
    interp = GPUInterpreter(min_world)
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)

    # inference populates the cache, so we don't need to jl_get_method_inferred
    wvc = WorldView(CI_CACHE, min_world, max_world)
    @assert Core.Compiler.haskey(wvc, mi)

    # if src is rettyp_const, the codeinfo won't cache ci.inferred
    # (because it is normally not supposed to be used ever again).
    # to avoid the need to re-infer, set that field here.
    ci = Core.Compiler.getindex(wvc, mi)
    if ci !== nothing && ci.inferred === nothing
        ci.inferred = src
    end

    return
end

function ci_cache_lookup(mi, min_world, max_world)
    wvc = WorldView(CI_CACHE, min_world, max_world)
    return Core.Compiler.get(wvc, mi, nothing)
end

## interface

function compile_method_instance(@nospecialize(job::CompilerJob), method_instance::MethodInstance, world)
    # set-up the compiler interface
    debug_info_kind = if Base.JLOptions().debug_level == 0
        LLVM.API.LLVMDebugEmissionKindNoDebug
    elseif Base.JLOptions().debug_level == 1
        LLVM.API.LLVMDebugEmissionKindLineTablesOnly
    elseif Base.JLOptions().debug_level >= 2
        LLVM.API.LLVMDebugEmissionKindFullDebug
    end
    if job.target isa PTXCompilerTarget && !job.target.debuginfo
        debug_info_kind = LLVM.API.LLVMDebugEmissionKindNoDebug
    end
    params = Base.CodegenParams(;
        track_allocations  = false,
        code_coverage      = false,
        prefer_specsig     = true,
        gnu_pubnames       = false,
        debug_info_kind    = Cint(debug_info_kind),
        lookup             = @cfunction(ci_cache_lookup, Any, (Any, UInt, UInt)))

    # popoulate the cache
    if ci_cache_lookup(method_instance, world, world) === nothing
        ci_cache_populate(method_instance, world, world)
    end

    # generate IR
    native_code = ccall(:jl_create_native, Ptr{Cvoid},
                        (Vector{MethodInstance}, Base.CodegenParams, Cint),
                        [method_instance], params, #=extern policy=# 1)
    @assert native_code != C_NULL
    llvm_mod_ref = ccall(:jl_get_llvm_module, LLVM.API.LLVMModuleRef,
                         (Ptr{Cvoid},), native_code)
    @assert llvm_mod_ref != C_NULL
    llvm_mod = LLVM.Module(llvm_mod_ref)

    # get the top-level code
    code = ci_cache_lookup(method_instance, world, world)

    # get the top-level function index
    llvm_func_idx = Ref{Int32}(-1)
    llvm_specfunc_idx = Ref{Int32}(-1)
    ccall(:jl_get_function_id, Nothing,
          (Ptr{Cvoid}, Any, Ptr{Int32}, Ptr{Int32}),
          native_code, code, llvm_func_idx, llvm_specfunc_idx)
    @assert llvm_func_idx[] != -1
    @assert llvm_specfunc_idx[] != -1

    # get the top-level function)
    llvm_func_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                          (Ptr{Cvoid}, UInt32), native_code, llvm_func_idx[]-1)
    @assert llvm_func_ref != C_NULL
    llvm_func = LLVM.Function(llvm_func_ref)
    llvm_specfunc_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                              (Ptr{Cvoid}, UInt32), native_code, llvm_specfunc_idx[]-1)
    @assert llvm_specfunc_ref != C_NULL
    llvm_specfunc = LLVM.Function(llvm_specfunc_ref)

    # configure the module
    triple!(llvm_mod, llvm_triple(job.target))
    if llvm_datalayout(job.target) !== nothing
        datalayout!(llvm_mod, llvm_datalayout(job.target))
    end

    return llvm_specfunc, llvm_mod
end

Base.empty!(cache::CodeCache) = empty!(cache.dict)


text:
function cluster_points(centers::Matrix{Float64}, grd::Grid)

    Nc = size(centers)[2]
    M = zeros(2, grd.Nk, Nc)
    J = zeros(Int64, (grd.Nk, Nc))
    ks = zeros(Int64, Nc)

    for i=1:grd.Nk
        j = sortperm([norm(grd.points[:,i] - centers[:,k]) for k=1:Nc])[1]
        ks[j] += 1
        M[:,ks[j],j] = grd.points[:,i]
        J[ks[j],j] = i
    end

    points = zeros(2, grd.Nk)
    is = Int64[1:grd.Nk...]
    js = zeros(Int64, grd.Nk)
    dat = ones(Int64, grd.Nk)
    k = 0
    for i=1:Nc
        for j=1:ks[i]
            k += 1
            points[:,k] = M[:,j,i]
            js[k] = J[j,i]
        end
    end

    return points, sparse(is, js, dat, grd.Nk, grd.Nk)
end

function sort_by_angle(grd::Grid, r0)

    thetas = Float64[atan(grd.points[2,i]-r0[2], grd.points[1,i]-r0[1]) for i=1:grd.Nk]

    is = Int64[1:grd.Nk...]
    js = sortperm(thetas)
    dat = ones(Int64, grd.Nk)

    return grd.points[:,js], sparse(is, js, dat, grd.Nk, grd.Nk)
end


text:
using LAR

V = [rand(2) for k=1:10000]
p.VIEW(p.STRUCT(p.AA(p.MK)(V)))


text:
module ImageContrast

# using Graphics
using Cairo
using Tk
using Winston
using Images
import Colors.Fractional
using Compat; import Compat.String

type ContrastSettings
    min
    max
end

type ContrastData
    imgmin
    imgmax
    phist::FramedPlot
    chist::Canvas
end

# The callback should have the syntax:
#    callback(cs)
# The callback's job is to replot the image with the new contrast settings
function contrastgui{T}(img::AbstractArray{T}, cs::ContrastSettings, callback::Function)
    win = Toplevel("Adjust contrast", 500, 300, true)
    contrastgui(win, img, cs, callback)
end

function contrastgui{T}(win::Tk.TTk_Container, img::AbstractArray{T}, cs::ContrastSettings, callback::Function)
    # Get initial values
    dat = img[:,:]
    immin = minfinite(dat)
    immax = maxfinite(dat)
    if is(cs.min, nothing)
        cs.min = immin
    end
    if is(cs.max, nothing)
        cs.max = immax
    end
    cs.min = oftype(immin, cs.min)
    cs.max = oftype(immax, cs.max)

    # Set up GUI
    fwin = Frame(win)
    w = width(win.w)
    h = height(win.w)
    pack(fwin, expand=true, fill="both")

    dimg = immax-immin
    rmin = isa(immin, Fractional) ? immin - dimg/10^4 : floor(Integer, immin)
    rmax = isa(immax, Fractional) ? immax + dimg/10^4 : ceil(Integer, immax)
    max_slider = Slider(fwin, rmin, rmax)
    set_value(max_slider, cs.max)
    chist = Canvas(fwin, 2w/3, h)
    min_slider = Slider(fwin, rmin, rmax)
    set_value(min_slider, cs.min)

    Tk.grid(max_slider, 1, 1, sticky="ew", padx=5)
    Tk.grid(chist, 2, 1, sticky="nsew", padx=5)
    Tk.grid(min_slider, 3, 1, sticky="ew", padx=5)
    grid_columnconfigure(fwin, 1, weight=1)
    grid_rowconfigure(fwin, 2, weight=1)

    emax = Entry(fwin, width=10)
    emin = Entry(fwin, width=10)
    set_value(emax, string(float32(cs.max)))
    set_value(emin, string(float32(cs.min)))
#    emax[:textvariable] = max_slider[:variable]
#    emin[:textvariable] = min_slider[:variable]

    fbuttons = Frame(fwin)
    zoom = Button(fbuttons, "Zoom")
    full = Button(fbuttons, "Full range")
    Tk.grid(zoom, 1, 1, sticky="we")
    Tk.grid(full, 2, 1, sticky="we")

    Tk.grid(emax, 1, 2, sticky="nw")
    Tk.grid(fbuttons, 2, 2, sticky="nw")
    Tk.grid(emin, 3, 2, sticky="nw")

    # Prepare the histogram
    nbins = ceil(Int, min(sqrt(length(img)), 200))
    p = prepare_histogram(dat, nbins, immin, immax)

    # Store data we'll need for updating
    cdata = ContrastData(immin, immax, p, chist)

    function rerender()
        pcopy = deepcopy(cdata.phist)
        bb = Winston.limits(cdata.phist)
        add(pcopy, Curve([cs.min, cs.max], [bb.y0, bb.y1], "linewidth", 10, "color", "white"))
        add(pcopy, Curve([cs.min, cs.max], [bb.y0, bb.y1], "linewidth", 5, "color", "black"))
        Winston.display(chist, pcopy)
        reveal(chist)
        callback(cs)
        Tk.update()
    end
    # If we have a image sequence, we might need to generate a new histogram.
    # So this function will be returned to the caller
    function replaceimage(newimg, minval = min(newimg), maxval = max(newimg))
        p = prepare_histogram(newimg, nbins, minval, maxval)
        cdata.imgmin = minval
        cdata.imgmax = maxval
        cdata.phist = p
        rerender()
    end

    # Set initial histogram scale
    setrange(cdata.chist, cdata.phist, cdata.imgmin, cdata.imgmax, rerender)

    # All bindings
    bind(emin, "<Return>") do path
        try
            my_min = float32(get_value(emin))
            my_max = float32(get_value(emax))
            # Don't let values cross
            my_max = my_max < my_min ? my_min + 0.01f0 : my_max # offset is arbitrary
            cs.min = convertsafely(typeof(cs.min), my_min)
            cs.max = convertsafely(typeof(cs.max), my_max)
            set_value(emax, string(my_max)) # caution: each widget gets it's own type
            set_value(min_slider, my_min)
            set_value(max_slider, my_max)
            rerender()
        catch
            set_value(emin, string(cs.min))
        end
    end
    bind(emax, "<Return>") do path
        try
            my_min = float32(get_value(emin))
            my_max = float32(get_value(emax))
            # Don't let values cross
            my_min = my_min > my_max ? my_max - 0.01f0 : my_min # offset is arbitrary
            cs.min = convertsafely(typeof(cs.min), my_min)
            cs.max = convertsafely(typeof(cs.max), my_max)
            set_value(emin, string(my_min))
            set_value(min_slider, my_min)
            set_value(max_slider, my_max)
            rerender()
        catch
            set_value(emax, string(cs.max))
        end
    end
    bind(min_slider, "command") do path
        my_min = float32(min_slider[:value])
        my_max = float32(max_slider[:value])
        # Don't let values cross
        my_max = my_max < my_min ? my_min + 0.01f0 : my_max # offset is arbitrary
        cs.min = convertsafely(typeof(cs.min), my_min)
        cs.max = convertsafely(typeof(cs.max), my_max)
        set_value(max_slider, my_max)
        set_value(emax, string(my_max))
        set_value(emin, string(my_min))
        try
            rerender()
        catch
            set_value(max_slider, my_max)
        end
    end
    bind(max_slider, "command") do path
        my_min = float32(min_slider[:value])
        my_max = float32(max_slider[:value])
        # Don't let values cross
        my_min = my_min > my_max ? my_max - 0.01f0 : my_min # offset is arbitrary
        cs.min = convertsafely(typeof(cs.min), my_min)
        cs.max = convertsafely(typeof(cs.max), my_max)
        set_value(min_slider, my_min)
        set_value(emin, string(my_min))
        set_value(emax, string(my_max))
        try
            rerender()
        catch
            set_value(min_slider, my_min)
        end
    end
    bind(zoom, "command", path -> setrange(cdata.chist, cdata.phist, cdata.imgmin, cdata.imgmax, rerender))
    bind(full, "command", path -> setrange(cdata.chist, cdata.phist, min(cdata.imgmin, cs.min), max(cdata.imgmax, cs.max), rerender))

    replaceimage
end

convertsafely{T<:Integer}(::Type{T}, val) = convert(T, round(val))
convertsafely{T}(::Type{T}, val) = convert(T, val)

function prepare_histogram(img, nbins, immin, immax)
    w = (immax-immin)/10^5
    e = immin-w:(immax-immin)/(nbins-1):immax+w
    dat = img[:]
    e, counts = hist(dat[isfinite(dat)], e)
    counts .+= 1   # because of log scaling
    x, y = stairs(e, counts)
    p = FramedPlot()
    setattr(p, "ylog", true)
    setattr(p.y, "draw_nothing", true)
    setattr(p.x2, "draw_nothing", true)
    setattr(p.frame, "tickdir", 1)
    add(p, FillBetween(x, ones(length(x)), x, y, "color", "black"))
    p
end

function stairs(xin::AbstractVector, yin::Vector)
    nbins = length(yin)
    if length(xin) != nbins+1
        error("Pass edges for x, and bin values for y")
    end
    xout = zeros(0)
    yout = zeros(0)
    sizehint!(xout, 2nbins)
    sizehint!(yout, 2nbins)
    push!(xout, xin[1])
    for i = 2:nbins
        xtmp = xin[i]
        push!(xout, xtmp)
        push!(xout, xtmp)
    end
    push!(xout, xin[end])
    for i = 1:nbins
        ytmp = yin[i]
        push!(yout, ytmp)
        push!(yout, ytmp)
    end
    xout, yout
end

function setrange(c::Canvas, p, minval, maxval, render::Function)
    setattr(p, "xrange", (minval, maxval))
    render()
end

end


text:
using Gym
using StatsBase: sample, Weights
using DataStructures: CircularBuffer
using Statistics: mean

env = make("Copy-v0")
reset!(env)

max_episodes = 75000
α = 3f-1
γ = 9f-1
ϵ = 1f-1
ϵ_schedule = 10000
# Copy is considered solved when the average reward for 100 continuous episodes >=25
goal = 25f0

num_dims = [length(space) for space in env.action_space.spaces]
num_actions = prod(num_dims)
num_states = length(env.observation_space)

Q = zeros(Float32, num_states, num_actions)
P = zeros(Float32, num_actions)

function policy(state)
    fill!(P, ϵ/num_actions)
    P[argmax(Q[state[1], :])] += 1 - ϵ
    return sample(1:num_actions, Weights(P))
end

function actioninttotup(int)
  @assert int <= prod(num_dims)
  res = []
  int -= 1
  for dim in reverse(num_dims)
    push!(res, (int%dim)+1)
    int = div(int, dim)
  end
  return reverse(res)
end

total_reward = CircularBuffer{Float32}(100)
for e in 1:max_episodes
    global α, γ, ϵ_schedule, ϵ, goal
    if e % 50 == 0
        println("Episode: $e | Average reward of the last 100 episodes: $(mean(total_reward))")
        if mean(total_reward) > goal
            println("Copy Solved  at episode = $e")
            break
        end
    end
    e % ϵ_schedule == 0 && (ϵ /= 2)
    current_state = reset!(env)
    done = false
    push!(total_reward, 0f0)
    while !done
        a = policy(current_state)
        decoded_a = actioninttotup(a)
        next_state, reward, done, _ = step!(env, decoded_a)
        #println("a = $a\tcurrent_state[1] = $(current_state[1])\tnext_state[1] = $(next_state[1])")
        Q[current_state[1], a] += α * ((reward + γ * maximum(Q[next_state[1], :])) - Q[current_state[1], a])
        current_state = next_state
        total_reward[end] += reward
    end
end


text:
####################################################################################
####################################################################################
############################       CHECK RHS      ##################################
####################################################################################
####################################################################################

function check_lhs_rhs(model)
    new_model = deepcopy(model)
    lhs_states, lhs_parameters, lhs_forcings, lhs_constants, lhs_equations = get_lhs(new_model)
    all_rhs, rhs_to_lhs, lhs_to_rhs = get_rhs(new_model.Equations)
    for i in all_rhs
        i ∉ [lhs_states, lhs_parameters, lhs_forcings, lhs_constants, lhs_equations] &&
          error("The input $i is missing from the model")
    end
    for i in lhs_parameters
        if i ∉ all_rhs
            warn("The output is not dependent on the parameter $i. It has been removed from the model.")
            delete!(new_model.Parameters, i)
        end
    end
    for i in lhs_constants
        if i ∉ all_rhs
            warn("The output is not dependent on the constant $i. It has been removed from the model.")
            delete!(new_model.Constants, i)
        end
    end
    for i in lhs_forcings
        if i ∉ all_rhs
            warn("The output is not dependent on the forcing $i. It has been removed from the model.")
            delete!(new_model.Forcings, i)
        end
    end
    for i in lhs_states
        if i ∉ all_rhs
            warn("The output is not dependent on the state $i, but it has NOT been removed.")
        end
    end
    for i in lhs_equations
        if (i ∉ all_rhs) && (i ∉ [model.NamesObserved, model.NamesDerivatives])
            #warn("The output is not dependent on the intermediate variable $i. It has been removed from the model.")
            delete!(new_model.Equations, i)
        end
    end
    return new_model, rhs_to_lhs, lhs_to_rhs
end

####################################################################################
####################################################################################
###########################       CHECK UNITS      #################################
####################################################################################
####################################################################################

function calculate_units(exp::Expr, units::Dict{Symbol, dt.Dimension}, c::Int64)
  ex = deepcopy(exp)
  for i in 1:length(ex.args)
    if isa(ex.args[i], Expr)
      ex.args[i], c = calculate_units(ex.args[i], units, c)
    elseif isa(ex.args[i], Symbol) && haskey(units, symbol(ex.args[i]))
      c += 1
      name = parse("__unit__$(c)")
      @eval global $name = ($(units[ex.args[i]]))
      ex.args[i] = name
    end
  end
  ex, c
end

function check_units(sorted_model::dt.OdeSorted)
  given_units = Dict{Symbol, dt.Dimension}()
  for (key,val) in sorted_model.Constants  given_units[symbol(key)] = val.Units.d end
  for (key,val) in sorted_model.Parameters  given_units[symbol(key)] = val.Units.d end
  for (key,val) in sorted_model.States  given_units[symbol(key)] = val.Units.d end
  for (key,val) in sorted_model.Forcings  given_units[symbol(key)] = val.Units.d end
  for level in sorted_model.SortedEquations
    for (key,val) in level  given_units[symbol(key)] = val.Units.d end
  end
  # Add units of time (as input that is always available)
  given_units[symbol("time")] = dt.s.d
  for level in sorted_model.SortedEquations[2:end]
    for (key,val) in level
      expected = given_units[symbol(key)]
      if isa(val.Expr, Symbol)
            infered = eval(given_units[val.Expr])
      else
        infered_expr, c = calculate_units(val.Expr, given_units, 0)
        infered = try eval(infered_expr) catch
            error("Error when calculating units: The rhs of equation $key is not dimensionally homogeneous.
The right hand side expression was $(val.Expr)") end
      end
      expected != infered && error("Error when calculating units: Expected and given units for $key did not coincide.
I infered the units $infered but you assign it $expected.
The right hand side expression was $(val.Expr)")
    end
  end
end

####################################################################################
####################################################################################
########################       GENERATE JACOBIAN      ##############################
####################################################################################
####################################################################################

# Calculation Jacobian matrix of the model
function generate_jacobian_matrix(compressed_model::dt.OdeSorted)
  names_states = collect(keys(compressed_model.States))
  Jacobian = Array(Union(Expr, Symbol, Number),(length(names_states), length(names_states)))
  cs = 1
    for j in names_states
      cd = 1
      for i in compressed_model.NamesDerivatives
        @inbounds Jacobian[cd,cs] =  differentiate(compressed_model.SortedEquations[2][i].Expr, parse(j))
        cd += 1
      end
    cs += 1
  end
  return Jacobian
end

####################################################################################
####################################################################################
###################       GENERATE EXTENDED SYSTEM      ############################
####################################################################################
####################################################################################

# Calculate array of sensitivities
function generate_sensitivity_array(compressed_model::dt.OdeSorted)
  names_states = collect(keys(compressed_model.States))
  names_derivatives = ["d_"*x*"_dt" for x in names_states]
  Sensitivity = Array{Union(Expr, Symbol, Number), 1}[]
  names_parameters = collect(keys(compressed_model.Parameters))
  for i in names_parameters
    sens = Array(Union(Expr, Symbol, Number), length(names_states))
    c = 1
    for j in names_derivatives
      sens[c] = differentiate(compressed_model.SortedEquations[2][j].Expr, parse(i))
      c += 1
    end
    push!(Sensitivity, sens)
  end
  return Sensitivity
end


text:
function (mowa::MovingWindowSolver)(solution::MovingWindowSolution, problem::AbstractInitialValueProblem; saveiterates::Bool=false)
    @↓ parallelsolver, budget, adaptive, τ, Δτ = mowa
    @↓ finesolver, coarsolver, tolerance, P, K = parallelsolver
    K₀ = K
    @↓ finestepsize ← stepsize = finesolver
    @↓ hF ← h = finestepsize
    hF₀ = hF
    @↓ coarstepsize ← stepsize = coarsolver
    @↓ hG ← h = coarstepsize
    hG₀ = hG
    @↓ ϵ = tolerance
    if budget !== nothing
        @↓ Kmax = budget
        @↑ parallelsolver = K ← Kmax
    end
    if adaptive !== nothing
        @↓ δτ⁻, δτ⁺, δΔτ⁻, δΔτ⁺, Nᵣ, fineupdate = adaptive
    end
    @↓ windows, restarts = solution
    @↓ u0, (t0, tN) ← tspan = problem

    τ0 = t0
    τN = τ0 + τ
    τJ = zero(τ0)
    # τ0, τN, τJ = windowparams(problem)

    cache = TimeParallel.TimeParallelCache(problem, parallelsolver)
    @↓ U, G, T = cache

    m = 1
    M = M₀ = length(windows)
    counter = counter2 = 0
    while τ0 < tN
        if τN ≥ tN
            counter += 1
        end
        # if still not done, δτ⁺ease solution.windows length
        if m > M
            append!(windows, Vector{AbstractTimeParallelSolution}(undef, M₀))
            append!(restarts, zeros(M₀))
            M += M₀
        end
        # start with a coarse guess
        if m == 1
            # define windowproblem
            windowproblem = subproblemof(problem, u0, τ0, τN)
            # update cache
            coarseguess!(cache, windowproblem, parallelsolver)
        elseif m > 1
            𝑢 = restarts[m] == 0 ? windows[m-1] : windows[m]
            # recompute T
            for n = 1:P+1
                T[n] = (P - n + 1) / P * τ0 + (n - 1) / P * τN
            end
            # recompute U
            for n = 1:P
                # transfer initial conditions from (m-1)-th window (or m-th window at previous "restart")
                if T[n] ≤ τJ
                    U[n] = 𝑢(T[n])
                    # compute G for Parareal
                    if n == 1
                        G[n] = U[n]
                    else
                        windowchunkproblem = subproblemof(problem, U[n-1], T[n-1], T[n])
                        G[n] = coarsolver(windowchunkproblem).u[end]
                    end
                # compute the remaining initial conditions with coarse solver
                else
                    windowchunkproblem = subproblemof(problem, U[n-1], T[n-1], T[n])
                    G[n] = U[n] = coarsolver(windowchunkproblem).u[end]
                end
            end
            windowproblem = subproblemof(problem, U[1], τ0, τN)
            # @↑ cache = U, G, T
        end
        windows[m] = parallelsolver(cache, windowproblem; saveiterates=saveiterates)
        # println("$(restarts[m]) restarts, $(length(solution[m].errors)) Parareal iterations")
        # println("——————————————————————————————————")
        # save τJ for next window
        τJ = τN
        # adaptive check
        if adaptive !== nothing
            if (windows[m].errors[end] > ϵ) && (restarts[m] < Nᵣ) && (hG > hF)
                # δτ⁻₂ = max(δτ⁻, P * hG / τ) # P * hG is minimum possible
                # τ *= δτ⁻₂
                # Δτ *= δτ⁻₂
                # if fineupdate
                #     hF *= δτ⁻₂
                # end
                # hG = max(hG * δτ⁻₂, hF)
                Δτ = max(δΔτ⁻ * Δτ, τ / P)
                restarts[m] += 1
                counter2 += 1
            else
                counter2 = 0
                # τ *= δτ⁺
                # Δτ *= δτ⁺
                # if fineupdate
                #     hF *= δτ⁺
                # end
                # hG *= δτ⁺
                Δτ = min(δΔτ⁺ * Δτ, τ) # avoids overshoot
                τ0 = min(τ0 + Δτ, τJ)
                m += 1
            end
            if fineupdate
                @↑ finestepsize = h ← hF
            end
            @↑ coarstepsize = h ← hG
        else
            τ0 += Δτ
            m += 1
        end
        τN = τ0 + τ
        if counter ≥ 1 && counter2 == 0
            break
        end
    end
    m -= 1 # to correct m += 1 at end of loop
    @↑ parallelsolver = K ← K₀
    @↑ finestepsize = h ← hF₀
    @↑ coarstepsize = h ← hG₀
    resize!(windows, m)
    resize!(restarts, m)
    return solution
end

function (mowa::MovingWindowSolver)(problem::AbstractInitialValueProblem; saveiterates::Bool=false)
    solution = MovingWindowSolution(problem, mowa)
    mowa(solution, problem; saveiterates=saveiterates)
    return solution
end


text:
# main program structure construction

function solveMain(T, Δt, N, allGen = false, qpopt = false, iterMin = 100,
    iterMax = 1000, cutDict = Dict(), ubGen = false, ubM = 200, hardened = [], simuRule = 0, simuSample = Dict())
    # readin data and execute the SDDP algorithm

    UB = 9999999999;
    LB = -9999999999;
    iterNo = 0;
    keepIter = true;
    LBHist = [];
    UBHist = [];
    UBuHist = [];
    UBlHist = [];
    timeHist = [];
    # initialize the cutDict
    for j in procs()
        remotecall_fetch(cutIni,j,cutDict);
    end
    sampleCounter = 0;

    # while the termination criteria is not met
    while (keepIter)&(iterNo <= iterMax)
        iterNo += 1;
        iterStart = time();
        if simuRule == 0
            if simuSample == Dict()
                trialPaths,currentLB,currentUBDict = exeForward(T, Δt, N, qpopt, Dict(), hardened);
            else
                pathIter = Dict();
                for sNo in 1:N
                    pathIter[sNo] = simuSample[sampleCounter+sNo];
                end
                trialPaths,currentLB,currentUBDict = exeForward(T, Δt, N, qpopt, pathIter, hardened);
                sampleCounter += N;
            end
        else
            trialPaths,currentLB,currentUBDict = exeForward_simuOpt(T, Δt, N, iterNo, qpopt, hardened);
        end
        push!(LBHist,currentLB);

        # forward pass: obtain the trial paths
        if !(ubGen)
            currentUBList = [currentUBDict[ubkey] for ubkey in keys(currentUBDict)];
        else
            trialPathsUB,currentLBUB,currentUBDict = exeForward(T, Δt, ubM, qpopt, Dict(), hardened);
            currentUBList = [currentUBDict[ubkey] for ubkey in keys(currentUBDict)];
        end
        push!(UBHist,mean(currentUBList));
        currentUBu = mean(currentUBList) + 1.96*std(currentUBList);
        push!(UBuHist,currentUBu);
        currentUBl = mean(currentUBList) - 1.96*std(currentUBList);
        push!(UBlHist,currentUBl);

        if (currentLB <= currentUBu)&(currentLB >= currentUBl)&(iterNo >= iterMin)
            keepIter = false;
        else
            if allGen
                exeBackwardAll(T, Δt, trialPaths, qpopt, hardened);
            else
                exeBackward(T, Δt, trialPaths, qpopt, hardened);
            end
        end
        iterElapsed = time() - iterStart;
        push!(timeHist,iterElapsed);
        println("========= Iteration $(iterNo) Finished, LB = $(round(currentLB, digits = 2)),
            UB = [$(round(currentUBl,digits = 2)),$(round(currentUBu, digits = 2))], Time = $(iterElapsed) sec. =========")
    end
    return cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist;
end

function preGen(T, Δt, N, iterMax, qpopt = false, cutDict = Dict(), hardened = [])
    # pregenerate cuts

    UB = 9999999999;
    LB = -9999999999;
    iterNo = 0;
    keepIter = true;
    LBHist = [];
    UBHist = [];
    UBuHist = [];
    UBlHist = [];
    timeHist = [];
    # initialize the cutDict
    for j in procs()
        remotecall_fetch(cutIni,j,cutDict);
    end

    while iterNo <= iterMax
        iterNo += 1;
        trialPaths,currentLB,currentUBDict = exeForward_last(T, Δt, N, qpopt, hardened);
        exeBackward_last(T, Δt, trialPaths, qpopt, hardened);
    end

    return cutDict;
end