text:
# # Horizontal convection example
#
# In "horizontal convection", a non-uniform buoyancy is imposed on top of an initially resting fluid.
#
# This example demonstrates:
#
#   * How to use `ComputedField`s for output.
#   * How to post-process saved output using `FieldTimeSeries`.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## Horizontal convection
#
# We consider here two-dimensional horizontal convection of an incompressible flow ``\boldsymbol{u} = (u, w)``
# on the ``(x, z)``-plane (``-L_x/2 \le x \le L_x/2`` and ``-H \le z \le 0``). The flow evolves
# under the effect of gravity. The only forcing on the fluid comes from a prescribed, non-uniform
# buoyancy at the top-surface of the domain.
#
# We start by importing `Oceananigans` and `Printf`.

using Oceananigans
using Printf

# ### The grid

H = 1.0          # vertical domain extent
Lx = 2H          # horizontal domain extent
Nx, Nz = 128, 64 # horizontal, vertical resolution

grid = RegularRectilinearGrid(size = (Nx, Nz),
                                 x = (-Lx/2, Lx/2),
                                 z = (-H, 0),
                              halo = (3, 3),
                          topology = (Bounded, Flat, Bounded))

# ### Boundary conditions
#
# At the surface, the imposed buoyancy is
# ```math
# b(x, z = 0, t) = - b_* \cos (2 \pi x / L_x) \, ,
# ```
# while zero-flux boundary conditions are imposed on all other boundaries. We use free-slip
# boundary conditions on ``u`` and ``w`` everywhere.

b★ = 1.0

@inline bˢ(x, y, t, p) = - p.b★ * cos(2π * x / p.Lx)

b_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(bˢ, parameters=(; b★, Lx)))

# ### Non-dimensional control parameters and Turbulence closure
#
# The problem is characterized by three non-dimensional parameters. The first is the domain's
# aspect ratio, ``L_x / H`` and the other two are the Rayleigh (``Ra``) and Prandtl (``Pr``)
# numbers:
#
# ```math
# Ra = \frac{b_* L_x^3}{\nu \kappa} \, , \quad \text{and}\, \quad Pr = \frac{\nu}{\kappa} \, .
# ```
#
# The Prandtl number expresses the ratio of momentum over heat diffusion; the Rayleigh number
# is a measure of the relative importance of gravity over viscosity in the momentum equation.
#
# For a domain with a given extent, the nondimensional values of ``Ra`` and ``Pr`` uniquely
# determine the viscosity and diffusivity, i.e.,
#
# ```math
# \nu = \sqrt{\frac{Pr b_* L_x^3}{Ra}} \quad \text{and} \quad \kappa = \sqrt{\frac{b_* L_x^3}{Pr Ra}} \, .
# ```
#
# We use isotropic viscosity and diffusivities, `ν` and `κ` whose values are obtain from the
# prescribed ``Ra`` and ``Pr`` numbers. Here, we use ``Pr = 1`` and ``Ra = 10^8``:

Pr = 1.0    # Prandtl number
Ra = 1e8    # Rayleigh number

ν = sqrt(Pr * b★ * Lx^3 / Ra)  # Laplacian viscosity
κ = ν * Pr                     # Laplacian diffusivity
nothing # hide

# ## Model instantiation
#
# We instantiate the model with the fifth-order WENO advection scheme, a 3rd order
# Runge-Kutta time-stepping scheme, and a `BuoyancyTracer`.

model = NonhydrostaticModel(
           architecture = CPU(),
                   grid = grid,
              advection = WENO5(),
            timestepper = :RungeKutta3,
                tracers = :b,
               buoyancy = BuoyancyTracer(),
                closure = IsotropicDiffusivity(ν=ν, κ=κ),
    boundary_conditions = (b=b_bcs,)
    )

# ## Simulation set-up
#
# We set up a simulation that runs up to ``t = 40`` with a `JLD2OutputWriter` that saves the flow
# speed, ``\sqrt{u^2 + w^2}``, the buoyancy, ``b``, andthe vorticity, ``\partial_z u - \partial_x w``.
#
# ### The `TimeStepWizard`
#
# The TimeStepWizard manages the time-step adaptively, keeping the Courant-Freidrichs-Lewy
# (CFL) number close to `0.75` while ensuring the time-step does not increase beyond the
# maximum allowable value for numerical stability.

max_Δt = 1e-1
wizard = TimeStepWizard(cfl=0.75, Δt=1e-2, max_change=1.2, max_Δt=max_Δt)

# ### A progress messenger
#
# We write a function that prints out a helpful progress message while the simulation runs.

CFL = AdvectiveCFL(wizard)

start_time = time_ns()

progress(sim) = @printf("i: % 6d, sim time: % 1.3f, wall time: % 10s, Δt: % 1.4f, CFL: %.2e\n",
                        sim.model.clock.iteration,
                        sim.model.clock.time,
                        prettytime(1e-9 * (time_ns() - start_time)),
                        sim.Δt.Δt,
                        CFL(sim.model))
nothing # hide

# ### Build the simulation
#
# We're ready to build and run the simulation. We ask for a progress message and time-step update
# every 50 iterations,

simulation = Simulation(model, Δt = wizard, iteration_interval = 50,
                                                     stop_time = 40.0,
                                                      progress = progress)

# ### Output
#
# We use `ComputedField`s to diagnose and output the total flow speed, the vorticity, ``\zeta``,
# and the buoyancy, ``b``. Note that `ComputedField`s take "AbstractOperations" on `Field`s as
# input:

u, v, w = model.velocities # unpack velocity `Field`s
b = model.tracers.b        # unpack buoyancy `Field`

## total flow speed
s = ComputedField(sqrt(u^2 + w^2))

## y-component of vorticity
ζ = ComputedField(∂z(u) - ∂x(w))

outputs = (s = s, b = b, ζ = ζ)
nothing # hide

# We create a `JLD2OutputWriter` that saves the speed, and the vorticity.
# We then add the `JLD2OutputWriter` to the `simulation`.

saved_output_prefix = "horizontal_convection"
saved_output_filename = saved_output_prefix * ".jld2"

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                  field_slicer = nothing,
                                                      schedule = TimeInterval(0.5),
                                                        prefix = saved_output_prefix,
                                                         force = true)
nothing # hide

# Ready to press the big red button:

run!(simulation)

# ## Load saved output, process, visualize
#
# We animate the results by opening the JLD2 file, extracting data for the iterations we ended
# up saving at, and ploting the saved fields. From the saved buoyancy field we compute the
# buoyancy dissipation, ``\chi = \kappa |\boldsymbol{\nabla} b|^2``, and plot that also.
#
# To start we load the saved fields are `FieldTimeSeries` and prepare for animating the flow by
# creating coordinate arrays that each field lives on.

using JLD2, Plots
using Oceananigans
using Oceananigans.Fields
using Oceananigans.AbstractOperations: volume

saved_output_prefix = "horizontal_convection"
saved_output_filename = saved_output_prefix * ".jld2"

## Open the file with our data
file = jldopen(saved_output_filename)

κ = file["closure/κ/b"]

s_timeseries = FieldTimeSeries(saved_output_filename, "s")
b_timeseries = FieldTimeSeries(saved_output_filename, "b")
ζ_timeseries = FieldTimeSeries(saved_output_filename, "ζ")

times = b_timeseries.times

## Coordinate arrays
xs, ys, zs = nodes(s_timeseries[1])
xb, yb, zb = nodes(b_timeseries[1])
xζ, yζ, zζ = nodes(ζ_timeseries[1])
nothing # hide

# Now we're ready to animate.

@info "Making an animation from saved data..."

anim = @animate for i in 1:length(times)

    ## Load fields from `FieldTimeSeries` and compute χ
    t = times[i]

    s_snapshot = interior(s_timeseries[i])[:, 1, :]
    ζ_snapshot = interior(ζ_timeseries[i])[:, 1, :]

    b = b_timeseries[i]
    χ = ComputedField(κ * (∂x(b)^2 + ∂z(b)^2))
    compute!(χ)

    b_snapshot = interior(b)[:, 1, :]
    χ_snapshot = interior(χ)[:, 1, :]

    ## determine colorbar limits and contour levels
    slim = 0.6
    slevels = vcat(range(0, stop=slim, length=31), [slim])

    blim = 0.6
    blevels = vcat([-1], range(-blim, stop=blim, length=51), [1])

    ζlim = 9
    ζlevels = range(-ζlim, stop=ζlim, length=41)

    χlim = 0.025
    χlevels = range(0, stop=χlim, length=41)

    @info @sprintf("Drawing frame %d:", i)

    kwargs = (      xlims = (-Lx/2, Lx/2),
                    ylims = (-H, 0),
                   xlabel = "x / H",
                   ylabel = "z / H",
              aspectratio = 1,
                linewidth = 0)

    s_plot = contourf(xs, zs, s_snapshot';
                      clims = (0, slim), levels = slevels,
                      color = :speed, kwargs...)
    s_title = @sprintf("speed √[(u²+w²)/(b⋆H)] @ t=%1.2f", t)

    b_plot = contourf(xb, zb, b_snapshot';
                      clims = (-blim, blim), levels = blevels,
                      color = :thermal, kwargs...)
    b_title = @sprintf("buoyancy, b/b⋆ @ t=%1.2f", t)

    ζ_plot = contourf(xζ, zζ, clamp.(ζ_snapshot', -ζlim, ζlim);
                      clims=(-ζlim, ζlim), levels = ζlevels,
                      color = :balance, kwargs...)
    ζ_title = @sprintf("vorticity, (∂u/∂z - ∂w/∂x) √(H/b⋆) @ t=%1.2f", t)

    χ_plot = contourf(xs, zs, χ_snapshot';
                      clims = (0, χlim), levels = χlevels,
                      color = :dense, kwargs...)
    χ_title = @sprintf("buoyancy dissipation, κ|∇b|² √(H/b⋆⁵) @ t=%1.2f", t)

    plot(s_plot, b_plot, ζ_plot, χ_plot,
           size = (700, 1200),
           link = :x,
         layout = Plots.grid(4, 1),
          title = [s_title b_title ζ_title χ_title])
end

mp4(anim, "horizontal_convection.mp4", fps = 16) # hide

# At higher Rayleigh numbers the flow becomes much more vigorous. See, for example, an animation
# of the voricity of the fluid at ``Ra = 10^{12}`` on [vimeo](https://vimeo.com/573730711).

# ### The Nusselt number
#
# Often we are interested on how much the flow enhances mixing. This is quantified by the
# Nusselt number, which measures how much the flow enhances mixing compared if only diffusion
# was in operation. The Nusselt number is given by
#
# ```math
# Nu = \frac{\langle \chi \rangle}{\langle \chi_{\rm diff} \rangle} \, ,
# ```
#
# where angle brackets above denote both a volume and time average and  ``\chi_{\rm diff}`` is
# the buoyancy dissipation that we get without any flow, i.e., the buoyancy dissipation associated
# with the buoyancy distribution that satisfies
#
# ```math
# \kappa \nabla^2 b_{\rm diff} = 0 \, ,
# ```
#
# with the same boundary conditions same as our setup. In this case we can readily find that
#
# ```math
# b_{\rm diff}(x, z) = b_s(x) \frac{\cosh \left [2 \pi (H + z) / L_x \right ]}{\cosh(2 \pi H / L_x)} \, ,
# ```
# which implies ``\langle \chi_{\rm diff} \rangle = \kappa b_*^2 \pi \tanh(2 \pi Η /Lx)``.
#
# We use the loaded `FieldTimeSeries` to compute the Nusselt number from buoyancy and the volume
# average kinetic energy of the fluid.
#
# First we compute the diffusive buoyancy dissipation, ``\chi_{\rm diff}`` (which is just a
# scalar):

H  = file["grid/Lz"]
Lx = file["grid/Lx"]
κ  = file["closure/κ/b"]

χ_diff = κ * b★^2 * π * tanh(2π * H / Lx)
nothing # hide

# We then create two `ReducedField`s to store the volume integrals of the kinetic energy density
# and the buoyancy dissipation. We need the `grid` to do so; the `grid` can be recoverd from
# any `FieldTimeSeries`, e.g.,

grid = b_timeseries.grid

∫ⱽ_s² = ReducedField(Nothing, Nothing, Nothing, CPU(), grid, dims=(1, 2, 3))
∫ⱽ_mod²_∇b = ReducedField(Nothing, Nothing, Nothing, CPU(), grid, dims=(1, 2, 3))

# We recover the time from the saved `FieldTimeSeries` and construct two empty arrays to store
# the volume-averaged kinetic energy and the instantaneous Nusselt number,

t = b_timeseries.times

kinetic_energy, Nu = zeros(length(t)), zeros(length(t))
nothing # hide

# Now we can loop over the fields in the `FieldTimeSeries`, compute `KineticEnergy` and `Nu`,
# and plot.

for i = 1:length(t)
    s = s_timeseries[i]
    sum!(∫ⱽ_s², s^2 * volume)
    kinetic_energy[i] = 0.5 * ∫ⱽ_s²[1, 1, 1]  / (Lx * H)

    b = b_timeseries[i]
    sum!(∫ⱽ_mod²_∇b, (∂x(b)^2 + ∂z(b)^2) * volume)
    Nu[i] = (κ *  ∫ⱽ_mod²_∇b[1, 1, 1]) / χ_diff
end

p1 = plot(t, kinetic_energy,
          xlabel = "time",
          ylabel = "KE / (b⋆H)",
       linewidth = 3,
          legend = :none)

p2 = plot(t, Nu,
          xlabel = "time",
          ylabel = "Nu",
       linewidth = 3,
          legend = :none)

plot(p1, p2, layout = Plots.grid(2, 1))


text:
const NUM_DIRECTIONS = 4
const RIGHT = 0
const UP = 1
const LEFT = 2
const DOWN = 3

turn_left(direction::Integer) = mod(direction + 1, NUM_DIRECTIONS)
turn_right(direction::Integer) = mod(direction - 1, NUM_DIRECTIONS)

move_up(position::CartesianIndex{2}) = CartesianIndex(position[1] - 1, position[2])
move_down(position::CartesianIndex{2}) = CartesianIndex(position[1] + 1, position[2])
move_left(position::CartesianIndex{2}) = CartesianIndex(position[1], position[2] - 1)
move_right(position::CartesianIndex{2}) = CartesianIndex(position[1], position[2] + 1)

function move_forward(position::CartesianIndex{2}, direction::Integer)
    if direction == UP
        return move_up(position)
    elseif direction == DOWN
        return move_down(position)
    elseif direction == LEFT
        return move_left(position)
    else direction == RIGHT
        return move_right(position)
    end
end

function move_backward(position::CartesianIndex{2}, direction::Integer)
    if direction == UP
        return move_down(position)
    elseif direction == DOWN
        return move_up(position)
    elseif direction == LEFT
        return move_right(position)
    else direction == RIGHT
        return move_left(position)
    end
end

ext: jl


text:
# This file contains some code that is chared across several tests

using NeighbourLists, NBodyIPs, JuLIP, StaticArrays, Test,
      ForwardDiff, Test
using NBodyIPs: NBodyFunction
import NBodyIPs: descriptor
import JuLIP.Potentials: evaluate, cutoff, evaluate_d!

using LinearAlgebra, Printf

# --------- MANY BODY CODE THAT IS SHARED ACROSS TESTS ------------

struct TestBL{N} <: NBodyFunction{N, NBodyIPs.NullSiteDesc}
   r0::Float64
   rcut::Float64
   vN::Val{N}
end

@pot TestBL

function edgelengths1(Rs, J)
   r = [ norm(Rs[j])         for j in J ] |> collect
   S = [ Rs[j] / norm(Rs[j]) for j in J ] |> collect
   for i = 1:length(J)-1, j = i+1:length(J)
      push!(r, norm(Rs[J[j]]-Rs[J[i]]))
      push!(S, (Rs[J[i]]-Rs[J[j]]) / norm(Rs[J[i]]-Rs[J[j]]))
   end
   return r, S
end

# function edgelengths(Rs, J)
#    r = NBodyIPs.edge_lengths(Rs, J)

descriptor(::TestBL) = NBodyIPs.NullSiteDesc()

evaluate(V::TestBL{N}, Rs, ii, J) where {N} =
      fnbody(edgelengths1(Rs, J)[1], V.r0, V.rcut) / N

function evaluate_d!(dVsite, V::TestBL{N}, Rs, ii, J) where {N}
   r, S = edgelengths1(Rs, J)
   dV = (fnbody_d(r, V.r0, V.rcut) / N) .* S
   idx = 0
   for i = 1:length(J)
      idx += 1
      dVsite[J[i]] += dV[idx]
   end
   for i = 1:length(J)-1, j = i+1:length(J)
      idx += 1
      dVsite[J[i]] += dV[idx]
      dVsite[J[j]] -= dV[idx]
   end
   return dVsite
end

evaluate(V::TestBL, rs::Vector{Float64}) = fnbody(rs, V.r0, V.rcut)

cutoff(V::TestBL) = V.rcut

sqrt1p(x) = expm1( 0.5 * log1p( x ) )   # sqrt(1 + x) - 1
sqrt1p_d(x) = 0.5 / sqrt(1+x)

function fnbody(r, r0, rcut)
   E = sum( exp.(1.0 .- r./r0) )
   F = prod( (r./rcut .- 1.0) .* (r .< rcut) )
   return sqrt1p(E) * F^2
end

function fnbody_d(r, r0, rcut)
   E = sum( exp.(1.0 .- r./r0) )
   F = prod( (r./rcut .- 1.0) .* (r .< rcut) )
   ∇f = (- sqrt1p_d(E) * F^2 / r0) .* (exp.(1.0 .- r./r0))
   ∇f += (sqrt1p(E) * 2.0 / rcut * F^2) .* (r./rcut .- 1.0 .- 1e-14).^(-1)
   return ∇f
end

# Generate a MODEL N-Body function
gen_fnbody(N; r0=1.0, rcut = N < 5 ? 2.5 * r0 : 1.7 * r0) =
   TestBL(r0, rcut, Val(N))

println("Checking that `fnbody` is correct...")
fnbody_ad(rs, r0, rcut) = ForwardDiff.gradient( t -> fnbody(t, r0, rcut), rs )
for n = 1:10
   print(".")
   r = SVector(rand(n)...)
   @test fnbody_d(r, 1.0, 2.0) ≈ fnbody_ad(r, 1.0, 2.0)
end
println()

function naive_n_body(X::Vector{JVec{T}}, V::NBodyFunction{N}) where {T, N}
   nX = length(X)
   E = 0.0
   cnt = 0
   fact_N = factorial(N)
   for j in CartesianIndices(ntuple(_->1:nX, N))
      s = zeros(T, (N*(N-1))÷2); n = 0
      for a = 1:N-1, b = a+1:N
         n += 1
         s[n] = norm(X[j[a]] - X[j[b]])
      end
      if 0 < minimum(s) && maximum(s) <= V.rcut
         # each of these terms occur factorial(N) times, i.e. for
         # every permutation of the indices!
         E += V(s) / fact_N
      end
   end
   return E
end

println("--------------------------------------")
println("    Testing EdgeLength Implementation")
println("--------------------------------------")
for N in [2,2,3,3,4,4,5]
   Rs = rand(JVecF, 20)
   J = unique(rand(1:20, 30))
   if length(J) < N
      continue
   end
   J = SVector(J[1:N]...)
   r1, S1 = edgelengths1(Rs, J)
   r = NBodyIPs.edge_lengths(Rs, J)
   (@test r1 ≈ Vector(r)) |> print; print(" ")
   # (@test S1 ≈ Vector(S)) |> print; print(" ")
end
println()

println("--------------------------------------")
println("    Testing NBodyIterator")
println("--------------------------------------")
println("Check that the energy is consistent with a naive implementation")
MM = [2,2,3,3,3,4,4,5]  # body orders
println("   N     Nat    =>   |Emr-Enaive|")
for N in MM
   # create N-body
   VN = gen_fnbody(N; r0 = rnn(:Cu))
   if N < 5
      # create a not-too-large copper cell
      at = bulk(:Cu, cubic=true, pbc = false) * 2
   else
      at = bulk(:Cu, cubic=true, pbc = false) * (1,2,2)
   end

   rattle!(at, 0.1)
   nat = length(at)

   # assemble energy via neighbourlist and map-reduce
   Emr = energy(VN, at)
   # assemble energy naively
   Enaive = naive_n_body(positions(at), VN)

   print("   $N      $nat    =>   $(abs(Emr - Enaive)): ")
   (@test Emr ≈ Enaive) |> println
end

println("----------------------------------------------")
println("Finite-difference tests for forces and stress")
MM = [2,3,4,5] # [2,3,4,5]  # body orders
for N in MM
   # create N-body
   VN = gen_fnbody(N; r0 = rnn(:Cu))
   # create a not-too-large copper cell
   at = bulk(:Cu, cubic=true, pbc = true)
   rattle!(at, 0.03)
   set_constraint!(at, VariableCell(at))
   set_calculator!(at, VN)

   nat = length(positions(at))
   dE = JuLIP.gradient(at)
   dEh = nothing
   E = energy(at)
   println("[N = $N]")
   @printf("    h    | error \n")
   @printf("---------|----------- \n")
   x = dofs(at)
   errors = []
   for p = 2:11
      h = 0.1^p
      dEh = copy(dE)
      for n = 1:length(x)
         x[n] += h
         set_dofs!(at, x)
         dEh[n] = (energy(at) - E) / h
         x[n] -= h
      end
      push!(errors, norm(dE - dEh, Inf))
      @printf(" %1.1e | %4.2e  \n", h, errors[end])
   end
   (@test minimum(errors) <= 1e-3 * maximum(errors)) |> println
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using OnlinePackage
using Documenter: deploydocs, makedocs

makedocs(sitename = "OnlinePackage.jl", modules = [OnlinePackage], doctest = false)
deploydocs(repo = "github.com/bramtayl/OnlinePackage.jl.git")


text:
if OS_NAME == :Darwin
    const unw_getcontext = :unw_getcontext
    const unw_init = :unw_init_local_dwarf
    const unw_step = :unw_step
    const unw_get_reg = :unw_get_reg
    const UNW_REG_IP = -1
else
    const unw_getcontext = :jl_unw_getcontext
    const unw_init = :jl_unw_init_local
    const unw_step = :jl_unw_step
    const unw_get_reg = :jl_unw_get_reg
    const UNW_X86_64_RAX   =  0
    const UNW_X86_64_RDX   =  1
    const UNW_X86_64_RCX   =  2
    const UNW_X86_64_RBX   =  3
    const UNW_X86_64_RSI   =  4
    const UNW_X86_64_RDI   =  5
    const UNW_X86_64_RBP   =  6
    const UNW_X86_64_RSP   =  7
    const UNW_X86_64_R8    =  8
    const UNW_X86_64_R9    =  9
    const UNW_X86_64_R10   = 10
    const UNW_X86_64_R11   = 11
    const UNW_X86_64_R12   = 12
    const UNW_X86_64_R13   = 13
    const UNW_X86_64_R14   = 14
    const UNW_X86_64_R15   = 15
    const UNW_X86_64_RIP   = 16
    const UNW_REG_IP       = UNW_X86_64_RIP

    const UC_MCONTEXT_GREGS_RSP = 0xa0
    const UC_MCONTEXT_GREGS_RIP = 0xa8
end

function get_reg(cursor, reg)
    res = Ref{UInt64}()
    ccall(unw_get_reg, Cint, (Ptr{Void}, Cint, Ref{UInt64}), cursor, reg, res)
    res[]
end
get_ip(cursor) = get_reg(cursor, UNW_REG_IP)

# The first step is hooking specific, since the frame chain has not yet been
# established. It's also very simple since all we need to do is pop the return
# address from the stack and store it as the new IP
function step_first!(RC)
    RC.data[RegisterMap[:rip]] = unsafe_load(convert(Ptr{UInt},RC.data[RegisterMap[:rsp]]))
    RC.data[RegisterMap[:rsp]] += sizeof(Ptr{Void})
end

function rec_backtrace(callback, RC)
    cursor = Array(UInt8, 1000)
    ccall(unw_init, Void, (Ptr{Void}, Ptr{Void}), cursor, RC.data)
    callback(cursor)
    step_first!(RC)
    ccall(unw_init, Void, (Ptr{Void}, Ptr{Void}), cursor, RC.data)
    callback(cursor)
    while ccall(unw_step, Cint, (Ptr{Void},), cursor) > 0
        callback(cursor)
    end
end

function rec_backtrace(RC)
    ips = Array(UInt64, 0)
    rec_backtrace(cursor->push!(ips,get_ip(cursor)), RC)
    ips
end

function local_RC()
    data = Array(UInt, 100)
    ccall(unw_getcontext,Void,(Ptr{Void},),data)
    RegisterContext(data)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libportaudio"
version = v"19.7.0"

# Collection of sources required to build libportaudio. Not all of these
# are used for all platforms.
sources = [
    ArchiveSource("http://files.portaudio.com/archives/pa_stable_v190700_20210406.tgz",
                  "47efbf42c77c19a05d22e627d42873e991ec0c1357219c0d74ce6a2948cb2def"),

    # This includes a patch
    DirectorySource("./bundled"),

    # uncomment the following lines to include ASIO support. To distribute the
    # resulting binaries you'll need to sign the licence agreement included with
    # the SDK at: https://www.steinberg.net/en/company/developers.html

    # FileSource("http://www.steinberg.net/sdk_downloads/ASIOSDK2.3.1.zip",
    #            "31074764475059448a9b7a56f103f4723ed60465e0e9d1a9446ca03dcf840f04")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/portaudio

# move the ASIO SDK to where CMake can find it
if [ -d "asiosdk2.3.1" ]; then
    mv "asiosdk2.3.1 svnrev312937/ASIOSDK2.3.1" asiosdk2.3.1
fi

# apply the patch
atomic_patch -p1 ../patches/win_ds_fix_warning.patch
atomic_patch -p1 ../patches/wasapi-32-bit-windows.patch

# First, build libportaudio
mkdir build-bb
cd build-bb
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DCMAKE_DISABLE_FIND_PACKAGE_PkgConfig=ON \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    ..
make
make install
install_license "../LICENSE.txt"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libportaudio", :libportaudio),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Build against `ALSA` on Linux
    Dependency("alsa_jll"; platforms=filter(Sys.islinux, platforms)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

## Step 7: Final before & after

temp = SimpleSDMPredictor(WorldClim, BioClim, 1);
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0);

## Before
# getindex
l1 = temp[coords];
size(l1)
(l1.bottom, l1.top)
(l1.left, l1.right)
stride(l1)

# geotiff reading
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...);
l2.grid == l1.grid
size(l2)
(l2.bottom, l2.top)
(l2.left, l2.right)
stride(l2)

# geotiff writing
l2_backup = copy(l2);
tempfile = tempname();
geotiff(tempfile, l2);
l3 = replace(geotiff(SimpleSDMPredictor, tempfile), NaN => nothing);
l2_backup.grid == l2.grid
l2_backup.grid == l3.grid

# Inverse bounds
SimpleSDMPredictor(temp.grid, 180.0, -180.0, 90.0, -90.0)

# hcat/vcat
l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0);
l5 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0);
l6 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=20.0, top=30.0);
vcat(l4, l5)
vcat(l4, l6)


text:
cd(@__DIR__); include("setups/simpletree.jl");
plotlyjs(dpi = 200)

## (a)
D = natural_eigdist(𝚽, 𝛌, Q; α = 1.0, input_format = :pmf1, distance = :ROT)
E = transform(fit(MDS, D, maxoutdim=3, distances=true))
E[1, :] .*= 1; E[2, :] .*= -1; E[3, :] .*= 1;
plt = simpletree_mds_plot(E, ieb1, ieb2, ieb3, ieb4, iejc); plot!(camera=(210,-37), xlims = [-18, 15])
# savefig(plt, "../figs/SimpleTree_MDS_ROT1_pmf1_alpha1.html")
savefig(plt, "../figs/SimpleTree_MDS_ROT1_pmf1_alpha1.pdf")

## (b)
D = dist_sROT
E = transform(fit(MDS, D, maxoutdim=3, distances=true))
plt = simpletree_mds_plot(E, ieb1, ieb2, ieb3, ieb4, iejc)
# savefig(plt, "../figs/SimpleTree_MDS_sROT_pmf1_alpha1.html")
savefig(plt, "../figs/SimpleTree_MDS_sROT_pmf1_alpha1.pdf")


text:
module PartitionModel

using DataFrames
using JuMP
using Cbc

export Solve

type PartitionData
  n::Int                                    # Number of entities
  df::DataFrame                             # Dataframe containing entity data
  Ids::DataArray                            # Array of entities
  vNum::Vector{Symbol}                      # Vector containing numeric variables
  vCat::Vector{Symbol}                      # Vector containing categoric variables
  nMet::Dict{Symbol, Dict{String, Float64}} # Numerical metrics
  cMet::Dict{Symbol, Dict{String, Float64}} # Categorical metrics
end

function Solve(fname::String, n_groups::Int, timelimit::Float64=999.9, vars=[])
  return Solve(CreateData(fname, vars), n_groups, timelimit)
end

function GroupSize(m::Int, n::Int)
  k = floor(n / m)
  y = Array(Int, m)
  for g = 1:(n % m)
    y[g] = k + 1
  end
  for g = (n % m + 1) : m
    y[g] = k
  end
  return y
end

function Solve(data::PartitionData, n_groups::Int, timelimit::Float64=999.9)
  groups = GroupSize(n_groups, data.n)

  m = Model(solver=CbcSolver(seconds=timelimit))

  @defVar(m, x[1:length(data.Ids), 1:n_groups], Bin, string("x", i, j))
  @defVar(m, mmin[data.vNum])           # Minimum mean of group
  @defVar(m, mmax[data.vNum])           # Maximum mean of group
  @defVar(m, vmin[data.vNum])           # Minimum variance of group
  @defVar(m, vmax[data.vNum])           # Maximum variance of group
  @defVar(m, cvio[data.vCat] >= 0)           # Violation of categorical constraint

  k1 = 1.0
  k2 = 1.0
  k3 = 1000.0
  @setObjective(m, :Min,
                k1 * sum([(mmax[v] - mmin[v]) / data.nMet[v]["mean"] for v in data.vNum])
                + k2 * sum([(vmax[v] - vmin[v]) / data.nMet[v]["var"] for v in data.vNum])
                + k3 * sum(cvio)
                )

  # Each entity can only be allocated to single group
  @addConstraint(m, c1[i=1:data.n], sum([x[i,j] for j=1:n_groups]) == 1)

  # Set number of people per group
  @addConstraint(m, c2[j=1:n_groups], sum([x[i,j] for i=1:data.n]) == groups[j])

  # Bound numerical variables
  @addConstraint(m, cMeanMax[v=data.vNum, j=1:n_groups], sum([data.df[v][i] * x[i, j] for i=1:data.n]) >= groups[j] * mmin[v])
  @addConstraint(m, cMeanMin[v=data.vNum, j=1:n_groups], sum([data.df[v][i] * x[i, j] for i=1:data.n]) <= groups[j] * mmax[v])
  @addConstraint(m, cVarMax[v=data.vNum, j=1:n_groups], sum([(data.df[v][i] - data.nMet[v]["mean"])^2 * x[i, j] for i=1:data.n]) >= groups[j] * vmin[v])
  @addConstraint(m, cVarMax[v=data.vNum, j=1:n_groups], sum([(data.df[v][i] - data.nMet[v]["mean"])^2 * x[i, j] for i=1:data.n]) <= groups[j] * vmax[v])

  status = solve(m)

  if status == :Optimal || status == :UserLimit
    solution = get_solution(x, data)
    data.df[:Group] = solution
    return data.df
  elseif status == :Unbounded
    error("Unable to solve. Model is unbounded, which means something is wrong with your data")
  elseif status == :Error
    error("Cbc couldn't solve this model.")
  end
end

function get_solution(x, data::PartitionData)
  solution = Array(Int, data.n)
  for i = 1:data.n
    for j = 1:size(x, 2)
      if getValue(x[i, j]) > 0.5
        solution[i] = j
      end
    end
  end
  return solution
end

function CreateData(fname::String, vars::Vector=[])
  df = readtable(fname)
  Ids = df[names(df)[1]]
  vNum = Array(Symbol, 0)
  vCat = Array(Symbol, 0)
  vMet = Dict{Symbol, Dict{String, Float64}} ()
  cMet = Dict{Symbol, Dict{String, Float64}} ()
  for v in names(df)[2:end]
  	if length(vars) > 0
	   	if !(string(v) in vars)
  			# Excluded variable
  			continue
  		end
  	end
    if typeof(df[v]) == DataArray{UTF8String, 1}
      push!(vCat, v)
      cMet[v] = categorical_metrics(df[v])
    else
      push!(vNum, v)
      vMet[v] = numerical_metrics(df[v])
    end
  end
  return PartitionData(length(Ids), df, Ids, vNum, vCat, vMet, cMet)
end

function numerical_metrics(c::DataArray{Float64, 1})
#   return Dict("mean" => mean(c), "var" => var(c))
  return {string("mean") => mean(c), string("var") => var(c)}

end

function categorical_metrics(c::DataArray{UTF8String, 1})
  metrics = Dict{String, Float64}()
  for i in c
    if haskey(metrics, string(i))
      metrics[string(i)] += 1.0
    else
      metrics[string(i)] = 1.0
    end
  end

  for (k, v) in metrics
    metrics[k] /= length(c)
  end
  return metrics
end

end


text:
# Neighbor finders

export
    NoNeighborFinder,
    find_neighbors!,
    DistanceNeighborFinder,
    TreeNeighborFinder

"""
    NoNeighborFinder()

Placeholder neighbor finder that returns no neighbors.
When using this neighbor finder, ensure that `nl_only` for the interactions is
set to `false`.
"""
struct NoNeighborFinder <: NeighborFinder end

"""
    find_neighbors!(simulation, neighbor_finder, step_n; parallel=true)

Obtain a list of close atoms in a system.
Custom neighbor finders should implement this function.
"""
function find_neighbors!(s::Simulation,
                            ::NoNeighborFinder,
                            ::Integer;
                            kwargs...)
    return
end

"""
    DistanceNeighborFinder(nb_matrix, n_steps, dist_cutoff)
    DistanceNeighborFinder(nb_matrix, n_steps)

Find close atoms by distance.
"""
struct DistanceNeighborFinder{T} <: NeighborFinder
    nb_matrix::BitArray{2}
    n_steps::Int
    dist_cutoff::T
end

function DistanceNeighborFinder(nb_matrix::BitArray{2},
                                 n_steps::Integer)
    return DistanceNeighborFinder(nb_matrix, n_steps, 1.2)
end

function find_neighbors!(s::Simulation,
                            nf::DistanceNeighborFinder,
                            step_n::Integer;
                            parallel::Bool=true)
    !iszero(step_n % nf.n_steps) && return

    neighbors = s.neighbors
    empty!(neighbors)
    sqdist_cutoff = nf.dist_cutoff ^ 2

    if parallel && nthreads() > 1
        nl_threads = [Tuple{Int, Int}[] for i in 1:nthreads()]

        @threads for i in 1:length(s.coords)
            nl = nl_threads[threadid()]
            ci = s.coords[i]
            nbi = nf.nb_matrix[:, i]
            for j in 1:(i - 1)
                r2 = sum(abs2, vector(ci, s.coords[j], s.box_size))
                if r2 <= sqdist_cutoff && nbi[j]
                    push!(nl, (i, j))
                end
            end
        end

        for nl in nl_threads
            append!(neighbors, nl)
        end
    else
        for i in 1:length(s.coords)
            ci = s.coords[i]
            nbi = nf.nb_matrix[:, i]
            for j in 1:(i - 1)
                r2 = sum(abs2, vector(ci, s.coords[j], s.box_size))
                if r2 <= sqdist_cutoff && nbi[j]
                    push!(neighbors, (i, j))
                end
            end
        end
    end
end

"""
    TreeNeighborFinder(nb_matrix, n_steps, dist_cutoff)
    TreeNeighborFinder(nb_matrix, n_steps)

Find close atoms by distance (using a tree search).
"""
struct TreeNeighborFinder{T} <: NeighborFinder
    nb_matrix::BitArray{2}
    n_steps::Int
    dist_cutoff::T
end

function TreeNeighborFinder(nb_matrix::BitArray{2},
                                 n_steps::Integer)
    return TreeNeighborFinder(nb_matrix, n_steps, 1.2)
end

function find_neighbors!(s::Simulation,
                          nf::TreeNeighborFinder,
                          step_n::Integer;
                          parallel::Bool=true)
    !iszero(step_n % nf.n_steps) && return

    neighbors = s.neighbors
    empty!(neighbors)

    bv = SVector{3}(s.box_size, s.box_size, s.box_size)
    btree = BallTree(s.coords, PeriodicEuclidean(bv))

    if parallel && nthreads() > 1
        nl_threads = [Tuple{Int, Int}[] for i in 1:nthreads()]

        @threads for i in 1:length(s.coords)
            nl = nl_threads[threadid()]
            ci = s.coords[i]
            nbi = nf.nb_matrix[:, i]
            idxs = inrange(btree, ci, nf.dist_cutoff, true)
            for j in idxs
                if nbi[j] && i > j
                    push!(nl, (i, j))
                end
            end
        end

        for nl in nl_threads
            append!(neighbors, nl)
        end
    else
        for i in 1:length(s.coords)
            ci = s.coords[i]
            nbi = nf.nb_matrix[:, i]
            idxs = inrange(btree, ci, nf.dist_cutoff, true)
            for j in idxs
                if nbi[j] && i > j
                    push!(neighbors, (i, j))
                end
            end
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LCIO_Julia_Wrapper_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LCIO_Julia_Wrapper")
JLLWrappers.@generate_main_file("LCIO_Julia_Wrapper", UUID("1fe59b3e-ae39-56b2-a121-c4c5e6ddab56"))
end  # module LCIO_Julia_Wrapper_jll


text:
using Images, FileIO

ima = load("data/lenna50.jpg")
imb = Gray.(ima)

medcol = median(imb)
imb[imb .≤ medcol] = Gray(0.0)
imb[imb .> medcol] = Gray(1.0)
save("data/lennaGray.jpg", imb)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
ArithmeticBlocks
"""

export MathFunctionBlock

mutable struct MathFunctionBlock <: AbstractArithmeticBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    op::Symbol

    function MathFunctionBlock(op::Symbol)
        @createblock new(Vector{InPort}(), Vector{OutPort}(), op) 1 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::MathFunctionBlock)
    println(io, "MathFunctionBlock($(blk.op))")
end

"""
to expr
"""

function _toexpr(blk::MathFunctionBlock)
    if blk.op == :reciprocal
        Expr(:call, :/, 1, _toexpr(blk.inport[1]))
    elseif blk.op == :square
        Expr(:call, :^, _toexpr(blk.inport[1]), 2)
    end
end


text:
"""
    loglik_drill_lease!(grad, lease, theta, sim, dtv, dograd)

Return log lik of `lease` history. Optionally *adds* to gradient vector `grad`.

Uses temp vector `ubv` from thread-specific `dtv::DrillingTmpVarsThread`
"""
function loglik_drill_lease!(grad, lease, theta, sim, dtv, dograd)

    LL = azero(theta)
    ubv = _ubv(dtv)
    dubv = _dubv(dtv)

    for obs in lease
        actions = actionspace(obs)
        n = length(actions)
        resize!(ubv, n)
        dubv_vw = view(dubv, :, OneTo(n))

        @inbounds @simd for d in actions
            dp1 = d+1
            dubv_slice = view(dubv, :, dp1)
            ubv[dp1] = full_payoff!(dubv_slice, d, obs, theta, sim, dograd)
        end

        dp1 = _y(obs)+1
        LL += ubv[dp1] - logsumexp!(ubv)

        if dograd
            ubv[dp1] -= 1
            BLAS.gemv!('N', -1.0, dubv_vw, ubv, 1.0, grad)
        end
    end

    return LL
end

"""
    loglik_drill_unit!(grad, unit, theta, sim, dtv, dograd)

Return log lik of `unit` history, integrating over `J`. possible leases.

If `dograd`, then `grad .+= ∇(log Lᵢ)` using thread-specific temp variables
`dtv::DrillingTmpVarsThread`. This overwrites `_gradJ(dtv)`
"""
function loglik_drill_unit!(grad, unit, theta, sim, dtv, dograd)

    nJ    = num_initial_leases(unit)
    gradJ = view(_gradJ(dtv), :, OneTo(nJ))
    LLj = _llj(dtv)
    resize!(LLj, nJ)

    fill!(gradJ, 0)

    LL = azero(theta)
    if nJ > 0
        LLj .= log.(j1chars(InitialDrilling(unit)))
        for (ji,lease) in enumerate(InitialDrilling(unit))
            gradj = uview(gradJ, :, ji)
            LLj[ji] += loglik_drill_lease!(gradj, lease, theta, sim, dtv, dograd)
        end
        LL += logsumexp!(LLj)

        dograd && BLAS.gemv!('N', 1.0, gradJ, LLj, 1.0, grad)

    end

    for lease in DevelopmentDrilling(unit)
        LL += loglik_drill_lease!(grad, lease, theta, sim, dtv, dograd)
    end

    return LL
end

"""
    simloglik!(grad, unit::DrillUnit, theta, sims, dtv, dograd)

Threaded computation of `_llm(sims)[m] += log lik(DrillUnit, sims[m])`
"""
function simloglik!(grad, unit::DrillUnit, theta, sims::SimulationDrawsVector, dograd; kwargs...)

    dtv = DrillingTmpVars(unit)
    M = _num_sim(sims)
    llm = _llm(sims)
    gradM = _drillgradm(sims)
    fill!(gradM, 0)

    length(grad) == length(theta) || throw(DimensionMismatch("grad,theta incompatible"))
    size(gradM, 1) == length(grad) || throw(DimensionMismatch("size(gradM) = $(size(gradM)) vs size(theta) = $(size(theta)) not OK"))

    model = _model(_data(unit))
    solve_vf_and_update_itp!(model, theta, ichars(unit), dograd; kwargs...)

    @inbounds for m in OneTo(M)
        sims_m = sims[m]        # get 1 particular simulation
        gradm_i = uview(gradM, :, m)
        llm[m] = loglik_drill_unit!(gradm_i, unit, theta, sims_m, dtv, dograd)
    end
end

function grad_simloglik!(grad, unit::DrillUnit, theta, sims::SimulationDrawsVector)
    M = _num_sim(sims)
    llm = _llm(sims)
    gradM = _drillgradm(sims)

    BLAS.gemv!('N', 1.0, gradM, llm, 1.0, grad)
end

"""
    simloglik_drill_unit!(grad, unit, theta, sim, dtv, dograd)

Compute *simulated* log Lᵢ of `unit` using `sims::SimulationDrawsVector`. If
`dograd`, update cols of `_drillgradm(sims) .+= ∇(log Lᵢₘ)`.

Uses set of thread-specific `dtv::DrillingTmpVarsAll`
"""
function simloglik_drill_unit!(grad, unit, theta, sims, dograd)
    simloglik!(grad, unit, theta, sims, dograd)
    LL = logsumexp!(_llm(sims)) - log(_num_sim(sims))
    dograd && grad_simloglik!(grad, unit, theta, sims)
    return LL
end

function simloglik_drill_data!(grad, hess, data, theta, sim::SimulationDrawsMatrix, dograd=false, dohess=false)

    dohess == true && dograd == false && throw(error("can't dohess without dograd"))
    checksquare(hess) == length(grad) == length(theta) || throw(DimensionMismatch("grad, theta incompatible"))

    θρ = theta_drill_ρ(_model(data), theta)
    update!(sim, θρ)
    g = dohess ? similar(grad) : grad

    LL = azero(theta)
    for (i,unit) in enumerate(data)
        dohess && fill!(g, 0)
        LL += simloglik_drill_unit!(g, unit, theta, view(sim, i), dograd)
        if dohess
            BLAS.ger!(1.0, g, g, hess)
            grad .+= g
        end
    end

    return LL
end


text:
using Sockets
using OpenSSL
using Nghttp2
using Test

function read_all(io::IO)::Vector{UInt8}
    # Create IOBuffer and copy chunks until we read eof.
    result_stream = IOBuffer()

    while !eof(io)
        buffer_chunk = read(io)
        write(result_stream, buffer_chunk)
    end

    seekstart(result_stream)
    return read(result_stream)
end

function read_all_by_byte(io::IO)::Vector{UInt8}
    iob = IOBuffer()
    while !eof(io)
        c = read(io, UInt8)
        write(iob, c)
    end

    seek(iob, 0)

    return read(iob)
end

const DEFAULT_STATUS_200 = [":status" => "200"]
const DEFAULT_TRAILER = ["grpc-status" => "0"]

const DEFAULT_REQUEST_HEADERS = [
    ":method" => "POST",
    ":path" => "/default",
    ":authority" => "localhost:5000",
    ":scheme" => "http",
    "content-type" => "application/text"]

const INVALID_REQUEST_HEADERS = [
    ":method" => "POST",
    ":path" => "/default",
    ":scheme" => "http",
    "content-type" => "application/text",
    ":authority" => "localhost:5000"]

function test_server(socket::Sockets.TCPServer)
    accepted_socket = accept(socket)

    server_session = Nghttp2.from_accepted(accepted_socket)
    @test Nghttp2.is_server_session(server_session.session)

    local request_stream::Http2Stream

    try
        request_stream = recv(server_session)
    catch ex
        close(server_session)
        throw(ex)
    end

    @test isopen(request_stream) == true

    println("==> test_server recv a stream")
    request_data = read_all(request_stream)
    println("test_server, received length: $(length(request_data))")

    send_buffer = IOBuffer(request_data)

    submit_response(request_stream, send_buffer, DEFAULT_STATUS_200)

    close(server_session)

    return close(socket)
end

function test_server()
    socket = listen(5000)
    return test_server(socket)
end

function test_client(request_io::IO=IOBuffer(), headers=DEFAULT_REQUEST_HEADERS)
    tcp_connection = connect(5000)

    # Read all the input data to the buffer, use it later for comparision.
    request_data = read(request_io)
    request_iob = IOBuffer(request_data)

    client_session = Nghttp2.open(tcp_connection)

    stream1 = submit_request(client_session, request_iob, headers)
    response_data = read_all(stream1)

    @test length(response_data) == length(request_data)
    @test response_data == request_data

    return response_data == request_data
end

max_stars_repo_name: grlap/nghttpjl

max_issues_repo_name: grlap/nghttpjl

max_forks_repo_name: grlap/nghttpjl


text:
@inline weighted_sum_mixing(s, v) = matmul(v, s)

@inline mixing(f, v, g, args...) = f(attention_score(g, args...), v)


text:
#=
Run this file to regenerate `geos_c.jl` and `geos_common.jl`.

It expects a GEOS install in the deps folder, run `build LibGEOS` in Pkg mode
if these are not in place.

The wrapped GEOS version and provided GEOS version should be kept in sync.
So when updating the GEOSBuilder provided version, also rerun this wrapper.
This way we ensure that the provided library has the same functions available
as the wrapped one. Furthermore this makes sure constants in `geos_common.jl`
like `GEOS_VERSION`, which are just strings, are correct.
=#

using Clang

includedir = normpath(joinpath(@__DIR__, "..", "deps", "usr", "include"))
headerfiles = [joinpath(includedir, "geos_c.h")]

wc = init(; headers = headerfiles,
            output_file = joinpath(@__DIR__, "geos_c.jl"),
            common_file = joinpath(@__DIR__, "geos_common.jl"),
            clang_includes = [includedir, CLANG_INCLUDE],
            clang_args = ["-I", includedir],
            header_wrapped = (root, current) -> root == current,
            header_library = x -> "libgeos",
            clang_diagnostics = true,
            )

run(wc)


text:
using Musical, DSP
sawwave(x) = x - floor(x)
mywave(x) = 0.5*sin(2*pi*x) + 0.5*sign(sin(2*pi*x*2)) + 0.2*sin(2*pi*x*4+1)+ 0.2*sin(2*pi*x*4.02+1)
#mywave(x) = 1.1*sin(2*pi*x) + 0.3*sign(sin(2*pi*x))+sawwave(x*2)*0.7
expchirpf(x, k::Real=1.1, f0::Real = 1f-1) = x.*(f0.*k.^(x))
mywave(x) = sinewave(expchirpf(x))
env = ADSREnvelope(0.05,0.05,0.9,1,1)
wf = WaveForm(wave=x->sinewave(expchirpf(x)), env=env, nchannels=2, phase=[0, 0.2])

w = waveform(16, 7, 44100, wf)

w2 = waveform(16, 1, 44100, wf)
addto!(w, 6.5, w2)

exportsound(w,"test44b.wav")
#w = importsound("asmz.wav")

#run(`vlc test44b.wav`)
plotall(stereo(w))
#plotall(addto!(w, 0, (-6dB)*waveform(2000, 1, 44100, WaveForm(nchannels=2))))

#f, ax = plotall(stereo(w))
#ax[3][:grid](axis="y", which="major")

wf = WaveForm(env=ENVELOPE_STANDARD)
w = waveform(notefreq("E4"), 5, 44100, wf)+0.5*waveform(notefreq("E5"), 5, 44100, wf)+0.3*waveform(notefreq("E6"), 5, 44100, wf)
w = w/2
w = waveform(notefreq("E2"), 5, 44100, wf)+0.5*waveform(notefreq("E3"), 5, 44100, wf)+0.3*waveform(notefreq("E4"), 5, 44100, wf)
w = w/2
w = waveform(notefreq("A2"), 5, 44100, wf)+0.5*waveform(notefreq("A3"), 5, 44100, wf)+0.3*waveform(notefreq("A4"), 5, 44100, wf)
w = w/2
w = waveform(notefreq("D3"), 5, 44100, wf)+0.5*waveform(notefreq("D4"), 5, 44100, wf)+0.3*waveform(notefreq("D5"), 5, 44100, wf)
w = w/2
w = waveform(notefreq("G3"), 5, 44100, wf)+0.5*waveform(notefreq("G4"), 5, 44100, wf)+0.3*waveform(notefreq("G5"), 5, 44100, wf)
w = w/2
w = waveform(notefreq("B3"), 5, 44100, wf)+0.5*waveform(notefreq("B4"), 5, 44100, wf)+0.3*waveform(notefreq("B5"), 5, 44100, wf)
w = w/2
exportsound(w, "bobo.wav")
wf = WaveForm(wave=sawtooth, env=ENVELOPE_STANDARD, nchannels=2, phase=[0, 0.1])
w = waveform(notefreq("G2"), 5, 44100, wf)+1*waveform(1.01*notefreq("D3"), 5, 44100, wf)+1.3*waveform(notefreq("G3"), 5, 44100, wf)+1*waveform(1.05*notefreq("A#3"), 5, 44100, wf)
w = w
exportsound(w, "bobo.wav")


text:
my_I(T, N) = Matrix{T}(I, N, N)

Zygote._pullback(::AContext, ::typeof(my_I), args...) = my_I(args...), nograd_pullback

function lgssm_components(k::Separable, x::SpaceTimeGrid, storage)

    # Compute spatial covariance, and temporal GaussMarkovModel.
    t = get_times(x)
    r = get_space(x)

    kr, kt = k.l, k.r
    Kr = kernelmatrix(kr, r)
    As_t, as_t, Qs_t, emission_proj_t, x0_t = lgssm_components(kt, t, storage)

    # Compute components of complete LGSSM.
    Nr = length(r)
    ident = my_I(eltype(storage), Nr)
    As = map(A -> kron(ident, A), As_t)
    as = map(a -> repeat(a, Nr), as_t)
    Qs = map(Q -> kron(Kr + ident_eps(1e-12), Q), Qs_t)
    emission_proj = _build_st_proj(emission_proj_t, Nr, ident)
    x0 = Gaussian(repeat(x0_t.m, Nr), kron(Kr, x0_t.P))
    return As, as, Qs, emission_proj, x0
end

function _build_st_proj((Hs, hs)::Tuple{AbstractVector, AbstractVector}, Nr::Integer, ident)
    return (map(H -> kron(ident, H), Hs), map(h -> fill(h, Nr), hs))
end

function build_prediction_obs_vars(
    pr_indices::AbstractVector{<:Integer},
    r_full::AbstractVector{<:AbstractVector},
    σ²s_pr::AbstractVector{<:Diagonal{T}},
) where {T<:Real}
    σ²s_pr_full = map(r -> Diagonal(zeros(T, length(r))), r_full)
    σ²s_pr_full[pr_indices] .= σ²s_pr
    return σ²s_pr_full
end


text:
#import Pkg; Pkg.add("DifferentialEquations")
#import Pkg; Pkg.add("ForwardDiff")
#import Pkg; Pkg.add("DiffResults")
#import Pkg; Pkg.add("DiffEqSensitivity")
#import Pkg; Pkg.add("DataFrames")
#import Pkg; Pkg.add("CSV")
#import Pkg; Pkg.add("Dierckx")
#import Pkg; Pkg.add("StatsPlots")
#import Pkg; Pkg.add("LaTeXStrings")
#import Pkg; Pkg.add("Plots")

using DifferentialEquations
using ForwardDiff
using DiffResults
using DiffEqSensitivity
using Statistics
using Random
using LinearAlgebra
using DataFrames
using CSV
using Dierckx
using Plots
using StatsPlots
using LaTeXStrings

function f_apop(dx, x, p, t)
    x1, x2, x3, x4, x5, x6 = x
    E, tp, ti, td, k, Rm, a1, C1, C2, C4, C5, Ub, U0, Um, Rg, alpha, beta = p
    Vp = 3
    Vi = 11
    Vg = 10
    C3 = 100/100
    meal_t = 300, 650, 1100
    meal_q = 60e3, 40e3, 50e3
    f1 = Rm / (1 + exp(x3 / (Vg * C1) - a1))
    f2 = Ub * (1 - exp(-x3 / (Vg * C2)))
    kappa = (1 / Vi + 1 / (E * ti)) / C4
    f3 = (U0 + Um / (1 + (kappa * x2)^(-beta))) / (Vg * C3)
    f4 = Rg / (1 + exp(alpha * (1 - x6 / (Vp * C5))))
    dt1 = t - meal_t[1]
    dt2 = t - meal_t[2]
    dt3 = t - meal_t[3]
    IG1 = 0.5 * meal_q[1] * k * exp(-k * dt1) * (sign(dt1) + 1)
    IG2 = 0.5 * meal_q[2] * k * exp(-k * dt2) * (sign(dt2) + 1)
    IG3 = 0.5 * meal_q[3] * k * exp(-k * dt3) * (sign(dt3) + 1)
    IG = IG1 + IG2 + IG3
    tmp = E * (x1 / Vp - x2 / Vi)
    dx[1] = (f1 - tmp - x1 / tp)
    dx[2] = (tmp - x2 / ti)
    dx[3] = (f4 + IG - f2 - f3 * x3)
    dx[4] = (x1 - x4) / td
    dx[5] = (x4 - x5) / td
    dx[6] = (x5 - x6) / td
end

cscale = 1; tscale = 1.;
p = [0.200853979,5.986360211,101.2036365,11.97776077,0.00833154,208.6221286,6.592132806,301.2623884,37.65196466,78.75855574,25.93801618,71.32611075,4.063671595,88.97766871,179.85672,7.536059566,1.783340558]

x0 = [36., 44., 11000., 0., 0., 0.] / cscale
tspan = (0.0, 1800.0)
prob_apop = ODELocalSensitivityProblem(f_apop, x0, tspan, p)

sol_apop = solve(prob_apop, alg_hints=[:stiff], saveat=0.1)
x_apop, dp_apop = extract_local_sensitivities(sol_apop)

lab = [L"E", L"tp", L"ti", L"td", L"k", L"Rm", L"a1", L"C1", L"C2", L"C4", L"C5", L"Ub", L"U0", L"Um", L"Rg", L"alpha", L"beta"]
σ = 0.01 * std(x_apop, dims=2)
cov_ϵ = σ[4]
dp = dp_apop
cols = 4:4

#plot(sol_apop.t, x_apop[4,:], lw=2)

Nt = length(dp[1][1,:])
Nstate = length(dp[1][:,1])
Nparam = length(dp[:,1])
F = zeros(Float64, Nparam, Nparam)

perm = vcat(1, sort(rand(2:Nt-1, Nt÷5)), Nt)

for i in perm
    S = reshape(dp[1][:,i], (Nstate,1))
    for j = 2:Nparam
        S = hcat(S, reshape(dp[j][:,i], (Nstate,1)))
    end
    global F += S[cols,:]' * inv(cov_ϵ) * S[cols,:]
end

C = inv(F)
R = ones(size(C))
R = [C[i,j]/sqrt(C[i,i]*C[j,j]) for i = 1:size(C)[1], j = 1:size(C)[1]]
#heatmap(R, xlims=(0.5,size(R)[1]+0.5), aspect_ratio = 1, color = :inferno, clims = (-1, 1),
#        xticks = (1:1:size(C)[1], lab), xtickfont = font(14, "Times"),
#        yticks = (1:1:size(C)[1], lab), ytickfont = font(14, "Times"), fmt = :pdf, dpi=300)
#savefig("correlation_matrix")

abs.(R) .> 0.99

lowerbound = sqrt.(diag(inv(F))) / tscale
lowerbound[1:3:7] = lowerbound[1:3:7] / cscale
for i = 1:length(lab)
    println(lab[i], '\t', lowerbound[i])
end

for i = 1:Nparam
    println(eigvals(F)[i])
    println(eigvecs(F)[:,i])
    println('\n')
end

#bar(eigvecs(F)[:,1:9], ylabel = "FIM null eigenvector coefficients", ytickfont = font(12, "Times"),
#    xticks = (1:1:size(C)[1], lab), xtickfont = font(12, "Courier"),
#    legendfontsize = 10, label = [L"null_1" L"null_2" L"null_3" L"null_4" L"null_5" L"null_6" L"null_7" L"null_8" L"null_9"], fmt = :png,
#    legend=:topright, dpi=300)
#savefig("nulleigen_apop")

max_stars_repo_name: lu-group/sbinn

max_issues_repo_name: lu-group/sbinn

max_forks_repo_name: lu-group/sbinn


text:
@testset "filter" begin

# test filter
sc = SparkContext(master="local")

nums1 = parallelize(sc, 1:30)
n2 = filter(nums1, x->(x>10))
n3 = filter(n2, x->(x<15))

@test count(n2) == 20
@test count(n3) == 4

pnums1 = cartesian(nums1, nums1)
pn2 = filter(pnums1, x->(x[1]>10 && x[2]<10))
pn3 = filter(pnums1, x->(x[1]<15 && x[2]>5))
@test count(pn2) == 180
@test count(pn3) == 350

close(sc)

end


text:
using Pkg
 Pkg.activate(".")

verbose = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectML
 include("shared/scripts.jl")
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics, Dates

include("read_expres_data_101501.jl")

order_list_timeseries = extract_orders(all_spectra,pipeline_plan,recalc=true)

line_list_df = prepare_line_list_pass1(linelist_for_ccf_filename, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true)
 (ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_scale_factor=10.0, ccf_mid_velocity=ccf_mid_velocity, recalc=true)
 line_width = RvSpectML.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.05)
line_list_df = prepare_line_list_pass1(linelist_for_ccf_filename, all_spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = line_width)

(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_scale_factor=2.0, range_no_mask_change=line_width, ccf_mid_velocity=ccf_mid_velocity, v_step=100, recalc=true)
alg_fit_rv = RVFromCCF.MeasureRvFromCCFGaussian(frac_of_width_to_fit=3)
rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)

 if make_plot(pipeline_plan, :ccf_total)
   using Plots
   t_idx = 20
   #plt = plot(v_grid,ccfs[:,t_idx]./maximum(ccfs[:,t_idx],dims=1),label=:none)
   #scatter!(plt,v_grid,ccfs_expr[:,t_idx]./maximum(ccfs_expr[:,t_idx],dims=1),markersize=1.2,label=:none)
   plt = plot(v_grid,ccfs[:,t_idx],label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
   if save_plot(pipeline_plan,:ccf_total)   savefig(plt,joinpath(output_dir,target_subdir * "_ccf_sum.png"))   end
   display(plt)
end

need_to!(pipeline_plan, :template)
if need_to(pipeline_plan, :template)  # Compute order CCF's & measure RVs
   if verbose println("# Making template spectra.")  end
   @assert !need_to(pipeline_plan,:extract_orders)
   GC.gc()   # run garbage collector for deallocated memory
   map(i->order_list_timeseries.metadata[i][:rv_est] = 0.0, 1:length(order_list_timeseries) )
   @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = make_template_spectra(order_list_timeseries)
   if save_data(pipeline_plan, :template)
      using JLD2, FileIO
      save(joinpath(output_dir,target_subdir * "_matrix.jld2"), Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
   dont_need_to!(pipeline_plan, :template);
end

if make_plot(pipeline_plan,:template)
   using Plots
   chunkid = 10
   idx = spectral_orders_matrix.chunk_map[chunkid]
   plt = scatter(spectral_orders_matrix.λ[idx],(f_mean[idx].-1.0)./maximum(abs.((f_mean[idx].-1.0))),markersize=1.0,label=:none)
   scatter!(plt,spectral_orders_matrix.λ[idx],deriv[idx]./maximum(abs.(deriv[idx])),markersize=1.0,label=:none)
   scatter!(plt,spectral_orders_matrix.λ[idx],deriv2[idx]./maximum(abs.(deriv2[idx])),markersize=1.0,label=:none)
   xlabel!("λ (Å)")
   ylabel!("f(λ), f'(λ), f''(λ), all standardized")
   title!("Template spectrum for chunk " * string(chunkid) )
end

need_to!(pipeline_plan,:fit_lines)
if need_to(pipeline_plan,:fit_lines)
   if verbose println("# Performing fresh search for lines in template spectra.")  end
   cl = ChunkList(map(grid->ChunkOfSpectrum(spectral_orders_matrix.λ,f_mean, var_mean, grid), spectral_orders_matrix.chunk_map))
   # We're done with the spectral_orders_matrix, so we can release the memory now
   spectral_orders_matrix = nothing
   GC.gc()
   need_to!(pipeline_plan,:template)
   lines_in_template = LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(min_deriv2=0.5))  # TODO: Automate threshold for finding a line

   if verbose println("# Finding above lines in all spectra.")  end
   @time fits_to_lines = LineFinder.fit_all_lines_in_chunklist_timeseries(order_list_timeseries, lines_in_template )

   if save_data(pipeline_plan,:fit_lines)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_linefinder_lines.csv"), lines_in_template )
      CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits.csv"), fits_to_lines )
      #CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits_clean.csv"), lines_to_try )
   end
   dont_need_to!(pipeline_plan,:fit_lines);
end

using Printf
function select_line_fits_with_good_std_v(line_fits_df::DataFrame, quantile_threshold::Real; verbose::Bool = false, write_csv::Bool = false)
   fit_distrib = line_fits_df |> @groupby(_.line_id) |>
            @map( { median_a=median(_.fit_a), median_b=median(_.fit_b), median_depth=median(_.fit_depth), median_σ²=median(_.fit_σ²), median_λc=median(_.fit_λc),
                   std_a=std(_.fit_a), std_b=std(_.fit_b), std_depth=std(_.fit_depth), std_σ²=std(_.fit_σ²), std_λc=std(_.fit_λc),
                   line_id=first(_.line_id),  frac_converged=mean(_.fit_converged)  } ) |>
            @filter(_.frac_converged == 1.0 ) |> DataFrame

   std_v_threshold =  quantile(fit_distrib.std_λc./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps ,quantile_threshold)
   median_σ_width_treshold = 2000 # quantile(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,1-quantile_threshold)
   good_lines_std_v = fit_distrib |>
            @filter( _.std_λc./_.median_λc.*RvSpectML.speed_of_light_mps <= std_v_threshold ) |>
            @filter( sqrt.(_.median_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps >= median_σ_width_treshold ) |>
                  DataFrame
   if verbose
      println("# Found ", size(good_lines_std_v,1), " good lines (std_v), rejected ", size(fit_distrib,1)-size(good_lines_std_v,1), " lines.")
   end
   if write_csv
      val_str = Printf.@sprintf("%1.2f",quantile_threshold)
      CSV.write(joinpath(output_dir,target_subdir * "_good_lines_stdv_quant=" * val_str * ".csv"), good_lines_std_v )
   end
   return good_lines_std_v
end

function select_line_fits_with_good_depth_width_slope(line_fits_df::DataFrame, quantile_threshold::Real; verbose::Bool = false, write_csv::Bool = false )
   fit_distrib = line_fits_df |> @groupby(_.line_id) |>
            @map( { median_a=median(_.fit_a), median_b=median(_.fit_b), median_depth=median(_.fit_depth), median_σ²=median(_.fit_σ²), median_λc=median(_.fit_λc),
                   std_a=std(_.fit_a), std_b=std(_.fit_b), std_depth=std(_.fit_depth), std_σ²=std(_.fit_σ²), std_λc=std(_.fit_λc),
                   line_id=first(_.line_id),  frac_converged=mean(_.fit_converged)  } ) |>
            @filter(_.frac_converged == 1.0 ) |> DataFrame

   std_depth_treshold = quantile(fit_distrib.std_depth,quantile_threshold)
   median_σ_width_treshold = 2000 # quantile(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,1-quantile_threshold)
   std_σ_width_treshold = quantile(sqrt.(fit_distrib.std_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,quantile_threshold)
   std_b_treshold = quantile(fit_distrib.std_b,quantile_threshold)
   std_a_treshold = quantile(fit_distrib.std_a,quantile_threshold)
   good_lines_alt = fit_distrib |>
      @filter( 0.05 <= _.median_depth <= 0.95 ) |>
      #@filter( sqrt.(_.median_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps >= median_σ_width_treshold ) |>
      @filter( sqrt.(_.std_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps <= std_σ_width_treshold ) |>
      @filter( _.std_b < std_b_treshold) |>

      DataFrame
   if verbose
      println("# Found ", size(good_lines_alt,1), " good lines (std_depth_width_slope), rejected ", size(fit_distrib,1)-size(good_lines_alt,1), " lines.")
   end
   if write_csv
      val_str = Printf.@sprintf("%1.2f",quantile_threshold)
      CSV.write(joinpath(output_dir,target_subdir * "_good_lines_fit_quant=" * val_str * ".csv"), good_lines_alt )
   end
   return good_lines_alt
end

lines_in_template_vs_threshold = map(thresh->LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(min_deriv2=thresh)),[0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.25, 1.5, 2.0])
  map(x->size(x,1),lines_in_template_vs_threshold)
lines_in_template = lines_in_template_vs_threshold[4]
  @time fits_to_lines = LineFinder.fit_all_lines_in_chunklist_timeseries(order_list_timeseries, lines_in_template )

function calc_rvs_from_line_by_line_fits_stdv(fits_to_lines::DataFrame, threshold_stdv::Real)
   good_lines_stdv = select_line_fits_with_good_std_v(fits_to_lines,threshold_stdv,write_csv=false)
   df = fits_to_lines |> @join(good_lines_stdv, _.line_id, _.line_id, {obs_id=_.obs_idx, line_id=_.line_id, Δv=(_.fit_λc.-__.median_λc).*3e8./__.median_λc, weight=(__.median_λc/(3e8*__.std_λc))^2 }) |>
         @groupby(_.obs_id) |> @map({obs_id=first(_.obs_id), Δv= sum(_.Δv.*_.weight)./sum(_.weight) }) |> DataFrame
   rms_stdv = std(df.Δv)
end
function calc_rvs_from_line_by_line_fits_alt(fits_to_lines::DataFrame, threshold_alt::Real)
   good_lines_alt = select_line_fits_with_good_depth_width_slope(fits_to_lines,threshold_alt)
   df2 = fits_to_lines |> @join(good_lines_alt, _.line_id, _.line_id, {obs_id=_.obs_idx, line_id=_.line_id, Δv=(_.fit_λc.-__.median_λc).*3e8./__.median_λc, weight=(__.median_λc/(3e8*__.std_λc))^2 }) |>
            @groupby(_.obs_id) |> @map({obs_id=first(_.obs_id), Δv= sum(_.Δv.*_.weight)./sum(_.weight) }) |> DataFrame
   rms_alt = std(df2.Δv)
end

thresholds = 0.05:0.05:1 #[0.05, 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95]
 plt = plot()
 scatter!(plt,thresholds,map(t->calc_rvs_from_line_by_line_fits_stdv(fits_to_lines,t),thresholds),color=1,label="σ_λ")
 plot!(plt,thresholds,map(t->calc_rvs_from_line_by_line_fits_stdv(fits_to_lines,t),thresholds),color=1,label=:none)
 scatter!(plt,thresholds,map(t->calc_rvs_from_line_by_line_fits_alt(fits_to_lines,t),thresholds),color=2,label="σ_other")
 plot!(plt,thresholds,map(t->calc_rvs_from_line_by_line_fits_alt(fits_to_lines,t),thresholds),color=2,label=:none)
 xlabel!("Threshold for accepting lines")
 ylabel!("RMS v (m/s)")
 savefig("rms_rv_vs_line_accept_threshold.png")

map(c->(typeof(c)),fits_to_lines[1,:])
map(c->(typeof(c)),good_lines_stdv[1,:])
good_lines_stdv

println(names(df))
for row in eachrow(df)
   obs_id = row.obs_id
   good_lines_stdv.line_id .== row.
   mean_λc = row.mean_λ
   #line_weights = 1.0 ./ row.std_λ.^2
   Δvs = map(λ->(λ.-mean_λc).*3e8./mean_λc,row.all_λ)
   ave_Δv = sum(Δvs.*lines_weights) / sum(line_weights)
   println(obs_id,ave_Δv)
end

@map( { obs_id=first(_.obs_id),
        median_a=median(_.fit_a), median_b=median(_.fit_b), median_depth=median(_.fit_depth), median_σ²=median(_.fit_σ²), median_λc=median(_.fit_λc),
               std_a=std(_.fit_a), std_b=std(_.fit_b), std_depth=std(_.fit_depth), std_σ²=std(_.fit_σ²), std_λc=std(_.fit_λc),
               line_id=first(_.line_id),  frac_converged=mean(_.fit_converged)  } ) |>
        @filter(_.frac_converged == 1.0 ) |> DataFrame

if make_plot(pipeline_plan,:template)
   using Plots
   chunkid = 31
   idx = spectral_orders_matrix.chunk_map[chunkid]
   plt = [plot() for i in 1:3]
   for i in 1:3
      local line_plt_idx = findall(x->minimum(spectral_orders_matrix.λ[idx])<=x<=maximum(spectral_orders_matrix.λ[idx]) , good_lines_alt.median_λc)
      map(x->plot!(plt[i],[x,x],[-0.5,0], color=:black,label=:none), good_lines_alt.median_λc[line_plt_idx])
      line_plt_idx = findall(x->minimum(spectral_orders_matrix.λ[idx])<=x<=maximum(spectral_orders_matrix.λ[idx]) , good_lines_stdv.median_λc)
      map(x->plot!(plt[i],[x,x],[ 0,0.5], color=:black,label=:none), good_lines_stdv.median_λc[line_plt_idx])
      scatter!(plt[i],spectral_orders_matrix.λ[idx],(f_mean[idx].-1.0)./maximum(abs.((f_mean[idx].-1.0))),markersize=1.0,color=2,label=:none)
      scatter!(plt[i],spectral_orders_matrix.λ[idx],deriv[idx]./maximum(abs.(deriv[idx])),markersize=1.0,color=3,label=:none)
      scatter!(plt[i],spectral_orders_matrix.λ[idx],deriv2[idx]./maximum(abs.(deriv2[idx])),markersize=1.0,color=4,label=:none)
      ylims!(plt[i],-0.5,0.5)
   end
   λ_div1 = spectral_orders_matrix.λ[idx][floor(Int,length(idx)/3)]
   λ_div2 = spectral_orders_matrix.λ[idx][floor(Int,length(idx)*2/3)]
   xlims!(plt[1],minimum(spectral_orders_matrix.λ[idx]),λ_div1)
   xlims!(plt[2],λ_div1,λ_div2)
   xlims!(plt[3],λ_div2,maximum(spectral_orders_matrix.λ[idx]))
   xlabel!(plt[3],"λ (Å)")
   ylabel!(plt[2],"f(λ), f'(λ), f''(λ), all standardized")
   title!(plt[1],"Template spectrum for chunk " * string(chunkid) )
   plt_combo = plot(plt[1], plt[2], plt[3], layout = (3,1))
   display(plt_combo)
end

# Reset steps of pipeline_plan to rerun with new linelist.
need_to!(pipeline_plan, :ccf_total);
 need_to!(pipeline_plan, :rvs_ccf_total);
 need_to!(pipeline_plan, :ccf_orders);
 need_to!(pipeline_plan, :rvs_ccf_orders);
 need_to!(pipeline_plan, :template);
 need_to!(pipeline_plan, :dcpca);

lines_to_try = good_lines_stdv[5] |> @map({ lambda=_.median_λc, weight=_.median_depth} ) |> DataFrame
 need_to!(pipeline_plan, :ccf_total);
 need_to!(pipeline_plan, :rvs_ccf_total);
 sort!(lines_to_try,:lambda)
 (ccfs2, v_grid2) = ccf_total(order_list_timeseries, lines_to_try, pipeline_plan,  mask_scale_factor=1.0, range_no_mask_change=line_width, ccf_mid_velocity=ccf_mid_velocity, v_step=100, recalc=true)
 alg_fit_rv = RVFromCCF.MeasureRvFromCCFGaussian(frac_of_width_to_fit=3)
 rvs_ccf2 = calc_rvs_from_ccf_total(ccfs2, pipeline_plan, v_grid=v_grid2, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)

lines_to_try = good_lines_alt[4] |> @map({ lambda=_.median_λc, weight=_.median_depth} ) |> DataFrame
  need_to!(pipeline_plan, :ccf_total);
  need_to!(pipeline_plan, :rvs_ccf_total);
  sort!(lines_to_try,:lambda)
  (ccfs2, v_grid2) = ccf_total(order_list_timeseries, lines_to_try, pipeline_plan,  mask_scale_factor=1.0, range_no_mask_change=line_width, ccf_mid_velocity=ccf_mid_velocity, v_step=100, recalc=true)
  alg_fit_rv = RVFromCCF.MeasureRvFromCCFGaussian(frac_of_width_to_fit=3)
  rvs_ccf2 = calc_rvs_from_ccf_total(ccfs2, pipeline_plan, v_grid=v_grid2, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)

names(good_lines_stdv)

if make_plot(pipeline_plan, :ccf_total)
   using Plots
   t_idx = 20
   using Plots
   plt = plot(v_grid,ccfs[:,t_idx]./maximum(ccfs[:,t_idx],dims=1),label=:none)
   scatter!(plt,v_grid2,ccfs2[:,t_idx]./maximum(ccfs2[:,t_idx],dims=1),markersize=1.2,label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
   if save_plot(pipeline_plan,:ccf_total)   savefig(plt,joinpath(output_dir,target_subdir * "_ccf2_sum.png"))   end
   display(plt)
end

# Random plots that may be useful
histogram(log10.(fit_distrib.std_λc./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps ),nbins=100, alpha=0.5, label="All converged fits")
 good_lines = select_line_fits_with_good_depth_width_slope(fits_to_lines,0.8)
 histogram!(log10.(good_lines.std_λc./good_lines.median_λc.*RvSpectML.speed_of_light_mps),nbins=50,alpha=0.5,label="Consistent Fits")
 good_lines = select_line_fits_with_good_std_v(fits_to_lines,0.7)
 histogram!(log10.(good_lines.std_λc./good_lines.median_λc.*RvSpectML.speed_of_light_mps),nbins=50,alpha=0.5, label="Low σ v Fits")

scatter(log10.(fit_distrib.median_depth),(fit_distrib.median_σ²),markersize=2.0)
 scatter!(log10.(bad_lines.median_depth),(bad_lines.median_σ²),markersize=3.0)

quantile(fit_distrib.std_λc./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps , 0.8)

histogram(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,nbins=200)
 histogram!(sqrt.(good_lines.median_σ²)./good_lines.median_λc.*RvSpectML.speed_of_light_mps,nbins=200,alpha=0.5)

quantile(sqrt.(good_lines.median_σ²)./good_lines.median_λc.*RvSpectML.speed_of_light_mps,0.05)

scatter( sqrt.(fit_distrib.std_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,good_lines.std_λc,markersize=2,label=:none)
 scatter!( sqrt.(good_lines.std_σ²)./good_lines.median_λc.*RvSpectML.speed_of_light_mps,good_lines.std_λc,markersize=3,label=:none)

std_depth_treshold = quantile(fit_distrib.std_depth,0.90)
scatter(good_lines.std_depth,good_lines.std_λc,markersize=2,label=:none)

histogram(good_lines.median_depth,nbins=100)

sqrt.(good_lines.median_σ²)./good_lines.median_λc.*RvSpectML.speed_of_light_mps
#scatter(good_lines.median_depth,good_lines.std_λc,markersize=2,label=:none)

scatter(sqrt.(good_lines.median_σ²)./good_lines.median_λc.*RvSpectML.speed_of_light_mps,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.median_depth,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.median_a,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.median_b,good_lines.std_λc,markersize=2,label=:none)

scatter(sqrt.(good_lines.std_σ²)./good_lines.median_λc.*RvSpectML.speed_of_light_mps,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.std_depth,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.std_a,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.std_b,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.std_b,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.std_b,good_lines.std_λc,markersize=2,label=:none)

names(lines_in_template)

println(names(lines_in_template_vs_threshold[1]))


text:
@testset "4.2.7 (d trig)^m (a+b (c cos)^n)^p" begin
    (a, b, n, x, ) = @variables a b n x

    #= ::Package:: =#

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^2)^p*when*a+b=0=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [sin(x)^6/(a - a*cos(x)^2), x, 4, (3*x)/(8*a) - (3*cos(x)*sin(x))/(8*a) - (cos(x)*sin(x)^3)/(4*a)]
    @test_int [sin(x)^5/(a - a*cos(x)^2), x, 3, -(cos(x)/a) + cos(x)^3/(3*a)]
    @test_int [sin(x)^4/(a - a*cos(x)^2), x, 3, x/(2*a) - (cos(x)*sin(x))/(2*a)]
    @test_int [sin(x)^3/(a - a*cos(x)^2), x, 2, -(cos(x)/a)]
    @test_int [sin(x)^2/(a - a*cos(x)^2), x, 2, x/a]
    @test_int [sin(x)^1/(a - a*cos(x)^2), x, 2, -(atanh(cos(x))/a)]
    @test_int [csc(x)^1/(a - a*cos(x)^2), x, 3, -(atanh(cos(x))/(2*a)) - (cot(x)*csc(x))/(2*a)]
    @test_int [csc(x)^2/(a - a*cos(x)^2), x, 3, -(cot(x)/a) - cot(x)^3/(3*a)]
    @test_int [csc(x)^3/(a - a*cos(x)^2), x, 4, -((3*atanh(cos(x)))/(8*a)) - (3*cot(x)*csc(x))/(8*a) - (cot(x)*csc(x)^3)/(4*a)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^2)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [sin(x)^7/(a + b*cos(x)^2), x, 4, -(((a + b)^3*atan((sqrt(b)*cos(x))/sqrt(a)))/(sqrt(a)*b^(7/2))) + ((a^2 + 3*a*b + 3*b^2)*cos(x))/b^3 - ((a + 3*b)*cos(x)^3)/(3*b^2) + cos(x)^5/(5*b)]
    @test_int [sin(x)^5/(a + b*cos(x)^2), x, 4, -(((a + b)^2*atan((sqrt(b)*cos(x))/sqrt(a)))/(sqrt(a)*b^(5/2))) + ((a + 2*b)*cos(x))/b^2 - cos(x)^3/(3*b)]
    @test_int [sin(x)^3/(a + b*cos(x)^2), x, 3, -(((a + b)*atan((sqrt(b)*cos(x))/sqrt(a)))/(sqrt(a)*b^(3/2))) + cos(x)/b]
    @test_int [sin(x)^1/(a + b*cos(x)^2), x, 2, -(atan((sqrt(b)*cos(x))/sqrt(a))/(sqrt(a)*sqrt(b)))]
    @test_int [csc(x)^1/(a + b*cos(x)^2), x, 4, -((sqrt(b)*atan((sqrt(b)*cos(x))/sqrt(a)))/(sqrt(a)*(a + b))) - atanh(cos(x))/(a + b)]
    @test_int [csc(x)^3/(a + b*cos(x)^2), x, 5, -((b^(3/2)*atan((sqrt(b)*cos(x))/sqrt(a)))/(sqrt(a)*(a + b)^2)) - ((a + 3*b)*atanh(cos(x)))/(2*(a + b)^2) - (cot(x)*csc(x))/(2*(a + b))]
    @test_int [csc(x)^5/(a + b*cos(x)^2), x, 6, -((b^(5/2)*atan((sqrt(b)*cos(x))/sqrt(a)))/(sqrt(a)*(a + b)^3)) - ((3*a^2 + 10*a*b + 15*b^2)*atanh(cos(x)))/(8*(a + b)^3) - ((3*a + 7*b)*cot(x)*csc(x))/(8*(a + b)^2) - (cot(x)*csc(x)^3)/(4*(a + b))]

    @test_int [sin(x)^6/(a + b*cos(x)^2), x, 6, -(((8*a^2 + 20*a*b + 15*b^2)*x)/(8*b^3)) - ((a + b)^(5/2)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(sqrt(a)*b^3) + ((4*a + 7*b)*cos(x)*sin(x))/(8*b^2) + (cos(x)*sin(x)^3)/(4*b)]
    @test_int [sin(x)^4/(a + b*cos(x)^2), x, 5, -(((2*a + 3*b)*x)/(2*b^2)) - ((a + b)^(3/2)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(sqrt(a)*b^2) + (cos(x)*sin(x))/(2*b)]
    @test_int [sin(x)^2/(a + b*cos(x)^2), x, 4, -(x/b) - (sqrt(a + b)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(sqrt(a)*b)]
    @test_int [sin(x)^0/(a + b*cos(x)^2), x, 2, -(atan((sqrt(a + b)*cot(x))/sqrt(a))/(sqrt(a)*sqrt(a + b)))]
    @test_int [csc(x)^2/(a + b*cos(x)^2), x, 3, -((b*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(sqrt(a)*(a + b)^(3/2))) - cot(x)/(a + b)]
    @test_int [csc(x)^4/(a + b*cos(x)^2), x, 4, -((b^2*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(sqrt(a)*(a + b)^(5/2))) - ((a + 2*b)*cot(x))/(a + b)^2 - cot(x)^3/(3*(a + b))]
    @test_int [csc(x)^6/(a + b*cos(x)^2), x, 4, -((b^3*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(sqrt(a)*(a + b)^(7/2))) - ((a^2 + 3*a*b + 3*b^2)*cot(x))/(a + b)^3 - ((2*a + 3*b)*cot(x)^3)/(3*(a + b)^2) - cot(x)^5/(5*(a + b))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^3)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^3)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [sin(x)/(4 - 3*cos(x)^3), x, 7, -(atan((1 + 6^(1/3)*cos(x))/sqrt(3))/(2*2^(1/3)*3^(5/6))) + log(2^(2/3) - 3^(1/3)*cos(x))/(6*6^(1/3)) - log(2*2^(1/3) + 2^(2/3)*3^(1/3)*cos(x) + 3^(2/3)*cos(x)^2)/(12*6^(1/3))]

    #= ::Section:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^4)^p=#

    #= ::Section:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^2)^p*when*a+b=0=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/(1 - cos(x)^2), x, 3, -cot(x)]
    @test_int [1/(1 - cos(x)^2)^2, x, 3, -cot(x) - cot(x)^3/3]
    @test_int [1/(1 - cos(x)^2)^3, x, 3, -cot(x) - (2*cot(x)^3)/3 - cot(x)^5/5]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^2)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [cos(x)^7/(a + b*cos(x)^2), x, 4, -((a^3*atanh((sqrt(b)*sin(x))/sqrt(a + b)))/(b^(7/2)*sqrt(a + b))) + ((a^2 - a*b + b^2)*sin(x))/b^3 + ((a - 2*b)*sin(x)^3)/(3*b^2) + sin(x)^5/(5*b)]
    @test_int [cos(x)^5/(a + b*cos(x)^2), x, 4, (a^2*atanh((sqrt(b)*sin(x))/sqrt(a + b)))/(b^(5/2)*sqrt(a + b)) - ((a - b)*sin(x))/b^2 - sin(x)^3/(3*b)]
    @test_int [cos(x)^3/(a + b*cos(x)^2), x, 3, -((a*atanh((sqrt(b)*sin(x))/sqrt(a + b)))/(b^(3/2)*sqrt(a + b))) + sin(x)/b]
    @test_int [cos(x)^1/(a + b*cos(x)^2), x, 2, atanh((sqrt(b)*sin(x))/sqrt(a + b))/(sqrt(b)*sqrt(a + b))]
    @test_int [sec(x)^1/(a + b*cos(x)^2), x, 4, atanh(sin(x))/a - (sqrt(b)*atanh((sqrt(b)*sin(x))/sqrt(a + b)))/(a*sqrt(a + b))]
    @test_int [sec(x)^3/(a + b*cos(x)^2), x, 5, ((a - 2*b)*atanh(sin(x)))/(2*a^2) + (b^(3/2)*atanh((sqrt(b)*sin(x))/sqrt(a + b)))/(a^2*sqrt(a + b)) + (sec(x)*tan(x))/(2*a)]
    @test_int [sec(x)^5/(a + b*cos(x)^2), x, 6, ((3*a^2 - 4*a*b + 8*b^2)*atanh(sin(x)))/(8*a^3) - (b^(5/2)*atanh((sqrt(b)*sin(x))/sqrt(a + b)))/(a^3*sqrt(a + b)) + ((3*a - 4*b)*sec(x)*tan(x))/(8*a^2) + (sec(x)^3*tan(x))/(4*a)]

    @test_int [cos(x)^6/(a + b*cos(x)^2), x, 6, ((8*a^2 - 4*a*b + 3*b^2)*x)/(8*b^3) + (a^(5/2)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(b^3*sqrt(a + b)) - ((4*a - 3*b)*cos(x)*sin(x))/(8*b^2) + (cos(x)^3*sin(x))/(4*b)]
    @test_int [cos(x)^4/(a + b*cos(x)^2), x, 5, -(((2*a - b)*x)/(2*b^2)) - (a^(3/2)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(b^2*sqrt(a + b)) + (cos(x)*sin(x))/(2*b)]
    @test_int [cos(x)^2/(a + b*cos(x)^2), x, 3, x/b + (sqrt(a)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(b*sqrt(a + b))]
    @test_int [cos(x)^0/(a + b*cos(x)^2), x, 2, -(atan((sqrt(a + b)*cot(x))/sqrt(a))/(sqrt(a)*sqrt(a + b)))]
    @test_int [sec(x)^2/(a + b*cos(x)^2), x, 3, (b*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(a^(3/2)*sqrt(a + b)) + tan(x)/a]
    @test_int [sec(x)^4/(a + b*cos(x)^2), x, 4, -((b^2*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(a^(5/2)*sqrt(a + b))) + ((a - b)*tan(x))/a^2 + tan(x)^3/(3*a)]
    @test_int [sec(x)^6/(a + b*cos(x)^2), x, 4, (b^3*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(a^(7/2)*sqrt(a + b)) + ((a^2 - a*b + b^2)*tan(x))/a^3 + ((2*a - b)*tan(x)^3)/(3*a^2) + tan(x)^5/(5*a)]

    @test_int [1/(a + b*cos(x)^2)^2, x, 4, -(((2*a + b)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(2*a^(3/2)*(a + b)^(3/2))) - (b*cos(x)*sin(x))/(2*a*(a + b)*(a + b*cos(x)^2))]

    @test_int [1/(a + b*cos(x)^2)^3, x, 5, -(((8*a^2 + 8*a*b + 3*b^2)*atan((sqrt(a + b)*cot(x))/sqrt(a)))/(8*a^(5/2)*(a + b)^(5/2))) - (b*cos(x)*sin(x))/(4*a*(a + b)*(a + b*cos(x)^2)^2) - (3*b*(2*a + b)*cos(x)*sin(x))/(8*a^2*(a + b)^2*(a + b*cos(x)^2))]

    @test_int [1/(1 + cos(x)^2)^1, x, 2, x/sqrt(2) - atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2))/sqrt(2)]
    @test_int [1/(1 + cos(x)^2)^2, x, 4, (3*x)/(4*sqrt(2)) - (3*atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2)))/(4*sqrt(2)) - (cos(x)*sin(x))/(4*(1 + cos(x)^2))]
    @test_int [1/(1 + cos(x)^2)^3, x, 5, (19*x)/(32*sqrt(2)) - (19*atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2)))/(32*sqrt(2)) - (cos(x)*sin(x))/(8*(1 + cos(x)^2)^2) - (9*cos(x)*sin(x))/(32*(1 + cos(x)^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^2)^(p/2)*when*a+b=0=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sqrt(1 - cos(x)^2), x, 3, -(cot(x)*sqrt(sin(x)^2))]
    @test_int [sqrt(-1 + cos(x)^2), x, 3, (-cot(x))*sqrt(-sin(x)^2)]

    @test_int [(1 - cos(x)^2)^(3/2), x, 4, (-(2/3))*cot(x)*sqrt(sin(x)^2) - (1/3)*cot(x)*(sin(x)^2)^(3/2)]
    @test_int [(-1 + cos(x)^2)^(3/2), x, 4, (2/3)*cot(x)*sqrt(-sin(x)^2) - (1/3)*cot(x)*(-sin(x)^2)^(3/2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/sqrt(1 - cos(x)^2), x, 3, -((atanh(cos(x))*sin(x))/sqrt(sin(x)^2))]
    @test_int [1/sqrt(-1 + cos(x)^2), x, 3, -((atanh(cos(x))*sin(x))/sqrt(-sin(x)^2))]

    @test_int [1/(1 - cos(x)^2)^(3/2), x, 4, -(cot(x)/(2*sqrt(sin(x)^2))) - (atanh(cos(x))*sin(x))/(2*sqrt(sin(x)^2))]
    @test_int [1/(-1 + cos(x)^2)^(3/2), x, 4, cot(x)/(2*sqrt(-sin(x)^2)) + (atanh(cos(x))*sin(x))/(2*sqrt(-sin(x)^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sqrt(1 + cos(x)^2), x, 1, Elliptic.E(pi/2 + x, -1)]
    @test_int [sqrt(-1 - cos(x)^2), x, 2, (sqrt(-1 - cos(x)^2)*Elliptic.E(pi/2 + x, -1))/sqrt(1 + cos(x)^2)]
    @test_int [sqrt(a + b*cos(x)^2), x, 2, (sqrt(a + b*cos(x)^2)*Elliptic.E(pi/2 + x, -(b/a)))/sqrt(1 + (b*cos(x)^2)/a)]

    @test_int [(1 + cos(x)^2)^(3/2), x, 4, 2*Elliptic.E(pi/2 + x, -1) - (2/3)*Elliptic.F(pi/2 + x, -1) + (1/3)*cos(x)*sqrt(1 + cos(x)^2)*sin(x)]
    @test_int [(-1 - cos(x)^2)^(3/2), x, 6, -((2*sqrt(-1 - cos(x)^2)*Elliptic.E(pi/2 + x, -1))/sqrt(1 + cos(x)^2)) - (2*sqrt(1 + cos(x)^2)*Elliptic.F(pi/2 + x, -1))/(3*sqrt(-1 - cos(x)^2)) - (1/3)*cos(x)*sqrt(-1 - cos(x)^2)*sin(x)]
    @test_int [(a + b*cos(x)^2)^(3/2), x, 6, (2*(2*a + b)*sqrt(a + b*cos(x)^2)*Elliptic.E(pi/2 + x, -(b/a)))/(3*sqrt(1 + (b*cos(x)^2)/a)) - (a*(a + b)*sqrt(1 + (b*cos(x)^2)/a)*Elliptic.F(pi/2 + x, -(b/a)))/(3*sqrt(a + b*cos(x)^2)) + (1/3)*b*cos(x)*sqrt(a + b*cos(x)^2)*sin(x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/sqrt(1 + cos(x)^2), x, 1, Elliptic.F(pi/2 + x, -1)]
    @test_int [1/sqrt(-1 - cos(x)^2), x, 2, (sqrt(1 + cos(x)^2)*Elliptic.F(pi/2 + x, -1))/sqrt(-1 - cos(x)^2)]
    @test_int [1/sqrt(a + b*cos(x)^2), x, 2, (sqrt(1 + (b*cos(x)^2)/a)*Elliptic.F(pi/2 + x, -(b/a)))/sqrt(a + b*cos(x)^2)]

    @test_int [1/(1 + cos(x)^2)^(3/2), x, 3, (1/2)*Elliptic.E(pi/2 + x, -1) - (cos(x)*sin(x))/(2*sqrt(1 + cos(x)^2))]
    @test_int [1/(-1 - cos(x)^2)^(3/2), x, 4, (sqrt(-1 - cos(x)^2)*Elliptic.E(pi/2 + x, -1))/(2*sqrt(1 + cos(x)^2)) + (cos(x)*sin(x))/(2*sqrt(-1 - cos(x)^2))]
    @test_int [1/(a + b*cos(x)^2)^(3/2), x, 4, (sqrt(a + b*cos(x)^2)*Elliptic.E(pi/2 + x, -(b/a)))/(a*(a + b)*sqrt(1 + (b*cos(x)^2)/a)) - (b*cos(x)*sin(x))/(a*(a + b)*sqrt(a + b*cos(x)^2))]

    @test_int [cos(x)/sqrt(1 + cos(x)^2), x, 2, asin(sin(x)/sqrt(2))]
    @test_int [cos(5 + 3*x)/sqrt(3 + cos(5 + 3*x)^2), x, 2, asin(sin(5 + 3*x)/2)/3]
    @test_int [cos(x)/sqrt(4 - cos(x)^2), x, 2, asinh(sin(x)/sqrt(3))]

    #= ::Section:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^3)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^4)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^4)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/(a + b*cos(x)^4), x, 10, ((sqrt(a) + sqrt(a + b))*atan((a^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b)) - sqrt(2)*(a + b)^(3/4)*cot(x))/(a^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b)))))/(2*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b))) - ((sqrt(a) + sqrt(a + b))*atan((a^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b)) + sqrt(2)*(a + b)^(3/4)*cot(x))/(a^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b)))))/(2*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b + sqrt(a)*sqrt(a + b))) - ((sqrt(a) - sqrt(a + b))*log(sqrt(a)*(a + b)^(1/4) - sqrt(2)*a^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b))*cot(x) + (a + b)^(3/4)*cot(x)^2))/(4*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b))) + ((sqrt(a) - sqrt(a + b))*log(sqrt(a)*(a + b)^(1/4) + sqrt(2)*a^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b))*cot(x) + (a + b)^(3/4)*cot(x)^2))/(4*sqrt(2)*a^(3/4)*(a + b)^(1/4)*sqrt(a + b - sqrt(a)*sqrt(a + b)))]

    @test_int [1/(a - b*cos(x)^4), x, 4, -(atan((sqrt(sqrt(a) - sqrt(b))*cot(x))/a^(1/4))/(2*a^(3/4)*sqrt(sqrt(a) - sqrt(b)))) - atan((sqrt(sqrt(a) + sqrt(b))*cot(x))/a^(1/4))/(2*a^(3/4)*sqrt(sqrt(a) + sqrt(b)))]

    @test_int [1/(1 + cos(x)^4), x, 10, x/(2*sqrt(-1 + sqrt(2))) + atan(((-2 + sqrt(2))*cos(x)*sin(x) + sqrt(-1 + sqrt(2))*(1 - 2*sin(x)^2))/(2 + sqrt(1 + sqrt(2)) + 2*sqrt(-1 + sqrt(2))*cos(x)*sin(x) + (-2 + sqrt(2))*sin(x)^2))/(4*sqrt(-1 + sqrt(2))) + atan(((-2 + sqrt(2))*cos(x)*sin(x) + sqrt(-1 + sqrt(2))*(-1 + 2*sin(x)^2))/(2 + sqrt(1 + sqrt(2)) - 2*sqrt(-1 + sqrt(2))*cos(x)*sin(x) + (-2 + sqrt(2))*sin(x)^2))/(4*sqrt(-1 + sqrt(2))) + (1/8)*sqrt(-1 + sqrt(2))*log(sqrt(2) - 2*sqrt(-1 + sqrt(2))*cot(x) + 2*cot(x)^2) - (1/8)*sqrt(-1 + sqrt(2))*log(1 + sqrt(2*(-1 + sqrt(2)))*cot(x) + sqrt(2)*cot(x)^2)]

    @test_int [1/(1 - cos(x)^4), x, 3, x/(2*sqrt(2)) - atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2))/(2*sqrt(2)) - cot(x)/2]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cos(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    @test_int [1/(a + b*cos(x)^5), x, 12, (2*atan((sqrt(a^(1/5) - b^(1/5))*tan(x/2))/sqrt(a^(1/5) + b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - b^(1/5))*sqrt(a^(1/5) + b^(1/5))) + (2*atan((sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))) + (2*atan((sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))) + (2*atan((sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))) + (2*atan((sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5)))]
    @test_int [1/(a + b*cos(x)^6), x, 7, -(atan((sqrt(a^(1/3) + b^(1/3))*cot(x))/a^(1/6))/(3*a^(5/6)*sqrt(a^(1/3) + b^(1/3)))) - atan((sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))*cot(x))/a^(1/6))/(3*a^(5/6)*sqrt(a^(1/3) - (-1)^(1/3)*b^(1/3))) - atan((sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3))*cot(x))/a^(1/6))/(3*a^(5/6)*sqrt(a^(1/3) + (-1)^(2/3)*b^(1/3)))]
    @test_int [1/(a + b*cos(x)^8), x, 9, atan((sqrt((-a)^(1/4) - b^(1/4))*cot(x))/(-a)^(1/8))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) - b^(1/4))) + atan((sqrt((-a)^(1/4) - im*b^(1/4))*cot(x))/(-a)^(1/8))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) - im*b^(1/4))) + atan((sqrt((-a)^(1/4) + im*b^(1/4))*cot(x))/(-a)^(1/8))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) + im*b^(1/4))) + atan((sqrt((-a)^(1/4) + b^(1/4))*cot(x))/(-a)^(1/8))/(4*(-a)^(7/8)*sqrt((-a)^(1/4) + b^(1/4)))]

    @test_int [1/(a - b*cos(x)^5), x, 12, (2*atan((sqrt(a^(1/5) + b^(1/5))*tan(x/2))/sqrt(a^(1/5) - b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - b^(1/5))*sqrt(a^(1/5) + b^(1/5))) + (2*atan((sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(1/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(1/5)*b^(1/5))) + (2*atan((sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(2/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(2/5)*b^(1/5))) + (2*atan((sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(3/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(3/5)*b^(1/5))) + (2*atan((sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5))*tan(x/2))/sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))))/(5*a^(4/5)*sqrt(a^(1/5) - (-1)^(4/5)*b^(1/5))*sqrt(a^(1/5) + (-1)^(4/5)*b^(1/5)))]
    @test_int [1/(a - b*cos(x)^6), x, 7, -(atan((sqrt(a^(1/3) - b^(1/3))*cot(x))/a^(1/6))/(3*a^(5/6)*sqrt(a^(1/3) - b^(1/3)))) - atan((sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))*cot(x))/a^(1/6))/(3*a^(5/6)*sqrt(a^(1/3) + (-1)^(1/3)*b^(1/3))) - atan((sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3))*cot(x))/a^(1/6))/(3*a^(5/6)*sqrt(a^(1/3) - (-1)^(2/3)*b^(1/3)))]
    @test_int [1/(a - b*cos(x)^8), x, 9, -(atan((sqrt(a^(1/4) - b^(1/4))*cot(x))/a^(1/8))/(4*a^(7/8)*sqrt(a^(1/4) - b^(1/4)))) - atan((sqrt(a^(1/4) - im*b^(1/4))*cot(x))/a^(1/8))/(4*a^(7/8)*sqrt(a^(1/4) - im*b^(1/4))) - atan((sqrt(a^(1/4) + im*b^(1/4))*cot(x))/a^(1/8))/(4*a^(7/8)*sqrt(a^(1/4) + im*b^(1/4))) - atan((sqrt(a^(1/4) + b^(1/4))*cot(x))/a^(1/8))/(4*a^(7/8)*sqrt(a^(1/4) + b^(1/4)))]

    @test_int [1/(1 + cos(x)^5), x, 11, (2*atan(sqrt((1 - (-1)^(2/5))/(1 + (-1)^(2/5)))*tan(x/2)))/(5*sqrt(1 - (-1)^(4/5))) + (2*atan(sqrt((1 - (-1)^(4/5))/(1 + (-1)^(4/5)))*tan(x/2)))/(5*sqrt(1 + (-1)^(3/5))) - (2*atanh(tan(x/2)/sqrt(-((1 - (-1)^(1/5))/(1 + (-1)^(1/5))))))/(5*sqrt(-1 + (-1)^(2/5))) - (2*sqrt(-((1 + (-1)^(3/5))/(1 - (-1)^(3/5))))*atanh(sqrt(-((1 + (-1)^(3/5))/(1 - (-1)^(3/5))))*tan(x/2)))/(5*(1 + (-1)^(3/5))) + sin(x)/(5*(1 + cos(x)))]
    @test_int [1/(1 + cos(x)^6), x, 7, atan(tan(x)/sqrt(2))/(3*sqrt(2)) + atan(tan(x)/sqrt(1 - (-1)^(1/3)))/(3*sqrt(1 - (-1)^(1/3))) + atan(tan(x)/sqrt(1 + (-1)^(2/3)))/(3*sqrt(1 + (-1)^(2/3))), x/(3*sqrt(2)) - atan(sqrt(1 - (-1)^(1/3))*cot(x))/(3*sqrt(1 - (-1)^(1/3))) - atan(sqrt(1 + (-1)^(2/3))*cot(x))/(3*sqrt(1 + (-1)^(2/3))) - atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2))/(3*sqrt(2))]
    @test_int [1/(1 + cos(x)^8), x, 9, -(atan(sqrt(1 - (-1)^(1/4))*cot(x))/(4*sqrt(1 - (-1)^(1/4)))) - atan(sqrt(1 + (-1)^(1/4))*cot(x))/(4*sqrt(1 + (-1)^(1/4))) - atan(sqrt(1 - (-1)^(3/4))*cot(x))/(4*sqrt(1 - (-1)^(3/4))) - atan(sqrt(1 + (-1)^(3/4))*cot(x))/(4*sqrt(1 + (-1)^(3/4)))]

    @test_int [1/(1 - cos(x)^5), x, 11, (2*atan(sqrt((1 - (-1)^(1/5))/(1 + (-1)^(1/5)))*tan(x/2)))/(5*sqrt(1 - (-1)^(2/5))) + (2*atan(sqrt((1 - (-1)^(3/5))/(1 + (-1)^(3/5)))*tan(x/2)))/(5*sqrt(1 + (-1)^(1/5))) - (2*atanh(tan(x/2)/sqrt(-((1 - (-1)^(2/5))/(1 + (-1)^(2/5))))))/(5*sqrt(-1 + (-1)^(4/5))) + (2*atanh(sqrt(-((1 + (-1)^(4/5))/(1 - (-1)^(4/5))))*tan(x/2)))/(5*sqrt(-1 - (-1)^(3/5))) - sin(x)/(5*(1 - cos(x)))]
    @test_int [1/(1 - cos(x)^6), x, 8, -(atan(sqrt(1 + (-1)^(1/3))*cot(x))/(3*sqrt(1 + (-1)^(1/3)))) - atan(sqrt(1 - (-1)^(2/3))*cot(x))/(3*sqrt(1 - (-1)^(2/3))) - cot(x)/3]
    @test_int [1/(1 - cos(x)^8), x, 10, x/(4*sqrt(2)) - atan(sqrt(1 - im)*cot(x))/(4*sqrt(1 - im)) - atan(sqrt(1 + im)*cot(x))/(4*sqrt(1 + im)) - atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2))/(4*sqrt(2)) - cot(x)/4]

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^2)^p=#

    @test_int [tan(x)/(1 + cos(x)^2), x, 4, -log(cos(x)) + (1/2)*log(1 + cos(x)^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [tan(x)*sqrt(a + b*cos(x)^2), x, 4, sqrt(a)*atanh(sqrt(a + b*cos(x)^2)/sqrt(a)) - sqrt(a + b*cos(x)^2)]

    @test_int [tan(x)*sqrt(1 - cos(x)^2), x, 5, atanh(sqrt(sin(x)^2)) - sqrt(sin(x)^2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tan(x)/sqrt(a + b*cos(x)^2), x, 3, atanh(sqrt(a + b*cos(x)^2)/sqrt(a))/sqrt(a)]
    @test_int [tan(x)/sqrt(1 + cos(x)^2), x, 3, atanh(sqrt(1 + cos(x)^2))]

    @test_int [tan(x)/sqrt(1 - cos(x)^2), x, 4, atanh(sqrt(sin(x)^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^3)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^3)^p=#

    @test_int [tan(x)^3/(a + b*cos(x)^3), x, 11, -((b^(2/3)*atan((a^(1/3) - 2*b^(1/3)*cos(x))/(sqrt(3)*a^(1/3))))/(sqrt(3)*a^(5/3))) + log(cos(x))/a + (b^(2/3)*log(a^(1/3) + b^(1/3)*cos(x)))/(3*a^(5/3)) - (b^(2/3)*log(a^(2/3) - a^(1/3)*b^(1/3)*cos(x) + b^(2/3)*cos(x)^2))/(6*a^(5/3)) - log(a + b*cos(x)^3)/(3*a) + sec(x)^2/(2*a)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^3)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [tan(x)*sqrt(a + b*cos(x)^3), x, 5, (2/3)*sqrt(a)*atanh(sqrt(a + b*cos(x)^3)/sqrt(a)) - (2/3)*sqrt(a + b*cos(x)^3)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tan(x)/sqrt(a + b*cos(x)^3), x, 4, (2*atanh(sqrt(a + b*cos(x)^3)/sqrt(a)))/(3*sqrt(a))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^4)^p=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^4)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^4)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [tan(x)*sqrt(a + b*cos(x)^4), x, 5, (1/2)*sqrt(a)*atanh(sqrt(a + b*cos(x)^4)/sqrt(a)) - (1/2)*sqrt(a + b*cos(x)^4)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tan(x)/sqrt(a + b*cos(x)^4), x, 4, (2*atanh(sqrt(a + b*cos(x)^4)/sqrt(a)))/(4*sqrt(a))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^n)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^m*(a+b*cos(e+f*x)^n)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [tan(x)*sqrt(a + b*cos(x)^n), x, 5, (2*sqrt(a)*atanh(sqrt(a + b*cos(x)^n)/sqrt(a)))/n - (2*sqrt(a + b*cos(x)^n))/n]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tan(x)/sqrt(a + b*cos(x)^n), x, 4, (2*atanh(sqrt(a + b*cos(x)^n)/sqrt(a)))/(sqrt(a)*n)]
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.7 (d trig)^m (a+b (c cos)^n)^p.jl


text:
using Colors, Highlights
using Highlights.Format
using Highlights.Tokens
using Highlights.Themes

if !isdefined(Highlights.Themes, :has_fg)
    has_fg(style) = style.fg.active
    css2color(c) = RGBA{Float32}(c.r/255, c.g/255, c.b/255, 1)
else
    import Highlights.Themes: has_fg
    css2color(str) = parse(RGBA{Float32}, string("#", str))
end

function style2color(style, default)
    if has_fg(style)
        css2color(style.fg)
    else
        default
    end
end

function render_str(
        ctx::Format.Context, theme::Format.Theme
    )
    defaultcolor = if has_fg(theme.base)
        css2color(theme.base.fg)
    else
        RGBA(0f0, 0f0, 0f0, 1f0)
    end
    colormap = map(s-> style2color(s, defaultcolor), theme.styles)
    tocolor = Dict(zip(Tokens.__TOKENS__, colormap))
    colors = RGBA{Float32}[]
    io = IOBuffer()
    for token in ctx.tokens
        t = Tokens.__TOKENS__[token.value.value]
        str = SubString(ctx.source, token.first, token.last)
        print(io, str)
        append!(colors, fill(tocolor[t], length(str)))
    end
    String(take!(io)), colors
end
function highlight_text(src::AbstractString, T = Themes.DefaultTheme)
    io = IOBuffer()
    render_str(
        Highlights.Compiler.lex(src, Lexers.JuliaLexer),
        Themes.theme(T)
    )
end


text:
using LabConnections
using Base.Test

# This flag is enabled if a dummy filesystem should be used for testing (for online testing)
# disabling the flag allows the BBB to be run in the loop, in this case blinking LEDS
RUNNING_TESTS = true

include("BeagleBone/SYS_LED_test.jl")
include("BeagleBone/PWM_test.jl")
include("BeagleBone/GPIO_test.jl")


text:
#=
Given the root to a binary search tree, find the second largest node in the tree.
=#

struct Node
    key::Int
    left::Union{Node, Nothing}
    right::Union{Node, Nothing}
end

function node_size(node::Node)
    size = 1
    if !isnothing(node.left)
        size = size + node_size(node.left)
    end
    if !isnothing(node.right)
        size = size + node_size(node.right)
    end
    return size
end

function secondlargestnode(node::Node)
    if node_size(node.right) >= node_size(node.left)
        return node.right
    else
        return node.left
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
function test_rabit()
    rabit_init()
    @test rabit_get_rank() == 0
    @test rabit_get_world_size() == 1
    @test !rabit_is_distributed()
    @test rabit_get_version_number() == 0
    rabit_finalize()
end

test_rabit()


text:
Base.:+(p::Pol) = p
Base.:-(p::Pol) = Pol([-c for c in p.coeffs])

function Base.:+(p1::Pol, p2::Pol)
    maxdeg = max(deg(p1), deg(p2))
    coeffs = [p1[i]+p2[i] for i in 1:maxdeg+1]
    return Pol(coeffs)
end

Base.:+(p1::Union{Pol, Number}, p2::Union{Pol, Number}) = +(promote(p1, p2)...)
Base.:-(p1::Union{Pol, Number}, p2::Union{Pol, Number}) = p1 + (-p2)

function Base.:*(p1::Pol{T}, p2::Pol{T}) where {T<:Number}
    coeffs = zeros(T, deg(p1)+deg(p2)+1)
    @inbounds for (i, c1) in enumerate(p1.coeffs)
        @inbounds for (j, c2) in enumerate(p2.coeffs)
            coeffs[i+j-1] += c1*c2
        end
    end
    return Pol(coeffs)
end

Base.:*(p1::Union{Pol, Number}, p2::Union{Pol, Number}) = *(promote(p1, p2)...)

Base.:/(p::Pol, n::Number)  = Pol([c/n for c in p.coeffs])

function Base.:^(p::Pol, n::Integer)
    iszero(n) && return Pol([1])
    if iseven(n)
        return (p*p)^(n÷2)
    else
        return p*(p*p)^((n-1)÷2)
    end
end

function Base.divrem(a::Pol, b::Pol)
    r = a
    q = zero(a)
    d = deg(b)
    c = lc(b)
    while !iszero(r) && deg(r) >= d
        s = lc(r)/c
        coeffs = zeros(typeof(s), deg(r)-d+1)
        coeffs[end] = s
        s = Pol(coeffs)
        q, r = q + s, r - s*b
    end
    return q, r
end

function Base.div(a::Pol, b::Pol)
    q, _ = divrem(a, b)
    return q
end

function Base.rem(a::Pol, b::Pol)
    _, r = divrem(a, b)
    return r
end

function Base.gcd(a::Pol, b::Pol)
    r0 = a
    r1 = b
    while !iszero(r1)
        r1, r0 = rem(r0, r1), r1
    end
    return r0/lc(r0)
end

# comparisons

function Base.:(==)(p1::Pol, p2::Pol)
    deg(p1) == deg(p2) || return false
    @inbounds for i in 1:length(p1.coeffs)
        p1.coeffs[i] == p2.coeffs[i] || return false
    end
    return true
end

for op in (:>, :<, :>=, :<=)
    @eval Base.$op(p1::Pol, p2::Pol) = $op(deg(p1), deg(p2))
end

Base.iszero(p::Pol) = iszero(p.coeffs)