text:
mutable struct FilePathDialog
  dialog::GtkDialog
  signal_insert_id::Culong
  signal_delete_id::Culong
end

mutable struct FilesPanel <: GtkScrolledWindow
    handle::Ptr{Gtk.GObject}
    list::GtkTreeStore
    tree_view::GtkTreeView
    paste_action
    menu
    current_iterator
    path_dialog::FilePathDialog
    main_window::MainWindow
    form_builder::GtkBuilder

    function FilesPanel(main_window::MainWindow)
        sc = GtkScrolledWindow()
        (tv,list,cols) = files_tree_view(["Icon","Name"])
        push!(sc,tv)

        form_builder = Gtk.GtkBuilderLeaf(filename=joinpath(HOMEDIR,"forms/forms.glade"))

        t = new(sc.handle,list,tv,nothing,nothing,nothing,FilePathDialog(form_builder),main_window,form_builder)
        t.menu = filespanel_context_menu_create(t)

        signal_connect(filespanel_treeview_clicked_cb,tv, "button-press-event",
                       Cint, (Ptr{Gtk.GdkEvent},), false,t)
        signal_connect(filespanel_treeview_keypress_cb,tv, "key-press-event",
                       Cint, (Ptr{Gtk.GdkEvent},), false,t)
        signal_connect(filespanel_treeview_row_expanded_cb,tv, "test-expand-row",
                       Cint, (Ptr{Gtk.TreeIter},Ptr{Gtk.TreePath}))

        Gtk.gobject_move_ref(t,sc)
    end
end

global const type_folder = 1
global const type_file = 2
global const type_placeholder = 3 #FIXME: add in GtkExtensions ?

function files_tree_view(rownames)
    n  = length(rownames)
    t = (Gtk.GdkPixbuf,AbstractString, AbstractString, Bool, Int)
    list = GtkTreeStore(t...)
    tv = GtkTreeView(GtkTreeModel(list))
    set_gtk_property!(tv,"level-indentation",4)
    cols = GtkTreeViewColumn[]

    r1 = GtkCellRendererPixbuf()
    c1 = GtkTreeViewColumn(rownames[1], r1, Dict([("pixbuf",0)]))
    Gtk.G_.sort_column_id(c1,0)
    push!(cols,c1)
    #Gtk.G_.max_width(c1,Int(200/n))
    push!(tv,c1)

    r2 = GtkCellRendererText()
    c2 = GtkTreeViewColumn(rownames[2], r2, Dict([("text",1)]))
    Gtk.G_.sort_column_id(c2,1)
    push!(cols,c2)
    #Gtk.G_.max_width(c2,Int(200/n))
    push!(tv,c2)
    return (tv,list,cols)
end

function FilePathDialog(form_builder)
    w = GAccessor.object(form_builder,"DialogCreateFile")
    Gtk.GAccessor.transient_for(w,main_window)#TODO global
    Gtk.GAccessor.modal(w,true)
    btn_create_file = GAccessor.object(form_builder,"btnCreateFile")
    signal_connect(close_file_path_dialog,btn_create_file, "clicked", Nothing, (), false,w)
    return FilePathDialog(w,0,0)

end
 function close_file_path_dialog(btn::Ptr,dialog)
    response(dialog,Gtk.GConstants.GtkResponseType.ACCEPT)
    return nothing
end

function copy_up(model::GtkTreeModel,iter::Gtk.GtkTreeIter)
  parent_iterator = Gtk.mutable(Gtk.GtkTreeIter)
  ret = Gtk.iter_parent(model,parent_iterator, iter)
  if (!ret)
    return nothing
  else
    return parent_iterator[]
  end
end

function filespanel_context_menu_create(t::FilesPanel)
    menu = GtkMenu() |>
    (changeDirectoryItem = GtkMenuItem("Change Directory")) |>
    (addToPathItem = GtkMenuItem("Add to Path")) |>
    GtkSeparatorMenuItem() |>
    (newFileItem = GtkMenuItem("New File")) |>
    (newFolderItem = GtkMenuItem("New Folder")) |>
    GtkSeparatorMenuItem() |>
    (deleteItem = GtkMenuItem("Delete")) |>
    (renameItem = GtkMenuItem("Rename")) |>
    (copyItem = GtkMenuItem("Copy")) |>
    (cutItem = GtkMenuItem("Cut")) |>
    (pasteItem = GtkMenuItem("Paste")) |>
    GtkSeparatorMenuItem() |>
    (copyFullPathItem = GtkMenuItem("Copy Full Path"))

    #FIXME disable until the dialog bugs are fixed
    set_gtk_property!(newFileItem,:sensitive,false)
    set_gtk_property!(newFolderItem,:sensitive,false)
    set_gtk_property!(deleteItem,:sensitive,false)
    set_gtk_property!(renameItem,:sensitive,false)
    set_gtk_property!(copyItem,:sensitive,false)
    set_gtk_property!(pasteItem,:sensitive,false)
    set_gtk_property!(cutItem,:sensitive,false)

    signal_connect(filespanel_changeDirectoryItem_activate_cb,
    changeDirectoryItem, "activate",Nothing, (),false,t)
    signal_connect(filespanel_addToPathItem_activate_cb, addToPathItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_newFileItem_activate_cb, newFileItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_newFolderItem_activate_cb, newFolderItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_deleteItem_activate_cb, deleteItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_renameItem_activate_cb, renameItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_copyItem_activate_cb, copyItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_cutItem_activate_cb, cutItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_pasteItem_activate_cb, pasteItem,
    "activate",Nothing, (),false,t)
    signal_connect(filespanel_copyFullPathItem_activate_cb, copyFullPathItem,
    "activate",Nothing, (),false,t)
    return menu
end

function get_sorted_files(path)
    sort_paths = (x,y)->
    begin
        try
            if isdir(x) && isdir(y)
                return x < y
            elseif isdir(x)
                return true
            elseif isdir(y)
                return false
            else
                return x < y
            end
        catch err
            return true
        end
    end
    sort(readdir(path),lt=sort_paths, by=(x)->return joinpath(path,x))
end

function create_treestore_file_item(path::AbstractString, filename::AbstractString)
    pixbuf = Gtk.icon_theme_load_icon_for_scale(Gtk.icon_theme_get_default(),"code",16,1,0)
    return (pixbuf,filename, joinpath(path,filename), true, type_file)
end
function create_treestore_placeholder_item()
    pixbuf = Gtk.icon_theme_load_icon_for_scale(Gtk.icon_theme_get_default(),"code",24,1,0)
    return (pixbuf,"", "", true, type_placeholder)
end
function add_placeholder(list::GtkTreeStore, parent=nothing)
    return push!(list,create_treestore_placeholder_item(),parent)
end
function add_file(list::GtkTreeStore,path::AbstractString,filename::AbstractString, parent=nothing)
    return push!(list,create_treestore_file_item(path,filename),parent)
end
function insert_file(list::GtkTreeStore,path::AbstractString,filename::AbstractString, iter)
    return insert!(list,iter,create_treestore_file_item(path,filename), how=:sibling, where=:before)
end
function create_treestore_folder_item(path::AbstractString)
  pixbuf = Gtk.icon_theme_load_icon_for_scale(Gtk.icon_theme_get_default(),"folder",24,1,0)
  return (pixbuf,basename(path),path,false,type_folder )
end
function add_folder(list::GtkTreeStore,path::AbstractString, parent=nothing)
    return push!(list,create_treestore_folder_item(path),parent)
end
function insert_folder(list::GtkTreeStore,path::AbstractString, iter)
    return insert!(list,iter,create_treestore_folder_item(path), how=:sibling, where=:before)
end

function populate_folder(w::GtkTreeStore,folder::GtkTreeIter)
    w[folder,4] = true #mark folderr as loaded
    path        = w[folder,3]
    n           = get_sorted_files(path)
    for el in n
        try #some operations are not allowed
            full_path = joinpath(path,string(el))
            if isdir(full_path)
                child = add_folder(w,full_path,folder)
                add_placeholder(w,child)
            else
                file_parts = splitext(el)
                if  (file_parts[2]==".jl" || file_parts[2]==".md")#FIXME need something more general
                    add_file(w,path,el,folder)
                end
            end
        catch err
            @warn "error in populate_folder: " err
        end
    end
end

function where_insert_folder(w::GtkTreeStore,path::AbstractString, parent)
    iter = Gtk.mutable(Gtk.GtkTreeIter)
    iter_valid = Gtk.iter_nth_child(Gtk.GtkTreeModel(w),iter,parent,1)
    while iter_valid && (w[iter,5] == type_folder && w[iter,2] < basename(path))
        iter_valid = Gtk.get_iter_next(Gtk.GtkTreeModel(w), iter)
    end
    if iter_valid
        return iter
    else
        return nothing
    end
end

function where_insert_file(w::GtkTreeStore,p::AbstractString, parent)
  iter = Gtk.mutable(Gtk.GtkTreeIter)

  iter_valid = Gtk.iter_nth_child(Gtk.GtkTreeModel(w),iter,parent,1)
  #Skip folders
  while iter_valid && w[iter,5] == type_folder
    iter_valid = Gtk.get_iter_next(Gtk.GtkTreeModel(w), iter)
  end
  while iter_valid &&  w[iter,2] < basename(p)
      iter_valid = Gtk.get_iter_next(Gtk.GtkTreeModel(w), iter)
  end
  if iter_valid
      return iter
  else
      return nothing
  end
end

function update!(w::GtkTreeStore, p::AbstractString, parent=nothing)
    if isdir(p)
        iter = where_insert_folder(w,p,parent)
        if iter!=nothing
            folder = insert_folder(w,p,iter)
        else
            folder = add_folder(w,p,parent)
        end
        populate_folder(w,folder)
    else
      iter = where_insert_file(w,p,parent)
      if iter!=nothing
          insert_file(w,dirname(p),basename(p),iter)
      else
          add_file(w,dirname(p),basename(p),parent)
      end
    end
end

function update!(w::FilesPanel, path = pwd(current_console(w.main_window)))
    empty!(w.list)
    update!(w.list,path)
    expand_root(w.tree_view)
end

#### FILES PANEL

function get_selected_path(treeview::GtkTreeView,list::GtkTreeStore)
    v = selected(treeview, list)
    if v != nothing
        return v[3]
    else
        return nothing
    end
end

function get_selected_file(treeview::GtkTreeView,list::GtkTreeStore)
    path = get_selected_path(treeview,list)
    if isfile(path)
        return path
    else
        return nothing
    end
end

function file_tab_index(filename,editor::Editor)
    for i=1:length(editor)
        t = get_tab(editor,i)
        if typeof(t) == EditorTab && t.filename == filename #in case we want to have something else in the editor
            return i
        end
    end
    return -1
end

function open_file(treeview::GtkTreeView,list::GtkTreeStore, editor::Editor)
    file = get_selected_file(treeview,list)
    file_idx = file_tab_index(file,editor)
    if file != nothing && file_idx==-1
        open_in_new_tab(file,editor)
    else
        index(editor,file_idx)
    end
end

function nearest_folder_from_current_iterator(filespanel)
    if filespanel.list[filespanel.current_iterator,5] == type_file
        #Get parent folder
        iterator = copy_up(GtkTreeModel(filespanel.list),filespanel.current_iterator)
    else
        iterator = filespanel.current_iterator
    end
    return iterator
end

#=File path menu =#
function path_dialog_create_file(filespanel)
    te_filename = GAccessor.object(filespanel.form_builder,"filename")
    #TODO check overwrite
    iterator = nearest_folder_from_current_iterator(filespanel)
    filename = get_gtk_property(te_filename, :text, AbstractString)
    touch(filename)
    update!(filespanel.list, filename, iterator)
    open_in_new_tab(filename, _editor(filespanel.main_window) )

    return nothing
end

function path_dialog_create_directory(filespanel)
    te_filename = GAccessor.object(filespanel.form_builder,"filename")
    iterator = nearest_folder_from_current_iterator(filespanel)
    #TODO check overwrite
    filename = get_gtk_property(te_filename, :text, AbstractString)
    mkdir(filename)

    update!(filespanel.list,filename,iterator)

    return nothing
end

function path_dialog_rename_file(filespanel)
    te_filename = GAccessor.object(filespanel.form_builder,"filename")
    #TODO check overwrite

    current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
    filename = get_gtk_property(te_filename, :text, AbstractString)
    mv(current_path,filename)

    update!(filespanel.list, filename,copy_up(GtkTreeModel(filespanel.list),filespanel.current_iterator))
    delete!(filespanel.list, filespanel.current_iterator)

    return nothing
end

@guarded function path_dialog_filename_inserted_text(text_entry_buffer_ptr::Ptr,
                                            cursor_pos,
                                            new_text::Cstring,
                                            n_chars,
                                            data)
    path              = data[1]
    delete_signal_id  = data[2]
    text_entry_buffer = convert(GtkEntryBuffer, text_entry_buffer_ptr)
    if (cursor_pos < length(path))
        Gtk.signal_handler_block(text_entry_buffer, delete_signal_id)
        delete_text(text_entry_buffer,cursor_pos,n_chars)
        Gtk.signal_handler_unblock(text_entry_buffer, delete_signal_id)
    end
    return nothing
end

@guarded function path_dialog_filename_deleted_text(text_entry_buffer_ptr::Ptr, cursor_pos,n_chars,data)
    path = data[1]
    insert_signal_id = data[2]
    text_entry_buffer = convert(GtkEntryBuffer, text_entry_buffer_ptr)
    if (cursor_pos < length(path))
        Gtk.signal_handler_block(text_entry_buffer, insert_signal_id)
        insert_text(text_entry_buffer,
                    cursor_pos,
                    path[cursor_pos+1:min(end,cursor_pos+n_chars+1)],
                    n_chars)
        Gtk.signal_handler_unblock(text_entry_buffer, insert_signal_id)
    end
    return nothing
end

#FIXME return mutable struct ?
@guarded function configure_text_entry_fixed_content(te, fixed, nonfixed="")
    set_gtk_property!(te, :text,"");
    set_gtk_property!(te, :text,string(fixed,nonfixed));
    te = buffer(te)
    id_signal_insert = Culong(0)
    id_signal_delete = Culong(0)
    id_signal_insert = signal_connect(path_dialog_filename_inserted_text,
        te,
        "inserted-text",
        Nothing,
        (Cuint,Cstring,Cuint),false,(fixed,id_signal_delete))
    id_signal_delete = signal_connect(path_dialog_filename_deleted_text,
        te,
        "deleted-text",
        Nothing,
        (Cuint,Cuint),false,(fixed,id_signal_insert))
    return (id_signal_insert, id_signal_delete)
end

function configure(dialog::FilePathDialog,
                   filespanel::FilesPanel,
                   path::AbstractString,
                   filename::AbstractString="")
    #Disconnect the Text Entry
    te_filename = GAccessor.object(filespanel.form_builder,"filename")

    te = buffer(te_filename)
    if (dialog.signal_insert_id > 0)
        signal_handler_disconnect(te,dialog.signal_insert_id)
    end
    if (dialog.signal_delete_id > 0)
        signal_handler_disconnect(te,dialog.signal_delete_id)
    end

    if (!isdirpath(path))
	    path = joinpath(path,"")
    end
    (dialog.signal_insert_id, dialog.signal_delete_id) = configure_text_entry_fixed_content(te_filename,path,filename)

end

function file_path_dialog_set_button_caption(filespanel, caption::AbstractString)
    btn_create_file = GAccessor.object(filespanel.form_builder,"btnCreateFile")
    set_gtk_property!(btn_create_file,:label,caption)
end

@guarded (PROPAGATE) function filespanel_treeview_row_expanded_cb(treeviewptr::Ptr,
                                             iterptr::Ptr{Gtk.TreeIter},
                                             path::Ptr{Gtk.TreePath},
                                             data)
    treeview        = convert(GtkTreeView, treeviewptr)
    iter            = unsafe_load(iterptr)
    tree_view_model = Gtk.model(treeview)
    if (tree_view_model[iter,5]==type_folder) && (!tree_view_model[iter,4])
        child_iter =Gtk.mutable(GtkTreeIter)
        if Gtk.iter_nth_child(GtkTreeModel(tree_view_model),child_iter,iter,1)
            delete!(tree_view_model, child_iter[])
        end
        populate_folder(tree_view_model,iter)
    end
    return PROPAGATE
end

@guarded (PROPAGATE) function filespanel_treeview_clicked_cb(widgetptr::Ptr, eventptr::Ptr, filespanel)

    treeview = convert(GtkTreeView, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    list = filespanel.list
    menu = filespanel.menu
    mod = get_default_mod_mask()

    macbook_right_click = Int(event.button) == 1 && (event.state & mod) == GdkModifierType.CONTROL

    if event.button == 3 || macbook_right_click
        (ret,current_path) = Gtk.path_at_pos(treeview,round(Int,event.x),round(Int,event.y));
        if ret
            (ret,filespanel.current_iterator) = Gtk.iter(
                                                      Gtk.GtkTreeModel(filespanel.list),
                                                      current_path)
            showall(menu)
            popup(menu,event)
        end
    else
        if event.event_type == Gtk.GdkEventType.DOUBLE_BUTTON_PRESS
            open_file(treeview,list,_editor(filespanel.main_window))
        end
    end
    return PROPAGATE
end

@guarded (PROPAGATE) function filespanel_treeview_keypress_cb(widgetptr::Ptr, eventptr::Ptr, filespanel)
    treeview = convert(GtkTreeView, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    list = filespanel.list
    if event.keyval == GdkKeySyms.Return
        open_file(treeview,list,_editor(filespanel.main_window))
    end
    return PROPAGATE
end

@guarded function filespanel_newFileItem_activate_cb(widgetptr::Ptr,filespanel)

    if (filespanel.current_iterator != nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        if isfile(current_path)
            current_path = dirname(current_path)
        end
        configure(filespanel.path_dialog ,
                  filespanel,
                  current_path )

        file_path_dialog_set_button_caption(filespanel,"+")

        ret = run(filespanel.path_dialog.dialog)
        visible(filespanel.path_dialog.dialog,false)
        if ret == Gtk.GConstants.GtkResponseType.ACCEPT
            path_dialog_create_file(filespanel)
        end

    end
    return nothing
end

@guarded function filespanel_deleteItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        delete!(filespanel.list,filespanel.current_iterator)
        rm(current_path,recursive=true)
    end
    return nothing
end

@guarded function filespanel_renameItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        base_path = dirname(current_path)
        resource  = basename(current_path)
        #TODO check overwrite
        configure(filespanel.path_dialog ,
                  filespanel,
                  base_path,
                  resource )

        file_path_dialog_set_button_caption(filespanel,"Rename it")
        ret = run(filespanel.path_dialog.dialog)
        visible(filespanel.path_dialog.dialog,false)
        if ret == Gtk.GConstants.GtkResponseType.ACCEPT
            path_dialog_rename_file(filespanel)
        end
    end
    return nothing
end

@guarded function filespanel_changeDirectoryItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        if isfile(current_path)
            current_path = dirname(current_path)
        end
        try
            cd(current_path)
            on_path_change(filespanel.main_window, current_path)
        catch err
        end
    end
    return nothing
end

@guarded function filespanel_addToPathItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        if isfile(current_path)
            current_path = dirname(current_path)
        end
        push!(LOAD_PATH,current_path)
        try
            cd(current_path)
            on_path_change(filespanel.main_window, current_path)
        catch err
        end
    end
    return nothing
end

@guarded function filespanel_newFolderItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        if isfile(current_path)
            current_path = dirname(current_path)
        end
        configure(filespanel.path_dialog ,
                  filespanel,
                  current_path )

        file_path_dialog_set_button_caption(filespanel,"+")
        ret = run(filespanel.path_dialog.dialog)
        visible(filespanel.path_dialog.dialog,false)
        if ret ==Gtk.GConstants.GtkResponseType.ACCEPT
            path_dialog_create_directory(filespanel)
        end

    end
    return nothing
end

@guarded function filespanel_copyItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        filespanel.paste_action = ("copy", filespanel.current_iterator)
    end
    return nothing
end

@guarded function filespanel_cutItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        filespanel.paste_action = ("cut", filespanel.current_iterator)
    end
    return nothing
end

@guarded nothing function filespanel_pasteItem_activate_cb(widgetptr::Ptr,filespanel)
    #TODO check overwrite
    if (filespanel.paste_action != nothing) && (filespanel.current_iterator!=nothing)
        orig_path    = filespanel.list[filespanel.paste_action[2],3]
        iterator     = nearest_folder_from_current_iterator(filespanel)
        current_path = filespanel.list[iterator,3]
        destination  = joinpath(current_path,basename(orig_path))

        if (filespanel.paste_action[1]=="copy")
            cp(orig_path,destination)
        else
            mv(orig_path,destination)
            delete!(filespanel.list, filespanel.paste_action[2])
        end
        update!(filespanel.list, destination,iterator)

    end
    return nothing
end

@guarded function filespanel_copyFullPathItem_activate_cb(widgetptr::Ptr,filespanel)
    if (filespanel.current_iterator!=nothing)
        current_path =  Gtk.getindex(filespanel.list,filespanel.current_iterator,3)
        try
            clipboard(current_path)
        catch e
            println(e)
        end
    end
    return nothing
end

#this is called from MainWindow
function on_path_change(filespanel::FilesPanel, path)
    update!(filespanel, path)
end

# Callback called from Gtk main loop to update the panel
# this is a bad idea because it updates the panel while we try to interact with it
function filespanel_update_cb(user_data)
    filespanel = unsafe_pointer_to_objref(user_data)::FilesPanel
    update!(filespanel)
    return Cint(true)
end


text:
## =========================== ## ## =========================== ##
## 	   Utility functions for CollectiveDynamics statistics		 ##
##	   Martin Zumaya Hernandez 						             ##
##     21 / 03 / 2017									         ##
## =========================== ## ## =========================== ##

### ============== ### ============== ### ============== ###
##                   MAKE DIR FROM PATH                   ##
### ============== ### ============== ### ============== ###
"""
    make_dir_from_path(path)
Try to make dir at path, if already exists doesn't and continues
"""
function make_dir_from_path(path)

    try
        mkdir(path)
    catch error
        println("Folder already exists")
    end

end

### ============== ### ============== ### ============== ###
##           COMPUTE RELATIVE DISTANCES VECTOR            ##
### ============== ### ============== ### ============== ###
"""
    calc_rij_3D_vect(vect)
Compute distance between particles and store them in a 1-D array
# Arguments
* vect -> array of particles' positions
"""
function calc_rij_3D_vect(vect)

    n = div(length(vect),3)

    vec_rij = zeros(Float64, div(n*(n-1),2) )

    k = 1

    for i in 1:3:3n, j in (i+3):3:3n
        vec_rij[k] = norm([vect[i] - vect[j], vect[i+1] - vect[j+1], vect[i+2] - vect[j+2]])
        k += 1
    end

    return vec_rij
end

### ================================== ###

"""
    calc_rij_2D_vect(vect)
Compute distance between particles and store them in a 1-D array
# Arguments
* vect -> array of particles' positions
"""
function calc_rij_2D_vect(vect)

    n = div(length(vect), 2)

    vec_rij = zeros(Float64, div(n*(n-1),2) )

    k = 1

    for i in 1:2:2n, j in (i+2):2:2n
        vec_rij[k] = norm([vect[i] - vect[j], vect[i+1] - vect[j+1]])
        k += 1
    end

    return vec_rij
end

### ============== ### ============== ### ============== ###
##  COMPUTE POSITIONS IN CENTER OF MASS REFERENCE FRAME   ##
### ============== ### ============== ### ============== ###
"""
    calc_vect_3D_cm(vals)
Compute positions in center of mass referenece frame
"""
function calc_vect_3D_cm(vals)

    N = size(vals, 1)

    x_cm = mean(vals[1:3:N, :], 1)
    y_cm = mean(vals[2:3:N, :], 1)
    z_cm = mean(vals[3:3:N, :], 1)

    for i in 1:size(vals,2)
        vals[1:3:N, i] -= x_cm[i]
        vals[2:3:N, i] -= y_cm[i]
        vals[3:3:N, i] -= z_cm[i]
    end
end

### ================================== ###
"""
    calc_vect_2D_cm(vals)
Compute positions in center of mass referenece frame
"""
function calc_vect_2D_cm(vals)

    N = size(vals, 1)

    x_cm = mean(vals[1:2:N, :], 1)
    y_cm = mean(vals[2:2:N, :], 1)

    for i in 1:size(vals,2)
        vals[1:2:N, i] -= x_cm[i]
        vals[2:2:N, i] -= y_cm[i]
    end
end

### ================================== ###

### ============== ### ============== ### ============== ###
##                 COMPUTE SAMPLING TIMES                 ##
### ============== ### ============== ### ============== ###
"""
    get_times(T)
Compute sampling times for 10^T time steps
"""
function get_times(T)

    times = [convert(Int, exp10(i)) for i in 0:T]

    tau = Int64[]

    for i in 1:(length(times) - 1)

        if i > 1

            for t in (times[i]+1):times[i+1]

                if t % times[i] == 0 || t % times[i-1] == 0
                    push!(tau, t)
                end
            end

        else

            for t in (times[i]+1):times[i+1]

                if t % times[i] == 0
                    push!(tau, t)
                end
            end

        end

    end

    return tau
end

### ============== ### ============== ### ============== ###
##             FINITE DIFFERENCES DERIVATIVES             ##
### ============== ### ============== ### ============== ###
function deriv_interpol(vals, time)

    n = length(vals)
    der = zeros(n)

    for i in 1:3:n

        h_1 = vals[i+1] - vals[i]
        h_2 = vals[i+2] - vals[i]

        der[i] = ((2*h_1 + h_2)/(h_1*(h_1+h_2))) * vals[i] + ((h_1 + h_2)/(h_1 * h_2)) * vals[i+1] - ((h_1)/(h_2 * (h_1 + h_2))) * vals[i+2]
        der[i + 1] = -(h_2/(h_1*(h_1+h_2))) * vals[i] - ((h_1 - h_2)/(h_1 * h_2)) * vals[i+1] - (h_1/(h_2 * (h_1 + h_2))) * vals[i+2]
        der[i + 2] = (h_2/(h_1*(h_1+h_2))) * vals[i] - ((h_1 + h_2)/(h_1 * h_2)) * vals[i+1] + ((h_1 + 2*h_2)/(h_2 * (h_1 + h_2))) * vals[i+2]

    end

    return der
end

ext: jl


text:
# Open a URL in a browser
function openurl(url::AbstractString)
    @static if is_apple() run(`open $url`) end
    @static if is_windows() run(`cmd /c start $url`) end
    @static if is_linux()   run(`xdg-open $url`) end
end

#Jupyter Notebook display
import Base.show

asset(url...) = readall(joinpath(dirname(@__FILE__), "..", "assets", "bower_components", url...))

function writemime(io::IO, ::MIME"text/html", v::VegaVisualization)
    divid = "vg" * randstring(3)
    script_contents = scriptstr(v, divid)
    display("text/html", """

    """)
end

import Patchwork: Elem

function patchwork_repr(v::VegaVisualization)
    divid = "vg" * randstring(3)
    script_contents = scriptstr(v, divid)
    Elem(:div, [
        Elem(:div, "") & Dict(:id=>divid, :style=>Dict("min-height"=>"$(v.height + 110)px")),
        Elem(:script, script_contents) & Dict(:type=>"text/javascript")
    ])
end

function show(io::IO, m::MIME"text/html", v::VegaVisualization)
    writemime(io, m, patchwork_repr(v))
end

function scriptstr(v::VegaVisualization, divid)
    spec = JSON.json(tojs(v))
    return """
        require.config({
          paths: {
            d3: ["https://vega.github.io/vega-editor/vendor/d3.min", "/files/vegajl/assets/bower_components/d3/d3.min"],
            vega: ["https://vega.github.io/vega/vega.min", "/files/vegajl/assets/bower_components/vega/vega.min"],
            cloud: ["https://vega.github.io/vega-editor/vendor/d3.layout.cloud", "/files/vegajl/assets/bower_components/d3-cloud/build/d3.layout.cloud"],
            topojson: ["https://vega.github.io/vega-editor/vendor/topojson", "/files/vegajl/assets/bower_components/topojson/topojson"]
          }
        });

        require(["d3", "topojson", "cloud"], function(d3, topojson, cloud){

            window.d3 = d3;
            window.topojson = topojson;
            window.d3.layout.cloud = cloud;

                  require(["vega"], function(vg) {

                  vg.parse.spec($spec, function(chart) { chart({el:\"#$divid\"}).update(); });

                  window.setTimeout(function() {
                    var pnglink = document.getElementById(\"$divid\").getElementsByTagName(\"canvas\")[0].toDataURL(\"image/png\")
                    document.getElementById(\"$divid\").insertAdjacentHTML('beforeend', '
Save as PNG')

                  }, 20);

              }); //vega require end

        }); //d3 require end
    """
end

function body(v::VegaVisualization)
    divid = "vg" * randstring(3)

    "

    "
end

function genhtml(v::VegaVisualization; title="Vega.jl Visualization")
    d3 = asset("d3","d3.min.js")
    topojson = asset("topojson","topojson.js")
    cloudlayout = asset("d3-cloud", "build", "d3.layout.cloud.js")
    vega = asset("vega", "vega.js")

    "

        <title>$title</title>

      $(body(v))

   "
end

#Vega Scaffold: https://github.com/vega/vega/wiki/Runtime
function writehtml(io::IO, v::VegaVisualization; title="Vega.jl Visualization")
    println(io, genhtml(v, title=title))
end

function Base.show(io::IO, v::VegaVisualization)

    if displayable("text/html")
        v
    elseif BLINK_INSTALLED
        blink_show(v)
    else
        # create a temporary file
        tmppath = string(tempname(), ".vega.html")
        io = open(tmppath, "w")
        writehtml(io, v)
        close(io)

        # Open the browser
        openurl(tmppath)

    end

    return
end

max_stars_repo_name: jpthor/vega


text:
# Generate random circuits and count the logic functions that are computed.
# Similar to results of Raman and Wagner 2011 to generate their Figure 2b.
# Run after including ../src/CGP.jl or uncomment the following line
#include("../src/CGP.jl")
# The important functions used are in ../src/Execute.jl

# Note that the parameters are set below
# The default functions are defined in ../src/Func.jl
# The types  MyInt and MyFunc are defined in CGP.jl
# Note that the printed output would be long, output is truncated in the functions in Execute.jl.

#using Base.Test
using Printf
using Main.CGP

# The following lines allow this file to included from CGP.jl, CGP.jl/test, or CGP.jl/src.
cwd = pwd()
if cwd[end-5:end] == "CGP.jl"
  csvfile = "$cwd/test/data/test_inf_alleles.csv"
  outfile = "$cwd/test/data/test_generate_random_functions.txt"
else
  csvfile = "../test/data/test_inf_alleles.csv"
  outfile = "../test/data/test_generate_random_functions.txt"
end
f = open(outfile,"w")
numinputs = 2
numoutputs = 2
numinteriors = 4
numlevelsback = numinteriors
funcs = Main.CGP.default_funcs(numinputs)
context = Main.CGP.construct_context(numinputs)
p = Parameters( numinputs=numinputs, numoutputs=numoutputs, numinteriors=numinteriors, numlevelsback=numlevelsback )

num_iterations = 50   # Use a value <= 50 to see full untruncated output
#num_iterations = 10000
println(f,"num_iterations: ",num_iterations)
print_parameters(f,p)
counts = create_count_function_hash( p.numinputs, p.numoutputs )
for _ = 1:num_iterations
    c = random_chromosome(p, funcs)
    #print_chromosome( c )
    outputs = execute_chromosome(c, context)
    if num_iterations <= 50
      print(f,"outputs: ",outputs)
      if length(outputs) > 1
        println(f,"   hamming: ",hamming(outputs[1],outputs[2]))
      else
        println(f)
      end
    end
    concat_outputs = concatenate_outputs(p.numinputs, outputs)
    #Printf.@printf(f,"concat_outputs: 0x%2x\n",concat_outputs)
    increment_count_function_hash(concat_outputs, counts )
    outputs
end
print_count_function_hash(f,counts)
print_count_function_summary(f,counts)
close(f)


text:
"""
function chroma(stft::Matrix{Complex{Cdouble}},
	fs::Cdouble, wsize::Cuint, γ::Cdouble, ϵ::Cdouble
)::Matrix{Cdouble}

Compute the chroma features of a signal from its Short Time Fourier Transform.

The transform must be provided as a matrix in which each column represents a
frame and with time evolving along the second dimension. The return will be a
matrix of size 12xN where N is the number of frames of the STFT and each column
is a vector of normalized chroma features.

**Parameters:**

`stft`: The Short Time Fourier Transform of the signal from which the features
must be computed.

`fs`: The sampling rate of the signal.

`wsize`: The number of samples of the window used in the STFT.

`γ`: The compression factor to be used in the chroma features extraction.

`ϵ`: The smallest norm that a chroma vector needs to have to not be replaced by
a stub one.
"""
function chroma(
	stft::Matrix{Complex{Cdouble}},
	fs::Cdouble,
	wsize::Cuint,
	γ::Cdouble,
	ϵ::Cdouble
)::Matrix{Cdouble}

	# The number of chroma vectors is defined by the frames in the STFT
	nofframes = size(stft)[2]
	pitchspectrum = Matrix{Cdouble}(undef, 128, nofframes)

	# Relation of a pitch with a frequency defined by the MIDI standard
	pitch2freq(p) = round(Cint, (wsize/fs)*2^((p - 69)/12)*440)

	# The bound frequencies of each one of the pitch classes
	bounds = [1 + pitch2freq(p - 0.5) for p in 0:128]

	for n in 1:nofframes

		# The pitch spectrum is a log spaced spectrum where a component is
		# proportional to the sum of the power of all stft components inside the
		# bounds of a pitch class
		pitchspectrum[:, n] =
			[sum(abs2.(stft[bounds[i]:bounds[i + 1],n])) for i in 1:128]

	end

	chroma = Matrix{Cdouble}(undef, 12, nofframes)

	for n in 1:nofframes

		# A chroma vector is given by the sum of the pitch classes of the same
		# note in all octaves
		chroma[:, n] =
			[sum(pitchspectrum[(class % 12 + 1):12:end, n]) for class in 0:11]

	end

	# Compresses the chroma vectors
	chroma = log10.(1 .+ γ .* chroma)

	stubFeature = ones(12)./√12

	# Evaluates the norm of each chroma vector...
	for i in 1:nofframes

		# ...and if it's smaller than the ϵ parameter normalizes the vector
		if ((norm = √(sum(chroma[:, i].^2))) > ϵ)
			chroma[:, i] ./= norm
		# Otherwise replaces the vector with a stub one
		else
			chroma[:, i] = stubFeature
		end

	end

	return chroma

end

max_stars_repo_name: gabriel-ss/interpretation-switcher

max_issues_repo_name: gabriel-ss/interpretation-switcher

max_forks_repo_name: gabriel-ss/interpretation-switcher


text:
using BAT
using DensityInterface
using ValueShapes
using Distributions
using Plots
using StatsBase
using Statistics
using SpecialFunctions
using ArraysOfArrays
using TypedTables
using CSV
import Cuba
using AHMI

function log_pdf_poisson(λ::T, k::U) where {T<:Real,U<:Real}
    R = float(promote_type(T,U))
    if λ >= 0 && k >= 0 && isinteger(k)
        a = iszero(k) ? R(k) : R(log(λ) * k)
        lg = logabsgamma(k + 1)[1]
        R(a - λ - lg)
    else
        R(-Inf)
    end
end

function bkg_signal_pdf(
    ν_B::Real,
    ν_S::Real,
    λ::Real,
    S_μ::Real,
    S_σ::Real,
    E::Real)

    res = (ν_B * pdf(Exponential(λ), E) + ν_S * pdf(Normal(S_μ, S_σ), E)) / (ν_B + ν_S)
    return res
end

function log_bkg_signal_pdf(
    ν_B::Real,
    ν_S::Real,
    λ::Real,
    S_μ::Real,
    S_σ::Real,
    E::Real)

    return log(bkg_signal_pdf(ν_B, ν_S, λ, S_μ, S_σ, E))
end

function μ_log_normal(
    γ::Real,
    σ::Real)

    return log(γ) - σ^2 / 2.0
end

function fit_function_sum_all(p, E)

    ν_B = summary_dataset_table.exposure .* p.B
    ν_S = summary_dataset_table.exposure .* summary_dataset_table.efficiency .* p.S

    return sum(summary_dataset_table.exposure .* bkg_signal_pdf.(ν_B, ν_S, p.λ, p.S_μ, p.S_σ, E)) / sum(summary_dataset_table.exposure)

end

# We can define a likelihood via logfuncdensity
function make_likelihood_bkg(summary_dataset_table, sample_table)
    logfuncdensity(function(p)
        ll_a = sum(log_pdf_poisson.(p.B .* summary_dataset_table.exposure, summary_dataset_table.event_counts))
        ll_b = sum(-log(p.λ) .- sample_table.E ./ p.λ)
        return ll_a + ll_b
    end)
end

# We can also define a likelihood as a struct that supports `densitykind` and `logdensityof`:
struct SignalBkgLikelihood{DS,ST}
    summary_dataset_table::DS
    sample_table::ST
end

@inline DensityInterface.DensityKind(::SignalBkgLikelihood) = IsDensity()

function DensityInterface.logdensityof(likelihood::SignalBkgLikelihood, p)
    summary_dataset_table =  likelihood.summary_dataset_table
    sample_table =  likelihood.sample_table

    ν_B = summary_dataset_table.exposure .* p.B
    ν_S = summary_dataset_table.exposure .* summary_dataset_table.efficiency .* p.S
    help_B = ν_B[sample_table.dataset]
    help_S = ν_S[sample_table.dataset]

    ll_a = sum(log_pdf_poisson.(ν_B .+ ν_S, summary_dataset_table.event_counts))
    ll_b = sum(log_bkg_signal_pdf.(help_B, help_S, p.λ, p.S_μ, p.S_σ, sample_table.E))

    ll_a + ll_b
end

N_datasets = 5
E_min = 0.0
E_max = 200.0
signal_events = 3
true_λ = 50.0

ΔE = E_max - E_min

summary_dataset_table = CSV.read("summary_dataset_table.csv", Table)
sample_table = CSV.read("sample_table.csv", Table)

function make_child_prior(N)
    v -> begin
        return NamedTupleDist(B = fill(LogNormal(μ_log_normal(v.m_B, v.σ_B), v.σ_B), N))
    end
end

parent_prior_bkg = NamedTupleDist(
    σ_B = Uniform(0.1, 1.0),
    m_B = Uniform(10^-10, 1e-1 * ΔE),
    λ = Uniform(10^-10, 100.0)
)

parent_prior_bkg_signal = NamedTupleDist(
    S = Uniform(0.0, 10.0),
    S_μ = 100.0,
    S_σ = 2.0,
    σ_B = Uniform(0.1, 1.0),
    m_B = Uniform(10^-10, 1e-1 * ΔE),
    λ = Uniform(10^-10, 100.0)
)

prior_bkg =  HierarchicalDistribution(make_child_prior(length(summary_dataset_table)), parent_prior_bkg)

prior_bkg_signal = HierarchicalDistribution(make_child_prior(length(summary_dataset_table)), parent_prior_bkg_signal)

posterior_bkg = PosteriorDensity(make_likelihood_bkg(summary_dataset_table, sample_table), prior_bkg)

posterior_bkg_signal = PosteriorDensity(SignalBkgLikelihood(summary_dataset_table, sample_table), prior_bkg_signal)

nchains = 4
nsteps = 10^5

algorithm = MCMCSampling(mcalg = HamiltonianMC(), nchains = nchains, nsteps = nsteps)

samples_bkg = bat_sample(posterior_bkg, algorithm).result

@show evidence_bkg_ahmi = bat_integrate(samples_bkg, AHMIntegration()).result
@show evidence_bkg_cuba = bat_integrate(posterior_bkg, VEGASIntegration(log_density_shift = -maximum(samples_bkg.logd), maxevals = 10^6, rtol = 0.005)).result

samples_bkg_signal = bat_sample(posterior_bkg_signal, algorithm).result

@show evidence_bkg_signal_ahmi = bat_integrate(samples_bkg_signal, AHMIntegration()).result
@show evidence_bkg_signal_cuba = bat_integrate(posterior_bkg_signal, VEGASIntegration(log_density_shift = -maximum(samples_bkg_signal.logd), maxevals = 10^6, rtol = 0.005)).result

@show BF_exponential_ahmi = evidence_bkg_signal_ahmi / evidence_bkg_ahmi
@show BF_exponential_cuba = evidence_bkg_signal_cuba / evidence_bkg_cuba

@show bkg_sig_marginal_modes = bat_marginalmode(samples_bkg_signal).result

p_1 = plot(size=(800,500), layout=(2,2), labelfontsize=12, tickfontsize=10, legendfontsize=7)
#upper left
p_1 = plot!(samples_bkg_signal, :S, subplot=1, label = "Posterior")
p_1 = plot!(parent_prior_bkg_signal, :S, subplot=1, label = "Prior", linecolor = "blue")
#lower right
p_1 = plot!(samples_bkg_signal, :λ, subplot=4, label = "Posterior", legend=false)
p_1 = plot!(parent_prior_bkg_signal, :λ, subplot=4, label = "Prior", linecolor = "blue")
#upper right
p_1 = plot!(samples_bkg_signal, (:S, :λ), subplot=2, st = :histogram, legend=false, colorbar=false)
#lower left
p_1 = plot!(samples_bkg_signal, (:S,:λ), subplot=3, legend=true)

savefig(p_1, "prior_posterior.pdf")
savefig(p_1, "prior_posterior.png")

p_2 = plot(size=(800,500), layout=(2,2), labelfontsize=12, tickfontsize=10, legendfontsize=7)
#upper left
p_2 = plot!(samples_bkg_signal, :m_B, subplot=1, label = "Posterior")
p_2 = plot!(parent_prior_bkg_signal, :m_B, subplot=1, label = "Prior", linecolor = "blue")
#lower right
p_2 = plot!(samples_bkg_signal, :σ_B, subplot=4, label = "Posterior", legend=false)
p_2 = plot!(parent_prior_bkg_signal, :σ_B, subplot=4, label = "Prior", linecolor = "blue")
#upper right
p_2 = plot!(samples_bkg_signal, (:m_B, :σ_B), subplot=2, st = :histogram, legend=false, colorbar=false)
#lower left
p_2 = plot!(samples_bkg_signal, (:m_B,:σ_B), subplot=3, legend=true)

savefig(p_2, "prior_posterior_hierarchical.pdf")
savefig(p_2, "prior_posterior_hierarchical.png")

p_hist = plot(size=(800,500), layout=(1,1), labelfontsize=12, tickfontsize=10, legendfontsize=7)
p_hist = histogram!(sample_table.E, bins = range(0.0, stop=maximum(sample_table.E)+20., length=100), title = "", xlabel = "Energy [keV]", ylabel = "Counts", label = "", box = :on, grid = :off)

savefig(p_hist, "total_hist.pdf")
savefig(p_hist, "total_hist.png")

p_fit_sum = plot(size=(800,500), layout=(1,1), labelfontsize=12, tickfontsize=10, legendfontsize=7)
p_fit_sum = plot!(range(0.0, (maximum(sample_table.E)+20), length=500), fit_function_sum_all, samples_bkg_signal, box = :on, grid = :off, xlabel = "Energy [keV]", ylabel = "Background distribution", legend = :topright)
p_fit_sum = histogram!(twinx(), xticks=([], []), sample_table.E, bins = range(0.0, stop=maximum(sample_table.E)+20., length=100), ylabel = "Counts", label = "Binned data", box = :on, grid = :off, fillalpha = 0.4, linealpha = 0.4, legend = :topleft)

savefig(p_fit_sum, "detector_sum.pdf")
savefig(p_fit_sum, "detector_sum.png")


text:
"""
    A version of E7 for testing linearization
"""
module E7A

using ModelBaseEcon

model = Model()
model.options.substitutions = false

@parameters model begin
    delta = 0.1000000000000000
    p_dlc_ss = 0.0040000000000000
    p_dlinv_ss = 0.0040000000000000
    p_growth = 0.0040000000000000
end

@variables model begin
    dlc; dlinv; dly; lc; linv;
    lk; ly;
end

@shocks model begin
    dlc_shk; dlinv_shk;
end

@autoexogenize model begin
    lc = dlc_shk
    linv = dlinv_shk
end

@equations model begin
    dlc[t] = (1 - 0.2 - 0.2) * p_dlc_ss + 0.2 * dlc[t - 1] + 0.2 * dlc[t + 1] + dlc_shk[t]
    dlinv[t] = (1 - 0.5) * p_dlinv_ss + 0.1 * dlinv[t - 2] + 0.1 * dlinv[t - 1] + 0.1 * dlinv[t + 1] + 0.1 * dlinv[t + 2] + 0.1 * dlinv[t + 3] + dlinv_shk[t]
    lc[t] = lc[t - 1] + dlc[t]
    linv[t] = linv[t - 1] + dlinv[t]
    ly[t] = log(exp(lc[t]) + exp(linv[t]))
    dly[t] = ly[t] - ly[t - 1]
    lk[t] = log((1 - delta) * exp(lk[t - 1]) + exp(linv[t]))
end

@initialize model

@steadystate model linv = lc - 7;
@steadystate model lc = 14;

end


text:
################################################################################
#
#   DEFINITION OF 2D UNITCELLS
#
#   Code is included in subfiles
#
################################################################################

# SQUARE / RECTANGULAR LIKE

# Square lattice
include("unitcells_2d/square.jl")

# Square octagon lattice
include("unitcells_2d/square_octagon.jl")

# Shastry Sutherland lattice
include("unitcells_2d/shastry_sutherland.jl")

# TRIANGULAR LIKE

# Triangular lattice
include("unitcells_2d/triangular.jl")

# Honeycomb lattice
include("unitcells_2d/honeycomb.jl")

# Kagome lattice
include("unitcells_2d/kagome.jl")


text:
#src Generate README.md by running `using Literate; Literate.markdown("docs/src/README.jl", "."; documenter=false, credit=false, execute=true)`

# <img align="right" width="180px" src="https://avatars.githubusercontent.com/u/58918144?s=200&v=4">

# <!-- DO NOT EDIT README.md directly, instead edit docs/README.jl and generate the markdown-->

# # Probabilistic<wbr>Circuits<wbr>.jl

# [![Unit Tests](https://github.com/Juice-jl/ProbabilisticCircuits.jl/workflows/Unit%20Tests/badge.svg)](https://github.com/Juice-jl/ProbabilisticCircuits.jl/actions?query=workflow%3A%22Unit+Tests%22+branch%3Amaster)  [![codecov](https://codecov.io/gh/Juice-jl/ProbabilisticCircuits.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/Juice-jl/ProbabilisticCircuits.jl) [![](https://img.shields.io/badge/docs-stable-green.svg)](https://juice-jl.github.io/ProbabilisticCircuits.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://juice-jl.github.io/ProbabilisticCircuits.jl/dev)

# This package provides functionalities for learning/constructing probabilistic circuits and using them to compute various probabilistic queries. It is part of the [Juice package](https://github.com/Juice-jl) (Julia Circuit Empanada).

# ## Example usage

include("usage.jl")

# ## Testing

# To make sure everything is working correctly, you can run our test suite as follows. The first time you run the tests will trigger a few slow downloads of various test resources.

# ```bash
# julia --color=yes -e 'using Pkg; Pkg.test("ProbabilisticCircuits")'
# ```


text:
module Fourier

export mdft, imdft, fftrange

include("mdft.jl")

end

max_stars_repo_name: JuliaOptics/Fourier

max_issues_repo_name: JuliaOptics/Fourier

max_forks_repo_name: JuliaOptics/Fourier


text:
module TimeScales

using ItemGraphs: ItemGraph, add_edge!, add_vertex!, items

import Dates

export
    TimeScale,
    find_path,
    link_scales!,
    register_scale!

"""
All timescales are subtypes of the abstract type `TimeScale`.
The following timescales are defined:

* [`UT1`](@ref) — Universal Time
* [`TAI`](@ref) — International Atomic Time
* [`TT`](@ref) — Terrestrial Time
* [`TCG`](@ref) — Geocentric Coordinate Time
* [`TCB`](@ref) — Barycentric Coordinate Time
* [`TDB`](@ref) — Barycentric Dynamical Time
"""
abstract type TimeScale end

const NAMES = (
    :UniversalTime,
    :InternationalAtomicTime,
    :TerrestrialTime,
    :GeocentricCoordinateTime,
    :BarycentricCoordinateTime,
    :BarycentricDynamicalTime,
)

const ACRONYMS = (
    :UT1,
    :TAI,
    :TT,
    :TCG,
    :TCB,
    :TDB,
)

for (acronym, name) in zip(ACRONYMS, NAMES)
    acro_str = String(acronym)
    name_str = String(name)
    name_split = join(split(name_str, r"(?=[A-Z])"), " ")
    wiki = replace(name_split, " "=>"_")
    @eval begin
        """
            $($name_str)

        A type representing the $($name_split) ($($acro_str)) time scale.

        # References

        - [Wikipedia](https://en.wikipedia.org/wiki/$($wiki))
        """
        struct $name <: TimeScale end

        """
            $($acro_str)

        The singleton instance of the [`$($name_str)`](@ref) type representing
        the $($name_split) ($($acro_str)) time scale.

        # References

        - [Wikipedia](https://en.wikipedia.org/wiki/$($wiki))
        """
        const $acronym = $name()

        export $name, $acronym

        Base.show(io::IO, ::$name) = print(io, "$($acro_str)")
        tryparse(::Val{Symbol($acro_str)}) = $acronym
    end
end

const SCALES = ItemGraph{TimeScale}()

function register_scale!(s)
    add_vertex!(SCALES, s)
end

function link_scales!(s1, s2; oneway=false)
    add_edge!(SCALES, s1, s2)
    oneway || add_edge!(SCALES, s2, s1)
end

link_scales!(TAI, TT)
link_scales!(TAI, UT1)
link_scales!(TT, TCG)
link_scales!(TT, TDB)
link_scales!(TCB, TDB)

find_path(from, to) = items(SCALES, from, to)

struct NotATimeScale <: TimeScale end

tryparse(s::AbstractString) = tryparse(Val(Symbol(s)))
tryparse(::Any) = nothing

@inline function Dates.tryparsenext(d::Dates.DatePart{'t'}, str, i, len, locale)
    next = Dates.tryparsenext_word(str, i, len, locale, Dates.max_width(d))
    next === nothing && return nothing
    word, i = next
    val = tryparse(word)
    val === nothing && throw(ArgumentError("'$word' is not a recognized time scale."))
    return val, i
end

end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
# # This is a title

# This is a description


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Colors, ColorSchemes, Random
Random.seed!(123)
function custom_palette(colormap=:viridis)
    cmap = get(colorschemes[colormap], LinRange(0,1,256))
    ctmp = "0 '#$(hex(cmap[1]))',"
    for i in 2:256; ctmp = ctmp*"$(i-1) '#$(hex(cmap[i]))'," end;
    "set palette defined("*ctmp[1:end-1]*")"
end
@gp " " :-  # this is just to reset session
@gp :- "set multiplot layout 3,3; set key off" :-
cmaps = [:magma, :viridis, :plasma, :inferno, :berlin, :leonardo,
    :devon, :spring, :ice]
for i in 1:9
    @gp :- i title = "$(cmaps[i]), pt $(i)" "set size square" :-
    @gp(:-, rand(15), rand(15), rand(15),
    "w p pt $(i) ps 3 lc palette", custom_palette(cmaps[i]),
    "set cbtics out nomirror", "set xtics 0,1", "set ytics 0,1",
    "set cbtics 0.0,0.2,1.0", xrange= (-0.1,1.1),
    yrange= (-0.1,1.1), :-)
end
@gp
save(term="pngcairo size 900,800", output="plt_ex12.png") # hide


text:
using DashBootstrapComponents, DashHtmlComponents

fluid_jumbotron = dbc_jumbotron(
    [
        dbc_container(
            [
                html_h1("Fluid jumbotron", className="display-3"),
                html_p(
                    "This jumbotron occupies the entire horizontal " *
                    "space of its parent.",
                    className="lead",
                ),
                html_p(
                    "You will need to embed a fluid container in " * "the jumbotron.",
                    className="lead",
                ),
            ],
            fluid=true,
        ),
    ],
    fluid=true,
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
struct BitsUInt8TupleInt32Int64
    a::UInt8
    b::Tuple{Int32, Int64}
end

struct BitsUInt8TupleInt32TupleInt16UInt16
    a::UInt8
    b::Tuple{Int32, Tuple{Int16, UInt16}}
end

@testset "Bits types with tuple fields" begin
    @test begin
        b = JlrsReflect.reflect([BitsUInt8TupleInt32Int64])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsUInt8TupleInt32Int64] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsUInt8TupleInt32Int64")]
        pub struct BitsUInt8TupleInt32Int64 {
            pub a: u8,
            pub b: ::jlrs::wrappers::inline::tuple::Tuple2<i32, i64>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsUInt8TupleInt32TupleInt16UInt16])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsUInt8TupleInt32TupleInt16UInt16] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsUInt8TupleInt32TupleInt16UInt16")]
        pub struct BitsUInt8TupleInt32TupleInt16UInt16 {
            pub a: u8,
            pub b: ::jlrs::wrappers::inline::tuple::Tuple2<i32, ::jlrs::wrappers::inline::tuple::Tuple2<i16, u16>>,
        }"""
    end
end


text:
using Luxor, Colors

Drawing(640, 640, "pattern.png")
origin()

pt1 = Point(0, 12)
pt2 = Point(25, 12)
pt3 = Point(25, 37)
pt4 = Point(0, 37)
pt5 = Point(9, 0)

function draw(size)
    size < 5 && return
    randomhue()
    poly([pt1, pt2, pt3, pt4], :fill)
    randomhue()
    poly([pt1, pt5, pt2], :fill)
    axes()
    textcentred(@sprintf "%d" size)
    gsave()
    rotate(atan2(4, -3))
    translate(pt5)
    scale(3/5, 3/5)
    draw(size*3/5)
    grestore()
    gsave()
    rotate(-atan2(3, 4))
    translate(pt2)
    scale(4/5, 4/5)
    draw(size*4/5)
    grestore()
end

scale(10, 10)
draw(10)
finish()

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
function metric_curvature(xn, vn, task::TaskGDS{<:TaskMapT{M,N}}, C::Chart{I,N}) where {M,N,I}
    n = dim(N)
    ∂g_∂v = ForwardDiff.jacobian(vn -> reshape(metric_chart(xn, vn, task, C), Size(n^2)), vn)
    Ξ = SMatrix{n,n,eltype(xn)}(0.5*sum(∂g_∂v[(i-1)*n+1:i*n,:].*vn[i] for i in 1:n)')
end

function force_curvature(xn, vn, task::TaskGDS{<:TaskMapT{M,N}}, C::Chart{I,N}) where {M,N,I}
    n = dim(N)
    ∂g_∂x = SMatrix{n^2,n,eltype(xn)}(ForwardDiff.jacobian(xn -> reshape(metric_chart(xn, vn, task, C), n^2), xn))
    ξ1 = reshape(∂g_∂x*vn, Size(n,n))*vn
    ξ2 = ForwardDiff.gradient(xn -> (vn'*metric_chart(xn, vn, task, C)*vn)[1], xn)
    ξ = ξ1 - 0.5*ξ2
end

function single_task_components(xm, vm, task::TaskGDS{<:TaskMapT{M,N}}, CM::Chart{I,M},
        CN::Chart{J,N}) where {M,N,I,J}
    n = dim(N)
    xn = base_task_map_chart(xm, task, CM, CN)
    Jf = base_task_jacobian_chart(xm, task, CM, CN)
    vn = Jf*vm
    Jfdot = base_task_jacobian_chart_dot(xm, vm, task, CM, CN)

    G = metric_chart(xn, vn, task, CN)
    Ξ = metric_curvature(xn, vn, task, CN)
    ξ = force_curvature(xn, vn, task, CN)
    ℱ_pot = potential_force_chart(xn, task, CN)
    B = dissipative_term_chart(xn, vn, task, CN)
    A = Ξ + G

    xnddot_des = SVector{n,eltype(xm)}(pinv(A)*(ℱ_pot - B*vn - ξ))

    JftAJf = Jf'*A*Jf
    JftA_a = Jf'*A*(xnddot_des - Jfdot*vm)

    JftAJf, JftA_a
end

function single_task_acceleration(xm, vm, task::TaskGDS{<:TaskMapT{M,N}}, CM::Chart{I,M},
        CN::Chart{J,N}, robot_coord_rep=ChartRep()) where {M,N,I,J}
    m, n = dim(M), dim(N)
    (xm, vm) = robot_coord_rep == EmbRep() ? emb_to_chart_differential(xm, vm, CM) : (xm, vm)
    CN = isglobal(CN) ? CN : choose_chart_chart(xm, task, CM, CN)
    JftAJf, JftA_a = single_task_components(xm, vm, task, CM, CN)
    σxddot = SMatrix{m,m,eltype(xm)}(pinv(Matrix(JftAJf)))*JftA_a
    if robot_coord_rep == EmbRep()
        xme, vme, σxddot = chart_to_emb_differential(xm, vm, σxddot, CM)
    end
    σxddot, CN
end

function multiple_task_acceleration(xm, vm, tasks::TaskGDSList, CM::Chart{I,M}, CNs::ChartList,
        robot_coord_rep=ChartRep(); log_task_chart=false) where {M,I}
    m = dim(M)
    JftAJf_sum = zeros(m,m)
    JftA_a_sum = zeros(m)
    CNs_out = ChartList()

    (xm, vm) = robot_coord_rep == EmbRep() ? emb_to_chart_differential(xm, vm, CM) : (xm, vm)
    for i in 1:length(tasks)
        CN = isglobal(CNs[i]) ? CNs[i] : choose_chart_chart(xm, tasks[i], CM, CNs[i])
        JftAJf, JftA_a = single_task_components(xm, vm, tasks[i], CM, CN)
        JftAJf_sum += JftAJf
        JftA_a_sum += JftA_a
        log_task_chart && push!(CNs_out, CN)
    end
    σxddot = SMatrix{m,m,eltype(xm)}(pinv(Matrix(JftAJf_sum)))*JftA_a_sum
    if robot_coord_rep == EmbRep()
        xme, vme, σxddot = chart_to_emb_differential(xm, vm, σxddot, CM)
    end
    σxddot, CNs_out
end


text:
export vertices, bonds, distance, unit, isconnected
export viz

"""
    AbstractSites

Abstrat type for all lattice or non-lattice atom collections.
The interface includes

* `vertices`
* `bonds`
* `Base.getindex`
"""
abstract type AbstractSites end
abstract type AbstractLattice <: AbstractSites end

"""
    vertices(sites)

Get the vertices for the site collection.
"""
function vertices end

"""
    isconnected(lattice, i, j; kwargs...)

return true if sites `i` and `j` are connected.
"""
function bonds end

"""
    bond(sites, i, j)

Get the bond (a pair of tuples) for the site `i` and `j`.
One can also access it with `sites[i; j]`.
"""
function bond(lt::AbstractSites, loc1::Tuple, loc2::Tuple)
    lt[loc1...], lt[loc2...]
end
function bond(lt::AbstractSites, loc1, loc2)
    lt[loc1], lt[loc2]
end
function bond(lt::AbstractSites, loc1::CartesianIndex, loc2::CartesianIndex)
    lt[loc1.I...], lt[loc2.I...]
end

Base.typed_vcat(lt::AbstractSites, loc1, loc2) = bond(lt, loc1, loc2)

distance(i::Tuple{T,T}, j::Tuple{T,T}) where T<:Real = sqrt((i[1] - j[1])^2 + (i[2] - j[2])^2)
Base.length(lt::AbstractSites) = length(vertices(lt))
Base.CartesianIndices(lt::AbstractLattice) = CartesianIndices(size(lt))
Base.LinearIndices(lt::AbstractLattice) = LinearIndices(size(lt))

"""
    bonds(atoms)

Get the bonds for the atom collection.
"""
function bonds(ud::AbstractSites)
    edges = Tuple{Int,Int}[]
    n = length(ud)
    for i = 1:n
        for j = i+1:n
            if isconnected(ud, i, j)
                push!(edges, (i,j))
            end
        end
    end
    return edges
end

function viz(lt::AbstractSites; line_style=bondstyle(:default, stroke("black"), linewidth(3*unit(lt)*mm)),
             node_style=nodestyle(:default, r=unit(lt)*0.2, stroke("black"), fill("white"), linewidth(3*unit(lt)*mm)),
             text_style=textstyle(:default, fontsize(unit(lt)*100pt)), labels=vertices(lt))
    empty_cache!()
    for (i,node) in enumerate(vertices(lt))
        node_style >> lt[node]
        text_style >> (lt[node], "$(labels[i])")
    end
    for bond in bonds(lt)
        line_style >> lt[bond[1]; bond[2]]
    end
    flush!()
end

function Base.show(io::IO, mime::MIME"text/html", lt::AbstractSites)
    Base.show(io, mime, viz(lt))
end

include("unitdisk.jl")
include("parallelogram.jl")
include("squarelattice.jl")
include("chimera.jl")


text:
#EDAData: TR0 Reader Interface
#-------------------------------------------------------------------------------

#==Main data structures
===============================================================================#
mutable struct Tr0Reader
	reader::SpiceData.DataReader
end

#==Open/read/close functions
===============================================================================#
function Base.open(::Type{Tr0Reader}, path::String)
	return Tr0Reader(SpiceData._open(path))
end

function Base.read(r::Tr0Reader, signame::String)
	y = read(r.reader, signame)
	return DataF1(r.reader.sweep, y)
end

Base.close(r::Tr0Reader) = close(r.reader)

#==Accessors
===============================================================================#
Base.names(r::Tr0Reader) = r.reader.signalnames

#==Generate friendly show functions
===============================================================================#
function Base.show(io::IO, r::Tr0Reader)
	print(io, "Tr0Reader:", r.reader)
end

#==High-level interface
===============================================================================#
open_tr0(path::String) = open(Tr0Reader, path) #Use Tr0 reader from this module.
open_tr0(fn::Function, path::String) = open(fn, Tr0Reader, path) #For do/end method

#Last line


text:
@actions app begin
    @get (ctx) begin
        @text "awesome"
    end
end


text:
"""
    FrankCopula{d,T}

Fields:
  - θ::Real - parameter

Constructor

    FrankCopula(d, θ)

The [Frank](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Most_important_Archimedean_copulas) copula in dimension ``d`` is parameterized by ``\\theta \\in [0,\\infty)``. It is an Archimedean copula with generator :

```math
\\phi(t) = -\\frac{\\log\\left(1+e^{-t}(e^{-\\theta-1})\\right)}{\theta}
```
"""
struct FrankCopula{d,T} <: ArchimedeanCopula{d}
    θ::T
end
FrankCopula(d,θ) = θ >= 0 ? FrankCopula{d,typeof(θ)}(θ) : @error "Theta must be positive"
ϕ(  C::FrankCopula,       t) = -log(1+exp(-t)*(exp(-C.θ)-1))/C.θ
ϕ⁻¹(C::FrankCopula,       t) = -log((exp(-t*C.θ)-1)/(exp(-C.θ)-1))

D₁ = GSL.sf_debye_1 # sadly, this is C code.
# could be replaced by :
# using QuadGK
# D₁(x) = quadgk(t -> t/(exp(t)-1), 0, x)[1]/x
# to make it more general. but once gain, it requires changing the integrator at each evlauation,
# which is problematic.
# Better option is to try to include this function into SpecialFunctions.jl.

τ(C::FrankCopula) = 1+4(D₁(C.θ)-1)/C.θ
function τ⁻¹(::Type{FrankCopula},τ)
    if τ == zero(τ)
        return τ
    end
    x₀ = (1-τ)/4
    return Roots.fzero(x -> (1-D₁(x))/x - x₀, 10)
end

radial_dist(C::FrankCopula) = Logarithmic(1-exp(-C.θ))


text:
function get_kegg_organism_codes()::Some

    try

        # initialize -
        organism_array = Array{KEGGOrganism,1}()

        # formulate the URL -
        url_string = "$(_KEGG_LIST_URL)/organism"

        # make the call -
        http_body = http_get_call_with_url(url_string) |> check
        if (isnothing(http_body) == true)
            return Some(nothing)
        end

        # split along the newline -
        tokenized_body = split(http_body,"\n")
        for record in tokenized_body

            # record -
            # {Txxx}\t{organism code}\t{description}\t{species_taxonomy}

            # ok, we need to cut again, along the \t -
            field_array = split(record,"\t")

            # check, do we have 4 elements -
            if (length(field_array) == 4)

                # create new organism -
                organism_wrapper = KEGGOrganism()
                organism_wrapper.organism_id = field_array[1]
                organism_wrapper.species_description = field_array[3]
                organism_wrapper.species_taxonomy = field_array[4]

                # # grab the organism_code -
                organism_code = field_array[2]
                organism_wrapper.organism_code = organism_code

                # cache -
                push!(organism_array, organism_wrapper)
            end
        end

        # return -
        return Some(organism_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_genes_in_organism_pathway(organism_code::String, pathway_code::String)::Some

    try

        # initialize -
        gene_list = Array{String,1}()

        # setup url -
        url_string = "$(_KEGG_LINK_URL)/$(organism_code)/$(pathway_code)"
        http_body = http_get_call_with_url(url_string) |> check
        if (isnothing(http_body) == true)
            return Some(nothing)
        end

        # split along the newline -
        tokenized_body = split(http_body, "\n")
        for token in tokenized_body
            fragment_array = split(token, "\t")
            if (length(fragment_array) > 1)
                gene_value = string(fragment_array[2])
                push!(gene_list, gene_value)
            end
        end

        # return -
        return Some(gene_list)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_pathways_for_organism(organism_code::String)::Some

    try

        # initialize -
        pathway_array = Array{KEGGPathway,1}()

        # formulate the url -
        url_string = "$(_KEGG_LIST_URL)/pathway/$(organism_code)"

        # get http -
        http_body = http_get_call_with_url(url_string) |> check
        if (isnothing(http_body) == true)
            return Some(nothing)
        end

        # split along the newline -
        tokenized_body = split(http_body, "\n")
        for token in tokenized_body

            # split along the \t -
            fragment_array = split(string(token), "\t")
            if (length(fragment_array) == 2)

                # build -
                pathway_object = KEGGPathway()
                pathway_object.pathway_id = string(fragment_array[1])
                pathway_object.pathway_description = string(fragment_array[2])

                # cache -
                push!(pathway_array, pathway_object)
            end
        end

        # return -
        return Some(pathway_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end


text:
# Helper for writing an AST to a string! Extended by e.g. Transpiler!

@compat abstract type ASTIO <: IO end

Base.print(io::ASTIO, args...) = print(io.io, args...)
Base.print(io::ASTIO, args::TypedSlot) = print(io.io, args)
Base.print(io::ASTIO, arg::String) = print(io.io, arg)
Base.print(io::ASTIO, arg::Char) = print(io.io, arg)
Base.print(io::ASTIO, arg::Symbol) = print(io.io, arg)
Base.print(io::ASTIO, arg::Float32) = print(io.io, arg)

Base.write(io::ASTIO, args...) = write(io.io, args...)
Base.write(io::ASTIO, arg::UInt8) = write(io.io, arg)
Base.write(io::ASTIO, arg::String) = write(io.io, arg)
Base.write(io::ASTIO, arg::Char) = write(io.io, arg)

immutable ExprNotSupported
    message::String
    line::Int
end

unsupported_expr(message, line) = throw(ExprNotSupported(message, line))


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

using Markdown
using InteractiveUtils

# ╔═╡ 1389c991-523c-46ee-8d54-7c24f2a8ab8b
md"""
# Health Standards Relevant to this Project

## [Discharge Diagnosis Codes](https://www.who.int/standards/classifications/classification-of-diseases)

>ICD is the foundation for the identification of health trends and statistics globally, and the international standard for reporting diseases and health conditions. It is the diagnostic classification standard for all clinical and research purposes. ICD defines the universe of diseases, disorders, injuries and other related health conditions, listed in a comprehensive, hierarchical fashion.

Upon discharge from a hospitalization, patients are assigned diagnostic codes by professional coders who read the discharge summary and based on the contained information assign a primary diagnosis as well as a number of associated diagnoses using ICD-10 codes.

These codes have the benefit of being widely used in hospitals and thus widely available, but they suffer from variability in who generates/enters the codes and, more importantly, suffer from up coding (purposelly selecting a code that reimburses at a higher rate rather than the most clinically accurate code), as the ICD codes are used to determine hospital funding in Australia (and elsewhere). (See, for example, [Victoria case-mix funding](https://www2.health.vic.gov.au/hospitals-and-health-services/funding-performance-accountability/activity-based-funding/casemix-funding).) Conversely, since the codes are not used for GP funding in Australia, they are not widely available from GP clinics.

## Lab Tests

In dominant standard for coding lab tests and observations is [LOINC](https://bioportal.bioontology.org/ontologies/LOINC).

## Prescriptions

In Australia prescription drugs available on the [Pharmaceutical Benefit Scheme (PBS)](https://www.pbs.gov.au/info/general/about-this-site) are coded using the WHO [Anatomical Therapeutic Chemical (ATC) codes combined with the Defined Daily Dose (DDD) codes](https://www.whocc.no/atc_ddd_index/). The WHO provides a searchable interface for the ATC/DDD codes. You can browse the ATC code heirarchy at [Bioportal](https://bioportal.bioontology.org/ontologies/ATC), but without the DDD codes.

"""

# ╔═╡ 43f379e2-133c-4ec0-a410-8432ebc670f1
html"""

 Limitations of PBS data

Monash University epidemiologist Jenni Ilomaki provided the following limitations to PBS data:

Apologies for my late reply. There are a couple of limitations in the Australian PBS data that I would like to have (in the ideal world)

 PBS claim records do not have instructions which makes it difficult to estimate the duration of prescription. In the US Medicare, and also in some Asian countries, the prescribed duration is available. In Australia, we have to estimate this based on the dispensed amount and the usual administration frequency and dose.

 I would like to have data which includes all dispensed medications, not only those dispensed via the PBS. There are lots of prescriptions dispensed as 'private'. Discount pharmacies give big discounts in these cases, for example, waiving the dispensing fee.

Indication for the medication would be also important, but this does not occur in any claims data that I am aware of so totally requires the 'ideal world'.

 Finally, it would be great to have a flag when medication is stopped. Or if a new, similar medication is dispensed, it would be great to know if the new one is a switch or an add-on. Again, I know I am dreaming....

"""

# ╔═╡ Cell order:
# ╠═1389c991-523c-46ee-8d54-7c24f2a8ab8b
# ╠═43f379e2-133c-4ec0-a410-8432ebc670f1


text:
"""
mutable struct InverseDistanceSurrogate{X,Y,P,C,L,U} <: AbstractSurrogate

The inverse distance weighting model is an interpolating method and the
unknown points are calculated with a weighted average of the sampling points.
p is a positive real number called the power parameter.
p > 1 is needed for the derivative to be continuous.

"""
mutable struct InverseDistanceSurrogate{X,Y,L,U,P} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    p::P
end

function InverseDistanceSurrogate(x,y,lb,ub;p::Number=1.0)
    return InverseDistanceSurrogate(x,y,lb,ub,p)
end

function (inverSurr::InverseDistanceSurrogate)(val)
    if val in inverSurr.x
        return inverSurr.y[findfirst(x->x==val,inverSurr.x)]
    else
        if length(inverSurr.lb) == 1
            num = sum(inverSurr.y[i]*(norm(val .- inverSurr.x[i]))^(-inverSurr.p) for i = 1:length(inverSurr.x))
            den = sum(norm(val .- inverSurr.x[i])^(-inverSurr.p) for i = 1:length(inverSurr.x))
            return num/den
        else
            βᵢ = [norm(val .- inverSurr.x[i])^(-inverSurr.p) for i = 1:length(inverSurr.x)]
            num = sum(inverSurr.y[i] * βᵢ[i] for i = 1:length(inverSurr.y))
            den = sum(βᵢ)
            return num/den
        end
    end
end

function add_point!(inverSurr::InverseDistanceSurrogate,x_new,y_new)
    if eltype(x_new) == eltype(inverSurr.x)
        #1D
        append!(inverSurr.x,x_new)
        append!(inverSurr.y,y_new)
    else
        #ND
        push!(inverSurr.x, x_new)
        push!(inverSurr.y, y_new)
    end
    nothing
end


text:
export RandIntNumSet

"""
    curveGrad(i::Int64, f::Array{Float64,1}; h::Int64=5, di::Int64=2)
    ->
    [grad,[iLB, iHB, dx]]
Numerical differentiation (5-point method) of the input function f which is in terms of an discrete array.
    \n== Arguments ==
    \n`i::Int64`: The point (index) for derivative calcualtion.
    \n`f::Array{Float64,1}`: The array of function values.
    \n`h::Int64`: The interval h. Default value is 5.
    \n`di::Int64`: An addtional interval for mean value at each point to compenstate the data fluctuation. Default value is 2, which means the mean of 2*2+1 sequential values.
    \n== Output variables ==
    \n`grad::Float64`: Derivative of f at point (index) i.
    \n`iLB::Int64`: Lower bound requirement for the index.
    \n`iHB::Int64`: Upper bound requirement for the index.
    \n`dx::Int64`: half interval length. iHB = i+dx iLB = i-dx
"""
function curveGrad(i::Int64, f::Array{Float64,1}; h::Int64=5, di::Int64=2)
    dx = 2h+di
    iHB = i+dx
    iLB = i-dx
    f_p2 = f[i+2h-di:iHB] |> mean
    f_p1 = f[i+h-di:i+h+di] |> mean
    f_m1 = f[i-h-di:i-h+di] |> mean
    f_m2 = f[iLB:i-2h+di] |> mean
    grad = ( - f_p2 + 8*f_p1 - 8*f_m1 + f_m2 ) / 12*h
    [grad,[iLB, iHB, dx]]
end

"""
    HofNPP(numSet::Array{Int64,1}, Coeff::Float64=1.0)
    ->
    H
Generating the corresponding Hamiltonian given a number set.
    \n== Arguments ==
    \n`numSet::Array{Int64,1}`: The target number set.
    \n`Coeff::Float64`: An overall coefficent. Default value is 1.0.
    \n== Output variables ==
    \n`H:CompositeBlock`: The corresponding Hamiltonian.
"""
function HofNPP(numSet::Array{Int64,1}, Coeff::Float64=1.0)
    nbit = length(numSet)
    sz = i->put(nbit, i=>Z)
    H = map(1:nbit) do i
            numSet[i]*sz(i)
        end |> sum
    H = Coeff*H^2
end

"""
    RandIntNumSet(size::Int64, sum::Int64)
    ->
    set::Array{Int64,1}
Randomly generating number set with equi-partition soluttions.
    \n== Arguments ==
    \n`size::Int64`: The size (number of elements) of the target number set.
    \n`sum::Int64`: The sum of elements inside the target number set.
    \n`isShuffled::Bool`: Shuffle the elements inside the set to eliminate a trivial solution for the equi-partition. Default value is true.
    \n== Output variables ==
    \n`set::Array{Int64,1}`: The generated number set.
"""
function RandIntNumSet(size::Int64, sum::Int64; isShuffled::Bool=true)
    indexSet = [1:sum;]
    middle = Int(round(sum/2))
    deleteat!(indexSet, (middle, sum))
    selection = rand(indexSet, size-2)
    append!(selection, [0,middle,sum]) |> sort!
    set = map(2:size+1) do i
            selection[i]-selection[i-1]
        end
    if isShuffled == true
        set |> shuffle!
    end
    set
end


text:
using Plots
import Statistics
using CausalTrees

"""
logistic
"""
logistic(α) = exp(α) / (1 + exp(α))

# TEST 1 ######################################################################
# generate data
x = randn(100)
y = rand(100)
y[x .> 0] .+=  2
# plot
scatter(x,y,label="",legend=:bottomright)
# test split finding algorithm: split #1
@time split_i,thres,loss = CausalTrees.evaluate_possible_splits(x[:,1],y,10)
#@time pos_thres = CausalTrees.compute_possible_splits(x[:,1],10)
vline!([thres], label = "split #1")
# test split finding algorithm: split #2
msk = x[:,1] .< 0.0
@time split_i,thres,loss = CausalTrees.evaluate_possible_splits(x[msk,1],y[msk],10)
vline!([thres], label = "split #2")

# TEST 2 #######################################################################
# generate data
x = randn(1000,1)
y = rand(1000)
y[x[:,1] .> 0] .+=  2
y[x[:,1] .> 1] .+=  1
# plot
scatter(x[:,1],y,label="")
# fit tree and predict
@time t1 = CausalTrees.decision_tree(x,y,mse_tau,min_leaf_size=100)
yhat = predict(t1,x)
scatter!(x[:,1],yhat,label="")

# TEST 3 #######################################################################
# generate data
x = randn(1000,5)
e = 6 .* rand(1000)
y = 2 .* x[:,1] .+ e
y[x[:,1] .> 0] .+=  1.5
# plot
scatter(x[:,1],y,label="")
# fit tree and predict
@time t1 = CausalTrees.decision_tree(x,y,mse_tau,min_leaf_size=100)
yhat = predict(t1,x)
scatter!(x[:,1],yhat,label="")

# TEST 4 #######################################################################
# generate data
x = randn(1000,1)
y = rand(1000)
y[x[:,1] .> 0] .+=  2
y[x[:,1] .> 1] .+=  1
# plot
scatter(x[:,1],y,label="",legend=:bottomright)
# fit tree and predict
@time t1 = CausalTrees.decision_tree(x,y,mse_tau,min_leaf_size=50) # not honest
@time t2 = CausalTrees.decision_tree(x,y,mse_tau,min_leaf_size=50,honesty=true) # honest
yhat1 = predict(t1,x)
yhat2 = predict(t2,x)
scatter!(x[:,1],yhat1,label="not honest")
scatter!(x[:,1],yhat2,label="honest")

# TEST 5 #######################################################################
# generate data
x = randn(1000,5)
e = 6 .* rand(1000)
y = 2 .* x[:,1] .+ e
y[x[:,1] .> 0] .+=  1.5
# plot
scatter(x[:,1],y,label="")
# fit random forest and predict
@time rf1 = CausalTrees.random_forest(x,y,mse_tau,n_trees=500,row_subsamp_p=0.8,
    min_leaf_size=30,honesty=true)
yhat = predict(rf1,x)
scatter!(x[:,1],yhat,label="")

# TEST 6 #######################################################################
# binary data
n = 2000
# covariates
x = randn((n,5))
# treatment group
w = rand([0,1],n)
# part of treatment effect
α = 1.8 .* x[:,1] + 0.8 .* x[:,2]
yp = logistic.(α + 0.01 .* randn(n))
y = zeros(n)
rr = rand(n)
y[yp .> rr] .= 1.0
y = reshape(y,n)
# plot
scatter(x[:,1],y,label="")
# fit random forest and predict
@time rf1 = CausalTrees.random_forest(x,y,mse_tau,n_trees=1000,row_subsamp_p=0.8,
    min_leaf_size=30,honesty=false)

yhat = predict(rf1,x)
scatter!(x[:,1],yhat,label="")
scatter!(x[:,1],logistic.(0.6 .* x[:,1] + 0.5 .* x[:,2]),label="")

scatter(yp,yhat,label="")
plot!(Shape([(0,0),(1,1)]),label="")


text:
## 2D cases
#==========================================================================================#

"""
    solve_case(foil :: Foil, uniform :: Uniform2D; sources = false, wake_length = 1e3, num_panels :: Integer = 60)

Evaluate a doublet-source case given a `Foil` with a `Uniform2D`, with optional named arguments to specify whether the source terms are non-zero, the length of the wake, and the number of panels for the analysis.
"""
function solve_case(foil :: Foil, freestream :: Uniform2D; viscous = false, sources = false, wake_length = 1e5, num_panels :: Integer = 60, wake_panels = 15)
    panels = paneller(foil, num_panels)
    cls, cms, cps, cl_wake = solve_problem(panels, velocity(freestream), freestream.angle, ifelse(viscous, wake_panels, sources), wake_length)
    # cdp, cl = sincos(freestream.angle) .* cl_wake
    # cm      = cl * r_ref
    # coeffs    = (cdp, cl, cm)

    cl_wake, cls, cms, cps, panels
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# Autogenerated wrapper script for Ninja_jll for armv7l-linux-gnueabihf-cxx11
export ninja

JLLWrappers.@generate_wrapper_header("Ninja")
JLLWrappers.@declare_executable_product(ninja)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        ninja,
        "bin/ninja",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ClusterManager for Slurm

export SlurmManager, addprocs_slurm

import Logging.@warn

struct SlurmManager <: ClusterManager
    np::Integer
end

function launch(manager::SlurmManager, params::Dict, instances_arr::Array,
                c::Condition)
    try
        exehome = params[:dir]
        exename = params[:exename]
        exeflags = params[:exeflags]

        stdkeys = keys(Distributed.default_addprocs_params())

	println(stdkeys)
        p = filter(x->(!(x[1] in stdkeys) && x[1] != :job_file_loc), params)
	println(p)

        srunargs = []
        for k in keys(p)
            if length(string(k)) == 1
                push!(srunargs, "-$k")
                val = p[k]
                if length(val) > 0
                    push!(srunargs, "$(p[k])")
                end
            else
                k2 = replace(string(k), "_"=>"-")
                val = p[k]
                if length(val) > 0
                    push!(srunargs, "--$(k2)=$(p[k])")
                else
                    push!(srunargs, "--$(k2)")
                end
            end
        end

        # Get job file location from parameter dictionary.
        job_file_loc = joinpath(exehome, get(params, :job_file_loc, "."))

        # Make directory if not already made.
        if !isdir(job_file_loc)
            mkdir(job_file_loc)
        end

        println("removing old files")
        # cleanup old files
	map(f->rm(joinpath(job_file_loc, f)), filter(t -> occursin(r"job(.*?).out", t), readdir(job_file_loc)))
        println("removing old Setting up srun commands")

        np = manager.np
        jobname = "julia-$(getpid())"
        srun_cmd = `srun -J $jobname -n $np -o "$(joinpath(job_file_loc, "job%4t.out"))" -D $exehome $(srunargs) $exename $exeflags $(worker_arg())`
        srun_proc = open(srun_cmd)
        for i = 0:np - 1
            print("connecting to worker $(i + 1) out of $np\r")
            local w=[]
            fn = "$(joinpath(exehome, job_file_loc))/job$(lpad(i, 4, "0")).out"
            t0 = time()
            while true
                if time() > t0 + 60 + np
                    @warn "dropping worker: file not created in $(60 + np) seconds"
                    break
                end
                sleep(0.001)
                if isfile(fn) && filesize(fn) > 0
                    w = open(fn) do f
                        return split(split(readline(f), ":")[2], "#")
                    end
                    break
                end
            end
            if length(w) > 0
                config = WorkerConfig()
                config.port = parse(Int, w[1])
                config.host = strip(w[2])
                # Keep a reference to the proc, so it's properly closed once
                # the last worker exits.
                config.userdata = srun_proc
                push!(instances_arr, config)
                notify(c)
            end
        end
    catch e
        println("Error launching Slurm job:")
        rethrow(e)
    end
end

function manage(manager::SlurmManager, id::Integer, config::WorkerConfig,
                op::Symbol)
    # This function needs to exist, but so far we don't do anything
end

addprocs_slurm(np::Integer; kwargs...) = addprocs(SlurmManager(np); kwargs...)


max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.9 trig^m (a+b cot^n+c cot^(2 n))^p.jl


text:
######################################################################
# domain.jl
#
# TODO: allow creation from tuples as well
######################################################################
export Domain

type Domain
	# user provided
	mins::Vector{Float64}			# minimum value per dimension
	maxes::Vector{Float64}			# maximum value per dimension
	cells::Vector{Int}				# number of cells per dimension

	# calculated and used internally
	lengths::Vector{Float64}		# maxes - mins
	cell_lengths::Vector{Float64}	# length of cell in each dimension
	num_dims::Int					# number of dimensions
	cell_size::Float64				# size (area, volume, etc) of a cell

	function Domain(mins, maxes, cells::Vector{Int})
		num_dims = length(mins)
		@assert num_dims == length(maxes)
		@assert num_dims == length(cells)
		@assert num_dims > 0

		d = new()

		# these things stay the same
		d.mins = deepcopy(mins)
		d.maxes = deepcopy(maxes)
		d.cells = deepcopy(cells)

		# these can be easily computed
		d.num_dims = num_dims
		d.lengths = zeros(num_dims)
		d.cell_lengths = zeros(num_dims)
		d.cell_size = 1.0
		for i = 1:num_dims
			d.lengths[i] = d.maxes[i] - d.mins[i]
			d.cell_lengths[i] = d.lengths[i] / d.cells[i]
			d.cell_size *= d.cell_lengths[i]
		end
		return d
	end

	# can just pass an int for num_cells too
	function Domain(mins, maxes, num_cells::Int)
		n = length(mins)
		return Domain(mins, maxes, num_cells*ones(Int, n))
	end

	# Allow user to only give maxes
	function Domain(maxes, bins::Vector{Int})
		return Domain(zeros(length(maxes)), maxes, bins)
	end

	function Domain(maxes, num_cells::Int)
		return Domain(zeros(length(maxes)), maxes, num_cells)
	end

	# allow user to only specify discretization (unit cell is assumed)
	Domain(num_cells::Int) = Domain([0,0], [1,1], [num_cells,num_cells])

end

x_min(domain::Domain) = domain.mins[1]
y_min(domain::Domain) = domain.mins[2]
z_min(domain::Domain) = domain.mins[3]

x_max(domain::Domain) = domain.maxes[1]
y_max(domain::Domain) = domain.maxes[2]
z_max(domain::Domain) = domain.maxes[3]

x_size(domain::Domain) = domain.cell_lengths[1]
y_size(domain::Domain) = domain.cell_lengths[2]
z_size(domain::Domain) = domain.cell_lengths[3]

x_cells(domain::Domain) = domain.cells[1]
y_cells(domain::Domain) = domain.cells[2]
z_cells(domain::Domain) = domain.cells[3]


text:
using Test
import QHull

@testset "chull" begin
    pts = [-1.0 0;
           1 1;
           3 0;
           3 3;
           2 2;
           0 3;
           -1 2]

    hull = QHull.chull(pts)
    @test hull.vertices == [1, 3, 4, 6, 7]
    @test size(hull.points) == size(pts)
    @test hull.simplices == Array{Int,1}[[3,1], [4,3], [6,4], [7,1], [7,6]]

    # multi-dim
    x = randn(1000, 5)
    hull = QHull.chull(x)
end

@testset "Polyhedra tests" begin
    include("polyhedron.jl")
end


text:
# Tests for mixed arithmetic
@testset "Mixed" begin

@testset "Mixed Arithmetic" begin
    @test 1.01USD + Monetary(:USD; storage=BigInt) == 2.01USD
    @test 1.01USD ÷ Monetary(:USD; storage=BigInt) == 1
    @test 1.01USD % Monetary(:USD; storage=BigInt) == 0.01USD
end

@testset "Mixed Comparison" begin
    @test 1.01USD > Monetary(:USD; storage=BigInt)
    @test 1.01USD > Monetary(:USD; storage=BigInt, precision=4)
    @test 0.99USD < Monetary(:USD; storage=BigInt)
    @test 0.99USD < Monetary(:USD; storage=BigInt, precision=4)
    @test USD == Monetary(:USD; storage=BigInt)
    @test USD == Monetary(:USD; storage=BigInt, precision=4)

    @test USD ≠ CAD
    @test 0USD == 0CAD  # NB: unexpected?
end

@testset "Type Safety" begin
    @test_throws MethodError 1.01USD > CAD
    @test_throws MethodError 1.01USD + CAD
end

end  # testset Monetary


text:
using TensorToolbox
using LinearAlgebra
using InvertedIndices
using StatsBase

"""
This function gives the best rank-1 approximation of the non-negative tensor `T`
"""
function LT1R(T)
    input_tensor_depth = ndims(T)
    input_tensor_shape = size(T)

    partial_sums = []
    one_to_N = [1:input_tensor_depth;]
    for k=1:input_tensor_depth
        partial_sum = vec(sum(T, dims=one_to_N[Not(k)]))
        push!(partial_sums, partial_sum)
    end

    P = ttt(partial_sums[1], partial_sums[2])
    for n=3:input_tensor_depth
        P = ttt(P, partial_sums[n])
    end

    dev = sum(T)^(input_tensor_depth-1)
    P .= P ./ dev
    return P
end

function get_basis_positions(tensor_size, reqrank)
    # Non bingo position become basis for reconstration tensor.
    # Basis position is non-bingo posision
    D = length(reqrank)
    basis_positions = []
    for d = 1 : D
        basis_position = [1; sample(2:tensor_size[d], reqrank[d]-1, replace=false, ordered=true)]
        push!(basis_positions, basis_position)
    end
    return basis_positions
end

"""
Legendre Tucker-rank Reduction,
proposed by Kazu Ghalamkari and Mahito Sugiyama.

# Aruguments
- `X`: input non-negative tensor, multi dimensional array.
- `reqrank` : target Tucker rank, array.
"""
function LTR(X, reqrank)
    X = copy(X)
    @assert ndims(X) == length(reqrank)
    D = ndims(X)
    tensor_size = size(X)
    basis_positions = get_basis_positions(tensor_size, reqrank)

    for d=1:D
        for b_number = 1 : reqrank[d]
            rank1_subtensor_begin = basis_positions[d][b_number]
            if b_number != reqrank[d]
                rank1_subtensor_end = basis_positions[d][b_number+1] - 1
            else
                rank1_subtensor_end = tensor_size[d]
            end

            if rank1_subtensor_begin == rank1_subtensor_end
                continue
            end
            idx = ntuple( l -> l == d ? (rank1_subtensor_begin:rank1_subtensor_end) : (:), D)

            X[idx...] = LT1R( X[idx...] )
        end
    end
    return X
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
export gradient_numerical

_expval(psi::AbstractVector, Ĥ, op::AbsLinearOperator) =
    psi'*Ĥ*psi/(psi'psi)

_expval(psi::AbstractVector, Ĥ, op::KLocalLiouvillian) =
    sum(abs2.((Ĥ*psi)./psi))

function gradient_numerical(net, op, h=0.001)
    Ĥ = Matrix(op)
    cnet = cached(net)
    hilb = basis(op)
    psi = Vector(cnet, hilb)

    ∇g = grad_cache(net)
    ∇gv = vec_data(∇g)[1]
    Δ = grad_cache(eltype(trainable_first(net)), net)
    Δv = vec_data(Δ)[1]

    optimizer = Optimisers.Descent(1)

    if eltype(Δv) isa Complex
        h = h+h*1im
    end

    Δv .= 0.0
    for i=1:length(Δv)
        Δv[i] = -h
        Optimisers.update!(optimizer, net, Δ)
        E_p = _expval(Vector(cnet, hilb), Ĥ, op)

        Δv[i] = 2h
        Optimisers.update!(optimizer, net, Δ)
        E_m = _expval(Vector(cnet, hilb), Ĥ, op)

        Δv[i] = -h
        Optimisers.update!(optimizer, net, Δ)
        Δv[i] = 0

        ∇gv[i] = (E_p - E_m) / (2h)
    end
    return ∇g
end