text:
# MAT.jl
# Tools for reading MATLAB v5 files in Julia
#
# Copyright (C) 2012   Timothy E. Holy and Simon Kornblith
#
# 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.

VERSION >= v"0.4.0-dev+6521" && __precompile__()

module MAT

using HDF5, Compat

include("MAT_HDF5.jl")
include("MAT_v5.jl")

using .MAT_HDF5, .MAT_v5

export matopen, matread, matwrite, names, exists, @read, @write

# Open a MATLAB file
const HDF5_HEADER = UInt8[0x89, 0x48, 0x44, 0x46, 0x0d, 0x0a, 0x1a, 0x0a]
function matopen(filename::AbstractString, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool)
    # When creating new files, create as HDF5 by default
    fs = filesize(filename)
    if cr && (tr || fs == 0)
        return MAT_HDF5.matopen(filename, rd, wr, cr, tr, ff)
    elseif fs == 0
        error("File \"$filename\" does not exist and create was not specified")
    end

    # Test whether this is a MAT file
    if fs < 128
        error("File \"$filename\" is too small to be a supported MAT file")
    end
    rawfid = open(filename, "r")

    # Check for MAT v4 file
    magic = read!(rawfid, Vector{UInt8}(4))
    for i = 1:length(magic)
        if magic[i] == 0
            close(rawfid)
            error("\"$filename\" is not a MAT file, or is an unsupported (v4) MAT file")
        end
    end

    # Check for MAT v5 file
    seek(rawfid, 124)
    version = read(rawfid, UInt16)
    endian_indicator = read(rawfid, UInt16)
    if (version == 0x0100 && endian_indicator == 0x4D49) ||
       (version == 0x0001 && endian_indicator == 0x494D)
        if wr || cr || tr || ff
            error("creating or appending to MATLAB v5 files is not supported")
        end
        return MAT_v5.matopen(rawfid, endian_indicator)
    end

    # Check for HDF5 file
    for offset = 512:512:fs-8
        seek(rawfid, offset)
        if read!(rawfid, Vector{UInt8}(8)) == HDF5_HEADER
            close(rawfid)
            return MAT_HDF5.matopen(filename, rd, wr, cr, tr, ff)
        end
    end

    close(rawfid)
    error("\"$filename\" is not a MAT file")
end

function matopen(fname::AbstractString, mode::AbstractString)
    mode == "r"  ? matopen(fname, true , false, false, false, false) :
    mode == "r+" ? matopen(fname, true , true , false, false, false) :
    mode == "w"  ? matopen(fname, false, true , true , true , false) :
#     mode == "w+" ? matopen(fname, true , true , true , true , false) :
#     mode == "a"  ? matopen(fname, false, true , true , false, true ) :
#     mode == "a+" ? matopen(fname, true , true , true , false, true ) :
    error("invalid open mode: ", mode)
end

matopen(fname::AbstractString) = matopen(fname, "r")

function matopen(f::Function, args...)
    fid = matopen(args...)
    try
        f(fid)
    finally
        close(fid)
    end
end

"""
    matopen(filename [, mode]) -> handle
    matopen(f::Function, filename [, mode]) -> f(handle)

Mode defaults to "r" for read.  It can also be "w" for write, or "r+" for
read or write without creation or truncation.

Use with `read`, `write`, `close`, `names`, and `exists`.
"""
matopen

# Read all variables from a MATLAB file
"""
    matread(filename) -> Dict

Return a dictionary of all the variables and values in a Matlab file,
opening and closing it automatically.
"""
function matread(filename::AbstractString)
    file = matopen(filename)
    local vars
    try
        vars = read(file)
    finally
        close(file)
    end
    vars
end

# Write a dict to a MATLAB file
"""
    matwrite(filename, d::Dict)

Write a dictionary containing variable names as keys and values as values
to a Matlab file, opening and closing it automatically.
"""
function matwrite{S, T}(filename::AbstractString, dict::Associative{S, T})
    file = matopen(filename, "w")
    try
        for (k, v) in dict
            local kstring
            try
                kstring = ascii(convert(String, k))
            catch x
                error("matwrite requires a Dict with ASCII keys")
            end
            write(file, kstring, v)
        end
    finally
        close(file)
    end
end
end


text:
module MonteCarlo

using LinearAlgebra: AbstractMatrix
using Reexport
# Loading the RNG will fail if Random is nto exported
@reexport using MonteCarloObservable, Random
import MonteCarloObservable.AbstractObservable
using Parameters, Requires
using TimerOutputs, StructArrays
using Printf, SparseArrays, LinearAlgebra, Dates, Statistics, Random, Distributed
import ProgressMeter

if get(ENV, "MONTECARLO_USE_LOOPVECTORIZATION", "true") == "true"
    import LoopVectorization
    using LoopVectorization: @turbo
else
    printstyled(
        "Using MonteCarlo.jl without LoopVectorization. This should only be done for tests.",
        color = :red
    )
    macro turbo(code)
        esc(quote @inbounds @fastmath $code end)
    end
end

import JLD, JLD2
using CodecZlib

# Because we fully load all data directly for JLD we lose access to the path
# This is supposed to keep track of path information
struct FileWrapper{T}
    file::T
    path::String

    FileWrapper(file::T, path) where T = new{T}(file, abspath(path))
end

function Base.getindex(fw::FileWrapper{T}, k) where T
    out = getindex(fw.file, k)
    if out isa Union{JLD2.JLDFile, JLD2.Group, Dict} # JLD generates a nested dict :(
        return FileWrapper(out, fw.path)
    else
        out
    end
end
Base.setindex!(fw::FileWrapper, k, v) = setindex!(fw.file, k, v)
Base.haskey(fw::FileWrapper, k) = haskey(fw.file, k)
Base.write(fw::FileWrapper, x) = write(fw.file, x)
Base.write(fw::FileWrapper, k, x) = write(fw.file, k, x)
Base.close(fw::FileWrapper) = close(fw.file)
Base.get(fw::FileWrapper, k, default) = haskey(fw, k) ? fw[k] : default
Base.keys(fw::FileWrapper) = keys(fw.file)

# To allow switching between JLD and JLD2:
const UnknownType = Union{JLD.UnsupportedType, JLD2.UnknownType}
const JLDFile = Union{FileWrapper{<: JLD.JldFile}, FileWrapper{<: JLD2.JLDFile}, JLD.JldFile, JLD2.JLDFile}

include("helpers.jl")
export enable_benchmarks, disable_benchmarks, print_timer, reset_timer!
include("flavors/abstract.jl")
include("models/abstract.jl")
include("lattices/abstract.jl")

include("configurations.jl")
export Discarder, ConfigRecorder, BufferedConfigRecorder, RelativePath, AbsolutePath
include("Measurements.jl")
export measurements, observables

include("lattices/lattice_iterators.jl")
include("lattices/square.jl")
include("lattices/chain.jl")
include("lattices/cubic.jl")
include("lattices/honeycomb.jl")
include("lattices/triangular.jl")
include("lattices/ALPS.jl")
include("lattices/deprecated.jl")
# export directions, RawMask, DistanceMask # maybe getorder?
export AbstractLattice, Chain, SquareLattice, CubicLattice, TriangularLattice, ALPSLattice
export EachSite, EachSiteAndFlavor, OnSite, EachSitePair, EachSitePairByDistance,
        EachLocalQuadByDistance, EachLocalQuadBySyncedDistance,
        Sum, ApplySymmetries
export neighbors, directions

include("flavors/MC/MC.jl")
include("flavors/DQMC/main.jl")
export Greens, GreensAt, CombinedGreensIterator
export GreensMatrix, swapop
export boson_energy_measurement, greens_measurement, occupation, magnetization
export charge_density, charge_density_correlation, charge_density_susceptibility
export spin_density, spin_density_correlation, spin_density_susceptibility
export pairing, pairing_correlation, pairing_susceptibility
export current_current_susceptibility, superfluid_density
export noninteracting_energy, interacting_energy, total_energy

export EmptyScheduler, SimpleScheduler, AdaptiveScheduler
export Adaptive, NoUpdate, LocalSweep
export GlobalFlip, GlobalShuffle, SpatialShuffle, TemporalShuffle,
        Denoise, DenoiseFlip, StaggeredDenoise
export ReplicaExchange, ReplicaPull, connect, disconnect
# export mask, uniform_fourier, structure_factor, SymmetryWrapped, swave, eswave

include("models/Ising/IsingModel.jl")
include("models/HubbardModel.jl")
export IsingEnergyMeasurement, IsingMagnetizationMeasurement

include("FileIO.jl")
export save, load, resume!

export reset!
export run!, resume!, replay!
export Model, IsingModel
# export RepulsiveGHQHubbardModel, AttractiveGHQHubbardModel
export HubbardModel, HubbardModelAttractive, HubbardModelRepulsive
export MonteCarloFlavor, MC, DQMC
export greens, greens!, lattice, model, parameters
export DensityHirschField, MagneticHirschField, MagneticGHQField

# For extending
export AbstractMeasurement, Model

function __init__()
    # @require LatPhysBase="eec5c15a-e8bd-11e8-0d23-6799ca40c963" include("lattices/LatPhys.jl")
    @require LatticePhysics = "53011200-ee7a-11e8-39f1-5f3e57afe4fd" include("lattices/LatPhys.jl")
    @require LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" include("lattices/ALPS.jl")
    @require MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" begin
        include("mpi.jl")
        export mpi_queue

        include("flavors/DQMC/updates/mpi_updates.jl")
        export MPIReplicaExchange, MPIReplicaPull
    end
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include("DataFrames.jl")
end

end # module

ext: jl


text:
using Parameters

@with_kw struct AD_Para
    # mesh
    Nintervals::Int = 50        # number of intervals for your mesh
    basisorder::Int = 2         # element basis function order (1 or 2)
    xmin::Float64 = -1.0        # minimum x-value for mesh
    xmax::Float64 = 1.0         # minimum x-value for mesh
    ymin::Float64 = -1.0        # minimum x-value for mesh
    ymax::Float64 = 1.0         # minimum x-value for mesh

    # boundary conditions
    dirVal::Float64 = 1.0       # (dirichlet) value at left boundary
    Pe::Float64 = 20.0           # Peclet number (advection / diffusion)

    # visualization
    solFolder::String = "solutions"         # folder for solution to go in
    solFile::String = "ad_example"     # name of vtk solutions output
    solName::String = "chemical"         # name of unknown
end

max_stars_repo_name: pseastham/eFEMpart


text:
using Documenter

makedocs(
	sitename = "MimiPAGE2009.jl",
	pages = [
		"Home" => "index.md",
		"Getting started" => "gettingstarted.md",
		"Model Structure" => "model-structure.md",
		"Technical User Guide" => "technicaluserguide.md",
		"Model Validation" => "validation.md"]
)

deploydocs(
    repo = "github.com/anthofflab/MimiPAGE2009.jl.git"
)


text:
struct SubProblem{F <: AbstractFeasibility, T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    id::Int
    π::T

    solver::S
    feasibility_solver::S

    h::Tuple{A,A}
    x::A
    y::A
    masterterms::Vector{Tuple{Int,Int,T}}

    function SubProblem(model::JuMP.Model,
                        id::Integer,
                        π::AbstractFloat,
                        x::AbstractVector,
                        y₀::AbstractVector,
                        masterterms::Vector{Tuple{Int,Int,R}},
                        optimsolver::MPB.AbstractMathProgSolver,
                        ::Type{F}) where {R <: AbstractFloat, F <: AbstractFeasibility}
        T = promote_type(eltype(x), eltype(y₀), R, Float32)
        x_ = convert(AbstractVector{T}, x)
        y₀_ = convert(AbstractVector{T}, y₀)
        masterterms_ = convert(Vector{Tuple{Int,Int,T}}, masterterms)
        A = typeof(x_)
        solver = LQSolver(model, optimsolver)
        subproblem = new{F,T,A,typeof(solver)}(id,
                                               π,
                                               solver,
                                               FeasibilitySolver(model, optimsolver, F),
                                               (convert(A, MPB.getconstrLB(solver.lqmodel)),
                                                convert(A, MPB.getconstrUB(solver.lqmodel))),
                                               x_,
                                               y₀_,
                                               masterterms_)
        return subproblem
    end

    function SubProblem(model::JuMP.Model,
                        parent::JuMP.Model,
                        id::Integer,
                        π::AbstractFloat,
                        x::AbstractVector,
                        y₀::AbstractVector,
                        optimsolver::MPB.AbstractMathProgSolver,
                        ::Type{F}) where F <: AbstractFeasibility
        T = promote_type(eltype(x), eltype(y₀), Float32)
        x_ = convert(AbstractVector{T}, x)
        y₀_ = convert(AbstractVector{T}, y₀)
        A = typeof(x_)

        solver = LQSolver(model, optimsolver)

        subproblem = new{F,T,A,typeof(solver)}(id,
                                               π,
                                               solver,
                                               FeasibilitySolver(model, optimsolver, F),
                                               (convert(A, MPB.getconstrLB(solver.lqmodel)),
                                                convert(A, MPB.getconstrUB(solver.lqmodel))),
                                               x_,
                                               y₀_,
                                               Vector{Tuple{Int,Int,T}}())
        parse_subproblem!(subproblem, model, parent)
        return subproblem
    end
end

function FeasibilitySolver(model::JuMP.Model, optimsolver::MPB.AbstractMathProgSolver, ::Type{IgnoreFeasibility})
    return LQSolver(model, optimsolver; load = false)
end

function FeasibilitySolver(model::JuMP.Model, optimsolver::MPB.AbstractMathProgSolver, ::Type{<:HandleFeasibility})
    solver = LQSolver(model, optimsolver)
    feasibility_problem!(solver)
    return solver
end

function parse_subproblem!(subproblem::SubProblem, model::JuMP.Model, parent::JuMP.Model)
    for (i, constr) in enumerate(model.linconstr)
        for (j, var) in enumerate(constr.terms.vars)
            if var.m == parent
                # var is a first stage variable
                push!(subproblem.masterterms, (i,var.col,-constr.terms.coeffs[j]))
            end
        end
    end
end

function update_subproblem!(subproblem::SubProblem, x::AbstractVector)
    lb = MPB.getconstrLB(subproblem.solver.lqmodel)
    ub = MPB.getconstrUB(subproblem.solver.lqmodel)
    for i in [term[1] for term in unique(term -> term[1], subproblem.masterterms)]
        lb[i] = subproblem.h[1][i]
        ub[i] = subproblem.h[2][i]
    end
    for (i,j,coeff) in subproblem.masterterms
        lb[i] += coeff*x[j]
        ub[i] += coeff*x[j]
    end
    MPB.setconstrLB!(subproblem.solver.lqmodel, lb)
    MPB.setconstrUB!(subproblem.solver.lqmodel, ub)
    update_feasibility_solver!(subproblem, lb, ub)
    subproblem.x[:] = x
    return nothing
end
function update_feasibility_solver!(::SubProblem{IgnoreFeasibility}, ::AbstractVector, ::AbstractVector)
    return nothing
end
function update_feasibility_solver!(subproblem::SubProblem{<:HandleFeasibility}, lb::AbstractVector, ub::AbstractVector)
    MPB.setconstrLB!(subproblem.feasibility_solver.lqmodel, lb)
    MPB.setconstrUB!(subproblem.feasibility_solver.lqmodel, ub)
end
update_subproblems!(subproblems::Vector{<:SubProblem}, x::AbstractVector) = map(prob -> update_subproblem!(prob,x), subproblems)

function get_solution(subproblem::SubProblem)
    return copy(subproblem.y), getredcosts(subproblem.solver), getduals(subproblem.solver), getobjval(subproblem.solver)
end

function solve(subproblem::SubProblem)
    subproblem.solver(subproblem.y)
    solvestatus = status(subproblem.solver)
    if solvestatus == :Optimal
        subproblem.y[:] = getsolution(subproblem.solver)
        return OptimalityCut(subproblem)
    elseif solvestatus == :Infeasible
        return Infeasible(subproblem)
    elseif solvestatus == :Unbounded
        return Unbounded(subproblem)
    else
        error(@sprintf("Subproblem %d was not solved properly, returned status code: %s", subproblem.id,string(solvestatus)))
    end
end

function (subproblem::SubProblem{<:HandleFeasibility})()
    subproblem.feasibility_solver(vcat(subproblem.y, rand(2*MPB.numconstr(subproblem.feasibility_solver.lqmodel))))
    w = getobjval(subproblem.feasibility_solver)
    if w > 0
        return FeasibilityCut(subproblem)
    end
    return solve(subproblem)
end
function (subproblem::SubProblem{IgnoreFeasibility})()
    return solve(subproblem)
end

function (subproblem::SubProblem)(x::AbstractVector)
    update_subproblem!(subproblem, x)
    subproblem.solver(subproblem.y)
    solvestatus = status(subproblem.solver)
    if solvestatus == :Optimal
        subproblem.y[:] = getsolution(subproblem.solver)
        return getobjval(subproblem.solver)
    elseif solvestatus == :Infeasible
        return Inf
    elseif solvestatus == :Unbounded
        return -Inf
    else
        error(@sprintf("Subproblem %d was not solved properly, returned status code: %s", subproblem.id, string(solvestatus)))
    end
end

function OptimalityCut(subproblem::SubProblem)
    λ = getduals(subproblem.solver)
    π = subproblem.π
    cols = zeros(length(subproblem.masterterms))
    vals = zeros(length(subproblem.masterterms))
    for (s,(i,j,coeff)) in enumerate(subproblem.masterterms)
        cols[s] = j
        vals[s] = -π*λ[i]*coeff
    end
    δQ = sparsevec(cols, vals, length(subproblem.x))
    q = π*getobjval(subproblem.solver)+δQ⋅subproblem.x

    return OptimalityCut(δQ, q, subproblem.id)
end

function FeasibilityCut(subproblem::SubProblem)
    λ = getduals(subproblem.feasibility_solver)
    cols = zeros(length(subproblem.masterterms))
    vals = zeros(length(subproblem.masterterms))
    for (s, (i,j,coeff)) in enumerate(subproblem.masterterms)
        cols[s] = j
        vals[s] = -λ[i]*coeff
    end
    G = sparsevec(cols,vals,length(subproblem.x))
    g = getobjval(subproblem.feasibility_solver)+G⋅subproblem.x

    return FeasibilityCut(G, g, subproblem.id)
end

Infeasible(subprob::SubProblem) = Infeasible(subprob.id)
Unbounded(subprob::SubProblem) = Unbounded(subprob.id)

function fill_submodel!(submodel::JuMP.Model, subproblem::SubProblem)
    fill_submodel!(submodel, get_solution(subproblem)...)
    return nothing
end

function fill_submodel!(submodel::JuMP.Model, x::AbstractVector, μ::AbstractVector, λ::AbstractVector, C::AbstractFloat)
    submodel.colVal = x
    submodel.redCosts = μ
    submodel.linconstrDuals = λ
    submodel.objVal = C
    submodel.objVal *= submodel.objSense == :Min ? 1 : -1
    return nothing
end


text:
function build()
    p = pwd()
    cd(dirname(@__FILE__))
    if Sys.iswindows()
        try
            run(`mingw32-make`)
        catch
            if Sys.WORD_SIZE == 32
                url = "https://github.com/sunoru/RandomNumbers.jl/releases/download/deplib-0.1/librandom123-32.dll"
            else
                url = "https://github.com/sunoru/RandomNumbers.jl/releases/download/deplib-0.1/librandom123.dll"
            end
            @info("You don't have MinGW32 installed, so it's now downloading the library binary from github.")
            download(url, "librandom123.dll")
        end
    elseif Sys.isbsd() && !Sys.isapple()  # e.g. FreeBSD
        run(`gmake`)
    else
        run(`make`)
    end
    cd(p)
end

function have_aesni()
    @static if VERSION < v"0.5-" || Sys.WORD_SIZE != 64
        return false
    else
        ecx = Base.llvmcall(
            """%1 = call { i32, i32, i32, i32 } asm "xchgq  %rbx,\${1:q}\\0A  cpuid\\0A  xchgq  %rbx,\${1:q}",
            "={ax},=r,={cx},={dx},0,~{dirflag},~{fpsr},~{flags}"(i32 1)
            %2 = extractvalue { i32, i32, i32, i32 } %1, 2
            ret i32 %2""", UInt32, Tuple{})
        return (ecx >> 25) & 1 == 1
    end
end

check_compiler() = Sys.iswindows() ? true : success(`gcc --version`)

if have_aesni() && check_compiler()
    build()
else
    @warn("AES-NI will not be compiled.")
end


text:
using SnpArrays, Distributions, Test, Statistics, VarianceComponentModels, Distributions, Random
using DataFrames, LinearAlgebra, StatsBase, TraitSimulation, Distributions
Random.seed!(1234)

# demo genetic data can be found in Snparrays package directory
filename = "EUR_subset"
EUR = SnpArray(SnpArrays.datadir(filename * ".bed"));
GRM = grm(EUR, minmaf = 0.05);
EUR_data = SnpData(SnpArrays.datadir(filename));
bimfile = EUR_data.snp_info # store the snp_info with the snp names
snpid  = bimfile[:, :snpid] # store the snp names in the snpid vector
causal_snp_index = findfirst(x -> x ==  "rs62057050" , snpid) # find the index of the snp of interest by snpid

locus = convert(Vector{Float64}, @view EUR[:, causal_snp_index])
famfile = EUR_data.person_info
sex = map(x -> strip(x) == "F" ? 0.0 : 1.0, famfile[!, :sex])
X = [sex locus]
β_full =
 [0.52
0.05]

I_n = Matrix(I, size(GRM))
vc = @vc [0.01][:,:] ⊗ (GRM + I_n) + [0.9][:, :] ⊗ I_n
trait = VCMTrait(X, β_full[:, :], vc)
X = trait.X
V1 = trait.vc[1].V
V2 = trait.vc[2].V

nsim = 5
γs = collect(0.0:0.05:0.4)
Y = zeros(size(trait.μ))
X = trait.X
V1 = trait.vc[1].V
V2 = trait.vc[2].V
@test eltype(Y) == Float64

# allocate variance component objects for null and alternative hypothesis
vcm_null, vcrot_null, vcm_alt, vcrot_alt = @time null_and_alternative_vcm_and_rotate(Y, X, V1, V2)
pvalue = power_simulation(trait, γs, nsim, vcm_null, vcrot_null, vcm_alt, vcrot_alt)

@test size(pvalue) == (nsim, length(γs))

alpha = 0.0000005
powers = power(pvalue, alpha)
@test powers[1] < powers[end]


text:
# ScalarField related functions loaded when Interpolations.jl is available.

"""
    $SIGNATURES

Load gridded elevation data `dem` into a factor graph `fg` as a collection
of Point3 variables. Each variable is connected to its 4-neighborhood by
relative Point3Point3 MvNormal constraints with mean defined by their
relative position on the grid (`x`, `y`) and covariance `meshEdgeSigma`.
"""
function _buildGraphScalarField!( fg::AbstractDFG,
                                  dem::AbstractMatrix, # assume grayscale image for now
                                  x::AbstractVector,
                                  y::AbstractVector;
                                  solvable::Int=0,
                                  marginalized::Bool=true,
                                  meshEdgeSigma=diagm([1;1;1]),
                                  refKey::Symbol = :simulated )
  #
  # assume regular grid
  dx, dy = x[2]-x[1], y[2]-y[1]
  for i in 1:length(x)
    for j in 1:length(y)
      s = Symbol("pt$(i)_$(j)") # unique identifier
      pt = addVariable!(fg, s, Point3, solvable=solvable)
      setMarginalized!(pt, marginalized) # assume solveKey=:default

      # ...
      refVal = [x[i];y[j];dem[i,j]]
      simPPE = DFG.MeanMaxPPE(refKey, refVal, refVal, refVal)
      setPPE!(pt, refKey, typeof(simPPE), simPPE)

      # Regular grid triangulation:
      #  add factor to (i-1,j)   |
      #  add factor to (i, j-1)  -
      #  add factor to (i-1, j-1)  \

      # no edges to prev row on first row
      if i>1
        dVal1 = dem[i,j]-dem[i-1,j]
        f = Point3Point3(MvNormal([dx, 0, dVal1], meshEdgeSigma))
        addFactor!(fg, [Symbol("pt$(i-1)_$(j)"), s], f, solvable=solvable, graphinit=false)
      end

      # no edges to prev column on first column
      if j>1
        dVal2 = dem[i,j]-dem[i,j-1]
        f = Point3Point3(MvNormal([0, dy, dVal2], meshEdgeSigma))
        addFactor!(fg, [Symbol("pt$(i)_$(j-1)"),s], f, solvable=solvable, graphinit=false)
      end

      # no edges to add on first element
      if i>1 && j>1
        dVal3 = dem[i,j]-dem[i-1,j-1]
        f = Point3Point3(MvNormal([dx, dy, dVal3], meshEdgeSigma))
        addFactor!(fg,[Symbol("pt$(i-1)_$(j-1)"),s], f, solvable=solvable, graphinit=false)
      end

    end
  end

  nothing
end

#


text:
module JuliaGridSynth

using LinearAlgebra
using SparseArrays
# Write your package code here.
include("Gates.jl")
export
    AbstractGate,
    X,
    Y,
    Z,
    H,
    S,
    Sdag,
    T,
    Tdag,
    XGate,
    YGate,
    ZGate,
    HGate,
    SGate,
    TGate,
    SdagGate,
    TdagGate,
    RZGate,
    matrix

end

max_stars_repo_name: isolatedinformation/JuliaGridSynth

max_issues_repo_name: isolatedinformation/JuliaGridSynth

max_forks_repo_name: isolatedinformation/JuliaGridSynth


text:
const ENVIRONMENTS_DIR = DFC.config_path("environments")

"""
    Environment(MPI_command::String, scheduler_flags::Vector{String}, exports::Vector{String})

Environment to run a [`Job`](@ref) in. When running on a server with a scheduler `scheduler_flags` holds what these should be. e.g. `#SBATCH -N 2`. `MPI_command` and `MPI_processes` will be used to prefix executables that should be ran in parallel, e.g. if `MPI_command = "mpirun -np 4 --bind-to core"` a parallel executable will be translated into a script line as `mpirun -np 4 --bind-to core exec`.
"""
@with_kw mutable struct Environment
    MPI_command::String
    scheduler_flags::Vector{String} = String[]
    exports::Vector{String} = String[]
end

StructTypes.StructType(::Type{Environment}) = StructTypes.Struct()
function Base.:(==)(e1::Environment, e2::Environment)
    if e1.MPI_command != e2.MPI_command
        return false
    else
        for fn in (:scheduler_flags, :exports)
            for f1 in getfield(e1, fn)
                if !(f1 ∈ getfield(e2, fn))
                    return false
                end
            end
        end
    end
    return true
end

"Reads the script file and tries to understand what the runtime environment is."
function environment_from_jobscript(scriptpath::String)
    lines = filter(!isempty, strip.(readlines(scriptpath)))
    id = findfirst(l -> any(x -> occursin(x, l), Calculations.RUN_EXECS), lines)
    if id !== nothing
        sline = split(replace(lines[id], "#" => ""))

        MPI_command = join(sline[1:findnext(x -> splitdir(x)[end] ∈ Calculations.allexecs() , sline, 2)-1], " ")
    else
        MPI_command = ""
    end
    #TODO only works with slurm!!!
    scheduler_flags = map(y -> strip(split(y, "SBATCH")[end]), filter(x -> x[1:7] == "#SBATCH", lines))
    deleteat!(scheduler_flags, findall(x -> occursin("-J", x) || occursin("--job-name", x), scheduler_flags))
    exports = map(y -> strip(split(y, "export")[end]), filter(x -> x[1:6] == "export", lines))
    return Environment(MPI_command, scheduler_flags, exports)
end

function load_environment(name::String)
    fullpath = joinpath(ENVIRONMENTS_DIR, name * ".json")
    return ispath(fullpath) ? JSON3.read(read(fullpath, String), Environment) : nothing
end

function save(env::Environment, name::String)
    if !ispath(ENVIRONMENTS_DIR)
        mkpath(ENVIRONMENTS_DIR)
    end
    return JSON3.write(joinpath(ENVIRONMENTS_DIR, name * ".json"), env)
end

function environment_name(env::Environment)
    if ispath(ENVIRONMENTS_DIR)
        for f in readdir(ENVIRONMENTS_DIR)
            name = splitext(f)[1]
            env2 = load_environment(name)
            if env2 == env
                return name
            end
        end
    end
end

function Base.write(f::IO, env::Environment)
    for flag in env.scheduler_flags
        write(f, "#SBATCH $flag\n")
    end
    for flag in env.exports
        write(f, "export $flag\n")
    end
end

function rm_environment!(name::AbstractString)
    fullpath = joinpath(ENVIRONMENTS_DIR, name * ".json")
    if ispath(fullpath)
        rm(fullpath)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

import .CuArrays

# Aovid `mat + a * I` with CuArrays which involes scalar operations and is slow
safe_diagm(mat::CuArrays.CuArray, a::T) where {T} = CuArrays.CuArray{T}(diagm(0 => fill(a, size(mat, 1))))


text:
# Helper data for color parsing
include("names_data.jl")

# Color Parsing
# -------------

const col_pat_hex1 = r"(#|0x)([[:xdigit:]])([[:xdigit:]])([[:xdigit:]])"
const col_pat_hex2 = r"(#|0x)([[:xdigit:]]{2})([[:xdigit:]]{2})([[:xdigit:]]{2})"
const col_pat_rgb  = r"rgb\((\d+%?),(\d+%?),(\d+%?)\)"
const col_pat_hsl  = r"hsl\((\d+%?),(\d+%?),(\d+%?)\)"
const col_pat_rgba = r"rgba\((\d+%?),(\d+%?),(\d+%?),(\d+(?:\.\d*)?%?)\)"
const col_pat_hsla = r"hsla\((\d+%?),(\d+%?),(\d+%?),(\d+(?:\.\d*)?%?)\)"

# Parse a number used in the "rgb()" or "hsl()" color.
function parse_rgb_hsl_num(num::AbstractString)
    if num[end] == '%'
        return parse(Int, num[1:end-1], 10) / 100
    else
        return parse(Int, num, 10) / 255
    end
end

# Parse a number used in the alpha field of "rgba()" and "hsla()".
function parse_alpha_num(num::AbstractString)
    if num[end] == '%'
        return parsefloat(num[1:end-1]) / 100
    else
        return parsefloat(num)
    end
end

# Parse a color description.
#
# This parses subset of HTML/CSS color specifications. In particular, everything
# is supported but: "currentColor".
#
# It does support named colors (though it uses X11 named colors, which are
# slightly different than W3C named colors in some cases), "rgb()", "hsl()",
# "#RGB", and "#RRGGBB' syntax.
#
# Args:
#   desc: A color name or description.
#
# Returns:
#   An RGB color, unless "hsl()" was used, in which case an HSL color.
#
function color(desc::AbstractString)
    desc_ = replace(desc, " ", "")
    mat = match(col_pat_hex2, desc_)
    if mat != nothing
        return RGB(parse(Int, mat.captures[2], 16) / 255,
                   parse(Int, mat.captures[3], 16) / 255,
                   parse(Int, mat.captures[4], 16) / 255)
    end

    mat = match(col_pat_hex1, desc_)
    if mat != nothing
        return RGB((16 * parse(Int, mat.captures[2], 16)) / 255,
                   (16 * parse(Int, mat.captures[3], 16)) / 255,
                   (16 * parse(Int, mat.captures[4], 16)) / 255)
    end

    mat = match(col_pat_rgb, desc_)
    if mat != nothing
        return RGB(parse_rgb_hsl_num(mat.captures[1]),
                   parse_rgb_hsl_num(mat.captures[2]),
                   parse_rgb_hsl_num(mat.captures[3]))
    end

    mat = match(col_pat_hsl, desc_)
    if mat != nothing
        return HSL(parse_rgb_hsl_num(mat.captures[1]),
                   parse_rgb_hsl_num(mat.captures[2]),
                   parse_rgb_hsl_num(mat.captures[3]))
    end

    mat = match(col_pat_rgba, desc_)
    if mat != nothing
        return RGBA(parse_rgb_hsl_num(mat.captures[1]),
                    parse_rgb_hsl_num(mat.captures[2]),
                    parse_rgb_hsl_num(mat.captures[3]),
                    parse_alpha_num(mat.captures[4]))
    end

    mat = match(col_pat_hsla, desc_)
    if mat != nothing
        return HSLA(parse_rgb_hsl_num(mat.captures[1]),
                    parse_rgb_hsl_num(mat.captures[2]),
                    parse_rgb_hsl_num(mat.captures[3]),
                    parse_alpha_num(mat.captures[4]))
    end

    desc_ = lowercase(desc_)

    if desc_ == "transparent"
        return RGBA(0.0, 0.0, 0.0, 0.0)
    end

    if !haskey(color_names, desc_)
        error("Unknown color: ", desc)
    end

    c = color_names[desc_]
    return RGB(c[1] / 255, c[2] / 255, c[3] / 255)
end

color(c::ColorValue) = c
color(c::AlphaColorValue) = c


text:
@rule Transition(:a, Marginalisation) (q_out::Any, q_in::Categorical) = begin
    return MatrixDirichlet(probvec(q_out)*probvec(q_in)' .+ 1.0)
end

@rule Transition(:a, Marginalisation) (q_out_in::Contingency, ) = begin
    return MatrixDirichlet(contingency_matrix(q_out_in) .+ 1.0)
end


text:
"""
writepheno
writepheno will simulate the traits under the given model for GLM and OrderedMultinomial traits (univariate) and write the simulation results to tab separated PLINK formatted phenotype file titled "pheno.txt".
The first two columns of the phenotype file will be the Family ID and the Individual ID the working directory `file_path`, and will find the family information stored in the ".fam" file in this data directory.
"""
function writepheno(plink_file::AbstractString, file_path::AbstractString, traitobject::Union{GLMTrait, OrderedMultinomialTrait}, nsim::Int; famname::AbstractString=plink_file*".fam")
    # simulate phenotypes
    pheno = DataFrame(simulate(traitobject, nsim), :auto)
    rename!(pheno, [Symbol("Trait$i") for i in 1:nsim])

    # load person info
    dlmdata = readdlm(file_path*"/"*famname, AbstractString)
    famfile = DataFrame(Tables.table(dlmdata, header = Symbol.(:x, axes(dlmdata,2))))
    personinfo = famfile[:, 1:2]
    outfile = hcat(personinfo, pheno)
    open("pheno.txt"; write=true) do f
        write(f, "#FID, IID, $nsim simulations replications\n")
        writedlm(f, Matrix(outfile))
    end
end

"""
writepheno
writepheno will simulate the traits under the given model for Variance Component Model traits and GLMM traits (multivariate) and write the simulation results to tab separated PLINK formatted phenotype file titled "pheno.txt".
The first two columns of the phenotype file will be the Family ID and the Individual ID the working directory `file_path`, and will find the family information stored in the ".fam" file in this data directory.
"""
function writepheno(plink_file::AbstractString, file_path::AbstractString, traitobject::Union{VCMTrait, GLMMTrait}, nsim::Int; famname::AbstractString=plink_file*".fam")
    # simulate phenotypes
    n_traits = ntraits(traitobject)
    phenotypes = simulate(traitobject, nsim)
    pheno = [phenotypes[i][:, j] for j in 1:n_traits, i in 1:nsim]
    pheno = DataFrame(hcat(pheno...), :auto)
    rename!(pheno[:,1:n_traits:n_traits*nsim], [Symbol("Simulation$i") for i in 1:nsim])

    # load person info
    dlmdata = readdlm(file_path*"/"*famname, AbstractString)
    famfile = DataFrame(Tables.table(dlmdata, header = Symbol.(:x, axes(dlmdata,2))))
    personinfo = famfile[:, 1:2]
    rename!(personinfo, [Symbol("FID"), Symbol("IID")])

    outfile = hcat(personinfo, pheno)

    open("pheno.txt"; write=true) do f
        write(f, "#FID, IID, $nsim simulations replications for $n_traits traits\n")
        writedlm(f, Matrix(outfile))
    end
end


text:
using Documenter, SignalBase

DocMeta.setdocmeta!(SignalBase, :DocTestSetup, :(using SignalBase; using SignalBase.Units); recursive=true)

makedocs(;
    modules=[SignalBase],
    format=Documenter.HTML(),
    pages=[
        "Reference" => "index.md",
    ],
    repo="https://github.com/haberdashpi/SignalBase.jl/blob/{commit}{path}#L{line}",
    sitename="SignalBase.jl",
    authors="David F. Little",
)

deploydocs(;
    repo="github.com/haberdashPI/SignalBase.jl",
)


text:
# %% markdown
# # selfplay.jl
# %%
include("selfPlay.jl")
# %%
selfPlay()
# %%
using ProgressMeter
# %%
@showprogress
# %%
using Distributed

using ProgressMeter
using Printf: @sprintf
using Random

addprocs(8)

@everywhere begin
    include("./mcts.jl")
    include("./utils.jl")
    using SharedArrays
    using ParallelDataTransfer
end
# %%
function catAll(x, y)
    @sync begin
        @async x[1] = cat(dims=4, x[1], y[1])
        @async x[2] = cat(dims=2, x[2], y[2])
        @async x[3] = cat(dims=1, x[3], y[3])
    end
    return x
end

function progress_meter_string(p::Progress)
    if p.counter < p.n
        percentage_complete = round(Int, p.counter / p.n *100)
    elseif p.counter > p.n
        percentage_complete = 99
    else
        percentage_complete = 100
    end
    bar = ProgressMeter.barstring(60, percentage_complete, barglyphs=p.barglyphs)
    @sprintf "%s%3u%%  %s " p.desc percentage_complete bar
end

function pregress_to_string(id::Int, steps::Int, statu::Bool)
    pre = @sprintf "|    %d   |  %2d   " id steps
    pre *= '|' * "❌⚪" ^ (steps÷2) * ('❌'^(steps % 2)) * ('⭐'^statu)
    @sprintf "%-80s" pre
end

function showProgressBar(sha::SharedArray, status::SharedArray)
    mainPro = Progress(BATCH_SIZE, dt=0.1, color=:black, output=stdout)
    mainPro.desc = "|  main  | "
    total = 0
    desc = "+-------------------------------- Self Playing --------------------------------+\n\n| Gameid | Steps | Time: "
    open("./logs/selfPlay.log", "a+") do io
        println(io, desc)
        print(io, Dates.format(now(), "mm-dd-HH:MM"))
    end
    finished_game = 0
    while finished_game < 8
        finished_game, total = sum(status), sum(sha)
        print_string = desc
        IJulia.clear_output(true)

        if sum(status) == 7
            gap = " " ^ ((total - BATCH_SIZE) % 7)
            print_string *= gap*"saving data..."
        end

        print_string *= ProgressMeter.durationstring(time() - mainPro.tfirst) * "\n"

        if finished_game != 8
            mainPro.counter = total
            print_string *= progress_meter_string(mainPro) * "\n"
        else
            mainPro.counter = BATCH_SIZE
            print_string *= progress_meter_string(mainPro) * " Done!!\n"
        end

        for i=1:7
            print_string *= pregress_to_string(i, sha[i], status[i]) *"\n"
        end
        print_string *= "\n+-------------------------------- Self Playing --------------------------------+\n"
        println(print_string)
        sleep(0.2)
    end
    open("./logs/selfPlay.log", "a+") do io
        println(io, "Time: "*Dates.format(now(), "mm-dd-HH:MM"))
    end
end
# %%
@everywhere function playonce(nn::NeuralNet, sha::SharedArray{Int}, status::SharedArray{Bool}, gameid::Int)
    "play game once with neural net for generating training data"
    sum(sha) > BATCH_SIZE && return []
    chunk = []
    player = MCTS()
    game = Game()
    randomPlace(game, 2)
    while isnothing(game.winner)
        move, p = predict(game, player, nn, test=false)
        push!(chunk, [getGameState(game), p, game.nextPlayer])
        place(game, move)
        sha[gameid] +=1
    end
    @async save_log("selfPlay", game)
    for i in eachindex(chunk)
        chunk[i][3] *= game.winner
    end
    status[gameid] = true
    return chunk
end
# %%

function selfPlay()
    "self play n batches and save the data"
    @sync begin
        if length(readdir(MODEL_PATH))<2
            @everywhere nn = NeuralNet(9, 9, 4)
            save_model(nn)
        else
            @everywhere nn = load_model()
        end
        sha = SharedArray{Int}(7); status = SharedArray{Bool}(8)
        @async showProgressBar(sha, status);

        examples = @distributed append! for gameid=1:8
            playonce(nn, sha, status, gameid)
        end
        Random.shuffle!(examples)
        onebatch = reduce(catAll, examples[1:BATCH_SIZE])
        save_data(onebatch)
        status .= true
        nothing
    end
end
# %%
Main.selfPlay()
# %%
@everywhere nn = NeuralNet(9, 9, 4)
# %%
sha = SharedArray{Int}(7); status = SharedArray{Bool}(8)
@async showProgressBar(sha, status);
examples = @distributed append! for gameid=1:7
    playonce(nn, sha, status, gameid)
end
# %%
status .= true
# %%
status .= true
# %%
nn= NeuralNet(9, 9, 4)
# %%
sha = zeros(Int8, 7)
status = zeros(Bool, 8)
# %%
examples = @distributed append! for gameid=1:7
    playonce(nn, sha, status, gameid)
end
# %%
status .= true
# %%
status .= true
# %%
examples = pmap(1:7) do gameid
    playonce(nn, sha, status, gameid)
end
# %%
examples = []
for work in workers()
    r = @spawn work playonce(nn, sha, status)
    append!(examples, fetch(r))
end
# %%
nn = NeuralNet(9, 9, 4)
# %%
examples = []
for work in workers()
    append!(examples, fetch(@spawnat worker playonce(nn, sha, status)))
end
# %%
status .= true
# %%
@everywhere using InteractiveUtils
# %%
varinfo()
# %%


text:
module Mueller

using StaticArrays

export rotation, rotate, linear_polarizer, hwp, qwp, mirror, waveplate

include("plots.jl")

# Various formulae for polarization components
"""
    rotation([T=Float64], θ)

Generate a rotation matrix with the given angle, in radians, `θ`. This can be used to rotate the axes of polarization components arbitrarily. For convenience, the [`rotate`](@ref) method will rotate a component, without having to generate this matrix, itself.

# Examples

Rotate a linear polarizer by 90 degrees counter-clockwise
```jldoctest
julia> M = linear_polarizer();

julia> r = rotation(π/4);

julia> Mr = r' * M * r;

julia> Mr ≈ linear_polarizer(π/4)
true
```

# See also
[`rotate`](@ref)
"""
function rotation(T, θ)
    sin2t, cos2t = sincos(2 * θ)
    return SA{T}[1 0 0 0
              0 cos2t sin2t 0
              0 -sin2t cos2t 0
              0 0 0 1]
end
rotation(θ) = rotation(Float64, θ)

"""
    rotate(M, θ)

Rotates the component represented by the Mueller matrix `M` counter-clockwise by angle `θ` (in radians).

# Examples
```jldoctest
julia> M = linear_polarizer();

julia> Mr = rotate(M, π/2);

julia> Mr ≈ linear_polarizer(π/2)
true
```

# See also
[`Mueller.rotation`](@ref)
"""
function rotate(mat::AbstractMatrix{T}, θ) where T
    r = rotation(T, θ)
    return r' * mat * r
end

"""
    linear_polarizer([T=Float64], θ=0; p=1)

A linear polarizer with the throughput axis given by `θ`, in radians, by default horizontal. The partial polarization can be given with the `p` keyword argument, which changes the intensity by a factor of `p^2/2`.

# Examples

```jldoctest
julia> M = linear_polarizer()
4×4 StaticArrays.SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 0.5  0.5  0.0  0.0
 0.5  0.5  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> S = [1, 0, 0, 0]; # I, Q, U, V

julia> M * S # only horizontal component (+Q) remains
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
 0.5
 0.5
 0.0
 0.0
```
"""
function linear_polarizer(T::Type, θ=0; p=1)
    I = T(p^2 / 2)
    sin2t, cos2t = sincos(2 * θ)
    M = I * SA{T}[1 cos2t sin2t 0
                 cos2t cos2t^2 cos2t * sin2t 0
                 sin2t cos2t * sin2t sin2t^2 0
                 0 0 0 0]
    return M
end
linear_polarizer(θ=0; kwargs...) = linear_polarizer(Float64, θ; kwargs...)

# wave plates

"""
    waveplate([T=Float64], δ=0, θ=0)

A generic phase retarder (waveplate) with fast axis aligned with angle `θ`, in radians, and phase delay of `δ`, in radians, along the slow axis.

# Examples
```jldoctest
julia> M = waveplate(π);

julia> M ≈ hwp()
true

julia> M = waveplate(π/2, π/4);

julia> M ≈ qwp(π/4)
true
```

# See also
[`hwp`](@ref), [`qwp`](@ref), [`mirror`](@ref)
"""
function waveplate(T::Type, δ=0, θ=0)
    sin2t, cos2t = sincos(2 * θ)
    sind, cosd = sincos(δ)
    return SA{T}[1 0 0 0
                 0 cos2t^2 + sin2t^2*cosd cos2t * sin2t * (1 - cosd) sin2t * sind
                 0 cos2t*sin2t*(1-cosd) cos2t^2*cosd + sin2t^2 -cos2t * sind
                 0 -sin2t * sind cos2t * sind cosd]
end
waveplate(δ=0, θ=0) = waveplate(Float64, δ, θ)

"""
    hwp([T=Float64], θ=0)

A half-wave plate (HWP) with fast axis oriented at angle `θ`, in radians.

# Examples
```jldoctest
julia> M = hwp()
4×4 StaticArrays.SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 1.0   0.0   0.0           0.0
 0.0   1.0   0.0           0.0
 0.0   0.0  -1.0          -1.22465e-16
 0.0  -0.0   1.22465e-16  -1.0

julia> S = [1, 1, 0, 0]; # I, Q, U, V

julia> M * S # allow +Q through unchanged
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
 1.0
 1.0
 0.0
 0.0

julia> rotate(M, π/8) * S # switch +Q to +U
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
  1.0
  1.9967346175427393e-16
  1.0
 -8.659560562354932e-17

```

# See also
[`waveplate`](@ref), [`qwp`](@ref)
"""
hwp(T::Type, θ=0) = waveplate(T, π, θ)
hwp(θ=0) = hwp(Float64, θ)

"""
    qwp([T=Float64], θ=0)

A quarter-wave plate (QWP) with fast axis oriented at angle `θ`, in radians.

# Examples
```jldoctest
julia> M = qwp()
4×4 StaticArrays.SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 1.0   0.0  0.0           0.0
 0.0   1.0  0.0           0.0
 0.0   0.0  6.12323e-17  -1.0
 0.0  -0.0  1.0           6.12323e-17

julia> S = [1, 1, 0, 0]; # I, Q, U, V

julia> M * S # allow +Q through unchanged
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
 1.0
 1.0
 0.0
 0.0

julia> qwp(-π/4) * S # switch +Q to +V
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
  1.0
  6.123233995736766e-17
 -6.123233995736765e-17
  1.0
```

# See also
[`waveplate`](@ref), [`hwp`](@ref)
"""
qwp(T::Type, θ=0) = waveplate(T, π/2, θ)
qwp(θ=0) = qwp(Float64, θ)

"""
    mirror([T=Float64], r=1, δ=π, θ=0)

A reflective mirror with reflectance `r`, oriented at angle `θ`, in radians, compared to the reference frame of the light, and with phase shift `δ`, in radians. An ideal mirror will have perfect reflectance (1) and a half-wave phase shift (π).

# Examples
```jldoctest
julia> M = mirror()
4×4 StaticArrays.SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 1.0  0.0   0.0           0.0
 0.0  1.0   0.0          -0.0
 0.0  0.0  -1.0           1.22465e-16
 0.0  0.0  -1.22465e-16  -1.0

julia> S = [1, 1, 0, 0]; # I, Q, U, V

julia> M * S # no change
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
 1.0
 1.0
 0.0
 0.0

julia> mirror(1, π, π/4) * S # rotates polarized light
4-element StaticArrays.SVector{4, Float64} with indices SOneTo(4):
  1.0
 -1.0
  1.2246467991473532e-16
  1.2246467991473532e-16
```
"""
function mirror(T::Type, r=1, δ=π, θ=0)
    a = 0.5 * (r + 1)
    b = 0.5 * (r - 1)
    sin2t, cos2t = sincos(2 * θ)
    sin4t = sin(4 * θ)
    sind, cosd = sincos(δ)
    sqrm = sqrt(r)
    M = SA{T}[a          b * cos2t                            b * sin2t                           0
              b * cos2t  a * cos2t^2 + sqrm * cosd * sin2t^2 (a - sqrm * cosd) * sin4t * 0.5     -sqrm * sind * sin2t
              b * sin2t (a - sqrm * cosd) * sin4t * 0.5       a * sin2t^2 + sqrm * cosd * cos2t^2 sqrm * sind * cos2t
              0          sqrm * sind * sin2t                 -sqrm * sind * cos2t                 sqrm * cosd]
    return M
end
mirror(r=1, δ=π, θ=0) = mirror(Float64, r, δ, θ)

end # module


text:
# Profiler interface

export
    @cuprofile

initialized = false

start_profiler() = @cucall(:cuProfilerStart, (Ptr{Void},), 0)
stop_profiler() = @cucall(:cuProfilerStop, (Ptr{Void},), 0)

macro cuprofile(ex)
    quote
        start_profiler()
        try
            $(esc(ex))
        finally
            stop_profiler()
        end
    end
end