text:
"""
    plot(arg1::Array; kwargs...)

reads (x,y) pairs from files [or standard input] and generates PostScript code that will plot lines,
polygons, or symbols at those locations on a map.

Full option list at [`psxy`]($(GMTdoc)plot.html)

Parameters
----------

- **A** | **steps** | **straight_lines** :: [Type => Str]

	By default, geographic line segments are drawn as great circle arcs.
	To draw them as straight lines, use this option.
    ($(GMTdoc)plot.html#a)
- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** | **shift** | **offset** :: [Type => Str]

    Offset the plot symbol or line locations by the given amounts dx/dy in cm, inch or points.
    ($(GMTdoc)plot.html#d)
- **E** | **error** | **error_bars** :: [Type => Str]

    Draw symmetrical error bars.
    ($(GMTdoc)plot.html#e)
- **F** | **conn** | **connection** :: [Type => Str]

    Alter the way points are connected
    ($(GMTdoc)plot.html#f)
- **G** | **fill** | **markerfacecolor** | **MarkerFaceColor** | **markercolor** | **mc** :: [Type => Str]

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    ($(GMTdoc)plot.html#g)
- **I** | **intens** :: [Type => Str | number]

    Use the supplied intens value (in the [-1 1] range) to modulate the fill color by simulating illumination.
    ($(GMTdoc)plot.html#i)
- **L** | **close** | **polygon** :: [Type => Str]

    Force closed polygons.
    ($(GMTdoc)plot.html#l)
- **N** | **no_clip** | **noclip** :: [Type => Str or []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)plot.html#n)
- $(GMT.opt_P)
- **S** | **symbol** | **marker** | **Marker** :: [Type => Str]

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    ($(GMTdoc)plot.html#s)
    Alternatively select a sub-set of symbols using the aliases: **marker** or **Marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
    + **y**, **y_dash**

    and select their sizes with the **markersize** or **size** keyword [default is 7p].
    The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
    points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT="c")*

- **W** | **pen** | **markeredgecolor** | **mec** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **pen** for the polygons and
    **markeredgecolor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)

- **Z** | **level** :: [Type => Str | NamedTuple]	`Arg = value|file[+f|+l] | (data=Array|Number, outline=_, fill=_)`

    Paint polygons after the level given as a cte or a vector with same size of number of polygons. Needs a color map.
    ($(GMTdoc)plot.html#z)

- **axis** | **aspect** :: [Type => Str]

    When equal to "equal" makes a square plot.
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_w)
- $(GMT.opt_swap_xy)
"""
function plot(arg1; first=true, kw...)
	common_plot_xyz("", cat_1_arg(arg1), "plot", first, false, kw...)
end
plot!(arg1; kw...) = plot(arg1; first=false, kw...)

function plot(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	common_plot_xyz("", cat_2_arg2(rang, [f(x) for x in rang]), "plot", first, false, kw...)
end
plot!(f::Function, rang=nothing; kw...) = plot(f, rang; first=false, kw...)

function plot(f1::Function, f2::Function, range_t=nothing; first=true, kw...)	# Parametric version
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; is3D=false, kw...), "plot", first, false, kw...)
end
plot!(f1::Function, f2::Function, range_t=nothing; kw...) = plot(f1, f2, range_t; first=false, kw...)

plot(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, arg1, "plot", true, false, kw...)
plot!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "plot", false, false, kw...)
plot(arg1, arg2; kw...)  = common_plot_xyz("", cat_2_arg2(arg1, arg2), "plot", true, false, kw...)
plot!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2), "plot", false, false, kw...)
plot(arg1::Number, arg2::Number; kw...)  = common_plot_xyz("", cat_2_arg2([arg1], [arg2]), "plot", true, false, kw...)
plot!(arg1::Number, arg2::Number; kw...) = common_plot_xyz("", cat_2_arg2([arg1], [arg2]), "plot", false, false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
plotyy(arg1, arg2; kwargs...)

Example:
```julia
plotyy([1 1; 2 2], [1.5 1.5; 3 3], R="0.8/3/0/5", title="Ai", ylabel=:Bla, xlabel=:Ble, seclabel=:Bli, show=1)
```
"""
function plotyy(arg1, arg2; first=true, kw...)
	d = KW(kw)
	(haskey(d, :xlabel)) ? (xlabel = string(d[:xlabel]);	delete!(d, :xlabel)) : xlabel = ""	# Only to used at the end
	(haskey(d, :seclabel)) ? (seclabel = string(d[:seclabel]);	delete!(d, :seclabel)) : seclabel = ""
	((val = find_in_dict(d, [:fmt])[1]) !== nothing) ? (fmt = arg2str(val)) : fmt = FMT[1]
	((val = find_in_dict(d, [:savefig :name])[1]) !== nothing) ? (savefig = arg2str(val)) : savefig = nothing
	Vd = ((val = find_in_dict(d, [:Vd])[1]) !== nothing) ? val : 0

	cmd, opt_B = parse_B(d, "", " -Baf -BW")
	if (opt_B != " -Baf -BW")
		if (occursin(" -Bx", opt_B) || occursin(" -By", opt_B) || occursin("+t", opt_B))
			# OK, so here's the problem. Both title and label maybe multi-words, case in which they will have the
			# form: "aa bb" but since they wont be fully parsed again by parse_B we loose the info about the sapces.
			# The trick is to find and replace the sapces inside the " " chunks like parse_B does and let it undo
			# the replacing without inserting extras -B, like the "aa -Bbb" that would otherwise result.
			t = findall(isequal('"'), opt_B)	# Find indices of the " and hope they came in pairs.
			if (!isempty(t))					# Have a title and maybe an ylabel (or vice-versa)
				opt_B = opt_B[1:t[1]] * replace(opt_B[t[1]+1:t[2]-1], ' '=>'\x7f') * opt_B[t[2]:end]
				if (length(t) == 4)				# Have also an ylabel
					opt_B = opt_B[1:t[3]] * replace(opt_B[t[3]+1:t[4]-1], ' '=>'\x7f') * opt_B[t[4]:end]
				end
			end
			d[:B] = replace(opt_B, "-B" => "")
		else
			d[:B] = " af W"
		end
	else
		d[:B] = " af W"
	end
	(Vd != 0) && (d[:Vd] = Vd)
	d[:lc] = "#0072BD"
	(haskey(d, :show)) ? (delete!(d, :show);  do_show = true) : do_show = false
	d[:par] = (MAP_FRAME_PEN="#0072BD", MAP_TICK_PEN="#0072BD", FONT_ANNOT_PRIMARY="#0072BD", FONT_LABEL="#0072BD")
	r1 = common_plot_xyz("", cat_1_arg(arg1), "plotyy", first, false, d...)

	(Vd != 0) && (d[:Vd] = Vd)
	(seclabel != "" && occursin(" ", seclabel)) && (seclabel = "\"" * seclabel * "\"")
	(seclabel != "") && (seclabel = " y+l" * seclabel)
	d[:B] = " af E"	* seclabel		# Also remember that previous -B was consumed in first call
	d[:lc]  = "#D95319"
	d[:par] = (MAP_FRAME_PEN="#D95319", MAP_TICK_PEN="#D95319", FONT_ANNOT_PRIMARY="#D95319", FONT_LABEL="#D95319")
	r2 = common_plot_xyz("", cat_1_arg(arg2), "plotyy", false, false, d...)

	(xlabel != "" && occursin(" ", xlabel)) && (xlabel = "\"" * xlabel * "\"")
	opt_B = (xlabel != "") ? "af Sn x+l" * xlabel : "af Sn"
	r3 = basemap!(J="", R="", B=opt_B, Vd=Vd, fmt=fmt, name=savefig, show=do_show)
	return (Vd == 2) ? [r1;r2;r3] : nothing
end
# ------------------------------------------------------------------------------------------------------

"""
plot3d(arg1::Array; kwargs...)

reads (x,y,z) triplets and generates PostScript code that will plot lines,
polygons, or symbols at those locations in 3-D.

Full option list at [`plot3d`]($(GMTdoc)plot3d.html)

Parameters
----------

- **A** | **steps** | **straight_lines** :: [Type => Str]

    By default, geographic line segments are drawn as great circle arcs. To draw them as straight lines, use this option.
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **C** | **color** :: [Type => Str]

    Give a CPT or specify -Ccolor1,color2[,color3,...] to build a linear continuous CPT from those colors automatically.
    ($(GMTdoc)plot3d.html#c)
- **D** | **offset** :: [Type => Str]

    Offset the plot symbol or line locations by the given amounts dx/dy.
    ($(GMTdoc)plot3d.html#d)
- **E** | **error_bars** :: [Type => Str]

    Draw symmetrical error bars.
    ($(GMTdoc)plot3d.html#e)
- **F** | **conn** | **connection** :: [Type => Str]

    Alter the way points are connected
    ($(GMTdoc)plot3d.html#f)
- **G** | **fill** | **markerfacecolor** | **MarkerFaceColor** | **markercolor** | **mc** :: [Type => Str]

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    ($(GMTdoc)plot3d.html#g)
- **I** | **intens** :: [Type => Str or number]

    Use the supplied intens value (in the [-1 1] range) to modulate the fill color by simulating illumination.
    ($(GMTdoc)plot3d.html#i)
- **L** | **closed_polygon** :: [Type => Str]

    Force closed polygons.
    ($(GMTdoc)plot3d.html#l)
- **N** | **no_clip** :: [Type => Str | []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)plot3d.html#n)
- $(GMT.opt_P)
- **S** | **symbol** | **marker** | **Marker** :: [Type => Str]

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    ($(GMTdoc)plot3d.html#s)
    Alternatively select a sub-set of symbols using the aliases: **marker** or **Marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
    + **y**, **y_dash**
- **W** | **pen** | **line_attribs** | **markeredgecolor** | **MarkerEdgeColor** | **mec**:: [Type => Str]
    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot3d.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **line_attribs** for the polygons and
    **markeredgecolor** or **MarkerEdgeColor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)

- **Z** | **level** :: [Type => Str | NamedTuple]	`Arg = value|file[+f|+l] | (data=Array|Number, outline=_, fill=_)`

    Paint polygons after the level given as a cte or a vector with same size of number of polygons. Needs a color map.
    ($(GMTdoc)plot3d.html#z)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_w)

### Example:

    plot3d(x -> sin(x)*cos(10x), y -> sin(y)*sin(10y), z -> cos(z), 0:pi/100:pi, show=true, aspect3=:equal)
"""
plot3d(arg1; kw...)  = common_plot_xyz("", cat_1_arg(arg1), "plot3d", true, true, kw...)
plot3d!(arg1; kw...) = common_plot_xyz("", cat_1_arg(arg1), "plot3d", false, true, kw...)

plot3d(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, arg1, "plot3d", true, true, kw...)
plot3d!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "plot3d", false, true, kw...)

# ------------------------------------------------------------------------------------------------------
function plot3d(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; first=true, kw...)
	common_plot_xyz("", hcat(arg1[:], arg2[:], arg3[:]), "plot3d", first, true, kw...)
end
plot3d!(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; kw...) = plot3d(arg1, arg2, arg3; first=false, kw...)

function plot3d(f1::Function, f2::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; kw...), "plot3d", first, true, kw...)
end
plot3d!(f1::Function, f2::Function, range_t=nothing; kw...) = plot3d(f1, f2, range_t; first=false, kw...)

function plot3d(f1::Function, f2::Function, f3::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_3f(f1, f2, f3, range_t; kw...), "plot3d", first, true, kw...)
end
plot3d!(f1::Function, f2::Function, f3::Function, range_t=nothing; kw...) = plot3d(f1, f2, f3, range_t; first=false, kw...)

const plot3  = plot3d			# Alias
const plot3! = plot3d!
# ------------------------------------------------------------------------------------------------------

"""
    scatter(cmd0::String="", arg1=nothing; kwargs...)

Reads (x,y) pairs and plot symbols at those locations on a map.
This module is a subset of ``plot`` to make it simpler to draw scatter plots. So many of
its (fine) controling parameters are not listed here. For a finer control, user should
consult the ``plot`` module.

Parameters
----------

- **G** | **fill** | **markerfacecolor** :: [Type => Str]

    Select color or pattern for filling of symbols or polygons.
    ($(GMTdoc)plot.html#g)
- **N** | **no_clip** :: [Type => Str | []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)plot.html#n)
- $(GMT.opt_P)
- **S** :: [Type => Str]

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    ($(GMTdoc)plot.html#s)

    Alternatively select a sub-set of symbols using the aliases: **symbol** or **marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
	+ **y**, **y_dash**

    and select their sizes with the **markersize** or **size** keyword [default is 8p].
    The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
    points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT=:c,)*
- **W** | **pen** | **markeredgecolor** | **mec** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)

[`Full man page`](https://genericmappingtools.github.io/GMT.jl/latest/scatter/)
[`GMT man page`]($(GMTdoc)plot.html)
"""
function scatter(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	common_plot_xyz("", cat_2_arg2(rang, [f(x) for x in rang]), "scatter",  first, false, kw...)
end
scatter!(f::Function, rang=nothing; kw...) = scatter(f, rang; first=false, kw...)

function scatter(f1::Function, f2::Function, range_t=nothing; first=true, kw...)	# Parametric version
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; is3D=false, kw...), "scatter",  first, false, kw...)
end
scatter!(f1::Function, f2::Function, range_t=nothing; kw...) = scatter(f1, f2, range_t; first=false, kw...)

scatter(cmd0::String="", arg1=nothing; first=true, kw...)  = common_plot_xyz(cmd0, arg1, "scatter",  first, false, kw...)
scatter!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "scatter",  false, false, kw...)

scatter(arg; kw...)  = common_plot_xyz("", cat_1_arg(arg), "scatter", true, false, kw...)
scatter!(arg; kw...) = common_plot_xyz("", cat_1_arg(arg), "scatter", false, false, kw...)

scatter(arg1, arg2; kw...)  = common_plot_xyz("", cat_2_arg2(arg1, arg2), "scatter", true, false, kw...)
scatter!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2), "scatter", false, false, kw...)

# ------------------------------------------------------------------------------------------------------
scatter3(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, arg1, "scatter3",  true, true, kw...)
scatter3!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "scatter3",  false, true, kw...)

scatter3(arg; kw...)  = common_plot_xyz("", cat_1_arg(arg), "scatte3", true, false, kw...)
scatter3!(arg; kw...) = common_plot_xyz("", cat_1_arg(arg), "scatte3", false, false, kw...)

function scatter3(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; kw...)
	common_plot_xyz("", hcat(arg1, arg2, arg3), "scatter3", true, true, kw...)
end
function scatter3!(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; kw...)
	common_plot_xyz("", hcat(arg1, arg2, arg3), "scatter3", false, true, kw...)
end

function scatter3(f1::Function, f2::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; kw...), "scatter3", first, true, kw...)
end
scatter3!(f1::Function, f2::Function, range_t=nothing; kw...) = scatter3(f1, f2, range_t; first=false, kw...)

function scatter3(f1::Function, f2::Function, f3::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_3f(f1, f2, f3, range_t; kw...), "scatter3", first, true, kw...)
end
scatter3!(f1::Function, f2::Function, f3::Function, range_t=nothing; kw...) = scatter3(f1, f2, f3, range_t; first=false, kw...)

const scatter3d  = scatter3			# Alias
const scatter3d! = scatter3!
# ------------------------------------------------------------------------------------------------------

"""
    bar(cmd0::String="", arg1=nothing; kwargs...)

Reads a file or (x,y) pairs and plots vertical bars extending from base to y.

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **fill** :: [Type => Str --

    Select color or pattern for filling the bars
    ($(GMTdoc)plot.html#g)
- **base** | **bottom** :: [Type => Str | Num]		``key=value``

    By default, base = ymin. Use this option to change that value. If base is not appended then we read it.
    from the last input data column.
- **size** | **width** :: [Type => Str | Num]		``key=value``

    The size or width is the bar width. Append u if size is in x-units. When *width* is used the default is plot-distance units.

Example:

    bar(sort(randn(10)), fill=:black, axis=:auto, show=true)
"""
function bar(cmd0::String="", arg=nothing; first=true, kw...)
	d = KW(kw)
	do_cat = (haskey(d, :stack) || haskey(d, :stacked) && isvector(arg) && length(arg) > 2) ? false : true
	if (cmd0 == "" && do_cat) arg = cat_1_arg(arg)  end	# If ARG is a vector, prepend it with a 1:N x column
	GMT.common_plot_xyz(cmd0, arg, "bar", first, false, kw...)
end
bar!(cmd0::String="", arg=nothing; kw...) = bar(cmd0, arg; first=false, kw...)

function bar(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	bar("", cat_2_arg2(rang, [f(x) for x in rang]); first=first, kw...)
end
bar!(f::Function, rang=nothing; kw...) = bar(f, rang; first=false, kw...)

bar(arg1, arg2; first=true, kw...)  = common_plot_xyz("", cat_2_arg2(arg1, arg2), "bar", first, false, kw...)
bar!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2), "bar", false, false, kw...)
bar(arg; kw...)  = bar("", arg; kw...)
bar!(arg; kw...) = common_plot_xyz("", cat_1_arg(arg), "bar", false, false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    bar3(cmd0::String="", arg1=nothing; kwargs...)

Read a grid file, a grid or a MxN matrix and plots vertical bars extending from base to z.

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **fill** :: [Type => Str]		``key=color``

    Select color or pattern for filling the bars
    ($(GMTdoc)plot.html#g)
- **base** :: [Type => Str | Num]		``key=value``

    By default, base = ymin. Use this option to change that value. If base is not appended then we read it.
- $(GMT.opt_p)

Example:

    G = gmt("grdmath -R-15/15/-15/15 -I0.5 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
    bar3(G, lw=:thinnest, show=true)
"""
function bar3(cmd0::String="", arg=nothing; first=true, kwargs...)
	# Contrary to "bar" this one has specific work to do here.
	d = KW(kwargs)
	opt_z = ""

	if (isa(arg, Array{<:GMTdataset,1}))  arg1 = arg[1]	# It makes no sense accepting > 1 datasets
	else                                  arg1 = arg	# Make a copy that may or not become a new thing
	end

	if (isa(arg1, Array))
		ny, nx = size(arg1)
		if ((nx >= 3 && ny > 3))  arg1 = mat2grid(arg1)  end		# Assume it is a 'bare grid'
	elseif (cmd0 != "")
		if ((val = find_in_dict(d, [:grd :grid])[1]) !== nothing)
			arg1 = gmtread(cmd0, grd=true)
		elseif ((val = find_in_dict(d, [:dataset :table])[1]) !== nothing)
			arg1 = gmtread(cmd0, dataset=true);		arg1 = arg1[1]
		else
			error("BAR3: When first arg is a name, must also state its type. e.g. grd=true or dataset=true")
		end
	end

	opt_base = add_opt(d, "", "", [:base])	# No need to purge because base is not a psxy option

	if (isa(arg1, GMTgrid))
		if (haskey(d, :bar))
			opt_S = parse_bar_cmd(d, :bar, "", "So")[1]
		else
			# 0.85 is the % of inc width of bars
			w1::Float64, w2::Float64 = arg1.inc[1]*0.85, arg1.inc[2]*0.85
			opt_S = " -So$(w1)u/$(w2)u"
			if     (haskey(d, :nbands))  opt_z = string("+z", d[:nbands]);	delete!(d, :nbands)
			elseif (haskey(d, :Nbands))  opt_z = string("+Z", d[:Nbands]);	delete!(d, :Nbands)
			end
		end
		opt, = parse_R(d, "", !first)
		if (opt == "" || opt == " -R")			# OK, no R but we know it here so put it in 'd'
			if (arg1.registration == 1)			# Fine, grid is already pixel reg
				push!(d, :R => arg1.range)
			else								# Need to get a pix reg R
				range = deepcopy(arg1.range)
				range[1] -= arg1.inc[1] / 2;	range[2] += arg1.inc[1] / 2;
				range[3] -= arg1.inc[2] / 2;	range[4] += arg1.inc[2] / 2;
				push!(d, :R => range)
			end
			z_min = arg1.range[5]
		elseif (opt_base == "")					# Shit, need to get zmin out of the opt_R string
			t = split(opt, '/')
			(length(t) == 6) ? z_min = t[5] : error("For 3D cases, region must have 6 selements")
		end
		if (opt_base == "")  push!(d, :base => z_min)  end	# Make base = z_min
		arg1 = gmt("grd2xyz", arg1)[1]			# Now arg1 is a GMTdataset
	else
		opt_S = parse_I(d, "", [:S :width], "So", true)
		if (opt_S == "")
			opt_S = parse_bar_cmd(d, :bar, "", "So", true)[1]
		end
		if (opt_S == "")
			if ((isa(arg1, Array) && size(arg1,2) < 5) || (isa(arg1, GMTdataset) && size(arg1.data,2) < 5))
				error("BAR3: When NOT providing *width* data must contain at least 5 columns.")
			end
		end
		#if (opt_S != "" && !isletter(opt_S[end]))   opt_S = opt_S * 'u'  end
		if     (haskey(d, :nbands))  opt_z = string("+z", d[:nbands]);	delete!(d, :nbands)
		elseif (haskey(d, :Nbands))  opt_z = string("+Z", d[:Nbands]);	delete!(d, :Nbands)
		end
	end

	opt_base = add_opt(d, "", "", [:base])		# Do this again because :base may have been added above
	if (opt_base == "")
		_z_min::Float32 = (isa(arg1, Array)) ? minimum(view(arg1, :, 3)) : minimum(view(arg1.data, :, 3))
		opt_S *= "+b$_z_min"
	else
		opt_S *= "+b" * opt_base
	end

	common_plot_xyz("", arg1, "bar3|" * opt_S * opt_z, first, true, d...)
end

bar3(arg1; kw...) = bar3("", arg1; first=true, kw...)
bar3!(cmd0::String="", arg1=nothing; kw...) = bar3(cmd0, arg1; first=false, kw...)
bar3!(arg1; kw...) = bar3("", arg1; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    arrows(cmd0::String="", arg1=nothing; arrow=(...), kwargs...)

Plots an arrow field. When the keyword *arrow=(...)* or *vector=(...)* is used, the direction (in degrees
counter-clockwise from horizontal) and length must be found in columns 3 and 4, and size, if not specified
on the command-line, should be present in column 5. The size is the length of the vector head. Vector stem
width is set by option *pen* or *line_attrib*.

The *vecmap=(...)* variation is similar to above except azimuth (in degrees east of north) should be
given instead of direction. The azimuth will be mapped into an angle based on the chosen map projection.
If length is not in plot units but in arbitrary user units (e.g., a rate in mm/yr) then you can use the
*input_col* option to scale the corresponding column via the +sscale modifier.

The *geovec=(...)* or *geovector=(...)* keywords plot geovectors. In geovectors azimuth (in degrees east from north) and geographical length must be found in columns 3 and 4. The size is the length of the vector head. Vector width is set by *pen* or *line_attrib*. Note: Geovector stems are drawn as thin filled polygons and hence pen attributes like dashed and dotted are not available. For allowable geographical units, see the *units=()* option.

The full *arrow* options list can be consulted at [Vector Attributes](@ref)

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)

Example:

	arrows([0 8.2 0 6], limits=(-2,4,0,9), arrow=(len=2,stop=1,shape=0.5,fill=:red), axis=:a, pen="6p", show=true)
"""
function arrows(cmd0::String="", arg1=nothing; first=true, kwargs...)
	# A arrows plotting method of plot
	d = KW(kwargs)
	cmd = helper_arrows(d, true)	# Have to delete to avoid double parsing in -W
	cmd = (cmd == "") ? " -Sv0.5+e+h0.5" : " -S" * cmd
	GMT.common_plot_xyz(cmd0, arg1, cmd, first, false, d...)
end

function helper_arrows(d::Dict, del::Bool=true)
	# Helper function to set the vector head attributes
	(show_kwargs[1]) && return print_kwarg_opts([:arrow :vector :arrow4 :vector4 :vecmap :geovec :geovector], "NamedTuple | String")

	cmd = ""
	val, symb = find_in_dict(d, [:arrow :vector :arrow4 :vector4 :vecmap :geovec :geovector], del)
	if (val !== nothing)
		code = 'v'
		if (symb == :geovec || symb == :geovector)
			code = '='
		elseif (symb == :vecmap)	# Uses azimuth and plots angles taking projection into account
			code = 'V'
		end
		if (isa(val, String))		# An hard core GMT string directly with options
			cmd = (val[1] != code) ? code * val : val	# In last case the GMT string already has vector flag char
		elseif (isa(val, Number))  cmd = code * "$val"
		elseif (symb == :arrow4 || symb == :vector4)  cmd = code * vector4_attrib(val)
		else                       cmd = code * vector_attrib(val)
		end
	end
	return cmd
end

arrows!(cmd0::String="", arg1=nothing; kw...) = arrows(cmd0, arg1; first=false, kw...)
arrows(arg1; kw...)  = arrows("", arg1; first=true, kw...)
arrows!(arg1; kw...) = arrows("", arg1; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    lines(cmd0::String="", arg1=nothing; decorated=(...), kwargs...)

Reads a file or (x,y) pairs and plots a collection of different line with decorations

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)

Examples:

    lines([0, 10]; [0, 20], limits=(-2,12,-2,22), proj="M2.5", pen=1, fill=:red,
		  decorated=(dist=(val=1,size=0.25), symbol=:box), show=true)

    lines(x -> cos(x) * x, y -> sin(y) * y, linspace(0,2pi,100), region=(-4,7,-5.5,2.5), lw=2, lc=:sienna,
          decorated=(quoted=true, const_label=" In Vino Veritas  - In Aqua, Rãs & Toads", font=(25,"Times-Italic"),
                     curved=true, pen=(0.5,:red)), aspect=:equal, fmt=:png, show=true)
"""
lines(cmd0::String="", arg1=nothing; first=true, kwargs...) = common_plot_xyz(cmd0, arg1, "lines", first, false, kwargs...)
lines!(cmd0::String="", arg=nothing; kw...) = lines(cmd0, arg; first=false, kw...)

function lines(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	lines("", cat_2_arg2(rang, [f(x) for x in rang]); first=first, kw...)
end
lines!(f::Function, rang=nothing; kw...) = lines(f, rang; first=false, kw...)

function lines(f1::Function, f2::Function, range_t=nothing; first=true, kw...)	# Parametric version
	lines("", help_parametric_2f(f1, f2, range_t; is3D=false, kw...); first=first, kw...)
end
lines!(f1::Function, f2::Function, range_t=nothing; kw...) = lines(f1, f2, range_t; first=false, kw...)

lines(arg1, arg2; kw...)  = lines("", cat_2_arg2(arg1, arg2); first=true, kw...)
lines!(arg1, arg2; kw...) = lines("", cat_2_arg2(arg1, arg2); first=false, kw...)
lines(arg; kw...)  = lines("", cat_1_arg(arg); first=true, kw...)
lines!(arg; kw...) = lines("", cat_1_arg(arg); first=false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    hlines(arg; decorated=(...), kwargs...)

Plots one or a collection of horizontal lines with eventual decorations

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for the horizontal lines
    ($(GMTdoc)plot.html#w)

Example:

    plot(rand(5,3))
    hlines!([0.2, 0.6], pen=(1, :red), show=true)
"""
function hlines(arg1=nothing; first=true, kwargs...)
	# A lines plotting method of plot
	d = KW(kwargs)
	(arg1 === nothing && ((arg1 = find_in_dict(d, [:data])[1]) === nothing)) && error("No input data")
	mat = ones(2, length(arg1))
	[mat[1,k] = mat[2,k] = arg1[k] for k = 1:length(arg1)]
	x = ((opt_R = parse_R(d, "")[2]) != "") ? vec(opt_R2num(opt_R)[1:2]) : [-1e50, 1e50]
	D = mat2ds(mat, x=x, multi=true)

	common_plot_xyz("", D, "lines", first, false, d...)
end
hlines!(arg=nothing; kw...) = hlines(arg; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------------------
"""
    vlines(arg; decorated=(...), kwargs...)

Plots one or a collection of vertical lines with eventual decorations

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for the horizontal lines
    ($(GMTdoc)plot.html#w)

Example:

    plot(rand(5,3), region=[0,1,0,1])
    vlines!([0.2, 0.6], pen=(1, :red), show=true)
"""
function vlines(arg1=nothing; first=true, kwargs...)
	# A lines plotting method of plot
	d = KW(kwargs)
	(arg1 === nothing && ((arg1 = find_in_dict(d, [:data])[1]) === nothing)) && error("No input data")
	mat = ones(2, length(arg1))
	mat[1,:] = mat[2,:] = arg1
	x = ((opt_R = parse_R(d, "")[2]) != "") ? vec(opt_R2num(opt_R)[3:4]) : [-1e50, 1e50]
	D = mat2ds(mat, x=x, multi=true)
	# Now we need tp swapp x / y columns because the vlines case is more complicated to implement.
	for k = 1:length(arg1)
		D[k].data[:,1], D[k].data[:,2] = D[k].data[:,2], D[k].data[:,1]
	end

	common_plot_xyz("", D, "lines", first, false, d...)
end
vlines!(arg=nothing; kw...) = vlines(arg; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------------------
"""
    ternary(cmd0="", arg1=nothing; image=false, clockwise=false, kwargs...)

Reads (a,b,c[,z]) records from table [or file] and plots image and symbols at those locations on a ternary diagram.

- **B** | **frame** :: [Type => NamedTuple | Str] --

    For ternary diagrams the three sides are referred to as a (bottom), b (right), and c (left). The default is to
    annotate and draw grid lines but without labeling the axes. But since labeling is a very important feature, you
    can use the `labels` option that take as argument a 3 elements Tuple with the labels of the 3 axes. Further
    control on annotations and grid spacing (on/off) is achieved by using the `frame=(annot=?, grid=?, alabel=?, blabel=?,
    clabel=?, suffix=?)` form. Note that not all options of the general `frame` options are accepted in this module and for more
    elaborated frame option selection you will have to resort to the pure GMT syntax in the form `frame="<arg> <arg> <arg>"`
    ($(GMTdoc)psternary.html#b)
- $(GMT.opt_C)
- **G** | **fill** :: [Type => Str] --

    Select color or pattern for filling the bars
    ($(GMTdoc)psternary.html#c)
- **L** | **vertex_labels** :: [Type => Str | Tuple of strings] --		`Arg = a/b/c`

    Set the labels for the three diagram vertices where the component is 100% [none].
    ($(GMTdoc)psternary.html#l)
- **M** | **dump** :: [Type => Str]

    Dumps the converted input (a,b,c[,z]) records to Cartesian (x,y,[,z]) records, where x, y
    are normalized coordinates on the triangle (i.e., 0–1 in x and 0–sqrt(3)/2 in y). No plotting occurs.
    ($(GMTdoc)coast.html#m)
- **N** | **no_clip** | **noclip** :: [Type => Str or []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)psternary.html#n)
- **R** | **region** | **limits** :: [Type => Tuple | Str]

    Give the min and max limits for each of the three axis a, b, and c. Default is (0,100,0,100,0,100)
- $(GMT.opt_P)
- **S** | **symbol** :: [Type => Str]

    Plot individual symbols in a ternary diagram. If `S` is not given then we will instead plot lines
    (requires `pen`) or polygons (requires `color` or `fill`).
    ($(GMTdoc)psternary.html#s)

    Alternatively select a sub-set of symbols using the aliases: **symbol** or **marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
    + **y**, **y_dash**

    and select their sizes with the **markersize** or **size** keyword [default is 8p].
    The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
    points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT=:c,)*
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **pen** :: [Type => Str | Number]

    Sets the attributes for the particular line.
    ($(GMTdoc)psternary.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_q)
- $(GMT.opt_t)

Other than the above options, the `kwargs` input accepts still the following options
- `image`: - Fills the ternary plot with an image computed automatically with `grdimage` from a grid interpolated with `surface`
- `contour`: - This option works in two different ways. If used together with `image` it overlays a contour
               by doing a call to `grdcontour`. However, if used alone it will call `contour` to do the contours.
               The difference is important because this option can be used in *default mode* with `contour=true`
               where the number and annotated contours is picked automatically, or the use can exert full control
               by passing as argument a NamedTuple with all options appropriated to that module. *e.g.*
               `contour=(cont=10, annot=20, pen=0.5)`
- `contourf`: - Works a bit like the _standalone_ `contour`. If used with `contourf=true` call make a filled contour
                using automatic parameters. The form `contourf=(...)` let us selects options of the contourf module.
- `clockwise`: - Set it to `true` to indicate that positive axes directions be clock-wise
                 [Default lets the a, b, c axes be positive in a counter-clockwise direction].
"""
function ternary(cmd0::String="", arg1=nothing; first::Bool=true, image::Bool=false, kwargs...)
	# A wrapper for psternary
	(cmd0 == "" && arg1 === nothing) && (arg1 = [0.0 0.0 0.0])	# No data in, just a kind of ternary basemap
	(cmd0 != "") && (arg1 = gmtread(cmd0))
	d = init_module(first, kwargs...)[1]
	opt_J = parse_J(d, "", " -JX" * split(def_fig_size, '/')[1] * "/0", true, false, false)[2]
	opt_R = parse_R(d, "")[1]
	d[:R] = (opt_R ==  "") ? "0/100/0/100/0/100" : opt_R[4:end]
	parse_B4ternary!(d)
	clockwise = haskey(d, :clockwise)
	if (image || haskey(d, :contourf) || haskey(d, :contour))
		t = tern2cart(isa(arg1, GMTdataset) ? arg1.data : isa(arg1, Vector{<:GMTdataset}) ? arg1[1].data : arg1, clockwise)
		!endswith(opt_J, "/0") && (opt_J *= "/0")			# Need the "/0". Very important.
		if (haskey(d, :contourf))
			contourf(t, R=(0.0,1.0,0,sqrt(3)/2), B=:none, J=opt_J[4:end], backdoor=d[:contourf])
			delete!(d, :contourf)
		else
			G = surface("-R0/1/0/0.865 -I0.005 -T0.5 -Vq", t)
			Gmask = grdmask("-R0/1/0/0.865 -I0.005 -NNaN/1/1", [0.0 0; 0.5 0.865; 1 0; 0 0])
			G *= Gmask
			if (image)			# grdimage plus eventual contours
				grdimage(G, B=:none, J=opt_J[4:end])
				if (haskey(d, :contour))
					grdcontour!(G, backdoor=d[:contour])
					delete!(d, :contour)
				end
			else				# Only contours
				pscontour(t, R=(0.0,1.0,0,sqrt(3)/2), B=:none, J=opt_J[4:end], backdoor=d[:contour])
				delete!(d, :contour)
			end
		end
		first = false
	end
	if (clockwise)
		endswith(opt_J, "/0") && (opt_J = opt_J[1:end-2])		# Strip the trailing "/0". Very important.
		d[:J] = "X-" * opt_J[5:end]
		del_from_dict(d, [:proj :projection])		# To avoid non-consumed warnings
		delete!(d, :clockwise)
	end
	if ((val = find_in_dict(d, [:par :conf :params], false)[1]) === nothing)
		d[:par] = (MAP_GRID_PEN_PRIMARY="thinnest,gray",)
	end
	if (GMTver <= v"6.2.0" && (val = find_in_dict(d, CPTaliases, false)[1]) !== nothing && isa(val, GMTcpt))
		_name = joinpath(tempdir(), "GMTjl_tmp.cpt");	gmtwrite(_name, val);	d[:C] = _name	# Workaround a bug in 6.2.0
	end
	(API == C_NULL) && gmt_restart()	# Force having a valid API. We can't afford otherwise here.
	(GMTver <= v"6.2.0") && gmtlib_setparameter(API, "MAP_FRAME_AXES", "WESNZ")	# Because of a bug in 6.2.0 modern theme
	r = common_plot_xyz("", arg1, "ternary", first, false, d...)
	(GMTver <= v"6.2.0") && gmtlib_setparameter(API, "MAP_FRAME_AXES", "auto")
	# With the following trick we leave the -R history in 0/1/0/1 and so we can append with plot, text, etc
	gmt("psxy -Scp -R0/1/0/1 -JX -O -Vq > " * joinpath(tempdir(), "lixo.ps"), [0. 0.])
	return r
end

function parse_B4ternary!(d::Dict)
	# Ternary accepts only a special brand of -B. Try to parse and/or build -B option
	opt_B = parse_B(d, "", " -Bafg")[2]
	if ((val = find_in_dict(d, [:labels])[1]) !== nothing)		# This should be the easier way
		!(isa(val,Tuple) && length(val) == 3) && error("The `labels` option must be Tuple with 3 elements.")
		opt_Bs = split(opt_B)							# This drops the leading ' '
		x = (opt_Bs[1][3] == 'p') ? opt_Bs[1][4:end] : opt_Bs[1][3:end]
		d[:B] = " -Ba$(x)+l" * string(val[1]) * " -Bb$(x)+l" * string(val[2]) * " -Bc$(x)+l" * string(val[3])
		[d[:B] *= " " * opt_Bs[k] for k = 2:length(opt_Bs)]		# Append the remains, if any.
	else		# Ui, try to parse a string like this: " -Bpag8+u\" %\" -Ba+la -Bb+lb -Bc+lc"
		opt_Bs = split(opt_B, " -B")[2:end]				# 2:end because surprisingly the first is = ""
		if (length(opt_Bs) == 1)  d[:B] = opt_B			# Accept whatever was selected
		else											# User may have used frame=(annot=?,grid=?, alabel=?,...)
			if (length(opt_Bs) == 3) d[:B] = opt_B		# OK, silly no annotations,ticks,grid
			else
				x = opt_Bs[1][2:end]
				d[:B] = " -Ba$(x)" * opt_Bs[2][2:end] * " -Bb$(x)" * opt_Bs[3][2:end] * " -Bc$(x)" * opt_Bs[4][2:end]
				[d[:B] *= " -B" * opt_Bs[k] for k = 5:length(opt_Bs)]		# Append the remains, if any.
			end
		end
	end
end

function tern2cart(abcz::Matrix{<:Real}, reverse::Bool=false)
	# coverts ternary to cartesian units. ABCZ is either a Mx3 (a,b,c) or Mx4 (a,b,c,z) matrix
	a,b,c = !reverse ? (3,1,2) : (1,2,3)
	s = view(abcz, :, a) + view(abcz, :, b) + view(abcz, :, c)	# s = (a + b + c)
	x = 0.5 .* (2.0 .* view(abcz, :, b) + view(abcz, :, c)) ./ s
	y = 0.5 .* sqrt(3.) .* view(abcz, :, c) ./ s
	return (size(abcz,2) == 3) ? [x y] : [x y abcz[:,4]]
end

function dict_auto_add!(d::Dict)
	# If the Dict 'd' has a 'backdoor' member that is a NamedTuple, add its contents to the Dict
	if ((val = find_in_dict(d, [:backdoor])[1]) !== nothing && isa(val, NamedTuple))
		key = keys(val)
		[d[key[n]] = val[n] for n = 1:length(val)]
		delete!(d, :backdoor)
	end
end

ternary!(cmd0::String="", arg1=nothing; kw...) = ternary(cmd0, arg1; first=false, kw...)
ternary(arg1;  kw...) = ternary("", arg1; first=true, kw...)
ternary!(arg1; kw...) = ternary("", arg1; first=false, kw...)
ternary(kw...) = ternary("", nothing; first=true, kw...)
const psternary  = ternary            # Aliases
const psternary! = ternary!           # Aliases

# ----------------------------------------------------------------------------------------------
"""
    events(cmd0::String, arg1=nothing; kwargs...)

Plot event symbols and labels for a moment in time

Full option list at [`events`]($(GMTdoc)events.html)

Parameters
----------

- **T** | **now** :: [Type => Int | Str]

    Set the current plot time. If absolute times are used you must also use -fT.
    ($(GMTdoc)events.html#t)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** | **shift** | **offset** :: [Type => Str]	``Arg = [j|J]dx[/dy][+v[pen]]``

    Offset the text from the projected (x,y) point by dx,dy [0/0].
    ($(GMTdoc)events.html#d)
- **E** | **knots** :: [Type => Str]	``Arg = s|t[+o|Odt][+rdt][+pdt][+ddt][+fdt][+ldt]``

    Set the time knots for the symbol or text time-functions.
    ($(GMTdoc)events.html#e)
- **G** | **fill** :: [Type => Str | Int | Touple]

    Set constant shade or color for all symbols.
    ($(GMTdoc)events.html#g)
- $(GMT.opt_J)
- **L** | **duration** :: [Type => Bool | Number | Str]		``Arg = [length|t]``

    Specify the length (i.e., duration) of the event.
    ($(GMTdoc)events.html#l)
- **M** | **rise** :: [Type => Str]		``Arg = i|s|t[val1][+cval2]``

    Modify the initial intensity, size magnification, or transparency of the symbol during the rise interval.
    ($(GMTdoc)events.html#m)
- **Q** | **save** :: [Type => Number]

    Save the intermediate event symbols and labels to permanent files instead of removing them when done.
    ($(GMTdoc)events.html#q)
- $(GMT.opt_R)
- **W** | **pen** | **markeredgecolor** | **mec** :: [Type => Str]

    Specify symbol outline pen attributes [Default is no outline].
    ($(GMTdoc)events.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_swap_xy)
"""
# ------------------------------------------------------------------------------------------------------
function events(cmd0::String="", arg1=nothing; kwargs...)
	# events share a lot of options with plot
	d = KW(kwargs)
	cmd = add_opt(d, "", "T", [:T :now])
	if (!occursin("-T", cmd))  error("The 'now' (T) option is mandatory")  end
	cmd = add_opt(d, cmd, "E", [:E :knots],
		(symbol=("s", nothing, 1), text=("t", nothing, 1), shift_startEnd = "+o", shift_start="+O", raise="+r", plateau="+p", decay="+d", fade="+f", text_duration="+l"))
	cmd = add_opt(d, cmd, "M", [:M :rise],
		(intensity=("i", arg2str, 1), size=("s", arg2str, 1), transparency=("t", arg2str, 1), coda="+c"))
	cmd = add_opt(d, cmd, "L", [:L :duration])
	cmd = add_opt(d, cmd, "Q", [:Q :save])
	cmd = add_opt(d, cmd, 'D', [:D :offset],
		(away=("j", nothing, 1), corners=("J", nothing, 1), shift="", line=("+v",add_opt_pen)))
	cmd = add_opt(d, cmd, 'F', [:F :attrib],
		(angle="+a", Angle="+A", font=("+f", font), justify="+j", region_justify="+c", header="_+h", label="_+l", rec_number="+r", text="+t", zvalues="+z"), false)
	common_plot_xyz(cmd0, arg1, "events|" * cmd, true, false, d...)
end
const psevents = events            # Alias

# ------------------------------------------------------------------------------------------------------
cat_1_arg(arg::GMTdataset) = return arg				# Miserable attempts to force type stability
cat_1_arg(arg::Vector{<:GMTdataset}) = return arg
function cat_1_arg(arg)
	# Add a first column with 1:n to all args that are not GMTdatasets
	#(isa(arg, Vector{<:GMTdataset}) || isa(arg, GMTdataset))  &&  return arg
	if (isa(arg, Vector) || typeof(arg) <: AbstractRange)
		#arg = hcat(collect(1:size(arg,1)), arg)
		arg = hcat(collect(eltype(arg), 1:size(arg,1)), arg)
	#elseif (isa(arg, Array) && size(arg,1) == 1)		# Accept also row arrays. CAN'T IT BREAKS plot([1 1])
		#arg = hcat(collect(eltype(arg), 1:length(arg)), arg')
	elseif (isa(arg, NTuple))
		arg = hcat(collect(eltype(arg), 1:length(arg)), collect(arg))
	end
	return arg
end

# ------------------------------------------------------------------------------------------------------
function cat_2_arg2(arg1, arg2)
	# Cat two vectors (or tuples) or a vector (or tuple) and a matrix in a Mx2 matrix

	arg2 === nothing && return arg1
	isa(arg1, Real) && isa(arg2, Real) && return [arg1 arg2]
	!((isa(arg1, Vector) || typeof(arg1) <: AbstractRange || isa(arg1, NTuple) || isa(arg1, Matrix)) && (isa(arg2, Vector) || typeof(arg2) <: AbstractRange || isa(arg2, NTuple) || isa(arg2, Matrix)) ) &&
		error("Unknown types ($(typeof(arg1))) and ($(typeof(arg2))) in cat_2_arg2() function")

	if (isa(arg1, NTuple))  arg1 = collect(arg1)  end
	if (isa(arg2, NTuple))  arg2 = collect(arg2)  end
	if (size(arg1,1) == 1 && size(arg1,2) != 1)  arg1 = arg1[:]  end
	if (size(arg2,1) == 1 && size(arg2,2) != 1)  arg2 = arg2[:]  end
	arg = hcat(arg1, arg2)
	if (size(arg,2) > 2)  global multi_col[1] = true  end
	return arg
end

# ------------------------------------------------------------------------------------------------------
function gen_coords4funs(rang=nothing, axis="x"; kw...)
	# Generate axes coordenates to use when plot functions

	if (rang === nothing)
		symb = (axis == "x") ? [:xlim] : [:ylim]
		if ((val = find_in_dict(KW(kw), symb)[1]) !== nothing)
		(!isa(val, Tuple) && !isa(val, Array{<:Real})) && error("$(string(symb[1])) must be a tuple or array, not '$(typeof(val))'")
			(length(val) != 2) && error("$(string(symb[1])) must have 2 elements")
			rang = linspace(val[1], val[2], 200)
		else
			rang = linspace(-5, 5)
		end
	elseif (isa(rang, Real))
		rang = linspace(-rang, rang)
	end
	return rang
end

# ------------------------------------------------------------------------------------------------------
# Common code shared by the functions that accept parametric equations
function help_parametric_2f(f1::Function, f2::Function, range_t=nothing; is3D=true, kw...)
	# This function is shared by both the 2D & 3D cases
	t = collect(gen_coords4funs(range_t, "x"; kw...));
	x = [f1(x) for x in t];		y = [f2(y) for y in t]
	out = (is3D) ? hcat(x[:], y[:], t[:]) : hcat(x[:], y[:])
	return out
end

function help_parametric_3f(f1::Function, f2::Function, f3::Function, range_t=nothing; kw...)
	t = collect(gen_coords4funs(range_t, "x"; kw...));
	x = [f1(x) for x in t];		y = [f2(y) for y in t];		z = [f3(z) for z in t]
	return hcat(x[:], y[:], z[:])
end
# ------------------------------------------------------------------------------------------------------


text:
using Test
using Random
using VectorCalculus

Random.seed!(1234)

@testset "VectorCalculus" begin
    using LinearAlgebra
    include("algebraic.jl")
    include("differential.jl")
end

ext: jl


text:
"""
Surface area scaling rules
"""
abstract type AbstractScaling end

"""
    scaling_correction(f::AbstractScaling, V)

Calculate the shape/scaling correction coefficient
from the current mass `V`.
"""
function scaling_correction end
scaling_correction(f::Nothing, V) = 1

"""
    Isomorph()
"""
struct Isomorph <: AbstractScaling end

scaling_correction(f::Isomorph, V) = 1

"""
    V0morph(Vd)

$(FIELDDOCTABLE)
"""
@columns struct V0morph{Mo} <: AbstractScaling
#   Field    | Def | Unit | Bounds        | Log | Description
    Vd::Mo   | 4.0 | mol  | (0.0, 1000.0) | _   | "reference"
end

scaling_correction(f::V0morph, V) = (V / f.Vd)^(-2//3)

"""
    V1morph(Vd)

$(FIELDDOCTABLE)
"""
@columns struct V1morph{Mo} <: AbstractScaling
#   Field    | Def | Unit | Bounds        | Log | Description
    Vd::Mo   | 4.0 | mol  | (0.0, 1000.0) | _   | "reference"
end

scaling_correction(f::V1morph, V) = (V / f.Vd)^(1//3)

"""
    V1V0morph(Vd, Vmax, β)

$(FIELDDOCTABLE)
"""
@columns struct V1V0morph{Mo,B} <: AbstractScaling
#   Field    | Def | Unit | Bounds        | Log | Description
    Vd::Mo   | 4.0 | mol  | (0.0, 1000.0) | _   | "reference"
    Vmax::Mo | 4.0 | mol  | (0.0, 1000.0) | _   | "reference"
    β::B     | 4.0 | mol  | (0.0, 10.0)   | _   | "reference"
end

scaling_correction(f::V1V0morph, V) = (V / f.Vd)^(1//3 - (V/f.Vmax)^f.β)

"""
    Plantmorph(M_Vref, M_Vscaling)

Plant morph formulation from DEBtool.

$(FIELDDOCTABLE)
"""
@columns struct Plantmorph{Mo} <: AbstractScaling
#   Field          | Def | Unit | Bounds           | Log  | Description
    M_Vref::Mo     | 0.1 | mol  | (0.00001, 20.0)  | true | "Scaling reference"
    M_Vscaling::Mo | 1.0 | mol  | (0.0001, 2000.0) | true | "Scaling mass"
end

scaling_correction(f::Plantmorph, V) = (V / f.M_Vref)^(-V / f.M_Vscaling)

update_scaling!(o, u) = set_scaling!(o, scaling_correction(scaling_pars(o), u[:V]))


text:
module Recursion

using Mixtape
import Mixtape: CompilationContext, transform, allow
using CodeInfoTools
using MacroTools

module Factorial

f(x::Int64) = x <= 1 ? 1 : x * f(x - 1)

end

struct MyMix <: CompilationContext end
allow(ctx::MyMix, m::Module) = m == Factorial

swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        isexpr(s, :call) || return s
        s.args[1] == Base.:(*) || return s
        return Expr(:call, Base.:(+), e.args[2:end]...)
    end
    return new
end

function transform(::MyMix, src)
    b = CodeInfoTools.Builder(src)
    for (v, st) in b
        b[v] = swap(st)
    end
    return CodeInfoTools.finish(b)
end

function postopt!(::MyMix, ir)
    display(ir)
    ir
end

Mixtape.@load_abi()
display(call(Factorial.f, 10; ctx = MyMix()))

end # module


text:
module HerbProlog

using Base: UInt16, String
include("structs.jl")
include("swipl-bindings.jl")
include("swipl-interface.jl")

export Prolog, Swipl
export start, stop, cleanup, is_initialised, asserta, assertz, query, has_solution

module LowLevel

export swipl_start, swipl_is_initialised, swipl_start, swipl_cleanup, swipl_new_atom, swipl_atom_chars, swipl_new_functor, swipl_functor_name, swipl_functor_arity
export swipl_term_type, swipl_is_variable, swipl_is_ground, swipl_is_atom, swipl_is_string, swipl_is_integer, swipl_is_float, swipl_is_callable, swipl_is_compound, swipl_is_list, swipl_is_pair, swipl_is_atomic, swipl_is_number
export swipl_get_integer, swipl_get_bool, swipl_get_atom_chrs, swipl_get_float, swipl_get_name_arity, swipl_get_compound_name_arity, swipl_get_arg
export swipl_new_term_ref, swipl_new_term_refs, swipl_put_variable, swipl_put_atom, swipl_put_bool, swipl_put_atom_chars, swipl_put_string_chars, swipl_put_integer, swipl_put_float, swipl_put_functor, swipl_put_list, swipl_put_nil, swipl_put_term, swipl_cons_functor, swipl_copy_term_ref
export swipl_unify, swipl_unify_atom, swipl_unify_bool, swipl_unify_atom_chars, swipl_unify_string_chars, swipl_unify_integer, swipl_unify_float, swipl_unify_functor, swipl_unify_compound, swipl_unify_list, swipl_unify_nil, swipl_unify_arg
export swipl_predicate, swipl_open_query, swipl_next_solution, swipl_cut_query, swipl_close_query

end

end # module


text:
using DiffEqSensitivity, OrdinaryDiffEq, Zygote, Test
function lv!(du, u, p, t)
    x,y = u
    a, b, c, d = p
    du[1] = a*x - b*x*y
    du[2] = -c*y + d*x*y
end
function test(u0,p)
    tspan = [0.,1.]
    prob = ODEProblem(lv!, u0, tspan, p)
    sol = solve(prob,Tsit5())
    return sol.u[end][1]
end
function test2(u0,p)
    tspan = [0.,1.]
    prob = ODEProblem(lv!, u0, tspan, p)
    sol = solve(prob,Tsit5())
    return Array(sol)[1,end]
end

u0 = [1.,1.]
p = [1.,1.,1.,1.]
@test Zygote.gradient(test,u0,p) == Zygote.gradient(test2,u0,p)


text:
# Passing Julia callback functions to Python routines.
#
# Note that this will typically involve two functions: the
# desired Julia function/closure, and a top-level C-callable
# wrapper function used with PyCFunction_NewEx -- the latter
# is called from Python and calls the former as needed.

################################################################

# Define a Python method/function object from f(PyPtr,PyPtr)::PyPtr.
# Requires f to be a top-level function.
function pymethod(f::Function, name::String, flags::Integer)
    # Python expects the PyMethodDef structure to be a *constant*,
    # and the strings pointed to therein must also be constants,
    # so we define anonymous globals to hold these
    def = gensym("PyMethodDef")
    defname = gensym("PyMethodDef_ml_name")
    @eval const $defname = bytestring($name)
    @eval const $def = PyMethodDef(convert(Ptr{Uint8}, $defname),
                                   $(cfunction(f, PyPtr, (PyPtr,PyPtr))),
                                   convert(Cint, $flags),
                                   convert(Ptr{Uint8}, C_NULL))
    PyObject(@pycheckn ccall((@pysym :PyCFunction_NewEx), PyPtr,
                             (Ptr{PyMethodDef}, Ptr{Void}, Ptr{Void}),
                             &eval(def), C_NULL, C_NULL))
end

################################################################

# To pass an arbitrary Julia function to Python, we wrap it
# in a jl_Function Python class, where jl_Function.__call__
# executes the jl_Function_callback function in Julia, which
# in turn fetches the actual Julia function from the "f" attribute
# and calls it.

function jl_Function_call(self_::PyPtr, args_::PyPtr, kw_::PyPtr)
    ret_ = convert(PyPtr, C_NULL)
    args = PyObject(args_)
    try
        if kw_ != C_NULL
            throw(ArgumentError("keywords not yet supported in callbacks"))
        end
        f = unsafe_pyjlwrap_to_objref(self_)::Function
        ret = PyObject(f(convert(PyAny, args)...))
        ret_ = ret.o
        ret.o = convert(PyPtr, C_NULL) # don't decref
    catch e
        pyraise(e)
    finally
        args.o = convert(PyPtr, C_NULL) # don't decref
    end
    return ret_::PyPtr
end
const jl_Function_call_ptr = cfunction(jl_Function_call,
                                       PyPtr, (PyPtr,PyPtr,PyPtr))

jl_FunctionType = PyTypeObject()

function pycallback_initialize()
    global jl_FunctionType
    if (jl_FunctionType::PyTypeObject).tp_name == C_NULL
        jl_FunctionType::PyTypeObject =
         pyjlwrap_type("PyCall.jl_Function",
                       t -> t.tp_call = jl_Function_call_ptr)
    end
    return
end

function pycallback_finalize()
    global jl_FunctionType
    jl_FunctionType::PyTypeObject = PyTypeObject()
end

function pycallback(f::Function)
    global jl_FunctionType
    if (jl_FunctionType::PyTypeObject).tp_name == C_NULL
        pycallback_initialize()
    end
    pyjlwrap_new(jl_FunctionType::PyTypeObject, f)
end

PyObject(f::Function) = pycallback(f)


text:
#!/usr/bin/julia
using PyPlot, QuantumOptics, LinearAlgebra, FFTW, Colors,
        Printf, DelimitedFiles

# make sure to set script directory as pwd()
cd(@__DIR__)

# include my custom cmds module
if Sys.iswindows()
    include("..\\cmds.jl")
    fn = "01_Output\\"
else
    include("../cmds.jl")
    fn = "01_Output/"
end

# to use PyPlot in GUI mode under Linux
pygui(true)
# run once with calc_2d = false to initialize functions
calc_2d = true

cmp = cmds.create_colormap("bright");

function draw_dipole(x,y,α,d)
    quiver(x- d/2*cos(α), y - d/2*sin(α), d*cos(α), d*sin(α), angles="xy",
            scale_units="xy",scale=1,color="r")
    scatter(x,y)
end

function angle2rad(α)
    A = α / 360 * 2 * π
    return A
end

## geometry of coupled dimer
# here, I use (0,0) as origin for point dipole 1 and an angle of the dipole of 0, such
# that it lies along the x-axis and has a length d. The point dipole is located at
# a distance r under the angle β, and has an dipole angle α and a length d. d
# controls the coupling strength.
#
# position of dipole 1 and angle wrt to x-axis (random)
x₁, y₁, α₁, d₁ = 0, 0, angle2rad(0), 0.25
# r: distance, β: slip angle
r , β          = 1.155, angle2rad(23.147495) # slip angle for TDBC 0.404 rad https://arxiv.org/pdf/1202.5712.pdf (theory)
# position of dipole 2 (calculated)
x₂, y₂         = x₁ + r*cos(β), y₁ + r*sin(β)
# angle and moment of dipole 2
α₂, d₂         = angle2rad(0), 0.25

## set up energies and couplings
# homo dimer    : E1 == E2 and J1 = J2
# hetero dimer  : E1 != E2
E₁ = 2.396
E₂ = 2.396

## draw dimer
#
figure(figsize=(10,6));
ax1 = subplot(331);
grid(b=1,which="both")
plot([x₁,x₂], [y₁,y₂], linestyle="--", color="k")
draw_dipole(x₁,y₁,α₁,d₁); draw_dipole(x₂,y₂,α₂,d₂)
xlim([-0.5, x₂+.5]), ylim([-0.5, y₂+0.5])
title("Dimer geometry")
xlabel("X"); ylabel("Y")

## determine x and y-components of tdm 1 and 2
d₁xy = [cos(α₁)*d₁, sin(α₁)*d₁] # x and y component of point dipole 1
d₂xy = [cos(α₂)*d₂, sin(α₂)*d₂] # x and y component of point dipole 2

## coupling strength will depend on geometry, first cos is out-of-plane rotation
# WHAT IS THE FACTOR (here 0.4) in κ ???
κ = 1 * (cos(0) - 3 * cos(β) * cos(β))
J = κ * (d₁*d₂) / norm([x₁,y₁]-[x₂,y₂])^3

# the following from Cho's book (p.170)
R12 = [x₁,y₁]-[x₂,y₂]
J = dot(d₁xy,d₂xy) / norm(R12)^3 -
        3 * dot(d₁xy,R12) * dot(R12,d₂xy) / norm(R12)^5

## make Hilbert-space
b_mon = NLevelBasis(2)  # Hilbert-space of monomer
b_dim = b_mon ⊗ b_mon   # Hilbert-space of dimer

## create transition operators in monomer basis
j12 = transition(b_mon,1,2) # == g_e, or a / σ₋ (annihilation)
j21 = transition(b_mon,2,1) # == e_g, or a† /σ⁺ (creation)
# g_e / e_g notation from Constantin script

## convert transition operators into dimer basis
mu1hat = j21 ⊗ one(b_mon) + j12 ⊗ one(b_mon)
mu2hat = one(b_mon) ⊗ j21 + one(b_mon) ⊗ j12

## add geometry and magnitude of tdm and calculate total tdm in x and y
μx = d₁xy[1] * mu1hat + d₂xy[1] * mu2hat
μy = d₁xy[2] * mu1hat + d₂xy[2] * mu2hat

## visualize combined tdm (offset in y a bit for clarity)
quiver((x₁+x₂)/2-(d₁xy[1]+d₂xy[1])/2, (y₁+y₂)/2-(d₁xy[2]+d₂xy[2])/2+0.05,
        d₁xy[1]+d₂xy[1], d₁xy[2]+d₂xy[2], angles="xy", scale_units="xy",
        scale=1, color="b", alpha=0.75)

quiver((x₁+x₂)/2-(d₁xy[1]-d₂xy[1])/2, (y₁+y₂)/2-(d₁xy[2]-d₂xy[2])/2+0.05,
        d₁xy[1]-d₂xy[1], d₁xy[2]-d₂xy[2], angles="xy", scale_units="xy",
        scale=1, color="r")

scatter((x₁+x₂)/2,(y₁+y₂)/2)

## define Hamiltonian of the dimer
# H₁ and H₂ are the Hamiltonians of the individual dipoles
H₁ = E₁ * j21 * j12
H₂ = E₂ * j21 * j12
# Hₓ describes the exchange between dipoles 1 and 2
Hₓ = J * (j12 ⊗ j21 + j21 ⊗ j12)
# use the following with TDMs between ground- and biexciton state, this modifies
# the energies of the ground- and doubly excited state slightly
#Hₓ = J * ((j21+j12) ⊗ (j21+j12)) #???????????
# H is the total Hamiltonian of the system
H_dim = one(b_mon) ⊗ H₁ + H₂ ⊗ one(b_mon) + Hₓ

## ADD CAVITY

wcav = 2.2 # cavity frequency
g  = 0.15   # coupling strength
κ = 0.5         # cavity dissipation rate
N = 2               # number of cavity fock states
b_fock = FockBasis(N)
a = destroy(b_fock) ⊗ one(b_dim)
at = sparse(dagger(a))

J12a = j12 ⊗ one(j12)
J12b = one(j12) ⊗ j12
J21a = j21 ⊗ one(j21)
J21b = one(j21) ⊗ j21

sm = one(b_fock) ⊗ (j12 ⊗ j21)
sp = sparse(dagger(sm))

A = at * (one(b_fock) ⊗ J12a) + at * (one(b_fock) ⊗ J12b)
B = a  * (one(b_fock) ⊗ J21a) + a  * (one(b_fock) ⊗ J21b)

H_cav = wcav * at * a

H = H_cav + one(b_fock) ⊗ H_dim + g * (A + B)

## get eigenstates and eigenvalues of H
#  using Julia LinearAlgebra functions instead of QuantumOptics requires to use H.data
states_dim   = eigvecs(dense(H_dim).data)
energies_dim = eigvals(dense(H_dim).data)

states_cav   = eigvecs(dense(H_cav).data)
energies_cav = eigvals(dense(H_cav).data)

states   = eigvecs(dense(H).data)
energies = eigvals(dense(H).data)

#states = eigenstates(dense(H))

## perform bases transformation of μ
#  for this I need to express the matrices as dense, data matrices
#μx =  states * μx * states
#μy =  states' * dense(μy) * states

## plot site basis and exciton basis energy diagrams
ax2 = subplot(132)
hlines([0, E₁], -1.25, -0.75)
hlines([0, E₂],-0.25, 0.25)
hlines(energies_dim[:], 0.75, 1.25)
hlines(energies[:], 1.75, 2.25)
hlines(energies_cav[:], 2.75, 3.25)
xticks([-1, 0, 1, 2, 3], ["monomer 1", "monomer 2", "dimer", "dimer in cavity", "cavity"],
        rotation=45, ha="right")
title("Energy levels")

## make initial density matrix
Psi0 = fockstate(b_fock,0) ⊗ nlevelstate(b_mon,1) ⊗ nlevelstate(b_mon,1)
rho0 = dm(Psi0)
println("rho0")
display(rho0)

## time list for evaluation of master equation
tlist = [0:0.2:150;]

## make collapse operator
ccc = diagonaloperator(b_mon,[1, 1])
#ccc = one(b_mon)
#c_ops = [sqrt(0.25)*one(b_mon)⊗j12, sqrt(0.25)*j12⊗one(b_mon)]
L = [one(b_fock)⊗ccc⊗j12, one(b_fock)⊗j12⊗ccc, a]

Γ_L = [sqrt(0.005), sqrt(0.005),               sqrt(0.02)]
Γ_a = [sqrt(0.2), sqrt(0.2), sqrt(0.2), sqrt(0.2)]

L = Γ_L .* L
## make transition dipole operator FOR DIMER SYSTEM !!!!!
# clearly now μ has a direction dependence, but if sample is isotropic  take:
μ   = (μx + μy)
# and the tdm with the ground state:
#μ12 = one(b_fock) ⊗ (μ .* (j12*j21⊗j21 + j12*j21⊗j12 + j21⊗(j12*j21) + j12⊗(j12*j21)))
#display(dense(μ12))
# ... and with the doubly excited state:
#μ23 = one(b_fock) ⊗ (μ .* (j21*j12⊗j21 + j21*j12⊗j12 + j21⊗(j21*j12) + j12⊗(j21*j12)))
#display(dense(μ23))
CC = ccc ⊗ ccc
## make transition dipole operator FOR CAVITY !!!!!
μ12 = at * (one(b_fock) ⊗ (J12a * J21a * CC)) + a * (one(b_fock) ⊗ (J12a * J21a * CC))
μ23 = at * (one(b_fock) ⊗ (J21a * J12a * CC)) + a * (one(b_fock) ⊗ (J21a * J12a * CC)) +
      at * (one(b_fock) ⊗ (J21b * J12b * CC)) + a * (one(b_fock) ⊗ (J21b * J12b* CC))
μ = μ12 + μ23

rho1 = μ12 * rho0 * μ12
rho1 = μ12 * rho0 * μ12

rho2 = μ23 * rho1 * μ23
μ23  = μ23 / sqrt(tr(rho1))
rho2 = μ23 * rho1 * μ23
μ12  = μ12 / sqrt(tr(rho1))

# calculate magnitude of tdm with singly excited states e1 and e2
# HOT TO ORDER THEM ???? ... depending on J>0 (H) or J<0 (J) ? Something more
# fundamental possible ?
if J < 0
        # tdm magnitude <0,0|μ|0,1>
        Ja = states[1,:]' * dense(μ).data * states[2,:]
        # tdm magnitude <0,0|μ|1,0>
        Jb = states[1,:]' * dense(μ).data * states[3,:]
elseif J > 0
        # everything reversed
        Ja = states[1,:]' * dense(μ).data * states[3,:]
        Jb = states[1,:]' * dense(μ).data * states[2,:]
end

# plot
ax3 = subplot(334)
bar(0, real(Ja))
bar(1, real(Jb))
xticks([0, 1], ["μ_g-e1", "μ_g-e2"])
xlim([-1, 2]); ylim([0, 1.5 * maximum(real([Ja, Jb]))])
ylabel("tdm strength")

## calculate and plot expectation values of population
#=
tout, rhot = timeevolution.master(tlist, rho0, H, F)
nc = real(expect(j12 ⊗ j21, rhot))
na = real(expect(j21 ⊗ j12, rhot))
figure(figsize=(6,3));
plot(tout,nc)
plot(tout,na)
=#
#rho0.data[1,1] = 0.9
#rho0.data[2,2] = .0
#rho0.data[3,3] = .1

## calculated 1st order correlation function and associated spectrum
# STILL UNSURE ABOUT FREQUENCY AXIS/CONVENTION
corr = timecorrelations.correlation(tlist, rho0, H, L, μ12, μ12)

## zeropad corr and extrapolate tlist
zp = 0
corr = cmds.zeropad(corr,zp)
tnew, ~ = cmds.interpt(tlist,zp)
ω,spec = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true)

## add the results to previous figure
ax4 = subplot(336)
plot(tnew,real(corr),"gray")
title("System correlation function")
xlabel("time"); ylabel("Corr. func. ⟨ ... | ... ⟩")

ax5 = subplot(339)
plot(-[E₁, E₁], [0, .5])
plot(-[E₂, E₂], [0, .5])
plot(ω,spec,"gray")
plot(-[energies[2], energies[2]], [0, .5],"k--")
plot(-[energies[3], energies[3]], [0, .5],"k--")
title("System abs. spectrum")
xlabel("ω"); ylabel("Absorption")

ωC = 1
S = 0.1
# check out: physics.stackexchange.com/questions/474313/ohmic-spectral-density
function spectral_density(ω)
    if ω == 0.0
        return 0
    else
        return #complex(ω)^S / ωC^(S-1) * exp(-ω/ωC);
    end
end

γ1 = 10
function noise_power(ω)
    if ω == 0.0
        return 0
    else
        return 1 / π * γ1 * ωC ./ (ω.^2 .+ ωC^2)
    end
end

λ = 2.2; γ = 22.1;
#check: home.uchicago.edu/~tokmakoff/TDQMS/Notes/7._Fluctuations_3-09.pdf
function debye(ω)
    if ω == 0.0
        return 0
    else
        return 2 * λ * ω * γ ./ (ω.^2 .+ γ^2)
    end
end

### HOW TO GET RIGHT OPERATORS ??????????
### ????????????????????????????????????
a_ops = [Γ_a[1] * one(b_fock) ⊗ one(b_mon) ⊗ (j12*j21) , noise_power,
         Γ_a[2] * one(b_fock) ⊗ (j12*j21) ⊗ one(b_mon) , noise_power,
         Γ_a[3] * one(b_fock) ⊗ (j21*j12) ⊗ (j21*j12)  , noise_power,
         Γ_a[4] * (at*a)                                 , noise_power]

#a_ops = [*(Γ_a[i], a_ops[2*i-1]) for i in 1:length(Γ_a)]

R, ekets = timeevolution.bloch_redfield_tensor(H, [a_ops], J=L)
# different elements in R ... https://www.pnas.org/content/pnas/108/52/20908.full.pdf
# changing R[1,1] let's 2D signal disappear after a certain time T
#R.data[1,1] = -0.01
# however ρ[1,1] goes to zero ... shouldn't it go back to one ???

tout, rhot = timeevolution.master_bloch_redfield(tlist,μ12*rho0,R,H)
corr = expect(μ,rhot)

zp = 11
corr = cmds.zeropad(corr,zp)
tnew, ~ = cmds.interpt(tout,zp)

#subplot(336)
ax4.plot(tnew,real(corr),"g",linewidth=1,label="RF")

ω,spec = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true)
ax5.plot(ω,spec)

#subplot(236)
#ax5.plot(ω,spec,"g",linewidth=1,label="RF")

tout, rhot = timeevolution.master_bloch_redfield(tlist,μ12*rho0*μ12,R,H)

es1 = expect(one(b_fock) ⊗ (j21*j12)  ⊗ one(b_mon), rhot)
es2 = expect(one(b_fock) ⊗ one(b_mon) ⊗ (j21*j12),  rhot)
cav = expect(create(b_fock) * destroy(b_fock) ⊗ one(b_mon) ⊗ one(b_mon), rhot)
gs  = expect(rho0, rhot)

subplot(333)
plot(tout,real(es1),tout,real(es2),tout,cav,tout,gs)
legend(["ES1", "ES2", "cav", "GS"],loc="right")
tight_layout()

"""
figure(figsize=(6,4))
ww = [0:.1:100;]
plot(ww,debye(ww));
"""

method = "redfield"
if method == "redfield"
    F       = R
    use_sub = false
elseif method == "lindblad"
    F       = L
    use_sub = true
end

if calc_2d

        ## calculate (complex) 3rd order corr function (with T=0)
        zp = 11 # zeropad up to 2^zp

        ## calculate 2D spectra at
        T = [0,20] #fs

        out2d = Array{cmds.out2d}(undef, length(T))

        Threads.@threads for i = 1:length(T)
            out2d[i] = cmds.make2Dspectra(tlist,rho0,H,F,μ12,μ23,T[i],
                                                method;debug=true,use_sub=use_sub,zp=zp);
        end

        ## crop 2D data and increase dw
        out2d = [cmds.crop2d(out2d[i],1.3;w_max=3.5,step=1) for i = 1:length(T)]

        ## plot 2D spectra for each(?) T
        # what to plot
        rep="absorptive"
        scal="lin"

        ## make  subplot layout
        nplots = length(T);
        ncols = Int32(ceil(sqrt(nplots)));
        nrows = Int32(ceil(nplots / ncols));

        # determine maximum value in dataset out2d[:].full2d[:,:]
        maxi = maximum([maximum(real(out2d[i].full2d)) for i in 1:length(out2d)])

        # plot 2D spectra
        fig, ax = subplots(nrows,ncols,sharex=true,sharey=true,figsize=(ncols*3.2,nrows*3),squeeze=false)
        fig.suptitle(rep * " 2D spectrum (" * scal * ". scaling)")
        k = 0
        for i = 1:ncols
            for j = 1:nrows
                global k += 1
                if k > nplots
                    continue
                end
                sca(ax[k])
                ax[k].set_aspect = "equal"
                cmds.plot2d(out2d[k].ω,round.(out2d[k].full2d,digits=1);repr=rep,scaling=scal,norm=maxi)
                ax[k].set_title("2D spectrum at $(T[k]) fs")
            end
        end
        tight_layout()
        subplots_adjust(top=0.88)

        ω = out2d[1].ω
        ## plot additional things, like energy levels of states
        plot([E₁, E₁], [ω[1], ω[end]],"k--",linewidth=1,alpha=0.25)
        plot([E₂, E₂], [ω[1], ω[end]],"k--",linewidth=1,alpha=0.25)
        plot([ω[1], ω[end]],[E₁, E₁],"k--",linewidth=1,alpha=0.25)
        plot([ω[1], ω[end]],[E₂, E₂],"k--",linewidth=1,alpha=0.25)
        plot([energies[2:3], energies[2:3]], [ω[1], ω[end]],"g--",linewidth=1,alpha=0.25)
        plot([ω[1], ω[end]],[energies[2:3], energies[2:3]],"g--",linewidth=1,alpha=0.25)

        ## Save data
        cmds.save_2d([round.(real(out2d[i].full2d),digits=1) for i = 1:length(T)],T,fn)

        ## plot TA (summed 2D spectrum)
        figure()
        ta = [sum(real(out2d[i].full2d),dims=1) for i in 1:length(out2d)] ./ length(out2d)
        plot(ω,vcat(ta...)')
        plot(ω,zeros(size(ω)),linestyle = "dashed")
        xlabel("Energy/frequency")
        ylabel("Difference absorption")
        tight_layout()

end

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
# This program solves the Laplace-Beltrami PDE on a unit sphere.
# In strong form, the problem can be stated as:
# find u s.t.
#   -Δ_S(u) = f on \Omega
# Where Δ_S is the spherical Laplacian defined as Δ_S = ∇_S⋅(∇_S(u)),
# with ∇_S() and  ∇_S⋅ being the spherical gradient and spherical
# divergence, resp.

# Given a known analytical solution u(x), the program manufactures f(x)
# by applying -Δ_S to u, and then performs an error convergence analysis.
# For the discretization of the sphere, it uses the so-called cubed sphere
# mesh

function u(x)
  sin(π*x[1])*cos(π*x[2])*exp(x[3])
end
function uθϕ(θϕ)
  u(θϕ2xyz(θϕ))
end
function f(x)
  -(laplacian_unit_sphere(uθϕ)(xyz2θϕ(x)))
end
function Gridap.∇(::typeof(u))
  gradient_unit_sphere(uθϕ)
end

function solve_laplace_beltrami(model,order,degree,ls=BackslashSolver())
    V = FESpace(model,ReferenceFE(lagrangian,Float64,order); conformity=:H1)
    U = TrialFESpace(V)

    Ω = Triangulation(model)
    dΩ = Measure(Ω,degree)

    a(u,v) = ∫(∇(v)⋅∇(u))dΩ
    b(v)   = ∫(v*f)dΩ

    op = AffineFEOperator(a,b,U,V)
    fels = LinearFESolver(ls)
    uh = solve(fels,op)

    e=u-uh
    ∇e = ∇(uh)-∇(u)∘xyz2θϕ

    # H1-norm
    # sqrt(sum(∫( e*e + (∇e)⋅(∇e) )dΩ))
    sqrt(sum(∫( (∇e)⋅(∇e) )dΩ))#,uh
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
module Refactoring
##
import ..opt

collect_skip_one(ex::Expr) = any( ex.head .== [:call, :(=), :kw, :macrocall])
dont_collect(ex::Expr) = any( ex.head .== [:line, :using])

#function remove_type_annotation(ex::Expr)
#    ex.head == :(::) && return ex.args[1]
#    return :nothing
#end
#remove_type_annotation(s::Symbol) = s

#get the variable name from the rhs
function assigned_var(ex::Expr)

    if ex.head == :call #e.g. f(x, y) = x + y -> []
        return [ex.args[1]] #map(remove_type_annotation, ex.args[2:end])
    end

    if ex.head == :tuple # x, y = 1, 2
        out = Symbol[]
        for i=1:length(ex.args)
            push!(out, tuple_assigment(ex.args[i]))
            typeof(ex.args[i]) == Expr && ex.args[i].head == :(=) && return out
        end
        return out
    end
    typeof(ex.args[1]) == Symbol && return [ex.args[1]]

    ex.args[1]
end
assigned_var(s::Symbol) = [s]

tuple_assigment(s::Symbol) = s
tuple_assigment(ex::Expr) = ex.args[1]

is_assignement(ex::Expr) = ex.head == :(=)

function is_tuple_assigment(ex::Expr)
    if ex.head == :tuple
        for i=1:length(ex.args)
             typeof(ex.args[i]) == Expr && ex.args[i].head == :(=) && return true
        end
    end
    false
end

function arguments(ex::Expr, out, assigned)

    is_assignement(ex) && push!(assigned, assigned_var(ex.args[1])...)
    is_tuple_assigment(ex) && push!(assigned, assigned_var(ex)...)

    if collect_skip_one(ex)
        for i = 2:length(ex.args)
            arguments(ex.args[i], out, assigned)
        end
    elseif !dont_collect(ex)
        for i = 1:length(ex.args)
            arguments(ex.args[i], out, assigned)
        end
    end
end

"Filter out types and modules."
function filter_types(s::Symbol)
    try
        t = eval(Main, :(typeof($s)))
        t <: DataType && return false
        t <: Module && return false
        t <: Type && return false
    catch
    end
    true
end

function arguments(ex::Expr)
    out, assigned = Symbol[], Symbol[]
    arguments(ex, out, assigned)
    filter(filter_types, unique(out))
end
arguments(s::Symbol, out, assigned) = !any(s .== assigned) && push!(out, s)
arguments(s::Any, out, assigned) = nothing

#ex = quote
#    f(ex::Expr, b) = a
#    ex = 1
#end
#dump(ex)
#
#arguments(ex)

## Formatting

function indent_body(body)

    lines = split(body, "\n")
    ident_length = [length(l) - length(lstrip(l)) for l in lines]
    line_length = [length(l) for l in lines]

    min_ident = minimum( ident_length[line_length.>0] )

    tabw = opt("Editor", "tab_width")
    pre = " "^tabw

    lines = [pre*l[(1+min_ident):end]  for l in lines]

    join(lines, "\n")
end

function extract_method(body::AbstractString)

    ex  = try
        ex = Base.parse_input_line(body)
    catch err
        println(err)
        return ""
    end

    args = arguments(ex)
    if !isempty(args)
        sargs = string(args[1])
        for i=2:length(args)
            sargs = string(sargs, ", ", string(args[i]))
        end
    else
        sargs = ""
    end

    body = indent_body(body)

"
function ($sargs)
$body
end
"
end

##
end


text:
using Distributions
#This file contains the definition of the type VariableWeight, representing different
#style of weights that you can create
abstract VariableWeight
export VariableWeight
getWeight(w::VariableWeight,icomp::Integer)=error("getWeight method not defined for $(typeof(w))")

immutable LaplaceWeight
  w::Vector{Float64}
end
export LaplaceWeight
function LaplaceWeight(Ncomp)
  w = rand(Laplace(),Ncomp)
  scale!(w,inv(sqrt(sumabs2(w))))
  LaplaceWeight(w)
end
getWeight(w::LaplaceWeight,icomp::Integer,kw,d)=w.w[icomp]

"""
Type to represent a weight vector where the weights can shift during an event. the
field w contains the normal events and the field wd contains the disturbed events.
"""
immutable DisturbedLaplaceWeight
  w::Vector{Float64}
  wd::Vector{Float64}
end
export DisturbedLaplaceWeight
function DisturbedLaplaceWeight(Ncomp)
    w = rand(Laplace(),Ncomp)
    wd = rand(Laplace(),Ncomp)
    scale!(w,inv(sqrt(sumabs2(w))))
    #Make sure disturbance is orthogonal to original w
    wdort = wd - dot(wd,w)*w
    scale!(wdort,inv(sqrt(sumabs2(wdort))))
    DisturbedLaplaceWeight(w,wdort)
end
getWeight(w::DisturbedLaplaceWeight,icomp::Integer,kw,d)= d>0.0 ? kw * w.wd[icomp] + (1.0-kw) * w.w[icomp] : w.w[icomp]

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
module Masks

using NetCDF
using NCDatasets
using Interpolations

function CreateMask(settings)
    println("  Creating region masks...")

    # ----------------------------------------------------------------------------
    # Make a 0.5° mask of US coastal sea level to combine with CDS/CMEMS altimetry
    # 8 basins:
    # 1. NE Northeast Coast
    # 2. SE southeast coast
    # 3. EGOM Eastern Gulf of Mexico
    # 4. WGOM Gulf coast west
    # 5. SW Southwest coast
    # 6. NW Northwest coast
    # 7. PAC Pacific Islands (Hawaii)
    # 8. CAR Carribean Islands
    # 9. NAL North Alaska
    # 10. SAL South Alaska
    # ----------------------------------------------------------------------------
    fn_slm   = homedir()*"/Data/GRACE/JPL_mascon/mask.nc"

    # boundaries
    NE = [282.5 292;35.2  43.0]
    SE = [278.5 286; 25.0  35.2]
    EGOM = [ 270.5  278.5; 25.0  31.0]
    WGOM = [261.0  270.5; 25.0  31.0]
    SW = [234.0  242.0; 31.0  42.0]
    NW = [234.0  242.0; 42.0  48.0]
    PAC = [180.0  210.0; 18.0  30.0]
    CAR = [280.0  300.0; 14.0  25.0]
    NAL = [187.0  219.0; 51.0  72.0]
    SAL1 = [172.0  213.0; 47.0  60.0]
    SAL2 = [200.0  231.0; 54.75  62.5]

    # NE- Lat: [35.2 43]; Lon: [282.5 292]
    # SE- Lat: [25 35.2]; Lon: [278.5 286]
    # GC- Lat: [25 31]; Lon: [261 278.5]
    # WC- Lat: [31 47.5]; Lon: [234, 242]

    ϕ = reshape(ncread(fn_slm,"lon"),(:,1))
    θ = reshape(ncread(fn_slm,"lat"),(1,:))
    slm = convert.(Bool,1 .- ncread(fn_slm,"land"))
    depth = ncread(settings["fn_bathymetry"],"z")
    basin = ncread(settings["fn_basin_codes"],"basin")

    NE_mask = @. (ϕ>=NE[1,1]) & (ϕ<=NE[1,2]) & (θ>=NE[2,1]) & (θ<=NE[2,2]) & slm & (depth > -500)
    SE_mask = @. (ϕ>=SE[1,1]) & (ϕ<=SE[1,2]) & (θ>=SE[2,1]) & (θ<=SE[2,2]) & slm & (depth > -700)
    EGOM_mask = @. (ϕ>=EGOM[1,1]) & (ϕ<=EGOM[1,2]) & (θ>=EGOM[2,1]) & (θ<=EGOM[2,2]) & slm & (depth > -2000)
    WGOM_mask = @. (ϕ>=WGOM[1,1]) & (ϕ<=WGOM[1,2]) & (θ>=WGOM[2,1]) & (θ<=WGOM[2,2]) & slm & (depth > -2000)
    SW_mask = @. (ϕ>=SW[1,1]) & (ϕ<=SW[1,2]) & (θ>=SW[2,1]) & (θ<=SW[2,2]) & slm & (depth > -2500)
    NW_mask = @. (ϕ>=NW[1,1]) & (ϕ<=NW[1,2]) & (θ>=NW[2,1]) & (θ<=NW[2,2]) & slm & (depth > -2500)
    CAR_mask = @. (ϕ>=CAR[1,1]) & (ϕ<=CAR[1,2]) & (θ>=CAR[2,1]) & (θ<=CAR[2,2]) & slm & (depth > -2500)
    PAC_mask = @. (ϕ>=PAC[1,1]) & (ϕ<=PAC[1,2]) & (θ>=PAC[2,1]) & (θ<=PAC[2,2]) & slm & (depth > -5000)
    NAL_mask = @. (ϕ>=NAL[1,1]) & (ϕ<=NAL[1,2]) & (θ>=NAL[2,1]) & (θ<=NAL[2,2]) & slm & (depth > -500) & (basin !=1) & (basin !=0)

    SAL1_mask = @. (ϕ>=SAL1[1,1]) & (ϕ<=SAL1[1,2]) & (θ>=SAL1[2,1]) & (θ<=SAL1[2,2]) & slm & (depth > -500) & (basin ==1)
    SAL2_mask = @. (ϕ>=SAL2[1,1]) & (ϕ<=SAL2[1,2]) & (θ>=SAL2[2,1]) & (θ<=SAL2[2,2]) & slm & (depth > -500) & (basin ==1)
    SAL_mask = @. SAL1_mask | SAL2_mask
    USA_mask = @. NE_mask | SE_mask | EGOM_mask | WGOM_mask | SW_mask | NW_mask

    # Mask without land removed (for tide-gauge selection)
    NE_unc = @. (ϕ>=NE[1,1]) & (ϕ<=NE[1,2]) & (θ>=NE[2,1]) & (θ<=NE[2,2])
    SE_unc = @. (ϕ>=SE[1,1]) & (ϕ<=SE[1,2]) & (θ>=SE[2,1]) & (θ<=SE[2,2])
    EGOM_unc = @. (ϕ>=EGOM[1,1]) & (ϕ<=EGOM[1,2]) & (θ>=EGOM[2,1]) & (θ<=EGOM[2,2])
    WGOM_unc = @. (ϕ>=WGOM[1,1]) & (ϕ<=WGOM[1,2]) & (θ>=WGOM[2,1]) & (θ<=WGOM[2,2])
    SW_unc = @. (ϕ>=SW[1,1]) & (ϕ<=SW[1,2]) & (θ>=SW[2,1]) & (θ<=SW[2,2])
    NW_unc = @. (ϕ>=NW[1,1]) & (ϕ<=NW[1,2]) & (θ>=NW[2,1]) & (θ<=NW[2,2])
    CAR_unc = @. (ϕ>=CAR[1,1]) & (ϕ<=CAR[1,2]) & (θ>=CAR[2,1]) & (θ<=CAR[2,2])
    PAC_unc = @. (ϕ>=PAC[1,1]) & (ϕ<=PAC[1,2]) & (θ>=PAC[2,1]) & (θ<=PAC[2,2])
    NAL_unc = @. (ϕ>=NAL[1,1]) & (ϕ<=NAL[1,2]) & (θ>=NAL[2,1]) & (θ<=NAL[2,2])
    SAL1_unc = @. (ϕ>=SAL1[1,1]) & (ϕ<=SAL1[1,2]) & (θ>=SAL1[2,1]) & (θ<=SAL1[2,2])
    SAL2_unc = @. (ϕ>=SAL2[1,1]) & (ϕ<=SAL2[1,2]) & (θ>=SAL2[2,1]) & (θ<=SAL2[2,2])
    SAL_unc = @. SAL1_unc | SAL2_unc

    USA_unc = @. NE_unc | SE_unc | EGOM_unc | WGOM_unc | SW_unc | NW_unc

    # Make numerical mask for plot
    mask_num = zeros(Int16,size(NE_mask))
    mask_num[NE_mask] .= 1
    mask_num[SE_mask] .= 2
    mask_num[EGOM_mask] .= 3
    mask_num[WGOM_mask] .= 4
    mask_num[SW_mask] .= 5
    mask_num[NW_mask] .= 6
    mask_num[PAC_mask] .= 7
    mask_num[CAR_mask] .= 8
    mask_num[NAL_mask] .= 9
    mask_num[SAL_mask] .= 10

    # Save
    println("   Saving...")

    fh = Dataset(settings["fn_region_mask"],"c")
    defDim(fh,"lon", length(ϕ[:]))
    defDim(fh,"lat", length(θ[:]))
    defVar(fh,"lon",ϕ[:],("lon",),deflatelevel=5)
    defVar(fh,"lat",θ[:],("lat",),deflatelevel=5)
    defVar(fh,"NE_mask",convert.(UInt8,NE_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"SE_mask",convert.(UInt8,SE_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"EGOM_mask",convert.(UInt8,EGOM_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"WGOM_mask",convert.(UInt8,WGOM_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"SW_mask",convert.(UInt8,SW_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"NW_mask",convert.(UInt8,NW_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"PAC_mask",convert.(UInt8,PAC_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"CAR_mask",convert.(UInt8,CAR_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"NAL_mask",convert.(UInt8,NAL_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"SAL_mask",convert.(UInt8,SAL_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"USA_mask",convert.(UInt8,USA_mask),("lon","lat",),deflatelevel=5)
    defVar(fh,"mask_num",mask_num,("lon","lat",),deflatelevel=5)

    defVar(fh,"NE_unc",convert.(UInt8,NE_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"SE_unc",convert.(UInt8,SE_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"EGOM_unc",convert.(UInt8,EGOM_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"WGOM_unc",convert.(UInt8,WGOM_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"SW_unc",convert.(UInt8,SW_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"NW_unc",convert.(UInt8,NW_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"CAR_unc",convert.(UInt8,CAR_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"PAC_unc",convert.(UInt8,PAC_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"NAL_unc",convert.(UInt8,NAL_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"SAL_unc",convert.(UInt8,SAL_unc),("lon","lat",),deflatelevel=5)
    defVar(fh,"USA_unc",convert.(UInt8,USA_unc),("lon","lat",),deflatelevel=5)
    close(fh)
    println("  Creating region done")

    return nothing
end

function ReadMask(settings)
    mask = Dict()
    mask["ϕ"] = ncread(settings["fn_region_mask"],"lon")
    mask["θ"] = ncread(settings["fn_region_mask"],"lat")
    for region ∈ settings["regions"]
        mask[region] = convert.(Bool,ncread(settings["fn_region_mask"],region*"_mask"))
        mask[region*"_unc"] = convert.(Bool,ncread(settings["fn_region_mask"],region*"_unc"))
    end
    return mask
end

function RegridBasinCodes(settings)
    # --------------------------------------------------------------
    # Regrid the basin code map from Eric Leuliette/NOAA
    # Data provided by the NOAA Laboratory for Satellite Altimetry."
    # --------------------------------------------------------------
    fn_in = homedir()*"/Data/Basins/basin_codes.nc"
    # Read
    ncinfo(fn_in)
    θ_in = ncread(fn_in,"lat")
    ϕ_in = ncread(fn_in,"lon")
    B_in = 1.0f0 .* ncread(fn_in,"basin")

    # Interpolate
    ϕ = [0.25f0:0.5f0:359.75f0...]
    θ = [-89.75f0:0.5f0:89.75f0...]
    B = interpolate((ϕ_in,θ_in), B_in, Gridded(Constant()))(ϕ,θ)

    # Write
    fh = Dataset(settings["fn_basin_codes"],"c")
    defDim(fh,"lon", length(ϕ))
    defDim(fh,"lat", length(θ))
    defVar(fh,"lon",ϕ,("lon",),deflatelevel=5)
    defVar(fh,"lat",θ,("lat",),deflatelevel=5)
    defVar(fh,"basin",convert.(UInt8,B),("lon","lat",),deflatelevel=5)
    fh.attrib["title"] = ["BASIN5 - Global 5-minute Basin Grid regridded to 0.5 degree"]
    fh.attrib["data_source"] = ["Data provided by the NOAA Laboratory for Satellite Altimetry, regridded afterwards"]
    fh.attrib["author"] = ["Original data: Eric Leuliette NOAA, regridded by Thomas Frederikse NASA JPL/Caltech"]
    close(fh)
    return nothing
end

end


text:
# shows some data analysis using DataFrames
using CSV, DataFrames, GLM
dir = dirname(dirname(pathof(Econometrics)))
nerlove = CSV.read(dir*"/Examples/Data/nerlove.csv")
nerlove[:lnC] = log.(nerlove[:cost])
nerlove[:lnQ] = log.(nerlove[:output])
nerlove[:lnPL] = log.(nerlove[:labor])
nerlove[:lnPF] = log.(nerlove[:fuel])
nerlove[:lnPK] = log.(nerlove[:capital])
f = @formula(lnC ~ 1 + lnQ + lnPL + lnPF + lnPK)
lm(f, nerlove)

max_stars_repo_name: manuelmontesinos/Econometrics

max_issues_repo_name: manuelmontesinos/Econometrics

max_forks_repo_name: manuelmontesinos/Econometrics


text:
s₁ = Scenario(ξ₁ = 6., ξ₂ = 8., probability = 0.5)
s₂ = Scenario(ξ₁ = 4., ξ₂ = 4., probability = 0.5)

infeasible = StochasticProgram([s₁,s₂])

@first_stage infeasible = begin
    @variable(model, x₁ >= 0)
    @variable(model, x₂ >= 0)
    @objective(model, Min, 3*x₁ + 2*x₂)
end

@second_stage infeasible = begin
    @decision x₁ x₂
    @uncertain ξ₁ ξ₂
    @variable(model, 0.8*ξ₁ <= y₁ <= ξ₁)
    @variable(model, 0.8*ξ₂ <= y₂ <= ξ₂)
    @objective(model, Min, -15*y₁ - 12*y₂)
    @constraint(model, 3*y₁ + 2*y₂ <= x₁)
    @constraint(model, 2*y₁ + 5*y₂ <= x₂)
end

push!(problems, (infeasible,"Infeasible"))


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

# Wavelength range [nm]
λ = 200:1100
# Angle of incidence [degrees]
θ = [0.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.fusedsilicauv2(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:tauclorentz),
    LayerTMMO(emergent),
]

# Measured absolute transmittance
Texp = SpectraDB.HafniaSpectrum(beam.λ)
plot(Spectrum1D(), beam.λ, Texp)
gui()

## Optimization using Optim

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=true,
);

# One oscillator
seed = [vcat(175.0, # thickness
             [3.0, # high-frequency index
               2.5, 2.5, 2.37, 5.0], # osc 1
        ),
]

solOptim = FitTMMOptics(
    Transmittance(), seed, beam, Texp, layers;
    options=options, alg=SAMIN(), lb=0.2.*seed, ub=10.0.*seed,
)

plot(FitSpectrum(),
    solOptim.Beam.λ, solOptim.spectrumExp, solOptim.spectrumFit,
    xaxis=("Wavelength [nm]"), yaxis=("Transmittance"),
)
gui()

# julia> solOptim.optParams
# 1-element Array{Array{Float64,1},1}:
#  [185.57495457185868, 1.4538249665571539, 24.99997160640067, 6.801747171200166, 0.47400001243358125, 4.579667310612569]
#
# julia> solOptim.objfunMin
# 0.0004221112894584093

# Two oscillators
seed2 = [vcat(185.0, # thickness
              [2.0, # high-frequency permittivity
                25.0, 7.0, 0.5, 5.0, # osc 1
                10.0, 7.0, 0.5, 5.0, # osc 2
              ]
        ),
]

solOptim2 = FitTMMOptics(
    Transmittance(), seed2, beam, Texp, layers;
    options=options, alg=SAMIN(), lb=0.1.*seed2, ub=10.0.*seed2,
)

plot(FitSpectrum(),
    solOptim2.Beam.λ, solOptim2.spectrumExp, solOptim2.spectrumFit,
    xaxis=("Wavelength [nm]"), yaxis=("Transmittance"),
)
gui()

# julia> solOptim2.optParams
# 1-element Array{Array{Float64,1},1}:
#  [189.09445508077323, 0.39609259415335923, 158.43421436032554, 8.610532933615291, 3.4744910306832337, 5.015700019476466, 15.249954668687487, 0.7000170470355342, 0.06250726388751021, 1.0341375083174082]
#
# julia> solOptim2.objfunMin
# 0.0001506570743847816


text:
module MergeTypesModule

using CoordsModule: Coords
using GridModule: Grid,for_all
import Base.zero

export MergeTypes
export SimpleMergeTypes
export no_merge_mtype
export connection_merge_as_primary_flood_merge_as_primary
export connection_merge_as_primary_flood_merge_as_secondary
export connection_merge_as_primary_flood_merge_not_set
export connection_merge_as_primary_flood_merge_as_both
export connection_merge_as_secondary_flood_merge_as_primary
export connection_merge_as_secondary_flood_merge_as_secondary
export connection_merge_as_secondary_flood_merge_not_set
export connection_merge_as_secondary_flood_merge_as_both
export connection_merge_not_set_flood_merge_as_primary
export connection_merge_not_set_flood_merge_as_secondary
export connection_merge_not_set_flood_merge_as_both
export connection_merge_as_both_flood_merge_as_primary
export connection_merge_as_both_flood_merge_as_secondary
export connection_merge_as_both_flood_merge_not_set
export connection_merge_as_both_flood_merge_as_both
export convert_to_simple_merge_type_connect
export convert_to_simple_merge_type_flood
export no_merge
export null_mtype
export primary_merge
export secondary_merge
export double_merge
export null_simple_mtype
export zero

@enum MergeTypes begin
no_merge_mtype = 0
connection_merge_as_primary_flood_merge_as_primary = 1
connection_merge_as_primary_flood_merge_as_secondary = 2
connection_merge_as_primary_flood_merge_not_set = 3
connection_merge_as_primary_flood_merge_as_both = 4
connection_merge_as_secondary_flood_merge_as_primary = 5
connection_merge_as_secondary_flood_merge_as_secondary = 6
connection_merge_as_secondary_flood_merge_not_set = 7
connection_merge_as_secondary_flood_merge_as_both = 8
connection_merge_not_set_flood_merge_as_primary = 9
connection_merge_not_set_flood_merge_as_secondary = 10
connection_merge_not_set_flood_merge_as_both = 11
connection_merge_as_both_flood_merge_as_primary = 12
connection_merge_as_both_flood_merge_as_secondary = 13
connection_merge_as_both_flood_merge_not_set = 14
connection_merge_as_both_flood_merge_as_both = 15
null_mtype = 16
end

function zero(::Type{MergeTypes})
  return no_merge_mtype::MergeTypes
end

@enum SimpleMergeTypes begin
no_merge = 0
primary_merge = 1
secondary_merge = 2
double_merge = 3
null_simple_mtype = 4
end

const convert_to_simple_merge_type_connect =
  SimpleMergeTypes[no_merge primary_merge primary_merge #=
                =# primary_merge primary_merge secondary_merge #=
                =# secondary_merge secondary_merge secondary_merge #=
                =# no_merge no_merge no_merge double_merge #=
                =# double_merge double_merge double_merge #=
                =# null_simple_mtype]
const convert_to_simple_merge_type_flood =
  SimpleMergeTypes[no_merge primary_merge secondary_merge #=
                =# no_merge double_merge primary_merge #=
                =# secondary_merge no_merge double_merge #=
                =# primary_merge secondary_merge double_merge #=
                =# primary_merge secondary_merge no_merge #=
                =# double_merge null_simple_mtype]

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.SymmetricLoop
using GeometricExamples

const Δt     = 5.
const ntime  = 10000
const nloop  = 2000
const nsave  = 10
const nplot  = 5
const run_id = "poincare_1st_symmetric_vprk_dt5"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_symmetric_poincare_invariant_1st.jl")


text:
function test_line_search()
    fcn(x) = x[1]^2 + x[2]^2

    x0 = [1.0, 2]
    step_direction = [-1.0, -1.0]
    f0 = fcn(x0)
    step_size = opti.line_search(fcn, x0, step_direction, f0)
    println(step_size)
    x1 = x0 + step_direction * step_size
    println(x1)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Reproduce
using FileIO, JLD2

function arg_parse_settings(as::ArgParseSettings = ArgParseSettings())

    @add_arg_table as begin
        "--opt1"
        arg_type=Int64
        "--opt2"
        arg_type=Int64
        "--steps"
        arg_type=Int64
    end
    return as
end

function main_experiment(args::Vector{String}, saveloc::String="default_save_loc")

    arg_settings = arg_parse_settings()
    parsed = parse_args(args, arg_settings)
    create_info!(parsed, saveloc)
    j = 0
    if parsed["opt1"] == 2
        throw("Oh No!!!")
    end
    # sleep(0.1*(parsed["opt1"]^4))
#    for i in 1:parsed["opt1"]^3
#        j += i
#    end
    @save joinpath(parsed["_SAVE"], "data.jld2") args

    return j
end


text:
function LVGNA_max_rank_experiments(testing=true)
    path = "../data/sparse_tensors/LVGNA/"
    files = readdir(path)
    if testing
        files = files[6:7]
        files =  ["worm_PHY1.ssten","worm_Y2H1.ssten"]
        alphas = [.5]
        betas = [1e1]
        result_path = all_results_path*"TAME_iterate_max_rank/test/"
    else
        alphas = [.5,1.0]
        betas = [0.0,1e1,1e2,1e3]
        result_path = all_results_path*"TAME_iterate_max_rank/"
    end

    if !isdir(result_path)
        mkpath(result_path)
    end

    X = distributed_pairwise_alignment(files,path,method=LowRankTAME_M(),iter=15;
                                       alphas,betas,profile=true,no_matching=true)
    filename=filter(x-> !isspace(x),"LVGNA_pairwise_alignments_alpha:$(alphas)_beta:$(betas)_iter:15_noMatching_profiled_results.json")
    save_results(result_path*filename,X)

    #=  # Additional drivers for running TAME
    X = distributed_pairwise_alignment(files,path,method=TAME_M(),iter=15,
                                       alphas,betas,profile=true,no_matching=true)
                                    profile=true,no_matching=true)
    =#
    #TODO: add in saving routines and filenames
end

function LVGNA_shift_experiment(testing=true)

    tensor_files = readdir(sparse_tensors_path)

    if testing
        samples = 5
        tol=1e-8
        max_iter=5
        tensor_files = tensor_files[6:7]
        result_path = all_results_path*"TAME_iterate_max_rank/test/"
    else

        samples = 3000
        tol=1e-16
        max_iter=30

        result_path = all_results_path*"TAME_iterate_max_rank/"
    end
    alphas = [1.0]
    shift_percentage=1.0

    if !isdir(result_path)
        mkpath(result_path)
    end

    foundMaxEigs = []
    for file in tensor_files
        A = load_ThirdOrderSymTensor(sparse_tensors_path*file)
        push!(foundMaxEigs,LambdaTAME.SSHOPM_sample(A,samples,0.0,max_iter,tol)[2])
    end

    rank_profiles = []
    for i = 1:length(tensor_files)
        A = load_ThirdOrderSymTensor(sparse_tensors_path*tensor_files[i])

        for j = i+1:length(tensor_files)
            B = load_ThirdOrderSymTensor(sparse_tensors_path*tensor_files[j])
            β = foundMaxEigs[i]*foundMaxEigs[j]*shift_percentage

            output = align_tensors_profiled(A,B;method = LowRankTAME_M(),alphas=alphas,betas=[β],iter=max_iter,tol=tol,no_matching=true)
            push!(rank_profiles,(tensor_files[i],tensor_files[j],foundMaxEigs[i],foundMaxEigs[j],output.profile))
        end
    end

    filename=filter(x-> !isspace(x),"LVGNAMaxEigShiftRanks_alphas:$(alphas)_iter:$(max_iter)_SSHOPMSamples:$(samples)_shiftPercentage:$(shift_percentage)_tol:$(tol)_results.json")
    save_results(result_path*filename,rank_profiles)
end

function random_graph_experiments(testing=true)

    if testing
        n_sizes = [100]
        alphas = [.5]
        betas = [100.0]
        trials = 2
        result_path = all_results_path*"TAME_iterate_max_rank/test/"
    else
        n_sizes = [100, 500, 1000, 2000, 5000, 10000]
        trials = 50
        alphas = [.5,1.0]
        betas = [0.0,1.0,10.0,100.0]
        result_path = all_results_path*"TAME_iterate_max_rank/"
    end

    p_remove = [.05]  #ER param
    edge_inclusion_p = [.5]
    step_percentage = [.25]

    if !isdir(result_path)
        mkpath(result_path)
    end

    X = distributed_random_trials(trials,ErdosRenyiNoise(),true;method=LowRankTAME_M(),graph=RandomGeometric(),iter=15,
                                    degreedist=LogNormal(log(5),1), p_remove,n_sizes,profile=true,no_matching=true,
                                    postProcessing=noPostProcessing(),tol=1e-12,alphas,betas)
    filename=filter(x-> !isspace(x),"LRTAME_noMatch_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_p:$(p_remove)_noiseModel:ER_trials:$(trials)_MaxRankResults.json")
    save_results(result_path*filename,X)

    X = distributed_random_trials(trials,DuplicationNoise(),true;method=LowRankTAME_M(),graph=RandomGeometric(),iter=15,
                                    degreedist=LogNormal(log(5),1), edge_inclusion_p, step_percentage,n_sizes,
                                    profile=true,no_matching=true, postProcessing=noPostProcessing(),
                                    tol=1e-12,alphas,betas)
    filename=filter(x-> !isspace(x),"LRTAME_noMatch_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_p:$(edge_inclusion_p)_noiseModel:Duplication_sp:$(step_percentage)_trials:$(trials)_MaxRankResults.json")
    save_results(result_path*filename,X)

    #=    # Additional drivers for running TAME
    X = distributed_random_trials(trials,ErdosRenyiNoise(),true;method=TAME_M(),graph=RandomGeometric(),iter=15,
                                    degreedist=LogNormal(log(5),1), p_remove,n_sizes,profile=true,no_matching=true,
                                    postProcessing=noPostProcessing(),tol=1e-12,alphas,betas)
    filename=filter(x-> !isspace(x),"TAME_noMatch_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_p:$(p_remove)_noiseModel:ER_trials:$(trials)_MaxRankResults.json")
    save_results(result_path*filename,X)

    X = distributed_random_trials(trials,DuplicationNoise(),true;method=TAME_M(),graph=RandomGeometric(),iter=15,
                                    degreedist=LogNormal(log(5),1), edge_inclusion_p, step_percentage,n_sizes,
                                    profile=true,no_matching=true, postProcessing=noPostProcessing(),
                                    tol=1e-12,alphas,betas)
    filename=filter(x-> !isspace(x),"TAME_noMatch_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_p:$(edge_inclusion_p)_noiseModel:Duplication_sp:$(step_percentage)_trials:$(trials)_MaxRankResults.json")
    save_results(result_path*filename,X)
    =#

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME