text:
"""
Generate figures 4, 5, .... which summarize a test run of the Gibbs.
=====================================================
This file generates figures based on the run from `scriptParallel`.
"""

"""
Load modules and set load and save links
------------------------------------------
"""
simdir   =  "/Users/ethananderes/Dropbox/BayesLense/simulations/scriptParallel_2657077506/" # contains the simulation
srcpath  =  "/Users/ethananderes/Dropbox/BayesLense/src/"
savepath = "/Users/ethananderes/Dropbox/BayesLense/paper/"

include(srcpath*"Interp.jl")
include(srcpath*"cmb.jl")
include(srcpath*"fft.jl")

using PyPlot, Interp

"""
Set the parameters  of the simulation
-----------------------------------------
"""
const percentNyqForC = 0.5 # used for T l_max
const numofparsForP  = 1500  # used for P l_max
const hrfactor = 2.0
const pixel_size_arcmin = 2.0
const n = 2.0^9
const beamFWHM = 0.0
const nugget_at_each_pixel = (4.0)^2
begin  #< ---- dependent run parameters
	local deltx =  pixel_size_arcmin * π / (180 * 60) #rads
	local period = deltx * n # side length in rads
	local deltk =  2 * π / period
	local nyq = (2 * π) / (2 * deltx)
	const maskupP  = √(deltk^2 * numofparsForP / π)  #l_max for for phi
	const maskupC  = min(9000.0, percentNyqForC * (2 * π) / (2 * pixel_size_arcmin * π / (180*60))) #l_max for for phi
end
const scale_grad =  2.0e-3
const scale_hmc  =  2.0e-3

d = 2
parlr = setpar(
	pixel_size_arcmin,
	n,
	beamFWHM,
	nugget_at_each_pixel,
	maskupC,
	maskupP,
	srcpath
	);
dirac_0 = 1 / parlr.grd.deltk^d
elp = parlr.ell[10:int(maskupP)]
elt = parlr.ell[10:1000]
cpl = parlr.CPell2d[10:int(maskupP)]
ctl = parlr.CTell2dLen[10:1000]
ct  = parlr.CTell2d[10:1000]
r2  = parlr.grd.r.^2
bin_mids_P = (parlr.grd.deltk*1.5):(parlr.grd.deltk):maskupP
bin_mids_T = (parlr.grd.deltk*1.5):(parlr.grd.deltk):1000

"""
Get Set the burning/thinning and detect the number of parallel runs
------------------------------------------
"""
krang = 551:100:2501  #  burning/thinning

jobs = Int[] # detect the number of parallel runs
if isdir("$simdir/job1")
	push!(jobs, 1)
else
	push!(jobs, 2)
	while isdir("$simdir/job$(jobs[end]+1)")
		push!(jobs,jobs[end]+1 )
	end
end
jobs

"""
Define functions which compute power and average over l - bins
------------------------------------------
"""
function binave(fk::Matrix, kmag::Matrix, bin_mids::Range)
fpwr = Array(Complex{Float64}, length(bin_mids))
fill!(fpwr, -1.0)
rtcuts  = collect(bin_mids +  step(bin_mids) / 2)
lftcuts = collect(bin_mids -  step(bin_mids) / 2)
	lftcuts[1] = 0.1 # extend the left boundary all the way down, but not zero
	for i in 1:length(bin_mids)
		ibin = lftcuts[i] .<= kmag .< rtcuts[i]
		fpwr[i] = fk[ibin] |> mean
	end
	fpwr
end
function binpower(fk::Matrix, kmag::Matrix, bin_mids::Range)
	fpwr = Array(Float64, length(bin_mids))
	fill!(fpwr, -1.0)
	rtcuts  = collect(bin_mids +  step(bin_mids) / 2)
	lftcuts = collect(bin_mids -  step(bin_mids) / 2)
	lftcuts[1] = 0.1 # extend the left boundary all the way down, but not zero
	for i in 1:length(bin_mids)
		ibin = lftcuts[i] .<= kmag .< rtcuts[i]
		fpwr[i] = fk[ibin] |> abs2 |> mean
	end
	fpwr
end

"""
Load the simulation runs
------------------------------------------"""
#-- simulation truths
phix   = readcsv("$simdir/phix.csv")
tildex = readcsv("$simdir/tildetx_lr.csv")
tx_lr  = readcsv("$simdir/tx_lr.csv")
phik   = fft2(phix, parlr)
tildek = fft2(tildex, parlr)
tk_lr  = fft2(tx_lr, parlr)
ytx    = readcsv("$simdir/ytx.csv")
qex    = readcsv("$simdir/qex_lr.csv")
x      = readcsv("$simdir/x.csv")
maskx    = readcsv("$simdir/maskvarx.csv")
maskbool = maskx .== Inf

# initialize  cross correlation
phik_cov_smpl = Array{Float64,1}[]
tildek_cov_smpl = Array{Float64,1}[]
tk_cov_smpl = Array{Float64,1}[]

# initialize  spectral coverage
phik_pwr_smpl = Array{Float64,1}[]
tildek_pwr_smpl = Array{Float64,1}[]
tk_pwr_smpl = Array{Float64,1}[]

# initialize  1-d slices
propslice = 0.75 # btwn 0 and 1
phix_slice_samples = Array{Float64,1}[]
tildex_slice_samples = Array{Float64,1}[]
tx_slice_samples = Array{Float64,1}[]

# initialize  average
phix_sum = zero(phix)
tildetx_lr_sum = zero(phix)
tx_lr_sum = zero(phix)

runs = 0
for r in jobs, k in krang
	if isfile("$simdir/job$r/phix_curr_$k.csv")

		phix_curr = readcsv("$simdir/job$r/phix_curr_$k.csv")
		phik_curr = fft2(phix_curr, parlr)
		tildetx_curr = readcsv("$simdir/job$r/tildetx_lr_curr_$k.csv")
		tildetk_curr = fft2(tildetx_curr, parlr)
		tx_curr = readcsv("$simdir/job$r/tx_lr_curr_$k.csv")
		tk_curr = fft2(tx_curr, parlr)

		# spectral coverage
		push!(phik_pwr_smpl,   binpower(r2    .* phik_curr, parlr.grd.r, bin_mids_P))
		push!(tildek_pwr_smpl, binpower(√(r2) .* tildetk_curr, parlr.grd.r, bin_mids_T))
		push!(tk_pwr_smpl,     binpower(√(r2) .* tk_curr, parlr.grd.r, bin_mids_T))

		# cross correlation
		pnxt   = binave(  r2 .* phik_curr .* conj(phik), parlr.grd.r, bin_mids_P)
		pnxt ./= binpower(sqrt(r2) .* phik_curr, parlr.grd.r, bin_mids_P) |> sqrt
		pnxt ./= binpower(sqrt(r2) .* phik, parlr.grd.r, bin_mids_P) |> sqrt
		push!(phik_cov_smpl, real(pnxt))

		tnxt   = binave(  r2 .* tildetk_curr .* conj(tildek), parlr.grd.r, bin_mids_T)
		tnxt ./= binpower(sqrt(r2) .* tildetk_curr, parlr.grd.r, bin_mids_T) |> sqrt
		tnxt ./= binpower(sqrt(r2) .* tildek, parlr.grd.r, bin_mids_T) |> sqrt
		push!(tildek_cov_smpl, real(tnxt))

		tnxt   = binave(  r2 .* tk_curr .* conj(tk_lr), parlr.grd.r, bin_mids_T)
		tnxt ./= binpower(sqrt(r2) .* tk_curr, parlr.grd.r, bin_mids_T) |> sqrt
		tnxt ./= binpower(sqrt(r2) .* tk_lr, parlr.grd.r, bin_mids_T) |> sqrt
		push!(tk_cov_smpl, real(tnxt))

		# 1-d slices
		push!(phix_slice_samples, phix_curr[int(end*propslice),:][:])
		push!(tildex_slice_samples, tildetx_curr[int(end*propslice),:][:])
		push!(tx_slice_samples, tx_curr[int(end*propslice),:][:])

		# averages
		phix_sum += phix_curr
		tildetx_lr_sum += tildetx_curr
		tx_lr_sum += tx_curr

	runs += 1
	end
end

# 1-d slices
phix_slice_samples   = hcat(phix_slice_samples...)
tildex_slice_samples = hcat(tildex_slice_samples...)
tx_slice_samples     = hcat(tx_slice_samples...)
phix_true_slice      = phix[int(end*propslice),:][:]
tildex_true_slice    = tildex[int(end*propslice),:][:]
tx_true_slice        = tx_lr[int(end*propslice),:][:]
phix_sum_slice       = phix_sum[int(end*propslice),:][:]
tildetx_lr_sum_slice = tildetx_lr_sum[int(end*propslice),:][:]
tx_lr_sum_slice      = tx_lr_sum[int(end*propslice),:][:]
qex_slice  = qex[int(end*propslice),:][:]
x_slice    = x[int(end*propslice),:][:]
varx_slice = maskx[int(end*propslice),:][:]
isempty(x_slice[varx_slice.==Inf]) || (maskmin=minimum(x_slice[varx_slice.==Inf]); maskmax = maximum(x_slice[varx_slice.==Inf]))

# cross correlation
# squish to matrix
phik_cov_smpl   = hcat(phik_cov_smpl...)
tildek_cov_smpl = hcat(tildek_cov_smpl...)
tk_cov_smpl     = hcat(tk_cov_smpl...)
# then split by the rows
phb_cov          = Array{Float64,1}[vec(phik_cov_smpl[k,:])   for k=1:size(phik_cov_smpl,   1)]
tib_cov          = Array{Float64,1}[vec(tildek_cov_smpl[k,:]) for k=1:size(tildek_cov_smpl, 1)]
tib_unlensed_cov = Array{Float64,1}[vec(tk_cov_smpl[k,:]) for k=1:size(tk_cov_smpl, 1)]

# spectral coverage
# squish to matrix
phik_pwr_smpl   = hcat(phik_pwr_smpl...)
tildek_pwr_smpl = hcat(tildek_pwr_smpl...)
tk_pwr_smpl     = hcat(tk_pwr_smpl...)
# then split by the rows
phb_pwr          = Array{Float64,1}[phik_pwr_smpl[k,:][:]   for k=1:size(phik_pwr_smpl,   1)]
tib_pwr          = Array{Float64,1}[tildek_pwr_smpl[k,:][:] for k=1:size(tildek_pwr_smpl, 1)]
tib_unlensed_pwr = Array{Float64,1}[tk_pwr_smpl[k,:][:] for k=1:size(tk_pwr_smpl, 1)]

# here is the truth
phb_truth          = binpower(r2 .* phik, parlr.grd.r, bin_mids_P)
tib_truth          = binpower(sqrt(r2) .* tildek, parlr.grd.r, bin_mids_T)
tib_unlensed_truth = binpower(sqrt(r2) .* tk_lr, parlr.grd.r, bin_mids_T)

"""
Cooling schedule version 2
------------------------------------------
"""
function repstretch(start, stop, rep::Int, total::Int)
    tdr = int(total / rep)
    return vec(linspace(start, stop, tdr)' .* ones(rep, tdr))
end
lcool = repstretch(2.0, 1.2maskupC, 75, 1000) |> int
λcool = Float64[]
d2x   = parlr.grd.deltx * parlr.grd.deltx
d2k   = parlr.grd.deltk * parlr.grd.deltk
delt0 = 1 / d2k
barNx = 0.99 * parlr.nugget_at_each_pixel
barNk = barNx * d2x # spectral density in fourier

ls = 1:4000
semilogy(ls, parlr.CTell2d[ls], color = "k",  linewidth=1.5, linestyle="-", label=L"$C_l^{TT}$")
annotate(L"$C_l^{TT}$",
        xy=(150, parlr.CTell2d[150]), xycoords="data",
        xytext=(35, 35), textcoords="offset points", fontsize=16,
        arrowprops={:arrowstyle=>"->"})
for cntr in [1, 100, 300, 500, 700, 900]
	λσ2dy = (cntr > 8000.0) ? barNk : max(barNk, parlr.CTell2d[lcool[cntr]])
	semilogy(ls, zeros(ls) + λσ2dy, color = "k",  linewidth=1.5, linestyle=":")
	annotate(latexstring("\$\\lambda_{$cntr}\\bar{\\sigma}^2dy\$"),
        xy=(lcool[cntr], λσ2dy), xycoords="data",
        xytext=(20, 7), textcoords="offset points", fontsize=16)
end
semilogy(ls, zeros(ls) + barNk, color = "k",  linewidth=1.5, linestyle="--", label="noise level")
annotate(L"homogeneous noise spectral density = $\bar{\sigma}^2dy$",
        xy=(1540, barNk), xycoords="data",
        xytext=(-100, -40), textcoords="offset points", fontsize=14,
        arrowprops={:arrowstyle=>"->"})
xlabel("wavenumber", fontsize=14)
ylabel("spectral density", fontsize=14)
xlim(0, 3700)
ylim(1e-7, 1e+4)
savefig(joinpath(savepath, "figure4.pdf"), dpi=300, bbox_inches="tight", pad_inches=0.1)

"""
True phix
------------------------------------------
"""
imshow(
  phix,
  interpolation = "nearest",
  vmin=minimum(phix),
  vmax=maximum(phix),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
plt.xlabel("degrees")
plt.ylabel("degrees")
savefig(joinpath(savepath, "figure5a.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Posterior mean for phix
------------------------------------------
"""
imshow(
  phix_sum/runs,
  interpolation = "nearest",
  vmin=minimum(phix),
  vmax=maximum(phix),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
plt.xlabel("degrees")
plt.ylabel("degrees")
savefig(joinpath(savepath, "figure5b.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Quadratic estimate
------------------------------------------
"""
imshow(
	qex,
  interpolation = "nearest",
  vmin=minimum(phix),
  vmax=maximum(phix),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
plt.xlabel("degrees")
plt.ylabel("degrees")
savefig(joinpath(savepath, "figure5c.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Data
------------------------------------------
"""
imshow(
  ytx,
  interpolation = "nearest",
  vmin=minimum(tildex),
  vmax=maximum(tildex),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
xlabel("degrees")
ylabel("degrees")
title("Data")
savefig(joinpath(savepath, "figure6a.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Poster mean tx
------------------------------------------
"""
imshow(
  tx_lr_sum/runs,
  interpolation = "nearest",
  vmin=minimum(tildex),
  vmax=maximum(tildex),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
plt.xlabel("degrees")
plt.ylabel("degrees")
colorbar()
title(L"$E(T(x)|data)$")
savefig(joinpath(savepath, "figure6b.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
tx - tildex
------------------------------------------
"""
imshow(
  (tx_lr - tildex) .* (1- maskbool)  ,
  interpolation = "nearest",
  vmin=0.8*minimum( (tx_lr - tildex) .* (1- maskbool) ),
  vmax=0.8*maximum( (tx_lr - tildex) .* (1- maskbool) ),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
plt.xlabel("degrees")
plt.ylabel("degrees")
colorbar(format="%.e")
title(L"$T(x) - \widetilde T(x)$")
savefig(joinpath(savepath, "figure6d.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
tx - E(tx|data)
------------------------------------------
"""
imshow(
   (tx_lr - tx_lr_sum/runs) .* (1- maskbool) ,
  interpolation = "nearest",
  vmin=0.8*minimum( (tx_lr - tildex) .* (1- maskbool) ),
  vmax=0.8*maximum( (tx_lr - tildex) .* (1- maskbool) ),
  origin="lower",
  extent=(180/pi)*[minimum(x), maximum(x),minimum(x), maximum(x)]
)
#title(latexstring("\$abs(T(x) - E(T(x)|data))\$"))
plt.xlabel("degrees")
plt.ylabel("degrees")
title(L"$T(x) - E(T(x)|data)$")
savefig(joinpath(savepath, "figure6c.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Slices of phix
------------------------------------------
"""
#=
The slices are done at  degree mark 12.7
(180/pi) * parlr.grd.y[int(end*propslice),:]
=#
plot((180/pi)*x_slice, phix_slice_samples[:,1], color="blue", alpha=0.15, label=L"posterior samples of $\phi(x)$")
plot((180/pi)*x_slice, phix_slice_samples[:,2:end], color="blue", alpha=0.15)
plot((180/pi)*x_slice, phix_true_slice, color = "red",  linewidth=1.5, linestyle="-", label=L"simulation truth $\phi(x)$")
annotate("mask",
         xy=((180/pi)*maskmin, 0), xycoords="data",
         xytext=(-50, 50), textcoords="offset points", fontsize=16,
         arrowprops={:arrowstyle=>"->"})
xlabel("degrees")
axis("tight")
axvspan((180/pi)*maskmin, (180/pi)*maskmax,  facecolor="0.5", alpha=0.3, label="mask region")
legend(loc="upper left", fontsize = 12)
ticklabel_format(style="sci", axis="y", scilimits=(0,0))
savefig(joinpath(savepath, "figure7a.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Slice of tx
------------------------------------------
"""
lper = 0.2
rper = 0.5
zoom = round(lper*length(x_slice)):round(rper*length(x_slice))
plot((180/pi)*x_slice[zoom], tx_slice_samples[zoom,1], color="blue", alpha=0.15, label=L"posterior samples $T(x)$")
plot((180/pi)*x_slice[zoom], tx_slice_samples[zoom,2:end], color="blue", alpha=0.15)
plot((180/pi)*x_slice[zoom], tx_true_slice[zoom], color = "r",  linewidth=1.4, linestyle="-", label=L"simulation truth $T(x)$")
axis("tight")
axvspan((180/pi)*maskmin, (180/pi)*maskmax,  facecolor="0.5", alpha=0.3, label="mask region")
annotate("mask", xy=((180/pi)*maskmin, 0), xycoords="data", xytext=(-50, -50), textcoords="offset points", fontsize=16, arrowprops={:arrowstyle=>"->"})
xlabel("degrees")
axis("tight")
legend(loc="upper left", fontsize = 12)
ticklabel_format(style="sci", axis="y", scilimits=(0,0))
savefig(joinpath(savepath, "figure7b.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Spectral coverage for phi
------------------------------------------
"""
plot(elp, elp.^4 .* cpl / 4, "-k", label = L"$l^4C_l^{\phi\phi}/4$")
plot(bin_mids_P, phb_truth / dirac_0 / 4, "or", label = L"$l^4 |\phi_l|^2/(4\delta_0)$")
rtcuts  = collect(bin_mids_P +  step(bin_mids_P) / 2)
lftcuts = collect(bin_mids_P -  step(bin_mids_P) / 2)
lftcuts[1] = 0.1 # extend the left boundary all the way down, but not zero
errorbar(
	bin_mids_P,
	map(median,  phb_pwr / dirac_0 / 4),
	xerr = Array{Float64,1}[bin_mids_P-lftcuts, rtcuts-bin_mids_P],
	yerr = Array{Float64,1}[map(x-> median(x)-quantile(x,0.025), phb_pwr / dirac_0 / 4), map(x-> quantile(x,0.975)-median(x), phb_pwr / dirac_0 / 4)],
	fmt="*b",
	label = L"95% posterior region for $l^4 |\phi_l|^2/ (4 \delta_0)$"
)
plot(collect(bin_mids_P), zero(bin_mids_P), ":k")
xlabel("wavenumber")
legend()
ticklabel_format(style="sci", axis="y", scilimits=(0,0))
axis("tight")
savefig(joinpath(savepath, "figure8a.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Spectral coverage for T
------------------------------------------
"""
plot(elt,  elt.^2 .* ct, "-k", label = L"$l^2C_l^{ TT}$")
plot(bin_mids_T, tib_unlensed_truth / dirac_0 , "or", label =  L"$l^2 |T_l|^2/ \delta_0$")
rtcuts  = collect(bin_mids_T +  step(bin_mids_T) / 2)
lftcuts = collect(bin_mids_T -  step(bin_mids_T) / 2)
lftcuts[1] = 0.1 # extend the left boundary all the way down, but not zero
errorbar(
	bin_mids_T,
	map(median,  tib_unlensed_pwr /  dirac_0),
	xerr = Array{Float64,1}[bin_mids_T-lftcuts, rtcuts-bin_mids_T],
	yerr = Array{Float64,1}[map(x-> median(x)-quantile(x,0.025), tib_unlensed_pwr /  dirac_0), map(x-> quantile(x,0.975)-median(x), tib_unlensed_pwr /  dirac_0)],
	fmt="*b",
	label = L"95\% posterior region for $l^2 |T_l|^2/ \delta_0$"
)
plot(collect(bin_mids_T), zero(bin_mids_T), ":k")
xlabel("wavenumber")
legend()
axis("tight")
ticklabel_format(style="sci", axis="y", scilimits=(0,0))
savefig(joinpath(savepath, "figure8b.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
Cross correlation for phi and T
------------------------------------------
"""
rtcuts  = collect(bin_mids_P +  step(bin_mids_P) / 2)
lftcuts = collect(bin_mids_P -  step(bin_mids_P) / 2)
lftcuts[1] = 0.1 # extend the left boundary all the way down, but not zero
errorbar(
	bin_mids_P,
	map(median,  phb_cov),
	xerr = Array{Float64,1}[bin_mids_P-lftcuts, rtcuts-bin_mids_P],
	yerr = Array{Float64,1}[map(x-> median(x)-quantile(x,0.025), phb_cov),
							map(x-> quantile(x,0.975)-median(x), phb_cov)],
	fmt="*b",
	label = L"empirical cross correlation with $\phi_l$"
)
annotate(L"corr between $\phi_l$ and $P(\phi_l|data)$",
         xy=(bin_mids_P[20], map(median,  phb_cov)[20]+0.02), xycoords="data",
         xytext=(-10, 80), textcoords="offset points", fontsize=16,
         arrowprops={:arrowstyle=>"->"})
rtcuts  = collect(bin_mids_T +  step(bin_mids_T) / 2)
lftcuts = collect(bin_mids_T -  step(bin_mids_T) / 2)
lftcuts[1] = 0.1 # extend the left boundary all the way down, but not zero
errorbar(
	bin_mids_T,
	map(median,  tib_unlensed_cov),
	xerr = Array{Float64,1}[bin_mids_T-lftcuts, rtcuts-bin_mids_T],
	yerr = Array{Float64,1}[map(x-> median(x)-quantile(x,0.025), tib_unlensed_cov),
							map(x-> quantile(x,0.975)-median(x), tib_unlensed_cov)],
	fmt="*b",
	label = L"empirical cross correlation with $T_l$"
)
annotate(L"corr between $T_l$ and $P(T_l|data)$",
         xy=(bin_mids_T[30], map(median,  tib_unlensed_cov)[30]-0.02), xycoords="data",
         xytext=(-170, -60), textcoords="offset points", fontsize=16,
         arrowprops={:arrowstyle=>"->"})
#plot(collect(bin_mids_P), zero(bin_mids_P), ":k")
xlabel("wavenumber")
#legend()
axis("tight")
ylim(-0.05,1)
savefig(joinpath(savepath, "figure9a.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

"""
trace plots
------------------------------------------"""
simdir_trc   =  "/Users/ethananderes/Dropbox/BayesLense/simulations/scriptTrack_2657077506/" # contains the simulation
trace_trc   = map(include_string, readcsv("$(simdir_trc)/trace.csv"))

# ----This is the degree reference for trace_trc[3, :]
# [(180/pi) * parlr.grd.x[400,300], (180/pi) * parlr.grd.y[400,300]] #<--- [9.9,13.2]
# ----This is the degree reference for trace_trc[4, :]
# [(180/pi) * parlr.grd.x[50,50], (180/pi) *  parlr.grd.y[50,50]] #<---------- [1.6, 1.6]
plot(vcat(0.0, real(trace_trc[3,2:5001])...),     color = "blue", alpha = 0.5,  linewidth=1.5, linestyle="-", label = L"\phi(x),\, x = (9.9^o, 13.2^o)")
plot(zeros(size(trace_trc,2)) + trace_trc[3,1],  color = "blue", alpha = 1.0,  linewidth=2.5, linestyle="--")
plot(vcat(0.0, real(trace_trc[4,2:5001])...),     color = "green", alpha = 0.5,  linewidth=1.5, linestyle="-",  label = L"\phi(x),\, x = (1.6^o, 1.6^o)")
plot(zeros(size(trace_trc,2)) + trace_trc[4,1],  color = "green", alpha = 1.0,  linewidth=2.5, linestyle="--")
xlabel("Gibbs iteration")
ylabel(L"Gibbs chain values of $\phi(x)$")
ticklabel_format(style="sci", axis="y", scilimits=(0,0))
title("Assessing Gibbs chain convergence")
legend(loc="best", fontsize = 12)
axis("tight")
savefig(joinpath(savepath, "figure10a.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

# ----This is the degree reference for trace_trc[6, :]
# k = [parlr.grd.k1[4, 7],  parlr.grd.k2[4, 7]] #<----  [126.56,63.28]
# wavenumber = norm(k) #<--- 141.50
plot(vcat(0.0, real(trace_trc[6,2:5001])...),     color = "blue", alpha = 0.5,  linewidth=1.5, linestyle="-", label = L"real(\phi_l),\,l = (126.56,63.28)")
plot(zeros(size(trace_trc,2)) + real(trace_trc[6,1]), color = "blue", alpha = 1.0,  linewidth=2.5, linestyle="--")
plot(vcat(0.0, imag(trace_trc[6,2:5001])...),  color = "green", alpha = 0.5,  linewidth=1.5, linestyle="-", label = L"imag(\phi_l),\,l = (126.56,63.28)")
plot(zeros(size(trace_trc,2)) + imag(trace_trc[6,1]),  color = "green", alpha = 1.0,  linewidth=2.5, linestyle="--")
xlabel("Gibbs iteration")
ylabel(L"Gibbs chain values of $\phi_l$")
title("Assessing Gibbs chain convergence")
ticklabel_format(style="sci", axis="y", scilimits=(0,0))
legend(loc="best", fontsize = 12)
axis("tight")
savefig(joinpath(savepath, "figure10b.pdf"), dpi=300, bbox_inches="tight",  pad_inches=0.1)

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
"""
The submodule `MLLabelUtils.LabelEnc` contains a selection of
popular label encodings:

Strictly binary label encodings:

- [`LabelEnc.FuzzyBinary`](@ref)
- [`LabelEnc.TrueFalse`](@ref)
- [`LabelEnc.ZeroOne`](@ref)
- [`LabelEnc.MarginBased`](@ref)

Multiclass label encodings:

- [`LabelEnc.Indices`](@ref)
- [`LabelEnc.NativeLabels`](@ref)
- [`LabelEnc.OneOfK`](@ref)

Multiclass to binary:

- [`LabelEnc.OneVsRest`](@ref)
"""
module LabelEnc
    import ..MLLabelUtils.LabelEncoding
    import ..MLLabelUtils.BinaryLabelEncoding

    """
        LabelEnc.FuzzyBinary <: LabelEncoding{Any,2,1}

    A vector-based binary label encoding without a specific labeltype.

    Primarily used for fuzzy comparision of binary true targets
    and predicted targets. It basically assumes that the encoding
    is either `TrueFalse`, `ZeroOne`, or `MarginBased` by treating
    all non-negative values as positive outputs.
    """
    struct FuzzyBinary <: BinaryLabelEncoding{Any,1} end

    """
        LabelEnc.TrueFalse <: LabelEncoding{Bool,2,1}

    A vector-based binary label encoding with labeltype `Bool`.

    Represents binary classification-labels as boolean, in which `true`
    represents the positive class, and `false` the negative class.
    """
    struct TrueFalse <: BinaryLabelEncoding{Bool,1} end

    """
        LabelEnc.ZeroOne{T<:Number} <: LabelEncoding{T,2,1}

    A vector-based binary label encoding with a numeric labeltype.

    Represents binary classification-labels as numbers, in which
    `one(T)` represents the positive class, and `zero(T)` the
    negative class.
    """
    struct ZeroOne{T<:Number,R<:Number} <: BinaryLabelEncoding{T,1}
        cutoff::R
        function ZeroOne{T,R}(cutoff::R = R(0.5)) where {T,R}
            @assert 0 <= cutoff <= 1
            new{T,R}(cutoff)
        end
    end
    ZeroOne(t::Type{T} = Float64, cutoff::R = 0.5) where {T<:Number,R<:Number} = ZeroOne{T,R}(cutoff)
    ZeroOne(cutoff::R) where {R<:Number} = ZeroOne(R, cutoff)

    """
        LabelEnc.MarginBased{T<:Number} <: LabelEncoding{T,2,1}

    A vector-based binary label encoding with a numeric labeltype.

    Represents binary classification-labels as numbers, in which
    `one(T)` represents the positive class, and `-one(T)` the
    negative class.
    """
    struct MarginBased{T<:Number} <: BinaryLabelEncoding{T,1} end
    MarginBased(::Type{T} = Float64) where {T<:Number} = MarginBased{T}()

    """
        LabelEnc.OneVsRest{T} <: LabelEncoding{T,2,1}

    A vector-based binary label encoding with a arbitrary labeltype.

    A special label-encoding in that it only uses a positive
    label to determine which class some element belongs to.
    If some value matches the specified positive label, then it is
    considered positive, otherwise it is considered negative.

    That said `OneVsRest` requires a negative label to be specified
    in order to be able to denote it somehow. Note, however,
    that the specified negative label is purely for asthetic
    reasons and is not used to determine class membership
    """
    struct OneVsRest{T} <: BinaryLabelEncoding{T,1}
        poslabel::T
        neglabel::T
    end
    OneVsRest(poslabel::Bool) = OneVsRest{Bool}(poslabel, !poslabel)
    OneVsRest(poslabel::String) = OneVsRest{String}(poslabel, "not_$(poslabel)")
    OneVsRest(poslabel::Symbol) = OneVsRest{Symbol}(poslabel, Symbol(:not_, poslabel))
    OneVsRest(poslabel::T) where {T<:Number} = OneVsRest{T}(poslabel, ifelse(poslabel == 0, poslabel+T(1), T(0)))

    """
        LabelEnc.Indices{T,K<:Number} <: LabelEncoding{T,K,1}

    A vector-based multi-label encoding with a numeric labeltype.

    Represents the class-labels as indices starting from `T(1)` up to
    (including) `T(K)`.
    In a binary setting `T(1)` corresponds to te positive class and
    `T(2)` corresponds to the negative class.
    """
    struct Indices{T<:Number,K} <: LabelEncoding{T,K,1}
        function Indices{T,K}() where {T,K}
            typeof(K) <: Int || throw(TypeError(:Indices,"constructor when checking typeof(K)",Type{Int},typeof(K)))
            new{T,K}()
        end
    end
    Indices(::Type{T}, ::Type{Val{K}}) where {T,K} = Indices{T,K}()
    Indices(::Type{Val{K}}) where {K} = Indices(Int,Val{K})
    Indices(::Type{T}, K::Number) where {T} = Indices(T,Val{Int(K)})
    Indices(K::Number) = Indices(Int,Val{Int(K)})

    """
        LabelEnc.OneOfK{T,K<:Number} <: LabelEncoding{T,K,2}

    A matrix-based multi-label encoding with a numeric labeltype.

    Represents the class labels in a one-hot encoding scheme.
    That means a matrix in which for every observation one of
    `K` elements is set to 1.
    """
    struct OneOfK{T<:Number,K} <: LabelEncoding{T,K,2}
        function OneOfK{T,K}() where {T,K}
            typeof(K) <: Int || throw(TypeError(:OneOfK,"constructor when checking typeof(K)",Type{Int},typeof(K)))
            new{T,K}()
        end
    end
    OneOfK(::Type{T}, ::Type{Val{K}}) where {T,K} = OneOfK{T,K}()
    OneOfK(::Type{Val{K}}) where {K} = OneOfK(Int,Val{K})
    OneOfK(::Type{T}, K::Number) where {T} = OneOfK(T,Val{Int(K)})
    OneOfK(K::Number) = OneOfK(Int,Val{Int(K)})

    """
        LabelEnc.NativeLabels{T,K,F} <: LabelEncoding{T,K,1}

    A vector-based multi-label encoding with a arbitrary labeltype.

    Represents arbitrary class labels by storing the possible values
    as a member variable.

    Constructor `LabelEnc.NativeLabels([fallbacklabel], classlabels)`
     - `classlabels` is a collection of labels that are allowed by the
       encoding. Each one representing a class.
     - `fallbacklabel` is an optional first argument.
         - It is either:
            - a value from `classlabels` which is used to encode any observed
              label that does not occur in `classlabels`
            - or; a function which takes as its input a label which does not
              occur in `classlabels` and returns one that does.
         - This can be used for example to handle out of vocabulary words in NLP
           as in `LabelEnc.NativeLabels("<OOV>", [red", "green", blue", ..., "<OOV>"])`
         - Or to standardize inputs: as in
           ```
            LabelEnc.NativeLabels([..., v"0.5.2", v"0.6.0", v"0.6.1", ....]) do oov
                # Discard build and prerelease fields, `v"0.7.0-beta.18"`->`v"0.7.0"`
                VersionNumber(oov.major, oov.minor, oov.patch)
            end
           ```

    In a binary setting the first element of the stored class labels
    represents the positive label and the second element the negative
    label.
    In a multiclass setting, the labels will `convertenc` to
    `LabelEnc.Indices`  using their position in the `label` argument
    to the constructor.
    """
    struct NativeLabels{T,K,F} <: LabelEncoding{T,K,1}
        getfallbacklabel::F
        label::Vector{T}
        invlabel::Dict{T,Int}
        function NativeLabels{T,K,F}(getfallbacklabel::F, label::Vector{T}) where {T,K,F<:Function}
            typeof(K) <: Int || throw(TypeError(:NativeLabels,"constructor when checking typeof(K)",Type{Int},typeof(K)))
            @assert length(label) == length(unique(label)) == K
            new{T,K,F}(getfallbacklabel, label, Dict(zip(label,1:K)))
        end
    end
    const default_getfallbacklabel = identity # by default leave label as is, so parent can change

    NativeLabels{T,K}(getfallbacklabel::F, label::Vector{T}) where {T,K,F<:Function} = NativeLabels{T,K,F}(getfallbacklabel, label)
    NativeLabels{T,K}(fallbacklabel, label::Vector{T}) where {T,K} = NativeLabels{T,K}(oov->fallbacklabel, label)
    NativeLabels{T,K}(label::Vector{T}) where {T,K} = NativeLabels{T,K}(default_getfallbacklabel, label)

    NativeLabels{T,K}(getfallbacklabel::F, label::AbstractVector{T}) where {F,T,K} = NativeLabels{T,K}(getfallbacklabel, collect(label))
    NativeLabels{T,K}(fallbacklabel::T, label::AbstractVector{T}) where {T,K} = NativeLabels{T,K}(oov->fallbacklabel, label)
    NativeLabels{T,K}(label::AbstractVector{T}) where {T,K} = NativeLabels{T,K}(default_getfallbacklabel, label)

    NativeLabels(getfallbacklabel::Function, label::AbstractVector{T}, ::Type{Val{K}})  where {T,K} = NativeLabels{T,K}(getfallbacklabel, label)
    NativeLabels(fallbacklabel, label::AbstractVector{T}, ::Type{Val{K}})  where {T,K} = NativeLabels(oov->fallbacklabel, label, Val{K})
    NativeLabels(label::AbstractVector, ::Type{Val{K}})  where {K} = NativeLabels(default_getfallbacklabel, label, Val{K})

    NativeLabels(getfallbacklabel, label) = NativeLabels(getfallbacklabel, label, Val{length(label)})
    NativeLabels(label) = NativeLabels(default_getfallbacklabel, label)

    Base.hash(a::NativeLabels, h::UInt) = hash(a.getfallbacklabel, hash(a.label, hash(:NativeLabels, h)))
    Base.:(==)(a::NativeLabels, b::NativeLabels) = isequal(a.label, b.label) && isequal(a.getfallbacklabel, b.getfallbacklabel)

end # submodule

standardize_label(lbl, lm::LabelEnc.NativeLabels) = haskey(lm.invlabel, lbl) ? lbl : lm.getfallbacklabel(lbl)

_ambiguous() = throw(ArgumentError("Can't infer the label meaning because argument types or values are ambiguous. Please specify the desired LabelEncoding manually."))

# Query the index
label2ind(lbl, lm::BinaryLabelEncoding) = ifelse(isposlabel(lbl, lm), 1, 2)
label2ind(lbl::Union{Number,T}, lm::LabelEnc.Indices{T}) where {T} = Int(lbl)
label2ind(lbl::Union{Number,T}, lm::LabelEnc.OneOfK{T}) where {T} = Int(lbl)
label2ind(lbl::AbstractVector, lm::LabelEnc.OneOfK) = argmax(lbl)
function label2ind(lbl, lm::LabelEnc.NativeLabels)
    std_lbl = standardize_label(lbl, lm)
    Int(lm.invlabel[std_lbl])
end

# Query the label
ind2label(i::Integer, lm::BinaryLabelEncoding) = ifelse(i == 1, ind2label(Val{1},lm), ind2label(Val{2},lm))
ind2label(i::Integer, ::LabelEnc.Indices{T}) where {T} = T(i)
ind2label(i::Integer, ::LabelEnc.OneOfK{T,K}) where {T,K} = (x = zeros(T,K); x[i] = T(1); x)
ind2label(i::Integer, lm::LabelEnc.NativeLabels) = lm.label[Int(i)]

ind2label(::Type{Val{1}}, ::LabelEnc.TrueFalse) = true
ind2label(::Type{Val{2}}, ::LabelEnc.TrueFalse) = false
ind2label(::Type{Val{1}}, ::LabelEnc.ZeroOne{T}) where {T} = T(1)
ind2label(::Type{Val{2}}, ::LabelEnc.ZeroOne{T}) where {T} = T(0)
ind2label(::Type{Val{1}}, ::LabelEnc.MarginBased{T}) where {T} = T(1)
ind2label(::Type{Val{2}}, ::LabelEnc.MarginBased{T}) where {T} = -T(1)
ind2label(::Type{Val{1}}, ovr::LabelEnc.OneVsRest) = ovr.poslabel
ind2label(::Type{Val{2}}, ovr::LabelEnc.OneVsRest) = ovr.neglabel
ind2label(::Type{Val{K}}, ::LabelEnc.Indices{T}) where {T,K} = T(K)
ind2label(::Type{Val{I}}, ::LabelEnc.OneOfK{T,K}) where {T,I,K} = (x = zeros(T,K); x[I] = T(1); x)
ind2label(::Type{Val{K}}, lm::LabelEnc.NativeLabels) where {K} = lm.label[K]

poslabel(lm::LabelEnc.BinaryLabelEncoding) = ind2label(Val{1}, lm)
neglabel(lm::LabelEnc.BinaryLabelEncoding) = ind2label(Val{2}, lm)

label(lm::LabelEnc.BinaryLabelEncoding{T,1}) where {T} = [poslabel(lm), neglabel(lm)]
label(::LabelEnc.Indices{T,K}) where {T,K} = collect(T(1):T(K))
label(::LabelEnc.OneOfK{T,K}) where {T,K} = collect(1:K)
label(lm::LabelEnc.NativeLabels) = lm.label

labeltype(::Type{LabelEnc.ZeroOne}) = Number
labeltype(::Type{LabelEnc.MarginBased}) = Number
labeltype(::Type{LabelEnc.Indices}) = Number
labeltype(::Type{LabelEnc.OneOfK}) = Number

# What it means to be a positive label
isposlabel(value, ::LabelEnc.FuzzyBinary) = _ambiguous()
isposlabel(value, ovr::LabelEnc.OneVsRest) = (value == ovr.poslabel)
isposlabel(value::Bool, ::LabelEnc.FuzzyBinary) = value
isposlabel(value::Bool, ::LabelEnc.TrueFalse)   = value
isposlabel(value::Bool, ::LabelEnc.MarginBased) = throw(MethodError(isposlabel,(value,)))
isposlabel(value::Bool, ::LabelEnc.Indices)     = throw(MethodError(isposlabel,(value,)))
isposlabel(value::Bool, ::LabelEnc.OneOfK)      = throw(MethodError(isposlabel,(value,)))
isposlabel(value::T, ::LabelEnc.FuzzyBinary) where {T<:Number} = (value > T(0))
isposlabel(value::T, zo::LabelEnc.ZeroOne)   where {T<:Number} = (value >= zo.cutoff)
isposlabel(value::T, ::LabelEnc.MarginBased) where {T<:Number} = (sign(value) >= T(0))
isposlabel(value::Number, lm::LabelEnc.Indices{T,2}) where {T} = value == poslabel(lm)
isposlabel(value::Number, lm::LabelEnc.OneOfK{T,2})  where {T} = value == T(1)
isposlabel(value::AbstractVector{<:Number}, lm::LabelEnc.OneOfK{T,2}) where {T} = argmax(value) == 1
isposlabel(value, lm::LabelEnc.NativeLabels{T,2}) where {T} = standardize_label(value, lm) == poslabel(lm)

# What it means to be a negative label
isneglabel(value, ::LabelEnc.FuzzyBinary) = _ambiguous()
isneglabel(value, ovr::LabelEnc.OneVsRest) = (value != ovr.poslabel)
isneglabel(value::Bool, ::LabelEnc.FuzzyBinary) = !value
isneglabel(value::Bool, ::LabelEnc.TrueFalse)   = !value
isneglabel(value::Bool, ::LabelEnc.MarginBased) = throw(MethodError(isneglabel,(value,)))
isneglabel(value::Bool, ::LabelEnc.Indices)     = throw(MethodError(isneglabel,(value,)))
isneglabel(value::Bool, ::LabelEnc.OneOfK)      = throw(MethodError(isneglabel,(value,)))
isneglabel(value::T, ::LabelEnc.FuzzyBinary) where {T<:Number} = (value <= T(0))
isneglabel(value::T, zo::LabelEnc.ZeroOne)   where {T<:Number} = (value < zo.cutoff)
isneglabel(value::T, ::LabelEnc.MarginBased) where {T<:Number} = (sign(value) == -1)
isneglabel(value::Number, lm::LabelEnc.Indices{T,2}) where {T} = value == neglabel(lm)
isneglabel(value::Number, lm::LabelEnc.OneOfK{T,2})  where {T} = value == T(2)
isneglabel(value::AbstractVector{<:Number}, lm::LabelEnc.OneOfK{T,2}) where {T} = argmax(value) == 2
isneglabel(value, lm::LabelEnc.NativeLabels{T,2}) where {T} = standardize_label(value, lm) == neglabel(lm)

# Check if the encoding is appropriate
islabelenc(targets::AbstractArray, args...) = false
islabelenc(targets::AbstractVector{T}, ::LabelEnc.FuzzyBinary)       where {T<:Union{Number,Bool}} = true
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.FuzzyBinary}) where {T<:Union{Number,Bool}} = true
islabelenc(targets::AbstractVector{T}, ::LabelEnc.ZeroOne{T})        where {T<:Union{Number,Bool}} = all(x == 0 || x == 1 for x in targets)
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.ZeroOne{T}})  where {T<:Union{Number,Bool}} = all(x == 0 || x == 1 for x in targets)
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.ZeroOne})     where {T<:Union{Number,Bool}} = all(x == 0 || x == 1 for x in targets)
islabelenc(targets::AbstractVector{T}, ::LabelEnc.MarginBased{T})       where {T<:Number} = all(x == -1 || x == 1 for x in targets)
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.MarginBased{T}}) where {T<:Number} = all(x == -1 || x == 1 for x in targets)
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.MarginBased})    where {T<:Number} = all(x == -1 || x == 1 for x in targets)
islabelenc(targets::AbstractVector{Bool}, ::LabelEnc.TrueFalse) = true
islabelenc(targets::AbstractVector{Bool}, ::Type{LabelEnc.TrueFalse}) = true
islabelenc(targets::AbstractVector{T}, lm::LabelEnc.OneVsRest{T})   where {T}           = any(x == lm.poslabel for x in targets)
islabelenc(targets::AbstractVector{T}, ::LabelEnc.Indices{T,K})     where {T<:Number,K} = all(0 < x <= T(K) && isinteger(x) for x in targets)
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.Indices{T}}) where {T<:Number}   = all(0 < x && isinteger(x) for x in targets)
islabelenc(targets::AbstractVector{T}, ::Type{LabelEnc.Indices})    where {T<:Number}   = all(0 < x && isinteger(x) for x in targets)

function islabelenc(targets::AbstractVector, lm::LabelEnc.NativeLabels; strict=true)
    if strict
        all(haskey(lm.invlabel, x) for x in targets)
    else
        all(haskey(lm.invlabel, standardize_label(x, lm)) for x in targets)
    end
end

function islabelenc(targets::AbstractMatrix{<:Union{Bool,Number}}, lm; obsdim = LearnBase.default_obsdim(targets))
    islabelenc(targets, lm, convert(LearnBase.ObsDimension,obsdim))
end

function islabelenc(targets::AbstractMatrix{<:Union{Bool,Number}}, lm::LabelEnc.OneOfK, obsdim)
    islabelenc(targets, typeof(lm), obsdim)
end

function islabelenc(targets::AbstractMatrix{T}, ::Type{LabelEnc.OneOfK{T}}, obsdim) where {T<:Union{Bool,Number}}
    islabelenc(targets, LabelEnc.OneOfK, obsdim)
end

function islabelenc(targets::AbstractMatrix{T}, ::Type{LabelEnc.OneOfK{T,K}}, ::Union{ObsDim.Last,ObsDim.Constant{2}}) where {T<:Union{Bool,Number},K}
    k, n = size(targets)
    ifelse(k != K, false, islabelenc(targets, LabelEnc.OneOfK, ObsDim.Last()))
end

function islabelenc(targets::AbstractMatrix{T}, ::Type{LabelEnc.OneOfK{T,K}}, ::ObsDim.First) where {T<:Union{Bool,Number},K}
    n, k = size(targets)
    ifelse(k != K, false, islabelenc(targets, LabelEnc.OneOfK, ObsDim.First()))
end

function islabelenc(targets::AbstractMatrix{T}, ::Type{LabelEnc.OneOfK}, ::Union{ObsDim.Last,ObsDim.Constant{2}}) where {T<:Union{Bool,Number}}
    k, n = size(targets)
    @inbounds for i in 1:n
        found = false
        for j in 1:k
            tcur = targets[j,i]
            if tcur == T(1)
                if found
                    return false
                end
                found = true
            elseif tcur == T(0)
                # this is fine
            else
                return false
            end
        end
        if !found
            return false
        end
    end
    return true
end

function islabelenc(targets::AbstractMatrix{T}, ::Type{LabelEnc.OneOfK}, ::ObsDim.First) where {T<:Union{Bool,Number}}
    n, k = size(targets)
    @inbounds for i in 1:n
        found = false
        for j in 1:k
            tcur = targets[i,j]
            if tcur == T(1)
                if found
                    return false
                end
                found = true
            elseif tcur == T(0)
                # this is fine
            else
                return false
            end
        end
        if !found
            return false
        end
    end
    return true
end

# Automatic determination of label mode
labelenc(target) = _ambiguous()
labelenc(targets::AbstractVector{Bool}) = LabelEnc.TrueFalse()

function labelenc(targets::AbstractVector)
    lbls = label(targets)
    LabelEnc.NativeLabels(lbls)
end

function labelenc(targets::AbstractVector{T}) where {T<:Number}
    lbls = label(targets)
    if length(lbls) == 1
        if lbls[1] == 0
            LabelEnc.ZeroOne(T)
        elseif lbls[1] == -1
            LabelEnc.MarginBased(T)
        else
            _ambiguous() # could be ZeroOne or MarginBased
        end
    elseif length(lbls) == 2
        mi, ma = extrema(lbls)
        if mi == 0 && ma == 1
            LabelEnc.ZeroOne(T)
        elseif mi == -1 && ma == 1
            LabelEnc.MarginBased(T)
        elseif mi == 1 && ma == 2
            LabelEnc.Indices(T, Val{2})
        else
            LabelEnc.NativeLabels([ma, mi])
        end
    elseif all(x > 0 && isinteger(x) for x in lbls)
        LabelEnc.Indices(T, maximum(lbls))
    elseif length(lbls) > 10 && any(!isinteger(x) for x in lbls)
        @warn("The number of distinct floating point numbers (including at least one that is non-integer!) in the label vector is quite large. Are you sure you want to perform classification and not regression?")
        LabelEnc.NativeLabels(lbls)
    else
        LabelEnc.NativeLabels(lbls)
    end
end

# TODO: Multilabel (Matrix as targets)
function labelenc(targets::AbstractMatrix{<:Number}; obsdim = LearnBase.default_obsdim(targets))
    labelenc(targets, convert(LearnBase.ObsDimension,obsdim))
end

function labelenc(targets::AbstractMatrix{T}, ::Union{ObsDim.Last,ObsDim.Constant{2}}) where {T<:Number}
    LabelEnc.OneOfK(T,size(targets,1))
end

function labelenc(targets::AbstractMatrix{T}, ::ObsDim.First) where {T<:Number}
    LabelEnc.OneOfK(T,size(targets,2))
end

ext: jl


text:
# Modified from Shreyas Kowshik's implementation.

"""
Categorical policy, for discrete action spaces.
"""
mutable struct CategoricalPolicy
    π # Neural network for the policy
    value_net # Neural network for the value function
    solver # Reference to the TRPOSolver or PPOSolver
end

function CategoricalPolicy(solver)
    policy_net = Chain(Dense(solver.state_size, solver.hidden_layer_size, relu; initW=_random_normal, initb=constant_init),
                       Dense(solver.hidden_layer_size, solver.action_size; initW=_random_normal, initb=constant_init),
                       x -> softmax(x))

    value_net = Chain(Dense(solver.state_size, solver.hidden_layer_size ,relu; initW=_random_normal),
                      Dense(solver.hidden_layer_size, solver.hidden_layer_size, relu; initW=_random_normal),
                      Dense(solver.hidden_layer_size, 1; initW=_random_normal))

    return CategoricalPolicy(policy_net, value_net, solver)
end

"""
Diagonal Gaussian policy, for continuous action space.
"""
mutable struct DiagonalGaussianPolicy
    μ # Neural network for the mean of the Gaussian distribution
    logΣ # Neural network for the log standard deviation of the Gaussian distribution
    value_net # Neural network for the value function
    solver # Reference to the TRPOSolver or PPOSolver
end

function DiagonalGaussianPolicy(solver, log_std)
    μ = Chain(Dense(solver.state_size, solver.hidden_layer_size, tanh; initW=_random_normal, initb=constant_init),
              Dense(solver.hidden_layer_size, solver.action_size; initW=_random_normal, initb=constant_init),
              x->tanh.(x),
              x->param(2.0) .* x) # TODO. Make this vector of environment STDs (i.e. scale it to [-x, +x] action bound)

    value_net = Chain(Dense(solver.state_size, solver.hidden_layer_size, tanh; initW=_random_normal),
                      Dense(solver.hidden_layer_size, solver.hidden_layer_size, tanh; initW=_random_normal),
                      Dense(solver.hidden_layer_size, 1; initW=_random_normal))

    logΣ = param(ones(solver.action_size) * log_std)

    return DiagonalGaussianPolicy(μ, logΣ, value_net, solver)
end

"""
Takes in the policy and returns an action based on the current state.
"""
function get_action(policy::CategoricalPolicy, state)
    action_probs = policy.π(state).data
    action_probs = reshape(action_probs, policy.solver.action_size)
    return Distributions.sample(1:policy.solver.action_size, Distributions.Weights(action_probs))
end

function get_action(policy::DiagonalGaussianPolicy, state)
    # Our policy outputs the parameters of a Normal distribution
    μ = reshape(policy.μ(state), policy.solver.action_size)
    log_std = policy.logΣ

    σ² = (exp.(log_std)).^2
    Σ = diagm(0=>σ².data)
    distribution = MvNormal(μ.data, Σ)

    return rand(distribution, policy.solver.action_size)
end

test_action(policy::CategoricalPolicy, state) = get_action(policy, state)
test_action(policy::DiagonalGaussianPolicy, state) = policy.μ(state).data # Use only the mean for prediction

"""
Translate action from neural network policy to either an ASTSeedAction or an ASTSampleAction
"""
function translate_ast_action(sim::GrayBox.Simulation, action, ::Type{ASTSeedAction})
    # Convert to seed: NOTE, not every useful when using PPO---effectively random. You've been warned (try MCTS)!
    return ASTSeedAction(hash_uint32(action...))
end

function translate_ast_action(sim::GrayBox.Simulation, action, ::Type{ASTSampleAction})
    gray_environment = GrayBox.environment(sim)
    environment_sample = GrayBox.EnvironmentSample()
    for (i,k) in enumerate(keys(gray_environment))
        # log-probability from the environment's distributions (not the log_prob from the NN policy)
        logp = logpdf(gray_environment[k], action[i])
        environment_sample[k] = GrayBox.Sample(action[i], logp)
    end
    return ASTSampleAction(environment_sample)
end

"""
Return the log-probability of an action under the current policy parameters.
"""
function log_prob(policy::CategoricalPolicy, states::Array, actions::Array)
    action_probs = policy.π(states)
    actions_one_hot = zeros(policy.solver.action_size, size(action_probs)[end])
    for i in 1:size(action_probs)[end]
        actions_one_hot[actions[:,i][1],i] = 1.0
    end
    return log.(sum((action_probs .+ 1f-5) .* actions_one_hot, dims=1))
end

function log_prob(policy::DiagonalGaussianPolicy, states::Array, actions::Array)
    μ = policy.μ(states)
    σ = exp.(policy.logΣ)
    σ² = σ.^2
    log_probs = broadcast(-, ((actions .- μ).^2)./(2 .* σ²)) .- 0.5*log.(sqrt(2π)) .- log.(σ)
    return log_probs
end

"""
Return the entropy of the policy distribution.
"""
function entropy(policy::CategoricalPolicy, states::Array)
    action_probs = policy.π(states)
    return sum(action_probs .* log.(action_probs .+ 1f-10), dims=1)
end

entropy(policy::DiagonalGaussianPolicy, states::Array) = 0.5 + 0.5 * log(2π) .+ policy.logΣ

"""
kl_params:

`old_log_probs` : CategoricalPolicy

`Array([μ, logΣ])` : DiagonalGaussianPolicy
"""
function kl_divergence(policy::CategoricalPolicy, kl_params, states::Array)
    old_log_probs = hcat(cat(kl_params..., dims=1)...)

    action_probs = policy.π(states)
    log_probs = log.(action_probs)

    log_ratio = log_probs .- old_log_probs
    kl_div = (exp.(old_log_probs)) .* log_ratio

    return -1.0f0 .* sum(kl_div, dims=1)
end

function kl_divergence(policy::DiagonalGaussianPolicy, kl_params, states::Array)
    μ0 = policy.μ(states)
    logΣ0 = policy.logΣ
    μ1 = mus = hcat([kl_params[i][1] for i in 1:length(kl_params)]...)
    logΣ1 = hcat([kl_params[i][2] for i in 1:length(kl_params)]...)

    var0 = exp.(2 .* logΣ0)
    var1 = exp.(2 .* logΣ1)
    pre_sum = 0.5 .* (((μ0 .- μ1).^2 .+ var0) ./ (var1 .+ 1e-8) .- 1.0f0) .+ logΣ1 .- logΣ0
    kl = sum(pre_sum, dims=1)

    return kl
end

get_policy_params(policy::CategoricalPolicy) = Flux.params(policy.π)
get_policy_params(policy::DiagonalGaussianPolicy) = Flux.params(Flux.params(policy.μ)..., Flux.params(policy.logΣ)...)

get_policy_net(policy::CategoricalPolicy) = [policy.π]
get_policy_net(policy::DiagonalGaussianPolicy) = [policy.μ, policy.logΣ]

get_value_params(policy::Union{CategoricalPolicy, DiagonalGaussianPolicy}) = Flux.params(policy.value_net)

get_value_net(policy::Union{CategoricalPolicy, DiagonalGaussianPolicy}) = [policy.value_net]

function get_policy(solver)
    if solver.policy_type == :discrete
        return CategoricalPolicy(solver)
    elseif solver.policy_type == :continuous
        return DiagonalGaussianPolicy(solver, solver.log_std)
    else
        error("Policy type not supported, $(solver.policy_type)")
    end
end

# TODO. Generalize for discrete spaces.
function set_action_size!(solver, mdp::MDP)
    if actiontype(mdp)== ASTSeedAction
        solver.action_size = 1
    elseif actiontype(mdp) == ASTSampleAction
        solver.action_size = length(GrayBox.environment(mdp.sim))
    else
        @warn "solver.action_size not set for actiontype: $(actiontype(mdp))"
    end
    return solver
end


text:
# Performance_measurement.jl
module Performance_measurement

# Performance_measurement related code

end # module


text:
"""
绘制bubble贡献的图
"""

using Plots

"""
bubb的贡献
"""
function bubb(x, y)
    eps_k = x
    neps_kp = y
    lamb = 1e-2
    if abs(eps_k - neps_kp) < 1.e-10
        #如果两个数值比较接近, Pi^{-}和Pi^{+}的公式完全一样，就是第二个能量要加个负号
        # lim (eps_k -> -eps_kp) Pi^{-} =
        # 1/T (e^{eps/T} (-eps/T*e^{eps/T} + eps/T + e^{eps/T} + 1))
        #/ (e^{eps/T} + 1)^3
        bval = eps_k / lamb
        if bval > 25
            #@info "数值不稳定"
            return 0.
        end
        expb = exp(bval)
        num = expb * (-bval * expb + bval + expb + 1)
        den = (1+expb)^3
        d_val = num / den / lamb
    else
        if (eps_k / lamb) > 25
            #@info "数值不稳定"
            num_left = 0.
        else
            #e^{epsilon_k / T}
            exp_k_t = exp(eps_k / lamb)
            num_left = eps_k / lamb * exp_k_t / ((1 + exp_k_t)^2)
        end
        if (neps_kp / lamb) > 25
            #@info "数值不稳定"
            num_righ = 0.
        else
            #e^{-epsilon_{-k+q} / T}
            exp_nkp_t = exp(neps_kp / lamb)
            num_righ = neps_kp / lamb * exp_nkp_t /
            ((1 + exp_nkp_t)^2)
        end
        d_val = (num_left - num_righ) / (eps_k - neps_kp)
    end # end if 小区域的贡献计算完
    return d_val
end

surface(-4:0.1:4, -4:0.1:4, bubb)

savefig("bubb_contrib2.png")


text:
using PEGParser
import PEGParser.uncached_parse
using Compat

type RepeatedRule <: Rule
  name::String
  rule
  action

  RepeatedRule(name::String, rule) = new(name, rule, no_action)
end

function repeated(name::String, pdata::ParserData, args::Array)
  rule = parseDefinition("$(name)_rule", args[1], pdata)
  return RepeatedRule(name, rule)
end

function uncached_parse(grammar::Grammar, rule::RepeatedRule, text::String, pos::Int64, cache)
  (ast, pos, error) = parse(grammar, rule.rule, text, pos, cache)
  if error !== nothing
    return (nothing, pos, ParseError("No match (RepeatedRule)", pos))
  end

  count = 1
  while error === nothing
    (ast, npos, error) = parse(grammar, rule.rule, text, pos, cache)
    if error === nothing
      count += 1
      pos = npos
    end
  end

  return (count, pos, nothing)
end

macro grammar(name, definitions)
  parsers = [:+, :*, :?, :|, :-, :^, :list, :integer, :float, :repeated]

  mapped_parsers = Dict{Symbol, Function}()
  for sym in parsers
    mapped_parsers[sym] = eval(sym)
  end

  return parseGrammar(name, definitions, ParserData(mapped_parsers))
end

@grammar custom begin
  start = repeated("abc") + -SPACE + integer
  SPACE = r"[ \t\n]+"
end

data = "abcabcabc 3"
(ast, pos, error) = parse(custom, data)

# value of our new rule should evaluate to the number of times we
# saw a repeat of the given text
@assert parse(Int, ast.children[1]) == parse(Int, ast.children[2].value)


text:
using Replay

instructions = """
using TestImages
using Sixel
sixel_encode(testimage("mandril_color"))
"""

replay(instructions, julia_project=@__DIR__)


text:
function geocode(address::String)
    url = "https://maps.googleapis.com/maps/api/geocode/json?address=$address&key=$(ENV["GOOGLE_MAPS_KEY"])"
    response = HTTP.get(url)
    return(JSON.parse(String(response.body)))
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct Float <: AbstractDelayed
    arg::AbstractDelayed
end

Base.float(arg::AbstractDelayed) = Float(arg)


text:
"""
`BufferedInputStream{T}` provides buffered reading from a source of type `T`.

Any type `T` wrapped in a `BufferedInputStream` must implement:

    BufferedStreams.readbytes!(source::T, buffer::Vector{UInt8}, from::Int, to::Int)

This function should:

* refill `buffer` starting at `from` and not filling past `to`.
* return the number of bytes read.

Failure to read any new data into the buffer is interpreted as eof.
"""
type BufferedInputStream{T} <: IO
    source::T
    buffer::Vector{UInt8}

    # Position of the next byte to be read in buffer;
    # `position ≤ 0` indicates that the stream is closed.
    position::Int

    # Number of bytes available in buffer.
    # I.e. buffer[1:available] is valid data.
    available::Int

    # If positive, preserve and move buffer[anchor:available]
    # when refilling the buffer.
    anchor::Int
end

function BufferedInputStream{T}(source::T, bufsize::Integer=default_buffer_size)
    if bufsize ≤ 0
        throw(ArgumentError("buffer size must be positive"))
    end
    return BufferedInputStream{T}(source, Vector{UInt8}(bufsize), 1, 0, 0)
end

function Base.show{T}(io::IO, stream::BufferedInputStream{T})
    bufsize = length(stream.buffer)
    filled = stream.available - stream.position + 1
    if isopen(stream)
        print(io,
            summary(stream), "(<",
            _datasize(bufsize), " buffer, ",
            round(Int, filled / bufsize * 100), "% filled>)")
    else
        print(io, summary(stream), "(<closed>)")
    end
end

"""
Refill the buffer, optionally moving and retaining part of the data.
"""
function fillbuffer!(stream::BufferedInputStream)
    if eof(stream.source)
        return 0
    end

    if ismarked(stream)
        keepfrom = stream.anchor
        keeplen = stream.available - keepfrom + 1
        # expand the buffer if we are attempting to keep most of it
        if keeplen > div(length(stream.buffer), 2)
            resize!(stream.buffer, length(stream.buffer) * 2)
        end
        stream.anchor = 1
        stream.position -= keepfrom - 1
    else
        keepfrom = stream.position
        keeplen = stream.available - keepfrom + 1
        stream.position = 1
    end

    copy!(stream.buffer, 1, stream.buffer, keepfrom, keeplen)
    nbytes = readbytes!(
        stream.source,
        stream.buffer,
        keeplen + 1,
        endof(stream.buffer))
    stream.available = keeplen + nbytes

    return nbytes
end

@inline function Base.eof(stream::BufferedInputStream)
    if stream.position > stream.available
        return fillbuffer!(stream) == 0
    else
        return false
    end
end

@inline function Base.skip(stream::BufferedInputStream, n_::Integer)
    n0 = n = convert(Int, n_)
    if n < 0
        throw(ArgumentError("n must be non-negative in skip(::BufferedInputStream, n)"))
    end

    while stream.position + n > stream.available + 1
        n -= stream.available - stream.position + 1
        stream.position = stream.available + 1
        if fillbuffer!(stream) < 1
            throw(EOFError())
        end
    end
    stream.position += n
    return n0
end

@inline function checkopen(stream::BufferedInputStream)
    if !isopen(stream)
        error("buffered input stream is already closed")
    end
end

"""
Return the next byte from the input stream without advancing the position.
"""
@inline function peek(stream::BufferedInputStream)
    checkopen(stream)
    if stream.position > stream.available
        if fillbuffer!(stream) < 1
            throw(EOFError())
        end
    end
    @inbounds c = stream.buffer[stream.position]
    return c
end

"""
Fills `buffer` with bytes from `stream`'s buffer without advancing the
position.

Unless the buffer is empty, we do not re-fill it. Therefore the number of bytes
read is limited to the minimum of `nb` and the remaining bytes in the buffer.
"""
function peekbytes!(stream::BufferedInputStream,
                    buffer::AbstractArray{UInt8},
                    nb=length(buffer))
    checkopen(stream)
    if stream.position > stream.available
        if fillbuffer!(stream) < 1
            throw(EOFError())
        end
    end
    nb = min(nb, stream.available - stream.position + 1)
    copy!(buffer, 1, stream.buffer, stream.position, nb)
    return nb
end

@inline function Base.read(stream::BufferedInputStream, ::Type{UInt8})
    checkopen(stream)
    if stream.position > stream.available
        if fillbuffer!(stream) < 1
            throw(EOFError())
        end
    end
    @inbounds c = stream.buffer[stream.position]
    stream.position += 1
    return c
end

# fast multi-byte data readers
for T in [Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
    @eval begin
        @inline function Base.read(stream::BufferedInputStream, ::Type{$(T)})
            checkopen(stream)
            if !ensurebuffered!(stream, $(sizeof(T)))
                throw(EOFError())
            end
            ptr::Ptr{$(T)} = pointer(stream)
            ret = unsafe_load(ptr)
            stream.position += $(sizeof(T))
            return ret
        end
    end
end

# Special purpose readuntil for plain bytes.
function Base.readuntil(stream::BufferedInputStream, delim::UInt8)
    checkopen(stream)
    anchor!(stream)
    while true
        p0 = pointer(stream.buffer, stream.position)
        p1 = ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Cint, Csize_t),
                   p0, delim, stream.available - stream.position + 1)
        if p1 != C_NULL
            stream.position += p1 - p0
            break
        else
            stream.position = stream.available + 1
            nb = fillbuffer!(stream)
            if nb == 0
                chunk = takeanchored!(stream)
                return chunk
            end
        end
    end
    chunk = stream.buffer[upanchor!(stream):stream.position]
    stream.position += 1
    return chunk
end

function readbytes!(stream::BufferedInputStream,
                    buffer::AbstractArray{UInt8},
                    nb=length(buffer))
    return readbytes!(stream, buffer, 1, nb)
end

function readbytes!(stream::BufferedInputStream,
                    buffer::AbstractArray{UInt8},
                    from::Int,
                    to::Int)
    checkopen(stream)
    oldbuflen = buflen = length(buffer)
    nb = to - from + 1
    while !eof(stream) && from <= to
        if stream.position > stream.available && fillbuffer!(stream) < 1
            break
        end

        num_chunk_bytes = min(to - from + 1, stream.available - stream.position + 1)
        if from + num_chunk_bytes > buflen
            buflen = max(buflen + num_chunk_bytes, 2*buflen)
            resize!(buffer, buflen)
        end

        copy!(buffer, from, stream.buffer, stream.position, num_chunk_bytes)
        stream.position += num_chunk_bytes
        from += num_chunk_bytes
    end

    if buflen > oldbuflen
        resize!(buffer, nb - (to - from + 1))
    end

    return nb - (to - from + 1)
end

function Base.ismarked(stream::BufferedInputStream)
    return stream.anchor > 0
end

function Base.mark(stream::BufferedInputStream)
    stream.anchor = stream.position
    return stream.anchor
end

function Base.unmark(stream::BufferedInputStream)
    if !ismarked(stream)
        return false
    end
    stream.anchor = 0
    return true
end

function Base.reset(stream::BufferedInputStream)
    if !ismarked(stream)
        error("buffered stream is not marked")
    end
    anchor = stream.anchor
    stream.position = anchor
    unmark(stream)
    return anchor
end

"""
Return true if the stream is anchored.
"""
function isanchored(stream::BufferedInputStream)
    return stream.anchor > 0
end

"""
Set the buffer's anchor to its current position.
"""
function anchor!(stream::BufferedInputStream)
    stream.anchor = stream.position
end

"""
Remove and return a buffer's anchor.
"""
function upanchor!(stream::BufferedInputStream)
    anchor = stream.anchor
    stream.anchor = 0
    return anchor
end

"""
Copy and return a byte array from the anchor up to, but not including the
current position, also removing the anchor.
"""
function takeanchored!(stream::BufferedInputStream)
    if stream.position - 1 > stream.available
        throw(EOFError())
    end
    chunk = stream.buffer[stream.anchor:stream.position - 1]
    upanchor!(stream)
    return chunk
end

function Base.position(stream::BufferedInputStream)
    return position(stream.source) - stream.available + stream.position - 1
end

function Base.seek{T}(stream::BufferedInputStream{T}, pos::Integer)
    if applicable(seek, stream.source, pos)
        upanchor!(stream)
        source_position = position(stream.source)
        # is the new position within the buffer?
        if source_position - stream.available <= pos <= source_position
            stream.position = 1 + pos - (source_position - stream.available)
        else
            seek(stream.source, pos)
            stream.position = 1
            stream.available = 0
        end
    else
        throw(ArgumentError(
            string("Can't seek in input stream with source of type ", T)))
        # TODO: Allow seeking forwards by just reading and discarding input
    end
end

function Base.isopen(stream::BufferedInputStream)
    return stream.position > 0
end

function Base.close(stream::BufferedInputStream)
    if !isopen(stream)
        return
    end
    if applicable(close, stream.source)
        close(stream.source)
    end
    stream.position = 0
    empty!(stream.buffer)
    return
end

function Base.pointer(stream::BufferedInputStream, index::Integer=1)
    return pointer(stream.buffer, stream.position + index - 1)
end

@inline function available_bytes(stream::BufferedInputStream)
    return stream.available - stream.position + 1
end

@inline function ensurebuffered!(stream::BufferedInputStream, nb::Integer)
    if available_bytes(stream) < nb
        fillbuffer!(stream)
        if available_bytes(stream) < nb
            return false
        end
    end
    return true
end


text:
using ITensors

let
  Ny = 6
  Nx = 12

  N = Nx * Ny

  sites = siteinds("S=1/2", N; conserve_qns=true)

  lattice = square_lattice(Nx, Ny; yperiodic=false)

  ampo = AutoMPO()
  for b in lattice
    ampo .+= 0.5, "S+", b.s1, "S-", b.s2
    ampo .+= 0.5, "S-", b.s1, "S+", b.s2
    ampo .+= "Sz", b.s1, "Sz", b.s2
  end
  H = MPO(ampo, sites)

  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  # Initialize wavefunction to a random MPS
  # of bond-dimension 10 with same quantum
  # numbers as `state`
  psi0 = randomMPS(sites, state, 20)

  sweeps = Sweeps(10)
  setmaxdim!(sweeps, 20, 60, 100, 100, 200, 400, 800)
  setcutoff!(sweeps, 1E-8)
  @show sweeps

  energy, psi = dmrg(H, psi0, sweeps)

  return nothing
end


text:
using Distributed

using Test

@everywhere begin
    using MLJBase
    using MLJTuning
    using ..Models
    import ComputationalResources: CPU1, CPUProcesses, CPUThreads
end

using Random
Random.seed!(1234*myid())
using .TestUtilities

N = 30
x1 = rand(N);
x2 = rand(N);
x3 = rand(N);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.4*rand(N);

m(K) = KNNRegressor(K=K)
r = [m(K) for K in 13:-1:2]

# TODO: replace the above with the line below and post an issue on
# the failure (a bug in Distributed, I reckon):
# r = (m(K) for K in 13:-1:2)

@testset "constructor" begin
    @test_throws(MLJTuning.ERR_SPECIFY_RANGE,
                 TunedModel(model=first(r), tuning=Grid(), measure=rms))
    @test_throws(MLJTuning.ERR_SPECIFY_RANGE,
                 TunedModel(model=first(r), measure=rms))
    @test_throws(MLJTuning.ERR_BOTH_DISALLOWED,
                 TunedModel(model=first(r),
                            models=r, tuning=Explicit(), measure=rms))
    tm = TunedModel(models=r, tuning=Explicit(), measure=rms)
    @test tm.tuning isa Explicit && tm.range ==r && tm.model == first(r)
    @test input_scitype(tm) == Unknown
    @test TunedModel(models=r, measure=rms) == tm
    @test_logs (:info, r"No measure") @test TunedModel(models=r) == tm

    @test_throws(MLJTuning.ERR_SPECIFY_MODEL,
                 TunedModel(range=r, measure=rms))
    @test_throws(MLJTuning.ERR_MODEL_TYPE,
                 TunedModel(model=42, tuning=Grid(),
                            range=r, measure=rms))
    @test_logs (:info, MLJTuning.INFO_MODEL_IGNORED) tm =
        TunedModel(model=42, tuning=Explicit(), range=r, measure=rms)
    @test_logs (:info, r"No measure") tm =
        TunedModel(model=first(r), range=r)
    @test_throws(MLJTuning.ERR_SPECIFY_RANGE_OR_MODELS,
                TunedModel(tuning=Explicit(), measure=rms))
    @test_throws(MLJTuning.ERR_NEED_EXPLICIT,
                 TunedModel(models=r, tuning=Grid()))
    tm = TunedModel(model=first(r), range=r, measure=rms)
    @test tm.tuning isa Grid
    @test input_scitype(tm) == Table(Continuous)
end

results = [(evaluate(model, X, y,
                     resampling=CV(nfolds=2),
                     measure=rms,
                     verbosity=0,)).measurement[1] for model in r]

@testset "measure compatibility check" begin
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=cross_entropy)
    @test_logs((:error, r"Problem"),
               (:info, r""),
               (:info, r""),
               @test_throws ArgumentError fit(tm, 0, X, y))
end

@testset_accelerated "basic fit (CPU1)" accel begin
    printstyled("\n Testing progressmeter basic fit with $(accel) and CPU1 resampling \n", color=:bold)
    best_index = argmin(results)
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=[rms, l1], acceleration=accel)
    verbosity = accel isa CPU1 ? 2 : 1
    fitresult, meta_state, _report = fit(tm, verbosity, X, y);
    history, _, state = meta_state;
    results2 = map(event -> event.measurement[1], history)
    @test results2 ≈ results
    @test fitresult.model == collect(r)[best_index]
    @test _report.best_model == collect(r)[best_index]
    @test _report.history[5] == MLJTuning.delete(history[5], :metadata)

    # training_losses:
    losses = training_losses(tm, _report)
    @test all(eachindex(losses)) do i
        minimum(results[1:i]) == losses[i]
    end
    @test MLJBase.iteration_parameter(tm) == :n

end

@static if VERSION >= v"1.3.0-DEV.573"
@testset_accelerated "Basic fit (CPUThreads)" accel begin
    printstyled("\n Testing progressmeter basic fit with $(accel) and CPUThreads resampling \n", color=:bold)
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=[rms, l1], acceleration= CPUThreads(),
                    acceleration_resampling=accel)
    fitresult, meta_state, report = fit(tm, 1, X, y);
    history, _, state = meta_state;
    results3 = map(event -> event.measurement[1], history)
    @test results3 ≈ results
end
end
@testset_accelerated "Basic fit (CPUProcesses)" accel begin
    printstyled("\n Testing progressmeter basic fit with $(accel) and CPUProcesses resampling \n", color=:bold)
    best_index = argmin(results)
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, resampling=CV(nfolds=2),
                    measures=[rms, l1], acceleration=CPUProcesses(),
                    acceleration_resampling=accel)
    fitresult, meta_state, report = fit(tm, 1, X, y);
    history, _, state = meta_state;
    results4 = map(event -> event.measurement[1], history)
    @test results4 ≈ results
end

@testset_accelerated("under/over supply of models", accel, begin
    tm = TunedModel(model=first(r), tuning=Explicit(),
                    range=r, measures=[rms, l1],
                    acceleration=accel,
                    resampling=CV(nfolds=2),
                    n=4)
    mach = machine(tm, X, y)
    fit!(mach, verbosity=0)
    history = MLJBase.report(mach).history
    @test map(event -> event.measurement[1], history) ≈ results[1:4]

    tm.n += 2
    @test_logs((:info, r"^Updating"),
               (:info, r"^Attempting to add 2.*to 6"),
               fit!(mach, verbosity=1))
    history = MLJBase.report(mach).history
    @test map(event -> event.measurement[1], history) ≈ results[1:6]

    tm.n=100
    @test_logs (:info, r"Only 12") fit!(mach, verbosity=0)
    history = MLJBase.report(mach).history
    @test map(event -> event.measurement[1], history) ≈ results
end)

@everywhere begin

    # variation of the Explicit strategy that annotates the models
    # with metadata
    mutable struct MockExplicit <: MLJTuning.TuningStrategy end

    annotate(model) = (model, params(model)[1])

    _length(x) = length(x)
    _length(::Nothing) = 0
    MLJTuning.models(tuning::MockExplicit,
                     model,
                     history,
                     state,
                     n_remaining,
                     verbosity) =
                         annotate.(state)[_length(history) + 1:end], state

    function default_n(tuning::Explicit, range)
        try
            length(range)
        catch MethodError
            DEFAULT_N
        end

    end

end

@testset_accelerated("passing of model metadata", accel,
                  begin
                     tm = TunedModel(model=first(r), tuning=MockExplicit(),
                                     range=r, resampling=CV(nfolds=2),
                                     measures=[rms, l1], acceleration=accel)
                     fitresult, meta_state, report = fit(tm, 0, X, y);
                     history, _, state = meta_state;
                     @test all(history) do event
                         event.metadata == event.model.K
                     end
                     end)

@testset "data caching" begin
    X = (x1= ones(5),);
    y = coerce(collect("abcaa"), Multiclass);

    m(b) = ConstantClassifier(testing=true, bogus=b)
    r = [m(b) for b in 1:5]

    tuned_model = TunedModel(model=first(r), tuning=Explicit(),
                             range=r, resampling=Holdout(fraction_train=0.8),
                             measure=log_loss, cache=true)

    # There is reformatting and resampling of X and y for training of
    # first model, and then resampling on the prediction side only
    # thereafter. Then, finally, for training on all supplied data, there
    # is reformatting and resampling again:
    @test_logs(
        (:info, "reformatting X, y"), # fit 1
        (:info, "resampling X, y"),   # fit 1
        (:info, "resampling X"),      # predict 1
        (:info, "resampling X"),      # predict 2
        (:info, "resampling X"),      # predict 3
        (:info, "resampling X"),      # predict 4
        (:info, "resampling X"),      # predict 5
        (:info, "reformatting X, y"), # fit on all data
        (:info, "resampling X, y"),   # fit one all data
        fit(tuned_model, 0, X, y)
    );

    # Otherwise, resampling and reformatting happen for every model
    # trained, and we get a reformat on the prediction side every
    # time:
    tuned_model.cache = false
    @test_logs(
        (:info, "reformatting X, y"), # fit 1
        (:info, "resampling X, y"),   # fit 1
        (:info, "reformatting X"),    # predict 1
        (:info, "reformatting X, y"), # fit 2
        (:info, "resampling X, y"),   # fit 2
        (:info, "reformatting X"),    # predict 2
        (:info, "reformatting X, y"), # fit 3
        (:info, "resampling X, y"),   # fit 3
        (:info, "reformatting X"),    # predict 3
        (:info, "reformatting X, y"), # fit 4
        (:info, "resampling X, y"),   # fit 4
        (:info, "reformatting X"),    # predict 4
        (:info, "reformatting X, y"), # fit 5
        (:info, "resampling X, y"),   # fit 5
        (:info, "reformatting X"),    # predict 5
        (:info, "reformatting X, y"), # fit on all data
        (:info, "resampling X, y"),   # fit on all data
        fit(tuned_model, 0, X, y)
    );
end

@testset "issue #128" begin
    X, y = make_regression(10, 2)
    dtc = DecisionTreeRegressor()
    r   = range(dtc, :max_depth, lower=1, upper=50);

    tmodel = TunedModel(model=dtc, ranges=[r, ],
                        tuning=Grid(resolution=50),
                        measure=mae,
                        n=48);
    mach = machine(tmodel, X, y)
    fit!(mach, verbosity=0);

    @test length(report(mach).history) == 48

    tmodel.n = 49
    fit!(mach, verbosity=0);

    @test length(report(mach).history) == 49
end

@testset_accelerated "Resampling reproducibility" accel begin
    X, y = make_regression(100, 2)
    dcr = DeterministicConstantRegressor()

    # Hold out reproducibility
    homodel = TunedModel(tuning=Explicit(),
                         models=fill(dcr, 10),
                         resampling=Holdout(rng=StableRNG(1234)),
                         acceleration_resampling=accel,
                         measure=mae)
    homach = machine(homodel, X, y)
    fit!(homach, verbosity=0);
    horep = report(homach)
    measurements = getproperty.(horep.history, :measurement)
    @test all(==(measurements[1]), measurements)

    # Cross-validation reproducibility
    cvmodel = TunedModel(tuning=Explicit(),
                         models=fill(dcr, 10),
                         resampling=CV(nfolds=5, rng=StableRNG(1234)),
                         acceleration_resampling=accel,
                         measure=mae)
    cvmach = machine(cvmodel, X, y)
    fit!(cvmach, verbosity=0);
    cvrep = report(cvmach)
    per_folds = getproperty.(cvrep.history, :per_fold)
    @test all(==(per_folds[1]), per_folds)
end

@testset "deterministic metrics for probabilistic models" begin

    # https://github.com/JuliaAI/MLJBase.jl/pull/599 allows mix of
    # deterministic and probabilistic metrics:
    X, y = MLJBase.make_blobs()
    model = DecisionTreeClassifier()
    range = MLJBase.range(model, :max_depth, values=[1,2])
    tmodel = TunedModel(model=model,
                        range=range,
                        measures=[MisclassificationRate(),
                                  LogLoss()])
    mach = machine(tmodel, X, y)
    @test_logs fit!(mach, verbosity=0)

end

true


text:
function maketype(name,param_dict,origex,funcs,syms,fex;
                  pex=:(), vector_ex = :(), vector_ex_return = :(),
                  tgradex = :(),
                  Jex=:(),
                  expJex=:(),
                  invJex=:(),invWex=:(),invWex_t=:(),
                  Hex = :(),
                  invHex = :(),
                  params = Symbol[],
                  pfuncs=Vector{Expr}(undef,0),
                  d_pfuncs = Vector{Expr}(undef,0),
                  symjac = Matrix{Expr}(undef,0,0),
                  symtgrad = Vector{Expr}(undef,0),
                  param_Jex=:(),
                  f_expr=:(),
                  tgrad_expr=:(),
                  jac_expr=:(),
                  invjac_expr=:(),
                  invW_expr=:(),
                  invW_t_expr=:(),
                  param_jac_expr=:())

    typeex = :(mutable struct $name{F,J,T,W,Wt,PJ,TT1,TT2} <: DiffEqBase.AbstractParameterizedFunction{true}
        f::F
        analytic::Nothing
        jac::J
        tgrad::T
        invW::W
        invW_t::Wt
        paramjac::PJ
        mass_matrix::LinearAlgebra.UniformScaling{Bool}
        jac_prototype::Nothing
        origex::Expr
        funcs::Vector{Expr}
        pfuncs::Vector{Expr}
        d_pfuncs::Vector{Expr}
        syms::Vector{Symbol}
        symjac::Matrix{TT1} # https://github.com/symengine/SymEngine.jl/issues/122
        symtgrad::Vector{TT2}
        tgradex::Expr
        Jex::Expr
        expJex::Expr
        param_Jex::Expr
        invJex::Expr
        invWex::Expr
        invWex_t::Expr
        Hex::Expr
        invHex::Expr
        fex::Expr
        pex::Expr
        vector_ex::Expr
        vector_ex_return::Expr
        params::Vector{Symbol}
        colorvec::Nothing
    end)

    # Make the default constructor
    new_ex = Meta.quot(origex)
    tgradex_ex = Meta.quot(tgradex)
    Jex_ex = Meta.quot(Jex)
    expJex_ex = Meta.quot(expJex)
    invJex_ex = Meta.quot(invJex)
    invWex_ex = Meta.quot(invWex)
    invWex_t_ex = Meta.quot(invWex_t)
    Hex_ex = Meta.quot(Hex)
    invHex_ex = Meta.quot(invHex)
    fex_ex = Meta.quot(fex)
    pex_ex = Meta.quot(pex)
    vector_ex_ex = Meta.quot(vector_ex)
    vector_ex_return_ex = Meta.quot(vector_ex_return)
    param_Jex_ex = Meta.quot(param_Jex)

    constructorex = :($(name)() =
                  $(name)($f_expr,nothing,
                  $jac_expr,$tgrad_expr,$invW_expr,$invW_t_expr,$param_jac_expr,
                  $(LinearAlgebra.I),nothing,$new_ex,$funcs,$pfuncs,$d_pfuncs,
                  $syms,$symjac,$symtgrad,
                  $tgradex_ex,$Jex_ex,$expJex_ex,$param_Jex_ex,
                  $invJex_ex,$invWex_ex,$invWex_t_ex,
                  $Hex_ex,$invHex_ex,$fex_ex,$pex_ex,$vector_ex_ex,
                  $vector_ex_return_ex,$params,nothing)) |> esc

    callex = :(((f::$name))(args...) = f.f(args...)) |> esc
    callex2 = :(((f::$name))(u,p,t::Number) = (du=similar(u);f.f(du,u,p,t);du)) |> esc

    # Make the type instance using the default constructor
    typeex,constructorex,callex,callex2
end

DiffEqBase.has_Wfact(::DiffEqBase.AbstractParameterizedFunction) = false
DiffEqBase.has_Wfact_t(::DiffEqBase.AbstractParameterizedFunction) = false


text:
using Plots
# plots the contours of the OLS objective function,
# with the true and estimated coefficients
# generate data
n = 20
x = [ones(n) randn(n)]
beta = randn(2)
y = x*beta+randn(n)
# ols estimate and objective function
betahat = x\y
s = (a,b)-> (1/n)*sum((y - x*[a,b]).^2)
# plot the contours and the points
closeall()
b1 = range(-4.0,stop=4.0,length=100)
b2 = range(-4.0,stop=4.0,length=100)
p = contour(b1,b2,(b1,b2)->s(b1,b2),fill=true, c=:viridis)
scatter!([beta[1]],[beta[2]], markersize=10, label="true")
scatter!([betahat[1]],[betahat[2]], markersize=10, label="estimated")

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
module FileSessionAdapter

using Nullables
using Genie.Sessions, Genie, Genie.Loggers, Genie.Configuration
using Serialization

"""
    write(session::Sessions.Session) :: Sessions.Session

Persists the `Session` object to the file system, using the configured sessions folder and returns it.
"""
function write(session::Sessions.Session) :: Sessions.Session
  if ! isdir(joinpath(Genie.SESSIONS_PATH))
    log("The configured sessions folder $(Genie.SESSIONS_PATH) does not exist -- switching to using a temporary folder.")
    Core.eval(Genie, :(Genie.SESSIONS_PATH = mktempdir()))
  end

  try
    open(joinpath(Genie.SESSIONS_PATH, session.id), "w") do (io)
      serialize(io, session)
    end
  catch ex
    log("Error when serializing session in $(@__FILE__):$(@__LINE__)", :err)
    log(string(ex), :err)
    log("$(@__FILE__):$(@__LINE__)", :err)

    Genie.Configuration.isprod() && rethrow(ex)
  end

  session
end

"""
    read(session_id::Union{String,Symbol}) :: Nullable{Sessions.Session}
    read(session::Sessions.Session) :: Nullable{Sessions.Session}

Attempts to read from file the session object serialized as `session_id`.
"""
function read(session_id::Union{String,Symbol}) :: Nullable{Sessions.Session}
  try
    isfile(joinpath(Genie.SESSIONS_PATH, session_id)) || return Nullable{Sessions.Session}(write(Session(session_id)))
  catch ex
    log("Can't check session file", :err)
    log(string(ex), :err)
    log("$(@__FILE__):$(@__LINE__)", :err)

    Nullable{Sessions.Session}(write(Session(session_id)))
  end

  try
    session = open(joinpath(Genie.SESSIONS_PATH, session_id), "r") do (io)
      deserialize(io)
    end

    Nullable{Sessions.Session}(session)
  catch ex
    log("Can't read session", :err)
    log(string(ex), :err)
    log("$(@__FILE__):$(@__LINE__)", :err)

    Nullable{Sessions.Session}(write(Session(session_id)))
  end
end
function read(session::Sessions.Session) :: Nullable{Sessions.Session}
  read(session.id)
end

end


text:
DATAPATH = joinpath(@__DIR__, "data")

"""
    load_iris() -> X, y, names

Loads the first 150 observations from the
Iris flower data set introduced by Ronald Fisher (1936).
The 4 by 150 matrix `X` contains the numeric measurements,
in which each individual column denotes an observation.
The vector `y` contains the class labels as strings.
The vector `names` contains the names of the features (i.e. rows of `X`)

[1] Fisher, Ronald A. "The use of multiple measurements in taxonomic problems." Annals of eugenics 7.2 (1936): 179-188.
"""
function load_iris()
    path = joinpath(DATAPATH, "iris.csv")
    raw_csv = readdlm(path, ',')
    X = convert(Matrix{Float64}, raw_csv[:, 1:4]')
    y = convert(Vector{String}, raw_csv[:, 5])
    vars = ["Sepal length", "Sepal width", "Petal length", "Petal width"]
    X, y, vars
end


text:
using Documenter, MQTT

makedocs()

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
if !isdefined(Reactive, :number)
    number() = rand(0:100)
end

@testset "Call counting" begin
    a = Signal(0)
    b = Signal(0)

    c = map(+, a, b)
    d = merge(a, b)
    e = map(+, a, map(x->2x, a)) # Both depend on a
    f = map(+, a, b, c, e)

    count = s -> foldp((x, y) -> x+1, 0, s)

    ca = count(a)
    cb = count(b)
    cc = count(c)
    cd = count(d)
    ce = count(e)
    cf = count(f)

    for i in 1:100
        push!(a, number())
        step()
        push!(b, number())
        step()

        @test (ca.value) == (i)
        @test (cb.value) == (i)
        @test (cc.value) == (2i)
        @test (cd.value) == (2i)
        @test (ce.value) == (i)
        @test (cf.value) == (2i)
    end

end


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

module ClassifyModel

using CSV
using DataFrames
using DataStructures

function loadlabel(fname::AbstractString)
    df = CSV.read(fname, delim = ",", header = 0)
    list = filter(s -> !startswith(coalesce(s, "#"), "#"), df[1])
    OrderedDict(enumerate(list))
end

Labels = collect(loadlabel("data/class_list.txt"))

mutable struct ClassifyData
    queries::Vector{String}
    labels::Vector{String}
    size::Integer
    offset::Integer
    ClassifyData(queryset::Vector{String}) = new(queryset, fill("", size(queryset, 1)), size(queryset, 1), 0)
    function ClassifyData(queryset::Vector{String}, labelset::Vector{String})
        @assert size(queryset, 1) >= size(labelset, 1)
        new(queryset, vcat(labelset, fill("", size(queryset, 1) - size(labelset, 1))),
            size(queryset, 1), size(labelset, 1))
    end
end

function getquery(data::ClassifyData)
    (data.offset < data.size) ? data.queries[data.offset + 1] : nothing
end

function update!(data::ClassifyData, label::AbstractString)
    data.offset += 1
    data.labels[data.offset] = label
end

function loadquery(fquery::AbstractString, fcache::AbstractString = "")
    dfq = CSV.read(fquery, delim = ",", header = 0)
    queryset = collect(skipmissing(dfq[1]))
    if fcache != ""
        dfc = CSV.read(fcache, delim = ",", header = 0)
        index = minimum(filter(k -> ismissing(dfc[k, 2]) || dfc[k, 2] == "", 1:size(df, 1)))
        if index > 1
            @assert dfc[1:(index - 1), 1] == queryset[1:(index - 1)]
            labelset = collect(skipmissing(dfc[1:(index - 1), 2]))
            return ClassifyData(queryset, labelset)
        end
    end
    return ClassifyData(queryset)
end

function dumpdata(data::ClassifyData)
    df = DataFrame(query = data.queries, label = data.labels)
    CSV.write("data/classify_out.csv", df, delim = ",", writeheader = false)
end

DataSet = loadquery("data/demo_query.csv")

end


text:
include("../src/opti_module.jl")
include("../opti_test_functions/testfunctions.jl")

module gtest
using BenchmarkTools
using ..opti
using ..optitestfun

function test_steepest_gradient_descent()
    # fcn(x) = x[1]^2 + x[2]^2
    fcn(x) = optitestfun.rosenbrock(x, 2)
    # x0 = [1.0, 2]
    x0 = [-1.8, -2]
    stat = opti.steepest_gradient_descent(fcn, x0, tol_options = opti.ToleranceOptions(max_iter = 1000))

    # stat = opti.steepest_gradient_descent_original(fcn, x0, tol_options = opti.ToleranceOptions(max_iter = 1000))

    show(stat)

end

function test_steepest_gradient_descent_v2()

    fcn(x) = optitestfun.rosenbrock(x, 2)

    x0 = [-1.8, -2]
    (stat, opti_search_path) = opti.steepest_gradient_descent(fcn, x0, tol_options = opti.ToleranceOptions(max_iter = 1000), log_path = true)

    show(stat)

    optitestfun.plot_rosenbrock([-2, 5], [-2, 5])
    PyPlot.plot( opti_search_path.x[1, :], opti_search_path.x[2, :], marker = :., markersize = 15 )
end

function benchmark_steepestgrad()

    fcn(x) = optitestfun.rosenbrock(x, 2)
    x0 = [-1.8, -2]
    b = @benchmark opti.steepest_gradient_descent($fcn, $x0, tol_options = opti.ToleranceOptions(max_iter = 1000))

    println("=== steepest_gradient_descent ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_steepestgrad_orig()

    fcn(x) = optitestfun.rosenbrock(x, 2)
    x0 = [-1.8, -2]
    b = @benchmark opti.steepest_gradient_descent_original($fcn, $x0, tol_options = opti.ToleranceOptions(max_iter = 1000))

    println("=== steepest_gradient_descent_original ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_cg_pr()

    fcn(x) = optitestfun.rosenbrock(x, 2)
    x0 = [-1.8, -2]
    b = @benchmark opti.cg_pr($fcn, $x0)

    println("=== cg_pr ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function test_cg_pr()
    # fcn(x) = x[1]^2 + x[2]^2
    fcn(x) = optitestfun.rosenbrock(x, 4)
    x0 = [-1.8, -2]

    stat = opti.cg_pr(fcn, x0; tol_options = opti.ToleranceOptions(f_abs_tol = 1e-10))

    show(stat)

end

function test_cg_pr_quad()
    fcn(x) = optitestfun.nd_quadratic(x)
    n_dim = 100
    x0 = 10 * rand(n_dim)

    stat = opti.cg_pr(fcn, x0; tol_options = opti.ToleranceOptions(f_abs_tol = 1e-10))

    show(stat)
end

function test_cg_pr_quadwithl1()

    n_dim = 10
    x_sol = ones(n_dim)
    fcn(x) = optitestfun.quad_with_norm_l1(x, x_sol, 0.8)

    x0 = 10 * rand(n_dim)

    stat = opti.cg_pr(fcn, x0; tol_options = opti.ToleranceOptions(f_abs_tol = 1e-10))

    show(stat)

end

function test_cg_pr_v2()
    # fcn(x) = x[1]^2 + x[2]^2
    fcn(x) = optitestfun.rosenbrock(x, 2)
    x0 = [-1.8, -2]

    (stat, opti_search_path) = opti.cg_pr(fcn, x0, log_path = true)

    show(stat)

    optitestfun.plot_rosenbrock([-2, 5], [-2, 5])
    PyPlot.plot( opti_search_path.x[1, :], opti_search_path.x[2, :], marker = :., markersize = 15 )
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license
export evaluate, selection_string

"""
Create a ``Selection`` from a selection string
"""
function Selection(selection::AbstractString)
    handle = lib.chfl_selection(pointer(selection))
    return Selection(handle)
end

"""
Get the selection string used to create a given ``selection``.
"""
function selection_string(selection::Selection)
    return _call_with_growing_buffer(
        (buffer, size) -> _check(lib.chfl_selection_string(selection.handle, buffer, size))
    )
end

"""
Get the size of a ``selection``, *i.e.* the number of atoms we are selecting
together.
"""
function Base.size(selection::Selection)
    result = Ref{UInt64}(0)
    _check(
        lib.chfl_selection_size(selection.handle, result)
    )
    return result[]
end

"""
Evaluate a ``selection`` on a given ``frame``. This function return a list of
indexes or tuples of indexes of atoms in the frame matching the selection.
"""
function evaluate(selection::Selection, frame::Frame)
    matching = Ref{UInt64}(0)
    _check(
        lib.chfl_selection_evaluate(selection.handle, frame.handle, matching)
    )
    matching = matching[]
    matches = Array{lib.chfl_match}(matching)
    _check(
        lib.chfl_selection_matches(selection.handle, pointer(matches), matching)
    )
    result = []
    selection_size = size(selection)
    for match in matches
        assert(match.size == selection_size)
        if selection_size == 1
            push!(result, match.atoms_1)
        elseif selection_size == 2
            push!(result, (match.atoms_1, match.atoms_2))
        elseif selection_size == 3
            push!(result, (match.atoms_1, match.atoms_2, match.atoms_3))
        elseif selection_size == 4
            push!(result, (match.atoms_1, match.atoms_2, match.atoms_3, match.atoms_4))
        end
    end
    return result
end

"""
Free the allocated memory for the ``Selection`` object.
"""
function _free(selection::Selection)
    lib.chfl_selection_free(selection.handle)
end

"""
Make a deep copy of a ``selection``.
"""
function Base.deepcopy(selection::Selection)
    handle = lib.chfl_selection_copy(selection.handle)
    return Selection(handle)
end


text:
using Plots
function export_plot(
        token_stream,
        code;
        path::S,
        filename::S = "./scope_per_flag.png",
    ) where {S}
    ei = eachindex(code)
    plot()
    @inbounds for (flag, scope) in token_stream
        plot!([i for i in ei], scope, label=replace(flag, "\n" => ""))
    end
    scope_sum = [0 for i in ei]
    @inbounds for (flag, scope) in token_stream
        scope_sum .+= scope
    end
    plot!([i for i in ei], scope_sum, label="total scope", legend=:bottomleft)
    savefig(joinpath(path,filename))
end


text:
srcs=[[162,237,1] [220,979,1] [1237,861,1] [1254,301,1.]];

tgts=[[193,209,1] [224,854,1] [1308,867,1] [1359,227,1.]];

src_img = load("/home/somvt/work/DiffImages.jl/examples/source.jpeg")'

tgt_img = load("/home/somvt/work/DiffImages.jl/examples/target.jpeg")'

src2 = imresize(src_img,ratio=1/8)
tgt2 = imresize(tgt_img,ratio=1/8)
src3 = map(x->Gray(x).val, src2)
tgt3 = map(x->Gray(x).val,tgt2)
function warp_me2(img::AbstractArray{T}, tform) where T
    img = ImageTransformations.box_extrapolation(img, RGB(.0,.0,.0))
    # inds = DiffImages.custom_autorange(img, inv(tform))
    inds = map(x->SVector{length(x),Float64}(x.I),CartesianIndices(img))
    # tgts = H * srcs
    # img = map(x -> ImageTransformations._getindex(img, tform(x)), inds)
    img = map(x->img(inv(tform)(x)...),inds)
    return img
end

h=homography2d(srcs,tgts);
h=DiffImages.homography{Float64}(h);
function L1_float(y,ŷ)
    loss = 0
    for i in CartesianIndices(y)
        if checkbounds(Bool,ŷ,i)
            s = ŷ[i]-y[i]
            if s > 0
                loss += s
            else
                loss += -s
            end
        else
            if y[i] > 0
                loss += y[i]
            else
                loss += -y[i]
            end
        end
    end
    loss
end

function L1_color(y,ŷ)
    # ŷ => warped image, y => target image
    loss = 0
    for i in CartesianIndices(y)
        if checkbounds(Bool,ŷ,i)
            r = ŷ[i]-y[i]
            loss += abs(r.r+r.g+r.b)
        else
            c = y[i]
            loss += abs(c.r+c.g+c.b)
        end
    end
    loss
end
# L(y,ȳ) = sum(abs.(ȳ-begin
#                    y = map(x->x.val,y)
#                    y
#                    end)) # => y in Gray
# for i in 1:num_iters
#     ∇ = Zygote.gradient(src, h) do img, tform
#         c = sum(warp_me2(img, tform))
#         gray(c)
#         end[2].H
#     h = map((y,ȳ)->map(L(y,ȳ)),m,∇)
#     h = homography(SMatrix{3,3,Float64,9}(h))
# end

max_stars_repo_path: kornia vs julia/mytests.jl

max_stars_repo_name: SomTambe/reproduction

max_issues_repo_path: kornia vs julia/mytests.jl

max_issues_repo_name: SomTambe/reproduction

max_forks_repo_path: kornia vs julia/mytests.jl

max_forks_repo_name: SomTambe/reproduction


text:
cd(dirname(@__FILE__)) do
    run(`make`)
end


text:
# Tests of Adjacency List types

using Graphs
using Test

gd = simple_adjlist(3)
gu = simple_adjlist(3; is_directed=false)

# concept test

@test implements_vertex_list(gd)    == true
@test implements_edge_list(gd)      == false
@test implements_vertex_map(gd)     == true
@test implements_edge_map(gd)       == false

@test implements_adjacency_list(gd) == true
@test implements_incidence_list(gd) == false
@test implements_bidirectional_adjacency_list(gd) == false
@test implements_bidirectional_incidence_list(gd) == false
@test implements_adjacency_matrix(gd) == false

# graph without edges

@test vertex_type(gd) == Int
@test edge_type(gd) == Edge{Int}

@test num_vertices(gd) == 3
@test num_edges(gd) == 0
@test vertices(gd) == 1:3
@test is_directed(gd) == true

@test num_vertices(gu) == 3
@test num_edges(gu) == 0
@test vertices(gu) == 1:3
@test is_directed(gu) == false

for i = 1 : 3
    @test out_degree(i, gd) == 0
    @test isempty(out_neighbors(i, gd))
    @test out_degree(i, gu) == 0
    @test isempty(out_neighbors(i, gu))
end

# graph with edges

add_edge!(gd, 1, 2)
add_edge!(gd, 1, 3)
add_edge!(gd, 2, 3)

@test num_edges(gd) == 3

nbs = Any[[2, 3], [3], Int[]]
for i = 1 : 3
    @test out_degree(i, gd) == length(nbs[i])
    @test out_neighbors(i, gd) == nbs[i]
end

add_edge!(gu, 1, 2)
add_edge!(gu, 2, 3)

@test num_edges(gu) == 2

nbs = Any[[2], [1, 3], [2]]
for i = 1 : 3
    @test out_degree(i, gu) == length(nbs[i])
    @test out_neighbors(i, gu) == nbs[i]
end

# adjacency list with key vertices

for g in [adjlist(KeyVertex{String}), adjlist(String)]

    vs = [  add_vertex!(g, "a"),
            add_vertex!(g, "b"),
            add_vertex!(g, "c") ]

    @test num_vertices(g) == 3

    for i = 1 : 3
        v = vs[i]
        @test vertices(g)[i] == v
        @test !isa(v, KeyVertex) || v.index == i
        @test out_degree(v, g) == 0
    end

    add_edge!(g, vs[1], vs[2])
    add_edge!(g, vs[1], vs[3])
    add_edge!(g, vs[2], vs[3])

    @test num_edges(g) == 3

    @test out_degree(vs[1], g) == 2
    @test out_degree(vs[2], g) == 1
    @test out_degree(vs[3], g) == 0

    @test out_neighbors(vs[1], g) == [vs[2], vs[3]]
    @test out_neighbors(vs[2], g) == [vs[3]]
    @test isempty(out_neighbors(vs[3], g))

end

# Construct an AdjacencyList via a Vector of vertices

g = adjlist(collect(1:3))
add_edge!(g, 2, 1)
add_edge!(g, 2, 3)

@test num_vertices(g) == 3
@test num_edges(g) == 2

@test out_degree(g.vertices[1], g) == 0
@test out_degree(g.vertices[2], g) == 2
@test out_degree(g.vertices[3], g) == 0

@test isempty(out_neighbors(g.vertices[1], g))
@test out_neighbors(g.vertices[2], g) == [1, 3]
@test isempty(out_neighbors(g.vertices[3], g))

# # construct via adjacency matrix
# A = [true true true; false false true; false false true]
# g = simple_adjlist(A)
# @test adjacency_matrix(g) == A


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgRuntimeRawExtension <: SwaggerModel

    function IoK8sApimachineryPkgRuntimeRawExtension(;)
        o = new()
        o
    end
end # type IoK8sApimachineryPkgRuntimeRawExtension

const _property_map_IoK8sApimachineryPkgRuntimeRawExtension = Dict{Symbol,Symbol}()
const _property_types_IoK8sApimachineryPkgRuntimeRawExtension = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }) = collect(keys(_property_map_IoK8sApimachineryPkgRuntimeRawExtension))
Swagger.property_type(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgRuntimeRawExtension[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgRuntimeRawExtension[property_name]

function check_required(o::IoK8sApimachineryPkgRuntimeRawExtension)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, name::Symbol, val)
end