text:
# HistoryRecorder
# maintained by @zsunberg

"""
A simulator that records the history for later examination

The simulation will be terminated when either
1) a terminal state is reached (as determined by `isterminal()` or
2) the discount factor is as small as `eps` or
3) max_steps have been executed

Keyword Arguments:
    - `rng`: The random number generator for the simulation
    - `capture_exception::Bool`: whether to capture an exception and store it in the history, or let it go uncaught, potentially killing the script
    - `show_progress::Bool`: show a progress bar for the simulation
    - `eps`
    - `max_steps`
    - `sizehint::Int`: the expected length of the simulation (for preallocation)

Usage (optional arguments in brackets):

    hr = HistoryRecorder()
    history = simulate(hr, pomdp, policy, [updater [, init_belief [, init_state]]])
"""
mutable struct HistoryRecorder <: Simulator
    rng::AbstractRNG

    # options
    capture_exception::Bool
    show_progress::Bool

    # optional: if these are null, they will be ignored
    max_steps::Union{Nothing,Any}
    eps::Union{Nothing,Any}
    sizehint::Union{Nothing,Integer}
end

# This is the only stable constructor
function HistoryRecorder(;rng=MersenneTwister(rand(UInt32)),
                          eps=nothing,
                          max_steps=nothing,
                          sizehint=nothing,
                          capture_exception=false,
                          show_progress=false)
    return HistoryRecorder(rng, capture_exception, show_progress, max_steps, eps, sizehint)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy) begin
    @req updater(::typeof(policy))
    up = updater(policy)
    @subreq simulate(sim, pomdp, policy, up)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater) begin
    @req initialstate_distribution(::typeof(pomdp))
    dist = initialstate_distribution(pomdp)
    @subreq simulate(sim, pomdp, policy, bu, dist)
end

function simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater=updater(policy))
    dist = initialstate_distribution(pomdp)
    return simulate(sim, pomdp, policy, bu, dist)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater, dist::Any) begin
    P = typeof(pomdp)
    S = statetype(pomdp)
    A = actiontype(pomdp)
    O = obstype(pomdp)
    @req initialize_belief(::typeof(bu), ::typeof(dist))
    @req isterminal(::P, ::S)
    @req discount(::P)
    @req generate_sor(::P, ::S, ::A, ::typeof(sim.rng))
    b = initialize_belief(bu, dist)
    B = typeof(b)
    @req action(::typeof(policy), ::B)
    @req update(::typeof(bu), ::B, ::A, ::O)
end

function simulate(sim::HistoryRecorder,
                           pomdp::POMDP{S,A,O},
                           policy::Policy,
                           bu::Updater,
                           initialstate_dist::Any,
                           initialstate::Any=get_initialstate(sim, initialstate_dist)
                  ) where {S,A,O}

    initial_belief = initialize_belief(bu, initialstate_dist)
    if sim.max_steps == nothing
        max_steps = typemax(Int)
    else
        max_steps = sim.max_steps
    end
    if sim.eps != nothing
        max_steps = min(max_steps, ceil(Int,log(sim.eps)/log(discount(pomdp))))
    end
    if sim.sizehint == nothing
        sizehint = min(max_steps, 1000)
    else
        sizehint = sim.sizehint
    end

    # aliases for the histories to make the code more concise
    sh = sizehint!(Vector{S}(undef, 0), sizehint)
    ah = sizehint!(Vector{A}(undef, 0), sizehint)
    oh = sizehint!(Vector{O}(undef, 0), sizehint)
    bh = sizehint!(Vector{typeof(initial_belief)}(undef, 0), sizehint)
    rh = sizehint!(Vector{Float64}(undef, 0), sizehint)
    ih = sizehint!(Vector{Any}(undef, 0), sizehint)
    aih = sizehint!(Vector{Any}(undef, 0), sizehint)
    uih = sizehint!(Vector{Any}(undef, 0), sizehint)
    exception = nothing
    backtrace = nothing

    push!(sh, initialstate)
    push!(bh, initial_belief)

    if sim.show_progress
        if (sim.max_steps == nothing) && (sim.eps == nothing)
            error("If show_progress=true in a HistoryRecorder, you must also specify max_steps or eps.")
        end
        prog = Progress(max_steps, "Simulating..." )
    end

    disc = 1.0
    step = 1

    try
        while !isterminal(pomdp, sh[step], get(ah,step-1,nothing), get(oh,step-1,nothing)) && step <= max_steps
            a, ai = action_info(policy, bh[step])
            push!(ah, a)
            push!(aih, ai)

            sp, o, r, i = generate_sori(pomdp, sh[step], ah[step], sim.rng)

            push!(sh, sp)
            push!(oh, o)
            push!(rh, r)
            push!(ih, i)

            bp, ui = update_info(bu, bh[step], ah[step], oh[step])
            bh = push_belief(bh, bp)
            push!(uih, ui)

            step += 1

            if sim.show_progress
                next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            exception = ex
            backtrace = catch_backtrace()
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return POMDPHistory(sh, ah, oh, bh, rh, ih, aih, uih, discount(pomdp), exception, backtrace)
end

@POMDP_require simulate(sim::HistoryRecorder, mdp::MDP, policy::Policy) begin
    init_state = initialstate(mdp, sim.rng)
    @subreq simulate(sim, mdp, policy, init_state)
end

@POMDP_require simulate(sim::HistoryRecorder, mdp::MDP, policy::Policy, initialstate::Any) begin
    P = typeof(mdp)
    S = statetype(mdp)
    A = actiontype(mdp)
    @req isterminal(::P, ::S)
    @req action(::typeof(policy), ::S)
    @req generate_sr(::P, ::S, ::A, ::typeof(sim.rng))
    @req discount(::P)
end

function simulate(sim::HistoryRecorder,
                  mdp::MDP{S,A}, policy::Policy,
                  init_state::S=get_initialstate(sim, mdp)) where {S,A}

    if sim.max_steps == nothing
        max_steps = typemax(Int)
    else
        max_steps = sim.max_steps
    end
    if sim.eps != nothing
        max_steps = min(max_steps, ceil(Int,log(sim.eps)/log(discount(mdp))))
    end
    if sim.sizehint == nothing
        sizehint = min(max_steps, 1000)
    else
        sizehint = sim.sizehint
    end

    # aliases for the histories to make the code more concise
    sh = sizehint!(Vector{S}(undef, 0), sizehint)
    ah = sizehint!(Vector{A}(undef, 0), sizehint)
    rh = sizehint!(Vector{Float64}(undef, 0), sizehint)
    ih = sizehint!(Vector{Any}(undef, 0), sizehint)
    aih = sizehint!(Vector{Any}(undef, 0), sizehint)
    exception = nothing
    backtrace = nothing

    if sim.show_progress
        prog = Progress(max_steps, "Simulating..." )
    end

    push!(sh, init_state)

    disc = 1.0
    step = 1

    try
        while !isterminal(mdp, sh[step], get(ah,step-1,nothing)) && step <= max_steps
            a, ai = action_info(policy, sh[step])
            push!(ah, a)
            push!(aih, ai)

            sp, r, i = generate_sri(mdp, sh[step], ah[step], sim.rng)

            push!(sh, sp)
            push!(rh, r)
            push!(ih, i)

            disc *= discount(mdp)
            step += 1

            if sim.show_progress
                next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            exception = ex
            backtrace = catch_backtrace()
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return MDPHistory(sh, ah, rh, ih, aih, discount(mdp), exception, backtrace)
end

function get_initialstate(sim::Simulator, initialstate_dist)
    return rand(sim.rng, initialstate_dist)
end

function get_initialstate(sim::Simulator, mdp::Union{MDP,POMDP})
    return initialstate(mdp, sim.rng)
end

# this is kind of a hack in cases where the belief isn't stable
push_belief(bh::Vector{T}, b::T) where T = push!(bh, b)
function push_belief(bh::Vector{T}, b::B) where {B, T}
    if !(T isa Union) # if T is not already a Union, try making a Union of the two types; don't jump straight to Any
        new = Vector{Union{T,B}}(undef, length(bh)+1)
    else
        new = Vector{promote_type(T, B)}(undef, length(bh)+1)
    end
    new[1:end-1] = bh
    new[end] = b
    return new
end


text:
module test_h3_api_hierarchy

using Test
using H3.API # h3ToParent h3ToChildren maxH3ToChildrenSize compact uncompact maxUncompactSize
using H3.Lib

res = 1
sf = GeoCoord(0.659966917655, 2pi - 2.1364398519396)
child = geoToH3(sf, res)
@test child == 0x081283ffffffffff
@test h3GetResolution(child) == res

parent = h3ToParent(child, res)
@test parent == child

sfHex8 = geoToH3(sf, 8)
sfHex9s = h3ToChildren(sfHex8, 9)
center = h3ToGeo(sfHex8)
sfHex9_0 = geoToH3(center, 9)
@test sfHex9s[1] == sfHex9_0

@test maxH3ToChildrenSize(sfHex8, 9) == 7

hexCount = Lib.NUM_BASE_CELLS
res0Hexes = fill(0x0000_0000_0000_0001, hexCount)
compressed = compact(res0Hexes)
@test length(compressed) == hexCount

maxHexes = maxUncompactSize(compressed, 0)
@test maxHexes == hexCount

decompressed = uncompact(compressed, 0)
@test length(decompressed) == 122

decompressed = uncompact(compressed, 1)
@test length(decompressed) == 854

end # module test_h3_api_hierarchy

ext: jl


text:
using Documenter
using OpenCV

makedocs(
    modules = [OpenCV],
    clean   = false,
    format   = Documenter.Formats.HTML,
    sitename = "OpenCV.jl",
    pages = Any[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Getting started" => "getting-started.md",
        "Modules" => Any[
            "Core" => "modules/core.md"
            "ImgProc" => "modules/imgproc.md"
            "ImgCodecs" => "modules/imgcodecs.md"
            "VideoIO" => "modules/videoio.md"
            "HighGUI" => "modules/highgui.md"
            "Calib3d" => "modules/calib3d.md"
            "LibOpenCV" => "modules/libopencv.md"
            ]
        ],
)


text:
using Documenter
using TFT

makedocs(
    modules     = [TFT],
    format      = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename    = "TFT.jl",
    authors     = "Tom Van Acker",
    pages       = [ "Home"              => "index.md",
                    "DSL Manual"        =>
                         ["DTFT"                       => "dtft.md",
                          "Utilities"                  => "util.md"],
                    "Developer"         =>
                         ["Developer's Guide"          => "guide.md",
                          "Mathematical Background"    => "math.md"]
                  ]
)

deploydocs(
     repo = "github.com/timmyfaraday/TFT.jl.git"
)


text:
function likelihood(k, model::AbstractBinomialModel, observation)
    return mean(
                exp.(-0.5f0 .* ((observation .- model.q .* k) ./ model.σ).^2)
                ./ (sqrt(2*Float32(pi)) .* model.σ)
           , dims = 2
           )[:,1]
end

function kernel_likelihood_indices!(u, v, idxT, kT, q, σ, observation, r, randstates)
    # use column-index first (gives a 3x speedup)
    # kT stands for the transpose of k
    # idxT stands for the transpose of idx
    j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    M_in, M_out = size(kT)
    @inbounds if j <= M_out
        vsum = 0f0
        CurMax = 1f0
        for i in 1:M_in
            # omitting normalization constant here; it is only needed for u
            vi      = CUDA.exp(-0.5f0 *((observation - q[j] * kT[i,j]) / σ[j])^2)
            vsum   += vi
            v[i, j] = vsum
            # sample descending sequence of sorted random numbers
            # r[i,M_in] >= ... >= r[i,2] >= r[i,1]
            # Algorithm by:
            # Bentley & Saxe, ACM Transactions on Mathematical Software, Vol 6, No 3
            # September 1980, Pages 359--364
            mirrori = M_in - i + 1 # mirrored index i
            CurMax *= CUDA.exp(CUDA.log(GPUArrays.gpu_rand(Float32, CUDA.CuKernelContext(), randstates)) / mirrori)
            r[mirrori, j] = CurMax
        end
        # compute average likelihood across inner particles
        # (with normalization constant that was omitted from v for speed)
        u[j] = vsum / (M_in * CUDA.sqrt(2*Float32(pi)) * σ[j])
        # O(n) binning algorithm for sorted samples
        bindex = 1 # bin index
        @inbounds for i in 1:M_in
            # scale random numbers (this is equivalent to normalizing v)
            rsample = r[i, j] * vsum
            # checking bindex <= M_in - 1 is redundant since
            # v[M_in, j] = vsum
            while rsample > v[bindex, j]
                bindex += 1
            end
            idxT[i, j] = bindex
        end
    end
    return nothing
end

function likelihood_indices(k, model::AbstractBinomialModel, observation)
    M_out, M_in = size(k)
    r           = CuArray{Float32}(undef, M_in, M_out)
    u           = CuArray{Float32}(undef, M_out)
    v           = CuArray{Float32}(undef, M_in, M_out)
    idx         = CuArray{Int}(undef, M_out, M_in)

    rng = GPUArrays.default_rng(CuArray)

    kernel  = @cuda launch=false kernel_likelihood_indices!(
                u, v,
                idx', k',
                model.q, model.σ,
                Float32(observation),
                r,
                rng.state
              )
    config  = launch_configuration(kernel.fun)
    threads = Base.min(M_out, config.threads, 256)
    blocks  = cld(M_out, threads)
    kernel(
        u, v,
        idx', k', # pass transposes for column-major indexing (gives a 3x speedup)
        model.q, model.σ,
        observation,
        r,
        rng.state
        ;
        threads=threads, blocks=blocks
    )
    return u, idx
end

function inner_resample_helper!(in, out, idx)
    function kernel(in, out, idx)
        j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
        @inbounds if j <= size(out, 2)
            @inbounds for i in 1:size(out, 1)
                out[i, j] = in[i, idx[i, j]]
            end
        end#if
        return nothing
    end
    kernel  = @cuda launch=false kernel(in, out, idx)
    config  = launch_configuration(kernel.fun)
    threads = Base.min(size(out, 2), config.threads, 256)
    blocks  = cld(size(out, 2), threads)
    kernel(in, out, idx; threads=threads, blocks=blocks)
    return out
end

function inner_resample_helper!(in, idx)
    out = similar(in)
    inner_resample_helper!(in, out, idx)
    in .= out
    return in
end

function likelihood_resample!(state::BinomialState, model, observation::BinomialObservation)
    u, idx = likelihood_indices(state.k, model, observation.EPSP)
    M_out  = size(state.n, 1)
    inner_resample_helper!(state.n, idx)
    inner_resample_helper!(state.k, idx)
    return u
end


text:
using MOONS.Grids
using MOONS.Fields
using Plots

abstract type AbstractGridPlot end
struct Isosurfaces3D <: AbstractGridPlot end
struct Streamlines3D <: AbstractGridPlot end
struct GridVis{slice} <: AbstractGridPlot end
struct Contour2D{slice} <: AbstractGridPlot end
struct Profile1D{axis} <: AbstractGridPlot end

plot_grid(g::Grid, filename::AbstractString, grid_viz::GridVis{slice}) where {F, slice} = nothing
plot_func(f::F, g::Grid, filename::AbstractString, ::GridVis{slice}) where {F, slice} = nothing

function plot_grid(g::Grid, filename::AbstractString, grid_viz::GridVis{slice}) where {F, slice}
  return plot_func((x,y) -> 1, g, filename, grid_viz, Plots.wireframe)
end

function plot_func(f::F,
                   g::Grid,
                   filename::AbstractString,
                   ::GridVis{slice},
                   viz_type::G=contourf) where {G, F, slice}
    i_orth = orthogs(Val(slice))
    _x = g.c[i_orth[1]]
    _y = g.c[i_orth[2]]
    x = _x.n
    y = _y.n
    data = [f(x.h[i], y.h[i]) for i in over_points(_x, Interior(), Vertex1D()),
                                  j in over_points(_y, Interior(), Vertex1D())]
    X = [x.h[i] for i in over_points(_x, Interior(), Vertex1D())]
    Y = [y.h[j] for j in over_points(_y, Interior(), Vertex1D())]
    viz_type(X,Y,data, camera=(90, 90), size=(2*320,2*300), aspect_ratio=:equal)
    viz_type(X,Y,data)
    savefig(filename)
end

# import Plots: plot

# function plot(grid::Grid, f::AbstractField)
#     plot()
# end


text:
"""
# Stipple.Elements

The `Elements` module provides utility methods for interfacing between Julia and Vue.js.
"""
module Elements

import Genie
using Stipple

import Genie.Renderer.Html: HTMLString, normal_element

export root, elem, vm, @iif, @elsiif, @els, @recur, @text, @bind, @data, @on, @showif
export stylesheet

#===#

"""
    `function root(app::M)::String where {M<:ReactiveModel}`

Generates a valid JavaScript object name to be used as the name of the Vue app -- and its respective HTML container.
"""
function root(app::M)::String where {M<:ReactiveModel}
  Genie.Generator.validname(typeof(app) |> string)
end

function root(app::Type{M})::String where {M<:ReactiveModel}
  Genie.Generator.validname(app |> string)
end

"""
    `function elem(app::M)::String where {M<:ReactiveModel}`

Generates a JS id `#` reference to the DOM element containing the Vue app template.
"""
function elem(app::M)::String where {M<:ReactiveModel}
  "#$(root(app))"
end

const vm = root

#===#

"""
    `function vue_integration(model::M; vue_app_name::String, endpoint::String, channel::String, debounce::Int)::String where {M<:ReactiveModel}`

Generates the JS/Vue.js code which handles the 2-way data sync between Julia and JavaScript/Vue.js.
It is called internally by `Stipple.init` which allows for the configuration of all the parameters.
"""
function vue_integration(m::Type{M}; vue_app_name::String = "StippleApp", core_theme::Bool = true,
                          channel::String = Genie.config.webchannels_default_route,
                          debounce::Int = Stipple.JS_DEBOUNCE_TIME,
                          transport::Module = Genie.WebChannels)::String where {M<:ReactiveModel}
  model = Base.invokelatest(m)

  vue_app = replace(json(model |> Stipple.render), "\"{" => " {")
  vue_app = replace(vue_app, "}\"" => "} ")
  vue_app = replace(vue_app, "\"$channel\"" => "CHANNEL")

  output =
  string(
    "
    Stipple.init($( core_theme ? "{theme: 'stipple-blue'}" : "" ));
    var $vue_app_name = new Vue($( replace(vue_app, "'$(Stipple.UNDEFINED_PLACEHOLDER)'"=>Stipple.UNDEFINED_VALUE) ));
    "

    ,

    join(
      [Stipple.watch(vue_app_name, field, "CHANNEL", debounce, model) for field in fieldnames(m) if Stipple.ispublic(field, model)]
    )

    ,

    """

  window.parse_payload = function(payload){
    if (payload.key) {
      window.$(vue_app_name).revive_payload(payload)
      window.$(vue_app_name).updateField(payload.key, payload.value);
    }
  }

  function app_ready() {
    if (
      (document.readyState === 'complete' || document.readyState === 'interactive')
        $(transport == Genie.WebChannels ? " && (Genie.WebChannels.socket.readyState === 1)" : "")
      ) {

      if (Genie.Settings.env == 'dev') {
        console.info('Waiting ' + $vue_app_name.isreadydelay + 'ms');
      }

      setTimeout(function(){
        $vue_app_name.isready = true;
        if (Genie.Settings.env == 'dev') {
          console.info('App ready');
        }
      }, $vue_app_name.isreadydelay); // let's give it a bit to process server side events

      $(transport == Genie.WebChannels &&
      "
      try {
        if (Genie.Settings.webchannels_keepalive_frequency > 0) {
          setInterval(keepalive, Genie.Settings.webchannels_keepalive_frequency);
        }
      } catch (e) {
        if (Genie.Settings.env == 'dev') {
          console.error('Error setting WebSocket keepalive interval: ' + e);
        }
      }
      ")

    } else {
      if (Genie.Settings.env == 'dev') {
        console.info('App starting');
      }
      setTimeout(app_ready, Genie.Settings.webchannels_timeout);
    }
  };

  window.onload = function() {
    if (Genie.Settings.env == 'dev') {
      console.info('Loading completed');
    }
    app_ready();
  }
  """
  )

  output = repr(output)
  output[2:prevind(output, lastindex(output))]
end

#===#

"""
    `@iif(expr)`

Generates `v-if` Vue.js code using `expr` as the condition.
<https://vuejs.org/v2/api/#v-if>

### Example

```julia
julia> span("Bad stuff's about to happen", class="warning", @iif(:warning))
"<span class=\"warning\" v-if='warning'>Bad stuff's about to happen</span>"
```
"""
macro iif(expr)
  :( "v-if='$($(esc(expr)))'" )
end

"""
    `@elsiif(expr)`

Generates `v-else-if` Vue.js code using `expr` as the condition.
<https://vuejs.org/v2/api/#v-else-if>

### Example

```julia
julia> span("An error has occurred", class="error", @elsiif(:error))
"<span class=\"error\" v-else-if='error'>An error has occurred</span>"
```
"""
macro elsiif(expr)
  :( "v-else-if='$($(esc(expr)))'" )
end

"""
    `@els(expr)`

Generates `v-else` Vue.js code using `expr` as the condition.
<https://vuejs.org/v2/api/#v-else>

### Example

```julia
julia> span("Might want to keep an eye on this", class="notice", @els(:notice))
"<span class=\"notice\" v-else='notice'>Might want to keep an eye on this</span>"
```
"""
macro els(expr)
  :( "v-else='$($(esc(expr)))'" )
end

"""
Generates `v-for` directive to render a list of items based on an array.
<https://vuejs.org/v2/guide/list.html#Mapping-an-Array-to-Elements-with-v-for>

### Example

```julia
julia> p(" {{todo}} ", class="warning", @recur(:"todo in todos"))
"<p v-for='todo in todos'>\n {{todo}} \n</p>\n"
```

"""
macro recur(expr)
  :( "v-for='$($(esc(expr)))'" )
end

"""
    `@text(expr)`

Creates a `v-text` or a `text-content.prop` Vue biding to the element's `textContent` property.
<https://vuejs.org/v2/api/#v-text>

### Example

```julia
julia> span("", @text("abc | def"))
"<span :text-content.prop='abc | def'></span>"

julia> span("", @text("abc"))
"<span v-text='abc'></span>"
```
"""
macro text(expr)
  quote
    directive = occursin(" | ", string($(esc(expr)))) ? ":text-content.prop" : "v-text"
    "$(directive)='$($(esc(expr)))'"
  end
end

"""
    `@bind(expr, [type])`

Binds a model parameter to a Vue component, generating a `v-model` property, optionally defining the parameter type.
<https://vuejs.org/v2/api/#v-model>

### Example

```julia
julia> input("", placeholder="Type your name", @bind(:name))
"<input placeholder=\"Type your name\"  v-model='name' />"

julia> input("", placeholder="Type your name", @bind(:name, :identity))
"<input placeholder=\"Type your name\"  v-model.identity='name' />"
```
"""
macro bind(expr)
  :( "v-model='$($(esc(expr)))'" )
end

macro bind(expr, type)
  :( "v-model.$($(esc(type)))='$($(esc(expr)))'" )
end

"""
    `@data(expr)`

Creates a Vue.js data binding for the elements that expect it.

### Example

```julia
julia> plot(@data(:piechart), options! = "plot_options")
"<template><apexchart :options=\"plot_options\" :series=\"piechart\"></apexchart></template>"
```
"""
macro data(expr)
  quote
    x = $(esc(expr))
    if typeof(x) <: Union{AbstractString,Symbol}
      Symbol(x)
    else
      strx = strip("$x")
      startswith(strx, "Any[") && (strx = strx[4:end])

      JSONText(string(":", strx))
    end
  end
end

"""
    `on(action, expr)`

Defines a js routine that is called by the given `action` of the Vue component, e.g. `:click`, `:input`

### Example

```julia
julia> input("", @bind(:input), @on("keyup.enter", "process = true"))
"<input  v-model='input' v-on:keyup.enter='process = true' />"
```
"""
macro on(args, expr)
  :( "v-on:$(string($(esc(args))))='$(replace($(esc(expr)),"'" => raw"\'"))'" )
end

"""
    `@showif(expr, [type])`

v-show will always be rendered and remain in the DOM; v-show only toggles the display CSS property of the element.
<https://vuejs.org/v2/guide/conditional.html#v-show>

Difference between @showif and @iif when to use either

v-if has higher toggle costs while v-show has higher initial render costs

### Example

```julia
julia> h1("Hello!", @showif(:ok))
"<h1 v-show="ok">Hello!</h1>"
```
"""
macro showif(expr)
  :( "v-show='$($(esc(expr)))'" )
end

#===#

"""
    `function stylesheet(href::String; args...) :: String`

Generates the corresponding HTML `link` tag to reference the CSS stylesheet at `href`.

### Example

```julia
julia> stylesheet("https://fonts.googleapis.com/css?family=Material+Icons")
"<link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" />"
```
"""
function stylesheet(href::String; args...) :: ParsedHTMLString
  Genie.Renderer.Html.link(href=href, rel="stylesheet", args...)
end

end


text:
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

library_path = joinpath("..", "src")
(library_path in LOAD_PATH) || push!(LOAD_PATH, library_path)

using Documenter, Stripeline

DocMeta.setdocmeta!(Stripeline, :DocTestSetup, :(using Stripeline); recursive=true)

ENV["GKSwstype"] = "100"  # Disable display of plots

makedocs(
    modules = [Stripeline],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    sitename = "Stripeline.jl",
    pages = [
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "Basic functions" => "basic.md",
        "Instrument database" => "instrumentdb.md",
        "Time-ordered data" => "tod.md",
        "Scanning strategy" => "scanning.md",
        "Data acquisition" => "acquisition.md",
        "Map making" => "mapmaking.md",
	"Full-scale simulation tutorial" => "simulation_tutorial.md",
    ],
    debug = true,
)

deploydocs(
    repo = "github.com/lspestrip/Stripeline.jl.git",
)


text:
using Gadfly
using Compose
using DataFrames
using StatsBase

abstract type Gridworld end

mutable struct simpleGrid <: Gridworld
    width::Number
    endstates::Any
    discount::Number
end

mutable struct volcanoGrid <: Gridworld
    width::Number
    penaltyStates::Any
    penalty::Number
    rewardStates::Any
    reward::Number
    probSlip::Number
    discount::Number
end

"Check if current state is in End State"
function isEnd(mdp::simpleGrid, state)
    return state in mdp.endstates
end

"Check if current state is in End State"
function isEnd(mdp::volcanoGrid, state)
    if state ∈ mdp.penaltyStates
        return true
    elseif state ∈ mdp.rewardStates
        return true
    else
        return false
    end
end

"Retunrs all possible actions given current state"
function actions(mdp::Gridworld, state)
    actions = ["left", "right", "up", "down"]
    if isEnd(mdp, state)
        return []
    end
    if state in 1:mdp.width
        filter!(x->x != "up", actions)
    end
    if state in mdp.width^2-mdp.width + 1:mdp.width^2
        filter!(x->x != "down", actions)
    end
    if mod(state, mdp.width) == 0
        filter!(x->x != "right", actions)
    end
    if mod(state, mdp.width) == 1
        filter!(x->x !="left", actions)
    end
    return actions
end

function newState(mdp::Gridworld, state::Number, action::String)
    if action == "right"
        return state + 1
    elseif action == "left"
        return state - 1
    elseif action == "down"
        return state + mdp.width
    elseif action == "up"
        return state - mdp.width
    end
end

"Helper function to return all states"
function states(mdp::Gridworld)
    return 1:mdp.width^2
end

"Returns tuple of (New State, Probability, Reward)"
function succProbReward(mdp::simpleGrid, state, action)
    result = []
    newS = newState(mdp, state, action)
    reward = isEnd(mdp, newS) ? 0 : -1
    push!(result, [newS, 1, reward])
    return result
end

"Returns tuple of (New State, Probability, Reward)"
function succProbReward(mdp::volcanoGrid, state, action)
    result = []
    n = length(actions(mdp, state)) - 1
    for i in actions(mdp, state)
        if newState(mdp, state, i) ∈ mdp.penaltyStates
            reward = mdp.penalty
        elseif newState(mdp, state, i) ∈ mdp.rewardStates
            reward = mdp.reward
        else
            reward = -1
        end
        if i == action
            push!(result, [newState(mdp, state, i), 1 - mdp.probSlip, reward])
        else
            push!(result, [newState(mdp, state, i), mdp.probSlip / n, reward])
        end
    end
    return result
end

"Sample from actions according to probabilities and return new state, reward"
function sampleAction(mdp::Gridworld, state, action)
    newStates = []
    probs = []
    rewards = []
    for (newState, prob, reward) in succProbReward(mdp, state, action)
        push!(newStates, newState)
        push!(probs, prob)
        push!(rewards, reward)
    end
    idx = sample(1:length(probs), ProbabilityWeights(convert(Array{Float64,1}, probs)))
    return newStates[idx], rewards[idx]
end

function TDIteraation(mdp::Gridworld, episodes, exploration_prob=0.2, startState=1, α=.001)
    value =  zeros(length(states(mdp)))
    delta = []
    function Q(mdp::Gridworld, state, action)
        q = sum(
                prob * (reward + mdp.discount * value[Int(newState)])
                for (newState, prob, reward) in succProbReward(mdp, state, action)
            )
        return q
    end
    function getAction(mdp::Gridworld, state)
        if rand() < exploration_prob
            return rand(actions(mdp, state))
        else
            return maximum((Q(mdp, state, action), action) for action in actions(mdp, state))[2]
        end
    end
    # For each trial
    for i ∈ 1:episodes
        state = startState
        # Run each episode
        while true
            action = getAction(mdp, state)
            newState, reward = sampleAction(mdp, state, action)
            before = value[Int(state)]
            value[Int(state)] += α .* (reward .+ mdp.discount .* value[Int(newState)] .- value[Int(state)])
            push!(delta, before - value[Int(state)])
            if isEnd(mdp, newState) break end
            state = newState
        end
    end
    #plot(1:length(delta), delta, Geom.line)
    return DataFrame(value=value, state = states(mdp))
end

"Impliments Value iteration on simpleGrid Struct"
function valueIteration(mdp::Gridworld)
    # Initialize value and policy arrays
    value =  Array{Float64}(undef, length(states(mdp)))
    p = Array{String}(undef, length(states(mdp)))

    for (i, state) in enumerate(states(mdp))
        value[i] = 0
    end
    # Q Function
    function Q(mdp::Gridworld, state, action)
        q = sum(
                prob * (reward + mdp.discount * value[Int(newState)])
                for (newState, prob, reward) in succProbReward(mdp, state, action)
            )
        return q
    end
    # This should be until convergence
    for j in 1:100
        vNew =  Array{Float64}(undef, length(states(mdp)))
        for (i, state) in enumerate(states(mdp))
            if isEnd(mdp, state)
                vNew[i] = 0
            else
                vNew[i] = maximum([Q(mdp, state, action) for action in actions(mdp, state)])
            end
        end
        value = vNew

        # Get policy
        for (i, state) in enumerate(states(mdp))
            if isEnd(mdp, state)
                p[i] = "None"
            else
                p[i] = maximum((Q(mdp, state, action), action) for action in actions(mdp, state))[2]
            end
        end
    end
    out = []
    for (i, state) in enumerate(states(mdp))
        push!(out, (state, value[i], p[i]))
    end
    return out
end

"Generates plot with V*(s) values for each state in Grid"
function plotValueGrid(mdp::simpleGrid)
    result = valueIteration(mdp)
    df = DataFrame(
        [(y=y, x=x,)
            for y in -(1:mdp.width)
            for x in 1:mdp.width
        ])
    df[!, :label] = [string(result[2]) for result in result]
    df[!, :endState] = [isEnd(mdp, state) for state in states(mdp)]
    end_coords = filter(row -> row.endState == true, df)
    x = mdp.width + .5
    plot(df, x=:x, y=:y, label=:label, Geom.label,
        Guide.xticks(ticks=0.5:x, label=false),
        Guide.yticks(ticks=-(0.5:x), label=false),
        Guide.title("Simple Gridworld V(s) ∀ s ∈ S"),
        Guide.xlabel(nothing), Guide.ylabel(nothing),
        Guide.Annotation(compose(context(),
        (context(),
            rectangle(
                end_coords[!,:x] .- .5,
                end_coords[!,:y] .- .5,
                ones(size(end_coords)[1]),
                ones(size(end_coords)[1])),
            fill("green"),fillopacity(0.1), stroke("green")),
        )),
        Theme(alphas = [.5], grid_color="black",
            grid_line_width=2pt, grid_line_style=:solid))
end

function plotValueGrid(mdp::volcanoGrid)
    result = valueIteration(mdp)
    df = DataFrame(
        [(y=y, x=x,)
            for y in -(1:mdp.width)
            for x in 1:mdp.width
        ])
    df[!, :label] = [string(round(result[2], digits=2)) for result in result]
    df[!, :rewardStates] = [state ∈ mdp.rewardStates for state in states(mdp)]
    df[!, :penaltyStates] = [state ∈ mdp.penaltyStates for state in states(mdp)]
    reward_coords = filter(row -> row.rewardStates == true, df)
    penalty_coords = filter(row -> row.penaltyStates == true, df)
    x = mdp.width + .5
    plot(df, x=:x, y=:y, label=:label, Geom.label,
        Guide.xticks(ticks=0.5:x, label=false),
        Guide.yticks(ticks=-(0.5:x), label=false),
        Guide.title("Volcano Gridworld V(s) ∀ s ∈ S"),
        Guide.xlabel(nothing), Guide.ylabel(nothing),
        Guide.Annotation(compose(context(),
        (context(),
            rectangle(
                reward_coords[!,:x] .- .5,
                reward_coords[!,:y] .- .5,
                ones(size(reward_coords)[1]),
                ones(size(reward_coords)[1])),
            fill("green"),fillopacity(0.1), stroke("green")),
        (context(),
            rectangle(
                penalty_coords[!,:x] .- .5,
                penalty_coords[!,:y] .- .5,
                ones(size(penalty_coords)[1]),
                ones(size(penalty_coords)[1])),
            fill("red"),fillopacity(0.1), stroke("red")),
        )),
        Theme(alphas = [.5], grid_color="black",
            grid_line_width=2pt, grid_line_style=:solid))
end

# grid = simpleGrid(5, [1, 25], 1)
# grid2 = volcanoGrid(4, [3,7], -100, [4], 10, 0.3, 1)
# plotValueGrid(grid)
# plotValueGrid(grid2)

max_stars_repo_name: blazickjp/ReinforcementLearning

max_issues_repo_name: blazickjp/ReinforcementLearning

max_forks_repo_name: blazickjp/ReinforcementLearning


text:
# it take Z as Gauss DCA
# and returns
# - fi_tens (q,N), fij_tens(q,q,N,N)
# - fi (q*N),fij (q*N*q*N)
# - Meff (effective numb of seq)
# - W (the vector of M weights)

# [normal gaussDCA returns N*(q-1), frequencies]

q = 21

function compute_new_fi_fij(Z::Matrix{Int8}, q, theta)

	#G.C
	#q+1 to get sequences which are 21X21

    W, Meff = GaussDCA.compute_weights(Z, q+1, theta)
    Pi_matr, Pij_matr = compute_all_freqs(Z, W, Meff,q+1)
    Pi_tens, Pij_tens = freq_tens_form(Pi_matr, Pij_matr, Z)

    return Pi_tens, Pij_tens, Pi_matr, Pij_matr, Meff, W
end

function compute_all_freqs(Z::Matrix{Int8}, W::Vector{Float64}, Meff::Float64,q::Int)
    N, M = size(Z)
    s = q - 1

    Ns = N * s

    Pij = zeros(Ns, Ns)
    Pi = zeros(Ns)

    ZZ = Vector{Int8}[vec(Z[i,:]) for i = 1:N]

    i0 = 0
    for i = 1:N
        Zi = ZZ[i]
        for k = 1:M
            a = Zi[k]
            a == q && continue
            Pi[i0 + a] += W[k]
        end
        i0 += s
    end
    Pi /= Meff

    i0 = 0
    for i = 1:N
        Zi = ZZ[i]
        j0 = i0
        for j = i:N
            Zj = ZZ[j]
            for k = 1:M
                a = Zi[k]
                b = Zj[k]
                (a == q || b == q) && continue
                Pij[i0+a, j0+b] += W[k]
            end
            j0 += s
        end
        i0 += s
    end
    for i = 1:Ns
        Pij[i,i] /= Meff
        for j = i+1:Ns
            Pij[i,j] /= Meff
            Pij[j,i] = Pij[i,j]
        end
    end

    return Pi, Pij
end

function freq_tens_form(fi, fij, Z::Matrix{Int8})

    N, M = size(Z)

	#adapt the frequencies
	fi_tens = zeros(q,N)

	for i = 1:N
		fi_tens[:,i] = fi[(1+(q)*(i-1)):((q)+(q)*(i-1))]
	end

	fij_tens= zeros(q,q,N,N);
	for i = 1:N
		for j = 1:N
			fij_tens[:,:,i,j] = fij[(1+(q)*(i-1)):((q)+(q)*(i-1)),(1+(q)*(j-1)):((q)+(q)*(j-1))]
		end
	end

	return fi_tens, fij_tens
end

max_stars_repo_name: juan-rodriguez-rivas/covmut

max_issues_repo_name: juan-rodriguez-rivas/covmut

max_forks_repo_name: juan-rodriguez-rivas/covmut


text:
@testset "pointer" begin

@testset "unsafe_load & unsafe_store!" begin

@eval struct LoadableStruct
    a::Int64
    b::UInt8
end
Base.one(::Type{LoadableStruct}) = LoadableStruct(1,1)
Base.zero(::Type{LoadableStruct}) = LoadableStruct(0,0)

@testset for T in (Int8, UInt16, Int32, UInt32, Int64, UInt64, Int128,
                   Float32, Float64,
                   LoadableStruct),
             cached in (false, true)
    d_a = Mem.upload(ones(T))
    d_b = Mem.upload(zeros(T))

    ptr_a = CUDAnative.DevicePtr{T,AS.Global}(Base.unsafe_convert(CuPtr{T}, d_a))
    ptr_b = CUDAnative.DevicePtr{T,AS.Global}(Base.unsafe_convert(CuPtr{T}, d_b))
    @test Mem.download(T, d_a) != Mem.download(T, d_b)

    let ptr_a=ptr_a, ptr_b=ptr_b #JuliaLang/julia#15276
        if cached && capability(dev) >= v"3.2"
            @on_device unsafe_store!(ptr_b, unsafe_cached_load(ptr_a))
        else
            @on_device unsafe_store!(ptr_b, unsafe_load(ptr_a))
        end
    end
    @test Mem.download(T, d_a) == Mem.download(T, d_b)
end

@testset "indexing" begin
    function kernel(src, dst)
        unsafe_store!(dst, CUDAnative.unsafe_cached_load(src, 4))
        return
    end

    T = Complex{Int8}

    src = Mem.upload([T(1) T(9); T(3) T(4)])
    dst = Mem.upload([0])

    @cuda kernel(
        CUDAnative.DevicePtr{T,AS.Global}(CuPtr{T}(src.ptr)),
        CUDAnative.DevicePtr{T,AS.Global}(CuPtr{T}(dst.ptr))
    )

    @test Mem.download(T, src, 4)[4] == Mem.download(T, dst)[1]
end

end

end


text:
@testset "single op" begin
    @test_throws MethodError augment!(rand(2,2), (rect,), Rotate90())
    @test_throws MethodError augment!((rand(2,2),), rect, Rotate90())
    @test_throws MethodError augment!((rand(2,2),rand(2,2)), (rect,), Rotate90())
    @test_throws BoundsError augment!(rand(2,2), rect, Rotate90())
    for pl in (Augmentor.ImmutablePipeline(Rotate90()), (Rotate90(),))
        img = @inferred Augmentor._augment(rect, pl)
        @test typeof(img) <: Array
        @test img == @inferred(augment(rect, Rotate90()))
        @test img == @inferred(augment(rect, (Rotate90(),)))
        @test eltype(img) <: eltype(rect)
        @test img == rotl90(rect)
        img1, img2 = @inferred augment((square2, rect), pl)
        @test img1 == rotl90(square2)
        @test img2 == img
        @test typeof(img2) == typeof(img)
        img3 = @inferred augment(rect, Rotate90())
        @test img2 == img3
        @test typeof(img3) == typeof(img2)
        out = similar(img)
        @test @inferred(augment!(out, rect, pl)) == img
        outs = (similar(square2), similar(img))
        res = @inferred(augment!(outs, (square2, rect), pl))
        @test res === outs
        @test outs[1] == img1
        @test outs[2] == img2
        out = similar(img)
        @test @inferred(augment!(out, rect, Rotate90())) == img
        @test_throws BoundsError augment!(rand(2,2), rect, pl)
    end
end

ops = Augmentor.ImmutablePipeline(Rotate(90),Rotate(-90)) # forces affine
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) === typeof(invwarpedview(camera, Augmentor.toaffinemap(NoOp(),rect), Flat()))
    @test wv == camera
    wv1, wv2 = @inferred Augmentor._augment((rgb_rect,camera), ops)
    @test typeof(wv1) === typeof(invwarpedview(rgb_rect, Augmentor.toaffinemap(NoOp(),rect), Flat()))
    @test typeof(wv2) === typeof(invwarpedview(camera, Augmentor.toaffinemap(NoOp(),rect), Flat()))
    @test wv1 == rgb_rect
    @test wv2 == camera
end

ops = (CacheImage(),CacheImage()) # forces eager
@testset "$(str_showcompact(ops))" begin
    img = @inferred Augmentor._augment(camera, ops)
    @test img === camera
    img1, img2 = @inferred Augmentor._augment((camera, square2), ops)
    @test img1 === camera
    @test img2 === square2
end

ops = (ElasticDistortion(4,4),CacheImage()) # forces lazy then eager
@testset "$(str_showcompact(ops))" begin
    img = @inferred Augmentor._augment(camera, ops)
    @test size(img) == size(camera)
    @test typeof(img) == typeof(camera)
    img1, img2 = @inferred Augmentor._augment((camera, N0f8.(camera)), ops)
    @test img1 isa Array{Gray{N0f8}}
    @test img2 isa Array{N0f8}
    @test img1 == img2
end

ops = (ShearX(45),ShearX(-45)) # forces affine
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) === typeof(invwarpedview(rect, Augmentor.toaffinemap(NoOp(),rect), Flat()))
    @test view(wv,1:512,1:512) == camera
end

ops = (ShearY(45),ShearY(-45)) # forces affine
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) === typeof(invwarpedview(rect, Augmentor.toaffinemap(NoOp(),rect), Flat()))
    @test view(wv,1:512,1:512) == camera
end

ops = (ShearY(45),ShearX(-2),CacheImage()) # forces affine then eager
@testset "$(str_showcompact(ops))" begin
    img = @inferred Augmentor._augment(camera, ops)
    @test typeof(img) <: OffsetArray
    @test indices(img) == (-255:768, 0:512)
end

ops = (Resize(2,2),Rotate90()) # forces affine
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(rect, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === rect
    @test round.(Float64.(wv),1) == round.(Float64.(rotl90(imresize(rect,2,2))),1)
end

ops = (Resize(2,2),Rotate90(),CacheImage()) # forces affine then eager
@testset "$(str_showcompact(ops))" begin
    img = @inferred Augmentor._augment(rect, ops)
    @test typeof(img) <: OffsetArray
    @test round.(Float64.(img),1) == round.(Float64.(rotl90(imresize(rect,2,2))),1)
end

buf = rand(Gray{N0f8}, 2, 2)
ops = (Resize(2,2),Rotate90(),CacheImage(buf)) # forces affine then eager
@testset "$(str_showcompact(ops))" begin
    img = @inferred Augmentor._augment(rect, ops)
    @test typeof(img) <: OffsetArray
    @test round.(Float64.(img),1) == round.(Float64.(rotl90(imresize(rect,2,2))),1)
    @test img == ops[3].buffer
    @test parent(img) === ops[3].buffer
end

ops = (Rotate180(),Crop(5:200,200:500),Rotate90(1),Crop(1:250, 1:150))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/rot_crop_either_crop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rot_crop_either_crop.txt" img
    out = similar(img)
    @test @inferred(augment!(out, camera, ops)) === out
    @test_reference "reference/rot_crop_either_crop.txt" out
    @test @allocated(augment!(out, camera, ops)) < @allocated(augment(camera, ops))
end

ops = Augmentor.ImmutablePipeline(Rotate180(),Crop(5:200,200:500),Rotate90(),Crop(50:300, 50:195))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(parent(wv)) <: Base.PermutedDimsArrays.PermutedDimsArray
    @test_reference "reference/rot_crop_rot_crop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rot_crop_rot_crop.txt" img
    out = similar(img)
    @test @inferred(augment!(out, camera, ops)) === out
    @test_reference "reference/rot_crop_rot_crop.txt" out
    @test @allocated(augment!(out, camera, ops)) < @allocated(augment(camera, ops))
end

ops = (Rotate180(),Crop(5:200,200:500),Rotate90(),Crop(50:300, 50:195),Resize(25,15))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test eltype(wv) <: eltype(camera)
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/rot_crop_rot_crop_resize.txt" wv
    img = @inferred Augmentor.augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rot_crop_rot_crop_resize.txt" img
    out = similar(img)
    @test @inferred(augment!(out, camera, ops)) === out
    @test_reference "reference/rot_crop_rot_crop_resize.txt" out
    @test @allocated(augment!(out, camera, ops)) < @allocated(augment(camera, ops))
end

ops = (Rotate(45),CropNative(1:512,1:512))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/rot45_crop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rot45_crop.txt" img
end

ops = (Rotate(45),CropSize(512,512))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/rot45_crop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rot45_crop.txt" img
end

ops = (Rotate(-45),CropSize(256,256))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/rotr45_cropsize.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rotr45_cropsize.txt" img
end

ops = (Scale(.1,.2),NoOp())
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: InvWarpedView
    @test parent(wv).itp.coefs === camera
    @test_reference "reference/scale_noop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/scale_noop.txt" img
end

ops = (Scale(.1,.2),CropRatio())
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/scale_cropratio.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/scale_cropratio.txt" img
end

ops = (ShearX(45),NoOp())
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: InvWarpedView
    @test parent(wv).itp.coefs === camera
    @test_reference "reference/shearx_noop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/shearx_noop.txt" img
end

ops = (ShearY(45),CropNative(1:512,1:512))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/sheary_crop.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/sheary_crop.txt" img
end

ops = (Crop(101:200,201:350),Scale(.2,.4))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/crop_scale.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/crop_scale.txt" img
end

ops = (Crop(101:200,201:350),Zoom(1.3))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/crop_zoom.txt" wv
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/crop_zoom.txt" img
end

ops = (Rotate(45),Zoom(2))
@testset "$(str_showcompact(ops))" begin
    wv = @inferred Augmentor._augment(camera, ops)
    @test typeof(wv) <: SubArray
    @test typeof(wv.indexes) <: Tuple{Vararg{IdentityRange}}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)).itp.coefs === camera
    @test_reference "reference/rot45_zoom.txt" wv
    wv2 = Augmentor._augment(camera, (ops..., NoOp()))
    @test_reference "reference/rot45_zoom.txt" wv2
    img = @inferred augment(camera, ops)
    @test img == parent(copy(wv))
    @test typeof(img) <: Array
    @test eltype(img) <: eltype(camera)
    @test_reference "reference/rot45_zoom.txt" img
    img1, img2 = @inferred augment((N0f8.(camera),camera), ops)
    @test img1 == img2
    @test img2 == img
    @test img1 isa Array{N0f8}
    @test img2 isa Array{Gray{N0f8}}
    @test_reference "reference/rot45_zoom.txt" img2
end

ops = (Rotate(45),CropSize(200,200),Zoom(1.1),ConvertEltype(RGB{Float64}),SplitChannels())
@testset "$(str_showcompact(ops))" begin
    wv1 = @inferred Augmentor._augment(camera, ops[1:3])
    wv2 = @inferred Augmentor._augment(camera, ops[1:4])
    wv3 = @inferred Augmentor._augment(camera, ops)
    _, wv4 = @inferred Augmentor._augment((rotl90(camera), camera), ops)
    @test wv3 == wv4
    @test typeof(wv3) == typeof(wv4)
    img = colorview(RGB{Float64}, wv3)
    @test RGB{Float64}.(copy(wv1)) ≈ wv2
    @test collect(wv1) ≈ img
    @test_reference "reference/rot45_crop_zoom_convert.txt" wv2
end

# just for code coverage
@test typeof(@inferred(Augmentor.augment_impl(Rotate90()))) <: Expr


text:
function get_plus_inside(obj::BoundaryAttribute)
    return jcall(obj, "getPlusInside", SubHyperplane, ())
end

function get_plus_outside(obj::BoundaryAttribute)
    return jcall(obj, "getPlusOutside", SubHyperplane, ())
end

function get_splitters(obj::BoundaryAttribute)
    return jcall(obj, "getSplitters", NodesSet, ())
end


text:
using Distributed

function calc_napier_parallel()
    n_trial =1e9
    total_count = @distributed (+) for i in 1:n_trial
        counter = 0
        accumulated = 0.0
        while true
            accumulated += rand()
            counter += 1
            if accumulated >= 1.0
                break
            end
        end
        counter
    end
    println(total_count/n_trial)
end

addprocs()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using JSON

function formatSolveBioResponse(res, raw=false)
    body = res.data
    if raw
        return body
    else
        res = JSON.parse(body)
        return res
    end
end


text:
# This module implements the MMA Algorithm in Julia
# as described in:
# AU  - Svanberg, Krister
# TI  - The method of moving asymptotes—a new method for structural optimization
# JO  - International Journal for Numerical Methods in Engineering
# JA  - Int. J. Numer. Meth. Engng.
module MMA

using Parameters, StructArrays, Setfield, TimerOutputs, GPUUtils, CuArrays

using LinearAlgebra

import Optim
import Optim: OnceDifferentiable, Fminbox, GradientDescent, update!,
                MultivariateOptimizationResults, OptimizationTrace, maxdiff,
                LineSearches, ConjugateGradient, LBFGS, AbstractOptimizer
import Base: min, max, show

export MMAModel, box!, ineq_constraint!, optimize

struct MMA87 <: AbstractOptimizer end
struct MMA02 <: AbstractOptimizer end
abstract type TopOptAlgorithm end

include("utils.jl")
include("model.jl")
include("primal.jl")
include("dual.jl")
include("lift.jl")
include("trace.jl")
include("workspace.jl")

const μ = 0.1
const ρmin = 1e-5

default_dual_caps(::MMA87, ::Type{T}) where T = (T(0.9), T(1.1))
#default_dual_caps(::MMA87, ::Type{T}) where T = (T(0.0), T(Inf))

#default_dual_caps(::MMA02, ::Type{T}) where T = (T(1), T(100))
default_dual_caps(::MMA02, ::Type{T}) where T = (T(1e6), T(1e6))

function optimize(model::MMAModel{T,TV}, x0::TV, optimizer=MMA02(), suboptimizer=Optim.ConjugateGradient(); s_init=T(0.5), s_incr=T(1.2), s_decr=T(0.7), dual_caps=default_dual_caps(optimizer, T)) where {T, TV}
    check_error(model, x0)
    workspace = MMAWorkspace(model, x0, optimizer, suboptimizer; s_init=s_init,
        s_incr=s_incr, s_decr=s_decr, dual_caps=dual_caps)
    optimize!(workspace)
end

struct MMAResult{TO, TX, T}
    optimizer::TO
    initial_x::TX
    minimizer::TX
    minimum::T
    iter::Int
    iteration_converged::Bool
    x_converged::Bool
    x_tol::T
    x_abschange::T
    f_converged::Bool
    f_tol::T
    f_abschange::T
    g_converged::Bool
    g_tol::T
    g_residual::T
    f_increased::Bool
    f_calls::Int
    g_calls::Int
    h_calls::Int
end

function optimize!(#=to, =#workspace::MMAWorkspace{T, TV, TM}) where {T, TV, TM}
    @unpack model, optimizer, suboptimizer, suboptions, x0, x, x1, x2, λ, l, u, ∇f_x, g,
        ng_approx, ∇g, f_x, f_calls, g_calls, f_x_previous, primal_data, tr, tracing,
        converged, x_converged, f_converged, gr_converged, f_increased, x_residual,
        f_residual, gr_residual, asymptotes_updater, variable_bounds_updater,
        cvx_grad_updater, lift_updater, lift_resetter, x_updater, dual_obj,
        dual_obj_grad, dual_caps, outer_iter, iter = workspace

    TSubOptions = typeof(workspace.suboptions)
    n_i = length(constraints(model))
    n_j = dim(model)
    maxiter = model.maxiter[]
    while !converged && iter < maxiter
        outer_iter += 1
        asymptotes_updater(Iteration(outer_iter))

        # Track trial points two steps back
        copyto!(x2, x1)
        copyto!(x1, x)

        # Update convex approximation
        ## Update bounds on primal variables
        variable_bounds_updater()

        ## Computes values and updates gradients of convex approximations of objective and constraints
        cvx_grad_updater()

        if optimizer isa MMA02
            lift_resetter(Iteration(outer_iter))
        end
        lift = true
        _suboptions::TSubOptions = @set suboptions.outer_iterations = model.maxiter[]
        _suboptions = @set suboptions.iterations = model.maxiter[]
        while lift && iter < model.maxiter[]
            iter += 1
            # Solve dual
            λ .= min.(dual_caps[2], max.(λ, dual_caps[1]))
            d = OnceDifferentiable(dual_obj, dual_obj_grad, λ)
            _minimizer = #=@timeit to "Fminbox"=# Optim.optimize(d, l, u, λ, Optim.Fminbox(suboptimizer), _suboptions).minimizer
            copyto!(λ, _minimizer)
            dual_obj_grad(ng_approx, λ)

            # Evaluate the objective function and its gradient
            f_x_previous, f_x = f_x, eval_objective(model, x, ∇f_x)
            f_calls, g_calls = f_calls + 1, g_calls + 1
            # Correct for functions whose gradients go to infinity at some points, e.g. √x
            while mapreduce((x)->(isinf(x) || isnan(x)), or, ∇f_x, init=false)
                map!((x1,x)->(T(0.01)*x1 + T(0.99)*x), x, x1, x)
                f_x = eval_objective(model, x, ∇f_x)
                f_calls, g_calls = f_calls + 1, g_calls + 1
            end
            primal_data.f_val[] = f_x

            # Evaluate the constraints and their Jacobian
            map!((i)->eval_constraint(model, i, x, @view(∇g[:,i])), g, 1:n_i)

            if optimizer isa MMA87
                lift = false
            else
                lift = lift_updater()
            end
        end

        # Assess convergence
        x_converged, f_converged, gr_converged,
        x_residual, f_residual, gr_residual,
        f_increased, converged = assess_convergence(x, x1, f_x, f_x_previous, ∇f_x,
            xtol(model), ftol(model), grtol(model))

        converged = converged && all((x)->(x<=ftol(model)), g)
        # Print some trace if flag is on
        @mmatrace()
    end
    h_calls = 0

    @pack! workspace = model, optimizer, suboptimizer, x0, x, x1, x2, λ, l, u, ∇f_x, g,
        ng_approx, ∇g, f_x, f_calls, g_calls, f_x_previous, primal_data, tr, tracing,
        converged, x_converged, f_converged, gr_converged, f_increased, x_residual,
        f_residual, gr_residual, asymptotes_updater, variable_bounds_updater,
        cvx_grad_updater, lift_updater, lift_resetter, x_updater, dual_obj,
        dual_obj_grad, dual_caps, outer_iter, iter

    results = MMAResult(optimizer,
            x0,
            x,
            f_x,
            iter,
            iter == model.maxiter[],
            x_converged,
            xtol(model),
            x_residual,
            f_converged,
            ftol(model),
            f_residual,
            gr_converged,
            grtol(model),
            gr_residual,
            f_increased,
            f_calls,
            g_calls,
            h_calls)
    return results
end

end # module


text:
function cos(obj::SinCos)
    return jcall(obj, "cos", jdouble, ())
end

function sin(obj::SinCos)
    return jcall(obj, "sin", jdouble, ())
end


text:
using Test, BitBasis

@test bsizeof(ind) == sizeof(Int) * 8
@test onehot(ComplexF64, 2, 2) == [0, 0, 1, 0]
@test bdistance(1, 7) == 2
@test onehot(ComplexF64, bit"01", 2) == transpose(ComplexF64[0 1 0 0; 0 1 0 0])
@test log2dim1(rand(4, 4)) == log2i(4)

@testset "log2i" begin
    for itype in [Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128]
        @test log2i(itype(2^5)) == 5
        @test typeof(log2i(itype(2^5))) == Int
    end
end

@testset "bit length" begin
    @test bit_length(8) == 4
    @test bit_length(Int32(8)) == 4
end

@testset "indices with" begin
    # indices_with
    nbit = 5
    poss, vals = [4, 2, 3], [1, 0, 1]
    @test indices_with(nbit, poss, vals) == filter(x -> readbit.(x, poss) == vals, basis(nbit))
end


text: Dict{String, Any}("g" => Dict{String, Any}("h" => Dict{String, Any}("i" => Dict{String, Any}())), "j" => Dict{String, Any}("ʞ" => Dict{String, Any}("l" => Dict{String, Any}())), "x" => Dict{String, Any}("1" => Dict{String, Any}("2" => Dict{String, Any}())), "a" => Dict{String, Any}(" x " => Dict{String, Any}(), "d.e" => Dict{String, Any}(), "b" => Dict{String, Any}("c" => Dict{String, Any}()), "b.c" => Dict{String, Any}()), "d" => Dict{String, Any}("e" => Dict{String, Any}("f" => Dict{String, Any}())))


text:
using Test
using LIKWID
using LinearAlgebra
using CUDA

@assert CUDA.functional()

# julia must have been started with `likwid-perfctr -G ... -W ... -m`
@assert GPUMarker.isactive()

# init
@test isnothing(GPUMarker.init())

N = 100_000_000
a = 3.141f0
z_gpu = CUDA.zeros(Float32, N)
x_gpu = CUDA.rand(Float32, N)
y_gpu = CUDA.rand(Float32, N)

saxpy_gpu!(z,a,x,y) = CUDA.@sync z .= a .* x .+ y

# regular workflow
saxpy_gpu!(z_gpu,a,x_gpu,y_gpu)
@test GPUMarker.startregion("saxpy_gpu!")
saxpy_gpu!(z_gpu,a,x_gpu,y_gpu)
@test GPUMarker.stopregion("saxpy_gpu!")

# nextgroup
@test isnothing(GPUMarker.nextgroup())

# registerregion
A = CUDA.rand(100,100)
B = CUDA.rand(100,100)
@test_broken GPUMarker.registerregion("mul")
@test GPUMarker.startregion("mul")
for _ in 1:10
    A * B
end
@test GPUMarker.stopregion("mul")

# TODO: fix!
# # getregion
# LIKWID.NvMon.init([0])
# x = GPUMarker.getregion("mul")
# @test typeof(x) == Tuple{Int32, Vector{Float64}, Float64, Int32}

# # resetregion
# @test GPUMarker.resetregion("mul")
# y = GPUMarker.getregion("mul")
# @test x != y
# @test y[3] == 0
# @test y[4] == 0

# close
@test isnothing(GPUMarker.close())


text:
using ModelingToolkit, OrdinaryDiffEq, Test

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

lorenz1 = ODESystem(eqs,name=:lorenz1)
lorenz2 = ODESystem(eqs,name=:lorenz2)

@variables a, α
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a*γ,
               α ~ 2lorenz1.x + a*γ]
connected = structural_simplify(ODESystem(connections,t,[a,α],[γ],systems=[lorenz1,lorenz2]))

u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0,
      a => 2.0]

p  = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8/3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8/3,
      γ => 2.0]

tspan = (0.0,100.0)
prob = ODEProblem(connected,u0,tspan,p)
sol = solve(prob,Rodas4())

@test sol[lorenz1.x] isa Vector
@test sol[lorenz1.x,2] isa Float64
@test sol[lorenz1.x,:] isa Vector
@test length(sol[lorenz1.x,1:5]) == 5
@test sol[α] isa Vector
@test sol[α,3] isa Float64
@test length(sol[α,5:10]) == 6

#=
using Plots
plot(sol,vars=(lorenz2.x,lorenz2.z))
plot(sol,vars=(α,lorenz2.z))
plot(sol,vars=(lorenz2.x,α))
plot(sol,vars=α)
plot(sol,vars=(t,α))
=#


text:
using Test, DLM

# Source: "Robustness of the subsonic doublet lattice method" by van Zyl
@testset "AR 20 - 10 X 40 - Van Zyl" begin

    AR = 20

    nchord = 10
    nspan = 40

    # define geometry
    b = 10
    semispan = b/2
    chord = b/AR
    symmetric = true

    # generate geometry
    xyz = zeros(3, nchord+1, nspan+1)
    for i = 1:nchord+1
        for j = 1:nspan+1
            xyz[1, i, j] = chord*(i-1)/nchord
            xyz[2, i, j] = semispan*(j-1)/nspan
        end
    end

    # set flow conditions
    Vinf = 10
    M = 0.0

    # set frequency range
    nfreq = 12
    kr = [0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

    # initialize matrices
    AIC = zeros(Complex{Float64}, nchord, nspan, nchord, nspan)
    w = zeros(Complex{Float64}, nchord, nspan)
    CL = zeros(Complex{Float64}, nfreq)

    for kernel in [ParabolicKernel(), QuarticKernel()]

        # loop through frequencies
        for k = 1:nfreq

            # reduced freqency is based on entire chord
            ω = Vinf*kr[k]/(chord)

            # get AIC
            influence_matrix!(AIC, ω, Vinf, M, xyz, true, kernel=kernel)

            # get downwash
            for i = 1:nchord
                for j = 1:nspan
                    # get control point x-location
                    x1 = (1/2)*xyz[1,i,j] + (1/2)*xyz[1,i,j+1]
                    x2 = (1/2)*xyz[1,i+1,j] + (1/2)*xyz[1,i+1,j+1]
                    x = (1/4)*x1 + (3/4)*x2
                    # compute downwash
                    # -1/U*(dh/dt + U*dh/dx), where h = (x - 0.5c)*e^{j*omega*t}
                    w[i,j] = -1.0 - 1im*(ω/Vinf)*(x-0.5*chord)
                end
            end

            # get pressure coefficients
            Cp = pressure_coefficients(AIC, w)

            # solve for lift coefficient
            CL[k] = sum(Cp)/(nspan*nchord)
        end

        # # 10 x 10 - extracted from plot
        # CL_van_Zyl = [5.5+0.0im, 4.4+0.2im, 3.9+1.6im, 3.5+4.1im, 3.2+6.5im, 3.0+8.5im,
        #     2.6+10.3im, 2.2+11.9im, 1.7+13.0im, 1.2+13.7im, 0.6+13.9im, -0.1+13.3im]

        # # 10 x 40 - extracted from plot
        # CL_van_Zyl = [5.5+0.0im, 4.3+0.2im, 3.8+1.5im, 3.4+4.0im, 3.2+6.2im, 2.9+8.2im,
        #     2.6+9.9im, 2.3+11.3im, 2.0+12.4im, 1.7+13.2im, 1.4+13.8im, 1.1+14.0im]

        # # extrapolated solution (grid independent) - extracted from plot
        # CL_van_Zyl = [5.45+0.0im, 4.27+0.2im, 3.80+1.6im, 3.54+4.1im, 3.46+6.4im,
        #     3.42+8.8im, 3.41+1.1im, 3.39+13.4im, 3.38+15.6im, 3.38+17.9im, 3.37+20.2im, 3.34+22.4im]

        if kernel === ParabolicKernel()
            # current outputs (validate new values before changing these)
            CL_test = [5.5+0.0im, 4.3+0.3im, 3.8+1.6im, 3.5+4.1im, 3.2+6.4im, 3.0+8.5im,
                2.7+10.3im, 2.4+11.8im, 2.0+13.1im, 1.6+14.0im, 1.2+14.6im, 0.8+14.7im]
        elseif kernel === QuarticKernel()
            # current outputs (identical to those produced by van Zyl)
            CL_test= [5.5+0.0im, 4.3+0.3im, 3.8+1.6im, 3.4+4.1im, 3.2+6.3im, 2.9+8.3im,
                2.6+10.0im, 2.3+11.4im, 2.0+12.5im, 1.7+13.3im, 1.4+13.8im, 1.1+14.1im]
        end

        # test that validated results have not changed
        for i = 1:nfreq
            @test isapprox(real(CL[i]), real(CL_test[i]), atol=0.1)
            @test isapprox(imag(CL[i]), imag(CL_test[i]), atol=0.1)
        end
    end
end


text:
#===============================================================================
# Practice using julia's type system along with dynamic, polymophic properties
# allowing julia to have its cool multiple dispatch system for specific types
#
# Author: Carl Cortright
# Date: 12/28/2016
#
===============================================================================#

# Note: Julia uses a type graph to determine inheritance

(1+2)::Int64

# In a variable declaration, it is liked typed c or java
function foo()
  x::Int8 = 100
  x
end
typeof(foo())

# This will always return a 64-bit float
function Sinc(x)::Float64
  if x == 0
    return 1
  end
  return sin(pi * x)/(pi * x)
end
Sinc(2)

# A type that is a subtype of real numbers
abstract subreal <: Real

# Note: when no subtype is given, an abstract type inherits from any, like object in java

# Playing with bitstypes
bitstype 128 genome
# Julia's type system is normative, meaning types of the same size are not interchangable

# Composite types - the equivalent to structs in c
type kungfoo
  bar
  baz::Int64
  qux::Float64
end

# A default constructor is generated when foo is created
myFoo = kungfoo("Hello world", 1, 1.1234)
typeof(myFoo)

# Fields of immutable types can't be changed
immutable complexnum
  real::Float64
  imag::Float64
end

# Unions can be used to specify a union of types
IntorString = Union{Int, AbstractString}

# Parameterized types can be used to create whole new classes of types
type point{T}
  x::T
  y::T
end

# This method can only be applied to subtypes
function norm1(p::point{Real})
  sqrt(p.x^2 + p.y^2)
end

# This is the right way to do it, so that we can apply to subtypes
function norm2{T<:Real}(p::point{T})
  sqrt(p.x^2 + p.y^2)
end

a = point{Float64}(1.0, 0.0)
typeof(a)

abstract Pointy{T}

type Pointer{T} <: Pointy{T}
  x::T
  y::T
end

Pointer{Float64} <: Pointy{Float64}

#
# DiagPoints only exist on the diagonal line y=x
#
type DiagPoint{T} <: Pointy{T}
  x::T
end

function CalcY{T <: Real}(p::DiagPoint{T})::Pointer{T}
  x = p.x
  y = x
  return Pointer(x, y)
end

# Typealiases can also be useful
typealias counter UInt

# There are also lots of functions that operate on types
supertype(Float16)
supertype(AbstractFloat)
supertype(Real)
supertype(Number)
# Finally got to Any...

# You can also put plain values in for types
myMatrix = Array{Float64, 2}

# This will allow you to dispatch on true/false (thinking indicator functions...)
immutable Val{T}
end
firstlast(::Type{Val{true}}) = "first"
firstlast(::Type{Val{false}}) = "last"
firstlast(Val{true})
firstlast(Val{false})

# Nullable objects are a way to interact with objects that could not exist
x1 = Nullable{Int64}()
isnull(x1)
# Throws a NullException
get(x1)

max_issues_repo_name: carlcortright/Playground


text:
push!(LOAD_PATH,"../src/")

using SplittingOperators
using Documenter

makedocs(modules=[SplittingOperators],
         doctest = false,
         format = :html,
         sitename = "SplittingOperators.jl",
         pages = ["Documentation"    => "index.md"])

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/JuliaVlasov/SplittingOperators.jl.git",
 )


text:
# SCIP typedefs

typealias SCIP_LINEARCONSTYPE SCIP_LinearConsType
typealias SCIP_SETPPCTYPE SCIP_SetppcType
typealias SCIP_CLOCKTYPE SCIP_ClockType
typealias SCIP_DISPSTATUS SCIP_DispStatus
typealias SCIP_EXPROP SCIP_ExprOp
typealias SCIP_EXPRCURV SCIP_ExprCurv
typealias SCIP_BRANCHDIR SCIP_BranchDir
typealias SCIP_LPSOLSTAT SCIP_LPSolStat
typealias SCIP_BOUNDTYPE SCIP_BoundType
typealias SCIP_SIDETYPE SCIP_SideType
typealias SCIP_ROWORIGINTYPE SCIP_RowOriginType
typealias SCIP_LPALGO SCIP_LPAlgo
typealias SCIP_VERBLEVEL SCIP_VerbLevel
typealias SCIP_NLPPARAM SCIP_NlpParam
typealias SCIP_NLPSOLSTAT SCIP_NlpSolStat
typealias SCIP_NLPTERMSTAT SCIP_NlpTermStat
typealias SCIP_PARAMTYPE SCIP_ParamType
typealias SCIP_PARAMSETTING SCIP_ParamSetting
typealias SCIP_PARAMEMPHASIS SCIP_ParamEmphasis
typealias SCIP_OBJSENSE SCIP_Objsense
typealias SCIP_RESULT SCIP_Result
typealias SCIP_RETCODE SCIP_Retcode
typealias SCIP_EFFICIACYCHOICE SCIP_Efficiacychoice
typealias SCIP_STAGE SCIP_Stage
typealias SCIP_SETTING SCIP_Setting
typealias SCIP_SOLORIGIN SCIP_SolOrigin
typealias SCIP_STATUS SCIP_Status
typealias SCIP_NODETYPE SCIP_NodeType
typealias SCIP_VARSTATUS SCIP_Varstatus
typealias SCIP_VARTYPE SCIP_Vartype
typealias SCIP_DOMCHGTYPE SCIP_DomchgType
typealias SCIP_BOUNDCHGTYPE SCIP_BoundchgType
typealias SCIP_VBCCOLOR SCIP_VBCColor


text:
using FixedLengthStrings
using Base.Test

test_string = "foo"

test_tiny = @tiny_string "foo"
test_short = @short_string "foo"
test_mid = @mid_string "foo"

@test test_tiny == test_string
@test test_tiny == test_short
@test test_tiny == test_mid

tiny_length = 4
short_length = 16
mid_length = 64

@test length(test_tiny.sv) == tiny_length
@test length(test_short.sv) == short_length
@test length(test_mid.sv) == mid_length

@test_throws ErrorException try @eval @tiny_string("foobar") catch err; throw(err.error) end


text:
module AlgebraOverFields

using Printf: @printf,@sprintf
using ...Interfaces: dimension
using ...Prerequisites: atol,rtol,Float,rawtype
using ...Prerequisites.NamedVectors: NamedVector
using ...Prerequisites.TypeTraits: efficientoperations
using ..Combinatorics: AbstractCombinatorics

import ...Interfaces: rank,add!,sub!,mul!,div!,⊗,⋅,sequence,permute

export SimpleID,ID
export idpropertynames,idisless
export Element,ScalarElement,Elements
export scalartype,idtype

"""
    SimpleID <: NamedVector

A simple id is the building block of the id system of an algebra over a field.
"""
abstract type SimpleID <: NamedVector end

"""
    ID{I<:SimpleID,N}

The (composite) id system of an algebra over a field.

Type alias for `NTuple{N,I} where {N,I<:SimpleID}`
"""
const ID{I<:SimpleID,N}=NTuple{N,I}
"""
    ID(ids::SimpleID...)
    ID(ids::NTuple{N,SimpleID}) where N

Get the composite id from simple ids.
"""
ID(ids::SimpleID...)=ids
ID(ids::NTuple{N,SimpleID}) where N=ids
"""
    ID(::Type{SID},attrs::Vararg{NTuple{N},M}) where {SID<:SimpleID,N,M}

Get the composite id from the components of simple ids.
"""
@generated function ID(::Type{SID},attrs::Vararg{NTuple{N,Any},M}) where {SID<:SimpleID,N,M}
    exprs=[]
    for i=1:N
        args=[:(attrs[$j][$i]) for j=1:M]
        push!(exprs,:(SID($(args...))))
    end
    return :(ID($(exprs...)))
end

"""
    idpropertynames(::Type{I}) where I<:ID{SimpleID} -> Tuple{Vararg{Symbol}}

Get the property names of a composite id.
"""
@generated function idpropertynames(::Type{I}) where I<:ID{SimpleID}
    exprs=[QuoteNode(Symbol(name,'s')) for name in I|>eltype|>fieldnames]
    return Expr(:tuple,exprs...)
end

"""
    getproperty(cid::ID{SimpleID},name::Symbol)

Get the property of a composite id.
"""
Base.getproperty(cid::ID{SimpleID},name::Symbol)=idgetproperty(cid,Val(name),Val(cid|>typeof|>idpropertynames))
@generated function idgetproperty(cid::ID{SimpleID},::Val{name},::Val{names}) where {name,names}
    index=findfirst(isequal(name),names)::Int
    exprs=[:(getfield(cid[$i],$index)) for i=1:fieldcount(cid)]
    return Expr(:tuple,exprs...)
end

"""
    show(io::IO,cid::Tuple{SimpleID,Vararg{SimpleID}})

Show a composite id.
"""
Base.show(io::IO,cid::Tuple{SimpleID,Vararg{SimpleID}})=@printf io "ID(%s)" join(cid,",")

"""
    promote_type(::Type{Tuple{}},::Type{Tuple{}})
    promote_type(::Type{Tuple{}},I::Type{<:ID{SimpleID,N}}) where N
    promote_type(I::Type{<:ID{SimpleID,N}},::Type{Tuple{}}) where N

Define the promote rule for ID types.
"""
Base.promote_type(::Type{Tuple{}},::Type{Tuple{}})=Tuple{}
Base.promote_type(::Type{Tuple{}},I::Type{<:ID{SimpleID,N}}) where N=ID{I|>eltype}
Base.promote_type(I::Type{<:ID{SimpleID,N}},::Type{Tuple{}}) where N=ID{I|>eltype}

"""
    idisless(cid1::ID{SimpleID},cid2::ID{SimpleID}) -> Bool

Compare two ids and judge whether the first is less than the second.

We define that ids with smaller ranks are always less than those with higher ranks. If two ids are of the same rank, the comparison goes just like that between tuples.
"""
function idisless(cid1::ID{SimpleID},cid2::ID{SimpleID})
    r1,r2=cid1|>rank,cid2|>rank
    r1<r2 ? true : r1>r2 ? false : isless(cid1,cid2)
end

"""
    rank(id::ID{SimpleID}) -> Int
    rank(::Type{<:ID{SimpleID}}) -> Any
    rank(::Type{<:ID{SimpleID,N}}) where N -> Int

Get the rank of a composite id.
"""
rank(id::ID{SimpleID})=id|>typeof|>rank
rank(::Type{<:ID{SimpleID}})=Any
rank(::Type{<:ID{SimpleID,N}}) where N=N

"""
    *(sid1::SimpleID,sid2::SimpleID) -> ID{SimpleID}
    *(sid::SimpleID,cid::ID{SimpleID}) -> ID{SimpleID}
    *(cid::ID{SimpleID},sid::SimpleID) -> ID{SimpleID}
    *(cid1::ID{SimpleID},cid2::ID{SimpleID}) -> ID{SimpleID}

Get the product of the id system.
"""
Base.:*(sid1::SimpleID,sid2::SimpleID)=ID(sid1,sid2)
Base.:*(sid::SimpleID,cid::ID{SimpleID})=ID(sid,cid...)
Base.:*(cid::ID{SimpleID},sid::SimpleID)=ID(cid...,sid)
Base.:*(cid1::ID{SimpleID},cid2::ID{SimpleID})=ID(cid1...,cid2...)

"""
    Element{V,I<:ID{SimpleID}}

An element of an algebra over a field.

The first and second attributes of an element must be
- `value`: the coefficient of the element
- `id::ID{SimpleID}`: the id of the element
"""
abstract type Element{V,I<:ID{SimpleID}} end

"""
    ScalarElement{V}

Identity element.
"""
const ScalarElement{V}=Element{V,Tuple{}}
scalartype(::Type{M}) where M<:Element=rawtype(M){M|>valtype,Tuple{}}

"""
    valtype(m::Element)
    valtype(::Type{<:Element})
    valtype(::Type{<:Element{V}}) where V

Get the type of the value of an element.

The result is also the type of the field over which the algebra is defined.
"""
Base.valtype(m::Element)=m|>typeof|>valtype
Base.valtype(::Type{<:Element})=Any
Base.valtype(::Type{<:Element{V}}) where V=V

"""
    idtype(m::Element)
    idtype(::Type{<:Element{V,I} where V}) where I<:ID{SimpleID}
    idtype(::Type{<:Element})

The type of the id of an element.
"""
idtype(m::Element)=m|>typeof|>idtype
idtype(::Type{<:Element{V,I} where V}) where I<:ID{SimpleID}=I
idtype(::Type{<:Element})=ID{SimpleID}

"""
    rank(m::Element) -> Int
    rank(::Type{M}) where M<:Element -> Int

Get the rank of an element.
"""
rank(m::Element)=m|>typeof|>rank
rank(::Type{M}) where M<:Element=M|>idtype|>rank

"""
    ==(m1::Element,m2::Element) -> Bool
    isequal(m1::Element,m2::Element) -> Bool

Compare two elements and judge whether they are equal to each other.
"""
Base.:(==)(m1::Element,m2::Element) = ==(efficientoperations,m1,m2)
Base.isequal(m1::Element,m2::Element)=isequal(efficientoperations,m1,m2)

"""
    isapprox(m1::Element,m2::Element;atol::Real=atol,rtol::Real=rtol) -> Bool

Compare two elements and judge whether they are inexactly equivalent to each other.
"""
Base.isapprox(m1::Element,m2::Element;atol::Real=atol,rtol::Real=rtol)=isapprox(efficientoperations,Val((:value,)),m1,m2;atol=atol,rtol=rtol)

"""
    replace(m::Element;kwargs...) -> typeof(m)

Return a copy of a concrete `Element` with some of the field values replaced by the keyword arguments.
"""
Base.replace(m::Element;kwargs...)=replace(efficientoperations,m;kwargs...)

"""
    promote_rule(::Type{M1},::Type{M2}) where {M1<:Element,M2<:Element}

Define the promote rule for Element types.
"""
function Base.promote_rule(::Type{M1},::Type{M2}) where {M1<:Element,M2<:Element}
    M1<:M2 && return M2
    M2<:M1 && return M1
    v1,i1,r1=M1|>valtype,M1|>idtype,M1|>rank
    v2,i2,r2=M2|>valtype,M2|>idtype,M2|>rank
    V,I=promote_type(v1,v2),promote_type(i1,i2)
    M=r2==0 ? rawtype(M1) : r1==0 ? rawtype(M2) : typejoin(rawtype(M1),rawtype(M2))
    return isconcretetype(I) ? M{V,I} : M{V,<:I}
end

"""
    mul_promote_type(::Type{M},::Type{V}) where {M<:Element,V<:Number}
    mul_promote_type(::Type{V},::Type{M}) where {M<:Element,V<:Number}

Define the promote rule for the multiplication between an Element and a scalar.
"""
mul_promote_type(::Type{M},::Type{V}) where {M<:Element,V<:Number}=mul_promote_type(V,M)
function mul_promote_type(::Type{V},::Type{M}) where {M<:Element,V<:Number}
    C,I=promote_type(M|>valtype,V),M|>idtype
    return isconcretetype(I) ? rawtype(M){C,I} : rawtype(M){C,<:I}
end

"""
    one(::Type{M}) where {M<:Element}

Get the identity operator.
"""
function Base.one(::Type{M}) where {M<:Element}
    rtype=rawtype(M)
    vtype=isconcretetype(valtype(M)) ? valtype(M) : Int
    @assert fieldnames(rtype)==(:value,:id) "one error: not supproted type($(nameof(rtype)))."
    return rtype(one(vtype),ID())
end

"""
    convert(::Type{M},m::ScalarElement) where {M<:ScalarElement}
    convert(::Type{M},m::Number) where {M<:ScalarElement}
    convert(::Type{M},m::Element) where {M<:Element}

1) Convert a scalar element from one type to another;
2) Convert a scalar to a scalar element;
3) Convert an element from one type to another.
"""
function Base.convert(::Type{M},m::ScalarElement) where {M<:ScalarElement}
    typeof(m)<:M && return m
    @assert fieldnames(M)==(:value,:id) "convert error: not supported type($(nameof(M)))."
    return rawtype(M)(convert(M|>valtype,m.value),m.id)
end
function Base.convert(::Type{M},m::Number) where {M<:ScalarElement}
    @assert fieldnames(M)==(:value,:id) "convert error: not supported type($(nameof(M)))."
    return rawtype(M)(convert(M|>valtype,m),ID())
end
function Base.convert(::Type{M},m::Element) where {M<:Element}
    typeof(m)<:M && return m
    @assert idtype(m)<:idtype(M) "convert error: dismatched ID type."
    @assert rawtype(typeof(m))<:rawtype(M) "convert error: dismatched raw Element type."
    return replace(m,value=convert(valtype(M),m.value))
end

"""
    Elements{I<:ID{SimpleID},M<:Element} <: AbstractDict{I,M}

An set of elements of an algebra over a field.

Type alias for `Dict{I<:ID{SimpleID},M<:Element}`.
Similar iterms are automatically merged thanks to the id system.
"""
const Elements{I<:ID{SimpleID},M<:Element}=Dict{I,M}
"""
    Elements(ms)
    Elements(ms::Pair{I,M}...) where {I<:ID{SimpleID},M<:Element}
    Elements(ms::Element...)

Get the set of elements with similar items merged.
"""
Elements(ms)=Base.dict_with_eltype((K,V)->Dict{K,V},ms,eltype(ms))
function Elements(ms::Pair{I,M}...) where {I<:ID{SimpleID},M<:Element}
    result=Elements{I,M}()
    for (id,m) in ms result[id]=m end
    return result
end
function Elements(ms::Element...)
    result=Elements{ms|>eltype|>idtype,ms|>eltype}()
    for m in ms add!(result,m) end
    return result
end

"""
    show(io::IO,ms::Elements)

Show a set of elements.
"""
function Base.show(io::IO,ms::Elements)
    @printf io "Elements with %s entries:\n" length(ms)
    for m in values(ms)
        @printf io "  %s\n" m
    end
end

"""
    repr(ms::Elements) -> String

Get the repr representation of a set of elements.
"""
function Base.repr(ms::Elements)
    cache=[@sprintf("Elements with %s entries:",length(ms))]
    for m in ms|>values
        push!(cache,@sprintf("  %s",repr(m)))
    end
    return join(cache,"\n")
end

"""
    zero(ms::Elements) -> Nothing
    zero(::Type{<:Elements}) -> Nothing

Get a zero set of elements.

A zero set of elements is defined to be the one with no elements.
"""
Base.zero(ms::Elements)=ms|>typeof|>zero
Base.zero(::Type{<:Elements})=nothing

"""
    ==(ms::Elements,::Nothing) -> Bool
    ==(::Nothing,ms::Elements) -> Bool
    isequal(ms::Elements,::Nothing) -> Bool
    isequal(::Nothing,ms::Elements) -> Bool
"""
Base.:(==)(ms::Elements,::Nothing)=length(ms)==0
Base.:(==)(::Nothing,ms::Elements)=length(ms)==0
Base.isequal(ms::Elements,::Nothing)=length(ms)==0
Base.isequal(::Nothing,ms::Elements)=length(ms)==0

"""
    add!(ms::Elements) -> typeof(ms)
    add!(ms::Elements,::Nothing) -> typeof(ms)
    add!(ms::Elements,m::Number) -> typeof(ms)
    add!(ms::Elements,m::Element) -> typeof(ms)
    add!(ms::Elements,mms::Elements) -> typeof(ms)

Get the inplace addition of elements to a set.
"""
add!(ms::Elements)=ms
add!(ms::Elements,::Nothing)=ms
add!(ms::Elements,m::Number)=add!(ms,convert(scalartype(ms|>valtype),m))
function add!(ms::Elements,m::Element)
    m=convert(ms|>valtype,m)
    old=get(ms,m.id,nothing)
    new=old===nothing ? m : replace(old,value=old.value+m.value)
    abs(new.value)==0.0 ? delete!(ms,m.id) : ms[m.id]=new
    return ms
end
add!(ms::Elements,mms::Elements)=(for m in mms|>values add!(ms,m) end; ms)

"""
    sub!(ms::Elements) -> typeof(ms)
    sub!(ms::Elements,::Nothing) -> typeof(ms)
    add!(ms::Elements,m::Number) -> typeof(ms)
    sub!(ms::Elements,m::Element) -> typeof(ms)
    sub!(ms::Elements,mms::Elements) -> typeof(ms)

Get the inplace subtraction of elements from a set.
"""
sub!(ms::Elements)=ms
sub!(ms::Elements,::Nothing)=ms
sub!(ms::Elements,m::Number)=add!(ms,convert(scalartype(ms|>valtype),-m))
function sub!(ms::Elements,m::Element)
    m=convert(ms|>valtype,m)
    old=get(ms,m.id,nothing)
    new=old==nothing ? -m : replace(old,value=old.value-m.value)
    abs(new.value)==0.0 ? delete!(ms,m.id) : ms[m.id]=new
    return ms
end
sub!(ms::Elements,mms::Elements)=(for m in mms|>values sub!(ms,m) end; ms)

"""
    mul!(ms::Elements,factor::ScalarElement) -> Elements
    mul!(ms::Elements,factor::Number) -> Elements

Get the inplace multiplication of elements with a scalar.
"""
mul!(ms::Elements,factor::ScalarElement)=mul!(ms,factor.value)
function mul!(ms::Elements,factor::Number)
    @assert isa(one(ms|>valtype|>valtype)*factor,ms|>valtype|>valtype) "mul! error: dismatched type, $(ms|>valtype) and $(factor|>typeof)."
    for m in values(ms)
        ms[m.id]=replace(m,value=m.value*factor)
    end
    return ms
end

"""
    div!(ms::Elements,factor::ScalarElement) -> Elements
    div!(ms::Elements,factor::Number) -> Elements

Get the inplace division of element with a scalar.
"""
div!(ms::Elements,factor::ScalarElement)=mul!(ms,1/factor.value)
div!(ms::Elements,factor::Number)=mul!(ms,1/factor)

"""
    +(m::Element) -> typeof(m)
    +(m::Element,::Nothing) -> typeof(m)
    +(::Nothing,m::Element) -> typeof(m)
    +(m::Element,factor::Number) -> Elements
    +(factor::Number,m::Element) -> Elements
    +(m1::Element,m2::Element) -> Elements
    +(ms::Elements) -> typeof(ms)
    +(ms::Elements,::Nothing) -> typeof(ms)
    +(::Nothing,ms::Elements) -> typeof(ms)
    +(ms::Elements,factor::Number) -> Elements
    +(factor::Number,ms::Elements) -> Elements
    +(ms::Elements,m::Element) -> Elements
    +(m::Element,ms::Elements) -> Elements
    +(ms1::Elements,ms2::Elements) -> Elements

Overloaded `+` operator between elements of an algebra over a field.
"""
Base.:+(m::Element)=m
Base.:+(ms::Elements)=ms
Base.:+(m::Element,::Nothing)=m
Base.:+(::Nothing,m::Element)=m
Base.:+(ms::Elements,::Nothing)=ms
Base.:+(::Nothing,ms::Elements)=ms
Base.:+(factor::Number,m::Element)=m+rawtype(m|>typeof)(factor,ID())
Base.:+(m::Element,factor::Number)=m+rawtype(m|>typeof)(factor,ID())
Base.:+(factor::Number,m::ScalarElement)=replace(m,value=m.value+factor)
Base.:+(m::ScalarElement,factor::Number)=replace(m,value=m.value+factor)
Base.:+(m1::ScalarElement,m2::ScalarElement)=replace(m1,value=m1.value+m2.value)
Base.:+(factor::Number,ms::Elements)=ms+rawtype(ms|>valtype)(factor,ID())
Base.:+(ms::Elements,factor::Number)=ms+rawtype(ms|>valtype)(factor,ID())
function Base.:+(m1::Element,m2::Element)
    I=promote_type(m1|>idtype,m2|>idtype)
    M=promote_type(typeof(m1),typeof(m2))
    return add!(Elements{I,M}(m1.id=>m1),m2)
end
Base.:+(m::Element,ms::Elements)=ms+m
function Base.:+(ms::Elements,m::Element)
    I=promote_type(ms|>keytype,m|>idtype)
    M=promote_type(valtype(ms),typeof(m))
    return add!(Elements{I,M}(ms),m)
end
function Base.:+(ms1::Elements,ms2::Elements)
    I=promote_type(ms1|>keytype,ms2|>keytype)
    M=promote_type(valtype(ms1),valtype(ms2))
    return add!(Elements{I,M}(ms1),ms2)
end

"""
    *(m::Element,::Nothing) -> Nothing
    *(::Nothing,m::Element) -> Nothing
    *(factor::Number,m::Element) -> Element
    *(m::Element,factor::Number) -> Element
    *(m1::Element,m2::Element) -> Element
    *(ms::Elements,::Nothing) -> Nothing
    *(::Nothing,ms::Elements) -> Nothing
    *(factor::Number,ms::Elements) -> Elements
    *(ms::Elements,factor::Number) -> Elements
    *(m::Element,ms::Elements) -> Elements
    *(ms::Elements,m::Element) -> Elements
    *(ms1::Elements,ms2::Elements) -> Elements

Overloaded `*` operator for element-scalar multiplications and element-element multiplications of an algebra over a field.
"""
Base.:*(m::Element,::Nothing)=nothing
Base.:*(::Nothing,m::Element)=nothing
Base.:*(ms::Elements,::Nothing)=nothing
Base.:*(::Nothing,ms::Elements)=nothing
Base.:*(factor::ScalarElement,m::Element)=m*factor.value
Base.:*(m::Element,factor::ScalarElement)=m*factor.value
Base.:*(m1::ScalarElement,m2::ScalarElement)=replace(m1,value=m1.value*m2.value)
Base.:*(factor::Number,m::Element)=m*factor
Base.:*(m::Element,factor::Number)=replace(m,value=factor*m.value)
Base.:*(factor::ScalarElement,ms::Elements)=ms*factor.value
Base.:*(ms::Elements,factor::ScalarElement)=ms*factor.value
Base.:*(factor::Number,ms::Elements)=ms*factor
function Base.:*(ms::Elements,factor::Number)
    abs(factor)==0 && return zero(Elements)
    result=Elements{ms|>keytype,mul_promote_type(ms|>valtype,factor|>typeof)}()
    for (id,m) in ms result[id]=m*factor end
    return result
end
Base.:*(m::Element,ms::Elements)=Elements((m*mm for mm in ms|>values)...)
Base.:*(ms::Elements,m::Element)=Elements((mm*m for mm in ms|>values)...)
Base.:*(ms1::Elements,ms2::Elements)=Elements((m1*m2 for m1 in ms1|>values for m2 in ms2|>values)...)
function Base.:*(m1::Element,m2::Element)
    @assert(    m1|>typeof|>nameof==m2|>typeof|>nameof && m1|>typeof|>fieldcount==m2|>typeof|>fieldcount==2,
                "\"*\" error: not implemented between $(m1|>typeof|>nameof) and $(m2|>typeof|>nameof)."
                )
    rawtype(typeof(m1))(m1.value*m2.value,m1.id*m2.id)
end

"""
    -(m::Element) -> typeof(m)
    -(m::Element,::Nothing) -> typeof(m)
    -(::Nothing,m::Element) -> typeof(m)
    -(m::Element,factor::Number) -> Elements
    -(factor::Number,m::Element) -> Elements
    -(m1::Element,m2::Element) -> Elements
    -(ms::Elements) -> typeof(ms)
    -(ms::Elements,::Nothing) -> typeof(ms)
    -(::Nothing,ms::Elements) -> typeof(ms)
    -(ms::Elements,factor::Number) -> Elements
    -(factor::Number,ms::Elements) -> Elements
    -(m::Element,ms::Elements) -> Elements
    -(ms::Elements,m::Element) -> Elements
    -(ms1::Elements,ms2::Elements) -> Elements

Overloaded `-` operator between elements of an algebra over a field.
"""
Base.:-(m::Element)=m*(-1)
Base.:-(ms::Elements)=ms*(-1)
Base.:-(m::Element,::Nothing)=m
Base.:-(::Nothing,m::Element)=-m
Base.:-(ms::Elements,::Nothing)=ms
Base.:-(::Nothing,ms::Elements)=-ms
Base.:-(factor::Number,m::Element)=rawtype(m|>typeof)(factor,ID())-m
Base.:-(m::Element,factor::Number)=m+rawtype(m|>typeof)(-factor,ID())
Base.:-(factor::Number,m::ScalarElement)=replace(m,value=factor-m.value)
Base.:-(m::ScalarElement,factor::Number)=replace(m,value=m.value-factor)
Base.:-(m1::ScalarElement,m2::ScalarElement)=replace(m1,value=m1.value-m2.value)
Base.:-(factor::Number,ms::Elements)=rawtype(ms|>valtype)(factor,ID())-ms
Base.:-(ms::Elements,factor::Number)=ms+rawtype(ms|>valtype)(-factor,ID())
function Base.:-(m1::Element,m2::Element)
    I=promote_type(m1|>idtype,m2|>idtype)
    M=promote_type(typeof(m1),typeof(m2))
    return sub!(Elements{I,M}(m1.id=>m1),m2)
end
function Base.:-(m::Element,ms::Elements)
    I=promote_type(m|>idtype,ms|>keytype)
    M=promote_type(typeof(m),valtype(ms))
    return sub!(Elements{I,M}(m.id=>m),ms)
end
function Base.:-(ms::Elements,m::Element)
    I=promote_type(ms|>keytype,m|>idtype)
    M=promote_type(valtype(ms),typeof(m))
    return sub!(Elements{M|>idtype,M}(ms),m)
end
function Base.:-(ms1::Elements,ms2::Elements)
    I=promote_type(ms1|>keytype,ms2|>keytype)
    M=promote_type(valtype(ms1),valtype(ms2))
    return sub!(Elements{I,M}(ms1),ms2)
end

"""
    /(m::Element,factor::Number) -> Element
    /(m::Element,factor::ScalarElement) -> Element
    /(ms::Elements,factor::Number) -> Elements
    /(ms::Elements,factor::ScalarElement) -> Elements

Overloaded `/` operator for element-scalar division of an algebra over a field.
"""
Base.:/(m::Element,factor::Number)=m*(1/factor)
Base.:/(m::Element,factor::ScalarElement)=m*(1/factor.value)
Base.:/(ms::Elements,factor::Number)=ms*(1/factor)
Base.:/(ms::Elements,factor::ScalarElement)=ms*(1/factor.value)

"""
    ^(m::Element,n::Integer) -> Element
    ^(ms::Elements,n::Integer) -> Elements

Overloaded `^` operator for element-integer power of an algebra over a field.
"""
Base.:^(m::Element,n::Integer)=(@assert n>0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->m,Val(n))))
Base.:^(ms::Elements,n::Integer)=(@assert n>0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->ms,Val(n))))

"""
    ⊗(m::Element,ms::Elements) -> Elements
    ⊗(ms::Elements,m::Element) -> Elements
    ⊗(ms1::Elements,ms2::Elements) -> Elements

Overloaded `⊗` operator for element-element multiplications of an algebra over a field.
"""
⊗(m::Element,ms::Elements)=Elements((m⊗mm for mm in ms|>values)...)
⊗(ms::Elements,m::Element)=Elements((mm⊗m for mm in ms|>values)...)
⊗(ms1::Elements,ms2::Elements)=Elements((m1⊗m2 for m1 in ms1|>values for m2 in ms2|>values)...)

"""
    ⋅(m::Element,ms::Elements) -> Elements
    ⋅(ms::Elements,m::Element) -> Elements
    ⋅(ms1::Elements,ms2::Elements) -> Elements

Overloaded `⋅` operator for element-element multiplications of an algebra over a field.
"""
⋅(m::Element,ms::Elements)=Elements((m⋅mm for mm in ms|>values)...)
⋅(ms::Elements,m::Element)=Elements((mm⋅m for mm in ms|>values)...)
⋅(ms1::Elements,ms2::Elements)=Elements((m1⋅m2 for m1 in ms1|>values for m2 in ms2|>values)...)

"""
    sequence(m::Element,table=nothing) -> NTuple{rank(m),Int}

Get the sequence of the ids of an element according to a table.
"""
@generated sequence(m::Element,table)=Expr(:tuple,[:(get(table,m.id[$i],nothing)) for i=1:rank(m)]...)

"""
    split(m::Element) -> Tuple{Any,Vararg{Element}}

Split an element into the coefficient and a sequence of rank-1 elements.
"""
@generated function Base.split(m::Element)
    @assert m|>fieldnames==(:value,:id) "split error: not supported split of $(nameof(typeof(m)))."
    exprs=[:(m.value)]
    for i=1:rank(m)
        push!(exprs,:(rawtype(typeof(m))(one(m|>valtype),ID(m.id[$i]))))
    end
    return Expr(:tuple,exprs...)
end

"""
    replace(m::Element,pairs::Pair{<:SimpleID,<:Union{Element,Elements}}...) -> Element/Elements
    replace(ms::Elements,pairs::Pair{<:SimpleID,<:Union{Element,Elements}}...) -> Elements

Replace the rank-1 components of an element with new element/elements.
"""
function Base.replace(m::Element,pairs::Pair{<:SimpleID,<:Union{Element,Elements}}...)
    @assert m|>typeof|>fieldnames==(:value,:id) "replace error: not supported replacement of $(nameof(typeof(m)))."
    replacedids=NTuple{length(pairs),idtype(m)|>eltype}(pair.first for pair in pairs)
    ms=split(m)
    result=ms[1]
    for i=1:rank(m)
        index=findfirst(isequal(m.id[i]),replacedids)
        result=result*(isa(index,Int) ? pairs[index].second : ms[i+1])
    end
    return result
end
function Base.replace(ms::Elements,pairs::Pair{<:SimpleID,<:Union{Element,Elements}}...)
    result=elementstype(pairs...)()
    for m in values(ms)
        add!(result,replace(m,pairs...))
    end
    return result
end
@generated function elementstype(ms::Vararg{Pair{<:SimpleID,<:Union{Element,Elements}},N}) where N
    ms=ntuple(i->(fieldtype(ms[i],2)<:Elements) ? fieldtype(ms[i],2)|>valtype : fieldtype(ms[i],2),Val(N))
    @assert mapreduce(m->(m|>fieldnames)==(:value,:id),&,ms) "elementstype error: not supported."
    I=mapreduce(idtype,promote_type,ms)|>eltype
    V=mapreduce(valtype,promote_type,ms)
    M=reduce(typejoin,ms)|>rawtype
    isconcretetype(V) && return Elements{ID{I},M{V,<:ID{I}}}
end

"""
    permute!(result::Elements,m::Element,table=nothing) -> Elements
    permute!(result::Elements,ms::Elements,table=nothing) -> Elements

Permute the ids of an-element/a-set-of-elements to the descending order according to a table, and store the permuted elements in result.
"""
function Base.permute!(result::Elements,m::Element,table=nothing)
    cache=valtype(result)[m]
    while length(cache)>0
        current=pop!(cache)
        pos=elementcommuteposition(sequence(current,table))
        if isa(pos,Nothing)
            add!(result,current)
        else
            left,right=elementleft(current,pos),elementright(current,pos)
            for middle in permute(typeof(current),current.id[pos],current.id[pos+1],table)
                temp=left*middle*right
                temp===nothing || push!(cache,temp)
            end
        end
    end
    return result
end
function Base.permute!(result::Elements,ms::Elements,table=nothing)
    for m in values(ms)
        permute!(result,m,table)
    end
    return result
end
function elementcommuteposition(seqs)
    pos=1
    while pos<length(seqs)
        seqs[pos]<seqs[pos+1] && return pos
        pos+=1
    end
    return nothing
end
elementleft(m::Element,i::Int)=rawtype(typeof(m))(m.value,m.id[1:i-1])
elementright(m::Element,i::Int)=rawtype(typeof(m))(one(m.value),m.id[i+2:end])

"""
    permute(m::Element,table=nothing) -> Elements
    permute(ms::Elements,table=nothing) -> Elements

Permute the ids of an-element/a-set-of-elements to the descending order according to a table.
"""
function permute(m::Element,table=nothing)
    M,V,I=m|>typeof|>rawtype,m|>valtype,m|>idtype|>eltype
    permute!(Elements{ID{I},M{V,<:ID{I}}}(),m,table)
end
function permute(ms::Elements,table=nothing)
    M,V,I=ms|>valtype|>rawtype,ms|>valtype|>valtype,ms|>keytype|>eltype
    permute!(Elements{ID{I},M{V,<:ID{I}}}(),ms,table)
end

end #module


text:
# Get Model properties from data in HDF5 format ?

function getModelProperties()

    Dict("prop" => "value")
end


text:
__precompile__()

module FilePaths

using Reexport, URIParser, MacroTools
@reexport using FilePathsBase

include("uri.jl")
include("compat.jl")

end # end of module


text:
####
#### Defines list of developer documents
####

dev_docs = Any[
    "Coding style" => "DevDocs/CodeStyle.md",
    "Acceptable Unicode" => "DevDocs/AcceptableUnicode.md",
    "Variable list" => "DevDocs/VariableList.md",
    "Diagnostic variable list" => "DevDocs/DiagnosticVariables.md",
]


text:
## Some Polynomial utilities

## conveniences for AbstractAlgebra.jl
# get value of GF(p) in (-p/2, p/2)
#function value(x::AbstractAlgebra.gfelem)
function value(x)
    # reach in to get a, q
    if hasfield(typeof(parent(x)), :p)
        q = parent(x).p
        a = x.d
        a <= q ÷ 2 ? a : a - q
    else
        x
    end
end

#value(x) = x

function poly_coeffs(f::AbstractAlgebra.Generic.Poly)
#    println("poly coef f=$f d=$(degree(f))")
    d = degree(f)
    d < 0 && return [value(coeff(f, 0))]
    [value(coeff(f, i)) for i in 0:degree(f)]
end

poly_coeffs(as) = as

variable(f::AbstractAlgebra.Generic.Poly) = gen(parent(f))

function maxnorm(f::AbstractAlgebra.Generic.Poly)
    norm(poly_coeffs(f), Inf)
end

function onenorm(f::AbstractAlgebra.Generic.Poly)
    norm(poly_coeffs(f), 1)
end

## Conversion to Generic.Poly
# If f = a0 + a1*y + ... an*y^n then compute
# a0*x + ... an*x^n for x a poly. (Could be different ring, ...)
function as_poly(f::AbstractAlgebra.Generic.Poly, x::AbstractAlgebra.Generic.Poly)
    degree(f) < 0 && return zero(x)
    sum(value(coeff(f, i)) * x^i for i in 0:degree(f))
end

function as_poly(as, x::AbstractAlgebra.Generic.Poly)
    sum( ai * x^(i-1) for  (i,ai) in enumerate(as))
end

# make poly over Z?
function as_poly(f::AbstractAlgebra.Generic.Poly{T}, x::String="x") where {T}
    R,y = ZZ[x]
    sum(value(coeff(f, i)) * y^i for i in 0:degree(f))
end

function as_poly(as, x0::String="x")
    R, x = ZZ[x0]
    as_poly(as, x)
end

# Z[x] with coefficients in (-p/2, p/2)
function _modp(x, p)
    a = mod(x,p)
    a > p ÷ 2 ? a - p : a
end

# make poly over x with coefficients in (-p/2...p/2)
function as_poly_modp(f, p, x::AbstractAlgebra.Generic.Poly)
    as = poly_coeffs(f)
    as_poly(_modp.(as, p), x)
end

# make poly in GP from coefficiens of x
# !!! NOTE
#     ?GF When check == false, no check is made, but the behaviour of the
#     resulting object is undefined if p is composite.
function as_poly_modp(f, p, x0::String="x")
    R, x = GF(p, check=false)[x0]
    as = poly_coeffs(f)
    as_poly(_modp.(as, p), x)
end

# Create a polynomial over Zp[x] from the coefficients
# as is (a0,a1, ..., an) as an interable
function as_poly_Zp(f::AbstractAlgebra.Generic.Poly, p, x)
    as_poly_Zp(poly_coeffs(f), p, x)
end

function as_poly_Zp(as, p::S, var="x") where {S}
    T = eltype(as)
    R = promote_type(T, S)
    q = convert(R, p)

    F = GF(q)
    R,x = PolynomialRing(F, var)
    sum(ai*x^(i-1) for (i,ai) in enumerate(as))
end