text:
#!/usr/local/bin/julia
# last tested Julia version: 1.6.1
# Kuramoto model on a 2D lattice
# FvW 03/2018

using Interpolations
using NPZ
using PyCall
using PyPlot
using Statistics
using VideoIO
@pyimport matplotlib.animation as anim

function kuramoto2d(N, T, t0, dt, s, K0, K1, sd_omega)
    s_sqrt_dt = s*sqrt(dt)
	omega = 1.0 .+ sd_omega * (rand(N,N) .- 0.5)
    phi = 2pi .* rand(N,N) .- pi
	X = zeros(Float64,T,N,N)
	dphi = zeros(Float64,N,N)
	phi_up = zeros(Float64,N,N)
	phi_down = zeros(Float64,N,N)
	phi_left = zeros(Float64,N,N)
    phi_right = zeros(Float64,N,N)
    # coupling constant time course
    u = K0*ones(t0)
    v = [_ for _ in range(K0,stop=K1,length=T)]
    K = vcat(u,v) # 1-dim time course of K
    K4 = K .* 0.25
    # iterate
    for t in range(1, stop=t0+T, step=1)
        (t%100 == 0) && print("    t = ", t, "/", t0+T, "\r")
		phi_up = circshift(phi, [-1 0])
        phi_down = circshift(phi, [1 0])
        phi_left = circshift(phi, [0 -1])
        phi_right = circshift(phi, [0 1])
        dphi = zeros(Float64,N,N)
	    dphi += (sin.(phi_up-phi) + sin.(phi_down-phi))
	    dphi += (sin.(phi_left-phi) + sin.(phi_right-phi))
        dphi = omega .+ K4[t] .* dphi
		# stochastic integration
		phi += (dphi*dt + s_sqrt_dt*randn(N,N))
		#(t > t0) && (X[t-t0,:,:] = cos.(phi))
        (t > t0) && (X[t-t0,:,:] = angle.(exp.(im.*phi))) # angle
    end
    println("\n")
    return X
end

function animate_pyplot(fname, data)
    """
    Animate 3D array as .mp4 using PyPlot, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    vmin = minimum(data)
    vmax = maximum(data)
    # setup animation image
    fig = figure(figsize=(6,6))
    axis("off")
    t = imshow(data[1,:,:], origin="lower", cmap=ColorMap("gray"),
			   vmin=vmin, vmax=vmax)
    tight_layout()
    # frame generator
    println("[+] animate")
    function animate(i)
        (i%100 == 0) && print("    t = ", i, "/", nt, "\r")
        t.set_data(data[i+1,:,:])
    end
    # create animation
    ani = anim.FuncAnimation(fig, animate, frames=nt, interval=10)
    println("\n")
    # save animation
    ani[:save](fname, bitrate=-1,
               extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    show()
end

function animate_video(fname, data, downsample=10)
    """
    Animate 3D array as .mp4 using VideoIO, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    n1 = nt
    if (downsample > 0)
        println("[+] Downsampling ", downsample)
        n1 = Int(nt/downsample) # number of samples after downsampling
        data1 = zeros(Complex64,n1,nx,ny)
        t0 = range(1,stop=nt)
        t1 = range(1,stop=nt,length=n1)
        for i in range(1,stop=nx)
            for j in range(1,stop=ny)
                y = data[:,i,j]
                f_ip = LinearInterpolation(t0, y)
                data1[:,i,j] = f_ip(t1)
            end
        end
    else
        data1 = copy(data)
    end
    # BW
    y = UInt8.(round.(255*(data1 .- minimum(data1)) /
                           (maximum(data1)-minimum(data1))))
    # BW inverted
    #y = UInt8.(round.(255 .- 255*(data1 .- minimum(data1)) /
    #                              (maximum(data1)-minimum(data1))))
    encoder_options = (color_range=2, crf=0, preset="medium")
    framerate=30
    T = size(data1,1)
    println("[+] Animate:")
    open_video_out(fname, y[1,end:-1:1,:], framerate=framerate,
                   encoder_options=encoder_options) do writer
        for i in range(2,stop=n1,step=1)
            (i%100 == 0) && print("    i = ", i, "/", T, "\r")
            write(writer, y[i,end:-1:1,:])
        end
    end
    println("\n")
end

function main()
    println("Kuramoto lattice model\n")
    N = 128
    T = 3000
    t0 = 0
    dt = 0.05
    s = 0.2
    # Kuramoto parameters
    K0 = 0.5
    K1 = 4.0
    sd_omega = 0.25
    println("[+] Lattice size N: ", N)
    println("[+] Time steps T: ", T)
    println("[+] Warm-up steps t0: ", t0)
    println("[+] Integration time step dt: ", dt)
    println("[+] Noise intensity s: ", s)
	println("[+] Coupling constant K0...K1: ", K0, " ... ", K1)
    println("[+] Frequency distribution std. dev.: ", sd_omega)

    # run simulation
    data = kuramoto2d(N, T, t0, dt, s, K0, K1, sd_omega)
    println("[+] Data dimensions: ", size(data))

    # plot mean signal
    m = mean(reshape(data, (T,N*N)), dims=2)
    plot(m, "-k"); show()

    # save data
    K0_str = rpad(K0, 4, '0') # coupl. param. 4-char string
    K1_str = rpad(K1, 4, '0') # coupl. param. 4-char string
    s_str = rpad(s, 4, '0') # noise as 4-char string
    fname1 = string("kuramoto2d_K_", K0_str, "_", K1_str, "_s_", s_str, ".npy")
    #npzwrite(fname1, data)
    #println("[+] Data saved as: ", fname1)

	# video
    fname2 = string("kuramoto2d_K_", K0_str, "_", K1_str, "_s_", s_str, ".mp4")
    #animate_pyplot(fname2, data) # slow
    animate_video(fname2, data, 5) # fast
    println("[+] Data saved as: ", fname2)
end

main()


text:
export make_offline_gui

mutable struct Analog_Viewer
    c::Gtk.GtkCanvasLeaf
    slider::Gtk.GtkScaleLeaf
    adj::Gtk.GtkAdjustmentLeaf
    pos::Int64
    changed::Bool
end

mutable struct Offline_GUI
    win::Gtk.GtkWindowLeaf
    sc::Single_Channel
    sortview::SortView
    v::Array{Int16,1} #Should be anything
    av::Analog_Viewer
end

function make_offline_gui(s)

    grid = Grid()

    vbox1_2=Grid()

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

    sb2=SpinButton(1:1000)
    setproperty!(sb2,:value,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")
    setproperty!(button_gain,:active,false)
    vbox1_2_1[1,2]=button_gain

    #THRESHOLD
    frame1_3=Frame("Threshold")
    vbox1_2[1,3]=frame1_3
    vbox1_3_1=Box(:v)
    push!(frame1_3,vbox1_3_1)

    #sb=SpinButton(-300:300)
    #setproperty!(sb,:value,0)
    sb=Label("0")
    push!(vbox1_3_1,sb)

    button_thres_all = CheckButton()
    add_button_label(button_thres_all,"All Channels")
    setproperty!(button_thres_all,:active,false)
    push!(vbox1_3_1,button_thres_all)

    button_thres = CheckButton()
    add_button_label(button_thres,"Show")
    setproperty!(button_thres,:active,false)
    push!(vbox1_3_1,button_thres)

    frame_hold=Frame("Spike")
    #grid[1,1]=frame_hold
    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

    button_rb=Array(RadioButton,3)
    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)
    setproperty!(adj_sort,:value,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 |> showall

    #COLUMN 2 - Threshold slider
    vbox_slider=Box(:v)
    thres_slider = Scale(true, -300,300,1)
    adj_thres = Adjustment(thres_slider)
    setproperty!(adj_thres,:value,0)

    c_thres=Canvas(10,200)
    setproperty!(c_thres,:vexpand,false)

    Gtk.GAccessor.inverted(thres_slider,true)
    Gtk.GAccessor.draw_value(thres_slider,false)

    setproperty!(thres_slider,:vexpand,true)
    push!(vbox_slider,thres_slider)
    push!(vbox_slider,c_thres)

    c_grid=Grid()

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

    show(c2)
    c_grid[1,1]=c2
    setproperty!(c2,:hexpand,true)
    setproperty!(c2,:vexpand,true)

    c3=Canvas(-1,200)
    @guarded draw(c3) do widget
        ctx = Gtk.getgc(c3)
        clear_c3(c3,1)
    end
    show(c3)
    c_grid[1,2]=c3
    setproperty!(c3,:hexpand,true)

    grid[2,2]=vbox1_2
    grid[3,2]=vbox_slider
    grid[4,2]=c_grid

    c_analog=Canvas(-1,200)
    grid[4,3]=c_analog
    analog_slider = Scale(false, 1,30000,1)
    adj_analog = Adjustment(analog_slider)
    setproperty!(adj_analog,:value,1)
    grid[4,4]=analog_slider

    analog_handles = Analog_Viewer(c_analog,analog_slider,adj_analog,1,false)

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

    win = Window(grid,"SpikeSorting.jl") |> showall

    sleep(5.0)

    thres_widgets=Thres_Widgets(sb,thres_slider,adj_thres,button_thres_all,button_thres)
    gain_widgets=Gain_Widgets(sb2,gain_checkbox,button_gain)

    sc_widgets=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,adj_thres,thres_slider,false,
    thres_widgets,gain_widgets)

    id = signal_connect(canvas_press_win,c2,"button-press-event",Void,(Ptr{Gtk.GdkEventButton},),false,(sc_widgets,))

    #=
    Pause, Restore, and Clear Callbacks
    =#
    id = signal_connect(pause_cb,button_pause,"toggled",Void,(),false,(sc_widgets,))
    id = signal_connect(clear_button_cb,button_clear,"clicked",Void,(),false,(sc_widgets,))
    id = signal_connect(restore_button_cb,button_restore,"clicked",Void,(),false,(sc_widgets,))

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

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

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

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

    id = signal_connect(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(check_cb_template,check_sort1,"clicked",Void,(),false,(sc_widgets,))

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

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

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

    #=
    Threshold callbacks
    =#

    id = signal_connect(thres_cb,thres_slider,"value-changed",Void,(),false,(sc_widgets,))
    id = signal_connect(thres_show_cb,button_thres,"clicked",Void,(),false,(sc_widgets,))

    #=
    Gain Callbacks
    =#

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

    v=round.(Int16,100 .* rand(60000)-50.0)

    myoffline = Offline_GUI(win,sc_widgets,sortview_handles,v,analog_handles)

    #=
    Analog Callbacks
    =#

    id = signal_connect(analog_slider_cb,analog_slider,"value-changed",Void,(),false,(myoffline,))
    id = signal_connect(analog_gain_cb, sb2, "value-changed",Void,(),false,(myoffline,))

    myoffline
end

function offline_loop(han,sorting,buf,nums)

    while true

        if han.sc.buf.c_changed

            clus = han.sc.buf.selected_clus

            if clus>0

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

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

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

            end

            #Send cluster information to sorting
            send_clus(sorting,han.sc)
            han.sc.buf.c_changed=false
        end
        if han.sc.buf.replot
            replot_all_spikes(han.sc)
            if visible(han.sortview.win)
                if !han.sc.pause

                else
                    #Refresh Screen
                    SpikeSorting.replot_sort(han.sortview)
                end
            end
            han.sc.buf.replot=false
        end
        if han.sc.thres_changed
            offline_thres_changed(han.sc,sorting)
        end
        if han.sc.show_thres
            SpikeSorting.plot_thres(han.sc)
        end
        if han.sc.rb_active
            SpikeSorting.draw_rb(han.sc)
        end
        if han.av.changed
            replot_analog(han)
        end

        reveal(han.sc.c2)
        reveal(han.sc.c3)

        sleep(0.01)
    end

end

function offline_thres_changed(sc,s)

    mythres=getproperty(sc.adj_thres,:value,Int)
    sc.thres=mythres

    #update sorting
    if (getproperty(sc.thres_widgets.all,:active,Bool))|(getproperty(sc.gain_widgets.all,:active,Bool))
        for i=1:length(s)
            s[i].thres=-1*sc.thres/sc.s #This isn't right, need a scale thing
        end
    else
        s[sc.spike].thres = -1 * sc.thres / sc.s
    end

    sc.thres_changed=false

    nothing
end

function analog_slider_cb(widget::Ptr,user_data::Tuple{Offline_GUI})

    han, = user_data

    han.av.pos = getproperty(han.av.adj,:value,Int)

    han.av.changed = true

    nothing
end

function analog_gain_cb(widget::Ptr,user_data::Tuple{Offline_GUI})

    han, = user_data

    han.av.changed = true

    nothing
end

function replot_analog(han)
    ctx = Gtk.getgc(han.av.c)

    set_source_rgb(ctx,0.0,0.0,0.0)
    paint(ctx)

    h=height(ctx)
    w=width(ctx)

    gain = han.sc.s

    set_source_rgb(ctx,1.0,1.0,1.0)
    move_to(ctx,1.0,h/2+han.v[han.av.pos]*gain)
    for i=2:round(Int64,w)
        line_to(ctx,i,h/2+han.v[han.av.pos+i-1]*gain)
    end
    stroke(ctx)

    reveal(han.av.c)

    han.av.changed = false
    nothing
end

ext: jl


text:
function equations_to_matrix_form(eqs, vars)
    expand_eqs = expand_derivatives.(eqs)
    A = Symbolics.jacobian(expand_eqs, vars; simplify = true)
    # # Ensure none of x is in A
    # vars_A = Symbolics.get_variables(A)
    # if !isempty(vars_A)
    #     for v in vars
    #         @assert all(vars_A .!= v)
    #     end
    # end

    b = simplify.(-substitute(expand_eqs, Dict(v => 0 for v in vars)))
    A, b
end

function generate_equations_of_motion()
    @variables t
    @variables u(t) w(t) θ′(t)

    @variables θ_t(t) θ_k(t) γ_OA(t)

    LT = symbolic_kinetic_energy(u, w, θ′, t, θ_t = θ_t, θ_k = θ_k, γ_OA = γ_OA)

    D = Differential(t)
    # Need to generate A and b for Ax = b where x is D(u), D(w), D(θ′)

    # To me A = jacobian of LT w.r.t x assuming LT is linear in x
    x = [D(u), D(w), D(θ′)]
    A, b = equations_to_matrix_form(LT, x)
    println("Variables in A and b after equations_to_matrix_form")
    println(Symbolics.get_variables.(A))
    println(Symbolics.get_variables.(b))

    generated_dir = joinpath(@__DIR__, "..", "generated")
    funA, _ = build_function(A, (θ_t, θ_k, γ_OA))
    write(joinpath(generated_dir, "inertia_lhs.jl"), string(funA))

    funb, _ = build_function(b, (u, w, θ′, #=t,=# θ_t, D(θ_t), D(D(θ_t)), θ_k, D(θ_k), D(D(θ_k)), γ_OA, D(γ_OA), D(D(γ_OA))))
    write(joinpath(generated_dir, "inertia_rhs.jl"), string(funb))
    (funA, funb)
end


text:
# Discretization methods as structures
# Full Discretization
abstract type DiscretizationMethod{fT} end
struct FullDiscretization{fT} <: DiscretizationMethod{fT}
    Δt::fT; # Time resolution of the discretisation
end
terms(::FullDiscretization) = 1
methodorder(::FullDiscretization) = 0

# Semi-Discretization
abstract type SemiDiscretizationMethodOrder end
struct NumericSD{N} <: SemiDiscretizationMethodOrder end
struct AnalyticZero <: SemiDiscretizationMethodOrder end
struct AnalyticZeroImproved <: SemiDiscretizationMethodOrder end
struct AnalyticFirst <: SemiDiscretizationMethodOrder end

struct SemiDiscretization{T,fT<:AbstractFloat} <: DiscretizationMethod{fT}
    Δt::fT; # Time resolution of the discretisation
end
SemiDiscretization(Δt::fT) where fT<:AbstractFloat = SemiDiscretization{NumericSD{1},fT}(Δt)
SemiDiscretization(SDorder::Int,Δt::fT) where fT<:AbstractFloat = SemiDiscretization{NumericSD{SDorder},fT}(Δt)
# n-th order semi discretisation for analytical calculation of the submatrices
SemiDiscretization(SDorder::Type{<:SemiDiscretizationMethodOrder},Δt::fT) where fT<:AbstractFloat = SemiDiscretization{SDorder,fT}(Δt)

_terms(::Type{AnalyticZero}) = 1
_terms(::Type{AnalyticZeroImproved}) = 2
_terms(::Type{AnalyticFirst}) = 2
_terms(::Type{NumericSD{N}}) where N = N+1

_methodorder(::Type{AnalyticZero}) = 0
_methodorder(::Type{AnalyticZeroImproved}) = 0
_methodorder(::Type{AnalyticFirst}) = 1
_methodorder(::Type{NumericSD{N}}) where N = N

terms(::SemiDiscretization{T,fT}) where {T<:SemiDiscretizationMethodOrder,fT} = _terms(T)
methodorder(::SemiDiscretization{T,fT}) where {T<:SemiDiscretizationMethodOrder,fT} = _methodorder(T)
# n-th order semi discretisation for numerical calculation of the submatrices


text:
module fenkar_pde

export fenkar!, p, u0, tspan, N

const N = Int(round(3/0.015 * sqrt(2)))
const D1= 0.0010/(0.015*0.015)
const D2= 0.0002/(0.015*0.015)

# define fenkar system
function H(x;k=100.0)
	return 0.5*(1.0 + tanh(k*x))
end
function fenkar!(dx, x, p, t)

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi,t0,TI,IA,θ = p[1:17]

	@inbounds for n=1:N

		# fenkar dynamics
		dx[n,1] =-(x[n,1]*H(uc-x[n,1])/to + H(x[n,1]-uc)/tr - x[n,2]*H(x[n,1]-uc)*(1.0-x[n,1])*(x[n,1]-uc)/td - x[n,3]*H(x[n,1]-ucsi;k=xk)/tsi)
		dx[n,2] = H(uc-x[n,1])*(1.0-x[n,2])/(tv1m*H(uv-x[n,1]) + tv2m*H(x[n,1]-uv)) - H(x[n,1]-uc)*x[n,2]/tvp
		dx[n,3] = H(uc-x[n,1])*(1.0-x[n,3])/twm - H(x[n,1]-uc)*x[n,3]/twp

		# no-flux with angular dependence on two diffusions
		if n==1
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(2*x[2,1] - 2*x[n,1])
		elseif n==N
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(2*x[N-1,1] - 2*x[n,1])
		else
			dx[n,1] = dx[n,1] + (D1*cos(θ)^2 + D2*sin(θ)^2)*(x[n-1,1] + x[n+1,1] - 2*x[n,1])
		end
	end

	@inbounds for n=1:10
		dx[n,1] = dx[n,1] + H(t-t0;k=1.0).*IA*sin(pi*(t-t0)/TI)^500
	end

	return nothing
end

# fenkar system parameters in p[1:13] (BR; from https://doi.org/10.1063/1.166311)
# and then p[14:14] = [fiber_angle]
# and then p[15:17] = [offset_time, stimulus_period, stimulus_amplitude]
   p = [29.0, 19.6, 1250.0, 3.33, 41.0, 870.0, 0.25, 12.5, 33.3, 10.0, 0.13, 0.04, 0.85, 0.0,   125.0, 0.1, 0.0]
   tspan = (0.0, 3000.0)
   u0 = zeros(Float64,N,3)
end


text:
# use BigFloats to determine the correct multipart significand values

setprecision(BigFloat, 1024-32)

if VERSION < v"0.7.0-"
    function randbf(n::Int, minpow2::Int=-30, maxpow2::Int=30)
        map(BigFloat, [
  "1.07667362842634302875882485034815716810629194710565656562969414207416239876804138664200409610113371132280829829629832263256714528995802234040849813704082672156927957863749387076287840301964688430754091065235913161629140290884101760484736726473976079401401147883477923675027520351735075797152855836416e+07",
  "7.60748557393455339828540112453332418619638860305803783802730707022859084587294339903810351152073795054157487281375223707799628501192534101282529017424792511150868032222132483546842951599762092688574875921051939584637424674000900494888159449760594598774242024183935370025516955967800612807273226203056e+03",
 "1.01009783883214796939239044340748787035775744586468925658997362182283343256995923331447855073801654281923815059649089988442055669719130377931811669893494753785172259260668401827608830443258153559955378038315877091487515417840929051078749844294237197559146430830186543200881912955210567697159734161746e+08",
 "4.1245214426363132106952925327575152946135655968460257705895100313409437154706034721072332923408364147795689183110749245398246334615815545526637837203883552374655798754154685643020718698704988357219715592219001016139798058056217494048704150413911515348170194880658304635415197358315798372069765192874e-08" ,
 "1.36992681634497157710823757218212462020923420764304266195733632579164888111419245630311779279748305397282610606960088714843091998438863291788182556221333846649556649520326334458000186369228808757627428735457286282838275315266930380123121254688334794848348755290800846374821276652967840753347718635029e-01",
 "2.94858317848686133293541130602062359010070703806418479610839760302786976368324411900869196286003468069815366567021751663397576055163808821822822097711517381468473106784563064985606637230351326570627909216412824526296346290474319021779700377267907116012148016739873785794441273223196937356937973298659e+02",
 "5.95055797578493207515924145850964017684622369150554788389091946765998252814685042346088932099000365451928691359946290283888488224325538350111090826576924308172152944045983054525124209794731824676188705811652720593410165071875328719162103666763195568281097309793543738523335250616935076596622638347407e+03",
 "3.15054948211191356314913169340959072089219610008297916033091017761704445253825251502165172463341604996908310454459355144649818087578185104293609662661839457745607545427806251391024523966288380842934561506016775114907752505720033212254076073607340345519878814733638203885238055887611157301459041318178e-05",
 "5.77174439114390945273223780323544620768413138937749434744711987997569297056682895255977202079242364047162306732488881077734196226219633554035583634986289034666025443917870332034794257064486275571281320317507018965620420900645995060677533252133506595149704703780132544755752158806504979663672033758911e-04",
 "4.25077457651057377879234217343886357354756446406779144818039083060523697061489328310373443441491281048311726114644057659947258050134289901253358963322154534095999369474744143070279038366242112550530139824644617376434795666898486123044574657401741178206605053191495847432137431639158207948239741121763e-07",
 "1.55682034004501874420112164333526085835012389073718535377079595939390816942886016384758068266677734217122812533071474821647671109235784675493592896971362614284072642114295321294088531594288947767273292078177131845531557139280374027532884522406602667434847367669450902383906530030725605802346291757329e+04",
 "1.14813889416809952394952796336233935453039197815740472314590980781511494565162088199218683863701886151913585902440224716610351397944534001776642148933018925885943822294877766583906263567114594697876001150704382796054810214591847356712635386930139623431463531365908887564165299926472475458986694936992"    ,
 "1.09736700932017264891502504112916215884978496395147385378024211548209007879430197010449264738980720175433779837070866624585059512848665987861336438350266374066628713795624485874888516535689087707302657229083418647724760413522167159656620818880631162538127686367509190385195187904208733155733917198565e+01",
 "4.73456459274572569812070971230583369816815626277683307633994713202237481923273540676481031080298232805624099440914043825555011082176793145710677507438923792398697743761679419998731749761221257578904492581543137680577929159687872144076997699538721156806148086106579722001805270108584960594025581846152e-01",
 "6.74968482387614057550330826710477012042628707503997951014418292047670613959235449126623600309141233050914872706417075420240044574773500648039320707780236936126504702480341633353248328084487055334556920416593820697678075832526726443921438576709327325827313454783411276976602317623583743430806741116158e-07",
 "3.21791225864409900109924218062718769949404895765446161061667938870057501681046915676247941906034526745941307181121408292398527357985090726690049720310000127458829866199757280416321334786265833264976108938476321783958201575418485860696556727064916445558253027138906003513871617900705218919073889971882e-08"
  ])
    end
else
    function randbf(n::Int, minpow2::Int=-30, maxpow2::Int=30)
        pow2s = rand(minpow2:maxpow2, n)
        bfs = rand(BigFloat, n)
        bfls = [ldexp(frexp(bfs[i])[1], pow2s[i]) for i=1:n]
        return bfls
    end
end

randfl(::Type{T}, n::Int, minpow2::Int=-30, maxpow2::Int=30)  where {T<:AbstractFloat} =
    map(T, randbf(n, minpow2, maxpow2))

const Floats1 = NTuple{1,T} where {T<:AbstractFloat}
const Floats2 = NTuple{2,T} where {T<:AbstractFloat}
const Floats3 = NTuple{3,T} where {T<:AbstractFloat}
const Floats4 = NTuple{4,T} where {T<:AbstractFloat}
const Floats5 = NTuple{5,T} where {T<:AbstractFloat}
const Floats6 = NTuple{6,T} where {T<:AbstractFloat}
const Floats7 = NTuple{7,T} where {T<:AbstractFloat}
const Floats8 = NTuple{8,T} where {T<:AbstractFloat}
const Floats9 = NTuple{9,T} where {T<:AbstractFloat}

function Floats1(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   return (fl1, )
end

function Floats2(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   return fl1, fl2
end

function Floats3(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   return fl1, fl2, fl3
end

function Floats4(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   fl4 = T(bf - fl1 - fl2 - fl3)
   return fl1, fl2, fl3, fl4
end

function Floats5(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   fl4 = T(bf - fl1 - fl2 - fl3)
   fl5 = T(bf - fl1 - fl2 - fl3 - fl4)
   return fl1, fl2, fl3, fl4, fl5
end

function Floats6(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   fl4 = T(bf - fl1 - fl2 - fl3)
   fl5 = T(bf - fl1 - fl2 - fl3 - fl4)
   fl6 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5)
   return fl1, fl2, fl3, fl4, fl5, fl6
end

function Floats7(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   fl4 = T(bf - fl1 - fl2 - fl3)
   fl5 = T(bf - fl1 - fl2 - fl3 - fl4)
   fl6 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5)
   fl7 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5 - fl6)
   return fl1, fl2, fl3, fl4, fl5, fl6, fl7
end

function Floats8(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   fl4 = T(bf - fl1 - fl2 - fl3)
   fl5 = T(bf - fl1 - fl2 - fl3 - fl4)
   fl6 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5)
   fl7 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5 - fl6)
   fl8 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5 - fl6 - fl7)
   return fl1, fl2, fl3, fl4, fl5, fl6, fl7, fl8
end

function Floats9(::Type{T}, bf::BigFloat) where {T<:AbstractFloat}
   fl1 = T(bf)
   fl2 = T(bf - fl1)
   fl3 = T(bf - fl1 - fl2)
   fl4 = T(bf - fl1 - fl2 - fl3)
   fl5 = T(bf - fl1 - fl2 - fl3 - fl4)
   fl6 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5)
   fl7 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5 - fl6)
   fl8 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5 - fl6 - fl7)
   fl9 = T(bf - fl1 - fl2 - fl3 - fl4 - fl5 - fl6 - fl7 - fl8)
   return fl1, fl2, fl3, fl4, fl5, fl6, fl7, fl8, fl9
end


text:
module ObservationSchemes

    using LinearAlgebra, StaticArrays
    using GaussianDistributions, Random
    using Distributions, OrderedCollections

    import GaussianDistributions: whiten, unwhiten

    # these are used throughout the suite, they need to be put into one place only
    # and imported from there.
    ismutable(el) = ismutable(typeof(el))
    ismutable(::Type) = Val(false)
    ismutable(::Type{<:Array}) = Val(true)

    include("starting_point_priors.jl")
    include("first_passage_times.jl")
    include("observations_common.jl")
    include("observations_linear_gaussian.jl")
    include("observations_general.jl")
    include("all_observations.jl")
    include("all_observations_utilities.jl")
    include("obs_scheme.jl")
    include("convenience_functions.jl")

    # starting_point_priors.jl
    export KnownStartingPt, GsnStartingPt

    # first_passage_times.jl
    export NoFirstPassageTimes, FirstPassageTimeInfo

    # observations_lienar_gaussian.jl
    export LinearGsnObs, fpt_info

    # observations_general.jl
    export GeneralObs

    # all_observations.jl
    export AllObservations, add_recording!, add_recordings!, add_dependency!
    export initialize

    # obs_scheme.jl
    export ObsScheme, load_data

    # convenience_functions.jl
    export package, build_recording, num_recordings, num_obs

    export print_parameters

end # module


text:
module BeaconK8sUtilities

using ProgressMeter, JSON3, kubectl_jll
using Preferences
using Mustache
using RelocatableFolders
using Compat
using FilePathsBase
using LoggingExtras
using LoggingFormats
using Logging, Dates
using Distributed # for `myid`

const TEMPLATES = @path joinpath(@__DIR__, "..", "templates")

export kubectl

export get_status, last_condition, wait_until_pod_ready, watch_logs, port_forward,
       get_pod_names, get_current_namespace, follow
include("utilities.jl")

export default_ecr, default_service_account
include("preferences.jl")

export setup_tensorboard, setup_follow
include("templates.jl")

export json_logger
include("logger.jl")

end # module


text:
using Bokeh
# simple example with debug true so we can read the json more easily
Bokeh.debug(true)
plotfile("debug.html")
y = collect(1:10).^2
# here we don't define x so it's auto generated as 0:(length(y) - 1)
plot(y, autoopen=true)


text:
mutable struct RegenerateContext{C <: CallSite,
                                 T <: Trace,
                                 K <: UnconstrainedSelection,
                                 P <: Parameters,
                                 D <: Diff} <: ExecutionContext
    prev::C
    tr::T
    select::K
    weight::Float64
    score::Float64
    discard::HierarchicalTrace
    visited::Visitor
    params::P
    argdiffs::D
    RegenerateContext(cl::C, select::K, argdiffs::D) where {C <: CallSite, K <: UnconstrainedSelection, D <: Diff} = new{C, typeof(cl.trace), K, EmptyParameters, D}(cl, typeof(cl.trace)(), select, 0.0, 0.0, Trace(), Visitor(), Parameters(), argdiffs)
    RegenerateContext(cl::C, select::K, params::P, argdiffs::D) where {C <: CallSite, K <: UnconstrainedSelection, P <: Parameters, D <: Diff} = new{C, typeof(cl.trace), K, P, D}(cl, typeof(cl.trace)(), select, 0.0, 0.0, Trace(), Visitor(), params, argdiffs)
end
Regenerate(cl, select, argdiffs) = RegenerateContext(cl, select, argdiffs)

# Regenerate has a special dynamo.
@dynamo function (mx::RegenerateContext)(a...)
    ir = IR(a...)
    ir == nothing && return
    recur!(ir)
    return ir
end

# ------------ Convenience ------------ #

function regenerate(ctx::RegenerateContext, bbcs::HierarchicalCallSite, new_args...)
    ret = ctx(bbcs.fn, new_args...)
    return ret, HierarchicalCallSite(ctx.tr, ctx.score, bbcs.fn, new_args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function regenerate(sel::L, bbcs::HierarchicalCallSite) where L <: UnconstrainedSelection
    argdiffs = NoChange()
    ctx = RegenerateContext(bbcs, sel, argdiffs)
    return regenerate(ctx, bbcs, bbcs.args...)
end

function regenerate(sel::L, params, bbcs::HierarchicalCallSite) where L <: UnconstrainedSelection
    argdiffs = NoChange()
    ctx = RegenerateContext(bbcs, sel, params, argdiffs)
    return regenerate(ctx, bbcs, bbcs.args...)
end

function regenerate(sel::L, bbcs::HierarchicalCallSite, argdiffs::D, new_args...) where {L <: UnconstrainedSelection, D <: Diff}
    ctx = RegenerateContext(bbcs, sel, argdiffs)
    return regenerate(ctx, bbcs, new_args...)
end

function regenerate(sel::L, params, bbcs::HierarchicalCallSite, argdiffs::D, new_args...) where {L <: UnconstrainedSelection, D <: Diff}
    ctx = RegenerateContext(bbcs, sel, params, argdiffs)
    return regenerate(ctx, bbcs, new_args...)
end

# TODO: fix for dispatch with params.
function regenerate(sel::L, vcs::VectorizedCallSite{typeof(plate)}) where {L <: UnconstrainedSelection, D <: Diff}
    argdiffs = NoChange()
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = RegenerateContext(vcs, v_sel, argdiffs)
    ret = ctx(plate, addr, vcs.fn, vcs.args)
    return ret, VectorizedCallSite{typeof(plate)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function regenerate(sel::L, vcs::VectorizedCallSite{typeof(markov)}) where {L <: UnconstrainedSelection, D <: Diff}
    argdiffs = NoChange()
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = RegenerateContext(vcs, v_sel, argdiffs)
    ret = ctx(markov, addr, vcs.fn, vcs.args[1], vcs.args[2]...)
    return ret, VectorizedCallSite{typeof(markov)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function regenerate(sel::L, vcs::VectorizedCallSite{typeof(markov)}, d::NoChange, len::Int) where {L <: UnconstrainedSelection, D <: Diff}
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = RegenerateContext(vcs, v_sel, d)
    ret = ctx(markov, addr, vcs.fn, len, vcs.args[2]...)
    return ret, VectorizedCallSite{typeof(markov)}(ctx.tr, ctx.score, vcs.fn, vcs.args, ret), ctx.weight, UndefinedChange(), ctx.discard
end

function regenerate(sel::L, vcs::VectorizedCallSite{typeof(markov)}, len::Int) where {L <: UnconstrainedSelection, D <: Diff}
    return regenerate(sel, vcs, NoChange(), len)
end

# ------------ includes ------------ #

include("hierarchical/regenerate.jl")
include("plate/regenerate.jl")
include("markov/regenerate.jl")
include("factor/regenerate.jl")

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct RegenerateContext{T <: Trace,
                                 L <: UnconstrainedSelection,
                                 P <: Parameters} <: ExecutionContext
    prev::T
    tr::T
    select::L
    weight::Float64
    score::Float64
    discard::T
    visited::Visitor
    params::P
end
```

Inner constructors:

```julia
function RegenerateContext(tr::T, sel::Vector{Address}) where T <: Trace
    un_sel = selection(sel)
    new{T, typeof(un_sel), EmptyParameters}(tr, Trace(), un_sel, 0.0, Trace(), Visitor(), Parameters())
end
function RegenerateContext(tr::T, sel::L) where {T <: Trace, L <: UnconstrainedSelection}
    new{T, L, EmptyParameters}(tr, Trace(), sel, 0.0, Trace(), Visitor(), Parameters())
end
```

Outer constructors:

```julia
Regenerate(tr::Trace, sel::Vector{Address}) = RegenerateContext(tr, sel)
Regenerate(tr::Trace, sel::UnconstrainedSelection) = RegenerateContext(tr, sel)
```

The `RegenerateContext` is used for MCMC algorithms, to propose new choices for addresses indicated by an `UnconstrainedSelection` in the `select` field.
""", RegenerateContext)

@doc(
"""
```julia
ret, cl = regenerate(sel::L, bbcs::HierarchicalCallSite, new_args...) where L <: UnconstrainedSelection
ret, cl = regenerate(sel::L, bbcs::HierarchicalCallSite) where L <: UnconstrainedSelection
```
`regenerate` is an API to the `RegenerateContext` execution context. `regenerate` requires that users provide an `UnconstrainedSelection`, an original call site, and possibly a set of new arguments to be used in the regeneration step. This context internally keeps track of the bookkeeping required to increment likelihood weights, as well as prune off parts of the trace which are invalid if a regenerated choice changes the shape of the trace (e.g. control flow), and returns a new return value `ret` as well as the modified call site `cl`.
""", regenerate)


text:
"""
constraints.jl
"""

using Altro
using LinearAlgebra

@inline bound_constraint_eval(x, u, x_max, x_min, u_max, u_min) = ([
    x - x_max;
    x_min - x;
    u - u_max;
    u_min - u;
])

function bound_constraint_jacobian(x, u, x_max, x_min, u_max, u_min)
    n = length(x)
    m = length(u)
    p = 2 * n + 2 * m
    Cx = zeros(p, n)
    for i = 1:n
        Cx[i, i] = 1.
        Cx[i + n, i] = -1.
    end
    Cu = zeros(p, m)
    for i = 1:m
        Cu[i + 2 * n, i] = 1.
        Cu[i + 2 * n + m, i] = -1.
    end
    return Cx, Cu
end

function test_bound_constraint()
    big = 1
    big_ = 1000
    n = 2
    m = 2
    M(x) = x
    V(x) = x
    for i = 1:big
        x_max = rand(n)
        x_min = -rand(n)
        u_max = rand(m)
        u_min = -rand(m)
        bc = BoundConstraint(x_max, x_min, u_max, u_min, n, m, M, V)
        p = 2 * n + 2 * m
        c = zeros(p)
        Cx = zeros(p, n)
        Cu = zeros(p, m)
        for j = 1:big_
            x = 2 .* rand(n) .- 1
            X = [x]
            u = 2 .* rand(m) .- 1
            U = [u]
            k = 1
            Altro.evaluate!(c, bc, X, U, k)
            c_ = bound_constraint_eval(x, u, x_max, x_min, u_max, u_min)
            @assert c ≈ c_
            Altro.jacobian!(Cx, Cu, bc, X, U, k)
            Cx_, Cu_ = bound_constraint_jacobian(x, u, x_max, x_min, u_max, u_min)
            @assert Cx ≈ Cx_
            @assert Cu ≈ Cu_
        end
    end
end


text:
using Revise
using SIREN
using Test
using Flux: params

function test_uniform()
    a = []

    for i in SIREN.uniform(100,100)
        if i < -1 || i > 1
            push!(a, i)
        end
    end
    @test length(a) == 0
end

function test_dense()
    fan_in = 5
    fan_out = 3
    W = uniform(fan_out, fan_in)
    x = uniform(fan_in)
    b = uniform(fan_out)
    omega_0 = 30

    d = SIREN.Dense(omega_0, W, b)

    @test params(d)[1] == W
    @test params(d)[2] == b

    @test d(x) == sin.(omega_0 * (W * x .+ b))
end

function test_conv()
    m = SIREN.Conv((3,2), 4=>5)
    d = rand(28,28,4,2)
    m(d)
end

function test_CuArrays()

end

@testset "SIREN.jl" begin
    test_uniform()
    test_dense()
    test_conv()
end

max_issues_repo_name: MacKenzieHnC/SIREN

max_forks_repo_name: MacKenzieHnC/SIREN


text:
using GLVisualize, GLAbstraction, GeometryTypes, Reactive

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end

description = """
Surface mesh with arbitrary x,y,z coordinates.
"""

function mgrid(dim1, dim2)
    X = [i for i in dim1, j in dim2]
    Y = [j for i in dim1, j in dim2]
    return X,Y
end

rotation_angle = const_lift(*, timesignal, 2f0*pi)
rotation = map(rotationmatrix_z, rotation_angle)

dphi, dtheta = pi/200.0f0, pi/200.0f0

phi,theta = mgrid(0f0:dphi:(pi+dphi*1.5f0), 0f0:dtheta:(2f0*pi+dtheta*1.5f0));
m0 = 4f0; m1 = 3f0; m2 = 2f0; m3 = 3f0; m4 = 6f0; m5 = 2f0; m6 = 6f0; m7 = 4f0;
a = sin.(m0*phi).^m1;
b = cos.(m2*phi).^m3;
c = sin.(m4*theta).^m5;
d = cos.(m6*theta).^m7;
r = a + b + c + d;
x = r.*sin.(phi).*cos.(theta);
y = r.*cos.(phi);
z = r.*sin.(phi).*sin.(theta);

surface = visualize((x,y,z), :surface, model=rotation)
_view(surface, window)

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


text:
# This file replicates the testsuite provided by MINPACK for the hybrj
# function. It consists in 14 functions with their analytical Jacobian.
#
# We run exactly the same examples than in MINPACK, with the following
# exceptions:
# - Chebyquad in dimension 8 is not run. It does not converge (it is not either
#   with MINAPCK).
# - Only one initial point is tested for each problem. For some problems,
#   MINPACK would test also other initial points which are 10 or 100 times the
#   original point.

# If the results should be printed to a
PRINT_FILE = false

@testset "minpack" begin

function rosenbrock()
    function f!(fvec, x)
        fvec[1] = 1 - x[1]
        fvec[2] = 10(x[2]-x[1]^2)
    end
    function j!(fjac, x)
        fjac[1,1] = -1
        fjac[1,2] = 0
        fjac[2,1] = -20x[1]
        fjac[2,2] = 10
    end
    (OnceDifferentiable(f!, j!, [-1.2, 1.0], [-1.2, 1.0]), [-1.2; 1.0], "Rosenbrock")
end

function powell_singular()
    function f!(fvec, x)
        fvec[1] = x[1] + 10x[2]
        fvec[2] = sqrt(5)*(x[3] - x[4])
        fvec[3] = (x[2] - 2x[3])^2
        fvec[4] = sqrt(10)*(x[1] - x[4])^2
    end
    function j!(fjac, x)
        fill!(fjac, 0)
        fjac[1,1] = 1
        fjac[1,2] = 10
        fjac[2,3] = sqrt(5)
        fjac[2,4] = -fjac[2,3]
        fjac[3,2] = 2(x[2] - 2x[3])
        fjac[3,3] = -2fjac[3,2]
        fjac[4,1] = 2sqrt(10)*(x[1] - x[4])
        fjac[4,4] = -fjac[4,1]
    end
    (OnceDifferentiable(f!, j!, rand(4), [3.; -1.; 0.; 1.]), [3.; -1.; 0.; 1.], "Powell singular")
end

function powell_badly_scaled()
    c1 = 1e4
    c2 = 1.0001
    function f!(fvec, x)
        fvec[1] = c1*x[1]*x[2] - 1
        fvec[2] = exp(-x[1]) + exp(-x[2]) - c2
    end
    function j!(fjac, x)
        fjac[1,1] = c1*x[2]
        fjac[1,2] = c1*x[1]
        fjac[2,1] = -exp(-x[1])
        fjac[2,2] = -exp(-x[2])
    end
    (OnceDifferentiable(f!, j!, rand(2), [0.; 1.]), [0.; 1.], "Powell badly scaled")
end

function wood()
    c3 = 2e2
    c4 = 2.02e1
    c5 = 1.98e1
    c6 = 1.8e2

    function f!(fvec, x)
        temp1 = x[2] - x[1]^2
        temp2 = x[4] - x[3]^2
        fvec[1] = -c3*x[1]*temp1 - (1 - x[1])
        fvec[2] = c3*temp1 + c4*(x[2] - 1) + c5*(x[4] - 1)
        fvec[3] = -c6*x[3]*temp2 - (1 - x[3])
        fvec[4] = c6*temp2 + c4*(x[4] - 1) + c5*(x[2] - 1)
    end

    function j!(fjac, x)
        fill!(fjac, 0)
        temp1 = x[2] - 3x[1]^2
        temp2 = x[4] - 3x[3]^2
        fjac[1,1] = -c3*temp1 + 1
        fjac[1,2] = -c3*x[1]
        fjac[2,1] = -2*c3*x[1]
        fjac[2,2] = c3 + c4
        fjac[2,4] = c5
        fjac[3,3] = -c6*temp2 + 1
        fjac[3,4] = -c6*x[3]
        fjac[4,2] = c5
        fjac[4,3] = -2*c6*x[3]
        fjac[4,4] = c6 + c4
    end
    (OnceDifferentiable(f!, j!, rand(4), [-3.; -1; -3; -1]), [-3.; -1; -3; -1], "Wood")
end

function helical_valley()
    tpi = 8*atan(1)
    c7 = 2.5e-1
    c8 = 5e-1

    function f!(fvec, x)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        fvec[1] = 10(x[3] - 10*temp1)
        fvec[2] = 10(temp2 - 1)
        fvec[3] = x[3]
    end

    function j!(fjac, x)
        temp = x[1]^2 + x[2]^2
        temp1 = tpi*temp
        temp2 = sqrt(temp)
        fjac[1,1] = 100*x[2]/temp1
        fjac[1,2] = -100*x[1]/temp1
        fjac[1,3] = 10
        fjac[2,1] = 10*x[1]/temp2
        fjac[2,2] = 10*x[2]/temp2
        fjac[2,3] = 0
        fjac[3,1] = 0
        fjac[3,2] = 0
        fjac[3,3] = 1
    end
    (OnceDifferentiable(f!, j!, rand(3), [-1.; 0; 0]), [-1.; 0; 0], "Helical Valley")
end

function watson(n::Integer)
    c9 = 2.9e1

    function f!(fvec, x)
        fill!(fvec, 0)
        for i = 1:29
            ti = i/c9
            sum1 = 0.0
            temp = 1.0
            for j = 2:n
                sum1 += (j-1)*temp*x[j]
                temp *= ti
            end
            sum2 = 0.0
            temp = 1.0
            for j = 1:n
                sum2 += temp*x[j]
                temp *= ti
            end
            temp1 = sum1 - sum2^2 - 1
            temp2 = 2*ti*sum2
            temp = 1/ti
            for k = 1:n
                fvec[k] += temp*(k-1-temp2)*temp1
                temp *= ti
            end
        end
        temp = x[2] - x[1]^2 - 1
        fvec[1] += x[1]*(1-2temp)
        fvec[2] += temp
    end

    function j!(fjac, x)
        fill!(fjac, 0)
        for i = 1:29
            ti = i/c9
            sum1 = 0.0
            temp = 1.0
            for j = 2:n
                sum1 += (j-1)*temp*x[j]
                temp *= ti
            end
            sum2 = 0.0
            temp = 1.0

            for j = 1:n
                sum2 += temp*x[j]
                temp *= ti
            end
            temp1 = 2(sum1 - sum2^2 - 1)
            temp2 = 2*sum2
            temp = ti^2
            tk = 1.0

            for k = 1:n
                tj = tk
                for j = k:n
                    fjac[k,j] += tj*((Float64(k-1)/ti - temp2)*(Float64(j-1)/ti - temp2) - temp1)
                    tj *= ti
                end
                tk *= temp
            end
        end
        fjac[1,1] += 6x[1]^2 - 2x[2] + 3
        fjac[1,2] -= 2x[1]
        fjac[2,2] += 1
        for k = 1:n
            for j = k:n
                fjac[j,k] = fjac[k,j]
            end
        end
    end
    (OnceDifferentiable(f!, j!, zeros(n), zeros(n)), zeros(n), "Watson")
end

function chebyquad(n::Integer)
    tk = 1/n

    function f!(fvec, x)
        fill!(fvec, 0)
        for j = 1:n
            temp1 = 1.0
            temp2 = 2x[j]-1
            temp = 2temp2
            for i = 1:n
                fvec[i] += temp2
                ti = temp*temp2 - temp1
                temp1 = temp2
                temp2 = ti
            end
        end
        iev = -1.0
        for k = 1:n
            fvec[k] *= tk
            if iev > 0
                fvec[k] += 1/(k^2-1)
            end
            iev = -iev
        end
    end

    function j!(fjac, x)
        for j = 1:n
            temp1 = 1.
            temp2 = 2x[j] - 1
            temp = 2*temp2
            temp3 = 0.0
            temp4 = 2.0
            for k = 1:n
                fjac[k,j] = tk*temp4
                ti = 4*temp2 + temp*temp4 - temp3
                temp3 = temp4
                temp4 = ti
                ti = temp*temp2 - temp1
                temp1 = temp2
                temp2 = ti
            end
        end
    end
    (OnceDifferentiable(f!, j!, collect(1:n)/(n+1), collect(1:n)/(n+1)), collect(1:n)/(n+1), "Chebyquad")
end

function brown_almost_linear(n::Integer)
    function f!(fvec, x)
        sum1 = sum(x) - (n+1)
        for k = 1:(n-1)
            fvec[k] = x[k] + sum1
        end
        fvec[n] = prod(x) - 1
    end

    function j!(fjac, x)
        fill!(fjac, 1)
        fjac[diagind(fjac)] .= 2
        prd = prod(x)
        for j = 1:n
            if x[j] == 0.0
                fjac[n,j] = 1.
                for k = 1:n
                    if k != j
                        fjac[n,j] *= x[k]
                    end
                end
            else
                fjac[n,j] = prd/x[j]
            end
        end
    end
    (OnceDifferentiable(f!, j!, 0.5*ones(n), 0.5*ones(n)), 0.5*ones(n), "Brown almost-linear")
end

function discrete_boundary_value(n::Integer)
    h = 1/(n+1)

    function f!(fvec, x)
        for k = 1:n
            temp = (x[k] + k*h + 1)^3
            if k != 1
                temp1 = x[k-1]
            else
                temp1 = 0.0
            end
            if k != n
                temp2 = x[k+1]
            else
                temp2 = 0.0
            end
            fvec[k] = 2x[k] - temp1 - temp2 + temp*h^2/2
        end
    end

    function j!(fjac, x)
        for k = 1:n
            temp = 3*(x[k]+k*h+1)^2
            for j = 1:n
                fjac[k,j] = 0
            end
            fjac[k,k] = 2 + temp*h^2/2
            if k != 1
                fjac[k,k-1] = -1
            end
            if k != n
                fjac[k,k+1] = -1
            end
        end
    end
    initial_x = collect(1:n)*h
    initial_x = initial_x .* (initial_x .- 1)

    (OnceDifferentiable(f!, j!, initial_x, initial_x), initial_x, "Discrete boundary value")
end

function discrete_integral_equation(n::Integer)
    h = 1/(n+1)

    function f!(fvec, x)
        for k = 1:n
            tk = k*h
            sum1 = 0.0
            for j = 1:k
                tj = j*h
                sum1 += tj*(x[j] + tj + 1)^3
            end
            sum2 = 0.0
            kp1 = k+1
            if n >= kp1
                for j = kp1:n
                    tj = j*h
                    sum2 += (1-tj)*(x[j] + tj + 1)^3
                end
            end
            fvec[k] = x[k] + h*((1-tk)*sum1 + tk*sum2)/2
        end
    end

    function j!(fjac, x)
        for k = 1:n
            tk = k*h
            for j = 1:n
                tj = j*h
                fjac[k,j] = h*min(tj*(1-tk), tk*(1-tj))*3(x[j] + tj + 1)^2/2
            end
            fjac[k,k] += 1
        end
    end

    initial_x = collect(1:n)*h
    initial_x = initial_x .* (initial_x .- 1)

    (OnceDifferentiable(f!, j!, initial_x, initial_x), initial_x, "Discrete integral equation")
end

function trigonometric(n::Integer)
    function f!(fvec, x)
        for j = 1:n
            fvec[j] = cos(x[j])
        end
        sum1 = sum(fvec)
        for k = 1:n
            fvec[k] = n+k-sin(x[k]) - sum1 - k*fvec[k]
        end
    end

    function j!(fjac, x)
        for j = 1:n
            temp = sin(x[j])
            for k = 1:n
                fjac[k,j] = temp
            end
            fjac[j,j] = (j+1)*temp - cos(x[j])
        end
    end
    (OnceDifferentiable(f!, j!,ones(n)/n, ones(n)/n), ones(n)/n, "Trigonometric")
end

function variably_dimensioned(n::Integer)
    function f!(fvec, x)
        sum1 = 0.0
        for j = 1:n
            sum1 += j*(x[j]-1)
        end
        temp = sum1*(1+2sum1^2)
        for k = 1:n
            fvec[k] = x[k] - 1 + k*temp
        end
    end

    function j!(fjac, x)
        sum1 = 0.0
        for j = 1:n
            sum1 += j*(x[j]-1)
        end
        temp = 1 + 6sum1^2
        for k = 1:n
            for j = k:n
                fjac[k,j] = k*j*temp
                fjac[j,k] = fjac[k,j]
            end
            fjac[k,k] += 1
        end
    end
    (OnceDifferentiable(f!, j!, 1 .- collect(1:n)/n, 1 .- collect(1:n)/n), 1 .- collect(1:n)/n, "Variably dimensioned")
end

function broyden_tridiagonal(n::Integer)
    function f!(fvec, x)
        for k = 1:n
            temp = (3-2x[k])*x[k]
            if k != 1
                temp1 = x[k-1]
            else
                temp1 = 0.0
            end
            if k != n
                temp2 = x[k+1]
            else
                temp2 = 0.0
            end
            fvec[k] = temp - temp1 - 2temp2 + 1
        end
    end

    function j!(fjac, x)
        fill!(fjac, 0)
        for k = 1:n
            fjac[k,k] = 3-4x[k]
            if k != 1
                fjac[k,k-1] = -1
            end
            if k != n
                fjac[k,k+1] = -2
            end
        end
    end
    (OnceDifferentiable(f!, j!, -ones(n), -ones(n)), -ones(n), "Broyden tridiagonal")
end

function broyden_banded(n::Integer)
    ml = 5
    mu = 1

    function f!(fvec, x)
        for k = 1:n
            k1 = max(1, k-ml)
            k2 = min(k+mu, n)
            temp = 0.0
            for j = k1:k2
                if j != k
                    temp += x[j]*(1+x[j])
                end
            end
            fvec[k] = x[k]*(2+5x[k]^2) + 1 - temp
        end
    end

    function j!(fjac, x)
        fill!(fjac, 0)
        for k = 1:n
            k1 = max(1, k-ml)
            k2 = min(k+mu, n)
            for j = k1:k2
                if j != k
                    fjac[k,j] = -(1+2x[j])
                end
            end
            fjac[k,k] = 2+15x[k]^2
        end
    end
    (OnceDifferentiable(f!, j!, -ones(n), -ones(n)), -ones(n), "Broyden banded")
end

alltests = [rosenbrock();
            powell_singular();
            powell_badly_scaled();
            wood();
            helical_valley();
            watson(6);
            watson(9);
            chebyquad(5);
            chebyquad(6);
            chebyquad(7);
            #chebyquad(8);
            chebyquad(9);
            brown_almost_linear(10);
            brown_almost_linear(30);
            brown_almost_linear(40);
            discrete_boundary_value(10);
            discrete_integral_equation(1);
            discrete_integral_equation(10);
            trigonometric(10);
            variably_dimensioned(10);
            broyden_tridiagonal(10);
            broyden_banded(10);
           ]

TESTS_FAIL_NEWTON = ["Trigonometric", "Chebyquad", "Brown almost-linear"]
TESTS_FAIL_BROYDEN = ["Powell singular", "Wood", "Trigonometric", "Watson"]# ["Rosenbrock", "Powell singular", "Powell badly scaled", "Trigonometric", "Watson"]

if PRINT_FILE; f_out = open("minpack_results.dat", "w"); end

@printf("%-45s   %5s   %5s   %5s   %14s     %10s\n", "Function", "Dim", "NFEV",
        "NJEV", "Final inf-norm", "total time")
println("-"^86)

if PRINT_FILE
    @printf(f_out, "%-45s   %5s   %5s   %5s   %14s\n", "Function", "Dim", "NFEV",
            "NJEV", "Final inf-norm")
    println(f_out, "-"^86)
end

for (df, initial, name) in alltests
    for method in (:trust_region, :newton, :broyden)
        if method == :broyden
            obj = NonDifferentiable(df.f, initial, initial)
        else
            obj = OnceDifferentiable(df.f, initial, initial)
        end
        if method == :newton && name in TESTS_FAIL_NEWTON
            continue
        end
        if method == :broyden
            if name in TESTS_FAIL_BROYDEN
                continue
            end
            tot_time = @elapsed r = nlsolve(obj, initial, method = method)
        else
            tot_time = @elapsed r = nlsolve(obj, initial, method = method)
        end
        @printf("%-45s   %5d   %5d   %5d   %14e   %10e\n", name*"-"*string(method), length(initial),
        r.f_calls, r.g_calls, r.residual_norm, tot_time)
        @test converged(r)
        # with autodiff
        if !(method == :broyden)
            tot_time2 = @elapsed r_AD = nlsolve(obj.f, initial, method = method, autodiff = :forward)
            @printf("%-45s   %5d   %5d   %5d   %14e   %10e\n", name*"-"*string(method)*"-AD",
                    length(initial), r_AD.f_calls, r_AD.g_calls, r_AD.residual_norm, tot_time)
            @test converged(r_AD)
        end
        # Fix this to work with ndf
        if PRINT_FILE
            @printf(f_out, "%-45s   %5d   %5d   %5d   %14e\n", name*"-"*string(method)*"-AD", length(initial),
                    r_AD.f_calls, r_AD.g_calls, r_AD.residual_norm)
            @printf(f_out, "%-45s   %5d   %5d   %5d   %14e\n", name*"-"*string(method), length(initial),
                    r.f_calls, r.g_calls, r.residual_norm)
        end
    end
end

if PRINT_FILE; close(f_out); end

end


text:
using SafeTestsets, Pkg

const GROUP = get(ENV, "GROUP", "All")

function dev_subpkg(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.develop(PackageSpec(path=subpkg_path))
end

function activate_subpkg_env(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.activate(subpkg_path)
    Pkg.develop(PackageSpec(path=subpkg_path))
    Pkg.instantiate()
end

@time begin
if GROUP == "All" || GROUP == "Core"
    dev_subpkg("GalacticOptimJL")
    dev_subpkg("GalacticFlux")
    @safetestset "AD Tests" begin
        include("ADtests.jl")
    end
    @safetestset "Mini batching" begin
        include("minibatch.jl")
    end
    @safetestset "DiffEqFlux" begin
        include("diffeqfluxtests.jl")
    end
elseif GROUP == "GPU"
    dev_subpkg("GalacticOptimJL")
    dev_subpkg("GalacticFlux")
    activate_downstream_env()
    @safetestset "DiffEqFlux GPU" begin
        include("downstream/gpu_neural_ode.jl")
    end
else
    dev_subpkg(GROUP)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", GROUP)
    Pkg.test(PackageSpec(name=GROUP, path=subpkg_path))
end
end


text:
module Map_removeEQ

export map_removeEQ

using Dates, JLD2, DSP, PlotlyJS, Printf, ORCA, FileIO

using SeisIO

include("get_kurtosis.jl")
include("remove_eq.jl")

using .Get_kurtosis, .Remove_eq

"""
    ParallelEQremoval(dlid, InputDict::Dict)
    remove earthquake and save it into jld2 file.

"""
function map_removeEQ(dlid, InputDict::Dict)

    #store data

    finame                 =InputDict["finame"]
    IsKurtosisRemoval      =InputDict["IsKurtosisRemoval"]
    max_edgetaper_duration =InputDict["max_edgetaper_duration"]
    kurtosis_timewindow    =InputDict["kurtosis_timewindow"]
    kurtosis_threshold     =InputDict["kurtosis_threshold"]
    IsSTALTARemoval        =InputDict["IsSTALTARemoval"]
    stalta_longtimewindow  =InputDict["stalta_longtimewindow"]
    stalta_threshold       =InputDict["stalta_threshold"]
    invert_tukey_α         =InputDict["invert_tukey_α"]
    max_wintaper_duration  =InputDict["max_wintaper_duration"]
    removal_shorttimewindow=InputDict["removal_shorttimewindow"]
    overlap                =InputDict["overlap"]
    plot_kurtosis_α        =InputDict["plot_kurtosis_α"]
    plot_boxheight         =InputDict["plot_boxheight"]
    plot_span              =InputDict["plot_span"]
    fodir                  =InputDict["fodir"]
    foname                 =InputDict["foname"]
    fopath                 =InputDict["fopath"]
    IsSaveFig              =InputDict["IsSaveFig"]

    DLtimestamplist        =InputDict["DLtimestamplist"]
    stationlist            =InputDict["stationlist"]
    NumofTimestamp         =InputDict["NumofTimestamp"]

    tstamp = DLtimestamplist[dlid]

    if mod(dlid, round(0.1*NumofTimestamp)+1) == 0
        println(@sprintf("start process %s", tstamp))
    end

    SRall = SeisData(length(stationlist))

    icount = 0

    bt_getkurtosis = 0.0
    bt_removeeq = 0.0

    for st = stationlist
        #S = t[joinpath(tstamp, st)]
        st1 = replace(st, "-"=>"")
        S = FileIO.load(finame, joinpath(tstamp, st1))

        if S.misc["dlerror"] == 0
            dt = 1/S.fs
            #tvec = collect(0:S.t[2,1]-1) * dt ./ 60 ./ 60
            tvec = collect(0:length(S.x)-1) * dt ./ 60 ./ 60

            #tapering to avoid instrumental edge artifacts
            SeisIO.taper!(S,  t_max = max_edgetaper_duration, α=0.05)
            S1 = deepcopy(S)

            #set long window length to user input since last window of previous channel will have been adjusted
            S1.misc["eqtimewindow"] = fill(true, length(S1.x))

            if IsKurtosisRemoval
                # compute kurtosis and detect earthqukes
                bt_1 = @elapsed S1 = Get_kurtosis.get_kurtosis(S1, float(kurtosis_timewindow))
                bt_2 = @elapsed S1 = Remove_eq.detect_eq_kurtosis(S1, tw=float(removal_shorttimewindow), kurtosis_threshold=float(kurtosis_threshold), overlap=float(overlap))

                btsta_1 = 0

                if IsSTALTARemoval
                    # detect earthquake and tremors by STA/LTA
                    btsta_1 = @elapsed S1 = Remove_eq.detect_eq_stalta(S1, float(stalta_longtimewindow), float(removal_shorttimewindow),
                                        float(stalta_threshold), float(overlap))
                end

                bt_3 = @elapsed S1 = Remove_eq.remove_eq(S1, S, float(invert_tukey_α), plot_kurtosis_α, max_wintaper_duration,
                                plot_boxheight, trunc(Int, plot_span), fodir, tstamp, tvec, IsSaveFig)

                bt_getkurtosis += bt_1
                bt_removeeq += bt_2 + bt_3 + btsta_1

                #if mod(dlid, round(0.1*NumofTimestamp)+1) == 0
                #    println([bt_1, bt_2, btsta_1, bt_3])
                #end

                #remove kurtosis for reduce size
                S1.misc["kurtosis"] = []
                S1.misc["eqtimewindow"] = []

            else
                #only STA/LTA
                if IsSTALTARemoval

                    bt_2 = @elapsed S1 = Remove_eq.detect_eq_stalta(S1, float(stalta_longtimewindow), float(removal_shorttimewindow),
                                        float(stalta_threshold), float(overlap))
                    bt_3 = @elapsed S1 = Remove_eq.remove_eq(S1, S, float(invert_tukey_α), plot_kurtosis_α, max_wintaper_duration,
                                    plot_boxheight, trunc(Int, plot_span), fodir, tstamp, tvec, IsSaveFig)

                    #remove kurtosis for reduce size
                    S1.misc["kurtosis"] = []
                    S1.misc["eqtimewindow"] = []

                    bt_getkurtosis += 0.0
                    bt_removeeq += bt_2 + bt_3

                else
                    #no removal process.
                    println("Both 'IsKurtosisRemoval' and 'IsKurtosisRemoval' are false. No removal process is executed.")
                    exit(0)
                end

            end

        else
            #download error found: save as it is.
            S1 = S
        end

        icount += 1
        SRall[icount] = S1

    end

    return (SRall, bt_getkurtosis, bt_removeeq)
end

end


text:
"""
---
Quadrature Amplitude Modulation (QAM) hard decoding function
Return the hard decoded constellation with voronoi baseds decision. The difference with bitDeMapping is that bitDeMapping returns the decoded bit sequence whereas hardConstellation returns the closest constellation point. This can be use to compute raw EVM estimation (assuming a sufficiently high SNR to avoid errors).
# --- Syntax
	hardConstellation!(qamDec,M,qamMat)
# --- Input parameters
- qamDec	: Vector to populate [Array{Complex{Float64},N}] with N = length(qamMat)
- M			: Constellation size (i.e 4 to 256)
- qamMat	: Vector to decode
# --- Output parameters
- []
"""
function hardConstellation!(qamThres,M, qamMat)
	# ----------------------------------------------------
	# --- Overall parameters
	# ----------------------------------------------------
	scalingFactor	= sqrt(2/3*(M-1));
	#qamThres = zeros(Complex{Float64},length(qamMat));
	# ----------------------------------------------------
	# --- Switch on modulation order
	# ----------------------------------------------------
	# Here all Voronoi region are defined with normalized constellation
	# Symbol remapper are recall from bitMappingQAM (without rescaling factor)
	if M == 4
		## --- QPSK demodulator
		# --- Gray encoding
		#  -1  1
		#   0  1
		# --- Voronoi region based on sign
		# < 0 --> -1
		# > 0 --> 1
		# --- Iterative hard demapping
		alphabet  = [-1 1];
		bounds	  = 0;
		for n = 1 : 1 : length(qamMat)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamMat[n]);
			# Getting area of interest: from area 1 to 4
			decI	= 1+ (realI > bounds[1]);
			# --- Decision for imag part
			# Getting only imag part
			imagP	= imag(qamMat[n]);
			# Getting area of interest: from area 1 to 4
			decQ	= 1+ (imagP > bounds[1]) ;
			# --- TO hard constellation
			qamThres[n] = alphabet[decI] + 1im*alphabet[decQ];
		end
	elseif M == 16
		## 16-QAM hard demapper
		# --- Gray encoding
		#  -3     -1     1     3
		#  10     11    01    00
		# --- Voronoi region definition
		bounds	  = [-2 0 2]/(sqrt(2/3*15));
		alphabet  = [-3 -1 1 3];
		#bounds	= [-sqrt(2)/2 0 sqrt(2)/2];
		# FIXME Voronoi bug in decision
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamMat)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamMat[n]);
			# Getting area of interest: from area 1 to 4
			decI	= 1+ (realI > bounds[1]) + (realI>bounds[2]) + (realI>bounds[3]);
			# --- Decision for imag part
			# Getting only imag part
			imagP	= imag(qamMat[n]);
			# Getting area of interest: from area 1 to 4
			decQ	= 1+ (imagP > bounds[1]) + (imagP>bounds[2]) + (imagP>bounds[3]);
			# --- TO hard constellation
			qamThres[n] = alphabet[decI] + 1im*alphabet[decQ];
		end
	elseif M == 64
		## 64-QAM hard demapper
		# --- Gray encoding
		#   -7     -5    -3    -1    +1    +3     +5     +7
		#  100    101   111   110   010   011    001    000
		# --- Voronoi region definition
		bounds = [-6 -4 -2 0 2 4 6]/(sqrt(2/3*63));
		alphabet  = [-7 -5 -3 -1 1 3 5 7];
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamMat)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamMat[n]);
			# Getting area of interest: from area 1 to 4
			decI	= 1+ (realI > bounds[1]) + (realI>bounds[2]) + (realI>bounds[3])+ (realI>bounds[4])+ (realI>bounds[5])+ (realI>bounds[6])+ (realI>bounds[7]);
			# --- Decision for imag part
			# Getting only imag part
			realQ	= imag(qamMat[n]);
			# Getting area of interest: from area 1 to 8
			decQ	= 1+ (realQ > bounds[1]) + (realQ>bounds[2]) + (realQ>bounds[3])+ (realQ>bounds[4])+ (realQ>bounds[5])+ (realQ>bounds[6])+(realQ>bounds[7]);
			# --- TO hard constellation
			qamThres[n] = alphabet[decI] + 1im*alphabet[decQ];
		end

	elseif M == 256
		## 64-QAM hard demapper
		# --- Gray encoding
		# ------------------------------------------------#
		#	   |       |         |        |        |        |        |        |
		#	 -15     -13 .     -11 .     -9 .     -7 .     -5 .     -3 .     -1
		#	0100    0101 .    0111 .   0110 .   0010 .   0011 .   0001 .   0000
		#	   |       |         |        |        |        |        |        |
		#     +1 .     +3 .     +5 .     +7 .     +9 .    +11 .    +13      +15.
		#   1000 .   1001 .   1011 .   1010 .   1110 .   1111 .   1101 .   1100
		# --- Voronoi region definition
		bounds = [-14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14]/(sqrt(2/3*255));
		alphabet  = collect(-15:2:15);
		# --- Iterative hard demapping
		for n = 1 : 1 : length(qamMat)
			# --- Decision for real part
			# Getting only real part
			realI	= real(qamMat[n]);
			# Getting area of interest: from area 1 to 15
			decI	= 1+ (realI > bounds[1]) + (realI>bounds[2]) + (realI>bounds[3])+ (realI>bounds[4])+ (realI>bounds[5])+ (realI>bounds[6])+ (realI>bounds[7])+ (realI>bounds[8])+(realI>bounds[9])+(realI>bounds[10])+(realI>bounds[11])+(realI>bounds[12])+(realI>bounds[13])+ (realI>bounds[14])+ (realI>bounds[15]) ;
			# --- Decision for imag part
			# Getting only imag part
			realQ	= imag(qamMat[n]);
			# Getting area of interest: from area 1 to 15
			decQ	= 1+ (realQ > bounds[1]) + (realQ>bounds[2]) + (realQ>bounds[3])+ (realQ>bounds[4])+ (realQ>bounds[5])+ (realQ>bounds[6])+ (realQ>bounds[7])+ (realQ>bounds[8])+(realQ>bounds[9])+(realQ>bounds[10])+(realQ>bounds[11])+(realQ>bounds[12])+(realQ>bounds[13])+ (realQ>bounds[14])+ (realQ>bounds[15]) ;
			# --- TO hard constellation
			qamThres[n] = alphabet[decI] + 1im*alphabet[decQ];
		end
	end
	qamThres .= qamThres ./ scalingFactor;
	# ---  Output on decoded sequence
	return qamThres;
end

"""
---
Quadrature Amplitude Modulation (QAM) hard decoding function
Return the hard decoded constellation with voronoi baseds decision. The difference with bitDeMapping is that bitDeMapping returns the decoded bit sequence whereas hardConstellation returns the closest constellation point. This can be use to compute raw EVM estimation (assuming a sufficiently high SNR to avoid errors).
# --- Syntax
qamDec = hardConstellation!(qamDec,M,qamMat)
# --- Input parameters
- M			: Constellation size (i.e 4 to 256)
- qamMat	: Vector to decode
# --- Output parameters
- qamDec	: Vector to populate [Array{Complex{Float64},N}] with N = length(qamMat)
"""
function hardConstellation(M,qamMat)
	qamThres  = zeros(Complex{Float64},length(qamMat));
	hardConstellation!(qamThres,M,qamMat);
	return qamThres;
end

# ----------------------------------------------------
# --- Multiple dispatch handling
# ----------------------------------------------------
# --- MD: String case for modulation order
function hardConstellation!(qamThres,M::String, qamMat)
	# --- Casting modulation order to int8
	if M == "QPSK" || M == "4-QAM" || M == "QAM-4" || M=="4QAM" || M == "QAM4"
		hardConstellation!(qamThres,4,qamMat);
	elseif M == "16-QAM" || M == "QAM-16" || M=="16QAM" || M == "QAM16"
		hardConstellation!(qamThres,16,qamMat);
	elseif M == "64-QAM" || M == "QAM-64" || M=="64QAM" || M == "QAM64"
		hardConstellation!(qamThres,64,qamMat);
	end
end
function hardConstellation!(M::String, qamMat)
	# --- Casting modulation order to int8
	if M == "QPSK" || M == "4-QAM" || M == "QAM-4" || M=="4QAM" || M == "QAM4"
		hardConstellation(4,qamMat);
	elseif M == "16-QAM" || M == "QAM-16" || M=="16QAM" || M == "QAM16"
		hardConstellation(16,qamMat);
	elseif M == "64-QAM" || M == "QAM-64" || M=="64QAM" || M == "QAM64"
		hardConstellation(64,qamMat);
	end

end


text:
addprocs(5)

using RoME, Distributions
using KernelDensityEstimate
using IncrementalInference
using TransformUtils, Colors

using RoMEPlotting, Gadfly
using KernelDensityEstimatePlotting
using JLD2

# for development
importall Caesar

include(joinpath(dirname(@__FILE__),"slamUtils.jl"))
Gadfly.push_theme(:default)

saswindow = 8
N = 100

beacon = :l1
# addFactor!(fg, [:l1], Prior(MvNormal(zeros(2), diagm([100;100]) )) )

symbolstart = 1;
factorgap = 5;

windowstart = 55;
nfactors = 40;
posegap = 2;

poses = Dict{Int,Array}();
nav = Dict{Int,Array}();
errorinds = [];
allframes = Dict{Int,Array}();
saveBF = Dict();

#datadir = "/media/data1/data/kayaks/07_18_parsed/";
savedirheader = "0813/0813";

fg = initfg()
addVariable!(fg, beacon, Point2)

for fitr in 1:nfactors
    navchecked = false;
    while ~navchecked
        sasframes = collect(windowstart:1:(windowstart+saswindow-1))
        navtemp = importdata_nav(sasframes);
        navchecked, errorind = sanitycheck_nav(navtemp);
        push!(errorinds,errorind);
        if navchecked
            nav[fitr] = navtemp;
            poses[fitr] = [Symbol("x$i") for i in symbolstart:(symbolstart+saswindow-1)]
            addsascluster_only_gps!(fg, poses[fitr], beacon, sasframes, N=N)
            allframes[fitr] = sasframes
        else
            windowstart+=errorind
        end
    end

    if !isInitialized(fg,:l1)
        IncrementalInference.doautoinit!(fg, :l1)
    end

    jldname2 = "$(pwd())/" * savedirheader * "_nav_$(fitr).jld"
    save(jldname2,"nav",nav,"poses",poses,"frames",allframes)

    tree = wipeBuildNewTree!(fg, drawpdf=false)
    inferOverTree!(fg, tree, N=N)

    PL = [];

    beacongt = [17.0499;1.7832];
    push!(PL,layer(x=[beacongt[1];],y=[beacongt[2];],label=String["Beacon Ground Truth";],Geom.point,Geom.label(hide_overlaps=false), order=2));

    l1fit = getKDEMean(getBelief(fg,:l1))
    push!(PL,layer(x=[l1fit[1];],y=[l1fit[2];],label=String["KDE Mean";],Geom.point,Geom.label(hide_overlaps=false),order=1));

    l1max = getKDEMax(getBelief(fg,:l1))
    push!(PL,layer(x=[l1max[1];],y=[l1max[2];],label=String["KDE Max";],Geom.point,Geom.label(hide_overlaps=false),order=1));

    push!(PL, plotKDE(getBelief(fg, :l1),levels=4, layers=true)... )
    L1 = getVal(fg, :l1)
    push!(PL, layer(x=L1[1,:],y=L1[2,:], Geom.histogram2d))

    for i in 1:fitr
        push!(PL,layer(x=nav[i][:,1],y=nav[i][:,2],Geom.point, Geom.path));
        saveBF[(fitr,i)] = approxConvFwdBFRaw(fg, poses[i], :l1, nav[i][1,:])
    end

    push!(PL,Coord.Cartesian(xmin=-100,xmax=150,ymin=-150,ymax=100))
    pl = plot(PL...);

    plname = savedirheader * "_$(fitr).png"
    Gadfly.draw(PNG(plname,1.5*20cm,1.5*20cm),pl)

    jldname = savedirheader * "_fg_$(fitr).jld"
    savejld(fg,file=jldname)

    symbolstart += (saswindow + factorgap)
    windowstart += (saswindow + posegap)
end

#poses1 = [Symbol("x$i") for i in 1:(1+saswindow-1)]
#windowstart = 100;
#sasframes = collect(windowstart:1:windowstart+saswindow-1)
#nav1 = importdata_nav(sasframes);
#addsascluster_only_gps!(fg, poses1, beacon, sasframes, N=N)

# writeGraphPdf(fg)
# @async run(`evince bt.pdf`)

#PLOTTING TOOLS

Gadfly.push_theme(:default)

PL = [];

beacongt = [17.0499;1.7832];
push!(PL,layer(x=[beacongt[1];],y=[beacongt[2];],label=String["Beacon Ground Truth";],Geom.point,Geom.label(hide_overlaps=false), order=2));

l1fit = fit(MvNormal,getVal(fg,:l1))
push!(PL,layer(x=[l1fit.μ[1];],y=[l1fit.μ[2];],label=String["Mean-Fit";],Geom.point,Geom.label(hide_overlaps=false),order=1));

l1max = getKDEMax(getBelief(fg,:l1))
push!(PL,layer(x=[l1max[1];],y=[l1max[2];],label=String["Max-Fit";],Geom.point,Geom.label(hide_overlaps=false),order=1));

push!(PL, plotKDE(getBelief(fg, :l1),levels=4, layers=true)... )
#L1 = getVal(fg, :l1)
#push!(PL, layer(x=L1[1,:],y=L1[2,:], Geom.histogram2d))

saveBF = Dict()
for i in 1:nfactors
    #push!(PL,layer(x=nav[i][:,1],y=nav[i][:,2],Geom.point, Geom.path, Theme(default_color=colorant"red")));
    #push!(PL, plotKDE(getBelief(fg, Symbol("x$(i)")),levels=4, layers=true, dimLbls=["",""])... )
    saveBF[i] = approxConvFwdBFRaw(fg, poses[i], :l1, nav[i][1,:])
end

jldnameBF = "$(pwd())/0809/0809_BF_all.jld"
save(jldnameBF,"saveBF",saveBF)

push!(PL,Coord.Cartesian(xmin=150,xmax=250,ymin=-200,ymax=-50))

pl = plot(PL...)

frames = collect(15:425);
dataOut = zeros(length(frames),2);
count = 0
for i in frames
    count+=1
    navfile = datadir*"/nav$(i).csv"
    dataOut[count,:] = readdlm(navfile,',',Float64,'\n')
end

ll = ls(fg,:l1)

stuff = IncrementalInference.localProduct(fg,:l1)

allv = collect(1:12)
lo =[4,7,9,10,11]
# keep = [1;2;3;5;6;8;12]
keep = setdiff(allv, lo)
plotKDE([getBelief(fg,:l1);stuff[2][keep]],levels=1,c=["red";["blue" for i in 1:12]])

something, = predictbelief(fg, :l1, ll[keep])

plotKDE([getBelief(fg,:l1);kde!(something);stuff[2][keep]],levels=1,c=["red";"cyan";["blue" for i in 1:12]])

plotKDE([getBelief(fg,:l1);kde!(something);stuff[2][keep]],levels=1,c=["red";"cyan";["blue" for i in 1:12]], dims=[1])

navfile = datadir*"inav.csv"
inav = readdlm(navfile,',',Float64,'\n')

push!(PL,layer(x=dataOut[:,1],y=dataOut[:,2],Geom.point, Geom.path,Theme(default_color=colorant"green")),
     layer(x=inav[:,1],y=inav[:,2]))

Gadfly.draw(SVG("stored/07-27-530-2.svg",1.5*20cm,1.5*20cm),pl)
savejld(fg,file="stored/07-27-530.jld")

## PREVIOUS DEVELOPMENT AND TESTING CODE BELOW

# backwards solve

# backwards leave one out case
setVal!(fg, :l1, 0.01*randn(2,N))
reset!(sas2d.dbg)
predX2 = IncrementalInference.approxConv(fg, :l1x1x2x3x4x5x6x7x8x9f1, :x9, N=N)

pl = plot(
layer(x=predX2[1,:], y=predX2[2,:], Geom.histogram2d(ybincount=200)),
layer(x=posData[1:9,1],y=posData[1:9,2], Geom.point)
)

pl.coord = Coord.Cartesian(xmin=0,xmax=20,ymin=-55,ymax=-35)
pl

## VISUALIZATION AND BEAMS DEBUGGIN

pl = plotKDE(kde!(predL1))
# pl.coord = Gadfly.Coord.Cartesian(xmin=-200, xmax=100, ymin=-100, ymax=200);
pl

# Look at sas2d.dbg
# @show length(sas2d.dbg)
plot(x=sas2d.dbg.azi_smpls, Geom.histogram)

n=200; pl = plot(x=linspace(0,2pi,length(sas2d.dbg.beams[n])), y=sas2d.dbg.beams[n], Geom.path())

n=200; pl = plot(x=linspace(0,2pi,length(sas2d.dbg.beams_cut[n])), y=sas2d.dbg.beams_cut[n], Geom.path())

pl = plot(x=linspace(0,2pi,length(avg_beam)), y=avg_beam, Geom.path(), Coord.Cartesian(ymin=0))
pl = plot(x=linspace(0,2pi,length(avg_beam_cut)), y=avg_beam_cut, Geom.path(), Coord.Cartesian(ymin=0))

#


text:
##############################################################################
##
## Iterate on terms
##
##############################################################################

eachterm(@nospecialize(x::AbstractTerm)) = (x,)
eachterm(@nospecialize(x::NTuple{N, AbstractTerm})) where {N} = x
TermOrTerms = Union{AbstractTerm, NTuple{N, AbstractTerm} where N}
hasintercept(@nospecialize(t::TermOrTerms)) =
    InterceptTerm{true}() ∈ terms(t) ||
    ConstantTerm(1) ∈ terms(t)
omitsintercept(@nospecialize(f::FormulaTerm)) = omitsintercept(f.rhs)
omitsintercept(@nospecialize(t::TermOrTerms)) =
    InterceptTerm{false}() ∈ terms(t) ||
    ConstantTerm(0) ∈ terms(t) ||
    ConstantTerm(-1) ∈ terms(t)
##############################################################################
##
## Parse IV
##
##############################################################################

function parse_iv(@nospecialize(f::FormulaTerm))
	for term in eachterm(f.rhs)
		if term isa FormulaTerm
            both = intersect(eachterm(term.lhs), eachterm(term.rhs))
            endos = setdiff(eachterm(term.lhs), both)
            exos = setdiff(eachterm(term.rhs), both)
            # otherwise empty collection. Later, there will be check
            isempty(endos) && throw("There are no endogeneous variables")
            length(exos) < length(endos) && throw("Model not identified. There must be at least as many instrumental variables as endogeneneous variables")
			formula_endo = FormulaTerm(ConstantTerm(0), tuple(ConstantTerm(0), endos...))
			formula_iv = FormulaTerm(ConstantTerm(0), tuple(ConstantTerm(0), exos...))
            formula_exo = FormulaTerm(f.lhs, tuple((term for term in eachterm(f.rhs) if !isa(term, FormulaTerm))..., both...))
            return formula_exo, formula_endo, formula_iv
		end
	end
	return f, nothing, nothing
end

##############################################################################
##
## Parse FixedEffect
##
##############################################################################
struct FixedEffectTerm <: AbstractTerm
    x::Symbol
end
StatsModels.termvars(t::FixedEffectTerm) = [t.x]
fe(x::Term) = FixedEffectTerm(Symbol(x))
fe(s::Symbol) = FixedEffectTerm(s)

has_fe(::FixedEffectTerm) = true
has_fe(::FunctionTerm{typeof(fe)}) = true
has_fe(t::InteractionTerm) = any(has_fe(x) for x in t.terms)
has_fe(::AbstractTerm) = false
has_fe(@nospecialize(t::FormulaTerm)) = any(has_fe(x) for x in eachterm(t.rhs))

fesymbol(t::FixedEffectTerm) = t.x
fesymbol(t::FunctionTerm{typeof(fe)}) = Symbol(t.args_parsed[1])

"""
    parse_fixedeffect(data, formula::FormulaTerm)
    parse_fixedeffect(data, ts::NTuple{N, AbstractTerm})

Construct any `FixedEffect` specified with a `FixedEffectTerm`.

# Returns
- `Vector{FixedEffect}`: a collection of all `FixedEffect`s constructed.
- `Vector{Symbol}`: names assigned to the fixed effect estimates (can be used as column names).
- `Vector{Symbol}`: names of original fixed effects.
- `FormulaTerm` or `NTuple{N, AbstractTerm}`: `formula` or `ts` without any term related to fixed effects (an intercept may be explicitly omitted if necessary).
"""
function parse_fixedeffect(data, @nospecialize(formula::FormulaTerm))
    fes = FixedEffect[]
    ids = Symbol[]
    fekeys = Symbol[]
     for term in eachterm(formula.rhs)
        result = _parse_fixedeffect(data, term)
        if result !== nothing
            push!(fes, result[1])
            push!(ids, result[2])
            append!(fekeys, result[3])
        end
    end
    if !isempty(fes)
        if any(fe.interaction isa UnitWeights for fe in fes)
            formula = FormulaTerm(formula.lhs, (InterceptTerm{false}(), (term for term in eachterm(formula.rhs) if !isa(term, Union{ConstantTerm,InterceptTerm}) && !has_fe(term))...))
        else
            formula = FormulaTerm(formula.lhs, Tuple(term for term in eachterm(formula.rhs) if !has_fe(term)))
        end
    end
    return fes, ids, unique(fekeys), formula
end

# Method for external packages
function parse_fixedeffect(data, @nospecialize(ts::NTuple{N, AbstractTerm})) where N
    fes = FixedEffect[]
    ids = Symbol[]
    fekeys = Symbol[]
    for term in eachterm(ts)
        result = _parse_fixedeffect(data, term)
        if result !== nothing
            push!(fes, result[1])
            push!(ids, result[2])
            append!(fekeys, result[3])
        end
    end
    if !isempty(fes)
        if any(fe.interaction isa UnitWeights for fe in fes)
            ts = (InterceptTerm{false}(), (term for term in eachterm(ts) if !isa(term, Union{ConstantTerm,InterceptTerm}) && !has_fe(term))...)
        else
            ts = Tuple(term for term in eachterm(ts) if !has_fe(term))
        end
    end
    return fes, ids, unique(fekeys), ts
end

# Construct FixedEffect from a generic term
function _parse_fixedeffect(data, @nospecialize(t::AbstractTerm))
    if has_fe(t)
        st = fesymbol(t)
        return FixedEffect(Tables.getcolumn(data, st)), Symbol(:fe_, st), [st]
    end
end

# Construct FixedEffect from an InteractionTerm
function _parse_fixedeffect(data, @nospecialize(t::InteractionTerm))
    fes = (x for x in t.terms if has_fe(x))
    interactions = (x for x in t.terms if !has_fe(x))
    if !isempty(fes)
        # x1&x2 from (x1&x2)*id
        fe_names = [fesymbol(x) for x in fes]
        v1 = _multiply(data, Symbol.(interactions))
        fe = FixedEffect((Tables.getcolumn(data, fe_name) for fe_name in fe_names)...; interaction = v1)
        interactions = string.(interactions)
        s = vcat(["fe_" * string(fe_name) for fe_name in fe_names], interactions)
        return fe, Symbol(reduce((x1, x2) -> x1*"&"*x2, s)), fe_names
    end
end

function _multiply(data, ss::AbstractVector)
    if isempty(ss)
        return uweights(size(data, 1))
    elseif length(ss) == 1
        # in case it has missing (for some reason *(missing) not defined) iv VERSION < 1.6
        # do NOT use ! since it would modify the vector
        return convert(AbstractVector{Float64}, replace(Tables.getcolumn(data, ss[1]), missing => 0))
    else
        return convert(AbstractVector{Float64}, replace!(.*((Tables.getcolumn(data, x) for x in ss)...), missing => 0))
    end
end


text:
using SimpleFanControl
using Test

@testset "SimpleFanControl.jl" begin
    # Write your tests here.
end


text:
"""
Returns (N1, N2) with `N1 x N2 == n`. N2 might become 1
"""
function close2square(n::Real)
    # a cannot be greater than the square root of n
    # b cannot be smaller than the square root of n
    # we get the maximum allowed value of a
    amax = floor(Int, sqrt(n));
    if 0 == rem(n, amax)
        # special case where n is a square number
        return (amax, div(n, amax))
    end
    # Get its prime factors of n
    primeFactors  = factor(n);
    # Start with a factor 1 in the list of candidates for a
    candidates = [1]
    for (f, _) in primeFactors
        # Add new candidates which are obtained by multiplying
        # existing candidates with the new prime factor f
        # Set union ensures that duplicate candidates are removed
        candidates = union(candidates, f .* candidates)
        # throw out candidates which are larger than amax
        filter!(x-> x <= amax, candidates)
    end
    # Take the largest factor in the list d
    (candidates[end], div(n, candidates[end]))
end

same_length_array(array, value::NativeFont) = Iterators.repeated(value, length(array))
function extrema_nan(x::ClosedInterval)
    (minimum(x), maximum(x))
end

# TODO upgrade GeometryTypes to do these kind of things.
# Problem: I don't really want to introduce a depedency on intervals
function Base.in(point::StaticVector{N}, rectangle::HyperRectangle{N}) where N
    mini, maxi = minimum(rectangle), maximum(rectangle)
    for i = 1:N
        point[i] in (mini[i] .. maxi[i]) || return false
    end
    return true
end

to_range(x::ClosedInterval) = (minimum(x), maximum(x))
to_range(x::VecTypes{2}) = x
to_range(x::AbstractRange) = (minimum(x), maximum(x))
to_range(x::AbstractVector) = (minimum(x), maximum(x))
function to_range(x::AbstractArray)
    if length(x) in size(x) # assert that just one dim != 1
        to_range(vec(x))
    else
        error("Can't convert to a range. Please supply a range/vector/interval or a tuple (min, max)")
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 03c25682-4e46-11eb-1b97-33699faa080f
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Plots
	using Flux
	using Statistics
end

# ╔═╡ 1290b1fa-4e4d-11eb-111b-37ade59c8f30
md"""
# Chapter06 Temporal-Difference Learning (Cliff Walking)
"""

# ╔═╡ 0f4a2ba0-4e4d-11eb-0134-5d4fb9897b33
md"""
In Example 6.6, a gridworld example of **Cliff Walking** is introduced to compare the difference between on-policy (SARSA) and off-policy (Q-learning). Although there's a package of [GridWorlds.jl](https://github.com/JuliaReinforcementLearning/GridWorlds.jl) dedicated to 2-D environments, we decide to write an independent implementation here as a showcase.
"""

# ╔═╡ 63a013f0-4e4e-11eb-2737-ff5104b993f2
begin
	const NX = 4
	const NY = 12
	const Start = CartesianIndex(4, 1)
	const Goal = CartesianIndex(4, 12)
	const LRUD = [
		CartesianIndex(0, -1),  # left
		CartesianIndex(0, 1),   # right
		CartesianIndex(-1, 0),  # up
		CartesianIndex(1, 0),   # down
	]
	const LinearInds = LinearIndices((NX, NY))
end

# ╔═╡ 89c9010e-4e4e-11eb-2d97-05ee771fd5f4
function iscliff(p::CartesianIndex{2})
    x, y = Tuple(p)
    x == 4 && y > 1 && y < NY
end

# ╔═╡ eb5e16fc-4e4e-11eb-0ce7-5b5510c7924a
# take a look at the wordmap
heatmap((!iscliff).(CartesianIndices((NX, NY))); yflip = true)

# ╔═╡ 4e8e67ac-4e51-11eb-0742-ad46e685a921
begin
	Base.@kwdef mutable struct CliffWalkingEnv <: AbstractEnv
		position::CartesianIndex{2} = Start
	end
	function (env::CliffWalkingEnv)(a::Int)
		x, y = Tuple(env.position + LRUD[a])
		env.position = CartesianIndex(min(max(x, 1), NX), min(max(y, 1), NY))
	end
end

# ╔═╡ af4da300-4e51-11eb-3821-0522243eea6a
RLBase.state(env::CliffWalkingEnv) = LinearInds[env.position]

# ╔═╡ 94e72e6e-4e51-11eb-25e4-81c253969de6
RLBase.state_space(env::CliffWalkingEnv) = Base.OneTo(length(LinearInds))

# ╔═╡ be36d904-4e51-11eb-1ead-1d89513b68e5
RLBase.action_space(env::CliffWalkingEnv) = Base.OneTo(length(LRUD))

# ╔═╡ e087f574-4e51-11eb-1abf-795bdc85f7b1
RLBase.reward(env::CliffWalkingEnv) = env.position == Goal ? 0.0 : (iscliff(env.position) ? -100.0 : -1.0)

# ╔═╡ f97f149a-4e51-11eb-22d9-f9ec9ed837a8
RLBase.is_terminated(env::CliffWalkingEnv) = env.position == Goal || iscliff(env.position)

# ╔═╡ 8e26fbea-4e51-11eb-3ae8-312009b7c7c1
RLBase.reset!(env::CliffWalkingEnv) = env.position = Start

# ╔═╡ 3d8ee702-4e52-11eb-031a-bf64267362d3
world = CliffWalkingEnv()

# ╔═╡ 30b3604c-4e52-11eb-277f-d5c354e9ec66
md"""
Now we have a workable environment. Next we create several factories to generate different policies for comparison.
"""

# ╔═╡ fd5eb880-4e52-11eb-0c30-618e1e9f8f17
begin
	NS = length(state_space(world))
	NA = length(action_space(world))
end

# ╔═╡ 5da833fc-4e52-11eb-3f96-1d129eb12a00
create_agent(α, method) = Agent(
	policy = QBasedPolicy(
		learner=TDLearner(
			approximator=TabularQApproximator(
				;n_state=NS,
				n_action=NA,
				opt=Descent(α),
			),
			method=method,
			γ=1.0,
			n=0
		),
		explorer=EpsilonGreedyExplorer(0.1)
	),
	trajectory=VectorSARTTrajectory()
)

# ╔═╡ 95bf52b0-4e53-11eb-1740-2fff3b8a9fbe
function repeated_run(α, method, N, n_episode, is_mean=true)
	env = CliffWalkingEnv()
	rewards = []
	for _ in 1:N
		h = TotalRewardPerEpisode(;is_display_on_exit=false)
		run(
			create_agent(α, method),
			env,
			StopAfterEpisode(n_episode;is_show_progress=false),
			h
		)
		push!(rewards, is_mean ? mean(h.rewards) : h.rewards)
	end
	mean(rewards)
end

# ╔═╡ 33ab57f2-4e55-11eb-3b53-89492ca4b5f6
begin
	p = plot(legend=:bottomright, xlabel="Episodes", ylabel="Sum of rewards during episode")
	plot!(p, repeated_run(0.5, :SARS, 1000, 500, false), label="QLearning")
	plot!(p, repeated_run(0.5, :SARSA, 1000, 500, false), label="SARSA")
	p
end

# ╔═╡ a9f81264-4e5b-11eb-2634-cb2179ea9e6b
begin
	A = 0.1:0.05:0.95
	fig_6_3 = plot(;legend=:bottomright, xlabel="α", ylabel="Sum of rewards per episode")

	plot!(fig_6_3, A, [repeated_run(α, :SARS, 100, 100) for α in A], linestyle=:dash ,markershape=:rect, label="Interim Q")
	plot!(fig_6_3, A, [repeated_run(α, :SARSA, 100, 100) for α in A], linestyle=:dash, markershape=:dtriangle, label="Interim SARSA")
	plot!(fig_6_3, A, [repeated_run(α, :ExpectedSARSA, 100, 100) for α in A], linestyle=:dash, markershape=:cross, label="Interim ExpectedSARSA")

	plot!(fig_6_3, A, [repeated_run(α, :SARS, 10, 5000) for α in A], linestyle=:solid ,markershape=:rect, label="Asymptotic interim Q")
	plot!(fig_6_3, A, [repeated_run(α, :SARSA, 10, 5000) for α in A], linestyle=:solid, markershape=:dtriangle, label="Asymptotic SARSA")
	plot!(fig_6_3, A, [repeated_run(α, :ExpectedSARSA, 10, 5000) for α in A], linestyle=:solid, markershape=:cross, label="Asymptotic ExpectedSARSA")
	fig_6_3
end

# ╔═╡ Cell order:
# ╟─1290b1fa-4e4d-11eb-111b-37ade59c8f30
# ╠═03c25682-4e46-11eb-1b97-33699faa080f
# ╟─0f4a2ba0-4e4d-11eb-0134-5d4fb9897b33
# ╠═63a013f0-4e4e-11eb-2737-ff5104b993f2
# ╠═89c9010e-4e4e-11eb-2d97-05ee771fd5f4
# ╠═eb5e16fc-4e4e-11eb-0ce7-5b5510c7924a
# ╠═4e8e67ac-4e51-11eb-0742-ad46e685a921
# ╠═af4da300-4e51-11eb-3821-0522243eea6a
# ╠═94e72e6e-4e51-11eb-25e4-81c253969de6
# ╠═be36d904-4e51-11eb-1ead-1d89513b68e5
# ╠═e087f574-4e51-11eb-1abf-795bdc85f7b1
# ╠═f97f149a-4e51-11eb-22d9-f9ec9ed837a8
# ╠═8e26fbea-4e51-11eb-3ae8-312009b7c7c1
# ╠═3d8ee702-4e52-11eb-031a-bf64267362d3
# ╟─30b3604c-4e52-11eb-277f-d5c354e9ec66
# ╠═fd5eb880-4e52-11eb-0c30-618e1e9f8f17
# ╠═5da833fc-4e52-11eb-3f96-1d129eb12a00
# ╠═95bf52b0-4e53-11eb-1740-2fff3b8a9fbe
# ╠═33ab57f2-4e55-11eb-3b53-89492ca4b5f6
# ╠═a9f81264-4e5b-11eb-2634-cb2179ea9e6b


text:
using FDM: grad, jacobian, _jvp, _j′vp, jvp, j′vp, to_vec

# Dummy type where length(x::DummyType) ≠ length(first(to_vec(x)))
struct DummyType{TX<:Matrix}
    X::TX
end

function FDM.to_vec(x::DummyType)
    x_vec, back = to_vec(x.X)
    return x_vec, x_vec -> DummyType(back(x_vec))
end

Base.:(==)(x::DummyType, y::DummyType) = x.X == y.X
Base.length(x::DummyType) = size(x.X, 1)

@testset "grad" begin

    @testset "grad" begin
        rng, fdm = MersenneTwister(123456), central_fdm(5, 1)
        x = randn(rng, 2)
        xc = copy(x)
        @test grad(fdm, x->sin(x[1]) + cos(x[2]), x) ≈ [cos(x[1]), -sin(x[2])]
        @test xc == x
    end

    function check_jac_and_jvp_and_j′vp(fdm, f, ȳ, x, ẋ, J_exact)
        xc = copy(x)
        @test jacobian(fdm, f, x, length(ȳ)) ≈ J_exact
        @test jacobian(fdm, f, x) == jacobian(fdm, f, x, length(ȳ))
        @test _jvp(fdm, f, x, ẋ) ≈ J_exact * ẋ
        @test _j′vp(fdm, f, ȳ, x) ≈ J_exact' * ȳ
        @test xc == x
    end

    @testset "jacobian / _jvp / _j′vp" begin
        rng, P, Q, fdm = MersenneTwister(123456), 3, 2, central_fdm(5, 1)
        ȳ, A, x, ẋ = randn(rng, P), randn(rng, P, Q), randn(rng, Q), randn(rng, Q)
        Ac = copy(A)

        check_jac_and_jvp_and_j′vp(fdm, x->A * x, ȳ, x, ẋ, A)
        @test Ac == A
        check_jac_and_jvp_and_j′vp(fdm, x->sin.(A * x), ȳ, x, ẋ, cos.(A * x) .* A)
        @test Ac == A
    end

    function test_to_vec(x)
        x_vec, back = to_vec(x)
        @test x_vec isa Vector
        @test x == back(x_vec)
        return nothing
    end

    @testset "to_vec" begin
        test_to_vec(1.0)
        test_to_vec(1)
        test_to_vec(randn(3))
        test_to_vec(randn(5, 11))
        test_to_vec(randn(13, 17, 19))
        test_to_vec(randn(13, 0, 19))
        test_to_vec(UpperTriangular(randn(13, 13)))
        test_to_vec(Symmetric(randn(11, 11)))
        test_to_vec(Diagonal(randn(7)))
        test_to_vec(DummyType(randn(2, 9)))

        @testset "$T" for T in (Adjoint, Transpose)
            test_to_vec(T(randn(4, 4)))
            test_to_vec(T(randn(6)))
            test_to_vec(T(randn(2, 5)))
        end

        @testset "Tuples" begin
            test_to_vec((5, 4))
            test_to_vec((5, randn(5)))
            test_to_vec((randn(4), randn(4, 3, 2), 1))
            test_to_vec((5, randn(4, 3, 2), UpperTriangular(randn(4, 4)), 2.5))
            test_to_vec(((6, 5), 3, randn(3, 2, 0, 1)))
            test_to_vec((DummyType(randn(2, 7)), DummyType(randn(3, 9))))
            test_to_vec((DummyType(randn(3, 2)), randn(11, 8)))
        end
    end

    @testset "jvp" begin
        rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
        x, y = randn(rng, N), randn(rng, M)
        ẋ, ẏ = randn(rng, N), randn(rng, M)
        xy, ẋẏ = vcat(x, y), vcat(ẋ, ẏ)
        ż_manual = _jvp(fdm, (xy)->sum(sin, xy), xy, ẋẏ)[1]
        ż_auto = jvp(fdm, x->sum(sin, x[1]) + sum(sin, x[2]), ((x, y), (ẋ, ẏ)))
        ż_multi = jvp(fdm, (x, y)->sum(sin, x) + sum(sin, y), (x, ẋ), (y, ẏ))
        @test ż_manual ≈ ż_auto
        @test ż_manual ≈ ż_multi
    end

    @testset "j′vp" begin
        rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
        x, y = randn(rng, N), randn(rng, M)
        z̄ = randn(rng, N + M)
        xy = vcat(x, y)
        x̄ȳ_manual = j′vp(fdm, xy->sin.(xy), z̄, xy)
        x̄ȳ_auto = j′vp(fdm, x->sin.(vcat(x[1], x[2])), z̄, (x, y))
        x̄ȳ_multi = j′vp(fdm, (x, y)->sin.(vcat(x, y)), z̄, x, y)
        @test x̄ȳ_manual ≈ vcat(x̄ȳ_auto...)
        @test x̄ȳ_manual ≈ vcat(x̄ȳ_multi...)
    end
end


text:
#!/usr/bin/env julia

# Tests against astropy.

using AstroAngles
using DelimitedFiles
using Printf
using SkyCoords
using Statistics
using Test

import SkyCoords: lat, lon

datapath = joinpath(dirname(@__FILE__), "data")

rad2arcsec(r) = 3600 * rad2deg(r)

# input coordinates
fname = joinpath(datapath, "input_coords.csv")
indata, inhdr = readdlm(fname, ','; header = true)

# Float32 has a large tolerance compared to Float64 and BigFloat, but here we
# are more interested in making sure that the infrastructure works for different
# floating types.
@testset "Testing $F" for (F, TOL) in (
    (Float32, 0.2),
    (Float64, 0.0001),
    (BigFloat, 0.0001),
)
    for (insys, T) in (
        ("icrs", ICRSCoords{F}),
        ("fk5j2000", FK5Coords{2000,F}),
        ("fk5j1975", FK5Coords{1975,F}),
        ("gal", GalCoords{F}),
    )
        c_in = T[T(indata[i, 1], indata[i, 2]) for i = 1:size(indata, 1)]
        for c in c_in
            @test convert(T, c) == c
        end

        for (outsys, S) in (
            ("icrs", ICRSCoords{F}),
            ("fk5j2000", FK5Coords{2000,F}),
            ("fk5j1975", FK5Coords{1975,F}),
            ("gal", GalCoords{F}),
        )
            (outsys == insys) && continue

            c_out = S[convert(S, c) for c in c_in]

            # Test pipe and constructor conversion
            @test c_out == S[S(c) for c in c_in]
            @test c_out == S[c |> S for c in c_in]

            # Read in reference answers.
            ref_fname = joinpath(datapath, "$(insys)_to_$(outsys).csv")
            refdata, hdr = readdlm(ref_fname, ','; header = true)
            c_ref = S[S(refdata[i, 1], refdata[i, 2]) for i = 1:size(refdata, 1)]

            # compare
            sep = separation.(c_out, c_ref)
            maxsep = rad2arcsec(maximum(sep))
            meansep = rad2arcsec(mean(sep))
            minsep = rad2arcsec(minimum(sep))
            @printf "%8s --> %8s : max=%6.4f\"  mean=%6.4f\"   min=%6.4f\"\n" insys outsys maxsep meansep minsep
            @test maxsep < TOL
        end
    end
end

# Test separation between coordinates and conversion with mixed floating types.
@testset "Separation" begin
    c1 = ICRSCoords(ℯ, pi / 2)
    c5 = ICRSCoords(ℯ, 1 + pi / 2)
    @test separation(c1, c5) ≈ separation(c5, c1) ≈ separation(c1, convert(GalCoords, c5)) ≈
          separation(convert(FK5Coords{1980}, c5), c1) ≈ 1
    for T in (GalCoords, FK5Coords{2000})
        c2 = convert(T{Float32}, c1)
        c3 = convert(T{Float64}, c1)
        c4 = convert(T{BigFloat}, c1)
        @test typeof(c2) === T{Float32}
        @test typeof(c3) === T{Float64}
        @test typeof(c4) === T{BigFloat}
        @test isapprox(lat(c2), lat(c3), rtol = sqrt(eps(Float32)))
        @test isapprox(lat(c3), lat(c4), rtol = sqrt(eps(Float64)))
        @test isapprox(lon(c2), lon(c3), rtol = sqrt(eps(Float32)))
        @test isapprox(lon(c3), lon(c4), rtol = sqrt(eps(Float64)))
        c6 = convert(T, c5)
        @test separation(c3, c6) ≈ separation(c6, c3) ≈ 1
    end
end

@testset "string construction" for C in [
    ICRSCoords,
    GalCoords,
    FK5Coords{2000},
    FK5Coords{1970},
]
    @test C(hms"0h0m0", dms"0d0m0") == C(0.0, 0.0)
    @test C(hms"12h0.0m0.0s", dms"90:0:0") == C(π, π / 2)
    @test C(hms"18h0:0", dms"90:0:0") == C(3π / 2, π / 2)
    @test C(hms"12:0:0", dms"90:0:0") == C(π, π / 2)
end

# Test separation between coordinates and conversion with mixed floating types.
@testset "Position Angles" begin
    c1 = ICRSCoords(0, 0)
    c2 = ICRSCoords(deg2rad(1), 0)

    # interface
    @test @inferred position_angle(c1, c2) ≈ @inferred position_angle(c1, c2 |> GalCoords)
    @test position_angle(c1, c2) ≈ position_angle(c1, c2 |> GalCoords)

    # accuracy
    @test position_angle(c1, c2) ≈ π / 2

    c3 = ICRSCoords(deg2rad(1), deg2rad(0.1))
    @test position_angle(c1, c3) < π / 2

    c4 = ICRSCoords(0, deg2rad(1))
    @test position_angle(c1, c4) ≈ 0

    # types
    for T in [ICRSCoords, GalCoords, FK5Coords{2000}]
        c1 = T(0, 0)
        c2 = T(deg2rad(1), 0)
        @test position_angle(c1, c2) ≈ π / 2
    end
end

@testset "Offset ($T1, $T2)" for T1 in [ICRSCoords, GalCoords, FK5Coords{2000}], T2 in [ICRSCoords, GalCoords, FK5Coords{2000}]
    # simple integration tests, depend that separation and position_angle are accurate
    c1s = [
        T1(0, -π/2), # south pole
        T1(0, π/2), # north pole
        T1(deg2rad(1), deg2rad(2))
    ]
    c2 = T2(deg2rad(5), deg2rad(10))

    for c1 in c1s
        sep, pa = @inferred offset(c1, c2)
        test_c2 = @inferred offset(c1, sep, pa)
        @test test_c2 isa T1
        test_c2 = T2(test_c2)
        @test lon(test_c2) ≈ lon(c2)
        @test lat(test_c2) ≈ lat(c2)
    end

    # specific cases to cover special cases.
    c1 = T1(0, deg2rad(89))
    for (pa, sep) in [(0, 2), (180, 358)]
        sep = deg2rad(sep)
        pa = deg2rad(pa)
        c2 = offset(c1, sep, pa)
        @test lon(c2) |> rad2deg ≈ 180
        @test lat(c2) |> rad2deg ≈ 89

        c2 = offset(c1, 2sep, pa)
        @test lon(c2) |> rad2deg ≈ 180
        @test lat(c2) |> rad2deg ≈ 87
    end

    # verify antipode
    c1 = T1(deg2rad(10), deg2rad(47))
    for pa in range(0, stop=377, length=10)
        c2 = offset(c1, deg2rad(180), deg2rad(pa))
        @test lon(c2) |> rad2deg ≈ 190
        @test lat(c2) |> rad2deg ≈ -47

        c2 = offset(c1, deg2rad(360), deg2rad(pa))
        @test lon(c2) |> rad2deg ≈ 10
        @test lat(c2) |> rad2deg ≈ 47
    end

    c1 = T1(deg2rad(10), deg2rad(60))
    c2 = offset(c1, deg2rad(1), deg2rad(90))
    @test 11.9 < lon(c2) |> rad2deg < 12.0
    @test 59.9 < lat(c2) |> rad2deg < 60.0
end