text:
# This file is designed to test the visibility interpolation algorithm.
using Test
import PyPlot
import PyPlot.plt
using LaTeXStrings

using DiskJockey.constants
using DiskJockey.gridding
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gauss

# Because of the flipped nature of the sky (but not flipped nature of the UV plane)
# there are some tricky conventions about how to pack the array.

# Test to see if we get the convolutional interpolation correct by using a 2D
# Elliptical Gaussian.
# Realistic Gaussian will have scale dimensions (fatter in x direction)
const mu_RA = 2.0 # [arcsec]
const mu_DEC = -0.5 # [arcsec]
const s_x = 1.2 # [arcsec]
const s_y = 1.0 # [arcsec]
const rho = 0.5
const p0 = [mu_RA, mu_DEC, s_x, s_y, rho] # [arcsec]
const p_center = [0.0, 0.0, s_x, s_y, rho] # [arcsec]

# If the synthesized beam is 0.3", to oversample this by a factor of 10 would require something like 512x512 pixels
nx = 128
ny = 128

# full span of the image
ra = fftspace(10.0, nx) # [arcsec]
dec = fftspace(10.0, ny) # [arcsec]

# convert ra and dec in [arcsec] to radians, and then take the sin to convert to ll, mm
ll = sin.(ra * arcsec) # direction cosines
mm = sin.(dec * arcsec)

# The natural, shifted Gaussian image
img = imageGauss(ll, mm, p0, 1)

# A centered Gaussian image which we will then shift via the phase theorem
img_plain1 = imageGauss(ll, mm, p_center, 1)
img_plain2 = imageGauss(ll, mm, p_center, 1)
img_plain3 = imageGauss(ll, mm, p_center, 1)

lam0 = lam0s["12CO2-1"] # [microns] 12CO 2-1

skim = SkyImage(img, ra, dec, lam0)

skim_plain1 = SkyImage(img_plain1, ra, dec, lam0)
skim_plain2 = SkyImage(img_plain2, ra, dec, lam0)
skim_plain3 = SkyImage(img_plain3, ra, dec, lam0)

# Apply a centered correction function on the offset image and Fourier transform the image
corrfun!(skim, 0.0, 0.0)
shift_fft = transform(skim)

# Now, create images that we will use as test cases for doing the interpolation before or after
# sampling

# Apply a centered correction function on the centered image (to phase shift downsampled
# visibilities later)
# This is the recommended approach.
corrfun!(skim_plain1, 0, 0)
vis_fft_center = transform(skim_plain1)

# Apply an offset correction function on the centered image (then  phase shift the full visibilities)
corrfun!(skim_plain2, mu_RA, mu_DEC)
vis_fft_shift = transform(skim_plain2)
phase_shift!(vis_fft_shift, mu_RA, mu_DEC)

# Apply an offset correction function on the centered image, but later we will phase shift the visibilities *after* they have been downsampled.
# NOTE, this DOES NOT WORK, and is simply included here for full testing.
corrfun!(skim_plain3, mu_RA, mu_DEC)
vis_fft_offset = transform(skim_plain3)

# This right here is not actually a fair comparison, because we haven't yet *interpolated* the FFTed values using the gridding convolution functions.

# Complex subtraction
# println("Maximum FFT discrepancy: ", maximum(abs(plain_fft.VV - vis_analytic)))
# println("Maximum FFT discrepancy: ", maximum(abs(shift_fft.VV - vis_analytic)))

# Return a normalized instance that is symmetric about 0
function scale(data)
    s = maximum(abs.(data))
    return norm = plt.Normalize(vmin = -s, vmax = s, clip = false)
end

# Basic plot of the shifted image

@testset "plot image array" begin
    # Because the sky convention is different than the way the SkyImage is stored,
    # we need to flip the array for plotting
    fig, ax = plt.subplots(nrows = 2, figsize = (5, 8))

    ext = (skim.ra[end], skim.ra[1], skim.dec[1], skim.dec[end])
    ax[1].imshow(reverse(skim.data[:,:,1], dims = 2), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
    ax[1].contour(reverse(skim.data[:,:,1], dims = 2), origin = "lower", extent = ext)
    ax[1].set_title("Sky Projection")
    ax[1].set_xlabel(L"$\alpha$ [arcsec]")
    ax[1].set_ylabel(L"$\delta$ [arcsec]")

    ext = (ll[1], ll[end], mm[1], mm[end])
    ax[2].imshow(skim.data[:,:,1], interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
    ax[2].contour(skim.data[:,:,1], origin = "lower", extent = ext)
    ax[2].set_title("Raw Array")
    ax[2].set_xlabel(L"$ll$")
    ax[2].set_ylabel(L"$mm$")

    fig.subplots_adjust(left = 0.15, right = 0.85, hspace = 0.25)
    plt.savefig(joinpath(outputdir, "gaussian_img_array.png"))
end

# Basic plot of the *centered* Gaussian image.
fig, ax = plt.subplots(nrows = 1, figsize = (5, 5))
# Real, analytic Gaussian
ext = (skim.ra[end], skim.ra[1], skim.dec[1], skim.dec[end])
aximg = ax.imshow(reverse(skim_plain1.data[:,:,1], dims = 2), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext) # , norm = scale(img))
ax.set_title("image")
ax.set_xlabel(L"$\alpha$ [arcsec]")
ax.set_ylabel(L"$\delta$ [arcsec]")
# [left, bottom, width, height]
cax = fig.add_axes([0.84, 0.25, 0.03, 0.45])
cb = fig.colorbar(aximg, cax = cax)

fig.subplots_adjust(left = 0.15, right = 0.85, hspace = 0.25)
plt.savefig(joinpath(outputdir, "gaussian_img_center.png"))

# It's nice to see these plots and what they look like, but it doesn't really make sense
# to compare them, unless we are comparing *interpolated* visibilities to the analytic formula.
# Just comparing the image is WRONG.

# THIS (following) IS THE ONLY THING WE SHOULD CARE ABOUT.

# Next, choose some (u,v) points within the bounds and see how the
# interpolated values compare to what the FTGauss would return.

nu = length(shift_fft.uu)
zer = zeros(nu)

function plot_1d(analytic, approx, fname)
    fig, ax = plt.subplots(nrows = 2, figsize = (5, 8))
    ax[1].plot(shift_fft.uu, zer, ".k", label = "Grid spacing")
    ax[1].plot(uu, real(approx), "ob", label = "Interp")
    ax[1].plot(uu, real(analytic), ".r", label = "Analytic")
    # ax[1][:plot](vis_fft.uu, real(analytic_u), "or", label="Analytic")
    ax[1].set_xlim(-100, 100)
    ax[1].set_title("Real")
    ax[1].set_xlabel(L"u [k $\lambda$]")
    ax[1].legend()

    ax[2].plot(shift_fft.uu, zer, ".k", label = "Grid spacing")
    ax[2].plot(uu, imag(approx), "ob", label = "Interp")
    ax[2].plot(uu, imag(analytic), ".r", label = "Analytic")
    # ax[2][:plot](vis_fft.uu, imag(analytic_u), "or", label="Analytic")
    ax[2].set_xlim(-100, 100)
    ax[2].set_title("Imag")
    ax[2].set_xlabel(L"u [k $\lambda$]")

    plt.savefig(joinpath(outputdir, fname))
end

n = 100
uu = LinRange(-100, 100, n) # [kλ]
approx = Array{ComplexF64}(undef, n)
analytic = Array{ComplexF64}(undef, n)

# First, let's see how the interpolated points, with *no shift*, correspond to the analytic form
# with no shift.

for i = 1:n
    u = uu[i]
    v = 0.0
    # println("Interpolating at $u, $v")
    approx[i] = interpolate_uv(u, v, vis_fft_center)
    analytic[i] = FTGauss(u, v, p_center, 1)
end

plot_1d(analytic, approx, "interpolation.png")

# Next, let's see how the image shifted, then interpolated points compare with the
# analytic form

for i = 1:n
    u = uu[i]
    v = 0.0
    # println("Interpolating at $u, $v")
    approx[i] = interpolate_uv(u, v, shift_fft)
    analytic[i] = FTGauss(u, v, p0, 1)
end

plot_1d(analytic, approx, "interpolation_raw_image_offset.png")

# Next,  let's see how the FFT'ed, then interpolated, then shifted points compare with the
# analytic form
mu = Float64[mu_RA, mu_DEC] * arcsec # [radians]

for i = 1:n
    u = uu[i]
    v = 0.0

    # Now shift the interpolated points according to the phase theorem
    R = Float64[u, v] * 1e3 # [λ]
    # Not actually in polar phase form
    shift = exp(-2pi * 1.0im * (R' * mu)[1])

    approx[i] = shift * interpolate_uv(u, v, vis_fft_center)
    analytic[i] = FTGauss(u, v, p0, 1)
end

plot_1d(analytic, approx, "interpolation_shift_visibilities.png")

# Finally, let's check the way we've been doing things, which is to phase shift the full visibilities first, then downsample.

for i = 1:n
    u = uu[i]
    v = 0.0
    # println("Interpolating at $u, $v")
    approx[i] = interpolate_uv(u, v, vis_fft_shift)
    analytic[i] = FTGauss(u, v, p0, 1)
end

plot_1d(analytic, approx, "interpolation_shift_image.png")

# As a final, final check that the resampling is in fact linear, let's use the centered image, offset corrfun, downsample this, then apply the phase shift at the end. This is where all my confusion started. If we were just able to do the previous two tests OK, it seems like this *shouldn't work*, because if in fact everything is linear, then the only difference is that this has the corrfun applied w/ a shift.

for i = 1:n
    u = uu[i]
    v = 0.0

    # Now shift the interpolated points according to the phase theorem
    R = Float64[u, v] * 1e3 # [λ]
    # Not actually in polar phase form
    shift = exp(-2pi * 1.0im * (R' * mu)[1])

    approx[i] = shift * interpolate_uv(u, v, vis_fft_offset)
    analytic[i] = FTGauss(u, v, p0, 1)
end

plot_1d(analytic, approx, "interpolation_linear_operations.png")

# Ok, this *does not work*. What is the reason? Is the resampling operation not linear apparently?
# I don't think it's that, rather, I think that in order for the gridding correction function/ gridding convolution function pair to work properly, and properly cancel each other out, then they need to correspond to the same phase center. So, either you do both fixed to 0,0, or you do both offset to some value. Once the interpolation is done, you can phase shift as much as you want, but don't mix order of operations.

# Now, we can check things on an actual 2D image if we so please.

# Create analytic function on a smaller grid
n = 256
uu = LinRange(-150, 150, n)
vv = LinRange(-150, 150, n)

function plot_2d(analytic::Matrix{ComplexF64}, approx::Matrix{ComplexF64}, fname::AbstractString)

    println("Analyzing $fname")

    # Plot the reals first

    fig, ax = plt.subplots(nrows = 3, figsize = (5, 11))

    ext = (uu[1], uu[end], vv[1], vv[end])

    axan = ax[1].imshow(real(analytic), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
    ax[1].set_title("Analytic FT")
    ax[1].set_xlabel(L"uu [k$\lambda$]")
    ax[1].set_ylabel(L"vv [k$\lambda$]")

    cax = fig.add_axes([0.84, 0.70, 0.03, 0.25])
    cb = fig.colorbar(axan, cax = cax)

    axfft = ax[2].imshow(real(approx), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
    ax[2].set_title("Interpolated Visibilites from FFT")
    ax[2].set_xlabel(L"uu [k$\lambda$]")
    ax[2].set_ylabel(L"vv [k$\lambda$]")

    cax = fig.add_axes([0.84, 0.40, 0.03, 0.25])
    cb = fig.colorbar(axfft, cax = cax)

    diff = real(analytic - approx)

    # Measure the max value in the analytic array, just as an opportunity to set scale.
    println("Maximum analytic value ", maximum(abs, analytic))

    # Then, measure the total Mod squared difference across the whole image
    println("Total sqrt(modsquared) difference across image ", sqrt(sum(abs2, diff)))

    axdif = ax[3].imshow(diff, interpolation = "none", origin = "lower", cmap = plt.get_cmap("bwr"), extent = ext, norm = scale(diff))
    ax[3].set_title("Difference")
    ax[3].set_xlabel(L"uu [k$\lambda$]")
    ax[3].set_ylabel(L"vv [k$\lambda$]")

    cax = fig.add_axes([0.84, 0.10, 0.03, 0.25])
    cb = fig.colorbar(axdif, cax = cax)

    fig.subplots_adjust(hspace = 0.25, top = 0.97, bottom = 0.06)

    plt.savefig(joinpath(outputdir, fname * "_real.png"))

    fig, ax = plt.subplots(nrows = 3, figsize = (5, 11))

    axan = ax[1].imshow(imag(analytic), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
    ax[1].set_title("Analytic FT")
    ax[1].set_xlabel(L"uu [k$\lambda$]")
    ax[1].set_ylabel(L"vv [k$\lambda$]")

    cax = fig.add_axes([0.84, 0.70, 0.03, 0.25])
    cb = fig.colorbar(axan, cax = cax)

    axfft = ax[2].imshow(imag(approx), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
    ax[2].set_title("Interpolated Visibilites from FFT")
    ax[2].set_xlabel(L"uu [k$\lambda$]")
    ax[2].set_ylabel(L"vv [k$\lambda$]")

    cax = fig.add_axes([0.84, 0.40, 0.03, 0.25])
    cb = fig.colorbar(axfft, cax = cax)

    diff = imag(analytic - approx)

    axdif = ax[3].imshow(diff, interpolation = "none", origin = "lower", cmap = plt.get_cmap("bwr"), extent = ext, norm = scale(diff))
    ax[3].set_title("Difference")
    ax[3].set_xlabel(L"uu [k$\lambda$]")
    ax[3].set_ylabel(L"vv [k$\lambda$]")

    cax = fig.add_axes([0.84, 0.10, 0.03, 0.25])
    cb = fig.colorbar(axdif, cax = cax)

    fig.subplots_adjust(hspace = 0.25, top = 0.97, bottom = 0.06)

    plt.savefig(joinpath(outputdir, fname * "_imag.png"))

    println()

end

# First, let's see how the interpolated points, with *no shift*, correspond to the analytic form
# with no shift.
approx = Array{ComplexF64}(undef, n, n)
analytic = FTGauss(uu, vv, p_center, 1) # can take in full arrays

for i = 1:n
    for j = 1:n
        approx[j,i] = interpolate_uv(uu[i], vv[j], vis_fft_center)
    end
end

plot_2d(analytic, approx, "2D_interpolation")

# Next, let's see how the image shifted, then interpolated points compare with the
# analytic form

analytic = FTGauss(uu, vv, p0, 1) # can take in full arrays

for i = 1:n
    for j = 1:n
        approx[j,i] = interpolate_uv(uu[i], vv[j], shift_fft)
    end
end

plot_2d(analytic, approx, "2D_interpolation_raw_image_offset")

# Next,  let's see how the FFT'ed, then interpolated, then shifted points compare with the
# analytic form

for i = 1:n
    for j = 1:n
        R = Float64[uu[i], vv[j]] * 1e3 # [λ]
        # Not actually in polar phase form
        shift = exp(-2pi * 1.0im * (R' * mu)[1])

        approx[j,i] = shift * interpolate_uv(uu[i], vv[j], vis_fft_center)
    end
end

plot_2d(analytic, approx, "2D_interpolation_shift_visibilities")

# Finally, let's check the way we've been doing things, which is to phase shift the full visibilities first, then downsample.

for i = 1:n
    for j = 1:n
        approx[j,i] = interpolate_uv(uu[i], vv[j], vis_fft_shift)
    end
end

plot_2d(analytic, approx, "2D_interpolation_shift_image")

fig, ax = plt.subplots(nrows = 2, figsize = (5, 8))

axan = ax[1].imshow(abs.(analytic), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
ax[1].set_title("Amplitude [Analytic]")
ax[1].set_xlabel(L"uu [k$\lambda$]")
ax[1].set_ylabel(L"vv [k$\lambda$]")

cax = fig.add_axes([0.84, 0.70, 0.03, 0.25])
cb = fig.colorbar(axan, cax = cax)

axfft = ax[2].imshow(angle.(analytic), interpolation = "none", origin = "lower", cmap = plt.get_cmap("Greys"), extent = ext)
ax[2].set_title("Phase [Analytic]")
ax[2].set_xlabel(L"uu [k$\lambda$]")
ax[2].set_ylabel(L"vv [k$\lambda$]")

cax = fig.add_axes([0.84, 0.20, 0.03, 0.25])
cb = fig.colorbar(axfft, cax = cax)

fig.subplots_adjust(hspace = 0.25, top = 0.95, bottom = 0.1)

plt.savefig(joinpath(outputdir, "interpolation_phase.png"))


text:
# ------------------------------------------------------------------------------------------
#
#
# # JuliaDB is an analytical database that...
#
# - **Loads multi-file datasets**
# - **Sorts data by index variables for fast filter, aggregation, sort, and join
# operations**
# - **Compiles queries** (it's Julia all the way down)
# - **Stores ANY data type** (again, it's Julia)
# - **Saves tables to disk for fast reloading**
# - **Uses Julia's built-in parallel computing features to fully utilize any machine or
# cluster**
# - **Integrates with [OnlineStats](https://github.com/joshday/OnlineStats.jl) for big (or
# small) data analytics**
#
# # Helpful links
#
# For additional documentation, please refer to
#
# - [JuliaDB Docs](http://juliadb.org/latest/api/index.html)
# - [OnlineStats Docs](http://joshday.github.io/OnlineStats.jl/stable/)
#
# # Getting Started
#
# In this notebook, we'll introduce **JuliaDB**'s two main data structures (`Table` and
# `NDSparse`) by
#
# 1. Creating tables from vectors
# 1. Accessing data from `Table` and `NDSparse`
# 1. Loading tables from CSVs
# 1. Saving tables into binary format
# 1. Reloading a saved table
# 1. Using selectors
# ------------------------------------------------------------------------------------------

using JuliaDB

# ------------------------------------------------------------------------------------------
# # Creating tables from Julia Vectors
#
# - First we'll make some data vectors
# ------------------------------------------------------------------------------------------

x = [false, true, false, true]
y = ['B', 'B', 'A', 'A']
z = [.1, .3, .2, .4]
x, y, z

# ------------------------------------------------------------------------------------------
# ## `Table`: Tabular Data Sorted by Primary Key(s)
#
# - We can create a table in **JuliaDB** using the `table` function.
# ------------------------------------------------------------------------------------------

# x, y, z become columns of the table, with default numbering 1, 2, 3.

t1 = table(x, y, z)

# The keyword argument `names` lets us label the columns

t2 = table(x,  y, z, names = [:x, :y, :z])

# ------------------------------------------------------------------------------------------
# ### Sorting by Index Variables
# - Furthermore, we can **sort** a table by a primary key (or keys), which can be set with
# the keyword argument `pkey`.
# - Below, the rows of `t3` are sorted first by the values in column `x` and second by the
# values in column `y`.
# ------------------------------------------------------------------------------------------

t3 = table(x, y, z, names = [:x, :y, :z], pkey = (:x, :y))

# ------------------------------------------------------------------------------------------
# ### Tuple vs. NamedTuple
#
# - We can also build a `table` from a `NamedTuple` of vectors.
# - A `NamedTuple` is created with the `@NT` macro.  Values can be accessed by position or
# name.
# ------------------------------------------------------------------------------------------

# Tuple

a = ("John", "Doe")

@show a[1];

# NamedTuple

b = @NT(firstname = "John", lastname = "Doe")

@show b[1]
@show b[:firstname]
@show b.firstname;

# Column names are taken from the NamedTuple

t4 = table(@NT(x=x, y=y, z=z), pkey = [:x, :y])

# ------------------------------------------------------------------------------------------
# ### `NDSparse`: N-dimensional Sparse Array with Arbitrary Indexes
#
# - Compare the following two data structures:
# ------------------------------------------------------------------------------------------

# maps (tuple of integers) -> value

sparse(reshape(z, 2, 2))

# maps (tuple of arbitrary index types) -> value

ndsparse((x, y), z)

# ------------------------------------------------------------------------------------------
# - Like `table`, we can use `ndsparse` with NamedTuples:
# ------------------------------------------------------------------------------------------

nd2 = ndsparse(@NT(x=x, y=y), z)

# ------------------------------------------------------------------------------------------
# ---
#
# # Accessing data from `Table` and `NDSparse`
#
# ## Index into `Table`
#
# - In the last section, we created the `Table`, `t3`:
# ------------------------------------------------------------------------------------------

t3

# ------------------------------------------------------------------------------------------
# - We can get the first row of `t3` by indexing:
# ------------------------------------------------------------------------------------------

t3[1]

# ------------------------------------------------------------------------------------------
# - We can get an individual element in the row by using the column name:
# ------------------------------------------------------------------------------------------

t3[1].z

# ------------------------------------------------------------------------------------------
# ## Index into `NDSparse`
#
# - Since `NDSparse` acts like an array, accessing data is slightly different.
# - We must look up a value by using the index columns:
# ------------------------------------------------------------------------------------------

nd2[false, 'A']

# ------------------------------------------------------------------------------------------
# - Since `NDSparse` acts like an array, we can look up all the rows where `x` is `false`:
# ------------------------------------------------------------------------------------------

nd2[false, :]

# ------------------------------------------------------------------------------------------
# ## Iteration
#
# - `Table` and `NDSparse` also differ in how they iterate over the data.
# ------------------------------------------------------------------------------------------

# Table: iterate over Tuples/NamedTuples of rows

for row in t3
    println(row)
end

# NDSparse: iterate over values

for item in nd2
    println(item)
end

# ------------------------------------------------------------------------------------------
# ---
#
# # Loading tables from CSVs
#
# - Now we'll load data from multiple CSV files into a table.
# - Our example data is 8 different stocks' OHLC data (each stock in a separate file).
#     - We can start by looking at all the available files in the `stocksample` directory.
#     - `;` allows us to run shell commands from Jupyter or a Julia REPL.
#
# ------------------------------------------------------------------------------------------

;ls stocksample

# ------------------------------------------------------------------------------------------
# - Each txt file has the structure:
# ------------------------------------------------------------------------------------------

;head stocksample/aapl.us.txt

# ------------------------------------------------------------------------------------------
# - We can now use the `loadtable` function to load all files in the `stocksample` directory
# into one table:
# ------------------------------------------------------------------------------------------

stocks = loadtable("stocksample")

# ------------------------------------------------------------------------------------------
# - All data from the 8 files is now in one table, but we don't have the stock ticker
# information!
# - `loadtable` has a `filenamecol` option that will add a column for us, which we will call
# `:Ticker`:
# ------------------------------------------------------------------------------------------

stocks = loadtable("stocksample"; filenamecol = :Ticker)

# ------------------------------------------------------------------------------------------
# - Furthermore, we can use the `indexcols` to specify how data is sorted.
#   - First we'll sort by `:Ticker`
#   - Then sort by `:Date`
# ------------------------------------------------------------------------------------------

stocks = loadtable("stocksample"; filenamecol = :Ticker, indexcols = [:Ticker, :Date])

# ------------------------------------------------------------------------------------------
# - Notice the printing style has changed.  A summary is printed when the display width is
# too narrow to print all the columns.  In Jupyter, the width is hardcoded as 80 characters,
# so the table actually fits in this case.  We can override this behavior with:
# ------------------------------------------------------------------------------------------

IndexedTables.set_show_compact!(false)

stocks

# ------------------------------------------------------------------------------------------
# ## Motivation for ND Sparse
#
# ### What was Apple's closing price on 1986-02-10?
#
# - For a `Table`, this requires a query
# - With `NDSparse`, this is just `getindex`
#     - Load data as `NDSparse` and use arbitrary indexes (`String` and `Date`) to look up
# the closing price:
# ------------------------------------------------------------------------------------------

# Load data as NDSparse:

stocksnd = loadndsparse("stocksample", filenamecol=:Ticker, indexcols = [:Ticker, :Date])

# Get the value associated with Apple and 1986-02-10:

stocksnd["aapl.us.txt", Date(1986, 2, 10)].Close

# ------------------------------------------------------------------------------------------
# # Saving Tables to Disk
#
# - **JuliaDB** can save tables into a binary format that can be loaded efficiently in
# future Julia sessions.
#     - `save(table, destination)`
# ------------------------------------------------------------------------------------------

save(stocks, "stocks.jdb")

# ------------------------------------------------------------------------------------------
# # Reloading a Saved Table
#
# - To load a saved table, we use the `load` function rather than `loadtable`.
# - This is typically **much** faster than reloading from the CSVs.
# ------------------------------------------------------------------------------------------

@time reloaded_stocks = load("stocks")

stocks == reloaded_stocks

# ------------------------------------------------------------------------------------------
# ---
#
# # Using  Selectors
#
# ### Selectors are powerful ways to select and manipulate data
#
# 1. `Integer`: column at position
# 2. `Symbol`: column by name
# 3. `Array`: itself
# 4. `Pair{Selection => Function}`: function mapped to selection
# 5. `Tuple` of selections: table of each selection
#
# ### Selectors show up in many places (everything in green)
#
# select(t, which)
# map(f, t; select)
# reduce(f, t; select)
# filter(f, t; select)
# groupby(f, t, by; select)
# groupreduce(f, t, by; select)
# join(f, l, r; how, lkey, rkey, lselect, rselect)
# groupjoin(f, l, r; how, lkey, rkey, lselect, rselect)
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Selector Examples
#
# - A single selection returns a Vector:
# ------------------------------------------------------------------------------------------

select(stocks, 1)[1:5]

select(stocks, :Date)[1:5]

# ------------------------------------------------------------------------------------------
# - Multiple selections return a table:
# ------------------------------------------------------------------------------------------

select(stocks, (1, :Date))[1:5]

ext: jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials