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

using Markdown
using InteractiveUtils

# ╔═╡ 61e06d99-e14d-462f-a83a-bffd734408f6
using Plots

# ╔═╡ 2b1cdf75-d93d-413a-b3f7-37eb9ea54910
md"""
In the additional file ("PS14\_MDP\_spt.jl"), you can find a generic implentation for dealing with MDP's.
"""

# ╔═╡ ee9e8098-5bb9-4598-bc7d-095412b8cd73
"""
	ingredients(path::String)

Helper function to reload external code. The external file is processed and
returned as a module.

### Example
```
# load a external file
mymodule = ingredients("path/to/external_file.jl")
# use the external functions
mymodule.foo(args)
```
"""
function ingredients(path::String)
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	return m
end

# ╔═╡ 9c9e40ce-6ed6-4dd5-951e-074e6bd3a2da
MDP = ingredients("./PS14_MDP_spt.jl")

# ╔═╡ eac907a6-b230-11eb-04f1-5df27e70da8b
md"""
# Complex decisions
A sequential decision problem for a fully observable, stochastic environment with a Markovian transition model and additive rewards is called a Markov decision process (MDP). It consists of the following items:
* a set of states (with its initial state $s_0$)
* a set of actions $\text{ACTION}(s)$ of actions in each state
* a transition model $P(s'|s,a)$
* a reward function $R(s)$

In MDP chance node represent uncertainty about what might happen based on (state,action). A Q-state (s,a) is when you were in a state and took an action.

"""

# ╔═╡ 981096f5-6948-4355-93b2-0fcf453213a5
md"""
## Transforming a problem into an MDP
Below we have a set of problem that we want to transform into the MDP formalism. Think about the states, transition probabilties and rewards. Provide the components of the MDP formalism in a graphical way.
1. a cleaning robot needs to collect trash. The robot has limited autonomy. Possible activities include include searching for trash, waiting  or recharging itself.
2. a human wants to maximize its earnings. The human can be tired, energetic or healthier. in function of how the human is feeling, he can sleep, work (=earn money) or do a workout. Working always fatigues the humen. Doing a workout makes him feel more energetic when tired and makes him healthier when energetic.
3. consider the example of the lectures with the robot in the 4x3 grid. provide the MDP formalism.

"""

# ╔═╡ 46e7fb69-8bd0-4a78-a7df-83e71cde2abc
md"""
## Solving an MDP
We can define a MDP from the example from the course by defining the initial and terminal states as well as the rewards for each state.
"""

# ╔═╡ 06eb94f5-18b8-420b-bb37-c4bed65672ec
begin
	# environment setup
	initial = (1, 1)
	terminal_states = Set([(2, 4), (3, 4)])
	grid = [-0.04 -0.04 -0.04 -0.04;
			-0.04 nothing -0.04 -1;
			-0.04 -0.04 -0.04 +1]
	# environment creation
	myMDP = MDP.GridMarkovDecisionProcess(initial, terminal_states, grid)
end

# ╔═╡ 344b2ed1-c4a5-427d-811d-4a2757c9677d
md"""
If you want to visualize a specific set of actions, you can use the `show_grid` function. If you have a set of actions that leads to a solution and describes what to do for each state, you have a policy for the MDP.
"""

# ╔═╡ b8e0a622-5f7d-4ebc-8758-7907cee24f34
begin
	# Illustration of limited actions (maps: state => action)
	myactions = Dict((1,1) => (0,1), (1,2) => (1,0) )
	# show results
	gridstring = MDP.show_grid(myMDP,myactions)
	println("Grid layout and actions:")
	println(gridstring)
end

# ╔═╡ 22070e43-650c-41b4-b4de-050da37f83fe
md"""
We can verify the behavior of our model and the different function for some small examples
"""

# ╔═╡ 7c933a68-a618-4df2-9e73-bff1c3faf4b6
begin
	pos = (2,3)
	action=(1,0)
	# all possible actions are shown
	MDP.actions(myMDP, pos)
end

# ╔═╡ a0ab6856-ce14-4f18-aee1-efb137a89ca8
md"""
Altough all actions are identified, invalid postions and action are taken into account in the transition model. You can clearly see this in the example below. You would expect to find (2,2) as a valid case, but as this state does not exist, you get staying in place instead.
"""

# ╔═╡ 3b7213c7-c782-4686-9dc2-a0d522587d3a
tt = MDP.transition_model(myMDP,pos,action)

# ╔═╡ beff15ad-ea24-44a4-912d-5101759ac937
md"""
You can observe something similar when you try an action that does not make sense e.g. trying to quit the space in a corner:
"""

# ╔═╡ f7f02056-6f87-4bbf-b157-a44693b47a6b
let
	state = (4,4)
	[MDP.transition_model(myMDP,state,action) for action in MDP.actions(myMDP,state)]
end

# ╔═╡ 03712063-cbdf-4680-ae77-994d37970674
md"""
Our transition model works as intended.
"""

# ╔═╡ c7c601d8-408e-4227-9c56-0e30f429eb50
md"""
### Value iteration
Starting from the Bellman equation you can solve the $n$ equations for the $n$ states in an iterative fashion.

```math
V_{k+1}(s) = \max_{a}\sum_{s'}T(s,a,s')[R(s,a,s') + \gamma V_k(s') ]
```
where
```math
\begin{eqnarray*}
T(s,a,s')&:&\text{transition probility of going to state $s'$ when doing action $a$ in state $s$}\\
R(s,a,s') &:&\text{reward}\\
\gamma&:&\text{discount factor}
\end{eqnarray*}
```
This algorithm will converge to the unique optimal values $V^{*}(s)$ after a number of iterations after being initialized to zero.
"""

# ╔═╡ 612b57ce-549e-430c-b040-42aa13830283
VI = MDP.value_iteration(myMDP)

# ╔═╡ 6ce9ef18-2342-4cab-86e7-9e97aeb670f8
md"""
Now that we have the final values, we can illustrate it
"""

# ╔═╡ 6e90ab1e-dd41-4228-9e01-d3bfada8c387
md"""
#### Obtaining a policy from value iteration
Once you have the solutions of the problem, you can obtain the optimal policy by doing policy extraction:
```math
\pi^{*}(s)=  \underset{a}{\text{argmax}} \sum_{s'}T(s,a,s')[R(s,a,s') + \gamma V^{*}(s,a,s')]
```
"""

# ╔═╡ 8f7ca271-87da-4df5-a80b-ebc59ed121db
pol = MDP.policy_extraction(myMDP, VI)

# ╔═╡ d1734b14-9e9f-4208-9d9b-8be19161f3dc
let
	F = zeros(size(myMDP.grid))
	for (key,val) in VI
		F[key...] = val
	end
	p = heatmap(F,seriescolor=:RdYlGn_11,yflip=true)
	title!("value iteration results")
	yticks!(1:3)
	for (key,val) in pol
		dir = text("""$(get(MDP.arrow_characters,val,"?"))""", :center, 28)
		annotate!((key[2],key[1]-1/6, dir))
		annotate!((key[2],key[1]+1/4, """$(round(get(VI,key,0),digits=4))"""))
	end
	p
end

# ╔═╡ 29226fd6-893e-4ca2-8bfb-6e834a9530bf
md"""
### Policy iteration
In the previous section we have used the non-linear $\max$ operator. If you impose a policy on the system, you no longer have to do this and this reduces the computation from $\mathcal{O}(S^2A)$ to $\mathcal{O}(S^2)$. The update is given by

```math
V_{k+1}(s) = \sum_{s'}T(s,\pi(s),s')[R(s,\pi(s),s') + \gamma V_k^{\pi}(s') ]
```
We end up with a linear system that can be solved exactly $\mathcal{O}(n^3)$ or in an iterative, approximate way.

Now we set a fixed policy to start and solve the problem again.
"""

# ╔═╡ 1521e183-e15e-4ac1-ac5b-6a058914bd05
polbis,ubis = MDP.policy_iteration(myMDP)

# ╔═╡ 40aa0d72-894f-473f-b644-15fcba096aba
let
	F = zeros(size(myMDP.grid))
	for (key,val) in ubis
		F[key...] = val
	end
	p = heatmap(F,seriescolor=:RdYlGn_11,yflip=true)
	title!("policy iteration results")
	yticks!(1:3)
	for (key,val) in polbis
		dir = text("""$(get(MDP.arrow_characters,val,"?"))""", :center, 28)
		annotate!((key[2],key[1]-1/6, dir))
		annotate!((key[2],key[1]+1/4, """$(round(get(VI,key,0),digits=4))"""))
	end
	p
end

# ╔═╡ b008fbdc-10c3-4f67-b124-9870fc7eaa58
md"""
**Note:**

In the previous example, the rewards of a state depends only on the state itself, i.e. $R(s)$. In practice however, you can encounter a reward that also depends on the action taken, i.e. $R(s,a)$, or even on the action taken and the outcome state, i.e. $R(s,a,s')$.
"""

# ╔═╡ 6dd8d501-ee40-4409-a5b8-5835297fdde0
md"""
### Application 1
Consider the human life cycle from before.

* What is different with respect to the previous application?
* Implement it and analyze the impact of the parameters in the model on the policy you would apply. Do your results make sense?

"""

# ╔═╡ Cell order:
# ╠═61e06d99-e14d-462f-a83a-bffd734408f6
# ╟─2b1cdf75-d93d-413a-b3f7-37eb9ea54910
# ╟─ee9e8098-5bb9-4598-bc7d-095412b8cd73
# ╠═9c9e40ce-6ed6-4dd5-951e-074e6bd3a2da
# ╟─eac907a6-b230-11eb-04f1-5df27e70da8b
# ╟─981096f5-6948-4355-93b2-0fcf453213a5
# ╟─46e7fb69-8bd0-4a78-a7df-83e71cde2abc
# ╠═06eb94f5-18b8-420b-bb37-c4bed65672ec
# ╟─344b2ed1-c4a5-427d-811d-4a2757c9677d
# ╠═b8e0a622-5f7d-4ebc-8758-7907cee24f34
# ╟─22070e43-650c-41b4-b4de-050da37f83fe
# ╠═7c933a68-a618-4df2-9e73-bff1c3faf4b6
# ╟─a0ab6856-ce14-4f18-aee1-efb137a89ca8
# ╠═3b7213c7-c782-4686-9dc2-a0d522587d3a
# ╟─beff15ad-ea24-44a4-912d-5101759ac937
# ╠═f7f02056-6f87-4bbf-b157-a44693b47a6b
# ╟─03712063-cbdf-4680-ae77-994d37970674
# ╟─c7c601d8-408e-4227-9c56-0e30f429eb50
# ╠═612b57ce-549e-430c-b040-42aa13830283
# ╟─6ce9ef18-2342-4cab-86e7-9e97aeb670f8
# ╟─6e90ab1e-dd41-4228-9e01-d3bfada8c387
# ╠═8f7ca271-87da-4df5-a80b-ebc59ed121db
# ╠═d1734b14-9e9f-4208-9d9b-8be19161f3dc
# ╟─29226fd6-893e-4ca2-8bfb-6e834a9530bf
# ╠═1521e183-e15e-4ac1-ac5b-6a058914bd05
# ╠═40aa0d72-894f-473f-b644-15fcba096aba
# ╟─b008fbdc-10c3-4f67-b124-9870fc7eaa58
# ╟─6dd8d501-ee40-4409-a5b8-5835297fdde0


text:
# Purpose: identify components of phenotype neutral network.

# steps is the number of random walk steps.
# maxsteps is the number of steps in mut_evolve().
# maxtries is the number of attempts to find the next chromosome in the random_neutral_walk.
# Assumes a single output goal
# Returns circuit_code_list of neutral circuits discovered on the neutral walk.
# Includes both walk circuits and neutral neighbors of walk circuits.
function neutral_walk( g::Goal, p::Parameters, steps::Int64, maxsteps::Int64, maxtries::Int64 )
  n_repeats = 10  # maximum number of tries to find c that maps to g
  res = mut_evolve_repeat(n_repeats, p, [g], funcs, maxsteps )
  c = res[1]
  neutral_walk( c, p, steps, maxsteps, maxtries )
end

# neutral walk starting with a given chromosome instead of a given circuit.
function neutral_walk( c::Chromosome, p::Parameters, steps::Int64, maxsteps::Int64, maxtries::Int64 )
  @assert p.numoutputs == 1
  funcs = default_funcs( p.numinputs )
  circuit_ints = Int64[]
  n_repeats = 10  # maximum number of tries to find c that maps to g
  res = mut_evolve_repeat(n_repeats, p, [g], funcs, maxsteps )
  c = res[1]
  println("complexity(c): ",complexity5(c))
  #@assert sort(output_values(c)) == sort(g)
  @assert output_values(c) == g    # Assumes p.numoutputs==1
  push!( circuit_ints, circuit_int(c) )
  for i = 1:steps
    (new_c,active) = mutate_chromosome!( deepcopy(c), funcs )
    outputs = output_values(new_c)
    attempts = 1
    while attempts < maxtries && outputs != g  # try to find a c to extend walk
      (new_c,active) = mutate_chromosome!( deepcopy(c), funcs )
      outputs = output_values(new_c)
      #println("attempts: ",attempts,"  outputs: ",outputs,"  goal: ",g)
      attempts += 1
    end
    if attempts == maxtries
      # Alternatively, try all mutations of c in our attempt to extend walk
      all_chromes = map( x->x[2], mutate_all( deepcopy(c), funcs, output_chromosomes=true ))
      filter!( x->output_values(x)==g, all_chromes )  # only chromosomes that map to g
      #println(" attempts == maxtries  len(all_chromes): ",length(all_chromes))
      if length(all_chromes) == 0
        error("unable to extend random_neutral_walk in function neutral_walk() at step: ",i)
        #continue
      else
        @assert output_values(all_chromes[1]) == g
        c = rand(all_chromes)
      end
    else
      c = new_c   # neutral walk extended
    end
    @assert output_values(c) == g
    # mutate_all() returns a list of pairs where the second element of the pair is the chromosome
    all_chromes = map( x->x[2], mutate_all( c, funcs, output_chromosomes=true ))
    filter!( x->output_values(c)==g, all_chromes )
    for ch in all_chromes
      push!( circuit_ints, circuit_int(ch) )
    end
    circuit_ints = unique(circuit_ints)
  end # for loop
  circuit_ints
end

# Does multiple random walks and combines the returned circuit code lists if they have circuits in common.
function run_neutral_walk( g::Goal, p::Parameters, n_walks::Int64, steps::Int64, maxsteps::Int64, maxtries::Int64 )
  walk_list = Int64[]
  circuit_int_list = Set{Int64}[]
  for w = 1:n_walks
    cclist = Set(neutral_walk( g, p, steps, maxsteps, maxtries))
    println("length(cclist): ",length(cclist))
    to_combine = Int64[]   # indices of circuit_int list to combine
    for i = 1:length(circuit_int_list)
      if length(intersect(cclist,circuit_int_list[i])) > 0
        push!(to_combine,i)
      end
    end
    println("to_combine: ",to_combine)
    if length(to_combine) > 0
      combined_list = cclist
      for i = length(to_combine):-1:1
        println("combining circuit_int_list[",to_combine[i],"]")
        combined_list = union(combined_list,circuit_int_list[to_combine[i]])
        if i > 1
          deleteat!( circuit_int_list, to_combine[i] )
          println("removing index: ",to_combine[i]," from circuit_int_list")
          println("walk_list before: ",walk_list)
          #walk_list = setdiff(walk_list,to_combine[i])
          deleteat!( walk_list, to_combine[i] )
          println("removing index: ",to_combine[i]," from walk list")
          println("walk_list after: ",walk_list)
        end
      end
      circuit_int_list[to_combine[1]] = combined_list
    else
      push!(walk_list,w)
      push!(circuit_int_list,cclist)
      #println("walk list: ",walk_list)
      #println("   lengths circuit_int_list: ", [length(circuit_int_list[k]) for k = 1:length(circuit_int_list)])
    end
    #println("w: ",w,"  walk_list after add: ",walk_list)
    #println("   lengths circuit_int_list: ", [length(circuit_int_list[k]) for k = 1:length(circuit_int_list)])
    print("w: ",w,"  length(circuit_int_list) after add: ",length(circuit_int_list))
    println("  lengths circuit_int_list: ",[length(ccl) for ccl in circuit_int_list])
    for i = 1:length(circuit_int_list)
      for j = 1:(i-1)
        if !isempty(intersect(circuit_int_list[i],circuit_int_list[j]))
          println("XXXXX i: ",i,"  j: ",j,"  len intersect: ",length(intersect(circuit_int_list[i],circuit_int_list[j])))
        end
      end
    end
  end
  walk_list
end

# Extract the list of discovered component sizes from an *_ints.txt file created by neutral_walks() function in this file
# Return a dataframe with two fields:  goal, and comp_list.
function extract_component_lists( intsfile::String )
  csdf = DataFrame()
  csdf.goal = String[]
  csdf.comp_list = String[]
  open( intsfile, "r" ) do f
    for line in readlines(f)
      m = match(r"(\[0x.*\]).*(\[.*\])",line)
      push!(csdf,m.captures)
    end
  end
  csdf
end

# Add the component sizes lists returned by function extract_component_sizes() to the dataframe of csvfile.
#  The resulting dataframe is written to a csvfile with "cs" appended to the file name (before the ".csv")
function add_component_sizes_to_df( csvfile::String )
  df = read_dataframe(csvfile)
  ints_file = string( csvfile[1:(end-4)], "_ints.txt" )
  csdf = extract_component_lists( ints_file )
  size_lists = [ eval(Meta.parse(csdf.comp_list[i])) for i = 1:size(csdf)[1] ]
  size_lists = map(x->sort(x,rev=true),size_lists)
  df.size_sum = map( sum, size_lists )
  df.size_ratio = [ (df.size_sum[i] - size_lists[i][1])/df.size_sum[i] for i = 1:size(df)[1] ]
  df.comp_list= csdf.comp_list
  dfcs_file = string( csvfile[1:(end-4)], "ccs.csv" )
  println("dfcs_file: ",dfcs_file )
  write_dataframe_with_comments(df,csvfile,dfcs_file)
end

ext: jl


text:
###################################################
## network/solvers/mip.jl
## MIP to select minimums
###################################################

"""
    `mipTimes` : optimize travel times to minimize L1 error from data with given paths
"""
function mipTimes(s::IterativeState; args...)
    g = s.data.network.graph
    paths = s.paths
    tripData = s.trips

    # Big M = max time of paths -- right now: 5km/h speed
    minTimes = uniformTimes(s.data.network, 5.)
    M = Float64[maximum(pathTime(minTimes,p) for p in l) for l in paths]

    #Create the model (will be changed to avoid hard-coded parameters)
    # !BarConvTol needs to be changed
    m = Model(solver = GurobiSolver(TimeLimit=10000, Method=3, BarConvTol=1e-6; args...))

    # DECISION VARIABLES
    # Road times
    @variable(m, t[i=vertices(g), j=outneighbors(g,i)] >= s.data.minTimes[i,j])
    # Absolute difference between tripData times and computed times
    @variable(m, epsilon[d=eachindex(tripData)] >= 0)
    # Time of the trips
    @variable(m, T[d=eachindex(tripData)] >= 0)
    # Integer variables to decide which path is the shortest
    @variable(m, minP[d=eachindex(tripData),p=1:length(paths[d])], Bin)

    # OBJECTIVE
    @objective(m, Min, sum(tripData[d].weight/tripData[d].time * epsilon[d] for d in eachindex(tripData)))

    # CONSTRAINTS
    # absolute values contraints (define epsilon), equal to time of first path
    @constraint(m, epsLower[d=eachindex(tripData)], T[d] - tripData[d].time >= - epsilon[d])
    @constraint(m, epsUpper[d=eachindex(tripData)], T[d] - tripData[d].time <=   epsilon[d])

    # inequality constraints
    @constraint(m, inequalityPath[d=eachindex(tripData), p=eachindex(paths[d])],
        sum(paths[d][p][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][p])) >= T[d] )

    # "minimum equality" contraints
    @constraint(m, equalityPath[d=eachindex(tripData), p=eachindex(paths[d])],
        sum(paths[d][p][edge] * t[src(edge), dst(edge)] for edge in keys(paths[d][p])) - M[d]*(1 - minP[d,p]) <= T[d])

    # integer constraints
    @constraint(m, equalityPath[d=eachindex(tripData)], sum(minP[d,p] for p in eachindex(paths[d])) == 1)

    # SOLVE LP
    status = solve(m)
    times = getvalue(t)

    # Export result as sparse matrix
    result = spzeros(Float64, nv(g), nv(g))
    for i in vertices(g), j in outneighbors(g,i)
        result[i,j] = times[i,j]
    end

    return result
end


text:
function kmr_markov_matrix_sequential{T<:Real}(n::Integer, p::T, ε::T)
    """
    Generate the markov matrix for the KMR model with *sequential* move

    n: number of players
    p: level of p-dominance for action 2
       = the value of p such that action 2 is the BR for (1-q, q) for any q > p,
         where q (1-q, resp.) is the prob that the opponent plays action 2 (1, resp.)
    ε: mutation probability

    References:
        kmr_markov_matrix_sequential is contributed from https://github.com/oyamad
    """
    P = zeros(T, n+1, n+1)

    P[1, 1], P[1, 2] = 1 - ε/2, ε/2
    @inbounds for i = 1:n-1
        P[i+1, i] = (i/n) * (ε/2 + (1 - ε) *
                             (((i-1)/(n-1) < p) + ((i-1)/(n-1) == p)/2))
        P[i+1, i+2] = ((n-i)/n) * (ε/2 + (1 - ε) *
                                 ((i/(n-1) > p) + (i/(n-1) == p)/2))
        P[i+1, i+1] = 1 - P[i+1, i] - P[i+1, i+2]
    end
    P[end, end-1], P[end, end] = ε/2, 1 - ε/2
    return P
end

function Base.isapprox{T<:Real,S<:Real}(x::Vector{Vector{T}},
                                        y::Vector{Vector{S}})
    n = length(x)
    length(y) == n || return false
    for i in 1:n
        isapprox(x[i], y[i])::Bool || return false
    end
    return true
end

@testset "Testing mc_tools.jl" begin
    # Matrix with two recurrent classes [1, 2] and [4, 5, 6],
    # which have periods 2 and 3, respectively
    Q = [0 1 0 0 0 0
         1 0 0 0 0 0
         1//2 0 0 1//2 0 0
         0 0 0 0 1 0
         0 0 0 0 0 1
         0 0 0 1 0 0]
    Q_stationary_dists = Vector{Rational{Int}}[
        [1//2, 1//2, 0, 0, 0, 0], [0, 0, 0, 1//3, 1//3, 1//3]
    ]
    Q_dict_Rational = Dict(
        "P" => Q,
        "stationary_dists" => Q_stationary_dists,
        "comm_classes" => Vector{Int}[[1, 2], [3], [4, 5, 6]],
        "rec_classes" => Vector{Int}[[1, 2], [4, 5, 6]],
        "is_irreducible" => false,
        "period" => 6,
        "is_aperiodic" => false,
    )
    Q_dict_Float64 = copy(Q_dict_Rational)
    Q_dict_Float64["P"] = convert(Matrix{Float64}, Q)
    Q_dict_Float64["stationary_dists"] =
        convert(Vector{Vector{Float64}}, Q_stationary_dists)

    # examples from
    # Graph-Theoretic Analysis of Finite Markov Chains by J.P. Jarvis & D. R. Shier

    fig1_p = zeros(Rational{Int64}, 5, 5)
    fig1_p[[3, 4, 9, 10, 11, 13, 18, 19, 22, 24]] =
        [1//2, 2//5, 1//10, 1, 1, 1//5, 3//10, 1//5, 1, 3//10]
    fig2_p = zeros(Rational{Int64}, 5, 5)
    fig2_p[[3, 10, 11, 13, 14, 17, 18, 19, 22]] =
        [1//3, 1, 1, 1//3, 1//2, 1//2, 1//3, 1//2, 1//2]

    fig1_dict_Rational = Dict(
        "P" => fig1_p,
        "stationary_dists" => Vector{Rational{Int}}[[0, 1//2, 0, 0, 1//2]],
        "comm_classes" => Vector{Int}[[1, 3, 4], [2, 5]],
        "rec_classes" => Vector{Int}[[2, 5]],
        "is_irreducible" => false,
        "period" => 2,
        "is_aperiodic" => false,
    )
    fig2_dict_Rational = Dict(
        "P" => fig2_p,
        "stationary_dists" => Vector{Rational{Int}}[[1//6, 0, 3//6, 2//6, 0]],
        "comm_classes" => Vector{Int}[[1, 3, 4], [2, 5]],
        "rec_classes" => Vector{Int}[[1, 3, 4]],
        "is_irreducible" => false,
        "period" => 1,
        "is_aperiodic" => true,
    )
    fig1_dict_Float64 = copy(fig1_dict_Rational)
    fig2_dict_Float64 = copy(fig2_dict_Rational)
    for (d_R, d_F) in zip((fig1_dict_Rational, fig2_dict_Rational),
                          (fig1_dict_Float64, fig2_dict_Float64))
        d_F["P"] = convert(Matrix{Float64}, d_R["P"])
        d_F["stationary_dists"] =
            convert(Vector{Vector{Float64}}, d_R["stationary_dists"])
    end

    testcases_Float64 = [
        Q_dict_Float64,
        Dict(
            "P" => [0.4 0.6; 0.2 0.8],
            "stationary_dists" => Vector{Float64}[[0.25, 0.75]],
            "comm_classes" => Vector{Int}[collect(1:2)],
            "rec_classes" => Vector{Int}[collect(1:2)],
            "is_irreducible" => true,
            "period" => 1,
            "is_aperiodic" => true,
            # "cyclic_classes" => Vector[collect(1:2)],
        ),
        Dict(
            "P" => [0. 1.; 1. 0.],
            "stationary_dists" => Vector{Float64}[[0.5, 0.5]],
            "comm_classes" => Vector{Int}[collect(1:2)],
            "rec_classes" => Vector{Int}[collect(1:2)],
            "is_irreducible" => true,
            "period" => 2,
            "is_aperiodic" => false,
            # "cyclic_classes" => Vector[[1], [2]],
        ),
        Dict(
            "P" => eye(2),
            "stationary_dists" => Vector{Float64}[[1, 0], [0, 1]],
            "comm_classes" => Vector{Int}[[1], [2]],
            "rec_classes" => Vector{Int}[[1], [2]],
            "is_irreducible" => false,
            "period" => 1,
            "is_aperiodic" => true,
        ),
        # Reducible mc with a unique recurrent class,
        # where n-1 is a transient state
        Dict(
            "P" => [1. 0.; 1. 0.],
            "stationary_dists" => Vector{Float64}[[1, 0]],
            "comm_classes" => Vector{Int}[[1], [2]],
            "rec_classes" => Vector{Int}[[1]],
            "is_irreducible" => false,
            "period" => 1,
            "is_aperiodic" => true,
        ),
        # these matrices come from RMT4 section 2.2.1
        Dict(
            "P" => [1 0 0; .2 .5 .3; 0 0 1],
            "stationary_dists" => Vector{Float64}[[1, 0, 0], [0, 0, 1]],
            "comm_classes" => Vector{Int}[[1], [2], [3]],
            "rec_classes" => Vector{Int}[[1], [3]],
            "is_irreducible" => false,
            "period" => 1,
            "is_aperiodic" => true,
        ),
        Dict(
            "P" => [.7 .3 0; 0 .5 .5; 0 .9 .1],
            "stationary_dists" => Vector{Float64}[[0, 9/14, 5/14]],
            "comm_classes" => Vector{Int}[[1], [2, 3]],
            "rec_classes" => Vector{Int}[[2, 3]],
            "is_irreducible" => false,
            "period" => 1,
            "is_aperiodic" => true,
        ),
        fig1_dict_Float64,
        fig2_dict_Float64,
    ]

    @testset "test MarkovChain with Float64" begin
        for test_dict in testcases_Float64
            mc = MarkovChain(test_dict["P"])
            stationary_dists = @inferred stationary_distributions(mc)
            @test isapprox(stationary_dists, test_dict["stationary_dists"])
            @test isequal(
                sort(communication_classes(mc), by=(x->x[1])),
                test_dict["comm_classes"]
                )
            @test isequal(
                sort(recurrent_classes(mc), by=(x->x[1])),
                test_dict["rec_classes"]
                )
            @test is_irreducible(mc) == test_dict["is_irreducible"]
            @test period(mc) == test_dict["period"]
            @test is_aperiodic(mc) == test_dict["is_aperiodic"]
        end
    end

    testcases_Rational = [
        Q_dict_Rational,
        Dict(
            "P" => [2//3 1//3; 1//4 3//4],
            "stationary_dists" => Vector{Rational{Int}}[[3//7, 4//7]],
            "comm_classes" => Vector{Int}[collect(1:2)],
            "rec_classes" => Vector{Int}[collect(1:2)],
            "is_irreducible" => true,
            "period" => 1,
            "is_aperiodic" => true,
            # "cyclic_classes" => Vector[collect(1:2)],
        ),
        fig1_dict_Rational,
        fig2_dict_Rational,
    ]

    @testset "test MarkovChain with Rational" begin
        for test_dict in testcases_Rational
            mc = MarkovChain(test_dict["P"])
            stationary_dists = @inferred stationary_distributions(mc)
            @test isequal(stationary_dists, test_dict["stationary_dists"])
            @test isequal(
                sort(communication_classes(mc), by=(x->x[1])),
                test_dict["comm_classes"]
                )
            @test isequal(
                sort(recurrent_classes(mc), by=(x->x[1])),
                test_dict["rec_classes"]
                )
            @test is_irreducible(mc) == test_dict["is_irreducible"]
            @test period(mc) == test_dict["period"]
            @test is_aperiodic(mc) == test_dict["is_aperiodic"]
        end
    end

    testcases_Int = [
        Dict(
            "P" => [0 1; 1 0],
            "stationary_dists" => Vector{Rational{Int}}[[1//2, 1//2]],
            "comm_classes" => Vector{Int}[collect(1:2)],
            "rec_classes" => Vector{Int}[collect(1:2)],
            "is_irreducible" => true,
            "period" => 2,
            "is_aperiodic" => false,
            # "cyclic_classes" => Vector[[1], [2]],
        ),
        Dict(
            "P" => eye(Int, 2),
            "stationary_dists" => Vector{Rational{Int}}[[1, 0], [0, 1]],
            "comm_classes" => Vector{Int}[[1], [2]],
            "rec_classes" => Vector{Int}[[1], [2]],
            "is_irreducible" => false,
            "period" => 1,
            "is_aperiodic" => true,
        ),
        # Reducible mc with a unique recurrent class,
        # where n-1 is a transient state
        Dict(
            "P" => [1 0; 1 0],
            "stationary_dists" => Vector{Rational{Int}}[[1, 0]],
            "comm_classes" => Vector{Int}[[1], [2]],
            "rec_classes" => Vector{Int}[[1]],
            "is_irreducible" => false,
            "period" => 1,
            "is_aperiodic" => true,
        ),
        Dict(
            "P" => [0 1 0; 1 0 0; 0 0 1],
            "stationary_dists" =>
                Vector{Rational{Int}}[[1//2,1//2,0//1], [0//1,0//1,1//1]],
            "comm_classes" => Vector{Int}[[1, 2], [3]],
            "rec_classes" => Vector{Int}[[1, 2], [3]],
            "is_irreducible" => false,
            "period" => 2,
            "is_aperiodic" => false,
        ),
    ]

    @testset "test MarkovChain with Int" begin
        for test_dict in testcases_Int
            mc = MarkovChain(test_dict["P"])
            stationary_dists = @inferred stationary_distributions(mc)
            @test isequal(stationary_dists, test_dict["stationary_dists"])
            @test isequal(
                sort(communication_classes(mc), by=(x->x[1])),
                test_dict["comm_classes"]
                )
            @test isequal(
                sort(recurrent_classes(mc), by=(x->x[1])),
                test_dict["rec_classes"]
                )
            @test is_irreducible(mc) == test_dict["is_irreducible"]
            @test period(mc) == test_dict["period"]
            @test is_aperiodic(mc) == test_dict["is_aperiodic"]
        end
    end

    @testset "test Sparse MarkovChain" begin
        # Test 2 testcases from each of testcases_*
        for test_dict in [testcases[i] for i in 1:2,
                                           testcases in (testcases_Float64,
                                                         testcases_Rational,
                                                         testcases_Int)]
            mc = MarkovChain(sparse(test_dict["P"]))
            stationary_dists = @inferred stationary_distributions(mc)
            @test isapprox(stationary_dists, test_dict["stationary_dists"])
            @test isequal(
                sort(communication_classes(mc), by=(x->x[1])),
                test_dict["comm_classes"]
                )
            @test isequal(
                sort(recurrent_classes(mc), by=(x->x[1])),
                test_dict["rec_classes"]
                )
            @test is_irreducible(mc) == test_dict["is_irreducible"]
            @test period(mc) == test_dict["period"]
            @test is_aperiodic(mc) == test_dict["is_aperiodic"]
        end
    end

    tol = 1e-15
    kmr_matrices = (
        kmr_markov_matrix_sequential(27, 1/3, 1e-2),
        kmr_markov_matrix_sequential(3, 1/3, 1e-14)
    )

    @testset "test gth_solve with KMR matrices" begin
        for P in kmr_matrices
            x = gth_solve(P)

            # Check elements sum to one
            @test isapprox(sum(x), 1; atol=tol)

            # Check elements are nonnegative
            for i in 1:length(x)
                @test x[i] >= -tol
            end

            # Check x is a left eigenvector of P
            @test isapprox(vec(x'*P), x; atol=tol)
        end
    end

    @testset "test gth_solve with generator matrices" begin
        P_dict_Int = Dict(
            "P" => [-3 3; 4 -4],
            "stationary_dist" => [4//7, 3//7],
        )

        P_dict_Float64 = Dict(
            "P" => [-3. 3.; 4. -4.],
            "stationary_dist" => [4/7, 3/7],
        )

        x = gth_solve(P_dict_Int["P"])
        @test isequal(x, P_dict_Int["stationary_dist"])

        x = gth_solve(P_dict_Float64["P"])
        @test isapprox(x, P_dict_Float64["stationary_dist"])
    end

    @testset "test gth_solve with reducible matrix" begin
        P_dict = Dict(
            "P" => [0.4 0 0.6 0; 0 0 0 1; 0.2 0 0.8 0; 0 1 0 0],
            # Stationary dist for the recurrent class {1, 3}
            "stationary_dist" => [0.25, 0, 0.75, 0],
        )

        x = gth_solve(P_dict["P"])
        @test isapprox(x, P_dict["stationary_dist"])
    end

    @testset "test MarkovChain with KMR matrices" begin
        for P in kmr_matrices
            mc = MarkovChain(P)
            stationary_dists = stationary_distributions(mc)
            for x in stationary_dists
                # Check elements sum to one
                @test isapprox(sum(x), 1; atol=tol)

                # Check elements are nonnegative
                for i in 1:length(x)
                    @test x[i] >= -tol
                end

                # Check x is a left eigenvector of P
                @test isapprox(vec(x'*P), x; atol=tol)
            end
        end
    end

    @testset "test MarkovChain throws errors" begin
        # not square
        @test_throws DimensionMismatch MarkovChain(rand(4, 5))

        # state_values to long
        @test_throws DimensionMismatch MarkovChain([0.5 0.5; 0.2 0.8], rand(3))

        # first row doesn't sum to 1
        @test_throws ArgumentError MarkovChain([0.0 0.5; 0.2 0.8])

        # negative element, but sums to 1
        @test_throws ArgumentError MarkovChain([-1 1; 0.2 0.8])
    end

    mc3 = MarkovChain([0.4 0.6; 0.2 0.8])

    @testset "test simulate_indices shape" begin

        for mc in (mc3, MarkovChain(sparse(mc3.p)))
            ts_length = 10
            init = 2
            nums_sims = [3, 1]

            @test size(simulate_indices(mc, ts_length)) == (ts_length,)
            @test size(simulate_indices(mc, ts_length; init=init)) ==
                       (ts_length, )
            for num_sims in nums_sims
                X = Array(Int64, ts_length, num_sims)
                @test size(simulate_indices!(X, mc)) ==
                           (ts_length, num_sims)
                @test size(simulate_indices!(X, mc; init=init)) ==
                           (ts_length, num_sims)
            end
        end
    end  # testset

    @testset "test simulate_indices init array" begin
        for mc in (mc3, MarkovChain(sparse(mc3.p)))
            ts_length = 10
            init = 2
            Y = simulate_indices(mc, ts_length; init=init)
            @test Y[1] == init

            num_sims = 5
            X = Array(Int64, ts_length, num_sims)
            simulate_indices!(X, mc; init=init)
            @test vec(X[1, :]) == init .* ones(Int, num_sims)

            init = [2, 1]
            simulate_indices!(X, mc; init=init)
            @test vec(X[1, :]) == collect(take(cycle(init), num_sims))
        end
    end  # testset

    @testset "Sparse Matrix" begin
        # use fig1_p from above because we already checked the graph theory
        # algorithms for this stochastic matrix.
        p = fig1_p
        p_s = sparse(p)

        @testset "constructors" begin
            mc_s = MarkovChain(p_s)
            @test isa(mc_s, MarkovChain{eltype(p), typeof(p_s)})
        end

        mc_s = MarkovChain(p_s)
        mc = MarkovChain(p)

        @testset "basic correspondence with dense version" begin
            @test n_states(mc) == n_states(mc_s)
            @test maxabs(mc.p - mc_s.p) == 0.0
            @test recurrent_classes(mc) == recurrent_classes(mc_s)
            @test communication_classes(mc) == communication_classes(mc_s)
            @test is_irreducible(mc) == is_irreducible(mc_s)
            @test is_aperiodic(mc) == is_aperiodic(mc_s)
            @test period(mc) == period(mc_s)
            @test isapprox(stationary_distributions(mc_s),
                           stationary_distributions(mc))
        end
    end

    @testset "simulate" begin
        for _mc in (mc3, MarkovChain(sparse(mc3.p)))
            for T in [Float16, Float32, Float64, Int8, Int16, Int32, Int64]
                mc = MarkovChain(_mc.p, rand(T, size(_mc.p, 1)))
                ts_length = 10
                init = 2
                nums_sims = [3, 1]

                @test size(@inferred(simulate(mc, ts_length))) == (ts_length,)
                for num_sims in nums_sims
                    X = Array(T, ts_length, num_sims)
                    @test size(simulate!(X, mc)) ==
                           (ts_length, num_sims)
                    @test size(simulate!(X, mc; init=init)) ==
                           (ts_length, num_sims)
                end
            end  # state_values eltypes
        end

        @testset "test simulate init array" begin
            for mc in (mc3, MarkovChain(sparse(mc3.p)))
                ts_length = 10
                init = 2
                Y = simulate(mc, ts_length; init=init)
                @test Y[1] == mc.state_values[init]

                num_sims = 5
                X = Array(eltype(mc.state_values), ts_length, num_sims)
                simulate!(X, mc; init=init)
                @test vec(X[1, :]) == mc.state_values[init .* ones(Int, num_sims)]

                init = [2, 1]
                simulate!(X, mc; init=init)
                @test vec(X[1, :]) == mc.state_values[collect(take(cycle(init), num_sims))]
            end
        end  # testset
    end

    @testset "simulate iterators" begin
        p = [0.0 1.0 0.0 0.0
             0.0 0.0 1.0 0.0
             0.0 0.0 0.0 1.0
             1.0 0.0 0.0 0.0]
        mc = MarkovChain(p, [10.0, 20.0, 30.0, 40.0])

        mcis = MCIndSimulator(mc, 50, 1)

        want = collect(take(cycle(1:4), 50))
        have = Array(Int, 50)
        for (ix, i) in enumerate(mcis)
            have[ix] = i
        end
        @test have == want

        @test start(mcis) == (1, 0)
        for i in 1:49
            @test !done(mcis, (1, i))
        end
        @test done(mcis, (1, 50))
        @test length(mcis) == 50
        @test eltype(mcis) == Int

        mcs = MCSimulator(mc, 50, 1)
        want = collect(take(cycle([10, 20, 30, 40]), 50))
        have = zeros(Float64, 50)
        for (ix, i) in enumerate(mcs)
            have[ix] = i
        end

        @test have == want

        @test start(mcs) == (1, 0)
        for i in 1:49
            @test !done(mcs, (1, i))
        end
        @test done(mcs, (1, 50))
        @test length(mcs) == 50
        @test eltype(mcs) == Float64

        if isdefined(Base, :iteratorsize)
            @test Base.iteratorsize(mcs) == Base.HasLength()
            @test Base.iteratorsize(mcis) == Base.HasLength()
        end
    end

end  # testset


text:
function julia_to_llvm(@nospecialize x)
    isboxed = Ref{UInt8}()
    # LLVMType(ccall(:jl_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Ref{UInt8}), x, isboxed))    # noserialize
    LLVMType(ccall(:jl_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Ref{UInt8}), x, isboxed))  # julia v1.1.1
end

const jl_value_t_ptr = julia_to_llvm(Any)
const jl_value_t = eltype(jl_value_t_ptr)
# const jl_value_t_ptr_ptr = LLVM.PointerType(jl_value_t_ptr)
# # cheat on these for now:
# const jl_datatype_t_ptr = jl_value_t_ptr
# const jl_unionall_t_ptr = jl_value_t_ptr
# const jl_typename_t_ptr = jl_value_t_ptr
# const jl_sym_t_ptr = jl_value_t_ptr
# const jl_svec_t_ptr = jl_value_t_ptr
# const jl_module_t_ptr = jl_value_t_ptr
# const jl_array_t_ptr = jl_value_t_ptr
#
# const bool_t  = julia_to_llvm(Bool)
# const int8_t  = julia_to_llvm(Int8)
# const int16_t = julia_to_llvm(Int16)
# const int32_t = julia_to_llvm(Int32)
# const int64_t = julia_to_llvm(Int64)
# const uint8_t  = julia_to_llvm(UInt8)
# const uint16_t = julia_to_llvm(UInt16)
# const uint32_t = julia_to_llvm(UInt32)
# const uint64_t = julia_to_llvm(UInt64)
# const float_t  = julia_to_llvm(Float32)
# const double_t = julia_to_llvm(Float64)
# const float32_t = julia_to_llvm(Float32)
# const float64_t = julia_to_llvm(Float64)
# const void_t    = julia_to_llvm(Nothing)
# const size_t    = julia_to_llvm(Int)
#
# const int8_t_ptr  = LLVM.PointerType(int8_t)
# const void_t_ptr  = LLVM.PointerType(void_t)

function module_setup(mod::LLVM.Module)
#    triple!(mod, "wasm32-unknown-unknown-wasm")
#    datalayout!(mod, "e-m:e-p:32:32-i64:64-n32:64-S128")
end

llvmmod(native_code) =
    LLVM.Module(ccall(:jl_get_llvm_module, LLVM.API.LLVMModuleRef,
                      (Ptr{Cvoid},), native_code.p))

function Base.write(mod::LLVM.Module, path::String)
    open(io -> write(io, mod), path, "w")
end

walk(f, x) = true
# walk(f, x::Instruction) = foreach(c->walk(f,c), operands(x))
# walk(f, x::Instruction) = f(x) || foreach(c->walk(f,c), operands(x))
walk(f, x::ConstantExpr) = f(x) || foreach(c->walk(f,c), operands(x))


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements rsubsolve function
# 3. Exports Type

#using LinearAlgebra

"""
    LinSysVecColProjStd <: LinSysVecColProjection

A mutable structure that represents a standard column projection method.

See Luo, Zhi-Quan, and Paul Tseng. "On the convergence of the coordinate descent method
    for convex differentiable minimization." Journal of Optimization Theory and
    Applications 72.1 (1992): 7-35.

For a generalization, see Patel, Vivak, Mohammad Jahangoshahi, and Daniel Adrian Maldonado.
    "An Implicit Representation and Iterative Solution of Randomly Sketched Linear Systems."
    SIAM Journal on Matrix Analysis and Applications (2021) 42:2, 800-831.
    https://doi.org/10.1137/19M1259481

# Aliases
- `CoordinateDescent`
- `GaussSeidel`

# Fields
- `alpha::Float64`, a relaxation parameter that should be set between `0.0` and `2.0`

Calling `LinSysVecColProjStd()` defaults the relaxatoin parameter to `1.0`.

"""
mutable struct LinSysVecColProjStd <: LinSysVecColProjection
    α::Float64
end
CoordinateDescent = LinSysVecColProjStd
GaussSeidel = LinSysVecColProjStd

LinSysVecColProjStd() = LinSysVecColProjStd(1.0)

# Common rsubsolve interface for linear systems
function rsubsolve!(
    type::LinSysVecColProjStd,
    x::AbstractVector,
    samp::Tuple{U,V,W} where {U<:AbstractVector,V<:AbstractArray,W<:Real},
    iter::Int64,
)
    # samp[1] is the search direction
    # samp[2] is the coefficient matrix A
    # samp[3] is the residual of normal system in A * samp[1], samp[1]' A' * (A * x - b)

    x .= x - samp[1] * (type.α * samp[3] / sum((samp[2] * samp[1]).^2))

    return nothing
end


text:
# included in handshaketests.jl
import Base.BufferStream
function templaterequests()
    chromeheaders = Dict{String, String}( "Connection"=>"Upgrade",
                                            "Upgrade"=>"websocket")
    firefoxheaders = Dict{String, String}("Connection"=>"keep-alive, Upgrade",
                                            "Upgrade"=>"websocket")
    chromerequest = HTTP.Request("GET", "/", collect(chromeheaders))
    firefoxrequest = HTTP.Request("GET", "/", collect(firefoxheaders))
    return [chromerequest, firefoxrequest]
end
convert(::Type{HTTP.Header}, pa::Pair{String,String}) = Pair(SubString(pa[1]), SubString(pa[2]))
sethd(r::HTTP.Request, pa::Pair) = sethd(r, convert(HTTP.Header, pa))
sethd(r::HTTP.Request, pa::HTTP.Header) = HTTP.setheader(r, pa)

takefirstline(buf::IOBuffer) = strip(split(buf |> take! |> String, "\r\n")[1])
takefirstline(buf::BufferStream) = strip(split(buf |> read |> String, "\r\n")[1])
"""
The dummy websocket don't use TCP. Close won't work, but we can manipulate the contents
using otherwise the same functions as TCP sockets.
"""
dummyws(server::Bool)  = WebSocket(BufferStream(), server)

function dummywshandler(req, dws::WebSocket{BufferStream})
    close(dws.socket)
    close(dws)
end
function handshakeresponse(request::HTTP.Request)
    buf = BufferStream()
    c = HTTP.Connection(buf)
    t = HTTP.Transaction(c)
    s = HTTP.Stream(request, t)
    WebSockets.upgrade(dummywshandler, s)
    close(buf)
    takefirstline(buf)
end
nothing


text:
struct BlackBoxOptimization <: QuerySynthesisOptimizer
    method::Symbol
    eps::Float64
    params::Dict{Symbol, Any}
    function BlackBoxOptimization(method::Symbol; eps::Float64=0.1, params...)
        new(method, eps, Dict(params))
    end
end

function query_synthesis_optimize(f::Function, optimizer::BlackBoxOptimization, data::Array{T, 2}, labels::Vector{Symbol})::Array{T, 2} where T <: Real
    lb, ub = vec.(data_boundaries(data[:, labels .!= :Lout], optimizer.eps))
    res = BlackBoxOptim.bboptimize(x -> -first(f(reshape(x, length(x), 1)));
        NumDimensions=size(data, 1), Method=optimizer.method,
        SearchRange=[(lb[i], ub[i]) for i in eachindex(lb)],
        TraceMode=:silent, optimizer.params...)
    return reshape(BlackBoxOptim.best_candidate(res), size(data, 1), 1)
end


text:
module Excel

using DataFrames
using XLSX

function read_sheet(filename, sheetname)
   #reading excel sheet in entirety
    try
        _workbook = XLSX.openxlsx(filename)
        _sheet    = _workbook[sheetname]
        _array    = _sheet[:]
        return DataFrame(_array)
    catch
        println("Cannot retrieve data.")
    end

end

end #module


text:
@testset "lines" begin
    p1 = Point([11.343, 44.495])
    @test linearc(p1, 5., 20., 60.).coordinates ≈ [[
						11.364576,
						44.537252
					],
					[
						11.370282,
						44.53554
					],
					[
						11.375724,
						44.533437
					],
					[
						11.380851,
						44.530964
					],
					[
						11.385612,
						44.528145
					],
					[
						11.389963,
						44.525006
					],
					[
						11.39386,
						44.521578
					],
					[
						11.397267,
						44.517894
					],
					[
						11.397614,
						44.51747]]
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module HJBSolver

include("types.jl")

include("solver.jl")

include("solver2d.jl")

end # module

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
#
# FineShift.jl -
#
# Direct and adjoint operators for sub-sampling separable shift of
# multi-dimensional arrays and for the discrete correlation and
# convolution.
#
module FineShift

export
    convolve!,
    convolve,
    correlate!,
    correlate,
    fineshift,
    fineshift!

using ArrayTools
using InterpolationKernels

include("impl.jl")

"""
    fineshift!(dst, src, ker, t, d=1, adj=false) -> dst

overwrites `dst` with array `src` shifted along its `d`-th dimension.  The
shift `t` is given in sampling units and may have a fractional part.  The shift
is performed by interpolating array `src` along its `d`-th dimension with
kernel function `ker`.  Arrays `src` and `dst` must have the same
non-interpolated dimensions (all dimensions but the `d`-th one).

The result of the operation amounts to doing:

    dst[i,j,k] = sum_jp src[i,jp,k])*ker(j - jp - t)

with `ker` the interpolation kernel and for all possible indices.  Indices `i`
and `k` may be multi-dimensional (including 0-dimensional).

For now, only *flat* boundary conditions are implemented.

If optional argument `adj` is true, the adjoint of the linear operator
implemented by this method is applied instead.

See also [`fineshift`](@ref), [`correlate!`](@ref),

"""
function fineshift!(dst::AbstractArray{T,N},
                    src::AbstractArray{T,N},
                    ker::Kernel{T,S},
                    t::Real,
                    d::Int = 1,
                    opt::Union{Bool,Val} = false) where {T<:AbstractFloat,N,S}
    off, wgt = shiftparams(ker, T(t))
    return correlate!(dst, src, off, wgt, d, opt)
end

"""
    fineshift([len=size(arr,d)], arr, ker, t, d=1, adj=false)

yields array `arr` shifted along its `d`-th dimension.  The shift `t` is given
in sampling units and may have a fractional part.  The shift is performed by
interpolating array `arr` along its `d`-th dimension with kernel function
`ker`.  The result is an array whose `d`-th dimension has length `len` and
whose other dimensions have the same length as those of `arr`.

Schematically, assuming `res` is the result of the interpolation and
unidimensional arrays:

    res[i] ≈ arr[i - t]    for i = 1, 2, ..., len

If optional argument `adj` is true, the adjoint of the linear operator
implemented by this method is applied instead.

For now, only *flat* boundary conditions are implemented.

See also [`fineshift!`](@ref)

"""
function fineshift(arr::AbstractArray{T,N},
                   ker::Kernel{T,S},
                   t::Real,
                   args...) where {T<:AbstractFloat,N,S}
    return fineshift!(Array{T,N}(undef, size(arr)), arr, ker, t, args...)
end

function fineshift(len::Int,
                   arr::AbstractArray{T,N},
                   ker::Kernel{T,S},
                   t::Real,
                   d::Int = 1,
                   args...) where {T<:AbstractFloat,N,S}
    dims = ntuple(i -> (i == d ? len : size(arr, i)), Val(N))
    return fineshift!(Array{T,N}(undef, dims), arr, ker, t, d, args...)
end

"""
    shiftparams(ker, t) -> off, wgt

yields the parameters for fine-shifting an array along one of its dimension by
amount `t` and by means of interpolation with kernel `ker`.  The returned
values are `off`, an integer offset, and `wgt`, a `S`-tuple of interpolation
weights with `S` the size of the support of `ker`.

To perform the fine-shifting a source vector `src`, the destination vector
`dst` should be computed as:

    dst[i] = wgt[1]*src[i-off+1] + ... + wgt[k]*src[i-off+k] +
             ... + wgt[S]*src[i-off+S]

that is the sum of `wgt[k]*src[i-off+k]` for `k = 1, ..., S` with `S` the
number of interpolation weights.

!!! note
    For out of range indices `k = i-off+k`, the value taken by `src[i-off+k]`
    depends on the chosen boundary conditions.

See the note `fast-interpolation.md` for detailed explanations.

"""
function shiftparams(ker::Kernel{T}, t::Real) where {T}
    off, wgt = InterpolationKernels.compute_offset_and_weights(ker, -T(t))
    return safe_int(-off), wgt
end

# Convert `x` to an integer of type `Int` taking care of overflows but not
# of "inexact errors".
safe_int(x::Int) = x
safe_int(x::Real) = (x ≤ typemin(Int) ? typemin(Int) :
                     x ≥ typemax(Int) ? typemax(Int) : Int(x))

"""
    convolve([len=size(arr,d)], arr, off, wgt, d=1, adj=false, opt=Val(0))

yields the discrete convolution of array `arr` by weights `wgt` with offset
`off` along the `d`-th dimension.  The result is an array whose `d`-th
dimension has length `len` and whose other dimensions have the same length as
those of `arr`.

If optional argument `adj` is true, the adjoint of the linear operator
implemented by this method is applied instead.

For now, only *flat* boundary conditions are implemented.

See also [`convolve!`](@ref)

"""
function convolve(arr::AbstractArray{T,N},
                  off::Int,
                  wgt::Tuple{Vararg{T}},
                  args...) where {T<:AbstractFloat,N,S}
    return convolve!(Array{T,N}(undef, size(arr)), arr, off, wgt, args...)
end

function convolve(len::Int,
                  arr::AbstractArray{T,N},
                  off::Int,
                  wgt::Tuple{Vararg{T}},
                  d::Int = 1,
                  args...) where {T<:AbstractFloat,N,S}
    dims = ntuple(i -> (i == d ? len : size(arr, i)), Val(N))
    return convolve!(Array{T,N}(undef, dims), arr, off, wgt, d, args...)
end

"""
    convolve!(dst, src, off, wgt, d=1, adj=false, opt=Val(0)) -> dst

overwrites the contents of `dst` with the discrete convolution of `src` by a
kernel whose coefficients are given by `wgt` and with an offset `off` along the
`d`-th dimension.  Flat boundary conditions are assumed.

The call is equivalent to computing:

    dst[i] = sum_{k ∈ 1:S} wgt[k]*src[clamp(i-off-k,1,n)]

with `S = length(wgt)` and `n = length(x)` for all `i ∈ 1:lenght(dst)`.

If optional argument `adj` is true, the adjoint of the linear operator
implemented by this method is applied instead.

See also: [`convolve`](@ref), [`correlate!`](@ref).

"""
function convolve!(dst::AbstractArray{T,N},
                   src::AbstractArray{T,N},
                   off::Int,
                   wgt::NTuple{S,T},
                   args...) where {T<:AbstractFloat,N,S}
    # We really want to compute a discrete correlation because addressing array
    # entries in increasing order is the key for vectorization.
    return correlate!(dst, src, off, reverse(wgt), args...)
end

"""
    correlate([len=size(arr,d)], arr, off, wgt, d=1, adj=false, opt=Val(0))

yields the discrete correlation of array `arr` by weights `wgt` with offset
`off` along the `d`-th dimension.  The result is an array whose `d`-th
dimension has length `len` and whose other dimensions have the same length as
those of `arr`.

If optional argument `adj` is true, the adjoint of the linear operator
implemented by this method is applied instead.

For now, only *flat* boundary conditions are implemented.

See also [`correlate!`](@ref)

"""
function correlate(arr::AbstractArray{T,N},
                   off::Int,
                   wgt::Tuple{Vararg{T}},
                   args...) where {T<:AbstractFloat,N,S}
    return correlate!(Array{T,N}(undef, size(arr)), arr, off, wgt, args...)
end

function correlate(len::Int,
                   arr::AbstractArray{T,N},
                   off::Int,
                   wgt::Tuple{Vararg{T}},
                   d::Int = 1,
                   args...) where {T<:AbstractFloat,N,S}
    dims = ntuple(i -> (i == d ? len : size(arr, i)), Val(N))
    return correlate!(Array{T,N}(undef, dims), arr, off, wgt, d, args...)
end

"""
    correlate!(dst, src, off, wgt, d=1, adj=false, opt=Val(0))

overwrite the contents of `dst` with the discrete correlation of `src` by a
kernel whose coefficients are given by `wgt` and with an offset `off` along the
`d`-th dimension.  See [`shiftparams`](@ref) for a description of these latter
arguments.  Flat boundary conditions are assumed.

The call is equivalent to computing:

    dst[i] = sum_{k ∈ 1:S} wgt[k]*src[clamp(i-off+k,1,n)]

with `S = length(wgt)` and `n = length(x)` for all `i ∈ 1:lenght(dst)`.

If optional argument `adj` is true, the adjoint of the linear operator
implemented by this method is applied instead.

See also: [`correlate`](@ref), [`fineshift!`](@ref),
[`shiftparams!`](@ref), [`correlate!`](@ref).

"""
function correlate!(dst::AbstractArray{T,N},
                    src::AbstractArray{T,N},
                    off::Int,
                    wgt::NTuple{S,T},
                    d::Int = 1,
                    adj::Bool = false,
                    opt::Int = Impl.JUMP) where {T<:AbstractFloat,N,S}
    dims, m, n = checkindices(dst, src, d)
    if adj
        fill!(dst, zero(T))
        #@assert m ≥ 1
        Impl.correlate!(Impl.Adjoint(),
                        opt, dst, src, off, wgt, d, dims, n, m)
    else
        #@assert n ≥ 1
        Impl.correlate!(Impl.Direct(),
                        opt, dst, src, off, wgt, d, dims, m, n)
    end
    return dst
end

#------------------------------------------------------------------------------
# UTILITIES

"""
    checkindices(A, B, d) -> dims, m, n

checks the indices of arguments `A` and `B` for a separable operation along
their `d`-th dimension and return `dims = size(A)`, `m = size(A,d)` and
`n = size(B,d)`.

"""
@inline function checkindices(A::AbstractArray{<:Any,N},
                              B::AbstractArray{<:Any,N},
                              d::Int) where {N}
    1 ≤ d ≤ N || throw_out_of_range_dimension_index()
    Adims = standard_size(A)
    Bdims = standard_size(B)
    @inbounds begin
        for k in 1:N
            k == d || Adims[k] == Bdims[k] ||
                throw_incompatible_dimensions()
        end
        return (Adims, Adims[d], Bdims[d])
    end
end

#------------------------------------------------------------------------------
# ERRORS

@noinline throw_non_standard_indices() =
    error("array have non standard indices")

@noinline throw_out_of_range_dimension_index() =
    error("out of range dimension index")

@noinline throw_incompatible_dimensions() =
    throw(DimensionMismatch("arrays have incompatible dimensions"))

end # module


text:
# ---
# title: 1337. The K Weakest Rows in a Matrix
# id: problem1337
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/the-k-weakest-rows-in-a-matrix/description/>
# hidden: true
# ---
#
# Given a `m * n` matrix `mat` of _ones_  (representing soldiers) and _zeros_
# (representing civilians), return the indexes of the `k` weakest rows in the
# matrix ordered from the weakest to the strongest.
#
# A row _**i**_ is weaker than row _**j**_ , if the number of soldiers in row
# _**i**_ is less than the number of soldiers in row _**j**_ , or they have the
# same number of soldiers but _**i**_ is less than _**j**_. Soldiers are
# **always** stand in the frontier of a row, that is, always _ones_  may appear
# first and then _zeros_.
#
#
#
# **Example 1:**
#
#
#
#     Input: mat =
#     [[1,1,0,0,0],
#      [1,1,1,1,0],
#      [1,0,0,0,0],
#      [1,1,0,0,0],
#      [1,1,1,1,1]],
#     k = 3
#     Output: [2,0,3]
#     Explanation:
#     The number of soldiers for each row is:
#     row 0 -> 2
#     row 1 -> 4
#     row 2 -> 1
#     row 3 -> 2
#     row 4 -> 5
#     Rows ordered from the weakest to the strongest are [2,0,3,1,4]
#
#
# **Example 2:**
#
#
#
#     Input: mat =
#     [[1,0,0,0],
#      [1,1,1,1],
#      [1,0,0,0],
#      [1,0,0,0]],
#     k = 2
#     Output: [0,2]
#     Explanation:
#     The number of soldiers for each row is:
#     row 0 -> 1
#     row 1 -> 4
#     row 2 -> 1
#     row 3 -> 1
#     Rows ordered from the weakest to the strongest are [0,2,3,1]
#
#
#
#
# **Constraints:**
#
#   * `m == mat.length`
#   * `n == mat[i].length`
#   * `2 <= n, m <= 100`
#   * `1 <= k <= m`
#   * `matrix[i][j]` is either 0 **or** 1.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using WebIO
using Test
using Sockets
import WebIO: DOM, instanceof

# Always bundle first so we're testing against latest version of WebIO JS.
include("bundling.jl")

@testset "node" begin

    @test dom"div"().children == []

    n = dom"div"("hello")
    @test n.children == ["hello"]
    @test instanceof(n) == DOM(:html, :div)

    n = dom"div"("hello", prop1=1, prop2=Dict(:x=>1, :y=>2))
    @test n.props[:prop1] == 1
    @test n.props[:prop2] == Dict(:x=>1, :y=>2)

    # test immutability
    n2 = append(n, ["world"])
    @test n.children == ["hello"]
    @test n2.children == ["hello", "world"]

    n3 = mergeprops(n, Dict(:x=>:x, :prop1=>2, :prop2=>Dict(:x=>2)))
    @test n.children == ["hello"]
    @test n3.children == ["hello"]
    @test n3.props == Dict(:x=>:x, :prop1=>2, :prop2=>Dict(:x=>2, :y=>2))

    n4 = withchild(n2, 2) do x
        @test x == "world"
        "www"
    end

    @test n4.children[2] == "www"
end

@testset "syntax for extending Nodes" begin

    n = dom"div"()
    @test children(n("hello", "world")) == ["hello", "world"]
    @test props(n(prop="x")) == Dict(:prop=>"x")

    n1 = n("hello", "world", prop="x")
    @test children(n1) == ["hello", "world"]
    @test props(n1) == Dict(:prop=>"x")
end

function capture_plaintext(x)
    io = IOBuffer()
    show(io, MIME("text/plain"), x)
    seek(io, 0)
    read(io, String)
end

@testset "plaintext printing" begin
    @test capture_plaintext(node(:div)) == "(div)"
    @test capture_plaintext(node(:div, node(:ul))) == "(div\n  (ul))"
    @test capture_plaintext(node(:div, node(:ul, node(:li, "hello")))) == "(div\n  (ul\n    (li\n      \"hello\")))"
    @test capture_plaintext(node(:div, id="foo", class="bar")) == "(div { id=\"foo\" class=\"bar\" })"
    @test capture_plaintext(dom"ul"(dom"li"("hello"), dom"li"("world"))) == "(ul\n  (li\n    \"hello\")\n  (li\n    \"world\"))"
end

@testset "@js_str" begin
    @test js"x=y".s == "x=y"
    y = 1
    @test js"x=$y".s == "x=1"

    y = js"f(42)"
    @test js"x=$y" == js"x=f(42)"

    y = Dict()
    @test js"x=$y" == js"x={}"
end

include("communication.jl")
include("node.jl")
include("http-tests.jl")
include("mux-tests.jl")
include("blink-tests.jl")
include("util-tests.jl")
include("ijulia-tests.jl")
include("syntax.jl")
include("deprecations.jl")


text:
abstract type AbstractExecution end
"""
    RawExecutionParameter

An optimizer attribute used for raw parameters of the execution. Defers to `RawParameter`.
"""
struct RawExecutionParameter <: ExecutionParameter
    name::Any
end

function MOI.get(executor::AbstractExecution, param::RawExecutionParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(executor)))
        error("Unrecognized parameter name: $(name) for executor $(typeof(executor)).")
    end
    return getfield(executor, name)
end

function MOI.set(executor::AbstractExecution, param::RawExecutionParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(executor)))
        error("Unrecognized parameter name: $(name) for executor $(typeof(executor)).")
    end
    setfield!(executor, name, value)
    return nothing
end
"""
    get_execution_attribute(stochasticprogram::StochasticProgram, name::String)

Return the value associated with the execution-specific attribute named `name` in `stochasticprogram`.

See also: [`set_execution_attribute`](@ref), [`set_execution_attributes`](@ref).
"""
function get_execution_attribute(stochasticprogram::StochasticProgram, name::String)
    return return MOI.get(optimizer(stochasticprogram), RawExecutionParameter(name))
end
"""
    set_execution_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)

Sets the execution-specific parameter identified by `name` to `value`.

"""
function set_execution_attribute(stochasticprogram::StochasticProgram, name::Union{Symbol, String}, value)
    return set_optimizer_attribute(stochasticprogram, RawExecutionParameter(name), value)
end
"""
    set_execution_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)

Given a list of `attribute => value` pairs or a collection of keyword arguments, calls
`set_execution_attribute(stochasticprogram, attribute, value)` for each pair.

"""
function set_execution_attributes(stochasticprogram::StochasticProgram, pairs::Pair...)
    for (name, value) in pairs
        set_execution_attribute(stochasticprogram, name, value)
    end
end
function set_execution_attributes(stochasticprogram::StochasticProgram; kw...)
    for (name, value) in kw
        set_execution_attribute(stochasticprogram, name, value)
    end
end

"""
    Serial

Factory object for [`LShaped.SerialExecution`](@ref)/[`ProgressiveHedging.SerialExecution`](@ref). Pass through the `Execution` attribute.

"""
struct Serial <: AbstractExecution end

function str(::Serial)
    return ""
end

"""
    Synchronous

Factory object for [`LShaped.SynchronousExecution`](@ref)/[`ProgressiveHedging.SynchronousExecution`](@ref). Pass through the `Execution` attribute.

"""
struct Synchronous <: AbstractExecution end

function str(::Synchronous)
    return "Synchronous "
end

"""
    Asynchronous

Factory object for [`LShaped.AsynchronousExecution`](@ref)/[`ProgressiveHedging.AsynchronousExecution`](@ref). Pass through the `Execution` attribute.

...
# Parameters
- `max_active::Int = 3`: Maximum number of active iterations that run asynchronously.
- `κ::T = 0.5`: Relative amount of finished subproblems required to start a new iterate. Governs the amount of asynchronicity.
...
"""
mutable struct Asynchronous <: AbstractExecution
    max_active::Int
    κ::Float64

    function Asynchronous(; max_active::Int = 3, κ::AbstractFloat = 0.5)
        return new(max_active, κ)
    end
end

function str(::Asynchronous)
    return "Asynchronous "
end


text:
# Below, there are two ways to define a function with the same name as the one defined in
# another package (e.g., MaxwellFDFD): extending the function defined in the other package,
# or calling the function defined in the other package by qualifying it with the package
# name within the definition of the function in the present package.
#
# For example, when defining a function that delegates to a function fun_XXX!() in
# MaxwellFDFD, we would want to use the same function name.  This can be done as either
#     MaxwellFDFD.fun_XXX!(mdl::Model, ...) = fun_XXX!(mdl.mx, ...)
# or
#     fun_XXX!(mdl::Model, ...) = MaxwellFDFD.fun_XXX!(mdl.mx, ...)
#
# Both definitions works.  If we reexport MaxwellFDFD from the present package, we do not
# even need to export the function, because the function name is already exported by
# MaxwellFDFD.  (If MaxwellFDFD were not reexported, we would need to export fun_XXX!
# explicitly in this file.)
#
# Which one is the right approach?
#
# If we are extending a function that is already exported by the present package (e.g.,
# using reexport), I think the former is better, because it explicitly specifies that we are
# extending the capabilities of an already exported function.
#
# On the other hand, if we are "extending" a function that is exported by another package
# (e.g., SALTBase) that is not reexported by the present package, then we need to export the
# function anyway to make it available to the users.  Then, we are actually not extending
# the capability of some function that is already available to the users.  Therefore, in
# that case I think using the latter makes more sense.
#
# Currently, I am reexporting MaxwellFDFD from the present package, but not SALTBase.  (The
# users may want to solve a driven Maxwell's equations, so I reexport MaxwellFDFD, but they
# do not need to know the internals of SALTBase, so I do not reexport SALTBase.)  Therefore,
# when redefining functions already defined in MaxwellFDFD and SALTBase below, I extend
# MaxwellFDFD's functions (the former approach), but call SALTBase's functions with
# qualification (the latter approach).

export GainProfile
export get_εcold, get_dblcurl, add_gainobj!
export init_salt!, simulate!

# - Do not export Model; quality it with the package name MaxwellSALT, because I would have
# similar types in other packages such as MaxwellFDFD and MaxwellGuide.
# - It might be possible to define mx for other Maxwell solvers (e.g., COMSOL).  Then MX
# doesn't have to be <:MaxwellFDFD.Model.  Let's consider the possibility of this
# generalization later.
mutable struct Model{K,Kₑ,VF,VC,GP<:GainProfile{VF},MX<:MaxwellFDFD.Model{K,Kₑ},MXS<:AbstractMaxwellSolver}
    # Model components
    mx::MX  # MaxwellFDFD's model
    gp::GP  # gain profile
    gobj_vec::Vector{GainObject{K,Kₑ}}  # gain objects

    # Solver
    mxsolver::MXS

    # Initial guess
    ωguess::VecComplexF
    Ψguess::Vector{VC}

    # Solution placeholders
    nlsol::NonlasingSol{VC}
    nlvar::NonlasingVar{MXS,VC}
    lsol::LasingSol{VF,VC}
    lvar::LasingVar{MXS,VF,VC}

    function Model{K,Kₑ,VF,VC,GP,MX}(mx::MaxwellFDFD.Model, gp::GainProfile) where {
            K,Kₑ,VF,VC,GP<:GainProfile{VF},MX<:MaxwellFDFD.Model{K,Kₑ}}
        Ps = create_paramops(mx)
        Cs = create_curls(mx)
        CC = create_A(EE, 0.0, Ps, Cs)  # ω = 0.0

        Me2c, _ = create_Mcs(mx)

        # Construct Mc2e.
        grid = mx.grid
        boundft = mx.boundft
        e⁻ⁱᵏᴸ = create_e⁻ⁱᵏᴸ(mx)

        N = grid.N
        isbloch = grid.isbloch

        # Arguments of create_mean:
        # - isfwd is set to average the fields at the voxel corners to get the fields at the
        # solution field locations (but with polarization not meant for that locations).
        # (For boundft[w] = ft, voxel-corner field is forward-averaged along w.)
        # - ∆l and ∆l′⁻¹ are not supplied; use unweighted arithmetic averaging.
        isfwdₑ = boundft.==EE
        Mₑ = create_mean(isfwdₑ, N, isbloch, e⁻ⁱᵏᴸ; mx.order_cmpfirst)

        permute = SVec(ntuple(identity,Val(Kₑ)))
        Mc2e = Vector{typeof(CC)}(undef, Kₑ-1)
        for k = 1:Kₑ-1
            permute = circshift(permute, 1)  # currently return MVector; see StaticArrays.jl/issues/745
            P = create_πcmp(N, permute; mx.order_cmpfirst)
            Mc2e[k] = Mₑ * P
        end

        mxsolver = DirectMaxwellSolver(CC, Me2c, tuple(Mc2e...))

        gobj_vec = GainObject{K,Kₑ}[]

        return new{K,Kₑ,VF,VC,GP,MX,typeof(mxsolver)}(mx, gp, gobj_vec, mxsolver)
    end
end

function Model(mx::MX, gp::GP; Vtype::Type=Vector) where {K,Kₑ,GP<:GainProfile,MX<:MaxwellFDFD.Model{K,Kₑ}}
    return Model{K,Kₑ,Vtype{Float},Vtype{ComplexF},GP,MX}(mx, gp)
end

# Operators needed for calculating the initial guess solution
get_dblcurl(mdl::Model) = mdl.mxsolver.CC

function get_εcold(mdl::Model)
    Pε, _ = create_paramops(mdl.mx)

    return diag(Pε)
end

function init_salt!(mdl::Model{K,Kₑ}, ωguess::AbsVecNumber, Ψguess::AbsMatComplexF) where {K,Kₑ}
    mxsolver = mdl.mxsolver
    gp = mdl.gp

    grid = mdl.mx.grid
    N = Kₑ * prod(grid.N)
    M = length(ωguess)

    mdl.nlsol = NonlasingSol(ωguess, Ψguess)
    mdl.nlvar = NonlasingVar(mxsolver, N, M)
    mdl.lsol = LasingSol(N, M)
    mdl.lvar = LasingVar(mxsolver, N, length(gp), M)

    return nothing
end

#= Gain objects =#
add_gainobj!(mdl::Model, shapes::Shape...; kwargs...) = add_gainobj!(mdl.gobj_vec, shapes...; kwargs...)
add_gainobj!(mdl::Model, shapes::AbsVec{<:Shape}; kwargs...) = add_gainobj!(mdl.gobj_vec, shapes; kwargs...)

function simulate!(mdl::Model{K,Kₑ,VF},
                   dvec::AbsVecReal;  # trajectory of pump strength parameter to follow
                   outωaψ::NTuple{3,Bool}=(true,true,false),  # true to output ω, a, ψ
                   doutvec::AbsVecReal=dvec,  # output results when dvec[i] ∈ doutvec
                   kwargs...  # keyword arguments of SALTBase > simulation.jl/simulate!()
                   ) where {K,Kₑ,VF}
    mx = mdl.mx
    grid, boundft, ise˔shp, order_cmpfirst = getproperty.(Ref(mx), (:grid, :boundft, :ise˔shp, :order_cmpfirst))

    gt₀ = ft2gt.(EE, boundft)
    setD₀!(D₀::VF, d::Real) = assign_pumpstr!(D₀, mdl.gobj_vec, d, gt₀, grid.l, ise˔shp, order_cmpfirst)

    gp, lsol, lvar, nlsol, nlvar = getproperty.(Ref(mdl), (:gp, :lsol, :lvar, :nlsol, :nlvar))

    # Return, as functions of d, the followings:
    # eigenfrequencies, modal amplitudes, normalized mode profiles, numbers of AA steps, times taken for convergence of AA
    ωout, aout, ψout, nAA, tAA = SALTBase.simulate!(lsol, lvar, nlsol, nlvar, gp, get_εcold(mdl),
                                                    dvec, setD₀!; outωaψ, doutvec, kwargs...)

    return ωout, aout, ψout, nAA, tAA
end


text:
using SparseWrappers
using Test
using Random
using LinearAlgebra
using SparseArrays

import SparseWrappers:  Conjugate

@testset "linalg"     begin include("linalg.jl") end
@testset "sparsify"   begin include("sparsify.jl") end
@testset "views"      begin include("views.jl") end
@testset "fallback"   begin include("fallback.jl") end
@testset "combine"    begin include("combine.jl") end
nothing


text:
function initialize!(integrator, cache::DP5ThreadedCache)
  integrator.kshortsize = 4
  resize!(integrator.k, integrator.kshortsize)
  integrator.k .= [cache.update,cache.bspl,cache.dense_tmp3,cache.dense_tmp4]
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::DP5ThreadedCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache.tab
  dp5threaded_loop1(dt,tmp,uprev,a21,k1,uidx)
  f(k2, tmp, p, t+c1*dt)
  dp5threaded_loop2(dt,tmp,uprev,a31,k1,a32,k2,uidx)
  f(k3, tmp, p, t+c2*dt)
  dp5threaded_loop3(dt,tmp,uprev,a41,k1,a42,k2,a43,k3,uidx)
  f(k4, tmp, p, t+c3*dt)
  dp5threaded_loop4(dt,tmp,uprev,a51,k1,a52,k2,a53,k3,a54,k4,uidx)
  f(k5, tmp, p, t+c4*dt)
  dp5threaded_loop5(dt,tmp,uprev,a61,k1,a62,k2,a63,k3,a64,k4,a65,k5,uidx)
  f(k6, tmp, p, t+dt)
  dp5threaded_loop6(dt,u,uprev,a71,k1,a73,k3,a74,k4,a75,k5,a76,k6,update,uidx)
  f(integrator.fsallast,u,p,t+dt)
  integrator.destats.nf += 6
  if integrator.opts.adaptive
    dp5threaded_adaptiveloop(dt,utilde,btilde1,k1,btilde3,k3,btilde4,k4,btilde5,k5,btilde6,k6,btilde7,k7,uidx)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  dp5threaded_denseloop(bspl,update,k1,k3,k4,k5,k6,k7,integrator.k,d1,d3,d4,d5,d6,d7,uidx)
end

@noinline @muladd function dp5threaded_denseloop(bspl,update,k1,k3,k4,k5,k6,k7,k,d1,d3,d4,d5,d6,d7,uidx)
  Threads.@threads for i in uidx
    bspl[i] = k1[i] - update[i]
    k[3][i] = update[i] - k7[i] - bspl[i]
    k[4][i] = d1*k1[i]+d3*k3[i]+d4*k4[i]+d5*k5[i]+d6*k6[i]+d7*k7[i]
  end
end

@noinline @muladd function dp5threaded_loop1(dt,tmp,uprev,a21,k1,uidx)
  a = dt*a21
  Threads.@threads for i in uidx
    tmp[i] = uprev[i]+a*k1[i]
  end
end

@noinline @muladd function dp5threaded_loop2(dt,tmp,uprev,a31,k1,a32,k2,uidx)
  Threads.@threads for i in uidx
    tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
  end
end

@noinline @muladd function dp5threaded_loop3(dt,tmp,uprev,a41,k1,a42,k2,a43,k3,uidx)
  Threads.@threads for i in uidx
    tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
  end
end

@noinline @muladd function dp5threaded_loop4(dt,tmp,uprev,a51,k1,a52,k2,a53,k3,a54,k4,uidx)
  Threads.@threads for i in uidx
    tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
  end
end

@noinline @muladd function dp5threaded_loop5(dt,tmp,uprev,a61,k1,a62,k2,a63,k3,a64,k4,a65,k5,uidx)
  Threads.@threads for i in uidx
    tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
end

@noinline @muladd function dp5threaded_loop6(dt,u,uprev,a71,k1,a73,k3,a74,k4,a75,k5,a76,k6,update,uidx)
  Threads.@threads for i in uidx
    update[i] = a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]
    u[i] = uprev[i]+dt*update[i]
  end
end

@noinline @muladd function dp5threaded_adaptiveloop(dt,utilde,btilde1,k1,btilde3,k3,btilde4,k4,btilde5,k5,btilde6,k6,btilde7,k7,uidx)
  Threads.@threads for i in uidx
    utilde[i] = dt*(btilde1*k1[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i])
  end
end


text:
push!(LOAD_PATH, homedir()*"/Projects/QuIPS/src")

using QuIPS

quip = [
    (:H, 1),
    (:CX, (1, 2)),
    (:MEASURE, 1),
    (:MEASURE, 2)
]

QC = QCircuit(quip, 2)

for i = 1:5
    run!(QC)
    println(QC.out)
    reset!(QC)
end

println()

N = 1000
gs = 0
for i = 1:N
    run!(QC, false)
    if QC.out[1] == 1
        global gs += 1
    end
    reset!(QC)
end

println("% of correlated outcomes = 1 for $N runs: ", gs/N)

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
using Combinatorics

function run_intcode(program::Array{Int,1}, phase::Int, input::Channel, output::Channel)
    pointer = 1
    ret = nothing
    while pointer <= length(program)
        # println(program)
        instruction = reverse(digits(program[pointer], pad=5))
        # println(instruction)
        opcode = instruction[end] + instruction[end-1]*10
        modes = reverse(instruction[1:end-2])
        # println(join([opcode,": ",modes]))
        if opcode == 1 # add
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            val2 = modes[2] == 0 ? program[pt2+1] : pt2
            program[pt3+1] = val1 + val2
            pointer = pointer + 4
        elseif opcode == 2 # multiply
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            val2 = modes[2] == 0 ? program[pt2+1] : pt2
            # println(join([modes, pt1, pt2, val1, val2, val1*val2], ","))
            program[pt3+1] = val1 * val2
            pointer = pointer + 4
        elseif opcode == 3 # store input
            pt1 = program[pointer+1]
            if phase >= 0
                program[pt1+1] = phase
                phase = -1
            else
                program[pt1+1] = take!(input)
            end
            pointer = pointer + 2
        elseif opcode == 4 # output
            pt1 = program[pointer+1]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            pointer = pointer + 2
            ret = val1
            put!(output, val1)
        elseif opcode == 5 # jump if true
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            val2 = modes[2] == 0 ? program[pt2+1] : pt2
            if val1 != 0
                pointer = val2 + 1
            else
                pointer = pointer + 3
            end
        elseif opcode == 6 # jump if false
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            val2 = modes[2] == 0 ? program[pt2+1] : pt2
            if val1 == 0
                pointer = val2 + 1
            else
                pointer = pointer + 3
            end
        elseif opcode == 7 # less than
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            val2 = modes[2] == 0 ? program[pt2+1] : pt2
            program[pt3+1] = val1 < val2 ? 1 : 0
            pointer = pointer + 4
        elseif opcode == 8 # equal
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = modes[1] == 0 ? program[pt1+1] : pt1
            val2 = modes[2] == 0 ? program[pt2+1] : pt2
            program[pt3+1] = val1 == val2 ? 1 : 0
            pointer = pointer + 4
        elseif opcode == 99
            return ret
        else
            throw(ErrorException("Invalid opcode: " * string(opcode)))
        end
    end
end

function find_max_signal(program::Array{Int,1}, phases::Array{Int,1})::Int
    max = 0
    settings = []
    for phase in permutations(phases)
        a_in = Channel(1)
        b_in = Channel(1)
        c_in = Channel(1)
        d_in = Channel(1)
        e_in = Channel(1)
        e = 0
        @sync begin
            @async run_intcode(deepcopy(program), phase[1], a_in, b_in)
            @async run_intcode(deepcopy(program), phase[2], b_in, c_in)
            @async run_intcode(deepcopy(program), phase[3], c_in, d_in)
            @async run_intcode(deepcopy(program), phase[4], d_in, e_in)
            @async e = run_intcode(deepcopy(program), phase[5], e_in, a_in)
            put!(a_in, 0)
        end
        if e > max
            max = e
            settings = phase
        end
    end
    println(settings)
    println(max)
    return max
end

program = open("thrusters.csv") do f
    parse.(Int, split(readlines(f)[1], ","))
end
# program = [3,15,3,16,1002,16,10,16,1,16,15,15,4,15,99,0,0]
# phases = [4,3,2,1,0]
# program = [3,23,3,24,1002,24,10,24,1002,23,-1,23,101,5,23,23,1,24,23,23,4,23,99,0,0]
# phases = [0,1,2,3,4]
# program = [3,31,3,32,1002,32,10,32,1001,31,-2,31,1007,31,0,33,1002,33,7,33,1,33,31,31,1,32,31,31,4,31,99,0,0,0]
# phases = [1,0,4,3,2]
# find_max_signal(program, [0,1,2,3,4])
# program = [3,26,1001,26,-4,26,3,27,1002,27,2,27,1,27,26, 27,4,27,1001,28,-1,28,1005,28,6,99,0,0,5]
# program = [3,52,1001,52,-5,52,3,53,1,52,56,54,1007,54,5,55,1005,55,26,1001,54, -5,54,1105,1,12,1,53,54,53,1008,54,0,55,1001,55,1,55,2,53,55,53,4, 53,1001,56,-1,56,1005,56,6,99,0,0,0,0,10]
find_max_signal(program, [5,6,7,8,9])


text:
module TestShapeDiagrams
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra.ShapeDiagrams

A, B, C = Ob(FreeCategory, :A, :B, :C)

# Spans.
f, g = Hom(:f, C, A), Hom(:g, C, B)
span = Span(f,g)
@test apex(span) == C
@test left(span) == f
@test right(span) == g
f = Hom(:f, A, B)
@test_throws Exception Span(f,g)

# Cospans.
f, g = Hom(:f, A, C), Hom(:g, B, C)
cospan = Cospan(f,g)
@test base(cospan) == C
@test left(cospan) == f
@test right(cospan) == g

end


text:
export logrange

# FIXME: how about an easy, *performant*, lazy version ?
logrange(itr) = collect(exp10(x) for x in itr)


text:
#using Compat.Test
#using Compat

using ImmersedLayers
using Literate
using Test
##using TestSetExtensions

const GROUP = get(ENV, "GROUP", "All")

notebookdir = "../examples"
docdir = "../docs/src/manual"
litdir = "./literate"

if GROUP == "All" || GROUP == "Auxiliary"
  include("tools.jl")
  include("surface_ops.jl")
end

if GROUP == "All" || GROUP == "Literate"
  for (root, dirs, files) in walkdir(litdir)
    for file in files
      endswith(file,".jl") && @testset "$file" begin include(joinpath(root,file)) end
    end
  end
end

#=
if GROUP == "Notebooks"
  for (root, dirs, files) in walkdir(litdir)
    for file in files
      #endswith(file,".jl") && startswith(file,"3") && Literate.notebook(joinpath(root, file),notebookdir)
      endswith(file,".jl") && Literate.notebook(joinpath(root, file),notebookdir)
    end
  end
end
=#

if GROUP == "Documentation"
  for (root, dirs, files) in walkdir(litdir)
    for file in files
      endswith(file,".jl") && Literate.markdown(joinpath(root, file),docdir)
    end
  end
end


text:
@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::PointMass) = PointMass(mean(m_out) - mean(m_in1))

@rule typeof(+)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass) = NormalMeanVariance(mean(m_out) - mean(m_in1), var(m_out))
@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::UnivariateNormalDistributionsFamily) = NormalMeanVariance(mean(m_out) - mean(m_in1), var(m_in1))

@rule typeof(+)(:in2, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in1::PointMass) = MvNormalMeanCovariance(mean(m_out) - mean(m_in1), cov(m_out))
@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::MultivariateNormalDistributionsFamily) = MvNormalMeanCovariance(mean(m_out) - mean(m_in1), cov(m_in1))

# TODO: Do we need it?
# Specific case for precision parametrisation
@rule typeof(+)(:in2, Marginalisation) (m_out::NormalMeanPrecision, m_in1::PointMass) = NormalMeanPrecision(mean(m_out) - mean(m_in1), precision(m_out))
@rule typeof(+)(:in2, Marginalisation) (m_out::PointMass, m_in1::NormalMeanPrecision) = NormalMeanPrecision(mean(m_out) - mean(m_in1), precision(m_in1))

@rule typeof(+)(:in2, Marginalisation) (m_out::NormalMeanPrecision, m_in1::NormalMeanPrecision) = begin
    mout, pout = mean_precision(m_out)
    min1, pin1 = mean_precision(m_in1)
    return NormalMeanPrecision(mout - min1, pout * pin1 / (pout + pin1))
end

@rule typeof(+)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily) = begin
    mout, vout = mean_var(m_out)
    min1, vin1 = mean_var(m_in1)
    return NormalMeanVariance(mout - min1, vout + vin1)
end

@rule typeof(+)(:in2, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in1::MultivariateNormalDistributionsFamily) = begin
    mout, vout = mean_cov(m_out)
    min1, vin1 = mean_cov(m_in1)
    return MvNormalMeanCovariance(mout - min1, vout + vin1)
end


text:
function minmeanmedianmax(vec)
    (minimum(vec),round(mean(vec),2),median(vec),maximum(vec))
end
function statcols(data)
    num_col = size(data,2)
    for ind in 1:num_col
        tmp = data[:,ind]
        print("min mean median max for $ind class: ")
        println(minmeanmedianmax(tmp))
    end
end

function model_flow()
    train_path = joinpath(data_dir, "tr.jld")
    val_path   = joinpath(data_dir, "val.jld")
    test_path  = joinpath(data_dir, "te.jld")

    train = JLD.load(train_path, "tr")
    val   = JLD.load(val_path, "val")
    test  = JLD.load(test_path, "te")

    train = map(x->parse(Float64,x), train)
    val = map(x->parse(Float64,x), val)
    test = map(x->parse(Float64,x), test)
    @show size(train)
    @show size(val)
    @show size(test)

    #println("minmeanmedianmax(sum(vcat(train,val,test).== 1.0,2)):")
    #println(minmeanmedianmax(sum(vcat(train,val,test)[:,1:end-1].== 1.0,2)))
    #statcols(train)
    sleep(1)
    info("training model using Gradient Boosting decision tree")
    num_iter = 5000
    bst = model_xgboost(train,val,test,num_iter)
end

function sparse_stat()
    train_path = joinpath(data_dir, "train.jld")
    val_path   = joinpath(data_dir, "val.jld")
    test_path  = joinpath(data_dir, "test.jld")
    train = JLD.load(train_path, "train")
    val   = JLD.load(val_path, "val")
    test  = JLD.load(test_path, "test")

    data = vcat(train,val,test)
end

function apply_model(bst,vector::Array{Float32,1})
    test = DMatrix(vector')
    label = XGBoost.predict(bst,test)
    label
end
function apply_model(mdl::ASCIIString,vector::Array{Float32,1})
    bst = Booster(model_file = mdl)
    test = DMatrix(vector')
    label = XGBoost.predict(bst,test)
    label
end

@doc """ modeling data with Gradient Boosting decesion tree
""" ->
function model_xgboost(data)
    model_xgboost(data[1],data[2],data[3])
end

@doc """ xgboost for train,val,test dataset without gridtune
         eval with weighted F1-score.
""" ->
function model_xgboost(train, val, test, num_iter::Int64)
    train_X,train_Y = train[:,1:end-1],train[:,end]
    val_X,val_Y = val[:,1:end-1],val[:,end]
    test_X,test_Y = test[:,1:end-1],test[:,end]

    num_class = maximum(train_Y)+1
    m_tr,n_tr   = size(train_X)
    m_val,n_val = size(val_X)
    info("in the model, there are $num_class classes")
    info("train_X size is :$m_tr,$n_tr")
    info("val_X size is : $m_val, $n_tr")
    #numcancersample_data(train,cancer_names())
    #sleep(3)

    num_round = num_iter

    dtrain = DMatrix(train_X, label = train_Y)
    dval   = DMatrix(val_X,   label = val_Y)
    dtest  = DMatrix(test_X,  label = test_Y)

    watch_list = [(dtrain,"train"), (dval,"val")]
    param      = Dict{ASCIIString,Any}("max_depth"=>7,"eta"=>0.01,"nthread"=>50,
                                       "objective"=>"multi:softprob","silent"=>1,
                                       #"alpha"=>0.7,

                                       "sub_sample"=>0.9) #,"num_class"=>num_class)

    bst = xgboost(dtrain,num_round,watchlist=watch_list,param=param,num_class=num_class,
                  metrics=["mlogloss"], seed=2015)

    #XGBoost.save(bst, cancer33_mdl)
    test_preds = XGBoost.predict(bst, test_X)

    #info("on test dataset mean-f1-score is $mean_f1")
    yprob = reshape(test_preds, Int64(num_class), size(test_X,1))
    ind_prob = map(ind->indmax(yprob[:,ind]), 1:size(test_X,1))
    test_preds = ind_prob - 1 # for index start different

    #test_preds = map(x->convert(Int64,x), test_preds)
    #test_Y = map(x->convert(Int64,x), test_Y)
    #num_class = convert(Int64, num_class)
    #ROC_print(test_Y.+1,test_preds.+1,num_class = num_class, challenge_list=cancer_names())
    #numcancersample_data(train,cancer_names())
    return bst
end

### Model Tuning useful in MLBase
function estfun(max_depth,eta)
    param["max_depth"] = max_depth
    param["eta"] = eta

    if CV
        bst = nfold_cv(dtrain,num_round,nfold,param=param,feval=evaluerror_softmax,show_stdv=false)
    else
        bst = xgboost(dtrain,num_round,watchlist=watchlist,param=param,feval=evaluerror_softmax,seed=0, early_stopping_rounds=80)
    end
    return bst
end

function evalfun(bst)
    test_preds = XGBoost.predict(bst,test_X)
    mean_f1 = evaluerror_softmax(test_preds,dtest)
    return mean_f1
end

function gridtune()
    if debug
        r = gridtune(estfun,evalfun,("max_depth",[2,10]),("eta",[0.1]))
    else
        r = gridtune(estfun,evalfun,("max_depth",[2,4,6,10,30,50,100]),("eta",[0.01,0.05,0.1,0.3,0.5,0.8,1]))
    end

    best_model, best_cfg, best_score = r
    @printf "===============================================================================================\n"
    @printf "softmax_best %10s  is %8f \n" best_score[1] best_score[2]
    @printf "===============================================================================================\n"
    f = open("softmax_best_scores.txt","a")
    @printf f "softmax_best %10s is %8f \n"  best_score[1] best_score[2]
    close(f)

    #save(best_model,"../models/softmax.model")

    #bst = Booster(model_file = "../models/softmax.model")
    #p = XGBoost.predict(bst,test[:,1:end-1]) .+ 1

    p = XGBoost.predict(best_model,test[:,1:end-1]) .+ 1

    test_Y_gt,train_pre = convert(Array{Int64},test[:,end]),convert(Array{Int64},p)

    ROC_print(test_Y_gt,train_pre,num_class=num_class)
    #   if WIN
    p = XGBoost.predict(best_model,test_I[:,1:end-1]) .+ 1

    test_Y_gt,train_pre = convert(Array{Int64},test_I[:,end]),convert(Array{Int64},p)

    ROC_print(test_Y_gt.+1,train_pre,num_class=num_class)
    #  end

end

@doc """eval ROC_value of model
"""->
function ROC_print(gt,pred;num_class=2,challenge_list=cancer_list)

    C = confusmat(num_class,gt,pred)
    println(C)

    output = Array{Pair{ASCIIString,Float64},1}(num_class)

    for i=1:num_class
        gt_tmp = (gt .== i)
        pred_tmp = (pred .== i)
        r_tmp = MLBase.roc(gt_tmp,pred_tmp)
        f1_tmp = MLBase.f1score(r_tmp)
        output[i] = Pair{ASCIIString,Float64}(challenge_list[i],f1_tmp)

    end
    output = sort(output,by=x->x[2],rev=true)
    for i=1:num_class
        @printf "%35s\t %8f\n" output[i][1] output[i][2]
    end
end

@doc """ shuffle dataset
         random indexes of row
""" ->
function shuffle(data::SparseMatrixCSC{Float64,Int64})
    m,n = size(data)
    inds = randperm(m)
    data[inds,:]
end

function split_cancer_data(data::SparseMatrixCSC{Float64,Int64})
    m,n = size(data)
    point1 = round(Int64, 0.6*m)
    point2 = round(Int64, 0.8*m)
    train = data[1:point1,:]
    val   = data[point1+1:point2,:]
    test  = data[point2+1:end,:]

    train,val,test
end

@doc """ split each cancer dataset
""" ->
function split_data(data::SparseMatrixCSC{Float64,Int64})
    num_label = maximum(data[:,end])

    function reducer(data1, data2)
        (vcat(data1[1],data2[1]), vcat(data1[2], data2[2]), vcat(data1[3], data2[3]))
    end

    data = @parallel (reducer) for i = 1:num_label
        split_cancer_data(data[data[:,end] .== i,:])
    end

    data
end


text:
import Flux: softmax, σ
using Random
using Zygote

rng = MersenneTwister(234)

xs = sort(Float32.([-1000, -10.0, -1.5, 0.0, 0.1, 1.0, 15.0]))
@testset "oneplus" for i in 1:length(xs)
    # Outputs correct range
    @test typeof(DNC.oneplus(xs[i])) == Float32
    @test DNC.oneplus(xs[i]) >= 1.0
    # Maintains order
    if i>1
        @test DNC.oneplus(xs[i]) > DNC.oneplus(xs[i-1])
    end
end

@testset "clip" begin
    arr = [-3, -2.5, -1, 1.5, 3.5]
    @test DNC.clip(arr, nothing) == arr
    @test DNC.clip(arr, 2) == [-2, -2, -1, 1.5, 2]
end

@testset "Calc output" begin
    outsize = 5
    N, W, R, B = 3, 5, 2, 1
    readvectors = rand(rng, Float32, W*R, B)
    Wr = rand(rng, Float32, outsize, R*W, B)
    v = rand(rng, Float32, outsize, B)
    res = DNC.calcoutput(v, readvectors, Wr)
    @test size(res) == (outsize, B)
    @test eltype(res) == Float32
end

@testset "Split ξ $(R) read head" for R in 1:2
    W, B = 5, 1
    ξ = rand(rng, Float32, 10, B)
    transforms = DNC.inputmappings(10, R, W)
    inputs = DNC.split_ξ(ξ, transforms)

    @testset "Type Float32" for inp in inputs
        eltype(inp) == Float32
    end

    @testset "Dimensions" begin
        @test size(inputs.kr) == (W, R, B)
        @test size(inputs.kw) == (W, 1, B)
        @test size(inputs.βr) == (R,B)
        @test size(inputs.βw) == (1,B)
        @test size(inputs.ga) == (1,B)
        @test size(inputs.gw) == (1,B)
        @test size(inputs.v) == (W,B)
        @test size(inputs.e) == (W,B)
        @test size(inputs.f) == (R,B)
        @test size(inputs.readmode) == (3, R, B)
    end

    @testset "Domain" begin
        # β ∈ [1, ∞)
        @test !any(x->x<1, inputs.βr)
        @test !any(x->x<1, inputs.βw)
        # g ∈ [0, 1]
        @test 0 <= inputs.ga[1] <= 1
        @test 0 <= inputs.gw[1] <= 1
        # e ∈ [0, 1]^W
        @test !any(x -> x<0 || x>1, inputs.e)
        # f ∈ [0, 1]^R
        @test !any(x -> x<0 || x>1, inputs.f)
        # π_i ∈ S_3
        @test round.(sum(inputs.readmode; dims=1)[1,:]; digits=5) == ones(R)
    end

    @testset "Differentiability" begin
        g = gradient(ξ, transforms) do ξ, transforms
            inputs = DNC.split_ξ(ξ, transforms)
            tot = 0
            for v in inputs
                tot += sum(v)
            end
            tot
        end
        @test !isnothing(g)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module Biplots

using Tables
using LinearAlgebra
using Statistics
using Printf

import Makie

# -----------------------
# TYPES OF NORMALIZATION
# -----------------------

center(X) = X .- mean(X, dims=1)

function logcenter(X)
  L  = log.(X)
  μn = mean(L, dims=1)
  μp = mean(L, dims=2)
  μ  = mean(L)
  L .- μn .- μp .+ μ
end

formtransform(X)  = (Z=center(X),    α=1.0, κ=√(size(X,1)-1))
covtransform(X)   = (Z=center(X),    α=0.0, κ=1.0           )
rformtransform(X) = (Z=logcenter(X), α=1.0, κ=1.0           )
rcovtransform(X)  = (Z=logcenter(X), α=0.0, κ=√size(X,2)    )

biplotof = Dict(
  :form => formtransform,
  :cov  => covtransform,
  :rform => rformtransform,
  :rcov => rcovtransform
)

"""
    biplot(table; kind=:form, dim=2, [aesthetics...])

Biplot of `table` of given `kind` in `dim`-dimensional space.

There are four kinds of biplots:

* `:form`  - standard biplot with shape parameter `α = 1`
* `:cov`   - standard biplot with shape parameter `α = 0`
* `:rform` - relative variation biplot with `α = 1`
* `:rcov`  - relative variation biplot with `α = 0`

# Biplot attributes

* `kind` - kind of biplot (`:form`, `:cov`, `:rform` or `:rcov`)
* `dim`  - number of dimensions `dim ∈ {2,3}` (default to `2`)

# Aesthetics attributes

* `axescolormap` - colormap of principal axes (default to theme colormap)
* `dotcolormap`  - colormap of sample dots (default to theme colormap)
* `textsize`     - text size in axes and dots (default to `16`)
* `axesbody`     - size of principal axes' body (depends on `dim`)
* `axeshead`     - size of principal axes' head (depends on `dim`)
* `axescolor`    - color of principal axes (default to `:black`)
* `axeslabel`    - names of principal axes (default to `x1,x2,...`)
* `dotsize`      - size of sample dots (depends on `dim`)
* `dotcolor`     - color of sample dots (default to `:black`)
* `dotlabel`     - names of sample dots (default to `1:nobs`)
* `showdots`     - show names of dots (default to `true`)
* `showlinks`    - show links between principal axes (default to`false`)

See https://en.wikipedia.org/wiki/Biplot.

## References

* Gabriel, K. 1971. [The Biplot Graphic Display of Matrices
  with Application to Principal Component Analysis]
  (https://www.jstor.org/stable/2334381)

* Aitchison, J. & Greenacre, M. 2002. [Biplots of Compositional data]
  (https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9876.00275)
"""
@Makie.recipe(Biplot, table) do scene
  Makie.Attributes(;
    # biplot attributes
    kind = :form,
    dim  = 2,

    # aesthetic attributes
    axescolormap = Makie.theme(scene, :colormap),
    dotcolormap  = Makie.theme(scene, :colormap),
    textsize     = 16,
    axesbody     = nothing,
    axeshead     = nothing,
    axescolor    = :black,
    dotsize      = nothing,
    dotcolor     = :black,
    dotlabel     = nothing,
    showdots     = true,
    showlinks    = false,
  )
end

function Makie.plot!(plot::Biplot{<:Tuple{Any}})
  # biplot attributes
  table = plot[:table][]
  kind  = plot[:kind][]
  dim   = plot[:dim][]

  # aesthetics attributes
  axescolormap = plot[:axescolormap][]
  dotcolormap  = plot[:dotcolormap][]
  textsize     = plot[:textsize][]
  axesbody     = plot[:axesbody][]
  axeshead     = plot[:axeshead][]
  axescolor    = plot[:axescolor][]
  dotsize      = plot[:dotsize][]
  dotcolor     = plot[:dotcolor][]
  dotlabel     = plot[:dotlabel][]
  showdots     = plot[:showdots][]
  showlinks    = plot[:showlinks][]

  # design matrix
  X = Tables.matrix(table)
  n, p = size(X)

  # defaults differ on 2 or 3 dimensions
  if isnothing(axesbody)
    axesbody = dim == 2 ? 2 : 0.01
  end
  if isnothing(axeshead)
    axeshead = dim == 2 ? 6 : 0.03
  end
  if axescolor isa AbstractVector{<:Number}
    min, max = extrema(axescolor)
    axesscale(x) = x / (max-min) - min / (max - min)
    axescolor = Makie.cgrad(axescolormap)[axesscale.(axescolor)]
  end
  if isnothing(dotsize)
    dotsize = dim == 2 ? 4 : 10
  end
  if isnothing(dotlabel)
    dotlabel = string.(1:n)
  end
  if dotcolor isa AbstractVector{<:Number}
    min, max = extrema(dotcolor)
    dotscale(x) = x / (max-min) - min / (max - min)
    dotcolor = Makie.cgrad(dotcolormap)[dotscale.(dotcolor)]
  end

  # sanity checks
  @assert dim ∈ [2,3] "dim must be 2 or 3"
  @assert kind ∈ [:form,:cov,:rform,:rcov] "$kind is not a valid kind of biplot"
  @assert length(dotlabel) == n "dotlabel must have length $n"

  # transformation
  Z, α, κ = biplotof[kind](X)

  # singular value decomposition
  U, σ, V = svd(Z)

  # variance explained
  σ² = σ[1:dim] .^ 2 / sum(σ .^ 2)

  # matrix factors X ≈ F*G'
  F = U[:,1:dim] .* (σ[1:dim] .^ α)'
  G = V[:,1:dim] .* (σ[1:dim] .^ (1 - α))'

  # plot principal axes
  points = fill(Makie.Point(ntuple(i->0., dim)), n)
  direcs = [Makie.Vec{dim}(v) for v in eachrow(G)] ./ κ
  Makie.arrows!(plot, points, direcs,
    linewidth  = axesbody,
    arrowsize  = axeshead,
    arrowcolor = axescolor,
    linecolor  = axescolor,
  )

  # plot axes names
  position = Tuple.(direcs)
  names = Tables.columnnames(table)
  Makie.text!(plot, collect(string.(names)),
    position = position,
    textsize = textsize,
    color = axescolor,
  )

  # plot links between axes
  if showlinks
    links = Makie.Vec{dim}[]
    for i in 1:p
      for j in i+1:p
        push!(links, direcs[i])
        push!(links, direcs[j])
      end
    end
    position = Tuple.(links)
    Makie.linesegments!(plot, position,
      color = :gray,
      linestyle = :dash,
      linewidth = 0.5,
    )
  end

  # plot samples
  points = [Makie.Point{dim}(v) for v in eachrow(F)]
  Makie.scatter!(plot, points,
    markersize = dotsize,
    color = dotcolor,
  )

  if showdots
    # plot samples names
    position = Tuple.(points)
    Makie.text!(plot, dotlabel,
      position = position,
      textsize = textsize,
      color = dotcolor,
    )
  end

  # plot variance explained
  minpos = fill(Inf, dim)
  for i in 1:dim
    for direc in direcs
      if direc[i] < minpos[i]
        minpos[i] = direc[i]
      end
    end
    for point in points
      if point[i] < minpos[i]
        minpos[i] = point[i]
      end
    end
  end
  textdim = [(@sprintf "Dim %d (%.01f " i 100*v)*"%)" for (i, v) in enumerate(σ²)]
  Makie.text!(plot, join(textdim, "\n"),
    position = Tuple(minpos),
  )
end

end