text:
# Represents types of Query Results
immutable QueryResultType name::String end
const VertexSetQueryResult = QueryResultType("VertexSetQueryResult")
const EdgeSetQueryResult = QueryResultType("EdgeSetQueryResult")
const InitialQueryResult = QueryResultType("InitialQueryResult")

# Represents options for the selection of query results
immutable QueryResultOption name::String end
const IncomingQueryResultOption = QueryResultOption("IncomingQueryResultOption")
const OutgoingQueryResultOption = QueryResultOption("OutgoingQueryResultOption")
const HeadQueryResultOption = QueryResultOption("HeadQueryResultOption")
const TailQueryResultOption = QueryResultOption("TailQueryResultOption")
const NoneQueryResultOption = QueryResultOption("NoneQueryResultOption")
const DistinctResultOption = QueryResultOption("DistinctResultOption")
const MoveToVerticesResultOption = QueryResultOption("MoveToVerticesResultOption")
const MoveToEdgesResultOption = QueryResultOption("MoveToEdgesResultOption")
const StoreResultOption = QueryResultOption("StoreResultOption")

type CompositeQuery
	# Represents a Query in its entirety, or a portion of a larger Query

	graph::Graph
	previous::CompositeQuery
	associatedfunction::Function
	associatedexpression::Expr
	inputtype::QueryResultType
	outputtype::QueryResultType
	option::QueryResultOption
	depth::Integer
	result::Any
	isrealised::Bool
	isgrouped::Bool
	repeatcount::Integer
	operationstepback::Integer
	containspartiallyloadedresults::Bool

	# loader declared Any instead of GraphLoader
	# in order to avoid mutually circular type declaration
	# see: https://github.com/JuliaLang/julia/issues/269
	loader::Any

	function CompositeQuery(g::Graph, outputtype::QueryResultType, option::QueryResultOption)
		# Constructs a CompositeQuery for the specified Graph

		cq = new()
		cq.graph = g
		cq.inputtype = InitialQueryResult
		cq.outputtype = outputtype
		cq.option = option
		cq.depth = 1
		cq.isrealised = false
		cq.isgrouped = false
		cq.result = Set{CompositeQueryResultItem}()
		cq.repeatcount = 0
		cq.operationstepback = 0
		cq.containspartiallyloadedresults = false

		if isdefined(g, :loader)
			cq.loader = g.loader
		end

		return cq
	end

	function CompositeQuery(g::Graph, outputtype::QueryResultType, option::QueryResultOption, associatedfunction::Function)
		# Constructs a CompositeQuery for the specified Graph with a specified where filter

		cq = CompositeQuery(g, outputtype, option)
		cq.associatedfunction = associatedfunction

		return cq
	end

	function CompositeQuery(g::Graph, outputtype::QueryResultType, option::QueryResultOption, associatedexpression::Expr)
		# Constructs a CompositeQuery for the specified Graph with a specified where filter

		cq = CompositeQuery(g, outputtype, option)
		cq.associatedexpression = associatedexpression

		return cq
	end

	function CompositeQuery(previous::CompositeQuery, outputtype::QueryResultType, option::QueryResultOption)
		# Constructs a CompositeQuery which extends a previous query

		cq = CompositeQuery(previous.graph, outputtype, option)
		cq.previous = previous
		cq.inputtype = previous.outputtype
		cq.outputtype = outputtype
		cq.option = option
		cq.depth = previous.depth + 1
		cq.repeatcount = previous.repeatcount

		if isdefined(previous, :loader)
			cq.loader = previous.loader
		end

		return cq
	end

	function CompositeQuery(previous::CompositeQuery, outputtype::QueryResultType, option::QueryResultOption, associatedfunction::Function)
		# Constructs a CompositeQuery which extends a previous query and has an initial where filter function

		cq = CompositeQuery(previous, outputtype, option)
		cq.associatedfunction = associatedfunction

		return cq
	end

	function CompositeQuery(previous::CompositeQuery, outputtype::QueryResultType, option::QueryResultOption, associatedexpression::Expr)
		# Constructs a CompositeQuery which extends a previous query and has an initial where filter expression

		cq = CompositeQuery(previous, outputtype, option)
		cq.associatedexpression = associatedexpression

		return cq
	end
end

function vertices(g::Graph)
	# Operation that selects the vertices of a graph.  A CompositeQuery which begins with this operation is constructed.

	return CompositeQuery(g, VertexSetQueryResult, NoneQueryResultOption)
end

function vertices(g::Graph, where::Function)
	# Operation that selects the vertices of a graph that match a supplied where condition.  A CompositeQuery which begins with this operation is constructed.

	return CompositeQuery(g, VertexSetQueryResult, NoneQueryResultOption, where)
end

function edges(g::Graph)
	# Operation that selects the edges of a graph.  A CompositeQuery which begins with this operation is constructed.

	return CompositeQuery(g, EdgeSetQueryResult, NoneQueryResultOption)
end

function edges(g::Graph, where::Function)
	# Operation that selects the edges of a graph that match a specified where condition.  A CompositeQuery which begins with this operation is constructed.

	return CompositeQuery(g, EdgeSetQueryResult, NoneQueryResultOption, where)
end

function outgoing(cq::CompositeQuery)
	# Operation that selects the outgoing edges of the Nodes resulting from the supplied query.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, EdgeSetQueryResult, OutgoingQueryResultOption)
end

function outgoing(cq::CompositeQuery, where::Function)
	# Operation that selects the outgoing edges of the Nodes resulting from the supplied query that match the supplied where condition.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, EdgeSetQueryResult, OutgoingQueryResultOption, where)
end

function incoming(cq::CompositeQuery)
	# Operation that selects the incoming edges of the Nodes resulting from the supplied query.  A CompositeQuery which extends the previous query to include this operation is returned.
	return CompositeQuery(cq, EdgeSetQueryResult, IncomingQueryResultOption)
end

function incoming(cq::CompositeQuery, where::Function)
	# Operation that selects the incoming edges of the Nodes resulting from the supplied query that match the supplied where function.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, EdgeSetQueryResult, IncomingQueryResultOption, where)
end

function head(cq::CompositeQuery)
	# Operation that selects the Vertices at the head of the edges resulting from the supplied query.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, VertexSetQueryResult, HeadQueryResultOption)
end

function head(cq::CompositeQuery, where::Function)
	# Operation that selects the Vertices at the head of the edges resulting from the supplied query that match the supplied where function.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, VertexSetQueryResult, HeadQueryResultOption, where)
end

function tail(cq::CompositeQuery)
	# Operation that selects the Vertices at the tail of the edges resulting from the supplied query.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, VertexSetQueryResult, TailQueryResultOption)
end

function tail(cq::CompositeQuery, where::Function)
	# Operation that selects the Vertices at the tail of the edges resulting from the supplied query that match the supplied where function.  A CompositeQuery which extends the previous query to include this operation is returned.

	return CompositeQuery(cq, VertexSetQueryResult, TailQueryResultOption, where)
end

function distinct(cq::CompositeQuery)
	# Operation that filters the current set of results to a distinct set
	return CompositeQuery(cq, cq.outputtype, DistinctResultOption)
end

function store(cq::CompositeQuery, select::Function)
	# Operation that stores data to be carried forward in the query
	return CompositeQuery(cq, cq.outputtype, StoreResultOption, select)
end

function movetoedges(cq::CompositeQuery, select::Function)
	# Operation that moves to a previously stored set of edges
	return CompositeQuery(cq, EdgeSetQueryResult, MoveToEdgesResultOption, select)
end

function movetovertices(cq::CompositeQuery, select::Function)
	# Operation that moves to a previously stored set of edges
	return CompositeQuery(cq, VertexSetQueryResult, MoveToVerticesResultOption, select)
end

function select(cq::CompositeQuery, select::Function)
	# executes the query and returns a set of mapped results
	if !cq.isrealised
		realise!(cq)
	end

	results = Set{Any}()

	if isspecified(cq.result)
		for r in cq.result
			push!(results, select(r))
		end
	end

	return results
end

function reduce(reducer::Function, cq::CompositeQuery, mapper::Function)
	# Performs a map reduce operation on the Edges or Vertices resulting from a query

	if !cq.isrealised
		realise!(cq)
	end

	if !isdefined(cq, :result)
		throw (CompositeQueryResultNotRealisedException())
	end

	currentvalue::Any = UnspecifiedValue
	if isspecified(cq.result)
		for r in cq.result
			currentvalue = reducer(currentvalue, mapper(r))
		end
	end

	return currentvalue
end

function sum(cq::CompositeQuery, map::Function)
	# Sums the values mapped from the Vertices or Edges resulting from the supplied query by using the supplied map function

	return reduce(cq, map) do x,y
			if isspecified(x) && isspecified(y)
				x + y
			else
				y
			end
		end
end

function maximum(cq::CompositeQuery, map::Function)
	# Returns the maximum of the values mapped from the Vertices or Edges resulting from the supplied query using the supplied map function

	return reduce(cq, map) do x,y
			if isspecified(x) && isspecified(y) && x > y
				x
			else
				y
			end
		end
end

function minimum(cq::CompositeQuery, map::Function)
	# Returns the minimum of the values mapped from the Vertices or Edges resulting from the supplied query by using the supplied map function

	return reduce(cq, map) do x,y
			if isspecified(x) && isspecified(y) && x < y
				x
			else
				y
			end
		end
end

function count(cq::CompositeQuery)
	# Returns the count of the  Vertices or Edges resulting from the supplied query

	result = reduce(cq, n -> 1) do x, y
			if isspecified(x) && isspecified(y)
				x + y
			else
				y
			end
		end

	if !isspecified(result)
		result = 0
	end

	return result
end

function average(cq::CompositeQuery, map::Function)
	# Returns the average of the  values mapped from the Vertices or Edges resulting from the supplied query by using the supplied map function

	result = UnspecifiedValue

	countresult = count(cq)

	if countresult > 0
		sum_result = sum(cq, map)
		result = sum_result / countresult
	end

	return result
end

function group(cq::CompositeQuery, groupby::Function)
	# groups the results of a composite query
	if !cq.isrealised
		realise!(cq)
	end

	if !isdefined(cq, :result)
		throw (CompositeQueryResultNotRealisedException())
	end

	groupedresult = Dict{Any, CompositeQuery}()

	if isspecified(cq.result)
		for r in cq.result
			key = groupby(r)

			groupset = Base.get(groupedresult, key, UnspecifiedValue)

			if groupset == UnspecifiedValue
				groupset = CompositeQuery(cq, cq.outputtype, NoneQueryResultOption, groupby)
				groupedresult[key] = groupset
				groupset.isrealised = true
			end

			push!(groupset.result, r)
		end
	end

	return groupedresult
end

function isspecified(value::Any)
	# Tests whether a value is specified (whether it matches the UnspecifiedValue)

	return value != UnspecifiedValue
end

function getdistinctresults(cq)
	# An operation that selects the distinct set of Vertices or Edges resulting from the previous query.  A new CompositeQuery is constructed which includes the distinct operation.
	if !cq.isrealised
		realise!(cq)
	end

	if !isdefined(cq, :result)
		throw (CompositeQueryResultNotRealisedException())
	end

	# the implementation of this method based on the implementation of unique in base library
	distinctset = Set{Container}()
	newresults = Set{CompositeQueryResultItem}()

	if isspecified(cq.result)
		for r in cq.result
			if !in(r.item, distinctset)
				push!(distinctset, r.item)
				push!(newresults, CompositeQueryResultItem(r.item, r))
			end
		end
	end

	return newresults
end

function loop(cq::CompositeQuery, stepfunction::Function)
	realise!(cq)

	newresults = Set{CompositeQueryResultItem}()

	if isspecified(cq.result)
		for r in cq.result
			newitem = CompositeQueryResultItem(r.item, r)
			# copy forward the isnew flag on repeat
			newitem.isnew = r.isnew
			push!(newresults, newitem)
		end
	end

	cqnew = CompositeQuery(cq, cq.outputtype, NoneQueryResultOption)
	cqnew.operationstepback = stepfunction(cqnew)
	cqnew.result = newresults
	cqnew.isrealised = true
	cqnew.containspartiallyloadedresults = cq.containspartiallyloadedresults

	if cqnew.operationstepback > 0
		cqnew.repeatcount = cqnew.repeatcount + 1
	else
		cqnew.repeatcount = 0
	end

	return cqnew
end

function mergedistinct(cq::CompositeQuery, stepfunction::Function)
	return merge(cq, stepfunction, true)
end

function merge(cq::CompositeQuery, stepfunction::Function)
	return merge(cq, stepfunction, false)
end

function merge(cq::CompositeQuery, stepfunction::Function, distinct::Bool)
	realise!(cq)

	stepcount = stepfunction(cq)

	distinctset = Set{Container}()
	newresults = Set{CompositeQueryResultItem}()
	cqnew = CompositeQuery(cq, cq.outputtype, NoneQueryResultOption)
	cqnew.containspartiallyloadedresults = cq.containspartiallyloadedresults

	if stepcount > 0
		source = cqnew
		stepped = 0
		for i in 1:stepcount
			if isdefined(source, :previous)
				source = source.previous
				stepped = stepped + 1
			else
				# throw
				throw(InvalidStepCountForMergeException())
			end
		end

		cqnew.containspartiallyloadedresults = cqnew.containspartiallyloadedresults || source.containspartiallyloadedresults

		if isspecified(source.result)
			for r in source.result
				add = !distinct
				if distinct
					if !in(r.item, distinctset)
						push!(distinctset, r.item)
						add = true
					end
				end

				if add
					mergeitem = CompositeQueryResultItem(r.item, r)
					mergeitem.isnew = false
					push!(newresults, mergeitem)
				end
			end
		end
	end

	if isspecified(cq.result)
		for r in cq.result
			add = !distinct
			if distinct
				if !in(r.item, distinctset)
					push!(distinctset, r.item)
					add = true
				end
			end

			if add
				newitem = CompositeQueryResultItem(r.item, r)
				push!(newresults, newitem)
			end
		end
	end

	cqnew.result = newresults
	cqnew.isrealised = true

	return cqnew
end

function filter(cq::CompositeQuery, filterfunction::Function)
	# Operation that filters results without traversing
	return CompositeQuery(cq, cq.outputtype, NoneQueryResultOption, filterfunction)
end

function newonly(cq::CompositeQuery)
	# Operation that filters results without traversing
	return filter(cq, i->i.isnew)
end

function resultcount(cq::CompositeQuery)
	realise!(cq)

	return length(cq.result)
end

function previousresultcount(cq::CompositeQuery, steps::Integer)
	source = cq
	stepped = 0
	for i in 1:steps
		if isdefined(source, :previous)
			source = source.previous
			stepped = stepped + 1
		else
			# throw
			throw(InvalidStepCountForResultCountException())
		end
	end
	realise!(source)
	return length(source.result)
end

# need methods
#   mergedistinct
#   latestmergecount

function realise!(cq::CompositeQuery)
	ispartial = false
	sourceisprefiltered = false

	# realise each part of the composite query in turn
	# there is alot of scope for optimisation here to combine parts of the composite query
	if !cq.isrealised
		# ensure previous parts of the query are realised first
		if isdefined(cq, :previous) && !cq.previous.isrealised
			realise!(cq.previous)
		end

		if cq.outputtype != EdgeSetQueryResult && cq.outputtype != VertexSetQueryResult
			# throw
			throw(InvalidOutputQueryResultTypeException())
		end

		source = Set{CompositeQueryResultItem}()

		if cq.inputtype == InitialQueryResult
			items = Set{Container}()

			if isdefined(cq, :loader)
				loadfunction = loadedges
				if cq.outputtype == EdgeSetQueryResult
					# do nothing
				elseif cq.outputtype == VertexSetQueryResult
					loadfunction = loadvertices
				else
					throw(InvalidInputQueryResultTypeException())
				end

				if isdefined(cq, :associatedfunction)
					sourceisprefiltered = true
					items = loadfunction(cq.loader, cq.associatedfunction)
				elseif isdefined(cq, :associatedexpression)
					sourceisprefiltered = true
					items = loadfunction(cq.loader, cq.associatedexpression)
				else
					items = loadfunction(cq.loader)
				end
			else
				if cq.outputtype == EdgeSetQueryResult
					items = values(cq.graph.edges)
				elseif cq.outputtype == VertexSetQueryResult
					items = values(cq.graph.vertices)
				else
					throw(InvalidInputQueryResultTypeException())
				end
			end

			for i in items
				ispartial = ispartial || i.partiallyloaded
				push!(source, CompositeQueryResultItem(i))
			end
		else
			if !isdefined(cq, :previous)
				# throw
				throw(MissingPreviousQueryContextException())
			end

			if !isdefined(cq.previous, :result)
				# throw
				throw(MissingPreviousQueryContextResultException())
			end

			source = cq.previous.result

			if cq.previous.containspartiallyloadedresults
				if isdefined(cq, :loader)
					# load missing data from the previous results
					loadmissingdata(cq.loader, source)
					cq.previous.containspartiallyloadedresults = false
				end
			end
		end

		if cq.option == StoreResultOption
			for r in source
				newresultitem = CompositeQueryResultItem(r.item, r)
				push!(cq.result, newresultitem)

				values = cq.associatedfunction(newresultitem)
				if isspecified(values)
					merge!(newresultitem.storedvalues, values)
				end
			end
		elseif cq.option == MoveToVerticesResultOption
			for r in source
				item = cq.associatedfunction(r)
				if isspecified(item)
					if isa(item,Vertex)
						ispartial = ispartial || item.partiallyloaded
						newresultitem = CompositeQueryResultItem(item, r)
						push!(cq.result, newresultitem)
					else
						throw(UnexpectedTypeOnMoveToVertices())
					end
				end
			end
		elseif cq.option == MoveToEdgesResultOption
			for r in source
				item = cq.associatedfunction(r)
				if isspecified(item)
					if isa(item,Edge)
						ispartial = ispartial || item.partiallyloaded
						newresultitem = CompositeQueryResultItem(item, r)
						push!(cq.result, newresultitem)
					else
						throw(UnexpectedTypeOnMoveToEdges())
					end
				end
			end
		elseif cq.option == DistinctResultOption
			cq.result = getdistinctresults(cq.previous)
		else
			# determine whether a further check for missing data is required
			loadcheckneeded = isdefined(cq, :loader) && (isdefined(cq, :associatedfunction) || isdefined(cq, :associatedexpression) && !sourceisprefiltered)

			# default the getitem function to returning the item given
			getitem = r -> r.item
			getitemsubset = UnspecifiedValue

			if cq.inputtype == cq.outputtype || cq.inputtype == InitialQueryResult
				# nothing to do
				loadcheckneeded = false
			elseif cq.inputtype == EdgeSetQueryResult && cq.outputtype == VertexSetQueryResult
				if cq.option == HeadQueryResultOption
					getitem = (r -> r.item.head)
				elseif cq.option == TailQueryResultOption
					getitem = (r -> r.item.tail)
				else
					# throw
					throw(InvalidQueryOptionOutputTypeCombinationException())
				end
			elseif cq.inputtype == VertexSetQueryResult && cq.outputtype == EdgeSetQueryResult
				if cq.option == IncomingQueryResultOption
					getitemsubset = (r->r.incomingedges)
				elseif cq.option == OutgoingQueryResultOption
					getitemsubset = (r->r.outgoingedges)
				else
					# throw
					throw(InvalidQueryOptionOutputTypeCombinationException())
				end
			else
				#throw
				throw(InvalidInputOutputQueryOptionCombinationException())
			end

			if loadcheckneeded
				# as a function or expression is being used, ensure the related items are fully loaded before being tested
				# build up a list of the items that must be fully loaded
				loadrequired = false
				itemstoload = Set{Container}()

				if getitemsubset != UnspecifiedValue
					for r in source
						item = getitem(r)
						subset = getitemsubset(item)
						for si in subset
							if si.partiallyloaded
								push!(itemstoload, si)
								loadrequired = true
							end
						end
					end
				else
					for r in source
						item = getitem(r)
						if item.partiallyloaded
							push!(itemstoload, item)
							loadrequired = true
						end
					end
				end

				if loadrequired
					# load the missing data now if any was detected
					loadmissingdata(cq.loader, itemstoload)
				end
			end

			for r in source
				item = getitem(r)
				if getitemsubset != UnspecifiedValue
					subset = getitemsubset(item)
					for si in subset
						cqi = CompositeQueryResultItem(si, r)
						if doesitemmatchquery(cqi, cq)
							ispartial = ispartial || item.partiallyloaded
							push!(cq.result, cqi)
						end
					end
				else
					cqi = CompositeQueryResultItem(item, r)
					if sourceisprefiltered || doesitemmatchquery(cqi, cq)
						ispartial = ispartial || item.partiallyloaded
						push!(cq.result, cqi)
					end
				end
			end
		end

		cq.containspartiallyloadedresults = ispartial
		cq.isrealised = true
	end
end

function doesitemmatchquery(resultitem::CompositeQueryResultItem, cq::CompositeQuery)
	# test whether a CompositeQueryResultItem matches the conditions for inclusion in a query
	match = true

	if isdefined(cq, :associatedfunction)
		match = cq.associatedfunction(resultitem)
	end

	return match
end

# A union of the types that may be used as a source of queries
QuerySource = Union(Graph,GraphLoader,CompositeQuery,Dict{Any,CompositeQuery})

# A union of the types that may be used as operations of a query
QueryOperation = Union(Function,(Function, Function),(Function, Expr))

function query(source::QuerySource,operations::QueryOperation...)
	# Builds a query from a source and series of operations

	current = source
	isgrouped = false
	operationindex = 1

	while operationindex <= length(operations) && !isgrouped
		operation = operations[operationindex]

		# if this is a group operation break out
		if isa(current, Dict{Any,CompositeQuery})
			isgrouped = true
		else
			if isa(operation,Function)
				current = operation(current)
			elseif isa(operation, (Function,Function))
				current = operation[1](current, operation[2])
			elseif isa(operation,(Function,Expr))
				current = operation[1](current, operation[2])
			end

			if isa(current, CompositeQuery) && current.operationstepback > 0
				operationindex = operationindex - current.operationstepback
			else
				operationindex = operationindex + 1
			end
		end
	end

	if isgrouped && length(operations) >= operationindex
		remainingoperations = operations[operationindex:length(operations)]

		groupedresult = Dict{Any, Any}()

		for (k,v) in current
			push!(groupedresult, k, query(v, remainingoperations...))
		end

		return groupedresult
	else
		return current
	end
end


text:
using Nabla: checkpoint

function _checkpoint_foo(x)
    y = checkpoint(sin, (x,))
    return cos(y)
end
_foo(x) = cos(sin(x))

_bar(x; y=5.0) = sin(x * y)

@testset "checkpointing" begin

    # Positional args testing.
    @test checkpoint(sin, (5.0,)) == sin(5.0)
    @test ∇(x->checkpoint(sin, (x,)))(5.0) == (cos(5.0),)

    @test _checkpoint_foo(5.0) == _foo(5.0)
    @test ∇(_checkpoint_foo)(5.0) == ∇(_foo)(5.0)

    # kwargs testing.
    @test checkpoint(_bar, (5.0,), (y=4.5,)) == _bar(5.0; y=4.5)
    @test ∇(x->checkpoint(_bar, (x,), (y=4.5,)))(5.0) == ∇(x->_bar(x; y=4.5))(5.0)

    # Check that closures error.
    y = 5.0
    bar(x) = x + y
    @test_throws ErrorException checkpoint(bar, (5.0,))
end

ext: jl


text:
using Graphene
using Test
using DataFrames
using Images: Gray
using FixedPointNumbers
using HDF5: h5read

test_csv_data = import_csv("test_files/test_data_flower.csv")
@test isa(test_csv_data, DataFrame)
@test ncol(test_csv_data) == 2

@test isa(read_nnimage("test_files/test_nnimage_flower.jpeg"), Array{Gray{FixedPointNumbers.Normed{UInt8,8}},2})

test_nnResult_stack = import_stack("test_files/graphene_defect_nnOutput.h5", range=1:10)
@test size(test_nnResult_stack) == (256, 256, 10)
test_nnResult_stack = import_stack("test_files/graphene_defect_nnOutput.h5")
@test size(test_nnResult_stack) == (256, 256, 100)

test_centroids = make_centroids(test_nnResult_stack[:,:,1])
@test isa(test_centroids, Array{Tuple{Float64,Float64},1})

@test isa(centroid2xy(test_centroids), Array{Float64,2})
@test size(centroid2xy(test_centroids)) == (2, 837)

@test isa(stack2xy(test_nnResult_stack), Array{Array{Float64,2},1})


text:
using Gen
import Random
using FunctionalCollections: PersistentVector

include("static_model.jl")
include("dataset.jl")

@gen (static) function is_outlier_proposal(trace, i::Int)
    prev_z::Bool = trace[:data => i => :z]
    @trace(bernoulli(prev_z ? 0.0 : 1.0), :data => i => :z)
end

Gen.load_generated_functions()

line_selection = StaticSelection(select(:slope, :intercept))
std_selection = StaticSelection(select(:log_inlier_std, :log_outlier_std))

function do_inference(xs, ys, num_iters)

    observations = choicemap()
    for (i, y) in enumerate(ys)
        observations[:data => i => :y] = y
    end
    observations[:log_inlier_std] = 0.
    observations[:log_outlier_std] = 0.

    # initial trace
    (trace, _) = generate(model, (xs,), observations)

    scores = Vector{Float64}(undef, num_iters)
    for i=1:num_iters
        trace, = hmc(trace, line_selection; eps=1e-2, L=10)
        trace, = hmc(trace, std_selection; eps=1e-2, L=10)
        trace, = mala(trace, line_selection, 0.001)
        trace, = mala(trace, std_selection, 0.001)

        # step on the outliers
        for j=1:length(xs)
            (trace, _) = metropolis_hastings(trace, is_outlier_proposal, (j,))
        end

        score = get_score(trace)
        scores[i] = score

        # print
        slope = trace[:slope]
        intercept = trace[:intercept]
        inlier_std = exp(trace[:log_inlier_std])
        outlier_std = exp(trace[:log_outlier_std])
        println("score: $score, slope: $slope, intercept: $intercept, inlier_std: $inlier_std, outlier_std: $outlier_std")
    end
    return scores
end

(xs, ys) = make_data_set(200)
do_inference(xs, ys, 10)
@time scores = do_inference(xs, ys, 100)


text:
############################################################
# judiPDE ####################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiPDE, judiPDEexception, judiPDEadjoint, judiPDEadjointException, subsample

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

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections

# Forward
struct judiPDE{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    geometry::Geometry
    options::Options
    fop::Function              # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiPDEexception <: Exception
    msg :: String
end

# Adjoint
struct judiPDEadjoint{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    geometry::Geometry
    options::Options
    fop::Function              # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiPDEadjointException <: Exception
    msg :: String
end

############################################################
## outer constructors

function judiPDE(name::String,info::Info,model::Model, geometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
    (DDT == Float32 && RDT == Float32) || throw(judiPDEexception("Domain and range types not supported"))
    m = n_samples(geometry, info)
    n = info.n * sum(info.nt)
    srcnum = 1:info.nsrc

    return F = judiPDE{Float32,Float32}(name, m, n, info, model, geometry, options,
                              args -> time_modeling(model, args[1], args[2], args[3], args[4], args[5], srcnum, 'F', 1, options),
                              args_T -> time_modeling(model, args_T[1], args_T[2], args_T[3], args_T[4], args_T[5], srcnum, 'F', -1, options)
                              )
end

function judiPDEadjoint(name::String,info::Info,model::Model, geometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
    (DDT == Float32 && RDT == Float32) || throw(judiPDEadjointException("Domain and range types not supported"))
    m = n_samples(geometry, info)
    n = info.n * sum(info.nt)

    srcnum = 1:info.nsrc

    return F = judiPDEadjoint{Float32,Float32}(name, m, n, info, model, geometry, options,
                                args_T -> time_modeling(model, args_T[1], args_T[2], args_T[3], args_T[4], args_T[5], srcnum, 'F', -1, options),
                                args -> time_modeling(model, args[1], args[2], args[3], args[4], args[5], srcnum, 'F', 1, options),
                                )
end

############################################################
## overloaded Base functions

# conj(judiPDE)
conj(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        A.fop,
        A.fop_T
        )

# transpose(jo)
transpose(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
    judiPDEadjoint{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
        A.fop_T,
        A.fop
        )

# adjoint(jo)
adjoint(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
	judiPDEadjoint{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
		A.fop_T,
		A.fop
		)

# conj(jo)
conj(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT}=
    judiPDEadjoint{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        A.fop,
        A.fop_T
        )

# transpose(jo)
transpose(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
        A.fop_T,
        A.fop
        )
# adjoint(jo)
adjoint(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
        A.fop_T,
        A.fop
        )

############################################################
## overloaded Base *(...judi...)

# *(judiPDE,judiRHS)
function *(A::judiPDE{ADDT,ARDT},v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiPDEexception("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiRHS):",A.name,typeof(A),vDT]," / "))
	args = (v.geometry, v.data, A.geometry, nothing, nothing)
    V = A.fop(args)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

function *(A::judiPDEadjoint{ADDT,ARDT},v::judiRHS{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiPDEadjointException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiRHS):",A.name,typeof(A),vDT]," / "))
	args = (A.geometry,nothing,v.geometry,v.data,nothing)
    V = A.fop(args)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(judiPDE,judiWavefield)
function *(A::judiPDE{ADDT,ARDT},v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiPDEexception("Shape mismatch: A:$(size(A)), v: $(size(v))"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiWavefield):",A.name,typeof(A),vDT]," / "))
	args = (nothing,v.data,A.geometry,nothing,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

function *(A::judiPDEadjoint{ADDT,ARDT},v::judiWavefield{vDT}) where {ADDT,ARDT,vDT}
	A.n == size(v,1) || throw(judiPDEadjointException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDE,judiWavefield):",A.name,typeof(A),vDT]," / "))
	args = (A.geometry,nothing,nothing,v.data,nothing)
	V = A.fop(args)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiVector):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(judiPDE,judiProjection)
function *(A::judiPDE{CDT,ARDT},B::judiProjection{BDDT,CDT}) where {ARDT,BDDT,CDT}
    A.n == size(B,1) || throw(judiPDEexception("Shape mismatch: A:$(size(A)), B: $(size(B))"))
    return judiModeling(A.info,A.model,B.geometry,A.geometry;options=A.options,DDT=CDT,RDT=ARDT)
end

function *(A::judiPDEadjoint{CDT,ARDT},B::judiProjection{BDDT,CDT}) where {ARDT,BDDT,CDT}
    A.n == size(B,1) || throw(judiPDEadjointException("Shape mismatch: A:$(size(A)), B: $(size(B))"))
    return adjoint(judiModeling(A.info,A.model,A.geometry,B.geometry,options=A.options,DDT=CDT,RDT=ARDT))
end

# *(judiPDE,judiLRWF)
function *(A::judiPDE{CDT,ARDT},B::judiLRWF{BDDT,CDT}) where {ARDT,BDDT,CDT}
    A.n == size(B,1) || throw(judiPDEexception("Shape mismatch: A:$(size(A)), B: $(size(B))"))
    return judiPDEextended(A.info,A.model,B.wavelet,A.geometry;options=A.options,DDT=CDT,RDT=ARDT)
end

function *(A::judiPDEadjoint{CDT,ARDT},B::judiLRWF{BDDT,CDT}) where {ARDT,BDDT,CDT}
    A.n == size(B,1) || throw(judiPDEadjointException("Shape mismatch: A:$(size(A)), B: $(size(B))"))
    return adjoint(judiPDEextended(A.info,A.model,B.wavelet,A.geometry,options=A.options,DDT=CDT,RDT=ARDT))
end

#multiplications with judiExtendedSource
# *(judiPDE,judiExtendedSource)
function *(A::judiPDE{ADDT,ARDT}, B::judiExtendedSource{vDT}) where {ADDT, ARDT,vDT}
    A.n == size(B,1) || throw(judiPDEexception("Shape mismatch: A:$(size(A)), B: $(size(B))"))
    jo_check_type_match(ADDT, vDT, join(["DDT for *(judiPDE,judiExtendedSource):",A.name,typeof(A),vDT]," / "))
    F = judiPDEextended(A.info, A.model, B.wavelet, A.geometry; options=A.options,DDT=vDT, RDT=ARDT)
    V = F*B.weights
    jo_check_type_match(ARDT, eltype(V),join(["RDT from *(judiPDE, judiExtendedSource):",A.name,typeof(A), eltype(V)]," / "))
    return V
end

# *(judiPDEadjoint, judiExtendedSource)
function *(A::judiPDEadjoint{ADDT,ARDT}, B::judiExtendedSource{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(B, 1) || throw(judiPDEadjointException("Shape mismatch: A:$(size(A)), B: $(size(B))"))
    jo_check_type_match(ADDT, vDT, join(["DDT for *(judiPDEadjoint,judiExtendedSource):",A.name,typeof(A),vDT]," / "))
    F = adjoint(judiPDEextended(A.info, A.model, B.wavelet, A.geometry; options=A.options,DDT=vDT,RDT=ARDT))
    V = F*B.weights
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDE,judiExtendedSource):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(num,judiPDE)
function *(a::Number,A::judiPDE{ADDT,ARDT}) where {ADDT,ARDT}
    return judiPDE{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1 -> jo_convert(ARDT,a, false)*A.fop(v1),
        v2 -> jo_convert(ARDT,a, false)*A.fop_T(v2)
        )
end

function *(a::Number,A::judiPDEadjoint{ADDT,ARDT}) where {ADDT,ARDT}
    return judiPDEadjoint{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1 -> jo_convert(ARDT,a, false)*A.fop(v1),
        v2 -> jo_convert(ADDT,a, false)*A.fop_T(v2)
        )
end

############################################################
## overloaded Basees +(...judiPDE...), -(...judiPDE...)

# -(judiPDE)
-(A::judiPDE{DDT,RDT}) where {DDT,RDT} =
    judiPDE{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1->-A.fop(v1),
        v2->-A.fop_T(v2)
        )

-(A::judiPDEadjoint{DDT,RDT}) where {DDT,RDT} =
    judiPDEadjoint{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        v1->-A.fop(v1),
        v2->-A.fop_T(v2)
        )

############################################################
## Additional overloaded functions

# Subsample Modeling operator
function subsample(F::judiPDE{ADDT,ARDT}, srcnum) where {ADDT,ARDT}
    geometry = subsample(F.geometry,srcnum)     # Geometry of subsampled data container
    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiPDE(F.name, info, F.model, geometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

function subsample(F::judiPDEadjoint{ADDT,ARDT}, srcnum) where {ADDT,ARDT}
    geometry = subsample(F.geometry,srcnum)     # Geometry of subsampled data container
    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiPDEadjoint(F.name, info, F.model, geometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

getindex(F::judiPDE,a) = subsample(F,a)
getindex(F::judiPDEadjoint,a) = subsample(F,a)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Missing, missing and ismissing are defined in essentials.jl

show(io::IO, x::Missing) = print(io, "missing")

"""
    MissingException(msg)

Exception thrown when a [`missing`](@ref) value is encountered in a situation
where it is not supported. The error message, in the `msg` field
may provide more specific details.
"""
struct MissingException <: Exception
    msg::AbstractString
end

showerror(io::IO, ex::MissingException) =
    print(io, "MissingException: ", ex.msg)

"""
    nonmissingtype(T::Type)

If `T` is a union of types containing `Missing`, return a new type with
`Missing` removed.

# Examples
```jldoctest
julia> nonmissingtype(Union{Int64,Missing})
Int64

julia> nonmissingtype(Any)
Any
```

!!! compat "Julia 1.3"
    This function is exported as of Julia 1.3.
"""
nonmissingtype(::Type{T}) where {T} = typesplit(T, Missing)

function nonmissingtype_checked(T::Type)
    R = nonmissingtype(T)
    R >: T && error("could not compute non-missing type")
    return R
end

promote_rule(T::Type{Missing}, S::Type) = Union{S, Missing}
promote_rule(T::Type{Union{Nothing, Missing}}, S::Type) = Union{S, Nothing, Missing}
function promote_rule(T::Type{>:Union{Nothing, Missing}}, S::Type)
    R = nonnothingtype(T)
    R >: T && return Any
    T = R
    R = nonmissingtype(T)
    R >: T && return Any
    T = R
    R = promote_type(T, S)
    return Union{R, Nothing, Missing}
end
function promote_rule(T::Type{>:Missing}, S::Type)
    R = nonmissingtype(T)
    R >: T && return Any
    T = R
    R = promote_type(T, S)
    return Union{R, Missing}
end

convert(::Type{T}, x::T) where {T>:Missing} = x
convert(::Type{T}, x::T) where {T>:Union{Missing, Nothing}} = x
convert(::Type{T}, x) where {T>:Missing} = convert(nonmissingtype_checked(T), x)
convert(::Type{T}, x) where {T>:Union{Missing, Nothing}} = convert(nonmissingtype_checked(nonnothingtype_checked(T)), x)

# Comparison operators
==(::Missing, ::Missing) = missing
==(::Missing, ::Any) = missing
==(::Any, ::Missing) = missing
# To fix ambiguity
==(::Missing, ::WeakRef) = missing
==(::WeakRef, ::Missing) = missing
isequal(::Missing, ::Missing) = true
isequal(::Missing, ::Any) = false
isequal(::Any, ::Missing) = false
<(::Missing, ::Missing) = missing
<(::Missing, ::Any) = missing
<(::Any, ::Missing) = missing
isless(::Missing, ::Missing) = false
isless(::Missing, ::Any) = false
isless(::Any, ::Missing) = true
isapprox(::Missing, ::Missing; kwargs...) = missing
isapprox(::Missing, ::Any; kwargs...) = missing
isapprox(::Any, ::Missing; kwargs...) = missing

# Unary operators/functions
for f in (:(!), :(~), :(+), :(-), :(*), :(&), :(|), :(xor),
          :(zero), :(one), :(oneunit),
          :(isfinite), :(isinf), :(isodd),
          :(isinteger), :(isreal), :(isnan),
          :(iszero), :(transpose), :(adjoint), :(float), :(complex), :(conj),
          :(abs), :(abs2), :(iseven), :(ispow2),
          :(real), :(imag), :(sign), :(inv))
    @eval ($f)(::Missing) = missing
end
for f in (:(Base.zero), :(Base.one), :(Base.oneunit))
    @eval ($f)(::Type{Missing}) = missing
    @eval function $(f)(::Type{Union{T, Missing}}) where T
        T === Any && throw(MethodError($f, (Any,)))  # To prevent StackOverflowError
        $f(T)
    end
end
for f in (:(Base.float), :(Base.complex))
    @eval $f(::Type{Missing}) = Missing
    @eval function $f(::Type{Union{T, Missing}}) where T
        T === Any && throw(MethodError($f, (Any,)))  # To prevent StackOverflowError
        Union{$f(T), Missing}
    end
end

# Binary operators/functions
for f in (:(+), :(-), :(*), :(/), :(^), :(mod), :(rem))
    @eval begin
        # Scalar with missing
        ($f)(::Missing, ::Missing) = missing
        ($f)(::Missing, ::Number)  = missing
        ($f)(::Number,  ::Missing) = missing
    end
end

div(::Missing, ::Missing, r::RoundingMode) = missing
div(::Missing, ::Number, r::RoundingMode) = missing
div(::Number, ::Missing, r::RoundingMode) = missing

min(::Missing, ::Missing) = missing
min(::Missing, ::Any)     = missing
min(::Any,     ::Missing) = missing
max(::Missing, ::Missing) = missing
max(::Missing, ::Any)     = missing
max(::Any,     ::Missing) = missing

# Rounding and related functions
round(::Missing, ::RoundingMode=RoundNearest; sigdigits::Integer=0, digits::Integer=0, base::Integer=0) = missing
round(::Type{>:Missing}, ::Missing, ::RoundingMode=RoundNearest) = missing
round(::Type{T}, ::Missing, ::RoundingMode=RoundNearest) where {T} =
    throw(MissingException("cannot convert a missing value to type $T: use Union{$T, Missing} instead"))
round(::Type{T}, x::Any, r::RoundingMode=RoundNearest) where {T>:Missing} = round(nonmissingtype_checked(T), x, r)
# to fix ambiguities
round(::Type{T}, x::Rational{Tr}, r::RoundingMode=RoundNearest) where {T>:Missing,Tr} = round(nonmissingtype_checked(T), x, r)
round(::Type{T}, x::Rational{Bool}, r::RoundingMode=RoundNearest) where {T>:Missing} = round(nonmissingtype_checked(T), x, r)

# Handle ceil, floor, and trunc separately as they have no RoundingMode argument
for f in (:(ceil), :(floor), :(trunc))
    @eval begin
        ($f)(::Missing; sigdigits::Integer=0, digits::Integer=0, base::Integer=0) = missing
        ($f)(::Type{>:Missing}, ::Missing) = missing
        ($f)(::Type{T}, ::Missing) where {T} =
            throw(MissingException("cannot convert a missing value to type $T: use Union{$T, Missing} instead"))
        ($f)(::Type{T}, x::Any) where {T>:Missing} = $f(nonmissingtype_checked(T), x)
        # to fix ambiguities
        ($f)(::Type{T}, x::Rational) where {T>:Missing} = $f(nonmissingtype_checked(T), x)
    end
end

# to avoid ambiguity warnings
(^)(::Missing, ::Integer) = missing

# Bit operators
(&)(::Missing, ::Missing) = missing
(&)(a::Missing, b::Bool) = ifelse(b, missing, false)
(&)(b::Bool, a::Missing) = ifelse(b, missing, false)
(&)(::Missing, ::Integer) = missing
(&)(::Integer, ::Missing) = missing
(|)(::Missing, ::Missing) = missing
(|)(a::Missing, b::Bool) = ifelse(b, true, missing)
(|)(b::Bool, a::Missing) = ifelse(b, true, missing)
(|)(::Missing, ::Integer) = missing
(|)(::Integer, ::Missing) = missing
xor(::Missing, ::Missing) = missing
xor(a::Missing, b::Bool) = missing
xor(b::Bool, a::Missing) = missing
xor(::Missing, ::Integer) = missing
xor(::Integer, ::Missing) = missing

*(d::Missing, x::AbstractString) = missing
*(d::AbstractString, x::Missing) = missing

function float(A::AbstractArray{Union{T, Missing}}) where {T}
    U = typeof(float(zero(T)))
    convert(AbstractArray{Union{U, Missing}}, A)
end
float(A::AbstractArray{Missing}) = A

"""
    skipmissing(itr)

Return an iterator over the elements in `itr` skipping [`missing`](@ref) values.
The returned object can be indexed using indices of `itr` if the latter is indexable.
Indices corresponding to missing values are not valid: they are skipped by [`keys`](@ref)
and [`eachindex`](@ref), and a `MissingException` is thrown when trying to use them.

Use [`collect`](@ref) to obtain an `Array` containing the non-`missing` values in
`itr`. Note that even if `itr` is a multidimensional array, the result will always
be a `Vector` since it is not possible to remove missings while preserving dimensions
of the input.

# Examples
```jldoctest
julia> x = skipmissing([1, missing, 2])
skipmissing(Union{Missing, Int64}[1, missing, 2])

julia> sum(x)
3

julia> x[1]
1

julia> x[2]
ERROR: MissingException: the value at index (2,) is missing
[...]

julia> argmax(x)
3

julia> collect(keys(x))
2-element Vector{Int64}:
 1
 3

julia> collect(skipmissing([1, missing, 2]))
2-element Vector{Int64}:
 1
 2

julia> collect(skipmissing([1 missing; 2 missing]))
2-element Vector{Int64}:
 1
 2
```
"""
skipmissing(itr) = SkipMissing(itr)

struct SkipMissing{T}
    x::T
end
IteratorSize(::Type{<:SkipMissing}) = SizeUnknown()
IteratorEltype(::Type{SkipMissing{T}}) where {T} = IteratorEltype(T)
eltype(::Type{SkipMissing{T}}) where {T} = nonmissingtype(eltype(T))

function iterate(itr::SkipMissing, state...)
    y = iterate(itr.x, state...)
    y === nothing && return nothing
    item, state = y
    while item === missing
        y = iterate(itr.x, state)
        y === nothing && return nothing
        item, state = y
    end
    item, state
end

IndexStyle(::Type{<:SkipMissing{T}}) where {T} = IndexStyle(T)
eachindex(itr::SkipMissing) =
    Iterators.filter(i -> @inbounds(itr.x[i]) !== missing, eachindex(itr.x))
keys(itr::SkipMissing) =
    Iterators.filter(i -> @inbounds(itr.x[i]) !== missing, keys(itr.x))
@propagate_inbounds function getindex(itr::SkipMissing, I...)
    v = itr.x[I...]
    v === missing && throw(MissingException("the value at index $I is missing"))
    v
end

function show(io::IO, s::SkipMissing)
    print(io, "skipmissing(")
    show(io, s.x)
    print(io, ')')
end

# Optimized mapreduce implementation
# The generic method is faster when !(eltype(A) >: Missing) since it does not need
# additional loops to identify the two first non-missing values of each block
mapreduce(f, op, itr::SkipMissing{<:AbstractArray}) =
    _mapreduce(f, op, IndexStyle(itr.x), eltype(itr.x) >: Missing ? itr : itr.x)

function _mapreduce(f, op, ::IndexLinear, itr::SkipMissing{<:AbstractArray})
    A = itr.x
    local ai
    inds = LinearIndices(A)
    i = first(inds)
    ilast = last(inds)
    while i <= ilast
        @inbounds ai = A[i]
        ai === missing || break
        i += 1
    end
    i > ilast && return mapreduce_empty(f, op, eltype(itr))
    a1::eltype(itr) = ai
    i += 1
    while i <= ilast
        @inbounds ai = A[i]
        ai === missing || break
        i += 1
    end
    i > ilast && return mapreduce_first(f, op, a1)
    # We know A contains at least two non-missing entries: the result cannot be nothing
    something(mapreduce_impl(f, op, itr, first(inds), last(inds)))
end

_mapreduce(f, op, ::IndexCartesian, itr::SkipMissing) = mapfoldl(f, op, itr)

mapreduce_impl(f, op, A::SkipMissing, ifirst::Integer, ilast::Integer) =
    mapreduce_impl(f, op, A, ifirst, ilast, pairwise_blocksize(f, op))

# Returns nothing when the input contains only missing values, and Some(x) otherwise
@noinline function mapreduce_impl(f, op, itr::SkipMissing{<:AbstractArray},
                                  ifirst::Integer, ilast::Integer, blksize::Int)
    A = itr.x
    if ifirst == ilast
        @inbounds a1 = A[ifirst]
        if a1 === missing
            return nothing
        else
            return Some(mapreduce_first(f, op, a1))
        end
    elseif ifirst + blksize > ilast
        # sequential portion
        local ai
        i = ifirst
        while i <= ilast
            @inbounds ai = A[i]
            ai === missing || break
            i += 1
        end
        i > ilast && return nothing
        a1 = ai::eltype(itr)
        i += 1
        while i <= ilast
            @inbounds ai = A[i]
            ai === missing || break
            i += 1
        end
        i > ilast && return Some(mapreduce_first(f, op, a1))
        a2 = ai::eltype(itr)
        i += 1
        v = op(f(a1), f(a2))
        @simd for i = i:ilast
            @inbounds ai = A[i]
            if ai !== missing
                v = op(v, f(ai))
            end
        end
        return Some(v)
    else
        # pairwise portion
        imid = (ifirst + ilast) >> 1
        v1 = mapreduce_impl(f, op, itr, ifirst, imid, blksize)
        v2 = mapreduce_impl(f, op, itr, imid+1, ilast, blksize)
        if v1 === nothing && v2 === nothing
            return nothing
        elseif v1 === nothing
            return v2
        elseif v2 === nothing
            return v1
        else
            return Some(op(something(v1), something(v2)))
        end
    end
end

"""
    filter(f, itr::SkipMissing{<:AbstractArray})

Return a vector similar to the array wrapped by the given `SkipMissing` iterator
but with all missing elements and those for which `f` returns `false` removed.

!!! compat "Julia 1.2"
    This method requires Julia 1.2 or later.

# Examples
```jldoctest
julia> x = [1 2; missing 4]
2×2 Matrix{Union{Missing, Int64}}:
 1         2
  missing  4

julia> filter(isodd, skipmissing(x))
1-element Vector{Int64}:
 1
```
"""
function filter(f, itr::SkipMissing{<:AbstractArray})
    y = similar(itr.x, eltype(itr), 0)
    for xi in itr.x
        if xi !== missing && f(xi)
            push!(y, xi)
        end
    end
    y
end

"""
    coalesce(x, y...)

Return the first value in the arguments which is not equal to [`missing`](@ref),
if any. Otherwise return `missing`.

See also [`something`](@ref).

# Examples

```jldoctest
julia> coalesce(missing, 1)
1

julia> coalesce(1, missing)
1

julia> coalesce(nothing, 1)  # returns `nothing`

julia> coalesce(missing, missing)
missing
```
"""
function coalesce end

coalesce() = missing
coalesce(x::Missing, y...) = coalesce(y...)
coalesce(x::Any, y...) = x

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
using WAVEFORM
using JOLI
using Test

# Run Tests
include("test_operators.jl")
include("test_krylov.jl")
include("test_utils.jl")
include("modeling.jl")
include("test_3D.jl")


text:
# indicator of a halfspace

"""
  IndHalfspace(a::Array{Real}, b::Real)

Returns the function `g = ind{x : <a,x> ⩽ b}`.
"""

immutable IndHalfspace{T <: AbstractArray, R <: Real} <: IndicatorConvex
  a::T
  b::R
  function IndHalfspace(a::T, b::R)
    norma = vecnorm(a)
    new(a/norma, b/norma)
  end
end

IndHalfspace{T <: AbstractArray, R <: Real}(a::T, b::R) = IndHalfspace{T, R}(a, b)

function (f::IndHalfspace){T <: Real}(x::AbstractArray{T})
  s = vecdot(f.a,x)-f.b
  if s <= 1e-14
    return 0.0
  end
  return +Inf
end

function prox!{T <: Real}(f::IndHalfspace, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  s = vecdot(f.a,x)-f.b
  if s <= 0
    y[:] = x
    return 0.0
  end
  for k in eachindex(x)
    y[k] = x[k] - s*f.a[k]
  end
  return 0.0
end

fun_name(f::IndHalfspace) = "indicator of a halfspace"
fun_dom(f::IndHalfspace) = "AbstractArray{Real}"
fun_expr(f::IndHalfspace) = "x ↦ 0 if <a,x> ⩽ b, +∞ otherwise"
fun_params(f::IndHalfspace) =
  string( "a = ", typeof(f.a), " of size ", size(f.a), ", ",
          "b = $(f.b)")

function prox_naive{T <: Real}(f::IndHalfspace, x::AbstractArray{T}, gamma::Real=1.0)
  s = vecdot(f.a,x)-f.b
  if s <= 0
    return x, 0.0
  end
  return x - s*f.a, 0.0
end


text:
# Nonlinear functions, in order to encourage inlining, these sadly are not an array of lambdas
function Round0(b,c,d)
    return UInt32((b & c) | (~b & d))
end

function Round1And3(b,c,d)
    return UInt32(b ⊻ c ⊻ d)
end

function Round2(b,c,d)
    return UInt32((b & c) | (b & d) | (c & d))
end

function transform!(context::SHA1_CTX)
    # Buffer is 16 elements long, we expand to 80
    buffer = reinterpret(eltype(context.state), context.buffer)
    for i in 1:16
        context.W[i] = bswap(buffer[i])
    end

    # First round of expansions
    for i in 17:32
        @inbounds begin
            context.W[i] = lrot(1, context.W[i-3] ⊻ context.W[i-8] ⊻ context.W[i-14] ⊻ context.W[i-16], 32)
        end
    end

    # Second round of expansions (possibly 4-way SIMD-able)
    for i in 33:80
        @inbounds begin
            context.W[i] = lrot(2, context.W[i-6] ⊻ context.W[i-16] ⊻ context.W[i-28] ⊻ context.W[i-32], 32)
        end
    end

    # Initialize registers with the previous intermediate values (our state)
    a = context.state[1]
    b = context.state[2]
    c = context.state[3]
    d = context.state[4]
    e = context.state[5]

    # Run our rounds, manually separated into the four rounds, unfortunately using an array of lambdas
    # really kills performance and causes a huge number of allocations, so we make it easy on the compiler
    for i = 1:20
        @inbounds begin
            temp = UInt32(lrot(5, a, 32) + Round0(b,c,d) + e + context.W[i] + K1[1])
            e = d
            d = c
            c = lrot(30, b, 32)
            b = a
            a = temp
        end
    end

    for i = 21:40
        @inbounds begin
            temp = UInt32(lrot(5, a, 32) + Round1And3(b,c,d) + e + context.W[i] + K1[2])
            e = d
            d = c
            c = lrot(30, b, 32)
            b = a
            a = temp
        end
    end

    for i = 41:60
        @inbounds begin
            temp = UInt32(lrot(5, a, 32) + Round2(b,c,d) + e + context.W[i] + K1[3])
            e = d
            d = c
            c = lrot(30, b, 32)
            b = a
            a = temp
        end
    end

    for i = 61:80
        @inbounds begin
            temp = UInt32(lrot(5, a, 32) + Round1And3(b,c,d) + e + context.W[i] + K1[4])
            e = d
            d = c
            c = lrot(30, b, 32)
            b = a
            a = temp
        end
    end

    context.state[1] += a
    context.state[2] += b
    context.state[3] += c
    context.state[4] += d
    context.state[5] += e
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# InputStream struct for reading SceneLang scripts

"""
    InputStream

Type wrapping the `IO` from the source code of a SceneLang script.

# Fields

- `stream::IO`: IO stream from the source code
- `location::SourceLocation`: a [`SourceLocation`](@ref) that keeps track of the reading position in the source file
- `saved_char::Union{Char, Nothing}`: stores a character from [`unread_char!`](@ref) or nothing
- `saved_location::SourceLocation`: a [`SourceLocation`](@ref) storing the previous reading position
- `saved_token::Union{Token, Nothing}`: stores an unreaded [`Token`](@ref)
- `tabulations::Int`: how many columns a `<tab>` character is worth
"""
mutable struct InputStream
    stream::IO
    location::SourceLocation
    saved_char::Union{Char, Nothing}
    saved_location::SourceLocation
    saved_token::Union{Token, Nothing}
    tabulations::Int

    """
        InputStream(stream::IO, file_name::String; tabulations::Int = 8, line_num::Int = 1, col_num::Int = 1)

    Construct an instance of [`InputStream`](@ref) with location at the beginning of the file and initialize
    `saved_char` and `saved_token` to `nothing`. The initial state of the location can be set using the `line_num` and `col_num`
    keyword arguments, but remember this will not affect the starting position of the stored stream, only how it will be represented.
    """
    function InputStream(stream::IO, file_name::String; tabulations::Int = 8, line_num::Int = 1, col_num::Int = 1)
        loc = SourceLocation(file_name=file_name, line_num = line_num, col_num = col_num)
        new(stream, loc, nothing, loc, nothing, tabulations)
    end
end

function Base.iterate(stream::InputStream, state::Int = 1)
    isa((token = read_token(stream)).value, StopToken) ? nothing : (token, state + 1)
end

Base.IteratorSize(::Type{InputStream}) = Base.SizeUnknown()

"""
    open_stream(f::Function, file_name::String; tabulations::Int = 8)

Open read-only a file named `file_name` as an [`InputStream`](@ref) and apply `f` to it.
"""
function open_stream(f::Function, file_name::String; tabulations::Int = 8)
    open(file_name, "r") do io
        InputStream(io, file_name, tabulations=tabulations) |> f
    end
end

"""
    eof(stream::InputStream)

Check if the stream has reached the end-of-file.
"""
Base.eof(stream::InputStream) = eof(stream.stream)


text:
#=

    distributions.jl

    Contains code to create alternatively parameterised distributions, including copulae.

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#

import Base: rand, minimum, maximum
import Distributions: pdf, cdf
import Statistics: mean, var, std, quantile
import StatsBase: skewness, sample

abstract type Copula end

##############################################################
## Alternative parameterisation of the gamma distribution
##############################################################
struct GammaAlt{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    ω::T
    d::ContinuousUnivariateDistribution  # Underlying (truncated) Gamma distribution
    GammaAlt{T}(μ,σ,ω,d) where {T} = new{T}(μ,σ,ω,d)
end
struct GammaAltNegative{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    ω::T
    d::ContinuousUnivariateDistribution  # Underlying (truncated) Gamma distribution
    GammaAltNegative{T}(μ,σ,ω,d) where {T} = new{T}(μ,σ,ω,d)
end

"""
    GammaAlt(μ,σ,ω)

Construct a truncated (x > 0) Gamma distribution `d` where the mean, standard deviation
and skewness of the untruncated distribution are given by μ, σ and ω.

"""
function GammaAlt(μ::T,σ::T,ω::T) where {T <: Real}
    μ > 0.0 || error("Mean must be positive.")
    σ > 0.0 || error("Standard deviation must be positive.")
    if ω < 0
        α = 4/ω^2
        θ = -σ * ω / 2
        d = Truncated(Gamma(α, θ) - α*θ - μ,-Inf,0.0)
        GammaAltNegative{T}(μ,σ,ω,d)
    else
        α = 4/ω^2
        θ = σ * ω / 2
        d = Truncated(Gamma(α, θ) - α*θ + μ,0.0,Inf)
        GammaAlt{T}(μ,σ,ω,d)
    end
end

#### Evaluation
rand(rng::AbstractRNG, d::GammaAltNegative) = -rand(rng,d.d)
pdf(d::GammaAltNegative,x::Real) = pdf(d.d,-x)
logpdf(d::GammaAltNegative,x::Real) = logpdf(d.d,-x)
cdf(d::GammaAltNegative,x::Real) = 1 .- cdf(d.d,-x)
quantile(d::GammaAltNegative,p::AbstractArray) = -quantile(d.d,1.0 .- p)
quantile(d::GammaAltNegative,p::Number) = -quantile(d.d,1.0 .- p)
minimum(d::GammaAltNegative) = 0.0
maximum(d::GammaAltNegative) = Inf

rand(rng::AbstractRNG, d::GammaAlt) = rand(rng,d.d)
pdf(d::GammaAlt,x::Real) = pdf(d.d,x)
logpdf(d::GammaAlt,x::Real) = logpdf(d.d,x)
cdf(d::GammaAlt,x::Real) = cdf(d.d,x)
quantile(d::GammaAlt,p::AbstractArray) = quantile(d.d,p)
quantile(d::GammaAlt,p::Number) = quantile(d.d,p)
minimum(d::GammaAlt) = 0.0
maximum(d::GammaAlt) = Inf

##############################################################
## Alternative parameterisation of the Log Normal distribution
##############################################################
"""
    LogNormalAlt(μ,μ₁,σ₁)

Construct a truncated (x > 0) Log-Normal distribution `d` where the mean of the
untruncated distribution is `μ`. `μ₁` and `σ₁` are the standard Log Normal parameters
in the unshifted distribution.
"""
LogNormalAlt(μ,μ₁,σ₁) = Truncated(LogNormal(μ₁,σ₁) - exp(μ₁ + σ₁^2/2) + μ,0,Inf)

##############################################################
## Gaussian copula
##############################################################
"""
    GaussianCopula(P)

Construct a Gaussian copula with correlation matrix `P`.
"""
struct GaussianCopula{N} <: Copula
    P::Matrix           # Covariance matrix
    L::LowerTriangular  # Cholesky decomposition of `P`
end
GaussianCopula(P::Matrix) = GaussianCopula{size(P,1)}(P,cholesky(P).U')
GaussianCopula(p) = (P = vec_to_cor(p); GaussianCopula{size(P,1)}(P,cholesky(P).U'))

rand(rng::AbstractRNG,C::GaussianCopula{N}) where {N} = normcdf.(C.L * randn(rng,N))
rand(rng::AbstractRNG,C::GaussianCopula{N},n::Int) where {N} = normcdf.(C.L * randn(rng,N,n))
rand(C::GaussianCopula{N},n::Int) where {N} = normcdf.(C.L * randn(N,n))
sample(C::GaussianCopula{N},n::Int=1) where {N} = normcdf.(C.L * randn(N,n))

pdf(C::GaussianCopula{N},u::Vector) where {N} = (x = norminvcdf.(u); 1 / sqrt(det(C.P)) * exp.(-0.5 * x' * (inv(C.P) - I) * x))

marginalize(C::GaussianCopula{N};dims=[1,2]) where {N} = GaussianCopula(C.P[dims,dims])

##############################################################
## t copula
##############################################################
"""
    TCopula(P,ν=1)

Construct a t-copula with correlation matrix `P` and degrees of freedom `ν`.

WARNING: CURRENTLY UNTESTED

"""
struct TCopula{N} <: Copula
    P::Matrix           # Covariance matrix
    ν::Real             # df parameter
    L::LowerTriangular  # Cholesky decomposition of `P`
end
TCopula(P::Matrix,ν::Real=1) = TCopula{size(P,1)}(P,ν,cholesky(P).U')
TCopula(p,ν::Real=1) = (P = vec_to_cor(p); TCopula{size(P,1)}(P,ν,cholesky(P).U'))

rand(rng::AbstractRNG,C::TCopula{N}) where {N} = tdistcdf.(C.ν,C.L * randn(rng,N) / sqrt.(rand(Chisq(C.ν))' / C.ν))
rand(rng::AbstractRNG,C::TCopula{N},n::Int) where {N} = tdistcdf.(C.ν,C.L * randn(rng,N,n) ./ sqrt.(rand(Chisq(C.ν),n)' / C.ν))
rand(C::TCopula{N},n::Int) where {N} = tdistcdf.(C.ν,C.L * randn(N,n) ./ sqrt.(rand(Chisq(C.ν),n)' / C.ν))
sample(C::TCopula{N},n::Int=1) where {N} = tdistcdf.(C.ν,C.L * randn(N,n) ./ sqrt.(rand(Chisq(C.ν),n)' / C.ν))

Γ = gamma
function pdf(C::TCopula{N},u::Vector) where {N}
    x = tdistinvcdf.(C.ν,u);
    Γ((C.ν + N)/2) /
    (Γ(C.ν/2) * C.ν^(N/2) * π^(N/2) * sqrt(det(C.P))) *
    (1 + x' * inv(C.P) * x / C.ν)^(-(C.ν + N) / 2) /
    prod(tdistpdf.(C.ν,x))
end

##############################################################
## Copula-MVDistribution
##############################################################
"""
    MvDependent(C,v)

Construct dependent multivariate distribution where the marginals are described by
elements of `v` (a vector), and the dependence structure is described a copula `C`.
"""
struct MvDependent{N}
    C::Copula
    v::Vector
    q::Vector{Bool}
end
MvDependent(C::Copula,v::Vector;q=falses(length(v))) = MvDependent{length(v)}(C,v,q)

function sample(d::MvDependent{N},n=1) where {N}
    U = rand(d.C,n)
    for i = 1:N
        U[i,:] = d.q[i] ? quantile_interpolated(d.v[i],U[i,:]) : quantile(d.v[i],U[i,:])
    end
    return U
end
rand(d::MvDependent,n::Int=1) = sample(d,n)

marginalize(d::MvDependent;dims=[1,2]) = MvDependent(
    marginalize(d.C;dims),
    d.v[dims],
    q=d.q[dims]
)

pdf(d::MvDependent,x::Vector) = pdf(d.C,max.(0.0,cdf.(d.v,x))) * prod(pdf.(d.v,x))

##############################################################
## Helpful functions
##############################################################
"""
    vec_to_cor(p)

Convert number of vector `p` to a correlation matrix.

Example:

    vec_to_cor([0.1,0.2,0.3])
    3×3 Matrix{Float64}:
     1.0  0.1  0.2
     0.1  1.0  0.3
     0.2  0.3  1.0
"""
function vec_to_cor(p::Vector)
    n = 0.5 * (1 + sqrt(1 + 8length(p)))
    isinteger(n) || error("Check length of input.")
    n = Int(n)
    P = zeros(n,n)
    P[(triu(ones(n,n)) - I) .== 1.0] = p
    P += P' + I
end
vec_to_cor(p::Real) = vec_to_cor([p])

"""
    quantile_interpolated(d,u)

Approximate quantile function using interpolation. Output is comparable to
`quantile(d,u)`, with a considerable speed-up if `length(u)` is large and
the quantile function of `d` is expensive.
"""
function quantile_interpolated(d::Distribution,u::Vector;pts::Int=100)
    length(u) > 2 || "`u` must contain more than 2 elements."
    x = range(quantile(d,[extrema(u)...])...,length=pts)
    y = cdf(d,x)
    LinearInterpolation(y,x,extrapolation_bc=Line()).(u)
end

Distributions.cdf(x::Vector) = invperm(sortperm(x)) / length(x) .- 1 / 2length(x)

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using Turing, ReverseDiff, Memoization
using DifferentialEquations
using Plots
using Statistics
using StatsPlots
using DataFrames
using CSV
Turing.setadbackend(:reversediff)
Turing.setrdcache(true)
Turing.turnprogress(true)
using Distributions
using LinearAlgebra

# here we are simulating the signal of two coupled oscillators
# with amplitudes A1 and A2 of kBT/k and kBT/(k+c) where c is the coupling constant
# D = kBT/gamma
# in terms of the standard parameters of a OU process:
# sigma = sqrt(2*D)
# Theta = D/A

function corr_osc(c)
    # function that returns two time series that are correlated through a coupling coefficient c
    μ = 0.0 # mean is zero
    σ = sqrt(2) # D=1
    Θ1 = 1.0
    Θ2 = 1.0+abs(c)

    W1 = OrnsteinUhlenbeckProcess(Θ1,μ,σ,0.0,1.0)
    W2 = OrnsteinUhlenbeckProcess(Θ2,μ,σ,0.0,1.0)
    prob1 = NoiseProblem(W1,(0.0,100.0))
    prob2 = NoiseProblem(W2,(0.0,100.0))
    sol1 = solve(prob1;dt=0.1)
    sol2 = solve(prob2;dt=0.1)

    # creating the two correlated
    x1 = (sol1.u .+ sol2.u)/2
    if c>0
        x2 = (sol1.u .- sol2.u)/2
    else
        x2 = (sol2.u .- sol1.u)/2
    end
    return x1,x2
end

# Ornstein-Uhlenbeck process
@model ou(rn,T,delta_t) = begin
    ampl ~ Uniform(0.0,5.0)
    b ~ beta(5.0,1.0)

    rn[1] ~ Normal(0,sqrt(ampl))

    for i=2:T
        rn[i] ~ Normal(rn[i-1]*b,sqrt(ampl*(1-b^2)))
    end
end

# Ornstein-Uhlenbeck process with added Gaussian noise
@model oupn(rn,T,delta_t,::Type{R}=Vector{Float64}) where {R} = begin
    ampl ~ Uniform(0.0,5.0)
    b ~ Beta(5.0,1.0)
    noise_ampl ~ Uniform(0.0,1)

    b = exp(-delta_t/tau)
    r = R(undef, T)

    r[1] ~ Normal(0,sqrt(ampl))

    for i=2:T
        r[i] ~ Normal(r[i-1]*b,sqrt(ampl*(1-b^2)))
    end
    rn ~ MvNormal(r,sqrt(noise_ampl))
end

@model ou_corr(rn1,rn2,T,delta_t) = begin
    ampl1 ~ Uniform(0.0,5.0)
    ampl2 ~ Uniform(0.0,5.0)
    d ~ Uniform(0.0,5.0)
    b1 = exp(-delta_t*d/ampl1)
    b2 = exp(-delta_t*d/ampl2)

    rn1[1] ~ Normal(0,sqrt(ampl1))
    rn2[1] ~ Normal(0,sqrt(ampl2))

    for i=2:T
        rn1[i] ~ Normal(rn1[i-1]*b1,sqrt(ampl1*(1-b1^2)))
        rn2[i] ~ Normal(rn2[i-1]*b2,sqrt(ampl2*(1-b2^2)))
    end
end

# here we would like to run MCMC on several artificial data sets with varying correlation coefficients
# we should test rho from -0.9 to 0.9
rho_list = -0.9:0.1:0.9
results = DataFrame(rho = Float64[], coupling = Float64[], pearson = Float64[], c = Float64[], dc = Float64[])
tries = 5
for rho in rho_list
    coupling = 2*abs(rho)/(1-abs(rho))*sign(rho)
    println("coupling: ",coupling," rho: ",rho)
    # do rho tries times
    for i in 1:tries
        x1, x2 = corr_osc(coupling)
        pearson = Statistics.cor(x1,x2)
        println("Pearson: ",pearson)

        # lets see whether we can estimate c from the data
        y1 = x1 .+ x2
        y2 = x1 .- x2

        # Ornstein-Uhlenbeck process of two coupled oscillators
        model_fct = ou_corr(y1,y2,length(y1),0.1)
        chn = sample(model_fct, NUTS(0.65), 5000)
#        Turing.emptyrdcache()

        print(describe(chn))
        println(chn)
        ampl1 = Array(chn[:ampl1])
        ampl2 = Array(chn[:ampl2])

        a1 = std(ampl1)
        a2 = std(ampl2)

        if a1>a2
            c = (ampl1 .- ampl2)./ampl2
        else
            c = (ampl1 .- ampl2)./ampl1
        end
        c_mean = mean(c)
        dc = std(c)

        println("A1,A2: ",a1,",",a2)
        println("c estimate: ",c_mean,"std: ",dc)

        push!(results,[rho,coupling,pearson,c_mean,dc])
    end
end
println(results)
CSV.write("results2.csv",results)


text:
using SignedDistanceFunctions
using Test
using Plots
using Meshes

function Plots.plot(sdf::SignedDistanceFunctions.SDF)
    n = 100
    x = y = range(-10, 10, length=n)
    p = contourf(x, y, (x...) -> sdf(Point(x...)), levels=30, clims=(-10., 10.))
    contour!(p, x, y, (x...) -> sdf(Point(x...)), levels=[0.], c=:cyan)
end

plot(sdf(HalfSpace(Vec(1., 1.), 1.)))
plot(sdf(Sphere((1., 1.), 3.)))
plot(sdf(Ellipse((1., 1.),(0.4, 0.9), 6.)))
plot(sdf(HalfSpace(Vec(0., 1.), 1.)) ∩ sdf(HalfSpace(Vec(1., 0.), 1.)))
plot(sdf(Box((0., 0.), (4., 4.))))
plot(sdf(Box((-6., -6.), (-1., 4.))) ∪ sdf(Box((-2., -6.), (2., 0.))))
plot(sdf(Box((-6., -6.), (-1., 4.))) ∪ sdf(Box((-2., -6.), (2., 0.))) ∪ sdf(Sphere((0., 0.), 3.)))
plot(sdf(Box((-6., -6.), (-1., 4.))) ∪ sdf(Box((-2., -6.), (2., 0.))) - sdf(Sphere((0., 0.), 3.)))
plot(sdf(Sphere((1., 2.), 3.)) ∩ sdf(Sphere((0., 0.), 2.)))
plot(onion(0.2, Sphere((1.,1.),1.)))
plot(round(2., onion(0., sdf(Box((-5.,-5.),(5.,5.))))))
plot(smooth_union(sdf(Box((-6., -6.), (-1., 4.))), sdf(Box((0., -6.), (2., 0.))), 3.))
plot(smooth_union(sdf(Box((-6., -6.), (-1., 4.))), sdf(Sphere((2.,0.), 2.)), 5.))


text:
using Test
using Random
using Statistics

using EnsembleKalmanProcesses.Observations

@testset "Observations" begin

    # Generate samples as vector of vectors
    sample_dim = 3  # sample dimension, i.e., number of elements per sample
    n_samples = 5   # number of samples
    samples = [vec(i * ones(sample_dim)) for i in 1:n_samples]
    data_names = ["d$(string(i))" for i in 1:sample_dim]
    obs = Obs(samples, data_names)
    @test obs.data_names == data_names
    @test obs.mean == [3.0, 3.0, 3.0]
    @test obs.obs_noise_cov == 2.5 * ones(3, 3)

    # Generate samples as vector of vectors, pass a covariance to Obs
    covar = ones(sample_dim, sample_dim)
    obs = Obs(samples, covar, data_names)
    @test obs.obs_noise_cov == covar
    covar_wrong_dims = ones(sample_dim, sample_dim - 1)
    @test_throws AssertionError Obs(samples, covar_wrong_dims, data_names)

    # Generate samples as a 2d-array (each row corresponds to 1 sample)
    samples_t = vcat([i * ones(sample_dim)' for i in 1:n_samples]...)
    samples = permutedims(samples_t, (2, 1))
    obs = Obs(samples, data_names)
    @test obs.mean == [3.0, 3.0, 3.0]
    @test obs.obs_noise_cov == 2.5 * ones(3, 3)

    # Generate samples as a 2d-array (each row corresponds to 1 sample),
    # pass a covariance to Obs
    obs = Obs(samples, covar, data_names)
    @test obs.obs_noise_cov == covar
    @test_throws AssertionError Obs(samples, covar_wrong_dims, data_names)
    @test_throws AssertionError Obs(samples_t, covar, data_names) #as default is data are columns

    ## Edge cases:
    # Generating a single sample (a column vector)
    sample = reshape([1.0, 2.0, 3.0], 3, 1)
    data_name = "d1"
    obs = Obs(sample, data_name)
    @test obs.mean == vec(sample)
    @test obs.obs_noise_cov == nothing

    #Generating a single sample (vector of vector)
    sample = [[1.0, 2.0, 3.0]]
    data_name = "d1"
    obs = Obs(sample, data_name)
    @test obs.mean == sample[1]
    @test obs.obs_noise_cov == nothing

    # Generating a single sample with cov (a column vector)
    sample = reshape([1.0, 2.0, 3.0], 3, 1)
    obs_noise_cov = ones(3, 3)
    data_name = "d1"
    obs = Obs(sample, obs_noise_cov, data_name)
    @test obs.mean == vec(sample)
    @test obs.obs_noise_cov == obs_noise_cov

    #Generating a single sample with cov (vector of vector)
    sample = [[1.0, 2.0, 3.0]]
    obs_noise_cov = ones(3, 3)
    data_name = "d1"
    obs = Obs(sample, obs_noise_cov, data_name)
    @test obs.mean == sample[1]
    @test obs.obs_noise_cov == obs_noise_cov

    # Generate 1D-samples (data are columns of the row vector) -- this should result in scalar values for the mean and obs_noise_cov
    sample = reshape([1.0, 2.0, 3.0], 1, 3)
    data_name = "d1"
    obs = Obs(sample, data_name)
    @test obs.mean == 2.0
    @test obs.obs_noise_cov == 1.0

    # Generate 1D-samples (Vector of vector)
    sample = [[1.0], [2.0], [3.0]]
    data_name = "d1"
    obs = Obs(sample, data_name)
    @test obs.mean == 2.0
    @test obs.obs_noise_cov == 1.0

    # Generate 1D-samples, but also passing a known "1x1 covarince matrix" (stored as a variance)
    # (data are columns of the row vector)
    # this should result in scalar values for the mean and obs_noise_cov
    sample = reshape([1.0, 2.0, 3.0], 1, 3)
    obs_noise_cov = reshape([3.0], 1, 1) #needs to be of type Array{,2}
    data_name = "d1"
    obs = Obs(sample, obs_noise_cov, data_name)
    @test obs.mean == 2.0
    @test obs.obs_noise_cov == 3.0

    # Generate 1D-samples, but also passing a known "1x1 covarince matrix"
    # (data are vector of vectors)
    sample = [[1.0], [2.0], [3.0]]
    obs_noise_cov = reshape([3.0], 1, 1) #needs to be of type Array{,2}
    data_name = "d1"
    obs = Obs(sample, obs_noise_cov, data_name)
    @test obs.mean == 2.0
    @test obs.obs_noise_cov == 3.0

end


text:
using Parameters
using ModelParameters

abstract type HyModNode <: NetworkNode end

"""
Simple implementation of HyMod - does not include snow melt processes (see [1]).

Adapted with kind permission from: https://github.com/jdherman/GRA-2020-SALib

# References
1. Gharari, S., Hrachowitz, M., Fenicia, F., Savenije, H.H.G., 2013.
    An approach to identify time consistent model parameters: sub-period calibration.
    Hydrology and Earth System Sciences 17, 149–161.
    https://doi.org/10.5194/hess-17-149-2013
"""
Base.@kwdef mutable struct SimpleHyModNode{P, A<:Real} <: HyModNode
    @network_node

    # parameters
    Sm_max::P = Param(250.0, bounds=(1.0, 500.0))
    B::P = Param(1.0, bounds=(0.0, 2.0))
    alpha::P = Param(0.2, bounds=(0.0, 1.0))
    Kf::P = Param(0.5, bounds=(0.1, 0.9999))
    Ks::P = Param(0.05, bounds=(0.001, 0.1))

    # stores
    Sm::Array{A} = [0.0]
    Sf1::Array{A} = [0.0]
    Sf2::Array{A} = [0.0]
    Sf3::Array{A} = [0.0]
    Ss1::Array{A} = [0.0]

    outflow::Array{A} = []
end

function SimpleHyModNode(name::String, spec::Dict)
    n = create_node(SimpleHyModNode, name, spec["area"])
    node_params = spec["parameters"]
    for (p_name, p_val) in node_params
        sym = Symbol(p_name)
        p = getfield(n, sym)
        setfield!(n, sym, Param(p_val, bounds=p.bounds))
    end

    return n
end

function SimpleHyModNode(name::String, area::Float64, sm_max::Float64, B::Float64,
                         alpha::Float64, Kf::Float64, Ks::Float64)
    n = create_node(SimpleHyModNode, name, area)
    update_params!(n, sm_max, B, alpha, Kf, Ks)

    # stores/output
    n.Sm = Float64[0.0]
    n.Sf1 = Float64[0.0]
    n.Sf2 = Float64[0.0]
    n.Sf3 = Float64[0.0]
    n.Ss1 = Float64[0.0]
    n.outflow = Float64[]

    return n
end

# function prep_state!(node::HyModNode, sim_length::Int64)
#     node.Sm = fill(0.0, sim_length+1)
#     node.Sf1 = fill(0.0, sim_length+1)
#     node.Sf3 = fill(0.0, sim_length+1)
#     node.Ss1 = fill(0.0, sim_length+1)

#     node.outflow = fill(0.0, sim_length)
# end

"""
    run_node!(node::HyModNode, climate::Climate;
              inflow=nothing, extraction=nothing, exchange=nothing)

Run given HyMod node for entire simulation period.
"""
function run_node!(node::HyModNode, climate::Climate;
                   inflow=nothing, extraction=nothing, exchange=nothing)
    timesteps = sim_length(climate)
    # prep_state!(node, timesteps)

    for ts in 1:timesteps
        run_timestep!(node, climate, ts;
                      inflow=inflow, extraction=extraction, exchange=exchange)
    end
end

"""
    run_timestep!(node::HyModNode, climate::Climate, timestep::Int,
                  inflow::Float64, extraction::Float64, exchange::Float64)

Run given HyMod node for a time step.
"""
function run_timestep!(node::SimpleHyModNode, climate::Climate, timestep::Int;
                       inflow=nothing, extraction=nothing, exchange=nothing)::Float64
    ts = timestep
    P, ET = climate_values(node, climate, ts)
    return run_timestep!(node, P, ET, ts; inflow=inflow, extraction=extraction, exchange=exchange)
end

function run_timestep!(node::SimpleHyModNode, rain, et, ts; inflow=nothing, extraction=nothing, exchange=nothing)
    in_flow = timestep_value(ts, node.name, "inflow", inflow)
    ext = timestep_value(ts, node.name, "extraction", extraction)
    flux = timestep_value(ts, node.name, "exchange", exchange)

    Sm = node.Sm[ts]
    Sf1 = node.Sf1[ts]
    Sf2 = node.Sf2[ts]
    Sf3 = node.Sf3[ts]
    Ss1 = node.Ss1[ts]

    return run_hymod!(node, rain, et, Sm, Sf1, Sf2, Sf3, Ss1, in_flow, ext, flux)
end

function run_hymod!(node::SimpleHyModNode, P, PET, Sm, Sf1, Sf2, Sf3, Ss1, inflow, ext, flux)
    Sm_max = node.Sm_max
    B = node.B
    alpha = node.alpha
    Kf = node.Kf
    Ks = node.Ks

    # Calculate fluxes
    Peff = P*(1 - max(1.0 - Sm/Sm_max, 0.0)^B)  # PDM model Moore 1985
    ETa = min(PET*(Sm/Sm_max), Sm)

    Qf1 = Kf*Sf1
    Qf2 = Kf*Sf2
    Qf3 = Kf*Sf3
    Qs1 = Ks*Ss1

    # update state variables
    Sm_t1 = max(Sm + P - Peff - ETa, 0.0)
    Sf1_t1 = Sf1 + alpha*Peff - Qf1
    Sf2_t1 = Sf2 + Qf1 - Qf2
    Sf3_t1 = Sf3 + Qf2 - Qf3
    Ss1_t1 = Ss1 + (1-alpha)*Peff - Qs1

    Qtmp = (Qs1 + Qf3) * node.area
    Q_t1 = inflow + Qtmp - ext + flux
    update_states(node, Sm_t1, Sf1_t1, Sf2_t1, Sf3_t1, Ss1_t1, Q_t1)

    return Q_t1
end

function update_params!(node::HyModNode, Sm_max, B, alpha, Kf, Ks)
    node.Sm_max = Param(Sm_max, bounds=node.Sm_max.bounds::Tuple)
    node.B = Param(B, bounds=node.B.bounds::Tuple)
    node.alpha = Param(alpha, bounds=node.alpha.bounds::Tuple)
    node.Kf = Param(Kf, bounds=node.Kf.bounds::Tuple)
    node.Ks = Param(Ks, bounds=node.Ks.bounds::Tuple)
end

function update_states(node::HyModNode, Sm, Sf1, Sf2, Sf3, Ss1, Q)
    append!(node.Sm, Sm)
    append!(node.Sf1, Sf1)
    append!(node.Sf2, Sf2)
    append!(node.Sf3, Sf3)
    append!(node.Ss1, Ss1)
    append!(node.outflow, Q)
end

function reset!(node::HyModNode)
    node.Sm = Float64[node.Sm[1]]
    node.Sf1 = Float64[node.Sf1[1]]
    node.Sf2 = Float64[node.Sf2[1]]
    node.Sf3 = Float64[node.Sf3[1]]
    node.Ss1 = Float64[node.Ss1[1]]
    node.outflow = Float64[]
end


text:
using HTTP
using Test
using Pluto
using Pluto: ServerSession, ClientSession, SessionActions, WorkspaceManager
using Pluto.Configuration
using Pluto.Configuration: notebook_path_suggestion, from_flat_kwargs, _convert_to_flags
using Pluto.WorkspaceManager: WorkspaceManager, poll

@testset "Web server" begin

@testset "Exports" begin
    @inferred Pluto.port_serversocket(Sockets.ip"0.0.0.0", nothing)

    port = 13432
    host = "localhost"
    local_url(suffix) = "http://$host:$port/$suffix"

    server_running() = HTTP.get(local_url("favicon.ico")).status == 200 && HTTP.get(local_url("edit")).status == 200

    # without notebook at startup
    options = Pluto.Configuration.from_flat_kwargs(; port, launch_browser=false, workspace_use_distributed=true, require_secret_for_access=false, require_secret_for_open_links=false)
    🍭 = Pluto.ServerSession(; options)
    server_task = @async Pluto.run(🍭)

    # FYI, you should normally use a PlutoEvent for things we do in this test instead of polling! Don't use this as an example.
    @test poll(10) do
        server_running()
    end

    @test isempty(🍭.notebooks)

    HTTP.get(local_url("sample/JavaScript.jl"); retry=false)

    # wait for the notebook to be added to the session
    @test poll(10) do
        length(🍭.notebooks) == 1
    end

    notebook = only(values(🍭.notebooks))

    # right now, the notebook was only added to the session and assigned an ID. Let's wait for it to get a process:
    @test poll(60) do
        haskey(WorkspaceManager.workspaces, notebook.notebook_id)
    end
    sleep(1)

    # Note that the notebook is running async right now! It's not finished yet. But we can already run these tests:

    fileA = download(local_url("notebookfile?id=$(notebook.notebook_id)"))
    fileB = tempname()
    write(fileB, sprint(Pluto.save_notebook, notebook))

    @test Pluto.only_versions_or_lineorder_differ(fileA, fileB)

    export_contents = read(download(local_url("notebookexport?id=$(notebook.notebook_id)")), String)

    @test occursin(string(Pluto.PLUTO_VERSION), export_contents)
    @test occursin("</html>", export_contents)

    SessionActions.shutdown(🍭, notebook; keep_in_session=false)

    schedule(server_task, InterruptException(); error=true)
    # wait for the server task to finish
    # normally this `wait` would rethrow the above InterruptException, but Pluto.run should catch for InterruptExceptions and not bubble them up.
    wait(server_task)
end

end # testset


text:
using AbstractTrees

Base.show(io::IO, p::AbstractPolicy) =
    AbstractTrees.print_tree(io, StructTree(p), maxdepth=get(io, :max_depth, 10))

is_expand(::AbstractEnv) = false


text:
#=
This file contains definitions (input output stream, registers and addresses) required to understand the data in the harp binaries.

Message protocol
https://github.com/harp-tech/protocol/blob/master/Binary%20Protocol%201.0%201.1%2020180223.pdf

Register and Addresses
https://bitbucket.org/fchampalimaud/device.behavior/src/master/Firmware/Behavior/app_ios_and_regs.h
=#

payloadtypes = Dict([
    (1 , UInt8),
    (2 , UInt16),
    (4 , UInt32),
    (8 , UInt64),
    (129 , Int8),
    (130 , Int16),
    (132 , Int32),
    (136 , Int64),
    (68  , Float32)
])

registerbits_A = Dict([
    (0, "PORT0_DO"),
    (1, "PORT1_DO"),
    (2, "PORT2_DO"),
    (3, "PORT0_12V"),
    (4, "PORT1_12V"),
    (5, "PORT2_12V"),
    (6, "LED0"),
    (7, "LED1"),
    (8, "RGB0"),
    (9, "RGB1"),
    (10,"DO0"),
    (11,"DO1"),
    (12,"DO2"),
    (13,"DO3")
])

registeradds_A = Dict([
    (46, "PORT0_DO"),
    (47, "PORT1_DO"),
    (48, "PORT2_DO"),
    (49, "PORT0_12V"),
    (50, "PORT1_12V"),
    (51, "PORT2_12V"),
    (52, "LED0"),
    (53, "LED1"),
    (54, "RGB0"),
    (55, "RGB1"),
    (56,"DO0"),
    (57,"DO1"),
    (58,"DO2"),
    (59,"DO3")
])

registeradds_PWM = Dict([
    (60, "FREQ_PWM_DO0"),
    (61, "FREQ_PWM_DO1"),
    (62, "FREQ_PWM_DO2"),
    (63, "FREQ_PWM_DO3"),
    (64, "DCYC_PWM_DO0"),
    (65, "DCYC_PWM_DO0"),
    (66, "DCYC_PWM_DO0"),
    (67, "DCYC_PWM_DO0"),
    (68, "PWM_START"),
    (69, "PWM_STOP")
])

max_stars_repo_name: Peter-Vincent/Harpist

max_issues_repo_name: Peter-Vincent/Harpist

max_forks_repo_name: Peter-Vincent/Harpist


text:
function mystars(x)
    if x < 0.001
        y = "***"
    elseif x < 0.05
        y = "**"
    elseif x < 0.10
        y =  "*"
    else
        y = ""
    end
end

# mystars.(table[:,2])

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
function nextslice(ind::Tuple{Int64, Int64, Int64})
    fromobs, toobs, slicelength = ind
    sliceend = min(toobs, fromobs + slicelength - 1)
    if fromobs <= sliceend
        ((fromobs, sliceend)::Union{Tuple{Int64, Int64}, Nothing}), (sliceend + 1, toobs, slicelength)
    else
        nothing, ind
    end
end

function nextdatachunk(state::Tuple{Vector{T}, IOStream, Integer, Integer, Integer}) where {T}
    buffer, iostream, fromobs, toobs, slicelength = state
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if eof(iostream) || slicelength == 0
        close(iostream)
        nothing, state
    else
        buffer = resize!(buffer, slicelength)
        read!(iostream, buffer)
        (buffer::Union{Vector{T}, Nothing}), (buffer, iostream, fromobs + slicelength, toobs, slicelength)
    end
end

function nextline(iostream::IOStream)
    if eof(iostream)
        close(iostream)
        nothing, iostream
    else
        (readline(iostream)::Union{String, Nothing}), iostream
    end
end

function nextlines(state::Tuple{Vector{String}, IOStream, Integer})
    buffer, iostream, maxlen = state
    if eof(iostream)
        close(iostream)
        nothing, state
    else
        for i = 1:maxlen
            line = readline(iostream)
            buffer[i] = line
            if eof(iostream) && i < maxlen
                buffer = resize!(buffer, i)
                break
            end
        end
        (buffer::Union{Vector{String}, Nothing}), (buffer, iostream, maxlen)
    end
end

function verifyslicelength(fromobs::Integer, toobs::Integer, slicelength::Integer)
    toobs < fromobs || slicelength <= 0 ? 0 : min(slicelength, toobs - fromobs + 1)
end

function slice(data::BitArray{1}, fromobs::Integer, toobs::Integer, slicelength::Integer)
    if fromobs > toobs
        EmptySeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(data))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do rng
                view(data, rng[1]:rng[2])
        end
    end
end

function slice(data::Vector{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T}
    if fromobs > toobs
        EmptySeq{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(data))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do rng
                view(data, rng[1]:rng[2])
        end
    end
end

function slice(data::SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T}
    if fromobs > toobs
        EmptySeq{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(data))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do rng
                view(data, rng[1]:rng[2])
        end
    end
end

function slice(data::SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, fromobs::Integer, toobs::Integer, slicelength::Integer)
    if fromobs > toobs
        EmptySeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(data))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do rng
                view(data, rng[1]:rng[2])
        end
    end
end

function mapslice(f::Function, slices::EmptySeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}, slicelength::Integer, ::Type{S}) where {S}
    if S == Bool
        EmptySeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}()
    else
        EmptySeq{SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}}()
    end
end

function mapslice(f::Function, slices::ConsSeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}, slicelength::Integer, ::Type{S}) where {S}
    if S == Bool
        buffer = BitArray{1}(undef, slicelength)
        map(slices, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do slice
            v = view(buffer, 1:length(slice))
            v .= f.(slice)
            return v
        end
    else
        buffer = Vector{S}(undef, slicelength)
        map(slices, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}) do slice
            v = view(buffer, 1:length(slice))
            v .= f.(slice)
            return v
        end
    end
end

function mapslice(f::Function, slices::EmptySeq{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}}, slicelength::Integer, ::Type{S}) where {T} where {S}
    EmptySeq{SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}}()
end

function mapslice(f::Function, slices::ConsSeq{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}}, slicelength::Integer, ::Type{S}) where {T} where {S}
    if S == Bool
        buffer = BitArray{1}(undef, slicelength)
        map(slices, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do slice
            v = view(buffer, 1:length(slice))
            v .= f.(slice)
            return v
        end
    else
        buffer = Vector{S}(undef, slicelength)
        map(slices, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}) do slice
            v = view(buffer, 1:length(slice))
            v .= f.(slice)
            return v
        end
    end
end

function mapslice2(f::Function, slices::EmptySeq{Tuple{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{U}) where {U}
    if U == Bool
        EmptySeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}()
    else
        EmptySeq{SubArray{U,1,Array{U,1},Tuple{UnitRange{Int64}},true}}()
    end
end

function mapslice2(f::Function, slices::ConsSeq{Tuple{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{U}) where {U}
    if U == Bool
        buffer = BitArray{1}(undef, slicelength)
        map(slices, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do slice
            slice1, slice2 = slice
            v = view(buffer, 1:length(slice1))
            v .= f.(slice1, slice2)
            return v
        end
    else
        buffer = Vector{U}(undef, slicelength)
        map(slices, SubArray{U,1,Array{U,1},Tuple{UnitRange{Int64}},true}) do slice
            slice1, slice2 = slice
            v = view(buffer, 1:length(slice1))
            v .= f.(slice1, slice2)
            return v
        end
    end
end

function mapslice2(f::Function, slices::EmptySeq{Tuple{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{U}) where {T} where {S} where {U}
    EmptySeq{SubArray{U,1,Array{U,1},Tuple{UnitRange{Int64}},true}}()
end

function mapslice2(f::Function, slices::ConsSeq{Tuple{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{U}) where {T} where {S} where {U}
    if U == Bool
        buffer = BitArray{1}(undef, slicelength)
        map(slices, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do slice
            slice1, slice2 = slice
            v = view(buffer, 1:length(slice1))
            v .= f.(slice1, slice2)
            return v
        end
    else
        buffer = Vector{U}(undef, slicelength)
        map(slices, SubArray{U,1,Array{U,1},Tuple{UnitRange{Int64}},true}) do slice
            slice1, slice2 = slice
            v = view(buffer, 1:length(slice1))
            v .= f.(slice1, slice2)
            return v
        end
    end
end

function mapslice3(f::Function, slices::EmptySeq{Tuple{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer)
    EmptySeq{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}()
end

function mapslice3(f::Function, slices::ConsSeq{Tuple{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer)
    buffer = BitArray{1}(undef, slicelength)
    map(slices, SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}) do slice
        slice1, slice2, slice3 = slice
        v = view(buffer, 1:length(slice1))
        v .= f.(slice1, slice2, slice3)
        return v
    end
end

function mapslice3(f::Function, slices::EmptySeq{Tuple{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{V}) where {T} where {S} where {V}
    EmptySeq{SubArray{V,1,Array{V,1},Tuple{UnitRange{Int64}},true}}()
end

function mapslice3(f::Function, slices::ConsSeq{Tuple{SubArray{Bool,1,BitArray{1},Tuple{UnitRange{Int64}},true}, SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{V}) where {T} where {S} where {V}
    buffer = Vector{V}(undef, slicelength)
    map(slices, SubArray{V,1,Array{V,1},Tuple{UnitRange{Int64}},true}) do slice
        slice1, slice2, slice3 = slice
        v = view(buffer, 1:length(slice1))
        v .= f.(slice1, slice2, slice3)
        return v
    end
end

function mapslice3(f::Function, slices::EmptySeq{Tuple{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}, SubArray{U,1,Array{U,1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{V}) where {T} where {S} where {U} where {V}
    EmptySeq{SubArray{V,1,Array{V,1},Tuple{UnitRange{Int64}},true}}()
end

function mapslice3(f::Function, slices::ConsSeq{Tuple{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}, SubArray{S,1,Array{S,1},Tuple{UnitRange{Int64}},true}, SubArray{U,1,Array{U,1},Tuple{UnitRange{Int64}},true}}}, slicelength::Integer, ::Type{V}) where {T} where {S} where {U} where {V}
    buffer = Vector{V}(undef, slicelength)
    map(slices, SubArray{V,1,Array{V,1},Tuple{UnitRange{Int64}},true}) do slice
        slice1, slice2, slice3 = slice
        v = view(buffer, 1:length(slice1))
        v .= f.(slice1, slice2, slice3)
        return v
    end
end

function getnfolds(n::Integer, stratified::Bool, len::Integer)
    n = UInt8(n)
    #len = length(label)
    perm = randperm(len)
    res = Vector{UInt8}(undef, len)
    if stratified
        #stats = getstats(label)

    else
        space = map((x -> Int64(floor(x))), LinSpace(1, len, n + 1))
        obs = 1
        for i in 1:n
            from = i == 1 ? space[i] : space[i] + 1
            to = space[i + 1]
            for j in from:to
                res[perm[obs]] = i
                obs += 1
            end
        end
    end
    res
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
module ForwardDiff

using DiffRules, DiffResults
using DiffResults: DiffResult, MutableDiffResult, ImmutableDiffResult
using StaticArrays
if VERSION >= v"1.6"
    using Preferences
end
using Random
using LinearAlgebra
import SIMD: Vec

import Printf
import NaNMath
import SpecialFunctions
import LogExpFunctions
import CommonSubexpressions

const SIMDFloat = Union{Float64, Float32}
const SIMDInt = Union{
                       Int128, Int64, Int32, Int16, Int8,
                       UInt128, UInt64, UInt32, UInt16, UInt8,
                     }
const SIMDType = Union{SIMDFloat, SIMDInt}

include("prelude.jl")
include("partials.jl")
include("dual.jl")
include("config.jl")
include("apiutils.jl")
include("derivative.jl")
include("gradient.jl")
include("jacobian.jl")
include("hessian.jl")

export DiffResults

end # module


text:
function args_apply(s)
  @add_arg_table s begin
    "--label", "-l"
      help = "should label files be generated?"
      action = :store_true
    "--density", "-d"
      help = "should density images be generated?"
      action = :store_true
    "--density-postfix"
      help = "postfix added to saved density images"
      arg_type = String
      default = "-dens"
    "--density-prefix"
      help = "prefix added to saved density images"
      arg_type = String
      default = ""
    "--no-gpu", "-g"
      help = "prevent usage of gpu even if support is detected"
      action = :store_true
    "--quiet", "-q"
      help = "do not print output"
      action = :store_true
    "--acceptance-level", "-a"
      help = "level used to determine which density peaks are selected for annotation"
      arg_type = Float64
      default = 50.
    "--assessment-level", "-A"
      help = "additional density level for which counts are given; only affects logging"
      arg_type = Float64
      default = 30.
    "--patchsize", "-p"
      help = "maximal patchsize to use for applying the network"
      arg_type = Int
      default = 256
    "--overlap"
      help = "overlap used for dividing the image in patches"
      arg_type = Int
      default = 24
    "model"
      required = true
      arg_type = String
      help = "DCellC model (.dccm file) to be used"
    "images"
      nargs='+'
      required = true
      arg_type = String
      help = "image files that the model is applied to"
  end
end

# first command can be "synth" or "data"
function cmd_apply(args)

  # load the model
  model = modelload(args["model"])

  # check if gpu support should be disabled
  if args["no-gpu"] || gpu() < 0
    if !args["quiet"] println("# CPU mode") end
    at = Array{Float32}
  else
    if !args["quiet"] println("# GPU mode") end
    at = KnetArray{Float32}
  end

  # iterate through all images
  for imgfile in args["images"]

    img = imgload(imgfile)

    # calculate density maps
    dens = density_patched(model, img,
                           patchsize = args["patchsize"],
                           overlap = args["overlap"], at = at)

    # calculate labels and candidate levels
    aclevel = args["acceptance-level"]
    aslevel = args["assessment-level"]
    lbl  = label(dens, level = aclevel)
    clbl = label(dens, level = aslevel)

    # print out classification information
    if !args["quiet"]
      @printf("Image %s: %d (%.0f) %d (%.0f)", imgfile,
              length(lbl), aclevel, length(clbl), aslevel)
    end

    # write a label file if so desired
    if args["label"]
      lblsave(splitext(imgfile)[1], lbl)
    end

    # write a density image if so desired
    if args["density"]
      base = splitext(imgfile)[1]
      dname = args["density-prefix"] * base * args["density-postfix"]
      imgsave(dname, GreyscaleImage(dens))
    end
  end
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
module FourVectors

using Parameters
using StaticArrays
using LinearAlgebra

export FourVector
export Particle, SParticle

export psq, invmasssq, mass
export sphericalangles, boostfactor
export dot
include("structs.jl")

export Rx, Rx!
export Ry, Ry!
export Rz, Rz!
export Bz, Bz!
include("transformations.jl")

end # module