text:
using Reexport

@reexport using DataFrames
using MLStyle
using Distributions

export sourceXform

function sourceXform(m::Model)
    m = canonical(m)
    pars = parameters(m)
    @gensym t result

    proc(m, st::Let)        = :($(st.x) = $(st.rhs))
    proc(m, st::Return)     = nothing
    proc(m, st::LineNumber) = nothing

    function proc(m, st::Follows)
        if st.x ∈ pars
            return @q begin
                $(st.x) = rand($(st.rhs))
                $t = xform($(st.rhs))

                $result = merge($result, ($(st.x)=$t,))
            end
        else
            return @q begin
                $(st.x) = rand($(st.rhs))
            end
        end
    end

    body = buildSource(m, proc) |> striplines

    argsExpr = Expr(:tuple,arguments(m)...)

    @gensym rand

    flatten(@q (
        function $rand(args...;kwargs...)
            @unpack $argsExpr = kwargs
            $result = NamedTuple()
            $body
            as($result)
        end
    ))

end

export makeXform
function makeXform(m :: Model)
    fpre = @eval $(sourceXform(m))
    f(;kwargs...) = Base.invokelatest(fpre; kwargs...)
end

export xform
xform(m::Model; kwargs...) = makeXform(m)(;kwargs...)

function xform(d)
    if hasmethod(support, (typeof(d),))
        return asTransform(support(d))
    end
end

function asTransform(supp:: RealInterval)
    (lb, ub) = (supp.lb, supp.ub)

    (lb, ub) == (-Inf, Inf) && (return asℝ)
    (lb, ub) == (0.0,  Inf) && (return asℝ₊)
    (lb, ub) == (0.0,  1.0) && (return as𝕀)
    error("asTransform($supp) not yet supported")
end

# export xform
# xform(::Normal)       = asℝ
# xform(::Cauchy)       = asℝ
# xform(::Flat)         = asℝ

# xform(::HalfCauchy)   = asℝ₊
# xform(::HalfNormal)   = asℝ₊
# xform(::HalfFlat)     = asℝ₊
# xform(::InverseGamma) = asℝ₊
# xform(::Gamma)        = asℝ₊
# xform(::Exponential)  = asℝ₊

# xform(::Beta)         = as𝕀
# xform(::Uniform)      = as𝕀

function xform(d::For)
    # allequal(d.f.(d.θs)) &&
    return as(Array, xform(d.f(d.θs[1])), size(d.θs)...)

    # TODO: Implement case of unequal supports
    @error "xform: Unequal supports not yet supported"
end

function xform(d::iid)
    as(Array, xform(d.dist), d.size...)
end

ext: jl


text:
#======================================================================================================================
Reach-Avoid Specification
======================================================================================================================#
dfa_states = collect(1:3)
dfa_props = ["a", "b"]
dfa_transitions = [(1, "!a∧!b", 1),
                    (1, "a∧!b", 2),
                    (2, "true", 2),
                    (1, "!a∧b", 3),
                    (3, "true", 3)]
dfa_accepting_state = 2
dfa_sink_state = 3
dfa_initial_state = 1


text:
# # Continuous SIMP example
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`csimp.ipynb`](@__NBVIEWER_ROOT_URL__/examples/csimp.ipynb)
#-
# ## Commented Program
#
# What follows is a program spliced with comments.
#md # The full program, without comments, can be found in the next [section](@ref csimp-plain-program).

using TopOpt

# ### Define the problem
E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force

problems = Any[
            PointLoadCantilever(Val{:Linear}, (60, 20, 20), (1.0, 1.0, 1.0), E, v, f),
            PointLoadCantilever(Val{:Linear}, (160, 40), (1.0, 1.0), E, v, f),
            HalfMBB(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
            LBeam(Val{:Linear}, Float64, force = f),
            TieBeam(Val{:Quadratic}, Float64)
            ]
problem_names = [
    "3d cantilever beam",
    "cantilever beam",
    "half MBB beam",
    "L-beam",
    "tie-beam",
    ]

i = 2
println(problem_names[i])
problem = problems[i]

# ### Parameter settings
V = 0.5 # volume fraction
xmin = 0.001 # minimum density
rmin = 3.0
steps = 40 # maximum number of penalty steps, delta_p0 = 0.1
reuse = true # adaptive penalty flag

convcriteria = Nonconvex.GenericCriteria()
penalty = TopOpt.PowerPenalty(1.0)
pcont = Continuation(penalty, steps = steps, xmin = xmin, pmax = 5.0)

# NOTE: non-convexity + computational error lead to different solutions that satisfy the KKT tolerance
mma_options = options = Nonconvex.MMAOptions(maxiter=1000)
maxtol = 0.01 # maximum tolerance
mintol = 0.0001 # minimum tolerance
b = log(mintol / maxtol) / steps
a = maxtol / exp(b)
mma_options_gen = TopOpt.MMAOptionsGen(
    steps = steps,
    initial_options = mma_options,
    ftol_gen = ExponentialContinuation(a, b, 0.0, steps + 1, mintol),
)
csimp_options = TopOpt.CSIMPOptions(
    steps = steps,
    options_gen = mma_options_gen,
    p_gen = pcont,
    reuse = reuse,
)

# ### Define a finite element solver
solver = FEASolver(Displacement, Direct, problem, xmin = xmin, penalty = penalty)

# ### Define compliance objective
comp = Compliance(problem, solver)
filter = if problem isa TopOptProblems.TieBeam
    identity
else
    DensityFilter(solver, rmin = rmin)
end
obj = x -> comp(filter(x))

# ### Define volume constraint
volfrac = TopOpt.Volume(problem, solver)
constr = x -> volfrac(filter(x)) - V

# ### Define subproblem optimizer
x0 = fill(V, length(solver.vars))
optimizer = Optimizer(
    obj, constr, x0, Nonconvex.MMA87(),
    options = mma_options, convcriteria = convcriteria,
)

# ### Define continuation SIMP optimizer
simp = SIMP(optimizer, solver, penalty.p)
cont_simp = ContinuationSIMP(simp, steps, csimp_options)

# ### Solve
result = cont_simp(x0)

@show result.convstate
@show optimizer.workspace.iter
@show result.objval

# ### (Optional) Visualize the result using Makie.jl
# Need to run `using Pkg; Pkg.add(Makie)` first
# ```julia
# using TopOpt.TopOptProblems.Visualization: visualize
# fig = visualize(problem; topology = result.topology,
#     problem; topology = result.topology, default_exagg_scale = 0.07,
#     scale_range = 10.0, vector_linewidth = 3, vector_arrowsize = 0.5,
# )
# Makie.display(fig)
# ```

#md # ## [Plain Program](@id csimp-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [csimp.jl](csimp.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
using
  Documenter,
  Literate,
  Plots,   # to not capture precompilation output
  PassiveTracerFlows

# Gotta set this environment variable when using the GR run-time on Travis CI.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")

examples = [
    "cellularflow.jl",
]

for example in examples
  withenv("GITHUB_REPOSITORY" => "FourierFlows/PassiveTracerFlowsDocumentation") do
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
    Literate.notebook(example_filepath, OUTPUT_DIR)
    Literate.script(example_filepath, OUTPUT_DIR)
  end
end

#####
##### Build and deploy docs
#####

# Set up a timer to print a space ' ' every 240 seconds. This is to avoid Travis CI
# timing out when building demanding Literate.jl examples.
Timer(t -> println(" "), 0, interval=240)

format = Documenter.HTML(
  collapselevel = 2,
     prettyurls = get(ENV, "CI", nothing) == "true",
      canonical = "https://fourierflows.github.io/PassiveTracerFlowsDocumentation/stable/"
)

makedocs(
   modules = [PassiveTracerFlows],
   doctest = false,
     clean = true,
 checkdocs = :all,
    format = format,
   authors = "Gregory L. Wagner and Navid C. Constantinou",
  sitename = "PassiveTracerFlows.jl",
     pages = Any[
           "Home" => "index.md",
           "Examples" => Any[
             "TracerAdvectionDiffusion" => Any[
               "literated/cellularflow.md",
               ]
           ],
           "Modules" => Any[
             "modules/traceradvectiondiffusion.md",
           ],
           "DocStrings" => Any[
             "man/types.md",
             "man/functions.md"
             ]
           ]
)

withenv("GITHUB_REPOSITORY" => "FourierFlows/PassiveTracerFlowsDocumentation") do
  deploydocs(       repo = "github.com/FourierFlows/PassiveTracerFlowsDocumentation.git",
                versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
            push_preview = false,
               devbranch = "main"
            )
end

max_stars_repo_name: FourierFlows/PassiveTracerFlows

max_issues_repo_name: FourierFlows/PassiveTracerFlows

max_forks_repo_name: FourierFlows/PassiveTracerFlows


text:
type Neighbours6{T} <: FieldVector{6, T}
  left::Union{T,Void}
  right::Union{T,Void}
  front::Union{T,Void}
  back::Union{T,Void}
  top::Union{T,Void}
  bottom::Union{T,Void}

  Neighbours6{T}() where T = new(nothing,nothing,nothing,nothing,nothing,nothing)
end

abstract type INode end
abstract type IBlock end
abstract type IChunk end
#RChunk = Union{Void,IChunk}

const NodeNeighbours6 = Neighbours6{INode}
const BlockNeighbours = Neighbours6{IBlock}

type Node6{T} <: INode
  value::Union{T,Void}
  next::NodeNeighbours6

  Node6{T}() where T = new(nothing,NodeNeighbours6())
  Node6{T}(value::T) where T = new(value,NodeNeighbours6())
end

type Block <: IBlock
  typ::Int32
  pos::Vec3f
  flags::UInt32
  active::Bool
  visible::Bool
  surrounded::Bool
  sides::Array{UInt32,1}

  parent::IChunk
  next::BlockNeighbours

  Block(parent::IChunk, pos=Vec3f,typ=0) = new(typ,pos,0,true,true,false,resetSides(),parent,BlockNeighbours()) #,zeros(Mat4x4f),zeros(Mat4x4f))
end

#Base.isless{T}(a::Ref{T}, b::Ref{T}) = a.x < b.x
Base.isless(a::Block, b::Block) = a.pos < b.pos

type Chunk <: IChunk
  active::Bool
  len::UInt32
  childs::Array{Block,3}
  filtered::Array{Block,1}
  count::UInt32
  fileredCount::UInt32

  function Chunk(len::Integer)
    this = new(true,len,Array{Block,3}(len,len,len),Block[],0,0)
    createBlocks(this)
    linkBlockNeighbours(this)
    this
  end
end

function clean(this::Union{Void,Chunk})
  if this == nothing return end
  this.childs = Array{Block,3}(0,0,0)
  this.filtered = Block[]
  gc() # force garbage collection, free memory
end

isType(this::Block, typ) = this.typ == typ

isValid(this::Chunk) = this.count > 0 && this.fileredCount > 0

isSeen(this::Block) = isActive(this) && isVisible(this) && this.typ > 0

i=0
i+=1; const LEFT_SIDE = i #0x1
i+=1; const RIGHT_SIDE = i #0x2
i+=1; const TOP_SIDE = i #0x4
i+=1; const BOTTOM_SIDE = i #0x8
i+=1; const FRONT_SIDE = i #0x10
i+=1; const BACK_SIDE = i #0x20
const MAX_SIDES = i

resetSides() = fill(UInt32(1),MAX_SIDES)
resetSides(this::Block) = this.sides=resetSides()

function createBlocks(this::Chunk)
  len = this.len
  for x=1:len; for y=1:len; for z=1:len; this.childs[x,y,z] = Block(this,Vec3f(x,y,z)); end; end; end
end

function linkBlockNeighbours(this::Chunk)
  a = this.childs
  len = this.len
  for x=1:len; for y=1:len; for z=1:len;
    pos = (x,y,z)
    b = getBlock(this,pos)
    b.next.left = getBlock(this,pos,LEFT_SIDE)
    b.next.right = getBlock(this,pos,RIGHT_SIDE)
    b.next.bottom = getBlock(this,pos,BOTTOM_SIDE)
    b.next.top = getBlock(this,pos,TOP_SIDE)
    b.next.back = getBlock(this,pos,BACK_SIDE)
    b.next.front = getBlock(this,pos,FRONT_SIDE)
  end; end; end
end

function getBlock(this::Chunk, pos::Tuple{Integer,Integer,Integer}, side=0)
  a = this.childs
  len = this.len
  (x,y,z) = pos
  b = nothing
  if side == LEFT_SIDE
    b=x>1?a[x-1,y,z]:nothing
  elseif side == RIGHT_SIDE
    b=x<len?a[x+1,y,z]:nothing
  elseif side == BOTTOM_SIDE
    b=y>1?a[x,y-1,z]:nothing
  elseif side == TOP_SIDE
    b=y<len?a[x,y+1,z]:nothing
  elseif side == BACK_SIDE
    b=z>1?a[x,y,z-1]:nothing
  elseif side == FRONT_SIDE
    b=z<len?a[x,y,z+1]:nothing
  else
    b=x>=1&&x<=len&&y>=1&&y<=len&&z>=1&&z<=len?a[x,y,z]:nothing
  end
  b
end

function hideUnseen(this::Chunk)
  for b in this.childs
    if !isSeen(b) continue end # skip invisible

    blocked=0

    if (next=b.next.left) != nothing
      r=b.sides[LEFT_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.right) != nothing
      r=b.sides[RIGHT_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.bottom) != nothing
      r=b.sides[BOTTOM_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.top) != nothing
      r=b.sides[TOP_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.back) != nothing
      r=b.sides[BACK_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.front) != nothing
      r=b.sides[FRONT_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    setSurrounded(b,blocked >= 6)
  end
end

function setFlag(this::Block, flag::Unsigned, add::Bool)
  if add this.flags |= flag
  else this.flags ⊻= flag #⊻ = xor and is not shown on windows...
  end
end

isActive(this::Block) = this.active
isVisible(this::Block) = this.visible
isSurrounded(this::Block) = this.surrounded

isValid(this::Block) = isActive(this) && isVisible(this) && !isSurrounded(this) && this.typ > 0

setActive(this::Block, active::Bool) = this.active = active
setVisible(this::Block, visible::Bool) = this.visible = visible
setSurrounded(this::Block, surrounded::Bool) = this.surrounded = surrounded

hideType(this::Chunk, typ::Integer) = for b in this.childs; if b.typ == typ; setVisible(b,false); end; end
removeType(this::Chunk, typ::Integer) = for b in this.childs; if b.typ == typ; setActive(b,false); end; end

showAll(this::Chunk) = for b in this.childs
  setVisible(b,true)
  setSurrounded(b,false)
  resetSides(b)
end

const FrustumNode = Node6{Tuple{Block,Any}}

printList(list,n=1) = begin
  n = round(n)
  if n<1 error("n < 1 is invalid!") end
  println(typeof(list)," with ",length(list)," entries shows each ",n,n==1?"st":(n==2?"nd":(n==3?"rd":"th"))," element:")
  i=0; for (k,v) in list; i+=1; if i == n i=0; println(" ",k) end; end
  println()

  #a = collect(keys(list))[1:X:end]
  #Base.showarray(STDOUT,a,false)
  #show(STDOUT, "text/plain", a)
  #display(a)
  #whos()
end

function checkInFrustum(this::Chunk, fstm::Frustum)
  pos = getNearPosition(fstm)
  list = SortedDict()

  open("frustum.txt", "w") do f

  for b in this.childs
    result = checkSphere(fstm, b.pos, 1.5)
    clist = result[2]
    visible = result[1] != :FRUSTUM_OUTSIDE
    setVisible(b,visible)

    # sort
    if visible
      (x,y,z) = (clist[:X],clist[:Y],clist[:Z])
      write(f, string("[",x,", ",y,", ",z,"]")*"\n")

      rx=round(x)
      ry=round(y)
      rz=round(z)

      if rx % 2 != 0
        rx+=1;
      end

      k = (rx,ry)

      if (abs(x) <= 1 && abs(y) <= 1 && abs(z) <= 1) || z < 0
        setVisible(b,false)
      else
        kv=z=>b
        if haskey(list,k)

          p = first(list[k])
          list[k][z] = b
          q = first(list[k])

          if p[1] == q[1] kv=p[1]=>p[2] else b=p[2] end
          setVisible(b,false)
        end
        list[k] = SortedDict(kv)
      end
    end

    #=
    if visible && false
      current = FrustumNode((b,clist))
      parent = nothing
      if parent == nothing parent = current
      else
        nlist = parent.value[2]
        while true
          lt = clist[:FRUSTUM_LEFT] > clist[:FRUSTUM_TOP]
          l = nlist[:FRUSTUM_LEFT] > clist[:FRUSTUM_LEFT]
          r = nlist[:FRUSTUM_RIGHT] < clist[:FRUSTUM_RIGHT]
          t = nlist[:FRUSTUM_TOP] > clist[:FRUSTUM_TOP]
          b = nlist[:FRUSTUM_BOTTOM] < clist[:FRUSTUM_BOTTOM]
          n = nlist[:FRUSTUM_NEAR] < clist[:FRUSTUM_NEAR]
          f = nlist[:FRUSTUM_FAR] > clist[:FRUSTUM_FAR]

          #lr = l>r?l:r
          #tb = t>b?t:b
          #nf = n<f?n:f

          if lt
            if l
              next = parent.next.left
              if next == nothing parent.next.left = current
              else parent = next
              end
            elseif r
              next = parent.next.right
              if next == nothing parent.next.right = current
              else parent = next
              end
            end
          else
            if t
              next = parent.next.top
              if next == nothing parent.next.top = current
              else parent = next
              end
            elseif b
              next = parent.next.bottom
              if next == nothing parent.next.bottom = current
              else parent = next
              end
            end
          end
        end

      end
    end
    =#

  end
  end

  #printList(list,length(list)/10)

  #for (k,z) in list
    #p = k1 + Vec3f(0,0,0)
    #l = k1 + Vec3f(-1,-1,-1)
    #r = k1 + Vec3f(1,1,1)

    #skip = true
    #for (k2,b) in z
      #if skip skip=false; continue end
      #if k2 == nothing continue end
      #if p.z - 10 >= k2.z continue end
      #if p.y + 1 >= k2.y continue end
      #if p.y - 1 <= k2.y continue end
      #if l.x <= k2.x continue end
      #if r.x >= k2.x continue end
      #setVisible(b,false)
    #end
  #end

end

#------------------------------------------------------------------------------------

setFilteredChilds(this::Chunk, r::Array{Block,1}) = begin this.filtered = r; this.fileredCount = length(r); r end
getFilteredChilds(this::Chunk) = this.filtered

getActiveChilds(this::Chunk) = filter(b->isActive(b),this.childs)
getVisibleChilds(this::Chunk) = filter(b->isVisible(b),this.childs)
getValidChilds(this::Chunk) = filter(b->isValid(b),this.childs)

function getData(this::Block)
  # get visible sides
  i=0; sides=0
  for side in this.sides
    if side > 0 sides |= (0x1 << i) end
    i+=1
  end

  #if this.sides[LEFT_SIDE] > 0 sides |= 0x1 end
  #if this.sides[RIGHT_SIDE] > 0 sides |= 0x2 end
  #if this.sides[TOP_SIDE] > 0 sides |= 0x4 end
  #if this.sides[BOTTOM_SIDE] > 0 sides |= 0x8 end
  #if this.sides[FRONT_SIDE] > 0 sides |= 0x10 end
  #if this.sides[BACK_SIDE] > 0 sides |= 0x20 end

  SVector(Float32[this.pos...,this.typ,sides]...)
end

getData(this::Chunk) = isValid(this)?vec((b->getData(b)).(getFilteredChilds(this))):Float32[]

function update(this::Chunk)
  this.count = length(this.childs)
  setFilteredChilds(this, getValidChilds(this))

  #fill copies references
  #refblocks = Ptr{Void}[]
  #for b in this.childs push!(refblocks, pointer_from_objref(b)) end
end

function createSingle(this::Chunk)
  b=this.childs[1,1,1]
  b.pos=Vec3f(0, 0, -10)
  b.typ=1
end

function createExample(this::Chunk)
  const DIST = Vec3f(2,2,2); #Vec3f(3,5,3) #r = 1f0/30 #0.005f0
  const START = Vec3f(-(this.len*DIST.x) / 2.0f0, -(this.len*DIST.y) / 2.0f0, (this.len*DIST.z) / 2.0f0)

  translate(x,y,z) = Vec3f(START.x+x*DIST.x, START.y+y*DIST.y, START.z-z*DIST.z)

  x=1; y=1; z=1;

  for i=1:(this.len^3)
    b = this.childs[i];
    b.pos=translate(x-1,y-1,z-1)
    b.typ=rand(1:16)
    #model = Mat4x4f(translation(Array(b.pos)))
    #push!(refblocks, pointer_from_objref(b))

    x += 1; if x > this.len
      y += 1; x=1;
      if y > this.len z += 1; y=1;
        if z > this.len break; end
      end
    end
  end
end

function createLandscape(this::Chunk)
  #Texture* heightTexture = m_pRenderer->GetTexture(m_pChunkManager->GetHeightMapTexture());

  rgbImage=Images.load("heightmap.jpg")
  (imgwidth, imgheight) = size(rgbImage)

  grayImage = Gray.(rgbImage)
  imgMatrix = reinterpret.(channelview(grayImage)) #convert(Array{Int32},raw(grayImage))

  w = imgwidth/this.len
  h = imgheight/this.len

  const DIST = Vec3f(2,2,2)
  const START = Vec3f(-(this.len*DIST.x) / 2.0f0, -(this.len*DIST.y) / 2.0f0, (this.len*DIST.z) / 2.0f0)

  translate(x,y,z) = Vec3f(START.x+x*DIST.x, START.y+y*DIST.y, START.z-z*DIST.z)

  for x=1:this.len
    for z=1:this.len
      # Use the height map texture to get the height value of x, z
      height = (imgMatrix[UInt32(trunc(z*h)), UInt32(trunc(x*w))] / 0xFF)* 1.25 * this.len
      if height >= this.len height=this.len end

      level_air = height * 0.95
      level_grass = height * 0.9
      level_dirt = height * 0.8
      level_stonebricks = height * 0.7
      level_stone = height * 0.4
      level_lava = height * 0

      for y=1:UInt32(round(height))
        #id = y*this.len^2+z*this.len+x
        #c = this.childs[id]
        #c.active = true
        b = this.childs[x,y,z]
        b.pos = translate(x-1,y-1,z-1)
        #if y >= level_air b.typ = 0 # air or nothing
        if y >= level_grass  b.typ = 2 #grass
        elseif y >= level_dirt b.typ = 1 #dirt
        elseif y >= level_stonebricks b.typ = 5 #stonebricks
        elseif y >= level_stone b.typ = 4 #stone
        elseif y >= level_lava b.typ = 15 #lava
        end
      end
    end
  end

  h=1 #UInt32(trunc(this.len/2))
  b = this.childs[h,h,h]
  b.pos=Vec3f(0,0,0)
  b.typ=7
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
mutable struct Metropolis{T <: AbstractFloat, S}
  state::S
  logprob_state::T
end

function update!(m::Metropolis{T, S},
                 logprob::Function,
                 logjump::Function,
                 randjump::Function;
                 verbose::Int=0) where {T <: AbstractFloat, S}
  # Draw candidate from jump (given current state)
  cand = randjump(m.state)

  # Compute log density of candidate draw
  logprob_cand = logprob(cand)

  # Compute log porposal density of prev -> candidate
  logjump_cand_from_prev = logjump(cand, m.state)

  # Compute log porposal density of candidate -> prev
  logjump_prev_from_cand = logjump(m.state, cand)

  # Compute log acceptance ratio
  log_accept_ratio = ((logprob_cand - m.logprob_state) -
                      (logjump_cand_from_prev - logjump_prev_from_cand))

  # Log of rand uniform
  log_u = log(rand())

  # With some probability, update current state with proposed state.
  if log_accept_ratio > log_u
    if verbose > 0
      alpha = round(log_accept_ratio, digits=3)
      println("Accepted proposal with log_accept_ratio = $(alpha)")
    end
    m.state = deepcopy(cand)
    m.logprob_state = logprob_cand
  elseif verbose > 0
    alpha = round(log_accept_ratio, digits=3)
    println("Rejected proposal with log_accept_ratio = $(alpha)")
  end

  return m
end


text:
using FastTransforms, LowRankApprox, Compat
using Compat.Test, Compat.LinearAlgebra

if VERSION < v"0.7-"
    const mul! = Base.A_mul_B!
    const lmul! = Base.scale!
    const rmul! = Base.scale!
else
    Ac_mul_B!(C, A, B) = mul!(C, A', B)
    At_mul_B!(C, A, B) = mul!(C, transpose(A), B)
end

import FastTransforms: Butterfly

@testset "Butterfly algorithm" begin
    kernel = (x,y) -> exp(im*x*y)

    function randfft(m,n,σ)
        x = (range(0,stop=m-1,length=m)+σ*rand(m))
        y = (range(0,stop=n-1,length=n)+σ*rand(n))*2π/n
        Complex{Float64}[kernel(x,y) for x in x, y in y]
    end

    function randnfft(m,n,σ)
        x = (range(0,stop=m-1,length=m)+σ*randn(m))
        y = (range(0,stop=n-1,length=n)+σ*randn(n))*2π/n
        Complex{Float64}[kernel(x,y) for x in x, y in y]
    end

    println("Testing the butterflied FFT")

    N = 10
    A = Vector{Matrix{Complex{Float64}}}(undef,N)
    for n in 1:N
        A[n] = randnfft(2^n,2^n,0.0)
    end

    for n in 7:N
        @time B = Butterfly(A[n], n-5)
        nb = size(A[n], 2)
        b = rand(Complex{Float64}, nb)./(1:nb)
        u = zero(b)
        @time uf = A[n]*b
        @time mul!(u, B, b)
        w = zero(b)
        @time Ac_mul_B!(w, B, u)
        lmul!(inv(2^n), w)
        println(norm(u-uf)/2^n)
        println(norm(w-b))
        println(norm(w-A[n]\u))
    end

    println("Testing the butterflied Hilbert matrix")

    N = 10
    A = Vector{Matrix{Float64}}(undef,N)
    for n in 1:N
        A[n] = Float64[1/(i+j-1) for i = 1:2^n+50,j=1:2^n+50]
    end

    for n in 7:N
        @time B = Butterfly(A[n], n-5)
        nb = size(A[n], 2)
        b = rand(Float64, nb)./(1:nb)
        u = zero(b)
        @time uf = A[n]*b
        @time mul!(u, B, b)
        w = zero(b)
        @time At_mul_B!(w, B, b)
        println(norm(u-uf)/nb)
        println(norm(w-A[n]'b))
        println(norm(u-w))
    end

    println("Testing the butterflied NUFFT")

    N = 10
    A = Vector{Matrix{Complex{Float64}}}(undef,N)
    B = Vector{Butterfly{Complex{Float64}}}(undef,N)
    for n in 7:N
        A[n] = randnfft(2^n+50,2^n+50,0.1)
        @time B[n] = Butterfly(A[n], n-6)
    end

    for n in 7:N
        nb = size(A[n], 2)
        b = rand(Complex{Float64}, nb)./(1:nb)
        uf = zero(b)
        u = zero(b)
        @time mul!(uf, A[n], b)
        @time mul!(u, B[n], b)
        println(norm(u-uf))
    end
end


text:
@testset "MC: IsingModel Simulation" begin
    Random.seed!(123)
    m = IsingModel(dims=2, L=8);
    mc = MC(m, beta=0.35);
    run!(mc, sweeps=1000, thermalization=10, verbose=false);

    # Check measurements
    measured = measurements(mc)[:ME]
    @test   25.47  ≈ measured[:Magn].M |> mean         atol=0.01
    @test    0.82  ≈ measured[:Magn].M |> std_error    atol=0.01
    @test  887.    ≈ measured[:Magn].M2 |> mean        atol=1.0
    @test   46.    ≈ measured[:Magn].M2 |> std_error   atol=1.0
    @test    0.398 ≈ measured[:Magn].m |> mean         atol=0.001
    @test    0.013 ≈ measured[:Magn].m |> std_error    atol=0.001
    @test    1.300 ≈ measured[:Magn].chi |> mean       atol=0.001

    @test  -55.10  ≈ measured[:Energy].E |> mean       atol=0.01
    @test    0.88  ≈ measured[:Energy].E |> std_error  atol=0.01
    @test 3342.    ≈ measured[:Energy].E2 |> mean      atol=1.0
    @test  104.    ≈ measured[:Energy].E2 |> std_error atol=1.0
    @test   -0.861 ≈ measured[:Energy].e |> mean       atol=0.001
    @test    0.014 ≈ measured[:Energy].e |> std_error  atol=0.001
    @test    0.585 ≈ measured[:Energy].C |> mean       atol=0.001

    @test [
        -0.028 -0.02 0.022 0.018 0.002 -0.008 -0.026 -0.058; -0.016 -0.028 0.004 0.014 0.024 0.024 -0.026 -0.044; -0.048 -0.018 0.004 0.02 0.02 0.002 -0.016 -0.018; -0.028 -0.016 0.016 -0.02 -0.018 -0.024 -0.026 -0.03; -0.02 -0.016 -0.024 -0.044 -0.028 -0.006 -0.026 -0.05; -0.008 -0.006 0.004 -0.038 -0.038 -0.016 -0.03 -0.022; -0.028 -0.018 0.014 -0.01 -0.018 0.008 0.008 -0.05; -0.038 -0.028 0.016 0.026 -0.018 -0.022 -0.018 -0.054
    ] ≈ measured[:conf].obs |> mean                    atol=0.001
    @test last(measured[:conf].obs) == mc.conf
end

@testset "DQMC: HubbardModel Simulation" begin
    Random.seed!(123)
    m = HubbardModelAttractive(dims=2, L=4);
    mc = DQMC(m, beta=1.0);
    run!(mc, sweeps=1000, thermalization=10, verbose=false);

    # Check measurements
    measured = measurements(mc)[:ME]
    @test [
        0.477228 -0.148092 0.000650671 -0.152739 -0.156428 0.0050096 0.0343558 0.00164354 0.000686254 0.0314209 -0.00315709 0.0363441 -0.149735 9.70741e-5 0.0335553 0.00358477; -0.153774 0.454577 -0.160642 0.000648542 0.0126497 -0.159687 0.00669076 0.0333576 0.037061 0.00569243 0.0358175 -0.00404491 0.00355269 -0.16096 0.0045801 0.0340806; 0.00137329 -0.149875 0.525756 -0.156449 0.0352624 0.000278655 -0.152092 -0.00291045 0.000989166 0.0322977 -0.00177604 0.0358907 0.0322688 -0.00133419 -0.151118 0.00029095; -0.145036 0.000208847 -0.148504 0.520469 0.00232344 0.0322997 -0.00403395 -0.141013 0.0340725 -0.00152847 0.0315167 -0.0043278 -0.0022816 0.032141 -0.00106526 -0.149073; -0.153225 0.0116195 0.0344808 0.00197823 0.459859 -0.154248 0.00721501 -0.148024 -0.157576 0.00551732 0.0349254 0.00233362 0.00219748 0.0371333 -0.00627411 0.03385; 0.00507538 -0.151558 -0.000391151 0.0329975 -0.155022 0.479573 -0.150905 0.00226697 0.00295499 -0.148607 0.00301816 0.0347823 0.0315428 0.00392446 0.0326532 -0.00249034; 0.0351943 0.00306724 -0.159841 -0.00715043 0.00575413 -0.155764 0.491007 -0.148957 0.0353336 0.00576441 -0.153046 -0.00141562 -0.000898893 0.0362212 0.00138285 0.0337357; 0.00118721 0.0348373 -0.00172386 -0.166127 -0.168602 0.000460718 -0.16182 0.508259 -0.000832228 0.036381 0.000702338 -0.165049 0.0363685 0.000530091 0.0375223 3.16583e-7; -0.000531546 0.031598 0.000765815 0.0349322 -0.150373 0.00357689 0.032698 -0.00108294 0.507991 -0.14935 0.00185427 -0.154108 -0.147832 -0.00104978 0.032203 -0.000816359; 0.0354625 0.00767388 0.0356703 0.000856574 0.00576076 -0.152291 0.0023841 0.0328979 -0.157559 0.475779 -0.156894 -0.00115001 -0.00177241 -0.158056 0.00557598 0.0338655; -0.00215687 0.0327828 -0.00226684 0.0351669 0.0352777 0.00517752 -0.149678 -0.000851672 0.00202061 -0.151299 0.480858 -0.157663 0.0346624 6.83759e-5 -0.155145 0.000856897; 0.0346939 -0.00253293 0.0333441 -0.00604537 -0.000972427 0.034324 -0.000290598 -0.147195 -0.160874 -0.00191725 -0.157674 0.52118 -0.00359494 0.0360121 0.00140374 -0.153626; -0.156798 0.00415648 0.0344931 -0.00219987 0.00366139 0.0365817 -0.00100057 0.0324253 -0.163216 -0.00268529 0.0365863 -0.00464487 0.526158 -0.162021 0.00374309 -0.151599; -0.000494114 -0.148824 -0.00554303 0.0321295 0.0336398 0.00246292 0.0314112 -0.00022372 -0.00163466 -0.149309 0.00169644 0.0332757 -0.149364 0.494472 -0.146979 -0.001565; 0.0343884 0.00516688 -0.153669 -0.00249342 -0.00596883 0.0360426 0.00222142 0.0311754 0.0350533 0.00385871 -0.154668 0.00290686 4.46472e-5 -0.155343 0.498268 -0.146554; 0.00117938 0.0336081 0.000884403 -0.159131 0.0366746 -0.00326657 0.0348769 0.00402202 -0.00051693 0.0347719 0.00395268 -0.160173 -0.154205 0.000746075 -0.154492 0.47354
    ] ≈ measured[:Greens].obs |> mean           atol=0.0001
    @test [
        0.0205191 0.00466986 0.00319713 0.00433101 0.0063448 0.00298264 0.00188447 0.00270596 0.00279388 0.00147149 0.00223948 0.00191749 0.00544694 0.0034801 0.00166324 0.00279062; 0.00581179 0.0173809 0.00615989 0.00298333 0.00264239 0.00544051 0.00288374 0.00180361 0.00208151 0.00286832 0.0019582 0.0017646 0.00267498 0.00597742 0.00256961 0.00180452; 0.00304324 0.00535371 0.0171345 0.00603167 0.00188977 0.00215002 0.00496448 0.00224604 0.0019083 0.00161783 0.00226727 0.00188134 0.00147511 0.00266665 0.00549195 0.00243711; 0.00409947 0.00231981 0.00510606 0.0172372 0.00209903 0.00155016 0.00232761 0.00370052 0.0018203 0.00180111 0.00127228 0.00238768 0.00236927 0.00137936 0.00228196 0.00505507; 0.0059045 0.0025362 0.00187802 0.00270614 0.0191752 0.00584846 0.00234163 0.00566209 0.00677845 0.00217692 0.001687 0.00235139 0.00203264 0.00186346 0.00162346 0.00230803; 0.00271169 0.00545808 0.00263302 0.0016385 0.00582896 0.0177921 0.00586849 0.00208916 0.00277508 0.00461105 0.00264704 0.00144147 0.00176448 0.00242271 0.0019095 0.00155428; 0.00196409 0.00206585 0.00588655 0.00288169 0.00244775 0.00538384 0.0147285 0.00616923 0.00185163 0.00234098 0.00470289 0.00209398 0.00170172 0.00177406 0.00201058 0.00195253; 0.00291793 0.00171317 0.00288986 0.00509084 0.00732493 0.00205185 0.00623073 0.0158797 0.00265401 0.00204915 0.00334973 0.00524786 0.00180304 0.00206516 0.00213687 0.00277024; 0.00228005 0.00158917 0.00184363 0.00168331 0.00600594 0.0023704 0.00159086 0.00215774 0.0179357 0.00458805 0.00265191 0.0057394 0.0048724 0.00230553 0.00160043 0.00256597; 0.00188113 0.00260856 0.00165067 0.00188637 0.00246927 0.00535851 0.00277197 0.00174793 0.00452975 0.0186795 0.00587128 0.00293908 0.00256491 0.00550431 0.00254789 0.00182041; 0.00199224 0.00198167 0.00239671 0.00148865 0.0019963 0.00225782 0.0048999 0.00253693 0.00227267 0.00495478 0.0167915 0.00544362 0.00193831 0.00261395 0.00605124 0.00213894; 0.00201428 0.00189837 0.00177325 0.00217569 0.00219486 0.00157229 0.00262323 0.00533245 0.00639645 0.00236338 0.00582814 0.0141972 0.00214637 0.00170952 0.00242774 0.00652343; 0.00551784 0.00314772 0.0018408 0.00244508 0.00239146 0.00211826 0.00165435 0.00167652 0.00526101 0.00247197 0.00198295 0.00232103 0.0166794 0.00479359 0.00235695 0.0051415; 0.00299558 0.00716651 0.00302339 0.00157788 0.00177353 0.00249744 0.00176033 0.00164171 0.00242951 0.00543594 0.00212568 0.00147208 0.00535451 0.0167751 0.00525319 0.00226052; 0.00194073 0.0021867 0.00667736 0.00247544 0.00198623 0.00230368 0.00222494 0.00156614 0.00178054 0.00234639 0.00585425 0.00259014 0.00238915 0.00565968 0.0169088 0.00539021; 0.00288388 0.00177993 0.00292207 0.00446213 0.00215593 0.00206059 0.00171043 0.00265896 0.00261172 0.00179292 0.00246213 0.00612522 0.00551238 0.00288897 0.00570215 0.0156827
    ] ≈ measured[:Greens].obs |> std_error      atol=0.0001

    @test  0.784 ≈ measured[:BosonEnergy].obs |> mean           atol=0.001
    @test  0.430 ≈ measured[:BosonEnergy].obs |> std_error      atol=0.001

    @test [
        0.16 -0.06 0.1 0.04 0.16 -0.1 -0.06 0.16 0.04 0.0; -0.04 0.1 0.14 0.16 0.22 0.32 0.06 0.12 0.06 -2.22045e-18; -0.16 0.14 -0.1 -0.02 0.02 0.1 -0.22 -0.16 -0.04 -0.08; 0.04 -0.06 -0.08 -0.06 -0.08 -0.16 -0.08 -0.1 -0.04 0.06; 0.1 0.12 0.14 0.18 0.14 0.08 -2.22045e-18 0.18 -0.02 0.04; 0.22 0.04 0.22 -0.16 0.04 -0.08 -0.04 0.06 0.12 9.99201e-18; -0.08 0.06 0.16 0.08 0.04 -0.06 -0.08 0.18 0.0 -0.04; -0.2 0.04 0.0 0.06 -0.06 0.04 -2.10942e-17 -0.06 0.18 -0.16; -0.1 0.06 0.1 -0.12 -0.08 0.08 -8.88178e-18 -0.12 -0.08 0.1; 0.16 -0.04 0.06 -0.02 -0.1 0.1 0.16 -0.06 0.16 0.04; 0.08 0.1 -0.12 0.06 -0.06 0.2 0.06 -0.06 0.04 0.04; 0.08 -0.2 -0.04 -0.02 -0.12 -0.08 -0.1 0.14 -0.14 -0.02; -0.06 -0.1 -0.02 0.06 -0.08 -0.08 -0.02 0.1 -0.14 -0.16; 0.1 -0.08 -0.08 0.12 0.04 -0.02 -0.22 0.08 -0.06 0.12; -0.16 0.04 0.16 0.12 0.08 0.1 0.22 -0.16 -0.12 -0.02; -2.22045e-18 0.06 0.02 0.04 0.1 0.06 0.22 0.16 -0.06 -0.02
    ] ≈ measured[:conf].obs |> mean                    atol=0.01
    @test last(measured[:conf].obs) == mc.conf
end


text:
using Test

@testset "Test GraphOrdering.jl" begin
    @testset "Test bandwidth calculation" begin
        include("test_bw.jl")
    end
    @testset "Test symrcm" begin
        include("test_symrcm.jl")
    end
end


text:
@everywhere using IterativeSolvers.lsqr

# Pulls K2vec and sparsify_codes
include("utils.jl")

# Update a dimension of a codebook using LSQR or LSMR
@everywhere function updatecb!(
  K::SharedMatrix{Float32},
  C::SparseMatrixCSC{Int32,Int32},
  X::Matrix{Float32},
  IDX::UnitRange{Int64},
  codebook_upd_method::AbstractString="lsqr")   # choose the codebook update method out of lsqr or lsmr

  if codebook_upd_method == "lsqr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsqr( C, vec(X[i, :]) );
    end
  elseif codebook_upd_method == "lsmr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsmr( C, vec(X[i, :]) );
    end
  else
    error("Codebook update method unknown: ", codebook_upd_method)
  end
end

@everywhere function updatecb!(
  K::SharedMatrix{Float32},
  C::SparseMatrixCSC{Float32,Int32},
  X::Matrix{Float32},
  IDX::UnitRange{Int64},
  codebook_upd_method::AbstractString="lsqr")   # choose the codebook update method out of lsqr or lsmr

  if codebook_upd_method == "lsqr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsqr( C, vec(X[i, :]) );
    end
  elseif codebook_upd_method == "lsmr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsmr( C, vec(X[i, :]) );
    end
  else
    error("Codebook update method unknown: ", codebook_upd_method)
  end
end

###############################################
## Codebook update for a fully-connected MRF ##
###############################################

function update_codebooks(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false,      # whether to print progress
  codebook_upd_method::AbstractString="lsqr")   # choose the codebook update method out of lsqr or lsmr

  if !(codebook_upd_method in ["lsmr", "lsqr"]); error("Codebook update method unknown"); end

  if V print("Doing " * uppercase(codebook_upd_method) * " codebook update... "); st=time(); end

  d, n   = size(X)
  m, _   = size(B)
  C = sparsify_codes( B, h )

  K = SharedArray{Float32, 2}((d, size(C,2)))
  if nworkers() == 1
    updatecb!( K, C, X, 1:d )
  else
    paridx = splitarray( 1:d, nworkers() )
    @sync begin
      for (i, wpid) in enumerate( workers() )
        @async begin
          remotecall_wait( updatecb!, wpid, K, C, X, paridx[i], codebook_upd_method )
        end #@async
      end #for
    end #@sync
  end

  new_C = K2vec( K, m, h )

  if V @printf("done in %.3f seconds.\n", time()-st); end

  return new_C
end # function update_codebooks

function get_cbdims_chain(
  d::Integer, # The number of dimensions
  m::Integer) # The number of codebooks

  subdims = splitarray(1:d, m-1);
  odims = Vector{UnitRange{Integer}}(m);

  odims[1] =  subdims[1];
  for i = 2:m-1;
    odims[i] = subdims[i-1][1]:subdims[i][end];
  end
  odims[end] = subdims[end];

  return odims
end

# Update a dimension of a codebook using LSQR
@everywhere function updatecb_struct!(
  K::SharedMatrix{Float32},
  C::SparseMatrixCSC{Float32,Int32},
  X::Matrix{Float32},
  dim2C, #::Vector{Bool},
  subcbs,
  IDX::UnitRange{Int64} )

  for i = IDX
    rcbs      = cat( 1, subcbs[find(dim2C[i,:])]... );
    K[i,rcbs] = lsqr( C[:,rcbs], vec(X[i, :]) );
  end
end

function update_codebooks_generic(
  X::Matrix{Float32},  # d-by-n. The data that was encoded.
  B::Union{Matrix{Int16},SharedArray{Int16,2}}, # d-by-m. X encoded.
  h::Integer,          # number of entries per codebooks
  odimsfunc::Function, # Function that says which dimensions each codebook has
  V::Bool=false)       # whether to print progress

  if V print("Doing LSQR codebook update... "); st=time(); end

  d, n   = size( X );
  m, _   = size( B );
  C = sparsify_codes( B, h );

  odims = odimsfunc(d, m);

  # Make a map of dimensions to codebooks
  dim2C = zeros(Bool, d, m);
  for i = 1:m; dim2C[ odims[i], i ] = true; end

  K = SharedArray{Float32,2}((d, h*m));
  subcbs = splitarray(1:(h*m), m);

  if nworkers() == 1
    updatecb_struct!( K, C, X, dim2C, subcbs, 1:d );
  else
    paridx = splitarray( 1:d, nworkers() );
    @sync begin
      for (i, wpid) in enumerate( workers() )
        @async begin
          remotecall_wait(updatecb_struct!, wpid, K, C, X, dim2C, subcbs, paridx[i] );
        end #@async
      end #for
    end #@sync
  end

  new_C = K2vec( K, m, h );
  if V @printf("done in %.3f seconds.\n", time()-st); end

  return new_C
end

# Update codebooks in a chain
function update_codebooks_chain(
  X::Matrix{Float32}, # d-by-n. The data that was encoded.
  B::Union{Matrix{Int16},SharedArray{Int16,2}}, # d-by-m. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false)      # whether to print progress

  return update_codebooks_generic(X, B, h, get_cbdims_chain, V);

end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
function FmeanSquare(x...)
    tmp=0.0
    for i in 1:length(x)
        tmp+=mean(x[i].*x[i])
    end
    return ([tmp/length(x)],nothing) # must always return an array for the value
end

function FmeanSquare_inplace(value,auxvalue,x...) # inplace
    tmp=0.0
    for i in 1:length(x)
        tmp+=mean(x[i].*x[i])
    end
    copy!(value,tmp/length(x))
end

DmeanSquare(derivativeIDX,f_c,faux_c,grad_c,grad_n,x...)=axpy!(2.0/(length(x)*length(x[derivativeIDX])),grad_c.*x[derivativeIDX],grad_n)

if PROC=="GPU"
    function FmeanSquare()
    end

    function FmeanSquare(x::CudaArray{Float64}...)
        tmp=CudaArray(Float64,zeros(1))
        for i in 1:length(x)
            CUBLAS.gemv!('T',1.0/(length(x)*length(x[i])),flatten(Float64,x[i]),vec(x[i]),1.0,tmp)
        end
    return (tmp,nothing) # must always return an array for the value
    end

    function FmeanSquare(x::CudaArray{Float32}...)
        tmp=CudaArray(Float32,zeros(1))
        for i in 1:length(x)
            CUBLAS.gemv!('T',Float32(1.0/(length(x)*length(x[i]))),flatten(Float32,x[i]),vec(x[i]),Float32(1.0),tmp)
        end
    return (tmp,nothing) # must always return an array for the value
    end

    function FmeanSquare_inplace(value::CudaArray{Float64},auxvalue,x::CudaArray{Float64}...) # inplace
        fill!(value,0.0)
        for i in 1:length(x)
            CUBLAS.gemv!('T',1.0/(length(x)*length(x[i])),flatten(Float64,x[i]),vec(x[i]),1.0,value);
        end
    end

    function FmeanSquare_inplace(value::CudaArray{Float32},auxvalue,x::CudaArray{Float32}...) # inplace
        fill!(value,Float32(0.0))
        for i in 1:length(x)
            CUBLAS.gemv!('T',Float32(1.0/(length(x)*length(x[i]))),flatten(Float32,x[i]),vec(x[i]),Float32(1.0),value);
        end
    end

    function DmeanSquare(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray...)
        alphaaxpy!(2.0/(length(x)*length(x[derivativeIDX])),grad_c,x[derivativeIDX],grad_n)
    end

end

Inplace[FmeanSquare]=FmeanSquare_inplace
Derivative[FmeanSquare]=DmeanSquare

meanSquare(n::ADnode)=ADnode(FmeanSquare,n)

function meanSquare(n::ArrayADnode)
    return ADnode(FmeanSquare,n)
end

#meanSquare(A::ADtrans)=ADnode(FmeanSquare, ftranspose(node[A.parent]))
meanSquare(A::ADtrans)=ADnode(FmeanSquare, node[A.parent]) # meanSq(A')=meanSq(A)

export meanSquare

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function inverse(side::HorizonSide)
    new_side=HorizonSide(mod(Int(side)+2,4))
    return new_side
end

function movements!(r,side,n::Int)
    for n in 1:n+1
        move!(r,side)
    end
end

function  shuttle!(r, side)
    n=0
    side=Ost
    while isborder(r,Nord)
        n += 1
        movements!(r, side, n)
        side = inverse(side)
    end
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
using Gadfly
draw(SVG("gadfly.svg",6inch,3inch), plot([x -> x^2],0, 25))
pl = plot([x -> cos(x)/x], 5, 25)
draw(PNG("gadfly.png", 300, 100), pl)

max_stars_repo_path: Module 1/Chapter10/gadfly.jl

max_issues_repo_path: Module 1/Chapter10/Old/gadfly.jl

max_forks_repo_path: Module 1/Chapter10/Old/gadfly.jl


text:
module PositionTypes

__precompile__(true)

export DTE, IV, Options, Underlying, Position
export to_τ, get_ivol, set_DTE, set_ivol

mutable struct DTE{T<:Real}
    days::T
end

mutable struct IV{T<:Real}
    imp_vol::T
end

struct Options{T<:Real}
    isCall::Bool
    strike::T
    expiry::DTE
    ivol::IV
end

struct Underlying{T<:Real}
    price::T
end

struct Position{T<:Real, S}
    holding::Union{Options{T}, Underlying{T}}
    quantity::S
end

to_τ(x::Options) = x.expiry.days / 365.0 # We use ACT/365 here
get_ivol(x::Options) = x.ivol.imp_vol

function set_DTE(x::Options, new_days)
    @assert new_days != 0 "For Greek calculation DTE cannot be 0"
    x.expiry.days = new_days
end

function set_ivol(x::Options, new_ivol)
    @assert new_ivol != 0 "For Greek calculation IVol cannot be 0"
    x.ivol.imp_vol = new_ivol
end

end

max_stars_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_issues_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_forks_repo_name: Stochastic-Adventure/BlackScholesGreeks


text:
module Custom

using MLJ, MLJScikitLearnInterface, DataFrames

export MLJ, MLJScikitLearnInterface, transform

function transform(data, model)
    select!(data, Not(:class))
end

end

max_stars_repo_name: amperie/user-models

max_issues_repo_name: amperie/user-models

max_forks_repo_name: amperie/user-models


text:
module MatchingTests
using  ExpressionPatterns.Matching
using  ExpressionPatterns.Helper
using  Base.Test
import Base.Meta.quot

matchval(exval) = :(@test match($(exval[1])) == $(exval[2]))

macro testmatch(pattern, block)

  arg_statements = (x->esc.(x.args[2:end])).(linesof(block))
  res_statements = matchval.(arg_statements)

  :(let match = $matcher($(esc(pattern)))
      $(res_statements...)
    end)
end

# anything

@testmatch (:x) begin
   1000   => true
  :(x,y)  => true
  :(x+10) => true
end

# tuple

@testmatch :(x,y,z) begin
  :(a,b,c) => true
  :(1,2,3) => true
  :[a,b,c] => false
  :(1,2)   => false
end

# vector / literal

@testmatch :[1,a] begin
  :[1,2]        => true
  :[1,:(a,b,c)] => true
  :[1,b,c]      => false
  :[a,a]        => false
  :[[1],a]      => false
  :(1,2)        => false
end

# call

@testmatch :(f(a,b)) begin
  :(1+2)                  => true
  :(reduce(*, [1,2,3,4])) => true
  :(g[1,2])               => false
  :(t{a,b})               => false
end

# ref

@testmatch :((:T{Symbol})[:T{Integer}]) begin
  :(A[1])      => true
  :(arr[1000]) => true
  :(A[1,2])    => false
  :(1[1])      => false
  :(A[a])      => false
end

# curly

@testmatch (:(T{*{a}})) begin
  :(Vector{Integer}) => true
  :(A{1,2,3,4})      => true
  :(T{})             => true
end

# block

@testmatch (quote *{:T{:P{s->all(islower,s)}, AbstractString}} end) begin
  (quote "abc"; "onetwothree"; "solami" end) => true
  (quote "abc";     "123";     "solami" end) => false
  (quote "abc"; "onetwothree"; "solami" end) => true
  (quote "aBc"; "onetwothree"; "solami" end) => false
  (quote "abc";      123;      "solami" end) => false
end

# forloop / nested slurps

@testmatch :(for i in [*{[*{1,2}]}]; *{body} end) begin
  :(for x in [[1,2],[],[1,2,1,2]]; print(x) end) => true
  :(for x in [[1,2],[],[1,2,1]];   print(x) end) => false
  :(for x in [[1,2],[3],[1,2,1]];  print(x) end) => false
end

# while loop / lazy slurps

@testmatch :(while x < t; :?{p(:?{args})}; x+=1 end) begin

  :(while i <   10; println(io, i);             i+=1 end) => true
  :(while k < t-10; splice!(A,2,k); push!(B,k); k+=1 end) => true
  :(while x < g(x); println(x);     push!(B,x); x+=1 end) => true
  :(while x < g(x); println(i);     push!(B,x); x+=1 end) => true
  :(while x < g(x); println(x,y);   push!(B,x); i+=1 end) => false
  :(while x > g(x); println(x);     push!(B,x); x+=2 end) => false

end

# macrocall / quote

@testmatch :(@m :x y) begin
  :(@transport :car   (john, peter))    => true
  :(@transport plane  (maria, isabela)) => false
end

# binding

@testmatch (:(:B{+}(1+2))) begin
  :(f(1+2)) => true
  :(f(1-2)) => false
end

# literal

@testmatch :(:L{x+y}) begin
  :(x+y) => true
  :(1+2) => false
end

# autobinding

@testmatch :(:L{:A{x+y}}) begin
  :(1+2) => true
  :(x+y) => true
  :(x-y) => false
end

# consistency

@testmatch :(x, x) begin
  :(1,1) => true
  :(1,2) => false
end

# equality

@testmatch :(:EQ{x, 10}) begin
  10 => true
  11 => false
end

# type check

@testmatch :(:T{x, Symbol}) begin
  :x => true
  10 => false
end

# predicate check

@testmatch :(:P{x, iseven}) begin
  2 => true
  1 => false
end

# iterable

@testmatch :(:I{x,y,z}) begin
  (1,2,3) => true
  [4,5,6] => true
 :(7,8,9) => false
 :[a,b,c] => false
end

# raw

@testmatch :(:R{:export, x, y, z}) begin
  :(export a,b,c) => true
  :(import a.b.c) => false
end

# line number nodes

@testmatch LineNumberNode(1) begin
  LineNumberNode(1) => true
  LineNumberNode(2) => true
  Expr(:line, 10)   => true
end

end


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

struct PolicyAssignmentsApi <: SwaggerApi
    client::Swagger.Client
end

"""
Creates or updates a policy assignment.
 This operation creates or updates a policy assignment with the given scope and name. Policy assignments apply to all resources contained within their scope. For example, when you assign a policy at resource group scope, that policy applies to all resources in the group.
Param: scope::String (required)
Param: policyAssignmentName::String (required)
Param: parameters::PolicyAssignment (required)
Param: api_version::String (required)
Return: PolicyAssignment
"""
function _swaggerinternal_policyAssignmentsCreate(_api::PolicyAssignmentsApi, scope::String, policyAssignmentName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicyAssignment, "/{scope}/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.path, "policyAssignmentName", policyAssignmentName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsCreate(_api::PolicyAssignmentsApi, scope::String, policyAssignmentName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsCreate(_api, scope, policyAssignmentName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsCreate(_api::PolicyAssignmentsApi, response_stream::Channel, scope::String, policyAssignmentName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsCreate(_api, scope, policyAssignmentName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Creates or updates a policy assignment.
This operation creates or updates the policy assignment with the given ID. Policy assignments made on a scope apply to all resources contained in that scope. For example, when you assign a policy to a resource group that policy applies to all resources in the group. Policy assignment IDs have this format: '{scope}/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}'. Valid scopes are: management group (format: '/providers/Microsoft.Management/managementGroups/{managementGroup}'), subscription (format: '/subscriptions/{subscriptionId}'), resource group (format: '/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}', or resource (format: '/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/{resourceProviderNamespace}/[{parentResourcePath}/]{resourceType}/{resourceName}'.
Param: policyAssignmentId::String (required)
Param: parameters::PolicyAssignment (required)
Param: api_version::String (required)
Return: PolicyAssignment
"""
function _swaggerinternal_policyAssignmentsCreateById(_api::PolicyAssignmentsApi, policyAssignmentId::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PolicyAssignment, "/{policyAssignmentId}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "policyAssignmentId", policyAssignmentId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsCreateById(_api::PolicyAssignmentsApi, policyAssignmentId::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsCreateById(_api, policyAssignmentId, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsCreateById(_api::PolicyAssignmentsApi, response_stream::Channel, policyAssignmentId::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsCreateById(_api, policyAssignmentId, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a policy assignment.
This operation deletes a policy assignment, given its name and the scope it was created in. The scope of a policy assignment is the part of its ID preceding '/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}'.
Param: scope::String (required)
Param: policyAssignmentName::String (required)
Param: api_version::String (required)
Return: PolicyAssignment
"""
function _swaggerinternal_policyAssignmentsDelete(_api::PolicyAssignmentsApi, scope::String, policyAssignmentName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", PolicyAssignment, "/{scope}/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.path, "policyAssignmentName", policyAssignmentName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsDelete(_api::PolicyAssignmentsApi, scope::String, policyAssignmentName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsDelete(_api, scope, policyAssignmentName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsDelete(_api::PolicyAssignmentsApi, response_stream::Channel, scope::String, policyAssignmentName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsDelete(_api, scope, policyAssignmentName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Deletes a policy assignment.
This operation deletes the policy with the given ID. Policy assignment IDs have this format: '{scope}/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}'. Valid formats for {scope} are: '/providers/Microsoft.Management/managementGroups/{managementGroup}' (management group), '/subscriptions/{subscriptionId}' (subscription), '/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}' (resource group), or '/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/{resourceProviderNamespace}/[{parentResourcePath}/]{resourceType}/{resourceName}' (resource).
Param: policyAssignmentId::String (required)
Param: api_version::String (required)
Return: PolicyAssignment
"""
function _swaggerinternal_policyAssignmentsDeleteById(_api::PolicyAssignmentsApi, policyAssignmentId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", PolicyAssignment, "/{policyAssignmentId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyAssignmentId", policyAssignmentId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsDeleteById(_api::PolicyAssignmentsApi, policyAssignmentId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsDeleteById(_api, policyAssignmentId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsDeleteById(_api::PolicyAssignmentsApi, response_stream::Channel, policyAssignmentId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsDeleteById(_api, policyAssignmentId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves a policy assignment.
This operation retrieves a single policy assignment, given its name and the scope it was created at.
Param: scope::String (required)
Param: policyAssignmentName::String (required)
Param: api_version::String (required)
Return: PolicyAssignment
"""
function _swaggerinternal_policyAssignmentsGet(_api::PolicyAssignmentsApi, scope::String, policyAssignmentName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyAssignment, "/{scope}/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "scope", scope)  # type String
    Swagger.set_param(_ctx.path, "policyAssignmentName", policyAssignmentName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsGet(_api::PolicyAssignmentsApi, scope::String, policyAssignmentName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsGet(_api, scope, policyAssignmentName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsGet(_api::PolicyAssignmentsApi, response_stream::Channel, scope::String, policyAssignmentName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsGet(_api, scope, policyAssignmentName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves the policy assignment with the given ID.
The operation retrieves the policy assignment with the given ID. Policy assignment IDs have this format: '{scope}/providers/Microsoft.Authorization/policyAssignments/{policyAssignmentName}'. Valid scopes are: management group (format: '/providers/Microsoft.Management/managementGroups/{managementGroup}'), subscription (format: '/subscriptions/{subscriptionId}'), resource group (format: '/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}', or resource (format: '/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/{resourceProviderNamespace}/[{parentResourcePath}/]{resourceType}/{resourceName}'.
Param: policyAssignmentId::String (required)
Param: api_version::String (required)
Return: PolicyAssignment
"""
function _swaggerinternal_policyAssignmentsGetById(_api::PolicyAssignmentsApi, policyAssignmentId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyAssignment, "/{policyAssignmentId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "policyAssignmentId", policyAssignmentId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsGetById(_api::PolicyAssignmentsApi, policyAssignmentId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsGetById(_api, policyAssignmentId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsGetById(_api::PolicyAssignmentsApi, response_stream::Channel, policyAssignmentId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsGetById(_api, policyAssignmentId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves all policy assignments that apply to a subscription.
This operation retrieves the list of all policy assignments associated with the given subscription that match the optional given $filter. Valid values for $filter are: 'atScope()' or 'policyDefinitionId eq '{value}''. If $filter is not provided, the unfiltered list includes all policy assignments associated with the subscription, including those that apply directly or from management groups that contain the given subscription, as well as any applied to objects contained within the subscription. If $filter=atScope() is provided, the returned list includes all policy assignments that apply to the subscription, which is everything in the unfiltered list except those applied to objects contained within the subscription. If $filter=policyDefinitionId eq '{value}' is provided, the returned list includes all policy assignments of the policy definition whose id is {value}.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Return: PolicyAssignmentListResult
"""
function _swaggerinternal_policyAssignmentsList(_api::PolicyAssignmentsApi, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyAssignmentListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Authorization/policyAssignments", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsList(_api::PolicyAssignmentsApi, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsList(_api, api_version, subscriptionId; filter=filter, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsList(_api::PolicyAssignmentsApi, response_stream::Channel, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsList(_api, api_version, subscriptionId; filter=filter, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves all policy assignments that apply to a management group.
This operation retrieves the list of all policy assignments applicable to the management group that match the given $filter. Valid values for $filter are: 'atScope()' or 'policyDefinitionId eq '{value}''. If $filter=atScope() is provided, the returned list includes all policy assignments that are assigned to the management group or the management group's ancestors. If $filter=policyDefinitionId eq '{value}' is provided, the returned list includes all policy assignments of the policy definition whose id is {value} that apply to the management group.
Param: managementGroupId::String (required)
Param: filter::String (required)
Param: api_version::String (required)
Return: PolicyAssignmentListResult
"""
function _swaggerinternal_policyAssignmentsListForManagementGroup(_api::PolicyAssignmentsApi, managementGroupId::String, filter::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PolicyAssignmentListResult, "/providers/Microsoft.Management/managementgroups/{managementGroupId}/providers/Microsoft.Authorization/policyAssignments", ["azure_auth"])
    Swagger.set_param(_ctx.path, "managementGroupId", managementGroupId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsListForManagementGroup(_api::PolicyAssignmentsApi, managementGroupId::String, filter::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsListForManagementGroup(_api, managementGroupId, filter, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsListForManagementGroup(_api::PolicyAssignmentsApi, response_stream::Channel, managementGroupId::String, filter::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsListForManagementGroup(_api, managementGroupId, filter, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves all policy assignments that apply to a resource.
This operation retrieves the list of all policy assignments associated with the specified resource in the given resource group and subscription that match the optional given $filter. Valid values for $filter are: 'atScope()' or 'policyDefinitionId eq '{value}''. If $filter is not provided, the unfiltered list includes all policy assignments associated with the resource, including those that apply directly or from all containing scopes, as well as any applied to resources contained within the resource. If $filter=atScope() is provided, the returned list includes all policy assignments that apply to the resource, which is everything in the unfiltered list except those applied to resources contained within the resource. If $filter=policyDefinitionId eq '{value}' is provided, the returned list includes all policy assignments of the policy definition whose id is {value} that apply to the resource. Three parameters plus the resource name are used to identify a specific resource. If the resource is not part of a parent resource (the more common case), the parent resource path should not be provided (or provided as ''). For example a web app could be specified as ({resourceProviderNamespace} == 'Microsoft.Web', {parentResourcePath} == '', {resourceType} == 'sites', {resourceName} == 'MyWebApp'). If the resource is part of a parent resource, then all parameters should be provided. For example a virtual machine DNS name could be specified as ({resourceProviderNamespace} == 'Microsoft.Compute', {parentResourcePath} == 'virtualMachines/MyVirtualMachine', {resourceType} == 'domainNames', {resourceName} == 'MyComputerName'). A convenient alternative to providing the namespace and type name separately is to provide both in the {resourceType} parameter, format: ({resourceProviderNamespace} == '', {parentResourcePath} == '', {resourceType} == 'Microsoft.Web/sites', {resourceName} == 'MyWebApp').
Param: resourceGroupName::String (required)
Param: resourceProviderNamespace::String (required)
Param: parentResourcePath::String (required)
Param: resourceType::String (required)
Param: resourceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Return: PolicyAssignmentListResult
"""
function _swaggerinternal_policyAssignmentsListForResource(_api::PolicyAssignmentsApi, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "policyAssignmentsListForResource", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "policyAssignmentsListForResource", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", PolicyAssignmentListResult, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/providers/{resourceProviderNamespace}/{parentResourcePath}/{resourceType}/{resourceName}/providers/Microsoft.Authorization/policyAssignments", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "parentResourcePath", parentResourcePath)  # type String
    Swagger.set_param(_ctx.path, "resourceType", resourceType)  # type String
    Swagger.set_param(_ctx.path, "resourceName", resourceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsListForResource(_api::PolicyAssignmentsApi, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsListForResource(_api, resourceGroupName, resourceProviderNamespace, parentResourcePath, resourceType, resourceName, api_version, subscriptionId; filter=filter, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsListForResource(_api::PolicyAssignmentsApi, response_stream::Channel, resourceGroupName::String, resourceProviderNamespace::String, parentResourcePath::String, resourceType::String, resourceName::String, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsListForResource(_api, resourceGroupName, resourceProviderNamespace, parentResourcePath, resourceType, resourceName, api_version, subscriptionId; filter=filter, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""
Retrieves all policy assignments that apply to a resource group.
This operation retrieves the list of all policy assignments associated with the given resource group in the given subscription that match the optional given $filter. Valid values for $filter are: 'atScope()' or 'policyDefinitionId eq '{value}''. If $filter is not provided, the unfiltered list includes all policy assignments associated with the resource group, including those that apply directly or apply from containing scopes, as well as any applied to resources contained within the resource group. If $filter=atScope() is provided, the returned list includes all policy assignments that apply to the resource group, which is everything in the unfiltered list except those applied to resources contained within the resource group. If $filter=policyDefinitionId eq '{value}' is provided, the returned list includes all policy assignments of the policy definition whose id is {value} that apply to the resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Return: PolicyAssignmentListResult
"""
function _swaggerinternal_policyAssignmentsListForResourceGroup(_api::PolicyAssignmentsApi, resourceGroupName::String, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "policyAssignmentsListForResourceGroup", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "policyAssignmentsListForResourceGroup", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", PolicyAssignmentListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Authorization/policyAssignments", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function policyAssignmentsListForResourceGroup(_api::PolicyAssignmentsApi, resourceGroupName::String, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsListForResourceGroup(_api, resourceGroupName, api_version, subscriptionId; filter=filter, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function policyAssignmentsListForResourceGroup(_api::PolicyAssignmentsApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; filter=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_policyAssignmentsListForResourceGroup(_api, resourceGroupName, api_version, subscriptionId; filter=filter, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export policyAssignmentsCreate, policyAssignmentsCreateById, policyAssignmentsDelete, policyAssignmentsDeleteById, policyAssignmentsGet, policyAssignmentsGetById, policyAssignmentsList, policyAssignmentsListForManagementGroup, policyAssignmentsListForResource, policyAssignmentsListForResourceGroup


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct MultilineChartContent <: ProtoType
    tag::Base.Vector{AbstractString}
    MultilineChartContent(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MultilineChartContent

mutable struct MarginChartContent_Series <: ProtoType
    value::AbstractString
    lower::AbstractString
    upper::AbstractString
    MarginChartContent_Series(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MarginChartContent_Series

mutable struct MarginChartContent <: ProtoType
    series::Base.Vector{MarginChartContent_Series}
    MarginChartContent(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MarginChartContent

mutable struct Chart <: ProtoType
    title::AbstractString
    multiline::MultilineChartContent
    margin::MarginChartContent
    Chart(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Chart
const __oneofs_Chart = Int[0,1,1]
const __oneof_names_Chart = [Symbol("content")]
meta(t::Type{Chart}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_Chart, __oneof_names_Chart, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Category <: ProtoType
    title::AbstractString
    chart::Base.Vector{Chart}
    closed::Bool
    Category(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Category

mutable struct Layout <: ProtoType
    version::Int32
    category::Base.Vector{Category}
    Layout(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Layout

export Chart, MultilineChartContent, MarginChartContent_Series, MarginChartContent, Category, Layout


text:
import TOML: parse
import Downloads: download
import JLFzf: inter_fzf
import Pkg: add, dependencies
import JSON3: read

#reglist = parse(seekstart(download("https://raw.githubusercontent.com/JuliaRegistries/General/master/Registry.toml", IOBuffer())))["packages"]
#k = [c for (b,c) in [a for (a, b) in values(reglist)]]
#v = collect(values(Dict(values(reglist))))
#v = [v[i][2] for i in 1:length(v)]

"""
Searches and returns the information of the selected package
"""
function search()
	x = inter_fzf(k, "--read0")
	if(x!="")
		z = uppercase(x[1])
		str = parse(seekstart(download("https://raw.githubusercontent.com/JuliaRegistries/General/master/$(z)/$(x)/Package.toml", IOBuffer())))["repo"][20:end-4]
		if(str[end-6:end-3] != "_jll")
			y = read(seekstart(download("https://api.github.com/repos/$(str)", IOBuffer())))[:description]
			println(x*".jl :", "\t", y)
		else
			println(x*".jll :", "\t", "A jll package")
		end
		return x
	end
end

"""
Searches and adds the selected package into the current environment
"""
function seadd()
	x = search()
	println("Do you want to add this package to the current environment? (y/n)")
	i = readline()
	if i == "y"
		Pkg.add(x)
	end
end

function installed()
    deps = dependencies()
    installs = Dict{String, VersionNumber}()
    for (uuid, dep) in deps
        dep.is_direct_dep || continue
        dep.version === nothing && continue
        installs[dep.name] = dep.version
    end
    return installs
end

"""
Search and get the latest version of a package
"""
function gelat(x::String = String(search()))
	y = "$(x[1])"
	y = uppercase(y)
	z = keys(parse(seekstart(download("https://raw.githubusercontent.com/JuliaRegistries/General/master/$(y)/$(x)/Versions.toml", IOBuffer()))))
	return maximum(VersionNumber.([i for i in z]))
end

"""
Compare the latest version of a package to the current version
"""
function chklat(x::String = String(inter_fzf(collect(keys(installed())), "--read0")))
	a = installed()[x]
	y = gelat(x)
	println("Installed version => $(a) \nLatest version => $(y)")
end

"""
Check if all the installed packages are at the latest version or not
"""
function chkall()
	a = installed()
	flag = false
	for (i, j) in a
		b = ""
		try
			b = gelat(i)
		catch e
			continue
		end
		if(j != b)
			flag = true
			println("$(i)       Latest version = $(b) \t Current version = $(j)\n")
		end
	end
	if(!flag)
		println("All packages are on the latest version")
	end
end


text:
function BPR(b, C, F, p)
    f(x) = F*(1 + b*(x/C)^p)
    df(x) = F*b*(p/C)*(x/C)^(p-1)
    ddf(x) = F*b*(p*(p-1)/(C^2))*(x/C)^(p-2)

    (f, df, ddf)
end
BPR(data::Dict) = BPR(data[:b], data[:capacity], data[:free_flow_time], data[:power])


text:
# This file was generated, do not modify it. # hide
baseurl = "https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data/glm/"

dfX = DataFrame(urldownload(baseurl * "X3.csv"))
dfYbinary = DataFrame(urldownload(baseurl * "Y3.csv"))
dfX1 = DataFrame(urldownload(baseurl * "X1.csv"))
dfY1 = DataFrame(urldownload(baseurl * "Y1.csv"));


text:
"""
Computation of Scattering Matrices

Accelerated if `nprocs()>1`.
"""
module SMatrices

export smatrix

dimension_files = (
    "1D/SMatrices1D.jl",
    )

using ..Common
import ..Common.AbstractComplexBL

using ..Scattering
using Distributed
using Interpolations
using LaTeXStrings
using LinearAlgebra
using RecipesBase
using SharedArrays
using Statistics

import ..Scattering.ScatteringSolution
import ..Scattering.EquivalentSource

"""
    ScatteringMatrix{N}

`N`×`N` Scattering Matrix container with fields `S` (scattering matrix), `ω` (frequencies)
"""
struct ScatteringMatrix{N,TW}
    S::Array{ComplexF64,3}
    ω::TW
end

"""
    SMatrixContainer{N}

Auxilliary for computation of `N`×`N` scattering matrix.
"""
struct SMatrixContainer{N,TSIM,TW,TLU} <: Function
    sim::TSIM
    ω::TW
    lupack::TLU
    ky::Float64
    kz::Float64
end

foreach(include,dimension_files)

Base.ndims(smc::SMatrixContainer{N}) where N = N

"""
    smatrix(::Simulation, ω; [kx, ky, kz, lupack=$(DEFAULT_LUPACK)]) -> S

Frequency `ω` is a (real) number, vector, or range.
Scattering matrix `S` is N-by-N-by-M, where M=`length(ω)` and N is the number
of channels.

Automatically parallelizes over frequency `ω` if `nprocs()>1`.

Noteworthy: `S` is easily plotted with Plots.jl: `plot(S)`.
"""
smatrix(sim::Simulation, ω::Real; kwargs...) = smatrix(sim, [ω]; kwargs...)

# core smatrix computation taking place during SharedMatrix initialization
function (spc::SMatrixContainer)(S::SharedArray)
    if indexpids(S)==0
        return nothing
    else
        idx = indexpids(S)
        nchunks = length(procs(S))
        splits = [round(Int, s) for s in range(0, stop=size(S,3), length=nchunks+1)]
        indices = splits[idx]+1:splits[idx+1]
    end
    dimensional_smatrix!(S,spc,indices)
    return nothing
end

################################################################################
# Pretty Printing
import ..Common.PRINTED_COLOR_LIGHT
import ..Common.PRINTED_COLOR_NUMBER

function Base.show(io::IO,s::ScatteringMatrix{N}) where N
    printstyled(io,"ScatteringMatrix ",color=PRINTED_COLOR_LIGHT)
    print(io,"($N-by-$N) @ $(length(s.ω)) frequencies")
end

################################################################################
# Plotting

# first smooth version when ω is LinRange
@recipe f(s::ScatteringMatrix{N,TW};n=1001) where {N,TW<:LinRange} = s,n
@recipe function f(s::ScatteringMatrix{1,TW},n::Integer) where TW<:LinRange
    ωs = LinRange(s.ω[1],s.ω[end],n)
    @series begin
        label := L"|r|^2"
        RL = CubicSplineInterpolation(s.ω,abs2.(s.S[1,1,:]))
        ωs, RL(ωs)
    end
end
@recipe function f(s::ScatteringMatrix{2,TW},n::Integer) where TW<:LinRange
    ωs = LinRange(s.ω[1],s.ω[end],n)
    @series begin
        label := L"|t|^2"
        T = CubicSplineInterpolation(s.ω,abs2.(s.S[1,2,:]))
        ωs, T(ωs)
    end
    @series begin
        label := L"|r_L|^2"
        RL = CubicSplineInterpolation(s.ω,abs2.(s.S[1,1,:]))
        ωs, RL(ωs)
    end
    @series begin
        label := L"|r_R|^2"
        RR = CubicSplineInterpolation(s.ω,abs2.(s.S[2,2,:]))
        ωs, RR(ωs)
    end
end

# next jagged version when ω is Array
@recipe function f(s::ScatteringMatrix{1,Vector})
    @series begin
        label := L"|r|^2"
        s.ω, abs2.(s.S[1,1,:])
    end
end
@recipe function f(s::ScatteringMatrix{2,Vector})
    @series begin
        label := L"|t|^2"
        s.ω, abs2.(s.S[1,2,:])
    end
    @series begin
        label := L"|r_L|^2"
        s.ω, abs2.(s.S[1,1,:])
    end
    @series begin
        label := L"|r_R|^2"
        s.ω, abs2.(s.S[2,2,:])
    end
end

end # module


text:
"""
    CubicPeriodicSemidiscretisation{T,Derivative,Dissipation}

A semidiscretisation of the cubic conservation law
    \$\\partial_t u(t,x) + \\partial_x u(t,x)^3 = 0\$
with periodic boundary conditions.
"""
struct CubicPeriodicSemidiscretisation{T,Derivative<:AbstractDerivativeOperator{T},
                                        Dissipation,
                                        SplitForm<:Union{Val{false}, Val{true}}} <: AbstractSemidiscretisation
    derivative::Derivative
    dissipation::Dissipation
    tmp1::Vector{T}
    tmp2::Vector{T}
    split_form::SplitForm

    function CubicPeriodicSemidiscretisation(derivative::Derivative, dissipation::Dissipation, split_form::SplitForm=Val{false}()) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}}
        if dissipation != nothing
            @argcheck size(derivative) == size(dissipation) DimensionMismatch
            @argcheck grid(derivative) == grid(dissipation) ArgumentError
        end
        N = size(derivative, 2)
        tmp1 = Array{T}(undef, N)
        tmp2 = Array{T}(undef, N)
        new{T,Derivative,Dissipation,SplitForm}(derivative, dissipation, tmp1, tmp2, split_form)
    end
end

function Base.show(io::IO, semidisc::CubicPeriodicSemidiscretisation)
    print(io, "Semidiscretisation of the cubic conservation law\n")
    print(io, "  \$ \\partial_t u(t,x) + \\partial_x u(t,x)^3 = 0 \$ \n")
    print(io, "with periodic boundaries using")
    if semidisc.split_form == Val{true}()
        print(io, " a split form and: \n")
    else
        print(io, " no split form and: \n")
    end
    print(io, semidisc.derivative)
    print(io, semidisc.dissipation)
end

function (disc::CubicPeriodicSemidiscretisation)(du, u, p, t)
    @unpack tmp1, tmp2, derivative, dissipation, split_form = disc
    @boundscheck begin
        @argcheck length(u) == length(tmp1)
        @argcheck length(u) == length(tmp2)
        @argcheck length(u) == length(du)
    end

    # volume terms
    if typeof(split_form) <: Val{true}
        mhalf = -one(eltype(u)) / 2

        ## u^2 * D * u
        @. tmp2 = u^2
        mul!(tmp1, derivative, u)
        @. du = mhalf * tmp2 * tmp1
        ## u * D * u^2
        mul!(tmp1, derivative, tmp2)
        @. du += mhalf * u * tmp1
        ## D * u^3
        @. tmp2 *= u
        mul!(tmp1, derivative, tmp2)
        @. du += mhalf * tmp1
    else
        @. tmp2 = u^3
        mul!(tmp1, derivative, tmp2)
        @. du = -tmp1
    end

    # dissipation
    if dissipation != nothing
        mul!(tmp1, dissipation, u)
        @. du += tmp1
    end

    nothing
end

"""
    CubicNonperiodicSemidiscretisation

A semidiscretisation of the cubic conservation law
    \$\\partial_t u(t,x) + \\partial_x u(t,x)^3 = 0\$
with nonperiodic boundary conditions.
"""
struct CubicNonperiodicSemidiscretisation{T,Derivative<:AbstractDerivativeOperator{T},
                                            Dissipation,
                                            SplitForm<:Union{Val{false}, Val{true}},
                                            LeftBC, RightBC} <: AbstractSemidiscretisation
    derivative::Derivative
    dissipation::Dissipation
    tmp1::Vector{T}
    tmp2::Vector{T}
    split_form::SplitForm
    left_bc::LeftBC
    right_bc::RightBC

    function CubicNonperiodicSemidiscretisation(derivative::Derivative, dissipation::Dissipation, split_form::SplitForm, left_bc::LeftBC, right_bc::RightBC) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}, LeftBC, RightBC}
        if dissipation != nothing
            @argcheck size(derivative) == size(dissipation) DimensionMismatch
            @argcheck grid(derivative) == grid(dissipation) ArgumentError
        end
        N = size(derivative, 2)
        tmp1 = Array{T}(undef, N)
        tmp2 = Array{T}(undef, N)
        new{T,Derivative,Dissipation,SplitForm,LeftBC,RightBC}(derivative, dissipation, tmp1, tmp2, split_form, left_bc, right_bc)
    end
end

function Base.show(io::IO, semidisc::CubicNonperiodicSemidiscretisation)
    print(io, "Semidiscretisation of the cubic conservation law\n")
    print(io, "  \$ \\partial_t u(t,x) + \\partial_x u(t,x)^3 = 0 \$ \n")
    print(io, "with nonperiodic boundaries using")
    if semidisc.split_form == Val{true}()
        print(io, " a split form and: \n")
    else
        print(io, " no split form and: \n")
    end
    print(io, semidisc.derivative)
    print(io, semidisc.dissipation)
end

function godunov_flux_cubic(uₗ::T, uᵣ::T) where {T<:Real}
    uₗ^3
end

function (disc::CubicNonperiodicSemidiscretisation)(du, u, p, t)
    @unpack tmp1, tmp2, derivative, dissipation, split_form, left_bc, right_bc = disc
    @boundscheck begin
        @argcheck length(u) == length(tmp1)
        @argcheck length(u) == length(tmp2)
        @argcheck length(u) == length(du)
    end

    # volume terms
    if typeof(split_form) <: Val{true}
        mhalf = -one(eltype(u)) / 2

        ## u^2 * D * u
        @. tmp2 = u^2
        mul!(tmp1, derivative, u)
        @. du = mhalf * tmp2 * tmp1
        ## u * D * u^2
        mul!(tmp1, derivative, tmp2)
        @. du += mhalf * u * tmp1
        ## D * u^3
        @. tmp2 *= u
        mul!(tmp1, derivative, tmp2)
        @. du += mhalf * tmp1
    else
        @. tmp2 = u^3
        mul!(tmp1, derivative, tmp2)
        @. du = -tmp1
    end

    # dissipation
    if dissipation != nothing
        mul!(tmp1, dissipation, u)
        @. du += tmp1
    end

    # boundary conditions via Godunov's flux
    @inbounds fnum_left = godunov_flux_cubic(left_bc(t), u[1])
    @inbounds du[1] += (fnum_left - u[1]^3) / left_boundary_weight(derivative)
    @inbounds fnum_right = godunov_flux_cubic(u[end], right_bc(t))
    @inbounds du[end] -= (fnum_right - u[end]^3) / right_boundary_weight(derivative)

    nothing
end


text:
struct Habitat{T, N}
    # Store habitats as a Tuple because there probably won't be very many
    # habitats in most cases, and it allows disparate element types (e.g. some
    # Float64 and some Bool).
    habs::T

    function Habitat(habs::T) where {T<:Tuple}
        N = length(habs)
        chksz = size(habs[1])
        for i in 2:N
            _check_habsize(habs[i], chksz) || error("Habitat $i size does not match")
        end
        new{T, N}(habs)
    end
end

function _check_habsize(hab, chksz)
    sz = size(hab)
    all(sz .== chksz)
end

function Habitat(habs::Vector{<:Matrix})
    Habitat(tuple(habs...))
end

function Habitat(habs...)
    Habitat(tuple(habs...))
end

length(hab::Habitat{T, N}) where {T, N} = N
size(hab::Habitat) = size(hab.habs[1])

getindex(hab::Habitat, i) = getindex(hab.habs, i)

struct HabitatPreference{F, N}
    pref::F

    function HabitatPreference(pref::F) where F<:Tuple
        N = length(pref)
        new{F, N}(pref)
    end
end

function HabitatPreference(prefs...)
    pref = tuple(prefs...)
    HabitatPreference(pref)
end

function HabitatPreference(prefvec::V) where V<:Vector
    pref = tuple(prefvec...)
    HabitatPreference(pref)
end

function (habpref::HabitatPreference{F, N} where F)(hab::Habitat{T, N} where T; normalize = true) where {N}
    pref = mapreduce(.*, 1:N) do n
        habpref.pref[n].(hab[n])
    end
    if (normalize)
        pref ./= sum(pref)
    end
    pref
end