text:
module pplot_

export
	pplot

using IUP
using IUP_IM
using IUP_CD

type Handles
	figure1::Ptr{Ihandle}
	iup_canvas::Ptr{Ihandle}
	cd_canvas::Ptr{cdCanvas}	# cdCanvas is a composite type
end

global const MAXPLOT = 6
global plot = Array(Ptr{Ihandle}, MAXPLOT)

#=
global dial1::Ptr{Ihandle}, dial2::Ptr{Ihandle},          # dials for zooming */
        tgg1::Ptr{Ihandle}, tgg2::Ptr{Ihandle},            # auto scale on|off toggles */
        tgg3::Ptr{Ihandle}, tgg4::Ptr{Ihandle},            # grid show|hide toggles */
        tgg5::Ptr{Ihandle},                   		# legend show|hide toggle */
        tabs::Ptr{Ihandle}
=#
global dial1, dial2
global tgg1, tgg2
global tgg3, tgg4
global tgg5
global tabs

include("../src/handlegraphics.jl")

# ----------------------------------------------------------------------------------
function pplot()

global dial1, dial2
global tgg1, tgg2
global tgg3, tgg4
global tgg5
global tabs

	IupOpen()	#Initializes IUP

	vboxr = Array(Ptr{Ihandle}, MAXPLOT+1) 	# tabs containing the plots

	IupControlsOpen();  # init the addicional controls library (we use IupTabs)
	IupPlotOpen();     # init IupPlot library

	#  cdInitGdiPlus();

	# create plots
	for (ii = 1:MAXPLOT)
		plot[ii] = IupPlot()
	end

	# left panel: plot control
	# Y zooming
	dial1 = IupDial("VERTICAL");
	lbl1 = IupLabel("+");
	lbl2 = IupLabel("-");
	boxinfo = IupVbox(lbl1, IupFill(), lbl2, C_NULL);
	boxdial1 = IupHbox(boxinfo, dial1, C_NULL);

	IupSetAttribute(boxdial1, "ALIGNMENT", "ACENTER");
	IupSetAttribute(boxinfo, "ALIGNMENT", "ACENTER");
	IupSetAttribute(boxinfo, "SIZE", "20x52");
	IupSetAttribute(boxinfo, "GAP", "2");
	IupSetAttribute(boxinfo, "MARGIN", "4");
	IupSetAttribute(boxinfo, "EXPAND", "YES");
	IupSetAttribute(lbl1, "EXPAND", "NO");
	IupSetAttribute(lbl2, "EXPAND", "NO");

	IupSetAttribute(dial1, "ACTIVE", "NO");
	IupSetAttribute(dial1, "SIZE", "20x52");
	IupSetCallback(dial1, "BUTTON_PRESS_CB", cfunction(dial1_btndown_cb, Int, (Ptr{Ihandle},)))
	IupSetCallback(dial1, "MOUSEMOVE_CB", cfunction(dial1_btnup_cb, Int, (Ptr{Ihandle}, Cdouble)) )
	IupSetCallback(dial1, "BUTTON_RELEASE_CB", cfunction(dial1_btnup_cb, Int, (Ptr{Ihandle}, Cdouble)) )

	tgg1 = IupToggle("Y Autoscale");
	IupSetCallback(tgg1, "ACTION", cfunction(tgg1_cb, Int, (Ptr{Ihandle}, Cint)) )
	IupSetAttribute(tgg1, "VALUE", "ON");

	f1 = IupFrame( IupVbox(boxdial1, tgg1) );
	IupSetAttribute(f1, "TITLE", "Y Zoom");

	# X zooming
	dial2 = IupDial("HORIZONTAL");
	lbl1 = IupLabel("-");
	lbl2 = IupLabel("+");
	boxinfo = IupHbox(lbl1, IupFill(), lbl2, C_NULL);
	boxdial2 = IupVbox(dial2, boxinfo, C_NULL);

	IupSetAttribute(boxdial2, "ALIGNMENT", "ACENTER");
	IupSetAttribute(boxinfo, "ALIGNMENT", "ACENTER");
	IupSetAttribute(boxinfo, "SIZE", "64x16");
	IupSetAttribute(boxinfo, "GAP", "2");
	IupSetAttribute(boxinfo, "MARGIN", "4");
	IupSetAttribute(boxinfo, "EXPAND", "HORIZONTAL");

	IupSetAttribute(lbl1, "EXPAND", "NO");
	IupSetAttribute(lbl2, "EXPAND", "NO");

	IupSetAttribute(dial2, "ACTIVE", "NO");
	IupSetAttribute(dial2, "SIZE", "64x16");
	IupSetCallback(dial2, "BUTTON_PRESS_CB", cfunction(dial2_btndown_cb, Int, (Ptr{Ihandle},)))
	IupSetCallback(dial2, "MOUSEMOVE_CB", cfunction(dial2_btnup_cb, Int, (Ptr{Ihandle}, Cdouble)) )
	IupSetCallback(dial2, "BUTTON_RELEASE_CB", cfunction(dial2_btnup_cb, Int, (Ptr{Ihandle}, Cdouble)) )

	tgg2 = IupToggle("X Autoscale");
	IupSetCallback(tgg2, "ACTION", cfunction(tgg2_cb, Int, (Ptr{Ihandle}, Cint)) )

	f2 = IupFrame( IupVbox(boxdial2, tgg2) )
	IupSetAttribute(f2, "TITLE", "X Zoom")

	lbl1 = IupLabel("")
	IupSetAttribute(lbl1, "SEPARATOR", "HORIZONTAL")

	tgg3 = IupToggle("Vertical Grid")
	IupSetCallback(tgg3, "ACTION", cfunction(tgg3_cb, Int, (Ptr{Ihandle}, Cint)) )
	tgg4 = IupToggle("Horizontal Grid");
	IupSetCallback(tgg4, "ACTION", cfunction(tgg4_cb, Int, (Ptr{Ihandle}, Cint)) )

	lbl2 = IupLabel("")
	IupSetAttribute(lbl2, "SEPARATOR", "HORIZONTAL");

	tgg5 = IupToggle("Legend");
	IupSetCallback(tgg5, "ACTION", cfunction(tgg5_cb, Int, (Ptr{Ihandle}, Cint)) )

	lbl3 = IupLabel("");
	IupSetAttribute(lbl3, "SEPARATOR", "HORIZONTAL");

	bt1 = IupButton("Export PDF");
	IupSetCallback(bt1, "ACTION", cfunction(bt1_cb, Int, (Ptr{Ihandle}, )) )

	vboxl = IupVbox(f1, f2, lbl1, tgg3, tgg4, lbl2, tgg5, lbl3, bt1, C_NULL);
	IupSetAttribute(vboxl, "GAP", "4");
	IupSetAttribute(vboxl, "EXPAND", "NO");

	# right panel: tabs with plots
	for (ii = 1:MAXPLOT)
		vboxr[ii] = IupVbox(plot[ii]); # each plot a tab
		s = @sprintf("Plot %d", ii)
		IupSetAttribute(vboxr[ii], "TABTITLE", s); # name each tab
	end
	vboxr[MAXPLOT+1] = C_NULL 	# mark end of vector

	tabs = IupTabsv(pointer(vboxr)) # create tabs
	IupSetCallback(tabs, "TABCHANGE_CB", cfunction(tabs_tabchange_cb, Int, (Ptr{Ihandle}, Ptr{Ihandle})) )

	# dialog
	hbox = IupHbox(vboxl, tabs);
	IupSetAttribute(hbox, "MARGIN", "4x4");
	IupSetAttribute(hbox, "GAP", "10");

	dlg = IupDialog(hbox);
	IupSetAttributes(dlg, "SIZE=500x240" );
	IupSetAttribute(dlg, "TITLE", "IupPlot Example");

	InitPlots(); # It must be able to be done independent of dialog Mapping

	tabs_tabchange_cb(tabs, vboxr[1]);

	IupShowXY(dlg, IUP_CENTER, IUP_CENTER);
	IupSetAttribute(dlg, "SIZE", C_NULL);

	IupMainLoop()
	IupClose()

	return 0

end

# --------------------------------------------------------------------------------
function delete_cb(ih::Ptr{Ihandle}, index::Cint, sample_index::Cint, x::Cdouble, y::Cdouble)
	s = @sprintf("DELETE_CB(%d, %d, %f, %f)\n", index, sample_index, x, y);
	println(s)
	return IUP_DEFAULT;
end

# --------------------------------------------------------------------------------
function select_cb(ih::Ptr{Ihandle}, index::Cint, sample_index::Cint, x::Cdouble, y::Cdouble, select::Cint)
	s = @sprintf("SELECT_CB(%d, %d, %f, %f, %d)\n", index, sample_index, x, y, select)
	println(s)
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function draw_cb(ih::Ptr{Ihandle}, index::Cint, sample_index::Cint, x::Cdouble, y::Cdouble, select::Cint)
	s = @sprintf("DRAWSAMPLE_CB(%d, %d, %f, %f, %d)\n", index, sample_index, x, y, select)
#	println(s)
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function edit_cb(ih::Ptr{Ihandle}, index::Cint, sample_index::Cint, x::Cdouble, y::Cdouble, new_x::Ptr{Float32},  new_y::Ptr{Float32})
	s = @sprintf("EDIT_CB(%d, %d, %f, %f, %f, %f)\n", index, sample_index, x, y, unsafe_load(new_x), unsafe_load(new_y))
	println(s)
	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function postdraw_cb(ih::Ptr{Ihandle}, cnv::Ptr{cdCanvas})
	x = pointer([0.0])
	y = pointer([0.0])

	IupPlotTransform(ih, 0.003, 0.02, x, y);
	x = int32(unsafe_load(x))
	y = int32(unsafe_load(y))

	CD_BOLD = 1 		# Very strange shit is happening here. If I use those below, it blow Julia right away.
	CD_SOUTH = 1
	cdCanvasFont(cnv, C_NULL, CD_BOLD, 10)
	cdCanvasTextAlignment(cnv, CD_SOUTH)
	cdCanvasText(cnv, x, y, pointer("My Inline Legend"))
	#println("POSTDRAW_CB()\n")

	return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function predraw_cb(ih::Ptr{Ihandle}, cnv::Ptr{cdCanvas})
	#println("PREDRAW_CB()\n");
	return IUP_DEFAULT;
end

# ---------------------------------------------------------------------------------------
function InitPlots()

	# PLOT 0 - MakeExamplePlot1
	IupSetAttribute(plot[1], "TITLE", "AutoScale");
	IupSetAttribute(plot[1], "MARGINTOP", "40")
	IupSetAttribute(plot[1], "MARGINLEFT", "40")
	IupSetAttribute(plot[1], "MARGINBOTTOM", "50")
	IupSetAttribute(plot[1], "TITLEFONTSIZE", "16")
	IupSetAttribute(plot[1], "AXS_XLABEL", "gnu (Foo)")
	IupSetAttribute(plot[1], "AXS_YLABEL", "Space (m^3)")
	IupSetAttribute(plot[1], "AXS_YFONTSIZE", "7");
	IupSetAttribute(plot[1], "AXS_YTICKFONTSIZE", "7");
	IupSetAttribute(plot[1], "LEGENDSHOW", "YES");
	IupSetAttribute(plot[1], "AXS_XFONTSIZE", "10");
	IupSetAttribute(plot[1], "AXS_YFONTSIZE", "10");
	IupSetAttribute(plot[1], "AXS_XLABELCENTERED", "NO");
	IupSetAttribute(plot[1], "AXS_YLABELCENTERED", "NO");

	#  IupSetAttribute(plot[1], "USE_IMAGERGB", "YES");
	#  IupSetAttribute(plot[1], "USE_GDI+", "YES");

	theFac = 1.0 / (100*100*100)
	IupPlotBegin(plot[1], 0);
	for (theI = -100:100)
		x = float64(theI+50);
		y = theFac*theI*theI*theI;
		IupPlotAdd(plot[1], x, y)
	end
	IupPlotEnd(plot[1]);
	IupSetAttribute(plot[1], "DS_LINEWIDTH", "3");
	IupSetAttribute(plot[1], "DS_LEGEND", "Line");

	theFac = 2.0/100;
	IupPlotBegin(plot[1], 0);
	for (theI = -100:100)
		x = float64(theI)
		y = -theFac*theI;
		IupPlotAdd(plot[1], x, y);
	end
	IupPlotEnd(plot[1]);
	IupSetAttribute(plot[1], "DS_LEGEND", "Curve 1");

	IupPlotBegin(plot[1], 0);
	for (theI = -100:100)
		x = (0.01*theI*theI-30);
		y = 0.01*theI;
		IupPlotAdd(plot[1], x, y);
	end
	IupPlotEnd(plot[1]);
	IupSetAttribute(plot[1], "DS_LEGEND", "Curve 2");

	# PLOT 1 - MakeExamplePlot2
	IupSetAttribute(plot[2], "TITLE", "No Autoscale+No CrossOrigin");
	IupSetAttribute(plot[2], "TITLEFONTSIZE", "16");
	IupSetAttribute(plot[2], "MARGINTOP", "40");
	IupSetAttribute(plot[2], "MARGINLEFT", "55");
	IupSetAttribute(plot[2], "MARGINBOTTOM", "50");
	IupSetAttribute(plot[2], "BGCOLOR", "0 192 192");
	IupSetAttribute(plot[2], "AXS_XLABEL", "Tg (X)");
	IupSetAttribute(plot[2], "AXS_YLABEL", "Tg (Y)");
	IupSetAttribute(plot[2], "AXS_XAUTOMIN", "NO");
	IupSetAttribute(plot[2], "AXS_XAUTOMAX", "NO");
	IupSetAttribute(plot[2], "AXS_YAUTOMIN", "NO");
	IupSetAttribute(plot[2], "AXS_YAUTOMAX", "NO");
	IupSetAttribute(plot[2], "AXS_XMIN", "10");
	IupSetAttribute(plot[2], "AXS_XMAX", "60");
	IupSetAttribute(plot[2], "AXS_YMIN", "-0.5");
	IupSetAttribute(plot[2], "AXS_YMAX", "0.5");
	IupSetAttribute(plot[2], "AXS_XCROSSORIGIN", "NO");
	IupSetAttribute(plot[2], "AXS_YCROSSORIGIN", "NO");
	IupSetAttribute(plot[2], "AXS_XFONTSTYLE", "BOLD");
	IupSetAttribute(plot[2], "AXS_YFONTSTYLE", "BOLD");
	IupSetAttribute(plot[2], "AXS_XREVERSE", "YES");
	IupSetAttribute(plot[2], "GRIDCOLOR", "128 255 128");
	IupSetAttribute(plot[2], "GRIDLINESTYLE", "DOTTED");
	IupSetAttribute(plot[2], "GRID", "YES");
	IupSetAttribute(plot[2], "LEGENDSHOW", "YES");

	theFac = 1.0/(100*100*100);
	IupPlotBegin(plot[2], 0);
	for (theI = 0:100)
		x = float64(theI);
		y = theFac*theI*theI*theI;
		IupPlotAdd(plot[2], x, y);
	end
	IupPlotEnd(plot[2]);

	theFac = 2.0/100;
	IupPlotBegin(plot[2], 0);
	for (theI = 0:100)
		x = float64(theI);
		y = -theFac*theI;
		IupPlotAdd(plot[2], x, y);
	end
	IupPlotEnd(plot[2]);

	# PLOT 2 - MakeExamplePlot4
	IupSetAttribute(plot[3], "TITLE", "Log Scale");
	IupSetAttribute(plot[3], "TITLEFONTSIZE", "16");
	IupSetAttribute(plot[3], "MARGINTOP", "40");
	IupSetAttribute(plot[3], "MARGINLEFT", "70");
	IupSetAttribute(plot[3], "MARGINBOTTOM", "50");
	IupSetAttribute(plot[3], "GRID", "YES");
	IupSetAttribute(plot[3], "AXS_XSCALE", "LOG10");
	IupSetAttribute(plot[3], "AXS_YSCALE", "LOG2");
	IupSetAttribute(plot[3], "AXS_XLABEL", "Tg (X)");
	IupSetAttribute(plot[3], "AXS_YLABEL", "Tg (Y)");
	IupSetAttribute(plot[3], "AXS_XFONTSTYLE", "BOLD");
	IupSetAttribute(plot[3], "AXS_YFONTSTYLE", "BOLD");

	theFac = 100.0/(100*100*100);
	IupPlotBegin(plot[3], 0);
	for (theI = 0:100)
		x = (0.0001+theI*0.001);
		y = (0.01+theFac*theI*theI*theI);
		IupPlotAdd(plot[3], x, y);
	end
	IupPlotEnd(plot[3]);
	IupSetAttribute(plot[3], "DS_COLOR", "100 100 200");

	# PLOT 3 - MakeExamplePlot5
	IupSetAttribute(plot[4], "TITLE", "Bar Mode");
	IupSetAttribute(plot[4], "TITLEFONTSIZE", "16");
	IupSetAttribute(plot[4], "MARGINTOP", "40");
	IupSetAttribute(plot[4], "MARGINLEFT", "30");
	IupSetAttribute(plot[4], "MARGINBOTTOM", "30");
	kLables = ["jan","feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"]
	kData = [1,2,3,4,5,6,7,8,9,0,1,2]

	IupPlotBegin(plot[4], 1);
	for (theI = 1:12)
		IupPlotAddStr(plot[4], kLables[theI], float64(kData[theI]))
	end
	IupPlotEnd(plot[4]);
	IupSetAttribute(plot[4], "DS_COLOR", "100 100 200");
	IupSetAttribute(plot[4], "DS_MODE", "BAR");

	#PLOT 4 - MakeExamplePlot6
	IupSetAttribute(plot[5], "TITLE", "Marks Mode");
	IupSetAttribute(plot[5], "TITLEFONTSIZE", "16");
	IupSetAttribute(plot[5], "MARGINTOP", "40");
	IupSetAttribute(plot[5], "MARGINLEFT", "45");
	IupSetAttribute(plot[5], "MARGINBOTTOM", "40");
	IupSetAttribute(plot[5], "AXS_XAUTOMIN", "NO");
	IupSetAttribute(plot[5], "AXS_XAUTOMAX", "NO");
	IupSetAttribute(plot[5], "AXS_YAUTOMIN", "NO");
	IupSetAttribute(plot[5], "AXS_YAUTOMAX", "NO");
	IupSetAttribute(plot[5], "AXS_XMIN", "0");
	IupSetAttribute(plot[5], "AXS_XMAX", "0.011");
	IupSetAttribute(plot[5], "AXS_YMIN", "0");
	IupSetAttribute(plot[5], "AXS_YMAX", "0.22");
	IupSetAttribute(plot[5], "AXS_XCROSSORIGIN", "NO");
	IupSetAttribute(plot[5], "AXS_YCROSSORIGIN", "NO");
	IupSetAttribute(plot[5], "AXS_XTICKFORMAT", "%1.3f");
	IupSetAttribute(plot[5], "LEGENDSHOW", "YES");
	IupSetAttribute(plot[5], "LEGENDPOS", "BOTTOMRIGHT");

	theFac = 100.0/(100*100*100);
	IupPlotBegin(plot[5], 0);
	for (theI = 0:10)
		x = (0.0001+theI*0.001);
		y = (0.01+theFac*theI*theI);
		IupPlotAdd(plot[5], x, y);
	end
	IupPlotEnd(plot[5]);
	IupSetAttribute(plot[5], "DS_MODE", "MARKLINE");
	IupSetAttribute(plot[5], "DS_SHOWVALUES", "YES");

	IupPlotBegin(plot[5], 0);
	for (theI = 0:10)
		x = (0.0001+theI*0.001);
		y = (0.2-theFac*theI*theI);
		IupPlotAdd(plot[5], x, y);
	end
	IupPlotEnd(plot[5]);
	IupSetAttribute(plot[5], "DS_MODE", "MARK");
	IupSetAttribute(plot[5], "DS_MARKSTYLE", "HOLLOW_CIRCLE");

	# PLOT 5 - MakeExamplePlot8
	IupSetAttribute(plot[6], "TITLE", "Data Selection and Editing");
	IupSetAttribute(plot[6], "TITLEFONTSIZE", "16");
	IupSetAttribute(plot[6], "MARGINTOP", "40");

	theFac = 100.0/(100*100*100);
	IupPlotBegin(plot[6], 0);
	for (theI = -10:10)
		x = (0.001*theI);
		y = (0.01+theFac*theI*theI*theI);
		IupPlotAdd(plot[6], x, y);
	end
	IupPlotEnd(plot[6]);
	IupSetAttribute(plot[6], "DS_COLOR", "100 100 200");
#	IupSetAttribute(plot[6], "DS_EDIT", "YES");
	IupSetCallback(plot[6], "DELETE_CB", cfunction(delete_cb, Int, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble)) )
	IupSetCallback(plot[6], "SELECT_CB", cfunction(select_cb, Int, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble, Cint)) )
	IupSetCallback(plot[6], "CLICKSAMPLE_CB", cfunction(select_cb, Int, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble, Cint)) )
	IupSetCallback(plot[6], "POSTDRAW_CB", cfunction(postdraw_cb, Int, (Ptr{Ihandle}, Ptr{cdCanvas})) )
	IupSetCallback(plot[6], "PREDRAW_CB", cfunction(predraw_cb, Int, (Ptr{Ihandle}, Ptr{cdCanvas})) )
#	IupSetCallback(plot[6], "EDIT_CB", cfunction(edit_cb, Int, (Ptr{Ihandle}, Cint, Cint, Cfloat, Cfloat, Ptr{Cfloat},  Ptr{Cfloat})) )
	IupSetCallback(plot[6], "DRAWSAMPLE_CB", cfunction(draw_cb, Int, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble, Cint)) )
end

# ------------------------------------------------------------------------------
# Some processing required by current tab change: the controls at left
# will be updated according to current plot props
function tabs_tabchange_cb(self::Ptr{Ihandle}, new_tab::Ptr{Ihandle})

global dial1, dial2
global tgg1, tgg2
global tgg3, tgg4
global tgg5
global tabs

  ss = bytestring(IupGetAttribute(new_tab, "TABTITLE"))
  #ss += 5; # Skip "Plot "
  #ii = atoi(ss);
  ii = parseint(ss[6]) 		# Where the hell is atoi() ????

  # autoscaling
  # X axis
  if (IupGetInt(plot[ii], "AXS_XAUTOMIN") != 0 && IupGetInt(plot[ii], "AXS_XAUTOMAX") != 0)
    IupSetAttribute(tgg2, "VALUE", "ON");
    IupSetAttribute(dial2, "ACTIVE", "NO");
  else
    IupSetAttribute(tgg2, "VALUE", "OFF");
    IupSetAttribute(dial2, "ACTIVE", "YES");
  end
  # Y axis
  if ( IupGetInt(plot[ii], "AXS_YAUTOMIN") != 0 && IupGetInt(plot[ii], "AXS_YAUTOMAX") != 0)
    IupSetAttribute(tgg1, "VALUE", "ON");
    IupSetAttribute(dial1, "ACTIVE", "NO");
  else
    IupSetAttribute(tgg1, "VALUE", "OFF");
    IupSetAttribute(dial1, "ACTIVE", "YES");
  end

  # grid
  if (IupGetInt(plot[ii], "GRID") != 0)
    IupSetAttribute(tgg3, "VALUE", "ON");
    IupSetAttribute(tgg4, "VALUE", "ON");
  else
    # X axis
    if (IupGetAttribute(plot[ii], "GRID") == 'V')
      IupSetAttribute(tgg3, "VALUE", "ON");
    else
      IupSetAttribute(tgg3, "VALUE", "OFF")
    end
    # Y axis
    if (IupGetAttribute(plot[ii], "GRID") == 'H')
      IupSetAttribute(tgg4, "VALUE", "ON");
    else
      IupSetAttribute(tgg4, "VALUE", "OFF")
    end
  end

  # legend
  if (IupGetInt(plot[ii], "LEGENDSHOW") != 0)
    IupSetAttribute(tgg5, "VALUE", "ON");
  else
    IupSetAttribute(tgg5, "VALUE", "OFF")
  end

  return IUP_DEFAULT
end

# show/hide V grid
function tgg3_cb(self::Ptr{Ihandle}, v::Cint)
	ii = tabs_get_index()

	if (v != 0)
		if (IupGetInt(tgg4, "VALUE") != 0)
			IupSetAttribute(plot[ii], "GRID", "YES");
		else
			IupSetAttribute(plot[ii], "GRID", "VERTICAL")
		end
	else
		if (IupGetInt(tgg4, "VALUE") != 0)
			IupSetAttribute(plot[ii], "GRID", "NO");
		else
			IupSetAttribute(plot[ii], "GRID", "HORIZONTAL")
		end
	end

	IupSetAttribute(plot[ii], "REDRAW");

	return IUP_DEFAULT
end

# show/hide H grid
function tgg4_cb(self::Ptr{Ihandle}, v::Cint)
	ii = tabs_get_index()

	if (v != 0)
		if (IupGetInt(tgg3, "VALUE") != 0)
			IupSetAttribute(plot[ii], "GRID", "YES")
		else
			IupSetAttribute(plot[ii], "GRID", "HORIZONTAL")
		end
	else
		if (IupGetInt(tgg3, "VALUE") != 0)
			IupSetAttribute(plot[ii], "GRID", "NO")
		else
			IupSetAttribute(plot[ii], "GRID", "VERTICAL")
		end
	end

	IupSetAttribute(plot[ii], "REDRAW")

	return IUP_DEFAULT
end

# ---------------------------------------------------------------------------
function tabs_get_index()
	curr_tab = IupGetHandle(IupGetAttribute(tabs, "VALUE"));
	ss = bytestring(IupGetAttribute(curr_tab, "TABTITLE"))
	ii = parseint(ss[6]) 		# Where the hell is atoi() ????
	return ii
end

# Some processing required by current tab change: the controls at left
# will be updated according to current plot props
function tabs_tabchange_cb(self::Ptr{Ihandle}, new_tab::Ptr{Ihandle})

	ss = bytestring(IupGetAttribute(new_tab, "TABTITLE"))
	ii = parseint(ss[6]) 		# Where the hell is atoi() ????

	# autoscaling
	# X axis
	if (IupGetInt(plot[ii], "AXS_XAUTOMIN") != 0 && IupGetInt(plot[ii], "AXS_XAUTOMAX") != 0)
		IupSetAttribute(tgg2, "VALUE", "ON");
		IupSetAttribute(dial2, "ACTIVE", "NO");
	else
		IupSetAttribute(tgg2, "VALUE", "OFF");
		IupSetAttribute(dial2, "ACTIVE", "YES");
	end
  # Y axis
  if (IupGetInt(plot[ii], "AXS_YAUTOMIN") != 0 && IupGetInt(plot[ii], "AXS_YAUTOMAX") != 0)
    IupSetAttribute(tgg1, "VALUE", "ON");
    IupSetAttribute(dial1, "ACTIVE", "NO");
  else
    IupSetAttribute(tgg1, "VALUE", "OFF");
    IupSetAttribute(dial1, "ACTIVE", "YES");
  end

  # grid
  if (IupGetInt(plot[ii], "GRID") != 0)
    IupSetAttribute(tgg3, "VALUE", "ON");
    IupSetAttribute(tgg4, "VALUE", "ON");
  else
    # X axis
    s = bytestring(IupGetAttribute(plot[ii], "GRID"))
    if (s[1] == 'V')
      IupSetAttribute(tgg3, "VALUE", "ON");
    else
      IupSetAttribute(tgg3, "VALUE", "OFF")
    end
    # Y axis
    if (s[1] == 'H')
      IupSetAttribute(tgg4, "VALUE", "ON");
    else
      IupSetAttribute(tgg4, "VALUE", "OFF")
    end
  end

  # legend
  if (IupGetInt(plot[ii], "LEGENDSHOW") != 0)
    IupSetAttribute(tgg5, "VALUE", "ON");
  else
    IupSetAttribute(tgg5, "VALUE", "OFF")
end

  return IUP_DEFAULT
end

# show/hide V grid
# ----------------------------------------------------------------------------------
function tgg3_cb(self::Ptr{Ihandle}, v)
  ii = tabs_get_index();

  if (v != 0)
    if (IupGetInt(tgg4, "VALUE"))
      IupSetAttribute(plot[ii], "GRID", "YES");
    else
      IupSetAttribute(plot[ii], "GRID", "VERTICAL")
    end
  else
    if (!IupGetInt(tgg4, "VALUE"))
      IupSetAttribute(plot[ii], "GRID", "NO");
    else
      IupSetAttribute(plot[ii], "GRID", "HORIZONTAL")
    end
  end

  IupSetAttribute(plot[ii], "REDRAW");

  return IUP_DEFAULT
end

# show/hide H grid
# ----------------------------------------------------------------------------------
function tgg4_cb(self::Ptr{Ihandle}, v)
  ii = tabs_get_index();

  if (v != 0)
    if (IupGetInt(tgg3, "VALUE"))
      IupSetAttribute(plot[ii], "GRID", "YES");
    else
      IupSetAttribute(plot[ii], "GRID", "HORIZONTAL");
    end
  else
    if (!IupGetInt(tgg3, "VALUE"))
      IupSetAttribute(plot[ii], "GRID", "NO");
    else
      IupSetAttribute(plot[ii], "GRID", "VERTICAL")
    end
  end

  IupSetAttribute(plot[ii], "REDRAW");

  return IUP_DEFAULT
end

# show/hide legend
# ----------------------------------------------------------------------------------
function tgg5_cb(self::Ptr{Ihandle}, v::Cint)
  ii = tabs_get_index()

  if (v != 0)
    IupSetAttribute(plot[ii], "LEGENDSHOW", "YES");
  else
    IupSetAttribute(plot[ii], "LEGENDSHOW", "NO")
  end

  IupSetAttribute(plot[ii], "REDRAW");

  return IUP_DEFAULT;
end

# autoscale Y
# ----------------------------------------------------------------------------------
function tgg1_cb(self::Ptr{Ihandle}, v::Cint)
  ii = tabs_get_index();

  if (v != 0)
    IupSetAttribute(dial1, "ACTIVE", "NO");
    IupSetAttribute(plot[ii], "AXS_YAUTOMIN", "YES");
    IupSetAttribute(plot[ii], "AXS_YAUTOMAX", "YES");
  else
    IupSetAttribute(dial1, "ACTIVE", "YES");
    IupSetAttribute(plot[ii], "AXS_YAUTOMIN", "NO");
    IupSetAttribute(plot[ii], "AXS_YAUTOMAX", "NO");
  end

  IupSetAttribute(plot[ii], "REDRAW");

  return IUP_DEFAULT;
end

# autoscale X
# ----------------------------------------------------------------------------------
function tgg2_cb(self::Ptr{Ihandle}, v::Cint)
  ii = tabs_get_index();

  if (v != 0)
    IupSetAttribute(dial2, "ACTIVE", "NO");
    IupSetAttribute(plot[ii], "AXS_XAUTOMIN", "YES");
    IupSetAttribute(plot[ii], "AXS_XAUTOMAX", "YES");
  else
    IupSetAttribute(dial2, "ACTIVE", "YES");
    IupSetAttribute(plot[ii], "AXS_XAUTOMIN", "NO");
    IupSetAttribute(plot[ii], "AXS_XAUTOMAX", "NO");
  end

  IupSetAttribute(plot[ii], "REDRAW");

  return IUP_DEFAULT;
end

# Y zoom
# ----------------------------------------------------------------------------------
function dial1_btndown_cb(self::Ptr{Ihandle})
	ii = tabs_get_index()
	IupStoreAttribute(plot[ii], "OLD_YMIN", IupGetAttribute(plot[ii], "AXS_YMIN"));
	IupStoreAttribute(plot[ii], "OLD_YMAX", IupGetAttribute(plot[ii], "AXS_YMAX"));
	return IUP_DEFAULT;
end

# ----------------------------------------------------------------------------------
function dial1_btnup_cb(self::Ptr{Ihandle}, angle::Cdouble)
	ii = tabs_get_index();

	x1 = IupGetFloat(plot[ii], "OLD_YMIN");
	x2 = IupGetFloat(plot[ii], "OLD_YMAX");

	ss = IupGetAttribute(plot[ii], "AXS_YMODE");
	if ( ss != C_NULL && ss[3] == '2' )
		# LOG2:  one circle will zoom 2 times
		xm = 4.0 * fabs(angle) / 3.141592;
		if (angle>0.0)
			x2 /= xm; x1 *= xm;
		else
			x2 *= xm; x1 /= xm
		end
	end
	if ( ss != C_NULL && ss[3] =='1' )
		# LOG10:  one circle will zoom 10 times
		xm = 10.0 * fabs(angle) / 3.141592;
		if (angle > 0.0)
			x2 /= xm; x1 *= xm
		else
			x2 *= xm; x1 /= xm
		end
	else
		# LIN: one circle will zoom 2 times
		xm = (x1 + x2) / 2.0;
		x1 = xm - (xm - x1)*(1.0-angle*1.0/3.141592);
		x2 = xm + (x2 - xm)*(1.0-angle*1.0/3.141592);
	end

	if (x1 < x2)
		s = @sprintf("%f", x1)
		IupSetAttribute(plot[ii], "AXS_YMIN", s);
		s = @sprintf("%f", x2)
		IupSetAttribute(plot[ii], "AXS_YMAX", s);
	end

	IupSetAttribute(plot[ii], "REDRAW");

	return IUP_DEFAULT;
end

# X zoom
# ----------------------------------------------------------------------------------
function dial2_btndown_cb(self::Ptr{Ihandle})
	ii = tabs_get_index();
	IupStoreAttribute(plot[ii], "OLD_XMIN", IupGetAttribute(plot[ii], "AXS_XMIN"));
	IupStoreAttribute(plot[ii], "OLD_XMAX", IupGetAttribute(plot[ii], "AXS_XMAX"));
	return IUP_DEFAULT;
end

# ----------------------------------------------------------------------------------
function dial2_btnup_cb(self::Ptr{Ihandle}, angle::Cdouble)
	ii = tabs_get_index();

	x1 = IupGetFloat(plot[ii], "OLD_XMIN");
	x2 = IupGetFloat(plot[ii], "OLD_XMAX");
	xm = (x1 + x2) / 2.0;
	x1 = xm - (xm - x1)*(1.0-angle*1.0/3.141592);		# one circle will zoom 2 times
	x2 = xm + (x2 - xm)*(1.0-angle*1.0/3.141592);

	s = @sprintf("%f", x1)
	IupSetAttribute(plot[ii], "AXS_XMIN", s);
	s = @sprintf("%f", x2)
	IupSetAttribute(plot[ii], "AXS_XMAX", s);
	IupSetAttribute(plot[ii], "REDRAW");

	return IUP_DEFAULT
end

# ----------------------------------------------------------------------------------
function bt1_cb(self::Ptr{Ihandle})
	ii = tabs_get_index()
	#cnv = cdCreateCanvas(cdContextPDF(), "pplot.pdf -o");
	#IupPlotPaintTo(plot[ii], cnv);
	#cdKillCanvas(cnv);
	return IUP_DEFAULT;
end

end 	# end module


text:
@testset "Check identifiability from given IO-equations" begin
    R, (x, y, z, t) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z", "t"])

    @test check_identifiability(
        x * (y + z) * t^4 - x^3 * y * z * t^5 + 7 + t * x,
        [x, y, z],
        0.99
    ) == [true, false, false]

    @test check_identifiability(
        [z + t * (x + y) + z^3, z - t^2 * x * y *5 + t^2],
        [x, y],
        [x, y, x + y, x - y, x^2 + y^2],
        0.99
    ) == [false, false, true, false, true]

end

ext: jl


text:
using DashBootstrapComponents, DashHtmlComponents

toast = dbc_toast(
    [html_p("This is the content of the toast", className="mb-0")],
    header="This is the header",
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using InteractionNetworkModels, Distributions, BenchmarkTools, Plots
using Distances, StructuredDistances
# The Model(s)
model_mode = Hollywood(-3.0, Poisson(7), 10)
S = sample(model_mode, 10)
S = [[1,1,1,1], [2,2,2,2], [3,3,3,3,3,3]]
V = 1:20

d_lcs = MatchingDistance(FastLCS(100))
d_lcs_t = ThreadedMatchingDistance(LCS())
d_lsp = MatchingDistance(FastLSP(100))
# d_f = FastMatchingDist(FastLCS(100), 51)
d_lsp_fp = FpMatchingDist(FastLSP(100), 100.0)

d_c = CouplingDistance(FastLSP(100))

K_inner = DimensionRange(1, 50)
K_outer = DimensionRange(1, 50)

d_lsp_as = AvgSizeMatchingDist(FastLSP(100), 0.5)
d_lsp_sc = SizeConstrainedMatchingDist(FastLSP(100), 1.0, 4)

model = SIM(S, 5.0, d_lcs, V, K_inner, K_outer)
model_t = SIM(S, 5.0, d_lcs_t, V, K_inner, K_outer)

@time d_lsp_as([[1,2]], [[1,2], [1,3,3,4,5]])
@time d_lsp_sc([[1,2]], [[1,2], [1,3,3,4,5,12,3]])
@time d_lsp([[1,2]], [[1,2], [1,3,3,4,5,12,3]])

# model_f = SIM(S, 4.0, d_f, V, 50, 50)
mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=5, β=0.6, ν_td=3,
    len_dist=TrGeometric(0.8, K_inner.l, K_inner.u),
    lag=1,
    K=200,
    burn_in=1000
)
mcmc_sampler_sp = SimMcmcInsertDeleteSubpath(
    ν_ed=5, β=0.6, ν_td=3,
    len_dist=TrGeometric(0.1, model.K_inner.l, model.K_inner.u),
    lag=1,
    K=200
)

mcmc_sampler_prop = SimMcmcInsertDeleteProportional(
    ν_ed=1, β=0.6, ν_td=3,
    len_dist=TrGeometric(0.1, model.K_inner.l, model.K_inner.u),
    lag=1,
    K=200, burn_in=1000
)

mcmc_sampler_len = SimMcmcInsertDeleteLengthCentered(
    ν_ed=1, β=0.6, ν_td=3,
    lag=1,
    K=200, burn_in=1000
)

@time out=mcmc_sampler(
    model,
    lag=20,
    init=model.mode,
    burn_in=0,
    desired_samples=200
)

@time out=mcmc_sampler(
    model_t,
    lag=20,
    init=model.mode,
    burn_in=0,
    desired_samples=200
)

plot(out)
summaryplot(out)
out.sample[100]

d_lsp_sc(, model.mode)

d_lsp_sc([[1,1,1], [2,2,2]], [[1,1,1],[2,2,2],[3,4,4,4]])

@time out=mcmc_sampler_len(
    model,
    lag=1,
    init=model.mode,
    burn_in=0,
    desired_samples=10000
)
plot(out)
summaryplot(out)
out.sample

@time out_prop=mcmc_sampler_prop(
    model,
    lag=1,
    init=model.mode,
    desired_samples=10000
)
plot(out_prop)
summaryplot(out_prop)
out_prop.sample

sample_frechet_var(out.sample, d_lcs, with_memory=true)

n = 1000
m = 4
samples = [draw_sample(mcmc_sampler, model, desired_samples=n) for i in 1:m]

mean_dists_summary(samples, d_lcs)

@time out_sp=mcmc_sampler_sp(
    model,
    lag=1,
    init=model.mode,
    desired_samples=10000,
    burn_in=0
)
plot!(out_sp)

summaryplot(out)
out_sp.sample[1:100]

@btime out=mcmc_sampler(
    model_f,
    lag=1,
    init=model.mode,
    desired_samples=2000,
    burn_in=0
)

plot(out)
summaryplot(out)
S
out.sample

tmp = [1,2,1,2,1,4]
V = 1:10

p = zeros(length(V))
p .+= [i ∈ tmp for i in V]
p ./= sum(p)


text:
using Flux
using Plots
using DifferentialEquations

# solve 2nd order ode numerically
k = 1.0
force(dx,x,k,t) = -k*x + 0.1sin(x)
prob = SecondOrderODEProblem(force,1.0,0.0,(0.0,10.0),k)
sol = solve(prob)
plot(sol,label=["Velocity" "Position"])

# Generate the dataset
plot_t = 0:0.01:10
data_plot = sol(plot_t)
positions_plot = [state[2] for state in data_plot]
force_plot = [force(state[1],state[2],k,t) for state in data_plot]

t = 0:3.3:10
dataset = sol(t)
position_data = [state[2] for state in sol(t)]
force_data = [force(state[1],state[2],k,t) for state in sol(t)]

plot(plot_t,force_plot,xlabel="t",label="True Force")
scatter!(t,force_data,label="Force Measurements")

# create NN
NNForce = Chain(x -> [x],
                Dense(1,32,tanh),
                Dense(32,1),
                first)

loss() = sum(abs2,NNForce(position_data[i]) - force_data[i] for i in 1:length(position_data))
loss()

opt = Flux.Descent(0.01)
data = Iterators.repeated((), 5000)
iter = 0
cb = function () #callback function to observe training
    global iter += 1
    if iter % 500 == 0
        display(loss())
    end
end
display(loss())

# train the NN
Flux.train!(loss, Flux.params(NNForce), data, opt; cb=cb)

# visualize the result
learned_force_plot = NNForce.(positions_plot)
plot(plot_t,force_plot,xlabel="t",label="True Force")
plot!(plot_t,learned_force_plot,label="Predicted Force")
scatter!(t,force_data,label="Force Measurements")

# solve SHO problem to justify PINN
force2(dx,x,k,t) = -k*x
prob_simplified = SecondOrderODEProblem(force2,1.0,0.0,(0.0,10.0),k)
sol_simplified = solve(prob_simplified)
plot(sol,label=["Velocity" "Position"])
plot!(sol_simplified,label=["Velocity Simplified" "Position Simplified"])

# create Physics Loss Term to force solution to satisfy SHO ode as well
random_positions = [2rand()-1 for i in 1:100] # random values in [-1,1]
loss_ode() = sum(abs2,NNForce(x) - (-k*x) for x in random_positions)
loss_ode()

λ = 0.1  # factor controlling physics regularization
composed_loss() = loss() + λ*loss_ode()

opt = Flux.Descent(0.01)
data = Iterators.repeated((), 5000)
iter = 0
cb = function () #callback function to observe training
    global iter += 1
    if iter % 500 == 0
        display(composed_loss())
    end
end
display(composed_loss())

# train the new model
Flux.train!(composed_loss, Flux.params(NNForce), data, opt; cb=cb)

# visualize the final result!
learned_force_plot = NNForce.(positions_plot)

plot(plot_t,force_plot,xlabel="t",label="True Force")
plot!(plot_t,learned_force_plot,label="Predicted Force")
scatter!(t,force_data,label="Force Measurements")


text:
# Autogenerated wrapper script for libwebp_jll for x86_64-apple-darwin
export cwebp, dwebp, gif2webp, img2webp, libwebp, libwebpdecoder, libwebpdemux, libwebpmux, webpinfo, webpmux

using Giflib_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("libwebp")
JLLWrappers.@declare_executable_product(cwebp)
JLLWrappers.@declare_executable_product(dwebp)
JLLWrappers.@declare_executable_product(gif2webp)
JLLWrappers.@declare_executable_product(img2webp)
JLLWrappers.@declare_library_product(libwebp, "@rpath/libwebp.7.dylib")
JLLWrappers.@declare_library_product(libwebpdecoder, "@rpath/libwebpdecoder.3.dylib")
JLLWrappers.@declare_library_product(libwebpdemux, "@rpath/libwebpdemux.2.dylib")
JLLWrappers.@declare_library_product(libwebpmux, "@rpath/libwebpmux.3.dylib")
JLLWrappers.@declare_executable_product(webpinfo)
JLLWrappers.@declare_executable_product(webpmux)
function __init__()
    JLLWrappers.@generate_init_header(Giflib_jll, JpegTurbo_jll, libpng_jll, Libtiff_jll, Libglvnd_jll)
    JLLWrappers.@init_executable_product(
        cwebp,
        "bin/cwebp",
    )

    JLLWrappers.@init_executable_product(
        dwebp,
        "bin/dwebp",
    )

    JLLWrappers.@init_executable_product(
        gif2webp,
        "bin/gif2webp",
    )

    JLLWrappers.@init_executable_product(
        img2webp,
        "bin/img2webp",
    )

    JLLWrappers.@init_library_product(
        libwebp,
        "lib/libwebp.7.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpdecoder,
        "lib/libwebpdecoder.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpdemux,
        "lib/libwebpdemux.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpmux,
        "lib/libwebpmux.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        webpinfo,
        "bin/webpinfo",
    )

    JLLWrappers.@init_executable_product(
        webpmux,
        "bin/webpmux",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module CURAND

using ..CuArrays
using ..CuArrays: libcurand, active_context

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

import CUDAnative

using CEnum

include("libcurand_common.jl")
include("error.jl")

version() = VersionNumber(curandGetProperty(CUDAapi.MAJOR_VERSION),
                          curandGetProperty(CUDAapi.MINOR_VERSION),
                          curandGetProperty(CUDAapi.PATCH_LEVEL))

include("libcurand.jl")
include("wrappers.jl")
include("random.jl")

const _generators = Dict{CuContext,RNG}()
const _generator = Ref{Union{Nothing,RNG}}(nothing)

function generator()
    if _generator[] == nothing
        CUDAnative.maybe_initialize("CURAND")
        _generator[] = get!(_generators, active_context[]) do
            context = active_context[]
            RNG()
        end
    end

    return _generator[]::RNG
end

end

const seed! = CURAND.seed!
const rand = CURAND.rand
const randn = CURAND.randn
const rand_logn = CURAND.rand_logn
const rand_poisson = CURAND.rand_poisson

@deprecate curand CuArrays.rand
@deprecate curandn CuArrays.randn


text:
normalize(l::Lens) = l
normalize(l::ComposedLens) = _compose(l.outer, l.inner)

"""
    _compose(lens1, lens2)

Like `∘` but fixes the associativity to match with the default one in
Setfield.
"""
_compose(l1::Lens, l2::Lens) = l1 ∘ l2
_compose(l1::Lens, l2::ComposedLens) =
    _compose(_compose(normalize(l1), l2.outer), l2.inner)


text:
using ExactReals
using Test

@testset "ExactReals.jl" begin

    @testset "Constructors" begin
        a = ExactReal(3//5)
        @test value_range(a) == (0.6, nextfloat(0.6))

        b = ExactReal(3)
        @test value_range(b) == (prevfloat(3.0), nextfloat(3.0))
    end

    @testset "Arithmetic" begin

        a = ExactReal(3//5)

        b = a + a
        @test value_range(b) == (1.2, nextfloat(1.2))

        c = a * a
        @test value_range(c) == (0.36, nextfloat(0.36))

        d = -a
        @test value_range(d) == (prevfloat(-0.6), -0.6)

        e = sqrt(a)
        s = sqrt(0.6)
        @test value_range(e) == (prevfloat(s), s)
    end

    @testset "Comparisons" begin
        a = ExactReal(3//5)
        @test ExactReals.compare(1 - a, a) == -1

        @test a > ExactReal(1//2)
    end

end


text:
#
# SpherePlots.jl –
#
# import RecipesBase: plot
# import Plots: scatter3d, scatter3d!, surface!, surface, current
# export plot
# """
#     plot(M::Sphere, signals) – plot signals on the two sphere
# """
# function plot(M::Sphere,signals::Array{Array{SnPoint,1},1})
#     if manifoldDimension(M) > 2
#         error("Plot only works for functions on S2")
#     end
#     plotlyjs()
#     n = 180
#     u = linspace(0,2*π,n);
#     v = linspace(0,π,n);
#     x = cos.(u) * sin.(v)';
#     y = sin.(u) * sin.(v)';
#     z = ones(n) * cos.(v)';
#     surface(x,y,z,opacity=.95,
#     f = :white,
#     axis=nothing,colorbar=nothing,legend=nothing,border=:none)
#     for signal in signals
#         if any(manifoldDimension.(signal) .> 2)
#             error("This plot method only works for points on S2.");
#         end
#         xS = [ getValue(p)[1] for p in signal]
#         yS = [ getValue(p)[2] for p in signal]
#         zS = [ getValue(p)[3] for p in signal]
#         scatter3d!(xS,yS,zS)
#     end
#     return current()
# end
# plot(M::Sphere, signal::Array{SnPoint,1}) = plot(M,[signal])


text:
# Pseudocounts
# ============

"Parametric abstract type to define pseudocount types"
abstract type Pseudocount{T<:Real} end

"""
You can use `NoPseudocount()` to avoid pseudocount corrections where a
`Pseudocount` type is needed.
"""
struct NoPseudocount <: Pseudocount{Float64} end

"""
**Additive Smoothing** or fixed pseudocount `λ` for `ResidueCount`
(in order to estimate probabilities when the number of samples is low).

Common values of `λ` are:

- `0` : No cell frequency prior, gives you the maximum likelihood estimator.
- `0.05` is the optimum value for `λ` found in Buslje et. al. 2009, similar results was obtained for `λ` in the range [0.025, 0.075].
- `1 / p` : Perks prior (Perks, 1947) where `p` the number of parameters (i.e. residues, pairs of residues) to estimate. If `p` is the number of residues (`20` without counting gaps), this gives you `0.05`.
- `sqrt(n) / p` : Minimax prior (Trybula, 1958) where `n` is the number of samples and `p` the number of parameters to estimate. If the number of samples `n` is 400 (minimum number of sequence clusters for achieve good performance in Buslje et. al. 2009) for estimating 400 parameters (pairs of residues without counting gaps) this gives you `0.05`.
- `0.5` : Jeffreys prior (Jeffreys, 1946).
- `1` : Bayes-Laplace uniform prior, aka. Laplace smoothing.
"""
struct AdditiveSmoothing{T} <: Pseudocount{T}
    λ::T
end

Base.zero(::Type{AdditiveSmoothing{T}}) where {T} = AdditiveSmoothing(zero(T))
Base.one(::Type{AdditiveSmoothing{T}}) where {T} = AdditiveSmoothing(one(T))


text:
using LabelledArrays

export ar1_bidir

"""
    eom_ar1_bidir(x, p, n) → SVector{2}

Equations of motion for a system consisting of two mutually
coupled first order autoregressive processes.

## Equations of motion

```math
\\begin{aligned}
x(t+1) &= a_{1}x + c_{yx}y + \\epsilon_{x} \\\\
y(t+1) &= b_{1}y + c_{xy}x + \\epsilon_{y}
\\end{aligned}
```

where ``\\epsilon_{x}`` and ``\\epsilon_{y}`` are drawn independently
at each time step from normal distributions with zero mean and standard
deviations `σx` and `σy`.
"""
function eom_ar1_bidir(x, p, n)
    a₁, b₁, c_xy, c_yx, σx, σy = (p...,)
    x, y = (x...,)

    ϵx = rand(Normal(0, σx))
    ϵy = rand(Normal(0, σy))

    dx = a₁*x + c_yx*y + ϵx
    dy = b₁*y + c_xy*x + ϵy
    return SVector{2}(dx, dy)
end

function ar1_bidir(u₀,a₁, b₁, c_xy, c_yx, σx, σy)
    p = @LArray [a₁, b₁, c_xy, c_yx, σx, σy] (:a₁, :b₁, :c_xy, :c_yx, :σx, :σy)
    logistic_system = DiscreteDynamicalSystem(eom_ar1_bidir, u₀, p)

    return logistic_system
end

"""
    ar1_bidir(;u₀ = rand(2), a₁ = 0.5, b₁ = 0.7, c_xy = 0, c_yx = 0.2,
        σx = 0.3, σy = 0.3) → DiscreteDynamicalSystem

A system consisting of two mutually coupled first order autoregressive processes.

## Equations of motion

```math
\\begin{aligned}
x(t+1) &= a_{1}x + c_{yx}y + \\epsilon_{x} \\\\
y(t+1) &= b_{1}y + c_{xy}x + \\epsilon_{y}
\\end{aligned}
```

where ``\\epsilon_{x}`` and ``\\epsilon_{y}`` are drawn independently
at each time step from normal distributions with zero mean and standard
deviations `σx` and `σy`.
"""
ar1_bidir(;a₁ = 0.5, b₁ = 0.7, u₀ = rand(2), c_xy = 0, c_yx = 0.2, σx = 0.3, σy = 0.3) =
    ar1_bidir(u₀, a₁, b₁, c_xy, c_yx, σx, σy)


text:
let user_data_dir
    ENV["PATH"] = JULIA_HOME*";"*joinpath(JULIA_HOME,"..","Git","bin")*";"*ENV["PATH"]
    haskey(ENV,"JULIA_EDITOR") || (ENV["JULIA_EDITOR"] = "start")
    user_data_dir = abspath(ENV["AppData"],"julia")
    isdir(user_data_dir) || mkdir(user_data_dir)
    ispty(s::Base.AsyncStream) = bool(ccall(:jl_ispty, Cint, (Ptr{Void},), s))
    ispty(s) = false
    if ispty(STDIN)
        function Base.Terminals.raw!(t::Base.Terminals.UnixTerminal,raw::Bool)
            if raw
                run(`stty raw -echo`,t.in_stream,t.out_stream,t.err_stream)
            else
                run(`stty cooked echo`,t.in_stream,t.out_stream,t.err_stream)
            end
            true
        end
    end
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
module DiffResults

using StaticArrays

#########
# Types #
#########

abstract type DiffResult{O,V,D<:Tuple} end

struct ImmutableDiffResult{O,V,D<:Tuple} <: DiffResult{O,V,D}
    value::V
    derivs::D # ith element = ith-order derivative
    function ImmutableDiffResult(value::V, derivs::NTuple{O,Any}) where {O,V}
        return new{O,V,typeof(derivs)}(value, derivs)
    end
end

mutable struct MutableDiffResult{O,V,D<:Tuple} <: DiffResult{O,V,D}
    value::V
    derivs::D # ith element = ith-order derivative
    function MutableDiffResult(value::V, derivs::NTuple{O,Any}) where {O,V}
        return new{O,V,typeof(derivs)}(value, derivs)
    end
end

################
# Constructors #
################

"""
    DiffResult(value::Union{Number,AbstractArray}, derivs::Tuple{Vararg{Number}})
    DiffResult(value::Union{Number,AbstractArray}, derivs::Tuple{Vararg{AbstractArray}})

Return `r::DiffResult`, with output value storage provided by `value` and output derivative
storage provided by `derivs`.

In reality, `DiffResult` is an abstract supertype of two concrete types, `MutableDiffResult`
and `ImmutableDiffResult`. If all `value`/`derivs` are all `Number`s or `StaticArray`s,
then `r` will be immutable (i.e. `r::ImmutableDiffResult`). Otherwise, `r` will be mutable
(i.e. `r::MutableDiffResult`).

Note that `derivs` can be provide in splatted form, i.e. `DiffResult(value, derivs...)`.
"""
DiffResult

DiffResult(value::Number, derivs::Tuple{Vararg{Number}}) = ImmutableDiffResult(value, derivs)
DiffResult(value::Number, derivs::Tuple{Vararg{StaticArray}}) = ImmutableDiffResult(value, derivs)
DiffResult(value::StaticArray, derivs::Tuple{Vararg{StaticArray}}) = ImmutableDiffResult(value, derivs)
DiffResult(value::Number, derivs::Tuple{Vararg{AbstractArray}}) = MutableDiffResult(value, derivs)
DiffResult(value::AbstractArray, derivs::Tuple{Vararg{AbstractArray}}) = MutableDiffResult(value, derivs)
DiffResult(value::Union{Number,AbstractArray}, derivs::Union{Number,AbstractArray}...) = DiffResult(value, derivs)

"""
    GradientResult(x::AbstractArray)

Construct a `DiffResult` that can be used for gradient calculations where `x` is the input
to the target function.

Note that `GradientResult` allocates its own storage; `x` is only used for type and
shape information. If you want to allocate storage yourself, use the `DiffResult`
constructor instead.
"""
GradientResult(x::AbstractArray) = DiffResult(first(x), similar(x))
GradientResult(x::StaticArray) = DiffResult(first(x), x)

"""
    JacobianResult(x::AbstractArray)

Construct a `DiffResult` that can be used for Jacobian calculations where `x` is the input
to the target function. This method assumes that the target function's output dimension
equals its input dimension.

Note that `JacobianResult` allocates its own storage; `x` is only used for type and
shape information. If you want to allocate storage yourself, use the `DiffResult`
constructor instead.
"""
JacobianResult(x::AbstractArray) = DiffResult(similar(x), similar(x, length(x), length(x)))
JacobianResult(x::StaticArray) = DiffResult(x, zeros(StaticArrays.similar_type(typeof(x), Size(length(x),length(x)))))

"""
    JacobianResult(y::AbstractArray, x::AbstractArray)

Construct a `DiffResult` that can be used for Jacobian calculations where `x` is the
input to the target function, and `y` is the output (e.g. when taking the Jacobian
of `f!(y, x)`).

Like the single argument version, `y` and `x` are only used for type and
shape information and are not stored in the returned `DiffResult`.
"""
JacobianResult(y::AbstractArray, x::AbstractArray) = DiffResult(similar(y), similar(y, length(y), length(x)))
JacobianResult(y::StaticArray, x::StaticArray) = DiffResult(y, zeros(StaticArrays.similar_type(typeof(x), Size(length(y),length(x)))))

"""
    HessianResult(x::AbstractArray)

Construct a `DiffResult` that can be used for Hessian calculations where `x` is the
input to the target function.

Note that `HessianResult` allocates its own storage; `x` is only used for type and
shape information. If you want to allocate storage yourself, use the `DiffResult`
constructor instead.
"""
HessianResult(x::AbstractArray) = DiffResult(first(x), zeros(eltype(x), size(x)), similar(x, length(x), length(x)))
HessianResult(x::StaticArray) = DiffResult(first(x), x, zeros(StaticArrays.similar_type(typeof(x), Size(length(x),length(x)))))

#############
# Interface #
#############

@generated function tuple_eltype(x::Tuple, ::Type{Val{i}}) where {i}
    return quote
        $(Expr(:meta, :inline))
        return $(x.parameters[i])
    end
end

@generated function tuple_setindex(x::NTuple{N,Any}, y, ::Type{Val{i}}) where {N,i}
    T = x.parameters[i]
    new_tuple = Expr(:tuple, [ifelse(i == n, :(convert($T, y)), :(x[$n])) for n in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        return $new_tuple
    end
end

Base.eltype(r::DiffResult) = eltype(typeof(r))

Base.eltype(::Type{D}) where {O,V,D<:DiffResult{O,V}} = eltype(V)

Base.:(==)(a::DiffResult, b::DiffResult) = a.value == b.value && a.derivs == b.derivs

Base.copy(r::DiffResult) = DiffResult(copy(r.value), map(copy, r.derivs))

# value/value! #
#--------------#

"""
    value(r::DiffResult)

Return the primal value stored in `r`.

Note that this method returns a reference, not a copy.
"""
value(r::DiffResult) = r.value

"""
    value!(r::DiffResult, x)

Return `s::DiffResult` with the same data as `r`, except for `value(s) == x`.

This function may or may not mutate `r`. If `r::ImmutableDiffResult`, a totally new
instance will be created and returned, whereas if `r::MutableDiffResult`, then `r` will be
mutated in-place and returned. Thus, this function should be called as `r = value!(r, x)`.
"""
value!(r::MutableDiffResult, x::Number) = (r.value = x; return r)
value!(r::MutableDiffResult, x::AbstractArray) = (copyto!(value(r), x); return r)
value!(r::ImmutableDiffResult{O,V}, x::Union{Number,AbstractArray}) where {O,V} = ImmutableDiffResult(convert(V, x), r.derivs)

"""
    value!(f, r::DiffResult, x)

Equivalent to `value!(r::DiffResult, map(f, x))`, but without the implied temporary
allocation (when possible).
"""
value!(f, r::MutableDiffResult, x::Number) = (r.value = f(x); return r)
value!(f, r::MutableDiffResult, x::AbstractArray) = (map!(f, value(r), x); return r)
value!(f, r::ImmutableDiffResult{O,V}, x::Number) where {O,V} = value!(r, convert(V, f(x)))
value!(f, r::ImmutableDiffResult{O,V}, x::AbstractArray) where {O,V} = value!(r, convert(V, map(f, x)))

# derivative/derivative! #
#------------------------#

"""
    derivative(r::DiffResult, ::Type{Val{i}} = Val{1})

Return the `ith` derivative stored in `r`, defaulting to the first derivative.

Note that this method returns a reference, not a copy.
"""
derivative(r::DiffResult, ::Type{Val{i}} = Val{1}) where {i} = r.derivs[i]

"""
    derivative!(r::DiffResult, x, ::Type{Val{i}} = Val{1})

Return `s::DiffResult` with the same data as `r`, except `derivative(s, Val{i}) == x`.

This function may or may not mutate `r`. If `r::ImmutableDiffResult`, a totally new
instance will be created and returned, whereas if `r::MutableDiffResult`, then `r` will be
mutated in-place and returned. Thus, this function should be called as
`r = derivative!(r, x, Val{i})`.
"""
function derivative!(r::MutableDiffResult, x::Number, ::Type{Val{i}} = Val{1}) where {i}
    r.derivs = tuple_setindex(r.derivs, x, Val{i})
    return r
end

function derivative!(r::MutableDiffResult, x::AbstractArray, ::Type{Val{i}} = Val{1}) where {i}
    copyto!(derivative(r, Val{i}), x)
    return r
end

function derivative!(r::ImmutableDiffResult, x::Union{Number,StaticArray}, ::Type{Val{i}} = Val{1}) where {i}
    return ImmutableDiffResult(value(r), tuple_setindex(r.derivs, x, Val{i}))
end

function derivative!(r::ImmutableDiffResult, x::AbstractArray, ::Type{Val{i}} = Val{1}) where {i}
    T = tuple_eltype(r.derivs, Val{i})
    return ImmutableDiffResult(value(r), tuple_setindex(r.derivs, T(x), Val{i}))
end

"""
    derivative!(f, r::DiffResult, x, ::Type{Val{i}} = Val{1})

Equivalent to `derivative!(r::DiffResult, map(f, x), Val{i})`, but without the implied
temporary allocation (when possible).
"""
function derivative!(f, r::MutableDiffResult, x::Number, ::Type{Val{i}} = Val{1}) where {i}
    r.derivs = tuple_setindex(r.derivs, f(x), Val{i})
    return r
end

function derivative!(f, r::MutableDiffResult, x::AbstractArray, ::Type{Val{i}} = Val{1}) where {i}
    map!(f, derivative(r, Val{i}), x)
    return r
end

function derivative!(f, r::ImmutableDiffResult, x::Number, ::Type{Val{i}} = Val{1}) where {i}
    return derivative!(r, f(x), Val{i})
end

function derivative!(f, r::ImmutableDiffResult, x::StaticArray, ::Type{Val{i}} = Val{1}) where {i}
    return derivative!(r, map(f, x), Val{i})
end

function derivative!(f, r::ImmutableDiffResult, x::AbstractArray, ::Type{Val{i}} = Val{1}) where {i}
    T = tuple_eltype(r.derivs, Val{i})
    return derivative!(r, map(f, T(x)), Val{i})
end

# special-cased methods #
#-----------------------#

"""
    gradient(r::DiffResult)

Return the gradient stored in `r`.

Equivalent to `derivative(r, Val{1})`.
"""
gradient(r::DiffResult) = derivative(r)

"""
    gradient!(r::DiffResult, x)

Return `s::DiffResult` with the same data as `r`, except `gradient(s) == x`.

Equivalent to `derivative!(r, x, Val{1})`; see `derivative!` docs for aliasing behavior.
"""
gradient!(r::DiffResult, x) = derivative!(r, x)

"""
    gradient!(f, r::DiffResult, x)

Equivalent to `gradient!(r::DiffResult, map(f, x))`, but without the implied temporary
allocation (when possible).

Equivalent to `derivative!(f, r, x, Val{1})`; see `derivative!` docs for aliasing behavior.
"""
gradient!(f, r::DiffResult, x) = derivative!(f, r, x)

"""
    jacobian(r::DiffResult)

Return the Jacobian stored in `r`.

Equivalent to `derivative(r, Val{1})`.
"""
jacobian(r::DiffResult) = derivative(r)

"""
    jacobian!(r::DiffResult, x)

Return `s::DiffResult` with the same data as `r`, except `jacobian(s) == x`.

Equivalent to `derivative!(r, x, Val{1})`; see `derivative!` docs for aliasing behavior.
"""
jacobian!(r::DiffResult, x) = derivative!(r, x)

"""
    jacobian!(f, r::DiffResult, x)

Equivalent to `jacobian!(r::DiffResult, map(f, x))`, but without the implied temporary
allocation (when possible).

Equivalent to `derivative!(f, r, x, Val{1})`; see `derivative!` docs for aliasing behavior.
"""
jacobian!(f, r::DiffResult, x) = derivative!(f, r, x)

"""
    hessian(r::DiffResult)

Return the Hessian stored in `r`.

Equivalent to `derivative(r, Val{2})`.
"""
hessian(r::DiffResult) = derivative(r, Val{2})

"""
    hessian!(r::DiffResult, x)

Return `s::DiffResult` with the same data as `r`, except `hessian(s) == x`.

Equivalent to `derivative!(r, x, Val{2})`; see `derivative!` docs for aliasing behavior.
"""
hessian!(r::DiffResult, x) = derivative!(r, x, Val{2})

"""
    hessian!(f, r::DiffResult, x)

Equivalent to `hessian!(r::DiffResult, map(f, x))`, but without the implied temporary
allocation (when possible).

Equivalent to `derivative!(f, r, x, Val{2})`; see `derivative!` docs for aliasing behavior.
"""
hessian!(f, r::DiffResult, x) = derivative!(f, r, x, Val{2})

###################
# Pretty Printing #
###################

Base.show(io::IO, r::ImmutableDiffResult) = print(io, "ImmutableDiffResult($(r.value), $(r.derivs))")

Base.show(io::IO, r::MutableDiffResult) = print(io, "MutableDiffResult($(r.value), $(r.derivs))")

end # module


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Level Set Implementation First Stage

# Reference paper: https://gmd.copernicus.org/articles/4/591/2011/gmd-4-591-2011.pdf

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux, BenchmarkTools
using Plots #GalacticOptim
using PyPlot

@parameters t x y θ
@variables u(..) #u corresponds to ψ in reference paper, z is the terrain height
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dy'~y
@derivatives Dxx''~x
@derivatives Dyy''~y

#operators
z    = x^2 + y^2
gn   = (Dx(u(t,x,y,θ))^2 + Dy(u(t,x,y,θ))^2)^(1/2) #gradient norm
∇u   = [Dx(u(t,x,y,θ)), Dy(u(t,x,y,θ))]
#Δu   = Dxx(u(t,x,y,θ)) + Dyy(u(t,x,y,θ))
Dxz  = 2*x
Dyz  = 2*y
#∇z   = [Dx(z),Dy(z)]

#fuel dependent coefficients
n    = ∇u/gn #normal to the fire region
tanϕ = sum([Dxz,Dyz].*n)
#f = open("namelist.fire.txt")

#Fuel parameters given by file (convert to Imperial unit system!!!!)
windrf = [0.36, 0.36, 0.44,  0.55,  0.42,  0.44,  0.44, 0.36, 0.36, 0.36,  0.36,  0.43,  0.46, 1e-7]
fgi    = [0.166, 0.897, 1.076, 2.468, 0.785, 1.345, 1.092, 1.121, 0.780, 2.694, 2.582, 7.749, 13.024, 1.e-7]
fueldepthm = [0.305, 0.305, 0.762, 1.829, 0.61,  0.762, 0.762, 0.061, 0.061, 0.305, 0.305, 0.701, 0.914, 0.305]
savr = [3500., 2784., 1500., 1739., 1683., 1564., 1562., 1889., 2484., 1764., 1182., 1145., 1159., 3500.]
fuelmce = [0.12, 0.15, 0.25, 0.20, 0.20, 0.25, 0.40, 0.30, 0.25, 0.25, 0.15, 0.20, 0.25, 0.12]
fueldens = [32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.,32.] #from namelist.fire: "! 32 if solid, 19 if rotten"
st = [0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555]
se = [0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010, 0.010]
cmbcnst  = 17.433e+06
fuelmc_g = 0.08  #or 0.09
weight = [7.,  7.,  7., 180., 100., 100., 100., 900., 900., 900., 900., 900., 900., 7.]

a     = windrf[1]
zf    =
z0    =
w     = weight[1]
wl    = fgi[1]
δm    = fueldepthm[1]
sigma = savr[1]
Mx    = fuelmce[1]
ρP    = fueldens[1]
ST    = st[1]
SE    = se[1]
h     = cmbcnst[1]
Mf    = fuelmc_g[1]

#Fire spread rate equations
βop  = 0.1 #dummy
U    = 2 #wind vector (dummy)
w0   = wl/(1 + Mf)
ρb   = w0/δm #different from paper for units reasons
β    = ρb/ρP
ξ    = exp((0.792 + 0.618*sigma^0.5)*(β+0.1))
ηs   = 0.174*SE^-0.19 #probably mistaken in the paper
ηM   = 1 - 2.59*Mf/Mx + 5.11*(Mf/Mx)^2 - 3.52*(Mf/Mx)^3
wn   = w0/(1 + ST)
Γmax = sigma^(1.5)/(495 + 0.594*sigma^(1.5))
A    = 1/(4.77*sigma^(0.1) - 7.27)
Γ    = Γmax*(β/βop)^A*exp(A*(1 - β/βop))
ϵ    = exp(-138/sigma)
Qig  = 250*β + 1116*Mf
C    = 7.47*exp(-0.133*sigma^0.55)
Ua   = a*U #U is the wind vector
E    = 0.715*exp(-0.000359*sigma)
IR   = Γ*wn*h*ηM*ηs
R0   = IR*ξ/(ρb*ϵ*Qig)             #spread rate without wind
ϕw   = C*max(Ua^β, (β/βop)^E)      #wind factor
ϕS   = 5.275*β^-0.3*tanϕ^2       #slope factor (problem is we need to give tanϕ as a number)

S = 2 #R0*(1 + ϕw + ϕS) #fire spread rate [max(S0, R0 + c*min(e,max(0,U))^b + d*max(0,tanϕ)^2)]

#level set equation
eq = Dt(u(t,x,y,θ)) + S*gn ~ 0
#eq2 = z(x,y) - x^2 - y^2 ~ 0

#eqs = [eq1,eq2]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           y ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)] # Discretization

dx = 0.0125; dy = 0.0125; dt = 0.0125;

#ts,xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

# Initial and boundary conditions
bcs = [u(0,x,y,θ) ~ 0]

# LOW LEVEL
# Neural network
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

#strategy = GridTraining()
strategy = StochasticTraining()
discretization = PhysicsInformedNN([dt,dx,dy],chain,strategy=strategy)

indvars = [t,x,y]
depvars = [u]
dim = length(domains)

expr_pde_loss_function = build_loss_function(eq,indvars,depvars)
expr_bc_loss_functions = [build_loss_function(bc,indvars,depvars) for bc in bcs]
train_sets = generate_training_sets(domains,[dt,dx,dy],bcs,indvars,depvars)

train_domain_set, train_bound_set, train_set = train_sets

phi = discretization.phi
autodiff = discretization.autodiff
derivative = discretization.derivative
initθ = discretization.initθ

pde_loss_function = get_loss_function(eval(expr_pde_loss_function),
                                      train_domain_set,
                                      phi,
                                      derivative,
                                      strategy)
bc_loss_function = get_loss_function(eval.(expr_bc_loss_functions),
                                     train_bound_set,
                                     phi,
                                     derivative,
                                     strategy)

function loss_function(θ,p)
    return pde_loss_function(θ) + bc_loss_function(θ)
end

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

f = NeuralPDE.OptimizationFunction(loss_function, GalacticOptim.AutoZygote())

prob = GalacticOptim.OptimizationProblem(f, initθ)

# optimizer
opt = GalacticOptim.ADAM(0.01)
res = GalacticOptim.solve(prob, opt; cb = cb, maxiters = 30)
phi = discretization.phi

ts,xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

#u_predict_contourf = reshape([first(phi([t,x,y],res.minimizer)) for t in ts for x in xs for y in ys] ,length(ys),length(xs),length(ts))
#Plots.plot(ts, xs, ys, u_predict_contourf, linetype=:contourf,title = "predict")

#u_predict = [[first(phi([t,x,y],res.minimizer)) for x in xs for y in ys] for t in ts ]
u_predict = [reshape([first(phi([t,x,y],res.minimizer)) for x in xs for y in ys], (length(xs),length(ys))) for t in ts]
p1= Plots.plot(xs, u_predict[2],title = "t = 0.1");
p2= Plots.plot(xs, u_predict[6],title = "t = 0.5");
p3= Plots.plot(xs, u_predict[end],title = "t = 1");
Plots.plot(p1,p2,p3)

#HIGH LEVEL
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN([dt,dx,dy],
                                             chain,
                                             strategy = NeuralPDE.StochasticTraining(include_frac=0.9))
pde_system = PDESystem(eq,bcs,domains,[t,x,y],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.04), progress = false; cb = cb, maxiters=300)
phi = discretization.phi

xs,ys,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

analytic_sol_func(x,y,t) = []
u_real  = [[analytic_sol_func(t,x)[i] for t in ts for x in xs] for i in 1:3]
u_predict  = [[phi([x,y,t],res.minimizer)[i] for x in xs for y in ys for t in ts] for i in 1:4]
diff_u = [abs.(u_real[i] .- u_predict[i] ) for i in 1:3]

for i in 1:3
    p1 = plot(xs, ts, u_real[i], st=:surface,title = "u$i, analytic");
    p2 = plot(xs, ts, u_predict[i], st=:surface,title = "predict");
    p3 = plot(xs, ts, diff_u[i],linetype=:contourf,title = "error");
    plot(p1,p2,p3)
    savefig("sol_u$i")
end

#Another Trial
@parameters t x y
@variables u(..), X(..), Y(..), z(..) #u corresponds to ψ in reference paper, z is the terrain height
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dy'~y

#operators
∇u = [Dt(u(X(t),Y(t)))/Dt(X(t)),Dt(u(X(t),Y(t)))/Dt(Y(t))]
gn = ((Dt(u(X(t),Y(t)))/Dt(X(t)))^2 + (Dt(u(X(t),Y(t)))/Dt(X(t)))^2)^(1/2) #gradient norm
∇z = [Dx(z(x,y)),Dy(z(x,y))]

#fuel dependent coefficients
n    = ∇u/gn #normal to the fire region
S0   = 1
R0   = 3
b    = 2
c    = 1
d    = 1
e    = 1
U    = 4#wind vector
tanϕ = ∇z.*n

#definitions
S = 1#max(S0, R0 + c*min(e,max(0,U))^b + d*max(0,tanϕ)^2)

#level set equation
eq = Dt(u(X(t),Y(t))) + gn ~ 0

# Initial and boundary conditions

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           y ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)]# Discretization

dx = 0.125; dy= 0.125; dt = 0.125;

# Neural network
chain = FastChain(FastDense(7,32,Flux.σ),FastDense(32,32,Flux.σ),FastDense(32,1))

#strategy = GridTraining()
strategy = StochasticTraining()
discretization = PhysicsInformedNN([dt,dx,dy],chain,strategy=strategy)

indvars = [t,x,y]
depvars = [u]
dim = length(domains)

expr_pde_loss_function = build_loss_function(eq,indvars,depvars)
expr_bc_loss_functions = [build_loss_function(bc,indvars,depvars) for bc in bcs]
train_sets = generate_training_sets(domains,dx,bcs,indvars,depvars)

train_domain_set, train_bound_set, train_set = train_sets

phi = discretization.phi
autodiff = discretization.autodiff
derivative = discretization.derivative
initθ = discretization.initθ

pde_loss_function = get_loss_function(eval(expr_pde_loss_function),
                                      train_domain_set,
                                      phi,
                                      derivative,
                                      strategy)
bc_loss_function = get_loss_function(eval.(expr_bc_loss_functions),
                                     train_bound_set,
                                     phi,
                                     derivative,
                                     strategy)

function loss_function(θ,p)
    return pde_loss_function(θ) + bc_loss_function(θ)
end

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

f = NeuralPDE.OptimizationFunction(loss_function, GalacticOptim.AutoZygote())

prob = GalacticOptim.OptimizationProblem(f, initθ)

# optimizer
opt = GalacticOptim.ADAM(0.01)
res = GalacticOptim.solve(prob, opt; cb = cb, maxiters=30)
prob
phi = discretization.phi

ts,xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
#u_predict_contourf = reshape([first(phi([t,x,y],res.minimizer)) for t in ts for x in xs for y in ys] ,length(ys),length(xs),length(ts))
#Plots.plot(ts, xs, ys, u_predict_contourf, linetype=:contourf,title = "predict")

#u_predict = [[first(phi([t,x,y],res.minimizer)) for x in xs for y in ys] for t in ts ]
u_predict = [reshape([first(phi([t,x,y],res.minimizer)) for x in xs for y in ys], (length(xs),length(ys))) for t in ts ]
p1= Plots.plot(xs, u_predict[2],title = "t = 0.1");
p2= Plots.plot(xs, u_predict[6],title = "t = 0.5");
p3= Plots.plot(xs, u_predict[end],title = "t = 1");
Plots.plot(p1,p2,p3)


text:
import Pkg; Pkg.add("InfrastructureSystems")

using DelimitedFiles
using Documenter

# Include SLiDE modules.
using SLiDE

if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

DocMeta.setdocmeta!(SLiDE, :DocTestSetup, :(using SLiDE); recursive=true)

# Now, generate the documentation.
makedocs(clean = true,
    modules = [SLiDE],
    format = Documenter.HTML(
        mathengine = Documenter.MathJax(),
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    sitename = "SLiDE.jl",
    authors = "Jonathon Becker, Maxwell Brown, Caroline L. Hughes",
    # workdir = "../",
    pages = [
        "Home" => "index.md",
        "Introduction" => Any[
            "Data" => Any[
                "Overview" => "man/data/overview.md",
                "Preparation" => "man/data/preparation.md",
            ],
            "Build" => Any[
                "Overview" => "man/build/overview.md",
                "Partition" => "man/build/partition.md",
                "Share" => "man/build/share.md",
                "Disaggregate" => "man/build/disagg.md",
            ],
            # "EEM" => Any[
            #     "SEDS" => "man/eem/seds.md",
            # ],
            "Scaling" => "man/scaling.md",
            "Parameters" => "man/parameters.md"
        ],
        "API" => [
            "Types" => map(
                s -> "api/types/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/api/types")))),
            "Functions" => map(
                s -> "api/functions/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/api/functions")))),
            "Internals" => map(
                s -> "api/internals/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/api/internals"))))
            ],
        "Model" => "api/model.md",
    ]
)

# deploydocs(
#     repo = "https://github.com/NREL/SLiDE.git",
#     target = "build",
#     branch = "gh-pages",
#     devbranch = "docs",
#     devurl = "dev",
#     versions = ["stable" => "v^", "v#.#"],
# )

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using KhepriMeshCat
using Test

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


text:
"""
    mul!(C, A, B)
    mul!(C, A, B; kwargs...)

Multiply A×B and store the result in C (overwriting the contents of C).

This function uses multithreading.

This function is part of the public API of Gaius.
"""
function mul! end

"""
    mul(A, B)
    mul(A, B; kwargs...)

Multiply A×B and return the result.

This function uses multithreading.

This function is part of the public API of Gaius.
"""
function mul end

"""
    mul_serial(A, B)
    mul_serial(A, B; kwargs...)

Multiply A×B and return the result.

This function will run single-threaded.

This function is part of the public API of Gaius.
"""
function mul_serial end

"""
    mul_serial!(C, A, B)
    mul_serial!(C, A, B; kwargs...)

Multiply A×B and store the result in C (overwriting the contents of C).

This function will run single-threaded.

This function is part of the public API of Gaius.
"""
function mul_serial! end

function mul(A::MatTypes, B::MatTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Matrix{T}(undef, size(A,1), size(B,2))
    mul!(C, A, B)
    C
end

function mul_serial(A::MatTypes, B::MatTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Matrix{T}(undef, size(A,1), size(B,2))
    mul_serial!(C, A, B)
    C
end

function mul(A::StructArray{Complex{T}, 2}, B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Matrix{T}(undef, size(A, 1), size(B,2)),
                                 Matrix{T}(undef, size(A, 1), size(B,2))))
    mul!(C, A, B)
    C
end

function mul_serial(A::StructArray{Complex{T}, 2}, B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Matrix{T}(undef, size(A, 1), size(B,2)),
                                 Matrix{T}(undef, size(A, 1), size(B,2))))
    mul_serial!(C, A, B)
    C
end

function mul!(C::AbstractArray{T}, A::AbstractArray{T}, B::AbstractArray{T};
              block_size = nothing, sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C, A, B)

    _block_size = choose_block_size(C, A, B, block_size)

    _mul!(multithreaded, C, A, B, _block_size)
    C
end

function mul_serial!(C::AbstractArray{T}, A::AbstractArray{T}, B::AbstractArray{T};
                              block_size = nothing, sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C, A, B)

    _block_size = choose_block_size(C, A, B, block_size)

    _mul!(singlethreaded, C, A, B, _block_size)
    C
end

function mul!(C::StructArray{Complex{T}}, A::StructArray{Complex{T}}, B::StructArray{Complex{T}};
              block_size = default_block_size(), sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C.re, A.re, B.re)

    _block_size = choose_block_size(C, A, B, block_size)

    GC.@preserve C A B begin
        Cre, Cim = C.re, C.im
        Are, Aim = A.re, A.im
        Bre, Bim = B.re, B.im
        _mul!(    multithreaded,     Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(multithreaded,     Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    multithreaded,     Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(multithreaded,     Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    C
end

function mul_serial!(C::StructArray{Complex{T}}, A::StructArray{Complex{T}}, B::StructArray{Complex{T}};
                              block_size = default_block_size(), sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C.re, A.re, B.re)

    _block_size = choose_block_size(C, A, B, block_size)

    GC.@preserve C A B begin
        Cre, Cim = (C.re), (C.im)
        Are, Aim = (A.re), (A.im)
        Bre, Bim = (B.re), (B.im)
        _mul!(    singlethreaded,     Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(singlethreaded,     Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    singlethreaded,     Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(singlethreaded,     Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    C
end

function mul!(C::Adjoint{Complex{T}, <:StructArray{Complex{T}}},
              A::Adjoint{Complex{T}, <:StructArray{Complex{T}}},
              B::StructArray{Complex{T}};
              block_size = default_block_size(), sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C.parent.re', A.parent.re', B.re)

    _block_size = choose_block_size(C, A, B, block_size)
    A.parent.im .= (-).(A.parent.im) #ugly hack
    GC.@preserve C A B begin
        Cre, Cim = (C.parent.re'), (C.parent.im')
        Are, Aim = (A.parent.re'), (A.parent.im')
        Bre, Bim = (B.re), (B.im)
        _mul!(    multithreaded,     Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(multithreaded,     Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    multithreaded,     Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(multithreaded,     Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    A.parent.im .= (-).(A.parent.im)
    C.parent.im .= (-).(C.parent.im) # ugly hack
    C
end

function mul_serial!(C::Adjoint{Complex{T}, <:StructArray{Complex{T}}},
                              A::Adjoint{Complex{T}, <:StructArray{Complex{T}}},
                              B::StructArray{Complex{T}};
                              block_size = default_block_size(), sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C.parent.re', A.parent.re', B.re)

    _block_size = choose_block_size(C, A, B, block_size)
    A.parent.im .= (-).(A.parent.im) #ugly hack
    GC.@preserve C A B begin
        Cre, Cim = (C.parent.re'), (C.parent.im')
        Are, Aim = (A.parent.re'), (A.parent.im')
        Bre, Bim = (B.re), (B.im)
        _mul!(    singlethreaded,     Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(singlethreaded,     Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    singlethreaded,     Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(singlethreaded,     Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    A.parent.im .= (-).(A.parent.im)
    C.parent.im .= (-).(C.parent.im) # ugly hack
    C
end

function mul!(C::Transpose{Complex{T}, <:StructArray{Complex{T}}},
              A::Transpose{Complex{T}, <:StructArray{Complex{T}}},
              B::StructArray{Complex{T}};
              block_size = default_block_size(), sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C.parent.re |> transpose, A.parent.re |> transpose, B.re)

    _block_size = choose_block_size(C, A, B, block_size)

    GC.@preserve C A B begin
        Cre, Cim = (C.parent.re |> transpose), (C.parent.im |> transpose)
        Are, Aim = (A.parent.re |> transpose), (A.parent.im |> transpose)
        Bre, Bim = (B.re), (B.im)
        _mul!(    multithreaded,     Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(multithreaded,     Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    multithreaded,     Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(multithreaded,     Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    C
end

function mul_serial!(C::Transpose{Complex{T}, <:StructArray{Complex{T}}},
                              A::Transpose{Complex{T}, <:StructArray{Complex{T}}},
                              B::StructArray{Complex{T}};
                              block_size = default_block_size(), sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C.parent.re |> transpose, A.parent.re |> transpose, B.re)

    _block_size = choose_block_size(C, A, B, block_size)

    GC.@preserve C A B begin
        Cre, Cim = (C.parent.re |> transpose), (C.parent.im |> transpose)
        Are, Aim = (A.parent.re |> transpose), (A.parent.im |> transpose)
        Bre, Bim = (B.re), (B.im)
        _mul!(    singlethreaded,     Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(singlethreaded,     Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    singlethreaded,     Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(singlethreaded,     Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    C
end

function mul(A::MatTypes, B::VecTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Vector{T}(undef, size(A,1))
    mul!(C, A, B)
    C
end

function mul_serial(A::MatTypes, B::VecTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Vector{T}(undef, size(A,1))
    mul_serial!(C, A, B)
    C
end

function mul(A::StructArray{Complex{T}, 2}, B::StructArray{Complex{T}, 1}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Vector{T}(undef, size(A, 1)),
                                 Vector{T}(undef, size(A, 1))))
    mul!(C, A, B)
    C
end

function mul_serial(A::StructArray{Complex{T}, 2}, B::StructArray{Complex{T}, 1}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Vector{T}(undef, size(A, 1)),
                                 Vector{T}(undef, size(A, 1))))
    mul_serial!(C, A, B)
    C
end

function mul(A::CoVecTypes, B::MatTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Vector{T}(undef, size(B,2))'
    mul!(C, A, B)
    C
end

function mul_serial(A::CoVecTypes, B::MatTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Vector{T}(undef, size(B,2))'
    mul_serial!(C, A, B)
    C
end

function mul(A::Adjoint{Complex{T}, <:StructArray{Complex{T}, 1}},
             B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Vector{T}(undef, size(B, 2)),
                                 Vector{T}(undef, size(B, 2))))'
    mul!(C, A, B)
    C
end

function mul_serial(A::Adjoint{Complex{T}, <:StructArray{Complex{T}, 1}},
                             B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Vector{T}(undef, size(B, 2)),
                                 Vector{T}(undef, size(B, 2))))'
    mul_serial!(C, A, B)
    C
end

function mul(A::Transpose{Complex{T}, <:StructArray{Complex{T}, 1}},
             B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Vector{T}(undef, size(B, 2)),
                                 Vector{T}(undef, size(B, 2)))) |> transpose
    mul!(C, A, B)
    C
end

function mul_serial(A::Transpose{Complex{T}, <:StructArray{Complex{T}, 1}},
                             B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Vector{T}(undef, size(B, 2)),
                                 Vector{T}(undef, size(B, 2)))) |> transpose
    mul_serial!(C, A, B)
    C
end


text:
# Fallback: Collection of selections
select(t::CTable, sel) = isempty(sel) ? CTable(NamedTuple()) : mapreduce(x -> select(t, x), merge, sel)

#-----------------------------------------------------------------------------# field name/index
field_name(t::CTable, sel) = error("Selection $sel returns multiple field names.")
field_name(t::CTable, sel::Symbol) = sel
field_name(t::CTable, sel::String) = Symbol(sel)
field_name(t::CTable, sel::Integer) = names(t)[sel]

field_index(t::CTable, sel) = error("Selection $sel returns multiple field indexes.")
field_index(t::CTable, sel::Symbol) = findfirst(==(sel), names(t))
field_index(t::CTable, sel::String) = field_index(t, Symbol(sel))
field_index(t::CTable, sel::Integer) = sel

# selections that return a single column
select(t::CTable, sel::Union{Symbol, String, Integer}) = CTable(select(cols(t), [field_name(t, sel)]))

select(t::CTable, sel::Type) = select(t, names(t)[collect(map(x -> eltype(x) <: sel, cols(t)))])

select(t::CTable, sel::Regex) = select(t, filter(x -> occursin(sel, string(x)), names(t)))

#-----------------------------------------------------------------------------# Col
select(t::CTable, c::Col) = cols(t)[field_name(t, c.sel)]

#-----------------------------------------------------------------------------# Not
select(t::CTable, not::Not) = CTable(delete(cols(t), names(select(t, not.sel))))

#-----------------------------------------------------------------------------# Before
field_index(t::CTable, b::Before) = field_index(t, b.sel) - 1
function select(t::CTable, b::Before)
    i = field_index(t, b)
    i < 1 ? error("Attempted to select 0-th field.") : select(t, i)
end

#-----------------------------------------------------------------------------# After
field_index(t::CTable, a::After) = field_index(t, a.sel) + 1
function select(t::CTable, a::After)
    i = field_index(t, a)
    i > length(ncols(t)) ?
        error("Attempted to select $i-th field when only $(ncols(t)) fields exist.") :
        select(t, i)
end

#-----------------------------------------------------------------------------# FromTo
select(t::CTable, b::FromTo) = select(t, field_index(t, b.first):field_index(t, b.last))

#-----------------------------------------------------------------------------# And
And(args...) = And(args)
select(t::CTable, a::And) = mapreduce(x -> select(t, x), (a,b) -> select(a, intersect(names(a), names(b))), a.sel)


text:
export recursive_ols, recursive_ols!

recursive_ols_overview =
"""Calculates the OLS estimators for the model y ~ x recursively and returns
them as an k by n - R matrix. The jth column of `estimates` is the OLS
estimator estimated by regressing y[1:(R+j-1)] on x[1:(R+j-1),:] and is used to
forecast y[R+j] rom x[R+j,:].

* `recursive_ols` creates a new matrix containing the estimates and does basic
  checks that the arrays have the right dimensions. The actual calculations are
  done by `recursive_ols!`

* `recursive_ols!` fills its first argument with the OLS estimates and does not
  check the dimensions of its arguments. It is intended to be used inside
  simulations, etc., where the arguments are known to be safe.

There are also a methods that return the OLS estimators and the implied
prediction errors.
"""
Index = Union(AbstractVector{Int},Range{Int},UnitRange{Int})

@doc recursive_ols_overview ->
function recursive_ols!(estimates::AbstractMatrix{Float64},
                        y::AbstractVector{Float64},
                        x::AbstractMatrix{Float64},
                        order::Index)
    k, P = size(estimates)
    R = length(y) - P
    ## QR decomposition used to construct the (hypothetical) period 0
    ## estimate
    Rmat = vcat(qrfact(hcat(x[order[1:R-1],:], y[order[1:R-1]]))[:R], zeros(1, k+1))
    for j in 1:P
        ols_update!(sub(estimates,:,j), Rmat, x, y, order[j+R-1])
    end
end

function recursive_ols!(estimates::AbstractMatrix{Float64},
                        errors::AbstractVector{Float64},
                        y::AbstractVector,
                        x::AbstractMatrix,
                        order::Index)
    k, P = size(estimates)
    R = length(y) - P
    ## QR decomposition used to construct the (hypothetical) period 0
    ## estimate
    Rmat = vcat(qrfact(hcat(x[order[1:R-1],:], y[order[1:R-1]]))[:R], zeros(1, k+1))
    for j in 1:P
        ols_update!(sub(estimates,:,j), Rmat, x, y, order[j+R-1])
        errors[j] = y[order[R+j]] - scalar(x[order[R+j],:] * estimates[:,j])
    end
end

@doc recursive_ols_overview ->
function recursive_ols(y::AbstractVector, x::AbstractMatrix, R::Integer, order::Index)
    length(y) == size(x,1) ||
      error("`y` should be a vector with the same number of elements as rows of `x`")
    estimates = Array(Float64, size(x, 2), length(y) - R)
    recursive_ols!(estimates, y, x, order)
    return estimates
end

@doc "Function that actually updates the OLS coefficients." ->
function ols_update!(b::AbstractVector{Float64}, R::AbstractMatrix{Float64},
                     x::AbstractMatrix{Float64}, y::AbstractVector{Float64},
                     j::Integer)
    k = size(x, 2)
    R[k+2,1:k] = x[j,:]
    R[k+2,k+1] = y[j]
    for i in 1:k
        Gi, _ = givens(R, i, k+2, i)
        A_mul_B!(Gi, R)
    end
    b[1:k] = copy(R[1:k,k+1])
    Base.LinAlg.LAPACK.trtrs!('U', 'N', 'N', sub(R, 1:k, 1:k), b)
end

## Additional convenience methods
recursive_ols!(estimates, y, x) = recursive_ols!(estimates, y, x, 1:length(y))
recursive_ols!(estimates, errors, y, x) =
    recursive_ols!(estimates, errors, y, x, 1:length(y))
recursive_ols(y, x, R) = recursive_ols(y, x, R, 1:length(y))


text:
using SafeTestsets

@safetestset "Networks Types" begin
    include("networks.jl")
end
@safetestset "Neuron Types" begin
    include("learning.jl")
end
#TODO:
# @safetestset "Neuron Types" begin include("neuron_type.jl") end


text:
export MultiContrastFile
import Base: getindex, length, iterate

struct MultiContrastFile <: MPIFile
  files::Vector{MPIFile}

  function MultiContrastFile(filenames::Vector{String})
    return new([MPIFile(f) for f in filenames])
  end

  function MultiContrastFile(f::Vector{T}) where T<:MPIFile
    return new(f)
  end
end

getindex(f::MultiContrastFile, index::Integer) = f.files[index]
length(f::MultiContrastFile) = length(f.files)
start_(f::MultiContrastFile) = 1
next_(f::MultiContrastFile, state) = (f[state],state+1)
done_(f::MultiContrastFile, state) = state > length(f.files)
iterate(f::MultiContrastFile, s=start_(f)) = done_(f, s) ? nothing : next_(f, s)

function Base.show(io::IO, f::MultiContrastFile)
  print(io, "Multi Contrast File: ", f.files)
end

for op in [:filepath, :version, :uuid, :time, :studyName, :studyNumber, :studyUuid, :studyDescription,
            :experimentName, :experimentNumber, :experimentUuid, :experimentDescription,
            :experimentSubject, :experimentHasMeasurement,
            :experimentIsSimulation, :experimentIsCalibration, :experimentHasProcessing,
            :tracerName, :tracerBatch, :tracerVendor, :tracerVolume, :tracerConcentration,
            :tracerSolute, :tracerInjectionTime,
            :scannerFacility, :scannerOperator, :scannerManufacturer, :scannerName,
            :scannerTopology, :acqNumBGFrames, :acqGradient,
            :acqStartTime, :acqNumFrames, :acqNumPeriodsPerFrame,
            :dfNumChannels, :dfBaseFrequency, :dfDivider, :dfStrength, :dfPhase,
            :dfCycle, :dfWaveform, :rxNumChannels, :acqNumAverages, :rxBandwidth,
            :rxNumSamplingPoints, :rxTransferFunction, :rxHasTransferFunction, :rxTransferFunctionFileName, :rxInductionFactor, :rxUnit,
            :rxDataConversionFactor,
            :calibSNR, :calibFov, :calibFovCenter, #:calibSize,
            :calibOrder, :calibPositions, :calibOffsetField,
            :calibDeltaSampleSize, :calibMethod]
  @eval MPIFiles.$op(f::MultiContrastFile) = MPIFiles.$op(f.files[1])
end

measPath(bs::MultiContrastFile) = [measPath(b) for b in bs]
gridSize(bs::MultiContrastFile) = [gridSize(b) for b in bs]
calibSize(bs::MultiContrastFile) = [calibSize(b) for b in bs]
gridSizeCommon(bs::MultiContrastFile) = gridSize(bs[1])
fov(b::MultiContrastFile) = fov(b[1])

function filterFrequencies(bSFs::MultiContrastFile; kargs...) where {T<:MPIFile}
  return intersect([filterFrequencies(bSF; kargs...) for bSF in bSFs]...)
end

### deprecated ###
measPath(bs::Vector{T}) where {T<:MPIFile} = [measPath(b) for b in bs]
gridSize(bs::Vector{T}) where {T<:MPIFile} = [gridSize(b) for b in bs]
calibSize(bs::Vector{T}) where {T<:MPIFile} = [calibSize(b) for b in bs]
gridSizeCommon(bs::Vector{T}) where {T<:MPIFile} = gridSize(bs[1])
fov(b::Vector{T}) where {T<:MPIFile} = fov(b[1])
calibFov(b::Vector{T}) where {T<:MPIFile} = calibFov(b[1])

function filterFrequencies(bSFs::Vector{T}; kargs...) where {T<:MPIFile}
  return intersect([filterFrequencies(bSF; kargs...) for bSF in bSFs]...)
end


text:
"""
    TestFunctionRunner.run(module::Module)
    TestFunctionRunner.run(modules::AbstractVector)

Run tests defined in a `module`.

This function discovers "test functions" in `module` and its "test modules"
recursively.  Test functions are nullary functions named `test` or nullary
functions whose name start with `test_`.  Test modules are the sub-modules
whose name starts with "Test" followed by a capital letter.
"""
TestFunctionRunner.run

function _run_options(;
    prepare_distributed::Bool = true,
    timeout::Union{Nothing,Real} = default_global_timeout(),
    failfast::Bool = default_fastfail(),
    kwargs...,
)
    runtests_options = (; timeout = timeout, failfast = failfast, kwargs...)
    return (prepare_distributed, failfast, runtests_options)
end

function TestFunctionRunner.run(m::Module; kwargs...)
    (prepare_distributed, failfast, runtests_options) = _run_options(; kwargs...)
    @debug "`TestFunctionRunner.run`" options = sort!(collect(kwargs))
    with_project(m) do
        if prepare_distributed
            load_everywhere(m)
        else
            prepare_local(m)
        end
        testset = testset_iter([m]; failfast = failfast) do m
            "$(nameof(m))" => function ()
                runtests(m; runtests_options...)
            end
        end
        Success(testset)
    end
end

function TestFunctionRunner.run(modules::AbstractVector; kwargs...)
    for m in modules
        TestFunctionRunner.run(m::Module; kwargs...)
    end
end

function default_global_timeout()
    timeout = get(ENV, "TEST_FUNCTION_RUNNER_JL_TIMEOUT", nothing)
    timeout === nothing && return timeout
    return parse(Float64, timeout)
end

function default_fastfail()
    str = get(ENV, "TEST_FUNCTION_RUNNER_JL_FASTFAIL", nothing)
    str === nothing && return false
    return lowercase(str) in ("true", "yes", "1")
end

struct Success
    testset::Any
end

function Base.show(io::IO, ::MIME"text/plain", ::Success)
    print(io, "Success: 🎉 All test passed")
    # TODO: print more info
end

function find_project_file(dir)
    p = joinpath(dir, "JuliaProject.toml")
    isfile(p) && return p
    p = joinpath(dir, "Project.toml")
    isfile(p) && return p
    return nothing
end

function root_module(m::Module)
    while true
        p = parentmodule(m)
        m === p && return m
        m = p
    end
end

function projectfile(m::Module)
    dir = dirname(dirname(pathof(root_module(m))))
    p = find_project_file(dir)
    p === nothing && error("cannot find the project for module $m")
    return p
end

function is_in_path(m::Module)
    project = projectfile(m)
    paths = Base.load_path()
    project in paths && return true
    realproject = realpath(project)
    realproject in paths && return true
    matches(path) = path == project || path == realproject
    return any(paths) do path
        matches(path) || matches(realpath(path))
    end
end

function with_project(f, m::Module)
    is_in_path(m) && return f()
    load_path = copy(LOAD_PATH)
    push!(LOAD_PATH, projectfile(m))
    try
        f()
    finally
        append!(empty!(LOAD_PATH), load_path)
    end
end

function load_everywhere(m::Module)
    pkgid = Base.PkgId(m)
    prepare_local(m)
    nprocs() == 1 && return
    @everywhere append!(empty!(LOAD_PATH), $(copy(LOAD_PATH)))
    @everywhere Base.require($pkgid)
    @everywhere $prepare_local($m)
end

function prepare_local(m::Module)
    LoadAllPackages.loadall(projectfile(m))
end

function should_test(m::Module)
    f = try
        m.should_test_module
    catch
        return true
    end
    return f()::Bool
end

function run_before_test_module_hook(m::Module)
    f = try
        m.before_test_module
    catch
        return
    end
    f()
    return
end

function run_after_test_module_hook(m::Module)
    f = try
        m.after_test_module
    catch
        return
    end
    f()
    return
end

function with_module_context(f, m::Module)
    ctx = try
        m.module_context
    catch
        nothing
    end
    if ctx === nothing
        return f()
    else
        return ctx(f)
    end
end

default_timeout_of(@nospecialize(_)) = nothing

function get_timeout_of(m::Module)
    try
        return m.timeout_of
    catch
        return default_timeout_of
    end
end

function testset_iter(f, iter; failfast::Bool)
    tests = Iterators.map(f, iter)
    if failfast
        for (_label, test) in tests
            test()
        end
    else
        @testset "$label" for (label, test) in tests
            test()
        end
    end
end

function runtests(
    m::Module;
    recursive::Bool = true,
    timeout::Union{Nothing,Real} = nothing,
    failfast::Bool = false,
)
    should_test(m) || return
    timeout === nothing ||
        timeout > 0 ||
        error("`timeout` should be positive; got: $timeout")
    timeout_of = get_timeout_of(m)
    run_before_test_module_hook(m)
    try
        with_module_context(m) do
            @debug "Testing module: `$m`"
            testset_iter(test_functions(m); failfast = failfast) do f
                "$f" => function ()
                    label = "$m.$f"
                    @debug "Testing function: `$label`"
                    tout = something(timeout_of(f), timeout, Some(nothing))
                    if tout === nothing
                        f()
                    else
                        Terminators.withtimeout(f, tout; label = label)
                    end
                end
            end
            recursive || return
            testset_iter(test_modules(m); failfast = failfast) do sub
                "$(nameof(sub))" => function ()
                    runtests(sub; recursive = recursive, timeout = timeout, failfast = failfast)
                end
            end
        end
    finally
        run_after_test_module_hook(m)
    end
end

is_test_function(name::Symbol) = name === :test || startswith(string(name), "test_")
is_test_module(name::Symbol) = match(r"^Test[A-Z]", string(name)) !== nothing

function test_functions(m::Module)
    tests = map(names(m, all = true)) do n
        is_test_function(n) || return nothing
        f = getproperty(m, n)
        f !== m || return nothing
        if f isa Type
            t = f
        else
            t = typeof(f)
        end
        parentmodule(t) === m || return nothing
        applicable(f) || return nothing  # removed by Revise?
        return f
    end
    filter!(!isnothing, tests)
    sort!(tests, by = string ∘ nameof)
    return tests
end

function test_modules(root::Module)
    modules = Module[]
    for n in names(root, all = true)
        m = getproperty(root, n)
        m isa Module || continue
        m === root && continue
        parentmodule(m) === root || continue
        is_test_module(nameof(m)) || continue
        push!(modules, m)
    end
    sort!(modules, by = string ∘ nameof)
    return modules
end

"""
    TestFunctionRunner.@run()

Run tests defined in the package `@__DIR__/\$TestPackage/src/\$TestPackage.jl`.

# Keyword Arguments

- `paths::AbstractVector{<:AbstractString}`: Additional load paths. Relative
  paths are resolved with respect to the parent direcotry of the current file.
  Example: `TestFunctionRunner.@run(paths = ["../benchmark/MyBenchmarks/"])`.

- `packages::AbstractVector{<:AbstractString}`: Specify packages to be tested by
  relative paths.
"""
macro run(options...)
    kwargs = map(assignment_as_kw, options)
    esc(:($at_run_impl($(QuoteNode(__source__)), $(QuoteNode(__module__)); $(kwargs...))))
end

function assignment_as_kw(ex)
    if Meta.isexpr(ex, :(=), 2)
        return Expr(:kw, ex.args...)
    end
    error("a key-value pair (`k = v`) is expected; got: $ex")
end

function at_run_impl(
    __source__::LineNumberNode,
    __module__::Module;
    paths::AbstractVector{<:AbstractString} = String[],
    packages::Union{AbstractVector{<:AbstractString},Nothing} = nothing,
    prepare_distributed::Bool = true,
    kwargs...,
)
    testdir = dirname(string(__source__.file))
    paths = joinpath.(testdir, String.(paths))
    if packages === nothing
        pkgid, project = find_test_package(testdir)
        pushfirst!(paths, project)
        m = load_test_package(pkgid, paths, prepare_distributed)
    else
        m = map(packages) do pkgpath
            local project = joinpath(testdir, pkgpath)
            local pkgid = pkgid_from_project_path(project)
            load_test_package(pkgid, [project; paths], prepare_distributed)
        end
    end
    Base.invokelatest(
        TestFunctionRunner.run,
        m;
        prepare_distributed = prepare_distributed,
        kwargs...,
    )
end

function load_test_package(
    pkgid::PkgId,
    paths::Vector{String},
    prepare_distributed::Bool,
)::Module
    missing_paths = setdiff!(map(abspath, paths), LOAD_PATH)
    try
        return Base.require(pkgid)
    catch
        isempty(missing_paths) && rethrow()
        append!(LOAD_PATH, missing_paths)
    end
    if prepare_distributed
        @everywhere append!(empty!(LOAD_PATH), $(copy(LOAD_PATH)))
    end
    @debug "Try loading `$pkgid` again after `LOAD_PATH` hack..." LOAD_PATH
    return Base.require(pkgid)
end

function find_test_package(dir::AbstractString)
    projects = Tuple{Base.PkgId,String}[]
    for n in readdir(dir)
        subdir = joinpath(dir, n)
        isdir(subdir) || continue

        path = find_project_file(subdir)
        path === nothing && continue

        toml = TOML.parsefile(path)
        name = get(toml, "name", nothing)
        uuid = get(toml, "uuid", nothing)
        name === nothing && continue
        uuid === nothing && continue
        isfile(subdir, "src", name * ".jl") || continue

        pkgid = PkgId(Base.UUID(uuid), name)
        push!(projects, (pkgid, path))
    end
    if isempty(projects)
        error("test package not found at: $dir")
    elseif length(projects) > 1
        error("multiple packages found at: $dir")
    end
    return projects[1]
end

function pkgid_from_project_path(path)
    if isdir(path)
        project = find_project_file(path)
    else
        project = path
    end
    toml = TOML.parsefile(project)
    name = toml["name"]::String
    uuid = toml["uuid"]::String
    return PkgId(Base.UUID(uuid), name)
end

#=
abstract type Testable end

struct TestModule <: Testable
    m::Module
    tests::Vector{Testable}
end

struct TestFunction <: Testable
    f::Any
end

testname(t::TestModule) = nameof(t.m)
testname(t::TestFunction) = nameof(t.f)

function TestFunctionRunner.run(testable::TestModule)
    @testset "$(testname(t))" for t in testable.tests
        TestFunctionRunner.run(t)
    end
end

TestFunctionRunner.run(testable::TestFunction) = testable.f()

TestFunctionRunner.collect(m::Module) = TestModule(m, discover(m))

function discover(m::Module)
    for n in names(root, all = true)
        m = getproperty(root, n)
        m isa Module || continue
        m === root && continue
        startswith(string(nameof(m)), "Test") || continue
        push!(modules, m)
    end
end
=#


text:
function mixing_length(param_set, ml_model::MinDisspLen{FT}) where {FT}
    c_m::FT = CPEDMF.c_m(param_set)
    c_d::FT = CPEDMF.c_d(param_set)
    smin_ub::FT = CPEDMF.smin_ub(param_set)
    smin_rm::FT = CPEDMF.smin_rm(param_set)
    l_max::FT = ICP.l_max(param_set)
    c_b::FT = ICP.static_stab_coeff(param_set)
    g::FT = CPP.grav(param_set)
    molmass_ratio::FT = CPP.molmass_ratio(param_set)
    vkc::FT = CPSGS.von_karman_const(param_set)
    ustar = ml_model.ustar
    z = ml_model.z
    tke_surf = ml_model.tke_surf
    ∂b∂z = ml_model.∇b.∂b∂z
    tke = ml_model.tke

    # kz scale (surface layer)
    if ml_model.obukhov_length < 0.0 #unstable
        l_W =
            vkc * z / (sqrt(tke_surf / ustar / ustar) * c_m) *
            min((1 - 100.0 * z / ml_model.obukhov_length)^0.2, 1 / vkc)
    else # neutral or stable
        l_W = vkc * z / (sqrt(tke_surf / ustar / ustar) * c_m)
    end

    # compute l_TKE - the production-dissipation balanced length scale
    a_pd = c_m * (ml_model.Shear² - ∂b∂z / ml_model.Pr) * sqrt(tke)
    # Dissipation term
    c_neg = c_d * tke * sqrt(tke)
    # Subdomain exchange term
    b_exch = ml_model.b_exch

    if abs(a_pd) > eps(FT) && 4 * a_pd * c_neg > -b_exch * b_exch
        l_TKE = max(-b_exch / 2 / a_pd + sqrt(b_exch * b_exch + 4 * a_pd * c_neg) / 2 / a_pd, 0)
    elseif abs(a_pd) < eps(FT) && abs(b_exch) > eps(FT)
        l_TKE = c_neg / b_exch
    else
        l_TKE = FT(0)
    end

    # compute l_N - the effective static stability length scale.
    N_eff = sqrt(max(∂b∂z, 0))
    if N_eff > 0.0
        l_N = min(sqrt(max(c_b * tke, 0)) / N_eff, l_max)
    else
        l_N = l_max
    end

    # add limiters
    l = SA.SVector(
        (l_N < eps(FT) || l_N > l_max) ? l_max : l_N,
        (l_TKE < eps(FT) || l_TKE > l_max) ? l_max : l_TKE,
        (l_W < eps(FT) || l_W > l_max) ? l_max : l_W,
    )

    # get soft minimum
    min_len, min_len_ind = findmin(l)
    mix_len = lamb_smooth_minimum(l, smin_ub, smin_rm)
    ml_ratio = mix_len / min_len
    return MixLen(min_len_ind, mix_len, ml_ratio)
end