text:
"""
    RandomDotOrg

Use various functions from https://random.org (q.v.).

# Currently available
- getQuota(): obtain current bit quota for your IP.
- checkQuota(): check if quota is non-zero.
- random_numbers(): obtain integers.
- random_sequences(): obtain randomized sequences of integers 1..N
- random_strings(): obtain random strings of characters (upper/lower case letters, digits)
- random_gaussian(): obtain random Gaussian numbers
- random_decimal_fractions(): obtain random numbers on the interval (0,1)
- random_bytes(): obtain random bytes in various formats
- random_bitmap(): obtain a random bitmap of size up to 300 x 300 as GIF or PNG.

Github repository at: https://github.com/edwcarney/RandomDotOrg

"""
module RandomDotOrg

using HTTP, Printf
export  getQuota, checkQuota, random_numbers, random_sequences, random_strings, random_gaussian,
        random_decimal_fractions, random_bytes, random_bitmap

"""
    Get the current bit quota from Random.org
"""
function getQuota()
    r = HTTP.get("https://www.random.org/quota/?format=plain");
    return parse(Int64, rstrip(String(r.body)))
end;

"""

    checkQuota(minimum = 500::Number)

    Test for sufficient quota to insure response. This should be set to match
    user's needs.
"""
function checkQuota(minimum = 500::Number)
    return getQuota() >= minimum
end;

"""

    random_numbers(n = 100::Number, min = 1, max = 20, base = 10, parse=true, check = true, col = 5)

Get `n` random integers on the interval `[min, max]` as strings
in one of 4 `base` values--[binary (2), octal (8), decimal (10), or hexadecimal (16)]
All numbers are returned as strings (as Random.org sends them).

# Arguments
- `max`,`min` : [-1e9, 1e9]
- `base::Integer`: retrieved Integer format [2, 8, 10, 16]
- `numeric`: return numbers instead of strings (base = 10, only)
- `check::Bool`: perform a call to `checkQuota` before making request
- `col::Integer`: used to fulfill parameter requirments of random.org

# Examples
```
julia> random_numbers(5, max=50, base=16)
5-element Array{SubString{String},1}:
 "1f"
 "27"
 "12"
 "05"
 "2c"

julia> random_numbers(5, max=4096, base=2, min=2048)
5-element Array{SubString{String},1}:
 "0110011101101"
 "0110101001011"
 "0100101001110"
 "0101000001110"
 "0101101101100"

julia> random_numbers(5, max=200)
5-element Array{Int64,1}:
  69
 107
  69
 155
  72

julia> random_numbers(5, max=200, numeric=false)
5-element Array{SubString{String},1}:
 "121"
 "178"
 "155"
 "12"
 "71"
```
"""
function random_numbers(n = 100::Number; min = 1, max = 20, base = 10, numeric = true, check = true, col = 5)
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (min < -1f+09 || max > 1f+09 || min > max) && return "Range must be between -1e9 and 1e9"

    (!(base in [2, 8, 10, 16])) && return "Base has to be one of 2, 8, 10 or 16"

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    urlbase = "https://www.random.org/integers/"
    urltxt = @sprintf("%s?num=%d&min=%d&max=%d&col=%d&base=%d&format=plain&rnd=new",
                    urlbase, n, Int(min), Int(max), col, base)
    response = HTTP.get(urltxt)
    (numeric && base == 10) && return [parse(Int, x) for x in split(rstrip(String(response.body)))]
    return split(rstrip(String(response.body)))
end;

"""
Get a randomized interval `[min, max]`. Returns (max - min + 1) randomized integers
All numbers are returned as strings (as Random.org sends them).

# Arguments
- `min` : no less than 1
- `max` : must be [-1e9, 1e9]
- `col::Integer` : used to fulfill parameter requirments of random.org
- `check::Bool`: perform a call to `checkQuota` before making request

# Example
```
julia> random_sequences(max=8)
8-element Array{Int64,1}:
 4
 1
 5
 3
 6
 7
 8
 2
```
"""
function random_sequences(;min = 1::Number, max = 20::Number, col = 1, check = true)
    (min < -1f+09 || max > 1f+09 || min > max) && return "Range must be between -1e9 and 1e9"

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    urlbase = "https://www.random.org/sequences/"
    urltxt = @sprintf("%s?min=%d&max=%d&col=%d&format=plain&rnd=new",
                    urlbase, Int(min), Int(max), col)
    response = HTTP.get( urltxt)
    return [parse(Int64, x) for x in split(rstrip(String(response.body)))]
end;

"""

    random_strings(n=10::Number, len=5; digits=true, upperalpha=true, loweralpha=true, unique=true, check=true)

Get `n` random strings of length `len`

# Arguments
- `len` : length of strings [1, 20]
- `digits::Bool`: include digits in strings
- `upperalpha::Bool`: include upper case letters in strings
- `loweralpha::Bool`: include lowercase letters in strings
- `unique::Bool`: strings must be unique if `true`
- `check::Bool`: perform a call to `checkQuota` before making request

# Examples
```
julia> random_strings(5, 8, loweralpha=false)
5-element Array{SubString{String},1}:
 "8LPPSHYL"
 "Y3I7ILXJ"
 "0O9C41EK"
 "Y2N7U9WZ"
 "OH9EZ3ST"

julia> random_strings(3, 4, digits=false)
3-element Array{SubString{String},1}:
 "cKIG"
 "HUir"
 "gAzX"
```
"""
function random_strings(n=10::Number, len=5; digits=true, upperalpha=true, loweralpha=true, unique=true, check=true)
    (n < 1 || n > 10000) && return "1 to 10,000 requests only"

    (len < 1 || len > 20) && return "Length must be between 1 and 20"

    (typeof(digits) != Bool || typeof(upperalpha) != Bool || typeof(loweralpha) != Bool || typeof(unique) != Bool) && return "The 'digits', '(lower|upper)alpha' and 'unique' arguments have to be logical"

    (!digits && !upperalpha && !loweralpha) && return "The 'digits', 'loweralpha' and 'upperalpha' cannot all be false"

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    urlbase = "https://www.random.org/strings/"
    urltxt = @sprintf("%s?num=%d&len=%d&digits=%s&upperalpha=%s&loweralpha=%s&unique=%s&format=plain&rnd=new",
                urlbase, n, len, ifelse(digits, "on", "off"),
                ifelse(upperalpha, "on", "off"), ifelse(loweralpha, "on", "off"),
                ifelse(unique, "on", "off"))
    response = HTTP.get( urltxt)
    split(rstrip(String(response.body)))
end;

"""

    random_gaussian(n=10::Number, mean=0.0, stdev=1.0; dec=10, col=1, notation="scientific", check=true)

Get n numbers from a Gaussian distribution with `mean` and `stdev`.
Returns strings in `dec` decimal places.
Scientific notation only for now.

# Arguments
- `mean`, `stdev` : between [-1e6, 1e6]
- `dec` : decimal places [2,20]
- `col` : unique value
- `scientific` : unique value
- `check::Bool`: perform a call to `checkQuota` before making request

# Example
```
julia> random_gaussian(5, dec=5)
5-element Array{Float64,1}:
 -0.36174
  0.95992
  1.2278
  1.0429
  0.31885
```
"""
function random_gaussian(n=10::Number, mean=0.0, stdev=1.0; dec=10, col=1, notation="scientific", check=true)
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (mean < -1f+06 || mean > 1f+06) && return "Mean must be between -1e6 and 1e6"

    (stdev < -1f+06 || stdev > 1f+06) && return "Std dev must be between -1e6 and 1e6"

    (dec < 2 || dec > 20) && return "Decimal places must be between 2 and 20"

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    urlbase = "https://www.random.org/gaussian-distributions/"
    urltxt = @sprintf("%s?num=%d&mean=%f&stdev=%f&dec=%d&col=%d&notation=%s&format=plain&rnd=new",
                    urlbase, n, mean, stdev, dec, col, notation)
    # print(urltxt)
    response = HTTP.get( urltxt)
    return [parse(Float64, x) for x in split(rstrip(String(response.body)))]
end;

"""
    random_decimal_fractions(n=10::Number; dec=10, col=2, check=true)

Get n decimal fractions on the interval (0,1).
Returns strings in `dec` decimal places.

# Arguments
- `dec` : decimal places [2, 20]
- `col` : unique value
- `check::Bool`: perform a call to `checkQuota` before making request

# Example
```
julia> random_decimal_fractions(5, dec=3)
5-element Array{Float64,1}:
 0.788
 0.949
 0.773
 0.875
 0.552
```
"""
function random_decimal_fractions(n=10::Number; dec=10, col=2, check=true)
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (dec < 2 || dec > 20) && return "Decimal places must be between 2 and 20"

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    urlbase = "https://www.random.org/decimal-fractions/"
    urltxt = @sprintf("%s?num=%d&dec=%d&col=%d&format=plain&rnd=new",
                        urlbase, n, dec, col)
    # print(urltxt)
    response = HTTP.get( urltxt)
    return [parse(Float64, x) for x in split(rstrip(String(response.body)))]
end;

"""
    random_bytes(n=10::Number; format="o", check=true)

Get n random bytes.
Returns strings in binary, decimal, octal, or hexadecimal.
The request will also download a DMS file.

# Arguments
- `format` : format of bytes; one of b, d, o, h, or file (Int in hex)
- `check::Bool`: perform a call to `checkQuota` before making request

# Examples
```
julia> random_bytes(5, format='o')
5-element Array{SubString{String},1}:
 "075"
 "030"
 "317"
 "277"
 "356

julia> random_bytes(5, format="file")
5-element Array{UInt8,1}:
 0xda
 0xdf
 0xd9
 0xb3
 0x97
```
"""
function random_bytes(n=10::Number; format='o', check=true)
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (!(format in ['b', 'd', 'o', 'h', "file"])) && return "Format must be one of b, d, o, h, or file."

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    urlbase = "https://www.random.org/cgi-bin/randbyte"
    urltxt = @sprintf("%s?nbytes=%d&format=%s",
                        urlbase, n, format)
    # print(urltxt)
    # response = HTTP.get( urltxt)
    response = HTTP.get(urltxt)
    if (format == "file")
        return response.body
    else
        return split(rstrip(String(response.body)))
    end
end

"""
    random_bitmap([format="png"; width=64, height=64, save="", overwrite='n', check=true])

Get a random bitmap as a PNG or GIF.

# Arguments
- `format::String`: png or gif
- `height, width`: 1-300 pixels in each dimension
- `save`: filename/filepath to write file (extension added by default)
- `overwrite`: will not overwrite by default; set to 'y' to overwrite
- `check::Bool`: perform a call to `checkQuota` before making request

"""
function random_bitmap(format="png"; width=64, height=64, save="", overwrite='n', check=true)
    (!(format in ["png", "gif"])) && return "Format must be png or gif."

    (width < 1 || width > 300 || height < 1 || height > 300) && return "Height/Width must be no more than 300."

    (check && !checkQuota()) && return "random.org suggests to wait until tomorrow"

    full_path = @sprintf("%s.%s", save, format)
    (isfile(full_path) && overwrite=='n') && return @sprintf("File %s exists. Set overwrite to 'y' to save.", full_path)

    urlbase = "https://www.random.org/bitmaps"
    urltxt = @sprintf("%s?format=%s&height=%d&width=%d&zoom=1",
                        urlbase, format, height, width)
    response = HTTP.get(urltxt)
    if save == ""
        return(response.body)
    else
        # full_path = @sprintf("%s.%s", save, format)
        # if isfile(full_path) && overwrite == 'n'
        #     @printf("File not saved; %s exists. Set overwrite to \"y\".", full_path)
        # else
        open(full_path, "w") do outfile
            write(outfile, response.body)
        end
        @printf("File saved as: %s\n", full_path)
    # end
    end
end

end;  # module


text:
export precom

# precompiling information
struct PreInfo
    routines::Vector{Any}
end
PreInfo() = PreInfo([])

"""
    precom(module, ex)

Precompile a function, returns a tuple of (macros, function name, arguments, type parameters, function body).
"""
function precom(m::Module, ex)
    mc, fname, args, ts, body = match_function(ex)
    vars = Symbol[]
    newargs = map(args) do arg
        @smatch arg begin
            :(::$tp)=>Expr(:(::), gensym(), tp)
            _ => arg
        end
    end
    for arg in newargs
        pushvar!(vars, arg)
    end
    info = PreInfo()
    body_out = precom_body(m, body, info)
    if !isempty(info.routines)
        error("`@routine` and `~@routine` must appear in pairs, mising `~@routine`!")
    end
    st = SymbolTable(vars, Symbol[], Symbol[])
    st_after = copy(st)
    variable_analysis_ex.(body_out, Ref(st_after))
    checksyms(st_after, st)
    mc, fname, newargs, ts, body_out
end

function precom_body(m::Module, body::AbstractVector, info)
    Any[precom_ex(m, ex, info) for ex in body]
end

# precompile `+=`, `-=`, `*=` and `/=`
function precom_opm(f, out, arg2)
    if f in [:(+=), :(-=), :(*=), :(/=)]
        @smatch arg2 begin
            :($x |> $view) => Expr(f, out, :(identity($arg2)))
            :($subf($(subargs...))) => Expr(f, out, arg2)
            _ => Expr(f, out, :(identity($arg2)))
        end
    elseif f in [:(.+=), :(.-=), :(.*=), :(./=)]
        @smatch arg2 begin
            :($x |> $view) || :($x .|> $view) => Expr(f, out, :(identity.($arg2)))
            :($subf.($(subargs...))) => Expr(f, out, arg2)
            :($subf($(subargs...))) => Expr(f, out, arg2)
            _ => Expr(f, out, :(identity.($arg2)))
        end
    end
end

# precompile `⊻=`
function precom_ox(f, out, arg2)
    if f == :(⊻=)
        @smatch arg2 begin
            :($x |> $view) => Expr(f, out, :(identity($arg2)))
            :($subf($(subargs...))) ||
            :($a || $b) || :($a && $b) => Expr(f, out, arg2)
            _ => Expr(f, out, :(identity($arg2)))
        end
    elseif f == :(.⊻=)
        @smatch arg2 begin
            :($x |> $view) || :($x .|> $view) => Expr(f, out, :(identity.($arg2)))
            :($subf.($(subargs...))) => Expr(f, out, arg2)
            :($subf($(subargs...))) => Expr(f, out, arg2)
            _ => Expr(f, out, :(identity.($arg2)))
        end
    end
end

"""
    precom_ex(module, ex, info)

Precompile a single statement `ex`, where `info` is a `PreInfo` instance.
"""
function precom_ex(m::Module, ex, info)
    @smatch ex begin
        :($x ← $val) || :($x → $val) => ex
        :($x ↔ $y) => ex
        :($(xs...), $y ← $val) => precom_ex(m, :(($(xs...), $y) ← $val), info)
        :($(xs...), $y → $val) => precom_ex(m, :(($(xs...), $y) → $val), info)
        :($a += $b) => precom_opm(:+=, a, b)
        :($a -= $b) => precom_opm(:-=, a, b)
        :($a *= $b) => precom_opm(:*=, a, b)
        :($a /= $b) => precom_opm(:/=, a, b)
        :($a ⊻= $b) => precom_ox(:⊻=, a, b)
        :($a .+= $b) => precom_opm(:.+=, a, b)
        :($a .-= $b) => precom_opm(:.-=, a, b)
        :($a .*= $b) => precom_opm(:.*=, a, b)
        :($a ./= $b) => precom_opm(:./=, a, b)
        :($a .⊻= $b) => precom_ox(:.⊻=, a, b)
        Expr(:if, _...) => precom_if(m, copy(ex), info)
        :(while ($pre, $post); $(body...); end) => begin
            post = post == :~ ? pre : post
            info = PreInfo()
            Expr(:while, :(($pre, $post)), Expr(:block, precom_body(m, body, info)...))
        end
        :(@from $line $post while $pre; $(body...); end) => precom_ex(m, Expr(:while, :(($pre, !$post)), ex.args[4].args[2]), info)
        :(begin $(body...) end) => begin
            Expr(:block, precom_body(m, body, info)...)
        end
        # TODO: allow ommit step.
        :(for $i=$range; $(body...); end) ||
        :(for $i in $range; $(body...); end) => begin
            info = PreInfo()
            Expr(:for, :($i=$(precom_range(range))), Expr(:block, precom_body(m, body, info)...))
        end
        :(@safe $line $subex) => ex
        :(@cuda $line $(args...)) => ex
        :(@launchkernel $line $(args...)) => ex
        :(@inbounds $line $subex) => Expr(:macrocall, Symbol("@inbounds"), line, precom_ex(m, subex, info))
        :(@simd $line $subex) => Expr(:macrocall, Symbol("@simd"), line, precom_ex(m, subex, info))
        :(@threads $line $subex) => Expr(:macrocall, Symbol("@threads"), line, precom_ex(m, subex, info))
        :(@avx $line $subex) => Expr(:macrocall, Symbol("@avx"), line, precom_ex(m, subex, info))
        :(@invcheckoff $line $subex) => Expr(:macrocall, Symbol("@invcheckoff"), line, precom_ex(m, subex, info))
        :(@routine $line $expr) => begin
            precode = precom_ex(m, expr, info)
            push!(info.routines, precode)
            precode
        end
        :(~(@routine $line)) => begin
            if isempty(info.routines)
                error("`@routine` and `~@routine` must appear in pairs, mising `@routine`!")
            end
            precom_ex(m, dual_ex(m, pop!(info.routines)), info)
        end
        # 1. precompile to expand macros
        # 2. get dual expression
        # 3. precompile to analyze vaiables
        :(~$expr) => precom_ex(m, dual_ex(m, precom_ex(m, expr, PreInfo())), info)
        :($f($(args...))) => :($f($(args...)))
        :($f.($(args...))) => :($f.($(args...)))
        :(nothing) => ex
        Expr(:macrocall, _...) => precom_ex(m, macroexpand(m, ex), info)
        ::LineNumberNode => ex
        ::Nothing => ex
        _ => error("unsupported statement: $ex")
    end
end

precom_range(range) = @smatch range begin
    _ => range
end

function precom_if(m, ex, exinfo)
    _expand_cond(cond) = @smatch cond begin
        :(($pre, ~)) => :(($pre, $pre))
        :(($pre, $post)) => :(($pre, $post))
        :($pre) => :(($pre, $pre))
    end
    if ex.head == :if
        ex.args[1] = _expand_cond(ex.args[1])
    elseif ex.head == :elseif
        ex.args[1].args[2] = _expand_cond(ex.args[1].args[2])
    end
    info = PreInfo()
    ex.args[2] = Expr(:block, precom_body(m, ex.args[2].args, info)...)
    if length(ex.args) == 3
        if ex.args[3].head == :elseif
            ex.args[3] = precom_if(m, ex.args[3], exinfo)
        elseif ex.args[3].head == :block
            info = PreInfo()
            ex.args[3] = Expr(:block, precom_body(m, ex.args[3].args, info)...)
        else
            error("unknown statement following `if` $ex.")
        end
    end
    ex
end

export @code_preprocess

"""
    @code_preprocess ex

Preprocess `ex` and return the symmetric reversible IR.

```jldoctest; setup=:(using NiLangCore)
julia> NiLangCore.rmlines(@code_preprocess if (x < 3, ~) x += exp(3.0) end)
:(if (x < 3, x < 3)
      x += exp(3.0)
  end)
```
"""
macro code_preprocess(ex)
    QuoteNode(precom_ex(__module__, ex, PreInfo()))
end

precom_ex(m::Module, ex) = precom_ex(m, ex, PreInfo())

# push a new variable to variable set `x`, for allocating `target`
function pushvar!(x::Vector{Symbol}, target)
    @smatch target begin
        ::Symbol => begin
            if target in x
                throw(InvertibilityError("Symbol `$target` should not be used as the allocation target, it is an existing variable in the current scope."))
            else
                push!(x, target)
            end
        end
        :(($(tar...),)) => begin
            for t in tar
                pushvar!(x, t)
            end
        end
        :($tar = _) => pushvar!(x, tar)
        :($tar...) => pushvar!(x, tar)
        :($tar::$tp) => pushvar!(x, tar)
        Expr(:parameters, targets...) => begin
            for tar in targets
                pushvar!(x, tar)
            end
        end
        Expr(:kw, tar, val) => begin
            pushvar!(x, tar)
        end
        _ => error("unknown variable expression $(target)")
    end
    nothing
end

ext: jl


text:
module samplers

using Distributions, LinearAlgebra
using StatsBase
using Printf
using CSV
using DataFrames
using DataStructures

include("outFiles.jl")
include("misc.jl")

export runSampler

#main sampler
function runSampler(rowID,A,Y,X,Z,chainLength,burnIn,outputFreq,priorVCV,M,paths2maps,rS)

	println("priorVCV $(priorVCV)")
	#output settings
	these2Keep  = collect((burnIn+outputFreq):outputFreq:chainLength) #print these iterations

        #some info
	##This is not really nFix, but the "blocks" of fixed effects
        nFix  = length(X)
	nRand = length(Z)
	nData = length(Y)
	nMarkerSets = length(M)
	nMarkers    = [size(M[m],2) for m in keys(M)]
	println("nMarkers: $nMarkers")

        #initial computations and settings
	ycorr = deepcopy(Y)

	##make iXpX, Z', zpz (for uncor)
        iXpX = deepcopy(X)
        for x in keys(X)
                iXpX[x] = inv(X[x]'X[x])
        end

	Zp  = deepcopy(Z) #similar(Z')
#	zpz = Array{Array{Float64, 1},1}(undef,0)
	zpz = OrderedDict{Any,Any}()
	for z in keys(Z)
                zpz[z] = diag(Z[z]'Z[z])
		Zp[z]  = Z[z]'
        end

        #key positions for speed
        XKeyPos = OrderedDict{Any,Int64}()
        [XKeyPos[collect(keys(X))[i]]=i for i in 1:length(keys(X))]
	println("XKeyPos: $XKeyPos")

	ZKeyPos = OrderedDict{Any,Int64}()
	[ZKeyPos[collect(keys(Z))[i]]=i for i in 1:length(keys(Z))]
        println("ZKeyPos: $ZKeyPos")

        ##make b and u arrays
        b = Array{Array{Float64, 1},1}(undef,0)
        ##counts columns per effect
        nColEachX = []
        for xSet in keys(X)
                println(xSet)
                nCol = size(X[xSet],2)
                push!(b,fill(0.0,nCol))
                nColEachX = push!(nColEachX,nCol)
        end

        u = Array{Array{Float64, 1},1}(undef,0)
        ##counts columns per effect
        nColEachZ = []
	##get priors per effect
	iVarStr = Dict{Any,Array{Float64,2}}() #inverses will be computed
	varU_prior = OrderedDict{Any,Any}()
        for zSet in keys(Z)
                nCol = size(Z[zSet],2)
                push!(u,fill(0.0,nCol))
                nColEachZ = push!(nColEachZ,nCol)
		#var structures and priors
		if haskey(priorVCV,zSet)
			if isempty(priorVCV[zSet][1])
				println("priorVCV structure for $zSet is empty, an identity matrix will be used")
				iVarStr[zSet] = Matrix(1.0I,nCol,nCol)
			elseif priorVCV[zSet][1]=="A"
				iVarStr[zSet] = inv(A)
				println("priorVCV structure for $zSet is A, computed A matrix will be used")
			else 	iVarStr[zSet] = inv(priorVCV[zSet][1])
			end
			varU_prior[zSet] = priorVCV[zSet][2]
		else	println("priorVCV structure for $zSet is empty, an identity matrix will be used with an arbitrary variance of 100")
			iVarStr[zSet] = Matrix(1.0I,nCol,nCol)
			varU_prior[zSet] = 100
		end
        end

	#set up for E
	isempty(priorVCV["e"][1]) ? strE = Matrix(1.0I,nData,nData) : strE = priorVCV["e"][1]
	varE_prior = priorVCV["e"][2]

	#parameters for priors
        dfE = 4.0
	dfDefault = 4.0

	if varE_prior==0.0
		varE_prior  = 0.0005
       		scaleE     = 0.0005
        else
       		scaleE    = varE_prior*(dfE-2.0)/dfE
   	end

	##no correlated random effects
	scaleU = Dict{Any,Any}()
	for zSet in keys(Z)
		scaleU[zSet] = varU_prior[zSet]*(dfDefault-2.0)/dfDefault
	end

	#pre-computations using priors, not relevant for correlated random effects
   	νS_E = scaleE*dfE

	#ADD MARKERS
	# read map file and make regions

	#key positions for each effect in beta, for speed. Order of matrices in M are preserved here.
        BetaKeyPos = OrderedDict{String,Int64}()
        println("keysM: $(keys(M))")
        for mSet in keys(M)
                pos = findall(mSet.==collect(keys(M)))[]
                BetaKeyPos[mSet] = pos
        end
        println("BetaKeyPos: $BetaKeyPos")

	#make mpm

	Mp = OrderedDict{Any,Any}()
       	mpm = OrderedDict{Any,Any}()

	corM = OrderedDict{Any,Any}()
	corMPos = OrderedDict{Any,Any}()
	regionArray = OrderedDict{Any,Array{UnitRange{Int64},1}}()

	for pSet ∈ keys(priorVCV)
		corEffects = []
		corPositions = []
		if typeof(pSet)==String
			println("$pSet is univariate")
			if pSet ∈ keys(M)
				println("univariate mpm for $pSet")
				tempmpm = []
				nowM = M[pSet]
				for c in eachcol(nowM)
					push!(tempmpm,BLAS.dot(c,c))
				end
				mpm[pSet] = tempmpm
				theseRegions = prep2RegionData(paths2maps[pSet],rS[pSet])
		                regionArray[pSet] = theseRegions
			end
		else println("$pSet will be correlated")
			correlate = collect(pSet)
			for pSet in correlate
				println(pSet)
				push!(corEffects,pSet)
				push!(corPositions,findall(pSet.==keys(M))[])
			end
			if issubset(corEffects,collect(keys(M)))
				corMPos[pSet] = corPositions
				corM[pSet] = corEffects
				tempM = hcat.(eachcol.(getindex.(Ref(M), (pSet)))...)
				for d in corEffects
                       			delete!(M,d)
               			end
				M[pSet]   = tempM
				mpm[pSet] = MatByMat.(tempM)
				Mp[pSet]  = transpose.(tempM)
				tempM = 0
				nowMap = first(pSet)					 #should throw out error if sets have different lengths! implement it here!
				theseRegions = prep2RegionData(paths2maps[nowMap],rS[nowMap]) ###first data
                		regionArray[pSet] = theseRegions
			end
		end
	end

#	nRegions  = [mSet => length(regionArray[mSet]) for mSet in keys(regionArray)] #per marker set
	nRegions  = OrderedDict(mSet => length(regionArray[mSet]) for mSet in keys(regionArray))
	println("nRegions: $(values(nRegions))")

	println("corM: $corM")
	println("corMPos: $corMPos")

	dfM = Dict{Any,Any}()
	for mSet ∈ keys(mpm)
		dfM[mSet] = 3.0+size(priorVCV[mSet],1)
	end

	println("dfM $dfM")

	scaleM = Dict{Any,Any}()
        for mSet in keys(mpm)
                nMComp = size(priorVCV[mSet],1)
                nMComp > 1 ? scaleM[mSet] = priorVCV[mSet].*(dfM[mSet]-nMComp-1.0)  : scaleM[mSet] = priorVCV[mSet]*(dfM[mSet]-2.0)/dfM[mSet] #I make float and array of float
        end

	println("scaleM $scaleM")

	println("keysM: $(keys(M))")

	#storage

	varU = varU_prior #for storage

	beta = zeros(Float64,nMarkerSets,maximum(vcat([0,nMarkers]...))) #zero is for max to work when no SNP data is present #can allow unequal length! Remove tail zeros for printing....
#	vcovBeta = fill(Matrix(Diagonal(varM)),maximum(nRegions)) #can allow unequal length! Remove tail zeros for printing....

        varBeta = OrderedDict{Any,Any}()
        for mSet in keys(mpm)
                varBeta[mSet] = [priorVCV[mSet] for i in 1:length(regionArray[mSet])] #later, direct reference to key when varM_prior is a dictionary
        end
        println("keys of varBeta: $(keys(varBeta))")
        println("varBeta: $(varBeta)")

	#Start McMC
        for iter in 1:chainLength
		#sample residual variance
	       	varE = sampleVarE(νS_E,ycorr,dfE,nData)

		#sample fixed effects
        	#always returns corrected Y and new b
	        sampleX!(X,b,iXpX,nFix,nColEachX,XKeyPos,ycorr,varE)

		#sample random effects
		# always returns corrected Y and new u
		sampleZ!(iVarStr,Z,Zp,zpz,nRand,ZKeyPos,varE,varU,u,ycorr)

		#sample variances
		sampleRanVar!(varU,nRand,scaleU,dfDefault,u,ZKeyPos,iVarStr)

		#sample marker effects and variances
	        sampleMandMVar_view!(M,Mp,corM,corMPos,beta,mpm,nMarkerSets,BetaKeyPos,regionArray,nRegions,ycorr,varE,varBeta,scaleM,dfM)

        	#print
		if iter in these2Keep
			IO.outMCMC(pwd(),"b",vcat(b...)') ### currently no path is provided!!!!
			IO.outMCMC(pwd(),"u",vcat(u...)')
			IO.outMCMC(pwd(),"varE",varE)
			IO.outMCMC(pwd(),"varU",hcat([varU[k] for k in keys(varU)]...))
			for mSet in keys(BetaKeyPos)
                                IO.outMCMC(pwd(),"beta$mSet",beta[BetaKeyPos[mSet]])
                        end
			for pSet in keys(mpm)
				IO.outMCMC(pwd(),"var".*pSet,varBeta[pSet])
			end
	#		if onScreen==true
            			println("b, $(vcat(b...))") #i always vectorize b. maybe better to make it vector initially
        #		end
		end
	end
end

#Sampling fixed effects

function sampleX!(X,b,iXpX,nFix,nColEachX,keyX,ycorr,varE)
        #block for each effect
        for xSet in keys(X)
		pos = keyX[xSet]
                ycorr    .+= X[xSet]*b[pos]
                rhs      = X[xSet]'*ycorr
                meanMu   = iXpX[xSet]*rhs
                if nColEachX[pos] == 1
                        b[pos] .= rand(Normal(meanMu[],sqrt((iXpX[xSet]*varE))[]))
                else b[pos] .= rand(MvNormal(vec(meanMu),convert(Array,Symmetric(iXpX[xSet]*varE))))
                end
                ycorr    .-= X[xSet]*b[pos]
        end
end

#Sampling random effects
function sampleZ!(iStrMat,Zmat,ZpMat,zpzMat,nRand,ZKeyPos,varE,varU,u,ycorr)
	#block for each effect
	for zSet in keys(Zmat)
		pos = ZKeyPos[zSet]
		uVec = deepcopy(u[pos])
		iMat = iStrMat[zSet]
		tempzpz = zpzMat[zSet] ###added
		λz = varE/(varU[zSet])
	        ycorr .+= Zmat[zSet]*uVec
	        Yi = ZpMat[zSet]*ycorr #computation of Z'ycorr for ALL  rhsU
		nCol = length(zpzMat[zSet])
	        for i in 1:nCol
        	        uVec[i] = 0.0 #also excludes individual from iMat! Nice trick.
			rhsU = Yi[i] - λz*dot(iMat[:,i],uVec)
                	lhsU = tempzpz[i] + (iMat[i,i]*λz)[1]
			invLhsU = 1.0/lhsU
                	meanU = invLhsU*rhsU
                	uVec[i] = rand(Normal(meanU,sqrt(invLhsU*varE)))
        	end
		u[pos] = uVec
        	ycorr .-= Zmat[zSet]*uVec
	end
end

#Sampling marker effects
function sampleM!(MMat,beta,mpmMat,nMSet,keyM,regionsMat,regions,ycorr,varE,varBeta)
        #for each marker set
        for mSet in keys(MMat)
                for r in 1:regions[keyM[mSet]]
                        theseLoci = regionsMat[keyM[mSet]][r]
                        regionSize = length(theseLoci)
                        lambda = varE/(varBeta[keyM[mSet]][r])
                        for locus in theseLoci
                                BLAS.axpy!(beta[keyM[mSet],locus],MMat[mSet][:,locus],ycorr)
                                rhs = BLAS.dot(MMat[mSet][:,locus],ycorr)
                                lhs = mpmMat[mSet][locus] + lambda
                                meanBeta = lhs\rhs
                                beta[keyM[mSet],locus] = sampleBeta(meanBeta, lhs, varE)
                                BLAS.axpy!(-1.0*beta[keyM[mSet],locus],MMat[mSet][:,locus],ycorr)
                        end
                end
        end
end

function sampleMandMVar!(MMat,beta,mpmMat,nMSet,keyM,regionsMat,regions,ycorr,varE,varBeta,scaleM,dfM)
        #for each marker set
        for mSet in keys(MMat)
		pos = keyM[mSet]
                for r in 1:regions[pos]
                        theseLoci = regionsMat[pos][r]
                        regionSize = length(theseLoci)
                        lambda = varE/(varBeta[mSet][r])
                        for locus in theseLoci
                                BLAS.axpy!(beta[pos,locus],MMat[mSet][:,locus],ycorr)
                                rhs = BLAS.dot(MMat[mSet][:,locus],ycorr)
                                lhs = mpmMat[mSet][locus] + lambda
                                meanBeta = lhs\rhs
                                beta[pos,locus] = sampleBeta(meanBeta, lhs, varE)
                                BLAS.axpy!(-1.0*beta[pos,locus],MMat[mSet][:,locus],ycorr)
                        end
			varBeta[mSet][r] = sampleVarBeta(scaleM[pos],dfM[pos],beta[pos,theseLoci],regionSize)
                end
        end
end

function sampleMandMVar_view!(MMat,MpMat,correlatedM,keyCorM,beta,mpmMat,nMSet,keyBeta,regionsMat,regions,ycorr,varE,varBeta,scaleM,dfM)
        #for each marker set
        for mSet in keys(mpmMat)
		if mSet in keys(correlatedM)
			println("$mSet in corM")
			betaPos = keyCorM[mSet]
			nowMp = MpMat[mSet] ###
			for r in 1:regions[mSet]
				theseLoci = regionsMat[mSet][r]
				regionSize = length(theseLoci)
				invB = inv(varBeta[mSet][r])
				for locus in theseLoci
					RHS = zeros(size(invB,1))###
				#	for m in correlatedM[mSet]
				#		BLAS.axpy!(beta[keyBeta[m],locus],MMat[m][:,locus],ycorr) #beta pos is different than pos
				#	end

					ycorr .+= MMat[mSet][locus]*beta[betaPos,locus]

					RHS = (nowMp[locus]*ycorr)./varE ### FASTEST
				#	RHS = zeros(size(invB,1)) ### for mul!
				#	mul!(RHS,nowMp[locus],ycorr./varE) ### LESS MEMORY ALLOCATION
				#	RHS = [BLAS.dot(MMat[m][:,locus],ycorr)/varE for m in correlatedM[mSet]]
					invLHS::Array{Float64,2} = inv((mpmMat[mSet][locus]./varE) .+ invB)
					meanBeta::Array{Float64,1} = invLHS*RHS
					beta[betaPos,locus] = rand(MvNormal(meanBeta,convert(Array,Symmetric(invLHS))))
				#	for m in correlatedM[mSet]
                                #                BLAS.axpy!(-1.0*beta[keyBeta[m],locus],MMat[m][:,locus],ycorr)
                                #       end
					ycorr .-= MMat[mSet][locus]*beta[betaPos,locus]
				end
				varBeta[mSet][r] = sampleVarCovBeta(scaleM[mSet],dfM[mSet],beta[betaPos,theseLoci],regionSize)
			end
		else
			nowM = MMat[mSet]
                	betaPos = keyBeta[mSet]
                	for r in 1:regions[mSet]
                        	theseLoci = regionsMat[mSet][r]
                        	regionSize = length(theseLoci)
                        	lambda = varE/(varBeta[mSet][r])
                        	for locus in theseLoci
                                	BLAS.axpy!(beta[betaPos,locus],view(nowM,:,locus),ycorr)
                                	rhs::Float64 = BLAS.dot(view(nowM,:,locus),ycorr)
                               		lhs::Float64 = mpmMat[mSet][locus] + lambda
                                	meanBeta::Float64 = lhs\rhs
                                	beta[betaPos,locus] = sampleBeta(meanBeta, lhs, varE)
                                	BLAS.axpy!(-1.0*beta[betaPos,locus],view(nowM,:,locus),ycorr)
                       		end
                        	varBeta[mSet][r] = sampleVarBeta(scaleM[mSet],dfM[mSet],beta[betaPos,theseLoci],regionSize)
                	end
       		end
	 end
end

#sample marker effects
function sampleBeta(meanBeta, lhs, varE)
    return rand(Normal(meanBeta,sqrt(lhs\varE)))
end

#sample random effects' variances
function sampleRanVar!(varU,nRand,scale_ranVar,df_ranVar,effVec,keyZ,iStrMat)
	for z in keys(scale_ranVar)
		pos = keyZ[z]
		n = size(iStrMat[z],2)
		varU[z] = (scale_ranVar[z]*df_ranVar + effVec[pos]'*iStrMat[z]*effVec[pos])/rand(Chisq(df_ranVar + n))
	end
end

function sampleMarkerVar!(beta,varBeta,nMSet,keyM,regions,regionsMat,scaleM,dfM)
        #for each marker set
        for mSet in keys(varBeta)
		pos = keyM[mSet]
                for r in 1:regions[pos] #dont have to compute 1000000 times, take it out
                        theseLoci = regionsMat[pos][r]
                        regionSize = length(theseLoci)
                        varBeta[mSet][r] = sampleVarBeta(scaleM[pos],dfM[pos],beta[pos,theseLoci],regionSize)
                end
        end
end

#sample marker variances
function sampleVarBeta(scalem,dfm,whichLoci,regionSize)
	return (scalem*dfm + BLAS.dot(whichLoci,whichLoci)) / rand(Chisq(dfm + regionSize))
end

function sampleVarCovBeta(scalem,dfm,whichLoci,regionSize)
	Sb = whichLoci*whichLoci'
	return rand(InverseWishart(dfm + regionSize, scalem + Sb))
end

#Sample residual variance
function sampleVarE(νS_e,yCorVec,df_e,nRecords)
    return (νS_e + BLAS.dot(yCorVec,yCorVec))/rand(Chisq(df_e + nRecords))
end

end


text:
export @cast, @cast!, @reduce, @reduce!

"""
    @cast Z[i,j,...] := f(A[j,k,...])  options

Macro for broadcasting, reshaping, and slicing of arrays in index notation.
Understands the following things:
* `A[i,j,k]` is a three-tensor with these indices.
* `B[(i,j),k]` is the same thing, reshaped to a matrix. Its first axis (the bracket) is indexed
  by `n = i + (j-1) * N` where `i ∈ 1:N`. This may also be written `B[i\\j,k]` or `B[i⊗j,k]`.
* `C[k][i,j]` is a vector of matrices.
* `D[j,k]{i}` is an ordinary matrix of `SVector`s, which may be reinterpreted from `A[i,j,k]`.
* `E[i,_,k]` has two nontrivial dimensions, and `size(E,2)==1`. On the right hand side
  (or when writing to an existing array) you may also write `E[i,3,k]` meaning `view(E, :,3,:)`,
  or `E[i,\$c,j]` to use a variable `c`. Fixing inner indices, like `C[k][i,_]`, is not allowed.
* `f(x)[i,j,k]` is allowed; `f(x)` must return a 3-tensor (and will be evaluated only once).
* `g(H[:,k])[i,j]` is a generalised `mapslices`, with `g` mapping columns of `H`
    to matrices, which are glued into a 3-tensor `[i,j,k]`.
  Similarly `h(I[k])[i,j]` expects an `h` which maps scalars to matrices.
* `J[i,-j,k]` means roughly `reverse(J, dims=2)`.
* `K[i,j]'` conjugates each element, equivalent to `K'[j,i]` which is the
  conjugate-transpose of the matrix.
* `M[i,i]` means `diag(M)[i]`, but only for matrices: `N[i,i,k]` is an error.

The left and right hand sides must have all the same indices,
with the sole exception of `M[i,i]`.
See `@reduce` and `@mul` for related macros which can sum over things.

If a function of one or more tensors appears on the right hand side,
then this represents a broadcasting operation,
and the necessary re-orientations of axes are automatically inserted.

The following actions are possible:
* `=` writes into an existing array, `copyto!(Z, ...)`.
* `:=` creates a new object. This need not be named: `Z = @cast [i,j,k] := ...` is allowed.
* `==` insists on a view of the old object (error if impossible), and `|=` insists on a copy.
* `=>` creates an anonymous function, for which the output must be on the right.
   For example `@cast A[i] + B[j]' => Z[i,j]` gives `(A,B) -> A .+ B'`.

Re-ordering of indices `Z[k,j,i]` is done lazily with `PermutedDimsArray(A, ...)`.
Reversing of an axis `F[i,-j,k]` is also done lazily, by `Reverse{2}(F)` which makes a `view`.
Using `|=` (or broadcasting) will produce a simple `Array`.

Options can be specified at the end (if several, separated by `,` i.e. `options::Tuple`)
* `i:3` supplies the range of index `i`. Variables and functions like `j:Nj, k:length(K)`
  are allowed.
* `assert` or `!` will turn on explicit dimension checks of the input.
  (Providing ranges may also turn these on, but imperfectly, confirm with `@pretty @cast ...`.)
* `cat` will glue slices by things like `hcat(A...)` instead of the default `reduce(hcat, A)`,
  and `lazy` will instead make a `VectorOfArrays` container.
* `nolazy` disables `PermutedDimsArray` and `Reverse` in favour of `permutedims` and `reverse`,
  and `Diagonal` in favour of `diagm` for `Z[i,i]` output.

These options only work if you load other packages:
* `strided` will place `@strided` in front of broadcasting operations,
  and use `@strided permutedims(A, ...)` instead of `PermutedDimsArray(A, ...)`.
  For this you need `using Strided`.
* Static slices `D[j,k]{i}` need `using StaticArrays`, and to create them you should
  give all slice dimensions explicitly. You may write `D[k]{i:2,j:2}` to specify
  `Size(2,2)` slices. They are made most cleanly from the first indices of the input,
  i.e. this `D` from `A[i,j,k]`.
"""
macro cast(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@cast", exs...))
    _macro(exs...; reduce=false, where=where)
end

"""
    @cast! Z[i...] := A[j...] opt

Variant of `@cast` which effectively runs `@check!()` on each tensor.
"""
macro cast!(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@cast!", exs...))
    _macro(exs...; reduce=false, where=where, icheck=true)
end

"""
    @reduce A[i] := sum(j,k) B[i,j,k]             # A = vec(sum(B, dims=(2,3)))
    @reduce A[i] := prod(j) B[i] + ε * C[i,j]     # A = vec(prod(B .+ ε .* C, dims=2))
    @reduce A[i] = sum(j) exp( C[i,j] / D[j] )    # sum!(A, exp.(C ./ D') )

Tensor reduction macro:
* The reduction function can be anything which works like `sum(B, dims=(1,3))`,
  for instance `prod` and `maximum` and `Statistics.mean`.
* In-place operations `Z[j] = sum(...` will construct the banged version of the given function's name,
  which must work like `sum!(Z, A)`.
* The tensors can be anything that `@cast` understands, including gluing of slices `B[i,k][j]`
  and reshaping `B[i\\j,k]`. See `? @cast` for the complete list.
* If there are several tensors (or scalars) on the right, then this is a broadcasting operation.
* Index ranges may be given afterwards (as for `@cast`) or inside the reduction `sum(i:3, k:4)`.
* All indices appearing on the right must appear either within `sum(...)` etc, or on the left.

    F = @reduce sum(i,j)  B[i] + γ * D[j]         # sum(B .+ γ .* D')
    @reduce G[] := sum(i,j)  B[i] + γ * D[j]      # F == G[]

Complete reduction to a scalar output `F`, or a zero-dim array `G`.
`G[]` involves `sum(A, dims=(1,2))` rather than `sum(A)`.

    @reduce Z[k] := sum(i,j) A[i] * B[j] * C[k]  lazy, i:N, j:N, k:N

The option `lazy` replaces the broadcast expression with a `BroadcastArray`,
to avoid `materialize`ing the entire array before summing. In the example this is of size `N^3`.

The option `strided` will place `@strided` in front of the broadcasting operation.
You need `using Strided` for this to work.

    @reduce sum(i) A[i] * log(@reduce [i] := sum(j) A[j] * exp(B[i,j]))
    @cast W[i] := A[i] * exp(- @reduce S[i] = sum(j) exp(B[i,j]) lazy)

Recursion like this is allowed, inside either `@cast` or `@reduce`.
The array being created need not have a name, like `[i]`,
unless writing into an existing array, like `S[i]` here.
"""
macro reduce(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@reduce", exs...)) # wtf?
    _macro(exs...; reduce=true, where=where)
end

"""
    @reduce! Z[j] := sum(i,k) A[i,j,k]

Variant of `@reduce` which effectively runs `@check!()` on each tensor.
"""
macro reduce!(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@reduce!", exs...))
    _macro(exs...; reduce=true, where=where, icheck=true)
end

#=

This much is now per RHS term. None need exist outside the function:

nameA, indA, indAsub
sizeA      -- known & used for fixing ranges
getB, numB -- for view(A, get), and number of fixed indices
sizeC      -- after reshaping outer indices, in terms of sz_i
codeD      -- for gluing
indEflat   -- after gluing
negF       -- done before permutedims
shiftG     -- ditto, but not yet written
permH      --
codeI      -- to orient, if lacking some indices.

Of these things from LHS, only really store + canon need to be passed to RHS function inputex()
via walker(). But many need to be passed from readleft() to output, packaged into:
outUZ = (redUind, negV, codeW, sizeX, getY, numY, indZ, sizeZ, nameZ)

flags      -- todo, to-done, and options
store      -- information from parse! with sizes etc.
canon      -- canonical list of indices
canonsize  -- filled with size(A,2) etc, from store, at most one (:)

redUind    -- if reducing, this is done before slicing
negV       -- directions to reverse, done here.
codeW      -- for slicing. sizeWstatic is Size() of the slice.
sizeX      -- after slicing, container size, in sz_i
getY, numY -- for view(Z, get) in-place, and number fixed
sizeZ      -- for final reshape, in sz_i
nameZ, indZ, indZsub, indZred -- given LHS

=#

function _macro(exone, extwo=nothing, exthree=nothing;
    reduce=false, icheck=false, where=nothing, recurse=false)

    flags = Set{Symbol}()

    if reduce
        #===== parse @reduce input =====#

        if @capture(exone, left_ := redfun_(redind__) )     # Z[i] := sum(j) A[i] * B[j]
            right = extwo
            options = exthree
            push!(flags, :reduce)
            V && @info "partial reduce" left right redfun Tuple(redind) options
        elseif @capture(exone, left_ = redfun_(redind__) )
            right = extwo
            options = exthree
            push!(flags, :reduce)
            push!(flags, :inplace)
        elseif @capture(exone, left_ |= redfun_(redind__) )
            right = extwo
            options = exthree
            push!(flags, :reduce)
            push!(flags, :mustcopy)

        elseif @capture(exone, redfun_(redind__) )          # sum(i,j) A[i] * B[j]
            right = extwo
            options = exthree
            push!(flags, :reduce)
            push!(flags, :scalar)
            V && @info "full reduce" right redfun Tuple(redind)

        else
            throw(MacroError("don't know what to do with $exone", where))
        end

    else
        #===== parse @cast input =====#

        if @capture(exone, left_ := right_ )                # Z[i,j] := A[i] * B[j]
        elseif @capture(exone, left_ = right_ )
            push!(flags, :inplace)
        elseif @capture(exone, left_ |= right_ )
            push!(flags, :mustcopy)
        elseif @capture(exone, left_ == right_ )
            push!(flags, :mustview)
                                                            # A[i,j] -> Z[j][i]
        elseif @capture(exone, right_ -> left_ ) ||  @capture(exone, right_ => left_ )
            push!(flags, :anonfunc)
        else
            throw(MacroError("@cast doesn't know what to do with $exone", where))
        end
        V && @info "no reduction" left right
        options = extwo
        redind = []
        redfun = identity
        exthree == nothing || throw(MacroError("@cast doesn't know what to do with $exthree", where))

    end

    #===== parse LHS to get canonical list =====#

    store = SizeDict()
    canon, outUZ, nameZ, checkZ = readleft(left, redind, flags, store, icheck, where)

    # parse options both to look for keywords and sizes
    @capture(options, (optvec__,)) || (optvec = Any[options])
    optind, _,_,_ = parse!(store, nothing, [], optvec; allowranges=true, flags=flags)

    if count(i -> i != nothing, setdiff(optind, canon)) > 0
        str = join(something.(setdiff(optind, canon), "nothing"), ", ")
        m_error("attempting to ignore unrecognised options: $str", where)
    end

    #===== parse and process RHS =====#

    outex = quote end

    if @capture(right, AA_[ii__] ) || @capture(right, AA_{ii__} )
        newright = walker(outex, right, canon, flags, store, icheck, where)
    else
        newright = MacroTools.prewalk(
            x -> walker(outex, x, canon, flags, store, icheck, where), right)
        push!(flags, :broadcast)
    end

    notseen = setdiff(canon, unique(store.seen))
    isempty(notseen) || throw(MacroError("did not see index $(join(notseen, ", ")) on the right", where))

    #===== almost done =====#

    packagecheck(flags, where) # disabled for issue #2

    canonsize = sizeinfer(store, canon, where, true)

    V && @info "before in/out choice" store flags Tuple(canonsize)

    if :inplace in flags

        #===== in-place output  =====#

        if checkZ != nothing
            push!(outex.args, checkZ)
        end

        inout = outputinplace(newright, outUZ, redfun, canonsize, canon, flags, store, nameZ, where)
        push!(outex.args, inout)

        if :finalres in flags
            push!(outex.args, nameZ)
        end

    else
        #===== out-of-place output  =====#

        if :broadcast in flags
            if (:reduce in flags) && (:lazy in flags)
                newright = :( Base.@__dot__ TensorCast.lazy($newright) ) # really LazyArrays.lazy
            elseif :strided in flags
                dotted = Broadcast.__dot__(newright) # @strided does not work on @.
                newright = :( Strided.@strided $dotted )
            else
                newright = :( Base.@__dot__ $newright ) # prettier in @pretty
            end
        end

        finalright = outputnew(newright, outUZ, redfun, canonsize, canon, flags, store, where)
        push!(outex.args, :( $nameZ =  $finalright ) )

        if checkZ != nothing
            push!(outex.args, checkZ)
            push!(outex.args, nameZ)
        end
    end

    #===== finalise =====#

    if :needsize in flags
        szcanon = Any[ Symbol(:sz_,i) for i in canon ]
        pushfirst!(outex.args, :( local ($(szcanon...),) = ($(canonsize...),) ) )
    end

    if :assert in flags || :(!) in flags || check_options.size
        for ch in store.checks
            pushfirst!(outex.args, ch)
        end
    end

    for tex in store.topex
        pushfirst!(outex.args, tex)
    end

    if recurse==true # only for @reduce inside something
        indZ = outUZ[end-1]
        return outex, indZ
    end

    if :anonfunc in flags
        leftex = anoninput(store.rightnames, where)
        outex = :( $leftex -> $outex )
        outex.args[2] = MacroTools.unblock(outex.args[2])
        return esc(outex)
    else
        outex = MacroTools.unblock(outex) # removes begin end if it can
        return esc(outex)
    end
end

#="""
    canon, outUZ, nameZ, checkZ = readleft(left, redind, flags, store, icheck, where)

outUZ = (redUind, negV, codeW, sizeX, getY, numY, indZ, sizeZ)
are things passed to output construction.
"""=#
function readleft(left, redind, flags, store, icheck, where)

    if @capture(left, Z_[outer__][inner__]) ||  @capture(left, [outer__][inner__])
        push!(flags, :slice)
    elseif @capture(left, Z_[outer__]{inner__}) || @capture(left, [outer__]{inner__})
        push!(flags, :staticslice)
    elseif @capture(left, Z_[outer__]) || @capture(left, [outer__])
        inner = []
    elseif left==nothing
        @assert :scalar in flags
        inner = []
        outer = []
        Z = nothing
    else
        throw(MacroError("readleft doesn't know what to do with $left", where))
    end

    if Z == nothing
        nameZ = gensym(:Z) # no output name
        if :inplace in flags
            throw(MacroError("can't write in-place into nowhere!", where))
        end
    else
        nameZ = Z
    end

    if !(:inplace in flags)
        Z = nothing # tells parse! not to think about size(Z,...)
    end
    flat12, getY, sizeZ, negV = parse!(store, Z, outer, inner; allowranges=true) # true allows A[i]{j:3}
    redUind, _,_,_ = parse!(store, nothing, [], redind; allowranges=true) # true means sum(i:3) allowed

    if length(flat12)-length(inner) == 2 && flat12[end] == flat12[end-1] # Z[i,i] special case
        length(getY) == 2 || throw(MacroError("can't fix output index with Diagonal output", where)) # because I'm lazy, this would not be impossible
        push!(flags, :diagleft)
        pop!(flat12) # remove index from flat list
        pop!(getY)
        pop!(sizeZ)
    end

    canon = vcat(flat12, redUind) # order here is [inner, outer, reduced]

    checkrepeats(flat12, " on left hand side", where)
    checkrepeats(canon, " in reduction function", where)

    codeW = repeat(Any[*], length(canon) - length(redUind))
    codeW[1:length(inner)] .= (:)

    indW = canon[1:length(inner)] # inner without minuses etc

    indX = flat12[length(inner)+1:end] # flattened outer without minuses, also without fixed
    sizeX = Any[ Symbol(:sz_,i) for i in indX ] # only for in-place

    numY = count(!isequal(:), getY) # the number of fixed indices in output

    if length(getY) - numY != length(sizeX)
        push!(flags, :backshape) # whether reshaping of view(Z, getY) is needed, only in-place
    end

    if (length(sizeZ) + numY) != (length(canon) - length(inner) - length(redUind))
        push!(flags, :outshape) # whether final reshaping is needed, for := case
    end

    indZ = [] # like indX, but with fixed indices re-inserted, for :named / :axis
    iz = 1
    for g in getY
        if g == (:)
            push!(indZ, indX[iz])
            iz += 1
        else
            push!(indZ, g)
        end
    end

    checkZ = nothing
    if icheck
        checknow = check_one(:($nameZ[$(outer...)]), where)
        if check_options.size
            checkZ = checknow # check!(...) to be inserted
        end
    end

    outUZ = (redUind, negV, codeW, indW, sizeX, getY, numY, indZ, sizeZ)
    V && @info "readleft" Tuple(canon) outUZ nameZ
    return canon, outUZ, nameZ, checkZ
end

#="""
    walker(outex, x, canon, flags, store, icheck, where)

Called by `MacroTools.prewalk` on RHS, finds tensors & pushes `:( sym = inputex(this) )` into
`outex.args`, and then replaces them with `sym`.
"""=#
function walker(outex, ex, canon, flags, store, icheck, where)

    # allow @reduce inside RHS, by simply calling the entire macro first
    if @capture(ex, @reduce(redex__))
        Rsym = gensym(:R)
        Rval, Rind = _macro(redex...; reduce=true, where=where, recurse=true)

        # sew the result of that into ... not outex but topex, as may need its size:
        push!(store.topex, :(local $Rsym = $Rval ) )

        # construct name for resulting tensor, for outer @cast to further process
        ex =  :( $Rsym[$(Rind...)] )

    elseif @capture(ex, @mul(mulex__))
        Msym = gensym(:M)
        Mval, Mind = _mul(mulex...; where=where, recurse=true)

        push!(store.topex, :(local $Msym = $Mval ) )

        ex =  :( $Msym[$(Mind...)] )
    end

    # find any indexed tensor, process it, and replace it
    if @capture(ex, A_[ij__][kl__] ) || @capture(ex, A_[ij__]{kl__} ) || @capture(ex, A_[ij__] )

        push!(store.rightnames, A) # used for @cast A[i,j] -> B[i\j]

        if isa(A, Symbol) || @capture(A, Astruct_.field_)

        # allow something like mapslices:
        # @pretty @cast Z[i,j] := fun(X[:,i,:],42)[j] + X[j,i,3] + Y[j]^2
        elseif @capture(A, f_(B_[i_colons__],rest__) )

            kl == nothing || throw(MacroError("can't nest this deeply: $ex", where))

            # i_colons = vcat(rvec(ii1), rvec(i1), Any[:(:)], rvec(i2), rvec(ii2)) # NB contains the symbol :(:) not the function (:)

            A, ex = slicemap(f, B, i_colons, rest, ij, store)

        # if we have f(x)[i,j] then we should evaluate f just once, e.g.
        # @pretty @cast A[i]{j:5} |= rand(15)[i\j]
        else
            Atop = gensym(:A)
            push!(store.topex, :(local $Atop = $A ) )
            A = Atop
        end

        V && @info "walker before inputex" A ex canon store

        Aval = inputex(A, ex, canon, flags, store, icheck, where)

        if isa(Aval, Symbol)
            ex = Aval
        else
            Asym = gensym(:A)
            push!(outex.args, :(local $Asym = $Aval) )
            ex =  Asym
        end
        V && @info "walker" ex Aval

    # catch A[i,j]' as element-wise conj
    elseif @capture(ex, arg_')
        ex = :( Base.conj($arg) )

    # try to protect log(2) etc from @. , but not log(A[i])
    elseif @capture(ex, f_(x_Symbol) ) || @capture(ex, f_(x_Int) ) || @capture(ex, f_(x_Float64) )
        fxsym = gensym(:fx)
        push!(outex.args, :( local $fxsym = $ex ))
        ex = fxsym
    end

    return ex
end

#="""
    slicemap(f, :B, [:i, :(:), :j], rest, kl, store)

This processes `f(B[i,:mj],rest...)[k,l]`, such as in these:
```
@pretty @cast Z[i,j] := fun(X[:,i,:],42)[j] + X[j,i,3] + Y[j]^2
@pretty @cast Z[i,j\k] := fun(X[:,i\k,:])[j] + Y[k]^2
```

* Note that `f(X[i,:], Y[i,:])` is not allowed (yet).

* `g(X[i,:,_])` should work, but notices _ after map,
and thus `g(X[i,:,3])` will evaluate g on many slices which aren't required.

These could be fixed by upgrading this function -- TODO maybe?
"""=#
function slicemap(f, B, i_colons, rest, kl, store)

    if isa(B, Symbol) || @capture(B, Bstruct_.field_)
    else # pull out B = B(x) so that it's evaluated only once
        Btop = gensym(:B)
        push!(store.topex, :(local $Btop = $B ) )
        B = Btop
    end

    # assertions, and size information for resising etc.
    Bdims = length(i_colons)
    str = "expected a $Bdims-tensor $B[" * join(i_colons, ", ") * "]"
    push!(store.checks, :(TensorCast.@assert_ ndims($B)==$Bdims $str) )
    for (dout, ind) in enumerate(i_colons)
        if ind != :(:) && !isconstant(ind)
            savesize!(store, ind, :( size($B, $dout) ) )
        end
    end

    # construct map(f,eachslice(B)) expression
    Atop = gensym(:slicemap)
    i_code = map(i -> i==:(:) ? (:) : (*), Tuple(i_colons))
    if count(isequal(:), i_code) > 0
        sliceex = :( TensorCast.sliceview($B, $i_code) )
    else
        sliceex = B
    end
    if length(rest) == 0
        push!(store.topex, :(local $Atop = map($f, $sliceex) ) )
    else
        Xsymb = gensym(:X)
        push!(store.topex, :(local $Atop = map($Xsymb -> $f($Xsymb, $(rest...)), $sliceex) ) )
    end

    # construct ex for inputex later
    ij = filter(!isequal(:(:)), i_colons)
    ex = :( $Atop[$(ij...)][$(kl...)] )

    return Atop, ex
end

#="""
    inputex(:A, :( A[i,j][k] ), target, flags, store, icheck, where)

Figures out all the steps needed to transform the given tensor to a boring one,
aligned with the given target, and returns the necessary expression.
Writes sizes which can be read from `A` into `store`, and necessary reshapes in terms of `sz_i`.

Now `target` need not be `== canon`, since `sz_i` is independent of that.

Now needs explicitly the name of tensor `A`,
so that when `inex = rand(2,3)[i,j]` this is not evaluated twice.
"""=#
function inputex(A, inex, target, flags, store, icheck, where)

    if @capture(inex, Aa_[outer__][inner__])
        glue = :yes
    elseif @capture(inex, Aa_[outer__]{inner__})
        glue = :static
    elseif @capture(inex, Aa_[outer__])
        inner = []
        glue = :no
    else error("inputex should not have been called")
    end

    flatE, getB, _, negF = parse!(store, A, outer, inner)

    append!(store.seen, flatE)

    ex = A
    if icheck
        ex = check_one(:($A[$(outer...)]), where)           # @check!
    end

    numB = count(!isequal(:), getB)
    if numB > 0                                             # A[_,i]
        if numB == length(getB) || :nolazy in flags
            needview!(getB) # replace _ with 1 if any
            ex = :( $ex[$(getB...)] )
        elseif needview!(getB) # then at least one index fixed and not _
            ex = :(view($ex, $(getB...) ))
        else
            ex = :(TensorCast.rview($ex, $(getB...) )) # really a reshape
        end
    end

    sizeC = Any[ Symbol(:sz_, i) for i in flatE[length(inner)+1 : end] ]

    if (length(sizeC) + numB) != length(outer)              # A[i\j]
        sizeCex = :(($(sizeC...) ,))
        ex = :( reshape($ex, $sizeCex) )
        push!(flags, :needsize)
    end

    if length(flatE)-length(inner) == 2 && flatE[end] == flatE[end-1]  # A[i,i] special case
        if :nolazy in flags
            ex = :( TensorCast.diag($ex) ) # LinearAlgebra might not be loaded by caller
            push!(flags, :havecopied)
        else
            ex = :( TensorCast.diagview($ex) )
        end
        pop!(flatE) # remove repeated index from the end
    end

    checkrepeats(flatE, " in term $inex", where) # after dealing with diag story

    dirs = [ findcheck(i, target, where) for i in flatE ]

    if glue == :yes                                         # A[i][k]
        codeD = repeat(Any[*],length(flatE))
        codeD[1:length(inner)] .= (:)

        if codeD == [:,*] && dirs[1] > dirs[2] # then we can avoid a transpose
            codeD = [*,:]
            dirs = [dirs[2], dirs[1]]
        end
        # you could perform more elaborate versions of that, e.g. for this:
        # @pretty @reduce A[i\j,_] = sum(k) B[i,j][k]
        # however only copy_glue and julienne_glue understand arbitrary codes

        if :lazy in flags || :mustview in flags
            ex = :( TensorCast.lazy_glue($ex, $(codeD...,))  )

        elseif :cat in flags
            ex = :( TensorCast.cat_glue($ex, $(codeD...,))  )
            push!(flags, :havecopied)
        elseif :glue in flags
            ex = :( TensorCast.copy_glue($ex, $(codeD...,))  )
            push!(flags, :havecopied)
        elseif :julienne in flags
            ex = :( TensorCast.julienne_glue($ex, $(codeD...,))  )
            push!(flags, :havecopied)
        else
            ex = :( TensorCast.red_glue($ex, $(codeD...,))  )
            push!(flags, :havecopied)
        end

    elseif glue == :static                                  # A[i]{k}
        ex = :( TensorCast.static_glue($ex)  )
        push!(flags, :staticglue) # for packagecheck
    end

    perm = ntuple(identity, length(dirs))
    if dirs != sort(dirs)                                   # A[j,i]
        perm = Tuple(sortperm(dirs))
        if :nolazy in flags
            ex = :( permutedims($ex, $perm) )
            push!(flags, :havecopied)
        elseif :strided in flags
            ex = :( Strided.@strided permutedims($ex, $perm) )
        elseif perm == (2,1)
            ex = :( TensorCast.PermuteDims($ex) ) # avoids transpose, except on numbers, because it's recursive
        else
            ex = :( PermutedDimsArray($ex, $perm) )
        end
        dirs = sort(dirs) # avoids unnecessary orient(PermuteDims(M), (:,:,*))
    end

    for i in negF                                           # A[-i,j]
        d = invperm(perm)[findcheck(i, flatE, where)]
        if :nolazy in flags
            ex = :( reverse($ex, dims=$d) )
            push!(flags, :havecopied)
        else
            ex = :( TensorCast.Reverse{$d}($ex) )
        end
    end

    if length(flatE) != length(target) && dirs != 1:length(dirs) # A[i] + B[j]
        codeH = repeat(Any[*],length(target))
        codeH[dirs] .= (:)
        while codeH[end] == *
            pop!(codeH)
        end
        ex = :( TensorCast.orient($ex, $(codeH...,)) )
    end

    return ex
end

#="""
     outputnew(newright, outUZ, redfun, canonsize, canon, flags, store, where)

For the case of `:=`, this constructs the expression to do reduction if needed,
and slicing/reshaping/reversing for LHS.

outUZ = (redUind, negV, codeW, indW, sizeX, getY, numY, indZ, sizeZ)
"""=#
function outputnew(newright, (redUind, negV, codeW, indW, sizeX, getY, numY, indZ, sizeZ),
        redfun, canonsize, canon, flags, store, where)

    ex = newright

    for ri in negV                                          # Z[-i,j]
        d = findcheck(ri, canon, where)
        ex = :( reverse($ex, dims=$d) )
    end

    if :reduce in flags && :scalar in flags                 # Z = @reduce sum(i) ...
        ex = :( $redfun($ex) )
    elseif :reduce in flags                                 # Z[i] := sum(j) ...
        rdims = Tuple([findcheck(i, canon, where) for i in redUind])
        if length(rdims)==1
            rdims = first(rdims)
        end
        ex = :( $redfun($ex, dims=$rdims) )
        if :outshape in flags && !(:slice in flags || :staticslice in flags)
            # then we need not dropdims, as reshape will handle it
        else
            ex = :( dropdims($ex, dims=$rdims) )
        end
    end

    if :slice in flags                                      # Z[i][k] :=
        if :mustcopy in flags # && !(:havecopied in flags) # make |= slightly stronger
            ex = :( TensorCast.slicecopy($ex, $(codeW...,)) )
            push!(flags, :havecopied)
        elseif :julienne in flags
            ex = :( TensorCast.julienne_slice($ex, $(codeW...,)) )
        else
            ex = :( TensorCast.sliceview($ex, $(codeW...,)) )
        end
    elseif :staticslice in flags                            # Z[i]{k} :=
        # codeW worked out already, but sizeWstatic must be done here
        # TODO sizeinfer!(store, canon, where, false)
        sizeWstatic = :( StaticArrays.Size($([store.dict[i] for i in indW]...)) )
        if :outshape in flags
            ex = :( TensorCast.static_slice($ex, $sizeWstatic, false) )
        else
            ex = :( TensorCast.static_slice($ex, $sizeWstatic) )
        end
    end

    if :outshape in flags                                   # Z[i\j, _, k]
        sizeZex = :(($(sizeZ...) ,))
        ex = :( reshape($ex, $sizeZex) )
        push!(flags, :needsize)
        for n in filter(!isequal(:), getY)
            n == 1 || n == :_ || throw(MacroError("can't fix output index to $n != 1, when creating a new array", where))
        end
    end

    if :mustcopy in flags                                   # Z[i] |= ...
        if !(:havecopied in flags) && !(:broadcast in flags)
            ex = :( collect($ex) )
        elseif (:slice in flags) && (:strided in flags)     # Z[i][j] |= ... strided
            ex = :( collect.($ex) ) # because slicecopy actually makes views here
        elseif (:staticslice in flags) || (:strided in flags) # i.e. then copy even if there was a broadcast etc.
            ex = :( collect($ex) )
        end

    elseif :mustview in flags                               # Z[i] == ...
        if :havecopied in flags
            m_error("can't do what you ask without copying, sorry", where)
        elseif :broadcast in flags
            m_error("can't broadcast without copying, sorry", where)
        end
    end

    if :diagleft in flags                                   # Z[i,i] := ...
        if :nolazy in flags
            ex = :( TensorCast.diagm(0 => $ex) )
        else
            ex = :( TensorCast.Diagonal($ex) )
        end
    end

    if :namedarray in flags                                 # Z[i] := ... named
        ex = :( TensorCast.namedarray($ex, $(indZ...,) ) )
    elseif :nameddims in flags || :named in flags
        ex = :( NamedDims.NamedDimsArray{$(indZ...,)}($ex) )
    # elseif :axis in flags
    #     ex = :( TensorCast.axisarray($ex, $(indZ...,) ) )
    end

    return ex
end

#="""
    outputinplace(newright, outUZ, redfun, canonsize, canon, flags, store, nameZ, where)

For the case of `=` this figures out how to write RHS into LHS, in one of three ways:
* reduction `sum!(Z, newright)`
* broadcasting `@. Z[...] = newright`
* neither, `copyto!(Z, newright)`

No longer attempts to write `permutedims!(Z, A, ...)`, now just `copyto!(Z, PermutedDimsArray(A, ...))`.
Doesn't really need so many arguments...
"""=#
function outputinplace(newright, (redUind, negV, codeW, indW, sizeX, getY, numY, indZ, sizeZ),
        redfun, canonsize, canon, flags, store, nameZ, where)

    if :slice in flags
        throw(MacroError("can't write to sliced arrays in-place, for now", where))
    end
    if length(negV) > 0
        m_error("can't reverse axes of in-place output, try moving -$(negV[1]) to right hand side", where)
    end

    if :reduce in flags                                     # sum!(revleft, newright)

        if :broadcast in flags
            if :lazy in flags
                newright = :( Base.@__dot__ TensorCast.lazy($newright) ) # really LazyArrays.lazy
            elseif :strided in flags
                dotted = Broadcast.__dot__(newright) # @strided does not work on @.
                newright = :( Strided.@strided $dotted )
            else
                newright = :( Base.@__dot__ $newright ) # prettier in @pretty
            end
        end

        # working backwards
        revleft = nameZ

        if numY > 0
            if needview!(getY)
                revleft = :(view($revleft, $(getY...) ))
            else
                revleft = :(TensorCast.rview($revleft, $(getY...) )) # really a reshape
            end
            push!(flags, :finalres)
        end

        if :diagleft in flags
            revleft = :( TensorCast.diagview($revleft) )
            push!(flags, :finalres)
        end

        if :backshape in flags
            sizeXex = :(($(sizeX...) ,))
            revleft = :( reshape($revleft, $sizeXex) )
            push!(flags, :needsize)
            push!(flags, :finalres)
        end

        if !endswith(string(redfun), '!')
            redfun = Symbol(redfun, '!')
        end

        ex = :( $redfun($revleft, $newright) )

    elseif :broadcast in flags                              # @. revleft[...] = newright

        # working backwards
        revleft = nameZ

        if numY > 0 # when getY has only : and 1, and backshape, then you could skip this
            if needview!(getY)
                revleft = :(view($revleft, $(getY...) ))
            else
                revleft = :(TensorCast.rview($revleft, $(getY...) )) # really a reshape
            end
            push!(flags, :finalres)
        end

        if :diagleft in flags
            revleft = :( TensorCast.diagview($revleft) )
            push!(flags, :finalres)
        end

        if :backshape in flags
            sizeXex = :(($(sizeX...) ,))
            revleft = :( reshape($revleft, $sizeXex) )
            push!(flags, :needsize)
            push!(flags, :finalres)
        end

        if :strided in flags
            dotted = Broadcast.__dot__(newright)
            ex = :( Strided.@strided $revleft .= $dotted )
        else
            ex = :( $revleft .= Base.@__dot__ $newright ) # don't apply @. to reshape(diagview(Z))
        end

    else                                                    # copyto!(revleft, newright)

        if :strided in flags
            newright = :( Strided.@strided $newright )
        end

        # working backwards
        revleft = nameZ

        if numY > 0
            if needview!(getY)
                revleft = :(view($revleft, $(getY...) ))
            else
                revleft = :(TensorCast.rview($revleft, $(getY...) )) # really a reshape
            end
            push!(flags, :finalres)
        end

        if :diagleft in flags
            revleft = :( TensorCast.diagview($revleft) )
            push!(flags, :finalres)
        end

        ex = :( $copyto!($revleft, $newright) )

    end

    return ex
end

#="""
    anoninput(store.rightnames)

Given a vector of names captured from `A_[i...]`, returns ex needed for `(A,B) -> ...`
"""=#
function anoninput(rightnames, where)
    for A in rightnames
        isa(A,Symbol) || throw(MacroError("can't use $A as anonymous function input", where))
        # TODO make this understand that $B means interpolate not arg
    end
    if length(rightnames) == 1
        return rightnames[1]
    else
        return :( ($(rightnames...),) )
    end
end

using LinearAlgebra  # for diag()

function packagecheck(flags, where)
    where === nothing && return
    # now check in caller's scope?
    if :staticslice in flags || :staticglue in flags
        isdefined(where.mod, :StaticArrays) || m_error("can't use static arrays without using StaticArrays", where)
    end
    if :strided in flags
        isdefined(where.mod, :Strided) || m_error("can't use option strided without using Strided", where)
    end
    if :julienne in flags
        isdefined(where.mod, :JuliennedArrays) || m_error("can't use option julienne without using JuliennedArrays", where)
    end
    if :namedarray in flags
        isdefined(where.mod, :NamedArrays) || m_error("can't use option namedarray without using NamedArrays", where)
    end
    if :named in flags || :nameddims in flags
        isdefined(where.mod, :NamedDims) || m_error("can't use option nameddims without using NamedDims", where)
    end
    # if :axis in flags
    #     isdefined(where.mod, :AxisArrays) || m_error("can't use option axis without using AxisArrays", where)
    # end
end


text:
# Exact Search
# ============
#
# Exact sequence search tools.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Reference Implementation
# ------------------------
#
# `searchindex` and `rsearchindex` should work exactly the same way as the
# following reference implementation:
#=
# Return the index of the first occurrence of `pat` in `seq[start:stop]`.
function searchindex(seq, pat, start=1, stop=lastindex(seq))
    m = length(pat)
    n = length(seq)
    for s in max(start-1, 0):min(stop, n)-m
        if occurs_with_shift(pat, seq, s)
            return s+1  # found
        end
    end
    return 0  # not found
end

# Return the index of the last occurrence of `pat` in `seq[stop:start]`.
function rsearchindex(seq, pat, start=lastindex(seq), stop=1)
    n = length(seq)
    m = length(pat)
    for s in min(start, n)-m:-1:max(stop-1, 0)
        if occurs_with_shift(pat, seq, s)
            return s+1  # found
        end
    end
    return 0  # not found
end

function occurs_with_shift(pat, seq, shift)
    for i in 1:lastindex(pat)
        if !iscompatible(pat[i], seq[shift+i])
            return false
        end
    end
    return true
end
=#

# Actual Implementation
# ---------------------

"""
Query type for exact sequence search.
"""
struct ExactSearchQuery{S<:BioSequence}
    seq::S         # query sequence
    cbits::UInt32  # compatibility bits
    fshift::Int    # shift length for forward search
    bshift::Int    # shift length for backward search
end

function ExactSearchQuery(query::BioSequence)
    cbits, fshift, bshift = preprocess(query)
    return ExactSearchQuery(query, cbits, fshift, bshift)
end

function preprocess(query)
    if length(query) == 0
        return UInt32(0), 0, 0
    end
    m = length(query)
    first = query[1]
    last = query[end]
    cbits::UInt32 = 0
    fshift = bshift = m
    for i in 1:lastindex(query)
        x = query[i]
        cbits |= compatbits(x)
        if iscompatible(x, last) && i < m
            fshift = m - i
        end
    end
    for i in lastindex(query):-1:1
        x = query[i]
        if iscompatible(x, first) && i > 1
            bshift = i - 1
        end
    end
    return cbits, fshift, bshift
end

function checkeltype(seq1, seq2)
    if eltype(seq1) != eltype(seq2)
        throw(ArgumentError("the element type of two sequences must match"))
    end
end

# Forward
# -------

"""
    findfirst(pat, seq::BioSequence, [, start=1[, stop=lastindex(seq)]])

Return the index of the first occurrence of `pat` in `seq[start:stop]`; symbol
comparison is done using `BioSequences.iscompatible`.
"""
function Base.findfirst(pat::BioSequence, seq::BioSequence,
                        start::Integer=1, stop::Integer=lastindex(seq))
    return findfirst(ExactSearchQuery(pat), seq, start, stop)
end

function Base.findfirst(query::ExactSearchQuery, seq::BioSequence,
                        start::Integer=1, stop::Integer=lastindex(seq))
    checkeltype(seq, query.seq)
    i = quicksearch(query, seq, start, stop)
    if i == 0
        return nothing
    else
        return i:i+length(query.seq)-1
    end
end

# This algorithm is borrowed from base/strings/search.jl, which looks similar to
# Sunday's Quick Search algorithm.
function quicksearch(query, seq, start, stop)
    pat = query.seq
    m = length(pat)
    n = length(seq)
    stop′ = min(stop, n) - m
    s::Int = max(start - 1, 0)

    if m == 0  # empty query
        if s ≤ stop′
            return s + 1  # found
        else
            return 0  # not found
        end
    end

    while s ≤ stop′
        if iscompatible(pat[m], seq[s+m])
            i = m - 1
            while i > 0
                if !iscompatible(pat[i], seq[s+i])
                    break
                end
                i -= 1
            end
            if i == 0
                return s + 1  # found
            elseif s < stop′ && (query.cbits & compatbits(seq[s+m+1]) == 0)
                s += m + 1
            elseif isambiguous(seq[s+m])
                s += 1
            else
                s += query.fshift
            end
        elseif s < stop′ && (query.cbits & compatbits(seq[s+m+1]) == 0)
            s += m + 1
        else
            s += 1
        end
    end

    return 0  # not found
end

# Backward
# --------

"""
    findlast(pat, seq::BioSequence, [, start=1[, stop=lastindex(seq)]])

Return the index of the last occurrence of `pat` in `seq[start:stop]`; symbol
comparison is done using `BioSequences.iscompatible`.
"""
function Base.findlast(pat::BioSequence, seq::BioSequence,
                       start::Integer=lastindex(seq), stop::Integer=1)
    return findlast(ExactSearchQuery(pat), seq, start, stop)
end

function Base.findlast(query::ExactSearchQuery, seq::BioSequence,
                       start::Integer=lastindex(seq), stop::Integer=1)
    checkeltype(seq, query.seq)
    i = quickrsearch(seq, query, start, stop)
    if i == 0
        return nothing
    else
        return i:i+length(query.seq)-1
    end
end

function quickrsearch(seq, query, start, stop)
    pat = query.seq
    m = length(pat)
    n = length(seq)
    stop′ = max(stop - 1, 0)
    s::Int = min(start, n) - m

    if m == 0  # empty query
        if s ≥ stop′
            return s + 1  # found
        else
            return 0  # not found
        end
    end

    while s ≥ stop′
        if iscompatible(pat[1], seq[s+1])
            i = 2
            while i < m + 1
                if !iscompatible(pat[i], seq[s+i])
                    break
                end
                i += 1
            end
            if i == m + 1
                return s + 1  # found
            elseif s > stop′ && (query.cbits & compatbits(seq[s]) == 0)
                s -= m + 1
            elseif isambiguous(seq[s+1])
                s -= 1
            else
                s -= query.bshift
            end
        elseif s > stop′ && (query.cbits & compatbits(seq[s]) == 0)
            s -= m + 1
        else
            s -= 1
        end
    end

    return 0  # not found
end

"""
occursin(x:BioSequence, y::BioSequence)
occursin(x:ExactSearchQuery, y::BioSequence)

Return Bool indicating presence of exact match of x in y.
"""
Base.occursin(x::BioSequence, y::BioSequence) = occursin(ExactSearchQuery(x), y)
Base.occursin(x::ExactSearchQuery, y::BioSequence) = quicksearch(x, y, 1, length(y)) != 0


text:
@testset "Statistics" begin

	set = BitsMNIST.Datasets.mnist()
	sx, sy = BitsMNIST.ZeroOne.sample(set["train_x"], set["train_y"], 0.01)
	model = Chain(BitDense(784, 2, true, sigmoid))
	error = BitsMNIST.Statistics.error(model, sx, sy)

	@test error >= 0
	@test error <= 1

end


text:
function test_energy(mesh, sbp, eqn, opts)
  energy_final = calcNorm(eqn, eqn.q_vec)
  q_initial = zeros(eqn.q_vec)
  ICfunc = AdvectionEquationMod.ICDict[opts["IC_name"]]
  ICfunc(mesh, sbp, eqn, opts, q_initial)
  energy_initial = calcNorm(eqn, q_initial)

  @test  abs(energy_initial - energy_final)  < 1e-12
end

function test_energy_parallel()
  @testset "----- Testing Energy Stability -----" begin

    start_dir = pwd()

    fname = "input_vals_periodic.jl"
    cd("./2dp4")
    mesh, sbp, eqn, opts = solvePDE(fname)
    test_energy(mesh, sbp, eqn, opts)

    cd("../3dp4")
    mesh, sbp, eqn, opts = solvePDE(fname)
    test_energy(mesh, sbp, eqn, opts)

  end

  return nothing
end

test_energy_parallel()


text:
#__precompile__()

"Main module for `LTspice.jl` - a Julia interface to LTspice"
module LTspice
using StringEncodings
using Dates: DateTime
import IterTools.chain

include("open_with_unknown_encoding.jl")
include("specialarrays.jl")
include("LTspiceSimulation.jl")
include("ParseCircuitFile.jl")
include("ParseLogFile.jl")
include("perlineiterator.jl")
include("utility.jl")

end # module


text:
using Images

img=load("julia.jpeg")
@show img

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


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

__precompile__(true)

module SolidStateDetectors

using LinearAlgebra
using Random
using Statistics

using ArraysOfArrays
using FillArrays
using Formatting
using Interpolations
using IntervalSets
using JSON
using LaTeXStrings
using ParallelProcessingTools
using ProgressMeter
using RadiationDetectorSignals
using RecipesBase
using Requires
using Rotations
using StaticArrays
using StatsBase
using Unitful
using YAML

import Clustering
import Distributions
import Tables
import TypedTables

import Base: size, sizeof, length, getindex, setindex!, axes, getproperty, broadcast,
             range, ndims, eachindex, enumerate, iterate, IndexStyle, eltype, in
import Base: show, print, println, display, +, -, &
import Base.convert

export SolidStateDetector
export SSD_examples

export Grid, CylindricalPoint, CartesianPoint

export ElectricPotential, PointTypes, ChargeDensity, DielectricDistribution, WeightingPotential, ElectricField
export apply_initial_state!
export calculate_electric_potential!, calculate_weighting_potential!, calculate_electric_field!
export update_till_convergence!, refine!
export set_charge_drift_model!, calculate_drift_fields!
export get_active_volume
export generate_charge_signals, generate_charge_signals!
export ElectricFieldChargeDriftModel, ADLChargeDriftModel
export Simulation, simulate!
export Event, drift_charges!

const SSDFloat = Union{Float16, Float32, Float64}

struct ConfigFileError <: Exception
    msg::AbstractString
end
Base.showerror(io::IO, e::ConfigFileError) = print(io, "ConfigFileError: ", e.msg)

include("Geometries/Geometries.jl")

include("Axes/DiscreteAxis.jl")
include("World/World.jl")
include("Grids/Grids.jl")

include("Types/Types.jl")

include("MaterialProperties/MaterialProperties.jl")
include("Config/Config.jl")
include("ChargeDensityModels/ChargeDensityModels.jl")
include("SolidStateDetector/DetectorGeometries.jl")
include("GeometryRounding.jl")

include("PotentialSimulation/PotentialSimulation.jl")

include("ElectricField/ElectricField.jl")

include("ChargeDriftModels/ChargeDriftModels.jl")
include("ChargeCloudModels/ChargeCloudModels.jl")
include("ChargeDrift/ChargeDrift.jl")
include("SignalGeneration/SignalGeneration.jl")

include("ChargeStatistics/ChargeStatistics.jl")
include("ChargeClustering/ChargeClustering.jl")

include("Simulation/Simulation.jl")
include("Event/Event.jl")
include("MCEventsProcessing/MCEventsProcessing.jl")

include("IO/IO.jl")

include("examples.jl")

include("PlotRecipes/PlotRecipes.jl")

function __init__()
    @require HDF5="f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" begin
        @require LegendHDF5IO="c9265ca6-b027-5446-b1a4-febfa8dd10b0" begin
            include("IO/hdf5_specific.jl")
        end
        include("MCEventsProcessing/MCEventsProcessing_hdf5.jl")
    end
end

end # module


text:
module SpringCollab2020CrystalBombDetonator

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/crystalBombDetonator" CrystalBombDetonator(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight)

const placements = Ahorn.PlacementDict(
    "Crystal Bomb Detonation Field (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CrystalBombDetonator,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::CrystalBombDetonator) = 8, 8
Ahorn.resizable(entity::CrystalBombDetonator) = true, true

function Ahorn.selection(entity::CrystalBombDetonator)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CrystalBombDetonator, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.45, 0.0, 0.45, 0.8), (0.0, 0.0, 0.0, 0.0))
end

end


text:
# info("Running Test: testmonitor.jl")
# include("testmonitor.jl")

info("Running Test: testTAP.jl")
include("testTAP.jl")


text:
export EthHdr, IpFlags, IpHdr,
       UdpHdr, TcpFlags, TcpHdr,
       IcmpHdr, DecPkt, decode_pkt

mutable struct EthHdr
    dest_mac::AbstractString
    src_mac::AbstractString
    ptype::UInt16
    EthHdr() = new("","",0)
end # struct EthHdr

mutable struct VLAN
	priority::UInt8
	dei::UInt8
	id::Int16
	ptype::UInt16
	VLAN() = new(0,0,0,0)
end # strut 802.1Q Virtual LAN

mutable struct IpFlags
    reserved::Bool
    dont_frag::Bool
    more_frags::Bool
    IpFlags() = new(false,false,false)
end # struct IpFlags

mutable struct IpHdr
    version::UInt8
    length::UInt8
    services::UInt8
    totlen::UInt16
    id::UInt16
    flags::IpFlags
    frag_offset::UInt16
    ttl::UInt8
    protocol::UInt8
    checksum::Bool
    src_ip::AbstractString
    dest_ip::AbstractString
    IpHdr() = new(0,0,0,0,0,IpFlags(),0,0,0,false,"","")
end # struct IpHdr

mutable struct TcpFlags
    reserved::Bool
    nonce::Bool
    cwr::Bool
    ecn::Bool
    urgent::Bool
    ack::Bool
    push::Bool
    reset::Bool
    syn::Bool
    fin::Bool
    TcpFlags() = new(false,false,false,false,false,
                     false,false,false,false,false)
end # struct TcpFlags

mutable struct TcpHdr
    src_port::UInt16
    dest_port::UInt16
    seq::UInt32
    ack::UInt32
    offset::UInt8
    flags::TcpFlags
    window::UInt16
    checksum::UInt16
    uptr::UInt16
    data::Vector{UInt8}
    TcpHdr() = new(0,0,0,0,0,TcpFlags(),0,0,0, Vector{UInt8}(undef, 0))
end # struct TcpHdr

mutable struct UdpHdr
    src_port::UInt16
    dest_port::UInt16
    length::UInt16
    checksum::UInt16
    data::Vector{UInt8}
    UdpHdr() = new(0,0,0,0,Vector{UInt8}(undef, 0))
end # struct UdpHdr

mutable struct IcmpHdr
    ptype::UInt8
    code::UInt8
    checksum::UInt16
    identifier::UInt16
    seqno::UInt16
    IcmpHdr() = new(0,0,0,0,0)
end # struct IcmpHdr

mutable struct DecPkt
    datalink::EthHdr
	vlan::Any
    network::IpHdr
    protocol::Any
    DecPkt() = new(EthHdr(), VLAN(), IpHdr(), nothing)
end # struct DecPkt

@inline function getindex_he(::Type{T}, b::Vector{UInt8}, i) where {T}
    # When 0.4 support is dropped, add @boundscheck
    checkbounds(b, i + sizeof(T) - 1)
    return unsafe_load(Ptr{T}(pointer(b, i)))
end

@inline getindex_be(::Type{T}, b::Vector{UInt8}, i) where {T} = hton(getindex_he(T, b, i))

#----------
# decode ethernet header
#----------
function decode_eth_hdr(d::Array{UInt8})
	hex(n) = string(n, base = 16, pad = 2) # Added due to hex being deprecated
    eh = EthHdr()
    eh.dest_mac = join(hex.(d[1:6]), ":") # apply hex to first 6 elements and join result with ":"
    eh.src_mac  = join(hex.(d[7:12]), ":") # do it for next 6 elements
    eh.ptype    = getindex_be(UInt16, d, 13)
    eh
end # function decode_eth_hdr

#----------
# decode 802.1q VLAN
#----------
function decode_vlan(d::Array{UInt8})
	vlan = VLAN()
    vlan_temp::UInt16 = (0 + d[1]) << 8 + d[2] # converts first 2 bytes of UInt8 array to UInt16
	vlan.priority = (vlan_temp & 0xe000) >> 13
	vlan.dei = (vlan_temp & 0x1000) >> 12
	vlan.id = (vlan_temp & 0x0fff)
	vlan.ptype = getindex_be(UInt16, d, 3)
#	vlan.payload = d[4:end]
	vlan
end

#----------
# calculate IP checksum
#----------
function ip_checksum(buf::Array{UInt8})
    sum::UInt64 = 0
    for i in 1:2:length(buf)
        pair = getindex_he(UInt16, buf, i)
        sum += pair
        if (sum & 0x80000000) != 0
            sum = (sum & 0xFFFF) + (sum >> 16)
        end
    end

    while ((sum >> 16) != 0)
        sum = (sum & 0xFFFF) + (sum >> 16)
    end
    ~sum
end # function ip_checksum

#----------
# decode IP header
#----------
function decode_ip_hdr(d::Array{UInt8})
    iph = IpHdr()
    iph.version     = (d[1] & 0xf0) >> 4
    iph.length      = (d[1] & 0x0f) * 4
    if ip_checksum(d[1:iph.length]) == 0xFFFFFFFFFFFF0000
        iph.checksum = true
    end
    iph.services    = d[2]
    iph.totlen      = getindex_be(UInt16, d, 3)
    iph.id          = getindex_be(UInt16, d, 5)

    # set flags
    flags = IpFlags()
    flags.reserved   = (d[7] & (1 << 7)) > 0
    flags.dont_frag  = (d[7] & (1 << 6)) > 0
    flags.more_frags = (d[7] & (1 << 5)) > 0
    iph.flags        = flags

    iph.frag_offset = getindex_be(UInt16, d, 7) & 0x7ff
    iph.ttl         = d[9]
    iph.protocol    = d[10]
    iph.src_ip      = join(Int.(d[13:16]), ".")
    iph.dest_ip     = join(Int.(d[17:20]), ".")
    iph
end # function decode_ip_hdr

#----------
# decode TCP header
#----------
function decode_tcp_hdr(d::Array{UInt8})
    tcph = TcpHdr()
    tcph.src_port  = getindex_be(UInt16, d, 1)
    tcph.dest_port = getindex_be(UInt16, d, 3)
    tcph.seq       = getindex_be(UInt32, d, 5)
    tcph.ack       = getindex_be(UInt32, d, 9)
    tcph.offset    = (d[13] & 0xf0) >> 4

    # set flags
    flags = TcpFlags()
    flags.reserved = ((d[13] & 0x0e) >> 1) > 0
    flags.nonce    = (d[13] & 1) > 0
    flags.cwr      = (d[14] & (1 << 7)) > 0
    flags.ecn      = (d[14] & (1 << 6)) > 0
    flags.urgent   = (d[14] & (1 << 5)) > 0
    flags.ack      = (d[14] & (1 << 4)) > 0
    flags.push     = (d[14] & (1 << 3)) > 0
    flags.reset    = (d[14] & (1 << 2)) > 0
    flags.syn      = (d[14] & (1 << 1)) > 0
    flags.fin      = (d[14] & 1) > 0
    tcph.flags     = flags

    tcph.window    = getindex_be(UInt16, d, 15)
    tcph.checksum  = getindex_be(UInt16, d, 17)
    tcph.uptr      = getindex_be(UInt16, d, 19)
    tcph.data      = d[tcph.offset * 4 + 1:end]
    tcph
end # function decode_tcp_hdr

#----------
# decode UDP header
#----------
function decode_udp_hdr(d::Array{UInt8})
    udph = UdpHdr()
    udph.src_port  = getindex_be(UInt16, d, 1)
    udph.dest_port = getindex_be(UInt16, d, 3)
    udph.length    = getindex_be(UInt16, d, 5)
    udph.checksum  = getindex_be(UInt16, d, 7)
    udph.data      = d[9:end]
    udph
end # function decode_udp_hdr

#----------
# decode ICMP header
#----------
function decode_icmp_hdr(d::Array{UInt8})
    icmph = IcmpHdr()
    icmph.ptype      = d[1]
    icmph.code       = d[2]
    icmph.checksum   = getindex_be(UInt16, d, 3)
    icmph.identifier = getindex_be(UInt16, d, 5)
    icmph.seqno      = getindex_be(UInt16, d, 7)
    icmph
end # function decode_icmp_hdr

#----------
# decode ethernet packet
#----------
function decode_pkt(pkt::Array{UInt8})
    decoded           = DecPkt()
    decoded.datalink  = decode_eth_hdr(pkt)

	vlans = []
	byte_loc = 13
	temp_vlan_check = getindex_be(UInt16, pkt, byte_loc)
	while temp_vlan_check in Set([0x8100, 0x88a8, 0x9100])
		vlan = decode_vlan(pkt[byte_loc + 2:byte_loc + 5])
		push!(vlans, vlan)
		byte_loc += 4
		temp_vlan_check = getindex_be(UInt16, pkt, byte_loc)
	end

	decoded.vlan = vlans
	proto = nothing

    if temp_vlan_check == 0x0800
		iphdr = decode_ip_hdr(pkt[byte_loc + 2:byte_loc + 21])
		decoded.network = iphdr
		if (iphdr.protocol == 1)
        	proto = decode_icmp_hdr(pkt[byte_loc + 2 + iphdr.length:end])
    	elseif (iphdr.protocol == 6)
        	proto = decode_tcp_hdr(pkt[byte_loc + 2 + iphdr.length:end])
    	elseif (iphdr.protocol == 17)
        	proto = decode_udp_hdr(pkt[byte_loc + 2 + iphdr.length:end])
		end
	end

    decoded.protocol = proto
    decoded
end # function decode_pkt


text:
using Turing, BenchmarkTools, BenchmarkHelper
using LinearAlgebra

## Dummny benchmarks

BenchmarkSuite["dummy"] = BenchmarkGroup(["dummy"])

data = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

@model function constrained_test(obs)
    p ~ Beta(2,2)
    for i = 1:length(obs)
        obs[i] ~ Bernoulli(p)
    end
    p
end

BenchmarkSuite["dummy"]["dummy"] = @benchmarkable sample($(constrained_test(data)), $(HMC(0.01, 2)), 2000)

## gdemo

BenchmarkSuite["gdemo"] = BenchmarkGroup(["gdemo"])

@model function gdemo(x, y)
    s² ~ InverseGamma(2, 3)
    m ~ Normal(0, sqrt(s²))
    x ~ Normal(m, sqrt(s²))
    y ~ Normal(m, sqrt(s²))
    return s², m
end

BenchmarkSuite["gdemo"]["hmc"] = @benchmarkable sample($(gdemo(1.5, 2.0)), $(HMC(0.01, 2)), 2000)

##
## MvNormal
##

BenchmarkSuite["mnormal"] = BenchmarkGroup(["mnormal"])

# Define the target distribution and its gradient

@model function target(dim)
   Θ = Vector{Real}(undef, dim)
   θ ~ MvNormal(zeros(dim), I)
end

# Sampling parameter settings
dim = 10
n_samples = 100_000
n_adapts = 2_000

BenchmarkSuite["mnormal"]["hmc"] = @benchmarkable sample($(target(dim)), $(HMC(0.1, 5)), $n_samples)

## MvNormal: ForwardDiff vs BackwardDiff (Tracker)

@model function mdemo(d, N)
    Θ = Vector(undef, N)
   for n=1:N
      Θ[n] ~ d
   end
end

dim2 = 250
A    = rand(Wishart(dim2, Matrix{Float64}(I, dim2, dim2)));
d    = MvNormal(zeros(dim2), A)

# ForwardDiff
Turing.setadbackend(:forwarddiff)
BenchmarkSuite["mnormal"]["forwarddiff"] = @benchmarkable sample($(mdemo(d, 1)), $(HMC(0.1, 5)), 5000)

# BackwardDiff
Turing.setadbackend(:reversediff)
BenchmarkSuite["mnormal"]["reversediff"] = @benchmarkable sample($(mdemo(d, 1)), $(HMC(0.1, 5)), 5000)