text:
mutable struct Console{T<:GtkTextView, B<:GtkTextBuffer} <: GtkScrolledWindow

    handle::Ptr{Gtk.GObject}
    view::T
    buffer::B
    prompt_position::Integer
    prompt::String
    stdout_buffer::IOBuffer
    worker::TCPSocket
    worker_port::Int
    worker_idx::Int
    busy::Bool
    history::HistoryProvider
    main_window
    eval_in::Module
    mode::REPLMode

    function Console{T, B}(w_idx::Int, main_window, worker::TCPSocket,
        init_fun=init_view_buffer!, b_args=()) where {T<:GtkTextView, B<:GtkTextBuffer}

        b = B(b_args...)
        v = T(b)

        prompt = "Main>"
        set_gtk_property!(b, :text, prompt)

        init_fun(v, b)

        sc = GtkScrolledWindow()
        set_gtk_property!(sc, :hscrollbar_policy, 1)
        push!(sc, v)

        history = setup_history(w_idx)
        worker = isnothing(worker) ? w_idx : worker

        n = new(sc.handle, v, b, length(prompt)+1, prompt, IOBuffer(), worker, 0, w_idx,
            false, history, main_window, Main, NormalMode()
        )
        Gtk.gobject_move_ref(n, sc)
    end
end
Console(w_idx::Int, main_window, worker::TCPSocket) = Console{GtkTextView, GtkTextBuffer}(w_idx, main_window, worker)

function init_view_buffer!(v, b)
    set_gtk_property!(v, :margin_bottom, 10)
    set_gtk_property!(v, :margin_left, 4)
end

console_manager(c::Console) = console_manager(c.main_window)
worker(c::Console) = c.worker_idx == 1 ? c.worker_idx : c.worker

Base.pwd(c::Console) = remotecall_fetch(pwd, worker(c))

import Base.write
function write(c::Console, str::AbstractString)
    insert!(c.buffer, end_iter(c.buffer), str)
    place_cursor(c.buffer, end_iter(c.buffer))
end

function switch_mode(c::Console, mode::REPLMode)
    previous_prompt = prompt(c)
    c.mode = mode
    switch_prompt(c, mode, previous_prompt)
end

function check_switch_mode(c::Console, event)
    for mode in REPLModes
        #if doing(switch_key(mode), event)
        if (switch_key(mode) == event.keyval) && (mode != c.mode)
            switch_mode(c, mode)
            return true
        end
    end
    false
end

function switch_prompt(c::Console, mode::REPLMode, previous_prompt::String)

    its = GtkTextIter(c.buffer, c.prompt_position-length(previous_prompt))
    ite = GtkTextIter(c.buffer, c.prompt_position)
    replace_text(c.buffer, its, ite, prompt(c))

    diff = length(prompt(c))-length(previous_prompt)
    c.prompt_position += diff
end

prompt(c::Console) = prompt(c, c.mode)

# used to write completions before the prompt when we have more than one
# here if I use c.buffer instead it leads to crashes on v1.6
function write_before_prompt(c::Console, buffer, str::AbstractString)

    it = GtkTextIter(buffer, c.prompt_position-length(prompt(c)))
    it = Gtk.mutable(it)
    #insert!(buffer, it, str)

    ccall((:gtk_text_buffer_insert, Gtk.libgtk), Nothing,
        (Ptr{Gtk.GObject}, Ptr{GtkTextIter}, Ptr{UInt8}, Cint), buffer, it, Gtk.bytestring(str), sizeof(str))

    c.prompt_position += length(str)

    it = GtkTextIter(buffer, c.prompt_position-length(prompt(c)))
    if get_text_left_of_iter(it) != "\n"
        it = GtkTextIter(buffer, c.prompt_position-length(prompt(c)))
        insert!(buffer, it, "\n")
        c.prompt_position += 1
    end
end

function new_prompt(c::Console)
    insert!(c.buffer, end_iter(c.buffer), "\n$(prompt(c))")
    c.prompt_position = length(c.buffer)+1
    place_cursor(c.buffer, end_iter(c.buffer))
end

function clear(c::Console)
    set_gtk_property!(c.buffer, :text, "")
end

##
function check_worker(c::Console)
    c.worker_idx == 1 && return true
    s = worker(c)
    if s.status < Sockets.StatusConnecting
        @info "worker not connect, try to reconnect on port $(c.worker_port)"
        try
            c.worker = connect(c.worker_port)
        catch err
            @warn err
            return false
        end
    end
    return true
end

function on_return(c::Console, cmd::AbstractString)

    cmd = string(strip(cmd))#avoid getting a substring
    buffer = c.buffer

    write(c, "\n")
    if !check_worker(c)
        write(c, "Couldn't connect to worker\n")
        return
    end

    push!(c.history, cmd)
    seek_end(c.history)

    found = false
    if c.mode == NormalMode()
        (found, t) = check_console_commands(cmd, c)
    end

    if !found
        @async on_return(c, c.mode, cmd)#we don't need the result right away
    else
        @async begin
            time = @elapsed wait(t)
            write_output_to_console(c.worker_idx, t.result, time)
        end
    end
    c.busy = true
    push!(c.main_window.statusBar, "console", "Busy")

    #g_timeout_add(() -> write_output_to_console(c, found), 50)
    nothing
end

function kill_current_task(c::Console)
    try #otherwise this makes the callback fail in some versions
        #@async Base.throwto(c.run_task, InterruptException())
        interrupt_task(c)
    catch
    end
end

interrupt_task(c::Console) = remotecall_fetch(RemoteGtkREPL.interrupt_task, worker(c))

function RemoteGtkREPL.process_message(m::RemoteGtkREPL.EvalDone)
    # to avoid writing at random times in the console we call the function
    # from Gtk's main loop
    g_idle_add(() -> write_output_to_console(m.console_idx, m.data, m.time))
end

RemoteGtkREPL.process_message(m::RemoteGtkREPL.StdOutData) = print_to_console_remote(m.data, m.console_idx)

@guarded (PROPAGATE) function write_output_to_console(console_idx::Int, result, time)

    c = get_console(main_window.console_manager, console_idx)

    if typeof(result) <: Tuple #console commands can return just a string
        str, v = result
    else
        str, v = (result, nothing)
    end

    !isnothing(v) && display(v)

    finalOutput = isnothing(str) ? "" : str

    if str == InterruptException()
        finalOutput = string(str) * "\n"
    end

    write(c, finalOutput)
    new_prompt(c)
    on_command_done(c.main_window, c)
    c.busy = false

    t = @sprintf("%4.6f", time)
    push!(c.main_window.statusBar, "console", "Run time $(t)s")

    return PROPAGATE#this will remove the function from the list of event sources
end

"Get the text after the prompt"
function command(c::Console)
    its = GtkTextIter(c.buffer, c.prompt_position)
    ite = GtkTextIter(c.buffer, length(c.buffer)+1)
    cmd = (its:ite).text[String]
    return cmd
end
function command(c::Console, str::AbstractString, offset::Integer)

    its = GtkTextIter(c.buffer, c.prompt_position)
    ite = GtkTextIter(c.buffer, length(c.buffer)+1)
    replace_text(c.buffer, its, ite, str)
    if offset >= 0 && c.prompt_position+offset <= length(c.buffer)
        place_cursor(c.buffer, c.prompt_position+offset)
    end
end
command(c::Console, str::AbstractString) = command(c, str, -1)

function move_cursor_to_end(c::Console)
    place_cursor(c.buffer, end_iter(c.buffer))
end
function move_cursor_to_prompt(c::Console)
    place_cursor(c.buffer, c.prompt_position)
end

"return cursor position in the prompt text"
function cursor_position(c::Console)
    a = c.prompt_position
    b = cursor_position(c.buffer)
    b-a+1
end
cursor_position(b::GtkTextBuffer) = get_gtk_property(b, :cursor_position, Int)

function select_on_ctrl_shift(direction, c::Console)

    buffer = c.buffer
    (found, its, ite) = selection_bounds(buffer)

    if direction == :start
        ite, its = its, ite
    end

    its = found ? nonmutable(buffer, its) : get_text_iter_at_cursor(buffer)

    direction == :start && move_cursor_to_prompt(c)
    direction == :end && move_cursor_to_sentence_end(buffer)

    ite = get_text_iter_at_cursor(buffer)
    select_range(buffer, ite, its)#invert here so the cursor end up on the far right
end

##
ismodkey(event::Gtk.GdkEvent, mod::Integer) =
    any(x -> x == event.keyval, [
        GdkKeySyms.Control_L, GdkKeySyms.Control_R,
        GdkKeySyms.Meta_L, GdkKeySyms.Meta_R,
        GdkKeySyms.Hyper_L, GdkKeySyms.Hyper_R,
        GdkKeySyms.Shift_L, GdkKeySyms.Shift_R
    ]) ||
    any(x -> x == event.state & mod, [
        GdkModifierType.CONTROL, GdkKeySyms.Meta_L, GdkKeySyms.Meta_R,
        PrimaryModifier, SHIFT, GdkModifierType.GDK_MOD1_MASK,
        SecondaryModifer, PrimaryModifier+SHIFT, PrimaryModifier+GdkModifierType.META])

before_prompt(console, pos::Integer) = pos < console.prompt_position
before_prompt(console) = before_prompt(console, get_gtk_property(console.buffer, :cursor_position, Int)+1)

before_or_at_prompt(console, pos::Integer) = pos <= console.prompt_position
before_or_at_prompt(console) = before_or_at_prompt(console, get_gtk_property(console.buffer, :cursor_position, Int)+1)
at_prompt(console, pos::Integer) = pos == console.prompt_position

function iters_at_console_prompt(console)
    its = GtkTextIter(console.buffer, console.prompt_position)
    ite = nonmutable(console.buffer, end_iter(console.buffer) )
    (its, ite)
end

#FIXME disable drag and drop text above cursor
@guarded (PROPAGATE) function console_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    textview = convert(GtkTextView, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    console = user_data
    buffer = console.buffer

    cmd = command(console)
    pos = cursor_position(console)
    prefix = length(cmd) >= pos ? cmd[1:pos] : ""

    mod = get_default_mod_mask()

    #put back the cursor after the prompt
    if before_prompt(console)
        #check that we are not trying to copy or something of the sort
        if !ismodkey(event, mod)
            move_cursor_to_end(console)
        end
    end

    if pos == 0#the cursor is at the prompt
        check_switch_mode(console, event) && return INTERRUPT
    end

    (found, it_start, it_end) = selection_bounds(buffer)

    #prevent deleting text before prompt
    if event.keyval == GdkKeySyms.BackSpace ||
       event.keyval == GdkKeySyms.Delete ||
       event.keyval == GdkKeySyms.Clear ||
       doing(Actions["cut"], event)

        if found

            before_prompt(console, char_offset(it_start)) && return INTERRUPT
        else
            before_or_at_prompt(console) && return INTERRUPT
        end
    end
    if doing(Actions["move_to_line_start"], event) ||
        doing(Action(GdkKeySyms.Left, PrimaryModifier), event)
        move_cursor_to_prompt(console)
        return INTERRUPT
    end
    if doing(Actions["move_to_line_end"], event) ||
       doing(Action(GdkKeySyms.Right, PrimaryModifier), event)
        move_cursor_to_end(console)
        return INTERRUPT
    end
    if doing(Actions["clear_console"], event)
        clear(console)
        new_prompt(console)
        return INTERRUPT
    end
    if doing(Action(GdkKeySyms.Right, PrimaryModifier+GdkModifierType.SHIFT), event)
        select_on_ctrl_shift(:end, console)
        return INTERRUPT
    end
    if doing(Action(GdkKeySyms.Left, PrimaryModifier+GdkModifierType.SHIFT), event)
        select_on_ctrl_shift(:start, console)
        return INTERRUPT
    end

    if doing(Action(GdkKeySyms.Left, NoModifier), event)
        if found
            at_prompt(console, char_offset(it_start)) && return INTERRUPT
        else
            before_or_at_prompt(console) && return INTERRUPT
        end
        return PROPAGATE
    end
    if doing(Action(GdkKeySyms.Left, GdkModifierType.SHIFT), event)

        at_prompt(console, char_offset(it_start)) && return INTERRUPT

        return PROPAGATE
    end

    if event.keyval == GdkKeySyms.Up
        if found
            if !before_prompt(console, offset(it_start))
                select_range(buffer, GtkTextIter(buffer, console.prompt_position), nonmutable(buffer, it_end))
                return INTERRUPT
            end
            return PROPAGATE
        end
        !history_up(console.history, prefix, cmd) && return convert(Cint, true)
        command(console, history_get_current(console.history), length(prefix))
        return INTERRUPT
    end

    if event.keyval == GdkKeySyms.Down
        hasselection(buffer) && return PROPAGATE
        history_down(console.history, prefix, cmd)
        command(console, history_get_current(console.history), length(prefix))

        return INTERRUPT
    end
    if event.keyval == GdkKeySyms.Tab
        #convert cursor position into index
        autocomplete(console, cmd, pos)
        return INTERRUPT
    end
    if doing(Actions["select_all"], event)
        #select all
        before_prompt(console) && return PROPAGATE
        #select only prompt
        its, ite = iters_at_console_prompt(console)
        select_range(buffer, its, ite)
        return INTERRUPT
    end
    if doing(Actions["interrupt_run"], event)
        kill_current_task(console)
        return INTERRUPT
    end
    if doing(Actions["copy"], event)
        auto_select_prompt(found, console, buffer)
        signal_emit(textview, "copy-clipboard", Nothing)
        return INTERRUPT
    end
    if doing(Actions["paste"], event)
        signal_emit(textview, "paste-clipboard", Nothing)
        return INTERRUPT
    end
    if doing(Actions["cut"], event)
        auto_select_prompt(found, console, buffer)
        signal_emit(textview, "cut-clipboard", Nothing)
        return INTERRUPT
    end

    return PROPAGATE
end

"""
Auto select the prompt text when nothing is selected
and we are trying to copy or cut.
"""
function auto_select_prompt(found, console, buffer)
    if !found && !before_prompt(console)
        its, ite = iters_at_console_prompt(console)
        select_range(buffer, its, ite)
    end
end

function _callback_only_for_return(widgetptr::Ptr, eventptr::Ptr, user_data)
    event = convert(Gtk.GdkEvent, eventptr)
    console = user_data
    buffer = console.buffer

    cmd = command(console)

    if event.keyval == GdkKeySyms.Return
        if !console.busy
            on_return(console, cmd)
        end
        return INTERRUPT
    end
    return PROPAGATE
end

## MOUSE CLICKS

@guarded (INTERRUPT) function _console_button_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    textview = convert(GtkTextView, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    buffer = get_gtk_property(textview, :buffer, GtkTextBuffer)
    console = user_data
    main_window = console.main_window

    if event.event_type == Gtk.GdkEventType.DOUBLE_BUTTON_PRESS
        select_word_double_click(textview, buffer, Int(event.x), Int(event.y))
        return INTERRUPT
    end

    mod = get_default_mod_mask()
    if Int(event.button) == 1 && Int(event.state & mod) == Int(PrimaryModifier)
        open_method(textview) && return INTERRUPT
    end

    if rightclick(event)
        menu = buildmenu([
            MenuItem("Close Console", remove_console_cb),
            MenuItem("Add Console", add_console_cb),
            MenuItem("Clear Console", clear_console_cb),
            MenuItem("Toggle Wrap Mode", toggle_wrap_mode_cb)
            #GtkSeparatorMenuItem,
            #MenuItem("Toggle Wrap Mode", kill_current_task_cb),
            ],
            (console_manager(main_window), console, main_window)
        )
        popup(menu, event)
        return INTERRUPT
    end

    return PROPAGATE
end

global const console_mousepos = zeros(Int, 2)
global const console_mousepos_root = zeros(Int, 2)

#FIXME replace this by the same thing at the window level ?
#or put this as a field of the type.
@guarded (PROPAGATE) function console_motion_notify_event_cb(widget::Ptr,  eventptr::Ptr, user_data)
    event = convert(Gtk.GdkEvent, eventptr)

    console_mousepos[1] = round(Int, event.x)
    console_mousepos[2] = round(Int, event.y)
    console_mousepos_root[1] = round(Int, event.x_root)
    console_mousepos_root[2] = round(Int, event.y_root)
    return PROPAGATE
end

##

## auto-scroll the textview
function console_scroll_cb(widgetptr::Ptr, rectptr::Ptr, user_data)

    c = user_data
    adj = get_gtk_property(c, :vadjustment, GtkAdjustment)
    set_gtk_property!(adj, :value,
        get_gtk_property(adj, :upper, AbstractFloat) -
        get_gtk_property(adj, :page_size, AbstractFloat)
    )
    adj = get_gtk_property(c, :hadjustment, GtkAdjustment)
    set_gtk_property!(adj, :value, 0)

    nothing
end

## Auto-complete

function complete_additional_symbols(str, S)
    comp = Array{String}(undef, 0)
    for s in S
        startswith(s, str) && push!(comp, s)
    end
    comp
end

# TODO use the same code than in the editor
function autocomplete(c::Console, cmd::AbstractString, pos::Integer)

    isempty(cmd) && return
    pos > length(cmd) && return

    scmd = JuliaWordsUtils.CharArray(cmd)

    (ctx, m) = console_commands_context(cmd)

    if c.mode == ShellMode()
        ctx, m = :file, nothing
    end

    lastpart = pos < length(scmd) ? scmd[pos+1:end] : ""
    cmd = scmd[1:pos]

    if ctx == :normal
        isempty(cmd) && return
        comp, dotpos = completions_in_module(cmd, c)
    end
    if ctx == :file
        (comp, dotpos) = remotecall_fetch(REPL.shell_completions, worker(c), cmd, lastindex(cmd))
        comp = REPL.completion_text.(comp)
    end

    update_completions(c, comp, dotpos, cmd, lastpart)
end

function completions_in_module(cmd, c::Console)
    prefix = string(c.eval_in, ".")
    @info prefix * cmd
    comp, dotpos = remotecall_fetch(RemoteGtkREPL.remote_completions, worker(c), prefix * cmd)
    @info comp
    dotpos = dotpos .- lastindex(prefix)
    comp, dotpos
end

##

# cmd is the word, including dots we are trying to complete
function update_completions(c::Console, comp, dotpos, cmd, lastpart)

    isempty(comp) && return

    dotpos = dotpos.start
    prefix = dotpos > 1 ? cmd[1:dotpos-1] : ""

    if(length(comp)>1)

        maxLength = maximum(map(length, comp))
        w = width(c.view)
        nchar_to_width(x) = 0.9*x*fontsize #TODO pango_font_metrics_get_approximate_char_width
        n_per_line = max(1, round(Int, w/nchar_to_width(maxLength)))

        out = "\n"
        for i = 1:length(comp)
            spacing = repeat(" ", maxLength-length(comp[i]))
            out = "$out $(comp[i]) $spacing"
            if mod(i, n_per_line) == 0
                out = out * "\n"
            end
        end

        write_before_prompt(c, c.buffer, out)
        out = prefix * REPL.LineEdit.common_prefix(comp)
    else
        out = prefix * comp[1]
    end

    offset = length(out)#place the cursor after the newly inserted piece
    #update entry
    out = out * lastpart
    #out = remove_filename_from_methods_def(out)

    command(c, out, offset)
end

function init!(c::Console)
    signal_connect(console_key_press_cb, c.view, "key-press-event",
    Cint, (Ptr{Gtk.GdkEvent}, ), false, c)
    signal_connect(_callback_only_for_return, c.view, "key-press-event",
    Cint, (Ptr{Gtk.GdkEvent}, ), false, c)
    signal_connect(_console_button_press_cb, c.view, "button-press-event",
    Cint, (Ptr{Gtk.GdkEvent}, ), false, c)
    signal_connect(console_motion_notify_event_cb, c, "motion-notify-event",
    Cint, (Ptr{Gtk.GdkEvent}, ), false)
    signal_connect(console_scroll_cb, c.view, "size-allocate", Nothing,
    (Ptr{Gtk.GdkRectangle}, ), false, c)

    # Note that due to historical reasons, GtkNotebook refuses to switch to a page unless the child widget is visible.
    # Therefore, it is recommended to show child widgets before adding them to a notebook.
    show(c)
    push!(console_manager(c), c)
    Gtk.GAccessor.tab_label_text(console_manager(c), c, "C" * string(length(console_manager(c))))
    g_timeout_add(()->print_to_console(c), 100)
end

"Run from the main Gtk loop, and print to console
the content of stdout_buffer"
# TODO : still needed ?
function print_to_console(console)

    s = String(take!(console.stdout_buffer))
    if !isempty(s)
        s = translate_colors(s)
        write(console, s)
    end
    if is_running
        return INTERRUPT
    else
        return Cint(false)
    end
end
#cfunction(print_to_console, Cint, Ptr{Console})

#FIXME dirty hack?
function translate_colors(s::AbstractString)

    s = replace(s, "\e[1m\e[31m" => "* ")
    s = replace(s, "\e[1m\e[32m" => "* ")
    s = replace(s, "\e[1m\e[31" => "* ")
    s = replace(s, "\e[0m" => "")
    s
end

@guarded (nothing) function toggle_wrap_mode_cb(btn::Ptr, user_data)
    ntbook, tab, main_window = user_data
    tab.view.wrap_mode[Int] = !Bool(tab.view.wrap_mode[Int])
    nothing
end
@guarded (nothing) function clear_console_cb(btn::Ptr, user_data)
    ntbook, tab, main_window = user_data
    clear(tab)
    new_prompt(tab)
    nothing
end
@guarded (nothing) function kill_current_task_cb(btn::Ptr, user_data)
    ntbook, tab, main_window = user_data
    kill_current_task(tab)
    nothing
end
@guarded (nothing) function add_console_cb(btn::Ptr, user_data)
    console_manager, tab, main_window = user_data
    #add_console(main_window)
    add_remote_console(main_window, console_manager.top_module)
    nothing
end
@guarded (nothing) function remove_console_cb(btn::Ptr, user_data)
    ntbook, tab, main_window = user_data
    idx = index(ntbook, tab)
    if idx != 1#can't close the main console
        c = ntbook[idx]
#        remotecall_fetch(info, worker(c), "Goodbye.")
        splice!(ntbook, idx)
        index(ntbook, max(idx-1, 0))
    end
    nothing
end

get_current_console(console_mng::GtkNotebook) = console_mng[index(console_mng)]

function get_console(cm::ConsoleManager, idx)
    for i = 1:length(cm)
        c = get_tab(cm, i)
        c.worker_idx == idx && return c
    end
end

#this is called by remote workers
function print_to_console_remote(s, idx::Integer)
    #copy the output to the right console buffer
    c = get_console(main_window.console_manager, idx)
    if !isnothing(c)
        write(c.stdout_buffer, s)
    else
        @warn "Failed to print $s for Console $(idx)."
    end
end

## REDIRECT_STDOUT for main console

function send_stream(rd::IO, stdout_buffer::IO)
    nb = bytesavailable(rd)
    if nb > 0
        d = read(rd, nb)
        s = String(copy(d))

        if !isempty(s)
            write(stdout_buffer, s)
        end
    end
end

function watch_stream(rd::IO, c::Console)
    while !eof(rd) && is_running # blocks until something is available
        send_stream(rd, c.stdout_buffer)
        sleep(0.01) # a little delay to accumulate output
    end
end

#


text:
"""
    sphere(n)

A minimal triangulation of the `n`-dimensional sphere.
"""
function sphere(n::Integer)
    s = Simplex(collect(1:n+1))
    return SimplicialComplex(s)
end

"""
    torus()

A minimal triangulation of the torus.

This is a simplicial complex with 7 vertices, 21 edges and
14 faces. It is the unique triangulation of the torus with
7 vertices, and has been found by Möbius in 1861.
"""
function torus()
    cplx = SimplicialComplex{Simplex{Int}}()
    for s in [[1, 2, 3], [2, 3, 5], [2, 4, 5], [2, 4, 7], [1, 2, 6],
              [2, 6, 7], [3, 4, 6], [3, 5, 6], [3, 4, 7], [1, 3, 7],
              [1, 4, 5], [1, 4, 6], [5, 6, 7], [1, 5, 7]]
        addsimplices!(cplx, Simplex(s))
    end
    return cplx
end

"""
    kleinbottle()

A minimal triangulation of the Klein bottle, as presented for example
in Davide Cervone's thesis.
"""
function kleinbottle()
    cplx = SimplicialComplex{Simplex{Int}}()
    for s in [[3, 4, 8], [2, 3, 4], [2, 4, 6], [2, 6, 8], [2, 5, 8],
              [3, 5, 7], [2, 3, 7], [2, 7, 1], [2, 5, 1], [3, 5, 1],
              [4, 5, 8], [4, 5, 7], [4, 6, 7], [6, 7, 1], [3, 6, 1], [3, 6, 8]]
        addsimplices!(cplx, Simplex(s))
    end
    return cplx
end

"""
    projectiveplane()

A minimal triangulation of the real projective plane.
"""
function projectiveplane()
    cplx = SimplicialComplex{Simplex{Int}}()
    for s in [[1, 2, 3], [1, 3, 4], [1, 2, 6], [1, 5, 6], [1, 4, 5],
              [2, 3, 5], [2, 4, 5], [2, 4, 6], [3, 4, 6], [3, 5, 6]]
        addsimplices!(cplx, Simplex(s))
    end
    return cplx
end

"""
    randomcomplex(n, d)

A random `d`-dimensional simplicial complex on `n` vertices.
"""
function randomcomplex(n, d)
    cplx = SimplicialComplex(map(Simplex, 1:n))
    if d+1 >= n
        addsimplices!(cplx, Simplex(collect(1:n)))
    else
        pts = collect(1:n)
        simsset = Set{Set{Int}}()
        while length(simsset) < n
            shuffle!(pts)
            push!(simsset, Set(pts[1:d+1]))
        end
        for s in simsset
            addsimplices!(cplx, Simplex(s...))
        end
    end
    return cplx
end

ext: jl


text:
using Oceananigans.Advection
using Oceananigans.Buoyancy
using Oceananigans.Coriolis
using Oceananigans.Operators
using Oceananigans.SurfaceWaves
using Oceananigans.TurbulenceClosures: ∂ⱼ_2ν_Σ₁ⱼ, ∂ⱼ_2ν_Σ₂ⱼ, ∂ⱼ_2ν_Σ₃ⱼ, ∇_κ_∇c

"""
    u_velocity_tendency(i, j, k, grid,
                        advection,
                        coriolis,
                        surface_waves,
                        closure,
                        background_fields,
                        velocities,
                        tracers,
                        diffusivities,
                        forcings,
                        hydrostatic_pressure,
                        clock)

Return the tendency for the horizontal velocity in the x-direction, or the east-west
direction, ``u``, at grid point `i, j, k`.

The tendency for ``u`` is called ``G_u`` and defined via

    ``∂_t u = G_u - ∂_x ϕ_n``

where ∂_x ϕ_n is the non-hydrostatic pressure gradient in the x-direction.

`coriolis`, `surface_waves`, and `closure` are types encoding information about Coriolis
forces, surface waves, and the prescribed turbulence closure.

`background_fields` is a `NamedTuple` containing background velocity and tracer
`FunctionFields`.

The arguments `velocities`, `tracers`, and `diffusivities` are `NamedTuple`s with the three
velocity components, tracer fields, and precalculated diffusivities where applicable.
`forcings` is a named tuple of forcing functions. `hydrostatic_pressure` is the hydrostatic
pressure anomaly.

`clock` keeps track of `clock.time` and `clock.iteration`.
"""
@inline function u_velocity_tendency(i, j, k, grid,
                                     advection,
                                     coriolis,
                                     surface_waves,
                                     closure,
                                     background_fields,
                                     velocities,
                                     tracers,
                                     diffusivities,
                                     forcings,
                                     hydrostatic_pressure,
                                     clock)

    return ( - div_Uu(i, j, k, grid, advection, velocities, velocities.u)
             - div_Uu(i, j, k, grid, advection, background_fields.velocities, velocities.u)
             - div_Uu(i, j, k, grid, advection, velocities, background_fields.velocities.u)
             - x_f_cross_U(i, j, k, grid, coriolis, velocities)
             - ∂xᶠᵃᵃ(i, j, k, grid, hydrostatic_pressure)
             + ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, clock, closure, velocities, diffusivities)
             + x_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, velocities, clock.time)
             + ∂t_uˢ(i, j, k, grid, surface_waves, clock.time)
             + forcings.u(i, j, k, grid, clock, merge(velocities, tracers)))
end

"""
    v_velocity_tendency(i, j, k, grid,
                        advection,
                        coriolis,
                        surface_waves,
                        closure,
                        background_fields,
                        velocities,
                        tracers,
                        diffusivities,
                        forcings,
                        hydrostatic_pressure,
                        clock)

Return the tendency for the horizontal velocity in the y-direction, or the north-south
direction, ``v``, at grid point `i, j, k`.

The tendency for ``v`` is called ``G_v`` and defined via

    ``∂_t v = G_v - ∂_y ϕ_n``

where ∂_y ϕ_n is the non-hydrostatic pressure gradient in the y-direction.

`coriolis`, `surface_waves`, and `closure` are types encoding information about Coriolis
forces, surface waves, and the prescribed turbulence closure.

`background_fields` is a `NamedTuple` containing background velocity and tracer
`FunctionFields`.

The arguments `velocities`, `tracers`, and `diffusivities` are `NamedTuple`s with the three
velocity components, tracer fields, and precalculated diffusivities where applicable.
`forcings` is a named tuple of forcing functions. `hydrostatic_pressure` is the hydrostatic
pressure anomaly.

`clock` keeps track of `clock.time` and `clock.iteration`.
"""
@inline function v_velocity_tendency(i, j, k, grid,
                                     advection,
                                     coriolis,
                                     surface_waves,
                                     closure,
                                     background_fields,
                                     velocities,
                                     tracers,
                                     diffusivities,
                                     forcings,
                                     hydrostatic_pressure,
                                     clock)

    return ( - div_Uv(i, j, k, grid, advection, velocities, velocities.v)
             - div_Uv(i, j, k, grid, advection, background_fields.velocities, velocities.v)
             - div_Uv(i, j, k, grid, advection, velocities, background_fields.velocities.v)
             - y_f_cross_U(i, j, k, grid, coriolis, velocities)
             - ∂yᵃᶠᵃ(i, j, k, grid, hydrostatic_pressure)
             + ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, clock, closure, velocities, diffusivities)
             + y_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, velocities, clock.time)
             + ∂t_vˢ(i, j, k, grid, surface_waves, clock.time)
             + forcings.v(i, j, k, grid, clock, merge(velocities, tracers)))
end

"""
    w_velocity_tendency(i, j, k, grid,
                        advection,
                        coriolis,
                        surface_waves,
                        closure,
                        background_fields,
                        velocities,
                        tracers,
                        diffusivities,
                        forcings,
                        clock)

Return the tendency for the vertical velocity ``w`` at grid point `i, j, k`.
The tendency for ``w`` is called ``G_w`` and defined via

    ``∂_t w = G_w - ∂_z ϕ_n``

where ∂_z ϕ_n is the non-hydrostatic pressure gradient in the z-direction.

`coriolis`, `surface_waves`, and `closure` are types encoding information about Coriolis
forces, surface waves, and the prescribed turbulence closure.

`background_fields` is a `NamedTuple` containing background velocity and tracer
`FunctionFields`.

The arguments `velocities`, `tracers`, and `diffusivities` are `NamedTuple`s with the three
velocity components, tracer fields, and precalculated diffusivities where applicable.
`forcings` is a named tuple of forcing functions.

`clock` keeps track of `clock.time` and `clock.iteration`.
"""
@inline function w_velocity_tendency(i, j, k, grid,
                                     advection,
                                     coriolis,
                                     surface_waves,
                                     closure,
                                     background_fields,
                                     velocities,
                                     tracers,
                                     diffusivities,
                                     forcings,
                                     clock)

    return ( - div_Uw(i, j, k, grid, advection, velocities, velocities.w)
             - div_Uw(i, j, k, grid, advection, background_fields.velocities, velocities.w)
             - div_Uw(i, j, k, grid, advection, velocities, background_fields.velocities.w)
             - z_f_cross_U(i, j, k, grid, coriolis, velocities)
             + ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, clock, closure, velocities, diffusivities)
             + z_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, velocities, clock.time)
             + ∂t_wˢ(i, j, k, grid, surface_waves, clock.time)
             + forcings.w(i, j, k, grid, clock, merge(velocities, tracers)))
end

"""
    tracer_tendency(i, j, k, grid,
                    val_tracer_index::Val{tracer_index},
                    advection,
                    closure,
                    buoyancy,
                    background_fields,
                    velocities,
                    tracers,
                    diffusivities,
                    forcing,
                    clock)

Return the tendency for a tracer field with index `tracer_index`
at grid point `i, j, k`.

The tendency is called ``G_c`` and defined via

    ``∂_t c = G_c``

where `c = C[tracer_index]`.

`closure` and `buoyancy` are types encoding information about the prescribed
turbulence closure and buoyancy model.

`background_fields` is a `NamedTuple` containing background velocity and tracer
`FunctionFields`.

The arguments `velocities`, `tracers`, and `diffusivities` are `NamedTuple`s with the three
velocity components, tracer fields, and precalculated diffusivities where applicable.
`forcings` is a named tuple of forcing functions.

`clock` keeps track of `clock.time` and `clock.iteration`.
"""
@inline function tracer_tendency(i, j, k, grid,
                                 val_tracer_index::Val{tracer_index},
                                 advection,
                                 closure,
                                 buoyancy,
                                 background_fields,
                                 velocities,
                                 tracers,
                                 diffusivities,
                                 forcing,
                                 clock) where tracer_index

    @inbounds c = tracers[tracer_index]
    @inbounds background_fields_c = background_fields.tracers[tracer_index]

    return ( - div_Uc(i, j, k, grid, advection, velocities, c)
             - div_Uc(i, j, k, grid, advection, background_fields.velocities, c)
             - div_Uc(i, j, k, grid, advection, velocities, background_fields_c)
             + ∇_κ_∇c(i, j, k, grid, clock, closure, c, val_tracer_index, diffusivities, tracers, buoyancy)
             + forcing(i, j, k, grid, clock, merge(velocities, tracers)))
end


text:
using InteractionNetworkModels, Plots, StatsPlots, LinearAlgebra

S = [[1,2],[2,1],[1,2],[2,3],[3,4],[1,2],[2,1]]

K_in_ub = 9
K_out_ub = 19
K_inner = DimensionRange(2,K_in_ub)
K_outer = DimensionRange(1,K_out_ub)
γ_s = 3.8
γ_m = 3.0
d_s = FastEditDistance(FastLCS(20), 20)
d_m = MatchingDist(FastLCS(20))

V = 1:9
sis_model = SIS(S, γ_s, d_s, V, K_inner, K_outer)
sim_model = SIM(S, γ_m, d_m, V, K_inner, K_outer)

sis_mcmc = SisMcmcInsertDelete(
    ν_ed=3, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.8,1,K_in_ub),
    burn_in=2000, lag=50, desired_samples=500
)

sim_mcmc = SimMcmcInsertDelete(
    ν_ed=3, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.8,1,K_in_ub),
    burn_in=2000, lag=100, desired_samples=500
)

sis_out = sis_mcmc(sis_model)
sim_out = sim_mcmc(sim_model)

plot(sis_out)
plot!(sim_out)

summaryplot(sis_out)
summaryplot!(sim_out)

# Order info

order_inf_sis = [d_s(x,y)==0.0 ? 1.0 : d_m(x,y)/d_s(x,y) for x in sis_out.sample, y in sis_out.sample]
order_inf_sis = order_inf_sis[tril!(trues(size(order_inf_sis)),-1)]

order_inf_sim = [d_s(x,y)==0.0 ? 1.0 : d_m(x,y)/d_s(x,y) for x in sim_out.sample, y in sim_out.sample]
order_inf_sim = order_inf_sim[tril!(trues(size(order_inf_sim)),-1)]

density(order_inf_sis)
density!(order_inf_sim)


text:
using Test

include("evaluation/test_diagnostics.jl")
include("evaluation/test_roc_auc.jl")


text:
"""
Metadata for an [`Invite`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/invite#invite-metadata-object).
"""
struct InviteMetadata
    inviter::User
    uses::Int
    max_uses::Int
    max_age::Int
    temporary::Bool
    created_at::DateTime
    revoked::Bool
end
@boilerplate InviteMetadata :constructors :docs :lower :merge :mock


text:
import DifferentialEquations
const DE = DifferentialEquations

# This is curiously slow
function sde_exercise()
    f = (u,p,t) -> 1.01*u
    g = (u,p,t) -> 0.87*u
    u0 = 0.5
    tspan = (0.0,1.0)
    prob = DE.SDEProblem(f,g,u0,tspan)
    sol = DE.solve(prob,reltol=1e-3,abstol=1e-3)
    return nothing
end

function ode_exercise()
    f = (u,p,t) -> -u
    u0 = 0.5
    tspan = (0., 1.)
    prob = DE.ODEProblem(f, u0, tspan)
    sol = DE.solve(prob)
    return nothing
end

function ode_exercise2()
    f = function(u,p,t)
        x, y, z = u
        sigma, rho, beta = p
        return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
    end
    u0 = [1.0,0.0,0.0]
    tspan = (0., 100.)
    p = [10.0,28.0,8/3]
    prob = DE.ODEProblem(f, u0, tspan, p)
    sol = DE.solve(prob,saveat=0.01)
    return nothing
end

# From ODE docs
function ode_exercise3()
    f(u,p,t) = 1.01*u
    u0 = 1/2
    tspan = (0.0,1.0)
    prob = ODEProblem(f,u0,tspan)
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
    return nothing
end

ode_exercise()
ode_exercise2()
ode_exercise3()
sde_exercise()


text:
const minstepsize = 1.5*(pi/180.0) #1 degree minimum stepsize for largest L

"""Generates a random start structure for the rotational search"""
function randomStartStructure(volume::SphericalHarmonicsVolume, K::Int64, L::Int64)
    new = deepcopy(volume)
    #Initialize with random unitary matrices
    for l = 2:2:L
        rot = complex(random_rotation(2*l+1))
        ctr = Umat(l) #complex to real transformation
        rtc = ctranspose(ctr) #real to complex transformation
        rot = rtc*rot*ctr
        for k = 1:K
            cvec_set(new,k,l, rot*cvec_get(volume,k,l))
        end
    end
    return new
end

""""Randomly pertubates a given structure by applying a single random rotations to different l's."""
function randomPertubation(volume::SphericalHarmonicsVolume, K::Int64, L::Int64, stepsize=0.1)
    new = deepcopy(volume)
    #Initialize with random unitary matrices
    for l = 2:2:L
        rot = complex(random_rotation_step(2*l+1, stepsize))
        ctr = Umat(l) #complex to real transformation
        rtc = ctranspose(ctr) #real to complex transformation
        rot = rtc*rot*ctr
        for k = 1:K
            cvec_set(new, k,l, rot*cvec_get(volume, k,l))
        end
    end
    return new
end

#--------------------------------------------------------------------------------
# Iterative projections approach
#--------------------------------------------------------------------------------
# function project_C(coeff)
#     newCoeff = deepcopy(coeff)
#     L = Int64(floor((KMAX-1)/2))
#     for l = 2:2:L
#         u,s,v = svd((c[l]*transpose(b[l])) * diagm([1:KMAX;])^2 * coefficients_to_matix(coeff,l)')
#         newCoeffMatrix = b[l]*c[l]*u*v'
#
#         for k = 1:KMAX
#           start = KMAX - (2*l)
#             cvec_set(newCoeff[k], l, newCoeffMatrix[k,start:KMAX])
#         end
#
#     end
#     return real_to_comp(newCoeff)
# end

# function three_photon_walk(coeff, klow = 2, K=10, stepsize=pi/180.0*33.0)
#     coeff = project_C(coeff)
#     bestprob = energy(coeff,ablists, basisindices,basis,K,L)
#     step = deepcopy(coeff)
#     for i = 1:4
# #         surf = map(backtransform, coeff)
# #         surf = map(abs, surf)
# #         coeff = map(transform, surf)
#         for l = 2:2:L
#             rot = complex(random_rotation_step(2*l+1, gaussian(l*stepsize, l*stepsize/2)))
#             rot =  rtc[l]*rot*ctr[l]#rtc*rot*ctr
#             for k = 1:KMAX cvec_set(step[k],l, rot*cvec_get(coeff[k],l)) end
#         end
#         prob = energy(step, K, L)
#         if prob > bestprob
#             bestprob = prob
#             coeff = deepcopy(step)
#         end
#     end
#     # return deleteTerms(coeff,K,L)
#     return coeff
# end

#Projects real electron density shells onto the contraints given by the 2 and 3 photon correlation
# function B2(density,klow=2, K=10)
#     Alm = sphBesselTransform(map(transform,density), true)
#     Almsurf = map(backtransform, Alm)
#
#     #Get amplitudes from projection
#     Ilm = getIntensityCoefficientsIndirect(Alm)
#     pIlm = three_photon_walk(Ilm, klow, K)
#     # pIlm = project_C(Ilm)
#     # pIlm = deleteTerms(pIlm, KMAX, L)
#     amplitudes = map((x) -> sqrt(max(real(backtransform(x)),0)),pIlm)
#
#     #Apply amplitudes
#     for k = 1:KMAX
#         for i = 1:length(Almsurf[k])
#             Almsurf[k][i] = Almsurf[k][i] / abs(Almsurf[k][i]) * amplitudes[k][i]
#         end
#     end
#
#     #go back to real space
#     newrlm = sphBesselTransform(map(transform, Almsurf), false)
#     return map(backtransform, newrlm)
# end

# function projection_refinement(;ablistsfilename="expdata/ablists32.dat", L = 8,klow=2, K = 10, number = 166)
#
#     #loading triplet data
#     loadHistograms(; klow=klow, K=K )
#
#     #precalc rotations
#     global ctr = Dict([l => Umat(l) for l=2:2:L])
#     global rtc = Dict([l => ctranspose(ctr[l]) for l=2:2:L])
#
#     #start with random spherical harmonics coefficients representing a function
#     rlm = getDensityCoefficients( (x) -> x <= rmax/2.0 ? ((p,t) -> 1 + rand()) : ((p,t) -> 0.0) )
#     d = map(backtransform,rlm)
#
#     iterations = 1000
#     beta0 = 0.75
#     beta_max = 0.75
#     tau= 350
#     global Koff = KMAX / 2
#
#   	tmp1 = 2* B2(d, klow, K) - d #--reflection
#     tmp2 = 0.0
#     # diffa = deepcopy(tmp1)
#
#     progress = Progress(iterations, 1, "Structure determination...", 50)
#   	for n = 1:iterations
#       	beta = exp((-n/tau)^3.0)*beta0 + (1-exp((-n/tau)^3.0))*beta_max
#     		tmp3 = A(tmp1)
#     		tmp_u = 0.5*(beta*(2.0*tmp3 -tmp1) + (1-beta)*tmp1 + d)
#     		tmp2 = B2(tmp_u, klow, K)
#
#       	# for k = 1:length(d)
#       	# 	diffa[k] = d[k] - tmp_u[k]
#       	# end
#       	# diff = sum(abs, map(sum(abs, ,diffa))/sum(abs, map(sum(abs, ,d))
#         diff = 0.0
#
#     		d = deepcopy(tmp_u)
#     		tmp1 = 2*tmp2 - tmp_u
#
#       	progress.desc = "$diff: "
#       	next!(progress)
#   	end
#   	d = deepcopy(tmp2)
#   	print("Phasing complete...")
#     coeff = map(transform, d)
#     saveCube(getCube(coeff, cubesize), cubesize, "parallel/result$(number).mrc", qmax)
# end

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

"""Main rotation search method
`params` - rotation search parameters
`state` - starting state"""
function rotation_search(params = Dict("reference_pdb_path"=>"crambin.pdb","stepsizefactor"=>1.02, "initial_stepsize" => pi/180.0 * 180.0, "L"=>8, "K3_range" => 1:8, "N"=>32, "histograms"=>"expdata/correlations_N32_K25.dat", "optimizer"=>rotate_all_at_once, "initial_temperature_factor"=>1.0, "measure"=>"Bayes", "temperature_decay"=>0.99, "LMAX"=>25, "K2_range"=>1:35, "qmax"=>1.0, "lambda"=>0.0, "include_negativity"=>false), state = Dict{Any,Any}("newRun"=>true) )

    #Don't start a finished run
    if haskey(state, "state") && state["state"] == "finished_structure"
        return
    end

    #load the histogrammed dataset into the global scope
    c2ref_full, c2ref, c3ref_full, c3ref = loadHistograms(maximum(params["K3_range"]), maximum(params["K3_range"]), params["histograms"], false)

    #Open output file for logging
    out = open("det.out", state["newRun"] ? "w+" : "a")
    println("Starting on machine $(gethostname())")
    flush(STDOUT)

    if state["newRun"]

        #Mark it as "non"-new run for future continuations
        state["newRun"] = false

        #Lets start with the whole range by default (not applied for hierarchical approach)
        state["K"] = maximum(params["K3_range"])

        #Save reference intensity for later use
        if haskey(params, "reference_pdb_path") && params["reference_pdb_path"] != ""
            density,fourier,intensity = createSphericalHarmonicsStructure(params["reference_pdb_path"], params["LMAX"], maximum(params["K2_range"]), qmax(maximum(params["K2_range"]), params["qmax"]))
            state["reference_intensity"] = intensity
        end

        #Retrieve initial structure from 2p-correlation if not provided
        if !haskey(state,"intensity")
            state["intensity"] = retrieveSolution(c2ref_full/sum(abs, c2ref_full), params["L"], params["LMAX"], params["K2_range"], params["qmax"], params["lambda"])

            #added this to check determination run, becaus retrieveSolution was unstable
            # c2_theo = twoPhotons(state["reference_intensity"], BasisType(params["N"], params["L"], params["LMAX"], maximum(params["K3_range"]), params["lambda"], dq(state["reference_intensity"])), maximum(params["K2_range"]), true, false)
            #
            # state["intensity"] = retrieveSolution(c2_theo, params["L"], params["LMAX"], params["K2_range"],  params["qmax"], params["lambda"])

            state["intensity"] = randomStartStructure(state["intensity"], state["intensity"].KMAX, state["intensity"].LMAX)
        end

        state["stepsizes"] = Dict()
        state["L"] = 2
        state["i"] = 0
        state["state"] = "running"

        #Make some logging things
        write(out,"""
        #Starting structure determination with the following parameters:
        #$params
        """)

        saveState(params, state)
    end

    try
        CUDA_init()
    catch x
        println(x)
        println("!!! Init of CUDA failed")
        println("!!! Fallback to CPU-backend.")
    end
    flush(STDOUT)

    #begin optimization with given optimizer
    params["optimizer"](out, params, state, c3ref)

    flush(out)
    close(out)

    #Save the final state
    state["state"] = "finished_structure"
    saveState(params, state)
    return params,state
end

"""Creates a Monte Carlo step with given stepsize"""
function get_MC_step(intensity::SphericalHarmonicsVolume, basis::AbstractBasisType, stepsizes::Dict)
    step = deepcopy(intensity)
    for l = 2:2:basis.L
        rot =  basis.rtc[l]*complex(random_rotation_step(2*l+1, gaussian(0.0, stepsizes[l])))*basis.ctr[l]
        for k = 1:intensity.KMAX
            cvec_set(step,k,l, rot*cvec_get(intensity,k,l))
        end
    end
    return step
end

"""Given a particular configuration, calculates the average energy of changes to the structure"""
function calculate_temperature(state::Dict, params::Dict, basis::AbstractBasisType, c3ref::C3, iterations::Int64=50)
    if params["initial_temperature_factor"] > 0.0
        println("Explore energy landscape for temperature... ")
        #Measure average energy variation for 100 steps with the initial_stepsize
        energy_sample = Float64[]

        for i = 1:iterations
            step = get_MC_step(state["intensity"], basis, state["stepsizes"])
            e = energy(step, basis, c3ref, params["K3_range"], params["measure"])
            println("Energy calculation step $i: $e")
            push!(energy_sample,e)
        end
        return std(energy_sample)*params["initial_temperature_factor"]
    else
        return 0.0
    end

end

"""Saves current state and paramters to a file."""
function saveState(params, state)
    serializeToFile("params.dat", params)
    serializeToFile("state.dat", state)
end

"""Regularly do something"""
function regular_action(out, state, params)
    #save out information every 100 iterations
    if mod(state["i"], 500) == 0
        @everywhere gc()
        saveState(params, state)
        flush(out)
        flush(STDOUT)
    end
end

"""Optimizes structure by changing all rotation matrices at once and comparing the structure's three photon correlation with the experimental one"""
function rotate_all_at_once(out, params::Dict, state::Dict, c3ref::C3)

    state["L"] = params["L"]
    d_basis = complexBasis_choice(state["L"], params["LMAX"], params["N"], maximum(params["K3_range"]), params["lambda"], dq(state["intensity"]))

    #Check if this is a continuation or a new run by checking stepsizes
    if !haskey(state["stepsizes"], state["L"])

        #gradually increase K
        state["K"] = maximum(params["K3_range"])

        #set stepsize for this resolution
        #This makes sure that the initial stepsize for L=2 stays as intended and the other L stepsizes are hierarchical higher
        # state["stepsizes"] = Dict(l=>params["initial_stepsize"]*2^(l/2.0-1.0) for l=2:2:state["L"])
        state["stepsizes"] = Dict(l=>params["initial_stepsize"]*(l-1) for l=2:2:state["L"])

        #calculate the initial energy with that resolution
        state["E"] = energy(state["intensity"], d_basis, c3ref, params["K3_range"], params["measure"])

        #Calculate reference energy
        state["reference_energy"] = haskey(state, "reference_intensity") ? energy(state["reference_intensity"], d_basis, c3ref, params["K3_range"], params["measure"]) : 0.0

        println("Reference energy: ", state["reference_energy"])

        #reset the temperature for that resolution
        state["T"] = calculate_temperature(state, params, d_basis, c3ref)

        if params["include_negativity"]
            state["negativity_factor"] = log((state["E"]+state["T"])/state["E"])/0.3 #0.3 is expected initial negativity

            #With new negativity_factor, recalculate energy
            state["E"] = energy(state["intensity"], d_basis, c3ref, params["K3_range"], params["measure"], state["negativity_factor"])
        else
            state["negativity_factor"] = 0.0
        end

        #Let's acknowledge that higher expansion limits L require slower temperature decay
        state["temperature_decay"] = params["temperature_decay"]
    end

    #run until we have converged with the stepsize (stepsize for largest L must below threshold)
    while state["stepsizes"][state["L"]] > minstepsize

        #increase overall step counter
        state["i"] += 1

        #adjust temperature each step with T_new = decay * T_old
        state["T"] *= state["temperature_decay"]

        #Monte Carlo step
        state["step"] = get_MC_step(state["intensity"], d_basis, state["stepsizes"])

        #evaluate step
        @time evaluate_step(out, params, state, d_basis, c3ref)

        regular_action(out, state, params)
    end
end

# "Iterativly optimizing a structure by increasing the angular resolution and adding more shells that are taken into consideration"
# function rotate_hierarchical(out, params::Dict, state::Dict, c3ref::C3)
#
#     for state["L"] = state["L"]:2:params["L"]
#
#         d_basis = complexBasis_choice(state["L"], params["LMAX"], params["N"], maximum(params["K3_range"]), params["lambda"], dq(state["intensity"]))
#
#         #Check if this is a continuation or a new run by checking stepsizes
#         if !haskey(state["stepsizes"], state["L"])
#
#             #gradually increase K
#             state["K"] = maximum(params["K3_range"]) - params["L"] + state["L"]
#
#             #set stepsize for this resolution
#             state["stepsizes"] = Dict(l=>params["initial_stepsize"]/(state["L"]-l+1) for l=2:2:state["L"])
#
#             #calculate the initial energy with that resolution
#             state["E"] = energy(state["intensity"], d_basis, c3ref, params["measure"])
#
#             #Calculate reference energy
#             state["reference_energy"] = haskey(state, "reference_intensity") ? energy(state["reference_intensity"], d_basis, c3ref, params["measure"]) : 0.0
#
#             #reset the temperature for that resolution
#             state["T"] = calculate_temperature(state, params, d_basis, c3ref)
#
#             #Let's acknowledge that higher expansion limits L require slower temperature decay
#             state["temperature_decay"] = params["temperature_decay"]
#         end
#
#         #run until we have converged with the stepsize (stepsize for largest L must below threshold)
#         while state["stepsizes"][state["L"]] > minstepsize
#
#             #increase overall step counter
#             state["i"] += 1
#
#             #adjust temperature each step with T_new = decay * T_old
#             state["T"] *= state["temperature_decay"]
#
#             #Monte Carlo step
#             state["step"] = get_MC_step(state["intensity"], d_basis, state["stepsizes"])
#
#             #evaluate step
#             @time evaluate_step(out, params, state, d_basis, c3ref)
#
#             #perform regular action such as logging
#             regular_action(out, state, params)
#         end
#     end
# end

# "MC sampling based on integration of the Hamilton equations stemming from negative log-probability"
# function rotate_HMC(new, step, out, initial_stepsize, iterations, stepsizecheck, Tstart, L, klow, K, N, E, number)
#
#   #variables for MC annealing
#   # num_accepted_steps = 0
#   # num_MC_steps = 0
#   # stepsize = initial_stepsize
#   dim = 10
#   q = rand(dim)
#   p = rand(dim)
#
#   grad_p = function(p) return p end
#   grad_q = function(q)
#     return q
#   end
#
#   # global coeff = deepcopy(new)
#   # global L_optim = L
#   # numangles = sum([(2*l+1)*l for l=2:2:L])
#   # q = rand(numangles)
#   # p = rand(numangles)
#   # eps = initial_stepsize
#   # E = opt_energy(q)
#   #
#   # grad_p = function(p) return p end
#   # grad_q = function(q)
#   #     derivative = deepcopy(q)
#   #     opt_energy_derivative!(q, derivative)
#   #     return derivative
#   # end
#
#   propose_step = function(q,p, iterations)
#     p = p - 0.5 * eps * grad_q(q)
#
#     for j = 1:iterations
#         q = q + eps * grad_p(p)
#         p = p - eps * grad_q(q)
#     end
#
#     q = q + eps * grad_p(p)
#     p = p - 0.5 * eps * grad_q(q)
#     return q,p
#   end
#
#   for i = 1:iterations
#
#       # # Monitoring acceptance rate and modify stepsize accordingly
#       # if num_MC_steps == stepsizecheck
#       #     rate = num_accepted_steps/num_MC_steps
#       #     num_MC_steps = 0
#       #     num_accepted_steps = 0
#       #     if rate < 0.25
#       #         stepsize /= 2
#       #     # end
#       #     elseif rate > 0.75
#       #         stepsize *= 2
#       #     end
#       # end
#
#       new_q,new_p = propose_step(q,p,100)
#
#       E_new = opt_energy(new_q)
#       delta_E = E_new - E
#       # T =  Tstart*exp(-i/(iterations/10))
#       ap  = exp(-(delta_E)) #acceptance energy
#       accepted = false
#       sign = 0.0
#       #Check if the step improved the triple correlation energy
#       if E_new < E || ( E_new > E && rand() <  ap)
#
#           sign = E_new < E ? 1.0 : -1.0
#           accepted = true
#
#           #save the rotational information
#           E = E_new
#           q = deepcopy(new_q)
#       end
#
#       println("$i\t$sign\t$E_new\t$E\t$delta_E\t$T\t$ap\t$eps")
#
#       if accepted num_accepted_steps += 1 end
#
#       num_MC_steps += 1
#   end
#   return new
# end

"""Increases the stepsize for all L"""
function increase_stepsize(state::Dict, params::Dict)
    for l in keys(state["stepsizes"])
        state["stepsizes"][l] *= params["stepsizefactor"]
    end
end

"""Decreases the stepsize for all L"""
function decrease_stepsize(state::Dict, params::Dict)
    for l in keys(state["stepsizes"])
        state["stepsizes"][l] /= params["stepsizefactor"]
    end
end

function evaluate_step(out, params::Dict, state::Dict, basis::AbstractBasisType, c3ref::C3)
    accepted = false
    sign = 0.0

    #Calculate the triple correlation and return difference to reference triple correlation
    E_new = energy(state["step"], basis, c3ref, params["K3_range"], params["measure"], state["negativity_factor"])
    delta_E = E_new - state["E"]
    ap  = exp(-(delta_E)/state["T"]) #acceptance energy

    #Check if the step improved the triple correlation energy
    if delta_E < 0.0 || ( delta_E > 0.0 && rand() <  ap)
        sign = delta_E < 0.0 ? 1.0 : -1.0
        increase_stepsize(state, params)
        #save the rotational information
        state["E"] = E_new
        state["intensity"] = deepcopy(state["step"])
    else
        decrease_stepsize(state, params)
    end

    #Logging current step
    line = "$(state["i"])\t$sign\t$E_new\t$(state["E"])\t$delta_E\t$(state["reference_energy"])\t$(state["E"]-state["reference_energy"])\t$(state["T"])\t$ap\t$(state["stepsizes"])\t$(state["L"])\t$(state["K"])\n"
    write(out, line)
    println(line)
end

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

function checkRotationSearch(reference::SphericalHarmonicsVolume, K::Int64, L::Int64, basis::AbstractBasisType; histogram="../expdata/correlations_N32_K25_P2048000.dat", save_structures::Bool=false)
    c2full,c2,c3full,c3 = loadHistograms(K, K, histogram)
    start = retrieveSolution(c2full, L, reference.LMAX, reference.KMAX, reference.rmax)
    start = randomStartStructure(deleteTerms(start, K, L), K, L)
    checkRotationSearch(start, deleteTerms(reference, K, L), 1:K, L, basis, save_structures=save_structures)
end

"""Tries out the rotational search scheme by comparing given structures with the known original structure
This yields convergence within the bounds of the scheme"""
function checkRotationSearch(start::SphericalHarmonicsVolume, reference::SphericalHarmonicsVolume, K_range::UnitRange{Int64}, L::Int64, basis::AbstractBasisType; iterations=1.0e4, reduce_stepsize=1000, plotting=false, include_negativity=false, energy=(x)->0.0, save_structures::Bool=false)

    new = deepcopy(start)
    step = deepcopy(new)

    FSC_list = Float64[]
    E_list = Float64[]

    FSC = similarity(reference, start, K_range)
    stepsizes = Dict(l=>pi*(l-1) for l=2:2:L)

    E = energy(new)

    for i = 0:round(Int64,iterations)

        if mod(i, reduce_stepsize) == 0
            E = energy(new)
            push!(E_list, E)
            push!(FSC_list, FSC)

            for l in keys(stepsizes)
                stepsizes[l] /= 2.0
            end
            println("$i\tFSC=$FSC\tE=$E")
            if save_structures
                saveCube(new, "output/rotation_check_$i.mrc")
            end
        end

        #Modify the structure a little bit
        step = get_MC_step(new, basis, stepsizes)

        #Calculate FSC of this new structure
        FSC_new = similarity(reference, step, K_range)

        #Check if we accept this step
        if FSC_new > FSC
            FSC = FSC_new
            new = deepcopy(step)
        end
    end

    return new, FSC
end

# """Takes a set of coefficients and applies rotations to it"""
# function rotate_coefficients(coeff::SphericalHarmonicsVolume, L::Int64, angles::Array{Float64})
#   rotated_coeff = deepcopy(coeff)
#
#   #Rotate the coefficients by the given angles
#   for l = 2:2:L
#
#       num_l_angles = (2*l+1)*l
#       #Unpacking the angles from the vector and calculating the rotation matri
#       R = SON_parametrized(2*l+1,angles[ai:ai+num_l_angles-1])
#
#       #Rotating the individual coefficients for all shells
#       for k = 1:KMAX cvec_set(rotated_coeff,k,l, R*cvec_get(coeff, k,l)) end
#   end
#
#   return rotated_coeff
# end

# end


text:
#=
    x / 2pi
            x * tripleprecision(inv(2pi))
            intpart, fracpart
=#
function rem2pi(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}}
    hi, md, lo = mul233(HILO(x), inv_pi_2o1_t64)

    intfrac_hi = modf(hi)
    intfrac_md = modf(md)
    intfrac_lo = modf(lo)
    fracparts = intfrac_hi[1], intfrac_md[1], intfrac_lo[1]
    # fracpart is the signed (directed) portion of a circle
    hi, md, lo = mul333(fracparts, pi_2o1_t64)
    return T(hi, md)
end

function rem1pi(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}}
    hi, md, lo = mul233(HILO(x), inv_pi_1o1_t64)

    intfrac_hi = modf(hi)
    intfrac_md = modf(md)
    intfrac_lo = modf(lo)
    fracparts = intfrac_hi[1], intfrac_md[1], intfrac_lo[1]
    # fracpart is the signed (directed) portion of a circle
    hi, md, lo = mul333(fracparts, pi_1o1_t64)
    return T(hi, md)
end

"""
    negpi_pospi(x)

    x --> [-pi..+pi)
"""
function negpi_pospi(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}}
    result = rem2pi(x)
    if result >= T(pi)
        result = -rem2pi(-x)
    end
    return result
end

# in -pi..+pi
# determine nearest multiple of pi/2
# within (pi/2) quadrant determine nearest multiple of pi/16

function whichquadrant(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}}
    if signbit(x) # quadrant -2 or -1
        quadrant = -x < T(pi)/2 ? -2 : -1
    else          # quadrant  1 or  2
        quadrant =  x < T(pi)/2 ? 1 :  2
    end
    return quadrant
end

function which64th(x::DoubleFloat{Float64}, quadrant::Int)
    y = T(mul232(HILO(x), nv_pi_1o2_t64))
    z = trunc(Int, y*4)
    return z
end
function which64th(x::DoubleFloat{Float32}, quadrant::Int)
    y = T(mul232(HILO(x), nv_pi_1o2_t32))
    z = trunc(Int, y*4)
    return z
end
function which64th(x::DoubleFloat{Float16}, quadrant::Int)
    y = T(mul232(HILO(x), nv_pi_1o2_t16))
    z = trunc(Int, y*4)
    return z
end

# ========================================== ^^^^^ =========================================

const pi_1o1_t64_hi = pi_1o1_t64[1]
const pi_1o1_t64_md = pi_1o1_t64[2]
const pi_1o1_t64_lo = pi_1o1_t64[3]

const pi_1o4_t64_hi = pi_1o4_t64[1]
const pi_1o4_t64_md = pi_1o4_t64[2]
const pi_1o4_t64_lo = pi_1o4_t64[3]

function value_minus_pi(x::DoubleFloat{Float64})
    hi, lo = sub232(HI(x), LO(x), pi_1o1_t64_hi, pi_1o1_t64_md, pi_1o1_t64_lo)
    return DoubleFloat{Float64}(hi, lo)
end

function value_minus_qrtrpi(x::DoubleFloat{Float64})
    hi, lo = sub232(HI(x), LO(x), pi_1o4_t64_hi, pi_1o4_t64_md, pi_1o4_t64_lo)
    return DoubleFloat{Float64}(hi, lo)
end

function mod2pi(x::DoubleFloat{Float64})
    signbit(x) && return mod2pi_neg(x)
    x < twopi_d64 && return x
    w1 = mul323(inv_pi_2o1_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_2o1_t64, w2)
    z = Double64(y)
    return z
end

function mod2pi_neg(x::DoubleFloat{Float64})
    m = mod2pi(-x)
    return Double64(sub322(pi_2o1_t64, HILO(m)))
end

mod2pi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(mod2pi(DoubleFloat{Float64}(x)))
mod2pi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(mod2pi(DoubleFloat{Float64}(x)))

function modpi(x::DoubleFloat{Float64})
    signbit(x) && return modpi_neg(x)
    x < onepi_d64 && return x
    w1 = mul323(inv_pi_1o1_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_1o1_t64, w2)
    z = Double64(y)
    return z
end

function modpi_neg(x::DoubleFloat{Float64})
    m = modpi(-x)
    return Double64(sub322(pi_1o1_t64, HILO(m)))
end
modpi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modpi(DoubleFloat{Float64}(x)))
modpi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modpi(DoubleFloat{Float64}(x)))

function modpipm(x::DoubleFloat{Float64})
    abs(x) < onepi_d64 && return x
    w1 = mul323(inv_pi_1o1_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_1o1_t64, w2)
    z = Double64(y)
    return z
end

modpipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modpipm(DoubleFloat{Float64}(x)))
modpipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modpipm(DoubleFloat{Float64}(x)))

function modhalfpi(x::DoubleFloat{Float64})
    signbit(x) && return modhalfpi_neg(x)
    x < halfpi_d64 && return x
    w1 = mul323(inv_pi_1o2_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_1o2_t64, w2)
    z = Double64(y)
    return z
end

function modhalfpi_neg(x::DoubleFloat{Float64})
    m = modhalfpi(-x)
    return Double64(sub322(pi_1o2_t64, HILO(m)))
end

modhalfpi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modhalfpi(DoubleFloat{Float64}(x)))
modhalfpi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modhalfpi(DoubleFloat{Float64}(x)))

function modhalfpipm(x::DoubleFloat{Float64})
    abs(x) < halfpi_d64 && return x
    w1 = mul323(inv_pi_1o2_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_1o2_t64, w2)
    z = Double64(y)
    return z
end

modhalfpipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modhalfpipm(DoubleFloat{Float64}(x)))
modhalfpipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modhalfpipm(DoubleFloat{Float64}(x)))

function modqrtrpi(x::DoubleFloat{Float64})
    signbit(x) && return modqrtrpi_neg(x)
    x < qrtrpi_d64 && return x
    w1 = mul323(inv_pi_1o4_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_1o4_t64, w2)
    z = Double64(y)
    return z
end

function modqrtrpi_neg(x::DoubleFloat{Float64})
    m = modqrtrpi(-x)
    return Double64(sub322(pi_1o4_t64, HILO(m)))
end

modqrtrpi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modqrtrpi(DoubleFloat{Float64}(x)))
modqrtrpi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modqrtrpi(DoubleFloat{Float64}(x)))

function modqrtrpipm(x::DoubleFloat{Float64})
    abs(x) < qrtrpi_d64 && return x
    w1 = mul323(inv_pi_1o4_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_1o4_t64, w2)
    z = Double64(y)
    return z
end

modqrtrpipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modqrtrpipm(DoubleFloat{Float64}(x)))
modqrtrpipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modqrtrpipm(DoubleFloat{Float64}(x)))

#=
    rem2pi(x) =  x - 2π*round(x/(2π),r)
    rem2pi(x, RoundDown) == mod2pi(x)
    •    if r == RoundNearest, then the result is in the interval [-π, π]. This will
        generally be the most accurate result.
    •    if r == RoundToZero, then the result is in the interval [0, 2π] if x is positive,.
        or [-2π, 0] otherwise.
    •    if r == RoundDown, then the result is in the interval [0, 2π].
    •    if r == RoundUp, then the result is in the interval [-2π, 0].
=#

rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Down}) = mod2pi(x)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Up}) = -rem2pi(-x, RoundDown)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest}) = modpipm(x)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:ToZero}) =
    signbit(x) ? rem2pi(x, RoundUp) : rem2pi(x, RoundDown)

rem2pi(x::DoubleFloat{Float32}, rounding::RoundingMode) =
    DoubleFloat{Float32}(rem2pi(DoubleFloat{Float64}(x), rounding))
rem2pi(x::DoubleFloat{Float16}, rounding::RoundingMode) =
    DoubleFloat{Float16}(rem2pi(DoubleFloat{Float64}(x), rounding))

rempi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Down}) = modpi(x)
rempi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Up}) = -rempi(-x, RoundDown)
rempi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest}) = modhalfpipm(x)
rempi(x::DoubleFloat{Float64}, rounding::RoundingMode{:ToZero}) =
    signbit(x) ? rempi(x, RoundUp) : rempi(x, RoundDown)

rempi(x::DoubleFloat{Float32}, rounding::RoundingMode) =
    DoubleFloat{Float32}(rempi(DoubleFloat{Float64}(x), rounding))
rempi(x::DoubleFloat{Float16}, rounding::RoundingMode) =
    DoubleFloat{Float16}(rempi(DoubleFloat{Float64}(x), rounding))

remhalfpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Down}) = modhalfpi(x)
remhalfpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Up}) = -remhalfpi(-x, RoundDown)
remhalfpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest}) = modqrtrpipm(x)
remhalfpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:ToZero}) =
    signbit(x) ? remhalfpi(x, RoundUp) : remhalfpi(x, RoundDown)

remhalfpi(x::DoubleFloat{Float32}, rounding::RoundingMode) =
    DoubleFloat{Float32}(remhalfpi(DoubleFloat{Float64}(x), rounding))
remhalfpi(x::DoubleFloat{Float16}, rounding::RoundingMode) =
    DoubleFloat{Float16}(remhalfpi(DoubleFloat{Float64}(x), rounding))

remqrtrpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Down}) = modqrtrpi(x)
remqrtrpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Up}) = -remqrtrpi(-x, RoundDown)
# remqrtrpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest}) = modeighthpipm(x)
remqrtrpi(x::DoubleFloat{Float64}, rounding::RoundingMode{:ToZero}) =
    signbit(x) ? remqrtrpi(x, RoundUp) : remqrtrpi(x, RoundDown)

remqrtepi(x::DoubleFloat{Float32}, rounding::RoundingMode) =
    DoubleFloat{Float32}(remqrtpi(DoubleFloat{Float64}(x), rounding))
remqrtrpi(x::DoubleFloat{Float16}, rounding::RoundingMode) =
DoubleFloat{Float16}(remqrtrpi(DoubleFloat{Float64}(x), rounding))


text:
#=
global:
- Julia version:
- Author: Dan
- Date: 2021-06-02
=#
#not the same implementation as in C++BALL

export
	Byte, Position, Handle, ASCII, Amino_Acids,
	BALL_QSAR_RINGPERCEEPTIONPROCESSOR_RUN_COUNT,
	BALL_QSAR_RINGPERCEEPTIONPROCESSOR_MAX_RUNS,
	BALL_HALF_OF_MAX_RING_SIZE, BALL_Properties

const BALL_SIZE_TYPE = Int64
const BALL_INDEX_TYPE = Int64

const Byte = UInt8
const Position = BALL_SIZE_TYPE
const Handle = BALL_SIZE_TYPE
#const Size = BALL_SIZE_TYPE
const Distance = BALL_INDEX_TYPE
const Index = BALL_INDEX_TYPE

#this being an enum on one hand makes sense, because we want to store this info in bytes (UInt8s)
#but it also doesnt because we want it to be interpreted as a Char - requiring to convert it every time youw ant to use it
@enum ASCII begin

		ASCII__BACKSPACE        = convert(UInt8,'\b')
		ASCII__BELL             = convert(UInt8,'\a')
		ASCII__CARRIAGE_RETURN  = convert(UInt8,'\r')
		ASCII__HORIZONTAL_TAB   = convert(UInt8,'\t')
		ASCII__NEWLINE          = convert(UInt8,'\n')
		ASCII__SPACE            = convert(UInt8,' ')
		ASCII__VERTICAL_TAB     = convert(UInt8,'\v')

		ASCII__COLON            = convert(UInt8,':')
		ASCII__COMMA            = convert(UInt8,',')
		ASCII__EXCLAMATION_MARK = convert(UInt8,'!')
		ASCII__POINT            = convert(UInt8,'.')
		ASCII__QUESTION_MARK    = convert(UInt8,'?')
		ASCII__SEMICOLON        = convert(UInt8,';')
	end
		ASCII__RETURN           = ASCII__NEWLINE
		ASCII__TAB              = ASCII__HORIZONTAL_TAB
		#whenever one tries to convert an ascii to any type, it will actually convert to a character
		Base.convert(::Type{Any}, x::ASCII) = convert(Char, x)

const INVALID_Distance = typemax(Distance)
const Distance_MIN = typemin(Distance) + 1
const Distance_MAX = typemax(Distance)

const INVALID_Handle = typemax(Handle)
const Handle_MIN = 0 ;
const Handle_MAX = typemax(Handle) - 1

const INVALID_Index = -1;
const Index_MIN = 0;
const Index_MAX = typemax(BALL_INDEX_TYPE)

const INVALID_Position = typemax(Position)
const Position_MIN = 0
const Position_MAX = typemax(BALL_INDEX_TYPE) - 1

#const INVALID_Size = typemax(Size)
#const Size_MIN = 0;
#const Size_MAX = typemax(Size) - 1;
"A String-Vector of all Amino Acid abbreviations."
const Amino_Acids = String["ALA", "CYS", "ASP", "GLU", "PHE", "GLY", "HIS", "ILE", "LYS", "LEU",
 							"MET", "ASN", "PRO", "GLN", "ARG", "SER", "THR", "VAL", "TRP", "TYR"]
BALL_QSAR_RINGPERCEEPTIONPROCESSOR_RUN_COUNT = 1
BALL_QSAR_RINGPERCEEPTIONPROCESSOR_MAX_RUNS = 20
"A String-Vector of all available Properties. Note: Properties can be added at runtime. (Needs testing tho)"
BALL_Properties = String["amino_acid", "InRing","fromNMR"]
"Ball can only percieve Rings of Size 2*`BALL_HALF_OF_MAX_RING_SIZE`"
const BALL_HALF_OF_MAX_RING_SIZE = 20		#Only detect Rings in SSSR with maximally 2*BALL_HALF_OF_MAX_RING_SIZE atoms


text:
SYNSET_TYPES = Dict{Char, String}(
    'n' => "noun", 'v' => "verb", 'a' => "adj", 'r' => "adv"
)

NOUN_POINTERS = Dict{String, String}(
    "-c" => "Member of this domain - TOPIC",
    "+" => "Derivationally related form",
    "%p" => "Part meronym",
    "~i" => "Instance Hyponym",
    "@" => "Hypernym",
    ";r" => "Domain of synset - REGION",
    "!" => "Antonym",
    "#p" => "Part holonym",
    "%s" => "Substance meronym",
    ";u" => "Domain of synset - USAGE",
    "-r" => "Member of this domain - REGION",
    "#s" => "Substance holonym",
    "=" => "Attribute",
    "-u" => "Member of this domain - USAGE",
    ";c" => "Domain of synset - TOPIC",
    "%m" => "Member meronym",
    "~" => "Hyponym",
    "@i" => "Instance Hypernym",
    "#m" => "Member holonym"
)
VERB_POINTERS = Dict{String, String}(
    "+" => "Derivationally related form",
    "@" => "Hypernym",
    ";r" => "Domain of synset - REGION",
    "!" => "Antonym",
    ";u" => "Domain of synset - USAGE",
    "\$" => "Verb Group",
    ";c" => "Domain of synset - TOPIC",
    ">" => "Cause",
    "~" => "Hyponym",
    "*" => "Entailment"
)

ADJECTIVE_POINTERS = Dict{String, String}(
    ";r" => "Domain of synset - REGION",
    "!" => "Antonym",
    "\\" => "Pertainym (pertains to noun)",
    "<" => "Participle of verb",
    "&" => "Similar to",
    "=" => "Attribute",
    ";c" => "Domain of synset - TOPIC"
)

ADVERB_POINTERS = Dict{String, String}(
    ";r" => "Domain of synset - REGION",
    "!" => "Antonym",
    ";u" => "Domain of synset - USAGE",
    "\\" => "Derived from adjective",
    ";c" => "Domain of synset - TOPIC"
)

MEMBER_OF_THIS_DOMAIN_TOPIC = "-c"
DERIVATIONALLY_RELATED_FORM = "+"
PART_MERONYM = "%p"
InstanceHyponym = "~i"
HYPERNYM = "@"
DOMAIN_OF_SYNSET_REGION = ";r"
ANTONYM = "!"
PART_HOLONYM = "#p"
SUBSTANCE_MERONYM = "%s"
VERB_GROUP = "\$"
DOMAIN_OF_SYNSET_USAGE = ";u"
MEMBER_OF_THIS_DOMAIN_REGION = "-r"
SUBSTANCE_HOLONYM = "#s"
DERIVED_FROM_ADJECTIVE = "\\"
PARTICIPLE_OF_VERB = "<"
SIMILAR_TO = "&"
ATTRIBUTE = "="
ALSO_SEE = "^"
CAUSE = ">"
MEMBER_OF_THIS_DOMAIN_USAGE = "-u"
DOMAIN_OF_SYNSET_TOPIC = ";c"
MEMBER_MERONYM = "%m"
HYPONYM = "~"
INSTANCE_HYPERNYM = "@i"
ENTAILMENT = "*"
MEMBER_HOLONYM = "#m"


text:
import LinearAlgebra: qr, lu, eigen, diag

function feast!(X::AbstractMatrix, A::AbstractMatrix;
				nodes::Integer=8, iter::Integer=10, c=complex(0.0,0.0), r=1.0, ϵ=1e-12,
				debug=false, store=false, mixed_prec=false, factorizer=lu, left_divider=ldiv!)
	contour = circular_contour_trapezoidal(c, r, nodes)
	feast!(X, A, contour; iter=iter, debug=debug, ϵ=ϵ, store=store, mixed_prec=mixed_prec, factorizer=factorizer, left_divider=left_divider)
end

finalize!(x::Any) = nothing

function feast!(X::AbstractMatrix, A::AbstractMatrix, contour::Contour;
					 iter::Integer=10, ϵ=1e-12, debug=false, store=false, mixed_prec=false, factorizer=lu, left_divider=ldiv!)
	N, m₀ = size(X)
	if size(A, 1) != size(A, 2)
		 error("Incorrect dimensions of A, must be square")
	elseif size(A,1) != N
		 error("Incorrect dimensions of X, must match A")
	end

	Ctype = if mixed_prec ComplexF32 else ComplexF64 end

	Λ, resolvent, res = zeros(ComplexF64, m₀), zeros(ComplexF64, m₀), zeros(m₀)
	temp, R, Q = zeros(Ctype, N, m₀), similar(X, ComplexF64), deepcopy(X)
	Aq, Xq = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)

	ZmA = similar(A, Ctype)
	nodes = size(contour.nodes, 1)

	if store
		ZmA .= (A - I*contour.nodes[1])
		facts1 = factorizer(ZmA)
		facts = Array{typeof(facts1)}(undef, nodes)
		facts[1] = facts1

		Threads.@threads for i=2:nodes
			  ZmA .= (A - I*contour.nodes[i])
			  facts[i] = factorizer(ZmA)
		end
	end

	for nit=0:iter
		Q .= Matrix(qr(Q).Q)
		mul!(R, A, Q) ## why does this one allocate?
		mul!(Aq, Q', R) ### Aq = Q' * A * Q
		# mul!(Bq, Q', Q) ### Bq = Q' * Q = I
		F = eigen!(Aq)
		Λ .= F.values
		Xq .= F.vectors
		mul!(X, Q, Xq) ### Recover eigenvectors from Ritz vectors ( X = Q * Xq )
		update_R!(X, R, Λ, A) ### compute residual vectors R for RII update
		residuals!(res, R, Λ, A) ### compute actual residuals
		contour_nonempty = reduce(|, in_contour(Λ, contour))
		if debug iter_debug_print(nit, Λ, res, contour, 1e-5) end
		if contour_nonempty && maximum(res[in_contour(Λ, contour)]) < ϵ
			if debug println("converged in $nit iteration") end
			break
		end
		if nit < iter ### Do not solve linear systems / form Q on last iteration
			Q .= 0.00
			for i=1:nodes
				resolvent .= 1.0 ./(contour.nodes[i] .- Λ)
				if store
					left_divider(temp, facts[i], R)
				else
					ZmA .= A - I*contour.nodes[i]
					linsolve!(temp, ZmA, R, factorizer, left_divider)
				end

				temp .= X - temp
				rmul!(temp, Diagonal(resolvent .* contour.weights[i]))
				Q .+= temp
			end
		end
	end
	if store
		foreach(finalize!, facts)
	end
	contour_nonempty = reduce(|, in_contour(Λ, contour))
	if !contour_nonempty println("no eigenvalues found in contour!") end
	Λ[in_contour(Λ, contour)], X[:,in_contour(Λ, contour)], res[in_contour(Λ, contour)]
end

function gen_feast!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix;
					nodes::Integer=8, iter::Integer=10, c=complex(0.0,0.0), r=1.0,
					debug=false, store=false, ϵ=1e-12, factorizer=lu, left_divider=ldiv!)
	contour = circular_contour_trapezoidal(c, r, nodes)
	gen_feast!(X, A, B, contour; iter=iter, debug=debug, ϵ=ϵ, factorizer=factorizer, left_divider=ldiv!)
end

function gen_feast!(X::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, contour::Contour;
					iter::Integer=10, debug=false, store=false, ϵ=1e-12, factorizer=lu, left_divider=ldiv!)
	N, m₀ = size(X)
	if size(A, 1) != size(A, 2)
		error("Incorrect dimensions of A, must be square")
	elseif size(A,1) != N
		error("Incorrect dimensions of X, must match A")
	end

	Λ, resolvent, res = zeros(ComplexF64, m₀), zeros(ComplexF64, m₀), zeros(m₀)
	temp, R, Q = zeros(ComplexF64, N, m₀), similar(X, ComplexF64), copy(X)
	Aq, Bq, Xq = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)
	ZmA = similar(A, ComplexF64)
	nodes = size(contour.nodes, 1)

	if store
		 ZmA .= (A - B*contour.nodes[1])
		 facts1 = factorizer(ZmA)
		 facts = Array{typeof(facts1)}(undef, nodes)
		 facts[1] = facts1

		 Threads.@threads for i=2:nodes
				ZmA .= (A - B*contour.nodes[i])
				facts[i] = factorizer(ZmA)
		 end
	end

	for nit=0:iter
		Q .= Matrix(qr(Q).Q)
		mul!(R, A, Q) ## why does this one allocate?
		mul!(Aq, Q', R) ### Aq = Q' * A * Q
		mul!(R, B, Q)
		mul!(Bq, Q', R) ### Bq = Q' * Q = I
		F = eigen!(Aq, Bq)
		Λ .= F.values
		Xq .= F.vectors
		mul!(X, Q, Xq) ### Recover eigenvectors from Ritz vectors ( X = Q * Xq )
		update_R!(X, R, Λ, A, B) ### compute residual vectors R for RII update
		residuals!(res, R, Λ, A) ### compute actual residuals
		contour_nonempty = reduce(|, in_contour(Λ, contour))
		if debug iter_debug_print(nit, Λ, res, contour, 1e-5) end
		if contour_nonempty && maximum(res[in_contour(Λ, contour)]) < ϵ
			  if debug println("converged in $nit iteration") end
			  break
		end
		if nit < iter ### Do not solve linear systems / form Q on last iteration
			Q .= 0.00
			for i=1:nodes
				resolvent .= (1.0 ./(contour.nodes[i] .- Λ))
				if store
					left_divider(temp, facts[i], R)
				else
					ZmA .= (A - B*contour.nodes[i])
					linsolve!(temp, ZmA, R, factorizer, left_divider)
				end
				temp .= X - temp
				rmul!(temp, Diagonal(resolvent .* contour.weights[i]))
				Q .+= temp
			end
		end
	end
	if store
		foreach(finalize!, facts)
	end
	contour_nonempty = reduce(|, in_contour(Λ, contour))
	if !contour_nonempty println("no eigenvalues found in contour!") end
	Λ[in_contour(Λ, contour)], X[:,in_contour(Λ, contour)], res[in_contour(Λ, contour)]
end

function dual_gen_feast!(Xr::AbstractMatrix, Xl::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix;
					nodes::Integer=8, iter::Integer=10, c=complex(0.0,0.0), r=1.0,
					debug=false, store=false, ϵ=1e-12, factorizer=lu, left_divider=ldiv!)
	contour = circular_contour_trapezoidal(c, r, nodes)
	dual_gen_feast!(Xr, Xl, A, B, contour; iter=iter, debug=debug, ϵ=ϵ, factorizer=factorizer, left_divider=ldiv!)
end

function dual_gen_feast!(Xr::AbstractMatrix, Xl::AbstractMatrix, A::AbstractMatrix, B, contour::Contour;
					iter::Integer=10, debug=false, store=false, ϵ=1e-12, factorizer=lu, left_divider=ldiv!, resize_subspace=true)
	N, m₀ = size(Xl)
	if size(A, 1) != size(A, 2)
		error("Incorrect dimensions of A, must be square")
	elseif size(A,1) != N
		error("Incorrect dimensions of X, must match A")
	end

	Λ, resolvent, resr, resl, temp, Rr, Rl, Ql, Qr, Aq, Bq, Xq = init(Xr, Xl, N, m₀)
	ZmA = similar(A, ComplexF64)
	nodes = size(contour.nodes, 1)

	if store
		if debug print("Precomputing factorizations: ") end
			ZmA .= (A - B*contour.nodes[1])
			rfacts1 = factorizer(ZmA)
			rfacts = Array{typeof(rfacts1)}(undef, nodes)
			rfacts[1] = rfacts1
			if debug print("1R ") end
			ZmA .= (A - B*contour.nodes[1])'
			lfacts1 = factorizer(ZmA)
			lfacts = Array{typeof(lfacts1)}(undef, nodes)
			lfacts[1] = lfacts1
			if debug print("1L ") end
			for i=2:nodes
				ZmA .= (A - B*contour.nodes[i])
				rfacts[i] = factorizer(ZmA)
				if debug print("$(i)R ") end
				ZmA .= (A - B*contour.nodes[i])'
				lfacts[i] = factorizer(ZmA)
				if debug print("$(i)L ") end
			end
		if debug println() end
	end

	try
	nin = 0
	for nit=0:iter
		if resize_subspace
			new_m₀ = new_subspace_size(nin, m₀)
			if new_m₀ > 0 && new_m₀ != m₀
				if debug println("Resizing subspace from $(m₀) to $(new_m₀)") end
				m₀ = new_m₀
				Λ, Xr, Xl, resolvent, resr, resl, temp, Rr, Rl, Ql, Qr, Aq, Bq, Xq = reinit(Λ, Xr, Xl, resolvent, resr, resl, temp, Rr, Rl, Ql, Qr, Aq, Bq, Xq, new_m₀)
			end
		end
		S = svd!(Ql'*B*Qr)
		Qr .= Qr*S.V*Diagonal(1.0/S.S)
		Ql .= Ql*S.U*Diagonal(1.0/S.S)
		mul!(Rr, A, Qr) ## why does this one allocate?
		mul!(Aq, Ql', Rr) ### Aq = Q' * A * Q
		mul!(Rr, B, Qr)
		mul!(Bq, Ql', Rr) ### Bq = Q' * Q = I
		F = eigen(Aq, Bq)
		Λ .= F.values
		Xq .= F.vectors
		mul!(Xr, Qr, Xq) ### Recover eigenvectors from Ritz vectors ( X = Q * Xq )
		F = eigen(Aq', Bq')
		Xq .= F.vectors
		mul!(Xl, Ql, Xq) ### Recover eigenvectors from Ritz vectors ( X = Q * Xq )
		update_R!(Xr, Rr, Λ, A, B) ### compute residual vectors R for RII update
		update_R!(Xl, Rl, Λ, A', B') ### compute residual vectors R for RII update
		residuals!(resr, Rr, Λ, A) ### compute actual residuals
		incntr = in_contour(Λ, contour)
		nin = sum(resr[incntr] .< 1e-5)
		contour_nonempty = reduce(|, incntr)
		if debug iter_debug_print(nit, Λ, resr, contour, 1e-5) end
		if contour_nonempty && maximum(resr[in_contour(Λ, contour)]) < ϵ
			if debug println("converged in $nit iteration") end
			break
		end
		if nit < iter ### Do not solve linear systems / form Q on last iteration
			Qr .= 0.00
			Ql .= 0.00
			for i=1:nodes
				resolvent .= (1.0 ./(contour.nodes[i] .- Λ))
				if store
					left_divider(temp, rfacts[i], Rr)
				else
					ZmA .= (A - B*contour.nodes[i])
					linsolve!(temp, ZmA, Rr, factorizer, left_divider)
				end
				temp .= Xr - temp
				rmul!(temp, Diagonal(resolvent .* contour.weights[i]))
				Qr .+= temp

				resolvent .= (1.0 ./(conj(contour.nodes[i]) .- conj.(Λ)))
				if store
					left_divider(temp, lfacts[i], Rl)
				else
					ZmA .= (A - B*contour.nodes[i])'
					linsolve!(temp, ZmA, Rl, factorizer, left_divider)
				end
				temp .= Xl - temp
				rmul!(temp, Diagonal(resolvent .* conj(contour.weights[i])))
				Ql .+= temp
			end
		end
	end

	contour_nonempty = reduce(|, in_contour(Λ, contour))
	if !contour_nonempty
		println("no eigenvalues found in contour!")
	end
	output = Λ[in_contour(Λ, contour)], Xr[:,in_contour(Λ, contour)], Xl[:,in_contour(Λ, contour)], resr[in_contour(Λ, contour)]
	return output
	finally
		if store
			foreach(finalize!, rfacts)
			foreach(finalize!, lfacts)
		end
	end
end

function linsolve!(Y, C, X, factorizer, left_divider)
	F = factorizer(C)
	left_divider(Y, F, X)
	finalize!(F)
end

function init(Xr, Xl, N::Int, m₀::Int)
	Λ, resolvent, resr, resl = zeros(ComplexF64, m₀), zeros(ComplexF64, m₀), zeros(m₀), zeros(m₀)
	temp, Rr, Rl, Ql, Qr = zeros(ComplexF64, N, m₀), similar(Xr, ComplexF64), similar(Xr, ComplexF64), copy(Xl), copy(Xr)
	Aq, Bq, Xq = zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀), zeros(ComplexF64, m₀, m₀)
	return Λ, resolvent, resr, resl, temp, Rr, Rl, Ql, Qr, Aq, Bq, Xq
end

function reinit(Λ, Xr, Xl, resolvent, resr, resl, temp, Rr, Rl, Ql, Qr, Aq, Bq, Xq, new_m₀, fillf::F=x->randn(typeof(x))) where F
	N, m₀ = size(temp,1), length(Λ)
	Xrn, Xln = zeros(ComplexF64, N, new_m₀), zeros(ComplexF64, N, new_m₀)
	map!(fillf, Xrn, Xrn)
	map!(fillf, Xln, Xln)
	copyto!(Xrn, Xr)
	copyto!(Xln, Xl)
	Λn, resolventn, resrn, resln, tempn, Rrn, Rln, Qln, Qrn, Aqn, Bqn, Xqn = init(Xrn, Xln, N, new_m₀)
	copyto!(Λn, Λ)
	copyto!(resolventn, resolvent)
	copyto!(resrn, resr)
	copyto!(resln, resl)
	copyto!(tempn, temp)
	copyto!(Rrn, Rr)
	copyto!(Rln, Rl)
	copyto!(Qln, Ql)
	copyto!(Qrn, Qr)
	copyto!(Aqn, Aq)
	copyto!(Bqn, Bq)
	copyto!(Xqn, Xqn)
	return Λn, Xrn, Xln, resolventn, resrn, resln, tempn, Rrn, Rln, Qln, Qrn, Aqn, Bqn, Xqn
end

function new_subspace_size(nin, m₀, fac_test=0.75, fac_new=2.0)
	if nin > ceil(Int, fac_test*m₀)
		return ceil(Int, fac_new*nin)
	else
		return 0
	end
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 5dc3b4c6-f4b9-11ea-3245-6d726e7e013f
begin
	using Images
	philip = load("philip.jpg")
end

# ╔═╡ f9893e7c-f4b8-11ea-22ab-7b9845fbe9f3
url = "https://i.imgur.com/VGPeJ6s.jpg"

# ╔═╡ 4e582dbe-f4b9-11ea-29c6-8d62ded32a99
download(url, "philip.jpg")

# ╔═╡ a28945a2-f4ba-11ea-1ae0-abcdbc83b171
typeof(philip)

# ╔═╡ cffcfb96-f4ba-11ea-1cd7-35a30f8077a7
RGBX(0.4, 0.6, 0.6)

# ╔═╡ e59de58c-f4ba-11ea-157e-096eed92ecfb
size(philip)

# ╔═╡ f35c83fe-f4ba-11ea-08c1-7515d7146923
philip[1:1000,1:400]

# ╔═╡ 1c3d9b0c-f4bb-11ea-1638-3d2076c23093
begin
	(h, w) = size(philip)
	head = philip[(h ÷ 2):h, (w ÷ 10): (9w ÷ 10)]
end

# ╔═╡ 567b6ffe-f4bb-11ea-35c7-d5e123d6cb40
size(head)

# ╔═╡ 644b4988-f4bb-11ea-1d5d-73900abda3d2
typeof([head head])

# ╔═╡ 7745ad30-f4bb-11ea-107d-b16cb31cfcff
[
	head 				  reverse(head, dims=2)
	reverse(head, dims=1) reverse(reverse(head,dims=1), dims=2)
]

# ╔═╡ d1f9c798-f4bb-11ea-0b58-09ef4a2ad157
new_phil = copy(head)

# ╔═╡ e7d2f878-f4bb-11ea-20be-91a92d4b1174
red = RGB(1,0,0)

# ╔═╡ f0b55c9c-f4bb-11ea-0300-c9fae6038d52
for i in 1:100
	for j in 1:300
		new_phil[i,j] = red
	end
end

# ╔═╡ 083eb4d0-f4bc-11ea-23c3-53c7e93e3d2c
new_phil

# ╔═╡ 0aaffd78-f4bc-11ea-0c9d-07264a6d052b
begin
	new_phil2 = copy(new_phil)
	new_phil2[100:200, 1:100] .= RGB(0,1,0)
	new_phil2
end

# ╔═╡ 3a6bf530-f4bc-11ea-391e-8772d38656c0
function redify(color)
	return RGB(color.r, 0, 0)
end

# ╔═╡ 49c8e54c-f4bc-11ea-0ce9-7fa4ffd7d818
begin
	color = RGB(0.5, 0.5, 0.2)

	[color, redify(color)]
end

# ╔═╡ 5cef3af4-f4bc-11ea-1d40-cde5ec741881
redify.(philip)

# ╔═╡ 674acdf6-f4bc-11ea-1bfd-9fbfa76c82c5
begin
	poor_phil = decimate(new_phil, 5)
	size(poor_phil)
end

# ╔═╡ 02fcc358-f4bd-11ea-300f-d77c789739fc
convolve(poor_phil, blur(2))

# ╔═╡ Cell order:
# ╠═f9893e7c-f4b8-11ea-22ab-7b9845fbe9f3
# ╠═4e582dbe-f4b9-11ea-29c6-8d62ded32a99
# ╠═5dc3b4c6-f4b9-11ea-3245-6d726e7e013f
# ╠═a28945a2-f4ba-11ea-1ae0-abcdbc83b171
# ╠═cffcfb96-f4ba-11ea-1cd7-35a30f8077a7
# ╠═e59de58c-f4ba-11ea-157e-096eed92ecfb
# ╠═f35c83fe-f4ba-11ea-08c1-7515d7146923
# ╠═1c3d9b0c-f4bb-11ea-1638-3d2076c23093
# ╠═567b6ffe-f4bb-11ea-35c7-d5e123d6cb40
# ╠═644b4988-f4bb-11ea-1d5d-73900abda3d2
# ╠═7745ad30-f4bb-11ea-107d-b16cb31cfcff
# ╠═d1f9c798-f4bb-11ea-0b58-09ef4a2ad157
# ╠═e7d2f878-f4bb-11ea-20be-91a92d4b1174
# ╠═f0b55c9c-f4bb-11ea-0300-c9fae6038d52
# ╠═083eb4d0-f4bc-11ea-23c3-53c7e93e3d2c
# ╠═0aaffd78-f4bc-11ea-0c9d-07264a6d052b
# ╠═3a6bf530-f4bc-11ea-391e-8772d38656c0
# ╠═49c8e54c-f4bc-11ea-0ce9-7fa4ffd7d818
# ╠═5cef3af4-f4bc-11ea-1d40-cde5ec741881
# ╠═674acdf6-f4bc-11ea-1bfd-9fbfa76c82c5
# ╠═02fcc358-f4bd-11ea-300f-d77c789739fc

max_stars_repo_name: ciriarte/learning

max_issues_repo_name: ciriarte/learning

max_forks_repo_name: ciriarte/learning


text:
module NonDimensional

using LinearAlgebra: norm

"""
    dynamic_pressure(ρ, V)

Compute the dynamic pressure given density ``ρ`` and speed ``V``.
"""
dynamic_pressure(ρ, V) = 1/2 * ρ * V^2

"""
    force_coefficient(force, q, S)

Compute the non-dimensional force coefficient given force, dynamic pressure ``q``, and area ``S``.
"""
force_coefficient(force, q, S) = force / (q * S)

"""
    moment_coefficient(moment, q, S, L)

Compute the non-dimensional moment coefficient given force, dynamic pressure ``q``, area ``S``, and reference length ``L``.
"""
moment_coefficient(moment, q, S, L) = force_coefficient(moment, q, S) / L

"""
    rate_coefficient(Ω, V, L)

Compute the non-dimensional angular velocity coefficient given angular speed ``Ω``, reference speed ``V``, and reference length ``L``.
"""
rate_coefficient(Ω, V, L) = Ω * L / 2V

"""
    pressure_coefficient(force, ρ, V, S)

Compute the pressure coefficient given force, reference density ``ρ``, speed ``V`` and area ``S``.
"""
pressure_coefficient(force, ρ, V, S) = force_coefficient(force, dynamic_pressure(ρ, V), S)

pressure_coefficient(mag, vels) = 1 - (norm(vels) / mag)^2

moment_coefficient(moment, q, S, b, c) = moment_coefficient.(moment, q, S, [b, c, b])

rate_coefficient(Ω, V, b, c) = rate_coefficient.(Ω, V, [b, c, b])

"""
    aerodynamic_coefficients(force, moment, Ω, V, S, b, c, ρ)

Compute the relevant aerodynamic coefficients given a net force, moment and angular rates ``Ω`` with reference speed ``V``, area ``S``, span ``b``, chord ``c``, and density ``ρ``.
"""
function aerodynamic_coefficients(force, moment, V, S, b, c, ρ)
    q = dynamic_pressure(ρ, V)

    [ force_coefficient(force, q, S)         ;
      moment_coefficient(moment, q, S, b, c) ]
end

reynolds_number(ρ, V, c, μ) = ρ * V * c / μ
reynolds_number(V, c, ν) = V * c / ν

force(CF, q, S) = CF * q * S
moment(CM, q, S, c) = CM * q * S * c
moment(CM, q, S, b, c) = moment.(CM, q, S, [b, c, b])

end


text:
function fizzbuzz(n)
    fiz(x) = "Fizz"
    buzzz(x::String) = string(x, "Buzz")
    buzzz(x) = "Buzz"
    g = identity
    foreach(f -> (iszero∘%)(n, (length∘string)(f)) && (g=f∘g), (fiz, buzzz))
    g(n)
end

(println ∘ fizzbuzz).(1:20)


text:
module TestMmap
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO: Format, FormatEGS, FormatIAEA

PS = [
    EGSParticle{Nothing},
    EGSParticle{Float32},
    IAEAParticle{0,0},
    IAEAParticle{1,0},
    IAEAParticle{0,1},
    IAEAParticle{rand(0:3), rand(0:3)},
]
@testset "PhspVector $P" for P in PS
    mktempdir() do dir
        for len in [1,2,3,4,10,100,rand(1000:2000)]
            ps = [arbitrary(P) for _ in 1:len]
            ext = if Format(first(ps)) isa FormatEGS
                ".egsphsp"
            else
                ".IAEAheader"
            end
            path = joinpath(dir, "phsp" * ext)
            phsp_write(path, ps)

            ps1 = phsp_iterator(collect, path)
            ps2 = PhspVector(path)
            @test all(ps1 .≈ ps)
            @test all(ps2 .≈ ps)

            p1 = first(ps)
            p2 = first(ps2)
            @test p1.w == p2.w
            @test length(ps2) == length(ps)
            @test first(ps) == first(ps2)
            @test last(ps) == last(ps2)
            @test ps1 == ps2

        end
    end
end

@testset "MultiPhspVector $P" for P in PS
    mktempdir() do dir
        # MultiPhspVector
        n = rand(10:100)
        ps = [arbitrary(P) for _ in 1:n]
        paths = String[]
        ext = if Format(first(ps)) isa FormatEGS
            ".egsphsp"
        else
            ".IAEAheader"
        end
        for (i, psi) in enumerate(Iterators.partition(ps, 5))
            path = joinpath(dir, "phsp$i" * ext)
            push!(paths, path)
            phsp_write(path, psi)
        end
        phsp = MultiPhspVector(paths)
        @test phsp isa MultiPhspVector
        @test phsp == ps
    end
end

end#module


text:
# Computes graph Stein discrepancy based on the Langevin operator.
# Args:
#  points - n x d array of sample points
#  weights - n x 1 array of real-valued weights associated with sample points
#    (default: equal weights)
#  gradlogdensity - the gradlogdensity of the target distribution
#  solver - optimization program solver supported by JuMP
function langevin_graph_discrepancy(; points=[],
                                    weights=fill(1/size(points,1), size(points,1)),
                                    gradlogdensity=nothing,
                                    solver=nothing,
                                    kwargs...)
    # make sure solver is defined
    isa(solver, AbstractMathProgSolver) ||
        error("Must specify solver of type AbstractMathProgSolver")
    # Aggregate repeated samples
    sample = SteinDiscrete(points, weights)
    points = sample.support
    weights = sample.weights
    n = length(weights)
    d = size(points,2)
    # Objective coefficients for each g(x_i)
    gradlogdensities = zeros(n, d)
    for i in 1:n
        gradlogdensities[i,:] = gradlogdensity(points[i,:])
    end
    gobjcoefficients = broadcast(*, gradlogdensities, weights)
    # Objective coefficients for each grad g(x_i)
    gradgobjcoefficients = zeros(n, d, d)
    for j = 1:d
        gradgobjcoefficients[:,j,j] = weights
    end
    # solve the optimization problem
    result = affine_graph_discrepancy(sample,
                                      gobjcoefficients,
                                      gradgobjcoefficients,
                                      solver,
                                      "langevin";
                                      kwargs...)
    result
end


text:
# test simple policies on the crying baby problem

using Test

using POMDPModels
using POMDPs

problem = BabyPOMDP(-5, -10, 0.1, 0.8, 0.1, 0.9)

# starve policy
# when the baby is never fed, the reward for starting in the hungry state should be -100
sim = RolloutSimulator()
sim.eps = 0.0001
sim.initialstate = true
ib = EmptyBelief()
policy = Starve()
r = simulate(sim, problem, policy, updater(policy), ib)
@test r ≈ -100.0 atol=0.01

# when the baby is never fed the average reward for starting in the full state should be -47.37
# should take ~5 seconds
# by the Hoeffding inequality there is less than a 1.3% chance that the average will be off by 0.5
n = 100000
r_sum = @parallel (+) for i in 1:n
    sim = RolloutSimulator(MersenneTwister(i),
                           false,
                           0.001,
                           nothing)
    policy = Starve()
    simulate(sim, problem, policy, updater(policy), EmptyBelief())
end
@test r_sum/n ≈ -47.37 atol=0.5

# always feed policy
# when the baby is always fed the reward for starting in the full state should be -50
sim = RolloutSimulator(MersenneTwister(), false, 0.0001, nothing)
policy = AlwaysFeed()
r = simulate(sim, problem, policy, updater(policy), EmptyBelief())
@test r ≈ -50.0 atol=0.01

# when the baby is always fed the reward for starting in the hungry state should be -60
sim = RolloutSimulator(MersenneTwister(), true, 0.0001, nothing)
policy = AlwaysFeed()
r = simulate(sim, problem, policy, updater(policy), EmptyBelief())
@test r ≈ -60.0 atol=0.01

# println("finished easy tests")

# good policy - feed when the last observation was crying - this is *almost* optimal
# from full state, reward should be -17.14
n = 100000
r_sum = @parallel (+) for i in 1:n
    rng = MersenneTwister(i)
    init_state = false
    od = observation(problem, init_state, false, init_state)
    obs = rand(rng, od)
    sim = RolloutSimulator(rng=rng, initialstate=init_state, eps=0.0001)
    policy = FeedWhenCrying()
    simulate(sim, problem, policy, updater(policy), PreviousObservation(obs))
    # println(i)
end
@test r_sum/n ≈ -17.14 atol=0.1

# from hungry state, reward should be -32.11
n = 100000
r_sum = @parallel (+) for i in 1:n
    rng = MersenneTwister(i)
    init_state = true
    od = observation(problem, init_state, false, init_state)
    obs = rand(rng, od)
    sim = RolloutSimulator(rng=rng, initialstate=init_state, eps=0.0001)
    policy = FeedWhenCrying()
    simulate(sim, problem, policy, updater(policy), PreviousObservation(obs))
end
@test r_sum/n ≈ -32.11 atol=0.1


text:
abstract type AbstractVectorFunctionSpace{dim,T,FE<:FiniteElement,MM,N1,N2,N3,fdim} <: DiscreteFunctionSpace{dim,T,FE} end

# VectorFunctionSpace
struct VectorFunctionSpace{dim,T<:Real,FE<:FiniteElement,MM,N1,N2,N3,fdim} <: AbstractVectorFunctionSpace{dim,T,FE,MM,N1,N2,N3,fdim}
    n_dof::Int
    ssp::ScalarFunctionSpace{dim,T,FE,MM,N1,N2,N3,fdim}
end

#Constructor
function VectorFunctionSpace(mesh::PolygonalMesh, felem::FiniteElement{dim,shape,order,1};
    update_face_data = true, quad_degree = order+1) where {dim, shape, order}
    quad_rule = QuadratureRule{dim,shape}(DefaultQuad(), quad_degree)
    f_quad_rule = QuadratureRule{dim-1,shape}(DefaultQuad(), quad_degree)
    fs = _scalar_fs(Float64, mesh, quad_rule, f_quad_rule, felem, update_face_data)
    n_func_basefuncs = getnbasefunctions(felem)
    dof = n_func_basefuncs*dim
    return VectorFunctionSpace(dof, fs)
end

#Data
@inline getnquadpoints(fs::AbstractVectorFunctionSpace) = length(fs.ssp.qr_weights)
@inline getdetJdV(fs::AbstractVectorFunctionSpace, q_point::Int) = getdetJdV(fs.ssp, q_point)
@inline getnbasefunctions(fs::AbstractVectorFunctionSpace) = fs.n_dof
@inline getnlocaldofs(fs::AbstractVectorFunctionSpace) = fs.n_dof
@inline getfiniteelement(fs::AbstractVectorFunctionSpace) = fs.ssp.fe
@inline getmesh(fs::AbstractVectorFunctionSpace) = fs.ssp.mesh

function shape_value(fs::AbstractVectorFunctionSpace{dim,T}, q_point::Int, base_func::Int) where {dim,T}
    @assert 1 <= base_func <= fs.n_dof "invalid base function index: $base_func"
    N_comp = zeros(T, dim)
    n = size(fs.ssp.N,1)
    N_comp[div(base_func,n+1)+1] = fs.ssp.N[mod1(base_func,n),q_point]
    return N_comp
end

function shape_gradient(fs::VectorFunctionSpace{dim,T}, q_point::Int, base_func::Int)  where {dim,T}
    @assert 1 <= base_func <= fs.n_dof "invalid base function index: $base_func"
    dN_comp = zeros(T, dim, dim)
    n = size(fs.ssp.N,1)
    dN_comp[div(base_func,n+1)+1, :] = fs.ssp.dNdξ[mod1(base_func,n), q_point]
    return Tensor{2,dim,T}((dN_comp...,)) ⋅ fs.ssp.Jinv[]
end

function shape_divergence(fs::AbstractVectorFunctionSpace{dim,T}, q_point::Int, base_func::Int)  where {dim,T}
    @assert 1 <= base_func <= fs.n_dof "invalid base function index: $base_func"
    return tr(shape_gradient(fs, q_point, base_func))
end

#Face Data
function shape_value(fs::AbstractVectorFunctionSpace{dim,T}, face::Int, q_point::Int, base_func::Int, orientation::Bool=true) where {dim,T}
    @assert 1 <= base_func <= size(fs.ssp.E,1)*dim "invalid base function index: $base_func"
    N_comp = zeros(T, dim)
    n = size(fs.ssp.E,1)
    q_p = orientation ? q_point : getnfacequadpoints(fs.ssp)-q_point + 1
    N_comp[div(base_func,n+1)+1] = fs.ssp.E[mod1(base_func,n), q_p][face]
    return N_comp
end

@inline get_normal(fs::AbstractVectorFunctionSpace, face::Int) = get_normal(fs.ssp, face)

function reinit!(fs::AbstractVectorFunctionSpace{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    reinit!(fs.ssp, x)
end


text:
using Circuits
using Documenter

makedocs(;
    modules=[Circuits],
    authors="David Gustavsson <david.e.gustavsson@gmail.com> and contributors",
    repo="https://github.com/gustaphe/Circuits.jl/blob/{commit}{path}#L{line}",
    sitename="Circuits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gustaphe.github.io/Circuits.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/gustaphe/Circuits.jl")


text:
@testset "Gabor" begin
    v1 = rand(3)
    v2 = rand(3)
    ell = abs(rand())
    p = abs(rand())
    k = GaborKernel(ell=ell, p=p)
    @test k.ell ≈ ell atol=1e-5
    @test k.p ≈ p atol=1e-5

    k_manual = exp(-sqeuclidean(v1, v2) / (2 * k.ell^2)) * cospi(euclidean(v1, v2) / k.p)
    @test k(v1,v2) ≈ k_manual atol=1e-5

    lhs_manual = transform(SqExponentialKernel(), 1/k.ell)(v1,v2)
    rhs_manual = transform(CosineKernel(), 1/k.p)(v1,v2)
    @test k(v1,v2) ≈ lhs_manual * rhs_manual atol=1e-5

    k = GaborKernel()
    @test k.ell ≈ 1.0 atol=1e-5
    @test k.p ≈ 1.0 atol=1e-5
    @test repr(k) == "Gabor Kernel (ell = 1.0, p = 1.0)"

    test_interface(k, Vector{Float64})

    test_ADs(x -> GaborKernel(ell = x[1], p = x[2]), [ell, p], ADs = [:Zygote])

    # Tests are also failing randomly for ForwardDiff and ReverseDiff but randomly
end