text:
import NonUniformRandomVariateGeneration.sampleCategorical

# basic routines used by both smcSerial and smcParallel

# adapted from https://discourse.julialang.org/t/how-to-copy-all-fields-without-changing-the-referece/945/5
@inline @generated function particleCopy!(dest, src)
  fieldNames = fieldnames(dest)
  fieldTypes = dest.types
  numFields = length(fieldNames)
  expressions = Array{Expr}(undef, numFields)

  for i = 1:numFields
    fieldName = fieldNames[i]
    fieldType = fieldTypes[i]
    type_mutable = false
    if VERSION < v"1.7"
      type_mutable = fieldType.mutable
    else
      type_mutable = ismutabletype(fieldType)
    end

    @assert !type_mutable || hasmethod(copyto!, (fieldType,
      fieldType)) "$fieldName::$fieldType : copyto! must exist for mutable Particle fields"
    if type_mutable
      @inbounds expressions[i] = :(copyto!(dest.$fieldName, src.$fieldName))
    else
      @inbounds expressions[i] = :(dest.$fieldName = src.$fieldName)
    end
  end
  body = Expr(:block, expressions...)

  quote
    # $(Expr(:meta, :inline))
    $body
    return
  end
end

## better not to inline ; enough work should be done in the loop
function _copyParticles!(out::Vec, in::Vec) where {Particle,
  Vec<:AbstractVector{Particle}}
  for i in eachindex(out)
    @inbounds particleCopy!(out[i], in[i])
  end
end

## better not to inline ; enough work should be done in the loop
function _copyParticles!(zetaAncs::Vec1, zetas::Vector{Particle},
  as::Vec2) where {Particle, Vec1<:AbstractVector{Particle},
    Vec2<:AbstractVector{Int64}}
  for i in eachindex(zetaAncs)
    @inbounds particleCopy!(zetaAncs[i], zetas[as[i]])
  end
end

## better not to inline ; enough work should be done in the loop
function _setEves!(eves::Vec1, oldEves::Vec2, as::Vec1) where
  {Vec1<:AbstractVector{Int64}, Vec2<:AbstractVector{Int64}}
  for i in eachindex(eves)
    @inbounds eves[i] = oldEves[as[i]]
  end
end

## better not to inline ; enough work should be done in the loop
function _mutateParticles!(zetas::Vec, engine::RNG, p::Int64, M!::F,
  zetaAncs::Vec, pScratch::ParticleScratch) where {Particle,
  Vec<:AbstractVector{Particle}, F<:Function, ParticleScratch}
  for j in eachindex(zetas)
    @inbounds M!(zetas[j], engine, p, zetaAncs[j], pScratch)
  end
end

## better not to inline ; enough work should be done in the loop
function _mutateParticles!(zetas::Vec, engine::RNG, p::Int64, M!::F,
  zetaAncs::Vec, pScratch::ParticleScratch, xref::Particle) where {Particle,
  Vec<:AbstractVector{Particle}, F<:Function, ParticleScratch}
  particleCopy!(zetas[1], xref)
  for j in 2:length(zetas)
    @inbounds M!(zetas[j], engine, p, zetaAncs[j], pScratch)
  end
end

## better not to inline ; enough work should be done in the loop
function _logWeightParticles!(lws::Vec1, p::Int64, lG::F,
  zetas::Vec2, pScratch::ParticleScratch, resample::Bool, oldWs::Vec1) where
  {Vec1<:AbstractVector{Float64}, F<:Function, Particle,
  Vec2<:AbstractVector{Particle}, ParticleScratch}
  for j in eachindex(zetas)
    @inbounds lws[j] = lG(p, zetas[j], pScratch)
  end
  !resample && (lws .+= log.(oldWs))
end

@inline function pickParticle!(path::Vector{Particle}, smcio::SMCIO{Particle}) where Particle
  @assert smcio.fullOutput
  @assert length(path) == smcio.n
  n::Int64 = smcio.n
  allAs::Vector{Vector{Int64}} = smcio.allAs
  allZetas::Vector{Vector{Particle}} = smcio.allZetas
  rng::RNG = getRNG()
  k::Int64 = sampleCategorical(smcio.ws, rng)

  @inbounds particleCopy!(path[n], allZetas[n][k])
  for p = n-1:-1:1
    @inbounds k = allAs[p][k]
    @inbounds particleCopy!(path[p], allZetas[p][k])
  end
  return
end

@inline function pickParticleBS!(path::Vector{Particle},
  smcio::SMCIO{Particle, ParticleScratch}, lM::F) where {Particle,
  ParticleScratch, F<:Function}
  @assert smcio.fullOutput
  @assert length(path) == smcio.n
  n::Int64 = smcio.n
  N::Int64 = smcio.N
  allAs::Vector{Vector{Int64}} = smcio.allAs
  allZetas::Vector{Vector{Particle}} = smcio.allZetas
  allWs::Vector{Vector{Float64}} = smcio.allWs
  bws::Vector{Float64} = smcio.internal.scratch1
  pScratch::ParticleScratch = smcio.internal.particleScratch

  rng = getRNG()
  k::Int64 = sampleCategorical(smcio.ws, rng)

  @inbounds particleCopy!(path[n], allZetas[n][k])
  for p = n-1:-1:1
    @inbounds bws .= log.(allWs[p])
    for j in 1:N
      @inbounds bws[j] += lM(p+1, allZetas[p][j], path[p+1], pScratch)
    end
    m::Float64 = maximum(bws)
    bws .= exp.(bws .- m)
    k = sampleCategorical(bws, rng)
    @inbounds particleCopy!(path[p], allZetas[p][k])
  end
end


text:
using SystemBenchmark
using CSV
using DataFrames
using Plots
using Statistics
gr()

function make_relative!(df)
    for col in 12:size(df,2)
        df[!,col] = df[!,col] ./ df[1,col]
    end

    df[!,:mean_cpu] = map(x->mean([x.FloatMul,
        # x.FusedMulAdd, # results unreliable before 0.2.1
        x.FloatSin,
        x.VecMulBroad,
        #x.CPUMatMul,
        x.MatMulBroad,
        x[Symbol("3DMulBroad")]]),eachrow(df))

    df[!,:mean_diskio] = map(x->mean([x.DiskWrite1KB,
        x.DiskWrite1MB,
        x.DiskRead1KB,
        x.DiskRead1MB]),eachrow(df))
    return df
end

function plotreport(df, figurepath; scale=2.0)
    platforms = ["Windows", "macOS", "Linux (x86", "Linux (aarch"]
    colors = [:blue,:orange,:green,:purple]

    p1 = plot(dpi=300)
    plot!(0:100,0:100,color=:gray)
    i = 1
    for plat = platforms
        df2 = df[occursin.(plat,df.OS),:]
        scatter!(df2[!,:mean_cpu], df2[!,:compilecache], label=replace(plat,"("=>""), leg=false,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time (relative to ref)")
        ylabel!("compilecache time\n(relative to ref)")
        i += 1
    end
    xlims!(0,7); ylims!(0,7)

    p2 = plot()
    plot!(0:100,0:100,color=:gray,label="parity")
    i = 1
    for plat in platforms
        df2 = df[occursin.(plat,df.OS),:]
        scatter!(df2[!,:mean_cpu],df2[!,:JuliaLoad],label=replace(plat,"("=>""), legend=:best, bg_legend = :transparent, fg_legend = :transparent,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time (relative to ref)")
        ylabel!("Julia startup time\n(relative to ref)")
        i += 1
    end
    xlims!(0,7); ylims!(0,7)

    p3 = plot()
    plot!(0:100,0:100,color=:gray)
    i = 1
    for plat in platforms
        df2 = df[occursin.(plat,df.OS),:]
        scatter!(df2[!,:mean_cpu],df2[!,:mean_diskio],label=replace(plat,"("=>""), fontsize=4,leg=false,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time (relative to ref)")
        ylabel!("Mean disk IO time\n(relative to ref)")
        i += 1
    end
    xlims!(0,7); ylims!(0,45)

    p4 = plot()
    plot!(0:100,0:100,color=:gray)
    i = 1
    for plat in platforms
        df2 = df[occursin.(plat,df.OS),:]
        scatter!(df2[!,Symbol("compilecache")],df2[!,Symbol("FFMPEGH264Write")],label=replace(plat,"("=>""), fontsize=4,leg=false,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("compilecache time\n(relative to ref)")
        ylabel!("FFMPEGH264Write time\n(relative to ref)")
        i += 1
    end
    xlims!(0,12); ylims!(0,12)

    plot(p1,p2,p3,p4,dpi=300,size=(400*scale,300*scale))

    savefig(figurepath)
end

function memoryreport(df_in, figurepath; scale=2.0)
    df = deepcopy(df_in)
    platforms = ["Windows", "macOS", "Linux (x86", "Linux (aarch"]
    colors = [:blue,:orange,:green,:purple]

    p1 = plot(dpi=300)
    i = 1
    for plat = platforms
        df2 = df[occursin.(plat,df.OS),:]
        x = df2[!,:mean_cpu]
        y = df2[!,:Bandwidth10kB]
        x = x[.!ismissing.(y)]
        y = y[.!ismissing.(y)]
        scatter!(x, y, label=replace(plat,"("=>""), leg=false,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time")
        ylabel!("Bandwidth10kB")
        i += 1
    end
    #xlims!(0,10); ylims!(0,15)

    p2 = plot()
    i = 1
    for plat in platforms
        df2 = df[occursin.(plat,df.OS),:]
        x = df2[!,:mean_cpu]
        y = df2[!,:Bandwidth100kB]
        x = x[.!ismissing.(y)]
        y = y[.!ismissing.(y)]
        scatter!(x, y, label=replace(plat,"("=>""), legend=:best, bg_legend = :transparent, fg_legend = :transparent,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time")
        ylabel!("Bandwidth100kB")
        i += 1
    end
    #xlims!(0,10); ylims!(0,15)

    # p3 = plot()
    # i = 1
    # for plat in platforms
    #     df2 = df[occursin.(plat,df.OS),:]
    #     x = df2[!,:mean_cpu]
    #     y = df2[!,:Bandwidth1MB]
    #     x = x[.!ismissing.(y)]
    #     y = y[.!ismissing.(y)]
    #     scatter!(x, y, label=replace(plat,"("=>""), leg=false,ms=4,markerstrokewidth=0,color=colors[i])
    #     xlabel!("Mean CPU time")
    #     ylabel!("Bandwidth1MB")
    #     i += 1
    # end

    p3 = plot()
    i = 1
    for plat in platforms
        df2 = df[occursin.(plat,df.OS),:]
        x = df2[!,:mean_cpu]
        y = df2[!,:Bandwidth10MB]
        x = x[.!ismissing.(y)]
        y = y[.!ismissing.(y)]
        scatter!(x, y, label=replace(plat,"("=>""), leg=false,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time")
        ylabel!("Bandwidth10MB")
        i += 1
    end

    p4 = plot()
    i = 1
    for plat in platforms
        df2 = df[occursin.(plat,df.OS),:]
        x = df2[!,:mean_cpu]
        y = df2[!,:DeepCopy]
        x = x[.!ismissing.(y)]
        y = y[.!ismissing.(y)]
        scatter!(x, y, label=replace(plat,"("=>""), leg=false,ms=4,markerstrokewidth=0,color=colors[i])
        xlabel!("Mean CPU time")
        ylabel!("DeepCopy")
        i += 1
    end

    plot(p1,p2,p3,p4,dpi=300,size=(400*scale,300*scale))

    savefig(figurepath)
end

using SystemBenchmark
df = getsubmittedbenchmarks()
make_relative!(df)
savebenchmark(joinpath(@__DIR__,"all.csv"), df)
plotreport(df, joinpath(@__DIR__,"summary_report.png"))
memoryreport(df, joinpath(@__DIR__,"memory_report.png"))

ext: jl


text:
#
# generator macro that translates the body of a generator to a Julia implementation.
#
# example syntax:
#
# @generator GenName(subgen1, subgen2, ...) begin
#   generates: ["..", ".."]
#		method1 = begin ... end
#		method2(param) = begin ... end
#		function method3()
#      ...
#		end
# end
#
# note: we generally want to avoid the standard hygiene changes performed by the macro expander
# as the rule code should execute in the context in which the macro is called.  However, code added
# to handle the constructs (such as choice points) should operate in the context of the DataGenerators
# module.
# The ideal approach would be to apply hygiene to DataGenerators code, and then escape any user-defined rule code.  However
# since DataGenerators code is mixed very closely in AST with user code (for instance when a GT construct appears in
# argument values passed to a function), this difficult to achieve correctly. For example, escaping all nodes
# apart from those that are GT construct will work in most cases, but fails if the user defines an anonymous
# function such i->... (e.g. in a call to map(..., )) as the escaped variable gives rise to an error (invalid
# assignment target)
#
#
#
macro generator(gensig, genbody)

	# macro is passed:
	# (1) the 'signature' of the generator: generator name and (optionally) subgenerator arguments
	# (2) the 'body' of the generator: a block containing methods and metadata

	# extract and check the name of the generator and (zero or more) subgen parameters
	genname, subgenargs = extractfuncsig(gensig, true) # true: allows permits function signature without parenthesis when there are no parameters

	# extractfuncsig returns nothing if there is a problem
	if genname == nothing
		error("The generator name and arguments are not valid: $(gensig)")
	end

	# check that all arguments are simple symbols: currently do not support typed arguments
	if any(arg -> (typeof(arg) != Symbol), subgenargs)
		error("Not all arguments to the generator are valid - remove any type specifications: $(genargs)")
	end

	# extract metadata and rules from the body of the generator
	metadata, generatorrules = extractmetadataandrules(genbody)

	# rewrite rules as functions and transform constructs
	gencontext = GeneratorContext(genname, subgenargs)
	transformrules!(generatorrules, gencontext)

	# construct and return code for generator type and methods that implement the rules
	typeblock = constructtype(genname, subgenargs, metadata, gencontext)
	methodsblock = constructmethods(generatorrules)

	tidyfinish = quote $(genname); end # avoid outputting the last method's name as result of macro; output type instead

	typeandmethods = mergeexprs(typeblock, methodsblock)

	esc(mergeexprs(typeandmethods, tidyfinish))

end

#
# context during rule transformation (e.g. generator name, choice points)
#
type GeneratorContext
	# Random offset for all choice point numbers for this generator. A random offset is used
	# to make it very unlikely that choice points from different generators would have the
	# same number.
	choicepointoffset::UInt
	numchoicepoints::UInt
	choicepointinfo::Dict{UInt, Dict{Symbol, Any}} # Map from choice point num to info dict
	genname::Symbol
	subgenargs::Vector{Symbol}
	rulemethodnames::Dict{Symbol, Symbol} # mapping of user rule name to actual methods name to avoid issues when methods of same name already exist and can't be extended
	genparam::Symbol
	genarg::Expr
	stateparam::Symbol
	statearg::Expr
	function GeneratorContext(genname, subgenargs)
		choicepointoffset = rand63bitint() # all choice points numbers will come in sequence after this number
		genparam = gensym(:g)
		stateparam = gensym(:s)
		genarg = :( $(genparam)::$(genname) )
		statearg = :( $(stateparam)::$(THIS_MODULE).DerivationState )
		new(choicepointoffset, 0, Dict{Int, Dict{Symbol, Any}}(), genname, subgenargs, Dict{Symbol, Symbol}(), genparam, genarg, stateparam, statearg)
	end
end

#
# holds a generator rule during transformation
#
type GeneratorRule
	rulename::Symbol
	args::Any
	body::Any
	methodname::Any
	function GeneratorRule(rulename, args, body, methodname=nothing)
		new(rulename, args, body, methodname)
	end
end

#
# extract metadata and rules from the body of the generator
#
# metadata syntax:
#		tag: value
#
# rule syntax:
#		method = begin ... end
#		method(param) = begin ... end
#		function method() ... end
#
function extractmetadataandrules(genbody::Expr)

	metadata = Dict{Symbol, Any}()  	# key is tag, value is value
	generatorrules = Vector{GeneratorRule}()

	extractmetadataandrulesfromblock(genbody, metadata, generatorrules)

	metadata, generatorrules

end

#
# extract metadata and rules from the body of the generator
# separate from extractmetadataandrules function above so it may be called recursively should the
# metadata and rules be enclosed in further begin ... end blocks at the top level
#
function extractmetadataandrulesfromblock(block::Expr, metadata::Dict{Symbol, Any}, generatorrules::Vector{GeneratorRule})

	# metadata and rules are defined at the top level
	for node in removelinenodes(block.args)

		if typeof(node) == Expr

			if node.head == :block
				extractmetadataandrulesfromblock(node, metadata, generatorrules)
				continue
			end

			# metadata syntax tree pattern
			#		:(:)
			#			:key
			#			expr
			if node.head == :(:)
				colonargs = removelinenodes(node.args)
				if length(colonargs) == 2 && typeof(colonargs[1]) == Symbol
					metadata[colonargs[1]] = eval(colonargs[2])  # TODO could eval have unintended side effects?
					continue
				end
			end

			# function call syntax tree pattern tested by extractfuncdef - it returns nothing if not a function definition
			(rulename, ruleargs, rulebody) = extractfuncdef(node)
			if rulename != nothing
				push!(generatorrules, GeneratorRule(rulename, ruleargs, rulebody))
				continue
			end

			# TODO recursively permit blocks within blocks to still be 'top level'
		end

		error("Unrecognised statement at the top level of the generator body: $(node)")

	end

end

function constructtype(genname, subgenargs, metaInfo, gencontext::GeneratorContext)

	# code for the generator type
	# note that macro hygeine will ensure that variables/functions/types not explicitly esc'ed will be transformed
	# into the context of the DataGenerators module, which is what we require here

	# we store also the current module at the time of calling this macro as the owning module of this new type:
	# this is so that rules are executed in the correct context when the generator is run since the same context
	# as this type is where the method corresponding to generator rules will be created.
	# it would be possible to derive the owning module from the fully specified type of the generator,
	# but there is currently no built-in Julia function to do this cleanly; instead we would need to
	# perform some custom string-manipulation on the type name, and this wouldn't be very robust (e.g.
	# may not survive type-aliasing etc.)

	# Further note: Module is stored as a symbol in the type since deepcopy does not support fields of type Module,
	# and to recreate the Module type, we also need its parent

	quote
		type $(genname) <: $(THIS_MODULE).Generator
			meta::Dict{Symbol, Any}
			choicepointinfo::Dict{UInt, Dict{Symbol, Any}}
			rulemethodnames::Dict{Symbol,Symbol}
			subgens::Vector{$(THIS_MODULE).Generator}
			choicemodel::$(THIS_MODULE).ChoiceModel
			evalfn::Function
			function $(genname)(subgens::Vector = [])
				if length(subgens) != $(length(subgenargs))
					error("Incorrect number of sub generators $(length(subgens))")
				end

				if !all([typeof(sg) <: $(THIS_MODULE).Generator for sg in subgens])
					error("Not all subgenerators are of type $(THIS_MODULE).Generator $(subgens)")
				end
				cpinfo = $(gencontext.choicepointinfo)
				new($metaInfo, cpinfo, $(gencontext.rulemethodnames), subgens, $(THIS_MODULE).SamplerChoiceModel(cpinfo),
				 		ex->eval($(current_module()),ex))
			end

			$(genname)(subgens...) = $(genname)(collect(subgens))
		end
	end

end

function constructmethods(generatorrules::Vector{GeneratorRule})

	# construct methods for generator rules (as short form function since body is guaranteed to be a block)

	methodexprs = map(gm -> Expr(:(=), Expr(:call, gm.methodname, gm.args...), gm.body), generatorrules)
	Expr(:block, methodexprs...)

end

#
# main entry point -- transforms rules written with Data Generator constructs into standard Julia methods
#
function transformrules!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)

	# (1) choose(Type) where Type is a not a directly supported type, nor generator
	transformchoosetypes!(genrules, gencontext)

	# (2) choose(<: String, [regex])
	transformchoosestrings!(genrules, gencontext)

	# up to here is processing of syntactic sugar: expansion adds rules that use Data Generators constructs
	# (which will be processed in subsequent steps, so this sugar processing must occur first)
	# however, subsequent transforms change Data Generators constructs into standard Julia

	# (3) choose(<: Number, [min, [max]])
	transformchoosenumbers!(genrules, gencontext)

	# (4) choose(Gen) where Gen is a generator type (means include choice points into this generator)
	# transformchoosegentypes!(genrules, gencontext)

	# (5) choose(subgen) where subgen is a generator instance (means call gen but don't include choice points)
	transformchoosesubgens!(genrules, gencontext)

	# (6) sequence choice point
	transformsequencechoicepoints!(genrules, gencontext)

    # (7) transform signatures (method names and arguments)
	transformsignatures!(genrules, gencontext)

	# (8) create umbrella methods, and transform rule choice points
	createumbrellamethods!(genrules, gencontext)

	# (9) calls to other rules
	transformrulecalls!(genrules, gencontext)

	# (10) indirect calls to other rules - feature used to avoid long-running compilation as a result of (seemingly) recursive type inference
	transformindirectrulecalls!(genrules, gencontext)

	# (11) direct call to subgenerators (no longer supported)
	# transformdirectsubgencalls!(genrules, gencontext)

end

#
# transform all type choice points: choose(Type)
#
function transformchoosetypes!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, ischoosetype, transformchoosetype)
	end
end

# is a call to choose(Type) point
function ischoosetype(callname, callparams, gencontext)
	if (callname == :choose)
		if (length(callparams) >= 1)
			datatype = converttodatatype(callparams[1])
			# type must not be a directly supported type (e.g. Int64, String) nor a Generator
			if (typeof(datatype) <: Type) && !(datatype in GENERATOR_SUPPORTED_CHOOSE_TYPES) && !(datatype <: Generator)
				return true
			end
		end
	end
	false
end

# expand Type choice points of the using translator for Types
function transformchoosetype(callname, callparams, genrules, gencontext, matchfn, transformfn)

	datatype = converttodatatype(callparams[1])

	if length(callparams) > 2
		error("choose($(datatype)) must have no further parameters")
	end

	# attempt to make rulenames unique for additional type rules by adding a prefix containing a random 64-bit values
	# in particular, it is unlikely any user-defined (or auto-generated) rules will have names the same as rule generated
	# for this type
	rulenameprefix = "choosetype" * hex(rand(UInt64))
	# create rules for regex
	typerules = type_rules(datatype, rulenameprefix)
	# add rules to generator
	addrulesources(genrules, typerules)
	# replace choose with call to entry point for regex rules (the first rule)
	:( $(typerules[1].rulename)() )

end

#
# transform all string choice points: choose(StringType, [regex])
#
function transformchoosestrings!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, ischoosestring, transformchoosestring)
	end
end

# is a call to choose(String) point
function ischoosestring(callname, callparams, gencontext)
	if (callname == :choose)
		if (length(callparams) >= 1)
			datatype = converttodatatype(callparams[1])
			if (typeof(datatype) <: Type) && (datatype in GENERATOR_SUPPORTED_CHOOSE_STRING_TYPES)
				return true
			end
		end
	end
	false
end

# expand String choice points of the form:
#		choose(type,[regex])
# by creating additional rules to emit strings satisfying the regex
function transformchoosestring(callname, callparams, genrules, gencontext, matchfn, transformfn)

	datatype = converttodatatype(callparams[1])

	if length(callparams) > 2
		error("choose($(datatype),...) must have at most one further parameter")
	end

	regex = "" # interpreted as wildcard
	if length(callparams) >= 2
		if !(typeof(regex) <: AbstractString)
			error("regex in choose($(datatype),...) must be a literal string")
		end
		regex = callparams[2]
	end

	# attempt to make rulenames unique for additional regex rules by adding a prefix containing a random 64-bit values
	# in particular, it is unlikely any user-defined (or auto-generated) rules will have names the same as rule generated
	# for this regex
	rulenameprefix = "choosestring" * hex(rand(UInt64))
	# create rules for regex
	regexrules = DataGeneratorTranslators.regex_rules(regex, datatype, rulenameprefix)
	# add rules to generator
	addrulesources(genrules, regexrules)
	# replace choose with call to entry point for regex rules (the first rule)
	:( $(regexrules[1].rulename)() )

end

#
# transform value choice points with simple numeric (inc. Boolean) types
#

function transformchoosenumbers!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, ischoosenumber, transformchoosenumber)
	end
end

# is a call to value choice point
function ischoosenumber(callname, callparams, gencontext)
	if (callname == :choose)
		if (length(callparams) >= 1)
			datatype = converttodatatype(callparams[1])
			if (typeof(datatype) <: Type) && (datatype in GENERATOR_SUPPORTED_CHOOSE_NUMBER_TYPES)
				return true
			end
		end
	end
	false
end

# transform numeric value choice points of the form:
#		choose(type, [min, [max]])
# (where parameters after the datatype constrain the range of the type)
# to:
#		DataGenerators.choosenumber(s, cpid, datatype, minval, maxval, rangeisliteral)
function transformchoosenumber(callname, callparams, genrules, gencontext, matchfn, transformfn)

	datatype = converttodatatype(callparams[1])
	# thanks to ischoosenumber matching function, we know we that there is at least one param, and that it is a symbol

	if datatype <: Bool

		# follows same pattern as other 'numeric' types, which is possible since 0~false 1~true
		# we don't allow any possibility to restrict this range

		if length(callparams) > 1
			error("choose($(datatype)) must have no further parameters")
		end

		minval = false
		maxval = true
		rangeisliteral = true
		cpinfo = Dict{Symbol,Any}(:datatype=>datatype, :min=>minval, :max=>maxval)
		cpid = recordchoicepoint(gencontext, :value, cpinfo)
		chooseexpr = :( $(THIS_MODULE).choosenumber($(gencontext.stateparam), $(cpid), $(datatype), $(minval), $(maxval), $(rangeisliteral)) )

	else

		# here parameters can be used to define a range of possible values
		# this can either be done via literals or expression - in the former case, the literal values are recorded and passed to the choice model
		# since knowing the bound(s) of the valid range can enable a better model than one that must potentially varying ranges

		if length(callparams) > 3
			error("choose($(datatype),...) must have at most two further parameters")
		end

		# parameters themselves could be nested choose(Number, ...) expression, so transform these recursively:
		for i in 2:length(callparams)
			callparams[i] = transformfunccall(callparams[i], genrules, gencontext, matchfn, transformfn)
		end

		if length(callparams) >= 2
			minval, minisliteral = processpossiblyliteralparam(callparams[2], datatype, gencontext)
		else
			minval, minisliteral = typemin(datatype), true
		end

		if length(callparams) >= 3
			maxval, maxisliteral = processpossiblyliteralparam(callparams[3], datatype, gencontext)
		else
			maxval, maxisliteral = typemax(datatype), true
		end

		cpinfo = Dict{Symbol,Any}(:datatype=>datatype)

		# record literal values in choice point info as an indicator to choice model that limit on range will not change
		if minisliteral
			cpinfo[:min] = minval
		end
		if maxisliteral
			cpinfo[:max] = maxval
		end

		# rangeisliteral parameter will avoid a further runtime check on type validity if both limits are literal
		rangeisliteral = minisliteral && maxisliteral

		cpid = recordchoicepoint(gencontext, :value, cpinfo)
		chooseexpr = :( $(THIS_MODULE).choosenumber($(gencontext.stateparam), $(cpid), $(datatype), $(minval), $(maxval), $(rangeisliteral)) )

	end

	chooseexpr

end

#
# transform choose(subgen)
#

function transformchoosesubgens!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, ischoosesubgen, transformchoosesubgen)
	end
end

# is a call to value choice point
ischoosesubgen(callname, callparams, gencontext) = (callname == :choose) && (length(callparams) >= 1) && (callparams[1] in gencontext.subgenargs)

# call to a sub-generator becomes: DataGenerators.subgen(g, s, i) where i it the index of the sub-generators in the arguments
function transformchoosesubgen(callname, callparams, genrules, gencontext, matchfn, transformfn)
	if length(callparams) > 1
		error("choose($(callparams[1])) must have no further parameters")
	end
	i = findfirst(gencontext.subgenargs, callparams[1])
	:( $(THIS_MODULE).subgen($(gencontext.genparam), $(gencontext.stateparam), $(i)) )
end

#
# transform all sequence choice points
#

function transformsequencechoicepoints!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, issequencechoicepoint, transformsequencechoicepoint)
	end
end

# is a call to sequence choice point
issequencechoicepoint(callname, callparams, gencontext) = (callname in [:mult, :plus, :reps,])

# transform sequence choice points of the form:
#		reps(:rule, min, max)
#		mult(:rule)
#		plus(:rule)
#	to:
#		DataGenerators.choosereps(s, cpid, ()->rule(g,s), minreps, maxreps, rangeisliteral)
function transformsequencechoicepoint(callname, callparams, genrules, gencontext, matchfn, transformfn)

	if length(callparams) < 1
		error("$(callname) must specify an expression to sequence")
	end

	# allow for nested sequence choice points by transforming the expression to sequence
	functocallexpr = transformfunccall(callparams[1], genrules, gencontext, matchfn, transformfn)

	# synatic sugar: intepret a symbol of a rule name as a no-parameter call to that rule
	functocallexpr = convertrulenametorulecall(functocallexpr, genrules)

	if callname == :mult

		if length(callparams) > 1
			error("$(callname) must have no parameters other than function to call")
		end

		minreps, minisliteral = 0, true
		maxreps, maxisliteral = typemax(Int), true

	elseif callname == :plus

		if length(callparams) > 1
			error("$(callname) must have no parameters other than function to call")
		end

		minreps, minisliteral = 1, true
		maxreps, maxisliteral = typemax(Int), true

	elseif callname == :reps

		if length(callparams) > 3
			error("$(callname) must have at most two parameters other than function to call")
		end

		if length(callparams) >= 2
			minreps, minisliteral = processpossiblyliteralparam(callparams[2], Int, gencontext)
		else
			minreps, minisliteral = 0, true
		end

		if length(callparams) >= 3
			maxreps, maxisliteral = processpossiblyliteralparam(callparams[3], Int, gencontext)
		else
			maxreps, maxisliteral = typemax(Int), true
		end

	else

		error("unrecognised sequence choice point construct $(callname)") # shouldn't happen

	end

	cpinfo = Dict{Symbol, Any}()

	# record literal values in choice point info as an indicator to choice model that limit on range will not change
	if minisliteral
		cpinfo[:min] = minreps
	end
	if maxisliteral
		cpinfo[:max] = maxreps
	end

	# rangeisliteral parameter will avoid a further runtime check on type validity if both limits are literal
	rangeisliteral = minisliteral && maxisliteral

	cpid = recordchoicepoint(gencontext, :sequence, cpinfo)
	idxvar = gensym("idx")
	:( [ $(functocallexpr) for $(idxvar) in 1:($(THIS_MODULE).choosereps($(gencontext.stateparam), $(cpid), $(minreps), $(maxreps), $(rangeisliteral))) ] )

end

#
# transform rule signatures (name and args)
# rule name is converted to a unique method name; internal arguments are added
#

function transformsignatures!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)

	for genrule in genrules

		genrule.methodname = uniquemethodname(genrule.rulename)

		# if rule body is not already a block, make it so (so can be written as short form function)
		if (typeof(genrule.body) != Expr) || (genrule.body.head != :block)
			genrule.body = Expr(:block, genrule.body)
		end

		genrule.args = [gencontext.genarg; gencontext.statearg; genrule.args]

	end

end

#
# create an 'umbrella' method for each unique rule name
# the umbrella method is the entry point for the rule, and records rule entry and exit during generation
# for rule choice points (where there is more than one rule with the same name), it handles the selection of which rule is called
#

function createumbrellamethods!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)

	# group rules with the same name together (so we can process rule choice constructs)
	groupedrules = Dict{Symbol, Vector{GeneratorRule}}()
	for genrule in genrules
		rules = get!(groupedrules, genrule.rulename, Vector{GeneratorRule}())
		push!(rules, genrule)
		# TODO: check that arguments also match?
	end

	# TODO sort(collect(rules)) to ensure predictable order?
	for (rulename, rules) in groupedrules

		# assign a unique method name
		umbrellamethodname = uniquemethodname(rulename)
		# and store this since it is the method that must be called when the rule is executed
		gencontext.rulemethodnames[rulename] = umbrellamethodname

		#
		# call (one of) the method(s) for the rule(s)
		# if there is more than one rule, then create a rule choice point that determines which to call
		#

		# assume that all individual rules have the same set of arguments, and take these from the first rule
		umbrellaargs = rules[1].args

		# extract parameters from these args for calls made within the new function to the rule methods
		umbrellaparams = [extractparamfromarg(arg) for arg in umbrellaargs]

		if length(rules) > 1

			# if more than one rule with this name, then a rule choice point is required

			# ensure unique name for chosenidx in case of other variable defined in context
			chosenidxvar = gensym("chosenidx")

			condexpr = nothing
			for idx in 1:length(rules)

				# call method with new name
				callexpr = Expr(:call, rules[idx].methodname, umbrellaparams...)

				# build conditional expr
				if condexpr == nothing
					condexpr = callexpr
				else
					condexpr = :( ($(chosenidxvar) == $(idx)) ? $(callexpr) : ($condexpr) )
				end

			end

			cpinfo = Dict{Symbol,Any}(:rulename => rulename, :min => 1, :max => length(rules))
			cpid = recordchoicepoint(gencontext, :rule, cpinfo)

			umbrellabody = Expr(:block, :( $(chosenidxvar) = $(THIS_MODULE).chooserule($(gencontext.stateparam), $(cpid), $(length(rules))) ), condexpr)

		else

			# if only one method for this rule, simply call the method

			callexpr = Expr(:call, rules[1].methodname, umbrellaparams...)

			umbrellabody = Expr(:block, callexpr)

		end

		#
		# now 'wrap' rule body with calls to indicate the start and end of the rule:
		#
		# begin
		#	recordstartofrule(state, rulemethodname)
		# 	result = <rulebody block created in code above>
		#	recordendofrule(state, rulemethodname)
		#	result
		# end
		#
		# note: we use the umbrella method name (which might be somewhat obfuscated) rather than the rulename,
		# since the rulename need not be unique across a generator and its subgenerators

		resultvar = gensym("result")
		rulenameexpr = QuoteNode(umbrellamethodname)
		umbrellabody = Expr(:block,
			:( $(THIS_MODULE).recordstartofrule($(gencontext.stateparam), $(rulenameexpr)) ),
			:( $(resultvar) = $(umbrellabody) ),
			:( $(THIS_MODULE).recordendofrule($(gencontext.stateparam)) ),
			:( $(resultvar) )
		)

		push!(genrules, GeneratorRule(rulename, umbrellaargs, umbrellabody, umbrellamethodname))

	end

end

#
# transform all call to rules by calling the umbrella method
#

function transformrulecalls!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, isrulecall, transformrulecall)
	end
end

# is a call to a generator rule if called function is a key to the list of rulemethods
isrulecall(callname, callparams, gencontext) = haskey(gencontext.rulemethodnames, callname)

# to call a method, call the umbrella method
function transformrulecall(callname, callparams, genrules, gencontext, matchfn, transformfn)
	umbrellamethodname = gencontext.rulemethodnames[callname]
	rewrittenparams = [gencontext.genparam; gencontext.stateparam; callparams]
	Expr(:call, umbrellamethodname, rewrittenparams...)
end

#
# transform all indirect call to rules by calling the umbrella method
#

function transformindirectrulecalls!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, isindirectrulecall, transformindirectrulecall)
	end
end

# is a call to a generator rule if called function is a key to the list of rulemethods
isindirectrulecall(callname, callparams, gencontext) = (callname == :rule) && length(callparams) >= 1

# to call a method, call a helper method in generation code
function transformindirectrulecall(callname, callparams, genrules, gencontext, matchfn, transformfn)
	rewrittenparams = [gencontext.genparam; gencontext.stateparam; callparams]
	:( $(THIS_MODULE).indirectrulecall($(rewrittenparams...)) )
end

#
# transform all direct calls to rules by calling the umbrella method (deprecated)
#

function transformdirectsubgencalls!(genrules::Vector{GeneratorRule}, gencontext::GeneratorContext)
	for genrule in genrules
		genrule.body = transformfunccall(genrule.body, genrules, gencontext, isdirectsubgencall, transformdirectsubgencall)
	end
end

# is a call to a subgen of the generator
isdirectsubgencall(callname, callparams, gencontext) = (callname in gencontext.subgenargs)

# call to a sub-generator becomes: DataGenerators.subgen(g, s, i)
# where i it the index of the sub-generators in the arguments
function transformdirectsubgencall(callname, callparams, gencontext, matchfn, transformfn)
	warn("Direct calls to sub-generators are deprecated; please use: choose(<subgenname>)")
	i = findfirst(gencontext.subgenargs, callname)
	:( $(THIS_MODULE).subgen($(gencontext.genparam), $(gencontext.stateparam), $(i)) )
end

#
# Utility methods
#

# We use 63 bit offset to ensure there is room to add choice points after the offset
# value. If we used 64 bit we might get too close to the end of the range...
function rand63bitint()
	candidate = rand(UInt64)
	while candidate > (2^63-1)
		candidate = rand(UInt64)
	end
	candidate
end

function recordchoicepoint(gencontext::GeneratorContext, cptype::Symbol, cpinfo::Dict)
	cpid = nextchoicepointnum(gencontext)
	cpinfo[:type] = cptype
	gencontext.choicepointinfo[cpid] = cpinfo
	cpid
end

# get a unique choice point num
function nextchoicepointnum(gencontext::GeneratorContext)
	gencontext.numchoicepoints += 1
	gencontext.choicepointoffset + gencontext.numchoicepoints
end

# generate a unique rule name
uniquemethodname(rulename) = gensym(string(rulename))

# if node is a just rule name, convert it to a call to rule with no parameters
function convertrulenametorulecall(node, genrules::Vector{GeneratorRule})
	if isa(node, Symbol) && (node in map(genrule->genrule.rulename, genrules))
		return Expr(:call, node)
	end
	node
end

# generic method for parsing node tree to identify and then transform particular types of function calls
# matchfn should return true if the desired type of function call is identified
# transformfn should return the tree to replace the node containing the call
# note: transformfn must handle recursively checking any nodes in its argument list itself (which is why matchfn is passed)
function transformfunccall(node, genrules::Vector{GeneratorRule}, gencontext::GeneratorContext, matchfn::Function, transformfn::Function)

	if islinenode(node)
		return node
	end

	# if node is a reference to an explicit DataGenerator method, then don't transform and go no deeper
	if (typeof(node) == Expr) && (node.head == Symbol(".")) && (node.args[1] == THIS_MODULE)
		return node
	end

	callname, callparams = extractfunccall(node)

	if (callname != nothing) && matchfn(callname, callparams, gencontext)
		return transformfn(callname, callparams, genrules, gencontext, matchfn, transformfn)
	end

	if typeof(node) == Expr
		node.args = map(arg -> transformfunccall(arg, genrules, gencontext, matchfn, transformfn), node.args)
		return node
	end

	return node

end

# process a param that should be a specified literal datatype or an expression
# note: cannot simply check for Expr or not since call to GT-specific no-arg func would be non-Expr but is non-literal
function processpossiblyliteralparam(param, datatype, gencontext::GeneratorContext)
	value = param
	paramisliteral = true
	try
		value = convert(datatype, param) # see if value can be interpreted as a literal of specified datatype
	catch
		paramisliteral = false # assume instead is not a literal
	end
	value, paramisliteral
end

# extract an parameter name from the function argument expression
# the argument *may*
#		- have type, e.g. x::Int
#		- have default value e.g. x=2
#		- be a named parameter
# TODO vargs (i.e. ...)
function extractparamfromarg(argexpr)
	if typeof(argexpr) == Symbol
		return argexpr
	end
	if typeof(argexpr) == Expr
		filteredargs = removelinenodes(argexpr.args)
		if (argexpr.head == :(::)) && (length(filteredargs) == 2)
			return extractparamfromarg(filteredargs[1])
		end
	end
	if (argexpr.head == :kw) && (length(filteredargs) == 2) # :kw indicates parameter defaulting
		return extractparamfromarg(filteredargs[1])
	end
	if argexpr.head == :parameter # :parameter indicates parameter defaulting
		# TODO - difficulty here is that these named parameters come first, and when using them need to use a special syntax, so disallow for the moment
		error("currently cannot use named parameters in a rule choice rule")
	end
	error("cannot extract parameter name from argument definition")
end

# attempt to convert to DataType (otherwise returns nothing)
function converttodatatype(s)
	datatype = nothing
	if (typeof(s) == Symbol) || ((typeof(s) == Expr) && (s.head == :curly))
		# :curly expression for Union{...}, Vector{...} etc.
		try
			datatype = eval(s)
			# TODO could eval perform side-effects here? should be OK when we know it is a symbol, but what about curly expression?
		catch
			datatype = nothing
		end
		if !(typeof(datatype) <: Type) # Type rather than datatype as Union{} is a subtype of Type, but not a DataType
			datatype = nothing
		end
	end
	datatype
end

# test whether syntax tree node in an expression is a line node (i.e. a file and line location)
# (as of v0.3, LineNumberNode does not seem to occur)
islinenode(node) = ((typeof(node) == Expr) && (node.head == :line)) || (typeof(node) == LineNumberNode)

# remove nodes in array that are line nodes
removelinenodes(args) = filter(arg->!islinenode(arg),args)

# extract method signature into tuple (name, arguments)
# (1) standard long or short form: foo(bar) or foo(bar) =
# (2) GT-specific no args form: foo (but only when allownoparenthesisform is true)
function extractfuncsig(node, allownoparenthesisform = false)
	if typeof(node) == Expr
		if node.head == :call # standard long or short form
			callargs = removelinenodes(node.args)
			if length(callargs) >= 1
				return (callargs[1], callargs[2:end])
			end
		end
	elseif allownoparenthesisform && (typeof(node) == Symbol) # GT-specific no parenthesis form when no args
		return (node, [])
	end
	return (nothing,nothing)
end

# extract function definition into tuple (name, arguments, body)
# (1) long form: function foo(bar) ... end
# (2) short form: foo(bar) = ...
# (3) GT-specific no args form: foo = ... (but only when allownoparenthesisform is true)
function extractfuncdef(node, allownoparenthesisform = false)
	if typeof(node) == Expr
		if (node.head == :function) || (node.head == :(=) )
			eqargs = removelinenodes(node.args)
			if length(eqargs) == 2
				(name, args) = extractfuncsig(eqargs[1], allownoparenthesisform)
				if name != nothing
					return (name, args, eqargs[2])
				end
			end
		end
	end
	return (nothing,nothing,nothing)
end

# parse function call into tuple (name, params)
# (1) standard form: foo(bar) = ...
# (2) GT-specific no args form: foo (but only when allownoparenthesisform is true)
function extractfunccall(node, allownoparenthesisform = false)
	if typeof(node) == Expr
		if node.head == :call # standard long or short form
			callargs = removelinenodes(node.args)
			if length(callargs) >= 1
				return (callargs[1], callargs[2:end])
			end
		end
	elseif allownoparenthesisform && (typeof(node) == Symbol)
		return(node, [])
	end
	return (nothing,nothing)
end

# merge two expression into a block
# if both are already blocks, then combine contents as a single block
function mergeexprs(ex1::Expr, ex2::Expr)
	if ex1.head == :block && ex2.head == :block
		Expr(:block, ex1.args..., ex2.args...)
	else
		Expr(:block, ex1, ex2)
	end
end

# parses rule sources (as returned by translator) and appends them generator rules
function addrulesources(genrules::Vector{GeneratorRule}, rulesources::Vector{DataGeneratorTranslators.RuleSource})
	for rulesource in rulesources
		genrule = GeneratorRule(rulesource.rulename, rulesource.args, parse(join(rulesource.source, "\n")))
		push!(genrules, genrule)
	end
end


text:
const QuasiArrayMulArray{p, q, T, V} = Mul{<:Any, <:Any, <:AbstractQuasiArray{T,p}, <:AbstractArray{V,q}}
const ArrayMulQuasiArray{p, q, T, V} = Mul{<:Any, <:Any, <:AbstractArray{T,p}, <:AbstractQuasiArray{V,q}}
const QuasiArrayMulQuasiArray{p, q, T, V} = Mul{<:Any, <:Any, <:AbstractQuasiArray{T,p}, <:AbstractQuasiArray{V,q}}
####
# Matrix * Vector
####
const QuasiMatMulVec{T, V} = QuasiArrayMulArray{2, 1, T, V}
const QuasiMatMulMat{T, V} = QuasiArrayMulArray{2, 2, T, V}
const QuasiMatMulQuasiMat{T, V} = QuasiArrayMulQuasiArray{2, 2, T, V}

combine_mul_styles(::MulStyle, ::QuasiArrayApplyStyle) = QuasiArrayApplyStyle()
combine_mul_styles(::QuasiArrayApplyStyle, ::MulStyle) = QuasiArrayApplyStyle()
combine_mul_styles(::DefaultArrayApplyStyle, ::QuasiArrayApplyStyle) = QuasiArrayApplyStyle()
combine_mul_styles(::QuasiArrayApplyStyle, ::DefaultArrayApplyStyle) = QuasiArrayApplyStyle()

ApplyStyle(::typeof(*), ::Type{<:AbstractQuasiArray}) = QuasiArrayApplyStyle()
ApplyStyle(::typeof(*), ::Type{<:AbstractArray}, ::Type{<:AbstractQuasiArray}) = MulStyle()
ApplyStyle(::typeof(*), ::Type{<:AbstractQuasiArray}, ::Type{<:AbstractArray}) = MulStyle()
ApplyStyle(::typeof(*), ::Type{<:AbstractQuasiArray}, ::Type{<:AbstractQuasiArray}) = MulStyle()

*(A::AbstractQuasiMatrix) = copy(A)
*(A::AbstractQuasiArray, B::AbstractQuasiArray) = mul(A, B)
*(A::AbstractArray, B::AbstractQuasiArray) = mul(A, B)
*(A::AbstractQuasiArray, B::AbstractArray) = mul(A, B)

for op in (:pinv, :inv)
    @eval $op(A::AbstractQuasiArray) = apply($op,A)
end

@inline axes(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractQuasiMatrix}) = (axes(L.A, 2),axes(L.B,2))
@inline axes(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractQuasiVector}) = (axes(L.A, 2),)

@inline \(A::AbstractQuasiArray, B::AbstractQuasiArray) = ldiv(A,B)
@inline \(A::AbstractQuasiArray, B::AbstractArray) = ldiv(A,B)
@inline \(A::AbstractArray, B::AbstractQuasiArray) = ldiv(A,B)

@inline /(A::AbstractQuasiArray, B::AbstractQuasiArray) = rdiv(A,B)
@inline /(A::AbstractQuasiArray, B::AbstractArray) = rdiv(A,B)
@inline /(A::AbstractArray, B::AbstractQuasiArray) = rdiv(A,B)

####
# MulQuasiArray
#####

const MulQuasiArray{T, N, Args<:Tuple} = ApplyQuasiArray{T, N, typeof(*), Args}

const MulQuasiVector{T, Args<:Tuple} = MulQuasiArray{T, 1, Args}
const MulQuasiMatrix{T, Args<:Tuple} = MulQuasiArray{T, 2, Args}

const Vec = MulQuasiVector

_ApplyArray(F, factors...) = ApplyQuasiArray(F, factors...)
_ApplyArray(F, factors::AbstractArray...) = ApplyArray(F, factors...)

MulQuasiOrArray = Union{MulArray,MulQuasiArray}

_factors(M::MulQuasiOrArray) = M.args
_factors(M) = (M,)

@inline _flatten(A::MulQuasiArray, B...) = _flatten(Applied(A), B...)
@inline flatten(A::MulQuasiArray) = ApplyQuasiArray(flatten(Applied(A)))
@inline flatten(A::SubQuasiArray{<:Any,2,<:MulQuasiArray}) = materialize(flatten(Applied(A)))

adjoint(A::MulQuasiArray) = ApplyQuasiArray(*, reverse(adjoint.(A.args))...)
transpose(A::MulQuasiArray) = ApplyQuasiArray(*, reverse(transpose.(A.args))...)

function similar(A::MulQuasiArray)
    B,a = A.args
    B*similar(a)
end

function similar(A::QuasiArrayMulArray)
    B,a = A.args
    applied(*, B, similar(a))
end

function copyto!(dest::MulQuasiArray, src::MulQuasiArray)
    d = last(dest.args)
    s = last(src.args)
    copyto!(IndexStyle(d), d, IndexStyle(s), s)
    dest
end

Array(M::Mul{<:Any,<:Any,<:Any,<:AbstractQuasiMatrix}) = eltype(M)[M[k,j] for k in axes(M)[1], j in axes(M)[2]]
Array(M::Mul{<:Any,<:Any,<:Any,<:AbstractQuasiVector}) = eltype(M)[M[k] for k in axes(M)[1]]
_quasi_mul(M, _) = QuasiArray(M)
_quasi_mul(M, ::NTuple{N,OneTo{Int}}) where N = Array(M)
copy(M::Mul{QuasiArrayLayout,QuasiArrayLayout}) = _quasi_mul(M, axes(M))
copy(M::Mul{QuasiArrayLayout}) = _quasi_mul(M, axes(M))
copy(M::Mul{<:Any,QuasiArrayLayout}) = _quasi_mul(M, axes(M))
copy(M::Mul{<:AbstractLazyLayout,QuasiArrayLayout}) = ApplyQuasiArray(M)
copy(M::Mul{ApplyLayout{typeof(\)},QuasiArrayLayout}) = ApplyQuasiArray(M)
copy(M::Mul{QuasiArrayLayout,<:AbstractLazyLayout}) = ApplyQuasiArray(M)
for op in (:+, :-)
    @eval begin
        copy(M::Mul{BroadcastLayout{typeof($op)},QuasiArrayLayout}) = simplify(M)
        copy(M::Mul{QuasiArrayLayout,BroadcastLayout{typeof($op)}}) = simplify(M)
    end
end
@inline copy(M::Mul{BroadcastLayout{typeof(*)},QuasiArrayLayout}) = copy(Mul{BroadcastLayout{typeof(*)},UnknownLayout}(M.A,M.B))

LazyArrays._vec_mul_view(a::AbstractQuasiVector, kr, ::Colon) = view(a, kr)

###
# Scalar special case, simplifies x * A and A * x
# Find an AbstractArray to latch on to by commuting
###

_lmul_scal_reduce(x::Number, A) = (x * A,)
_lmul_scal_reduce(x::Number, A::AbstractArray, B...) = (x * A, B...)
_lmul_scal_reduce(x::Number, A, B...) = (A, _lmul_scal_reduce(x, B...)...)

_rmul_scal_reduce(x::Number, Z) = (Z * x,)
_rmul_scal_reduce(x::Number, Z::AbstractArray, Y...) = (Y..., Z*x)
_rmul_scal_reduce(x::Number, Z, Y...) = (_rmul_scal_reduce(x, Y...)..., Z)

_ldiv_scal_reduce(x::Number, A) = (x \ A,)
_ldiv_scal_reduce(x::Number, A::AbstractArray, B...) = (x \ A, B...)
_ldiv_scal_reduce(x::Number, A, B...) = (A, _ldiv_scal_reduce(x, B...)...)

_rdiv_scal_reduce(x::Number, Z) = (Z / x,)
_rdiv_scal_reduce(x::Number, Z::AbstractArray, Y...) = (Y..., Z/x)
_rdiv_scal_reduce(x::Number, Z, Y...) = (_rdiv_scal_reduce(x, Y...)..., Z)

broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(*), x::Number, A::MulQuasiArray{<:Any,N}) where N =
    ApplyQuasiArray(*, _lmul_scal_reduce(x, arguments(A)...)...)
broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(*), A::MulQuasiArray{<:Any,N}, x::Number) where N =
    ApplyQuasiArray(*, _rmul_scal_reduce(x, reverse(arguments(A))...)...)

broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(\), x::Number, A::MulQuasiArray{<:Any,N}) where N =
    ApplyQuasiArray(*, _ldiv_scal_reduce(x, arguments(A)...)...)
broadcasted(::LazyQuasiArrayStyle{N}, ::typeof(/), A::MulQuasiArray{<:Any,N}, x::Number) where N =
    ApplyQuasiArray(*, _ldiv_scal_reduce(x, arguments(A)...)...)


text:
##############################################################################
##
## Iterate on terms
##
##############################################################################
eachterm(x::AbstractTerm) = (x,)
eachterm(x::NTuple{N, AbstractTerm}) where {N} = x

##############################################################################
##
## Parse IV
##
##############################################################################

function decompose_iv(f::FormulaTerm)
	formula_endo = nothing
	formula_iv = nothing
	for term in eachterm(f.rhs)
		if isa(term, FormulaTerm)
			if formula_endo != nothing
				throw("There can only be one instrumental variable specification")
			end
			formula_endo = FormulaTerm(ConstantTerm(0), tuple(ConstantTerm(0), eachterm(term.lhs)...))
			formula_iv = FormulaTerm(ConstantTerm(0), tuple(ConstantTerm(0), eachterm(term.rhs)...))
		end
	end
	return FormulaTerm(f.lhs, tuple((term for term in eachterm(f.rhs) if !isa(term, FormulaTerm))...)), formula_endo, formula_iv
end

##############################################################################
##
## Parse FixedEffect
##
##############################################################################
struct FixedEffectTerm <: AbstractTerm
    x::Symbol
end
StatsModels.termvars(t::FixedEffectTerm) = [t.x]
fe(x::Term) = FixedEffectTerm(Symbol(x))

has_fe(::FixedEffectTerm) = true
has_fe(::FunctionTerm{typeof(fe)}) = true
has_fe(t::InteractionTerm) = any(has_fe(x) for x in t.terms)
has_fe(::AbstractTerm) = false
has_fe(t::FormulaTerm) = any(has_fe(x) for x in eachterm(t.rhs))

fesymbol(t::FixedEffectTerm) = t.x
fesymbol(t::FunctionTerm{typeof(fe)}) = Symbol(t.args_parsed[1])

function parse_fixedeffect(df::AbstractDataFrame, formula::FormulaTerm)
    fe = FixedEffect[]
    id = Symbol[]
    for term in eachterm(formula.rhs)
        result = parse_fixedeffect(df, term)
        if result != nothing
            push!(fe, result[1])
            push!(id, result[2])
        end
    end
    formula = FormulaTerm(formula.lhs, tuple((term for term in eachterm(formula.rhs) if !has_fe(term))...))
    return fe, id, formula
end

# Constructors from dataframe + Term
function parse_fixedeffect(df::AbstractDataFrame, t::AbstractTerm)
    if has_fe(t)
        st = fesymbol(t)
        return FixedEffect(df[!, st]), Symbol(:fe_, st)
    end
end

# Constructors from dataframe + InteractionTerm
function parse_fixedeffect(df::AbstractDataFrame, t::InteractionTerm)
    fes = (x for x in t.terms if has_fe(x))
    interactions = (x for x in t.terms if !has_fe(x))
    if !isempty(fes)
        # x1&x2 from (x1&x2)*id
        fe_names = [fesymbol(x) for x in fes]
        fe = FixedEffect(group((df[!, fe_name] for fe_name in fe_names)...); interaction = _multiply(df, Symbol.(interactions)))
        interactions = setdiff(Symbol.(terms(t)), fe_names)
        s = vcat(["fe_" * string(fe_name) for fe_name in fe_names], string.(interactions))
        return fe, Symbol(reduce((x1, x2) -> x1*"&"*x2, s))
    end
end

function _multiply(df, ss::Vector)
    if isempty(ss)
        out = Ones(size(df, 1))
    else
        out = ones(size(df, 1))
        for j in eachindex(ss)
            _multiply!(out, df[!, ss[j]])
        end
    end
    return out
end
function _multiply!(out, v)
    if v isa CategoricalVector
        throw("Fixed Effects cannot be interacted with Categorical Vector. Use fe(x)&fe(y)")
    end
    for i in eachindex(out)
        if v[i] === missing
            # may be missing when I remove singletons
            out[i] = 0.0
        else
            out[i] = out[i] * v[i]
        end
    end
end


text:
function ueberpruefeLGS(Gleichungssystem,f1,f2,f3,f4)

if Gleichungssystem(f1,f2,f3,f4) == [0;0;0;0]
    println("Ihr habt das lineare Gleichungssystem korrekt aufgestellt!")
else
    println("Überprüft noch einmal eure Eingabe.")
end

end


text:
using Plots, LaTeXStrings, Statistics,  Distributions,StatsPlots, JLD
TotalNum=10000
NList=[5,10,100,1000]
Data=[[] for i in 1:length(NList)]
for N in 1:length(NList)
    for i in 1:TotalNum
        push!(Data[N], sum(rand(NList[N])))
    end
end
plots=[]
colors=[:silver, :lightblue,:pink2,:mediumpurple3]
for i in 1:length(NList)
    plt=begin
        histogram(Data[i],color=colors[i],normalize = true, label=L"Histogram\ of\ Samples", alpha=0.9)
        plot!(Normal(mean(Data[i]),std(Data[i])), c = :black, label= L"Normal(\mu = %$(round(mean(Data[i]),digits = 2)), \sigma = %$(round(std(Data[i]),digits = 2)))")
    end
    push!(plots, plt)
end
plot(plots..., size = (1000,1000), plot_title = "Normalized Distribution of Samples", dpi=600)
save("../../chapter6/6.3/samples.jld","Samples", Data)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter6\\6.3\\Figs\\N_Distribution.png")


text:
"""
Generalized Gaussian Distribution package.
"""
module GGD

using Distributions, SpecialFunctions, Statistics, Random
using SpecialFunctions: gamma
using Distributions: Bernoulli, Gamma, sampler, rand
import Base: rand
import Random: GLOBAL_RNG

### source files
include("generalizedgaussian.jl")

### methods
export
   GeneralizedGaussian,     # principle distribution type.
   params,                  # returns shape, location, and scale parameters.
   shape,                   # shape parameter β.
   scale,                   # scale parameter α, which is related to variance via the shape parameter and gamma functions.
   location,                # location parameter μ, the same as the mean/median/mode.
   std,                     # standard deviation σ.
   mean,                    # mean of the distribution (equal to μ)
   median,                  # median of distribution
   mode,                    # mode of distribution
   var,                     # variance of distribution, (σ²)
   skewness,                # skewness of Generalized Gaussian Distribution is zero.
   kurtosis,                # Excess kurtosis of distribution, (this subtracts 3.)
   entropy,                 # entropy of distribution
   logpdf,                  # log probability density
   cdf,                     # cummulative distribution function
   pdf,                     # probability density function
   rand                     # random sample from the distribution

"""
A Julia Package for the Generalized Gaussian Distribution (GGD).

API:
- `d = GGD(params...)` creates the GGD Type, containing the relevant
    distribution parameters, including mu, sigma, alpha, and beta.
- `z = rand(rng, d)` will generate a sample from the specified GGD
    distribution "d".
- `pdf(d, x)` computes the probability density function at x for GGD
    distribution "d".
- `cdf(d, x)` computes the cumulative distribution function at x.

"""

end # module


text:
JSON.lower(a::Aggregator) = non_nothing_dict(a)

"""
    Count(output_name)

Count number of rows.
"""
struct Count <: Aggregator
    type::String
    name::String
    Count(name) = new("count", name)
end

"""
    SingleField(type, field, output_name)

Apply `type` aggregator to the `field` and get result as `output_name`.

Available `type`s: "longSum", "doubleSum", "floatSum", "longMin", "doubleMin",
"floatMin", "longMax", "doubleMax", "floatMax", "longFirst", "doubleFirst",
"floatFirst", "longLast", "doubleLast", "floatLast", "longAny", "doubleAny",
"floatAny", "doubleMean"
"""
struct SingleField <: Aggregator
    type::String
    fieldName::String
    name::String
    function SingleField(type, of, as)
        type ∈ [
            "longSum", "doubleSum", "floatSum",
            "longMin", "doubleMin", "floatMin",
            "longMax", "doubleMax", "floatMax",
            "longFirst", "doubleFirst", "floatFirst",
            "longLast", "doubleLast", "floatLast",
            "longAny", "doubleAny", "floatAny",
            "doubleMean"
        ] || error("Unknown SingleField aggregation")
        new(type, of, as)
    end
end

"""
    StringAgg(type, field, output_name; maxStringBytes=nothing)

Apply `type` string aggregator to `field` and get result as `output_name`.

Available `type`s: "first", "last", "any"

See Druid documentation for `maxStringBytes`.
"""
struct StringAgg <: Aggregator
    type::String
    fieldName::String
    name::String
    maxStringBytes
    function StringAgg(type, of, as; maxStringBytes=nothing)
        type = lowercase(type)
        type ∈ ["first", "last", "any"] || error("Unknown StringAgg aggregation")
        maxStringBytes === nothing || (isa(maxStringBytes, Integer) && maxStringBytes >= 0) || error("maxStringBytes must be a positive integer")
        new("string" * titlecase(type), of, as, maxStringBytes)
    end
end

"""
    Grouping(groupings, name)

Grouping aggregator
"""
struct Grouping <: Aggregator
    type::String
    groupings::Vector{String}
    name::String
    Grouping(groupings, name) = new("grouping", groupings, name)
end

"""
    Filtered(filter, aggregator)

An aggregator that combines a filter with the `aggregator`.
"""
struct Filtered <: Aggregator
    type::String
    filter::Filter
    aggregator::Aggregator
    Filtered(filter, aggregator) = new("filtered", filter, aggregator)
end


text:
"""
    function store_data(metrics::BasicMetrics{O, H}, title::String) where{O<:AbstractTakeObjective, H<:ValueSelection}

Store useful results from consecutive search in `.csv` file.
"""
function storedata(metrics::AbstractMetrics; filename::String="")
    df = DataFrame(
        Episode = Int[],
        Instance = Missing[],
        Solution = Int[],
        Nodes = Int[],
        Time = Float64[],
        Score = Union{Missing, Float64}[],
        Reward = Union{Missing, Float64}[],
        Loss = Union{Missing, Float64}[]
    )
    for i = 1:metrics.nbEpisodes
        episodeData = Dict(
            :Episode => i,
            :Instance => missing,
            :Solution => 0,
            :Nodes => metrics.meanNodeVisitedUntilEnd[i],
            :Time => metrics.timeneeded[i],
            :Score => missing,
            :Reward => isnothing(metrics.totalReward) ? missing : metrics.totalReward[i],
            :Loss => isnothing(metrics.loss) ? missing : metrics.loss[i]
        )
        push!(df, episodeData)
        for j = 1:length(metrics.nodeVisited[i])
            solutionData = copy(episodeData)
            solutionData[:Solution] = j
            solutionData[:Nodes] = metrics.nodeVisited[i][j]
            solutionData[:Score] = isnothing(metrics.scores) | isnothing(metrics.scores[i][j]) ? missing : metrics.scores[i][j]
            push!(df, solutionData)
        end
    end
    CSV.write(filename*".csv", df)
    return df
end

function storedata(metrics::Vector{<:AbstractMetrics}; filename::String="")
    df = DataFrame(
        Episode = Int[],
        Instance = Int[],
        Solution = Int[],
        Nodes = Int[],
        Time = Float64[],
        Score = Union{Missing, Float64}[],
        Reward = Union{Missing, Float64}[],
        Loss = Union{Missing, Float64}[]
    )
    nbInstances = length(metrics)
    for j in 1:nbInstances, i = 1:metrics[j].nbEpisodes
        episodeData = Dict(
            :Episode => i,
            :Instance => j,
            :Solution => 0,
            :Nodes => metrics[j].meanNodeVisitedUntilEnd[i],
            :Time => metrics[j].timeneeded[i],
            :Score => missing,
            :Reward => isnothing(metrics[j].totalReward) ? missing : metrics[j].totalReward[i],
            :Loss => isnothing(metrics[j].loss) ? missing : metrics[j].loss[i]
        )
        push!(df, episodeData)
        for k = 1:length(metrics[j].nodeVisited[i])
            solutionData = copy(episodeData)
            solutionData[:Solution] = k
            solutionData[:Nodes] = metrics[j].nodeVisited[i][k]
            solutionData[:Score] = isnothing(metrics[j].scores) | isnothing(metrics[j].scores[i][k]) ? missing : metrics[j].scores[i][k]
            push!(df, solutionData)
        end
    end
    CSV.write(filename*".csv", df)
    return df
end


text:
function acf(Z1)
    fZ1 = fft(Z1)
    pZ1 = abs2(fZ1)
    fftshift(real(ifft(pZ1)))
end


text:
using Plots, StatsPlots
using Plots.Measures
using DataFrames, Dates, Statistics, Distributions, LaTeXStrings
import Bootstrap: bootstrap, BalancedSampling

import .Analysis: TemporalCrossSection

function quickplot(node::NetworkNode)
    fig = plot(node.outflow)
    return fig
end

function quickplot(node::NetworkNode, climate::Climate)
    date = timesteps(climate)

    @assert length(date) == length(node.outflow) || "Date length and result lengths do not match!"

    fig = plot(date, node.outflow)

    return fig
end

function quickplot(obs, node::NetworkNode, climate::Climate, label="", log=true; burn_in=1, limit=nothing, metric=Streamfall.mKGE)
    return quickplot(obs, node.outflow, climate, label, log; burn_in=burn_in, limit=limit, metric=metric)
end

function quickplot(obs::Array, sim::Array, climate::Climate, label="", log=true; burn_in=1, limit=nothing, metric=Streamfall.mKGE)
    date = timesteps(climate)
    last_e = !isnothing(limit) ? limit : lastindex(obs)
    show_range = burn_in:last_e
    return quickplot(obs[show_range], sim[show_range], date[show_range], label, log; metric=metric)
end

function quickplot(obs::Array, sim::Array, xticklabels::Array, label="Modeled", log=true; metric=Streamfall.mKGE)
    @assert length(xticklabels) == length(obs) || "x-axis tick label length and observed lengths do not match!"
    @assert length(xticklabels) == length(sim) || "x-axis tick label length and simulated lengths do not match!"

    score = round(metric(obs, sim), digits=4)
    metric_name = String(Symbol(metric))

    if log
        # Add small constant in case of 0-flow
        obs = obs .+ 1e-2
        sim = sim .+ 1e-2
    end

    label = "$(label) ($(metric_name): $(score))"
    fig = plot(xticklabels, obs,
                label="Observed",
                legend=:best,
                ylabel="Streamflow",
                xlabel="Date",
                fg_legend=:transparent,
                bg_legend=:transparent)
    plot!(xticklabels, sim, label=label, alpha=0.7)

    if log
        # modify yaxis
        yaxis!(fig, :log10)
    end

    qqfig = qqplot(obs, sim, legend=false, markerstrokewidth=0, alpha=0.7, xlabel="Observed", ylabel="Modeled")

    if log
        xaxis!(qqfig, :log10)
        yaxis!(qqfig, :log10)
    end

    combined = plot(fig, qqfig, size=(800, 400), left_margin=10mm, layout=(1,2))

    return combined
end

"""
    plot_residuals(obs::Array, sim::Array; xlabel="", ylabel="", title="")

Plot residual between two sequences.

# Arguments
- x : x-axis data
- y : y-axis data
- xlabel : x-axis label
- ylabel : y-axis label
- title : title text
"""
function plot_residuals(x::Array, y::Array; xlabel="", ylabel="", title="")
    # 1:1 Plot
    fig_1to1 = scatter(x, y, legend=false,
                       markerstrokewidth=0, markerstrokealpha=0, alpha=0.2)
    plot!(x, y, color=:red, markersize=0.1, markerstrokewidth=0,
    xlabel=xlabel, ylabel=ylabel, title=title)

    return fig_1to1
end

"""Symmetrical log values.

https://kar.kent.ac.uk/32810/2/2012_Bi-symmetric-log-transformation_v5.pdf
https://discourse.julialang.org/t/symmetrical-log-plot/45709/3
"""
function symlog(y)
    return sign.(y) .* log10.(1.0 .+ abs.(y))
end

"""
    temporal_cross_section(dates, obs; ylabel=nothing, period::Function=month)

Provides indication of temporal variation and uncertainty across time, grouped by `period`.

Notes:
Assumes daily data.
Filters out leap days.

# Arguments
- dates : Date of each observation
- obs : observed data
- ylabel : Optional replacement ylabel. Uses name of `func` if not provided.
- `period::Function` : Method from `Dates` package to group (defaults to `monthday`)
"""
function temporal_cross_section(dates, obs;
                                title="", ylabel="Median Error", label=nothing,
                                period::Function=monthday,
                                kwargs...)  # show_extremes::Bool=false,
    if isnothing(label)
        label = ylabel
    end

    arg_keys = keys(kwargs)
    format_func = y -> y
    logscale = [:log, :log10]
    tmp = nothing

    xsect_res = TemporalCrossSection(dates, obs, period)

    if :yscale in arg_keys || :yaxis in arg_keys
        tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis]

        if tmp in logscale
            log_obs = symlog(copy(obs))

            # Format function for y-axis tick labels (e.g., 10^x)
            format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0"

            log_xsect_res = TemporalCrossSection(dates, log_obs, period)
            target = log_xsect_res.cross_section
        else
            target = xsect_res.cross_section
        end
    end

    x_section = target.median
    lower_75 = target.lower_75
    upper_75 = target.upper_75
    lower_95 = target.lower_95
    upper_95 = target.upper_95

    sp = target.subperiod
    xlabels = join.(sp, "-")

    if !isnothing(tmp) & (tmp in logscale)
        # Remove keys so later plotting methods do not error
        kwargs = Dict(kwargs)
        delete!(kwargs, :yscale)
        delete!(kwargs, :yaxis)
    end

    xsect_df = xsect_res.cross_section

    # Display indicator values using original data instead of log-transformed data
    med = xsect_df.median
    m_ind = round(mean(med), digits=2)
    sd_ind = round(std(med), digits=2)

    wr75_m_ind = round(xsect_res.mean_75, digits=2)
    wr75_sd_ind = round(xsect_res.std_75, digits=2)
    wr95_m_ind = round(xsect_res.mean_95, digits=2)
    wr95_sd_ind = round(xsect_res.std_95, digits=2)

    fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0)
    plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0)
    plot!(fig, xlabels, x_section,
            label="$(label) μ: $(m_ind), σ: $(sd_ind)",
            color="black",
            xlabel=nameof(period),
            ylabel=ylabel,
            legend=:bottomleft,
            legendfont=Plots.font(10),
            fg_legend=:transparent,
            bg_legend=:transparent,
            left_margin=5mm,
            bottom_margin=5mm,
            title=title,
            yformatter=format_func;
            kwargs...)

    # if show_extremes
    #     scatter!(fig, xlabels, min_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...)
    #     scatter!(fig, xlabels, max_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...)
    # end

    return fig
end

"""
    temporal_cross_section(dates, obs, sim; ylabel=nothing, period::Function=month)

Provides indication of temporal variation and uncertainty across time, grouped by `period`.

Notes:
Assumes daily data.
Filters out leap days.

# Arguments
- dates : Date of each observation
- obs : observed data
- ylabel : Optional replacement ylabel. Uses name of `func` if not provided.
- `period::Function` : Method from `Dates` package to group (defaults to `month`)
"""
function temporal_cross_section(dates, obs, sim;
                                title="", ylabel="Median Error", label=nothing,
                                period::Function=monthday,
                                kwargs...)  # show_extremes::Bool=false,
    if isnothing(label)
        label = ylabel
    end

    arg_keys = keys(kwargs)
    format_func = y -> y
    logscale = [:log, :log10]
    tmp = nothing

    xsect_res = TemporalCrossSection(dates, obs, sim)
    target = xsect_res.cross_section

    if :yscale in arg_keys || :yaxis in arg_keys
        tmp = (:yscale in arg_keys) ? kwargs[:yscale] : kwargs[:yaxis]

        if tmp in logscale
            log_obs = symlog(copy(obs))
            log_sim = symlog(copy(sim))

            # Format function for y-axis tick labels (e.g., 10^x)
            format_func = y -> (y != 0) ? L"%$(Int(round(sign(y)) * 10))^{%$(round(abs(y), digits=1))}" : L"0"

            log_xsect_res = TemporalCrossSection(dates, log_obs, log_sim)
            target = log_xsect_res.cross_section
        end
    end

    x_section = target.median
    lower_75 = target.lower_75
    upper_75 = target.upper_75
    lower_95 = target.lower_95
    upper_95 = target.upper_95

    sp = target.subperiod
    xlabels = join.(sp, "-")

    if !isnothing(tmp) & (tmp in logscale)
        # Remove keys so later plotting methods do not error
        kwargs = Dict(kwargs)
        delete!(kwargs, :yscale)
        delete!(kwargs, :yaxis)
    end

    xsect_df = xsect_res.cross_section

    # Display indicator values using original data instead of log-transformed data
    med = xsect_df.median
    m_ind = round(mean(med), digits=2)
    sd_ind = round(std(med), digits=2)

    wr75_m_ind = round(xsect_res.mean_75, digits=2)
    wr75_sd_ind = round(xsect_res.std_75, digits=2)
    wr95_m_ind = round(xsect_res.mean_95, digits=2)
    wr95_sd_ind = round(xsect_res.std_95, digits=2)

    fig = plot(xlabels, lower_95, fillrange=upper_95, color="lightblue", alpha=0.3, label="CI₉₅ μ: $(wr95_m_ind), σ: $(wr95_sd_ind)", linealpha=0)
    plot!(fig, xlabels, lower_75, fillrange=upper_75, color="lightblue", alpha=0.5, label="CI₇₅ μ: $(wr75_m_ind), σ: $(wr75_sd_ind)", linealpha=0)
    plot!(fig, xlabels, x_section,
            label="$(label) μ: $(m_ind), σ: $(sd_ind)",
            color="black",
            xlabel=nameof(period),
            ylabel=ylabel,
            legend=:bottomleft,
            legendfont=Plots.font(10),
            fg_legend=:transparent,
            bg_legend=:transparent,
            left_margin=5mm,
            bottom_margin=5mm,
            title=title,
            yformatter=format_func;
            kwargs...)

    # if show_extremes
    #     scatter!(fig, xlabels, min_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...)
    #     scatter!(fig, xlabels, max_section, label="", alpha=0.5, color="lightblue", markerstrokewidth=0; kwargs...)
    # end

    return fig
end


text:
module TestDataFrame
    using Base.Test
    using DataFrames

    #test_group("Operations on DataFrames that have column groupings")

    x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])
    y = DataFrame(c = [1, 2, 3], d = [4, 5, 6])

    #
    # Equality
    #
    @test isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2]), b=@data([4, 5])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), c=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(b=@data([4, 5, 6]), a=@data([1, 2, 3])))
    @test !isequal(DataFrame(a=@data([1, 2, 2]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test isequal(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])))

    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2]), b=@data([4, 5])) != DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]), c=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(b=@data([4, 5, 6]), a=@data([1, 2, 3]))
    @test DataFrame(a=@data([1, 2, 2]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 3, NA]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6]))
    @test isna(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])))
    @test isna(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))

    z = deepcopy(x)

    z = hcat(x, y)

    v = DataFrame(
        a = [5,6,7],
        b = [8,9,10]
    )
    z = vcat(DataFrame[v, x])

    z = vcat(v, x)

    # Deleting columns removes any mention from groupings
    delete!(x, :a)
    @test names(x) == [:b]

    ## del calls ref, which properly deals with groupings
    z2 = z[:,[1,1,2]]
    @test names(z2) == [:a, :a_1, :b]

    #test_group("DataFrame assignment")
    # Insert single column
    x[symbol("c")] = 1:3
    @test_throws ErrorException x[symbol("\u212b")] = 1:3
    @test_throws ErrorException x[symbol("end")] = 1:3
    @test_throws ErrorException x[symbol("1a")] = 1:3

    x0 = x[[], :]
    @test_throws ErrorException x0[:d] = [1]
    @test_throws ErrorException x0[:d] = 1:3

    # Insert single value
    x[:d] = 3
    @test x[:d] == [3, 3, 3]

    x0[:d] = 3
    @test x0[:d] == Int[]

    #test_group("Empty DataFrame constructors")
    df = DataFrame(Int, 10, 3)
    @test size(df, 1) == 10
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Int}
    @test typeof(df[:, 3]) == DataVector{Int}

    df = DataFrame({Int, Float64, ASCIIString}, 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Float64}
    @test typeof(df[:, 3]) == DataVector{ASCIIString}

    df = DataFrame({Int, Float64, ASCIIString}, [:A, :B, :C], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Float64}
    @test typeof(df[:, 3]) == DataVector{ASCIIString}

    df = convert(DataFrame, zeros(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, ones(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, eye(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    #test_group("Other DataFrame constructors")
    df = DataFrame([{:a=>1, :b=>'c'}, {:a=>3, :b=>'d'}, {:a=>5}])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == DataVector{Int}
    @test typeof(df[:,:b]) == DataVector{Char}

    df = DataFrame([{:a=>1, :b=>'c'}, {:a=>3, :b=>'d'}, {:a=>5}], [:a, :b])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == DataVector{Int}
    @test typeof(df[:,:b]) == DataVector{Char}

    data = {:A => [1, 2], :C => [:1, :2], :B => [3, 4]}
    df = DataFrame(data)
    # Specify column_names
    df = DataFrame(data, [:C, :A, :B])

    # This assignment was missing before
    df = DataFrame(Column = [:A])
    df[1, :Column] = "Testing"

    # zero-row dataframe and subdataframe test
    df = DataFrame(x=[], y=[])
    @test nrow(df) == 0
    df = DataFrame(x=[1:3], y=[3:5])
    sdf = sub(df, df[:x] .== 4)
    @test size(sdf, 1) == 0

    @test hash(convert(DataFrame, [1 2; 3 4])) == hash(convert(DataFrame, [1 2; 3 4]))
    @test hash(convert(DataFrame, [1 2; 3 4])) != hash(convert(DataFrame, [1 3; 2 4]))

    # push!(df, row)
    df=DataFrame( first=[1,2,3], second=["apple","orange","pear"] )

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, {3,"pear"})
    @test df==dfb

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, (3,"pear"))
    @test df==dfb

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, (33.33,"pear"))

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, ("coconut",22))

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, [ :first=>3,:second=>"pear" ])
    @test df==dfb

    df=DataFrame( first=[1,2,3], second=["apple","orange","banana"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, [ "first"=>3,"second"=>"banana" ])
    @test df==dfb

    df0= DataFrame( first=[1,2], second=["apple","orange"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, [ :first=>true,:second=>false ])
    @test df0==dfb

    df0= DataFrame( first=[1,2], second=["apple","orange"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, ["first"=>"chicken", "second"=>"stuff" ])
    @test df0==dfb

end


text:
module GeometryTests

using Test

@testset "Triangulations" begin include("TriangulationsTests.jl") end

@testset "Grids" begin include("GridsTests.jl") end

@testset "UnstructuredGrids" begin include("UnstructuredGridsTests.jl") end

@testset "CartesianGrids" begin include("CartesianGridsTests.jl") end

@testset "GridTopologies" begin include("GridTopologiesTests.jl") end

@testset "UnstructuredGridTopologies" begin include("UnstructuredGridTopologiesTests.jl") end

@testset "FaceLabelings" begin include("FaceLabelingsTests.jl") end

@testset "DiscreteModels" begin include("DiscreteModelsTests.jl") end

@testset "UnstructuredDiscreteModels" begin include("UnstructuredDiscreteModelsTests.jl") end

@testset "CartesianDiscreteModels" begin include("CartesianDiscreteModelsTests.jl") end

@testset "PeriodicBC" begin include("PeriodicBCTests.jl") end

@testset "RestrictedTriangulations" begin include("RestrictedTriangulationsTests.jl") end

@testset "BoundaryTriangulations" begin include("BoundaryTriangulationsTests.jl") end

@testset "SkeletonTriangulations" begin include("SkeletonTriangulationsTests.jl") end

@testset "GridPortions" begin include("GridPortionsTests.jl") end

@testset "DiscreteModelPortions" begin include("DiscreteModelPortionsTests.jl") end

@testset "RestrictedDiscreteModels" begin include("RestrictedDiscreteModelsTests.jl") end

@testset "AppendedTriangulations" begin include("AppendedTriangulationsTests.jl") end

#
#@testset "CellFields" begin include("CellFieldsTests.jl") end
#
#@testset "QPointCellFields" begin include("QPointCellFieldsTests.jl") end
#
#
#@testset "SkeletonPairs" begin include("SkeletonPairsTests.jl") end
#
#@testset "CellQuadratures" begin include("CellQuadraturesTests.jl") end
#

end # module


text:
immutable Binomial <: DiscreteUnivariateDistribution
    size::Int
    prob::Float64
    function Binomial(n::Real, p::Real)
    	if n <= 0
	    	error("size must be positive")
	    else
	    	if 0.0 <= p <= 1.0
	    		new(int(n), float64(p))
	    	else
	    		error("prob must be in [0, 1]")
			end
	    end
	end
end

Binomial(size::Integer) = Binomial(size, 0.5)
Binomial() = Binomial(1, 0.5)

@_jl_dist_2p Binomial binom

entropy(d::Binomial) = d.size * (-xlogx(1.0 - d.prob) - xlogx(d.prob))

insupport(d::Binomial, x::Number) = isinteger(x) && 0 <= x <= d.size

kurtosis(d::Binomial) = (1.0 - 2.0 * d.prob * (1.0 - d.prob)) / var(d)

mean(d::Binomial) = d.size * d.prob

median(d::Binomial) = iround(d.size * d.prob)

# TODO: May need to subtract 1 sometimes
modes(d::Binomial) = iround((d.size + 1.0) * d.prob)

function mgf(d::Binomial, t::Real)
	n, p = d.size, d.prob
	return (1.0 - p + p * exp(t))^n
end

function cf(d::Binomial, t::Real)
	n, p = d.size, d.prob
	return (1.0 - p + p * exp(im * t))^n
end

modes(d::Binomial) = iround([d.size * d.prob])

# TODO: rand() is totally screwed up

skewness(d::Binomial) = (1.0 - 2.0 * d.prob) / std(d)

var(d::Binomial) = d.size * d.prob * (1.0 - d.prob)


text:
module PoissonLagrangeMultiplierTests

using Test
using Gridap
using Gridap.Geometry

function is_left(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  ( x[1] <= 0.5 )
end

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
bgface_to_mask = get_face_mask(labels,"boundary",1)
Γface_to_bgface = findall(bgface_to_mask)
model_Γ = BoundaryDiscreteModel(Polytope{1},model,Γface_to_bgface)

Γ = Triangulation(model_Γ)
Γface_coords = get_cell_coordinates(Γ)
Γface_mask = lazy_map(is_left,Γface_coords)
Γlface_Γface = findall(Γface_mask)
Γrface_Γface = findall(!,Γface_mask)
Γl = BoundaryTriangulation(model,view(Γface_to_bgface,Γlface_Γface))
Γr = BoundaryTriangulation(model,view(Γface_to_bgface,Γrface_Γface))

Λ = SkeletonTriangulation(model_Γ)

order = 2
reffe_u = ReferenceFE(lagrangian,Float64,order)
reffe_λ = ReferenceFE(lagrangian,Float64,order-1)
V = TestFESpace(model,reffe_u,conformity=:H1)
S = TestFESpace(model_Γ,reffe_λ,conformity=:L2)
U = TrialFESpace(V)
L = TrialFESpace(S)
Y = MultiFieldFESpace([V,S])
X = MultiFieldFESpace([U,L])

degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
dΓl = Measure(Γl,degree)
dΓr = Measure(Γr,degree)
dΛ = Measure(Λ,degree)
n = get_normal_vector(Γ)
nΛ = get_normal_vector(Λ)

uₑ(x) = x[1]^2 + x[2]^2
f(x) = -Δ(uₑ)(x)

# Weak form. Additional non needed terms are added for testing purposes
a((u,λ),(v,η)) = ∫( ∇(u)⋅∇(v) )dΩ +
               ∫( (u+λ)*(η+v) - u*v - η*λ + 0.0*(n⋅∇(u)-λ)*(n⋅∇(v)-η) )dΓl +
               ∫( (u+λ)*(η+v) - u*v - η*λ + 0.0*(n⋅∇(u)-λ)*(n⋅∇(v)-η) )dΓr +
               ∫( 0.0*jump(n⋅∇(u))*mean(η) )dΛ +
               ∫( 0.0*jump(nΛ⋅∇(u))*mean(η) )dΛ
l((v,η)) = ∫( f*v )dΩ + ∫( uₑ*η )dΓl + ∫( uₑ*η )dΓr

op = AffineFEOperator(a,l,X,Y)
uh,λh = solve(op)

e = uh-uₑ
l2(v) = √(∑(∫(v*v)dΩ))

tol = 1.0e-10
@test l2(e) < tol

end


text:
facts("Testing creation of fixed accuracy agents pool") do
    agents = TcrCore.setupAgentsWithFixedAccuracy(3, 50)
    @fact length(agents) --> 3
    @fact agents[1].accuracy --> 50
    @fact agents[2].accuracy --> 50
    @fact agents[3].accuracy --> 50
end

facts("Testing creation of low diversified accuracy agents pool") do
    srand(1234)
    agents = TcrCore.setupRandomAgents(5, 50, 10)
    @fact length(agents) --> 5
    @fact agents[1].accuracy --> 58
    @fact agents[2].accuracy --> 40
    @fact agents[3].accuracy --> 45
    @fact agents[4].accuracy --> 40
    @fact agents[5].accuracy --> 58
end

facts("Testing creation of high diversified accuracy agents pool") do
    srand(1234)
    agents = TcrCore.setupRandomAgents(5, 50, 20)
    @fact length(agents) --> 5
    @fact agents[1].accuracy --> 67
    @fact agents[2].accuracy --> 31
    @fact agents[3].accuracy --> 40
    @fact agents[4].accuracy --> 31
    @fact agents[5].accuracy --> 67
end

facts("Testing candidate evaluation by a perfect agent") do
    srand(1234)
    agent = Agent(100)
    @fact TcrCore.evaluateCandidateByAgent(30, agent) --> 30
    @fact TcrCore.evaluateCandidateByAgent(50, agent) --> 50
    @fact TcrCore.evaluateCandidateByAgent(80, agent) --> 80
end

facts("Testing candidate evaluation by a high accuracy agent") do
    srand(1234)
    agent = Agent(80)
    @fact TcrCore.evaluateCandidateByAgent(30, agent) --> 35.78231467967497
    @fact TcrCore.evaluateCandidateByAgent(50, agent) --> 56.011625439045446
    @fact TcrCore.evaluateCandidateByAgent(80, agent) --> 83.29652502336155
end

facts("Testing candidate evaluation by a medium accuracy agent") do
    srand(1234)
    agent = Agent(50)
    @fact TcrCore.evaluateCandidateByAgent(30, agent) --> 44.45578669918743
    @fact TcrCore.evaluateCandidateByAgent(50, agent) --> 65.02906359761361
    @fact TcrCore.evaluateCandidateByAgent(80, agent) --> 88.2413125584039
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
function FieldLine(arg0::FieldLine)
    return FieldLine((FieldLine,), arg0)
end

function FieldLine(arg0::FieldVector3D, arg1::FieldVector3D, arg2::jdouble)
    return FieldLine((FieldVector3D, FieldVector3D, jdouble), arg0, arg1, arg2)
end

function closest_point(obj::FieldLine, arg0::FieldLine)
    return jcall(obj, "closestPoint", FieldVector3D, (FieldLine,), arg0)
end

function contains(obj::FieldLine, arg0::FieldVector3D)
    return jcall(obj, "contains", jboolean, (FieldVector3D,), arg0)
end

function contains(obj::FieldLine, arg0::Vector3D)
    return jcall(obj, "contains", jboolean, (Vector3D,), arg0)
end

function distance(obj::FieldLine, arg0::FieldLine)
    return jcall(obj, "distance", RealFieldElement, (FieldLine,), arg0)
end

function distance(obj::FieldLine, arg0::FieldVector3D)
    return jcall(obj, "distance", RealFieldElement, (FieldVector3D,), arg0)
end

function distance(obj::FieldLine, arg0::Vector3D)
    return jcall(obj, "distance", RealFieldElement, (Vector3D,), arg0)
end

function get_abscissa(obj::FieldLine, arg0::FieldVector3D)
    return jcall(obj, "getAbscissa", RealFieldElement, (FieldVector3D,), arg0)
end

function get_abscissa(obj::FieldLine, arg0::Vector3D)
    return jcall(obj, "getAbscissa", RealFieldElement, (Vector3D,), arg0)
end

function get_direction(obj::FieldLine)
    return jcall(obj, "getDirection", FieldVector3D, ())
end

function get_origin(obj::FieldLine)
    return jcall(obj, "getOrigin", FieldVector3D, ())
end

function get_tolerance(obj::FieldLine)
    return jcall(obj, "getTolerance", jdouble, ())
end

function intersection(obj::FieldLine, arg0::FieldLine)
    return jcall(obj, "intersection", FieldVector3D, (FieldLine,), arg0)
end

function is_similar_to(obj::FieldLine, arg0::FieldLine)
    return jcall(obj, "isSimilarTo", jboolean, (FieldLine,), arg0)
end

function point_at(obj::FieldLine, arg0::RealFieldElement)
    return jcall(obj, "pointAt", FieldVector3D, (RealFieldElement,), arg0)
end

function point_at(obj::FieldLine, arg0::jdouble)
    return jcall(obj, "pointAt", FieldVector3D, (jdouble,), arg0)
end

function reset(obj::FieldLine, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(obj, "reset", void, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function revert(obj::FieldLine)
    return jcall(obj, "revert", FieldLine, ())
end


text:
include("../src/components/climatedynamics_component.jl")

@testset "climatedynamics" begin

    Precision = 1.0e-11
    T = length(MimiDICE2010.model_years)
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, climatedynamics, :climatedynamics)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :climatedynamics, :FORC, read_params(f, "B122:BI122", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :climatedynamics, :fco22x, p[:fco22x])
    set_param!(m, :climatedynamics, :t2xco2, p[:t2xco2])
    set_param!(m, :climatedynamics, :tatm0, p[:tatm0])
    set_param!(m, :climatedynamics, :tatm1, p[:tatm1])
    set_param!(m, :climatedynamics, :tocean0, p[:tocean0])
    set_param!(m, :climatedynamics, :c1, p[:c1])
    set_param!(m, :climatedynamics, :c3, p[:c3])
    set_param!(m, :climatedynamics, :c4, p[:c4])

# Run the one-component model
    run(m)

# Extract the generated variables
    TATM    = m[:climatedynamics, :TATM]
    TOCEAN  = m[:climatedynamics, :TOCEAN]

# Extract the true values
    True_TATM   = read_params(f, "B121:BI121", T)
    True_TOCEAN = read_params(f, "B123:BI123", T)

# Test that the values are the same
    @test maximum(abs, TATM .- True_TATM) ≈ 0. atol = Precision
    @test maximum(abs, TOCEAN .- True_TOCEAN) ≈ 0. atol = Precision

end


text:
"""
Types derived from `B1pMultiplier` can be used to multiply a 1-p basis
with another function. They are implemented as another 1-p basis, but with just
a single basis function attached to them. In particular they don't get
a basis function index attached to them. The abstract type
`abstract type B1pMultiplier` implements the generic functionality
to enable this.

Derived types must implement
* `_inner_evaluate`
* `_inner_evaluate_d`
"""
abstract type B1pMultiplier{T} <: OneParticleBasis{T} end

valtype(::B1pMultiplier{T}, args...)  where {T} = T

function _inner_evaluate end
function _inner_evaluate_d end

function evaluate!(A, mult::B1pMultiplier, X::AbstractState)
    A[1] = _inner_evaluate(mult, X)
    return A
end

symbols(::B1pMultiplier) = Tuple{}()

indexrange(::B1pMultiplier) = Dict()

Base.length(basis::B1pMultiplier) = 1

isadmissible(b, ::B1pMultiplier) = true

degree(b, ::B1pMultiplier, args...) = 0

get_index(::B1pMultiplier, b) = 1

get_spec(::B1pMultiplier, n::Integer) = NamedTuple()


text:
using EffectSizes, Statistics
using Test
using EffectSizes: AbstractEffectSize, correction, pooledstd1, pooledstd2, _effectsize

@testset "$T constructor" for (T, v) = [(EffectSize, :d), (CohenD, :d), (HedgeG, :g),
                                        (GlassΔ, :Δ)]
    e = 0.
    l = -1.
    u = 1.
    q = .95
    es = T(e, ConfidenceInterval(l, u, q))
    @test es isa AbstractEffectSize
    @test effectsize(es) == getfield(es, v)
    @test confint(es) == es.ci
    @test quantile(es) == es.ci.quantile
    @test lower(es.ci) == es.ci.ci[1]
    @test upper(es.ci) == es.ci.ci[2]

    xs = randn(90)
    ys = randn(110)

    @test effectsize(T(xs .+ 1, xs)) > 0
    @test effectsize(T(xs, xs .+ 1)) < 0

    @testset "constructors" begin
        # Normal
        es = T(xs, ys)
        @test quantile(es) == 0.95
        @test typeof(effectsize(es)) == eltype(xs)
        @test es == T(xs, ys)

        @test effectsize(T(xs, xs)) == 0

        es2 = T(xs, ys, quantile=0.8)
        @test lower(confint(es2)) > lower(confint(es))
        @test upper(confint(es2)) < upper(confint(es))
        @test quantile(es2) == 0.8
        @test effectsize(es2) == effectsize(es)

        # bootstrap
        es3 = T(xs, ys, 100)
        @test effectsize(es3) == effectsize(es)
        @test quantile(es3) == 0.95

        es4 = T(xs, ys, 100, quantile=0.1)
        @test quantile(es4) == 0.1
        @test lower(confint(es3)) < lower(confint(es4))
        @test upper(confint(es3)) > upper(confint(es4))

        @test effectsize(T(xs, xs, 100)) == 0
    end

    @testset "Integers" begin
        # Smoke tests for `Int`s.
        xs = round.(Int, randn(90) .* 10)
        ys = round.(Int, randn(110) .* 10)
        T(xs, ys, quantile=0.95)
    end
end

@testset "correction" begin
    @test_throws DomainError correction(1)
    @test correction(2) == 0.
    @test correction(20) > .9
    @test correction(10^9) ≈ 1
end

@testset "_effectsize" begin
    # from http://staff.bath.ac.uk/pssiw/stats2/page2/page14/page14.html
    @test round(_effectsize(20., 24., 4.53), digits=3) == -0.883

    xs = 10:110
    ys = 1:100
    mx = mean(xs)
    my = mean(ys)
    s = pooledstd1(xs, ys)
    @test _effectsize(mx, my, s, 100) > _effectsize(mx, my, s, 50) > _effectsize(mx, my, s, 10)
end


text: # import Merlin


text:
module ParallelUtils

export @ensureworkers, sendto, getfrom

macro ensureworkers(numWorkers)
    :(if nworkers() < $numWorkers
        if nprocs() == 1
            addprocs($numWorkers)
            @everywhere import ParallelUtils
        else
            addprocs($numWorkers-nworkers())
            @everywhere import ParallelUtils
        end
    end)
end

# from http://stackoverflow.com/questions/27677399/julia-how-to-copy-data-to-another-processor-in-julia
function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, eval(Main, Expr(:(=), nm, val)))
    end
end
function sendto(ps::Vector{Int}; args...)
    for p in ps
        sendto(p; args...)
    end
end
getfrom(p::Int, nm::Symbol; mod=Main) = fetch(@spawnat(p, getfield(mod, nm)))

"Save the type and data of the passed dictionary as text."
function writedict(f::IOStream, dict::Dict)
    println(f, join(match(r"[A-z]+{([^,]+),([^}]+)}", string(typeof(dict))).captures, '\t'))
    for (k,v) in dict
        println(f, string(k)*"\t"*string(v))
    end
end
function writedict(fileName::AbstractString, dict::Dict)
    f = open(fileName, "w")
    writedict(f, dict)
    close(f)
end

"""
Read the type and data of the passed dictionary.

Numeric types are converted automatically, and string type are
preserved, but other types may fail to load.
"""
function readdict(f::IOStream; reversed=false)
    parts = split(strip(readline(f)), '\t')
    type1 = eval(parse(parts[1]))
    type2 = eval(parse(parts[2]))
    d = reversed ? Dict{type2,type1}() : Dict{type1,type2}()
    for line in eachline(f)
        k,v = split(strip(line), '\t', limit=2)
        if type1 <: Number
            k = parse(type1, k)
        end
        if type2 <: Number
            v = parse(type2, v)
        end

        reversed ? d[v] = k : d[k] = v
    end
    d
end
function readdict(fileName::AbstractString; reversed=false)
    f = open(fileName)
    d = readdict(f, reversed=reversed)
    close(f)
    d
end

end # module

max_stars_repo_name: slundberg/ParallelUtils

max_issues_repo_name: slundberg/ParallelUtils

max_forks_repo_name: slundberg/ParallelUtils


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

using BAT
using Test

using Random, Distributions, StatsBase

@testset "bat_sample" begin
    @testset "RandSampling" begin
        dist = MvNormal([0.4, 0.6], [2.0 1.2; 1.2 3.0])

        @test length(@inferred(bat_sample(dist, 1)).samples) == 1
        @test length(@inferred(bat_sample(dist, 10^3)).samples) == 10^3

        @test length(@inferred(bat_sample(Random.GLOBAL_RNG, dist, 10^2)).samples) == 10^2
        @test length(@inferred(bat_sample(dist, 10^2, BAT.RandSampling())).samples) == 10^2
        @test length(@inferred(bat_sample(Random.GLOBAL_RNG, dist, 10^2, BAT.RandSampling())).samples) == 10^2

        samples, stats = bat_sample(dist, 10^5)
        @test isapprox(mean(samples.params), [0.4, 0.6]; rtol = 0.05)
        @test isapprox(cov(samples.params), [2.0 1.2; 1.2 3.0]; rtol = 0.05)
        @test all(isequal(1), samples.weight)
        @test stats isa NamedTuple{(:mode,:mean,:cov)}
        @test stats == bat_stats(samples)
    end
end


text:
"""
    boruvka_mst(g, distmx = weights(g); minimize = true)

Return a tuple `(mst, weights)` where `mst` is a vector of edges representing the
optimum (minimum, by default) spanning tree of a connected, undirected graph
`g` with optional matrix `distmx` that provides distinct edge weights, and
`weights` is the sum of all the edges in the solution by using
[Boruvka's algorithm](https://en.wikipedia.org/wiki/Bor%C5%AFvka%27s_algorithm).
The algorithm requires that all edges have different weights to correctly generate a minimun/maximum spanning tree
### Optional Arguments
- `minimize=true`: if set to `false`, calculate the maximum spanning tree.
"""
function boruvka_mst end

@traitfn function boruvka_mst(g::AG::(!IsDirected),
        distmx::AbstractMatrix{T} = weights(g);
        minimize = true) where {T<:Real, U, AG<:AbstractGraph{U}}

    djset = IntDisjointSets(nv(g))
    # maximizing Z is the same as minimizing -Z
    # mode will indicate the need for the -1 multiplication
    mode = minimize ? 1 : -1
    mst = Vector{edgetype(g)}()
    sizehint!(mst, nv(g) - 1)
    weight = zero(T)

    while true
        cheapest = Vector{Union{edgetype(g), Nothing}}(nothing, nv(g))
        # find cheapest edge that connects two components
        found_edge = false
        for edge in edges(g)
            set1 = find_root(djset, src(edge))
            set2 = find_root(djset, dst(edge))
            if set1 != set2
                found_edge = true
                e1 = cheapest[set1]
                if e1 === nothing || distmx[src(e1), dst(e1)] * mode > distmx[src(edge), dst(edge)] * mode
                    cheapest[set1] = edge
                end
                e2 = cheapest[set2]
                if e2===nothing || distmx[src(e2), dst(e2)] * mode > distmx[src(edge), dst(edge)] * mode
                    cheapest[set2] = edge
                end
            end
        end
        #no more edges between two components
        !found_edge && break
        # add cheapest edges to the tree
        for v in vertices(g)
            if cheapest[v] !== nothing
                edge = cheapest[v]
                if !in_same_set(djset, src(edge), dst(edge))
                    weight += distmx[src(edge), dst(edge)]
                    union!(djset, src(edge), dst(edge))
                    push!(mst, edge)
                end
            end
        end
    end
    return (mst=mst, weight=weight)
end