text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 78b36c68-895b-11eb-3754-2f78ee415f6c
begin
	using PlutoUI
end

# ╔═╡ 7a65ebb8-8958-11eb-3c1f-cf1daed9200e
md"""
# Adversarial search
"""

# ╔═╡ 88a95bfa-895b-11eb-2336-b3e894079c4f
md"""
## Application 1 (Ex 5.3)

Suppose two friends live in different cities on a map. Imagine that due to a contagious virus, one of the friends wants to avoid the other. Let's call the two persons $P$, the pursuer and $E$, the evader. Suppose that the players take turns moving. The game ends only when the players are on the same node; the terminal payoff to the pursuer is minus the total time taken (expressed in number of turns). The evader "wins" by never losing. An example is shown below (P moves first):
"""

# ╔═╡ 0a3f61fa-895c-11eb-3fae-f3cef63c52c4
PlutoUI.LocalResource("./img/adversarialmap.png", :width => 500, :align => "middle")

# ╔═╡ 6783207e-895a-11eb-2d87-bf236e515f97
md"""
A partial partial game tree for this map is shown in the image below. Each node is labeled by its position xy, e.g. bd corresponds with player $P$ on location b and player $E$ on location d. This is the starting position on the image above.

"""

# ╔═╡ 615c38a0-895c-11eb-26bf-2f78b04eed50
PlutoUI.LocalResource("./img/adversarialtree.png", :width => 500, :align => "middle")

# ╔═╡ 544bdd62-895c-11eb-14fa-f15a675494f8
md"""

1. Copy the game tree and mark the values of the terminal nodes (point of view of the pursuer).
2. Next to each internal node, write the strongest fact you can infer about its value (a number, one or more inequalities such as “≥ 14”, or a “?”).
3. Beneath each question mark, write the name of the node reached by that branch.
4. Explain how a bound on the value of the nodes in (3) can be derived from consideration of shortest-path lengths on the map, and derive such bounds for these nodes. Remember the cost to get to each leaf as well as the cost to solve it.
5. Now suppose that the tree as given, with the leaf bounds from (4), is evaluated from left to right. Circle those “?” nodes that would not need to be expanded further, given the bounds from part (4), and cross out those that need not be considered at all.
6. Can you claim anything in general about who wins the game on a map that is a tree? (Suppose that the tree is finite)
"""

# ╔═╡ 4dac7eea-895d-11eb-1ece-7bbfdf60aa4e
md"""
## Application 2 (Ex 5.8)
Consider the game illustrated below:

$(PlutoUI.LocalResource("./img/game2.png", :width => 300, :align => "middle"))
"""

# ╔═╡ b1e4afa4-895d-11eb-3f54-a32c30da203c
md"""
Player A moves first. The two players take turns moving, and each player must move his token to an open adjacent space in either direction. If the opponent occupies an adjacent space, then a player may jump over the opponent to the next open space if any. (For example, if A is on 3 and B is on 2, then A may move back to 1.) The game ends when one player reaches the opposite end of the board. If player A reaches space 4 first, then the value of the game to A is +1; if player B reaches space 1 first, then the value of the game to A is −1.

1. Draw the complete game tree, using the following conventions:
    * Write each state as $(s_A, s_B )$, where $s_A$ and $s_B$ denote the token locations.
    * Put each terminal state in a square box and write its game value in a circle.
    * Put loop states (states that already appear on the path to the root) in double square boxes. Since their value is unclear, annotate each with a “?” in a circle.
2. Now mark each node with its backed-up minimax value (also in a circle). Explain how you handle the “?” values and why. You should provide small reasoning.
3. Explain why the standard minimax algorithm would fail on this game tree and briefly sketch how you might fix it, drawing on your answer to (2). Does your modified algorithm give optimal decisions for all games with loops?
"""

# ╔═╡ 6a252a9e-895e-11eb-0410-29ab7456d6f9
md"""
## Application 3 - minimax algorithm
We will implement the minimax algorithm and apply it on a game of tic-tac-toe.Before going into the details of the algorithm implementation, we need to think about what we need and how we represent a game. From the course we know that we need the following:
* \\( S_0\\): The initial state, which specifies how the game is set up at the start.
* \\( \texttt{player}(s)\\): Defines which player has the move in a state.
* \\( \texttt{actions}(s)\\): Returns the set of legal moves in a state.
* \\( \texttt{result}(s,a)\\): The transition model, which defines the result of a move.
* \\( \texttt{terminal-test}(s)\\): A terminal test, which is true when the game is over and false otherwise.
* \\( \texttt{utility}(s, p)\\): A utility function (also called an objective function or payoff function), defines the final numeric value for a game that ends in terminal state s for a player p
* \\( \texttt{display}(g)\\): Display a game in way that is easy to interpret for us humans.

Below you find an implementation of tic-tac-toe. As with previous search methods, we provide a generic framework you can use to implement other games as well.
"""

# ╔═╡ 5c0a3b24-895f-11eb-33aa-8d74ee1eaacd
begin
	abstract type AbstractGame end
	abstract type AbstractGameState end

	"""
    	TicTacToeState

	Representation of a TicTacToe state.

	# Arguments
	- `turn::String`: current player ("X" or "O")
	- `utility::Int64`: +1/-1/0 (in case of a draw)
	- `board::Dict`: store game evolution (e.g. (1,3) => "X")
	- `board::Vector{Tuple{Int64,Int64}}`: holds all possible moves
	"""
	struct TicTacToeState <: AbstractGameState
		turn::String
		utility::Int64
		board::Dict
		moves::Vector{Tuple{Int64,Int64}}
	end

	"""
    	TicTacToeGame

	Implementation of the Tic-tac-toe game.

	# Arguments
	- `initial::TicTacToeState`: initial state
	- `h::Int64`: number of horizontal rows
	- `v::Int64`: number of vertical rows
	- `k::Int64`: number of aligned values required to win
	"""
	struct TicTacToeGame <: AbstractGame
		initial::TicTacToeState
		h::Int64
		v::Int64
		k::Int64
	end

	"""
		TicTacToeGame(initial::TicTacToeState)

	From an initial state generate a 3x3 game.
	"""
	function TicTacToeGame(initial::TicTacToeState)
		return TicTacToeGame(initial, 3, 3, 3)
	end

	"""
		TicTacToeGame()

	Generate a 3x3 empty TicTacToe game with player "X" starting
	"""
	function TicTacToeGame()
		init = TicTacToeState("X", 0, Dict(), collect((x, y) for x in 1:3 for y in 1:3))
		return TicTacToeGame(init, 3, 3, 3)
	end

	"""
    	actions(game, state)

	Retrieve possible actions from the current state of the game.
	"""
	function actions(game::TicTacToeGame, state::TicTacToeState)
		return state.moves
	end

	"""
    	result(game, state, move)

	Return the state of a game after applying a move.
	"""
	function result(game::TicTacToeGame, state::TicTacToeState, move::Tuple{Int64,Int64})
		# only allowed moves
		if move ∉ state.moves
			@warn "Illegal move"
			return state
		end

		# copy board in order not to modify previous state & store move
		board = copy(state.board)
		board[move] = state.turn

		# remaining moves
		moves = filter(x-> !isequal(x, move), state.moves)

		return TicTacToeState(state.turn == "X" ? "O" : "X",
   							  compute_utility(game, board, move, state.turn),
							  board,
							  moves)

	end

	"""
    	utility(game, state, player::String)

	For a game, quantify the utility of a state given a player's turn
	"""
	function utility(game::TicTacToeGame, state::TicTacToeState, player::String)
    	return player == "X" ? state.utility : -state.utility
	end

	"""
    	terminal_test(game, state)

	For a game, check if it's finished
	"""
	function terminal_test(game::TicTacToeGame, state::TicTacToeState)
		return (state.utility != 0) || (length(state.moves) == 0)
	end

	"""
    	to_move(game, state)

	Determine whose round it is for a game in a state.
	"""
	function to_move(game::TicTacToeGame, state::TicTacToeState)
    	return state.turn
	end

	"""
    	compute_utility(game, board, move, player)

	Compute the utility for a terminal state, dealing with all possible cases (rows, columns, diagonals)
	"""
	function compute_utility(game::TicTacToeGame, board::Dict,
							 move::Tuple{Int64, Int64}, player::String)
		# if three in a row is found, return matching utility value
		if k_in_row(game, board, move, player, (0, 1)) ||  # check rows
		   k_in_row(game, board, move, player, (1, 0)) ||  # check columns
		   k_in_row(game, board, move, player, (1, -1))||  # check main diagonal
		   k_in_row(game, board, move, player, (1, 1))     # check other diagonal

			return player == "X" ? 1 : -1
		else
			return 0
		end
	end

	"""
    	k_in_row(game, board, move, player, delta)

	Supporting function for `compute_utility`. Determine if a player obtains three in 	  a row after applying a move on the board. Checks in a direction determined by
	delta.
	"""
	function k_in_row(game::TicTacToeGame, board::Dict, move::Tuple{Int64, Int64}, 						  player::String, delta::Tuple{Int64, Int64})
		(delta_x, delta_y) = delta
		(x, y) = move
		n = 0
		# go from the current move in one direction (adding delta_x, delta_y)
		while get(board, (x,y), nothing) == player
			n = n + 1;
			x = x + delta_x
			y = y + delta_y
		end
		# reset to current position
		(x, y) = move
		# go from current move to the other direction (subtracting delta_x, delta_y)
		while get(board, (x,y), nothing) == player
			n = n + 1
			x = x - delta_x
			y = y - delta_y
		end
		# remove the duplicate check on the position of the move
		n = n - 1;
		return n >= game.k
	end

	"""
    	display(game, state)

	Display the TicTacToe game in a human-readable way.
	"""
	function Base.display(game::TicTacToeGame, state::TicTacToeState)
		io = IOBuffer()
		print(io, """Current game layout (X's score: $(state.utility), current player: $(state.turn)):\n\n""")
		for x in 1:game.h
			for y in 1:game.v
				print(io, get(state.board, (x, y), "."))
			end
			print(io, "\n")
    	end

		return println(stdout,String(take!(io)))
	end
end

# ╔═╡ 1e045506-8966-11eb-2c2f-f39e9d812b3b
md"""
### Illustration of a small game
"""

# ╔═╡ 298bb4d2-8966-11eb-2d9e-e15e24682397
begin
	# Initialise the game and show it's layout
	G = TicTacToeGame()
	state = G.initial
	display(G, state)
	# make random moves
	for _ = 1:9
		state = result(G, state, rand(actions(G, state)))
		display(G,state);println()
	end
end

# ╔═╡ 38c21f1c-8967-11eb-32ed-8bc5633cefd0
md"""
Note: for this illustration we do not check if a player wins during the (random) game,  we just fill the entire board to illustrate the game is working as intended.
"""

# ╔═╡ 4c136544-8967-11eb-322c-952f52cbf62f
@show state

# ╔═╡ 7711bc0c-8967-11eb-0fdf-e30e182b06b1
@assert terminal_test(G, state)

# ╔═╡ 8ae142e6-8967-11eb-32fb-7f214a71d8b5
md"""
### Minimax implementation
Now that we have evaluated that our game that is fully functional, we implement the minimax algorithm.
"""

# ╔═╡ b6ed262a-8967-11eb-3326-0398ad5931d5
begin
	"""
    	minimax_decision(state, game)

	Given all the possible actions of the current state of the game, returns the one that maximizes the utility. This is done by searching through the moves all the way to the terminal states.

	Makes us of a helper function inside (called `foo`) that determines the utility of a specific action).

	Note: argmax returns the index of the largest value.
	"""
	function minimax_decision(state::M, game::T) where {M <:AbstractGameState, T <: AbstractGame}
		# detemermine what player needs to make a move
		player = to_move(game, state);

		# supporting function for argmax
		function foo(action::Tuple{Int64,Int64};
					 relevant_game::AbstractGame=game,
					 relevant_player::String=player,
				     relevant_state::M=state)
			return minimax_min_value(relevant_game,
									 relevant_player,
									 result(relevant_game, relevant_state, action));
		end

		bestactionindex = argmax([foo(action) for action in actions(game, state)])

		return actions(game, state)[bestactionindex]
	end;

	"""
		minimax_min_value(game, player, state)

	return the minimal utility value of all possible actions for a game in a state
	"""
	function minimax_min_value(game::T, player::String, state::M) where {T <: AbstractGame, M <: AbstractGameState}
		if (terminal_test(game, state))
			return utility(game, state, player);
		end

		outcomes= [minimax_max_value(game, player, result(game, state, action)) for action in actions(game, state)]

		return minimum(vcat(Inf,outcomes))
	end

	"""
    	minimax_max_value(game, player, state)

	return the maximal utility value of all possible actions for a game in a state
	"""
	function minimax_max_value(game::T, player::String, state::M) where {T <: AbstractGame, M <: AbstractGameState}
		if (terminal_test(game, state))
			return utility(game, state, player)
		end

		outcomes = [minimax_min_value(game, player, result(game, state, action)) for action in actions(game, state)]

		return maximum(vcat(-Inf, outcomes))
	end
end

# ╔═╡ abb41398-8969-11eb-29d0-05fe2438b3a8
md"""
### Illustration of the workings of minimax
Below we have a game that is near completion (only two cases left). Player $X$ has the turn. It is clear that the best possible move is to choose $(3,3)$, as this will lead to a victory for $X$. We see that the minimax algorithm correctly identifies this move.
"""

# ╔═╡ c6a5b54e-8969-11eb-056e-9d4dd769f7c4
let
	G = TicTacToeGame()
	state = TicTacToeState("X",0, Dict((1,1)=>"X", (1,2)=>"O", (1,3)=>"X", (2,1)=>"O", (2,2)=>"X",(3,1)=>"X",(3,2)=>"O"), Tuple{Int64,Int64}[(2, 3), (3, 3)])
	display(G,state)
	move = minimax_decision(state, G);
	println("Player $(state.turn)'s best possible move: $(move)")
end

# ╔═╡ a65acd26-896c-11eb-21f0-3b8479d17e1c
md"""
### Actual game playing:
Try to realise the following:

* have two "intelligent" agents play against eachother
* have one "intelligent" agent play against a "dumb" agent i.e. an agent that makes random moves

What are the outcomes? Is this to be expected?
"""

# ╔═╡ 15b8bc4a-8b08-11eb-1991-95df70d5d7ee

# ╔═╡ beaa55c2-896c-11eb-30b8-7bfafa9b7441
begin
	"""
		playintelligentgame(;method::Function=minimax_decision)

	Play a TicTacToe game between two intelligent players.

	method: `minimax_decision` or `alphabeta_full_search`
	"""
	function playintelligentgame(;method::Function=minimax_decision)
		G = TicTacToeGame()
		state = G.initial
		for _ in 1:9
			move = method(state, G)
			state = result(G, state, move)
			if terminal_test(G,state)
				if state.utility == 0
					println("It's a draw!")
				else
					println("winner: $(state.turn == "X" ? "O" : "X")")
				end
				break
			end
		end
		display(G,state)
	end

	"""
		playdumbgame(;method::Function=minimax_decision)

	Play a TicTacToe game between a dumb (random moves) and an intelligent player.

	method: `minimax_decision` or `alphabeta_full_search`
	"""
	function playdumbgame(;method::Function=minimax_decision)
		G = TicTacToeGame()
		state = G.initial
		for _ in 1:9
			move = state.turn == "X" ? rand(actions(G, state)) : method(state, G)
			state = result(G, state, move)
			if terminal_test(G,state)
				if state.utility == 0
					println("It's a draw!")
				else
					println("winner: $(state.turn == "X" ? "O" : "X")")
				end
				break
			end
		end
		display(G,state)
	end
end

# ╔═╡ e6309ebc-896c-11eb-2a7e-2fdb4660393c
@time playintelligentgame()

# ╔═╡ 709fa20a-896d-11eb-3d35-0909bf0124c3
@time playdumbgame()

# ╔═╡ f859e3d4-8a23-11eb-2f36-65ccdeb9a782

# ╔═╡ ff0b3db8-8a23-11eb-261e-5b3d3a549a19
md"""
## Application 4 - α-β pruning
On the illustration below, you see the state of a two-player, zero-sum game with the values of the terminal states. Which branches won’t be explored when using α − β pruning? Suppose nodes are evaluated from the left to the right.

$(PlutoUI.LocalResource("./img/abtree.png", :width => 300, :align => "middle"))

"""

# ╔═╡ bb6e0a6a-896d-11eb-13da-adbd86f43b93
md"""
## Application 5 - $\alpha\text{-}\beta$ pruning
In the previous application you have used and studied the minimax algorithm. Below you find an implementation of the  $\alpha\text{-}\beta$ pruning algorithm (still exploring the game tree to the bottom). Use it in the same setting as before (i.e. TicTacToe) and compare its performance with the standard minimax algorithm.
"""

# ╔═╡ ee38ca70-896d-11eb-1256-97a4b1b3f321
begin
	"""
 	   alphabeta_full_search(state, game)

	Given all the possible actions of the current state of the game, returns the one that maximizes the utility. This is done by searching through the moves all the way to the terminal states, making use of α-β pruning (Fig 5.7)

	Makes us of a helper function inside (called `foo`) that determines the utility of a specific action).
	"""
	function alphabeta_full_search(state::M, game::T) where {M <:AbstractGameState, T <: AbstractGame}
		# detemermine what player needs to make a move
		player = to_move(game, state);

		# supporting function for argmax
		function foo(action::Tuple{Int64,Int64};
					 relevant_game::AbstractGame=game,
					 relevant_player::String=player,
					 relevant_state::M=state)
			 return alphabeta_full_search_min_value(relevant_game,
													relevant_player,
						result(relevant_game, relevant_state, action), -Inf, Inf);
		end

		bestactionindex = argmax([foo(action) for action in actions(game, state)])
		return actions(game, state)[bestactionindex]
	end

	function alphabeta_full_search_max_value(game::T,
											 player::String,
											 state::M,
											 α::Number, β::Number) where {T <: AbstractGame, M<:AbstractGameState}
		if (terminal_test(game, state))
			return utility(game, state, player)
		end
    	v = -Inf
    	for action in actions(game, state)
        	v = max(v, alphabeta_full_search_min_value(game, player, result(game, state, action), α, β))
        	if (v >= β)
            	return v
	        end
    	    α = max(α, v)
	    end
    	return v
	end

	function alphabeta_full_search_min_value(game::T,
											 player::String,
											 state::M,
											 α::Number, β::Number) where {T <: AbstractGame, M<:AbstractGameState}
    	if (terminal_test(game, state))
			return utility(game, state, player)
		end

		v = Inf
		for action in actions(game, state)
			v = min(v, alphabeta_full_search_max_value(game, player, result(game, state, action), α, β))
			if (v <= α)
				return v
			end
			β = min(β, v)
		end
		return v
	end

end

# ╔═╡ cd81ae40-8a22-11eb-04a3-1bb77668a99b
let
	# three remaining values
	G = TicTacToeGame()
	state = TicTacToeState("O",
						   0,
						   Dict((1,1)=>"X", (1,2)=>"O", (2,1)=>"O", (2,2)=>"X",(3,1)=>"X",(3,2)=>"O"),
						   Tuple{Int64,Int64}[(2, 3), (3, 3),(1,3)])
	display(G,state)
	@time minimax_decision(state, G)
	@time alphabeta_full_search(state,G)
end

# ╔═╡ Cell order:
# ╟─7a65ebb8-8958-11eb-3c1f-cf1daed9200e
# ╠═78b36c68-895b-11eb-3754-2f78ee415f6c
# ╟─88a95bfa-895b-11eb-2336-b3e894079c4f
# ╟─0a3f61fa-895c-11eb-3fae-f3cef63c52c4
# ╟─6783207e-895a-11eb-2d87-bf236e515f97
# ╟─615c38a0-895c-11eb-26bf-2f78b04eed50
# ╟─544bdd62-895c-11eb-14fa-f15a675494f8
# ╟─4dac7eea-895d-11eb-1ece-7bbfdf60aa4e
# ╟─b1e4afa4-895d-11eb-3f54-a32c30da203c
# ╟─6a252a9e-895e-11eb-0410-29ab7456d6f9
# ╠═5c0a3b24-895f-11eb-33aa-8d74ee1eaacd
# ╟─1e045506-8966-11eb-2c2f-f39e9d812b3b
# ╠═298bb4d2-8966-11eb-2d9e-e15e24682397
# ╟─38c21f1c-8967-11eb-32ed-8bc5633cefd0
# ╠═4c136544-8967-11eb-322c-952f52cbf62f
# ╠═7711bc0c-8967-11eb-0fdf-e30e182b06b1
# ╟─8ae142e6-8967-11eb-32fb-7f214a71d8b5
# ╠═b6ed262a-8967-11eb-3326-0398ad5931d5
# ╟─abb41398-8969-11eb-29d0-05fe2438b3a8
# ╠═c6a5b54e-8969-11eb-056e-9d4dd769f7c4
# ╟─a65acd26-896c-11eb-21f0-3b8479d17e1c
# ╠═15b8bc4a-8b08-11eb-1991-95df70d5d7ee
# ╠═beaa55c2-896c-11eb-30b8-7bfafa9b7441
# ╠═e6309ebc-896c-11eb-2a7e-2fdb4660393c
# ╠═709fa20a-896d-11eb-3d35-0909bf0124c3
# ╠═f859e3d4-8a23-11eb-2f36-65ccdeb9a782
# ╟─ff0b3db8-8a23-11eb-261e-5b3d3a549a19
# ╟─bb6e0a6a-896d-11eb-13da-adbd86f43b93
# ╠═ee38ca70-896d-11eb-1256-97a4b1b3f321
# ╠═cd81ae40-8a22-11eb-04a3-1bb77668a99b

max_stars_repo_path: PS/PS06 - adversarial search.jl

max_issues_repo_path: PS/PS06 - adversarial search.jl

max_forks_repo_path: PS/PS06 - adversarial search.jl


text:
module Datasets
import ..Cubes.Axes: axsym, axname, CubeAxis, findAxis, CategoricalAxis, RangeAxis, caxes
import ..Cubes: YAXArray, concatenatecubes, CleanMe
using ...YAXArrays: YAXArrays, YAXDefaults
using DataStructures: OrderedDict, counter
using Dates: Day,Hour,Minute,Second,Month,Year, Date, DateTime, TimeType
using IntervalSets: Interval, (..)
using CFTime: timedecode, timeencode, DateTimeNoLeap, DateTime360Day, DateTimeAllLeap
using YAXArrayBase
using YAXArrayBase: iscontdimval
using DiskArrayTools: CFDiskArray, ConcatDiskArray
using Glob: glob

struct Dataset
    cubes::OrderedDict{Symbol,YAXArray}
    axes::Dict{Symbol,CubeAxis}
end
function Dataset(;cubesnew...)
    axesall = Set{CubeAxis}()
    foreach(values(cubesnew)) do c
        ax = caxes(c)
        foreach(a->push!(axesall,a),ax)
    end
    axesall = collect(axesall)
    axnameall = axsym.(axesall)
    axesnew = Dict{Symbol,CubeAxis}(axnameall[i]=>axesall[i] for i in 1:length(axesall))
    Dataset(OrderedDict(cubesnew), axesnew)
end

function Base.show(io::IO,ds::Dataset)
    println(io,"YAXArray Dataset")
    println(io,"Dimensions: ")
    foreach(a->println(io,"   ", a),values(ds.axes))
    print(io,"Variables: ")
    foreach(i->print(io,i," "),keys(ds.cubes))
end
function Base.propertynames(x::Dataset, private::Bool=false)
    if private
        Symbol[:cubes; :axes; collect(keys(x.cubes)); collect(keys(x.axes))]
    else
        Symbol[collect(keys(x.cubes)); collect(keys(x.axes))]
    end
end
function Base.getproperty(x::Dataset,k::Symbol)
    if k === :cubes
        return getfield(x,:cubes)
    elseif k === :axes
        return getfield(x,:axes)
    else
        x[k]
    end
end
Base.getindex(x::Dataset,i::Symbol) = haskey(x.cubes, i) ? x.cubes[i] : haskey(x.axes,i) ? x.axes[i] : throw(ArgumentError("$i not found in Dataset"))
function Base.getindex(x::Dataset,i::Vector{Symbol})
  cubesnew = [j=>x.cubes[j] for j=i]
    Dataset(;cubesnew...)
end

function fuzzyfind(s::String,comp::Vector{String})
    sl = lowercase(s)
    f = findall(i->startswith(lowercase(i),sl),comp)
    if length(f) != 1
        throw(KeyError("Name $s not found"))
    else
        f[1]
    end
end
function Base.getindex(x::Dataset,i::Vector{String})
    istr = string.(keys(x.cubes))
    ids  = map(name->fuzzyfind(name,istr),i)
    syms   = map(j->Symbol(istr[j]),ids)
    cubesnew = [Symbol(i[j])=>x.cubes[syms[j]] for j=1:length(ids)]
    Dataset(;cubesnew...)
end
Base.getindex(x::Dataset,i::String)=getproperty(x,Symbol(i))
function subsetcube(x::Dataset; var=nothing, kwargs...)
    if var ===nothing
        cc = x.cubes
        Dataset(;map(ds->ds=>subsetcube(cc[ds];kwargs...),collect(keys(cc)))...)
    elseif isa(var,String) || isa(var, Symbol)
        subsetcube(getproperty(x,Symbol(var));kwargs...)
    else
        cc = x[var].cubes
        Dataset(;map(ds->ds=>subsetcube(cc[ds];kwargs...),collect(keys(cc)))...)
    end
end
function collectdims(g)
  dlist = Set{Tuple{String,Int,Int}}()
  varnames = get_varnames(g)
  foreach(varnames) do k
    d = get_var_dims(g,k)
    v = get_var_handle(g,k)
    for (len,dname) in zip(size(v),d)
      if !occursin("bnd",dname) && !occursin("bounds",dname)
        datts = get_var_attrs(g,dname)
        offs = get(datts,"_ARRAY_OFFSET",0)
        push!(dlist,(dname,offs,len))
      end
    end
  end
  outd = Dict(d[1] => (ax = toaxis(d[1],g,d[2],d[3]), offs = d[2]) for d in dlist)
  length(outd)==length(dlist) || throw(ArgumentError("All Arrays must have the same offset"))
  outd
end

function toaxis(dimname,g,offs,len)
    axname = dimname
    if !haskey(g,dimname)
      return RangeAxis(dimname, 1:len)
    end
    ar = g[dimname]
    aratts = get_var_attrs(g,dimname)
    if uppercase(axname)=="TIME" && haskey(aratts,"units")
        tsteps = timedecode(ar[:],aratts["units"],get(aratts,"calendar","standard"))
        RangeAxis(dimname,tsteps[offs+1:end])
    elseif haskey(aratts,"_ARRAYVALUES")
      vals = identity.(aratts["_ARRAYVALUES"])
      CategoricalAxis(axname,vals)
    else
      axdata = testrange(ar[offs+1:end])
      if eltype(axdata)<:AbstractString || (!issorted(axdata) && !issorted(axdata,rev=true))
        CategoricalAxis(axname,axdata)
      else
        RangeAxis(axname,axdata)
      end
    end
end
propfromattr(attr) = Dict{String,Any}(filter(i->i[1]!="_ARRAY_DIMENSIONS",attr))

"Test if data in x can be approximated by a step range"
function testrange(x)
  r = range(first(x),last(x),length=length(x))
  all(i->isapprox(i...),zip(x,r)) ? r : x
end

testrange(x::AbstractArray{<:AbstractString}) = x

_glob(x) = startswith(x,"/") ? glob(x[2:end], "/") : glob(x)

open_mfdataset(g::AbstractString; kwargs...) = open_mfdataset(_glob(g); kwargs...)
open_mfdataset(g::Vector{<:AbstractString};kwargs...) =
  merge_datasets(map(i->open_dataset(i;kwargs...), g))

function open_dataset(g; driver = :all)
  g = to_dataset(g, driver=driver)
  isempty(get_varnames(g)) && throw(ArgumentError("Group does not contain datasets."))
  dimlist = collectdims(g)
  dnames  = string.(keys(dimlist))
  varlist = filter(get_varnames(g)) do vn
    upname = uppercase(vn)
    !occursin("BNDS",upname) && !occursin("BOUNDS",upname) && !any(i->isequal(upname,uppercase(i)),dnames)
  end
  allcubes = OrderedDict{Symbol,YAXArray}()
  for vname in  varlist
    vardims = get_var_dims(g,vname)
    iax = [dimlist[vd].ax for vd in vardims]
    offs = [dimlist[vd].offs for vd in vardims]
    subs = if all(iszero,offs)
      nothing
    else
      ntuple(i->(offs[i]+1):(offs[i]+length(iax[i])),length(offs))
    end
    ar = get_var_handle(g,vname)
    att = get_var_attrs(g,vname)
    if subs !== nothing
      ar = view(ar,subs...)
    end
    allcubes[Symbol(vname)] = YAXArray(iax,ar,propfromattr(att), cleaner=CleanMe[])
  end
  sdimlist = Dict(Symbol(k)=>v.ax for (k,v) in dimlist)
  Dataset(allcubes,sdimlist)
end
Base.getindex(x::Dataset;kwargs...) = subsetcube(x;kwargs...)
YAXDataset(;kwargs...) = Dataset(YAXArrays.YAXDefaults.cubedir[];kwargs...)

function Cube(ds::Dataset; joinname="Variable")

  dl = collect(keys(ds.axes))
  dls = string.(dl)
  length(ds.cubes)==1 && return first(values(ds.cubes))
  #TODO This is an ugly workaround to merge cubes with different element types,
  # There should bde a more generic solution
  eltypes = map(eltype,values(ds.cubes))
  majtype = findmax(counter(eltypes))[2]
  newkeys = Symbol[]
  for k in keys(ds.cubes)
    c = ds.cubes[k]
    if all(axn->findAxis(axn,c)!==nothing,dls) && eltype(c)==majtype
      push!(newkeys,k)
    end
  end
  if length(newkeys)==1
    return ds.cubes[first(newkeys)]
  else
    varax = CategoricalAxis(joinname, string.(newkeys))
    return concatenatecubes([ds.cubes[k] for k in newkeys], varax)
  end
end

"""
    function createdataset(DS::Type,axlist; kwargs...)

Creates a new datacube with axes specified in `axlist`. Each axis must be a subtype
of `CubeAxis`. A new empty Zarr array will be created and can serve as a sink for
`mapCube` operations.

### Keyword arguments

* `folder=tempname()` location where the new cube is stored
* `T=Union{Float32,Missing}` data type of the target cube
* `chunksize = ntuple(i->length(axlist[i]),length(axlist))` chunk sizes of the array
* `chunkoffset = ntuple(i->0,length(axlist))` offsets of the chunks
* `persist::Bool=true` shall the disk data be garbage-collected when the cube goes out of scope?
* `overwrite::Bool=false` overwrite cube if it already exists
* `properties=Dict{String,Any}()` additional cube properties
* `fillvalue= T>:Missing ? defaultfillval(Base.nonmissingtype(T)) : nothing` fill value
* `datasetaxis="Variable"` special treatment of a categorical axis that gets written into separate zarr arrays
"""
function createdataset(DS, axlist;
  path="",
  persist = nothing,
  T=Union{Float32,Missing},
  chunksize = ntuple(i->length(axlist[i]),length(axlist)),
  chunkoffset = ntuple(i->0,length(axlist)),
  overwrite::Bool=false,
  properties=Dict{String,Any}(),
  datasetaxis = "Variable",
  kwargs...
  )
  if persist === nothing
    persist = !isempty(path)
  end
  path = getsavefolder(path, persist)
  check_overwrite(path, overwrite)
  splice_generic(x::AbstractArray,i) = [x[1:(i-1)];x[(i+1:end)]]
  splice_generic(x::Tuple,i)         = (x[1:(i-1)]...,x[(i+1:end)]...)
  myar = create_empty(DS, path)
  finalperm = nothing
  idatasetax = datasetaxis === nothing ? nothing : findAxis(datasetaxis,axlist)
  if idatasetax !== nothing
    groupaxis = axlist[idatasetax]
    axlist = splice_generic(axlist,idatasetax)
    chunksize = splice_generic(chunksize,idatasetax)
    chunkoffset = splice_generic(chunkoffset,idatasetax)
    finalperm = ((1:idatasetax-1)..., length(axlist)+1, (idatasetax:length(axlist))...)
  else
    groupaxis = nothing
  end
  foreach(axlist,chunkoffset) do ax,co
    arrayfromaxis(myar,ax,co)
  end
  attr = properties
  s = map(length,axlist) .+ chunkoffset
  if all(iszero,chunkoffset)
    subs = nothing
  else
    subs = ntuple(length(axlist)) do i
      (chunkoffset[i]+1):(length(axlist[i].values)+chunkoffset[i])
    end
  end
  if groupaxis===nothing
    cubenames = ["layer"]
  else
    cubenames = groupaxis.values
  end
  cleaner = CleanMe[]
  persist || push!(cleaner,CleanMe(path,false))
  allcubes = map(cubenames) do cn
    axnames = map(axname,axlist)
    axlengths = map(i->length(get_var_handle(myar,i)),axnames)
    v = if allow_missings(myar) || !(T>:Missing)
      add_var(myar, T, cn, axlengths, axnames, attr; chunksize = chunksize, kwargs...)
    else
      S = Base.nonmissingtype(T)
      if !haskey(attr,"missing_value")
        attr["missing_value"] = YAXArrayBase.defaultfillval(S)
      end
      v = add_var(myar, S, cn, axlengths, axnames, attr; chunksize = chunksize, kwargs...)
      CFDiskArray(v,attr)
    end
    if subs !== nothing
      v = view(v,subs...)
    end
    YAXArray(axlist,v,propfromattr(attr),cleaner=cleaner)
  end
  if groupaxis===nothing
    return allcubes[1],allcubes[1]
  else
    cube = concatenatecubes(allcubes,groupaxis)
    return permutedims(cube, finalperm), cube
  end
end

function getsavefolder(name,persist)
  if isempty(name)
    name = persist ? split(tempname(),"/")[end] : tempname()[2:end]
    joinpath(YAXDefaults.workdir[],name)
  else
    occursin("/",name) ? name : joinpath(YAXDefaults.workdir[],name)
  end
end

function check_overwrite(newfolder, overwrite)
  if isdir(newfolder) || isfile(newfolder)
    if overwrite
      rm(newfolder, recursive=true)
    else
      error("$(newfolder) already exists, please pick another name or use `overwrite=true`")
    end
  end
end

function arrayfromaxis(p,ax::CubeAxis,offs)
    data, attr = dataattfromaxis(ax,offs)
    attr["_ARRAY_OFFSET"]=offs
    za = add_var(p, data, axname(ax), (axname(ax),), attr)
    za
end

prependrange(r::AbstractRange,n) = n==0 ? r : range(first(r)-n*step(r),last(r),length=n+length(r))
function prependrange(r::AbstractVector,n)
  if n==0
    return r
  else
    step = r[2]-r[1]
    first = r[1] - step*n
    last = r[1] - step
    radd = range(first,last,length=n)
    return [radd;r]
  end
end

defaultcal(::Type{<:TimeType}) = "standard"
defaultcal(::Type{<:DateTimeNoLeap}) = "noleap"
defaultcal(::Type{<:DateTimeAllLeap}) = "allleap"
defaultcal(::Type{<:DateTime360Day}) = "360_day"

datetodatetime(vals::AbstractArray{<:Date}) = DateTime.(vals)
datetodatetime(vals) = vals

function dataattfromaxis(ax::CubeAxis,n)
    prependrange(ax.values,n), Dict{String,Any}()
end
# function dataattfromaxis(ax::CubeAxis,n)
#     prependrange(1:length(ax.values),n), Dict{String,Any}("_ARRAYVALUES"=>collect(ax.values))
# end
function dataattfromaxis(ax::CubeAxis{T},n) where T<:TimeType
    data = timeencode(datetodatetime(ax.values),"days since 1980-01-01",defaultcal(T))
    prependrange(data,n), Dict{String,Any}("units"=>"days since 1980-01-01","calendar"=>defaultcal(T))
end

#The good old Cube function:
Cube(s::String;kwargs...) = Cube(open_dataset(s);kwargs...)
function Cube(;kwargs...)
  if !isempty(YAXArrays.YAXDefaults.cubedir[])
    Cube(YAXArrays.YAXDefaults.cubedir[];kwargs...)
  else
    error("A path should be specified")
  end
end

#Defining joins of Datasets
abstract type AxisJoin end
struct AllEqual <: AxisJoin
     ax
end
struct SortedRanges <: AxisJoin
     axlist
     perm
end
blocksize(x::AllEqual) = 1
blocksize(x::SortedRanges) = length(x.axlist)
getperminds(x::AllEqual) = 1:1
getperminds(x::SortedRanges) = x.perm
wholeax(x::AllEqual) = x.ax
wholeax(x::SortedRanges) = reduce(vcat,x.axlist[x.perm])
struct NewDim <: AxisJoin
     newax
end
#Test for a range of categorical axes how to concatenate them
function analyse_axjoin_ranges(dimvallist)
    firstax = first(dimvallist)
    if all(isequal(firstax), dimvallist)
        return AllEqual(firstax)
    end
    revorder = if all(issorted,dimvallist)
        false
    elseif all(i->issorted(i,rev=true),dimvallist)
        true
    else
        error("Dimension values are not sorted")
    end
    function ltfunc(ax1,ax2)
        min1,max1 = extrema(ax1)
        min2,max2 = extrema(ax2)
        if max1 < min2
            return true
        elseif min1 > max2
            return false
        else
            error("Dimension ranges overlap")
        end
    end
    sp = sortperm(dimvallist, rev = revorder, lt = ltfunc)
    SortedRanges(dimvallist, sp)
end
using YAXArrayBase: YAXArrayBase, getdata, getattributes, yaxcreate
function create_mergedict(dimvallist)
    allmerges = Dict{Symbol,Any}()

    for (axn,dimvals) in dimvallist
        iscont = iscontdimval.(dimvals)
        if all(iscont)
            allmerges[axn] = analyse_axjoin_ranges(dimvals)
        elseif any(iscont)
            error("Mix of continous and non-continous values")
        else
            allmerges[axn] = analyse_axjoin_categorical(dimvals)
        end
    end
    allmerges
end
function merge_datasets(dslist)
    allaxnames = counter(Symbol)
    for ds in dslist, k in keys(ds.axes)
        push!(allaxnames,k)
    end
    dimvallist = Dict(ax=>map(i->i.axes[ax].values,dslist) for ax in keys(allaxnames))
    allmerges = create_mergedict(dimvallist)
    repvars = counter(Symbol)
    for ds in dslist, v in keys(ds.cubes)
        push!(repvars,v)
    end
    tomergevars = filter(i->i[2]==length(dslist), repvars)
    mergedvars=Dict{Symbol,Any}()
    for v in keys(tomergevars)
        dn = YAXArrayBase.dimnames(first(dslist)[v])
        howmerge = getindex.(Ref(allmerges), dn)
        sizeblockar = map(blocksize,howmerge)
        perminds = map(getperminds,howmerge)
        @assert length(dslist) == prod(sizeblockar)
        vcol = map(i->getdata(i[v]), dslist)
        allatts = mapreduce(i->getattributes(i[v]),merge,dslist, init=Dict{String,Any}())
        aa = [vcol[i] for (i,_) in enumerate(Iterators.product(perminds...))]
        dvals = map(wholeax, howmerge)
        mergedvars[v] = yaxcreate(YAXArray, ConcatDiskArray(aa), dn, dvals, allatts)
    end
    Dataset(;mergedvars...)
end

end

ext: jl


text:
# This file holds the definition of the report retrieval functions pertaining to
#   the SimulationReport type.

export getCountReport,
       getFluxInReport,
       getFluxOutReport,
       getFluxOutBreakdown,
       getAgeDistributionReport,
       getAgeStatsReport

"""
```
getCountReport( mpSim::ManpowerSimulation,
                timeRes::T )
    where T <: Real
```
This function retrieves a report of the number of active persons in the manpower
simulation `mpSim` for a time grid with time resolution `timeRes`. If the
current simulation time is zero, or the time resolution is negative, no report
is retrieved.

This function returns a `Tuple{Vector{Float64}, Vector{Int}}` where the first
element is the time grid of the report, and the second element is the number of
active personnel at each time step. If no report is available for whatever
reason, the function returns `nothing`.
"""
# Functions used
# --------------
# From same file:
# - safeInitialiseReport
# - isReportGenerated
# From simulationReport.jl:
# - generateCountReport
function getCountReport( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    # Try to initialise a report if needed. If it didn't work, don't continue.
    if !safeInitialiseReport( mpSim, timeRes )
        return
    end  # if !isReportInitialised( mpSim, timeRes )

    # Create the report if it doesn't exist yet.
    if !isReportGenerated( mpSim, timeRes, :count )
        generateCountReport( mpSim, timeRes )
    end  # if !isReportGenerated( mpSim, timeRes, :count )

    simReport = mpSim.simReports[ timeRes ]
    return ( simReport.timeGrid, simReport.activeCount )

end  # getCountReport( mpSim, timeRes )

"""
```
getFluxInReport( mpSim::ManpowerSimulation,
                 timeRes::T )
    where T <: Real
```
This function retrieves a report of the number of persons in the manpower
simulation `mpSim` who became active in each time interval of a time grid with
time resolution `timeRes`. If the current simulation time is zero, or the time
resolution is negative, no report is retrieved.

This function returns a `Tuple{Vector{Float64}, Vector{Int}}` where the first
element is the time grid of the report, and the second element is the number of
personnel who became active in each time interval. If no report is available for
whatever reason, the function returns `nothing`.
"""
# Functions used
# --------------
# From same file:
# - safeInitialiseReport
# - isReportGenerated
# From simulationReport.jl:
# - generateFluxInReport
function getFluxInReport( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    # Try to initialise a report if needed. If it didn't work, don't continue.
    if !safeInitialiseReport( mpSim, timeRes )
        return
    end  # if !isReportInitialised( mpSim, timeRes )

    # Create the report if it doesn't exist yet.
    if !isReportGenerated( mpSim, timeRes, :fluxIn )
        generateFluxInReport( mpSim, timeRes )
    end  # if !isReportGenerated( mpSim, timeRes, :fluxIn )

    simReport = mpSim.simReports[ timeRes ]
    return ( simReport.timeGrid[ 2:end ], simReport.fluxIn )

end  # getFluxInReport( mpSim, timeRes )

"""
```
getFluxOutReport( mpSim::ManpowerSimulation,
                  timeRes::T )
    where T <: Real
```
This function retrieves a report of the number of persons in the manpower
simulation `mpSim` who became inactive in each time interval of a time grid with
time resolution `timeRes`. If the current simulation time is zero, or the time
resolution is negative, no report is retrieved.

This function returns a `Tuple{Vector{Float64}, Vector{Int}}` where the first
element is the time grid of the report, and the second element is the number of
personnel who became inactive in each time interval. If no report is available
for whatever reason, the function returns `nothing`.
"""
# Functions used
# --------------
# From same file:
# - safeInitialiseReport
# - isReportGenerated
# From simulationReport.jl:
# - generateFluxOutReport
function getFluxOutReport( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    # Try to initialise a report if needed. If it didn't work, don't continue.
    if !safeInitialiseReport( mpSim, timeRes )
        return
    end  # if !isReportInitialised( mpSim, timeRes )

    # Create the report if it doesn't exist yet.
    if !isReportGenerated( mpSim, timeRes, :fluxOut )
        generateFluxOutReport( mpSim, timeRes )
    end  # if !isReportGenerated( mpSim, timeRes, :fluxOut )

    simReport = mpSim.simReports[ timeRes ]
    return ( simReport.timeGrid[ 2:end ], simReport.fluxOut )

end  # getFluxOutReport( mpSim, timeRes )

"""
```
getFluxOutBreakdown( mpSim::ManpowerSimulation,
                     timeRes::T )
    where T <: Real
```
This function retrieves a report of the number of persons in the manpower
simulation `mpSim` who became inactive in each time interval of a time grid with
time resolution `timeRes`, broken down by cause of inactivity. If the current
simulation time is zero, or the time resolution is negative, no report is
retrieved.

This function returns a `Tuple{Vector{Float64}, Dict{String, Vector{Int}}}`
where the first element is the time grid of the report, and the second element
is the number of personnel who became inactive in each time interval per cause
(= keys of the dictionary). If no report is available for whatever reason, the
function returns `nothing`.
"""
# Functions used
# --------------
# From same file:
# - safeInitialiseReport
# - isReportGenerated
# From simulationReport.jl:
# - generateFluxOutReport
function getFluxOutBreakdown( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    # Try to initialise a report if needed. If it didn't work, don't continue.
    if !safeInitialiseReport( mpSim, timeRes )
        return
    end  # if !isReportInitialised( mpSim, timeRes )

    # Create the report if it doesn't exist yet.
    if !isReportGenerated( mpSim, timeRes, :fluxOut )
        generateFluxOutReport( mpSim, timeRes )
    end  # if !isReportGenerated( mpSim, timeRes, :fluxOut )

    simReport = mpSim.simReports[ timeRes ]
    return ( simReport.timeGrid[ 2:end ], simReport.fluxOutBreakdown )

end  # getFluxOutBreakdown( mpSim, timeRes )

"""
```
getAgeDistributionReport( mpSim::ManpowerSimulation,
                          timeRes::T1,
                          ageRes::T2 )
    where T1 <: Real where T2 <: Real
```
This function retrieves a report on the distribution of the ages of the active
persons for the manpower simulation `mpSim` for a time grid with time resolution
`timeRes`. The age points for the age distribution form a grid with resolution
`ageRes`.

This function returns a `Tuple{Vector{Float64}, Vector{Float64}, Array{Int, 2}}`
where the first element is the time grid of the report, the second element is
the age grid, and the last element is the matrix of the age distribution at
every point in the time grid, corresponding to the rows of the matrix. If no
report is available for whatever reason, the function returns `nothing`.
"""
# Functions used
# --------------
# From same file:
# - safeInitialiseReport
# - isReportGenerated
# - getAgeResolution
# From simulationReport.jl:
# - generateAgeReports
function getAgeDistributionReport( mpSim::ManpowerSimulation, timeRes::T1,
    ageRes::T2 ) where T1 <: Real where T2 <: Real

    # Try to initialise a report if needed. If it didn't work, don't continue.
    if !safeInitialiseReport( mpSim, timeRes )
        return
    end  # if !isReportInitialised( mpSim, timeRes )

    # Do nothing if the age resolution isn't positive.
    if ageRes <= 0
        warn( "Negative age resolution entered. No report generated." )
        return
    end  # if ageRes <= 0

    # Create the report if it doesn't exist yet, and if the age resolution is
    #   correct.
    if !isReportGenerated( mpSim, timeRes, :age ) ||
        ( getAgeResolution( mpSim, timeRes ) != ageRes )
        generateAgeReports( mpSim, timeRes, ageRes )
    end  # if !isReportGenerated( mpSim, timeRes, :age ) || ...

    simReport = mpSim.simReports[ timeRes ]
    return ( simReport.timeGrid, simReport.ageGrid, simReport.ageDist )

end  # getAgeDistributionReport( mpSim, timeRes, ageRes )

"""
```
getAgeStatsReport( mpSim::ManpowerSimulation,
                   timeRes::T )
    where T <: Real
```
This function retrieves a report on the statistics of the distribution of the
ages of the active persons for the manpower simulation `mpSim` for a time grid
with time resolution `timeRes`.

This function returns a `Tuple{Vector{Float64}, Array{Int, 2}}` where the first
element is the time grid of the report, and the second element is the matrix
with the statistics with each column corresponding to one of these statistics,
in that order: mean, standard deviation, median, minimum, maximum. If no report
is available for whatever reason, the function returns `nothing`.
"""
# Functions used
# --------------
# From same file:
# - safeInitialiseReport
# - isReportGenerated
# From simulationReport.jl:
# - generateAgeReports
function getAgeStatsReport( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    # Try to initialise a report if needed. If it didn't work, don't continue.
    if !safeInitialiseReport( mpSim, timeRes )
        return
    end  # if !isReportInitialised( mpSim, timeRes )

    # Create the report if it doesn't exist yet.
    if !isReportGenerated( mpSim, timeRes, :age )
        generateAgeReports( mpSim, timeRes )
    end  # if !isReportGenerated( mpSim, timeRes, :age )

    simReport = mpSim.simReports[ timeRes ]
    return ( simReport.timeGrid, simReport.ageStats )

end  # getAgeStatsReport( mpSim, timeRes )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

"""
```
isReportInitialised( mpSim::ManpowerSimulation,
                     timeRes::T )
    where T <: Real
```
This function checks whether the manpower simulation `mpSim` has a report
initialised for time resolution `timeRes`.

This function returns a `Bool`.
"""
# Functions used: none
function isReportInitialised( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    return haskey( mpSim.simReports, timeRes )

end  # isReportInitialised( mpSim, timeRes )

"""
```
safeInitialiseReport( mpSim::ManpowerSimulation,
                      timeRes::T )
    where T <: Real
```
This function initialises a report for the manpower simulation `mpSim` for time
resolution `timeRes` if necessary, and checks if such a report has actually been
initialised.

This function returns a `Bool`.
"""
# Functions used
# --------------
# From same file:
# - isReportInitialised
# From simulationReport.jl:
# - initialiseReport
function safeInitialiseReport( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    # Check if a report has been initialised for this time resolution. If not,
    #   initialise one.
    if !isReportInitialised( mpSim, timeRes )
        initialiseReport( mpSim, timeRes )
    end  # if !isReportInitialised( mpSim, timeRes )

    return isReportInitialised( mpSim, timeRes )

end  # safeInitialiseReport( mpSim, timeRes )

"""
```
isReportGenerated( mpSim::ManpowerSimulation,
                   timeRes::T,
                   report::Symbol )
    where T <: Real
```
This function checks whether the manpower simulation `mpSim` has a report for
time resolution `timeRes` of the type `report`. Allowed types are:
* `:count`
* `:fluxIn`
* `:fluxOut`
* `:age`

This function returns a `Bool`.
"""
# Functions used
# --------------
# From same file:
# - isReportInitialised
function isReportGenerated( mpSim::ManpowerSimulation, timeRes::T,
    report::Symbol ) where T <: Real

    if report ∉ [ :count, :fluxIn, :fluxOut, :age ]
        warn( "Unknown report type. This warning should not appear." )
        return false
    end  # if report ∉ [ :count, :fluxIn, :fluxOut, :age ]

    if !isReportInitialised( mpSim, timeRes )
        return false
    end  # if !isReportInitialised( mpSim, timeRes )

    simReport = mpSim.simReports[ timeRes ]

    if ( report === :count ) && !isempty( simReport.activeCount )
        return true
    elseif ( report === :fluxIn ) && !isempty( simReport.fluxIn )
        return true
    elseif ( report === :fluxOut ) && !isempty( simReport.fluxOut )
        return true
    elseif ( report === :age ) && !isempty( simReport.ageDist ) &&
        !isempty( simReport.ageStats )
        return true
    end  # if ( report === :count ) && ...

    return false

end  # isReportGenerated( mpSim, timeRes, report )

"""
```
getAgeResolution( mpSim,
                  timeRes )
    where T <: Real
```
This function retrieves the age resolution of the simulation reports for the
manpower simulation `mpSim`, with time resolution `timeRes`.

This function returns a `Float64`. If there is no age distribution report
initialised for the given time grid, or if the age grid has only one entry, the
function returns `NaN`.
"""
# Functions used
# --------------
# From same file:
# - isReportInitialised
function getAgeResolution( mpSim::ManpowerSimulation, timeRes::T ) where T <: Real

    if !isReportInitialised( mpSim, timeRes )
        return NaN
    end  # if !isReportInitialised( mpSim, timeRes )

    simReport = mpSim.simReports[ timeRes ]

    return length( simReport.ageGrid ) < 2 ? NaN :
        simReport.ageGrid[ 2 ] - simReport.ageGrid[ 1 ]

end  # getAgeResolution( mpSim, timeRes )


text:
# startup script for solving an equation

push!(LOAD_PATH, "/users/creanj/julialib_fork/PUMI.jl")
using PumiInterface # pumi interface
using PdePumiInterface  # common mesh interface - pumi
using SummationByParts  # SBP operators
include("../equation/Equation.jl")  # equation types
include("../rk4/rk4.jl")  # timestepping
include("./euler/euler.jl")  # solver functions
include("./euler/ic.jl")  # initial conditions functions
include("./euler/adaptfuncs.jl")  # mesh adapt function
include("./euler/output.jl")  # print results to file
#=
# timestepping parameters
delta_t = 0.5
t_max = 1.00
=#

# create operator
sbp = TriSBP{Float64}()  # create linear sbp operator

# create mesh
dmg_name = ".null"
#smb_name = "tri2l.smb"
#smb_name = "tri18l.smb"
smb_name = "adapt_big.smb"
mesh = PumiMesh2(dmg_name, smb_name, 1; dofpernode=4)  #create linear mesh with 1 dof per node

# create euler equation
eqn = EulerEquation(sbp)

# create vectors to hold solution at current, previous timestep
u0 = zeros(mesh.numDof)  # solution at previous timestep
u = zeros(mesh.numDof) # solution at current timestep

# populate u0 with initial condition
ICsmoothHeaviside(mesh, sbp, eqn, u0)
saveSolutionToMesh(mesh, u0)

# create C-callible function
cfunc3 = cfunction(shockRefine2, Void, (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Void}, Ptr{TriSBP{Float64}}))

# pass the function to C
createAnisoFunc(mesh.m_ptr, cfunc3, mesh.f_ptr, sbp)

writeVtkFiles("output_pre", mesh.m_ptr)

# run mesh adaptation
runAnisoAdapt(mesh.m_ptr)

writeVtkFiles("output_post", mesh.m_ptr)

# reinitialize mesh
reinitPumiMesh2(mesh)
u =zeros(mesh.numDof)
retrieveSolutionFromMesh(mesh, u);

printSolution(mesh, u)
printCoordinates(mesh)

#println("u = \n", u)


text:
#copy&paste into REPL:

using BenchmarkTools

module M

using GeneralizedGenerated

mutable struct Ef
    e::Expr
    f1::Function
    f2::Function
    f3::Function
    f4
    function Ef()
        e = :(a + b)
        f1 = Base.eval( :( (a::Int,b::Int)->$e ) )
        f2 = Core.eval( M , :( (a::Int,b::Int)->$e ) )
        f3 = (a,b)->a+b
        f4 = mk_function( :(  (a, b) -> $e ) )
        new(e,f1,f2,f3,f4)
    end
end

function eval_global(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        global a=3
        global b=5
        try
            r += Core.eval(M, ef.e)
        catch e
        end
    end
    return r
end

function eval_let(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        val1=3
        val2=5
        ex=ef.e
        try
            r += @eval begin
                    let
                        a=$val1
                        b=$val2
                        $ex
                    end
                end
        catch e
        end
    end
    return r
end

interpolate_from_dict(ex::Expr, dict) = Expr(ex.head, interpolate_from_dict.(ex.args, Ref(dict))...)
interpolate_from_dict(ex::Symbol, dict) = get(dict, ex, ex)
interpolate_from_dict(ex::Any, dict) = ex
function eval_interpolate(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        try
            r += Base.eval(interpolate_from_dict(ef.e,Dict( :a => 3, :b => 5 )))
        catch e
        end
    end
    return r
end

function eval_invokelatest(@nospecialize(a_ef::Array{Ef,1}))
    r=0
    for ef in a_ef
        r += Base.invokelatest(ef.f1,3,5)
    end
    return r
end

function eval_applylatest(@nospecialize(a_ef::Array{Ef,1}))
    r=0
    for ef in a_ef
        r += Core._apply_latest(ef.f1,3,5)
    end
    return r
end

function eval_f1(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        try
            r += ef.f1(3,5)
        catch e
        end
    end
    return r
end

function eval_f2(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        try
            r += ef.f2(3,5)
        catch e
        end
    end
    return r
end

function eval_f3(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        try
            r += ef.f3(3,5)
        catch e
        end
    end
    return r
end

function eval_f4(a_ef::Array{Ef,1})
    r=0
    for ef in a_ef
        try
            r += ef.f4(3,5)
        catch e
        end
    end
    return r
end

function f(a,b)
    try
        a+b
    catch e
       0
    end
end

end

a_ef=[ M.Ef() for i in 1:1000 ]

#1
@btime M.eval_global($a_ef)

#2
@btime M.eval_let($a_ef)

#3
@btime M.eval_interpolate($a_ef)

#4
@btime M.eval_f1($a_ef)

#5
@btime M.eval_f2($a_ef)

#6
@btime M.eval_f3($a_ef)

#7
@btime sum( [ M.f(3,5) for i in 1:1000 ] )

#8
@btime M.eval_invokelatest($a_ef)

#8
@btime M.eval_applylatest($a_ef)

#9
@btime M.eval_f4($a_ef)


text:
const RECURSIVE_MAX = 5

function run(s)
    # Your code here
    lines = eachline(IOBuffer(s))
    rules = Dict{SubString,SubString{String}}()
    for line in lines
        if isempty(line)
            break
        end
        (idx, rule) = split(line, ": ", limit=2)
        if idx == "8"
            rule = "42+"[1:end]
        elseif idx == "11"
            rule = "42_31"[1:end]
        end
        rules[idx] = rule
    end
    memory = Dict{SubString{String},String}()
    function resolve_rule(rule::SubString{String})::String
        # memoize here
        get!(memory, rule) do
            if '+' in rule
                return "(?:$(resolve_rule(rule[1:end - 1])))+"
            elseif '_' in rule
                (left, right) = split(rule, "_", limit=2)
                left_pat = resolve_rule(left)
                right_pat = resolve_rule(right)
                tmp = ("(?:$(left_pat){$(i)}$(right_pat){$(i)})" for i in 1:RECURSIVE_MAX)
                return "(?:$(join(tmp, "|")))"
            elseif '"' in rule
                return rule[2:2]
            elseif '|' in rule
                (left, right) = split(rule, " | ", limit=2)
                return "(?:$(resolve_rule(left))|$(resolve_rule(right)))"
            else
                parts = split(rule, " ")
                return join((resolve_rule(rules[i]) for i in parts), "")
            end
        end
    end

    regex_str = resolve_rule(rules["0"])
    regex = Regex("^$(regex_str)\$")
    return sum(occursin.(regex, lines))
end

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

function main()
    run(ARGS[1])
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
module BDisposal

using JuMP, Ipopt, GLPK, AmplNLWriter

export efficiencyScores, prodIndex, dmuEfficiency, dmuEfficiencyDual
# for debug only:
# export problem, convexProblem, nonConvexProblem, dmuPass2

include("EfficiencyScores.jl")
include("ProdIndex.jl")
include("IndividualDMUProblem.jl")

end # module


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import Base.push!

"""
    KVStore(kv_type = :local)

For single machine training, there are two commonly used types:

- `local`: Copies all gradients to CPU memory and updates weights there.

- `device`: Aggregates gradients and updates weights on GPU(s).
  With this setting, the `KVStore` also attempts to use GPU peer-to-peer
  communication, potentially accelerating the communication.

For distributed training, `KVStore` also supports a number of types:

- `dist_sync`: Behaves similarly to `local` but with one major difference.
  With `dist_sync`, batch-size now means the batch size used on each machine.
  So if there are `n` machines and we use batch size ``b``,
  then `dist_sync` behaves like `local` with batch size `n * b`.

- `dist_device_sync`: Identical to `dist_sync` with the difference similar
  to `device` vs `local`.

- `dist_async`: Performs asynchronous updates.
  The weights are updated whenever gradients are received from any machine.
  No two updates happen on the same weight at the same time.
  However, the order is not guaranteed.
"""
mutable struct KVStore
  handle    :: MX_KVStoreHandle
  updater_c :: Ptr{Cvoid}
  updater   :: Function

  KVStore(hdr::MX_KVStoreHandle) = new(hdr, Ptr{Cvoid}(0))
end

function KVStore(kv_type::Symbol = :local)
  @assert kv_type ∈ (:local, :device, :dist_sync, :dist_device_sync, :dist_async)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXKVStoreCreate, (char_p, Ref{MX_handle}), dump_mx_param(kv_type), ref_hdr)
  KVStore(MX_KVStoreHandle(ref_hdr[]))
end

Base.unsafe_convert(::Type{MX_handle}, obj::KVStore) =
  Base.unsafe_convert(MX_handle, obj.handle)
Base.convert(t::Type{MX_handle}, obj::KVStore) = Base.unsafe_convert(t, obj)
Base.cconvert(t::Type{MX_handle}, obj::KVStore) = Base.unsafe_convert(t, obj)

Base.show(io::IO, kv::KVStore) =
    print(io, "mx.KVStore @ $(get_type(kv))")

function _flatten_kvlist(keys::Vector{Int}, vals::Vector{<:Vector{<:NDArray}})
  @assert length(keys) == length(vals)
  keys_flt = Int[]
  vals_flt = NDArray[]
  for (k,v) in zip(keys, vals)
    append!(keys_flt, Base.ones(Int, length(v))*k)
    append!(vals_flt, v)
  end
  return (keys_flt, vals_flt)
end

"""
    init!(kv::KVStore, key::Int, val::NDArray)
    init!(kv::KVStore, keys, vals)

Initializes a single or a sequence of key-value pairs into the store.

For each key, one must `init!` it before calling `push!` or `pull!`.
When multiple workers invoke `init!` for the same key, only
the value supplied by worker with rank `0` is used. This function returns
after data has been initialized successfully.

```jldoctest
julia> kv = KVStore(:local)
mx.KVStore @ local

julia> init!(kv, 42, mx.rand(2, 3))
```
"""
init!(kv::KVStore, key::Int, val::NDArray) = init!(kv, [key], [val])
init!(kv::KVStore, key::Int, vals::Vector{<:NDArray}) =
  init!(kv, Base.ones(Int, length(vals)) * key, vals)
init!(kv::KVStore, keys::Vector{Int}, vals::Vector{<:Vector{<:NDArray}}) =
  init!(kv, _flatten_kvlist(keys, vals)...)

function init!(kv::KVStore, keys::Vector{Int}, vals::VecOfNDArray)
  @assert length(keys) == length(vals)
  keys = Cint[keys...]
  vals = MX_handle[vals...]
  @mxcall(:MXKVStoreInit, (MX_handle, MX_uint, Ptr{Cint}, Ptr{MX_handle}),
          kv, length(keys), keys, vals)
end

"""
    push!(kv::KVStore, key,  val;  priority = 0)
    push!(kv::KVStore, key,  vals; priority = 0)
    push!(kv::KVStore, keys, vals; priority = 0)

Pushes a single or a sequence of key-value pairs into the store.

This function returns immediately after adding an operator to the engine.
The actual operation is executed asynchronously. If there are consecutive
pushes to the same key, there is no guarantee on the serialization of pushes.
The execution of a push does not guarantee that all previous pushes are
finished. There is no synchronization between workers by default.
One can use ``barrier()`` to sync all workers.

`push!` and `pull!` single `NDArray`:
```jldoctest
julia> kv = KVStore(:local)
mx.KVStore @ local

julia> x = NDArray(undef, 2, 3);

julia> init!(kv, 3, x)

julia> push!(kv, 3, mx.ones(2, 3) * 8)

julia> pull!(kv, 3, x)

julia> x
2×3 mx.NDArray{Float32,2} @ CPU0:
 8.0  8.0  8.0
 8.0  8.0  8.0
```

Aggregate values and `push!`:
```jldoctest
julia> vals = [mx.ones((2, 3), gpu(0)) * 3, mx.ones((2, 3), gpu(1)) * 4];

julia> push!(kv, 3, vals)

julia> pull!(kv, 3, x)

julia> x
2×3 mx.NDArray{Float32,2} @ CPU0:
 7.0  7.0  7.0
 7.0  7.0  7.0
```

`push!` a list of key to single device:

```jldoctest
julia> keys = [4, 5];

julia> init!(kv, keys, [NDArray(undef, 2, 3), NDArray(undef, 2, 3)])

julia> push!(kv, keys, [x, x])

julia> y, z = NDArray(undef, 2, 3), NDArray(undef, 2, 3);

julia> pull!(kv, keys, [y, z])
```
"""
push!(kv::KVStore, key::Int, val::NDArray; priority::Int = 0) =
  push!(kv, [key], [val]; priority = priority)
push!(kv::KVStore, key::Int, vals::Vector{<:NDArray}; priority::Int = 0) =
  push!(kv, Base.ones(Int, length(vals)) * key, vals; priority = priority)
push!(kv:: KVStore, keys::Vector{Int}, vals::Vector{<:Vector{<:NDArray}};
      priority::Int = 0) =
  push!(kv, _flatten_kvlist(keys, vals)...; priority = priority)

function push!(kv::KVStore, keys::Vector{Int}, vals::Vector{<:NDArray}; priority::Int = 0)
  @assert length(keys) == length(vals)
  keys = Cint[keys...]
  vals = MX_handle[vals...]
  @mxcall(:MXKVStorePush, (MX_handle, MX_uint, Ptr{Cint}, Ptr{MX_handle}, Cint),
          kv, length(keys), keys, vals, priority)
end

""" Pulls a single value or a sequence of values from the store.

This function returns immediately after adding an operator to the engine.
Subsequent attempts to read from the `out` variable will be blocked until the
pull operation completes.

`pull` is executed asynchronously after all previous `pull` calls and only
the last `push` call for the same input key(s) are finished.

The returned values are guaranteed to be the latest values in the store.

See [`pull!`](@ref) for more examples.
"""
pull!(kv::KVStore, key::Int, out::NDArray; priority::Int = 0) =
  pull!(kv, [key], [out], priority = priority)
pull!(kv::KVStore, key::Int, outs::Vector{<:NDArray}; priority::Int = 0) =
  pull!(kv, Base.ones(Int, length(outs))*key, outs; priority = priority)
pull!(kv::KVStore, keys::Vector{Int}, outs::Vector{<:Vector{<:NDArray}};
      priority::Int = 0) =
  pull!(kv, _flatten_kvlist(keys, outs)...; priority = priority)

function pull!(kv::KVStore, keys::Vector{Int}, outs::Vector{<:NDArray}; priority::Int = 0)
  @assert length(keys) == length(outs)
  keys = Cint[keys...]
  outs = MX_handle[outs...]
  @mxcall(:MXKVStorePull, (MX_handle, MX_uint, Ptr{Cint}, Ptr{MX_handle}, Cint),
          kv, length(keys), keys, outs, priority)
end

function get_type(kv::KVStore)
  type_ref = Ref{char_p}(0)
  @mxcall(:MXKVStoreGetType, (MX_handle, Ref{char_p}), kv, type_ref)
  return Symbol(unsafe_string(type_ref[]))
end

function get_num_workers(kv::KVStore)
  ref_size = Ref{Cint}(0)
  @mxcall(:MXKVStoreGetGroupSize, (MX_handle, Ref{Cint}), kv, ref_size)
  return Int(ref_size[])
end

function get_rank(kv::KVStore)
  ref_rank = Ref{Cint}(0)
  @mxcall(:MXKVStoreGetRank, (MX_handle, Ref{Cint}), kv, ref_rank)
  return Int(ref_rank[])
end

"""
    barrier(kv::KVStore)

Invokes global barrier among all worker nodes.

For example, assume there are `n` machines. We would like machine `0` to first
`init` the values and then have all the workers `pull` the initialized value.
Before pulling, we can place invoke `barrier(kv)` to guarantee that the
initialization is finished.
"""
barrier(kv::KVStore) = @mxcall(:MXKVStoreBarrier, (MX_handle,), kv)

# TODO: Currently Julia does not support closure in c-callbacks, so we are making use of the
# extra handle parameter of the API to pass the updater object around. Fix this when someday
# full closure cfunction is supported in Julia.
function _kvstore_update_wrapper(key::Cint, nd_recv::MX_handle, nd_local::MX_handle,
                                 updater::Ptr{Cvoid})
  updater_func = unsafe_pointer_to_objref(updater)
  updater_func(Int(key), NDArray(MX_NDArrayHandle(nd_recv)),
               NDArray(MX_NDArrayHandle(nd_local)))
  nothing
end

"""
    setupdater!(kv, updater)

Sets a `push!` updater into the store.

This function only changes the local store.
When running on multiple machines one must use `set_optimizer`.

```jldoctest
julia> update(key, val, orig) = mx.@inplace orig += val .* .2
update (generic function with 1 method)

julia> kv = KVStore(:local)
mx.KVStore @ local

julia> mx.setupdater!(kv, update)

julia> init!(kv, 42, mx.ones(2, 3))

julia> push!(kv, 42, mx.ones(2, 3))

julia> x = NDArray(undef, 2, 3);

julia> pull!(kv, 42, x)

julia> x
2×3 mx.NDArray{Float32,2} @ CPU0:
 1.2  1.2  1.2
 1.2  1.2  1.2
```
"""
function setupdater!(kv::KVStore, updater)
  kv.updater = updater # keep a reference to the julia object so that updater_c is kept valid
  kv.updater_c = @cfunction(_kvstore_update_wrapper, Cvoid,
                            (Cint,MX_handle,MX_handle,Ptr{Cvoid}))
  @mxcall(:MXKVStoreSetUpdater, (MX_handle, Ptr{Cvoid}, Any),
          kv, kv.updater_c, updater)
end

"""
    setoptimizer!(kv::KVStore, opt)

Registers an optimizer with the kvstore.

When using a single machine, this function updates the local optimizer.
If using multiple machines and this operation is invoked from a worker node,
it will serialized the optimizer with pickle and send it to all servers.
The function returns after all servers have been updated.

```jldoctest
julia> kv = KVStore()
mx.KVStore @ local

julia> W = mx.zeros(2, 3)  # 2×3 weight matrix
2×3 mx.NDArray{Float32,2} @ CPU0:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> init!(kv, 42, W)

julia> setoptimizer!(kv, SGD(η = .2))  # SGD with .2 as learning rate

julia> ∇W = mx.ones(2, 3)  # assume it's the gradient
2×3 mx.NDArray{Float32,2} @ CPU0:
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> push!(kv, 42, ∇W)

julia> pull!(kv, 42, W)  # fetch weight and write back to `W`

julia> W
2×3 mx.NDArray{Float32,2} @ CPU0:
 -0.2  -0.2  -0.2
 -0.2  -0.2  -0.2
```
"""
function setoptimizer!(kv::KVStore, opt::AbstractOptimizer)
  if occursin(r"dist", string(get_type(kv))) && _isworker()
    # TODO
    error("not implemented")
  else
    setupdater!(kv, getupdater(opt))
  end
end

function _isworker()::Bool
  ref = Ref{Cint}(0)
  @mxcall(:MXKVStoreIsWorkerNode, (Ref{Cint},), ref)
  ref_is_worker[]
end

# TODO: sparse support?


text:
using Arpack, SparseArrays, LinearAlgebra
using Plots

# define functions used here
include("ades.jl")
include("acre.jl")
include("pascaltriangle.jl")
include("in2b.jl")
# include("b2in.jl")
include("epsilon.jl")
include("cutMsizeEne.jl")
include("Hintfunccutoff.jl")
include("delta.jl")

function Hsocfunccutoff!(indvec::Vector{Int64}, Msize0::Int64, Np::Int64, matp::Matrix{Int64}, energyijHOmat::Vector{Float64}, energyijSOCmat::Matrix{ComplexF64}, energyijWmat::Matrix{Float64}, ksoc::Float64, Omega::Float64, matHO::SparseMatrixCSC{Float64}, matSOC::SparseMatrixCSC{ComplexF64}, matW::SparseMatrixCSC{Float64})

    Msize = Msize0*2
    maxmatpcut = length(indvec)

    # matp = zeros(Int,Msize+1,Np+1);
    # pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # maxmatp = matp[Msize+1,Np+1] # the indices are m+1 and n+1 for N^m_ns

    # defines vectors and matrices
    # Hsoc .= 0; #Hsoc = spzeros(ComplexF64,maxmatpcut,maxmatpcut); #spzeros(ComplexF64,maxmatp,maxmatp);
    matHO .= 0.0
    matSOC .= 0.0
    matW .= 0.0

    # define a matrix for the Hamiltonian
    for nn = 1:maxmatpcut #maxmatp

        vecmbnn = spzeros(Int64,Msize+1)
        in2b!(indvec[nn],Msize,Np,matp,vecmbnn) #vecmbnn = in2b(indvec[nn],Msize,Np,matp) #in2b(nn,Msize,Np)
        energyij = 0.

        # free terms
        for jj = 1:Msize

            vecmbnnj = spzeros(Int64,Msize+1)
            ades!(jj,vecmbnn,vecmbnnj) #vecmbnnj = ades(jj,vecmbnn)
            if vecmbnnj[Msize+1] == 0
               continue
            end

            for ii = 1:Msize

                vecmbnnij = spzeros(Int64,Msize+1)
                acre!(ii,vecmbnnj,vecmbnnij) #vecmbnnij = acre(ii,vecmbnnj)

                # energyij = epsilon(ii,jj,ksoc,Omega)
                # energyij = energyijmat[ii,jj]
                energyij = energyijHOmat[ii]*delta(ii-jj) + energyijSOCmat[ii,jj] + energyijWmat[ii,jj]
                if isapprox(energyij,0) #energyij == 0 # energy is Int and zero if ii==jj
                   continue #vecmbnnij[1+Msize] = 0
                end

                for mm = 1:nn

                    vecmbmm = spzeros(Int64,Msize+1)
                    in2b!(indvec[mm],Msize,Np,matp,vecmbmm) #vecmbmm = in2b(indvec[mm],Msize,Np)

                    if vecmbnnij[1:Msize] == vecmbmm[1:Msize]

                       # Hsoc[mm,nn] = Hsoc[mm,nn] + energyij*sqrt(vecmbnnij[Msize+1])
                       matHO[mm,nn] = matHO[mm,nn] + energyijHOmat[ii]*delta(ii-jj)*sqrt(vecmbnnij[Msize+1])
                       matSOC[mm,nn] = matSOC[mm,nn] + energyijSOCmat[ii,jj]*sqrt(vecmbnnij[Msize+1])
                       matW[mm,nn] = matW[mm,nn] + energyijWmat[ii,jj]*sqrt(vecmbnnij[Msize+1])

                    end

                end

            end

        end

        # # Interactions
        # for ll = 1:Msize
        #
        # end

    end

    # since the Hamiltonian is helmitian
    # Hsoc .= Hsoc + Hsoc' - spdiagm(diag(Hsoc))
    # matHO .=
    matSOC .= matSOC + matSOC' - spdiagm(diag(matSOC))
    matW .= matW + matW' - spdiagm(diag(matW))

    # return Hsoc

end

function main(gdown::Float64, gup::Float64, gdu::Float64, ksoc::Float64, Omega0::Float64, Omega1::Float64, NOmega::Int64, Msize0::Int64, Np::Int64, Ene0minumhalf::Int64, specnum::Int64)

    # define cutoff of energy in Fock states
    Msize = Msize0*2
    matp = zeros(Int64,Msize+1,Np+1)
    pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # note the indices are m+1 and n+1 for N^m_n
    indvec = cutMsizeEne(Msize0,Np,matp,Ene0minumhalf)
    maxmatpcut = length(indvec)

    # interaction Hamiltonian
    # mat1, mat2, mat3 = Hintfunccutoff(indvec,Msize0,Np)
    mat1 = spzeros(Float64,maxmatpcut,maxmatpcut)
    mat2 = spzeros(Float64,maxmatpcut,maxmatpcut)
    mat3 = spzeros(Float64,maxmatpcut,maxmatpcut)
    Hintfunccutoff!(indvec,Msize0,Np,matp,mat1,mat2,mat3)
    # lambda, phi = eigs(mat0+gdown*mat1+gup*mat2+gdu*mat3,nev=specnum,which=:SR)

    # single-particle Hamiltonian
    # calculate and save energyij
    energyijmat = zeros(Float64,Msize,Msize);
    energyijHOmat = zeros(Float64,Msize);
    energyijSOCmat = zeros(ComplexF64,Msize,Msize);
    energyijWmat = zeros(Float64,Msize,Msize);
    for ii = 1:Msize
        for jj = 1:ii

            # energyijmat[ii,jj] = epsilon(ii,jj,ksoc,Omega)

            if ii == jj
               energyijHOmat[ii] = epsilon(ii,jj,0.0,0.0)
            elseif iseven(ii+jj)
               energyijSOCmat[ii,jj] = epsilon(ii,jj,1.0,0.0)
            elseif abs(ii-jj) == 1 && ceil(Int64,ii/2) == ceil(Int64,jj/2) #nii == njj # && isodd(ii+jj)
               energyijWmat[ii,jj] = epsilon(ii,jj,0.0,1.0)
            end

        end
    end
    # energyijmat .= energyijmat + energyijmat' + spdiagm(diag(energyijmat))
    energyijSOCmat .= energyijSOCmat + energyijSOCmat' + spdiagm(diag(energyijSOCmat))
    energyijWmat .= energyijWmat + energyijWmat' + spdiagm(diag(energyijWmat))

    matHO = spzeros(Float64,maxmatpcut,maxmatpcut)
    matSOC = spzeros(ComplexF64,maxmatpcut,maxmatpcut)
    matW = spzeros(Float64,maxmatpcut,maxmatpcut)
    Hsocfunccutoff!(indvec,Msize0,Np,matp,energyijHOmat,energyijSOCmat,energyijWmat,1.0,1.0,matHO,matSOC,matW)

    # make an array for different Omega
    arrayOmega = LinRange(Omega0,Omega1,NOmega)
    lambdajj = zeros(Float64,specnum,NOmega)

    # different Omega
    for jj = 1:NOmega

        Omegajj = arrayOmega[jj]

        # diagonalise
        lambda, phi = eigs(matHO+matSOC*ksoc+matW*Omegajj+gdown*mat1+gup*mat2+gdu*mat3,nev=specnum,which=:SR)
        lambdajj[:,jj] .= real(lambda[:])

    end

    return lambdajj, arrayOmega

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
using DictUtils

x = ["bar"=>0.0,"baz"=>1.0]
y = ["foo"=>1.0,"bar"=>1.0,"baz"=>1.0]
u = ["foo"=>1.0,"bar"=>2.0,"baz"=>3.0]
v = ["foo"=>1.0,"bar"=>2.0,"baz"=>3.0]

@assert coalesce(1,null,2) == 1
@assert coalesce(null,1,2) == 1
@assert coalesce(null,null,2) == 2

addDicts!(x, y, 1.0)
@assert x == ["foo"=>1,"bar"=>1,"baz"=>2]

normalizeDict!(x)
@assert x == ["foo"=>0.25,"bar"=>0.25,"baz"=>0.5]

sortedHead!(u, 4)
@assert u == ["foo"=>1,"bar"=>2,"baz"=>3]
sortedHead!(u, 2)
@assert u == ["foo"=>1,"bar"=>2]

sortedTail!(v, 4)
@assert v == ["foo"=>1,"bar"=>2,"baz"=>3]
sortedTail!(v, 2)
@assert v == ["bar"=>2,"baz"=>3]

@assert showCompact(v) == "bar:2.0,baz:3.0"
@assert parsePrimitive(String, "3") == "3"
@assert parsePrimitive(Float64, "3") == 3.0
@assert parsePrimitive(Int64, "3") == 3
@assert parseCompact("bar:2,baz:3", String, Float64) == v


text:
abstract type Order end

"""
Trait container for dimension and array ordering in AlignedGrid.

The default is `Ordered(Forward()`, `Forward())`

All combinations of forward and reverse order for data and indices seem to occurr
in real datasets, as strange as that seems. We cover these possibilities by specifying
the order of both explicitly.

Knowing the order of indices is important for using methods like `searchsortedfirst()`
to find indices in sorted lists of values. Knowing the order of the data is then
required to map to the actual indices. It's also used to plot the data later - which
always happens in smallest to largest order.

Base also defines Forward and Reverse, but they seem overly complicated for our purposes.
"""
struct Ordered{D,A,R} <: Order
    index::D
    array::A
    relation::R
end
Ordered(; index=Forward(), array=Forward(), relation=Forward()) =
    Ordered(index, array, relation)

indexorder(order::Ordered) = order.index
arrayorder(order::Ordered) = order.array
relationorder(order::Ordered) = order.relation

"""
Trait indicating that the array or dimension has no order.
"""
struct Unordered{R} <: Order
    relation::R
end
Unordered() = Unordered(Forward())

indexorder(order::Unordered) = order
arrayorder(order::Unordered) = order
relationorder(order::Unordered) = order.relation

"""
Trait indicating that the array or dimension is in the normal forward order.
"""
struct Forward <: Order end

"""
Trait indicating that the array or dimension is in the reverse order.
Selector lookup or plotting will be reversed.
"""
struct Reverse <: Order end

Base.reverse(::Reverse) = Forward()
Base.reverse(::Forward) = Reverse()
# Base.reverse(o::Ordered) =
    # Ordered(indexorder(o), reverse(relationorder(o)), reverse(arrayorder(o)))
# Base.reverse(o::Unordered) =
    # Unordered(reverse(relationorder(o)))

reverseindex(o::Unordered) =
    Unordered(reverse(relationorder(o)))
reverseindex(o::Ordered) =
    Ordered(reverse(indexorder(o)), arrayorder(o), reverse(relationorder(o)))

reversearray(o::Unordered) =
    Unordered(reverse(relationorder(o)))
reversearray(o::Ordered) =
    Ordered(indexorder(o), reverse(arrayorder(o)), reverse(relationorder(o)))

isrev(::Forward) = false
isrev(::Reverse) = true

"""
Indicates wether the cell value is specific to the locus point
or is related to the whole the span.

The span may contain a value if the distance between locii if known.
This will often be identical to the distance between any two sequential
cell values, but may be distinct due to rounding errors in a vector index,
or context-dependent spans such as `Month`.
"""
abstract type Sampling end

"""
Each cell value represents a siegle discrete sample taken at the index location.
"""
struct SingleSample <: Sampling end

"""
Multiple samples from the span combined using method `M`,
where `M` is `typeof(mean)`, `typeof(sum)` etc.
"""
struct MultiSample{M} <: Sampling end
MultiSample() = MultiSample{Nothing}()

"""
The sampling method is unknown.
"""
struct UnknownSampling <: Sampling end

"""
Indicate the position of index values in grid cells.

This is frequently `Start` for time series, but may be `Center`
for spatial data.
"""
abstract type Locus end

"""
Indicates dimensions that are defined by their center coordinates/time/position.
"""
struct Center <: Locus end

"""
Indicates dimensions that are defined by their start coordinates/time/position.
"""
struct Start <: Locus end

"""
Indicates dimensions that are defined by their end coordinates/time/position.
"""
struct End <: Locus end

struct UnknownLocus <: Locus end

"""
Traits describing the grid type of a dimension.
"""
abstract type Grid end

dims(g::Grid) = nothing
arrayorder(grid::Grid) = arrayorder(order(grid))
indexorder(grid::Grid) = indexorder(order(grid))
relationorder(grid::Grid) = relationorder(order(grid))

Base.reverse(g::Grid) = rebuild(g, reverse(order(g)))
reversearray(g::Grid) = rebuild(g, reversearray(order(g)))
reverseindex(g::Grid) = rebuild(g, reverseindex(order(g)))

slicegrid(grid::Grid, index, I) = grid

"""
Fallback grid type
"""
struct UnknownGrid <: Grid end

order(::UnknownGrid) = Unordered()

# TODO handle unordered in all grids
bounds(grid::UnknownGrid, dim) = first(dim), last(dim)

"""
A grid dimension that is independent of other grid dimensions.
"""
abstract type IndependentGrid{O} <: Grid end

"""
A grid dimension aligned exactly with a standard dimension, such as lattitude or longitude.
"""
abstract type AbstractAlignedGrid{O,L,Sa} <: IndependentGrid{O} end

order(g::AbstractAlignedGrid) = g.order
locus(g::AbstractAlignedGrid) = g.locus
sampling(g::AbstractAlignedGrid) = g.sampling

"""
An [`AlignedGrid`](@ref) grid without known regular spacing. These grids will generally be paired
with a vector of coordinates along the dimension, instead of a range.

Bounds are given as the first and last points, which omits the span of one cell,
as it is not known. To fix this use either a [`BoundedGrid`](@ref) with specified
starting bounds or a [`RegularGrid`](@ref) with a known constand cell span.

## Fields
- `order::Order`: `Order` trait indicating array and index order
- `locus::Locus`: `Locus` trait indicating the position of the indexed point within the cell span
- `sampling::Sampling`: `Sampling` trait indicating wether the grid cells are single samples or means
"""
struct AlignedGrid{O<:Order,L<:Locus,Sa<:Sampling} <: AbstractAlignedGrid{O,L,Sa}
    order::O
    locus::L
    sampling::Sa
end
AlignedGrid(; order=Ordered(), locus=Start(), sampling=UnknownSampling()) =
    AlignedGrid(order, locus, sampling)

rebuild(g::AlignedGrid, order=order(g), locus=locus(g), sampling=sampling(g)) =
    AlignedGrid(order, locus, sampling)

bounds(grid::AlignedGrid, dim) = first(dim), last(dim)

"""
An alligned grid without known regular spacing and tracked bounds.
These grids will generally be paired with a vector of coordinates along the
dimension, instead of a range.

As the size of the cells is not known, the bounds must be actively tracked.

## Fields
- `order::Order`: `Order` trait indicating array and index order
- `locus::Locus`: `Locus` trait indicating the position of the indexed point within the cell span
- `sampling::Sampling`: `Sampling` trait indicating wether the grid cells are single samples or means
- `bounds`: the outer edges of the grid (different to the first and last coordinate).
"""
struct BoundedGrid{O<:Order,L<:Locus,Sa<:Sampling,B} <: AbstractAlignedGrid{O,L,Sa}
    order::O
    locus::L
    sampling::Sa
    bounds::B
end
BoundedGrid(; order=Ordered(), locus=Start(), sampling=UnknownSampling(), bounds=nothing) =
    BoundedGrid(order, locus, sampling, bounds)

bounds(g::BoundedGrid) = g.bounds

rebuild(g::BoundedGrid, order=order(g), locus=locus(g), sampling=sampling(g), bounds=bounds(g)) =
    BoundedGrid(order, locus, sampling, bounds)

# TODO: deal with unordered AbstractArray
slicegrid(g::BoundedGrid, index, I) =
    rebuild(g, order(g), locus(g), sampling(g), slicebounds(locus(g), bounds(g), index, I))

slicebounds(loci::Start, bounds, index, I) =
    index[first(I)],
    last(I) >= lastindex(index)  ? bounds[2] : index[last(I) + 1]
slicebounds(loci::End, bounds, index, I) =
    first(I) <= firstindex(index) ? bounds[1] : index[first(I) - 1],
    index[last(I)]
slicebounds(loci::Center, bounds, index, I) =
    first(I) <= firstindex(index) ? bounds[1] : (index[first(I) - 1]   + index[first(I)]) / 2,
    last(I)  >= lastindex(index)  ? bounds[2] : (index[last(I) + 1] + index[last(I)]) / 2

abstract type AbstractRegularGrid{O,L,Sa,Sp} <: AbstractAlignedGrid{O,L,Sa} end

"""
An [`AlignedGrid`](@ref) where all cells are the same size and evenly spaced.

## Fields
- `order::Order`: `Order` trait indicating array and index order
- `locus::Locus`: `Locus` trait indicating the position of the indexed point within the cell span
- `sampling::Sampling`: `Sampling` trait indicating wether the grid cells are single samples or means
- `span::Number`: the size of a grid step, such as 1u"km" or `Month(1)`
"""
struct RegularGrid{O<:Order,L<:Locus,Sa<:Sampling,Sp} <: AbstractRegularGrid{O,L,Sa,Sp}
    order::O
    locus::L
    sampling::Sa
    span::Sp
end
RegularGrid(; order=Ordered(), locus=Start(), sampling=UnknownSampling(), span=nothing) =
    RegularGrid(order, locus, sampling, span)

span(grid::RegularGrid) = grid.span

rebuild(g::RegularGrid, order=order(g), locus=locus(g), sampling=sampling(g), span=span(g)) =
    RegularGrid(order, locus, sampling, span)

bounds(grid::RegularGrid, dim) = bounds(indexorder(grid), locus(grid), grid, dim)
bounds(::Forward, ::Start, grid, dim) = first(dim), last(dim) + span(grid)
bounds(::Reverse, ::Start, grid, dim) = last(dim), first(dim) + span(grid)
bounds(::Forward, ::Center, grid, dim) = first(dim) - span(grid) / 2, last(dim) + span(grid) / 2
bounds(::Reverse, ::Center, grid, dim) = last(dim) - span(grid) / 2, first(dim) + span(grid) / 2
bounds(::Forward, ::End, grid, dim) = first(dim) - span(grid), last(dim)
bounds(::Reverse, ::End, grid, dim) = last(dim) - span(grid), first(dim)

abstract type AbstractCategoricalGrid{O} <: IndependentGrid{O} end

"""
A grid dimension where the values are categories.

## Fields
- `order`: `Order` trait indicating array and index order
"""
struct CategoricalGrid{O<:Order} <: AbstractCategoricalGrid{O}
    order::O
end
CategoricalGrid(; order=Ordered()) = CategoricalGrid(order)

order(g::CategoricalGrid) = g.order

rebuild(g::CategoricalGrid, order=order(g)) = CategoricalGrid(order)

bounds(grid::CategoricalGrid, dim) = bounds(indexorder(grid), grid, dim)
bounds(::Forward, grid, dim) = first(dim), last(dim)
bounds(::Reverse, grid, dim) = last(dim), first(dim)
bounds(::Unordered, grid, dim) = error("Cannot call `bounds` on an unordered grid")

"""
Traits describing a grid dimension that is dependent on other grid dimensions.

Indexing into a dependent dimension must provide all other dependent dimensions.
"""
abstract type DependentGrid <: Grid end

locus(g::DependentGrid) = g.locus
dims(g::DependentGrid) = g.dims
sampling(g::DependentGrid) = g.sampling

"""
Grid type using an affine transformation to convert dimension from
`dim(grid)` to `dims(array)`.

## Fields
- `dims`: a tuple containing dimenension types or symbols matching the order
          needed by the transform function.
- `sampling`: a `Sampling` trait indicating wether the grid cells are sampled points or means
"""
struct TransformedGrid{D,L,Sa<:Sampling} <: DependentGrid
    dims::D
    locus::L
    sampling::Sa
end
TransformedGrid(dims=(), locus=Start(), sampling=UnknownSampling()) =
    TransformedGrid(dims, locus, sampling)

rebuild(g::TransformedGrid, dims=dims(g), locus=locus(g), sampling=sampling(g)) =
    TransformedGrid(dims, locus, sampling)

"""
A grid dimension that uses an array lookup to convert dimension from
`dim(grid)` to `dims(array)`.

## Fields
- `dims`: a tuple containing dimenension types or symbols matching the order
          needed to index the lookup matrix.
- `sampling`: a `Sampling` trait indicating wether the grid cells are sampled points or means
"""
struct LookupGrid{D,L,Sa<:Sampling} <: DependentGrid
    dims::D
    locus::L
    sampling::Sa
end
LookupGrid(dims=(), locus=Start(), sampling=UnknownSampling()) =
    LookupGrid(dims, locus, sampling)

rebuild(g::LookupGrid; dims=dims(g), locus=locus(g), sampling=sampling(g)) =
    LookupGrid(dims, locus, sampling)


text:
@testset "File splitting" begin

    W = addprocs(8)
    @everywhere using GigaSOM

    # re-use the PBMC8 data:
    # try load it by parts of different size and scrutinize the parts

    for splitSize = 1:8
        di = loadFCSSet(:test, md[:, :file_name], W[1:splitSize])
        dselect(di, fcsAntigens, antigens)
        cols = Vector(1:length(antigens))
        dtransform_asinh(di, cols, 5)
        dscale(di, cols)
        @test isapprox(pbmc8_data, gather_array(di), atol = 1e-4)

        splits = dmapreduce(di, d -> size(d, 1), vcat)
        @test sum(splits) == size(pbmc8_data, 1)
        @test minimum(splits) + 1 >= maximum(splits)

        sizes = loadFCSSizes(md[:, :file_name])
        dis = distributeFCSFileVector(:testF, md[:, :file_name], W[1:splitSize])
        @test dcount(length(sizes), dis) == sizes
        @test dcount_buckets(length(sizes), dis, length(sizes), dis) ==
              Matrix{Int}(LinearAlgebra.Diagonal(sizes))
    end

    rmprocs(W)
end


text:
immutable Exponential <: ContinuousUnivariateDistribution
    β::Float64 		# note: scale not rate

    function Exponential(β::Real)
        β > zero(β) || error("scale must be positive")
        new(float64(β))
    end

    Exponential() = new(1.0)
end

@distr_support Exponential 0.0 Inf

#### Parameters

scale(d::Exponential) = d.β
rate(d::Exponential) = 1.0 / d.β

params(d::Exponential) = (d.β,)

#### Statistics

mean(d::Exponential) = scale(d)

median(d::Exponential) = logtwo * scale(d)

mode(d::Exponential) = 0.0

var(d::Exponential) = scale(d)^2

skewness(d::Exponential) = 2.0

kurtosis(d::Exponential) = 6.0

entropy(d::Exponential) = 1.0 + log(scale(d))

#### Evaluation

zval(d::Exponential, x::Float64) = x / d.β
xval(d::Exponential, z::Float64) = z * d.β

pdf(d::Exponential, x::Float64) = (λ = rate(d); x < 0.0 ? 0.0 : λ * exp(-λ * x))
logpdf(d::Exponential, x::Float64) =  (λ = rate(d); x < 0.0 ? -Inf : log(λ) - λ * x)

cdf(d::Exponential, x::Float64) = x > 0.0 ? -expm1(-zval(d, x)) : 0.0
ccdf(d::Exponential, x::Float64) = x > 0.0 ? exp(-zval(d, x)) : 0.0
logcdf(d::Exponential, x::Float64) = x > 0.0 ? log1mexp(-zval(d, x)) : -Inf
logccdf(d::Exponential, x::Float64) = x > 0.0 ? -zval(d, x) : 0.0

quantile(d::Exponential, p::Float64) = -xval(d, log1p(-p))
cquantile(d::Exponential, p::Float64) = -xval(d, log(p))
invlogcdf(d::Exponential, lp::Float64) = -xval(d, log1mexp(lp))
invlogccdf(d::Exponential, lp::Float64) = -xval(d, lp)

gradlogpdf(d::Exponential, x::Float64) = x > 0.0 ? -rate(d) : 0.0

mgf(d::Exponential, t::Real) = 1.0/(1.0 - t * scale(d))
cf(d::Exponential, t::Real) = 1.0/(1.0 - t * im * scale(d))

#### Sampling

rand(d::Exponential) = xval(d, randexp())

#### Fit model

immutable ExponentialStats <: SufficientStats
    sx::Float64   # (weighted) sum of x
    sw::Float64   # sum of sample weights

    ExponentialStats(sx::Real, sw::Real) = new(float64(sx), float64(sw))
end

suffstats(::Type{Exponential}, x::Array) = ExponentialStats(sum(x), length(x))
suffstats(::Type{Exponential}, x::Array, w::Array) = ExponentialStats(dot(x, w), sum(w))

fit_mle(::Type{Exponential}, ss::ExponentialStats) = Exponential(ss.sx / ss.sw)