text:
function optimize!(structure::DeterministicEquivalent, optimizer::MOI.AbstractOptimizer, x₀::AbstractVector)
    # Sanity check
    backend(structure.model) === optimizer || error("Stochastic program optimizer has not been connected to the deterministically equivalent problem.")
    # Crash if supported
    if MOI.supports(optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)
        for (i, idx) in enumerate(all_decisions(structure))
            MOI.set(optimizer, MOI.VariablePrimalStart(), idx, x₀[i])
        end
    end
    # Run standard MOI optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function supports_structure(::MOI.AbstractOptimizer, ::DeterministicEquivalent)
    return true
end

function check_loadable(::MOI.AbstractOptimizer, ::DeterministicEquivalent)
    return nothing
end

function set_master_optimizer!(structure::DeterministicEquivalent, optimizer)
    # Ensure decision bridges are added
    for bridge_type in structure.model.bridge_types
        JuMP._moi_add_bridge(structure.model.moi_backend, bridge_type)
    end
    return nothing
end

function set_master_optimizer_attribute!(structure::DeterministicEquivalent, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(backend(structure.model), attr, value)
    return nothing
end

function set_subproblem_optimizer!(structure::DeterministicEquivalent, optimizer)
    return nothing
end

function set_subproblem_optimizer_attribute!(::DeterministicEquivalent, ::MOI.AbstractOptimizerAttribute, value)
    return nothing
end

function optimize!(structure::DeterministicEquivalent, optimizer::AbstractStructuredOptimizer)
    # Sanity check
    supports_structure(optimizer, typeof(structure)) || throw(UnsupportedStructure{typeof(optimizer), typeof(structure)}())
    # Ensure that no decisions are fixed in the first stage
    untake_decisions!(structure.model, structure.decision_variables[1])
    # Load structure
    load_structure!(optimizer, structure)
    # Run structure-exploiting optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end


text:
using Test, DelayEmbeddings

println("\nTesting utils...")

@testset "utils" begin

    A = orthonormal(50, 50)
    @test size(A) == (50, 50)
    @test A isa Matrix

    B = orthonormal(3, 3)
    @test B isa SMatrix
end

ext: jl


text:
# solver file to be loaded into eFEMpart

################################
###### General Framework #######
################################

import LinearAlgebra: lu

# General solver for scalar problems
function solve(prob::Problem,mesh::S,param::T) where
                {S<:AbstractMesh,T<:AbstractParameter}
  if !(prob.OperatorType in [:Poisson2D,:AdvDiff2D,:Stokes2D,:Darcy2D,:PoissonAS,:StokesAS,:AdvDiffAS,:Brinkman2D,:BrinkmanMP2D])
    error("operator type not recognized")
  end

  # generate linear operator matrix
  LinOp = GenerateSystem(mesh,prob,param)

  # Apply boundary conditions
  ApplyBC!(LinOp,mesh,prob,param,prob.OperatorType)

  # apply solver
  sol = LinearSolve(mesh,LinOp,prob,param,prob.OperatorType)

  return sol
end
# using new boundary condition application during matrix assembly
function solve_darcyvar_2b(prob,mesh,param)
    LinOp = LinearOperator(mesh,prob)

    # generate linear operator matrix & apply BC!
    start = time()
    LinOp.Op,LinOp.rhs = Darcy2DMatrix_2b(mesh,prob.bcval[:forcing],param,prob)
    #println("linear system generation & bc apply (together):")
    println("linear system generation & bc apply:")
    println("  ",time() - start," seconds")

    #start = time()
    #ApplyBC!(LinOp,mesh,prob,param,prob.OperatorType)
    #println("linear system generation & bc apply (together):")
    #println("  ",time() - start," seconds")

    # apply solver
    start = time()
    #stemp = @elapsed(sol = LinearSolve(mesh,LinOp,prob,param,prob.OperatorType))
    sol = LinearSolve(mesh,LinOp,prob,param,prob.OperatorType)
    println("solver time:")
    println("  ",time() - start," seconds")

    return sol
end

# solve function without parameter included (for :Poisson2D)
function solve(prob::Problem,mesh::S) where S<:AbstractMesh
  param = PoissonParam(1.0);  return solve(prob,mesh,param)
end

# General solver for scalar problems
function solve!(sol::R,prob::Problem,mesh::S,param::T) where
                {R<:AbstractSolution,S<:AbstractMesh,T<:AbstractParameter}
  # generate linear operator matrix
  LinOp = GenerateSystem(mesh,prob,param)

  # Apply boundary conditions
  ApplyBC!(LinOp,mesh,prob,param,prob.OperatorType)

  # apply solver
  LinearSolve!(sol,mesh,LinOp,prob,param,prob.OperatorType)

  return nothing
end

##############################
###### Generate Matrix #######
##############################

function GenerateSystem(mesh::R,prob::S,param::T) where
          {R<:AbstractMesh,S<:AbstractProblem,T<:AbstractParameter}
  LinOp = LinearOperator(mesh,prob)

  # POISSON ✓
  if prob.OperatorType == :Poisson2D
    LinOp.Op,LinOp.rhs = Laplace2DMatrix(mesh,prob.bcval[:forcing],param)
    LinOp.Op *= param.κ

  # DARCY -- NOT TESTED YET
  elseif prob.OperatorType == :Darcy2D
    LinOp.Op,LinOp.rhs = Darcy2DMatrix(mesh,prob.bcval[:forcing],param)

  # ADVECTION-DIFFUSION 2D ✓
  elseif prob.OperatorType == :AdvDiff2D
    LinOp.Op,LinOp.rhs = AdvDiff2DMatrix(mesh,prob.bcval[:forcing],param)

  # Stokes 2D ✓
  elseif prob.OperatorType == :Stokes2D
    LinOp.Op,LinOp.rhs = Stokes2DMatrix(mesh,prob,param)

  # BRINKMAN 2D  -- NOT TESTED YET
  elseif prob.OperatorType == :Brinkman2D
    LinOp.Op,LinOp.rhs = Brinkman2DMatrix(mesh,prob,param)

  # BRINKMAN MULTIPHASE 2D -- NOT TESTED YET
  elseif prob.OperatorType == :BrinkmanMP2D
    LinOp.Op,LinOp.rhs = BrinkmanMP2DMatrix(mesh,prob,param)

  # ADVECTION-DIFFUSION AXISYMMETRIC
  elseif prob.OperatorType == :AdvDiffAS
    LinOp.Op,LinOp.rhs = AdvDiffASMatrix(mesh,prob.bcval[:forcing],param)

  # STOKES AXISYMMETRIC
  elseif prob.OperatorType == :StokesAS
    LinOp.Op,LinOp.rhs = StokesASMatrix(mesh,prob,param)

  elseif prob.OperatorType == :PoissonAS
    LinOp.Op,LinOp.rhs = LaplaceASMatrix(mesh,prob.bcval[:forcing],param)
  end

  return LinOp
end

########################################
###### Apply Boundary Conditions #######
########################################

# solve boundary condition
function ApplyBC!(LinOp::LinearOperator,mesh::S,prob::Problem,param::T,
                  OpType::Symbol) where {S<:AbstractMesh,T<:AbstractParameter}
  # POISSON 2D or DARCY 2D
  if OpType == :Poisson2D || OpType == :Darcy2D
    if :dNodes in bc(prob)
      scalarDirichlet!(prob.bcid[:dNodes],prob.bcval[:dBC],LinOp.Op,LinOp.rhs)
    end
    if :nNodes in bc(prob)
      scalarNeumann!(mesh.xy,mesh.cm,mesh.order,
                     prob.bcid[:nNodes],prob.bcval[:nBC],LinOp.rhs)
    end

  # ADVECTION-DIFFUSION 2D
  elseif OpType == :AdvDiff2D
    if :dNodes in bc(prob)
      scalarDirichlet!(prob.bcid[:dNodes],prob.bcval[:dBC],LinOp.Op,LinOp.rhs)
    end

    if :nNodes in bc(prob)
      scalarNeumann!(mesh.xy,mesh.cm,mesh.order,
                     prob.bcid[:nNodes],prob.bcval[:nBC],LinOp.rhs)
    end

    if :rNodes in bc(prob)
      AdvDiffRobin!(mesh.xy,mesh.cm,prob.bcid[:rNodes],
                    LinOp.Op,LinOp.rhs,hcat(param.u,param.v))
    end

  # STOKES 2D or BRINKMAN 2D or BRINKMAN MULTIPHASE 2D
  elseif OpType == :Stokes2D || OpType == :Brinkman2D || OpType == :BrinkmanMP2D
    if :dUNodes in bc(prob)
      fluidDirichlet!(mesh,prob,LinOp)
    end

    if !(:nNodes in bc(prob))
      LinOp.Op,LinOp.rhs = fluidNeumann(mesh,prob,LinOp)
    end

  # ADVECTION-DIFFUSION AXISYMMETRIC
  elseif OpType == :AdvDiffAS || OpType == :PoissonAS
    if :dNodes in bc(prob)
      scalarDirichlet!(prob.bcid[:dNodes],prob.bcval[:dBC],LinOp.Op,LinOp.rhs)
    end

    if :nNodes in bc(prob)
      scalarNeumann!(mesh.xy,mesh.cm,mesh.order,
                    prob.bcid[:nNodes],prob.bcval[:nBC],LinOp.rhs)
    end

    if :rNodes in bc(prob)
      AdvDiffRobin!(mesh.xy,mesh.cm,prob.bcid[:rNodes],
                    LinOp.Op,LinOp.rhs,prob.wind)
    end

  # STOKES AXISYMMETRIC
  elseif OpType == :StokesAS
    if :dUNodes in bc(prob)
      fluidDirichlet!(mesh,prob,LinOp)
    end

    if !(:nNodes in bc(prob))
      LinOp.Op,LinOp.rhs = fluidNeumann(mesh,prob,LinOp)
    end
  end
end

##############################
######## Solve System ########
##############################

function LinearSolve!(sol::AbstractSolution,mesh::AbstractMesh,
                      LinOp::LinearOperator,prob::AbstractProblem,
                      param::AbstractParameter,OpType::Symbol)
  # POISSON 2D
  if OpType == :Poisson2D || OpType == :PoissonAS
    sol.u = GaussElimSolve(LinOp)

  # DARCY 2D
  elseif OpType == :Darcy2D
    println("LinearSolve! not defined for Darcy")
    return FluidSolution(u,v,p)

  # ADVECTION-DIFFUSION 2D
  elseif OpType == :AdvDiff2D || OpType == :AdvDiffAS
    sol.u = GaussElimSolve(LinOp)

  # STOKES 2D or BRINKMAN 2D or BRINKMAN MULTIPHASE 2D
  elseif ( OpType == :Stokes2D     || OpType == :Brinkman2D ||
           OpType == :BrinkmanMP2D)
    fluidSolve!(sol,mesh,prob,LinOp)

  # STOKES AS requires row-normalization preconditioner!
  elseif OpType == :StokesAS
    #rowNormalizePreconditioner!(LinOp)
    fluidSolve!(sol,mesh,prob,LinOp)
  end

  return nothing
end

function LinearSolve(mesh::AbstractMesh,
                      LinOp::LinearOperator,prob::AbstractProblem,
                      param::AbstractParameter,OpType::Symbol)
  # POISSON 2D
  if OpType == :Poisson2D || OpType == :PoissonAS
    u = GaussElimSolve(LinOp)
    #u = cg(LinOp.Op,LinOp.rhs)

  return ScalarSolution(u)

  # DARCY 2D
  elseif OpType == :Darcy2D
    p = GaussElimSolve(LinOp)

    # preconditioner=
    #pc = AMGPreconditioner{SmoothedAggregation}(LinOp.Op)

    #p = IterativeSolvers.cg(LinOp.Op, LinOp.rhs;verbose=false,Pl=pc)

    #p = GaussElimSolve(LinOp)
    #u,v = DarcyVelocity(mesh,param.α,p)
  #return FluidSolution(u,v,p)
  return ScalarSolution(p)

  # ADVECTION-DIFFUSION 2D
  elseif OpType == :AdvDiff2D || OpType == :AdvDiffAS
    #pc = DiagonalPreconditioner(LinOp.Op)

    #u = gmres(LinOp.Op, LinOp.rhs;verbose=true,Pl=pc)
    u = GaussElimSolve(LinOp)
  return ScalarSolution(u)

  # STOKES 2D or BRINKMAN 2D or BRINKMAN MULTIPHASE 2D
  elseif ( OpType == :Stokes2D     || OpType == :Brinkman2D ||
    OpType == :BrinkmanMP2D)
    return fluidSolve(mesh,prob,LinOp)

  # STOKES AS requires row-normalization preconditioner!
  elseif OpType == :StokesAS
    #rowNormalizePreconditioner!(LinOp)
    return fluidSolve(mesh,prob,LinOp)
  end
end

function GaussElimSolve(LinOp::LinearOperator)
  U = LinOp.Op\LinOp.rhs

  return U
end

"""
  fluidSolve(xy,cm,nNodes,Stiff,F;PRINT=false) -> u,v,p

Computes and decomposes the array-solution U into vector velocity [u,v] and scalar pressure p
"""
function fluidSolve(mesh::FluidMesh,prob::Problem,
                    LinOp::AbstractLinearOperator;PRINT=false)
  # preconditioner
  #pc = IterativeSolvers.Identity()
  #N,_  = size(LinOp.Op) #fit_candidates
  #ml = AlgebraicMultigrid.ruge_stuben(LinOp.Op)
  #P  = AlgebraicMultigrid.aspreconditioner(ml)
  #pc   = MyLinOp(Pmap)
  #pc   = Preconditioners.AMGPreconditioner{SmoothedAggregation}(LinOp.Op)
  #soln = IterativeSolvers.gmres(LinOp.Op, LinOp.rhs;maxiter=100,verbose=false,Pl=P)

  #soln = qr(LinOp.Op)\(LinOp.rhs)  # qr factorization does better on MMS
  #rowNormalizePreconditioner!(LinOp)
  soln = lu(LinOp.Op)\(LinOp.rhs)   # lu factorization does better on external domains
                                    # I HAVE NO IDEA WHY

  # decompose solution
  NumUnodes = length(mesh.xy)
  u = soln[1:NumUnodes]
  v = soln[NumUnodes+1:2*NumUnodes]
  p = soln[2*NumUnodes+1:end]

  # adjust for enclosed boundary condition
  if !(:nNodes in bc(prob))
    if PRINT
      println("adjusting for enclosed boundary condition")
    end
    p=p[1:end-1]
  end

  return FluidSolution(u,v,p)
end

"""
  fluidSolve(xy,cm,nNodes,Stiff,F;PRINT=false) -> u,v,p

Computes and decomposes the array-solution U into vector velocity [u,v] and scalar pressure p
"""
function fluidSolve!(sol::FluidSolution,
                     mesh::FluidMesh,prob::Problem,
                     LinOp::AbstractLinearOperator;PRINT=false)
  soln = lu(LinOp.Op)\LinOp.rhs

  # decompose solution
  NumUnodes = length(mesh.xy)
  sol.u = soln[1:NumUnodes]
  sol.v = soln[NumUnodes+1:2*NumUnodes]
  sol.p = soln[2*NumUnodes+1:end]

  # adjust for enclosed boundary condition
  if !(:nNodes in bc(prob))
    if PRINT
      println("adjusting for enclosed boundary condition")
    end
    sol.p=sol.p[1:end-1]
  end

  return nothing
end

function DarcyVelocity(mesh,alpha::Vector{Float64},p::Vector{Float64})
  Nel    = length(mesh.cm)
  Nnodes = length(mesh.xy)

  nodeCounter = zeros(Int,Nnodes)

  u = zeros(Float64,Nnodes)
  v = zeros(Float64,Nnodes)

  xNodes = zeros(Float64,9)
  yNodes = zeros(Float64,9)
  αNodes = zeros(Float64,9)
  pNodes = zeros(Float64,9)

  w,s,t = GaussQuadPoints2D(3)

  # for each element
  for el=1:Nel
    c = mesh.cm[el].NodeList
    for i=1:9
      xNodes[i] = mesh.xy[c[i]].x
      yNodes[i] = mesh.xy[c[i]].y
      αNodes[i] = alpha[c[i]]
      pNodes[i] = p[c[i]]
    end

    # for each node in element (assumes quadratic basis)
    for i=1:9
      phi,dphidx,dphidy,jac = derivShape2D(s[i],t[i],xNodes,yNodes,2)
      dpdx = shapeEval(pNodes,dphidx)
      dpdy = shapeEval(pNodes,dphidy)
      αg = shapeEval(αNodes,phi)

      u[c[i]] += -αg*dpdx/jac
      v[c[i]] += -αg*dpdy/jac
    end
    nodeCounter[c] .+= 1
  end

  # take average
  for i=1:Nnodes
    u[i] /= nodeCounter[i]
    v[i] /= nodeCounter[i]
  end

  return u,v
end
function DarcyVelocity(mesh,alpha::Float64,p::Vector{Float64})
  return DarcyVelocity(mesh,alpha*ones(Float64,length(mesh.xy)),p)
end

function rowNormalizePreconditioner!(LinOp)
    # apply preconditioner
    rows,cols,_ = findnz(LinOp.Op)
    nrows = length(rows)

    # compute row norms
    counter = 0
    alreadyDone = []
    for ti=1:nrows
      if rows[ti] in alreadyDone

      else
        normval = maximum(LinOp.Op[rows[ti],:])

        # normalize Matrix
        LinOp.Op[rows[ti],:] ./= normval

        # normalize forcing -- this is going twice!!!
        LinOp.rhs[rows[ti]] /= normval

        push!(alreadyDone,rows[ti])
        counter += 1
      end
    end

  return nothing
end


text:
module Models

  using ITensors

  export Model, critical_point, local_boltzmann_weight

  struct Model{model} end
  Model(s::AbstractString) = Model{Symbol(s)}()
  macro Model_str(s)
    :(Model{$(Expr(:quote, Symbol(s)))})
  end

  #mpo_itensor(m::String, args...; kwargs...) =
  #  mpo_itensor(Model(m), args...; kwargs...)

  local_boltzmann_weight(m::String, args...; kwargs...) =
    local_boltzmann_weight(Model(m), args...; kwargs...)

  f(λ₊, λ₋) = [
    (λ₊ + λ₋)/2 (λ₊ - λ₋)/2
    (λ₊ - λ₋)/2 (λ₊ + λ₋)/2
  ]

  function sqrt_bond_matrix(; β::Real, J::Real=1.0)
    # Alternative method
    #Q = [exp(β * J) exp(-β * J); exp(-β * J) exp(β * J)]
    #return √Q
    λ₊ = √(exp(β * J) + exp(-β * J))
    λ₋ = √(exp(β * J) - exp(-β * J))
    return f(λ₊, λ₋)
  end

  # The local Boltzmann weight for the Ising
  # model in 1 dimension
  function local_boltzmann_weight(
    ::Model"ising",
    ::Val{1};
    β::Real,
    J::Real=1.0,
    sz::Bool=false,
  )
    d = 2 # local dimension of the Ising local Boltzmann factor
    s, s′ = Index.((d, d))
    T = ITensor(s, s′)
    for i in 1:d
      T[i, i] = 1.0
    end
    if sz
      T[1, 1] = -T[1, 1]
    end
    s̃, s̃′ = sim.((s, s′))
    T̃ = T * δ(s, s̃) * δ(s′, s̃′)
    sqrtQ = sqrt_bond_matrix(; β=β, J=J)
    @show sqrtQ
    X = itensor(vec(sqrtQ), s̃, s)
    X′ = itensor(vec(sqrtQ), s̃′, s′)
    return array(permute(T̃ * X′ * X, s, s′))
  end

  # The local Boltzmann weight for the Ising
  # model in 2 dimensions
  function local_boltzmann_weight(
    ::Model"ising",
    ::Val{2};
    β::Real,
    J::Real=1.0,
    sz::Bool=false,
  )
    d = 2 # local dimension of the Ising local Boltzmann factor
    sₕ, sₕ′ = Index.((d, d))
    sᵥ, sᵥ′ = Index.((d, d))
    @assert dim(sₕ) == dim(sᵥ)
    d = dim(sₕ)
    T = ITensor(sₕ, sₕ′, sᵥ, sᵥ′)
    for i in 1:d
      T[i, i, i, i] = 1.0
    end
    if sz
      T[1, 1, 1, 1] = -T[1, 1, 1, 1]
    end
    s̃ₕ, s̃ₕ′, s̃ᵥ, s̃ᵥ′ = sim.((sₕ, sₕ′, sᵥ, sᵥ′))
    T̃ = T * δ(sₕ, s̃ₕ) * δ(sₕ′, s̃ₕ′) * δ(sᵥ, s̃ᵥ) * δ(sᵥ′, s̃ᵥ′)
    X = sqrt_bond_matrix(; β, J)
    Xₕ = itensor(vec(X), s̃ₕ, sₕ)
    Xₕ′ = itensor(vec(X), s̃ₕ′, sₕ′)
    Xᵥ = itensor(vec(X), s̃ᵥ, sᵥ)
    Xᵥ′ = itensor(vec(X), s̃ᵥ′, sᵥ′)
    return array(permute(T̃ * Xₕ′ * Xᵥ′ * Xₕ * Xᵥ, sₕ, sₕ′, sᵥ, sᵥ′))
  end

  function mpo_itensor(
    m::Model,
    pair_sₕ::Pair{<:Index,<:Index},
    pair_sᵥ::Pair{<:Index,<:Index};
    kwargs...
  )
    sₕ, sₕ′ = pair_sₕ
    sᵥ, sᵥ′ = pair_sᵥ
    return itensor(mpo_array(m; kwargs...), sₕ, sₕ′, sᵥ, sᵥ′)
  end

  ## function ising_mpo(sₕ::Index, sᵥ::Index, args...; kwargs...)
  ##   return ising_mpo(sₕ => sₕ', sᵥ => sᵥ', args...; kwargs...)
  ## end

  critical_point(::Model"ising") = 0.5 * log(√2 + 1)
  #const βc = critical_point(Model("ising"))

  function free_energy(::Model"ising"; β::Real, J::Real=1.0)
    k = β * J
    c = cosh(2 * k)
    s = sinh(2 * k)
    xmin = 0.0
    xmax = π
    integrand(x) = log(c^2 + √(s^4 + 1 - 2 * s^2 * cos(x)))
    integral, err = quadgk(integrand, xmin, xmax)::Tuple{Float64,Float64}
    return -(log(2) + integral / π) / (2 * β)
  end

  function magnetization(::Model"ising"; β::Real)
    β > βc && return (1 - sinh(2 * β)^(-4))^(1 / 8)
    return 0.0
  end

end


text:
@safetestset "isboolean" begin
    using Validator

    # loose=false
    @test Validator.isboolean("true") === true
    @test Validator.isboolean("false") === true
    @test Validator.isboolean("0") === true
    @test Validator.isboolean("1") === true

    @test Validator.isboolean("1.0") === false
    @test Validator.isboolean("0.0") === false
    @test Validator.isboolean("true ") === false
    @test Validator.isboolean("False") === false
    @test Validator.isboolean("True") === false
    @test Validator.isboolean("yes") === false

    # loose=true
    @test Validator.isboolean("true", true) === true
    @test Validator.isboolean("True", true) === true
    @test Validator.isboolean("TRUE", true) === true
    @test Validator.isboolean("false", true) === true
    @test Validator.isboolean("False", true) === true
    @test Validator.isboolean("FALSE", true) === true
    @test Validator.isboolean("1", true) === true
    @test Validator.isboolean("0", true) === true
    @test Validator.isboolean("yes", true) === true
    @test Validator.isboolean("Yes", true) === true
    @test Validator.isboolean("YES", true) === true
    @test Validator.isboolean("no", true) === true
    @test Validator.isboolean("No", true) === true
    @test Validator.isboolean("NO", true) === true

    @test Validator.isboolean("1.0", true) === false
    @test Validator.isboolean("0.0", true) === false
    @test Validator.isboolean("true ", true) === false
    @test Validator.isboolean(" false", true) === false
end


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

mutable struct TrafficAnalyticsProperties <: SwaggerModel
    networkWatcherFlowAnalyticsConfiguration::Any # spec type: Union{ Nothing, TrafficAnalyticsConfigurationProperties } # spec name: networkWatcherFlowAnalyticsConfiguration

    function TrafficAnalyticsProperties(;networkWatcherFlowAnalyticsConfiguration=nothing)
        o = new()
        validate_property(TrafficAnalyticsProperties, Symbol("networkWatcherFlowAnalyticsConfiguration"), networkWatcherFlowAnalyticsConfiguration)
        setfield!(o, Symbol("networkWatcherFlowAnalyticsConfiguration"), networkWatcherFlowAnalyticsConfiguration)
        o
    end
end # type TrafficAnalyticsProperties

const _property_map_TrafficAnalyticsProperties = Dict{Symbol,Symbol}(Symbol("networkWatcherFlowAnalyticsConfiguration")=>Symbol("networkWatcherFlowAnalyticsConfiguration"))
const _property_types_TrafficAnalyticsProperties = Dict{Symbol,String}(Symbol("networkWatcherFlowAnalyticsConfiguration")=>"TrafficAnalyticsConfigurationProperties")
Base.propertynames(::Type{ TrafficAnalyticsProperties }) = collect(keys(_property_map_TrafficAnalyticsProperties))
Swagger.property_type(::Type{ TrafficAnalyticsProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TrafficAnalyticsProperties[name]))}
Swagger.field_name(::Type{ TrafficAnalyticsProperties }, property_name::Symbol) =  _property_map_TrafficAnalyticsProperties[property_name]

function check_required(o::TrafficAnalyticsProperties)
    true
end

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


text:
module Hulls

export ConicHull, verify, validate, create_hull, init_hull!, ftype

using ..Common
using ..Dets
using ..Verify
import ..Common: nconic, gtype, get_canonical_winding, indexof
import ..Common: dot
import ..Common: add_facet!, del_facet!, mark_facet!, ismarked
import ..Common: nbof, nbsof, generatorof, generatorsof, set_nb!

# --------------------------------- AFacet -----------------------------------

global numfacets = 0

immutable AFacet{G} <: Facet
    generators::Vector{G}
    links::Vector{Union(AFacet{G},Nothing)}
    positive::Bool
    id::Int

    function AFacet(generators, positive)
        global numfacets
        NF = length(generators)
        f = new(G[generators...], fill!(Array(Union(AFacet{G},Nothing), NF), nothing),
                                        positive, numfacets += 1)
        f
    end
end
function AFacet{F<:AFacet}(face::Face{F}, generator::Generator)
    facet = face.parent
    generators = copy(facet.generators)
    generators[face.k] = generator
    F(generators, facet.positive)
end

gtype{G}( ::Type{AFacet{G}}) = G
nconic(facet::AFacet) = length(facet.generators) + 1

nbof(f::AFacet, k::Int)                    = f.links[k]
set_nb!(f::AFacet, new_nb::AFacet, k::Int) = (f.links[k] = new_nb; nothing)
generatorof(f::AFacet, k::Int)             = f.generators[k]
nbsof(f::AFacet)        = f.links
generatorsof(f::AFacet) = f.generators

# Base.show(io::IO, f::AFacet) = print(io, "Facet$(f.id)")
function Base.show(io::IO, f::AFacet)
    print(io, "Facet$(f.id)(", f.generators, ", [")
    for (k, nb) in enumerate(f.links)
        print(io, "Facet$(nb.id)")
        if k < nfacet(f); print(io, ", "); end
    end
    print(io, "])");
end

indexof{F<:AFacet}(f::F, link::F)    = indexof(f.links, link)
indexof{G}(f::AFacet{G}, g::G) = indexof(f.generators, g)

function get_canonical_winding(facet::AFacet)
    perm = sortperm(facet.generators, by=g->g.id)
    generators = facet.generators[perm]
    even = facet.positive $ isoddperm(perm)
    (generators, even)
end

# -------------------------------- Dominance ---------------------------------

dot(f::AFacet, g::Generator)   = xorsign(det(f.generators..., g), !f.positive)
xorsign(x::Number, flip::Bool) = flip ? -x : x

hassign(x::Number, positive::Bool) = positive ? x > 0 : x < 0
function antidominates_replaced(g::Generator, face::Face, replacement::Generator)
    gs = copy(face.parent.generators)
    gs[face.k] = replacement
    hassign(det(gs..., g), !face.parent.positive)
end

# -------------------------------- ConicHull ---------------------------------

type ConicHull{F<:AFacet,G<:Generator} <: Hull
    nc::Int
    generators::Vector{G}
    facets::Set{F}
    marked::Set{F}

#    ConicHull() = (@assert nconic(F) == nconic(G); new())
    ConicHull(nc::Int) = new(nc)
    function ConicHull(nc::Int, generators, facets)
        hull = ConicHull{F,G}(nc)
        init_hull!(hull, generators, facets)
        hull
    end
end

function init_hull!{F,G}(hull::ConicHull{F,G}, generators, facets)
    hull.generators = G[generators...]
    hull.facets = Set{F}(facets)
    hull.marked = Set{F}()
end

hulltype{G<:Generator}(::Type{G}) = ConicHull{AFacet{G},G}

ftype{F,G}(::Type{ConicHull{F,G}}) = F
ftype{F,G}(::ConicHull{F,G})       = F

gtype{F,G}(::Type{ConicHull{F,G}}) = G
gtype{F,G}(::ConicHull{F,G})       = G

nconic(hull::ConicHull) = hull.nc

function add_facet!(h::ConicHull, face::Face, g::Generator)
    f = AFacet(face, g)
    push!(h.facets, f) # todo: add in when it doesn't signify marked status anymore
    f
end
del_facet!(h::ConicHull, facet::Facet) = (pop!(h.facets, facet); nothing)
function mark_facet!(h::ConicHull, facet::Facet, mark::Bool)
    mark ? push!(h.marked, facet) : delete!(h.marked, facet)
    nothing
end
ismarked(h::ConicHull, facet::Facet) = facet in h.marked

function verify(hull::ConicHull)
    NF = nfacet(hull)
    for facet in hull.facets
        @assert length(facet.links) == length(facet.generators) == NF
        if length(Set(facet.links)) != NF
            error("facet $facet has duplicate neighbors")
        end
        if length(Set(facet.generators)) != NF
            error("facet $facet has duplicate generators")
        end
        for (nb, generator) in zip(facet.links, facet.generators)
            # @assert !is(nb, nothing)
            if !(nb in hull.facets)
                error("facet $facet links to facet $nb, which is not in the hull.")
            end
            # Check that neighbor is not dominated by the opposing generator
            @assert !dominates(generator, nb)

            # Check that there is a link back
            nb_g = generatorof(nb, facet)
            # Check that the linked facets have all generators in common except the opposite ones
            link_gs = Set(facet.generators); pop!(link_gs, generator)
            link_nb_gs = Set(nb.generators); pop!(link_nb_gs, nb_g)
            @assert link_gs == link_nb_gs
            # Check that the opposite generators are different
            @assert nb_g != generator
        end
    end
end

function validate(hull::ConicHull)
    links = create_nb_dict(ftype(hull), hull.facets)
    for facet in hull.facets
        for (g,nb) in zip(facet.generators, facet.links)
            @assert nb === get_nb(links, facet, g)
        end
    end
end

create_simplex(NC, F, G, gs::Matrix) = create_simplex(NC, F, G, [gs[:,k] for k in 1:size(gs,2)])
function create_simplex(NC, F, G, gs::Vector)
    @assert length(gs) == NC
    gs = [G(x) for x in gs]
    for g in gs; (@assert nconic(g) == NC); end
    d = det(gs...)
    @assert d != 0
    odd = d < 0

    fs = F[F(except_index(gs, kf), isodd(NC-kf) $ odd) for kf in 1:NC]

    for (kf, f) in enumerate(fs)
        links = except_index([1:NC...], kf)
        for (kn, nb) in enumerate(links)
            f.links[kn] = fs[nb]
        end
    end

    gs, fs
end

create_hull{F,G}(nc::Int, H::Type{ConicHull{F,G}}) = H(nc)
function create_hull{F,G}(nc::Int, H::Type{ConicHull{F,G}}, generators)
    hull = create_hull(nc, H)
    init_hull!(hull, generators)
    hull
end

init_hull!(hull::ConicHull, gs::Matrix) = init_hull!(hull, [gs[:,k] for k in 1:size(gs,2)])
function init_hull!(hull::ConicHull, generators::Vector)
    NC = nconic(hull)
    gs, fs = create_simplex(NC, ftype(hull), gtype(hull), generators)
    init_hull!(hull, gs, fs)
    copy(gs)
end

end # module


text:
using MatFun, Base.Test

@testset "aaa" begin
	Z = collect(linspace(-3.1, 3.1, 200))
	r, pol, res, zer, z, f, w, errvec = aaa(gamma, Z)
	@test r(Z) ≈ gamma.(Z)
	@test any(pol .≈ -1) && any(pol .≈ -2) && any(pol .≈ -3)
	@test norm(gamma.(zer), Inf) <= 0.1

	srand(789)
	poles = randn(10) + im*randn(10)
	residues = randn(10) + im*randn(10)
	func = (x) -> sum(residues./(x .- poles))
	Z = MatFun.lpdisk(sqrt(2.0), 400)
	F = func.(Z)
	r, pol, res, zer, z, f, w, errvec = aaa(F, Z)
	@test r(Z) ≈ F
	@test sort(real(pol)) ≈ sort(real(poles)) && sort(imag(pol)) ≈ sort(imag(poles))
	@test sort(real(res)) ≈ sort(real(residues)) && sort(imag(res)) ≈ sort(imag(residues))
	@test norm(func.(zer), Inf) <= 0.1
end


text:
module Alea

include("internal/Internal.jl")
include("math/Math.jl")
include("statistics/Statistics.jl")
include("quadrature/Quadrature.jl")
include("bases/Bases.jl")
include("gpc/GPC.jl")

using Reexport
@reexport using .Math
@reexport using .Statistics
@reexport using .Quadrature
@reexport using .Bases
@reexport using .GPC

end


text:
## Low level gemm! call with pointers
## Borrowed from Knet.jl

using LinearAlgebra
using LinearAlgebra.BLAS: libblas, BlasInt,  @blasfunc

# C := alpha*op(A)*op(B) + beta*C, where:
# op(X) is one of op(X) = X, or op(X) = XT, or op(X) = XH,
# alpha and beta are scalars,
# A, B and C are matrices:
# op(A) is an m-by-k matrix,
# op(B) is a k-by-n matrix,
# C is an m-by-n matrix.

for (gemm, elty) in ((:dgemm_,:Float64), (:sgemm_,:Float32))
    @eval begin
        function gemm!(transA::Char, transB::Char, M::Int, N::Int, K::Int, alpha::($elty), A::Ptr{$elty}, B::Ptr{$elty}, beta::($elty), C::Ptr{$elty})
            if transA=='N'; lda=M; else; lda=K; end
            if transB=='N'; ldb=K; else; ldb=N; end
            ldc = M;
            ccall((@blasfunc($(gemm)), libblas), Nothing,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                   Ref{BlasInt}),
                  transA, transB, M, N, K,
                  alpha, A, lda, B, ldb, beta, C, ldc)
        end
    end
end


text:
using CSV
using DataFrames
using WordCloud

df = CSV.File(pkgdir(WordCloud)*"/res/guxiang_frequency.txt", header=false)|> DataFrame;
words = df[!, "Column2"]
weights = df[!, "Column3"]

wc = wordcloud(words, weights, density=0.65)
gifdirectory = "guxiang_animation"
#eval# try rm("guxiang_animation", force=true, recursive=true) catch end
generate_animation!(wc, 100, outputdir=gifdirectory)
println("results are saved in guxiang_animation")
wc
#eval# runexample(:animation)
#md# ![](guxiang_animation/result.gif)


text:
#
# Options
#
@doc raw"""
    ParticleSwarmOptions{P,T} <: Options

Describes a particle swarm optimizing algorithm, with

# Fields

* `x0` – a set of points (of type `AbstractVector{P}`) on a manifold as initial particle positions
* `velocity` – a set of tangent vectors (of type `AbstractVector{T}`) representing the velocities of the particles
* `inertia` – (`0.65`) the inertia of the patricles
* `social_weight` – (`1.4`) a social weight factor
* `cognitive_weight` – (`1.4`) a cognitive weight factor
* `stopping_criterion` – (`[`StopAfterIteration`](@ref)`(500) | `[`StopWhenChangeLess`](@ref)`(1e-4)`)
  a functor inheriting from [`StoppingCriterion`](@ref) indicating when to stop.
* `retraction_method` – (`default_retraction_method(M)`) the rectraction to use
* `inverse_retraction_method` - (`default_inverse_retraction_method(M)`) an inverse retraction to use.
* `vector_transport_method` - (`default_vector_transport_method(M)`) a vector transport to use

# Constructor

    ParticleSwarmOptions(M, x0, velocity; kawrgs...)

construct a particle swarm Option for the manifold `M` starting at point `x0` with velocities `x0`,
where the manifold is used within the defaults of the other fields mentioned above,
which are keyword arguments here.

# See also
[`particle_swarm`](@ref)
"""
mutable struct ParticleSwarmOptions{
    TX<:AbstractVector,
    TG,
    TVelocity<:AbstractVector,
    TParams<:Real,
    TStopping<:StoppingCriterion,
    TRetraction<:AbstractRetractionMethod,
    TInvRetraction<:AbstractInverseRetractionMethod,
    TVTM<:AbstractVectorTransportMethod,
} <: Options
    x::TX
    p::TX
    g::TG
    velocity::TVelocity
    inertia::TParams
    social_weight::TParams
    cognitive_weight::TParams
    stop::TStopping
    retraction_method::TRetraction
    inverse_retraction_method::TInvRetraction
    vector_transport_method::TVTM
    @deprecate ParticleSwarmOptions(
        x0::AbstractVector,
        velocity::AbstractVector,
        inertia=0.65,
        social_weight=1.4,
        cognitive_weight=1.4,
        stopping_criterion::StoppingCriterion=StopWhenAny(
            StopAfterIteration(500), StopWhenChangeLess(10.0^(-4))
        ),
        retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
        inverse_retraction_method::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
        vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
    ) ParticleSwarmOptions(
        DefaultManifold(2),
        x0,
        velocity;
        inertia=inertia,
        social_weight=social_weight,
        cognitive_weight=cognitive_weight,
        stopping_criterion=stopping_criterion,
        retraction_method=retraction_method,
        inverse_retraction_method=inverse_retraction_method,
        vector_transport_method=vector_transport_method,
    )
    function ParticleSwarmOptions(
        M::AbstractManifold,
        x0::AbstractVector,
        velocity::AbstractVector;
        inertia=0.65,
        social_weight=1.4,
        cognitive_weight=1.4,
        stopping_criterion::StoppingCriterion=StopAfterIteration(500) |
                                              StopWhenChangeLess(1e-4),
        retraction_method::AbstractRetractionMethod=default_retraction_method(M),
        inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method(
            M
        ),
        vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(
            M
        ),
    )
        o = new{
            typeof(x0),
            eltype(x0),
            typeof(velocity),
            typeof(inertia + social_weight + cognitive_weight),
            typeof(stopping_criterion),
            typeof(retraction_method),
            typeof(inverse_retraction_method),
            typeof(vector_transport_method),
        }()
        o.x = x0
        o.p = deepcopy(x0)
        o.velocity = velocity
        o.inertia = inertia
        o.social_weight = social_weight
        o.cognitive_weight = cognitive_weight
        o.stop = stopping_criterion
        o.retraction_method = retraction_method
        o.inverse_retraction_method = inverse_retraction_method
        o.vector_transport_method = vector_transport_method
        return o
    end
end


text:
module Unidecode
    using REPL.REPLCompletions: emoji_symbols, latex_symbols
    using Unicode: graphemes
    const emoji_dict = Dict((v => k[2:end]) for (k, v) in emoji_symbols)
    const latex_dict = Dict((v => k[2:end]) for (k, v) in latex_symbols)
    function decode_emoji(x)
        string((get(emoji_dict, c, c) for c in graphemes(x))...)
    end
    function decode_latex(x)
        string((get(latex_dict, c, c) for c in graphemes(x))...)
    end
    function unidecode(x)
        decode_emoji(decode_latex(x))
    end
    export decode_emoji, decode_latex, unidecode
end


text:
using PackageCompiler

userimg = Pkg.dir("Makie", "test", "makie_precompile.jl")
image_path = PackageCompiler.sysimg_folder()
PackageCompiler.build_sysimg(image_path, userimg, compilecache = "yes")
isfile(backup) || mv(syspath, backup)
cp(imgfile, syspath)


text:
using Reinforce.MountainCarEnv

@testset "MountainCarEnv" begin

@info "Reinforce.MountainCarEnv"

env = MountainCar()

i = 0
for sars in Episode(env, RandomPolicy())
    i += 1
end

@test i > 0

end  # @testset "MountainCarEnv"


text:
module CImGuiExtensions

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGui: ImVec2, ImVec4, IM_COL32, ImU32
using CSV
using DataFrames
using Printf

include("abstract_types.jl")
include("util.jl")
include("status.jl")
include("path.jl")
include("context/export_context.jl")
include("context/import_context.jl")

include("control/control.jl")

include("layout.jl")
include("axis.jl")
include("tickmark.jl")

include("view/plotlines_display_properties.jl")
include("model/plotlines_model.jl")
include("control/plotlines_control.jl")

include("view/file_dialog_display_properties.jl")
include("model/file_dialog_model.jl")
include("control/file_dialog_control.jl")

include("view/nested_interval_display_properties.jl")
include("model/nested_interval_model.jl")
include("control/nested_interval_control.jl")

include("view/label_intervals_display_properties.jl")
include("model/labelled_interval_model.jl")
include("control/label_interval_control.jl")

include("schemas.jl")

include("operation/csv_import.jl")
include("operation/csv_export.jl")
include("operation/labelled_intervals_import.jl")
include("operation/labelled_intervals_export.jl")

include("context/plotlines_context.jl")

include("context/nested_interval_context.jl")
include("context/label_interval_context.jl")

include("control/truncated_plot_control.jl")
include("context/truncated_plot_context.jl")

#include("displayproperties.jl")
# include("tickmark.jl")
# include("axis.jl")
# include("schemas.jl")
# include("image_importer.jl")
# include("filedialog.jl")
# include("control.jl")
# include("plotlines.jl")
# include("nestedinterval.jl")
# include("labelintervals.jl")
# include("truncatedplot.jl")
# include("import.jl")
# include("export.jl")

        # Abstract types
export  AbstractModel,
        AbstractDialogModel,
        AbstractStatus,
        AbstractDisplayProperties,
        AbstractControl,
        AbstractOperation,
        AbstractImporter,
        AbstractExporter,
        AbstractVendor,
        AbstractProduct,
        AbstractDialogControl,
        AbstractData,
        AbstractLayout,
        AbstractContext,
        AbstractPlotContext,
        AbstractSchema,
        AbstractModelViewControl,
        # Concrete types
        Path,
        FileDialogModel,
        EnabledStatus,
        DisabledStatus,
        ConfirmedStatus,
        UnconfirmedStatus,
        FileDialogDisplayProperties,
        NestedIntervalDisplayProperties,
        FileDialogControl,
        GenericVendor,
        GenericProduct,
        CSVSchema,
        CSVImporter,
        ImageImporter,
        CSVExporter,
        RectangularLayout,
        NestedIntervalControl,
        NestedIntervalContext,
        PlotlinesControl,
        PlotlinesModel,
        PlotlinesDisplayProperties,
        ImportContext,
        ExportContext,
        PlotContext,
        NestedInterval,
        Tickmark,
        Axis,
        TruncatedPlotContext,
        LabelledIntervals,
        LabelledInterval,
        LabelIntervalContext,
        LabelIntervalControl,
        LabelIntervalDisplayProperties,
        IntervalLabels,
        # Functions
        get_path,
        set_path!,
        get_description,
        set_description!,
        get_labelled_intervals,
        set_labelled_intervals!,
        get_label,
        set_label!,
        get_nested_interval,
        set_nested_interval!,
        get_start,
        set_start!,
        get_stop,
        set_stop!,
        get_interval,
        set_interval!,
        copy,
        get_data,
        set_data!,
        isenabled,
        enable!,
        disable!,
        isrunning,
        load_dataframe,
        get_caption,
        set_caption!,
        get_action,
        set_action!,
        get_width,
        set_width!,
        get_height,
        set_height!,
        get_id,
        set_id!,
        get_caption,
        set_caption!,
        get_col,
        set_col!,
        get_textcol,
        set_textcol!,
        get_layout,
        set_layout!,
        get_padding,
        set_padding!,
        get_plotcontex,
        set_plotcontext!,
        get_xtick,
        set_xtick!,
        get_ytick,
        set_ytick!,
        get_yaxis,
        set_yaxis!,
        get_xaxis,
        set_xaxis!,
        get_thickness,
        set_thickness!,
        get_pos,
        set_pos!,
        get_directory,
        set_directory!,
        get_filename,
        set_filename!,
        get_interpreter,
        set_interpreter!,
        get_length,
        set_length!,
        get_spacing,
        set_spacing!,
        get_interpreter,
        set_interpreter!,
        get_captioner,
        set_captioner!,
        stretch_linearly
end # module


text:
function estimator_extrema(x::MC{N,T}, y::MC{N,T}, s, dP) where {N,T}
    xcv = x.cv;   xcvg = x.cv_grad
    ycv = y.cv;   ycvg = y.cv_grad
    xccn = -x.cc; xccgn = -x.cc_grad
    yccn = -y.cc; yccgn = -y.cc_grad
    t3 = affine_expand_del(dP, xcv,  xcvg,  s)
    t4 = affine_expand_del(dP, ycv,  ycvg,  s)
    s3 = affine_expand_del(dP, xccn, xccgn, s)
    s4 = affine_expand_del(dP, yccn, yccgn, s)
    return t3, t4, s3, s4
end

function estimator_under(xv, yv, x::MC{N,T}, y::MC{N,T}, s, dp, dP, p_rel, p_diam) where {N,T}
    xcv = x.cv;   xcvg = x.cv_grad
    ycv = y.cv;   ycvg = y.cv_grad
    u1cv = affine_expand_del(dp, xcv, xcvg, s)
    u2cv = affine_expand_del(dp, ycv, ycvg, s)
    return u1cv, u2cv, xcvg, ycvg
end

function estimator_over(xv, yv, x::MC{N,T}, y::MC{N,T}, s, dp, dP, p_rel, p_diam) where {N,T}
    xccn = -x.cc; xccgn = -x.cc_grad
    yccn = -y.cc; yccgn = -y.cc_grad
    v1ccn = affine_expand_del(dp, xccn, xccgn, s)
    v2ccn = affine_expand_del(dp, yccn, yccgn, s)
    return v1ccn, v2ccn, xccgn, yccgn
end

include(joinpath(@__DIR__, "enumeration.jl"))
include(joinpath(@__DIR__, "affine_arithmetic.jl"))

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
synth_paletted_imgs = [
    ("RGB_paletted", IndirectArray(rand(1:100, 127, 257), rand(RGB{N0f8}, 100))),
    ("BRG_paletted", IndirectArray(rand(1:100, 127, 257), rand(BGR{N0f8}, 100))),
    ("ARGB_paletted", IndirectArray(rand(1:100, 127, 257), rand(ARGB{N0f8}, 100))),
    ("ABGR_paletted", IndirectArray(rand(1:100, 127, 257), rand(ABGR{N0f8}, 100))),
    ("RGBA_paletted", IndirectArray(rand(1:100, 127, 257), rand(RGBA{N0f8}, 100))),
    ("RGB_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(RGB{Float64}, 100))),
    ("BRG_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(BGR{Float64}, 100))),
    ("ARGB_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(ARGB{Float64}, 100))),
    ("ABGR_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(ABGR{Float64}, 100))),
    ("RGBA_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(RGBA{Float64}, 100))),
]

expected_img(x::Matrix{<:TransparentColor}) = convert(Array{RGBA{N0f8}}, x)
expected_img(x::Matrix{<:AbstractRGB}) = convert(Array{RGB{N0f8}}, x)

@testset "synthetic paletted images" begin
    for (case, image) in synth_paletted_imgs
        @testset "$(case)" begin
            expected = expected_img(collect(image))
            fpath = joinpath(PNG_TEST_PATH, "test_img_$(case).png")
            @testset "write" begin
                open(io->PNGFiles.save(io, image), fpath, "w")
                @test PNGFiles.save(fpath, image) == 0
            end
            @testset "read" begin
                global read_in_pngf = PNGFiles.load(fpath)
                @test typeof(read_in_pngf) <: AbstractMatrix
            end
            @testset "compare" begin
                @test all(expected .≈ read_in_pngf)
            end
            global read_in_immag = _standardize_grayness(ImageMagick.load(fpath))
            @testset "$(case): ImageMagick read type equality" begin
                # The lena image is Grayscale saved as RGB...
                @test eltype(_standardize_grayness(read_in_pngf)) == eltype(read_in_immag)
            end
            @testset "$(case): ImageMagick read values equality" begin
                imdiff_val = imdiff(read_in_pngf, read_in_immag)
                onfail(@test imdiff_val < 0.01) do
                    PNGFiles._inspect_png_read(fpath)
                    _add_debugging_entry(fpath, case, imdiff_val)
                end
            end
            path, ext = splitext(fpath)
            newpath = path * "_new" * ext
            PNGFiles.save(newpath, read_in_pngf)
            @testset "$(case): IO is idempotent" begin
                @test all(read_in_pngf .≈ PNGFiles.load(newpath))
                @test all(read_in_pngf .≈ open(io->PNGFiles.load(io), newpath))
            end
        end
    end
end


text:
"""
	legend(cmd0::String="", arg1=nothing; kwargs...)

Make legends that can be overlaid on maps. It reads specific legend-related information from input or file file.

Full option list at [`legend`]($(GMTdoc)legend.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **C** | **clearance** :: [Type => Str]

    Sets the clearance between the legend frame and the internal items [4p/4p].
    ($(GMTdoc)legend.html#c)
- **D** | **refpoint** | **pos** | **position** :: [Type => Str]  `Arg=[g|j|J|n|x]refpoint+wwidth[/height][+jjustify][+lspacing][+odx[/dy]]`

    Defines the reference point on the map for the legend using one of four coordinate systems.
    ($(GMTdoc)legend.html#d)
- **F** | **box** :: [Type => Str | Number]   `Arg=[+cclearances][+gfill][+i[[gap/]pen]][+p[pen]][+r[radius]][+s[[dx/dy/][shade]]]`

    Without further options, draws a rectangular border around the legend using *MAP_FRAME_PEN*.
    ($(GMTdoc)legend.html#f)
- $(GMT.opt_Jz)
- $(GMT.opt_P)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function legend(cmd0::String="", arg1=nothing; first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("pslegend", cmd0, arg1)

	d = KW(kwargs)
    K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd = parse_common_opts(d, cmd, [:F :c :p :t :JZ :UVXY :params], first)
	#if (!first && current_view !== nothing && !occursin("-p", cmd))  cmd *= current_view  end

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

	#cmd = add_opt(cmd, 'D', d, [:D :pos :position :refpoint],
	#			  (map=("g", nothing, 1), inside=("j", nothing, 1), paper=("x", nothing, 1),
	#			  anchor=("", arg2str, 2), width="+w", justify="+j", spacing="+l", offset="+o"))
	cmd = parse_type_anchor(d, cmd, [:D :pos :position :refpoint])
	cmd = add_opt(cmd, 'C', d, [:C :clearance])

	r = finish_PS_module(d, "pslegend " * cmd, "", K, O, true, arg1)
	gmt("destroy")
	return r
end

# ---------------------------------------------------------------------------------------------------
legend!(cmd0::String="", arg1=nothing; first=false, kw...) = legend(cmd0, arg1;first=false, kw...)
legend(arg1; first=true, kw...)   = legend("", arg1; first=first, kw...)
legend!(arg1; first=false, kw...) = legend("", arg1; first=first, kw...)

const pslegend  = legend			# Alias
const pslegend! = legend!			# Alias


text:
using Documenter, XxxxYyyy

makedocs(;
    modules=[XxxxYyyy],
    format = Documenter.HTML(
        prettyurls = false
    ),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/ollecram/XxxxYyyy.jl/blob/{commit}{path}#L{line}",
    sitename="XxxxYyyy.jl",
    authors="Marcello Vitaletti",
)

deploydocs(;
    repo="github.com/ollecram/XxxxYyyy.jl",
)


text:
using Primes: isprime
import Base.divgcd

const primetable =
    [2,3,5]
const factortable =
    [UInt8[], UInt8[1], UInt8[0,1], UInt8[2], UInt8[0,0,1]]
const factorialtable =
    [UInt32[], UInt32[], UInt32[1], UInt32[1,1], UInt32[3,1], UInt32[3,1,1]]
const bigprimetable =
    [[big(2)], [big(3)], [big(5)]]
const bigone = Ref{BigInt}(big(1))

# Make a prime iterator
struct PrimeIterator
end
primes() = PrimeIterator()

Base.IteratorSize(::Type{PrimeIterator}) = Base.IsInfinite()
Base.IteratorEltype(::Type{PrimeIterator}) = Base.HasEltype()
Base.eltype(::PrimeIterator) = Int

# Get the `n`th prime; store all primes up to the `n`th if not yet available
function prime(n::Int)
    p = last(primetable)
    while length(primetable) < n
        p = p + 2
        while !isprime(p)
            p += 2
        end
        push!(primetable, p)
        push!(bigprimetable, [big(p)])
    end
    @inbounds return primetable[n]
end

Base.iterate(::PrimeIterator, n = 1) = prime(n), n+1

# get primes and their powers as `BigInt`, also cache all results
function bigprime(n::Integer, e::Integer=1)
    e == 0 && return bigone[]
    p = prime(n) # triggers computation of prime(n) if necessary
    @inbounds l = length(bigprimetable[n])
    @inbounds while l < e
        # compute next prime power as approximate square of existing results
        k = (l+1)>>1
        push!(bigprimetable[n], bigprimetable[n][k]*bigprimetable[n][l+1-k])
        l += 1
    end
    @inbounds return bigprimetable[n][e]
end

# A custom `Integer` subtype to store an integer as its prime factorization
struct PrimeFactorization{U<:Unsigned} <: Integer
    powers::Vector{U}
    sign::Int8
end
PrimeFactorization(powers::Vector{U}) where {U<:Unsigned} =
    PrimeFactorization{U}(powers, one(Int8))

# define our own factor function, returning an instance of PrimeFactorization
function primefactor(n::Integer)
    iszero(n) && return PrimeFactorization(UInt8[], zero(Int8))
    sn = n < 0 ? -one(Int8) : one(Int8)
    n = abs(n)
    m = length(factortable)
    while m < abs(n)
        m += 1
        powers = UInt8[] # should be sufficient for all integers up to 2^255
        a = m
        for p in primes()
            f = 0
            anext, r = divrem(a, p)
            while r == 0
                f += 1
                a = anext
                anext, r = divrem(a, p)
            end
            push!(powers, f)
            a == 1 && break
        end
        push!(factortable, powers)
    end
    @inbounds return PrimeFactorization(copy(factortable[n]), sn)
end

function primefactorial(n::Integer)
    n < 0 && throw(DomainError(n))
    m = length(factorialtable)-1
    @inbounds while m < n
        prevfactorial = factorialtable[m+1]
        m += 1
        f = primefactor(m).powers
        powers = copy(prevfactorial)
        if length(f) > length(powers) # can at most be 1 larger
            push!(powers, 0)
        end
        for k = 1:length(f)
            powers[k] += f[k]
        end
        push!(factorialtable, powers)
    end
    @inbounds return PrimeFactorization(copy(factorialtable[n+1]))
end

# Methods for PrimeFactorization:
Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign)

Base.one(::Type{PrimeFactorization{U}}) where {U<:Unsigned} =
    PrimeFactorization(Vector{U}())
Base.zero(::Type{PrimeFactorization{U}}) where {U<:Unsigned} =
    PrimeFactorization(Vector{U}(), zero(Int8))

Base.promote_rule(P::Type{<:PrimeFactorization},::Type{<:Integer}) = P
Base.promote_rule(P::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt

Base.convert(P::Type{<:PrimeFactorization}, n::Integer) = convert(P, primefactor(n))
function Base.convert(::Type{BigInt}, a::PrimeFactorization)
    A = one(BigInt)
    for (n, e) in enumerate(a.powers)
        if !iszero(e)
            MPZ.mul!(A, bigprime(n, e))
        end
    end
    return a.sign < 0 ? MPZ.neg!(A) : A
end
Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization{U}) where {U<:Unsigned} =
    a
Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization) where {U<:Unsigned} =
    PrimeFactorization(convert(Vector{U}, a.powers), a.sign)

Base.:(==)(a::PrimeFactorization, b::PrimeFactorization) =
    a.powers == b.powers && a.sign == b.sign
function Base.:<(a::PrimeFactorization, b::PrimeFactorization)
    if a.sign != b.sign
        return a.sign < b.sign
    elseif a.sign < 0
        return <(-b, -a)
    else
        ag, bg = divgcd(a, b)
        if length(ag.powers) <= length(bg.powers) &&
                all(k->ag.powers[k]<bg.powers[k], 1:length(ag.powers))
            return true
        else
            return convert(BigInt, ag) < convert(BigInt, bg)
        end
    end
end

# Methods for PrimeFactorization: only fast multiplication, and lcm and gcd.
# Addition and subtraction will require conversion to BigInt
Base.sign(a::PrimeFactorization) = a.sign

Base.:-(a::PrimeFactorization) = PrimeFactorization(a.powers, -a.sign)
function Base.:*(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T}
    if a.sign == 0
        return a
    elseif b.sign ==0
        return b
    else
        return PrimeFactorization(_vadd!(copy(a.powers), b.powers), a.sign*b.sign)
    end
end
function Base.gcd(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T}
    if a.sign == 0
        return b
    elseif b.sign ==0
        return a
    else
        return PrimeFactorization(_vmin!(copy(a.powers), b.powers))
    end
end
function Base.lcm(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T}
    if a.sign == 0
        return a
    elseif b.sign ==0
        return b
    else
        return PrimeFactorization(_vmax!(copy(a.powers), b.powers))
    end
end
Base.divgcd(a::PrimeFactorization, b::PrimeFactorization) = divgcd!(copy(a), copy(b))
function divgcd!(a::PrimeFactorization, b::PrimeFactorization)
    af, bf = a.powers, b.powers
    for k = 1:min(length(af), length(bf))
        gk = min(af[k], bf[k])
        af[k] -= gk
        bf[k] -= gk
    end
    while length(af) > 0 && iszero(last(af))
        pop!(af)
    end
    while length(bf) > 0 && iszero(last(bf))
        pop!(bf)
    end
    return a, b
end

# split `a::PrimeFactorization` into a square `s` and a remainder `r`, such that
# `a = s^2 * r` and the powers in the prime factorization of `r` are zero or one
function splitsquare(a::PrimeFactorization)
    r = PrimeFactorization(map(p->convert(UInt8, isodd(p)), a.powers), a.sign)
    while length(r.powers) > 0 && iszero(last(r.powers))
        pop!(r.powers)
    end
    s = PrimeFactorization(map(p->(p>>1), a.powers))
    while length(s.powers) > 0 && iszero(last(s.powers))
        pop!(s.powers)
    end
    return s, r
end

# given a list of numerators and denominators, compute the common denominator and
# the rescaled numerator after putting all fractions at the same common denominator
function commondenominator!(nums::Vector{P}, dens::Vector{P}) where {P<:PrimeFactorization}
    isempty(nums) && return one(P)
    # accumulate lcm of denominator
    den = PrimeFactorization(copy(dens[1].powers))
    for i = 2:length(dens)
        _vmax!(den.powers, dens[i].powers)
    end
    # rescale numerators
    for i = 1:length(nums)
        _vsub!(_vadd!(nums[i].powers, den.powers), dens[i].powers)
    end
    return den
end

# auxiliary function to compute sums of a list of PrimeFactorizations as quickly as possible
function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list))
    # first compute gcd to take out common factors
    g = PrimeFactorization(copy(list[ind[1]].powers))
    for k in ind
        _vmin!(g.powers, list[k].powers)
    end
    for k in ind
        _vsub!(list[k].powers, g.powers)
    end
    L = length(ind)
    if L > 32
        l = L >> 1
        s = sumlist!(list, first(ind).+(0:l-1)) + sumlist!(list, first(ind).+(l:L-1))
    else
        # do sum
        s = big(0)
        for k in ind
            MPZ.add!(s, convert(BigInt, list[k]))
        end
    end
    return MPZ.mul!(s, convert(BigInt, g))
end

# Mutating vector methods that also grow and shrink as required
function _vmin!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned}
    while length(af) > length(bf)
        pop!(af)
    end
    @inbounds for k = 1:length(af)
        af[k] = min(af[k], bf[k])
    end
    while length(af) > 0 && iszero(last(af))
        pop!(af)
    end
    return af
end
function _vmax!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned}
    while length(bf) > length(af)
        push!(af, zero(U))
    end
    @inbounds for k = 1:length(bf)
        af[k] = max(af[k], bf[k])
    end
    return af
end
function _vadd!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned}
    while length(bf) > length(af)
        push!(af, zero(U))
    end
    @inbounds for k = 1:length(bf)
        af[k] = +(af[k], bf[k])
    end
    return af
end
function _vsub!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned}
    if length(bf) > length(af)
        throw(OverflowError())
    end
    @inbounds for k = 1:length(bf)
        bf[k] > af[k] && throw(OverflowError())
        af[k] -= bf[k]
    end
    while length(af) > 0 && iszero(last(af))
        pop!(af)
    end
    return af
end


text:
"""
This struct contains the REC market parameters
    for price forcasts using CEM and endogeneous Economic Dispatch.
The REC market assumes a percentage RPS requirment of annual energy consumption
which is parameterized within this struct.
"""

struct RECMarket

    rec_req::Float64    # $/MW/investment period
    price_cap::Float64        # MW

    function RECMarket(req::Number, price_cap::Number)
        @assert 0 <= req <= 1
        @assert price_cap >= 0
        new(req, price_cap)
    end

end


text:
@testset "ad.jl" begin
    @testset "logp" begin
        # Hand-written log probabilities for vector `x = [s, m]`.
        function logp_gdemo_default(x)
            s = x[1]
            m = x[2]
            dist = Normal(m, sqrt(s))

            return logpdf(InverseGamma(2, 3), s) +
                   logpdf(Normal(0, sqrt(s)), m) +
                   logpdf(dist, 1.5) +
                   logpdf(dist, 2.0)
        end

        test_model_ad(gdemo_default, logp_gdemo_default)

        @model function wishart_ad()
            return v ~ Wishart(7, [1 0.5; 0.5 1])
        end

        # Hand-written log probabilities for `x = [v]`.
        function logp_wishart_ad(x)
            dist = Wishart(7, [1 0.5; 0.5 1])
            return logpdf(dist, reshape(x, 2, 2))
        end

        test_model_ad(wishart_ad(), logp_wishart_ad)
    end

    # https://github.com/TuringLang/Turing.jl/issues/1595
    @testset "dot_observe" begin
        function f_dot_observe(x)
            logp, _ = DynamicPPL.dot_observe(
                SampleFromPrior(), [Normal(), Normal(-1.0, 2.0)], x, VarInfo()
            )
            return logp
        end
        function f_dot_observe_manual(x)
            return logpdf(Normal(), x[1]) + logpdf(Normal(-1.0, 2.0), x[2])
        end

        # Manual computation of the gradient.
        x = randn(2)
        val = f_dot_observe_manual(x)
        grad = ForwardDiff.gradient(f_dot_observe_manual, x)

        @test ForwardDiff.gradient(f_dot_observe, x) ≈ grad

        y, back = Tracker.forward(f_dot_observe, x)
        @test Tracker.data(y) ≈ val
        @test Tracker.data(back(1)[1]) ≈ grad

        y, back = Zygote.pullback(f_dot_observe, x)
        @test y ≈ val
        @test back(1)[1] ≈ grad
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
using RCall
using DataFrames, DataFramesMeta
using CSV

max_stars_repo_name: ppalmes/AMLPipelineR

max_issues_repo_name: ppalmes/AMLPipelineR

max_forks_repo_name: ppalmes/AMLPipelineR


text:
using Test
using HKQM
using Logging
using ForwardDiff

disable_logging(Logging.Info)

@testset "Elements" begin
    ca = CubicElementArray(5u"Å", 3)
    @test size(ca) == (3,3,3)
    @test get_center(ca[2,2,2]) ≈ zeros(3).*u"bohr"

    ceg = CubicElementGrid(5u"Å", 4, 16)
    @test size(ceg) == (16,16,16,4,4,4)

    ev = ElementVector(0, 2, 5, 7)
    @test length(ev) == 3

    egv = ElementGridVector(ev, 6)
    @test size(egv) == (6,3)
    @test egv[begin, begin] >= 0
    @test egv[end, end] <= 7
end

@testset "Tensors" begin
    ceg = CubicElementGrid(5u"Å", 4, 16)
    ct = CoulombTransformation(ceg, 16)
    clog = CoulombTransformationLog(ceg, 16)
    ca = CoulombTransformationLocal(ceg, 16)
    cll  = CoulombTransformationLogLocal(ceg, 16)
    @test size(ct) == (16,16,4,4,16)
    @test size(ct) == size(clog)
    @test size(ct) == size(ca)
    @test size(ct) == size(cll)

    # Symmetry
    @test ct.w[1]*ct[1,3,2,4,1] ≈ ct.w[3]*ct[3,1,4,2,1]
    @test ct.w[1]*clog[1,3,2,4,2] ≈ ct.w[3]*clog[3,1,4,2,2]
    @test ct.w[1]*ca[1,3,2,4,3] ≈ ct.w[3]*ca[3,1,4,2,3]
    @test ct.w[1]*cll[1,3,2,4,4] ≈ ct.w[3]*cll[3,1,4,2,4]
end

@testset "Operators" begin
    ceg = CubicElementGrid(5u"Å", 4, 16)
    r = position_operator(ceg)
    p = momentum_operator(ceg)
    x = r[1]
    @test unit(x) == u"bohr"
    @test x == r.operators[1]
    @test x != r.operators[2]
    @test x != r.operators[3]
    @test size(x) == size(r)
    @test_throws AssertionError r + p
    @test sin(1u"bohr^-1"*x).vals == sin.(x.vals)
    @test cos(1u"bohr^-1"*x).vals == cos.(x.vals)
    @test exp(1u"bohr^-1"*x).vals == exp.(x.vals)
    @test (x + x).vals == 2.0 .* x.vals
    @test (2x).vals ==  (x + x).vals
    @test (x^2).vals == x.vals.^2
    @test sqrt(x^2).vals ≈ abs.(x.vals)
    @test (x^2).vals == (x*x).vals
    @test unit(x^2) == u"bohr^2"
    @test unit(x+x) == u"bohr"
    @test unit(sqrt(x^2)) == u"bohr"
end

@testset "Quantum States" begin
    ceg = CubicElementGrid(5u"Å", 4, 16)
    r = position_operator(ceg)
    r² = r⋅r
    gv = exp(-1u"bohr^-2"*r²)
    ψ = QuantumState(ceg, gv.vals)
    normalize!(ψ)
    @test bracket(ψ,ψ) ≈ 1
    ϕ = ψ + 2ψ
    @test bracket(ϕ,ψ) ≈ 3
    rr = bracket(ψ, r, ψ)
    @test sqrt( sum(x->x^2, rr) ) < 1u"bohr"*1E-12
end

@testset "Nuclear potential" begin
    tn = test_nuclear_potential(5u"Å", 4, 64, 64; mode="preset")
    @test abs( (tn["integral"] - tn["total reference"]) / tn["total reference"]) < 1e-6
end

@testset "Poisson equation" begin
    ec = test_accuracy(5u"Å",4,24,48; showprogress=false)  # Poor accuracy
    eref = gaussian_coulomb_integral()[1]*u"hartree"
    e = ec["calculated"]
    @test abs((e-eref)/eref) < 1e-3  # Calculation had poor accuracy
end

@testset "Forward mode AD" begin
    d = test_accuracy_ad(5u"Å", 4, 16, 96; α1=1.3, α2=0.8, d=0.5, showprogress=false)
    # Grid resolution is small so not very accurate
    @test all(abs.(d[1] .- d[2])./d[2] .< 1e-2)

    # Helmholtz equation derivative
    function f(k)
        ceg = CubicElementGrid(5u"Å",2,16)
        ct=optimal_coulomb_tranformation(ceg, 4; k=k)
        ρ = HKQM.density_tensor(ceg)
        V = poisson_equation(ρ, ct)
        return HKQM.integrate(ρ, ceg, V)
    end
    g(x) = ForwardDiff.derivative(f, x)
    g(1.)
end

@testset "Derivative and kinetic energy" begin
    a = test_kinetic_energy(5u"Å", 4, 32; ν=1, ω=3)
    @test abs(a[1]-a[2]) / a[2] < 1e-10
end


text:
module sssudoko
using Base.Dates

row1 = [0,0,0,0,0,0,0,7,0]
row2 = [0,4,0,2,9,1,8,6,5]
row3 = [8,5,0,0,0,0,0,0,9]
row4 = [4,0,8,6,0,0,2,0,0]
row5 = [3,2,0,0,4,0,0,8,6]
row6 = [0,0,1,0,0,9,4,0,3]
row7 = [1,0,0,0,0,0,0,9,8]
row8 = [5,8,4,9,6,3,0,2,0]
row9 = [0,6,0,0,0,0,0,0,0]
p1=[0,7,0,9,6,4,0,0,0]
p2=[0,3,0,1,0,0,0,4,6]
p3=[0,0,0,8,0,0,0,0,5]
p4=[0,2,0,0,0,0,0,0,3]
p5=[6,0,0,0,0,0,0,0,4]
p6=[4,0,0,0,0,0,0,9,0]
p7=[1,0,0,0,0,3,0,0,0]
p8=[3,9,0,0,0,5,0,1,0]
p9=[0,0,0,7,1,6,0,5,0]
#array = hvcat (9,row1,row2,row3,row4,row5,row6,row7,row8,row9)
array =hvcat(9,p1,p2,p3,p4,p5,p6,p7,p8,p9)
result = []
type CandidateSet{T}
    rx::Int
    ry::Int
    candidates::Array{T}
end

function howManySquares(matrix::Array)
    count = 0
    for thing in matrix
        if thing==0 count+=1
        end
    end
    return count
end

"Returns an instance of CandidateSet"
function generateCandidates (rx::Int,ry::Int, matrix::Array)
    starter=CandidateSet(rx,ry,Int[0])
    if matrix [rx,ry]!=0
        return starter #because we don't need a candidate
	#not good code - need to think of a way to just not have this value
    else
    allMembers = rowColumnBlock(rx,ry,matrix)
   # print (allMembers)
    elements=Int[]
    for i in 1:9
        if !in(i,allMembers)
            push!(elements,i)
        end
    end
    if (elements==[])
       # print ("rx=",rx,"ry=",ry)
        throw (UndefRefError()) # because this is an matrix with a contradiction not yet discovered
        end
          candidate=CandidateSet(rx,ry,elements)
    return candidate
    end
end

function createCandidateMatrix (things::Array)
    matrix = CandidateSet[]
    for n in 1:9
        for m in 1:9
            candidates = generateCandidates (m,n,things)
                push!(matrix,candidates)
            end
    end
    matrix2 = reshape (matrix,9,9)
   #@show(matrix2)
    return (matrix2)
end

"extract relevant elements for considering a square in sudoko"
"the row the item is in, the column it is in and the 3*3 local that its in"
function rowColumnBlock (rx::Int, ry::Int, matrix::Array)
   # matrix[rx,ry] = 0
  row = matrix[rx,:]
    column = matrix[:, ry]
    x1 = div(rx-1,3)*3
    y1 = div(ry-1,3)*3
    x2=x1+3
    y2=y1+3
    x1+=1
    y1+=1
   # @show x1,x2,y1,y2
    block = matrix[x1:x2,y1:y2]
    #=we've got the row, column and block so now
    get membership and see if there's a missing one =#
    allMembers = [reshape(block,1,9);reshape(column,1,9);reshape(row,1,9)]
    return allMembers
end
#=
function exclusiveRowColumnBlock(rx::Int, ry::Int, matrix::Array)
    value = matrix[rx,ry]
    matrix[rx,ry] = 0
    retval = rowColumnBlock(rx,ry,matrix)
    matrix[rx,ry] = value
    return retval
end
=#

function copytx (matrix::Array{Int})
    retval = Array (Int,length(matrix))
    for i in 1:length (matrix)
        retval[i] = matrix[i]
        end
        return reshape (retval,9,9)
        end

function notConsistent (matrix::Array)
    copyt = copytx(matrix)
    for n in 1:9
        for m in 1:9
        if (matrix[n,m]!=0)
          copyt[n,m]=10
          row = copyt[n,:]
          if (in(matrix[n,m],row))
            #print ("not consistent")
            return (true)
            end
          column = copyt[:, m]
          if (in(matrix[n,m],column))
            #print ("not consistent")
            return (true)
          end
         x1 = div(n-1,3)*3
        y1 = div(m-1,3)*3
        x2=x1+3
        y2=y1+3
        x1+=1
        y1+=1
   # @show x1,x2,y1,y2
        block = copyt[x1:x2,y1:y2]
        if (in(matrix[n,m],block))
            #print ("not consistent")
            return (true)
        end #if
        copyt [n,m] = matrix[n,m]
        end #for
        end #for
    end #if
    return false
    end

function solutionSearch( matrix::Array, depth::Int)
    # for every element in matrix which is a 0 (undefined) look
    # at all the other 0s to find any that contain exclusions
    # which can solve this square
    # check singles and pairs only for now !
    try
    touched = false
    newMatrix = Int[]
    candidateSet = createCandidateMatrix(matrix)
    if (notConsistent(matrix))
    return false
    else
    for i in 1:81
        if (matrix[i]==0 ) #something to do
           if (length(candidateSet[i].candidates)==1)
              #simple write in
        #      println("write in=",candidateSet[i].candidates[1]," at ",i)
              matrix[i]=candidateSet[i].candidates[1]
            if (notConsistent(matrix))
         #         println("found contradiction")
                  return false
              end #unwind as this is a dead attempt
              touched = true
            end #end if
        end #end if
    end #end for
    candidateSet=createCandidateMatrix(matrix) #in case a repair was needed
    tries = zeros(9,9)
    marktries(tries,matrix)
    while (!complete (tries))
    index = shortest(candidateSet,tries)
    println("depth=",depth," index=",index)
    for guess in candidateSet[index].candidates
              #= for every option try writing it in and then solving the matrix with that assumption
              if a contradiction arises then stop and unwind and try a different option =#
              newMatrix = copytx(matrix)
              newMatrix[index] =guess
              if(!notConsistent (newMatrix))
                  if(solutionSearch (newMatrix,depth+1))
                    touched = true
                    matrix = newMatrix
                   end
              else
               matrix[index] = 0
              end#if
        end #for
    return false
    end #else
    end
     if(rem(datetime2unix(now()),10)==0.0)
       println(matrix)
       end
      if (howManySquares(matrix)==0)
       throw (InterruptException())
       end

 catch (e)
        if (isa(e,InterruptException))
            global result = matrix
            return true
         end
  return false
  end
 end

function marktries(tries::Array,tester::Array)
    for i in 1:length(tester)
        if (tester[i]!=0)
            tries[i]=1
        end#if
    end#for
end

 function complete (tries::Array)
    if (in(0,tries))
        return false
    else
    return true
    end
end

 function shortest (candidates::Array{CandidateSet}, tries::Array)
    theOne =0
    for i in 1:length(candidates)
       if (tries[i]==0.0)
       if candidates[i].candidates[1]!=0 # there is a list
       if (theOne==0)
        theOne=i else
        if (length(candidates[i].candidates) < length(candidates[theOne].candidates))
            theOne = i
        end #if
        end #else
        end #if
    end #if
    end #for
    return theOne
end

cache= []
cache2 = []

 function cleanArray(thing::Int,target::Array)
    retVal =Int[]
    if (length(target)==1) return [0] end
  #  print ("starting>>", target)
    for i = 1:length(target)
        if target[i]!=thing
            push!(retVal,target[i])
        end
    end
   #         print ("done >> ",retVal)
    return retVal
    end

function hypothesis (guess::Int, tester::Array{CandidateSet})
  #print (tester)
  for i in tester
    if in (guess,i.candidates)
    return (false)
    end
    end
    return (true)
end

@time(solutionSearch(array,0))

println("result is ",result)

end #module


text:
# In sequential case one can drive the system with for loop:

while !(learner.loss() < step)
    xi = ask!(learner,true)
    yi = f(xi)
    tell!(learner,(xi,yi))
end

# However when paralelism is taken into account it gets a bit more complex

tasks = Channel()
results = Channel()

unresolved = 0
while true

    if !(learner.loss() < step)

        xi = ask!(learner,true)
        put!(tasks,xi)
        unresolved += 1

        if unresolved < N
            continue
        end
    end

    if unresolved == 0
        break
    end

    yi = take!(results)
    tell!(learner,xi)
    unresolved -= 1

end

using Adaptive
using TaskMaster

@everywhere f(p) = exp(-p[1]^2 - p[2]^2)

master = WorkMaster(f)
learner2d = AdaptiveLearner2D([(-3,+3),(-3,+3)])
loop = Loop(master,learner2d)
evaluate!(loop,learner->learner.loss()<0.05)