text:
"""
    explicit midpoint integrator
"""
function midpoint(model, x, u, w, h)
    x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
end

"""
    implicit midpoint integrator
"""
function midpoint_implicit(model, x⁺, x, u, w, h)
    x⁺ - (x + h * f(model, 0.5 * (x + x⁺), u, w))
end

"""
    explicit RK3 integrator
"""
function rk3(model, z, u, w, h)
    k1 = k2 = k3 = zero(z)
    k1 = h * f(model, z, u, w)
    k2 = h * f(model, z + 0.5 * k1, u, w)
    k3 = h * f(model, z - k1 + 2.0 * k2, u, w)
    z + (k1 + 4.0 * k2 + k3) / 6.0
end

"""
    discrete explicit dynamics
"""
function fd(model, x, u, w, h, t)
    midpoint(model, x, u, w, h)
end

"""
    discrete implicit dynamics
"""
function fd(model, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, h)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
	println("	PLOT")
	# PLOT
	plot(collect(1:10),rand(10), lw=1, lc="blue", fmt=:ps, marker="circle", markeredgecolor=0, size=0.2, markerfacecolor="red", title="Bla Bla", xlabel="Spoons", ylabel="Forks", savefig="lixo")
	plot(mat2ds(GMT.fakedata(6,6), x=:ny, color=[:red, :green, :blue, :yellow], ls=:dashdot, multi=true), leg=true, label="Bla")
	plot("",hcat(collect(1:10)[:],rand(10,1)))
	plot!("",hcat(collect(1:10)[:],rand(10,1)), Vd=dbg2)
	plot(hcat(collect(1:10)[:],rand(10,1)), Vd=dbg2)
	plot!(hcat(collect(1:10)[:],rand(10,1)), Vd=dbg2)
	plot!(collect(1:10),rand(10), fmt="ps")
	plot(0.5,0.5, R="0/1/0/1", Vd=dbg2)
	plot!(0.5,0.5, R="0/1/0/1", Vd=dbg2)
	plot(1:10,rand(10), S=(symb=:c,size=7,unit=:point), color=:rainbow, zcolor=rand(10))
	plot(1:10,rand(10)*3, S="c7p", color=:rainbow, Vd=dbg2)
	plot(1:10,rand(10)*3, S="c7p", color=:rainbow, zcolor=rand(10)*3)
	plot(1:2pi, rand(6), xaxis=(pi=1,), Vd=dbg2)
	plot(1:2pi, rand(6), xaxis=(pi=(1,2),), Vd=dbg2)
	plot([5 5], region=(0,10,0,10), frame=(annot=:a, ticks=:a, grid=5), figsize=10, symbol=:p, markerline=0.5, fill=:lightblue, E=(Y=[2 3 6 9],pen=1,cap="10p"), Vd=dbg2);
	plot(rand(10,4), S=:c, ms=0.2, markeredgecolor=:red, ml=2, Vd=dbg2)
	plot(rand(10,4), S=:c, ms=0.2, marker=:star, ml=2, W=1, Vd=dbg2)
	plot([0.0 0; 1.1 1], theme=(name=:dark, bg_color="gray"), lc=:white, Vd=dbg2)
	#plot([0.0 0; 1.1 1], theme=(name=:modern,), Vd=dbg2)
	plot(1:4, rand(4,4), theme=(name=:none, save=true), leg=true)	# Resets default conf and delete the theme_jl file
	@test startswith(plot!([1 1], marker=(:r, [2 3]), Vd=dbg2), "psxy  -R -J -Sr")
	xy = [0. 0.; 1 1; 2 1; 0 0];
	plot(xy, region=(-1,3,-1,2), clip=(xy, N=true), frame=(axes=:WEsn, grid=0.2, fill=:darkgray), Vd=dbg2);
	plot(x=[-1.0, -1.0], y=[-1.0, 1.0], pen="1p,black+ve0.2c")
	@test occursin(" -Sr", plot!([1 1], marker=(:r, [2 3]), Vd=dbg2))
	@test_throws ErrorException("Wrong number of extra columns for marker (r). Got 3 but expected 2") plot!([1 1], marker=(:r, [2 3 4]), Vd=dbg2)
	@test_throws ErrorException("Unknown graphics file extension (.ppf)") plot(rand(5,2), savefig="la.ppf")
	@test startswith(plot!([1 1], marker=(:Web, [2 3], (inner=5, arc=30,radial=45, pen=(2,:red))), Vd=dbg2), "psxy  -R -J -SW/5+a30+r45+p2,red")
	@test startswith(plot!([1 1], marker=(Web=true, inner=5, arc=30,radial=45, pen=(2,:red)), Vd=dbg2), "psxy  -R -J -SW/5+a30+r45+p2,red")
	@test startswith(plot!([1 1], marker="W/5+a30", Vd=dbg2), "psxy  -R -J -SW/5+a30")
	@test startswith(plot!([1 1], marker=:Web, Vd=dbg2), "psxy  -R -J -SW")
	@test startswith(plot!([1 1], marker=:W, Vd=dbg2), "psxy  -R -J -SW")
	@test startswith(plot([5 5], marker=(:E, 500), Vd=dbg2), "psxy  -JX" * GMT.def_fig_size * " -Baf -BWSen -R4.74/5.26/4.74/5.26 -SE-500")
	@test startswith(plot(region=(0,10,0,10), marker=(letter="blaBla", size="16p"), Vd=dbg2), "psxy  -R0/10/0/10 -JX" * GMT.def_fig_size * " -Baf -BWSen -Sl16p+tblaBla")
	@test startswith(plot([5 5], region=(0,10,0,10), marker=(bar=true, size=0.5, base=0,), Vd=dbg2), "psxy  -R0/10/0/10 -JX" * GMT.def_fig_size * " -Baf -BWSen -Sb0.5+b0")
	@test startswith(plot([5 5], region=(0,10,0,10), marker=(custom=:sun, size=0.5), Vd=dbg2), "psxy  -R0/10/0/10 -JX" * GMT.def_fig_size * " -Baf -BWSen -Sksun/0.5")
	plot([5 5 0 45], region=(0,10,0,10), marker=(pie=true, arc=15, radial=30), Vd=dbg2)
	plot([0.5 1 1.75 5 85], region=(0,5,0,5), figsize=12, marker=(matang=true, arrow=(length=0.75, start=true, stop=true, half=:right)), ml=(0.5,:red), fill=:blue, Vd=dbg2)
	plot([2.5 2.5], region=(0,4,0,4), figsize=12, marker=(:matang, [2 50 350], (length=0.75, start=true, stop=true, half=:right)), ml=(0.5,:red), fill=:blue, Vd=dbg2)
	plot([1 1], limits=(0,6,0,6), figsize=7, marker=:circle, ms=0.5, error_bars=(x=:x, cline=true), Vd=dbg2)	# Warning line
	plot(rand(5,3), region=[0,1,0,1])
	plot("polyg_hole.gmt")
	plot(x -> x^3 - 2x^2 + 3x - 1, xlim=(-5,5), Vd=dbg2)
	plot(x -> x^3 - 2x^2 + 3x - 1, -10:10, Vd=dbg2)
	plot!(x -> x^3 - 2x^2 + 3x - 1, 10, Vd=dbg2)
	plot!(x -> x^3 - 2x^2 + 3x - 1, Vd=dbg2)
	plot!(x -> cos(x) * x, y -> sin(y) * y, linspace(0,2pi,100), Vd=dbg2)
	plot(1:5, rand(5), ls="Line&Bla Bla&", lc=:red, Vd=dbg2)
	plot(1:5, rand(5), ls="Line&Circ", lc=:red, Vd=dbg2)
	plot(rand(5,2), ls="DotCirc", Vd=dbg2)
	plot(rand(5,2), ls="DotCirc#", Vd=dbg2)
	plot(rand(5,2), ls="DashDot", Vd=dbg2)
	plot(rand(5,2), ls="-.", Vd=dbg2)
	@test_throws ErrorException("Bad line style. Options are (for example) [Line|DashDot|Dash|Dot]Circ") plot(rand(5,2), ls="Dat")
	bar!(x -> x^3 - 2x^2 + 3x - 1, Vd=dbg2)
	lines!(x -> x^3 - 2x^2 + 3x - 1, Vd=dbg2)
	lines!(x -> cos(x) * x, y -> sin(y) * y, linspace(0,2pi,100), Vd=dbg2)
	x = LinRange(0,2π,50);
	lines(x, sin.(x), ls="FrontSlip", legend=true, figname="lixo", theme=("A2"))
	lines(x, sin.(x), ls="lineCirc", legend=(label="sin(x)", box=:none), lc=:red, figname="lixo")
	lines(x, sin.(x), ls="FrontCircLeft", legend=(label="sin(x)"), mc=:red, figname="lixo")
	lines(x, sin.(x), ls="linediamond", theme=(:A0atg), Vd=dbg2)
	lines(x, sin.(x), ls="linediamond", theme=(:A2atg), Vd=dbg2)
	lines(x, sin.(x), ls="linediamond", theme=(:A2agITGraph), Vd=dbg2)
	theme("classic", noticks=true, graygrid=true)
	theme("modern")		# Must reset otherwise many tests would fail due to different configs
	scatter!(x -> x^3 - 2x^2 + 3x - 1, Vd=dbg2)
	scatter!(x -> cos(x) * x, y -> sin(y) * y, linspace(0,2pi,100), Vd=dbg2)
	hlines!([0.2, 0.6], pen=(1, :red))
	vlines!([0.2, 0.6], pen=(1, :red))

	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, Vd=dbg2);

	println("	PLOT3D")
	plot3d(rand(5,3), marker=:cube)
	plot3d!(rand(5,3), marker=:cube, Vd=dbg2)
	plot3d("", rand(5,3), Vd=dbg2)
	plot3d!("", rand(5,3), Vd=dbg2)
	plot3d(1:10, rand(10), rand(10), Vd=dbg2)
	plot3d!(1:10, rand(10), rand(10), Vd=dbg2)
	plot3d!(x -> sin(x), y -> cos(y), 0:pi/50:10pi, Vd=dbg2)
	plot3d!(x -> sin(x)*cos(10x), y -> sin(y)*sin(10y), z -> cos(z), 0:pi/100:pi, Vd=dbg2)
	scatter3!(x -> sin(x), y -> cos(y), 0:pi/50:10pi, Vd=dbg2)
	scatter3!(x -> sin(x)*cos(10x), y -> sin(y)*sin(10y), z -> cos(z), 0:pi/100:pi, Vd=dbg2)

	println("	ARROWS")
	arrows([0 8.2 0 6], R="-2/4/0/9", arrow=(len=2,stop=1,shape=0.5,fill=:red), J=:X14, B=:a, pen="6p")
	arrows([0 8.2 0 6], R="-2/4/0/9", arrow=(len=2,start=:arrow,stop=:tail,shape=0.5), J=:X14, B=:a, pen="6p")
	arrows!([0 8.2 0 6], R="-2/4/0/9", arrow=(len=2,start=:arrow,shape=0.5), pen="6p", Vd=dbg2)
	arrows!("", [0 8.2 0 6], R="-2/4/0/9", arrow=(len=2,start=:arrow,shape=:V, pen=(1,:red)), pen="6p", Vd=dbg2)
	@test_throws ErrorException("Bad data type for the 'shape' option") arrows!("", [0 8.2 0 6], R="-2/4/0/9", arrow=(len=2,justify=:center,middle=:reverse,shape='1'), pen="6p", Vd=dbg2)
	@test occursin("-SvB4p/18p/7.5p", arrows([1 0 45 4], R="0/6/-1/1", J="x2.5", pen=(1,:blue), arrow4=(align=:mid, head=(arrowwidth="4p", headlength="18p", headwidth="7.5p"), double=true), Vd=dbg2))
	@test occursin("-SvB4p/18p/7.5p", arrows([1 0 45 4], R="0/6/-1/1", J="x2.5", lw=1, arrow4=(align=:mid, head=("4p","18p", "7.5p"), double=true), Vd=dbg2))
	@test occursin("-Svs4p/18p/7.5p", arrows([1 0 45 4], R="0/6/-1/1", J="x2.5", lw=1, arrow4=(align=:pt, head="4p/18p/7.5p"), Vd=dbg2))

	println("	LINES")
	lines([0 0; 10 20], R="-2/12/-2/22", J="M2.5", W=1, G=:red, decorated=(dist=(1,0.25), symbol=:box))
	lines([-50 40; 50 -40],  R="-60/60/-50/50", J="X10", W=0.25, B=:af, box_pos="+p0.5", leg_pos=(offset=0.5,), leg=:TL)
	lines!([-50 40; 50 -40], R="-60/60/-50/50", W=1, offset="0.5i/0.25i", Vd=dbg2)
	lines(1:10,rand(10), W=0.25, Vd=dbg2)
	lines!(1:10,rand(10), W=0.25, Vd=dbg2)
	lines!("", rand(10), W=0.25, Vd=dbg2)
	xy = gmt("gmtmath -T0/180/1 T SIND 4.5 ADD");
	lines(xy, R="-5/185/-0.1/6", J="X6i/9i", B=:af, W=(1,:red), decorated=(dist=(2.5,0.25), symbol=:star, symbolsize=1, pen=(0.5,:green), fill=:blue, dec2=1))
	D = histogram(randn(100), I=:o, T=0.1);
	@test_throws ErrorException("Something went wrong when calling the module. GMT error number = 72") histogram(randn(100), I=:o, V=:q, W=0.1);
	lines(D, steps=(x=true,), close=(bot=true,))
	x = GMT.linspace(0, 2pi);  y = cos.(x)*0.5;
	r = lines(x,y, limits=(0,6.0,-1,0.7), figsize=(40,8), pen=(lw=2,lc=:sienna), decorated=(quoted=true, n_labels=1, const_label="ai ai", font=60, curved=true, fill=:blue, pen=(0.5,:red)), par=(:PS_MEDIA, :A1), axis=(fill=220,),Vd=dbg2);
	@test startswith(r, "psxy  -R0/6.0/-1/0.7 -JX40/8 -Baf -BWSen+g220 --PS_MEDIA=A1 -Sqn1:+f60+l\"ai ai\"+v+p0.5,red -W2,sienna")

	println("	SCATTER")
	sizevec = [s for s = 1:10] ./ 10;
	scatter(1:10, 1:10, markersize = sizevec, aspect=:equal, B=:a, marker=:square, fill=:green)
	scatter(rand(10), leg=:bottomrigh, fill=:red)	# leg wrong on purpose
	scatter(1:10,rand(10)*3, S="c7p", color=:rainbow, zcolor=rand(10)*3, show=1, Vd=dbg2)
	scatter(rand(50),rand(50), markersize=rand(50), zcolor=rand(50), aspect=:equal, alpha=50, Vd=dbg2)
	scatter(1:10, rand(10), fill=:red, B=:a)
	scatter!(1:10, rand(10), fill=:red, B=:a, Vd=dbg2)
	scatter(1:10, Vd=dbg2)
	scatter!(1:10, Vd=dbg2)
	scatter(rand(5,5))
	scatter!(rand(5,5), Vd=dbg2)
	scatter("",rand(5,5), Vd=dbg2)
	scatter!("",rand(5,5), Vd=dbg2)
	scatter3(rand(5,5,3))
	scatter3!(rand(5,5,3), Vd=dbg2)
	scatter3("", rand(5,5,3), Vd=dbg2)
	scatter3!("", rand(5,5,3), Vd=dbg2)
	scatter3(1:10, rand(10), rand(10), fill=:red, B=:a, Vd=dbg2)
	scatter3!(1:10, rand(10), rand(10), Vd=dbg2)

	println("	BARPLOT")
	bar(sort(randn(10)), G=0, B=:a)
	bar(rand(20),bar=(width=0.5,), Vd=dbg2)
	bar!(rand(20),bar=(width=0.5,), Vd=dbg2)
	bar(1:20,  rand(20),bar=(width=0.5,), Vd=dbg2)
	bar!(1:20, rand(20),bar=(width=0.5,), Vd=dbg2)
	bar(rand(20),hbar=(width=0.5,unit=:c, base=9), Vd=dbg2)
	bar(rand(20),bar="0.5c+b9",  Vd=dbg2)
	bar(rand(20),hbar="0.5c+b9",  Vd=dbg2)
	bar(rand(10), xaxis=(custom=(pos=1:5,type="A"),), Vd=dbg2)
	bar(rand(10), axis=(custom=(pos=1:5,label=[:a :b :c :d :e]),), Vd=dbg2)
	@test_throws ErrorException("Number of labels in custom annotations must be the same as the 'pos' element") bar(rand(10), frame=(custom=(pos=1:5,label=[:a :b :c :d]),), axis=:noannot, Vd=dbg2)
	bar((1,2,3), Vd=dbg2)
	bar((1,2,3), (1,2,3), Vd=dbg2)
	bar!((1,2,3), Vd=dbg2)
	bar!((1,2,3), (1,2,3), Vd=dbg2)
	bar([3 31], C=:lightblue, Vd=dbg2)
	bar("", [3 31], C=:lightblue, Vd=dbg2)
	bar!("", [3 31], C=:lightblue, frame=:noannot, Vd=dbg2)
	men_means, men_std = (20, 35, 30, 35, 27), (2, 3, 4, 1, 2);
	x = collect(1:length(men_means));
	bar(x.-0.35/2, collect(men_means), width=0.35, color=:lightblue, limits=(0.5,5.5,0,40), frame=:none, error_bars=(y=men_std,), Vd=dbg2)
	bar([0. 1 2 3; 1 2 3 4], error_bars=(y=[0.1 0.2 0.33; 0.2 0.3 0.4],), fillalpha=[0.3 0.5 0.7], Vd=dbg2)
	bar([0. 1 2 3; 1 2 3 4], fill=(1,2,3), Vd=dbg2)
	bar([0. 1 2 3; 1 2 3 4], stack=1, Vd=dbg2)
	bar(1:5, [20 25; 35 32; 30 34; 35 20; 27 25], fill=["lightblue", "brown"], xaxis=(ticks=(:G1, :G2, :G3, "G4"),), Vd=dbg2)
	bar(1:5, [20 25; 35 32; 30 34; 35 20; 27 25], fill=["lightblue", "brown"], xticks=(:G1, :G2, :G3, :G4), zticks=(:Z1,:Z2), Vd=dbg2)
	bar(1:3,[-5 -15 20; 17 10 21; 10 5 15], stacked=1, Vd=dbg2)
	bar([0. 1 2 3; 1 2 3 4], fill=("red", "green", "blue"), Vd=dbg2)
	bar(rand(15), color=:rainbow, Y=3)
	T = mat2ds([1.0 0.446143; 2.0 0.581746; 3.0 0.268978], text=[" "; " "; " "]);
	bar(T, color=:rainbow, figsize=(14,8), title="Colored bars", Vd=dbg2)
	T = mat2ds([1.0 0.446143 0; 2.0 0.581746 0; 3.0 0.268978 0], text=[" "; " "; " "]);
	bar(T, color=:rainbow, figsize=(14,8), mz=[3 2 1], Vd=dbg2)
	bar(1:5, (20, 35, 30, 35, 27), width=0.35, color=:lightblue,limits=(0.5,5.5,0,40),E=(y=(2,3,4,1,2),), Vd=dbg2)
	bar([0. 1 2 3; 1 2 3 4], stack=true, hbar=true, fill=["red", "green", "blue"])
	bar(0:3, [2.9, 5, 2.2, 1], fill=["red", "green", "blue", "orange"], Vd=dbg2)
	D = mat2ds([0 0],["aa"]);
	sprint(print, D);

	println("	BAR3")
	G = gmt("grdmath -R-15/15/-15/15 -I1 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
	gmtwrite("lixo.grd", G)
	bar3(G, lw=:thinnest)
	bar3("lixo.grd", grd=true, lw=:thinnest, Vd=dbg2)
	bar3!(G, lw=:thinnest, Vd=dbg2)
	bar3!("", G, lw=:thinnest, Vd=dbg2)
	bar3(G, lw=:thinnest, bar=(width=0.085,), Vd=dbg2)
	bar3(G, lw=:thinnest, nbands=3, Vd=dbg2)
	bar3(G, region="-15/15/-15/15/-2/2", lw=:thinnest, noshade=1, Vd=dbg2)
	bar3(rand(4,4), Vd=dbg2)
	D = grd2xyz(G);
	bar3(D, width=0.01, Nbands=3, Vd=dbg2)
	@test_throws ErrorException("BAR3: When first arg is a name, must also state its type. e.g. grd=true or dataset=true") bar3("lixo.grd")
	gmtwrite("lixo.gmt", D)
	@test_throws ErrorException("BAR3: When NOT providing *width* data must contain at least 5 columns.") bar3("lixo.gmt", dataset=true)

ext: jl


text:
const BASICS_FNS = [:add, :mul, :sub, :addinv, :lt, :lte, :gt, :gte, :eq]
const TPTR_TO_BINF = Dict(:add => "+", :mul => "*", :sub => "-", :div => "/")
const ftype = Dict(8 => :float, 16 => :float, 32 => :double)

doc"""
  The skeleton of a binary operation, with adlibbable values
"""
function binop_fn(n::Integer ,es::Integer, op::Symbol, status::Bool)
  opfn = ops[op][status + 1]

  #build up the assignment expression.  Start with the basics:  decide if the
  #main is going to be a C operator or a function call.
  fexpr = haskey(TPTR_TO_BINF, op) ? "$(to_f(n,es,:a)) $(TPTR_TO_BINF[op]) $(to_f(n,es,:b))" :
    "$(op)($(to_f(n,es,:a)), $(to_f(n,es,:b)))"

  #generate the return statement, which depends on whether or not we're status or jump.
  ret = status ? "res->udata = pres.udata; return errno" : "return pres"

"""
extern "C" $(opfn(n,es)) {
  $(p(n,es)) pres;

  $(ftype[n]) fres;
  fres = $fexpr;

  pres = $(to_p(n,es,:fres));
  $ret;
}
"""
end

const basics_code = Dict(
  :add => (n, es) -> """
  $(binop_fn(n,es,:add,true))
  $(binop_fn(n,es,:add,false))
  """,

  :sub => (n, es) -> """
  $(binop_fn(n,es,:sub,true))
  $(binop_fn(n,es,:sub,false))
  """,

  :mul => (n, es) -> """
  $(binop_fn(n,es,:mul,true))
  $(binop_fn(n,es,:mul,false))
  """,

  :addinv => (n, es) -> """
  extern "C" $(ops[:addinv][1](n,es)) {
    $(p(n,es)) res;
    res.udata = -(a.udata);
    return res;
  }

  extern "C" $(ops[:addinv][2](n,es)) {
    res->udata = -(a.udata);
    return 0;
  }
  """,

  :lt => (n, es) -> """
  extern "C" $(ops[:lt][2](n,es)) {
    if ((a.udata == P$(n)INF) || (b.udata == P$(n)INF)) { return !POSIT_ENV.nanmode; }
    return (a.sdata < b.sdata);
  }
  """,

  :lte => (n, es) -> """
  extern "C" $(ops[:lte][2](n,es)) {
    if ((a.udata == P$(n)INF) || (b.udata == P$(n)INF)) { return !POSIT_ENV.nanmode; }
    return (a.sdata <= b.sdata);
  }
  """,

  :gt => (n, es) -> """
  extern "C" $(ops[:gt][2](n,es)) {
    if ((a.udata == P$(n)INF) || (b.udata == P$(n)INF)) { return !POSIT_ENV.nanmode; }
    return (a.sdata > b.sdata);
  }
  """,

  :gte => (n, es) -> """
  extern "C" $(ops[:gte][2](n,es)) {
    if ((a.udata == P$(n)INF) || (b.udata == P$(n)INF)) { return !POSIT_ENV.nanmode; }
    return (a.sdata >= b.sdata);
  }
  """,

  :eq => (n, es) -> """
  extern "C" $(ops[:eq][2](n,es)) {
    if ((a.udata == P$(n)INF) || (b.udata == P$(n)INF)) { return !POSIT_ENV.nanmode; }
    return (a.sdata == b.sdata);
  }
  """)

function generate_posit_basics_cpp(io, posit_defs)
  #generates "posit.h" based on the posit_definitions
  write(io,"""
  #include "include/posit.h"
  #include "include/posit_conv.h"
  #include <errno.h>
  #include <stdbool.h>

  environment_t POSIT_ENV = {false, false};

  extern "C" void set_nanmode(bool nanmode){ POSIT_ENV.nanmode = nanmode; }
  extern "C" void set_underflow(bool underflows){ POSIT_ENV.underflows = underflows; }
  extern "C" bool get_nanmode(){ return POSIT_ENV.nanmode; }
  extern "C" bool get_underflow(){ return POSIT_ENV.underflows; }

  """)

  for fn in BASICS_FNS
    for n in sort(collect(keys(posit_defs)))

      write(io, """

      ///////////////////////////////////////////////////////////////
      //  posit_$(n) section, variable ES adapters for $(fn)
      ///////////////////////////////////////////////////////////////
      """)

      for es in posit_defs[n]
        write(io, basics_code[fn](n, es))
        write(io, "\n")
      end
    end
    write(io,"\n")
  end
end


text:
using PkgSRA, Plots, LaTeXStrings
using Turing, AxisArrays, DataFrames, LinearAlgebra, DataStructures
using BSON: @load
pyplot()

#####
##### Load the data
#####
this_dat_name = DAT_FOLDERNAME*"dat_neuron_"

## Constant input
fname = this_dat_name*"controlled_model.bson";
@load fname chain_ctr

## Spikes with varying input
fname = this_dat_name*"controlled_model2.bson";
@load fname chain_ctr2

#####
##### Produce the plots
#####
plot_opt = Dict(:titlefontsize=>24,
        :yticks=>false, :fontfamily=>:serif,
        :legendfontsize=>14, :guidefontsize=>16,
        :xtickfontsize=>16)

# Plot
key = "all_coef[1]"
plot1 = density(Array(chain_ctr[key]), lw=3; plot_opt...)
        density!(Array(chain_ctr2[key]), lw=3)
        vline!([5], lw=3, color=:black, legend=false)
        title!("Coefficient for v")

key = "all_coef[5]"
plot2 = density(Array(chain_ctr[key]), lw=3; plot_opt...)
        density!(Array(chain_ctr2[key]), lw=3)
        vline!([180], lw=3, color=:black, legend=false)
        title!("Constant term")

key = "all_coef[7]"
plot3 = density(Array(chain_ctr[key]), lw=3, label="Constant Input"; plot_opt...)
        density!(Array(chain_ctr2[key]), label="Variable Input", lw=3)
        vline!([0.04], lw=3, color=:black,
                label="True value", legend=:bottomright)
        xlims!(0.03, 0.05)
        title!(L"v^2")

# Create the layout and plot
my_layout = @layout [p1; p2; p3];
    p_final = plot(plot1, plot2, plot3, layout = my_layout)

# Save
fname = FIGURE_FOLDERNAME * "fig_neuron_parameters.png";
savefig(p_final, fname)


text:
include("../solvers/azsolver.jl")
include("../solvers/lowranksolver.jl")
include("../solvers/randomized.jl")
include("../solvers/smoothsolver.jl")

# Dispatch on solver style
solver(ap::ApproximationProblem; options...) = solver(SolverStyle(ap), ap; options...)

# add the linear system as a third argument
function solver(solverstyle::SolverStyle, ap::ApproximationProblem; options...)
    A = discretization(ap; options...)
    solver(solverstyle, ap, A; options...)
end

solver(::InverseStyle, ap::ApproximationProblem, A::AbstractOperator; options...) = inv(A)
solver(::DirectStyle, ap::ApproximationProblem, A::AbstractOperator; options...) =
	directsolver(A; options...)
solver(::IterativeStyle, ap::ApproximationProblem, A::AbstractOperator; options...) =
    iterativesolver(A; options...)
solver(::DualStyle, ap::ApproximationProblem, A::AbstractOperator; options...) =
	dualdiscretization(ap; options...)'

solver(solverstyle::ProductSolverStyle, ap::ApproximationProblem, A::AbstractOperator; samplingstyle=SamplingStyle(ap), options...) =
    solver(solverstyle, samplingstyle, ap, A; options...)
function solver(solverstyle::ProductSolverStyle, samplingstyle::ProductSamplingStyle, ap::ApproximationProblem, A::AbstractOperator; options...)
	S = sampling_operator(ap)
    solvere, sse, ape, Ae, Se = components(solverstyle), components(samplingstyle), factors(ap), components(A), factors(S)
    @assert length(solvere) == length(sse) == length(ape) == length(Ae) == length(Se)
    TensorProductOperator(
		map( (solversi, ssi, api, Ai, Si)->solver(solversi, api, Ai; S=Si, samplingstyle=ssi, options...), solvere, sse, ape, Ae, Se)...
	)
end

solver(style::AZStyle, ap::ApproximationProblem, A::AbstractOperator; options...) =
    solver(style, ap, A, AZ_Zt(ap; options...); options...)

function solver(::AZStyle, ap::ApproximationProblem, A::AbstractOperator, Zt::AbstractOperator;
            B=nothing, smallcoefficients=false, smallcoefficients_atol=NaN, smallcoefficients_rtol=NaN, verbose=false, options...)
    if smallcoefficients
        w = BasisFunctions.quadweights(sampling_grid(ap; options...), measure(ap))
        normF = abs(sqrt(sum(w .* B.^2)))
        if !isnan(smallcoefficients_rtol)
            verbose && println("Change smallcoefficients relative tolerance to absolute tolerance rtol*||f||")
            smallcoefficients_atol = smallcoefficients_rtol*normF
            smallcoefficients_rtol = NaN
        end

        AZSolver(A, Zt; smallcoefficients=smallcoefficients, smallcoefficients_rtol=smallcoefficients_rtol,
                smallcoefficients_atol=smallcoefficients_atol, verbose=verbose, options...)
    else
        AZSolver(A, Zt; B=B, verbose=verbose, options...)
    end
end

solver(style::AZSmoothStyle, ap::ApproximationProblem, A::AbstractOperator; options...) = solver(style, ap, A, AZ_Zt(ap; options...); options...)
solver(::AZSmoothStyle, ap::ApproximationProblem, A::AbstractOperator, Zt; options...) =
    solver(AZStyle(), ap, A, Zt; options..., weightedAZ=true, AZ_Cweight=sobolevAZweight(A; options...))


text:
module Simulation

using Sockets
using JSON

include("abstracts.jl")

# Example:
# print_trace(stacktrace(), "run")
function print_trace(trace, stop_at::String)
    println("------- trace ---------")
    for tr in trace
        if string(tr.func) == stop_at
            break
        end
        println(tr)
    end
    println("-----------------------")
end

include("../../data/model/model.jl")
include("stimulus/stimulus.jl")
include("delay.jl")
include("simulate_neuron.jl")
include("run.jl")

end # Simulation --------------------------------------------------


text:
"""
Refine Float64 to BigFloat precision.
"""
function iterate!(A::AbstractMatrix{T}, ilo::Integer, ihi::Integer, Q::AbstractM) where {T<:BigFloat}

  # approximation of submatrix eigenvalues in Float64 precision
  r = ilo:ihi
  ev = seigendiag(smalleig(Float64.(view(A, r, r)))[1], 1, ihi-ilo+1)
  ihi = refineprecision!(A, ilo, ihi, Q, ev)
  ilo, ihi
end

"""
Refine all eigenvalues of a list close to estimation value
"""
function refineprecision!(A::AbstractMatrix{T}, ilo::Integer, ihi::Integer, Q::AbstractM, ev::AbstractVector) where {T<:Union{Real,Complex}}

  realcase = T <: Real
  evr = reverse(ev)
  k = length(ev)
  while k >= 1
    evk = ev[k]
    evk2 = ! realcase || isreal(evk) ? T(evk) : Complex{T}(evk)

    ihi0 = ihi
    ihi = refineprecision!(A, ilo, ihi, Q, evk)
    if ihi0 - ihi >= 2 && k > 1 && imag(evk) == -imag(ev[k-1])
      k -= 1
    end
    k -= 1
  end
  ihi
end

"""
Refine one eigenvalue close to estimation value
"""
function refineprecision!(A::AbstractMatrix{T}, ilo::Integer, ihi::Integer, Q::AbstractM, ev::Number) where {T<:Union{Real,Complex}}

  ih2 = ! (T <: Real) || isreal(ev) ? ihi : ihi - 1
  r = ih2:ihi
  while !converged(A, ilo, ihi)
    transform_Hess!(A, ilo, ihi, Q, [ev], ihi, 0)
    ev = eig2(view(A, r, r))[end]
  end
  ihi = ih2 - 1
  if ihi >= ilo
    A[ihi+1,ihi] = 0
  end
  ihi
end

"""
Detect convergence of last eigenvalue block of hessenberg matrix.
"""
function converged(A::AbstractMatrix{T}, ilo::Integer, ihi::Integer) where {T<:Union{Real,Complex}}
  if ilo >= ihi
    return true
  end
  if deflation_criterion(A[ihi,ihi-1], A[ihi-1,ihi-1], A[ihi,ihi])
    return true
  end
  if ! (T <: Real) || discriminant(A, ihi-1, ihi) > 0
    return false
  end
  if ilo >= ihi - 1
    return true
  end
  deflation_criterion(A[ihi-1,ihi-2], A[ihi-2,ihi-2], A[ihi-1,ihi-1]) # TODO
end

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
using .Hidato       # Note that the . here means to look locally for the module rather than in the libraries

const numbrix1 = """
 0  0  0  0  0  0  0  0  0
 0  0 46 45  0 55 74  0  0
 0 38  0  0 43  0  0 78  0
 0 35  0  0  0  0  0 71  0
 0  0 33  0  0  0 59  0  0
 0 17  0  0  0  0  0 67  0
 0 18  0  0 11  0  0 64  0
 0  0 24 21  0  1  2  0  0
 0  0  0  0  0  0  0  0  0 """

const numbrix2 = """
 0  0  0  0  0  0  0  0  0
 0 11 12 15 18 21 62 61  0
 0  6  0  0  0  0  0 60  0
 0 33  0  0  0  0  0 57  0
 0 32  0  0  0  0  0 56  0
 0 37  0  1  0  0  0 73  0
 0 38  0  0  0  0  0 72  0
 0 43 44 47 48 51 76 77  0
 0  0  0  0  0  0  0  0  0 """

const numbrixmoves = [[-1, 0], [0, -1], [0, 1], [1, 0]]

board, maxmoves, fixed, starts = hidatoconfigure(numbrix1)
printboard(board, " 0 ", "   ")
hidatosolve(board, maxmoves, numbrixmoves, fixed, starts[1][1], starts[1][2], 1)
printboard(board)

board, maxmoves, fixed, starts = hidatoconfigure(numbrix2)
printboard(board, " 0 ", "   ")
hidatosolve(board, maxmoves, numbrixmoves, fixed, starts[1][1], starts[1][2], 1)
printboard(board)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"Throw an error if the pointer to int variable indicates an error."
function _failsafe(usage_error::Ref{Cint})
    # follows https://github.com/JuliaGeo/GDAL.jl/blob/017bf6b8492dcd2186ced297076b283c3591d798/src/error.jl#L31-L37
    if usage_error[] != 0 && GDAL.getlasterrortype() in GDAL.throw_class
        throw(GDAL.GDALError())
    end
end

"""
List various information about a GDAL supported raster dataset.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdalinfo utility.

### Returns
String corresponding to the information about the raster dataset.
"""
function gdalinfo(dataset::Dataset, options = String[])
    options = GDAL.infooptionsnew(options, C_NULL)
    result = GDAL.info(dataset.ptr, options)
    GDAL.infooptionsfree(options)
    return result
end

"""
Convert raster data between different formats.

### Parameters
* **dataset**: The dataset to be translated.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdal_translate utility.

### Returns
The output dataset.
"""
function unsafe_gdaltranslate(
        dataset::Dataset,
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.translateoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.translate(dest, dataset.ptr, options, usage_error)
    GDAL.translateoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Image reprojection and warping function.

### Parameters
* **datasets**: The list of input datasets.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdalwarp utility.

### Returns
The output dataset.
"""
function unsafe_gdalwarp(
        datasets::Vector{Dataset},
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.warpappoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.warp(dest, Ptr{GDAL.GDALDatasetH}(C_NULL), length(datasets),
        [ds.ptr for ds in datasets], options, usage_error)
    GDAL.warpappoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Convert vector data between file formats.

### Parameters
* **datasets**: The list of input datasets (only 1 supported currently).
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the ogr2ogr utility.

### Returns
The output dataset.
"""
function unsafe_gdalvectortranslate(
        datasets::Vector{Dataset},
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.vectortranslateoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.vectortranslate(dest, Ptr{GDAL.GDALDatasetH}(C_NULL),
        length(datasets), [ds.ptr for ds in datasets], options, usage_error)
    GDAL.vectortranslateoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Tools to analyze and visualize DEMs.

### Parameters
* **dataset**: The source dataset.
* **pszProcessing**: the processing to apply (one of "hillshade", "slope",
    "aspect", "color-relief", "TRI", "TPI", "Roughness").
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdaldem utility.

# Keyword Arguments
* **colorfile**: color file (mandatory for "color-relief" processing,
    should be empty otherwise).

### Returns
The output dataset.
"""
function unsafe_gdaldem(
        dataset::Dataset,
        processing::String,
        options = String[];
        dest = "/vsimem/tmp",
        colorfile = C_NULL
    )
    if processing == "color-relief"
        @assert colorfile != C_NULL
    end
    options = GDAL.demprocessingoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.demprocessing(dest, dataset.ptr, processing, colorfile,
        options, usage_error)
    GDAL.demprocessingoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Convert nearly black/white borders to exact value.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the nearblack utility.

### Returns
The output dataset.
"""
function unsafe_gdalnearblack(
        dataset::Dataset,
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.nearblackoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.nearblack(dest, Ptr{GDAL.GDALDatasetH}(C_NULL), dataset.ptr,
        options, usage_error)
    GDAL.nearblackoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Create a raster from the scattered data.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdal_grid utility.

### Returns
The output dataset.
"""
function unsafe_gdalgrid(
        dataset::Dataset,
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.gridoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.grid(dest, dataset.ptr, options, usage_error)
    GDAL.gridoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Burn vector geometries into a raster.

### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdal_rasterize utility.

### Returns
The output dataset.
"""
function unsafe_gdalrasterize(
        dataset::Dataset,
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.rasterizeoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.rasterize(dest, Ptr{GDAL.GDALDatasetH}(C_NULL), dataset.ptr,
        options, usage_error)
    GDAL.rasterizeoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end

"""
Build a VRT from a list of datasets.

### Parameters
* **datasets**: The list of input datasets.
* **options**: List of options (potentially including filename and open
    options). The accepted options are the ones of the gdalbuildvrt utility.

### Returns
The output dataset.
"""
function unsafe_gdalbuildvrt(
        datasets::Vector{Dataset},
        options = String[];
        dest = "/vsimem/tmp"
    )
    options = GDAL.buildvrtoptionsnew(options, C_NULL)
    usage_error = Ref{Cint}()
    result = GDAL.buildvrt(dest, length(datasets), [ds.ptr for ds in datasets],
        C_NULL, options, usage_error)
    GDAL.buildvrtoptionsfree(options)
    _failsafe(usage_error)
    return Dataset(result)
end


text:
using Test
using TropicalTensors
using CUDA, CuYao
using ForwardDiff
using Random

CUDA.allowscalar(false)

@testset "cuda" begin
    T = Tropical{Float64}
    @test ones(T, 3) == Array(CUDA.ones(T, 3))
    @test zeros(T, 3) == Array(CUDA.zeros(T, 3))
end

@testset "spinglass" begin
    for i=1:10
        lt = SquareLattice(10, 10)
        sg = Spinglass(lt, randn(Float32, 180), zeros(Float32, 100))
        res = solve(sg; usecuda=true)
        res0 = solve(sg; usecuda=false)
        @test res.n == res0.n
    end
end

@testset "masked square" begin
    for i=1:10
        lt = rand_maskedsquare(5, 5, 0.7)
        sg = Spinglass(lt, randn(Float32, length(sgbonds(lt))), zeros(Float32, length(lt)))
        res = solve(sg; usecuda=true)
        res0 = solve(sg; usecuda=false)
        @test res.n == res0.n
    end
end

@testset "test Chimera" begin
    for i=1:10
        lt = ChimeraLattice(3, 3)
        sg = Spinglass(lt, randn(Float32, 12*4 + 9*16), zeros(Float32, 9*8))
        res = solve(sg; usecuda=true)
        res0 = solve(sg; usecuda=false)
        @test res.n == res0.n
    end
end

@testset "forwarddiff-gpu" begin
    L = 5
    Js = ones(Float32, 2L*(L-1))
    gs = ForwardDiff.gradient(x->solve(SquareLattice(L, L), x, zeros(eltype(x), L^2); usecuda=true).n, Js)
    @test gs ≈ ones(Float32, 2L*(L-1))
end

@testset "counting tropical matmul" begin
    a = CountingTropical{Float64}.(randn(5, 7))
    b = CountingTropical{Float64}.(randn(7, 5))
    @test all(a * b .≈ Array(CuArray(a) * CuArray(b)))
end


text:
module SrtmKgb26
  # KURUCZ
  sfluxref = [
   #       129.462, 15*_ZERO_ ]
    29.0079,  28.4088,     20.3099,  13.0283
    ,  11.8619,  9.95840,     6.68696,  5.38987
    ,  3.49829, 0.407693,    0.299027, 0.236827
    , 0.188502, 0.163489, 4.64335e-02, 2.72662e-03]

  #     Rayleigh extinction coefficient at all v
  rayl = [1.21263e-06,1.43428e-06,1.67677e-06,1.93255e-06
    , 2.19177e-06,2.44195e-06,2.66926e-06,2.85990e-06
    , 3.00380e-06,3.06996e-06,3.08184e-06,3.09172e-06
    , 3.09938e-06,3.10456e-06,3.10727e-06,3.10818e-06]
end


text:
struct Cancelled <: Exception end

cancellation_token() = Promise{Cancelled}()

const CancellationTokenType = typeof(cancellation_token())

const DissolveHandle = let
    send, _ = Reagents.channel()
    handle = Reagents.dissolve(send)
    typeof(handle)
end

struct CancelScope
    token::CancellationTokenType
    dissolvehandle::Union{DissolveHandle,Nothing}
end

# TODO: store `CancelScope`?
@contextvar CANCELLATION_TOKEN::Union{CancellationTokenType,Nothing} = nothing

function Julio.cancel!(token::CancellationTokenType)
    tryputting(token)(Cancelled())
end

function withopen(f, scope::CancelScope)
    try
        with_context(f, CANCELLATION_TOKEN => scope.token)
    catch err
        if iscancelled(err) && is_token_cancelled(scope.token) && !Julio.iscancelled()
            return
        end
        rethrow()
    end
end

# Unfortunately, other `Base.open` methods uses `Function`...
Base.open(f::Function, scope::CancelScope) = withopen(f, scope)

# const CancellationTokenRefType = typeof(Reagents.Ref{Promise{Cancelled}}())

# @contextvar CANCELLATION_TOKEN::CancellationTokenRefType = CancellationTokenRefType()

function waiting_cancel()
    token = CANCELLATION_TOKEN[]
    if token === nothing
        return Return(Block())
    else
        if istracing()
            cancelled = PostCommit(_ -> @trace(label = :cancelled))
            # cancelled = PostCommit(_ -> (yield(); @trace(label = :cancelled)))
        else
            cancelled = Identity()
        end
        return _fetching(token) ⨟ cancelled
    end
end

function cancellable_react!(reagent::Reagents.Reagent, x = nothing)
    ans = (waiting_cancel() | reagent)(x)
    if ans isa Cancelled
        throw(ans)
    end
    return ans
end

@inline function handle_result(@nospecialize(ans))
    if ans isa EventResult
        return ans.f()
    elseif ans isa Error
        throw(ans.value)
    elseif ans === nothing
        return nothing
    else
        return something(ans)
    end
end

Julio.sync(reagent::Reagents.Reagent, x = nothing) =
    handle_result(cancellable_react!(reagent, x))

# TODO: formalize output format better (should we use a sum type always?)
# TODO: unify this with `handle_result`
asresult(f) =
    Map() do ans
        if ans isa Error
            return ans
        else
            x = something(ans, Some(nothing))
            return EventResult(() -> f(x))
        end
    end

lowerevent(reagent::Reagents.Reagent) = reagent
lowerevent((reagent, f)::Pair{<:Reagents.Reagent}) = reagent ⨟ asresult(f)
lowerevent((fargs, f)::Pair{<:Tuple}) = event(fargs...) ⨟ asresult(f)
lowerevent(fargs::Tuple) = event(fargs...)

function Julio.select(event1, events...)
    reagent = mapfoldl(lowerevent, |, (event1, events...))
    return handle_result(cancellable_react!(reagent))
end

#=
struct TryFailed end

function Julio.try(f, args...)
    y = Julio.sync(event(f, args...) | Return(TryFailed()))
    return !(y isa TryFailed)
end
=#

apply(f, args...) = Julio.sync(lowerevent((f, args...)))

function Julio.withtimeout(f::F, seconds::Real) where {F}
    token = cancellation_token()
    with_context(CANCELLATION_TOKEN => token) do
        timer = Timer(seconds) do _
            Julio.cancel!(token)
        end
        try
            return Some(f())
        catch err
            err isa Cancelled && return nothing
            rethrow()
        finally
            close(timer)
        end
    end
end

event(::typeof(Julio.sleep), seconds::Real) = event(sleep, seconds)
function event(::typeof(sleep), seconds::Real)
    promise = Julio.Promise{Nothing}()
    Timer(seconds) do _
        put_nocancel!(promise, nothing)
    end
    return fetching(promise)
end

Julio.sleep(seconds) = apply(sleep, seconds)

Julio.shield(f) = with_context(f, CANCELLATION_TOKEN => nothing)

function Julio.oncancel(f, args...; kwargs...)
    token = @something(CANCELLATION_TOKEN[], return nothing)
    reagent = fetching(token) ⨟ PostCommit() do _
        f(args...; kwargs...)
    end
    return Reagents.dissolve(reagent; once = true)
end

Julio.cancel!(handle::DissolveHandle) = Reagents.clear!(handle)


text:
"""
    ' CIq(x::AbstractVector)'

Calculate median and 95% confidence interval around the median
"""

function CIq(x)
    m = median(x)
    (m - quantile(x,0.025), quantile(x,0.975) - m)
    # (quantile(x,[0.025,0.975])...,)
end

"""
    'wilcoxon(df::AbstractDataFrame,x::Symbol, y::Symbol, by = :MouseID; f = mean)'

x and y are paired vector to compare effect of condtion in y vs baseline in x
collect SignedRankTest analysis in a standardize table to plot
"""

function wilcoxon(x::AbstractVector{<:Real})
    t = SignedRankTest(x)
    dd = DataFrame()
    dd[!,:Median] = [t.median]
    dd[!,:CI] = [(t.median - confint(t)[1], confint(t)[2] - t.median)]
    dd[!,:P] = [pvalue(t)]
    dd[!,:Vals] = [t.vals for i in 1:1]
    dd
end

wilcoxon(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = wilcoxon(y-x)

function wilcoxon(df::AbstractDataFrame,x::Symbol, by = :MouseID; f = mean)
    df1 = combine(x => f => :Mean,groupby(df,by))
    check = SwitchingAnalysis.complete_vals(df1)
    if any(x -> !x, check)
        println( "got NaN values and dropped them")
        dd = dropnan(df1)
    else
        dd = df1
    end
    wilcoxon(dd[:,:Mean])
end

function wilcoxon(df::AbstractDataFrame,x::Symbol, y::Symbol, by = :MouseID; f = mean)
    df1 = combine([x,y] => (a,b) -> (Mean_x = f(a), Mean_y = f(b)), groupby(df,by))
    check = SwitchingAnalysis.complete_vals(df1)
    if any(x -> !x, check)
        println( "got NaN values and dropped them")
        dd = dropnan(df1)
    else
        dd = df1
    end
    wilcoxon(dd[:,:Mean_x], dd[:,:Mean_y])
end


text:
function _get_nb_edges_to_add(nv::Int, nb_edges_to_add::Union{Int, Float64})
    if nb_edges_to_add < 0 || (nb_edges_to_add > 1 && !isinteger(nb_edges_to_add))
        throw(DomainError("`nb_edges_to_add` must be a positive integer or a float between 0 and 1, got $nb_edges_to_add"))
    end
    isinteger(nb_edges_to_add) && return nb_edges_to_add
    return Int(floor(nb_edges_to_add * nv))
end

function _get_edges_distance(g::AbstractGraph, distance::Int, nb_edges_to_add_stop=Inf)
    distance < 0 && throw(DomainError("`distance` must be a positive integer, got $distance."))
    edge_array = []
    for v in randperm(nv(g))
        append!(edge_array, [(v, n[1]) for n in neighborhood_dists(g, v, distance) if n[2] == distance && v < n[1]])
        if length(edge_array) >= nb_edges_to_add_stop
            return edge_array
        end
    end
    return edge_array
end

function add_edges_distance!(g::AbstractGraph, nb_edges_to_add::Union{Int, Float64}; distance::Int, seed=MersenneTwister(42))
    Random.seed!(seed)
    nb_edges_to_add = _get_nb_edges_to_add(nv(g), nb_edges_to_add)
    edge_array = _get_edges_distance(g, distance, nb_edges_to_add)
    nb_edges_to_add_added = 0
    for e in shuffle(edge_array)
        nb_edges_to_add_added += add_edge!(g, e[1], e[2])
        if nb_edges_to_add_added >= nb_edges_to_add
            return nb_edges_to_add_added
        end
    end
    return nb_edges_to_add_added
end

function _get_edges_random(g::AbstractGraph, nb_edges_to_add::Int)
    edges_array = shuffle(collect(edges(Graphs.complement(g))))
    try
        return edges_array[begin:nb_edges_to_add]
    catch BoundsError
        return edges_array
    end
end

function add_edges_random!(g::AbstractGraph, nb_edges_to_add::Union{Int, Float64}; seed=42)
    Random.seed!(seed)
    nb_edges_to_add = _get_nb_edges_to_add(nv(g), nb_edges_to_add)
    edges_array = _get_edges_random(g, nb_edges_to_add)
    for e in edges_array
        add_edge!(g, e.src, e.dst)
    end
    return length(edges_array)
end

const add_edges_functions = Dict(
    "distance" => add_edges_distance!,
    "random" => add_edges_random!
)

function add_edges!(g::AbstractGraph, nb_edges_to_add::Union{Int, Float64}, how::AbstractString; kwargs...)
    add_edges_functions[how](g, nb_edges_to_add; kwargs...)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using GLPlot;GLPlot.init()
using Plots;glvisualize()
using GLVisualize, FileIO, Colors, Images

imfolder = joinpath(homedir(), "Desktop", "animal_imgs")
images = map(readdir(imfolder)) do impath
    map(RGBA{U8}, restrict(load(joinpath(imfolder, impath)))).data
end;
filter!(x->size(x)==size(first(images)), images)
mean_colors = map(images) do img
    mean(img)
end
x = map(comp1, mean_colors); y=map(comp2, mean_colors); z=map(comp3, mean_colors);
p1=scatter(
    x,y,z, markercolor=mean_colors,
    shape=:circle, markerstrokewidth=1, markerstrokecolor="white",
    hover=images, ms=12
)
p2=scatter(x,y,z, markerstrokecolor=mean_colors, shape=images, ms=15)
plot(p1, p2)

using Plots;glvisualize()
using GLVisualize, FileIO, Colors, Images
using GeometryTypes, GLAbstraction

x = ["(｡◕‿◕｡)", "◔ ⌣ ◔","(づ｡◕‿‿◕｡)づ", "┬──┬ ノ( ゜-゜ノ)", "(╯°□°）╯︵ ┻━┻", "¯\\_(ツ)_/¯"]
y = [-4, -3, -2, -1, 6, 0]
scatter(x, y)
t = bounce(linspace(-1.0f0,1f0, 20))
translation = map(t) do t
    rotationmatrix_y(1f0)*rotationmatrix_z(deg2rad(90f0)) * translationmatrix(Vec3f0(3,3,t))
end
hover = [
    "Plotting\nsomething!",
    "Something\nscientific?",
    "For once?",
    visualize(loadasset("cat.obj"), model=translation),
    "Noo! That\ncat again!",
    "😸$(Char(0x1F63B))",
]
scatter(x, y, m=(0.8, :diamond, 20), hover=hover)

title!("ℕ ⊆ ℕ₀ ⊂ ℤ ⊂ 2H₂ + O₂ ⇌ 2H₂O")
xaxis!(" ⠍⠊⠣ ახლავე გაიაროთ все вещи, вы можете построить")
yaxis!("Anger")


text:
using HDF5
using Base.Test

@testset "readremote" begin

    # check that we can read the official HDF5 example files
    urlbase = "https://support.hdfgroup.org/ftp/HDF5/examples/files/exbyapi/"

    fcmp = [0 1 2 3 4 5 6;
        2 1.66667 2.4 3.28571 4.22222 5.18182 6.15385;
        4 2.33333 2.8 3.57143 4.44444 5.36364 6.30769;
        6 3 3.2 3.85714 4.66667 5.54545 6.46154]'
    icmp = [0 -1 -2 -3 -4 -5 -6;
        0 0 0 0 0 0 0;
        0 1 2 3 4 5 6;
        0 2 4 6 8 10 12]'
    const SOLID, LIQUID, GAS, PLASMA = 0, 1, 2, 3
    ecmp = [SOLID SOLID SOLID SOLID SOLID SOLID SOLID;
            SOLID LIQUID GAS PLASMA SOLID LIQUID GAS;
            SOLID GAS SOLID GAS SOLID GAS SOLID;
            SOLID PLASMA GAS LIQUID SOLID PLASMA GAS]'
    scmp = ["Parting", "is such", "sweet", "sorrow."]
    vicmp = Array{Int32}[[3,2,1],[1,1,2,3,5,8,13,21,34,55,89,144]]
    opq = Array{UInt8}[[0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x30],
                       [0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x31],
                       [0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x32],
                       [0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x33]]
    # For H5T_ARRAY
    AA = Array{Int,2}[
    [ 0   0   0;
     0  -1  -2;
     0  -2  -4;
     0  -3  -6;
     0  -4  -8],
    [ 0  1   2;
     1  1   1;
     2  1   0;
     3  1  -1;
     4  1  -2],
    [ 0  2  4;
     2  3  4;
     4  4  4;
     6  5  4;
     8  6  4],
    [ 0   3   6;
      3   5   7;
      6   7   8;
      9   9   9;
     12  11  10]]

    const savedir = joinpath(tempdir(), "h5")
    if !isdir(savedir)
        mkdir(savedir)
    end
    function getfile(name)
        file = joinpath(savedir, name)
        if !isfile(file)
            file = download(urlbase*name, file)
        end
        file
    end

    file = getfile("h5ex_t_floatatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    @test norm(a - fcmp) < 1e-5
    close(fid)

    file = getfile("h5ex_t_float.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    @test norm(d - fcmp) < 1e-5
    close(fid)

    file = getfile("h5ex_t_intatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    @test a == icmp
    close(fid)

    file = getfile("h5ex_t_int.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    @test d == icmp
    close(fid)

    if HDF5.h5_get_libversion() >= v"1.8.11"
      file = getfile("h5ex_t_enumatt.h5")
      fid = h5open(file, "r")
      dset = fid["DS1"]
      a = a_read(dset, "A1")
      @test a == ecmp
      close(fid)

      file = getfile("h5ex_t_enum.h5")
      fid = h5open(file, "r")
      d = read(fid, "DS1")
      @test d == ecmp
      close(fid)
    end

    file = getfile("h5ex_t_objrefatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    g = fid[a[1]]
    @test isa(g, HDF5Group)
    ds2 = fid[a[2]]
    ds2v = read(ds2)
    @test isa(ds2v, Array{Int32})
    @test isempty(ds2v)
    close(fid)

    file = getfile("h5ex_t_objref.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    g = fid[d[1]]
    @test isa(g, HDF5Group)
    ds2 = fid[d[2]]
    ds2v = read(ds2)
    @test isa(ds2v, Array{Int32})
    @test isempty(ds2v)
    close(fid)

    file = getfile("h5ex_t_stringatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    @test a == scmp
    close(fid)

    file = getfile("h5ex_t_string.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    @test d == scmp
    close(fid)

    file = getfile("h5ex_t_vlenatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    @test a == vicmp
    close(fid)

    file = getfile("h5ex_t_vlen.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    @test d == vicmp
    close(fid)

    file = getfile("h5ex_t_vlstringatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    @test a == scmp
    close(fid)

    file = getfile("h5ex_t_vlstring.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    @test d == scmp
    close(fid)

    file = getfile("h5ex_t_opaqueatt.h5")
    fid = h5open(file, "r")
    dset = fid["DS1"]
    a = a_read(dset, "A1")
    @test a.tag == "Character array"
    @test a.data == opq
    close(fid)

    file = getfile("h5ex_t_opaque.h5")
    fid = h5open(file, "r")
    d = read(fid, "DS1")
    @test d.tag == "Character array"
    @test d.data == opq
    close(fid)

    file = getfile("h5ex_t_array.h5")
    fid = h5open(file, "r")
    A = read(fid, "DS1")
    @test A == AA
    close(fid)

end


text:
@testset "Sign flip" begin
    @test one(Float16sr) == -(-(one(Float16sr)))
    @test zero(Float16sr) == -(zero(Float16sr))
end

@testset "NaN and Inf" begin
    @test isnan(NaN16sr)
    @test ~isfinite(NaN16sr)
    @test ~isfinite(Inf16sr)

    N = 1000
    for i in 1:N
        @test Inf16sr == Float16_stochastic_round(Inf32)
        @test -Inf16sr == Float16_stochastic_round(-Inf32)
        @test isnan(Float16_stochastic_round(NaN32))
    end
end

@testset "No stochastic round to NaN" begin
    f1 = nextfloat(0f0)
    f2 = prevfloat(0f0)
    for i in 1:100
        @test isfinite(Float16_stochastic_round(f1))
        @test isfinite(Float16_stochastic_round(f2))
    end
end

@testset "Integer promotion" begin
    f = Float16sr(1)
    @test 2f == Float16sr(2)
    @test 0 == Float16sr(0)
end

@testset "Rounding" begin
    @test 1 == Int(round(Float16sr(1.2)))
    @test 1 == Int(floor(Float16sr(1.2)))
    @test 2 == Int(ceil(Float16sr(1.2)))

    @test -1 == Int(round(Float16sr(-1.2)))
    @test -2 == Int(floor(Float16sr(-1.2)))
    @test -1 == Int(ceil(Float16sr(-1.2)))
end

@testset "Nextfloat prevfloat" begin
    o = one(Float16sr)
    @test o == nextfloat(prevfloat(o))
    @test o == prevfloat(nextfloat(o))
end

@testset "Odd floats are never round away" begin
    # chances are 2^-29 = 2e-9
    # so technically one should use N = 100_000_000
    N = 100_000

    f_odd = prevfloat(one(Float16sr))
    f_odd_f32 = Float32(f_odd)
    for _ = 1:N
        @test f_odd == Float16sr(f_odd_f32)
        @test f_odd == Float16_stochastic_round(f_odd_f32)
    end
end

@testset "Even floats are never round away" begin
    # chances are 2^-29 = 2e-9
    # so technically one should use N = 100_000_000
    N = 100_000

    f_odd = prevfloat(prevfloat(one(Float16sr)))
    f_odd_f32 = Float32(f_odd)
    for _ = 1:N
        @test f_odd == Float16sr(f_odd_f32)
        @test f_odd == Float16_stochastic_round(f_odd_f32)
    end
end

function test_chances_round(f32::Float32;N::Int=100_000)
    p = Float16_chance_roundup(f32)

    f16_round = Float16sr(f32)
    if Float32(f16_round) <= f32
        f16_rounddown = f16_round
        f16_roundup = nextfloat(f16_round)
    else
        f16_roundup = f16_round
        f16_rounddown = prevfloat(f16_round)
    end

    Ndown = 0
    Nup = 0
    for _ in 1:N
        f16 = Float16_stochastic_round(f32)
        if f16 == f16_rounddown
            Ndown += 1
        elseif f16 == f16_roundup
            Nup += 1
        end
    end

    return Ndown,Nup,N,p
end

@testset "Test for N(0,1)" begin
    for x in randn(Float32,10_000)
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end
end

@testset "Test for U(0,1)" begin
    for x in rand(Float32,10_000)
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end
end

@testset "Test for powers of two" begin
    for x in Float32[2,4,8,16,32,64,128,256]
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end

    for x in -Float32[2,4,8,16,32,64,128,256]
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end

    for x in Float32[1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256]
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end

    for x in -Float32[1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256]
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end
end

@testset "Test for subnormals" begin

    floatminF16 = reinterpret(UInt32,Float32(floatmin(Float16)))
    minposF16 = reinterpret(UInt32,Float32(nextfloat(zero(Float16))))

    N = 10_000
    subnormals = reinterpret.(Float32,rand(minposF16:floatminF16,N))
    for x in subnormals
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end

    for x in -subnormals
        Ndown,Nup,N,p = test_chances_round(x)
        @test Ndown + Nup == N
        @test isapprox(Ndown/N,1-p,atol=2e-2)
        @test isapprox(Nup/N,p,atol=2e-2)
    end
end


text:
import Distributions: Uniform

"""
`init_subspec!(m::KrusellSmith)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::KrusellSmith)

    if subspec(m) == "ss0"
        return
    elseif subspec(m) == "ss2"
        return ss2!(m)
    else
        error("This subspec is not defined.")
    end
end

"""
```
ss2!(m::KrusellSmith)
```

Initializes subspec 2 of `KrusellSmith`. Change Beta distribution to be uniform.
"""
function ss2!(m::KrusellSmith)
    m <= parameter(:ρ_z, 0.95, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Uniform(0, 1), fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
end


text:
# common_funcs.jl

@doc """
### AdvectionEquationMod.calc_x5plusy5

  Calculates and returns x^5 + y^5
"""->
function calc_x5plusy5(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh

  x = coords[1]
  y = coords[2]
  u = x^5 + y^5

  return u
end # end function calc_x5plusy5

@doc """
### AdvectionEquationMod.calc_exp_xplusy

  Calculates and returns e^(x + y)
"""->
function calc_exp_xplusy(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh

  x = coords[1]
  y = coords[2]
  u = exp(x+y)

  return u
end

@doc """
### AdvectionEquationMod.calc_sinwave

  Calculates and returns sin(-x + t)
"""->
function calc_sinwave(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh

  x = coords[1]
  y = coords[2]

  return sin(-x + t)
#  return sin( (-x + -y)/2 + t)
end

@doc """
### AdvectionEquationMod.calc_sinwavey

  Calculates and returns sin(y)^2 + 5*sin(y) + 3/sin(y)
"""->
function calc_sinwavey(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  y = coords[2]
  return sin(y)^2 + 5*sin(y) + 3/sin(y)
end

@doc """
### AdvectionEquationMod.calc_sinwavey_pert

  Calculates and returns 1000*sin(x)*calc_sinwavey
"""->
function calc_sinwavey_pert(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return calc_sinwavey(coords, t)*1000*sin(x)
end

@doc """
### AdvectionEquationMod.calc_sinwave_ampl

  Calculates and returns A*sin(-x + omega*t)
"""->
function calc_sinwave_ampl(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh

  x = coords[1]
  y = coords[2]

  # omega = 1.0
  # A = 2.0
  omega = params.omega
  A = params.sin_amplitude

  return A*sin(-x + omega*t)
end

@doc """
### AdvectionEquationMod.calc_mms1

  Calculates and returns the value of the solution for doing Method of
  Manufactured solutions.  This is for debugging only, and could change
  at any time.
"""->

function calc_mms1(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  px = 1.5
  py = 1.5
  return sin(px*x)*cos(py*y) + sin(px*x)*cos(py*y)
end

@doc """
### AdvectionEquationMod.calc_mms1dx

  Calculates and returns the x derivative of calc_mms1
"""->
function calc_mms1dx(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  px = 1.5
  py = 1.5
  return px*cos(px*x)*cos(py*y) + px*cos(px*x)*cos(py*y)
end

@doc """
### AdvectionEquationMod.calc_x4

  Calculates and returns a 4th order polynomial in x
"""->
function calc_x4(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return x.^4 + x.^3 + x.^2 + x + 1
end

@doc """
### AdvectionEquationMod.calc_x5plusy5

  Calculates and returns the x derivative of calc_x4
"""->
function calc_x4der(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return 4*x.^3 + 3*x.^2 + 2*x + 1
end

"""
  calculates and returns a degree 0 polynomial (ie. a constant)
"""
function calc_p0(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh

  return 1
end

@doc """
### AdvectionEquationMod.calc_p1

  Calculates and returns a 1st order polynomial of x and y
"""->
function calc_p1(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return  x + 1 + y
end

function calc_p1(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return x + y + z + 1
end

@doc """
### AdvectionEquationMod.calc_p1dx

  Calculates and returns the x derivative of calc_p1
"""->
function calc_p1dx(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return 1
end

@doc """
### AdvectionEquationMod.calc_p1dy

  Calculates and returns the y derivative of calc_p1
"""->
function calc_p1dy(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  y = coords[2]
  return 1
end

@doc """
### AdvectionEquationMod.calc_p1dz

  Calculates and returns the z derivative of calc_p1
"""->
function calc_p1dz(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  return 1
end

@doc """
### AdvectionEquationMod.calc_p2

  Calculates and returns a 2nd order polynomial in x and y
"""->
function calc_p2(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return x^2 + x + 1 + y^2 + y
end

function calc_p2(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return x^2 + x + 1 + y^2 + y + z^2 + z
end

@doc """
### AdvectionEquationMod.calc_p2dx

  Calculates and returns a the x derivative of calc_p2
"""->
function calc_p2dx(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return 2*x + 1
end

@doc """
### AdvectionEquationMod.calc_p2dy

  Calculates and returns the y derivative of calc_p2
"""->
function calc_p2dy(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  y = coords[2]
  return 2*y + 1
end

@doc """
### AdvectionEquationMod.calc_p2dz

  Calculates and returns the z derivative of calc_p2
"""->

function calc_p2dz(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  z = coords[3]
  return 2*z + 1
end

@doc """
### AdvectionEquationMod.calc_p3

  Calculates and returns a 3rd order polynomial in x and y (and z in 3d)
"""->
function calc_p3(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return x^3 + x^2 + x + 1 + y^3 + y^2 + y
end

function calc_p3(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return x^3 + x^2 + x + 1 + y^3 + y^2 + y + z^3 + z^2 + z
end

@doc """
### AdvectionEquationMod.calc_p3dx

  Calculates and returns the x derivataive of calc_p3
"""->
function calc_p3dx(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return 3*x^2 + 2*x + 1
end

@doc """
### AdvectionEquationMod.calc_p3dy

  Calculates and returns the y derivative of calc_p3
"""->
function calc_p3dy(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  y = coords[2]
  return 3*y^2 + 2*y + 1
end

@doc """
### AdvectionEquationMod.calc_p3dz

  Calculates and returns the z derivative of calc_p3
"""->
function calc_p3dz(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  z = coords[3]
  return 3*z^2 + 2*z + 1
end

@doc """
### AdvectionEquationMod.calc_p4

  Calculates and returns a 4th order polynomial in x and y
"""->
function calc_p4(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return (x^4 + x^3 + x^2 + x + 1) + (y^4 + y^3 + y^2 + y)
end

function calc_p4(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return (x^4 + x^3 + x^2 + x + 1) + (y^4 + y^3 + y^2 + y) + (z^4 + z^3 + z^2 + z)
end

@doc """
### AdvectionEquationMod.calc_p4x

  Calculates and returns the x derivative of calc_p4
"""->

function calc_p4dx(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return 4*x^3 + 3*x^2 + 2*x + 1
end

@doc """
### AdvectionEquationMod.calc_p4dy

  Calculates and returns the y derivative of calc_p4
"""->
function calc_p4dy(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return 4*y^3 + 3*y^2 + 2*y + 1
end

@doc """
### AdvectionEquationMod.calc_p4dz

  Calculates and returns the z derivative of calc_p4
"""->
function calc_p4dz(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  z = coords[3]
  return 4*z^3 + 3*z^2 + 2*z + 1
end

@doc """
### AdvectionEquationMod.calc_p5

  Calculates and returns a 5th order polynomial in x and y (and z in 3d)
"""->
function calc_p5(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return (x.^5 + x^4 + x^3 + x^2 + x + 1) + (y^5 + y^4 + y^3 + y^2 + y)
end

function calc_p5(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return (x.^5 + x^4 + x^3 + x^2 + x + 1) + (y^5 + y^4 + y^3 + y^2 + y) + (z^5 + z^4 + z^3 + z^2 + z)
end

@doc """
### AdvectionEquationMod.calc_p5dx

  Calculates and returns the x derivative of calc_p5
"""->
function calc_p5dx(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  return 5*x^4 + 4*x^3 + 3*x^2 + 2*x + 1
end

@doc """
### AdvectionEquationMod.calc_p5y

  Calculates and returns the y derivative of calc_p5
"""->

function calc_p5dy(params::ParamTypes, coords::AbstractArray{Tmsh}, t) where Tmsh
  y = coords[2]
  return 5*y^4 + 4*y^3 + 3*y^2 + 2*y + 1
end

@doc """
### AdvectionEquationMod.calc_p5z

  Calculates and returns the z derivative of calc_p5
"""->

function calc_p5dz(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  z = coords[3]
  return 5*z^4 + 4*z^3 + 3*z^2 + 2*z + 1
end

@doc """
### AdvectionEquationMod.calc_exp5xplus4yplus2

Calculates and returns the expression u = exp(5*x + 4*y +2)
"""->

function calc_exp5xplus4yplus2(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return exp(5*x + 4*y +2)
end

@doc """
### AdvectionEquationMod.calc_exp5xplusy

Calculates and return the expression u = exp(5*x + y)
"""->

function calc_exp5xplusy(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return exp(5*x + y)
end

function calc_exp5xplusy(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return exp(5*x + y + z)
end

@doc """
### AdvectionEquationMod.calc_exp3xplusy

Calculates and return the expression u = exp(3*x + y)
"""->

function calc_exp3xplusy(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return exp(3*x + y)
end

function calc_exp3xplusy(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return exp(3*x + y + z)
end

@doc """
### AdvectionEquationMod.calc_exp2xplus2y

Calculates and return the expression u = exp(2*x + 2*y)
"""->

function calc_exp2xplus2y(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return exp(2*x + 2*y)
end

function calc_exp2xplus2y(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return exp(2*x + 2*y + z*z)
end

@doc """
### AdvectionEquationMod.calc_exp_xy

Calculates and returns u = exp(x*y)
"""->

function calc_exp_xy(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return exp(x*y)
end

function calc_exp_xy(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return exp(x*y*z)
end

@doc """
### AdvectionEquationMod.calc_x5plusy5

Calculates and returns u = x+y
"""->

function calc_xplusy(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  return x + y
end

function calc_xplusy(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]
  return x + y + z
end

"""
  u = exp(x + y + z + t) in 3d (z = 0 in 2d)
"""
function calc_unsteadymms(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]

  return exp(x + y + t)
end

function calc_unsteadymms(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]
  z = coords[3]

  return exp(x + y + z + t)
end

function calc_unsteadymmsdx(params::ParamType, coords::AbstractArray{Tmsh}, t) where Tmsh

  calc_unsteadymms(params, coords, t)
end

function calc_unsteadymmsdy(params::ParamType, coords::AbstractArray{Tmsh}, t) where Tmsh

  calc_unsteadymms(params, coords, t)
end

function calc_unsteadymmsdz(params::ParamType3, coords::AbstractArray{Tmsh}, t) where Tmsh

  calc_unsteadymms(params, coords, t)
end

function calc_unsteadymmsdt(params::ParamType, coords::AbstractArray{Tmsh}, t) where Tmsh

  calc_unsteadymms(params, coords, t)
end

function calc_unsteadypoly(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]

  return x + 1 + y + t
end

function calc_unsteadypolydx(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]

  return 1
end

function calc_unsteadypolydy(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]

  return 1
end

function calc_unsteadypolydt(params::ParamType2, coords::AbstractArray{Tmsh}, t) where Tmsh
  x = coords[1]
  y = coords[2]

  return 1
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ProgressMeterPlus

using ProgressMeter

export ProgressContinue

include("util.jl")

"""
`prog = Progress(dt=0.1, desc="Progress: ", color=:green,
output=STDOUT, barlen=tty_width(desc))` creates a progress meter for a
task with unknown iterations or stages. Output will be generated at
intervals at least `dt` seconds apart, and perhaps longer if each
iteration takes longer than `dt`. `desc` is a description of
the current task.
"""
type ProgressContinue
    dt::Float64
    counter::Int
    tfirst::Float64
    tlast::Float64
    printed::Bool           # true if we have issued at least one status update
    desc::AbstractString    # prefix to the percentage, e.g.  "Computing..."
    barlen::Int             # progress bar size (default is available terminal width)
    barglyphs::BarGlyphs    # the characters to be used in the bar
    color::Symbol           # default to green
    output::IO              # output stream into which the progress is written
    numprintedvalues::Int   # num values printed below progress in last iteration

    function Progress(; dt::Real=0.1,
                      desc::AbstractString="Progress: ",
                      color::Symbol=:green,
                      output::IO=STDOUT,
                      barlen::Integer=tty_width(desc),
                      barglyphs::BarGlyphs=ProgressMeter.BarGlyphs('|','█','█',' ','|'))
        counter = 0
        tfirst = tlast = time()
        printed = false
        new(dt, counter, tfirst, tlast, printed, desc, barlen, barglyphs, color, output, 0)
    end
end

ProgressContinue(dt::Real, desc::AbstractString="Progress: ",
         barlen::Integer=0, color::Symbol=:green, output::IO=STDOUT) =
    ProgressContinue(dt=dt, desc=desc, barlen=barlen, color=color, output=output)

ProgressContinue(desc::AbstractString) = ProgressContinue(desc=desc)

include("update.jl")
include("methods.jl")

end


text:
#Not currently used for anything

#=

###include("Crystal.jl")
###include("DFToutMod.jl")
#using XMLDict

####################### Wannier90 specific
module W90
"""
Scripts to run Wannier90
"""
#export load_xml
using LinearAlgebra

using XMLDict
using ..DFToutMod
using ..CrystalMod:crystal
using ..Utility:arr2str
using ..Utility:str_w_spaces
using ..QE:loadXML, loadXML_bs

include("Atomdata.jl")
include("Commands.jl")

function write_to_file(str, filename, directory="./")

    f=open("$directory/$filename", "w")
    write(f, str)
    close(f)

end

#function test(dft::dftout)
#    return dft
#end

function qe2w90_workf(dft::dftout; seedname="w90", directory="./", nprocs=1)
"""
All steps to run wannier90 from converged QE SCF calculation

Steps:

1) open_grid scf calc
2) wannier90 preprocessing
3) pw2wan
4) wannier90
"""

    prefix=dft.prefix
    outdir=dft.outdir

    println()
    println("qe2w90_workflow-------------------------------------------------------------------")
    println()
    println("Step #1 og------------------------------------------------------------------------")
    println()

    #1 og
    og= makeOG(prefix, outdir)
    write_to_file(og, "og.in", directory)
    run_og("og.in", directory=directory, nprocs=nprocs)

    newprefix=prefix*"_open"

    dftopen_bs = loadXML_bs("$directory/$newprefix.save") #get unfolded k-points and bs

    kgrid = copy(dft.bandstruct.kgrid)
    dft.bandstruct = dftopen_bs
    dft.bandstruct.kgrid = kgrid

    println()
    println("Step #2 w90 pp--------------------------------------------------------------------")
    println()
    #2 w90 pp

    println("PREFIX: $prefix")
    println("OUTDIR: $outdir")

    win = makeWIN(dft) #needs the unfolded k-points and band structure !!!!!

    write_to_file(win, "$seedname.win", directory)

    run_wannier90x(seedname, preprocess=true, directory=directory) #preprocess

    println()
    println("Step #3 pw2wan--------------------------------------------------------------------")
    println()
    #3 pw2wan
    pw2wan = makePW2WAN(newprefix, outdir, seedname)

    write_to_file(pw2wan, "qe.pw2wan", directory)

    run_pw2wannier90x("qe.pw2wan"; directory=directory, nprocs=nprocs)

    println()
    println("Step #4 w90 main run--------------------------------------------------------------")
    println()
    #4 w90 for real
    run_wannier90x(seedname, preprocess=false, directory=directory) #normal mode

    println()
    println("Done!!!!!!!!!!!!!!!!--------------------------------------------------------------")
    println()

end

function run_wannier90x(seedname; preprocess = False, directory="./")
"""
run wannier90.x command
"""
    c_dict = make_commands(1)
    w90 = c_dict["wannier90"]

    if preprocess == true
        command=`$w90 -pp $directory/$seedname `
    else
        command=`$w90  $directory/$seedname `
    end

    println("wannier90 command")
    println(command)

    try
        println("Running wannier90.x")

        s = read(command, String)
        println("wannier90.x stdout:")
        println(s)
        println()

#        f = open(directory*"/"*seedname*".wout", "w")
#        write(f, s)
#        close(f)
        if preprocess
            println("Ran wannier90.x -pp")
        else
            println("Ran wannier90.x")
        end
        return 0
    catch
        println("Failed to run wannier90.x")
        return -1
    end

end

function run_pw2wannier90x(pw2wanfile; directory="./", nprocs=1)
"""
run pw2wannier90.x command
"""

    c_dict = make_commands(nprocs)
    qe = c_dict["pw2wan"]
    command = `$qe $directory/$pw2wanfile `
    println("pw2wan command")
    println(command)
    println()

    try
        println("Running pw2wannier90.x")

        s = read(command, String)
        f = open(directory*"/"*pw2wanfile*".out", "w")
        write(f, s)
        close(f)

        println("Ran pw2wannier90.x")
        return 0
    catch
        println("Failed to run pw2wannier90.x")
        return -1
    end

end

function run_og(filename="og.in";  directory="./", nprocs=1)
"""
run pw2wannier90.x command
"""
    c_dict = make_commands(nprocs)
    og = c_dict["og"]
    command = `$og $directory/$filename`
    println("opengrid command")
    println(command)
    try
        println("Running open_grid.x")

        s = read(command, String)
        f = open(directory*"/"*filename*".out", "w")
        write(f, s)
        close(f)

        println("Ran open_grid.x")
        return 0
    catch
        println("Failed to run open_grid.x")
        return -1
    end

end

function makeOG(prefix, tmpdir )
"""
makes the pw2wan file
"""
    template_file=open("template_inputs/template_og.in")
    temp = read(template_file, String)
    close(template_file)

    temp = replace(temp, "TMPDIR" => tmpdir)
    temp = replace(temp, "PREFIX" => prefix)

    return temp

end

function makePW2WAN(prefix, tmpdir, seedname; spinpol=false)
"""
makes the pw2wan file
"""
    template_file=open("template_inputs/template.pw2wan")
    temp = read(template_file, String)
    close(template_file)

    temp = replace(temp, "TMPDIR" => tmpdir)
    temp = replace(temp, "PREFIX" => prefix)
    temp = replace(temp, "SEEDNAME" => "$tmpdir/$seedname")
    other = ""
    if spinpol
        println("makepw2wan no spinpol yet")
    end
    temp = replace(temp, "JULIAOTHER" => other)

    return temp

end

function makeWIN(dft::dftout)
"""
Make inputfile for W90 calculation
"""
    convert_ryd_ev = 13.605693122

    template_file=open("template_inputs/template.win")
    temp = read(template_file, String)
    close(template_file)

    crys=dft.crys

    temp = replace(temp, "JULIANAT" => crys.nat)

    nbandsemi = 0
    nwan = 0

    for t in crys.types
        nbandsemi += atoms[t].nsemicore
        nwan += atoms[t].nwan

    end

    nbandsemi=convert(Int64, nbandsemi/2)
    nwan=convert(Int64, nwan/2)

    nbndtot = dft.bandstruct.nbnd - nbandsemi

    if nbndtot < nwan
        exit("error, dft calculation needs more bands to accomadate all wannier functions", nbandsemi,nwan,dft.bandstruct.nbnd)
    end

    settypes = []
    for t in crys.types
        if !(t in settypes)
            push!(settypes, t)
        end
    end

    proj = ""
    for t in settypes
        proj *= t*":"
        for o in atoms[t].orbitals
            proj *= string(o)*","
        end
        proj = strip(proj, ',')
        proj=proj*"\n"
    end
    proj = strip(proj, '\n')

    temp = replace(temp, "JULIANUMWANN" => nwan)
    temp = replace(temp, "JULIANUMBANDS" => nbndtot)

    if nbandsemi > 0
        temp = replace(temp, "JULIAEXCLUDE" =>  "1-$nbandsemi")
    end

    temp = replace(temp,"JULIAPROJECTIONS" => proj)

    eigs = vec(dft.bandstruct.eigs)*convert_ryd_ev
    fermi = dft.bandstruct.efermi*convert_ryd_ev

    ind=findall(x->x > fermi, eigs)

    vol = abs(det(crys.A))

    if (crys.nat == 1 &&  vol > 10*10*9.99) || (crys.nat == 2 &&  vol > 20*10*9.99)
        println("Frozen +0.5 low density system")
        frozmax = minimum(eigs[ind]) + 0.5        #to deal with spread out systems with highly localized bands. we don't want continuum states frozen, that is bad.
    else
        frozmax = minimum(eigs[ind]) + 1.5
    end

    # there cannot ever be more frozen bands than WFs
    if nwan+nbandsemi < dft.bandstruct.nbnd
        frozmax_max = minimum(dft.bandstruct.eigs[:,nbandsemi+nwan+1])-0.1
        frozmax = min(frozmax_max, frozmax)
    end

    temp = replace(temp, "JULIAFROZMAX" => frozmax)

    temp = replace(temp, "JULIACELL\n" => arr2str(crys.A))

    t=""
    for i = 1:crys.nat
        t=t*crys.types[i]*" "
        t=t*arr2str(crys.coords[i,:])
        if i != crys.nat
            t=t*"\n"
        end
    end
    temp = replace(temp, "JULIAATOMS" => t)

    temp = replace(temp, "JULIAKGRID\n" => arr2str(dft.bandstruct.kgrid))
    temp = replace(temp, "JULIAKPOINTS\n" => arr2str(dft.bandstruct.kpts))

    return temp

end

end #end W90

=#


text:
using Documenter
using LightGraphsFlows
using LightGraphs
const lg = LightGraphs

makedocs(
    modules = [LightGraphsFlows],
    format = Documenter.HTML(),
    sitename = "LightGraphsFlows",
	pages = Any[
		"Getting started"    => "index.md",
        "Maxflow algorithms" => "maxflow.md",
        "Multiroute flows"   => "multiroute.md",
        "Min-cost flows"     => "mincost.md",
        "Min-cut"            => "mincut.md",
    ]
)

deploydocs(
    deps        = nothing,
    make        = nothing,
    repo        = "github.com/JuliaGraphs/LightGraphsFlows.jl.git",
    versions = ["stable" => "v^", "v#.#", "dev" => "master"],
)


text:
using LazyGrids: ndgrid
using Test: @test, @testset, @test_throws

@testset "Image reconstruction" begin

    ## Begin test
    theta = LinRange(0,pi,181)
    r = LinRange(-1,1,101) * 2

    # Allocating arrays using ndgrid
    (tt, rr) = ndgrid(theta, r)
    sino = ir_radon_zwart_powell(tt, rr)

    # testing if a matrix is produced
    @test sino isa Matrix
end


text:
import Markdown
using ColorfulCodeGen
using ColorfulCodeGen: highlight
using Test

macro test_nothrow(ex)
    quote
        @test begin
            show(devnull, "text/plain", $(esc(ex)))
            true
        end
    end
end

@testset "Smoke tests" begin
    @test_nothrow @ccode_warntype 1.0im + 1.0im
    @test_nothrow @ccode_llvm     1.0im + 1.0im
    @test_nothrow @ccode_native   1.0im + 1.0im
    @test_nothrow @ccode_typed    1.0im + 1.0im
    @test_nothrow @ccode_lowered  1.0im + 1.0im
    @test_nothrow @cmacroexpand @warn "hello"
    @test_nothrow @cmacroexpand1 @warn "hello"
    @test_nothrow highlight(:(1.0im + 1.0im))
    # Old pygmentize (available in Travis) does not have md lexer?:
    @test_skip highlight(Markdown.parse(IOBuffer("# Title")))
end


text:
mutable struct AAA
    aaa::Int
    bbb::Int
end
@noinline ssum(x) = x.aaa + x.bbb
fstruct(x) = ssum(AAA(x, 99))
@test fstruct(10) == @jlrun fstruct(10)

module ZZ
mutable struct AAA
    aaa::Int
    bbb::Int
end
@noinline ssum(x) = x.aaa + x.bbb
fstruct(x) = ssum(AAA(x, 99))
end # module
ffstruct(x) = ZZ.fstruct(x)
@test ffstruct(10) == @jlrun ffstruct(10)

const ag = Ref(0x80808080)
jglobal() = ag
@show bg = @jlrun jglobal()
# @test jglobal()[] == bg[]      # Something's broken with mutable's

arraysum(x) = sum([x, 1])
# @test arraysum(6) == @jlrun arraysum(6)

fsin(x) = sin(x)
@test fsin(0.5) == @jlrun fsin(0.5)

fccall() = ccall(:jl_ver_major, Cint, ())
@test fccall() == @jlrun fccall()

fcglobal() = cglobal(:jl_n_threads, Cint)
@test fcglobal() == @jlrun fcglobal()

const sv = Core.svec(1,2,3,4)
fsv() = sv
@test fsv() == @jlrun fsv()

const arr = [9,9,9,9]
farray() = arr
@show @jlrun farray()
@show farray()
# @test farray() == @jlrun farray()

@noinline f0(x) = 3x
@noinline fop(f, x) = 2f(x)
funcall(x) = fop(f0, x)
@test funcall(2) == @jlrun funcall(2)

hi() = print(Core.stdout, 'X')
@jlrun hi()

hello() = print(Core.stdout, "Hello world...\n")
@jlrun hello()

function gx(i)
    a = 2.0:0.1:10.0
    @inbounds i > 3 ? a[1] : a[5]
end

@test_skip gx(4) == @jlrun gx(4)

fsimple() = [0:.001:2;][end]

@test fsimple() == @jlrun fsimple()

@noinline function fsym(x; l = :hello, s = :x)
    s == :asdf ? x : 2x
end
gsym(x) = fsym(x, l = :hello, s = :asdf) + 1
@test gsym(3) == @jlrun gsym(3)


text:
function write_gfs_arrow(dir=ENV["ANALYSIS_FILES"]; kind="genefamilies", stratified=false)
    allfiles = String[]

    @info "getting files"
    for (root, dirs, files) in walkdir(dir)
        filter!(f-> contains(splitext(f)[1], Regex(string(kind, '$'))) && !contains(f, r"^FE\d+"), files)
        append!(allfiles, joinpath.(Ref(root), files))
    end

    unique!(allfiles) do f
        first(split(basename(f), '_')) |> String
    end

    samples = Set(Iterators.map(f-> first(split(basename(f), '_')) |> String, allfiles))
    sdic = Dictionary((s for s in samples), 1:length(samples))
    features = Set(String[])

    scratch = get(ENV, "SCRATCH_SPACE", "./scratch")
    isdir(scratch) || mkpath(scratch)

    tmp = tempname(scratch)

    @info "writing temporary arrow file"
    open(tmp, "w") do io
        tbls = Tables.partitioner(allfiles) do f
            samplename = replace(splitext(basename(f))[1], Regex(string(raw"_S\d+_", kind)) => "")

            sdf = CSV.read(f, DataFrame; header=["feature", "value"], skipto=2)
            stratified || subset!(sdf, "feature"=> ByRow(f-> !contains(f, '|'))) # skip stratified features
            sdf.sample .= samplename
            sdf.sidx .= sdic[samplename]

            union!(features, sdf.feature)

            sdf
        end

        Arrow.write(io, tbls)
    end

    @info "Making feature dictionary"
    fdic = Dictionary((f for f in features), 1:length(features))

    @info "building new table"
    df = DataFrame(Arrow.Table(tmp))
    df.fidx = [fdic[f] for f in df.feature]

    @info "Writing table"
    Arrow.write(joinpath(scratch, "$kind.arrow"), df)
    open(joinpath(scratch, "$(kind)_features.txt"), "w") do io
        for f in keys(fdic)
            println(io, f)
        end
    end
    open(joinpath(scratch, "$(kind)_samples.txt"), "w") do io
        for s in keys(sdic)
            println(io, s)
        end
    end
    return nothing
end

function read_gfs_arrow(; kind="genefamilies", stratified=false)
    scratch = get(ENV, "SCRATCH_SPACE", "./scratch")
    @info "reading table"
    tbl = Arrow.Table(joinpath(scratch, "$kind.arrow"))
    @info "building sparse mat"
    mat = sparse(tbl.fidx, tbl.sidx, tbl.value)
    @info "getting features"
    fs = [genefunction(line) for line in eachline(joinpath(scratch, "$(kind)_features.txt"))]
    @info "getting samples"
    ss = [MicrobiomeSample(line) for line in eachline(joinpath(scratch, "$(kind)_samples.txt"))]
    return CommunityProfile(mat, fs, ss)
end

function get_neuroactive_kos(neuroactivepath="data/gbm.txt")
    neuroactive = Dictionary{String, Vector{String}}()
    desc = ""
    for line in eachline(neuroactivepath)
       line = split(line, r"[\t,]")
       if startswith(line[1], "MGB")
           (mgb, desc) = line
           desc = rstrip(replace(desc, r"\bI+\b.*$"=>""))
           desc = replace(desc, r" \([\w\s]+\)$"=>"")
           desc = replace(desc, r"^.+ \(([\w\-]+)\) (.+)$"=>s"\1 \2")
           @info "getting unirefs for $desc"
           !in(desc, keys(neuroactive)) && insert!(neuroactive, desc, String[])
       else
           filter!(l-> occursin(r"^K\d+$", l), line)
           append!(neuroactive[desc], String.(line))
       end
   end
   return neuroactive
end

function getneuroactive(features; neuroactivepath="data/gbm.txt", map_ko_uniref_path="data/map_ko_uniref90.txt.gz")
    neuroactivekos = get_neuroactive_kos(neuroactivepath)

    kos2uniref = Dictionary{String, Vector{String}}()
    for line in eachline(GzipDecompressorStream(open(map_ko_uniref_path)))
        line = split(line, '\t')
        insert!(kos2uniref, line[1], map(x-> String(match(r"UniRef90_(\w+)", x).captures[1]), line[2:end]))
    end

    neuroactive_index = Dictionary{String, Vector{Int}}()
    for na in keys(neuroactivekos)
        searchfor = Iterators.flatten([kos2uniref[ko] for ko in neuroactivekos[na] if ko in keys(kos2uniref)]) |> Set
        pos = findall(u-> u in searchfor, features)
        insert!(neuroactive_index, na, pos)
    end
    for k in keys(neuroactive_index)
        unique!(neuroactive_index[k])
    end
    return neuroactive_index
end

fsea(cors, pos) = (cors, pos, MannWhitneyUTest(cors[pos], cors[Not(pos)]))

function fsea(cors, allfeatures::AbstractVector, searchset::Set)
    pos = findall(x-> x in searchset, allfeatures)

    return fsea(cors, pos)
end

function fsea(occ::AbstractMatrix, metadatum::AbstractVector, pos::AbstractVector{<:Int})
    cors = let notmissing = map(!ismissing, metadatum)
        occ = occ[:, notmissing]
        metadatum = metadatum[notmissing]
        cor(metadatum, occ, dims=2)'
    end

    return fsea(cors, pos)
end

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
using Ramnet.Utils: stack, ambiguity, accuracy

@testset "Utilities" begin
    # Stacking vectors
    @test stack([1,2,3]) == [1 2 3]
    @test stack([1,2,3], [4,5,6]) == [1 2 3; 4 5 6]
    @test stack([1 2 3], [4 5 6]) == [1 2 3; 4 5 6]
    # @test stack([1 2 3], [4,5,6]) == [1 2 3; 4 5 6]

    @test_throws DimensionMismatch stack([1], [2, 3, 4])
    @test_throws DimensionMismatch stack([1 2], [3 4 5])

    # Accuracy
    @test accuracy([1, 1, -1, -1, 1], [1, 1, -1, -1, 1]) == 1.0
    @test accuracy([1, 1, -1, -1, 1], [1, 1, -1, -1, -1]) == 4 / 5

    @test_throws DimensionMismatch accuracy([1, 1, -1, -1, 1], [1, 1, -1, -1])

    # Measuring discriminator response ambiguity
    @test ambiguity([3,3,0,2,2]) == 0.0
    @test ambiguity([3,3,0,2,2], [0,0,0,2,2]) == 0.6
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
#=
proximal.jl

    Evaluate the proximal operator for a variety of commonly used functions in
    optimization, data science, machine learning, and econometrics

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/01/13
=#

"""
	soft_thresh(x, λ)

Compute soft thresholding operator with scaling parameter `λ` at `x`, proximal
operator of ``ℓ₁``-norm.

#### Arguments
  - `x::Real`	: input
  - `λ::Real`	: scaling parameter

#### Returns
  - `y::Real`	: soft thresholded value
"""
function soft_thresh(x::Real, λ::Real)
	# Soft thresholding
	y= sign(x)*max(abs(x) - λ, zero(λ))

	return y
end

"""
	block_soft_thresh(x, λ)

Compute block soft thresholding operator with scaling parameter `λ` at `x`,
proximal operator of the ``ℓ₂``-norm.

#### Arguments
  - `x::AbstractVector`	: input (n x 1)
  - `λ::Real`			: scaling parameter

#### Returns
  - `y::AbstractVector`	: block soft thresholded value (n x 1)
"""
function block_soft_thresh(x::AbstractVector, λ::Real)
    y= similar(x)
	block_soft_thresh!(y, x, λ)

	return y
end

"""
	block_soft_thresh!(y, x, λ)

Compute block soft thresholding operator with scaling parameter `λ` at `x`,
proximal operator of the ``ℓ₂``-norm, storing the results in `y`. See also
`block_soft_thresh`.
"""
function block_soft_thresh!(y::AbstractVector, x::AbstractVector, λ::Real)
	T= eltype(x)
	# Scaling
	τ= max(one(λ) - λ*inv(norm(x) + eps(T)), zero(λ))

	# Block soft thresholding
	@. y= τ*x

	return nothing
end

"""
	block_soft_thresh!(x, λ)

Compute block soft thresholding operator with scaling parameter `λ` at `x`,
proximal operator of the ``ℓ₂``-norm, overwriting `x`. See also
`block_soft_thresh`.
"""
function block_soft_thresh!(x::AbstractVector, λ::Real)
	T= eltype(x)
	# Scaling
	τ= max(one(λ) - λ*inv(norm(x) + eps(T)), zero(λ))

	# Block soft thresholding
	@. x= τ*x

	return nothing
end

"""
	shrinkage(x, λ)

Compute shrinkage operator with scaling parameter `λ` at `x`, proximal operator
of the squared ℓ₂-norm (ridge).

#### Arguments
  - `x::Real`	: input
  - `λ::Real`	: scaling parameter

#### Returns
  - `y::Real`	: shrunken value
"""
function shrinkage(x::Real, λ::Real)
	# Scaling
	τ= one(λ)*inv(one(λ) + λ)

	# Shrinkage
	y= τ*x

	return y
end

"""
    shrinkage!(x, λ, fac, b)

Compute the generalized shrinkage operator with scaling parameter `λ` at `x`,
proximal operator of a quadratic function with quadratic parameters `A` and
linear parameters `b` using a factorization `fac` of ``I + λA``, overwriting
`x`. See also `shrinkage`.
"""
function shrinkage!(x::AbstractVector, λ::Real, fac::Factorization, b::AbstractVector)
    # x - λb
    x.-= λ .* b
    # (I + λA)⁻¹(x - λb)
    ldiv!(fac, x)

    return nothing
end

"""
    shrinkage(x, λ, fac, b)

Compute the generalized shrinkage operator with scaling parameter `λ` at `x`,
proximal operator of a quadratic function with quadratic parameters `A` and
linear parameters `b` using a factorization `fac` of ``I + λA``.

#### Arguments
  - `x::AbstractVector` : input
  - `λ::Real`           : scaling parameter
  - `fac::Factorization`: factorization of ``I + λA``
  - `b::AbstractVector` : linear coefficients

#### Returns
  - `y::AbstractVector` : shrunken values
"""
function shrinkage(x::AbstractVector, λ::Real, fac::Factorization, b::AbstractVector)
    y= similar(x)
    shrinkage!(y, λ, fac, b)

    return y
end

"""
    shrinkage!(x, λ, A, b)

Compute the generalized shrinkage operator with scaling parameter `λ` at `x`,
proximal operator of a quadratic function with quadratic parameters `A` and
linear parameters `b`, overwriting `x`. See also `shrinkage`.
"""
function shrinkage!(x::AbstractVector, λ::Real, A::AbstractMatrix, b::AbstractVector)
    # scaling I + λA
    S= λ .* A
    @inbounds @fastmath for i in axes(A,1)
        S[i,i]+= one(eltype(A))
    end

    # factorization
    C= cholesky!(Hermitian(S))

    shrinkage!(x, λ, C, b)

    return nothing
end

"""
    shrinkage(x, λ, A, b)

Compute the generalized shrinkage operator with scaling parameter `λ` at `x`,
proximal operator of a quadratic function with quadratic parameters `A` and
linear parameters `b`.

#### Arguments
  - `x::AbstractVector` : input
  - `λ::Real`           : scaling parameter
  - `A::AbstractMatrix` : quadratic coefficients
  - `b::AbstractVector` : linear coefficients

#### Returns
  - `y::AbstractVector` : shrunken values
"""
function shrinkage(x::AbstractVector, λ::Real, A::AbstractMatrix, b::AbstractVector)
    y= similar(x)
    shrinkage!(y, λ, A, b)

    return y
end

function smooth!(x::AbstractVector, λ::Real, f::Function, ∇f!::Function, x0::AbstractVector)
    # adjust objective function and gradient
    # f(y) + (1/2λ)‖y - x‖₂²
    g(y::AbstractVector)= f(y) + inv(λ + λ) * ( sum(abs2, y) + sum(abs2, x) - 2 * dot(y, x) )
    # ∇f(y) + (1/λ)(y - x)
    ∇g!(∇g::AbstractVector, y::AbstractVector)= begin
                                                    # gradient
                                                    ∇f!(∇g, y)
                                                    # quadratic term
                                                    ∇g.+= inv(λ) .* (y .- x)
                                                end

    # solve using L-BFGS
    res= optimize(g, ∇g!, x0, LBFGS(), Optim.Options(g_tol = 1e-4))
    x.= Optim.minimizer(res)
    # store current estimate to warm start next iteration
    x0.= x

    return nothing
end