text:
#!/usr/bin/env julia
################################################################################
# This script is a short program for everyone who is interested in molecular
# dynamics. It will simulate a small amount of gas (default is argon) in a fixed
# volume. The potential will be desribed by the Lennard-Jones potential. The
# simulation itself is implemented with the velocity Verlet algorithm.
# Feel free to change/distribute/whatever the script.
# Also, if you have questions feel free to contact me at:
#     <cerebrosus at gmx.net>
################################################################################

push!(LOAD_PATH, ".")

using PyCall

@pyimport matplotlib.pyplot as pyplot
@pyimport matplotlib.animation as animation

using Distributions
using Potentials

################################################################################
# Global variables
# Removed most global variables cause of performance issues as explained in:
# http://docs.julialang.org/en/release-0.4/manual/performance-tips/

################################################################################
# Types

type Model
  # Different properties for the simulated system.
  steps::Int64
  particles::Int64
  sideLength::Float64
  halfSideLength::Float64
  initTemp::Float64

  # Physical properties of the particles.
  diameter::Float64
  mass::Float64

  # Properties for the Lennard-Jones-Potential.
  sigma::Float64
  epsilon::Float64
  ljm::Potentials.LennardJonesModel

  timeStep::Float64
  timeStep2::Float64

  rcut::Float64

  positions::Array{Float64, 3}
  velocities::Array{Float64, 2}

  accelerations::Array{Float64, 2}
  forces::Array{Float64, 2}
  temperatures::Array{Float64, 1}
end

################################################################################
# Functions

function correctVelocities!(velocities::Array{Float64, 2}, mass::Float64,
    particles::Int64, initTemp::Int64)

  s::Float64 = 0.0
  for i::Int64 = 1:size(velocities, 1), j::Int64 = 1:size(velocities, 2)
    s += velocities[i, j]^2
  end

  lambda = ((3 * (particles - 1) * 1.38064852e-23 * initTemp) / (mass * s))^(1/2)

  for i::Int64 = 1:size(velocities, 1), j::Int64 = 1:size(velocities, 2)
    velocities[i, j] *= lambda
  end
end

function initPositions!(positions::Array{Float64, 3}, sideLength::Float64,
      halfSideLength::Float64, particles::Int64)

  # How many particles on a line.
  numPartSide::Int64 = Int64(round(particles^(1.0 / 3.0)))
  numPartSide2::Int64 = numPartSide^2

  #The distance between every particle and the box limit.
  distance::Float64 = sideLength / (numPartSide + 1.0)

  parPos::Int64 = 0
  for x = 1:numPartSide, y = 1:numPartSide, z = 1:numPartSide
    parPos = z + (y - 1) * numPartSide + (x - 1) * numPartSide^2
    positions[1, parPos, 1] = -halfSideLength + distance * z
    positions[2, parPos, 1] = -halfSideLength + distance * y
    positions[3, parPos, 1] = -halfSideLength + distance * x
  end
end

"""
Constructor function for type Model.

return: The contructed model object.
"""
function Model()
  # To high takes to long. :-D
  steps::Int64 = 1000

  # Please choose a number with an integer result at ^(1/3).
  particles::Int64 = 64

  # The side length is one complete side length of the box.
  sideLength::Float64 = 6.0e-10
  halfSideLength::Float64 = sideLength / 2.0

  # The starting temperature of the particles which correlates to the velocity
  # of the particles.
  initTemp::Int64 = 50

  # Properties related to the atoms/molecules you choose.
  diameter::Float64 = 2.6e-10
  mass::Float64 = 6.69e-26

  # Different properties for the Lennard-Jones potential. You should change them
  # only if you know hat you are doing as they are well choosen to the default
  # system.
  sigma::Float64 = 3.4e-10
  epsilon::Float64 = 1.65e-21
  ljm::Potentials.LennardJonesModel = Potentials.LennardJonesModel(sigma,
      epsilon)

  # Properties of the simulated system.
  timeStep::Float64 = 0.001 * sigma * sqrt(mass / epsilon)
  timeStep2::Float64 = timeStep^2
  rcut::Float64 = 2 * sigma

  # Init the initial positions.
  positions::Array{Float64, 3} = fill(0.0, 3, particles, steps)
  initPositions!(positions, sideLength, halfSideLength, particles)

  # Create the velocity array and init the coordinates to the normal
  # distribution.
  dist = Normal(0.0, initTemp^(0.5))
  velocities::Array{Float64, 2} = rand(dist, 3, particles)
  correctVelocities!(velocities, mass, particles, initTemp);

  # Create the accelerations, temperatures and forces array, which are zero at
  # the beginning.
  accelerations::Array{Float64, 2} = fill(0.0, 3, particles)
  forces::Array{Float64, 2} = fill(0.0, 3, particles)
  temperatures::Array{Float64, 1} = fill(0.0, steps)

  return Model(steps, particles, sideLength, halfSideLength, initTemp,
               diameter, mass, sigma, epsilon, ljm, timeStep, timeStep2, rcut,
               positions, velocities, accelerations, forces, temperatures)
end

"""
Calculating the force between two particles.

posIndexA: Index of the first particle.
posIndexB: Index of the second particle.
step:      Step of simulation.
"""
function potential!(m::Model, posIndexA::Int64, posIndexB::Int64, step::Int64)
  diff::Array{Float64, 1} = fill(0.0, 3)
  diff[1] = m.positions[1, posIndexA, step] - m.positions[1, posIndexB, step]
  diff[2] = m.positions[2, posIndexA, step] - m.positions[2, posIndexB, step]
  diff[3] = m.positions[3, posIndexA, step] - m.positions[3, posIndexB, step]

  r::Float64 = norm(diff)
  if r < m.rcut
    pot::Float64 = ljp(m, r)

    m.forces[1, posIndexA] += pot * diff[1]
    m.forces[2, posIndexA] += pot * diff[2]
    m.forces[3, posIndexA] += pot * diff[3]
    m.forces[1, posIndexB] -= pot * diff[1]
    m.forces[2, posIndexB] -= pot * diff[2]
    m.forces[3, posIndexB] -= pot * diff[3]
  end

  return 0.0
end

"""
Adjust the given position to fit in the given area. There may be a closed or
a periodic box.

position:     The positions to correct.
velocity:     The current velocity of the particle.
sideLength:   The length of a box with equal side length.
systemClosed: True if the system is a closed box and false if it is periodic.
"""
function adjustPosition!(position::Array{Float64, 1},
    velocities::Array{Float64, 2}, index::Int64, halfSideLength::Float64,
    systemClosed::Bool = true)

  # We suggest, that the center of the box has the coordinates [0, 0, 0].
  if systemClosed::Bool == false
    for i::Int64 = 1 : 3
      if position[i] > halfSideLength
        position[i] = -halfSideLength
      elseif position[i] < -halfSideLength
        position[i] = halfSideLengh
      end
    end
  else
    for i::Int64 = 1 : 3
      if position[i] > halfSideLength || position[i] < -halfSideLength
        velocities[i, index] *= -1
      end
    end
  end
end

"""
Calculate the temperature of the system based on the mass and velocity of
every particle.

temperatures: The temperatures object.
index:        The index for temperature step.
mass:         Mass of the particle.
velocities:   Velocities object of all particle as a vector with n-coordinates.
"""
function temperature!(temperatures::Array{Float64, 1}, index::Int64,
    mass::Float64, velocities::Array{Float64, 2}, particles::Int64)

  const kB::Float64 = 1.38064852e-23
  s::Float64 = 0.0
  for i::Int64 = 1:size(velocities, 1), j::Int64 = 1:size(velocities, 2)
    s += velocities[i, j]^2
  end

  temperatures[index] =  s * mass / (kB * 3.0 * (particles - 1))
end

"Main simulation function for running the system."
function simulation(m::Model)
  position::Array{Float64, 1} = fill(0.0, 3)

  # Update forces
  for j = 1:m.particles
    for k = j + 1:m.particles
      potential!(m, j, k, 1)
    end

    # Update accelerations
    m.accelerations[1, j] = m.forces[1, j] / m.mass
    m.accelerations[2, j] = m.forces[2, j] / m.mass
    m.accelerations[3, j] = m.forces[3, j] / m.mass
  end

  # Running main loop
  for i = 1 : m.steps - 1
    for j = 1 : m.particles
      # Update velocities + 0.5 timestep.
      m.velocities[1, j] = m.velocities[1, j] + 0.5 * m.accelerations[1, j] * m.timeStep
      m.velocities[2, j] = m.velocities[2, j] + 0.5 * m.accelerations[2, j] * m.timeStep
      m.velocities[3, j] = m.velocities[3, j] + 0.5 * m.accelerations[3, j] * m.timeStep

      # Update particle position followed by a correction of particle position
      # and velocitiy.
      position[1] = m.positions[1, j, i] + m.velocities[1, j] * m.timeStep
      position[2] = m.positions[2, j, i] + m.velocities[2, j] * m.timeStep
      position[3] = m.positions[3, j, i] + m.velocities[3, j] * m.timeStep
      adjustPosition!(position, m.velocities, j, m.halfSideLength)
      m.positions[1, j, i + 1] = position[1]
      m.positions[2, j, i + 1] = position[2]
      m.positions[3, j, i + 1] = position[3]
    end
    for j = 1 : m.particles
      # Update forces
      for k = j + 1 : m.particles
        potential!(m, j, k, i + 1)
      end

      # Update accelerations
      m.accelerations[1, j] = m.forces[1, j] / m.mass
      m.accelerations[2, j] = m.forces[2, j] / m.mass
      m.accelerations[3, j] = m.forces[3, j] / m.mass

      # Update velocities + full timestep.
      m.velocities[1, j] = m.velocities[1, j] + 0.5 * m.accelerations[1, j] * m.timeStep
      m.velocities[2, j] = m.velocities[2, j] + 0.5 * m.accelerations[2, j] * m.timeStep
      m.velocities[3, j] = m.velocities[3, j] + 0.5 * m.accelerations[3, j] * m.timeStep
    end

    # Calculate and save temperature for every step.
    temperature!(m.temperatures, i, m.mass, m.velocities, m.particles)
    fill!(m.forces, 0.0)
  end
end

function showAnimationPlot(m::Model)
  plots = [pyplot.plot([m.positions[1, i, j] for i in 1:size(m.positions, 2)],
                [m.positions[2, i, j] for i in 1:size(m.positions, 2)], "ro")
           for j in 1:size(m.positions, 3)]

  pyplot.axis([-6e-10, 6e-10, -6e-10, 6e-10])
  fig1 = pyplot.figure(1)
  animation.ArtistAnimation(fig1, plots, interval = 1, blit = false)
  pyplot.show()
end

function showTemperaturePlot(m::Model)
  pyplot.title("Molecular Dynamics Simulation", fontsize = 14)
  pyplot.plot([i for i in 1:m.steps], m.temperatures)
  pyplot.axis([0, m.steps, m.initTemp - 20, m.initTemp + 20])
  pyplot.xlabel("Timesteps")
  pyplot.ylabel("Temperature [K]")
  pyplot.show()
end

function writeAnimation(m::Model)
  for i = 1:m.steps
    file = open(string("result/sim-", i, ".csv"), "w")

    for j = 1:m.particles
      write(file, join(m.positions[:, j, i], ","), "\n")
    end

    close(file)
  end
end

m = Model()
@time simulation(m)

showTemperaturePlot(m)

max_stars_repo_name: Christian-Kr/simLJP-julia

max_issues_repo_name: Christian-Kr/simLJP-julia

max_forks_repo_name: Christian-Kr/simLJP-julia


text:
using Distributions
using StatsBase, Clustering
using MLBase # Just for cross validation

################################################################################
#####
##### Helper functions
#####
"""
Gets number of degrees of freedom, which is equal to the
number of nonzero terms + 1 for noise
"""
my_dof(m::DynamicalSystemModel) = length(get_nonzero_terms(m)) + 1

#####
##### Information theory
#####
"""
Akaike Information Criterion (AIC) for SINDy models
"""
my_aic(m::DynamicalSystemModel, dat, predictor; dist=Normal()) =
    -2loglikelihood(dist, (m(dat) .- predictor)') + 2my_dof(m)

"""
Corrected Akaike Information Criterion.
    Used for small sample sizes (Hurvich and Tsai 1989)
"""
function my_aicc(m::DynamicalSystemModel, dat, predictor; dist=Normal())
    k = my_dof(m)
    n = length(dat) # equal to numel(dat)
    correction = 2k*(k+1)/(n-k-1)
    return my_aic(m, dat, predictor, dist=dist) + correction
end

"""
Akaike Information Criterion (AIC), without being 'clever' and using an error
    distribution as in the above function
"""
simple_aic(m::DynamicalSystemModel, dat, predictor; dist=Normal()) =
    2log(rss_sindy_derivs(m, dat, predictor)) + 2my_dof(m)

################################################################################
#####
##### Cross validation
#####

"""
Computes the residual sum of squares error (derivatives) for a SINDy model
"""
rss_sindy_derivs(m::DynamicalSystemModel, dat, predictor) =
    sum((m(dat, 0) .- predictor).^2);

"""
Computes the residual sum of squares error (simulated) for a SINDy model
"""
rss_sindy_integrate(m::DynamicalSystemModel, dat, predictor) =
    sum((simulate_model(m, dat[:,1]) .- predictor).^2);

"""
Computes the k-folds cross validation of a SINDYc model
    From: https://mlbasejl.readthedocs.io/en/latest/crossval.html
"""
function sindyc_cross_validate(m::sindycModel,
                            dat, predictor;
                            dist=Normal())
    # BUG??
    n = size(dat, 2);
    @show n
    @show size(m.U)
    my_retrain(train_ind) =
        sindyc_retrain(m,
        dat[:,train_ind], predictor[:,train_ind], m.U[:,train_ind])

    scores = cross_validate(
        (ind)->sindyc_retrain(m, dat[:,ind], predictor[:,ind], m.U[:,ind]),
        # TODO: the rss shouldn't need the data itself
        (new_model, test_inds) -> rss_sindy_derivs(new_model,
                        dat[:, test_inds], predictor[:, test_inds]),  # evaluation function
        n,              # total number of samples
        Kfold(n, 3))    # cross validation plan: 3-fold
    return mean(scores)
end

"""
Computes the k-folds cross validation of a SINDY model
    From: https://mlbasejl.readthedocs.io/en/latest/crossval.html
"""
function sindy_cross_validate(m::sindyModel,
                            dat, predictor;
                            dist=Normal())
    n = size(dat, 2);
    my_retrain(train_ind) =
        sindy_retrain(m,
        dat[:,train_ind], predictor[:,train_ind])

    scores = cross_validate(
        my_retrain,
        # TODO: the rss shouldn't need the data itself
        (new_model, test_inds) -> rss_sindy_derivs(new_model,
                        dat[:, test_inds], predictor[:, test_inds]),  # evaluation function
        n,              # total number of samples
        Kfold(n, 3))    # cross validation plan: 3-fold
    return mean(scores)
end

"""
Calculates the error in coefficients between two models
"""
function calc_coefficient_error(m1::DynamicalSystemModel, m2::DynamicalSystemModel)
    err = sum((m1.A .- m2.A).^2)
    return err
end

# sindy_cross_validate(m::sindycModel, dat, predictor; dist=Normal()) =
#     sindy_cross_validate(m, dat, predictor)

################################################################################
#####
##### Ensemble testing functions
#####
"""
sindyc_ensemble(X, X_grad, library, val_list;
                U=nothing,
                selection_criterion=my_aic,
                selection_dist=Normal(),
                sparsification_mode="quantile",
                use_clustering_minimization=false)

Loop over several values of sparsity and choose the best model via
    MINIMIZING the 'selection_criterion' function
        selection_criterion() should have signature:
        selection_criterion(model, data, predictor)

    Also returns all the sparsity values and all the models
"""
function sindyc_ensemble(X, X_grad, library, val_list;
                U=nothing, ts=nothing,
                selection_criterion=my_aic,
                selection_dist=Normal(),
                sparsification_mode="quantile",
                use_clustering_minimization=false,
                var_names=nothing)
    if ts == nothing
        ts = 1:size(X,2)
    end
    n = length(val_list)
    all_models = Vector(undef, n)
    all_criteria = zeros(n)
    # Convenience function
    # TODO: Update with optimizer interface
    make_model(x) = if sparsification_mode == "quantile"
        sindyc(X, X_grad, U, ts; library=library,
                                use_lasso=true,
                                quantile_threshold=x,
                                var_names=var_names)
    elseif sparsification_mode == "num_terms"
        sindyc(X, X_grad, U, ts; library=library,
                                use_lasso=true,
                                quantile_threshold=nothing,
                                num_terms=x,
                                var_names=var_names)
    else
        error("Unknown sparsification mode")
    end

    for (i, val) in enumerate(val_list)
        m = make_model(val) # Defined above with if statement
        all_criteria[i] = selection_criterion(m,
                X, X_grad, dist=selection_dist)
        all_models[i] = m
    end

    if !use_clustering_minimization
        best_criterion, best_index = findmin(all_criteria)
    else
        # Use AIC, but get the sparsest of the similar models
        #   Very helpful when the noise is highly Non-Gaussian
        tmp = reshape(all_criteria, length(all_criteria),1)
        c = kmeans(tmp', 6) # TODO: learn the number of clusters
        ind = findmin(all_criteria)[2]
        best_model_clust = c.assignments[ind]
        clust_ind_in_all_models = findall(c.assignments.==best_model_clust)
        nnz = sum.(val_list)[clust_ind_in_all_models]
        best_criterion, model_ind_in_clust = findmin(nnz)
        best_index = clust_ind_in_all_models[model_ind_in_clust]
    end

    return (best_model=all_models[best_index],
            best_criterion=best_criterion,
            all_criteria=all_criteria,
            all_models=all_models,
            best_index=best_index)
end

"""
An options struct to make it easier to call sindyc_ensemble
"""
function get_sindyc_ensemble_parameters()
    return Dict(:U=>nothing,
        :selection_criterion=>my_aicc,
        :selection_dist=>Normal(),
        :sparsification_mode=>"num_terms", # NOTE: different from raw object above
        :use_clustering_minimization=>false)
end

################################################################################
"""
sindy_ensemble(model_template, X, X_grad, val_list;
                selection_criterion=rss_sindy_derivs)

Loop over several values of optimizer hyperparameters via
    MINIMIZING the 'selection_criterion' function
        selection_criterion() should have signature:
        selection_criterion(model, data, predictor)

    Also returns all the sparsity values and all the models
"""
function sindy_ensemble(model_template::sindyModel,
                X, X_grad, val_list;
                selection_criterion=rss_sindy_derivs,
                test_ind=1:500)
    # TODO: check these test indices
    n = length(val_list)
    all_models = Vector(undef, n)
    all_criteria = zeros(n)

    opt = model_template.optimizer
    for (i, val) in enumerate(val_list)
        # println("Testing optimizer $(typeof(model_template.optimizer)) with value:")
        # @show val
        # Update optimzer and retrain
        model_template.optimizer = copy_optimizer(opt, val)
        this_m = sindy_retrain(model_template, X, X_grad)
        all_criteria[i] = selection_criterion(
            this_m, X[:,test_ind], X_grad[:,test_ind])
        all_models[i] = this_m
    end

    best_criterion, best_index = findmin(all_criteria)

    return (best_model=all_models[best_index],
            best_criterion=best_criterion,
            all_criteria=all_criteria,
            all_models=all_models,
            best_index=best_index)
end

#
export sindyc_ensemble, sindy_ensemble, get_sindyc_ensemble_parameters,
        my_aic, my_aicc, my_dof, rss_sindy_derivs,
        sindy_cross_validate, sindyc_cross_validate,
        calc_coefficient_error

ext: jl


text:
using DataFrames;

"""
sw_align(seq1::String, seq2::String, smat::DataFrame, gapp::Float64, gapcp::Float64, out_fp::String)
Align FASTA sequences SEQ1 and SEQ2 using the Smith-Waterman algorithm with gap penalty
GP and extension penalty EP and according to the scoring matrix SMAT.
Stores the result in OUT_FP.

"""
function sw_align(seq1, seq2, smat, gp, ep, out_fp; write_align = true)

    M = length(seq1);
    N = length(seq2);

    H = zeros(M+1, N+1);
    Wk = zeros(M+1);
    Wm = zeros(N+1);

    for i in 2:(M+1)
        for j in 2:(N+1)

            diag = (H[i-1, j-1] + smat[smat[:, :AA] .== string(uppercase(seq1[i-1])),
                                            Symbol(uppercase(seq2[j-1]))])[1];

            # W_k = extension penalty * length of extension + gap penalty
            # We don't take into account extension length until after a gap is allowed
            vert = maximum([(H[i-k, j] - (k-1)*ep - gp) for k in 1:(i-1)]);
            horiz = maximum([(H[i, j-k] - (k-1)*ep - gp) for k in 1:(j-1)]);
            opt = max(diag, vert, horiz, 0);
            H[i, j] = opt;
        end
    end

    #println(DataFrame(H));
    if write_align
        sw_traceback(seq1, seq2, H, out_fp)
    else
        return findmax(H)[1], H
    end

end

"""
many_sw_align(seqs::Array, smat::DataFrame, gp::Float64, ep::Float64, out_fp::String)
Align a set of sequences in SEQS with the Smith-Waterman algorithm. Use the gap penalty
GP and extension penalty EP and score according to the scoring matrix SMAT.
Stores the result in OUT_FP.

"""
function many_sw_align(seqs, smat, gp, ep, out_fp; write_output=true, normalize=false,
                        write_alignments = true)

    println(gp)
    println(ep)

    scores = [];
    Hs = [];

    progress_i = 1

    for seq_pair in seqs

        println(string("Aligning pair ", progress_i, " of ", length(seqs)))

        s, H = sw_align(seq_pair[1], seq_pair[2], smat, gp, ep, out_fp; write_align=false)

        push!(Hs, H)

        if normalize
            s = s / min(length(seq_pair[1]), length(seq_pair[2]))
        end

        push!(scores, s)

        progress_i += 1

    end

    if write_output && !write_alignments
        open(out_fp, "w") do f
            for s in scores
                write(f, string(s, '\t'));
            end
            write(f, '\n');
        end
    elseif write_alignments && write_output
        sw_traceback_many(seqs, Hs, out_fp)
    else
        scores
    end

end

"""
sw_traceback(seq1::String, seq2::String, H::DataFrame, out_fp::String)
Reconstruct the optimal alignment given the alignment matrix H and the two
alignments SEQ1 and SEQ2. Write output to OUT_FP.
"""

function sw_traceback(seq1, seq2, H, out_fp)

    m = length(seq1)+1;
    n = length(seq2)+1;

    align1 = "";
    matches = "";
    align2 = "";

    mii = ind2sub(H, findmax(H)[2]);

    i, j = mii[1], mii[2];

    _end = false
    while !_end
        diag = H[i-1, j-1];
        up = H[i-1, j];
        left = H[i, j-1];

        if diag >= up && diag >= left
            if diag == 0
                _end = true
            else
                align1 = string(seq1[i-1], align1)
                align2 = string(seq2[j-1], align2)
                i -= 1
                j -= 1
            end
        elseif up > diag && up >= left
            if up == 0
                _end = true
            else
                align1 = string(seq1[i-1], align1)
                align2 = string("-", align2)
                i -= 1
            end
        else
            if left == 0
                _end = true
            else
                align1 = string("-", align1)
                align2 = string(seq2[j-1], align2)
                j -= 1
            end
        end
    end

    align1 = string(seq1[i-1], align1)
    align2 = string(seq2[j-1], align2)

    open(out_fp, "w") do f
        write(f, string(align1, '\n'))
        write(f, string(align2, '\n'))
        write(f, string("Score: ", findmax(H)[1], '\n'))
    end

    return

end

"""
sw_traceback_many(seqs::Array, Hs::Array, out_fp::String)
Reconstruct the optimal alignments for all the sequences pairs in SEQS,
where there is an accompanying alignment matrix H in Hs. Write output to OUT_FP.
"""

function sw_traceback_many(seqs, Hs, out_fp)

    open(out_fp, "w") do f
        for i in 1:length(seqs)

            seq1 = seqs[i][1];
            seq2 = seqs[i][2];
            H = Hs[i];

            m = length(seq1)+1;
            n = length(seq2)+1;

            align1 = "";
            matches = "";
            align2 = "";

            mii = ind2sub(H, findmax(H)[2]);

            i, j = mii[1], mii[2];

            _end = false
            while !_end
                diag = H[i-1, j-1];
                up = H[i-1, j];
                left = H[i, j-1];

                if diag >= up && diag >= left
                    if diag == 0
                        _end = true
                    else
                        align1 = string(seq1[i-1], align1)
                        align2 = string(seq2[j-1], align2)
                        i -= 1
                        j -= 1
                    end
                elseif up > diag && up >= left
                    if up == 0
                        _end = true
                    else
                        align1 = string(seq1[i-1], align1)
                        align2 = string("-", align2)
                        i -= 1
                    end
                else
                    if left == 0
                        _end = true
                    else
                        align1 = string("-", align1)
                        align2 = string(seq2[j-1], align2)
                        j -= 1
                    end
                end
            end

            align1 = string(seq1[i-1], align1)
            align2 = string(seq2[j-1], align2)

            write(f, string(align1, '\n'))
            write(f, string(align2, '\n'))
            write(f, string("Score: ", findmax(H)[1], '\n'))

        end
    end

end


text:
# Implement code blocks 3.11 through 3.19 in Statistical Rethinking, 2nd edition.

using Distributions
using Statistics
using StatsBase
using LinearAlgebra

grid_size = 1000
num_observations = 3
num_samples = 1000

# Code block 3.11
param_grid = range(0, 1, length=grid_size)
prior = ones((grid_size))
likelihood = [pdf(Binomial(3, p), num_observations) for p in param_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)

@assert abs(sum(posterior) - 1) < 1e-4 "Posterior does not sum to 1."
@assert minimum(posterior) >= 0
@assert maximum(posterior) <= 1

samples = sample(param_grid, Weights(posterior), num_samples; replace=true, ordered=false)

# Code block 3.12
function central_interval(samples, proba::AbstractFloat)
    @assert proba > 0
    @assert proba < 1
    quantile(samples, [0.5 - proba / 2, 0.5 + proba / 2])
end

central_interval(samples, 0.5)

# Code block 3.13

# UPDATE: all the prob. prog. packages have it right. Because they *sort* the samples
# first, they're really working with the ECDF. Subtracting in the domain of the ECDF
# such that the range always has at least `proba` mass.

# - I'm 85% sure the StatisticalRethinking.jl package has this wrong, it assumes
#   some mass on both sides of the tail, which is not correct.

function high_density_interval(samples, proba::AbstractFloat)
    @assert proba >= 0
    @assert proba <= 1

    count_threshold = proba * length(samples)
    hist = fit(Histogram, samples; closed=:right, nbins=100)

    interval_mass = 0
    start_index = 1
    min_interval_length = Inf
    candidate = (NaN, NaN)
    candidate_mass = 0

    for (i, w) in enumerate(hist.weights)
        interval_mass += w

        while interval_mass >= count_threshold
            interval_length = i - start_index

            if (interval_length < min_interval_length) ||
                    ((interval_length == min_interval_length) &&
                     interval_mass > candidate_mass)

                candidate = (start_index, i)
                candidate_mass = interval_mass
                min_interval_length = interval_length
            end

            interval_mass = interval_mass - hist.weights[start_index]
            start_index += 1
        end
    end

    interval = [hist.edges[1][candidate[1]], hist.edges[1][candidate[2]]]
end

high_density_interval(samples, 0.5)

# Code block 3.14 - MAP of the explicit posterior
_, i = findmax(posterior)
posterior_mode = param_grid[i]

# Code block 3.15 - MAP based on posterior samples
# - Based on a density estimate again

# Code block 3.16
mean(samples)
median(samples)

# Code block 3.17 - mean absolute error loss for a point estimate
function mae_loss(point_estimate, param_grid, posterior)
    sum(posterior .* abs.(point_estimate .- param_grid))
end

loss = mae_loss(0.5, param_grid, posterior)

# Code blocks 3.18 and 3.19 - mean absolute error for all param values
# - The book uses `sapply` in R. I don't know how to vectorize like that yet in Julia,
#   so let's just loop.
losses = [mae_loss(p, param_grid, posterior) for p in param_grid]

_, i = findmin(losses)

best_point_estimate = param_grid[i]


text:
"""
 Additive identity.
"""
struct Zero <: ChainRulesCore.AbstractZero end
@inline Base.:+(::Zero, a) = a
@inline Base.:+(a, ::Zero) = a
@inline Base.:+(::Zero, ::Zero) = Zero()
# @inline Base.:+(sptr::StackPointer, ::Zero, a) = (sptr, a)
# @inline Base.:+(sptr::StackPointer, a, ::Zero) = (sptr, a)
# @inline Base.:+(sptr::StackPointer, ::Zero, ::Zero) = (sptr, Zero())
@inline LoopVectorization.vadd(::Zero, a) = a
@inline LoopVectorization.vadd(a, ::Zero) = a
@inline LoopVectorization.vadd(::Zero, ::Zero) = Zero()
# @inline SIMDPirates.vadd(sptr::StackPointer, ::Zero, a) = (sptr, a)
# @inline SIMDPirates.vadd(sptr::StackPointer, a, ::Zero) = (sptr, a)
# @inline SIMDPirates.vadd(sptr::StackPointer, ::Zero, ::Zero) = (sptr, Zero())

@inline Base.:*(::Zero, a) = Zero()
@inline Base.:*(a, ::Zero) = Zero()
@inline Base.:*(::Zero, ::Zero) = Zero()
# @inline Base.:*(sptr::StackPointer, ::Zero, a) = (sptr, Zero())
# @inline Base.:*(sptr::StackPointer, a, ::Zero) = (sptr, Zero())
# @inline Base.:*(sptr::StackPointer, ::Zero, ::Zero) = (sptr, Zero())
@inline LoopVectorization.vmul(::Zero, a) = Zero()
@inline LoopVectorization.vmul(a, ::Zero) = Zero()
@inline LoopVectorization.vmul(::Zero, ::Zero) = Zero()
# @inline SIMDPirates.vmul(sptr::StackPointer, ::Zero, a) = (sptr, Zero())
# @inline SIMDPirates.vmul(sptr::StackPointer, a, ::Zero) = (sptr, Zero())
# @inline SIMDPirates.vmul(sptr::StackPointer, ::Zero, ::Zero) = (sptr, Zero())

"""
Multiplicative identity.
"""
struct One <: ChainRulesCore.AbstractDifferential end
@inline Base.:*(::One, a) = a
@inline Base.:*(a, ::One) = a
@inline Base.:*(::One, ::One) = One()
# @inline Base.:*(sptr::StackPointer, ::One, a) = (sptr, a)
# @inline Base.:*(sptr::StackPointer, a, ::One) = (sptr, a)
# @inline Base.:*(sptr::StackPointer, ::One, ::One) = (sptr, One())
@inline LoopVectorization.vmul(::One, a) = a
@inline LoopVectorization.vmul(a, ::One) = a
@inline LoopVectorization.vmul(::One, ::One) = One()
# @inline LoopVectorization.vmul(sptr::StackPointer, ::One, a) = (sptr, a)
# @inline LoopVectorization.vmul(sptr::StackPointer, a, ::One) = (sptr, a)
# @inline LoopVectorization.vmul(sptr::StackPointer, ::One, ::One) = (sptr, One())
@inline Base.inv(::One) = One()
@inline Base.FastMath.inv_fast(::One) = One()
@inline Base.:/(a, ::One) = a
@inline Base.FastMath.div_fast(a, ::One) = a

Base.size(::Zero) = Core.tuple()
Base.size(::One) = Core.tuple()


text: hello() = "Hello, World!"


text:
"""
1. Grey results are PRESENT_NOWHERE.
2. Yellow results are PRESENT_ELSEWHERE.
3. Green results are PRESENT_HERE.
"""
@enum ResultCode begin
    PRESENT_NOWHERE
    PRESENT_ELSEWHERE
    PRESENT_HERE
end

"""
Convert Wordle emojis to ResultCode objects.
"""
function ResultCode(c::Char)
    if c === '⬜'
        PRESENT_NOWHERE
    elseif c === '🟨'
        PRESENT_ELSEWHERE
    elseif c === '🟩'
        PRESENT_HERE
    else
        error("Invalid code character")
    end
end


text:
using DiffEqFlux, GeometricFlux, GraphSignals, OrdinaryDiffEq, Test

# Fully Connected Graph
adj_mat = FeaturedGraph(Float32[0 1 1 1
                                1 0 1 1
                                1 1 0 1
                                1 1 1 0])

features = [-10.0f0 -9.0f0 9.0f0 10.0f0
              0.0f0  0.0f0 0.0f0  0.0f0]

target = [1.0 1.0 0.0 0.0
          0.0 0.0 1.0 1.0]

model = Flux.Chain(
    NeuralODE(
        GCNConv(adj_mat, 2=>2),
        (0.f0, 1.f0), Tsit5(), save_everystep = false,
        reltol = 1e-3, abstol = 1e-3, save_start = false
    ),
    x -> reshape(cpu(x), size(x)[1:2])
)

ps = Flux.params(model)
opt = ADAM(0.1)

initial_loss = Flux.Losses.logitcrossentropy(model(features), target)

# for i in 1:100
#     gs = gradient(() -> Flux.Losses.logitcrossentropy(model(features), target), ps)
#     Flux.Optimise.update!(opt, ps, gs)
# end
updated_loss = Flux.Losses.logitcrossentropy(model(features), target)

@test_broken updated_loss < initial_loss


text:
push!(LOAD_PATH,"../src/Retriever.jl")
include("../src/Retriever.jl")

Pkg.add("Documenter")
using Documenter, DocumenterTools

makedocs(
    modules = [Retriever],
    clean = false,
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        # Or set prettyurls = false
        prettyurls = !("local" in ARGS)),
    build = "build",
    sitename = "Retriever.jl",
    authors = "Ethan White",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "intro.md",
        "Installation Guide" => "tutorial.md",
        "Developer's Guide" => "developer.md",
        "Source" => "index.md",
        "Command Documentation" => "lib/public.md",
        "Code of Conduct" => "CODE_OF_CONDUCT.md"
    ]
)

deploydocs(
    repo = "github.com/weecology/Retriever.jl.git",
    target = "build",
    deps = nothing,
    julia = "",
    make = nothing
)