text:
using GeoJSON, GeoInterface, JSON3
using Test

include(joinpath(@__DIR__, "geojson_samples.jl"))

@testset "GeoJSON" begin
@testset "A: CRS link" begin
    feature = GeoJSON.read(a)
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) == nothing
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 2
    @test GeoInterface.bbox(feature) == nothing
    @test GeoInterface.crs(feature) == feature.properties["crs"]
    @test GeoInterface.crs(feature)["properties"]["href"] == "data.crs"
    @test GeoInterface.crs(feature)["properties"]["type"] == "ogcwkt"
    @test GeoInterface.crs(feature)["type"] == "link"
    @test "Ã" in keys(GeoInterface.properties(feature))
    @test GeoInterface.properties(feature)["Ã"] == "Ã"
    dict = geo2dict(feature)
    @test dict["geometry"] == nothing
    @test dict["type"] == "Feature"
    @test dict["crs"]["properties"]["href"] == "data.crs"
    @test dict["crs"]["properties"]["type"] == "ogcwkt"
    @test dict["crs"]["type"] == "link"
end

@testset "B: Multipoint" begin
    feature = GeoJSON.read(b)
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) isa GeoInterface.MultiPoint
    coords = GeoInterface.coordinates(GeoInterface.geometry(feature))
    @test GeoInterface.coordinates(GeoInterface.geometry(feature)) == [[-155.52,19.61],[-156.22,20.74],[-157.97,21.46]]
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 2
    @test !("Ã" in keys(feature.properties))
    @test feature.properties["type"] == "é"
    @test feature.properties["featureid"] == 1
end

@testset "C: EPSG" begin
    feature = GeoJSON.read(c)
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) == nothing
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 3
    @test feature.properties["type"] == "meow"
    @test GeoInterface.crs(feature) == feature.properties["crs"]
    @test GeoInterface.crs(feature) isa JSON3.Object
    @test GeoInterface.crs(feature)["type"] == "name"
    @test GeoInterface.crs(feature)["properties"]["name"] == "urn:ogc:def:crs:EPSG::3785"
    dict = geo2dict(feature)
    @test dict["geometry"] == nothing
    @test dict["id"] == 1
    @test dict["type"] == "Feature"
    @test dict["crs"]["properties"]["name"] == "urn:ogc:def:crs:EPSG::3785"
    @test dict["crs"]["type"] == "name"
end

@testset "D: MultiLineString" begin
    feature = GeoJSON.read(d)
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) isa GeoInterface.MultiLineString
    @test GeoInterface.geometry(feature) == feature.geometry
    @test feature.geometry.coordinates == [[[3.75,9.25],[-130.95,1.52]],[[23.15,-34.25],[-1.35,-4.65],[3.45,77.95]]]
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 3
    @test feature.properties["title"] == "Dict 1"
    @test GeoInterface.bbox(feature) == feature.properties["bbox"]
    @test GeoInterface.bbox(feature) == [-180.0, -90.0, 180.0, 90.0]
    @test GeoInterface.properties(feature)["featureid"] == "1"
end

@testset "E: HTTP Links" begin
    feature = GeoJSON.read(e)
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) isa GeoInterface.Point
    @test GeoInterface.geometry(feature) == feature.geometry
    @test feature.geometry.coordinates == [53, -4]
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 4
    @test feature.properties["title"] == "Feature 1"
    @test feature.properties["summary"] == "The first feature"
    @test feature.properties["link"] == "http://example.org/features/1"
    dict = geo2dict(feature)
    @test dict["geometry"]["coordinates"] ≈ [53.0,-4.0]
    @test dict["geometry"]["type"] == "Point"
    @test dict["id"] == "1"
    @test dict["type"] == "Feature"
end

@testset "F: Null geometry" begin
    feature = GeoJSON.read(f)
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) == nothing
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 2
    @test feature.properties["foo"] == "bar"
    @test feature.properties["featureid"] == 12
end

@testset "G: MultiPolygon" begin
    featurecollection = GeoJSON.read(g)
    # printing to GeoJSON not yet implemented
    # @test JSON3.write(featurecollection) ==
    #     "{\"features\":[{\"geometry\":{\"coordinates\":[[[[-117.913883,33.96657],[-117.907767,33.967747],[-117.912919,33.96445],[-117.913883,33.96657]]]],\"type\":\"MultiPolygon\"},\"properties\":{\"addr2\":\"Rowland Heights\",\"cartodb_id\":46,\"addr1\":\"18150 E. Pathfinder Rd.\",\"park\":\"Pathfinder Park\"},\"type\":\"Feature\"}],\"bbox\":[100.0,0.0,105.0,1.0],\"type\":\"FeatureCollection\",\"crs\":{\"properties\":{\"name\":\"urn:ogc:def:crs:EPSG::3785\"},\"type\":\"name\"}}"
    @test featurecollection.bbox ≈ [100,0,105,1]
    @test featurecollection.crs == Dict("properties" => Dict("name" => "urn:ogc:def:crs:EPSG::3785"),
                                        "type" => "name")
end

@testset "H: Print" begin
    feature = GeoJSON.read(h)
    # printing to GeoJSON not yet implemented
    # @test JSON3.write(feature) ==
    #     "{\"geometry\":{\"coordinates\":[[[3.75,9.25],[-130.95,1.52]],[[23.15,-34.25],[-1.35,-4.65],[3.45,77.95]]],\"type\":\"MultiLineString\"},\"properties\":{\"title\":\"Dict 1\"},\"bbox\":[-180.0,-90.0,180.0,90.0],\"type\":\"Feature\"}"
end

@testset "Multipolygon" begin
    testobj = GeoJSON.read(multipolygon)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    features = GeoInterface.features(testobj)
    @test length(features) == 1
    feature = features[1]
    @test GeoInterface.bbox(feature) == feature.properties["bbox"]
    @test GeoInterface.bbox(feature) == [100.0, 0.0, 105.0, 1.0]
    @test GeoInterface.geometry(feature) isa GeoInterface.MultiPolygon
    @test GeoInterface.geometry(feature) == feature.geometry
    @test size(feature.geometry.coordinates) == (1,)
    @test length(feature.geometry.coordinates) == 1
    @test length(feature.geometry.coordinates[1]) == 1
    @test length(feature.geometry.coordinates[1][1]) == 4
    @test feature.geometry.coordinates[1][1][1] ≈ [-117.914, 33.9666] atol=1e-3
    @test feature.geometry.coordinates[1][1][2] ≈ [-117.908,33.9677] atol=1e-3
    @test feature.geometry.coordinates[1][1][3] ≈ [-117.913,33.9644] atol=1e-3
    @test feature.geometry.coordinates[1][1][4] ==
        feature.geometry.coordinates[1][1][4]
    @test feature.properties["addr2"] == "Rowland Heights"
    @test feature.properties["cartodb_id"] == 46
    @test feature.properties["addr1"] == "18150 E. Pathfinder Rd."
    @test feature.properties["park"] == "Pathfinder Park"
end

@testset "Realmultipolygon" begin
    testobj = GeoJSON.read(realmultipolygon)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    features = GeoInterface.features(testobj)
    @test length(features) == 1
    feature = features[1]
    @test GeoInterface.geometry(feature) isa GeoInterface.MultiPolygon
    @test GeoInterface.geometry(feature) == feature.geometry
    @test size(feature.geometry.coordinates) == (2,)
    @test feature.geometry.coordinates == [[[[102.0,2.0],[103.0,2.0],[103.0,3.0],[102.0,3.0],[102.0,2.0]]],[[[100.0,0.0],[101.0,0.0],[101.0,1.0],[100.0,1.0],[100.0,0.0]],[[100.2,0.2],[100.8,0.2],[100.8,0.8],[100.2,0.8],[100.2,0.2]]]]
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1]["geometry"]["type"] == "MultiPolygon"
    dict = geo2dict(feature)
    @test dict["type"] == "Feature"
    @test length(dict["geometry"]["coordinates"]) == 2
end

@testset "Polyline" begin
    testobj = GeoJSON.read(polyline)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    features = GeoInterface.features(testobj)
    @test length(features) == 1
    feature = features[1]
    @test feature.geometry isa GeoInterface.LineString
    @test feature.properties["InLine_FID"] == 0
    @test feature.properties["SimLnFLag"] == 0
    @test feature.geometry.coordinates == [[-89.0,43.0],[-88.0,44.0],[-88.0,45.0]]
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1]["geometry"]["type"] == "LineString"
    dict = geo2dict(feature)
    @test dict["type"] == "Feature"
    @test dict["id"] == "a73ws67n775q"
    @test length(dict["geometry"]["coordinates"]) == 3
end

@testset "Point" begin
    testobj = GeoJSON.read(point)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    feature = testobj.features[1]
    @test feature.properties["fax"] == "305-571-8347"
    @test feature.properties["phone"] == "305-571-8345"
    @test feature.geometry isa GeoInterface.Point
    @test feature.geometry.coordinates == [-89.0, 44.0]
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1]["geometry"]["type"] == "Point"
    dict = geo2dict(feature)
    @test dict["type"] == "Feature"
    @test dict["id"] == "a7vs0i9rnyyx"
    @test length(dict["geometry"]["coordinates"]) == 2
end

@testset "Pointnull" begin
    testobj = GeoJSON.read(pointnull)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    feature = testobj.features[1]
    @test feature.properties["featureid"] == "a7vs0i9rnyyx"
    @test feature.properties["fax"] == "305-571-8347"
    @test feature.properties["phone"] == "305-571-8345"
    @test feature.geometry isa Nothing
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1]["geometry"] == nothing
    dict = geo2dict(feature)
    @test dict["geometry"] == nothing
    @test dict["id"] == "a7vs0i9rnyyx"
    @test dict["type"] == "Feature"
end

@testset "Poly" begin
    testobj = GeoJSON.read(poly)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    feature = testobj.features[1]
    @test feature.properties["featureid"] == "a7ws7wldxold"
    @test feature.properties["DIST_NUM"] == 7.0
    @test feature.properties["PHONE"] == "686-3070"
    @test feature.properties["AREA_SQMI"] == 12.41643
    @test feature.properties["LOCATION"]  == "Bustleton Ave. & Bowler St"
    @test feature.properties["DIV_CODE"]  == "NEPD"
    @test feature.properties["DIST_NUMC"] == "07"
    @test feature.geometry isa GeoInterface.Polygon
    @test length(feature.geometry.coordinates) == 1
    @test length(feature.geometry.coordinates[1]) == 5
    @test feature.geometry.coordinates[1][1] ==
        feature.geometry.coordinates[1][end]
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1] == geo2dict(feature)
    dict = geo2dict(feature)
    @test dict["geometry"]["type"] == "Polygon"
    @test dict["id"] == "a7ws7wldxold"
    @test dict["type"] == "Feature"
end

@testset "Polyhole" begin
    testobj = GeoJSON.read(polyhole)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    feature = testobj.features[1]
    @test feature.properties["featureid"] == "a7ws7wldxold"
    @test feature.properties["DIST_NUM"] == 7.0
    @test feature.properties["PHONE"] == "686-3070"
    @test feature.properties["AREA_SQMI"] == 12.41643
    @test feature.properties["LOCATION"]  == "Bustleton Ave. & Bowler St"
    @test feature.properties["DIV_CODE"]  == "NEPD"
    @test feature.properties["DIST_NUMC"] == "07"
    @test feature.geometry isa GeoInterface.Polygon
    @test length(feature.geometry.coordinates) == 2
    @test length(feature.geometry.coordinates[1]) == 5
    @test feature.geometry.coordinates[1][1] ==
        feature.geometry.coordinates[1][end]
    @test length(feature.geometry.coordinates[2]) == 5
    @test feature.geometry.coordinates[2][1] ==
        feature.geometry.coordinates[2][end]
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1] == geo2dict(feature)
    dict = geo2dict(feature)
    @test dict["geometry"]["type"] == "Polygon"
    @test dict["id"] == "a7ws7wldxold"
    @test dict["type"] == "Feature"
end

@testset "Collection" begin
    testobj = GeoJSON.read(collection)
    @test testobj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    feature = testobj.features[1]
    @test feature.properties["STATE_ABBR"] == "ZZ"
    @test feature.properties["STATE_NAME"] == "Top"
    @test feature.geometry isa GeoInterface.GeometryCollection
    @test length(feature.geometry.geometries) == 3
    @test feature.geometry.geometries[1] isa GeoInterface.Polygon
    @test feature.geometry.geometries[2] isa GeoInterface.Polygon
    @test feature.geometry.geometries[3] isa GeoInterface.Point
    coords = feature.geometry.geometries[1].coordinates
    @test length(coords) == 1
    @test length(coords[1]) == 5
    @test coords[1][1] == coords[1][end]
    coords = feature.geometry.geometries[2].coordinates
    @test length(coords) == 1
    @test length(coords[1]) == 5
    @test coords[1][1] == coords[1][end]
    @test feature.geometry.geometries[3].coordinates ≈ [-94.0,46.0]
    dict = geo2dict(testobj)
    @test dict["type"] == "FeatureCollection"
    @test dict["features"][1] == geo2dict(feature)
    dict = geo2dict(feature)
    @test dict["geometry"]["type"] == "GeometryCollection"
    @test length(dict["geometry"]["geometries"]) == 3
    @test dict["geometry"]["geometries"][1]["type"] == "Polygon"
    @test dict["geometry"]["geometries"][2]["type"] == "Polygon"
    @test dict["geometry"]["geometries"][3]["type"] == "Point"
    @test dict["id"] == "a7xlmuwyjioy"
    @test dict["type"] == "Feature"
end

@testset "OSM buildings" begin
    buildings = GeoJSON.read(osm_buildings)
    @test buildings isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    @test length(buildings.features) == 4
    @test map(typeof,buildings.features) == fill(GeoInterface.Feature,4)

    feature = buildings.features[1]
    @test feature isa GeoInterface.Feature
    @test GeoInterface.geometry(feature) isa GeoInterface.Polygon
    @test GeoInterface.properties(feature) == feature.properties
    @test feature.properties isa Dict{String,Any}
    @test length(feature.properties) == 2
    @test feature.properties["height"] == 150
    @test feature.properties["color"] == "rgb(255,200,150)"

    building_dict = geo2dict(buildings)
    @test building_dict isa Dict{String,Any}

    # printing to GeoJSON not yet implemented
    # @test GeoJSON.write(buildings) ==
    # """{\"features\":[{\"geometry\":{\"coordinates\":[[[13.42634,52.49533],[13.4266,52.49524],[13.42619,52.49483],[13.42583,52.49495],[13.4259,52.49501],[13.42611,52.49494],[13.4264,52.49525],[13.4263,52.49529],[13.42634,52.49533]]],\"type\":\"Polygon\"},\"properties\":{\"height\":150,\"color\":\"rgb(255,200,150)\"},\"type\":\"Feature\"},{\"geometry\":{\"coordinates\":[[[13.42706,52.49535],[13.42745,52.4952],[13.42745,52.4952],[13.42741,52.49516],[13.42717,52.49525],[13.42692,52.49501],[13.42714,52.49494],[13.42686,52.49466],[13.4265,52.49478],[13.42657,52.49486],[13.42678,52.4948],[13.42694,52.49496],[13.42675,52.49503],[13.42706,52.49535]]],\"type\":\"Polygon\"},\"properties\":{\"height\":130,\"color\":\"rgb(180,240,180)\"},\"type\":\"Feature\"},{\"geometry\":{\"coordinates\":[[[[13.42746,52.4944],[13.42794,52.49494],[13.42799,52.49492],[13.42755,52.49442],[13.42798,52.49428],[13.42846,52.4948],[13.42851,52.49478],[13.428,52.49422],[13.42746,52.4944]]],[[[13.42803,52.49497],[13.428,52.49493],[13.42844,52.49479],[13.42847,52.49483],[13.42803,52.49497]]]],\"type\":\"MultiPolygon\"},\"properties\":{\"height\":120,\"color\":\"rgb(200,200,250)\"},\"type\":\"Feature\"},{\"geometry\":{\"coordinates\":[[[13.42857,52.4948],[13.42918,52.49465],[13.42867,52.49412],[13.4285,52.49419],[13.42896,52.49465],[13.42882,52.49469],[13.42837,52.49423],[13.42821,52.49428],[13.42863,52.49473],[13.42853,52.49476],[13.42857,52.4948]]],\"type\":\"Polygon\"},\"properties\":{\"height\":140,\"color\":\"rgb(150,180,210)\"},\"type\":\"Feature\"}],\"type\":\"FeatureCollection\"}"""
    # @test GeoJSON.write(buildings) == JSON3.write(JSON3.read(osm_buildings))
    # @test building_dict == JSON3.read(osm_buildings)
end

@testset "Tech Square: read file" begin
    path = joinpath(@__DIR__, "tech_square.geojson")
    obj = GeoJSON.read(read(path))
    @test obj isa GeoInterface.FeatureCollection{GeoInterface.Feature}
    @test length(GeoInterface.features(obj)) == 171
    feature = obj.features[1]
    @test feature.properties["featureid"] == "relation/2119819"
    @test feature.geometry isa GeoInterface.Polygon
    @test length(feature.geometry.coordinates) == 1
    @test length(feature.geometry.coordinates[1]) == 38
    @test feature.geometry.coordinates[1][1] == feature.geometry.coordinates[1][end]
end

@testset "Deprecations" begin
    path = joinpath(@__DIR__, "tech_square.geojson")
    @test (@test_deprecated repr(GeoJSON.parsefile(path))) == repr(GeoJSON.read(read(path)))
    dict = GeoJSON.read(g)
    @test (@test_deprecated repr(GeoJSON.parse(g))) == repr(dict)
    @test (@test_deprecated GeoJSON.geojson(dict)) == GeoJSON.write(dict)
end

end # testset "GeoJSON"


text:
# main file for evalResidual_revm

import PDESolver.evalResidual_revm

@doc """
### EulerEquationMod.evalrevm_transposeproduct

Reverse mode of evalResidual with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""->

function evalResidual_revm(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData,
                     opts::Dict, t::Number=0.0)

  time = eqn.params.time
  eqn.params.t = t  # record t to params
  myrank = mesh.myrank

  # Forward sweep
  time.t_dataprep += @elapsed dataPrep_for_revm(mesh, sbp, eqn, opts)

  time.t_volume += @elapsed if opts["addVolumeIntegrals"]
    evalVolumeIntegrals_revm(mesh, sbp, eqn, opts)
  end

  if opts["use_GLS"]
    error("GLS not supported for revm product")
  end

  time.t_bndry += @elapsed if opts["addBoundaryIntegrals"]
    evalBoundaryIntegrals_revm(mesh, sbp, eqn, opts)
  end

   time.t_stab += @elapsed if opts["addStabilization"]
     addStabilization_revm(mesh, sbp, eqn, opts)
   end

  time.t_face += @elapsed if mesh.isDG && opts["addFaceIntegrals"]
    evalFaceIntegrals_revm(mesh, sbp, eqn, opts)
  end

  time.t_sharedface += @elapsed if mesh.commsize > 1
    evalSharedFaceIntegrals_revm(mesh, sbp, eqn, opts)
  end

  time.t_shock += @elapsed if opts["addShockCapturing"]
    evalShockCapturing_revm(mesh, sbp, eqn, opts)
  end

  time.t_source += @elapsed evalSourceTerm_revm(mesh, sbp, eqn, opts)

  # apply inverse mass matrix to eqn.res, necessary for CN
  if opts["use_Minv"]
    error("use_Minv not supported for revm product")
  end

  time.t_dataprep += @elapsed dataPrep_revm(mesh, sbp, eqn, opts)

  return nothing
end  # end evalResidual

"""
  Dataprep-type function that should be called at the beginning of a reverse
  mode product, rather than the regular `dataPrep` function.
"""
function dataPrep_for_revm(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                   eqn::AbstractEulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  # apply filtering to input
  if eqn.params.use_filter
    applyFilter(mesh, sbp, eqn, eqn.q, opts)
  end

  getAuxVars(mesh, eqn)

  if opts["use_staggered_grid"]
    error("staggered grid not supported for revm product")
    aux_vars = zeros(Tres, 1, mesh.mesh2.numNodesPerElement)
    for i=1:mesh.numEl
      qs = ro_sview(eqn.q, :, :, i)
      qf = sview(eqn.q_flux, :, :, i)
      interpolateElementStaggered(eqn.params, mesh, qs, aux_vars, qf)
    end
  end

  # only eqn.q_bndry is required, all other methods are no-precompute
  if mesh.isDG
    interpolateBoundary(mesh, sbp, eqn, opts, eqn.q, eqn.q_bndry, eqn.aux_vars_bndry)
  end

  if eqn.params.use_edgestab
    stabscale(mesh, sbp, eqn)
  end

  return nothing
end

@doc """
### EulerEquationMod.dataPrep_revm

Reverse mode of dataPrep w.r.t mesh metrics

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""->
function dataPrep_revm(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                   eqn::AbstractEulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  # nothing to do here

  return nothing
end

@doc """
### EulerEquationMod.evalVolumeIntegrals_revm

Reverse mode of evalVolumeIntegrals with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""->
function evalVolumeIntegrals_revm(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh,  Tsol, Tres, Tdim}

  integral_type = opts["volume_integral_type"]
  if integral_type == 1
    if opts["Q_transpose"] == true
      calcVolumeIntegrals_nopre_revm(mesh, sbp, eqn, opts)
    else
      error("Q_transpose = false not supported for revm product")
    end  # end if
  elseif integral_type == 2
    calcVolumeIntegralsSplitForm_revm(mesh, sbp, eqn, opts, eqn.volume_flux_func,
                                 eqn.volume_flux_func_revm)
  else
    throw(ErrorException("Unsupported volume integral type = $integral_type"))
  end

  return nothing
end  # end evalVolumeIntegrals

@doc """
### EulerEquationMod.evalBoundaryIntegrals_revm

Reverse mode of evalBoundaryIntegrals with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""->

function evalBoundaryIntegrals_revm(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh, Tsol, Tres, Tdim}

  getBCFluxes_revm(mesh, sbp, eqn, opts)

  return nothing

end  # end evalBoundaryIntegrals

function addStabilization_revm(mesh::AbstractMesh, sbp::AbstractOperator,
                               eqn::EulerData, opts)

  if opts["use_lps"]
    applyLPStab_revm(mesh, sbp, eqn, opts)
  end

  return nothing
end

@doc """
### EulerEquationMod.evalFaceIntegrals_revm

Reverse mode of evalFaceIntegrals with respect to the mesh metrics ∂ξ/∂x

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""->

function evalFaceIntegrals_revm(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator, eqn::EulerData{Tsol}, opts) where {Tmsh, Tsol}

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1
    calcFaceIntegral_nopre_revm(mesh, sbp, eqn, opts, eqn.flux_func_revm,
                         mesh.interfaces)

  elseif face_integral_type == 2

    getFaceElementIntegral_revm(mesh, sbp, eqn, eqn.face_element_integral_func,
                           eqn.flux_func, mesh.sbpface, mesh.interfaces)

  else
    throw(ErrorException("Unsupported face integral type = $face_integral_type"))
  end

  # do some output here?
  return nothing
end

@doc """

Reverse mode evalSharedFaceIntegrals with respect to the mesh metrics

"""

function evalSharedFaceIntegrals_revm(mesh::AbstractDGMesh, sbp, eqn, opts)

  if getParallelData(eqn.shared_data) != PARALLEL_DATA_ELEMENT
    error("""parallel data setting must be PARALLEL_DATA_ELEMENT """)
  end

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1

    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceIntegrals_element_revm)

  elseif face_integral_type == 2

    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceElementIntegrals_element_revm)
#    getSharedFaceElementIntegrals_element(mesh, sbp, eqn, opts, eqn.face_element_integral_func,  eqn.flux_func)
  else
    throw(ErrorException("unsupported face integral type = $face_integral_type"))
  end

  return nothing
end

#=
#TODO: what is this and why is it here?
"""
### EulerEquationMod.calcSharedFaceIntegrals_revm

Reverse mode of calcSharedFaceIntegrals w.r.t mesh metrics, ∂ξ/∂x.

"""

function calcSharedFaceIntegrals_revm( mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol},
                            opts, functor_revm::FluxType) where {Tmsh, Tsol}
  # calculate the face flux and do the integration for the shared interfaces

  #TODO: update this to use the new parallel communication primatives
  if opts["parallel_data"] != PARALLEL_DATA_FACE
    throw(ErrorException("cannot use calcSharedFaceIntegrals without parallel face data"))
  end

  params = eqn.params

  npeers = mesh.npeers
  val = sum(mesh.recv_waited)
  if val !=  mesh.npeers && val != 0
    throw(ErrorException("Receive waits in inconsistent state: $val / $npeers already waited on"))
  end

  for i=1:mesh.npeers
    if val == 0
      params.time.t_wait += @elapsed idx, stat = MPI.Waitany!(mesh.recv_reqs)
      mesh.recv_stats[idx] = stat
      mesh.recv_reqs[idx] = MPI.REQUEST_NULL  # don't use this request again
      mesh.recv_waited[idx] = true
    else
      idx = i
    end

    # calculate the flux
    interfaces = mesh.shared_interfaces[idx]
    qL_arr = eqn.q_face_send[idx]
    qR_arr = eqn.q_face_recv[idx]
    aux_vars_arr = eqn.aux_vars_sharedface[idx]
#    dxidx_arr = mesh.dxidx_sharedface[idx]
    nrm_arr = mesh.nrm_sharedface[idx]
    flux_arr_bar = eqn.flux_sharedface_bar[idx]

    # permute the received nodes to be in the elementR orientation
    permuteinterface!(mesh.sbpface, interfaces, qR_arr)
    for j=1:length(interfaces)
      interface_i = interfaces[j]
      for k=1:mesh.numNodesPerFace
        eL = interface_i.elementL
        fL = interface_i.faceL

        qL = sview(qL_arr, :, k, j)
        qR = sview(qR_arr, :, k, j)
        aux_vars = sview(aux_vars_arr, :, k, j)
        nrm_xy = sview(nrm_arr, :, k, j)
        flux_j = sview(flux_arr_bar, :, k, j)
        functor_revm(params, qL, qR, aux_vars, nrm_xy, flux_j)
      end
    end
    # end flux calculation

    # do the integration
    boundaryintegrate_rev!(mesh.sbpface, mesh.bndries_local[idx], flux_arr_bar, eqn.res_bar, SummationByParts.Subtract())
  end  # end loop over npeers

  @debug1 sharedFaceLogging(mesh, sbp, eqn, opts, eqn.q_face_send, eqn.q_face_recv)

  return nothing
end
=#

"""
### EulerEquationMod.evalSourceTerm_revm

Reverse mode of evalSourceTerm w.r.t mesh metrics.

*This function has not been written properly. The developer must update this
documentation as and when this code is developed*

"""
function evalSourceTerm_revm(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  # placeholder for multiple source term functionality (similar to how
  # boundary conditions are done)
  if opts["use_src_term"]
    error("source terms not supported for revm product")
  end

  return nothing
end  # end function

"""

  Reverse mode of `evalShockCapturing` wrt metrics

  **Inputs**:

   * mesh : Abstract mesh type
   * sbp  : Summation-by-parts operator
   * eqn  : Euler equation object
   * opts : options dictonary
"""
function evalShockCapturing_revm(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  # currently there is only one choice for the shock capturing scheme.
  # If there are more, need something more sophisiticated to choose.
  # Perhaps add an abstract typed field to eqn containing the structs

  capture = eqn.shock_capturing
  applyShockCapturing_revm(mesh, sbp, eqn, opts, capture)

  return nothing
end

#=
@doc """
###EulerEquationMod.updatePumiMesh

Given the volume nodes in the MeshMovement `volNodes` data structure, update the
Pumi mesh object

* mesh  : a mesh object
* sbp   : SBP operator object
* volNodes : A 2D array containing coordinates of all the mesh vertices/nodes
             owned by a perticular rank. There are no repeated vertices, i.e.
             all the coordinates are unique. size(volNodes) = (3, numVert)

"""->

function updatePumiMesh(mesh::AbstractDGMesh{Tmsh}, sbp::AbstractOperator,
                        volNodes::AbstractArray{Tmsh, 2}) where Tmsh

  for i = 1:mesh.numEl
    for j = 1:size(mesh.vert_coords,2)
      # Get the vertex numbering on the portion of mesh owned by the processor
      local_vertnum = mesh.element_vertnums[j,i]
      for k = 1:mesh.dim
        mesh.vert_coords[k,j,i] = volNodes[k,local_vertnum]
      end
    end
  end

  for i = 1:mesh.numEl
    update_coords(mesh, i, mesh.vert_coords[:,:,i])
  end
  commit_coords(mesh, sbp)

  return nothing
end
=#

ext: jl


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Actors, Test, .Threads
import Actors: spawn, newLink

t = Ref{Task}()
a = Ref{Int}()
a[] = 1

inca(a, b) = a[] = a[] + b

me = newLink()

act = spawn(inca, a, taskref=t)

@test t[].state == :runnable
A = request(act, Actors.Diag, :act)
@test A.mode == :default

send(act, 1)
sleep(0.1)
@test a[] == 2
@test request(act, 1) == 3
@test a[] == 3
become!(act, threadid)
@test request(act) in 1:nthreads()
@test A.bhv == threadid
become!(act, (lk, x, y) -> send(lk, x^y))
send(act, me, 123, 456)
@test receive(me) == 2409344748064316129

act1 = spawn(threadid, sticky=true)
@test request(act1) == 1
act2 = spawn(threadid, thrd=2)
@test request(act2) == 2

println(info(act2))


text:
using Distributions

"""
    bubbleKernel(x::Float64, r::Float64)::Float64

Return 1.0 if dist <= r, otherwise 0.0.
"""
function bubbleKernel(x::Float64, r::Float64)::Float64
    return x <= r ? 1.0 : 0.0
end

"""
    gaussianKernel(x::Float64, r::Float64)::Float64

Return Gaussian(x) for μ=0.0 and σ = r/3.
(a value of σ = r/3 makes the training results comparable between different kernels
for same values of r).
"""
function gaussianKernel(x::Float64, r::Float64)::Float64
    return Distributions.pdf.(Distributions.Normal(0.0,r/3), x)
end


text:
function get_grid_bots(base_path::String=Pinguu.base_dir, api_path::String=ENV["api"])
    config_files = glob(joinpath(base_path, "grid_bots/user/*/*/config.json"))
    api = JSON.parsefile(api_path)
    telegram = JSON.parsefile(joinpath(base_path, "telegram.json"))
    dev_settings = JSON.parsefile(joinpath(base_path, "dev.json"))
    config_paths = []
    clients = []
    user = []
    for file in config_files
        tmp = JSON.parsefile(file)
        if !(tmp["status"] in ["dead", "finished"])
            bot_path = joinpath("/", dirname(file))
            user_id = string(tmp["user_id"])
            if haskey(api, user_id)
                push!(clients, binance.Client(api[user_id]["api"], api[user_id]["api_secret"]))
                push!(config_paths, bot_path)
            end
            tmp_usr = Dict()
            tmp_usr["user_id"] = user_id
            if haskey(telegram, user_id)
                tmp_usr["telegram"] = telegram[user_id]
            else
                tmp_usr["telegram"] = ""
            end
            push!(user, tmp_usr)
        end
    end
    return config_paths, clients, user
end
export get_grid_bots


text:
export i_magtensor!, i_gvb_tensor!, i_bondtensor!, i_g16_tensor!
export i_magtensor, i_gvb_tensor, i_bondtensor, i_g16_tensor
export apply_Ghb!, apply_Gvb!, apply_G16!, apply_Gh!, apply_Gcp!, apply_Greset!, apply_Gcut!

(f::Type{Inv{TropicalMatrixBlock}})(m::TropicalMatrixBlock) = m.mat

"""
    i_magtensor!(v, Jij)

`v` should be a length 2 vector.
"""
@i function i_magtensor!(v::AbstractVector{TT}, h::Real) where TT<:TropicalTypes
    @safe @assert size(v) == (2,)
    Tropical(h)
    v[1] *= h
    v[2] /= h
    (~Tropical)(h)
end

"""
    i_bondtensor!(mat, Jij)

`mat` should be a `one` tensor.
"""
@i function i_bondtensor!(mat::AbstractMatrix{TT}, Jij::Real) where TT<:TropicalTypes
    @safe @assert size(mat) == (2,2)
    Tropical(Jij)
    mat[1,1] *= Jij
    mat[2,2] *= Jij
    mat[1,2] /= Jij
    mat[2,1] /= Jij
    (~Tropical)(Jij)
end

@i function i_gvb_tensor!(mat::Diagonal{TT}, Jij::Real) where TT<:TropicalTypes
    @safe @assert size(mat) == (4,4)
    Tropical(Jij)
    mat[1,1] *= Jij
    mat[2,2] /= Jij
    mat[3,3] /= Jij
    mat[4,4] *= Jij
    (~Tropical)(Jij)
end

@i function i_g16_tensor!(out!::AbstractMatrix{T}, Js) where T<:TropicalTypes
    @safe @assert length(Js) == 16
    @safe @assert size(out!) == (16, 16)
    @routine begin
        y ← ones(T,2,2,2,2,2,2,2,2)
        xs ← ([ones(T,2,2) for i=1:16]...,)
        for i = 1:length(Js)
            i_bondtensor!(xs |> tget(i), Js[i])
        end
        ixs ← ([(ix...,) for ix in split("aα,aβ,aγ,aδ,bα,bβ,bγ,bδ,cα,cβ,cγ,cδ,dα,dβ,dγ,dδ", ',')]...,)
        iy ← ("abcdαβγδ"...,)
        LogLikeNumbers.i_einsum!(ixs, xs, iy, y)
    end
    for i=1:length(out!)
        out![i] *= y[i]
    end
    ~@routine
end

i_bondtensor(::Type{TT}, Jij::T) where {T<:Real, TT} = i_bondtensor!(ones(TT, 2, 2), Jij)[1]
i_magtensor(::Type{TT}, h::T) where {T<:Real, TT} = i_magtensor!(ones(TT, 2), h)[1]
i_gvb_tensor(::Type{TT}, Jij::T) where {T<:Real, TT} = i_gvb_tensor!(Diagonal(ones(TT, 4)), Jij)[1]
i_g16_tensor(::Type{TT}, Js::Vector{T}) where {T<:Real, TT} = i_g16_tensor!(ones(TT, 16, 16), Js)[1]

i_bondtensor(Jij::T) where T<:Real = i_bondtensor(Tropical{T}, Jij)
i_magtensor(h::T) where T<:Real = i_magtensor(Tropical{T}, h)
i_gvb_tensor(Jij::T) where T<:Real = i_gvb_tensor(Tropical{T}, Jij)
i_g16_tensor(Js::Vector{T}) where T<:Real = i_g16_tensor(Tropical{T}, Js)

@i function YaoBlocks.apply!(reg::ArrayReg{1}, pb::PutBlock{N,C,<:TropicalMatrixBlock}, REG_STACK) where {N, C}
    i_instruct!(reg.state |> vec, pb.content.mat, pb.locs, (), (), REG_STACK)
end
@i function YaoBlocks.apply!(reg::ArrayReg{1}, cb::ControlBlock{N,<:TropicalMatrixBlock}, REG_STACK) where {N, C}
    i_instruct!(reg.state |> vec, cb.content.mat, cb.locs, cb.ctrl_locs, cb.ctrl_config, REG_STACK)
end
@i function YaoBlocks.apply!(reg::ArrayReg{1,<:Tropical}, cb::ControlBlock{N,<:XGate}, REG_STACK) where {N, C}
    i_instruct!(reg.state |> vec, Val(:X), cb.locs, cb.ctrl_locs, cb.ctrl_config, REG_STACK)
end
YaoBlocks.apply!(reg::ArrayReg{B}, b::TropicalMatrixBlock, REG_STACK) where B = throw(NotImplementedError(:apply!, typeof((reg, b))))

"""
    apply_Ghb!(reg::ArrayReg, i::Int, J::Real, REG_STACK)

Apply a horizontal  spinglassbond tensor (parametrized by coupling `J`) on site `i`.
This instruct will increase the stack top of `REG_STACK` by 1.
"""
@i function apply_Ghb!(reg::ArrayReg{1,T}, i::Int, J::Real, REG_STACK) where T<:Tropical
    @routine @invcheckoff begin
        nbit ← nqubits(reg)
        blk ← put(nbit, i=>tropicalblock(MMatrix{2,2}(ones(T, 2, 2))))
        i_bondtensor!(blk.content.mat, J)
    end
    apply!(reg, blk, REG_STACK)
    ~@routine
end

"""
    apply_Gh!(reg::ArrayReg, i::Int, h::Real, REG_STACK)

Apply a spinglass magnetic tensor (parametrized by coupling `h`) on site `i`.
This instruct will increase the stack top of `REG_STACK` by 0.
"""
@i function apply_Gh!(reg::ArrayReg{1,T}, i::Int, h::Real, REG_STACK) where T<:Tropical
    @routine @invcheckoff begin
        nbit ← nqubits(reg)
        blk ← put(nbit, i=>tropicalblock(Diagonal(MVector{2}(ones(T, 2)))))
        i_magtensor!(blk.content.mat.diag, h)
    end
    apply!(reg, blk, REG_STACK)
    ~@routine
end

"""
    apply_Gvb!(reg::ArrayReg, i::NTuple{2,Int}, J::Real, REG_STACK)

Apply a vertical spinglass bond tensor (parametrized by coupling `J`) on sites `i[1]` and `i[2]`.
This instruct will increase the stack top of `REG_STACK` by 0.
"""
@i function apply_Gvb!(reg::ArrayReg{1,T}, i::NTuple{2,Int}, J::Real, REG_STACK) where T<:Tropical
    @routine @invcheckoff begin
        nbit ← nqubits(reg)
        blk ← put(nbit, i=>tropicalblock(Diagonal(MVector{4}(ones(T, 4)))))
        i_gvb_tensor!(blk.content.mat, J)
    end
    apply!(reg, blk, REG_STACK)
    ~@routine
end

"""
    apply_G16!(reg::ArrayReg, i::NTuple{4,Int}, Js::AbstractVector, REG_STACK)

Apply a spin glass tensor of K(4,4) graph (the intra-block coupling term in the Chimera graph that parametrized by 16 coupling terms `Js`) on sites `i[1:4]`.
This instruct will increase the stack top of `REG_STACK` by 1.
"""
@i function apply_G16!(reg::ArrayReg{1,T}, i::NTuple{4,Int}, Js::AbstractVector{<:Real}, REG_STACK) where T<:Tropical
    @routine @invcheckoff begin
        nbit ← nqubits(reg)
        blk ← put(nbit, i=>tropicalblock(MMatrix{16,16}(ones(T, 16, 16))))
        i_g16_tensor!(blk.content.mat, Js)
    end
    apply!(reg, blk, REG_STACK)
    ~@routine
end

"""
    apply_Gcp!(reg::ArrayReg, i::NTuple{2,Int}, REG_STACK)

Apply a copy gate (or CNOT).
This instruct will increase the stack top of `REG_STACK` by 1.
"""
@i function apply_Gcp!(reg::ArrayReg{1,T}, i::NTuple{2,Int}, REG_STACK) where T<:Tropical
    @routine @invcheckoff g ← put(nqubits(reg), i=>tropicalblock(Diagonal([one(T), T(-999999), T(-999999), one(T)])))
    apply!(reg, g, REG_STACK)
    ~@routine
end

"""
    apply_Greset!(reg::ArrayReg, i::Int, REG_STACK)

Apply a reset gate (or SUM).
This instruct will increase the stack top of `REG_STACK` by 1.
"""
@i function apply_Greset!(reg::ArrayReg{1,T}, i::Int, REG_STACK) where T<:Tropical
    @routine @invcheckoff g ← put(nqubits(reg), i=>tropicalblock(MMatrix{2,2}(resettensor(T))))
    apply!(reg, g, REG_STACK)
    ~@routine
end

"""
    apply_Gcut!(reg::ArrayReg, i::Int, REG_STACK)

Apply a cut gate.
This instruct will increase the stack top of `REG_STACK` by 1.
"""
@i function apply_Gcut!(reg::ArrayReg{1,T}, i::Int, REG_STACK) where T<:Tropical
    @routine @invcheckoff g ← put(nqubits(reg), i=>tropicalblock(MMatrix{2,2}(one(T), one(T), one(T), one(T))))
    apply!(reg, g, REG_STACK)
    ~@routine
end


text:
module ComponentPackages
    import JSON, JSON2, HTTP

    include("ComponentMetas.jl")
    using .ComponentMetas

    const ROOT_PATH = (@__DIR__) * "/.."
    const META_FILENAME = ROOT_PATH * "/components_meta.json"

    struct ComponentPackage
        name ::String
        package_name ::String
        source_path ::String
        scripts ::NamedTuple{(:dev, :prod), Tuple{Vector{String}, Vector{String}}}
        has_components ::Bool
    end

    macro _make_packages()
        raw_meta = open(META_FILENAME, "r") do f
            JSON2.read(f, Vector{ComponentPackage})
        end

        package_pairs = map(raw_meta) do p
            return :($(Symbol(p.name)) = $(p))
        end
        return esc(:($(package_pairs...),))
    end

    const _components_packages = @_make_packages

    function components_js_include(prefix = ""; debug=false)
        type = debug ? :dev : :prod
        join(
            map(values(_components_packages)) do p
                join(
                    map(p.scripts[type]) do script
                        """<script src="$(prefix)_dash-component-suites/$(p.package_name)/$(script)"></script>"""
                    end, "\n"
                )
            end
        )
    end

    macro register_js_sources(path, prefix)
        registers = []
        for package in values(_components_packages)
            package_uri = "_dash-component-suites/$(package.package_name)/"
            package_path = ROOT_PATH * package.source_path
            part = quote
                full_package_uri = $(prefix) * $(package_uri)
                if startswith($(path), full_package_uri)
                    script_name = replace($(path), full_package_uri=>"")
                    filename = $(package_path) * script_name
                    try
                        return HTTP.Response(200, ["Content-Type" => "application/javascript"], body = read(filename))
                    catch
                        return HTTP.Response(404)
                    end
                end
            end
            push!(registers, part)
        end
        return esc(quote
            $(registers...)
        end)
    end

    function package_namespace(package_name)
        return _components_packages[Symbol(package_name)].package_name
    end

    function load_package_components(package_name)
        package_symbol = Symbol(package_name)
        if !haskey(_components_packages, package_symbol)
            error("undefined package $(package_symbol)")
        end
        package = _components_packages[Symbol(package_name)]

        meta_filename = "$(ROOT_PATH)/$(package.source_path)metadata.json"
        return load_components_meta(meta_filename)
        raw_meta = JSON.parsefile(meta_filename, dicttype = OrderedDict{String, Any})
        components = Dict{Symbol, ComponentMeta}()
        for (script, meta) in raw_meta
            name = match(r"([a-zA-Z0-9_]+)\.react\.js", script).captures[1]
            components[Symbol(name)] = ComponentMeta(name, meta["description"], parse_props(meta["props"]))
            break
        end
        return components
    end

    function component_doc_list()
        result = ""
        for package in values(_components_packages)
            if package.has_components
                result *= "## Component package: `$(package.name)`:\n"
                components = load_package_components(Symbol(package.name))
                names = sort(collect(keys(components)))
                names = map(names) do name
                    lowercase("`$(package.name)_$(name)`")
                end
                result *= join(names, ", ") * "\n"
            end
        end
        return result
    end

    function make_component(package, component)
        maker_name = Symbol(lowercase(package.name),"_", lowercase(component.name))
        props = collect(keys(component.properties))
        add_signs = ""
        if any(x->x==:children, props)
            add_signs = """\n    $(maker_name)(children::Any;kwags...)\n    $(maker_name)(children_maker::Function;kwags...)\n"""
        end

        avaible_props = Set(Symbol.(getproperty.(values(component.properties), :name)))

        docstr = """    $(maker_name)(;kwags...)
        $(add_signs)

        $(component.description)

        # Arguments
        $(join(map(values(component.properties)) do prop
            "- `$(prop.name)` - $(prop.description)"
        end,"\n"))
        """
        makers = Vector{Expr}()
        push!(makers,esc(quote
            export $(maker_name)
            function $(maker_name) end
            @doc $(docstr) $(maker_name)
            function $(maker_name)(;kwargs...)
                avaible_props = $(avaible_props)
                result = Component($(component.name), $(package.package_name), Dict{Symbol, Any}(), avaible_props)
                for (prop, value) in pairs(kwargs)
                    if !(prop in avaible_props)
                        throw(ArgumentError("Invalid property $(string(prop)) for component " * $(string(maker_name))))
                    end
                    push!(result.props, prop=>value)
                end
                return result
            end
        end))
        if any(x->x==:children, props)
            props_nochildren = filter(x->x!=:children, props)
            push!(makers,esc(quote
                function $(maker_name)(children::Any;  kwargs...)
                    avaible_props = $(avaible_props)
                    result = Component($(component.name), $(package.package_name), Dict{Symbol, Any}(), avaible_props)
                    for (prop, value) in pairs(kwargs)
                        if !(prop in avaible_props)
                            throw(ArgumentError("Invalid property $(string(prop)) for component " * $(string(maker_name))))
                        end
                        push!(result.props, prop=>value)
                    end
                    push!(result.props, :children=>children)
                    return result
                end
                function $(maker_name)(children_maker::Function; kwargs...)
                    avaible_props = $(avaible_props)
                    result = Component($(component.name), $(package.package_name), Dict{Symbol, Any}(), avaible_props)
                    for (prop, value) in pairs(kwargs)
                        if !(prop in avaible_props)
                            throw(ArgumentError("Invalid property $(string(prop)) for component " * $(string(maker_name))))
                        end
                        push!(result.props, prop=>value)
                    end
                    push!(result.props, :children=>children_maker())
                    return result
                end
        end))
        end

        return quote
            $(makers...)
        end
    end

    macro reg_components()
        components_code::Vector{Expr} = []
        for package in values(_components_packages)
            if package.has_components
                components = ComponentPackages.load_package_components(package.name)
                for component in values(components)
                    push!(components_code, make_component(package, component))
                end
            end
        end
        return quote
            $(components_code...)
        end

    end
end

max_stars_repo_name: EricForgy/Dashboards

max_issues_repo_name: EricForgy/Dashboards

max_forks_repo_name: EricForgy/Dashboards


text:
export getIndicesOfCellsArray;
function getIndicesOfCellsArray(Mesh::RegularMesh, overlap::Array{Int64,1},numDomains::Array{Int64,1},getIndicesOfCell::Function)
	# Mesh 				= param.Mesh;
	# overlap 			= param.overlap;
	# numDomains 			= param.numDomains
	# getIndicesOfCell	= param.getIndicesOfCell;
	ncells				= Mesh.n;
	## This is just a tryout:
	Idxs = getIndicesOfCell(numDomains,overlap,div.(ncells,2).+1,ncells);
	lenIdxs = length(Idxs);
	ArrIdxs = zeros(UInt32,lenIdxs,prod(numDomains));
	for ic = 1:prod(numDomains);
		icloc = cs2loc(ic,numDomains);
		IIp = getIndicesOfCell(numDomains,overlap,icloc,ncells);
		ArrIdxs[1:length(IIp),ic] = IIp;
	end
	return ArrIdxs;
end

export loc2cs3D
function loc2cs3D(loc::Array{Int64,1},n::Array{Int64,1})
@inbounds cs = loc[1] + (loc[2]-1)*n[1] + (loc[3]-1)*n[1]*n[2];
return cs;
end

export loc2cs
function loc2cs(loc::Array{Int64,1},n::Array{Int64,1})
if length(n)==2
	@inbounds cs = loc[1] + (loc[2]-1)*n[1];
else
	@inbounds cs = loc[1] + (loc[2]-1)*n[1] + (loc[3]-1)*n[1]*n[2];
end
return cs;
end

export cs2loc
function cs2loc(cs_loc::Int64,n::Array{Int64,1})
	if length(n)==3
		@inbounds loc1 = mod(cs_loc-1,n[1])+1;
		@inbounds loc2 = div(mod(cs_loc-1,n[1]*n[2]),n[1]) + 1;
		@inbounds loc3 = div(cs_loc-1,n[1]*n[2])+1;
		return [loc1;loc2;loc3];
	else
		@inbounds loc1 = mod(cs_loc-1,n[1]) + 1;
		@inbounds loc2 = div(cs_loc-1,n[1])+1;
		return [loc1;loc2];
	end
end


text:
#TODO: Make inmmutable later. Requires refactor of ThreePhase Fault callbacks
mutable struct SimulationInputs
    sys::PSY.System
    injectors_data::Vector{<:PSY.StaticInjection}
    counts::Base.ImmutableDict{Symbol, Int}
    Ybus::SparseMatrixCSC{Complex{Float64}, Int64}
    dyn_lines::Bool
    voltage_buses_ix::Vector{Int}
    current_buses_ix::Vector{Int}
    global_index::Dict{String, Dict{Symbol, Int64}}
    total_shunts::Dict{Int64, Float64}
    global_vars::Dict{Symbol, Number}
    lookup::Dict{Int, Int}
    DAE_vector::Vector{Bool}
    aux_arrays::Dict{Int, Vector}
    tspan::NTuple{2, Float64}
end

function SimulationInputs(;
    sys::PSY.System,
    dynamic_injectors = Vector{PSY.StaticInjection}(),
    counts::Base.ImmutableDict{Symbol, Int} = Base.ImmutableDict{Symbol, Int}(),
    Ybus::SparseMatrixCSC{Complex{Float64}, Int64} = SparseMatrixCSC{
        Complex{Float64},
        Int64,
    }(zeros(1, 1)),
    dyn_lines::Bool = false,
    voltage_buses_ix::Vector{Int} = Vector{Int}(),
    current_buses_ix::Vector{Int} = Vector{Int}(),
    global_index::Dict{String, Dict{Symbol, Int64}} = Dict{String, Dict{Symbol, Int64}}(),
    total_shunts::Dict{Int64, Float64} = Dict{Int64, Float64}(),
    global_vars::Dict{Symbol, Number} = Dict{Symbol, Number}(),
    lookup::Dict{Int, Int} = Dict{Int, Int}(),
    DAE_vector::Vector{Bool} = Vector{Bool}(),
    aux_arrays::Dict{Int, Vector} = Dict{Int, Vector}(),
    tspan::NTuple{2, Float64} = (0.0, 0.0),
)
    return SimulationInputs(
        sys,
        dynamic_injectors,
        counts,
        Ybus,
        dyn_lines,
        voltage_buses_ix,
        current_buses_ix,
        global_index,
        total_shunts,
        global_vars,
        lookup,
        DAE_vector,
        aux_arrays,
        tspan,
    )
end

function build!(inputs::SimulationInputs)
    sys = get_system(inputs)
    n_buses = length(PSY.get_components(PSY.Bus, sys))
    DAE_vector = collect(falses(n_buses * 2))
    global_state_index = MAPPING_DICT()
    state_space_ix = [n_buses * 2]
    current_buses_ix = collect(1:n_buses)
    static_bus_var_count = 2 * length(current_buses_ix)
    voltage_buses_ix = Vector{Int}()
    total_states = 0
    first_dyn_branch_point = -1
    branches_n_states = 0
    global_vars = Dict{Symbol, Number}(
        :ω_sys => 1.0,
        :ω_sys_index => -1, #To define 0 if infinite source, bus_number otherwise,
    )
    total_shunts = Dict{Int, Float64}()
    found_ref_bus = false
    sys_basepower = PSY.get_base_power(sys)

    inputs.Ybus, inputs.lookup = _get_Ybus(sys)
    dyn_branches = PSY.get_components(PSY.DynamicBranch, sys)

    if !(isempty(dyn_branches))
        inputs.dyn_lines = true
        first_dyn_branch_point = state_space_ix[1] + 1
        for br in dyn_branches
            arc = PSY.get_arc(br)
            n_states = PSY.get_n_states(br)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            bus_ix_from = get_lookup(inputs)[from_bus_number]
            bus_ix_to = get_lookup(inputs)[to_bus_number]
            merge!(
                +,
                total_shunts,
                Dict(
                    bus_ix_from => 1 / PSY.get_b(br).from,
                    bus_ix_to => 1 / PSY.get_b(br).to,
                ),
            )
            push!(voltage_buses_ix, bus_ix_from, bus_ix_to)
            DAE_vector[bus_ix_from] = DAE_vector[bus_ix_from + n_buses] = true
            DAE_vector[bus_ix_to] = DAE_vector[bus_ix_to + n_buses] = true
            DAE_vector = push!(DAE_vector, collect(trues(n_states))...)
            total_states += n_states
            _add_states_to_global!(global_state_index, state_space_ix, br)
        end

        for (ix, val) in enumerate(DAE_vector[1:n_buses])
            if val
                #This V_ix should be V_number.
                global_state_index["V_$(ix)"] = Dict(:R => ix, :I => ix + n_buses)
                total_states += 2
                static_bus_var_count -= 2
                push!(voltage_buses_ix, ix)
                @assert static_bus_var_count >= 0
            end
        end
        branches_n_states = state_space_ix[1] - n_buses * 2
    else
        inputs.dyn_lines = false
        @debug("System doesn't contain Dynamic Branches")
    end
    unique!(voltage_buses_ix)
    sources = PSY.get_components(PSY.Source, sys)
    for s in sources
        btype = PSY.get_bustype(PSY.get_bus(s))
        if (btype == PSY.BusTypes.REF) && found_ref_bus
            throw(IS.ConflictingInputsError("The system can't have more than one source or generator in the REF Bus"))
        end
        btype != PSY.BusTypes.REF && continue
        global_vars[:ω_sys_index] = 0 #To define 0 if infinite source, bus_number otherwise,
        found_ref_bus = true
    end

    dynamic_injection = PSY.get_components(
        PSY.StaticInjection,
        sys,
        x -> PSY.get_dynamic_injector(x) !== nothing,
    )

    if isempty(dynamic_injection)
        error("System doesn't contain any DynamicInjection devices")
    end

    for d in dynamic_injection
        @debug PSY.get_name(d)
        dynamic_device = PSY.get_dynamic_injector(d)
        if !(:states in fieldnames(typeof(dynamic_device)))
            continue
        end
        device_bus = PSY.get_bus(d)
        btype = PSY.get_bustype(device_bus)
        if (btype == PSY.BusTypes.REF) && found_ref_bus
            throw(IS.ConflictingInputsError("The system can't have more than one source or generator in the REF Bus"))
        end
        _make_device_index!(d)
        device_n_states = PSY.get_n_states(dynamic_device)
        DAE_vector = push!(DAE_vector, collect(trues(device_n_states))...)
        total_states += device_n_states
        _add_states_to_global!(global_state_index, state_space_ix, dynamic_device)
        push!(inputs.injectors_data, d)

        btype != PSY.BusTypes.REF && continue
        if typeof(dynamic_device) <: PSY.DynamicGenerator
            ω_ix = global_state_index[PSY.get_name(d)][:ω]
        elseif typeof(dynamic_device) <: PSY.DynamicInverter
            #TO DO: Make it general for cases when ω is not a state (droop)!
            ω_ix = global_state_index[PSY.get_name(d)][:ω_oc]
        end
        global_vars[:ω_sys_index] = ω_ix #To define 0 if infinite source, bus_number otherwise,
        found_ref_bus = true
    end
    injection_n_states = state_space_ix[1] - branches_n_states - n_buses * 2
    @assert total_states == state_space_ix[1] - static_bus_var_count
    @debug total_states
    setdiff!(current_buses_ix, voltage_buses_ix)

    inputs.global_vars = global_vars
    inputs.DAE_vector = DAE_vector
    inputs.global_index = global_state_index
    inputs.voltage_buses_ix = voltage_buses_ix
    inputs.current_buses_ix = current_buses_ix
    inputs.total_shunts = total_shunts
    inputs.counts = Base.ImmutableDict(
        :total_states => total_states,
        :injection_n_states => injection_n_states,
        :branches_n_states => branches_n_states,
        :first_dyn_injection_pointer => 2 * n_buses + branches_n_states + 1,
        :first_dyn_branch_point => first_dyn_branch_point,
        :total_variables => total_states + static_bus_var_count,
        :bus_count => n_buses,
    )

    @assert get_ω_sys(inputs) != -1
    return
end

get_system(inputs::SimulationInputs) = inputs.sys
get_injectors_data(inputs::SimulationInputs) = inputs.injectors_data
get_counts(inputs::SimulationInputs) = inputs.counts
get_voltage_buses_ix(inputs::SimulationInputs) = inputs.voltage_buses_ix
get_current_buses_ix(inputs::SimulationInputs) = inputs.current_buses_ix
get_global_index(inputs::SimulationInputs) = inputs.global_index
get_Ybus(inputs::SimulationInputs) = inputs.Ybus
get_total_shunts(inputs::SimulationInputs) = inputs.total_shunts
get_global_vars(inputs::SimulationInputs) = inputs.global_vars
get_dyn_lines(inputs::SimulationInputs) = inputs.dyn_lines
get_lookup(inputs::SimulationInputs) = inputs.lookup
get_DAE_vector(inputs::SimulationInputs) = inputs.DAE_vector
get_aux_arrays(inputs::SimulationInputs) = inputs.aux_arrays
get_tspan(inputs::SimulationInputs) = inputs.tspan

get_injection_pointer(inputs::SimulationInputs) =
    get_counts(inputs)[:first_dyn_injection_pointer]
get_branches_pointer(inputs::SimulationInputs) = get_counts(inputs)[:first_dyn_branch_point]
get_n_injection_states(inputs::SimulationInputs) = get_counts(inputs)[:injection_n_states]
get_n_branches_states(inputs::SimulationInputs) = get_counts(inputs)[:branches_n_states]
get_system_state_count(inputs::SimulationInputs) = get_counts(inputs)[:total_states]
get_variable_count(inputs::SimulationInputs) = get_counts(inputs)[:total_variables]
get_device_index(inputs::SimulationInputs, device::D) where {D <: PSY.DynamicInjection} =
    get_global_index(inputs)[device.name]
get_bus_count(inputs::SimulationInputs) = get_counts(inputs)[:bus_count]
get_ω_sys(inputs::SimulationInputs) = get_global_vars(inputs)[:ω_sys]


text:
#=
# MNIST

Trained without centering.
=#

import Makie
import CairoMakie
import MLDatasets
import Flux
import RestrictedBoltzmannMachines as RBMs
using Statistics: mean, std, var
using Random: bitrand
using ValueHistories: MVHistory
using RestrictedBoltzmannMachines: visible, BinaryRBM, sample_from_inputs
using RestrictedBoltzmannMachines: initialize!, log_pseudolikelihood, pcd!, minibatch_count
nothing #hide

# Useful function to plot grids of MNIST digits.

"""
    imggrid(A)

Given a four dimensional tensor `A` of size `(width, height, ncols, nrows)`
containing `width x height` images in a grid of `nrows x ncols`, this returns
a matrix of size `(width * ncols, height * nrows)`, that can be plotted in a heatmap
to display all images.
"""
imggrid(A::AbstractArray{<:Any,4}) =
    reshape(permutedims(A, (1,3,2,4)), size(A,1)*size(A,3), size(A,2)*size(A,4))

#=
Load the MNIST dataset.
=#

Float = Float32
train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float}(train_x[:, :, train_y .== 0] .≥ 0.5)
nothing #hide

# Let's visualize some random digits.

nrows, ncols = 10, 15
fig = Makie.Figure(resolution=(40ncols, 40nrows))
ax = Makie.Axis(fig[1,1], yreversed=true)
idx = rand(1:size(train_x,3), nrows * ncols) # random indices of digits
digits = reshape(train_x[:,:,idx], 28, 28, ncols, nrows)
Makie.image!(ax, imggrid(digits), colorrange=(1,0))
Makie.hidedecorations!(ax)
Makie.hidespines!(ax)
fig

# Initialize an RBM with 400 hidden units.

rbm = initialize!(BinaryRBM(Float, (28,28), 400), train_x)
batchsize = 256
batchcount = minibatch_count(train_x; batchsize)
epochs = 500
history = MVHistory()
time_0 = time()
@time pcd!(
    rbm, train_x; epochs, batchsize, center=false,
    callback = function(; epoch, batch_idx, _...)
        push!(history, :t, time() - time_0)
        if batch_idx == batchcount && epoch % 5 == 0
            lpl = log_pseudolikelihood(rbm, train_x)
            push!(history, :lpl, mean(lpl))
        end
    end
)
nothing #hide

# Plot of log-pseudolikelihood of trian data during learning.

fig = Makie.Figure(resolution=(500,300))
ax = Makie.Axis(fig[1,1], xlabel = "train time", ylabel="pseudolikelihood")
Makie.lines!(ax, get(history, :lpl)...)
fig

# Sample digits from the RBM starting from a random condition.

nsteps = 3000
fantasy_F = zeros(nrows*ncols, nsteps)
fantasy_x = bitrand(28,28,nrows*ncols)
fantasy_F[:,1] .= RBMs.free_energy(rbm, fantasy_x)
@time for t in 2:nsteps
    fantasy_x .= RBMs.sample_v_from_v(rbm, fantasy_x)
    fantasy_F[:,t] .= RBMs.free_energy(rbm, fantasy_x)
end
nothing #hide

# Check equilibration of sampling

fig = Makie.Figure(resolution=(400,300))
ax = Makie.Axis(fig[1,1], xlabel="sampling time", ylabel="free energy")
fantasy_F_μ = vec(mean(fantasy_F; dims=1))
fantasy_F_σ = vec(std(fantasy_F; dims=1))
Makie.band!(ax, 1:nsteps, fantasy_F_μ - fantasy_F_σ/2, fantasy_F_μ + fantasy_F_σ/2)
Makie.lines!(ax, 1:nsteps, fantasy_F_μ)
fig

# Plot the sampled digits.

fig = Makie.Figure(resolution=(40ncols, 40nrows))
ax = Makie.Axis(fig[1,1], yreversed=true)
Makie.image!(ax, imggrid(reshape(fantasy_x, 28, 28, ncols, nrows)), colorrange=(1,0))
Makie.hidedecorations!(ax)
Makie.hidespines!(ax)
fig


text:
import Optim, LineSearches

struct GradientAscent <: Algorithm end

function train(learner::Learner{GradientAscent})
    num_samples = sample_count(learner.train_data)
    β_init = _unif_random(learner.param_lb, learner.param_ub)
    offset_init = rand(Uniform(learner.offset_lb, learner.offset_ub))

    l1_coeff = 0.0
    l2_coeff = 0.0
    for penalty in learner.objective_penalties
        if penalty isa L1Regularizer
            @assert l1_coeff == 0.0
            l1_coeff = penalty.coefficient
        elseif penalty isa L2Regularizer
            @assert l2_coeff = 0.0
            l2_coeff = penalty.coefficient
        else
            error("Unrecognized type of regularizer")
        end
    end
    function l(x)
        β = x[1:(end-1)]
        offset = x[end]
        model = LinearModel(β, offset)
        return -1 * (1 / num_samples * sum(_reward_function(data_point, model) for data_point in learner.train_data.samples) +
            l1_coeff * sum(abs.(β)) +
            l2_coeff * sum(β.^2))
    end
    function ∇l(x)
        β = x[1:(end-1)]
        offset = x[end]
        model = LinearModel(β, offset)
        return -1 * (1 / num_samples * sum(vcat(_grad_reward_function(data_point, model)...) for data_point in learner.train_data.samples) +
            l1_coeff * vcat(sign.(β), 0.0) +
            2l2_coeff * vcat(β, 0.0))
    end

    results = Optim.optimize(l, ∇l, vcat(β_init, offset_init), Optim.GradientDescent(linesearch=LineSearches.StrongWolfe()), Optim.Options(show_trace=false, allow_f_increases=true, extended_trace=false), inplace=false)

    β_values = results.minimizer[1:(end-1)]
    offset_value = results.minimizer[end]

    return TrainingResults(learner, LinearModel(β_values, offset_value))
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
"""
    @option [alias::String] <struct def>

Define an option struct type. This will auto-generate methods that parse a given `Dict{String}`
object (the keys must be of type `String`) into an instance of the struct type you defined. One
can use `alias` string to distinguish multiple possible option type for the same field.

# Example

One can define option type via `@option` macro with or without an alias.

```julia-repl
julia> "Option A"
       @option "option_a" struct OptionA
           name::String
           int::Int = 1
       end

julia> "Option B"
       @option "option_b" struct OptionB
           opt::OptionA = OptionA(;name = "Sam")
           float::Float64 = 0.3
       end
```

and convert a dict to an option type via [`from_dict`](@ref).

```julia-repl
julia> d = Dict{String, Any}(
           "opt" => Dict{String, Any}(
               "name" => "Roger",
               "int" => 2,
           ),
           "float" => 0.33
       );

julia> option = from_dict(OptionB, d)
OptionB(;
    opt = OptionA(;
        name = "Roger",
        int = 2,
    ),
    float = 0.33,
)
```

when there are multiple possible option type for one field,
one can use the alias to distinguish them

```julia-repl
julia> @option struct OptionD
           opt::Union{OptionA, OptionB}
       end

julia> d1 = Dict{String, Any}(
               "opt" => Dict{String, Any}(
                   "option_b" => d
               )
           );

julia> from_dict(OptionD, d1)
OptionD(;
    opt = OptionB(;
        opt = OptionA(;
            name = "Roger",
            int = 2,
        ),
        float = 0.33,
    ),
)
```
"""
macro option(ex)
    esc(option_m(ex))
end

macro option(alias::String, ex)
    esc(option_m(ex, alias))
end

function option_m(ex, type_alias=nothing)
    def = JLKwStruct(ex, type_alias)
    return codegen_option_type(def)
end

"""
    codegen_option_type(def::JLKwStruct)

Generate the `Configurations` option type definition from
a given `JLKwStruct` created by [`Expronicon`](https://github.com/Roger-luo/Expronicon.jl).
"""
function codegen_option_type(def::JLKwStruct)
    quote
        $(codegen_ast(def))
        Core.@__doc__ $(def.name)
        $(codegen_create(def))
        $(codegen_is_option(def))
        $(codegen_convert(def))
        $(codegen_field_default(def))
        $(codegen_type_alias(def))
        $(codegen_isequal(def))
        # pretty printings
        $(codegen_show(def))
        nothing
    end
end

"""
    codegen_is_option(x::JLKwStruct)

Generate the [`is_option`](@ref) method.
"""
function codegen_is_option(def::JLKwStruct)
    quote
        $Configurations.is_option(::$(def.name)) = true
        $Configurations.is_option(::Type{<:$(def.name)}) = true
    end
end

"""
    codegen_convert(x::JLKwStruct)

Generate `Base.convert` from `AbstractDict{String}` to the given option type.
"""
function codegen_convert(def::JLKwStruct)
    quote
        function $Base.convert(::Type{<:$(def.name)}, d::AbstractDict{String})
            $Configurations.from_dict($(def.name), d)
        end
    end
end

function codegen_field_default(def::JLKwStruct)
    obj = gensym(:x)
    msg = Expr(:string, "type $(def.name) does not have field ", obj)
    err = :(error($msg))
    ret = JLIfElse()
    ret.otherwise = err

    isconst = Dict{Symbol, Bool}()
    default = Dict{Symbol, Any}()
    prev_field_names = Symbol[]

    for (k, field) in enumerate(def.fields)
        vars = Symbol[]

        for name in prev_field_names
            if has_symbol(field.default, name)
                push!(vars, name)
            end
        end

        cond = :($obj == $(QuoteNode(field.name)))

        if isempty(vars) # const default
            ret[cond] = field.default
        else
            jlfn = JLFunction(;
                args=vars,
                body=field.default,
            )
            fn = gensym(:fn)
            ret[cond] = quote
                $fn = $(codegen_ast(jlfn))
                $PartialDefault($fn, $vars, $(Expr(:quote, field.default)))
            end
        end

        push!(prev_field_names, field.name)
    end

    type = gensym(:type)
    typevars = name_only.(def.typevars)
    ub = isempty(def.typevars) ? def.name : Expr(:curly, def.name, typevars...)

    return codegen_ast(
        JLFunction(;
            name=:($Configurations.field_default),
            args=[:(::Type{$type}), :($obj::Symbol)],
            body=codegen_ast(ret),
            whereparams=[typevars..., :($type <: $ub)],
        )
    )
end

# function codegen_field_alias(def::JLKwStruct)
# end

"""
    codegen_isequal(x::JLKwStruct)

Generate `Base.:(==)` to overload comparison operator to [`compare_options`](@ref)
for given option type.
"""
function codegen_isequal(x::JLKwStruct)
    return :(Base.:(==)(a::$(x.name), b::$(x.name)) = $compare_options(a, b))
end

"""
    codegen_type_alias(def::JLKwStruct)

Generate type alias method [`type_alias`](@ref).
"""
function codegen_type_alias(def::JLKwStruct)
    quote
        $Configurations.type_alias(::Type{<:$(def.name)}) = $(def.typealias)
    end
end

function codegen_create(def::JLKwStruct)
    codegen_ast_kwfn(def, :($Configurations.create))
end

function codegen_show(def::JLKwStruct)
    quote
        function $Base.show(io::IO, m::MIME"text/plain", x::$(def.name))
            $show_option(io, m, x)
        end

        function $Base.show(io::IO, m::MIME"text/html", x::$(def.name))
            $show_option(io, m, x)
        end

        function $Base.show(io::IO, ::MIME"application/toml", x::$(def.name))
            return print(io, to_toml(x))
        end
    end
end


text:
using Documenter, SOFA

# Generate documents
makedocs(
    modules   = [SOFA],
    doctest   = false,
    clean     = true,
    linkcheck = true,
    format    = Documenter.HTML(),
    sitename  = "SOFA.jl",
    authors   = "Duncan Eddy",
    pages     = Any[
        "Home" => "index.md",
        "SOFA Library" => "sofa.md",
    ]
)

deploydocs(
    repo = "github.com/sisl/SOFA.jl",
    devbranch = "master",
    devurl = "latest",
)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Distributions
using Plots

Distributions.Normal{T}(; μ, σ) where T = Normal{T}(μ, σ)
normals = [
    Normal{Float64}(μ=-0.42117172154916505, σ=1.2253594150160463)
    Normal{Float64}(μ=1.4244177594878171, σ=1.2253594150160463)
    Normal{Float64}(μ=-0.7646882557011875, σ=1.2253594150160463)
    Normal{Float64}(μ=-1.4994608515495296, σ=1.2253594150160463)
    Normal{Float64}(μ=0.1588834413196244, σ=1.2253594150160463)
    Normal{Float64}(μ=0.6906462431515679, σ=1.2253594150160463)
    Normal{Float64}(μ=-0.8610106914769075, σ=1.2253594150160463)
    Normal{Float64}(μ=0.39571181451559956, σ=1.2253594150160463)
    Normal{Float64}(μ=-0.7506648963208741, σ=1.2253594150160463)
    Normal{Float64}(μ=1.0347673088571538, σ=1.2253594150160463)
] |> x -> sort(x; by=mean)

ys = range(-10, 10; length=300)
m, n = length(normals), length(ys)
plot(legend=false)
for (x, normal) in zip(range(-1.5, 1.5; length=m), normals)
    plot3d!(fill(x, n), ys, pdf.(normal, ys))
end
plot!()

# %%
using LinearAlgebra
using Distributions
using Plots

dist_true = MvNormal([50, 50], 25^2*[1 1; 1 1.3])
n = 30
sample = rand(dist_true, n)
x, y = sample[1,:], sample[2,:]
X = x .^ (0:1)'
b̂ = X\y # b\hat TAB -> b̂
ŷ = X*b̂
u = norm(y - ŷ)/√(n - 2)
f(x) = evalpoly(x, b̂)
s(x) = u*√(1 + dot([1, f(x)], X'X\[1, f(x)]))
#disty(x) = Normal(f(x), s(x))
disty(x) = LocationScale(f(x), s(x), TDist(n-2))

plot(legend=false, xtick=-100:50:200, ytick=-100:50:200)
for xs in 0:10:100
    ys = range(-50, 150; length=200)
    plot3d!(fill(xs, length(ys)), ys, pdf.(disty(xs), ys))
end
scatter3d!(x, y, zeros(n); ms=2, msc=:auto, color=:red)
xs = range(-30, 130; length=100)
plot3d!(xs, f.(xs), zero(xs); color=:blue)

# %%

max_stars_repo_path: 0015/plot3d example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/plot3d example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/plot3d example.jl

max_forks_repo_name: genkuroki/public


text:
"""
    enforce_no_flow_through!(b::ConformalBody, motion, elements, t)

Update the body, `b`, to enforce the no-flow-through condition given ambient vortex elements, `elements`, and while moving with kinematics specified by `motion`.

# Example

```jldoctest
julia> p = Bodies.Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0])
Polygon with 4 vertices at
             (-1.0,-1.0) (0.2,-1.0) (1.0,0.5) (-1.0,1.0)
             interior angles/π = [0.5, 0.656, 0.422, 0.422]

julia> b = Bodies.ConformalBody(p)
Body generated by: Schwarz-Christoffel map of unit circle to exterior of polygon with 4 vertices

  centroid at 0.0 + 0.0im
  angle 0.0

julia> motion = RigidBodyMotion(1.0, 0.0);

julia> point = Vortex.Point(0.0 + 2im, 1.0);

julia> Bodies.enforce_no_flow_through!(b, motion, point, 0.0)

julia> b.img
1-element Array{Element,1}:
 Vortex.Point(0.0 + 0.5im, -1.0)
```
"""
function enforce_no_flow_through!(b::ConformalBody, ṗ, elements, t)
    @get ṗ (ċ, α̇)

    # should set up images here

    b.ċ = ċ
    b.α̇ = α̇

    get_image!(b,elements)

    nothing
end

Elements.image(z::ComplexF64,b::ConformalBody) = 1.0/conj(z)

Elements.image(s::T,b::ConformalBody) where T <: Union{Blob,Point} = Elements.image(s.z,b)

# TODO: The `get_image` functions are causing some type instability
# might have to fix it later
function get_image!(b::ConformalBody, sources::T) where T <: Union{Tuple, AbstractArray}
    b.img = Elements.Element[]  #Points.Point[]
    for source in sources
        get_image!(b.img,source,b)
    end
    nothing
end

function get_image!(b::ConformalBody, source::T) where T <: Union{Blob,Point}
  get_image!(b.img,source,b)
end

get_image!(tmp,src,b::ConformalBody) = get_image!(tmp,Elements.unwrap_src(src), b, kind(Elements.unwrap_src(src)))

function get_image!(tmp,src,b::ConformalBody,::Type{Singleton})
  push!(tmp,get_image(src,b))
end

function get_image!(tmp,src,b::ConformalBody,::Type{Group})
  for i in eachindex(src)
      push!(tmp,get_image(src[i], b))
  end
end

function get_image(src::Union{Blob{T},Point{T}}, b::ConformalBody) where T <: Complex
    Point{T}(Elements.image(src.z,b),src.S)
end

function get_image(src::Union{Blob{T},Point{T}}, b::ConformalBody) where T <: Real
    Point{T}(Elements.image(src.z,b),-src.S)
end

function get_image(src::Freestream, b::ConformalBody)
    # it is assumed that freestream is given in the circle plane form. Note that
    # src.U is already the conjugate of the free stream (due to the convention
    # in this code), so there is no need to take the conjugate of it here
    Doublet{ComplexF64}(0.0,src.U*π)
end

"""
    vorticity_flux(b::ConformalBody, edge, sys, v, t,
                   tesp = 0.0)

Return strength of a new vortex element that satisfies edge suction parameter on an edge of a conformally-mapped body.
For a given edge, if the current suction parameter is less than the criticial suction parameter, then no vorticity is released.  If it is higher, however, vorticity will be released so that the suction parameter equals the critical value.

# Arguments

- `b`: the body
- `edge`: index of the vertex on the body corresponding to designated edge
- `sys`: current system of body and fluid vorticity
- `v`: the vortex element (with unit circulation) that the vorticity flux is going into
- `t`: the current time (for evaluating body motion)
- `tesp`: the critical trailing edge suction parameter we want to enforce.  By default, the parameters is set to 0.0 to enforce the Kutta condition on the edge.  We can disable vortex shedding from an edge by setting the its critical suction parameter to `Inf`

# Returns

- `Γ`: the strength that the vortex element should have in order to satisfy the edge suction parameters

"""
function vorticity_flux(b::Bodies.ConformalBody, edge::Integer, sys, v, t, tesp = 0.0)
    σ̃ = suction_parameter(edge,b,sys,t)

    # enforce boundary conditions for the elements v on a stationary airfoil
    db = deepcopy(b)
    motion = RigidBodyMotion(0.0, 0.0)
    Bodies.enforce_no_flow_through!(db, motion, v, 0)

    dσ = suction_parameter(edge,b,(db,v),t)

    Γ = circulation(v)

    if (abs2(tesp) > abs2(σ̃))
        K = 0.0
    else
        K = (sign(σ̃)*tesp - σ̃)/dσ
    end
    return K*Γ
end

"""
    vorticity_flux(b::ConformalBody, edge₁, edge₂, sys, v₁, v₂, t,
                   sp₁ = 0.0, sp₂ = 0.0)

Return strength of two new vortex elements that satisfies edge suction parameter on the specified edges of a conformally-mapped body.
For a given edge, if the current suction parameter is less than the criticial suction parameter, then no vorticity is released.
If it is higher, however, vorticity will be released so that the suction parameter equals the critical value.

# Arguments

- `b`: the body
- `edge₁`: index of the vertex on the body corresponding to first designated edge
- `edge₂`: index of the vertex on the body corresponding to second designated edge
- `sys`: current system of body and fluid vorticity
- `v₁`: the first vortex element (with unit circulation) that the vorticity flux is going into
- `v₂`: the second vortex element (with unit circulation) that the vorticity flux is going into
- `t`: the current time (for evaluating body motion)
- `sp₁`: the critical edge suction parameter at edge₁ we want to enforce. By default, the parameters is set to 0.0 to enforce the Kutta condition on the edge.  We can disable vortex shedding from an edge by setting the its critical suction parameter to `Inf`
- `sp₂`: the critical edge suction parameter at edge₂ we want to enforce. By default, the parameters is set to 0.0 to enforce the Kutta condition on the edge.  We can disable vortex shedding from an edge by setting the its critical suction parameter to `Inf`

# Returns

- `Γ₁`: the strength that the vortex element `v₁` should have in order to satisfy the edge suction parameters
- `Γ₂`: the strength that the vortex element `v₂` should have in order to satisfy the edge suction parameters

"""
function vorticity_flux(b::Bodies.ConformalBody, edge₁::Integer, edge₂::Integer, sys, v₁, v₂, t, sp₁ = 0.0, sp₂ = 0.0)

    # existing suction parameters
    σ̃₁ = Bodies.suction_parameter(edge₁,b,sys,t)
    σ̃₂ = Bodies.suction_parameter(edge₂,b,sys,t)

    # enforce boundary conditions for the elements v on a stationary airfoil
    db₁ = deepcopy(b)
    db₂ = deepcopy(b)
    motion = RigidBodyMotion(0.0, 0.0)
    Bodies.enforce_no_flow_through!(db₁, motion, v₁, 0)
    Bodies.enforce_no_flow_through!(db₂, motion, v₂, 0)

    # Unit suction parameters
    dσ₁₁ = suction_parameter(edge₁,b,(db₁,v₁),t)
    dσ₁₂ = suction_parameter(edge₂,b,(db₁,v₁),t)

    dσ₂₁ = suction_parameter(edge₁,b,(db₂,v₂),t)
    dσ₂₂ = suction_parameter(edge₂,b,(db₂,v₂),t)

    Γ₁ = circulation(v₁)
    Γ₂ = circulation(v₂)

    if (abs2(sp₁) > abs2(σ̃₁)) && (abs2(sp₂) ≤ abs2(σ̃₂))
        K₁, K₂ = 0.0, (sign(σ̃₂)*sp₂ - σ̃₂)/dσ₂₂
    elseif (abs2(sp₁) ≤ abs2(σ̃₁)) && (abs2(sp₂) > abs2(σ̃₂))
        K₁, K₂ = (sign(σ̃₁)*sp₁ - σ̃₁)/dσ₁₁, 0.0
    elseif (abs2(sp₁) > abs2(σ̃₁)) && (abs2(sp₂) > abs2(σ̃₂))
        K₁ = K₂ = 0.0
    else
        σ̃₁ = sign(σ̃₁)*sp₁ - σ̃₁
        σ̃₂ = sign(σ̃₂)*sp₂ - σ̃₂
        detA = dσ₁₁*dσ₂₂ - dσ₁₂*dσ₂₁

        @assert (detA != 0) "Cannot enforce suction parameters"

        K₁ = (dσ₂₂*σ̃₁ - dσ₂₁*σ̃₂)/detA
        K₂ = (dσ₁₁*σ̃₂ - dσ₁₂*σ̃₁)/detA

    end
    return K₁*Γ₁, K₂*Γ₂
end

function suction_parameter_factor(k::Integer,m::ExteriorMap)
    # Return the factor in front of the velocity tangent to the circle for vertex k
    beta = 1 .- m.angle
    zeta = m.preprev
    fact = (1+beta[k])^beta[k]*abs(m.constant)
    for j = 1:m.N
        if j == k
            continue
        end
        fact *= abs(zeta[k]-zeta[j])^beta[j]
    end
    return -fact^(-1/(1+beta[k]))
end
function suction_parameter(edge::Integer,b::Bodies.ConformalBody,sys,t)
    w̃ = induce_velocity(b,sys,t)
    σ = real(-im*conj(b.zetas[edge])*w̃[edge])*suction_parameter_factor(edge,b.m)
    return σ
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie, FileIO
GLMakie.activate!() # HIDE
let
    moon = load(download("https://svs.gsfc.nasa.gov/vis/a000000/a004600/a004675/phases.0001_print.jpg"))
    n = 1024 # 2048
    θ = LinRange(0,pi,n)
    φ = LinRange(0,2*pi,2*n)
    x = [cos(φ)*sin(θ) for θ in θ, φ in φ]
    y = [sin(φ)*sin(θ) for θ in θ, φ in φ]
    z = [cos(θ) for θ in θ, φ in φ]
    fig = Figure(resolution = (900,900), backgroundcolor = "#748AA6")
    ax = Axis3(fig, aspect = :data, viewmode = :fitzoom, #perspectiveness = 0.5,
        azimuth = 0.01π, elevation = 0.85π,)
    surface!(ax, x, y, z, color = moon, shading = true,
        lightposition = Vec3f0(-2, -3, -3), ambient = Vec3f0(0.8, 0.8, 0.8),
        backlight = 1.5f0)
    hidedecorations!(ax)
    hidespines!(ax)
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "moon.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["GLMakie","FileIO", "ImageMagick", "QuartzImageIO"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module EntropyStableEuler

using StaticArrays
using UnPack

export Euler
Base.@kwdef struct Euler{d}
    γ::Float64 = 1.4 # default value of γ
end

export logmean
include("logmean.jl")

export betafun, pfun, prim_to_cons, cons_to_prim_beta
export Sfun, u_vfun, v_ufun
include("variables.jl")

export fS, fS_prim, fS_prim_log
include("fluxes.jl")

end


text:
using GeneralizedNullspaces
using GenericLinearAlgebra
using Test
using Random
using LinearAlgebra

Random.seed!(1234)

include("testutils.jl")

quiet = true

# example from appendix B of the Guglielmi et al. paper
let tol=1e-10
    a = 233/896
    b = 248/896
    c = 171/896
    d = 29/896
    e = 15/896
    z = 69/448
    s = 2101/9632
    t = 295/19264
    u = 1403/28896

    A = [z    s      s      s      t     u     t     u     t     u;
         a    b      c      c      e     d     0     0     0     d;
         a    c      b      c      0     d     e     d     0     0;
         a    c      c      b      0     0     0     d     e     d;
         3/32 7/16   3/32   3/32   3/32  3/32  0     0     0     3/32;
         9/64 39/128 39/128 3/64   3/128 9/64  3/128 1/128 0     1/128;
         3/32 3/32   7/16   3/32   0     3/32  3/32  3/32  0     0;
         9/64 3/64   39/128 39/128 0     1/128 3/128 9/64  3/128 1/128;
         3/32 3/32   3/32   7/16   0     0     0     3/32  3/32  3/32;
         9/64 39/128 3/64   39/128 3/128 1/128 0     1/128 3/128 9/64]

    B, V, ν, μ = gnsd(A,tol)

    # check destructuring
    F = gnsd(A,tol)
    @test F isa GNSD{Float64}
    @test F.B == B
    @test F.V == V
    @test F.ν == ν
    @test F.μ == μ

    @test ν == 2
    @test μ == [3,1]
    @test gnsdcheck(A,B,V,ν,μ)
end

let tol=1e-9
    ν=1
    μx = [3,2]
    AA = makegnb(1,μx)
    n = size(AA,1)
    Q,R = qr(randn(n,n))
    A = Q' * AA * Q
    B,V,ν1,μ = gnsd(A,tol)
    @test ν == 1
    @test μ == μx[1:ν]
    @test gnsdcheck(A,B,V,ν,μ)
end

let tol=1e-9
    ν = 2
    μx = [3,2,1]
    AA = makegnb(ν,μx)
    n = size(AA,1)
    Q,R = qr(randn(n,n))
    A = Q' * AA * Q
    B,V,ν1,μ = gnsd(A,tol)
    quiet || println("input ν=$ν μ'=$μx")
    quiet || println("ν=$ν1 μ=$μ")
    @test ν == ν1
    @test gnsdcheck(A,B,V,ν1,μ; quiet=quiet)
end

for TT in (Float64,Complex{Float64})
    let tol=1e-9
        ν=2
        μx = [3,2,1]
        AA = makegnb(ν,μx)
        n = size(AA,1)
        Q,R = qr(randn(n,n))
        A = Matrix{TT}(Q' * AA * Q)
        B,V,ν1,μ = gnsd(A,tol)
        quiet || println("input ν=$ν μ'=$μx")
        quiet || println("ν=$ν1 μ=$μ")
        @test ν == ν1
        @test gnsdcheck(A,B,V,ν1,μ; quiet=quiet)
    end
end

let tol=1e-9
    ν=2
    μx=[3,2,1]
    AA = makegnb(ν,μx;cplx=true)
    n = size(AA,1)
    Q,R = qr(randn(n,n))
    A = Q' * AA * Q
    B,V,ν1,μ = gnsd(A,tol)
    quiet || println("input ν=$ν μ'=$μx")
    quiet || println("ν=$ν1 μ=$μ")
    @test ν == ν1
    @test gnsdcheck(A,B,V,ν1,μ; quiet=quiet)
end

let tol=1e-29
    ν=1
    μx = [3,2]
    AA = BigFloat.(makegnb(1,μx))
    n = size(AA,1)
    Q,R = qr(BigFloat.(randn(n,n)))
    Ab = Q' * AA * Q
    B,V,ν1,μ = gnsd(Ab,tol)
    @test ν == 1
    @test μ == μx[1:ν]
    @test gnsdcheck(Ab,B,V,ν,μ)
end


text:
module OptimBase

using Reexport
using Printf
@reexport using NLSolversBase

import Base.summary

export OptimizationOptions,
       OptimizationState,
       OptimizationTrace,
       Optimizer,
       UnivariateOptimizationResults,
       MultivariateOptimizationResults,
       # API
       minimizer,
       converged

include("types.jl")
include("api.jl")

end # module


text:
using ADCME
using ADSeismic
using PyPlot
using DelimitedFiles
using JLD2

scale = 201

param = AcousticPropagatorParams(NX=scale, NY=scale,
    NSTEP=1000, DELTAT=1e-4,  DELTAX=1.0, DELTAY=1.0,
    PropagatorKernel = 0, Rcoef = 1e-8)

rc = Ricker(param, 30.0, 200.0, 1e6)
srci = [div(param.NX,2)]
srcj = [div(param.NY,5)]
srcv = reshape(rc, :, 1)
src = AcousticSource(srci, srcj, srcv)

C = Variable(3300*ones(param.NX+2, param.NY+2))
model = AcousticPropagatorSolver(param, src, C)
@load "data.jld2" u

U = model.u
loss = sum((u[:, :, 40] - U[:,:,40])^2)
g = gradients(loss, C)

sess = Session(); init(sess)
G = run(sess, g)

close("all")
G[:,1:49] .= NaN
pcolormesh(G')
gca().invert_yaxis()
savefig("gradient.png")


text:
export Map_Affine
export eval_map

# ****************************************************************

@doc raw"""
Defines an affine transformation:

```math
X(\hat{X}) = L_A(\hat{X}) X_A θ_A +  L_B(\hat{X}) X_B θ_B
```
where ``L_A`` and ``L_B`` are the Lagrange basis.

# Constructors

## Scaled interpolation

```julia
Map_Affine(X_hat_A => X_A, X_hat_B => X_B)
```

for ``(θ_A,θ_B) = (1,1)`` we have an affine map such that
``X(\hat{X}_A) = X_A`` and ``X(\hat{X}_B) = X_B``.

## Usual interpolation

The second constructor
```julia
Map_Affine(X_hat_A, X_hat_B)
```

is a shortcut for

```julia
Map_Affine(X_hat_A => 1, X_hat_B => 1)
```
that defines an **usual** interpolation (without scaling):
``(θ_A,θ_B)`` are simply the values at ``(X_A,X_B)``
"""
struct Map_Affine{T} <: Abstract_Map
    _hat_to_A::Pair{T,T}
    _hat_to_B::Pair{T,T}
end

function Map_Affine(x_hat_A::T,x_hat_B::T) where {T}
    Map_Affine(x_hat_A=>one(T),x_hat_B=>one(T))
end

parameter_size(m::Map_Affine) = 2

function eval_map(m::Map_Affine{T},
                X_hat::AbstractVector{T_HAT},
                θ::AbstractVector{T_θ}) where {T_HAT,T,T_θ}

    @assert length(θ) == parameter_size(m)

    X_T = promote_type(T_HAT,T,T_θ)

    θ_A = θ[1]
    θ_B = θ[2]

    X_hat_A = first(m._hat_to_A)
    X_A = last(m._hat_to_A)
    X_hat_B = first(m._hat_to_B)
    X_B = last(m._hat_to_B)

    Δ = 1/(X_hat_B-X_hat_A)

    basis_A(x) = (X_hat_B-x)*Δ
    basis_B(x) = (x-X_hat_A)*Δ
    transform(x) = basis_A(x)*X_A*θ_A + basis_B(x)*X_B*θ_B

    X=similar(X_hat,X_T)
    @. X = transform(X_hat)

    X
end

# note: maybe better to use:
# @article{article,
# author = {Murray, Kevin and Müller, Samuel and Turlach, Berwin},
# year = {2016},
# month = {02},
# pages = {1-21},
# title = {Fast and flexible methods for monotone polynomial fitting},
# volume = {86},
# journal = {Journal of Statistical Computation and Simulation},
# doi = {10.1080/00949655.2016.1139582}
# }


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, Debugger
path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/ANNBN.jl"))

# This is a simple example in 1 dimension. You may set the unknown function below,
# as well as the level of noise, and number of observations
ff(x)=0.3sin(exp(3x))+0.5; noise_multiplier=0.0; i_train=100

i_test=i_train-1
vars=1
xx_train=range(0.05,0.95,length=i_train)
yy_train_init=ff.(xx_train)
yy_train=yy_train_init .+ (rand(i_train).-1/2)*noise_multiplier
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=0.8;yy_train.+=0.1
scatter(xx_train,yy_train,label="input yy_train",color=:black)
xx_test=(xx_train[2:end].+xx_train[1:end-1])./2
yy_test=ff.(xx_test)
yy_test.-=mi;yy_test./=ma;yy_test.*=0.8;yy_test.+=0.1
scatter!(xx_test,yy_test,label="initial yy_test, without noise",color=:black,markershape=:diamond,legend=:bottomleft)
xx_train=collect(xx_train)[:,:];xx_test=collect(xx_test)[:,:]

##### Clustering
neurons=Int64(floor(i_train/(vars+1))) # this automatically generates number of neurons which equals to the number of clusters (§3.1)
inds_all,n_per_part=ANNBN.___clustering(neurons,xx_train,200)

##### Train Sigmoid
a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
predl1,layer1_train=ANNBN.predict_new(a_all,a_layer1,xx_train,i_train,neurons)
predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
scatter!(xx_test,predl1_test,label=@sprintf("ANNBN-%1d-neurons",neurons),color=:black,markershape=:rect)
# savefig("1D-noise.pdf")

##### Train RBF
# Select 1 Gaussian from calc_phi.jl
cc1=0.01
a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1)
predl1,layer1_train=ANNBN.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))


text:
# Translating
# https://proplot.readthedocs.io/en/latest/subplots.html#Subplot-numbers-and-a-b-c-labels
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

f, axs = proplot.subplots(nrows=8, ncols=8, axwidth=0.7, space=0)
[ax.format(
    abc=true, abcloc="ur", xlabel="x axis", ylabel="y axis",
    xticks=[], yticks=[], suptitle="Subplot labels demo"
) for ax in axs]
f.savefig("figs/Subplot-numbers-and-a-b-c-labels.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
function ComplexODEState(arg0::jdouble, arg1::Vector{JComplex})
    return ComplexODEState((jdouble, Vector{JComplex}), arg0, arg1)
end

function ComplexODEState(arg0::jdouble, arg1::Vector{JComplex}, arg2::Vector{Vector{JComplex}})
    return ComplexODEState((jdouble, Vector{JComplex}, Vector{Vector{JComplex}}), arg0, arg1, arg2)
end

function get_complete_state(obj::ComplexODEState)
    return jcall(obj, "getCompleteState", Vector{JComplex}, ())
end

function get_complete_state_dimension(obj::ComplexODEState)
    return jcall(obj, "getCompleteStateDimension", jint, ())
end

function get_number_of_secondary_states(obj::ComplexODEState)
    return jcall(obj, "getNumberOfSecondaryStates", jint, ())
end

function get_primary_state(obj::ComplexODEState)
    return jcall(obj, "getPrimaryState", Vector{JComplex}, ())
end

function get_primary_state_dimension(obj::ComplexODEState)
    return jcall(obj, "getPrimaryStateDimension", jint, ())
end

function get_secondary_state(obj::ComplexODEState, arg0::jint)
    return jcall(obj, "getSecondaryState", Vector{JComplex}, (jint,), arg0)
end

function get_secondary_state_dimension(obj::ComplexODEState, arg0::jint)
    return jcall(obj, "getSecondaryStateDimension", jint, (jint,), arg0)
end

function get_time(obj::ComplexODEState)
    return jcall(obj, "getTime", jdouble, ())
end


text:
using AttractorNetworksTests
using Test

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


text:
# Render the widgets from WidgetsBase.jl
function jsrender(session::Session, button::Button)
    return jsrender(session, DOM.input(
        type = "button",
        value = button.content,
        onclick = js"JSServe.update_obs($(button.value), true);";
        button.attributes...
    ))
end

function jsrender(session::Session, tf::TextField)
    return jsrender(session, DOM.input(
        type = "textfield",
        value = tf.value,
        onchange = js"JSServe.update_obs($(tf.value),  this.value);";
        tf.attributes...
    ))
end

function jsrender(session::Session, ni::NumberInput)
    return jsrender(session, DOM.input(
        type = "number",
        value = ni.value,
        onchange = js"JSServe.update_obs($(ni.value), parseFloat(this.value));";
        ni.attributes...
    ))
end

function jsrender(session::Session, slider::Slider)
    return jsrender(session, DOM.input(
        type = "range",
        min = map(first, slider.range),
        max = map(last, slider.range),
        value = slider.value,
        step = map(step, slider.range),
        oninput = js"JSServe.update_obs($(slider.value), parseFloat(value))";
        slider.attributes...
    ))
end

function JSServe.jsrender(session::Session, tb::Checkbox)
    return jsrender(session, DOM.input(
        type = "checkbox",
        checked = tb.value,
        onchange = js"JSServe.update_obs($(tb.value), this.checked);";
        tb.attributes...
    ))
end

const noUiSlider = Dependency(
    :noUiSlider,
    [dependency_path("nouislider.min.js"), dependency_path("nouislider.min.css")]
)

function jsrender(session::Session, slider::RangeSlider)
    args = (slider.range, slider.connect, slider.orientation,
            slider.tooltips, slider.ticks)
    style = map(session, args...) do range, connect, orientation, tooltips, ticks
        value = slider.value[]
        return Dict(
            "range" => Dict(
                "min" => fill(range[1], length(value)),
                "max" => fill(range[end], length(value)),
            ),
            "step" => step(range),
            "start" => value,
            "connect" => connect,
            "tooltips" => tooltips,
            "direction" => "ltr",
            "orientation" => string(orientation),
            "pips" => isempty(ticks) ? nothing : ticks
        )
    end
    rangediv = DOM.div()
    create_slider = js"""function create_slider(style){
        var range = $(rangediv);
        range.noUiSlider.updateOptions(JSServe.deserialize_js(style), true);
    }"""
    onload(session, rangediv, js"""function onload(range){
        var style = $(style[]);
        $(noUiSlider).create(range, style);
        range.noUiSlider.on('update', function (values, handle, unencoded, tap, positions){
            JSServe.update_obs($(slider.value), [parseFloat(values[0]), parseFloat(values[1])]);
        });
    }""")
    onjs(session, style, create_slider)
    return rangediv
end

"""
A simple wrapper for types that conform to the Tables.jl Table interface,
which gets rendered nicely!
"""
struct Table
    table
    class::String
    row_renderer::Function
end

render_row_value(x) = x
render_row_value(x::Missing) = "n/a"
render_row_value(x::AbstractString) = string(x)
render_row_value(x::String) = x

function Table(table; class="", row_renderer=render_row_value)
    return Table(table, class, row_renderer)
end

function JSServe.jsrender(session::Session, table::Table)
    names = string.(Tables.schema(table.table).names)
    header = DOM.thead(DOM.tr(DOM.th.(names)...))
    rows = []
    for row in Tables.rows(table.table)
        push!(rows, DOM.tr(DOM.th.(table.row_renderer.(values(row)))...))
    end
    body = DOM.tbody(rows...)
    return DOM.table(header, body; class=table.class)
end

const ace = JSServe.Dependency(
    :ace,
    ["https://cdn.jsdelivr.net/gh/ajaxorg/ace-builds/src-min/ace.js"]
)

struct CodeEditor
    theme::String
    language::String
    options::Dict{Symbol, Any}
    onchange::Observable{String}
    element::Hyperscript.Node{Hyperscript.HTMLSVG}
end

"""
    CodeEditor(language::String; initial_source="", theme="chrome", editor_options...)

Defaults for `editor_options`:
```
(
    autoScrollEditorIntoView = true,
    copyWithEmptySelection = true,
    wrapBehavioursEnabled = true,
    useSoftTabs = true,
    enableMultiselect = true,
    showLineNumbers = false,
    fontSize = 16,
    wrap = 80,
    mergeUndoDeltas = "always"
)
```
"""
function CodeEditor(language::String; initial_source="", theme="chrome", editor_options...)
    defaults = Dict(
        :autoScrollEditorIntoView => true,
        :copyWithEmptySelection => true,
        :wrapBehavioursEnabled => true,
        :useSoftTabs => true,
        :enableMultiselect => true,
        :showLineNumbers => false,
        :fontSize => 16,
        :wrap => 80,
        :mergeUndoDeltas => "always"
    )
    options = Dict{Symbol, Any}(merge(Dict(editor_options), defaults))
    onchange = Observable(initial_source)
    element = DOM.div("", id="editor")
    return CodeEditor(theme, language, options, onchange, element)
end

function jsrender(session::Session, editor::CodeEditor)
    theme = "ace/theme/$(editor.theme)"
    language = "ace/mode/$(editor.language)"
    JSServe.onload(session, editor.element, js"""
        function (element){
            var editor = $ace.edit(element);
            editor.setTheme($theme);
            editor.session.setMode($language);
            editor.resize();
            // use setOptions method to set several options at once
            editor.setOptions($(editor.options));

            editor.session.on('change', function(delta) {
                JSServe.update_obs($(editor.onchange), editor.getValue());
            });

            editor.session.setValue($(editor.onchange[]));
        }
    """)
    return editor.element
end

# Ok, this is bad piracy, but I donno how else to make the display nice for now!
function Base.show(io::IO, m::MIME"text/html", widget::WidgetsBase.AbstractWidget)
    show(io, m, App(widget))
end


text:
# LTI primitives
#
# Procedure for efficient numerical evaluation of convolution, deconvolution, and related
# operators.

# FIXME: The `convolve` function should work even if `g` is shorter than `u`.

"""
    convolve(g, u)

Convolution of signals `g` and `u`.
"""
function convolve(g::Vector{T}, u::Vector{T}) where {T<:Number}
    @assert length(g) == length(u) "Romeo.LTI.convolve: Input signals must be of the same length."
    y = similar(u)
    if length(y) > 0
        @inbounds @simd for t = 1:length(u)
            val = zero(y[t])
            @simd for τ = 1:t
                val += g[t-τ+1] * u[τ]
            end # for
            y[t] = val
        end # for
    end # if
    return y
end # function

"""
    deconvolve(y, u)

Deconvolution of signals `y` and `u`.
Compute signal `g` such that `g ⋆ u = y`.
"""
function deconvolve(y::Vector{T}, u::Vector{T}) where {T<:Number}
    @assert length(y) == length(u) "Romeo.LTI.deconvolve: Input signals must be of the same length."
    g = Vector{typeof(one(T) / one(T))}(undef, size(u))
    if length(u) > 0
        m = one(T) / u[1]
        g[1] = y[1] * m
        @inbounds @simd for t = 2:length(u)
            val = y[t]
            @simd for τ = 2:t
                val -= g[t-τ+1] * u[τ]
            end # for
            g[t] = val * m
        end # for
    end # if
    return g
end # function

"""
    convroot(g)

Convolutional root of `h`.
Compute signal `g` such that `h = g ⋆ g`.

For signals having negative initial values, the convolutive root is multivalued
(similar to the algebraic root). The convolutive root will return the solution
obtained using the primary branch of the algebraic root.
"""
function convroot(h::Vector{T}) where {T<:Number}
    if length(h) == 0
        g = Vector{typeof(one(T) / one(T))}(undef, size(h))
    else
        @assert h[1] != 0 "Romeo.LTI.convroot: The initial value of the input signal cannot be zero."
        if T <: Real
            if h[1] > 0
                g = Vector{typeof(one(T) / one(T))}(undef, size(h))
                h1 = h[1]
            else
                g = Vector{typeof(one(Complex{T}) / one(Complex{T}))}(undef, size(h))
                h1 = complex(h[1])
            end
        else
            g = Vector{typeof(one(T) / one(T))}(undef, size(h))
            h1 = h[1]
        end
        if length(g) > 0
            g[1] = sqrt(h1)
            if length(g) >= 2
                @inbounds @simd for k = 2:length(h)
                    val = h[k]
                    @simd for i = 2:(k-1)
                        val -= g[k-i+1] * g[i]
                    end
                    g[k] = val / 2 / g[1]
                end
            end
        end
    end
    return g
end

"""
    convid([T=typeof(1.0)], N)

Convolutional identity operator.

Kernel length is `N`. Element type is determined by `T`.
"""
function convid(T::Type{<:Number}, N::Integer)
    δ = zeros(T, N)
    δ[1] = one(T)
    return δ
end # function

convid(N::Integer) = convid(typeof(1.0), N)

"""
    convinv(g)

Convolutional inverse of the signal `g`.
"""
function convinv(g::Vector{T}) where {T<:Number}
    id = convid(T, length(g))
    return deconvolve(id, g)
end # end

export convolve, deconvolve, convroot, convid, convinv