text:
@testset "Result" begin
    @testset "Result+PathResult" begin
        R = solve(
            equations(heart()),
            save_all_paths = true,
            seed = 13186,
            show_progress = false,
        )
        @test R isa Result
        @test nat_infinity(R) ≤ 572
        @test nfinite(R) == 4
        @test length(collect(R)) == 576
        @test finite(R) isa Vector{<:PathResult}
        @test seed(R) == 13186

        @test length(finite(R, only_nonsingular = false)) == 4
        @test length(finite(R, only_nonsingular = true)) == 4
        @test length(finite(R, only_singular = true)) == 0
        @test 572 - length(failed(R)) == nat_infinity(R)
        @test length(real(R, tol = 1e-6)) == 2
        @test nreal(R, tol = 1e-6) == 2
        @test length(real_solutions(R)) == 2
        @test_deprecated realsolutions(R)
        @test length(at_infinity(R)) ≤ 572
        @test length(results(R, only_real = true, real_tol = 1e-8)) == 2
        @test length(results(R, only_nonsingular = true, singular_tol = 1e9)) == 4
        @test length(finite(results(R, only_real = true))) == 2
        @test nresults(R, only_nonsingular = true, singular_tol = 1e9) == 4
        @test length(results(R, onlyfinite = false)) == 576
        @test nresults(R, onlyfinite = false) == 576
        @test nnonsingular(R) == 4
        @test length(nonsingular(R)) == 4
        @test length(singular(R)) == 0
        @test length(singular(real(R))) == 0

        @test_nowarn mapresults(solution, R)
        @test_nowarn mapresults(start_solution, R)
        @test_nowarn mapresults(residual, R)
        @test_nowarn mapresults(is_singular, R)
        @test count(mapresults(is_real, R)) == 2
        # test fallback
        @test count(results(is_real, R)) == 2

        @test length(solutions(R)) == 4
        @test solutions(R) isa Vector{Vector{ComplexF64}}
        @test real_solutions(R, real_tol = 1e-8) isa Vector{Vector{Float64}}
        @test length(real_solutions(R, real_tol = 1e-8)) == 2

        @test_nowarn string(R)
        @test_nowarn string(R[end])
        @test_nowarn show(IOContext(devnull, :compact => true), R)
        @test_nowarn show(IOContext(devnull, :compact => true), R[end])

        test_treeviews(R)

        @polyvar x y z
        R = solve([(x - 3)^3, (y - 2)], affine_tracking = true)
        @test R isa Result
        test_treeviews(R)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 1
        @test nsingular(R, counting_multiplicities = true) == 3
        @test_nowarn sprint(show, R)

        R = solve([(x - 3)^3, (y - 2)], affine_tracking = true, system_scaling = nothing)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 1
        @test nsingular(R, counting_multiplicities = true) == 3
        @test_nowarn sprint(show, R)

        @polyvar x y z
        R = solve([(x - 3z), (y - 2z)])
        @test R isa Result{<:ProjectiveVectors.PVector}
        @test_nowarn string(R)
        test_treeviews(R)
        @test nnonsingular(R) == 1

        @polyvar x y z
        R = solve([(x - 3z)^3, (y - 2z)])
        @test R isa Result{<:ProjectiveVectors.PVector}
        test_treeviews(R)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 1
        @test nsingular(R, counting_multiplicities = true) == 3

        @polyvar x y z
        V = x^2 + y^2 - z^2
        L = randn(1, 3) * [x, y, z]
        R = solve([V; L])
        @test_nowarn sprint(show, R)
    end

    @testset "singular result" begin
        @polyvar x
        f = (x - 3)^3 * (x - 2)
        R = solve([f])
        @test length(solutions(R)) == 2
        @test sort(results(multiplicity, R)) == [1, 3]

        @polyvar x y
        f = (x - 3y)^3 * (x - 2y)
        R = solve([f])
        @test nsolutions(R) == 2
        @test sort(results(multiplicity, R)) == [1, 3]
    end

    @testset "multiplicities" begin
        @polyvar x y z
        f = (x - 1 * y) * (x - 1.01 * y)
        g = x - y - z
        S = solve([f, g], threading = false)
        @test length(multiplicities(S, tol = 1e-1)) == 1
        @test length(multiplicities(S, tol = 1e-5)) == 0
        @test length(multiplicities(S.pathresults, tol = 1e-1)) == 1
        @test length(multiplicities(S.pathresults, tol = 1e-5)) == 0

        @polyvar x z
        f = (x - 1) * (x - 1.01)
        g = x - 1 - z
        S = solve([f, g], threading = false)
        @test length(multiplicities(S, tol = 1e-1)) == 1
        @test length(multiplicities(S, tol = 1e-5)) == 0
        @test length(multiplicities(S.pathresults, tol = 1e-1)) == 1
        @test length(multiplicities(S.pathresults, tol = 1e-5)) == 0
        @test norm(solution(S[1]) - solution(S[2])) < 1e-1
    end

    @testset "real projective solutions" begin
        @polyvar x y v w
        res = solve([x * y - 6 * v * w, x^2 - 5 * v^2], variable_groups = [(x, v), (y, w)])
        @test all(is_real, res)
        @test sum(fubini_study.(real_solutions(res), solutions(res))) ≈ 0.0 atol = 1e-6
    end
end


text:
@testset "Testing Main deconvolution code" begin
    Random.seed!(42)

    img = abs.(randn((3, 3)))
    psf = abs.(randn((3, 3)))

    res = [4.933040319558362 1.640892951656587e-9 1.941347101658743e-8; 1.4783685597711925e-8 4.723499796776938e-7 3.668727889064141e-9; 2.765081805017183e-9 3.280055085010056e-9 2.583774989955379]

    #= @show deconvolution(img, psf, λ=0.01)[1] =#
    @test all(≈(res, deconvolution(img, psf, λ=0.01)[1], rtol=0.1))
    @test all(≈(res, deconvolution(img, psf, plan_fft=false, λ=0.01)[1], rtol=0.1))
    @test all(≈(res, deconvolution(img, psf, plan_fft=false, λ=0.01, iterations=20)[1], rtol=0.1))

    # testing regularizer
    res2 = [4.188270526990536 5.999388400251461e-10 2.8299849680327642e-8; 1.725273124171714e-7 2.54195544512864 2.0216187854619135e-9; 9.594324085846738e-10 1.2000166997002865e-8 0.7863126081711094]
    @test all(≈(res2, deconvolution(img, psf, regularizer=nothing, padding=0.0)[1], rtol=0.1))

    # testing padding
    img = abs.(20 .* randn((3, 3)))
    psf = abs.(randn((3, 3)))
    #= @show deconvolution(img, psf, regularizer=nothing, padding=0.1)[1] =#
    res3 = [2.7358233029549654e-13 91.81808945246226 4.8779432830104286e-17; 2.5136104031366187e-13 89.00977286151101 3.0110340818892347e-13; 3.244894460054766e-12 18.817221108658277 89.77201822876616]
    @test all(≈(res3, deconvolution(img, psf, regularizer=nothing, padding=0.1)[1], rtol=1e-2))

    # test without mapping
    res4 = [-5.332320754919289 40.63633433786456 -18.173096688407522; -15.597390722842547 59.01033285248438 -4.122592175577954; 13.402949467681973 6.4622595721116465 42.82143693411391]
    #= @show deconvolution(img, psf, regularizer=nothing, padding=0.1, mapping=nothing)[1] =#
    @test all(≈(res4, deconvolution(img, psf, regularizer=nothing, padding=0.1, mapping=nothing)[1], rtol=0.1))

    # test broadcasting with image having more dimensions
    img = zeros((3, 3, 2))
    imgc = abs.(randn((3, 3, 1)))
    img[:, :, 1] = imgc
    img[:, :, 2] = imgc
    psf = zeros((3, 3))
    psf[1,1] = 1

    res = deconvolution(img, psf, regularizer=GR(num_dims=3, sum_dims=[1,2]))[1]
    @test all(res[:, :, 1] .≈ res[:, :, 2])

    img = zeros((3, 3, 2, 1))
    imgc = abs.(randn((3, 3, 1, 1)))
    img[:, :, 1, 1] = imgc
    img[:, :, 2, 1] = imgc
    psf = zeros((3, 3))
    psf[1,1] = 1

    res = deconvolution(img, psf, regularizer=GR(num_dims=4, sum_dims=[1,2]))[1]
    @test all(≈(res[:, :, 1, :], res[:, :, 2, :], rtol=0.1))

end

ext: jl


text:
# --------------------------------------------------------------------------------------
# Serial connection and stepper motor initialization
# --------------------------------------------------------------------------------------

mutable struct StepperSystem{N}
    con
    pos::MVector{N,Float64}
    id::MVector{N,String}
    step2coord::MVector{N,Any}
    coord2step::MVector{N,Any}
    depend::MVector{N,Any}
end

"""
    stepper_open(dof::Int; port=nothing, baud=9600)

## Description
Function to open a generic connection and returns a StepperSystem type to control a system of stepper motors
## Arguments
- dof: degrees of freedom
- port: port path. If nothing it'll get the first port of SerialPort.list_serialports
- baud: baud rate (default = 9600)
- testnocon: true for testing purposes. Makes dev.con equals missing
## Observations
The output is a StepperSystem type including the following elements with fixed size (dof):
- con: port connection
- pos: position
- id: stepper motors IDs
- step2coord: array of functions to convert steps into coordinates
- coord2step: array of functions to convert coordinates into steps
- depend: attribute with the dependecies. To be used with the conversion functions. Allows dependent systems.
Except from con, all the other attributes must be configured (see StepperControl.stepper_config)
"""
function stepper_open(dof::Int; port=nothing, baud=9600, testnocon=false)

    if testnocon
        con = missing
    else
        if isnothing(port)
            port = list_serialports()[1]
        end
        con = SerialPort(port, baud)
    end

    pos = repeat([0.0], dof)
    id = "m" .* string.(1:dof)
    f(x) = float(x)
    g(x) = Int(round(x))

    dev = StepperSystem{dof}(con, pos, id, repeat([f], dof), repeat([g], dof), 1:dof)

    return dev
end

"""
    linear_step2coord(; spr::Int, r=1.0)

## Description
Linear function to find coordinates from a number of steps
## Arguments
- spr: steps per revolution
- r: bell crank radius. The unit can be used to express angular displacements, in radians.
## Example
```jldoctest
julia> x = linear_step2coord(spr=2048);

julia> x(512)
1.5707963267948966

julia> x(512)*180/π
90.0
```
"""
function linear_step2coord(;spr::Int, r=1.0)

    function f(steps::Real)
        steps = Int(round(steps))
        r*2π*steps/spr
    end

    return f
end

"""
    linear_coord2step(; spr::Int, r=1.0)

## Description
Linear function to find steps from a displacement
## Arguments
- spr: steps per revolution
- r: bell crank radius. The unit can be used to express angular displacements, in radians.
## Example
```jldoctest
julia> x = linear_coord2step(spr=2048);

julia> x(π/2)
512
```
"""
function linear_coord2step(;spr::Int, r=1.0)
    g(coord::Real) = Int(round(coord*spr/(r*2π)))
    return g
end

"""
    stepper_config!(dev::StepperSystem; motorID::AbstractVector=dev.id, step2coord::AbstractVector=dev.step2coord, coord2step::AbstractVector=dev.coord2step, depend::AbstractVector=dev.depend)

## Description
Configure stepper system.
## Arguments
- dev: element of StepperSystem type
- motorID: array with motor IDs
- step2coord: function to convert steps into coordinates
- coord2step: function to convert coordinates into steps
-
## Examples
```jldoctest
julia> r = stepper_open(2);

julia> r.id
2-element MArray{Tuple{2},String,1,2} with indices SOneTo(2):
 "m1"
 "m2"

julia> r.step2coord[1](10)
10
julia> r.coord2step[1](10)
10

julia> f = linear_step2coord(spr=2048, r=2);

julia> g = linear_coord2step(spr=2048, r=2);

julia> stepper_config!(r, motorID=["x","y"], step2coord=[f], coord2step=[g]);

julia> r.id
2-element MArray{Tuple{2},String,1,2} with indices SOneTo(2):
 "x"
 "y"

 julia> r.step2coord[1](10)
 0.06135923151542565
 julia> r.step2coord[1](2048)/(2π)
 2.0

 julia> r.coord2step[1](10)
 1630
 julia> r.coord2step[1](4π)
 2048
```
"""
function stepper_config!(dev::StepperSystem; motorID::AbstractVector = dev.id,
                        step2coord::AbstractVector = dev.step2coord, coord2step::AbstractVector = dev.coord2step,
                        depend::AbstractVector = dev.depend)

    n = length(dev.id)
    dev.id = motorID

    if length(step2coord) == 1
        step2coord = repeat(step2coord, n)
    end
    dev.step2coord = step2coord

    if length(coord2step) ==1
        coord2step = repeat(coord2step, n)
    end
    dev.coord2step = coord2step

    dev.depend = depend

end


text:
"""
    O = gd2gmt(dataset; band=0, bands=[], sds=0, pad=0)

Convert a GDAL raster dataset into either a GMTgrid (if type is Int16 or Float) or a GMTimage type
Use `band` to select a single band of the dataset. When you know that the dataset contains several
bands of an image, use the kwarg `bands` with a vector the wished bands. By default it reads all
bands of the image or grid object.

When DATASET is a string it may contain the file name or the name of a subdataset. In former case
you can use the kwarg `sds` to selec the subdataset numerically. Alternatively, provide the full `sds` name.
For files with `sds` with a scale_factor (e.g. MODIS data), that scale is applyied automaticaly.

    Examples:
       G = gd2gmt("AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc", sds=1);
    or
       G = gd2gmt("SUBDATASET_1_NAME=NETCDF:AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc:sst");
    or
       G = gd2gmt("NETCDF:AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc:sst");
"""
function gd2gmt(_dataset; band::Int=0, bands=Vector{Int}(), sds::Int=0, pad::Int=0, layout::String="")

	(isa(_dataset, GMTgrid) || isa(_dataset, GMTimage) || isGMTdataset(_dataset)) && return _dataset

	# For some bloody reason it would print annoying (& false?) warning messages. Have to use brute force
	#Gdal.CPLPushErrorHandler(@cfunction(Gdal.CPLQuietErrorHandler, Cvoid, (UInt32, Cint, Cstring)))
	dataset, scale_factor, add_offset, got_fill_val, fill_val = gd2gmt_helper(_dataset, sds)
	#Gdal.CPLPopErrorHandler();
	(!isa(_dataset, String) && _dataset.ptr == C_NULL) && error("NULL dataset sent in")

	n_dsbands = Gdal.nraster(dataset)
	xSize, ySize, nBands = Gdal.width(dataset), Gdal.height(dataset), n_dsbands
	dType = Gdal.pixeltype(getband(dataset, 1))
	is_grid = (sizeof(dType) >= 4 || dType == Int16) ? true : false		# Simple (too simple?) heuristic
	in_bands = (!isempty(bands)) ? bands : ((band == 0) ? collect(1:n_dsbands) : [band])
	(maximum(in_bands) > n_dsbands) && error("One selected band is larger then number of bands in this dataset")
	ncol, nrow = xSize+2pad, ySize+2pad
	(scale_factor != 1 || add_offset != 0) && (dType = Float32)
	mat = (dataset isa Gdal.AbstractRasterBand) ? zeros(dType, ncol, nrow) : zeros(dType, ncol, nrow, length(in_bands))
	n_colors = 0
	desc = Vector{String}(undef,0)
	if (isa(dataset, Gdal.AbstractRasterBand))
		Gdal.rasterio!(dataset, mat, in_bands, 0, 0, xSize, ySize, Gdal.GF_Read, 0, 0, C_NULL, pad)
		colormap, n_colors = get_cpt_from_colortable(dataset)
		append!(desc, [Gdal.GDALGetDescription(dataset.ptr)])
	else
		ds = (isa(dataset, Gdal.RasterDataset)) ? dataset.ds : dataset
		Gdal.rasterio!(ds, mat, in_bands, 0, 0, xSize, ySize, Gdal.GF_Read, 0, 0, 0, C_NULL, pad)
		colormap, n_colors = get_cpt_from_colortable(ds)
		for bd in in_bands
			append!(desc, [Gdal.GDALGetDescription(Gdal.GDALGetRasterBand(ds.ptr, bd))])
		end
	end

	(!isa(mat, Matrix) && size(mat,3) == 1) && (mat = reshape(mat, size(mat,1), size(mat,2)))	# Fck pain
	if (layout != "")		# From GDAL it always come as a TR but not sure about the interleave
		if     (startswith(layout, "BR"))  mat = reverse(mat, dims=1)		# Just flipUD
		elseif (startswith(layout, "TC"))  mat = collect(mat')
		elseif (startswith(layout, "BC"))  mat = reverse(mat', dims=1)
		else   @warn("Unsuported layout ($(layout)) change")
		end
		layout = layout[1:2] * "B"		# Till further knowledge, assume it's always Band interleaved
	end

	# If we found a scale_factor above, apply it
	(scale_factor != 1) && (mat = gd2gmt_helper_scalefac(mat, scale_factor, add_offset, got_fill_val, fill_val))

	try
		global gt = getgeotransform(dataset)
	catch
		global gt = [0.5, 1.0, 0.0, ySize+0.5, 0.0, 1.0]	# Resort to no coords
	end

	x_inc, y_inc = gt[2], abs(gt[6])
	x_min, y_max = gt[1], gt[4]
	(is_grid) && (x_min += x_inc/2;	 y_max -= y_inc/2)	# Maitain the GMT default that grids are gridline reg.
	x_max = x_min + (xSize - 1*is_grid - 2pad) * x_inc
	y_min = y_max - (ySize - 1*is_grid - 2pad) * y_inc
	z_min, z_max = (is_grid) ? extrema_nan(mat) : extrema(mat)
	hdr = [x_min, x_max, y_min, y_max, z_min, z_max, Float64(!is_grid), x_inc, y_inc]
	prj = getproj(dataset)
	(prj != "" && !startswith(prj, "+proj")) && (prj = toPROJ4(importWKT(prj)))
	(prj == "") && (prj = seek_wkt_in_gdalinfo(gdalinfo(dataset)))
	if (is_grid)
		(eltype(mat) == Float64) && (mat = Float32.(mat))
		O = mat2grid(mat; hdr=hdr, proj4=prj, names=desc)
		O.layout = (layout == "") ? "TRB" : layout
	else
		O = mat2img(mat; hdr=hdr, proj4=prj, noconv=true, names=desc)
		O.layout = (layout == "") ? "TRBa" : layout * "a"
		if (n_colors > 0)
			O.colormap = colormap;	O.n_colors = n_colors
			((nodata = Gdal.getnodatavalue(Gdal.getband(dataset))) !== nothing) && (O.nodata = nodata)
		end
	end
	O.inc = [x_inc, y_inc]		# Reset because if pad != 0 they were recomputed inside the mat2? funs
	O.pad = pad
	return O
end

# ---------------------------------------------------------------------------------------------------
function gd2gmt_helper_scalefac(mat, scale_factor, add_offset, got_fill_val, fill_val)
	# Apply a scale + offset
	(got_fill_val) && (nodata = isnodata(mat, fill_val))
	if (eltype(mat) <: Integer)
		mat = mat .* scale_factor .+ add_offset		# Also promotes (and pay) the array to float32
	else
		@inbounds Threads.@threads for k = 1:length(mat)  mat[k] = mat[k] * scale_factor + add_offset  end
	end
	(got_fill_val) && (mat[nodata] .= NaN32)
	mat
end

# ---------------------------------------------------------------------------------------------------
function gd2gmt_helper(input, sds)
	# Deal with datasets as file names. Extract SUBDATASETS if wished.
	# But `input` can also be gdal dataset, case in which the `sds` option will probably not work.
	scale_factor, add_offset, got_fill_val, fill_val = Float32(1), Float32(0), false, Float32(0)
	if (sds > 0)					# Must fish the SUBDATASET name in in gdalinfo
		# Let's fish the SDS 1 name in: SUBDATASET_1_NAME=NETCDF:"AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc":sst
		info = gdalinfo(input)
		ind = findall("SUBDATASET_", info)
		nn = 0
		for k = 1:2:length(ind)
			ind2 = findfirst('_', info[ind[1][end]+1 : ind[1][end]+3])	# Allow for up to 999 SDSs
			n = tryparse(Int, info[ind[k][end]+1 : ind[k][end]+ind2[1]-1])
			if (n == sds)			# Ok, found it, but we still need to find the EOL
				ind3 = findfirst('\n', info[ind[k][1]:end])
				input = info[ind[k][1] : ind[k][1] + ind3[1]-2]
				nn = n
				break
			end
		end
		(nn == 0) && error("SUBDATASET $(sds) not found in " * input)
	end
	if (isa(input, AbstractString))
		sds_name = trim_SUBDATASET_str(input)
		((dataset = Gdal.unsafe_read(sds_name)) == C_NULL) && error("GDAL failed to read " * sds_name)
	else
		dataset = input
	end

	# Hmmm, check also for scale_factor, add_offset, _FillValue
	info = gdalinfo(dataset)
	(info === nothing) && error("GDAL failed to find " * sds_name)
	if (occursin("Metadata:", info))
		if ((ind = findfirst("scale_factor=", info)) !== nothing)	# OK, found one
			ind2 = findfirst('\n', info[ind[1]:end])
			scale_factor = tryparse(Float32, info[ind[1]+13 : ind[1] + ind2[1]-2])
			add_offset = Float32(0)
			if ((ind = findfirst("add_offset=", info)) !== nothing)
				ind2 = findfirst('\n', info[ind[1]:end])
				add_offset = tryparse(Float32, info[ind[1]+11 : ind[1] + ind2[1]-2])
			end
			fill_val, got_fill_val = get_FillValue(info)
		end
	end
	return dataset, scale_factor, add_offset, got_fill_val, fill_val
end

# ---------------------------------------------------------------------------------------------------
function gd2gmt(geom::Gdal.AbstractGeometry, proj::String="")::Vector{<:GMTdataset}
	# Convert a geometry into a single GMTdataset
	gmtype = Gdal.getgeomtype(geom)
	if (gmtype == Gdal.wkbPolygon)		# getx() doesn't work for polygons
		geom = Gdal.getgeom(geom,0)
	elseif (gmtype == wkbMultiPolygon || gmtype == wkbMultiLineString)
		n_pts = Gdal.ngeom(geom)
		D = Vector{GMTdataset}(undef, n_pts)
		[D[k] = gd2gmt(Gdal.getgeom(geom,k-1), "")[1] for k = 1:n_pts]
		(proj != "") && (D[1].proj4 = proj)
		return D
	elseif (gmtype == wkbMultiPoint)
		n_dim, n_pts = Gdal.getcoorddim(geom), Gdal.ngeom(geom)
		mat = Array{Float64,2}(undef, n_pts, n_dim)
		[mat[k,1] = Gdal.getx(Gdal.getgeom(geom,k-1), 0) for k = 1:n_pts]
		[mat[k,2] = Gdal.gety(Gdal.getgeom(geom,k-1), 0) for k = 1:n_pts]
		(n_dim == 3) && ([mat[k,2] = Gdal.getz(Gdal.getgeom(geom,k-1), 0) for k = 1:n_pts])
		return [GMTdataset(mat, String[], "", String[], proj, "", gmtype)]
	end

	n_dim, n_pts = Gdal.getcoorddim(geom), Gdal.ngeom(geom)
	n = (n_dim == 2) ? 2 : 3
	mat = Array{Float64,2}(undef, n_pts, n)
	[mat[k,1] = Gdal.getx(geom, k-1) for k = 1:n_pts]
	[mat[k,2] = Gdal.gety(geom, k-1) for k = 1:n_pts]
	(n_dim == 3) && ([mat[k,3] = Gdal.getz(geom, k-1) for k = 1:n_pts])
	[GMTdataset(mat, String[], "", String[], proj, "", gmtype)]
end

# ---------------------------------------------------------------------------------------------------
function gd2gmt(dataset::Gdal.AbstractDataset)
	# This method is for OGR formats only
	(Gdal.OGRGetDriverByName(Gdal.shortname(getdriver(dataset))) == C_NULL) && return gd2gmt(dataset; pad=0)

	D, ds = Vector{GMTdataset}(undef, Gdal.ngeom(dataset)), 1
	for k = 1:Gdal.nlayer(dataset)
		layer = getlayer(dataset, 0)
		Gdal.resetreading!(layer)
		proj = ((p = getproj(layer)) != C_NULL) ? toPROJ4(p) : ""
		while ((feature = Gdal.nextfeature(layer)) !== nothing)
			for j = 0:Gdal.ngeom(feature)-1
				geom = Gdal.getgeom(feature, j)
				_D = gd2gmt(geom, proj)
				gt = Gdal.getgeomtype(geom)
				# Maybe when there nlayers > 1 or other cases, starting allocated size is not enough
				(length(_D) + ds >= length(D)) && append!(D, Vector{GMTdataset}(undef, round(Int, 0.5 * length(D))))
				for d in _D
					D[ds] = d
					D[ds].geom = gt
					ds += 1
				end
			end
		end
	end
	(length(D) != ds-1) && (D = deleteat!(D,ds:length(D)))
	return D
end

# ---------------------------------------------------------------------------------------------------
function get_cpt_from_colortable(dataset)
	# Extract the color info from a GDAL colortable and put it in a row vector for GMTimage.colormap
	band = (!isa(dataset, Gdal.AbstractRasterBand)) ? Gdal.getband(dataset) : dataset
	ct = Gdal.getcolortable(band)
	(ct.ptr == C_NULL) && return Vector{Clong}(), 0
	n_colors = Gdal.ncolorentry(ct)
	cmap, n = Vector{Clong}(undef, 4 * n_colors), 1
	for k = 0:n_colors-1
		c = Gdal.getcolorentry(ct, k)
		cmap[n] = c.c1;	n += 1; cmap[n] = c.c2;	n += 1; cmap[n] = c.c3;	n += 1; cmap[n] = c.c4;	n += 1;
	end
	return cmap, n_colors
end

# ---------------------------------------------------------------------------------------------------
function trim_SUBDATASET_str(sds::String)
	# If present, trim the initial "SUBDATASET_X_NAME=" from a subdataset string name
	ind = 0
	(startswith(sds, "SUBDATASET_") && (ind = findfirst('=', sds)) === nothing) && error("Badly formed SUBDATASET string")
	sds_name = sds[ind+1:end]
end

# ---------------------------------------------------------------------------------------------------
function get_FillValue(str::String)
	# Search for a _FillValue in str. Return it if found or NaN otherwise. We use the second return
	# value to tell between a true NaN as _FillValue and one that's only use for type stability
	((ind = findfirst("_FillValue=", str)) === nothing) && return NaN, false
	ind2 = findfirst('\n', str[ind[1]:end])
	fill_val = (ind2 !== nothing) ? tryparse(Float32, str[ind[1]+11 : ind[1] + ind2[1]-2]) : tryparse(Float32, str[ind[1]+11:end])
	return fill_val, true
end

# ---------------------------------------------------------------------------------------------------
function seek_wkt_in_gdalinfo(info::String)
	# Seek for a SRS string in gdalinfo info. Seems that files can have no explicit projinfo but still
	# expose them with gdalinfo.
	((ind = findfirst("SRS=GEO", info)) === nothing) && return ""
	ind2 = findfirst('\n', view(info, ind[5]:length(info)))		# Find next EOL
	proj4 = toPROJ4(importWKT(info[ind[5] : ind[5]+ind2[1]-2]))
end

# ---------------------------------------------------------------------------------------------------
"""
    ds = gmt2gd(GI)

Create GDAL dataset from the contents of GI that can be either a Grid or an Image

    ds = gmt2gd(D, save="", geometry="")

Create GDAL dataset from the contents of D, which can be a GMTdataset, a vector of GMTdataset ir a MxN array.
The `save` keyword instructs GDAL to save the contents as an OGR file. Format is determined by file estension.
`geometry` can be a string with "polygon", where file will be converted to polygon/multipolygon depending
on `D` is a single or a multi-segment object, or "point" to convert to a multipoint geometry.
"""
function gmt2gd(GI)
	width, height = (GI.layout != "" && GI.layout[2] == 'C') ? (size(GI,2), size(GI,1)) : (size(GI,1), size(GI,2))
	ds = Gdal.create("", driver=getdriver("MEM"), width=width, height=height, nbands=size(GI,3), dtype=eltype(GI[1]))
	if (isa(GI, GMTgrid))
		if (GI.layout != "" && GI.layout[2] == 'C')
			indata = (GI.layout[1] == 'B') ? collect(reverse(GI.z, dims=1)') : collect(GI.z')
		else
			indata = GI.z
		end
		Gdal.write!(ds, indata, isa(GI.z, Array{<:Real, 3}) ? Cint.(collect(1:size(GI,3))) : 1)
	elseif (isa(GI, GMTimage))
		if (GI.layout != "" && GI.layout[2] == 'C')
			indata = (GI.layout[1] == 'B') ? collect(reverse(GI.image, dims=1)') : collect(GI.image')
		else
			indata = GI.image
		end
		Gdal.write!(ds, indata, isa(GI.image, Array{<:Real, 3}) ? Cint.(collect(1:size(GI,3))) : 1)

		if (GI.n_colors > 0)
			ct = Gdal.createcolortable(UInt32(1))	# RGB
			n = 1
			for k = 0:GI.n_colors-1
				c1, c2, c3, c4 = GI.colormap[n],GI.colormap[n+1],GI.colormap[n+2],GI.colormap[n+3]
				#Gdal.createcolorramp!(ct, k, Gdal.GDALColorEntry(c1,c2,c3,c4), k, Gdal.GDALColorEntry(c1,c2,c3,c4))
				Gdal.setcolorentry!(ct, k, Gdal.GDALColorEntry(c1,c2,c3,c4));
				n += 4
			end
			Gdal.setcolortable!(Gdal.getband(ds), ct)
			(!isnan(GI.nodata)) && (setnodatavalue!(Gdal.getband(ds), GI.nodata))
		end
	end

	if (!isempty(GI.names))			# Set bands description
		[Gdal.GDALSetDescription(Gdal.GDALGetRasterBand(ds.ptr, k), GI.names[k]) for k = 1:min(size(GI,3), length(GI.names))]
	end

	x_min, y_max = GI.range[1], GI.range[4]
	(GI.registration == 0) && (x_min -= GI.inc[1]/2;  y_max += GI.inc[2]/2)
	setgeotransform!(ds, [x_min, GI.inc[1], 0.0, y_max, 0.0, -GI.inc[2]])
	if     (GI.wkt != "")    setproj!(ds, GI.wkt)
	elseif (GI.proj4 != "")  setproj!(ds, toWKT(importPROJ4(GI.proj4), true))
	end
	return ds
end

# ---------------------------------------------------------------------------------------------------
gmt2gd(D::Array{<:Real,2}; save::String="", geometry::String="") = gmt2gd(mat2ds(D); save=save, geometry=geometry)
gmt2gd(D::GMTdataset; save::String="", geometry::String="") = gmt2gd([D]; save=save, geometry=geometry)
function gmt2gd(D::Vector{<:GMTdataset}; save::String="", geometry::String="")
	# ...
	n_cols = size(D[1].data, 2)
	(n_cols < 2) && error("GMTdataset must have at least 2 columns")

	geometry = lowercase(geometry)
	ismulti  = (length(D) > 1) || (D[1].geom == wkbPoint && size(D[1].data, 1) > 1)
	ispolyg  = occursin("poly", geometry);
	isline   = occursin("line", geometry);
	ispoint  = occursin("point", geometry);
	(geometry != "" && !isline && !ispoint && !ispolyg) && error("Geometry $(geometry) not yet implemented")
	if (D[1].geom == 0 && !isline && !ispoint && !ispolyg)	# If all multi-segments are closed create a Polygon/MultiPolygon
		ispolyg = true
		for k = 1:length(D)
			(D[k].data[1,1:2] != D[k].data[end,1:2]) && (ispolyg = false; break)
		end
		isline = !ispolyg						# Otherwise make a Line/MultiLine
	end

	ds = Gdal.create(getdriver("MEMORY"));
	#ds = Gdal.create(getdriver("ESRI Shapefile"), filename="/vsimem/mem.shp")
	if     (D[1].proj4 != "")  sr = Gdal.importPROJ4(D[1].proj4)
	elseif (D[1].wkt   != "")  sr = Gdal.importWKT(D[1].wkt)
	else                       sr = Gdal.ISpatialRef(C_NULL)
	end

	if (ispolyg || D[1].geom == wkbPolygon || D[1].geom == wkbMultiPolygon)	# If guessed or in Dataset
		geom_code, geom_cmd = (!ismulti) ? (wkbPolygon, Gdal.createpolygon()) :
		                                   (wkbMultiPolygon, Gdal.createmultipolygon())
	elseif (isline || D[1].geom == wkbLineString || D[1].geom == wkbMultiLineString)
		geom_code, geom_cmd = (!ismulti) ? (wkbLineString, Gdal.createlinestring()) :
		                                   (wkbMultiLineString, Gdal.createmultilinestring())
	elseif (D[1].geom == wkbMultiPoint || (ispoint && !ismulti))
		geom_code, geom_cmd = wkbMultiPoint, Gdal.createmultipoint()
	else
		geom_code, geom_cmd = (!ismulti) ? (wkbPoint, Gdal.createpoint()) :
		                                   (wkbMultiPoint, Gdal.createmultipoint())
	end

	layer = Gdal.createlayer(name="layer1", dataset=ds, geom=geom_code, spatialref=sr);
	feature = Gdal.unsafe_createfeature(layer)
	geom = geom_cmd

	if (ispolyg || D[1].geom == wkbPolygon || D[1].geom == wkbMultiPolygon)
		if (ismulti)
			for k = 1:length(D)
				poly = Gdal.creategeom(Gdal.wkbPolygon)
				Gdal.addgeom!(geom, Gdal.addgeom!(poly, makering(D[k].data)))
			end
		else			# Polygons with islands
			[Gdal.addgeom!(geom, makering(D[k].data)) for k = 1:length(D)]
		end
		Gdal.setgeom!(feature, geom)
	elseif (isline || D[1].geom == wkbLineString || D[1].geom == wkbMultiLineString)
		if (ismulti)
			for k = 1:length(D)
				line = Gdal.creategeom(wkbLineString)
				x,y,z = helper_gmt2gd_xyz(D[k], n_cols)
				(n_cols == 2) ? Gdal.OGR_G_SetPoints(line.ptr, size(D[k].data, 1), x, 8, y, 8, C_NULL, 8) :
				                Gdal.OGR_G_SetPoints(line.ptr, size(D[k].data, 1), x, 8, y, 8, z, 8)
        		Gdal.addgeom!(geom, line)
			end
		else
			x,y,z = helper_gmt2gd_xyz(D[1], n_cols)
			(n_cols == 2) ? Gdal.OGR_G_SetPoints(geom.ptr, size(D[1].data, 1), x, 8, y, 8, C_NULL, 8) :
			                Gdal.OGR_G_SetPoints(geom.ptr, size(D[1].data, 1), x, 8, y, 8, z, 8)
		end
		Gdal.setgeom!(feature, geom)
	elseif (ispoint || D[1].geom == wkbPoint || D[1].geom == wkbMultiPoint)
		if (D[1].geom == wkbMultiPoint || ismulti)
			for k = 1:length(D)
				x,y,z = helper_gmt2gd_xyz(D[k], n_cols)
				if (n_cols == 2)
					[Gdal.addgeom!(geom, Gdal.createpoint(x[n], y[n])) for n = 1:length(x)]
				else
					[Gdal.addgeom!(geom, Gdal.createpoint(x[n], y[n], z[n])) for n = 1:length(x)]
				end
			end
		else
			x,y,z = helper_gmt2gd_xyz(D[1], n_cols)
			(n_cols == 2) ? Gdal.OGR_G_SetPoints(geom.ptr, size(D[1].data, 1), x, 8, y, 8, C_NULL, 8) :
			                Gdal.OGR_G_SetPoints(geom.ptr, size(D[1].data, 1), x, 8, y, 8, z, 8)
		end
		Gdal.setgeom!(feature, geom)
	else
		@warn("Geometries with geometry code $(D[1].geom) are not yet implemented")
	end

	Gdal.setfeature!(layer, feature)
	Gdal.destroy(feature)

	if (save != "")
		ogr2ogr(ds, dest=save)
		Gdal.destroy(ds)
		return nothing
	end
	return ds
end

function helper_gmt2gd_xyz(D::GMTdataset, n_cols::Int)
	# Helper funtion to split a matrix in x,y[z] and make sure doubles are returned
	if (eltype(D) == Float64)
		x, y = D.data[:,1], D.data[:,2]
		(n_cols > 2) && (z = D.data[:,3])
	else
		x, y = Float64.(D.data[:,1]), Float64.(D.data[:,2])
		(n_cols > 2) && (z = Float64.(D.data[:,3]))
	end
	return (n_cols == 3) ? (x, y, z) : (x, y, 0.0)
end

	# This is possibly a wasting solution implying a data copy but it's bloody simpler
	#if (conv2polyg)
		#geom.ptr = (length(D) == 1) ? Gdal.OGR_G_ForceToPolygon(geom.ptr) : Gdal.OGR_G_ForceToMultiPolygon(geom.ptr)
	#end

function makering(data)
	ring = Gdal.creategeom(Gdal.wkbLinearRing)
	if (size(data,2) == 2)
		[Gdal.addpoint!(ring, data[k,1], data[k,2]) for k = 1:size(data,1)]
	else
		[Gdal.addpoint!(ring, data[k,1], data[k,2], data[k,3]) for k = 1:size(data,1)]
	end
	ring
end

# ---------------------------------------------------------------------------------------------------
function R_inc_to_gd(inc::Vector{Float64}, opt_R::String="", BB::Vector{Float64}=Vector{Float64}())
	# Convert an opt_R string or BB vector + inc vector in the equivalent set of options for GDAL
	if (opt_R != "")
		ind = findfirst("R", opt_R)
		BB = (ind !== nothing) ? tryparse.(Float64, split(opt_R[ind[1]+1:end], '/')) : tryparse.(Float64, split(opt_R, '/'))
	end
	nx = round(Int32, (BB[2] - BB[1]) / inc[1])
	inc_y = (length(inc) == 1) ? inc[1] : inc[2]
	ny = round(Int32, (BB[4] - BB[3]) / inc_y)
	return ["-txe", "$(BB[1])", "$(BB[2])", "-tye", "$(BB[3])", "$(BB[4])", "-outsize", "$(nx)", "$(ny)"]
end

# ---------------------------------------------------------------------------------------------------
"""
	gdalshade(filename; kwargs...)

Create a shaded relief with the GDAL method (color image blended with shaded intensity).

- kwargs hold the keyword=value to pass the arguments to `gdaldem hillshade`

Example:
    I = gdalshade("hawaii_south.grd", C="faa.cpt", zfactor=4);

### Returns
A GMT RGB Image
"""
function gdalshade(fname; kwargs...)
	d = KW(kwargs)
	band = ((val = find_in_dict(d, [:band], false)[1]) !== nothing) ? string(val) : "1"
	cmap = find_in_dict(d, CPTaliases)[1]	# The color cannot be passed to second call to gdaldem

	A = gdaldem(fname, "color-relief", ["-b", band], C=cmap)
	B = gdaldem(fname, "hillshade"; d...)
	blendimg!(A, B)
end

# ---------------------------------------------------------------------------------------------------
"""
    gdalread(fname::AbstractString, opts=String[]; gdataset=false, kwargs...)

Read a raster or a vector file from a disk file and return the result either as a GMT type (the default)
or a GDAL dataset.

- `fname`: Input data. It can be a file name, a GMTgrid or GMTimage object or a GDAL dataset
- `opts`:  List of options. The accepted options are the ones of the gdal_translate utility.
           This list can be in the form of a vector of strings, or joined in a simgle string.
- `gdataset`: If set to `true` forces the return of a GDAL dataset instead of a GMT type.
- `kwargs`: This options accept the GMT region (-R) and increment (-I)

### Returns
A GMT grid/image or a GDAL dataset
"""
function gdalread(fname::AbstractString, optsP=String[]; opts=String[], gdataset=false, kw...)
	(fname == "") && error("Input file name is missing.")
	(isempty(optsP) && !isempty(opts)) && (optsP = opts)		# Accept either Positional or KW argument
	ressurectGDAL();
	ds_t = Gdal.read(fname, I=false)
	if (Gdal.OGRGetDriverByName(Gdal.shortname(getdriver(ds_t))) == C_NULL)
		ds = gdaltranslate(ds_t, optsP; gdataset=gdataset, kw...)
	else
		optsP = (isempty(optsP)) ? ["-overwrite"] : append!(optsP, "-overwrite")
		ds = ogr2ogr(ds_t, optsP; gdataset=true, kw...)
		Gdal.deletedatasource(ds, "/vsimem/tmp")		# WTF I need to do this?
	end
	Gdal.GDALClose(ds_t.ptr)			# WTF it needs explicit close?
	return (gdataset) ? ds : gd2gmt(ds)
end

# ---------------------------------------------------------------------------------------------------
const ncType = Dict(
	Int8    => 1,
	UInt8   => 7,
	Int16   => 3,
	UInt16  => 8,
	Int32   => 4,
	UInt32  => 9,
	Int64   => 10,
	UInt64  => 11,
	Float32 => 5,
	Float64 => 6,
	Char    => 2,
	String  => 12)
"""
    gdalwrite(data, fname::AbstractString, opts=String[]; kwargs...)
or

    gdalwrite(fname::AbstractString, data, opts=String[]; kwargs...)

Write a raster or a vector file to disk

- `fname`: Output file name. If not explicitly selected via `opts` the used driver will be picked from the file extension.
- `data`:  The data to be saved in file. It can be a GMT type or a GDAL dataset.
- `opts`:  List of options. The accepted options are the ones of the gdal_translate or ogr2ogr utility.
           This list can be in the form of a vector of strings, or joined in a simgle string.
- `kwargs`: This options accept the GMT region (-R) and increment (-I)

or

    gdalwrite(cube, fname::AbstractString, v=nothing; dim_name::String="time", dim_units::String="")

Write a MxNxP `cube` object to disk as a multilayered file.

- `cube`: A GMTgrid or GMTimage cube
- `fname`: The file name where to save the `cube`
- `v`: A vector ith the coordinates of the Z layers (if omitted create one as 1:size(cube,3))
- `dim_name`: The name of the variable of the ``vertical`` dimension.
- `dim_units`: The units of the `v` vector. If not provided, use the `cube.z_units` if exist (GMTgrid only)
"""
gdalwrite(fname::AbstractString, data, optsP=String[]; opts=String[], kw...) = gdalwrite(data, fname, optsP, true; opts=opts, kw...)
function gdalwrite(data, fname::AbstractString, optsP=String[], pickme::Bool=true; opts=String[], kw...)
	(fname == "") && error("Output file name is missing.")
	(isempty(optsP) && !isempty(opts)) && (optsP = opts)		# Accept either Positional or KW argument
	ds, = Gdal.get_gdaldataset(data, optsP)
	if (Gdal.OGRGetDriverByName(Gdal.shortname(getdriver(ds))) == C_NULL)
		gdaltranslate(ds, optsP; dest=fname, kw...)
	else
		ogr2ogr(ds, optsP; dest=fname, kw...)
	end
end

function gdalwrite(cube::GItype, fname::AbstractString, v=nothing; dim_name::String="time", dim_units::String="")
	nbands = size(cube,3)
	(v === nothing && !isempty(cube.v)) && (v = cube.v)
	_v = (v === nothing) ? "$(collect(1:nbands))" : "$(v)"
	ds = gmt2gd(cube)
	Gdal.setmetadataitem(ds, "NETCDF_DIM_EXTRA", dim_name)
	code = (v === nothing) ? 4 : ncType[eltype(v)]		# Pick the v data type.
	Gdal.setmetadataitem(ds, "NETCDF_DIM_" * dim_name * "_DEF", "{$(nbands),$(code)}")
	Gdal.setmetadataitem(ds, "NETCDF_DIM_" * dim_name * "_VALUES", "{" *_v[2:end-1] * "}")
	Gdal.setmetadataitem(ds, dim_name * "#axis", string(dim_name[1]))
	#Gdal.setmetadataitem(ds, "Band1#actual_range","{0,50000}")
	if (dim_units != "")                              Gdal.setmetadataitem(ds, dim_name * "#units", dim_units)
	elseif (isa(cube, GMTgrid) && cube.z_unit != "")  Gdal.setmetadataitem(ds, dim_name * "#units", cube.z_unit)
	end
	crs = Gdal.getproj(cube, wkt=true)
	(crs != "" ) && Gdal.setproj!(ds, crs)
	Gdal.unsafe_copy(ds, filename=fname, driver=getdriver("netCDF"), options=["FORMAT=NC4", "COMPRESS=DEFLATE", "ZLEVEL=4"])
	return nothing
end

# TODO ==> Input as vector. EPSGs. Test input t_srs
# ---------------------------------------------------------------------------------------------------
"""
    lonlat2xy(lonlat::Matrix{<:Real}, t_srs::String; s_srs=::String="+proj=longlat +datum=WGS84")
or

    lonlat2xy(D::GMTdataset, t_srs::String; s_srs=::String="+proj=longlat +datum=WGS84")

Computes the forward projection from LatLon to XY in the given projection. The input is assumed to be in WGS84.
If it isn't, pass the appropriate projection info via the `s_srs` option.

### Parameters
* `lonlat`: The input data. It can be a Matrix, or a GMTdataset (or vector of it)
* `t_srs`:  The destiny projection system. This can be a PROJ4 or a WKT string

### Returns
A Matrix if input is a Matrix or a GMTdadaset if input had that type
"""
lonlat2xy(xy::Vector{<:Real}, t_srs::String; s_srs::String=prj4WGS84) = vec(lonlat2xy(reshape(xy[:],1,length(xy)), t_srs; s_srs=s_srs))
function lonlat2xy(lonlat::Matrix{<:Real}, t_srs::String; s_srs::String=prj4WGS84)
	D = ogr2ogr(lonlat, ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite"])
	return D[1].data		# Return only the array because that's what was sent in
end

lonlat2xy(D::GMTdataset, t_srs::String; s_srs::String=prj4WGS84) = lonlat2xy([D], t_srs; s_srs=s_srs)
function lonlat2xy(D::Vector{<:GMTdataset}, t_srs::String; s_srs::String=prj4WGS84)
	(startswith(D[1].proj4, "+proj=longl") || startswith(D[1].proj4, "+proj=latlon")) && (s_srs = D[1].proj4)
	o = ogr2ogr(D, ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite"])
	(isa(o, Gdal.AbstractDataset)) && (o = gd2gmt(o))
	o
end

# ---------------------------------------------------------------------------------------------------
"""
    xy2lonlat(xy::Matrix{<:Real}, s_srs::String; t_srs=::String="+proj=longlat +datum=WGS84")
or

    xy2lonlat(D::GMTdataset, s_srs::String; t_srs=::String="+proj=longlat +datum=WGS84")

Computes the inverse projection from XY to LonLat in the given projection. The output is assumed to be in WGS84.
If that isn't right, pass the appropriate projection info via the `t_srs` option.

### Parameters
* `xy`: The input data. It can be a Matrix, or a GMTdataset (or vector of it)
* `s_srs`:  The data projection system. This can be a PROJ4 or a WKT string
* `t_srs`:  The target SRS. If the default is not satisfactory, provide a new projection info (PROJ4 or WKT)

### Returns
A Matrix if input is a Matrix or a GMTdadaset if input had that type
"""
xy2lonlat(xy::Vector{<:Real}, s_srs::String; t_srs::String=prj4WGS84) = vec(xy2lonlat(reshape(xy[:],1,length(xy)), s_srs; t_srs=t_srs))
function xy2lonlat(xy::Matrix{<:Real}, s_srs::String; t_srs::String=prj4WGS84)
	D = ogr2ogr(xy, ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite"])
	return D[1].data		# Return only the array because that's what was sent in
end

xy2lonlat(D::GMTdataset, s_srs::String=""; t_srs::String=prj4WGS84) = xy2lonlat([D], s_srs; t_srs=t_srs)
function xy2lonlat(D::Vector{<:GMTdataset}, s_srs::String=""; t_srs::String=prj4WGS84)
	(D[1].proj4 == "" && D[1].wkt == "" && s_srs == "") && error("No projection information whatsoever on the input data.")
	if (s_srs != "") _s_srs = s_srs
	else             _s_srs = (D[1].proj4 != "") ? D[1].proj4 : D[1].wkt
	end
	o = ogr2ogr(D, ["-s_srs", _s_srs, "-t_srs", t_srs, "-overwrite"])
	(isa(o, Gdal.AbstractDataset)) && (o = gd2gmt(o))
	o
end


text:
"""
    jointoperator_bc(operators, Q::Array)

Returns a discretized operator that solves systems of differential equations defined by
`operators` with transitions by `Q` where `operators` is an N-length collection of
discretized operators with boundary conditions applied and `Q` is N by N intensity matrix.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> L1bc = L₁₋bc(x̄, (Reflecting(), Reflecting()))
4×4 LinearAlgebra.Tridiagonal{Float64,Array{Float64,1}}:
  0.0   0.0    ⋅    ⋅
 -1.0   1.0   0.0   ⋅
   ⋅   -1.0   1.0  0.0
   ⋅     ⋅   -1.0  1.0

julia> L2bc = L₂bc(x̄, (Reflecting(), Reflecting()))
4×4 LinearAlgebra.Tridiagonal{Float64,Array{Float64,1}}:
 -1.0   1.0    ⋅     ⋅
  1.0  -2.0   1.0    ⋅
   ⋅    1.0  -2.0   1.0
   ⋅     ⋅    1.0  -1.0

julia> Q = [-0.5 0.5; 0.3 -0.3]
2×2 Array{Float64,2}:
 -0.5   0.5
  0.3  -0.3

julia> jointoperator_bc((L1bc, L2bc), Q)
2×2-blocked 8×8 BlockBandedMatrices.BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}}}:
 -0.5   0.0    ⋅    ⋅   │   0.5   0.0    ⋅     ⋅
 -1.0   0.5   0.0   ⋅   │   0.0   0.5   0.0    ⋅
   ⋅   -1.0   0.5  0.0  │    ⋅    0.0   0.5   0.0
   ⋅     ⋅   -1.0  0.5  │    ⋅     ⋅    0.0   0.5
 ───────────────────────┼────────────────────────
  0.3   0.0    ⋅    ⋅   │  -1.3   1.0    ⋅     ⋅
  0.0   0.3   0.0   ⋅   │   1.0  -2.3   1.0    ⋅
   ⋅    0.0   0.3  0.0  │    ⋅    1.0  -2.3   1.0
   ⋅     ⋅    0.0  0.3  │    ⋅     ⋅    1.0  -1.3
```
"""
function jointoperator_bc(operators, Q::Array)
    M = size(operators[1],1)
    N = length(operators)

    # check if all operators are square
    @assert all(operator->(size(operator,1) == size(operator,2)),
                operators)

    # check if all operators have same size
    @assert all(operator->(size(operator) == size(operators[1])),
                operators)

    # check if the size of transition matrix is
    # same as the number of operators
    @assert size(Q,1) == size(Q,2) == N

    # extract operators and append them to form a diagonal block tridiagonal
    Ls = blockdiag(sparse.(operators)...)
    Ls = BandedBlockBandedMatrix(Ls, (M*ones(Int64, N), M*ones(Int64, N)), (0,0), (1,1))

    # construct a kronecker product of Q times I_M
    Qs = BandedBlockBandedMatrix(Kron(Q, Eye(M)))

    return (Ls + Qs)
end


text:
function burgularSimple()

    burgular,earthquake,alarm,radio = 1,2,3,4 #variables

    yes,no=1,2 #states

    #useful for printing reaults
    dictVariable=Dict{Integer,DiscreteVariable}()
    dictVariable[burglar]=DiscreteVariable("burglar",["yes","no"])
    dictVariable[earthquake]=DiscreteVariable("earthquake",["yes","no"])
    dictVariable[alarm]=DiscreteVariable("alarm",["yes","no"])
    dictVariable[radio]=DiscreteVariable("radio",["yes","no"])

    #dicrete probability tables:
    potBurglar=PotArray(burgular,[0.01 0.99])
    potEarthquake=PotArray(earthquake,[0.000001 , 1-0.000001])

    table = zeros(2,2)
    table[yes,yes]=1
    table[no,no]=0
    table[yes,no]=0
    table[no,no]=1
    potRadio=PotArray([radio earthquake],table)

    table=zeros(2,2,2)
    table[yes,yes,yes]=0.9999
    table[yes,yes,no]=0.99
    table[yes,no,yes]=0.99
    table[yes,no,no]=0.0001

    table[no,:,:] = 1-table[yes,:,:]
    potAlarm=PotArray([radio burgular earthquake], table)

    jointPot=prod([potBurglar potRadio potAlarm potEarthquake])

    println("p(burglar|alarm=yes):")
    show(condpot(setpot(jointpot,alarm,yes),burglar),dictVariable)
end


text:
##
using DelimitedFiles
using Plots
include("minimalESN.jl")

##
#esn = ESN(inSize=1, outSize=1, resSize=1000, leaking_rate=0.5, sparsity=0.6, spectral_radius=0.9)
esn = ESN() # ESN com valores padrões

##
data = readdlm("MackeyGlass_t17.txt")

##
fit!(esn, data[1:2000])

##
preds = predict!(esn, 1000, 0)
preds = preds'

##
plot(data[2001:3000])
plot!(preds)

##
rmse = sqrt( sum((data[2001:3000] - preds).^2) ./ length(preds))
println(rmse)


text:
# Run tests for the dep_distributed example

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

dep_distributed_exact_eigvals = [
   -0.400236388049641 + 0.970633098237807im
   -0.400236388049641 - 0.970633098237807im
   2.726146249832675 + 0.000000000000000im
   -1.955643591177653 + 3.364550574688863im
   -1.955643591177653 - 3.364550574688863im
   4.493937056300693 + 0.000000000000000im
   -1.631513006819252 + 4.555484848248613im
   -1.631513006819252 - 4.555484848248613im
   -1.677320660400946 + 7.496870451838560im
   -1.677320660400946 - 7.496870451838560im]

myerrmeasure=(λ,v) -> minimum(abs.(λ .- dep_distributed_exact_eigvals) ./ abs(λ))

@testset "gallery(dep_distributed" begin
    dep=nep_gallery("dep_distributed");

    n=size(dep,1);

    M1=compute_Mder(dep, dep_distributed_exact_eigvals[1]);
    @test minimum(svdvals(M1))<eps()*100;

    x=Array{Float64,1}(1:n)
    λ=1+1.0im;
    #
    @test norm(compute_Mder(dep,λ)*x-compute_Mlincomb(dep,λ,x))<eps()*10

    # Check derivative computation with finite difference
    zder=compute_Mlincomb(dep,λ,x,[1],1)
    ee=1e-4;
    zder_approx=(compute_Mlincomb(dep,λ+ee,x)-compute_Mlincomb(dep,λ-ee,x))/(2ee)
    @test norm(zder-zder_approx)<ee^2*10

    (λ,V)=iar(dep,σ=3,neigs=5,errmeasure=myerrmeasure, v=ones(n),
              logger=displaylevel,maxit=100,tol=eps()*100)

    @info λ

    @testset "IAR eigval[$i]" for i in 1:length(λ)
        @test estimate_error(myerrmeasure,λ[i],V[:,i])<1e-10
    end

    @bench @testset "Quasinewton eigval[$i]" for i in 1:length(dep_distributed_exact_eigvals[1:3])
        λ0=round(dep_distributed_exact_eigvals[i]*10)/10
        λ,v=quasinewton(ComplexF64,dep,v=ones(n),λ=λ0,
                        #displaylevel=displaylevel,
                        armijo_factor=0.5,maxit=200,
                        errmeasure=myerrmeasure,
                        tol=eps()*100)
        @test abs((dep_distributed_exact_eigvals[i]-λ)/λ)<eps()*100
    end

    @testset "Errors thrown" begin
        @test_throws MethodError nep_gallery("dep_distributed", 15)
        @test_throws MethodError nep_gallery("dep_distributed", t=15)
    end
end


text:
typealias TuningIntervals Vector{Real}

immutable Tuning
  intervals :: TuningIntervals
  base :: Real
  name :: Nullable{Symbol}

  function Tuning(n::Symbol, base::Real)
    if n == :edo_12
      new(map(n -> 2 ^ (n / 12), 1:12), base, n)
    end
  end
end

Tuning(n) = Tuning(n, 440)

@delegate Tuning.intervals [
  Base.get, Base.getindex, Base.haskey,
  Base.getkey, Base.start, Base.done, Base.next,
  Base.isempty, Base.length, Base.hash
]


text:
# ---
# title: 105. Construct Binary Tree from Preorder and Inorder Traversal
# id: problem105
# author: zhwang
# date: 2022-03-02
# difficulty: Medium
# categories: Array, Tree, Depth-first Search
# link: <https://leetcode.com/problems/construct-binary-tree-from-preorder-and-inorder-traversal/description/>
# hidden: true
# ---
#
# Given preorder and inorder traversal of a tree, construct the binary tree.
#
# **Note:**
# You may assume that duplicates do not exist in the tree.
#
# For example, given
#
#
#
#     preorder = [3,9,20,15,7]
#     inorder = [9,3,15,20,7]
#
# Return the following binary tree:
#
#
#
#         3
#        / \
#       9  20
#         /  \
#        15   7
#
#
## @lc code=start
using LeetCode

## using @view macro
function build_tree_105(preorder::AbstractArray, inorder::AbstractArray)::TreeNode
    root = TreeNode(first(preorder))
    pos = findfirst(==(root.val), inorder)
    pos != 1 && (root.left = build_tree_105(@view(preorder[2:pos]), @view(inorder[1:(pos - 1)])))
    pos != length(preorder) &&
     (root.right = build_tree_105(@view(preorder[(pos + 1):end]), @view(inorder[(pos + 1):end])))
    root
end
## @lc code=end

## @lc test=start
@testset "105.Construct-Binary-Tree-from-Preorder-and-Inorder-Traversal.jl" begin
    @test build_tree_105([3,9,20,15,7], [9,3,15,20,7]) == TreeNode{Int}([3,9,20,nothing,nothing,15,7])
    @test build_tree_105([-1], [-1]) == TreeNode(-1)
    @test build_tree_105([3, 9, 8, 5, 4, 10, 20, 15, 7], [4, 5, 8, 10, 9, 3, 15, 20, 7]) == TreeNode{Int}([3,9,20,8,nothing,15,7,5,10,nothing,nothing,nothing,nothing,4])
end
## @lc test=end

## @lc info=start
# link: [solution to 105](https://leetcode-cn.com/problems/construct-binary-tree-from-preorder-and-inorder-traversal/solution/105-cong-qian-xu-yu-zhong-xu-bian-li-xu-7lzvc/)
## @lc info=end


text:
function ju2py_obj(ju_obj, py_obj)
    py_attr = keys(py_obj)
    n_attr = length(py_attr)
    for i = 1:n_attr
        # skip base class attributes
        if !contains(String(py_attr[i]), "__")
            if py_obj[py_attr[i]] != nothing
                error("Cannot assign to non-None attribute $(py_obj[py_attr[i]]) in Python object.")
            end
            py_attr_name = py_attr[i]
            ju_key_name = String(py_attr_name)
            if ju_obj[ju_key_name] != nothing
                ju_eltype = eltype(ju_obj[ju_key_name])
                # assign Julia array of arrays to Python list of numpy arrays
                if ju_eltype == Array{Int64, 1} || ju_eltype == Array{Float64,1} || ju_eltype == Array{Int64,2} || ju_eltype == Array{Float64, 2}
                    try py_obj[py_attr_name] = PyVector(ju_obj[ju_key_name])
                    catch
                        error("Could not assign attribute $(py_attr[i]).")
                    end
                # assign Julia array to numpy array
                elseif (ju_eltype == Int64 || ju_eltype == Float64) && size(ju_obj[ju_key_name]) != ()
                    try py_obj[py_attr_name] = np.array(ju_obj[ju_key_name])
                    catch
                        error("Could not assign attribute $(py_attr[i]).")
                    end
                # assign scalar
                elseif size(ju_obj[ju_key_name]) == ()
                    try py_obj[py_attr_name] = ju_obj[ju_key_name]
                    catch
                        error("Could not assign attribute $(py_attr[i]).")
                    end
                # the Julia Dictionary containts something that it should not contain
                else
                    error("Julia dictionary contains something not an Array{Int/Float64} or Int/Float64.")
                end
            end
        end
    end
    return py_obj
end


text:
#= Custom calibration configuration file. =#

using Distributions
using StatsBase
using LinearAlgebra
using CalibrateEDMF
using CalibrateEDMF.ModelTypes
using CalibrateEDMF.DistributionUtils
using CalibrateEDMF.ReferenceModels
using CalibrateEDMF.ReferenceStats
using CalibrateEDMF.LESUtils
using CalibrateEDMF.TurbulenceConvectionUtils
const src_dir = dirname(pathof(CalibrateEDMF))
using CalibrateEDMF.HelperFuncs
# Import EKP modules
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions
using JLD2

namelist_args = [
    ("time_stepping", "t_max", 2.0 * 3600),
    ("time_stepping", "dt_max", 30.0),
    ("time_stepping", "dt_min", 20.0),
    ("stats_io", "frequency", 120.0),
    ("grid", "dz", 150.0),
    ("grid", "nz", 20),
]

# Cases defined as structs for quick access to default configs
struct Bomex end

function get_config()
    config = Dict()
    # Flags for saving output data
    config["output"] = get_output_config()
    # Define regularization of inverse problem
    config["regularization"] = get_regularization_config()
    # Define reference used in the inverse problem
    config["reference"] = get_reference_config(Bomex())
    # Define the parameter priors
    config["prior"] = get_prior_config()
    # Define the kalman process
    config["process"] = get_process_config()
    # Define the SCM static configuration
    config["scm"] = get_scm_config()
    return config
end

function get_output_config()
    config = Dict()
    config["outdir_root"] = mktempdir()
    config["overwrite_scm_file"] = false # Flag for overwritting SCM input file
    return config
end

function get_regularization_config()
    config = Dict()
    config["perform_PCA"] = true # Performs PCA on data
    config["variance_loss"] = 1.0e-2 # Variance truncation level in PCA
    config["normalize"] = true  # whether to normalize data by pooled variance
    config["tikhonov_mode"] = "relative" # Tikhonov regularization
    config["tikhonov_noise"] = 1.0e-4 # Tikhonov regularization
    config["dim_scaling"] = false # Dimensional scaling of the loss
    return config
end

function get_process_config()
    config = Dict()
    config["N_iter"] = 2
    config["N_ens"] = 5
    config["algorithm"] = "Inversion" # "Sampler", "Unscented"
    return config
end

function get_reference_config(::Bomex)
    config = Dict()
    config["case_name"] = ["Bomex"]
    # Flag to indicate whether reference data is from a perfect model (i.e. SCM instead of LES)
    config["y_reference_type"] = SCM()
    config["Σ_reference_type"] = SCM()
    config["y_names"] = [["thetal_mean", "qt_mean"]]
    ref_root_dir = mktempdir()
    config["y_dir"] = [joinpath(ref_root_dir, "Output.Bomex.ref")]
    config["scm_suffix"] = ["scm"]
    config["scm_parent_dir"] = [ref_root_dir]
    config["t_start"] = [0.0]
    config["t_end"] = [2.0 * 3600]
    config["namelist_args"] = repeat([namelist_args], 2)
    return config
end

function get_prior_config()
    config = Dict()
    config["constraints"] =
        Dict("entrainment_factor" => [bounded(0.0, 0.5)], "detrainment_factor" => [bounded(0.0, 0.5)])
    config["unconstrained_σ"] = 0.5
    return config
end

function get_scm_config()
    config = Dict()
    config["namelist_args"] = nothing
    return config
end


text:
export activation!
export
  CUDNN_ACTIVATION_SIGMOID, CUDNN_ACTIVATION_RELU,
  CUDNN_ACTIVATION_TANH, CUDNN_ACTIVATION_CLIPPED_RELU

function activation_desc(mode, relu_nanopt, relu_ceiling)
  p = Ptr{Void}[0]
  cudnnCreateActivationDescriptor(p)
  cudnnSetActivationDescriptor(p[1], mode, relu_nanopt, relu_ceiling)
  p[1]
end

"""
reluNanOpt: whether propagates NaN or not
reluCeiling: floating point number to specify the clipping threashod
when the activation mode is set to CUDNN_ACTIVATION_CLIPPED_RELU.
"""
function activation!{T}(mode, x::CuArray{T}, y::CuArray{T};
  relu_nanopt=CUDNN_NOT_PROPAGATE_NAN, relu_ceiling=1.0, alpha=1.0, beta=0.0)

  h = gethandle(x.dev)
  adesc = activation_desc(mode, relu_nanopt, relu_ceiling)
  xdesc = tensor_desc(x)
  ydesc = tensor_desc(y)
  cudnnActivationForward(h, adesc, T[alpha], xdesc, x, T[beta], ydesc, y)

  cudnnDestroyActivationDescriptor(adesc)
  cudnnDestroyTensorDescriptor(xdesc)
  cudnnDestroyTensorDescriptor(ydesc)
  y
end

function ∇activation!{T}(mode, y::CuArray{T}, dy::CuArray{T}, x::CuArray{T}, dx::CuArray{T};
  relu_nanopt=CUDNN_NOT_PROPAGATE_NAN, relu_ceiling=1.0, alpha=1.0, beta=0.0)

  h = gethandle(x.dev)
  adesc = activation_desc(mode, relu_nanopt, relu_ceiling)
  ydesc = tensor_desc(y)
  dydesc = tensor_desc(y)
  xdesc = tensor_desc(x)
  dxdesc = tensor_desc(dx)
  cudnnActivationBackward(h, adesc, T[alpha], ydesc, y, dydesc, dy, xdesc, x,
    T[beta], dxdesc, dx)

  cudnnDestroyActivationDescriptor(adesc)
  cudnnDestroyTensorDescriptor(ydesc)
  cudnnDestroyTensorDescriptor(dydesc)
  cudnnDestroyTensorDescriptor(xdesc)
  cudnnDestroyTensorDescriptor(dxdesc)
  dx
end


text:
#!/usr/bin/env julia
include("rf_cfchange.jl")
imgdir = joinpath("..", "..", "images")
function pltse8015(age, sex, indyr=1980)
	close()
	plotsex(age, 93, sex, 1980, 2015, indyr)
	savefig(joinpath(imgdir, "Sv8015_$(sex)$(age)RfInd$(indyr).svg"))
	close()
end
pltse8015(:a55_64, 2)
pltse8015(:a55_64, 1)
pltse8015(:a55_64, 2, 2015)
pltse8015(:a75_84, 2)
pltse8015(:a75_84, 1)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
using Test
import GDAL
import ArchGDAL as AG

@testset "test_rasterband.jl" begin
    @testset "Test methods for rasterband" begin
        AG.read("data/utmsmall.tif") do dataset
            ds_result = """
            GDAL Dataset (Driver: GTiff/GeoTIFF)
            File(s):
              data/utmsmall.tif

            Dataset (width x height): 100 x 100 (pixels)
            Number of raster bands: 1
              [GA_ReadOnly] Band 1 (Gray): 100 x 100 (UInt8)
            """
            @test sprint(print, dataset) == ds_result
            rb = AG.getband(dataset, 1)
            sprint(print, rb) == """
            [GA_ReadOnly] Band 1 (Gray): 100 x 100 (UInt8)
            blocksize: 100×81, nodata: -1.0e10, units: 1.0px + 0.0
            overviews: """
            @test sprint(print, AG.getdataset(rb)) == ds_result

            @test AG.getunittype(rb) == ""
            AG.setunittype!(rb, "ft")
            @test AG.getunittype(rb) == "ft"
            AG.setunittype!(rb, "")
            @test AG.getunittype(rb) == ""

            @test AG.getoffset(rb) == 0
            AG.setoffset!(rb, 10)
            @test AG.getoffset(rb) ≈ 10
            AG.setoffset!(rb, 0)
            @test AG.getoffset(rb) ≈ 0

            @test AG.getscale(rb) == 1
            AG.setscale!(rb, 0.5)
            @test AG.getscale(rb) ≈ 0.5
            AG.setscale!(rb, 2)
            @test AG.getscale(rb) ≈ 2
            AG.setscale!(rb, 1)
            @test AG.getscale(rb) ≈ 1

            @test isnothing(AG.getnodatavalue(rb))
            AG.setnodatavalue!(rb, -100)
            @test AG.getnodatavalue(rb) ≈ -100
            AG.deletenodatavalue!(rb)
            @test isnothing(AG.getnodatavalue(rb))

            AG.copy(dataset) do dest
                destband = AG.getband(dest, 1)
                AG.copywholeraster!(rb, destband)
                @test sprint(print, destband) == """
                [GA_Update] Band 1 (Gray): 100 x 100 (UInt8)
                    blocksize: 100×81, nodata: nothing, units: 1.0px + 0.0
                    overviews: """
                @test AG.noverview(destband) == 0
                AG.buildoverviews!(dest, Cint[2, 4, 8])
                @test AG.noverview(destband) == 3
                @test sprint(print, destband) == """
                [GA_Update] Band 1 (Gray): 100 x 100 (UInt8)
                    blocksize: 100×81, nodata: nothing, units: 1.0px + 0.0
                    overviews: (0) 50x50 (1) 25x25 (2) 13x13
                               """
                @test AG.getcolorinterp(destband) == AG.GCI_GrayIndex
                AG.setcolorinterp!(destband, AG.GCI_RedBand)
                @test AG.getcolorinterp(destband) == AG.GCI_RedBand

                @test sprint(print, AG.sampleoverview(destband, 100)) == """
                [GA_Update] Band 1 (Gray): 13 x 13 (UInt8)
                    blocksize: 128×128, nodata: nothing, units: 1.0px + 0.0
                    overviews: """
                @test sprint(print, AG.sampleoverview(destband, 200)) == """
                [GA_Update] Band 1 (Gray): 25 x 25 (UInt8)
                    blocksize: 128×128, nodata: nothing, units: 1.0px + 0.0
                    overviews: """
                @test sprint(print, AG.sampleoverview(destband, 500)) == """
                [GA_Update] Band 1 (Gray): 25 x 25 (UInt8)
                    blocksize: 128×128, nodata: nothing, units: 1.0px + 0.0
                    overviews: """
                AG.sampleoverview(destband, 1000) do result
                    @test sprint(print, result) == """
                    [GA_Update] Band 1 (Gray): 50 x 50 (UInt8)
                        blocksize: 128×128, nodata: nothing, units: 1.0px + 0.0
                        overviews: """
                end
                @test sprint(print, AG.getmaskband(destband)) == """
                [GA_ReadOnly] Band 0 (Undefined): 100 x 100 (UInt8)
                    blocksize: 100×81, nodata: nothing, units: 1.0px + 0.0
                    overviews: """
                @test AG.maskflags(destband) == 1
                @test AG.maskflaginfo(rb) == (
                    all_valid = true,
                    per_dataset = false,
                    alpha = false,
                    nodata = false,
                )
                AG.createmaskband!(destband, 3)
                AG.getmaskband(destband) do maskband
                    @test sprint(print, maskband) == """
                    [GA_Update] Band 1 (Gray): 100 x 100 (UInt8)
                        blocksize: 100×81, nodata: nothing, units: 1.0px + 0.0
                        overviews: """
                end
                @test AG.maskflags(destband) == 3
                @test AG.maskflaginfo(destband) == (
                    all_valid = true,
                    per_dataset = true,
                    alpha = false,
                    nodata = false,
                )
                AG.fillraster!(destband, 3)
                AG.setcategorynames!(destband, ["foo", "bar"])
                @test AG.getcategorynames(destband) == ["foo", "bar"]

                AG.getoverview(destband, 0) do overview
                    return AG.regenerateoverviews!(
                        destband,
                        [overview, AG.getoverview(destband, 2)],
                    )
                end

                AG.createRAT() do rat
                    AG.setdefaultRAT!(destband, rat)
                    @test AG.getdefaultRAT(destband).ptr != rat.ptr
                end
                @test AG.getdefaultRAT(destband).ptr !=
                      GDAL.GDALRasterAttributeTableH(C_NULL)

                AG.getcolortable(destband) do ct
                    @test ct.ptr == GDAL.GDALColorTableH(C_NULL)
                end
                AG.createcolortable(AG.GPI_RGB) do ct
                    @test ct.ptr != GDAL.GDALColorTableH(C_NULL)
                    return AG.setcolortable!(destband, ct)
                end
                AG.clearcolortable!(destband)
                AG.getcolortable(destband) do ct
                    @test ct.ptr == GDAL.GDALColorTableH(C_NULL)
                end
            end
        end
    end
end


text:
using Gadfly, Compose

set_default_plot_size(12inch,3inch)

p=plot(y=[1,2,3], Geom.line)
c=Compose.context()
hstack(Union{Plot,Compose.Context}[p,c,p,c,p])


text:
@require "../"

suite("should parse milliseconds") do
  parsed = parse(100)
  @test parsed["milliseconds"] == 100
end

suite("should parse seconds") do
  parsed = parse(1000)
  @test parsed["seconds"] == 1
end

suite("should parse minutes") do
  parsed = parse(60000)
  @test parsed["minutes"] == 1
end

suite("should parse hours") do
  parsed = parse(3600000)
  @test parsed["hours"] == 1
end

suite("should parse days") do
  parsed = parse(86400000)
  @test parsed["days"] == 1
end

max_stars_repo_name: queckezz/ms-parse

max_issues_repo_name: queckezz/ms-parse

max_forks_repo_name: queckezz/ms-parse


text:
using Turing
using Test
using Random
Random.seed!(125)

x = [1.5 2.0]

@model mhtest(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  for i in 1:length(x)
    x[i] ~ Normal(m, sqrt(s))
  end
  s, m
end

# MH with prior as proposal
alg = MH(2000)
chain = sample(mhtest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.1
@test mean(chain[:m]) ≈ 7/6 atol=0.1

# MH with Gaussian proposal
GKernel(var) = (x) -> Normal(x, sqrt.(var))
alg = MH(5000, (:s, GKernel(5)), (:m, GKernel(1.0)))
chain = sample(mhtest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.2
@test mean(chain[:m]) ≈ 7/6 atol=0.1

# MH within Gibbs
alg = Gibbs(1000, MH(5, :m), MH(5, :s))
chain = sample(mhtest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.1
@test mean(chain[:m]) ≈ 7/6 atol=0.1


text:
using AutoGrid
using Test

@testset "AutoGrid.jl" begin
    # Write your tests here.
end


text:
function Subplot{T<:AbstractBackend}(::T; parent = RootLayout())
    Subplot{T}(
        parent,
        Series[],
        (20mm, 5mm, 2mm, 10mm),
        defaultbox,
        defaultbox,
        KW(),
        nothing,
        nothing
    )
end

plotarea(sp::Subplot) = sp.plotarea
plotarea!(sp::Subplot, bbox::BoundingBox) = (sp.plotarea = bbox)

Base.size(sp::Subplot) = (1,1)
Base.length(sp::Subplot) = 1
Base.getindex(sp::Subplot, r::Int, c::Int) = sp

leftpad(sp::Subplot)   = sp.minpad[1]
toppad(sp::Subplot)    = sp.minpad[2]
rightpad(sp::Subplot)  = sp.minpad[3]
bottompad(sp::Subplot) = sp.minpad[4]

get_subplot(plt::Plot, sp::Subplot) = sp
get_subplot(plt::Plot, i::Integer) = plt.subplots[i]
get_subplot(plt::Plot, k) = plt.spmap[k]
get_subplot(series::Series) = series.d[:subplot]

get_subplot_index(plt::Plot, idx::Integer) = Int(idx)
get_subplot_index(plt::Plot, sp::Subplot) = findfirst(_ -> _ === sp, plt.subplots)

series_list(sp::Subplot) = sp.series_list # filter(series -> series.d[:subplot] === sp, sp.plt.series_list)

function should_add_to_legend(series::Series)
    series.d[:primary] && series.d[:label] != "" &&
        !(series.d[:seriestype] in (
            :hexbin,:histogram2d,:hline,:vline,
            :contour,:contourf,:contour3d,:surface,:wireframe,
            :heatmap, :pie, :image
        ))
end

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


text:
"""
Displays and image in GIMP

# Example:
view_img(rand(1:10,10,10))
"""
function view_img(img)
    save_img("/tmp/image.tif", img)
    run(`gimp /tmp/image.tif`)
end

## More will be added to this library. Missinngs can be marked with red.
## Boundary layers can be added.


text:
# syntax: proto2
using ProtoBuf
import ProtoBuf.meta

mutable struct Blob <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Blob(; kwargs...)
        obj = new(meta(Blob), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Blob
const __meta_Blob = Ref{ProtoMeta}()
function meta(::Type{Blob})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Blob)
            __meta_Blob[] = target = ProtoMeta(Blob)
            allflds = Pair{Symbol,Union{Type,String}}[:raw => Vector{UInt8}, :raw_size => Int32, :zlib_data => Vector{UInt8}, :lzma_data => Vector{UInt8}, :OBSOLETE_bzip2_data => Vector{UInt8}]
            meta(target, Blob, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Blob[]
    end
end
function Base.getproperty(obj::Blob, name::Symbol)
    if name === :raw
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :raw_size
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :zlib_data
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :lzma_data
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :OBSOLETE_bzip2_data
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct BlobHeader <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BlobHeader(; kwargs...)
        obj = new(meta(BlobHeader), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct BlobHeader
const __meta_BlobHeader = Ref{ProtoMeta}()
function meta(::Type{BlobHeader})
    ProtoBuf.metalock() do
        if !isassigned(__meta_BlobHeader)
            __meta_BlobHeader[] = target = ProtoMeta(BlobHeader)
            req = Symbol[:_type,:datasize]
            allflds = Pair{Symbol,Union{Type,String}}[:_type => AbstractString, :indexdata => Vector{UInt8}, :datasize => Int32]
            meta(target, BlobHeader, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_BlobHeader[]
    end
end
function Base.getproperty(obj::BlobHeader, name::Symbol)
    if name === :_type
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :indexdata
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :datasize
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export Blob, BlobHeader


text:
sum = 0
for i in 1:1000000000
    sum += i
end
print(sum)


text:
function add_expressions!(
    container::OptimizationContainer,
    ::Type{T},
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <: ExpressionType,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    add_expression_container!(container, T(), D, names, time_steps)
    return
end

function add_expressions!(
    container::OptimizationContainer,
    ::Type{T},
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <:
    Union{ComponentActivePowerRangeExpressionUB, ComponentActivePowerRangeExpressionLB},
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.HybridSystem}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    subcomp_keys = string.([PSY.ThermalGen, PSY.RenewableGen])
    add_expression_container!(
        container,
        T(),
        D,
        names,
        subcomp_keys,
        time_steps;
        sparse=true,
    )
    return
end

function add_expressions!(
    container::OptimizationContainer,
    ::Type{T},
    devices::U,
    model::ServiceModel{V, W},
) where {
    T <: ExpressionType,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    V <: PSY.Reserve,
    W <: AbstractReservesFormulation,
} where {D <: PSY.Component}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    add_expression_container!(container, T(), D, names, time_steps)
    return
end

# Note: add_to_jump_expression! are legacy when more control was needed over the calls to
# add_to_expression. These might be removed post JuMP 1.0 release.
function _add_to_jump_expression!(
    expression::T,
    var::JuMP.VariableRef,
    multiplier::Float64,
) where {T <: JuMP.AbstractJuMPScalar}
    JuMP.add_to_expression!(expression, multiplier, var)
    return
end

function _add_to_jump_expression!(
    expression::T,
    value::Float64,
) where {T <: JuMP.AbstractJuMPScalar}
    JuMP.add_to_expression!(expression, value)
    return
end

function _add_to_jump_expression!(
    expression::T,
    parameter::PJ.ParameterRef,
    multiplier::Float64,
) where {T <: JuMP.AbstractJuMPScalar}
    PJ.add_to_expression!(expression, multiplier, parameter)
    return
end

function _add_to_jump_expression!(
    expression::T,
    var::Union{JuMP.VariableRef, PJ.ParameterRef},
    multiplier::Float64,
    constant::Float64,
) where {T <: JuMP.AbstractJuMPScalar}
    _add_to_jump_expression!(expression, constant)
    _add_to_jump_expression!(expression, var, multiplier)
    return
end

function _add_to_jump_expression!(
    expression::T,
    parameter::Float64,
    multiplier::Float64,
) where {T <: JuMP.AbstractJuMPScalar}
    _add_to_jump_expression!(expression, parameter * multiplier)
    return
end

"""
Default implementation to add parameters to SystemBalanceExpressions
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: SystemBalanceExpressions,
    U <: TimeSeriesParameter,
    V <: PSY.Device,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    parameter = get_parameter_array(container, U(), V)
    multiplier = get_parameter_multiplier_array(container, U(), V)
    for d in devices, t in get_time_steps(container)
        bus_number = PSY.get_number(PSY.get_bus(d))
        name = PSY.get_name(d)
        _add_to_jump_expression!(
            get_expression(container, T(), X)[bus_number, t],
            parameter[name, t],
            multiplier[name, t],
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: OnStatusParameter,
    V <: PSY.ThermalGen,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    parameter = get_parameter_array(container, U(), V)

    for d in devices, t in get_time_steps(container)
        bus_number = PSY.get_number(PSY.get_bus(d))
        name = PSY.get_name(d)
        mult = get_expression_multiplier(U(), T(), d, W())
        _add_to_jump_expression!(
            get_expression(container, T(), X)[bus_number, t],
            parameter[name, t],
            mult,
        )
    end
    return
end

"""
Default implementation to add variables to SystemBalanceExpressions
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: SystemBalanceExpressions,
    U <: VariableType,
    V <: PSY.Device,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    variable = get_variable(container, U(), V)
    expression = get_expression(container, T(), X)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        bus_number = PSY.get_number(PSY.get_bus(d))
        _add_to_jump_expression!(
            expression[bus_number, t],
            variable[name, t],
            get_variable_multiplier(U(), V, W()),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: SystemBalanceExpressions,
    U <: OnVariable,
    V <: PSY.ThermalGen,
    W <: Union{AbstractCompactUnitCommitment, ThermalCompactDispatch},
    X <: PM.AbstractPowerModel,
}
    variable = get_variable(container, U(), V)
    expression = get_expression(container, T(), X)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        bus_number = PSY.get_number(PSY.get_bus(d))
        _add_to_jump_expression!(
            expression[bus_number, t],
            variable[name, t],
            get_variable_multiplier(U(), d, W()),
        )
    end
    return
end

"""
Default implementation to add parameters to SystemBalanceExpressions
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: SystemBalanceExpressions,
    U <: TimeSeriesParameter,
    V <: PSY.StaticInjection,
    W <: AbstractDeviceFormulation,
    X <: CopperPlatePowerModel,
}
    parameter = get_parameter_array(container, U(), V)
    multiplier = get_parameter_multiplier_array(container, U(), V)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        _add_to_jump_expression!(
            get_expression(container, T(), X)[t],
            parameter[name, t],
            multiplier[name, t],
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: OnStatusParameter,
    V <: PSY.ThermalGen,
    W <: AbstractDeviceFormulation,
    X <: CopperPlatePowerModel,
}
    parameter = get_parameter_array(container, U(), V)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        mult = get_expression_multiplier(U(), T(), d, W())
        _add_to_jump_expression!(
            get_expression(container, T(), X)[t],
            parameter[name, t],
            mult,
        )
    end
    return
end

"""
Default implementation to add variables to SystemBalanceExpressions
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: VariableType,
    V <: PSY.StaticInjection,
    W <: AbstractDeviceFormulation,
    X <: CopperPlatePowerModel,
}
    variable = get_variable(container, U(), V)
    expression = get_expression(container, T(), X)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        _add_to_jump_expression!(
            expression[t],
            variable[name, t],
            get_variable_multiplier(U(), V, W()),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: OnVariable,
    V <: PSY.ThermalGen,
    W <: Union{AbstractCompactUnitCommitment, ThermalCompactDispatch},
    X <: CopperPlatePowerModel,
}
    variable = get_variable(container, U(), V)
    expression = get_expression(container, T(), X)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        _add_to_jump_expression!(
            expression[t],
            variable[name, t],
            get_variable_multiplier(U(), d, W()),
        )
    end
    return
end

"""
Default implementation to add parameters to SystemBalanceExpressions
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: SystemBalanceExpressions,
    U <: TimeSeriesParameter,
    V <: PSY.StaticInjection,
    W <: AbstractDeviceFormulation,
    X <: Union{PTDFPowerModel, StandardPTDFModel},
}
    parameter = get_parameter_array(container, U(), V)
    multiplier = get_parameter_multiplier_array(container, U(), V)
    sys_expr = get_expression(container, T(), PSY.System)
    nodal_expr = get_expression(container, T(), PSY.Bus)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        bus_no = PSY.get_number(PSY.get_bus(d))
        _add_to_jump_expression!(sys_expr[t], parameter[name, t], multiplier[name, t])
        _add_to_jump_expression!(
            nodal_expr[bus_no, t],
            parameter[name, t],
            multiplier[name, t],
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: OnStatusParameter,
    V <: PSY.ThermalGen,
    W <: AbstractDeviceFormulation,
    X <: Union{PTDFPowerModel, StandardPTDFModel},
}
    parameter = get_parameter_array(container, U(), V)
    sys_expr = get_expression(container, T(), PSY.System)
    nodal_expr = get_expression(container, T(), PSY.Bus)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        bus_no = PSY.get_number(PSY.get_bus(d))
        mult = get_expression_multiplier(U(), T(), d, W())
        _add_to_jump_expression!(sys_expr[t], parameter[name, t], mult)
        _add_to_jump_expression!(nodal_expr[bus_no, t], parameter[name, t], mult)
    end
    return
end

"""
Default implementation to add variables to SystemBalanceExpressions
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: VariableType,
    V <: PSY.StaticInjection,
    W <: AbstractDeviceFormulation,
    X <: Union{PTDFPowerModel, StandardPTDFModel},
}
    variable = get_variable(container, U(), V)
    sys_expr = get_expression(container, T(), PSY.System)
    nodal_expr = get_expression(container, T(), PSY.Bus)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        bus_no = PSY.get_number(PSY.get_bus(d))
        _add_to_jump_expression!(
            sys_expr[t],
            variable[name, t],
            get_variable_multiplier(U(), V, W()),
        )
        _add_to_jump_expression!(
            nodal_expr[bus_no, t],
            variable[name, t],
            get_variable_multiplier(U(), V, W()),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: OnVariable,
    V <: PSY.ThermalGen,
    W <: Union{AbstractCompactUnitCommitment, ThermalCompactDispatch},
    X <: Union{PTDFPowerModel, StandardPTDFModel},
}
    variable = get_variable(container, U(), V)
    sys_expr = get_expression(container, T(), PSY.System)
    nodal_expr = get_expression(container, T(), PSY.Bus)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        bus_no = PSY.get_number(PSY.get_bus(d))
        _add_to_jump_expression!(
            sys_expr[t],
            variable[name, t],
            get_variable_multiplier(U(), d, W()),
        )
        _add_to_jump_expression!(
            nodal_expr[bus_no, t],
            variable[name, t],
            get_variable_multiplier(U(), d, W()),
        )
    end
    return
end

"""
Implementation of add_to_expression! for lossless branch/network models
"""
function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: FlowActivePowerVariable,
    V <: PSY.Branch,
    W <: AbstractBranchFormulation,
    X <: PM.AbstractActivePowerModel,
}
    var = get_variable(container, U(), V)
    expression = get_expression(container, T(), X)
    for d in devices
        for t in get_time_steps(container)
            flow_variable = var[PSY.get_name(d), t]
            _add_to_jump_expression!(
                expression[PSY.get_number(PSY.get_arc(d).from), t],
                flow_variable,
                -1.0,
            )
            _add_to_jump_expression!(
                expression[PSY.get_number(PSY.get_arc(d).to), t],
                flow_variable,
                1.0,
            )
        end
    end
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ActivePowerBalance,
    U <: HVDCTotalPowerDeliveredVariable,
    V <: PSY.DCBranch,
    W <: AbstractBranchFormulation,
    X <: PM.AbstractPowerModel,
}
    var = get_variable(container, U(), V)
    expression = get_expression(container, T(), X)
    for d in devices
        for t in get_time_steps(container)
            flow_variable = var[PSY.get_name(d), t]
            _add_to_jump_expression!(
                expression[PSY.get_number(PSY.get_arc(d).to), t],
                flow_variable,
                1.0,
            )
        end
    end
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: Union{ActivePowerRangeExpressionUB, ActivePowerRangeExpressionLB},
    U <: VariableType,
    V <: PSY.Device,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    variable = get_variable(container, U(), V)
    if !has_container_key(container, T, V)
        add_expressions!(container, T, devices, model)
    end
    expression = get_expression(container, T(), V)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        _add_to_jump_expression!(expression[name, t], variable[name, t], 1.0)
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <:
    Union{ComponentActivePowerRangeExpressionUB, ComponentActivePowerRangeExpressionLB},
    U <: ComponentActivePowerVariable,
    V <: PSY.HybridSystem,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    variable = get_variable(container, U(), V)
    if !has_container_key(container, T, V)
        add_expressions!(container, T, devices, model)
    end
    expression = get_expression(container, T(), V)
    for d in devices,
        t in get_time_steps(container),
        sub_comp in [PSY.ThermalGen, PSY.RenewableGen]

        sub_comp_key = string(sub_comp)
        name = PSY.get_name(d)
        _add_to_jump_expression!(
            expression[name, sub_comp_key, t],
            variable[name, sub_comp_key, t],
            1.0,
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}},
    model::ServiceModel{X, W},
) where {
    T <: Union{ActivePowerRangeExpressionUB, ReserveRangeExpressionUB},
    U <: VariableType,
    V <: PSY.Component,
    X <: PSY.Reserve{PSY.ReserveUp},
    W <: AbstractReservesFormulation,
}
    service_name = get_service_name(model)
    variable = get_variable(container, U(), X, service_name)
    if !has_container_key(container, T, V)
        add_expressions!(container, T, devices, model)
    end
    expression = get_expression(container, T(), V)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        _add_to_jump_expression!(expression[name, t], variable[name, t], 1.0)
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::Union{Vector{V}, IS.FlattenIteratorWrapper{V}},
    model::ServiceModel{X, W},
) where {
    T <: Union{ActivePowerRangeExpressionLB, ReserveRangeExpressionLB},
    U <: VariableType,
    V <: PSY.Component,
    X <: PSY.Reserve{PSY.ReserveDown},
    W <: AbstractReservesFormulation,
}
    service_name = get_service_name(model)
    variable = get_variable(container, U(), X, service_name)
    if !has_container_key(container, T, V)
        add_expressions!(container, T, devices, model)
    end
    expression = get_expression(container, T(), V)
    for d in devices, t in get_time_steps(container)
        name = PSY.get_name(d)
        _add_to_jump_expression!(expression[name, t], variable[name, t], -1.0)
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::U,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {
    T <: Union{ActivePowerRangeExpressionUB, ActivePowerRangeExpressionLB},
    U <: OnStatusParameter,
    V <: PSY.Device,
    W <: AbstractDeviceFormulation,
}
    parameter_array = get_parameter_array(container, U(), V)
    if !has_container_key(container, T, V)
        add_expressions!(container, T, devices, model)
    end
    expression = get_expression(container, T(), V)
    for d in devices, mult in get_expression_multiplier(U(), T(), d, W())
        for t in get_time_steps(container)
            name = PSY.get_name(d)
            _add_to_jump_expression!(
                expression[name, t],
                parameter_array[name, t],
                -mult,
                mult,
            )
        end
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::U,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {
    T <: Union{ActivePowerRangeExpressionUB, ActivePowerRangeExpressionLB},
    U <: OnStatusParameter,
    V <: PSY.ThermalGen,
    W <: AbstractThermalDispatchFormulation,
}
    parameter_array = get_parameter_array(container, U(), V)
    if !has_container_key(container, T, V)
        add_expressions!(container, T, devices, model)
    end
    expression = get_expression(container, T(), V)
    for d in devices, mult in get_expression_multiplier(U(), T(), d, W())
        for t in get_time_steps(container)
            name = PSY.get_name(d)
            _add_to_jump_expression!(expression[name, t], parameter_array[name, t], -mult)
        end
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{U},
    model::ServiceModel{V, W},
    devices_template::Dict{Symbol, DeviceModel},
) where {U <: VariableType, V <: PSY.Reserve, W <: AbstractReservesFormulation}
    contributing_devices_map = get_contributing_devices_map(model)
    for (device_type, devices) in contributing_devices_map
        device_model = get(devices_template, Symbol(device_type), nothing)
        device_model === nothing && continue
        expression_type = get_expression_type_for_reserve(U(), device_type, V)
        add_to_expression!(container, expression_type, U, devices, model)
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    ::PSY.System,
    ::NetworkModel{W},
    ::Type{W},
) where {
    T <: ActivePowerBalance,
    U <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    W <: Union{CopperPlatePowerModel, StandardPTDFModel},
}
    variable = get_variable(container, U(), PSY.System)
    expression = get_expression(container, T(), PSY.System)
    for t in get_time_steps(container)
        _add_to_jump_expression!(
            expression[t],
            variable[t],
            get_variable_multiplier(U(), PSY.System, W()),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    sys::PSY.System,
    ::NetworkModel{W},
    ::Type{W},
) where {
    T <: ActivePowerBalance,
    U <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    W <: PM.AbstractActivePowerModel,
}
    variable = get_variable(container, U(), PSY.Bus)
    expression = get_expression(container, T(), PSY.Bus)
    bus_numbers = PSY.get_number.(PSY.get_components(PSY.Bus, sys))
    for t in get_time_steps(container), n in bus_numbers
        _add_to_jump_expression!(
            expression[n, t],
            variable[n, t],
            get_variable_multiplier(U(), PSY.Bus, W),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    sys::PSY.System,
    ::NetworkModel{W},
    ::Type{W},
) where {
    T <: ActivePowerBalance,
    U <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    W <: PM.AbstractPowerModel,
}
    variable = get_variable(container, U(), PSY.Bus, "P")
    expression = get_expression(container, T(), PSY.Bus)
    bus_numbers = PSY.get_number.(PSY.get_components(PSY.Bus, sys))
    for t in get_time_steps(container), n in bus_numbers
        _add_to_jump_expression!(
            expression[n, t],
            variable[n, t],
            get_variable_multiplier(U(), PSY.Bus, W),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    sys::PSY.System,
    ::NetworkModel{W},
    ::Type{W},
) where {
    T <: ReactivePowerBalance,
    U <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    W <: PM.AbstractPowerModel,
}
    variable = get_variable(container, U(), PSY.Bus, "Q")
    expression = get_expression(container, T(), PSY.Bus)
    bus_numbers = PSY.get_number.(PSY.get_components(PSY.Bus, sys))
    for t in get_time_steps(container), n in bus_numbers
        _add_to_jump_expression!(
            expression[n, t],
            variable[n, t],
            get_variable_multiplier(U(), PSY.Bus, W),
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{S},
    cost_expression::JuMP.AbstractJuMPScalar,
    component::T,
    time_period::Int,
) where {S <: CostExpressions, T <: PSY.Component}
    if has_container_key(container, S, T)
        device_cost_expression = get_expression(container, S(), T)
        component_name = PSY.get_name(component)
        JuMP.add_to_expression!(
            device_cost_expression[component_name, time_period],
            cost_expression,
        )
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    areas::IS.FlattenIteratorWrapper{V},
    model::ServiceModel{PSY.AGC, W},
) where {
    T <: Union{EmergencyUp, EmergencyDown},
    U <:
    Union{AdditionalDeltaActivePowerUpVariable, AdditionalDeltaActivePowerDownVariable},
    V <: PSY.Area,
    W <: AbstractServiceFormulation,
}
    names = PSY.get_name.(areas)
    time_steps = get_time_steps(container)
    if !has_container_key(container, T, V)
        expression = add_expression_container!(container, T(), V, names, time_steps)
    end
    expression = get_expression(container, T(), V)
    variable = get_variable(container, U(), V)
    for n in names, t in time_steps
        _add_to_jump_expression!(expression[n, t], variable[n, t], 1.0)
    end
    return
end

function add_to_expression!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    services::IS.FlattenIteratorWrapper{V},
    model::ServiceModel{V, W},
) where {
    T <: RawACE,
    U <: SteadyStateFrequencyDeviation,
    V <: PSY.AGC,
    W <: AbstractServiceFormulation,
}
    names = [PSY.get_name(PSY.get_area(s)) for s in services]
    time_steps = get_time_steps(container)
    if !has_container_key(container, T, V)
        expression = add_expression_container!(container, T(), PSY.Area, names, time_steps)
    end
    expression = get_expression(container, T(), PSY.Area)
    variable = get_variable(container, U(), PSY.Area)
    for s in services, t in time_steps
        name = PSY.get_name(PSY.get_area(s))
        _add_to_jump_expression!(
            expression[name, t],
            variable[t],
            get_variable_multiplier(U(), s, W()),
        )
    end
    return
end


text:
using KLIEPInference
using Random
using Distributions
using JLD

Random.seed!(325)

## chain (1)
idx0 = KLIEPInference.trimap(5, 6)
idx1 = KLIEPInference.trimap(4, 5)
idx2 = KLIEPInference.trimap(6, 7)
idx3 = KLIEPInference.trimap(4, 6)
idx4 = KLIEPInference.trimap(5, 7)

for m in [25, 50]
    γx = chain(m; lenC=m, sgn=0)
    γy = deepcopy(γx)
    γy[idx0] = abs(γy[idx0]) > 0.2 ? γy[idx0] - sign(γy[idx0])*0.1 : γy[idx0] + sign(γy[idx0])*0.1
    γy[idx1] *= -1.0
    γy[idx2] *= -1.0
    γy[idx3] += (2. * rand(Bernoulli()) - 1.) * 0.1
    γy[idx4] += (2. * rand(Bernoulli()) - 1.) * 0.1
    @save "graphs/params_exp1_chain1_$(m).jld" γx γy
end

## chain (2)
idx0 = KLIEPInference.trimap(5, 6)
idx1 = KLIEPInference.trimap(3, 4)
idx2 = KLIEPInference.trimap(7, 8)
idx3 = KLIEPInference.trimap(4, 5)
idx4 = KLIEPInference.trimap(4, 6)

for m in [25, 50]
    γx = chain(m; lenC=m, sgn=0)
    γy = deepcopy(γx)
    γy[idx0] = abs(γy[idx0]) > 0.2 ? γy[idx0] - sign(γy[idx0])*0.1 : γy[idx0] + sign(γy[idx0])*0.1
    γy[idx1] *= -1.0
    γy[idx2] *= -1.0
    γy[idx3] = abs(γy[idx3]) > 0.2 ? γy[idx3] - sign(γy[idx3])*0.1 : γy[idx3] + sign(γy[idx3])*0.1
    γy[idx4] += (2. * rand(Bernoulli()) - 1.) * 0.1
    @save "graphs/params_exp1_chain2_$(m).jld" γx γy
end

Random.seed!(324)

## tree (1)
idx0 = KLIEPInference.trimap(1, 3)
idx1 = KLIEPInference.trimap(1, 2)
idx2 = KLIEPInference.trimap(3, 8)
idx3 = KLIEPInference.trimap(1, 9)
idx4 = KLIEPInference.trimap(3, 4)

for m in [25, 50]
    γx = kary_tree(m, sgn=0)
    γy = deepcopy(γx)
    γy[idx0] = abs(γy[idx0]) > 0.2 ? γy[idx0] - sign(γy[idx0])*0.1 : γy[idx0] + sign(γy[idx0])*0.1
    γy[idx1] *= -1.0
    γy[idx2] *= -1.0
    γy[idx3] += (2. * rand(Bernoulli()) - 1.) * 0.1
    γy[idx4] = abs(γy[idx4]) > 0.2 ? γy[idx4] - sign(γy[idx4])*0.1 : γy[idx4] + sign(γy[idx4])*0.1
    @save "graphs/params_exp1_tree1_$(m).jld" γx γy
end

## tree (2)
idx0 = KLIEPInference.trimap(1, 3)
idx1 = KLIEPInference.trimap(2, 5)
idx2 = KLIEPInference.trimap(8, 23)
idx3 = KLIEPInference.trimap(2, 3)
idx4 = KLIEPInference.trimap(3, 4)

for m in [25, 50]
    γx = kary_tree(m, sgn=0)
    γy = deepcopy(γx)
    γy[idx0] = abs(γy[idx0]) > 0.2 ? γy[idx0] - sign(γy[idx0])*0.1 : γy[idx0] + sign(γy[idx0])*0.1
    γy[idx1] *= -1.0
    γy[idx2] *= -1.0
    γy[idx3] += (2. * rand(Bernoulli()) - 1.) * 0.1
    γy[idx4] += (2. * rand(Bernoulli()) - 1.) * 0.1
    @save "graphs/params_exp1_tree2_$(m).jld" γx γy
end


text:
# Ilija Spasić PR 34-2020

# Predmetni projekat - Bucket Sort algoritam

# Zadatak: Neka je dat niz brojeva iz intervala [1..100].
# Sortirati ovaj niz Bucket Sort metodom i ispisati tako dobijeni rezultat na ekran.
#
#     Skripta je dužna da kreira dati niz sa bar 30 brojeva,
# pozove funkciju BucketSort() koja sortira ovaj niz i ispiše ga na ekran.
# Funkcija BucketSort() sortira niz Bucket Sort metodom, i vraća sortirani niz.

using Printf

include("bucketSort.jl")

a = 1
b = 100

X = generateArray(a,b)

sortedX = bucketSort(X,"r")

println("Niz pre sortiranja:\t Niz posle sortiranja: \n")

for i in 1:length(X)
    @printf("Niz[%2i] = %3i \t\t Niz[%2i] = %3i\n",
            i, X[i], i, sortedX[i])
end