text:
import Pkg; Pkg.instantiate();

using GeometryBasics, Colors
import Observables
import Observables: on, off, notify!

import GLMakie
import GLFW: SetWindowTitle
using AbstractPlotting
using AbstractPlotting.MakieLayout

include("../src/mcmc.jl");
include("../src/visualization.jl")

### GUI

scene, layout = layoutscene(resolution=(768,1024))

state_plot = layout[1,1:2] = LScene(scene, raw=true, camera = cam2d!, autolimitaspect=1.0)
sub_plot_energy = layout[3,1] = LAxis(scene, title="Energy")
sub_plot_mag = layout[3,2] = LAxis(scene, title="Magnetization")

g_T = labelslider!(scene, "temp", 0f0:0.005f0:5f0)
T = g_T.slider.value
g_h = labelslider!(scene, "field", -5f0:0.1f0:5f0)
h = g_h.slider.value
g_mult = labelslider!(scene, "speed", 0.0:0.05:1.0, sliderkw=Dict(:startvalue=>0.5))
mult = g_mult.slider.value
sub_controls = layout[2,:] = GridLayout()
sub_controls[1,1] = g_T.layout
sub_controls[1,2] = g_h.layout
sub_controls[2,1:2] = g_mult.layout

rowsize!(layout, 3, Relative(1/6))

color1 = colorant"black"
color2 = colorant"cornflowerblue"

### Logic

algorithm = Node(:metropolis) # either :metropolis or :wolff

# Timey-Wimey
frame_node = Node(0)
run_signal = Node(false)
BASE_FPS = Node(60)
BUFFER_LENGTH = 1024 # no. of stored observations

L = Node(128) # linear system size

pos_and_marker = lift(get_primitives, L)
positions = lift(x->reshape(x[1], length(x[1])), pos_and_marker)
square = lift(x->x[2], pos_and_marker)
config0 = lift(frustratedConfiguration, L)
cluster = lift(L->zeros(Bool, L,L), L)

color_node = lift(L->[color1 for j in 1:L^2], L)

size_change = lift(L) do L
    run_signal[] = false

    update_cam!(state_plot.scene, FRect(-1,-1,2,2))
    display(scene);

    nothing
end

color_update = lift(config0) do config0
    # TODO: recalculating all colors is wasteful
    # -> return flipped spins -> adjust color
    color_node[] = reshape(color_gen(config0, color1, color2),L[]^2)
end

# Plot of the lattice
meshscatter!(state_plot, positions; color=color_node,
    markersize=1, camera=cam2d!, marker=square,
    raw=true, shading=false
    )
state_plot.scene[end][:light] = Vec{3,Float32}[[1.0, 1.0, 1.0], [0.1, 0.1, 0.1], [0.9, 0.9, 0.9], [0.0, 0.0, -20.0]]

## Events

toggle_run!() = run_signal[] = !run_signal[]
toggle_algo!() = algorithm[] = algorithm[] == :metropolis ? :wolff : :metropolis

function set_title!()
    title_string = "Ising - $(titlecase(string(algorithm[])))"
    if mult[] == 0.0 || !run_signal[]
        title_string *= " - PAUSED"
    end
    SetWindowTitle(GLMakie.gl_screens[1], title_string)
end

# Keys:
#   - r: reset state and plots
#   - c: reset camera
#   - +/- (QWERTY): adjust simulation speed
#   - s: short sweep: 1 time
#   - l: long sweep 1_000 times
#   - a: switch algorithm Metropolis/Wolff
on(scene.events.keyboardbuttons) do buttons
    if ispressed(scene, Keyboard.p)
        # if the renderloop ended for some reason
        # it is restarted here
        if istaskdone(render_task)
            global render_task = renderloop()
        end
        toggle_run!()
    elseif ispressed(scene, Keyboard.c)
        update_cam!(state_plot.scene, FRect(-1,-1,2,2))
    elseif ispressed(scene, Keyboard.r)
        config0[] = frustratedConfiguration(L[]);
        cluster[] = zeros(Bool, L[],L[]);
        empty!(time_buffer[])
        empty!(E_buffer[])
        empty!(M_buffer[])
        notify!(config0)
    elseif ispressed(scene, Keyboard.s)
        run_signal[] = false
        _sweep!(config0, L[]^2)
        _measure!()
    elseif ispressed(scene, Keyboard.l)
        run_signal[] = false
        _sweep!(config0, 1_000*L[]^2)
    elseif ispressed(scene, Keyboard.equal)
        mult[] = min(mult[]+0.05, 1.0)
    elseif ispressed(scene, Keyboard.minus)
        mult[] = max(mult[]-0.05, 0.0)
    elseif ispressed(scene, Keyboard.a)
        toggle_algo!()
    end
    set_title!()
end
on(scene.events.window_area) do wa
    update_cam!(state_plot.scene, FRect(-1,-1,2,2))
end

## Physical observables
time_buffer = Node([0.0])
E_buffer = Node([Point2(0.0,0.0)]) # pairs of (time, value)
M_buffer = Node([Point2(0.0,0.0)])

E_plot = lines!(sub_plot_energy, E_buffer)
ylims!(sub_plot_energy, (-2.1, 0.1))
M_plot = lines!(sub_plot_mag, M_buffer)
ylims!(sub_plot_mag, (-1.1, 1.1))

# recalculate the plot limits when observables change.
# MakieLayout doesn't  adjust them automatically when
# new points get pushed to the plot.
plot_limits = Observables.on(E_buffer) do a
    trange = (E_buffer[][1][1], E_buffer[][end][1]+mult[])
    xlims!(sub_plot_energy, trange)
    xlims!(sub_plot_mag, trange)
    nothing
end

## Render- and simulation loop

function _sweep!(config0, n)
    if algorithm[] == :metropolis
        metropolis_sweep!(config0[], n, 1/T[], h[])
    elseif algorithm[] == :wolff
        wolff_sweep!(config0[], cluster[], ceil(Int, n/L[]^2), 1/T[], h[])
    end
    notify!(config0)
    nothing
end

function _measure!()
    t = if length(time_buffer[])>0
            time_buffer[][end]
        else
             0.0
        end + mult[]
    push!(time_buffer[], t)
    push!(E_buffer[], Point2(t, H(config0[], 1.0, h[])/L[]^2))
    push!(M_buffer[], Point2(t, m(config0[])))
    if length(E_buffer[]) > BUFFER_LENGTH
        deleteat!(E_buffer[], 1)
        deleteat!(M_buffer[], 1)
        deleteat!(time_buffer[], 1)
    end
    Observables.notify!(E_buffer)
    Observables.notify!(M_buffer)
    nothing
end

function renderloop()
    @async begin
        while scene.events.window_open[]
            t1 = time()
            if run_signal[] && mult[] > 0.0
                _sweep!(config0, round(Int, mult[] * L[]^2))
                # Measure observables
                # and push them to the relevant buffers.
                # Prune buffers if necessary.
                _measure!()
            end
            t2 = time()
            if (t2-t1) < 1/BASE_FPS[]
                sleep(1/BASE_FPS[] - (t2-t1))
            end
        end
        :ended
    end
end

render_task = renderloop()

display(scene)

max_stars_repo_name: skleinbo/Ising

max_issues_repo_name: Yaren-Github/Ising

max_forks_repo_name: Yaren-Github/Ising


text:
@doc """ multi-class logarithmic loss
""" ->
function mlogloss{T<:AbstractFloat}(pd::Array{T,2}, gt::Array{T,2})
    # T can't be Real since raised an integer x to a negative power -n.
    function max_min(x::T)
        max(min(x,1-10.0^(-15)),10.0^(-15))
    end
    pd = map(max_min, pd)
    log_pd = log(pd)
    -sum(sum(gt .* log_pd, 1))
end

ext: jl


text:
module SnoopBench

# Assignment of parcel to modules
struct A end
f3(::A) = 1
f2(a::A) = f3(a)
f1(a::A) = f2(a)

# Like map! except it uses push!
# With a single call site
mappushes!(f, dest, src) = (for item in src push!(dest, f(item)) end; return dest)
mappushes(@nospecialize(f), src) = mappushes!(f, [], src)
function mappushes3!(f, dest, src)
    # A version with multiple call sites
    item1 = src[1]
    push!(dest, item1)
    item2 = src[2]
    push!(dest, item2)
    item3 = src[3]
    push!(dest, item3)
    return dest
end
mappushes3(@nospecialize(f), src) = mappushes3!(f, [], src)

# Useless specialization
function spell_spec(::Type{T}) where T
    name = Base.unwrap_unionall(T).name.name
    str = ""
    for c in str
        str *= c
    end
    return str
end
function spell_unspec(@nospecialize(T))
    name = Base.unwrap_unionall(T).name.name
    str = ""
    for c in str
        str *= c
    end
    return str
end

end


text:
# Test mapping arrays directly to aesthetics without a data fram.

using Gadfly, DataArrays, DataFrames

plot(x=collect(1:100), y=sort(rand(100)),
     Guide.xlabel("index"), Guide.ylabel("position"), Guide.title("Gadfly Rules"))


text:
const METIS_NOPTIONS = 40
## Return codes
const METIS_OK           = @compat Int32(1)  # normal return
const METIS_ERROR_INPUT  = @compat Int32(-2) # erroneous inputs and/or options
const METIS_ERROR_MEMORY = @compat Int32(-3) # insufficient memory
const METIS_ERROR        = @compat Int32(-4) # Other errors
## Operation type codes
const METIS_OP_PMETIS = @compat Int32(0)
const METIS_OP_KMETIS = @compat Int32(1)
const METIS_OP_OMETIS = @compat Int32(2)
## (1-based) positions in options vector
const METIS_OPTION_PTYPE     = 1
const METIS_OPTION_OBJTYPE   = 2
const METIS_OPTION_CTYPE     = 3
const METIS_OPTION_IPTYPE    = 4
const METIS_OPTION_RTYPE     = 5
const METIS_OPTION_DBGLVL    = 6
const METIS_OPTION_NITER     = 7
const METIS_OPTION_NCUTS     = 8
const METIS_OPTION_SEED      = 9
const METIS_OPTION_NO2HOP    = 10
const METIS_OPTION_MINCONN   = 11
const METIS_OPTION_CONTIG    = 12
const METIS_OPTION_COMPRESS  = 13
const METIS_OPTION_CCORDER   = 14
const METIS_OPTION_PFACTOR   = 15
const METIS_OPTION_NSEPS     = 16
const METIS_OPTION_UFACTOR   = 17
const METIS_OPTION_NUMBERING = 18
const METIS_OPTION_HELP      = 19
const METIS_OPTION_TPWGTS    = 20
const METIS_OPTION_NCOMMON   = 21
const METIS_OPTION_NOOUTPUT  = 22
const METIS_OPTION_BALANCE   = 23
const METIS_OPTION_GTYPE     = 24
const METIS_OPTION_UBVEC     = 25
## Partitioning Schemes
const METIS_PTYPE_RB   = @compat Int32(0)
const METIS_PTYPE_KWAY = @compat Int32(1)
## Graph types for meshes
const METIS_GTYPE_DUAL  = @compat Int32(0)
const METIS_GTYPE_NODAL = @compat Int32(1)
## Coarsening Schemes
const METIS_CTYPE_RM   = @compat Int32(0)
const METIS_CTYPE_SHEM = @compat Int32(1)
## Initial partitioning schemes
const METIS_IPTYPE_GROW    = @compat Int32(0)
const METIS_IPTYPE_RANDOM  = @compat Int32(1)
const METIS_IPTYPE_EDGE    = @compat Int32(2)
const METIS_IPTYPE_NODE    = @compat Int32(3)
const METIS_IPTYPE_METISRB = @compat Int32(4)
## Refinement schemes
const METIS_RTYPE_FM        = @compat Int32(0)
const METIS_RTYPE_GREEDY    = @compat Int32(1)
const METIS_RTYPE_SEP2SIDED = @compat Int32(2)
const METIS_RTYPE_SEP1SIDED = @compat Int32(3)
## Debug levels (bit positions)
const METIS_DBG_INFO       = @compat Int32(1) # Shows various diagnostic messages
const METIS_DBG_TIME       = @compat Int32(2) # Perform timing analysis
const METIS_DBG_COARSEN    = @compat Int32(4) # Show the coarsening progress
const METIS_DBG_REFINE     = @compat Int32(8) # Show the refinement progress
const METIS_DBG_IPART      = @compat Int32(16) # Show info on initial partitioning
const METIS_DBG_MOVEINFO   = @compat Int32(32) # Show info on vertex moves during refinement
const METIS_DBG_SEPINFO    = @compat Int32(64) # Show info on vertex moves during sep refinement
const METIS_DBG_CONNINFO   = @compat Int32(128) # Show info on minimization of subdomain connectivity
const METIS_DBG_CONTIGINFO = @compat Int32(256) # Show info on elimination of connected components
const METIS_DBG_MEMORY     = @compat Int32(2048) # Show info related to wspace allocation
## Types of objectives
const METIS_OBJTYPE_CUT  = @compat Int32(0)
const METIS_OBJTYPE_VOL  = @compat Int32(1)
const METIS_OBJTYPE_NODE = @compat Int32(2)


text:
# This file is part of Kpax3. License is MIT.

"""
# User defined settings for a Kpax3 run

## Description

## Fields

"""
struct KSettings
  ifile::String
  ofile::String
  # Common parameters
  protein::Bool
  miss::Vector{UInt8}
  misscsv::Vector{String}
  l::Int
  α::Real
  θ::Real
  γ::Vector{Float64}
  r::Float64
  maxclust::Int
  maxunit::Int
  verbose::Bool
  verbosestep::Int
  # Genetic Algorithm parameters
  popsize::Int
  maxiter::Int
  maxgap::Int
  xrate::Float64
  mrate::Float64
  # MCMC parameters
  T::Int
  burnin::Int
  tstep::Int
  op::StatsBase.ProbabilityWeights
end

function KSettings(ifile::String,
                   ofile::String;
                   protein::Bool=true,
                   miss::Vector{UInt8}=zeros(UInt8, 0),
                   misscsv::Vector{String}=Array{String}(undef, 0),
                   l::Int=100000000,
                   alpha::Real=0.5,
                   theta::Real=-0.25,
                   gamma::Vector{Float64}=[0.6; 0.35; 0.05],
                   r::Float64=log(0.001) / log(0.95),
                   maxclust::Int=500,
                   maxunit::Int=500,
                   verbose::Bool=false,
                   verbosestep::Int=1000,
                   popsize::Int=20,
                   maxiter::Int=20000,
                   maxgap::Int=5000,
                   xrate::Float64=0.9,
                   mrate::Float64=0.005,
                   T::Int=100000,
                   burnin::Int=10000,
                   tstep::Int=1,
                   op::Vector{Float64}=[0.5; 0.0; 0.5])
  # open files and immediately close them. We do this to throw a proper Julia
  # standard exception if something is wrong
  f = open(ifile, "r")
  close(f)

  dirpath = dirname(ofile)
  if !isdir(dirpath)
    mkpath(dirpath)
  end

  #=
  f = open(string(ofile, "_settings.bin"), "a")
  close(f)

  f = open(string(ofile, "_row_partition.bin"), "a")
  close(f)

  f = open(string(ofile, "_col_partition.bin"), "a")
  close(f)
  =#

  if length(miss) == 0
    miss = if protein
             UInt8['?', '*', '#', '-', 'b', 'j', 'x', 'z']
           else
             UInt8['?', '*', '#', '-', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's',
                   'v', 'w', 'x', 'y', 'j', 'z']
           end
  else
    if length(miss) == 1
      if miss[1] != UInt8(0)
        if UInt8(0) < miss[1] < UInt8(128)
          if UInt8(64) < miss[1] < UInt8(91)
            # convert to lowercase
            miss[1] += UInt8(32)
          end
        else
          throw(KDomainError(string("Value 'miss[1]' is not in the range ",
                                    "[1, ..., 127]: ", Int(miss[1]), ".")))
        end
      end
    else
      for i in 1:length(miss)
        if UInt8(0) < miss[i] < UInt8(128)
          if UInt8(64) < miss[i] < UInt8(91)
            # convert to lowercase
            miss[i] += UInt8(32)
          end
        else
          throw(KDomainError(string("Value 'miss[", i, "]' is not in the ",
                                    "range [1, ..., 127]: ",
                                    Int(miss[i]), ".")))
        end
      end
    end
  end

  if length(misscsv) == 0
    misscsv = [""]
  elseif !("" in misscsv)
    insert!(misscsv, 1, "")
  end

  if l < 1
    throw(KDomainError("Argument 'l' is not positive."))
  end

  if length(gamma) != 3
    throw(KInputError("Argument 'gamma' does not have length 3."))
  elseif gamma[1] < 0
    throw(KDomainError("Argument 'gamma[1]' is negative."))
  elseif gamma[2] < 0
    throw(KDomainError("Argument 'gamma[2]' is negative."))
  elseif gamma[3] < 0
    throw(KDomainError("Argument 'gamma[3]' is negative."))
  end

  if r <= 0.0
    throw(KDomainError("Argument 'r' is not positive."))
  end

  if maxclust < 1
    throw(KDomainError("Argument maxclust is lesser than 1."))
  end

  if maxunit < 1
    throw(KDomainError("Argument maxunit is lesser than 1."))
  end

  if verbosestep < 0
    throw(KDomainError("Argument 'verbosestep' is negative."))
  end

  # disable status reports if verbosestep is not positive
  verbose = verbose && (verbosestep > 0)

  if popsize < 4
    throw(KDomainError("Argument 'popsize' is lesser than 4."))
  end

  if maxiter < 1
    throw(KDomainError("Argument 'maxiter' is lesser than 1."))
  end

  if maxgap < 0
    throw(KDomainError("Argument 'maxgap' is lesser than 1."))
  end

  if !(0 <= xrate <= 1)
    throw(KDomainError("Argument 'xrate' is not in the range [0, 1]."))
  end

  if !(0 <= mrate <= 1)
    throw(KDomainError("Argument 'mrate' is not in the range [0, 1]."))
  end

  if T < 1
    throw(KDomainError("Argument 'T' is lesser than 1."))
  end

  if burnin < 0
    throw(KDomainError("Argument 'burnin' is negative."))
  end

  if tstep < 0
    throw(KDomainError("Argument 'tstep' is negative."))
  end

  if length(op) != 3
    throw(KInputError("Argument 'op' does not have length 3."))
  elseif op[1] < 0
    throw(KDomainError("Argument 'op[1]' is negative."))
  elseif op[2] < 0
    throw(KDomainError("Argument 'op[2]' is negative."))
  elseif op[3] < 0
    throw(KDomainError("Argument 'op[3]' is negative."))
  end

  KSettings(ifile, ofile, protein, miss, misscsv, l, alpha, theta, gamma, r,
            maxclust, maxunit, verbose, verbosestep, popsize, maxiter, maxgap,
            xrate, mrate, T, burnin, tstep, StatsBase.ProbabilityWeights(op))
end


text:
# cconvert: Not exported
# size: interfaces.jl
# asize_from: Not exported
# length: interfaces.jl
# elsize: Not exported
# sizeof: TODO
# strides: interfaces.jl
# isassigned: interfaces.jl
# unsafe_copy!: overwriting
# copy!: overwriting
# copy: abstractarray
# reinterpret: TODO
# reshape: abstractarray
# similar: interfaces.jl
# getindex: interfaces.jl
# fill!: overwriting
# fill: TODO
# cell: TODO
# zeros: interfaces.jl
# zero: interfaces.jl
# ones: interfaces.jl
# one: interfaces.jl
# eye
# convert
# promote_rule
# collect
# start: interfaces.jl
# next: interfaces.jl
# done: interfaces.jl
# unsafe_getindex: Not exported
# setindex!: overwriting
# unsafe_setindex!: Not exported
# _growat!: Not exported
# _growat_beg!: Not exported
# _growat_end!: Not exported
# _deleteat!: Not exported
# _deleteat_beg!: Not exported
# _deleteat_end!: Not exported
# push!: overwriting
# append!: overwriting
# prepend!: overwriting
# resize!: overwriting
# sizehint!: overwriting
# pop!: overwriting
# unshift!: overwriting
# shift!: overwriting
# insert!: overwriting
# deleteat!: overwriting
# splice!: overwriting
# empty!: overwriting
# lexcmp: TODO
# reverse: TODO
# reverseind: TODO
# reverse!: overwriting
# vcat: abstractarray
# hcat: abstractarray
# findnext
# findfirst
# findprev
# findlast
# find
# findn
# findnz
# findmax
# findmin
# indmax
# indmin
# indexin
# findin
# indcopy: Not exported
# filter
# filter!: overwriting
# intersect
# union
# setdiff
# symdiff


text:
export PermGroup, perm, parity, elements, cycles, character, setpermstyle,
       order

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent_type(::Type{perm}) = PermGroup

elem_type(::Type{PermGroup}) = perm

###############################################################################
#
#   Basic manipulation
#
###############################################################################

doc"""
    parent(a::perm)
> Return the parent of the given permutation group element.
"""
parent(a::perm) = a.parent

function deepcopy_internal(a::perm, dict::ObjectIdDict)
   G = parent(a)
   return G(deepcopy(a.d))
end

function Base.hash(a::perm, h::UInt)
   b = 0x595dee0e71d271d0%UInt
   for i in 1:length(a.d)
      b = xor(b, xor(hash(a[i], h), h))
      b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
   end
   return b
end

doc"""
    parity(a::perm)
> Return the parity of the given permutation, i.e. the parity of the number of
> transpositions that compose it. The function returns $1$ if the parity is odd
> otherwise it returns $0$.
"""
# TODO: 2x slower than Flint
function parity(a::perm)
   to_visit = trues(a.d)
   parity = length(to_visit)
   k = 1
   while any(to_visit)
      parity -= 1
      k = findnext(to_visit, k)
      to_visit[k] = false
      next = a[k]
      while next != k
         to_visit[next] = false
         next = a[next]
      end
   end
   return parity%2
end

doc"""
    sign(a::perm)
> Returns the sign of the given permutations, i.e. `1` if `a` is even and `-1`
> if `a` is odd.
"""
sign(a::perm) = (-1)^parity(a)

function getindex(a::perm, n::Int)
   return a.d[n]
end

function setindex!(a::perm, d::Int, n::Int)
   a.d[n] = d
end

doc"""
    eye(G::PermGroup)
> Return the identity permutation for the given permutation group.
"""
eye(G::PermGroup) = G()

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, G::PermGroup)
   print(io, "Permutation group over $(G.n) elements")
end

mutable struct PermDisplayStyle
   format::Symbol
end

const perm_display_style = PermDisplayStyle(:array)

doc"""
    setpermstyle(format::Symbol)
> Nemo can display (in REPL or in general as string) permutations by either
> array of integers whose `n`-th position represents value at `n`, or as
> (an array of) disjoint cycles, where the value of permutation at `n` is
> represented as the entry immediately following `n` in the cycle.
> The style can be switched by calling `setpermstyle` with `:array` or
> `:cycles` argument.
"""
function setpermstyle(format::Symbol)
   format in (:array, :cycles) || throw("Permutations can be displayed
   only as :array or :cycles.")
   perm_display_style.format = format
   return format
end

function show(io::IO, a::perm)
   if perm_display_style.format == :array
      print(io, "[" * join(a.d, ", ") * "]")
   elseif perm_display_style.format == :cycles
      if a == parent(a)()
         print(io, "()")
      else
         print(io, join(["("*join(c, ",")*")" for c in cycles(a) if
            length(c)>1],""))
      end
   end
end

###############################################################################
#
#   Comparison
#
###############################################################################

doc"""
    ==(a::perm, b::perm)
> Return `true` if the given permutations are equal, otherwise return `false`.
"""
function ==(a::perm, b::perm)
   parent(a) == parent(b) || return false
   a.d == b.d || return false
   return true
end

###############################################################################
#
#   Binary operators
#
###############################################################################

doc"""
    *(a::perm, b::perm)
> Return the composition of the two permutations, i.e. $a\circ b$. In other
> words, the permutation corresponding to applying $b$ first, then $a$, is
> returned.
"""
function *(a::perm, b::perm)
   d = similar(a.d)
   for i in 1:length(d)
      d[i] = a[b[i]]
   end
   return parent(a)(d, false)
end

doc"""
    ^(a::perm, n::Int)
> Return the `n`-th power of a permutation `a`. By default Nemo computes powers
> by cycle decomposition of `a`. `Nemo.power_by_squaring` provides a different
> method for powering which may or may not be faster, depending on the
> particuar case. Due to caching of cycle structure repeatedly powering should
> be faster with the default method. NOTE: cycle structure is not computed
> for `n<4`.
"""
function ^(a::perm, n::Int)
   if n <0
      return inv(a)^-n
   elseif n == 0
      return parent(a)()
   elseif n == 1
      return deepcopy(a)
   elseif n == 2
      return parent(a)(a.d[a.d])
   elseif n == 3
      return parent(a)(a.d[a.d[a.d]])
   else
      new_perm = similar(a.d)
      cyls = cycles(a)
      for cycle in cyls
         k = n % length(cycle)
         shifted = circshift(cycle, -k)
         for (idx,j) in enumerate(cycle)
            new_perm[j] = shifted[idx]
         end
      end
      return parent(a)(new_perm, false)
   end
end

function power_by_squaring(a::perm, b::Int)
   if n <0
      return inv(a)^-n
   elseif n == 0
      return parent(a)()
   elseif n == 1
      return deepcopy(a)
   elseif n == 2
      return parent(a)(a.d[a.d])
   elseif n == 3
      return parent(a)(a.d[a.d[a.d]])
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      cache1 = deepcopy(a.d)
      cache2 = deepcopy(a.d)
      bit >>= 1
      while bit != 0
         cache2 = cache1[cache1]
         cache1 = cache2
         if (UInt(bit) & b) != 0
            cache1 = cache1[a.d]
         end
         bit >>= 1
      end
      return parent(a)(cache1, false)
   end
end

###############################################################################
#
#   Inversion
#
###############################################################################

doc"""
    inv(a::perm)
> Return the inverse of the given permutation, i.e. the permuation $a^{-1}$
> such that $a\circ a^{-1} = a^{-1}\circ a$ is the identity permutation.
"""
function inv(a::perm)
   d = similar(a.d)
   for i in 1:length(a.d)
      d[a[i]] = i
   end
   G = parent(a)
   return G(d, false)
end

# TODO: can we do that in place??
function inv!(a::perm)
   d = similar(a.d)
   for i in 1:length(a.d)
      d[a[i]] = i
   end
   a.d = d
end

###############################################################################
#
#   Iterating over all permutations
#
###############################################################################

Base.start(A::AllPerms) = (collect(1:A.n), 1, 1, ones(Int, A.n))
Base.next(A::AllPerms, state) = all_perms(state...)
Base.done(A::AllPerms, state) = state[2] > A.all
Base.eltype(::Type{AllPerms}) = Vector{Int}
length(A::AllPerms) = A.all

function all_perms(elts, counter, i, c)
   if counter == 1
      return (copy(elts), (elts, counter+1, i, c))
   end
   n = length(elts)
   while i <= n
      if c[i] < i
         if isodd(i)
            elts[1], elts[i] = elts[i], elts[1]
         else
            elts[c[i]], elts[i] = elts[i], elts[c[i]]
         end
         c[i] += 1
         i = 1
         return (copy(elts), (elts, counter+1, i, c))
      else
         c[i] = 1
         i += 1
      end
   end
end

doc"""
    elements(G::PermGroup)
> Returns an iterator over all elements in the group $G$. You may use
> `collect(elements(G))` to get an array of all elements.
"""
elements(G::PermGroup) = (G(p, false) for p in AllPerms(G.n))

###############################################################################
#
#   Misc
#
###############################################################################

doc"""
    order(G::PermGroup)
> Returns the order of the full permutation group.
"""
order(G::PermGroup) = factorial(G.n)

doc"""
    cycles(a::perm)
> Decomposes permutation into disjoint cycles.
"""
function cycles(a::perm)
   if isdefined(a, :cycles)
      return a.cycles
   else
      to_visit = trues(a.d)
      cycles = Vector{Vector{Int}}()
      k = 1
      while any(to_visit)
         cycle = Vector{Int}()
         k = findnext(to_visit, k)
         to_visit[k] = false
         push!(cycle, k)
         next = a[k]
         while next != k
            push!(cycle, next)
            to_visit[next] = false
            next = a[next]
         end
         push!(cycles, cycle)
      end
      a.cycles = cycles
      return cycles
   end
end

doc"""
    order(a::perm)
> Returns the order of permutation `a`.
"""
order(a::perm) = lcm([length(c) for c in cycles(a)])

doc"""
    permtype(a::perm, rev=true)
> Returns the type of permutation `a`, i.e. lengths of disjoint cycles building
> `a`. This fully determines the conjugacy class of `a`. The lengths are sorted
> in reverse order by default.
"""
permtype(a::perm) = sort([length(c) for c in cycles(a)], rev=true)

doc"""
    matrix_repr(a::perm)
> Return the permutation matrix representing `a` via natural embedding of the
> permutation group into general linear group over ZZ
"""
function matrix_repr(a::perm)
   A = eye(Int, parent(a).n)
   return A[a.d,:]
end

doc"""
    emb!(result::perm, p::perm, V::Vector{Int})
> Embedds permutation `p` into `result` on the indices given by `V`. This
> corresponds to natural embedding of $S_k$ into $S_n$ as the subgroup
> permuting points indexed by `V`.
"""
function emb!(result::perm, p::perm, V::Vector{Int})
    result.d[V] = (result.d[V])[p.d]
    return result
end

doc"""
    emb(G::PermGroup, V::Vector{Int})
> Returns the natural embedding of a permutation group into `G` as the
> subgroup permuting points indexed by `V`.
"""
function emb(G::PermGroup, V::Vector{Int}, check::Bool=true)
   if check
      @assert length(Base.Set(V)) == length(V)
      @assert all(V .<= G.n)
   end
   return p -> Nemo.emb!(G(), p, V)
end

doc"""
    rand(G::PermGroup)
> Returns a random element from group `G`.
"""
rand(G::PermGroup) = G(randperm(G.n), false)

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (G::PermGroup)()
   z = perm(G.n)
   z.parent = G
   return z
end

function (G::PermGroup)(a::Array{Int, 1}, check::Bool=true)
   length(a) != G.n && error("Unable to coerce to permutation: lengths differ")
   if check
      Base.Set(a) != Base.Set(1:length(a)) && error("Unable to coerce to
         permutation: non-unique elements in array")
   end
   z = perm(a)
   z.parent = G
   return z
end

function (G::PermGroup)(a::perm)
   parent(a) != G && error("Unable to coerce to permutation: wrong parent")
   return a
end

###############################################################################
#
#   PermGroup constructor
#
###############################################################################

# handled by inner constructors

##############################################################################
#
#   Irreducible Characters
#
##############################################################################

doc"""
    character(lambda::Partition)
> Returns the $\lambda$-th irreducible character of permutation group on
> `lambda.n` generators. The returned character function
>     `chi(p::perm, check::Bool=true)`
> can be evaluated on a permutation `p::perm`. Checking in $\chi$ if `p`
> belongs to the appropriate group can be switched off by calling
> `chi(p, false)`. The values computed by $\chi$ are stored in look-up table.
>
> The computation follows the Murnaghan-Nakayama formula:
> $$\chi_\lambda(\sigma) = \sum_{\text{rimhook }\xi\subset \lambda}
(-1)^{ll(\lambda\backslash\xi)} \chi_{\lambda \backslash\xi}(\tilde\sigma)$$
> where $\lambda\backslash\xi$ denotes partition associated with Young Diagram
> of $\lambda$ with $\xi$ removed, $ll$ denotes the leg-length (i.e. number of
> rows - 1) and $\tilde\sigma$ is permutation obtained from $\sigma$ by the
> removal of the longest cycle.
>
> For more details see e.g. Chapter 2.8 of _Group Theory and Physics_ by
> S.Sternberg.
"""
function character(lambda::Partition)
   R = partitionseq(lambda)

   function char(p::perm, check::Bool=true)
      if check
         lambda.n == length(p.d) || throw("Can't evaluate character on $p : lengths differ.")
      end
      return MN1inner(R, Partition(permtype(p)), 1, _charvalsTable)
   end

   return char
end

doc"""
    character(lambda::Partition, p::perm, check::Bool=true)
> Returns the value of `lambda`-th irreducible character on permutation `p`.
"""
function character(lambda::Partition, p::perm, check::Bool=true)
   if check
      lambda.n == length(p.d) || throw("lambda-th irreducible character can be evaluated only on permutations of length $(lambda.n).")
   end

   return MN1inner(partitionseq(lambda), Partition(permtype(p)), 1, _charvalsTable)
end

doc"""
    character(lambda::Partition, mu::Partition)
> Returns the value of `lambda-th` irreducible character on the conjugacy class
> represented by partition `mu`. Values of characters computed by this method
> are not cached _globally_.
"""
character(lambda::Partition, mu::Partition) = MN1inner(partitionseq(lambda), mu, 1)


text:
# User interaction with RBTree.

# (c) E. Adrian Garro S. Costa Rica Institute of Technology.

include("tree.jl")

type UI
	option::String
	tree_num::Tree{Int64, Int64}
	tree_str::Tree{String, String}
end

"""
`input(prompt::AbstractString="")`

Read a string from STDIN. The trailing newline is stripped.

The prompt string, if given, is printed to standard output without a
trailing newline before reading input.
"""
function input(prompt::AbstractString="")
	print(prompt)
	return chomp(readline())
end

"""
`int_input(prompt::AbstractString="")`

Use input function to get integers from user.
"""
function int_input(string_num::AbstractString)
	possible_int = input(string_num)
	parse(Int, possible_int)
end

function read_option(app::UI)
	println("")
	println("--------(Red Black Tree)---------")
	println("     1) Clear tree numbers.")
	println("     2) Clear tree strings.")
	println("     3) Insert number.")
	println("     4) Insert string.")
	println("     5) Delete number.")
	println("     6) Delete string.")
	println("     7) Search number.")
	println("     8) Search string.")
	println("     9) Print inorder tree numbers.")
	println("     10) Print inorder tree strings.")
	println("     11) Pretty print tree numbers.")
	println("     12) Pretty print tree strings.")
	println("     0) End.")
	println("----------------------------------")
	println("")
	app.option = input("Select an option to continue: ")
end

# Option 1
function clear_tree_num(app::UI)
	clear!(app.tree_num)
	println("The tree has been successfully cleaned!")
end

# Option 2
function clear_tree_str(app::UI)
	clear!(app.tree_str)
	println("The tree has been successfully cleaned!")
end

# Option 3
function insert_num(app::UI)
	try
		num = int_input("Please enter the number you want insert: ")
		insert!(app.tree_num, num, num)
		println("The number has been successfully inserted!")
	catch e
		println(e)
	end
end

# Option 4
function insert_str(app::UI)
	string = input("Please enter the string you want insert: ")
	insert!(app.tree_str, string, string)
	println("The string has been successfully inserted!")
end

# Option 5
function delete_num(app::UI)
	try
		# try read number.
		num_key = int_input("Please enter the number you want delete: ")
		try
			# try delete number.
			delete!(app.tree_num, num_key)
			println("The number has been successfully deleted!")
		catch e
			# case where number key doesn't exist.
			println(e)
		end
	catch e
		# case where user doesn't type a number.
		println(e)
	end
end

# Option 6
function delete_str(app::UI)
	str_key = input("Please enter the string you want delete: ")
	try
		# try delete string.
		delete!(app.tree_str, str_key)
		println("The string has been successfully deleted!")
	catch e
		# case where string key doesn't exist.
		println(e)
	end
end

# Option 7
function search_num(app::UI)
	try
		# try read number.
		num_key = int_input("Please enter the number you want search: ")
		try
			# try search number.
			num_val = search(app.tree_num, num_key)
			println(string("The number ", string(string(num_val), " was found!")))
		catch e
			# case where number doesn't exist.
			println(e)
		end
	catch e
		# case where user doesn't type a number.
		println(e)
	end
end

# Option 8
function search_str(app::UI)
	str_key = input("Please enter the string you want search: ")
	try
		# try search string.
		str_val = search(app.tree_str, str_key)
		println(string("The string ", string(string(str_val), " was found!")))
	catch e
		# case where string doesn't exist.
		println(e)
	end
end

# Option 9
function print_inorder_num(app::UI)
	for num in @task inorderwalk(app.tree_num)
		println(num)
	end
end

# Option 10
function print_inorder_str(app::UI)
	for string in @task inorderwalk(app.tree_str)
		println(string)
	end
end

# Option 11
function pretty_print_num(app::UI)
	prettyprint(app.tree_num)
end

# Option 12
function pretty_print_str(app::UI)
	prettyprint(app.tree_str)
end

function menu(app::UI)
	if app.option == ""
		true
	elseif app.option == "1"
		clear_tree_num(app::UI)
		true
	elseif app.option == "2"
		clear_tree_str(app::UI)
		true
	elseif app.option == "3"
		insert_num(app::UI)
		true
	elseif app.option == "4"
		insert_str(app::UI)
		true
	elseif app.option == "5"
		delete_num(app::UI)
		true
	elseif app.option == "6"
		delete_str(app::UI)
		true
	elseif app.option == "7"
		search_num(app::UI)
		true
	elseif app.option == "8"
		search_str(app::UI)
		true
	elseif app.option == "9"
		print_inorder_num(app::UI)
		true
	elseif app.option == "10"
		print_inorder_str(app::UI)
		true
	elseif app.option == "11"
		pretty_print_num(app::UI)
		true
	elseif app.option == "12"
		pretty_print_str(app::UI)
		true
	elseif app.option == "0"
		println("Thank you for your time. ;)")
		false
	else
		println("Pick a valid option from menu!")
		true
	end
end

function main()
	app = UI("", Tree{Int64, Int64}(), Tree{String, String}())
	while menu(app::UI)
		read_option(app::UI)
	end
end

main()

max_stars_repo_name: adriangarro/llrb

max_issues_repo_name: adriangarro/llrb

max_forks_repo_name: adriangarro/llrb


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Libdl

name = "libbitfield"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    DirectorySource(@__DIR__),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cmake -B build -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
cmake --build build
cmake --install build
touch LICENSE
"""

# Only need to build for host platform since we are not publishing
platforms = [HostPlatform()]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libbitfield"], :libbitfield)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
const ContactImplicitMPC = Main
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
s_sim = get_simulation("hopper_2D", "sine2_2D_lc", "sinusoidal")
s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
model = s.model
env = s.env

# get trajectory
ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_forward.jld2"),
    load_type = :joint_traj))
# time
H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 100*H*N_sample #500*H*N_sample

# Select the initial speed
v0 = 0.2
Tstance = 0.13 # measure using hop-in-place gait
Tflight = 0.62 # measure using hop-in-place gait
p = raibert_policy(s_sim.model, v0=v0, Tstance=Tstance, Tflight=Tflight, h=h)

off0 = SVector{model.dim.q,T}([0.0, 0.5, 0.0, 0.0])
off1 = SVector{model.dim.q,T}([0*v0*h_sim, 0.5, 0.0, 0.0])
q_ref = SVector{model.dim.q,T}([0.0, 0.5, 0.0, 0.5])
q0_sim = copy(q_ref) + off0
q1_sim = copy(q_ref) + off1

sim = ContactImplicitMPC.simulator(s_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = ContactImplicitMPC.InteriorPointOptions(
		γ_reg = 0.0,
		undercut = Inf,
		r_tol = 1.0e-8,
		κ_tol = 1.0e-8,),
    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

status = ContactImplicitMPC.simulate!(sim, verbose = true)

plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.(sim.traj.q)...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[1:model.dim.u] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
plot!(plt[3,1], hcat(Vector.([γ[1:model.dim.c] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)

plot_lines!(vis, model, sim.traj.q[1:10:end])
plot_surface!(vis, s_sim.env, n=200, xlims = [-1, 40])
anim = visualize_robot!(vis, model, sim.traj, sample=5)
anim = visualize_force!(vis, model, s_sim.env, sim.traj, anim=anim, h=h_sim, sample = 5)

plot(hcat([u[[1,2]] for u in Vector.(sim.traj.u)[1:end]]...)')

# filename = "hopper_ipopt"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)

# const ContactImplicitMPC = Main


text:
type Thermostat{F<:Function}
    name::String
    distribution::F
end

type Potential{F<:Function}
    name::String
    f::F
end

type Integrator{F<:Function}
    name::String
    f::F
end

type Parameters
    T::Float64
    Q::Float64
    dtsampling::Float64
    dt::Float64
    nsimulations::Int64
    nsteps::Int64
    thermo::Thermostat
    potential::Potential
    integrator::Integrator
end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
module Common

using StaticArrays

export Float, Coordinate, Coordinate2, AngMomArray

Float = Float64

"""
    Coordinate(x::NTuple{3, Float64})
    Coordinate(x::AbstractArray{3, Float64})
    Coordinate(x1, x2, x3)

Represents the cartesian coordinates of a point in 3D space.
"""
struct Coordinate <: StaticArray{Tuple{3}, Float, 1}
    coords::NTuple{3, Float}
    Coordinate(coords::NTuple{3, Float}) = new(coords)
    Coordinate(coords::NTuple{3, T}) where {T} = throw(TypeError(:Coordinate, "coords", Type{Float}, T))
end
Base.getindex(R::Coordinate, i::Int) = R.coords[i]
Tuple(R::Coordinate) = R.coords

#  struct AngMomArray{L, S, T}
    #  data::NTuple{S, T}
#  end

#  function AngMomArray{L, T}(data) where {L<:Tuple, T}
    #  S = (L.parameters[1]+1)^2
    #  AngMomArray{L, S, T}(data)
#  end
#  Base.getindex(a::AngMomArray{L, S, T}, l::Int, m::Int) where {L, S, T} = a.data[l^2+m+l+1]

end


text:
module ForecastData

using Pkg.Artifacts
using CSV, DataFrames, DataFramesMeta, Dates, Downloads, Inflate, SQLite, Tar

include("art_path.jl")

# Datasets
include("co2.jl")
export co2

include("seaborne.jl")
export seaborne

include("quakes.jl")
export quakes

include("london.jl")
export london

"""
List of available dataset functions:

| Name | Description |
|---|---
|**co2** | Atmospheric carbon dioxide
|**london** | Weather and crime in Greater London
|**quakes** | Earthquakes around the world
|**seaborne** | World seaborne trade
"""
ForecastData

end


text:
"""
    cgls(op, b; <keyword arguments>)
CGLS algorithms for operators whose input and output are files saved on disk

# Arguments
- `op::Function`: linear operator whose forward is called as op(d, m, 1; op_params...) and the adjoint op(m, d, 2; op_params...)
- `b ::String`  : path to oberved data.

# Keyword arguments
- `op_params::NamedTuple=NamedTuple()`: the parameters for the keyword arguments of the linear operator `op`.
- `d_axpby::Function`: compute `y = a*x + b*y` for variables in data space.
- `m_axpby::Function`: compute `y = a*x + b*y` for variables in model space.
- `d_norm::Function` : compute the L2 norm for variables in data space.
- `m_norm::Function` : compute the L2 norm for variables in model space.
- `x0::String="NULL"`      : initial guess of the model parameters.
- `dir_work::String="NULL"`: working directory for the files generated by the linear operator, the directory of pwd() is used if not provided.
- `shift      = 0.0` : trade-off parameter.
- `tol        = 1e-6`: tolerence for convergence.
- `maxIter    = 50`  : maximum number of iterations.
- `print_flag = true`: print the running information.
- `save_flag  = true`: save the intermidiate result of model parameters.

# Improvement
- return the history of iterations
"""
function cgls(op::Function, b::Ts; dir_work::Ts="NULL", op_params::NamedTuple=NamedTuple(),
              d_axpby::Function, m_axpby::Function,  d_norm::Function, m_norm::Function,
              x0::Ts="NULL", shift=0.0, tol=1e-6, maxIter=50, print_flag=true, save_flag=true) where {Ts<:String}

    # determine the working directory2`
    if dir_work == "NULL"
       dir_work =  pwd()
    end

    # create directory for variables in model space
    dir_model = joinpath(dir_work, "cgls_model_space")
    rm(dir_model, force=true, recursive=true)               # clean space
    mkdir(dir_model)
    if !isdir(dir_model)
       error("could not create directory for variables in model space")
    end
    x = joinpath(dir_model, "x.rsf")
    s = joinpath(dir_model, "s.rsf")
    p = joinpath(dir_model, "p.rsf")

    # create directory for variables in data space
    dir_data = joinpath(dir_work, "cgls_data_space")
    rm(dir_data, force=true, recursive=true)                # clean space
    mkdir(dir_data)
    r = joinpath(dir_data, "r"); mkdir(r);
    q = joinpath(dir_data, "q"); mkdir(q);
    if !isdir(dir_data) || !isdir(r) || !isdir(q)
       error("could not create directory of residue")
    end

    # create a directory save the intermidiate result
    if save_flag
       dir_iterations = joinpath(dir_work, "iterations")
       rm(dir_iterations, force=true, recursive=true)        # clean space
       mkdir(dir_iterations)
       if !isdir(dir_iterations)
          error("could not create directory for iterations")
       end
    end

    # initial guess provided
    if x0 != "NULL"

       cp(x0, x, force=true)
       op(r, x, 1; op_params...)
       d_axpby(1.0, b, -1.0, r)                            # r = b - A*x

       op(s, r, 2; op_params...)
       m_axpby(-shift, x, 1.0, s)                          # s = A'*r - shift*x

    # no initial guess
    else

       cp(b, r, force=true)                                # r = b
       op(s, r, 2; op_params...)                           # s = A'*r

       # create x
       hdr = read_RSheader(s)
       write_RSdata(x, hdr, zeros(hdr))                    # (x=0)
    end

    # compute residue
    data_fitting = (d_norm(r))^2                           # ||b-Ax||_2^2
    constraint   = 0.0
    cost0        = data_fitting + constraint
    convergence  = Float64[]; push!(convergence, 1.0);

    # initialize some intermidiate vectors
    cp(s, p, force=true)                                   # p = s

    norms0= m_norm(s)                                      # ||s||_2
    gamma = norms0^2                                       # ||s||_2^2
    normx = m_norm(x)                                      # ||x||_2
    xmax  = normx     # keep the initial one               # ||x||_2
    resNE = 1.0

    gamma0= copy(gamma)                                    # ||s||_2^2
    delta = 0.0

    # iteration counter and stop condition
    k = 0
    run_flag = true

    if print_flag
       header = "  k         data_fitting           constraint             normx                resNE"
       println(""); println(header);
       @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
    end

    while (k < maxIter) && run_flag

          k = k + 1
          op(q, p, 1; op_params...)                        # q = A * p

          delta = (d_norm(q))^2 + shift * (m_norm(p))^2    # ||q||_2^2 + shift*||p||_2^2
          indefinite = delta <= 0.0 ? true  : false
          delta      = delta == 0.0 ? eps() : delta

          alpha = gamma / delta

          m_axpby(alpha, p, 1.0, x)                        # x = x + alpha * p
          d_axpby(-alpha, q, 1.0, r)                       # r = r - alpha * q

          data_fitting = (d_norm(r))^2                     # ||r||_2^2
          constraint   = shift * (m_norm(x))^2             # shift*||x||_2^2
          cost         = data_fitting + constraint         # ||r||_2^2 + shift*||x||_2^2

          # save the intermidiate result
          path_iter = join([dir_iterations "/iteration" "_" "$k" ".rsf"])
          cp(x, path_iter, force=true)

          op(s, r, 2; op_params...)
          m_axpby(-shift, x, 1.0, s)                       # s = A' * r - shift * x

          norms  = m_norm(s)                               # ||s||_2
          gamma0 = gamma                                   # ||s||_2^2 previous iteration
          gamma  = norms^2                                 # ||s||_2^2
          beta   = gamma / gamma0

          m_axpby(1.0, s, beta, p)                         # p = s + beta * p

          # check the stopping crietia
          normx = m_norm(x)
          xmax  = normx > xmax ? normx : xmax
          if norms <= norms0 * tol || normx * tol >= 1.0
             run_flag = false
          end

          # print information
          resNE = cost / cost0
          if print_flag
             @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
          end
          push!(convergence, resNE);
    end

    return x, convergence
end


text:
using Random, Plots, Distributions, Statistics

# Example solution for example 6.5, exercise 6.9 and 6.10 using Julia 1.4.1
#
# By @burmecia at GitHub, 2nd May 2020
#

α, ε, γ = 0.5, 0.1, 1

# grid world map, origin (1,1) is at the bottom left
rows, cols = 7, 10
start, goal = [4, 1], [4, 8]

# full action list
actions = ("N", "E", "S", "W", "NE", "SE", "SW", "NW", "Z")

# wind
wind = zeros(Int64, cols)
wind[4:9] .= 1
wind[7:8] .= 2
wind_offset = 0:0

# global Q array
Q = Array{Float64}

episode_num = 300

function next_state(s, a)
    s′ = copy(s)
    if a == 1
        s′[1] += 1  # north
    elseif a == 2
        s′[2] += 1  # east
    elseif a == 3
        s′[1] -= 1  # south
    elseif a == 4
        s′[2] -= 1  # west
    elseif a == 5
        s′[1] += 1  # north east
        s′[2] += 1
    elseif a == 6
        s′[1] -= 1  # south east
        s′[2] += 1
    elseif a == 7
        s′[1] -= 1  # south west
        s′[2] -= 1
    elseif a == 8
        s′[1] += 1  # north west
        s′[2] -= 1
    elseif a == 9
        # do nothing
    end

    # row shift by wind
    w = wind[clamp(s[2], 1, length(wind))]
    s′[1] += w + (w > 0 ? rand(wind_offset) : 0)

    s′
end

# choose action using ε-greedy policy
function choose_action(s, act_len; ε = ε)
    if rand() >= ε
        argmax(Q[s[1], s[2], :])
    else
        rand(1:act_len)
    end
end

# take action, observe next state and get next action using ε-greedy policy
function next_state_action(s, a, act_len; ε = ε)
    # take action, get next state
    s′ = next_state(s, a)

    # if agent is off the grid, leave its position unchanged
    s′[1] = clamp(s′[1], 1, rows)
    s′[2] = clamp(s′[2], 1, cols)

    # choose A′ from S′
    a′ = choose_action(s′, act_len, ε = ε)

    (s′, a′)
end

# output an optimal tracjectory using greedy policy
function output_optimal_tractory()
    s, a = start, choose_action(start, act_len, ε = 0)
    while s != goal
        println("S: $(s), A: $(actions[a])")
        s, a = next_state_action(s, a, act_len, ε = 0)
    end
end

function game(act_len)
    global Q

    # initialise Q
    Q = zeros(Float64, rows, cols, act_len)

    hist = zeros(Int64, episode_num)
    r = -1

    for i = 1:episode_num
        step = 0
        s, a = start, choose_action(start, act_len)
        while s != goal
            s′, a′ = next_state_action(s, a, act_len)

            # upate Q[sa]
            sa = CartesianIndex(Tuple([s..., a]))
            s′a′ = CartesianIndex(Tuple([s′..., a′]))
            Q[sa] += α * (r + γ * Q[s′a′] - Q[sa])

            s, a = s′, a′
            step += 1
        end

        # cap extream values in the first several episodes
        hist[i] = min(step, 300)
    end

    reshape(hist, 1, :)
end

function main()
    global wind_offset

    hist_4 = Array{Float64,2}(undef, 0, episode_num)
    hist_8 = Array{Float64,2}(undef, 0, episode_num)
    hist_9 = Array{Float64,2}(undef, 0, episode_num)

    for i = 1:200
        hist_4 = vcat(hist_4, game(4))
        hist_8 = vcat(hist_8, game(8))
        hist_9 = vcat(hist_9, game(9))
    end

    # exercise 6.10
    wind_offset = -1:1
    hist_8_rand_wind = Array{Float64,2}(undef, 0, episode_num)
    for i = 1:200
        hist_8_rand_wind = vcat(hist_8_rand_wind, game(8))
    end

    hist_4 = mean(hist_4, dims = 1)
    hist_8 = mean(hist_8, dims = 1)
    hist_9 = mean(hist_9, dims = 1)
    hist_8_rand_wind = mean(hist_8_rand_wind, dims = 1)

    p = plot(
        title = "4 vs 8 vs 9 action space with random wind",
        xlabel = "Episodes (averaged over 200 experiments)",
        ylabel = "Steps",
        legend = :topright,
    )
    plot!(p, hist_4[1, :], label = "actions=4")
    plot!(p, hist_8[1, :], label = "actions=8")
    plot!(p, hist_9[1, :], label = "actions=9")
    plot!(p, hist_8_rand_wind[1, :], label = "actions=8, random wind")
end

# uncomment the main function below to start learning
#main()


text:
using Test, SummationByPartsOperators
using LinearAlgebra
using SpecialFunctions

function accuracy_test!(res, ufunc, dufunc, D)
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    mul!(res, D, u)
    maximum(abs, du-res) < 5*length(res)*eps(eltype(res))
end

# Accuracy Tests
for T in (Float32, Float64)
    xmin = -one(T)
    xmax = one(T)

    for N in 2 .^ (1:4)
        D = legendre_derivative_operator(xmin, xmax, N)
        println(devnull, D)
        @test SummationByPartsOperators.derivative_order(D) == 1
        @test issymmetric(D) == false
        u = compute_coefficients(zero, D)
        res = D*u
        for k in 1:N-1
            ufunc = x -> x^k/typeof(x)(gamma(k+1))
            dufunc = x -> x^(k-1)/typeof(x)(gamma(k))
            @test accuracy_test!(res, ufunc, dufunc, D)
            xplot, duplot = evaluate_coefficients(res, D)
            @test maximum(abs, duplot - dufunc.(xplot)) < 5N*eps(T)
            @test abs(integrate(u, D)) < N*eps(T)
        end
    end
end

# Modal Filtering
for T in (Float32, Float64), filter_type in (ExponentialFilter(),)
    xmin = -one(T)
    xmax = one(T)

    for N in 2 .^ (1:4)
        D = legendre_derivative_operator(xmin, xmax, N)
        filter! = ConstantFilter(D, filter_type)
        u = compute_coefficients(zero, D)
        res = D*u
        for k in 1:N-1
            compute_coefficients!(u, x->exp(sinpi(x)), D)
            norm2_u = integrate(u->u^2, u, D)
            filter!(u)
            @test integrate(u->u^2, u, D) <= norm2_u
            @test integrate(u, D) ≈ sum(mass_matrix(D) * u)
            @test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u)
        end
    end
end


text:
"""
    skan!(repo::Repo, pages::Vector{<:Page}; notify=true) -> Vector{PageScan}

Compare the `pages` to the stored state in `repo`.
Return a vector of changed `pages` and update the repo for each changed page.
When no pages changed, the vector is empty.
"""
function skan!(repo, pages::Vector{<:Page}; notify=true)::Vector{PageScan}
    pull_or_clone!(repo)
    state = retrieve(repo)

    changed = map(pages) do page
        oldscan = retrieve(state, page)::Union{PageScan,Nothing}
        newscan = scan(page)::PageScan
        if isnothing(oldscan) || oldscan.content != newscan.content
            return BeforeAfter(oldscan, newscan)
        else
            return nothing
        end
    end
    filter!(!isnothing, changed)
    # Without this, filtered may become a `Vector{Union{Nothing, Skans.BeforeAfter}}`.
    filtered = convert(Array{BeforeAfter}, changed)

    new_pages = getproperty.(filtered, :new)
    update!(repo, state, new_pages)

    if notify && !isempty(filtered)
        post_issue_comment!(repo, filtered)
    end

    return new_pages
end

function mail(changed::Vector{PageScan})
    error("Not implemented yet")
end


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

@doc raw"""A node selector requirement is a selector that contains values, a key, and an operator that relates the key and values.

    IoK8sApiCoreV1NodeSelectorRequirement(;
        key=nothing,
        operator=nothing,
        values=nothing,
    )

    - key::String : The label key that the selector applies to.
    - operator::String : Represents a key's relationship to a set of values. Valid operators are In, NotIn, Exists, DoesNotExist. Gt, and Lt.
    - values::Vector{String} : An array of string values. If the operator is In or NotIn, the values array must be non-empty. If the operator is Exists or DoesNotExist, the values array must be empty. If the operator is Gt or Lt, the values array must have a single element, which will be interpreted as an integer. This array is replaced during a strategic merge patch.
"""
mutable struct IoK8sApiCoreV1NodeSelectorRequirement <: SwaggerModel
    key::Any # spec type: Union{ Nothing, String } # spec name: key
    operator::Any # spec type: Union{ Nothing, String } # spec name: operator
    values::Any # spec type: Union{ Nothing, Vector{String} } # spec name: values

    function IoK8sApiCoreV1NodeSelectorRequirement(;key=nothing, operator=nothing, values=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeSelectorRequirement, Symbol("key"), key)
        setfield!(o, Symbol("key"), key)
        validate_property(IoK8sApiCoreV1NodeSelectorRequirement, Symbol("operator"), operator)
        setfield!(o, Symbol("operator"), operator)
        validate_property(IoK8sApiCoreV1NodeSelectorRequirement, Symbol("values"), values)
        setfield!(o, Symbol("values"), values)
        o
    end
end # type IoK8sApiCoreV1NodeSelectorRequirement

const _property_map_IoK8sApiCoreV1NodeSelectorRequirement = Dict{Symbol,Symbol}(Symbol("key")=>Symbol("key"), Symbol("operator")=>Symbol("operator"), Symbol("values")=>Symbol("values"))
const _property_types_IoK8sApiCoreV1NodeSelectorRequirement = Dict{Symbol,String}(Symbol("key")=>"String", Symbol("operator")=>"String", Symbol("values")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeSelectorRequirement }) = collect(keys(_property_map_IoK8sApiCoreV1NodeSelectorRequirement))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeSelectorRequirement }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeSelectorRequirement[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeSelectorRequirement }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeSelectorRequirement[property_name]

function check_required(o::IoK8sApiCoreV1NodeSelectorRequirement)
    (getproperty(o, Symbol("key")) === nothing) && (return false)
    (getproperty(o, Symbol("operator")) === nothing) && (return false)
    true
end

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


text:
using Plots, BenchmarkTools
Threads.nthreads()
function juliaset_pixel(z₀, c)
    z = z₀
    for i in 1:255
        abs2(z)> 4.0 && return (i - 1)%UInt8
        z = z*z + c
    end
    return UInt8(255)
end

function juliaset_column!(img, c, n, j)
    x = -2.0 + (j-1)*4.0/(n-1)
    for i in 1:n
        y = -2.0 + (i-1)*4.0/(n-1)
        @inbounds img[i,j] = juliaset_pixel(x+im*y, c)
    end
    nothing
end

function juliaset_single!(img, c, n)
    for j in 1:n
        juliaset_column!(img, c, n, j)
    end
    nothing
end

function juliaset(x, y, n=1000, method = juliaset_single!, extra...)
    c = x + y*im
    img = Array{UInt8,2}(undef,n,n)
    method(img, c, n, extra...)
    return img
end

# frac = juliaset(-0.79, 0.15)
# plot(heatmap(1:size(frac,1),1:size(frac,2), frac, color=:Spectral))

@btime juliaset(-0.79, 0.15);

function juliaset_static!(img, c, n)
    Threads.@threads for j in 1:n
        juliaset_column!(img, c, n, j)
    end
    nothing
end

@btime juliaset(-0.79, 0.15, 1000, juliaset_static!);

using Folds
function juliaset_folds!(img, c, n)
    Folds.foreach(j -> juliaset_column!(img, c, n, j), 1:n)
    nothing
end
julia> @btime juliaset(-0.79, 0.15, 1000, juliaset_folds!);
  16.267 ms (25 allocations: 978.20 KiB)

function juliaset_folds!(img, c, n, nt)
    parts = collect(Iterators.partition(1:n, cld(n, nt)))
    Folds.foreach(parts) do ii
        foreach(j ->juliaset_column!(img, c, n, j), ii)
    end
    nothing
end
julia> @btime juliaset(-0.79, 0.15, 1000, (args...) -> juliaset_folds!(args..., 16));
  16.716 ms (25 allocations: 978.61 KiB)

using FLoops, FoldsThreads
function juliaset_folds!(img, c, n)
    @floop ThreadedEx(basesize = 2) for j in 1:n
        juliaset_column!(img, c, n, j)
    end
    nothing
end
@btime juliaset(-0.79, 0.15, 1000, juliaset_folds!);


text:
export CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION

function convolution_desc(T::Type, padding, stride, mode)
    N = length(padding)
    p = Ptr{Void}[0]
    cudnnCreateConvolutionDescriptor(p)
    cpadding = Cint[padding[i] for i=N:-1:1]
    cstrides = Cint[stride[i] for i=N:-1:1]
    cupscale = fill(Cint(1), N)
    cudnnSetConvolutionNdDescriptor(p[1], N, cpadding, cstrides, cupscale, mode, datatype(T))
    p[1]
end

function convolution{T}(x::CuArray{T}, w::CuArray{T}, padding, stride;
    mode=CUDNN_CROSS_CORRELATION, alpha=1.0, beta=0.0)

    N = length(padding)
    outdims = ntuple(i -> (size(x,i)+2padding[i]-size(w,i)) ÷ stride[i] + 1, N)
    y = similar(x, outdims..., size(w,N+2), size(x,N+2))

    h = handle(x)
    xdesc = tensor_desc(x)
    wdesc = filter_desc(w)
    convdesc = convolution_desc(T, padding, stride, mode)
    ydesc = tensor_desc(y)

    algo_p = cudnnConvolutionFwdAlgo_t[0]
    cudnnGetConvolutionForwardAlgorithm(h, xdesc, wdesc, convdesc, ydesc,
        CUDNN_CONVOLUTION_FWD_PREFER_FASTEST, 0, algo_p)
    algo = algo_p[1]

    worksize_p = Cint[0]
    cudnnGetConvolutionForwardWorkspaceSize(h, xdesc, wdesc, convdesc, ydesc, algo, worksize_p)
    worksize = worksize_p[1]
    workspace = CuArray{Int8}(Int(worksize))

    cudnnConvolutionForward(h, T[alpha], xdesc, x, wdesc, w, convdesc,
        algo, workspace, worksize, T[beta], ydesc, y)

    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyFilterDescriptor(wdesc)
    cudnnDestroyConvolutionDescriptor(convdesc)
    cudnnDestroyTensorDescriptor(ydesc)
    y
end

function ∇convolution_bias!(dy, db; mode=CUDNN_CROSS_CORRELATION, alpha=1.0, beta=0.0)
    T = eltype(dy)
    h = handle(dy)
    dydesc = tensor_desc(dy)
    dbdesc = tensor_desc(db)
    cudnnConvolutionBackwardBias(h, T[alpha], dydesc, dy, T[beta], dbdesc, db)

    cudnnDestroyTensorDescriptor(dydesc)
    cudnnDestroyTensorDescriptor(dbdesc)
    db
end

function ∇convolution_filter!(x, padding, stride, dy, dw; mode=CUDNN_CROSS_CORRELATION, alpha=1.0, beta=0.0)
    T = eltype(dy)
    h = handle(dy)
    xdesc = tensor_desc(x)
    dydesc = tensor_desc(dy)
    convdesc = convolution_desc(T, padding, stride, mode)
    dwdesc = filter_desc(dw)

    algo_p = cudnnConvolutionBwdFilterAlgo_t[0]
    cudnnGetConvolutionBackwardFilterAlgorithm(h, xdesc, dydesc, convdesc, dwdesc,
        CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST, 0, algo_p)
    algo = algo_p[1]

    worksize_p = Cint[0]
    cudnnGetConvolutionBackwardFilterWorkspaceSize(h, xdesc, dydesc, convdesc, dwdesc, algo, worksize_p)
    worksize = worksize_p[1]
    workspace = CuArray{Int8}(Int(worksize))

    cudnnConvolutionBackwardFilter(h, T[alpha], xdesc, x, dydesc, dy, convdesc,
        algo, workspace, worksize, T[beta], dwdesc, dw)

    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(dydesc)
    cudnnDestroyConvolutionDescriptor(convdesc)
    cudnnDestroyFilterDescriptor(dwdesc)
    dw
end

function ∇convolution_data!(w, padding, stride, dy, dx; mode=CUDNN_CROSS_CORRELATION, alpha=1.0, beta=0.0)
    T = eltype(dy)
    h = handle(dy)
    wdesc = filter_desc(w)
    dydesc = tensor_desc(dy)
    convdesc = convolution_desc(T, padding, stride, mode)
    dxdesc = tensor_desc(dx)

    algo_p = cudnnConvolutionBwdDataAlgo_t[0]
    cudnnGetConvolutionBackwardDataAlgorithm(h, wdesc, dydesc, convdesc, dxdesc,
        CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST, 0, algo_p)
    algo = algo_p[1]

    worksize_p = Cint[0]
    cudnnGetConvolutionBackwardDataWorkspaceSize(h, wdesc, dydesc, convdesc,
        dxdesc, algo, worksize_p)
    worksize = worksize_p[1]
    workspace = CuArray{Int8}(Int(worksize))

    cudnnConvolutionBackwardData(h, T[alpha], wdesc, w, dydesc, dy, convdesc,
        algo, workspace, worksize, T[beta], dxdesc, dx)

    cudnnDestroyFilterDescriptor(wdesc)
    cudnnDestroyTensorDescriptor(dydesc)
    cudnnDestroyConvolutionDescriptor(convdesc)
    cudnnDestroyTensorDescriptor(dxdesc)
    dx
end


text:
export Filters

module Filters

using ..Ranger
using ..Ranger.Nodes
using ..Ranger.Rendering
using ..Ranger.Math

include("transform_filter.jl")
include("translate_filter.jl")

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module CommonUtils

using Caesar
using Images
using FileIO
using Cairo
using RoMEPlotting

export plotSLAM2D_KeyAndSim, plotHMSLevel

@deprecate buildDEMSimulated(w...;kw...) RoME.generateField_CanyonDEM(w...;kw...)
@deprecate getSampleDEM(w...;kw...) RoME.generateField_CanyonDEM(w...;kw...)
@deprecate loadDEM!(w...;kw...) RoME._buildGraphScalarField(w...;kw...)

@deprecate plotSLAM2D_KeyAndSim(w...;kw...) RoMEPlotting.plotSLAM2D_KeyAndRef(w...;kw...)

## plot some of the ROIs

function plotHMSLevel(fg::AbstractDFG,
                      lbl::Symbol;
                      coord=Coord.cartesian(xmin=-9000, xmax=9000, ymin=-9000,ymax=9000)  )
  #
  loc = getPPE(fg, lbl, :simulated).suggested
  plp = plot( x=[loc[1]], y=[loc[2];],
              Geom.point,
              Guide.title(string(lbl)),
              Theme(default_color=colorant"purple") )
  #
  fct = intersect(ls(fg, lbl), lsf(fg, tags=[:DEM;]))[1]
  plk = getFactorType(fg[fct]).Z.densityFnc |> plotKDE;

  #
  union!(plp.layers, plk.layers);
  plp.coord = coord
  plp
end

end


text:
using DelimitedFiles, Plots

"""
`StrongField.plot( comps::Array{StrongField.Computation,1}, results::Array{Dict{String,Any},1}, energyScale::String = "atomic", probabilityScaling::String = "linear", dataLabel::String = "StrongFieldData" )`
    ... generates a graphical representation of the observable (StrongField.SfaEnergyDistribution, StrongField.SfaMomentumDistribution, StrongField.SfaAzimuthalAngularDistribution or StrongField.SfaPolarAngularDistribution) with the results of all computations comps combined in one plot.
    - All comps.observable need to be equal!
    - The probabilities for all comps need to be at the same grid points (either energies, momenta or angles)
    - energyScale determines the scaling of the energy axis either in atomic units (energyScale = "atomic") or in units of hbar*omega (energyScale = "omega").
    - The y-axis is scaled either linearly (probabilityScaling = "linear") or logarithmically (probabilityScaling = "log")
"""
function plot( comps::Array{StrongField.Computation,1}, results::Array{Dict{String,Any},1}, energyScale::String = "atomic", probabilityScaling::String = "linear", dataLabel::String = "StrongFieldData" )

    #---Check if all observables are equal---
    observable = typeof(comps[1].observable)
    for comp in comps
        if  typeof(comp.observable) != observable
            error("Unequal types of observables for data given to StrongField.plot().")
            return nothing
        end
    end

    #Define colors and linestyles for plot of SfaEnergyDistribution and SfaAngularDistribution
    colors = [:black,:blue,:red,:green,:purple]
    styles = filter((s->begin
            s in Plots.supported_styles()
        end), [:solid, :dash, :dot, :dashdot, :dashdotdot])

    #---Photoelectron energy spectra---
    if  observable == StrongField.SfaEnergyDistribution

        for k = 1:size(comps,1) #Loop over all comps

            #prepare data and rescale the x-axis if neccessary
            if  energyScale == "atomic"
                energyLabel = "E (a.u.)"
                energies = results[k]["energy distribution"].energies
            elseif  energyScale == "omega"
                energyLabel = "E/w"
                energies = results[k]["energy distribution"].energies / comps[k].beam.omega
            end
            probabilities = results[k]["energy distribution"].probabilities
            gr() #sets the plotting backend to the package "GR"

            #set scaling of the y-axis
            scaling = :identity
            if  probabilityScaling == "log"
                        scaling = :log10
            end

            #generate the plot
            if  k == 1 #initial plot generation
                Plots.plot(energies, probabilities,
                            title = "Photoelectron energy spectrum",
                            xlabel = energyLabel,
                            ylabel = "P(E)",
                            xscale = :identity,
                            yscale = scaling,
                            framestyle = :box,
                            legend = :none,
                            #markershape = :circle,
                            #markershape = :none,
                            line = (2,styles[1]),
                            linecolor = colors[1],
                            gridlinewidth = 2,
                            tickfontsize = 10,
                            labelfontsize = 10,
                            labelfontfamily = "Latin Modern Roman",
                            titlefontfamily = "Latin Modern Roman"
                            )
            else    #add to the plot already generated for k = 1
                Plots.plot!(energies, probabilities,
                            title = "Photoelectron energy spectrum",
                            xlabel = energyLabel,
                            ylabel = "P(E)",
                            xscale = :identity,
                            yscale = scaling,
                            framestyle = :box,
                            legend = :none,
                            #markershape = :circle,
                            #markershape = :none,
                            line = (2,styles[k]),
                            linecolor = colors[k],
                            gridlinewidth = 2,
                            tickfontsize = 10,
                            labelfontsize = 10,
                            labelfontfamily = "Latin Modern Roman",
                            titlefontfamily = "Latin Modern Roman"
                            )
            end

        end

        #export the plot as png-file
        savefig(dataLabel * "-energy_spectrum.pdf")
    #-----------------------------------

    #---Photoelectron momentum distributions - DOES NOT WORK YET---
    elseif  observable == StrongField.SfaMomentumDistribution
            #prepare data
            #pList = results["momentum distribution"].momenta
            #probList = results["momentum distribution"].probabilities

            #generate the plot
            #pyplot()
            #r = range(0,stop=10,length=11)
            #theta = range(0,stop=360,length=361)
            #f(r,theta) = r^2
            #println("$(f.(r,theta'))")
            #Plots.plot( heatmap( f.(r,theta'), proj=:polar ) )
            #println("$(transpose(probList))")
            #Plots.plot( heatmap( pList, transpose(probList), proj=:polar ) )
            #imshow(pList,transpose(probList))
            #probList = probList[2:-1]
            #matplotlib.pyplot.pcolormesh(pList[:,1], pList[:,2], transpose(probList))
            contourf( pList[:,1], pList[:,2], probList )
            #matplotlib.pyplot.imshow(pList)
            #heatmap( pList[:,1], pList[:,2], probList, proj=:polar, legend=true
                            #heatmap( probList
            #                )

            #p = Plots.plot( #heatmap( pList[:,1], pList[:,2], probList
            #               Plots.GR.polarheatmap( probList
                        ##heatmap( probList
            #                 )
            #              )

            #export the plot as png-file
            #savefig(dataLabel * "-momentum_distribution.pdf")
    #-----------------------------------------

    #---Photoelectron angular distributions---
    elseif  observable == StrongField.SfaAzimuthalAngularDistribution || observable == StrongField.SfaPolarAngularDistribution
        gr() #sets the plotting backend to the package "GR"

        for k = 1:size(comps,1) #Loop over all comps

            #prepare data
            if  observable == StrongField.SfaAzimuthalAngularDistribution
                angleList = results[k]["angular distribution"].phis
                plotTitle = "Photoelectron angular distribution (azimuthal)"
            elseif observable == StrongField.SfaPolarAngularDistribution
                angleList = results[k]["angular distribution"].thetas
                plotTitle = "Photoelectron angular distribution (polar)"
            end
            probList = results[k]["angular distribution"].probabilities

            lineLabel = ""
            if energyScale == "atomic"
                lineLabel = string( round( results[k]["angular distribution"].energy, digits = 2) ) * " a.u."
                legendTitle = "Energy"
            elseif energyScale == "omega"
                lineLabel = string( round( results[k]["angular distribution"].energy / comps[k].beam.omega, digits = 1 ) ) * " w"
                legendTitle = "Energy"
            end

            if  k == 1 #initial plot generation
                GR.polar(angleList, probList)
                Plots.plot(angleList,probList,
                            title = plotTitle,
                            proj = :polar,
                            line = (2,styles[1]),
                            linecolor = colors[1],
                            labelfontfamily = "Latin Modern Roman",
                            titlefontfamily = "Latin Modern Roman",
                            gridlinewidth = 2,
                            tickfontsize = 10,
                            labelfontsize = 10,
                            label = lineLabel,
                            legendtitle = legendTitle,
                            legend = true

                            )
            else
                Plots.plot!(angleList,probList,
                            proj = :polar,
                            line = (2,styles[k]),
                            linecolor = colors[k],
                            gridlinewidth = 2,
                            tickfontsize = 10,
                            labelfontsize = 10,
                            label = lineLabel,
                            labelfontfamily = "Latin Modern Roman",

                            )
            end

        end

        #export the plot as pdf-file
        savefig(dataLabel * "-angular_distribution.pdf")
    #-----------------------------------------

    #---Not a valid obserable---
    else
            error("Undefined observable for strong-field computations in StrongField.plot().")
    end
    #---------------------------

end

"""
`StrongField.exportData( comps::Array{StrongField.Computation,1}, results::Array{Dict{String,Any},1}, dataLabel::String = "StrongFieldData" )`
    ... exports the results = [Array1 Array2 ...] returned by StrongField.perform with the StrongField computations comps = [computation1 computation2 ...]
        into files with name dataLabel-dataType-1.csv, dataLabel-dataType-2.csv, ... where dataType = energy_distribution, azimuthal_angular_distribution, etc.
"""
function exportData( comps::Array{StrongField.Computation,1}, results::Array{Dict{String,Any},1}, dataLabel::String = "StrongFieldData" )
    for j = 1:size(results)[1]
        w = results[j]
        if typeof(comps[j].observable) == StrongField.SfaEnergyDistribution
            energyDistribution = w["energy distribution"]
            writedlm(dataLabel * "-energy_distribution" * "-" * string(j) * ".csv",hcat(energyDistribution.energies,energyDistribution.probabilities))
        elseif typeof(comps[j].observable) == StrongField.SfaMomentumDistribution
            angularDistribution = w["momentum distribution"]
            writedlm(dataLabel * "-momentum_distribution" * "-" * string(j) * ".csv",hcat(angularDistribution.phis,angularDistribution.probabilities))
        elseif typeof(comps[j].observable) == StrongField.SfaAzimuthalAngularDistribution
            angularDistribution = w["angular distribution"]
            writedlm(dataLabel * "-azimuthal_angular_distribution" * "-" * string(j) * ".csv",hcat(angularDistribution.phis,angularDistribution.probabilities))
        elseif typeof(comps[j].observable) == StrongField.SfaPolarAngularDistribution
            angularDistribution = w["angular distribution"]
            writedlm(dataLabel * "-polar_angular_distribution" * "-" * string(j) * ".csv",hcat(angularDistribution.phis,angularDistribution.probabilities))
        end
    end
end

#------------------------TEST: Plot the radial wave functions (JAC + hydrogenic)------------------------------------
"""
`StrongField.exportRadialWavefunctions( comps::Array{StrongField.Computation,1}, dataLabel::String = "StrongFieldData", savePlot::Bool = false )`
    ... exports the radial wave functions (initial state) = [Array1 Array2 ...] that are used in the StrongField computations comps = [computation1 computation2 ...]
        into files with name dataLabel-radial_wavefunction-1.csv, dataLabel-radial_wavefunction-2.csv, etc.
        - If savePlot == true, the wave functions are also plotted and exported into a single figure dataLabel-radial_wave_function.pdf
"""
function exportRadialWavefunctions( comps::Array{StrongField.Computation,1}, dataLabel::String = "StrongFieldData", savePlot::Bool = false )

        minIonizationPotential = Float64

        for k = 1:size(comps)[1]

            #Extract the initial orbital of the active electron from the many-electron comp.initialLevel and set quantum numbers
            initialOrbitals = comps[k].initialLevel.basis.orbitals

            #Find highest lying orbital (smallest ionization potential)
            defaultSubshell = [sh for (sh,or) in initialOrbitals][1] #This is not nice; must be a better way to simply get a default element from a Dict
            o = initialOrbitals[defaultSubshell]
            minIonizationPotential = abs(o.energy)
            for (subshell,orbital) in initialOrbitals
                if   abs(orbital.energy) < minIonizationPotential
                    o = orbital
                    minIonizationPotential = abs(orbital.energy)
                end
            end

            ls = LevelSymmetry(o.subshell)
            n = o.subshell.n;      l = Int((ls.J.num+1)/2);    j = ls.J.num/2;

            if  (sign((-1)^l) == -1 && ls.parity == plus::Parity) || (sign((-1)^l) == 1 && ls.parity == minus::Parity)
                l = l - 1
            end
            l = floor(Int,l)

            if  comps[k].settings.hydrogenic
                if comps[k].settings.hydrogenic1s
                    P = StrongField.HydrogenPnl( o.energy, 1, 0, rGrid.r )
                else
                    P = StrongField.HydrogenPnl( o.energy, n, l, rGrid.r )
                end
            else
                P = o.P
            end

            writedlm(dataLabel * "-initial_radial_wavefunction" * string(k) * ".csv",hcat(comps[k].grid.r,P))

            if savePlot
                if k == 1
                    Plots.plot(comps[k].grid.r, P,
                                                title = "Radial wave functions",
                                                xlabel = "r (a.u.)",
                                                ylabel = "P(r)",
                                                markershape = :circle,
                                                gridlinewidth = 2,
                                                tickfontsize = 10,
                                                labelfontsize = 10,
                                                labelfontfamily = "Latin Modern Roman",
                                                titlefontfamily = "Latin Modern Roman"
                                )
                else
                    Plots.plot!(comps[k].grid.r, P)
                end

            end

        end

        if savePlot
            savefig(dataLabel * "-initial_radial_wave_function.pdf")
        end
end

#if false
#    minIonizationPotential = 0.
    #Extract the initial orbital of the active electron from the many-electron comp.initialLevel and set quantum numbers
#    initialOrbitals = initialLevel.basis.orbitals

    #Find highest lying orbital (smallest ionization potential)
#    defaultSubshell = [sh for (sh,or) in initialOrbitals][1] #This is not nice; must be a better way to simply get a default element from a Dict
#    o = initialOrbitals[defaultSubshell]
#    minIonizationPotential = abs(o.energy)
#    for (subshell,orbital) in initialOrbitals
#        if   abs(orbital.energy) < minIonizationPotential
#            global o = orbital
#            global minIonizationPotential = abs(orbital.energy)
#        end
#    end

#    ls = LevelSymmetry(o.subshell)
#    n = o.subshell.n;      l = Int((ls.J.num+1)/2);    j = ls.J.num/2;
#
#    if  (sign((-1)^l) == -1 && ls.parity == plus::Parity) || (sign((-1)^l) == 1 && ls.parity == minus::Parity)
#        l = l - 1
#    end
#    l = floor(Int,l)
#
#    hydrogenP = StrongField.HydrogenPnl( o.energy, n, l, rGrid.r )
#
#    Plots.plot(rGrid.r, [o.P hydrogenP],
#                                title = "Radial wave function: n=" * string(n) * ", l=" * string(l) * ", Ip=" * string(round(convertUnits("energy: from atomic to eV", o.energy),digits=2)) * " eV",
#                                xlabel = "r (a.u.)",
#                                ylabel = "P(r)",
#                                markershape = :circle,
#                                gridlinewidth = 2,
#                                tickfontsize = 10,
#                                labelfontsize = 10,
#                                labelfontfamily = "Latin Modern Roman",
#                                titlefontfamily = "Latin Modern Roman",
#                                label = ["JAC" "Hydrogen"]
#                              )
#
#    savefig("radial_wavefunction.pdf")
#
#    writedlm(dataName * "-radial_wavefunction.csv",hcat(rGrid.r,o.P,hydrogenP))
#end
##-----------------------------------------------------------------------------------------------------------------------

#if false
#epsilonp=2*omega
#kappa =-1
#lp=0

#nrContinuum = Continuum.gridConsistency(epsilonp, rGrid)
#contSettings = Continuum.Settings(false, nrContinuum)
#contSettings = Continuum.Settings(false, rGrid.NoPoints)

#newiLevel = Basics.generateLevelWithSymmetryReducedBasis(initialLevel, initialLevel.basis.subshells)
#newfLevel = Basics.generateLevelWithSymmetryReducedBasis(finalLevel, newiLevel.basis.subshells)
#newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, kappa), newiLevel)
#cOrbital, phase  = Continuum.generateOrbitalForLevel(epsilonp, Subshell(101, kappa), newfLevel, nuclearModel, rGrid, contSettings)

#Plots.plot( rGrid.r[1:size(cOrbital.P)[1]], real(cOrbital.P * exp(im*phase)),
#                            title = "Radial wave function continuum",
#                            xlabel = "r (a.u.)",
#                            ylabel = "P(r)",
#                            markershape = :circle,
#                            gridlinewidth = 2,
#                            tickfontsize = 10,
#                            labelfontsize = 10,
#                            labelfontfamily = "Latin Modern Roman",
#                            titlefontfamily = "Latin Modern Roman"
#          )
#
#    Z=1.0
#    cVolkov = StrongField.CoulombVolkovP( epsilonp, lp, Z, rGrid.r )
#    Plots.plot!( rGrid.r, real(cVolkov) )
#
#    savefig("radial_wavefunction_continuum.pdf")

#end


text:
using PyPlot
using LaTeXStrings

data = readdlm("error_vs_length",'\t')

error = data[:,1]
l = data[:,2]

species_indices = Dict()
for i in 1:length(labels)
	species_indices[labels[i]] = i
end

error /= error[1]
#l /= l[1]

m = 8
a = .3

no_glc = Array{Int64}(readdlm("no_glc"))
no_g6p = Array{Int64}(readdlm("no_g6p"))

figure(figsize=(12,8))
for i in 1:length(error)
	if in(i,no_glc)
		semilogx(error[i],l[i],color="orange","o",markersize=m,markeredgecolor="none",alpha=a)
	else
		semilogx(error[i],l[i],"go",markersize=m,markeredgecolor="none",alpha=a)
	end
end

semilogx(error[1],l[1],"k*",markersize=m*2,markeredgecolor="none")
xlabel("error / base case error")
ylabel("length")
savefig("New.pdf")