text:
using SeisPlot, SeisAcoustic

# ==============================================================================
#             read physical model
function lsrtm_BP()

  dir_work = joinpath(homedir(), "Desktop/lsrtm");

  path_rho = joinpath(dir_work, "physical_model/rho.rsf");
  path_vel = joinpath(dir_work, "physical_model/vel.rsf");

  # read the physical model
  (hdr_rho, rho) = read_RSdata(path_rho);
  (hdr_vel, vel) = read_RSdata(path_vel);

  # # cropped model for imaging
  vel = vel[1:2:1350,650:2:3700];
  rho = rho[1:2:1350,650:2:3700];
  vel = model_smooth(vel, 15);

  rho1 = copy(rho); rho1 .= minimum(rho);
  vel1 = copy(vel); vel1 .= vel[1];

  SeisPlotTX(vel, dx=0.0125, dy=0.0125, yticks=0:2:8, xticks=0:2:19, hbox=1.3*3, wbox=3.7*3, xlabel="X (km)", ylabel="Z (km)", cmap="rainbow", vmin=minimum(vel), vmax=maximum(vel));
  cbar = colorbar(); cbar.set_label("km/s"); tight_layout(); savefig("/Users/wenlei/Desktop/vel.pdf"); close();

  SeisPlotTX(rho, dx=0.0125, dy=0.0125, yticks=0:2:8, xticks=0:2:19, hbox=1.3*3, wbox=3.7*3, xlabel="X (km)", ylabel="Z (km)", cmap="rainbow", vmin=minimum(rho), vmax=maximum(rho));
  cbar = colorbar(); cbar.set_label("g/cm^3"); tight_layout(); savefig("/Users/wenlei/Desktop/rho.pdf"); close();

  # vertical and horizontal grid size
  dz = 6.25; dx = 6.25;

  # time step size and maximum modelling length
  dt = 0.0007; tmax = 6.0;

  # top boundary condition
  free_surface = false;
  data_format  = Float32;
  order        = 5;

  # tdparams for generating observations
  fidiff_hete = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
                    data_format=data_format, order=order);

  # tdparams for removing direct arrival
  fidiff_homo = TdParams(rho1, vel1, free_surface, dz, dx, dt, tmax;
                         data_format=data_format, order=order);

  # tdparams for imaging
  fidiff = TdParams(rho1, vel, free_surface, dz, dx, dt, tmax;
                    data_format=data_format, order=order);

  # initialize a source
  # isz = 2; isx = fidiff.nx;
  # src = Source(isz, isx, fidiff; amp=100000, fdom=20, type_flag="miniphase");

  # # vector of source
  isx = collect(40 : 12 : fidiff.nx-50); isz = 2*ones(length(isx));
  src = get_multi_sources(isz, isx, fidiff; amp=100000, fdom=20, type_flag="miniphase");

  # generate observed data
  irx = collect(1: 2 : fidiff.nx);
  irz = 2 * ones(Int64, length(irx));

  dir_obs        = joinpath(dir_work, "observations");
  dir_sourceside = joinpath(dir_work, "sourceside");

  # prepare observations and sourceside wavefield
  get_reflections(dir_obs, irz, irx, src, fidiff_hete, fidiff_homo);
  get_wavefield_bound(dir_sourceside, src, fidiff);

  born_params = (irz=irz, irx=irx, dir_sourceside=dir_sourceside, fidiff=fidiff, normalization_flag=true, mute_index=10);
  (x, his) = cgls(born_approximation, dir_obs; dir_work=dir_work, op_params=born_params,
                  d_axpby=recordings_axpby!, m_axpby=image_axpby!, d_norm=recordings_norm, m_norm=l2norm_rsf);
end

lsrtm_BP();

# download the data
scp wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/Desktop/lsrtm/iterations/iteration_1.rsf /Users/wenlei/Desktop/
scp wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/Desktop/lsrtm/sourceside/normalization.rsf /Users/wenlei/Desktop/

(hdr, m) = read_RSdata("/Users/wenlei/Desktop/iteration_1.rsf")
(hdr, s) = read_RSdata("/Users/wenlei/Desktop/normalization.rsf")
s = reshape(s, 675, 1526);
m = reshape(m, 675, 1526);
m .= m .* s;
m1= laplace_filter(m);

p2 = copy(rho)
for i2 = 1 : size(p2,2)
    for i1 = 1 : size(p2, 1)-1
        p2[i1,i2] = (p2[i1+1,i2]-p2[i1,i2]) / (p2[i1+1,i2]+p2[i1,i2])
    end
end

SeisPlotTX(m, cmap="gray", hbox=6.75, wbox=15.2);
SeisPlotTX(s, cmap="rainbow", hbox=6.75, wbox=15.2, vmax=maximum(s), vmin=minimum(s));
SeisPlotTX(m1, cmap="gray", hbox=6.75, wbox=15.2, pclip=96);
SeisPlotTX(p2, cmap="gray", hbox=6.75, wbox=15.2);
SeisPlotTX(vel, cmap="rainbow", hbox=6.75, wbox=15.2, vmax=maximum(vel), vmin=minimum(vel));

# # ==============================================================================
# #             convert physical model of SU to RSF
# work_dir = joinpath(homedir(), "Desktop/tmp_LSRTM");
#
# path_rho = joinpath(work_dir, "physical_model/rho.su");
# path_vel = joinpath(work_dir, "physical_model/vel.su");
#
# num_samples = 1911;
# num_traces  = 5395;
# trace = zeros(Float32, num_samples);
# vel   = zeros(Float32, num_samples, num_traces);
# rho   = zeros(Float32, num_samples, num_traces);
#
# fid_vel = open(path_vel, "r");
# fid_rho = open(path_rho, "r");
# idx = [1];
# for i = 1 : num_traces
#     skip(fid_vel, 240)
#     skip(fid_rho, 240)
#
#     read!(fid_vel, trace)
#     copyto!(vel, idx[1], trace, 1, num_samples)
#
#     read!(fid_rho, trace)
#     copyto!(rho, idx[1], trace, 1, num_samples)
#     idx[1] = idx[1] + num_samples
# end
# close(fid_vel);
# close(fid_rho);
#
# # save the physical model
# hdr_rho = RegularSampleHeader(rho; o1=0., d1=6.25, label1="Z", unit1="m",
#                                    o2=0., d2=6.25, label2="X", unit2="m",
#                                    title="density model");
#
# hdr_vel = RegularSampleHeader(vel; o1=0., d1=6.25, label1="Z", unit1="m",
#                                    o2=0., d2=6.25, label2="X", unit2="m",
#                                    title="velocity model");
#
# path_rho = joinpath(work_dir, "physical_model/rho.rsf");
# path_vel = joinpath(work_dir, "physical_model/vel.rsf");
# write_RSdata(path_rho, hdr_rho, rho);
# write_RSdata(path_vel, hdr_vel, vel);


text:
module UnstructuredDiscreteModelsTests

using Test
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Geometry: DiscreteModelMock
using Gridap.Io

m = DiscreteModelMock()
g = get_grid(m)

model = UnstructuredDiscreteModel(g)
test_discrete_model(model)

domain = (0,1,0,1)
partition = (2,3)
m = CartesianDiscreteModel(domain,partition)
model = UnstructuredDiscreteModel(m)
test_discrete_model(model)
@test model === UnstructuredDiscreteModel(model)

m = DiscreteModelMock()
model = UnstructuredDiscreteModel(m)
test_discrete_model(model)
@test model === UnstructuredDiscreteModel(model)

m = DiscreteModelMock()
model = UnstructuredDiscreteModel(m)

dict = to_dict(model)
model2 = from_dict(UnstructuredDiscreteModel,dict)
test_discrete_model(model2)

model2 = from_json(UnstructuredDiscreteModel,to_json(model))
test_discrete_model(model2)

d = mktempdir()
f = joinpath(d,"model.jld2")

to_jld2_file(model,f)
model2 == from_jld2_file(typeof(model),f)

@test model.face_labeling.d_to_dface_to_entity == model2.face_labeling.d_to_dface_to_entity
@test model.face_labeling.tag_to_entities == model2.face_labeling.tag_to_entities
@test model.face_labeling.tag_to_name == model2.face_labeling.tag_to_name
@test model.grid.cell_nodes == model2.grid.cell_nodes
@test model.grid.cell_types == model2.grid.cell_types
@test model.grid.node_coordinates == model2.grid.node_coordinates
@test model.grid.reffes == model2.grid.reffes
@test model.grid_topology.cell_type == model2.grid_topology.cell_type
@test model.grid_topology.polytopes == model2.grid_topology.polytopes
@test model.grid_topology.vertex_coordinates == model2.grid_topology.vertex_coordinates

rm(d,recursive=true)

end # module

ext: jl


text:
"""
`module  JAC.MultiPhotonDeExcitation`
    ... a submodel of JAC that contains all methods for computing multi-photon excitation and decay rates.
"""
module MultiPhotonDeExcitation

    using Printf, QuadGK, ..AngularMomentum, ..AtomicState, ..Basics, ..Defaults, ..ManyElectron, ..Nuclear,
                          ..PhotoEmission, ..Radial, ..TableStrings

    """
    `abstract type MultiPhotonDeExcitation.AbstractMultiPhotonProperty`
        ... defines an abstract type to distinguish different multi-photon properties; however, by far not all of these properties
            are usually relevant for all the process::MultiPhotonDeExcitation.AbstractMultiPhotonProcess; cf. the JAC User Guide:

        + struct EnergyDiffCs
                    ... Energy differential cross section for initially unpolarized atoms (TwoPhotonEmission, ...)
        + struct TotalCsLinear
                    ... Total cross section for linerarly-polarized incident radiation (TwoPhotonAbsorptionMonochromatic, ...)
        + struct TotalCsRightCircular
                    ... Total cross section for right-circularly polarized incident radiation (TwoPhotonAbsorptionMonochromatic, ...)
        + struct TotalCsUnpolarized
                    ... Total cross section for unpolarized incident radiation (TwoPhotonAbsorptionMonochromatic, ...)
        + struct TotalCsDensityMatrix
                    ... Total cross section for incident radiation with user-defined density matrix (not yet).
                        (TwoPhotonAbsorptionMonochromatic, ...)
    """
    abstract type  AbstractMultiPhotonProperty       end

    """
    `struct  MultiPhotonDeExcitation.EnergyDiffCs  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty`
        ... a struct to request for the calculation of energy-diffrential cross sections in a two-photon emission process.
    """
    struct   EnergyDiffCs  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty      end

    function Base.show(property::EnergyDiffCs)
        sa = "energy-differential cross sections for multi-photon emission";  print(io, sa, "\n")
    end

    """
    `struct  MultiPhotonDeExcitation.TotalCsLinear  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty`
        ... a struct to request for the calculation of the total (absorption) cross sections for linearly-polarized
            incident radiation.
    """
    struct   TotalCsLinear  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty      end

    function Base.show(property::TotalCsLinear)
        sa = "total (absorption) cross sections for linearly-polarized radiation";  print(io, sa, "\n")
    end

    """
    `struct  MultiPhotonDeExcitation.TotalCsRightCircular  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty`
        ... a struct to request for the calculation of the total (absorption) cross sections for right-circularly polarized
            incident radiation.
    """
    struct   TotalCsRightCircular  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty      end

    function Base.show(property::TotalCsRightCircular)
        sa = "total (absorption) cross sections for right-circularly polarized radiation";  print(io, sa, "\n")
    end

    """
    `struct  MultiPhotonDeExcitation.TotalCsUnpolarized  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty`
        ... a struct to request for the calculation of the total (absorption) cross sections for unpolarized
            incident radiation.
    """
    struct   TotalCsUnpolarized  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProperty      end

    function Base.show(property::TotalCsUnpolarized)
        sa = "total (absorption) cross sections for unpolarized radiation";  print(io, sa, "\n")
    end

    """
    `abstract type MultiPhotonDeExcitation.AbstractMultiPhotonProcess`
        ... defines an abstract type to distinguish different multi-photon processes; see also:

        + struct TwoPhotonAbsorptionMonochromatic
                    ... to specify a two-photon absorption process with monochromatic and equally-polarized photons,
                        usually from the same beam.
        + struct TwoPhotonAbsorptionBichromatic
                    ... to specify a two-photon absorption process by photons with well-defined frequency, propagation
                        direction and helicity, usually from different (oriented) beams.
        + struct TwoPhotonEmission
                    ... to specify a two-photon emission process, often for initially-unpolarized atoms.
    """
    abstract type  AbstractMultiPhotonProcess       end

    """
    `struct  MultiPhotonDeExcitation.TwoPhotonAbsorptionMonochromatic  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProcess`
        ... a struct to specify a two-photon absorption process with monochromatic and equally-polarized photons,
            usually from the same beam.

        + properties        ::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1}
            ... List of the multi-photon processes that are to be calculated for this absorption process.
    """
    struct   TwoPhotonAbsorptionMonochromatic  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProcess
        properties          ::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1}
    end

    """
    `MultiPhotonDeExcitation.TwoPhotonAbsorptionMonochromatic()`
        ... constructor for an 'default' instance of a MultiPhotonDeExcitation.TwoPhotonAbsorptionMonochromatic() process.
    """
    function TwoPhotonAbsorptionMonochromatic()
        TwoPhotonAbsorptionMonochromatic( [TotalCsLinear(), TotalCsUnpolarized()] )
    end

    function Base.string(process::MultiPhotonDeExcitation.TwoPhotonAbsorptionMonochromatic)
        sa = "Two-photon absorption process with mono-chromatic photons (from the same beam):"
        return( sa )
    end

    function Base.show(io::IO, process::MultiPhotonDeExcitation.TwoPhotonAbsorptionMonochromatic)
        sa = Base.string(process);               print(io, sa, "\n")
        println(io, "properties:                 $(process.properties)  ")
    end

    """
    `struct  MultiPhotonDeExcitation.TwoPhotonAbsorptionBichromatic  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProcess`
        ... a struct to specify a two-photon absorption process by photons with well-defined frequency, propagation direction
            and helicity, usually from different (oriented) beams.

        + properties        ::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1}
            ... List of the multi-photon processes that are to be calculated for this absorption process.
        + omegaLess         ::Float64
            ... energy of the photon with smaller frequency; the larger frequency is derived from energy conservation.
    """
    struct   TwoPhotonAbsorptionBichromatic  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProcess
        properties          ::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1}
        omegaLess           ::Float64
    end

    """
    `MultiPhotonDeExcitation.TwoPhotonAbsorptionBichromatic()`
        ... constructor for an 'default' instance of a MultiPhotonDeExcitation.TwoPhotonAbsorptionBichromatic() process.
    """
    function TwoPhotonAbsorptionBichromatic()
        TwoPhotonAbsorptionBichromatic( [TotalCsUnpolarized()], 0.05 )
    end

    function Base.string(process::MultiPhotonDeExcitation.TwoPhotonAbsorptionBichromatic)
        sa = "Two-photon absorption process with bi-chromatic photons (from different beams):"
        return( sa )
    end

    function Base.show(io::IO, process::MultiPhotonDeExcitation.TwoPhotonAbsorptionBichromatic)
        sa = Base.string(process);               print(io, sa, "\n")
        println(io, "properties:                 $(process.properties)  ")
        println(io, "omegaLess:                  $(process.omegaLess)  ")
    end

    """
    `struct  MultiPhotonDeExcitation.TwoPhotonEmission  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProcess`
        ... a struct to specify a two-photon emission process.

        + properties        ::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1}
            ... List of the multi-photon processes that are to be calculated for this emission process.
        + noSharings  ::Int64   ... number of energy sharings for atomic transition, using the zeros of a Gauß-Legendre integration
    """
    struct   TwoPhotonEmission  <:  MultiPhotonDeExcitation.AbstractMultiPhotonProcess
        properties          ::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1}
        noSharings          ::Int64
    end

    """
    `MultiPhotonDeExcitation.TwoPhotonEmission()`
        ... constructor for an 'default' instance of a MultiPhotonDeExcitation.TwoPhotonEmission() process.
    """
    function TwoPhotonEmission()
        TwoPhotonEmission( [EnergyDiffCs()], 4 )
    end

    function Base.string(process::MultiPhotonDeExcitation.TwoPhotonEmission)
        sa = "Two-photon emission process:"
        return( sa )
    end

    function Base.show(io::IO, process::MultiPhotonDeExcitation.TwoPhotonEmission)
        sa = Base.string(process);               print(io, sa, "\n")
        println(io, "properties:                 $(process.properties)  ")
    end

    """
    `struct  MultiPhotonDeExcitation.Settings  <:  AbstractProcessSettings`  ... defines a type for the settings in estimating multi-photon excitation and decay rates.

        + process           ::MultiPhotonDeExcitation.AbstractMultiPhotonProcess      ... considered multi-photon process.
        + multipoles        ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges            ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + green             ::Array{AtomicState.GreenChannel,1}
                              ... Green function channels for these multi-photon calculations which need to be generated independently.
        + NoEnergySharings  ::Int64                        ... Number of energy sharings that are used for all emission lines in the
                                                               computations for each line.
        + printBefore       ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + lineSelection     ::LineSelection                ... Specifies the selected levels, if any.

    """
    struct Settings  <:  AbstractProcessSettings
        process             ::MultiPhotonDeExcitation.AbstractMultiPhotonProcess
        multipoles          ::Array{EmMultipole}
        gauges              ::Array{UseGauge}
        green               ::Array{AtomicState.GreenChannel,1}
        NoEnergySharings    ::Int64
        printBefore         ::Bool
        lineSelection       ::LineSelection
    end

    """
    `MultiPhotonDeExcitation.Settings()`  ... constructor for the default values of multi-photon excitation and decay rates.
    """
    function Settings()
        Settings( TwoPhotonEmission(), EmMultipole[], UseGauge[], AtomicState.GreenExpansion(), false, LineSelection() )
    end

    # `Base.show(io::IO, settings::MultiPhotonDeExcitation.Settings)`
    #		... prepares a proper printout of the variable settings::MultiPhotonDeExcitation.Settings.
     function Base.show(io::IO, settings::MultiPhotonDeExcitation.Settings)
        println(io, "process:                  $(settings.process)  ")
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        print(io,   "green:                     (settings.green)  ... with symmetries:   ")
        for channel in settings.green  print(io, "   $(channel.symmetry)")  end;   println(io, " ")
        println(io, "NoEnergySharings:         $(settings.NoEnergySharings)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "lineSelection:            $(settings.lineSelection)  ")
    end

    include("module-MultiPhotonDeExc-inc-2pAbsMono.jl")
    include("module-MultiPhotonDeExc-inc-2pEmission.jl")

end # module


text:
#  test MultipleFeatures2D
using RoME
using Test

using Optim

# using Gadfly

x1 = [0.0;0;0]
x2 = [1.0;-1.0;pi/2]

l1 = [0.5;0.5]
l2 = [1.0;0.0]
l3 = [1.5;-0.5]
l3b = [1.5;0.0]
L = [l1';l2';l3']

# plot(x=L[:,1],y=L[:,2],Geom.point)

the11 = atan(l1[2],l1[1])
the12 = atan(l2[2],l2[1])
the13 = atan(l3[2],l3[1])

Xi = SE2(x1)
Xj = SE2(x2)
Δxij = se2vee(Xi\Xj)

tl1 = se2vee(Xj \ SE2([l1;0.0]))[1:2]
the21 = atan(tl1[2],tl1[1])

tl2 = se2vee(Xj \ SE2([l2;0.0]))[1:2]
the22 = atan(tl2[2],tl2[1])

tl3 = se2vee(Xj \ SE2([l3;0.0]))[1:2]
the23 = atan(tl3[2],tl3[1])

tl3b = se2vee(Xj \ SE2([l3b;0.0]))[1:2]
the23b = atan(tl3b[2],tl3b[1])

xir1 = Normal(the11,1e-4)
xir2 = Normal(the12,1e-4)
xir3 = Normal(the13,1e-4)
xjr1 = Normal(the21,1e-4)
xjr2 = Normal(the22,1e-4)
xjr3 = Normal(the23,1e-4)
xjr3b = Normal(the23b,1e-4)

# Not using second feature
mm2 = MultipleFeatures2D(
  xir1,
  xir2,
  xir3,
  xjr1,
  xjr2,
  xjr3,
  xjr3b,
  Categorical([1.0;0.0]),
  SE2(zeros(3))  )

global res = 0.1*randn(3)
idx = 1
meas = getSample(mm2)

wAbi = reshape(x1,3,1)
wAbj = reshape(x2,3,1)
wAo1 = reshape(l1,2,1)
wAo2 = reshape(l2,2,1)
wAo3 = reshape(l3,2,1)

@testset "test MultipleFeatures unit vector functions..." begin

res, rhat, resid = zeros(3),zeros(2), zeros(2)

for i in 1:3
  r1, rhathat1, α1 = getUvecScaleFeature2D(Xi, SE2(zeros(3)), SE2([L[i,:];0.0]))
  rhat[1] = cos(meas[1][i,idx,1])
  rhat[2] = sin(meas[1][i,idx,1])
  # @show rhat, rhathat1
  resid[1:2] = rhat-rhathat1
  res[1] += dot(resid, resid)
end

@test 0.0 <= res[1] < 1e-5

end

@testset "test residual function..." begin

global res

res = zeros(3)

wAbi[:] .= 0.0
# @show wAbi

mm2(res,
  nothing,
  idx,
  meas,
  wAbi,
  wAbj,
  wAo1,
  wAo2,
  wAo3  )

# @show res

function gg2(x::Float64,y::Float64)
  wAbi[2] = x
  wAbi[3] = y

  res = zeros(3)

  mm2(res, nothing, idx, meas,
    wAbi,
    wAbj,
    wAo1,
    wAo2,
    wAo3  )
  return sum(res)[1]
end
XX = range(-0.5,stop=0.5,length=200)
YY = range(-pi/4,stop=pi/4,length=100)
gg2(0.0,0.0)
# Gadfly.plot(z=gg2,x=XX,y=YY, Geom.contour)

end

function minmickey(x::Array{Float64,1})
  global res
  wAbi[1:3,1] = x[1:3]
  res = zeros(3)
  mm2(res, nothing, idx, meas,
    wAbi,
    wAbj,
    wAo1,
    wAo2,
    wAo3  )

  return sum(res[1:2])
end

# ar99
@testset "test as minimization problem to pose..." begin

global res

gg = (x,y) -> minmickey([x;y;0.0])
# gg = (x,y) -> x*y
# plot(z=gg, x=range(-2,stop=2,length=100),y=range(-2,stop=2,length=100), Geom.contour)

gg(0.0,0.1)

@time r = optimize(minmickey, randn(3))
@show r.minimizer

@test norm(res) < 1e-3

end

function minmickey(x::Array{Float64,1})
  wAo1[1:2,1] = x[1:2]
  res = zeros(3)
  mm2(res, nothing, idx, meas,
    wAbi,
    wAbj,
    wAo1,
    wAo2,
    wAo3  )

  return sum(res[1:2])
end

@testset "test as minimization problem to landmark..." begin

wAbi[:] .= 0.0

gg = (x,y) -> minmickey([x;y])
#plot(z=gg, x=range(1,stop=2,length=100),y=range(-1,stop=1,length=100), Geom.contour(levels=50))

# now build in factor graph form for further testing

N=100 # TODO return to 50
fg = initfg()

initCov = 0.05*Matrix{Float64}(LinearAlgebra.I, 3,3)
initCov[3,3] = 0.001
odoCov = deepcopy(initCov)

ipp = PriorPose2(MvNormal(zeros(3), initCov^2))

v1 = addVariable!(fg, :x1, Pose2, N=N)

f1  = addFactor!(fg,[:x1], ipp)

v2 = addVariable!(fg, :x2, Pose2, N=N)

px2 = zeros(3)
px2[1:3] = [1.1;-1.0;pi/2]
ipp2 = PriorPose2(MvNormal(px2, initCov^2))
f1  = addFactor!(fg,[:x2], ipp2)

vl1 = addVariable!(fg, :l1, Point2, N=N)
vl2 = addVariable!(fg, :l2, Point2, N=N)
vl3 = addVariable!(fg, :l3, Point2, N=N)

# # why is explicit call to autoinit required here? should not be necessary
# doautoinit!(fg, :x1)
# doautoinit!(fg, :x2)

f2 = addFactor!(fg, [v1;v2;vl1;vl2;vl3], mm2, threadmodel=SingleThreaded )

# getVal(fg, :l3)
#
# writeGraphPdf(fg)

data = getData(f2)
# fieldnames(data.fnc)
@test data.fnc.threadmodel == SingleThreaded

ef2pts = approxConv(fg, :x1x2l1l2l3f1, :l2)
# evalFactor(fg, f2, fg.IDs[:l2])

tree, smt, hist = solveTree!(fg)
# ensureAllInitialized!(fg)
# tree = wipeBuildNewTree!(fg)
# # spyCliqMat(tree.cliques[1])
#
# batchSolve!(fg)

end

# lets see what is happening during MCMC runs
# plotMCMC(tree, :l2, show=true, levels=3) # make true to see pictures, false for testing

# spyCliqMat(tree, :x2)

#
# pts, = getSample(ipp2, 100);
#
# plotKDE( kde!(pts[1:2,:]) )
#
# @show ipp2.Cov
# pts = rand(MvNormal(ipp2.Zi[:,1],ipp2.Cov),N);

# using KernelDensityEstimate
# plotKDE(kde!(ef2pts[1,:]))

# pl = plotKDE([kde!(randn(2,100));kde!(2+randn(2,100))],levels=1,c=["red";"green"],legend=["1";"2"])

# arrp = [getBelief(fg, :l1);getBelief(fg, :l2);getBelief(fg, :l3)]
# arrp1 = deepcopy(arrp)
# plotKDE([arrp;arrp1],levels=1,c=["red";"red";"red";"blue";"blue";"blue"])

# arrp=[marginal(getBelief(fg, :x1),[1;2]);marginal(getBelief(fg, :x2),[1;2]) ]
# plotKDE(arrp,levels=1)

# plotKDE(marginal(getBelief(fg,:x1),[3]))

# pp,arr,parts = localProduct(fg,:x2)
# cc = plotKDE([pp;arr],c=["red";"black";"blue"],levels=1);

# cc = drawLbl(fg,:l1);
# Gadfly.draw(PDF("/home/dehann/Desktop/test.pdf",20cm,15cm),cc)

@testset "test ambiguous bi-modal multifeature constraint operation" begin|

mm2 = MultipleFeatures2D(
  xir1,
  xir2,
  xir3,
  xjr1,
  xjr2,
  xjr3,
  xjr3b,
  Categorical([0.5;0.5]),
  SE2(zeros(3))  )

end

#


text:
while true
    str = readline()
    str == "-" && break
    m = parse(Int, readline())
    for i in 1:m
        h = parse(Int, readline())
        str = str[h+1:end] * str[1:h]
    end
    println(str)
end


text:
#using Compat.Test
#using Compat

using SchwarzChristoffel
using Test

include("Exterior.jl")
include("Polygons.jl")

#@testset ExtendedTestSet "All tests" begin
#  @includetests ARGS
#end

# if isempty(ARGS)
#     include("../docs/make.jl")
# end


text:
##############################################################################################
module ParallelExperimentDriver
##############################################################################################
# imports and exports
using Distributed
@everywhere push!(LOAD_PATH, "/DataAssimilationBenchmarks/src/")
@everywhere push!(LOAD_PATH, "/DataAssimilationBenchmarks/src/methods")
@everywhere push!(LOAD_PATH, "/DataAssimilationBenchmarks/src/models")
@everywhere push!(LOAD_PATH, "/DataAssimilationBenchmarks/src/experiments")
@everywhere using JLD, ..FilterExps, ..SmootherExps, ..EnsembleKalmanSchemes, ..DeSolvers,
                  ..L96, ..ParallelExperimentDriver
@everywhere export experiment
#@everywhere export wrap_exp

##############################################################################################
##############################################################################################
## Timeseries data
##############################################################################################
# observation timeseries to load into the experiment as truth twin, timeseries are named by
# the model, seed to initialize, the integration scheme used to produce, number of analyses,
# the spinup length, and the time length between observation points

ts1 = "../data/time_series/L96_time_series_seed_0000_dim_40_diff_0.00_F_08.0_tanl_0.05_nanl_50000_spin_5000_h_0.010.jld2"
ts2 = "../data/time_series/L96_time_series_seed_0000_dim_40_diff_0.00_F_08.0_tanl_0.10_nanl_50000_spin_5000_h_0.010.jld2"
ts3 = "../data/time_series/L96_time_series_seed_0000_dim_40_diff_0.00_F_08.0_tanl_0.15_nanl_50000_spin_5000_h_0.010.jld2"
ts4 = "../data/time_series/L96_time_series_seed_0000_dim_40_diff_0.00_F_08.0_tanl_0.20_nanl_50000_spin_5000_h_0.010.jld2"
ts5 = "../data/time_series/L96_time_series_seed_0000_dim_40_diff_0.00_F_08.0)tanl_0.25_nanl_50000_spin_5000_h_0.010.jld2"
##############################################################################################

##############################################################################################
## Experiment parameter generation
##############################################################################################
##############################################################################################

##############################################################################################
# Filters
##############################################################################################
##############################################################################################
# filter_state
##############################################################################################
## [time_series, scheme, seed, nanl, obs_un, obs_dim, N_ens, infl] = args
#
#schemes = ["enkf-n-primal", "enkf-n-primal-ls", "enkf-n-dual"]
#seed = 0
#obs_un = 1.0
#obs_dim = 40
#N_ens = 15:43
#infl = [1.0]#LinRange(1.0, 1.20, 21)
#nanl = 2500
#
## load the experiments
#args = Tuple[]
#for scheme in schemes
#    for N in N_ens
#        for α in infl
#            tmp = (time_series, scheme, seed, nanl, obs_un, obs_dim, N, α)
#            push!(args, tmp)
#        end
#    end
#end
#
#experiment = FilterExps.filter_state
#
#
##############################################################################################
# filter_param
##############################################################################################
## [time_series, scheme, seed, nanl, obs_un, obs_dim, param_err, param_wlk, N_ens,
##  state_infl, param_infl] = args
#
#schemes = ["enkf", "etkf"]
#seed = 0
#obs_un = 1.0
#obs_dim = 40
#param_err = 0.03
#param_wlk = [0.0000, 0.0001, 0.0010, 0.0100]
#N_ens = 14:41
#state_infl = LinRange(1.0, 1.20, 21)
#param_infl = LinRange(1.0, 1.00, 1)
#nanl = 2500
#
## load the experiments
#args = Tuple[]
#for scheme in schemes
#    for wlk in param_wlk
#        for N in N_ens
#            for s_infl in state_infl
#                for p_infl in param_infl
#                    tmp = (time_series, scheme, seed, nanl, obs_un, obs_dim,
#                           param_err, wlk, N, s_infl, p_infl)
#                    push!(args, tmp)
#                end
#            end
#        end
#    end
#end
#
#experiment = FilterExps.filter_param
#
#
##############################################################################################
##############################################################################################
# Classic smoothers
##############################################################################################
# classic_state parallel run, arguments are
# time_series, method, seed, nanl, lag, shift, obs_un, obs_dim, γ, N_ens, state_infl = args

schemes = ["etks"]
seed = 0
lags = 1:3:52
shifts = [1]
#lags = [1, 2, 4, 8, 16, 32, 64]
#gammas = Array{Float64}(1:11)
gammas = [1.0]
shift = 1
obs_un = 1.0
obs_dim = 40
#N_ens = 15:2:41
N_ens = [21]
#state_infl = [1.0]
state_infl = LinRange(1.0, 1.10, 11)
time_series = [ts1, ts2, ts3, ts4, ts5]
nanl = 2500

# load the experiments
args = Tuple[]
for ts in time_series
    for scheme in schemes
        for γ in gammas
            for l in 1:length(lags)
                # optional definition of shifts in terms of the current lag parameter for a
                # range of shift values
                lag = lags[l]
                #shifts = lags[1:l]
                for shift in shifts
                    for N in N_ens
                        for s_infl in state_infl
                            tmp = (ts, scheme, seed, nanl, lag, shift, obs_un, obs_dim,
                                   γ, N, s_infl)
                            push!(args, tmp)
                        end
                    end
                end
            end
        end
    end
end

## define the robust to failure wrapper
#function wrap_exp(arguments)
#    try
#        classic_state(arguments)
#    catch
#        print("Error on " * string(args) * "\n")
#    end
#end

experiment = SmootherExps.classic_state
#experiment = wrap_exp

##############################################################################################
## classic_param single run for debugging, arguments are
##  [time_series, method, seed, nanl, lag, shift, obs_un, obs_dim, param_err,
##   param_wlk, N_ens, state_infl, param_infl = args
#
#schemes = ["enks", "etks"]
#seed = 0
#lag = 1:5:51
#shift = 1
#obs_un = 1.0
#obs_dim = 40
#N_ens = 14:41
#param_err = 0.03
#param_wlk = [0.0000, 0.0001, 0.0010, 0.0100]
#state_infl = LinRange(1.0, 1.20, 21)
#param_infl = LinRange(1.0, 1.00, 1)
#nanl = 2500
#
## load the experiments
#args = Tuple[]
#for scheme in schemes
#    for l in lag
#        for N in N_ens
#            for wlk in param_wlk
#                for s_infl in state_infl
#                    for p_infl in param_infl
#                        tmp = (time_series, scheme, seed, nanl, l, shift, obs_un, obs_dim,
#                               param_err, wlk, N, s_infl, p_infl)
#                        push!(args, tmp)
#                    end
#                end
#            end
#        end
#    end
#end
#
#experiment = SmootherExps.classic_param
#
#
##############################################################################################

##############################################################################################
# Single iteration smoothers
##############################################################################################
## single iteration single run for degbugging, arguments are
## [time_series, method, seed, nanl, lag, shift, mda, obs_un, obs_dim,
##  N_ens, state_infl = args
#
#schemes = ["enks-n-primal"]
#seed = 0
##lags = 1:3:52
#lags = [1, 2, 4, 8, 16, 32, 64]
##gammas = Array{Float64}(1:11)
#gammas = [1.0]
##shift = 1
#obs_un = 1.0
#obs_dim = 40
##N_ens = 15:2:41
#N_ens = [21]
#state_infl = [1.0]
##state_infl = LinRange(1.0, 1.10, 11)
#time_series = [time_series_1, time_series_2]
#mdas = [false]
#nanl = 2500
#
## load the experiments
#args = Tuple[]
#for m in mdas
#    for ts in time_series
#        for γ in gammas
#            for scheme in schemes
#                for l in 1:length(lags)
#                    # optional definition of shift in terms of the current lag parameter
#                    # for a range of shift values
#                    lag = lags[l]
#                    shifts = lags[1:l]
#                    for shift in shifts
#                        for N in N_ens
#                            for s_infl in state_infl
#                                tmp = (ts, scheme, seed, nanl, lag, shift, m, obs_un,
#                                       obs_dim, γ, N, s_infl)
#                                push!(args, tmp)
#                            end
#                        end
#                    end
#                end
#            end
#        end
#    end
#end
#
## define the robust to failure wrapper
#function wrap_exp(arguments)
#    try
#        SmootherExps.single_iteration_state(arguments)
#    catch
#        print("Error on " * string(args) * "\n")
#    end
#end
#
#experiment = wrap_exp
#
##############################################################################################
# Run the experiments in parallel over the parameter values
##############################################################################################
pmap(experiment, args)

##############################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
using Refactoring
using Test

@testset "Refactoring.jl" begin

    ex = Refactoring.econvert(Meta.parse("x[i] = 1"))
    @test ex.lhs.lhs == :x
    @test ex.lhs.rhs == :i

    include("test_extract_method.jl")
    include("test_search_and_replace.jl")

end


text:
"""
    print_tree(tree, maxdepth=5; kwargs...)
    print_tree(io, tree, maxdepth=5; kwargs...)
    print_tree(f::Function, io, tree, maxdepth=5; kwargs...)

# Usage
Prints an ASCII formatted representation of the `tree` to the given `io` object.
By default all children will be printed up to a maximum level of 5, though this
value can be overriden by the `maxdepth` parameter. Nodes that are truncated are
indicated by a vertical ellipsis below the truncated node, this indication can be
turned off by providing `indicate_truncation=false` as a kwarg. The charset to use in
printing can be customized using the `charset` keyword argument.
You can control the printing of individual nodes by passing a function `f(io, node)`;
the default is [`AbstractTrees.printnode`](@ref).

# Examples
```julia
julia> print_tree(stdout, Dict("a"=>"b","b"=>['c','d']))
Dict{String,Any}("b"=>['c','d'],"a"=>"b")
├─ b
│  ├─ c
│  └─ d
└─ a
   └─ b

julia> print_tree(stdout, '0'=>'1'=>'2'=>'3', 2)
'0'
└─ '1'
    └─ '2'
        ⋮

julia> print_tree(stdout, Dict("a"=>"b","b"=>['c','d']);
        charset = TreeCharSet('+','\\','|',"--","⋮"))
Dict{String,Any}("b"=>['c','d'],"a"=>"b")
+-- b
|   +-- c
|   \\-- d
\\-- a
   \\-- b
```

"""
print_tree

"""
    printnode(io::IO, node)

Print a single node. The default is to show a compact representation of `node`.
Override this if you want nodes printed in a custom way in [`print_tree`](@ref),
or if you want your print function to print part of the tree by default.

# Example

```
struct MyNode{T}
    data::T
    children::Vector{MyNode{T}}
end
AbstractTrees.printnode(io::IO, node::MyNode) = print(io, "MyNode(\$(node.data))")
```
"""
printnode(io::IO, node) = show(IOContext(io, :compact => true), node)

struct TreeCharSet
    mid
    terminator
    skip
    dash
    trunc
end

# Default charset
TreeCharSet() = TreeCharSet('├','└','│','─','⋮')
TreeCharSet(mid, term, skip, dash) = TreeCharSet(mid, term, skip, dash, '⋮')

function print_prefix(io, depth, charset, active_levels)
    for current_depth in 0:(depth-1)
        if current_depth in active_levels
            print(io,charset.skip," "^(textwidth(charset.dash)+1))
        else
            print(io," "^(textwidth(charset.skip)+textwidth(charset.dash)+1))
        end
    end
end

function _print_tree(printnode::Function, io::IO, tree, maxdepth = 5; indicate_truncation = true,
                     depth = 0, active_levels = Int[], charset = TreeCharSet(), withinds = false,
                     inds = [], from = nothing, to = nothing, roottree = tree)
    nodebuf = IOBuffer()
    isa(io, IOContext) && (nodebuf = IOContext(nodebuf, io))
    if withinds
        printnode(nodebuf, tree, inds)
    else
        tree != roottree && isa(treekind(roottree), IndexedTree) ?
            printnode(nodebuf, roottree[tree]) :
            printnode(nodebuf, tree)
    end
    str = String(take!(isa(nodebuf, IOContext) ? nodebuf.io : nodebuf))
    for (i,line) in enumerate(split(str, '\n'))
        i != 1 && print_prefix(io, depth, charset, active_levels)
        println(io, line)
    end
    c = isa(treekind(roottree), IndexedTree) ?
        childindices(roottree, tree) : children(roottree, tree)
    if c !== ()
        if depth < maxdepth
            s = Iterators.Stateful(from === nothing ? pairs(c) : Iterators.Rest(pairs(c), from))
            while !isempty(s)
                ind, child = popfirst!(s)
                ind === to && break
                active = false
                child_active_levels = active_levels
                print_prefix(io, depth, charset, active_levels)
                if isempty(s)
                    print(io, charset.terminator)
                else
                    print(io, charset.mid)
                    child_active_levels = push!(copy(active_levels), depth)
                end
                print(io, charset.dash, ' ')
                print_tree(printnode, io, child, maxdepth;
                indicate_truncation=indicate_truncation, depth = depth + 1,
                active_levels = child_active_levels, charset = charset, withinds=withinds,
                inds = withinds ? [inds; ind] : [], roottree = roottree)
            end
        elseif indicate_truncation
            print_prefix(io, depth, charset, active_levels)
            println(io, charset.trunc)
            print_prefix(io, depth, charset, active_levels)
            println(io)
        end
    end
end
print_tree(f::Function, io::IO, tree, args...; kwargs...) = _print_tree(f, io, tree, args...; kwargs...)
print_tree(io::IO, tree, args...; kwargs...) = print_tree(printnode, io, tree, args...; kwargs...)
print_tree(tree, args...; kwargs...) = print_tree(stdout::IO, tree, args...; kwargs...)

"""
Get the string result of calling [`print_tree`](@ref) with the supplied arguments.

The `context` argument works as it does in `Base.repr`.
"""
function repr_tree(tree, args...; context=nothing, kw...)
    buf = IOBuffer()
    io = context === nothing ? buf : IOContext(buf, context)
    print_tree(io, tree, args...; kw...)
    return String(take!(buf))
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, LidJul

makedocs(
    modules = [LidJul],
    format = Documenter.HTML(prettyurls = false),
    checkdocs = :exports,
    sitename = "LidJul.jl",
    pages = Any["index.md"]
)

# deploydocs(
#     repo = "github.com/{GHUSER}/LidJul.jl.git",
# )


text:
# this program illustrates sample selection bias associated with
# dropping observations for which the dep. vbl. is <0. The resulting figure is
# saved in sampsel.ps
using Plots
pyplot()
n = 100
sig = 5
x = 10*rand(n,1)
y = x+sig*randn(n,1)
# this drops the rows for which y < 0
z = y.>0.0
yin = y[z]
xin = x[z]
# regression using selected sample
xin = [ones(size(xin)) xin]
b = xin\yin
# plots
xx=(0:0.05:10)
yy=xx
yhat = b[1]+b[2]*xx
scatter(x,y, label = "data")
plot!(xx,yy,label = "population regression line")
plot!(xx,yhat,label = "fitted regression line")
gui()
savefig("sampsel.png")

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
mutable struct Accumulate{V,T,R} <: State{V}
    value::V
    time::T
    rate::R
    reset::Bool
end

Accumulate(; unit, time, timeunit, _value, _type, _...) = begin
    U = value(unit)
    V = valuetype(_type, U)
    v = _value
    #V = promote_type(V, typeof(v))
    TU = value(timeunit)
    t = unitfy(value(time), TU)
    T = typeof(t)
    RU = rateunittype(U, TU)
    R = valuetype(_type, RU)
    Accumulate{V,T,R}(v, t, zero(R), false)
end

constructortags(::Val{:Accumulate}) = (:unit, :init, :time, :timeunit)

@generated rateunit(::Accumulate{V,T,R}) where {V,T,R} = unittype(R)

updatetags!(d, ::Val{:Accumulate}; _...) = begin
    N = d[:_type]
    U = d[:unit]
    !haskey(d, :init) && (d[:init] = @q zero($N))
    !haskey(d, :time) && (d[:time] = :(context.clock.time))
    #TODO: automatic inference without explicit `timeunit` tag
    !haskey(d, :timeunit) && (d[:timeunit] = @q $C.timeunittype($U, $C.timeunit(__Context__)))
    !haskey(d, :reset) && (d[:reset] = false)
end

genvartype(v::VarInfo, ::Val{:Accumulate}; N, U, V, _...) = begin
    TU = gettag(v, :timeunit)
    T = @q $C.valuetype(Float64, $TU)
    RU = @q $C.rateunittype($U, $TU)
    R = @q $C.valuetype($N, $RU)
    @q Accumulate{$V,$T,$R}
end

gendefault(v::VarInfo, ::Val{:Accumulate}) = begin
    i = gettag(v, :init)
    u = gettag(v, :unit)
    @q $C.unitfy($C.value($i), $C.value($u))
end

genupdate(v::VarInfo, ::Val{:Accumulate}, ::MainStep; kw...) = begin
    @gensym s a0 t t0 a
    @q let $s = $(symstate(v)),
           $a0 = $s.reset ? $(gendefault(v)) : $s.value,
           $t = $C.value($(gettag(v, :time))),
           $t0 = $s.time,
           $a = $a0 + $s.rate * ($t - $t0)
        $(genstore(v, a; unitfy=false, minmax=true, round=false, when=false))
    end
end

genupdate(v::VarInfo, ::Val{:Accumulate}, ::PostStep; kw...) = begin
    w = gettag(v, :when)
    f = isnothing(w) ? genbody(v) : @q $C.value($w) ? $(genbody(v)) : zero($(gettag(v, :_type)))
    @gensym s t r e
    @q let $s = $(symstate(v)),
           $t = $C.value($(gettag(v, :time))),
           $r = $C.unitfy($f, $C.rateunit($s)),
           $e = $C.value($(gettag(v, :reset)))
        $s.time = $t
        $s.rate = $r
        $s.reset = $e
    end
end


text:
using JuMP

m = Model()
@variable(m, 0<=x<=2)
@variable(m,0<=y<=30)

@objective(m, Max, 5x+3*y)
@constraint(m, 1x+5y>=3.0)

print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))


text:
using DataFrames, Random, BenchmarkTools
using Douglass
using CategoricalArrays

rng = MersenneTwister(12345)

# This is similar to Mathieu Gomez's setup for benchmarking FixedEffectModels
N = 1_000_000
K = 100
id1 = rand(rng,1:div(N, K), N)
id2 = rand(rng, 1:K, N)
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)

df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
allowmissing!(df)

df.x2[rand(rng, 1:10,N) .> 8] .= missing
sort!(df, [:id2, :id1])

set_active_df(:df)
@time d"gen :z = :x1 + :x2"
@time d"drop :z"

# all times are from my 2020 MacBook Pro:
# julia> versioninfo()
# Julia Version 1.6.2
# Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
# Platform Info:
#   OS: macOS (x86_64-apple-darwin18.7.0)
#   CPU: Intel(R) Core(TM) i9-9980HK CPU @ 2.40GHz
#   WORD_SIZE: 64
#   LIBM: libopenlibm
#   LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
# Environment:
#   JULIA_EDITOR = code
#   JULIA_NUM_THREADS =
# Stata is StataMP 15.1 (see code)

gen = @benchmark begin
    d"gen :z = :x1 + :x2"
    d"drop :z"
end

egen = @benchmark begin
    d"bysort :id2 (:id1): egen :z = mean(:x1)"
    d"drop :z"
end

egenif = @benchmark begin
    d"bysort :id2 (:id1): egen :z = mean(:x1) if :x1 .> 0.0"
    d"drop :z"
end

using Statistics
egenpair = @benchmark begin
    d"bysort :id2 (:id1): egen :z = cor(:x1,:x2) if :x1 .> 0.0"
    d"drop :z"
end

# RESHAPING:
rng = MersenneTwister(12345)

N = 1_000_000
K = 100
id1 = floor.(div.(collect(1:N),10))
id2 = mod.(collect(1:N),10)
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)

df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
allowmissing!(df)

sort!(df, [:id2, :id1])

df_cpy = deepcopy(df)
# warm up
d"reshape_wide :x1 :x2 :y, i(:id1) j(:id2)"
d"reshape_long :x1 :x2 :y, i(:id1) j(:id2)"
# now for real
df = deepcopy(df_cpy)
reshapewide = @elapsed d"reshape_wide :x1 :x2 :y, i(:id1) j(:id2)"
reshapelong = @elapsed d"reshape_long :x1 :x2 :y, i(:id1) j(:id2)"

df = deepcopy(df_cpy)
d"duplicates_drop :id1 :id2"
df = deepcopy(df_cpy)
duplicatesdrop = @elapsed d"duplicates_drop :id1 :id2"

# plot the results

using DataFrames, CSV, Gadfly
results = CSV.read("benchmark/resultStata.csv", DataFrame)

rename!(results, :result => :stata)

results[!,:julia] = [median(gen).time/1e9, median(egen).time/1e9, median(egenif).time/1e9, median(egenpair).time/1e9, reshapewide, reshapelong, duplicatesdrop]

results[!,:relative] = results.julia ./ results.stata
results[!,:test] = ["gen", "egen/by", "egen/by\n + if", "egen\n + corr", "reshape\nwide", "reshape\nlong", "duplicates\ndrop"]

p = plot(results, x = "test", y = "relative", Guide.ylabel("Relative Time (Douglass/Stata)"), Guide.xlabel("Command"), Guide.yticks(ticks= [0, 1, 2, 3, 4, 5 ]))

draw(SVG("benchmark/benchmark.svg", 8inch, 5inch), p)


text:
#!/usr/bin/env julia

using Printf

# --- Def ---
include("lexerparser.inc.jl")

mutable struct State
    prods
    transitions
    id
end

state_id = 0

function state_new(prods)
    global state_id

    return State(prods, Dict(), state_id += 1)
end

# Same prods
state_eq_shallow(a, b) = issetequal(a.prods, b.prods)

function Base.print(io::IO, s::State)
    for p in s.prods
        Base.println(io, p)
    end

    for (k, v) in s.transitions
        Base.println("$(k.id) ---> $(v.id)")
    end
end

# --- Algo ---
# Updates the state set (states) by solving state
function slr_graph_solve!(state, states, prods)
    # Get each transition token
    # transitions[tok] = all prods with tok after the dot (aka pos)
    transitions = Dict()

    for p in state.prods
        # Reduce case
        if p.pos == length(p.right) + 1
            continue
        end

        next_tok = p.right[p.pos]
        if !haskey(transitions, next_tok)
            transitions[next_tok] = [p]
        else
            push!(transitions[next_tok], p)
        end
    end

    # For each new state
    for (transition_tok, old_prods) in transitions
        new_prods = Set()

        # Move the dot / remove reduce
        for p in old_prods
            # Reduce
            if p.pos == length(p.right) + 1
                continue
            end

            p = deepcopy(p)
            p.pos += 1
            push!(new_prods, p)
        end

        # Add non terminal prods recursively
        while true
            # All prods that have a non terminal token at pos
            nt_prods = Set([p for p in new_prods
                            if p.pos <= length(p.right) && !p.right[p.pos].terminal])
            # println(nt_prods)
            update = false
            for p in nt_prods
                # Token we create
                tok = p.right[p.pos]

                # Get all prods that make this token
                productions = Set([deepcopy(prod)
                                   for prod in prods if tok_eq(prod.left, tok)])
                # Add them
                old_size = length(new_prods)
                union!(new_prods, productions)

                if length(new_prods) != old_size
                    update = true
                end
            end

            if !update
                break
            end
        end

        # Either new or already present state
        next_state = state_new(new_prods)
        found = false
        for s in states
            # Found same state
            if state_eq_shallow(s, next_state)
                next_state = s
                found = true
                break
            end
        end

        # Link this state
        state.transitions[transition_tok] = next_state
        if !found
            push!(states, next_state)
            slr_graph_solve!(next_state, states, prods)
        end
    end
end

# !!! kernel is not removed from states
function slr_graph(prods)
    kernel = state_new(prods)
    states = [kernel]

    slr_graph_solve!(kernel, states, prods)

    for (i, s) in enumerate(states)
        s.id = i
    end

    return states
end

# Action is one of Sn, Rn, acc, nothing
function slr_table(states, tokens, prods, follow)
    # table[state][tok] = action
    table = Matrix{Any}(nothing, length(states), length(tokens))

    tok2index = Dict()
    for (i, t) in enumerate(tokens)
        tok2index[t] = i
    end

    prod2index = Dict()
    for (i, p) in enumerate(prods)
        prod2index[p] = i
    end

    for s in states
        # Shifts / gotos
        for (tok, next_state) in s.transitions
            table[s.id, tok2index[tok]] = "S$(next_state.id)"
        end

        # Reduce
        for p in s.prods
            if p.pos == length(p.right) + 1
                # Set reduce action for all tokens in follow(left of production)
                for reduce_tok in follow[p.left]
                    toki = tok2index[reduce_tok]
                    if table[s.id, toki] != nothing && length(table[s.id, toki]) > 0 && table[s.id, toki][1] == 'R'
                        println("--- Error Info ---")
                        println(s)
                        error("Reduce-reduce conflict for state $(s.id), token $(tokens[toki])")
                    end

                    if table[s.id, toki] == nothing
                        # Reduce using the prod when pos = 1
                        start_prod = deepcopy(p)
                        start_prod.pos = 1
                        rule = prod2index[start_prod]
                        table[s.id, toki] = "R$(rule)"
                    end
                end

                # Set init prod (root of the AST)
                if p.init
                    table[s.id, length(tokens)] = "acc"
                end
            end
        end
    end

    return table
end

function printtable(table, tokens)
    for row in 0:size(table, 1)
        for col in 0:size(table, 2)
            if col == 0
                @printf("%3d |", row)
            elseif row == 0
                @printf("%3s |", tokens[col])
            else
                val = table[row, col]
                if val == nothing
                    print("    |")
                else
                    @printf("%3s |", val)
                end
            end
        end
        println()
    end
end

function _first!(token, prods, first)
    if haskey(first, token)
        return first[token]
    end

    if token.terminal
        first[token] = Set([token])
        return first[token]
    else
        set = Set()
        for p in prods
            if p.left == token
                # No epsilon supported
                tok = p.right[1]
                if tok != token
                    union!(set, _first!(tok, prods, first))
                end
            end
        end

        first[token] = set

        return set
    end
end

function first_sets(tokens, prods)
    # first[token] = First set for this token
    first = Dict()

    for tok in tokens
        _first!(tok, prods, first)
    end

    return first
end

# No epsilon supported
function follow_sets(tokens, prods, first)
    follow = Dict()
    end_tok = tok_new("\$", terminal=true)
    follow[prods[begin].left] = Set([end_tok])

    changes = true
    while changes
        changes = false

        for p in prods
            for i in 1:length(p.right)
                tok = p.right[i]
                if !haskey(follow, tok)
                    follow[tok] = Set()
                end

                set = follow[tok]
                old_len = length(set)

                if i == length(p.right)
                    # Merge follow of the left token
                    if haskey(follow, p.left)
                        union!(set, follow[p.left])
                    end
                elseif p.right[i + 1].terminal
                    # Add this terminal to follow
                    push!(set, p.right[i + 1])
                else
                    # Merge first of the next token
                    if haskey(first, p.right[i + 1])
                        union!(set, first[p.right[i + 1]])
                    end
                end

                if old_len != length(set)
                    changes = true
                end
            end
        end
    end

    return follow
end

function setup!(tokens, prods)
    # Tokens
    @assert length(tokens) >= 1 "No tokens"

    end_tok = tok_new("\$", terminal=true)
    push!(tokens, end_tok)

    @assert length(unique(tokens)) == length(tokens) "Duplicate tokens or \$ (end) token declared manually"

    # Productions
    @assert length(prods) >= 1 "No productions"

    root_tok = prods[begin].left

    @assert all(p -> p.left != root_tok, prods[2:end]) "There must be only one rule to produce the root token ($root_tok)"
    prods[begin].init = true

    for (i, p) in enumerate(prods)
        @assert p.left in tokens "Left token ($(p.left)) not declared for production #$i ($p)"

        for (ti, tok) in enumerate(p.right)
            @assert tok in tokens "Right token #$ti ($tok) not declared for production #$i ($p)"
        end
    end
end

# Main function
function generate_parser(io, tokens, prods, produce_rules)
    parserdir = dirname(@__FILE__)

    # Setup and verify everything
    setup!(tokens, prods)

    # First / follow
    first = first_sets(tokens, prods)
    follow = follow_sets(tokens, prods, first)

    # Generate automaton
    states = slr_graph(prods)
    table = slr_table(states, tokens, prods, follow)

    # Def table
    src_table = "table = $table"
    src_tokens = "tokens = $tokens"
    src_prods = "prods = $prods"
    src_produce_rules = "produce_rules = [$(join(produce_rules, ", "))]"

    # Parse functions
    src_parse = read("$parserdir/slr.parse_template.jl", String)

    srcs = [
            src_table,
            src_tokens,
            src_prods,
            src_produce_rules,
            src_parse,
        ]

    src = join(srcs, "\n")
    println(io, src)
end


text:
import jInv.ForwardShare.getSensMatSize

function  getSensMatSize(pFor::Union{UnconstrainedOEDParam, EqconstrainedOEDParam, IneqconstrainedOEDParam})
	 (m,n,nrhs) = getProblemSize(pFor.Design)
	 return n*nrhs,getNumberOfDesignParameters(pFor.Design)
end

import jInv.ForwardShare.getSensMatVec
function getSensMatVec(x::Vector,p::Vector,pFor::Union{UnconstrainedOEDParam,EqconstrainedOEDParam})
	M   = getM(p,pFor.Design,false)
	D   = getD(p,pFor.Design)

	(m,n,nrhs) = getProblemSize(pFor.Design)
	nc   = getNumberOfConstraints(pFor)
	mv   = zeros(n,nrhs)
	for k=1:nrhs
		# get derivatives of design
		dpMx   = getdpMx(p,pFor.Xr[:,k],pFor.Design)
		dpMTMx = getdpMTx(p,M*pFor.Xr[:,k],pFor.Design)
		dpMTx  = getdpMTx(p,D[:,k],pFor.Design)
	 	dpD    = getdpD(p,k,pFor.Design)

		dwHx = M'*(dpMx*x) + dpMTMx*x
	    dwb  = M'*(dpD*x)  + dpMTx*x
		mv[:,k] = (pFor.Ainv\[dwb-dwHx;zeros(nc)])[1:n]
	end
    return vec(mv)
end

import jInv.ForwardShare.getSensTMatVec
function getSensTMatVec(x::Vector,p::Vector,pFor::Union{UnconstrainedOEDParam,EqconstrainedOEDParam})

	M    = getM(p,pFor.Design,false)
	D    = getD(p,pFor.Design)
    n    = size(M,2)
    nrhs = size(D,2)
	nc   = getNumberOfConstraints(pFor)

    Hinvx = (pFor.Ainv\[reshape(x,n,nrhs);zeros(nc,nrhs)])[1:n,:]

	mv = zeros(getNumberOfDesignParameters(pFor.Design))
	for k=1:nrhs
		# get derivatives of design
		dpMx   = getdpMx(p,pFor.Xr[:,k],pFor.Design)
		dpMTMx = getdpMTx(p,M*pFor.Xr[:,k],pFor.Design)
		dpMTx  = getdpMTx(p,D[:,k],pFor.Design)
	 	dpD    = getdpD(p,k,pFor.Design)

    	dwHx = dpMx'*(M*Hinvx[:,k]) + dpMTMx'*Hinvx[:,k]
		dwb  = dpD'*(M*Hinvx[:,k]) + dpMTx'*Hinvx[:,k]
		mv  += dwb-dwHx
	end
    return mv
end


text:
# Tutorial 2.7. Spin textures
# ===========================
#
# Physics background
# ------------------
#  - Spin textures
#  - Skyrmions
#
# Kwant features highlighted
# --------------------------
#  - operators
#  - plotting vector fields

sigma_0 = [1 0; 0 1]
sigma_x = [0 1; 1 0]
sigma_y = [0 -1im; 1im 0]
sigma_z = [1 0; 0 -1]

# vector of Pauli matrices σ_αiβ where greek
# letters denote spinor indices
sigma = [sigma_x, sigma_y, sigma_z]

field_direction(pos; r0, delta) = field_direction(pos,r0,delta)
function field_direction(pos, r0, delta)
    x, y = pos
    r = norm(pos)
    r_tilde = (r - r0) / delta
    theta = (tanh(r_tilde) - 1) * (pi / 2)

    if r == 0
        m_i = [0, 0, -1]
    else
        m_i = [
            (x / r) * sin(theta),
            (y / r) * sin(theta),
            cos(theta),
        ]
    end

    return m_i
end

py"""
def lead_onsite(site, J):
    return J * $sigma_z
"""

lat = kwant.lattice.square(norbs=2)

function make_system(L=80)

    syst = kwant.Builder()

    function square(pos)
        return all(-L/2 < p < L/2 for p in pos)
    end

    py"""
    import numpy as np
    def scattering_onsite(site, r0, delta, J):
        return J*sigma_Z# m_i = [1,1,1]
        # return J + np.dot(m_i,$sigma)
    """
    # m_i = $field_direction(site.pos, r0, delta)
    # return J * (m_i[0]*$(sigma[1]) + m_i[1]*$(sigma[2]) + m_i[2]*$(sigma[3]))

    syst[lat.shape(square, (0, 0))] = py"scattering_onsite"#J*sigma_z#
    syst[lat.neighbors()] = -sigma_0

    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
                         conservation_law=-sigma_z)

    lead[lat.shape(square, (0, 0))] = sigma_z#py"lead_onsite"
    lead[lat.neighbors()] = -sigma_0

    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    return syst
end

#NOTE: for some reason s.tag becomes s[2] in converting to julia. might work on this.
function plot_vector_field(syst; params...)
    xmin, ymin = minimum(s[2] for s in syst.sites)
    xmax, ymax = maximum(s[2] for s in syst.sites)
    x = broadcast((x,y)->x,range(xmin,stop=xmax),range(ymin,stop=ymax)')
    y = broadcast((x,y)->y,range(xmin,stop=xmax),range(ymin,stop=ymax)')

    m_i = map((w,z)->field_direction((w,z); params...), x,y)
    m_i_x, m_i_y, m_i_z = map(w->map(z->z[w],m_i),(1,2,3))

    fig, ax = plt.subplots(1, 1)
    im = ax.quiver(x, y, m_i_x,m_i_y,m_i_z, pivot="mid", scale=75)
    fig.colorbar(im)
    plt.gcf()
end

function plot_densities(syst, densities)
    fig, axes = plt.subplots(1, length(densities))
    for (ax, (title, rho)) in zip(axes, densities)
        kwant.plotter.map(syst, rho, ax=ax, a=4)
        ax.set_title(title)
    end
    plt.gcf()
end

function plot_currents(syst, currents)
    fig, axes = plt.subplots(1, length(currents))
    # if !hasattr(axes, "__len__")
    if length(axes)>1
        axes = (axes,)
    end
    for (ax, (title, current)) in zip(axes, currents)
        kwant.plotter.current(syst, current, ax=ax, colorbar=false)
        ax.set_title(title)
    end
    plt.gcf()
end


text:
import SatSolver: dpll, brute_force, is_formula_satisfiable, SatSolverResult

@testset "DPLL algorithm" begin
  function dpll_test(fla::Formula, ref::SatSolverResult)
    dpll_res = dpll(fla)
    correct_result = is_formula_satisfiable(dpll_res) == is_formula_satisfiable(ref)
    if correct_result && is_formula_satisfiable(ref)
      @test correct_result && all((cl) -> is_clause_satisfied(cl, dpll_res.satisfying_assignment), fla)
    else
      @test correct_result
    end
  end

  @testset "Random formulas" begin
    for _ in 1:500
      fla = gen_rand_formula(;
                             max_clauses=Unsigned(50),
                             max_variables=Unsigned(20),
                             max_clause_len=Unsigned(4))
      dpll_test(fla, brute_force(fla))
    end
  end
end


text:
# This file contains the code for differentiating the Cholesky decomposition

function level2partition(A::AbstractMatrix, j::Int)
    @inbounds begin
        N = size(A, 1)
        r = view(A, j, 1:j-1)
        d = view(A, j, j)
        B = view(A, j+1:N , 1:j-1)
        c = view(A, j+1:N, j)
    end
    return r, d, B, c
end

# Derivative of Cholesky decomposition, see Murray(2016). Differentiation of the Cholesky decomposition. arXiv.1602.07527
function chol_unblocked_rev!(L::AbstractMatrix, A_bar::AbstractMatrix)
    N = size(A_bar, 1)
    @inbounds for j in N:-1:1
        r, d, B, c = level2partition(L, j)
        r_bar, d_bar, B_bar, c_bar = level2partition(A_bar, j)
        d_bar[1] -= dot(c_bar, c) / d[1]
        d_bar[1] /= d[1]

        for i in eachindex(c_bar)
            c_bar[i] /= d[1]
        end

        for i in eachindex(r_bar)
            r_bar[i] -= d_bar[1]*r[i]
        end

        # r_bar .-= B'c_bar
        BLAS.gemv!('T', -1.0, B, c_bar, 1.0, r_bar) # r_bar = r_bar - Bᵀ c_bar
        # B_bar .-= c_bar * r'
        BLAS.ger!(-1.0, c_bar, r, B_bar) # B_bar = B_bar - c_bar * r'
        d_bar[1] /= 2
    end

    return A_bar
end


text:
######### fitting an Ornstein–Uhlenbeck process  ###########

using DataFrames, MCMC

# generate serie
srand(1)
duration = 1000  # 1000 time steps
mu0 = 10.        # convergence value
tau0 = 20        # convergence time
sigma0 = 0.1     # noise term

x = Array(Float64, duration)
x[1] = 1.
for i in 2:duration
	x[i] = x[i-1]*exp(-1/tau0) + mu0*(1-exp(-1/tau0)) +  sigma0*randn()
end

# model definition
ex = quote
    tau ~ Uniform(0, 100)
    sigma ~ Uniform(0, 2)
    mu ~ Uniform(0, 20)

    fac = exp(- 1. / tau)
    resid = x[2:end] - x[1:end-1] * fac - mu * (1. - fac)
    resid ~ Normal(0, sigma)
end

m = model(ex, tau=0.05, sigma=1., mu=1., gradient=true)
m.scale = [1000., 1., 10.]  # scale hint for tau, sigma and mu, to help sampling

# run random walk metropolis (10000 steps, 1000 for burnin, setting initial values)
mcchain01 = run(m * RAM() * SerialMC(1000:10000))

describe(mcchain01)
acceptance(mcchain01) # acceptance rate

# run Hamiltonian Monte-Carlo (10000 steps, 1000 for burnin, 5 inner steps, 0.002 inner step size)
mcchain02 = run(m * HMC(5, 0.002) * SerialMC(1000:10000))

# run NUTS - HMC (1000 steps, 500 for burnin)
mcchain03 = run(m * NUTS() * SerialMC(500:1000))

describe(mcchain03.diagnostics["ndoublings"]) # check # of doublings in NUTS algo


text:
immutable Network

    layers::Vector{Vector{Table}}
    ties::TieBreak

    function Network(inputrange::Int, layersizes::Vector{Int}, inputs::Vector{Int}, ties::TieBreak=negative)

        @assert length(layersizes) >= 1
        @assert layersizes[end] == 1
        @assert all(inputs .<= 16)
        @assert length(layersizes) == length(inputs)
        layerranges = (length(layersizes) == 1) ? [inputrange] : vcat([inputrange], layersizes[1:end-1])
        @assert all([range >= k for (range, k) in zip(layerranges, inputs)])

        layers = Vector{Vector{Table}}(length(layersizes))

        positions = [Set{Int}() for ll in 1:length(layersizes)+1]
        map(n -> push!(positions[end], n), layersizes[end])

        for ll in reverse(1:length(layers))
            layers[ll] = Vector{Table}(layersizes[ll])
            for tt in 1:length(layers[ll])
                if tt in positions[ll+1]
                    randinputs = shuffle(1:layerranges[ll])[1:inputs[ll]]
                    layers[ll][tt] = Table(randinputs, zeros(Int, 2^inputs[ll]), zeros(Int, 2^inputs[ll]))
                    map(n -> push!(positions[ll], n), randinputs)
                end
            end
        end

        new(layers, ties)

    end

end

function (net::Network)(x::Vector{Bool})

    prev = x
    next = Vector{Bool}

    for layer in net.layers
        next = Vector{Bool}(length(layer))
        for tt in 1:length(layer)
            if isassigned(layer, tt)
                next[tt] = layer[tt](prev, net.ties)
            end
        end
        prev = next
    end

    next[end]

end

function (net::Network)(X::Matrix{Bool})

    [net(X[:,ss]) for ss in 1:size(X, 2)]

end

function (net::Network)(X::Matrix{Bool}, y::AbstractArray{Bool})

    @assert size(X, 2) == length(y)

    for layer in net.layers

        for tt in 1:length(layer)
            if isassigned(layer, tt)
                for ss in 1:size(X, 2)
                    layer[tt](X[:,ss], y[ss])
                end
            end
        end

        R = Matrix{Bool}(length(layer), size(X, 2))

        for ss in 1:size(X, 2)
            for tt in 1:length(layer)
                if isassigned(layer, tt)
                    R[tt,ss] = layer[tt](X[:,ss], net.ties)
                end
            end
        end

        X = R

    end

    X[:]

end


text:
#==============================================================================#
# SimpleDB.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/sdb-2009-04-15.normal.json
#==============================================================================#

__precompile__()

module SimpleDB

using AWSCore

"""
    using AWSSDK.SimpleDB.batch_delete_attributes
    batch_delete_attributes([::AWSConfig], arguments::Dict)
    batch_delete_attributes([::AWSConfig]; DomainName=, Item=)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "BatchDeleteAttributes", arguments::Dict)
    sdb([::AWSConfig], "BatchDeleteAttributes", DomainName=, Item=)

# BatchDeleteAttributes Operation

Performs multiple DeleteAttributes operations in a single call, which reduces round trips and latencies. This enables Amazon SimpleDB to optimize requests, which generally yields better throughput.

**Note**
> If you specify BatchDeleteAttributes without attributes or values, all the attributes for the item are deleted.

BatchDeleteAttributes is an idempotent operation; running it multiple times on the same item or attribute doesn't result in an error.

The BatchDeleteAttributes operation succeeds or fails in its entirety. There are no partial deletes. You can execute multiple BatchDeleteAttributes operations and other operations in parallel. However, large numbers of concurrent BatchDeleteAttributes calls can result in Service Unavailable (503) responses.

This operation is vulnerable to exceeding the maximum URL size when making a REST request using the HTTP GET method.

This operation does not support conditions using Expected.X.Name, Expected.X.Value, or Expected.X.Exists.

The following limitations are enforced for this operation:

*   1 MB request size
*   25 item limit per BatchDeleteAttributes operation

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain in which the attributes are being deleted.

## `Item = [[ ... ], ...]` -- *Required*
A list of items on which to perform the operation.
```
 Item = [[
        "ItemName" => <required> ::String,
        "Attribute" =>  [[
            "Name" => <required> ::String,
            "Value" =>  ::String
        ], ...]
    ], ...]
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/BatchDeleteAttributes)
"""
@inline batch_delete_attributes(aws::AWSConfig=default_aws_config(); args...) = batch_delete_attributes(aws, args)

@inline batch_delete_attributes(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "BatchDeleteAttributes", args)

@inline batch_delete_attributes(args) = batch_delete_attributes(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.batch_put_attributes
    batch_put_attributes([::AWSConfig], arguments::Dict)
    batch_put_attributes([::AWSConfig]; DomainName=, Item=)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "BatchPutAttributes", arguments::Dict)
    sdb([::AWSConfig], "BatchPutAttributes", DomainName=, Item=)

# BatchPutAttributes Operation

The `BatchPutAttributes` operation creates or replaces attributes within one or more items. By using this operation, the client can perform multiple [PutAttribute](@ref) operation with a single call. This helps yield savings in round trips and latencies, enabling Amazon SimpleDB to optimize requests and generally produce better throughput.

The client may specify the item name with the `Item.X.ItemName` parameter. The client may specify new attributes using a combination of the `Item.X.Attribute.Y.Name` and `Item.X.Attribute.Y.Value` parameters. The client may specify the first attribute for the first item using the parameters `Item.0.Attribute.0.Name` and `Item.0.Attribute.0.Value`, and for the second attribute for the first item by the parameters `Item.0.Attribute.1.Name` and `Item.0.Attribute.1.Value`, and so on.

Attributes are uniquely identified within an item by their name/value combination. For example, a single item can have the attributes `{ "first_name", "first_value" }` and `{ "first_name", "second_value" }`. However, it cannot have two attribute instances where both the `Item.X.Attribute.Y.Name` and `Item.X.Attribute.Y.Value` are the same.

Optionally, the requester can supply the `Replace` parameter for each individual value. Setting this value to `true` will cause the new attribute values to replace the existing attribute values. For example, if an item `I` has the attributes `{ 'a', '1' }, { 'b', '2'}` and `{ 'b', '3' }` and the requester does a BatchPutAttributes of `{'I', 'b', '4' }` with the Replace parameter set to true, the final attributes of the item will be `{ 'a', '1' }` and `{ 'b', '4' }`, replacing the previous values of the 'b' attribute with the new value.

**Note**
> You cannot specify an empty string as an item or as an attribute name. The `BatchPutAttributes` operation succeeds or fails in its entirety. There are no partial puts.

**Important**
> This operation is vulnerable to exceeding the maximum URL size when making a REST request using the HTTP GET method. This operation does not support conditions using `Expected.X.Name`, `Expected.X.Value`, or `Expected.X.Exists`.

You can execute multiple `BatchPutAttributes` operations and other operations in parallel. However, large numbers of concurrent `BatchPutAttributes` calls can result in Service Unavailable (503) responses.

The following limitations are enforced for this operation:

*   256 attribute name-value pairs per item
*   1 MB request size
*   1 billion attributes per domain
*   10 GB of total user data storage per domain
*   25 item limit per `BatchPutAttributes` operation

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain in which the attributes are being stored.

## `Item = [[ ... ], ...]` -- *Required*
A list of items on which to perform the operation.
```
 Item = [[
        "ItemName" => <required> ::String,
        "Attribute" => <required> [[
            "Name" => <required> ::String,
            "Value" => <required> ::String,
            "Replace" =>  ::Bool
        ], ...]
    ], ...]
```

# Exceptions

`DuplicateItemName`, `InvalidParameterValue`, `MissingParameter`, `NoSuchDomain`, `NumberItemAttributesExceeded`, `NumberDomainAttributesExceeded`, `NumberDomainBytesExceeded`, `NumberSubmittedItemsExceeded` or `NumberSubmittedAttributesExceeded`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/BatchPutAttributes)
"""
@inline batch_put_attributes(aws::AWSConfig=default_aws_config(); args...) = batch_put_attributes(aws, args)

@inline batch_put_attributes(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "BatchPutAttributes", args)

@inline batch_put_attributes(args) = batch_put_attributes(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.create_domain
    create_domain([::AWSConfig], arguments::Dict)
    create_domain([::AWSConfig]; DomainName=)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "CreateDomain", arguments::Dict)
    sdb([::AWSConfig], "CreateDomain", DomainName=)

# CreateDomain Operation

The `CreateDomain` operation creates a new domain. The domain name should be unique among the domains associated with the Access Key ID provided in the request. The `CreateDomain` operation may take 10 or more seconds to complete.

**Note**
> CreateDomain is an idempotent operation; running it multiple times using the same domain name will not result in an error response.

The client can create up to 100 domains per account.

If the client requires additional domains, go to [http://aws.amazon.com/contact-us/simpledb-limit-request/](http://aws.amazon.com/contact-us/simpledb-limit-request/).

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain to create. The name can range between 3 and 255 characters and can contain the following characters: a-z, A-Z, 0-9, '_', '-', and '.'.

# Exceptions

`InvalidParameterValue`, `MissingParameter` or `NumberDomainsExceeded`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/CreateDomain)
"""
@inline create_domain(aws::AWSConfig=default_aws_config(); args...) = create_domain(aws, args)

@inline create_domain(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "CreateDomain", args)

@inline create_domain(args) = create_domain(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.delete_attributes
    delete_attributes([::AWSConfig], arguments::Dict)
    delete_attributes([::AWSConfig]; DomainName=, ItemName=, <keyword arguments>)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "DeleteAttributes", arguments::Dict)
    sdb([::AWSConfig], "DeleteAttributes", DomainName=, ItemName=, <keyword arguments>)

# DeleteAttributes Operation

Deletes one or more attributes associated with an item. If all attributes of the item are deleted, the item is deleted.

**Note**
> If `DeleteAttributes` is called without being passed any attributes or values specified, all the attributes for the item are deleted.

`DeleteAttributes` is an idempotent operation; running it multiple times on the same item or attribute does not result in an error response.

Because Amazon SimpleDB makes multiple copies of item data and uses an eventual consistency update model, performing a [GetAttributes](@ref) or [Select](@ref) operation (read) immediately after a `DeleteAttributes` or [PutAttributes](@ref) operation (write) might not return updated item data.

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain in which to perform the operation.

## `ItemName = ::String` -- *Required*
The name of the item. Similar to rows on a spreadsheet, items represent individual objects that contain one or more value-attribute pairs.

## `Attribute = [[ ... ], ...]`
A list of Attributes. Similar to columns on a spreadsheet, attributes represent categories of data that can be assigned to items.
```
 Attribute = [[
        "Name" => <required> ::String,
        "Value" =>  ::String
    ], ...]
```

## `Expected = [ ... ]`
The update condition which, if specified, determines whether the specified attributes will be deleted or not. The update condition must be satisfied in order for this request to be processed and the attributes to be deleted.
```
 Expected = [
        "Name" =>  ::String,
        "Value" =>  ::String,
        "Exists" =>  ::Bool
    ]
```

# Exceptions

`InvalidParameterValue`, `MissingParameter`, `NoSuchDomain` or `AttributeDoesNotExist`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/DeleteAttributes)
"""
@inline delete_attributes(aws::AWSConfig=default_aws_config(); args...) = delete_attributes(aws, args)

@inline delete_attributes(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "DeleteAttributes", args)

@inline delete_attributes(args) = delete_attributes(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.delete_domain
    delete_domain([::AWSConfig], arguments::Dict)
    delete_domain([::AWSConfig]; DomainName=)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "DeleteDomain", arguments::Dict)
    sdb([::AWSConfig], "DeleteDomain", DomainName=)

# DeleteDomain Operation

The `DeleteDomain` operation deletes a domain. Any items (and their attributes) in the domain are deleted as well. The `DeleteDomain` operation might take 10 or more seconds to complete.

**Note**
> Running `DeleteDomain` on a domain that does not exist or running the function multiple times using the same domain name will not result in an error response.

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain to delete.

# Exceptions

`MissingParameter`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/DeleteDomain)
"""
@inline delete_domain(aws::AWSConfig=default_aws_config(); args...) = delete_domain(aws, args)

@inline delete_domain(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "DeleteDomain", args)

@inline delete_domain(args) = delete_domain(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.domain_metadata
    domain_metadata([::AWSConfig], arguments::Dict)
    domain_metadata([::AWSConfig]; DomainName=)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "DomainMetadata", arguments::Dict)
    sdb([::AWSConfig], "DomainMetadata", DomainName=)

# DomainMetadata Operation

Returns information about the domain, including when the domain was created, the number of items and attributes in the domain, and the size of the attribute names and values.

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain for which to display the metadata of.

# Returns

`DomainMetadataResult`

# Exceptions

`MissingParameter` or `NoSuchDomain`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/DomainMetadata)
"""
@inline domain_metadata(aws::AWSConfig=default_aws_config(); args...) = domain_metadata(aws, args)

@inline domain_metadata(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "DomainMetadata", args)

@inline domain_metadata(args) = domain_metadata(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.get_attributes
    get_attributes([::AWSConfig], arguments::Dict)
    get_attributes([::AWSConfig]; DomainName=, ItemName=, <keyword arguments>)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "GetAttributes", arguments::Dict)
    sdb([::AWSConfig], "GetAttributes", DomainName=, ItemName=, <keyword arguments>)

# GetAttributes Operation

Returns all of the attributes associated with the specified item. Optionally, the attributes returned can be limited to one or more attributes by specifying an attribute name parameter.

If the item does not exist on the replica that was accessed for this operation, an empty set is returned. The system does not return an error as it cannot guarantee the item does not exist on other replicas.

**Note**
> If GetAttributes is called without being passed any attribute names, all the attributes for the item are returned.

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain in which to perform the operation.

## `ItemName = ::String` -- *Required*
The name of the item.

## `AttributeName = [::String, ...]`
The names of the attributes.

## `ConsistentRead = ::Bool`
Determines whether or not strong consistency should be enforced when data is read from SimpleDB. If `true`, any data previously written to SimpleDB will be returned. Otherwise, results will be consistent eventually, and the client may not see data that was written immediately before your read.

# Returns

`GetAttributesResult`

# Exceptions

`InvalidParameterValue`, `MissingParameter` or `NoSuchDomain`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/GetAttributes)
"""
@inline get_attributes(aws::AWSConfig=default_aws_config(); args...) = get_attributes(aws, args)

@inline get_attributes(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "GetAttributes", args)

@inline get_attributes(args) = get_attributes(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.list_domains
    list_domains([::AWSConfig], arguments::Dict)
    list_domains([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "ListDomains", arguments::Dict)
    sdb([::AWSConfig], "ListDomains", <keyword arguments>)

# ListDomains Operation

The `ListDomains` operation lists all domains associated with the Access Key ID. It returns domain names up to the limit set by [MaxNumberOfDomains](#MaxNumberOfDomains). A [NextToken](#NextToken) is returned if there are more than `MaxNumberOfDomains` domains. Calling `ListDomains` successive times with the `NextToken` provided by the operation returns up to `MaxNumberOfDomains` more domain names with each successive operation call.

# Arguments

## `MaxNumberOfDomains = ::Int`
The maximum number of domain names you want returned. The range is 1 to 100. The default setting is 100.

## `NextToken = ::String`
A string informing Amazon SimpleDB where to start the next list of domain names.

# Returns

`ListDomainsResult`

# Exceptions

`InvalidParameterValue` or `InvalidNextToken`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/ListDomains)
"""
@inline list_domains(aws::AWSConfig=default_aws_config(); args...) = list_domains(aws, args)

@inline list_domains(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "ListDomains", args)

@inline list_domains(args) = list_domains(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.put_attributes
    put_attributes([::AWSConfig], arguments::Dict)
    put_attributes([::AWSConfig]; DomainName=, ItemName=, Attribute=, <keyword arguments>)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "PutAttributes", arguments::Dict)
    sdb([::AWSConfig], "PutAttributes", DomainName=, ItemName=, Attribute=, <keyword arguments>)

# PutAttributes Operation

The PutAttributes operation creates or replaces attributes in an item. The client may specify new attributes using a combination of the `Attribute.X.Name` and `Attribute.X.Value` parameters. The client specifies the first attribute by the parameters `Attribute.0.Name` and `Attribute.0.Value`, the second attribute by the parameters `Attribute.1.Name` and `Attribute.1.Value`, and so on.

Attributes are uniquely identified in an item by their name/value combination. For example, a single item can have the attributes `{ "first_name", "first_value" }` and `{ "first_name", second_value" }`. However, it cannot have two attribute instances where both the `Attribute.X.Name` and `Attribute.X.Value` are the same.

Optionally, the requestor can supply the `Replace` parameter for each individual attribute. Setting this value to `true` causes the new attribute value to replace the existing attribute value(s). For example, if an item has the attributes `{ 'a', '1' }`, `{ 'b', '2'}` and `{ 'b', '3' }` and the requestor calls `PutAttributes` using the attributes `{ 'b', '4' }` with the `Replace` parameter set to true, the final attributes of the item are changed to `{ 'a', '1' }` and `{ 'b', '4' }`, which replaces the previous values of the 'b' attribute with the new value.

**Note**
> Using `PutAttributes` to replace attribute values that do not exist will not result in an error response.

You cannot specify an empty string as an attribute name.

Because Amazon SimpleDB makes multiple copies of client data and uses an eventual consistency update model, an immediate [GetAttributes](@ref) or [Select](@ref) operation (read) immediately after a [PutAttributes](@ref) or [DeleteAttributes](@ref) operation (write) might not return the updated data.

The following limitations are enforced for this operation:

*   256 total attribute name-value pairs per item
*   One billion attributes per domain
*   10 GB of total user data storage per domain

# Arguments

## `DomainName = ::String` -- *Required*
The name of the domain in which to perform the operation.

## `ItemName = ::String` -- *Required*
The name of the item.

## `Attribute = [[ ... ], ...]` -- *Required*
The list of attributes.
```
 Attribute = [[
        "Name" => <required> ::String,
        "Value" => <required> ::String,
        "Replace" =>  ::Bool
    ], ...]
```

## `Expected = [ ... ]`
The update condition which, if specified, determines whether the specified attributes will be updated or not. The update condition must be satisfied in order for this request to be processed and the attributes to be updated.
```
 Expected = [
        "Name" =>  ::String,
        "Value" =>  ::String,
        "Exists" =>  ::Bool
    ]
```

# Exceptions

`InvalidParameterValue`, `MissingParameter`, `NoSuchDomain`, `NumberDomainAttributesExceeded`, `NumberDomainBytesExceeded`, `NumberItemAttributesExceeded` or `AttributeDoesNotExist`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/PutAttributes)
"""
@inline put_attributes(aws::AWSConfig=default_aws_config(); args...) = put_attributes(aws, args)

@inline put_attributes(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "PutAttributes", args)

@inline put_attributes(args) = put_attributes(default_aws_config(), args)

"""
    using AWSSDK.SimpleDB.select
    select([::AWSConfig], arguments::Dict)
    select([::AWSConfig]; SelectExpression=, <keyword arguments>)

    using AWSCore.Services.sdb
    sdb([::AWSConfig], "Select", arguments::Dict)
    sdb([::AWSConfig], "Select", SelectExpression=, <keyword arguments>)

# Select Operation

The `Select` operation returns a set of attributes for `ItemNames` that match the select expression. `Select` is similar to the standard SQL SELECT statement.

The total size of the response cannot exceed 1 MB in total size. Amazon SimpleDB automatically adjusts the number of items returned per page to enforce this limit. For example, if the client asks to retrieve 2500 items, but each individual item is 10 kB in size, the system returns 100 items and an appropriate `NextToken` so the client can access the next page of results.

For information on how to construct select expressions, see Using Select to Create Amazon SimpleDB Queries in the Developer Guide.

# Arguments

## `SelectExpression = ::String` -- *Required*
The expression used to query the domain.

## `NextToken = ::String`
A string informing Amazon SimpleDB where to start the next list of `ItemNames`.

## `ConsistentRead = ::Bool`
Determines whether or not strong consistency should be enforced when data is read from SimpleDB. If `true`, any data previously written to SimpleDB will be returned. Otherwise, results will be consistent eventually, and the client may not see data that was written immediately before your read.

# Returns

`SelectResult`

# Exceptions

`InvalidParameterValue`, `InvalidNextToken`, `InvalidNumberPredicates`, `InvalidNumberValueTests`, `InvalidQueryExpression`, `MissingParameter`, `NoSuchDomain`, `RequestTimeout` or `TooManyRequestedAttributes`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/sdb-2009-04-15/Select)
"""
@inline select(aws::AWSConfig=default_aws_config(); args...) = select(aws, args)

@inline select(aws::AWSConfig, args) = AWSCore.Services.sdb(aws, "Select", args)

@inline select(args) = select(default_aws_config(), args)

end # module SimpleDB

#==============================================================================#
# End of file
#==============================================================================#


text:
using EMIRT
using EMIRT.Images
using Test

@testset "test image" begin
    img = rand(UInt8, 128,128,256)

    println("test image normalization ...")
    println("serial version: ")
    @time img3 = Images.normalize2d_serial(img)
    println("parallel version: ")
    @time img2 = Images.normalize2d(img)
    @test all( img2 .== img3 )

    println("test image mask")
    @time Images.image2mask(img; threshold = UInt8(128), sizeThreshold = UInt32(400))
end # end of test set


text:
@testset "Twists of elliptic curves" begin
  K = GF(2,2)
  a = gen(K)
  E = elliptic_curve_from_j_invariant(K(0))
  L = @inferred twists(E)
  twistlist = [EllipticCurve(K, [0, 0, 1, 0, 0 ]),
               EllipticCurve(K, [0, 0, 1, a, 0 ]),
               EllipticCurve(K, [0, 0, 1, 0, a ]),
               EllipticCurve(K, [0, 0, a, 0, 0 ]),
               EllipticCurve(K, [0, 0, a^2, 0, 0 ]),
               EllipticCurve(K, [0, 0, a, 0, 1 ]),
               EllipticCurve(K, [0, 0, a^2, 0, a^2])]
  @test L == twistlist
  for tw in L
    @test @inferred is_twist(E, tw)
  end

  E = elliptic_curve_from_j_invariant(K(a+1))
  F = @inferred quadratic_twist(E)
  @test is_twist(E, F)

  K = GF(3,2)
  a = gen(K)
  E = elliptic_curve_from_j_invariant(K(0))
  L = @inferred twists(E)

  twistlist = [EllipticCurve(K, [ 0, 0, 0, 2, 0 ]),
               EllipticCurve(K, [ 0, 0, 0, a^5, 0 ]),
               EllipticCurve(K, [ 0, 0, 0, a^6, 0 ]),
               EllipticCurve(K, [ 0, 0, 0, a^7, 0 ]),
               EllipticCurve(K, [ 0, 0, 0, 2, 1 ]),
               EllipticCurve(K, [ 0, 0, 0, a^6, a^3 ])]
  @test L == twistlist

  E = elliptic_curve_from_j_invariant(K(a+1))
  F = @inferred quadratic_twist(E)
  @test is_twist(E, F)

  K = GF(7,2)
  a = gen(K)
  E = elliptic_curve_from_j_invariant(K(0))
  L = @inferred twists(E)

  twistlist = [EllipticCurve(K, [ 0, 0, 0, 0, 2 ]),
               EllipticCurve(K, [ 0, 0, 0, 0, a^17 ]),
               EllipticCurve(K, [ 0, 0, 0, 0, a^18 ]),
               EllipticCurve(K, [ 0, 0, 0, 0, a^19 ]),
               EllipticCurve(K, [ 0, 0, 0, 0, a^20 ]),
               EllipticCurve(K, [ 0, 0, 0, 0, a^21 ])]
  @test L == twistlist

  E = elliptic_curve_from_j_invariant(K(1728))
  L = @inferred twists(E)

  twistlist = [EllipticCurve(K, [ 0, 0, 0, 1, 0]),
               EllipticCurve(K, [ 0, 0, 0, a, 0]),
               EllipticCurve(K, [ 0, 0, 0, a^2,0]),
               EllipticCurve(K, [ 0, 0, 0, a^3, 0])]
  @test L == twistlist

  E = EllipticCurve([1,2,3,4,5])
  F = @inferred quadratic_twist(E, 5)
  @test !is_isomorphic(E, F)
  @test is_twist(E, F)

end


text:
# It is not possible to test the Performance structure, since the
# output depends on the instant the columns is created and when it is
# called. We can test the SI prefix shifting however.

using Unitful
import SolverTraces: shift_prefix, SI_prefix_convert

@testset "SI prefices" begin
    @test shift_prefix(u"MHz", 6) == u"THz"
    @test shift_prefix(u"MHz", 3) == u"GHz"
    @test shift_prefix(u"MHz", -3) == u"kHz"
    @test shift_prefix(u"MHz", -6) == u"Hz"
    @test shift_prefix(u"MHz", -9) == u"mHz"
    @test_throws ArgumentError shift_prefix(u"MHz", 2)

    @test SI_prefix_convert(1.0u"nm") == 1.0u"nm"
    @test 0.10u"nm" ≠ 100.0u"pm"
    @test SI_prefix_convert(0.10u"nm") == 100.0u"pm"
    @test SI_prefix_convert(100.0u"nm") == 100.0u"nm"
    @test SI_prefix_convert(999.0u"nm") == 999.0u"nm"
end


text:
include("../../src/required.jl")
include("lib/quantiles.jl")

## Conditional arguments
# outcome = "raw" # desired outcome (required)
# outcome = "bart"
# save_quantile_figures = true # should quantile figures be overwritten (optional)

# Make sure "outcome" is defined
if !(@isdefined outcome)
    @warn "'outcome' not defined, must be either 'raw', 'bio', 'rf, or 'bart'"
elseif !(outcome in ("raw", "bio", "rf", "bart"))
    @warn "'outcome' invalid, must be either 'raw', 'bio', 'rf', or 'bart'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Load presence-absence data for all species
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name
distributions = [
    geotiff(
        SimpleSDMPredictor, joinpath("data", "proc", "distributions_$(outcome).tif"), i
    ) for i in eachindex(spenames)
]

## Y matrix (site-by-species community data table)

# Create matrix Y
Y = Ymatrix(distributions; transform=false)

# Create matrix of observed sites only
Yobs = _Yobs(Y)

## Richness

# Get richness
richness = calculate_richness(Y, distributions[1])

# Plot richness quantiles
richness_qplot = plot_layer(
    quantiles(richness);
    c=:viridis,
    # title = "Richness quantiles ($outcome distributions)",
    colorbar_title="Species richness (quantiles)",
    size=(650, 400),
)

## LCBD

# Get LCBD values
lcbd_rel = calculate_lcbd(Y, distributions[1]; transform=true, relative=true)

# Plot quantile scores
lcbdtr_qplot = plot_layer(
    quantiles(lcbd_rel);
    c=:viridis,
    # title = "LCBD quantiles ($(outcome) distributions, hellinger transformed)",
    colorbar_title="LCBD value (quantiles)",
    size=(650, 400),
)

## Export figures

# save_quantile_figures = true # should quantile figures be overwritten (optional)
if (@isdefined save_quantile_figures) && save_quantile_figures == true
    @info "Quantile figures saved ($(outcome))"
    savefig(
        plot(richness_qplot; dpi=200),
        joinpath("archive", "fig", "quantiles", "04_$(outcome)_richness_quantiles.png"),
    )
    savefig(
        plot(lcbdtr_qplot; dpi=200),
        joinpath("archive", "fig", "quantiles", "04_$(outcome)_lcbd.png"),
    )
else
    @info "Quantile figures not saved ($(outcome) richness)"
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
hurst_indes(Ra)

pain_index(Ra)

smoothing_index(Ra) 0.9825174 # R's PerformanceAnalytics result

ulcer_index(Ra) => 1.128081 # R's PerformanceAnalytics result


text:
using SymPy
using Compat.Test

include("tests.jl")
include("test-math.jl")
include("test-matrix.jl")
include("test-ode.jl")
include("test-logical.jl")
include("test-specialfuncs.jl")  ## XXX NEEDS WORK
include("test-permutations.jl")
include("test-physics.jl")


text:
using Documenter, Literate, HydroPowerModels

examples_dir = joinpath(dirname(dirname(@__FILE__)), "examples/HydroValleys")
docs_dir = dirname(@__FILE__)
testcases_dir = joinpath(dirname(dirname(@__FILE__)), "testcases")

plot_bool = true

makedocs(
    modules = [HydroPowerModels],
    doctest  = false,
    clean    = true,
    format   = Documenter.HTML(),
    sitename = "HydroPowerModels.jl",
    authors = "Andrew Rosemberg",
    pages = [
        "Home"      => "index.md",
        "Manual"    => "getstarted.md",
        "Case data" => "inputfiles.md",
        "Examples"  => "examples/cases.md",
        "Reference" => "apireference.md"
    ]
)

include("make_examples.jl")

deploydocs(
    repo   = "github.com/andrewrosemberg/HydroPowerModels.jl.git",
)


text:
using Fermi.HartreeFock

import Base: show

export RMP2

"""
    Fermi.MollerPlesset.RMP2Algorithm

Abstract type for RMP2 implementations.
"""
abstract type RMP2Algorithm end

"""
    Fermi.MollerPlesset.RMP2

Wave function object for Second Order Restricted Moller-Plesset methods.

# High Level Interface
Run a RMP2 computation and return the RMP2 object:
```
julia> @energy rmp2
```
Equivalent to
```
julia> Fermi.MollerPlesset.RMP2()
```
This function calls a constructor that runs a MP2 computation based on the options found in `Fermi.Options`.

# Fields

| Name   |   Description     |
|--------|---------------------|
| `correlation` |   Computed RMP2 correlation energy |
| `energy`   |   Total wave function energy (Reference energy + Correlation energy)      |

# Relevant options

These options can be set with `@set <option> <value>`

| Option         | What it does                      | Type      | choices [default]     |
|----------------|-----------------------------------|-----------|-----------------------|
| `mp2_alg`      | Picks MP2 algorithm               | `Int`     | [1]                   |
| `basis`        | What basis set to use             | `String`  | ["sto-3g"]            |
| `df`           | Whether to use density fitting    | `Bool`    | `true` [`false`]      |
| `rifit`        | What aux. basis set to use for RI | `String`  | ["auto"]              |
| `drop_occ`    | Number of occupied electrons to be dropped | `Int`  | [0]              |
| `drop_vir`    | Number of virtual electrons to be dropped | `Int`  | [0]              |
"""
struct RMP2{T} <: AbstractMPWavefunction
    correlation::T
    energy::T
end

"""
    Fermi.MollerPlesset.get_rmp2_alg()

Returns a singleton type corresponding to a RMP2 implementation based on the options.
"""
function get_rmp2_alg(N::Int = Options.get("mp2_alg"))
    try
        return get_rmp2_alg(Val(N))
    catch MethodError
        throw(FermiException("implementation number $N not available for RMP2."))
    end
end

# For each implementation a singleton type must be create
struct RMP2a <: RMP2Algorithm end
include("RMP2a.jl")
# And a number is assigned to the implementation
get_rmp2_alg(x::Val{1}) = RMP2a()

struct RMP2b <: RMP2Algorithm end
include("RMP2b.jl")
# And a number is assigned to the implementation
get_rmp2_alg(x::Val{2}) = RMP2b()

function RMP2(x...)
    if !any(i-> i isa RMP2Algorithm, x)
        RMP2(x..., get_rmp2_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RMP2 method: $args"))
    end
end

## MISCELLANEOUS
# Pretty printing
function string_repr(X::RMP2)
    out = ""
    out = out*" ⇒ Fermi Restricted MP2 Wave function\n"
    out = out*" ⋅ Correlation Energy:     $(X.correlation)\n"
    out = out*" ⋅ Total Energy:           $(X.energy)"
    return out
end

function show(io::IO, ::MIME"text/plain", X::RMP2)
    print(io, string_repr(X))
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
module CellFishing

using Arpack: svds
using LinearAlgebra: lu!, qr!, svd
using Random: MersenneTwister
using Serialization: serialize, deserialize
using SparseArrays: SparseMatrixCSC
using Statistics: mean, std
using Mmap: Mmap
using Blosc: libblosc
using Distributions: NegativeBinomial, logcdf, logccdf

include("svd.jl")
include("bitvectors.jl")
include("HammingIndexes.jl")
include("cmatrix.jl")
include("ematrix.jl")
include("preprocessor.jl")
include("features.jl")
include("index.jl")
include("search.jl")
include("degenes.jl")

using .HammingIndexes: bitsof

end  # module CellFishing