text:
"""
The `Movie` and `Scene` types and the `animate()` function are designed to help you
create the frames that can be used to make an animated GIF or movie.

1 Provide width, height, title, and optionally a frame range to the Movie constructor:

    demo = Movie(400, 400, "test", 1:500)

2 Define one or more scenes and scene-drawing functions.

3 Run the `animate()` function, calling those scenes.

Example

```
bang = Movie(400, 100, "bang")

backdrop(scene, framenumber) =  background("black")

function frame1(scene, framenumber)
    background("white")
    sethue("black")
    eased_n = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    circle(O, 40 * eased_n, :fill)
end

animate(bang, [
    Scene(bang, backdrop, 0:200),
    Scene(bang, frame1, 0:200, easingfunction=easeinsine)],
    creategif=true,
    pathname="/tmp/animationtest.gif")

```
"""
struct Movie
    width::Float64
    height::Float64
    movietitle::AbstractString
    movieframerange::AbstractRange
end

"""
    Movie(width, height, movietitle)

Define a movie, specifying the width, height, and a title. The title will be used to
make the output file name. The range defaults to `1:250`.
"""
Movie(width, height, movietitle::AbstractString) = Movie(width, height, movietitle, 1:250)

"""default linear transition - no easing, no acceleration"""
function lineartween(t, b, c, d)
   return c * t/d + b
end

"""
The Scene type defines a function to be used to render a range of frames in a movie.

- the `movie` created by Movie()
- the `framefunction` is a function taking two arguments: the scene and the framenumber.
- the `framerange` determines which frames are processed by the function. Defaults to the entire movie.
- the optional `easingfunction` can be accessed by the framefunction to vary the transition speed
- the optional `opts` which is a single argument of an abstract type which can be accessed within the framefunction
"""
mutable struct Scene
    movie::Movie
    framefunction::Function
    framerange::AbstractRange
    easingfunction::Function
    opts::Any
end

"""
    Scene(movie, function, range;
        easingfunction=easinoutquad,
        optarg=nothing)

Use the Scene() constructor function to create a scene. Supply a movie, a function
to generate the scene, and a range of frames. Optionally you can supply an
easing function, and other information, in `optarg`, which can be accessed as
`scene.opts`.

### Example

```
function initial(scene, framenumber)
    balls = scene.opts
    ...
end

animate(poolmovie, [
    Scene(poolmovie, initial, optarg=balls,   1:20),
    ...
    ])
```

To use an easing function inside the frame-generating function, you can create a normalized
value with, for example:

```
eased_n = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
```

Or, if the scene doesn't start at frame 1, calculate normalized easing function like this:

```
eased_n = scene.easingfunction(framenumber - scene.framerange.start,
    0, 1, scene.framerange.stop - scene.framerange.start)
```
"""
Scene(movie::Movie, framefunction::Function;
    easingfunction=lineartween, optarg=nothing) =
    Scene(movie, framefunction, movie.movieframerange, easingfunction, optarg)
Scene(movie::Movie, framefunction::Function, framerange::AbstractRange;
    easingfunction=lineartween, optarg=nothing) =
    Scene(movie, framefunction, framerange, easingfunction, optarg)

"""
    animate(movie::Movie, scenelist::Array{Scene, 1};
        creategif=false,
        framerate=30,
        pathname="",
        tempdirectory="",
        usenewffmpeg=true)

Create the movie defined in `movie` by rendering the frames define in the array of scenes
in `scenelist`.

If `creategif` is `true`, the function attempts to call the `ffmpeg` utility on
the resulting frames to build a GIF animation. This will be stored in `pathname`
(an existing file will be overwritten; use a ".gif" suffix), or in
`(movietitle).gif` in a temporary directory. `ffmpeg` should be installed and
available, of course, if this is to work.

### Example

```
animate(bang, [
    Scene(bang, backdrop, 0:200),
    Scene(bang, frame1, 0:200, easingfunction=easeinsine)],
    creategif=true,
    pathname="/tmp/animationtest.gif")
```

The `usenewffmpeg` option, `true` by default, uses single-pass palette
generation and more complex filtering provided by recent versions of the
`ffmpeg` utility, mainly to cope with transparent backgrounds. If set to
`false`, the behavior is the same as in previous versions of Luxor.
"""
function animate(movie::Movie, scenelist::Array{Scene, 1};
        creategif=false,
        framerate=30,
        pathname="",
        tempdirectory="",
        usenewffmpeg=true)

    if tempdirectory == ""
        tempdirectory = mktempdir()
    else
        if !isdir(tempdirectory)
            @info "$(tempdirectory) does not exist, using an autogenerated one."
            tempdirectory = mktempdir()
        end
    end
    @info("Frames for animation \"$(movie.movietitle)\" are being stored in directory: \n\t $(tempdirectory)")
    filecounter = 1
    rangelist = 0:-1
    for scene in scenelist
        rangelist = vcat(rangelist, collect(scene.framerange))
    end
    rangelist = unique(rangelist) # remove shared frames
    if rangelist[end] < movie.movieframerange.stop
        @warn("Movie framerange is longer than scene frame range: \n\t $(movie.movieframerange) > $(rangelist[end])")
    end
    for currentframe in rangelist[1]:rangelist[end]
        Drawing(movie.width, movie.height, "$(tempdirectory)/$(lpad(filecounter, 10, "0")).png")
        origin()
        # this frame needs doing, see if each of the scenes defines it
        for scene in scenelist
            if currentframe in scene.framerange
                scene.framefunction(scene, currentframe)
            end
        end
        finish()
        filecounter += 1
    end
    @info("... $(filecounter-1) frames saved in directory:\n\t $(tempdirectory)")
    if creategif == true
        # these two commands create a palette and then create animated GIF from the resulting images
        if !usenewffmpeg
            # old version of ffmpeg up to 2.1.3
            # these two commands create a palette and then an animated GIF from the resulting images using the palette
            run(`ffmpeg -loglevel panic -f image2 -i $(tempdirectory)/%10d.png -vf palettegen -y $(tempdirectory)/$(movie.movietitle)-palette.png`)
            run(`ffmpeg -loglevel panic -framerate $(framerate) -f image2 -i $(tempdirectory)/%10d.png -i $(tempdirectory)/$(movie.movietitle)-palette.png -lavfi paletteuse -y $(tempdirectory)/$(movie.movietitle).gif`)
        else
            # the latest version of ffmpeg uses built-in palettes and allegedly does transparency using complex filters
            run(`ffmpeg -loglevel panic -framerate $(framerate) -f image2 -i $(tempdirectory)/%10d.png -filter_complex "[0:v] split [a][b]; [a] palettegen=stats_mode=full:reserve_transparent=on:transparency_color=FFFFFF [p]; [b][p] paletteuse=new=1:alpha_threshold=128" -y $(tempdirectory)/$(movie.movietitle).gif`)
        end

        if ! isempty(pathname)
            mv("$(tempdirectory)/$(movie.movietitle).gif", pathname, force=true)
            @info("GIF is: $pathname")
        else
            @info("GIF is: $(tempdirectory)/$(movie.movietitle).gif")
        end
    end
    return true
end

"""
    animate(movie::Movie, scene::Scene; creategif=false, framerate=30)

Create the movie defined in `movie` by rendering the frames define in `scene`.
"""
animate(movie::Movie, scene::Scene; kwargs...) = animate(movie, [scene]; kwargs...)

"""
    easingflat(t, b, c, d)

A flat easing function, same as `lineartween()`.

For all easing functions, the four parameters are:

- `t` time, ie the current framenumber
- `b` beginning position or bottom value of the range
- `c` total change in position or top value of the range
- `d` duration, ie a framecount

1. `t/d` or `t/=d` normalizes `t` to between 0 and 1
2. `... * c` scales up to the required range value
3. `... + b` adds the initial offset

"""
function easingflat(t, b, c, d)
    return c * t / d + b
end

"""
    easeinquad(t, b, c, d)

quadratic easing in - accelerating from zero velocity
"""
function easeinquad(t, b, c, d)
   t /= d
   return c * t * t + b
end

"""
    easeoutquad(t, b, c, d)

quadratic easing out - decelerating to zero velocity
"""
function easeoutquad(t, b, c, d)
   t /= d
   return -c * t * (t - 2) + b
end

"""
    easeinoutquad(t, b, c, d)

quadratic easing in/out - acceleration until halfway, then deceleration
"""
function easeinoutquad(t, b, c, d)
   t /= d/2
   if t < 1
      return  (c/2) * t * t + b
    end
   t -= 1
   return -(c/2) * (t * (t - 2) - 1) + b
end

"""
    easeincubic(t, b, c, d)

cubic easing in - accelerating from zero velocity
"""
function easeincubic(t, b, c, d)
   t /= d
   return c * t * t * t + b
end

"""
    easeoutcubic(t, b, c, d)

cubic easing out - decelerating to zero velocity
"""
function easeoutcubic(t, b, c, d)
   t /= d
   t -= 1
   return c * (t * t * t + 1) + b
end

"""
    easeinoutcubic(t, b, c, d)

cubic easing in/out - acceleration until halfway, then deceleration
"""
function easeinoutcubic(t, b, c, d)
   t /= d/2
   if t < 1
       return c/2 * t * t * t + b
   end
   t -= 2
   return c/2 * (t * t * t + 2) + b
end

"""
    easeinquart(t, b, c, d)

quartic easing in - accelerating from zero velocity
"""
function easeinquart(t, b, c, d)
   t /= d
   return c * t * t * t * t + b
end

"""
    easeoutquart(t, b, c, d)

quartic easing out - decelerating to zero velocity
"""
function easeoutquart(t, b, c, d)
   t /= d
   t -= 1
   return -c * (t * t * t * t - 1) + b
end

"""
    easeinoutquart(t, b, c, d)

quartic easing in/out - acceleration until halfway, then deceleration
"""
function easeinoutquart(t, b, c, d)
   t /= d/2
   if t < 1
       return c/2 * t * t * t * t + b
   end
   t -= 2
   return -c/2 * (t * t * t * t - 2) + b
end

"""
    easeinquint(t, b, c, d)

quintic easing in - accelerating from zero velocity
"""
function easeinquint(t, b, c, d)
   t /= d
   return c * t * t * t * t * t + b
end

"""
    easeoutquint(t, b, c, d)

quintic easing out - decelerating to zero velocity
"""
function easeoutquint(t, b, c, d)
   t /= d
   t -= 1
   return c * (t * t * t * t * t + 1) + b
end

"""
    easeinoutquint(t, b, c, d)

quintic easing in/out - acceleration until halfway, then deceleration
"""
function easeinoutquint(t, b, c, d)
   t /= d/2
   if t < 1
       return c/2 * t * t * t * t * t + b
   end
   t -= 2
   return c/2 * (t * t * t * t * t + 2) + b
end

"""
    easeinsine(t, b, c, d)

sinusoidal easing in - accelerating from zero velocity
"""
function easeinsine(t, b, c, d)
   return -c * cos(t/d * (pi/2)) + c + b
end

"""
    easeoutsine(t, b, c, d)

sinusoidal easing out - decelerating to zero velocity
"""
function easeoutsine(t, b, c, d)
   return c * sin(t/d * (pi/2)) + b
end

"""
    easeinoutsine(t, b, c, d)

sinusoidal easing in/out - accelerating until halfway, then decelerating
"""
function easeinoutsine(t, b, c, d)
   return -c/2 * (cos(pi * t/d) - 1) + b
end

"""
    easeinexpo(t, b, c, d)

exponential easing in - accelerating from zero velocity
"""
function easeinexpo(t, b, c, d)
    if isapprox(t, 0.0, atol=1e-3)
        return 0.0
    else
        return c * 2 ^ (10 * (t/d - 1)) + b
    end
end

"""
    easeoutexpo(t, b, c, d)

exponential easing out - decelerating to zero velocity
"""
function easeoutexpo(t, b, c, d)
   return c * (-(2 ^ (-10 * t/d)) + 1) + b
end

"""
    easeinoutexpo(t, b, c, d)

exponential easing in/out - accelerating until halfway, then decelerating
"""
function easeinoutexpo(t, b, c, d)
   t /= d/2
   if t < 1
       return c/2 * 2 ^ (10 * (t - 1)) + b
   end
   t -= 1
   return c/2 * (-(2 ^ (-10 * t)) + 2) + b
end

"""
    easeincirc(t, b, c, d)

circular easing in - accelerating from zero velocity
"""
function easeincirc(t, b, c, d)
   t /= d
   return -c * (sqrt(abs((1 - t * t))) - 1) + b
end

"""
    easeoutcirc(t, b, c, d)

circular easing out - decelerating to zero velocity
"""
function easeoutcirc(t, b, c, d)
   t /= d
   t -= 1
   return c * sqrt(abs((1 - t * t))) + b
end

"""
    easeinoutcirc(t, b, c, d)

circular easing in/out - acceleration until halfway, then deceleration
"""
function easeinoutcirc(t, b, c, d)
   t /= d/2
   if (t < 1)
       return -c/2 * ((sqrt(abs(1 - t * t))) - 1) + b
   end
   t -= 2
   return c/2 * (sqrt(abs(1 - t * t)) + 1) + b
end

"""
    easeinoutinversequad(t, b, c, d)

ease in, then slow down, then speed up, and ease out
"""
function easeinoutinversequad(t, b, c, d)
   t /= d/2
   if t <= 1
      return -(c/2) * (t * (t - 2)) + b
    end
   t -= 1
   return c/2 * (t * t * t + 1) + b
end

"""
    easeinoutbezier(t, b, c, d, cpt1, cpt2)

This easing function takes six arguments, the usual `t`, `b`, `c`, and `d`, but
also two points. These are the normalized control points of a Bezier curve drawn
between `Point(0, 0)` to `Point(1.0, 1.0)`. The `y` value of the Bezier is the
eased value for `t`.

In your `frame()` generating function, if a Scene specifies the `easeinoutbezier` easing function, you can use this:

```
...
lineareasing = rescale(framenumber, 1, scene.framerange.stop)
beziereasing = scene.easingfunction(lineareasing, 0, 1, 1,
    Point(0.25, 0.25), Point(0.75, 0.75))
...
```

These two control points lie on the line between `0/0` and `1/1`, so it's equivalent to a linear easing (`lineartween()` or `easingflat`).

However, in the next example, the two control points define a wave-like curve
that changes direction before changing back. When animating with this easing
function, an object will 'go retrograde' for a while.

```
lineareasing = rescale(framenumber, 1, scene.framerange.stop)
beziereasing = scene.easingfunction(lineareasing, 0, 1, 1,
    Point(0.01, 1.99), Point(0.99, -1.5))
```

"""
function easeinoutbezier(t, b, c, d, cpt1::Point=Point(0.25, 0.25), cpt2::Point=Point(0.75, 0.75))
    t /= d + b
    bez = bezier(t, Point(0, 0), cpt1, cpt2, Point(1.0, 1.0))
    return bez.y * c
end

easingfunctions = [lineartween,
                    easeinquad,
                    easeoutquad,
                    easeinoutquad,
                    easeincubic,
                    easeoutcubic,
                    easeinoutcubic,
                    easeinquart,
                    easeoutquart,
                    easeinoutquart,
                    easeinquint,
                    easeoutquint,
                    easeinoutquint,
                    easeinsine,
                    easeoutsine,
                    easeinoutsine,
                    easeinexpo,
                    easeoutexpo,
                    easeinoutexpo,
                    easeincirc,
                    easeoutcirc,
                    easeinoutcirc,
                    easingflat,
                    easeinoutinversequad,
                    easeinoutbezier]


text:
function default_samplemapper(rig_name::AbstractString, daq_chan_name::String; sample_rate = 10000s^-1)
    #is_digi_funcs = [iscam, islas, isstim, iscammonitor]
    #any(map(f->f(daq_chan_name, rig_name)), is_digi_funcs)
    if isdigital(daq_chan_name, rig_name)
        return ttl_samplemapper(; sample_rate = sample_rate)
    elseif ispos(daq_chan_name, rig_name) || isposmonitor(daq_chan_name, rig_name)
        return piezo_samplemapper(PIEZO_RANGES[rig_name]...; rawtype = Int16, sample_rate = sample_rate)
    elseif isanalog(daq_chan_name, rig_name)
        if isoutput(daq_chan_name, rig_name)
            return generic_ao_samplemapper(-10.0V..10.0V; rawtype = Int16, sample_rate = sample_rate)
        else
            return generic_ai_samplemapper(-10.0V..10.0V; rawtype = Int16, sample_rate = sample_rate)
        end
    else
        error("Unrecognized channel name")
    end
end

function generic_ao_samplemapper(v::AbstractInterval{TV};
                                 rawtype=Int16, sample_rate::HasInverseTimeUnits{Int, TU}=10000s^-1) where{TV<:HasVoltageUnits, TU}
    return SampleMapper(typemin(rawtype), typemax(rawtype), minimum(v), maximum(v), minimum(v), maximum(v), sample_rate)
end

generic_ai_samplemapper = generic_ao_samplemapper

function piezo_samplemapper(p::AbstractInterval{TL}, v::AbstractInterval{TV};
                            rawtype=Int16, sample_rate::HasInverseTimeUnits{Int, TU}=10000s^-1) where{TL<:HasLengthUnits,TV<:HasVoltageUnits, TU}
    return SampleMapper(zero(rawtype), typemax(rawtype), minimum(v), maximum(v), minimum(p), maximum(p), sample_rate)
end

function galvo_ctrl_samplemapper(rawtype=Int16, sample_rate::HasInverseTimeUnits{Int, TU}=10000s^-1) where TU
    rad_min = -10.0*deg2rad(1.0)*Unitful.rad #1.0 degrees per volt
    rad_max = 10.0*deg2rad(1.0)*Unitful.rad
    return SampleMapper(zero(rawtype), typemax(rawtype), -10.0*Unitful.V, 10.0*Unitful.V, rad_min, rad_max, sample_rate)
end

function galvo_mon_samplemapper(rawtype=Int16, sample_rate::HasInverseTimeUnits{Int, TU}=10000s^-1) where TU
    rad_min = -5.0*deg2rad(2.0)*Unitful.rad #2.0 degrees per volt
    rad_max = 5.0*deg2rad(2.0)*Unitful.rad
    return SampleMapper(zero(rawtype), typemax(rawtype), -5.0*Unitful.V, 5.0*Unitful.V, rad_min, rad_max, sample_rate)
end

#Shortcut for creating a generic TTL SampleMapper, assumes TTL level of 3.3V (though this doesn't matter to Imagine, only for visualizing in Julia)
#By default the raw samples are true/false values encoded as UInt8.  If using an analog channel for TTL signals these default limits should be changed
function ttl_samplemapper(rawmin=UInt8(false), rawmax=(UInt8(true)); sample_rate::HasInverseTimeUnits{Int, U}=10000s^-1) where U
    return SampleMapper(rawmin, rawmax, 0.0*Unitful.V, 3.3*Unitful.V, false, true, sample_rate)
end

#returns an array of empty ImagineSignals, one for each channel accessible to OCPI2 users
function rigtemplate(rig::AbstractString; sample_rate::HasInverseTimeUnits{Int,U} = 10000s^-1) where U
    if !in(rig, RIGS)
        error("Unsupported rig")
    end
    coms = ImagineSignal[]
    shared_dict = Dict()
    name_lookup = DEFAULT_DAQCHANS_TO_NAMES[rig]
    #analog outputs
    ao_sampmapper = 0
    for c in AO_CHANS[rig]
        if ispos(c, rig)
            ao_sampmapper = piezo_samplemapper(PIEZO_RANGES[rig]...; rawtype = Int16, sample_rate = sample_rate)
        elseif iscam(c, rig) #if using an AO channel to control cameras (not advised, mostly for testing)
            ao_sampmapper = ttl_samplemapper(zero(Int16), ceil(Int16, typemax(Int16)*3.3/10.0); sample_rate = sample_rate)
        elseif isgalvo(c, rig)
            ao_sampmapper = galvo_ctrl_samplemapper(Int16, sample_rate)
        else
            ao_sampmapper = generic_ao_samplemapper(AO_RANGE[rig]; rawtype = Int16, sample_rate = sample_rate)
        end
        ao_vectype = RLEVector{rawtype(ao_sampmapper)}
        push!(coms, ImagineSignal{ao_vectype}(name_lookup[c], c, rig, [], String[], shared_dict, Int[], ao_sampmapper))
    end

    #digital outputs (includes cameras, lasers, and stimulus channels)
    do_sampmapper = ttl_samplemapper(;sample_rate = sample_rate)
    do_vectype = RLEVector{rawtype(do_sampmapper)}
    for c in DO_CHANS[rig]
        push!(coms, ImagineSignal{do_vectype}(name_lookup[c], c, rig, [], String[], shared_dict, Int[], do_sampmapper))
    end

    #analog inputs
    ai_sampmapper = 0
    for c in AI_CHANS[rig]
        if isposmonitor(c, rig)
            ai_sampmapper = piezo_samplemapper(PIEZO_RANGES[rig]...; rawtype = Int16, sample_rate = sample_rate)
        elseif iscammonitor(c, rig) #If using an AI channel for TTL camera inputs
            ai_sampmapper = ttl_samplemapper(zero(Int16), ceil(Int16, typemax(Int16)*3.3/10.0); sample_rate = sample_rate)
        elseif isgalvomonitor(c, rig)
            ai_sampmapper = galvo_mon_samplemapper(Int16, sample_rate)
        else
            ai_sampmapper = generic_ai_samplemapper(AI_RANGE[rig]; rawtype = Int16, sample_rate = sample_rate)
        end
        ai_vectype = Vector{rawtype(ai_sampmapper)}
        push!(coms, ImagineSignal{ai_vectype}(name_lookup[c], c, rig, [], String[], shared_dict, Int[], ai_sampmapper))
    end

    #digital inputs (including cameras)
    di_sampmapper = do_sampmapper
    di_vectype = Vector{rawtype(di_sampmapper)}
    for c in DI_CHANS[rig]
        push!(coms, ImagineSignal{di_vectype}(name_lookup[c], c, rig, [], String[], shared_dict, Int[], di_sampmapper)) #TODO: handle bit-packing (.di file convention)
    end

    return coms
end

ext: jl

max_forks_repo_name: HolyLab/ImagineInterface


text:
module ROMEO

using Statistics
using StatsBase

const NBINS = 256

include("utility.jl")
include("priorityqueue.jl")
include("weights.jl")
include("seed.jl")
include("region_handling.jl")
include("algorithm.jl")
include("unwrapping.jl")

export unwrap, unwrap!, unwrap_individual, unwrap_individual!

end # module


text:
#*****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module NNGrammarExpt1

export circuit_andor

using TFTools
using Datasets
using TensorFlow
import TensorFlow: DT_FLOAT32
import TensorFlow.API: l2_loss, AdamOptimizer, cast, round_, maximum_, minimum_

function circuit_andor(;
    datname::AbstractString="bin_synth",
    featfile::AbstractString="feats",
    labelfile::AbstractString="labels",
    labelfield::AbstractString="x1_and_x3",
    learning_rate::Float64=0.005,
    training_epochs::Int64=100,
    batch_size::Int64=1000,
    hidden_units::Vector{Int64}=Int64[15, 10, 15],
    display_step::Int64=1,
    b_debug::Bool=false)

    Dfeats = dataset(datname, featfile)
    Dlabels = dataset(datname, labelfile)
    n_feats = ncol(Dfeats)
    @assert nrow(Dfeats) == nrow(Dlabels)

    data_set = TFDataset(Dfeats, Dlabels[Symbol(labelfield)])

    # Construct model

    # inputs
    muxselect = Placeholder(DT_FLOAT32, [-1, n_feats])
    inputs = Placeholder(DT_FLOAT32, [-1, n_feats])

    # x mux
    x_muxin = inputs
    x_mux = SoftMux(n_feats, n_feats, hidden_units, Tensor(x_muxin), Tensor(muxselect))
    x_muxout = out(x_mux)

    # y mux
    y_muxin = inputs
    y_mux = SoftMux(n_feats, n_feats, hidden_units, Tensor(y_muxin), Tensor(muxselect))
    y_muxout = out(y_mux)

    # op block
    ops_in = (x_muxout, y_muxout)
    ops_list = [ops_And, ops_Or]
    ops_blk = OpsBlock(ops_in, ops_list)
    ops_out = out(ops_blk)

    # op mux
    op_muxin = ops_out
    op_mux = SoftMux(num_ops(ops_blk), n_feats, hidden_units, Tensor(op_muxin), Tensor(muxselect))
    op_muxout = out(op_mux)

    # outputs
    pred = op_muxout
    labels = Placeholder(DT_FLOAT32, [-1])

    # Define loss and optimizer
    cost = l2_loss(pred - labels) # Squared loss
    optimizer = minimize(AdamOptimizer(learning_rate), cost) # Adam Optimizer

    # Initializing the variables
    init = initialize_all_variables()

    # Rock and roll
    sess = Session()
    try
        run(sess, init)

        #debug
        #fd = FeedDict(inputs => data_set.X, muxselect => data_set.X, labels => data_set.Y)
        #@bp
        #tmp=1
        #run(sess, x_muxout, fd)
        #run(sess, y_muxout, fd)
        #run(sess, ops_out, fd)
        #run(sess, op_muxout, fd)
        #run(sess, cost, fd)
        #/debug

        # Training cycle
        for epoch in 1:training_epochs
            avg_cost = 0.0
            total_batch = div(num_examples(data_set), batch_size)

            # Loop over all batches
            for i in 1:total_batch
                batch_xs, batch_ys = next_batch(data_set, batch_size)
                fd = FeedDict(inputs => batch_xs, muxselect => batch_xs, labels => batch_ys)
                # Fit training using batch data
                run(sess, optimizer, fd)
                # Compute average loss
                batch_average_cost = run(sess, cost, fd)
                avg_cost += batch_average_cost / (total_batch * batch_size)
            end

            # Display logs per epoch step
            if epoch % display_step == 0
                println("Epoch $(epoch)  cost=$(avg_cost)")
            end
        end
        println("Optimization Finished")

        # Test model
        correct_prediction = (round_(pred) == labels)
        # Calculate accuracy
        accuracy = mean(cast(correct_prediction, DT_FLOAT32))
        fd = FeedDict(inputs => data_set.X, muxselect => data_set.X, labels => data_set.Y)
        acc = run(sess, accuracy, fd)
        println("Accuracy:", acc)

        if b_debug
            #reload data_set to recover original order
            data_set = TFDataset(Dfeats, Dlabels[Symbol(labelfield)])
            db_x = data_set.X
            db_labels = data_set.Y
            db_xmux = run(sess, x_muxout, fd)
            db_xmux_nnout = run(sess, x_mux.nnout, fd)
            db_xmux_hardselect = run(sess, x_mux.hardselect, fd)
            db_ymux = run(sess, y_muxout, fd)
            db_ymux_nnout = run(sess, y_mux.nnout, fd)
            db_ymux_hardselect = run(sess, y_mux.hardselect, fd)
            db_opmux = run(sess, op_muxout, fd)
            db_opmux_nnout = run(sess, op_mux.nnout, fd)
            db_opmux_hardselect = run(sess, op_mux.hardselect, fd)
            db_ypred = run(sess, pred, fd)
            #combine these into a single call to prevent multiple forward passes
            #run(sess, [x_muxout, y_muxout, ops_out, op_muxout, cost], fd)

            NSHOW = 20
            @show db_x[1:NSHOW]
            @show db_xmux[1:NSHOW]
            @show db_xmux_nnout[1:NSHOW,1:7]
            @show db_ymux[1:NSHOW]
            @show db_ymux_nnout[1:NSHOW,1:7]
            @show db_ymux[1:NSHOW]
            @show db_opmux[1:NSHOW]
            @show db_opmux_nnout[1:NSHOW,1:2]
            @show db_ypred[1:NSHOW]
            @show db_labels[1:NSHOW]
            @show db_xmux_hardselect[1:NSHOW]
            @show db_ymux_hardselect[1:NSHOW]
            @show db_opmux_hardselect[1:NSHOW]
            println("Accuracy:", acc)
        end
    finally
        close(sess)
    end
end

function ops_And(x::Tensor, y::Tensor)
    minimum_(x, y) #alternatively, use tf.logical_and?
end

function ops_Or(x::Tensor, y::Tensor)
    maximum_(x, y) #alternatively, use tf.logical_or?
end

end #module


text:
function loopheader!(integ)
    integ.iter += 1
    fix_dt_at_bounds!(integ)
end

function fix_dt_at_bounds!(integ)
    integ.dt = min(integ.opts.dtmax, integ.dt)
    integ.opts.force_dtmin && (integ.dt = max(integ.opts.dtmin, integ.dt))

    if integ.dt^(1/integ.constants.q) < eps(eltype(integ.u))
        @warn "Small step size: h^q < eps(u)! Continuing, but this could lead to numerical issues" integ.dt
    end

    next_t = integ.t + integ.dt
    if next_t + integ.opts.dtmin^(1/2) > integ.prob.tspan[2]
        # Avoid having to make a step smaller than dtmin in the next step
        @debug "Increasing the step size now to avoid having to do a really small step next to hit t_max"
        integ.dt = integ.prob.tspan[2] - integ.t
    end
end

function loopfooter!(integ)

    integ.accept_step = integ.opts.adaptive ? integ.EEst < 1 : true

    integ.opts.adaptive && (integ.dt = propose_step!(integ.steprule, integ))

    integ.accept_step ? (integ.destats.naccept += 1) : (integ.destats.nreject += 1)

    any(isnan.(integ.cache.σ_sq)) && error("Estimated sigma is NaN")
    integ.opts.adaptive && isnan(integ.EEst) && error("Error estimate is NaN")

    # Accept or reject the step: Moved this here from loopheader!
    if integ.iter > 0 && integ.accept_step
        integ.success_iter += 1
        apply_step!(integ)
    end
end

function postamble!(integ)
    if isstatic(integ.sigma_estimator)
        calibrate!(integ)
        integ.sigmas = [integ.sigmas[end] for s in integ.sigmas]
    end
    integ.smooth && smooth_all!(integ)

    integ.sol = DiffEqBase.build_solution(
        integ.prob, integ.alg,
        integ.times,
        integ.state_estimates,
        integ.sigmas,
        integ.proposals, integ;
        retcode=integ.retcode,
        destats=integ.destats)
end

function apply_step!(integ)
    copy!(integ.cache.x, integ.cache.x_filt)
    copy!(integ.u, integ.cache.u_filt)
    integ.t = integ.t_new

    integ.cache.σ_sq_prev = integ.cache.σ_sq

    # For the solution
    push!(integ.state_estimates, copy(integ.cache.x))
    push!(integ.times, copy(integ.t))
    push!(integ.sigmas, copy(integ.cache.σ_sq))
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# ---
# title: 1376. Time Needed to Inform All Employees
# id: problem1376
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search
# link: <https://leetcode.com/problems/time-needed-to-inform-all-employees/description/>
# hidden: true
# ---
#
# A company has `n` employees with a unique ID for each employee from `0` to `n
# - 1`. The head of the company has is the one with `headID`.
#
# Each employee has one direct manager given in the `manager` array where
# `manager[i]` is the direct manager of the `i-th` employee, `manager[headID] =
# -1`. Also it's guaranteed that the subordination relationships have a tree
# structure.
#
# The head of the company wants to inform all the employees of the company of an
# urgent piece of news. He will inform his direct subordinates and they will
# inform their subordinates and so on until all employees know about the urgent
# news.
#
# The `i-th` employee needs `informTime[i]` minutes to inform all of his direct
# subordinates (i.e After informTime[i] minutes, all his direct subordinates can
# start spreading the news).
#
# Return _the number of minutes_ needed to inform all the employees about the
# urgent news.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1, headID = 0, manager = [-1], informTime = [0]
#     Output: 0
#     Explanation: The head of the company is the only employee in the company.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/27/graph.png)
#
#
#
#     Input: n = 6, headID = 2, manager = [2,2,-1,2,2,2], informTime = [0,0,1,0,0,0]
#     Output: 1
#     Explanation: The head of the company with id = 2 is the direct manager of all the employees in the company and needs 1 minute to inform them all.
#     The tree structure of the employees in the company is shown.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/28/1730_example_3_5.PNG)
#
#
#
#     Input: n = 7, headID = 6, manager = [1,2,3,4,5,6,-1], informTime = [0,6,5,4,3,2,1]
#     Output: 21
#     Explanation: The head has id = 6. He will inform employee with id = 5 in 1 minute.
#     The employee with id = 5 will inform the employee with id = 4 in 2 minutes.
#     The employee with id = 4 will inform the employee with id = 3 in 3 minutes.
#     The employee with id = 3 will inform the employee with id = 2 in 4 minutes.
#     The employee with id = 2 will inform the employee with id = 1 in 5 minutes.
#     The employee with id = 1 will inform the employee with id = 0 in 6 minutes.
#     Needed time = 1 + 2 + 3 + 4 + 5 + 6 = 21.
#
#
# **Example 4:**
#
#
#
#     Input: n = 15, headID = 0, manager = [-1,0,0,1,1,2,2,3,3,4,4,5,5,6,6], informTime = [1,1,1,1,1,1,1,0,0,0,0,0,0,0,0]
#     Output: 3
#     Explanation: The first minute the head will inform employees 1 and 2.
#     The second minute they will inform employees 3, 4, 5 and 6.
#     The third minute they will inform the rest of employees.
#
#
# **Example 5:**
#
#
#
#     Input: n = 4, headID = 2, manager = [3,3,-1,2], informTime = [0,0,162,914]
#     Output: 1076
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#   * `0 <= headID < n`
#   * `manager.length == n`
#   * `0 <= manager[i] < n`
#   * `manager[headID] == -1`
#   * `informTime.length == n`
#   * `0 <= informTime[i] <= 1000`
#   * `informTime[i] == 0` if employee `i` has no subordinates.
#   * It is **guaranteed** that all the employees can be informed.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@apisetup begin
    PyObject_GetAttrString = PySym(:PyObject_GetAttrString)
    PyCFunction_NewEx = PySym(:PyCFunction_NewEx)

    # object
    PyObject_Str = PySym(:PyObject_Str)
    PyObject_VectorcallDict = PySymMaybe(:PyObject_VectorcallDict)
    PyObject_CallFunctionObjArgs = PySym(:PyObject_CallFunctionObjArgs)
    PyObject_CallMethodObjArgs = PySym(:PyObject_CallMethodObjArgs)
    PyObject_GetItem = PySym(:PyObject_GetItem)
    PyObject_SetItem = PySym(:PyObject_SetItem)
    PyObject_RichCompare = PySym(:PyObject_RichCompare)
    PyObject_IsInstance = PySym(:PyObject_IsInstance)
    PyObject_SetAttr = PySym(:PyObject_SetAttr)
    PyObject_GetAttr = PySym(:PyObject_GetAttr)

    # sequence
    PySequence_Length = PySym(:PySequence_Length)
    PySequence_GetItem = PySym(:PySequence_GetItem)

    # number
    PyNumber_AsSsize_t = PySym(:PyNumber_AsSsize_t)
    PyNumber_Long = PySym(:PyNumber_Long)

    # unicode
    PyUnicode_FromString = PySym(:PyUnicode_FromString)

    # int/long.
    # Int-prefixed means it's not an exported symbol from libpython
    # TODO: refactor names.
    Py_IntAsNumber = unsafe_load(reinterpret(Ptr{PyTypeObject}, PyO.int)).tp_as_number
    PyInt_Compare = unsafe_load(reinterpret(Ptr{PyTypeObject}, PyO.int)).tp_richcompare
    Py_IntPow = unsafe_load(Py_IntAsNumber).nb_power
    Py_IntAddInt = unsafe_load(Py_IntAsNumber).nb_add
    PyLong_AsDouble = PySym(:PyLong_AsDouble)
    PyLong_AsSsize_t = PySym(:PyLong_AsSsize_t)
    PyLong_FromSsize_t = PySym(:PyLong_FromSsize_t)
    PyLong_FromLong = PySym(:PyLong_FromLong)
    PyLong_FromSize_t = PySym(:PyLong_FromSize_t)

    # float
    PyFloat_FromDouble = PySym(:PyFloat_FromDouble)

    # bytes
    _PyBytes_Join = PySym(:_PyBytes_Join)

    # dict
    _PyDict_GetItem_KnownHash = PySym(:_PyDict_GetItem_KnownHash)
    PyDict_GetItemWithError = PySym(:PyDict_GetItemWithError)

    # list
    PyList_Append = PySym(:PyList_Append)
    PyList_GetItem = PySym(:PyList_GetItem)
    PyList_SetItem = PySym(:PyList_SetItem)
    PyList_New = PySym(:PyList_New)

    # function
    PyFunction_GetGlobals = PySym(:PyFunction_GetGlobals)

    # module
    PyModule_GetDict = PySym(:PyModule_GetDict)

    # exc
    PyExc_KeyError = PySym(:PyExc_KeyError)
    PyExc_TypeError = PySym(:PyExc_TypeError)

    # err
    PyErr_Clear = PySym(:PyErr_Clear)
    PyErr_SetString = PySym(:PyErr_SetString)
    PyErr_Print = PySym(:PyErr_Print)
    PyErr_ExceptionMatches = PySym(:PyErr_ExceptionMatches)
    PyErr_Occurred = PySym(:PyErr_Occurred)
    PyErr_SetObject = PySym(:PyErr_SetObject)
    PyErr_SetNone = PySym(:PyErr_SetNone)
end


text:
######################
# This is SemiLinearHeat but using FFTW

using FFTW
x=range(0,stop=2π,length=100)
u(x) = sin(x)
freqs = fft(u.(x))[ 1 : length(x) ÷ 2 + 1]  ## ÷ this means integer division
c= 2*abs.(freqs/length(x))

######################
# instead of the above section
using ApproxFun
s = Fourier()
n = 100
T = ApproxFun.plan_transform(s,n)
Ti = ApproxFun.plan_itransform(s,n)

x = points(s,n)    # 0  to 2*pi
r = (T*(cos.(cos.(x.-0.1))))

der = Derivative(s,2)

####################
# PDE-derived ODE

using Sundials
using Plots; gr()
using ApproxFun
using Sundials
using LinearAlgebra
using LinearMaps

import Base.At_mul_B!

S = Fourier()
n=100
x = points(S,n)
T = ApproxFun.plan_transform(S,n)
Ti = ApproxFun.plan_itransform(S,n)

u0= T*cos.(cos.(x.-0.1))

D2 = Derivative(S,2)
L = D2[1:n,1:n]
tspan = (0.0,10.0)

function heat(du,u,L,t)
    mul!(du,L,u)
end

prob = ODEProblem(heat, u0,tspan,L)
sol = solve(prob)
plot(x,Ti*sol(0.0))
plot!(x,Ti*sol(0.5))
plot!(x,Ti*sol(2.0))
plot!(x,Ti*sol(3.0))
plot!(x,Ti*sol(10.0))

max_stars_repo_path: Missing around and learning DifferentialEquations.jl/demos.jl


text:
# db.jl

using SQLite
using DataFrames

# Database item
struct DatabaseItem{T}
   time::Float64
   agent::String
   action::String
   key::String
   item::T
end # DatabaseItem

# createDatabaseItem - internal function
const agent = "AB9F"
createDatabaseItem(item::Any; agent=agent, action="CREATE") = DatabaseItem(time(), agent, action, item._id, item)

# Connect with PATH_CSV
connect(path::String)::SQLite.DB = SQLite.DBInterface.connect(SQLite.DB, path)

# connect to in-memory database
connect()::SQLite.DB = SQLite.DBInterface.connect(SQLite.DB)

# close connection
disconnect(db::SQLite.DB) = begin
   SQLite.DBInterface.close!(db)
end # disconnect

# archive an item as DatabaseItem
archive(db, table::String, items::Array{T, 1} where {T <: Any}) = begin
   # save items
   dbi = [ createDatabaseItem( i ; action="CREATE" ) for i in items]
   # return as dataframe
   DataFrame( dbi ) |> SQLite.load!(db, table)
end # archive

#store the fields of an item
store(db, table::String, items::Array{T, 1} where {T <: Any}) = begin
   DataFrame( items ) |> SQLite.load!(db, table)
end # store

# retrieve all item from a table
retrieve(db, table::String)::DataFrame = begin
   r = SQLite.DBInterface.execute(db, "select * from $table") |> DataFrame
   return r
end

# retrieve item form a table based on a sql condition
retrieve(db, table::String, condition::String )::DataFrame = begin
   r = SQLite.Query( db, "select * from $table where $condition")  |> DataFrame
   return r
end

# run a custom function
runfunct(funct, x, y, z) = funct(x, y, z)


text:
using DiffEqFlux, Flux

"""
    Module that modifies parameters
"""
abstract type Modifier <: DiffEqFlux.FastLayer end

""" Required as a FastLayer, define a paramsum field """
paramlength(f::Modifier) = f.param_size

""" Required as a FastLayer """
initial_params(f::Modifier) = f.initial_params()

""" Nominal modifier does not change the parameters """
struct Nominal{P} <: Modifier
    name::String
    param_size::Int
    initial_params::P
    function Nominal(init = Flux.zeros)
        initial_params() = init(0)
        new{typeof(initial_params)}("Nominal", 0, initial_params)
    end
end

function (f::Nominal)(p_sub, p_link, p)
    return p_sub, p_link
end

function debug(f::Nominal, p_sub, p_link, p)
    return 0.0
end

include("iss_ctrnn0.jl")
include("iss_ctrnn1.jl")

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Documenter, FMIFlux

makedocs(sitename="FMIFlux.jl",
         format = Documenter.HTML(
            collapselevel = 1,
            sidebar_sitename = false
         ),
         pages= Any[
            "Introduction" => "index.md"
            "Examples" => [
                "Examples Overview" => "examples/overview.md"
                "Simple CS-NeuralFMU" => "examples/simple_hybrid_CS.md"
                "Simple ME-NeuralFMU" => "examples/simple_hybrid_ME.md"
                "Modelica Conference 2021" => "examples/modelica_conference_2021.md"
            ]
            "Library Functions" => "library/overview.md"
            "Related Publication" => "related.md"
            "Contents" => "contents.md"
            ]
         )

deploydocs(repo = "github.com/ThummeTo/FMIFlux.jl.git", devbranch = "main")


text:
# # Baroclinic adjustment
#
# In this example, we simulate the evolution and equilibration of a baroclinically
# unstable front.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, CairoMakie"
# ```

using Oceananigans
using Oceananigans.Units

# ## Grid
#
# We use a three-dimensional channel that is periodic in the `x` direction:

Lx = 1000kilometers # east-west extent [m]
Ly = 1000kilometers # north-south extent [m]
Lz = 1kilometers    # depth [m]

Nx = 64
Ny = 64
Nz = 40

grid = RectilinearGrid(size = (Nx, Ny, Nz),
                       x = (0, Lx),
                       y = (-Ly/2, Ly/2),
                       z = (-Lz, 0),
                       topology = (Periodic, Bounded, Bounded))

# ## Turbulence closures

# We prescribe the values of vertical viscocity and diffusivity according to the ratio
# of the vertical and lateral grid spacing.

Δx, Δz = Lx/Nx, Lz/Nz
𝒜 = Δz/Δx # Grid cell aspect ratio.

κh = 0.1    # [m² s⁻¹] horizontal diffusivity
νh = 0.1    # [m² s⁻¹] horizontal viscosity
κz = 𝒜 * κh # [m² s⁻¹] vertical diffusivity
νz = 𝒜 * νh # [m² s⁻¹] vertical viscosity

horizontal_diffusive_closure = HorizontalScalarDiffusivity(ν = νh, κ = κh)

vertical_diffusive_closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization();
                                                       ν = νz, κ = κz)
nothing #hide

# ## Model

# We built a `HydrostaticFreeSurfaceModel` with an `ImplicitFreeSurface` solver.
# Regarding Coriolis, we use a beta-plane centered at 45° South.

model = HydrostaticFreeSurfaceModel(; grid,
                                    coriolis = BetaPlane(latitude = -45),
                                    buoyancy = BuoyancyTracer(),
                                    tracers = :b,
                                    closure = (vertical_diffusive_closure, horizontal_diffusive_closure),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    free_surface = ImplicitFreeSurface())

# We want to initialize our model with a baroclinically unstable front plus some small-amplitude
# noise.

"""
    ramp(y, Δy)

Linear ramp from 0 to 1 between -Δy/2 and +Δy/2.

For example:
```
            y < -Δy/2 => ramp = 0
    -Δy/2 < y < -Δy/2 => ramp = y / Δy
            y >  Δy/2 => ramp = 1
```
"""
ramp(y, Δy) = min(max(0, y/Δy + 1/2), 1)
nothing #hide

# We then use `ramp(y, Δy)` to construct an initial buoyancy configuration of a baroclinically
# unstable front. The front has a buoyancy jump `Δb` over a latitudinal width `Δy`.

N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

Δy = 50kilometers # width of the region of the front
Δb = Δy * M²      # buoyancy jump associated with the front
ϵb = 1e-2 * Δb    # noise amplitude

bᵢ(x, y, z) = N² * z + Δb * ramp(y, Δy) + ϵb * randn()

set!(model, b=bᵢ)

# Let's visualize the initial buoyancy distribution.

using CairoMakie

x, y, z = 1e-3 .* nodes((Center, Center, Center), grid) # convert m -> km

b = model.tracers.b

fig, ax, hm = heatmap(y, z, interior(b)[1, :, :],
                      colormap=:deep,
                      axis = (xlabel = "y [km]",
                              ylabel = "z [km]",
                              title = "b(x=0, y, z, t=0)",
                              titlesize = 24))

Colorbar(fig[1, 2], hm, label = "[m s⁻²]")

save("initial_buoyancy.svg", fig); nothing # hide

# ![](initial_buoyancy.svg)

# Now let's built a `Simulation`.

Δt₀ = 5minutes
stop_time = 40days

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# We add a `TimeStepWizard` callback to adapt the siulation's time-step,

wizard = TimeStepWizard(cfl=0.2, max_change=1.1, max_Δt=20minutes)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# Also, we add a callback to print a message about how the simulation is going,

using Printf

wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(100))

# ## Diagnostics/Output

# Add some diagnostics. Here, we save the buoyancy, ``b``, at the edges of our domain as well as
# the zonal (``x``) averages of buoyancy and zonal velocity ``u``.

u, v, w = model.velocities

B = Field(Average(b, dims=1))
U = Field(Average(u, dims=1))

filename = "baroclinic_adjustment"
save_fields_interval = 0.5day

slicers = (west = (1, :, :),
           east = (grid.Nx, :, :),
           south = (:, 1, :),
           north = (:, grid.Ny, :),
           bottom = (:, :, 1),
           top = (:, :, grid.Nz))

for side in keys(slicers)
    indices = slicers[side]

    simulation.output_writers[side] = JLD2OutputWriter(model, (; b);
                                                       filename = filename * "_$(side)_slice",
                                                       schedule = TimeInterval(save_fields_interval),
                                                       indices)
end

simulation.output_writers[:zonal] = JLD2OutputWriter(model, (b=B,);
                                                     schedule = TimeInterval(save_fields_interval),
                                                     filename = filename * "_zonal_average")

# Now let's run!

@info "Running the simulation..."

run!(simulation)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

# ## Visualization

# Now we are ready to visualize our resutls! We use `CairoMakie` in this example.
# On a system with OpenGL `using GLMakie` is more convenient as figures will be
# displayed on the screen.

using CairoMakie

# We load the saved buoyancy output on the top, bottom, and east surface as `FieldTimeSeries`es.

filename = "baroclinic_adjustment"

sides = keys(slicers)

slice_filenames = NamedTuple(side => filename * "_$(side)_slice.jld2" for side in sides)

b_timeserieses = (east   = FieldTimeSeries(slice_filenames.east, "b"),
                  north  = FieldTimeSeries(slice_filenames.north, "b"),
                  bottom = FieldTimeSeries(slice_filenames.bottom, "b"),
                  top    = FieldTimeSeries(slice_filenames.top, "b"))

avg_b_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "b")

nothing #hide

# We build the coordinates. We rescale horizontal coordinates so that they correspond to kilometers.

x, y, z = nodes(b_timeserieses.east)

x = x .* 1e-3 # convert m -> km
y = y .* 1e-3 # convert m -> km

x_xz = repeat(x, 1, Nz)
y_xz_north = y[end] * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz_east = x[end] * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), grid.Ny, 1)

x_xy = x
y_xy = y
z_xy_top = z[end] * ones(grid.Nx, grid.Ny)
z_xy_bottom = z[1] * ones(grid.Nx, grid.Ny)
nothing #hide

# Then we create a 3D axis. We use `zonal_slice_displacement` to control where the plot of the instantaneous
# zonal average flow is located.

fig = Figure(resolution = (900, 520))

zonal_slice_displacement = 1.2

ax = Axis3(fig[2, 1], aspect=(1, 1, 1/5),
           xlabel="x (km)", ylabel="y (km)", zlabel="z (m)",
           limits = ((x[1], zonal_slice_displacement * x[end]), (y[1], y[end]), (z[1], z[end])),
           elevation = 0.45, azimuth = 6.8,
           xspinesvisible = false, zgridvisible=false,
           protrusions=40,
           perspectiveness=0.7)

nothing #hide

# We use Makie's `Observable` to animate the data. To dive into how `Observable`s work we
# refer to [Makie.jl's Documentation](https://makie.juliaplots.org/stable/documentation/nodes/index.html).

n = Observable(1)

# Now let's make a 3D plot of the buoyancy and in front of it we'll use the zonally-averaged output
# to plot the instantaneous zonal-average of the buoyancy.

b_slices = (east   = @lift(interior(b_timeserieses.east[$n], 1, :, :)),
            north  = @lift(interior(b_timeserieses.north[$n], :, 1, :)),
            bottom = @lift(interior(b_timeserieses.bottom[$n], :, :, 1)),
            top    = @lift(interior(b_timeserieses.top[$n], :, :, 1)))

avg_b = @lift interior(avg_b_timeseries[$n], 1, :, :)

clims = @lift 1.1 .* extrema(b_timeserieses.top[$n][:])

kwargs = (colorrange = clims, colormap = :deep)

surface!(ax, x_yz_east, y_yz, z_yz;    color = b_slices.east, kwargs...)
surface!(ax, x_xz, y_xz_north, z_xz;   color = b_slices.north, kwargs...)
surface!(ax, x_xy, y_xy, z_xy_bottom ; color = b_slices.bottom, kwargs...)
surface!(ax, x_xy, y_xy, z_xy_top;     color = b_slices.top, kwargs...)

sf = surface!(ax, zonal_slice_displacement .* x_yz_east, y_yz, z_yz; color = avg_b, kwargs...)

contour!(ax, y, z, avg_b; transformation = (:yz, zonal_slice_displacement * x[end]),
         levels = 15, linewidth = 2, color = :black)

Colorbar(fig[2, 2], sf, label = "m s⁻²", height = 200, tellheight=false)

# Finally, we add a figure title with the time of the snapshot and then record a movie.

times = avg_b_timeseries.times

title = @lift "Buoyancy at t = " * string(round(times[$n] / day, digits=1)) * " days"

fig[1, 1:2] = Label(fig, title; textsize = 24, tellwidth = false, padding = (0, 0, -120, 0))

frames = 1:length(times)

record(fig, filename * ".mp4", frames, framerate=8) do i
    msg = string("Plotting frame ", i, " of ", frames[end])
    print(msg * " \r")
    n[] = i
end
nothing #hide

# ![](baroclinic_adjustment.mp4)


text:
facts("[nl_linearity] check simplification of formulae") do
    context("ifelse") do
        # First term is true, we should choose `then`
        expr = :(ifelse(1 > 0, x[1] ^ 2, x[2] + 1))
        lin_expr = AmplNLWriter.LinearityExpr(expr)
        @fact AmplNLWriter.convert_formula(lin_expr.c) --> :(x[1] ^ 2)
        @fact lin_expr.linearity --> :nonlinear

        # First term is false, we should choose `else`
        expr = :(ifelse(1 < 0, x[1] ^ 2, x[2] + 1))
        lin_expr = AmplNLWriter.LinearityExpr(expr)
        @fact AmplNLWriter.convert_formula(lin_expr.c) --> :(x[2] + 1)
        @fact lin_expr.linearity --> :linear

        # First term isn't constant, we can't simplify
        expr = :(ifelse(1 > x[1], x[1] ^ 2, x[2] + 1))
        lin_expr = AmplNLWriter.LinearityExpr(expr)
        @fact AmplNLWriter.convert_formula(lin_expr.c) --> expr
        @fact lin_expr.linearity --> :nonlinear
    end
end


text:
export RobotAxisRangeError, RobotDeviceError, RobotDOFError, RobotReferenceError, RobotStateError, RobotTeachError, RobotExplicitCoordinatesError, RobotSafetyError

abstract type RobotException <: Exception end

struct RobotStateError <: RobotException
    robot::Robot
    state::Union{RobotState,Nothing}
end

struct RobotAxisRangeError <: RobotException
    robot::Robot
    pos::RobotCoords
end

struct RobotDOFError <: RobotException
    robot::Robot
    dof::Int
end

struct RobotReferenceError <: RobotException
    robot::Robot
end

struct RobotDeviceError <: RobotException
    robot::Robot
    exc::Exception
end

struct RobotTeachError <: RobotException
    robot::Robot
    pos_name::String
end

struct RobotExplicitCoordinatesError <: RobotException
    robot::Robot
end

struct RobotSafetyError <: Exception
    robot::Robot
    pos::RobotCoords
end

function Base.showerror(io::IO, ex::RobotStateError)
    if ex.state === nothing
        print(io, "RobotStateError: robot '$(deviceID(ex.robot))' is in state $(state(ex.robot))")
    else
        print(io, "RobotStateError: robot '$(deviceID(ex.robot))' has to be in state $(ex.state) to perform the desired action. It was in state $(state(ex.robot)) instead.")
    end
end

Base.showerror(io::IO, ex::RobotAxisRangeError) = (print(io, "RobotAxisRangeError: position ");
                                               show(IOContext(io, :typeinfo => typeof(ex.pos)), ex.pos);
                                               print(io, " (ScannerCoords: ");
                                               show(IOContext(io, :typeinfo => typeof(toScannerCoords(ex.robot, ex.pos))), toScannerCoords(ex.robot, ex.pos));
                                               print(io, ") is out of the axis range ");
                                               show(IOContext(io, :typeinfo => typeof(axisRange(ex.robot))), axisRange(ex.robot));
                                               print(io, " of robot '$(deviceID(ex.robot))'"))

Base.showerror(io::IO, ex::RobotDOFError) = print(io, "RobotDOFError: coordinates included $(ex.dof) axes, but robot '$(deviceID(ex.robot))' has $(dof(ex.robot)) degrees-of-freedom")

Base.showerror(io::IO, ex::RobotReferenceError) = print(io, "RobotReferenceError: robot '$(deviceID(ex.robot))' has to be referenced to perform this action")

Base.showerror(io::IO, ex::RobotDeviceError) = (println(io, "RobotDeviceError: during the communication with robot '$(deviceID(ex.robot))' the following error occured:"); showerror(io, ex.exc))

function Base.showerror(io::IO, ex::RobotTeachError)
    if haskey(namedPositions(ex.robot), ex.pos_name)
        print(io, "RobotTeachError: the desired position name $(ex.pos_name) is already defined for robot '$(deviceID(ex.robot))', to override the current value pass override=true")
    else
        print(io, "RobotTeachError: the desired position $(ex.pos_name) is not defined for robot '$(deviceID(ex.robot))', use one of $(keys(namedPositions(ex.robot)))")
    end
end

Base.showerror(io::IO, ex::RobotExplicitCoordinatesError) = print(io, "RobotExplicitCoordinatesError: robot '$(deviceID(ex.robot))' has defined a coordinate transform and therefore needs explicit declaration of the used coordinate system. Please pass RobotCoords(pos) or ScannerCoords(pos).")

Base.showerror(io::IO, ex::RobotSafetyError) = (print(io, "RobotSafetyError: position ");
                                                show(IOContext(io, :typeinfo => typeof(ex.pos)), ex.pos);
                                                print(io, " (ScannerCoords: ");
                                                show(IOContext(io, :typeinfo => typeof(toScannerCoords(ex.robot, ex.pos))), toScannerCoords(ex.robot, ex.pos));
                                                print(io, ") is not safe to move to!"); )


text:
"""
    gather!(dst, src, idx)

Reverse operation of [`scatter!`](@ref). Gathers data from source `src`
and writes it in destination `dst` according to the index array `idx`.
For each `k` in `CartesianIndices(idx)`, assign values to `dst` according to

    dst[:, ... , k] .= src[:, ... , idx[k]...]

Notice that if `idx` is a vector containing integers,
and both `dst` and `src` are matrices, previous expression simplifies to

    dst[:, k] .= src[:, idx[k]]

and `k` will run over `1:length(idx)`.

The elements of `idx` can be integers or integer tuples and may be repeated.
A single `src` column can end up being copied into zero, one,
or multiple `dst` columns.

See [`gather`](@ref) for an allocating version.
"""
function gather!(dst::AbstractArray, src::AbstractArray, idx::AbstractArray)
    dims = scatter_dims(src, dst, idx)
    colons = ntuple(i -> Colon(), dims)
    for k in CartesianIndices(idx)
        _view(dst, colons, k) .= _view(src, colons, idx[k])
    end
    return dst
end

"""
    gather(src, idx) -> dst

Reverse operation of [`scatter`](@ref). Gathers data from source `src`
and writes it in a destination `dst` according to the index
array `idx`.
For each `k` in `CartesianIndices(idx)`, assign values to `dst`
according to

    dst[:, ... , k] .= src[:, ... , idx[k]...]

Notice that if `idx` is a vector containing integers
and `src` is a matrix, previous expression simplifies to

    dst[:, k] .= src[:, idx[k]]

and `k` will run over `1:length(idx)`.

The elements of `idx` can be integers or integer tuples and may be repeated.
A single `src` column can end up being copied into zero, one,
or multiple `dst` columns.

See [`gather!`](@ref) for an in-place version.
"""
function gather(src::AbstractArray{Tsrc, Nsrc},
                idx::AbstractArray{Tidx, Nidx}) where
                    {Tsrc, Nsrc, Nidx, Tidx}

    M = typelength(Tidx)
    dstsize = (size(src)[1:Nsrc-M]..., size(idx)...)
    dst = similar(src, Tsrc, dstsize)
    return gather!(dst, src, idx)
end

# Simple implementation with getindex for integer array.
# Perf equivalent to the one above (which can also handle the integer case)
# leave it here to show the simple connection with getindex.
function gather(src::AbstractArray{Tsrc, Nsrc},
                idx::AbstractArray{<:Integer}) where {Tsrc, Nsrc}
    colons = ntuple(i -> Colon(), Nsrc-1)
    return src[colons..., idx]
end

# Gradient

∇gather_src(Δ, src_size, idx) = scatter!(+, fill!(similar(Δ, eltype(Δ), src_size), 0), Δ, idx)

function rrule(::typeof(gather!), dst::AbstractArray, src::AbstractArray, idx::AbstractArray)
    y = gather!(dst, src, idx)
    src_size = size(src)
    gather!_pullback(Δ) = (NO_FIELDS, NoTangent(), ∇gather_src(Δ, src_size, idx), NoTangent())
    y, gather!_pullback
end


text:
using Test
using TestSetExtensions

using Oscar

include("../experiments_bounds.jl")

#------- Auxiliary functions --------------------------------------------------

function diff_sol_Lie_derivatives(ode::ODE, outputs, params, ic, inputs, prec::Int)
    # creating a new ring with variables for the derivatives of u
    new_varnames = map(string, gens(ode.poly_ring))
    if length(ode.u_vars) > 0
        append!(new_varnames, ["$(u)_$i" for u in ode.u_vars for i in 1:(prec - 1)])
    end
    new_ring, vars = PolynomialRing(base_ring(ode.poly_ring), new_varnames)

    # mapping everything to the new ring
    eval_point = [str_to_var(string(v), new_ring) for v in gens(ode.poly_ring)]
    new_eqs = Dict()
    for (x, f) in ode.equations
        if applicable(f, numerator)
            new_eqs[str_to_var(string(x), new_ring)] = evaluate(numerator(f), eval_point) // evaluate(denominator(f), eval_point)
        else
            new_eqs[str_to_var(string(x), new_ring)] = evaluate(f, eval_point)
        end
    end
    outputs_new = []
    for g in outputs
        if applicable(g, numerator)
            push!(outputs_new, evaluate(numerator(g), eval_point) // evaluate(denominator(g), eval_point))
        else
            push!(outputs_new, evaluate(g, eval_point))
        end
    end
    outputs = outputs_new
    params, ic = map(d -> Dict(str_to_var(string(k), new_ring) => v for (k, v) in d), [params, ic])

    # computing Lie derivatives
    derivation = copy(new_eqs)
    for u in ode.u_vars
        derivation[str_to_var(string(u), new_ring)] = str_to_var("$(u)_1", new_ring)
        for i in 1:(prec - 2)
            derivation[str_to_var("$(u)_$i", new_ring)] = str_to_var("$(u)_$(i + 1)", new_ring)
        end
    end
    Lie_derivatives = []
    for g in outputs
        push!(Lie_derivatives, [g // new_ring(1)])
        for i in 1:prec
            push!(
                Lie_derivatives[end],
                sum([derivative(Lie_derivatives[end][end], v) * get(derivation, v, 0) for v in gens(new_ring)])
            )
        end
    end

    # producing the result
    eval_dict = merge(params, ic)
    for u in ode.u_vars
        eval_dict[str_to_var(string(u), new_ring)] = inputs[u][1]
        for i in 2:prec
            eval_dict[str_to_var("$(u)_$(i - 1)", new_ring)] = inputs[u][i] * factorial(i - 1)
        end
    end
    eval_vec = [eval_dict[v] for v in gens(new_ring)]

    result = []
    for i in 1:length(outputs)
        push!(result, Dict())
        for v in vcat(ode.x_vars, ode.parameters)
            result[end][v] = []
            for j in 1:prec
                push!(
                    result[end][v],
                    eval_frac(derivative(Lie_derivatives[i][j], str_to_var("$v", new_ring)), eval_vec)
                )
            end
        end
    end

    return result
end

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

function rand_poly(deg, vars)
    if deg == 0
        return parent(vars[1])(1)
    end
    result = 0
    indices = collect(1:length(vars))
    monomials = []
    for d in 0:deg
        for subs in IterTools.subsets(indices, d)
            push!(monomials, subs)
        end
    end

    for subs in monomials
        monom = rand(-50:50)
        for v_ind in subs
            monom *= vars[v_ind]
        end
        result += monom
    end

    return result
end

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

@info "Testing started"

@testset "All the tests" begin
    @includetests ARGS
end

max_stars_repo_name: pogudingleb/ExperimentsBound

max_issues_repo_name: pogudingleb/ExperimentsBound

max_forks_repo_name: pogudingleb/ExperimentsBound


text:
module Mesh

include("BrickMesh.jl")
include("Topologies.jl")
include("Metrics.jl")
include("Elements.jl")
include("Grids.jl")
include("Filters.jl")
include("Geometry.jl")

end # module

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# Classification Network
struct ClassifierNODE{L1,L2,L3,P1,P2,T}
    preode::P1
    node::L2
    postode::P2
    p1::T
    p2::T
    p3::T

    function ClassifierNODE(preode, node, postode)
        p1, re1 = Flux.destructure(preode)
        p2 = node.p
        p3, re3 = Flux.destructure(postode)
        return new{
            typeof(preode),
            typeof(node),
            typeof(postode),
            typeof(re1),
            typeof(re3),
            typeof(p2),
        }(
            re1,
            node,
            re3,
            p1,
            p2,
            p3,
        )
    end
end

Flux.trainable(m::ClassifierNODE) = (m.p1, m.p2, m.p3)

function (m::ClassifierNODE{L1,L2,L3})(
    x,
    p1 = m.p1,
    p2 = m.p2,
    p3 = m.p3;
    node_kwargs...,
) where {L1,L2,L3}
    preode = m.preode(p1)::L1
    x = preode(x)
    x, nfe, sv = m.node(x, p2; node_kwargs...)
    postode = m.postode(p3)::L3
    return postode(x), nfe, sv
end

# Classification Network
struct ClassifierNSDE{L1,L2,L3,P1,P2,T}
    presde::P1
    nsde::L2
    postsde::P2
    p1::T
    p2::T
    p3::T

    function ClassifierNSDE(presde, nsde, postsde)
        p1, re1 = Flux.destructure(presde)
        p2 = nsde.p
        p3, re3 = Flux.destructure(postsde)
        return new{
            typeof(presde),
            typeof(nsde),
            typeof(postsde),
            typeof(re1),
            typeof(re3),
            typeof(p2),
        }(
            re1,
            nsde,
            re3,
            p1,
            p2,
            p3,
        )
    end
end

Flux.trainable(m::ClassifierNSDE) = (m.p1, m.p2, m.p3)

function (m::ClassifierNSDE{L1,L2,L3})(
    x,
    p1 = m.p1,
    p2 = m.p2,
    p3 = m.p3;
    trajectories::Int = 10,
    nsde_kwargs...,
) where {L1,L2,L3}
    # Neural SDE right now works only with CPU Inputs
    bsize = size(x)[end]
    x = _expand(x, trajectories) |> track
    presde = m.presde(p1)::L1
    x = presde(x)
    x, nfe1, nfe2, sv = m.nsde(x, p2; nsde_kwargs...)
    postsde = m.postsde(p3)::L3
    z = postsde(x)
    z = reshape(mean(reshape(z, size(z, 1), :, trajectories), dims = 3), :, bsize)
    return z, nfe1, nfe2, sv
end

_expand(x::AbstractArray{<:Number,4}, d::Int) = repeat(x, 1, 1, 1, d)
_expand(x::AbstractArray{<:Number,2}, d::Int) = repeat(x, 1, d)


text:
module TwiliteTimeline

using OAuth, JSON

export Authentictor, ResourceParams, collect_tweets
export ParamsGetTweets, ParamsPostTweet
export get_tweets, post_tweet
export write_to_json, read_from_json

include("collector.jl")
include("persistor.jl")
include("parameters.jl")
include("callers.jl")

end


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

mutable struct BgpSettings <: SwaggerModel
    asn::Any # spec type: Union{ Nothing, Int64 } # spec name: asn
    bgpPeeringAddress::Any # spec type: Union{ Nothing, String } # spec name: bgpPeeringAddress
    peerWeight::Any # spec type: Union{ Nothing, Int32 } # spec name: peerWeight

    function BgpSettings(;asn=nothing, bgpPeeringAddress=nothing, peerWeight=nothing)
        o = new()
        validate_property(BgpSettings, Symbol("asn"), asn)
        setfield!(o, Symbol("asn"), asn)
        validate_property(BgpSettings, Symbol("bgpPeeringAddress"), bgpPeeringAddress)
        setfield!(o, Symbol("bgpPeeringAddress"), bgpPeeringAddress)
        validate_property(BgpSettings, Symbol("peerWeight"), peerWeight)
        setfield!(o, Symbol("peerWeight"), peerWeight)
        o
    end
end # type BgpSettings

const _property_map_BgpSettings = Dict{Symbol,Symbol}(Symbol("asn")=>Symbol("asn"), Symbol("bgpPeeringAddress")=>Symbol("bgpPeeringAddress"), Symbol("peerWeight")=>Symbol("peerWeight"))
const _property_types_BgpSettings = Dict{Symbol,String}(Symbol("asn")=>"Int64", Symbol("bgpPeeringAddress")=>"String", Symbol("peerWeight")=>"Int32")
Base.propertynames(::Type{ BgpSettings }) = collect(keys(_property_map_BgpSettings))
Swagger.property_type(::Type{ BgpSettings }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_BgpSettings[name]))}
Swagger.field_name(::Type{ BgpSettings }, property_name::Symbol) =  _property_map_BgpSettings[property_name]

function check_required(o::BgpSettings)
    true
end

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


text:
#Only really going to work for xy plots?
function layer(plot1::VegaVisualization, plot2::VegaVisualization)

    #Make copy of plots as to not modify the inputs
    baseplot = deepcopy(plot1)
    overlay = deepcopy(plot2)

    #Stash plot2 main "group" scale into inner Vega group mark scale
    #This allows for different colors for every layer, presuming you set the color prior to layering
    #Still need to write code to avoid overwriting value if it's already there
    i = findfirst([z.name == "group" for z in overlay.scales])
    for m in 1:length(overlay.marks)
        overlay.marks[m].scales = [overlay.scales[i]]
    end

    #Iterate over data tables in second plot, append to baseplot
    #This should be fine as-is
    for i in 1:length(overlay.data)
        push!(baseplot.data, overlay.data[i])
    end

    #Add all marks to baseplot
    #This should be fine as-is
    for i in 1:length(overlay.marks)
        push!(baseplot.marks, overlay.marks[i])
    end

    #Make joint x and y scales
    xloc1 = findfirst([z.name == "x" for z in baseplot.scales])
    xloc2 = findfirst([z.name == "x" for z in overlay.scales])
    yloc1 = findfirst([z.name == "y" for z in baseplot.scales])
    yloc2 = findfirst([z.name == "y" for z in overlay.scales])

    #Overwrite baseplot domain value
    if baseplot.scales[xloc1].domain.fields == nothing
        xdomain1=baseplot.scales[xloc1].domain
    else
        xdomain1=baseplot.scales[xloc1].domain.fields
    end
    if overlay.scales[xloc2].domain.fields == nothing
        xdomain2=overlay.scales[xloc2].domain
    else
        xdomain2=overlay.scales[xloc2].domain.fields
    end

    if baseplot.scales[yloc1].domain.fields == nothing
        ydomain1=baseplot.scales[yloc1].domain
    else
        ydomain1=baseplot.scales[yloc1].domain.fields
    end
    if overlay.scales[yloc2].domain.fields == nothing
        ydomain2=overlay.scales[yloc2].domain
    else
        ydomain2=overlay.scales[yloc2].domain.fields
    end

    baseplot.scales[xloc1].domain = VegaDataRef(fields =
        [xdomain1, xdomain2;])
    baseplot.scales[yloc1].domain = VegaDataRef(fields =
        [ydomain1, ydomain2;])

    return baseplot

end

max_stars_repo_name: jpthor/vega


text:
"""
    FromState

a basic thermodynamic model, used to extract thermodynamic specifications from a state without using a Equation of State.
"""
struct FromState <: ThermoModel end
const URVec = AbstractVector{T} where T<:Number

moles2(a::ThermodynamicState,mw) = moles3(amount_type(a),a,mw)
mass2(a::ThermodynamicState,mw) =  mass3(amount_type(a),a,mw)
kg_per_mol2(a::ThermodynamicState,mw) =  kg_per_mol3(amount_type(a),a,mw)

function mw_error(::Nothing)
    throw(error("molecular weight information is required for this combination of properties."))
end

function mw_error(mw)
    return nothing
end
##one mol cases

#those cases should rarely be called
function moles3(::Tuple{T,OneMol},st,mw::Nothing) where T
    return one(numtype(st))
end

function moles3(::Tuple{T,OneMol},st,mw::T2) where {T,T2 <:Number}
    return one(promote_type(T2,numtype(st)))
end

#error?
function moles3(::Tuple{T,OneMol},st,mw::T2) where T where T2<:AbstractVector{T3} where T3 <: Number
    return one(promote_type(T3,numtype(st)))
end

function moles3(::Tuple{SingleComponent,MaterialAmount{MOLAR}},st,mw::T) where T <:Number
    return normalize_units(value(get_spec(MaterialAmount{MOLAR}(),st)))*_one(st)
end

function moles3(::Tuple{SingleComponent,MaterialAmount{MASS}},st,mw::T) where T <:Number
    # g/(g/mol) = mol
    mw_error(mw)
    x = mw_div(value(get_spec(MaterialAmount{MASS}(),st)),mw)
    return normalize_units(x)*_one(st)
end

## total ammounts:
function moles3(::Tuple{T,T},st,mw) where T<:MaterialCompounds{MOLAR,TOTAL_AMOUNT}
    return normalize_units(sum(value(get_spec(T(),st))))*_one(st)
end

function moles3(::Tuple{T,T},st,mw::T2) where T<:MaterialCompounds{MASS,TOTAL_AMOUNT} where T2<: URVec
    ###### needs mol
    ### exists mass numbers
    #mw_i * xmi = gc* (molc/gc)= molc, molmix = sum(molc)
    mw_error(mw)
    compounds::Spec = get_spec(T(),st)
    v = compounds.val
    sum_n_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_n_mw += mw_div(v[i],mw[i])
    end
    #sum_n_mw = mapreduce(mw_div,+,compounds,mw)
    return normalize_units(sum_n_mw)
end

function moles3(::Tuple{T,MaterialAmount{MOLAR}},st,mw) where T
    return normalize_units(value(get_spec(MaterialAmount{MOLAR}(),st)))*_one(st)
end

#inverse operations with fractions
function moles3(::Tuple{MaterialCompounds{MOLAR,FRACTION},MaterialAmount{MASS}},st,mw::T) where T<:URVec
    #mol fraction and mass
    #xni * mwi = mc/mmix * (gc/molc)= gc/molmix (sum)-> gmix/molmix
    #molmix = gmix / (gmix/molmix)
    mw_error(mw)
    compounds::Spec = get_spec(MaterialCompounds{MOLAR,FRACTION}(),st)
    amount::Spec = get_spec(MaterialAmount{MASS}(),st)
    v = compounds.val
    #sum_x_mw = mapreduce(mw_mul,+,value(compounds),mw)
    sum_x_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_x_mw += mw_mul(v[i],mw[i])
    end
    return normalize_units(value(amount)/sum_x_mw)
end

function moles3(::Tuple{MaterialCompounds{MASS,FRACTION},MaterialAmount{MASS}},st,mw::T) where T<:URVec
    #mass fractions and mass
    #mw_i * xmi = gc/gmix / (molc/gc)= molc/gmix (sum)-> molmix/gmix
    # molmix = molmix/gmix * (gmix)
    mw_error(mw)
    compounds::Spec = get_spec(MaterialCompounds{MASS,FRACTION}(),st)
    amount::Spec  = get_spec(MaterialAmount{MASS}(),st)
    v = compounds.val
    sum_x_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_x_mw += mw_div(v[i],mw[i])
    end
    #sum_x_mw = mapreduce(mw_div,+,value(compounds),mw)
    return normalize_units(value(amount)*sum_x_mw)
end

function mass3(::Tuple{SingleComponent,OneMol},st,mw::T2) where T2 <:Number
    mw_error(mw)
    return normalize_units(mw_mul(one(T2),mw))*_one(st)
end

#mass defined
function mass3(::Tuple{T,MaterialAmount{MASS}},st,mw) where T
    return normalize_units(value(get_spec(MaterialAmount{MASS}(),st)))*_one(st)
end

function mass3(::Tuple{SingleComponent,MaterialAmount{MOLAR}},st,mw) where T
    mw_error(mw)
    mol = value(get_spec(MaterialAmount{MOLAR}(),st))
    return normalize_units(mw_mul(mol,mw))*_one(st)
end

## total ammounts:
function mass3(::Tuple{T,T},st,mw) where T<:MaterialCompounds{MASS,TOTAL_AMOUNT}
    return normalize_units(sum(value(get_spec(T(),st))))*_one(st)
end

function mass3(::Tuple{T,T},st::S,mw::T2) where {S<:ThermodynamicState,T<:MaterialCompounds{MOLAR,TOTAL_AMOUNT},T2<:URVec}
    ###### needs mass
    ### exists mol numbers
    #mw_i * ni = (gc/molc) * molc = gc, sum(gc) = gmix
    mw_error(mw)
    compounds = get_spec(T(),st)::Spec
    v = compounds.val
    sum_n_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_n_mw += mw_mul(v[i],mw[i])
    end
    return normalize_units(sum_n_mw)
end

#inverse operations with fractions
function mass3(::Tuple{MaterialCompounds{MOLAR,FRACTION},MaterialAmount{MOLAR}},st,mw::T2) where T2<:URVec
    #moles, mol fractions
    #mw_i * xi = (gc/molc) * molc/molmix = gc/molmix (*molmix)
    mw_error(mw)
    amount = get_spec(MaterialAmount{MOLAR}(),st)
    compounds::Spec = get_spec(MaterialCompounds{MOLAR,FRACTION}(),st)
    sum_x_mw = zero(numtype(st))
    v = compounds.val
    for i in 1:length(mw)
        sum_x_mw += mw_mul(v[i],mw[i])
    end
    return normalize_units(value(amount)*sum_x_mw)
end

function mass3(::Tuple{MaterialCompounds{MASS,FRACTION},MaterialAmount{MOLAR}},st,mw::T2) where T2<:URVec
    #mass fractions and moles
    #sum_x_mw = xmi/mwi = mc/mmix / mc/molc = mc/mmix * molc/mc = molc/mmix -> sum -> molmix/mmix
    #a = molmix/mmix
    #inv(a) = mmix/molmix
    # molmix = molmix/gmix * (gmix)
    mw_error(mw)
    compounds::Spec = get_spec(MaterialCompounds{MASS,FRACTION}(),st)
    v = compounds.val
    amount = get_spec(MaterialAmount{MOLAR}(),st)
    sum_x_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_x_mw += mw_div(v[i],mw[i])
    end
    return normalize_units(value(amount)/sum_x_mw)
end

#single component cases, all the same
function kg_per_mol3(::Tuple{SingleComponent,OneMol},st,mw::T2) where T2 <:Number
    return normalize_units(mw_mul(one(T2),mw))*_one(st)
end

function kg_per_mol3(::Tuple{SingleComponent,MaterialAmount},st,mw::T2) where {T2 <:Number}
    return normalize_units(mw_mul(one(T2),mw))*_one(st)
end

#total amounts:
function kg_per_mol3(::Tuple{T,T},st,mw) where T<:MaterialCompounds{MASS,TOTAL_AMOUNT}
    #mw from mass numbers
    #gc / (gc/molc) = molc sum-> molmix
    mw_error(mw)
    compounds::Spec = get_spec(T(),st)
    v = compounds.val
    molmix = zero(numtype(st))
    for i in 1:length(mw)
        molmix += mw_div(v[i],mw[i])
    end
    gmix = sum(v)
    return normalize_units(gmix/molmix)
end

function kg_per_mol3(::Tuple{T,T},st,mw) where T<:MaterialCompounds{MOLAR,TOTAL_AMOUNT}
    #mw from mol numbers
    #molc * gc/molc = gc sum-> gmix
    mw_error(mw)
    compounds::Spec = get_spec(T(),st)
    v = compounds.val
    gmix = zero(numtype(st))
    for i in 1:length(mw)
        gmix += mw_mul(v[i],mw[i])
    end
    #gmix = mapreduce(mw_mul,+,value(compounds),mw)
    molmix = sum(value(compounds))
    return normalize_units(gmix/molmix)*_one(st)
end

#fraction amounts
function kg_per_mol3(::Tuple{MaterialCompounds{MOLAR,FRACTION},MaterialAmount{T}},st,mw) where T
    #mol fraction and mass
    #xni * mwi = mc/mmix * (gc/molc)= gc/molmix (sum)-> gmix/molmix
    mw_error(mw)
    compounds::Spec = get_spec(MaterialCompounds{MOLAR,FRACTION}(),st)
    #sum_x_mw = mapreduce(mw_mul,+,value(compounds),mw)*_one(st)
    v = compounds.val
    sum_x_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_x_mw += mw_mul(v[i],mw[i])
    end
    return normalize_units(one(sum_x_mw)/sum_x_mw)
end

function kg_per_mol3(::Tuple{MaterialCompounds{MASS,FRACTION},MaterialAmount{T}},st,mw) where T
    mw_error(mw)
    #mass fractions and mass
    compounds::Spec = get_spec(MaterialCompounds{MASS,FRACTION}(),st)
    #xmi * mwi = gc/gmix / gc/molc = molc/gmix (sum)-> molmix/gmix
    #sum_x_mw = mapreduce(mw_div,+,value(compounds),mw)*_one(st)
    v = compounds.val
    sum_x_mw = zero(numtype(st))
    for i in 1:length(mw)
        sum_x_mw += mw_div(v[i],mw[i])
    end
    return normalize_units(one(sum_x_mw)/sum_x_mw)
end

#invariant
function to_spec(st,sp::Spec{SP},mw,::MOLAR) where {SP<:AbstractIntensiveSpec{MOLAR}}
    return normalize_units(value(sp))*_one(st)
end

function to_spec(st,sp::Spec{SP},mw,any_value) where {SP<:Union{Pressure,Temperature}}
    return normalize_units(value(sp))*_one(st)
end

function to_spec(st,sp::Spec{SP},mw,::MASS) where {SP<:AbstractIntensiveSpec{MASS}}
    return normalize_units(value(sp))*_one(st)
end

function to_spec(st,sp::Spec{SP},mw,::TOTAL) where {SP<:AbstractIntensiveSpec{TOTAL}}
    return normalize_units(value(sp))*_one(st)
end

#to mass

#from molar to mass
function to_spec(st,sp::Spec{SP},mw,::MASS) where {SP<:AbstractIntensiveSpec{MOLAR}}
    return normalize_units(value(sp))/kg_per_mol2(st,mw)*_one(st)
end

function to_spec(st,sp::Spec{SP},mw,::MASS) where {SP<:AbstractIntensiveSpec{TOTAL}}
    return normalize_units(value(sp))/mass2(st,mw)*_one(st)
end

#to mol
function to_spec(st,sp::Spec{SP},mw,::MOLAR) where {SP<:AbstractIntensiveSpec{MASS}}
    return normalize_units(value(sp))*kg_per_mol2(st,mw)
end

function to_spec(st,sp::Spec{SP},mw,::MOLAR) where {SP<:AbstractIntensiveSpec{TOTAL}}
    return normalize_units(value(sp))/moles2(st,mw)
end

#to total
function to_spec(st,sp::Spec{SP},mw,::TOTAL) where {SP<:AbstractIntensiveSpec{MOLAR}}
    return normalize_units(value(sp))*moles2(st,mw)
end

function to_spec(st,sp::Spec{SP},mw,::TOTAL) where {SP<:AbstractIntensiveSpec{MASS}}
    return normalize_units(value(sp))*mass2(st,mw)
end

#volume and density
#same value
function to_spec_vol(st,sp::Spec{T},mw,::T) where {T<:VolumeAmount}
    return normalize_units(value(sp))
end

#inversion

function to_spec_vol(st,sp::Spec{VolumeAmount{T,VOLUME}},mw,::VolumeAmount{T,DENSITY}) where T
    val = value(sp)
    return normalize_units(one(val)/val)
end

function to_spec_vol(st,sp::Spec{VolumeAmount{T,DENSITY}},mw,::VolumeAmount{T,VOLUME}) where T
    val = value(sp)
    return normalize_units(one(val)/val)
end

#same type, volume
function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,VOLUME}},mw,::VolumeAmount{MASS,VOLUME})
    val = value(sp)
    return normalize_units(val/kg_per_mol2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,VOLUME}},mw,::VolumeAmount{MOLAR,VOLUME})
    val = value(sp)
    return normalize_units(val*kg_per_mol2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,VOLUME}},mw,::VolumeAmount{TOTAL,VOLUME})
    val = value(sp)
    return normalize_units(val*moles2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,VOLUME}},mw,::VolumeAmount{MOLAR,VOLUME})
    val = value(sp)
    return normalize_units(val/moles2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,VOLUME}},mw,::VolumeAmount{TOTAL,VOLUME})
    val = value(sp)
    return normalize_units(val*mass2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,VOLUME}},mw,::VolumeAmount{MASS,VOLUME})
    val = value(sp)
    return normalize_units(val/mass2(st,mw))
end

#same type, density
function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,DENSITY}},mw,::VolumeAmount{MASS,DENSITY})
    val = value(sp)
    return normalize_units(val*kg_per_mol2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,DENSITY}},mw,::VolumeAmount{MOLAR,DENSITY})
    val = value(sp)
    return normalize_units(val/kg_per_mol2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,DENSITY}},mw,::VolumeAmount{TOTAL,DENSITY})
    val = value(sp)
    return normalize_units(val/moles2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,DENSITY}},mw,::VolumeAmount{MOLAR,DENSITY})
    val = value(sp)
    return  (val*moles2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,DENSITY}},mw,::VolumeAmount{TOTAL,DENSITY})
    val = value(sp)
    return normalize_units(val/mass2(st,mw))
end

function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,DENSITY}},mw,::VolumeAmount{MASS,DENSITY})
    val = value(sp)
    return normalize_units(val*mass2(st,mw))
end

#all different combinations
function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,VOLUME}},mw,::VolumeAmount{TOTAL, DENSITY})
    val = value(sp)
    totv = val*moles2(st,mw)
    return normalize_units(one(totv)/totv)
  end

function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,VOLUME}},mw,::VolumeAmount{MASS, DENSITY})
    val = value(sp)
    val2 = val/kg_per_mol2(st,mw)
    return normalize_units(one(val2)/val2)
end

function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,VOLUME}},mw,::VolumeAmount{MOLAR, DENSITY})
    val = value(sp)
    val2 = val/moles2(st,mw)
    return normalize_units(one(val2)/val2)
end
function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,VOLUME}},mw,::VolumeAmount{MASS, DENSITY})
    val = value(sp)
    val2 = val/mass2(st,mw)
    return normalize_units(one(val2)/val2)
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,VOLUME}},mw,::VolumeAmount{MOLAR, DENSITY})
    val = value(sp)
    val2 = val*kg_per_mol2(st,mw)
    return normalize_units(one(val2)/val2)
end
function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,VOLUME}},mw,::VolumeAmount{TOTAL, DENSITY})
    val = value(sp)
    val2 = val*mass2(st,mw)
    return normalize_units(one(val2)/val2)
end

function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL,DENSITY}},mw,::VolumeAmount{MOLAR, VOLUME})
    val = value(sp)
    val2 = val*moles2(st,mw)
    return normalize_units(one(val2)/val2)
end
function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,DENSITY}},mw,::VolumeAmount{MOLAR, VOLUME})
    val = value(sp)
    val2 = val/kg_per_mol2(st,mw)
    return normalize_units(one(val2)/val2)
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,DENSITY}},mw,::VolumeAmount{TOTAL, VOLUME})
    val = value(sp)
    val2 = val/moles2(st,mw)
    return normalize_units(one(val2)/val2)
end
function to_spec_vol(st,sp::Spec{VolumeAmount{MASS,DENSITY}},mw,::VolumeAmount{TOTAL, VOLUME})
    val = value(sp)
    val2 = val/mass2(st,mw)
    return normalize_units(one(val2)/val2)
end

function to_spec_vol(st,sp::Spec{VolumeAmount{MOLAR,DENSITY }},mw,::VolumeAmount{MASS, VOLUME})
    val = value(sp)
    val2 = val*kg_per_mol2(st,mw)
    return normalize_units(one(val2)/val2)
end
function to_spec_vol(st,sp::Spec{VolumeAmount{TOTAL, DENSITY }},mw,::VolumeAmount{MASS, VOLUME})
    val = value(sp)
    val2 = val*mass2(st,mw)
    return normalize_units(one(val2)/val2)
end

function to_spec(st,sp::Spec{SP1},mw,x::SP2) where {SP1<:VolumeAmount,SP2<:VolumeAmount}
    return to_spec_vol(st,sp,mw,x)
end

#material compounds part

#fraction <-> total part, same base
function to_spec_mat(st,sp::Spec{MaterialCompounds{MOLAR,FRACTION}},mw,::MaterialCompounds{MOLAR,TOTAL_AMOUNT})
    return moles2(st,mw) .* normalize_units(value(sp))
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MASS,FRACTION}},mw,::MaterialCompounds{MASS,TOTAL_AMOUNT})
    return mass2(st,mw) .* normalize_units(value(sp))
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MOLAR,TOTAL_AMOUNT}},mw,::MaterialCompounds{MOLAR,FRACTION})
    return normalize_units(value(sp)) ./ moles2(st,mw)
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MASS,TOTAL_AMOUNT}},mw,::MaterialCompounds{MASS,FRACTION})
    return normalize_units(value(sp)) ./ mass2(st,mw)
end

#molar total <-> mass total
function to_spec_mat(st,sp::Spec{MaterialCompounds{MOLAR,TOTAL_AMOUNT}},mw,::MaterialCompounds{MASS,TOTAL_AMOUNT})
    return normalize_units(map(mw_mul,value(sp),mw))
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MASS,TOTAL_AMOUNT}},mw,::MaterialCompounds{MOLAR,TOTAL_AMOUNT})
    return  normalize_units(map(mw_div,value(sp),mw))
end

#molar fraction <-> mass fraction
function to_spec_mat(st,sp::Spec{MaterialCompounds{MOLAR,FRACTION}},mw,::MaterialCompounds{MASS,FRACTION})
    n =  normalize_units(value(sp)) .* moles2(st,mw)
    m = normalize_units(map(mw_mul,n,mw))
    return m ./ sum(m)
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MASS,FRACTION}},mw,::MaterialCompounds{MOLAR,FRACTION})
    m =  normalize_units(value(sp)) .* mass2(st,mw)
    n = map(mw_div,m,mw)
    return n ./ sum(n)
end

#crossed properties:

function to_spec_mat(st,sp::Spec{MaterialCompounds{MOLAR,FRACTION}},mw,::MaterialCompounds{MASS,TOTAL_AMOUNT})
    n =  value(sp).* moles2(st,mw)
    m = normalize_units(map(mw_mul,n,mw))
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MASS,FRACTION}},mw,::MaterialCompounds{MOLAR,TOTAL_AMOUNT})
    m =  value(sp) .* mass2(st,mw)
    n = normalize_units(map(mw_div,m,mw))
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MOLAR,TOTAL_AMOUNT}},mw,::MaterialCompounds{MASS,FRACTION})
    mw_error(mw)
    m =  map(mw_mul,value(sp),mw)
    return normalize_units(m ./ sum(m))
end

function to_spec_mat(st,sp::Spec{MaterialCompounds{MASS,TOTAL_AMOUNT}},mw,::MaterialCompounds{MOLAR,FRACTION})
    mw_error(mw)
    n =  map(mw_div,value(sp),mw)
    return normalize_units(n ./ sum(n))
end

function to_spec_compounds(st,mw,x::SP2) where {SP2<:MaterialCompounds}
    return to_spec_compounds(amount_type(st),st,mw,x)
end

function to_spec_compounds(amount_type::Tuple{SingleComponent,T},st,mw,x) where {T}
    return single_component_to_spec(st,mw,x)
end

function to_spec_compounds(amount_type::Tuple{T1,T2},st,mw,x::T1) where {T1<:MaterialCompounds,T2}
    sp = get_spec(T1,st)
    return normalize_units(value(sp))
end

function to_spec_compounds(amount_type::Tuple{T1,T2},st,mw,x) where {T1<:MaterialCompounds,T2}
    sp = get_spec(T1,st)
    return to_spec_mat(st,sp,mw,x)
end

function single_component_to_spec(st,mw,x::MaterialCompounds{MOLAR,FRACTION})
    return [_one(st)]
end

function single_component_to_spec(st,mw,x::MaterialCompounds{MASS,FRACTION})
    return [_one(st)]
end

function single_component_to_spec(st,mw,x::MaterialCompounds{MOLAR,TOTAL_AMOUNT})
    return [moles2(st,mw)]
end

function single_component_to_spec(st,mw,x::MaterialCompounds{MASS,TOTAL_AMOUNT})
    return [mass2(st,mw)]
end


text:
@recipe function f(val::VerySimpleRaster)
    seriestype --> :heatmap
    aspect_ratio --> 1
    grid --> false
#    framestyle --> :none
    val.xs, val.ys, @view replace(val.mat, val.nodata => NaN)'[end:-1:1, :]
end


text:
function process_all_retained!(gen_fn::Unfold{T,U}, params::Tuple, argdiffs::Tuple,
                               choices_or_selection, prev_length::Int, new_length::Int,
                               retained_and_targeted::Set{Int}, state) where {T,U}

    len_diff = argdiffs[1]
    init_state_diff = argdiffs[2]
    param_diffs = argdiffs[3:end] # a tuple of diffs

    if any(diff != NoChange() for diff in param_diffs)

        # visit every retained kernel application
        state_diff = init_state_diff
        for key=1:min(prev_length,new_length)
            state_diff = process_retained!(gen_fn, params, choices_or_selection,
                key, (NoChange(), state_diff, param_diffs...), state)
        end

    else
        # every parameter diff is NoChange()

        # visit only certain retained kernel applications
        to_visit::Vector{Int} = sort(collect(retained_and_targeted))
        key = 0
        state_diff = init_state_diff
        if state_diff != NoChange()
            key = 1
            visit = true
            while visit && key <= min(prev_length, new_length)
                state_diff = process_retained!(gen_fn, params, choices_or_selection,
                    key, (NoChange(), state_diff, param_diffs...), state)
                key += 1
                visit = (state_diff != NoChange())
            end
        end
        for i=1:length(to_visit)
            if key > to_visit[i]
                # we have already visited it
                continue
            end
            key = to_visit[i]
            visit = true
            while visit && key <= min(prev_length, new_length)
                state_diff = process_retained!(gen_fn, params, choices_or_selection,
                    key, (NoChange(), state_diff, param_diffs...), state)
                key += 1
                visit = (state_diff != NoChange())
            end
        end
    end
end

"""
Process all new applications.
"""
function process_all_new!(gen_fn::Unfold{T,U}, params::Tuple, choices_or_selection,
                          prev_len::Int, new_len::Int, state) where {T,U}
    for key=prev_len+1:new_len
        process_new!(gen_fn, params, choices_or_selection, key, state)
    end
end


text:
using Documenter, Metida
#using DocumenterLaTeX

makedocs(
    modules = [Metida],
    sitename = "Metida.jl",
    authors = "Vladimir Arnautov",
    pages = [
        "Home" => "index.md",
        "User guide" => [
        "First step" => "instanduse.md",
        "Examples" => "examples.md",
        "Bioequivalence" => "bioequivalence.md",
        "NLopt" => "nlopt.md",
        "CUDA" => "cuda.md",
        "F.A.Q." => "faq.md",
        ],
        "Details" => "details.md",
        "Benchmark" => "bench.md",
        "Validation" => "validation.md",
        "API" => "api.md",
        "Citation & Reference" => "ref.md",
    ],
)
#=
makedocs(format = DocumenterLaTeX.LaTeX(),
modules = [Metida],
sitename = "Metida.jl",
authors = "Vladimir Arnautov",
pages = [
    "Home" => "index.md",
    "User guide" => [
    "First step" => "instanduse.md",
    "Examples" => "examples.md",
    "Bioequivalence" => "bioequivalence.md",
    "NLopt" => "nlopt.md",
    "CUDA" => "cuda.md",
    "F.A.Q." => "faq.md",
    ],
    "Details" => "details.md",
    "Benchmark" => "bench.md",
    "Validation" => "validation.md",
    "API" => "api.md",
    "Citation & Reference" => "ref.md",
],
)
=#
deploydocs(repo = "github.com/PharmCat/Metida.jl.git", push_preview = true,
)


text:
using LinearAlgebraicRepresentation, ViewerGL, SparseArrays
Lar = LinearAlgebraicRepresentation; GL = ViewerGL

n,m,p = 1,1,1
V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
cube = V,FV,EV

threecubes = Lar.Struct([ cube,
    Lar.t(.3,.4,.25), Lar.r(pi/5,0,0), Lar.r(0,0,pi/12), cube,
    Lar.t(-.2,.4,-.2), Lar.r(0,pi/5,0), Lar.r(0,pi/12,0), cube ])

V,FV,EV = Lar.struct2lar(threecubes)
GL.VIEW([ GL.GLGrid(V,FV,GL.COLORS[1],.5), GL.GLFrame ])

cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.space_arrangement( W, cop_EV, cop_FE)

#triangulated_faces = Lar.triangulate(V, [copEV, copFE])
W = convert(Lar.Points, V')
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)

GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
meshes = GL.GLExplode(V,CVs[1:end],8,4,6,99,0.5);
GL.VIEW( push!( meshes, GL.GLFrame) );


text:
using DifferentialEquations
using SimpleDiffEq
using DataFrames
using StatsPlots
using BenchmarkTools

@inline function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

function sir_map!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    infection = rate_to_proportion(β*c*I/N,δt)*S
    recovery = rate_to_proportion(γ,δt)*I
    @inbounds begin
        du[1] = S-infection
        du[2] = I+infection-recovery
        du[3] = R+recovery
    end
    nothing
end;

δt = 0.1
nsteps = 400
tmax = nsteps*δt
tspan = (0.0,nsteps)
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0];

p = [0.05,10.0,0.25,δt]; # β,c,γ,δt

prob_map = DiscreteProblem(sir_map!,u0,tspan,p);

sol_map = solve(prob_map,solver=FunctionMap);

df_map = DataFrame(sol_map')
df_map[!,:t] = t;

@df df_map plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_map,solver=FunctionMap)


text:
using LinearAlgebra: pinv
using DSP: conv

export savitzky_golay_filter

"""
	savitzky_golay_filter(y::AbstractVector, window_size::Integer, polynomial_order::Integer; deriv_order::Integer = 0, boundary_mode = :interpolation)

Apply Savitzky-Golay polynomial smoothing to input data `y` using a polynomial of order `polynomial_order` fit to a moving window `window_size` points wide. Optionally derivatives can be taken by specifying the `deriv_order` and the caller is responsible for the appropriate scaling by the point spacing. Handling of data within half the window size is specified by `boundary_mode`. When set to `:interpolation` the polynomial fit will be used; when set to `:nearest` the data will be padded using the edge values before convolution and the valid portion will then be returned.

# References
1. Savitzky, A., & Golay, M. J. E. (1964). Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 36(8), 1627–1639. https://doi.org/10.1021/ac60214a047
2. Steinier, J., Termonia, Y., & Deltour, J. (1972). Comments on Smoothing and differentiation of data by simplified least square procedure. Analytical Chemistry, 44(11), 1906–1909. https://doi.org/10.1021/ac60319a045
3. Press, W. H., & Teukolsky, S. A. (1990). Savitzky-Golay Smoothing Filters. Computers in Physics, 4(6), 669. https://doi.org/10.1063/1.4822961
"""
function savitzky_golay_filter(y::AbstractVector, window_size::Integer, polynomial_order::Integer; deriv_order::Integer = 0, boundary_mode = :interpolation)

	# input validity checks
   	@assert isodd(window_size) "Window size must be an odd integer, i.e. fitting 2m + 1 points around the current value."
	@assert polynomial_order < window_size "Polynomial order must be less than the window size."
	@assert boundary_mode in (:interpolation, :nearest) "boundary_mode must be one of :interpolation, :nearest"

	# window size is 2m + 1 points
   	m = (window_size - 1) ÷ 2

	# build the Vandermonde design matrix A. Each row corresponds to a point in the fitting window -m:m
	# and each columns correspond to powers in the range 0:polynomial_order
   	fitting_points = -m:m
   	A = Matrix{Float64}(undef, window_size, polynomial_order + 1)
   	for i in 1:window_size, j in 1:polynomial_order + 1
        A[i,j] = fitting_points[i]^(j - 1)
    end

	if boundary_mode == :interpolation
		# for interpolation we'll want the full pseudo-inverse so we can calculate all the fit values at the edges
		# Ap = y
		C = pinv(A)

		# the filter coefficients are the rows of `C`
		filter_coeffs = C[deriv_order + 1,:] * factorial(deriv_order)

		# convolve with the filter coefficients with a couple extra steps:
		# 1. because of convolution will reverse coefficients we flip before
		# 2. c = conv(a,b) will return a vector of length(c) = length(a) + length(b) - 1 so we chop off the first and last m points
		smoothed = conv(reverse(filter_coeffs), y)[m+1:end-m]

		# for interpolation edge handling calculate the full fits
		if deriv_order == 0
			# if we are just smoothing then we can use the design and coefficient matrix as is
			AC = A*C
			smoothed[1:m] = (AC*y[1:window_size])[1:m]
			smoothed[end-m+1:end] = (AC*y[end-window_size+1:end])[end-m+1:end]
		else
			# otherwise we need to differentiate the polynomial coefficients
			# first m points
			p = C * y[1:window_size]
			for _ in 1:deriv_order
				p = [(i-1)*p[i] for i in 2:length(p)]
			end
			smoothed[1:m] = A[1:m, 1:size(A,2)-deriv_order]*p
			# last m points
			p = C * y[end-window_size+1:end]
			for _ in 1:deriv_order
				p = [(i-1)*p[i] for i in 2:length(p)]
			end
			smoothed[end-m+1:end] = A[m+2:end, 1:size(A,2)-deriv_order]*p

		end

		return smoothed

	elseif boundary_mode == :nearest

		# here we only need a single set of coefficients and so we can least-squares solve AᵀCᵀ = I for for a single row by picking a single column out of I
		Icol = zeros(Float64, polynomial_order+1, 1)
		Icol[deriv_order + 1] = 1.0
		filter_coeffs = transpose(A) \ Icol

		# pad the signal with the endpoints
		padded_y = [y[1] * ones(m); vec(y); y[end] * ones(m)]

		# convolve with filter
		smoothed = conv(filter_coeffs[end:-1:1], padded_y)

		# and return the valid midsection
		return smoothed[window_size:end-2*m]

	end

end