text:
function test_annealing(target, init, cost, modifier, cooldown)
    system = CorrelationTracker(target |> init, target)
    furnace = Furnace(system, target; T0 = 0.0)
    c0 = euclid_directional(system, target)

    for i in 1:10000
        furnace = annealing_step(furnace;
                                 cost     = cost,
                                 modifier = modifier,
                                 cooldown = cooldown)
    end

    c1 = euclid_directional(system, target)

    @test c1 < c0
end

@testset "300x300 image" begin
    array = read_cuboid("image3d.json")[:,:,1]
    target = CorrelationTracker(array;
                                directions = [:x, :y, :xy, :yx],
                                periodic   = true)

    # Need to find starting costs for this
    #modifiers = (InterfaceFlipper(), InterfaceSwapper(),
    #             RandomSwapper(),    RandomFlipper())
    #costs = (euclid_mean, euclid_directional)
    #initializers = (initialize_random, initialize_spheres)
    #cooldowns = (exponential_cooldown(), aarts_korst_cooldown(), frost_heineman_cooldown())

    #alltogether = product(modifiers, costs, initializers, cooldowns)

    alltogether = product((Flipper(InterfaceSampler()), Swapper(InterfaceSampler())),
                          (euclid_directional, euclid_mean),
                          (initialize_spheres,),
                          (aarts_korst_cooldown(), exponential_cooldown()))

    foreach(alltogether) do prod
        modifier, cost, initializer, cooldown = prod
        test_annealing(target, initializer, cost, modifier, cooldown)
    end
end


text:
using Pidoh
using Documenter

makedocs(;
    modules = [Pidoh],
    authors = "Amirhossein Rajabi",
    repo = "https://github.com/Pidoh/Pidoh.jl/blob/{commit}{path}#L{line}",
    sitename = "Pidoh.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://Pidoh.github.io/Pidoh.jl",
        assets = String[],
        collapselevel = 1,
    ),
    pages = [
        "Pidoh" => "index.md",
        "Algorithms" => Any["simulatedannealing.md"],
        "Problems" => Any["bitstringproblems.md", "problems.md"],
        "Instances" => "instances.md",
    ],
)

deploydocs(; repo = "github.com/Pidoh/Pidoh.jl")

ext: jl


text:
### MHState

# MHState holds the internal state ("local variables") of the Metropolis-Hastings sampler

mutable struct MHState{S<:ValueSupport, F<:VariateForm} <: MHSamplerState{F}
  proposal::Distribution{F, S} # Proposal distribution
  pstate::ParameterState{S, F} # Parameter state used internally by MH
  tune::MCTunerState
  ratio::Real # Acceptance ratio
  diagnosticindices::Dict{Symbol, Integer}

  function MHState{S, F}(
    proposal::Distribution{F, S},
    pstate::ParameterState{S, F},
    tune::MCTunerState,
    ratio::Real,
    diagnosticindices::Dict{Symbol, Integer}
  ) where {S<:ValueSupport, F<:VariateForm}
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    new(proposal, pstate, tune, ratio, diagnosticindices)
  end
end

MHState(
  proposal::Distribution{F, S},
  pstate::ParameterState{S, F},
  tune::MCTunerState,
  ratio::Real
) where {S<:ValueSupport, F<:VariateForm} =
  MHState{S, F}(proposal, pstate, tune, ratio, Dict{Symbol, Integer}())

MHState(
  proposal::Distribution{F, S},
  pstate::ParameterState{S, F},
  tune::MCTunerState=BasicMCTune()
) where {S<:ValueSupport, F<:VariateForm} =
  MHState(proposal, pstate, tune, NaN)

### Metropolis-Hastings (MH) sampler

# In its most general case it accommodates an asymmetric proposal distribution
# For symmetric proposals, the proposal correction factor equals 1, so the logproposal field is set to nothing
# For non-normalised proposals, the lognormalise() method is used for calculating the proposal correction factor

struct MH <: MHSampler
  symmetric::Bool # If symmetric=true then the proposal distribution is symmetric, else it is asymmetric
  normalised::Bool # If normalised=true then the proposal distribution is normalised, else it is non-normalised
  setproposal::Function # Function for setting the proposal distribution
end

MH(setproposal::Function; signature::Symbol=:high, args...) = MH(setproposal, Val{signature}; args...)

MH(setproposal::Function, ::Type{Val{:low}}; symmetric::Bool=true, normalised::Bool=true) =
  MH(symmetric, normalised, setproposal)

MH(setproposal::Function, ::Type{Val{:high}}; symmetric::Bool=true, normalised::Bool=true) =
  MH(symmetric, normalised, _state -> setproposal(_state.value))

# Random-walk Metropolis, i.e. Metropolis with a normal proposal distribution

MH(σ::Matrix{N}) where {N<:Real} = MH(x::Vector{N} -> MvNormal(x, σ), signature=:high)
MH(σ::Vector{N}) where {N<:Real} = MH(x::Vector{N} -> MvNormal(x, σ), signature=:high)
MH(σ::N) where {N<:Real} = MH(x::N -> Normal(x, σ), signature=:high)
MH(::Type{N}=Float64) where {N<:Real} = MH(x::N -> Normal(x, 1.0), signature=:high)

### Initialize Metropolis-Hastings sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{S, F},
  parameter::Parameter{S, F},
  sampler::MH,
  outopts::Dict
) where {S<:ValueSupport, F<:VariateForm}
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array{Any}(length(pstate.diagnostickeys))
  end
end

## Initialize MHState

function sampler_state(
  parameter::Parameter{S, F},
  sampler::MH,
  tuner::MCTuner,
  pstate::ParameterState{S, F},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
) where {S<:ValueSupport, F<:VariateForm}
  sstate = MHState(sampler.setproposal(pstate), generate_empty(pstate), tuner_state(parameter, sampler, tuner))
  set_diagnosticindices!(sstate, [:accept], diagnostickeys)
  sstate
end

## Reset parameter state

function reset!(pstate::UnivariateParameterState, x::Real, parameter::UnivariateParameter, sampler::MH)
  pstate.value = x
  parameter.logtarget!(pstate)
end

function reset!(pstate::MultivariateParameterState, x::RealVector, parameter::MultivariateParameter, sampler::MH)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

reset!(
  sstate::MHState{S, F},
  pstate::ParameterState{S, F},
  parameter::Parameter{S, F},
  sampler::MCSampler,
  tuner::MCTuner
) where {S<:ValueSupport, F<:VariateForm} =
  reset!(sstate.tune, sampler, tuner)

function show(io::IO, sampler::MH)
  issymmetric = sampler.symmetric ? "symmetric" : "non-symmmetric"
  isnormalised = sampler.normalised ? "normalised" : "non-normalised"
  print(io, "MH sampler: $issymmetric $isnormalised proposal")
end


text:
import CommandLine
import ManagerTypes
#
# This manager is also a server for all other remote managers. These
# remote managers are clients for current and may use "Client" module
# to send commands and obtain results. So, this command line property
# will be used for current manager(server) for listening clients. For
# this you have to run Manager like this:
# >julia AppManager.jl port=2000
# Fast port is used for "fast" pooling mode.
# Next const is an IP of server. Is used in pair with ARG_SERVER_PORT
#
const ARG_SERVER_PORT      = "port"
const ARG_FAST_SERVER_PORT = "fastPort"
const ARG_SERVER_IP        = "ip"
#
# Current manager may be connected to the left-sided other server or
# instance. It means, that current world will be joined with that
# left-sided world and all organisms will be shared between both of
# them.
#
const ARG_LEFT_SERVER_PORT = "leftPort"
const ARG_LEFT_SERVER_IP   = "leftIp"
#
# Current manager may be connected to the right-sided other server or
# instance. It means, that current world will be joined with that
# right-sided world and all organisms will be shared between both of
# them.
#
const ARG_RIGHT_SERVER_PORT = "rightPort"
const ARG_RIGHT_SERVER_IP   = "rightIp"
#
# Current manager may be connected to the up-sided other server or
# instance. It means, that current world will be joined with that
# up-sided world and all organisms will be shared between both of
# them.
#
const ARG_UP_SERVER_PORT = "upPort"
const ARG_UP_SERVER_IP   = "upIp"
#
# Current manager may be connected to the down-sided other server or
# instance. It means, that current world will be joined with that
# down-sided world and all organisms will be shared between both of
# them.
#
const ARG_DOWN_SERVER_PORT = "downPort"
const ARG_DOWN_SERVER_IP   = "downIp"
#
# In this mode no terminal output will be produced
#
const ARG_QUIET            = "quiet"
#
# Returns port number set in command line or in Config object (default)
# @param man Manager data type
# @param cmdParam Name of command line parameter
# @param cfgParam Name of configuration parameter if cmdParam wasn't set
# @return {Int} Port number or 0 if no port available
#
function _getPort(man::ManagerTypes.ManagerData, cmdParam::String, cfgParam::Symbol)
  local port::String = CommandLine.val(man.params, cmdParam)
  port == "" ? getfield(man.cfg, cfgParam) : parse(Int, port)
end
#
# Returns IP address of current server. See ARG_SERVER_IP for details
# @param man Managerdata type
# @param cmdParam Name of command line parameter
# @param cfgParam Name of configuration parameter if cmdParam wasn't set
# @return {IPv4} IP address
#
function _getIp(man::ManagerTypes.ManagerData, cmdParam::String, cfgParam::Symbol)
  local ip::String = CommandLine.val(man.params, cmdParam)
  ip == "" ? getfield(man.cfg, cfgParam) : IPv4(ip)
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
Base.@kwdef mutable struct Parameters
    maxiter::Int = 3000     # maximum number of iterations
    ncuts_per_iter::Int = 1 # number of cuts added per iteration
    obj_limit::Float64 = -Inf # termination condition

    ϵ_s::Float64 = 1.e-5 # termination tolerance
    m_L::Float64 = 1.e-4 # serious step condition parameter (0, 0.5)

    # Proximal method
    u::Float64 = 1.e-2       # initial proximal penalty value
    u_min::Float64 = 1.e-6   # minimum proximal penalty value
    ϵ_float::Float64 = 1.e-8 # tolerance for floating point comparison
    m_R::Float64 = 0.5       # proximal term update parameter (0.5, 1)
    max_age::Int = 10        # maximum number of iterations before removing inactive cuts

    # Trust region method
    Δ_ub::Float64 = 1.e+3 # trust region bound upper limit
    Δ_lb::Float64 = 1.e-4 # trust region bound lower limit

end

function get_parameter(params::Parameters, pname::String)
    return getfield(params, Symbol(pname))
end

function set_parameter(params::Parameters, pname::String, val)
    setfield!(params, Symbol(pname), val)
    return nothing
end


text:
using ASP
using FactCheck

tests = ["as_nnls"]

for t in tests
    include("$(t).jl")
end


text:
# - - parameters - - - - - - - - - - -
system["saveU_dir"] = ""
system["verboselevel"] = 2
system["L"] = (4, 4, 4, 4)
system["Nwing"] = 1
system["Nsteps"] = 15
system["firstlearn"] = 10
system["quench"] = true
system["logfile"] = "SLMC_L04040404_beta5.7_Staggered_mass0.5.txt"
system["βeff"] = 5.7
system["initial"] = "hot"
system["Dirac_operator"] = "Staggered"
system["log_dir"] = "./logs"
system["Nthermalization"] = 10
system["update_method"] = "SLMC"
system["randomseed"] = 111
system["NC"] = 3
system["BoundaryCondition"] = [1, 1, 1, -1]
system["saveU_format"] = nothing
system["β"] = 5.7
actions["use_autogeneratedstaples"] = false
actions["couplingcoeff"] = Any[]
actions["couplinglist"] = Any[]
cg["eps"] = 1.0e-19
cg["MaxCGstep"] = 3000
staggered["Nf"] = 4
staggered["mass"] = 0.5
measurement["measurement_methods"] = Dict[]
measurement["measurement_methods"] = Dict[Dict{Any,Any}("eps" => 1.0e-19,"fermiontype" => "Staggered","Nf" => 4,"mass" => 0.5,"measure_every" => 5,"MaxCGstep" => 3000,"methodname" => "Chiral_condensate"), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Polyakov_loop"),
                                Dict{Any,Any}("eps_flow" => 0.01,"fermiontype" => nothing,"numflow" => 10,"measure_every" => 10,"methodname" => "Topological_charge","Nflowsteps" => 4),
                                        Dict{Any,Any}("eps" => 1.0e-19,"fermiontype" => "Wilson","r" => 1,"measure_every" => 5,"hop" => 0.141139,"methodname" => "Pion_correlator","MaxCGstep" => 3000), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Plaquette")]
#measurement["measurement_methods"] = Dict[Dict{Any,Any}("eps" => 1.0e-19,"fermiontype" => "Staggered","Nf" => 4,"mass" => 0.5,"measure_every" => 50,"MaxCGstep" => 3000,"methodname" => "Chiral_condensate"), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Polyakov_loop"), Dict{Any,Any}("fermiontype" => nothing,"numflow" => 10,"measure_every" => 50,"methodname" => "Topological_charge"), Dict{Any,Any}("eps" => 1.0e-19,"fermiontype" => "Wilson","r" => 1,"measure_every" => 50,"hop" => 0.141139,"methodname" => "Pion_correlator","MaxCGstep" => 3000), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Plaquette")]
measurement["measurement_dir"] = "SLMC_L04040404_beta5.7_Staggered_mass0.5"
measurement["measurement_basedir"] = "./measurements"
# - - - - - - - - - - - - - - - - - - -


text:
inf2nan(x) = x == Inf ? NaN : x

function loess_errors(mat)
    est = Array{Float64}(undef,0,size(mat)[2])
    qs = float.(1:41)
      for j in 1:4
          group = 4*(j-1)
          tot = vec(sum(mat[group+1:group+4,:], dims = 1))
          for i in 1:4
          if i != j
              errs  = mat[group+i,:]
              rlogp = log10.((errs .+1) ./tot)  # 1 psuedocount for each err, but if tot=0 will give NA
              rlogp .= inf2nan.(rlogp)
              model = loess(qs,rlogp)
              pred = Loess.predict(model, qs)
              R"df <- data.frame(q=$(qs), errs=$(errs), tot=$(tot), rlogp=$(rlogp))"
              R"mod.lo <- loess(rlogp ~ q, df, weights=tot)" ###!
          #        mod.lo <- loess(rlogp ~ q, df)
              pred = rcopy(R"predict(mod.lo, $(qs))")
              #println(ismissing.(pred)')
              pred[ismissing.(pred)] .= NaN
              maxrli = findlast(x-> !isnan(x),pred)
              minrli = findfirst(x-> !isnan(x),pred)
              pred[maxrli:end] .= pred[maxrli]
              pred[1:minrli] .= pred[minrli]

              est = vcat(est, 10 .^pred')
              end
          end
      end
      clamp!(est,1e-7,0.25) # editing out this line give closer output to R
      err = vcat(1 .-sum(est[1:3,:],dims = 1),
                 est[1:3,:],
                 est[4,:]',
                 1 .-sum(est[4:6,:],dims = 1),
                 est[5:6,:],
                 est[7:8,:],
                 1 .-sum(est[7:9,:],dims = 1),
                 est[9,:]',
                 est[10:12,:],
                 1 .-sum(est[10:12,:], dims = 1))
      return err
  end


text:
#Include KG_Module
include("../src/KG_Module.jl")

filepath1 = "./data/binary_run_1/binary_bh_rotating_radius_1.5_alpha_10.0_Rorbit_3.0.h5"
filepath2 = "./data/binary_run_1/binary_bh_rotating_radius_1.5_alpha_10.0_Rorbit_5.0.h5"

#Get field at certain positions
time_vector_1 , ene_vector_1 = get_energy( filepath1 , 10 , 15000.0)
time_vector_2 , ene_vector_2 = get_energy( filepath2 , 10 , 15000.0)

#Write to file
open("energy_bhb_rorbit_3.txt", "w") do io
    writedlm(io , [time_vector_1 ene_vector_1])
end

open("energy_bhb_rorbit_5.txt", "w") do io
    writedlm(io , [time_vector_2 ene_vector_2])
end


text:
### Some methods for velocity calculation for a system driven by a surface force
"""
Interface velocity in a viscous liquid where regions seperated by interface do have the same viscousity. If γ is passed the curvattureless algorithm is used to take into account force due to surface tension. Returns velocity field projected on vertex normals. Usefull for equilibrium calculations.
"""
function stokesvelocity(points,normals,faces,forcen,etaP,gammap)
    vareas = vertexareas(points,faces)
    velocityn = zeros(Float64,length(points))

    for xkey in 1:length(points)

        x = points[xkey]
        nx = normals[xkey]
        fx = forcen[xkey]

        s = 0
        for ykey in 1:length(points)
            if ykey==xkey
                continue
            end

            y = points[ykey]
            ny = normals[ykey]
            fy = forcen[ykey]

            ### I will need to check a missing 2
            s += vareas[ykey]*1 ./8/pi/etaP* dot(y-x,nx+ny)/norm(y-x)^3*(1-3*dot(y-x,nx)*dot(y-x,ny)/norm(y-x)^2) * gammap

            ### ????????
            s += vareas[ykey]*1 ./8/pi/etaP* ( dot(nx,ny)/norm(x-y) + dot(nx,x -y)*dot(ny,x-y)/norm(x-y)^3 )*(fy - fx)
        end
        velocityn[xkey] = s
    end
    return velocityn
end

function stokesvelocity(points,normals,faces,forcen,etaP)

    vareas = vertexareas(points,faces)
    velocityn = zeros(Float64,size(points,2))

    for xkey in 1:length(points)

        x = points[xkey]
        nx = normals[xkey]
        fx = forcen[xkey]

        s = 0
        for ykey in 1:length(points)
            if ykey==xkey
                continue
            end

            y = points[ykey]
            ny = normals[ykey]
            fy = forcen[ykey]

            s += vareas[ykey]*1 ./8/pi/etaP* ( dot(nx,ny)/norm(x-y) + dot(nx,x -y)*dot(ny,x-y)/norm(x-y)^3 )*(fy - fx)
        end

        velocityn[xkey] = s
    end

    return velocityn
end


text:
# global valiables are specified in this file
# See https://docs.julialang.org/en/v1/manual/variables-and-scoping/#Constants
# for more details about "const"

const n_mesh = 40
const emax = 2
const Nnmax= 20
const pmax_fm = 4
const chi_order = 3 #0:LO 1:NLO 2:NNLO 3:N3LO
const calc_NN = true
const calc_3N = false #density-dependent 3NF

"""
hbar-omega can be specify by hand
const hw = 20.0
or using formula
const hw = hw_formula(Anum,fnum)
#fnum=1: 41A^(-1/3) MeV
#funm=2: 45A^(-1/3)-25A^(-2/3) MeV (used for sd-shell)
"""
const Anum = 4
const hw = 20.0
#Anum = 28; const hw = hw_formula(Anum,2)
#const hw = hw_formula(Anum,1)

## SRG evolution (srg_lambda is in fm^{-1}
#const srg = false
const srg = true
const srg_lambda = 2.0

## file name and format for TBME
tx = "bare";if srg; tx ="srg"*string(srg_lambda);end
#const tbme_fmt = "snt"
const tbme_fmt = "myg"
const fn_tbme = "tbme_em500n3lo_"*tx*"hw"*string(round(Int64,hw))*"emax"*string(emax)*"."*tbme_fmt

const calc_monopole = false
const calc_std = false

## target {[n,l,j]} for TBMEs
const target_nlj=[] # no-core
#const target_nlj=[[1,0,1],[0,2,3],[0,2,5]] # sd-shell part
#const target_nlj=[[1,1,1],[1,1,3],[0,3,5],[0,3,7]] # pf-shell part

##for valence space operators
#-(v_chi_order=1): vsNLO
#-(v_chi_order=3): vsN3LO (not implemnted)
const v_chi_order = 0 # 0: free-space only
const n_mesh_P = 10
const Pmax_fm = 3

## constants (no need to modify basically...)
const jmax = 6
const lmax = jmax+1
const lcmax = jmax+1
const Mp = 938.272
const Mn = 939.5653
const Mm = (Mp+Mn)/2
const Ms = [Mp,Mm,Mn]
const Lambchi = 500 # cutoff Lambda
const itts = [-2,0,2]
const hc = 197.327053
const gA = -1.29
const Fpi = 92.4
const hc2 = hc^2
const hc3 = hc^3
const hc4 = hc^4
const gA2 = gA^2
const gA4 = gA^4
const mpis = [139.5702,134.9766,139.5702]
### Fermi momentum for density-dependent NN from 3NFs
const kF = 1.35


text:
abstract type EmissionShape end
abstract type EllipseEmission <: EmissionShape end
abstract type RectEmission <: EmissionShape end
abstract type PointEmission <: EmissionShape end

"""
    RayOriginGenerator{T<:Real}

Generates 3D points in world space which serve as origins for rays.
"""
abstract type RayOriginGenerator{T<:Real} end

abstract type AbstractRayGenerator{T<:Real} end

"""
    GeometricRayGenerator{T,O<:RayOriginGenerator{T}} <: AbstractRayGenerator{T}

Generates geometric [`Ray`](@ref)s according to the specific implementation of the subclass.
"""
abstract type GeometricRayGenerator{T,O<:RayOriginGenerator{T}} <: AbstractRayGenerator{T} end

"""
    OpticalRayGenerator{T} <: AbstractRayGenerator{T}

Generates [`OpticalRay`](@ref)s according to the specific implementation of the subclass.
"""
abstract type OpticalRayGenerator{T} <: AbstractRayGenerator{T} end
export AbstractRayGenerator, GeometricRayGenerator, OpticalRayGenerator

Base.iterate(a::AbstractRayGenerator, state = 1) = state > length(a) ? nothing : (generateray(a, state - 1), state + 1)
Base.getindex(a::AbstractRayGenerator, index) = generateray(a, index)
Base.firstindex(a::AbstractRayGenerator) = 0
Base.lastindex(a::AbstractRayGenerator) = length(a) - 1
Base.copy(a::AbstractRayGenerator) = a # most don't have any heap allocated stuff so don't really need copying

"""
    randinsolidangle(direction::SVector{3,T}, uvec::SVector{3,T}, vvec::SVector{3,T}, θmax::T)

Generates a unit vector pointing somewhere within the cone with half angle `θmax` around `direction`. `uvec` and `vvec` should be orthogonal to each other and `direction`.
"""
function randinsolidangle(direction::SVector{3,T}, uvec::SVector{3,T}, vvec::SVector{3,T}, θmax::T)::SVector{3,T} where {T<:Real}
    ϕ = rand(T) * 2π
    θ = NaNsafeacos(one(T) + rand(T) * (cos(θmax) - 1))
    return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * direction)
end

################################################################

function getuvvecs(direction::SVector{3,T}, rotationvec::SVector{3,T}) where {T<:Real}
    direction = normalize(direction)
    if abs(dot(rotationvec, direction)) == one(T)
        rotationvec = SVector{3,T}(1.0, 0.0, 0.0)
    end
    uvec = normalize(cross(normalize(rotationvec), direction))
    vvec = normalize(cross(direction, uvec))
    return direction, uvec, vvec
end

"""
    RandomRectOriginPoints{T} <: RayOriginGenerator{T}

Generates ray origins randomly within a rectangle.

```julia
RandomRectOriginPoints(numrays, halfsizeu, halfsizev; position = (0.0, 0.0, 0.0), direction = (0.0, 0.0, -1.0), rotationvec = (0.0, 1.0, 0.0))
```
"""
struct RandomRectOriginPoints{T} <: RayOriginGenerator{T}
    position::SVector{3,T}
    direction::SVector{3,T}
    uvec::SVector{3,T}
    vvec::SVector{3,T}
    halfsizeu::T
    halfsizev::T
    nrays::Int
    function RandomRectOriginPoints(numrays::Int, halfsizeu::T, halfsizev::T; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        new{T}(position, getuvvecs(direction, rotationvec)..., halfsizeu, halfsizev, numrays)
    end
end

"""
    GridRectOriginPoints{T} <: RayOriginGenerator{T}

Generates ray origins on a rectangular grid.

```julia
GridRectOriginPoints(numraysu, numraysv, halfsizeu, halfsizev; position = (0.0, 0.0, 0.0), direction = (0.0, 0.0, -1.0), rotationvec = (0.0, 1.0, 0.0))
```
"""
struct GridRectOriginPoints{T} <: RayOriginGenerator{T}
    position::SVector{3,T}
    direction::SVector{3,T}
    uvec::SVector{3,T}
    vvec::SVector{3,T}
    halfsizeu::T
    halfsizev::T
    nraysu::Int
    nraysv::Int
    function GridRectOriginPoints(numraysu::Int, numraysv::Int, halfsizeu::T, halfsizev::T; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        new{T}(position, getuvvecs(direction, rotationvec)..., halfsizeu, halfsizev, numraysu, numraysv)
    end
end

"""
    RandomEllipseOriginPoints{T} <: RayOriginGenerator{T}

Generates ray origins randomly across an ellipse.

```julia
RandomEllipseOriginPoints(numrays, halfsizeu, halfsizev; position = (0.0, 0.0, 0.0), direction = (0.0, 0.0, -1.0), rotationvec = (0.0, 1.0, 0.0))
```
"""
struct RandomEllipseOriginPoints{T} <: RayOriginGenerator{T}
    position::SVector{3,T}
    direction::SVector{3,T}
    uvec::SVector{3,T}
    vvec::SVector{3,T}
    halfsizeu::T
    halfsizev::T
    nrays::Int
    function RandomEllipseOriginPoints(numrays::Int, halfsizeu::T, halfsizev::T; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        new{T}(position, getuvvecs(direction, rotationvec)..., halfsizeu, halfsizev, numrays)
    end
end

"""
    HexapolarOriginPoints{T} <: RayOriginGenerator{T}

Generates ray origins in a hexapolar pattern.

```julia
HexapolarOriginPoints(nrings::Int, halfsizeu::T, halfsizev::T; position = (0.0, 0.0, 0.0), direction = (0.0, 0.0, -1.0), rotationvec = (0.0, 1.0, 0.0))
```
"""
struct HexapolarOriginPoints{T} <: RayOriginGenerator{T}
    position::SVector{3,T}
    direction::SVector{3,T}
    uvec::SVector{3,T}
    vvec::SVector{3,T}
    halfsizeu::T
    halfsizev::T
    nrings::Int
    function HexapolarOriginPoints(nrings::Int, halfsizeu::T, halfsizev::T; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        new{T}(position, getuvvecs(direction, rotationvec)..., halfsizeu, halfsizev, nrings)
    end
end

"""
    OriginPoint{T} <: RayOriginGenerator{T}

Single point origin for a source.

```julia
OriginPoint{T}(numrays; position = (0.0, 0.0, 0.0), direction = (0.0, 0.0, -1.0), rotationvec = (0.0, 1.0, 0.0))
```
"""
struct OriginPoint{T} <: RayOriginGenerator{T}
    position::SVector{3,T}
    direction::SVector{3,T}
    uvec::SVector{3,T}
    vvec::SVector{3,T}
    nrays::Int
    function OriginPoint{T}(numrays::Int; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
        new{T}(position, getuvvecs(direction, rotationvec)..., numrays)
    end
end

export RandomRectOriginPoints, GridRectOriginPoints, HexapolarOriginPoints, RandomEllipseOriginPoints, OriginPoint

Base.length(a::RayOriginGenerator) = a.nrays
Base.length(a::GridRectOriginPoints) = a.nraysu * a.nraysv
Base.length(a::HexapolarOriginPoints) = 1 + round(Int, (a.nrings * (a.nrings + 1) / 2) * 6)

direction(a::RayOriginGenerator) = a.direction
position(a::RayOriginGenerator) = a.position
uvec(a::RayOriginGenerator) = a.uvec
vvec(a::RayOriginGenerator) = a.vvec

"""
    genorigin(o::RayOriginGenerator{T}, n::Int) -> SVector{3,T}

Generate origin positions for rays based on the type of the generator, e.g., randomly within a rectangle or ellipse. `n` is the index of the point being generated, starting from 0.
This has little meaning for random generators, but is important for `HexapolarOriginPoints` and `GridRectOriginPoints`.
"""
genorigin(a::OriginPoint, ::Int) = a.position
function genorigin(a::RandomRectOriginPoints{T}, ::Int) where {T<:Real}
    u = 2 * rand(T) * a.halfsizeu - a.halfsizeu
    v = 2 * rand(T) * a.halfsizev - a.halfsizev
    return a.position + u * a.uvec + v * a.vvec
end
function genorigin(a::RandomEllipseOriginPoints{T}, ::Int) where {T<:Real}
    ϕ = rand(T) * 2π
    ρ = sqrt(rand(T))
    u = cos(ϕ) * a.halfsizeu
    v = sin(ϕ) * a.halfsizev
    return a.position + ρ * (u * a.uvec + v * a.vvec)
end
function genorigin(a::HexapolarOriginPoints{T}, n::Int) where {T<:Real}
    n = mod(n, length(a))
    if n == 0
        return a.position
    else
        t = 1
        ringi = 1
        for i in 1:(a.nrings)
            t += 6 * i
            if n < t
                ringi = i
                break
            end
        end
        ρ = ringi / a.nrings
        pind = n - (t - 6 * ringi)
        ϕ = (pind / (6 * ringi)) * 2π
        u = cos(ϕ) * a.halfsizeu
        v = sin(ϕ) * a.halfsizev
        return a.position + ρ * (u * a.uvec + v * a.vvec)
    end
end
function genorigin(a::GridRectOriginPoints{T}, n::Int) where {T<:Real}
    n = mod(n, length(a))
    v = a.nraysv == 1 ? zero(T) : 2 * Int(floor(n / a.nraysu)) / (a.nraysv - 1) - 1.0
    u = a.nraysu == 1 ? zero(T) : 2 * mod(n, a.nraysu) / (a.nraysu - 1) - 1.0
    return a.position + (a.halfsizeu * u * a.uvec) + (a.halfsizev * v * a.vvec)
end

################################################################

"""
    CollimatedSource{T,O} <: GeometricRayGenerator{T,O}

Source which generates collimated rays.

```julia
CollimatedSource(generator::RayOriginGenerator)
```
"""
struct CollimatedSource{T,O<:RayOriginGenerator{T}} <: GeometricRayGenerator{T,O}
    generator::O
    function CollimatedSource(generator::O) where {T<:Real,O<:RayOriginGenerator{T}}
        new{T,O}(generator)
    end
end

"""
    GridSource{T,O} <: GeometricRayGenerator{T,O}

Source which generates rays in directions which fall on an even grid on a rectangle subtended by angles halfangleu, halfanglev.

```julia
GridSource(generator::RayOriginGenerator, numraysu, numraysv, halfangleu, halfanglev)
```
"""
struct GridSource{T,O<:RayOriginGenerator{T}} <: GeometricRayGenerator{T,O}
    generator::O
    nraysu::Int
    nraysv::Int
    halfangleu::T
    halfanglev::T
    function GridSource(generator::O, numraysu::Int, numraysv::Int, halfangleu::T, halfanglev::T) where {T<:Real,O<:RayOriginGenerator{T}}
        new{T,O}(generator, numraysu, numraysv, halfangleu, halfanglev)
    end
end

"""
    RandomSource{T,O} <: GeometricRayGenerator{T,O}

Source which generates rays in directions sampled randomly from the solid angle ±θmax centred on the direction of the source.

```julia
RandomSource(generator::RayOriginGenerator, numrays = 1, θmax = T(π / 2))
```
"""
struct RandomSource{T,O<:RayOriginGenerator{T}} <: GeometricRayGenerator{T,O}
    generator::O
    nrays::Int
    θmax::T
    function RandomSource(generator::O, numrays::Int = 1, θmax::T = T(π / 2)) where {T<:Real,O<:RayOriginGenerator{T}}
        @assert θmax <= π
        new{T,O}(generator, numrays, θmax)
    end
end

export CollimatedSource, GridSource, RandomSource

Base.length(a::RandomSource) = a.nrays * length(a.generator)
Base.length(a::GridSource) = a.nraysu * a.nraysv * length(a.generator)
Base.length(a::GeometricRayGenerator) = length(a.generator)

direction(a::GeometricRayGenerator) = direction(a.generator)
position(a::GeometricRayGenerator) = position(a.generator)
uvec(a::GeometricRayGenerator) = uvec(a.generator)
vvec(a::GeometricRayGenerator) = vvec(a.generator)

"""
    gendirection(o::GeometricRayGenerator{T}, n::Int) -> SVector{3,T}

Generate directions for rays based on the type of the generator, e.g., randomly within a cone or collimated. `n` is the index of the point being generated, starting from 0.
This has little meaning for random generators, but is important for `GridSource`.
"""
gendirection(a::RandomSource, ::Int) = randinsolidangle(a.generator.direction, a.generator.uvec, a.generator.vvec, a.θmax)
gendirection(a::CollimatedSource, ::Int) = a.generator.direction
function gendirection(a::GridSource{T}, n::Int) where {T<:Real}
    # distributing evenly across the area of the rectangle which subtends the given angle (*not* evenly across the angle range)
    dindex = mod(n, a.nraysu * a.nraysv)
    v = a.nraysv == 1 ? zero(T) : 2 * Int(floor(dindex / a.nraysu)) / (a.nraysv - 1) - 1.0
    u = a.nraysu == 1 ? zero(T) : 2 * mod(dindex, a.nraysu) / (a.nraysu - 1) - 1.0
    θu = atan(u * tan(a.halfangleu) / 2) * a.halfangleu
    θv = atan(v * tan(a.halfanglev) / 2) * a.halfanglev
    dir = cos(θv) * (cos(θu) * direction(a) + sin(θu) * uvec(a)) + sin(θv) * vvec(a)
    return dir
end

"""
    generateray(o::GeometricRayGenerator{T}, n::Int) -> Ray{T,3}

Generate geometric rays distributed according to the type of the generator. `n` is the index of the point being generated, starting from 0.
This has little meaning for random generators, but is important for `GridSource`, for example.
"""
function generateray(a::GeometricRayGenerator{T,O}, n::Int) where {T<:Real,O}
    origin = genorigin(a.generator, n)::SVector{3,T}
    direction = gendirection(a, n)::SVector{3,T}
    return Ray(origin, direction)
end

origingen(a::GeometricRayGenerator) = a.generator

################################################################

"""
    UniformOpticalSource{T,O,P} <: OpticalRayGenerator{T}

Source of OpticalRays with uniform power.

```julia
UniformOpticalSource(generator::GeometricRayGenerator, centralwavelength, power = 1.0; sourcenum = 0)
```
"""
struct UniformOpticalSource{T,O,P<:GeometricRayGenerator{T,O}} <: OpticalRayGenerator{T}
    generator::P
    power::T
    centralwavelength::T
    sourcenum::Int

    function UniformOpticalSource(generator::P, centralwavelength::Unitful.Length, power::T = one(T); sourcenum::Int = 0) where {T<:Real,O<:RayOriginGenerator{T},P<:GeometricRayGenerator{T,O}}
        λ = Unitful.ustrip(Unitful.u"μm", centralwavelength, sourcenum)
        new{T,O,P}(generator, power, λ)
    end

    function UniformOpticalSource(generator::P, centralwavelength::T, power::T = one(T); sourcenum::Int = 0) where {T<:Real,O<:RayOriginGenerator{T},P<:GeometricRayGenerator{T,O}}
        new{T,O,P}(generator, power, centralwavelength, sourcenum)
    end
end

"""
    CosineOpticalSource{T,O,P} <: OpticalRayGenerator{T}

Source of OpticalRays with power defined by:
``I(\\theta) \\approx I_0(\\cos\\theta)^C``
Where ``\\theta`` is the angle of the ray to the central direction of the source, and ``C`` is the `cosineexp` parameter.

```julia
CosineOpticalSource(generator::GeometricRayGenerator, cosineexp, centralwavelength, power = 1.0; sourcenum = 0)
```
"""
struct CosineOpticalSource{T,O,P<:GeometricRayGenerator{T,O}} <: OpticalRayGenerator{T}
    generator::P
    power::T
    centralwavelength::T
    cosineexp::T
    sourcenum::Int

    function CosineOpticalSource(generator::P, cosineexp::T, centralwavelength::Unitful.Length, power::T = one(T); sourcenum::Int = 0) where {T<:Real,O<:RayOriginGenerator{T},P<:GeometricRayGenerator{T,O}}
        λ = Unitful.ustrip(Unitful.u"μm", centralwavelength)
        new{T,O,P}(generator, power, λ, cosineexp, sourcenum)
    end

    function CosineOpticalSource(generator::P, cosineexp::T, centralwavelength::T, power::T = one(T); sourcenum::Int = 0) where {T<:Real,O<:RayOriginGenerator{T},P<:GeometricRayGenerator{T,O}}
        new{T,O,P}(generator, power, centralwavelength, cosineexp, sourcenum)
    end
end

"""
    GaussianOpticalSource{T,O,P} <: OpticalRayGenerator{T}

Source of OpticalRays with power defined by:
``I(\\theta) \\approx I_0e^{-(G_ul^2 + G_vm^2)}``
Where ``l`` and ``m`` are the direction cosines in the u and v directions to the central direction of the source, and ``G_u`` and ``G_v`` and the `gaussianu` and `gaussianv` parameters.

```julia
GaussianOpticalSource(generator::GeometricRayGenerator, gaussianu, gaussianv, centralwavelength, power = 1.0; sourcenum = 0)
```
"""
struct GaussianOpticalSource{T,P<:GeometricRayGenerator{T}} <: OpticalRayGenerator{T}
    generator::P
    power::T
    centralwavelength::T
    gaussianu::T
    gaussianv::T
    sourcenum::Int

    function GaussianOpticalSource(generator::P, gaussianu::T, gaussianv::T, centralwavelength::Unitful.Length, power::T = one(T); sourcenum::Int = 0) where {T<:Real,P<:GeometricRayGenerator{T}}
        λ = Unitful.ustrip(Unitful.u"μm", centralwavelength)
        new{T,P}(generator, power, λ, gaussianu, gaussianv, sourcenum)
    end

    function GaussianOpticalSource(generator::P, gaussianu::T, gaussianv::T, centralwavelength::T, power::T = one(T); sourcenum::Int = 0) where {T<:Real,P<:GeometricRayGenerator{T}}
        new{T,P}(generator, power, centralwavelength, gaussianu, gaussianv, sourcenum)
    end
end

export UniformOpticalSource, CosineOpticalSource, GaussianOpticalSource

function generateray(a::GaussianOpticalSource{T,P}, n::Int) where {T<:Real,P}
    r = generateray(a.generator, n)::Ray{T,3}
    l = dot(direction(r), uvec(a.generator))
    m = dot(direction(r), vvec(a.generator))
    power = a.power * exp(-(a.gaussianu * l^2 + a.gaussianv * m^2))
    return OpticalRay(r, power, a.centralwavelength, sourcenum = a.sourcenum)
end

function generateray(a::CosineOpticalSource{T,P}, n::Int) where {T<:Real,P}
    r = generateray(a.generator, n)::Ray{T,3}
    cosanglebetween = dot(direction(r), direction(a.generator))
    power = a.power * cosanglebetween^a.cosineexp
    return OpticalRay(r, power, a.centralwavelength, sourcenum = a.sourcenum)
end

"""
    generateray(o::OpticalRayGenerator{T}, n::Int) -> OpticalRay{T,3}

Generate optical rays distributed according to the type of the generator. `n` is the index of the point being generated, starting from 0.
This has little meaning for random generators, but is important for generators using `GridSource` or `GridRectOriginPoints`, for example.
"""
function generateray(a::UniformOpticalSource{T,P}, n::Int) where {T<:Real,P}
    r = generateray(a.generator, n)::Ray{T,3}
    return OpticalRay(r, one(T), a.centralwavelength, sourcenum = a.sourcenum)
end

Base.length(a::OpticalRayGenerator) = length(a.generator)
origingen(a::OpticalRayGenerator) = origingen(a.generator)
direction(a::OpticalRayGenerator) = direction(a.generator)
position(a::OpticalRayGenerator) = position(a.generator)
uvec(a::OpticalRayGenerator) = uvec(a.generator)
vvec(a::OpticalRayGenerator) = vvec(a.generator)

################################################################################################################################

"""
    PixelSource{T,C} <: RayGenerator{T}

Ray generator which encapsulates a number of subpixels, rays are generated in each subpixel in an interleaved manner.
The subpixels must be positioned correctly relative to each other when input to the constructor.

All subpixels must be coplanar, the orientation of the pixel and any display made using the pixel is taken from the first subpixel.
All subpixels must have the same number of rays (if we make subpixels be handled sequentially we wouldn't need this).
All subpixels must be of the same type.

```julia
PixelSource(subpixels::Vector{OpticalRayGenerator}; colormap = nothing, position = (0, 0, 0))
PixelSource(subpixels::Vararg{OpticalRayGenerator}; colormap = nothing, position = (0, 0, 0))
```
"""
struct PixelSource{T,P<:OpticalRayGenerator{T},C} <: OpticalRayGenerator{T}
    subpixels::Vector{P}
    colormap::SVector{C,Int}
    position::SVector{3,T}

    function PixelSource(subpixels::Vararg{P}; colormap = nothing, position::SVector{3,T} = SVector{3,T}(0, 0, 0)) where {T<:Real,P<:OpticalRayGenerator{T}}
        PixelSource(collect(subpixels), colormap = colormap, position = position)
    end

    function PixelSource(subpixels::Vector{P}; colormap = nothing, position::SVector{3,T} = SVector{3,T}(0, 0, 0)) where {T<:Real,P<:OpticalRayGenerator{T}}
        C = length(subpixels)
        # we need to map from subpixel index to color in image (e.g. for RGB pentile where there are multiple subpixels for one color)
        if (colormap === nothing)
            colormap = SVector{C,Int}(1:C)
        end
        for s in subpixels[1:end]
            @assert uvec(s) == uvec(subpixels[1]) && vvec(s) == vvec(subpixels[1])
            @assert length(s) == length(subpixels[1])
        end
        new{T,P,C}(subpixels, colormap, position)
    end
end
export PixelSource

"""
    generateray(a::PixelSource{T}, n::Int) -> OpticalRay{T,3}

Generates optical rays from all subpixels in the pixel. One ray is generated from each subpixel sequentially before looping back to the start.
"""
function generateray(a::PixelSource{T,P,C}, n::Int) where {T<:Real,P,C}
    subpixelindex = mod(n, C) + 1
    subpixeliter = Int(floor(n / C))
    ray = generateray(a.subpixels[subpixelindex], subpixeliter)::OpticalRay{T,3}
    return OpticalRay(origin(ray) + a.position, direction(ray), power(ray), wavelength(ray), sourcenum = sourcenum(ray))
end

uvec(a::PixelSource) = uvec(a.subpixels[1])
vvec(a::PixelSource) = vvec(a.subpixels[1])
colorindex(a::PixelSource{T,P,C}, i::Int) where {T<:Real,P,C} = a.colormap[mod(i, C) + 1]
Base.length(a::PixelSource) = sum(length.(a.subpixels))

"""
    OpticalSourceArray{T} <: RayGenerator{T}

Generates rays from an array of the given source at the specified locations.

```julia
OpticalSourceArray(generator::OpticalRayGenerator, positions::Vector{SVector{3,T}})
```
"""
struct OpticalSourceArray{T,S<:OpticalRayGenerator{T}} <: OpticalRayGenerator{T}
    generator::S
    positions::Vector{SVector{3,T}}

    function OpticalSourceArray(generator::S, positions::Vector{SVector{3,T}}) where {T<:Real,S<:OpticalRayGenerator{T}}
        new{T,S}(generator, positions)
    end
end
export OpticalSourceArray

Base.copy(a::OpticalSourceArray) = OpticalSourceArray(copy(a.generator), copy(a.positions))
Base.length(a::OpticalSourceArray) = length(a.generator) * length(a.positions)
Base.show(io::IO, a::OpticalSourceArray) = print(io, "OpticalSourceArray($(a.generator), $(length(a.positions)) positions)")

"""
    generateray(a::OpticalSourceArray{T}, n::Int) -> OpticalRay{T,3}

Generates optical rays from all generators in the array. One ray is generated from each element sequentially before looping back to the start of the array.
"""
function generateray(a::OpticalSourceArray{T,S}, n::Int) where {T<:Real,S}
    position = mod(n, length(a.positions)) + 1
    positioniter = Int(floor(n / length(a.positions)))
    ray = generateray(a.generator::S, positioniter)::OpticalRay{T,3}
    return OpticalRay(origin(ray) + a.positions[position], direction(ray), power(ray), wavelength(ray), sourcenum = sourcenum(ray))
end

"""
    BasicDisplayPanel{T,C} <: RayGenerator{T}

Ray generator representing a simple panel display. The panel is flat and pixles are on a regular rectangular grid. Each pixel corresponds to one image pixel.

```julia
BasicDisplayPanel(pixel::PixelSource, pitchx, pitchy, image, position = (0, 0, 0))
BasicDisplayPanel(generator::OpticalSourceArray{T,PixelSource}, pixelvals::Vector{Union{T,Vector{T}}})
```
"""
struct BasicDisplayPanel{T,P<:PixelSource{T}} <: OpticalRayGenerator{T}
    generator::OpticalSourceArray{T,P}
    pixelvals::Vector{Union{T,Vector{T}}}

    function BasicDisplayPanel(pixel::P, pitchx::T, pitchy::T, image::Array, position::SVector{3,T} = SVector{3,T}(0, 0, 0), gamma::T = 2.2) where {T<:Real,P<:PixelSource{T}}
        @assert pitchx > 0 && pitchy > 0
        h, w = size(image)
        positions = Vector{SVector{3,T}}(undef, 0)
        pixelvals = []
        for colindex in 1:w
            for rowindex in 1:h
                if size(channelview(image))[1] == h
                    k = T(channelview(image)[rowindex, colindex])^gamma
                else
                    k = T.(channelview(image)[:, rowindex, colindex]) .^ gamma
                end
                if sum(k) > 0.001
                    o = position + (colindex - (w + 1) / 2) * pitchx * uvec(pixel) + (rowindex - (h + 1) / 2) * pitchy * vvec(pixel)
                    push!(positions, o)
                    push!(pixelvals, k)
                end
            end
        end
        new{T,P}(OpticalSourceArray(pixel, positions), pixelvals)
    end

    function BasicDisplayPanel(generator::OpticalSourceArray{T,P}, pixelvals::Vector{Union{T,Vector{T}}}) where {T<:Real,P<:PixelSource{T}}
        new{T,P}(generator, pixelvals)
    end
end

export BasicDisplayPanel

Base.copy(a::BasicDisplayPanel) = BasicDisplayPanel(copy(a.generator), copy(a.pixelvals))
Base.length(a::BasicDisplayPanel) = length(a.generator)
Base.show(io::IO, a::BasicDisplayPanel) = print(io, "BasicDisplayPanel($(a.generator))")

"""
    generateray(a::BasicDisplayPanel{T}, n::Int) -> OpticalRay{T,3}

Generates optical rays from all pixels in the display. One ray is generated from each pixel sequentially before looping back to the start of the display.
"""
function generateray(a::BasicDisplayPanel{T,P}, n::Int) where {T<:Real,P}
    r = generateray(a.generator, n)::OpticalRay{T,3}
    position = mod(n, length(a.pixelvals)) + 1
    positioniter = Int(floor(n / length(a.pixelvals)))
    k = a.pixelvals[position]
    if length(k) == 1
        pow = power(r) * k
    else
        pow = power(r) * k[colorindex(a.generator.generator::P, positioniter)]
    end
    return OpticalRay(ray(r), pow, wavelength(r), sourcenum = sourcenum(r))
end

uvec(a::BasicDisplayPanel) = uvec(a.generator)
vvec(a::BasicDisplayPanel) = vvec(a.generator)

"""
    OpticalSourceGroup{T} <: RayGenerator{T}

Wrapper to group a number of separate optical sources into a single iterable object.
If the sources of the same type then performance will be improved.

```julia
OpticalSourceGroup(generators::Vector{<:OpticalRayGenerator})
OpticalSourceGroup(generators::Vararg{OpticalRayGenerator})
```
"""
abstract type OpticalSourceGroup{T} <: OpticalRayGenerator{T} end

struct OpticalSourceGroupG{T} <: OpticalSourceGroup{T}
    generators::Vector{OpticalRayGenerator{T}}
    rayspergen::Vector{Int}
end

struct OpticalSourceGroupS{T,S<:OpticalRayGenerator{T}} <: OpticalSourceGroup{T}
    generators::Vector{S}
    rayspergen::Vector{Int}
end

Base.show(io::IO, ::Type{<:OpticalSourceGroup}) = print(io, "OpticalSourceGroup")
Base.show(io::IO, a::OpticalSourceGroupS{T}) where {T<:Real} = print(io, "OpticalSourceGroupS{$T}($(length(a.generators)), $(a.rayspergen))")
Base.show(io::IO, a::OpticalSourceGroupG{T}) where {T<:Real} = print(io, "OpticalSourceGroupG{$T}($(length(a.generators)), $(a.rayspergen))")

OpticalSourceGroup(generators::Vararg{OpticalRayGenerator{T}}) where {T<:Real} = OpticalSourceGroup(collect(generators))
function OpticalSourceGroup(generators::Vector{<:OpticalRayGenerator{T}}) where {T<:Real}
    if length(generators) == 1
        return generators[1]
    end
    S = typeof(generators[1])
    for a in generators[2:end]
        if typeof(a) != S
            return OpticalSourceGroupG{T}(generators, length.(generators))
        end
    end
    OpticalSourceGroupS{T,S}(generators, length.(generators))
end
export OpticalSourceGroup

Base.copy(a::OpticalSourceGroup) = OpticalSourceGroup(a.generators)
Base.length(a::OpticalSourceGroup) = sum(a.rayspergen)

"""
    generateray(a::OpticalSourceGroup{T}, n::Int) -> OpticalRay{T,3}

Generate optical rays for each source in the group. All rays are generated for the first source, then all for the second source and so on as `n` increases.
"""
function generateray(a::OpticalSourceGroup{T}, n::Int) where {T<:Real}
    t = 0
    for (i, k) in enumerate(a.rayspergen)
        if n < t + k
            return generateray(a.generators[i], n - t)::OpticalRay{T,3}
            break
        end
        t += k
    end
    return nothing
end

"""
    RayListSource{T} <: OpticalRayGenerator{T}

Ray generator constructed manually from a list of rays which just outputs those rays in order.

```julia
RayListSource(rays::Vararg{OpticalRay})
RayListSource(rays::Vector{OpticalRay})
```
"""
struct RayListSource{T} <: OpticalRayGenerator{T}
    rays::Vector{OpticalRay{T,3}}
    RayListSource(rays::Vararg{OpticalRay{T,3}}) where {T<:Real} = new{T}(collect(rays))
    RayListSource(rays::Vector{OpticalRay{T,3}}) where {T<:Real} = new{T}(rays)
end
export RayListSource

Base.length(a::RayListSource) = length(a.rays)
origingen(::RayListSource) = nothing

function generateray(a::RayListSource, n::Int)
    return a.rays[n + 1]
end


text:
## PART 2 -- Note: you have to look at p1.jl for this to make sense

# Find integrated autocorrelation times

include("../q1/autocorrelation.jl");

println("\nPart 2");
println("\n\tFirst, plot autocorrelation to estimate τint");

# Fuse md trials
md = [md_1069 md_1304]

# Number of quantities to look at
N = 6

# Same (almost) as q1p3:
lengthmd = size(md)[1];

# Create x and y arrays for plotting
xcorr = [0:1:100];
ycorr = zeros( length( xcorr ), N );
for a in [1:N]
  # Calculate autocorrelations
  ycorr[:,a] = autocorrelation(lengthmd, md, a, xcorr);
  # Normalize by dividing through by Ĉ[1,a]
  ycorr[:,a] /= ycorr[1,a];
end

# Make plots for each alpha
println("\n\tMaking autocorrelation plots and placing in 'mdautocorr' vector");
using Gadfly
Gadfly.set_default_plot_size(20Gadfly.cm, 12Gadfly.cm)
mdautocorr = Plot[];
quantity = ["Temperature at T = 1.069", "Potential Energy at T = 1.069",
            "Virial at T = 1.069", "Temperature at T = 1.304",
            "Potential Energy at T = 1.304", "Virial at T = 1.304"];

for i in [1:N]
  # Plot
  push!(mdautocorr, Gadfly.plot ( x = 10xcorr, y = ycorr[:,i], Geom.line, Geom.point,
      Guide.xlabel("Step Separation"), Guide.ylabel("Autocorrelation"),
      Guide.title(join(["Autocorrelation for ",quantity[i]]," ")),
      Guide.xticks(ticks=[0:100:1000])))
end

# Find integrated autocorrelation times τint (similar to q1p4)

# Pick ncut for temperature, potential energy, and the virial
#   based on observed graph of autocorrelation
ncutt_1069 = 250;
ncutp_1069 = 250;
ncutv_1069 = 250;
ncutt_1304 = 100;
ncutp_1304 = 100;
ncutv_1304 = 100;

τtemp_1069 = τ(lengthmd, temps_1069, ncutt_1069);
τpe_1069 = τ(lengthmd, pe_1069, ncutp_1069);
τvir_1069 = τ(lengthmd, vir_1069, ncutv_1069);
τtemp_1304 = τ(lengthmd, temps_1304, ncutt_1304);
τpe_1304 = τ(lengthmd, pe_1304, ncutp_1304);
τvir_1304 = τ(lengthmd, vir_1304, ncutv_1304);

τmd = [τtemp_1069, τpe_1069, τvir_1069, τtemp_1304, τpe_1304, τvir_1304]

println("\n\tAutocorrelation times:\n");
for i in [1:N]
  println("\t\tτ̂v",quantity[i]":\t",τmd[i]);
end


text:
module WiringDiagrams

using Reexport

include("Core.jl")
include("Layers.jl")
include("Algebraic.jl")
include("Algorithms.jl")
include("Expressions.jl")
include("Serialization.jl")
include("GraphML.jl")
include("JSON.jl")

@reexport using .WiringDiagramCore
@reexport using .WiringLayers
@reexport using .AlgebraicWiringDiagrams
@reexport using .WiringDiagramAlgorithms
@reexport using .WiringDiagramExpressions

using .WiringDiagramSerialization
export convert_from_graph_data, convert_to_graph_data
@reexport using .GraphMLWiringDiagrams
@reexport using .JSONWiringDiagrams

end


text:
#mutable struct definitions
abstract type Stem end

mutable struct CancerStemCell <: Stem
    mutations::Array{Int64,1}
    fitness::Int64
end

mutable struct CancerNonStemCell <: Stem
    mutations::Array{Int64,1}
    fitness::Int64
    Ndivs::Int64
end

mutable struct CancerCells
  stemcells::Array{CancerStemCell, 1}
  nonstemcells::Array{CancerNonStemCell, 1}
  ncells::Array{Int64, 1}
  α::Float64
  d::Float64
end

mutable struct RawResults
  cells::CancerCells
  Nvec::Array{Float64, 1}
  divisions::Array{Int64, 1}
end

mutable struct cancercell
    mutations::Array{Int64,1}
    fitness::Int64
end

mutable struct RawOutput
    Nvec::Array{Int64, 1}
    tvec::Array{Float64, 1}
    muts::Array{Int64, 1}
    cells::Array{cancercell, 1}
    birthrates::Array{Float64, 1}
    deathrates::Array{Float64, 1}
    clonetype::Array{Int64, 1}
    clonetime::Array{Float64, 1}
    subclonemutations::Array{Any, 1}
    cloneN::Array{Int64, 1}
    Ndivisions::Array{Int64, 1}
    aveDivisions::Array{Float64, 1}
end

mutable struct bdprocess
  N::Array{Int64, 1}
  t::Array{Float64, 1}
  clonefreq::Array{Float64, 1}
end

mutable struct SimResult
    clonefreq::Array{Float64,1}
    clonefreqp::Array{Float64,1}
    clonetime::Array{Float64,1}
    subclonemutations::Array{Int64,1}
    birthrates::Array{Float64,1}
    deathrates::Array{Float64,1}
    tend::Float64
    trueVAF::Array{Float64,1}
    cloneN::Array{Int64, 1}
    clonetype::Array{Int64, 1}
    Ndivisions::Array{Int64, 1}
    cells::Array{cancercell, 1}
    aveDivisions::Array{Float64, 1}
end

mutable struct InputParameters
    numclones::Int64
    Nmax::Int64
    detectionlimit::Float64
    ploidy::Int64
    read_depth::Float64
    clonalmutations::Int64
    selection::Array{Float64,1}
    μ::Float64
    b::Float64
    d::Float64
    tevent::Array{Float64,1}
    ρ::Float64
    cellularity::Float64
    fixedmu::Bool
    timefunction::Function
    maxclonefreq::Int64
end

mutable struct StemCellSimResult
    trueVAF::Array{Float64,1}
    cells::CancerCells
    N::Array{Float64, 1}
    divisions::Array{Float64, 1}
    stemcellfraction::Array{Float64,1}
end

###############################################################################

function newmutations(cancercell, μ, mutID)
    #function to add new mutations to cells based on μ
    numbermutations= 1
    cancercell.mutations = append!(cancercell.mutations, mutID:mutID+numbermutations-1)
    mutID = mutID + numbermutations
    return cancercell, mutID
end

function newmutationsinit(cancercell, μ, mutID)
    numbermutations = 0
    cancercell.mutations = append!(cancercell.mutations,mutID:mutID+numbermutations-1)
    mutID = mutID + numbermutations
    return cancercell, mutID
end

function initializesim(clonalmutations)

    #initialize time to zero
    t = 0.0
    tvec = Float64[]
    push!(tvec,t)

    #population starts with one cell
    N = 1
    Nvec = Int64[]
    push!(Nvec,N)

    #Initialize array of cell type that stores mutations for each cell and their fitness type
    #fitness type of 1 is the host population, lowest fitness
    cells = cancercell[]
    #sizehint!(cells, 10^5)
    push!(cells,cancercell([],1))

    #need to keep track of mutations, assuming infinite sites, new mutations will be unique,
    #we assign each new muation with a unique integer by simply counting up from one
    mutID = 1
    cells[1],mutID = newmutationsinit(cells[1],clonalmutations,mutID)

    return t,tvec,N,Nvec,cells,mutID
end

function birthdeathprocess(b, d, Nmax; nclones = 0, s = repeat([1.0], inner = nclones), tevent = collect(1.0:0.5:100.0)[1:nclones], timefunction::Function = exptime() =  - log(rand()))

  simout = tumourgrow_birthdeath(b, d, Nmax, 0.0; numclones = nclones, clonalmutations = 0.0, s = s, tevent = tevent, timefunction = timefunction)
  M, fitness = cellsconvert(simout.cells)
  pctfit=Float64[]
  for i in 1:nclones push!(pctfit,sum(fitness.==(i+1))/Nmax) end

  return bdprocess(simout.Nvec, simout.tvec, pctfit)
end

exptime() = - log(rand())
meantime() = 1

function copycell(cancercellold::cancercell)
  newcancercell::cancercell = cancercell(copy(cancercellold.mutations), copy(cancercellold.fitness))
end

function tumourgrow_birthdeath(b, d, Nmax, μ; numclones=1, clonalmutations = μ, s = [0.0], tevent=[0.0], maxclonefreq = 200, timefunction::Function = meantime)

    #set array of birthrates
    birthrates = [b]
    deathrates = [d]
    times = vcat(tevent, 0.0)
    #time is defined in terms of population doublings
    timesN = round.(Int64, vcat(exp.(log(2).*times[1:end-1]), 0.0))

    #depending on number of clones add birthrates to model, fitness is randomly distributed between death and birth rates
    for i in 1:numclones
        push!(deathrates, rand() * deathrates[1])
        push!(birthrates,(1 + s[i]) * (birthrates[1] - deathrates[1]) + deathrates[i + 1])
    end

    #Rmax starts with b + d and changes once a fitter mutant is introduced, this ensures that
    # b and d have correct units
    Rmax = b + d

    #initialize arrays and parameters
    t,tvec,N,Nvec,cells,mutID = initializesim(clonalmutations)
    muts = Int64[]
    push!(muts,mutID)

    #we only want to introduce mutant once so have variable that keeps track of how many mutants have been introduced, keep track of which type of which clone aquires new clone
    fitmutant = 1
    clonetype = Int64[]
    clonetime = Float64[]
    subclonemutations = Any[]
    cloneN = Int64[]
    Ndivisions = Int64[]
    aveDivisions = Float64[]

    clonefreq = zeros(Int64, numclones + 1)
    clonefreq[1] = 1

    executed = false
    changemutrate = broadcast(!, BitArray(undef, numclones + 1))

    while N < Nmax

        #pick a random cell
        randcell = rand(1:N)
        r = rand(Uniform(0,Rmax))
	     Nt = N

        #birth event if r<birthrate, access correct birthrate from cells array
        if r < birthrates[cells[randcell].fitness]

            #population increases by one
            N = N + 1
            #copy cell and mutations for cell that reproduces
            #push!(cells, deepcopy(cells[randcell]))
            push!(cells, copycell(cells[randcell]))
            #add new mutations to both new cells
            if μ > 0.0
              cells[randcell],mutID = newmutations(cells[randcell],μ,mutID)
              cells[end],mutID = newmutations(cells[end],μ,mutID)
            end
            push!(muts,mutID)
            clonefreq[cells[randcell].fitness] = clonefreq[cells[randcell].fitness] + 1
            push!(Nvec, N)
            Δt =  1/(Rmax * Nt) .* timefunction()
            t = t + Δt
            push!(tvec,t)

            #if population time is tevent, cell is mutated into fitter cell
            if N >= timesN[fitmutant]
                if fitmutant != numclones + 1
                    #one mutant turns into another "type" so decreases in frequency

                    clonefreq[cells[randcell].fitness] = clonefreq[cells[randcell].fitness] - 1
                    #keep track of how many clones
                    fitmutant += 1
                    push!(clonetype, cells[randcell].fitness)
                    #change one mutant to fitter mutant
                    cells[randcell].fitness = fitmutant
                    #new type increases in frequency
                    clonefreq[cells[randcell].fitness] = clonefreq[cells[randcell].fitness] + 1

                    #change Rmax given that we now have a new fitter mutant
                    Rmax = maximum(birthrates[1:fitmutant]) + maximum(deathrates[1:fitmutant])

                    push!(clonetime, t)
                    push!(subclonemutations, deepcopy(cells[randcell].mutations))
                    push!(cloneN, N)
                    push!(Ndivisions, length(cells[randcell].mutations))
                    divs = map(x -> length(x.mutations), cells)
                    push!(aveDivisions, mean(divs))

                end
            end
        end

        if (birthrates[cells[randcell].fitness] + deathrates[cells[randcell].fitness]) <= r
          push!(Nvec, N)
          Δt =  1/(Rmax * Nt) * timefunction()
          t = t + Δt
          push!(tvec,t)
        end

        #death event if b<r<b+d
        if (birthrates[cells[randcell].fitness] <= r < birthrates[cells[randcell].fitness] + deathrates[cells[randcell].fitness])

            #population decreases by 1
            N = N - 1
            #frequency of cell type decreases
            clonefreq[cells[randcell].fitness] = clonefreq[cells[randcell].fitness] - 1
            #remove deleted cell
            deleteat!(cells,randcell)
            push!(Nvec,N)
            Δt =  1/(Rmax * Nt) * timefunction()
            t = t + Δt
            push!(tvec,t)
        end

        #every cell dies reinitialize simulation
        if (N == 0)
            t,tvec,N,Nvec,cells,mutID = initializesim(clonalmutations)
            muts = Int64[]
            push!(muts,mutID)
        end

        if (executed == false) && ((clonefreq.>maxclonefreq) == changemutrate)
            #if population of all clones is sufficiently large no new mutations
            #are acquired, can use this approximation as only mutations above 1%
            # frequency can be reliably detected
            μ = 0
            executed = true
        end

    end
    return RawOutput(Nvec, tvec, muts, cells, birthrates, deathrates, clonetype, clonetime, subclonemutations, cloneN, Ndivisions, aveDivisions)
end

function cellsconvert(cells)
    #convert from array of cell types to one array with mutations and one array with cell fitness

    fitness = zeros(Int64,length(cells))
    mutations = Int64[]
    sizehint!(mutations, length(cells) * 10)

    for i in 1:length(cells)
        append!(mutations,cells[i].mutations)
        fitness[i] = cells[i].fitness
    end

    return mutations, fitness
end

function allelefreq(mutations, cellnum)
    #create dictionary that maps mutation ID to allele frequency
    f = counts(mutations,minimum(mutations):maximum(mutations))
    muts = collect(minimum(mutations):maximum(mutations))
    idx = f .> 0.01
    f = map(Float64, f[idx])
    muts = muts[idx]
    Dict{Int64, Float64}(muts[i]::Int64 => f[i]::Float64 for i in 1:length(f))
end

function getresults(tevent::Array{Float64, 1}, s::Array{Float64, 1}, b, d, μ, Nmax; ploidy = 2, clonalmutations = 100, nc = 0, timefunction = exptime, maxclonefreq = 200)

    #Nvec,tvec,mvec,cells,br,dr,ct,clonetime
    sresult = tumourgrow_birthdeath(b, d, Nmax, μ; numclones = nc, s = s, tevent = tevent, clonalmutations = 0, timefunction = timefunction, maxclonefreq = maxclonefreq);
    M,fitness = cellsconvert(sresult.cells)

    return M, fitness, sresult.tvec[end], sresult.clonetime, sresult.subclonemutations, sresult.birthrates, sresult.deathrates, sresult.cloneN, sresult.clonetype, sresult.Ndivisions, sresult.cells, sresult.aveDivisions

end

function allelefreqexpand(AFDict, μ, subclonemutations; fixedmu = false)

  #expand allele frequency given mutation rate and calculate number of mutations in the subclones
  #subclonemutations = convert(Array{Array{Int64,1},1}, subclonemutations)
  if fixedmu == false
    cmuts = zeros(Int64, length(subclonemutations))
    mutfreqs = collect(values(AFDict))
    mutids = collect(keys(AFDict))
    mutations = rand(Poisson(μ), length(mutfreqs))
    AFnew = zeros(Int64, sum(mutations))

    for i in 1:length(cmuts)
      idx = findall((in)(mutids), subclonemutations[i])
      cmuts[i] = sum(mutations[idx])
    end

    j = 0
    for f in 1:length(mutfreqs)
        AFnew[(j + 1): j + mutations[f]] = fill(mutfreqs[f], mutations[f])
        j = j + mutations[f]
    end

  else
    mutfreqs = collect(values(AFDict))
    mutids = collect(keys(AFDict))
    μint = round(Int64, μ)
    mutations = fill(Int64(μ), length(mutfreqs))
    AFnew = zeros(Int64, sum(mutations))
    cmuts = zeros(Int64, length(subclonemutations))

    for i in 1:length(cmuts)
      idx = findall((in)(mutids), subclonemutations[i])
      cmuts[i] = sum(mutations[idx])
    end

    j = 0
    for f in 1:length(mutfreqs)
      AFnew[(j + 1): j + mutations[f]] = fill(mutfreqs[f], mutations[f])
      j = j + mutations[f]
    end
  end

  return AFnew, cmuts
end

function calculateclonefreq(pctfit, cmuts, clonetype)

  cf = deepcopy(pctfit)
  parent = clonetype
  cloneid = collect(1:length(clonetype))
  prev = 0
  for i in reverse(cloneid)[1:end-1]
    if parent[i] == 0
        continue
    end
    cf[parent[i]] = cf[parent[i]] + cf[i]
    cmuts[i] = cmuts[i] - cmuts[parent[i]]
  end

  return cf, cmuts
end

function run1simulation(IP::InputParameters, minclonesize, maxclonesize)

    M, fitness, tend, clonetime, subclonemutations, br, dr, cloneN, clonetype, Ndivisions, cells, aveDivisions = getresults(IP.tevent, IP.selection, IP.b, IP.d, IP.μ, IP.Nmax; ploidy = IP.ploidy, clonalmutations = IP.clonalmutations, nc = IP.numclones, timefunction = IP.timefunction, maxclonefreq = IP.maxclonefreq)

    if length(clonetime)!= IP.numclones

        IP.numclones = length(clonetime)
        IP.tevent = IP.tevent[1:IP.numclones]
        IP.selection = IP.selection[1:IP.numclones]
        br = br[1:IP.numclones]
        dr = dr[1:IP.numclones]
        clonetype = clonetype[1:IP.numclones]

    end

    AF = allelefreq(M, IP.Nmax)
    AF, cmuts = allelefreqexpand(AF, IP.μ, subclonemutations, fixedmu = IP.fixedmu)
    #prepend!(AF, repeat([Float64(IP.Nmax)], inner = IP.clonalmutations))
    prepend!(AF, fill(Float64(IP.Nmax), IP.clonalmutations))

    pctfit=Float64[]
    for i in 1:IP.numclones push!(pctfit,sum(fitness.==(i+1))/IP.Nmax) end

    #remove clones that have frequency < detectionlimit
#    detectableclones = (pctfit.>(IP.detectionlimit)) & (pctfit.<0.95)
    if length(pctfit) > 1
      clonefreq, cmuts = calculateclonefreq(pctfit, cmuts, clonetype .- 1)
      !(sum(clonefreq.>1.0) > 0) || error("There is a clone with frequency greater than 1, this should be impossible ($(clonefreq)), $(clonetype), $(pctfit)")
    else
      clonefreq = pctfit
    end

    if VERSION < v"0.6-"
        detectableclones = (clonefreq.>minclonesize) & (clonefreq.<maxclonesize)  # Deprecated as of v0.6
    else
        detectableclones = (clonefreq.>minclonesize) .& (clonefreq.<maxclonesize)
    end

    clonefreq = clonefreq[detectableclones]

    if sum(detectableclones) != IP.numclones
        IP.numclones = sum(detectableclones)
        IP.tevent = IP.tevent[detectableclones]
        IP.selection = IP.selection[detectableclones]
        clonetype = clonetype[detectableclones]
        pushfirst!(detectableclones, true)
        detectableclones = detectableclones[1:length(br)]
        br = br[detectableclones]
        dr = dr[detectableclones]
    end

    #return SimResults object
    return SimResult(clonefreq, pctfit, clonetime, cmuts, br, dr, tend, AF, cloneN, clonetype .- 1, Ndivisions, cells, aveDivisions), IP
end

###############################################################################


text:
__precompile__()

module BackgroundSegmenter
export
    ActiveQueue,
    apply,
    apply!,
    apply_mrf!,
    bg_energy!,
    cluster,
    components,
    components!,
    delete!,
    dequeue!,
    enqueue!,
    Filter,
    filter_components,
    filter_components!,
    front,
    Gaussian,
    GridGraph,
    HalfPath,
    increment_frames,
    MarkovModel,
    MinCut,
    MixtureModel,
    MixtureParams,
    mincut,
    morphological_close,
    morphological_close!,
    neglog_likelihood,
    Path,
    reset!,
    scale,
    segment!

include("Path.jl")
include("Queue.jl")
include("Filter.jl")
include("Components.jl")
include("MinCut.jl")
include("Mixture.jl")
include("Markov.jl")
end # module


text:
# Created by inigo quilez - iq/2013
# License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

# Volumetric clouds. It performs level of detail (LOD) for faster rendering and antialiasing
using MacroTools
using Visualize, GeometryTypes
using Visualize.JLRasterization: Sampler, JLCanvas, JLRasterizer
using Transpiler: mix, fract

Base.one(::Type{Vec3f0}) = Vec3f0(0f0, 0f0, 0f0)

"""
Unroll macro!
There is an unregistered package doing exactly this, called Unroll.jl.
But where'd be the fun in using that
Turns:
```Julia
@unroll for i=1:3
    println(i)
end
```
Into:
```Julia
i = 1
println(i)
i = 2
println(i)
i = 3
println(i)
```
"""
macro unroll(expr)
    idx, range, body = @match expr begin
        for idx_ = range_
            body__
        end => (idx, range, body)
        matchal__ => error("Only for loops like `for i=1:3` are supported")
    end
    block = Expr(:block)
    # Here, eval only works for static expressions.
    # If this would be serious code, we should include some error handling
    for i = eval(range)
        # create a block expression, starting with the assignment of the (now constant)
        # integer to the idx variable, and splicing in the rest of the for body
        # $ splices in the actual value. idx -> a symbol that will turn into a variable
        # i -> the actual index integer
        push!(block.args, :($idx = $i), body...)
    end
    esc(block)
end

function noise(x::Vec3f0)
    p = floor(x)
    f = fract(x)
    f = broadcast(*, f, f)
    f = (3f0 - 2f0*f)
    f = broadcast(*, f, f)
    uv = (p[Vec(1,2)] + Vec2f0(37f0, 17f0) * p[3]) + f[Vec(1,2)]
    rg = Vec2f0(0)
    return -1f0 + 2f0 * mix(rg[1], rg[2], f[3])
end

#noise(Vec3f0(0.242)

@generated function LODaa{N}(p::Vec3f0, time::Float32, ::Val{N})
    result = []
    qf = (2.02f0, 2.03f0, 2.01f0, 2.02f0)
    for i=1:N
        push!(result, :(f += $(2f0^(-i)) * noise(q)))
        i != N && push!(result, :(q = q * $(qf[i])))
    end
    quote
        f = 0f0
        q = p - Vec3f0(0f0, 0f0, 1f0) * time
        $(result...) # splice in unrolled expressions
        res = (1.5f0 - p[2] - 2f0 + 1.75f0) * f
        clamp(res, 0f0, 1f0)
    end
end
#LODaa(rand(Vec3f0), 1f0, Val{5}())

function integrate(
        sum::Vec4f0, dif::Float32, den::Float32,
        bgcol::Vec3f0, t::Float32
    )
    # lighting
    lin = Vec3f0(0.65f0, 0.68f0, 0.7f0) * 1.3f0 + 0.5f0 *
        Vec3f0(0.7f0, 0.5f0, 0.3f0) * dif

    colo = mix(1.15f0 * Vec3f0(1f0, 0.95f0, 0.8f0), Vec3f0(0.65f0), den)
    colo = colo .* lin
    colo = mix(colo, bgcol, 1f0 - exp(-0.003f0 * t*t))
    # front to back blending
    alpha = den * 0.4f0
    colo *= alpha
    return sum .+ Vec4f0(colo[1], colo[2], colo[3], alpha) * (1f0 - sum[4])
end

#
# integrate(
#     Vec4f0(0.2, 0.5, 0.3, 0.8), 2f0, 0.4f0,
#     Vec3f0(0.5, 0.58, 0.9), 1f0
# )

function march{N}(ro, rd, bgcol, sundir, sumt, steps, n::Val{N})
    sum = sumt[1]
    t = sum[2]
    for i = 1:steps
        pos = ro + t * rd
        if (pos[2] < -3f0 || pos[2] > 2f0 || sum[4] > 0.99f0)
            break
        end
        den = LODaa(pos, 1f0, n)
        if den > 0.01f0
            dif = clamp(
                (den - LODaa(pos + 0.3f0 * sundir, 1f0, n)) / 0.6f0,
                0f0, 1f0
            )
            sum = integrate(sum, dif, den, bgcol, t)
        end
        t += max(0.1f0, 0.02f0 * t)
    end
    sum, t
end

# march(
#     Vec3f0(0), Vec3f0(0), Vec3f0(0),
#     normalize(Vec3f0(-1f0, 0f0, -1f0)),
#     Vec4f0(0f0), 0f0, 30, Val{5}()
# )

function raymarch(ro::Vec3f0, rd::Vec3f0, bgcol::Vec3f0, sundir::Vec3f0)
    sumt = (Vec4f0(0f0), 0f0)
    #@unroll for i = 5:-1:2
        sumt = march(ro, rd, bgcol, sundir, sumt, 30, Val{5}())
        sumt = march(ro, rd, bgcol, sundir, sumt, 30, Val{4}())
        sumt = march(ro, rd, bgcol, sundir, sumt, 30, Val{3}())
        sumt = march(ro, rd, bgcol, sundir, sumt, 30, Val{2}())
    #end
    sum = sumt[1]

    return Vec4f0(
        clamp(sum[1], 0f0, 1f0),
        clamp(sum[2], 0f0, 1f0),
        clamp(sum[3], 0f0, 1f0),
        clamp(sum[4], 0f0, 1f0),
    )
end

# raymarch(
#     Vec3f0(0), Vec3f0(0), Vec3f0(0),
#     normalize(Vec3f0(-1f0, 0f0, -1f0)),
# )
#
# @code_warntype(raymarch(
#     Vec3f0(0), Vec3f0(0), Vec3f0(0),
#     normalize(Vec3f0(-1f0, 0f0, -1f0)),
# ))

function render(ro::Vec3f0, rd::Vec3f0, sundir::Vec3f0)
    # background sky
    sun = clamp(dot(sundir, rd), 0f0, 1f0)
    col = Vec3f0(0.6f0, 0.71f0, 0.75f0) - rd[2] *
        0.2f0 * Vec3f0(1f0, 0.5f0, 1f0) + 0.15f0 * 0.5f0

    col += 0.2f0 * Vec3f0(1f0, .6f0, 0.1f0) * (sun ^ 8f0)

    # clouds
    res = raymarch(ro, rd, col, sundir)
    col = col * (1f0 - res[4]) + res[Vec(1,2,3)]

    # sun glare
    col += 0.1f0 * Vec3f0(1f0, 0.4f0, 0.2f0) * (sun ^ 3f0)

    return Vec4f0(col[1], col[2], col[3], 1f0)
end

#render(Vec3f0(0), Vec3f0(0), normalize(Vec3f0(-1f0, 0f0, -1f0)))
function camera(cr = 0f0)
    # camera
    ro = 4f0 * normalize(Vec3f0(sin(3f0 * 1f0), 0.4*1f0, cos(3.0 * 1f0)))
    ta = Vec3f0(0.0, -1.0, 0.0)
    cw = normalize(ta - ro)
    cp = Vec3f0(sin(cr), cos(cr), 0.0)
    cu = normalize(cross(cw, cp))
    cv = normalize(cross(cu, cw))
    return ro, Mat3f0(cu..., cv..., cw...)
end

function frag_main(p, roca)
    ro = roca[1]
    ca = roca[2]
    #p = ptupl[1]
    sundir = normalize(Vec3f0(-1f0, 0f0, -1f0))
    # ray
    rd = ca * normalize(Vec3f0(p[1], p[2], 1.5f0))
    (render(ro, rd, sundir),)
end

function vert_main(p, roca)
    Vec4f0(p[1], p[2], 0, 1), p
end

noise_sampler = Sampler(rand(Vec3f0, 512, 512))
typeof(noise_sampler)
ro, ca = camera()
canvas = JLCanvas(Area => IRect(0, 0, 512, 512))
mesh = uvmesh(SimpleRectangle(-1f0, -1f0, 2f0, 2f0))
draw = JLRasterizer(
    mesh, (ro, ca),
    vert_main, frag_main,
)

draw(canvas, mesh, (ro, ca))
using FileIO
fb = canvas[FieldTraits.Fields(canvas)[3]][1]
save("/home/s/Desktop/test.png", fb)

canvas = GLFWWindow()

rect = SimpleRectangle(-1f0, -1f0, 2f0, 2f0)
mesh = Base.view(decompose(Point2f0, rect), decompose(GLTriangle, rect))
ro, ca = camera()
vbo = VertexArray(mesh)
uniforms = UniformBuffer((ro, ca))
draw = GLRasterizer(
    mesh, (uniforms,),
    vert_main, frag_main,
)


text:
function fibonacci(n)
    answer = pow(n-1)
    return answer
end

function pow(b)
    k = b
    t = [1 1
        1 1]
    p = [1 1
         1 0]
    while k>0
        if k%2==0
            k ÷= 2
            p *= p
        else
            k -= 1
            t *= p
        end
    end
    return t[2,2]
end


text:
function minTest1()
  Rect3x2::Matrix{Float64} = [0.1 0.2; 0.3 0.4;0.5 0.6]
  Rect2x3::Matrix{Float64} = zeros(2,3)
  Rect2x3' .= Rect3x2
  println("Still Zeros: ", Rect2x3)
  Rect2x3 .= Rect3x2'
  println("Seems to Work: ", Rect2x3)
end

function foo(A, x, b=Vector{Float64}(length(x)))
  b .= A * x
  return b
end

println(foo([1 2; 3 4], [1,2]))
----------------
#this tests for a notice period negative bias to performance
#this version
function NNP_1D(HFLong::DataFrame; dependentVar::Symbol = :performance1D,
    titleText::String =  "Performance Near Withdawels: Impact of Notice Period")::String
    thresholdLS::RedemptionThresholds = RedemptionThresholds(SIG_REDEMPTION_THRESHOLDS)

    modelsNNP_1D::Vector{CTLM} = Vector{CTLM}()

    XNNP_1DSpecs::Vector{CTExpr} = Vector{CTExpr}()
    XNNP_1DNames::Vector{Vector{Symbol}}  = Vector{Vector{Symbol}}()

    push!(XNNP_1DSpecs, parse("noticePeriod"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod])

    push!(XNNP_1DSpecs, parse("noticePeriod + main_strategy"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod])

    push!(XNNP_1DSpecs, parse("noticePeriod + fRedemptionDate + main_strategy"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod])

    push!(XNNP_1DSpecs, parse("noticePeriod + fDate + main_strategy"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod])

    push!(XNNP_1DSpecs, parse("noticePeriod&negLFlows"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod_negLFlows])

    push!(XNNP_1DSpecs, parse("noticePeriod&negLFlows + main_strategy"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod_negLFlows])

    push!(XNNP_1DSpecs, parse("noticePeriod&negLFlows + fRedemptionDate + main_strategy"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod_negLFlows])

    push!(XNNP_1DSpecs, parse("noticePeriod&negLFlows + fDate + main_strategy"))
    push!(XNNP_1DNames, [:intercept; :noticePeriod_negLFlows])

    #push the specifications (use a common name to make a clean table)
    #=for s::Symbol ∈ thresholdLS.negLFlows

        push!(XNNP_1DSpecs, parse("noticePeriod&$s"))
        push!(XNNP_1DNames, [:intercept; :noticePeriod])

        push!(XNNP_1DSpecs, parse("noticePeriod&$s + fDate + main_strategy"))
        push!(XNNP_1DNames, [:intercept; :noticePeriod])

    end=#

    #run the model
    for iX ∈ 1:length(XNNP_1DSpecs)
        push!(modelsNNP_1D, CTLM(HFLong, XNNP_1DSpecs[iX],
            dependentVar, XNames=XNNP_1DNames[iX], YName = :performance))
    end

    ######IO Code for performance bias analysis

    descRowNamesNNP_1D::Vector{String} = ["Date F.E.",  "Redemp. Date F.E.","Strategy F.E.",#="Outflow \\% >", =# "N ('000s)"]

    #need to print the descriptive rows
    descContentNNP_1D::Vector{Vector{String}} =
        [Vector{String}(length(modelsNNP_1D)) for i∈ 1:length(descRowNamesNNP_1D)]

    #Determine the threshold levels (searches for the threshold level in each spec)
    thresholdVec::Vector{Int} = ((i::Int)->maximum(
        ((j::Int)->contains(XNNP_1DSpecs[i],"$(thresholdLS.negLFlows[j])")?
            Int(100.0*collect(SIG_REDEMPTION_THRESHOLDS)[j]):0).(1:length(thresholdLS.negLFlows))
        )).(1:length(XNNP_1DSpecs))

    for i ∈ 1:length(XNNP_1DSpecs)
        descContentNNP_1D[1][i] = "$(contains(XNNP_1DSpecs[i], "fDate")?"X":"")"
        descContentNNP_1D[2][i] = "$(contains(XNNP_1DSpecs[i], "fRedemptionDate")?"X":"")"
        descContentNNP_1D[3][i] = "$(contains(XNNP_1DSpecs[i], "main_strategy")?"X":"")"
        #descContentNNP_1D[3][i] = "$(contains(XNNP_1DSpecs[i], "negLFlows")?"X":"")"
        #descContentNNP_1D[4][i] = "$(thresholdVec[i])\\%"
        descContentNNP_1D[end][i] = "$(round(modelsNNP_1D[i].N/1000.0,1))"
    end

    #rows and variable names for the regression
    rowsNNP_1D::Vector{Symbol} = [:intercept; :noticePeriod; :noticePeriod_negLFlows]

    TableTextNNP_1D::String = texTable(modelsNNP_1D, getHomoskedΣ!, rowsNNP_1D,
        titleCaption = "Performance Flow-Bias Specifications",
        colNames = [["($i)" for i∈1:length(modelsNNP_1D)]],
        contentRowNames = ["intercept", "Notice Period", "Notice X Outflow"],
        descRowNames = descRowNamesNNP_1D,
        descContent = descContentNNP_1D,
        decimalDigits = 2,
        columnSepPt = -20,
        scaling = fill(100.0,length(rowsNNP_1D))::Vector{Float64},
        clearMem = USE_AGGRESSIVE_GC,
        caption = """NNP_1D specification compares performance prior to notice,
            with performance after notice. In some specifications, the notice
            period dummy is interacted with the amount of the outflow. Threshold
            refers to the minimum outflow included, thus some versions include
            only significant outflows. Data is from alive and dead hedge funds
            2000-2017. Interpret a value of 1.0 for the non-interacted
            specification as the discontinuity increasing performance by 1\\%.
            For the interacted specification, a 1.0 implies that outflows of
            10\\% corresponding to a 0.1\\% drop in monthly performance.""")

    return TableTextNNP_1D
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using JuMP
using Ipopt
using CPLEX
include("io.jl")

function create_and_solve_det_milp(data::Dict{AbstractString,Any}, ref::Dict{Symbol,Any}, attacks::Vector{Attack}; k = 2)

    ϵ = 1e-5
    m = Model()
    # setsolver(m, GurobiSolver(LogToConsole=0))
    # setsolver(m, GurobiSolver())
    setsolver(m, CplexSolver(CPX_PARAM_PREIND=0))

    notgen_bus_indexes = Int[] # includes plain buses, buses with loads and condensors
    gen_bus_indexes = Int[] # includes only generators and not condensors

    w_vertex = Dict{Int,Float64}()
    for i in keys(ref[:bus])
        buses = ref[:bus]
        if buses[i]["bus_type"] == 2 || buses[i]["bus_type"] == 3
            total_pmin = sum([abs(ref[:gen][g]["pmin"]) for g in ref[:bus_gens][i]])
            total_pmax = sum([abs(ref[:gen][g]["pmax"]) for g in ref[:bus_gens][i]])
            is_condensor = (total_pmin == 0 && total_pmax == 0) ? true : false
            if is_condensor
                w_vertex[i] = abs(buses[i]["pd"])
                push!(notgen_bus_indexes, i)
                continue
            else
                w_vertex[i] = 0.0
                push!(gen_bus_indexes, i)
                continue
            end
        else
            w_vertex[i] = abs(buses[i]["pd"])
            push!(notgen_bus_indexes, i)
        end
    end

    @variable(m, x[i in keys(ref[:branch])], Bin) # x[e] = 1 means the edge is removed
    @variable(m, y[v in keys(ref[:bus])], Bin) # y[v] = 1 means the vertex is isolated
    @variable(m, dummy_edge[v in gen_bus_indexes] == 0) # dummy edges are added between a source and a gen_bus

    w_edge = Dict(i => 1 for i in keys(ref[:branch]))
    @constraint(m, gen_off[v in gen_bus_indexes], y[v] == 0)
    @constraint(m, sum(w_edge[i]*x[i] for i in keys(ref[:branch])) == k)
    @objective(m, Max, sum(w_vertex[i]*y[i] for i in keys(ref[:bus])))

    index_map = IndexMap()
    bus_indexes = collect(keys(ref[:bus]))
    for i in 1:length(bus_indexes)
        index_map.bus_to_vertex[bus_indexes[i]] = i
        index_map.vertex_to_bus[i] = bus_indexes[i]
    end
    index_map.bus_to_vertex[0] = length(keys(ref[:bus])) + 1
    index_map.vertex_to_bus[length(keys(ref[:bus]))+1] = 0

    function lazycuts(cb)
        x_vals = Dict(i => getvalue(x[i]) for i in keys(ref[:branch]))
        y_vals = Dict(v => getvalue(y[v]) for v in notgen_bus_indexes)
        y_vals = filter((v, val) -> val >= (1-ϵ), y_vals)
        g = simple_graph(length(keys(index_map.vertex_to_bus)), is_directed=false)
        dists = Float64[]
        # create support graph edges
        for (key, value) in x_vals
            f_bus = ref[:branch][key]["f_bus"]
            t_bus = ref[:branch][key]["t_bus"]
            i = index_map.bus_to_vertex[f_bus]
            j = index_map.bus_to_vertex[t_bus]
            add_edge!(g, i, j)
            if !haskey(index_map.parallel_edges, (i, j))
                index_map.parallel_edges[(i, j)] = Dict{AbstractString, Any}("count" => 1, "index" => [key])
                index_map.parallel_edges[(j, i)] = Dict{AbstractString, Any}("count" => 1, "index" => [key])
            else
                index_map.parallel_edges[(i, j)]["count"] += 1
                push!(index_map.parallel_edges[(i, j)]["index"], key)
                index_map.parallel_edges[(j, i)]["count"] += 1
                push!(index_map.parallel_edges[(j, i)]["index"], key)
            end
            push!(dists, value)
        end
        # create support graph dummy edges
        for gen in gen_bus_indexes
            f_bus = 0
            t_bus = gen
            i = index_map.bus_to_vertex[f_bus]
            j = index_map.bus_to_vertex[t_bus]
            add_edge!(g, i, j)
            if !haskey(index_map.parallel_edges, (i, j))
                index_map.parallel_edges[(i, j)] = Dict{AbstractString, Any}("count" => 1, "index" => [0])
                index_map.parallel_edges[(j, i)] = Dict{AbstractString, Any}("count" => 1, "index" => [0])
            end
            push!(dists, 0.0)
        end
        r = dijkstra_shortest_paths(g, dists, index_map.bus_to_vertex[0])
        if length(keys(y_vals)) != 0
            paths = enumerate_indices(r.parent_indices, [index_map.bus_to_vertex[i] for i in keys(y_vals)])
            for path in paths
                target_bus = index_map.vertex_to_bus[path[length(path)]]
                @assert path[1] == length(keys(ref[:bus]))+1
                @assert target_bus in notgen_bus_indexes
                @assert target_bus in keys(y_vals)
                path_indexes = Int[]
                for i in 2:length(path)-1
                    path_index = 0.0
                    index = index_map.parallel_edges[(path[i], path[i+1])]
                    if index["count"] > 1
                        path_index = (x_vals[index["index"][1]] > 1-ϵ) ? index["index"][2] : index["index"][1]
                    else
                        path_index = index["index"][1]
                    end
                    push!(path_indexes, path_index)
                end

                if sum([x_vals[i] for i in path_indexes]) < y_vals[target_bus]
                    @lazyconstraint(cb, sum(x[i] for i in path_indexes) >= y[target_bus])
                end
            end
        end
    end

    addlazycallback(m, lazycuts)
    status = solve(m)
    x_vals = Dict(i => getvalue(x[i]) for i in keys(ref[:branch]))
    x_vals = filter((i, val) -> val > (1-ϵ), x_vals)
    y_vals = Dict(i => getvalue(y[i]) for i in notgen_bus_indexes)
    y_vals = filter((i, val) -> val > (1-ϵ), y_vals)
    attack = Attack()
    for bus_id in collect(keys(y_vals))
        push!(attack.bus_ids, bus_id)
    end
    for branch_id in collect(keys(x_vals))
        push!(attack.branch_ids, branch_id)
    end
    attack.total_load = getobjectivevalue(m)
    return attack
end

function create_and_solve_stoch_milp(data::Dict{AbstractString,Any}, ref::Dict{Symbol,Any}, attacks::Vector{Attack}, log_p; k = 2)

    ϵ = 1e-5
    m = Model()
    setsolver(m, CplexSolver(CPX_PARAM_PREIND=0))

    notgen_bus_indexes = Int[] # includes plain buses, buses with loads and condensors
    gen_bus_indexes = Int[] # includes only generators and not condensors

    w_vertex = Dict{Int,Float64}()
    for i in keys(ref[:bus])
        buses = ref[:bus]
        if buses[i]["bus_type"] == 2 || buses[i]["bus_type"] == 3
            total_pmin = sum([abs(ref[:gen][g]["pmin"]) for g in ref[:bus_gens][i]])
            total_pmax = sum([abs(ref[:gen][g]["pmax"]) for g in ref[:bus_gens][i]])
            is_condensor = (total_pmin == 0 && total_pmax == 0) ? true : false
            if is_condensor
                w_vertex[i] = abs(buses[i]["pd"])
                push!(notgen_bus_indexes, i)
                continue
            else
                w_vertex[i] = 0.0
                push!(gen_bus_indexes, i)
                continue
            end
        else
            w_vertex[i] = abs(buses[i]["pd"])
            push!(notgen_bus_indexes, i)
        end
    end

    @variable(m, x[i in keys(ref[:branch])], Bin) # x[e] = 1 means the edge is removed
    @variable(m, y[v in keys(ref[:bus])], Bin) # y[v] = 1 means the vertex is isolated
    @variable(m, dummy_edge[v in gen_bus_indexes] == 0) # dummy edges are added between a source and a gen_bus
    @variable(m, p <= 0)
    @variable(m, z <= 1e6)
    @variable(m, 1e-6 <= eta <= 1e6)

    w_edge = Dict(i => 1 for i in keys(ref[:branch]))
    @constraint(m, gen_off[v in gen_bus_indexes], y[v] == 0)
    @constraint(m, sum(w_edge[i]*x[i] for i in keys(ref[:branch])) == k)
    @constraint(m, p == sum(x[i]*log_p[i] for i in keys(ref[:branch])))
    @constraint(m, eta == sum(w_vertex[i]*y[i] for i in keys(ref[:bus])))
    @objective(m, Max, z)

    index_map = IndexMap()
    bus_indexes = collect(keys(ref[:bus]))
    for i in 1:length(bus_indexes)
        index_map.bus_to_vertex[bus_indexes[i]] = i
        index_map.vertex_to_bus[i] = bus_indexes[i]
    end
    index_map.bus_to_vertex[0] = length(bus_indexes) + 1
    index_map.vertex_to_bus[length(bus_indexes)+1] = 0

    function lazycuts(cb)
        eta_val = getvalue(eta)
        p_val = getvalue(p)
        if eta_val < ϵ
            eta_val = 0.01
        end
        if eta_val > 0.0
            inv_eta_val = 1/eta_val
            @lazyconstraint(cb, z <= p + log(eta_val) + inv_eta_val*(eta - eta_val))
        end

        x_vals = Dict(i => getvalue(x[i]) for i in keys(ref[:branch]))
        y_vals = Dict(v => getvalue(y[v]) for v in notgen_bus_indexes)
        y_vals = filter((v, val) -> val >= (1-ϵ), y_vals)
        g = simple_graph(length(keys(index_map.vertex_to_bus)), is_directed=false)
        dists = Float64[]
        # create support graph edges
        for (key, value) in x_vals
            f_bus = ref[:branch][key]["f_bus"]
            t_bus = ref[:branch][key]["t_bus"]
            i = index_map.bus_to_vertex[f_bus]
            j = index_map.bus_to_vertex[t_bus]
            add_edge!(g, i, j)
            if !haskey(index_map.parallel_edges, (i, j))
                index_map.parallel_edges[(i, j)] = Dict{AbstractString, Any}("count" => 1, "index" => [key])
                index_map.parallel_edges[(j, i)] = Dict{AbstractString, Any}("count" => 1, "index" => [key])
            else
                index_map.parallel_edges[(i, j)]["count"] += 1
                push!(index_map.parallel_edges[(i, j)]["index"], key)
                index_map.parallel_edges[(j, i)]["count"] += 1
                push!(index_map.parallel_edges[(j, i)]["index"], key)
            end
            push!(dists, value)
        end
        # create support graph dummy edges
        for gen in gen_bus_indexes
            f_bus = 0
            t_bus = gen
            i = index_map.bus_to_vertex[f_bus]
            j = index_map.bus_to_vertex[t_bus]
            add_edge!(g, i, j)
            if !haskey(index_map.parallel_edges, (i, j))
                index_map.parallel_edges[(i, j)] = Dict{AbstractString, Any}("count" => 1, "index" => [0])
                index_map.parallel_edges[(j, i)] = Dict{AbstractString, Any}("count" => 1, "index" => [0])
            end
            push!(dists, 0.0)
        end
        r = dijkstra_shortest_paths(g, dists, index_map.bus_to_vertex[0])
        if length(keys(y_vals)) != 0
            paths = enumerate_indices(r.parent_indices, [index_map.bus_to_vertex[i] for i in keys(y_vals)])
            for path in paths
                target_bus = index_map.vertex_to_bus[path[length(path)]]
                @assert path[1] == length(keys(ref[:bus]))+1
                @assert target_bus in notgen_bus_indexes
                @assert target_bus in keys(y_vals)
                path_indexes = Int[]
                for i in 2:length(path)-1
                    path_index = 0.0
                    index = index_map.parallel_edges[(path[i], path[i+1])]
                    if index["count"] > 1
                        path_index = (x_vals[index["index"][1]] > 1-ϵ) ? index["index"][2] : index["index"][1]
                    else
                        path_index = index["index"][1]
                    end
                    push!(path_indexes, path_index)
                end

                if sum([x_vals[i] for i in path_indexes]) < y_vals[target_bus]
                    @lazyconstraint(cb, sum(x[i] for i in path_indexes) >= y[target_bus])
                end
            end
        end
    end

    addlazycallback(m, lazycuts)
    status = solve(m)
    x_vals = Dict(i => getvalue(x[i]) for i in keys(ref[:branch]))
    x_vals = filter((i, val) -> val > (1-ϵ), x_vals)
    y_vals = Dict(i => getvalue(y[i]) for i in notgen_bus_indexes)
    y_vals = filter((i, val) -> val > (1-ϵ), y_vals)
    p_val = getvalue(p)
    actual_prob = exp(p_val)
    log_eta = getobjectivevalue(m) - p_val
    eta_val = exp(log_eta)
    exp_load_shed = actual_prob*eta_val
    attack = Attack()
    for bus_id in collect(keys(y_vals))
        push!(attack.bus_ids, bus_id)
    end
    for branch_id in collect(keys(x_vals))
        push!(attack.branch_ids, branch_id)
    end
    attack.total_load = exp_load_shed
    return attack
end

type IndexMap
    bus_to_vertex::Dict{Int,Int}
    vertex_to_bus::Dict{Int,Int}
    parallel_edges::Dict{Any,Any}
    function IndexMap()
        imap = new()
        imap.bus_to_vertex = Dict{Int,Int}()
        imap.vertex_to_bus = Dict{Int,Int}()
        imap.parallel_edges = Dict{Any,Any}()
        return imap
    end
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
"""
    Base.iszero(t::Tuple)

Extends `iszero` to `Tuple` types. Return `true` iff all elements of `t` are equal to zero.
"""
Base.iszero(t::Tuple) = all(iszero, t)

"""
    Base.occursin(re, char::Char)
Extends `occursin` to test if a regular expression `re` is a match with a `Char`.
"""
Base.occursin(re, char::Char) = occursin(re, string(char))

"""
    name_to_string(name::BibInternal.Name)

Convert a name in an `Entry` to a string.
"""
function name_to_string(name)
    str = "$(name.particle)"
    if str != "" != name.last
        str *= " "
    end
    str *= name.last
    str *= name.junior == "" ? "" : ", $(name.junior)"
    if name.first != ""
        str *= ", $(name.first)"
    end
    if name.middle != ""
        str *= " $(name.middle)"
    end
    return str
end


text:
module PlotShapefiles

##############################################################################
##
## Dependencies and Reexports
##
##############################################################################

using Reexport
@reexport using Shapefile
@reexport using Compose
@reexport using Measures
@reexport using Colors
@reexport using Images
using Cairo: CairoARGBSurface
using FixedPointNumbers: Normed

import Base.identity

##############################################################################
##
## Exported methods and types
##
##############################################################################

export open_shapefile,
       lonlat_to_webmercator, identity,
       plotshape,
       choropleth,
       google_overlay, #CairoARGBSurface,
       rootpath, testdatapath

##############################################################################
##
## Load source files
##
##############################################################################

include("coordinate_conversions.jl")
include("shp_point.jl")
include("shp_multipoint.jl")
include("shp_polygon.jl")
include("shp_polyline.jl")
include("utils.jl")
include("plotshape.jl")
include("choropleth.jl")
include("google_overlay.jl")

end # module


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

# Wavelength range [nm]
λ = 200:2100
# Angle of incidence [degrees]
θ = [60.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.fused_silica_uv(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:tauclorentz),
    LayerTMMO(emergent),
]

# Measured absolute transmittance
Ψexp, Δexp, = SpectraDB.hafnia_ellips(beam.λ)
plot(beam.λ, Ψexp)
plot!(beam.λ, Δexp)
gui()

# Spectrum type: ellipsometry with the input experimental spectra
spectype = Ellipsometry([Ψexp Δexp])

## Optimization using BBO

# One oscillator
seed = [vcat(100.0, # thickness
             [3.0, 5.5,
              545, 6.2, 13.2], # osc 1
        ),
]
lb = 0.01.*seed
ub = 2.0.*seed
lb[1][1] = 80
ub[1][1] = 200
lb[1][3] = 4.0
ub[1][3] = 6.0

solOptim = fit_tmm_optics(
    spectype, seed, beam, layers;
    alg=:BBO,
    SearchRange=Utils.unfold_bnd(lb,ub),
)

plot(FitSpectrumEllip(),
    solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit,
    xaxis=("Wavelength [nm]"), yaxis=(),
)
gui()


text:
# This file is a part of MJDSigGen, licensed under the MIT License (MIT).

function read_config!(setup::SigGenSetup, config_filename::AbstractString)
    result = ccall(
        @sgsym(:read_config), Cint,
        (Cstring, Ptr{SigGenSetup}),
        config_filename, Ref(setup)
    )

    result != 0 && error("read_config failed.")
    return setup
end

read_config(config_filename::AbstractString) =
    read_config!(SigGenSetup(), config_filename)

function signal_calc_init!(setup::SigGenSetup, config_filename::AbstractString)
    result = ccall(
        @sgsym(:signal_calc_init), Cint,
        (Cstring, Ptr{SigGenSetup}),
        config_filename, Ref(setup)
    )
    result != 0 && error("signal_calc_init failed.")
    finalizer(signal_calc_finalize!, setup)
    return setup
end

signal_calc_init(config_filename::AbstractString) =
    signal_calc_init!(SigGenSetup(), config_filename)

function signal_calc_finalize!(setup::SigGenSetup)
    result = ccall(
        @sgsym(:signal_calc_finalize), Cint,
        (Ptr{SigGenSetup},),
        Ref(setup)
    )
    result != 0 && error("signal_calc_finalize failed.")
    return setup
end

function get_signal!(signal::DenseArray{Float32, 1}, setup::SigGenSetup, location::NTuple{3, Real})
    sl = length(signal)
    nsteps = setup.ntsteps_out
    sl != nsteps && throw(
        ArgumentError("signal length ($sl) must equal number of time steps ($nsteps)"))

    pt = CartPoint{Cfloat}(location[1], location[2], location[3])

    result = ccall(
        @sgsym(:get_signal), Cint,
        (CartPoint{Cfloat}, Ptr{Float32}, Ptr{SigGenSetup}),
        pt, signal, Ref(setup)
    )
    result < 0 && throw(ArgumentError("point not in crystal or has no field: $pt"))

    return signal
end

get_signal!(setup::SigGenSetup, location::NTuple{3, Real}) =
    get_signal!(zeros(Float32, setup.ntsteps_out), setup, location)

function outside_detector(setup, location::NTuple{3, Real})
    pt = CartPoint{Cfloat}(location[1], location[2], location[3])

    r = ccall(
        @sgsym(:outside_detector), Cint,
        (CartPoint{Cfloat}, Ptr{SigGenSetup}),
        pt, Ref(setup)
    )

    return r != 0
end

"""
    nearest_field_grid_index(setup::SigGenSetup, location::NTuple{3, Real})

Returns:
* (:outside, i, j), if outside crystal or too far from a valid grid point
* (:interpol, i, j, if interpolation is okay
* (:extrapol, i, j), if we can find a point but extrapolation is needed
"""
function nearest_field_grid_index(setup::SigGenSetup, location::NTuple{3, Real})
    r, ϕ, z = cart2cyl(location ...)
    cyl_pos_val = CylPoint{Cfloat}(r, ϕ, z)
    cyl_idx_ref = Ref(CylPoint{Cint}(0, 0, 0))

    retcode = ccall(
        @sgsym(:nearest_field_grid_index), Cint,
        (CylPoint{Cfloat}, Ptr{CylPoint{Cint}}, Ptr{SigGenSetup}, ),
        cyl_pos_val, cyl_idx_ref, Ref(setup)
    )

    cyl_idx = cyl_idx_ref.x

    i,j = Int(cyl_idx.z + 1), Int(cyl_idx.r + 1)

    if retcode < 0
        (:outside, 0, 0)
    elseif retcode == 0
        (:interpol, i, j)
    elseif retcode == 1
        (:extrapol, i, j)
    else
        throw(ArgumentError(
            "Don't know how to interpret return code $retcode of
            nearest_field_grid_index unknown"
        ))
    end
end

const fieldgen_exe = joinpath(dirname(@__FILE__), "..", "deps", "usr", "bin", "mjd_fieldgen")

function fieldgen(config_filename::AbstractString)
    setup = read_config(config_filename)

    mkpath(dirname(joinpath(dirname(config_filename), setup.field_name)))
    mkpath(dirname(joinpath(dirname(config_filename), setup.wp_name)))

    run(`$fieldgen_exe -c $config_filename`)
end

function get_drift_velocity(setup::SigGenSetup, location::NTuple{3, Real}, t::Symbol)
	if t==:e
		q = -1;
	elseif t==:h
		q = +1;
	else
		error("Charge carrier type must be :e or :h")
	end

    pt	= CartPoint{Cfloat}(location[1], location[2], location[3])
	vel	= Ref(CartPoint{Cfloat}(0, 0, 0))

    ccall(
        @sgsym(:drift_velocity), Cint,
        (CartPoint{Cfloat}, Float32, Ptr{CartPoint{Cfloat}}, Ptr{SigGenSetup}),
        pt, q, vel, Ref(setup)
    ) < 0 && error("Point not in crystal or has no field: $pt")

	#asd = unsafe_load(vel);
    return vel.x;
end

function get_drift_velocity_w_Eadd(setup::SigGenSetup, location::NTuple{3, Real}, t::Symbol, Eadd_cart::NTuple{3, Real})
	if t==:e
		q = -1;
	elseif t==:h
		q = +1;
	else
		error("Charge carrier type must be :e or :h")
	end

    pt	= CartPoint{Cfloat}(location[1], location[2], location[3])
	vel	= Ref(CartPoint{Cfloat}(0, 0, 0))
	Er, Eϕ, Ez = cart2cyl(Eadd_cart ...);
	if(Eadd_cart[1]<0)
		Er =-Er;
		Eϕ = 0;
	end
	Eadd = CylPoint{Cfloat}(Er, Eϕ, Ez);

    ccall(
        @sgsym(:drift_velocity_w_Eadd), Cint,
        (CartPoint{Cfloat}, Float32, Ptr{CartPoint{Cfloat}},
        Ptr{SigGenSetup}, CylPoint{Cfloat}),
        pt, q, vel, Ref(setup), Eadd
    ) < 0 && error("Point not in crystal or has no field: $pt")

	#asd = unsafe_load(vel);
    return vel.x;
end

function get_drift_velocity_from_Efield(setup::SigGenSetup, location::NTuple{3, Real}, t::Symbol, Efield_cart::NTuple{3, Real})
	if t==:e
		q = -1;
	elseif t==:h
		q = +1;
	else
		error("Charge carrier type must be :e or :h")
	end

    pt	= CartPoint{Cfloat}(location[1], location[2], location[3])
	vel	= Ref(CartPoint{Cfloat}(0, 0, 0))
	Er, Eϕ, Ez = cart2cyl(Efield_cart...);
	if(Efield_cart[1]<0)
		Er =-Er;
		Eϕ = 0;
	end
	Efield = CylPoint{Cfloat}(Er, Eϕ, Ez);

    ccall(
        @sgsym(:drift_velocity_from_Efield), Cint,
        (CartPoint{Cfloat}, Float32, Ptr{CartPoint{Cfloat}}, Ptr{SigGenSetup}, CylPoint{Cfloat}),
        pt, q, vel, Ref(setup), Efield
    ) < 0 && error("Point not in crystal or has no field: $pt")

    return Tuple(vel.x);
end


text:
##############################################################################
##
## Weight
##
##############################################################################
function get_weights(df::AbstractDataFrame, esample::AbstractVector, weights::Symbol)
    # there are no NA in it. DataVector to Vector
    out = convert(Vector{Float64}, df[esample, weights])
    map!(sqrt, out, out)
    return out
end
get_weights(df::AbstractDataFrame, esample::AbstractVector, ::Nothing) = Ones{Float64}(sum(esample))

#  remove observations with missing or negative weights
function isnaorneg(a::Vector{T}) where {T}
    out = BitArray(undef, length(a))
    @inbounds @simd for i in 1:length(a)
        out[i] = !ismissing(a[i]) & (a[i] > zero(T))
    end
    return out
end