text:
#!/usr/bin/env julia

# This file tests generation of multiblock VTK files.
# It also tests the non-default values of the optional parameters of vtk_grid,
# i.e., compress=false and append=false.

using WriteVTK
import Compat: UTF8String, view
const FloatType = Float32
const vtm_filename_noext = "multiblock"

function first_block_data()
    Ni, Nj, Nk = 10, 15, 20

    x = zeros(FloatType, Ni, Nj, Nk)
    y = copy(x)
    z = copy(x)

    # Just to test subarrays:
    qall = zeros(FloatType, 3Ni, 2Nj, 2Nk)
    q = view(qall, 1:Ni, 1:Nj, 1:Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r::FloatType = 1 + (i - 1)/(Ni - 1)
        th::FloatType = 3*pi/4 * (j-1) / (Nj-1)
        x[i, j, k] = r * cos(th)
        y[i, j, k] = r * sin(th)
        z[i, j, k] = (k-1) / (Nk-1)
        q[i, j, k] = k * exp(j) * sqrt(i)
    end

    return x, y, z, q
end

function second_block_data()
    Ni, Nj, Nk = 20, 16, 12

    x = zeros(FloatType, Ni)
    y = zeros(FloatType, Nj)
    z = zeros(FloatType, Nk)

    [x[i] =  1 +   (i-1) / (Ni-1) for i = 1:Ni]
    [y[j] = -2 + 2*(j-1) / (Nj-1) for j = 1:Nj]
    [z[k] =        (k-1) / (Nk-1) for k = 1:Nk]

    q = zeros(FloatType, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        q[i, j, k] = k * exp(j) * sqrt(i)
    end

    return x, y, z, q
end

function third_block_data()
    # This is basically the same as in the unstructured example.
    Ni, Nj, Nk = 40, 50, 20
    dims = (Ni, Nj, Nk)

    # Create points and point data.
    pts = Array{FloatType}(3, Ni, Nj, Nk)
    pdata = Array{FloatType}(Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        r = 1 + (i - 1)/(Ni - 1)
        th = 3*pi/4 * (j - 1)/(Nj - 1)
        pts[1, i, j, k] = r * cos(th)
        pts[2, i, j, k] = r * sin(th)
        pts[3, i, j, k] = -2 + (k - 1)/(Nk - 1)
        pdata[i, j, k] = i*i + k*sqrt(j)
    end

    # Create cells (all hexahedrons in this case) and cell data.
    celltype = VTKCellTypes.VTK_HEXAHEDRON
    cells = MeshCell[]
    cdata = FloatType[]

    for k = 2:Nk, j = 2:Nj, i = 2:Ni
        # Define connectivity of cell.
        inds = Array{Int32}(undef, 8)
        inds[1] = CartesianIndices(dims, i-1, j-1, k-1)
        inds[2] = CartesianIndices(dims, i  , j-1, k-1)
        inds[3] = CartesianIndices(dims, i  , j  , k-1)
        inds[4] = CartesianIndices(dims, i-1, j  , k-1)
        inds[5] = CartesianIndices(dims, i-1, j-1, k  )
        inds[6] = CartesianIndices(dims, i  , j-1, k  )
        inds[7] = CartesianIndices(dims, i  , j  , k  )
        inds[8] = CartesianIndices(dims, i-1, j  , k  )

        # Define cell.
        c = MeshCell(celltype, inds)

        push!(cells, c)
        push!(cdata, i*j*k)
    end

    return pts, cells, pdata, cdata
end

function main()
    # Generate data.
    x1, y1, z1, q1 = first_block_data()
    x2, y2, z2, q2 = second_block_data()
    points3, cells3, q3, c3 = third_block_data()

    # Create multiblock file.
    @time outfiles = vtk_multiblock(vtm_filename_noext) do vtm
        # First block.
        vtk = vtk_grid(vtm, x1, y1, z1; compress=false, append=false)
        vtk_point_data(vtk, q1, "q_values")

        # Second block.
        vtk = vtk_grid(vtm, x2, y2, z2; compress=false)
        vtk_point_data(vtk, q2, "q_values")

        # Third block.
        vtk = vtk_grid(vtm, points3, cells3; append=false)
        vtk_point_data(vtk, q3, "q_values")
        vtk_cell_data(vtk, c3, "c_values")
    end
    println("Saved:  ", join(outfiles, "  "))

    return outfiles::Vector{UTF8String}
end

ProjDir = dirname(@__FILE__)
cd(ProjDir*"/output") do

  main()

end


text:
using AutoLandmarking
using Flux
using CuArrays # to run models on the GPU for massive performance increase
import ImageView

push!(LOAD_PATH, Base.@__DIR__)

# Import data with custom import functions
image_data, imname = load_imgs("C:/Users/immanueldiener/Desktop/Master/scripts/no_enamel", (128,128,128), false)
landmark_data = read_landmarks("C:/Users/immanueldiener/Desktop/Master/scripts/no_enamel", 10, "@1")

# X and Y coordinates are swapped in avizo files, wo I'll swap them back here.
#this step is essential for all the data augmentation to work.
landmark_data = swap_xy(landmark_data)

#some images have a higher value in the fill voxels than in the actual object voxels. This will fix the issue.
# The highest value must always be the filler value!
image_data[:,:,:,[109,130]] = image_data[:,:,:,[109,130]] .- 1
image_data[:,:,:,[109,130]] = image_data[:,:,:,[109,130]] .* -1

image_data2, imname2 = load_imgs("C:/Users/immanueldiener/Desktop/Master/master_data\\additional_teeth_volumes\\", (128,128,128), true)
landmark_data2 = read_landmarks("C:/Users/immanueldiener/Desktop/Master/master_data\\additional_teeth_volumes\\", 10, "@1")

#Swap X and Y again
landmark_data2 = swap_xy(landmark_data2)

#all these images have reversed fill and object values. Tha highest value must always be the filler value!
image_data2 = (image_data2 .- 1) .* -1

names = vcat(imname, imname2)
images = cat(image_data, image_data2, dims=4)
lms = cat(landmark_data, landmark_data2, dims=2)

# align along principal axes
aligned, aligned_lms, retro = align_principal(images, lms, 192)
resized, resized_lms, scales = resize_relevant(aligned, aligned_lms, 144)
ImageView.imshow(depth_map_all_sides(resized))
#define training and testing sets
X_train, X_test, y_train, y_test = regular_train_test_split_3d(resized, resized_lms)

# perform data augmentation. In case of RAM limitations the depth-map generation can be
# done right after every augmentation step, after which the augmented volumes can be deleted.
X_train, y_train = (mirror_vol(X_train, y_train))

X_train = Float32.(X_train)

flip1, lm_flip1 = flip_3D(X_train, y_train)
flip2, lm_flip2 = flip_volume_front(flip1, lm_flip1)
flip3, lm_flip3 = flip_volume_side(flip2, lm_flip2)

rot, lm_rot = rotate_volumes(flip2, lm_flip2, 10)
rot2, lm_rot2 = rotate_volumes(X_train, y_train, 20)
rot3, lm_rot3 = rotate_volumes(flip3, lm_flip3, -15)

jit, lm_jit = jitter_3D(flip1, lm_flip1, 10)
jit2, lm_jit2 = jitter_3D(X_train, y_train, 10)
jit3, lm_jit3 = jitter_3D(flip3, lm_flip3, 10)
jit4, lm_jit4 = jitter_3D(rot2, lm_rot2, 10)

complete1 = depth_map_all_sides(X_train)
X_train=nothing
complete2 = depth_map_all_sides(flip1)
flip1=nothing
complete3 = depth_map_all_sides(flip2)
flip2=nothing
complete4 = depth_map_all_sides(flip3)
flip3=nothing
complete5 = depth_map_all_sides(rot)
rot=nothing
complete6 = depth_map_all_sides(rot2)
rot2=nothing
complete7 = depth_map_all_sides(rot3)
rot3=nothing
complete8 = depth_map_all_sides(jit)
jit=nothing
complete9 = depth_map_all_sides(jit2)
jit2=nothing
complete10 = depth_map_all_sides(jit3)
jit3=nothing
complete11 = depth_map_all_sides(jit4)

X_train = cat(complete1, complete2, complete3, complete4, complete5, complete6,
    complete7, complete8, complete9, complete10, complete11, dims=4)
y_train = cat(y_train, lm_flip1, lm_flip2, lm_flip3, lm_rot, lm_rot2, lm_rot3, lm_jit, lm_jit2, lm_jit3, lm_jit4, dims=2)

X_train = image_gradients(X_train)

# make depthmaps of the testing data
X_test = depth_map_all_sides(X_test)
X_test = image_gradients(X_test)

# define the cost function
cost(x, y) = sum((model(x)-y).^2)|>gpu

#define the model
model = Flux.mapleaves(cu, AutoLandmarking.vgg19)

# define the trainingrate and optimiser
opt = Flux.ADAM(0.00003)

# redefinition of the dropout function to work with testmode!() since this
# functionality is not working at the moment. This step only needs to be done
# if you wish to perform uncertainty estimation using the dropout method (AutoLM.response_distribution())
# after training.
using Random
function Flux.dropout(x, p; dims = :)
    q = 1 - p
    y = rand!(similar(x, Flux._dropout_shape(x, dims)))
    y .= Flux._dropout_kernel.(y, p, q)
    x .* y
end

import Zygote
Zygote.@adjoint function Flux.dropout(x, p; dims = :)
   q = 1 - p
   y = rand!(similar(x, Flux._dropout_shape(x, dims)))
   y .= Flux._dropout_kernel.(y, p, q)
   return x .* y, Δ -> (Δ .* y, nothing)
end
Zygote.refresh()

function run_model(modell, X, y)
  train_data = Flux.mapleaves(cu, X)
  train_labels = Flux.mapleaves(cu, y)
  dataset = Flux.Data.DataLoader(train_data, train_labels, batchsize = 4, shuffle=true)
  Flux.train!(cost, params(model), dataset, opt)
  testmode!(model)
  cosima = cost_whole_data_2D(train_data, train_labels, cost)
  testmode!(model, false)
  return cosima
end

# create lists to store metrics
accs = []
costs = []

# run the training for 300 epochs. batches of 128 will be loaded onto the GPU,
# which will be further subdevided into minibatches of 4 (as defined in AutoML.run_model).
for i in 1:400
  costr = 0
  for j in 1:128:size(X_train, 4)
    if j+127 > size(X_train, 4)
      train_data = X_train[:,:,:,j:end]
      train_labels = y_train[:,j:end]
    else
      train_data = X_train[:,:,:,j:j+127]
      train_labels = y_train[:,j:j+127]
    end
    cosima = run_model(model, train_data, train_labels)
    costr = costr + cosima
    train_data = nothing
    train_labels = nothing
    # GC.gc()
  end
  push!(costs, costr)
  println("epoch ", i, " finished")
  println(costr)
  testmode!(model)
  acc2, max1, min1 = AutoLandmarking.avg_accuracy_per_point(model, gpu(X_test[:,:,:,:]), gpu(y_test[:,:]), 3)
  testmode!(model, false)
  println("median deviation per point on testing dataset: ", acc2, "with maximum", max1)
  println("median sum is ", sum(acc2))
  push!(accs, sum(acc2))
  if sum(acc2)[1]<7.6
    break
  end
end

#save model
using BSON
testmode!(model) #save the model in testmode so that the dropout layers are inactive!
model_cpu = cpu(model)
BSON.@save "C:/Users/immanueldiener/Desktop/Master/master_data/model_teeth.bson" model_cpu

#load model back in
BSON.@load "C:/Users/immanueldiener/Desktop/Master/master_data/model_teeth.bson" model_cpu
model = gpu(model_cpu)

# ================== After training ===================

#solve issue with arrays inside vector
aks = []
for i in accs
  push!(aks, i[1])
end

# check what the best accuracy is
minimum(aks)

#plot the development of cost and accuracy
import Plots
Plots.plot(aks, legend=:topright, label="sum of deviations per point", color= :red, xlabel="epochs", ylabel = "acuracy")
plt = Plots.plot!(Plots.twinx(), costs, label="training loss", legend=:topleft, ylabel = "loss")

# check for outliers (greatest distance from the procrustes mean).
# teh idea is to eliminate out all the predicted landmark sets that
# deviate a lot from the procrustes mean (reference), since these are
# most likely bad predictions. do get the procrustes mean, we can use training
# testing data together to get a more accureate mean. This is an ideal approch
# because otherwise we would not realize if all the predictions are biased in some
# way that is congruent over all predictions but still wrong.

function print_costs(y, resp)
  costs=[]
  for i in 1:size(resp,2)
    dev = sum((y[:,i] .- resp[:,i]) .^2)
    push!(costs, dev)
    println(i, "    ", dev)
  end
  return costs
end

X_train2, X_test2, y_train2, y_test2 = regular_train_test_split_3d(images, lms) #to get the volumes again

response = cpu(predict_set(gpu(X_test), model)) #predictions

# translate predictions back to the original volumes and get the names for training
# and testing set. This is necessary since we performed PC-alignment and resizing before training.
retro_test = zeros(3,5,48)
scales_test = zeros(1,48)
names_test = []
names_train = []
for i in 1:48
  push!(names_test, names[i*5])
  retro_test[:,:,i] = retro[:,:,i*5]
  scales_test[1,i] = scales[i*5]
end
for i in 1:243
  if i%5!=0
    push!(names_train, names[i])
  end
end

response_scaled = deepcopy(response)

for i in 1:48
  response_scaled[:,i] .= response[:,i]./scales_test[i]
end

accuracies = []
for i in 1:48
  acc2, max1, min1 = AutoLandmarking.avg_accuracy_per_point(model, gpu(X_test[:,:,:,i:i]), gpu(y_test[:,i:i]), 3)
  push!(accuracies, sum(acc2))
  println(i, "   ", sum(acc2))
end

# set negative values to 0 and values outside the volume to the maximum along each dimension
response_unaligned = relu.(translate_lms_back(response_scaled, retro_test))
AutoLandmarking.change_values!(response_unaligned, 12.8, 12.8, >)

# put landmarks near the surface onto the surface and look at the cost change
response_on_vol = landmark_to_surface(X_test2, response_unaligned, 10)
costs = print_costs(y_test2, response_on_vol) # using the landmark_to_surface function is useful

#outlier detection based on procrustes distances

three_d_train = AutoLandmarking.to_3d_array(y_train2)
three_d_test = AutoLandmarking.to_3d_array(response_on_vol)

aligned_train = AutoLandmarking.align_all(three_d_train)
aligned_test = AutoLandmarking.align_all(three_d_test)

mean_train, stdevs = AutoLandmarking.mean_shape(aligned_train)

distances = AutoLandmarking.proc_distance(mean_train, aligned_test, stdevs)
dists = []
for i in 1:48
  push!(dists, distances[i])
end

# look at distances from mean and respective costs
for i in 1:48
  println("ind $i cost is ", costs[i], "  and dist is ", dists[i])
end

# outlier detection based on model uncertainty using the response_distribution function
uncertainties = []
for i in 1:48
  means, stddev = AutoLandmarking.response_distribution(model, gpu(X_test[:,:,:,i:i]), gpu(y_test[:,i:i]), 200)
  println("$i has uncevertaty", sum(stddev), "  and cost  ", costs[i])
  push!(uncertainties, sum(stddev))
end

ext: jl

max_stars_repo_name: imi-diener/AutoLM

max_issues_repo_name: imi-diener/AutoLM

max_forks_repo_name: imi-diener/AutoLM


text:
export expected_policy_values, nash_conv

function expected_policy_values(π::AbstractPolicy, env::AbstractEnv)
    if get_terminal(env)
        [get_reward(env, p) for p in get_players(env) if p != get_chance_player(env)]
    elseif get_current_player(env) == get_chance_player(env)
        vals = [0. for p in get_players(env) if p != get_chance_player(env)]
        for a::ActionProbPair in get_legal_actions(env)
            vals .+= a.prob .* expected_policy_values(π, child(env, a))
        end
        vals
    else
        vals = [0. for p in get_players(env) if p != get_chance_player(env)]
        actions = get_actions(env)
        probs = get_prob(π, env)
        @assert length(actions) == length(probs)

        for (a, p) in zip(actions, probs)
            if p > 0 #= ignore illegal action =#
                vals .+= p .* expected_policy_values(π, child(env, a))
            end
        end
        vals
    end
end

function nash_conv(π, env;is_reduce=true,kw...)
    e = copy(env)
    RLBase.reset!(e)

    σ′ = [
        best_response_value(BestResponsePolicy(π, e, i;kw...), e)
        for i in get_players(e)
        if i != get_chance_player(e)
    ]

    σ = expected_policy_values(π, e)
    if is_reduce
        mapreduce(-, +, σ′, σ)
    else
        map(-, σ′, σ)
    end
end


text:
# Only 2D for now, 3 features from two poses
# will add multi-modality soon

# using Distributions, TransformUtils

function getUvecScaleFeature2D(wTb, bTc, wTo)
  cTo = (wTb * bTc) \ wTo
  r = cTo[1:2,3]
  α = norm(r)
  return r, r/α, α
end

function getUvecScaleBaseline2D(wTbi, wTbj, bTc)
  ciTcj = (wTbi * bTc) \ (wTbj * bTc)
  b = ciTcj[1:2,3]
  β = norm(b)
  return b, b/β, β
end

# function getEstFeatDirectii(wTbi, wTobjs::Vector, bTc)
#   xirvecs = Vector{Tuple}()
#   for wTo in wTobjs
#     push!(xirvecs, getUvecScaleFeature2D(wTbi, bTc, wTo))
#   end
#   return xirvecs
# end

function residruvec()
  nothing
end

"""
$(TYPEDEF)

two sets of three feature sightings and a body to camera lever arm transform
"""
mutable struct MultipleFeatures2D <: IIF.AbstractRelativeMinimize
  # treat these as angles positive X->Y, X forward, Y left and Z up
  xir1::Normal
  xir2::Normal
  xir3::Normal
  xjr1::Normal
  xjr2::Normal
  xjr3a::Normal
  xjr3b::Normal
  abC::Categorical # OBSOLETE, use addFactor!(..., multihypo=[...]...) instead, see Caesar docs
  bTc::Array{Float64,2}
  zDim::Tuple{Int,Int,Int,Int,Int}
  # Takes five
  MultipleFeatures2D(a,b,c,d,e,f,g,h,i ) = new(a,b,c,d,e,f,g,h,i, (3,3,2,2,2) )
end

function getSample(mm2d::MultipleFeatures2D, N::Int=1)

  meas = Array{Float64,3}(undef, 3,N,2)
  meas[1,:,1] = rand(mm2d.xir1,N)
  meas[2,:,1] = rand(mm2d.xir2,N)
  meas[3,:,1] = rand(mm2d.xir3,N)
  meas[1,:,2] = rand(mm2d.xjr1,N)
  meas[2,:,2] = rand(mm2d.xjr2,N)
  meas[3,:,2] = rand(mm2d.xjr3a,N)

  select = rand(mm2d.abC,N)
  beb = select .== 2
  numb = sum(beb)
  meas[3,beb,2] = rand(mm2d.xjr3b,numb)

  return (meas, beb, numb)
end

# redo with angles on sightings for minimization
# NOTE -- an inefficient implementation
function (mm2d::MultipleFeatures2D)(res::Array{Float64},
          userdata,
          idx::Int, meas::Tuple,
          wAbi::Array{Float64,2},
          wAbj::Array{Float64,2},
          wAo1::Array{Float64,2},
          wAo2::Array{Float64,2},
          wAo3::Array{Float64,2},
          wAo3b::Union{Array{Float64,2},Bool}=false  )
  #
  len3b = wAo3b==false ? 0 : size(wAo3b,2)
  dobimodal = len3b > 0
  meas[3] > 0 && !dobimodal ? error("mm2d::MultipleFeatures2D -- bi-modal mismatch") : nothing

  wTbi = SE2(wAbi[:,idx])
  wTbj = SE2(wAbj[:,idx])
  wTo1 = SE2([wAo1[:,idx];0.0])
  wTo2 = SE2([wAo2[:,idx];0.0])
  wTo3 = SE2([wAo3[:,idx];0.0])
  wTo3b = dobimodal ? SE2([wAo3b[:,idx];0.0]) : nothing

  ciTo1 = (wTbi * mm2d.bTc) \ wTo1
  ciTo2 = (wTbi * mm2d.bTc) \ wTo2
  ciTo3 = (wTbi * mm2d.bTc) \ wTo3

  cjTo1 = (wTbj * mm2d.bTc) \ wTo1
  cjTo2 = (wTbj * mm2d.bTc) \ wTo2
  cjTo3 = (wTbj * mm2d.bTc) \ wTo3
  cjTo3b = dobimodal ? (wTbj * mm2d.bTc) \ wTo3 : nothing

  the1 = atan(se2vee(ciTo1)[2], se2vee(ciTo1)[1])
  the2 = atan(se2vee(ciTo2)[2], se2vee(ciTo2)[1])
  the3 = atan(se2vee(ciTo3)[2], se2vee(ciTo3)[1])

  res[1] = 0.0
  res[1] += ( meas[1][1,idx,1] - the1 )^2
  res[1] += ( meas[1][2,idx,1] - the2 )^2
  res[1] += ( meas[1][3,idx,1] - the3 )^2

  the1 = atan(se2vee(cjTo1)[2], se2vee(cjTo1)[1])
  the2 = atan(se2vee(cjTo2)[2], se2vee(cjTo2)[1])
  the3 = atan(se2vee(cjTo3)[2], se2vee(cjTo3)[1])
  the3b = dobimodal ? atan(se2vee(cjTo3b)[2], se2vee(cjTo3b)[1]) : nothing

  # res[2] = 0.0
  tempres = 0.0
  tempres += ( meas[1][1,idx,2] - the1 )^2
  tempres += ( meas[1][2,idx,2] - the2 )^2

  # FIXME OUTDATED, multihypo= and in-place residual IIF v0.21
  if dobimodal
    if !meas[2][idx]
      # 3a case
      tempres += ( meas[1][3,idx,2] - the3 )^2
    else
      # 3b case
      tempres += ( meas[1][3,idx,2] - the3b )^2
    end
  else
    tempres += ( meas[1][3,idx,2] - the3 )^2
  end
  res[1]+tempres

  return nothing
end

#


text:
# A vector of fixed effects paired with a vector of interactions (empty if not interacted)
const FETerm = Pair{Vector{Symbol},Vector{Symbol}}

# Parse fixed effects from a generic term
function _parsefeterm(@nospecialize(t::AbstractTerm))
    if has_fe(t)
        s = fesymbol(t)
        return [s]=>Symbol[]
    end
end

# Parse fixed effects from an InteractionTerm
function _parsefeterm(@nospecialize(t::InteractionTerm))
    fes = (x for x in t.terms if has_fe(x))
    interactions = (x for x in t.terms if !has_fe(x))
    if !isempty(fes)
        fes = sort!([fesymbol(x) for x in fes])
        feints = sort!([Symbol(x) for x in interactions])
        return fes=>feints
    end
end

getfename(feterm::FETerm) = join(vcat("fe_".*string.(feterm[1]), feterm[2]), "&")

# Count the number of singletons dropped
function drop_singletons!(esample, fe::FixedEffect)
    cache = zeros(Int, fe.n)
    n = 0
    @inbounds for i in eachindex(esample)
        if esample[i]
            cache[fe.refs[i]] += 1
        end
    end
    @inbounds for i in eachindex(esample)
        if esample[i] && cache[fe.refs[i]] <= 1
            esample[i] = false
            n += 1
        end
    end
    return n
end

function Fstat(coef::Vector{Float64}, vcov_mat::AbstractMatrix{Float64}, has_intercept::Bool)
    length(coef) == has_intercept && return NaN
    if has_intercept
        coef = coef[1:end-1]
        vcov_mat = vcov_mat[1:end-1, 1:end-1]
    end
    return (coef' * (vcov_mat \ coef)) / length(coef)
end


text:
factorial_ij(i::Integer, j::Integer) = factorial(big(min(i,j))) / factorial(big(max(i,j)))

z(x::Real, p::Real) = sqrt(2.)*(x + p*im)
z(x::Vector{<:Real}, p::Vector{<:Real}) = z.(x, p')

α(m::Integer, n::Integer) = n - m
α(m::Vector{<:Integer}, n::Vector{<:Integer}) = α.(min.(m, n'), max.(m, n'))

gaussian_function(x::Real, p::Real) = exp(-0.5 * abs2(z(x,p))) / π
gaussian_function(x::Vector{<:Real}, p::Vector{<:Real}) = gaussian_function.(x,p')

neg_one_to_power_of(i::Integer) = (i % 2 == 0) ? 1 : -1

function coefficient_of_wave_function(m::Integer, n::Integer)
    if n ≥ m
        return neg_one_to_power_of(m) * sqrt(factorial_ij(m, n))
    else
        return neg_one_to_power_of(n) * sqrt(factorial_ij(m, n))
    end
end

coefficient_of_wave_function(m::Vector{<:Integer}, n::Vector{<:Integer}) = coefficient_of_wave_function.(m, n')

function z_to_power(m::Integer, n::Integer, x::Real, p::Real)
    if n ≥ m
        return conj(z(x, p'))^(n - m)
    else
        return z(x, p')^(m - n)
    end
end

function z_to_power(m::Vector{<:Integer}, n::Vector{<:Integer}, x::Vector{<:Real}, p::Vector{<:Real})
    x = reshape(x, 1, 1, length(x))
    p = reshape(p, 1, 1, 1, length(p))
    z_to_power.(m, n', x, p)
end

function z_to_power(m::Vector{<:Integer}, n::Vector{<:Integer})
    function z_to_power_xp(x::Vector{<:Real}, p::Vector{<:Real})
        x = reshape(x, 1, 1, length(x))
        p = reshape(p, 1, 1, 1, length(p))
        z_to_power.(m, n', x, p)
    end
    return z_to_power_xp
end

function z_to_power(x::Vector{<:Real}, p::Vector{<:Real})
    x = reshape(x, 1, 1, length(x))
    p = reshape(p, 1, 1, 1, length(p))
    function z_to_power_mn(m::Vector{<:Integer}, n::Vector{<:Integer})
        z_to_power.(m, n', x, p)
    end
    return z_to_power_mn
end


text:
ishorizontal(s::Segment) = s.p1[2]==s.p2[2]
isvertical(s::Segment) = s.p1[1]==s.p2[1]

function draw_segment!(field, s::Segment)
    start, goal = s.p1, s.p2
    field[start[2], start[1]] += 1
    dif = sign.(goal - start)
    current = start
    while current != goal
        current += dif
        field[current[2], current[1]] += 1
    end
end

function day05a(segments; include_diagonals=false)
    field = zeros(Int8, 1000, 1000)
    for s in segments
        if include_diagonals || ishorizontal(s) || isvertical(s)
            draw_segment!(field, s)
        end
    end
    return length(findall(field .> 1))
end

function day05b(segments)
    return day05a(segments, include_diagonals=true)
end

function day05(io)
    segments = parse_segments(io)
    return (day05a(segments), day05b(segments))
end


text:
using HorseML.NeuralNetwork
using HorseML.Preprocessing

@testset "NeuralNetwork" begin
    #Test for utils
    @test_nowarn Dense(10=>5, relu, set_w="He")
    @test_nowarn Dense(10=>5, relu, set_w=rand)
    @test_nowarn Dense(10=>5, relu, high_accuracy=true)
    init_w = Float32.(randn((5, 10)) ./ sqrt(10))
    init_b = Float32.(fill!(similar(init_w, size(init_w, 1)), 0))
    for layer in [Dense, Denseσ]
        @test_nowarn layer(10=>5, relu, initW=init_w, initb=init_b)
        @test_nowarn layer(10=>5, relu, initW=init_w)
    end
    @test_throws ArgumentError Dense(10=>5, relu, set_w=length)
    @test_nowarn Conv((2, 2), 2=>1, relu, set_w="He")
    @test_throws ArgumentError Conv((2, 2), 2=>1, relu, set_w=rand)
    @test_nowarn Conv((2, 2), 2=>1, relu, padding=KeepSize())
    @test_nowarn Conv((2, 2), 2=>1, relu, padding=KeepSize(), high_accuracy=true)
    @test_throws ArgumentError Dropout(-0.2)

    data = [(rand(Float64, 10), rand(Float64)) for i in 1 : 10]
    opt = Descent()
    NN = NetWork(Dense(10=>5, relu), Dense(5=>1, tanh))
    loss(x, y) = LossFunction.mse(NN(x), y)
    train!(NN, loss, data, opt)
    @test_nowarn @epochs 10 train!(NN, loss, data, opt)
    @test_nowarn @simple_epochs 10 train!(NN, loss, data, opt)
    @test typeof(NN[1]) <: Dense
    println(NN)
    @test_nowarn NN(rand(10, 10))

    #Test for activations
    ACTIVATIONS=[
        line, σ, hardσ, hardtanh, relu,
        leakyrelu, relu6, rrelu(0.01, 0.02), elu, gelu, swish, selu,
        celu, softplus, softsign, logσ, logcosh,
        mish, tanhshrink, softshrink, trelu, lisht,
        gaussian, GCU, SQU, NCU, SSU, DSU
    ]
    for f in ACTIVATIONS
        NN = NetWork(Dense(10=>1, f))
        @test_nowarn train!(NN, loss, data, opt)
        @test_nowarn NN(rand(10))
    end

    NN = NetWork(Denseσ(10=>1, prelu()))
    @test_nowarn train!(NN, loss, data, opt)
    @test_nowarn train!(NN, loss, data, opt)

    #Test for optimizers
    for opt in [Descent(), Momentum(), AdaGrad(), Adam()]
        NN = NetWork(Dense(10=>5, relu), Dense(5=>1, tanh))
        @test_nowarn train!(NN, loss, data, opt)
        @test_nowarn NN(rand(10))
    end

    #Test for Conv, Pooling, Flatten, Dropout
    datax = rand(Float32, 5, 10, 10, 2)
    datay = rand(Float32, 5)
    data = databuilder(datax, datay)
    NN = NetWork(MeanPool((2, 2)), Conv((2, 2), 2=>1, relu), MaxPool((2, 2)), Flatten(), Dropout(0.25), Dense(4=>1, tanh))
    opt = Descent()
    @test_nowarn train!(NN, loss, data, opt)
    println(NN)
    @test_nowarn NN(rand(10, 10, 2, 1))
end


text:
import ForwardDiff: Dual, Partials, value
"""
    inverse_flash_update!(storage, eos, c, V)

Update internal matrix of partial derivatives for a converged flash result.
"""
function inverse_flash_update!(storage, eos, c, V)
    x, y = storage.x, storage.y
    J_p, = update_primary_jacobian!(storage, eos, c, storage.forces, x, y, V)
    J_s, = update_secondary_jacobian!(storage, eos, c, storage.forces_secondary, x, y, V)
    invert_sp!(J_s, J_p)
    return J_s
end

function invert_sp!(J_s, J_p)
    F = lu!(J_p)
    ldiv!(F, J_s)
end

function invert_sp!(J_s::MMatrix, J_p::MMatrix)
    tmp = SMatrix(J_p)\SMatrix(J_s)
    J_s .= tmp
end

"""
    set_partials(v, storage, eos, c, index)

Modify a value `v` of ForwardDiff.Dual type to get the correct derivatives.

If the mixture is two-phase, the partial derivatives of phase molar fractions and vapor fraction
with respect to the flash conditions (pressure, temperature and overall mole fractions) are not
trivial. This function is the main gateway for setting these values.

# Notes
Experimental interface, subject to change. You most likely want to use either
[`set_partials_phase_mole_fractions!`](@ref) or [`set_partials_vapor_fraction`](@ref)

In order for this routine to work, the storage must be initialized using `flash_storage`
with the following options enabled:
- `inc_jac = true`
- `diff_externals = true`
and `inverse_flash_update!` must be called after a successful flash. The partial derivatives
with respect to p, T, z is then contained in `storage.buf_inv` with a negative sign. This
function then performs the requisite chain rule operations for the input.
"""
function set_partials(v, storage, eos, c, index)
    M = storage.J_inv
    p, z, T = c.p, c.z, c.T
    buf = storage.buf_inv
    return set_partials(v, M, buf, p, T, z, index)
end

function set_partials(v::Dual{T,V,N}, M, buf, p, temp, z, index) where {T,V,N}
    # Zero out buffer just in case
    @. buf = 0
    # ∂v/∂p
    if isa(p, Dual)
        set_partials_scalar!(buf, p, M, index, 1)
    end
    # ∂v/∂T
    if isa(temp, Dual)
        set_partials_scalar!(buf, temp, M, index, 2)
    end
    # ∂v/∂z_i for all i
    if eltype(z)<:Dual
        set_partials_vector!(buf, z, M, index, 2)
    end
    ∂ = Partials{N, V}(Tuple(buf))
    val = value(v)
    return Dual{T, V, N}(val, ∂)
end

set_partials_scalar!(buf, X::AbstractFloat, M, index, pos) = nothing

function set_partials_scalar!(buf, X::ForwardDiff.Dual{T,V,N}, M, index, pos) where {T, V, N}
    @inbounds @. buf -= M[index, pos]*X.partials
end

set_partials_vector!(buf, z, M, index, offset) = nothing

function set_partials_vector!(buf, z::AbstractVector{ForwardDiff.Dual{T,V,N}}, M, index, offset) where {T, V, N}
    for (ind, zi) = enumerate(z)
        @inbounds ∂z = -M[index, ind+offset]
        for i in eachindex(zi.partials)
            @inbounds buf[i] += ∂z*zi.partials[i]
        end
    end
end

"""
    set_partials_vapor_fraction(V, storage, eos, ∂c)

Set partial derivatives to a vapor mole fraction `ForwardDiff.Dual` instance with correct value, but
missing partial derivatives.
"""
set_partials_vapor_fraction(V, storage, eos, ∂c) = set_partials(V, storage, eos, ∂c, 2*number_of_components(eos) + 1)

"""
set_partials_phase_mole_fractions!(xy, storage, eos, ∂c, [phase_symbol])

Set partial derivatives to phase mole fraction vector with type `ForwardDiff.Dual` with correct values, but
missing partial derivatives.
"""
function set_partials_phase_mole_fractions!(xy, storage, eos, ∂c, phase = :liquid)
    n = length(xy)
    if phase == :liquid
        offset = 0
    else
        offset = n
    end
    @inbounds for i = 1:n
        xy[i] = set_partials(xy[i], storage, eos, ∂c, offset + i)
    end
    return xy
end


text:
module ObservationDims

using Compat: eachslice, hasproperty
using Distributions
using NamedDims
using Tables

export obs_arrangement, organise_obs
export SingleObs, IteratorOfObs, ArraySlicesOfObs
export MatrixRowsOfObs, MatrixColsOfObs

"""
    type ObsArrangement

Defines the orientation of data that is expected by a function
"""
abstract type ObsArrangement end

"""
    SingleObs <: ObsArrangement

The data consists of a single observation regardless of dimension.
"""
struct SingleObs <: ObsArrangement end

"""
    IteratorOfObs <: ObsArrangement

The data consists of an iterator over multiple observations.
"""
struct IteratorOfObs <: ObsArrangement end

"""
    ArraySlicesOfObs{D} <: ObsArrangement

The data consists of a multi-dimensional array where the observations are along dimension `D`.
"""
struct ArraySlicesOfObs{D} <: ObsArrangement end

"""
    MatrixRowsOfObs <: ObsArrangement

A special case of `ArraySlicesOfObs` where the observations are along the first dimension.
"""
const MatrixRowsOfObs = ArraySlicesOfObs{1}

"""
    MatrixColsOfObs <: ObsArrangement
A special case of `ArraySlicesOfObs` where the observations are along the second dimension.
"""
const MatrixColsOfObs = ArraySlicesOfObs{2}

"""
    obs_arrangement(f::Type{<:Function}) -> ObsArrangement

Specify the observation arrangement trait of a function `f`.
"""
function obs_arrangement end

"""
    _TableHolder

A internal wrapper type for dispatch purposes.
Since tables do not have a type, wrap them in this so we can dispatch on them.
"""
struct _TableHolder{T}
    data::T
end
"""
    organise_obs(f, data; obsdim=_default_obsdim(data))
    organise_obs(::ObsArrangement, data; obsdim=_default_obsdim(data))

Organise the `data` according to the `ObsArrangement` expected by `f`.

# Arguments
- `f`: the function or method needing the data in a certain orientation
- `data`: the data to transform
- `obsdim`: the dimension of the observations, resorts to `_default_obsdim` when not specified.
"""
function organise_obs(f, data; obsdim=_default_obsdim(data))
    return organise_obs(obs_arrangement(f), data; obsdim=obsdim)
end

# Specify arrangement based on type of data:

# Scalars have no "orientation" so no rearrangement required
for T in (Sampleable, Number, Symbol)
    @eval organise_obs(::SingleObs, data::$T; obsdim=nothing) = data
    @eval organise_obs(::IteratorOfObs, data::$T; obsdim=nothing) = data
    @eval organise_obs(::ArraySlicesOfObs, data::$T; obsdim=nothing) = data
end

## Vectors: obsdim is optional, we may or may not need it.
for T in (Any, AbstractVector)
    # Iterator -> IteratorOfObs
    @eval function organise_obs(::IteratorOfObs, obs_iter::$T; obsdim=nothing)
        if Tables.istable(obs_iter)
            return organise_obs(IteratorOfObs(), _TableHolder(obs_iter); obsdim=obsdim)
        else
            return obs_iter
        end
    end

    # Iterator -> ArraySlicesOfObs
    @eval function organise_obs(::ArraySlicesOfObs{D}, obs_iter::$T; obsdim=nothing) where D
        if Tables.istable(obs_iter)
            return organise_obs(ArraySlicesOfObs{D}(), _TableHolder(obs_iter); obsdim=obsdim)
        end
        # we assume all obs have same number of dimensions else nothing makes sense
        ndims_per_obs = ndims(first(obs_iter))

        # If collection of scalars then return the collection
        if ndims_per_obs == 0
            return collect(obs_iter)
        end

        # This should just be a mapreduce but that is slow
        # see https://github.com/JuliaLang/julia/issues/31137
        shaped_obs = Base.Generator(obs_iter) do obs
            new_shape = ntuple(ndims_per_obs + 1) do ii
                if ii < D
                    size(obs, ii)
                elseif ii > D
                    size(obs, ii-1)
                else  # ii = D
                    1
                end
            end
            # add singleton dim that we wil concatenate on
            Base.ReshapedArray(obs, new_shape, ())
        end
        return cat(shaped_obs...; dims=D)
    end
end

# Specify arrangement based on desired ObsArrangement:

# Any -> SingleObs: never any need to rearrage
organise_obs(::SingleObs, data; obsdim=nothing) = data

# Tables support
function organise_obs(arrangement::IteratorOfObs, holder::_TableHolder; obsdim=1)
    _warn_about_table_obsdim(obsdim)
    data = holder.data
    return (collect(obs) for obs in Tables.rows(data))
end

function organise_obs(arrangement::ArraySlicesOfObs, holder::_TableHolder; obsdim=1)
    _warn_about_table_obsdim(obsdim)
    data = Tables.matrix(holder.data)
    return organise_obs(arrangement, data; obsdim=1)
end

function _warn_about_table_obsdim(obsdim)
    if obsdim !== 1 && obsdim !== nothing
        @warn "Arraying a Table, obsdim not equal to 1 ignored" obsdim
    end
end

# Array -> IteratorOfObs or ArraySlicesOfObs: depends on obsdim
# Resorts to default obsdim which redispatches to the 3 arg form below.
for A in (IteratorOfObs, ArraySlicesOfObs)

    @eval function organise_obs(arrangement::$A, data::AbstractArray; obsdim=_default_obsdim(data))
        if obsdim isa Symbol && hasproperty(data, :axes)  # support AxisArrays via structural ducktyping
            obsdim = findfirst(x -> in(obsdim, typeof(x).parameters), data.axes)
        end
        return organise_obs(arrangement, data, obsdim)
    end

    @eval function organise_obs(arrangement::$A, data::NamedDimsArray; obsdim=_default_obsdim(data))
        obsdim = (obsdim isa Symbol) ? NamedDims.dim(data, obsdim) : obsdim
        return organise_obs(arrangement, data, obsdim)
    end
end

# 3 arg forms rearrange (non 1D) arrays according to the obsdim

# Slice up the array to get an iterator of observations
function organise_obs(::IteratorOfObs, data::AbstractArray, obsdim::Integer)
    return eachslice(data, dims=obsdim)
end

# If default obsdim is nothing we should gracefully dispatch on the default
function organise_obs(arrangement::IteratorOfObs, data, ::Nothing)
    return organise_obs(arrangement, data, _default_obsdim(data))
end

function organise_obs(arrangement::ArraySlicesOfObs, data, ::Nothing)
    return organise_obs(arrangement, data, _default_obsdim(data))
end

# Permute the array so the observations are arranged correctly
function organise_obs(
    ::ArraySlicesOfObs{D}, data::AbstractArray{<:Any, N}, obsdim::Integer
) where {D, N}

    if obsdim == D
        return data
    else
        # Swap around the obsdim with the dimension we want to assign it to
        perm = ntuple(N) do ii
            if ii == D
                return obsdim
            elseif ii == obsdim
                return D
            else
                return ii
            end
        end
        return PermutedDimsArray(data, perm)
    end
end

# Assign rows as observations by default
_default_obsdim(x) = 1
function _default_obsdim(x::NamedDimsArray{L}) where L
    obsnames = (:obs, :observations, :samples)

    # These obsnames specify the order of preference when returning your observation dimension
    # e.g. if :obs and :samples both exist in your NamedDimsArray then :obs is always returned
    used = findfirst(in(L), obsnames)
    if used === nothing
        throw(DimensionMismatch(string(
            "No observation dimension found. Provide one of the valid dimension names = $L."
        )))
    end
    return @inbounds obsnames[used]
end

end # module


text:
#===============================================================================

    amp_analysis.jl

    amplification analysis

    issue:
    the code assumes only 1 step/ramp because the current data format
    does not allow us to break the fluorescence data down by step_id/ramp_id

===============================================================================#

import DataStructures.OrderedDict
import StaticArrays: SVector, SMatrix
import Memento: debug, warn
using Ipopt

#===============================================================================
    function definitions >>
===============================================================================#

"Analyse amplification data via calls to `calibrate` and `get_fit_results`."
function amp_analysis(i ::AmpInput) # ; asrp ::AmpStepRampProperties)
    debug(logger, "at amp_analysis()")
    #
    ## deconvolute and normalize
    calibration_results =
        calibrate(i, i.calibration_data, i.calibration_args, i.raw, array)
    #
    ## initialize output
    o = AmpOutput(
        Val{i.amp_output},
        i,
        calibration_results...,
        # cq_method,
        default_ct_fluos(i))
    #
    ## fit amplification models and report results
    if i.num_cycles <= 2
		num_cycles = i.num_cycles
        warn(logger, "number of cycles $num_cycles <= 2: baseline subtraction " *
            "and Cq calculation will not be performed")
    else ## num_cycles > 2
        baseline_cyc_bounds = check_bl_cyc_bounds(i, DEFAULT_AMP_BL_CYC_BOUNDS)
        try
            set_ct_fluos!(o, i, baseline_cyc_bounds)
        catch e
            warn(logger, "set_ct_fluos catch error: " * sprint(showerror, e))
        end
        try
            set_output_fields!(o, i, get_fit_results(o, i, baseline_cyc_bounds))
        catch e
            warn(logger, "get_fit_results catch error: " * sprint(showerror, e))
        end
        set_qt_fluos!(o, i)
        set_report_cq!(o, i)
    end ## if
    return o
end ## amp_analysis()

#==============================================================================#

"Set the amplification output field `ct_fluos`."
function set_ct_fluos!(
    o                       ::AmpOutput,
    i                       ::AmpInput,
    baseline_cyc_bounds     ::AbstractArray,
)
    debug(logger, "at calc_ct_fluos()")
    (length(o.ct_fluos) > 0) && return nothing
    o.ct_fluos = default_ct_fluos(i)
    (i.cq_method != :ct)     && return nothing
    (i.amp_model != :SFC)    && return nothing
    ## else
    o.ct_fluos =
        map(1:i.num_channels) do channel_i
            fits =
                map(1:i.num_wells) do well_i
                    fluos = o.rbbs_ary3[:, well_i, channel_i]
                    amp_fit_model(
                        Val{SFCModel},
                        AmpCqFluoModelResults,
                        i,
                        fluos,
                        bl_cyc_bounds[well_i, channel_i],
                        DEFAULT_AMP_CT_FLUO_METHOD, ## cq_method
                        NaN_T) ## i.ct_fluo
                end ## do well_i
            fits |>
                their(:quant_status) |>
                find_idc_useful |>
                curry(getindex)(fits) |>
                their(:cq_fluo) |>
                median
        end ## do channel_i
    return nothing ## side effects only
end ## calc_ct_fluos()

## called by calc_ct_fluos() >>

@inline function find_idc_useful(postbl_stata ::AbstractVector)
    idc_useful = find(postbl_stata .== :Optimal)
    (length(idc_useful) > 0) && return idc_useful
    idc_useful = find(postbl_stata .== :UserLimit)
    (length(idc_useful) > 0) && return idc_useful
    return eachindex(postbl_stata)
end ## find_idc_useful()

default_ct_fluos(i ::AmpInput) =
    SVector{i.num_channels, Float_T}(fill(NaN_T, i.num_channels))

#==============================================================================#

"Fit amplification model to data for each well and channel."
function get_fit_results(
    o                       ::AmpOutput,
    i                       ::AmpInput,
    bl_cyc_bounds           ::AbstractArray,
)
    function do_model_fit(wi ::Int_T, ci ::Int_T)
        debug(logger, "at do_model_fit($wi, $ci)")
        if isa(solver, Ipopt.IpoptSolver) && length(prefix) > 0
            ipopt_file = string(join([prefix, ci, wi], '_')) * ".txt"
            push!(solver.options, (:output_file, ipopt_file))
        end
        fluos = o.rbbs_ary3[:, wi, ci]
        amp_fit_model(
            Val{i.amp_model},
            i.amp_model_results,
            i,
            fluos,
            bl_cyc_bounds[wi, ci],
            i.cq_method,
            o.ct_fluos[ci])
    end ## fit model()

    ## << end of function definition nested within set_fit_results!()

    debug(logger, "at set_fit_results!()")
    solver = i.solver
    prefix = i.ipopt_print_prefix
    fit_results =
        SMatrix{i.num_wells, i.num_channels, i.amp_model_results}([
            do_model_fit(wi, ci)
            for wi in 1:i.num_wells, ci in 1:i.num_channels])
end ## set_fit_results!()

#==============================================================================#

"Format the results of the amplification analyses."
@inline function set_output_fields!(
    o                       ::AmpOutput,
    i                       ::AmpInput,
    results                 ::AbstractArray,
)
    debug(logger, "at set_output_fields!()")
    foreach(fieldnames(first(results))) do fieldname
        output_field = getfield(o, fieldname)
        T = output_field |> eltype
        vector_output_field = ndims(output_field) == 3
        if vector_output_field
            setfield!(o, fieldname,
                results |>
                moose(bless(Vector{T}) ∘ field(fieldname), hcat) |>
                morph(:, i.num_wells, i.num_channels))
        else
            setfield!(o, fieldname,
                results |> mold(bless(T) ∘ field(fieldname)) |>
                bless(SMatrix{i.num_wells, i.num_channels, T}))
        end ## if
    end ## next fieldname
    return nothing ## side effects only
end ## set_output_fields!()

#==============================================================================#

"Calculate the amplification analysis fields `qt_fluos` and `max_qt_fluo`, when
the output format is `long`."
function set_qt_fluos!(
    o                       ::AmpLongOutput,
    i                       ::AmpInput,
)
    debug(logger, "at set_qt_fluos!()")
    o.qt_fluos =
        [   quantile(o.blsub_fluos[:, well_i, channel_i], i.qt_prob)
            for well_i in 1:i.num_wells, channel_i in 1:i.num_channels ]
    o.max_qt_fluo = maximum(o.qt_fluos)
    return nothing ## side effects only
end ## set_qt_fluos!()

"Do nothing when the output format is `short`."
set_qt_fluos!(
    o                       ::AmpShortOutput,
    i                       ::AmpInput,
) = nothing

#==============================================================================#

"Call `report_cq!` for each well and channel, when the output format is `long`."
function set_report_cq!(
    o                       ::AmpLongOutput,
    i                       ::AmpInput,
)
    debug(logger, "at set_report_cq!()")
    for well_i in 1:i.num_wells, channel_i in 1:i.num_channels
        report_cq!(i, o, well_i, channel_i)
    end
    return nothing ## side effects only
end ## set_report_cq!()

"Do nothing when the output format is `short`."
function set_report_cq!(
    o                       ::AmpShortOutput,
    i                       ::AmpInput,
)
    debug(logger, "at set_report_cq!() SHORT FORMAT")
    temp_cq=fill(NaN_T,size(o.cq)[1], size(o.cq)[2])

    max_dr1_lb, max_dr2_lb, max_bsf_lb = i.max_dr1_lb, i.max_dr2_lb, i.max_bsf_lb
    qt_fluos_ =
        [   quantile(o.blsub_fluos_flb[:, well_i, channel_i], i.qt_prob)
            for well_i in 1: i.num_wells, channel_i in 1: i.num_channels ]
    max_qt_fluo_ = maximum(qt_fluos_)
    for well_i in 1: i.num_wells, channel_i in 1: i.num_channels
        b_ = o.coefs[1, well_i, channel_i]
        postbl_status = o.quant_status[well_i, channel_i]
        max_dr1= maximum(o.dr1_pred[:, well_i, channel_i])
        max_dr2= maximum(o.dr2_pred[:, well_i, channel_i])
        max_bsf = maximum(o.blsub_fluos_flb[:, well_i, channel_i])
		max_bsf_ = maximum(o.blsub_fluos[(end - Int(floor((size(o.blsub_fluos)[1]) / 4))): end, well_i, channel_i])
        (scaled_max_dr1, scaled_max_dr2, scaled_max_bsf) =
        [max_dr1, max_dr2, max_bsf] ./ max_qt_fluo_

        min_reliable_cyc_ = length(i.baseline_cyc_bounds) == 2 ? i.baseline_cyc_bounds[1]: i.min_reliable_cyc

        if (max_bsf < i.max_bsf_lb || scaled_max_bsf < i.scaled_max_bsf_lb ||
            max_dr2 < max_dr2_lb || max_dr1 < max_dr1_lb  ||
            scaled_max_dr2 < i.scaled_max_dr2_lb ||
            scaled_max_dr1 < i.scaled_max_dr1_lb || o.cq[well_i, channel_i] > i.num_cycles ||
            o.cq[well_i, channel_i] < min_reliable_cyc_ || postbl_status == :Error ||
            b_ > 0)

            why_cq_failed = "Cq for W:$well_i, C:$channel_i failed due to: "
            why_cq_failed *= max_bsf < i.max_bsf_lb ? "max_bsf($max_bsf) < i.max_bsf_lb($(i.max_bsf_lb)) |" : ""
            why_cq_failed *= scaled_max_bsf < i.scaled_max_bsf_lb ? "scaled_max_bsf($scaled_max_bsf) < i.scaled_max_bsf_lb($(i.scaled_max_bsf_lb)) |" : ""
            why_cq_failed *= max_dr2 < max_dr2_lb ? "max_dr2($max_dr2) < max_dr2_lb($max_dr2_lb) |" : ""
            why_cq_failed *= max_dr1 < max_dr1_lb ? "max_dr1 < max_dr1_lb : " * "$max_dr1 < $max_dr1_lb |" : ""
            why_cq_failed *= scaled_max_dr2 < i.scaled_max_dr2_lb ? "scaled_max_dr2($scaled_max_dr2) < i.scaled_max_dr2_lb($(i.scaled_max_dr2_lb)) |" : ""
            why_cq_failed *= scaled_max_dr1 < i.scaled_max_dr1_lb ? "scaled_max_dr1($scaled_max_dr1) < i.scaled_max_dr1_lb($(i.scaled_max_dr1_lb)) |" : ""
            why_cq_failed *= o.cq[well_i, channel_i] > i.num_cycles ? "o.cq[well_i, channel_i]($(o.cq[well_i, channel_i])) > i.num_cycles($(i.num_cycles)) |" : ""
            why_cq_failed *= o.cq[well_i, channel_i] < min_reliable_cyc_ ? "o.cq[well_i, channel_i]($(o.cq[well_i, channel_i])) < min_reliable_cyc_($(min_reliable_cyc_)) |" : ""
            why_cq_failed *= postbl_status == :Error ? "postbl_status == :Error |" : ""
            why_cq_failed *= b_ > 0 ? "b_($b_) > 0 |" : ""
            debug(logger, why_cq_failed)

            o.blsub_fluos[:,well_i, channel_i]=o.blsub_fluos_flb[:,well_i, channel_i]
            o.dr1_pred[:,well_i, channel_i]=o.dr1_pred1[:,well_i, channel_i]
            o.dr2_pred[:,well_i, channel_i]=o.dr2_pred1[:,well_i, channel_i]

        else
            temp_cq[well_i, channel_i] = o.cq[well_i, channel_i]
            debug(logger, "Cq for W:$well_i, C:$channel_i >> $(o.cq[well_i, channel_i])")
        end
    end
    o.cq=convert(typeof(o.cq), temp_cq)
    return nothing
end

"Report amplification output fields relating to the calculation of `cq`."
function report_cq!(
    o                       ::AmpLongOutput,
    i                       ::AmpInput,
    well_i                  ::Int_T,
    channel_i               ::Int_T,
)
    if i.before_128x
        max_dr1_lb, max_dr2_lb, max_bsf_lb = [i.max_dr1_lb, i.max_dr2_lb, i.max_bsf_lb] ./ 128
    else
        max_dr1_lb, max_dr2_lb, max_bsf_lb = i.max_dr1_lb, i.max_dr2_lb, i.max_bsf_lb
    end
    #
    num_cycles = size(o.raw_data, 1)
    (postbl_status, cq_raw, max_dr1, max_dr2) =
        map([ :postbl_status, :cq_raw, :max_dr1, :max_dr2 ]) do fieldname
            fieldname -> getfield(o, fieldname)[well_i, channel_i]
        end
    max_bsf = maximum(o.blsub_fluos[:, well_i, channel_i])
    b_ = o.coefs[1, well_i, channel_i]
    (scaled_max_dr1, scaled_max_dr2, scaled_max_bsf) =
        [max_dr1, max_dr2, max_bsf] ./ o.max_qt_fluo
    why_NaN =
        if postbl_status == :Error
            "postbl_status == :Error"
        elseif b_ > 0
            "b > 0"
        elseif o.cq_method == ct && o.cq_raw == AMP_CT_VAL_DOMAINERROR
            "DomainError when calculating Ct"
        elseif o.cq_raw <= 0.1 || o.cq_raw >= num_cycles
            "cq_raw <= 0.1 || cq_raw >= num_cycles"
        elseif max_dr1 < max_dr1_lb
            "max_dr1 $max_dr1 < max_dr1_lb $max_dr1_lb"
        elseif max_dr2 < max_dr2_lb
            "max_dr2 $max_dr2 < max_dr2_lb $max_dr2_lb"
        elseif max_bsf < max_bsf_lb
            "max_bsf $max_bsf < max_bsf_lb $max_bsf_lb"
        elseif scaled_max_dr1 < i.scaled_max_dr1_lb
            "scaled_max_dr1 $scaled_max_dr1 < scaled_max_dr1_lb $(i.scaled_max_dr1_lb)"
        elseif scaled_max_dr2 < i.scaled_max_dr2_lb
            "scaled_max_dr2 $scaled_max_dr2 < scaled_max_dr2_lb $(i.scaled_max_dr2_lb)"
        elseif scaled_max_bsf < i.scaled_max_bsf_lb
            "scaled_max_bsf $scaled_max_bsf < scaled_max_bsf_lb $(i.scaled_max_bsf_lb)"
        else
            ""
        end ## why_NaN
    (why_NaN != "") && (o.cq[well_i, channel_i] = NaN_T)
    #
    for tup in (
        (:max_bsf,        max_bsf),
        (:scaled_max_dr1, scaled_max_dr1),
        (:scaled_max_dr2, scaled_max_dr2),
        (:scaled_max_bsf, scaled_max_bsf),
        (:why_NaN,        why_NaN))
        getfield(o, tup[1])[well_i, channel_i] = tup[2]
    end
    return nothing ## side effects only
end ## report_cq!

#===============================================================================
    helper functions >>
===============================================================================#

amp_init(i ::AmpInput, x) =
    SMatrix{i.num_wells, i.num_channels, typeof(x)}(
        fill(x, i.num_wells, i.num_channels))

function check_bl_cyc_bounds(
    i               ::AmpInput,
    bl_cyc_bounds   ::Union{Vector{I},Array{Vector{I},2}} where {I <: Int_T},
)
    debug(logger, "at check_bl_cyc_bounds()")
    (i.num_cycles <= 2) && return bl_cyc_bounds
    size_bcb = size(bl_cyc_bounds)
    if size_bcb == (0,) || size_bcb == (2,)
        return amp_init(i, bl_cyc_bounds)
    elseif size_bcb == (i.num_wells, i.num_channels) &&
        eltype(bl_cyc_bounds) <: AbstractVector ## final format of `baseline_cyc_bounds`
            return bl_cyc_bounds
    end
    throw(ArgumentError("`baseline_cyc_bounds` is not in the right format"))
end ## check_bl_cyc_bounds()

## baseline estimation parameters
# kwargs_bl(i ::AmpInput) =
#     Dict{Symbol,Any}(
#         :bl_method          => i.bl_method,
#         :bl_fallback_func   => i.bl_fallback_func,
#         :min_reliable_cyc   => i.min_reliable_cyc,
#     )

## quantitation parameters
# kwargs_quant(i ::AmpInput) =
#     Dict{Symbol,Any}(
#         :cq_method          => i.cq_method,
#         :denser_factor      => i.denser_factor,
#     )

## arguments for process_ad()
# kwargs_ad(i ::AmpInput) =
#     Dict{Symbol,Any}(
#         :ctrl_well_dict     => i.ctrl_well_dict,
#         # :cluster_method     => i.cluster_method,
#         # :norm_l             => i.norm_l,
#         # :encgr              => i.encgr,
#         # :categ_well_vec     => i.categ_well_vec
#     )

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using Base.Test
using RMB

# println("hi")

tol=1e-10

a=rand(2,5)
b=rand(5,2)
@test norm(I-RMB.explicitinv(a,b)*[a zeros(2,2);eye(5) b]) < tol


text:
using GeoStats
using Plots; gr(size=(600,400))
using VisualRegressionTests
using Test, Random, Pkg

# list of maintainers
maintainers = ["juliohm"]

# environment settings
istravis = "TRAVIS" ∈ keys(ENV)
ismaintainer = "USER" ∈ keys(ENV) && ENV["USER"] ∈ maintainers
datadir = joinpath(@__DIR__,"data")

# load data sets
data1D = readgeotable(joinpath(datadir,"data1D.tsv"), delim='\t', coordnames=[:x])
data2D = readgeotable(joinpath(datadir,"data2D.tsv"), delim='\t', coordnames=[:x,:y])

# list of tests
testfiles = [
  "solvers.jl",
  "comparisons.jl"
]

@testset "GeoStats.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
module PlotlyLight

using Random
using JSON3
using EasyConfig
using Cobweb

export Plot, Config

function __init__()
    global plotlyjs = joinpath(@__DIR__, "..", "deps", "plotly-2.8.3.min.js")
    !(isfile(plotlyjs)) && @warn "Cannot find plotly.js. PlotlyLight should be built again."
end

#-----------------------------------------------------------------------------# defaults
module Defaults
using EasyConfig: Config

src             = Ref(:cdn)
class           = Ref("")
style           = Ref("height: 100%;")
parent_class    = Ref("")
parent_style    = Ref("height: 100vh;")
config          = Ref(Config(displaylogo=false, responsive=true))
layout          = Ref(Config())

function reset!()
    src[]           = :cdn
    class[]         = ""
    style[]         = "height: 100%;"
    parent_class[]  = ""
    parent_style[]  = "height: 100vh;"
    config[]        = Config(displaylogo=false, responsive=true)
    layout[]        = Config()
end
end

#-----------------------------------------------------------------------------# src
src_opts = [:cdn, :local, :standalone, :none]
"""
    src!(x::Symbol) # `x` must be one of: $src_opts

- `:cdn` → Use PlotlyJS CDN.
- `:local` → Use local artifact.
- `:standalone` → Write JS into the HTML file directly (can be shared and viewed offline).
- `:none` → For when inserting into a page with Plotly.js already included.
"""
src!(x::Symbol) = (x in src_opts || error("src must be one of: $src_opts"); Defaults.src[] = x)

#-----------------------------------------------------------------------------# Plot
"""
    Plot(data, layout, config)

- A Plotly.js plot with components `data`, `layout`, and `config`.
    - `data = Config()`: A `Config` (single trace) or `Vector{Config}` (multiple traces).
    - `layout = Config()`.
    - `config = Config(displaylogo=false, responsive=true)`.
- Each of the three components are converted to JSON via `JSON3.write`.
- See the Plotly Javascript docs here: https://plotly.com/javascript/.

### Example

    p = Plot(Config(x=1:10, y=randn(10)))
    p.layout.title.text = "My Title!"
    p
"""
Base.@kwdef mutable struct Plot
    data::Vector{Config}    = Config[]
    layout::Config          = Defaults.layout[]
    config::Config          = Defaults.config[]
    id::String              = randstring(10)            # id of graphDiv
    js::Cobweb.Javascript   = Cobweb.Javascript("console.log('plot created!')")
end
function Plot(traces, layout=Defaults.layout[], config=Defaults.config[]; kw...)
    data = traces isa Config ? [traces] : traces
    Plot(; kw..., data, layout, config)
end

#-----------------------------------------------------------------------------# Display
Base.display(::Cobweb.CobwebDisplay, o::Plot) = display(Cobweb.CobwebDisplay(), Cobweb.Page(o))

function Base.show(io::IO, M::MIME"text/html", o::Plot)
    class, style = Defaults.class, Defaults.style
    parent_class, parent_style = Defaults.parent_class, Defaults.parent_style
    parent_style = if get(io, :is_pluto, false)
        s = replace(parent_style[], r"height.*;" => "")
        "height: 400px;" * s
    else
        parent_style[]
    end
    src = Defaults.src[]
    src in [:cdn, :standalone, :none, :local] || error("`src` must be :cdn, :standalone, :none, or :local")
    println(io, "
")
    println(io, "

")
    println(io, "
")

    if src === :cdn
        println(io, "")
    elseif src === :standalone
        print(io, "")
    elseif src === :local
        println(io, "")
    else
        # :none
    end

    println(io, "\n")
end

end # module


text:
@testset "learnedheuristic.jl" begin

    @testset "LearnedHeuristic" begin

        learnedheuristic = SeaPearl.LearnedHeuristic(agent)

        @test learnedheuristic.agent == agent
        @test isnothing(learnedheuristic.fitted_problem)
        @test isnothing(learnedheuristic.fitted_strategy)
        @test isnothing(learnedheuristic.action_space)
        @test isnothing(learnedheuristic.current_state)
        @test isnothing(learnedheuristic.reward)
        @test isnothing(learnedheuristic.search_metrics)

        learnedheuristic.fitted_problem = SeaPearl.LegacyGraphColoringGenerator
        learnedheuristic.fitted_strategy = SeaPearl.DFSearch

        @test learnedheuristic.fitted_problem == SeaPearl.LegacyGraphColoringGenerator
        @test learnedheuristic.fitted_strategy == SeaPearl.DFSearch

    end

    @testset "LearnedHeuristic in action" begin

    approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    target_approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    gnnlayers = 1
    approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
            [approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, 4),
    )
    target_approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
            [target_approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, 4),
    )

    agent = RL.Agent(
        policy = RL.QBasedPolicy(
            learner = RL.DQNLearner(
                approximator = RL.NeuralNetworkApproximator(
                    model = approximator_model,
                    optimizer = ADAM(0.001f0)
                ),
                target_approximator = RL.NeuralNetworkApproximator(
                    model = target_approximator_model,
                    optimizer = ADAM(0.001f0)
                ),
                loss_func = Flux.Losses.huber_loss,
                stack_size = nothing,
                γ = 0.99f0,
                batch_size = 32,
                update_horizon = 1,
                min_replay_history = 1,
                update_freq = 1,
                target_update_freq = 100,
            ),
            explorer = RL.EpsilonGreedyExplorer(
                ϵ_stable = 0.01,
                kind = :exp,
                ϵ_init = 1.0,
                warmup_steps = 0,
                decay_steps = 500,
                step = 1,
                is_break_tie = false,
                #is_training = true,
                rng = MersenneTwister(33)
            )
        ),
        trajectory = RL.CircularArraySLARTTrajectory(
            capacity = 500,
            state = SeaPearl.DefaultTrajectoryState[] => (),
            action = Int => (),
            legal_actions_mask = Vector{Bool} => (4, ),
        )
    )

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(1, 2, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 2, "x2", trailer)
        x3 = SeaPearl.IntVar(2, 3, "x3", trailer)
        x4 = SeaPearl.IntVar(1, 4, "x4", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)
        SeaPearl.addVariable!(model, x4)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x3, x4, trailer))

        lh = SeaPearl.LearnedHeuristic(agent)
        SeaPearl.update_with_cpmodel!(lh, model)
        false_x = first(values(model.variables))
        env = SeaPearl.get_observation!(lh, model, false_x)

        Flux.reset!(lh.agent)
        lh.agent(RL.PRE_EPISODE_STAGE, env)
        _, _ = SeaPearl.fixPoint!(model)
        env = SeaPearl.get_observation!(lh, model, x1)

        v1 = lh.agent(env)
        lh.agent(RL.PRE_ACT_STAGE, env, v1)
        SeaPearl.assign!(x1, v1)
        _, _ = SeaPearl.fixPoint!(model, SeaPearl.getOnDomainChange(x1))
        env = SeaPearl.get_observation!(lh, model, x2)
        lh.agent(RL.POST_ACT_STAGE, env)
        v2 = lh.agent(env)
        lh.agent(RL.PRE_ACT_STAGE, env, v2)
        SeaPearl.assign!(x2, v2)
        _, _ = SeaPearl.fixPoint!(model, SeaPearl.getOnDomainChange(x2))
        env = SeaPearl.get_observation!(lh, model, x3)
        lh.agent(RL.POST_ACT_STAGE, env)

        v3 = lh.agent(env)
        lh.agent(RL.PRE_ACT_STAGE, env, v3)
        SeaPearl.assign!(x3, v3)
        _, _ = SeaPearl.fixPoint!(model, SeaPearl.getOnDomainChange(x3))
        env = SeaPearl.get_observation!(lh, model, x4)
        lh.agent(RL.POST_ACT_STAGE, env)

        v4 = lh.agent(env)
        lh.agent(RL.PRE_ACT_STAGE, env, v4)
        SeaPearl.assign!(x4, v4)
        _, _ = SeaPearl.fixPoint!(model, SeaPearl.getOnDomainChange(x4))
        env = SeaPearl.get_observation!(lh, model, x4)
        lh.agent(RL.POST_ACT_STAGE, env)

        @test length(lh.agent.trajectory[:state]) == 4
        @test size(lh.agent.trajectory[:legal_actions_mask]) == (4,4)
        @test length(lh.agent.trajectory[:action]) == 4
        @test length(lh.agent.trajectory[:reward]) == 4
        @test length(lh.agent.trajectory[:terminal]) == 4

    end

    @testset "advanced test on learnedheuristic.jl" begin

    approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    target_approximator_GNN = SeaPearl.GraphConv(64 => 64, Flux.leakyrelu)
    gnnlayers = 1
    approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
            [approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, 4),
    )
    target_approximator_model = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            SeaPearl.GraphConv(3 => 64, Flux.leakyrelu),
            [target_approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu),
            Flux.Dense(32, 32, Flux.leakyrelu),
            Flux.Dense(32, 16, Flux.leakyrelu),
        ),
        outputChain = Flux.Dense(16, 4),
    )

    agent = RL.Agent(
        policy = RL.QBasedPolicy(
            learner = RL.DQNLearner(
                approximator = RL.NeuralNetworkApproximator(
                    model = approximator_model,
                    optimizer = ADAM(0.001f0)
                ),
                target_approximator = RL.NeuralNetworkApproximator(
                    model = target_approximator_model,
                    optimizer = ADAM(0.001f0)
                ),
                loss_func = Flux.Losses.huber_loss,
                stack_size = nothing,
                γ = 0.99f0,
                batch_size = 32,
                update_horizon = 1,
                min_replay_history = 1,
                update_freq = 1,
                target_update_freq = 100,
            ),
            explorer = RL.EpsilonGreedyExplorer(
                ϵ_stable = 0.01,
                kind = :exp,
                ϵ_init = 1.0,
                warmup_steps = 0,
                decay_steps = 500,
                step = 1,
                is_break_tie = false,
                #is_training = true,
                rng = MersenneTwister(33)
            )
        ),
        trajectory = RL.CircularArraySLARTTrajectory(
            capacity = 500,
            state = SeaPearl.DefaultTrajectoryState[] => (),
            action = Int => (),
            legal_actions_mask = Vector{Bool} => (4, ),
        )
    )

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(1, 2, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 2, "x2", trailer)
        x3 = SeaPearl.IntVar(2, 3, "x3", trailer)
        x4 = SeaPearl.IntVar(1, 4, "x4", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)
        SeaPearl.addVariable!(model, x4)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x3, x4, trailer))

        variableHeuristic = SeaPearl.MinDomainVariableSelection{false}()
        x = variableHeuristic(model)

        lh = SeaPearl.LearnedHeuristic(agent)
        SeaPearl.update_with_cpmodel!(lh, model)

        lh.firstActionTaken = true
        lh(SeaPearl.DecisionPhase, model, x)
        @test !isnothing(model.statistics.lastVar)
        @test lh.agent.trajectory[:reward][end] == -0.025f0

        lh(SeaPearl.EndingPhase, model, :Optimal)
        @test lh.agent.trajectory[:reward][end] == 0.975f0    #last DecisionReward (-0.025f0) + EndingReward (+1)

    end
end


text:
export lanford, doubling, tupling

# Lanford map
for (fn,dfn,n) in ((:lanford_v0,:lanford_dv0,25),(:lanford_v1,:lanford_dv1,17))
  @eval ($fn)(x) = (5-sqrt($n-8x))/2
  @eval ($dfn)(x) = 2/sqrt($n-8x)
end
lanford_brk(T=Float64) = lanford_v0(one(T))

"""
    lanford(T=Float64)

Return the Lanford map, with type encoding T.

See also: [`tupling`](@ref), [`doubling`](@ref)
"""
lanford(T=Float64) = MarkovMap([lanford_v0,lanford_v1],[0..lanford_brk(T),lanford_brk(T)..1];dir=Reverse,diff=[lanford_dv0,lanford_dv1])

# tupling map
"""
    tupling(k::Int, d = Interval(0,1.))

Returns the full-branch interval map on domain d with k equally-sized branches.

See also: [`doubling`](@ref), [`lanford`](@ref)
"""
tupling(k::Int,d=Interval(0,1.)) = modulomap(x->k*(x-leftendpoint(d))+leftendpoint(d),d)
tupling(k::Int,d::PeriodicSegment) = RevCircleMap(x->(x-leftendpoint(d))/k+leftendpoint(d),d,d,x->one(x)/k)

"""
    doubling(d = Interval(0,1.))

Returns the full-branch interval map on domain d with 2 equally-sized branches.

See also: [`tupling`](@ref), [`lanford`](@ref)
"""
doubling(d=Interval(0,1.)) = tupling(2,d)


text:
TN93(α1::F, α2::F, β::F, πA::F, πC::F, πG::F, πT::F, safe::Bool=true) where F <: Float64 = TN93abs(α1, α2, β, πA, πC, πG, πT, safe)

TN93(κ1::F, κ2::F, πA::F, πC::F, πG::F, πT::F, safe::Bool=true) where F <: Float64 = TN93rel(κ1, κ2, πA, πC, πG, πT, safe)

function TN93(θ_vec::A,
              π_vec::A,
              safe::Bool=true) where A <: AbstractArray
  if safe && length(π_vec) != 4
    error("Incorrect base frequency vector length")
  end
  if length(θ_vec) == 2
    return TN93rel(θ_vec[1], θ_vec[2],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe)
  elseif length(θ_vec) == 3
    return TN93abs(θ_vec[1], θ_vec[2], θ_vec[3],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe)
  else
    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
  end
end

function TN93rel(θ_vec::A,
                 π_vec::A,
                 safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 2
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return TN93rel(θ_vec[1], θ_vec[2], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
end

function TN93abs(θ_vec::A,
                 π_vec::A,
                 safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 3
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return TN93abs(θ_vec[1], θ_vec[2], θ_vec[3], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
end


text:
module Opt_PID_data

export My_Eval, create_My_Eval,   condEntropy, ∇f

using MathProgBase
# http://mathprogbasejl.readthedocs.io/en/latest/nlp.html

struct My_Eval <: MathProgBase.AbstractNLPEvaluator
    n_x     :: Int64
    n_y     :: Int64
    n_z     :: Int64

    n     :: Int64         # number of variables
    m     :: Int64         # number of marginal equations

    varidx  :: Array{Int64,3}                   # 0 if variable not present; otherwise idx of var in {1,...,n}
    xyz     :: Vector{Tuple{Int64,Int64,Int64}} # xyz-triple of varidx

    eqidx   :: Dict{ Tuple{String,Int64,Int64},  Int64} # first idx is "xy" or "xz". Index of eqn in {1,...,m}
    mr_eq   :: Vector{ Tuple{String,Int64,Int64} }      # ("xy",x,y) / ("yz",y,z) of an eqn
    prb_xyz :: Array{Float64,3}
    marg_xy :: Array{Float64,2}
    marg_xz :: Array{Float64,2}

    rhs   :: Vector{Float64}                # m-vector
    Gt    :: SparseMatrixCSC{Float64,Int64} # G^T n x m; transpose of constraint matrix
    Gt_K  :: Vector{Int64}                  # findn()-info of G^T: rows (vars)
    Gt_L  :: Vector{Int64}                  # findn()-info of G^T: columns (equations)

    TmpFloat       :: DataType
    bigfloat_nbits :: Int64
end

function condEntropy{TFloat_2,TFloat}(e::My_Eval, p::Vector{TFloat_2}, dummy::TFloat)   :: TFloat
    # m_yz = marg_yz(e,p,dummy)
    s::TFloat = TFloat(0.)
    for y = 1:e.n_y
        for z = 1:e.n_z
            P_yz = TFloat(0.)
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_yz += p[i]
                end
            end
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    p_xyz = p[i]
                    s  +=  (  (p_xyz ≤ 0 || P_yz  ≤ 0)  ?   TFloat(0.)   :   - p_xyz*log( p_xyz / P_yz )   )
                end
            end
        end
    end
    return s
    ;
end

function ∇f{TFloat,TFloat_2}(e::My_Eval, grad::Vector{TFloat_2}, p::Vector{TFloat_2}, dummy::TFloat) :: Void
    for y = 1:e.n_y
        for z = 1:e.n_z
            # make marginal P(*yz)
            P_yz::TFloat = TFloat(0.)
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_yz += p[i]
                end
            end
            # make log-expressions  log( P(xyz) / P(*yz) )
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_xyz::TFloat = TFloat( p[i] )
                    grad[i] = TFloat_2(   (P_xyz ≤ 0 || P_yz ≤ 0) ?  log(TFloat(e.n_x))  : -log( P_xyz / P_yz )  )
                end
            end
        end# for y
    end# for x
    ;
end # ∇f()

function create_My_Eval(q::Array{Float64,3}, tmpFloat::DataType, bigfloat_precision=256)
    const n_x::Int64 = size(q,1);
    const n_y::Int64 = size(q,2);
    const n_z::Int64 = size(q,3);

    # Create marginals
    prb_xyz::Array{Float64,3}  = zeros(n_x,n_y,n_z)
    marg_xy::Array{Float64,2}  = zeros(n_x,n_y)
    marg_xz::Array{Float64,2}  = zeros(n_x,n_z)

    for x in 1:n_x
        for y in 1:n_y
            for z in 1:n_z
                marg_xy[x,y] += q[x,y,z]
                marg_xz[x,z] += q[x,y,z]
            end
        end
    end
    # Find the variables
    varidx::Array{Int64,3}                   = zeros(Bool,size(q));
    xyz   = Vector{Tuple{Int64,Int64,Int64}}() # = [ (0,0,0) for i in 1:n_x*n_y*n_z ]
    sizehint!(xyz, n_x*n_y*n_z)
    for x in 1:n_x
        for y in 1:n_y
            for z in 1:n_z
                if marg_xy[x,y] > 0  &&  marg_xz[x,z] > 0
                    push!(xyz, (x,y,z))
                    varidx[x,y,z] = endof(xyz)
                    prb_xyz[x,y,z]= q[x,y,z]
                else
                    varidx[x,y,z] = 0
                end#if
            end
        end
    end
    const n::Int64 = length(xyz)

    # Find the equations
    eqidx = Dict{ Tuple{String,Int64,Int64},Int64}() # first idx is "xy" or "xz"
    mr_eq ::Vector{ Tuple{String,Int64,Int64} }   = [ ("",0,0)   for i in 1:n_x*(n_y+n_z) ]
    m::Int64 = 0
    for x in 1:n_x
        for y in 1:n_y
            if marg_xy[x,y] > 0
                m += 1
                eqidx["xy",x,y] = m
                mr_eq[m]        = ("xy",x,y)
            else
                eqidx["xy",x,y] = 0
            end#if
        end
        for z in 1:n_z
            if marg_xz[x,z] > 0
                m += 1
                eqidx["xz",x,z] = m
                mr_eq[m]        = ("xz",x,z)
            else
                eqidx["xz",x,z] = 0
            end#if
        end
    end #for x

    denseGt :: Array{Float64,2} = zeros(n,m)
    for l in 1:n
        (x,y,z) = xyz[l]
        for k in 1:m
            mr = mr_eq[k]
            if mr[1] == "xy"
                xy = mr[2:3]
                if xy[1]==x && xy[2]==y
                    denseGt[l,k] = 1.
                end
            elseif mr[1]=="xz"
                xz = mr[2:3]
                if xz[1]==x && xz[2]==z
                    denseGt[l,k] = 1.
                end
            else
                print("Fuck! Bug!")
                return;
            end
        end
    end

    Gt::SparseMatrixCSC{Float64,Int64} = sparse(denseGt)
    local Gt_K::Array{Int64,1}
    local Gt_L::Array{Int64,1}
    (Gt_K,Gt_L) = findn(Gt)

    TmpFloat       :: DataType  = tmpFloat
    bigfloat_nbits :: Int64     = bigfloat_precision

    rhs::Vector{Float64} = Vector{Float64}(m)

    return My_Eval(n_x,n_y,n_z, n,m, varidx,xyz, eqidx,mr_eq, prb_xyz, marg_xy,marg_xz, rhs, Gt,Gt_K,Gt_L,  TmpFloat,bigfloat_nbits)
    ;
end #^ create_My_eval()

function set_rhs!(ev::My_Eval, p::Array{Float64})::Void
    # Clean out old marginals
    for k in 1:ev.m
        mr = ev.mr_eq[k]
        if mr[1] == "xy"
            ev.marg_xy[ mr[2], mr[3] ] = 0.
        elseif mr[1]=="xz"
            ev.marg_xz[ mr[2], mr[3] ] = 0.
        else
            @assert false
        end
    end

    # Compute new marginals
    for (i,xyz) in enumerate(ev.xyz)
        x,y,z = xyz
        ev.marg_xy[ x,y ]   += max(0.,  p[ x,y,z ]  )
        ev.marg_xz[ x,z ]   += max(0.,  p[ x,y,z ]  )
    end

    # Fill RHS vector
    ev.rhs .= 0.
    for k in 1:ev.m
        mr = ev.mr_eq[k]
        if mr[1] == "xy"
            ev.rhs[k] = max(1.e-300,   ev.marg_xy[ mr[2], mr[3] ]   )  # the max is here because ...
        elseif mr[1]=="xz"
            ev.rhs[k] = max(1.e-300,   ev.marg_xz[ mr[2], mr[3] ]   )  # ... marginals must not be 0
        else
            @assert false
        end
    end #for all marg equations
    ;
end

end #module Opt_PID_data


text:
# Plotting the data structure for the state of an SRA analysis
#   See also: sra_model_object.jl for the saved states
#   See also: sra_model_function.jl for the analysis steps

import PkgSRA.plot_sindy_model

#####
##### Basic data plotting
#####
# function plot_data(m::sra_stateful_object)
#     names = m.sra_parameters
#     plot(m.ts,m.dat[1,:], label=names[1]);
#     for i = 2:size(m.dat,1)
#         plot!(m.ts,m.dat[i,:], label=names[i]);
#     end
#     title!("Raw data")
# end

function plot_data(m::sra_stateful_object, which_dim=1)
    names = m.sindy_model.variable_names
    plot(m.ts,m.dat[which_dim,:], label=names[which_dim]);
    title!("Raw data (variable $which_dim)")
end

function plot_derivative(m::sra_stateful_object, which_dim=1)
    names = m.sindy_model.variable_names
    plot(m.ts,m.numerical_grad[which_dim,:], label=names[which_dim]);
    title!("Numerical derivative (variable $which_dim)")
end

#####
##### Plotting against the truth
#####
function plot_true_grads(m::sra_stateful_object,
            t::sra_truth_object,
            which_dim=1)
    plot(m.ts, t.true_grad[which_dim,:], label="True gradient");
    plot!(m.ts, m.numerical_grad[which_dim,:], label="Numerical gradient")
end

# function plot_data_and_control(m::sra_stateful_object,
#             t::sra_truth_object)
#     plot_data(m)
#     plot!(m.ts, t.U_true, label="True control")
#     title!("Data and controller")
# end

function plot_data_and_control(m::sra_stateful_object,
            t::sra_truth_object; which_dim=1)
    plot_data(m, which_dim)
    plot!(m.ts, t.U_true[which_dim,:], label="True control")
    title!("Data and controller (Variable $which_dim)")
end

function plot_subsampled_points_and_control(m::sra_stateful_object,
            t::sra_truth_object, which_dim=1)
    plot_data_and_control(m, t, which_dim=which_dim)
    dat_sub = m.dat[:, m.subsample_ind]
    scatter!(m.ts[m.subsample_ind], dat_sub[which_dim,:], color=:blue)
    title!("Subsampled points plotted on data (Variable $which_dim)")
end

#####
##### Sindy equation printing
#####
function print_current_equations(m::sra_stateful_object; kwargs...)
    i = m.i
    if i > 1
        println("Current model equations for iteration $i:")
    else
        println("Naive model equations (initial iteration):")
    end
    print_equations(m.sindy_model; kwargs...)
end

function print_true_equations(t::sra_truth_object; kwargs...)
    println("True equations:")
    print_equations(t.core_dyn_true; kwargs...)
end

#####
##### Sindy equation simulations
#####
function plot_sindy_model(m::sra_stateful_object, which_dim=1)
    if m.sindy_dat == nothing
        simulate_model(m)
    end
    plot_data(m, which_dim)
    name = m.sindy_model.variable_names[which_dim]
    plot!(m.ts, m.sindy_dat[which_dim,:],
        label="Simulation ($name)");
    title!("Integration of current model")
end

function plot_sindy_derivatives(m::sra_stateful_object, which_dim=1)
    plot_derivative(m, which_dim)
    name = m.sindy_model.variable_names[which_dim]
    sindy_deriv = m.sindy_model(m.dat, 0)
    plot!(m.ts, sindy_deriv[which_dim,:],
        label="Simulation ($name)");
    title!("Derivatives of current model")
end

#####
##### Subsampling visualizations
#####
function plot_subsampled_points(m::sra_stateful_object, which_dim=1)
    plot_data(m, which_dim)
    dat_sub = m.dat[:, m.subsample_ind]
    scatter!(m.ts[m.subsample_ind], dat_sub[which_dim,:], color=:blue)
    title!("Subsampled points plotted on data (Variable $which_dim)")
end

function plot_residual(m::sra_stateful_object, which_dim=1)
    # TODO: This has to recalculate the residual
    residual = calc_residual(m)

    plot(residual[which_dim,:], label="Residual")
    hline!([m.noise_guess], label="Noise line");
    hline!([-m.noise_guess], label="Noise line")
end

"""
Plots the data and the integrated model.
"""
function plot_subsampled_simulation(m::sra_stateful_object, which_dim=1)
    plot_sindy_model(m, which_dim)
    dat_sub = m.dat[:, m.subsample_ind]
    ts_sub = m.ts[m.subsample_ind]
    scatter!(ts_sub, dat_sub[which_dim,:], color=:blue)
    xlims!((ts_sub[1], ts_sub[end]))
    title!("Subsampled points with integrated model (Variable $which_dim)")
end

"""
Plots the numerical derivative and the derivatives evaluated from the model.
    Note: no integration is performed
"""
function plot_subsampled_derivatives(m::sra_stateful_object, which_dim=1)
    plot_sindy_derivatives(m, which_dim)
    deriv_sub = m.numerical_grad[:, m.subsample_ind]
    ts_sub = m.ts[m.subsample_ind]
    scatter!(ts_sub, deriv_sub[which_dim,:], color=:blue)
    xlims!((ts_sub[1], ts_sub[end]))
    title!("Subsampled points with sindy derivatives (Variable $which_dim)")
end

##### Export
export plot_data, plot_data_and_control, plot_subsampled_points,
    plot_sindy_model, print_current_equations, print_true_equations,
    plot_subsampled_derivatives, plot_residual,
    plot_subsampled_points_and_control, plot_subsampled_simulation


text:
"""
	xy(cmd0::String="", arg1=[]; kwargs...)

reads (x,y) pairs and plot lines, polygons, or symbols at those locations on a map.

Full option list at [`psxy`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html)

Parameters
----------

- **A** : **straight_lines** : -- Str --

    By default, geographic line segments are drawn as great circle arcs. To draw them as straight
    lines, use the -A flag.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#a)
- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** : **offset** : -- Str --

    Offset the plot symbol or line locations by the given amounts dx/dy.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#d)
- **E** : **error_bars** : -- Str --

    Draw symmetrical error bars.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#e)
- **F** : **conn** : **connection** : -- Str --

    Alter the way points are connected
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#f)
- **G** : **fill** : **markerfacecolor** : -- Str --

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#g)
- **I** : **intens** : -- Str or number --

	Use the supplied intens value (in the [-1 1] range) to modulate the fill color by simulating
	shading illumination.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#i)
- **L** : **closed_polygon** : -- Str --

    Force closed polygons.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#l)
- **N** : **no_clip** : -- Str or [] --

    Do NOT clip symbols that fall outside map border
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#n)
- $(GMT.opt_P)
- **S** : -- Str --

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
	[`-S`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#s)

    Alternatively select a sub-set of symbols using the aliases: **symbol** or **marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
	+ **y**, **y_dash**

	and select their sizes with the **markersize** or **size** keyword [default is 8p].
	The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
	points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT="c")*
- **W** : **line_attrib** : **markeredgecolor** : -- Str --

    Set pen attributes for lines or the outline of symbols
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/psxy.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **line_attribs** for the polygons and
    **markeredgecolor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **axis** : **aspect** : -- Str --
    When equal to "equal" makes a square plot.
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function xy(cmd0::String="", arg1=[]; caller=[], K=false, O=false, first=true, kwargs...)
	if (isempty(cmd0) && isa(arg1, Array) && size(arg1,2) == 1 || size(arg1,1) == 1)	# y only
		arg1 = hcat(1:length(arg1), arg1[:])
	end
	common_plot_xyz(cmd0, arg1, caller, K, O, first, false, kwargs...)
end

# ---------------------------------------------------------------------------------------------------
function xyz(cmd0::String="", arg1=[]; caller=[], K=false, O=false, first=true, kwargs...)
	common_plot_xyz(cmd0, arg1, caller, K, O, first, true, kwargs...)
end

# ---------------------------------------------------------------------------------------------------
function common_plot_xyz(cmd0, arg1, caller, K, O, first, is3D, kwargs...)
	arg2 = []		# May be needed if GMTcpt type is sent in via C
	N_args = isempty_(arg1) ? 0 : 1

	if (is3D)	gmt_proggy = "psxyz"
	else		gmt_proggy = "psxy"
	end

	((isempty(cmd0) && isempty_(arg1)) || occursin(" -", cmd0)) && return monolitic(gmt_proggy, cmd0, arg1)

	if (occursin(" -", caller))
		cmd = caller
		caller = []			# It was piggy-backed by scatter or others
	else
		cmd = ""
	end

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	opt_J = " -JX12c/8c"
	for symb in [:axis :aspect]
		if (haskey(d, symb))
			if (d[symb] == "equal" || d[symb] == :equal)	# Need also a 'tight' option?
				opt_J = " -JX12c"
			end
			break
		end
	end
	cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, caller, O, opt_J)
	if (is3D)	cmd,opt_JZ = parse_JZ(cmd, d)	end
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_a(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_swap_xy(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, opt_i = read_data(d, cmd0, cmd, arg1, opt_R, opt_i, opt_bi, opt_di, is3D)

	if (is3D && isempty(opt_JZ) && length(collect(eachmatch(r"/", opt_R))) == 5)
		cmd = cmd * " -JZ6c"	# Default -JZ
	end

	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1, arg2)

	cmd = add_opt(cmd, 'A', d, [:A :straight_lines])
	cmd = add_opt(cmd, 'D', d, [:D :offset])
	cmd = add_opt(cmd, 'E', d, [:E :error_bars])
	cmd = add_opt(cmd, 'F', d, [:F :conn :connection])

	cmd = add_opt(cmd, 'G', d, [:G :fill])
	opt_Gsymb = add_opt("", 'G', d, [:G :markerfacecolor])		# Filling color for symbols

	opt_Wmarker = ""
	if (haskey(d, :markeredgecolor))
		opt_Wmarker = "0.5p," * arg2str(d[:markeredgecolor])	# 0.25p is so thin
	end

	cmd = add_opt(cmd, 'I', d, [:I :intens])
	cmd = add_opt(cmd, 'L', d, [:L :closed_polygon])
	cmd = add_opt(cmd, 'N', d, [:N :no_clip])

	opt_W = ""
	pen = build_pen(d)						# Either a full pen string or empty ("")
	if (!isempty(pen))
		opt_W = " -W" * pen
	else
		for sym in [:W :line_attrib]
			if (haskey(d, sym))
				if (isa(d[sym], Tuple))		# Like this it can hold the pen, not extended atts
					opt_W = " -W" * parse_pen(d[sym])
				else
					opt_W = " -W" * arg2str(d[sym])
				end
				break
			end
		end
	end

	opt_S = add_opt("", 'S', d, [:S :symbol])
	if (isempty(opt_S))			# OK, no symbol given via the -S option. So fish in aliases
		marca = get_marker_name(d, [:marker])
		if (!isempty(marca))
			ms = ""
			for symb in [:size :markersize]
				if (haskey(d, symb))
					if (isa(d[symb], AbstractArray))
						if (length(d[symb]) == size(arg1,1))
							arg1 = hcat(arg1, d[symb][:])
							ms = " "		# Just to defeat the empty test below
						else
							error("The size array must have the same number of elements rows in the data")
						end
					else
						marca = marca * arg2str(d[symb]);	ms = " "
					end
					break
				end
			end
			if (ms == "")  marca = marca * "8p"		end		# Default to 8p
		end
		if (!isempty(marca))  opt_S = " -S" * marca  end
	end

	if (!isempty(opt_S))
		opt_ML = ""
		if (haskey(d, :markerline))
			if (isa(d[:markerline], Tuple))			# Like this it can hold the pen, not extended atts
				opt_ML = " -W" * parse_pen(d[:markerline])
			else
				opt_ML = " -W" * arg2str(d[:markerline])
			end
			if (!isempty(opt_Wmarker))
				opt_Wmarker = ""
				@warn("markerline overrides markeredgecolor")
			end
		end
		if (!isempty(opt_W) && !isempty(opt_ML))
			@warn("You cannot use both markeredgecolor and W or line_attrib keys.")
		end
	end

	if (!isempty(opt_W) && isempty(opt_S)) 			# We have a line/polygon request
		cmd = [finish_PS(d, cmd * opt_W, output, K, O)]
	elseif (isempty(opt_W) && !isempty(opt_S))		# We have a symbol request
		if (!isempty(opt_Wmarker) && isempty(opt_W))
			opt_Gsymb = opt_Gsymb * " -W" * opt_Wmarker	# Piggy back in this option string
		end
		if (!isempty(opt_ML))						# If we have a symbol outline pen
			cmd = cmd * opt_ML
		end
		cmd = [finish_PS(d, cmd * opt_S * opt_Gsymb, output, K, O)]
	elseif (!isempty(opt_W) && !isempty(opt_S))		# We have both line/polygon and a symbol
		# that is not a vector (because Vector width is set by -W)
		if (opt_S[4] == 'v' || opt_S[4] == 'V' || opt_S[4] == '=')
			cmd = [finish_PS(d, cmd * opt_W * opt_S * opt_Gsymb, output, K, O)]
		else
			if (!isempty(opt_Wmarker))
				opt_Wmarker = " -W" * opt_Wmarker	# Set Symbol edge color
			end
			cmd1 = cmd * opt_W
			cmd2 = replace(cmd, opt_B => "") * opt_S * opt_Gsymb * opt_Wmarker	# Don't repeat option -B
			if (!isempty(opt_ML))					# If we have a symbol outline pen
				cmd1 = cmd1 * opt_ML
			end
			cmd = [finish_PS(d, cmd1, output, true, O)
			       finish_PS(d, cmd2, output, K, true)]
		end
	elseif (!isempty(opt_S) && !isempty(opt_ML))		# We have a symbol outline pen
		cmd = [finish_PS(d, cmd * opt_ML * opt_S * opt_Gsymb, output, K, O)]
	else
		cmd = [finish_PS(d, cmd, output, K, O)]
	end

    return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, gmt_proggy, arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
function get_marker_name(d::Dict, symbs, del=false)
	marca = ""
	for symb in symbs
		if (haskey(d, symb))
			t = d[symb]
			if (isa(t, Symbol))	t = string(t)	end
			if     (t == "-" || t == "x-dash")   marca = "-"
			elseif (t == "+" || t == "plus")     marca = "+"
			elseif (t == "a" || t == "*" || t == "star")     marca = "a"
			elseif (t == "c" || t == "circle")   marca = "c"
			elseif (t == "d" || t == "diamond")  marca = "d"
			elseif (t == "g" || t == "octagon")  marca = "g"
			elseif (t == "h" || t == "hexagon")  marca = "h"
			elseif (t == "i" || t == "v" || t == "inverted_tri")  marca = "i"
			elseif (t == "n" || t == "pentagon")  marca = "n"
			elseif (t == "p" || t == "." || t == "point")     marca = "p"
			elseif (t == "r" || t == "rectangle") marca = "r"
			elseif (t == "s" || t == "square")    marca = "s"
			elseif (t == "t" || t == "^" || t == "triangle")  marca = "t"
			elseif (t == "x" || t == "cross")     marca = "x"
			elseif (is3D && (t == "u" || t == "cube"))  marca = "u"
			elseif (t == "y" || t == "y-dash")    marca = "y"
			end
			if (del)
				delete!(d, symb)
			end
			break
		end
	end
	return marca
end

# ---------------------------------------------------------------------------------------------------
xy!(cmd0::String="", arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
	xy(cmd0, arg1; caller=caller, K=K, O=O,  first=first, kw...)
xy!(arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
	xy("", arg1; caller=caller, K=K, O=O,  first=first, kw...)

# ---------------------------------------------------------------------------------------------------
xyz!(cmd0::String="", arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
	xyz(cmd0, arg1; caller=caller, K=K, O=O,  first=first, kw...)
xyz!(arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
	xyz("", arg1; caller=caller, K=K, O=O,  first=first, kw...)

# ---------------------------------------------------------------------------------------------------
psxy   = xy				# Alias
psxy!  = xy!			# Alias
psxyz  = xyz			# Alias
psxyz! = xyz!			# Alias


text:
module DataFrameFuncs

export transpose_dataframe, first_rows, last_rows, find_rows, column_types,
       count_missing, replace_missing!, remove_cols_missing!, describe_number_cols,
       describe_string_cols, describe_date_cols

using DataFrames
using CSV
using Dates

"""
Transpose rows / columns of DataFrame `df`, starting at row `startrow` and using
`nrows` rows. Return value is a DataFrame.
"""
function transpose_dataframe(df::DataFrame,
                             startrow::Integer = 1,
                             nrows::Integer = 1)::DataFrame
    # Check input
    nrows_df::Int64 = size(df, 1)
    if startrow < 1 || nrows < 1
        throw(error("ERROR: invalid start row or number of rows specified, both should be >= 1."))
    elseif startrow > nrows_df
        throw(error("ERROR: start row exceeds number of rows in DataFrame."))
    end
    # Transpose a specific number of rows
    # Check if requested number of rows is available: otherwise limit number of rows
    dft::DataFrame = DataFrame(Column=names(df))
    for row in startrow:minimum([startrow + nrows - 1, nrows_df])
        dft[!, Symbol("Row", row)] = collect(df[row, :])
    end
    return dft
end

"""
Get first `n` rows of DataFrame `df`. Optionally `transpose` rows / columns or limit
the `columns` which are shown in the results. Return value is a DataFrame.
"""
function first_rows(df::DataFrame,
                    n::Integer;
                    transpose::Bool=false,
                    columns::Vector=[])::DataFrame
    dff::DataFrame = first(df, n)
    if !isempty(columns)
        dff = dff[:, columns]
    end
    transpose ? (return transpose_dataframe(dff, 1, n)) : return dff
end

"""
Get last `n` rows of DataFrame `df`. Optionally `transpose` rows / columns or limit
the `columns` which are shown in the results. Return value is a DataFrame.
"""
function last_rows(df::DataFrame,
                   n::Integer;
                   transpose::Bool=false,
                   columns::Vector{Symbol}=[])::DataFrame
    dfl::DataFrame = last(df, n)
    if !isempty(columns)
        dfl = dfl[:, columns]
    end
    transpose ? (return transpose_dataframe(dfl, 1, n)) : return dfl
end

"""
Find rows in DataFrame `df` for which the column `keycolumn` contains the value
`searchvalue`. Optionally limit the `columns` which are shown in the results.
Return value is a DataFrame.
"""
function find_rows(df::DataFrame,
                   keycolumn::Union{String, Symbol},
                   searchvalue::Any; transpose::Bool=false,
                   columns::Vector{Symbol}=[])::DataFrame
    if isempty(columns)
        dfr = df[df[:, keycolumn] .== searchvalue, :]
    else
        dfr = df[df[:, keycolumn] .== searchvalue, columns]
    end
    transpose ? (return transpose_dataframe(dfr, 1, nrow(dfr))) : return dfr
end

"""
Return the columns and their `type` of the columns in DataFrame `df`.
Columns are represented as symbols. Return value should be a sorted Dict(*).

(*) Note that sort of a Dict returns an OrderedCollections.OrderedDict{Symbol, Type}
    instead of a Dict{Symbol, Type} as one would expect. This is a known issue:

    https://github.com/JuliaCollections/OrderedCollections.jl/issues/25
"""
function column_types(df::DataFrame)
    return sort(Dict(Symbol.(names(df)) .=> eltype.(eachcol(df))))
end

"""
Return a `DataFrame` which contains the columns with `missing`
values of DataFrame `df` (and the number of missing values).
"""
function count_missing(df::DataFrame)::DataFrame
    return subset(describe(df, :nmissing), :nmissing => x -> x .> 0)
end

"""
Replace missing values in a `column` of DataFrame `df` with another
`value` (number / string). Changes the original DataFrame.
"""
function replace_missing!(df::DataFrame, column::Symbol, value::Union{Real, String})
    replace!(df[!, column], missing => value);
end

"""
Identify and remove columns with only `missing` values. Changes the original
DataFrame.
"""
function remove_cols_missing!(df::DataFrame)
    dfm::DataFrame = describe(df, :nmissing)
    cols_missing = dfm[dfm[!, :nmissing] .== nrow(df), :][:, 1]
    select!(df, Not(cols_missing))
end

"""
Describe various properties of the number columns in the DataFrame `df`. The output
columns are presented as strings with formatting for better readability.
"""
function describe_number_cols(df::DataFrame)::DataFrame
    # Select number columns from describe DataFrame
    dfn::DataFrame = describe(df, :eltype, :nmissing, :mean, :q25, :median, :q75, :min, :max, sum => :sum)
    number_cols = propertynames(select(df, findall(col -> all(value -> value isa Union{Missing, Number}, col), eachcol(df))))
    dfn = dfn[in(number_cols).(dfn.variable), :] |> sort
    # Change types of number columns (type -> Real)
    foreach(col -> dfn[!, col] = convert(Vector{Real}, dfn[!, col]),
            [:nmissing, :mean, :q25, :median, :q75, :min, :max, :sum])
    return dfn
end

"""
Describe various properties of the string columns in the DataFrame `df`. One of these properties
is a list of unique items, of which the number of items shown is controlled by `max_unique_items`
(default `6`).
"""
function describe_string_cols(df::DataFrame;
                              max_unique_items::Int64=6)::DataFrame
    # Select string columns from describe DataFrame
    dfs = describe(df, :eltype, :nmissing, :nunique)
    string_cols = propertynames(select(df, findall(col -> all(value -> value isa Union{Missing, AbstractString}, col), eachcol(df))))
    dfs = dfs[in(string_cols).(dfs.variable), :]
    # Change types of number columns (type -> Real)
    foreach(col -> dfs[!, col] = convert(Vector{Real}, dfs[!, col]),
            [:nmissing, :nunique])
    # Determine unique items
    vec_unique::Vector{String} = []
    for (idx, col) in enumerate(dfs[:, :variable])
        if dfs[idx, :nunique] <= max_unique_items
            push!(vec_unique, join( unique(df[:, col]) |> sort, ", ") )
        else
            push!(vec_unique, ">$max_unique_items")
        end
    end
    return hcat(dfs, Dict(:unique_items => vec_unique) |> DataFrame)
end

"""
Describe various properties of the date columns in the DataFrame `df`. Dates are formatted
as strings according to the specified `dateFormat` (default `dd-mm-yyyy`).
"""
function describe_date_cols(df::DataFrame;
                            dateFormat::String="dd-mm-yyyy")::DataFrame
    # Select relevant columns from describe DataFrame
    dfd = describe(df, :eltype, :nmissing, :nunique, :min, :max)
    date_cols = propertynames(select(df, findall(col -> all(value -> value isa Union{Missing, Dates.Date}, col), eachcol(df))))
    dfd = dfd[in(date_cols).(dfd.variable), :]
    # Change types of number columns (type -> Real)
    foreach(col -> dfd[!, col] = convert(Vector{Real}, dfd[!, col]),
            [:nmissing, :nunique])
    # Change types of date columns (type -> String)
    foreach(col -> dfd[!, col] = Dates.format.(dfd[!, col], dateFormat),
            [:min, :max])
    return dfd
end

end

max_stars_repo_name: bogersw/julia-snippets

max_issues_repo_name: bogersw/julia-snippets

max_forks_repo_name: bogersw/julia-snippets