text:
# # Using MLJ with Flux to train the iris dataset

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# **Julia version** is assumed to be 1.6.*

using MLJ
using Flux
import RDatasets

using Random
Random.seed!(123)

MLJ.color_off()

using Plots
pyplot(size=(600, 300*(sqrt(5)-1)));

# Following is a very basic introductory example, using a default
# builder and no standardization of input features.

# For a more advanced illustration, see the [MNIST images
# example](https://github.com/FluxML/MLJFlux.jl/blob/dev/examples/mnist).

# ## Loading some data and instantiating a model

iris = RDatasets.dataset("datasets", "iris");
y, X = unpack(iris, ==(:Species), colname -> true, rng=123);
NeuralNetworkClassifier = @load NeuralNetworkClassifier
clf = NeuralNetworkClassifier()

# ## Incremental training

import Random.seed!; seed!(123)
mach = machine(clf, X, y)
fit!(mach)

training_loss = cross_entropy(predict(mach, X), y) |> mean

# Increasing learning rate and adding iterations:

clf.optimiser.eta = clf.optimiser.eta * 2
clf.epochs = clf.epochs + 5

fit!(mach, verbosity=2);

#-

training_loss = cross_entropy(predict(mach, X), y) |> mean

# ## Accessing the Flux chain (model)

chain = fitted_params(mach).chain

# ##  Evolution of out-of-sample performance

r = range(clf, :epochs, lower=1, upper=200, scale=:log10)
curve = learning_curve(clf, X, y,
                       range=r,
                       resampling=Holdout(fraction_train=0.7),
                       measure=cross_entropy)
using Plots
plot(curve.parameter_values,
       curve.measurements,
       xlab=curve.parameter_name,
       xscale=curve.parameter_scale,
       ylab = "Cross Entropy")

#-

savefig("iris_history.png")

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src


text:
using SketchyCGAL
using Test

include("maxcut.jl")

@testset "maxcut examples" begin
	maxcut_example_1()
	maxcut_example_2()
	maxcut_example_3()
end

ext: jl


text:
using BinaryBuilder

include("../common.jl")

name = "Go"
version = v"1.16.3"

# https://golang.org/dl/
sources = [
    ArchiveSource(
        "https://golang.org/dl/go1.16.3.linux-amd64.tar.gz",
        "951a3c7c6ce4e56ad883f97d9db74d3d6d80d5fec77455c6ada6c1f7ac4776d2",
    )
]

# Bash recipe for building across all platforms
script = raw"""
mv ${WORKSPACE}/srcdir/go ${prefix}/
"""

# We only build for host platform: x86_64-linux-musl
platforms = [
    host_platform,
]

# Dependencies that must be installed before this package can be built
dependencies = []

# The products that we will ensure are always built
products = [
    ExecutableProduct("go", :go, "go/bin"),
    ExecutableProduct("gofmt", :gofmt, "go/bin"),
]

# Build the tarballs, and possibly a `build.jl` as well.
ndARGS, deploy_target = find_deploy_arg(ARGS)
build_info = build_tarballs(ndARGS, name, version, sources, script, platforms, products, dependencies; skip_audit=true)

if deploy_target !== nothing
    upload_and_insert_shards(deploy_target, name, version, build_info)
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    ParameterizedBoundaryConditionFunction(func, parameters)

A wrapper for boundary condition functions implemented as functions with parameters.
The boundary condition `func` is called with the signature

    `func(i, j, grid, clock, state, parameters)`

where `i, j` are the indices along the boundary, `grid` and `clock` are `model.grid` and
`model.clock`, and `state` is a `NamedTuple` with fields `velocities`, `tracers`,
and `diffusivities`, which are each `NamedTuple`s of `OffsetArray`s that reference
the data associated with their corresponding fields.

Example
=======

@inline linear_drag(i, j, grid, clock, state, parameters) =
    @inbounds - parameters.μ * state.velocities.u[i, j, 1]

u_boundary_condition = ParameterizedBoundaryCondition(linear_drag, (μ=π,))
"""
struct ParameterizedBoundaryConditionFunction{F, P} <: Function
    func :: F
    parameters :: P
end

@inline (bc::ParameterizedBoundaryConditionFunction)(args...) = bc.func(args..., bc.parameters)

"""
    ParameterizedBoundaryCondition(bctype, func, parameters)

Returns a `BoundaryCondition` of `bctype` with a `ParameterizedBoundaryConditionFunction`
with function `func` and `parameters`.
"""
ParameterizedBoundaryCondition(bctype, func, parameters) =
    BoundaryCondition(bctype, ParameterizedBoundaryConditionFunction(func, parameters))


text:
# file for determining MeshEntity orientations
# not generalizable (but may not be necessary for DG)
function getEntityOrientations(mesh::PumiMesh)
# get the offset for node access of each mesh entity of each element
# to read/write from a Pumi field, accessing the data stored on an edge looks like:
# apf.getComponents(f, e, abs(offset - node) - 1, vals)
# where f is the field, e is a mesh entity, node is the node index (one based),
# and vals is an array of values to be written to the field
# offset = 0 corresponds to the original ordering
# offset = number of nodes on the entity + 1 corresponds to reversing the
#  ordering
# other values of offset should corresponds to rotations in 3d
# for 2d, the offsets for all nodes on an entity would typically be the same,
# corresponding to either flipping or not flipping the edge, but in 3d they
# might be different for a rotation
  # create array of arrays
  # hold an offset for every node of every entity of every element
  offsets = zeros(UInt8, mesh.numNodesPerElement, mesh.numEl)
  # hold a bit telling if the entity is in the proper orientation for this element
  # ie. it is "owned" by this element, so the nodes can be accessed in order

  # the elementNodeOffsets array is internal, therefore should be
  # stored in Pumi node order, not SBP order
  flags = Array{BitArray{2}}(3)
  # create the arrays
  for i=1:3
    flags[i] = trues(mesh.numTypePerElement[i], mesh.numEl)
  end

  # now do edges

  edge_flags = flags[2]
  edges_i = Array{Ptr{Void}}(12)
  edges_start = mesh.typeOffsetsPerElement[2]
  nnodes_per_edge = mesh.numNodesPerType[2]
  for i=1:mesh.numEl
    apf.getDownward(mesh.m_ptr, mesh.elements[i], 1, edges_i)

    for j=1:mesh.numTypePerElement[2]  # loop over edges on element i
      # get orientation of edge j on element i
      # get global number of edge j
      edgenum_global = apf.getNumberJ(mesh.edge_Nptr, edges_i[j], 0, 0) + 1

      orient, edge_idx = getEdgeOrientation(mesh, i, edgenum_global)

      # figure out offset value
      # write n = mesh.numNodesPerType[2] + 1 if orientation = -1, or n=0 if orientation=1
      if orient == 1
	offset_val = 0
      else
        offset_val = mesh.numNodesPerType[2] + 1
      end

      # starting index for nodes on this edge in Pumi ordering
      start_idx = edges_start + (edge_idx - 1)*nnodes_per_edge

      # take into account Pumi to SBP mapping
      for k=1:mesh.numNodesPerType[2]  # loop over nodes on the edge
	# figure out the SBP index of this node
	node_idx = start_idx + k - 1  # pumi node_idx
        offsets[node_idx, i] = offset_val
      end

      # save value to flag array
      edge_flags[j, i] = div(orient + 1, 2)

    end  # end loop over edges
  end # end loop over elements

  # now do faces
#  getFaceOffsets(mesh, offsets, flags)

  return offsets, flags
end

function getEdgeOrientation(mesh::PumiMesh, elnum::Integer, edgenum::Integer)
# figure out what the orientation of the specified edge is relative ot the
# element by looking at the order of the vertices that define the edge
# if the edge is in the same orientation as the element, return 1, otherwise -1
# because we are dealing with simplex elements, the third edge has to be
# handled specially

#  println("\nEntered getEdgeOrientation")
  # get all the vertices
  el_verts, tmp = apf.getDownward(mesh.m_ptr, mesh.elements[elnum], 0)
  edge_verts, tmp = apf.getDownward(mesh.m_ptr, mesh.edges[edgenum], 0)
  el_edges, tmp = apf.getDownward(mesh.m_ptr, mesh.elements[elnum], 1)

  # find out which edge of the face it is
  edge_idx = 0
  for i=1:3
    edgenum_ret = apf.getNumberJ(mesh.edge_Nptr, el_edges[i], 0, 0) + 1
    if edgenum_ret == edgenum
      edge_idx = i
      break
    end
  end

  # an edge is defined by its two vertices
  # therefore, if the vertices are in the same order on the edge as on
  # the element, then the edge is in the same orientation as the element,
  # otherwise opposite orientation.
  # For simpleces, the 3rd edge has to be special cased because it uses
  # vertices 3 and 1
  pos1 = 0  # position of edge_verts[1] in el_verts
  pos2 = 0  # position of edge_verts[2] in el_verts

  # find the positions of the edge vertices on the element
  for i=1:3  # loop over el_verts
    if el_verts[i] == edge_verts[1]
      pos1 = i
    elseif el_verts[i] == edge_verts[2]
      pos2 = i
    end
  end

  @assert pos1 != 0
  @assert pos2 != 0

  # use the vertex positions to determine oridentation
  if edge_idx <= 2
    if pos2 - pos1 > 0  # positively oriented
      return 1, edge_idx
    elseif pos2 - pos1 < 0  # negatively oriented
      return -1, edge_idx
    else
      throw(ErrorException("bad orientation determination"))
    end
  elseif edge_idx == 3  # ordering is reversed for 3rd edge
    if -(pos2 - pos1) > 0  # positively oriented
      return 1, edge_idx
    elseif -(pos2 - pos1) < 0  # negatively oriented
      return -1, edge_idx
    else
      throw(ErrorException("bad orientation determination"))
    end
  else
    throw(ErrorException("bad orientation determination"))
  end

end  # end function

struct EntityOrientation
  which::Cint
  flip::Bool
  rotate::Cint
end

function getEntityOrientations(mesh::PumiMesh3DG)
# for DG (without subtriangulation), these shouldn't ever be needed

  offsets = zeros(UInt8, mesh.numNodesPerElement, mesh.numEl)
  flags = Array{BitArray{2}}(4)  # on array for Verts, edges, faces, regiosn
  for i=1:4
    flags[i] = falses(mesh.numTypePerElement[i], mesh.numEl)
  end

  return offsets, flags
end

function getRelRotate(mesh::PumiMesh3, elementL::Integer, elementR::Integer, facenum::Integer)
# calculate the rotation of faceR relative to faceL, where faceL and faceR
# are shared between the two element
# this uses getAligmnet, which describes how to transform the face from
# its current orientation to the canonical orientation as part of its parent
# element
# the difference in the required rotations is the relative rotation
# elementL : global element index
# elementR : global element index
# faceL : local face number of shared face on elementL
# faceR : local face number of shared face on elementR
# face : global face number of the shared face

eL = mesh.elements[elementL]
eR = mesh.elements[elementR]
face = mesh.faces[facenum]

# get rotations to bring faces into canonical orientation
# for each element
# we can ignore flip because one of the faces will always be flipped
whichL, flipL, rotateL = apf.getAlignment(mesh.m_ptr, eL, face)
whichR, flipR, rotateR = apf.getAlignment(mesh.m_ptr, eR, face)

r1 = EntityOrientation(whichL, flipL, rotateL)
r2 = EntityOrientation(whichR, flipR, rotateR)

rel_rotate = calcRelRotation(mesh, r1, r2)
return rel_rotate

end  # end function

# calculate the relative rotation of 2 faces
function calcRelRotation(mesh::PumiMesh, r1::EntityOrientation,
                         r2::EntityOrientation, localfaces=true)

  if r1.flip == r2.flip && localfaces
    throw(ErrorException("Both faces cannot be flipped"))
  end

  # sum rotations because they each rotate ccw in their parent
  # element's orientation, so they rotate in opposite directions
  rel_rotate = r1.rotate + r2.rotate
  #println("rel_rotate before wrapping = ", rel_rotate)
  rel_rotate = wrapNumber(rel_rotate, 0, 2)
  #println("rel_rotate after wrapping = ", rel_rotate)

  # add 1 so output is in range [1,3]
  return rel_rotate + 1
end

function wrapNumber(num::Integer, lower::Integer, upper::Integer)
# make num perioid, where upper and lower are the max and min values allowable
# ie. num can only be in the range [1 3], if num == 0, then it gets mapped to 3
# similarly, 4 gets mapped to 1

range = upper - lower + 1

if num < lower
  diff = lower - num
  return upper - (diff % range) + 1

elseif  num > upper
  diff = num - upper
  return lower + (diff % range) - 1
else
  return num
end

end
"""
  Type to hold data about 2 elements that share a face, along with some
  arrays needed by getRelativeOrientation

  **Fields**

   * elnumL: number of the left element
   * eL: the MeshEntity * of the left element
   * elnumR: the number of the right element
   * elR: the MeshEntity * of the right element
   * faceL: the local number of the face from the perspective of elL
   * faceR: the local number of the face from the perspective of elR
   * vertsL: array of length 4 holding the MeshEntity* for the vertices of
             elL
   * vertsR: same as vertsL, but for elR
   * facevertsL: array holding the MeshEntity* vertices on faceL of elL
   * facevertsR: same as facevertsL, but for elR
   * part_nums: array to hold part numbers retrieved from Pumi (length 400?)
   * matched_entities: array to hold MeshEntity * retrieves from
     Pumi (length 400?)
   * match2: bool indicating if faceL and faceR are matched

  When the type is constructed, the first 6 fields and the match2 field
  should have the correct values, all the other fields can be empty arrays
  of the correct size.

"""
struct FaceData
  elnumL::Int
  elL::Ptr{Void}
  elnumR::Int
  elR::Ptr{Void}
  faceL::Int  # local face numbers
  faceR::Int
  vertsL::Array{Ptr{Void}, 1}
  vertsR::Array{Ptr{Void}, 1}
  facevertsL::Array{Ptr{Void}, 1}
  facevertsR::Array{Ptr{Void}, 1}
  part_nums::Array{Cint, 1}
  matched_entities::Array{Ptr{Void}, 1}
  match2::Bool  # whether or not the second face is a Matched entity
end

"""
  Maps the position where vertex 1 of faceL is found in faceR to the relative
  rotation of the face.
"""
global const pos_to_rotation = [1, 2, 3]

"""
  This function takes the data from element that share a face and determines
  the relative orientation of the shared face (ie. how many rotations cw
  faceR is from faceL when looking from the perspective of faceL).  The
  user supplied topology information is used to determine the ordering of the
  vertices on both faces, therefore the relative orientation is consistent with
  the users definition of the faces
"""
function getRelativeOrientation(fdata::FaceData, mesh::PumiMesh3DG)

  faceL = fdata.faceL
  faceR = fdata.faceR
  vertmap = mesh.topo.face_verts

  # get all vertices of the tet
  apf.getDownward(mesh.m_ptr, fdata.elL, 0, fdata.vertsL)
  apf.getDownward(mesh.m_ptr, fdata.elR, 0, fdata.vertsR)
  # extract the face vertices
  # uses the user suppied topology information to order the face verts
  for i=1:3
    fdata.facevertsL[i] = fdata.vertsL[vertmap[i, faceL]]
    fdata.facevertsR[i] = fdata.vertsR[vertmap[i, faceR]]
  end

  if fdata.match2
    getMatchedVertsR(fdata, mesh)
  end

  return calcRelativeOrientation(fdata.facevertsL, fdata.facevertsR)

end

function getMatchedVertsR(fdata::FaceData, mesh::PumiMesh3DG)
# do an in place replacement of facevertsR with their matches

  # temporary arrays
  part_nums = fdata.part_nums
  matched_entities = fdata.matched_entities

  # if any have matches, replace with remote pointer
  for i=1:3
    n = apf.countMatches(mesh.m_ptr, fdata.facevertsR[i])
    apf.getMatches(part_nums, matched_entities)
    # make sure it is possible to correlate the vertices uniquely
    assertUnique(sview(part_nums, 1:Int(n)), sview(matched_entities, 1:Int(n)),
                 mesh.myrank, fdata.facevertsL)
    foundflag = false
    for j=1:n  # find the match on the same part
      nfound = countin(matched_entities[j], fdata.facevertsL)
      if part_nums[j] == mesh.myrank && ( nfound == 1)
        fdata.facevertsR[i] = matched_entities[j]
        foundflag = true
        break
      end
    end

    if !foundflag
      # we get here if either there were no matches or no unique matches
      @assert false
    end
  end

  return nothing
end

"""
  This function verifies that it is possible to determine the relative
  orientation of two faces that are matched.

  Specifically,
  Given two vectors of the same length, part_nums and matched entities, this
  function assures that for all entities on the same part as myrank, the
  matched entities appears no more than once in other_entities

  Inputs:
    part_nums: array of integers containing the MPI ranks of the processes
               on which the corresponding entries in matched_entities reside
    matched_entities: the entities with the associated MPI ranks
    myrank: the MPI rank of the current process
    other_entities: entities to compare matched_entities against

"""
function assertUnique(part_nums::AbstractVector,
                      matched_entities::AbstractVector, myrank::Integer,
                      other_entities::AbstractVector)

  n = length(part_nums)
  numfound = 0
  for i=1:n
    if part_nums[i] == myrank
      numfound += countin(matched_entities[i], other_entities)
    end
  end

  # for each vertex, it matches cannot contain more than one vertex of the
  # other face
  @assert numfound == 1

  return nothing
end

function calcRelativeOrientation(facevertsL::AbstractArray, facevertsR::AbstractArray)

  # find vert1 of faceL in faceR
  idx = 0
  v1 = facevertsL[1]
  for i=1:3
    if facevertsR[i] == v1
      idx = i
      break
    end
  end

  return pos_to_rotation[idx]
end

function getRelativeOrientation(fdata::FaceData, mesh::PumiMesh2DG)

  return 1
end


text:
#=
This is a script that computes the total rate for non-thermal
radiation chemistry processes using experimental and theoretical
data.
=#

ϕ = 10.0E0 #cosmic-ray flux in cm^-2 s^-1
data = [0.88,0.59,6.22,2.86,1.66] # Ralf's experimental data in molecules/eV
S_e = 3.837E7 # Electronic stopping power in ASW in units of eV*cm^2/g
μ = 7.30606E-26 # 44 amu in grams (i.e. mass of acetaldehyle/vinylalcohol
G = μ*data
println("The experimental data in units of grams/eV is $G")

σ = G*S_e
println("The cross-section for this process is $σ")

σ_tot = sum(σ)
J = σ_tot*ϕ
println("********************************************")
println("The rate constant for the process is $J s^-1")
println("********************************************")


text:
exp(p::T) where {T<:AbstractPosit} = T(exp(Float64(p)))
exp2(p::T) where {T<:AbstractPosit} = T(exp2(Float64(p)))
exp10(p::T) where {T<:AbstractPosit} = T(exp10(Float64(p)))

for P in (:Posit8, :Posit16, :Posit32, :Posit8_1, :Posit16_1, :Posit24_1,
        :Posit8_2, :Posit16_2, :Posit24_2)
        @eval begin
                sin(p::Complex{$P}) = Complex{$P}(sin(Complex{Float64}(p)))
                cos(p::Complex{$P}) = Complex{$P}(cos(Complex{Float64}(p)))

                exp(p::Complex{$P}) = cos(im*p) - im*sin(im*p)
        end
end

expm1(p::T) where {T<:AbstractPosit} = exp(p) - one(T)

log(p::T) where {T<:AbstractPosit} = T(log(Float64(p)))
log2(p::T) where {T<:AbstractPosit} = T(log2(Float64(p)))
log10(p::T) where {T<:AbstractPosit} = T(log10(Float64(p)))

log1p(p::T) where {T<:AbstractPosit} = log(p) + one(T)

sin(p::T) where {T<:AbstractPosit} = T(sin(Float64(p)))
cos(p::T) where {T<:AbstractPosit} = T(cos(Float64(p)))
tan(p::T) where {T<:AbstractPosit} = T(tan(Float64(p)))

sincos(p::AbstractPosit) = sin(p),cos(p)


text:
# Autogenerated wrapper script for OpenSpecFun_jll for armv6l-linux-musleabihf-libgfortran5
export libopenspecfun

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenSpecFun")
JLLWrappers.@declare_library_product(libopenspecfun, "libopenspecfun.so.1")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libopenspecfun,
        "lib/libopenspecfun.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Failed attempt to reproduce the Taylor coefficients in appendix B of Rančić et al. (1996).

function taylor_coefficients(A)
    A₁_expected = 1.47713062600964
    for _ in 1:10
        @printf("A₁ ≈ %.16f + %.16fim, ΔA₁ = %.16f + %.16fim\n", real(A[1]), imag(A[1]), real(A₁_expected - A[1]), imag(A₁_expected - A[1]))
        ϕ = range(-π/4, π/4, length=101)
        z = @. r * exp(im * ϕ)
        z′ = @. 1 - z
        Z = @. z′ .^ 4
        W = [sum(A[k] * Z[n]^k for k in 1:length(A)) for n in 1:length(Z)]
        g = [W[n] * r^4 * exp(im * θ) for (n, θ) in enumerate(range(-π, π, length=length(W)))]
        g̃ = fft(g)
        A = [g̃[k] / r^(4k) for k in 1:length(g̃)]
    end
end


text:
module PProf

export pprof, @pprof

using Profile
using ProtoBuf
using OrderedCollections

using Profile: clear

"""
    PProf.clear()

Alias for `Profile.clear()`
"""
clear

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("PProf not installed properly, run Pkg.build(\"PProf\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

include(joinpath("..", "lib", "perftools.jl"))

import .perftools.profiles: ValueType, Sample, Function, Location, Line
const PProfile = perftools.profiles.Profile

const proc = Ref{Union{Base.Process,Nothing}}(nothing)

"""
    _enter!(dict::OrderedDict{T, Int64}, key::T) where T

Resolves from `key` to the index (zero-based) in the dict.
Useful for the Strings table

NOTE: We must use Int64 throughout this package (regardless of system word-size) b/c the
proto file specifies 64-bit integers.
"""
function _enter!(dict::OrderedDict{T,Int64}, key::T) where T
    if haskey(dict, key)
        return dict[key]
    else
        l = Int64(length(dict))
        dict[key] = l
        return l
    end
end

using Base.StackTraces: lookup, StackFrame

# TODO:
# - Mappings

"""
    pprof(data, period;
            web = true, webhost = "localhost", webport = 57599,
            out = "profile.pb.gz", from_c = true, drop_frames = "", keep_frames = "")

Fetches the collected `Profile` data, exports to the `pprof` format, and (optionally) opens
a `pprof` web-server for interactively viewing the results.

If `web=true`, the web-server is opened in the background. Re-running `pprof()` will refresh
the web-server to use the new output.

If you manually edit the output file, `PProf.refresh()` will refresh the server without
overwriting the output file. `PProf.kill()` will kill the server.

# Arguments:
- `data::Vector{UInt}`: The data provided by `Profile.fetch` [optional].
- `period::UInt64`: The sampling period in nanoseconds [optional].

# Keyword Arguments
- `web::Bool`: Whether to launch the `go tool pprof` interactive webserver for viewing results.
- `webhost::AbstractString`: If using `web`, which host to launch the webserver on.
- `webport::Integer`: If using `web`, which port to launch the webserver on.
- `out::String`: Filename for output.
- `from_c::Bool`: If `false`, exclude frames that come from from_c. Defaults to `true`.
- `drop_frames`: frames with function_name fully matching regexp string will be dropped from the samples,
                 along with their successors.
- `keep_frames`: frames with function_name fully matching regexp string will be kept, even if it matches drop_functions.
"""
function pprof(
    data::Union{Nothing,Vector{UInt}} = nothing,
    period::Union{Nothing,UInt64} = nothing;
    web::Bool = true,
    webhost::AbstractString = "localhost",
    webport::Integer = 57599,
    out::AbstractString = "profile.pb.gz",
    from_c::Bool = true,
    drop_frames::Union{Nothing,AbstractString} = nothing,
    keep_frames::Union{Nothing,AbstractString} = nothing,
)
    if data === nothing
        data = copy(Profile.fetch())
    end
    if period === nothing
        period = ccall(:jl_profile_delay_nsec, UInt64, ())
    end

    string_table = OrderedDict{AbstractString,Int64}()
    enter!(string) = _enter!(string_table, string)
    enter!(::Nothing) = _enter!(string_table, "nothing")
    ValueType!(_type, unit) = ValueType(_type = enter!(_type), unit = enter!(unit))

    # Setup:
    enter!("")  # NOTE: pprof requires first entry to be ""
    # Functions need a uid, we'll use the pointer for the method instance
    seen_funcs = Set{UInt64}()
    funcs = Dict{UInt64,Function}()

    seen_locs = Set{UInt64}()
    locs = Dict{UInt64,Location}()
    locs_from_c = Dict{UInt64,Bool}()

    sample_type = [
        ValueType!("events", "count"), # Mandatory
        ValueType!("stack_depth", "count"),
    ]

    prof = PProfile(
        sample = [],
        location = [],
        _function = [],
        mapping = [],
        string_table = [],
        sample_type = sample_type,
        default_sample_type = 1, # events
        period = period,
        period_type = ValueType!("cpu", "nanoseconds"),
    )

    if drop_frames !== nothing
        prof.drop_frames = enter!(drop_frames)
    end
    if keep_frames !== nothing
        prof.keep_frames = enter!(keep_frames)
    end

    # start decoding backtraces
    location_id = Vector{eltype(data)}()
    lastwaszero = true

    for ip in data
        # ip == 0x0 is the sentinel value for finishing a backtrace, therefore finising a sample
        if ip == 0
            # Avoid creating empty samples
            if lastwaszero
                @assert length(location_id) == 0
                continue
            end

            # End of sample
            value = [
                1,                   # events
                length(location_id), # stack_depth
            ]
            push!(prof.sample, Sample(; location_id = location_id, value = value))
            location_id = Vector{eltype(data)}()
            lastwaszero = true
            continue
        end
        lastwaszero = false

        # A backtrace consists of a set of IP (Instruction Pointers), each IP points
        # a single line of code and `litrace` has the necessary information to decode
        # that IP to a specific frame (or set of frames, if inlining occured).

        # if we have already seen this IP avoid decoding it again
        if ip in seen_locs
            # Only keep C frames if from_c=true
            if (from_c || !locs_from_c[ip])
                push!(location_id, ip)
            end
            continue
        end
        push!(seen_locs, ip)

        # Decode the IP into information about this stack frame (or frames given inlining)
        location = Location(; id = ip, address = ip, line = [])
        location_from_c = true
        for frame in lookup(ip)
            # ip 0 is reserved
            frame.pointer == 0 && continue
            # if any of the frames is not from_c the entire location is not from_c
            location_from_c &= frame.from_c

            push!(location.line, Line(function_id = frame.pointer, line = frame.line))
            # Known function
            frame.pointer in seen_funcs && continue
            push!(seen_funcs, frame.pointer)

            # Store the function in our functions dict
            funcProto = Function()
            funcProto.id = frame.pointer
            file = nothing
            if frame.linfo !== nothing && frame.linfo isa Core.MethodInstance
                linfo = frame.linfo::Core.MethodInstance
                meth = linfo.def
                file = string(meth.file)
                funcProto.name = enter!(string(meth.module, ".", meth.name))
                funcProto.start_line = convert(Int64, meth.line)
            else
                # frame.linfo either nothing or CodeInfo, either way fallback
                file = string(frame.file)
                funcProto.name = enter!(string(frame.func))
                funcProto.start_line = convert(Int64, frame.line) # TODO: Get start_line properly
            end
            file = Base.find_source_file(file)
            funcProto.filename = enter!(file)
            funcProto.system_name = funcProto.name
            # Only keep C functions if from_c=true
            if (from_c || !frame.from_c)
                funcs[frame.pointer] = funcProto
            end
        end
        locs_from_c[ip] = location_from_c
        # Only keep C frames if from_c=true
        if (from_c || !location_from_c)
            locs[ip] = location
            push!(location_id, ip)
        end
    end

    # Build Profile
    prof.string_table = collect(keys(string_table))
    # If from_c=false funcs and locs should NOT contain C functions
    prof._function = collect(values(funcs))
    prof.location = collect(values(locs))

    # Write to disk
    open(out, "w") do io
        writeproto(io, prof)
    end

    if web
        refresh(webhost = webhost, webport = webport, file = out)
    end

    out
end

"""
    refresh(; webhost = "localhost", webport = 57599, file = "profile.pb.gz")

Start or restart the go pprof webserver.

- `webhost::AbstractString`: Which host to launch the webserver on.
- `webport::Integer`: Which port to launch the webserver on.
- `file::String`: Profile file to open.
"""
function refresh(
    ;
    webhost::AbstractString = "localhost",
    webport::Integer = 57599,
    file::AbstractString = "profile.pb.gz",
)

    if proc[] === nothing
        # The first time, register an atexit hook to kill the web server.
        atexit(PProf.kill)
    else
        # On subsequent calls, restart the pprof web server.
        Base.kill(proc[])
    end

    proc[] = open(pipeline(`$go_pprof -http=$webhost:$webport $file`))
end

"""
    pprof_kill()

Kills the pprof server if running.
"""
function kill()
    if proc[] !== nothing
        Base.kill(proc[])
        proc[] = nothing
    end
end

"""
    @pprof ex

Profiles the expression using `@profile` and starts or restarts `pprof`.
"""
macro pprof(ex)
    esc(quote
        $Profile.@profile $ex
        $(@__MODULE__).pprof()
    end)
end

end # module


text:
using Sockets
include("strings.jl")

function total_value(hand)
    total = 0
    aces = 0
    for card in hand
        if card == 14
            aces += 1
            total += 11
        elseif card > 10
            total += 10
        else
            total += card
        end
    end

    for i in 1:aces
        if(total > 21)
            total -= 10
        end
    end
    return total
end
function show_cards(hand, all = true)
    s = "** "
    if !all && size(hand,1) > 0
        if hand[1] == 11 s *= "Jack"
        elseif hand[1] == 12 s *= "Queen"
        elseif hand[1] == 13 s *= "King"
        elseif hand[1] == 14 s *= "Ace"
        else
            s = string(s, hand[1])
        end
        return s * " **"
    end
    for card in hand
        if card == 11 s *= "Jack"
        elseif card == 12 s *= "Queen"
        elseif card == 13 s *= "King"
        elseif card == 14 s *= "Ace"
        else
            s = string(s, card)
        end
        s *= " "
    end
    return string(s, "** ", "with a total sum of ", total_value(hand), "")
end
function hit(hand, deck)
    append!(hand, pop!(deck))
end
function natural(p :: Player, bet, dealer_hand, player_hand)
    dealer_total = total_value(dealer_hand)
    player_total = total_value(player_hand)

    if size(player_hand,1) == 2 && size(dealer_hand,1) == 2
        if player_total == 21 && dealer_total == 21
            print_dict(p, "black_jack_0")
            return true
        elseif player_total == 21
            print_dict(p, "black_jack_1")
            p.balance += bet ÷ 2
            return true
        elseif dealer_total == 21
            print_dict(p, "black_jack_2")
            p.balance -= bet
            return true
        end
    end
    return false
end

function evaluate(p :: Player, bet, dealer_hand, player_hand)
    dealer_total = total_value(dealer_hand)
    player_total = total_value(player_hand)

    if dealer_total > player_total
        write(p.socket, "The dealer wins with $(show_cards(dealer_hand)) against your $(show_cards(player_hand)).\nBetter luck next time!\n")
        p.balance -= bet
    elseif dealer_total < player_total
        write(p.socket, "You win with $(show_cards(player_hand)) against the dealers $(show_cards(dealer_hand)).\nCongratulations!\n")
        p.balance += bet
    else
        write(p.socket, "A standoff! The dealers $(show_cards(dealer_hand)) against your $(show_cards(player_hand)).\nWell played!\n")
    end
end

function play_black_jack(p::Player)
    bet = 0
    while true
        print_dict(p, "black_jack_3")
        printBalance(p)
        s = readline(p.socket)

        bet = tryparse(Int64, s)
        if bet == nothing || bet < 0
            print_dict(p, "repeat")
            continue
        elseif bet > p.balance
            print_dict(p, "black_jack_4")
            continue
        elseif bet == 0
            print_dict(p, "black_jack_5")
            p.status = reception
            return
        else
            break
        end
    end

    deck = repeat(collect(Int8, 2:14),4)

    shuffle!(deck)
    dealer_hand = [pop!(deck), pop!(deck)]
    player_hand = [pop!(deck), pop!(deck)]

    write(p.socket, "The dealer shows: $(show_cards(dealer_hand, false)).\n")
    write(p.socket, "You hold: $(show_cards(player_hand)).\n")

    if natural(p, bet, dealer_hand, player_hand)
        return
    end

    choice = "s"
    print_dict(p, "black_jack_6")
    choice = readline(p.socket)
    while choice != "s"
        if choice == "h"
            hit(player_hand, deck)
            write(p.socket, "You now hold: $(show_cards(player_hand))\n")
            if total_value(player_hand) > 21
                print_dict(p, "black_jack_7")
                p.balance -= bet
                return
            end
        else
            print_dict(p, "repeat")
        end

        print_dict(p, "black_jack_6")
        choice = readline(p.socket)
    end
    print_dict(p, "spacer")
    print_dict(p, "black_jack_8")

    write(p.socket, "The dealer shows: $(show_cards(dealer_hand))\n")
    dealer_total = total_value(dealer_hand)
    while dealer_total < 17
        hit(dealer_hand, deck)
        write(p.socket, "The dealer shows: $(show_cards(dealer_hand))\n")
        dealer_total = total_value(dealer_hand)
        if dealer_total > 21
            print_dict(p, "black_jack_9")
            p.balance += bet
            return
        end
    end
    evaluate(p, bet, dealer_hand, player_hand)
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct Attribute{S, T} <: AbstractAttribute
    iid::String
    type::AttributeType{S}
    value::T
end

function Attribute(t::Proto.Thing)
    iid = bytes2hex(t.iid)
    isempty(iid) && throw(TypeDBClientException(CONCEPT_MISSING_IID))

    attribute_type = instantiate(t._type)
    value_type = first(typeof(attribute_type).parameters)

    if value_type == Proto.AttributeType_ValueType.BOOLEAN
        return Attribute(iid, attribute_type, t.value.boolean)
    elseif value_type == Proto.AttributeType_ValueType.LONG
        return Attribute(iid, attribute_type, t.value.long)
    elseif value_type == Proto.AttributeType_ValueType.DOUBLE
        return Attribute(iid, attribute_type, t.value.double)
    elseif value_type == Proto.AttributeType_ValueType.STRING
        return Attribute(iid, attribute_type, t.value.string)
    elseif value_type == Proto.AttributeType_ValueType.DATETIME
        return Attribute(iid, attribute_type, t.value.date_time)
    else
        throw(TypeDBClientException(CONCEPT_BAD_VALUE_TYPE))
    end
end


text:
function FieldEventsLogger()
    return FieldEventsLogger(())
end

function clear_logged_events(obj::FieldEventsLogger)
    return jcall(obj, "clearLoggedEvents", void, ())
end

function get_logged_events(obj::FieldEventsLogger)
    return jcall(obj, "getLoggedEvents", List, ())
end

function monitor_detector(obj::FieldEventsLogger, arg0::FieldEventDetector)
    return jcall(obj, "monitorDetector", FieldEventDetector, (FieldEventDetector,), arg0)
end


text:
"""
    vlsrelief(data::Array{<:Real,2}, target::Array{<:Integer,1}, subset_size::Integer=-1,
                   num_subsets::Integer=100; rba::Any=relieff)::Array{Float64,1}

Compute feature weights using VLSRelief algorithm. The num_partitions_to_select argument specifies how many partitions
to select in each iteration. The num_subsets argument specifies the number of subset iterations to perform. The subset_size
argument specifies the size of a partition. The rba argument specifies a (partially applied) wrapped
RBA algorithm that should accept just the data and target values.

---
# Reference:
- Margaret Eppstein and Paul Haake. Very large scale ReliefF for genome-
wide association analysis. In 2008 IEEE Symposium on Computational
Intelligence in Bioinformatics and Computational Biology, CIBCB ’08,
2008.
"""
function vlsrelief(data::Array{<:Real,2}, target::Array{<:Integer,1}, subset_size::Integer=-1,
                   num_subsets::Integer=100; rba::Any=relieff)::Array{Float64,1}

    # If subset_size argument has signal value of -1, compute default value from data size.
    if subset_size == -1
        subset_size = Int64(ceil(size(data, 2) / 10))
    elseif subset_size > size(data, 2)  # If specified subset size greater than number of features, throw error.
        throw(DomainError(subset_size, "Feature subset size cannot exceed the number of features in the dataset."))
    end

    # Initialize feature weights vector.
    weights = zeros(Float64, size(data, 2))

    # Go over feature subsets.
    @inbounds for i = 1:num_subsets

        # Sample features.
        sample_idxs = StatsBase.sample(1:size(data, 2), subset_size, replace=false)

        # Use RBA on subset to obtain local weights.
        rba_weights = rba(data[:, sample_idxs], target)

        # Update weights using computed local weights.
        weights[sample_idxs] = vec(maximum(hcat(weights[sample_idxs], rba_weights), dims=2))
    end

    # Return computed feature weights.
    return weights
end


text:
# Tests for entropy function

arr = rand(1000)
sqrt_1000 = sqrt(1000)

# Test with 1 values array
@test get_entropy(arr) ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with 1 array passed.")

# Test with 2 values arrays
@test get_entropy(arr, arr) ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with 2 arrays passed.")

# Test with 3 values arrays
@test get_entropy(arr, arr, arr) ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with 3 arrays passed.")

# Test with uniform_width
@test get_entropy(arr, mode = "uniform_width") ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with uniform width passed.")

# Test with uniform_count
@test get_entropy(arr, mode = "uniform_count") ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with uniform count passed.")

# Test with bayesian_blocks
# @test get_entropy(arr, mode = "bayesian_blocks") ≈ log2(10) atol = 0.05
# println("Entropy with Bayesian blocks passed.")

# Test with number_of_bins
@test get_entropy(arr, number_of_bins = 5) ≈ log2(5) atol = 0.05
println("Entropy with number of bins passed.")

# Test with get_number_of_bins
@test get_entropy(arr, get_number_of_bins = function(values) return 2 end) ≈ log2(2) atol = 0.05
println("Entropy with get number of bins passed.")

# Test with maximum_likelihood
@test get_entropy(arr, estimator = "maximum_likelihood") ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with maximum likelihood passed.")

# Test with shrinkage
@test get_entropy(arr, estimator = "shrinkage") ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with shrinkage passed.")

# Test with shrinkage and lambda
@test get_entropy(arr, estimator = "shrinkage", lambda = 0) ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with shrinkage and lambda passed.")

# Test with dirichlet
@test get_entropy(arr, estimator = "dirichlet") ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with Dirichlet passed.")

# Test with dirichlet and prior
@test get_entropy(arr, estimator = "dirichlet", prior = 1) ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with Dirichlet and prior passed.")

# Test with miller_madow
@test get_entropy(arr, estimator = "miller_madow") ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with Miller-Madow passed.")

# Test with base e
@test get_entropy(arr, base = e) ≈ log(sqrt_1000) atol = 0.05
println("Entropy with change of base passed.")

# Test with discretized
@test get_entropy(discretize_values(arr), discretized = true) ≈ log2(sqrt_1000) atol = 0.05
println("Entropy with discretized input passed.")


text:
# main.jl runs on small systems sizes, to test that code is fully run-able

function maintest()
    dl = [4,6,8]
    cl = zeros(size(dl))

    for i in 1:length(dl)
    # define a model
    dims = (dl[i],dl[i]) # dimension of the lattice
    lattice = Lattices.get_SquareLattice(dims) # make a square lattice
    filling = 1 # setting filling ≢ 1 means there are holes.

    model = Dict(
        "dims" => dims,
        "lattice" => lattice,
        "filling" => filling,
        "hamiltonian" => Lattices.get_Tightbinding_Wavefunctions,
        "fermi_energy" => 0,
        )

    policy = FFG.SwapNeighborsPolicy()

    mc_spec = Dict(
        "mc_warmup_steps" => Int(1e1*prod(dims)),
        "mc_steps" => Int(1e2*prod(dims)),
        "sample_interval" => Int(prod(dims)*4),
        )

    gutzwiller_chain = FFG.GutzwillerChain()

    FFG.init_Chain!(gutzwiller_chain,
        model=model,observable=FFG.NeighborSzSz,
        policy=policy,mc_spec=mc_spec)

    FFG.get_Mc_Spec(gutzwiller_chain)

    FFG.runMC!(gutzwiller_chain)

    cl[i] = FFG.get_Data(gutzwiller_chain)["SzSz"]*3;
    end

    return cl
end


text:
const N = 624
const M = 397
const UPPER_MASK = 0x80000000
const LOWER_MASK = 0x7fffffff

mutable struct MT19937
    mt::Vector{UInt32}
    mti::Int
    function MT19937(x::Vector{UInt32}, i::Int)
        @assert length(x) == N
        new(x, i)
    end
end

@inline mt_magic(y) = ((y % Int32) << 31 >> 31) & 0x9908b0df
"Get a random `UInt32` number from a `MT19937` object."
@inline function mt_get(r::MT19937)
    mt = r.mt
    if r.mti > N
        @inbounds for i in 1:N-M
            y = (mt[i] & UPPER_MASK) | (mt[i+1] & LOWER_MASK)
            mt[i] = mt[i + M] ⊻ (y >> 1) ⊻ mt_magic(y)
        end
        @inbounds for i in N-M+1:N-1
            y = (mt[i] & UPPER_MASK) | (mt[i+1] & LOWER_MASK)
            mt[i] = mt[i + M - N] ⊻ (y >> 1) ⊻ mt_magic(y)
        end
        @inbounds begin
            y = (mt[N] & UPPER_MASK) | (mt[1] & LOWER_MASK)
            mt[N] = mt[M] ⊻ (y >> 1) ⊻ mt_magic(y)
        end
        r.mti = 1
    end
    k = mt[r.mti]
    k ⊻= (k >> 11)
    k ⊻= (k << 7) & 0x9d2c5680
    k ⊻= (k << 15) & 0xefc60000
    k ⊻= (k >> 18)

    r.mti += 1
    k
end

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
using Conda
using BinDeps
using Base.Test
using Compat

# Save previous state of the Conda installation
CURL_ALREADY_INSTALLED = false
PREVIOUS_CHANNELS = []
PREVIOUS_PACKAGES = []

"Save any state in the Conda installation"
function setup()
    global CURL_ALREADY_INSTALLED = "curl" in Conda._installed_packages()
    global PREVIOUS_PACKAGES = Conda._installed_packages_dict()
    global PREVIOUS_CHANNELS = Conda.channels()
    for channel in PREVIOUS_CHANNELS
        Conda.rm_channel(channel)
    end
    Conda.add_channel("defaults")
end

"Restore the state of the Conda installation"
function teardown()
    if CURL_ALREADY_INSTALLED
        Conda.add("curl")
    end

    for channel in PREVIOUS_CHANNELS
        Conda.add_channel(channel)
    end

    installed = Conda._installed_packages()

    for (name, (_, package)) in PREVIOUS_PACKAGES
        if !(name in installed)
            Conda.add(package)
        end
    end
end

function main()
    setup()

    @test Conda.exists("curl")
    Conda.add("curl")

    if is_unix()
        curl_path = joinpath(Conda.PREFIX, "bin", "curl-config")
    end
    if is_windows()
        manager = Conda.Manager(["curl"])
        curl_libpath = BinDeps.libdir(manager, "")
        curl_path = joinpath(curl_libpath, "curl.exe")
    end

    @test isfile(curl_path)

    @test isfile(joinpath(Conda.BINDIR, basename(curl_path)))

    Conda.rm("curl")
    if is_unix()
        @test !isfile(curl_path)
    end

    @test isfile(joinpath(Conda.SCRIPTDIR, "conda" * (is_windows() ? ".exe": "")))

    @test isfile(joinpath(Conda.PYTHONDIR, "python" * (is_windows() ? ".exe": "")))

    channels = Conda.channels()
    @test channels == ["defaults"]

    Conda.add_channel("foo")
    @test Conda.channels() == ["foo", "defaults"]
    # Testing that calling the function twice do not fail
    Conda.add_channel("foo")

    Conda.rm_channel("foo")
    channels = Conda.channels()
    @test channels == ["defaults"]

    teardown()
end

main()


text:
"""
# pairsplot

A simple version of StatsPlots' `cornerplot`.

$(SIGNATURES)

### Required arguments
```julia
* `df::DataFrame`                      : DataFrame containing the variables (as columns)
* `vars::Vector{Symbol}`               : Vector of variables to include in plot
```

"""
function pairsplot(df::DataFrame, vars::Vector{Symbol})

  l = length(vars)
  pmat = Matrix{Plots.Plot{Plots.GRBackend}}(undef, l, l)

  for i in 1:l
    for j in 1:l
      if i == j
        pmat[i, j] = plot(xlims=(0, 2), ylims=(0, 2), leg=false, framestyle=:none)
        plot!(pmat[i, j], annotations=(1, 1, Plots.text(String(vars[i]), :mid)),
          leg=false)
      else
        pmat[i,j] = scatter(df[:, vars[i]], df[:, vars[j]],
          markersize=5, color=:darkblue, leg=false)
      end
    end
  end
  plot(pmat..., layout=(l, l))

end

export
  pairsplot


text:
function cubeside_const_lift(_, id, mousehover)
    h_id, h_index = value(mousehover)
    if h_id == id && h_index >= 1 && h_index <= 6
        side =  CubeSides(mousehover.value[2]-1)
        side == TOP     && return Vec3f0(0,0,1)
        side == BOTTOM  && return Vec3f0(0,0,-1)
        side == FRONT   && return Vec3f0(0,1,0)
        side == BACK    && return Vec3f0(0,-1,0)
        side == LEFT    && return Vec3f0(1,0,0)
        side == RIGHT   && return Vec3f0(-1,0,0)
    end
    Vec3f0(1)
end
import Base: +
Base.clamp(x::T) where {T} = clamp(x, T(0),T(1))
+(a::RGBA{T}, b::Number) where {T} = RGBA{T}(clamp(comp1(a)+b), clamp(comp2(a)+b), clamp(comp3(a)+b), clamp(alpha(a)+b))

function cubeside_color(id, mousehover, startcolors, colortex)
    index = mousehover[2]
    update!(colortex, startcolors)
    if mousehover[1] == id
        (index >= 1 && index <= 6) && (colortex[index] = [startcolors[index] + 0.5])
    end
    nothing
end
function colored_cube()
    xdir    = Vec3f0(1f0,0f0,0f0)
    ydir    = Vec3f0(0f0,1f0,0f0)
    zdir    = Vec3f0(0f0,0f0,1f0)
    origin  = Vec3f0(-0.5)
    const quads = [
        (Quad(origin + zdir,   xdir, ydir), RGBA(1.0f0,0f0,0f0,1f0)), # Top
        (Quad(origin,          ydir, xdir), RGBA(0.5f0,0f0,0f0,1f0)), # Bottom
        (Quad(origin + ydir,   zdir, xdir), RGBA(0f0,1.0f0,0f0,1f0)), #Front
        (Quad(origin,          xdir, zdir), RGBA(0f0,0.5f0,0f0,1f0)), # Back
        (Quad(origin + xdir,   ydir, zdir), RGBA(0f0,0f0,1.0f0,1f0)), # Right
        (Quad(origin,          zdir, ydir), RGBA(0f0,0f0,0.5f0,1f0)), # Left
    ]
    cube_steering = merge(map(GLNormalMesh, quads))
end

Base.middle(r::SimpleRectangle{T}) where {T} = Point{2, T}(r.x+(r.w/T(2)), r.y+(r.mousehover/T(2)))

function get_rotation(m)
    idx = Vec3(1, 2, 3)
    xs = norm(m[1, idx])
    ys = norm(m[2, idx])
    zs = norm(m[3, idx])
    Mat4f0(
        m[1,1]/xs, m[1,2]/xs, m[1,3]/xs, 0,
        m[2,1]/ys, m[2,2]/ys, m[2,3]/ys, 0,
        m[3,1]/zs, m[3,2]/zs, m[3,3]/zs, 0,
        0     , 0     , 0     , 1,
    )
end

"""
Creates a camera which is steered by a cube for `window`.
"""
function cubecamera(
        window, projectiontype=Signal(GLAbstraction.PERSPECTIVE);
        cube_area   = Signal(SimpleRectangle(0,0,150,150)),
        eyeposition = Vec3f0(2),
        lookatv     = Vec3f0(0),
        trans       = Signal(Vec3f0(0)),
        theta       = Signal(Vec3f0(0))
    )
    T = Float32
    @materialize mouse_buttons_pressed, mouseposition, buttons_pressed = window.inputs
    dd = doubleclick(mouse_buttons_pressed, 0.3)
    mousehover = GLWindow.mouse2id(window)
    id = Signal(3)
    should_reset = filter(false, dd) do _
        x = value(mousehover).id == value(id)
    end
    cube_steering = colored_cube()
    new_eyeposition = preserve(const_lift(cubeside_const_lift,
        should_reset, id,
        Signal(mousehover)
    ))

    left_ctrl = Set([GLFW.KEY_LEFT_CONTROL])
    use_cam = map(buttons_pressed, window.inputs[:mouseinside]) do b, mi
        b == left_ctrl && mi
    end
    theta, trans = default_camera_control(
        window.inputs, Signal(0.02f0), Signal(1f0),
        use_cam
    )
    fov, upvector = map(Signal, (43f0, Vec3f0(0,0,1)))
    lookatvec, eyeposition =  Signal(lookatv), Signal(eyeposition)
    far = map(eyeposition, lookatvec) do a,b
        max(norm(b-a) * 5f0, 100f0)
    end
    near = map(eyeposition, lookatvec) do a,b
        norm(b-a) * 0.007f0
    end
    main_cam = PerspectiveCamera(
        theta ,trans, lookatvec, eyeposition, upvector,
        window.area, fov, near, far,
        projectiontype
    )
    preserve(map(new_eyeposition) do neweye_dir
        eypos, lookv = value(eyeposition), value(lookatvec)
        oldlength = norm(eypos - lookv)
        push!(eyeposition, lookv+neweye_dir*oldlength)
        for i=1:3
            unitvec = unit(Vec3f0, i)
            if dot(neweye_dir, unitvec) == 0f0 # find an orthogonal vector
                push!(upvector, unitvec)
            end
        end
    end)
    window.cameras[:perspective] = main_cam

    viewmatrix = map(main_cam.view) do m
        inv(get_rotation(m))
    end

    robj = visualize(cube_steering, model=viewmatrix)

    start_colors = cube_steering.attributes
    color_tex    = robj.children[][:attributes]
    preserve(const_lift(cubeside_color, id, mousehover, Signal(start_colors), Signal(color_tex)))
    preserve(id)
    push!(id, robj.children[].id)
    robj
end


text:
using JuDO
using Test
using DotEnv

import Dates: DateTime
import JuDO: getdata!, postdata!, putdata!, deletedata!, ENDPOINT, MAXOBJECTS
import JSON: parse

DotEnv.config()

DATADIR = joinpath(dirname(abspath(@__FILE__)), "data")

struct TestClient <: AbstractClient
    token::String
end

function getdata!(client::TestClient, uri::String)
    uri = replace(uri, ENDPOINT => "")
    # strip off query string
    idx = findfirst(isequal('?'), uri)
    if idx != nothing
        uri = uri[1:idx-1]
    end

    path = joinpath(DATADIR, uri, "get.json")
    data = open(path, "r") do io
        parse(read(io, String))
    end

    # get the only key which is neither meta or links
    pkey = pop!(setdiff(keys(data), ["meta", "links"]))
    data[pkey]
end

function postdata!(client::TestClient, uri::String, body::Dict{String})
    if endswith(uri, "actions")
        path = joinpath(DATADIR, replace(uri, ENDPOINT => ""),
                        "$(body["type"]).json")
    else
        path = joinpath(DATADIR, replace(uri, ENDPOINT => ""), "post.json")
    end
    data = open(path, "r") do io
        parse(read(io, String))
    end

    # data should have only one key
    pop!(data).second
end

function putdata!(client::TestClient, uri::String, body::Dict{String})
    path = joinpath(DATADIR, replace(uri, ENDPOINT => ""), "put.json")
    data = open(path, "r") do io
        parse(read(io, String))
    end

    # data should have only one key
    pop!(data).second
end

function deletedata!(client::TestClient, uri::String)
    true
end

client = TestClient("fake token")

# account tests
include("testaccount.jl")
# action tests
include("testactions.jl")
# block storage tests
include("testvolumes.jl")
# certificate tests
include("testcertificates.jl")
# domain tests
include("testdomains.jl")
# domain record tests
include("testrecords.jl")
# droplet tests
include("testdroplets.jl")

# real client test
if haskey(ENV, "DOTOKEN")
    client = Client(ENV["DOTOKEN"])
end


text:
using AirfoilDatabase
using Test
@testset "query" begin
    sd7003 = query_airfoil("7003")
    @test length(sd7003) == 1
    @test sd7003[1].name == "SD7003"
    sc2 = query_airfoil("NASA SC2")
    @test length(sc2) > 1
end

@testset "NACA" begin
    n0012 = query_airfoil("N0012")
    @test length(n0012) == 1
    x = n0012[1].x
    ile = argmin(x)
    xc = x[ile:-1:1]
    af = NACA("0012", xc)
    @test af.x[ile] == af.y[ile] == 0
    @test all(isapprox.(af.y, n0012[1].y, atol=1e-5))
    af = NACA("4412")
    @test length(af.x) == 65
    @test (af.x[end] + af.x[1])/2 == 1
    @test af.x[33] == af.y[33] == 0
    @test af.y[end] < af.y[1]
    af = NACA("4412", sharpTE=true)
    @test af.x[end] == af.x[1] == 1
    @test isapprox(af.y[end], af.y[1], atol=1e-6)
end


text:
using FrankWolfe
using LinearAlgebra
using LaTeXStrings
using Plots

# # FrankWolfe for scaled, shifted ``\ell^1`` and ``\ell^{\infty}`` norm balls

# In this example, we run the vanilla FrankWolfe algorithm on a scaled and shifted ``\ell^1`` and ``\ell^{\infty}`` norm ball, using the `ScaledBoundL1NormBall`
# and `ScaledBoundLInfNormBall` LMOs. We shift both onto the point ``(1,0)`` and then scale them by a factor of ``2`` along the x-axis. We project the point ``(2,1)`` onto the polytopes.

n = 2

k = 1000

xp = [2.0,1.0]

f(x) = norm(x-xp)^2

function grad!(storage,x)
    @. storage = 2 * (x - xp)
    return nothing
end

lower = [-1.0,-1.0]
upper = [3.0,1.0]

l1 = FrankWolfe.ScaledBoundL1NormBall(lower, upper)

linf = FrankWolfe.ScaledBoundLInfNormBall(lower, upper)

x1 = FrankWolfe.compute_extreme_point(l1, zeros(n))
gradient = collect(x1)

x_l1, v_1, primal_1, dual_gap_1, trajectory_1 = FrankWolfe.frank_wolfe(
    f,
    grad!,
    l1,
    collect(copy(x1)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=50,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

println("\nFinal solution: ", x_l1)

x2 = FrankWolfe.compute_extreme_point(linf, zeros(n))
gradient = collect(x2)

x_linf, v_2, primal_2, dual_gap_2, trajectory_2 = FrankWolfe.frank_wolfe(
    f,
    grad!,
    linf,
    collect(copy(x2)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=50,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

println("\nFinal solution: ", x_linf)

# We plot the polytopes alongside the solutions from above:

xcoord1 = [1,3,1,-1,1]
ycoord1 = [-1,0,1,0,-1]

xcoord2 = [3,3,-1,-1,3]
ycoord2 = [-1,1,1,-1,-1]

plot(xcoord1, ycoord1, title = "Visualization of scaled shifted norm balls", lw = 2, label = L"\ell^1 \textrm{ norm}")
plot!(xcoord2, ycoord2, lw = 2, label = L"\ell^{\infty} \textrm{ norm}")
plot!([x_l1[1]], [x_l1[2]], seriestype = :scatter, lw = 5, color = "blue", label = L"\ell^1 \textrm{ solution}")
plot!([x_linf[1]], [x_linf[2]], seriestype = :scatter, lw = 5, color = "orange", label = L"\ell^{\infty} \textrm{ solution}", legend = :bottomleft)


text:
abstract type LearningUpdate end
abstract type Learner end

export QLearner, LinearQLearner, SRLearner, LSTDLearner, GPI, predict
include("learners/value.jl")
include("learners/SR.jl")
include("learners/GPI.jl")
include("learners/LSTD.jl")

include("LearnedPolicy.jl")

update!(learner::Learner, args...) =
    update!(update(learner), learner, args...)

zero_eligibility_traces!(l::Learner) = zero_eligibility_traces!(l.update)

export TB, TBAuto, ESARSA, SR, update!, SARSA
abstract type AbstractTraceUpdate end
export AccumulatingTraces, ReplacingTraces
include("updates/update_utils.jl")
include("updates/TB.jl")
include("updates/SARSA.jl")
# include("updates/TB_Auto.jl")
include("updates/ESARSA.jl")
# include("updates/SR.jl")
include("updates/EmphESARSA.jl")
include("updates/ETB.jl")
include("updates/StateETB.jl")
include("updates/PriorESARSA.jl")
include("updates/PriorTB.jl")
include("updates/InterestTB.jl")

include("updates/SR_TB_Replay.jl")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# ------------------------------------------------------------------- #
# Allocations
# ------------------------------------------------------------------- #
"""
Computes `b2` given `a1` and `zhat` for model where `σ1 = σ2`

Utilitzes the planner's FOC wrt a1, a2, b1, and b2

The condition `σ1 = σ2` is not checked or enforced in this function. Thus,
unless you are certain that the σ's are the same, you should use `eval_b2(m, a1,
zh)` instead
"""
function analytical_b2(m::AbstractModel, a1, zh)
    ω1, ω2, σ1 = m.agent1.ω, m.agent2.ω, m.agent1.σ
    foo = (((1-ω1)*(1-ω2))/(ω1*ω2))^(1/(σ1-1))*(a1./(zh-a1))
    B = 1 ./ zh
    b2 = B ./ (1 + foo)
end

"Given `(a1, lzh)`, solve for b2"
function eval_b2(m::AbstractModel, a1, lzh)
    if is_same_sigma(m)
        analytical_b2(m, a1, exp(lzh))
    else
        # Could allow for possibility of different σ
        error("Need the values of sigma to be the same")
    end
end

"Given (lλ☆, lzh), solve for allocations `a1, b1, c1, a2, b2, c2`"
function get_allocation(m::AbstractModel, lλ☆, lzh)

    # Get a1 and b2 from interpolation
    a1 = m.a1_itp[lλ☆, lzh]
    b2 = m.b2_itp[lλ☆, lzh]

    # Now we can get everything else
    zh = exp(lzh)
    a2 = zh - a1
    b1 = 1./zh - b2
    c1 = _agg_c(m.agent1, a1, b1)
    c2 = _agg_c(m.agent2, b2, a2)

    return a1, b1, c1, a2, b2, c2
end

"Given `curr_state=(lλ☆, lzh)`, solve for allocations `a1, b1, c1, a2, b2, c2`"
get_allocations(m::AbstractModel, curr_state::Array{Float64, 1}) =
    get_allocations(m, curr_state[1], curr_state[2])

# ------------------------------------------------------------------- #
# Bounds
# ------------------------------------------------------------------- #
"""
Compute bounds on grids for exogenous variables, given their process

It is hard to pin down the exact bounds we care about for the exogenous
processes, so, instead of guessing, we will simulate for a large number
of periods and then add a little extra space on either end to get a
good idea of where the bounds should live for these processes.
"""
function get_exog_bounds(exog::AbstractExog)
    # Get number of exogenous states
    n_exog_states = n_exog(exog)

    # Use same seed and simulate for a long time
    srand(42)
    exog_sim = simulate(exog, init(exog), randn(n_exog_states, 2_500_000))

    # Get the exogenous state (es) bounds
    es_bounds = map(i->extrema(exog_sim[i, :]), collect(1:n_exog_states))

    # Expand bounds
    exp_es_bounds = Array(Tuple{Float64, Float64}, n_exog_states)
    for i=1:n_exog_states
        x_lb, x_ub = es_bounds[i][1], es_bounds[i][2]
        exp_bds = 0.025*abs(x_ub - x_lb)
        exp_es_bounds[i] = x_lb - exp_bds, x_ub + exp_bds
    end

    return exp_es_bounds
end

"""
Compute reasonable bounds to begin bisection on 'Lars k'

Note that

    log(Lars k) = (ρ2-α2) log(μ2) - (ρ1-α1) log(μ1)

Also notice that if

* (ρi - αi) > 0.
* μ2 is increasing in lλ☆ (μ1 is decreasing in lλ☆)
* μ2 is decreasing in lzh (μ1 is increasing in lzh) <- Don't think we use

We will have that  `Lars k` is increasing in lλ☆, which implies

* upper bound to lars k is when lλ☆ is big (and lzh is small)
* lower bound to lars k is when lλ☆ is small (and lzh is big)

## Parameters

- `i::Int` an integer specifying which state in `m.grid` we are at for time t
  variables
- `vf_itps::VFITP`: A scheme for interpolating over the time t+1 value functions

## Returns
- `lb, ub`: The lower and upper bounds for Lars' k at state `i`.
"""
function log_lars_k_bounds(m::AbstractModel, i::Int, vf_itps::VFITP)

    # unpack
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)
    Nϵ, Nshocks = size(m.ϵ)

    # current states
    lλ☆_i = m.grid[i, 1]

    # t+1 states
    lgp_i = get_exog_tomorrow(m, i)[end]
    gp_i = exp(lgp_i)

    # Build a function that gets the pieces we care about
    function get_k_K(lλ☆p)
        ljp, lup = Array(Float64, Nϵ), Array(Float64, Nϵ)

        for j=1:Nϵ
            # `curr_exog_statep` will be ordered `[lzh, others..., lg]`,
            # where others is empty for cosntant volatility model and equal to
            # `v` for the stochastic volatility model.
            curr_exog_statep = get_exog_tomorrow(m, i, j)
            ljp[j] = vf_itps[1][lλ☆p[j], curr_exog_statep[1:end-1]...]
            lup[j] = vf_itps[2][lλ☆p[j], curr_exog_statep[1:end-1]...]
        end
        # evaluate certainty equivalents
        lμ1 = eval_lμ1(m, gp_i, exp(ljp))
        lμ2 = eval_lμ2(m, gp_i, exp(lup))

        # now get the lars k implied from both the FOC wrt U' _and_ the
        # definition of lars' k from in the docstring
        k_upres = up_residual(m, ljp, lup, lλ☆p, lgp_i, lλ☆_i, 0.0)
        k_defn = (ρ2-α2)*lμ2 - (ρ1-α1)*lμ1

        return [k_upres; k_defn]
    end

    # TODO: Chase will you comment about what this stuff does? Maybe add it
    #       to docstring (please check the docstring)
    lλ_step = max(m.m_ls[1][end] / 12, .1)
    lλ☆_big, lλ☆_small = m.grid[i, 1] + lλ_step, m.grid[i, 1] - lλ_step

    div_val = 2.5
    bstep = round(Int, Nϵ/div_val)

    lb = sort(get_k_K(fill(lλ☆_small, Nϵ)))[end-bstep]
    ub = sort(get_k_K(fill(lλ☆_big, Nϵ)))[bstep]

    return lb, ub
end

# ------------------------------------------------------------------- #
# Residuals
# ------------------------------------------------------------------- #
"""
Evaluate the residual for Lars' k:

    lk = (ρ2-α2) log(μ2) - (ρ1-α1) log(μ1)
"""
function lk_residual(m::AbstractModel, lμ1, lμ2, lk)

    # Unpack
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)

    return ((ρ2-α2)*lμ2 - (ρ1-α1)*lμ1) - lk
end

"""
Compute the residual for the FOC wrt `U_{t+1}` at a particular t+1 state:

Note that here we take `log Lars k` as an argument, but when the argument is
set to 0.0 then the output of this function is the implied `log Lars k`.
"""
function up_residual(m::AbstractModel, ljp, lup, lλ☆p, lgp, lλ☆, lk)

    # Unpack
    ρ1, α1, β1, ω1, σ1 = _unpack(m.agent1)
    ρ2, α2, β2, ω2, σ2 = _unpack(m.agent2)

    log(β1/β2) + (α1-α2).*lgp + (α1-ρ1).*ljp - (α2-ρ2).*lup + lλ☆p - lλ☆ - lk
end


text:
"""
Nlayer_cylinder_DCS

Provides default parameters to use to compute the autocorrelation function g2 in a N-layered cylinder.
The arguments defined with a Vector should be the same length and same type.

# Keyword Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa::Vector{T}`: absorption coefficient (cm⁻¹)
- `μsp::Vector{T}`: reduced scattering coefficient (cm⁻¹)
- `n_med::Vector{T}`: medium's index of refraction
- `n_ext`: external medium's index of refraction (air or detector)
- `l`: the thicknesses of each layer (cm)
- `z`: the z-depth orthogonal from the boundary (cm)
- `a`: the radius of the cylinder (cm)
- `β`: constant in Siegert relation dependent on collection optics
- `BFi::Vector{T}`: Blood flow index ~αDb (cm²/s)
- `λ`: wavelength (nm)
- `bessels::AbstractVector{T}`: roots of J0

# Examples
```
julia> data = Nlayer_cylinder_DCS() # return default parameters
julia> data = Nlayer_cylinder_DCS(ρ = 1.5) # return ρ = 1.5 with the rest of the parameters given by defaults

# we can now define our correlation times τ
julia> τ = 10 .^(range(-10,stop=0,length=250))
julia> g2_DA_Nlay_cylinder_CW(τ, data)
```
"""
@with_kw struct Nlayer_cylinder_DCS{T <: Real}
    μsp::Vector{T} = [10.0, 10.0]                       # reduced scattering coefficient (1/cm)
    μa::Vector{T} = [0.1, 0.1]                          # absorption coefficient (1/cm)
    n_ext::T = 1.0                                      # surrounding index of refraction
    n_med::Vector{T} = [1.0, 1.0]                       # layers index of refraction

    l::Vector{T} = [0.5, 10.0]                          # length of cylinder layers (cm)
    ρ::Union{T, AbstractVector{T}} = 1.0                # source-detector separation (cm)
    a::T = 20.0                                         # radius of cylinder (cm)
    z::T = 0.0                                          # detector depth (cm)

    β::T = 1.0                                          # constant in Siegert relation dependent on collection optics
    BFi::Vector{T} = [2.0e-8, 2.0e-8]                   # Blood flow index ~αDb (cm²/s)
    λ::T = 750.0                                        # wavelength (nm)

    bessels::AbstractVector{T} = besselroots[1:1000]    # roots of J0

    N_quad::Int = 50                                    # number of nodes in gauss-legendre integration
    N_laplace::Int = 12                                 # number of laplace evaluations in time-domain inversion
end

"""
    g1_DA_Nlay_cylinder_CW(τ::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:500])

Compute the electric field autocorrelation function function g1 in a layered cylinder (CW).

# Arguments
- `τ`: time lags (s)
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `BFi`: Blood flow index ~αDb (cm²/s)
- `β`: constant in Siegert relation dependent on collection optics
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: length of cylinder layers (cm)
- `a`: cylinder radius (cm)
- `z`: the z-depth orthogonal from the boundary (cm)
- `λ`: wavelength (nm)
- `bessels`: roots of J0 (bessel function of first kind order zero)
"""
function g1_DA_Nlay_cylinder_CW(τ::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:500])
    g1 = similar(τ)
    μa_dynamic = similar(μa)

    k0 = 2 * π * n_med / (λ * 1e-7) # this should be in cm
    tmp = @. 2 * μsp * BFi * k0^2

    G0 = fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, bessels)

    Threads.@threads for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        G1 = fluence_DA_Nlay_cylinder_CW(ρ, μa_dynamic, μsp, n_ext, n_med, l, a, z, bessels)
        g1[ind] = G1 / G0
    end

    return g1
end

"""
    g1_DA_Nlay_cylinder_TD(τ::AbstractVector, t::AbstractFloat, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000], N_laplace = 8)
    g1_DA_Nlay_cylinder_TD(τ::AbstractVector, t::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000], N_laplace = 8, N_quad = 100)

Compute the time-domain electric field autocorrelation function function g1 in a layered cylinder.

If t is a AbstractVector then g1(tau) will be calculated by ∫P(t)g1(tau, t)dt with integration bounds given by t[1] and t[end].
This is needed for time-gated applications.

# Arguments
- `τ`: time lags (s)
- `ρ`: the source detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword arguments
- `BFi`: Blood flow index ~αDb (cm²/s)
- `β`: constant in Siegert relation dependent on collection optics
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: length of cylinder layers (cm)
- `a`: cylinder radius (cm)
- `z`: the z-depth orthogonal from the boundary (cm)
- `λ`: wavelength (nm)
- `bessels`: roots of J0 (bessel function of first kind order zero)
- `N_quad`: number of nodes in gauss-legendre integration
- `N_laplace`: number of laplace evaluations in inversion to time-domain
"""
function g1_DA_Nlay_cylinder_TD(τ::AbstractVector, t::AbstractFloat, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000], N_laplace = 8)
    g1 = similar(τ)
    μa_dynamic = similar(μa)

    k0 = 2 * π * n_med / (λ * 1e-7) # this should be in cm
    tmp = @. 2 * μsp * BFi * k0^2

    G0 = fluence_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, bessels, N = N_laplace)

    for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        G1 = fluence_DA_Nlay_cylinder_TD(t, ρ, μa_dynamic, μsp, n_ext, n_med, l, a, z, bessels, N = N_laplace)
        g1[ind] = G1 / G0
    end

    return g1
end

function g1_DA_Nlay_cylinder_TD(τ::AbstractVector, t::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:1000], N_laplace = 8, N_quad = 100)
    g1 = similar(τ)
    μa_dynamic = similar(μa)

    k0 = 2 * π * n_med / (λ * 1e-7) # this should be in cm
    tmp = @. 2 * μsp * BFi * k0^2
    x, w = gausslegendre(N_quad)

    tpsf_norm = integrate_compute_y_first(t -> fluence_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, bessels, N = N_laplace), x, w, t[1], t[2])

    for ind in eachindex(τ)
        μa_dynamic = μa + tmp * τ[ind]
        g1[ind] = integrate_compute_y_first(t -> (fluence_DA_Nlay_cylinder_TD(t, ρ, μa_dynamic, μsp, n_ext, n_med, l, a, z, bessels, N = N_laplace) / tpsf_norm), x, w, t[1], t[2])
    end

    return g1
end

function g2_DA_Nlay_cylinder_CW(τ::AbstractVector, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:500], β = 1.0)
    g1 = g1_DA_Nlay_cylinder_CW(τ, ρ, μa, μsp, BFi = BFi, n_ext = n_ext, n_med = n_med, l = l, a = a, z = z, λ = λ, bessels = bessels)
    for ind in eachindex(τ)
        g1[ind] = 1 + β * abs(g1[ind])^2
    end
    return g1
end
function g2_DA_Nlay_cylinder_CW(τ::AbstractVector, data::Nlayer_cylinder_DCS)
    return g2_DA_Nlay_cylinder_CW(τ, data.ρ, data.μa, data.μsp, BFi = data.BFi, n_ext = data.n_ext, n_med = data.n_med, l = data.l, a = data.a, λ = data.λ, bessels = data.bessels, β = data.β)
end

function g2_DA_Nlay_cylinder_TD(τ::AbstractVector, t::Union{AbstractVector, AbstractFloat}, ρ, μa, μsp; BFi = [2e-8, 2e-8], n_ext = 1.0, n_med = [1.0, 1.0], l = [1.0, 10.0], a = 20.0, z = 0.0, λ = 700.0, bessels = besselroots[1:500], N_laplace = 12, N_quad = 50, β = 1.0)
    if t isa AbstractVector
        g1 = g1_DA_Nlay_cylinder_TD(τ, t, ρ, μa, μsp, BFi = BFi, n_ext = n_ext, n_med = n_med, l = l, a = a, z = z, λ = λ, bessels = bessels, N_laplace = N_laplace, N_quad = N_quad)
    elseif t isa AbstractFloat
        g1 = g1_DA_Nlay_cylinder_TD(τ, t, ρ, μa, μsp, BFi = BFi, n_ext = n_ext, n_med = n_med, l = l, a = a, z = z, λ = λ, bessels = bessels, N_laplace = N_laplace)
    end

    for ind in eachindex(τ)
        g1[ind] = 1 + β * abs(g1[ind])^2
    end
    return g1
end
function g2_DA_Nlay_cylinder_TD(τ::AbstractVector, t, data::Nlayer_cylinder_DCS)
    return g2_DA_Nlay_cylinder_TD(τ, t, data.ρ, data.μa, data.μsp, BFi = data.BFi, n_ext = data.n_ext, n_med = data.n_med, l = data.l, a = data.a, λ = data.λ, bessels = data.bessels, N_laplace = data.N_laplace, N_quad = data.N_quad, β = data.β)
end


text:
module HWconstrained

greet() = print("Hello World!")

	using JuMP, NLopt, DataFrames, Ipopt
	using LinearAlgebra

	export data, table_NLopt, table_JuMP, obj, constr, max_NLopt, max_JuMP

	function data(a=0.5)
		n = 3
		p = [1, 1, 1]
		e = [2, 0, 0]
		ns = 4
		nss = ns^2
		z1 = [0.72, 0.92, 1.12, 1.32]
		z2 = [0.86, 0.96, 1.06, 1.16]
		za = [[1, i, j] for i in z1 for j in z2]
		z = vcat(za'...)
		pi = repeat([1/nss], nss)
		return Dict("a"=>a,"ns"=>ns,"nss"=>nss,"e"=>e,"p"=>p,"z"=>z,"pi"=>pi, "za"=>za)
	end

	function max_JuMP(a=0.5)
		d = data(a)
		m = Model(with_optimizer(Ipopt.Optimizer))
		@variable(m, c >= 0)
		@variable(m, x[1:3])
		@NLconstraint(m, c + sum(d["p"][i] * (x[i] - d["e"][i]) for i in 1:3) == 0)
		@NLobjective(m, Max, -exp(-a * c) + sum(-exp(-a * sum(d["za"][j][i] * x[i] for i in 1:3)) * d["pi"][j] for j in 1:16))
		JuMP.optimize!(m)
		return Dict("obj"=>objective_value(m),"c"=>value(c),"omegas"=>[value(x[i]) for i in 1:length(x)])
	end

	function table_JuMP()
		d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3))
		for i in 1:nrow(d)
			xx = max_JuMP(d[i,:a])
			d[i,:c] = xx["c"]
			d[i,:omega1] = xx["omegas"][1]
			d[i,:omega2] = xx["omegas"][2]
			d[i,:omega3] = xx["omegas"][3]
			d[i,:fval] = xx["obj"]
		end
		return d
	end

	function obj(x::Vector,grad::Vector,data::Dict)
		a = data["a"]
		z = data["z"]
		u(x) = -exp(-a * x)
		up(x) = a * exp(-a * x)
		if length(grad) > 0
			grad[1] = up(x[1])
			grad[2:end] = sum(data["pi"] .* z .* up.(z * x[2:end]), dims=1)
		end
		return u(x[1]) + sum(data["pi"] .* u.(z * x[2:end]))
	end

	function constr(x::Vector,grad::Vector,data::Dict)
		if length(grad) > 0
			grad[1] = 1
			grad[2:end] = data["p"]
		end
		return x[1] + sum(data["p"] .* (x[2:end] - data["e"]))
	end

	function max_NLopt(a=0.5)
		d = data(a)
		opt = Opt(:LD_MMA, 4)
		lower_bounds!(opt, [0., -Inf, -Inf, -Inf])
		max_objective!(opt, (x, g)->obj(x, g, d))
		inequality_constraint!(opt, (x, g)->constr(x, g, d))
		ftol_rel!(opt, 1e-9)
		NLopt.optimize(opt, vcat(0, d["e"]))
	end

	function table_NLopt()
		d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3))
		for i in 1:nrow(d)
			xx = max_NLopt(d[i,:a])
			for j in 2:ncol(d)-1
				d[i,j] = xx[2][j-1]
			end
			d[i,end] = xx[1]
		end
		return d
	end

end # module


text:
#import CUDA: _cuview, ViewIndex, NonContiguous
#using CUDA: genperm
# fallback to SubArray when the view is not contiguous

#=
function LinearAlgebra.permutedims!(dest::GPUArray, src::GPUArray, perm) where N
    perm isa Tuple || (perm = Tuple(perm))
    gpu_call(dest, (dest, src, perm)) do state, dest, src, perm
        I = @cartesianidx src state
        @inbounds dest[genperm(I, perm)...] = src[I...]
        return
    end
    return dest
end
=#

import CUDA: pow, abs, angle
for (RT, CT) in [(:Float64, :ComplexF64), (:Float32, :ComplexF32)]
    @eval cp2c(d::$RT, a::$RT) = CUDA.Complex(d*CUDA.cos(a), d*CUDA.sin(a))
    for NT in [RT, :Int32]
        @eval CUDA.pow(z::$CT, n::$NT) = CUDA.Complex((CUDA.pow(CUDA.abs(z), n)*CUDA.cos(n*CUDA.angle(z))), (CUDA.pow(CUDA.abs(z), n)*CUDA.sin(n*CUDA.angle(z))))
    end
end

@inline function bit_count(x::UInt32)
    x = ((x >> 1) & 0b01010101010101010101010101010101) + (x & 0b01010101010101010101010101010101)
    x = ((x >> 2) & 0b00110011001100110011001100110011) + (x & 0b00110011001100110011001100110011)
    x = ((x >> 4) & 0b00001111000011110000111100001111) + (x & 0b00001111000011110000111100001111)
    x = ((x >> 8) & 0b00000000111111110000000011111111) + (x & 0b00000000111111110000000011111111)
    x = ((x >> 16)& 0b00000000000000001111111111111111) + (x & 0b00000000000000001111111111111111)
    return x
end

@inline function bit_count(x::Int32)
    x = ((x >> 1) & Int32(0b01010101010101010101010101010101)) + (x & Int32(0b01010101010101010101010101010101))
    x = ((x >> 2) & Int32(0b00110011001100110011001100110011)) + (x & Int32(0b00110011001100110011001100110011))
    x = ((x >> 4) & Int32(0b00001111000011110000111100001111)) + (x & Int32(0b00001111000011110000111100001111))
    x = ((x >> 8) & Int32(0b00000000111111110000000011111111)) + (x & Int32(0b00000000111111110000000011111111))
    x = ((x >> 16)& Int32(0b00000000000000001111111111111111)) + (x & Int32(0b00000000000000001111111111111111))
    return x
end

bit_count(UInt32(0b11111))

# TODO
# support norm(view(reshape(A, m, n), :, 1))
using LinearAlgebra
import LinearAlgebra: norm
norm2(A::DenseCuArray; dims=1) = mapreduce(abs2, +, A, dims=dims) .|> CUDA.sqrt

export piecewise, cudiv
@inline function cudiv(x::Int)
    max_threads = 256
    threads_x = min(max_threads, x)
    threads_x, ceil(Int, x/threads_x)
end

# NOTE: the maximum block size is 65535
@inline function cudiv(x::Int, y::Int)
    max_threads = 256
    threads_x = min(max_threads, x)
    threads_y = min(max_threads ÷ threads_x, y)
    threads = (threads_x, threads_y)
    blocks = ceil.(Int, (x, y) ./ threads)
    threads, blocks
end

fix_cudiv(A::AbstractVector, firstdim::Int) = cudiv(firstdim)
fix_cudiv(A::AbstractMatrix, firstdim::Int) = cudiv(firstdim, size(A,2))

piecewise(state::AbstractVector, inds) = state
piecewise(state::AbstractMatrix, inds) = @inbounds view(state,:,inds[2])

import Base: kron, getindex
function kron(A::DenseCuArray{T1}, B::DenseCuArray{T2}) where {T1, T2}
    res = CUDA.zeros(promote_type(T1,T2), (size(A).*size(B))...)
    CI = Base.CartesianIndices(res)
    @inline function kernel(res, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        @inbounds inds = CI[state].I
        inds_A = (inds.-1) .÷ size(B) .+ 1
        inds_B = (inds.-1) .% size(B) .+ 1
        state <= length(res) && (@inbounds res[state] = A[inds_A...]*B[inds_B...])
        return
    end

    X, Y = cudiv(length(res))
    @cuda threads=X blocks=Y kernel(res, A, B)
    res
end

"""
    kron!(C::CuArray, A, B)

Computes Kronecker products in-place on the GPU.
The results are stored in 'C', overwriting the existing values of 'C'.
"""
function kron!(C::CuArray{T3}, A::DenseCuArray{T1}, B::DenseCuArray{T2}) where {T1, T2, T3}
    @boundscheck (size(C) == (size(A,1)*size(B,1), size(A,2)*size(B,2))) || throw(DimensionMismatch())
    CI = Base.CartesianIndices(C)
    @inline function kernel(C, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        @inbounds inds = CI[state].I
        inds_A = (inds.-1) .÷ size(B) .+ 1
        inds_B = (inds.-1) .% size(B) .+ 1
        state <= length(C) && (@inbounds C[state] = A[inds_A...]*B[inds_B...])
        return
    end

    X, Y = cudiv(length(C))
    @cuda threads=X blocks=Y kernel(C, A, B)
    C
end

function getindex(A::DenseCuVector{T}, B::DenseCuArray{<:Integer}) where T
    res = CUDA.zeros(T, size(B)...)
    @inline function kernel(res, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        state <= length(res) && (@inbounds res[state] = A[B[state]])
        return
    end

    X, Y = cudiv(length(B))
    @cuda threads=X blocks=Y kernel(res, A, B)
    res
end

function getindex(A::AbstractVector, B::DenseCuArray{<:Integer})
    A[Array(B)]
end

YaoBlocks.AD.as_scalar(x::DenseCuArray) = Array(x)[]


text:
using IOIndents
using Base.Test
using TestSetExtensions

function test_ioindent(input, output)
    iob = IOBuffer()
    io = IOIndent(iob)
    write(io, input...)
    @test String(take!(iob)) == output
end

@testset ExtendedTestSet "Tests" begin

@testset ExtendedTestSet "nops" begin
    test_ioindent("abc", "abc")
    test_ioindent("\nabc", "\nabc")
    test_ioindent("\nabc\n", "\nabc\n")
    test_ioindent("\nabc\n ", "\nabc\n ")

    test_ioindent(("abc", Indent()), "abc")
    test_ioindent(("abc", Align(), "\n"), "abc\n")
end

@testset ExtendedTestSet "Aligns" begin
    test_ioindent(("abc", Align(), "\nde"),
    """
    abc
       de""")

    test_ioindent(("αβγ", Align(), "\nϵκ"),
    """
    αβγ
       ϵκ""")

    test_ioindent(("abc", Align(), "\nde", Align(), "fg\nhi", "\nqr",
                   Dealign(), "\njk\n", Dealign(), "lm"),
    """
    abc
       defg
         hi
         qr
       jk
    lm""")

    test_ioindent(("abc", Align(), "\nde", Align(), "fg\nhi",
                   Dealign(), "\njk\n", Dealign(), "lm"),
    """
    abc
       defg
         hi
       jk
    lm""")

    test_ioindent((
    "\\begin{axis}[", Align(),
                   "\na = 5"),
    raw"""
    \begin{axis}[
                 a = 5""")
end

@testset ExtendedTestSet "Indents" begin
    test_ioindent(("abc\n", Indent(), Indent(), "cde", Dedent(),
                    "\nfgh\n", Dedent(), "ijk"),
    raw"""
    abc
            cde
        fgh
    ijk""")
end  # testset

@testset ExtendedTestSet "Mixed" begin
    test_ioindent((
    "\\begin{axis}[", Align(),
                   "a = 5,\nb = 5,", "\n",
                   "c = 5]",
        Dealign(), "\n",
        Indent(),
        "\\begin{plot}[", Align(),
                       "a = 5,\n", "b = 5,\n",
                        "c = 5]", Dealign(), "\n",
        "foo", Align(),
            "\nbar", Dealign(), Dedent(), Dedent(),
    "\nbase"),
    raw"""
    \begin{axis}[a = 5,
                 b = 5,
                 c = 5]
        \begin{plot}[a = 5,
                     b = 5,
                     c = 5]
        foo
           bar
    base""")
end

@testset ExtendedTestSet "README" begin

    test_ioindent(("[", Align(),
      "a = 1,",
      "\nb = ", Align(), "[",
      "c = 3,",
      "\n d = 4]", Dealign(),
      "\ne = 5]", Dealign()),
    raw"""
    [a = 1,
     b = [c = 3,
          d = 4]
     e = 5]""")

    test_ioindent(("function f(", Align(),
                          "a_long_argument,",
                          "\nanother_argument;",
                          "\na_keyword = :default", Dealign(),
                          "\nend"),
    raw"""
    function f(a_long_argument,
               another_argument;
               a_keyword = :default
    end""")

    test_ioindent(("function f(x)\n", Indent(),
                     "if x == 1\n", Indent(),
                     "return x\n", Dedent(),
                     "end\n", Dedent(),
                     "end"),
    raw"""
    function f(x)
        if x == 1
            return x
        end
    end""")
end

@testset ExtendedTestSet "IOContext" begin
    io = IOIndent(IOContext(IOBuffer(), :foo => true))
    @test get(io, :foo, false) == true
    @test get(io, :foobar, false) == false
    @test haskey(io, :foo) == true
    @test haskey(io, :foobar) == false
end

@testset ExtendedTestSet "nop normal IO" begin
    io = IOBuffer()
    print(io, "foo", Indent(), Dedent(), "bar\n", Indent(), "baz")
    @test String(take!(io)) == "foobar\nbaz"
end

@testset ExtendedTestSet "show" begin
    io = IOBuffer()
    show(io, IOIndent(IOBuffer()))
    @test String(take!(io)) == """
    IOIndent:
        IO: IOBuffer(data=UInt8[...], readable=true, writable=true, seekable=true, append=false, size=0, maxsize=Inf, ptr=1, mark=-1)
        Indent string: "    "
        Align char: " "
        Indent: 0
        Aligns: """
end

end


text:
# p33.jl - solve linear BVP u_xx = exp(4x), u'(-1)=u(1)=0

N = 16;
(D,x) = cheb(N);
D2 = D^2;
D2[N+1,:] = D[N+1,:];            # Neumann condition at x = -1
D2 = D2[2:N+1,2:N+1];
f = @. exp(4*x[2:N]);
u = D2\[f;0];
u = [0;u];
clf();
plot(x,u,".",markersize=10);
axis([-1,1,-4,0]);
xx = -1:.01:1;
uu = polyval(polyfit(x,u),xx);
plot(xx,uu); grid(true);
exact = @. (exp(4*xx) - 4*exp(-4)*(xx-1) - exp(4))/16;
title("max err = $(signif(norm(uu-exact,Inf),5))",fontsize=12);


text:
using InfiniteLinearAlgebra, InfiniteArrays, BandedMatrices, FillArrays, LazyBandedMatrices, LazyArrays, Test
import BandedMatrices: _BandedMatrix

@testset "∞-banded" begin
    @testset "Diagonal and BandedMatrix" begin
        D = Diagonal(Fill(2,∞))

        B = D[1:∞,2:∞]
        @test B isa BandedMatrix
        @test B[1:10,1:10] == diagm(-1 => Fill(2,9))
        @test B[1:∞,2:∞] isa BandedMatrix

        A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0,∞), -1 => Fill(3.0,∞))
        x = [1; 2; zeros(∞)]
        @test A*x isa Vcat
        @test (A*x)[1:10] == A[1:10,1:10]*x[1:10]

        @test InfBandCartesianIndices(0)[1:5] == CartesianIndex.(1:5,1:5)
        @test InfBandCartesianIndices(1)[1:5] == CartesianIndex.(1:5,2:6)
        @test InfBandCartesianIndices(-1)[1:5] == CartesianIndex.(2:6,1:5)

        @test A[band(0)][2:10] == 2:10
        @test D[band(0)] ≡ Fill(2,∞)
        @test D[band(1)] ≡ Fill(0,∞)

        @test B*A*x isa Vcat
        @test (B*A*x)[1:10] == [0; 10; 14; 12; zeros(6)]

        @test _BandedMatrix((1:∞)', ∞, -1, 1) isa BandedMatrix
    end

    @testset "∞-Toeplitz" begin
        A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
        @test A isa InfToeplitz
        @test MemoryLayout(typeof(A.data)) == ConstRows()
        @test MemoryLayout(typeof(A)) == BandedToeplitzLayout()
        V = view(A,:,3:∞)
        @test MemoryLayout(typeof(bandeddata(V))) == ConstRows()
        @test MemoryLayout(typeof(V)) == BandedToeplitzLayout()
        @test BandedMatrix(V) isa InfToeplitz
        @test A[:,3:end] isa InfToeplitz

        @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
        @test (A*A)[1:10,1:10] ≈ A[1:10,1:16]A[1:16,1:10]
        @test (A * Fill(2,∞))[1:10] ≈ 2A[1:10,1:16]*ones(16)
        @test (Fill(2,∞,∞)*A)[1:10,1:10] ≈ fill(2,10,13)A[1:13,1:10]

        @test Eye(∞) * A isa BandedMatrix
        @test A * Eye(∞) isa BandedMatrix

        @test A * [1; 2; Zeros(∞)] isa Vcat
        @test A * [1; 2; Zeros(∞)] == [A[1:5,1:2] * [1,2]; Zeros(∞)]

        @test MemoryLayout(Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
        @test MemoryLayout(Bidiagonal(Fill(1,∞), Fill(2,∞), :U)) isa BidiagonalToeplitzLayout
        @test MemoryLayout(SymTridiagonal(Fill(1,∞), Fill(2,∞))) isa TridiagonalToeplitzLayout
        @test MemoryLayout(LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
        @test MemoryLayout(LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)) isa BidiagonalToeplitzLayout
        @test MemoryLayout(LazyBandedMatrices.SymTridiagonal(Fill(1,∞), Zeros(∞))) isa TridiagonalToeplitzLayout

        T = LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))
        @test T[2:∞,3:∞] isa SubArray
        @test exp.(T) isa BroadcastMatrix
        @test exp.(T)[2:∞,3:∞] isa SubArray

        B = LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)
        @test B[2:∞,3:∞] isa SubArray
        @test exp.(B) isa BroadcastMatrix
        @test exp.(B)[2:∞,3:∞] isa SubArray

        @testset "algebra" begin
            T = Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))
            @test T isa InfiniteLinearAlgebra.TriToeplitz
            @test (T + 2I)[1:10,1:10] == (2I + T)[1:10,1:10] == T[1:10,1:10] + 2I
        end

        @testset "constant data" begin
            A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
            B = _BandedMatrix(Fill(2,4,∞), ∞, 1,2)
            @test (B*B)[1:10,1:10] ≈ B[1:10,1:14]B[1:14,1:10]
            @test (A*B)[1:10,1:10] ≈ A[1:10,1:14]B[1:14,1:10]
            @test (B*A)[1:10,1:10] ≈ B[1:10,1:14]A[1:14,1:10]
        end
    end

    @testset "Pert-Toeplitz" begin
        @testset "Inf Pert" begin
            A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
            @test A isa PertToeplitz
            @test MemoryLayout(typeof(A)) == PertToeplitzLayout()
            V = view(A,2:∞,2:∞)
            @test MemoryLayout(typeof(V)) == PertToeplitzLayout()
            @test BandedMatrix(V) isa PertToeplitz
            @test A[2:∞,2:∞] isa PertToeplitz

            @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I

            @test Eye(∞) * A isa BandedMatrix
            @test A * Eye(∞) isa BandedMatrix
        end

        @testset "TriPert" begin
            A = SymTridiagonal(Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞)))
            @test A isa InfiniteLinearAlgebra.SymTriPertToeplitz
            @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I

            A = Tridiagonal(Vcat([3.,4.], Fill.(0.5,∞)), Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞)))
            @test A isa InfiniteLinearAlgebra.TriPertToeplitz
            @test Adjoint(A) isa InfiniteLinearAlgebra.AdjTriPertToeplitz
            @test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
            @test (Adjoint(A) + 2I)[1:10,1:10] == (2I + Adjoint(A))[1:10,1:10] == Adjoint(A)[1:10,1:10] + 2I
        end

        @testset "InfBanded" begin
            A = _BandedMatrix(Fill(2,4,∞),ℵ₀,2,1)
            B = _BandedMatrix(Fill(3,2,∞),ℵ₀,-1,2)
            @test mul(A,A) isa PertToeplitz
            @test A*A isa PertToeplitz
            @test (A*A)[1:20,1:20] == A[1:20,1:23]*A[1:23,1:20]
            @test (A*B)[1:20,1:20] == A[1:20,1:23]*B[1:23,1:20]
        end
    end

    @testset "adjortrans" begin
        A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0+im,∞), -1 => Fill(3.0,∞))
        @test copy(A')[1:10,1:10] == (A')[1:10,1:10]
        @test copy(transpose(A))[1:10,1:10] == transpose(A)[1:10,1:10]
    end

    @testset "Eye subindex" begin
        @test Eye(∞)[:,1:3][1:5,:] == Eye(∞)[Base.Slice(oneto(∞)),1:3][1:5,:] == Eye(5,3)
        @test Eye(∞)[1:3,:][:,1:5] == Eye(∞)[1:3,Base.Slice(oneto(∞))][:,1:5] == Eye(3,5)
        @test Eye(∞)[:,:][1:5,1:3] == Eye(∞)[Base.Slice(oneto(∞)),Base.Slice(oneto(∞))][1:5,1:3] == Eye(5,3)
    end

    @testset "band(0) indexing" begin
        D = ApplyArray(*, Diagonal(1:∞), Diagonal(1:∞))
        @test D[band(0)][1:10] == (1:10).^2
    end

    @testset "Fill * Banded" begin
        A = _BandedMatrix(Ones(1,∞), ∞, 1,-1)
        B = _BandedMatrix(Fill(1.0π,1,∞), ∞, 0,0)
        @test (A*B)[1:10,1:10] ≈ BandedMatrix(-1 => Fill(1.0π,9))
    end

    @testset "Diagonal{Fill} * Bidiagonal" begin
        A, B = Diagonal(Fill(2,∞)) , LazyBandedMatrices.Bidiagonal(exp.(1:∞), exp.(1:∞), :L)
        @test (A*B)[1:10,1:10] ≈ (B*A)[1:10,1:10] ≈ 2B[1:10,1:10]
    end

    @testset "concat" begin
        H = ApplyArray(hvcat, 2, 1, [1 Zeros(1,∞)], [1; Zeros(∞)], Diagonal(1:∞))
        @test bandwidths(H) == (1,1)
        H = ApplyArray(hvcat, 2, 1, [0 Zeros(1,∞)], [0; Zeros(∞)], Diagonal(1:∞))
        @test bandwidths(H) == (0,0)
        H = ApplyArray(hvcat, (2,2), 1, [1 Zeros(1,∞)], [1; Zeros(∞)], Diagonal(1:∞))
        @test_broken bandwidths(H) == (1,1)
    end

    @testset "Banded * PaddedMatrix" begin
        A = Eye(∞)[2:∞,:]
        B = PaddedArray(randn(3,3),ℵ₀,ℵ₀)
        @test (A*B)[1:10,1:10] ≈ A[1:10,1:10] * B[1:10,1:10]
    end
end


text:
using Base.Test, DataFrames, MixedModels

contra = readtable(joinpath(dirname(@__FILE__), "data", "Contraception.csv.gz"))
for c in [:district, :use, :urban, :livch, :urbdist]
    contra[c] = pool(contra[c])
end
contra[:age2] = abs2(contra[:age])
gm1 = fit!(glmm(use01 ~ 1 + age + age2 + urban + livch + (1 | urbdist), contra, Binomial()));

@test isapprox(LaplaceDeviance(gm1), 2361.5457541; atol = 0.0001)
# There may be multiple optima here.
#@test isapprox(logdet(gm1), 75.7204822; atol = 0.0001)
#@test isapprox(sumabs2(gm1.u[1]), 48.47486965; atol = 0.0001)
#@test isapprox(sum(gm1.devresid), 2237.3504024; atol = 0.0001)

cbpp = readtable(joinpath(dirname(@__FILE__), "data", "cbpp.csv.gz"))
for c in [:herd, :period]
    cbpp[c] = pool(oftype(Int8[], cbpp[c]))
end
for c in [:incidence, :size]
  cbpp[c] = oftype(Int8[], cbpp[c])
end
cbpp[:prop] = cbpp[:incidence] ./ cbpp[:size]
gm2 = fit!(glmm(prop ~ 1 + period + (1 | herd), cbpp, Binomial(), LogitLink(); wt = cbpp[:size]));

@test isapprox(LaplaceDeviance(gm2), 100.095856; atol = 0.0001)
@test isapprox(sumabs2(gm2.u[1]), 9.72305; atol = 0.0001)
@test isapprox(logdet(gm2), 16.901054; atol = 0.0001)
@test isapprox(sum(gm2.devresid), 73.47143; atol = 0.001)
#VerbAgg = rcopy("lme4::VerbAgg")
#VerbAgg[:r201] = [Int8(x == "N" ? 0 : 1) for x in VerbAgg[:r2]]
#m3 = glmm(r201 ~ 1 + Anger + Gender + btype + situ + (1 | id) + (1 | item), VerbAgg, Bernoulli());
#f = r201 ~ 1 + Anger + Gender + btype + situ + (1 | id)

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
"""
phllf
"""

function phllf(cf1::CF1{Tv},
    data::AbstractPHSample, ::Type{MatT} = SparseMatrixCSC; eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01)) where {Tv,MatT}
    phllf(GPH(cf1, MatT), data, eps=eps, ufact=ufact)
end

function phllf(cf1::CF1{Tv}, cf1deriv::Dict{Symbol,CF1{Tv}},
    data::AbstractPHSample, ::Type{MatT} = SparseMatrixCSC; eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01)) where {Tv,MatT}
    gph = GPH(cf1, MatT)
    deriv = Dict((k, GPH(v, MatT)) for (k,v) in cf1deriv)
    phllf(gph, deriv, data, eps=eps, ufact=ufact)
end

function phllf(cf1::CF1{Tv}, cf1deriv::Dict{Symbol,CF1{Tv}}, cf1deriv2::Dict{Tuple{Symbol,Symbol},CF1{Tv}},
    data::AbstractPHSample, ::Type{MatT} = SparseMatrixCSC; eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01)) where {Tv,MatT}
    gph = GPH(cf1, MatT)
    deriv = Dict((k, GPH(v, MatT)) for (k,v) in cf1deriv)
    deriv2 = Dict((k, GPH(v, MatT)) for (k,v) in cf1deriv2)
    phllf(gph, deriv, deriv2, data, eps=eps, ufact=ufact)
end

function phllf(ph::GPH{Tv,MatT},
    data::WeightedSample{Tv};
    eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01)) where {Tv,MatT}
    dim, alpha, P, tau, qv = _phunif(ph, ufact)
    llf = Tv(0)

    right = rightbound(qv * data.maxtime, eps) + 1
    prob = Vector{Tv}(undef, right + 1)
    f0 = copy(alpha)
    f1 = similar(alpha)

    for k = eachindex(data.tdat)
        f1 .= Tv(0)
        unifforward!(f0, data.tdat[k], P, qv, prob, eps, f1)
        llf += data.wdat[k] * log(@dot(f1, tau))
        f0 .= f1
    end
    return llf
end

function phllf(ph::GPH{Tv,MatT}, phderiv::Dict{Symbol,GPH{Tv,MatT}},
    data::WeightedSample{Tv};
    eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01)) where {Tv,MatT}
    dim, alpha, P, tau, qv = _phunif(ph, ufact)
    Pdash = Dict((k,x.T/qv) for (k,x) in phderiv)
    llf = Tv(0)
    llfdash = Dict((k,Tv(0)) for (k,x) in phderiv)

    right = rightbound(qv * data.maxtime, eps) + 1
    prob = Vector{Tv}(undef, right + 1)
    f0 = copy(alpha)
    f1 = similar(alpha)
    f0dash = Dict((k,copy(x.alpha)) for (k,x) in phderiv)
    f1dash = Dict((k,similar(x.alpha)) for (k,x) in phderiv)

    for k = eachindex(data.tdat)
        f1 .= Tv(0)
        for j = eachindex(phderiv)
            f1dash[j] .= Tv(0)
        end
        unifforward!(f0, f0dash, data.tdat[k], P, Pdash, qv, prob, eps, f1, f1dash)
        tmp = @dot(f1, tau)
        llf += data.wdat[k] * log(tmp)
        for i = eachindex(phderiv)
            llfdash[i] += data.wdat[k] * (@dot(f1dash[i], tau) + @dot(f1, phderiv[i].tau)) / tmp
        end
        f0 .= f1
        for j = eachindex(phderiv)
            f0dash[j] .= f1dash[j]
        end
    end
    return llf, llfdash
end

function phllf(ph::GPH{Tv,MatT}, phderiv::Dict{Symbol,GPH{Tv,MatT}}, phderiv2::Dict{Tuple{Symbol,Symbol},GPH{Tv,MatT}},
    data::WeightedSample{Tv};
    eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01)) where {Tv,MatT}
    dim, alpha, P, tau, qv = _phunif(ph, ufact)
    Pdash = Dict((k,x.T/qv) for (k,x) in phderiv)
    Pdashdash = Dict((k,x.T/qv) for (k,x) in phderiv2)
    llf = Tv(0)
    llfdash = Dict((k,Tv(0)) for (k,x) in phderiv)
    tmp1 = Dict((k,Tv(0)) for (k,x) in phderiv)
    llfdashdash = Dict((k,Tv(0)) for (k,x) in phderiv2)

    right = rightbound(qv * data.maxtime, eps) + 1
    prob = Vector{Tv}(undef, right + 1)
    f0 = copy(alpha)
    f1 = similar(alpha)
    f0dash = Dict((k,copy(x.alpha)) for (k,x) in phderiv)
    f1dash = Dict((k,similar(x.alpha)) for (k,x) in phderiv)
    f0dashdash = Dict((k,copy(x.alpha)) for (k,x) in phderiv2)
    f1dashdash = Dict((k,similar(x.alpha)) for (k,x) in phderiv2)

    for k = eachindex(data.tdat)
        f1 .= Tv(0)
        for i = eachindex(f1dash)
            f1dash[i] .= Tv(0)
        end
        for i = eachindex(f1dashdash)
            f1dashdash[i] .= Tv(0)
        end
        unifforward!(f0, f0dash, f0dashdash, data.tdat[k],
            P, Pdash, Pdashdash, qv, prob, eps, f1, f1dash, f1dashdash)

        tmp = @dot(f1, tau)
        llf += data.wdat[k] * log(tmp)
        for i = eachindex(f1dash)
            tmp1[i] = (@dot(f1dash[i], tau) + @dot(f1, phderiv[i].tau)) / tmp
            llfdash[i] += data.wdat[k] * tmp1[i]
        end
        for i = eachindex(f1dashdash)
            k1 = i[1]
            k2 = i[2]
            llfdashdash[i] += data.wdat[k] * ((@dot(f1dashdash[i], tau) + @dot(f1dash[k1], phderiv[k2].tau) +
                            @dot(f1dash[k2], phderiv[k1].tau) + @dot(f1, phderiv2[i].tau)) / tmp - tmp1[k1] * tmp1[k2])
        end
        f0 .= f1
        for i = eachindex(f1dash)
            f0dash[i] .= f1dash[i]
        end
        for i = eachindex(f1dashdash)
            f0dashdash[i] .= f1dashdash[i]
        end
    end
    return llf, llfdash, llfdashdash
end

"""
Uniformization
"""

@origin prob => left function unifforward!(f0::Vector{Tv}, t::Tv,
    P::AbstractMatrix{Tv}, qv::Tv,
    prob::Vector{Tv}, eps::Tv, f1::Vector{Tv}) where Tv
    left = 0
    right = rightbound(qv*t, eps)
    weight = poipmf!(qv*t, prob, left=left, right=right)
    f1 .= prob[left] * f0
    for u = left+1:right
        f0 .= P' * f0
        @axpy(prob[u], f0, f1)
    end
    @scal(1/weight, f1)
    nothing
end

@origin prob => left function unifforward!(
    f0::Vector{Tv}, f0dash::Dict{Symbol,Vector{Tv}}, t::Tv,
    P::AbstractMatrix{Tv}, Pdash::Dict{Symbol,<:AbstractMatrix{Tv}}, qv::Tv,
    prob::Vector{Tv}, eps::Tv,
    f1::Vector{Tv}, f1dash::Dict{Symbol,Vector{Tv}}) where Tv
    left = 0
    right = rightbound(qv*t, eps)
    weight = poipmf!(qv*t, prob, left=left, right=right)
    f1 .= prob[left] * f0
    for k = eachindex(f1dash)
        f1dash[k] .= prob[left] * f0dash[k]
    end
    for u = left+1:right
        for k = eachindex(f0dash)
            f0dash[k] .= P' * f0dash[k] + Pdash[k]' * f0
        end
        f0 .= P' * f0
        @axpy(prob[u], f0, f1)
        for k = eachindex(f0dash)
            @axpy(prob[u], f0dash[k], f1dash[k])
        end
    end
    @scal(1/weight, f1)
    for k = eachindex(f1dash)
        @scal(1/weight, f1dash[k])
    end
    nothing
end

@origin prob => left function unifforward!(
    f0::Vector{Tv}, f0dash::Dict{Symbol,Vector{Tv}}, f0dashdash::Dict{Tuple{Symbol,Symbol},Vector{Tv}}, t::Tv,
    P::AbstractMatrix{Tv}, Pdash::Dict{Symbol,<:AbstractMatrix{Tv}}, Pdashdash::Dict{Tuple{Symbol,Symbol},<:AbstractMatrix{Tv}}, qv::Tv,
    prob::Vector{Tv}, eps::Tv,
    f1::Vector{Tv}, f1dash::Dict{Symbol,Vector{Tv}}, f1dashdash::Dict{Tuple{Symbol,Symbol},Vector{Tv}}) where Tv
    left = 0
    right = rightbound(qv*t, eps)
    weight = poipmf!(qv*t, prob, left=left, right=right)
    f1 .= prob[left] * f0
    for k = eachindex(f0dash)
        f1dash[k] .= prob[left] * f0dash[k]
    end
    for k = eachindex(f0dashdash)
        f1dashdash[k] .= prob[left] * f0dashdash[k]
    end
    for u = left+1:right
        for k = eachindex(f0dashdash)
            k1 = k[1]
            k2 = k[2]
            f0dashdash[k] .= P' * f0dashdash[k] + Pdash[k2]' * f0dash[k1] + Pdash[k1]' * f0dash[k2] + Pdashdash[k]' * f0
        end
        for k = eachindex(f0dash)
            f0dash[k] .= P' * f0dash[k] + Pdash[k]' * f0
        end
        f0 .= P' * f0

        @axpy(prob[u], f0, f1)
        for k = eachindex(f0dash)
            @axpy(prob[u], f0dash[k], f1dash[k])
        end
        for k = eachindex(f0dashdash)
            @axpy(prob[u], f0dashdash[k], f1dashdash[k])
        end
    end
    @scal(1/weight, f1)
    for k = eachindex(f1dash)
        @scal(1/weight, f1dash[k])
    end
    for k = eachindex(f1dashdash)
        @scal(1/weight, f1dashdash[k])
    end
    nothing
end