text:
using Plots, Latexify, LaTeXStrings

e = [2.15, 3.11, 4.81, 7.06, 10, 13.1, 17.1, 22.1, 25.1, 28.4, 31.6, 34.7, 36.8, 41.7, 46.6, 51.7, 56]

b₀ = [21280, 12260, 9312, 6937, 6219, 5105, 5363, 3960, 3063, 2967, 3113, 2449, 2237, 1979, 2266, 2266, 2266]

b₁ = [1231, 973.3, 776.8, 697.9, 589.8, 572.8, 415.7, 362.1, 366, 333.3, 281.1, 302, 314, 325.3, 302, 302, 302]

scatter(e,b₀,legend=true,
 xlabel="Thickness", ylabel = "Unit pressure drop", labels = "b₀", title = "NC 2733 foam", xlim=(0,60),ylim=(0,25000))
plot!(e, b₁, label = "b₁", seriestype = :scatter)

A = [e[i]^j for i in 1:17, j in -1:0]

ê₀ = A \ b₀

b̂₀(e) = ê₀' * [e^-1,1]

ê₁ = A \ b₁

b̂₁(e) = ê₁' * [e^-1,1]

plot!(b̂₀, label = "b₀(e)")
plot!(b̂₁, label = "b₁(e)")

E = [5.12, 10, 13.2, 15.6, 20.4, 21.7, 25.3, 30.7, 37.5, 41.8, 44.9, 47.2, 50.8, 55, 59.9, 62.1]

b₂ = [5665, 2689, 2144, 2236, 1592, 1872, 1151, 1196, 1334, 947, 1076, 1060, 1093, 888.7, 817.1, 949.2]

b₃ = [256.2, 271.7, 235.9, 208.6, 171.3, 139.5, 155.8, 127.5, 93.9, 104.2, 91.3, 86.1, 79.9, 87.6, 86.7, 82.4]

scatter(E,b₂,legend=true,
 xlabel="Thickness", ylabel = "Unit pressure drop", label = "b₂", xlim=(0,70), ylim=(0,6000))
plot!(E, b₃, label = "b₃",seriestype = :scatter)

B = [E[i]^j for i in 1:16, j in -1:0]

ê₂ = B \ b₂

b̂₂(E) = ê₂' * [E^-1,1]

ê₃ = B \ b₃

b̂₃(E) = ê₃' * [E^-1,1]

plot!(b̂₂, label = "b₂(e)")
plot!(b̂₃, label = "b₃(e)")

compₐ = []
erroₐ = []
a_C = [1623.500, 433.890, 39415.000, 25798.000, 284.200, 83.501, 2177.600, 1230.200]
a_AR = [ê₀[2], ê₂[2], ê₀[1], ê₂[1], ê₁[2], ê₃[2], ê₁[1], ê₃[1]]

for i in 1:8
    append!(compₐ, abs(a_AR[i]-a_C[i]))
    append!(erroₐ, round(compₐ[i]/a_C[i]*100, digits=2))
end

NC2733_AR = []
NC2733_C = []
NC2733_erro = []
NCX1116_AR = []
NCX1116_C = []
NCX1116_erro = []

for i in 1:4
    imp = 2i - 1
    par = 2i
    append!(NC2733_AR, a_AR[imp])
    append!(NC2733_C, a_C[imp])
    append!(NC2733_erro, erroₐ[imp])
    append!(NCX1116_AR, a_AR[par])
    append!(NCX1116_C, a_C[par])
    append!(NCX1116_erro, erroₐ[par])

end

using DataFrames

df = DataFrame(Coeficientes=["a₀", "a₁", "a₂", "a₃"],
    NC2733_Artigo=NC2733_AR, NC2733_Experimento=NC2733_C, NC2733_Erro=NC2733_erro,
    NCX1116_Artigo=NCX1116_AR, NCX1116_Experimento=NCX1116_C, NCX1116_Erro=NCX1116_erro)

show(df, allcols=true)
print("\n\n(Erros já em porcentagem)")
display(df, )

using PyCall
import CoolProp.CoolProp

CP = pyimport("CoolProp.CoolProp")

T = 90 + 273.15
P = 100

stt = CP.State("Air", {"T":T, "P":P})


text:
# Defines some commonly used states in quantum information

export bell, bellDiagState, isotropicState, rState, sState, sStateQutrit, rStatePhase, rStateCorrPhase, rStateCorrPhaseCopies;

""" `bell`

The Bell states (as vectors), bell[j,1:4] gives the j-th bell state.
"""
bell = zeros(4,4);
bell[1,1:4] = maxEntVec(2);
bell[2,1:4] = (kron(eVec(2,1),eVec(2,1)) -  kron(eVec(2,2),eVec(2,2)))/sqrt(2);
bell[3,1:4] = (kron(eVec(2,1),eVec(2,2)) +  kron(eVec(2,2),eVec(2,1)))/sqrt(2);
bell[4,1:4] = (kron(eVec(2,1),eVec(2,2)) -  kron(eVec(2,2),eVec(2,1)))/sqrt(2);

""" `epr` and `singlet`

Two of these have special names: we will use epr for the EPR pair, and singlet for the singlet.
"""
epr = bell[1,1:4];
singlet = bell[4,1:4];

""" `bellS`

An array of density matrices corresponding to the 4 Bell states.
"""
bS1 = epr*epr';
bS2 = bell[2,1:4]* bell[2,1:4]';
bS3 = bell[3,1:4]* bell[3,1:4]';
bS4 = bell[4,1:4]*bell[4,1:4]';
bellS = [bS1; bS2; bS3; bS4];

""" `PS`
Projector onto the symmetric subspace of 2 qubits.
"""
PS = bS1 + bS2 + bS3;

""" `PA`
Projector onto the antisymmetric subspace of 2 qubits.
"""
PA = bS4;

""" `rho = bellDiagState(p1,p2,p3)`

Returns a bell diagonal state that is a mixture between the 4 Bell states: p1 * phi^+ + p2 * psi^+ + p3 * phi^- + (1 - p1 - p2 - p3) * psi^-.
"""

function bellDiagState(p1::Number, p2::Number, p3::Number)

	@assert 0 <= p1 "Probilities must be positive."
	@assert 0 <= p2 "Probilities must be positive."
	@assert 0 <= p3 "Probilities must be positive."

	@assert p1 + p2 + p3 <= 1 "Probabilities cannot exceed 1."

	# Produce the desired mixture
	out = p1 * bS1 + p2 * bS3 + p3 * bS2 + (1-p1-p2-p3) * bS4;
	return out;
end

""" `rho = sState(p)`

Returns a state that is an a mixture between the EPR pair (with probability *p*), and the state |11><11|.
"""

function sState(p::Number)

	@assert 0 <= p "Probilities must be positive."
	@assert p <= 1 "Probabilities cannot exceed 1."

	# Generate the maximally entangled state
	epr = maxEnt(2);

	# Generate |11>
	e0 = [1 0];
	e1 = [0 1];
	v11 = kron(e1,e1);

	# Produce the desired mixture
	out = p * epr + (1-p) * v11'*v11;
	return out;
end

""" `rho = sStateQutrit(p)`

Returns a state that is a mixture between a 3 dimensional maximally entangled state (with probability *p*) and the state |00><00|
"""

function sStateQutrit(p::Number)

	@assert 0 <= p "Probilities must be positive."
	@assert p <= 1 "Probabilities cannot exceed 1."

	# Maximally entangled state
	epr = maxEnt(3);

	# Generate |00>
	e0 = [1 0 0];
	e1 = [0 0 1];
	v00 = kron(e0,e0);

	# Produce the desired mixture
	out = p * epr + (1-p) * v00'*v00;

  	return out;
end

""" `rho = rState(p)`

Returns a mixture between a state proportional to |01> + |10> (with probability *p*) and the state |11><11|.
"""

rState(p::Number) = rStatePhase(p, 0)

""" `rho = rStatePhase(p, phi)` or `rho = rStatePhase(p)`

Returns a mixture between a state proportional to |01> + e^(i *phi*) |10> (with probability *p*) and the state |11><11|. No value for *phi* defaults to *phi*=0.

"""

function rStatePhase(p::Number, phi::Number = 0.0)

	@assert 0 <= p "Probilities must be positive."
	@assert p <= 1 "Probabilities cannot exceed 1."

	# Produce the state |01> + e^(i phi) |10>
	e0 = eVec(2, 1);
	e1 = eVec(2, 2);
	vec = (kron(e0, e1) + e^(im*phi) * kron(e1,e0))/sqrt(2);

	if phi == 0 || phi == pi
		vec = real(vec);
	end

	# Produce a state orthogonal to the one above
	v11 = kron(e1,e1);

	# Construct the desire mixture
	out = p * vec*vec' + (1-p) * v11*v11';

  	return out
end

""" `rho = rStateCorrPhase(p, pd)`

Returns a state of the form integral phi [pd * r(p,phi) + (1-pd) * r(p,phi+pi)] tensor r(p,phi), where
r(p,phi) = rStatePhase(p,phi). Default for pd is 1.
"""
function rStateCorrPhase(p::Number, pd::Number = 1)

	@assert 0 <= p "Probilities must be positive."
	@assert p <= 1 "Probabilities cannot exceed 1."

    rStatePhaseDep(phi) = pd * rStatePhase(p, phi) + (1-pd) * rStatePhase(p, phi + pi)

	integrand(phi) = ( 1/(2*pi) ) * kron(rStatePhaseDep(phi), rStatePhase(p, phi));

	eps = 10.0^(-4)
	out = real( quadgk( integrand, 0, 2 * pi; reltol=sqrt(eps), abstol=0, maxevals=10^7, order=7, norm=vecnorm)[1])
  	return out
end

""" `rho = rStateCorrPhaseCopies(p)`

Returns a state of the form integral phi r(p,phi)^(tensor n) , where
r(p,phi) = rStatePhase(p,phi) and n is the number of desired copies.
"""
function rStateCorrPhaseCopies(p::Number, n::Int)

	@assert 0 <= p "Probilities must be positive."
	@assert p <= 1 "Probabilities cannot exceed 1."
	@assert n > 1 "Number of copies must be at least 1."

	integrand(phi) = ( 1/(2*pi) ) * copies(rStatePhase(p, phi), n);

	eps = 10.0^(-4)
	out = real( quadgk( integrand, 0, 2 * pi; reltol=sqrt(eps), abstol=0, maxevals=10^7, order=7, norm=vecnorm)[1])
  	return out
end

""" `rho = isotropicState(p)` or `rho = isotropicState(p,d)`

Returns an isotropic state, i.e., a mixture of a maximally entangled pair (with probability *p*) and the maximally mixed state in local dimension *d*. If no argument *d* is given the default is *d*=2, that is, we take the mixture of the EPR pair with the maximally mixed state of local dimension 2.
"""

function isotropicState(p::Number; d::Int = 2)

	@assert 0 <= p "Probilities must be positive."
	@assert p <= 1 "Probabilities cannot exceed 1."

        epr = maxEnt(d);
	out = p * epr*epr' + (1-p) * eye(d^2)/d^2;

	return out;
end

ext: jl


text:
#!/bin/bash
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
# below gets the directory name of the script, even if there is a symlink involved
# from https://stackoverflow.com/questions/59895/get-the-source-directory-of-a-bash-script-from-within-the-script-itself
SOURCE="${BASH_SOURCE[0]}"
while [ -h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
  DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
  SOURCE="$(readlink "$SOURCE")"
  [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
export JULIA_PROJECT=$DIR/..
export JULIA_LOAD_PATH=@:@stdlib  # exclude default environment
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "$SOURCE" "$@"
=#
# Parse command-line first, so a failure can be detected before the compilation
# or execution of any code unrelated to argument parsing.
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    s.description="""Analyze one or more LJH data files containing noise data,
to determine the noise autocorrelation, power-spectral density, and best-fit
ARMA model. Also store the results in an HDF5 file. By default, will fail if the
output HDF5 file already exists, though you can specify -r to replace an existing
file or -u to update it by addition of new channels (noise analysis for existing
channels cannot be replaced)."""

    @add_arg_table s begin
        "--outputfile", "-o"
            help = "store the results in OUTPUTFILE (an HDF5 file) instead of the inferred file or files"
            arg_type = String
        "--replaceoutput", "-r"
            help = "delete and replace any existing output files"
            action = :store_true
        "--nlags", "-n"
            help = "compute autocorrelation for this many lags (default: LJH record length)"
            arg_type = Int
            default = 0
        "--nfreq", "-f"
            help = "compute power spectrum for this many frequencies (default: LJH record length//2)"
            arg_type = Int
            default = 0
        "pulse_file"
            help = "name of the pulse containing single ljh noise file, all channels will be procssed"
            required = true
            arg_type = String
        "--dontcrash"
            help = "pass this to move on past any channel that fails"
            action = :store_true
        "--maxchannels"
            default = 10^4
            arg_type = Int
    end
    return parse_args(s)
end
parsed_args = parse_commandline()

using Pope
using Pope.LJH
if parsed_args["outputfile"]==nothing
    parsed_args["outputfile"] = Pope.outputname(parsed_args["pulse_file"],"noise")
end
display(parsed_args);println()
if !parsed_args["replaceoutput"] && isfile(parsed_args["outputfile"])
    println("intended output file $(parsed_args["outputfile"]) exists, pass --replaceoutput if you would like to overwrite it")
    exit(1) # anything other than 0 indicated process unsuccesful
end
using HDF5
using ARMA
using Pope.NoiseAnalysis
outputh5 = h5open(parsed_args["outputfile"],"w")
ljhdict = LJH.allchannels(parsed_args["pulse_file"],parsed_args["maxchannels"]) # ordered dict mapping channel number to filename
function analyze_one_file(filename::AbstractString, channum::Integer,
        outputh5::HDF5.HDF5File, nlags::Integer=0, nfreq::Integer=0;
        max_samples::Integer=50000000)

    # Open the LJH file for reading
    f = LJHFile(filename)
    frametime = LJH.frametime(f)
    nsamp = LJH.record_nsamples(f)
    nrec = LJH.ljh_number_of_records(f)
    if nsamp*nrec > max_samples
        nrec = max_samples // nsamp
    end
    rawdata = vcat([rec.data for rec in f[1:nrec]]...)
    samplesUsed = length(rawdata)

    if length(rawdata) == 0
        error("rawdata has length=0 for $filename")
    elseif all(rawdata.==rawdata[1])
        error("rawdata is all the same value $(rawdata[1]) for $filename")
    end

    if nlags <= 0
        nlags = nsamp
    end
    if nfreq <= 0
        nfreq = NoiseAnalysis.round_up_dft_length(div(nsamp, 2)) + 1
    end

    autocorr = compute_autocorr(rawdata, nlags, max_exc=1000)
    psd = compute_psd(rawdata, nfreq, frametime, max_exc=1000)
    freq = NoiseAnalysis.psd_freq(nfreq, frametime)
    freqstep = freq[2]-freq[1]

    max_ARMA_order = 5
    # outputh5["noise_analysis_autocorr_input/chan$channum"]=autocorr # save some data for debugging
    model = fitARMA(autocorr, max_ARMA_order, pmin=1)

    noise = NoiseResult(autocorr, psd, samplesUsed, freqstep, filename, model)
    NoiseAnalysis.hdf5save(outputh5, channum, noise)
end

if length(ljhdict) == 0
    println("Analyzed 0 files")
    exit()
end

for (channum, ljhname) in ljhdict
    println("Analyzing $ljhname")
    try
        analyze_one_file(ljhname, channum, outputh5, parsed_args["nlags"], parsed_args["nfreq"])
    catch ex
        if !parsed_args["dontcrash"]
            rethrow(ex)
        else
            println("$ljhname FAILED noise analysis")
            println(ex)
        end
    end
end


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear
Random.seed!(100)
prob = prob_sde_linear
f(u,p,t) = 2u
prob = SDEProblem{false}(f,prob.f,prob.u0,prob.tspan)
sol = solve(prob,EM(),dt=1//2^(4),saveat = 0.33)

sol.t == [0.0,.33,.33+.33,1]

sol = solve(prob,EM(),dt=1//2^(4),saveat = 0.33,save_start=false)

sol.t == [.33,.33+.33,1]
sol(0.4)

sol = solve(prob,EM(),dt=1//2^(4),saveat = 0.33,save_start=false,save_everystep=true)

@test 0.33 ∈ sol.t
@test 0.66 ∈ sol.t
sol.t != [.33,.33+.33,1]

prob = prob_sde_2Dlinear
f(du,u,p,t) = prob.f(du,u,p,t)
prob2 = SDEProblem(f,prob.g,vec(prob.u0),prob.tspan)
Random.seed!(200)
sol = solve(prob2,EM(),dt=1//2^(4),saveat = 0.33,save_start=false,save_everystep=true)
Random.seed!(200)
sol2= solve(prob2,EM(),dt=1//2^(4),saveat = 0.33,save_start=false,save_everystep=true,save_idxs=1:2:5)

@test sol(0.43)[1:2:5] == sol2(0.43)

# Test save_on switch
sol = solve(prob2, EM(), dt=1//2^(4), save_on=false, save_start=false, save_end=false)
@test isempty(sol.t) && isempty(sol.u)
sol = solve(prob2, EM(), dt=1//2^(4), saveat=0.33, save_on=false, save_start=false, save_end=false)
@test isempty(sol.t) && isempty(sol.u)


text:
#=
Example 2:
  Fit the following analytical model with 5 parameters:
  (p1  +  p2 * x  +  p3 * x^2  +  p4 * sin(p5 * x))  *  cos(x)

A `FuncWrap` component is used to wrap a common Julia function
=#
using Random, DataFitting

# Define the analytic model function
f(x, p1, p2, p3, p4, p5) = @. (p1  +  p2 * x  +  p3 * x^2  +  p4 * sin(p5 * x))  *  cos(x)

# "True" parameter values:
params = [1, 1.e-3, 1.e-6, 4, 5]

# "True" values
x = 1.:5:5000
y = f(x, params...);

# Add some noise to simulate a measurement process
noise = 1.
y .+= noise .* randn(MersenneTwister(0), length(x));

# Prepare:
# - domain for model evaluation
# - empirical data
data = Measures(y, noise)

# Build a model with appropriate components and parameter guess values
model = Model(:comp1 => FuncWrap(f, params...))

# Add an instrument (on the model domain), and an expression involving
# defined components:
add_dom!(model, x)
addexpr!(model, 1, :comp1)

# Fit the model to empirical data
result = fit!(model, data)

# Print best fit parameters
for i in 1:length(result.comp1.p)
    println("p[$i] = ", result.comp1.p[i].val, " ± ", result.comp1.p[i].unc)
end

# Plot data and best fit model
using Gnuplot
@gp    model[].domain data.val :-
@gp :- model[].domain model[].expr1 "w line"


text:
clf()
nmed = 10
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdq0[i,2:count[i],k,l])) > 0
#    diff1 = abs.(asinh.(dtdq0_big[i,2:count[i],k,l])-asinh.(dtdq0[i,2:count[i],k,l]));
#    diff1 = abs.(dtdq0_big[i,2:count[i],k,l]-dtdq0[i,2:count[i],k,l])/maximum(dtdq0[i,2:count[i],k,l]);
    diff1 = abs.(dtdq0_big[i,2:count[i],k,l]-dtdq0[i,2:count[i],k,l])
    diff2 = copy(diff1);
    diff3 = abs.(dtdelements_big[i,2:count[i],k,l]-dtdelements_mixed[i,2:count[i],k,l])
#    diff3 = abs.(asinh.(dtdelements_big[i,2:count[i],k,l])-asinh.(dtdelements_mixed[i,2:count[i],k,l]))
    diff4 = copy(diff3);
    ntt1 = size(diff1)[1];
    for it=nmed:ntt1-nmed
      # Divide by a median smoothed dtdq0:
      diff2[it] = diff1[it]/maximum(abs.(dtdq0[i,it+1-nmed+1:it+1+nmed,k,l]))
      maxit = maximum(abs.(dtdelements_mixed[i,it+1-nmed+1:it+1+nmed,k,l]))
      if maxit > 0
        diff4[it] = diff3[it]/maximum(abs.(dtdelements_mixed[i,it+1-nmed+1:it+1+nmed,k,l]))
      end
    end;
    diff2[1:nmed-1] .= diff2[nmed]; diff2[ntt1-nmed+1:ntt1] .= diff2[ntt1-nmed];
    diff4[1:nmed-1] .= diff4[nmed]; diff4[ntt1-nmed+1:ntt1] .= diff4[ntt1-nmed];
#    diff3 = abs.(asinh.(dtdq2_big[i,2:count[i],k,l])-asinh.(dtdq2[i,2:count[i],k,l]));
    loglog(tt[i,2:count[i]]-tt[i,1],diff2,linestyle=":");
    loglog(tt[i,2:count[i]]-tt[i,1],diff4,linestyle="--");
#    loglog(tt[i,2:count[i]]-tt[i,1],diff3);
    println("xvs: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff2))); #read(STDIN,Char);
    println("els: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff4))); #read(STDIN,Char);
#    println(i," ",k," ",l," asinh error h/2: ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
end
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5)

clf()
nmed = 10
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdq0[i,2:count[i],k,l])) > 0
    diff1 = abs.(asinh.(dtdq0_big[i,2:count[i],k,l])-asinh.(dtdq0[i,2:count[i],k,l]));
    diff3 = abs.(asinh.(dtdelements_big[i,2:count[i],k,l])-asinh.(dtdelements_mixed[i,2:count[i],k,l]))
    loglog(tt[i,2:count[i]]-tt[i,1],diff1,linestyle=":");
    loglog(tt[i,2:count[i]]-tt[i,1],diff3,linestyle="--");
    println("xvs: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff1))); #read(STDIN,Char);
    println("els: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
end
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5)

@inbounds for i=1:n, j=1:count[i]
  @inbounds for k=1:n, l=1:7; dtdelements_mixed[i,j,l,k] = 0.0
    @inbounds for p=1:n, q=1:7
      dtdelements_mixed[i,j,l,k] += dtdq0[i,j,q,p]*convert(Float64,jac_init_big[(p-1)*7+q,(k-1)*7+l])
    end
  end
end

dtdelements_mixed_big = copy(dtdelements_big)
@inbounds for i=1:n, j=1:count[i]
  @inbounds for k=1:n, l=1:7; dtdelements_mixed_big[i,j,l,k] = 0.0
    @inbounds for p=1:n, q=1:7
      dtdelements_mixed_big[i,j,l,k] += dtdq0_big[i,j,q,p]*jac_init_big[(p-1)*7+q,(k-1)*7+l]
    end
  end
end

clf()
nmed = 10
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdq0[i,2:count[i],k,l])) > 0
    diff1 = abs.(dtdq0_big[i,2:count[i],k,l]-dtdq0[i,2:count[i],k,l])/median(abs.(dtdq0[i,2:count[i],k,l]))
    loglog(tt[i,2:count[i]]-tt[i,1],diff1,linestyle=":");
    println("xvs: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff1))); #read(STDIN,Char);
  end
  if maximum(abs.(dtdelements[i,2:count[i],k,l])) > 0
    diff3 = abs.(dtdelements_big[i,2:count[i],k,l]-dtdelements[i,2:count[i],k,l])/median(abs.(dtdelements[i,2:count[i],k,l]));
    loglog(tt[i,2:count[i]]-tt[i,1],diff3,linestyle="--");
    println("els: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
end
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5)


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

export linesearch

linesearch(ls, sol::result, d::AbstractArray{T}, f::Function, g!::Function, fg!::Function,
           t, funRef, gtd, gvec::AbstractArray{T}) where T =
        linesearch(ls, sol, d, f, g!, fg!, T(t), T(funRef), T(gtd), gvec)

"""
    linesearch(ls, sol, d, f, g!, fg!, t, funRef, gtd, gvec)

Line search interface to LineSearches.jl

# Arguments

- `ls`: Line search structure (see LineSearches.jl documentation)
- `f`: Objective function, x -> f(x)
- `g!``: Gradient in place function, x-> (g.= gradient(x))
- `fg!``: Objective and in place gradient function, x-> (f = f(x); g.= gradient(x))
- `t`: Initial steplength gess
- `funRef`: Reference objective function value
- `gtd`: Reference direction inner product dot(g0, d0)
- `gvec`: prealocated array for thhe gradient

"""
function linesearch(ls, sol::result, d::AbstractArray{T}, f::Function, g!::Function, fg!::Function,
                    t::T, funRef::T, gtd::T, gvec::AbstractArray{T}) where T
    # Univariate line search functions
    ϕ(α) = f(sol.x .+ α.*d)

    function dϕ(α)
        g!(gvec, sol.x .+ α.*d)
        return dot(gvec, d)
    end

    function ϕdϕ(α)
        phi = fg!(gvec, sol.x .+ α.*d)
        dphi = dot(gvec, d)
        return (phi, dphi)
    end

    # Line search. Prevents it to throw error.
    try
        return ls(ϕ, dϕ, ϕdϕ, t, funRef, gtd)
    catch e
        @info "Line search failed"
        return 0, funRef
    end
end


text:
using LinearAlgebra

struct PDS
    ∇h
    proxF
    proxG
    L
    Lᵀ

    iterations::Integer
    τ::Real
    σ::Real
    λ::Real

    absTol::Real
end

function PDS(∇h; proxF=(t,a)->t, proxG=(t,a)->t, L=identity, adjL=identity, τ=1.0, σ=1.0, λ=1.0, iterations=nothing, absTol=1e2*eps())
    PDS(∇h, proxF, proxG, L, adjL, iterations, τ, σ, λ, absTol)
end

# FB-based primal-dual splitting algorithm
# N. Komodakis and J. Pesquet, "Playing with Duality," IEEE Signal Processing Magazine, vol. 32, no. 6, pp. 31--54, 2015.
# Algorithm 3
# argmin_{x} f(x) + g(Lx) + h(x)
# f() and g(): proxiable function
# h(x): smooth convex function having a Lipschitzian gradient
# ∇loss : gradient of loss function f(x)
# L: linear operator
# Lᵀ: adjoint of L
function (pds::PDS)(x0, v0=pds.L(x0); iterations=pds.iterations, verboseFunction::Function=(args...)->nothing)
    xₖ = x0
    v = v0
    cproxG = fmconj(pds.proxG)
    for k = 1:iterations
        xₖ₋₁ = xₖ
        p = pds.proxF(xₖ - pds.τ*(pds.∇h(xₖ) + pds.Lᵀ(v)), pds.τ)
        q = cproxG(v + pds.L(2.0*p - xₖ₋₁), pds.σ^-1)

        # (x, v) <- (x, v) + λ((p,q) - (x,v)) for　λ == 1
        xₖ, v = (x, v) .+ λ .* (p - x, q - v)
        # xₖ, v = p, q

        errx = norm(xₖ - xₖ₋₁)^2 /2
        verboseFunction(nItr, xₖ, errx)
        # TODO: break condition
    end
    (xₖ, ())
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
abstract type AbstractDirection end

struct Up <: AbstractDirection end
const UP = Up()
move(::Up, pos::CartesianIndex{2}) = pos + CartesianIndex(-1, 0)

struct Down <: AbstractDirection end
const DOWN = Down()
move(::Down, pos::CartesianIndex{2}) = pos + CartesianIndex(1, 0)

struct Left <: AbstractDirection end
const LEFT = Left()
move(::Left, pos::CartesianIndex{2}) = pos + CartesianIndex(0, -1)

struct Right <: AbstractDirection end
const RIGHT = Right()
move(::Right, pos::CartesianIndex{2}) = pos + CartesianIndex(0, 1)

const DIRECTIONS = (UP, DOWN, LEFT, RIGHT)


text:
module TouchOSC

export Layout, Page
export Led, Push, Toggle, Fader, Label, Encoder
export VERTICAL, HORIZONTAL
export IPHONE, IPAD, IPHONE5, IPHONE6, IPHONE6PLUS, IPADPRO
export RED, GREEN, BLUE, YELLOW, PURPLE, GRAY, ORANGE, BROWN, PINK
export RELATIVE, ABSOLUTE
export save

using LightXML
using ZipFile

const TOUCHOSC_FILE_VERSION = 15

# ascii(base64decode("LzEvdGVzdA=="))
# base64encode("/1/test")

@enum Orientation VERTICAL HORIZONTAL

# I haven't checked all these against the editor-generated XML, I'm assuming
# based on checking the first two
@enum LayoutMode IPHONE IPAD IPHONE5 IPHONE6 IPHONE6PLUS IPADPRO
@enum ControlColor RED GREEN BLUE YELLOW PURPLE GRAY ORANGE BROWN PINK

abstract Control

"The properties common to all Controls"
type ControlInfo
    name::ASCIIString
    position::Tuple{Int, Int}
    color::ControlColor
    size::Tuple{Int, Int}
    oscpath::ASCIIString
end

function ControlInfo(; name="", position=(0, 0), color=RED, size=(50, 50), oscpath="")
    ControlInfo(name, position, color, size, oscpath)
end

# create the accessor functions
for prop in [:name, :position, :color, :oscpath]
    @eval $prop(ctrl::Control) = ctrl.info.$prop
end

# size is special because we're extending Base.size
Base.size(ctrl::Control) = ctrl.info.size

type Page
    name::ASCIIString
    controls::Vector{Control}
end

function Page(;name="", controls=Control[])
    Page(name, controls)
end

type Layout
    orientation::Orientation
    mode::LayoutMode
    pages::Vector{Page}
end

function Layout(;orientation=VERTICAL, mode=IPHONE, pages=Page[])
    Layout(orientation, mode, pages)
end

function pushxml!(xmldoc, layout::Layout)
    root = create_root(xmldoc, "layout")
    set_attribute(root, "mode", string(Int(layout.mode)))
    # in the XML the string is reversed from what's displayed in the editor
    set_attribute(root, "orientation",
        layout.orientation == VERTICAL ? "horizontal" : "vertical")
    set_attribute(root, "version", TOUCHOSC_FILE_VERSION)

    for page in layout.pages
        # we have to give the layout because some of the XML generation
        # varies depending on the layout
        pushxml!(root, page, layout)
    end

    root
end

# needs to include the layout object so we can inspect it during XML generation
# for some of the controls
function pushxml!(layoutxml, page::Page, layout::Layout)
    child = new_child(layoutxml, "tabpage")
    set_attribute(child, "name", base64encode(page.name))
    set_attribute(child, "scalef", "0.0")
    set_attribute(child, "scalet", "1.0")

    for control in page.controls
        pushxml!(child, control, layout)
    end

    child
end

# transforms an XY position and size into the correct coordinates for the XML,
# which involves flipping X and Y axes in landscape mode, and inverting the X
# axis. This means that when we create the objects we can use locations that
# match what we'd see in the editor, and they'll end up in the right places
function transform(layout::Layout, pos, siz)
    if layout.orientation == VERTICAL
        return pos, siz
    end
    # take into account the space for the top bar
    screenheight = size(layout)[2] - 40

    (screenheight - pos[2] - siz[2], pos[1]), (siz[2], siz[1])
end

function pushxml!(pagexml, ctrl::Control, layout::Layout)
    landscape = layout.orientation == HORIZONTAL
    layoutpos, layoutsize = transform(layout, position(ctrl), size(ctrl))
    xml = new_child(pagexml, "control")
    set_attribute(xml, "name", base64encode(name(ctrl)))
    set_attribute(xml, "x", layoutpos[1])
    set_attribute(xml, "y", layoutpos[2])
    set_attribute(xml, "w", layoutsize[1])
    set_attribute(xml, "h", layoutsize[2])
    set_attribute(xml, "color", lowercase(string(color(ctrl))))
    if oscpath(ctrl) != ""
        set_attribute(xml, "osc_cs", base64encode(oscpath(ctrl)))
    end
    set_attribute(xml, "type", typestring(ctrl, landscape))
    # now add the type-specific info
    addxmlinfo!(xml, ctrl)
end

Base.push!(layout::Layout, page::Page) = push!(layout.pages, page)
Base.push!(page::Page, ctrl::Control) = push!(page.controls, ctrl)

function save(layout::Layout, filename)
    zip = ZipFile.Writer(filename)
    index = ZipFile.addfile(zip, "index.xml")

    xmldoc = XMLDocument()
    pushxml!(xmldoc, layout)
    print(xmldoc)
    write(index, string(xmldoc))
    close(zip)
end

const layoutsizes = [
    (480, 320), # IPHONE
    (1024, 768), # IPAD
    (568, 320), # IPHONE5
    (667, 375), # IPHONE6
    (736, 414), # IPHONE6PLUS
    (1366, 1024) # IPADPRO
]
function Base.size(layout::Layout)
    dims = layoutsizes[Int(layout.mode)+1]

    layout.orientation == HORIZONTAL ? dims : (dims[2], dims[1])
end

include("controls.jl")

end # module


text:
"""
    ngwp_analysis(G::GraphSig, wavelet_packet::Array{Float64,3})

For a GraphSig object `G`, generate the matrix of NGWP expansion coefficients.

# Input Arguments
- `G::GraphSig`: an input GraphSig object
- `wavelet_packet::Array{Float64,3}`: the varimax wavelets packet.

# Output Argument
- `dmatrix::Array{Float64,3}`: the expansion coefficients matrix.

"""
function ngwp_analysis(G::GraphSig, wavelet_packet::Array{Float64,3})
    f = G.f
    fcols = size(f, 2)
    (N, jmax, ) = Base.size(wavelet_packet)

    dmatrix = zeros(N, jmax, fcols)
    dmatrix[:, 1, :] = f

    for j = 2:jmax
        for i = 1:N
            dmatrix[i, j, :] = f' * wavelet_packet[i, j, :]
        end
    end

    return dmatrix
end

"""
    const_proj_wavelets(𝚽,vlist,elist; method = "Modified Gram-Schmidt with Lp Pivoting")

construct projection wavelets, i.e., project δ ∈ vlist onto span({φⱼ| j ∈ elist}).

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors 𝚽
- `vlist::Array{Int}`: the list of considered node indices.
- `elist::Array{Int}`: the list of considered eigenvector indices.
- `method::Symbol`: default is `:MGSLp`. other options: `:IP` (Iterative-Projection),
    `GS` (Gram Schmidt).

# Output Argument
- `Wav::Matrix{Float64}`: a matrix whose columns are projected wavelet vectors.

"""
function const_proj_wavelets(𝚽, vlist, elist; method = :MGSLp)
    if length(vlist) == 1
        return 𝚽[:, elist]
    end
    N = size(𝚽, 1)
    m = length(vlist)
    Wav = zeros(N, m)

    B = 𝚽[:, elist]

    if method == :IP
        for k in 1:length(vlist)
            wavelet = Proj(spike(vlist[k], N), B)
            Wav[:, k] .= wavelet ./ norm(wavelet)
            B = wavelet_perp_Matrix(wavelet, B)
        end
    elseif method == :GS
        P = B * B'
        for k in 1:length(vlist)
            wavelet = P * spike(vlist[k], N)
            Wav[:,k] .= wavelet ./ norm(wavelet)
        end
        Wav, complement_dim = gram_schmidt(Wav)
        if complement_dim != 0
            complement_space = B * nullspace(Wav' * B)
            Wav = hcat(Wav, complement_space)
        end
    elseif method == :MGSLp
        P = B * B'
        for k in 1:length(vlist)
            wavelet = P * spike(vlist[k], N)
            Wav[:,k] .= wavelet ./ norm(wavelet)
        end
        Wav, complement_dim = mgslp(Wav)
        if complement_dim != 0
            complement_space = B * nullspace(Wav' * B)
            Wav = hcat(Wav, complement_space)
        end
    else
        error("Do not support method ", method)
    end

    return Wav
end

"""
    function NGWP_jkl(GP_star::GraphPart, drow::Int, dcol::Int)

Generate the (j,k,l) indices for the NGWP basis vector corresponding to the coefficient dmatrix(drow,dcol)

### Input Arguments
* `GP_star::GraphPart`: a GraphPart object of the dual grpah
* `drow::Int`: the row of the expansion coefficient
* `dcol::Int`: the column of the expansion coefficient

### Output Argument
* `j`: the level index of the expansion coefficient
* `k`: the subregion in dual graph's index of the expansion coefficient
* `l`: the tag of the expansion coefficient
"""
function NGWP_jkl(GP_star::GraphPart, drow::Int, dcol::Int)
    (j, k, l) = GHWT_jkl(GP_star, drow, dcol)
    return j,k,l
end


text:
using Ripemd
if VERSION < v"1.0.0"
    using Base.Test
else
    using Test
end

files = readdir(".")

# use only files that are named test_*.jl and start with the newest

if VERSION < v"0.7.0"
    filter!(x -> contains(x, r"^test_"), files)
else
    filter!(x -> occursin(r"^test_", x), files)
end

sort!(files, by = x -> stat(x).mtime, rev = true)

for f in files
    include(f)
end


text:
module WhereTraits
export @traits, @traits_test, @traits_get_state, @traits_show_implementation

using Compat

@Base.kwdef mutable struct _Config
  # set to false because `@suppress` does currently not work with function definitions
  # for updates see https://github.com/JuliaIO/Suppressor.jl/issues/29
  suppress_on_traits_definitions::Bool = false
  auto_documentation::Bool = true
end
# TODO documentation of struct fields does not seem to work - hence we document the constant instead
"""
suppress_on_traits_definitions::Bool = false

    if true, all warnings, e.g. overwrite warnings are suppressed within @traits definitions

    this is useful, because in the current Julia, @traits is inherintly state-ful and edge-cases
    may need to overwrite previous definitions of the same @traits
"""
const CONFIG = _Config()

include("Utils/Utils.jl")
using .Utils

# this syntax is so complex that we need to store a state for each function
# which because of its complexity is defined in an extra submodule
include("InternalState.jl")
using .InternalState

include("Syntax/Syntax.jl")
using .Syntax

include("ExtraHelpers.jl")
using .ExtraHelpers

include("BasicTraits.jl")

end # module


text:
function FieldAbsoluteDate(arg0::AbsoluteDate, arg1::RealFieldElement)
    return FieldAbsoluteDate((AbsoluteDate, RealFieldElement), arg0, arg1)
end

function FieldAbsoluteDate(arg0::Field)
    return FieldAbsoluteDate((Field,), arg0)
end

function FieldAbsoluteDate(arg0::Field, arg1::AbsoluteDate)
    return FieldAbsoluteDate((Field, AbsoluteDate), arg0, arg1)
end

function FieldAbsoluteDate(arg0::Field, arg1::Date, arg2::TimeScale)
    return FieldAbsoluteDate((Field, Date, TimeScale), arg0, arg1, arg2)
end

function FieldAbsoluteDate(arg0::Field, arg1::DateComponents, arg2::TimeComponents, arg3::TimeScale)
    return FieldAbsoluteDate((Field, DateComponents, TimeComponents, TimeScale), arg0, arg1, arg2, arg3)
end

function FieldAbsoluteDate(arg0::Field, arg1::DateComponents, arg2::TimeScale)
    return FieldAbsoluteDate((Field, DateComponents, TimeScale), arg0, arg1, arg2)
end

function FieldAbsoluteDate(arg0::Field, arg1::DateTimeComponents, arg2::TimeScale)
    return FieldAbsoluteDate((Field, DateTimeComponents, TimeScale), arg0, arg1, arg2)
end

function FieldAbsoluteDate(arg0::Field, arg1::JString, arg2::TimeScale)
    return FieldAbsoluteDate((Field, JString, TimeScale), arg0, arg1, arg2)
end

function FieldAbsoluteDate(arg0::Field, arg1::jint, arg2::Month, arg3::jint, arg4::TimeScale)
    return FieldAbsoluteDate((Field, jint, Month, jint, TimeScale), arg0, arg1, arg2, arg3, arg4)
end

function FieldAbsoluteDate(arg0::Field, arg1::jint, arg2::Month, arg3::jint, arg4::jint, arg5::jint, arg6::jdouble, arg7::TimeScale)
    return FieldAbsoluteDate((Field, jint, Month, jint, jint, jint, jdouble, TimeScale), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function FieldAbsoluteDate(arg0::Field, arg1::jint, arg2::jint, arg3::jint, arg4::TimeScale)
    return FieldAbsoluteDate((Field, jint, jint, jint, TimeScale), arg0, arg1, arg2, arg3, arg4)
end

function FieldAbsoluteDate(arg0::Field, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jint, arg6::jdouble, arg7::TimeScale)
    return FieldAbsoluteDate((Field, jint, jint, jint, jint, jint, jdouble, TimeScale), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function FieldAbsoluteDate(arg0::FieldAbsoluteDate, arg1::RealFieldElement)
    return FieldAbsoluteDate((FieldAbsoluteDate, RealFieldElement), arg0, arg1)
end

function FieldAbsoluteDate(arg0::FieldAbsoluteDate, arg1::jdouble)
    return FieldAbsoluteDate((FieldAbsoluteDate, jdouble), arg0, arg1)
end

function FieldAbsoluteDate(arg0::FieldAbsoluteDate, arg1::jdouble, arg2::TimeScale)
    return FieldAbsoluteDate((FieldAbsoluteDate, jdouble, TimeScale), arg0, arg1, arg2)
end

function compare_to(obj::FieldAbsoluteDate, arg0::FieldAbsoluteDate)
    return jcall(obj, "compareTo", jint, (FieldAbsoluteDate,), arg0)
end

function create_besselian_epoch(::Type{FieldAbsoluteDate}, arg0::RealFieldElement)
    return jcall(FieldAbsoluteDate, "createBesselianEpoch", FieldAbsoluteDate, (RealFieldElement,), arg0)
end

function create_besselian_epoch(::Type{FieldAbsoluteDate}, arg0::RealFieldElement, arg1::TimeScales)
    return jcall(FieldAbsoluteDate, "createBesselianEpoch", FieldAbsoluteDate, (RealFieldElement, TimeScales), arg0, arg1)
end

function create_gps_date(::Type{FieldAbsoluteDate}, arg0::jint, arg1::RealFieldElement)
    return jcall(FieldAbsoluteDate, "createGPSDate", FieldAbsoluteDate, (jint, RealFieldElement), arg0, arg1)
end

function create_gps_date(::Type{FieldAbsoluteDate}, arg0::jint, arg1::RealFieldElement, arg2::TimeScale)
    return jcall(FieldAbsoluteDate, "createGPSDate", FieldAbsoluteDate, (jint, RealFieldElement, TimeScale), arg0, arg1, arg2)
end

function create_jd_date(::Type{FieldAbsoluteDate}, arg0::jint, arg1::RealFieldElement, arg2::TimeScale)
    return jcall(FieldAbsoluteDate, "createJDDate", FieldAbsoluteDate, (jint, RealFieldElement, TimeScale), arg0, arg1, arg2)
end

function create_julian_epoch(::Type{FieldAbsoluteDate}, arg0::RealFieldElement)
    return jcall(FieldAbsoluteDate, "createJulianEpoch", FieldAbsoluteDate, (RealFieldElement,), arg0)
end

function create_julian_epoch(::Type{FieldAbsoluteDate}, arg0::RealFieldElement, arg1::TimeScales)
    return jcall(FieldAbsoluteDate, "createJulianEpoch", FieldAbsoluteDate, (RealFieldElement, TimeScales), arg0, arg1)
end

function create_mjd_date(::Type{FieldAbsoluteDate}, arg0::jint, arg1::RealFieldElement, arg2::TimeScale)
    return jcall(FieldAbsoluteDate, "createMJDDate", FieldAbsoluteDate, (jint, RealFieldElement, TimeScale), arg0, arg1, arg2)
end

function duration_from(obj::FieldAbsoluteDate, arg0::AbsoluteDate)
    return jcall(obj, "durationFrom", RealFieldElement, (AbsoluteDate,), arg0)
end

function duration_from(obj::FieldAbsoluteDate, arg0::FieldAbsoluteDate)
    return jcall(obj, "durationFrom", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function equals(obj::FieldAbsoluteDate, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_arbitrary_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getArbitraryEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_ccsds_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getCCSDSEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_components(obj::FieldAbsoluteDate, arg0::TimeScale)
    return jcall(obj, "getComponents", DateTimeComponents, (TimeScale,), arg0)
end

function get_components(obj::FieldAbsoluteDate, arg0::TimeZone)
    return jcall(obj, "getComponents", DateTimeComponents, (TimeZone,), arg0)
end

function get_components(obj::FieldAbsoluteDate, arg0::TimeZone, arg1::TimeScale)
    return jcall(obj, "getComponents", DateTimeComponents, (TimeZone, TimeScale), arg0, arg1)
end

function get_components(obj::FieldAbsoluteDate, arg0::jint)
    return jcall(obj, "getComponents", DateTimeComponents, (jint,), arg0)
end

function get_components(obj::FieldAbsoluteDate, arg0::jint, arg1::TimeScale)
    return jcall(obj, "getComponents", DateTimeComponents, (jint, TimeScale), arg0, arg1)
end

function get_date(obj::FieldAbsoluteDate)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_field(obj::FieldAbsoluteDate)
    return jcall(obj, "getField", Field, ())
end

function get_fifties_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getFiftiesEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_future_infinity(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getFutureInfinity", FieldAbsoluteDate, (Field,), arg0)
end

function get_galileo_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getGalileoEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_gps_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getGPSEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_j2000_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getJ2000Epoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_java_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getJavaEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_julian_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getJulianEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_modified_julian_epoch(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getModifiedJulianEpoch", FieldAbsoluteDate, (Field,), arg0)
end

function get_past_infinity(::Type{FieldAbsoluteDate}, arg0::Field)
    return jcall(FieldAbsoluteDate, "getPastInfinity", FieldAbsoluteDate, (Field,), arg0)
end

function hash_code(obj::FieldAbsoluteDate)
    return jcall(obj, "hashCode", jint, ())
end

function is_after(obj::FieldAbsoluteDate, arg0::FieldTimeStamped)
    return jcall(obj, "isAfter", jboolean, (FieldTimeStamped,), arg0)
end

function is_after_or_equal_to(obj::FieldAbsoluteDate, arg0::FieldTimeStamped)
    return jcall(obj, "isAfterOrEqualTo", jboolean, (FieldTimeStamped,), arg0)
end

function is_before(obj::FieldAbsoluteDate, arg0::FieldTimeStamped)
    return jcall(obj, "isBefore", jboolean, (FieldTimeStamped,), arg0)
end

function is_before_or_equal_to(obj::FieldAbsoluteDate, arg0::FieldTimeStamped)
    return jcall(obj, "isBeforeOrEqualTo", jboolean, (FieldTimeStamped,), arg0)
end

function is_between(obj::FieldAbsoluteDate, arg0::FieldTimeStamped, arg1::FieldTimeStamped)
    return jcall(obj, "isBetween", jboolean, (FieldTimeStamped, FieldTimeStamped), arg0, arg1)
end

function is_between_or_equal_to(obj::FieldAbsoluteDate, arg0::FieldTimeStamped, arg1::FieldTimeStamped)
    return jcall(obj, "isBetweenOrEqualTo", jboolean, (FieldTimeStamped, FieldTimeStamped), arg0, arg1)
end

function is_close_to(obj::FieldAbsoluteDate, arg0::FieldTimeStamped, arg1::jdouble)
    return jcall(obj, "isCloseTo", jboolean, (FieldTimeStamped, jdouble), arg0, arg1)
end

function is_equal_to(obj::FieldAbsoluteDate, arg0::FieldTimeStamped)
    return jcall(obj, "isEqualTo", jboolean, (FieldTimeStamped,), arg0)
end

function offset_from(obj::FieldAbsoluteDate, arg0::FieldAbsoluteDate, arg1::TimeScale)
    return jcall(obj, "offsetFrom", RealFieldElement, (FieldAbsoluteDate, TimeScale), arg0, arg1)
end

function parse_ccsds_calendar_segmented_time_code(obj::FieldAbsoluteDate, arg0::jbyte, arg1::Vector{jbyte})
    return jcall(obj, "parseCCSDSCalendarSegmentedTimeCode", FieldAbsoluteDate, (jbyte, Vector{jbyte}), arg0, arg1)
end

function parse_ccsds_calendar_segmented_time_code(obj::FieldAbsoluteDate, arg0::jbyte, arg1::Vector{jbyte}, arg2::TimeScale)
    return jcall(obj, "parseCCSDSCalendarSegmentedTimeCode", FieldAbsoluteDate, (jbyte, Vector{jbyte}, TimeScale), arg0, arg1, arg2)
end

function parse_ccsds_day_segmented_time_code(::Type{FieldAbsoluteDate}, arg0::Field, arg1::jbyte, arg2::Vector{jbyte}, arg3::DateComponents)
    return jcall(FieldAbsoluteDate, "parseCCSDSDaySegmentedTimeCode", FieldAbsoluteDate, (Field, jbyte, Vector{jbyte}, DateComponents), arg0, arg1, arg2, arg3)
end

function parse_ccsds_day_segmented_time_code(::Type{FieldAbsoluteDate}, arg0::Field, arg1::jbyte, arg2::Vector{jbyte}, arg3::DateComponents, arg4::TimeScale)
    return jcall(FieldAbsoluteDate, "parseCCSDSDaySegmentedTimeCode", FieldAbsoluteDate, (Field, jbyte, Vector{jbyte}, DateComponents, TimeScale), arg0, arg1, arg2, arg3, arg4)
end

function parse_ccsds_unsegmented_time_code(::Type{FieldAbsoluteDate}, arg0::Field, arg1::jbyte, arg2::jbyte, arg3::Vector{jbyte}, arg4::FieldAbsoluteDate)
    return jcall(FieldAbsoluteDate, "parseCCSDSUnsegmentedTimeCode", FieldAbsoluteDate, (Field, jbyte, jbyte, Vector{jbyte}, FieldAbsoluteDate), arg0, arg1, arg2, arg3, arg4)
end

function parse_ccsds_unsegmented_time_code(::Type{FieldAbsoluteDate}, arg0::Field, arg1::jbyte, arg2::jbyte, arg3::Vector{jbyte}, arg4::FieldAbsoluteDate, arg5::FieldAbsoluteDate)
    return jcall(FieldAbsoluteDate, "parseCCSDSUnsegmentedTimeCode", FieldAbsoluteDate, (Field, jbyte, jbyte, Vector{jbyte}, FieldAbsoluteDate, FieldAbsoluteDate), arg0, arg1, arg2, arg3, arg4, arg5)
end

function shifted_by(obj::FieldAbsoluteDate, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldAbsoluteDate, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldAbsoluteDate, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldAbsoluteDate, (jdouble,), arg0)
end

function time_scales_offset(obj::FieldAbsoluteDate, arg0::TimeScale, arg1::TimeScale)
    return jcall(obj, "timeScalesOffset", RealFieldElement, (TimeScale, TimeScale), arg0, arg1)
end

function to_absolute_date(obj::FieldAbsoluteDate)
    return jcall(obj, "toAbsoluteDate", AbsoluteDate, ())
end

function to_date(obj::FieldAbsoluteDate, arg0::TimeScale)
    return jcall(obj, "toDate", Date, (TimeScale,), arg0)
end

function to_string(obj::FieldAbsoluteDate)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::FieldAbsoluteDate, arg0::TimeScale)
    return jcall(obj, "toString", JString, (TimeScale,), arg0)
end

function to_string(obj::FieldAbsoluteDate, arg0::TimeZone)
    return jcall(obj, "toString", JString, (TimeZone,), arg0)
end

function to_string(obj::FieldAbsoluteDate, arg0::TimeZone, arg1::TimeScale)
    return jcall(obj, "toString", JString, (TimeZone, TimeScale), arg0, arg1)
end

function to_string(obj::FieldAbsoluteDate, arg0::jint)
    return jcall(obj, "toString", JString, (jint,), arg0)
end

function to_string(obj::FieldAbsoluteDate, arg0::jint, arg1::TimeScale)
    return jcall(obj, "toString", JString, (jint, TimeScale), arg0, arg1)
end


text:
using LeastSquaresOptim, Printf, SparseArrays, Test

# simple factor model
# only problem with "real" optimization
# nice example because J'J is not invertible
# but cholfact in sparse handles this case
function factor()
    name = "factor"
    function f!(fvec, x)
        fvec[1] = 3.0 - x[1] * x[4]
        fvec[2] = 2.0 - x[1] * x[5]
        fvec[3] = 5.0 - x[1] * x[6]

        fvec[4] = 4.5 - x[2] * x[4]
        fvec[5] = 3.2 - x[2] * x[5]
        fvec[6] = 2.0 - x[2] * x[6]

        fvec[7] = 5.0 - x[3] * x[4]
        fvec[8] = 1.3 - x[3] * x[5]
        fvec[9] = 1.5 - x[3] * x[6]
    end

    function g!(J, x)
        fill!(J, 0.0)
        J[1, 1] = -x[4]
        J[1, 4] = -x[1]
        J[2, 1] = -x[5]
        J[2, 5] = -x[1]
        J[3, 1] = -x[6]
        J[3, 6] = -x[1]

        J[4, 2] = -x[4]
        J[4, 4] = -x[2]
        J[5, 2] = -x[5]
        J[5, 5] = -x[2]
        J[6, 2] = -x[6]
        J[6, 6] = -x[2]

        J[7, 3] = -x[4]
        J[7, 4] = -x[3]
        J[8, 3] = -x[5]
        J[8, 5] = -x[3]
        J[9, 3] = -x[6]
        J[9, 6] = -x[3]
    end

    function g!(J::SparseMatrixCSC, x)
        Jvals = nonzeros(J)
        i = 0
        i += 1
        Jvals[i] = -x[4]
        i += 1
        Jvals[i] = -x[1]
        i += 1
        Jvals[i] = -x[5]
        i += 1
        Jvals[i] = -x[1]
        i += 1
        Jvals[i] = -x[6]
        i += 1
        Jvals[i] = -x[1]
        i += 1
        Jvals[i] = -x[4]
        i += 1
        Jvals[i] = -x[2]
        i += 1
        Jvals[i] = -x[5]
        i += 1
        Jvals[i] = -x[2]
        i += 1
        Jvals[i] = -x[6]
        i += 1
        Jvals[i] = -x[2]
        i += 1
        Jvals[i] = -x[4]
        i += 1
        Jvals[i] = -x[3]
        i += 1
        Jvals[i] = -x[5]
        i += 1
        Jvals[i] = -x[3]
        i += 1
        Jvals[i] = -x[6]
        i += 1
        Jvals[i] = -x[3]
    end
    x = ones(6)
    return name, f!, g!, x
end
iter = 0
for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
    for (solver, solver_abbr) in ((LeastSquaresOptim.QR(), :qr), (LeastSquaresOptim.LSMR(), :iter))
        global iter += 1
        global name, f!, g!, x = factor()
        global fcur = ones(9)
        global J = ones(9, 6)
        g!(J, x)
        if solver == LeastSquaresOptim.LSMR()
            J = sparse(J)
        end
        global nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
        global r = optimize!(nls, optimizer(solver))
        if iter == 1
            show(r)
        end
        @printf("%4s %2s %30s %5d %5d   %5d   %10e\n", solver_abbr, optimizer_abbr, name, r.iterations, r.f_calls, r.g_calls, r.ssr)
        @test r.ssr <= 12
        @test r.converged
    end
end

max_stars_repo_name: matthieugomez/LeastSquares

max_issues_repo_name: matthieugomez/LeastSquares

max_forks_repo_name: matthieugomez/LeastSquares


text:
################################################################################
# Mechanisms
################################################################################
include("ant/methods/initialize.jl")
include("atlas/methods/initialize.jl")
include("block/methods/initialize.jl")
include("block2d/methods/initialize.jl")
include("cartpole/methods/initialize.jl")
include("dzhanibekov/methods/initialize.jl")
include("fourbar/methods/initialize.jl")
include("halfcheetah/methods/initialize.jl")
include("hopper/methods/initialize.jl")
include("humanoid/methods/initialize.jl")
include("orbital/methods/initialize.jl")
include("panda/methods/initialize.jl")
include("pendulum/methods/initialize.jl")
include("quadruped/methods/initialize.jl")
include("raiberthopper/methods/initialize.jl")
include("rexhopper/methods/initialize.jl")
include("slider/methods/initialize.jl")
include("sphere/methods/initialize.jl")
include("snake/methods/initialize.jl")
include("tennisracket/methods/initialize.jl")
include("tippetop/methods/initialize.jl")
include("twister/methods/initialize.jl")
include("walker/methods/initialize.jl")
include("tugbot/methods/initialize.jl")

################################################################################
# Environments
################################################################################
include("ant/methods/env.jl")
include("atlas/methods/env.jl")
include("block/methods/env.jl")
include("block2d/methods/env.jl")
include("cartpole/methods/env.jl")
include("halfcheetah/methods/env.jl")
include("hopper/methods/env.jl")
include("panda/methods/env.jl")
include("pendulum/methods/env.jl")
include("quadruped/methods/env.jl")
include("raiberthopper/methods/env.jl")
include("rexhopper/methods/env.jl")
include("walker/methods/env.jl")


text:
"""
    Grey-Wolf Optimizer

    params:
    -------
    * Xb := buffer
    * Xn := n-th solution
    * Xα := α-wolf solution
    * Xβ := β-wolf solution
    * Xδ := δ-wolf solution
"""
function gwo_move!(Xb::VecB, Xn::VecI, Xα::VecI, Xβ::VecI, Xδ::VecI, a::Real) # @code_warntype ✓
    @nrand 3 2 2.0
    @inbounds for i in eachindex(Xb)
        Xni = Xn[i]
        Xαi = Xα[i]
        Xβi = Xβ[i]
        Xδi = Xδ[i]
        tmp = Xαi + Xβi + Xδi
        tmp = tmp - a * (r11 - 1.0) * abs(r12 * Xαi - Xni)
        tmp = tmp - a * (r21 - 1.0) * abs(r22 * Xβi - Xni)
        tmp = tmp - a * (r31 - 1.0) * abs(r32 * Xδi - Xni)
        Xb[i] = tmp / 3.0
    end
end
"""
    Sine-Cosine Optimizer

    params:
    -------
    * Xb := buffer
    * Xn := n-th solution
    * Xr := referred solution
"""
function sco_move!(Xb::VecB, Xn::VecI, Xr::VecI, a::Real) # @code_warntype ✓
    r = 2.0 * rand()
    @inbounds for i in eachindex(Xb)
        Xb[i] = Xn[i] + a * abs(Xn[i] - Xr[i]) * ifelse(rand() < 0.5, sinpi(r), cospi(r))
    end
end
"""
    Arithmetic Search Algorithm Optimizer

    params:
    -------
    * Xb := buffer
    * Xn := n-th solution
    * lb := lower bounds
    * ub := upper bounds
    * p1 := tuning probability, p1 = 0.2 + 0.7 * (it / max. it)
    * p2 := tuning probability, p2 = 1 - exp(log(it / max. it) / 5)
"""
function asa_move!(Xb::VecB, Xn::VecI, lb::NTuple, ub::NTuple, p1::Real, p2::Real) # @code_warntype ✓
    @nrand 3
    if r1 < p1
        # exploitation:  addition / subtractioin
        r3 *= r2 < 0.5 ? 0.5 * p2 : -0.5 * p2
        @simd for i in eachindex(Xb)
            @inbounds Xb[i] = Xn[i] + r3 * (ub[i] - lb[i])
        end
    else
        # exploration:       division     / multiplication
        r3 *= r2 < 0.5 ? 0.5 / (p2 + 1.0) : 0.5 * p2
        @simd for i in eachindex(Xb)
            @inbounds Xb[i] = Xn[i] * r3 * (ub[i] - lb[i])
        end
    end
end
"""
    Water-Cycle Algorithm Optimizer

    params:
    -------
    * Xb    := buffer
    * Xbest := the best solution currently
    * lb    := lower bounds
    * ub    := upper bounds
"""
function wca_move!(Xb::VecB, Xbest::VecI) # @code_warntype ✓
    r = randn()
    @simd for i in eachindex(Xb)
        @inbounds Xb[i] = Xbest[i] + r * 0.31622776601683794 # sqrt(0.1)
    end
end


text:
# WARNING: This should not be used as-is.  It's unfinished experimentation.

using Transformations
using StochasticOptimization
using PenaltyFunctions
import OnlineStats: Mean, Variances, Weight, BoundedEqualWeight

mutable struct Actor{PHI<:Learnable, DIST<:MvNormalTransformation} <: AbstractPolicy
    ϕ::PHI      # map states to dist inputs. ∇logπ is grad(ϕ)
    D::DIST     # the N(ϕ) = N(μ,σ) from which we sample actions
    prep::PreprocessStep
    last_action
    testing::Bool
end
Actor(ϕ,D,prep=NoPreprocessing()) = Actor(ϕ, D, prep, zeros(D.n), false)

# TODO: specialize this on the distribution type
function Reinforce.action(actor::Actor, r, s′, A′)
    z = transform!(actor, s′)
    # @show map(extrema, (s′, z, params(actor), grad(actor)))

    # TODO: make this a general utility method... project_to_actions?
    # project our squashed sample onto into the action space to get our actions
    # a = (â --> [lo,hi])
    a = A′.lo .+ logistic.(z) .* (A′.hi .- A′.lo)
    if !(a in A′)
        warn("a=$a not in A=$(A′)")
        dump(actor)
        error()
        a = rand(A′)
    end
    copy!(actor.last_action, a)
    a
end

@with actor function grad!(actor::Actor)
    grad!(D)
    grad!(ϕ, input_grad(D))
end
@with actor function transform!(actor::Actor, xs::AbstractVector)
    x = if isa(prep, NoPreprocessing)
        xs
    else
        learn!(prep, xs)
        transform!(prep, xs)
    end
    phi = transform!(ϕ, x)
    if testing
        # if we're testing, return the mean exactly
        phi[1:D.n]
    else
        transform!(D, phi)
    end
end
params(actor::Actor) = params(actor.ϕ)
grad(actor::Actor) = grad(actor.ϕ)

#=
This is an implementation of the "AC" algorithm (Actor-critic Algorithm) from:
    Degris et al, "Model-Free Reinforcement Learning with Continuous Action in Practice"
=#

mutable struct OnlineActorCritic{ALGO, T, WGT<:Weight, PEN<:Penalty, ACTOR} <: AbstractPolicy
    δ::T            # last TD δ
    # r̄::T            # estimate of average return
    r̄::Mean{WGT}
    # svar::Variances{WGT}       # feaure whitener TODO: do actual whitening!
    penalty::PEN
    ns::Int         # number of inputs
    nv::Int         # 2ns+na because we do: vcat(s′, s′-s, a)
    na::Int         # the number of actions
    nu::Int         # the number of policy params (2na, since ϕ = vcat(μ,σ))
    x::Function     # a "feature mapping" function x(s)
    v::Vector{T}    # critic params
    eᵛ::Vector{T}   # eligibility trace for updating v (critic params)
    # ϕ::PHI            # u from the paper is really params(ϕ), ∇logπ is grad(ϕ)
    # D::DIST         # the N(ϕ) = N(μ,σ) from which we sample actions
    actor::ACTOR
    w::Vector{T}    # for INAC algos
    eᵘ::Vector{T}   # eligibility trace for updating u (actor params)
    γ::T            # δ decay
    λ::T            # e decay
    # αʳ::T           # learning rate for r̄
    # αᵛ::T           # learning rate for v
    # αᵘ::T           # learning rate for u
    gaᵛ
    gaᵘ
    # gaʷ
    # last_sars′
    xs
end

function OnlineActorCritic(s::AbstractVector, na::Int;
                     T::DataType = eltype(s),
                     algo::Symbol = :AC,
                     wgt_lookback::Int = 20000,
                     prep::PreprocessStep = Whiten(T,2length(s)+na,2length(s)+na,lookback=wgt_lookback),
                     penalty::Penalty = L2Penalty(1e-5),
                     ϕ::Learnable = nnet(2length(s)+na, 2na, [], :relu),
                     D::MvNormalTransformation = MvNormalTransformation(zeros(T,na),zeros(T,na)),
                     x::Function = identity,
                     γ::Number = 1.0,
                     λ::Number = 0.5,
                    #  αʳ::Number = 0.01,
                     αᵛ::Number = 0.01,
                     αᵘ::Number = 0.01,
                     αʷ::Number = 0.01,
                     gaᵛ = OnlineGradAvg(100, lr=αᵛ),
                     gaᵘ = OnlineGradAvg(100, lr=αᵘ)
                    #  gaʷ = OnlineGradAvg(100, lr=αʷ)
                    )
    @assert algo in (:AC, :INAC)
    T = eltype(s)
    ns = length(x(s))
    nv = 2ns+na
    nu = length(params(ϕ))
    wgt = BoundedEqualWeight(wgt_lookback)
    r̄ = Mean(wgt)
    # svar = Variances(nv, wgt)
    actor = Actor(ϕ,D,prep)
    link_nodes!(ϕ, D)

    setup!(gaᵛ, zeros(nv))
    setup!(gaᵘ, ϕ)
    # setup!(gaʷ, ϕ)

    OnlineActorCritic{algo,T,typeof(wgt),typeof(penalty),typeof(actor)}(
        zero(T),
        # zero(T),
        r̄,
        # svar,
        penalty,
        ns,
        nv,
        na,
        nu,
        x,
        zeros(T,nv),
        zeros(T,nv),
        # ϕ,
        # D,
        actor,
        zeros(T,nu),
        zeros(T,nu),
        γ,
        λ,
        # αʳ,
        # αᵛ,
        # αᵘ
        gaᵛ,
        gaᵘ,
        # gaʷ,
        # (s, zeros(na), 0.0, s)
        vcat(s, zeros(T,ns+na))
    )
end

function Reinforce.reset!(ac::OnlineActorCritic{A,T}) where {A,T}
    # reset!(ac.actor)
    fill!(ac.eᵛ, zero(T))
    fill!(ac.eᵘ, zero(T))
end

@with ac function Reinforce.action(ac::OnlineActorCritic, r, s′, A′)
    # ignore s′ and use latest xs (which should include s′)
    for i=1:ns
        xs[i+ns] = s′[i] - xs[i]
        xs[i] = s′[i]
    end
    xs[2ns+1:end] = actor.last_action
    a = action(actor, r, xs, A′)
    # if !isa(actor.prep, NoPreprocessing)
    #     xs[:] = output_value(actor.prep)
    # end
    a
    # # for i=1:ns
    # #     xs[i+ns] = s′[i] - xs[i]
    # #     xs[i] = s′[i]
    # # end
    # # xs[ns+1:2ns] = s′ -
    # # s′ = vcat(s′, s′-last_sars′[1], last_sars′[2])
    # # transform!(ϕ, xs)
    # # z = transform!(D)
    # z = transform!(actor)
    #
    # # TODO: make this a general utility method... project_to_actions?
    # # project our squashed sample onto into the action space to get our actions
    # # a = (â --> [lo,hi])
    # a = A′.lo .+ logistic.(z) .* (A′.hi .- A′.lo)
    # if !(a in A′)
    #     warn("a=$a not in A=$(A′)")
    #     dump(ac)
    #     error()
    #     a = rand(A′)
    # end
    # # xs[2ns+1:end] = a
    # a
end

#= Notes:
    - we can replace αᵛ/αᵘ with a GradientLearner and call update! instead
    - the @with macro (in StochasticOptimization) replaces variables with the
        dot versions if that object has a field of that name
=#

# function whiten(x::AbstractArray, vars::Variances)
#     return x
#     σ = vars.value
#     μ = vars.μ
#     @assert length(x) == length(σ)
#     out = zeros(x)
#     @inbounds for (i,j) in enumerate(eachindex(x))
#         out[j] = (x[j] - μ[i])
#         if σ[i] > 0
#             out[j] /= sqrt(σ[i])
#         end
#     end
#     out
# end

# This is similar to replacing traces for continuous values.
# If we're reversing the trace then we simply add, if we're increasing
# the magnitude then we take the larger of ei/xi.
# This should help with the stability of traces.
# Note: invented by @tbreloff, but I'm sure it exists somewhere already.
function update_eligibilty!(e::AbstractArray{T}, x::AbstractArray{T},
                           γλ::Number; clip::Number = 1e2) where T
    @assert length(e) == length(x)
    @inbounds for i=1:length(e)
        ei = γλ * e[i]
        xi = clamp(x[i], -clip, clip)
        e[i] = if ei < zero(T)
            xi < zero(T) ? min(ei, xi) : ei+xi
        else
            xi > zero(T) ? max(ei, xi) : ei+xi
        end
    end
end

@with ac function learn!(ac::OnlineActorCritic, s, a, r, s′)
    actor.testing && return
    # xs = whiten(x(s), svar)
    # xs′ = whiten(x(s′), svar)
    #
    # # update our input whitener
    # fit!(svar, x(s′))

    # xs = vcat(s, s-last_sars′[1], last_sars′[2])
    xs′ = vcat(s′, s′-s, a)

    prepped_xs = transform!(actor.prep, xs)
    prepped_xs′ = transform!(actor.prep, xs′)

    # compute TD delta
    δ = r - mean(r̄) + γ * dot(v, prepped_xs′) - dot(v, prepped_xs)

    # update average reward
    # r̄ += αʳ * δ
    # r̄ = αʳ * r + (one(αʳ) - αʳ) * r̄
    fit!(r̄, r)

    # update critic
    γλ = γ * λ
    update_eligibilty!(eᵛ, prepped_xs, γλ)
    chg = zeros(v)
    for i=1:nv
        # eᵛ[i] = γλ * eᵛ[i] + xs[i]
        chg[i] = -δ * eᵛ[i] + deriv(penalty, v[i])
        # v[i] += αᵛ * chg #+ (one(αᵛ) - αᵛ) * v[i]
    end
    learn!(v, gaᵛ, chg)

    # compute ∇logπ
    # TODO: add penalty?
    # grad!(D)
    # grad!(ϕ)

    # update actor eligibility trace
    grad!(actor)
    ∇logπ = grad(actor)
    # @show extrema(∇logπ), extrema(eᵘ)
    update_eligibilty!(eᵘ, ∇logπ, γλ)
    # for i=1:nu
    #     eᵘ[i] = γλ * eᵘ[i] + ∇logπ[i]
    # end

    # update the actor (different by algo)
    update_actor!(ac, params(actor), ∇logπ)
    return
end

@with ac function update_actor!(ac::OnlineActorCritic{:AC}, u, ∇logπ)
    # σ² = D.dist.Σ.diag
    chg = zeros(u)
    for i=1:nu
        # note: we multiply by σ² to reduce instabilities
        chg[i] = -δ * eᵘ[i] + deriv(penalty, u[i])
        if isnan(chg[i])
            @show i, δ, eᵘ[i], u[i]
        end
        # chg = δ * eᵘ[i] * σ²[mod1(i,na)] - deriv(penalty, u[i])
        # u[i] += αᵘ * chg
    end
    learn!(u, gaᵘ, chg)
end

# # BROKEN:
# @with ac function update_actor!{T}(ac::OnlineActorCritic{:INAC,T}, u, ∇logπ)
#     ∇ᵀw = dot(∇logπ, w)
#
#     # update w
#     chg = zeros(u)
#     for i=1:nu
#         chg[i] = -δ * eᵘ[i] + ∇logπ[i] * ∇ᵀw
#     end
#     learn!(w, gaʷ, chg)
#
#     # update u
#     for i=1:nu
#         chg[i] = w[i] - deriv(penalty, u[i])
#     end
#     learn!(u, gaᵘ, chg)
# end

# # BROKEN:
# @with ac function update_actor!{T}(ac::OnlineActorCritic{:INAC,T}, u, ∇logπ)
#     ∇ᵀw = dot(∇logπ, w)
#     # @show extrema(∇logπ), ∇ᵀw
#     chg = zeros(u)
#     for i=1:nu
#         w[i] += αᵘ * (δ * eᵘ[i] - ∇logπ[i] * ∇ᵀw) #+ (one(αᵘ) - αᵘ) * w[i]
#         if !isfinite(w[i])
#             @show w δ eᵘ ∇logπ ∇ᵀw
#             @show map(extrema, (w, δ, eᵘ, ∇logπ, ∇ᵀw))
#             error()
#         end
#         chg[i] = w[i] - deriv(penalty, u[i])
#         # u[i] += αᵘ * chg
#     end
#     learn!(u, gaᵘ, chg)
# end

# ----------------------------------------------------------------------
# ----------------------------------------------------------------------

mutable struct EpisodicActorCritic

end


text:
include("lib_for_robot.jl")
function task5(r::Robot)
    robot_path = []
    while (isborder(r, Sud) == false) | (isborder(r, West) == false)
        push!(robot_path, show_path_after_moving(r, West))
        push!(robot_path, show_path_after_moving(r, Sud))
    end
    for side in (Nord, Ost, Sud, West)
        while isborder(r, side) == false
            move!(r, side)
        end
        putmarker!(r)
    end
    back_by_path(r, robot_path)
end


text:
using MPI
using ClimateMachine
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using LinearAlgebra
using Printf
using Dates
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
import ClimateMachine.DGMethods.NumericalFluxes:
    normal_boundary_flux_second_order!

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output = parse(Bool, lowercase(get(ENV, "JULIA_CLIMA_OUTPUT", "false")))

include("advection_diffusion_model.jl")

struct HeatEqn{n, κ, A} <: AdvectionDiffusionProblem end

function init_velocity_diffusion!(
    ::HeatEqn{n},
    aux::Vars,
    geom::LocalGeometry,
) where {n}
    # diffusion of strength 1 in the n direction
    aux.diffusion.D = n * n'
end

# solution is such that
# u(1, t) = 1
# ∇u(0,t) = n
function initial_condition!(
    ::HeatEqn{n, κ, A},
    state,
    aux,
    localgeo,
    t,
) where {n, κ, A}
    ξn = dot(n, localgeo.coord)
    state.ρ = ξn + sum(A .* cos.(κ * ξn) .* exp.(-κ .^ 2 * t))
end
inhomogeneous_data!(::Val{0}, P::HeatEqn, x...) = initial_condition!(P, x...)

function normal_boundary_flux_second_order!(
    ::CentralNumericalFluxSecondOrder,
    bcs,
    ::AdvectionDiffusion{1, dim, HeatEqn{nd, κ, A}},
    fluxᵀn::Vars{S},
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    t,
    _...,
) where {S, dim, nd, κ, A}

    if any_isa(bcs, InhomogeneousBC{0})
        fluxᵀn.ρ = -diff⁻.σ' * n⁻
    elseif any_isa(bcs, InhomogeneousBC{1})
        # Get exact gradient of ρ
        x = aux⁻.coord
        ξn = dot(nd, x)
        ∇ρ = SVector(ntuple(
            i ->
                nd[i] * (1 - sum(A .* κ .* sin.(κ * ξn) .* exp.(-κ .^ 2 * t))),
            Val(3),
        ))

        # Compute flux value
        D = aux⁻.diffusion.D
        fluxᵀn.ρ = -(D * ∇ρ)' * n⁻
    end
end

function test_run(
    mpicomm,
    ArrayType,
    dim,
    topl,
    N,
    timeend,
    FT,
    direction,
    dt,
    n,
    κ = 10 * FT(π) / 2,
    A = 1,
)

    numberofsteps = convert(Int64, cld(timeend, dt))
    dt = timeend / numberofsteps
    @info "time step" dt numberofsteps dt * numberofsteps timeend

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    bcs = (InhomogeneousBC{1}(), InhomogeneousBC{0}())
    model = AdvectionDiffusion{dim}(HeatEqn{n, κ, A}(), bcs; advection = false)
    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
        direction = direction(),
    )

    Q = init_ode_state(dg, FT(0))

    lsrk = LSRK144NiegemannDiehlBusch(dg, Q; dt = dt, t0 = 0)

    eng0 = norm(Q)
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %.16e
                runtime = %s
                norm(Q) = %.16e""",
                gettime(lsrk),
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end
    callbacks = (cbinfo,)

    solve!(
        Q,
        lsrk;
        numberofsteps = numberofsteps,
        callbacks = callbacks,
        adjustfinalstep = false,
    )

    # Print some end of the simulation information
    engf = norm(Q)
    Qe = init_ode_state(dg, FT(timeend))

    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

using Test
let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    base_num_elem = 4

    expected_result = Dict()
    expected_result[2, 1, Float64, EveryDirection] = 0.005157483268127576
    expected_result[2, 2, Float64, EveryDirection] = 6.5687731035717e-5
    expected_result[2, 3, Float64, EveryDirection] = 1.6644861275185443e-6
    expected_result[2, 1, Float64, HorizontalDirection] = 0.020515449798977983
    expected_result[2, 2, Float64, HorizontalDirection] = 0.0005686256942296802
    expected_result[2, 3, Float64, HorizontalDirection] = 1.0132022682547854e-5
    expected_result[2, 1, Float64, VerticalDirection] = 0.02051544979897792
    expected_result[2, 2, Float64, VerticalDirection] = 0.0005686256942296017
    expected_result[2, 3, Float64, VerticalDirection] = 1.0132022682754848e-5
    expected_result[3, 1, Float64, EveryDirection] = 0.001260581018671256
    expected_result[3, 2, Float64, EveryDirection] = 2.214908522198975e-5
    expected_result[3, 3, Float64, EveryDirection] = 5.931735594156876e-7
    expected_result[3, 1, Float64, HorizontalDirection] = 0.005157483268127569
    expected_result[3, 2, Float64, HorizontalDirection] = 6.568773103570526e-5
    expected_result[3, 3, Float64, HorizontalDirection] = 1.6644861273865866e-6
    expected_result[3, 1, Float64, VerticalDirection] = 0.020515449798978087
    expected_result[3, 2, Float64, VerticalDirection] = 0.0005686256942297547
    expected_result[3, 3, Float64, VerticalDirection] = 1.0132022682817856e-5

    expected_result[2, 1, Float32, EveryDirection] = 0.005157135
    expected_result[2, 2, Float32, EveryDirection] = 6.5721644e-5
    expected_result[2, 3, Float32, EveryDirection] = 3.280845e-6
    expected_result[2, 1, Float32, HorizontalDirection] = 0.020514594
    expected_result[2, 2, Float32, HorizontalDirection] = 0.0005684704
    expected_result[2, 3, Float32, HorizontalDirection] = 1.02350195e-5
    expected_result[2, 1, Float32, VerticalDirection] = 0.020514673
    expected_result[2, 2, Float32, VerticalDirection] = 0.0005684843
    expected_result[2, 3, Float32, VerticalDirection] = 1.0227403e-5
    expected_result[3, 1, Float32, EveryDirection] = 0.0012602004
    expected_result[3, 2, Float32, EveryDirection] = 2.2415780e-5
    expected_result[3, 3, Float32, EveryDirection] = 1.1309192e-5
    expected_result[3, 1, Float32, HorizontalDirection] = 0.005157044
    expected_result[3, 2, Float32, HorizontalDirection] = 6.66792e-5
    expected_result[3, 3, Float32, HorizontalDirection] = 9.930429e-5
    expected_result[3, 1, Float32, VerticalDirection] = 0.020514654
    expected_result[3, 2, Float32, VerticalDirection] = 0.0005684157
    expected_result[3, 3, Float32, VerticalDirection] = 3.224683e-5

    @testset "$(@__FILE__)" begin
        for FT in (Float64, Float32)
            numlevels =
                integration_testing ||
                ClimateMachine.Settings.integration_testing ?
                3 : 1
            result = zeros(FT, numlevels)
            for dim in 2:3
                connectivity = dim == 2 ? :face : :full
                for direction in
                    (EveryDirection, HorizontalDirection, VerticalDirection)
                    if direction <: EveryDirection
                        n =
                            dim == 2 ?
                            SVector{3, FT}(1 / sqrt(2), 1 / sqrt(2), 0) :
                            SVector{3, FT}(
                                1 / sqrt(3),
                                1 / sqrt(3),
                                1 / sqrt(3),
                            )
                    elseif direction <: HorizontalDirection
                        n =
                            dim == 2 ? SVector{3, FT}(1, 0, 0) :
                            SVector{3, FT}(1 / sqrt(2), 1 / sqrt(2), 0)
                    elseif direction <: VerticalDirection
                        n =
                            dim == 2 ? SVector{3, FT}(0, 1, 0) :
                            SVector{3, FT}(0, 0, 1)
                    end
                    for l in 1:numlevels
                        Ne = 2^(l - 1) * base_num_elem
                        brickrange = ntuple(
                            j -> range(FT(0); length = Ne + 1, stop = 1),
                            dim,
                        )
                        periodicity = ntuple(j -> false, dim)
                        bc = ntuple(j -> (1, 2), dim)
                        topl = StackedBrickTopology(
                            mpicomm,
                            brickrange;
                            periodicity = periodicity,
                            boundary = bc,
                            connectivity = connectivity,
                        )
                        dt = 1 / (Ne * polynomialorder^2)^2

                        timeend = 0.01

                        @info (ArrayType, FT, dim, direction)
                        result[l] = test_run(
                            mpicomm,
                            ArrayType,
                            dim,
                            topl,
                            polynomialorder,
                            timeend,
                            FT,
                            direction,
                            dt,
                            n,
                        )

                        @test (
                            result[l] ≈
                            FT(expected_result[dim, l, FT, direction]) ||
                            result[l] <
                            FT(expected_result[dim, l, FT, direction])
                        )
                    end
                    @info begin
                        msg = ""
                        for l in 1:(numlevels - 1)
                            rate = log2(result[l]) - log2(result[l + 1])
                            msg *= @sprintf(
                                "\n  rate for level %d = %e\n",
                                l,
                                rate
                            )
                        end
                        msg
                    end
                end
            end
        end
    end
end

nothing


text:
using MixedModels, GLM, DataFrames

data = simulate_data_correlation([[0.3, -0.1], [0.1, -0.4], [0.5, -0.25]])
model = fit(LinearMixedModel, @formula(y ~ Var1 * Var2 + (1|Group)), data)
confint(model)
boot = bootstrap(1000, model)

a = describe(model)

fieldnames(typeof(model.L.data))
fieldnames(typeof(model))

model.formula
model.formula.ex_orig
model.formula.ex
model.formula.lhs
model.formula.rhs

model.trms
model.trms[1]
model.trms[1].refs
model.trms[1].levels

model.trms[2]
model.trms[3]

model.sqrtwts

model.A
model.A.blocks
model.A.block_sizes

model.L
model.L.data.blocks
model.L.data.block_sizes

model.optsum
model.optsum.initial
model.optsum.optimizer
model.optsum.final

using InteractiveUtils
methodswith(typeof(model))

vcov(model)

C = model.trms[3]
ModelMatrix(model)
methods(typeof(model.formula))
varinfo(MixedModels)

model.formula
function model_description(model::MixedModels.LinearMixedModel)

    modelname = "linear mixed model"
    outcome = model.mf.terms.eterms[1]
    predictors = model.mf.terms.eterms[2:end]
    formula = model.formula

    output = Dict(
        "Model" => modelname,
        "Outcome" => outcome,
        "Predictors" => predictors,
        "Formula" => formula,
        "text_description" => "We fitted a linear regression to predict $outcome with $(join(predictors, ", ", " and ")) ($formula)."
        )
    return Report(text=output["text_description"], values=output)
end


text:
"""
    Mlmnet_cvDeprecated(MLMNets, lambdas, alphas, data, rowFolds, colFolds)

Type for storing the results of running cross-validation for `mlmnet`

"""
mutable struct Mlmnet_cv

    # MlmnetDeprecated objects
    MLMNets::Array{Mlmnet,1}
    # Lambda penalties
    lambdas::Array{Float64,1}
    alphas::Array{Float64,1}

    # Response and predictor matrices
    data::RawData
    # Row folds
    rowFolds::Array{Array,1}
    # Column folds
    colFolds::Array{Array,1}

    # These will be generated and should not be supplied as input to the
    # constructor.
    # Test MSE for each of the CV folds for each lambdaL1 and lambdaL2
    mse::Array{Float64,3}
    # Proportion of zero interaction coefficients for each of the CV folds
    # for each lambdaL1 and lambdaL2
    propZero::Array{Float64,3}

    Mlmnet_cv(MLMNets, lambdas, alphas, data, rowFolds, colFolds, dig) =
        new(MLMNets, lambdas, alphas, data, rowFolds, colFolds,
            calc_mse(MLMNets, data, lambdas, alphas, rowFolds, colFolds),
            calc_prop_zero(MLMNets, lambdas, alphas; dig=dig))
end

"""
    mlmnet_cv(data, lambdas, alphas, rowFolds, colFolds;
              method="ista", isNaive=false,
              hasXIntercept=true, hasZIntercept=true,
              toXReg=trues(size(get_X(data), 2)),
              toZReg=trues(size(get_Z(data), 2)),
              toXInterceptReg=false, toZInterceptReg=false,
              toNormalize=true, isVerbose=true,
              stepsize=0.01, setStepsize=true,
              dig=12, funArgs...)

Performs cross-validation for `mlmnet` using row and column folds from user
input.

# Arguments

- data = RawData object
- lambdas = 1d array of floats consisting of the total penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2
- rowFolds = 1d array of arrays (one array for each fold), each containing
  the indices for a row fold; must be same length as colFolds. Can be
  generated with a call to `make_folds`, which is based on `Kfold` from the
  MLBase package.
- colFolds = 1d array of arrays (one array for each fold), each containing
  the indices for a column fold; must be same length as rowFolds. Can be
  generated with a call to `make_folds`, which is based on `Kfold` from the
  MLBase package

# Keyword arguments

- methods = function name that applies the Elastic-net penalty estimate method;
  default is `ista`, and the other methods are `fista`, `fista_bt`, `admm` and `cd`
- isNaive = boolean flag indicating whether to solve the Naive or non-Naive
  Elastic-net problem
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- toXReg = 1d array of bit flags indicating whether or not to regularize each
  of the `X` (row) effects. Defaults to 2d array of `true`s with length
  equal to the number of `X` effects (equivalent to `data.p`).
- toZReg = 1d array of bit flags indicating whether or not to regularize each
  of the `Z` (column) effects. Defaults to 2d array of `true`s with length
  equal to the number of `Z` effects (equivalent to `data.q`).
- toXInterceptReg = boolean flag indicating whether or not to regularize the
  `X` intercept Defaults to `false`.
- toZInterceptReg = boolean flag indicating whether or not to regularize the
  `Z` intercept. Defaults to `false`.
- toNormalize = boolean flag indicating if the columns of `X` and `Z`
  should be standardized (to mean 0, standard deviation 1). Defaults to `true`.
- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for coordinate
  descent and when `setStepsize` is set to `true` for `ista!` and `fista!`).
  Defaults to `0.01`.
- setStepsize = boolean flag indicating whether the fixed step size should be
  calculated (for `ista!` and `fista!`). Defaults to `true`.
- dig = integer; digits of precision for zero coefficients. Defaults to 12.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

An Mlmnet_cvDeprecated object.

# Some notes

This is the base `mlmnet_cv` function that all other variants call. Folds
are computed in parallel when possible.

"""
function mlmnet_cv(data::RawData,
                   lambdas::AbstractArray{Float64,1},
                   alphas::AbstractArray{Float64,1},
                   rowFolds::Array{Array{Int64,1},1},
                   colFolds::Array{Array{Int64,1},1};
                   method::String="ista", isNaive::Bool=false,
                   hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                   toXReg::BitArray{1}=trues(size(get_X(data), 2)),
                   toZReg::BitArray{1}=trues(size(get_Z(data), 2)),
                   toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                   toNormalize::Bool=true, isVerbose::Bool=true,
                   stepsize::Float64=0.01, setStepsize::Bool=true,
                   dig::Int64=12, funArgs...)
                   # tuningMode::String = "l1, l2"

                   # if(tuningMode == "alpha, lambda")
                   # lambdasL1 = lambdas.*alpha
                   # lambdasL2 = lambdas.*(1.-alpha)

    # Number of folds
    nFolds = length(rowFolds)
    # Ensure that the number of row folds is the same as the number of column
    # folds
    if nFolds != length(colFolds)
        error("Number of row folds must equal number of column folds. ")
    else
        println_verbose(string("Performing ", nFolds,
                               "-fold cross validation."), isVerbose)
    end

    # Generate list of RawData objects that hold out each of the folds
    dataFolds = Array{RawData}(undef, nFolds)
    for i in 1:nFolds
        dataFolds[i] = RawData(Response(get_Y(data)[rowFolds[i],
                                                    colFolds[i]]),
                               Predictors(get_X(data)[rowFolds[i],:],
                                          get_Z(data)[colFolds[i],:],
                                          data.predictors.hasXIntercept,
                                          data.predictors.hasZIntercept))
    end

    # Run mlmnet on each RawData object, in parallel when possible

    MLMNets = Distributed.pmap(data -> mlmnet(data, lambdas, alphas;
                                              method=method, isNaive=isNaive,
                                              hasXIntercept=hasXIntercept,
                                              hasZIntercept=hasZIntercept,
                                              toXReg=toXReg, toZReg=toZReg,
                                              toXInterceptReg=toXInterceptReg,
                                              toZInterceptReg=toZInterceptReg,
                                              toNormalize=toNormalize,
                                              isVerbose=isVerbose,
                                              stepsize=stepsize,
                                              setStepsize=setStepsize,
                                              funArgs...), dataFolds)

    return Mlmnet_cv(MLMNets, lambdas, alphas, data, rowFolds, colFolds, dig)
end

"""
    mlmnet_cv(data, lambdas, alphas, rowFolds, nColFolds;
    method="ista", isNaive=false,
    hasXIntercept=true, hasZIntercept=true,
    toXReg=trues(size(get_X(data), 2)),
    toZReg=trues(size(get_Z(data), 2)),
    toXInterceptReg=false, toZInterceptReg=false,
    toNormalize=true, isVerbose=true,
    stepsize=0.01, setStepsize=true,
    dig=12, funArgs...)

Performs cross-validation for `mlmnet` using row folds from user input and
non-overlapping column folds randomly generated using a call to `make_folds`.
Calls the base `mlmnet_cv` function.

# Arguments

- data = RawData object
- lambdas = 1d array of floats consisting of lambda penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2
- rowFolds = 1d array of arrays (one array for each fold), each containing
  the indices for a row fold. Can be generated with a call to `make_folds`,
  which is based on `Kfold` from the MLBase package.
- nColFolds = integer corresponding to the number of column folds to randomly
  generate. Can be either `length(rowFolds)` or `1` (to use all columns in
  every fold).

# Keyword arguments

- methods = function name that applies the Elastic-net penalty estimate method;
  default is `ista`, and the other methods are `fista`, `fista_bt`, `admm` and `cd`
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- toXReg = 1d array of bit flags indicating whether or not to regularize each
  of the `X` (row) effects. Defaults to 2d array of `true`s with length
  equal to the number of `X` effects (equivalent to `data.p`).
- toZReg = 1d array of bit flags indicating whether or not to regularize each
  of the `Z` (column) effects. Defaults to 2d array of `true`s with length
  equal to the number of `Z` effects (equivalent to `data.q`).
- toXInterceptReg = boolean flag indicating whether or not to regularize the
  `X` intercept Defaults to `false`.
- toZInterceptReg = boolean flag indicating whether or not to regularize the
  `Z` intercept. Defaults to `false`.
- toNormalize = boolean flag indicating if the columns of `X` and `Z`
  should be standardized (to mean 0, standard deviation 1). Defaults to `true`.
- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for coordinate
  descent and when `setStepsize` is set to `true` for `ista!` and `fista!`).
  Defaults to `0.01`.
- setStepsize = boolean flag indicating whether the fixed step size should be
  calculated (for `ista!` and `fista!`). Defaults to `true`.
- dig = integer; digits of precision for zero coefficients. Defaults to 12.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

An Mlmnet_cvDeprecated object.

# Some notes

Folds are computed in parallel when possible.

"""

function mlmnet_cv(data::RawData,
                   lambdas::AbstractArray{Float64,1},
                   alphas::AbstractArray{Float64,1},
                   rowFolds::Array{Array{Int64,1},1}, nColFolds::Int64;
                   method::String="ista", isNaive::Bool=false,
                   hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                   toXReg::BitArray{1}=trues(size(get_X(data), 2)),
                   toZReg::BitArray{1}=trues(size(get_Z(data), 2)),
                   toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                   toNormalize::Bool=true, isVerbose::Bool=true,
                   stepsize::Float64=0.01, setStepsize::Bool=true,
                   dig::Int64=12, funArgs...)

    # Generate random column folds
    colFolds = make_folds(data.m, nColFolds, length(rowFolds))

    # Pass in user input row folds and randomly generated column folds to the
    # base mlmnet_cv function
    mlmnet_cv(data, lambdas, alphas, rowFolds, colFolds;
              method=method, isNaive=isNaive,
              hasXIntercept=hasXIntercept, hasZIntercept=hasZIntercept,
              toXReg=toXReg, toZReg=toZReg,
              toXInterceptReg=toXInterceptReg,
              toZInterceptReg=toZInterceptReg,
              isVerbose=isVerbose, toNormalize=toNormalize,
              stepsize=stepsize, setStepsize=setStepsize, dig=dig, funArgs...)
end

"""
    mlmnet_cv(data, lambdas, alphas, nRowFolds, colFolds;
              method="ista", isNaive=false,
              hasXIntercept=true, hasZIntercept=true,
              toXReg=trues(size(get_X(data), 2)),
              toZReg=trues(size(get_Z(data), 2)),
              toXInterceptReg=false, toZInterceptReg=false,
              toNormalize=true, isVerbose=true,
              stepsize=0.01, setStepsize=true,
              dig=12, funArgs...)

Performs cross-validation for `mlmnet` using non-overlapping row folds
randomly generated using a call to `make_folds` and column folds from user
input. Calls the base `mlmnet_cv` function.

# Arguments

- data = RawData object
- lambdas = 1d array of floats consisting of lambda penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2
- nRowFolds = integer corresponding to the number of row folds to randomly
  generate. Can be either `length(colFolds)` or `1`
  (to use all rows in every fold).
- colFolds = 1d array of arrays (one array for each fold), each containing
  the indices for a column fold. Can be generated with a call to `make_folds`,
  which is based on `Kfold` from the MLBase package.

# Keyword arguments

- methods = function name that applies the Elastic-net penalty estimate method;
  default is `ista`, and the other methods are `fista`, `fista_bt`, `admm` and `cd`
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- toXReg = 1d array of bit flags indicating whether or not to regularize each
  of the `X` (row) effects. Defaults to 2d array of `true`s with length
  equal to the number of `X` effects (equivalent to `data.p`).
- toZReg = 1d array of bit flags indicating whether or not to regularize each
  of the `Z` (column) effects. Defaults to 2d array of `true`s with length
  equal to the number of `Z` effects (equivalent to `data.q`).
- toXInterceptReg = boolean flag indicating whether or not to regularize the
  `X` intercept Defaults to `false`.
- toZInterceptReg = boolean flag indicating whether or not to regularize the
  `Z` intercept. Defaults to `false`.
- toNormalize = boolean flag indicating if the columns of `X` and `Z`
  should be standardized (to mean 0, standard deviation 1). Defaults to `true`.
- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for coordinate
  descent and when `setStepsize` is set to `true` for `ista!` and `fista!`).
  Defaults to `0.01`.
- setStepsize = boolean flag indicating whether the fixed step size should be
  calculated (for `ista!` and `fista!`). Defaults to `true`.
- dig = integer; digits of precision for zero coefficients. Defaults to 12.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

An Mlmnet_cvDeprecated object.

# Some notes

Folds are computed in parallel when possible.

"""

function mlmnet_cv(data::RawData,
                   lambdas::AbstractArray{Float64,1},
                   alphas::AbstractArray{Float64,1},
                   nRowFolds::Int64, colFolds::Array{Array{Int64,1},1};
                   method::String="ista", isNaive::Bool=false,
                   hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                   toXReg::BitArray{1}=trues(size(get_X(data), 2)),
                   toZReg::BitArray{1}=trues(size(get_Z(data), 2)),
                   toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                   toNormalize::Bool=true, isVerbose::Bool=true,
                   stepsize::Float64=0.01, setStepsize::Bool=true,
                   dig::Int64=12, funArgs...)

    # Generate random row folds
    rowFolds = make_folds(data.n, nRowFolds, length(colFolds))

    # Pass in randomly generated row folds and user input column folds to the
    # base mlmnet_cv function
    mlmnet_cv(data, lambdas, alphas, rowFolds, colFolds;
              method=method, isNaive=isNaive,
              hasXIntercept=hasXIntercept, hasZIntercept=hasZIntercept,
              toXReg=toXReg, toZReg=toZReg,
              toXInterceptReg=toXInterceptReg,
              toZInterceptReg=toZInterceptReg,
              isVerbose=isVerbose, toNormalize=toNormalize,
              stepsize=stepsize, setStepsize=setStepsize, dig=dig, funArgs...)
end

"""
    mlmnet_cv(data, lambdas, alphas, nRowFolds=10, nColFolds=10;
    method="ista", isNaive=false,
    hasXIntercept=true, hasZIntercept=true,
    toXReg=trues(size(get_X(data), 2)),
    toZReg=trues(size(get_Z(data), 2)),
    toXInterceptReg=false, toZInterceptReg=false,
    toNormalize=true, isVerbose=true,
    stepsize=0.01, setStepsize=true,
    dig=12, funArgs...)

Performs cross-validation for `mlmnet` using non-overlapping row and column
folds randomly generated using calls to `make_folds`. Calls the base
`mlmnet_cv` function.

# Arguments

- data = RawData object
- lambdas = 1d array of floats consisting of lambda penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2
- nRowFolds = integer corresponding to the number of row folds to randomly
  generate. Can be either equal to `nColFolds` or `1`
  (to use all rows in every fold).
- nColFolds = integer corresponding to the number of column folds to randomly
  generate. Can be either equal to `nRowFolds` or `1` (to use all columns in
  every fold).

# Keyword arguments

- methods = function name that applies the Elastic-net penalty estimate method;
  default is `ista`, and the other methods are `fista`, `fista_bt`, `admm` and `cd`
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- toXReg = 1d array of bit flags indicating whether or not to regularize each
  of the `X` (row) effects. Defaults to 2d array of `true`s with length
  equal to the number of `X` effects (equivalent to `data.p`).
- toZReg = 1d array of bit flags indicating whether or not to regularize each
  of the `Z` (column) effects. Defaults to 2d array of `true`s with length
  equal to the number of `Z` effects (equivalent to `data.q`).
- toXInterceptReg = boolean flag indicating whether or not to regularize the
  `X` intercept Defaults to `false`.
- toZInterceptReg = boolean flag indicating whether or not to regularize the
  `Z` intercept. Defaults to `false`.
- toNormalize = boolean flag indicating if the columns of `X` and `Z`
  should be standardized (to mean 0, standard deviation 1). Defaults to `true`.
- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for coordinate
  descent and when `setStepsize` is set to `true` for `ista!` and `fista!`).
  Defaults to `0.01`.
- setStepsize = boolean flag indicating whether the fixed step size should be
  calculated (for `ista!` and `fista!`). Defaults to `true`.
- dig = integer; digits of precision for zero coefficients. Defaults to 12.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

An Mlmnet_cvDeprecated object.

# Some notes

Folds are computed in parallel when possible.

The default method for choosing the fixed step size for `fista!` or `ista!`
is to use the reciprocal of the product of the maximum eigenvalues of
`X*transpose(X)` and `Z*transpose(Z)`. This is computed when `fista!` or
`ista!` is passed into the `fun` argument and `setStepsize` is set to `true`.
If `setStepsize` is set to `false`, the value of the `stepsize` argument will
be used as the fixed step size. Note that obtaining the eigenvalues when `X`
and/or `Z` are very large may exceed computational limitations.

Specifying a good starting step size (`stepsize`) and multiplying factor
(`gamma`) when `fista_bt!` is passed into the `fun` argument can be difficult.
Shrinking the step size too gradually can result in slow convergence. Doing so
too quickly can cause the criterion to diverge. We have found that setting
`stepsize` to 0.01 often works well in practice; choice of `gamma` appears to
be less consequential.

"""

function mlmnet_cv(data::RawData,
                   lambdas::Array{Float64,1}, alphas::Array{Float64,1},
                   nRowFolds::Int64=10, nColFolds::Int64=10;
                   method::String="ista",
                   isNaive::Bool=false,
                   hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                   toXReg::BitArray{1}=trues(size(get_X(data), 2)),
                   toZReg::BitArray{1}=trues(size(get_Z(data), 2)),
                   toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                   toNormalize::Bool=true, isVerbose::Bool=true,
                   stepsize::Float64=0.01, setStepsize::Bool=true,
                   dig::Int64=12, funArgs...)

    # Generate random row and column folds
    rowFolds = make_folds(data.n, nRowFolds, max(nRowFolds, nColFolds))
    colFolds = make_folds(data.m, nColFolds, max(nRowFolds, nColFolds))

    # Pass in randomly generated row and column folds to the base mlmnet_cv
    # function
    mlmnet_cv(data, lambdas, alphas, rowFolds, colFolds;
              method=method, isNaive=isNaive,
              hasXIntercept=hasXIntercept, hasZIntercept=hasZIntercept,
              toXReg=toXReg, toZReg=toZReg,
              toXInterceptReg=toXInterceptReg,
              toZInterceptReg=toZInterceptReg,
              isVerbose=isVerbose, toNormalize=toNormalize,
              stepsize=stepsize, setStepsize=setStepsize, dig=dig, funArgs...)
end

"""
    mlmnet_cv(data, lambdas, nRowFolds=10, nColFolds=10;
              method="ista", isNaive=false,
              hasXIntercept=true, hasZIntercept=true,
              toXReg=trues(size(get_X(data), 2)),
              toZReg=trues(size(get_Z(data), 2)),
              toXInterceptReg=false, toZInterceptReg=false,
              toNormalize=true, isVerbose=true,
              stepsize=0.01, setStepsize=true,
              dig=12, funArgs...)

Performs cross-validation for `mlmnet` using non-overlapping row and column
folds randomly generated using calls to `make_folds`. Calls the base
`mlmnet_cv` function.

"""

function mlmnet_cv(data::RawData,
                      lambdas::Array{Float64,1},
                      nRowFolds::Int64=10, nColFolds::Int64=10;
                      method::String="ista",
                      isNaive::Bool=false,
                      hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                      toXReg::BitArray{1}=trues(size(get_X(data), 2)),
                      toZReg::BitArray{1}=trues(size(get_Z(data), 2)),
                      toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                      toNormalize::Bool=true, isVerbose::Bool=true,
                      stepsize::Float64=0.01, setStepsize::Bool=true,
                      dig::Int64=12, funArgs...)

    alphas = [1.0]

    # Pass in randomly generated row and column folds to the base mlmnet_cv
    # function
    mlmnet_cv(data, lambdas, alphas, nRowFolds, nColFolds;
                  method=method, isNaive=isNaive,
                  hasXIntercept=hasXIntercept, hasZIntercept=hasZIntercept,
                  toXReg=toXReg, toZReg=toZReg,
                  toXInterceptReg=toXInterceptReg,
                  toZInterceptReg=toZInterceptReg,
                  isVerbose=isVerbose, toNormalize=toNormalize,
                  stepsize=stepsize, setStepsize=setStepsize, dig=dig, funArgs...)

end