text:
mutable struct Befunge<:固体
	dire::UInt8
	coll::Bool
	error::Bool
	stack::Vector{Int}
	on::Union{Nothing,地毯}
	input::IO
	output::String
end
Befunge(dire::UInt8=0x1)=Befunge(dire,false,false,[-1],nothing,stdin,"")
function i_show(::Befunge,con::DrawContext)
	clear_rect(con,0,0;color=:slategray)
	so=get_loadedimg("Befunge0")
	fill_image(con,so,2,2,28,28)
end
function b_show(b::Befunge,con,x::Int,y::Int)
	if b.error
		so=get_loadedimg("Befunge1")
	else
		so::Matrix=get_loadedimg("Befunge0")
		if so.dire==0x2
			so=rotr90(so)
		elseif so.dire==0x3
			so=rot180(so)
		elseif so.dire==0x4
			so=rotr90(so)
		end
	end
	if b.on===nothing
		fill_image(con,so,x,y)
	else
		clear_rect(con,x,y;color=b.on.background_color)
		fill_image(con,so,x+2,y+2,28,28)
	end
end
function upd(b::Befunge,x::Int,y::Int)
	if b.on===nothing
		b.error=true
		return
	end
	tx=wtext(b.on::地毯)
	bs=b.stack
	if b.coll
		if tx=='"' b.coll=false
		else push!(bs,tx) end
		wmove(b,x,y)
		return
	end
	try
		tx=='+' ? begin fi=pop!(bs);se=pop!(bs);push!(bs,fi+se) end :
		tx=='-' ? begin fi=pop!(bs);se=pop!(bs);push!(bs,se-fi) end :
		tx=='*' ? begin fi=pop!(bs);se=pop!(bs);push!(bs,fi*se) end :
		tx=='/' ? begin fi=pop!(bs);se=pop!(bs);push!(bs,div(se,fi,RoundDown)) end :
		tx=='%' ? begin fi=pop!(bs);se=pop!(bs);push!(bs,mod(se,fi)) end :
		tx=='!' ? begin p=pop!(bs);push!(bs,p==0 ? 0 : 1) end :
		tx=='`' ? begin fi=pop!(bs);se=pop!(bs);push!(bs,se>fi ? 1 : 0) end :
		tx=='>' ? begin b.dire=0x2 end :
		tx=='<' ? begin b.dire=0x4 end :
		tx=='^' ? begin b.dire=0x1 end :
		tx=='v' ? begin b.dire=0x3 end :
		tx=='?' ? begin b.dire=crandom(1,4) end :
		tx=='_' ? begin p=pop!(bs);b.dire=p==0 ? 0x2 : 0x4 end :
		tx=='|' ? begin p=pop!(bs);b.dire=p==0 ? 0x3 : 0x1 end :
		# "
		tx==':' ? begin push!(bs,last(bs)) end :
		tx=='\\' ? begin fi=pop!(bs);se=pop!(bs);push!(bs.fi);push!(bs,se) end :
		tx=='$' ? begin pop!(bs) end :
		tx=='.' ? begin b.output*=string(pop!(bs)) end :
		tx==',' ? begin b.output*=Char(pop!(bs)) end :
		tx=='#' ? begin wmove(b,x,y);x+=direx[b.dire];y+=direy[b.dire] end :
		tx=='g' ? begin
			y=pop!(bs);x=pop!(bs)
			blk=getblk(x,y)
			if isa(blk,地毯) push!(bs,wtext(blk))
			else throw(true)
			end
		end :
		tx=='p' ? begin
			y=pop!(bs);x=pop!(bs);v=pop!(bs)
			blk=getblk(x,y)
			if isa(blk,地毯) blk.text=Char(v)
			else throw(true)
			end
		end :
		tx=='&' ? begin int=parse(Int,readline(b.input));push!(bs,int) end :
		tx=='~' ? begin char=read(b.input,Char);push!(bs,Int(char)) end :
		tx=='@' ? begin return end :
		'0'<=tx<='9' ? begin push!(bs,tx-'0') end : nothing
	catch
		b.error=true
		return
	end
	wmove(b,x,y)
end
function wmove(b::Befunge,x::Int,y::Int)
	x0=x+direx[b.dire]
	y0=y+direy[b.dire]
	tar=getblk(x0,y0)
	if isa(tar,地毯)
		setblk(x,y,b.on)
		b.on=tar
		setblk(x0,y0,b)
	else
		b.error=true
	end
end


text:
# 2D SLAM with velocity states

export   DynPose2, DynPose2VelocityPrior, PackedDynPose2VelocityPrior, DynPose2Pose2, PackedDynPose2Pose2

"""
$(TYPEDEF)
"""
mutable struct DynPose2VelocityPrior{T1,T2} <: IncrementalInference.AbstractPrior where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief}
  Zpose::T1
  Zvel::T2
  DynPose2VelocityPrior{T1,T2}() where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = new{T1,T2}()
  DynPose2VelocityPrior{T1,T2}(z1::T1,z2::T2) where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = new{T1,T2}(z1,z2)
end
DynPose2VelocityPrior(z1::T1,z2::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = DynPose2VelocityPrior{T1,T2}(z1,z2)

getSample(cf::CalcFactor{<:DynPose2VelocityPrior}, N::Int=1) = ([rand(cf.factor.Zpose,N);rand(cf.factor.Zvel,N)], )

function IIF.getParametricMeasurement(s::DynPose2VelocityPrior{<:MvNormal, <:MvNormal})

  meas = [mean(s.Zpose); mean(s.Zvel)]

  iΣp = invcov(s.Zpose)
  iΣv = invcov(s.Zvel)

  iΣ = zeros(eltype(iΣp), 5,5)

  iΣ[1:3,1:3] .= iΣp
  iΣ[4:5,4:5] .= iΣv

  return meas, iΣ
end

function (cf::CalcFactor{<:DynPose2VelocityPrior})(meas, X)
  #pose part, reused from PriorPose2
  iXihat = SE2(meas[1:3]) \ SE2(X[1:3])
  res_pose = se2vee(iXihat)

  #velocity part, normal prior
  res_vel = meas[4:5] .- X[4:5]

  return [res_pose; res_vel]
end

"""
$(TYPEDEF)
"""
mutable struct DynPose2Pose2{T <: IIF.SamplableBelief} <: IIF.AbstractRelativeRoots
  Zpose::Pose2Pose2{T} #Zpose::T1
  # reuseres::Vector{Float64}
  partial::Tuple{Int,Int,Int}
  DynPose2Pose2{T}() where {T <: IIF.SamplableBelief} = new{T}()
  DynPose2Pose2{T}(z1::T) where {T <: IIF.SamplableBelief} = new{T}(Pose2Pose2(z1), (1,2,3))
end
DynPose2Pose2(z1::T) where {T <: IIF.SamplableBelief} = DynPose2Pose2{T}(z1)

getSample(cf::CalcFactor{<:DynPose2Pose2}, N::Int=1) = (rand(cf.factor.Zpose.z,N), )

function (cf::CalcFactor{<:DynPose2Pose2})( meas,
                                            wXi,
                                            wXj  )
  #
  # cf.factor.Zpose(res, meas, wXi, wXj)
  wXjhat = SE2(wXi)*SE2(meas)
  jXjhat = SE2(wXj) \ wXjhat
  return se2vee(jXjhat)
  # z = meas[1][:,idx]
  # wxi, wxj = wXi[:,idx], wXj[:,idx]
  # dt = (userdata.variableuserdata[2].ut - userdata.variableuserdata[1].ut)*1e-6
  # wpj = ( wxi[1:2]+dt*wxi[4:5] + z[1:2] )
  # thetaj = se2vee(SE2([0;0;wxi[3]])*SE2([0;0;z[3]]))[3]
  # res[1:3] = se2vee( SE2(wxj[1:3])\SE2([wpj;thetaj]) )
  # # res[1:2] = (wXj[1:2] - (wXi[1:2]+dt*wXi[4:5])+z[1:2])
  # res[4:5] = z[4:5] - (wxj[4:5] -  wxi[4:5])
  # nothing
end

function compare(a::DynPose2VelocityPrior, b::DynPose2VelocityPrior)::Bool
  RoME.compareDensity(a.Zpose, b.Zpose) && RoME.compareDensity(a.Zvel, b.Zvel)
end

## Packing types

"""
$(TYPEDEF)
"""
mutable struct PackedDynPose2VelocityPrior <: IncrementalInference.PackedInferenceType
  strpose::AbstractString
  strvel::AbstractString
  PackedDynPose2VelocityPrior() = new()
  PackedDynPose2VelocityPrior(z1::AS, z2::AS) where {AS <: AbstractString} = new(z1, z2)
end

function convert(::Type{PackedDynPose2VelocityPrior}, d::DynPose2VelocityPrior)
  return PackedDynPose2VelocityPrior(convert(PackedSamplableBelief, d.Zpose),convert(PackedSamplableBelief, d.Zvel))
end
function convert(::Type{DynPose2VelocityPrior}, d::PackedDynPose2VelocityPrior)
  posedistr = convert(SamplableBelief, d.strpose)
  veldistr = convert(SamplableBelief, d.strvel)
  return DynPose2VelocityPrior(posedistr, veldistr)
end

"""
$(TYPEDEF)
"""
mutable struct PackedDynPose2Pose2 <: IncrementalInference.PackedInferenceType
  strpose::AbstractString
  PackedDynPose2Pose2() = new()
  PackedDynPose2Pose2(z1::AS) where {AS <: AbstractString} = new(z1)
end

function convert(::Type{PackedDynPose2Pose2}, d::DynPose2Pose2)
  return PackedDynPose2Pose2(string(d.Zpose.z))
end
function convert(::Type{DynPose2Pose2}, d::PackedDynPose2Pose2)
  posedistr = extractdistribution(d.strpose)
  return DynPose2Pose2(posedistr)
end

## ===========================================================================================

"""
$(TYPEDEF)
"""
mutable struct DynPose2DynPose2{T <: IIF.SamplableBelief} <: AbstractRelativeRoots
  Z::T
  reuseres::Vector{Vector{Float64}}
  DynPose2DynPose2{T}() where {T <: IIF.SamplableBelief} = new{T}()
  DynPose2DynPose2{T}(z1::T) where {T <: IIF.SamplableBelief} = new{T}(z1,[zeros(5) for i in 1:Threads.nthreads()])
end
DynPose2DynPose2(z1::T=MvNormal(zeros(5), diagm([0.01;0.01;0.001;0.1;0.1].^2))) where {T <: IIF.SamplableBelief} = DynPose2DynPose2{T}(z1)

# getManifolds(::Type{DynPose2DynPose2}) = (:Euclid,:Euclid,:Circular,:Euclid,:Euclid)
# getManifolds(::DynPose2DynPose2) = getManifolds(DynPose2DynPose2)

getSample(cf::CalcFactor{<:DynPose2DynPose2}, N::Int=1) = (rand(cf.factor.Z, N), )

function (cf::CalcFactor{<:DynPose2DynPose2})(meas,
                                              wXi,
                                              wXj  )
  #
  # vp2vp2.Zpose(res, userdata, idx, meas, Xi, Xj)
    # wXjhat = SE2(wxi[1:3,idx])*SE2(meas[1][1:3,idx])
    # jXjhat = SE2(wxj[1:3,idx]) \ wXjhat
    # se2vee!(res, jXjhat)
  z = meas
  wxi, wxj = wXi, wXj
  dt = Dates.value(cf.metadata.fullvariables[2].nstime - cf.metadata.fullvariables[1].nstime)*1e-9
  wpj = ( wxi[1:2]+dt*wxi[4:5] + z[1:2] )
  thetaj = se2vee(SE2([0;0;wxi[3]])*SE2([0;0;z[3]]))[3]
  res13 = se2vee( SE2(wxj[1:3])\SE2([wpj;thetaj]) )
  # res[1:2] = (wXj[1:2] - (wXi[1:2]+dt*wXi[4:5])+z[1:2])
  res45 = z[4:5] - (wxj[4:5] -  wxi[4:5])
  return [res13;res45]
end

function compare(a::DynPose2DynPose2, b::DynPose2DynPose2; tol::Float64=1e-10)::Bool
  TP = true
  TP = TP && RoME.compareDensity(a.Z, b.Z)
  TP = TP && norm(a.reuseres - b.reuseres) < tol
  return TP
end

"""
$(TYPEDEF)
"""
mutable struct PackedDynPose2DynPose2 <: IncrementalInference.PackedInferenceType
  Z::String
  PackedDynPose2DynPose2() = new()
  PackedDynPose2DynPose2(z1::AbstractString) = new(z1)
end

function convert(::Type{PackedDynPose2DynPose2}, d::DynPose2DynPose2)
  return PackedDynPose2DynPose2(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{DynPose2DynPose2}, d::PackedDynPose2DynPose2)
  return DynPose2DynPose2(convert(SamplableBelief, d.Z))
end

#

ext: jl


text:
using ConjugatePriors, Distributions
using ConjugatePriors: NormalInverseChisq, NormalInverseWishart
using Distributions: NormalStats, MvNormalStats
using Particles: add, sub, empty_suffstats

@testset "Component" begin
    @testset "add/sub from NormalStats" begin
        @test reduce(add, 1:10, init=NormalStats(0,0,0,0)) == suffstats(Normal, 1:10)
        @test reduce(sub, 1:5, init=suffstats(Normal, 1:10)) == suffstats(Normal, 6:10)
        @test sub(suffstats(Normal, [1]), 1) == NormalStats(0,0,0,0)
    end

    @testset "add/sub from MvNormalStats" begin
        Base.:(==)(a::MvNormalStats, b::MvNormalStats) =
            a.m == b.m &&
            a.s == b.s &&
            a.s2 == b.s2 &&
            a.tw == b.tw
        vecs = [[i,i+1] for i in 1.:10]
        @test reduce(add, vecs, init=MvNormalStats(zeros(2), zeros(2), zeros(2,2), 0)) ==
            suffstats(MvNormal, hcat(vecs...))
        @test reduce(sub, vecs[1:5], init=suffstats(MvNormal, hcat(vecs...))) ==
            suffstats(MvNormal, hcat(vecs[6:10]...))
        @test sub(suffstats(MvNormal, hcat(vecs[1])), vecs[1]) ==
            MvNormalStats(zeros(2), zeros(2), zeros(2,2), 0)
    end

    @testset "empty sufficient statistics" begin
        @test empty_suffstats(NormalInverseChisq()) == NormalStats(0, 0, 0, 0)
        ss1 = empty_suffstats(NormalInverseWishart(zeros(1), 1., ones(1,1), 1.))
        ss2 = MvNormalStats(zeros(1), zeros(1), zeros(1,1), 0.)
        @test all(getfield(ss1, n) == getfield(ss2, n) for n in fieldnames(typeof(ss1)))
        ss1 = empty_suffstats(NormalInverseWishart(zeros(3), 1., Matrix(1.0I,3,3), 1.))
        ss2 = MvNormalStats(zeros(3), zeros(3), zeros(3,3), 0.)
        @test all(getfield(ss1, n) == getfield(ss2, n) for n in fieldnames(typeof(ss1)))
    end

end


text:
"""
A Raster Tileset Provider typically follows the following conventions:
    * Tiles are 256 × 256 pixel PNG files
    * Each zoom level is a directory, each column is a subdirectory, and
        each tile in that column is a file
    * Filename(url) format is `/zoom/x/y.png`

The slippy map expects tiles to be served up at URLs following this scheme, so
all tile server URLs look pretty similar. (We currently only assume support for
raster tiles.)

Tiles are not always in these dimensions; for example there could be 64×64 pixel
images for mobile use, however 256×256 pixel images are a de facto standard.
515×512 pixel seems to be the usual size of high-resolution tiles.

Please obey the usage restrictions of each tile server if you use it heavily!

# Interface
In this interface, an AbstractProvider must implement the following method:
```
geturl(provider::AbstractProvider,x,y,z)
```
which returns a string representing the url to poll for tiles by that provider.

In addition, the AbstractProvider should contain the following attributes
    * `maxtiles`: maximum number of tiles to fetch, before issuing an error.

The AbstractProvider can override the default error message (when the max.
number of tiles is exceeded) by providing an implementation for
```
usagewarning(provider::AbstractProvider)
```

# Methods
`MapTiles.fetchtile(provider::AbstractProvider,x,y,z)` will return a
matrix of values corresponding to the raster tile at `/z/x/y.png`

# References
For more details see
* http://mapschool.io/
* http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
* http://wiki.openstreetmap.org/wiki/Tiles
* http://leaflet-extras.github.io/leaflet-providers/preview/
"""
abstract type AbstractProvider end

usagewarning(provider::AbstractProvider) =
    string("You are requesting a very large map, please consider ",
           "setting up your own tile server instead. For details, ",
           "please see https://switch2osm.org/serving-tiles/.")

function request(
        provider::AbstractProvider,
        x::Integer,
        y::Integer,
        z::Integer
    )
    url = geturl(provider,x,y,z)
    result = HTTP.request("GET", url)
    Logging.@info("Requesting $url (status: $(result.status))")
    result.body
end

function fetchrastertile(
        provider::AbstractProvider,
        x::Integer,
        y::Integer,
        z::Integer
    )
    data = request(provider, x, y, z)
    ImageMagick.readblob(data)
end

function fetchvectortile(
        provider::AbstractProvider,
        x::Integer,
        y::Integer,
        z::Integer
    )
    data = request(provider, x, y, z)
    ProtoBuf.readproto(IOBuffer(data), MapTiles.vector_tile.Tile())
end

"""
```
struct GBIFProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
    epsg::Int = 3857
    tileset::String = "omt"
    format::String = "pbf"
end
```

GBIF hosts map tiles in four projections, for use on GBIF.org and by users of
our API. We serve vector tiles based on the OpenMapTiles.org schema, and from
these generate raster tiles in several styles. Both the vector and raster tiles
are available for public use.

## Data schema and sources

The Web Mercator tiles are simply those from OpenMapTiles.org, see the
OpenMapTiles schema for details. The other projections have been generated to
fit the same layer schema, but with fewer layers — see https://tile.gbif.org/ui/
for details.

* Web Mercator tiles are from 2017-07-03.
* The other projections use OpenStreetMap and Natural Earth data. We do not plan
    to make frequent updates to the tiles, since the geological features of most
    interest to us do not often change.
* WGS84 tiles are from OpenStreetMap as of 2017-08-28, and Natural Earth version
    2.0.0.
* Polar tiles are from OpenStreetMap as of 2017-02-27, and Natural Earth version
    2.0.0. (They will soon be processed to add extra layers.)

## Vector Tiles

The schema for the vector tile layers follows the OpenMapTiles schema, but with
fewer layers included for the non-Mercator projections. In addition, the Arctic
LAEA projection includes a graticule layer, containing a polygon for the
Northern Hemisphere and lines for 60°N and the Arctic Circle.

## Format

* `format`: Either `.pbf` for a vector tile, or `@1x.png`, ..., `@4x.png` for
    PNG tiles at various pixel densities. Use higher pixel densities for high
    resolution screens, otherwise the tiles will appear fuzzy.
* `params`: Only for PNG tiles, specify a style e.g. `?style=gbif-classic`

## Feedback and issues

The process used by GBIF to generate the tiles in the three other projections
is on GitHub. Please log issues there, or use the GBIF API Users mailing list.

## Credits and copyright

* OpenStreetMap and Natural Earth for the map data,
* OpenMapTiles for the cartography and processing to produce vector tiles,
* OpenLayers for good support of non-Mercator projections.
* Base map tiles are © OpenMapTiles © OpenStreetMap contributors, and use
    Natural Earth data.

## References
https://tile.gbif.org/
https://www.gbif.org/developer/maps
"""
Parameters.@with_kw struct GBIFProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
    epsg::Int = 3857
    tileset::String = "omt"
    format::String = "pbf"
end
geturl(p::GBIFProvider, x::Integer, y::Integer, z::Integer) =
    "https://tile.gbif.org/$(p.epsg)/$(p.tileset)/$z/$x/$y.$(p.format)"

"""
```
struct OpenStreetMapProvider <: AbstractProvider
    variant::String = "standard"
    maxtiles::Int = 16
end
```
## Possible variants
    * `"standard"` (default)
    * `"blackwhite"`
    * `"deutschland"`
    * `"nolabels"`
    * `"france"`
    * `"humanitarian"`

## Where does the map data come from
OpenStreetMap is open data, created by OpenStreetMap contributors,
and available under the Open Data Commons Open Database License (ODbL).

## How do I fix an error on the map
https://www.openstreetmap.org/fixthemap

## Remarks
Apart from very limited testing purposes, you should not use the tiles supplied
by OpenStreetMap.org itself. OpenStreetMap is a volunteer-run non-profit body
and cannot supply tiles for large-scale commercial use. Rather, you should use
a third party provider that makes tiles from OSM data, or generate your own.

## References
See
    * https://switch2osm.org/using-tiles/
    * http://wiki.openstreetmap.org/wiki/Tile_servers
for details.
"""
Parameters.@with_kw struct OpenStreetMapProvider <: AbstractProvider
    variant::String = "standard"
    maxtiles::Int = 16
end

function geturl(provider::OpenStreetMapProvider, x::Integer, y::Integer, z::Integer)
    if provider.variant == "standard"
        "http://tile.openstreetmap.org/$z/$x/$y.png"
    elseif provider.variant == "blackwhite"
        "http://tiles.wmflabs.org/bw-mapnik/$z/$x/$y.png"
    elseif provider.variant == "nolabels"
        "https://tiles.wmflabs.org/osm-no-labels/$z/$x/$y.png"
    elseif provider.variant == "deutschland"
        "http://tile.openstreetmap.de/tiles/osmde/$z/$x/$y.png"
    elseif provider.variant == "france"
        "http://tile.openstreetmap.fr/osmfr/$z/$x/$y.png"
    elseif provider.variant == "humanitarian"
        "http://tile.openstreetmap.fr/hot/$z/$x/$y.png"
    else
        warn("Unknown provider: $(provider.variant), defaulting to standard")
        "http://tile.openstreetmap.org/$z/$x/$y.png"
    end
end
usagewarning(provider::OpenStreetMapProvider) =
    string("You are requesting a very large map, please see ",
           "OpenStreetMap tile usage policy ",
           "(https://operations.osmfoundation.org/policies/tiles/).")

"""
```
struct WikimediaMapsProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
end
```

## Terms of Use
Static images of maps can be used under a Creative Commons Attribution-ShareAlike
4.0 license, with attribution to https://wikimediafoundation.org/w/index.php?title=Maps_Terms_of_Use.

## Disclaimers
While we aim to provide a useful service for the Wikimedia community and beyond,
no software is perfect.

Wikimedia Maps service may undergo future development, so all your use is at
your sole risk. We provide Wikimedia Maps on an “as is” and “as available”
basis, and we expressly disclaim all warranties of all kinds, including implied
warranties of fitness for a particular purpose, merchantability, and
non-infringement. We make no warranty that Wikimedia Maps will meet your
requirements, be safe, secure, or uninterrupted.

## Where does the map data come from
The maps shown on Wikipedia, Wikivoyage and other Wikimedia projects use data
from OpenStreetMap. OpenStreetMap is open data, created by OpenStreetMap contributors,
and available under the Open Data Commons Open Database License (ODbL).

The map designs are based on the style OSM Bright for Mapbox Studio, available
under a Creative Commons Attribution 3.0 license.

## How do I fix an error on the map
https://www.openstreetmap.org/fixthemap

## References
See https://www.mediawiki.org/wiki/Maps for details.
"""
Parameters.@with_kw struct WikimediaMapsProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
end
geturl(provider::WikimediaMapsProvider, x::Integer, y::Integer, z::Integer) =
    "https://maps.wikimedia.org/osm-intl/$z/$x/$y.png"

"""
```
struct OpenCycleMapProvider <: AbstractProvider
    variant::String = "cycle"
    maxtiles::Int = typemax(Int)
end
```

The OpenCycleMap global cycling map is based on data from the OpenStreetMap
project. At low zoom levels it is intended for overviews of national cycling
networks; at higher zoom levels it should help with planning which streets to
cycle on, where you can park your bike and so on.

The map is updated every few days with the latest data from OpenStreetMap. News
about OpenCycleMap can be found in the OpenStreetMap archives on my blog.

## Possible variants
    * `cycle`
    * `transport`

## Credits
* Height Data is from NASA’s SRTM and is available in the public domain.
* Terrain Imagery is provided by True Marble by Unearthed Outdoors, LLC and is
    licensed under a Creative Commons Attribution 3.0 United States License.
* All other map data, including all cycling-related data, is from
    OpenStreetMap.org contributors and is licensed under the Open Data Commons
    Open Database License (ODbL).
* The map tiles are available for re-use under the Creative Commons Attribution
    Share-Alike 2.0 License. For full details see the full terms and conditions
    at http://thunderforest.com/terms/.

## References
See https://www.opencyclemap.org/docs/ for further details.
"""
Parameters.@with_kw struct OpenCycleMapProvider <: AbstractProvider
    variant::String = "cycle"
    maxtiles::Int = typemax(Int)
end
function geturl(provider::OpenCycleMapProvider, x::Integer, y::Integer, z::Integer)
    if provider.variant == "cycle"
        "http://tile.opencyclemap.org/cycle/$z/$x/$y.png"
    elseif provider.variant == "transport"
        "http://tile2.opencyclemap.org/transport/$z/$x/$y.png"
    else
        warn("Unknown provider: $(provider.variant), defaulting to cycle")
        "http://tile.opencyclemap.org/cycle/$z/$x/$y.png"
    end
end

"""
OpenSeaMap was created in 2009 in response to a great need for freely-accessible
seafaring maps. OpenSeaMap's goal is to add nautical and tourism information
that would interest sailors OSM, and to present it in a pleasing way. This
includes beacons, buoys and other seamarks, port information, repair shops, ship
supplies and much more, but also shops, restaurants and places of interest.
OpenSeaMap is part of OpenStreetMap and uses its database.

The basemap is rendered using the OpenStreetMap-data. This map is extended with
nautical data that is saved in the OSM-Database as well.

Visit http://openseamap.org/index.php and http://wiki.openseamap.org/wiki/Main_Page
for details.
"""
Parameters.@with_kw struct OpenSeaMapProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
end
geturl(provider::OpenSeaMapProvider, x::Integer, y::Integer, z::Integer) =
    "http://tiles.openseamap.org/seamark/$z/$x/$y.png"

"""
```
struct OpenTopoMapProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
end
```

OpenTopoMap is a project aiming at rendering topographic maps from OSM and SRTM
data. The map style should look familiar to Germans.

## Terms of Use
Please make sure to stick to the licence CC-BY-SA and always name opentopomap.org.
If you plan to use the tiles for bigger projects, please contact either

* http://www.openstreetmap.org/user/derstefan, or
* http://www.openstreetmap.org/user/mogstar

for more information.

## References
Visit http://wiki.openstreetmap.org/wiki/OpenTopoMap for further details.
"""
Parameters.@with_kw struct OpenTopoMapProvider <: AbstractProvider
    maxtiles::Int = typemax(Int)
end
geturl(provider::OpenTopoMapProvider, x::Integer, y::Integer, z::Integer) =
    "http://tile.opentopomap.org/$z/$x/$y.png"

"""
```
struct OpenWeatherMapProvider <: AbstractProvider
    apikey::String
    variant::String = "clouds"
    maxtiles::Int = typemax(Int)
end
```

OpenWeatherMap is an online service that provides weather data, including
current weather data, forecasts, and historical data to the developers of web
services and mobile applications. For data sources, it utilizes meteorological
broadcast services, raw data from airport weather stations, raw data from radar
stations, and raw data from other official weather stations. All data is processed
by OpenWeatherMap in a way that it attempts to provide accurate online weather
forecast data and weather maps, such as those for clouds or precipitation.

## Possible variants
    * `clouds` (default)
    * `clouds_cls`
    * `precipitation`
    * `precipitation_cls`
    * `rain`
    * `rain_cls`
    * `pressure`
    * `pressure_cntr`
    * `wind`
    * `temp`
    * `snow`

## Terms of Use
All data that OpenWeatherMap provides is used under the terms of a Creative
Commons Attribution-ShareAlike license. However, any access beyond 60 requests
of basic data per minute, as well as bulk downloading of weather readings and
forecasts, require a paid subscription.

## References
For further details, visit https://openweathermap.org/.
"""
Parameters.@with_kw struct OpenWeatherMapProvider <: AbstractProvider
    apikey::String
    variant::String = "clouds"
    maxtiles::Int = typemax(Int)
end
geturl(provider::OpenWeatherMapProvider, x::Integer, y::Integer, z::Integer) =
    "http://tile.openweathermap.org/map/$(provider.variant)/$z/$x/$y.png?appid=$(provider.apikey)"

"""
```
struct ThunderForestProvider <: AbstractProvider
    variant::String = "cycle"
    apikey::String
    maxtiles::Int = typemax(Int)
end
```
The Thunderforest Platform is the 4th major revision of the systems that powers
OpenCycleMap. With mapnik and postgis at its core, supported with mod_tile,
squid, and a host of other components, it provides robust, fast, reliable and
wonderful-looking maps.

## Possible variants
    * `"cycle"` (default)
    * `"transport"`
    * `"transport-dark"`
    * `"spinal-map"`
    * `"landscape"`
    * `"outdoors"`
    * `"pioneer"`
    * `"mobile-atlas"`
    * `"neighbourhood"`

## Attribution
Attribution must be given to both “Thunderforest” and “OpenStreetMap contributors”.
Users of your site/application must have a working links to www.thunderforest.com
and www.openstreetmap.org/copyright, for example:

    Maps © Thunderforest, Data © OpenStreetMap contributors

For media that doesn’t support links, such as posters, books and other printed
media, then you can provide attribution as follows:

    Maps © www.thunderforest.com, Data © www.osm.org/copyright

## References
For further details, visit https://www.thunderforest.com/docs/map-tiles-api/ and
https://www.thunderforest.com/terms/.
"""
Parameters.@with_kw struct ThunderForestProvider <: AbstractProvider
    variant::String = "cycle"
    apikey::String
    maxtiles::Int = typemax(Int)
end
geturl(provider::ThunderForestProvider, x::Integer, y::Integer, z::Integer) =
    "http://tile.thunderforest.com/$(provider.variant)/$z/$x/$y.png?apikey=$(provider.apikey)"

"""
```
struct OpenMapSurferProvider <: AbstractProvider
    variant::String = "roads"
    maxtiles::Int = typemax(Int)
end
```
OpenMapSurfer (http://korona.geog.uni-heidelberg.de/) - is a web map service
with a couple of maps rendered with MapSurfer.NET. This web map service is
developed by Maxim Rylov User:Runge and hosted by the GIScience (Geoinformatics)
Research Group, Heidelberg University.

## Possible variants
    * `"roads"` (default)
    * `"adminb"`
    * `"roadsg"`

## Terms of Use
Tiles that were rendered by this service can be used freely and without charge
by any individuals through the website http://korona.geog.uni-heidelberg.de

If you intend to use tiles from OpenMapSurfer services in your own applications
please contact them.

Commercial usage of the services provided by OpenMapSurfer does need approval!

OpenStreetMap data is available under the Open Database License.

Relief shading derived from CIAT-CSI SRTM (zoom levels 0-7): Users are
prohibited from any commercial, non-free resale, or redistribution without
explicit written permission from CIAT.

Original data by Jarvis A., H.I. Reuter, A. Nelson, E. Guevara, 2008, Hole-filled
seamless SRTM data V4, International Centre for Tropical Agriculture (CIAT),
available from http://srtm.csi.cgiar.org.

Relief shading derived from ASTER GDEM (zoom levels 8-18): ASTER GDEM is a
product of METI and NASA'.

Relief shading derived from ETOPO1 (zoom levels 0-4). Amante, C. and B. W. Eakins,
ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis.
NOAA Technical Memorandum NESDIS NGDC-24, 19 pp, March 2009.

## License
Content is available under Creative Commons Attribution-ShareAlike 2.0 license.

## References
For further details, visit http://wiki.openstreetmap.org/wiki/OpenMapSurfer.
"""
Parameters.@with_kw struct OpenMapSurferProvider <: AbstractProvider
    variant::String = "roads"
    maxtiles::Int = typemax(Int)
end
geturl(provider::OpenMapSurferProvider, x::Integer, y::Integer, z::Integer) =
    "http://korona.geog.uni-heidelberg.de/tiles/$(provider.variant)/x=$x&y=$y&z=$z"

"""
Hydda is a OpenStreetMap style, created for Swedish use, but usable in the rest
of the world too.

## Possible variants
    * `"full"` (default)
    * `"base"`
    * `"roads_and_labels"`

## Attribution
I don't know, maybe a reference to http://hydda.se or
https://github.com/karlwettin/tilemill-style-hydda?

I have to ask https://github.com/joakimfors and https://github.com/kodapan what they think.

Please visit http://openstreetmap.se/tjanster and
https://github.com/osmlab/editor-layer-index/issues/347 for further details.
"""
Parameters.@with_kw struct HyddaProvider <: AbstractProvider
    variant::String = "full"
    maxtiles::Int = typemax(Int)
end
geturl(provider::HyddaProvider, x::Integer, y::Integer, z::Integer) =
    "http://tile.openstreetmap.se/hydda/$(provider.variant)/$z/$x/$y.png"

"""
```
struct MapBoxProvider <: AbstractProvider
    accesstoken::String
    id::String = "streets"
    maxtiles::Int = typemax(Int)
end
```
The Mapbox Maps API supports reading raster tilesets, vector tilesets, and
Mapbox Editor project features. Tilesets can be retrieved as images, TileJSON,
or HTML slippy maps for embedding. Mapbox Editor project features can be
retrieved as GeoJSON or KML.

## Map IDs
The following map IDs are accessible to all accounts using a valid access token:

    * `streets` (default)
    * `light`
    * `dark`
    * `satellite`
    * `streets-satellite`
    * `wheatpaste`
    * `streets-basic`
    * `comic`
    * `outdoors`
    * `run-bike-hike`
    * `pencil`
    * `pirates`
    * `emerald`
    * `high-contrast`

## Access Token
Access to Mapbox web services? requires an access token that connects API
requests to your account. The requests in this package don't include an access
token: you will need to supply one by specifying the token.

## References
For further details, visit https://www.mapbox.com/api-documentation/.
"""
Parameters.@with_kw struct MapBoxProvider <: AbstractProvider
    accesstoken::String
    id::String = "streets"
    maxtiles::Int = typemax(Int)
end
geturl(provider::MapBoxProvider, x::Integer, y::Integer, z::Integer) =
    "https://api.tiles.mapbox.com/v4/mapbox.$(provider.id)/$z/$x/$y.png?access_token=$(provider.accesstoken)"

"""
```
StamenProvider <: AbstractProvider
    variant::String = "toner"
    maxtiles::Int = typemax(Int)
end
```

Stamen (stamen.com) is a San Francisco design and development studio focused on
data visualization and map-making. Stamen heavily uses OpenStreetMap data in
many of their map visualizations and it has provided three CC-BY OpenStreetMap
tilesets: Toner, Terrain, and Watercolor.

## Possible variants
    * `toner` (default)
    * `toner-background`
    * `toner-hybrid`
    * `toner-labels`
    * `toner-lite`
    * `watercolor`
    * `terrain`
    * `terrain-background`

## License
Except otherwise noted, each of these map tile sets are © Stamen Design, under a
Creative Commons Attribution (CC BY 3.0) license.

## Attribution
* For Toner and Terrain: Map tiles by Stamen Design, under CC BY 3.0.
    Data by OpenStreetMap, under ODbL.

* For Watercolor: Map tiles by Stamen Design, under CC BY 3.0.
    Data by OpenStreetMap, under CC BY SA.

* If you roll your own tiles from another source you will still need to credit
    “Map data by OpenStreetMap, under ODbL.”

## References
For further details, visit http://maps.stamen.com/.
"""
Parameters.@with_kw struct StamenProvider <: AbstractProvider
    variant::String = "toner"
    maxtiles::Int = typemax(Int)
end
geturl(provider::StamenProvider, x::Integer, y::Integer, z::Integer) =
    "http://tile.stamen.com/$(provider.variant)/$z/$x/$y.png"

"""
```
struct CARTOProvider <: AbstractProvider
    variant::String = "light_all"
    maxtiles::Int = typemax(Int)
end
```

## Possible variants
    * `light_all` (default)
    * `light_nolabels
    * `light_only_labels`
    * `dark_all`
    * `dark_nolabels`
    * `dark_only_labels`

## Remarks
For raster basemaps you can make use of them directly without an API KEY. For
vector basemaps you will need to get an API Key from Mapzen to start using them.
Remember that to use the basemaps you need to ensure you provide proper credit
and that in some cases you will need a commercial API KEY to use them.

## Terms of Service
Are based on https://carto.com/legal/, and detailed here:
https://drive.google.com/file/d/0B3OBExqwT6KJNHp3U3VUamx6U1U/view.

## References
For further details, visit https://carto.com/location-data-services/basemaps/ and
https://carto.com/docs/faqs/basemaps/.
"""
Parameters.@with_kw struct CARTOProvider <: AbstractProvider
    variant::String = "light_all"
    maxtiles::Int = typemax(Int)
end
geturl(provider::CARTOProvider, x::Integer, y::Integer, z::Integer) =
    "http://basemaps.cartocdn.com/$(provider.variant)/$z/$x/$y.png"

"""
```
GoogleMapsProvider <: AbstractProvider
    variant::String = ""
    maxtiles::Int = typemax(Int)
end
```

The Google Static Maps API lets you embed a Google Maps image on your web page
without requiring JavaScript or any dynamic page loading. The Google Static Maps
API service creates your map based on URL parameters sent through a standard
HTTP request and returns the map as an image you can display on your web page.

## Possible Variants
    * `roadmap` (default)
    * `terrain`
    * `satellite`
    * `hybrid`

For more details, visit https://developers.google.com/maps/documentation/static-maps/intro#MapTypes.

## Terms of Service
Please abide by the terms in https://developers.google.com/maps/terms?hl=en#section_10_5

## References
For further details, visit https://developers.google.com/maps/documentation/static-maps/
"""
Parameters.@with_kw struct GoogleMapsProvider <: AbstractProvider
    variant::String = ""
    maxtiles::Int = typemax(Int)
end
function geturl(provider::GoogleMapsProvider, x::Integer, y::Integer, z::Integer)
    layer = if provider.variant == ""
        ""
    elseif provider.variant == "roads-only"
        "lyrs=h"
    elseif provider.variant == "roadmap"
        "lyrs=m"
    elseif provider.variant == "terrain"
        "lyrs=p"
    elseif provider.variant == "altered-roadmap"
        "lyrs=r"
    elseif provider.variant == "satellite"
        "lyrs=s"
    elseif provider.variant == "terrain-only"
        "lyrs=t"
    elseif provider.variant == "hybrid"
        "lyrs=y"
    else
        warn("unknown map type: $(provider.variant), defaulting to standard")
        ""
    end
    "http://mt1.google.com/vt/$layer&x=$x&y=$y&z=$z"
end

"""
```
struct NASAGIBSProvider <: AbstractProvider
    variant::String = "MODIS_Terra_CorrectedReflectance_TrueColor"
    maxtiles::Int = typemax(Int)
end
```

The Global Imagery Browse Services (GIBS) are designed to deliver global,
full-resolution satellite imagery to users in a highly responsive manner,
enabling interactive exploration of the Earth.

## Imagery Layers & Endpoints
GIBS imagery layers are named and made available through a set of defined
endpoints based on the following characteristics:

* **Projection & Resolution**: Imagery layers are available in one or more
    projected coordinate systems (e.g. EPSG:4326 - "Geographic Lat/Lon") at a
    specific resolution (e.g. 2km/pixel)
* **Near Real-Time vs Standard Latency**: Imagery layers are available in a near
    real-time (e.g. within 3 hours of observation) or standard (e.g. within X
    days of observation) latency.
* **Data Version**: Imagery layers may be available for more than one version
    (e.g. MODIS v5 and v6) of the same science parameter.

## Populating the fields

https://gibs.earthdata.nasa.gov/wmts/epsg{EPSG:Code}/best/{ProductName}/default/{Time}/{TileMatrixSet}/{ZoomLevel}/{TileRow}/{TileCol}.png

with the desired projection, product, time, etc (Terra/MODIS Aerosol Optical
depth from 2014/04/09, in this case), GIBS products can be used by clients such
as ESRI's ArcGIS Online to add a "Tile Layer" by leaving the row, column, and
zoom level as parameters:

https://gibs.earthdata.nasa.gov/wmts/epsg3857/best/MODIS_Terra_Aerosol/default/2014-04-09/GoogleMapsCompatible_Level6/{level}/{row}/{col}.png

The time parameter supports a single day (YYYY-MM-DD) or repeating interval
(Rx/YYYY-MM-DD/PyYmMdD), as specified by the ISO 8601 specification, where:

    * x - number of repetitions (day frames)
    * YYYY - Year
    * MM - Month
    * DD - Day
    * y - Period years
    * m - period months
    * d - period days

## Possible variants
    * `MODIS_Terra_CorrectedReflectance_TrueColor` (default)
    * `VIIRS_CityLights_2012`
    * `MODIS_Terra_Land_Surface_Temp_Day`
    * `MODIS_Terra_Snow_Cover`
    * `MODIS_Terra_Aerosol`
    * `MODIS_Terra_Chlorophyll_A`

## References
Please visit https://wiki.earthdata.nasa.gov/display/GIBS/GIBS+API+for+Developers
for further details.
"""
Parameters.@with_kw struct NASAGIBSProvider <: AbstractProvider
    variant::String = "MODIS_Terra_CorrectedReflectance_TrueColor"
    time::String = ""
    tilematrixset::String = "GoogleMapsCompatible_Level"
    maxtiles::Int = 150_000
end
function geturl(provider::NASAGIBSProvider, x::Integer, y::Integer, z::Integer)
    ext = if provider.variant in ("MODIS_Terra_Land_Surface_Temp_Day", "MODIS_Terra_Snow_Cover", "MODIS_Terra_Aerosol", "MODIS_Terra_Chlorophyll_A")
        "png"
    else
        "jpg"
    end
    maxzoom = if provider.variant == "MODIS_Terra_Aerosol"
        6
    elseif provider.variant in ("MODIS_Terra_Chlorophyll_A", "MODIS_Terra_Land_Surface_Temp_Day")
        7
    elseif provider.variant in ("VIIRS_CityLights_2012", "MODIS_Terra_Snow_Cover")
        8
    else
        9
    end
    "https://map1.vis.earthdata.nasa.gov/wmts-webmerc/$(provider.variant)/default/$(provider.time)/$(provider.tilematrixset)$maxzoom/$z/$y/$x.$ext"
end
usagewarning(provider::NASAGIBSProvider) =
    string("You are requesting a very large map, please see the Bulk Downloading section of ",
           "https://wiki.earthdata.nasa.gov/display/GIBS/GIBS+API+for+Developers",
           " for details. Prior to beginning your bulk downloading activities, please ",
           "contact the GIBS support team at support@earthdata.nasa.gov.")

"""
Possible variants:
    * `World_Street_Map` (default)
    * `Specialty/DeLorme_World_Base_Map`
    * `World_Topo_Map`
    * `World_Imagery`
    * `World_Terrain_Base`
    * `World_Shaded_Relief`
    * `World_Physical_Map`
    * `Ocean_Basemap`
    * `NatGeo_World_Map`
    * `Canvas/World_Light_Gray_Base`
"""
Parameters.@with_kw struct ESRIProvider <: AbstractProvider
    variant::String = "World_Street_Map"
    maxtiles::Int = typemax(Int)
end
geturl(provider::ESRIProvider, x::Integer, y::Integer, z::Integer) =
    "http://server.arcgisonline.com/ArcGIS/rest/services/$(provider.variant)/MapServer/tile/$z/$x/$y"


text:
module test_82_AllExportedFunctions

import ModiaResult
ModiaResult.@usingModiaPlot

# Test the Figure operations
println("\n... test test_82_AllExportedFunctions.jl:\n")

t = range(0.0, stop=10.0, length=100)
result = Dict{String,Any}("time" => t, "phi" => sin.(t))
info   = ModiaResult.resultInfo(result)
ModiaResult.printResultInfo(result)

plot(result, "phi", figure=2)

showFigure(2)
saveFigure(2, "test_saveFigure.png")
closeFigure(2)
closeAllFigures()

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
"""
    train!(loss, params, data, opt; cb)
For each datapoint `d` in `data` computes the gradient of `loss(p,d...)` through
backpropagation and calls the optimizer `opt`.
Takes a callback as keyword argument `cb`. For example, this will print "training"
every 10 seconds:
```julia
DiffEqFlux.sciml_train(loss, params, data, opt,
            cb = throttle(() -> println("training"), 10))
```
The callback can call `Flux.stop()` to interrupt the training loop.
Multiple optimisers and callbacks can be passed to `opt` and `cb` as arrays.
"""
function sciml_train(loss, _θ, opt; cb = (args...) -> (false), maxiters)
  θ = copy(_θ)
  ps = Flux.params(θ)
  data = Iterators.repeated((), maxiters)
  t0 = time()
  # Flux is silly and doesn't have an abstract type on its optimizers, so assume
  # this is a Flux optimizer
  local x
  @progress for d in data
    gs = Flux.Zygote.gradient(ps) do
      x = loss(θ)
      first(x)
    end
    Flux.Optimise.update!(opt, ps, gs)
    cb_call = cb(θ,x...)
    if !(typeof(cb_call) <: Bool)
      error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
    elseif cb_call
      break
    end
  end
  _time = time()
  Optim.MultivariateOptimizationResults(opt,
                                        _θ,# initial_x,
                                        θ, #pick_best_x(f_incr_pick, state),
                                        first(x), # pick_best_f(f_incr_pick, state, d),
                                        maxiters, #iteration,
                                        maxiters >= maxiters, #iteration == options.iterations,
                                        false, # x_converged,
                                        0.0,#T(options.x_tol),
                                        0.0,#T(options.x_tol),
                                        NaN,# x_abschange(state),
                                        NaN,# x_abschange(state),
                                        false,# f_converged,
                                        0.0,#T(options.f_tol),
                                        0.0,#T(options.f_tol),
                                        NaN,#f_abschange(d, state),
                                        NaN,#f_abschange(d, state),
                                        false,#g_converged,
                                        0.0,#T(options.g_tol),
                                        NaN,#g_residual(d),
                                        false, #f_increased,
                                        nothing,
                                        maxiters,
                                        maxiters,
                                        0,
                                        true,
                                        NaN,
                                        _time-t0,)
end

decompose_trace(trace::Optim.OptimizationTrace) = last(trace)
decompose_trace(trace) = trace

function sciml_train(loss, θ, opt::Optim.AbstractOptimizer;
                      cb = (args...) -> (false), maxiters = 0, kwargs...)
  local x
  function _cb(trace)
    cb_call = cb(decompose_trace(trace).metadata["x"],x...)
    if !(typeof(cb_call) <: Bool)
      error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
    end
    cb_call
  end

  function optim_loss(θ)
    x = loss(θ)
    first(x)
  end

  function optim_loss_gradient!(g,θ)
    g .= Flux.Zygote.gradient(optim_loss,θ)[1]
    nothing
  end
  optimize(optim_loss, optim_loss_gradient!, θ, opt,
           Optim.Options(extended_trace=true,callback = _cb,
                         f_calls_limit = maxiters; kwargs...))
end


text:
const Ord = Base.Order.ForwardOrdering
const ObsData = SortedDict{Tuple, OrderedDict{String,Observable}, Ord}

mutable struct ObsTable
    data::ObsData
    par_names::OrderedSet{String}
end

## TODO deprecate, no table should be initialized
ObsTable() = ObsTable(ObsData(Ord()), OrderedSet{String}())

function ObsTable(::Type{Params}) where Params
    t = ObsTable(string.(fieldnames(Params)))
    t
end
ObsTable(params::T) where {T} = ObsTable(T)
ObsTable(pnames::Vector) = ObsTable(string.(pnames))
ObsTable(pnames::Tuple) = ObsTable(collect(pnames))

function ObsTable(pnames::Vector{String})
    t = ObsTable()
    set_params_names!(t, pnames)
    t
end

==(t1::ObsTable, t2::ObsTable) = t1.data==t2.data && t1.par_names==t2.par_names

#set_params_names!(t::ObsTable, a) = set_params_names!(fieldnames(a))
set_params_names!(t::ObsTable, keys...) = set_params_names!(t, collect(keys))
set_params_names!(t::ObsTable, keys::Vector{Symbol}) = set_params_names!(t, string.(keys))
set_params_names!(t::ObsTable, keys::Vector) = error("non valid keys")

function set_params_names!(t::ObsTable, keys::Vector{<:AbstractString})
    length(t.par_names) != 0 && Error("Can be done only once!")
    t.par_names = OrderedSet(String.(keys))
end

params_names(t::ObsTable) = t.par_names

function obs_names(t::ObsTable)
    names = OrderedSet{String}()
    for (k,v) in t.data
        push!(names, keys(v)...)
    end
    names
end

function nsamples(t::ObsTable, par::Tuple)
    haskey(t, par) || return 0
    return first(t[par])[2].t
end

# Indexing inteface
getindex(t::ObsTable, i...) = get!(t.data, to_index(t, i), OrderedDict{String,Observable}())
setindex!(t::ObsTable, v, i...) = setindex!(t.data, v, to_index(t, i))

getindex(t::ObsTable, i) = get!(t.data, to_index(t, i), OrderedDict{String,Observable}())
setindex!(t::ObsTable, v, i) = setindex!(t.data, v, to_index(t, i))

splat(a) = [getfield(a,f) for f in fieldnames(typeof(a))]

function to_index(t, i::Tuple)
    @assert length(params_names(t)) == length(i) "wrong number of params in indexing"
    i
end

to_index(t, i::Vector) = to_index(t, (i...,))
function to_index(t, i)
    if length(fieldnames(typeof(i))) == 0 #it is not a composite type
        return to_index(t, (i,))
    else
        return to_index(t, (splat(i)...,))
    end
end

# SOME type piracy here
getindex(d::OrderedDict{String,Observable}, i::String) = get!(d, i, Observable())
getindex(d::OrderedDict{String,Observable}, i::Symbol) = get!(d, string(i), Observable())
setindex!(d::OrderedDict{String,Observable}, x, i::Symbol) = setindex!(d, x, string(i))

endof(t::ObsTable) = endof(t.data)

haskey(t::ObsTable, k) = haskey(t.data, k)

# Iterable inteface
Base.iterate(t::ObsTable, x...) = iterate(t.data, x...)
eltype(::Type{ObsTable}) = eltype(typeof(t.data))
length(t::ObsTable) = length(t.data)

"""
    header(t::ObsTable)

Returns a string containg the header of the printed table.
"""
function header(t::ObsTable; lenpar=9, lenobs=18)
    pnames= t.par_names
    onames = obs_names(t)
    h = ""
    i = 1
    for k in pnames
        s = i == 1 ? "# $i:$k " : "$i:$k "
        h *= rpad(s, lenpar)
        i+=1
    end
    h *= rpad("$i:num ", lenpar)
    i+= 1
    for k in onames
        h *=  rpad("$i:$k ", lenobs)
        i +=2
    end
    return h
end

function Base.show(io::IO, t::ObsTable)
    lenpar = 9
    lenobs = 21
    println(io, header(t, lenpar=lenpar, lenobs=lenobs))
    onames = obs_names(t)
    for (par, obs) in t
        for p in par
            print(io, rpad("$p ", lenpar))
        end

        print(io, rpad("$(nsamples(t, par)) ", lenpar))
        for name in onames
            s = haskey(obs, name) ? "$(obs[name]) " : "NaN NaN "
            print(io, rpad(s, lenobs))
        end
        println(io)
    end
end

"""
    merge(t::ObsTables, ts::ObsTable...)
    merge!(t::ObsTables, ts::ObsTable...)

Fuse togheter the observations of many `ObsTable`s
as if the where independent measurements.
"""
function merge(t1::ObsTable, ts::ObsTable...)
    t = deepcopy(t1)
    return merge!(t, ts...)
end

function merge!(t::ObsTable, ts::ObsTable...)
    for t2 in ts
        t.par_names = union(t.par_names, t2.par_names)
        for (k,v) in t2
            for (n, o) in v
                t[k][n] &= o
            end
        end
    end
    return t
end

"""
    tomatrices(t::ObsTable) -> (params, means, errs)

Converts `t` into three matrices, `params`,`means` and `errs`.
"""
function tomatrices(t::ObsTable)
    n = length(t)
    m1 = length(params_names(t))
    pars = zeros(n,m1)
    onames = obs_names(t)
    m2 = length(obs_names(t))
    y = zeros(n, m2)
    yerr = zeros(n, m2)
    i = 1
    for (p, olist) in t
        pars[i,:] = [p...]
        for (j, oname) in enumerate(onames)
            obs = olist[oname]
            y[i,j] = mean(obs)
            yerr[i,j] = err(obs)
        end
        i += 1
    end
    return pars, y, yerr
end

"""
    tomatrix(t::ObsTable)

Converts `t` into a matrix corresponding to the printed table.
See [`header`](@ref) for the meaning of the columns.
"""
function tomatrix(t::ObsTable)
    n = length(t)
    m1 = length(params_names(t))
    onames = obs_names(t)
    m2 = length(obs_names(t))
    y = zeros(n, m1 + 1 + 2m2) #condidering also #samples
    i = 1
    for (p, olist) in t
        y[i,1:m1] = [p...]
        y[i,m1+1] = nsamples(t, p)
        for (j, oname) in enumerate(onames)
            obs = olist[oname]
            j = 2j+m1
            y[i,j] = mean(obs)
            y[i,j+1] = err(obs)
        end
        i += 1
    end
    return y
end

"""
    ObsTable(datfile::String)

Read an observable from a properly formatted data file
(i.e. the result of a `print(datfile,t::ObsTable)`).
"""
function ObsTable(datfile::String)
    f = open(datfile,"r")
    t = ObsTable()
    eof(f) && return t
    header = readline(f)
    header = split(header)
    @assert header[1] == "#"

    colnums = parse.(Int, map(x->split(x,':')[1], header[2:end]))
    names = map(x->split(x,':')[2], header[2:end])

    count = findlast(colnums .== 1:length(colnums))
    nparams = count - 2
    @assert names[nparams+1] ∈ ["num", "samples", "nsamples", "nsamp"]
    set_params_names!(t, names[1:nparams]...)
    eof(f) && return t
    res = readdlm(f)
    close(f)
    onames = String.(names[count:end])
    for i=1:size(res,1)
        pars = (res[i,1:nparams]...,)
        nsamp = res[i,nparams+1]
        for j=count:2:size(res,2)-1
            m, e = res[i,j], res[i,j+1]
            name = onames[(j-count)÷2 + 1]
            if isfinite(m)
                t[pars][name] = obs_from_mean_err_samp(m, e, nsamp)
            end
        end
    end
    return t
end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
function play(cur_player::AbstractPlayer)
  if cur_player.game.in_focus
    cur_player.state.is_playing || return
  else
    cur_player.state.is_playing = true
  end

  cur_js = ""

  if cur_player.state.is_done
    cur_player.state.is_done = false

    cur_js *= """
      \$(".js-tetris-container td").removeClass();
      \$(".js-level-text").text("Level $(cur_player.level)");
      \$(".js-score-text").text("$(lpad(cur_player.score, 8, "0"))");

      setTimeout(
        function() {
          \$(".js-splash-text h1").text("tetris.");
          \$(".js-splash-text div").text("double-click to play.");
        },
        500
      );
    """
  end

  cur_game = cur_player.game

  cur_game.in_focus = true

  drop_clock(cur_player)

  cur_js *= """
    \$(".js-tetris-container").removeClass("cs-has-splash");
    \$(".js-splash-screen").addClass("cs-disappear");
  """

  tetris_js(
    cur_game.scope,
    JSString(cur_js)
  )

  return
end


text:
## BIOCLIM model functions

# import SimpleSDMLayers: bioclim

# Small internal function for BIOCLIM. Ensures values are equal for both tails and between 0-1
_bcscore(x) = isnothing(x) ? nothing : (x > 0.5 ? 2 * (1 - x) : 2x)

# 1st part of BIOCLIM model for single environmental layer
function _bioclim_layer(
    occ::Union{G,D}, layer::T; training_layer::T=layer
) where {T<:SimpleSDMLayer,D<:DataFrame,G<:GBIFRecords}
    # occ: occurences of a single species as a DataFrame with latitude and longitude columns
    # layer: layer of environmental variables used for prediction
    # training_layer: optional, training environmental layer, can use a different resolution

    # Get observed environmental values (training values)
    obs = training_layer[occ]
    filter!(!isnothing, obs)
    # Create empty layer for predictions
    pred = similar(layer)
    # Make predictions
    if length(unique(obs)) <= 1 # special case if 1 unique value only
        # If 1 unique value only, set 1.0 as prediction for all sites with this value
        pred.grid = replace(x -> x == obs[1] ? 1.0 : 0.0, layer.grid)
    else
        # Create ECDF function to extract quantile values
        qfinder = ecdf(obs)
        # Get prediction score
        qfgrid = replace(x -> isnothing(x) ? nothing : qfinder(x), layer.grid)
        pred.grid = _bcscore.(qfgrid)
    end
    # Restore nothings from original layer
    for idx in findall(isnothing, layer.grid)
        pred.grid[idx] = nothing
    end
    return pred
end

# Complete BIOCLIM model on all environmental layers, including 1st part function
function bioclim(occ, layers; training_layers=layers, binary=true, threshold::Float64=0.0)
    # occ: occurences of a single species as a DataFrame with latitude and longitude columns
    # layers: layers of environmental variables used for prediction
    # training_layers: optional, training environmental layers, can use a different resolution
    # binary: optional, convert result to binary presence-absence values
    # with_threshold: optional, apply threshold to remove lower predictions
    # threshold: optional, percentage of lower predictions to remove (e.g. 0.05)

    @assert 0.0 ≤ threshold < 1.0 "Threshold must be between 0.0 and 1.0"

    # Apply 1st part of BIOCLIM on each environmental variable
    predictions = [
        _bioclim_layer(occ, layers[i]; training_layer=training_layers[i]) for
        i in eachindex(layers)
    ]
    # Reduce to single layer with minimum values
    prediction = reduce(min, predictions)
    # Apply threshold (if specified)
    if threshold > 0.0
        # Select non-nothing values only
        no_nothing = filter(!isnothing, prediction[occ])
        # Apply threshold only if there are non-nothing values
        if length(no_nothing) != 0
            # Get prediction value equivalent to threshold
            threshold_value = first(quantile(no_nothing, [threshold]))
            # Replace values smaller than threshold by nothing
            replace!(
                x -> !isnothing(x) || x < threshold_value ? nothing : x, prediction.grid
            )
        end
    end
    # Replace zeros by nothing
    replace!(x -> x == 0.0 ? nothing : x, prediction.grid)
    # Convert to binary presence-absence values (if requested, default is true)
    if binary == true
        replace!(x -> x > 0.0 ? 1.0 : x, prediction.grid)
    end
    return prediction
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames <: SwaggerModel
    categories::Any # spec type: Union{ Nothing, Vector{String} } # spec name: categories
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    listKind::Any # spec type: Union{ Nothing, String } # spec name: listKind
    plural::Any # spec type: Union{ Nothing, String } # spec name: plural
    shortNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: shortNames
    singular::Any # spec type: Union{ Nothing, String } # spec name: singular

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames(;categories=nothing, kind=nothing, listKind=nothing, plural=nothing, shortNames=nothing, singular=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames, Symbol("categories"), categories)
        setfield!(o, Symbol("categories"), categories)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames, Symbol("listKind"), listKind)
        setfield!(o, Symbol("listKind"), listKind)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames, Symbol("plural"), plural)
        setfield!(o, Symbol("plural"), plural)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames, Symbol("shortNames"), shortNames)
        setfield!(o, Symbol("shortNames"), shortNames)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames, Symbol("singular"), singular)
        setfield!(o, Symbol("singular"), singular)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames = Dict{Symbol,Symbol}(Symbol("categories")=>Symbol("categories"), Symbol("kind")=>Symbol("kind"), Symbol("listKind")=>Symbol("listKind"), Symbol("plural")=>Symbol("plural"), Symbol("shortNames")=>Symbol("shortNames"), Symbol("singular")=>Symbol("singular"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames = Dict{Symbol,String}(Symbol("categories")=>"Vector{String}", Symbol("kind")=>"String", Symbol("listKind")=>"String", Symbol("plural")=>"String", Symbol("shortNames")=>"Vector{String}", Symbol("singular")=>"String")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    (getproperty(o, Symbol("plural")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionNames }, name::Symbol, val)
end


text:
# Degrade resolution of maps generated by "gen_pysm_maps.py"
using Healpix
using Statistics
using Printf
using CFITSIO
ν = [27, 39, 93, 145, 225, 280, 350, 410, 850]  # in GHz
nν = length(ν)
nside = 128
for iν = 1:nν
    infile = @sprintf("map_equ_%03dghz_r8_nested_uKcmb.fits", ν[iν])
    fiin, fqin, fuin = (readMapFromFITS(infile, ic, Float32) for ic = 1:3)
    # Degrade and convert to RING order
    q = HealpixMap{Float32,RingOrder}(nside)
    u = HealpixMap{Float32,RingOrder}(nside)
    fact = Int((fiin.resolution.nside / nside)^2)
    for ip = 1:12*nside^2
        degrade = (ip-1)*fact+1:ip*fact
        jp = nest2ring(q.resolution, ip)
        q[jp] = mean(fqin[degrade])
        u[jp] = mean(fuin[degrade])
    end
    outfile = @sprintf("map_d1s1_equ_%03dghz_r7_uKcmb_qu.fits", ν[iν])
    extname = "Archive Map Table"
    typechar = "E"
    f = CFITSIO.fits_create_file(outfile)
    CFITSIO.fits_create_binary_tbl(
    f,
    0,
    [
        ("Stokes Q", "1$typechar", "uK, thermodynamic"),
        ("Stokes U", "1$typechar", "uK, thermodynamic"),
    ],
    extname,
    )
    saveToFITS(q, f, 1)
    saveToFITS(u, f, 2)
    CFITSIO.fits_close_file(f)
end


text:
# cpu info
# executed from command line
# julia> is the prompt

# cpu summary
julia> Sys.cpu_summary
cpu_summary (generic function with 3 methods)

julia> Sys.cpu_summary()
Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz:
       speed         user         nice          sys         idle          irq
#1  2496 MHz   70265296            0     52962390    1452947656      2625906  ticks
#2  2496 MHz   62610218            0     74936468    1438628203     22954593  ticks
#3  2496 MHz   72306906            0     52867046    1451000937       942187  ticks
#4  2496 MHz   67773812            0     43695703    1464705359      1063671  ticks

# cpu info
julia> Sys.cpu_info()
4-element Array{Base.Sys.CPUinfo,1}:
 Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz:
        speed         user         nice          sys         idle          irq
     2496 MHz   70256484            0     52955265    1452763812      2625656  ticks
 Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz:
        speed         user         nice          sys         idle          irq
     2496 MHz   62601796            0     74927593    1438445718     22952421  ticks
 Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz:
        speed         user         nice          sys         idle          irq
     2496 MHz   72298156            0     52859765    1450817187       942046  ticks
 Intel(R) Core(TM) i5-6300U CPU @ 2.40GHz:
        speed         user         nice          sys         idle          irq
     2496 MHz   67762500            0     43689281    1464523312      1063468  ticks

julia> Sys.CPU_THREADS
4

julia> Sys.CPU_NAME
"skylake"

julia> Sys.uptime()
7.2779536887234e6

julia> Sys.free_memory()
0x00000001d4ba9000

julia> Sys.ARCH
:x86_64

julia> Sys.MACHINE
"x86_64-w64-mingw32"

julia> Sys.total_memory()
0x00000002f7b02000

julia> Sys.KERNEL
:NT

julia> Sys.iswindows()
true

julia> Sys.islinux()
false

julia> Sys.isunix()
false

julia> Sys.loadavg()
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia>


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using RCall
using CSV
using DataFrames
using Dates
using Plots
using Underscores
using StatsBase
using ColorSchemes
plotlyjs()
## ========================================================================

@rlibrary var"data.table"
@rlibrary EpiNow2
@rlibrary covidregionaldata

import NetworkOptions
ENV["JULIA_SSL_CA_ROOTS_PATH"] = NetworkOptions.bundled_ca_roots()

begin
    R"""
    library(data.table)
    library(EpiNow2)
    library(covidregionaldata)
    country <- "france"
    options(mc.cores = 4)
    generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani")
    incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer")
    reporting_delay <- list(mean = convert_to_logmean(4, 1),
                            mean_sd = 0.1,
                            sd = convert_to_logsd(4, 1),
                            sd_sd = 0.1,
                            max = 15)
    # reported_cases <- read.csv(file = 'WHO-COVID-19-global-data.csv') #covidregionaldata::get_national_data(country)
    # reported_cases <- data.table::setDT(reported_cases)
    # reported_cases <- reported_cases[, .(date, confirm = cases_new)]
    """
end

regions = ["Hungary", "The United Kingdom", "France", "Sweden", "Denmark", "Germany", "Italy", "United States of America"]
old2new = Dict(
    "Hungary" => "Hungary",
    "The United Kingdom" => "UK",
    "France" => "France",
    "Sweden" => "Sweden",
    "Denmark" => "Denmark",
    "Germany" => "Germany",
    "Italy" => "Italy",
    "United States of America" => "USA"
)

function load_data(regions)
    fname = projectdir("data/epinow/WHO-COVID-19-global-data.csv")
    who = CSV.read(fname, DataFrame)
    observations = Dict{String, DataFrame}()
    for region in regions
        regional = @_ who |>
                filter(:Country => ==(region), __) |>
                select(__, :Date_reported => :date, :New_cases => :cases, :New_deaths => :deaths) #|>
        #__[Date("2020-08-01") <. __.date <. Date("2021-02-01") ]
        observations[old2new[region]] = regional
    end
    observations
end

obs = load_data(regions)

fname = projectdir("data/epinow/WHO-COVID-19-global-data.csv")
who = CSV.read(fname, DataFrame)
## ========================================================================
# for region in regions
function firefox(p; fname="tmp.html", fdir=normpath(homedir(), ".tmp"))
    # mkpath(dir)
    filepath = normpath(fdir, fname)
    savefig(p, filepath)
    run(`firefox $(filepath)`, wait=false)
    return p
end

function plot_regions(contacts, observations, mobility)
    regions = keys(contacts)
    region = first(regions)
    clr = palette(:tab10)
    for region in keys(contacts)
        o = observations[region]
        c = contacts[region]
        m = mobility[region]
        p1 = plot(o.date, o.cases, lab="cases", lw=3, xticks=:native, yticks=:native, c=clr[1])
        p2 = plot(o.date, o.deaths, lab="deaths", lw=3, xticks=:native, yticks=:native, c=clr[2])
        p3 = plot(c.date, c.total |> zscore, lab="total contacts", lw=3, xticks=:native, yticks=:native, c=clr[3])
        plot!(p3, m.date, m.apple |> zscore, lab="apple", lw=3, xticks=:native, yticks=:native, c=clr[4])
        plot!(p3, m.date, m.google |> zscore, lab="google", lw=3, xticks=:native, yticks=:native, c=clr[5])
        p = plot(p1, p2, p3, layout=(3,1), size=(1000, 3*250), sharex=true, link=:x)
        firefox(p, fname="$region.html")
    end
end

plot_regions(
    BSON.load(projectdir("data/contacts/int/mean_contact_rates_aggregation=14.bson")),
    obs,
    BSON.load(projectdir("data/mobility/", "mobility_international.bson"))
    )
## ========================================================================
fname = projectdir("data/epinow/rt.csv")
regions = ["Hungary", "United Kingdom", "France", "Sweden", "Denmark", "Germany", "Italy", "United States"]
old2new = Dict(
    "Hungary" => "Hungary",
    "United Kingdom" => "UK",
    "France" => "France",
    "Sweden" => "Sweden",
    "Denmark" => "Denmark",
    "Germany" => "Germany",
    "Italy" => "Italy",
    "United States" => "USA"
)

raw_rt = CSV.read(fname, DataFrame)
region = "Hungary"
rt = @_ raw_rt |>
        filter(:country => ==(region), __) |>
        select(__, :date, :mean)

all_regions = unique(raw_rt.country)

let
    c = contacts[:Hungary]
    plot(c.date, c.total ./ c.nrow |> zscore)
    plot!(rt.date, rt.mean)
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
using Test
using SparseArrays
import sbp

m = 2 << 5
n = 2 << 4

@testset "read/write vector" begin
        file = "x.bin"

        for type in [Int32, Float32, Float64]
                x = ones(type, n)
                sbp.lmv.write_vector(file, x, verbose=true)
                sbp.lmv.write_vector(file, x, verbose=false)
                y = sbp.lmv.read_vector(type, file)
                @test typeof(y) == Array{type, 1}
                @test isapprox(x, y)
        end

        rm(file)
end

@testset "read/write CSR matrix" begin
        file = "A.bin"
        A = spzeros(m, n)
        A[1,1] = 1.0
        A[1,2] = 2.0
        A[1,3] = 3.0
        A[2,1] = 4.0
        A[2,2] = 5.0
        A[2,3] = 6.0
        sbp.lmv.write_csr_matrix(file, A, verbose=true)
        sbp.lmv.write_csr_matrix(file, A, verbose=false)
        B = sbp.lmv.read_csr_matrix(Float64, Int64, file, verbose=true)
        B = sbp.lmv.read_csr_matrix(Float64, Int64, file, verbose=false)
        @test typeof(B) == SparseMatrixCSC{Float64, Int64}
        @test isapprox(A, B)
        rm(file)
end

@testset "read/write dense matrix" begin
        file = "A.bin"
        A = zeros(n, n)
        A[1,1] = 1.0
        A[1,2] = 2.0
        A[1,3] = 3.0
        A[2,1] = 4.0
        A[2,2] = 5.0
        A[2,3] = 6.0
        sbp.lmv.write_dense_matrix(file, A, verbose=true)
        sbp.lmv.write_dense_matrix(file, A, verbose=false)
        B = sbp.lmv.read_dense_matrix(file, verbose=true)
        B = sbp.lmv.read_dense_matrix(file, verbose=false)
        @test typeof(B) == Array{Float64, 2}
        @test isapprox(A, B)

        rm(file)
end

@testset "read/write config" begin
        file = "config.txt"
        cfg = Dict("nt" => Int64(1), "dt" => Float64(0.1))
        sbp.lmv.write_config(file, cfg, verbose=false, skip_types=true)
        sbp.lmv.write_config(file, cfg, verbose=false)
        sbp.lmv.write_config(file, cfg, verbose=true)
        cfg2 = sbp.lmv.read_config(file, verbose=false)
        cfg2 = sbp.lmv.read_config(file, verbose=true)
        @test cfg == cfg2

        rm(file)

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
# GCM-specific Sources
# This file contains helpers and lists currently available options

# Current options for GCM-specific sources:

"""
    HeldSuarezForcing <: TendencyDef{Source}

Defines a forcing that parametrises radiative and frictional effects using
Newtonian relaxation and Rayleigh friction, following Held and Suarez (1994)
"""
struct HeldSuarezForcing <: TendencyDef{Source} end

prognostic_vars(::HeldSuarezForcing) = (Momentum(), Energy())

function held_suarez_forcing_coefficients(bl, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed
    FT = eltype(state)

    # Parameters
    T_ref = FT(255)

    param_set = parameter_set(bl)
    _R_d = FT(R_d(param_set))
    _day = FT(day(param_set))
    _grav = FT(grav(param_set))
    _cp_d = FT(cp_d(param_set))
    _p0 = FT(MSLP(param_set))

    # Held-Suarez parameters
    k_a = FT(1 / (40 * _day))
    k_f = FT(1 / _day)
    k_s = FT(1 / (4 * _day))
    ΔT_y = FT(60)
    Δθ_z = FT(10)
    T_equator = FT(315)
    T_min = FT(200)
    σ_b = FT(7 / 10)

    # Held-Suarez forcing
    φ = latitude(bl, aux)
    p = air_pressure(ts)

    #TODO: replace _p0 with dynamic surface pressure in Δσ calculations to account
    #for topography, but leave unchanged for calculations of σ involved in T_equil
    σ = p / _p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor
    return (k_v = k_v, k_T = k_T, T_equil = T_equil)
end

function source(::Energy, s::HeldSuarezForcing, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    nt = held_suarez_forcing_coefficients(m, args)
    FT = eltype(state)
    param_set = parameter_set(bl)
    _cv_d = FT(cv_d(param_set))
    @unpack k_T, T_equil = nt
    T = air_temperature(ts)
    return -k_T * state.ρ * _cv_d * (T - T_equil)
end

function source(::Momentum, s::HeldSuarezForcing, m, args)
    nt = held_suarez_forcing_coefficients(m, args)
    return -nt.k_v * projection_tangential(m, args.aux, args.state.ρu)
end


text:
# Update_alloc.jl

export confusion_matrixζ;

function update_ζ!(y::Vector{T}, state::State_GPMTD,
    mixcomp_type::MixComponentNormal) where T <: Real

    n = length(y)
    L = length(state.mixcomps)
    lp = Vector{Float64}(undef, L+1)

    for i = 1:n

        lp[1] = state.lλ[1] - 0.5*log(state.intercept.σ2) -
            0.5 * (y[i] - state.intercept.μ)^2 / state.intercept.σ2

        for j in 1:L
            lp[j+1] = state.lλ[j+1] -
                0.5*log(state.mixcomps[j].σ2) -
                0.5 * (y[i] - state.mixcomps[j].μ - state.mixcomps[j].fx[i])^2 / state.mixcomps[j].σ2
        end

        w = exp.( lp .- maximum(lp) )
        state.ζ[i] = StatsBase.sample(Weights(w)) - 1 # so that the values are 0:L

    end

    for ℓ = 1:L
        state.mixcomps[ℓ].ζon_indx = findall( state.ζ .== ℓ )
    end

    return nothing
end

function confusion_matrixζ(ζtrue::Vector{Int}, ζest::Vector{Int}, L::Int)

    n = length(ζtrue)
    n2 = length(ζest)

    n == n2 || throw("Two vectors must have same length.")

    out = zeros(Int, L+1, L+1)

    for i = 1:(L+1)
        for j = 1:(L+1)
            out[i,j] = sum( [ζtrue[k] == (i-1) && ζest[k] == (j-1) for k = 1:n] )
        end
    end

    out
end


text:
using LinearAlgebra, MixedModels, StableRNGs, Test, SparseArrays

include("modelcache.jl")

@testset "Xymat" begin
    trm = MixedModels.FeTerm(hcat(ones(30), repeat(0:9, outer = 3)), ["(Intercept)", "U"])
    piv = trm.piv
    ipiv = invperm(piv)
    mat = MixedModels.FeMat(trm, Float64.(collect(axes(trm.x, 1))))
    @test size(mat) == (30, 3)
    @test length(mat) == 90
    @test size(mat') == (3, 30)
    @test eltype(mat) == Float64
    @test mat.xy === mat.wtxy
    prd = mat'mat
    @test typeof(prd) == Matrix{Float64}
    @test prd[ipiv, ipiv] == [30.0 135.0; 135.0 855.0]
    wts = rand(StableRNG(123454321), 30)
    MixedModels.reweight!(mat, wts)
    @test mul!(prd, mat', mat)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)

    # empty fixed effects
    trm = MixedModels.FeTerm(ones(10,0), String[])
    #@test size(trm) == (10, 0)  # There no longer are size and length methods for FeTerm
    #@test length(trm) == 0
    #@test size(trm') == (0, 10)
    @test eltype(trm) == Float64
    @test trm.rank == 0
end

@testset "XymatSparse" begin

    @testset "sparse and dense yield same fit" begin
        # deepcopy because we're going to modify
        m = last(models(:insteval))
        # this is kinda sparse:
        # julia> mean(first(m.feterm).x)
        # 0.10040140325753434

        fe = m.feterm
        X =  MixedModels.FeTerm(SparseMatrixCSC(fe.x), fe.cnames)
        @test typeof(X.x) <: SparseMatrixCSC
        @test X.rank == 28
        @test X.cnames == fe.cnames
        m1 = fit!(LinearMixedModel(collect(m.y), X, deepcopy(m.reterms), m.formula); progress=false)
        @test isapprox(m1.θ, m.θ, rtol = 1.0e-5)
    end

    @testset "rank deficiency in sparse FeTerm" begin
        trm = MixedModels.FeTerm(SparseMatrixCSC(hcat(ones(30),
                                                     repeat(0:9, outer = 3),
                                                     2repeat(0:9, outer = 3))),
                                ["(Intercept)", "U", "V"])
        # at present there is no attempt to evaluate the rank of a SparseMatrixCSC
        piv = trm.piv
        ipiv = invperm(piv)
        @test_broken rank(trm) == 2
        mat = MixedModels.FeMat(trm, Float64.(collect(axes(trm.x, 1))))
        @test size(mat) == (30, 4)
        @test length(mat) == 120
        @test size(mat') == (4, 30)
        @test eltype(mat) == Float64
        @test mat.xy === mat.wtxy
        prd = MixedModels.densify(mat'mat)
        @test typeof(prd) == typeof(prd)
        @test prd[ipiv, ipiv] == [30.0 135.0 270.0; 135.0 855.0 1710.0; 270.0 1710.0 3420.0]
        wts = rand(StableRNG(123454321), 30)
        MixedModels.reweight!(mat, wts)
        @test mul!(prd, mat', mat)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
    end

end


text:
using Test
using Pkg
Pkg.develop(path = joinpath(pwd(), "..", "..", "OmegaTest"))
Pkg.develop(path = joinpath(pwd(), "..", "..", "connectors", "OmegaDistributions"))

@testset "alltests" begin
  include("typeinfer.jl")
  include("ciid.jl")
  include("condition.jl")
  include("distributions.jl")
  include("core.jl")
  include("intervene.jl")
  include("higherintervene.jl")
  include("logpdf.jl")
  include("multivariate.jl")
  include("namedtuple.jl")
  # include("solution.jl")
  include("tagging.jl")
  include("typeinfer.jl")
  include("pointwise.jl")
  # include("var.jl")
  include("mem.jl")
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Documenter, HomotopyContinuation
import LinearAlgebra

makedocs(
    sitename = "Homotopy Continuation",
    pages = [
        "Introduction" => "index.md",
        "Solving general systems" => "solving.md",
        "Solving paremeterized systems with monodromy" => "monodromy.md",
        "Sorting arrays of solutions" => "sorting.md",
        "PathTracker" => "path_tracker.md",
        "CoreTracker" => "core_tracker.md",
        "Newton's method" => "newton.md",
        "Homotopies" => "homotopies.md",
        "Predictors and correctors" => "predictors-correctors.md",
        "Data structures for polynomial systems" => "systems.md",
        "Norms and distances" => "norms_distances.md",
        "Reference" => "reference.md"
        ],
    strict=true)

deploydocs(
    repo = "github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl.git"
)


text:
using CaimCore.ImageStack
using Colors, FixedPointNumbers

const STACKS_DIR = joinpath(@__DIR__, "data", "stacks")
const VIDEOS_DIR = joinpath(@__DIR__, "data", "videos")

@testset "Loading Frames Fails" begin
    @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "colormix"))
    @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "sizemix"))
    #  @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "unknown"))
    #  @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "unloadable"))
    @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "nonimage"))
    @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "fmtmix"))
    @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "fmtmix"); ext=".jpg")

    @test_throws ErrorException load(Frames, joinpath(STACKS_DIR, "withdots");
                                     ignoredots=false)

    @test_throws ErrorException load(Frames, String[])

    @test_throws ErrorException load(Frames, joinpath(VIDEOS_DIR, "empty.avi"))
    @test_throws ErrorException load(Frames, joinpath(VIDEOS_DIR, "garbage.avi"))
    @test_throws ErrorException load(Frames, joinpath(VIDEOS_DIR, "nonvideo.txt"))
end

@testset "Loading Frames from a Directory" begin
    let baseframe = repeat([RGB{N0f8}(1., 0., 0.) RGB{N0f8}(0., 1., 0.);
                            RGB{N0f8}(0., 0., 0.) RGB{N0f8}(0., 0., 1)];
                           inner=(10,10))
        frames = cat(CaimCore.scanl(rotr90, baseframe, 4)...; dims=3)

        let stack = load(Frames, joinpath(STACKS_DIR, "png"))
            @test size(stack) == size(frames)
            @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(STACKS_DIR, "tiff"))
            @test size(stack) == size(frames)
            @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(STACKS_DIR, "withdots"))
            @test size(stack) == size(frames)
            @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(STACKS_DIR, "withdots"); ext=".tiff",
                         ignoredots=false)
            @test size(stack) == size(frames)
            @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end
    end

    let baseframe = repeat([RGB{N0f8}(1., 0., 0.) RGB{N0f8}(0., 1., 0.);
                            RGB{N0f8}(0., 0., 0.) RGB{N0f8}(0., 0., 1)];
                           inner=(10,10))
        frames = cat(CaimCore.scanl(rot180, baseframe, 1)...; dims=3)

        stack = load(Frames, joinpath(STACKS_DIR, "fmtmix"); ext=".png")
        @test size(stack) == size(frames)
        @test ImageStack.frames(stack) == frames
        @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×2"
    end

    let baseframe = repeat([Gray{N0f8}(0.) Gray{N0f8}(0.33); Gray{N0f8}(1.) Gray{N0f8}(0.66)];
                           inner=(10,10))
        frames = cat(CaimCore.scanl(rotr90, baseframe, 4)...; dims=3)

        stack = load(Frames, joinpath(STACKS_DIR, "grayscale"))
        @test size(stack) == size(frames)
        @test ImageStack.frames(stack) == frames
        @test string(stack) == "Frames{Gray{N0f8}}\n  size: 20×20×5"
    end
end
#
@testset "Loading Frames from Video File" begin
    let baseframe = repeat([RGB{N0f8}(1., 0., 0.) RGB{N0f8}(0., 1., 0.);
                            RGB{N0f8}(0., 0., 0.) RGB{N0f8}(0., 0., 1)];
                           inner=(10,10))
        frames = cat(CaimCore.scanl(rotr90, baseframe, 4)...; dims=3)

        let stack = load(Frames, joinpath(VIDEOS_DIR, "png.avi"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "png.mov"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "png.mp4"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "tiff.avi"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "tiff.mov"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "tiff.mp4"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end
    end

    let baseframe = repeat([Gray{N0f8}(0.) Gray{N0f8}(0.33); Gray{N0f8}(1.) Gray{N0f8}(0.66)];
                           inner=(10,10))
        frames = cat(CaimCore.scanl(rotr90, baseframe, 4)...; dims=3)

        let stack = load(Frames, joinpath(VIDEOS_DIR, "grayscale.avi"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "grayscale.mov"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end

        let stack = load(Frames, joinpath(VIDEOS_DIR, "grayscale.mp4"))
            @test size(stack) == size(frames)
            #  @test ImageStack.frames(stack) == frames
            @test string(stack) == "Frames{RGB{N0f8}}\n  size: 20×20×5"
        end
    end
end


text:
using JuliaLSODA, Test

@testset "LSODA Example" begin
    function fex(du, u, p, t)
        du[1] = 1e4 * u[2] * u[3] - 0.04e0 * u[1]
        du[3] = 3e7 * u[2] * u[2]
        du[2] = - (du[1] + du[3])
    end

    prob = ODEProblem(fex, [1.0, 0, 0], [0.0,0.4e0])
    sol2 = solve(prob, LSODA())
    reference = [9.8517230200179651e-01, 3.3863992446169744e-05, 1.4793834005757239e-02]
    @test  sol2.u[end] ≈ reference rtol = 1e-15
end

@testset "generic_type" begin
    function fex(du, u, p, t)
        du[1] = 1e4 * u[2] * u[3] - 0.04e0 * u[1]
        du[3] = 3e7 * u[2] * u[2]
        du[2] = - (du[1] + du[3])
    end

    prob = ODEProblem(fex, Float32[1.0, 0, 0], Float64[0.0,0.4e0])
    sol2 = solve(prob, LSODA(),saveat=0.1)
    @show typeof(sol2.t)
    @show sol2
    reference = Float32[9.8517230200179651e-01, 3.3863992446169744e-05, 1.4793834005757239e-02]
    @test  sol2.u[end] ≈ reference rtol = 1e-5
end

@testset "LSODA Example tsaveat" begin
    function fex(du, u, p, t)
        du[1] = 1e4 * u[2] * u[3] - 0.04e0 * u[1]
        du[3] = 3e7 * u[2] * u[2]
        du[2] = - (du[1] + du[3])
    end

    prob = ODEProblem(fex, Float64[1.0, 0, 0], Float32[0.0,0.4e0])
    sol2 = solve(prob, LSODA(),saveat = [0.02,0.2,0.3], atol = 1e-5, rtol = 1e-3)
    reference = [9.8517230200179651e-01, 3.3863992446169744e-05, 1.4793834005757239e-02]
    @show sol2
    #@test  sol2 ≈ reference rtol = 1e-15
end
@testset "LSODA Example" begin
    function fex(du, u, p, t)
        du[1] = 1e4 * u[2] * u[3] - 0.04e0 * u[1]
        du[3] = 3e7 * u[2] * u[2]
        du[2] = - (du[1] + du[3])
    end

    prob = ODEProblem(fex, [1.0, 0, 0], [0.0,0.4e0])
    sol2 = solve(prob, LSODA(),saveat=[0.02,0.2,0.3], tstops=[0.1,0.25])
    reference = [9.8517230200179651e-01, 3.3863992446169744e-05, 1.4793834005757239e-02]
    @show sol2
end
#=
@testset "Rober" begin
    function rober(du,u,p,t)
      y₁,y₂,y₃ = u
      k₁,k₂,k₃ = p
      du[1] = -k₁*y₁+k₃*y₂*y₃
      du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
      du[3] =  k₂*y₂^2
      nothing
    end
    prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
    sol2 = solve(prob, LSODA())
end

@testset "Linear ODE" begin
    function f(du,u,p,t)
        du[1] = 1.01*u[1]
    end
    u0=[1/2]
    tspan = (0.0,1.0)
    prob = ODEProblem(f,u0,tspan)
end

@testset "Pendulum" begin
    l = 1.0                             # length [m]
    m = 1.0                             # mass[m]
    g = 9.81                            # gravitational acceleration [m/s²]

    function pendulum!(du,u,p,t)
        du[1] = u[2]                    # θ'(t) = ω(t)
        du[2] = -3g/(2l)*sin(u[1]) + 3/(m*l^2)*p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
    end

    θ₀ = 0.01                           # initial angular deflection [rad]
    ω₀ = 0.0                            # initial angular velocity [rad/s]
    u₀ = [θ₀, ω₀]                       # initial state vector
    tspan = (0.0,10.0)                  # time interval

    M = t->0.1sin(t)                    # external torque [Nm]

    prob = ODEProblem(pendulum!,u₀,tspan,M)
    sol = solve(prob, LSODA())
end
=#


text:
import Base.+;

#Strategies Implementation

const Portfolio = Dict{String, Position};

#Link an underlying name to a Position
function →(x::String, y::Position)
    out = Portfolio(x => y)

    return out
end

#Link an underlying name to a SingleNamePayoff
function →(x::String, y::SingleNamePayoff)
    sep = findfirst("_", x)
    if (!isnothing(sep))
        error("NO UNDERSCORE ALLOWED IN SINGLE NAME OPTIONS")
    end
    return x → (1.0 * y)
end

#Link an underlying name to a BasketPayoff
function →(x::String, y::BasketPayoff)
    sep = findfirst("_", x)[1]
    if (isnothing(sep))
        error("Bivariate payoff underlying must follow the format: INDEX1_INDEX2")
    end
    return x → (1.0 * y)
end

#Join Portfolio s
function +(x::Portfolio, y::Portfolio)
    out = copy(x)
    y_keys = keys(y)
    for y_key in y_keys
        if haskey(out, y_key)
            out[y_key] = out[y_key] + y[y_key]
        else
            out[y_key] = y[y_key]
        end
    end
    return out
end

#Overload for display
import Base.Multimedia.display;
function display(p::Portfolio)
    keys_ = keys(p)
    for key_ in keys_
        print(key_)
        print(" => ")
        print(p[key_])
    end
    println("")
end

#Overload for extraction, output needs particular shape.
function extract_(x::String, dict_::Portfolio)
    keys_ = keys(dict_)
    out = Position()
    out_vec = Position[]
    ap_keys = String[]
    for key_ in keys_
        if ((x == key_) || (any(z -> z == key_, split(x, "_"))))
            out = dict_[key_]
            push!(out_vec, out)
            push!(ap_keys, key_)
        end
    end

    if (length(out_vec) == 1) && (x == ap_keys[1])
        return out
    else
        str_v = split(x, "_")
        X = compute_indices(length(str_v))

        #tmp_map=Dict( str_v .=> collect(1:length(str_v)))
        tmp_map_rev = Dict(collect(1:length(str_v)) .=> str_v)

        out__ = Array{Position}(undef, length(X))

        for i_ = 1:length(X)
            idx_ = X[i_]
            str_underlv = [tmp_map_rev[idx_2] for idx_2 in idx_]
            str_underl = join(str_underlv, "_")

            if (haskey(dict_, str_underl))
                out__[i_] = dict_[str_underl]
            end
        end

        return out__
    end
end


text:
using CartesianGP
using Base.Test
include("../src/Mutate.jl")

# Utility functions

function mutate_test_parameters()
    mu = 1
    lambda = 4
    numperlevel = 2
    numlevels = 8
    numlevelsback = 8
    numinputs = 3
    numoutputs = 2
    mutrate = 0.3
    targetfitness = 1.0
    funcs = [AND, OR, XOR, NOT, ONE, ZERO ]
    fitfunc = hamming_max

    p =  Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
    return p
end

function parameters_half_adder()
    numperlevel = 1
    numlevels = 2
    numlevelsback = numlevels
    numinputs = 2
    numoutputs = 2

    mu = 1
    lambda = 4
    mutrate = 0.15
    targetfitness = 1.0
    funcs = [AND, OR, XOR]
    fitfunc = hamming_max

    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end

# This version has one node per level
function chromosome_half_adder(numperlevel=1, numlevels=2, numlevelsback=2)
    numinputs = 2
    numoutputs = 2

    #p = Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)
    p = parameters_half_adder()
    c = Chromosome(p)

    c.inputs = [InputNode(1), InputNode(2)]
    c.interiors = Array(InteriorNode,2,1)
    c.interiors[1,1] = InteriorNode(XOR, [(0, 1), (0, 2)])
    c.interiors[2,1] = InteriorNode(AND, [(0, 1), (0, 2)])
    c.outputs = [OutputNode((1, 1)), OutputNode((2, 1))]

    return c
end

function goal_half_adder()
    g = Goal(2, (0b0110, 0b1000))
    return g
end

c0 = chromosome_half_adder()
cm0 = mutate(c0)
#print_chromosome(c0)
#print_chromosome(cm0)

num_iter = 20
rlist = [rand(1:1000) for i in 1:num_iter]
p = mutate_test_parameters()
for i in 1:num_iter
    r = rlist[i]
    #println("r:",r)   # If there is a problem, print r and run with it as a seed
    srand(r)
    c0 = random_chromosome(p)
    cm0 = mutate(c0)
    #print_chromosome(c0)
    #print_chromosome(cm0)
end


text:
# tests for problems and ensemble problem types
@testset "Unitary X gate pulse in place w/ GRAPE" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = Uinit,
        Xt = Ufin,
        T = 1.0,
        n_controls = 2,
        guess = rand(2, 10),
        sys_type = UnitaryGate(),
    )

    sol = solve(prob, GRAPE(n_slices = 10, isinplace = true))

    @test sol.result.minimum - C1(Ufin, Ufin) < tol
end

@testset "Unitary X gate pulse static arrays" begin

    prob = Problem(
        B = [SSx, SSy],
        A = SSz,
        Xi = UinitS,
        Xt = UfinS,
        T = 1.0,
        n_controls = 2,
        guess = rand(2, 10),
        sys_type = UnitaryGate(),
    )

    sol = solve(prob, GRAPE(n_slices = 10, isinplace = false))
    @test sol.result.minimum - C1(UfinS, UfinS) < tol

end

@testset "Robust X gate pulse inplace" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = Uinit,
        Xt = Ufin,
        T = 5.0,
        n_controls = 2,
        guess = rand(2, 100),
        sys_type = UnitaryGate(),
    )

    ens = EnsembleProblem(
        prob = prob,
        n_ens = 5,
        A_g = A_gens,
        B_g = B_gens,
        XiG = U_init_gens,
        XtG = U_target_gens,
        wts = ones(5) / 5,
    )

    sol = solve(
        ens,
        GRAPE(
            n_slices = 100,
            isinplace = true,
            optim_options = Optim.Options(f_tol = 1e-3),
        ),
    )
    @test sol.result.minimum - C1(ρfin, ρfin) < tol

end

@testset "Robust X gate pulse inplace static arrays" begin

    prob = Problem(
        B = [SSx, SSy],
        A = SSz,
        Xi = UinitS,
        Xt = UfinS,
        T = 10.0,
        n_controls = 2,
        guess = rand(2, 100),
        sys_type = UnitaryGate(),
    )

    ens = EnsembleProblem(
        prob = prob,
        n_ens = 5,
        A_g = A_gens_static,
        B_g = B_gens_static,
        XiG = U_init_gens_static,
        XtG = U_target_gens_static,
        wts = ones(5) / 5,
    )

    sol = solve(
        ens,
        GRAPE(
            n_slices = 100,
            isinplace = false,
            optim_options = Optim.Options(f_tol = 1e-3),
        ),
    )
    @test sol.result.minimum - C1(UinitS, UfinS) < tol

end

@testset "Unitary X gate pulse ADGRAPE" begin

    prob = Problem(
        B = [Sx, Sy],
        A = 0.0 * Sz,
        Xi = Uinit,
        Xt = Ufin,
        T = 1.0,
        n_controls = 2,
        guess = rand(2, 25),
        sys_type = UnitaryGate(),
    )

    sol = solve(prob, ADGRAPE(n_slices = 25))

    pw_evolve(prob.A, prob.B, sol.opti_pulses, 2, 1.0 / 25, 25, Uinit)

    @test sol.result.minimum - C1(Ufin, Ufin) < tol
end

@testset "Robust X gate pulse ADGRAPE" begin

    prob = Problem(
        B = [Sx, Sy],
        A = Sz,
        Xi = Uinit,
        Xt = Ufin,
        T = 5.0,
        n_controls = 2,
        guess = rand(2, 100),
        sys_type = UnitaryGate(),
    )

    ens = EnsembleProblem(
        prob = prob,
        n_ens = 5,
        A_g = A_gens,
        B_g = B_gens,
        XiG = U_init_gens,
        XtG = U_target_gens,
        wts = ones(5) / 5,
    )

    sol = solve(ens, ADGRAPE(n_slices = 100))

    pw_evolve(prob.A, prob.B, sol.opti_pulses, 2, 5 / 100, 100, Uinit)

    @test sol.result.minimum - C1(ρfin, ρfin) < tol

end