text:
if Sys.iswindows()
    const glade_path = string(dirname(Base.source_path()),"\\interface.glade")
else
    const glade_path = string(dirname(Base.source_path()),"/interface.glade")
end

function makegui(r::RHD2000,s,task,fpga)

    b = Builder(filename=glade_path)

    #GUI ARRANGEMENT
    grid = Grid()

    #COLUMN 1 - control buttons

    #ROW 1

    #ROW 2
    vbox1_2=Grid()
    scroll_left=ScrolledWindow()
    Gtk.GAccessor.policy(scroll_left,Gtk.GConstants.GtkPolicyType.NEVER,Gtk.GConstants.GtkPolicyType.AUTOMATIC)
    push!(scroll_left,vbox1_2)
    grid[1,2]=scroll_left

    frame_control=Frame("Control")
    grid[1,1]=frame_control
    vbox_control = Grid()
    push!(frame_control,vbox_control)

    button_init = Button()
    add_button_label(button_init,"Init")
    vbox_control[1,1]=button_init

    button_run = ToggleButton()
    add_button_label(button_run,"Run")
    vbox_control[2,1]=button_run

    button_record = ToggleButton()
    add_button_label(button_record,"Record")
    vbox_control[3,1]=button_record

    #GAIN
    frame1_2=Frame("Gain")
    vbox1_2[1,2]=frame1_2
    vbox1_2_1=Grid()
    push!(frame1_2,vbox1_2_1)

    sb2=SpinButton(1:1000)
    Gtk.GAccessor.value(sb2,1)
    vbox1_2_1[1,1]=sb2

    gain_checkbox=CheckButton()
    add_button_label(gain_checkbox," x 10")
    vbox1_2_1[2,1]=gain_checkbox

    button_gain = CheckButton()
    add_button_label(button_gain,"All Channels")
    Gtk.GAccessor.active(button_gain,false)
    vbox1_2_1[1,2]=button_gain

    #THRESHOLD
    (thres_widgets,frame1_3,vbox_slider)=_make_thres_gui()
    vbox1_2[1,3]=frame1_3

    #SPIKE
    frame_hold=Frame("Spike")
    vbox1_2[1,4]=frame_hold
    vbox_hold=Grid()
    push!(frame_hold,vbox_hold)

    button_pause=ToggleButton()
    add_button_label(button_pause,"Pause")
    vbox_hold[2,2]=button_pause

    button_clear=Button()
    add_button_label(button_clear,"Refresh")
    vbox_hold[1,2]=button_clear

    button_restore=Button()
    add_button_label(button_restore,"Restore")
    vbox_hold[3,2]=button_restore

    if VERSION > v"0.7-"
        button_rb=Array{RadioButton}(undef, 3)
    else
        button_rb=Array{RadioButton}(3)
    end
    button_rb[1]=RadioButton(active=true)
    button_rb[2]=RadioButton(button_rb[1])
    button_rb[3]=RadioButton(button_rb[2])

    vbox_hold[1,3]=button_rb[1]
    Gtk.GAccessor.mode(button_rb[1],false)
    add_button_label(button_rb[1],"RubberBand")
    vbox_hold[2,3]=button_rb[2]
    Gtk.GAccessor.mode(button_rb[2],false)
    add_button_label(button_rb[2],"Draw")
    vbox_hold[3,3]=button_rb[3]
    Gtk.GAccessor.mode(button_rb[3],false)
    add_button_label(button_rb[3],"Selection")

    #CLUSTER
    frame1_4=Frame("Clustering")

    vbox1_3_2=Grid()
    push!(frame1_4,vbox1_3_2)

    button_sort1 = Button()
    button_sort2 = Button()
    button_sort3 = Button()

    button_sort4 = ToggleButton()
    button_sort5 = Button()

    check_sort1 = CheckButton()

    slider_sort = Scale(false, 0.0, 2.0,.02)
    adj_sort = Adjustment(slider_sort)
    Gtk.GAccessor.value(adj_sort,1.0)
    slider_sort_label=Label("Slider Label")

    sort_list=ListStore(Int32)
    push!(sort_list,(0,))
    sort_tv=TreeView(TreeModel(sort_list))
    sort_r1=CellRendererText()
    sort_c1=TreeViewColumn("Cluster",sort_r1, Dict([("text",0)]))
    Gtk.GAccessor.activate_on_single_click(sort_tv,1)

    push!(sort_tv,sort_c1)

    vbox1_3_2[1,3] = button_sort1
    vbox1_3_2[1,4] = button_sort2
    vbox1_3_2[1,5] = button_sort3
    vbox1_3_2[1,6] = button_sort4
    vbox1_3_2[1,7] = check_sort1
    #vbox1_3_2[1,7]=button_sort5
    vbox1_3_2[1,8] = slider_sort
    vbox1_3_2[1,9] = slider_sort_label

    myscroll=ScrolledWindow()
    Gtk.GAccessor.min_content_height(myscroll,150)
    Gtk.GAccessor.min_content_width(myscroll,100)
    push!(myscroll,sort_tv)
    vbox1_3_2[1,10]=myscroll
    vbox1_3_2[1,11]=Canvas(180,10)

    vbox1_2[1,5]=frame1_4 |> Gtk.showall

    #COLUMN 2 - Threshold slider

    grid[3,2]=vbox_slider

    #COLUMN 3 - MAXIMIZED CHANNEL PLOTTING

    #ROW 2
    c_grid=Grid()

    c2=Canvas()

    @guarded draw(c2) do widget
        ctx = Gtk.getgc(c2)
        SpikeSorting.clear_c2(c2,1)
    end

    show(c2)
    c_grid[1,1]=c2
    Gtk.GAccessor.hexpand(c2,true)
    Gtk.GAccessor.vexpand(c2,true)

    #ROW 2
    c3=Canvas(-1,200)
    @guarded draw(c3) do widget
        ctx = Gtk.getgc(c3)
        SpikeSorting.clear_c3(c3,1)
    end
    show(c3)
    c_grid[1,2]=c3
    Gtk.GAccessor.hexpand(c3,true)

    grid[4,2]=c_grid

    #ROW 3
    c2_slider=Scale(false, 1:16)
    adj2 = Adjustment(c2_slider)
    Gtk.GAccessor.value(adj2,1)
    grid[4,3]=c2_slider

    #COLUMN 3 - 16 CHANNEL DISPLAY

    #ROW 1 - Time
    (mytime_widgets,frame_time) = _make_clock()
    grid[5,1]=frame_time

    #ROW 2
#c=@Canvas(500,800)
c=Canvas(500)
    @guarded draw(c) do widget
        ctx = Gtk.getgc(c)
        set_source_rgb(ctx,0.0,0.0,0.0)
        set_operator(ctx,Cairo.OPERATOR_SOURCE)
        paint(ctx)
    end
    show(c)
grid[5,2]=c
    Gtk.GAccessor.vexpand(c,true)

    #ROW 3
    #Which 16 channels can be selected with a slider
    c_slider = Scale(false, 0:(div(length(r.nums)-1,16)+1))
    adj = Adjustment(c_slider)
    Gtk.GAccessor.value(adj,1)
    grid[5,3]=c_slider

    #COLUMN 4
    #ROW 2
    vbox_42=Box(:v)
    grid[6,2]=vbox_42

vbox_rb_upper=Box(:v)
    push!(vbox_42,vbox_rb_upper)

    push!(vbox_rb_upper,Label("Top Panel"))

    if VERSION > v"0.7-"
        rbs=Array{RadioButton}(undef,5)
    else
        rbs=Array{RadioButton}(5)
    end
    rbs[1]=RadioButton("16 Channel",active=true)
    rbs[2]=RadioButton(rbs[1],"32 Channel")
    rbs[3]=RadioButton(rbs[2],"64 Channel")
    rbs[4]=RadioButton(rbs[3],"64 Raster")
    rbs[5]=RadioButton(rbs[4],"Blank")

    push!(vbox_rb_upper,rbs[1])
    push!(vbox_rb_upper,rbs[2])
    push!(vbox_rb_upper,rbs[3])
    push!(vbox_rb_upper,rbs[4])
    push!(vbox_rb_upper,rbs[5])

c_rb=Canvas(40)
Gtk.GAccessor.vexpand(c_rb,true)

    push!(vbox_42,c_rb)

    vbox_rb_lower=Box(:v)
    push!(vbox_42,vbox_rb_lower)
    push!(vbox_rb_lower,Label("Lower Panel"))

    if VERSION > v"0.7-"
        rbs2=Array{RadioButton}(undef,8)
    else
        rbs2=Array{RadioButton}(8)
    end
    rbs2[1]=RadioButton("Events",active=true)
    rbs2[2]=RadioButton(rbs2[1],"16 Raster")
    rbs2[3]=RadioButton(rbs2[2],"32 Raster")
    rbs2[4]=RadioButton(rbs2[3],"Soft Scope")
rbs2[5]=RadioButton(rbs2[4],"64 Channel")
rbs2[6]=RadioButton(rbs2[5],"64 Raster")
rbs2[7]=RadioButton(rbs2[6],"Spectrogram")
rbs2[8]=RadioButton(rbs2[7],"Nothing")

    push!(vbox_rb_lower,rbs2[1])
    push!(vbox_rb_lower,rbs2[2])
    push!(vbox_rb_lower,rbs2[3])
    push!(vbox_rb_lower,rbs2[4])
    push!(vbox_rb_lower,rbs2[5])
    push!(vbox_rb_lower,rbs2[6])
push!(vbox_rb_lower,rbs2[7])
push!(vbox_rb_lower,rbs2[8])

    #MENU ITEMS

    #SORTING
    sortopts = MenuItem("_File")
    sortmenu = Menu(sortopts)
    load_sort_ = MenuItem("Load Sorting Parameters")
    push!(sortmenu,load_sort_)
    save_sort_ = MenuItem("Save Sorting Parameters")
push!(sortmenu,save_sort_)
load_backup_ = MenuItem("Recover Parameters from Backup")
push!(sortmenu,load_backup_)

saving_pref_=MenuItem("Saving Preferences")
push!(sortmenu,saving_pref_)

#Export

exopts = MenuItem("_Export")
push!(sortmenu,exopts)
    exmenu = Menu(exopts)
    export_plex_ = MenuItem("Plexon")
    push!(exmenu,export_plex_)
    export_klusta_ = MenuItem("KlustaFormat")
    push!(exmenu,export_klusta_)
    export_nwb_ = MenuItem("NWB")
    push!(exmenu,export_nwb_)
    export_jld_ = MenuItem("JLD")
    push!(exmenu,export_jld_)
export_mat_ = MenuItem("MAT")
push!(exmenu,export_mat_)

viewopts = MenuItem("_View")
viewmenu = Menu(viewopts)

    define_ref_ = MenuItem("Reference Configuration")
push!(viewmenu,define_ref_)

define_params = MenuItem("Parameter Table")
push!(viewmenu,define_params)

sv_open = MenuItem("Sort Viewer")
push!(viewmenu,sv_open)

define_ttls_ = MenuItem("TTL Configuration")
push!(viewmenu,define_ttls_)

#Options
opopts = MenuItem("_Options")
opmenu = Menu(opopts)

op_align = MenuItem("Alignment")
push!(opmenu,op_align)
op_align_menu = Menu(op_align)
op_align_min = MenuItem("Minimum")
push!(op_align_menu,op_align_min)
op_align_cross = MenuItem("Threshold Crossing")
push!(op_align_menu,op_align_cross)

op_band = MenuItem("Bandwidth")
push!(opmenu,op_band)

    mb = MenuBar()
    push!(mb,sortopts)
push!(mb,viewopts)
push!(mb,opopts)
grid[4,1]=mb

#SortView

sortview_handles = SpikeSorting.sort_gui(s[1].s.win+1)
visible(sortview_handles.win,false)

#POPUP MENUS

#Enable-Disable
    popupmenu = Menu()
    popup_enable = MenuItem("Enable")
    push!(popupmenu, popup_enable)
    popup_disable = MenuItem("Disable")
    push!(popupmenu, popup_disable)
Gtk.showall(popupmenu)

#Event
(event_handles,popup_event,popup_event_none) = _make_event_gui()

    setproperty!(grid, :column_spacing, 15)
    setproperty!(grid, :row_spacing, 15)
win = Window(grid, "Intan.jl GUI") |> Gtk.showall

#Spectrogram Menus
(spect_w_handles,spect_f_handles, popupmenu_spect) = _make_spectogram_gui()

#Soft Scope Menus
(scope_t_handles,scope_v_handles,scope_thres_handles,scope_signal_handles,popupmenu_scope) = _make_scope_gui()

prepare_save_folder(r)

    #Callback functions that interact with canvas depend on spike sorting method that is being used

    scales=ones(Float64,size(r.v,2),2) .* -.125
    scales[:,2]=scales[:,2].*.25
offs=zeros(Int64,size(r.v,2))

scope_mat=ones(Float64,500,3)

for i=1:500
    scope_mat[i,1]=550.0
    scope_mat[i,2]=650.0
    scope_mat[i,3]=750.0
end

sort_widgets=Sort_Widgets(button_sort1,button_sort2,button_sort3,button_sort4,check_sort1,false)

gain_widgets=SpikeSorting.Gain_Widgets(sb2,gain_checkbox,button_gain)
spike_widgets=Spike_Widgets(button_clear,button_pause)

#Make Filter Menu
band_widgets=_make_filter_gui()
band_widgets.lfp_en=falses(size(r.v,2))
setproperty!(band_widgets.sw_chan_sb,:upper,size(r.v,2))

#Make Table Menu
table_widgets=_make_table_gui()
for i=1:size(r.v,2)
    push!(table_widgets.list,(i,125,0,0,true))
end

#Reference Menu
ref_widgets = _make_reference_gui()
for i=1:size(r.v,2)
    push!(ref_widgets.list1,(i,))
    push!(ref_widgets.list2,(i,))
end

spect_widgets=Spectrogram(r.sr)

save_widgets=_make_save_gui()
Gtk.GAccessor.text(save_widgets.input,r.save.folder)

sleep(5.0)

sc_widgets=SpikeSorting.Single_Channel(c2,c3,Gtk.getgc(c2),copy(Gtk.getgc(c2)),false,RubberBand(Vec2(0.0,0.0),
Vec2(0.0,0.0),Vec2(0.0,0.0),[Vec2(0.0,0.0)],false,0),1,falses(500),falses(500),
false,false,button_pause,button_rb,1,(0.0,0.0),false,width(Gtk.getgc(c2)),
    height(Gtk.getgc(c2)),s[1].s.win,1.0,0.0,sortview_handles.buf,0.0,0.0,
    ClusterTemplate(convert(Int64,s[1].s.win)),0,1,false,sort_list,sort_tv,adj_sort,thres_widgets.adj,thres_widgets.slider,false,
    thres_widgets,gain_widgets)

    #Create type with handles to everything
handles=Gui_Handles(win,button_run,button_init,button_record,c_slider,adj,c2_slider,adj2,
                    c,1,1,scales,offs,(0.0,0.0),zeros(Int64,length(r.nums)),
                    0,-1 .*ones(Int64,6),
                    trues(length(r.nums)),mytime_widgets,
                    s[1].s.win,1,1,popupmenu,popup_event,popupmenu_spect,rbs,rbs2,scope_mat,
                    false,16,
                    false,slider_sort,
                    1,1,zeros(Int64,500),zeros(UInt32,20),
                    zeros(UInt32,500),zeros(Int64,50),SoftScope(r.sr,Gtk.getgc(c),SAMPLES_PER_DATA_BLOCK),
                    popupmenu_scope,sort_widgets,spike_widgets,
                    sortview_handles,band_widgets,table_widgets,spect_widgets,save_widgets,ref_widgets,
                    b,sc_widgets,sortview_handles.buf,rand(Int8,r.sr))

    handles.sc.s = -.125

Gtk.GAccessor.value(sb2,round(Int,scales[1,1]*-1000))

#=
Template Sorting Callbacks
=#

    id = signal_connect(SpikeSorting.canvas_release_template,c2,"button-release-event",Void,(Ptr{Gtk.GdkEventButton},),false,(sortview_handles.buf,sc_widgets))

    id = signal_connect(SpikeSorting.b1_cb_template,button_sort1,"clicked",Void,(),false,(sc_widgets,))
    add_button_label(button_sort1,"Delete Unit")

    id = signal_connect(SpikeSorting.b2_cb_template,button_sort2,"clicked",Void,(),false,(sc_widgets,))
    add_button_label(button_sort2,"Add Unit")

    id = signal_connect(SpikeSorting.b3_cb_template,button_sort3,"clicked",Void,(),false,(sc_widgets,sortview_handles))
    add_button_label(button_sort3,"Collect Templates")

    id = signal_connect(SpikeSorting.b4_cb_template,button_sort4,"clicked",Void,(),false,(sc_widgets,))
    add_button_label(button_sort4,"Show Template Bounds")

    setproperty!(check_sort1,:label,"Show Template")
    id = signal_connect(SpikeSorting.check_cb_template,check_sort1,"clicked",Void,(),false,(sc_widgets,))

    setproperty!(slider_sort_label,:label,"Tolerance")

id = signal_connect(SpikeSorting.unit_select_cb,sort_tv, "row-activated", Void, (Ptr{Gtk.GtkTreePath},Ptr{Gtk.GtkTreeViewColumn}), false, (handles.sc,))
id = signal_connect(SpikeSorting.canvas_press_win,c2,"button-press-event",Void,(Ptr{Gtk.GdkEventButton},),false,(handles.sc,))

#=
Window Callbacks
=#

id = signal_connect(SpikeSorting.win_resize_cb, win, "size-allocate",Void,(Ptr{Gtk.GdkRectangle},),false,(sc_widgets,))
id = signal_connect(close_cb,win, "destroy",Void,(),false,(handles,r,fpga))

#=
ISI canvas callbacks
=#
id = signal_connect(SpikeSorting.c3_press_win,c3,"button-press-event",Void,(Ptr{Gtk.GdkEventButton},),false,(sc_widgets,))

#=
Start Button Callbacks
=#

    id = signal_connect(run_cb, button_run, "clicked",Void,(),false,(handles,r,s,task,fpga))
    id = signal_connect(init_cb, button_init, "clicked", Void, (), false, (handles,r,task,fpga))
    id = signal_connect(record_cb,button_record,"clicked",Void,(),false,(handles,r))

#=
Slider callbacks
=#

    id = signal_connect(update_c1, c_slider, "value-changed", Void, (), false, (handles,))
id = signal_connect(update_c2_cb, c2_slider, "value-changed", Void, (), false, (handles,))

#=
Enable/Disable callbacks
=#

id = signal_connect(c_popup_select,c,"button-press-event",Void,(Ptr{Gtk.GdkEventButton},),false,(handles,))
id = signal_connect(popup_enable_cb,popup_enable,"activate",Void,(),false,(handles,r))
id = signal_connect(popup_disable_cb,popup_disable,"activate",Void,(),false,(handles,r))

#=
Save Load Callbacks
=#

id = signal_connect(export_plex_cb, export_plex_, "activate",Void,(),false,(handles,r))
id = signal_connect(export_jld_cb, export_jld_, "activate",Void,(),false,(handles,r))
id = signal_connect(export_mat_cb, export_mat_, "activate",Void,(),false,(handles,r))
id = signal_connect(save_config_cb, save_sort_, "activate",Void,(),false,(handles,r,s))
id = signal_connect(load_config_cb, load_sort_, "activate",Void,(),false,(handles,r,s))
id = signal_connect(load_backup_cb, load_backup_, "activate",Void,(),false,(handles,r,s))

#=
Sort Slider Callbacks
=#

id = signal_connect(template_slider, slider_sort, "value-changed", Void, (), false, (handles,))
id = signal_connect(slider_release_template,slider_sort,"button-release-event",Void,(Ptr{Gtk.GdkEventButton},),false,(handles,))

#=
Gain Callbacks
=#

id = signal_connect(sb2_cb,sb2, "value-changed",Void,(),false,(handles,r))
id = signal_connect(SpikeSorting.gain_check_cb,gain_checkbox, "clicked", Void,(),false,(sc_widgets,))

#=
Threshold callbacks
=#

add_thres_cb(sc_widgets)

#=
Pause, Restore, and Clear Callbacks
=#

id = signal_connect(SpikeSorting.pause_cb,button_pause,"toggled",Void,(),false,(handles.sc,))
id = signal_connect(SpikeSorting.clear_button_cb,button_clear,"clicked",Void,(),false,(handles.sc,))
id = signal_connect(SpikeSorting.restore_button_cb,button_restore,"clicked",Void,(),false,(handles.sc,))

#=
Rubberband, Draw, Stretch
=#

for i=1:3
    id = signal_connect(SpikeSorting.pause_state_cb,button_rb[i],"clicked",Void,(),false,(handles.sc,i))
end

#Event Viewer Callbacks
add_event_callbacks(event_handles,popup_event_none,handles)

#=
Radiobutton Callback
=#

for i=1:5
    id = signal_connect(rb1_cb,rbs[i],"clicked",Void,(),false,(handles,i))
end

for i=1:8
    id = signal_connect(rb2_cb,rbs2[i],"clicked",Void,(),false,(handles,i))
end

#=
Save Preferences Callbacks
=#

signal_connect(saving_pref_cb,saving_pref_,"activate",Void,(),false,(handles,r))
add_save_callbacks(save_widgets,handles,r,fpga)

#Soft Scope Callbacks
add_scope_callbacks(scope_v_handles,scope_t_handles,scope_thres_handles,scope_signal_handles,handles)

#Spectrogram Callbacks
add_spect_callbacks(spect_f_handles,spect_w_handles,handles)

#Filtering
signal_connect(band_adj_cb, op_band, "activate",Void,(),false,(handles,r))
add_filter_callbacks(band_widgets,handles,r,fpga)

#Reference
signal_connect(ref_cb, define_ref_, "activate",Void,(),false,(handles,r))
add_reference_callbacks(ref_widgets,handles,r,fpga)

#Parameter Table
signal_connect(table_cb, define_params, "activate",Void,(),false,(handles,r))
add_parameter_callbacks(table_widgets,handles,r,fpga)

signal_connect(ttl_cb,define_ttls_,"activate",Void,(),false,(handles,r))

#=
Sortview Callbacks
=#

id = signal_connect(sv_open_cb, sv_open, "activate",Void,(),false,(handles,))

signal_connect(sortview_handles.win, :delete_event) do widget, event
    visible(sortview_handles.win, false)
    true
end

#=
Backup
=#
f=open(string(r.save.backup,"enabled.bin"),"w")
write(f,handles.enabled)
close(f)

resize!(handles.win,1200,800)

handles
end

#Drawing
function run_cb(widgetptr::Ptr,user_data::Tuple{Gui_Handles,RHD2000,DArray{T,1,Array{T,1}},Task,DArray{I,1,Array{I,1}}}) where {T<:Sorting,I<:IC}

    widget = convert(ToggleButton, widgetptr)

    @async if getproperty(widget,:active,Bool)==true

        #unpack tuple
        han, rhd, s,task,fpga = user_data

	if rhd.debug.state==false
            map(runBoard,fpga)
        end
        while getproperty(widget,:active,Bool)==true
           main_loop_par(rhd,han,s,task,fpga)
        end
    end
    nothing
end

function run_cb(widgetptr::Ptr,user_data::Tuple{Gui_Handles,R,Array{S,1},T,Array{I,1}}) where {R<:RHD2000,S<:Sorting,T<:Task,I<:IC}

    widget = convert(ToggleButton, widgetptr)

    @async if getproperty(widget,:active,Bool)==true

        #unpack tuple
        han, rhd, s,task,fpga = user_data

	if rhd.debug.state==false
            map(runBoard,fpga)
        end
        while getproperty(widget,:active,Bool)==true
           main_loop_s(rhd,han,s,task,fpga)
        end
    end

    nothing
end

function main_loop_s(rhd::RHD2000,han::Gui_Handles,s::Array{T,1},task::Task,fpga::Array{I,1}) where {T<:Sorting, I<:IC}
    if rhd.debug.state==false
        myread=readDataBlocks(rhd,1,s,fpga)
    else
        myread=readDataBlocks(rhd,s)
    end
    main_loop(rhd,han,s,task,myread,fpga)
end

function main_loop_par(rhd::RHD2000,han::Gui_Handles,s::DArray{T,1,Array{T,1}},task::Task,fpga::DArray{I,1,Array{I,1}}) where {T<:Sorting, I<:IC}
    if rhd.debug.state==false
        if rhd.cal<3
            calibrate_parallel(fpga,s,rhd.v,rhd.buf,rhd.nums,rhd.time,rhd.cal)
        else
            onlinesort_parallel(fpga,s,rhd.v,rhd.buf,rhd.nums,rhd.time)
        end
        cal_update(rhd)
	myread=true
    else
        myread=readDataBlocks(rhd,s)
    end
    main_loop(rhd,han,s,task,myread,fpga)
end

function main_loop(rhd::RHD2000,han::Gui_Handles,s,task::Task,myread::Bool,fpga)

    #process and output (e.g. kalman, spike triggered stim calc, etc)
    do_task(task,rhd,myread,han,fpga)

    #plot spikes
    if myread
	if han.num16>0

            #top right
            if han.c_right_top==1
                draw_spike16(rhd,han)
            elseif han.c_right_top==2
                draw_spike32(rhd,han)
            elseif han.c_right_top==3
		draw_spike64(rhd,han)
            elseif han.c_right_top==4
            elseif han.c_right_top==5
            else
            end

            #bottom right
            if han.c_right_bottom==1
	        plot_events(fpga,han,han.draws)
            elseif han.c_right_bottom==2
                draw_raster16(rhd,han)
            elseif han.c_right_bottom==3
                draw_raster32(rhd,han)
            elseif han.c_right_bottom==4
                draw_scope(rhd,han)
            elseif han.c_right_bottom==5
		        draw_spike64(rhd,han)
            elseif han.c_right_bottom==6
                draw_raster64(rhd,han)
            elseif han.c_right_bottom==7
                draw_spectrogram(rhd,han)
            else
            end

            update_time(rhd,han)
	        reveal(han.c)

            if han.spike_changed
                new_single_channel(han,rhd,s,fpga)

                #Clear Sortview
            end
            if han.buf.c_changed

                clus = han.buf.selected_clus

                if clus>0

                    #Find Cluster characteristics from selected waveforms
                    (mymean,mystd)=SpikeSorting.make_cluster(han.buf.spikes,han.buf.mask,han.buf.ind,.!han.buf.selected)

                    #Apply cluster characteristics to handles cluster
                    SpikeSorting.change_cluster(han.sc.temp,mymean,mystd,clus)
                    setproperty!(han.sc.adj_sort, :value, 1.0)

                    if (han.buf.count>0)&(han.sc.pause)
                        SpikeSorting.template_cluster(han.sc,clus,mymean,mystd[:,2],mystd[:,1],1.0)
                    end

                end

                #Send cluster information to sorting
                SpikeSorting.send_clus(s,han.sc)
                backup_clus(han.sc.temp,han.sc.spike,rhd.save.backup)
                han.buf.c_changed=false
            end
            if han.buf.replot
                SpikeSorting.replot_all_spikes(han.sc)
                if visible(han.sortview_widgets.win)
                    if !han.sc.pause

                    else
                        #Refresh Screen
                        SpikeSorting.replot_sort(han.sortview_widgets)
                    end
                end
                han.buf.replot=false
            end
	    if (han.num>0)&(!han.sc.pause)
		   draw_spike(rhd,han)
	    end
            if han.sc.thres_changed
                thres_changed(han,s,fpga,rhd.save.backup)
            end
            if han.sc.show_thres
                SpikeSorting.plot_thres(han.sc)
            end
            if han.sc.rb_active
                SpikeSorting.draw_rb(han.sc)
            end
            draw_c3(rhd,han)
	end
	han.draws+=1
	if han.draws>500
	    han.draws=0
            if han.num16>0
	        clear_c(han)
            end
            if (!han.sc.hold)&(!han.sc.pause)
	            SpikeSorting.clear_c2(han.sc.c2,han.sc.spike)
                han.sc.ctx2=Gtk.getgc(han.sc.c2)
                han.sc.ctx2s=copy(han.sc.ctx2)
            end
            SpikeSorting.clear_c3(han.sc.c3,han.sc.spike)
            #Sort Button
            if han.sc.sort_cb
                SpikeSorting.draw_templates(han.sc)
            end
	end
        reveal(han.sc.c2)
        reveal(han.sc.c3)
	#write to disk, clear buffers
        queueToFile(rhd,task,fpga)
    end
    sleep(.00001)
    #sleep(.02) #debug
    nothing
end

function update_c1(widget::Ptr,user_data::Tuple{Gui_Handles})

    han, = user_data
    han.num16=getproperty(han.adj,:value,Int64) # 16 channels

    if han.num16>0
        clear_c(han)
        han.spike_changed=true
    end
    nothing
end

update_c2_cb(w::Ptr,d::Tuple{Gui_Handles})=update_c2(d[1])

function update_c2(han::Gui_Handles)

    han.num=getproperty(han.adj2, :value, Int64) # primary display

    if han.num16>0

        old_spike=rem(han.sc.spike-1,han.chan_per_display)+1
        han.spike_changed=true

        #Highlight channel
        highlight_channel(han,old_spike)
    end
    nothing
end

function get_cluster(han::Gui_Handles,s::Array{T,1}) where T<:Sorting
    han.sc.temp=deepcopy(s[han.sc.spike].c)
end

function get_cluster(han::Gui_Handles,s::DArray{T,1,Array{T,1}}) where T<:Sorting
    (nn,mycore)=get_thres_id(s,han.sc.spike)

    han.sc.temp=remotecall_fetch(((x,ss)->localpart(x)[ss].c),mycore,s,nn)
end

function backup_clus(myclus,chan,backup)

    f=open(string(backup,"/cluster/",chan,".bin"),"w")
    for i in fieldnames(typeof(myclus))
        write(f,getfield(myclus,i))
    end
    close(f)

    nothing
end

function init_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000,Task,Array{I,1}}) where I<:IC

    han, rhd,task,fpga = user_data
    if !rhd.initialized
        init_board!(rhd,fpga)
        init_task(task,rhd,han,fpga)
        rhd.initialized=true
    end

    nothing
end

function init_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000,Task,DArray{I,1,Array{I,1}}}) where I<:IC

    han, rhd,task,fpga = user_data
    if !rhd.initialized
        init_board!(rhd,fpga)
        init_task(task,rhd,han,fpga)
        rhd.initialized=true
    end

    nothing
end

function record_cb(widgetptr::Ptr, user_data::Tuple{Gui_Handles,R}) where R<:RHD2000

    han, rhd = user_data

    widget = convert(ToggleButton, widgetptr)

    if getproperty(widget,:active,Bool)
        rhd.save.record_mode = true
    else
        sleep(1.0)
        rhd.save.record_mode = false
    end
    nothing
end

#=
SortView Callbacks
=#

function sv_open_cb(widget::Ptr,user_data::Tuple{Gui_Handles})

    han, = user_data
    visible(han.sortview_widgets.win,true)
    nothing
end

#=
Right Canvas Callbacks
=#

function c_popup_select(widget::Ptr,param_tuple,user_data::Tuple{Gui_Handles})

    han, = user_data
    event = unsafe_load(param_tuple)

    ctx=Gtk.getgc(han.c)
    myheight=height(ctx)

    han.mim=(event.x,event.y)

    if event.y<(myheight-300) #top

        if han.c_right_top==1 #disable enable 16

            if event.button == 1 #left click

                (inmulti,channel_num)=check_multi(han,4,4,16,event.x,event.y)

                if inmulti
                    Gtk.GAccessor.value(han.adj2,channel_num)
                    update_c2(han)
                end

            elseif event.button == 3 #right click

                popup(han.popup_ed,event)
            end
        elseif han.c_right_top==2 #disable enable 32

            if event.button == 1 #left click

                (inmulti,channel_num)=check_multi(han,6,6,32,event.x,event.y)
                if inmulti
                    Gtk.GAccessor.value(han.adj2,channel_num)
                    update_c2(han)
                end

            elseif event.button == 3 #right click

                popup(han.popup_ed,event)
            end

        elseif han.c_right_top==3 #disable enable 64

            if event.button == 1 #left click

                (inmulti,channel_num)=check_multi(han,11,6,64,event.x,event.y)
                if inmulti
                    Gtk.GAccessor.value(han.adj2,channel_num)
                    update_c2(han)
                end

            elseif event.button == 3 #right click

                popup(han.popup_ed,event)
            end

        elseif han.c_right_top==4 #64 channel raster - nothing
        else
        end
    else #bottom

        if han.c_right_bottom==1 #event select

            if event.button ==3 #right click
                popup(han.popup_event,event)
            end

        elseif han.c_right_bottom==2 #16 raster - nothing
        elseif han.c_right_bottom==3 #32 raster - nothing
        elseif han.c_right_bottom==4 #select scope channels

            if event.button == 3 #right click
                popup(han.popup_scope,event)
            end

        elseif han.c_right_bottom==5 #disable enable 64

            if event.button == 1 #left click

                (inmulti,channel_num)=check_multi(han,11,6,64,event.x,event.y)
                if inmulti
                    Gtk.GAccessor.value(han.adj2,channel_num)
                    update_c2(han)
                end

            elseif event.button == 3 #right click

                popup(han.popup_ed,event)
            end

        elseif han.c_right_bottom==6 #64 channel raster - nothing
        elseif han.c_right_bottom==7 #spectrogram

            if event.button == 3 #right click
                popup(han.popup_spect,event)
            end

        else

        end
    end
    nothing
end

function check_multi(han::Gui_Handles,n_row::Int64,n_col::Int64,num_chan::Int64,x,y)

    (xbounds,ybounds)=get_multi_bounds(han,n_col,n_row,num_chan)

    count=1
    inmulti=false
    for j=2:length(xbounds), i=2:length(ybounds)
        if (x<xbounds[j])&(y<ybounds[i])
            inmulti=true
            break
        end
        count+=1
    end

    (inmulti,count)
end

function get_multi_bounds(han::Gui_Handles,n_col,n_row,num_chan)

    ctx=Gtk.getgc(han.c)

    mywidth=width(ctx)
    if num_chan<64
        myheight=height(ctx)-300.0
    else
        myheight=height(ctx)
    end

    xbounds=linspace(0.0,mywidth,n_col+1)
    ybounds=linspace(0.0,myheight,n_row+1)

    (xbounds,ybounds)
end

function get_multi_dims(han::Gui_Handles,n_col,n_row,num_chan,spike)

    (xbounds,ybounds)=get_multi_bounds(han,n_col,n_row,num_chan)

    y=rem(rem(spike-1,num_chan),n_row)+1
    x=div(rem(spike-1,num_chan),n_row)+1

    (xbounds[x],xbounds[x+1],ybounds[y],ybounds[y+1])
end

#(handles,r,fpga))
function close_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000,Array{I,1}}) where I<:IC

    han, rhd, fpgas = user_data

    #If we were in run loop, turn off

    for i=1:length(fpgas)
        resetBoard(fpgas[i])
        if (OPEN_EPHYS)
            enableBoardLeds(fpgas[i],false)
        end
    end

    println("Bye!")

    nothing
end


text:
module PlutoUI

import Base: show, get
import Markdown: htmlesc, withtag

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))

include("./Builtins.jl")
include("./Clock.jl")
include("./Resource.jl")
include("./Terminal.jl")
include("./TableOfContents.jl")
include("./RangeSlider.jl")
include("./DisplayTricks.jl")

end

ext: jl


text:
type PathGenerator{RSG <: AbstractRandomSequenceGenerator, S <: StochasticProcess1D}
  brownianBridge::Bool
  generator::RSG
  dimension::Int
  timeGrid::TimeGrid
  process::S
  nextSample::Sample
  temp::Vector{Float64}
  bb::BrownianBridge
end

function PathGenerator(process::StochasticProcess, tg::TimeGrid, generator::AbstractRandomSequenceGenerator, brownianBridge::Bool)
  generator.dimension == length(tg.times) - 1 || error("wrong dimensions")

  return PathGenerator(brownianBridge, generator, generator.dimension, tg, process, Sample(Path(tg), 1.0), zeros(generator.dimension), BrownianBridge(tg))
end

function PathGenerator(process::StochasticProcess, len::Float64, timeSteps::Int, generator::AbstractRandomSequenceGenerator, brownianBridge::Bool)
  dims = generator.dimension
  timeSteps == dims || error("sequence generator dimensionality error")

  tg = TimeGrid(len, timeSteps)

  return PathGenerator(brownianBridge, generator, dims, tg, process, Sample(Path(tg), 1.0), zeros(dims), BrownianBridge(tg))
end

get_next!(pg::PathGenerator) = get_next!(pg, false)
get_antithetic!(pg::PathGenerator) = get_next!(pg, true)

function get_next!(pg::PathGenerator, antithetic::Bool)
  sequenceVals, sequenceWeight = antithetic ? last_sequence(pg.generator) : next_sequence!(pg.generator)

  # TODO Brownian Bridge
  pg.temp = copy(sequenceVals)

  pg.nextSample.weight = sequenceWeight
  # path = pg.nextSample.value # TODO check to see if this is just a pass by reference

  pg.nextSample.value[1] = get_x0(pg.process)

  for i = 2:length(pg.nextSample.value)
    t = pg.timeGrid[i-1]
    dt = pg.timeGrid.dt[i - 1]
    pg.nextSample.value[i] = evolve(pg.process, t, pg.nextSample.value[i-1], dt, antithetic ? -pg.temp[i-1] : pg.temp[i-1])
  end

  return pg.nextSample
end


text:
# ---
# title: 1478. Allocate Mailboxes
# id: problem1478
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/allocate-mailboxes/description/>
# hidden: true
# ---
#
# Given the array `houses` and an integer `k`. where `houses[i]` is the location
# of the ith house along a street, your task is to allocate `k` mailboxes in the
# street.
#
# Return the **minimum** total distance between each house and its nearest
# mailbox.
#
# The answer is guaranteed to fit in a 32-bit signed integer.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/05/07/sample_11_1816.png)
#
#
#
#     Input: houses = [1,4,8,10,20], k = 3
#     Output: 5
#     Explanation: Allocate mailboxes in position 3, 9 and 20.
#     Minimum total distance from each houses to nearest mailboxes is |3-1| + |4-3| + |9-8| + |10-9| + |20-20| = 5
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/07/sample_2_1816.png)**
#
#
#
#     Input: houses = [2,3,5,12,18], k = 2
#     Output: 9
#     Explanation: Allocate mailboxes in position 3 and 14.
#     Minimum total distance from each houses to nearest mailboxes is |2-3| + |3-3| + |5-3| + |12-14| + |18-14| = 9.
#
#
# **Example 3:**
#
#
#
#     Input: houses = [7,4,6,1], k = 1
#     Output: 8
#
#
# **Example 4:**
#
#
#
#     Input: houses = [3,6,14,10], k = 4
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `n == houses.length`
#   * `1 <= n <= 100`
#   * `1 <= houses[i] <= 10^4`
#   * `1 <= k <= n`
#   * Array `houses` contain unique integers.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
diffper(n,xspan)

Construct 2nd-order differentiation matrices for functions with
periodic end conditions, using `n` unique nodes in the interval
`xspan`. Return a vector of nodes and the  matrices for the first
and second derivatives.
"""
function diffper(n,xspan)
    a,b = xspan
    h = (b-a)/n
    x = @. a + h*(0:n-1)   # nodes, omitting the repeated data

    # Construct Dx by diagonals, then correct the corners.
    dp = fill(0.5/h,n-1)        # superdiagonal
    dm = fill(-0.5/h,n-1)       # subdiagonal
    Dx = diagm(-1=>dm,1=>dp)
    Dx[1,n] = -1/(2*h)
    Dx[n,1] = 1/(2*h)

    # Construct Dxx by diagonals, then correct the corners.
    d0 =  fill(-2/h^2,n)        # main diagonal
    dp =  ones(n-1)/h^2         # superdiagonal and subdiagonal
    Dxx = diagm(-1=>dp,0=>d0,1=>dp)
    Dxx[1,n] = 1/(h^2)
    Dxx[n,1] = 1/(h^2)

    return x,Dx,Dxx
end

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
#############################################################################
# SCS.jl
# Wrapper around the SCS solver https://github.com/cvxgrp/scs
#############################################################################
# SCSSolverInterface.jl
# MathProgBase.jl interface for the SCS.jl solver wrapper
#############################################################################

using LinearAlgebra: dot
using MathProgBase.SolverInterface

import MathProgBase.SolverInterface: ConicModel, LinearQuadraticModel,
    getdual, getobjval, getsolution, getsolvetime, getvardual, loadproblem!,
    numconstr, numvar, optimize!, setbvec!, setwarmstart!, status,
    supportedcones

import Base.convert

# TODO: don't add to Base.convert!
function convert(x::Type{Int}, y::UnitRange{Int})
    if length(y) == 1
        return y[1]
    else
        error("convert` has no method matching convert(::Type{Int}, ::UnitRange{Int})")
    end
end
#############################################################################
# Define the MPB Solver and Model objects
export SCSSolver
struct SCSSolver <: AbstractMathProgSolver
    options
end
SCSSolver(;kwargs...) = SCSSolver(kwargs)

mutable struct SCSMathProgModel <: AbstractConicModel
    m::Int                            # Number of constraints
    n::Int                            # Number of variables
    input_numconstr::Int
    input_numvar::Int
    A::SparseMatrixCSC{Float64,Int}   # The A matrix (equalities)
    b::Vector{Float64}                # RHS
    c::Vector{Float64}                # The objective coeffs (always min)
    f::Int                            # number of zero cones
    l::Int                            # number of linear cones { x | x >= 0}
    q::Vector{Int}                    # Array of SOC sizes
    s::Vector{Int}                    # Array of SDP sizes
    ep::Int                           # Number of primal exponential cones
    ed::Int                           # Number of dual exponential cones
    orig_sense::Symbol                # Original objective sense
    # Post-solve
    solve_stat::Symbol
    solve_time::Float64
    obj_val::Float64
    primal_sol::Vector{Float64}
    dual_sol::Vector{Float64}
    slack::Vector{Float64}
    row_map_ind::Vector{Int}
    row_map_type::Vector{Symbol}
    col_map_ind::Vector{Int}          # map from MPB variables to rows
    col_map_type::Vector{Symbol}
    options
end

SCSMathProgModel(;kwargs...) = SCSMathProgModel(0, 0, 0, 0, spzeros(0, 0), Int[], Int[],
                                      0, 0, Int[], Int[], 0, 0,
                                      :Min, :NotSolved, 0.0, 0.0, Float64[], Float64[],
                                      Float64[], Int[], Symbol[],
                                      Int[], Symbol[], kwargs)

#############################################################################
# Begin implementation of the MPB low-level interface
# Implements
# - ConicModel
# - loadproblem!
# - optimize!
# - status
# - numvar
# - numconstr
# http://mathprogbasejl.readthedocs.org/en/latest/solverinterface.html

ConicModel(s::SCSSolver) = SCSMathProgModel(;s.options...)
LinearQuadraticModel(s::SCSSolver) = ConicToLPQPBridge(ConicModel(s))

#=
function setsense!(m::SCSMathProgModel, sns::Symbol)
    if m.orig_sense != sns
        sns == :Min || sns == :Max || error("Unrecognized sense $sns")
        m.orig_sense = sns
        m.c *= -1
    end
    nothing
end
=#

# TODO needs to be updated for newest constants
const status_map = Dict{Int, Symbol}(
    1 => :Optimal,
    -2 => :Infeasible,
    -1 => :Unbounded,
    -3 => :Indeterminate,
    -4 => :Error
)

function optimize!(m::SCSMathProgModel)
    linear_solver, options = sanitize_SCS_options(m.options)
    t = time()
    solution = SCS_solve(linear_solver, m.m, m.n, m.A, m.b, m.c, m.f, m.l, m.q,
                         m.s, m.ep, m.ed, Float64[],
                         m.primal_sol, m.dual_sol, m.slack; options...)
    m.solve_time = time() - t

    m.solve_stat = get(status_map, solution.ret_val, :UnknownError)
    m.primal_sol = solution.x

    m.dual_sol = solution.y

    # TODO: Get the right slack variables in the right order
    m.slack = solution.s

    m.obj_val = dot(m.c, m.primal_sol) * (m.orig_sense == :Max ? -1 : +1)
end

status(m::SCSMathProgModel) = m.solve_stat
getobjval(m::SCSMathProgModel) = m.obj_val
getsolution(m::SCSMathProgModel) = copy(m.primal_sol)

function invertsdconesize(p)
    return (sqrt(8*p+1) - 1) / 2
end

function isintegertol(n)
    return abs(n - convert(Int, n)) < 1e-16
end

#############################################################################
# Begin implementation of the MPB conic interface
# Implements
# - loadconicproblem!
# - supportedcones
# http://mathprogbasejl.readthedocs.org/en/latest/conic.html

function orderconesforscs(A_in, b_in, c_cones, v_cones)
    # Order the cones as:
    # Free, Zero, NonNeg (NonPos are converted), SOC, SDP, ExpPrimal, ExpDual
    #
    # Returns:
    # - scs_A (A ordered as needed), b
    # - num_free, num_zero, num_linear, soc_sizes, soc_sizes, sqrt_sdp_sizes,
    # - sqrt_sdp_size, num_expprimal, num_expdual

    m, n = size(A_in)
    A_in_t = A_in'
    A_t = spzeros(n,0)
    b = zeros(0)
    row_map_ind = zeros(Int, length(b_in))
    row_map_type = Array{Symbol}(undef, length(b_in))
    col_map_ind = zeros(Int, n)
    col_map_type = Array{Symbol}(undef, n)

    # First, count the total number of variables
    num_vars = 0
    for (cone, idxs) in v_cones
        col_map_type[idxs] .= cone
        num_vars += length(idxs)
    end
    @assert num_vars == n

    num_free = 0
    zeroidx = Int[]
    nonnegidx = Int[]
    nonposidx = Int[]
    socidx = Int[]
    soc_sizes = Int[]
    new_c_cones = Any[]

    for (cone, idxs) in c_cones
        if cone == :Free
            error("Why are you passing in a free constraint?")
        end
        # merge some cones for efficiency
        if cone == :Zero
            append!(zeroidx, idxs)
        elseif cone == :NonNeg
            append!(nonnegidx, idxs)
        elseif cone == :NonPos
            append!(nonposidx, idxs)
        elseif cone == :SOC
            append!(socidx, idxs)
            push!(soc_sizes, length(idxs))
        else
            push!(new_c_cones, (cone,idxs))
        end
    end
    length(zeroidx) > 0 && push!(new_c_cones, (:Zero, zeroidx))
    length(nonnegidx) > 0 && push!(new_c_cones, (:NonNeg, nonnegidx))
    length(nonposidx) > 0 && push!(new_c_cones, (:NonPos, nonposidx))
    length(socidx) > 0 && push!(new_c_cones, (:SOC, socidx))

    for (cone, idxs) in v_cones
        if cone == :Free
            num_free += length(idxs)
        end
    end

    num_zero = 0
    for (cone, idxs) in new_c_cones
        if cone == :Zero
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            A_t = [A_t A_in_t[:,idxs]]
            b = [b; b_in[idxs,:]]
            num_zero += length(idxs)
        end
    end
    for (cone, idxs) in v_cones
        if cone == :Zero
            nidx = length(idxs)
            col_map_ind[idxs] = (length(b)+1):(length(b)+nidx)
            A_t = [A_t sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]
            num_zero += nidx
        end
    end

    num_lin = 0
    for (cone, idxs) in new_c_cones
        if cone == :NonNeg
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            A_t = [A_t A_in_t[:,idxs]]
            b = [b; b_in[idxs,:]]
            num_lin += length(idxs)
        elseif cone == :NonPos
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            A_t = [A_t -A_in_t[:,idxs]]
            b = [b; -b_in[idxs,:]]
            num_lin += length(idxs)
        end
    end
    for (cone, idxs) in v_cones
        nidx = length(idxs)
        col_map_ind[idxs] = (length(b)+1):(length(b)+nidx)
        if cone == :NonNeg
            A_t = [A_t -sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]
            num_lin += nidx
        elseif cone == :NonPos
            A_t = [A_t sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]
            num_lin += nidx
        end
    end

    for (cone, idxs) in new_c_cones
        if cone == :SOC
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            A_t = [A_t A_in_t[:,idxs]]
            b = [b; b_in[idxs,:]]
        end
    end
    for (cone, idxs) in v_cones
        if cone == :SOC
            nidx = length(idxs)
            col_map_ind[idxs] = (length(b)+1):(length(b)+nidx)
            A_t = [A_t -sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]
            push!(soc_sizes, nidx)
        end
    end

    sqrt_sdp_sizes = Int[]
    for (cone, idxs) in new_c_cones
        if cone == :SDP
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            A_t = [A_t A_in_t[:,idxs]]
            b = [b; b_in[idxs,:]]
            # n must be a square integer
            n = length(idxs)
            isintegertol(invertsdconesize(n)) || error("number of SDP variables must be n*(n+1)/2")
            sqrt_n = convert(Int, invertsdconesize(n));
            push!(sqrt_sdp_sizes, sqrt_n)
        end
    end
    for (cone, idxs) in v_cones
        if cone == :SDP
            nidx = length(idxs)
            col_map_ind[idxs] = (length(b)+1):(length(b)+nidx)
            A_t = [A_t -sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]
             # n must be a square integer
            isintegertol(invertsdconesize(nidx)) || error("number of SDP variables must be n*(n+1)/2")
            sqrt_n = convert(Int, invertsdconesize(nidx));
            push!(sqrt_sdp_sizes, sqrt_n)
        end
    end

    num_expprimal = 0
    for (cone, idxs) in new_c_cones
        if cone == :ExpPrimal
            length(idxs) % 3 == 0 ||
                error("Number of ExpPrimal variables must be a multiple of 3")
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            A_t = [A_t A_in_t[:,idxs]]
            b = [b; b_in[idxs,:]]

            num_expprimal += div(length(idxs), 3)
        end
    end
    for (cone, idxs) in v_cones
        if cone == :ExpPrimal
            length(idxs) % 3 == 0 ||
                error("Number of ExpPrimal variables must be a multiple of 3")
            nidx = length(idxs)
            col_map_ind[idxs] = (length(b)+1):(length(b)+nidx)
            A_t = [A_t -sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]

            num_expprimal += div(length(idxs), 3)
        end
    end

    num_expdual = 0
    for (cone, idxs) in new_c_cones
        if cone == :ExpDual
            row_map_ind[idxs] = size(A_t, 2)+1:size(A_t, 2)+length(idxs)
            row_map_type[idxs] = [cone for i in 1:length(idxs)]

            length(idxs) % 3 == 0 ||
                error("Number of ExpDual variables must be a multiple of 3")
            A_t = [A_t A_in_t[:,idxs]]
            b = [b; b_in[idxs,:]]

            num_expdual += div(length(idxs), 3)
        end
    end
    for (cone, idxs) in v_cones
        if cone == :ExpDual
            length(idxs) % 3 == 0 ||
                error("Number of ExpDual variables must be a multiple of 3")
            nidx = length(idxs)
            col_map_ind[idxs] = (length(b)+1):(length(b)+nidx)
            A_t = [A_t -sparse(idxs, 1:nidx, ones(nidx), num_vars, nidx)]
            b = [b; zeros(nidx)]

            num_expdual += div(length(idxs), 3)
        end
    end

    return A_t', b, num_free, num_zero, num_lin, soc_sizes,
           sqrt_sdp_sizes, num_expprimal, num_expdual, col_map_ind, col_map_type, row_map_ind, row_map_type
end

loadproblem!(model::SCSMathProgModel, c, A, b, constr_cones, var_cones) =
    loadproblem!(model, c, sparse(A), b, constr_cones, var_cones)

function loadproblem!(model::SCSMathProgModel, c, A::SparseMatrixCSC, b, constr_cones, var_cones)
    # TODO: We should support SOCRotated
    bad_cones = [:SOCRotated]
    for cone_vars in constr_cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) not supported")
    end
    for cone_vars in var_cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) not supported")
    end

    # Convert idxs to an array
    c_cones = [(cone, [idxs...]) for (cone, idxs) in constr_cones]
    v_cones = [(cone, [idxs...]) for (cone, idxs) in var_cones]

    scs_A, scs_b, num_free, f, l, q, s, ep, ed, col_map_ind, col_map_type, row_map_ind, row_map_type =
        orderconesforscs(A, b, c_cones, v_cones)

    m, n = size(scs_A)

    model.n             = n
    model.m             = m # + rows_G
    model.A             = scs_A
    model.b             = scs_b[:]
    model.c             = c[:]
    model.q             = q
    model.s             = s
    model.ep            = ep
    model.ed            = ed
    model.orig_sense    = :Min
    model.f             = f
    model.l             = l
    model.col_map_ind   = col_map_ind
    model.col_map_type  = col_map_type
    model.row_map_ind   = row_map_ind
    model.row_map_type  = row_map_type
    model.input_numconstr = size(A,1)
    model.input_numvar    = size(A,2)

    return model
end

numvar(model::SCSMathProgModel) = model.input_numvar
numconstr(model::SCSMathProgModel) = model.input_numconstr

supportedcones(s::SCSSolver) = [:Free, :Zero, :NonNeg, :NonPos, :SOC, :SDP, :ExpPrimal, :ExpDual]

function getdual(m::SCSMathProgModel)
    dual = m.dual_sol[m.row_map_ind]
    # flip sign for NonPos since it's treated as NonNeg by SCS
    for i in 1:length(m.row_map_type)
        if m.row_map_type[i] == :NonPos
            dual[i] = -dual[i]
        end
    end
    return dual
end

function getvardual(m::SCSMathProgModel)
    dual = zeros(length(m.col_map_ind))
    for i in 1:length(m.col_map_type)
        if m.col_map_type[i] == :Free
            continue # dual is zero
        elseif m.col_map_type[i] == :NonPos
            # flip sign for NonPos since it's treated as NonNeg by SCS
            dual[i] = -m.dual_sol[m.col_map_ind[i]]
        else
            dual[i] = m.dual_sol[m.col_map_ind[i]]
        end
    end
    return dual
end

function addoption!(m::SCSMathProgModel, option::Symbol, value)
    nt = NamedTuple{(option,), Tuple{typeof(value)}}((value,))
    m.options = pairs(merge(m.options.data, nt))
    return m
end

# warmstart
# kwargs can be `primal_sol`, `dual_sol`, and `slack`
function setwarmstart!(m::SCSMathProgModel, primal_sol; kwargs...)
    addoption!(m, :warm_start, true)
    m.primal_sol = primal_sol
    for (k,v) in kwargs
        setfield!(m, k, v)
    end

    # check sizes to prevent segfaults
    nconstr, nvar = size(m.A)
    length(m.primal_sol) == nvar || (m.primal_sol = zeros(nvar))
    length(m.dual_sol) == nconstr || (m.dual_sol = zeros(nconstr))
    length(m.slack) == nconstr || (m.slack = zeros(nconstr))
    m
end

function setbvec!(m::SCSMathProgModel, b::Vector{Float64})
    m.b[m.row_map_ind] = b
end

getsolvetime(m::SCSMathProgModel) = m.solve_time


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

using LaTeXStrings
using Plots
pyplot()

# DeltaV max
println("Delta V max, ", sqrt(DeltaV2(0.1 * AU, 0.0, pars, mol)) * 1e-5)

# For a set of velocities, make channel maps to test whether we can mask appropriately.

# V4046 Sgr has: -3 to + 3 kms, in 0.4 km/s chunks.
nv = 22
vs = linspace(-4.0, 4.0, nv)

npix = 128
xs = linspace(-500 * AU, 500 * AU, npix)
ys = linspace(-500 * AU, 500 * AU, npix)

arr = Array{Bool}(npix, npix, nv)

# Fill the boolean mask
for k=1:nv
    for j=1:npix
        for i=1:npix
            arr[j,i,k] = verify_pixel(xs[i], ys[j], pars, vs[k], 0.5, 400*AU)
        end
    end
end

println("Fraction traced: ",sum(arr)/(npix * npix * nv))

function plot_channel_maps(arr)
    parr = [heatmap(arr[:,:,k]) for k=1:nv]
    p = plot(parr..., layout=(2, 11), size=(11 * 250,500))
    savefig(p, "mask.png")
end

plot_channel_maps(arr)


text:
@testset "Vertical interpolation" begin
    D = ClimArray([1.0:1.0:11.0 2.0:1.0:12.0 3.0:1.0:13.0], (Hei(0.0:2000.0:20000.0), Ti(1:3)))
    pressure_levels = [950.0,850.0,650.0,350.0,250.0,150.0] .* 100.0
    D_pre = interpolate_height2pressure(D, pressure_levels,extrapolation_bc=NaN)
    D_back = interpolate_pressure2height(D_pre, Vector(dims(D,Hei).val),extrapolation_bc=Line())
    my_dim = Dim{:My_Dim}
    E = ClimArray([1.0:1.0:10.0 2.0:1.0:11.0 3.0:2.0:21.0], (my_dim(1:10), Ti(1:3)))
    pressure = ClimArray([1000.0:-100.0:100.0 1000.0:-100.0:100.0 1000.0:-100.0:100.0] *100.0, (my_dim(1:10), Ti(1:3)))
    E_pre = interpolation2pressure(E, pressure, pressure_levels; vertical_coord=my_dim, extrapolation_bc=NaN )
    E_pre2 = interpolation2pressure(reverse(E,dims=my_dim), reverse(pressure,dims=my_dim), pressure_levels; vertical_coord=my_dim, extrapolation_bc=NaN, descending = false )

    @test hasdim(D_pre,Pre())
    @test gnv(dims(D_back,Hei)) == gnv(dims(D,Hei)) && gnv(dims(D_back,Ti)) == gnv(dims(D,Ti))
    @test E_pre.data == E_pre2.data
    @test E_pre[Pre(1)] < E_pre[Pre(2)] < E_pre[Pre(3)] < E_pre[Pre(4)] < E_pre[Pre(5)]
    @test D_pre[Pre(1)] < D_pre[Pre(2)] < D_pre[Pre(3)] < D_pre[Pre(4)] < D_pre[Pre(5)]
    @test D_back[Hei(1)] < D_back[Hei(2)] < D_back[Hei(3)] < D_back[Hei(4)] < D_back[Hei(5)]
    @test D[Hei(1)] < D_pre[Pre(3)] < D[Hei(11)]
end


text:
using GLVisualize, GeometryTypes, ModernGL, GLAbstraction, Colors, FileIO
if !isdefined(:runtests)
    window = glscreen()
end

description = """
Demonstrating how wireframes can be displayed.
"""

# loadasset is defined as `loadasset(name) = load(assetphat(name))`
# load comes from FileIO, which automatically detects the file format and loads it
# in general, load will try to load the mesh as is, which may return different
# meshtypes. This is why it's wrapped into GLNormalMesh, to guarantee some type
# stability.
# Formats supported are currently: obj, stl, ply, off and 2DM
msh = GLNormalMesh(loadasset("cat.obj"))
v = vertices(msh)
f = faces(msh)
colors = map(v) do _ # for each vertex
    RGBA(rand(RGB{Float32}), 0.4f0)
end

colored_mesh = GLNormalVertexcolorMesh(
    vertices = v, faces = f, color = colors
)

_view(visualize(colored_mesh), window, camera=:perspective)
_view(visualize(
    msh, :lines, thickness = 1f0,
    color = RGBA(1f0, 1f0, 1f0, 0.8f0)
), window, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
function save_model(model, prefix)
    m = cpu(model)
    w = weights = Tracker.data.(params(m))
    @save "$(prefix)-model.bson" m
    @save "$(prefix)-weights.bson" w
end

function load_model(prefix)
    @load "$(prefix)-model.bson" m
    @load "$(prefix)-weights.bson" w
    Flux.loadparams!(m, w)
    return m
end

max_stars_repo_name: darsnack/rl-experiments

max_issues_repo_name: darsnack/rl-experiments

max_forks_repo_name: darsnack/rl-experiments


text:
include("parallelLayerForProp.jl")

using ProgressMeter

try
    ProgressMeter.ijulia_behavior(:clear)
catch

end
using Random
using LinearAlgebra

###

@doc raw"""
    function chainForProp(
        X::AbstractArray{T,N},
        cLayer::Layer,
        cnt::Integer = -1;
        FCache = Dict{Layer,Dict{Symbol,AbstractArray}}(),
        kwargs...,
    ) where {T,N}

perform the chained forward propagation using recursive calls

# Arguments:

- `X::AbstractArray{T,N}` := input of the input layer

- `cLayer::Layer` := Input Layer

- `cnt::Integer` := an internal counter used to cache the layers was performed not to redo it again

# Returns

- `Cache::Dict{Layer, Dict{Symbol, Array}}` := the output each layer either A, Z or together As Dict of layer to dict of Symbols and Arrays for internal use, it set again the values of Z and A in each layer to be used later in back propagation and add one to the layer forwCount value when pass through it
"""
function chainForProp(
    X::AbstractArray{T,N},
    cLayer::Layer,
    cnt::Integer = -1;
    FCache = Dict{Layer,Dict{Symbol,AbstractArray}}(),
    kwargs...,
) where {T,N}
    if cnt < 0
        cnt = cLayer.forwCount + 1
    end

    if length(cLayer.nextLayers) == 0
        if cLayer.forwCount < cnt
            FCache[cLayer] = layerForProp(cLayer; FCache = FCache, kwargs...)
        end #if cLayer.forwCount < cnt
        return FCache
    elseif isa(cLayer, MILayer) #if typeof(cLayer)==AddLayer
        if all(
            i -> (i.forwCount == cLayer.prevLayer[1].forwCount),
            cLayer.prevLayer,
        )
            FCache[cLayer] = layerForProp(cLayer; FCache = FCache, kwargs...)
            for nextLayer in cLayer.nextLayers
                FCache =
                    chainForProp(X, nextLayer, cnt; FCache = FCache, kwargs...)
            end
        end #if all

        return FCache
    else #if cLayer.prevLayer==nothing
        if cLayer.forwCount < cnt
            if cLayer isa Input
                FCache[cLayer] =
                    layerForProp(cLayer, X; FCache = FCache, kwargs...)
            else
                FCache[cLayer] =
                    layerForProp(cLayer; FCache = FCache, kwargs...)
            end
            for nextLayer in cLayer.nextLayers
                FCache =
                    chainForProp(X, nextLayer, cnt; FCache = FCache, kwargs...)
            end

        end #if cLayer.forwCount < cnt
        return FCache
    end #if cLayer.prevLayer!=nothing

    return FCache

end #function chainForProp!

export chainForProp

@doc raw"""
    predictBatch(model::Model, X::AbstractArray, Y = nothing; kwargs...)

predict Y using the model and the input X and the labels Y

# Inputs

- `model::Model` := the trained model

- `X::AbstractArray` := the input `Array`

- `Y` := the input labels to compare with (optional)

# Output

- a `Tuple` of

    * `Ŷ` := the predicted values
    * `Ŷ_bool` := the predicted labels
    * `"accuracy"` := the accuracy of the predicted labels

"""
function predictBatch(model::Model, X::AbstractArray, Y = nothing; kwargs...)

    kwargs = Dict{Symbol, Any}(kwargs...)
    kwargs[:prediction] = getindex(kwargs, :prediction, default = true)

    FCache = chainForProp(X, model.inLayer; kwargs...)
    Ŷ = FCache[model.outLayer][:A]
    T = eltype(Ŷ)
    outLayer = model.outLayer
    actFun = outLayer.actFun
    lossFun = model.lossFun
    costs = nothing
    if Y != nothing
        costs = cost(eval(:($lossFun)), Ŷ, Y)
    end
    # if isbool(Y)
    return Ŷ, probToValue(eval(:($actFun)), Ŷ; labels = Y)..., costs

end #predict

export predictBatch

###
@doc raw"""
    predict(model::Model, X_In::AbstractArray, Y_In = nothing; kwargs...)

Run the prediction based on the trained `model`

# Arguments

- `model::Model` := the trained `Model` to predict on

- `X_In` := the input `Array`

- `Y_In` := labels (optional) to evaluate the model

## Key-word Arugmets

- `batchSize` := default `32`

- `useProgBar` := (`Bool`) where or not to shoe the prograss bar

# Return

- a `Dict` of:

    * `:YhatValue` := Array of the output of the integer prediction values
    * `:YhatProb` := Array of the output probabilities
    * `:accuracy` := the accuracy of prediction in case `Y_In` is given
"""
function predict(model::Model, X_In::AbstractArray, Y_In = nothing; kwargs...)

    kwargs = Dict{Symbol, Any}(kwargs...)
    batchSize = getindex(kwargs, :batchSize; default = 32)
    printAcc = getindex(kwargs, :printAcc; default = true)
    useProgBar = getindex(kwargs, :useProgBar; default = true)
    GCInt = getindex(kwargs, :GCInt, default = 5)
    noBool = getindex(kwargs, :noBool, default = false)
    kwargs[:prediction] = getindex(kwargs, :prediction, default = true)

    outLayer, lossFun, α = model.outLayer, model.lossFun, model.α
    Costs = []

    nX = ndims(X_In)
    T = eltype(X_In)
    m = size(X_In)[end]
    # c = size(Y_train)[end-1]
    nB = m ÷ batchSize
    axX = axes(X_In)[1:end-1]
    Y = nothing
    if Y_In != nothing
        axY = axes(Y_In)[1:end-1]
        costs =
            Array{AbstractFloat,1}(undef, nB + ((m % batchSize == 0) ? 0 : 1))
    end
    if useProgBar
        p = Progress((m % batchSize != 0 ? nB + 1 : nB), 0.1)
    end

    nY = length(model.outLayer.outputS)
    Ŷ_out =
        Array{AbstractArray{T,nY},1}(undef, nB + ((m % batchSize == 0) ? 0 : 1))
    Ŷ_prob_out =
        Array{AbstractArray{T,nY},1}(undef, nB + ((m % batchSize == 0) ? 0 : 1))
    accuracy =
        Array{AbstractFloat,1}(undef, nB + ((m % batchSize == 0) ? 0 : 1))
    # Threads.@threads for j = 1:nB
    @simd for j = 1:nB
        downInd = (j - 1) * batchSize + 1
        upInd = j * batchSize
        X = view(X_In, axX..., downInd:upInd)
        if Y_In != nothing
            Y = view(Y_In, axY..., downInd:upInd)
        end
        Ŷ_prob_out[j], Ŷ_out[j], accuracy[j], costs[j] =
            predictBatch(model, X, Y; kwargs...)

        if useProgBar
            update!(p, j, showvalues = [("Instances $m", j * batchSize)])
        end #if useProgBar
        # X = Y = nothing
        # Ŷ = nothing
        # if j%GCInt == 0
        #     Base.GC.gc()
        # end
    end

    if m % batchSize != 0
        downInd = (nB) * batchSize + 1
        X = view(X_In, axX..., downInd:m)
        if Y_In != nothing
            Y = view(Y_In, axY..., downInd:m)
        end
        Ŷ_prob_out[nB+1], Ŷ_out[nB+1], accuracy[nB+1], costs[nB+1] =
            predictBatch(model, X, Y; kwargs...)

        # X = Y = nothing
        # Base.GC.gc()
        if useProgBar
            update!(p, nB + 1, showvalues = [("Instances $m", m)])
        end
    end

    # accuracyM = nothing
    # if !(isempty(accuracy))
    accuracyM = mean(filter(x -> x != nothing, accuracy))
    costM = nothing
    if Y_In != nothing
        costM = mean(filter(x -> x != nothing, costs))
    end
    # end
    if noBool
        Ŷ_values = nothing
        Ŷ_prob = nothing
    else
        # Ŷ_values = Array{T,nY}(undef, repeat([0], nY)...)
        Ŷ_values = cat(Ŷ_out... ; dims = nY)
        Ŷ_prob = cat(Ŷ_prob_out...; dims = nY)
    end

    return Dict(
        :YhatValue => Ŷ_values,
        :YhatProb => Ŷ_prob,
        :accuracy => accuracyM,
        :cost => costM,
    )

end #function predict(
# model::Model,
# X_In::AbstractArray,
# Y_In=nothing;
# batchSize = 32,
# printAcc = true,
# useProgBar = false,
# )

export predict

### chainBackProp
include("parallelLayerBackProp.jl")
# include("layerUpdateParams.jl")

@doc raw"""
    function chainBackProp(
        X::AbstractArray{T1,N1},
        Y::AbstractArray{T2,N2},
        model::Model,
        FCache::Dict{Layer,Dict{Symbol,AbstractArray}},
        cLayer::L = nothing,
        BCache::Dict{Layer,Dict{Symbol,AbstractArray}}=Dict{Layer,Dict{Symbol,AbstractArray}}(),
        cnt = -1;
        tMiniBatch::Integer = -1, #can be used to perform both back and update params
        kwargs...,
    ) where {L<:Union{Layer,Nothing},T1,T2,N1,N2}

# Arguments

- `X` := train data

- `Y` := train labels

- `model` := is the model to perform the back propagation on

- `FCache` := the cached values of the forward propagation as `Dict{Layer, Dict{Symbol, AbstractArray}}`

- `cLayer` := is an internal variable to hold the current layer

- `BCache` := to hold the cache of the back propagtion (internal variable)

- `cnt` := is an internal variable to count the step of back propagation currently on to avoid re-do it

## Key-word Arguments

- `tMiniBatch` := to perform both the back prop and update trainable parameters in the same recursive call (if less than 1 update during back propagation is ditched)

- `kwargs` := other key-word arguments to be bassed to `layerBackProp` methods

# Return

- `BCache` := the cached values of the back propagation

"""
function chainBackProp(
    X::AbstractArray{T1,N1},
    Y::AbstractArray{T2,N2},
    model::Model,
    FCache::Dict{Layer,Dict{Symbol,AbstractArray}},
    cLayer::L = nothing,
    BCache::Dict{Layer,Dict{Symbol,AbstractArray}}=Dict{Layer,Dict{Symbol,AbstractArray}}(),
    cnt = -1;
    tMiniBatch::Integer = -1, #can be used to perform both back and update params
    kwargs...,
) where {L<:Union{Layer,Nothing},T1,T2,N1,N2}
    if cnt < 0
        cnt = model.outLayer.backCount + 1
    end

    if cLayer == nothing
        cLayer = model.outLayer
        BCache[cLayer] =
            layerBackProp(cLayer, model, FCache, BCache; labels = Y, kwargs...)

        if tMiniBatch > 0
            layerUpdateParams!(
                model,
                cLayer,
                cnt;
                tMiniBatch = tMiniBatch,
                kwargs...,
            )
        end

        BCache = chainBackProp(
            X,
            Y,
            model,
            FCache,
            cLayer.prevLayer,
            BCache,
            cnt;
            tMiniBatch = tMiniBatch,
            kwargs...,
        )

    elseif cLayer isa MILayer
        BCache[cLayer] = layerBackProp(cLayer, model, FCache, BCache; kwargs...)

        if tMiniBatch > 0
            layerUpdateParams!(
                model,
                cLayer,
                cnt;
                tMiniBatch = tMiniBatch,
                kwargs...,
            )
        end

        if cLayer.backCount >= cnt #in case layerBackProp did not do the back
            #prop becasue the next layers are not all
            #done yet
            for prevLayer in cLayer.prevLayer
                BCache = chainBackProp(
                    X,
                    Y,
                    model,
                    FCache,
                    prevLayer,
                    BCache,
                    cnt;
                    tMiniBatch = tMiniBatch,
                    kwargs...,
                )
            end #for
        end
    else #if cLayer==nothing
        BCache[cLayer] = layerBackProp(cLayer, model, FCache, BCache; kwargs...)

        if tMiniBatch > 0
            layerUpdateParams!(
                model,
                cLayer,
                cnt;
                tMiniBatch = tMiniBatch,
                kwargs...,
            )
        end

        if cLayer.backCount >= cnt #in case layerBackProp did not do the back
            #prop becasue the next layers are not all
            #done yet
            if !(cLayer isa Input)
                BCache = chainBackProp(
                    X,
                    Y,
                    model,
                    FCache,
                    cLayer.prevLayer,
                    BCache,
                    cnt;
                    tMiniBatch = tMiniBatch,
                    kwargs...,
                )
            end #if cLayer.prevLayer == nothing
        end
    end #if cLayer==nothing

    return BCache
end #backProp

export chainBackProp

###update parameters

include("layerUpdateParams.jl")

@doc raw"""
    chainUpdateParams!(model::Model,
                       cLayer::L=nothing,
                       cnt = -1;
                       tMiniBatch::Integer = 1) where {L<:Union{Layer,Nothing}}

Update trainable parameters using recursive call

# Arguments

- `model` := the model holds the training and update process

- `cLayer` := internal variable for recursive call holds the current layer

- `cnt` := an internal variable to hold the count of update in each layer not to re-do it

## Key-word Arguments

- `tMiniBatch` := the number of mini-batch of the total train collection

# Return

- `nothing`

"""
function chainUpdateParams!(model::Model,
                           cLayer::L=nothing,
                           cnt = -1;
                           tMiniBatch::Integer = 1) where {L<:Union{Layer,Nothing}}

    if cnt < 0
        cnt = model.outLayer.updateCount + 1
    end

    if cLayer==nothing
        layerUpdateParams!(model, model.outLayer, cnt, tMiniBatch=tMiniBatch)
        chainUpdateParams!(model, model.outLayer.prevLayer, cnt, tMiniBatch=tMiniBatch)

    elseif cLayer isa MILayer
        #update the AddLayer updateCounter
        if cLayer.updateCount < cnt
            cLayer.updateCount += 1
        end
        for prevLayer in cLayer.prevLayer
            chainUpdateParams!(model, prevLayer, cnt, tMiniBatch=tMiniBatch)
        end #for
    else #if cLayer==nothing
        layerUpdateParams!(model, cLayer, cnt, tMiniBatch=tMiniBatch)
        if !(cLayer isa Input)
            chainUpdateParams!(model, cLayer.prevLayer, cnt, tMiniBatch=tMiniBatch)
        end #if cLayer.prevLayer == nothing
    end #if cLayer==nothing

    return nothing
end #function chainUpdateParams!

export chainUpdateParams!

### train

@doc raw"""
    train(
          X_train,
          Y_train,
          model::Model,
          epochs;
          testData = nothing,
          testLabels = nothing,
          kwargs...,
          )

Repeat the trainging (forward/backward propagation and update parameters)

# Argument

- `X_train` := the training data

- `Y_train` := the training labels

- `model`   := the model to train

- `epochs`  := the number of repetitions of the training phase

# Key-word Arguments

- `testData` := to evaluate the training process over test data too

- `testLabels` := to evaluate the training process over test data too

- `batchSize` := the size of training when mini batch training

` `useProgBar` := (true, false) value to use prograss bar

- `kwargs` := other key-word Arguments to pass for the lower functions in hierarchy

# Return

- A `Dict{Symbol, Vector}` of:

    * `:trainAccuracies` := an `Array` of the accuracies of training data at each epoch
    * `:trainCosts` := an `Array` of the costs of training data at each epoch
    * In case `testDate` and `testLabels` are givens:

        + `:testAccuracies` := an `Array` of the accuracies of test data at each epoch
        + `:testCosts` := an `Array` of the costs of test data at each epoch

"""
function train(
               X_train,
               Y_train,
               model::Model,
               epochs;
               testData = nothing,
               testLabels = nothing,
               kwargs...,
               )

    kwargs = Dict{Symbol, Any}(kwargs...)
    batchSize = getindex(kwargs, :batchSize; default = 32)
    printCostsInterval = getindex(kwargs, :printCostsInterval; default = 0)
    useProgBar = getindex(kwargs, :useProgBar; default = true)
    embedUpdate = getindex(kwargs, :embedUpdate; default = true)
    metrics = getindex(kwargs, :metrics; default = [:accuracy, :cost])
    Accuracy = :accuracy in metrics
    Cost = :cost in metrics
    test = testData != nothing && testLabels != nothing

    inLayer, outLayer, lossFun, α = model.inLayer, model.outLayer, model.lossFun, model.α

    outAct = outLayer.actFun
    m = size(X_train)[end]
    c = size(Y_train)[end-1]
    nB = m ÷ batchSize

    N = ndims(X_train)
    axX = axes(X_train)[1:end-1]
    axY = axes(Y_train)[1:end-1]
    nMiniBatch = (nB + Integer(m % batchSize != 0))
    if useProgBar
        p = Progress(epochs*nMiniBatch, 0.1)
        showValues = Dict{String, Real}("Epoch $(epochs)"=>0,
                                        "Instances $(m)"=>0,
                                        )
    end

    if Accuracy
        showValues["Train Accuracy"] = 0.0
    end
    if Cost
        showValues["Train Cost"] = 0.0
    end

    if test
        testAcc = []
        testCost = []
        showValues["Test Accuracy"] = 0.0
        showValues["Test Cost"] = 0.0
    end

    Accuracies = []
    Costs = []
    for i=1:epochs
        shufInd = randperm(m)
        minCosts = [] #the costs of all mini-batches
        minAcc = []
        for j=1:nB
            downInd = (j-1)*batchSize+1
            upInd   = j * batchSize
            batchInd = shufInd[downInd:upInd]
            X = X_train[axX..., batchInd]
            Y = Y_train[axY..., batchInd]

            FCache = chainForProp(X, inLayer; kwargs...)
            a = FCache[outLayer][:A]
            if Accuracy
                _, acc = probToValue(eval(:($outAct)), a; labels = Y)
                push!(minAcc, acc)
                if useProgBar
                    tmpAcc = mean(minAcc)
                    showValues["Train Accuracy"] = round(tmpAcc ; digits=4)
                end
            end

            if Cost
                minCost = cost(eval(:($lossFun)), a, Y)
                push!(minCosts, minCost)
                if useProgBar
                    tmpCost = mean(minCosts)
                    showValues["Train Cost"] = round(tmpCost; digits = 4)
                end
            end

            if embedUpdate
                BCache = chainBackProp(X,Y,
                                       model,
                                       FCache;
                                       tMiniBatch = j,
                                       kwargs...)
            else
                BCache = chainBackProp(X,Y,
                                      model,
                                      FCache;
                                      kwargs...)

                chainUpdateParams!(model; tMiniBatch = j)
            end #if embedUpdate
            if useProgBar
                showValues["Epoch $(epochs)"] = i
                showValues["Instances $(m)"] = j*batchSize
                update!(p, (i-1)*nMiniBatch+j; showvalues=[("Epoch $(epochs)", pop!(showValues, "Epoch $(epochs)")), ("Instances $(m)", pop!(showValues, "Instances $(m)")), showValues...])
            end
                # if :accuracy in metrics && :cost in metrics
                #     update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+j; showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize), (:Accuracy, round(mean(minAcc); digits=4)), (:Cost, round(mean(minCosts); digits=4))])
                # elseif :accuracy in metrics
                #     update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+j; showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize), (:Accuracy, round(mean(minAcc); digits=4))])
                # elseif :cost in metrics
                #     update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+j; showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize), (:Cost, round(mean(minCosts); digits=4))])
                # else
                #     update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+j; showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize)])
                # end
            # end
            # chainUpdateParams!(model; tMiniBatch = j)

        end #for j=1:nB iterate over the mini batches

        if m%batchSize != 0
            downInd = (nB)*batchSize+1
            batchInd = shufInd[downInd:end]
            X = X_train[axX..., batchInd]
            Y = Y_train[axY..., batchInd]

            FCache = chainForProp(X, inLayer; kwargs...)
            if Accuracy
                _, acc = probToValue(eval(:($outAct)), FCache[outLayer][:A]; labels = Y)
                push!(minAcc, acc)
                if useProgBar
                    tmpAcc = mean(minAcc)
                    showValues["Train Cost"] = round(tmpAcc ; digits=4)
                end
            end

            a = FCache[outLayer][:A]
            if Cost
                minCost = cost(eval(:($lossFun)), a, Y)
                push!(minCosts, minCost)
                if useProgBar
                    tmpCost = mean(minCosts)
                    showValues["Train Cost"] = round(tmpCost; digits = 4)
                end
            end

            if embedUpdate
                BCache = chainBackProp(X,Y,
                                       model,
                                       FCache;
                                       tMiniBatch = nB + 1,
                                       kwargs...)
            else
                BCache = chainBackProp(X,Y,
                                      model,
                                      FCache;
                                      kwargs...)

                chainUpdateParams!(model; tMiniBatch = nB + 1)
            end #if embedUpdate
        end

        if Accuracy
            tmpAcc = mean(minAcc)
            push!(Accuracies, tmpAcc)
            if useProgBar
                showValues["Train Accuracy"] = round(tmpAcc; digits=4)
            end
        end

        if Cost
            tmpCost = mean(minCosts)
            push!(Costs, tmpCost)
            if useProgBar
                showValues["Train Cost"] = round(tmpCost; digits=4)
            end
        end

        if printCostsInterval>0 && i%printCostsInterval==0
            println("N = $i, Cost = $(Costs[end])")
        end

        if test
            testDict = predict(model, testData, testLabels; useProgBar = false, batchSize=batchSize)
            testcost = cost(eval(:($lossFun)), testDict[:YhatProb], testLabels)
            if Accuracy
                push!(testAcc, testDict[:accuracy])
                if useProgBar
                    showValues["Test Accuracy"] = round(testDict[:accuracy]; digits=4)
                end
            end
            if Cost
                push!(testCost, testcost)
                if useProgBar
                    showValues["Test Cost"] = round(testcost; digits=4)
                end
            end
        end

        if useProgBar
            showValues["Epoch $(epochs)"] = i
            showValues["Instances $(m)"] = m
            update!(p, (i-1)*nMiniBatch+(nB+1); showvalues=[("Epoch $(epochs)", pop!(showValues, "Epoch $(epochs)")), ("Instances $(m)", pop!(showValues, "Instances $(m)")), showValues...])
        end

        # if useProgBar
        #     if useProgBar
        #         if :accuracy in metrics && :cost in metrics
        #             update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+(nB + 1); showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", m), (:Accuracy, round(mean(minAcc); digits=4)), (:Cost, round(mean(minCosts); digits=4))])
        #         elseif :accuracy in metrics
        #             update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+(nB + 1); showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize), (:Accuracy, round(mean(minAcc); digits=4))])
        #         elseif :cost in metrics
        #             update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+(nB + 1); showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize), (:Cost, round(mean(minCosts); digits=4))])
        #         else
        #             update!(p, ((i-1)*(nB + ((m % batchSize == 0) ? 0 : 1)))+(nB + 1); showvalues=[("Epoch ($epochs)", i), ("Instances ($m)", j*batchSize)])
        #         end
        #     end
        # end
    end

    # model.W, model.B = W, B
    outDict = Dict{Symbol, Array{T,1} where {T}}()
    if Accuracy
        outDict[:trainAccuracies] = Accuracies
        if test
            outDict[:testAccuracies] = testAcc
        end
    end

    if Cost
        outDict[:trainCosts] = Costs
        if test
            outDict[:testCosts] = testCost
        end
    end
    return outDict
end #train

export train


text:
function Base.show(io::IO, csinfo::CSInfo)
    println("Info: ")
    for name in fieldnames(CSInfo)
        println(io, "$name = $(getfield(csinfo, name))")
    end
end

function compress_var_string(variable::CS.Variable)
    if CS.isfixed(variable)
        return string(CS.value(variable))
    end

    sorted_vals = sort(CS.values(variable))
    if sorted_vals[1] + length(CS.values(variable)) - 1 == sorted_vals[end]
        return "$(sorted_vals[1]):$(sorted_vals[end])"
    end
    return string(sort(CS.values(variable)))
end

function get_str_repr(variables::Array{CS.Variable})
    if length(size(variables)) == 1
        output = ""
        for i in 1:length(variables)
            if !CS.isfixed(variables[i])
                output *= "$(compress_var_string(variables[i])), "
            else
                output *= "$(CS.value(variables[i])), "
            end
        end
        return [output[1:(end - 2)]]
    elseif length(size(variables)) == 2
        max_length = 1
        for j in 1:size(variables)[2]
            for i in 1:size(variables)[1]
                if !CS.isfixed(variables[i, j])
                    len = length(compress_var_string(variables[i, j]))
                    if len > max_length
                        max_length = len
                    end
                else
                    len = length(string(CS.value(variables[i, j])))
                    if len > max_length
                        max_length = len
                    end
                end
            end
        end
        lines = String[]
        for j in 1:size(variables)[2]
            line = ""
            for i in 1:size(variables)[1]
                pstr = ""
                if !CS.isfixed(variables[i, j])
                    pstr = compress_var_string(variables[i, j])
                else
                    pstr = string(CS.value(variables[i, j]))
                end
                space_left = floor(Int, (max_length - length(pstr)) / 2)
                space_right = ceil(Int, (max_length - length(pstr)) / 2)
                line *= repeat(" ", space_left) * pstr * repeat(" ", space_right) * " "
            end
            push!(lines, line)
        end
        return lines
    else
        @warn "Currently not supported to print more than 2 dimensions. Maybe file an issue with your problem and desired output"
    end
end

function Base.show(io::IO, variables::Array{CS.Variable})
    lines = get_str_repr(variables)
    for line in lines
        println(line)
    end
end


text:
using StatsBase: mean

include("decisiontree.jl")

"""Gradient tree boosting algorithm.

The implmentation follows **Algorithm 10.3** in the *Elements of
Statistical Learning*. The original implementation of this algorithm called
MART for "multiple additive regression tree.
"""

# Loss functions
abstract type LossFunction end

struct MeanSquaredError <: LossFunction end

function (mse::MeanSquaredError)(y, pred)
    mean((y - pred).^2)
end

negativegradient(mse::MeanSquaredError, y, pred) = y - pred

"Initial prediction in case of Meansquared error is global mean."
predict(lossfunc::MeanSquaredError, y) = mean(y)

struct LogLoss <: LossFunction end

"""
LogLoss for binary classification problem.
The implmentation follows the article:
http://zpz.github.io/blog/gradient-boosting-tree-for-binary-classification/.
Here `pred` is the predicted log-odds and logistic/sigmmoid function
converted the log-odds to the predicted probability.
"""
function (logloss::LogLoss)(y, pred)
    sum(log.(1+exp.(pred)) - y.*pred)
end

sigmoid(pred) = 1 / (1 + exp(-pred))

gradient(logloss::LogLoss, y, pred) =  sigmoid.(pred) .- y

function hessian(logloss::LogLoss, y, pred)
    p = sigmoid.(pred)
    return p.*(1 .- p)
end
negativegradient(logloss::LogLoss, y, pred) =  -gradient(logloss, y, pred)

function predict(logloss::LogLoss, y)
    numobs = length(y)
    numpos = sum(y)
    return log(numpos/(numobs-numpos))
end

"""
Stopping rule for boosting model.
TODO:
    - Check the minimal increase between two consecutive model updates.
"""
struct BoostingStoppingRule
    numtrees::Int
    treestoppingrule::TreeStoppingRule
end

# Gradient boosting model
mutable struct GradientBoostingTrees
    lossfunc::LossFunction
    boostingstoppingrule::BoostingStoppingRule
    trees::Vector{DecisionTree}
    learningrate::Float64
    initialpred::Float64

    function GradientBoostingTrees(X, y, lossfunc, boostingstoppingrule; learningrate=0.1)
        numtrees = boostingstoppingrule.numtrees
        treestoppingrule = boostingstoppingrule.treestoppingrule
        trees = RegressionTree[]

        # Initialize boosting model: f0(x)
        initialpred = predict(lossfunc, y)
        preds = fill(initialpred, length(y))

        for m = 1:numtrees
            # generalized residual (negative gradient)
            residuals = negativegradient(lossfunc, y, preds)
            tree = RegressionTree(X, residuals, treestoppingrule)

            # update predictions on the terminal regions (leaves) in the newly fitted tree
            updateleaves!(lossfunc, tree, X, y, preds)
            newtreepreds = predict(tree, X)
            preds += learningrate * newtreepreds
            push!(trees, tree)
        end

        new(lossfunc, boostingstoppingrule, trees, learningrate, initialpred)
    end

end

"""
    updateleaves!(lossfunc, tree, X)

Update tree leaves with Meansquared error as the loss function.
For Meansquared error we have a closed form solution to minimize the loss
w.r.t the current model: the optimal point the mean of the residual in each
leave. See https://explained.ai/gradient-boosting/descent.html for details.

Because the regression tree has already fitted the residual, we don't need to
change the leaf value in this case.
"""
function updateleaves!(lossfunc, tree, X, y, preds) end

"""
    updateleaves!(lossfunc::LogLoss, tree, X)

Update tree leaves for log-loss. In general, we don't have a closed form solution.
We will conduct one step update using Newton's method. To use Newton's method,
we need gradient, Hessian and regularization terms.
"""
function updateleaves!(lossfunc::LogLoss, tree, X, y, preds)
    numregions = tree.numregions
    residuals = negativegradient(lossfunc, y, preds)
    hessians = hessian(lossfunc, y, preds)

    regionresidualsums = zeros(numregions)
    regionhessiansums = zeros(numregions)

    for i in 1:size(X)[1]
        regiondid = findregionid(lossfunc, tree, X[i, :])
        regionresidualsums[regiondid] += residuals[i]
        regionhessiansums[regiondid] += hessians[i]
    end

    regionupdates = regionresidualsums ./ regionhessiansums

    replacepreds!(tree, regionupdates)
end

function findregionid(lossfunc::LogLoss, tree, x)
    node = tree.rootnode
    while !node.isleaf
        featureindex = node.featureindex
        thresh = node.thresh
        if x[featureindex] < thresh
            node = node.leftnode
        else
            node = node.rightnode
        end
    end

    return node.regionid
end

function replacepreds!(tree::DecisionTree, regionupdates)
    stack = [tree.rootnode]

    while !isempty(stack)
        node = pop!(stack)
        if node.isleaf
            node.prediction = regionupdates[node.regionid]
        else
            @assert !isnothing(node.leftnode)
            @assert !isnothing(node.rightnode)
            push!(stack, node.leftnode)
            push!(stack, node.rightnode)
        end
    end

end

function predict(boostingmodel::GradientBoostingTrees, X)
    lossfunc = boostingmodel.lossfunc
    learningrate = boostingmodel.learningrate
    preds = fill(boostingmodel.initialpred, size(X)[1])

    for tree in boostingmodel.trees
        preds += learningrate*predict(tree, X)
    end

    preds
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
using Statistics
import Base.Meta: isexpr

"""
    ret = @freshexec [setup_ex] ex

Runs `ex` in an external process and gets back the final result (, which is supposed to be
  such a simple Julia object that we can restore it from its string representation).
Running  in external process can be useful for testing JET analysis, because:
- the first time analysis is not affected by native code cache and JET's global report cache
- we can do something that might break a Julia process without wondering it breaks the
  original test process later on
The optional positional argument `setup_ex` runs before each execution of `ex` and defaults
  to `JET_LOAD_EX`, which just loads JET into the process.
"""
macro freshexec(args...)
    args = map(a->Expr(:quote,a), args)
    return Expr(:escape, Expr(:call, GlobalRef(@__MODULE__, :freshexec), __module__, args...))
end
freshexec(mod, ex) = freshexec(mod, JET_LOAD_EX, ex)
freshexec(args...) = _freshexec(args..., collect_last_result)

function collect_last_result(exs)
    lines = string.(exs)
    lines[end] = "ret = $(lines[end])"
    return join(lines, '\n')
end

"""
    stats = @benchmark_freshexec [setup_ex] ex

Runs `ex` in an external process and collects execution statistics from [`@timed`](@ref).
This is particularly useful for testing the performance of first-time analysis, where
  the native code cache and JET's global report cache have no effect for its performance.
The optional positional argument `setup_ex` runs before each execution of `ex` and its
  execution statistis are not included in the benchmark result; it defaults to
  `JET_WARMUP_EX`, which loads JET and runs a warm up analysis `@profile identity(nothing)`.
"""
macro benchmark_freshexec(args...)
    args = map(a->Expr(:quote,a), args)
    return Expr(:escape, Expr(:call, GlobalRef(@__MODULE__, :benchmark_freshexec), __module__, args...))
end
benchmark_freshexec(mod, ex) = benchmark_freshexec(mod, JET_WARMUP_EX, ex)
benchmark_freshexec(args...) = _freshexec(args..., collect_statistics)

function collect_statistics(exs)
    return """
    ret = @timed begin
        $(join(exs, '\n'))
        nothing # ensure `stats` can be parsed
    end
    """
end

"""
    stats = @nbenchmark_freshexec [ntimes = 5] [setup_ex] ex

Runs `ex` in an external process multiple times (which can be configured by the optional
  keyword argument `ntimes`), and collects execution statistics from [`@timed`](@ref).
  The statistics are generated by taking the `mean` of all the trials.
This is particularly useful for benchmarking the performance of first-time analysis, where
  the native code cache and JET's global report cache have no effect for its performance.
The optional positional argument `setup_ex` runs before each execution of `ex` and its
  execution statistis are not included in the benchmark result; it defaults to
  `JET_WARMUP_EX`, which loads JET and runs a warm up analysis `@profile identity(nothing)`.
"""
macro nbenchmark_freshexec(args...)
    isn(x) = isexpr(x, :(=)) && first(x.args) === :ntimes
    i = findfirst(isn, args)
    ntimes = i === nothing ? 5 : Expr(:quote, last(args[i].args))
    args = map(a->Expr(:quote,a), filter(!isn, args))
    return Expr(:escape, Expr(:call, GlobalRef(@__MODULE__, :nbenchmark_freshexec), ntimes, __module__, args...))
end
nbenchmark_freshexec(ntimes, mod, ex) = nbenchmark_freshexec(ntimes, mod, JET_WARMUP_EX, ex)
function nbenchmark_freshexec(ntimes, args...)
    stats = [_freshexec(args..., collect_statistics) for _ in 1:ntimes]
    return (; time   = mean(getproperty.(stats, :time)),
              bytes  = mean(getproperty.(stats, :bytes)),
              gctime = mean(getproperty.(stats, :gctime)),
              )
end

function _freshexec(mod, setup_ex, ex, exs2script)
    # we need to flatten block expression into a toplevel expression to correctly handle
    # e.g. macro expansions
    setup_exs = isexpr(setup_ex, :block) ? setup_ex.args : [setup_ex]
    setup_script = join(string.(setup_exs), '\n')
    exs = isexpr(ex, :block) ? ex.args : [ex]
    script = exs2script(exs)

    prog = """
    old = stdout
    rw, wr = redirect_stdout()
    $(setup_script)
    $(script)
    redirect_stdout(old)
    close(rw); close(wr)
    println(stdout, repr(ret))
    """

    cmd = Cmd([JULIA_BIN, "-e", prog])
    io = IOBuffer()
    run(pipeline(cmd; stdout = io))

    return Core.eval(mod, Meta.parse(String(take!(io))))
end

const JET_LOAD_EX = :(using JET)
const JET_WARMUP_EX = quote
    using JET
    @profile_call identity(nothing) # warm-up for JET analysis
end
const JULIA_BIN = normpath(Sys.BINDIR, "julia")


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_VERBOSE)

Sparso.dprintln(0, 1, "should not print \n as the \n level is too low")
Sparso.dprintln(1, 0, "Testing tabbed display:")
Sparso.dprintln(1, 1, "1 tab")
Sparso.dprintln(1, 1, "OK")
Sparso.dprintln(1, 2, "2 tabs", "\nOK")
Sparso.dprintln(1, 3, "3 tabs", "\nOK")
flush(STDOUT)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
"""
check https://github.com/thautwarm/MLStyle.jl/blob/master/matrix_benchmark.jl
"""

include("matrix-benchmark/sampler.jl")
include("matrix-benchmark/utils.jl")
export ArbitrarySampler
export Utils

versus_items = ("datatype", "misc", "tuple", "array", "structfields", "vs-match")

function run_all()
    for item in versus_items
        run_one(item)
    end
end

function run_one(item)
    filename = "matrix-benchmark/bench-$item.jl"
    open(filename) do f
        include_string(Main, read(f, String), filename)
    end
end


text:
using PowerSystems
using NLsolve

include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
omib_file_dir = joinpath(dirname(@__FILE__), "OMIB.raw")
omib_sys = System(omib_file_dir, runchecks = false)
add_source_to_ref(omib_sys)
############### Data Dynamic devices ########################
function dyn_gen_classic(generator)
    return DynamicGenerator(
        name = get_name(generator),
        ω_ref = 1.0,
        machine = machine_classic(),
        shaft = shaft_damping(),
        avr = avr_none(),
        prime_mover = tg_none(),
        pss = pss_none(),
    )
end

#Attach dynamic generator. Currently use PSS/e format based on bus #.
gen = [g for g in get_components(Generator, omib_sys)][1]
case_gen = dyn_gen_classic(gen)
add_component!(omib_sys, case_gen, gen)

#Compute Y_bus after fault
fault_branch = deepcopy(collect(get_components(Branch, omib_sys))[1])
fault_branch.r = 0.00;
fault_branch.x = 0.1;
Ybus_fault = PSY.Ybus([fault_branch], get_components(Bus, omib_sys))[:, :]


text:
function update!(args...)
    dict = args[end]
    btns = ["btn1";"btn2";"btn3";"btn4"]
    txts = ["txt1";"txt2";"txt3";"txt4"]

    dict[btns[1]][] == 0 && dict[btns[2]][] == 0 && dict[btns[3]][] == 0 && dict[btns[4]][] == 0 && return

    for (i,btn) in enumerate(btns)
        if dict[btn][] == 1
            dict[btn]["is-loading"][] = true
            println("Starting simulation...")

            i==1 && run1!(dict[txts[i]][])
            i==2 && run2!(dict[txts[i]][])
            i==3 && run3!(dict[txts[i]][])
            i==4 && run4!(dict[txts[i]][])

            dict[btn]["is-loading"][] = false
            dict[btn][] = 0
        end
    end

    return
end

function friction(Fcart, Fpend, v, ω, k, c, Fs)
    if abs(v) < 1e-3 && abs(Fcart) < Fs
        Fcart = 0
    end
    Fcart = Fcart - k*v
    Fpend = -c*ω

    return Fcart, Fpend
end

function friction(Fcart, Fpend1, Fpend2, v, ω1, ω2, k, c, Fs)
    if abs(v) < 1e-3 && abs(Fcart) < Fs
        Fcart = 0
    end
    Fcart = Fcart - k*v
    Fpend1 = -c*ω1
    Fpend2 = -c*ω2

    return Fcart, Fpend1, Fpend2
end


text:
if VERSION >= v"0.4-"
  __precompile__()
end

module Qwt

# # VERSION >= v"0.4.0-dev+6521" && __precompile__()
# __precompile__(false)

export
  plot,
  oplot,
  subplot,
  scatter,
  heatmap,

  currentPlot,
  currentPlot!,

  getline,
  setdata,
  addRegressionLine,
  refresh,
  title,
  xlabel,
  ylabel,
  yrightlabel,
  windowtitle,
  hidelegend,
  showlegend,

  Widget,
  PlotWidget,
  Plot,

  hidewidget,
  showwidget,
  widgetpos,
  widgetsize,
  movewidget,
  resizewidget,
  move_resizewidget,
  screenCount,
  screenPosition,
  screenSize,
  moveToScreen,
  moveToLastScreen,

  Layout,
  Splitter,
  vbox,
  hbox,
  vsplitter,
  hsplitter,

  savepng,
  animation,
  saveframe,
  makegif,

  P2, # 2D point
  P3, # 3D point
  Point,
  ORIGIN,
  Scene,
  SceneItem,
  Ellipse,
  Circle,
  Rect,
  Square,
  SceneText,
  Line,

  currentScene,
  currentScene!,
  defaultBrush,
  defaultBrush!,
  defaultPen,
  defaultPen!,
  top,
  bottom,
  left,
  right,
  topleft,
  bottomright,
  background!,
  foreground!,

  makecolor,
  makebrush,
  makepen,
  position3d,
  position!,
  zvalue,
  zvalue!,
  rotation,
  rotation!,
  visible,
  visible!,
  parent,
  parent!,
  brush!,
  pen!,
  move!,
  movex!,
  movey!,
  scale!,
  scalex!,
  scaley!,
  rotate!,
  settext,

  ellipse!,
  circle!,
  rect!,
  square!,
  text!,
  line!

# ------------------------------------------------------
# Setup Qt
# ------------------------------------------------------

function addToPythonPath()
  # add the python dir to PYTHONPATH
  qwtPythonDir = "$(Pkg.dir("Qwt"))/src/python"
  try
      ENV["PYTHONPATH"] = ENV["PYTHONPATH"] * ":" * qwtPythonDir
  catch
      ENV["PYTHONPATH"] = qwtPythonDir
  end
end

addToPythonPath()

# # add the python dir to PYTHONPATH
# qwtPythonDir = "$(Pkg.dir("Qwt"))/src/python"
# try
#     ENV["PYTHONPATH"] = ENV["PYTHONPATH"] * ":" * qwtPythonDir
# catch
#     ENV["PYTHONPATH"] = qwtPythonDir
# end

using Colors
using PyCall

  # print("Initializing Qwt... ")
  # global QT, PLOT, WIDGETS, QAPP

  # # set the PYTHONPATH
  # unshift!(PyVector(pyimport("sys")["path"]), "")

  # # global QT, QWT, PLOT, WIDGETS
  # @pyimport PyQt4.Qt as QT
  # @pyimport PyQt4.Qwt5 as QWT
  # @pyimport BasicPlot as PLOT
  # @pyimport pythonwidgets as WIDGETS

# const QT = PyCall.PyNULL()
# const QWT = PyCall.PyNULL()
# const PLOT = PyCall.PyNULL()
# const WIDGETS = PyCall.PyNULL()

function __init__()

  # # add the python dir to PYTHONPATH
  # qwtPythonDir = "$(Pkg.dir("Qwt"))/src/python"
  # try
  #     ENV["PYTHONPATH"] = ENV["PYTHONPATH"] * ":" * qwtPythonDir
  # catch
  #     ENV["PYTHONPATH"] = qwtPythonDir
  # end

  addToPythonPath()

  # copy!(QT, pyimport("PyQt4"))
  global const QT = pywrap(pyimport("PyQt4.Qt"))
  global const QWT = pywrap(pyimport("PyQt4.Qwt5"))
  global const PLOT = pywrap(pyimport("BasicPlot"))
  global const WIDGETS = pywrap(pyimport("pythonwidgets"))

  # QT = QT2
  # PLOT = PLOT2
  # WIDGETS = WIDGETS2
  # unshift!(PyVector(pyimport("sys")["path"]), "")  # so you can load python files from the current directory
  # @pyimport FancyPlot as FPLOT
  pygui_start(:qt_pyqt4)
  global const QAPP = QT.QApplication([])
  # println("done.")

  global const CURRENT_SCENE = CurrentScene(Nullable{Scene}(), makebrush(:black), makepen(2, :black))
end
# ------------------------------------------------------

using ImmutableArrays

typealias P2 Vector2{Float64}
typealias P3 Vector3{Float64}
typealias Point Union{P2,P3}

const ORIGIN = P3(0,0,0)

P2(p::P3) = P2(p[1], p[2])
P3(p::P2, z::Real) = P3(p[1], p[2], z)
P3(p::P2) = P3(p, 0.0)

include("widgets.jl")
include("layout.jl")
include("plotitems.jl")
include("plot.jl")
include("subplot.jl")
include("animation.jl")
include("scene.jl")

end