text:
facts("Ranges") do

  context("range shatterinto range") do
    @fact shatterinto(r1, r1) => Range{Int}[]
    @fact shatterinto(r1, r2) => Range{Int}[]
    @fact shatterinto(r1, r3) => [Range(Box{Int}(Point(1,10), Point(12,12)), 3,3)
                                      Range(Box{Int}(Point(10,1), Point(12,10)), 3,3)]
    @fact shatterinto(r1, r4) => [r4]
    @fact shatterinto(r1, r5) => [r5]
    @fact shatterinto(r1, r6) => [Range(Box{Int}(Point(-1,10), Point(13,14)), 6,6),
                                      Range(Box{Int}(Point(-1,-2), Point(13,0)), 6,6),
                                      Range(Box{Int}(Point(10,0),  Point(13,10)), 6,6),
                                      Range(Box{Int}(Point(-1,0),  Point(0,10)), 6,6)]
  end

  context("push layout") do
    @fact push(Layout([r2]), r1) => Layout{Int64}([r1])
    @fact push(Layout([r1]), r2) => Layout{Int64}(
                    [Range{Int64}(Box{Int64}(Point{Int64}(0,2),Point{Int64}(10,10)),1,1),
                     Range{Int64}(Box{Int64}(Point{Int64}(0,0),Point{Int64}(10,1)),1,1),
                     Range{Int64}(Box{Int64}(Point{Int64}(2,1),Point{Int64}(10,2)),1,1),
                     Range{Int64}(Box{Int64}(Point{Int64}(0,1),Point{Int64}(1,2)),1,1),
                     Range{Int64}(Box{Int64}(Point{Int64}(1,1),Point{Int64}(2,2)),2,2)])
    @fact push(Layout([r2]), r1) => Layout{Int64}(
      [Range{Int64}(Box{Int64}(Point{Int64}(0,0),Point{Int64}(10,10)),1,1)])

    @fact push(push(Layout([r1]), r2), r5) => Layout{Int64}(
      [Range{Int64}(Box{Int64}(Point{Int64}(0,2),Point{Int64}(10,10)),1,1),
       Range{Int64}(Box{Int64}(Point{Int64}(0,0),Point{Int64}(10,1)),1,1),
       Range{Int64}(Box{Int64}(Point{Int64}(2,1),Point{Int64}(10,2)),1,1),
       Range{Int64}(Box{Int64}(Point{Int64}(0,1),Point{Int64}(1,2)),1,1),
       Range{Int64}(Box{Int64}(Point{Int64}(1,1),Point{Int64}(2,2)),2,2),
       Range{Int64}(Box{Int64}(Point{Int64}(5,10),Point{Int64}(12,12)),5,5)])
  end
end


text:
# This is data_types.jl
# contains type definitions to be used in GeophysicalModelGenerator

import Base: show

export  GeoData, ParaviewData, UTMData, CartData,
        LonLatDepthGrid, XYZGrid, Velocity_SphericalToCartesian!,
        Convert2UTMzone, Convert2CartData, ProjectionPoint

"""
    struct ProjectionPoint
        Lon     :: Float64
        Lat     :: Float64
        EW      :: Float64
        NS      :: Float64
        zone    :: Integer
        isnorth :: Bool
    end

Structure that holds the coordinates of a point that is used to project a data set from Lon/Lat to a Cartesian grid and vice-versa.
"""
struct ProjectionPoint
    Lat     :: Float64
    Lon     :: Float64
    EW      :: Float64
    NS      :: Float64
    zone    :: Integer
    isnorth :: Bool
end

"""
    ProjectionPoint(; Lat=49.9929, Lon=8.2473)

Defines a projection point used for map projections, by specifying latitude and longitude
"""
function ProjectionPoint(; Lat=49.9929, Lon=8.2473)
    # Default = Mainz (center of universe)
    x_lla = LLA(Lat, Lon, 0.0);    # Lat/Lon/Alt of geodesy package
    x_utmz = UTMZ(x_lla, wgs84)    # UTMZ of

    ProjectionPoint(Lat, Lon, x_utmz.x, x_utmz.y, Int64(x_utmz.zone), x_utmz.isnorth)
end

"""
    ProjectionPoint(EW::Float64, NS::Float64, Zone::Int64, isnorth::Bool)

Defines a projection point used for map projections, by specifying UTM coordinates (EW/NS), UTM Zone and whether you are on the northern hemisphere

"""
function ProjectionPoint(EW::Float64, NS::Float64, Zone::Int64, isnorth::Bool)

    x_utmz = UTMZ(EW,NS,0.0,Zone, isnorth)    # UTMZ of
    x_lla = LLA(x_utmz, wgs84);    # Lat/Lon/Alt of geodesy package

    ProjectionPoint(x_lla.lat, x_lla.lon, EW, NS, Zone, isnorth)
end

# data structure for a list of values - TO BE REMOVED
mutable struct ValueList
    name::String
    unit::String
    values::Vector{Float64}
end

"""
    GeoData(lon::Any, lat:Any, depth::GeoUnit, fields::NamedTuple)

Data structure that holds one or several fields with longitude, latitude and depth information.

- `depth` can have units of meter, kilometer or be unitless; it will be converted to km.
- `fields` should ideally be a NamedTuple which allows you to specify the names of each of the fields.
- In case you only pass one array we will convert it to a NamedTuple with default name.
- A single field should be added as `(DataFieldName=Data,)` (don't forget the comma at the end).
- Multiple fields  can be added as well. `lon`,`lat`,`depth` should all have the same size as each of the `fields`.
- In case you want to display a vector field in paraview, add it as a tuple: `(Velocity=(Veast,Vnorth,Vup), Veast=Veast, Vnorth=Vnorth, Vup=Vup)`; we automatically apply a vector transformation when transforming this to a `ParaviewData` structure from which we generate Paraview output. As this changes the magnitude of the arrows, you will no longer see the `[Veast,Vnorth,Vup]` components in Paraview which is why it is a good ideas to store them as separate Fields.
- Yet, there is one exception: if the name of the 3-component field is `colors`, we do not apply this vector transformation as this field is regarded to contain RGB colors.
- `Lat`,`Lon`,`Depth` should have the same size as the `Data` array. The ordering of the arrays is important. If they are 3D arrays, as in the example below, we assume that the first dimension corresponds to `lon`, second dimension to `lat` and third dimension to `depth` (which should be in km). See below for an example.

# Example
```julia-repl
julia> Lat         =   1.0:3:10.0;
julia> Lon         =   11.0:4:20.0;
julia> Depth       =   (-20:5:-10)*km;
julia> Lon3D,Lat3D,Depth3D = LonLatDepthGrid(Lon, Lat, Depth);
julia> Lon3D
3×4×3 Array{Float64, 3}:
[:, :, 1] =
 11.0  11.0  11.0  11.0
 15.0  15.0  15.0  15.0
 19.0  19.0  19.0  19.0

[:, :, 2] =
 11.0  11.0  11.0  11.0
 15.0  15.0  15.0  15.0
 19.0  19.0  19.0  19.0

[:, :, 3] =
 11.0  11.0  11.0  11.0
 15.0  15.0  15.0  15.0
 19.0  19.0  19.0  19.0
julia> Lat3D
 3×4×3 Array{Float64, 3}:
 [:, :, 1] =
  1.0  4.0  7.0  10.0
  1.0  4.0  7.0  10.0
  1.0  4.0  7.0  10.0

 [:, :, 2] =
  1.0  4.0  7.0  10.0
  1.0  4.0  7.0  10.0
  1.0  4.0  7.0  10.0

 [:, :, 3] =
  1.0  4.0  7.0  10.0
  1.0  4.0  7.0  10.0
  1.0  4.0  7.0  10.0
julia> Depth3D
  3×4×3 Array{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}, 3}:
  [:, :, 1] =
   -20.0 km  -20.0 km  -20.0 km  -20.0 km
   -20.0 km  -20.0 km  -20.0 km  -20.0 km
   -20.0 km  -20.0 km  -20.0 km  -20.0 km

  [:, :, 2] =
   -15.0 km  -15.0 km  -15.0 km  -15.0 km
   -15.0 km  -15.0 km  -15.0 km  -15.0 km
   -15.0 km  -15.0 km  -15.0 km  -15.0 km

  [:, :, 3] =
   -10.0 km  -10.0 km  -10.0 km  -10.0 km
   -10.0 km  -10.0 km  -10.0 km  -10.0 km
   -10.0 km  -10.0 km  -10.0 km  -10.0 km
julia> Data        =   zeros(size(Lon3D));
julia> Data_set    =   GeophysicalModelGenerator.GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Data,))
GeoData
  size      : (3, 4, 3)
  lon       ϵ [ 11.0 : 19.0]
  lat       ϵ [ 1.0 : 10.0]
  depth     ϵ [ -20.0 km : -10.0 km]
  fields    : (:DataFieldName,)
  attributes: ["note"]
```
"""
struct GeoData
    lon     ::  GeoUnit
    lat     ::  GeoUnit
    depth   ::  GeoUnit
    fields  ::  NamedTuple
    atts    ::  Dict

    # Ensure that the data is of the correct format
    function GeoData(lon,lat,depth,fields,atts=nothing)

        # check depth & convert it to units of km in case no units are given or it has different length units
        if unit.(depth)[1]==NoUnits
            depth = depth*km                # in case depth has no dimensions
        end
        depth = uconvert.(km,depth)         # convert to km
        depth = GeoUnit(depth,km)           # convert to GeoUnit structure with units of km

        # Check ordering of the arrays in case of 3D
        if sum(size(lon).>1)==3
            if maximum(abs.(diff(lon,dims=2)))>maximum(abs.(diff(lon,dims=1))) || maximum(abs.(diff(lon,dims=3)))>maximum(abs.(diff(lon,dims=1)))
                error("It appears that the lon array has a wrong ordering")
            end
            if maximum(abs.(diff(lat,dims=1)))>maximum(abs.(diff(lat,dims=2))) || maximum(abs.(diff(lat,dims=3)))>maximum(abs.(diff(lat,dims=2)))
                error("It appears that the lat array has a wrong ordering")
            end
        end

        # fields should be a NamedTuple. In case we simply provide an array, lets transfer it accordingly
        if !(typeof(fields)<: NamedTuple)
            if (typeof(fields)<: Tuple)
                if length(fields)==1
                    fields = (DataSet1=first(fields),)  # The field is a tuple; create a NamedTuple from it
                else
                    error("Please employ a NamedTuple as input, rather than  a Tuple")  # out of luck
                end
            else
                fields = (DataSet1=fields,)
            end
        end

        DataField = fields[1];
        if typeof(DataField)<: Tuple
            DataField = DataField[1];           # in case we have velocity vectors as input
        end

        if !(size(lon)==size(lat)==size(depth)==size(DataField))
            error("The size of Lon/Lat/Depth and the Fields should all be the same!")
        end

        if isnothing(atts)
            # if nothing is given as attributes, then we note that in GeoData
            atts = Dict("note" => "No attributes were given to this dataset")
        else
            # check if a dict was given
            if !(typeof(atts)<: Dict)
                error("Attributes should be given as Dict!")
            end
        end

        return new(lon,lat,depth,fields,atts)

     end

end

# Print an overview of the Geodata struct:
function Base.show(io::IO, d::GeoData)
    println(io,"GeoData ")
    println(io,"  size      : $(size(d.lon))")
    println(io,"  lon       ϵ [ $(minimum(d.lon.val)) : $(maximum(d.lon.val))]")
    println(io,"  lat       ϵ [ $(minimum(d.lat.val)) : $(maximum(d.lat.val))]")
    println(io,"  depth     ϵ [ $(minimum(d.depth.val)) : $(maximum(d.depth.val))]")
    println(io,"  fields    : $(keys(d.fields))")
    if any( propertynames(d) .== :atts)
        println(io,"  attributes: $(keys(d.atts))")
    end
end

"""
    ParaviewData(x::GeoUnit, y::GeoUnit, z::GeoUnit, values::NamedTuple)

Cartesian data in `x/y/z` coordinates to be used with Paraview.
This is usually generated automatically from the `GeoData` structure, but you can also invoke do this manually:

```julia-repl
julia> Data_set    =   GeophysicalModelGenerator.GeoData(1.0:10.0,11.0:20.0,(-20:-11)*km,(DataFieldName=(-20:-11),))
julia> Data_cart = convert(ParaviewData, Data_set)
```
"""
mutable struct ParaviewData
    x       ::  GeoUnit
    y       ::  GeoUnit
    z       ::  GeoUnit
    fields  ::  NamedTuple
end

# Print an overview of the ParaviewData struct:
function Base.show(io::IO, d::ParaviewData)
    println(io,"ParaviewData ")
    println(io,"  size  : $(size(d.x))")
    println(io,"  x     ϵ [ $(minimum(d.x.val)) : $(maximum(d.x.val))]")
    println(io,"  y     ϵ [ $(minimum(d.y.val)) : $(maximum(d.y.val))]")
    println(io,"  z     ϵ [ $(minimum(d.z.val)) : $(maximum(d.z.val))]")
    println(io,"  fields: $(keys(d.fields))")
end

# conversion function from GeoData -> ParaviewData
function Base.convert(::Type{ParaviewData}, d::GeoData)

    # Utilize the Geodesy.jl package & use the Cartesian Earth-Centered-Earth-Fixed (ECEF) coordinate system
    lon         =   Array(ustrip.(d.lon.val));
    lat         =   Array(ustrip.(d.lat.val));
    LLA_Data    =   LLA.(lat,lon, Array(ustrip.(d.depth.val))*1000);            # convert to LLA from Geodesy package
    X,Y,Z       =   zeros(size(lon)), zeros(size(lon)), zeros(size(lon));

    # convert to cartesian ECEF reference frame. Note that we use kilometers and the wgs84
    for i in eachindex(X)
        data_xyz = ECEF(LLA_Data[i], wgs84)
        X[i] = data_xyz.x/1e3;
        Y[i] = data_xyz.y/1e3;
        Z[i] = data_xyz.z/1e3;
    end

    # This is the 'old' implementation, which does not employ a reference ellipsoid
    # X = R .* cosd.( lon ) .* cosd.( lat );
    # Y = R .* sind.( lon ) .* cosd.( lat );
    # Z = R .* sind.( lat );

    # In case any of the fields in the tuple has length 3, it is assumed to be a vector, so transfer it
    field_names = keys(d.fields)
    for i=1:length(d.fields)
        if typeof(d.fields[i]) <: Tuple
            if length(d.fields[i]) == 3
                # the tuple has length 3, which is therefore assumed to be a velocity vector

                # If the field name contains the string "color" we do not apply a vector transformation as it is supposed to contain RGB colors
                if !occursin("color", string(field_names[i]))
                    println("Applying a vector transformation to field: $(field_names[i])")
                    Velocity_SphericalToCartesian!(d, d.fields[i])  # Transfer it to x/y/z format
                end
            end
        end
    end

    return ParaviewData(GeoUnit(X,km),GeoUnit(Y,km),GeoUnit(Z,km),d.fields)
end

"""
    UTMData(EW::Any, NS:Any, depth::GeoUnit, UTMZone::Int, NorthernHemisphere=true, fields::NamedTuple)

Data structure that holds one or several fields with UTM coordinates (east-west), (north-south) and depth information.

- `depth` can have units of meters, kilometer or be unitless; it will be converted to meters (as UTMZ is usually in meters)
- `fields` should ideally be a NamedTuple which allows you to specify the names of each of the fields.
- In case you only pass one array we will convert it to a NamedTuple with default name.
- A single field should be added as `(DataFieldName=Data,)` (don't forget the comma at the end).
- Multiple fields  can be added as well.
- In case you want to display a vector field in paraview, add it as a tuple: `(Velocity=(Veast,Vnorth,Vup), Veast=Veast, Vnorth=Vnorth, Vup=Vup)`; we automatically apply a vector transformation when transforming this to a `ParaviewData` structure from which we generate Paraview output. As this changes the magnitude of the arrows, you will no longer see the `[Veast,Vnorth,Vup]` components in Paraview which is why it is a good ideas to store them as separate Fields.
- Yet, there is one exception: if the name of the 3-component field is `colors`, we do not apply this vector transformation as this field is regarded to contain RGB colors.
- `Lat`,`Lon`,`Depth` should have the same size as the `Data` array. The ordering of the arrays is important. If they are 3D arrays, as in the example below, we assume that the first dimension corresponds to `lon`, second dimension to `lat` and third dimension to `depth` (which should be in km). See below for an example.

# Example
```julia-repl
julia> ew          =   422123.0:100:433623.0
julia> ns          =   4.514137e6:100:4.523637e6
julia> depth       =   -5400:250:600
julia> EW,NS,Depth =   XYZGrid(ew, ns, depth);
julia> Data        =   ustrip.(Depth);
julia> Data_set    =   UTMData(EW,NS,Depth,33, true, (FakeData=Data,Data2=Data.+1.))
UTMData
  UTM zone : 33-33 North
    size    : (116, 96, 25)
    EW      ϵ [ 422123.0 : 433623.0]
    NS      ϵ [ 4.514137e6 : 4.523637e6]
    depth   ϵ [ -5400.0 m : 600.0 m]
    fields  : (:FakeData, :Data2)
  attributes: ["note"]
```
If you wish, you can convert this from `UTMData` to `GeoData` with
```julia-repl
julia> Data_set1 =  convert(GeoData, Data_set)
GeoData
  size      : (116, 96, 25)
  lon       ϵ [ 14.075969111533457 : 14.213417764154963]
  lat       ϵ [ 40.77452227533946 : 40.86110443583479]
  depth     ϵ [ -5.4 km : 0.6 km]
  fields    : (:FakeData, :Data2)
  attributes: ["note"]
```
which would allow visualizing this in paraview in the usual manner:
```julia-repl
julia> Write_Paraview(Data_set1, "Data_set1")
1-element Vector{String}:
 "Data_set1.vts"
```
"""
struct UTMData
    EW       ::  GeoUnit
    NS       ::  GeoUnit
    depth    ::  GeoUnit
    zone     ::  Any
    northern ::  Any
    fields   ::  NamedTuple
    atts     ::  Dict

    # Ensure that the data is of the correct format
    function UTMData(EW,NS,depth,zone,northern,fields,atts=nothing)

        # check depth & convert it to units of km in case no units are given or it has different length units
        if unit.(depth)[1]==NoUnits
            depth = depth*m                # in case depth has no dimensions
        end
        depth = uconvert.(m,depth)         # convert to meters
        depth = GeoUnit(depth,m)           # convert to GeoUnit structure with units of meters

        # Check ordering of the arrays in case of 3D
        if sum(size(EW).>1)==3
            if maximum(abs.(diff(EW,dims=2)))>maximum(abs.(diff(EW,dims=1))) || maximum(abs.(diff(EW,dims=3)))>maximum(abs.(diff(EW,dims=1)))
                error("It appears that the EW array has a wrong ordering")
            end
            if maximum(abs.(diff(NS,dims=1)))>maximum(abs.(diff(NS,dims=2))) || maximum(abs.(diff(NS,dims=3)))>maximum(abs.(diff(NS,dims=2)))
                error("It appears that the NS array has a wrong ordering")
            end
        end

        # fields should be a NamedTuple. In case we simply provide an array, lets transfer it accordingly
        if !(typeof(fields)<: NamedTuple)
            if (typeof(fields)<: Tuple)
                if length(fields)==1
                    fields = (DataSet1=first(fields),)  # The field is a tuple; create a NamedTuple from it
                else
                    error("Please employ a NamedTuple as input, rather than  a Tuple")  # out of luck
                end
            else
                fields = (DataSet1=fields,)
            end
        end

        DataField = fields[1];
        if typeof(DataField)<: Tuple
            DataField = DataField[1];           # in case we have velocity vectors as input
        end

        if !(size(EW)==size(NS)==size(depth)==size(DataField))
            error("The size of EW/NS/Depth and the Fields should all be the same!")
        end

        if length(zone)==1
            zone = ones(Int64,size(EW))*zone
            northern = ones(Bool,size(EW))*northern
        end

        # take care of attributes
        if isnothing(atts)
            # if nothing is given as attributes, then we note that in GeoData
            atts = Dict("note" => "No attributes were given to this dataset")
        else
            # check if a dict was given
            if !(typeof(atts)<: Dict)
                error("Attributes should be given as Dict!")
            end
        end

        return new(EW,NS,depth,zone,northern, fields,atts)

     end

end

# Print an overview of the UTMData struct:
function Base.show(io::IO, d::UTMData)
    println(io,"UTMData ")
    if d.northern[1]
        println(io,"  UTM zone : $(minimum(d.zone))-$(maximum(d.zone)) North")
    else
        println(io,"  UTM zone : $(minimum(d.zone))-$(maximum(d.zone)) South")
    end
    println(io,"    size    : $(size(d.EW))")
    println(io,"    EW      ϵ [ $(minimum(d.EW.val)) : $(maximum(d.EW.val))]")
    println(io,"    NS      ϵ [ $(minimum(d.NS.val)) : $(maximum(d.NS.val))]")
    println(io,"    depth   ϵ [ $(minimum(d.depth.val)) : $(maximum(d.depth.val))]")
    println(io,"    fields  : $(keys(d.fields))")
    if any( propertynames(d) .== :atts)
        println(io,"  attributes: $(keys(d.atts))")
    end
end

"""
Converts a `UTMData` structure to a `GeoData` structure
"""
function Base.convert(::Type{GeoData}, d::UTMData)

    Lat = zeros(size(d.EW));
    Lon = zeros(size(d.EW));
    for i in eachindex(d.EW.val)

        # Use functions of the Geodesy package to convert to LLA
        utmz_i  = UTMZ(d.EW.val[i],d.NS.val[i],Float64(ustrip.(d.depth.val[i])),d.zone[i],d.northern[i])
        lla_i   = LLA(utmz_i,wgs84)
        lon = lla_i.lon;
       # if lon<0; lon = 360+lon; end # as GMT expects this

        Lat[i] = lla_i.lat
        Lon[i] = lon
    end

    # handle the case where an old GeoData structure is converted
    if any( propertynames(d) .== :atts)
        atts = d.atts;
    else
        atts = Dict("note" => "No attributes were given to this dataset") # assign the default
    end

    return GeoData(Lon,Lat,d.depth.val,d.fields,atts)

end

"""
Converts a `GeoData` structure to a `UTMData` structure
"""
function Base.convert(::Type{UTMData}, d::GeoData)

    EW = zeros(size(d.lon));
    NS = zeros(size(d.lon));
    depth = zeros(size(d.lon));
    zone = zeros(Int64,size(d.lon));
    northern = zeros(Bool,size(d.lon));
    for i in eachindex(d.lon.val)

        # Use functions of the Geodesy package to convert to LLA
        lla_i   = LLA(d.lat.val[i],d.lon.val[i],Float64(ustrip.(d.depth.val[i])*1e3))
        utmz_i  = UTMZ(lla_i, wgs84)

        EW[i] = utmz_i.x
        NS[i] = utmz_i.y
        depth[i] = utmz_i.z
        zone[i] = utmz_i.zone;
        northern[i] = utmz_i.isnorth
    end

    # handle the case where an old GeoData structure is converted
    if any( propertynames(d) .== :atts)
        atts = d.atts;
    else
        atts = Dict("note" => "No attributes were given to this dataset") # assign the default
    end

    return UTMData(EW,NS,depth,zone, northern, d.fields, atts)

end

"""
    Convert2UTMzone(d::GeoData, p::ProjectionPoint)

Converts a `GeoData` structure to fixed UTM zone, around a given `ProjectionPoint`
    This useful to use real data as input for a cartesian geodynamic model setup (such as in LaMEM). In that case, we need to project map coordinates to cartesian coordinates.
    One way to do this is by using UTM coordinates. Close to the `ProjectionPoint` the resulting coordinates will be rectilinear and distance in meters. The map distortion becomes larger the further you are away from the center.

"""
function Convert2UTMzone(d::GeoData, proj::ProjectionPoint)

    EW = zeros(size(d.lon));
    NS  = zeros(size(d.lon));
    zone        = zeros(Int64,size(d.lon));
    northern    = zeros(Bool,size(d.lon));
    trans       = UTMfromLLA(proj.zone, proj.isnorth, wgs84)
    for i in eachindex(d.lon.val)

        # Use functions of the Geodesy package to convert to LLA
        lla_i  =   LLA(d.lat.val[i],d.lon.val[i],Float64(ustrip.(d.depth.val[i])*1e3))
        utm_i  =   trans(lla_i)

        EW[i] = utm_i.x
        NS[i] = utm_i.y
        zone[i] = proj.zone;
        northern[i] = proj.isnorth
    end

    # handle the case where an old GeoData structure is converted
    if any( propertynames(d) .== :atts)
        atts = d.atts;
    else
        atts = Dict("note" => "No attributes were given to this dataset") # assign the default
    end

    return UTMData(EW,NS,d.depth.val,zone, northern, d.fields,atts)

end

"""
    CartData(x::Any, y::Any, z::GeoUnit, fields::NamedTuple)

Data structure that holds one or several fields with with Cartesian x/y/z coordinates. Distances are in kilometers

- `x`,`y`,`z` can have units of meters, kilometer or be unitless; they will be converted to kilometers
- `fields` should ideally be a NamedTuple which allows you to specify the names of each of the fields.
- In case you only pass one array we will convert it to a NamedTuple with default name.
- A single field should be added as `(DataFieldName=Data,)` (don't forget the comma at the end).
- Multiple fields  can be added as well.
- In case you want to display a vector field in paraview, add it as a tuple: `(Velocity=(Vx,Vnorth,Vup), Veast=Veast, Vnorth=Vnorth, Vup=Vup)`; we automatically apply a vector transformation when transforming this to a `ParaviewData` structure from which we generate Paraview output. As this changes the magnitude of the arrows, you will no longer see the `[Veast,Vnorth,Vup]` components in Paraview which is why it is a good ideas to store them as separate Fields.
- Yet, there is one exception: if the name of the 3-component field is `colors`, we do not apply this vector transformation as this field is regarded to contain RGB colors.
- `x`,`y`,`z` should have the same size as the `Data` array. The ordering of the arrays is important. If they are 3D arrays, as in the example below, we assume that the first dimension corresponds to `x`, second dimension to `y` and third dimension to `z` (which should be in km). See below for an example.

# Example
```julia-repl
julia> x        =   0:2:10
julia> y        =   -5:5
julia> z        =   -10:2:2
julia> X,Y,Z    =   XYZGrid(x, y, z);
julia> Data     =   Z
julia> Data_set =   CartData(X,Y,Z, (FakeData=Data,Data2=Data.+1.))
CartData
    size    : (6, 11, 7)
    x       ϵ [ 0.0 km : 10.0 km]
    y       ϵ [ -5.0 km : 5.0 km]
    z       ϵ [ -10.0 km : 2.0 km]
    fields  : (:FakeData, :Data2)
  attributes: ["note"]
```
`CartData` is particularly useful in combination with cartesian geodynamic codes, such as LaMEM, which require cartesian grids.
You can directly save your data to Paraview with
```julia-repl
julia> Write_Paraview(Data_set, "Data_set")
1-element Vector{String}:
 "Data_set.vts"
```

If you wish, you can convert this to `UTMData` (which will simply convert the )
```julia-repl
julia> Data_set1 =  convert(GeoData, Data_set)
GeoData
  size  : (116, 96, 25)
  lon   ϵ [ 14.075969111533457 : 14.213417764154963]
  lat   ϵ [ 40.77452227533946 : 40.86110443583479]
  depth ϵ [ -5.4 km : 0.6 km]
  fields: (:FakeData, :Data2)
```
which would allow visualizing this in paraview in the usual manner:

"""
struct CartData
    x       ::  GeoUnit
    y       ::  GeoUnit
    z       ::  GeoUnit
    fields  ::  NamedTuple
    atts    ::  Dict

    # Ensure that the data is of the correct format
    function CartData(x,y,z,fields,atts=nothing)

        # Check ordering of the arrays in case of 3D
        if sum(size(x).>1)==3
            if maximum(abs.(diff(x,dims=2)))>maximum(abs.(diff(x,dims=1))) || maximum(abs.(diff(x,dims=3)))>maximum(abs.(diff(x,dims=1)))
                error("It appears that the x-array has a wrong ordering")
            end
            if maximum(abs.(diff(y,dims=1)))>maximum(abs.(diff(y,dims=2))) || maximum(abs.(diff(y,dims=3)))>maximum(abs.(diff(y,dims=2)))
                error("It appears that the y-array has a wrong ordering")
            end
        end

        # check depth & convert it to units of km in case no units are given or it has different length units
        x = Convert!(x,km)
        y = Convert!(y,km)
        z = Convert!(z,km)

        # fields should be a NamedTuple. In case we simply provide an array, lets transfer it accordingly
        if !(typeof(fields)<: NamedTuple)
            if (typeof(fields)<: Tuple)
                if length(fields)==1
                    fields = (DataSet1=first(fields),)  # The field is a tuple; create a NamedTuple from it
                else
                    error("Please employ a NamedTuple as input, rather than a Tuple")  # out of luck
                end
            else
                fields = (DataSet1=fields,)
            end
        end

        DataField = fields[1];
        if typeof(DataField)<: Tuple
            DataField = DataField[1];           # in case we have velocity vectors as input
        end

        if !(size(x)==size(y)==size(z)==size(DataField))
            error("The size of x/y/z and the Fields should all be the same!")
        end

        # take care of attributes
        if isnothing(atts)
            # if nothing is given as attributes, then we note that
            atts = Dict("note" => "No attributes were given to this dataset")
        else
            # check if a dict was given
            if !(typeof(atts)<: Dict)
                error("Attributes should be given as Dict!")
            end
        end

        return new(x,y,z,fields,atts)

     end

end

# Print an overview of the UTMData struct:
function Base.show(io::IO, d::CartData)
    println(io,"CartData ")
    println(io,"    size    : $(size(d.x))")
    println(io,"    x       ϵ [ $(minimum(d.x.val)) : $(maximum(d.x.val))]")
    println(io,"    y       ϵ [ $(minimum(d.y.val)) : $(maximum(d.y.val))]")
    println(io,"    z       ϵ [ $(minimum(d.z.val)) : $(maximum(d.z.val))]")
    println(io,"    fields  : $(keys(d.fields))")
    if any( propertynames(d) .== :atts)
        println(io,"  attributes: $(keys(d.atts))")
    end
end

"""
    CartData(xyz::Tuple{Array,Array,Array})

This creates a `CartData` struct if you have a Tuple with 3D coordinates as input.
# Example
```julia
julia> data = CartData(XYZGrid(-10:10,-5:5,0))
CartData
    size    : (21, 11, 1)
    x       ϵ [ -10.0 km : 10.0 km]
    y       ϵ [ -5.0 km : 5.0 km]
    z       ϵ [ 0.0 km : 0.0 km]
    fields  : (:Z,)
  attributes: ["note"]
```
"""
CartData(xyz::Tuple) = CartData(xyz[1],xyz[2],xyz[3],(Z=xyz[3],))

"""
    Convert2UTMzone(d::CartData, proj::ProjectionPoint)

This transfers a `CartData` dataset to a `UTMData` dataset, that has a single UTM zone. The point around which we project is `ProjectionPoint`
"""
function Convert2UTMzone(d::CartData, proj::ProjectionPoint)

    return UTMData(ustrip.(d.x.val).*1e3 .+ proj.EW,ustrip.(d.y.val).*1e3 .+ proj.NS,
                   ustrip.(d.z.val).*1e3,proj.zone, proj.isnorth, d.fields, d.atts)

end

"""
    Convert2CartData(d::UTMData, proj::ProjectionPoint)
Converts a `UTMData` structure to a `CartData` structure, which essentially transfers the dimensions to km
"""
function Convert2CartData(d::UTMData, proj::ProjectionPoint)

        # handle the case where an old structure is converted
        if any( propertynames(d) .== :atts)
            atts = d.atts;
        else
            atts = Dict("note" => "No attributes were given to this dataset") # assign the default
        end

    return CartData( (ustrip.(d.EW.val) .- proj.EW)./1e3, (ustrip.(d.NS.val) .- proj.NS)./1e3,
                     ustrip.(d.depth.val)./1e3, d.fields,atts)
end

"""
    Convert2CartData(d::GeoData, proj::ProjectionPoint)
Converts a `GeoData` structure to a `CartData` structure, which essentially transfers the dimensions to km
"""
function Convert2CartData(d::GeoData, proj::ProjectionPoint)

    d_UTM = Convert2UTMzone(d,proj)
    return CartData( (ustrip.(d_UTM.EW.val) .- proj.EW)./1e3, (ustrip.(d_UTM.NS.val) .- proj.NS)./1e3,
                     ustrip.(d_UTM.depth.val)./1e3, d_UTM.fields,d_UTM.atts)
end

"""
    Lon, Lat, Depth = LonLatDepthGrid(Lon::Any, Lat::Any, Depth:Any)

Creates 3D arrays of `Lon`, `Lat`, `Depth` from 1D vectors or numbers

# Example 1: Create 3D grid
```julia-repl
julia> Lon,Lat,Depth =  LonLatDepthGrid(10:20,30:40,(-10:-1)km);
julia> size(Lon)
(11, 11, 10)
```

# Example 2: Create 2D lon/lat grid @ a given depth
```julia-repl
julia> Lon,Lat,Depth =  LonLatDepthGrid(10:20,30:40,-50km);
julia> size(Lon)
(11, 11)
```

# Example 3: Create 2D lon/depth grid @ a given lat
```julia-repl
julia> Lon,Lat,Depth =  LonLatDepthGrid(10:20,30,(-10:-1)km);
julia> size(Lon)
(11, 11)
```
# Example 4: Create 1D vertical line @ a given lon/lat point
```julia-repl
julia> Lon,Lat,Depth =  LonLatDepthGrid(10,30,(-10:-1)km);
julia> size(Lon)
(10, )
```

"""
function LonLatDepthGrid(Lon::Any, Lat::Any, Depth::Any)

    nLon    = length(Lon)
    nLat    = length(Lat)
    nDepth  = length(Depth)

    if nLon==nLat==nDepth==1
        error("Cannot use this routine for a 3D point (no need to create a grid in that case")
    end
    if maximum([length(size(Lon)), length(size(Lat)), length(size(Depth))])>1
        error("You can only give 1D vectors or numbers as input")
    end

    Lon3D   =   zeros(nLon,nLat,nDepth);
    Lat3D   =   zeros(nLon,nLat,nDepth);
    Depth3D =   zeros(nLon,nLat,nDepth);

    for i=1:nLon
        for j=1:nLat
            for k=1:nDepth
                Lon3D[i,j,k]    =   ustrip.(Lon[i]);
                Lat3D[i,j,k]    =   ustrip.(Lat[j]);
                Depth3D[i,j,k]  =   ustrip.(Depth[k]);
            end
        end
    end

    # Add dimensions back
    Lon3D   =   Lon3D*unit(  Lon[1])
    Lat3D   =   Lat3D*unit(  Lat[1])
    Depth3D = Depth3D*unit(Depth[1])

    return Lon3D, Lat3D, Depth3D
end

"""
    X,Y,Z = XYZGrid(X_vec::Any, Y_vec::Any, Z_vec::Any)

Creates a `X,Y,Z` grid. It works just as `LonLatDepthGrid` apart from the better suited name.

# Example 1: Create 3D grid
```julia-repl
julia> X,Y,Z =  XYZGrid(10:20,30:40,(-10:-1)km);
julia> size(X)
(11, 11, 10)
```

See `LonLatDepthGrid` for more examples.
"""
function XYZGrid(X_vec::Any, Y_vec::Any, Z_vec::Any)
    return X,Y,Z = LonLatDepthGrid(X_vec,Y_vec,Z_vec)
end

"""
    Velocity_SphericalToCartesian!(Data::GeoData, Velocity::Tuple)

In-place conversion of velocities in spherical velocities `[Veast, Vnorth, Vup]` to cartesian coordinates (for use in paraview).

NOTE: the magnitude of the vector will be the same, but the individual `[Veast, Vnorth, Vup]` components
will not be retained correctly (as a different `[x,y,z]` coordinate system is used in paraview).
Therefore, if you want to display or color that correctly in Paraview, you need to store these magnitudes as separate fields

"""
function Velocity_SphericalToCartesian!(Data::GeoData, Velocity::Tuple)
    # Note: This is partly based on scripts originally written by Tobias Baumann, Uni Mainz

    for i in eachindex(Data.lat.val)
        az  =   Data.lon.val[i];
        el  =   Data.lat.val[i];

        R           = [-sind(az) -sind(el)*cosd(az) cosd(el)*cosd(az);
                        cosd(az) -sind(el)*sind(az) cosd(el)*sind(az);
                        0.0       cosd(el)          sind(el)            ];

        V_sph       =   [Velocity[1][i]; Velocity[2][i]; Velocity[3][i] ];

        # Normalize spherical velocity
        V_mag       =  sum(sqrt.(V_sph.^2));        # magnitude
        V_norm      =  V_sph/V_mag

        V_xyz_norm  =  R*V_norm;
        V_xyz       =  V_xyz_norm.*V_mag;          # scale with magnitude

        # in-place saving of rotated velocity
        Velocity[1][i] = V_xyz[1];
        Velocity[2][i] = V_xyz[2];
        Velocity[3][i] = V_xyz[3];
    end
end

# Internal function that converts arrays to a GeoUnit with certain units
function Convert!(d,u)
    if unit.(d)[1]==NoUnits
        d = d*u                # in case it has no dimensions
    end
    d = uconvert.(u,d)         # convert to u
    d = GeoUnit(d,u)           # convert to GeoUnit structure with units of u

    return d
end

ext: jl


text:
# Greedy Coordinate Descent (GreedyCD) algorithm
#
# Author: Takehiro Sano
#
# Reference: Cho-Jui Hsieh and Inderjit S. Dhillon. "Fast coordinate descent methods
#  with variable selection for non-negative matrix factorization." In Proceedings of
#  the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining,
#  pp. 1064–1072, 2011.

mutable struct GreedyCD{T}
    maxiter::Int           # maximum number of iterations (in main procedure)
    verbose::Bool          # whether to show procedural information
    tol::T                 # tolerance of changes on W and H upon convergence
    update_H::Bool         # whether to update H
    lambda_w::T            # L1 regularization coefficient for W
    lambda_h::T            # L1 regularization coefficient for H

    function GreedyCD{T}(;maxiter::Integer=100,
                          verbose::Bool=false,
                          tol::Real=cbrt(eps(T)),
                          update_H::Bool=true,
                          lambda_w::Real=zero(T),
                          lambda_h::Real=zero(T)) where T

        maxiter > 1 || throw(ArgumentError("maxiter must be greater than 1."))
        tol > 0 || throw(ArgumentError("tol must be positive."))
        lambda_w >= 0 || throw(ArgumentError("lambda_w must be non-negative."))
        lambda_h >= 0 || throw(ArgumentError("lambda_h must be non-negative."))
        new{T}(maxiter, verbose, tol, update_H, lambda_w, lambda_h)
    end
end

solve!(alg::GreedyCD{T}, X, W, H) where T =
    nmf_skeleton!(GreedyCDUpd{T}(alg.update_H, alg.lambda_w, alg.lambda_h), X, W, H, alg.maxiter, alg.verbose, alg.tol)

struct GreedyCDUpd{T} <: NMFUpdater{T}
    update_H::Bool
    lambda_w::T
    lambda_h::T
end

struct GreedyCDUpd_State{T}
    WH::Matrix{T}
    PWW::Matrix{T}
    PHH::Matrix{T}
    PXW::Matrix{T}
    PXH::Matrix{T}
    GW::Matrix{T}
    GH::Matrix{T}
    Wnew::Matrix{T}
    Hnew::Matrix{T}
    SW::Matrix{T}
    SH::Matrix{T}
    DW::Matrix{T}
    DH::Matrix{T}
    qW::Vector{Int}
    qH::Vector{Int}

    function GreedyCDUpd_State{T}(X, W::Matrix{T}, H::Matrix{T}) where T
        p, n, k = nmf_checksize(X, W, H)
        new{T}(Matrix{T}(undef, p, n),
               Matrix{T}(undef, k, k),
               Matrix{T}(undef, k, k),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, n, k),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, n, k),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, n, k),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, n, k),
               Matrix{T}(undef, p, k),
               Matrix{T}(undef, n, k),
               Vector{Int}(undef, p),
               Vector{Int}(undef, n))
    end
end

prepare_state(::GreedyCDUpd{T}, X, W, H) where T = GreedyCDUpd_State{T}(X, W, H)

function evaluate_objv(u::GreedyCDUpd{T}, s::GreedyCDUpd_State{T}, X, W, H) where T
    mul!(s.WH, W, H)
    r = convert(T, 0.5) * sqL2dist(X, s.WH)
    if u.lambda_w > 0
        r += u.lambda_w * norm(W, 1)
    end
    if u.lambda_h > 0
        r += u.lambda_h * norm(H, 1)
    end
    return r
end

function _update_GreedyCD!(upd::GreedyCDUpd{T}, s::GreedyCDUpd_State{T}, X,
                           W::AbstractArray{T}, Ht::AbstractArray{T}, W_flag::Bool) where T
    if W_flag
        lambda = upd.lambda_w
        Wnew = s.Wnew
        P = s.PWW
        Z = s.PXW
        G = s.GW
        S = s.SW
        D = s.DW
        q = s.qW
    else
        lambda = upd.lambda_h
        Wnew = s.Hnew
        P = s.PHH
        Z = s.PXH
        G = s.GH
        S = s.SH
        D = s.DH
        q = s.qH
    end
    n_samples, n_components = size(W)

    mul!(P, transpose(Ht), Ht)
    mul!(Z, X, Ht)
    mul!(G, W, P)
    G .-= Z
    if lambda > zero(T)
        G .+= lambda
    end

    for r in 1:n_components
        for i in 1:n_samples
            S[i, r] = max(zero(T), W[i, r] - G[i, r] / (eps(T) + P[r, r])) - W[i, r]
            D[i, r] = - G[i, r] * S[i, r] - convert(T, 0.5) * P[r, r] * S[i, r]^2
        end
    end

    p_init = convert(T, -1.0)
    for i in 1:n_samples
        qi = argmax(view(D, i, :))
        q[i] = qi
        p_init = max(p_init, D[i, qi])
    end

    fill!(Wnew, zero(T))
    nu = convert(T, 0.001)

    for i in 1:n_samples
        qi = q[i]
        for _ in 1:n_components^2
            if D[i, qi] < nu * p_init
                break
            end

            Wnew[i, qi] += S[i, qi]

            for r in 1:n_components
                G[i, r] += S[i, qi] * P[qi, r]
            end

            for r in 1:n_components
                S[i, r] = max(zero(T), W[i, r] - G[i, r] / (eps(T) + P[r, r])) - W[i, r]
                D[i, r] = - G[i, r] * S[i, r] - convert(T, 0.5) * P[r, r] * S[i, r]^2
            end

            qi = argmax(view(D, i, :))
        end
        q[i] = qi
    end
    copyto!(W, W + Wnew)
    projectnn!(W) # countermeasure on rounding error
end

function update_wh!(upd::GreedyCDUpd{T}, s::GreedyCDUpd_State{T}, X, W::Matrix{T}, H::Matrix{T}) where T
    # update W
    Ht = transpose(H)
    _update_GreedyCD!(upd, s, X, W, Ht, true)

    # update H
    if upd.update_H
        Xt = transpose(X)
        _update_GreedyCD!(upd, s, Xt, Ht, W, false)
    end
end


text:
##############################################################################
##
## Combination behaves like [A B C ...] without forming it
##
##############################################################################

struct Combination{N}
    A::NTuple{N, Matrix{Float64}}
    cumlength::Vector{Int}
end

function Combination(A::Matrix{Float64}...)
    Combination(A, cumsum([size(x, 2) for x in A]))
end

function size(c::Combination, i)
    if i == 1
        size(c.A[1], 1)
    elseif i == 2
        c.cumlength[end]
    end
end

function view(c::Combination, ::Colon, j)
    index = searchsortedfirst(c.cumlength, j)
    newj = j
    if index > 1
        newj = j - c.cumlength[index-1]
    end
    view(c.A[index], :, newj)
end

##############################################################################
##
## Crossprod computes [A B C ...]' [A B C ...] without forming it
##
##############################################################################

# Construct [A B C]'[A B C] without generating [A B C]
function crossprod(c::Combination{N}) where {N}
    out = Array{Float64}(undef, size(c, 2), size(c, 2))
    for j in 1:size(c, 2)
        viewj = view(c, :, j)
        for i in j:size(c, 2)
            out[i, j] = dot(viewj, view(c, :, i))
        end
    end
    # make symmetric
    for j in 1:size(c, 2), i in 1:(j-1)
        out[i, j] = out[j, i]
    end
    return out
end
crossprod(A::Matrix{Float64}) = A'A
crossprod(A::Matrix{Float64}...) = crossprod(Combination(A...))

##############################################################################
##
## Returns base of [A B C ...]
##
## TODO: You could protect against roundoff error by using a controlled sum algorithm (similar to sum_kbn) to compute elements of X'X, then converting to BigFloat before factoring.
##
##
##############################################################################

# rank(A) == rank(A'A)
function basecol(X::Matrix{Float64}...; factorization = :Cholesky)
    cholm = cholesky!(Symmetric(crossprod(X...)), Val(true); tol = -1, check = false)
    r = 0
    if size(cholm, 1) > 0
        r = sum(diag(cholm.factors) .> size(X[1],1)^2 * eps())
        # used to be r = rank(cholm) but does not work wiht very high regressors at the same time as intercept
    end
    invpermute!(1:size(cholm, 1) .<= r, cholm.piv)
end

function getcols(X::Matrix{Float64},  basecolX::BitArray{1})
    sum(basecolX) == size(X, 2) ? X : X[:, basecolX]
end


text:
# Recurrent inference machines for seismic imaging
# Same notation as Putzky and Welling paper (2017)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using JUDI.TimeModeling, LinearAlgebra, InvertibleNetworks
using PyPlot, Random, JLD, Flux

Random.seed!(11)

# Training data
D = load("/data/pwitte3/models/overthrust_images_test.jld")
ntest = length(D["m"])

# Crop images and models to 400 x 120 (must be evenly dividable by 4)
for j=1:ntest
    D["m"][j] = D["m"][j][1:400, 1:120]
    D["m0"][j] = D["m0"][j][1:400, 1:120]
    D["dm"][j] = D["dm"][j][1:400, 1:120]
end

# Use one sample to set up operators
m = D["m"][1]
m0 = D["m0"][1]
dm = vec(D["dm"][1])

# Set up model structure
n = size(m0)
d = (25., 25.)
o = (0., 0.)

# Setup info and model structure
nsrc = 1
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 301
xrec = range(400f0, stop=9600f0, length=nxrec)
yrec = 0f0
zrec = range(250f0, stop=250f0, length=nxrec)

# receiver sampling and recording time
time = 2000f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
num_simsource = 21
xsrc = range(500f0, stop=9500f0, length=num_simsource)
ysrc = range(0f0, stop=0f0, length=num_simsource)
zsrc = range(20f0, stop=20f0, length=num_simsource)

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=time)

# setup wavelet
f0 = 0.015f0     # MHz
wavelet = zeros(Float32, srcGeometry.nt[1], num_simsource)
for j=1:num_simsource
    wavelet[:, j] = ricker_wavelet(time, dt, f0) * randn(1)[1]/sqrt(num_simsource)
end
q = judiVector(srcGeometry, wavelet)
q_sim = zeros(Float32, size(wavelet))

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model0)
info = Info(prod(n), nsrc, ntComp)

####################################################################################################

# Return data as julia array
opt = Options(return_array=true)

# Setup operators
F0 = judiModeling(info, model0; options=opt)
Pr = judiProjection(info, recGeometry)
Ps = judiProjection(info, srcGeometry)

# Modeling operator and Jacobian
F0 = Pr*F0*Ps'  # initial model
J = judiJacobian(F0, q)

####################################################################################################

# Dimensions
nx1 = n[1]
nx2 = n[2]
nx_in = 32
nx_hidden = 64
batchsize = 1

#ny1 = recGeometry.nt[1]
#ny2 = nxrec
#ny_in = 1
#ny_hidden = 64

maxiter = 4
Ψ(η) = identity(η)

SLIM = Array{AffineCouplingLayerSLIM}(undef, maxiter)
Params = Array{Parameter}(undef, 0)

# i-SLIM coupling layer
for j=1:maxiter
    SLIM[j] = AffineCouplingLayerSLIM(nx1, nx2, nx_in, nx_hidden, batchsize, Ψ; logdet=false, permute=true)
    #SLIM[j] = ConditionalLayerSLIM(nx1, nx2, nx_in, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type="affine")
    global Params = cat(Params, get_params(SLIM[j]); dims=1)
end

function forward(X, Y, J)
    for j=1:maxiter
        X = SLIM[j].forward(X, Y, J)
    end
    return X
end

# Load params
iter = 600
P_save = load(join(["../data/network_islim_params_iteration_", string(iter), ".jld"]))["P"]
for j=1:length(Params)
    Params[j].data = P_save[j].data
end

####################################################################################################
# Evaluate trained network for new data

# Step 1: Draw new m0, q_sim and d
i = 99 #randperm(ntest)[1]
x = D["dm"][i]
m0 = D["m0"][i]
for k=1:num_simsource
    q_sim[:,k] = ricker_wavelet(time, dt, f0) * randn(Float32, 1)[1]/sqrt(1f0*num_simsource)
end

# Set up Jacobian
J.model.m = m0
J.source[1] = q_sim
Y = J*vec(x)    # Observed data

# Step 2: Compute predicted image
X0 = zeros(Float32, nx1, nx2, nx_in, batchsize)
X = forward(X0, Y, J)
x_ = X[:,:,1,1]

# Comparison to RTM
rtm = J'*vec(Y)

# Plot
x_ = x_ / norm(x_, 2)
x = x / norm(x, 2)
rtm = rtm / norm(rtm, 2)

figure(figsize=(5, 7))
subplot(3,1,1); imshow(reshape(x, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("True image")
subplot(3,1,2); imshow(reshape(rtm, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("RTM")
subplot(3,1,3); imshow(reshape(x_, model0.n)', cmap="gray", vmin=-2e-2, vmax=2e-2); title("i-SLIM (affine)")

max_stars_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_issues_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_forks_repo_name: slimgroup/LoopUnrolledSeismicImaging


text:
# just trying an example

push!(LOAD_PATH, "/home/renee/phd/src/charge-distribution.project/ebit-ode-server/")

import EbitSolver


text:
using JLD2
using FileIO
using Glob
using Distances
using Statistics

files = glob("out/img_*.jld2")
ids = map(files) do f
    parse(Int, match(r"img_(\d*)\.jld", f).captures[1])
end

templates = load("out/template.jld2", "template_features")
alldata = load.(files, "res")

x = pairwise(CosineDist(), templates, templates, dims = 2)

T = 0.6

y = filter(!isnothing, map(zip(ids, alldata)) do (id, (bboxes, landmarks, features))
    likely = vec(minimum(pairwise(CosineDist(), features', templates, dims = 2), dims = 2))
    n = argmin(likely)
    likely[n] > T && return nothing
    (id = id, p = likely[n], bbox = bboxes[n, :], landmark = landmarks[n,:], feature = features[n,:])
end)

@save "t.jld2" y

max_stars_repo_name: BillKerman/FaceNetCustomized

max_issues_repo_name: BillKerman/FaceNetCustomized

max_forks_repo_name: BillKerman/FaceNetCustomized


text:
using ConstantArrays
using Test: @test, @testset, @test_throws, @inferred

@testset "ConstantArrays" begin
    dim = (3,4)
    val = 7
    x = @inferred ConstantArray(val, dim)

    show(isinteractive() ? stdout : devnull, MIME("text/plain"), x)

    # getindex:
    @test x[1] == val
    @test x[end] == val
    @test x[[2,3]] == [val,val]
    @test x[:,1] == fill(val, dim[1])
    @test x[:] == fill(val, prod(dim))

    @test all(x .== val)
    @test size(x) == dim
    @test size(x,1) == dim[1]
    @test size(x,2) == dim[2]
    @test eltype(x) == eltype(val)
    @test eltype(typeof(x)) == eltype(val)
    @test length(x) == prod(dim)
    @test axes(x) == axes(ones(dim))
    @test x .* ones(dim) == fill(val, dim)
    @test vec(x) == fill(val, prod(dim))
    @test [v for v in x] == fill(val, prod(dim)) # iterator

    @test Base.IteratorSize(typeof(x)) == Base.HasLength()
    @test Base.IteratorEltype(typeof(x)) == typeof(val)
    @test Base.IndexStyle(typeof(x)) == IndexLinear()
    @test Base.firstindex(x) == 1
    @test Base.lastindex(x) == prod(dim)

    z = @inferred copy(x)
    @test z === x # ! (because isbitstype(ConstantArray{T,D}) is true)
	@test isbits(x)

    @test_throws BoundsError x[0]
    @test_throws BoundsError x[end+1]

    mask = @inferred ConstantArray(dim)
    @test eltype(mask) == Bool

    z = rand(dim...)
    @test z[mask] == vec(z)

    dim = (2^9,2^8)
    mask = @inferred ConstantArray(dim)
    x = rand(dim...)
    y = similar(vec(x))
    @test getindex!(y, x, mask) == vec(x)
    @test 100 > @allocated getindex!(y, x, mask)
end


text:
using Test
using KernelFunctions
using Distances
using Random

@testset "KernelFunctions" begin
include("test_kernelmatrix.jl")
include("test_approximations.jl")
include("test_constructors.jl")
# include("test_AD.jl")
include("test_transform.jl")
include("test_distances.jl")
include("test_kernels.jl")
include("test_generic.jl")
include("test_adjoints.jl")
include("test_custom.jl")
    #include("types.jl")
end


text:
# discrete-time case
"""
    PeriodicMatrix(M, T) -> A::PeriodicMatrix

Discrete-time periodic matrix representation.

The discrete-time periodic matrix object `A` is built from a
`p`-vector `M` of real matrices, the associated time period `T` and
the number of subperiods specified via the keyword argument `nperiod = k`.

`M` contains the cyclic component matrices `M[i]`, `i = 1,..., p`,
where `M[i]` represents the value `M(Δ(i-1))` of a time periodic matrix `M(t)`
of period `T/k`, with `Δ := T/(kp)`, the associated sample time.
It is assumed that `M[i] := M[mod(i-1,p)+1]` for arbitrary `i`.
All component matrices are allowed to have arbitrary (time-varying) dimensions.
The component matrices `M`, the period `T`, the number of subperiods `k`, the discrete period `p`
and the sample time `Δ` can be accessed via `A.M`, `A.period`, `A.nperiod`, `A.dperiod` and `A.Ts`, respectively.
"""
struct  PeriodicMatrix{Domain,T} <: AbstractPeriodicArray{Domain,T}
    M::Vector{Matrix{T}}
    period::Float64
    nperiod::Int
end
# additional constructors
function  PeriodicMatrix(M::Vector{MT}, period::Real; nperiod::Int = 1) where {MT <: Array}
   period > 0 || error("period must be positive")
   nperiod > 0 || error("number of subperiods must be positive")
   any(ndims.(M) .> 2) && error("only vectors with vector or matrix elements supported")
   m = size.(M,2)
   T = promote_type(eltype.(M)...)
   return any(m .== 1) ?  PeriodicMatrix{:d,T}([T.(reshape(M[i],size(M[i],1),m[i])) for i in 1:length(M)], Float64(period), nperiod)  :
                   PeriodicMatrix{:d,T}([T.(M[i]) for i in 1:length(M)], Float64(period), nperiod)
end
PeriodicMatrix{:d,T}(A::Vector{Matrix{T}}, period::Real; nperiod::Int = 1) where {T} =
      PeriodicMatrix{:d,T}(A, period, nperiod)

# period change
function PeriodicMatrix{:d,T}(A::PeriodicMatrix{:d,T1}, period::Real) where {T,T1}
   period > 0 || error("period must be positive")
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      PeriodicMatrix{:d,T}([T.(A.M[i]) for i in 1:length(A)], Aperiod*d, A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      PeriodicMatrix{:d,T}([T.(A.M[i]) for i in 1:length(A)], Aperiod/n, nperiod)
   end
end
set_period(A::AbstractVecOrMat,period::Real) = A
function set_period(A::AbstractPeriodicArray{:d,T}, period::Real) where T
   period > 0 || error("period must be positive")
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      PeriodicMatrix{:d,T}(A.M, Aperiod*d, A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      PeriodicMatrix{:d,T}(A.M, Aperiod/n, nperiod)
   end
end

# function PeriodicMatrix{:d,T}(A::Vector{Matrix{T1}}, period::Real) where {T,T1}
#     PeriodicMatrix([T.(A[i]) for i in 1:length(A)], period)
# end
# PeriodicMatrix(M::Vector{MT}, period::Real) where {T <: Real, MT <: Array{T}} =
#            PeriodicMatrix{:d,T}(M, period)
PeriodicMatrix(M::Vector{Matrix{T}}, period::Real; nperiod::Int = 1) where {T <: Real} =
       PeriodicMatrix{:d,T}(M, period, nperiod)
PeriodicMatrix(M::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Real} =
   PeriodicMatrix{:d,T}([reshape(M,size(M,1),size(M,2))], period, nperiod)
function Base.getproperty(A::PeriodicMatrix, d::Symbol)
   if d === :dperiod
      return length(getfield(A, :M))
   elseif d === :Ts
      return A.period/A.dperiod/A.nperiod
   else
      getfield(A, d)
   end
end
Base.propertynames(A::PeriodicMatrix) = (:dperiod, :Ts, fieldnames(typeof(A))...)
isconstant(A::PeriodicMatrix) = (A.dperiod == 1)
Base.size(A::PeriodicMatrix) = (size.(A.M,1),size.(A.M,2))
Base.size(A::PeriodicMatrix, d::Integer) = size.(A.M,d)
Base.length(A::PeriodicMatrix) = A.dperiod
Base.eltype(A::PeriodicMatrix{:d,T}) where T = T
"""
    PeriodicArray(M, T; nperiod = k) -> A::PeriodicArray

Discrete-time periodic array representation.

The discrete-time periodic array object `A` is built from a `m×n×p` real array
`M`, the associated time period `T` and the number of subperiods specified via
the keyword argument `nperiod = k`.
`M` contains the cyclic component matrices `M[:,:,i]`, `i = 1,..., p`,
where `M[:,:,i]` represents the value `M(Δ(i-1))` of a time periodic matrix `M(t)`
of period `T/k`, with `Δ := T/(kp)`, the associated sample time.
It is assumed that  `M[:,:,k] := M[:,:,mod(k-1,p)+1]` for arbitrary `k`.
The component matrices `M`, the period `T`, the number of subperiods `k`, the discrete period `p`
and the sample time `Δ` can be accessed via `A.M`, `A.period`, `A.nperiod`, `A.dperiod` and `A.Ts`, respectively.
"""
struct PeriodicArray{Domain,T} <: AbstractPeriodicArray{Domain,T}
    M::Array{T,3}
    period::Float64
    nperiod::Int
end
# additional constructors
function  PeriodicArray{:d,T}(M::Array{T,3}, period::Real; nperiod::Int = 1) where {T <: Real}
   period > 0 || error("period must be positive")
   nperiod > 0 || error("number of subperiods must be positive")
   PeriodicArray{:d,T}(M, Float64(period), nperiod)
end
function PeriodicArray{:d,T}(A::PeriodicArray{:d,T1}, period::Real) where {T,T1}
   period > 0 || error("period must be positive")
   #isconstant(A) && (return PeriodicArray{:d,T}(convert(Array{T,3},A.M), period; nperiod = 1))
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      PeriodicArray{:d,T}(convert(Array{T,3},A.M), Aperiod*d; nperiod = A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      PeriodicArray{:d,T}(convert(Array{T,3},A.M), Aperiod/n; nperiod)
   end
end

#PeriodicArray{:d,T}(M::Array{T1,3}, period::Real; nperiod::Int = 1) where {T,T1} = PeriodicArray(T.(M), period; nperiod)
PeriodicArray(M::Array{T,3}, period::Real; nperiod::Int = 1) where {T <: Real} = PeriodicArray{:d,T}(M, period, nperiod)
PeriodicArray(M::VecOrMat{T}, period::Real; nperiod::Int = 1) where T = PeriodicArray(reshape(M,size(M,1),size(M,2),1), period; nperiod)
function Base.getproperty(A::PeriodicArray, d::Symbol)
   if d === :dperiod
      return size(getfield(A, :M), 3)
   elseif d === :Ts
      return A.period/A.dperiod/A.nperiod
   else
      getfield(A, d)
   end
end
Base.propertynames(A::PeriodicArray) = (:dperiod, :Ts, fieldnames(typeof(A))...)
isconstant(A::PeriodicArray) = (A.dperiod == 1)
Base.size(A::PeriodicArray) = (size(A.M,1),size(A.M,2))
Base.size(A::PeriodicArray, d::Integer) = size(A.M,d)
Base.length(A::PeriodicArray) = A.dperiod
Base.eltype(A::PeriodicArray{:d,T}) where T = T
iscontinuous(A::AbstractPeriodicArray) = typeof(A).parameters[1] == :c
#iscontinuous(A) = typeof(A).parameters[1] == :c
iscontinuous(A::Type) = A.parameters[1] == :c
"""
    PeriodicFunctionMatrix(f, T; nperiod = k) -> A::PeriodicFunctionMatrix

Continuous-time periodic function matrix representation.

The continuous-time periodic function matrix object `A` is built from a
time periodic real matrix function `f(t)` of real time variable `t`,
the associated time period `T` and the associated number of subperiods
specified via the keyword argument `nperiod = k`.
It is assumed that  `F(t) = F(t+T/k)` for any real time value `t`.
The function `f(t)`, the period `T`, the row and column dimensions
of `f(t)`, the number of subperiods `k` can be accessed via `A.f`, `A.period`,
the tuple `A.dims` and `A.nperiod`, respectively.
"""
struct PeriodicFunctionMatrix{Domain,T} <: AbstractPeriodicArray{Domain,T}
   f::Function
   period::Float64
   dims::Tuple{Int,Int}
   nperiod::Int
   _isconstant::Bool
end
# additional constructors
function PeriodicFunctionMatrix{:c,Tf}(f::Function, period::Real; isconst::Bool = false, nperiod::Int = 1) where {Tf}
   period > 0 || error("period must be positive")
   nperiod > 0 || error("number of subperiods must be positive")
   F0 = f(period)
   nd = ndims(F0)
   nd <= 2 || error("two-dimensional function array expected, got an $nd -dimensional array")
   m, n = size(F0,1),size(F0,2)
   eltype(F0) == Tf ? PeriodicFunctionMatrix{:c,Tf}(t -> n == 1 ? reshape(f(t),m,n) : f(t), Float64(period), (m,n), nperiod, isconst) :
                      PeriodicFunctionMatrix{:c,Tf}(t -> n == 1 ? convert(Matrix{Tf},reshape(f(Tf(t)),m,n)) : convert(Matrix{Tf},f(Tf(t))), Float64(period), (m,n), nperiod, isconst)
end
PeriodicFunctionMatrix(f::F, period::Real; isconst::Bool = false, nperiod::Int = 1) where {F<:Function}  =
             PeriodicFunctionMatrix{:c,eltype(f(period))}(f, period; isconst, nperiod)
function PeriodicFunctionMatrix(A::VecOrMat{T}, period::Real) where {T <: Real}
   if T == Num
      @variables t
      f = eval(build_function(reshape(A,size(A,1),size(A,2)), t, expression=Val{false})[1])
      PeriodicFunctionMatrix{:c,Float64}(t -> f(t), period; isconst = false)
   else
      PeriodicFunctionMatrix{:c,T}(t -> reshape(A,size(A,1),size(A,2)), period; isconst = true)
   end
end
function PeriodicFunctionMatrix{:c,T}(at::PeriodicFunctionMatrix, period::Real) where {T}
   period > 0 || error("period must be positive")
   Aperiod = at.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      PeriodicFunctionMatrix{:c,T}(at.f, Aperiod*d, at.dims, at.nperiod*d, at._isconstant)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      PeriodicFunctionMatrix{:c,T}(at.f, Aperiod/n, at.dims, at.nperiod, at._isconstant)
   end
end
#set_period(A::PeriodicFunctionMatrix, period::Real) = PeriodicFunctionMatrix{:c,eltype(A)}(A,period)

# function PeriodicFunctionMatrix(at::PeriodicFunctionMatrix, period::Real = at.period; isconst::Bool = isconstant(at))
#    # at.period = period
#    # at._isconstant = isconst
#    # return at
#    return PeriodicFunctionMatrix(at.f, period; isconst)
# end
# properties
isconstant(A::PeriodicFunctionMatrix) = A._isconstant
function isperiodic(f::Function, period::Real)
   t = rand(typeof(period))
   return f(t) ≈ f(t+period)
end
isperiodic(A::PeriodicFunctionMatrix) = isconstant(A) ? true : isperiodic(A.f,A.period/A.nperiod)
Base.size(A::PeriodicFunctionMatrix) = A.dims
Base.size(A::PeriodicFunctionMatrix, d::Integer) = d <= 2 ? size(A)[d] : 1
Base.eltype(A::PeriodicFunctionMatrix{:c,T}) where T = T

"""
    PeriodicSymbolicMatrix(F, T; nperiod = k) -> A::PeriodicSymbolicMatrix

Continuous-time periodic symbolic matrix representation.

The continuous-time periodic symbolic matrix object `A` is built from `F`, a
symbolic periodic real matrix or vector of symbolic variable `t`,
the associated time period `T` and the associated number of subperiods
specified via the keyword argument `nperiod = k`.
It is assumed that  `F(t) = F(t+T/k)` for any real time value `t`.
The symbolic matrix `F`, the period `T` and the number of subperiods `k`
can be accessed via `A.F`, `A.period` and `A.nperiod`, respectively.
"""
struct PeriodicSymbolicMatrix{Domain,T} <: AbstractPeriodicArray{Domain,T}
   F::Matrix{<:Num}
   period::Float64
   nperiod::Int
end
# additional constructors
function  PeriodicSymbolicMatrix{:c,T}(F::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Num}
   period > 0 || error("period must be positive")
   nperiod > 0 || error("number of subperiods must be positive")
   # check that array F is depending only on t
   tt = rand()
   @variables t
   Ft = substitute.(F, (Dict(t => tt),))
   m, n = size(Ft,1), size(Ft,2)
   any(length.(Symbolics.get_variables.(Ft)) .> 0 ) && error("t must be the only variable in F")
   PeriodicSymbolicMatrix{:c,T}(n == 1 ? reshape(F,m,n) : F, Float64(period), nperiod)
end
PeriodicSymbolicMatrix(F::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Union{Num,Real}} =
    PeriodicSymbolicMatrix{:c,Num}(Num.(F), period; nperiod)
# PeriodicSymbolicMatrix(F::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Real} =
#     PeriodicSymbolicMatrix{:c,Num}(Num.(F), Float64(period), nperiod)
# PeriodicSymbolicMatrix{:c,Num}(F::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Number} =
#     PeriodicSymbolicMatrix{:c,Num}(Num.(F), Float64(period), nperiod)
function PeriodicSymbolicMatrix{:c,T}(A::PeriodicSymbolicMatrix, period::Real) where {T}
   period > 0 || error("period must be positive")
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      PeriodicSymbolicMatrix{:c,T}(A.F, Aperiod*d, A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      PeriodicSymbolicMatrix{:c,T}(A.F, Aperiod/n, A.nperiod)
   end
end
set_period(A::PM, period::Real) where {T, PM <: AbstractPeriodicArray{:c,T}} = PM(A,period)

# properties
isconstant(A::PeriodicSymbolicMatrix) = all(length.(Symbolics.get_variables.(A.F)) .== 0)
function isperiodic(A::VecOrMat{T}, period::Real) where {T <: Num}
   tt = rand()
   @variables t
   At = substitute.(A, (Dict(t => tt),))
   return norm(At - substitute.(A, (Dict(t => tt+period),))) <= eps(10.)*max(1.,norm(At))
end
isperiodic(A::PeriodicSymbolicMatrix) = isconstant(A) ? true : isperiodic(A.F,A.period)
Base.size(A::PeriodicSymbolicMatrix) = size(A.F)
Base.size(A::AbstractPeriodicArray, d::Integer) = d <= 2 ? size(A)[d] : 1
Base.eltype(A::PeriodicSymbolicMatrix{:c,T}) where T = T

struct FourierFunctionMatrix{Domain,T} <: AbstractPeriodicArray{Domain,T}
   M::Fun
   period::Float64
   nperiod::Int
end
# additional constructors
"""
     FourierFunctionMatrix(Afun, T) -> A::FourierFunctionMatrix

Continuous-time Fourier function matrix representation.

The Fourier function matrix object `A` of period `T` is built using
the Fourier series representation of a periodic matrix `Afun(t)` of subperiod `T′ = T/k`,
where each entry of `Afun(t)` has the form

             p
      a_0 +  ∑ ( ac_i*cos(i*t*2*π/T′)+as_i*sin(i*2*π*t/T′) ) ,
            i=1

where `k ≥ 1` is the number of subperiods (default: `k = 1`).
The matrix `Afun` containing the Fourier representation, the period `T` and the
number of subperiods `k` can be accessed via `A.M`, `A.period` and `A.nperiod`, respectively.
"""
function FourierFunctionMatrix{:c,T}(A::Fun, period::Real) where {T}
   period > 0 || error("period must be positive")
   n, m = size(A,1), size(A,2)
   sint = domain(A)
   (sint.a == 0 && sint.b > 0) || error("the domain must be of the form 0..period")
   ti = rationalize(period/sint.b)
   denominator(ti) == 1 || error("only integer multiple periods supported")
   FourierFunctionMatrix{:c,eltype(domain(A))}(m == 1 ? reshape(A,n,m) : A, Float64(period), numerator(ti))
end
FourierFunctionMatrix(A::Fun, period::Real)  =
       FourierFunctionMatrix{:c,eltype(domain(A))}(A::Fun, period::Real)
FourierFunctionMatrix(A::Fun)  = FourierFunctionMatrix{:c,eltype(domain(A))}(A::Fun, domain(A).b)
function FourierFunctionMatrix{:c,T}(A::FourierFunctionMatrix, period::Real) where {T}
   period > 0 || error("period must be positive")
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      FourierFunctionMatrix{:c,T}(A.M, Aperiod*d, A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      FourierFunctionMatrix{:c,T}(A.M, Aperiod/n, A.nperiod)
   end
end
FourierFunctionMatrix(A0::VecOrMat{T}, period::Real) where {T <: Real}  =
    FourierFunctionMatrix{:c,float(T)}(Fun(t->float(T).(A0),Fourier(0..period)), period)
FourierFunctionMatrix{:c,T}(A0::VecOrMat, period::Real) where {T <: Real}  =
    FourierFunctionMatrix{:c,float(T)}(Fun(t->float(T).(A0),Fourier(0..period)), period)
function isconstant(A::FourierFunctionMatrix)
   for i = 1:size(A.M,1)
       for j = 1: size(A.M,2)
           ncoefficients(chop(A.M[i,j])) <= 1 || (return false)
       end
   end
   return true
end
isperiodic(A::FourierFunctionMatrix) = true
Base.size(A::FourierFunctionMatrix) = size(A.M)
Base.eltype(A::FourierFunctionMatrix{:c,T}) where T = T

struct HarmonicArray{Domain,T} <: AbstractPeriodicArray{Domain,T}
   values::Array{Complex{T},3}
   period::Float64
   nperiod::Int
end
# additional constructors
"""
     HarmonicArray(Ahr, T; nperiod = k) -> A::HarmonicArray

Continuous-time harmonic array representation.

The harmonic array object `A` of period `T` is built using
the harmonic representation of a periodic matrix `Ahr(t)` of subperiod `T′ = T/k` in the form

                     p
     Ahr(t) = A_0 +  ∑ ( Ac_i*cos(i*t*2*π/T′)+As_i*sin(i*2*π*t/T′) ) ,
                    i=1

where `k ≥ 1` is the number of subperiods (default: `k = 1`).
The `m×n×(p+1)` complex array `Ahr` contains the harmonic components as follows:
`Ahr[:,:,1]` contains the constant term `A_0` (the mean value) and
the real and imaginary parts of `Ahr[:,:,i+1]`
for `i = 1, ..., p` contain the coefficient matrices `Ac_i` and `As_i`, respectively.
The complex matrix `Ahr` containing the harmonic components, the period `T` and the
number of subperiods `k` can be accessed via `A.values`, `A.period` and `A.nperiod`, respectively.
"""
function HarmonicArray{:c,T}(Ahr::Array{<:Complex,3}, period::Real; nperiod::Int = 1) where {T}
   period > 0 || error("period must be positive")
   nperiod > 0 || error("number of subperiods must be positive")
   (size(Ahr,3) > 0 && iszero(imag(view(Ahr,:,:,1)))) || error("imaginary part of constant term must be zero")
   HarmonicArray{:c,T}(convert(Array{Complex{T},3},Ahr), Float64(period), nperiod)
end
HarmonicArray(Ahr::Array{Complex{T},3}, period::Real; nperiod::Int = 1) where T = HarmonicArray{:c,T}(Ahr, period; nperiod)
# change period and type
function HarmonicArray{:c,T}(A::HA, period::Real) where {HA <: HarmonicArray} where {T}
   period > 0 || error("period must be positive")
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      HarmonicArray{:c,T}(convert(Array{Complex{T},3},A.values), Aperiod*d; nperiod = A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      HarmonicArray{:c,T}(convert(Array{Complex{T},3},A.values), Aperiod/n; nperiod)
   end
end
"""
     HarmonicArray(A0, Ac, As, T) -> A::HarmonicArray

Construct a harmonic array representation from the harmonic components.

The harmonic array object `A` is built for
the harmonic representation `Ahr(t)` of a periodic matrix of period `T` in the form

                     p
     Ahr(t) = A_0 +  ∑ ( Ac_i*cos(i*t*2*π/T)+As_i*sin(i*2*π*t/T) ) ,
                    i=1

where the constant term `A_0` is contained in the real matrix `A0`, and `Ac` and `As` are
vectors of real matrices such that the `i`-th (cosinus) coefficient matrix
`Ac_i` is contained in `Ac[i]` and the `i`-th (sinus) coefficient matrix
`As_i` is contained in `As[i]`. `p` is the maximum of length of the vectors of matrices `Ac` and `As`.
If the length of `Ac` or `As` is less than `p`, then zero trailing matrices are assumed in the respective matrix.
All component matrices must have the same dimensions.
The complex matrix containing the harmonic components and the period `T`
can be accessed via `A.values` and `A.period`, respectively.
"""
function HarmonicArray(A0::MT, Acos::Union{Vector{MT},Nothing},
                               Asin::Union{Vector{MT},Nothing}, period::Real; nperiod::Int = 1) where {T <: Real, MT <: VecOrMat{T}}
   nc = isnothing(Acos) ? 0 : length(Acos)
   ns = isnothing(Asin) ? 0 : length(Asin)
   nmin = min(nc,ns)
   N = max(1,max(nc,ns)+1)
   ahr = Array{Complex{T},3}(undef, size(A0,1), size(A0,2), N)
   ahr[:,:,1] = A0
   [ahr[:,:,i+1] = complex.(Acos[i],Asin[i]) for i in 1:nmin]
   [ahr[:,:,i+1] = Acos[i] for i in nmin+1:nc]
   [ahr[:,:,i+1] = im*Asin[i] for i in nmin+1:ns]
   #HarmonicArray(ahr, period)
   HarmonicArray{:c,T}(ahr, period, nperiod)
end
HarmonicArray(A0::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Real}  =
          HarmonicArray(complex(reshape(A0,size(A0,1),size(A0,2),1)), period; nperiod)
HarmonicArray(A0::VecOrMat{T}, Acos::Vector{MT}, period::Real; nperiod::Int = 1) where {T <: Real, MT <: VecOrMat{T}}  =
          HarmonicArray(A0, Acos, nothing, period; nperiod)

# properties
isconstant(A::HarmonicArray) = size(A.values,3) <= 1
isperiodic(A::HarmonicArray) = true
Base.size(A::HarmonicArray) = (size(A.values,1),size(A.values,2))
Base.eltype(A::HarmonicArray{:c,T}) where T = T

"""
    PeriodicTimeSeriesMatrix(At, T; nperiod = k) -> A::PeriodicTimeSeriesMatrix

Continuous-time periodic time series matrix representation.

The continuous-time periodic time series matrix object `A` of period `T` is built from a
`p`-vector `At` of real matrices and the associated subperiod `T′ = T/k`, where
`k ≥ 1` is the number of subperiods (default: `k = 1`).
`At` contains the cyclic component matrices `At[i]`, `i = 1,..., p`,
where `At[i]` represents the value `A(Δ(i-1))` of a time periodic matrix `A(t)`
of period `T′`, with `Δ := T′/p`, the associated sampling time.
It is assumed that `At[i] := At[mod(i-1,p)+1]` for arbitrary `i`.
All component matrices must have the same dimensions.
The component matrices `At`, the period `T` and the number of subperiods `k`
can be accessed via `A.values`, `A.period`, and `A.nperiod`, respectively.
"""
struct PeriodicTimeSeriesMatrix{Domain,T} <: AbstractPeriodicArray{Domain,T}
   values::Vector{Array{T,2}}
   period::Float64
   nperiod::Int
end
# additional constructors
function PeriodicTimeSeriesMatrix{:c,T}(At::Union{Vector{Vector{T}},Vector{Matrix{T}}}, period::Real; nperiod::Int = 1) where {T <: Real}
   period > 0 || error("period must be positive")
   nperiod > 0 || error("number of subperiods must be positive")
   N = length(At)
   #N <= 1 && (return PeriodicTimeSeriesMatrix{:c,T}(At, Float64(period); nperiod) ) # the constant matrix case
   n1, n2 = size(At[1],1), size(At[1],2)
   (all(size.(At,1) .== n1) && all(size.(At,2) .== n2)) || error("all component matrices must have the same dimensions")

   # adjust final data to matrix type
   PeriodicTimeSeriesMatrix{:c,T}(n2 == 1 ? [reshape(At[j],n1,n2) for j in 1:N] : At, Float64(period), nperiod)
end
PeriodicTimeSeriesMatrix{:c,T}(A::Vector{Matrix{T1}}, period::Real; nperiod::Int = 1) where {T,T1} =
   PeriodicTimeSeriesMatrix([T.(A[i]) for i in 1:length(A)], period; nperiod)
PeriodicTimeSeriesMatrix(At::Union{Vector{Vector{T}},Vector{Matrix{T}}}, period::Real; nperiod::Int = 1) where {T <: Real} =
     PeriodicTimeSeriesMatrix{:c,T}(At, period; nperiod)
PeriodicTimeSeriesMatrix(At::VecOrMat{T}, period::Real; nperiod::Int = 1) where {T <: Real}  =
        PeriodicTimeSeriesMatrix([reshape(At,size(At,1),size(At,2))], period; nperiod)
# period change
function PeriodicTimeSeriesMatrix{:c,T}(A::PeriodicTimeSeriesMatrix{:c,T1}, period::Real) where {T,T1}
   Aperiod = A.period
   r = rationalize(Aperiod/period)
   n, d = numerator(r), denominator(r)
   min(n,d) == 1 || error("new period is incommensurate with the old period")
   if period >= Aperiod
      PeriodicTimeSeriesMatrix{:c,T}([T.(A.values[i]) for i in 1:length(A)], Aperiod*d; nperiod = A.nperiod*d)
   elseif period < Aperiod
      nperiod = div(A.nperiod,n)
      nperiod < 1 && error("new period is incommensurate with the old period")
      PeriodicTimeSeriesMatrix{:c,T}([T.(A.values[i]) for i in 1:length(A)], Aperiod/n; nperiod)
   end
end

# properties
isconstant(At::PeriodicTimeSeriesMatrix) = length(At.values) <= 1
isperiodic(At::PeriodicTimeSeriesMatrix) = true
Base.length(At::PeriodicTimeSeriesMatrix) = length(At.values)
Base.size(At::PeriodicTimeSeriesMatrix) = length(At) > 0 ? size(At.values[1]) : (0,0)
Base.eltype(At::PeriodicTimeSeriesMatrix{:c,T}) where T = T

# conversions to discrete-time PeriodicMatrix
Base.convert(::Type{PeriodicMatrix}, A::PeriodicArray{:d,T}) where T =
             PeriodicMatrix{:d,T}([A.M[:,:,i] for i in 1:size(A.M,3)],A.period; nperiod = A.nperiod)

# conversions to discrete-time PeriodicArray
# function Base.convert(::Type{PeriodicArray}, A::PeriodicMatrix{:d,T}) where T
#     N = length(A)
#     m, n = size.(A,1), size.(A,2)
#     N == 0 && PeriodicArray(Array{T,3}(undef,0,0,0),A.period)
#     if any(m .!= m[1]) || any(m .!= m[1])
#        @warn "Non-constant dimensions: the resulting component matrices padded with zeros"
#        t = zeros(T,maximum(m),maximum(n),N)
#        [copyto!(view(t,1:m[i],1:n[i],i),A[i]) for i in 1:N]
#        PeriodicArray{:d,T}(t,A.period)
#     else
#        t = zeros(T,m[1],n[1],N)
#        [copyto!(view(t,:,:,i),A[i]) for i in 1:N]
#        PeriodicArray{:d,T}(t,A.period)
#     end
# end
Base.convert(::Type{PeriodicArray}, A::PeriodicMatrix) where T = pm2pa(A)

# conversions to continuous-time PeriodicFunctionMatrix
function Base.convert(::Type{PeriodicFunctionMatrix{:c,T}}, A::PeriodicFunctionMatrix) where T
   return eltype(A) == T ? A : PeriodicFunctionMatrix{:c,T}(x -> T.(A.f(T(x))), A.period, A.dims, A.nperiod, A._isconstant)
end

function Base.convert(::Type{PeriodicFunctionMatrix{:c,T}}, A::FourierFunctionMatrix) where T
   return eltype(A) == T ? A : PeriodicFunctionMatrix{:c,T}(x -> T.(A.M(T(x))), A.period, size(A), A.nperiod, isconstant(A))
end
function Base.convert(::Type{PeriodicFunctionMatrix}, A::FourierFunctionMatrix)
   return PeriodicFunctionMatrix{:c,eltype(A)}(x -> A.M(x), A.period, size(A.M), A.nperiod, isconstant(A))
end

function Base.convert(::Type{PeriodicFunctionMatrix{:c,T}}, A::PeriodicSymbolicMatrix) where T
   @variables t
   f = eval(build_function(A.F, t, expression=Val{false})[1])
   PeriodicFunctionMatrix{:c,T}(x -> f(T(x)), A.period, size(A), A.nperiod, isconstant(A))
end
function Base.convert(::Type{PeriodicFunctionMatrix}, A::PeriodicSymbolicMatrix) where T
   @variables t
   f = eval(build_function(A.F, t, expression=Val{false})[1])
   PeriodicFunctionMatrix{:c,Float64}(x -> f(x), A.period, size(A), A.nperiod, isconstant(A))
end
# function Base.convert(::Type{PeriodicFunctionMatrix{:c,T}}, A::PeriodicSymbolicMatrix) where T
#    @variables t
#    f = eval(build_function(A.F, t, expression=Val{false})[1])
#    PeriodicFunctionMatrix{:c,T}(x -> f(x), A.period, size(A), A.nperiod, isconstant(A))
# end

# PeriodicFunctionMatrix(ahr::HarmonicArray, period::Real = ahr.period; exact = true)  =
#           PeriodicFunctionMatrix(t::Real -> hreval(ahr,t;exact)[1], period)
Base.convert(::Type{PeriodicFunctionMatrix}, ahr::HarmonicArray)  where T =
         PeriodicFunctionMatrix{:c,real(eltype(ahr.values))}(t::Real -> hreval(ahr,t), ahr.period, size(ahr), ahr.nperiod, isconstant(ahr))
Base.convert(::Type{PeriodicFunctionMatrix{:c,T}}, ahr::HarmonicArray)  where T =
         PeriodicFunctionMatrix{:c,T}(t::Real -> hreval(ahr,T(t)), ahr.period, size(ahr), ahr.nperiod, isconstant(ahr))
Base.convert(::Type{PeriodicFunctionMatrix}, At::PeriodicTimeSeriesMatrix)  where T =
    ts2pfm(At; method = "cubic")
# function PeriodicFunctionMatrix(A::PeriodicTimeSeriesMatrix; method = "linear")
#    N = length(A.values)
#    N == 0 && error("empty time array")
#    N == 1 && (return t -> A.values[1])
#    #dt = A.time[2]-A.time[1]
#    dt = A.period/N
#    ts = (0:N)*dt
#    if method == "linear"
#       itp = LinearInterpolation(ts, push!(copy(A.values),A.values[1]))
#       return PeriodicFunctionMatrix(t -> itp(mod(t, A.period)), A.period )
#    elseif method == "cubic"
#       n1, n2 = size(A.values[1])
#       intparray = Array{Any,2}(undef,n1, n2)
#       [intparray[i,j] = CubicSplineInterpolation(ts,[getindex.(A.values,i,j);A.values[1][i,j]],bc=Line(OnGrid())) for i in 1:n1, j in 1:n2]
#       return PeriodicFunctionMatrix(t -> [intparray[i,j](mod(t, A.period)) for i in 1:n1, j in 1:n2 ], A.period )
#    end
# end

# conversions to continuous-time Fourier function matrix
function Base.convert(::Type{FourierFunctionMatrix}, A::PeriodicFunctionMatrix)
   T = float(eltype(A))
   return FourierFunctionMatrix{:c,T}(Fun(x -> T.(A.f(x)), Fourier(0..A.period/A.nperiod)), Float64(A.period), A.nperiod)
end
function Base.convert(::Type{FourierFunctionMatrix}, A::HarmonicArray)
   tA = convert(PeriodicFunctionMatrix,A)
   return FourierFunctionMatrix{:c,eltype(tA)}(Fun(x -> tA.f(x), Fourier(0..tA.period/tA.nperiod)), Float64(tA.period), tA.nperiod)
end
function Base.convert(::Type{FourierFunctionMatrix}, A::PeriodicSymbolicMatrix)
   tA = convert(PeriodicFunctionMatrix,A)
   return FourierFunctionMatrix{:c,eltype(tA)}(Fun(x -> tA.f(x), Fourier(0..tA.period/tA.nperiod)), Float64(tA.period), tA.nperiod)
end
function Base.convert(::Type{FourierFunctionMatrix}, A::PeriodicTimeSeriesMatrix)
   tA = convert(PeriodicFunctionMatrix,A)
   return FourierFunctionMatrix{:c,eltype(tA)}(Fun(x -> tA.f(x), Fourier(0..tA.period/tA.nperiod)), Float64(tA.period), tA.nperiod)
end

# conversions to continuous-time PeriodicSymbolicMatrix
function Base.convert(::Type{PeriodicSymbolicMatrix}, A::PeriodicFunctionMatrix) where T
   @variables t
   PeriodicSymbolicMatrix(Num.(A.f(t)), A.period; nperiod = A.nperiod)
end
function Base.convert(::Type{PeriodicSymbolicMatrix{:c,T}}, A::PeriodicFunctionMatrix) where T
   @variables t
   PeriodicSymbolicMatrix(Num.(A.f(t)), A.period; nperiod = A.nperiod)
end
Base.convert(::Type{PeriodicSymbolicMatrix}, ahr::HarmonicArray)  where T =
   PeriodicSymbolicMatrix(hr2psm(ahr), ahr.period; nperiod = ahr.nperiod)
Base.convert(::Type{PeriodicSymbolicMatrix}, A::PeriodicTimeSeriesMatrix) =
   PeriodicSymbolicMatrix(hr2psm(convert(HarmonicArray,A)), A.period; nperiod = A.nperiod)
Base.convert(::Type{PeriodicSymbolicMatrix}, A::FourierFunctionMatrix) =
   convert(PeriodicSymbolicMatrix,convert(HarmonicArray,A))

# conversion to continuous-time HarmonicArray
Base.convert(::Type{HarmonicArray}, A::PeriodicTimeSeriesMatrix) = ts2hr(A)
Base.convert(::Type{HarmonicArray}, A::PeriodicFunctionMatrix) = pfm2hr(A)
Base.convert(::Type{HarmonicArray}, A::PeriodicSymbolicMatrix) = psm2hr(A)
#Base.convert(::Type{HarmonicArray}, A::FourierFunctionMatrix) = pfm2hr(convert(PeriodicFunctionMatrix,A))
Base.convert(::Type{HarmonicArray}, A::FourierFunctionMatrix) = ffm2hr(A)

# conversions to PeriodicTimeSeriesMatrix
function Base.convert(::Type{PeriodicTimeSeriesMatrix}, A::PeriodicSymbolicMatrix)
    tA = convert(PeriodicFunctionMatrix,A)
    PeriodicTimeSeriesMatrix(tA.f.((0:127)*tA.period/128/tA.nperiod), tA.period; nperiod = tA.nperiod)
end
Base.convert(::Type{PeriodicTimeSeriesMatrix}, A::PeriodicFunctionMatrix) =
         PeriodicTimeSeriesMatrix(A.f.((0:127)*A.period/128/A.nperiod), A.period; nperiod = A.nperiod)
Base.convert(::Type{PeriodicTimeSeriesMatrix}, A::HarmonicArray) =
    PeriodicTimeSeriesMatrix(tvmeval(A, collect((0:127)*A.period/128/A.nperiod)), A.period; nperiod = A.nperiod)
Base.convert(::Type{PeriodicTimeSeriesMatrix}, A::PeriodicMatrix) =
    convert(PeriodicTimeSeriesMatrix, pm2pa(A))
Base.convert(::Type{PeriodicTimeSeriesMatrix}, A::PeriodicArray) =
    PeriodicTimeSeriesMatrix([A.M[:,:,i] for i in 1:size(A.M,3)], A.period; nperiod = A.nperiod)
Base.convert(::Type{PeriodicTimeSeriesMatrix}, A::FourierFunctionMatrix) =
   convert(PeriodicTimeSeriesMatrix,convert(PeriodicFunctionMatrix,A))


text:
using BigArrays.S3Dicts
using Test

@testset "test s3 IO" begin
    as3  = S3Dict("s3://seunglab/jpwu/test/image/4_4_40")

    a = rand(UInt8, 50)

    println("save object")
    as3["test"] = a
    println("get object")
    b = as3["test"]
    @test all(a.==b)
    @test haskey(as3, "test") == true
    @test haskey(as3, "test2") == false
    @info("delete the file in s3")
    delete!(as3, "test")

    println("test no such key error...")
    try
        as3["test"]
    catch err
        @test isa(err, KeyError)
    end
end


text:
type DiscreteHopfieldNet <: HopfieldNet
    s::Vector{Float64} # State -- could have used Int's
    W::Matrix{Float64} # Weights
end

function DiscreteHopfieldNet(n::Integer)
    s = ones(n)
    W = zeros(n, n)
    DiscreteHopfieldNet(s, W)
end

# Perform one asynchronous update on randomly selected neuron
function update!(net::DiscreteHopfieldNet)
    i = rand(1:length(net.s))
    net.s[i] = dot(net.W[:, i], net.s) > 0 ? +1 : -1
    return
end

function Base.show(io::IO, net::DiscreteHopfieldNet)
    @printf io "A discrete Hopfield net with %d neurons\n" length(net.s)
end


text:
using VoronoiCells
using GeometryBasics
using Random
using Test

@testset "Cells" begin
    @testset "PointCollection from GeometryBasics" begin
        Random.seed!(2)
        points = [Point2(rand(), rand()) for _ in 1:5]
        rect = Rectangle(Point2(0, 0), Point2(1, 1))

        pc = VoronoiCells.PointCollection(points, rect)

        @test isa(pc, VoronoiCells.PointCollection{GeometryBasics.Point2{Float64}})
    end

    @testset "PointCollection from VoronoiDelaunay" begin
        Random.seed!(2)
        points = [VoronoiDelaunay.Point2D(1 + rand(), 1 + rand()) for _ in 1:5]
        rect = Rectangle(VoronoiDelaunay.Point2D(1.1, 1.1), VoronoiDelaunay.Point2D(1.9, 1.9))

        pc = VoronoiCells.PointCollection(points, rect)

        @test isa(pc, VoronoiCells.PointCollection{VoronoiDelaunay.Point2D})
    end

    @testset "Error making PointCollection with mix of VoronoiDelaunay and GeometryBasics" begin
        Random.seed!(2)
        points = [VoronoiDelaunay.Point2D(1 + rand(), 1 + rand()) for _ in 1:5]
        rect = Rectangle(Point2(1, 1), Point2(2, 2))

        @test_throws MethodError VoronoiCells.PointCollection(points, rect)
    end

    @testset "Simple point set far from corners" begin
        rect = Rectangle(Point2(0, 0), Point2(1, 1))

        points = [
            Point2(0.75, 0.75),
            Point2(0.25, 0.25),
            Point2(0.25, 0.75)
        ]

        tess = voronoicells(points, rect)

        @test tess.Cells[1] == [Point2(0.5, 1.0), Point2(0.5, 0.5), Point2(1.0, 0.0), Point2(1.0, 1.0)]
        @test tess.Cells[2] == [Point2(0.0, 0.5), Point2(0.0, 0.0), Point2(1.0, 0.0), Point2(0.5, 0.5)]
        @test tess.Cells[3] == [Point2(0.0, 1.0), Point2(0.0, 0.5), Point2(0.5, 0.5), Point2(0.5, 1.0)]
    end

    @testset "Point set on corners" begin
        rect = Rectangle(Point2(0, 0), Point2(1, 1))

        points = [
            Point2(0.0, 0.0),
            Point2(1.0, 0.0),
            Point2(1.0, 1.0),
            Point2(0.0, 1.0)
        ]

        tess = voronoicells(points, rect)

        @test tess.Cells[1] == [Point2(0.0, 0.5), Point2(0.0, 0.0), Point2(0.5, 0.0), Point2(0.5, 0.5)]
        @test tess.Cells[2] == [Point2(0.5, 0.5), Point2(0.5, 0.0), Point2(1.0, 0.0), Point2(1.0, 0.5)]
        @test tess.Cells[3] == [Point2(0.5, 1.0), Point2(0.5, 0.5), Point2(1.0, 0.5), Point2(1.0, 1.0)]
        @test tess.Cells[4] == [Point2(0.0, 1.0), Point2(0.0, 0.5), Point2(0.5, 0.5), Point2(0.5, 1.0)]
    end

    @testset "Random point set" begin
        Random.seed!(1)
        points = [Point2(rand(), rand()) for _ in 1:5]

        rect = Rectangle(Point2(0, 0), Point2(1, 1))

        tess = voronoicells(points, rect)

        @test tess.Cells[1] ≈ [
            Point2(0.0, 0.67254), Point2(0.0, 0.11508), Point2(0.31246, 0.18583), Point2(0.56624, 0.65872),
        ] atol = 1e-4
        @test tess.Cells[2] ≈ [
         Point2(0.0, 0.11508), Point2(0.0, 0.0), Point2(0.52699, 0.0), Point2(0.31246, 0.18583)
        ] atol = 1e-4
        @test tess.Cells[3] ≈ [
             Point2(0.60787, 0.67143), Point2(0.56624, 0.65872), Point2(0.31246, 0.18583), Point2(0.52699, 0.0), Point2(1.0, 0.0), Point2(1.0, 0.44116)
        ] atol = 1e-4
        @test tess.Cells[4] ≈ [
             Point2(0.60166, 1.0), Point2(0.60787, 0.67143), Point2(1.0, 0.44116), Point2(1.0, 1.0)
        ] atol = 1e-4
        @test tess.Cells[5] ≈ [
             Point2(0.0, 1.0), Point2(0.0, 0.67254), Point2(0.56624, 0.65872), Point2(0.60787, 0.67143), Point2(0.60166, 1.0)
        ] atol = 1e-4
    end
end


text:
module CCBladeOpenMDAOExample

using CCBlade: Rotor, Section, OperatingPoint

# convenience function to set fields within an array of structs
function Base.setproperty!(obj::Array{Rotor{TF, TI, TB, TF2}, N}, sym::Symbol, x) where {TF, TI, TB, TF2, N}
    setfield!.(obj, sym, x)
end

# convenience function to set fields within an array of structs
function Base.setproperty!(obj::Array{Section{TF1, TF2, TF3, TAF}, N}, sym::Symbol, x) where {TF1, TF2, TF3, TAF, N}
    setfield!.(obj, sym, x)
end

# convenience function to set fields within an array of structs
function Base.setproperty!(obj::Array{OperatingPoint{TF, TF2}, N}, sym::Symbol, x) where {TF, TF2, N}
    setfield!.(obj, sym, x)
end

include("derivatives.jl")
include("openmdao.jl")

end # module


text:
abstract type AbstractStorageFormulation <: AbstractDeviceFormulation end

struct BookKeeping <: AbstractStorageFormulation end

struct BookKeepingwReservation <: AbstractStorageFormulation end

#################################################Storage Variables#################################

function active_power_variables!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{St}) where {St<:PSY.Storage}

    add_variable(canonical,
                 devices,
                 Symbol("Psin_$(St)"),
                 false,
                 :nodal_balance_active,
                 -1.0;
                 lb_value = d -> 0.0,)
    add_variable(canonical,
                 devices,
                 Symbol("Psout_$(St)"),
                 false,
                 :nodal_balance_active;
                 lb_value = d -> 0.0,)

    return

end

function reactive_power_variables!(canonical::Canonical,
                                  devices::IS.FlattenIteratorWrapper{St}) where {St<:PSY.Storage}
    add_variable(canonical,
                 devices,
                 Symbol("Qst_$(St)"),
                 false,
                 :nodal_balance_reactive)

    return

end

function energy_storage_variables!(canonical::Canonical,
                                  devices::IS.FlattenIteratorWrapper{St}) where St<:PSY.Storage

    add_variable(canonical,
                 devices,
                 Symbol("Est_$(St)"),
                 false;
                 lb_value = d -> 0.0,)

    return

end

function storage_reservation_variables!(canonical::Canonical,
                                       devices::IS.FlattenIteratorWrapper{St}) where St<:PSY.Storage

    add_variable(canonical,
                 devices,
                 Symbol("Rst_$(St)"),
                 true)

    return

end

###################################################### output power constraints#################################

function active_power_constraints!(canonical::Canonical,
                                  devices::IS.FlattenIteratorWrapper{St},
                                  ::Type{BookKeeping},
                                  ::Type{S}) where {St<:PSY.Storage,
                                                                      S<:PM.AbstractPowerModel}

    range_data_in = [(PSY.get_name(s), PSY.get_inputactivepowerlimits(s)) for s in devices]
    range_data_out = [(PSY.get_name(s), PSY.get_outputactivepowerlimits(s)) for s in devices]

    device_range(canonical,
                 range_data_in,
                 Symbol("inputpower_range_$(St)"),
                 Symbol("Psin_$(St)"))

    device_range(canonical,
                range_data_out,
                Symbol("outputpower_range_$(St)"),
                Symbol("Psout_$(St)"))

    return

end

function active_power_constraints!(canonical::Canonical,
                                  devices::IS.FlattenIteratorWrapper{St},
                                  ::Type{BookKeepingwReservation},
                                  ::Type{S}) where {St<:PSY.Storage,
                                                                      S<:PM.AbstractPowerModel}

    range_data_in = [(PSY.get_name(s), PSY.get_inputactivepowerlimits(s)) for s in devices]
    range_data_out = [(PSY.get_name(s), PSY.get_outputactivepowerlimits(s)) for s in devices]

    reserve_device_semicontinuousrange(canonical,
                                       range_data_in,
                                       Symbol("inputpower_range_$(St)"),
                                       Symbol("Psin_$(St)"),
                                       Symbol("Rst_$(St)"))

    reserve_device_semicontinuousrange(canonical,
                                       range_data_out,
                                       Symbol("outputpower_range_$(St)"),
                                       Symbol("Psout_$(St)"),
                                       Symbol("Rst_$(St)"))

    return

end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactive_power_constraints!(canonical::Canonical,
                                   devices::IS.FlattenIteratorWrapper{St},
                                   ::Type{D},
                                   ::Type{S}) where {St<:PSY.Storage,
                                                                       D<:AbstractStorageFormulation,
                                                                       S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(s), PSY.get_reactivepowerlimits(s)) for s in devices]

    device_range(canonical,
                 range_data,
                 Symbol("reactiverange_$(St)"),
                 Symbol("Qst_$(St)"))

    return

end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(canonical::Canonical,
                            devices::IS.FlattenIteratorWrapper{St},
                            ::Type{D}) where {St<:PSY.Storage,
                                                                D<:AbstractStorageFormulation}

    storage_energy_init(canonical, devices)

return

end

###################################################### Energy Capacity constraints##########

function energy_capacity_constraints!(canonical::Canonical,
                                    devices::IS.FlattenIteratorWrapper{St},
                                    ::Type{D},
                                    ::Type{S}) where {St<:PSY.Storage,
                                                                        D<:AbstractStorageFormulation,
                                                                        S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(s), PSY.get_capacity(s)) for s in devices]

    device_range(canonical,
                 range_data,
                 Symbol("energy_capacity_$(St)"),
                 Symbol("Est_$(St)"))
    return

end

###################################################### book keeping constraints ############

function make_efficiency_data(devices::IS.FlattenIteratorWrapper{St}) where {St<:PSY.Storage}

    names = Vector{String}(undef, length(devices))
    in_out = Vector{InOut}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        names[ix] = PSY.get_name(d)
        in_out[ix] = PSY.get_efficiency(d)
    end

    return names, in_out

end

function energy_balance_constraint!(canonical::Canonical,
                                   devices::IS.FlattenIteratorWrapper{St},
                                   ::Type{D},
                                   ::Type{S}) where {St<:PSY.Storage,
                                                            D<:AbstractStorageFormulation,
                                                            S<:PM.AbstractPowerModel}

    key = ICKey(DeviceEnergy, St)

    if !(key in keys(canonical.initial_conditions))
        error("Initial Conditions for $(St) Energy Constraints not in the model")
    end

    efficiency_data = make_efficiency_data(devices)

    energy_balance(canonical,
                   canonical.initial_conditions[key],
                   efficiency_data,
                   Symbol("energy_balance_$(St)"),
                   (Symbol("Psout_$(St)"), Symbol("Psin_$(St)"), Symbol("Est_$(St)")))

    return

end


text:
include("eliminacao01.jl");
include("sr00.jl");

function resolve00(A,b)

	(m,n) = size(A);

	U = eliminacao01([A b]);
	b = U[:,n+1];
	U = U[:,1:n];

	x = sr00(U,b);

	return x;
end


text:
using ModelingToolkit
using Unitful
using LinearAlgebra
using ForwardDiff

# TODO: Check which of these are still in use and/or necessary.

# Avoid Plots compatibility issue with Unitful
Base.occursin(s::Any, q::Unitful.Quantity) = occursin(s, repr(q))
Base.pointer(q::Unitful.Quantity) = pointer(repr(q))

# Ensure that ForwardDiff values can get sent to C calls
Base.unsafe_convert(T::Type{<:Any}, x::ForwardDiff.Dual) = T(ForwardDiff.value.(x))

# Compute norms for arrays of symbolic variables (as needed by EphemerisNBP)
LinearAlgebra.norm(a::AbstractArray{<:ModelingToolkit.Num}) = sum(a .^ 2)^(1/2)


text:
using Test
using CartesianGeneticProgramming
using RDatasets
import Cambrian
import Random

cfg = get_config("../cfg/test.yaml")

function data_setup()
    iris = dataset("datasets", "iris")
    X = convert(Matrix, iris[:, 1:4])'
    X = X ./ maximum(X; dims=2)
    r = iris[:, 5].refs
    Y = zeros(maximum(r), size(X, 2))
    for i in 1:length(r)
        Y[r[i], i] = 1.0
    end
    X, Y
end

@testset "Lexicase selection Evolution" begin
    X, Y = data_setup()
    # TODO: check small example
    cfg = get_config("../cfg/iris.yaml")
    cfg["n_gen"] = 1000

    e = Cambrian.Evolution(CGPInd, cfg; id="iris")
    mutation = i::CGPInd->goldman_mutate(cfg, i)
    e.populate = x::Cambrian.Evolution->Cambrian.oneplus_populate!(
        x; mutation=mutation)
    e.evaluate = x::Cambrian.Evolution->Cambrian.lexicase_evaluate!(
        x, X, Y, CartesianGeneticProgramming.interpret)

    Cambrian.step!(e)
    @test length(e.population) == cfg["n_population"]
    best = sort(e.population)[end]
    @test best.fitness[1] >= 0.0
    @test best.fitness[1] <= size(X, 2)
    @test e.gen == 1

    Cambrian.run!(e)
    @test length(e.population) == cfg["n_population"]
    @test e.gen == cfg["n_gen"]
    println("Evolution step for lexicase selection symbolic regression")
    @timev Cambrian.step!(e)
    new_best = sort(e.population)[end]
    println("Final fitness: ", new_best.fitness[1])
    @test new_best.fitness[1] >= 0.0
end


text:
using Distributions, LinearAlgebra,Parameters

struct rkhs_reg
    alf
    sig
    x
end

function rkhs_gaussian(x::Array{<:Any,1},y,sig,lambda)
    N = size(x,1)
    #check dimensions
    if size(y,1) != N
        @error "Dimension Mismatch"
    end
    kernel_mat = zeros(N,N)
    kernel(x,z) = pdf(Normal(z,sig),x)
    for i in 1:N
        kernel_mat[:,i] = [kernel(x[j],x[i])/N for j in 1:N]
    end
    alf = pinv(kernel_mat + lambda*I(N))*y/sqrt(N)
    return rkhs_reg(alf,sig,x)
end

function rkhs_gaussian(x::Array{<:Any,2},y,sig,lambda)
    N = size(x,1)
    #check dimensions
    if size(y,1) != N
        @error "Dimension Mismatch"
    end
    kernel_mat = zeros(N,N)
    kernel(x,z) = pdf(MvNormal(z,sig),x)
    for i in 1:N
        kernel_mat[:,i] = [kernel(x[j,:],x[i,:])/N for j in 1:N]
    end
    alf = pinv(kernel_mat + lambda*I(N))*y/sqrt(N)
    return rkhs_reg(alf,sig,x)
end

function(model::rkhs_reg)(val)
    @unpack alf,sig,x = model
    N = size(x,1)
    kern_val = Array{Any}(undef,N)
    if size(x,2) == 1
        for i in 1:N
            kern_val[i] = pdf(Normal(x[i],sig),val)
        end
    else
        for i in 1:N
            kern_val[i] = pdf(MvNormal(x[i,:],sig),val)
        end
    end
    return 1/sqrt(N)*alf'*kern_val
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC