text:
module CoilsPlot

using Statistics
using LightGraphs
using ProgressMeter
using PyPlot
using PyCall

using Coils

export plot_vertex, plot_vertices, plot_edges, plot_system,
    plot_cell, plot_edge_current, plot_loop, plot_simpleloop,
    save_plots_simpleloops, save_report, plot_edge_currents,
    plot_cells, plot_loops_field, plot_deviation_histogram

# Would like to do:
# @pyimport matplotlib.patheffects as patheffects
# @pyimport arrow3d
# but annot use precompilation with @pyimport
# ref. https://github.com/JuliaPy/PyCall.jl#using-pycall-from-julia-modules
function __init__()
    global patheffects, arrow3d
    patheffects = pywrap(pyimport("matplotlib.patheffects"))
    pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(dirname(@__FILE__)))
    arrow3d = pywrap(pyimport("arrow3d"))

    using3D()
end

# To be able to write linspace.(starts, stops, lengths)
linspace(start, stop, length) = range(start, stop = stop, length = length)
signif(x, n) = round(x, sigdigits = n)

function plot_poi(poi; standalone = true)
    standalone && subplot(111, projection = "3d")

    # plot needs array of xs, then ys and then zs
    gca().plot(getindex.(poi, 1), getindex.(poi, 2), getindex.(poi, 3), ".")

    if standalone
        xlabel("x")
        ylabel("y")
        # tight_layout()
        equalize_3d_aspect()
    end
end

function plot_vertex(vertex_position; standalone = true, label = "")
    standalone && subplot(111, projection = "3d")

    gca().plot([vertex_position[1]], [vertex_position[2]], [vertex_position[3]],
        "o", color = "black", markersize = 1)
    if label != ""
        txt = gca().text(vertex_position...,  label,
            horizontalalignment = "center", verticalalignment = "center")
        txt.set_path_effects([patheffects.Stroke(linewidth = 3, foreground = "white"),
                           patheffects.Normal()])
    end

    if standalone
        xlabel("x")
        ylabel("y")
        # tight_layout()
        equalize_3d_aspect()
    end
end

function plot_vertices(vertex_positions; standalone = true, labels = true)
    standalone && subplot(111, projection = "3d")

    for (i, vertex_position) in enumerate(vertex_positions)
        plot_vertex(vertex_position, standalone = false, label = labels ? "$i" : "")
    end

    if standalone
        xlabel("x")
        ylabel("y")
        # tight_layout()
        equalize_3d_aspect()
    end
end

function plot_edges(g, vertex_positions; standalone = true, alpha = 1)
    standalone && subplot(111, projection = "3d")

    for edge in edges(g)
        pos1 = vertex_positions[src(edge)]
        pos2 = vertex_positions[dst(edge)]
        gca().plot([pos1[1], pos2[1]], [pos1[2], pos2[2]], [pos1[3], pos2[3]],
            color = "black", alpha = alpha)
    end

    if standalone
        xlabel("x")
        ylabel("y")
        # tight_layout()
        equalize_3d_aspect()
    end
end

function plot_system(g, vertex_positions, poi; standalone = true)
    standalone && subplot(111, projection = "3d")

    plot_poi(poi, standalone = false)
    plot_vertices(vertex_positions, standalone = false)
    plot_edges(g, vertex_positions, standalone = false)

    if standalone
        xlabel("x")
        ylabel("y")
        # tight_layout()
        equalize_3d_aspect()
    end
end

function plot_cell(cell, vertex_positions; shrink = 0.8, label = "")
    vrts = deepcopy(vertex_positions[cell])

    xmean = mean([ vrts[ivrts][1] for ivrts in eachindex(vrts) ])
    ymean = mean([ vrts[ivrts][2] for ivrts in eachindex(vrts) ])
    zmean = mean([ vrts[ivrts][3] for ivrts in eachindex(vrts) ])

    shrunkvrts = deepcopy(vrts)

    for ivrts in eachindex(vrts)
        shrunkvrts[ivrts][1] = (vrts[ivrts][1] - xmean) * shrink + xmean
        shrunkvrts[ivrts][2] = (vrts[ivrts][2] - ymean) * shrink + ymean
        shrunkvrts[ivrts][3] = (vrts[ivrts][3] - zmean) * shrink + zmean
    end

    for k in 1:length(shrunkvrts)-1
        for k in 1:length(shrunkvrts)-1
            x = [ shrunkvrts[k][1], shrunkvrts[k + 1][1] ]
            y = [ shrunkvrts[k][2], shrunkvrts[k + 1][2] ]
            z = [ shrunkvrts[k][3], shrunkvrts[k + 1][3] ]

            arrow3d.add_arrow3D(x, y, z, color = "C0")
        end
    end

    if label != ""
        txt = gca().text(xmean, ymean, zmean, label,
            horizontalalignment = "center", verticalalignment = "center", color = "C0")
        txt.set_path_effects([patheffects.Stroke(linewidth = 3, foreground = "white"),
                           patheffects.Normal()])
    end
end

function plot_cells(cells, vertex_positions, currents = [])
    for icell in eachindex(cells)
        label = "$icell"
        if length(currents) > 0
            label *= ": $(signif(currents[icell], 6))"
        end
        plot_cell(cells[icell], vertex_positions, label = label, shrink = 0.9)
    end
end

function plot_edge_current(edge, current, vertex_positions)
    pos1 = vertex_positions[src(edge)]
    pos2 = vertex_positions[dst(edge)]
    arrow3d.add_arrow3D([pos1[1], pos2[1]], [pos1[2], pos2[2]], [pos1[3], pos2[3]],
        color = "black")

    txt = gca().text(mean([pos1, pos2])...,  "$(signif(current, 4))",
        horizontalalignment = "center", verticalalignment = "center")
    txt.set_path_effects([patheffects.Stroke(linewidth = 3, foreground = "white"),
                       patheffects.Normal()])
end

function plot_edge_currents(g, currents, vertex_positions)
    for edge in edges(g)
        plot_edge_current(edge, currents[edge], vertex_positions)
    end
end

function plot_loop(loop, vertex_positions; color = nothing, standalone = true)
    for k in 1:length(loop)-1
        v1 = loop[k]
        v2 = loop[k + 1]

        x = [ vertex_positions[v1][1], vertex_positions[v2][1] ]
        y = [ vertex_positions[v1][2], vertex_positions[v2][2] ]
        z = [ vertex_positions[v1][3], vertex_positions[v2][3] ]

        arrow3d.add_arrow3D(x, y, z, color = color)
    end
end

function equalize_3d_aspect()
    limmin = minimum([ l[1] for l in [xlim(), ylim(), zlim()] ])
    limmax = maximum([ l[2] for l in [xlim(), ylim(), zlim()] ])
    xlim(limmin, limmax)
    ylim(limmin, limmax)
    zlim(limmin, limmax)
end

function decompose_string(x, coeffs, elements)
    "$(signif(Float64(x), 6)) ≈ $(join([ "$(c)×$e" for (c, e) in zip(coeffs, elements)], " + "))"
end

function plot_simpleloop(loop, current, vertex_positions,
        elemcurrents, simpleloopscurrent_decomp)
    figure(figsize = (15,15))

    plot_vertices(vertex_positions)
    plot_loop(loop, vertex_positions, color = "C0")

    dstr = decompose_string(current, simpleloopscurrent_decomp, elemcurrents)
    title("""current: $(signif(Float64(current), 6))
             vertices: $loop
             decomposition: $dstr
          """, fontsize = "x-large")
end

function save_plots_simpleloops(folder, simpleloops, simpleloopscurrents,
        vertex_positions, elemcurrents, simpleloopscurrents_decomp)
    @showprogress 1 for i in eachindex(simpleloops)
        plot_simpleloop(simpleloops[i], simpleloopscurrents[i],
            vertex_positions, elemcurrents, simpleloopscurrents_decomp[i])
        savefig(joinpath(folder, "simpleloop$(lpad(i, 3, 0)).png"), dpi = 150)
        close("all")
    end
end

"dirs - what to plot, eg. [3, 1, 2] will put z (the 3rd dimension) on the
    x-axis of the plot, x (the 1st dimension) on the y axis and will keep
    y (the 2nd dimension) fixed.
    spanA is the span of the first of the plotted dimension (along the x axis
    of the resulting plot).
    spanB is the span of the second of the plotted dimension (along the y axis
    of the resulting plot)."
function plot_loops_field(loops, currents, vertex_positions,
        dirs, valC;
        n = 50, Bref = x -> [0, 0, 0],
        spanA = extrema([ p[dirs[1]] for p in vertex_positions]),
        spanB = extrema([ p[dirs[2]] for p in vertex_positions]),
        levels = -5:5)

    As = linspace(spanA..., n)
    Bs = linspace(spanB..., n)

    B = zeros(length(As), length(Bs), 3)
    for iA in eachindex(As), iB in eachindex(Bs)
        p = [As[iA], Bs[iB], valC]
        Bp = field_loops(p, loops, currents, vertex_positions) .- Bref(p)
        B[iA, iB, :] .= Bp
    end

    f, axs = subplots(1, 4, figsize = (10,4), gridspec_kw = Dict("width_ratios" => [1, 1, 1, 0.1]))
    for i in 1:3
        sca(axs[i])
        contourf(As, Bs, B[:,:,i]' * 1e6, levels = levels, extend = "both", cmap="RdBu_r")
        gca().set_aspect("equal")
        xlabel("xyz"[dirs[1:1]])
        ylabel("xyz"[dirs[2:2]])
        dirletter = "xyz"[i]
        title("\$B_$dirletter\$")
        xlim(spanA...)
        ylim(spanB...)
    end

    colorbar(cax = axs[4], label = "field deviation from goal (μT)")
    # tight_layout()
end

function plot_deviation_histogram(poi, loops, loopscurrents, vertex_positions,
        Bgoal)

    deviations = vcat([
        field_loops(p, loops, loopscurrents, vertex_positions) .- Bgoal(p)
        for p in poi ]...)
    plt.hist(deviations * 1e6, bins = 100)
    yscale("log", nonposy="clip")
    xlabel("deviation from goal (μT)")
    ylabel("# POI")
end

function printlnflush(s)
    println(s)
    flush(STDOUT)
end

function save_report(folder, vertex_positions, g, poi, Bgoal, simpleloops,
        simpleloopscurrents, elemcurrents, simpleloopscurrents_decomp;
        extension = "png", dpi = 150,
        levels = -5:5)
    isdir(folder) || mkdir(folder)

    printlnflush("Plotting the geometry of the system...")

    figure(figsize = (15,15))
    plot_system(g, vertex_positions, poi)
    savefig(joinpath(folder, "system.$extension"), dpi = dpi)
    close("all")

    figure(figsize = (15,15))
    plot_system(g, vertex_positions, poi)
    gca().azim = -90
    gca().elev = 90
    savefig(joinpath(folder, "system_top.$extension"), dpi = dpi)
    close("all")

    figure(figsize = (15,15))
    plot_system(g, vertex_positions, poi)
    gca().azim = -90
    gca().elev = 0
    savefig(joinpath(folder, "system_front.$extension"), dpi = dpi)
    close("all")

    figure(figsize = (15,15))
    plot_system(g, vertex_positions, poi)
    gca().azim = 0
    gca().elev = 0
    savefig(joinpath(folder, "system_right.$extension"), dpi = dpi)
    close("all")

    simpleloopscurrents_real = real_decomposed_currents(
        simpleloopscurrents_decomp, elemcurrents)

    printlnflush("Plotting the deviation histogram...")
    figure()
    plot_deviation_histogram(poi, simpleloops, simpleloopscurrents_real,
        vertex_positions, Bgoal)
    savefig(joinpath(folder, "deviation_histogram.$extension"), dpi = dpi)
    close("all")

    printlnflush("Plotting the field of the solution...")

    zspan = extrema([ p[3] for p in vertex_positions])
    for z in linspace(zspan..., 5)[2 : end-1]
        plot_loops_field(simpleloops, simpleloopscurrents_real, vertex_positions,
            [1,2,3], z, n = 50, Bref = Bgoal, levels = levels)
        savefig(joinpath(folder, "field_XY_z$(signif(z,2)).$extension"),
            dpi = dpi)
        close("all")
    end

    xspan = extrema([ p[1] for p in vertex_positions])
    for x in linspace(xspan..., 5)[2 : end-1]
        plot_loops_field(simpleloops, simpleloopscurrents_real, vertex_positions,
            [3,2,1], x, n = 50, Bref = Bgoal, levels = levels)
        savefig(joinpath(folder, "field_ZY_x$(signif(x,2)).$extension"),
            dpi = dpi)
        close("all")
    end

    printlnflush("Plotting the simple loops...")

    save_plots_simpleloops(folder, simpleloops, simpleloopscurrents,
        vertex_positions, elemcurrents, simpleloopscurrents_decomp)
end

end


text:
"""
    plume_rise(scenario::Scenario, plumerise)
Implements the Briggs plume rise equations for buoyancy and momentum driven
plume rise as described in the ISC3 model guide EPA-454/B-95-003b
"""
function plume_rise(scenario, plumerise)

    if plumerise==false
        return x->zero(x)
    end

    # physics parameters
    g = 9.80616 #m/s^2

    required_params = [:jet_diameter, :jet_velocity, :release_temperature,
                       :windspeed, :ambient_temperature, :pasquill_gifford]
    if all(key -> !(ismissing(getproperty(scenario,key))), required_params)

        # parameters of the jet
        Dⱼ = scenario.jet_diameter
        uⱼ = scenario.jet_velocity
        Tᵣ = scenario.release_temperature

        # parameters of the environment
        u = scenario.windspeed
        Tₐ = scenario.ambient_temperature
        stability = scenario.pasquill_gifford

    else
        missing_params = [ String(i) for i in filter(key -> ismissing(getproperty(scenario,key)), required_params)]
        error_string = "These parameters cannot be missing: " * join(missing_params, ", ")
        e = MissingException(error_string)
        throw(e)
    end

    # buoyancy flux
    Fb = g * uⱼ * Dⱼ^2 * (Tᵣ - Tₐ) / (4Tᵣ)

    # momentum flux
    Fm = uⱼ^2 * Dⱼ^2 * Tₐ/(4Tᵣ)

    # stability check
    if stability ∈ Set(["A","B","C","D"])
        stable = false
    elseif stability == "E"
        stable = true
        Γ = 0.020    # default lapse rate K/m
        s = (g/Tₐ)*Γ # stability
    elseif stability == "F"
        stable = true
        Γ = 0.035    # default lapse rate K/m
        s = (g/Tₐ)*Γ # stability
    else
        err = "$stability is not a valid stability class"
        error(err)
    end

    if stable
        ΔTc = 0.019582*Tᵣ*uⱼ*√(s)
        if (Tᵣ - Tₐ) > ΔTc
            # buoyancy dominated plume rise
            final_rise = 2.6*(Fb/(u*s))^(1/3)
            xf = 2.0715*u/√(s)
            return x -> if (x < xf) min(1.60*(Fb*x^2/u^3)^(1/3), final_rise) else final_rise end
        else
            # momentum dominated plume rise
            stable_momentum_rise = 1.5*(Fm/(uⱼ*√(s)))^(1/3)
            unstable_momentum_rise = 3*Dⱼ*(uⱼ/u)
            final_rise = min(stable_momentum_rise, unstable_momentum_rise)
            xf = (π/2)*(u/√(s))
            β = (1/3) + (u/uⱼ)
            return x -> if (x < xf) min((3Fm*sin(x*√(s)/u)/(β^2*u*√(s)))^(1/3), final_rise) else final_rise end
        end
    else
        if Fb < 55
            ΔTc = 0.0297*Tᵣ*(uⱼ/Dⱼ)^(1/3)
            xf = 49*Fb^(5/8)
            buoyant_rise = 21.425*Fb^(3/4)/u
        else
            ΔTc = 0.00575*Tᵣ*(uⱼ^2/Dⱼ)^(1/3)
            xf = 119*Fb^(2/5)
            buoyant_rise = 38.71*Fb^(3/5)/u
        end

        if (Tᵣ - Tₐ) > ΔTc
            # buoyancy dominated plume rise
            final_rise = buoyant_rise
            return x -> if (x < xf) min(1.60*(Fb*x^2/u^3)^(1/3), final_rise) else final_rise end
        else
            # momentum dominated plume rise
            final_rise = 3Dⱼ*(uⱼ/u)
            xf = if (Fb<=0) 4Dⱼ*(uⱼ+3u)^2/(uⱼ*u) else xf end
            β = (1/3) + (u/uⱼ)
            return x -> if (x < xf) min((3Fm*x/(β*u)^2)^(1/3), final_rise) else final_rise end
        end
    end
end

ext: jl


text:
using Revise
using gtpmip
using westbrick

grounded_plan = solve_hpd(joinpath(dirname(pathof(gtpmip)), "..", "experiments/sim/task3/domain.hpd"),
                        joinpath(dirname(pathof(gtpmip)), "..", "experiments/sim/task3/problem.hpd");
                         max_levels=1000)
grounded_plan.plan

num_objects = 9
inits = grounded_plan.graph.initprops[:continuous][1]
xinits = [Symbol("xb"*string(i)) for i=1:num_objects]
yinits = [Symbol("yb"*string(i)) for i=1:num_objects]
init_poses = [[inits[x], inits[y]] for (x,y) in zip(xinits, yinits)]

objects = init_environment(num_objects, init_poses)
obs_dict, ax, fig = visualize_environment!(objects)
bobby = load_robot()
traj, obj_traj = init_visualization()

for i = 1:grounded_plan.size
    pose = grounded_plan.robot_poses[i]
    action = grounded_plan.plan[i]
    move!(bobby, pose, traj, obj_traj)
    if action[1] == :pick
        object_id = parse(Int, string(action[2][1])[end])
        pick!(bobby, -π/2., objects[object_id], traj, obj_traj)
    elseif action[1] == :place
        object_id = string(action[2][1])[end]
        place!(bobby, -π/2., traj, obj_traj)
    end
end

println("Execution time:")
@time visualize_trajectory!(bobby, traj, obj_traj, obs_dict, ax, fig;name="build_i.gif")


text:
add_verbose_scope(:MPolyGcd)

####################################################
# exported is the RecoCtx (and thus the rational_reconstruction functions)
# not exported are the helpers...

module RecoNF

using Hecke
import Nemo

function basis_matrix(d::fmpz, f::fmpz_poly, k::AnticNumberField)
  #assumes f is idl as above!!!
  #1st need to deconstruct f into the different degrees:
  #CRT of degree a>b and implies lead(b) = 0 mod q, hence gcd's are my friend
  #claim: in this situation, the "obvious" method will produce a Howell form
  #tries to compute the basis matrix for the ideal <d, f(a)> where a = gen(k)
  #assumes deg f < deg k, d coprime to the conductor/ index/ everything
  de = []
  g = d
  N = zero_matrix(FlintZZ, degree(k), degree(k))
  dN = fmpz(1)
  res = []
  f_orig = f
  d_orig = d
  for i=degree(f):-1:1
    if degree(f)<i
      continue
    end
    r = Base.gcd(coeff(f, i), g)
    #so I have <d/r, f> of degree i and
    #          <f, f mod r> of smaller degree
    n = div(g, r)
    c = invmod(lead(f), n)
    fn = mod(c*f, n)
    @assert ismonic(fn)
    @assert degree(fn) == i
    if degree(f) == degree(k)
      M = MatrixSpace(FlintZZ, degree(k), degree(k))(n)
    else
      M = zero_matrix(FlintZZ, degree(k), degree(k))
      for j=1:i
        M[j,j] = n
      end
      for j=1:degree(fn)+1
        M[i+1, j] = coeff(fn, j-1)
      end
      t = gen(parent(fn))^i-fn
      for j=i+2:degree(k)
        t = t*gen(parent(fn))
        t -= lead(t)*fn
        mod!(t, n)
        M[j,j] = 1
        for k=1:j-1
          M[j, k] = -coeff(t, k-1)
        end
      end
    end
    if dN == 1
      N = M
      dN = n
    else
      N, dN = induce_crt(N, dN, M, n)
    end
    f = mod(f, r)
    g = r
    if isone(g)
      break
    end
  end
  #TODO: implement the Fieker-Hofmann lifting step to avoid the hnf...
  N = Hecke._hnf_modular_eldiv(N, dN, :lowerleft)
  return N
end

export RecoCtx

mutable struct RecoCtx
  L::fmpz_mat  # should be LLL reduced, will do so on creation
  p1::fmpz     # the "prime": L is the basis matrix of an ideal, p1 is the
               # minimum
  f::fmpz_poly # the implicit ideal is <p1, f(gen(k))>
  LI::fmpz_mat #(Li, d) = pseudo_inv(L) - if set (integral = true)
  d::fmpz
  k::AnticNumberField
  new_data::Bool
  function RecoCtx(A::fmpz_mat, k::AnticNumberField)
    r= new()
    r.k = k
    r.L = lll(A)
    r.p1 = det(A)
    r.new_data = false
    return r
  end
  function RecoCtx(k::AnticNumberField)
    r = new()
    r.L = identity_matrix(FlintZZ, degree(k))
    r.p1 = fmpz(1)
    r.k = k
    r.new_data = false
    return r
  end
end

function Base.push!(R::RecoCtx, p::fmpz, f::fmpz_poly)
  @assert gcd(R.p1, p) == 1

  if R.p1 == 1
    R.f = f
    R.p1 = p
  else
    R.f, R.p1 = induce_crt(R.f, R.p1, f, p)
  end
  R.new_data = true
end

function data_assure(R::RecoCtx)
  R.new_data || return

  R.L = lll(basis_matrix(R.p1, R.f, R.k))
  if isdefined(R, :Li) #to keep stucture consistent
    R.Li, R.d = pseudo_inv(R.L)
  end
  R.new_data = false
  return R
end

function has_small_coeffs(a::nf_elem, B::fmpz)
  z = fmpz()
  for i=0:degree(parent(a))-1
    Nemo.num_coeff!(z, a, i)
    if cmpabs(z, B) >0
      return false
    end
  end
  return true
end

function Hecke.induce_rational_reconstruction(a::Generic.MPoly{nf_elem}, R::RecoCtx; integral::Bool = false)
  b = MPolyBuildCtx(parent(a))
  k = base_ring(a)
  d = k(2)
  if integral
    B = fmpz(1)
  else
    B = abs(det(R.L))
    B = fmpz(2)^div(nbits(B), 2*degree(k))
  end
  for i=1:length(a)
    if integral
      fl, c = rational_reconstruction(coeff(a, i), R, integral = integral)
      if !fl
        return fl, a
      end
    else
      #implicitly assumes elements have a common denominator
      fl, c = rational_reconstruction(coeff(a, i)*d, R, integral = true)
      if !fl || !has_small_coeffs(c, B)
        fl, c, dd = rational_reconstruction(coeff(a, i)*d, R, integral = false, split = true)
        !fl && return fl, a
        (has_small_coeffs(c, B) && has_small_coeffs(d*dd, B)) || return false, a
        d *= dd
      end
      c = c//d
    end
    push_term!(b, c, exponent_vector(a, i))
  end
  return true, finish(b)
end

#TODO: split needs to be a val-arg
function Hecke.rational_reconstruction(a::nf_elem, R::RecoCtx; integral::Bool = false, split::Bool = false)
  data_assure(R)
  if integral
    if !isdefined(R, :LI)
      R.LI, R.d = pseudo_inv(R.L)
    end
    t = zero_matrix(FlintZZ, 1, degree(R.k))
    z = fmpz()
    for i=1:degree(R.k)
      Nemo.num_coeff!(z, a, i-1)
      t[1, i] = z
    end
    s = t*R.LI
    for i=1:degree(R.k)
      s[1, i] = round(s[1, i]//R.d)
    end
    tt = s*R.L
    b = parent(a)()
    nb = div(3*nbits(R.d), 2)
    for i=1:degree(R.k)
      Hecke._num_setcoeff!(b, i-1, t[1, i]-tt[1, i])
      nb -= nbits(t[1, i] - tt[1, i])
    end
    return nb >= 0, b
  end
  n = degree(parent(a))
  Znn = MatrixSpace(FlintZZ, n, n)
  L = [ Znn(1) representation_matrix_q(a)[1] ; Znn(0) R.L]
  lll!(L)
  K = parent(a)
  d = Nemo.elem_from_mat_row(K, sub(L, 1:1, 1:n), 1, fmpz(1))
  n = Nemo.elem_from_mat_row(K, sub(L, 1:1, n+1:2*n), 1, fmpz(1))
  if split
    return true, n, d
  else
    return true, n//d
  end
end

end  #RecoNF module

using .RecoNF

module MPolyGcd

using Hecke
import Nemo, Nemo.nmod_mpoly, Nemo.NmodMPolyRing
import AbstractAlgebra
import Hecke.RecoCtx

function Hecke.gcd(f::Hecke.Generic.MPoly{nf_elem}, g::Hecke.Generic.MPoly{nf_elem})
  Hecke.check_parent(f, g)
  @vprint :MPolyGcd 1 "multivariate gcd of f with $(length(f)) and g with $(length(g)) terms over $(base_ring(f))\n"

  k = base_ring(f)
  ps = PrimesSet(Hecke.p_start, -1)
  fl, c = Hecke.iscyclotomic_type(k)
  if fl
    @vprint :MPolyGcd 2 "field is cyclotomic with conductor $c\n"
    ps = PrimesSet(Hecke.p_start, -1, c, 1)
  end
  fl, c = Hecke.isquadratic_type(k)
  if fl && abs(c) < typemax(Int)
    @vprint :MPolyGcd 2 "field is quadratic, using conductor $(4*c)\n"
    ps = PrimesSet(Hecke.p_start, -1, Int(4*c), 1)
  end
  return _gcd(f, g, ps)
end

function _gcd(f::Hecke.Generic.MPoly{nf_elem}, g::Hecke.Generic.MPoly{nf_elem}, ps::PrimesSet{Int})
#  @show "gcd start"
  p = iterate(ps)[1]
  K = base_ring(f)
  max_stable = 2
  stable = max_stable

  if iszero(f)
    return g
  end
  if iszero(g)
    return f
  end

  # compute deflation and deflate
  shifta, defla = Generic.deflation(f)
  shiftb, deflb = Generic.deflation(g)
  shiftr = min.(shifta, shiftb)
  deflr = broadcast(gcd, defla, deflb)
  f = deflate(f, shifta, deflr)
  g = deflate(g, shiftb, deflr)

  d = fmpz(1)
  gc = parent(f)()
  gd = parent(f)()
  Zx = Hecke.Globals.Zx
  R = RecoCtx(K)

  de = lcm(lcm(map(denominator, coefficients(f))), lcm(map(denominator, coefficients(g))))
  f*=de
  g*=de
  E = equation_order(K)
  lI = E*E(lead(f)) + E*E(lead(g))
  gl = Hecke.short_elem(lI)
  gl *= evaluate(derivative(K.pol), gen(K))  # use Kronnecker basis

  fl = true
  while true
    p = iterate(ps, p)[1]
    @vprint :MPolyGcd 2 "Main loop: using $p\n"
    @vtime :MPolyGcd 3 me = Hecke.modular_init(K, p, deg_limit = 1)
    if isempty(me)
      continue
    end

    @vtime :MPolyGcd 3 fp = Hecke.modular_proj(f, me)
    @vtime :MPolyGcd 3 gp = Hecke.modular_proj(g, me)
    glp = Hecke.modular_proj(gl, me)
    gcd_p = nmod_mpoly[]
    @vtime :MPolyGcd 3 for i=1:length(fp)
      _g = gcd(fp[i], gp[i])
      if length(_g) == 1 && iszero(exponent_vector(_g, 1))
        return inflate(one(parent(f)), shiftr, deflr)
      end
      push!(gcd_p, coeff(glp[i], 0)*_g)
    end
    #gcd_p = [coeff(glp[i], 0)*gcd(fp[i], gp[i]) for i=1:length(fp)]
    @vtime :MPolyGcd 3 tp = Hecke.modular_lift(gcd_p, me)
    if d==1
      d = fmpz(p)
      gc = tp
      idl = lift(Zx, me.ce.pr[end])
      push!(R, d, idl)
      fl, gd = induce_rational_reconstruction(gc, R, integral = true)
      if fl && divides(f, gd)[1] && divides(g, gd)[1]
#          @show "gcd stop", nbits(d), length(gd), gd
#          @time fl, q = divides(f, gd)
#          @time q = div(f, gd)
#          @time q*gd == f
          gd*=inv(gl)
          @assert isone(lead(gd))
          return inflate(gd, shiftr, deflr)
      end
      stable = max_stable
    else
      #TODO: instead of lifting "idl" and doing basis_matrix from
      #      scratch, do the basis_matrix for the new ideal and
      #      use CRT to combine them
      #TODO: explore combining LLL matrices to speed up LLL....
#TODO: deal with bad primes...

      push!(R, fmpz(p), lift(Zx, me.ce.pr[end]))
      if (!fl) || any(i->(parent(me.ce.pr[end])(coeff(tp, i) - coeff(gd, i))) % me.ce.pr[end] != 0, 1:length(tp))
        gc, d = induce_crt(gc, d, tp, fmpz(p), true)
        fl, gd = induce_rational_reconstruction(gc, R, integral = true)
        stable = max_stable
      else
        d *= p
        stable -= 1
      end
        if true || stable <= 0
          if divides(f, gd)[1] && divides(g, gd)[1]
#            @show "gcd stop", nbits(d), length(gd), gd
            gd*=inv(gl)
            @assert isone(lead(gd))
            return inflate(gd, shiftr, deflr)
          else
            stable = max_stable
          end
        end
      #before I forget: gc is "the" gcd modulo <d, idl>
    end
  end
end

function Hecke.induce_crt(a::Hecke.Generic.MPoly{nf_elem}, p::fmpz, b::Hecke.Generic.MPoly{nf_elem}, q::fmpz, signed::Bool = false)
  pi = invmod(p, q)
  mul!(pi, pi, p)
  pq = p*q
  if signed
    pq2 = div(pq, 2)
  else
    pq2 = fmpz(0)
  end
  z = zero(base_ring(a))

  #=
  c = (b-a)*pi+a
  mod!(c, pq)
  return c
  =#

  N = ngens(parent(a))

  ta = terms(a)
  tb = terms(b)
  c = MPolyBuildCtx(parent(a))
  aa, sa = iterate(ta)
  bb, sb = iterate(tb)
#  @assert length(a) == length(b)
#  @assert ==(aa, bb, true) # leading terms must agree or else...
  while !(aa === nothing) && !(bb === nothing)
    if ==(aa.exps, bb.exps) #monomial equality
      push_term!(c, Hecke.induce_inner_crt(coeff(aa, 1), coeff(bb, 1), pi, pq, pq2), exponent_vector(aa, 1))
      aa = iterate(ta, sa)
      bb = iterate(tb, sb)
      aa === nothing && break
      aa, sa = aa
      bb === nothing && break
      bb, sb = bb
    elseif Generic.monomial_isless(aa.exps, 1, bb.exps, 1, N, parent(aa), UInt(0)) #aa < bb
      push_term!(c, Hecke.induce_inner_crt(z, coeff(bb, 1), pi, pq, pq2), exponent_vector(bb, 1))
      bb, sb = iterate(tb, sb)
      bb === nothing && break
      bb, sb = bb
    else
      push_term!(c, Hecke.induce_inner_crt(coeff(aa, 1), z, pi, pq, pq2), exponent_vector(aa, 1))
      aa = iterate(ta, sa)
      aa === nothing && break
      aa, sa = aa
    end
  end
  while !(aa === nothing)
    push_term!(c, Hecke.induce_inner_crt(coeff(aa, 1), z, pi, pq, pq2), exponent_vector(aa, 1))
    aa = iterate(ta, sa)
    if !aa == nothing
      aa, sa = aa
    end
  end
  while !(bb === nothing)
    push_term!(c, Hecke.induce_inner_crt(z, coeff(bb, 1), pi, pq, pq2), exponent_vector(bb, 1))
    bb = iterate(tb, sb)
    if !(bb === nothing)
      bb, sb = bb
    end
  end
  return finish(c), pq
end

function Hecke.induce_crt(a::fmpz_mat, p::fmpz, b::fmpz_mat, q::fmpz, signed::Bool = false)
  pi = invmod(p, q)
  mul!(pi, pi, p)
  pq = p*q
  if signed
    pq2 = div(pq, 2)
  else
    pq2 = fmpz(0)
  end

  @assert size(a) == size(b)
  c = similar(a)
  for i=1:nrows(a)
    for j=1:ncols(a)
      c[i,j] = Hecke.inner_crt(a[i,j], b[i,j], pi, pq, pq2)
    end
  end
  return c, pq
end

function Hecke.modular_proj(f::Generic.MPoly{nf_elem}, me::Hecke.modular_env)
  if !isdefined(me, :Kxy)
    me.Kxy = parent(f)
  else
    @assert me.Kxy === parent(f)
  end
  if !isdefined(me, :Kpxy)
    p = characteristic(me.Fpx)
    me.Kpxy, _ = PolynomialRing(base_ring(me.Fpx), ["$(x)_$p" for x = me.Kxy.S])
  end
  fp = [MPolyBuildCtx(me.Kpxy) for x = me.fld]
  s = length(me.fld)
  for i=1:length(f)
    c = coeff(f, i)
    e = exponent_vector(f, i)
    cp = Hecke.modular_proj(c, me)
    R = base_ring(me.Fpx)
    for x = 1:s
      push_term!(fp[x], Hecke.nmod(coeff(cp[x], 0), R), e)
    end
  end
  return map(finish, fp)
end

function Hecke.modular_lift(g::Array{nmod_mpoly, 1}, me::Hecke.modular_env)

  #TODO: no dict, but do s.th. similar to induce_crt
  d = Dict{Array{Int, 1}, Array{Tuple{Int, Hecke.nmod}, 1}}()
  for i=1:length(g)
    for (c, e) = Base.Iterators.zip(Generic.MPolyCoeffs(g[i]), Generic.MPolyExponentVectors(g[i]))
      if Base.haskey(d, e)
        push!(d[e], (i, c))
      else
        d[e] = [(i, c)]
      end
    end
  end
  bt = MPolyBuildCtx(me.Kxy)

  for e = keys(d)
    for x=1:length(g)
      me.res[x] = zero!(me.res[x])
    end
    for (i, c) = d[e]
      me.res[i] = parent(me.res[i])(lift(c))
    end
    push_term!(bt, Hecke.modular_lift(me.res, me), e)
  end
  return finish(bt)

  bt = MPolyBuildCtx(me.Kxy)
  #TODO deal with different vectors properly (check induce_crt)
  @assert all(x->collect(exponent_vectors(g[1])) == collect(exponent_vectors(g[x])), 2:length(g))
  for i=1:length(g[1])
    for x=1:length(g)
      me.res[x] = parent(me.res[x])(lift(coeff(g[x], i)))
    end
    push_term!(bt, Hecke.modular_lift(me.res, me), exponent_vector(g[1], i))
  end
  return finish(bt)
end

function Hecke.mod!(f::fmpz_poly, p::fmpz)
  for i=0:degree(f)
    setcoeff!(f, i, mod(coeff(f, i), p))
  end
end

function Hecke.mod(f::fmpz_poly, p::fmpz)
  g = parent(f)()
  for i=0:degree(f)
    setcoeff!(g, i, mod(coeff(f, i), p))
  end
  return g
end

function Hecke.mod_sym!(f::fmpz_poly, p::fmpz)
  for i=0:degree(f)
    setcoeff!(f, i, Hecke.mod_sym(coeff(f, i), p))
  end
end

#=
import Base.//, Base.==

struct Term{T}
  f::T
  i::Int
  function Term(f::T, i::Int) where {T <: AbstractAlgebra.MPolyElem}
    return new{T}(f, i)
  end
end

function Base.show(io::IO, t::Term)
  print(io, "$(t.i)-th term of $(t.f)")
end

struct Terms{T}
  f::T
  function Terms(f::T) where {T <: AbstractAlgebra.MPolyElem}
    return new{T}(f)
  end
end

function Base.show(io::IO, t::Terms)
  print(io, "Iterator for the terms of $(t.f)")
end

function Base.iterate(T::Terms, st::Int = 0)
  st += 1
  if st > length(T.f)
    return nothing
  end
  return Term(T.f, st), st
end

Base.IteratorEltype(M::Terms) = Base.HasEltype()
Base.eltype(M::Terms{T}) where {T} = Term{T}

Base.IteratorSize(M::Terms) = Base.HasLength()
Base.length(M::Terms) = length(M.f)

function Base.lastindex(a::Terms)
  return length(a.f)
end

function Base.getindex(a::Terms, i::Int)
  return Term(a.f, i)
end

function Base.isless(f::Term, g::Term)
  R = parent(f.f)
  @assert R == parent(g.f)
  return AbstractAlgebra.Generic.monomial_isless(f.f.exps, f.i, g.f.exps, g.i, ngens(R), R, UInt(0))
end

function ==(f::Term, g::Term, monomial_only::Bool = false)
  R = parent(f.f)
  @assert R == parent(g.f)

  return AbstractAlgebra.Generic.monomial_cmp(f.f.exps, f.i, g.f.exps, g.i, ngens(R), R, UInt(0))==0 && (monomial_only || coeff(f.f, f.i) == coeff(g.f, g.i))
end

#=
function push_term!(M::MPolyBuildCtx{<:Generic.MPoly{T}}, t::Term{T}) where {T}
  push_term!(M, coeff(t.f, t.f.i), exponent_vector(t.f, t.f.i))
end
=#

function Hecke.coeff(t::Term)
  return coeff(t.f, t.i)
end

function Hecke.exponent_vector(t::Term)
  return exponent_vector(t.f, t.i)
end

function monomial(t::Term)
  m = parent(r.f)()
  set_exponent_vector!(m, 1, exponent_vector(t))
  setcoeff!(m, one(base_ring(m)))
  return m
end

function lead_term(f::AbstractAlgebra.MPolyElem)
  return Term(f, 1)
end

=#
#=TODO
  fit! for nmod_mpoly
  coeff(fq_nmod) -> UInt (should be nmod)
  nmod_mpoly -> gfp_mpoly? at least in Nemo
  set_coeff should accept UInt

  deal with bad primes (wrong expo vectors)
  reconstruction - and use it in the _hensel stuff elsewhere...
  deal with content

=#

end


text:
module Pda
using StaticArrays
import HDF5
using CSV
using DataFrames
using Statistics
using LinearAlgebra
import Polynomials: fit
"""
    polyfit(x,y,n)
polynomial fitting.
"""
polyfit(x,y,n)=fit(x,y,n)
"""
    linearfit(x,y)
Fit a straight line.
"""
linearfit(x,y)=polyfit(x,y,1)

include("io.jl")
include("coordinates.jl")
include("geometry.jl")

include("density.jl")
include("polarorder.jl")

end


text:
using Piecewise
using Test
using Measurements

include("algebraic_datatypes.jl")

@testset "Static polynomials" begin
    include("staticpolynomials.jl")
end

@testset "Piecewise functions" begin

#@testset "Piecewise polynomials" begin
polynomials = [
    StaticPolynomial(0.0, 0.0, 0.0),
    p"1.0-x^2",
    p"x^2 - 1.0",
    StaticPolynomial(0.0, 0.0, 0.0),
]

functions = [
    x -> 0.0,
    x -> 1.0 - x^2,
    x -> x^2 - 1.0,
    x -> 0.0
]

breakpoints = [-1.0, 0.0, 1.0]

mypiecewise = PiecewisePolynomial{3}(polynomials, breakpoints)
mypiecewise_funcs = OrderedPiecewiseFunction{3}(functions, breakpoints)

differentiated_polys = [zero(StaticPolynomial{Float64, 2}), p"-2.0x", p"2.0x", zero(StaticPolynomial{Float64, 2})]
integrated_polys = [zero(StaticPolynomial{Float64, 4}), p"x - x^3/3", p"x^3/3 - x", zero(StaticPolynomial{Float64, 4})]

mypiecewise_differentiated = PiecewisePolynomial{3}(differentiated_polys, breakpoints)
mypiecewise_integrated = PiecewisePolynomial{3}(integrated_polys, breakpoints)

@test mypiecewise_differentiated == differentiate(mypiecewise)
@test mypiecewise_integrated == integrate(mypiecewise)

@test mypiecewise(1.1) == mypiecewise(-1.1) == mypiecewise(-1) == 0 == mypiecewise_funcs(1.1) == mypiecewise_funcs(-1.1) == mypiecewise(-1)

op1 = Piecewise.@ordered_piecewise begin
    1 => x -> 0
    _ => x -> 1
end

op2 = Piecewise.@ordered_piecewise [x -> 0, x -> 1] [1]

pp1 = Piecewise.@piecewise_polynomial begin
    1 => p"0"
    _ => p"1"
end

pp2 = Piecewise.@piecewise_polynomial [p"0", p"1"] [1]

@test op1(1) == pp1(1) == op2(1) == pp2(1) == 1
@test op1(0.99) == pp1(0.99) == pp2(0.99) == 0

pp3 = Piecewise.@piecewise_polynomial begin
    1 => p"0"
    _ => p"x"
end

@test all(length.(pp3.polynomials) .== 2) && pp3.polynomials isa Vector{StaticPolynomial{Int64, 2}}

op_macro = Piecewise.@ordered_piecewise begin
    -1.0 => x -> 0.0
    0.0 => x -> 1.0 - x^2
    1.0 => x -> x^2- 1.0 - cos(x^2)
    _ => x -> 0.0
end

@test op_macro.expressions == ["x -> 0.0", "x -> 1.0 - x^2", "x -> (x^2 - 1.0) - cos(x^2)", "x -> 0.0"]

@test p"6*(2 - q)" == p"12 - 6q"
@test p"(2 - q) + (3 - q)*3" == p"11 - 4q"

_M6_1 =  @piecewise_polynomial begin
    0.0 => p"0.0"
    1.0 => p"(3 - q)^5 - 6*(2 - q)^5 + 15*(1 - q)^5"
    2.0 => p"(3 - q)^5 - 6*(2 - q)^5"
    3.0 => p"(3 - q)^5"
    _ => p"0.0"
end

# TODO: Write tests for printing

end

nothing


text:
module Metatheory

using RuntimeGeneratedFunctions
using Base.Meta

include("docstrings.jl")

RuntimeGeneratedFunctions.init(@__MODULE__)

# TODO document options
# Options
options = Dict{Symbol, Any}(
    :verbose => false,
    :printiter => false,
)

macro log(args...)
    quote options[:verbose] && @info($(args...)) end |> esc
end

export options

include("Util/Util.jl")
using .Util
export Util

# TODO document this interface
include("TermInterface.jl")
using .TermInterface
export TermInterface

include("rgf.jl")
export @metatheory_init

include("rule.jl")
export Rule
export @rule
export RHS_FUNCTION_CACHE
export getrhsfun

include("theory.jl")
export @theory
export Theory

include("Classic/Classic.jl")
using .Classic: gettheory
export Classic

include("EGraphs/EGraphs.jl")
export EGraphs

include("Library/Library.jl")
export Library

end # module


text:
using SurfaceTopology
using Documenter
using Literate

Literate.markdown(joinpath(@__DIR__, "../examples/features.jl"), joinpath(@__DIR__,"src/"); credit = false, name = "index")

makedocs(sitename="SurfaceTopology.jl",pages = ["index.md"])

deploydocs(
     repo = "github.com/akels/SurfaceTopology.jl.git",
 )


text:
@testset "AnnulusView" begin
    data = ones(10, 501, 501)

    av = AnnulusView(data; inner=10, outer=30)

    @test size(av) == size(data)
    @test 10 * π * (av.rmax^2 - av.rmin^2) ≈ sum(av) rtol=1/sqrt(sum(av))

    X = av()
    @test all(size(X) .== map(length, av.indices))
    @test all(isone, X)
    sub_cube = inverse(av, -ones(size(X)))
    @test sub_cube ≈ copyto!(copy(av), -ones(size(X)))

    X_view = av(true)
    @test X_view isa SubArray
    @test X_view ≈ X

    @test -10π * (av.rmax^2 - av.rmin^2) ≈ sum(sub_cube) rtol=1/sqrt(count(!iszero, sub_cube))

    data = randn(rng, 10, 501, 501)

    @test AnnulusView(data; inner=0, outer=Inf) ≈ data

    nanview = AnnulusView(data; fill=NaN)
    @test_throws BoundsError nanview[-1, 2, 3]
    @test isnan(nanview[1, 1, 1])
end

@testset "MultiAnnulusView" begin
    data = ones(10, 501, 501)
    fwhm = 5
    av = MultiAnnulusView(data, fwhm; inner=10, outer=30)

    @test size(av) == size(data)
    @test av.radii ≈ 12.5:5:27.5
    @test length(av.indices) == length(av.radii)
    @test av.width ≈ fwhm
    area = 10 * π * (30^2 - 10^2)
    @test area ≈ sum(av) rtol=1/sqrt(count(!iszero, av))

    X = av(1)
    @test all(size(X) .== map(length, av.indices[1]))
    @test all(isone, X)
    # only subs the first annulus
    @test_throws MethodError inverse(av, -X)
    @test_throws ErrorException copyto!(copy(av), -X)
    sub_cube = inverse(av, 1, -X)
    @test all(sub_cube[av.indices[1]...] .≈ -1)
    copy_cube = copyto!(copy(av), 1, -X)
    @test all(copy_cube[av.indices[1]...] .≈ -1)

    @test_throws DimensionMismatch inverse(av, 2, X)

    X_view = av(1, true)
    @test X_view isa SubArray
    @test X_view ≈ X

    av = MultiAnnulusView(data, fwhm; inner=10, outer=30)
    @test eachannulus(av) isa Base.Generator
    Xs = [-X for X in eachannulus(av)]
    sub_cube = inverse(av, Xs)
    @test sub_cube ≈ copyto!(copy(av), Xs)

    @test -area ≈ sum(sub_cube) rtol=1/sqrt(count(!iszero, sub_cube))

    data = randn(rng, 10, 501, 501)

    @test MultiAnnulusView(data, fwhm; inner=0, outer=Inf) ≈ data

    nanview = MultiAnnulusView(data, fwhm; fill=NaN)
    @test_throws BoundsError nanview[-1, 2, 3]
    @test isnan(nanview[1, 1, 1])
end


text:
@testset "NfAbsNS" begin

  Qx, x = FlintQQ["x"]
  K, (a, b) = @inferred number_field([x^2 - 2, x^3 - 3])

  @test K isa NfAbsNS
  @test elem_type(NfAbsNS) == NfAbsNSElem
  @test parent_type(NfAbsNSElem) == NfAbsNS
  @test show_minus_one(NfAbsNSElem) isa Bool
  @test !issimple(K)
  @test !issimple(NfAbsNS)

  @testset "Basics" begin
    @test FlintQQ == @inferred base_ring(K)
    @test 6 == @inferred degree(K)
    @test [2, 3] == @inferred degrees(K)
    @test 2 == @inferred ngens(K)
    @test [a, b] == @inferred gens(K)
    @test [one(K), a, b, a * b, b^2, a * b^2] == @inferred basis(K)
    @test string(K) isa String
    @test string(a) isa String
  end

  @testset "Elements" begin
    c = @inferred deepcopy(a)
    @test parent(c) === K
    d = @inferred rand(K, -10:10)
    @test parent(d) === K
    e = @inferred zero(K)
    @test @inferred iszero(e)
    e = @inferred one(K)
    @test @inferred isone(e)

    for i in 1:10
      z = @inferred rand(K, -10:10)
      @test z == @inferred K(data(z))
    end

    for i in 1:10
      z = rand(K, -10:10)

      @test @inferred ==(z, z)
      @test !(@inferred ==(z, z + 1))

      zm = @inferred -z
      @test iszero(zm + z)

      y = rand(K, -10:10)
      w = rand(K, -10:10)
      @test z * y + w * y == @inferred (z + w) *y

      while iszero(z)
        z = rand(K, -10:10)
      end
      zinv = @inferred inv(z)
      @test isone(@inferred z * zinv)

      for T in [Int, BigInt, fmpq, fmpz]
        u = rand(-10:10)
        while iszero(u)
          u = rand(-10:10)
        end
        @test y + u == y + u
        @test y - u == -(u - y)
        @test y * u == u * y
        @test y//u == inv(K(u)//y)
      end

      u = rand(-10:10)
      @test @inferred z^u == @inferred z^fmpz(u)

      M = zero_matrix(FlintQQ, 1, 6)
      Hecke.elem_to_mat_row!(M, 1, z)
      @test z == sum(M[1, j] * basis(K)[j] for j in 1:6)
      M = matrix(FlintQQ, 1, 6, [rand(-10:10) for j in 1:6])
      zz = Hecke.elem_from_mat_row(K, M, 1)
      @test zz == sum(M[1, j] * basis(K)[j] for j in 1:6)

      f = @inferred minpoly(z)
      @test length(factor(f)) == 1 && first(factor(f))[2] == 1
      @test ismonic(f)
      @test iszero(f(z))

      f = @inferred Hecke.minpoly_sparse(z)
      @test length(factor(f)) == 1 && first(factor(f))[2] == 1
      @test ismonic(f)
      @test iszero(f(z))

      f = @inferred Hecke.minpoly_dense(z)
      @test length(factor(f)) == 1 && first(factor(f))[2] == 1
      @test ismonic(f)
      @test iszero(f(z))

      f = @inferred charpoly(z)
      @test length(factor(f)) == 1
      @test ismonic(f)
      @test iszero(f(z))
      @test degree(f) == 6

      M = @inferred representation_matrix(z)
      for i in 1:nrows(M)
        @test Hecke.elem_from_mat_row(K, M, i) == z * basis(K)[i]
      end

      n = @inferred norm(z)
      @test n == det(M)
      t = @inferred tr(z)
      @test t == tr(M)
    end
  end

  @testset "Morphisms" begin
    L, m = @inferred simple_extension(K)
    LL,  = number_field(x^6-6*x^4-6*x^3+12*x^2-36*x+1)
    @test isisomorphic(L, LL)[1]
    for i in 1:100
      z = rand(L, -10:10)
      zz = rand(L, -10:10)
      while iszero(z)
        z = rand(L, -10:10)
      end

      @test @inferred m(z + zz) == m(z) + m(zz)
      @test @inferred m(z * zz) == m(z) * m(zz)
      @test @inferred m(-z) == -m(z)
      @test @inferred m(inv(z)) == inv(m(z))
      @test m\m(z) == z
    end

    for i in 1:100
      z = rand(K, -10:10)
      zz = rand(K, -10:10)
      while iszero(z)
        z = rand(K, -10:10)
      end

      @test @inferred m\(z + zz) == m\(z) + m\(zz)
      @test @inferred m\(z * zz) == m\(z) * (m\(zz))
      @test @inferred m\(-z) == -(m\(z))
      @test @inferred m\(inv(z)) == inv(m\z)
      @test m(m\z) == z
    end
  end

  @testset "Maximal order" begin
    K2,  = @inferred number_field([x^2 - 50, x^3 - 3])
    O = EquationOrder(K2)
    Omax = @inferred MaximalOrder(O)
    @test discriminant(Omax) == FlintZZ(30233088)
  end
end


text:
export UCEvo, UCInd, evaluate, populate, generation, mean_fitness

using Cambrian
using Statistics

import Cambrian.populate, Cambrian.evaluate, Cambrian.log_gen, Cambrian.save_gen

"""
Individuals definition
"""
mutable struct UCInd <: Cambrian.Individual
    e_chromosome::Vector{Float64}
    c_chromosome::Vector{Float64}
    fitnesses::Vector{Float64}
    lifetime::Int64
    reached_frames::Int64
end

function UCInd(e_chromosome::Vector{Float64}, c_chromosome::Vector{Float64})
    UCInd(e_chromosome, c_chromosome, Vector{Float64}(), 0, 0)
end

function UCInd(
	config::NamedTuple,
	e_chromosome::Vector{Float64},
	c_chromosome::Vector{Float64}
)
    UCInd(e_chromosome, c_chromosome)
end

"""
Evolution framework
"""
mutable struct UCEvo{T} <: Cambrian.AbstractEvolution
    config::NamedTuple
    logid::String
    logger::CambrianLogger
	resdir::String
    population::Vector{T}
    fitness::Function
    gen::Int64
    atari_games::Vector{Game}
end

function UCEvo{T}(
	config::NamedTuple,
    resdir::String,
	fitness::Function,
    init_population::Function,
    rom_name::String
) where T
    logid = config.id
    log_path = joinpath(resdir, logid, "logs/logs.csv")
    logger = CambrianLogger(log_path)
    population = init_population(T, config)
	# n_children = config.n_population - config.n_elite
	atari_games = [Game(rom_name, 0) for _ in 1:length(population)]
    UCEvo(config, logid, logger, resdir, population, fitness, 0, atari_games)
end

evaluate(e::UCEvo{T}) where T = fitness_evaluate(e, e.fitness)
populate(e::UCEvo{T}) where T = tournament_populate(e)

ucb(m::Float64, n_total::Int64, n_eval::Int64) = m + sqrt(2.0 * log(n_total) / n_eval)

function ucb(ind::UCInd)
	n_eval = length(ind.fitnesses)
	if n_eval > 0
		return ucb(mean(ind.fitnesses), ind.lifetime, n_eval)
	else
		return Inf
	end
end

function mean_fitness(ind::UCInd)
	if length(ind.fitnesses) > 0
		return mean(ind.fitnesses)
	else
		return -Inf
	end
end

function increment_lifetime!(e::UCEvo{T}, n::Int64) where T
	for ind in e.population
		ind.lifetime += n
	end
end

"""
Elite selection function
"""
function generation(e::UCEvo{T}) where T
    if e.config.greedy_elitism
        sort!(e.population, by=ind->mean_fitness(ind), rev=true)
    else
        sort!(e.population, by=ind->ucb(ind), rev=true)
    end
    e.population = e.population[1:e.config.n_elite]
end

"""
Evaluation function
"""
function fitness_evaluate(e::UCEvo{T}, fitness::Function) where T
	if e.config.surrogate
		if e.gen == 1
			exact_fitness_evaluate(e, fitness)
		else
			surrogate_fitness_evaluate(e, fitness)
		end
	else
		exact_fitness_evaluate(e, fitness)
	end
end

"""
An individual is not evaluated if enough evaluations have been made.
"""
function evaluation_utility(n_eval_max::Int64, ind::UCInd) where T
	if n_eval_max > 0
		if length(ind.fitnesses) < n_eval_max
			ucb(ind)
		else
			-Inf
		end
	else
		ucb(ind)
	end
end

function evaluation_utility(
	n_eval_max::Int64,
	m::Float64,
	n_total::Int64,
	n_eval::Int64
) where T
	if n_eval_max > 0
		if n_eval < n_eval_max
			ucb(m, n_total, n_eval)
		else
			-Inf
		end
	else
		ucb(m, n_total, n_eval)
	end
end

function exact_fitness_evaluate(e::UCEvo{T}, fitness::Function) where T
	# 1. Evaluate all new individuals in parallel
	n_children = e.config.n_population - e.config.n_elite
	sort!(e.population, by=ind->ucb(ind), rev=true) # Children are 1st
	n = e.gen == 1 ? length(e.population) : n_children
	#=
	for i in 1:n
		push!(e.population[i].fitnesses,
			  fitness(e.population[i], e.gen, e.atari_games[i]))
	end
	=#
	lck = ReentrantLock()
	@sync for i in 1:n
        Threads.@spawn begin
  			f = fitness(e.population[i], e.gen, e.atari_games[i])
  			lock(lck)
  			try
  				push!(e.population[i].fitnesses, f)
  			finally
  				unlock(lck)
  			end
        end
    end
	increment_lifetime!(e, n)
	@assert all([length(ind.fitnesses) > 0 for ind in e.population])
	# 2. Finish evaluation budget sequentially
	n_remaining_eval = e.config.n_eval - n
	for _ in 1:n_remaining_eval
		sort!(e.population, by=ind->evaluation_utility(e.config.n_eval_max, ind), rev=true)
		push!(e.population[1].fitnesses,
			  fitness(e.population[1], e.gen, e.atari_games[1]))
		increment_lifetime!(e, 1)
	end
end

function surrogate_fitness_evaluate(e::UCEvo{T}, fitness::Function) where T
	# 1. Evaluate all new individuals in parallel
	n_children = e.config.n_population - e.config.n_elite
	sort!(e.population, by=ind->ucb(ind), rev=true) # Children are 1st
	lck = ReentrantLock()
	@sync for i in 1:n_children
        Threads.@spawn begin
			# OUTDATED unsafe (I thought it was safe)
			#=push!(e.population[i].fitnesses,
				fitness(e.population[i], e.gen, e.atari_games[i]))=#
			# NEW safe
			f = fitness(e.population[i], e.gen, e.atari_games[i])
			lock(lck)
			try
				push!(e.population[i].fitnesses, f)
			finally
				unlock(lck)
			end
        end
    end
	increment_lifetime!(e, n_children)
	@assert all([length(ind.fitnesses) > 0 for ind in e.population])
	# 2. Determine individuals to evaluate next
	n_remaining_eval = e.config.n_eval - n_children
	f = [mean_fitness(ind) for ind in e.population]
	N = [ind.lifetime for ind in e.population]
	n = [length(ind.fitnesses) for ind in e.population]
	to_eval = Vector{Int64}()
	for i in 1:n_remaining_eval
		utilities = [evaluation_utility(e.config.n_eval_max, f[i], N[i], n[i]) for i in eachindex(f)]
		index = argmax(utilities)
		push!(to_eval, index)
		N .+= 1
		n[index] += 1
	end
	# 3. Evaluate remaining individuals
	lck = ReentrantLock()
	@sync for k in eachindex(to_eval)
        Threads.@spawn begin
			i = to_eval[k]
			# OUTDATED unsafe
			#=push!(e.population[i].fitnesses,
				fitness(e.population[i], e.gen, e.atari_games[i]))=#
			# NEW safe
			f = fitness(e.population[i], e.gen, e.atari_games[i])
			lock(lck)
			try
				push!(e.population[i].fitnesses, f)
			finally
				unlock(lck)
			end
		end
    end
	increment_lifetime!(e, n_remaining_eval)
end

"""
Populate function
"""
function tournament_populate(e::UCEvo{T}) where T
    children = Vector{T}()
    for _ in e.config.n_elite+1:e.config.n_population
        parent = ucb_tournament_selection(e.population, e.config.tournament_size)
        push!(children, mutate(parent))
    end
    push!(e.population, children...)
end

function ucb_tournament_selection(
    pop::Vector{T},
    tournament_size::Int64
) where T
    inds = shuffle!(collect(1:length(pop)))
    sort(pop[inds[1:tournament_size]], by=ind->ucb(ind))[end]
end

function fitness_tournament_selection(
    pop::Vector{T},
    tournament_size::Int64
) where T
    inds = shuffle!(collect(1:length(pop)))
    sort(pop[inds[1:tournament_size]], by=ind->ind.fitness)[end]
end

"""
Logging
"""
function log_gen(e::UCEvo{T}) where T
	sep = ";"
    if e.gen == 1
        f = open(joinpath(e.resdir, e.logid, "logs/logs.csv"), "w+")
        write(f, string("gen_number", sep, "fitnesses", sep,
			"reached_frames", sep, "dna_id\n"))
        close(f)
    end
	enco_path = joinpath(e.resdir, e.logid, Formatting.format("gens/encoder_{1:04d}", e.gen))
    cont_path = joinpath(e.resdir, e.logid, Formatting.format("gens/controller_{1:04d}", e.gen))
    mkpath(enco_path)
    mkpath(cont_path)
    for i in eachindex(e.population)
        dna_id = Formatting.format("{1:04d}", i)
		# Log results
		f = open(joinpath(e.resdir, e.logid, "logs/logs.csv"), "a+")
        write(f, Formatting.format(string(e.gen, sep,
			   e.population[i].fitnesses, sep,
			   e.population[i].reached_frames, sep, dna_id, "\n")))
        close(f)
		# Log individuals
		ind_fit = [mean_fitness(e.population[i])]
		enco = IPCGPInd(e.config.e_config, e.population[i].e_chromosome)
		enco.fitness .= ind_fit
		f = open(string(enco_path, "/", dna_id, ".dna"), "w+")
        write(f, string(enco))
        close(f)
        cont = CGPInd(e.config.c_config, e.population[i].c_chromosome)
		cont.fitness .= ind_fit
        f = open(string(cont_path, "/", dna_id, ".dna"), "w+")
        write(f, string(cont))
        close(f)
    end
end


text:
using ProjAssistant
@quickactivate

# ------------------------------------------------------
@time begin
    import Chemostat_InSilico
    const Dyn = Chemostat_InSilico.Dynamic
    import Chemostat_InSilico.Dynamic:
        run_simD3!, SimD3, hist,
        n_ave_conv!, tail_ave,
        Container, vec!, Vi,
        normalizeP!, batch_file,
        set_status, get_status,
        save_batch, save_bfiles, save_simdat,
        UNDONE_SIM_STATUS, DEAD_SIM_STATUS,
        EXPLODED_SIM_STATUS, NITERS_SIM_STATUS,
        STST_SIM_STATUS, reset!,
        STST, is_steady

    import SimTools
    import SimTools: grad_desc, gd_value

    using Base.Threads

end

## ------------------------------------------------------
function run_sim3D(simid, Ds, ϵs, cg;
        dosave_batches = true,
        dosave_simdat = true,
        docollect_Ps = true
    )

    iter = Iterators.product(Ds, ϵs, cg)
    tot_sims = length(iter)
    Ch = Channel() do Ch_
        for (simi, (D, ϵ, cg)) in enumerate(iter)
            put!(Ch_, (simi, (D, ϵ, cg)))
        end
    end

    # Space
    V = lglob(Dyn, :Vcell3D)
    Vz = Vi(V, :z)
    Vug = Vi(V, :ug)
    Vuo = Vi(V, :uo)

    @threads for _ in 1:nthreads()
        thid = threadid()

        for (simi, (D, ϵ, cg)) in Ch

            # ---------------------------------------------------------------
            # sim globals
            # SimD3
            S = SimD3(;
                # Space
                V,
                # Chemostat
                D, ϵ, cg,
                X = 0.5,
                sg = cg,
                Δt = 0.05,
                niters = 20000
            )

            P = S.P

            simparams = (;S.D, S.ϵ, S.cg)

            # batchs
            batch_files = []

            # Welcome & status
            status = get_status(simid, simparams)
            @info("At", simid, simparams, thid, status)

            (status != UNDONE_SIM_STATUS) && continue

            # exploded th
            z_tail_acc = 0.0
            z_tail_count = 0
            z_tail_len = 20
            exploded_Xth = 1e2
            dead_Xth = 5e-2
            z_D_diff_th = 0.05

            # time series
            push_frec = 10
            save_frec = 100 * push_frec

            Xts = Float64[]
            sgts = Float64[]
            z_avts = Float64[]
            ug_avts = Float64[]
            uo_avts = Float64[]
            cgD_Xts = Float64[]

            Pzts = Dict{Float64, Float64}[]
            Pugts = Dict{Float64, Float64}[]
            Puots = Dict{Float64, Float64}[]

            tseries = [Xts, sgts, z_avts, ug_avts, uo_avts, cgD_Xts, Pzts, Pugts]

            # stst params
            stst_w = 300
            stst_th = 0.01
            nbuffs = length((:X, :z, :sg, :ug, :uo))
            stst = STST(nbuffs, stst_w)

            # Info
            info_frec = 20
            _time = time()
            eltime = 0.0

            # ---------------------------------------------------------------
            # dobreak
            isbreak = false
            function dobreak()

                # niter
                isniter = (S.it >= S.niters)

                # stst
                push!(stst, S.X, S.z_av, S.sg, S.ug_av, S.uo_av)
                dostst1 = isniter
                dostst2 = (S.it >= 2 * stst_w) && iszero(rem(S.it, stst_w))
                dostst = dostst1 || dostst2
                isstst = dostst && is_steady(stst, stst_th)

                # z_tail
                iszero(rem(S.it, z_tail_len)) &&
                    (z_tail_count = 0, z_tail_acc = 0.0)
                z_tail_acc += S.z_av
                z_tail_count += 1
                z_tail = z_tail_acc / z_tail_count

                # dead
                isdead1 = (S.X < dead_Xth)
                isdead2 = (isstst || isniter) && ((S.D - z_tail) / S.D) > z_D_diff_th
                isdead = isdead1 || isdead2

                # explosion
                isexplosion1 = (S.X > exploded_Xth)
                isexplosion2 = (isstst || isniter) && ((z_tail - S.D) / S.D) > z_D_diff_th
                isexplosion = isexplosion1 || isexplosion2

                # ended
                if isdead
                    status = set_status(DEAD_SIM_STATUS, simid, simparams)
                elseif isexplosion
                    status = set_status(EXPLODED_SIM_STATUS, simid, simparams)
                elseif isstst
                    status = set_status(STST_SIM_STATUS, simid, simparams)
                elseif isniter
                    status = set_status(NITERS_SIM_STATUS, simid, simparams)
                end

                isbreak = isdead || isexplosion || isstst || isniter

                isbreak && @info("At isbreak", thid, status, isdead, isexplosion, isstst, isniter)

                return isbreak

            end

            # ---------------------------------------------------------------
            # print_info
            _round(v) = round(v; sigdigits = 4)
            function print_info(msg::String)
                @extract S: it cg D ϵ z_av X dXdt sg ug_av uo_av cgD_X

                # Info
                zdiff = abs(z_av - D)
                ugdiff = abs(ug_av - cgD_X)
                tug_relerr = ugdiff / cgD_X
                D, z_av, zdiff = _round.([D, z_av, zdiff])
                ug_av, uo_av, cgD_X, ugdiff = _round.([ug_av, uo_av, cgD_X, ugdiff])
                @info(msg,
                    progress = string(simi, " of ", tot_sims, " (", div(simi * 100, tot_sims), "%)"),
                    it, thid,
                    status,
                    cg, ϵ,
                    (D, z_av, zdiff),
                    (ug_av, cgD_X, ugdiff),
                    uo_av, X, dXdt,
                    sg, isbreak,
                    tfactor, tug_relerr,
                    eltime
                )
            end

            # ---------------------------------------------------------------
            # feedback
            function feedback()
                @extract S: it cg D ϵ z_av ug_av uo_av X dXdt sg cgD_X

                # info
                eltime = (time() - _time)
                _time = time()
                showinfo = (it == 1)
                showinfo |= iszero(rem(it, info_frec))
                showinfo && print_info("At iter")

                # collecting
                docollect = isbreak
                docollect |= (it == 1)
                docollect |= iszero(rem(it, push_frec))
                docollect &= dosave_batches
                if docollect

                    push!(Xts, X)
                    push!(z_avts, z_av)
                    push!(ug_avts, ug_av)
                    push!(uo_avts, uo_av)
                    push!(cgD_Xts, cgD_X)
                    push!(sgts, sg)

                    if docollect_Ps
                        push!(Pzts, hist(Vz, P))
                        push!(Pugts, hist(Vug, P))
                        push!(Puots, hist(Vuo, P))
                    end

                    # save batch
                    dosave = isbreak
                    dosave |= iszero(rem(it, save_frec))
                    if dosave
                        batchdat = (;
                            push_frec,
                            Xts, sgts, z_avts, ug_avts, uo_avts,
                            cgD_Xts, Pzts, Pugts, Puots
                        )

                        bfile_params = (simid, simparams, S.it)
                        bfile = save_batch(batchdat, bfile_params...)
                        push!(batch_files, bfile_params)

                        empty!.(tseries)
                    end
                end
            end

            # ---------------------------------------------------------------
            # P tranformation
            # globals
            Uug = maximum(Vug) * 1.1
            tkernel = (Vug ./ Uug)
            tauxP = similar(P)
            tfactor = 1.0
            min_tfactor = 1.001
            gdth = 1e-3
            gd_err = 0.0

            # gd fun
            function gd_up_fun(gdmodel)
                @extract S: ug_av cgD_X

                tfactor = gd_value(gdmodel)
                tfactor = max(min_tfactor, tfactor)
                tauxP .= S.P .* (tfactor .- tkernel)
                normalizeP!(tauxP)
                ug_av_ = sum(Vug .* tauxP)

                gd_err = ((cgD_X - ug_av_) / cgD_X)

                return ug_av_
            end

            gd_break_cond(gdmodel) = (gd_err > 0.0 && gd_err < gdth)

            function tranformation()
                @extract S: simit=it cg ug_av X cgD_X
                isinf(cg) && return

                # set up gradient descend
                tauxP .= P
                maxΔx = 5.0
                x0 = tfactor * 0.9
                x1 = tfactor
                target = cgD_X
                maxiter = 100
                verbose = false

                gdmodel_ = grad_desc(gd_up_fun;
                    break_cond = gd_break_cond,
                    target, x0, x1, maxΔx, gdth, maxiter, verbose
                )
                gdit = gdmodel_.iter

                # update P
                S.P .= tauxP

                # info
                showinfo = iszero(rem(simit, info_frec))
                showinfo && @info("At transformation",
                    thid, simit, gdit, tfactor, gd_err
                )

            end

            # ------------------------------------------------------
            # run sim
            run_simD3!(S;
                dobreak,
                tranformation,
                feedback
            )

            # ------------------------------------------------------
            # Save last batch
            feedback()

            # ------------------------------------------------------
            # save sim
            dosave_simdat && save_simdat(S, simid, simparams)
            dosave_batches && save_bfiles(batch_files, simid, simparams)

            # ------------------------------------------------------
            # info
            print_info("At end")

        end # for (D, ϵ, sg)
    end # for thid
end

## ------------------------------------------------------
# finite cg
let
    Ds = range(0.1, 0.5; length = 12)
    ϵs = [
        10.0.^range(log10(1e-3), log10(0.1); length=6);
        range(0.2, 1.0; length=6);
    ]
    cg = 15.0
    simid = "SimD3"

    params = (;Ds, ϵs, cg)
    sglob(Dyn, params, :SimD3, :params, :finite_cg)
    run_sim3D(simid, Ds, ϵs, cg;
        docollect_Ps = false
    )
end

## ------------------------------------------------------
# infinite cg
let
    Ds = range(0.1, 0.5; length = 50)
    ϵs = range(0.01, 1.0; length = 50)
    cg = Inf
    simid = "SimD3"
    params = (;Ds, ϵs, cg)
    sglob(Dyn, params, :SimD3, :params, :infinite_cg)
    run_sim3D(simid, Ds, ϵs, cg;
        dosave_batches = false,
        dosave_simdat = false,
        docollect_Ps = false
    )
end


text:
@with_kw struct GasPlasma{T}
    foo::Int = 6
    ns::Float64
    wp::Float64
    cn2::Float64
    npts::Int64
    Pe::Array{T,1} = zeros(Float64, npts)
    Pe0::Array{T,1} = zeros(Float64, npts)
    Na::Array{T,1} = zeros(Float64, npts)
    Ne::Array{T,1} = zeros(Float64, npts)
    WE::Array{T,1} = zeros(Float64, npts)
    wt::Array{T,1} = zeros(Float64, npts)
    Pe_l::Array{T,1} = zeros(Float64, npts)
    Pe_f::Array{T,1} = zeros(Float64, npts)
    temp1::Array{T,1} = zeros(Float64, npts)
    temp2::Array{T,1} = zeros(Float64, npts)
    temp3::Array{T,1} = zeros(Float64, npts)
    temp4::Array{T,1} = zeros(Float64, npts)
end

struct Gas{T, S}
    type::S #gas type
    p::Float64 #total pressure
    Tk::Float64
    pp::Float64 #partial pressure
    Ip::Float64 #ionization energy
    Reχ::Array{T,1} # refractive index^2 -1
    χ3::Float64
    n2::Float64
    gp::GasPlasma
    #raman to be added
end

function gas(type::Symbol, pp::Float64, Tk::Float64, p::Float64, grid::Grid)
    if type == :Ar
        Ip = 15.75962
    elseif type == :Kr
        Ip = 13.99961
    elseif type == :Ne
        Ip = 21.5646
    elseif type == :He
        Ip = 24.58741
    elseif type == :N2
        Ip = 15.581
    elseif type == :Xe
        Ip = 12.1298
    elseif type == :O2
        Ip = 12.0697
    elseif type == :O3
        Ip = 12.43
    end
    Ip = Ip*e

    Reχ = Realχ(type, p, pp, Tk) # not for O2
    χ3 = χ30(type, pp, Tk)*p
    n2 = n_2(χ3)
    gasplasma = gasPlasma(Ip, grid.npts)
    return Gas{Float64, Symbol}(type, p, Tk, pp, Ip, Reχ.(grid.W), χ3, n2, gasplasma)
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
module Plots

export Plot, Histogram, Histogram2, Linear, Linear3, ErrorBars, Image, Contour, Scatter, Quiver, Node, Circle, Ellipse, Command

using ..ColorMaps
using Compat
using Discretizers

typealias RealRange @compat Tuple{Real,Real}

include("ndgrid.jl")

abstract Plot

type Linear <: Plot
    data::AbstractArray{Real,2}
    mark
    markSize
    style
    legendentry
    onlyMarks
    Linear{T<:Real}(data::AbstractArray{T,2}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) = new(data, mark, markSize, style, legendentry, onlyMarks)
end

type Linear3 <: Plot
    data::AbstractArray{Real,2}
    mark
    markSize
    style
    legendentry
    onlyMarks
    Linear3{T<:Real}(data::AbstractArray{T,2}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) = new(data, mark, markSize, style, legendentry, onlyMarks)
end

const THRESHOLD_NSAMPLES_DISC_OURSELVES = 1000 # if we have more samples than this we discretize ourselves
function _construct_histogram_linear_data{Q<:Real,R<:Real}(
    data::Vector{Q},
    binedges::Vector{R},
    density::Bool, # If true, the bar height will be based on the probability density - otherwise directly on counts
    cumulative::Bool, # A cumulative histogram uses the sum of all previous bins and the current one as final value.
    )

    n = length(binedges)
    disc = LinearDiscretizer(binedges)
    counts = get_discretization_counts(disc, data)
    if cumulative
        cumsum!(counts, counts)
    end

    arr_x = convert(Vector{Float64}, binedges)
    arr_y = convert(Vector{Float64}, counts)
    if density
        arr_y ./= sum(counts)
        arr_y ./= binwidths(disc)
    end
    push!(arr_y, arr_y[end])

    Linear(hcat(arr_x, arr_y)', style="ybar interval,fill=blue!10, draw=blue", mark="none")
end
type Histogram <: Plot
    data::AbstractArray{Real,1}
    bins::Integer
    density::Bool
    cumulative::Bool
    style::AbstractString
    discretization::Symbol
    Histogram(data; bins=10, discretization=:default, density=false, cumulative=false, style="fill=blue!10") = new(data,bins,density,cumulative,style,discretization)
end

type Contour <: Plot
    data::AbstractMatrix{Real}
    xbins
    ybins
    style
    number
    levels
    labels
    Contour(data, xbins, ybins; style=nothing, number=nothing, levels=nothing,labels=nothing) = new(data, xbins, ybins, style, number, levels, labels)
    function Contour(f::Function, xrange::RealRange, yrange::RealRange; xbins=40, ybins=40, style=nothing, number=nothing, levels=nothing, labels=nothing)
        x = linspace(xrange[1], xrange[2], xbins)
        y = linspace(yrange[1], yrange[2], ybins)
        A = zeros(xbins, ybins)
        try
            A = Float64[f(xi, yi) for xi in x, yi in y]
        catch
            A = Float64[f([xi,yi]) for xi in x, yi in y]
        end
        new(A, x, y, style, number, levels, labels)
    end
end

function Contour(z::AbstractMatrix, x::Range, y::Range; style=nothing, number=nothing, levels=nothing, labels=nothing)
    (X, Y) = meshgrid(x, y)
    Contour([X[:]'; Y[:]'; z[:]'], length(x), length(y); style = style, number = number, levels = levels, labels=labels)
end

type Scatter <: Plot
    data::AbstractArray{Any,2}
    mark
    markSize
    style
    legendentry
    onlyMarks
    scatterClasses

    function Scatter{T<:Any}(data::AbstractArray{T,2}; mark=nothing, markSize=nothing, style=nothing, onlyMarks=true, legendentry=nothing, scatterClasses=nothing)
        if size(data,1) == 2
            return Linear(data, mark=mark, markSize=markSize, style=style, onlyMarks=onlyMarks, legendentry=legendentry)
        else
            return new(data, mark, markSize, style, legendentry, onlyMarks, scatterClasses)
        end
    end
end

type ErrorBars <: Plot
    data::AbstractArray{Real,2}
    mark
    style
    legendentry
    ErrorBars{T<:Real}(data::AbstractArray{T,2}; mark=nothing, style=nothing, legendentry=nothing) = new(data, mark, style, legendentry)
end

type Quiver <: Plot
    data::Matrix{Real}
    style
    legendentry
    Quiver{T<:Real}(data::Matrix{T}; style=nothing, legendentry=nothing) = new(data, style, legendentry)
end

type Node <: Plot
    data
    style
    x
    y
    axis # `nothing` will default to "axis cs", other options include "axis description cs", "xticklabel cs", etc.
    Node(data, x, y; style=nothing, axis=nothing) = new(data, style, x, y, axis)
end

type Circle <: Plot
    xc
    yc
    radius
    style
    Circle(xc=0,yc=0,radius=1;style=nothing) = new(xc,yc,radius,style)
end

type Ellipse <: Plot
    xc
    yc
    xradius
    yradius
    style
    Ellipse(xc=0,yc=0,xradius=1,yradius=1;style=nothing) = new(xc,yc,xradius,yradius,style)
end

type Command <: Plot
    cmd::AbstractString
    Command(cmd::AbstractString) = new(cmd)
end

function Quiver(f::Function, xrange::RealRange, yrange::RealRange; style=nothing, legendentry=nothing, samples=15, normalize=true)
    x = linspace(xrange[1], xrange[2], samples)
    y = linspace(yrange[1], yrange[2], samples)
    (X, Y) = meshgrid(x, y)
    n = length(X)
    U = zeros(n)
    V = zeros(n)
    for i = 1:n
        (U[i], V[i]) = f(X[i], Y[i])
    end
    if normalize
        r = max(maximum(U),maximum(V))
        r /= min(minimum(diff(x)),minimum(diff(y)))
        U /= r
        V /= r
    end
    Quiver(X[:], Y[:], U, V, style=style, legendentry=legendentry)
end

Quiver{A<:Real,B<:Real,C<:Real,D<:Real}(x::Vector{A}, y::Vector{B}, u::Vector{C}, v::Vector{D}; style=nothing, legendentry=nothing) = Quiver([x y u v]', style=style, legendentry=legendentry)

Linear{A<:Real, B<:Real}(x::AbstractArray{A,1}, y::AbstractArray{B,1}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) = Linear([x y]', mark=mark, markSize=markSize, style=style, legendentry=legendentry, onlyMarks=onlyMarks)
Linear{A<:Real}(data::AbstractArray{A,1}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) = Linear(collect(1:length(data)), data, mark=mark, markSize=markSize, style=style, legendentry=legendentry, onlyMarks=onlyMarks)

Linear3{A<:Real, B<:Real, C<:Real}(x::AbstractVector{A}, y::AbstractVector{B}, z::AbstractVector{C}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) = Linear3([x y z]', mark=mark, markSize=markSize, style=style, legendentry=legendentry, onlyMarks=onlyMarks)

ErrorBars{A<:Real, B<:Real, C<:Real, D<:Real, E<:Real, F<:Real}(x::AbstractArray{A,1}, y::AbstractArray{B,1}, xplus::AbstractArray{C,1}, yplus::AbstractArray{D,1},
                                                                xminus::AbstractArray{E,1}, yminus::AbstractArray{F,1}; mark=nothing, style=nothing, legendentry=nothing) = ErrorBars([x y xplus yplus xminus yminus]',mark=mark, style=style, legendentry=legendentry)
ErrorBars{A<:Real, B<:Real, C<:Real, D<:Real}(x::AbstractArray{A,1}, y::AbstractArray{B,1}, yplus::AbstractArray{C,1},yminus::AbstractArray{D,1}; mark=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) = ErrorBars([x y zeros(length(x)) yplus zeros(length(x)) yminus]', mark=mark, style=style, legendentry=legendentry)
ErrorBars{A<:Real, B<:Real, C<:Real}(x::AbstractArray{A,1}, y::AbstractArray{B,1}, yplusminus::AbstractArray{C,1}; mark=nothing, style=nothing, legendentry=nothing) = ErrorBars([x y zeros(length(x)) yplusminus zeros(length(x)) yplusminus]', mark=mark, style=style, legendentry=legendentry)

Scatter{A<:Real, B<:Real}(x::AbstractArray{A,1}, y::AbstractArray{B,1}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, scatterClasses=nothing) = Scatter([x y]', mark=mark, markSize=markSize, style=style, legendentry=legendentry, scatterClasses=scatterClasses)
Scatter{A<:Real, B<:Real, C<:Any}(x::AbstractArray{A,1}, y::AbstractArray{B,1}, f::AbstractArray{C,1}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, scatterClasses=nothing) = Scatter([x y f]', mark=mark, markSize=markSize, style=style, legendentry=legendentry, scatterClasses=scatterClasses)
Scatter{A<:Real, B<:Real}(x::A, y::B; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing) = Scatter([x y]', mark=mark, markSize=markSize, style=style, legendentry=legendentry)
Scatter{A<:Real, B<:Real}(x::A, y::B, f; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, scatterClasses=nothing) = Scatter([x y f]', mark=mark, markSize=markSize, style=style, legendentry=legendentry, scatterClasses=scatterClasses)

global _imgid = 1

type Image <: Plot
    filename::AbstractString
    xmin::Real
    xmax::Real
    ymin::Real
    ymax::Real
    zmin::Real
    zmax::Real
    colorbar::Bool
    colormap::ColorMaps.ColorMap
    style
    function Image{T <: Real}(A::Matrix{T}, xrange::RealRange, yrange::RealRange; filename=nothing, colorbar=true, colormap=ColorMaps.Gray(), zmin=nothing, zmax=nothing, style=nothing)
        global _imgid
        if filename == nothing
            id=myid()*10000000000000+_imgid
            filename = "tmp_$(id).png"
            _imgid += 1
        end
        if zmin == nothing
            zmin = minimum(A)
        end
        if zmax == nothing
            zmax = maximum(A)
        end
        if zmin == zmax
            zmin -= 1.
            zmax += 1.
        end
        A = clamp(A, zmin, zmax)
        A = A .- zmin
        A = A ./ (zmax - zmin)
        if isa(colormap, ColorMaps.ColorMap)
            write(colormap, A, filename)
        else
            write(ColorMaps.RGBArray(colormap), A, filename)
        end
        new(filename, xrange[1], xrange[2], yrange[1], yrange[2], zmin, zmax, colorbar, colormap, style)
    end
    function Image(f::Function, xrange::RealRange, yrange::RealRange; filename=nothing, colorbar=true, colormap=ColorMaps.Gray(), zmin=nothing, zmax=nothing, xbins=100, ybins=100, style=nothing)
        x = linspace(xrange[1], xrange[2], xbins)
        y = linspace(yrange[1], yrange[2], ybins)
        (X, Y) = meshgrid(x, y)
        A = map(f, X, Y)
        A = flipdim(A, 1)
        Image(A, xrange, yrange, filename=filename, colorbar=colorbar, colormap=colormap, zmin=zmin, zmax=zmax, style=style)
    end
end

function Histogram2{A<:Real, B<:Real}(x::Vector{A}, y::Vector{B}; xmin=minimum(x), xmax=maximum(x), ymin=minimum(y), ymax=maximum(y), xbins=50, ybins=50, density=false, filename=nothing, colorbar=true, colormap=ColorMaps.Gray(), zmin=nothing, zmax=nothing, style=nothing)
    ex = linspace(xmin, xmax, xbins+1)
    ey = linspace(ymin, ymax, ybins+1)
    ex, ey, M = hist2d(hcat(y, x), ey, ex)
    M = flipdim(M, 1)
    if density
        scale =  xbins * ybins / ((xmax-xmin) * (ymax-ymin) * sum(M))
        M = M * scale
    end
    Image(M, (xmin, xmax), (ymin, ymax), filename=filename, colorbar=colorbar, colormap=colormap, zmin=zmin, zmax=zmax, style=style)
end

end # end plot module


text:
using ESDL
using Test

newcubedir = mktempdir()
ESDLdir(newcubedir)
# Download Cube subset
c = S3Cube()
cgermany = c[
  region = "Germany",
  var = ["gross", "net_ecosystem", "air_temperature_2m", "terrestrial_ecosystem", "soil_moisture"],
  time = 2000:2010
]
saveCube(cgermany,"germanycube", chunksize=(20,20,92,1))
ESDL.ESDLDefaults.cubedir[] = joinpath(newcubedir,"germanycube")

include("access.jl")
include("axes.jl")
include("analysis.jl")
#include("artype.jl")
include("transform.jl")
include("remap.jl")
include("table.jl")
include("tabletocube.jl")


text:
using Distributions
using FileIO
using GeometryBasics
using ImageMagick
using Images
using OffsetArrays
using Plots
using Random
using RecipesPipeline
using RDatasets
using SparseArrays
using StaticArrays
using Statistics
using StatsPlots
using Test
using TestImages

# makie_test_dir = joinpath(@__DIR__, "test_makie")
# mkpath(makie_test_dir)

@testset "RecipesPipeline.jl" begin
    # @testset "Makie integration" begin cd(makie_test_dir) do
    #     include("makie.jl")
    # end end
    @testset "Plots tests" begin
    for i in eachindex(Plots._examples)
        if i ∉ Plots._backend_skips[:gr]
            @test Plots.test_examples(:gr, i, disp=false) isa Plots.Plot
        end
    end
    end
end


text:
__precompile__()
module StochasticPrograms

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using JuMP
using Distributions
using MathProgBase
using MacroTools
using MacroTools: @q, postwalk, prewalk
using Reexport
using ProgressMeter

import Base: getindex, length, in, issubset, show

const MPB = MathProgBase

export
    StochasticModel,
    StochasticProgram,
    Probability,
    AbstractScenario,
    Scenario,
    StageParameters,
    parameters,
    AbstractSampler,
    Sampler,
    ConfidenceInterval,
    StochasticSolution,
    AbstractStructuredSolver,
    AbstractStructuredModel,
    AbstractSampledSolver,
    AbstractSampledModel,
    StructuredModel,
    SampledModel,
    optimize_structured!,
    optimize_sampled!,
    fill_solution!,
    stochastic_solution,
    solverstr,
    internal_solver,
    Crash,
    CrashMethod,
    spsolver,
    spsolver_model,
    distributed,
    deferred,
    set_spsolver,
    internal_model,
    add_scenario!,
    add_scenarios!,
    scenarioproblems,
    decision_length,
    recourse_length,
    first_stage_nconstraints,
    first_stage_dims,
    nstages,
    scenario,
    scenarios,
    scenariotype,
    scenariotext,
    stage_parameters,
    stage_probability,
    probability,
    ExpectedScenario,
    expected,
    has_generator,
    generator,
    subproblem,
    subproblems,
    nsubproblems,
    parentmodel,
    masterterms,
    transfer_model!,
    nscenarios,
    AbstractSampler,
    sample,
    sampler,
    sample!,
    generate!,
    stage_one_model,
    stage_two_model,
    outcome_model,
    evaluate_decision,
    lower_bound,
    lower,
    upper,
    confidence,
    confidence_interval,
    lower_bound,
    upper_bound,
    gap,
    instantiate,
    optimize!,
    decision,
    optimal_decision,
    optimal_value,
    calculate_objective_value,
    @scenario,
    @sampler,
    @first_stage,
    @second_stage,
    @stage,
    @stochastic_model,
    @zero,
    @expectation,
    @sample,
    @decision,
    @parameters,
    @uncertain,
    WS,
    WS_decision,
    EWS,
    SAA,
    DEP,
    VRP,
    EVPI,
    EVP,
    EVP_decision,
    EV,
    EEV,
    VSS

macro exportJuMP()
    Expr(:export, names(JuMP)...)
end
@exportJuMP

include("types/types.jl")
include("methods/methods.jl")
include("spinterface.jl")
include("crash.jl")
include("solvers/solvers.jl")

end # module


text:
if basename(pwd()) == "aoc"
    cd("2018/18")
end

@enum Acre Open = UInt8('.') Trees = UInt8('|') Lumberyard = UInt8('#')

function loadlumber(filename::AbstractString)
    cs = [Acre(UInt8(char)) for line in eachline(filename) for char in line]
    permutedims(reshape(cs, isqrt(length(cs)), :))
end

function prettyprint(as::Matrix{Acre})
    for row in eachrow(as)
        println((Char(UInt8(c)) for c in row)...)
    end
end

function prettyprint(acres::Matrix{Acre}, steps::Integer)
    println("Initial state:")
    prettyprint(acres)
    if steps <= 0
        return
    end
    current, previous = takestep(acres), copy(acres)
    println("\nAfter 1 minute:")
    prettyprint(current)
    for i in 2:steps
        current, previous = step!(previous, current), current
        println("\nAfter ", i, " minutes:")
        prettyprint(current)
    end
    current
end

takestep(acres::AbstractMatrix{Acre}) = step!(similar(acres), acres)

function neighbours(matrix::AbstractMatrix, index::CartesianIndex)
    (
        neighbour
        for neighbour in (index + i for i in filter(!iszero, CartesianIndices((-1:1, -1:1))))
        if checkbounds(Bool, matrix, neighbour)
    )
end

function step!(dst::AbstractMatrix{Acre}, src::AbstractMatrix{Acre})
    copy!(dst, src)
    for (index, acre) in pairs(IndexCartesian(), src)
        ns = collect(src[n] for n in neighbours(src, index))
        if acre == Open && count(==(Trees), ns) >= 3
            dst[index] = Trees
        elseif acre == Trees && count(==(Lumberyard), ns) >= 3
            dst[index] = Lumberyard
        elseif acre == Lumberyard && !(any(==(Trees), ns) && any(==(Lumberyard), ns))
            dst[index] = Open
        end
    end
    dst
end

function info(acres::AbstractMatrix{Acre})
    woods = count(==(Trees), acres)
    yards = count(==(Lumberyard), acres)
    println("there are ", woods, " wooded acres and ", yards, " lumberyards.")
    woods * yards
end

part1(filename::AbstractString) = info(prettyprint(loadlumber(filename), 10))

@assert part1("example.txt") == 1147
part1("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
function WF_loglikelihood_from_adaptive_filtering(α, data, do_the_pruning::Function; silence = false)

    @assert length(α) == length(data[collect(keys(data))[1]])
    Δts = keys(data) |> collect |> sort |> diff |> unique
    if length(Δts) > 1
        test_equal_spacing_of_observations(data; override = false)
    end
    Δt = mean(Δts)

    sα = sum(α)
    times = keys(data) |> collect |> sort
    data, data_2D_array =  prepare_WF_dat_1D_2D(data)

    smmax = values(data) |> sum |> sum
    log_ν_ar = Array{Float64}(undef, smmax, smmax)
    log_Cmmi_ar = Array{Float64}(undef, smmax, smmax)
    log_binomial_coeff_ar_offset = Array{Float64}(undef,    smmax+1, smmax+1)
    μν_prime_im1 = Array{Float64,1}(undef, length(times))

    Λ_of_t = Dict()
    wms_of_t = Dict()

    μν_prime_im1[1] = logμπh_WF(α, zeros(Int64, length(α)), data[times[1]])

    filtered_Λ, filtered_wms = update_WF_params([1.], α, [repeat([0], inner = length(α))], data_2D_array[times[1]])
    Λ_pruned, wms_pruned = do_the_pruning(filtered_Λ, filtered_wms)

    Λ_of_t[times[1]] = filtered_Λ
    wms_of_t[times[1]] = filtered_wms
    new_sm_max = maximum(sum.(Λ_pruned))
    precompute_next_terms_ar!(0, new_sm_max, log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset, sα, Δt)
    sm_max_so_far = new_sm_max

    for k in 1:(length(times)-1)
        if (!silence)
            println("Step index: $k")
            println("Number of components: $(length(filtered_Λ))")
        end
        last_sm_max = maximum(sum.(Λ_pruned))
        new_sm_max = last_sm_max + sum(data[times[k+1]])

        if sm_max_so_far < new_sm_max
            precompute_next_terms_ar!(sm_max_so_far, new_sm_max, log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset, sα, Δt)
            sm_max_so_far = max(sm_max_so_far,new_sm_max)
        end

        filtered_Λ, filtered_wms, μν_prime_im1[k+1] = get_next_WF_filtering_distribution_and_loglikelihood_precomputed(Λ_pruned, wms_pruned, times[k], times[k+1], α, sα, data_2D_array[times[k+1]], log_ν_ar, log_Cmmi_ar, log_binomial_coeff_ar_offset)
        Λ_pruned, wms_pruned = do_the_pruning(filtered_Λ, filtered_wms)
        Λ_of_t[times[k+1]] = filtered_Λ
        wms_of_t[times[k+1]] = filtered_wms
    end

    # return Λ_of_t, wms_of_t
    log_lik_terms = Dict(zip(times, cumsum(μν_prime_im1)))
    return log_lik_terms

end

function logμπh_WF(α::AbstractArray{T, 1}, m::Union{AbstractArray{U, 1}, Tuple}, y::AbstractArray{U, 1}) where {T <: Real, U <: Integer}
    ## Needs to be written for multiple observations too
    sy = sum(y)
    sα = sum(α)
    return SpecialFunctions.logfactorial(sy) - sum(SpecialFunctions.logfactorial.(y)) + sum(log_pochammer.(α .+ m, y)) - log_pochammer(sα + sum(m), sy)
end

function WF_loglikelihood_from_adaptive_filtering_keep_fixed_number(α, data, fixed_number::Int64; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    WF_loglikelihood_from_adaptive_filtering(α, data, prune_keeping_fixed_number; silence = silence)

end

function WF_loglikelihood_from_adaptive_filtering_keep_fixed_fraction(α, data, fraction::Float64; silence = false)
    # println("filter_WF_mem2")
    function prune_keeping_fixed_fraction(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_fraction(Λ_of_t, wms_of_t, fraction)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    WF_loglikelihood_from_adaptive_filtering(α, data, prune_keeping_fixed_fraction; silence = silence)

end

function WF_loglikelihood_from_adaptive_filtering_keep_above_threshold(α, data, ε::Float64; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_above_threshold(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_above_threshold(Λ_of_t, wms_of_t, ε)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    WF_loglikelihood_from_adaptive_filtering(α, data, prune_keeping_above_threshold; silence = silence)

end

log_likelihood_WF_fixed_number = WF_loglikelihood_from_adaptive_filtering_keep_fixed_number

log_likelihood_WF_fixed_fraction = WF_loglikelihood_from_adaptive_filtering_keep_fixed_fraction

log_likelihood_WF_keep_above_threshold = WF_loglikelihood_from_adaptive_filtering_keep_above_threshold

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
VERSION >= v"0.4" && __precompile__()

module ProcessConvolutions

using Base
using Distributions
using PDMats

import PDMats: dim

export
    # Convolution Kernels
    AbstractConvolutionKernel,
    GaussianKernel,
    SquaredExponentialKernel,
    # Basic type
    ProcessConvolution,
    GaussianProcessSample,
    # Methods
    dim,		    # get dimensionality of process or kernel
    conv_wt,	    # Get convolution weights for a new location
    knot_wt,        # Return matrix of conv weights for new locations
    nknot,	        # Return number of knots
    predict        # Give value of GP at new locations

# include convolution kernels
include("ConvolutionKernels.jl")

#----------------------------------------------------------------------------
# Basic type
abstract type AbstractProcess <: Any end

immutable ProcessConvolution{F <: AbstractFloat} <: AbstractProcess
    knot_locs::Array{F, 2}
    knot_values::Vector{F}
    dim::Int
    nknot::Int

    function ProcessConvolution(knot_locs, knot_values)
      new(knot_locs,
          knot_values,
          size(knot_locs, 2),
          size(knot_locs, 1))
    end
end

function ProcessConvolution{F <: AbstractFloat}(knot_locs::Array{F, 2},
                                                dist::UnivariateDistribution)
     ProcessConvolution(knot_locs,
                        rand(dist, size(knot_locs, 1)))
end

ProcessConvolution{F <: AbstractFloat}(knot_locs::AbstractArray) =
    ProcessConvolution(knot_locs, Normal(0, 1))

knot_locs(pc::ProcessConvolution) = pc.knot_locs
knot_values(pc::ProcessConvolution) = pc.knot_values
nknot(pc::ProcessConvolution) = pc.nknot
dim(pc::ProcessConvolution) = pc.dim

immutable ZeroProcess <: AbstractProcess end

immutable ConstantProcess{F <: AbstractFloat} <: AbstractProcess
    v::F
end

#-----------------------------------------------------------------------------
# abstract PredictiveProcessConvolution
#
# type ContinuousPredictivePC{F <: AbstractFloat} <: PredictiveProcessConvolution
#     ProcConv::AbstractProcess
#     ConvKern::AbstractConvolutionKernel
#     PredLocs::Array{F}
#     KnotWt::Array{F}
#     Transform::Function
# end
#
# type DiscretePredictivePC{T <: Any} <: PredictiveProcessConvolution
#     KnotLocs::Array{F}
#     PredLocs::Array{Float}
#     Process::Array{T}
#     KnotValue::Dict{T, Array{Float}}
#     ConvKernel::Dict{T, AbstractConvolutionKernel}
#     KnotWt::Dict{AbstractConvolutionKernel, Array{Float}}
# end
#
# # Outer constructor for simplest case: no knot values specified, one
# # shared kernel for all processes.
# function DiscretePredictivePC(knotlocs::Array{Float},
#                               predlocs::Array{Float},
#                               processlist::Vector{T},
#                               kernel::AbstractConvolutionKernel)
#     nkts = size(knotlocs, 1)
#     npred = size(predlocs, 1)
#     nproc = length(processlist)
#
#     kern = Dict{T, AbstractConvolutionKernel}()
#     for proc in processlist
#         kern[proc] = kernel
#     end
#
#     kv = Dict{T, Vector{Float}}()
#     for proc in processlist
#         kv[proc] = randn(nkts)
#     end
#
#     kw = Dict{AbstractConvolutionKernel,
#               Array{Float}}(kernel => knot_wt(knotlocs, kernel, predlocs))
#
#     DiscretePredictivePC(knotlocs,
#                          predlocs,
#                          processlist,
#                          kv,
#                          kw)
# end
#
# # Outer constructor for different kernel case: no knot values specified, one
# # shared kernel for all processes.
# function DiscretePredictivePC(knotlocs::Array{Float},
#                               predlocs::Array{Float},
#                               kernel::Dict{T, AbstractConvolutionKernel})
#
#     processlist = collect(keys(kernel))
#     nkts = size(knotlocs, 1)
#     npred = size(predlocs, 1)
#     nproc = length(processlist)
#
#     kv = Dict{T, Vector{Float}}()
#     for proc in processlist
#         kv[proc] = randn(nkts)
#     end
#
#     kw = Dict{AbstractConvolutionKernel,
#               Array{Float64, 2}}()
#     for kern in unique(values(kernel))
#         kw[kern] = knot_wt(knotlocs, kern, predlocs)
#     end
#
#     DiscretePredictivePC(knotlocs,
#                          predlocs,
#                          processlist,
#                          kv,
#                          kw)
# end
#
#------------------------------------------------------------------------------
# Putting them together
function predict(pc::ProcessConvolution,
                 kern::AbstractConvolutionKernel,
                 new_loc::AbstractArray)

    nnew = size(new_loc, 1)
    new_val = zeros(nnew)

    for l in 1:nnew
      d = pc.knot_locs' .- new_loc'[:, l]
      new_val[l] = dot(conv_wt(kern, d), pc.knot_values)
    end
    new_val
end

function predict(pc::ProcessConvolution,
                 knot_wt::AbstractArray)
    knot_wt * pc.knot_values
end

## Calculate knot weight matrix
function knot_wt(pc::ProcessConvolution,
                 kern::AbstractConvolutionKernel,
                 new_locs::AbstractArray)
    nloc = size(new_locs, 1)
    nk = nknot(pc)

    k_wt = Array{Float64, 2}(nloc, nk)
    for l in 1:nloc
        k_wt[l, :] = conv_wt(kern, knot_locs(pc)' .- new_locs[l, :]')'
    end
    k_wt
end

function knot_wt(knot_locs::Array,
                 kern::AbstractConvolutionKernel,
                 new_locs::Array)
    nloc = size(new_locs, 1)
    nk = size(knot_locs, 1)

    k_wt = Array{Float64, 2}(nloc, nk)
    for l in 1:nloc
        k_wt[l, :] = conv_wt(kern, knot_locs' .- new_locs[l, :]')'
    end
    k_wt
end

# Include efficient sample storage type
include("sampling.jl")

end  # module


text:
include(joinpath("Duration", "decisionTree_test.jl"))
include(joinpath("Duration", "randomForest_test.jl"))
include(joinpath("Duration", "predByMean_test.jl"))
include(joinpath("Duration", "linReg_test.jl"))
include(joinpath("Duration", "gammaReg_test.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
module IdealPoints
    import StatsFuns: logistic
    import Distributions: logpdf, Normal
    import Optim
    import RollCallDataIO
    import RollCallDataIO: ORDFile, SparseRollCall

    include("nd_model.jl")
end


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

mutable struct ManagementPolicyBaseBlob <: SwaggerModel
    tierToCool::Any # spec type: Union{ Nothing, DateAfterModification } # spec name: tierToCool
    tierToArchive::Any # spec type: Union{ Nothing, DateAfterModification } # spec name: tierToArchive
    delete::Any # spec type: Union{ Nothing, DateAfterModification } # spec name: delete

    function ManagementPolicyBaseBlob(;tierToCool=nothing, tierToArchive=nothing, delete=nothing)
        o = new()
        validate_property(ManagementPolicyBaseBlob, Symbol("tierToCool"), tierToCool)
        setfield!(o, Symbol("tierToCool"), tierToCool)
        validate_property(ManagementPolicyBaseBlob, Symbol("tierToArchive"), tierToArchive)
        setfield!(o, Symbol("tierToArchive"), tierToArchive)
        validate_property(ManagementPolicyBaseBlob, Symbol("delete"), delete)
        setfield!(o, Symbol("delete"), delete)
        o
    end
end # type ManagementPolicyBaseBlob

const _property_map_ManagementPolicyBaseBlob = Dict{Symbol,Symbol}(Symbol("tierToCool")=>Symbol("tierToCool"), Symbol("tierToArchive")=>Symbol("tierToArchive"), Symbol("delete")=>Symbol("delete"))
const _property_types_ManagementPolicyBaseBlob = Dict{Symbol,String}(Symbol("tierToCool")=>"DateAfterModification", Symbol("tierToArchive")=>"DateAfterModification", Symbol("delete")=>"DateAfterModification")
Base.propertynames(::Type{ ManagementPolicyBaseBlob }) = collect(keys(_property_map_ManagementPolicyBaseBlob))
Swagger.property_type(::Type{ ManagementPolicyBaseBlob }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagementPolicyBaseBlob[name]))}
Swagger.field_name(::Type{ ManagementPolicyBaseBlob }, property_name::Symbol) =  _property_map_ManagementPolicyBaseBlob[property_name]

function check_required(o::ManagementPolicyBaseBlob)
    true
end

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


text:
fsum(x::Int,p...)=x+fsum(p[1],p[2:end]...)
fsum(x::Int) = x


text:
module ParsePSSE

function parse_skipsection(lines::Array{String}, pos::Int)
    while pos <= length(lines)
        line = strip(lines[pos])

        if line[1] in [ '0', 'Q' ]
            break
        end

        pos += 1
    end

    return pos
end

function parse_data(lines::Array{String}, pos::Int, info::Array;
                  delim=",")
    start = last = pos
    last = parse_skipsection(lines, start)
    ndata = last - start
    data = Array{Any}(undef, ndata, length(info))

    for i=1:ndata
        # Get rid of the single quotes.
        line = replace(lines[start+i-1], "'" => "")
        fields = split(line, delim)

        for j=1:length(info)
            if info[j][3] == String
                data[i,j] = strip(fields[info[j][2]])
            else
                data[i,j] = parse(info[j][3], fields[info[j][2]])
            end
        end
    end

    pos = last + 1
    return pos, data
end

include("parse_raw.jl")
include("parse_rop.jl")
include("parse_inl.jl")
include("parse_con.jl")

end


text:
using NQCDynamics: RingPolymers

function RingPolymerSimulation{Ehrenfest}(atoms::Atoms{S,T}, model::Model, n_beads::Integer; kwargs...) where {S,T}
    RingPolymerSimulation(atoms, model, Ehrenfest{T}(NQCModels.nstates(model)), n_beads; kwargs...)
end

function DynamicsMethods.motion!(du, u, sim::RingPolymerSimulation{<:Ehrenfest}, t)
    dr = DynamicsUtils.get_positions(du)
    dv = DynamicsUtils.get_velocities(du)
    dσ = DynamicsUtils.get_quantum_subsystem(du)

    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)
    σ = DynamicsUtils.get_quantum_subsystem(u)

    DynamicsUtils.velocity!(dr, v, r, sim, t)
    Calculators.update_electronics!(sim.calculator, r)
    acceleration!(dv, v, r, sim, t, σ)
    DynamicsUtils.apply_interbead_coupling!(dv, r, sim)
    DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim)
end

function DynamicsUtils.classical_hamiltonian(sim::RingPolymerSimulation{<:Ehrenfest}, u)
    v = DynamicsUtils.get_velocities(u)
    r = DynamicsUtils.get_positions(u)
    kinetic = DynamicsUtils.classical_kinetic_energy(sim, v)
    spring = RingPolymers.get_spring_energy(sim.beads, sim.atoms.masses, r)

    Calculators.evaluate_potential!(sim.calculator, r)
    Calculators.eigen!(sim.calculator)
    population = Estimators.adiabatic_population(sim, u)
    potential = sum([dot(population, eigs.values) for eigs in sim.calculator.eigen])

    return kinetic + potential + spring
end


text:
# MMMMMMMMMMMMMMMMMMMMMMMMMMMWNNKN0KMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
# MMMMMMMMMMMMMMMMMMMMMMMMMW0OMMMMX0MMMMMMMMMMXXNMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMNWMMMMNMMMMMMMMMMMMMMMMMMMMMMMM
# MMMMMMMMMMMMMMMMMMMMMMMMMKMMMMMOKWMMMMMMMN:lkkdXMMNXWMMMWMNNMMMMWXWMMMNNMMMWMNNMMMlXWM;.xMMW.MMMWXWMM0dWWMMMMMMMMMMMM
# MMMMMMMMMMMMMMMMMMMMMMMMNXMMMMMXXMMMMMMMM:oMMMMMX,xkl:MK.xkd'MN,dkxMx:Ox,NW dkd'N0 xOM;OccWW.MW;dOcoN:;kNMMMMMMMMMMMM
# MMMMMMMMMMMMMMMMMMWXXXXXkMMMMMMM0MMMMMMMM:lMMMMMocMMM.XK.MMM Mk;MMMM.lxxlKM MMM;xM.NMM;OMk,X.MO.xxdoMolMMMMMMMMMMMMMM
# MMMMMMMMMMMMMMXXXXNMMMMMoMMMMMMMOMMMMMMMMW:cdxdWN;dxlcMK.MMM MW;oxdMk;dxdWM dxo,NM;lkM;OMMX'.MW:oxd0MO'xXMMMMMMMMMMMM
# MMMMMMMMMMMMXXWMMMMMMMMM0KMMMMMOOMMMMMMMMMMMWWMMMMMWMMMMMMMMMMMMMWMMMMWWMMM MWWMMMMWWMMMMMMMMMMMMWWMMMMWMMMMMMMMMMMMM
# MMMMMMMMMMW0WMMMMMMMMMMMMNMMMM0:WMMMMMMMMXxxNMM0kMMMMMMMMMMMMMMMMMMMMXl0MMMxMMMMMMMMMMMMModMMMMMMMMMMMMW0WMMMMMMMMdoM
# MMMMMMMMMM0MMMMMMMMMMMMMMMMMMk'0MMMMMMMMM0cxlXMxoMNxKMWd0M0xOxdOkxd0MXcxkdkWM0xkx0MNd0xxMlokdxNMNxkxkNNxcxkMXxxxOMdlk
# MMMMMMMMMMKMMMMMMMMMNK0XNMMMM:'OWMMMMMMMM0cNOlKkoMXcOMWckMkcXMxcNMdlMXcOMXcOXcdOxcKNcxWMMloMWooMN0OkcOMOcNMWldMMMMdlW
# NKOkk0KK0OkMMMMMMWd;:kWMWXX0kdo;'OMMMMMMM0cNMKlloMNlxXOckMkcWMxcMMdlMXcxXOcKNlxXXKWNcOMMMllKKlxMloX0cOM0c0XMooKXKMdoM
# MMMMWX0kdl:0MMMMX;'lWMMMMMMWk;'0XWMMMMMMMN0WMMX0KMMX00XKXMX0WMX0MMKKMW0N00XMMWK00KMW0XMMMKXK0KWMX00NKXMWK00MWK00XMKKM
# MMMMMMMMMMMWOOoc'o0KWMMMMMMMMWK0MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM

# ConceptnetNumberbatch.jl - an interface for ConceptNetNumberbatch written in Julia at 0x0α Research,
#                            by Corneliu Cofaru, 2018
# Paper:
#   Robert Speer, Joshua Chin, and Catherine Havasi (2017).
#       "ConceptNet 5.5: An Open Multilingual Graph of General Knowledge." in proceedings of AAAI 2017.

module ConceptnetNumberbatch

using TranscodingStreams
using CodecZlib
using HDF5
using Languages

import Base: get, getindex, size, length, show, keys, values, in

export CONCEPTNET_MULTI_LINK,
       CONCEPTNET_EN_LINK,
       CONCEPTNET_HDF5_LINK,
       ConceptNet,
       download_embeddings,
       load_embeddings,
       tokenize_for_conceptnet,
       embed_document

include("defaults.jl")
include("interface.jl")
include("files.jl")
include("document_embeddings.jl")

end # module


text:
#=
This file demonstrates the framework of functions and data structures necessary to implement a task
=#

export Task_NoTask

#=
Data structure
should be a subtype of Task abstract type
=#

mutable struct Task_NoTask <: Task
end

#=
Constructors for data type
=#

#=
Initialization function
This will build all of the necessary elements before anything starts running
(initializing external boards, creating GUIs etc)
=#
function init_task(myt::Task_NoTask,rhd::RHD2000,han,fpga)
end

#=
Experimental Control Function

This will implement the control logic of the task
such as updating GUIs, modifying the data structure, talking to external boards
=#
function do_task(myt::Task_NoTask,rhd::RHD2000,myread,han,fpga)
end

#=
Logging Function

This will save the appropriate elements of the data structure, as well as specifying what
analog streams from either the Intan or other external DAQs
=#
function save_task(myt::Task_NoTask,rhd::RHD2000)
end