text:
# Copyright (c) 2015: Christof Stocker
#
# ==========================================================================
# Influenced by closely looking at the implementation of SVM.jl
#
# Copyright (c) 2013: John Myles White and other contributors.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
# DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
# OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
# USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# ==========================================================================
# Influenced by looking at the implementation of LIBLINEAR
#
# Copyright (c) 2007-2015 The LIBLINEAR Project.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither name of copyright holders nor the names of its contributors
# may be used to endorse or promote products derived from this software
# without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# ==========================================================================

# ==========================================================================
# Specify the custom linear SVM solver

"""
`DualCD <: SvmSolver`

Description
============

Implementation of the dual coordinate descent algorithm as it
was proposed in (Hsieh et al., 2008).

Usage
======

    model = svm(X, y, solver = DualCD(shrinking=false))

References
===========

- Hsieh, Cho-Jui, et al. "A dual coordinate descent method for large-scale linear SVM."
Proceedings of the 25th international conference on Machine learning (ICML-08).
ACM, 2008. DOI=10.1145/1390156.1390208, http://doi.acm.org/10.1145/1390156.1390208
"""
immutable DualCD{T} <: SvmSolver
  shrinking::Bool
  DualCD() = new(T)
end

DualCD(;shrinking::Bool = false) = DualCD{shrinking}()
DualCDWithShrinking() = DualCD{true}()

show(io::IO, ::DualCD{false}) = print(io, "DualCD(shrinking = false)")
show(io::IO, ::DualCD{true}) = print(io, "DualCD(shrinking = true)")

# ==========================================================================
# Implementation of dual coordinate descent for linear L1-, and L2-SVMs
# This particular method is specialized for
# - dense arrays
# - univariate prediction

function fit{TKernel<:ScalarProductKernel, TLoss<:L1orL2HingeLoss, TReg<:L2Penalty, INTERCEPT, T}(
    spec::CSVM{TKernel, TLoss, TReg},
    X::StridedMatrix{T}, y⃗::StridedVector,
    ::DualCD{false},
    predmodel::UvPredicton{INTERCEPT},
    problem::SvmProblem;
    iterations::Int = 1000,
    ftol::Real = 1.0e-4,
    show_trace::Bool = false,
    callback::Union{Function,Void} = nothing,
    nargs...)

  # Do some housekeeping: use shorter variable names for the options
  C = spec.C::Float64
  has_callback = typeof(callback) <: Function
  bias = INTERCEPT ? Float64(predmodel.bias) : zero(Float64)
  dual = typeof(problem) <: DualProblem

  # Print log header if show_trace is set
  show_trace && print_iter_head()

  # Get the size of the design matrix
  #   k ... number of features
  #   l ... number of observations
  k, l = size(X)

  # Initialize α and w to all zero.
  # Note: w has length k+1 to have space for the potential intercept.
  #       If no intercept is fit, the bias term will be ignored
  α  = zeros(l)
  w  = zeros(k + 1)

  # We generate two view to the a.ctual coefficients in w⃗ and the
  # bias in w₀.
  w⃗  = view(w, 1:k)
  w₀ = view(w, k+1)

  # We store the whole gradient if the user is interested in it
  ▽  = has_callback ? zeros(Float64, l) : zeros(Float64, 0)

  # U and D depend on the type of lossfunction used (which is the type of `loss`)
  # Note: Although D is conceptually a matrix, the diagonal Dᵢᵢ is constant and
  #       we can thus save memory by using just a scalar to store it's value
  U, Dᵢᵢ = if typeof(spec.loss) <: HingeLoss
    # L1-SVM: U = C,  Dᵢᵢ = 0
    (C, zero(Float64))
  else # typeof(loss) <: L2HingeLoss
    # L2-SVM: U = ∞,  Dᵢᵢ = (2⋅C)⁻¹
    (typemax(Float64), Float64(.5 / C))
  end

  # The concrete implementation of Q̄  depends on the bias being present or not.
  # This is because the bias is fit seperately instead of being part of the
  # design matrix X.
  # Note: We skip defining Q and go straight for Q̄  = Q + D,
  #       because Q itself is not really required
  Q̄ = if INTERCEPT
    vec(sum(X .^ 2, 1)  + bias ^ 2) + Dᵢᵢ
  else
    vec(sum(X .^ 2, 1)) + Dᵢᵢ
  end

  # Indicies into the observations of X (columns)
  # This array defines the order in which the observations of X will be iterated over
  # Note: It will improve convergence if we permutate this array each iteration.
  #       In order to achieve good performance we preallocate this array and shuffle
  #       its elements inplace each itertation
  indicies = collect(1:l)

  # Preallocate the vector-views into the design matrix to iterate over the observations
  # Note: We have to process each observation of X in every iteration.
  #       This means that even though arrayviews are cheap,
  #       it is still cheaper to preallocate them beforehand.
  X̄ = Array(ContiguousView{T,1,Array{T,2}}, l)
  for i in indicies
    @inbounds X̄[i] = view(X, :, i)
  end

  iteration = 0; G = zero(Float64); PG = zero(Float64)
  PGmin = zero(Float64); PGmax = zero(Float64)
  converged = false; stopped = false
  while !converged && iteration < iterations && !stopped

    # The bounds and difference of the projected gradient of the current iteration
    PGmin = typemax(Float64)
    PGmax = typemin(Float64)

    # Shuffle the indicies to improve convergence
    shuffle!(indicies)

    # loop over all observations
    @inbounds for i in indicies

      # Buffer αᵢ since we need it a few times
      ᾱᵢ = α[i]

      # The concrete computations of G depends on the intercept being present
      # The compiler will remove the if block
      if !INTERCEPT
        G = y⃗[i] .* dot(w⃗, X̄[i]) - 1 + Dᵢᵢ .* ᾱᵢ
      else
        G = y⃗[i] .* (dot(w⃗, X̄[i]) + bias * w₀[1]) - 1. + Dᵢᵢ .* ᾱᵢ
      end

      if has_callback
        # Keep track of the gradient vector if a callback is specified
        ▽[i] = G
      end

      PG = G
      if ᾱᵢ == zero(Float64)
        PG = min(G, zero(Float64))
      elseif ᾱᵢ == U
        PG = max(G, zero(Float64))
      end

      # Keep track of the current gradient bounds
      PGmin = min(PGmin, PG)
      PGmax = max(PGmax, PG)

      # if |PG| != 0 then update αᵢ and w accordingly
      if abs(PG) > 1e-13
        # This will update aᵢ but leave āᵢ alone (we need both to update w)
        α[i] = min(max(ᾱᵢ - G / Q̄[i], zero(Float64)), U)
        # The concrete update of w depends on the intercept being present.
        # The compiler will remove the if block
        #   w⃗ = w⃗ + (αᵢ − ᾱᵢ)yᵢ ⋅ x⃗ᵢ
        #   w⃗ = w⃗ + (αᵢ − ᾱᵢ)yᵢ ⋅ x⃗ᵢ; w₀ += (αᵢ − ᾱᵢ)yᵢ * bias
        if !INTERCEPT
          Base.BLAS.axpy!((α[i] - ᾱᵢ) * y⃗[i], X̄[i], w⃗)
        else
          tmp = (α[i] - ᾱᵢ) * y⃗[i]
          Base.BLAS.axpy!(tmp, X̄[i], w⃗)
          w₀[1] += tmp * bias
        end
      end
    end

    # If the distance of the gradient bounds is close to zero we assume convergence
    converged = (PGmax - PGmin) < ftol
    iteration += 1

    # In case the user requested to print the learning process and/or provided
    # callback function this code will provide just that.
    # Note: If no callback function is provided then the compiler will be able
    #       to optimize the empty callback call away.
    if typeof(callback) <: Function || show_trace
      f = dot(w, w)
      @simd for i in 1:l
        @inbounds f += α[i] * (Dᵢᵢ * α[i] - 2.)
      end
      f = f / 2.
      show_trace && print_iter(iteration, f, PGmax - PGmin)
      stopped = _docallback(predmodel, dual, callback, iteration, w, w⃗, α, f, ▽)
    end
  end

  # Compute the final objective value to store in the result
  f = dot(w, w)
  for i in 1:l
    @inbounds f += α[i] * (Dᵢᵢ * α[i] - 2.)
  end
  f = f / 2.

  #show_trace && print_final(iteration, f, converged)
  if dual
    DualSolution(α, w₀[1], f, iteration, converged)
  else
    # The size of the fitted w depends on if an intercept was fit or not
    PrimalSolution(INTERCEPT ? w : w[1:k], f, iteration, converged)
  end
end

# ==========================================================================
# Implementation of dual coordinate descent for linear L1-, and L2-SVMs
# This particular method is specialized for
# - sparse arrays
# - univariate prediction

function fit{TKernel<:ScalarProductKernel, TLoss<:L1orL2HingeLoss, TReg<:L2Penalty, INTERCEPT}(
    spec::CSVM{TKernel, TLoss, TReg},
    X::SparseMatrixCSC, y⃗::StridedVector,
    ::DualCD{false},
    predmodel::UvPredicton{INTERCEPT},
    problem::SvmProblem;
    iterations::Int = 1000,
    ftol::Real = 1.0e-4,
    show_trace::Bool = false,
    callback::Union{Function,Void} = nothing,
    nargs...)

  # Do some housekeeping: use shorter variable names for the options
  C = spec.C::Float64
  has_callback = typeof(callback) <: Function
  bias = INTERCEPT ? Float64(predmodel.bias) : zero(Float64)
  dual = typeof(problem) <: DualProblem

  # Print log header if show_trace is set
  show_trace && print_iter_head()

  # Get the size of the design matrix
  #   k ... number of features
  #   l ... number of observations
  k, l = size(X)

  # Initialize α and w to all zero.
  # Note: w has length k+1 to have space for the potential intercept.
  #       If no intercept is fit, the bias term will be ignored
  α  = zeros(l)
  w  = zeros(k + 1)

  # We generate two view to the actual coefficients in w⃗ and the
  # bias in w₀.
  w⃗  = view(w, 1:k)
  w₀ = view(w, k+1)

  # We store the gradient if the user is interested in it
  ▽  = has_callback ? zeros(l) : zeros(0)

  # U and D depend on the type of lossfunction used (which is the type of `loss`)
  # Note: Although D is conceptually a matrix, the diagonal Dᵢᵢ is constant and
  #       we can thus save memory by using just a scalar to store it's value
  U, Dᵢᵢ = if typeof(spec.loss) <: HingeLoss
    # L1-SVM: U = C,  Dᵢᵢ = 0
    (C, zero(Float64))
  else # typeof(loss) <: L2HingeLoss
    # L2-SVM: U = ∞,  Dᵢᵢ = (2⋅C)⁻¹
    (typemax(Float64), convert(Float64, .5 / C))
  end

  # Indicies into the observations of X (columns)
  # This array defines the order in which the observations of X will be iterated over
  # Note: It will improve convergence if we permutate this array each iteration.
  #       In order to achieve good performance we preallocate this array and shuffle
  #       its elements inplace each iertation
  indicies = collect(1:l)

  # We compute Q̄  as efficiently as possible by using the inner structure
  # of the SparseMatrixCSC directly
  # Note: We skip defining Q and go straight for Q̄  = Q + D,
  #       because Q itself is not really required
  Q̄ = zeros(l)
  @inbounds for i in indicies
    tstart = X.colptr[i]
    tstop  = X.colptr[i+1] - 1
    for j = tstart:tstop
      Q̄[i] += (X.nzval[j])^2
    end
    if INTERCEPT
      Q̄[i] += bias ^ 2
    end
    Q̄[i] += Dᵢᵢ
  end

  iteration = 0; G = zero(Float64); PG = zero(Float64)
  PGmin = zero(Float64); PGmax = zero(Float64)
  converged = false; stopped = false
  while !converged && iteration < iterations && !stopped

    # The bounds and difference of the projected gradient of the current iteration
    PGmin = typemax(Float64)
    PGmax = typemin(Float64)

    # Shuffle the indicies to improve convergence
    shuffle!(indicies)

    # loop over all observations
    @inbounds for i in indicies

      # These two variables define where the column of interest and how many
      # elements that column contains. This way we can efficiently iterate
      # over exactly those values that correspond to the non-zero features
      # of the current column (i.e. observation)
      tstart = X.colptr[i]
      tstop  = X.colptr[i+1] - 1

      # Buffer αᵢ since we need it a few times
      ᾱᵢ = α[i]

      # Compute G efficiently by working directly with the inner structure
      # of the SparseMatrixCSC
      #   G = yᵢ * (w⃗ ⋅ Xᵢ + bias * w₀) - 1 + Dᵢᵢ * αᵢ
      G = zero(Float64)
      # This does: "w⃗ ⋅ Xᵢ"
      for j = tstart:tstop
        G += w⃗[X.rowval[j]] * X.nzval[j]
      end
      # This does: "(... + bias * w₀)"
      if INTERCEPT
        G += bias * w₀[1]
      end
      # This does: "yᵢ * (...)"
      G *= y⃗[i]
      # This does: "... - 1 + Dᵢᵢ * αᵢ"
      G += Dᵢᵢ .* ᾱᵢ - 1.

      # Keep track of the gradient vector if a callback is specified
      if has_callback
        ▽[i] = G
      end

      PG = G
      if ᾱᵢ == zero(Float64)
        PG = min(G, zero(Float64))
      elseif ᾱᵢ == U
        PG = max(G, zero(Float64))
      end

      # Keep track of the current gradient bounds
      PGmin = min(PGmin, PG)
      PGmax = max(PGmax, PG)

      # if PḠ != 0 then update αᵢ and w accordingly
      if abs(PG) > 1e-13
        # This will update aᵢ but leave āᵢ alone (we need both to update w)
        α[i] = min(max(ᾱᵢ - G / Q̄[i], zero(Float64)), U)
        # Again, we use the structure of SparseMatrixCSC to compute w efficiently
        #   w⃗ = w⃗ + (αᵢ − ᾱᵢ)yᵢ ⋅ x⃗ᵢ
        tmp = (α[i] - ᾱᵢ) * y⃗[i]
        for j = tstart:tstop
          w⃗[X.rowval[j]] += tmp * X.nzval[j]
        end
        #   w₀ += (αᵢ − ᾱᵢ)yᵢ * bias
        if INTERCEPT
          w₀[1] += tmp * bias
        end
      end
    end

    # If the distance of the gradient bounds is close to zero we assume convergence
    converged = (PGmax - PGmin) < ftol
    iteration += 1

    # In case the user requested to print the learning process and/or provided
    # callback function this code will provide just that.
    # Note: If no callback function is provided then the compiler will be able
    #       to optimize the empty callback call away.
    if typeof(callback) <: Function || show_trace
      f = dot(w, w)
      for i in 1:l
        @inbounds f += α[i] * (Dᵢᵢ * α[i] - 2.)
      end
      f = f / 2.
      show_trace && print_iter(iteration, f, PGmax - PGmin)
      stopped = _docallback(predmodel, dual, callback, iteration, w, w⃗, α, f, ▽)
    end
  end

  # Compute the final objective value to store in the result
  f = dot(w, w)
  for i in 1:l
    @inbounds f += α[i] * (Dᵢᵢ * α[i] - 2.)
  end
  f = f / 2.

  #show_trace && print_final(iteration, f, converged)
  if dual
    DualSolution(α, w₀[1], f, iteration, converged)
  else
    # The size of the fitted w depends on if an intercept was fit or not
    PrimalSolution(INTERCEPT ? w : w[1:k], f, iteration, converged)
  end
end


text:
# ------------------------------------------------------------------------------------------
# # Multiple neural network layers with `Flux.jl`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In a previous notebook, we saw that one layer of neurons wasn't enough to distinguish
# between three types of fruit (apples, bananas *and* grapes), since the data is quite
# complex. To solve this problem, we need to use more layers, so heading into the territory
# of **deep learning**!
#
# By adding another layer between the inputs and the output neurons, a so-called "hidden
# layer", we will get our first serious **neural network**, looking something like this:
# ------------------------------------------------------------------------------------------

include("draw_neural_net.jl")
draw_network([2, 4, 3])

# ------------------------------------------------------------------------------------------
# We will continue to use two input data and try to classify into three types, so we will
# have three output neurons. We have chosen to add a single "hidden layer" in between, and
# have arbitrarily chosen to put 4 neurons there.
#
# Much of the *art* of deep learning is choosing a suitable structure for the neural network
# that will allow the model to be sufficiently complex to model the data well, but
# sufficiently simple to allow the parameters to be learned in a reasonable length of time.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Read in and process data
#
# As before, let's load some pre-processed data using code we've seen in the previous
# notebook.
# ------------------------------------------------------------------------------------------

using Flux
using Flux: onehot

using CSV

apples_1 = CSV.read("data/Apple_Golden_1.dat", delim='\t')
apples_2 = CSV.read("data/Apple_Golden_2.dat", delim='\t')
apples_3 = CSV.read("data/Apple_Golden_3.dat", delim='\t')
bananas = CSV.read("data/Banana.dat", delim='\t')
grapes_1 = CSV.read("data/Grape_White.dat", delim='\t')
grapes_2 = CSV.read("data/Grape_White_2.dat", delim='\t');

apples = vcat(apples_1, apples_2, apples_3)
grapes = vcat(grapes_1, grapes_2);

col1 = :red
col2 = :blue

x_apples  = [ [apples_1[i, col1], apples_1[i, col2]] for i in 1:size(apples_1)[1] ]
append!(x_apples, [ [apples_2[i, col1], apples_2[i, col2]] for i in 1:size(apples_2)[1] ])
append!(x_apples, [ [apples_3[i, col1], apples_3[i, col2]] for i in 1:size(apples_3)[1] ])

x_bananas = [ [bananas[i, col1], bananas[i, col2]] for i in 1:size(bananas)[1] ]

x_grapes = [ [grapes_1[i, col1], grapes_1[i, col2]] for i in 1:size(grapes_1)[1] ]
append!(x_grapes, [ [grapes_2[i, col1], grapes_2[i, col2]] for i in 1:size(grapes_2)[1] ])

xs = vcat(x_apples, x_bananas, x_grapes);

# ------------------------------------------------------------------------------------------
# We now we wish to classify the three types of fruit, so we again use one-hot vectors to
# represent the desired outputs $y^{(i)}$:
# ------------------------------------------------------------------------------------------

labels = [ones(length(x_apples)); 2*ones(length(x_bananas)); 3*ones(length(x_grapes))];

ys = [onehot(label, 1:3) for label in labels];  # onehotbatch(labels, 1:3)

# ------------------------------------------------------------------------------------------
# The input data is in `xs` and the one-hot vectors are in `ys`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Multiple layers in Flux
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's tell Flux what structure we want the network to have. We first specify the number of
# neurons in each layer, and then construct each layer as a `Dense` layer:
# ------------------------------------------------------------------------------------------

inputs = 2
hidden = 4
outputs = 3

layer1 = Dense(inputs, hidden, σ)
layer2 = Dense(hidden, outputs, σ)

# ------------------------------------------------------------------------------------------
# To stitch together multiple layers to make a multi-layer network, we use Flux's `Chain`
# function:
# ------------------------------------------------------------------------------------------

model = Chain(layer1, layer2)

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# What is the internal structure and sub-structure of this `model` object?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Training the model
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We have now set up a model and we have some training data.
# How do we train the model on the data?
#
# The amazing thing is that the rest of the code in `Flux` is **exactly the same as
# before**. This is possible thanks to the design of Julia itself, and of the `Flux`
# package.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Train the model as before, now using the popular `ADAM` optimizer. You may need to train
# the network for longer than before, since we have many more parameters.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Visualizing the results
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# What does this neural network represent? It is simply a more complicated function with two
# inputs and three outputs, i.e. a function $f: \mathbb{R}^2 \to \mathbb{R}^3$.
# Before, with a single layer, each component of the function $f$ basically corresponded to
# a hyperplane; now it will instead be a **more complicated nonlinear function** of the
# input data!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Visualize each component of the output separately as a heatmap and/or contours
# superimposed on the data. Interpret the results.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## What we have learned
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Adding an intermediate layer allows the network to start to deform the separating surfaces
# that it is learning into more complicated, nonlinear (curved) shapes. This allows it to
# separate data that were previously unable to be separated!
#
# However, using only two features means that data from different classes overlaps. To
# distinguish it we would need to use more features.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Exercise 4
#
# Use three features (red, green and blue) and build a network with one hidden layer. Does
# this help to distinguish the data better?
# ------------------------------------------------------------------------------------------

ext: jl

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/18. Tools - Multi-layer neural networks with Flux.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
export normlogpdf, mvnormlogpdf

using Flux: unsqueeze, stack

# watch https://github.com/JuliaStats/Distributions.jl/issues/1183
const log2π = log(2.0f0π)

"""
     normlogpdf(μ, σ, x; ϵ = 1.0f-8)

GPU automatic differentiable version for the logpdf function of normal distributions.
Adding an epsilon value to guarantee numeric stability if sigma is exactly zero
(e.g. if relu is used in output layer).
"""
function normlogpdf(μ, σ, x; ϵ = 1.0f-8)
    z = (x .- μ) ./ (σ .+ ϵ)
    -(z .^ 2 .+ log2π) / 2.0f0 .- log.(σ .+ ϵ)
end

"""
    mvnormlogpdf(μ::AbstractVecOrMat, L::AbstractMatrix, x::AbstractVecOrMat)

GPU automatic differentiable version for the logpdf function of multivariate
normal distributions.  Takes as inputs `mu` the mean vector, `L` the lower
triangular matrix of the cholesky decomposition of the covariance matrix, and
`x` a matrix of samples where each column is a sample.  Return a Vector
containing the logpdf of each column of x for the `MvNormal` parametrized by `μ`
and `Σ = L*L'`.
"""
function mvnormlogpdf(μ::AbstractVecOrMat, L::AbstractMatrix, x::AbstractVecOrMat)
    return -(
        (size(x, 1) * log2π + logdetLorU(L)) .+ vec(sum(abs2.(L \ (x .- μ)), dims = 1))
    ) ./ 2
end

"""
    mvnormlogpdf(μ::A, LorU::A, x::A; ϵ = 1f-8) where A <: AbstractArray

Batch version that takes 3D tensors as input where each slice along the 3rd
dimension is a batch sample.  `μ` is a (action_size x 1 x batch_size) matrix,
`L` is a (action_size x action_size x batch_size), x is a (action_size x
action_samples x batch_size).  Return a 3D matrix of size (1 x action_samples x
batch_size).
"""
function mvnormlogpdf(μ::A, LorU::A, x::A; ϵ = 1.0f-8) where {A<:AbstractArray}
    logp = [mvnormlogpdf(μ[:, :, k], LorU[:, :, k], x[:, :, k]) for k in 1:size(x, 3)]
    return unsqueeze(stack(logp, 2), 1) #returns a 3D vector
end


text:
# benchmarks for ISDF in BSE

# %%
# loading packages
using BenchmarkTools, JLD2, FileIO, LinearAlgebra, FFTW
using BSE_k_ISDF

BLAS.set_num_threads(1)
FFTW.set_num_threads(1)

# %% benchmark different values of N_k

# fixed parameters
example_path = "diamond/"
N_1d = 20
N_rs = (N_1d, N_1d, N_1d)
N_core = 0
N_v = 4
N_c = 10

N_μ_vv = 70
N_μ_cc = 220
N_μ_vc = 100

#variable parameters
N_ks_vec = [(2, 2, 2), (3, 3, 3), (4, 4, 4), (5, 5, 5), (7, 7, 7), (9, 9, 9), (13, 13, 13)]

for N_ks in N_ks_vec
    prob = BSE_k_ISDF.BSEProblemExciting(N_core, N_v, N_c, N_ks, N_rs, example_path * "$(N_ks...)_$(N_1d)/")

    isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)
    function f()
        BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)
    end
    t_isdf = @benchmark $f()

    function g()
        return BSE_k_ISDF.setup_H(prob, isdf)
    end
    t_H_setup = @benchmark $g()

    H = BSE_k_ISDF.setup_H(prob, isdf);
    x = rand(N_v * N_c * prod(N_ks)) + im * rand(N_v * N_c * prod(N_ks))
    function h()
        return H * x
    end
    t_H_x = @benchmark $h()

    save(example_path * "$(N_ks...)_$(N_1d)/benchmark_$(N_μ_vv)_$(N_μ_cc)_$(N_μ_vc).jld2", "t_isdf", t_isdf, "t_H_setup", t_H_setup, "t_H_x", t_H_x)
end

# %% save results

results = [load("diamond/$(N_ks...)_$(N_1d)/benchmark_$(N_μ_vv)_$(N_μ_cc)_$(N_μ_vc).jld2") for N_ks in N_ks_vec]
timings = ["t_isdf", "t_H_setup", "t_H_x"]

setup_times = sum(1e-9 .* time.(minimum.([res[t] for res in results, t in timings[1:2]])); dims = 2)
evaluation_times = 1e-9 .* time.(minimum.([res[t] for res in results, t in timings[3:3]]))

save(example_path * "/benchmark_$(N_μ_vv)_$(N_μ_cc)_$(N_μ_vc).jld2", "N_ks_vec", N_ks_vec, "setup_times", setup_times, "evaluation_times", evaluation_times)


text:
module TableIOInterface

export get_file_type, get_example_code, is_extension_supported, multiple_tables

## definition of file formats and extensions

abstract type AbstractFormat end

struct CSVFormat <: AbstractFormat end
struct ZippedFormat <: AbstractFormat end
struct JDFFormat <: AbstractFormat end
struct ParquetFormat <: AbstractFormat end
struct ExcelFormat <: AbstractFormat end
struct SQLiteFormat <: AbstractFormat end
struct StataFormat <: AbstractFormat end
struct SPSSFormat <: AbstractFormat end
struct SASFormat <: AbstractFormat end
struct JSONFormat <: AbstractFormat end
struct ArrowFormat <: AbstractFormat end
struct HDF5Format <: AbstractFormat end
struct JLD2Format <: AbstractFormat end
struct JuliaFormat <: AbstractFormat end

# data base only formats - not completely integrated yet
struct PostgresFormat <: AbstractFormat end

# define if copy of columns is required when converting Tables.jl interface to DataFrame
require_copycols(::AbstractFormat) = false
require_copycols(::JSONFormat) = true

# specify if mutliple tables could be inside a single file
multiple_tables(::AbstractFormat) = false

multiple_tables(::ZippedFormat) = true
multiple_tables(::ExcelFormat) = true
multiple_tables(::SQLiteFormat) = true
multiple_tables(::HDF5Format) = true
multiple_tables(::JLD2Format) = true
multiple_tables(::PostgresFormat) = true
multiple_tables(::JuliaFormat) = true

# mapping of file extensions to table file formats
# multiple extensions can be mapped to the same format
const FILE_EXTENSIONS = Dict{String, DataType}(
    "zip" => ZippedFormat,
    "csv" => CSVFormat,
    "jdf" => JDFFormat,
    "parquet" => ParquetFormat,
    "xlsx" => ExcelFormat,
    "xlsm" => ExcelFormat,
    "db" => SQLiteFormat,
    "sqlite" => SQLiteFormat,
    "sqlite3" => SQLiteFormat,
    "dta" => StataFormat,
    "sav" => SPSSFormat,
    "sas7bdat" => SASFormat,
    "json" => JSONFormat,
    "arrow" => ArrowFormat,
    "hdf" => HDF5Format,
    "hdf5" => HDF5Format,
    "jld2" => JLD2Format,
    "jl" => JuliaFormat,
)

_get_file_extension(filename:: AbstractString) = lowercase(splitext(filename)[2][2:end])

get_file_type(filename:: AbstractString) = FILE_EXTENSIONS[_get_file_extension(filename)]()

is_extension_supported(extension:: AbstractString) = lowercase(extension) ∈ keys(FILE_EXTENSIONS)

const import_table_io = "import TableIO"
const import_dataframes = "import DataFrames"

# `directory` is not used anymore because it corresponds to the notebook directory when inserting the data,
# not the current notebook directory.
get_example_code(t::T, directory:: AbstractString, filename:: AbstractString) where {T <: AbstractFormat} = get_example_code(t, filename)

function get_example_code(t::T, filename:: AbstractString) where {T <: AbstractFormat}
    if multiple_tables(t)
        """
        df_$(_get_varname(filename)) = let
            $import_table_io
            $import_dataframes
            # table_list = TableIO.list_tables($(_get_path_code(filename))) # uncomment to get a list of all tables in this file
            DataFrames.DataFrame(TableIO.read_table($(_get_path_code(filename)),
                # "tablename", # define the table name here, if not defined the alphabetically first table is loaded
                ); copycols=$(require_copycols(t)))
        end"""
    else
        """
        df_$(_get_varname(filename)) = let
            $import_table_io
            $import_dataframes
            DataFrames.DataFrame(TableIO.read_table($(_get_path_code(filename))); copycols=$(require_copycols(t)))
        end"""
    end
end

_get_path_code(filename) = """joinpath(split(@__FILE__, '#')[1] * ".assets", "$filename")"""

_get_varname(filename) = replace(filename, r"[\"\-,\.#@!\%\s+\;()\$&*\[\]\{\}'^]" => "_")

get_example_code(directory:: AbstractString, filename:: AbstractString, args...; kwargs...) = get_example_code(get_file_type(filename), directory, filename, args...; kwargs...)
end


text:
@diffusion_process JansenRitAux{K, R} begin
    :dimensions
    process --> 6
    wiener --> 1
    :parameters
    (a, b, σy) --> K
    :auxiliary_info
    t0 --> Float64
    T --> Float64
    vT --> R

    :additional
    constdiff --> true
    linear --> true
end
# θ
#param_names(P::JansenRitAux) = (:a, :b,, :σy)
#θ =[100.0, 50.0, 2000.0]
#P = JansenRitAux(θ...)

function B(t, P::JansenRitAux)
    @SMatrix [  0.0 0.0 0.0 1.0 0.0 0.0;
                0.0 0.0 0.0 0.0 1.0 0.0;
                0.0 0.0 0.0 0.0 0.0 1.0;
                -P.a^2 0.0 0.0 -2.0*P.a 0.0 0.0;
                0.0 -P.a^2 0.0 0.0 -2.0*P.a 0.0;
                0.0 0.0 -P.b^2 0.0 0.0 -2.0*P.b]
end

β(t, P::JansenRitAux) = @SVector [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

function σ(t, P::JansenRitAux)
    @SVector [0.0, 0.0, 0.0, 0.0, P.σy, 0.0]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpdateTrustedIdProviderWithAccountParameters <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, UpdateTrustedIdProviderProperties } # spec name: properties

    function UpdateTrustedIdProviderWithAccountParameters(;name=nothing, properties=nothing)
        o = new()
        validate_property(UpdateTrustedIdProviderWithAccountParameters, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(UpdateTrustedIdProviderWithAccountParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type UpdateTrustedIdProviderWithAccountParameters

const _property_map_UpdateTrustedIdProviderWithAccountParameters = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"))
const _property_types_UpdateTrustedIdProviderWithAccountParameters = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("properties")=>"UpdateTrustedIdProviderProperties")
Base.propertynames(::Type{ UpdateTrustedIdProviderWithAccountParameters }) = collect(keys(_property_map_UpdateTrustedIdProviderWithAccountParameters))
Swagger.property_type(::Type{ UpdateTrustedIdProviderWithAccountParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpdateTrustedIdProviderWithAccountParameters[name]))}
Swagger.field_name(::Type{ UpdateTrustedIdProviderWithAccountParameters }, property_name::Symbol) =  _property_map_UpdateTrustedIdProviderWithAccountParameters[property_name]

function check_required(o::UpdateTrustedIdProviderWithAccountParameters)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ UpdateTrustedIdProviderWithAccountParameters }, name::Symbol, val)
end


text:
export
    Creation,
    create!,
    create,
    Annihilation,
    annihilate!,
    annihilate,

    ComplexVec,
    α,
    ξ,

    Displacement,
    displace!,

    Squeezing,
    squeeze!

############
# a† and a #
############

"""
    Creation(; dim=DIM)

Creation operator in matrix representation

``\\hat{a}^{\\dagger}``
"""
Creation(T::Type{<:Number}; dim=DIM) = diagm(-1 => sqrt.(T.(1:dim-1)))
Creation(; dim=DIM) = Creation(ComplexF64, dim=dim)

"""
    create!(state::AbstractState)

Apply creation operator on the quantum state.

# Examples
```jldoctest
julia> state = VacuumState();

julia> create!(state);

julia> vec(state) == vec(SinglePhotonState())
true
```
"""
function create!(state::StateVector{T}) where {T}
    dim = state.dim
    state.v = Creation(T, dim=dim) * state.v

    return state
end

function create!(state::StateMatrix{T}) where {T}
    dim = state.dim
    𝐜 = Creation(T, dim=dim)
    state.𝛒 = 𝐜 * state.𝛒 * 𝐜'

    return state
end

"""
    create(state::AbstractState)

Apply creation operator on the new instance of the quantum state.

# Examples
```jldoctest
julia> state = VacuumState();

julia> new_state = create(state);

julia> vec(state) == vec(VacuumState())
true

julia> vec(new_state) == vec(SinglePhotonState())
true
```
"""
create(state::AbstractState) = create!(copy(state))

"""
    Annihilation(; dim=DIM)

Annihilation operator in matrix representation

``\\hat{a}``
"""
Annihilation(T::Type{<:Number}; dim=DIM) = diagm(1 => sqrt.(T.(1:dim-1)))
Annihilation(; dim=DIM) = Annihilation(ComplexF64, dim=dim)

"""
    annihilate!(state::AbstractState)

Apply annihilation operator on the quantum state.

# Examples
```jldoctest
julia> state = SinglePhotonState();

julia> annihilate!(state);

julia> vec(state) == vec(VacuumState())
true
```
"""
function annihilate!(state::StateVector{T}) where {T}
    dim = state.dim
    state.v = Annihilation(T, dim=dim) * state.v

    return state
end

function annihilate!(state::StateMatrix{T}) where {T}
    dim = state.dim
    𝐚 = Annihilation(T, dim=dim)
    state.𝛒 = 𝐚 * state.𝛒 * 𝐚'

    return state
end

"""
    annihilate!(state::AbstractState)

Apply annihilation operator on the new instance of the quantum state.

# Examples
```jldoctest
julia> state = SinglePhotonState();

julia> new_state = annihilate(state);

julia> vec(state) == vec(SinglePhotonState())
true

julia> vec(new_state) == vec(VacuumState())
true
```
"""
annihilate(state::AbstractState) = annihilate!(copy(state))

###########
# α and ξ #
###########

"""
    ComplexVec{T <: Real}(r::T, θ::T)

Vector in polar coordinate for complex plane.

``v = r e^{-i\\theta}``
"""
struct ComplexVec{T<:Real}
    r::T
    θ::T
end

function Base.show(io::IO, complexvec::ComplexVec{T}) where {T}
    print(io, "ComplexVec{$T}($(complexvec.r)exp(-$(complexvec.θ)im))")
end

z(complexvec::ComplexVec) = complexvec.r * exp(-im * complexvec.θ)

"""
    α(r::Real, θ::Real)

Eigenvalue of annihilation operator.

``\\hat{a} | \\alpha \\rangle = \\alpha | \\alpha \\rangle``

# Examples
```jldoctest
julia> α(1.5, π/4)
ComplexVec{Float64}(1.5exp(-0.7853981633974483im))
```
"""
α(T::Type{<:Real}, r::Real, θ::Real) = ComplexVec{T}(T(r), T(θ))
α(r, θ) = α(Float64, r, θ)

"""
    ξ(r::Real, θ::Real)

# Examples
```jldoctest
julia> ξ(1.5, π/4)
ComplexVec{Float64}(1.5exp(-0.7853981633974483im))
```
"""
const ξ = α

################
# displacement #
################

"""
    Displacement(α::ComplexVec{<:Real}; dim=DIM)

Displacement operator in matrix representation

``\\hat{D}(\\alpha) = exp(\\alpha \\hat{a}^{\\dagger} - \\alpha^{*} \\hat{a})``
"""
function Displacement(T::Type{<:Number}, α::ComplexVec; dim=DIM)
    return exp(z(α) * Creation(T, dim=dim) - z(α)' * Annihilation(T, dim=dim))
end

Displacement(α::ComplexVec; dim=DIM) = Displacement(ComplexF64, α, dim=dim)

"""
    displace!(state::AbstractState)

Apply displacement operator on the quantum state.

# Examples
```jldoctest
julia> state = VacuumState();

julia> displace!(state, α(5., π/4));

julia> vec(state) == vec(CoherentState(α(5., π/4)))
true
```
"""
function displace!(state::StateVector{T}, α::ComplexVec) where {T}
    dim = state.dim
    state.v = Displacement(T, α, dim=dim) * state.v

    return state
end

function displace!(state::StateMatrix{T}, α::ComplexVec) where {T}
    dim = state.dim
    𝐝 = Displacement(T, α, dim=dim)
    state.𝛒 = 𝐝 * state.𝛒 * 𝐝'

    return state
end

#############
# squeezing #
#############

"""
    Squeezing(ξ::ComplexVec{<:Real}; dim=DIM)

Squeezing operator in matrix representation

``\\hat{S}(\\xi) = exp(\\frac{1}{2} (\\xi^{*} \\hat{a}^{2} - \\xi \\hat{a}^{\\dagger 2}))``
"""
function Squeezing(T::Type{<:Number}, ξ::ComplexVec; dim=DIM)
    return exp(0.5 * z(ξ)' * Annihilation(T, dim=dim)^2 - 0.5 * z(ξ) * Creation(T, dim=dim)^2)
end

Squeezing(ξ::ComplexVec; dim=DIM) = Squeezing(ComplexF64, ξ, dim=dim)

"""
    squeeze!(state::AbstractState)

Apply squeezing operator on the quantum state.

# Examples
```jldoctest
julia> state = VacuumState();

julia> squeeze!(state, α(0.5, π/4));

julia> vec(state) == vec(SqueezedState(ξ(0.5, π/4)))
true
```
"""
function squeeze!(state::StateVector{T}, ξ::ComplexVec) where {T}
    dim = state.dim
    state.v = Squeezing(T, ξ, dim=dim) * state.v

    return state
end

function squeeze!(state::StateMatrix{T}, ξ::ComplexVec) where {T}
    dim = state.dim
    𝐬 = Squeezing(T, ξ, dim=dim)
    state.𝛒 = 𝐬 * state.𝛒 * 𝐬'

    return state
end

###################
# BHD measurement #
###################

##### for arb. state in intensity-to-measurement-phase quadrature coordinate #####

# |θ, x⟩ = ∑ₙ |n⟩ ⟨n|θ, x⟩ = ∑ₙ ψₙ(θ, x) |n⟩
# ⟨n|θ, x⟩ = ψₙ(θ, x) = exp(im n θ) (2/π)^(1/4) exp(-x^2) Hₙ(√2 x)/√(2^n n!)
function ψₙ(n::Integer, θ::Real, x::Real)
    return (2/π)^(1/4) * exp(im*n*θ - x^2) * hermiteh(n, sqrt(2)x) / sqrt(2^n * factorial(n))
end

function 𝛑̂!(result::Matrix{<:Complex}, θ::Real, x::Real; dim=DIM)
    view(result, :, 1) .= ψₙ.(big(0):big(dim-1), θ, x)
    result .= view(result, :, 1) * view(result, :, 1)'

    return result
end

function 𝛑̂(T::Type{<:Complex}, θ::Real, x::Real; dim=DIM)
    result = Matrix{T}(undef, dim, dim)

    return 𝛑̂!(result, θ, x, dim=dim)
end

𝛑̂(θ::Real, x::Real; dim=DIM) = 𝛑̂(ComplexF64, θ, x, dim=dim)

##### for Gaussian state in intensity-to-measurement-phase quadrature coordinate #####

# π̂ₓ = (â exp(-im θ) + â† exp(im θ)) / 2

tr_mul(𝐚, 𝐛) = sum(𝐚[i, :]' * 𝐛[:, i] for i in 1:size(𝐚, 1))
create_μ(state::StateMatrix{T}) where {T} = tr_mul(Creation(T, dim=state.dim), state.𝛒)
create²_μ(state::StateMatrix{T}) where {T} = tr_mul(Creation(T, dim=state.dim)^2, state.𝛒)
annihilate_μ(state::StateMatrix{T}) where {T} = tr_mul(Annihilation(T, dim=state.dim), state.𝛒)
annihilate²_μ(state::StateMatrix{T}) where {T} = tr_mul(Annihilation(T, dim=state.dim)^2, state.𝛒)
create_annihilate_μ(state::StateMatrix{T}) where {T} = tr_mul(
    Creation(T, dim=state.dim) * Annihilation(T, dim=state.dim),
    state.𝛒
)

# ⟨π̂ₓ²⟩ = ⟨ââ exp(-2im θ) + â†â† exp(2im θ) + ââ† + â†â⟩ / 4
# ⟨π̂ₓ²⟩ = (exp(-2im θ)⟨â²⟩ + exp(2im θ)⟨â†²⟩ + 1 + 2⟨ââ†⟩) / 4
# here, ⟨ââ† + â†â⟩ = 1 + 2⟨ââ†⟩ due to the commutation relation
function π̂ₓ²_μ(θs::AbstractVector{<:Number}, state::StateMatrix)
    return (
        exp.(-2im*θs) .* annihilate²_μ(state) .+
        exp.(2im*θs) .* create²_μ(state) .+
        1 .+ 2create_annihilate_μ(state)
    ) ./ 4
end

# ⟨π̂ₓ⟩ = ⟨â exp(-im θ) + â† exp(im θ)⟩ / 2
# ⟨π̂ₓ⟩ = (exp(-im θ)⟨â⟩ + exp(im θ)⟨â†⟩) / 2
function π̂ₓ_μ(θs::AbstractVector{<:Number}, state::StateMatrix)
    return (
        exp.(-im*θs) .* annihilate_μ(state) .+
        exp.(im*θs) .* create_μ(state)
    ) ./ 2
end


text:
const INDENT = 4

mutable struct IndentState
    indent::Int
    edits
end

function indent_pass(x, state)
    if headof(x) === :file
        if length(x) > 0
            for a in x
                check_indent(a, state)
                indent_pass(a, state)
            end
        end
    elseif (headof(x) === :block && headof(x[1]) == :BEGIN) || (headof(x) === :quote && headof(x[1]) === :QUOTE)
        state.offset += x[1].fullspan
        state.edits.indent += 1
        if length(x) > 2
            for i = 2:length(x) - 1
                a = x[i]
                check_indent(a, state)
                indent_pass(a, state)
            end
        end
        state.edits.indent -= 1
        check_indent(x[3], state)
        state.offset += x[3].fullspan
    elseif headof(x) == :struct
        state.offset += x[1].fullspan + x.args[2].fullspan
        if length(x.trivia) > 2 #mutable
            state.offset += x.trivia[2].fullspan
        end
        if headof(x.args[3]) === :block && length(x.args[3]) > 0
            state.edits.indent += 1
            for a in x.args[3]
                check_indent(a, state)
                indent_pass(a, state)
            end
            state.edits.indent -= 1
            check_indent(last(x.trivia), state)
            state.offset += last(x.trivia).fullspan
        end
    elseif headof(x) in (:function, :macro, :for, :while)
        state.offset += x[1].fullspan + x[2].fullspan
        if headof(x[3]) === :block
            state.edits.indent += 1
            if length(x[3]) > 0
                for a in x[3]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.edits.indent -= 1
            check_indent(x[4], state)
            state.offset += x[4].fullspan
        else
            check_indent(x[3], state)
            state.offset += x[3].fullspan
        end
    elseif headof(x) === :do
        state.offset += x[1].fullspan + x[2].fullspan + x[3][1].fullspan
        if headof(x[3][3]) === :block
            state.edits.indent += 1
            if length(x[3][3]) > 0
                for a in x[3][3]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.edits.indent -= 1
            check_indent(x[4], state)
            state.offset += x[4].fullspan
        else
            check_indent(x[4], state)
            state.offset += x[4].fullspan
        end
    elseif headof(x) === :macrocall
        if headof(x[1]) === :globalrefdoc
            state.offset += x[1].fullspan
            doc = x.args[3]
            doc_strs = split(str_value(doc), "\n")
            state.offset += 4
            for (i, s) in enumerate(doc_strs)
                # Skip indenting lines of "".
                # The final "" is associated with identing the
                # trailing docstring triple quote
                if s == "" && i != length(doc_strs)
                    state.offset += 1
                else
                    a = EXPR(:STRING, length(s) + 1, length(s), String(s))
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.offset += 3

            check_indent(x.args[4], state)
            indent_pass(x.args[4], state)
        elseif length(x) > 0
            for a in x
                indent_pass(a, state)
            end
        end
    elseif headof(x) === :mutable
        state.offset += x[1].fullspan + x[2].fullspan + x[3].fullspan
        if headof(x[4]) === :block
            state.edits.indent += 1
            if length(x[4]) > 0
                for a in x[4]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.edits.indent -= 1
            check_indent(x[5], state)
            state.offset += x[5].fullspan
        else
            check_indent(x[3], state)
            state.offset += x[4].fullspan
        end
    elseif headof(x) === :try && length(x) > 5
        state.offset += x[1].fullspan
        state.edits.indent += 1
        if length(x[2]) > 0
            for a in x[2]
                check_indent(a, state)
                indent_pass(a, state)
            end
        end
        state.edits.indent -= 1
        check_indent(x[3], state)
        state.offset += x[3].fullspan + x[4].fullspan

        state.edits.indent += 1
        if length(x[5]) > 0
            for a in x[5]
                check_indent(a, state)
                indent_pass(a, state)
            end
        end
        state.edits.indent -= 1
        check_indent(x[6], state)
        state.offset += x[6].fullspan

        if length(x) == 8
            state.edits.indent += 1
            if length(x[7]) > 0
                for a in x[7]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.edits.indent -= 1
            check_indent(x[8], state)
            state.offset += x[8].fullspan
        end
    elseif headof(x) === :if
        if headof(first(x)) == :IF
            state.offset += x[1].fullspan + x[2].fullspan
            state.edits.indent += 1
            if length(x[3]) > 0
                for a in x[3]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            if length(x) >= 4
                state.edits.indent -= 1
                check_indent(x[4], state)
                state.offset += x[4].fullspan
            end
            if length(x) == 5
                # state.edits.indent += 1
                check_indent(x[5], state)
                state.offset += x[5].fullspan
            elseif length(x) == 6
                state.edits.indent += 1
                for a in x[5]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
                state.edits.indent -= 1
                check_indent(x[6], state)
                state.offset += x[6].fullspan
            end
        else
            state.offset += x[1].fullspan
            if length(x[2]) > 0
                for a in x[2]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            if length(x) > 2
                state.edits.indent -= 1
                check_indent(x[3], state)
                state.offset += x[3].fullspan
                state.edits.indent += 1
                if length(x[4]) > 0
                    for a in x[4]
                        check_indent(a, state)
                        indent_pass(a, state)
                    end
                end
            end
        end
    elseif headof(x) === :let
        if length(x) > 3
            state.offset += x[1].fullspan + x[2].fullspan
            state.edits.indent += 1
            if length(x[3]) > 0
                for a in x[3]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.edits.indent -= 1
            check_indent(x[4], state)
            state.offset += x[4].fullspan
        else
            state.offset += x[1].fullspan
            state.edits.indent += 1
            if length(x[2]) > 0
                for a in x[2]
                    check_indent(a, state)
                    indent_pass(a, state)
                end
            end
            state.edits.indent -= 1
            check_indent(x[3], state)
            state.offset += x[3].fullspan
        end

    elseif x.args === nothing
        state.offset += x.fullspan
    else
        if length(x) > 0
            for a in x
                indent_pass(a, state)
            end
        end
    end
    state
end

function check_indent(x, state)
    for (l, i) in state.lines
        if state.offset == l + i
            if state.edits.indent * INDENT != i
                #= @info CSTParser.str_value(CSTParser.get_name(x)), state.edits.indent*INDENT, i, state.offset =#
                push!(state.edits.edits, (l, state.edits.indent * INDENT - i))
            end
        end
    end
end

function indents(text, opts)
    x = CSTParser.parse(text, true)

    state = indent_pass(x, State(0, IndentState(0, []), opts, text, get_lines(text)))

    sort!(state.edits.edits, lt=(a, b) -> a[1] < b[1], rev=true)
    for (l, d) in state.edits.edits
        if d > 0
            text = string(text[1:l], " "^d, text[l + 1:end])
        else
            text = string(text[1:l], text[l + 1 - d:end])
        end
    end
    return text
end


text:
#
# base.jl -
#
# Implements basic and low level methods for Julia interface to ActiveSilicon
# Phoenix (PHX) library.
#
#------------------------------------------------------------------------------
#
# This file is part of the `Phoenix.jl` package which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut (https://github.com/emmt/Phoenix.jl).
# Copyright (C) 2016, Éric Thiébaut & Jonathan Léger.
#

# Override bitwise operators for frame grabber parameters.
import Base: |, &, ~
(~)(x::Param{T,A}) where {T,A} = Param{T,A}(~x.ident)
(|)(x::Param{T,A}, y::Integer) where {T,A} = Param{T,A}(x.ident | y)
(&)(x::Param{T,A}, y::Integer) where {T,A} = Param{T,A}(x.ident & y)
Base.xor(x::Param{T,A}, y::Integer) where {T,A} = Param{T,A}(xor(x.ident, y))

# Conversion to integer values of frame grabber parameters.
convert(::Type{ParamValue}, x::Param) = x.ident
convert(::Type{T}, x::Param) where {T<:Integer} = convert(T, x.ident)

# Override `isreadable` and `iswritable` methods for frame grabber parameters
# and CoaXPress registers.  For the more general case, we do not want to
# pollute the other modules so we left the method definition be local to this
# module.
isreadable(::Any) = false
iswritable(::Any) = false
Base.isreadable(::Type{A}) where {A<:Readable} = true
Base.iswritable(::Type{A}) where {A<:Writable} = true
Base.isreadable(::Param) = false
Base.iswritable(::Param) = false
Base.isreadable(::Param{T,A}) where {T,A<:Readable} = true
Base.iswritable(::Param{T,A}) where {T,A<:Writable} = true
Base.isreadable(::Register) = false
Base.iswritable(::Register) = false
Base.isreadable(::RegisterValue{T,A}) where {T,A<:Readable} = true
Base.iswritable(::RegisterValue{T,A}) where {T,A<:Writable} = true
Base.isreadable(::RegisterAddress{T,A}) where {T,A<:Readable} = true
Base.iswritable(::RegisterAddress{T,A}) where {T,A<:Writable} = true
Base.isreadable(::RegisterString{N,A}) where {N,A<:Readable} = true
Base.iswritable(::RegisterString{N,A}) where {N,A<:Writable} = true

# Make a camera instance usable as an indexable object, which is useful
# to write configuration scripts.

Base.getindex(cam::Camera, param::Param) = getparam(cam, param)
Base.getindex(cam::Camera, reg::Register) = getparam(cam, reg)
Base.getindex(cam::Camera, key) = error("invalid key type `$(typeof(key))`")

Base.setindex!(cam::Camera, val, param::Param) =
    (setparam!(cam, param, val); return cam)
Base.setindex!(cam::Camera, val, reg::Register) =
    (setparam!(cam, reg, val); return cam)
Base.setindex!(cam::Camera, val, key) =
    error("invalid key type `$(typeof(key))`")

# Other overrides.
Base.eltype(cam::Camera) = eltype(eltype(cam.bufs))
Base.length(cam::Camera) = length(cam.bufs)
Base.getindex(cam::Camera, i::Integer) = getindex(cam.bufs, i)

@inline checkstatus(status::Status) =
    (status == PHX_OK || throw(PHXError(status)); nothing)

#-------------------------------------------------------------------------------
# Get/set Frame Grabber Parameters
# ================================

"""
    getparam(cam, key) -> val

yields the value `val` associated with `key` for the camera `cam`.  Argument
`key` can be a frame grabber parameter or a CoaXPress register (which must be
readable, *i.e.* `isreadable(key)` must be true).  Beware that the camera board
must be open before retrieving parameters (otherwise you get a
`PHX_ERROR_BAD_HANDLE` error).

This method implements the syntax:

    cam[key] -> val

A variant is available which does not throw an error but returns a status and
stores the value in `buf` (a pointer or a reference):

    _getparam(cam, key, buf) -> status

See also: [`setparam!`](@ref), [`isreadable`](@ref).

"""
function getparam(cam::Camera,
                  key::Param{T,A}) :: T where {T<:Integer,A<:Readable}
    buf = Ref{T}()
    checkstatus(_getparam(cam, key, buf))
    return buf[]
end

function getparam(cam::Camera,
                  key::Param{String,A}) :: String where {A<:Readable}
    buf = Ref{Ptr{UInt8}}()
    checkstatus(_getparam(cam, key, buf))
    return (buf[] == C_NULL ? "" : unsafe_string(buf[]))
end

function getparam(cam::Camera,
                  key::Param{Ptr{Cvoid},A}) :: Ptr{Cvoid} where {A<:Readable}
    buf = Ref{Ptr{Cvoid}}()
    checkstatus(_getparam(cam, key, buf))
    return buf[]
end

function getparam(cam::Camera,
                  key::RegisterString{N,A}) :: String where {N,A<:Readable}
    buf = Vector{UInt8}(undef, N)
    checkstatus(_readregister(cam, key, buf, N))
    return String(buf)
end

function getparam(cam::Camera,
                  key::RegisterValue{T,A}) :: T where {T<:Real,A<:Readable}
    buf = Ref{T}()
    checkstatus(_getparam(cam, key, buf))
    return buf[]
end

function getparam(cam::Camera,
                  key::RegisterAddress{T,A}) :: T where {T,A<:Readable}
    return getparam(cam, resolve(cam, key))
end

getparam(cam::Camera, key::Param) =
    error("unreadable parameter or undetermined parameter type")

getparam(cam::Camera, key::Register) =
    error("attempt to get an unreadable CoaXPress parameter")

function _getparam(cam::Camera,
                   key::Param{T,A},
                   buf::Union{Ptr,Ref}) where {T,A<:Readable}
    return ccall(_PHX_ParameterGet[], Status, (Handle, Cuint, Ptr{Cvoid}),
                 cam.handle, key.ident, buf)
end

function _getparam(cam::Camera,
                   key::RegisterValue{T,A},
                   buf::Ref{T}) where {T<:Real,A<:Readable}
    status = _readregister(cam, key, buf, sizeof(T))
    if cam.swap
        buf[] = bswap(buf[])
    end
    return status
end

@doc @doc(getparam) _getparam

"""

    resolve(cam, regaddr) -> regval

yields the CoaXPress register at indirect register address `regaddr` for camera
`cam`.

See also: [`getparam`](@ref), [`setparam!`](@ref).

"""
function resolve(cam::Camera,
                 reg::RegisterAddress{T,A}
                 ) :: RegisterValue{T,A} where {T,A<:AccessMode}
    buf = Ref{UInt32}()
    checkstatus(_readregister(cam, reg, buf, 4))
    addr = (cam.swap ? bswap(buf[]) : buf[])
    return RegisterValue{T,A}(addr)
end

"""
    setparam!(cam, key, val)

set the value associated with `key` to be `val` for the camera `cam`.  Argument
`key` can be a frame grabber parameter or a CoaXPress register (which must be
writable, *i.e.* `iswritable(key)` must be true).  This method implements the
syntax:

    cam[key] = val

A variant is available which does not throw an error but returns a status and
stores the value in `buf` (a pointer or a reference):

    _setparam(cam, key, buf) -> status

Beware that this low-level version does not check its arguments.

See also: [`getparam`](@ref), [`iswritable`](@ref).

"""
function setparam!(cam::Camera,
                   key::Param{T,A},
                   val::Integer) where {T<:Integer,A<:Writable}
    checkstatus(_setparam!(cam, key, Ref{T}(val)))
    return val
end

function setparam!(cam::Camera,
                   key::Param{String,A},
                   str::AbstractString) where {A<:Writable}
    checkstatus(_setparam!(cam, key, Ref(pointer(cstring(str)))))
    return str
end

function setparam!(cam::Camera,
                   key::Param{Nothing,A},
                   ::Nothing) where {A<:Writable}
    checkstatus(_setparam!(cam, key, C_NULL))
    return nothing
end

function setparam!(cam::Camera,
                   key::Param{Ptr{T},A},
                   buf::Union{Vector{T},Ptr{T},Ref{T}}) where {T,A<:Writable}
    checkstatus(_setparam!(cam, key, buf))
    return buf
end

function setparam!(cam::Camera,
                   key::Param{Ptr{Cvoid},A},
                   buf::Union{Vector,Ptr,Ref}) where {A<:Writable}
    checkstatus(_setparam!(cam, key, buf))
    return buf
end

function setparam!(cam::Camera,
                   key::RegisterString{N,A},
                   str::AbstractString) where {N,A<:Writable}
    buf = Array{UInt8}(undef, N)
    i = 0
    @inbounds for c in str
        if i ≥ N
            break
        end
        c != '\0' || error("string must not have embedded NUL characters")
        i += 1
        buf[i] = c
    end
    @inbounds while i < N
        i += 1
        buf[i] = zero(UInt8)
    end
    checkstatus(_writeregister(cam, key, buf, N))
    return str
end

function setparam!(cam::Camera,
                   key::RegisterValue{T,A},
                   val::Real) where {T<:Real,A<:Writable}
    setparam!(cam, key, convert(T, val))
    return val
end

function setparam!(cam::Camera,
                   key::RegisterValue{T,A},
                   val::T) where {T<:Real,A<:Writable}
    checkstatus(_setparam!(cam, key, Ref{T}(val)))
    return val
end

function setparam!(cam::Camera,
                   key::RegisterAddress{T,A},
                   val) where {T,A<:Writable}
    return setparam!(cam, resolve(cam, key), val)
end

setparam!(cam::Camera, key::Param, val) =
    error("unwritable parameter or undetermined parameter type")

setparam!(cam::Camera, key::Register, val) =
    error("attempt to set an unwritable CoaXPress parameter")

function _setparam!(cam::Camera,
                    key::Param{T,A},
                    buf::Union{Ptr,Ref,Vector}) where {T,A<:Writable}
    return ccall(_PHX_ParameterSet[], Status, (Handle, Cuint, Ptr{Cvoid}),
                 cam.handle, key.ident, buf)
end

function _setparam!(cam::Camera,
                    key::RegisterValue{T,A},
                    buf::Ref{T}) where {T<:Real,A<:Writable}
    if cam.swap
        buf[] = bswap(buf[])
    end
    return _writeregister(cam, key, buf, sizeof(T))
end

@doc @doc(setparam!) _setparam!

"""
    flushcache(cam)

makes sure that cached frame grabber parameters are effectively written to the
hardware.

"""
flushcache(cam::Camera) =
    setparam!(cam, PHX_DUMMY_PARAM|PHX_CACHE_FLUSH, nothing)

#------------------------------------------------------------------------------
# Configuration
# =============
#

CameraConfig(cam::Camera) = getconfig(cam)

"""
```julia
checkconfig(cam, cfg) -> boolean
```

yields whether the settings in configuration `cfg` are suitable for the camera
`cam`.  If keyword `throwerrors` is true, an exception is thrown if the
configuration has some invalid settings.

See also [`getconfig`](@ref), [`loadconfig`](@ref), [`saveconfig`](@ref) and
[`updateconfig!`](@ref).

"""
checkconfig(cam::Camera{M,C}, cfg::C; throwerrors::Bool=false) where {M,C} =
    error("not implemented")
function checkconfig(cam::Camera, cfg::CameraConfig; throwerrors::Bool=false)
    throwerrors && error("incompatible configuration type and camera model")
    return false
end

# Timings: getconfig  ~ 13 ns
#          getconfig! ~ 5.7 ns

"""
```julia
getconfig(cam) -> cfg
```

yields a copy of the current configuration of camera `cam`.

See also [`getconfig!`](@ref), [`setconfig`](@ref), [`loadconfig`](@ref),
[`saveconfig`](@ref) and [`updateconfig!`](@ref).

"""
getconfig(cam::Camera{M,C}) where {M,C} = getconfig!(C(), cam)

"""
```julia
getconfig!(cfg, cam) -> cfg
```

overwrites `cfg` with the current configuration of camera `cam` and returns
`cfg`.

 See also [`getconfig`](@ref), [`setconfig`](@ref), [`loadconfig`](@ref),
[`saveconfig`](@ref) and [`updateconfig!`](@ref).

"""
getconfig!(cfg::C, cam::Camera{M,C}) where {M<:CameraModel,C<:CameraConfig} =
    _fastcopy!(cfg, cam.config)
getconfig!(cfg::CameraConfig, cam::Camera) =
    error("incompatible configuration type and camera model")

# Private method to copy the contents of mutable data.  Beware that its use
# must be restricted to *simple* mutable structures which do not hold
# references.
function _fastcopy!(dst::T, src::T) where {T}
    # @assert isconcretetype(T) # not needed, sizeof() will complain if T is an
    #                           # abstract type
    ccall(:memcpy, Ptr{Cvoid}, (Ref{T}, Ref{T}, Csize_t), dst, src, sizeof(T))
    return dst
end

"""
```julia
setconfig!(cam, cfg) -> cam
```

applies the settings in configuration `cfg` to the camera `cam` and returns
`cam`.

See also [`getconfig`](@ref), [`loadconfig`](@ref), [`saveconfig`](@ref) and
[`updateconfig!`](@ref).

"""
setconfig!(cam::Camera{M,C}, cfg::C) where {M,C} =
    error("not implemented")
setconfig!(cam::Camera, cfg::CameraConfig) =
    error("incompatible configuration type and camera model")

"""
    saveconfig(cam, name, what = PHX_SAVE_ALL)

saves the actual configuration of camera `cam` in file `name`.  Optional
argument `what` specifies which parameters to save, it can be a combination
(*i.e.* bitwise or) of:

- `PHX_SAVE_CAM` to save the camera specific parameters.  These describe the
  camera.

- `PHX_SAVE_SYS` to save the system specific parameters.  These describe how
  the camera is connected to the Active Silicon board.

- `PHX_SAVE_APP` to save the application specific parameters.

- `PHX_SAVE_ALL` to save all three of the above types of parameters.

See also [`getconfig`](@ref), [`setconfig!`](@ref), [`loadconfig`](@ref) and
[`updateconfig!`](@ref).

"""
function saveconfig(cam::Camera, name::AbstractString,
                    what::Integer = PHX_SAVE_ALL)
    flushcache(cam)
    status = ccall(_PHX_Action[], Status,
                   (Handle, Action, ActionParam, Ptr{Cvoid}),
                   cam.handle, PHX_CONFIG_SAVE, what, cstring(name))
    return checkstatus(status)
end

"""
```julia
updateconfig!(cam; all=false) -> cam
```

updates the current configuration memorized by camera `cam`.  Normally only
parameters that might change whithout being explicitly modified (like the
temperature) are updated.  If keyword `all` is set true, all parameters are
read from the device.  This should not be needed in case if the caching
of parameters is done properly.

See also [`getconfig`](@ref), [`setconfig!`](@ref), [`loadconfig`](@ref) and
[`saveconfig`](@ref).

"""
updateconfig!(cam::Camera; all::Bool=false) =
    error("not implemented for this kind of camera")

#------------------------------------------------------------------------------
# Reading/Writing CoaXPress Registers
# ===================================
#

"""
`exec(cam, cmd)` executes CoaXPress command `cmd` for camera `cam`, throwing
exception in case of error.

`_exec(cam, cmd)` does the same but returns a status instead of throwing exceptions

"""
exec(cam::Camera, reg::RegisterCommand{T}) where {T} =
    checkstatus(_exec(cam, reg))

function _exec(cam::Camera, reg::RegisterCommand{T}) where {T}
    data = Ref{T}(cam.swap ? bswap(reg.value) : reg.value)
    return _writeregister(cam, reg, data, sizeof(T))
end

readstream(cam::Camera, args...) = checkstatus(_readstream(cam, args...))

_debug(cam::Camera, args...) = cam.debug && println(args...)

#-------------------------------------------------------------------------------
# Low Level Wrapper to Phoenix Dynamic Library
# ============================================
#

const _PHX_FUNCTIONS = (:_PHX_Create, :_PHX_Open, :_PHX_Close, :_PHX_Destroy,
                        :_PHX_StreamRead, :_PHX_ParameterGet,
                        :_PHX_ParameterSet, :_PHX_ControlRead,
                        :_PHX_ControlWrite, :_PHX_Action,
                        :_PHX_ErrCodeDecode, :_PHX_ErrHandlerDefault)
for sym in _PHX_FUNCTIONS
    @eval const $sym = Ref{Ptr{Cvoid}}(0)
end

function _errorhandler end
function _callback end

const _ERRORHANDLER_REF = Ref{Ptr{Cvoid}}(0)
const _CALLBACK_REF = Ref{Ptr{Cvoid}}(0)
function __init__()
    _ERRORHANDLER_REF[] = @cfunction(_errorhandler, Cvoid,
                                     (Ptr{Cchar}, Status, Ptr{Cchar}))
    _CALLBACK_REF[] = @cfunction(_callback, Cvoid,
                                 (Handle, UInt32, Ptr{Cvoid}))

    # Manage to load the dynamic library and its symbols with appropriate
    # flags.  It is still needed to start Julia with the correct dynamic library
    # search path (as checked above).
    handle = Libdl.dlopen(_PHXLIB, (Libdl.RTLD_LAZY |
                                    Libdl.RTLD_DEEPBIND |
                                    Libdl.RTLD_GLOBAL))
    for sym in _PHX_FUNCTIONS
        @eval $sym[] = Libdl.dlsym($handle, $(string(sym)[2:end]))
    end
end

"""
    open(model; configfile="",
                boardtype=0,
                boardnumber=PHX_BOARD_AUTO,
                channelnumber=PHX_CHANNEL_AUTO,
                boardmode=PHX_MODE_NORMAL)

yields a new camera instance.

    open(cam) -> cam

opens an ActiveSilicon board, configuring the hardware.  This operation must be
performed after creating the camera and setting some initial parameters.

See also: [`close`](@ref).

"""
Base.open(::Type{M}; kwds...) where {M<:CameraModel} =
    open(Camera(M); kwds...)

function Base.open(cam::Camera;
                   configfile::String = "",
                   boardtype::Integer = 0,
                   boardnumber::Integer = PHX_BOARD_AUTO,
                   channelnumber::Integer = PHX_CHANNEL_AUTO,
                   boardmode::Integer = PHX_MODE_NORMAL,
                   quiet::Bool = false)
    # Check state.
    if cam.state != 0
        1 ≤ cam.state ≤ 2 || error("camera structure corrupted")
        @warn "camera has already been opened"
        return cam
    end

    # Create references for parameter values and camera handle (these
    # references will persist until the final call to `PHX_Open`, i.e. we
    # do not assume that a copy of these settings is immediately made).
    ref_configfile = Ref{Ptr{UInt8}}(length(configfile) > 0 ?
                                     pointer(cstring(configfile)) :
                                     Ptr{UInt8}(0))
    ref_boardtype = Ref{ParamValue}(boardtype)
    ref_boardnumber = Ref{ParamValue}(boardnumber)
    ref_channelnumber = Ref{ParamValue}(channelnumber)
    ref_boardmode = Ref{ParamValue}(boardmode)

    # Set initial parameters.
    if boardtype != 0
        checkstatus(_setparam!(cam, PHX_BOARD_VARIANT, ref_boardtype))
    end
    checkstatus(_setparam!(cam, PHX_BOARD_NUMBER,   ref_boardnumber))
    checkstatus(_setparam!(cam, PHX_CHANNEL_NUMBER, ref_channelnumber))
    checkstatus(_setparam!(cam, PHX_CONFIG_MODE,    ref_boardmode))
    checkstatus(_setparam!(cam, PHX_CONFIG_FILE,    ref_configfile))

    # Open the camera.
    checkstatus(ccall(_PHX_Open[], Status, (Handle,), cam.handle))
    cam.state = 1

    # Discover whether we have a CoaXPress camera.
    coaxpress = ((cam[PHX_CXP_INFO] & PHX_CXP_CAMERA_DISCOVERED) != 0)
    if coaxpress != isa(cam, CoaXPressCamera)
        error("(non-)CoaXPress camera not properly discovered")
    end
    if coaxpress
        # Figure out byte order.
        magic = _read_magic(cam)
        if magic == bswap(CXP_MAGIC)
            cam.swap = ! cam.swap
        elseif magic != CXP_MAGIC
            error("unexpected magic number for CoaXPress camera")
        end

        # Get current width and height and initialize parameters.
        width  = _read_width(cam)
        height = _read_height(cam)
        cam[PHX_CAM_ACTIVE_XOFFSET] = 0
        cam[PHX_CAM_ACTIVE_YOFFSET] = 0
        cam[PHX_CAM_ACTIVE_XLENGTH] = width
        cam[PHX_CAM_ACTIVE_YLENGTH] = height

        if ! quiet
            pixelformat =
            println("  Vendor name:  ", _read_vendor_name(cam))
            println("  Model name:   ", _read_device_model_name(cam))
            println("  Image size:   $width × $height pixels")
            println("  Pixel format: ",_read_pixel_format(cam))
        end
    end

    # Apply specific post-open configuration.
    _openhook(cam)

    return cam
end

"""
    _openhook(cam)

performs specific operations just after camera `cam` has been opened.  This
function should return nothing but may throw exceptions to signal errors.

See also: [`open`](@ref)

"""
_openhook(cam::Camera) = nothing

"""
    close(cam)

closes an Active Silicon board, releasing the hardware.  Any running
acquisition is aborted.  The camera is automatically closed when finalized by
the garbage collector, so calling this method is optional.

See also: [`open`](@ref), [`abort`](@ref).

"""
function Base.close(cam::Camera)
    if cam.state == 2
        # Acquisition is running, abort it.
        abort(cam)
    end
    if cam.state == 0
        @warn "camera has already been closed"
    elseif cam.state == 1
        cam.state = 0 # avoid closing more than once
        status = _close(cam)
        status == PHX_OK || throw(PHXError(status))
    else
        error("camera structure corrupted")
    end
    return nothing
end

Base.isopen(cam::Camera) = (1 ≤ cam.state ≤ 2)

# Call PHX_Close, does no checkings but the validity of assumptions and return status.
function _close(cam::Camera)
    # Note that PHX_Close requires the address of the handle but left its
    # contents unchanged.
    ref = Ref(cam.handle)
    status = ccall(_PHX_Close[], Status, (Ptr{Handle},), ref)
    @assert cam.handle == ref[] # cf. note above
    return status
end

"""
    _readstream(cam, cmd, ptr) -> status

initiates and controls the acquisition of stream data such as images.  The
destination of the data can be system memory, graphics memory, or other user
addressable memory, such as slave PCI cards.  `cam` is a PHX camera, `cmd` is
the command to be performed and `ptr` is a command-specific parameter (ususally
a pointer).

"""
_readstream(cam::Camera, cmd::Acq, ptr::Union{Ref,Ptr}) =
    ccall(_PHX_StreamRead[], Status, (Handle, Acq, Ptr{Cvoid}),
          cam.handle, cmd, ptr)

"""
    _readregister(cam, reg, data, num, timeout = cam.timeout) -> status

reads register `reg` from camera `cam`.

See also: [`_readcontrol`](@ref), [`_writeregister`](@ref).

"""
function _readregister(cam::Camera, reg::Register, data, num::Integer,
                       timeout::Integer = cam.timeout)
    return _readcontrol(cam.handle, PHX_REGISTER_DEVICE,
                        Ref(reg.addr), data, Ref{UInt32}(num), timeout)
end

"""
    _readcontrol(handle, src, param, buf, num, timeout) -> status

This function receives `num[]` bytes from Active Silicon board `handle` via the
control port selected by `src` and `param` and stores them in buffer `buf`.  If
a data error occurs during reception of a byte, or a packet error occurs, or
the timeout period `timeout` (in milliseconds) is exceeded, no further bytes
will be stored in `buf`.  (Note: Depending on the protocol, the underlying
hardware may retry a read before reporting an error).  The actual number of
bytes received within the specified timeout period is stored in `num[]`.  Note
that registers in a CoaXPress device (e.g. camera) are big-endian.  Therefore
byte swapping (e.g. with `bswap`) may be needed when using this function to
access a CoaXPress device.  A standard status is returned.

See also: [`_readregister`](@ref), [`_writecontrol`](@ref).

"""
function _readcontrol(handle::Handle, src::ControlPort, param, buf,
                      num::Ref{UInt32}, timeout::Integer)
    return ccall(_PHX_ControlRead[], Status,
                 (Handle, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt32}, UInt32),
                 handle, src.port, param, buf, num, timeout)
end

"""
    _writeregister(cam, reg, data, num, timeout = cam.timeout) -> status

writes register `reg` to camera `cam`.

See also: [`_writecontrol`](@ref), [`_readregister`](@ref).

"""
function _writeregister(cam::Camera, reg::Register, data, num::Integer,
                        timeout::Integer = cam.timeout)
   return _writecontrol(cam.handle, PHX_REGISTER_DEVICE,
                        Ref(reg.addr), data, Ref{UInt32}(num), timeout)
end

"""
    _writecontrol(handle, src, param, buf, num, timeout) -> status

This function sends `num[]` bytes from buffer `buf` to Active Silicon board
`handle` via the control port selected by `src` and `param`.  If a data error
occurs during transmission of a byte, or a packet error occurs, or the timeout
period `timeout` (in milliseconds) is exceeded, no further bytes will be
transmitted.  (Note: Depending on the protocol, the underlying hardware may
retry a write before reporting an error).  The actual number of bytes
transmitted within the specified timeout period is stored in `num[]`.  Note
that registers in a CoaXPress device (e.g. camera) are big-endian.  Therefore
byte swapping (e.g. with `bswap`) may be needed when using this function to
access a CoaXPress device.  A standard status is returned.

See also: [`_writeregister`](@ref), [`_readcontrol`](@ref).

"""
function _writecontrol(handle::Handle, src::ControlPort, param, buf,
                       num::Ref{UInt32}, timeout::Integer)
    return ccall(_PHX_ControlWrite[], Status,
                 (Handle, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt32}, UInt32),
                 handle, src.port, param, buf, num, timeout)
end