text:
# Compatible with JuMP v0.18

mutable struct Stat
  count1::Int64
  count2::Int64
  count3::Int64
  maxmintime::Float64
  inter_computation_time::Float64
end

function compute_robust_sol_KP(data, Γ, heur, problem)::Float64
  opt = sum(data.hc) + sum(data.d)
  print("Preprocessing: solving robust models ... ")
  if problem == "KP"
    d = copy(data.d)
    push!(d,0.0)
    for dev in d
      model = create_model()
      @variable(model, x[1:data.n], Bin)
      @constraint(model, knp, sum(data.w[i] * x[i] for i in 1:data.n) ≥ data.C)
      @objective(model, Min, sum( ( data.hc[i] + max(0, data.d[i]-dev) )*x[i] for i in 1:data.n ) )
      print("$dev ")
      optimize!(model)
      profit = Γ * dev + objective_value(model)
      if profit < opt
        opt = profit
      end
    end
  else
    d = copy(data.d)
    push!(d,0.0)
    for dev in d
      g = DiGraph(data.nV)
      for i in 1:data.n
        add_edge!(g, data.from[i], data.to[i])
      end
      dist = zeros(data.nV,data.nV)
      for i in 1:data.n
        dist[data.from[i],data.to[i]] = data.hc[i] + max(0, data.d[i] - dev)
      end
      all_dist = dijkstra_shortest_paths(g, data.s, dist).dists

      profit = Γ * dev + all_dist[data.t]
      if profit < opt
        opt = profit
      end
    end
  end

  println("")
  println("Robust solution is $opt")
  println("Alternating heuristic is $heur")

  if heur < opt
    return convert(Float64, heur)
  else
    return convert(Float64, opt)
  end
end

function generate_solutions_KP(data, Γ, solution, heur)
  # compute the optimal solution of the min-max problem
  function gen_sol(x::BitArray,i::Int64,w::Int64,c::Float64,opt::Float64)
    if time() - STARTTIME > TIME_LIMIT
      return 0
    end
    if (i == data.n)
      push!(X,copy(x))
    else
      i = i+1

      x[i] = true
      if c + data.hc[i] < heur # filter by optimality
        c = c + data.hc[i]
        gen_sol(x,i,w,c,heur)
        c = c - data.hc[i]
      end

      x[i] = false
      if w - data.w[i] >= data.C  # filter by feasibility
        w = w - data.w[i]
        gen_sol(x,i,w,c,opt)
      end
    end
  end

  X = Vector{BitArray{1}}() # set of all feasible vectors for the knapsack
  x = falses(data.n)
  i = 0
  w = sum(data.w)
  c = 0.0

  gen_sol(x,i,w,c,heur)

  println("|X|=$(size(X)[1])")
  solution.n_sol = size(X)[1]

  #Y = zeros(0,data.n) # set of maximal feasible vectors for the knapsack
  #for i in 1:size(X)[1]
  #   y = X[i,1:data.n]
  #  dominated = ( (y'data.hc)[1] < opt) # remove immediately solutions whith deterministic cost below the robust optimal solution
  #  for j in 1:size(X)[1]
  #    if dominated
  #      break
  #    end
  #    x = X[j,1:data.n]
  #    dominated = dominate(y,x)
  #  end
  #  if !dominated
  #    Y = [Y; y']
  #  end
  #end
  #println("|Y|=$(size(Y)[1])")
  #solution.n_sol_red = size(Y)[1]

return X
end

function generate_solutions_SP_BB(data, Γ, solution, heur)
  m = length(data.hc)
  n = data.nV
  model = create_model()
  @variable(model, 0 <= f[1:m] <= 1)
  @constraint(model, equilibrium[v in setdiff(1:n,[data.s,data.t])], sum(f[i] for i in data.δ⁺[v]) - sum(f[i] for i in data.δ⁻[v]) == 0 )
  @constraint(model, sum(f[i] for i in data.δ⁺[data.s]) == 1)
  @objective(model, Min, sum(data.hc[i]*f[i] for i in 1:m))

  #print(model)

  Tree = [ Vector{BitArray}() for i in 1:(m+1) ]
  push!(Tree[1],BitArray(0))
  nnodes = 1
  X_bool = Vector{BitArray{1}}()
  starttime = time()

  while true && (time() - starttime < TIME_LIMIT)
    availablenodes = [ i for i in 1:(m+1) if !isempty(Tree[i]) ]
    deepest = 0
    if isempty(availablenodes)
      break
    else
      deepest = maximum(availablenodes)
    end
    currentnode = pop!(Tree[deepest])

    for i in 1:m
      set_lower_bound(f[i], 0.0)
      set_upper_bound(f[i], 1.0)
    end
    for i in 1:length(currentnode)
      if currentnode[i]
        set_lower_bound(f[i], 1.0)
      else
        set_upper_bound(f[i], 0.0)
      end
    end
    optimize!(model)

    relax = objective_value(model)
    if relax < heur
      if deepest == m + 1
        push!(X_bool,copy(currentnode))
      else
        leftbranch = copy(currentnode)
        push!(leftbranch, false)
        push!(Tree[deepest+1],leftbranch)

        rightbranch = copy(currentnode)
        push!(rightbranch, true)
        push!(Tree[deepest+1],rightbranch)

        nnodes = nnodes + 2
      end
    end
  end

  println("** Enumerating paths require $nnodes nodes")
  return X_bool
end

function generate_solutions_SP(data, Γ, solution, heur)
  g = DiGraph(data.nV)
  for i in 1:data.n
    add_edge!(g, data.from[i], data.to[i])
  end
  dist = zeros(data.nV,data.nV)
  for i in 1:data.n
    dist[data.from[i],data.to[i]] = data.hc[i]
  end
  sp = Vector(undef,data.nV)
  for v in 1:data.nV
    all_dist = dijkstra_shortest_paths(g, v, dist).dists
    sp[v] = all_dist[data.t]
  end

  print("Computing X ")
  @time X_bool = get_all_paths(data, heur, sp) # set of all feasible vectors for theSP
  println("|X|=$(length(X_bool))")

  #X = BitArray(length(X_bool),data.n)
  #print("creating matrix X ")
  #@time X = [ X_bool[s][i] for s in 1:length(X_bool), i in 1:data.n ]

return X_bool
end

function algo_combi_k_resist(heur,data)

    function populate(types_of_deviations,k,typeod)
        if k == K
            if length(typeod) > 0
                push!(types_of_deviations,typeod)
            end
        else
            k = k+1
            populate(types_of_deviations,k,typeod)
            typeod2 = copy(typeod)
            push!(typeod2,k)
            populate(types_of_deviations,k,typeod2)
        end
    end

    function remove(to_remove)
        all_but_s = setdiff(1:length(costs),to_remove)
        X = X[all_but_s]
        X_set = X_set[all_but_s]
        costs = costs[all_but_s]
        costs_dev = costs_dev[all_but_s]
        costs_Γ = costs_Γ[all_but_s,:]
        max_dev_index = max_dev_index[all_but_s]
    end

    function next(k::Int,next_dev::Int,available_dev::Vector{Float64})::Int
        next_index = next_dev + 1
        if next_index > length(max_dev_index[X_ω[ω[k]][σ[k]]])
            return 0
        elseif available_dev[max_dev_index[X_ω[ω[k]][σ[k]]][next_index]] < 1 - 1/granularity + 0.00001
            return next_index
        else
            next(k,next_index,available_dev)
        end
    end

    function maxmin_discrete(deviations::Array{Int,1},γ::Int)
        ############## body of the function where the actual maxmin is computed
        i = length(deviations)
        cost = Vector(undef,K)
        if i == TOD
            for kk in 1:K
                cost[kk] = costs[X_ω[ω[kk]][σ[kk]]]
            end
            for tod in 1:TOD
                for dev in 1:deviations[tod]
                    for kk in types_of_deviations[tod]
                        #println(types_of_deviations[tod])
                        #println(shared_dev)
                        cost[kk] = cost[kk] + shared_dev[tod,dev]
                    end
                end
            end
            mincost = minimum(cost)
            if mincost > maxmincost
                maxmincost = mincost
            end
            ################ the recursive part to build the possible vectors of deviations ##########
        elseif length(deviations) == TOD - 1
            push!(deviations,min(γ,bounds[i + 1]))
            maxmin_discrete(deviations,0)
        else
            for gamma in 0:min(γ,bounds[i + 1])
                deviations2 = copy(deviations)
                push!(deviations2,gamma)
                maxmin_discrete(deviations2,γ-gamma)
            end
        end
    end

    function maxmin_convex(ω::Vector{Int},σ::Vector{Int})
        dev_indexes = Array{Int,1}()
        for i in N
            for k in 1:K
                if X[X_ω[ω[k]][σ[k]]][i]
                    push!(dev_indexes,i)
                    break
                end
            end
        end

        I = 1:length(dev_indexes)
        Ik = Array{Set{Int},1}(K)
        for k in 1:K
            Ik[k] = Set{Int}()
            for i in 1:length(dev_indexes)
                if X[X_ω[ω[k]][σ[k]]][dev_indexes[i]]
                    push!(Ik[k],i)
                end
            end
        end

        slave = create_model()
        #slave = Model()

        @variable(slave, Ω ≥ 0)
        @variable(slave, c[I] ≥ 0)
        @variable(slave, 0 ≤ δ[I] ≤ 1)

        @constraint(slave, sum(δ[i] for i in I) ≤ Γ)
        @constraint(slave, [k in 1:K], Ω ≤ costs[X_ω[ω[k]][σ[k]]] + sum(data.d[dev_indexes[i]]*δ[i] for i in Ik[k]))
        @objective(slave, Max, Ω)

        optimize!(slave)

        return objective_value(slave)
    end

    function next_continuous(k::Int,current_dev::Int,dev_left::Array{Float64,1})
        next_index = current_dev + 1
        if dev_left[max_dev_index[X_ω[ω[k]][σ[k]]][next_index]] > 0.0001
            return next_index
        else
            next_continuous(k,next_index,dev_left)
        end
    end

    function enumerate_resist(k::Int,statistics::Stat) # enumerate all k-solutions between first and last
        ω_first[k] = Γ*granularity
        if k > 1
            ω_first[k] = ω[k-1]
        end
        ω_last[k] = 1
        if k == K
            ω_last[k] = max( 1,Γ*granularity + 1 - sum(ω[1:K-1]) )
        end

        #println("Before the ω loop for k=$k: ω_first[k]=$(ω_first[k]) and ω_last[k]=$(ω_last[k])")
        for ω[k] in ω_first[k]:-1:ω_last[k]
            #println("ω[$k] = $(ω[k])")
            σ_first[k] = 1
            if k > 1
                if ω[k] == ω[k-1]
                    σ_first[k] = σ[k-1]+1
                end
            end
            σ_last[k] = length(X_ω[ω[k]])

            #println("Before the σ loop for k=$k: σ_first[k]=$(σ_first[k]) and σ_last[k]=$(σ_last[k])")
            for σ[k] in σ_first[k]:σ_last[k]
                #println("σ       $k    $(σ_first[k])     $(σ[k])   $(σ_last[k])    $(length(X_ω[ω[k]]))")
                #println("σ[$k]=$(σ[k])   $BREAK")

                k==1 && BREAK && break
                time() - STARTTIME > TIME_LIMIT && break

                if k == 1 && discrete # allocate the last deviations for the heuristic
                    for γ in 1:rem(Γ,K)
                        costs_dev[X_ω[ω[1]][σ[1]]] = costs_dev[X_ω[ω[1]][σ[1]]] + data.d[max_dev_index[X_ω[ω[1]][σ[1]]][γ+div(Γ,K)]]
                    end
                end

                if k < K
                    enumerate_resist(k+1,statistics)

                    ################ STOP CALLING ITSELF, BODY OF THE FUNCTION
                else
                    start = time()
                    statistics.count1 = statistics.count1 + 1
                    if minimum( [ costs_dev[X_ω[ω[kk]][σ[kk]]] for kk in 1:K ] ) < minmaxmincost
                        # cheap lower bound on the max-min
                        statistics.count2 = statistics.count2 + 1

                        # computing the true max-min is expensive because it requires the interesection between the k solutions
                        # Thus, we compute a greedy lower bound on the max-min
                        c = [costs[X_ω[ω[k]][σ[k]]] for k in 1:K]
                        #if discrete
                        next_dev = fill(1,K)
                        available_dev = zeros(data.n)
                        granularity = 1 # remove this to experiment stronger (but slower) LB in the convex case
                        #number_deviations = [ 0 for kk in 1:K ]
                        for γ in 1:Γ*granularity
                          if sum(next_dev) > 0
                            k_min = 0
                            c_min = 999999999999
                            for kk in 1:K
                              if next_dev[kk] > 0 && c[kk] < c_min
                                k_min = kk
                                c_min = c[kk]
                              end
                            end
                            index = max_dev_index[X_ω[ω[k_min]][σ[k_min]]][next_dev[k_min]]
                            available_dev[index] = available_dev[index] + 1/granularity
                            for kk in 1:K
                              if X[X_ω[ω[kk]][σ[kk]]][index]
                                c[kk] = c[kk] + 1/granularity * data.d[index]
                                if available_dev[index] + 1/granularity > 1 - 1/granularity + 0.00001
                                  #number_deviations[kk] = number_deviations[kk] + 1
                                  next_dev[kk] = next(kk,next_dev[kk],available_dev)
                                end
                              end
                            end
                          else
                            break
                          end
                        end

                        if minimum(c) <= minmaxmincost # compute the true max-min only if the LB is low enough
                            if discrete
                                statistics.count3 = statistics.count3 + 1
                                start_inter_computation_time = time()
                                for tod in 1:TOD
                                    shared_set = []
                                    for k in types_of_deviations[tod]
                                        if length(shared_set) == 0
                                            shared_set = X_set[X_ω[ω[k]][σ[k]]]
                                        else
                                            shared_set = intersect(shared_set,X_set[X_ω[ω[k]][σ[k]]])
                                        end
                                    end
                                    for k in setdiff(1:K,types_of_deviations[tod])
                                        shared_set = setdiff(shared_set,X_set[X_ω[ω[k]][σ[k]]])
                                    end
                                    shared_dev[tod,1:length(shared_set)] = sort([ data.d[i] for i in shared_set ], rev=true)
                                    if (length(types_of_deviations[tod]) == 1)
                                        bounds[tod] = min( length(shared_set) , ω[types_of_deviations[tod][1]] )
                                    else
                                        bounds[tod] = min( length(shared_set) , maximum( [ω[types_of_deviations[tod][kk]] for kk in 1:length(types_of_deviations[tod])] ) )
                                    end
                                end
                                statistics.inter_computation_time = statistics.inter_computation_time + time() - start_inter_computation_time
                                maxmincost = 0
                                startmaxmintime = time()
                                maxmin_discrete(Vector{Int}(),Γ)
                                statistics.maxmintime = statistics.maxmintime + time() - startmaxmintime
                            elseif convex
                                statistics.count3 = statistics.count3 + 1
                                startmaxmintime = time()
                                maxmincost = maxmin_convex(ω,σ)
                                statistics.maxmintime = statistics.maxmintime + time() - startmaxmintime
                            end
                            if maxmincost < minmaxmincost
                                minmaxmincost = maxmincost
                                BREAK = true # BREAK ALL FOR-LOOPS
                            end
                        end
                    end
                end
                if k == 1 && discrete
                    for γ in 1:rem(Γ,K) # desallocate the last deviations for the heuristic
                        costs_dev[X_ω[ω[1]][σ[1]]] = costs_dev[X_ω[ω[1]][σ[1]]] - data.d[max_dev_index[X_ω[ω[1]][σ[1]]][γ+div(Γ,K)]]
                    end
                end
                if k == 1 && !BREAK
                    push!(to_remove,X_ω[ω[1]][σ[1]])
                end
            end
        end
    end

    ############################### some initialization ############################
    maxmincost = 0
    minmaxmincost = sum(data.hc) + sum(data.d)
    heur = compute_robust_sol_KP(data, Γ, heur, problem)
    if minmaxmincost > heur
        minmaxmincost = heur
    end
    sol = solution()
    LP = false
    N = 1:data.n
    if problem == "KP"
        X = generate_solutions_KP(data, Γ, sol, heur)
    elseif problem == "SP"
        X = generate_solutions_SP(data, Γ, sol, heur)
    end
    START = time()
    r = size(X)[1]
    number_initial_sol = r
    S = 1:r
    X_time = time() - STARTTIME
    cheap_order = sortperm(data.d,rev=true)
    original_order = invperm(cheap_order)
    data.d = data.d[cheap_order]
    data.hc = data.hc[cheap_order]

    BREAK = false
    X_set, X_ω, costs, max_dev_index, costs_dev, costs_Γ = 0,0,0,0,0,0
    types_of_deviations, TOD, shared_dev, bounds, pre_computed_bounds, to_remove, prep_time, = 0,0,0,0,0,0,0
    statistics = Stat(0,0,0,0.0,0.0)

    try
      print("Changing order of items in X ")
      @time X = [ X[s][cheap_order] for s in S ]
      print("Define X_set ")
      @time X_set = [ [i for i in N if X[s][i] ] for s in S ]
      print("Define costs ")
      @time costs = [ sum(data.hc[X_set[s]]) for s in S ]
      cheap_order = sortperm(costs)
      print("Reorganize X ")
      @time X = X[cheap_order]
      print("Reorganize X_set ")
      @time X_set = X_set[cheap_order]
      print("Compute costs ")
      @time costs = [ costs[s] for s in cheap_order ]
      print("Compute max_dev_index ")
      @time max_dev_index = [ X_set[s][sortperm(data.d[X_set[s]],rev=true)] for s in S ]
      print("Compute costs_dev ")
      @time costs_dev = [ costs[s] + sum(data.d[max_dev_index[s][1:min(div(Γ,K),length(max_dev_index[s]))]]) for s in S ]
      @time if convex
          for s in S
              if length(max_dev_index[s]) >= div(Γ,K)+1
                  costs_dev[s] = costs_dev[s] + rem(Γ,K)/K * data.d[max_dev_index[s][div(Γ,K)+1]]
              end
          end
      end
      costs_Γ = Matrix{Float64}(undef,r,granularity*Γ)
      print("Compute costs_Γ ")
      @time for s in S, g in 1:granularity*Γ
          γ = g/granularity
          floor_γ = convert(Int,floor(γ))
          frac_γ = γ - floor_γ
          costs_Γ[s,g] = costs[s] + sum(data.d[max_dev_index[s][1:min(floor_γ,length(max_dev_index[s]))]])
          if length(max_dev_index[s]) >= 1+floor_γ
            costs_Γ[s,g] = costs_Γ[s,g] + frac_γ * data.d[max_dev_index[s][1+floor_γ]]
          end
      end

      ############################### start the enumeration algorithm #####################################

      types_of_deviations = Vector{OrderedSet{Int}}()
      populate(types_of_deviations,0,OrderedSet{Int}())
      TOD = length(types_of_deviations)
      shared_dev = Matrix(undef,TOD,data.n)
      bounds = Vector{Int}(undef,TOD)
      pre_computed_bounds = [ length(X_set[s]) for s in 1:r ]
      to_remove = []

      prep_time = time() - START
      ############### Repeat until loop ################
      println("Starting enumeration part ...")
      X_ω = Vector{Vector{Int}}(undef,granularity*Γ)
      while true
          r = size(X)[1]
          S = 1:r
          for g in 1:granularity*Γ
              X_ω[g] = Vector{Int}()
          end
          γ = Vector{Int}(undef,r)
          for s in S
              if costs[s] > minmaxmincost
                  γ[s] = 0
              else
                  #println("$minmaxmincost      $(costs_Γ[s,:])")
                  γ[s] = minimum( [ g for g in 1:granularity*Γ if costs_Γ[s,g] >= minmaxmincost ] ) # Can be sped-up if need be?
                  push!(X_ω[γ[s]],s)
              end
          end
          println("Current best cost is $minmaxmincost")
          #for g in 1:granularity*Γ
          #    println("$g  $(length(X_ω[g]))")
          #end

          BREAK = false
          @time enumerate_resist(1,statistics)

          for s in S
              if γ[s] == 0
                  push!(to_remove,s)
              end
          end
          remove(to_remove)
          to_remove = []

          !BREAK && break
          println(time() - STARTTIME)
      end
    catch error
      println(error)
       if isa(error, OutOfMemoryError)
        statistics.maxmintime = -1
        statistics.inter_computation_time = -1
      end
    end
    ############# End of the repeat loop ###################

    println("optimal solution = $minmaxmincost")

    return (time() - STARTTIME, minmaxmincost, number_initial_sol, X_time, prep_time, statistics.count1, statistics.count2, statistics.count3, statistics.maxmintime, statistics.inter_computation_time)
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
## Fun.jl docs

# Constructors

"""
    Fun(s::Space,coefficients::AbstractVector)

returns a `Fun` with the specified `coefficients` in the space `s`
"""
Fun(::Space,::AbstractVector)

"""
    Fun(f,s::Space)

return a `Fun` representing the function, number, or vector `f` in the
space `s`.  If `f` is vector-valued, it returns a vector-valued analogue
of `s`.
"""
Fun(_,::Space)

"""
    Fun(f,d::Domain)

returns `Fun(f,Space(d))`, that is, it uses the default space for the specified
domain.
"""
Fun(_,::Domain)

"""
    Fun(s::Space)

returns `Fun(identity,s)`
"""
Fun(::Space)

"""
    Fun(f)

returns `Fun(f,Chebyshev())`
"""
Fun(f)

"""
    Fun()

returns `Fun(identity,Chebyshev())`.
"""
Fun()

"""
    ones(d::Space)

Return the `Fun` that represents the function one on the specified space.
"""
ones(::Space)

"""
    zeros(d::Space)

Return the `Fun` that represents the function one on the specified space.
"""
zeros(::Space)

# accessors

"""
    domain(f::Fun)

returns the domain that `f` is defined on
"""
domain(fun::Fun)

"""
    setdomain(f::Fun,d::Domain)

returns `f` projected onto `domain`
"""
setdomain(::Fun,::Domain)

"""
    space(f::Fun)

returns the space of `f`
"""
space(f::Fun)

"""
    values(f::Fun)

returns `f` evaluated at `points(f)`
"""
values(::Fun)

"""
    points(f::Fun)

returns a grid of points that `f` can be transformed into values
and back
"""
points(::Fun)

"""
    points(s::Space,n::Integer)

returns a grid of approximately `n` points, for which a transform exists
from values at the grid to coefficients in the space `s`.
"""
points(::Space,::Integer)

"""
    extrapolate(f::Fun,x)

returns an extrapolation of `f` from its domain to `x`.
"""
extrapolate(::Fun,x)

"""
    coefficients(f::Fun) -> Vector

returns the coefficients of `f`, corresponding to the space `space(f)`.
"""
coefficients(::Fun)

"""
    coefficients(f::Fun,s::Space) -> Vector

returns the coefficients of `f` in the space `s`, which
may not be the same as `space(f)`.
"""
coefficients(::Fun,::Space)

"""
    coefficients(cfs::AbstractVector,fromspace::Space,tospace::Space) -> Vector

converts coefficients in `fromspace` to coefficients in `tospace`
"""
coefficients(::AbstractVector,::Space,::Space)

"""
    ncoefficients(f::Fun) -> Integer

returns the number of coefficients of a fun
"""
ncoefficients(::Fun)

"""
    stride(f::Fun)

returns the stride of the coefficients, checked
numerically
"""
stride(::Fun)

## Modifiers

"""
   chop(f::Fun,tol) -> Fun

reduces the number of coefficients by dropping the tail that is below the specified tolerance.
"""
chop(::Fun,_)

"""
    reverseorientation(f::Fun)

return `f` on a reversed orientated contour.
"""
reverseorientation(::Fun)

## Spaces

"""
    canonicalspace(s::Space)

returns a space that is used as a default to implement missing functionality,
e.g., evaluation.  Implement a `Conversion` operator or override `coefficients` to support this.
"""
ApproxFun.canonicalspace(::Space)

"""
    transform(s::Space,vals::Vector)

Transform values on the grid specified by `points(s,length(vals))` to coefficients in the space `s`.
Defaults to `coefficients(transform(canonicalspace(space),values),canonicalspace(space),space)`
"""
transform(::Space,::Vector)

"""
    itransform(s::Space,coefficients::AbstractVector)

Transform coefficients back to values.  Defaults to using `canonicalspace` as in `transform`.
"""
itransform(::Space, ::AbstractVector)

"""
    evaluate(coefficients::AbstractVector, sp::Space, x)

Evaluates the expansion at a point `x`.
If `x` is in the domain, then this should return zero.
"""
evaluate(::AbstractVector, ::Space, _)

"""
    spacescompatible

Specifies equality of spaces while also supporting `AnyDomain`.
"""
spacescompatible(::Space,::Space)

"""
    conversion_type(a::Space,b::Space)

returns a `Space` that has a banded conversion operator to both `a` and `b`.
Override `ApproxFun.conversion_rule` when adding new `Conversion` operators.
"""
conversion_type(::Space,::Space)

"""
    dimension(s::Space)

returns the dimension of `s`, which is the maximum number of coefficients.
"""
dimension(::Space)

## Operator.jl docs

"""
    Operator{T}

is an abstract type to represent linear operators between spaces.
"""
Operator

"""
    domainspace(op::Operator)

gives the domain space of `op`.  That is, `op*f` will first convert `f` to
a `Fun` in the space `domainspace(op)` before applying the operator.
"""
domainspace(::Operator)

"""
    rangespace(op::Operator)

gives the range space of `op`.  That is, `op*f` will return a `Fun` in the
space `rangespace(op)`, provided `f` can be converted to a `Fun` in
`domainspace(op)`.
"""
rangespace(::Operator)

"""
    bandwidths(op::Operator)

returns the bandwidth of `op` in the form `(l,u)`, where `l ≥ 0` represents
the number of subdiagonals and `u ≥ 0` represents the number of superdiagonals.
"""
bandwidths(::Operator)

"""
    promotedomainspace(S::Operator,sp::Space)

returns the operator `S` but acting on the space `sp`.
"""
promotedomainspace(::Operator,::Space)

"""
    promoterangespace(S::Operator,sp::Space)

returns the operator `S` acting on the same space, but now return
functions in the specified range space `sp`
"""
promoterangespace(::Operator,::Space)

"""
    choosedomainspace(S::Operator,rangespace::Space)

returns a space `ret` so that `promotedomainspace(S,ret)` has the
specified range space.
"""
choosedomainspace(::Operator,::Space)

"""
    op[k,j]

returns the `k`th coefficient of `op*Fun([zeros(j-1);1],domainspace(op))`.
"""
getindex(::Operator,k,j)

"""
    op[f::Fun]

constructs the operator `op*Multiplication(f)`, that is, it multiplies on the right
by `f` first.  Note that `op*f` is different: it applies `op` to `f`.
"""
getindex(::Operator,::Fun)

"""
    Conversion(fromspace::Space,tospace::Space)

represents a conversion operator between `fromspace` and `tospace`, when available.
"""
Conversion(::Space,::Space)

ext: jl


text:
"""
```julia
m = arnoldi!(Hess, q, m, Ψ, H, dt=1.0; extended=true, norm_min=1e-15)
```

Calculate the Hessenberg matrix and Arnoldi vectors of `H dt`, from `Ψ`.

For a given order `m`, the `m×m` Hessemberg matrix is calculated and stored in
in the pre-allocated `Hess`. Further  an array of `m` normalized Arnoldi
vectors is stored in in the pre-allocated `q`, plus one additional unnormalized
Arnoldi vector.  The unnormalized `m+1`st vector could be used to easily
extend a given `m×m` Hessenberg matrix to a `(m+1)×(m+1)` matrix.

If the extended Hessenberg matrix is requested (`extended=true`, default), the
`m+1`st Arnoldi vector is also normalized, and it's norm will be stored in
`m+1, m` entry of the (extended) Hessenberg matrix, which is an `(m+1)×(m+1)`
matrix.

Return the size `m` of the calculated Hessenberg matrix. This will usually be
the input `m`, except when the Krylov dimension of `H` starting from `Ψ` is
less then `m`. E.g., if `Ψ` is an eigenstate of `H`, the returned `m` will
be 1.

See <http://en.wikipedia.org/wiki/Arnoldi_iteration> for a description of
the algorithm.

# Arguments

- `Hess::Matrix{ComplexF64}`: Pre-allocated storage for the Hessemberg matrix.
   Can be uninitialized on input. The matrix must be at least of size `m×m`, or
   `(m+1)×(m+1)` if `extended=true`. On output, the `m×m` sub-matrix of `Hess`
   (with the returned output `m`) will contain the Hessenberg matrix, and all
   other elements of `Hess` be be set to zero.
- `q`: Pre-allocated array of states similar to `Ψ`, as storage for the
  calculated Arnoldi vectors. These may be un-initialized on input. Must be at
  least of length `m+1`
- `m`: The requested dimensions of the output Hessenberg matrix.
- `Ψ`: The starting vector for the Arnoldi procedure. This can be of any type,
   as long as `Φ = H * Ψ` results in a vector similar to `Ψ`, there is an inner
   products of `Φ` and `Ψ` (`Ψ⋅Φ` is defined), and `norm(Ψ)` is defined.
- `H`: The operator (up to `dt`) for which to calculate the Arnoldi procedure.
  Can be of any type, as long as `H * Ψ` is defined.
- `dt`: The implicit time step; the total operator for which to calculate the
  Arnoldi procedure is `H * dt`
- `extended`: If `true` (default), calculate the extended Hessenberg matrix,
  and normalized the final Arnoldi vector
- `norm_min`: the minimum value of the norm of `Ψ` at which `Ψ` should be
   considered the zero vector
"""
function arnoldi!(Hess::Matrix{ComplexF64}, q::Array{T}, m::Int64, Ψ::T, H,
                  dt::Float64=1.0; extended=true, norm_min=1e-15) where T
    if extended
        dim_hess = m + 1
    else
        dim_hess = m
    end
    @assert size(Hess, 1) >= dim_hess && size(Hess, 2) >= dim_hess
    @assert length(q) >= m + 1
    fill!(Hess, 0)
    copyto!(q[1], Ψ)
    for j = 1:m
        mul!(q[j+1], H, q[j])
        for i = 1:j  # Orthogonalization with Gram-Schmidt
            Hess[i, j] = dt * (q[i] ⋅ q[j+1]) # = dt ⟨qᵢ|qⱼ₊₁⟩
            axpy!(-Hess[i,j] / dt, q[i], q[j+1])  # qⱼ₊₁ += -(Hessᵢⱼ/dt) qᵢ
        end
        if (j < m) || extended
            h = norm(q[j+1])
            Hess[j+1,j] = dt * h
            if h < norm_min
                # dimensionality exhausted. Returning reduced m
                m = j
                break
            end
            lmul!(1/h, q[j+1])
        end
    end
    return m
end

""" Extend dimension of Hessenberg matrix by one.

```julia
extend_arnoldi!(Hess, q, m, H, dt; norm_min=1e-15)
```

extends the entries in `Hess` from size (m-1)×(m-1) to size m×m, and the list
`q` of Arnoldi vectors from m to (m+1). It is assumed that the input `Hess` was
created by a call to [`arnoldi!`](@ref) with `extended=false` or a previous
call to `extend_arnoldi!`. Note that `Hess` itself is not resized, so it must
be allocated to size m×m or greater on input.
"""
function extend_arnoldi!(Hess, q, m, H, dt::Float64=1.0; norm_min=1e-15)
    h = norm(q[m])
    (h < norm_min) && return m
    Hess[m, m-1] = dt * h
    lmul!(1/h, q[m])
    mul!(q[m+1], H, q[m])
    for i in 1:m
        Hess[i, m] = dt * (q[i] ⋅ q[m+1])
        axpy!(-Hess[i,m] / dt, q[i], q[m+1])
    end
    # everything below the first sub-diagonal should be zero. We'll check the
    # last row (previous rows were checked in earlier extend_arnoldi!)
    @assert all(Hess[m,1:m-2] .== 0.0)
    return Hess
end

"""
```julia
diagonalize_hessenberg_matrix(Hess, m; accumulate=false)
```

Diagonalize the m × m top left submatrix of the given Hessenberg matrix.

If `accumulate` is `true`, return the concatenated eigenvalues for
`Hess[1:1,1:1]` to `Hess[1:m,1:m]`, that is, all sumatrices of size 1 through
`m`.
"""
function diagonalize_hessenberg_matrix(Hess, m; accumulate=false)
    j_min = m
    j_max = m
    if accumulate
        j_min = 1
        eigenvals = zeros(ComplexF64, (m * (m+1)) ÷ 2)
    else
        eigenvals = zeros(ComplexF64, m)
    end
    offset = 0
    for j = j_min:j_max
        if j == 1
            eigenvals[1] = Hess[1, 1]
        elseif j == 2
            a = Hess[1,1]
            c = Hess[2,1]
            b = Hess[1,2]
            d = Hess[2,2]
            s = sqrt(a^2 + 4 * b * c - 2*a*d + d^2)
            eigenvals[offset+1] = 0.5 * (a+d-s)
            eigenvals[offset+2] = 0.5 * (a+d+s)
        else
            # TODO: Wrap in UpperHessenberg?
            # https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Special-matrices
            # TODO: use a view?
            eigenvals[offset+1:offset+j] .= eigvals(Hess[1:j,1:j])
        end
        offset += j
    end
    return eigenvals
end


text:
typeof("abc")              # 返回 String， String是 AbstractString的子类型
isa("abc", AbstractString) # true
isa(1, Float64)            # false, 1 是整数而不是浮点数
isa(1.0, Float64)          # true
1.0 isa Number             # 另一种语法；true，Number （数）是一种抽象类型
supertype(Int64)           # Int64 的父类型
subtypes(Real)             # 抽象类型 Real 的子类型

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
module WaterMeanMeanForce

using LinearAlgebra: diagm, dot, eigen, norm, normalize, pinv
using Statistics: mean, std
using StatsBase: fit, Histogram

using PrettyTables: pretty_table
using ProgressMeter: Progress, ProgressWrapper

export
    list_models,
    get_model,

    list_integrators,
    get_integrator,

    MolecularDynamicsParameters,
    run_constrained_equilibration,
    run_constrained_production,
    run_restrained_equilibration,
    run_restrained_production,

    unbias_histograms,

    harmonic_exact_slope

include("misc.jl")
include("units.jl")
include("stats.jl")

include("state.jl")
include("model.jl")

include("classical.jl")
include("constraint.jl")

include("force.jl")

include("integrator.jl")
include("estimator.jl")
include("run.jl")

include("wham.jl")

end


text:
# # Wind over deep mixed layer
#
# In this example, we blow an inertially-oscillating atmospheric wind
# over a deep, pre-mixed boundary layer. We compare the evolution of the
# boundary layer with the K-Profile parameterization and the TKEMassFlux model.

using Plots, Printf

using OceanTurb

using OceanTurb.TKEMassFlux: ModelBoundaryConditions # This function works for TKE and KPP

# ## Resolution and domain size
#
# We use a grid with 256 points and 1 meter grid spacing.

grid = UniformGrid(N=256, H=256)

# ## Coriolis
constants = Constants(f=1e-4) # s⁻¹

# ## The inertial wind
#
# We blow an inertially-rotating wind with wind stress 1e-4 m² s⁻².
# This corresponds roughly to 6 m s⁻¹ atmospheric wind speed.

τ = 1e-4 # m² s⁻²

## Inertially rotating means Qᵘ + i Qᵛ = e⁻ⁱᶠᵗ τ:
Qᵘ_inertial(t, f) =   τ * cos(f * t)
Qᵛ_inertial(t, f) = - τ * sin(f * t)

## Boundary condition functions
Qᵘ(model) = Qᵘ_inertial(model.clock.time, model.constants.f)
Qᵛ(model) = Qᵛ_inertial(model.clock.time, model.constants.f)

U_boundary_conditions = BoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
V_boundary_conditions = BoundaryConditions(top = FluxBoundaryCondition(Qᵛ))

boundary_conditions = ModelBoundaryConditions(U = U_boundary_conditions,
                                              V = V_boundary_conditions)

# We are now ready to build both models.

tke_model = TKEMassFlux.Model(     grid = grid,
                                stepper = :BackwardEuler,
                                    bcs = boundary_conditions,
                              constants = constants)

kpp_model = ModularKPP.Model(     grid = grid,
                               stepper = :BackwardEuler,
                                   bcs = boundary_conditions,
                             constants = constants)

# ## Stratification
#
# Finally, we prescribe a linear stratification with buoyancy frequency

N² = 5e-6 # s⁻²

# The temperature gradient given the default thermal expansion coefficient `α`
# and gravitational acceleration `g` is

dTdz = N² / (constants.α * constants.g)

# ## Initial mixed layer depth
#
# We prescribe an initial mixed layer depth

h = 100 # m

# such that the temperature gradient is a step function,

T₀(z) = 20 + dTdz * min(z + h, 0)

# We then set initial conditions and the temperature boundary condition for both models,

for model in (tke_model, kpp_model)
    model.solution.T = T₀
    model.bcs.T.bottom = GradientBoundaryCondition(dTdz)
end

# Finally, a few shortcuts make plotting a little easier.
T_tke = tke_model.solution.T
U_tke = tke_model.solution.U
V_tke = tke_model.solution.V

T_kpp = kpp_model.solution.T
U_kpp = kpp_model.solution.U
V_kpp = kpp_model.solution.V

z = nodes(U_tke)

speed(U, V) = sqrt.(data(U).^2 .+ data(V).^2)

max_speed = 0.6

# ## Run the model
#
# We use a 1 minute time-step, which is quite small and safe for these models.

Δt = 1minute

# Animate!

anim = @animate for iplot = 1:201

    time_label = Plots.text(@sprintf("t = %.2f hours", tke_model.clock.time / hour))

    temperature_plot = plot(data(T_tke), z,
                            legend = :topleft,
                             label = "TKE",
                            xlabel = "Temperature (ᵒC)",
                            ylabel = "z (m)",
                             xlims = (19.7, 20.05),
                           )

    speed_plot = plot(speed(U_tke, V_tke), z,
                           legend = false,
                           xlabel = "Speed (m/s)",
                            xlims = (0, max_speed),
                      annotations = (0.4, -200, time_label),
                     )

    velocities_plot = plot(data(U_tke), z,
                                legend = :bottomright,
                                 label = "U (TKE)",
                                 color = :blue,
                                xlabel = "U, V (m/s)",
                                 xlims = (-max_speed, max_speed),
                          )

    plot!(velocities_plot, data(V_tke), z, label="V (TKE)", linecolor=:red)
    plot!(velocities_plot, data(U_kpp), z, label="U (KPP)", linecolor=:blue, linestyle=:dash)
    plot!(velocities_plot, data(V_kpp), z, label="V (KPP)", linecolor=:red, linestyle=:dash)

    plot!(temperature_plot, data(T_kpp), z, label = "KPP")
    plot!(speed_plot, speed(U_kpp, V_kpp), z)

    plot(temperature_plot, speed_plot, velocities_plot, layout = (1, 3), size=(800, 500))

    run_until!(tke_model, Δt, iplot * hour / 2)
    run_until!(kpp_model, Δt, iplot * hour / 2)
end

mp4(anim, "wind_over_deep_mixed_layer.mp4", fps = 15)


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary, RecursiveArrayTools
using Test

DDEProblemLibrary.importddeproblems()

const prob_ip = DDEProblemLibrary.prob_dde_constant_1delay_ip
const prob_scalar = DDEProblemLibrary.prob_dde_constant_1delay_scalar

const alg = MethodOfSteps(BS3(); constrained = false)

@testset for inplace in (true, false)
  prob = inplace ? prob_ip : prob_scalar

  @testset "integrator" begin
    integrator = init(prob, alg, dt= 0.01)
    solve!(integrator)

    u = recursivecopy(integrator.sol.u)
    t = copy(integrator.sol.t)

    reinit!(integrator)
    integrator.dt = 0.01
    solve!(integrator)

    @test u == integrator.sol.u
    @test t == integrator.sol.t
  end

  @testset "solution" begin
    integrator = init(prob, alg, dt= 0.01, tstops = [0.5], saveat = [0.33])
    sol = solve!(integrator)

    u = recursivecopy(sol.u)
    t = copy(sol.t)

    reinit!(integrator)
    integrator.dt = 0.01
    sol = solve!(integrator)

    @test u == sol.u
    @test t == sol.t
  end
end


text:
#####################################################################
# Example: integrate out the hyperparameters of a Gaussian process. #
#####################################################################

# Extra packages needed for this example
# If missing, they must be installed explicitly!

using AbstractGPs, PyPlot, Printf, LinearAlgebra

# Define log-likelihood

function logp(θ; x=x, y=y)

    local kernel = exp(θ[2]) * (Matern52Kernel() ∘ ScaleTransform(exp(θ[1])))

    local f = GP(kernel)

    local fx = f(x, exp(θ[3]))

    logpdf(fx, y)

end

# Generate some synthetic data

N = 25

σ = 0.1

x = rand(N)*10

y = sin.(x) .+ randn(N)*σ

xtest = collect(-1.0:0.1:11.0)

ytest = sin.(xtest)

# Approximate posterior with Gaussian

postθ, = VI( θ ->  logp(θ; x=x, y=y), [randn(3)*2 for i=1:10], S = 200, iterations = 100, show_every = 10)

# Draw samples from posterior and plot

for i in 1:3

    # draw hyperparameter sample from approximating Gaussian distribution
    local θ = rand(postθ)

    # instantiate kernel
    local sample_kernel = exp(θ[2]) * (Matern52Kernel() ∘ ScaleTransform(exp(θ[1])))

    # intantiate kernel, GP object and calculate posterior mean and covariance for the training data x, y generated above
    local f = GP(sample_kernel)
    local p_fx = AbstractGPs.posterior(f(x, exp(θ[3])), y)
    local μ, Σ = AbstractGPs.mean_and_cov(p_fx, xtest)

    figure()
    plot(x, y, "ko",label="Training data")
    plot(xtest, ytest, "b-", label="True curve")

    plot(xtest, μ, "r-")
    fill_between(xtest, μ.-sqrt.(diag(Σ)),μ.+sqrt.(diag(Σ)), color="r", alpha=0.3)

    title(@sprintf("GP posterior, sampled hyperparameters %.2f, %.2f, %.2f", exp(θ[1]),exp(θ[2]),exp(θ[3])))
    legend()

end


text:
include("graph.jl")

function contact_freq_to_distance_matrix(contact_file::AbstractString; remove_inf=true, scale=1, exponent=1/3)
    # turns a csv file of contact frequences into a distance matrix
    contact_data = readcsv(contact_file, Float64)
    # TODO: use the graphical approach... distances are 1/freq, use shortest
    # paths algorithm to get distances that don't exist

    # arbitrarily scale distances by some scale
    distances = scale./(contact_data).^exponent
    dist = distances
    # removing zero points
    if remove_inf
        dist_index = 1
        while dist_index <= size(dist)[1]
            if mapreduce(x->(isinf(x) || x==0), &, dist[dist_index, 1:end-1])
                dist = [dist[1:dist_index-1, :], dist[dist_index+1:end, :]]
                dist = [dist[:, 1:dist_index-1] dist[:,dist_index+1:end]]
            else
                dist_index += 1
            end
        end
    else

    end
    distances = dist
    n, n2 = size(distances)
    g = init_graph()
    for i in range(1, n)
        add_node!(g, i)
    end
    for i in range(1, n-1)
        for i2 in (i+1):n
            if distances[i, i2] != Inf
                add_edge!(g, i, i2, distances[i, i2])
            else
                #distances[i, i2] = 1.0
                #add_edge!(g, i, i2, distances[i, i2])
            end
        end
    end
    # done with constructing graph, now to find the distances
    dist = all_shortest_paths(g)
    # TODO: what if there's an "Inf" in the distance? replace by 1?
    # this is an incredibly hacky way of doing things, and it might mess
    # everything up entirely.
    # ... and it does seem to mess up everything entirely.
    # how else would this be done?
    #dist[dist.==Inf] = 1
    return dist
end

function mds_metric(dist::Array{Float64, 2}; dims = 3)
    # 'dist' is a matrix of distances (should be symmetric)
    # metric MDS - assumes that the data is metric
    # TODO: do some crazy matrix operations
    N, N = size(dist)
    d = zeros(N)
    b = 1/N^2*sum([sum([dist[j,k]^2 for k in j:N]) for j in range(1,N)])
    for i in range(1, N)
        d[i] = 1/N * sum([dist[i, j]^2 for j in range(1, N)]) - b
    end
    M = zeros(N, N)
    for i in range(1, N)
        for j in range(1, N)
            M[i,j] = 1/2*(d[i] + d[j] - dist[i,j]^2)
        end
    end
    lambda, E = eig(M)
    # largest 3 eigenvalues
    v = E[:, end-dims+1:end] .* sqrt(lambda[end-dims+1:end])'
    return fliplr(v[:,1:dims]) #[v[:,3] v[:,2] v[:,1]]
end

function test()
    m = "../test_data/combined_You_Pro_merged_chr18.csv"
    dist = contact_freq_to_distance_matrix(m)
    return dist, mds_metric(dist)
end


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
cd(@__DIR__)
using Revise
using ZigZagBoomerang
const ZZB = ZigZagBoomerang

using ZigZagBoomerang: Trace, sevent, waiting_time_ref, freezing_time!, ab,
     smove_forward!, grad_correct!, λ, sλ̄, reflect_sticky!, freezing_time,
     refresh_sticky_vel!, sep
using DataStructures
using LinearAlgebra
#using Observables
using Colors
using Random
using SparseArrays
using Test
using FileIO
using Dates
using Statistics
@time using Makie, AbstractPlotting
println(now())
using ForwardDiff
const ρ0 = 0.0
const d = 2
const dist = 1.
const R = 3.0
 r = -R:0.05:R
#=
mutable struct Level
    story
    gradϕ!
    logdensity
    potential
    ys
    minpoints
    T
    Target
    User
    θ0
    x0
end
=#
mutable struct Game
    kill_on_boundary
    pressed
    released
    ys
    canvas
    score
    speed
    timeleft
    auto
    potential
end
ZigZagBoomerang.λ(∇ϕx, θ, F::Union{BouncyParticle, Boomerang}) = ZZB.pos(dot(∇ϕx, θ)) + 0.0

function bouncy_inner!(game, X, Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, b, t′, f, θf, tfrez, tref, told, (acc, num),
    Flow::Union{ZigZag, BouncyParticle, Boomerang}, κ, args...; strong_upperbounds = false, factor=1.5, adapt=false)
    canvas = game.canvas
    factorized = Flow isa ZigZag
    while true
        ispressed(canvas.scene, Keyboard.q) && error("end")
        for (key, pressed) in game.pressed
            if pressed && !ispressed(canvas.scene, key)
                game.pressed[key] = false
            end
        end
        for (key, released) in game.released
            if released && ispressed(canvas.scene, key)
                game.released[key] = false
            end
        end

        finished = remove(game, x)
        finished |= ispressed(canvas.scene, Keyboard.enter)
        if ispressed(canvas.scene, Keyboard.a) && !game.pressed[Keyboard.a]
            game.auto[] = !game.auto[]
            game.pressed[Keyboard.a] = true
        end
        X[] = [point(POTENTIAL[], x)]
        sleep(0.002)
        yield()
        tᶠ, i = findmin(tfrez) # could be implemented with a queue
        tt, j = findmin([tref, tᶠ, t′...])
        τ = tt - t
        t, x, θ = smove_forward!(τ, t, x, θ, f, Flow)
        # move forward
        if j == 1 # refreshments of velocities
            if (!game.auto[] && ispressed(canvas.scene, Keyboard.c)) || (game.auto[] && !ispressed(canvas.scene, Keyboard.c))
                θ, θf = refresh_sticky_vel!(θ, θf, f, Flow)
            elseif ispressed(canvas.scene, Keyboard.equal)
                θ *= 1.1
                θf *= 1.1
            elseif ispressed(canvas.scene, Keyboard.minus)
                θ *= 0.9
                θf *= 0.9
            end
            tref = t + waiting_time_ref(Flow) # regenerate refreshment time
            b = ab(x, θ, c, Flow) # regenerate reflection time
            told = t
            t′ = t + poisson_times(b)
            tfrez = freezing_time!(tfrez, t, x, θ, f, Flow)
            for i in eachindex(f) # make function later...
                if !f[i]
                    tfrez[i] = t - log(rand())/(κ[i]*abs(θf[i]))
                end
            end
            game.speed[] = norm(θ)
        elseif j == 2 # get frozen or unfrozen in i
            if f[i] # if free
                if abs(x[i]) > 1e-8
                    tfrez[i] = t + freezing_time(x[i], θ[i], Flow.μ[i], Flow) # wrong zero of curve
                    error("x[i] = $(x[i]) !≈ 0 at $(tfrez[i])")
                end
                x[i] = -0*θ[i]
                θf[i], θ[i] = θ[i], 0.0 # stop and save speed
                f[i] = false # change tag
                tfrez[i] = t - log(rand())/(κ[i]*abs(θf[i])) # sticky time
                # tfrez[i] = t - log(rand()) # option 2
                if !(strong_upperbounds) #not strong upperbounds, draw new waiting time
                    b = ab(x, θ, c, Flow) # regenerate reflection time
                    told = t
                    t′ = t .+ poisson_times(b)
                end
                game.released[Keyboard.s] = false
            elseif !ispressed(canvas.scene, Keyboard.s)  && !game.released[Keyboard.s] # is frozen ->  unfreeze
                @assert x[i] == 0 && θ[i] == 0
                θ[i], θf[i] = θf[i], 0.0 # restore speed
                f[i] = true # change tag
                tfrez[i] = t + freezing_time(x[i], θ[i], Flow.μ[i], Flow)
                b = ab(x, θ, c, Flow) # regenerate reflection time
                told = t
                t′ = t .+ poisson_times(b)
                game.released[Keyboard.s] = true
            else
                tfrez[i] = t - log(rand())/(κ[i]*abs(θf[i])) # sticky time
            end
        elseif factorized
            i = j - 2

            ∇ϕx = ∇ϕ!(∇ϕx, x, args...)
            l, lb = ZZB.sλ(∇ϕx[i], i, x, θ, Flow), ZZB.sλ̄(b[i], t - told)
            num += 1

            left = (θ[1]*θ[2] > 0) ⊻ (i == 2)

            if ((game.auto[] || !inbounds(game, x)) && rand()*lb <= l ) ||  ( left && ispressed(canvas.scene, Keyboard.left) && !game.pressed[Keyboard.left]) ||
                ( !left && ispressed(canvas.scene, Keyboard.right) && !game.pressed[Keyboard.right])
                if ispressed(canvas.scene, Keyboard.right)
                    game.pressed[Keyboard.right] = true
                end
                if ispressed(canvas.scene, Keyboard.left)
                    game.pressed[Keyboard.left] = true
                end
                acc += 1
                if l > lb
                    !adapt && error("Tuning parameter `c` too small. l/lb = $(l/lb)")
                    acc = num = 0
                    adapt!(c, i, factor)
                end
                θ = ZZB.reflect!(i, ∇ϕx[i], x, θ, Flow)
                b = ab(x, θ, c, Flow)
                told = t
                t′ = t .+ poisson_times(b)
                tfrez = freezing_time!(tfrez, t, x, θ, f, Flow)
   #             if !inbounds(game, x)
    #                game.score[] = max(0, game.score[]-1)
     #           end
            else
                b = ab(x, θ, c, Flow)
                told = t
                t′ = t .+ poisson_times(b)
                finished || continue
            end
        else #   t′ usual bouncy particle / boomerang step
            ∇ϕx = ∇ϕ!(∇ϕx, x, args...)
            ∇ϕx = grad_correct!(∇ϕx, x, Flow)
            l, lb = λ(∇ϕx, θ, Flow), sλ̄(b, t - told) # CHECK if depends on f
            num += 1
            if ((game.auto[] || !inbounds(game, x)) && rand()*lb <= l ) ||  (ispressed(canvas.scene, Keyboard.space) && !game.pressed[Keyboard.space]) # reflect!
                if ispressed(canvas.scene, Keyboard.space)
                    game.pressed[Keyboard.space] = true
                end
                acc += 1
                if l > lb
                    !adapt && error("Tuning parameter `c` too small.")
                    c *= factor
                end
                θ = reflect_sticky!(∇ϕx, x, θ, f, Flow)
                b = ab(x, θ, c, Flow) # regenerate reflection time
                told = t
                t′ = t + poisson_times(b)
                tfrez = freezing_time!(tfrez, t, x, θ, f, Flow)

#                if !inbounds(game, x)
 #                   game.score[] = max(0, game.score[]-1)
  #              end
            else # nothing happened
                b = ab(x, θ, c, Flow)
                told = t
                t′ = t + poisson_times(b)
                finished || continue
            end
        end

        push!(Ξ, sevent(t, x, θ, f, Flow))
        return finished, t, x, θ, t′, tref, tfrez, told, f, θf, (acc, num), c, b
    end

end
ZZB.Trace(t0, x0, θ0, f, Flow::ZigZag) = ZZB.Trace(t0, x0, θ0, f, BouncyParticle(Flow.Γ, Flow.μ, 0.0))
function ZZB.sevent(t, x, θ, f, Z::ZigZag)
    t, copy(x), copy(θ), copy(f)
end
ZZB.freezing_time!(tfrez, t0, x0, θ0, f, Flow::ZigZag) = ZZB.freezing_time!(tfrez, t0, x0, θ0, f,  BouncyParticle(Flow.Γ, Flow.μ, 0.0))
function ZZB.ab(x, θ, c, Z::ZigZag)
    a = [c + (ZZB.idot(Z.Γ, i, x)  - ZZB.idot(Z.Γ, i, Z.μ))'*θ[i] for i in eachindex(x)]
    b = [c + θ[i]'*ZZB.idot(Z.Γ, i, θ) for i in eachindex(x)]
    a, b
end
poisson_times((a, b)::Tuple{Array{Float64,1},Array{Float64,1}}) = ZZB.poisson_time.(a, b, rand(length(a)))
poisson_times((a, b)) = ZZB.poisson_time(a, b, rand())

function bouncy(game, X, ∇ϕ!, t0, x0, θ0, T, c, Flow::Union{ZigZag, BouncyParticle, Boomerang},
    κ, args...;  strong_upperbounds = false, adapt=false, factor=2.0)
    factorized = Flow isa ZigZag
    X[] .= [point(x0)]
    t, x, θ, ∇ϕx = t0, deepcopy(x0), deepcopy(θ0), deepcopy(θ0)
    told = t0
    θf = 0*θ # tags
    f = [true for _ in eachindex(x)]
    Ξ = Trace(t0, x0, θ0, f, Flow)
    push!(Ξ, sevent(t, x0, θ0, f, Flow))
    tref = waiting_time_ref(Flow) #refreshment times
    tfrez = zero(x)
    tfrez = freezing_time!(tfrez, t0, x0, θ0, f, Flow) #freexing times
    num = acc = 0
    b = ab(x, θ, c, Flow)
    t′ = t .+ poisson_times(b) # reflection time
    finished = false
    while !finished && t < T  && (!game.kill_on_boundary || inbounds(game, 0.9x))
        finished, t, x, θ, t′, tref, tfrez, told, f, θf, (acc, num), c, b = bouncy_inner!(game, X, Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, b, t′, f, θf, tfrez, tref, told, (acc, num),
                Flow, κ, args...; strong_upperbounds = strong_upperbounds, factor = factor, adapt = adapt)
        game.timeleft[] = T - t
    end
    return Ξ, (t, x, θ), (acc, num), c
end

phi(x, y, rho = 0.0) =  1/(2*pi*sqrt(1-rho^2))*exp(-0.5*(x^2 + y^2 - 2x*y*rho)/(1-rho^2))
sphi(x, y, σ1, σ2 = σ1, ρ=0.0) = phi(x/σ1, y/σ2, ρ)/(σ1*σ2)
logdensity1(x, y) = log(sphi(x, y+0.2x, 0.9))
logdensity2(x, y) = log(0.5sphi(x - dist, y - dist, 0.6) + 0.5sphi(x + dist, y + dist, 0.6))
#potential(xy) = -0.1exp(logdensity(xy[1], xy[2])-log(phi(xy[1], xy[2])))
getpotential(logdensity, boom=false) = function (xy)
     -(2-1.75boom)*exp(logdensity(xy[1], xy[2]) - boom*log(phi(xy[1], xy[2])))
end
getpotentialxy(logdensity, boom=false) = function (x,y)
    -(2-1.75boom)*exp(logdensity(x,y) - boom*log(phi(x, y)))
end

potentials = [getpotentialxy(logdensity1), getpotentialxy(logdensity2)]
potential(xy) = POTENTIAL[](xy[1], xy[2])
#potential(x,y) = POTENTIAL[](x, x)

gradϕ(logdensity) = function (x,i)
    -ForwardDiff.partials(logdensity(ForwardDiff.Dual{}(x[1], 1.0*(i==1)), ForwardDiff.Dual{}(x[2], 1.0*(i==2))))[]
end

gradϕ!(logdensity, g = gradϕ(logdensity)) = function (y, x)
    y[1] = g(x,1)
    y[2] = g(x,2)
    y
end
function remove(game, x)
    ys = game.ys[]
    for i in eachindex(ys)
        if norm(x - ys[i]) < 0.15
            deleteat!(ys, i)
            #Makie.notify(game.ys)
            game.ys[] = game.ys[]
            length(ys) == 0 && return true
            game.score[] = game.score[] + 1
            return false
        end
    end
    return false
end

point(x) = Point3f0(x[1], x[2], 0.03+potential(x))
inbounds(x) = -R < x[1] < R && -R < x[2] < R

point(potential, x) = Point3f0(x[1], x[2], 0.03+potential(x[1],x[2]))
point(i::Int, x) = Point3f0(x[1], x[2], 0.03+potentials[i](x[1],x[2]))
inbounds(game, x) = -R < x[1] < R && -R < x[2] < R

κ = [100.0, 100.0]
t0 = 0.0
x0 = [1.0, 0.]
θ0 = [1.0, 1.0]
c = 100.0
T = 1000.0

trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ!(logdensity2), t0, x0, θ0, T, c, BouncyParticle(sparse(I(d)), 0*x0, 0.1), 0.01*κ; adapt=false)
ts, xs = sep(collect(discretize(trace, T/500)))

score = Node(0)
SCORE = lift(x->"Score "*string(x), score)
speed = Node(norm(θ0))
timeleft = Node(T)
SPEED = lift(x->"Speed "*string(round(x, digits=2)), speed)
TIMELEFT = lift(y->"Time left "*string(round(y, digits=0)),  timeleft)
auto = Node(false)
AUTO = lift(x -> x ? "Autopilot on" : "Manual", auto)
MESSAGE = Node(" "^10*"Get ready - SPACE to reflect"*" "^10*"\n Q to quit")
ys = Node([rand(2)])
POTENTIAL = Node{Any}(getpotentialxy(logdensity1))
X = Node([point(x0)])
M = lift(p -> [p(x, y) for x in r, y in r], POTENTIAL)
YS = lift((i, ys)->point.(i, ys), POTENTIAL, ys)
ELEMENTS = [TIMELEFT, SPEED, SCORE, AUTO]
MARKER = Node(:circle)
TITLE = Node(" "^20*"Bouncy Particle Sampler Game"*" "^20)

canvas = Figure(resolution=(1500,1000))
#canvas[2:2,1] =  ax0 = Axis(canvas, aspect=DataAspect())
#hidedecorations!(ax0)
canvas[3:10,1] =  ax1 = Axis(canvas, aspect=DataAspect())
hidedecorations!(ax1)
hidespines!(ax1)
canvas[1:2,1:5] =  ax = Axis(canvas, title = "RETURN - continue to next lvl, SPACE - reflect, Q - quit, S* - stick, C* - Crank-Nicolson, A - auto pilot, +-* - speed (*hold)", aspect=DataAspect())
text!(ax, TITLE, align=(:center, :bottom), textsize=.85)
text!(ax, MESSAGE, color=RGB(0.1, 0.3, 0.7), align=(:center, :bottom), textsize=0.65, position=(-0, -2.2))
hidedecorations!(ax)
hidespines!(ax)
lscene = LScene(canvas[3:10, 2:5], scenekw = (camera = cam3d!, raw = false))

using Colors
#surface!(lscene, r, r, [potential((x1, x2)) for x1 in r, x2 in r], colormap=:atlantic, show_axis=false)
contour3d!(lscene, r, r, M, levels=30, linewidth=3.0)

scatter!(lscene, X, marker=MARKER, color=:red, markersize=200)
meshscatter!(lscene, YS, color=:orange, markersize=0.05)
#text!(lscene.scene, SCORE, position=(1.5R, 0R, 1), strokecolor=:black, strokewidth=3.0, color=RGB(0.6, 0.2, 0.0), textsize=1.0, rotation= .75pi)
for i in eachindex(ELEMENTS)
    text!(ax1, ELEMENTS[i], strokecolor=:black, strokewidth=.0, color=RGB(0.6, 0.2, 0.0), position=(0, -2*i), textsize=2.0)
end

display(canvas)

println("Find your Makie window and reflect with the SPACE key, quit with Q")

level = [
    (courage="Good start!", title="Bouncy particle sampler (BPS)", message="Stear the particle through the posterior landscape and collect the sample. \nChange direction with a 'contour reflection' at the countour lines with SPACE. ", logdensity=logdensity1, F=BouncyParticle(sparse(I(d)), 0*x0, 2.0, 0.995), x0=[0.1,1.0], θ0=[1.0,0.3], κ = [100.0, 100.0], T=100.0, ys=[[2.0,2.0]])
    (courage="Improving your ESS!", title="Bouncy particle sampler 2", message="Collect some more posterior points. Randomize your energy to reach\n the center with the Crank-Nicolson scheme (Hold C and reflect with SPACE).", logdensity=logdensity1, F=BouncyParticle(sparse(I(d)), 0*x0, 2.0, 0.995), x0=[0.1,1.0], θ0=[1.0,0.3], T=100.0, ys=0.9*[randn(2) for _ in 1:50])
    (courage="Great mixing!",title="Bouncy particle sampler 3", message="A double well potential! Collect the points.\nContinue with RETURN or try the `autopilot' (the original sampler) with A.", logdensity=logdensity2, F=BouncyParticle(sparse(I(d)), 0*x0, 2.0, 0.995), x0=[0.1,1.0], θ0=[1.0,0.3], κ = [100.0, 100.0], T=100.0, ys=[[0.6*(randn(2)-[dist,dist]) for _ in 1:50];[ 0.6*(randn(2) + [dist,dist]) for _ in 1:50]])

    (courage="Too easy!", title="Zig-Zag", message="Try to collect the point with the Zig-Zag.\n It always moves on the diagonals. Change direction with the arrow keys. \n Continue with RETURN.", logdensity=logdensity1, F=ZigZag(sparse(I(d)), 0*x0), x0=[0.1,1.0], θ0=[1.0,1.0], κ = [100.0, 100.0], T=100.0, ys=[[2.0,2.0]])
    (courage="Almost ergodic!", title="Zig-Zag 2", message="A double well potential! Collect the points. Try also A(utopilot) to see\n how the actual Zig-Zag sampler explores this posterior landscape.", logdensity=logdensity2, F=ZigZag(sparse(I(d)), 0*x0),x0=[0.1,1.0], θ0=[1.0,1.0], κ = [100.0, 100.0], T=1000.0, ys=[[0.6*(randn(2)-[dist,dist]) for _ in 1:50];[ 0.6*(randn(2) + [dist,dist]) for _ in 1:50]])

    (courage="So Hamiltonian!", title="Boomerang", message="The Boomerang moves on circles and does 'contour reflections' with SPACE.\n Change your energy level with +/- to reach the center and the outskirts.", logdensity=logdensity2, F=Boomerang(sparse(I(d)), 0*x0, 2.0, 0.995), x0=[0.1,1.0], θ0=[1.0,0.5], κ = [100.0, 100.0], T=100.0, ys=[[1.4,1.4], [0.1,-.1]])
    (courage="Done. You achieved mastership!", title="Sticky Bouncy particle sampler", message="Some of the points are on the coordinate axes.\nSTICK to the axis with S to collect them all (release to unstick).", logdensity=logdensity2, F=BouncyParticle(sparse(I(d)), 0*x0, 2.0, 0.995), x0=[0.1,1.0], θ0=[1.0,0.5], κ = [100.0, 100.0],   T=100.0, ys=[[0.6*rand(Bool,2).*(randn(2)-[dist,dist]) for _ in 1:50];[ 0.6*rand(Bool,2).*(randn(2) + [dist,dist]) for _ in 1:50]])

    ]
sleep(1.0)
for i in 1:length(level)
    lvl = level[i]
    ∇ϕ! = gradϕ!(lvl.logdensity)
    x0 = lvl.x0
    θ0 = lvl.θ0
    ys[] = lvl.ys
    auto[] = false
    MESSAGE[] = lvl.message
    sleep(2.)
    game = Game(false, Dict(Keyboard.space => false, Keyboard.s => false, Keyboard.left => false, Keyboard.right => false, Keyboard.a => false), Dict(Keyboard.s => false), ys, canvas, score, speed, timeleft, auto, getpotential(lvl.logdensity))
    POTENTIAL[] = getpotentialxy(lvl.logdensity, lvl.F isa Boomerang)
    MARKER[] = :circle
    #sleep(1.5)

    TITLE[] = lvl.title

    trace, (tT, xT, θT), (acc, num) = bouncy(game, X, ∇ϕ!, t0, x0, θ0, lvl.T, c, lvl.F, κ; adapt=false)

   # MARKER[] = :xcross

    MESSAGE[] = "$(lvl.courage) \nCurrent score: $(score[]) posterior\nsamples."
    sleep(1.5)
end


text:
using NPZ
using Plots
using NearestNeighbors

include("io.jl")
include("functions.jl")
include("plot_utils.jl")
include("voronoi_utils.jl")
include("characteristics.jl")
include("irregular_ray_tracing.jl")

pyplot()

function searchlight_irregular()
    nx = ny = nz = 51

    n_sites = nz*nx*ny

    bounds = [[0 1]
              [0 1]
              [0 1]]

    temperature = ones(n_sites)*1u"K"
    electron_density = zeros(n_sites)*1u"m^-3"
    hydrogen_density = zeros(n_sites)*1u"m^-3"

    velocity_z = zeros(n_sites)u"m/s"
    velocity_x = zeros(n_sites)u"m/s"
    velocity_y = zeros(n_sites)u"m/s"

    println("---Computing grid---")
    positions = rand(3, n_sites)u"m"

    v0 = [0, 0.5, 0.5]
    R0 = 0.1u"m"

    n_sweeps = 3
    # positions = sample_beam(n_sites, bounds, beam, v0, R0, k)

    sites_file = "../data/searchlight_sites.txt"
    neighbours_file = "../data/searchlight_neighbours.txt"

    # write sites to file
    write_arrays(ustrip.(positions[2, :]),
                 ustrip.(positions[3, :]),
                 ustrip.(positions[1, :]),
                 sites_file)

    # compute neigbours
    run(`./voro.sh $sites_file $neighbours_file
            $(bounds[2,1]) $(bounds[2,2])
            $(bounds[3,1]) $(bounds[3,2])
            $(bounds[1,1]) $(bounds[1,2])`)

    # Voronoi grid
    sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                         temperature, electron_density, hydrogen_density,
                         velocity_z, velocity_x, velocity_y,
                         bounds[1,1]*1u"m", bounds[1,2]*1u"m",
                         bounds[2,1]*1u"m", bounds[2,2]*1u"m",
                         bounds[3,1]*1u"m", bounds[3,2]*1u"m",
                         n_sites)

    S = zeros(n_sites)u"kW*m^-2*nm^-1"
    α = zeros(n_sites)u"m^-1"

    tree = KDTree(ustrip(sites.positions))

    I_light = 1u"kW*m^-2*nm^-1"

    bottom_layer = sites.layers_up[2] - 1

    I_bottom = zeros(bottom_layer)u"kW*m^-2*nm^-1"
    for i in 1:bottom_layer
        idx = sites.perm_up[i]
        xi = sites.positions[2, idx]
        yi = sites.positions[3, idx]
        if sqrt((xi - 0.5u"m")^2 + (yi - 0.5u"m")^2) < R0
            I_bottom[i] = I_light
        end
    end

    top_layer = sites.layers_down[2] - 1

    I_top = zeros(top_layer)u"kW*m^-2*nm^-1"
    for i in 1:top_layer
        idx = sites.perm_down[i]
        xi = sites.positions[2, idx]
        yi = sites.positions[3, idx]
        if sqrt((xi - 0.5u"m")^2 + (yi - 0.5u"m")^2) < R0
            I_top[i] = I_light
        end
    end

    x = collect(LinRange(0,1,10*nx))
    y = collect(LinRange(0,1,10*ny))
    top_z = 1.0
    bottom_z = 0.0

    RES=500

    tot_time = 0.0

    weights, θ_array, ϕ_array, n_angles = read_quadrature("../quadratures/ul7n12.dat")

    for i in eachindex(θ_array)
        θ = θ_array[i]
        ϕ = ϕ_array[i]
        println("$(floor(Int,θ)), $(floor(Int,ϕ))")
        # Unit vector pointing in the direction of the ray
        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
        if θ > 90
            I, time = @timed Delaunay_upII(k, S, I_bottom, α, sites, n_sweeps)

            #=
            top_I = zeros(length(x), length(y))u"kW*nm^-1*m^-2"

            for i in 1:length(x)
                for j in 1:length(y)
                    position = [top_z, x[i], y[j]]
                    idx, dist = nn(tree, ustrip(position))
                    top_I[i, j] = I[idx]
                end
            end

            plot_searchlight(k, x, y, top_I, R0, "irregular_$(floor(Int,θ))_$(floor(Int,ϕ))")
            =#
        elseif θ < 90
            I, time = @timed Delaunay_downII(k, S, I_top, α, sites, n_sweeps)

            #=
            bottom_z = 0
            bottom_I = zeros(length(x), length(y))u"kW*nm^-1*m^-2"

            for i in 1:length(x)
                for j in 1:length(y)
                    position = [bottom_z, x[i], y[j]]
                    idx, dist = nn(tree, ustrip(position))
                    bottom_I[i, j] = I[idx]
                end
            end

            # npzwrite("../data/searchlight_data/I_$(θ)_$(ϕ)_voronoi.npy", ustrip.(bottom_I))
            # plot_searchlight(k, x, y, bottom_I, R0, "irregular_$(floor(Int,θ))_$(floor(Int,ϕ))")
            =#
        end
        tot_time += time
        # npzwrite("../data/searchlight_data/x_voronoi.npy", x)
        # npzwrite("../data/searchlight_data/y_voronoi.npy", y)
    end
    println("Total time $tot_time -- avg. time: $(tot_time/12)")
end

function searchlight_regular()
    nx = ny = nz = 51

    z = collect(LinRange(0,1,nz))u"m"
    x = collect(LinRange(0,1,nx))u"m"
    y = collect(LinRange(0,1,ny))u"m"

    # npzwrite("../data/searchlight_data/x_regular.npy", ustrip.(x))
    # npzwrite("../data/searchlight_data/y_regular.npy", ustrip.(y))

    temperature = ones(nz, nx, ny)u"K"
    electron_density = zeros(nz, nx, ny)u"m^-3"
    hydrogen_density = zeros(nz, nx, ny)u"m^-3"

    velocity_z = zeros(nz, nx, ny)u"m/s"
    velocity_x = zeros(nz, nx, ny)u"m/s"
    velocity_y = zeros(nz, nx, ny)u"m/s"

    atmos = Atmosphere(z, x, y, temperature, electron_density, hydrogen_density,
                       velocity_z, velocity_x, velocity_y)

    S_0 = zeros(nz, nx, ny)u"kW*m^-2*nm^-1"
    α = zeros(nz, nx, ny)u"m^-1"

    I_light = 1u"kW*m^-2*nm^-1"

    I_0 = zero(S_0[1,:,:])
    R0 = 0.1
    for i in 1:nx
        for j in 1:ny
            xi = i/nx
            yi = j/ny
            if sqrt((xi - 0.5)^2 + (yi - 0.5)^2) < R0
                I_0[i, j] = I_light
            end
        end
    end

    weights, θ_array, ϕ_array, n_angles = read_quadrature("../quadratures/ul7n12.dat")

    tot_time = 0.0

    for i in eachindex(θ_array)
        θ = θ_array[i]
        ϕ = ϕ_array[i]
        println("$(floor(Int,θ)), $(floor(Int,ϕ))")
        # Unit vector pointing in the direction of the ray
        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
        if θ > 90
            I, time = @timed short_characteristics_up(k, S_0, I_0, α, atmos;
                                         pt=true, n_sweeps=3)[:, 2:end-1, 2:end-1]

            I = I[end, :, :]
            # npzwrite("../data/searchlight_data/I_$(θ)_$(ϕ)_regular.npy", ustrip.(I))
            plot_searchlight(k, x[2:end-1], y[2:end-1], I, R0, "regular_$(floor(Int,θ))_$(floor(Int,ϕ))")
            println("Bottom: $(I_light*80), Top: $(sum(I))")
        elseif θ < 90
            I, time = @timed short_characteristics_down(k, S_0, I_0, α, atmos;
                                           pt=true, n_sweeps=3)[:, 2:end-1, 2:end-1]

            I = I[1, :, :]
            # npzwrite("../data/searchlight_data/I_$(θ)_$(ϕ)_regular.npy", ustrip.(I))
            plot_searchlight(k, x[2:end-1], y[2:end-1], I, R0, "regular_$(floor(Int,θ))_$(floor(Int,ϕ))")
            println("Top: $(I_light*80), Bottom: $(sum(I))")
        end
        tot_time += time
    end

    println("Total time $tot_time -- avg. time: $(tot_time/12)")

    print("")
end

searchlight_irregular()
searchlight_regular()

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
function _prepare_workspace!(ref::SimulationRef, base_name::AbstractString, folder::AbstractString)

    !isdir(folder) && error("Specified folder is not valid")
    global_path = joinpath(folder, "$(base_name)")
    !isdir(global_path) && mkpath(global_path)
    _sim_path = replace_chars("$(round(Dates.now(),Dates.Minute))", ":", "-")
    simulation_path = joinpath(global_path, _sim_path)
    raw_ouput = joinpath(simulation_path, "raw_output")
    mkpath(raw_ouput)
    models_json_ouput = joinpath(simulation_path, "models_json")
    mkpath(models_json_ouput)
    results_path = joinpath(simulation_path, "results")
    mkpath(results_path)

    ref.raw = raw_ouput
    ref.models = models_json_ouput
    ref.results = results_path

    return

end

function _validate_steps(stages::Dict{Int64, Stage},
                         steps::Int64,
                         stage_initial_times::Dict{Int64, Vector{Dates.DateTime}})

    for (stage_number, s) in stages
        forecast_count = length(stage_initial_times[stage_number])

        if steps*s.execution_count > forecast_count #checks that there are enough time series to run
            error("The number of available time series is not enough to perform the
                   desired amount of simulation steps.")
        end

    end

    return

end

function _get_dates(stages::Dict{Int64, Stage})
    k = keys(stages)
    k_size = length(k)
    range = Vector{Dates.DateTime}(undef, 2)
    @assert k_size == maximum(k)

    for i in 1:k_size
        initial_times = PSY.get_forecast_initial_times(stages[i].sys)
        i == 1 && (range[1] = initial_times[1])
        interval = PSY.get_forecasts_interval(stages[i].sys)
        for (ix,element) in enumerate(initial_times[1:end-1])
            if !(element + interval == initial_times[ix+1])
                error("The sequence of forecasts is invalid")
            end
        end
        (i == k_size && (range[end] = initial_times[end]))
    end

    return Tuple(range), true

end

function _populate_cache!(stage::_Stage)

    for (k, cache) in stage.cache
        build_cache!(cache, stage.canonical)
    end

    return
end

function _build_stages(sim_ref::SimulationRef,
                       stages::Dict{Int64, Stage},
                       verbose::Bool = true;
                       kwargs...)

    system_to_file = get(kwargs, :system_to_file, true)
    mod_stages = Vector{_Stage}(undef, length(stages))
    for (key, stage) in stages
        verbose && @info("Building Stage $(key)")
        canonical = Canonical(stage.model.transmission,
                                   stage.sys,
                                   stage.optimizer;
                                   use_parameters = true,
                                   initial_time = stage.initial_time,
                                   horizon = stage.horizon)
        mod_stages[key] = _Stage(key,
                                stage.model,
                                stage.op_problem,
                                stage.sys,
                                canonical,
                                stage.optimizer,
                                stage.execution_count,
                                stage.interval,
                                stage.chronology_ref,
                                stage.cache)
        _build!(mod_stages[key].canonical,
                stage.model,
                stage.sys;
                kwargs...)
        stage_path = joinpath(sim_ref.models,"stage_$(key)_model")
        mkpath(stage_path)
        _write_canonical(canonical, joinpath(stage_path, "optimization_model.json"))
        system_to_file && IS.to_json(stage.sys, joinpath(stage_path ,"sys_data.json"))
        _populate_cache!(mod_stages[key])
        sim_ref.date_ref[key] = PSY.get_forecast_initial_times(stage.sys)[1]
    end

    return mod_stages

end

function _feedforward_rule_check(::Type{T},
                              stage_number_from::Int64,
                              from_stage::Stage,
                              stage_number_to::Int64,
                              to_stage::Stage,) where T <: Chronology

    error("Feedforward Model $(T) not implemented")

    return

end

function _feedforward_rule_check(synch::Synchronize,
                                 stage_number_from::Int64,
                                 from_stage::Stage,
                                 stage_number_to::Int64,
                                 to_stage::Stage)

    #Don't check for same Stage.
    stage_number_from == stage_number_to && return

    from_stage_horizon = PSY.get_forecasts_horizon(from_stage.sys)
    to_stage_count = get_execution_count(to_stage)
    to_stage_synch = synch.to_steps
    from_stage_synch = synch.from_horizon

    if from_stage_synch > from_stage_horizon
        error("The lookahead length $(from_stage_horizon) in stage is insufficient to synchronize with $(from_stage_synch) feedforward steps")
    end

    if to_stage_synch*from_stage_synch != to_stage_count
        error("The execution total in stage is inconsistent with a chronology
                of $(from_stage_synch) feedforward steps and $(to_stage_synch) runs. The expected
                number of executions is $(to_stage_synch*from_stage_synch)")
    end

    if (from_stage_horizon % from_stage_synch) != 0
        error("The number of feedforward steps $(from_stage_horizon) in stage
               needs to be a mutiple of the horizon length $(from_stage_horizon)
               of stage to use Synchronize with parameters ($(from_stage_synch), $(to_stage_synch))")
    end

    return

end

function _feedforward_rule_check(sync::Sequential,
                              stage_number_from::Int64,
                              from_stage::Stage,
                              stage_number_to::Int64,
                              to_stage::Stage)

    return

end

function _feedforward_rule_check(sync::RecedingHorizon,
                                stage_number_from::Int64,
                                from_stage::Stage,
                                stage_number_to::Int64,
                                to_stage::Stage)

    return

end

function _check_chronology_ref(stages::Dict{Int64, Stage})

    for (stage_number,stage) in stages
        for (key, chron) in stage.chronology_ref
            key < 1 && continue
            _feedforward_rule_check(chron, key, stages[key], stage_number, stage)
        end
    end

    return

end

function _build_simulation!(sim_ref::SimulationRef,
                            steps::Int64,
                            stages::Dict{Int64, Stage};
                            verbose::Bool = false, kwargs...)

    stage_initial_times = Dict{Int64, Vector{Dates.DateTime}}()
    for (stage_number, stage) in stages
        PSY.check_forecast_consistency(stage.sys)
        if PSY.are_forecasts_contiguous(stage.sys)
            stage_initial_times[stage_number] = PSY.generate_initial_times(stage.sys,
                                                                        stage.interval,
                                                                        stage.horizon;
                                                                        initial_time = stage.initial_time)
        else
            stage_initial_times[stage_number] = get_forecast_initial_times(stage.sys)
        end
    end
    _validate_steps(stages, steps, stage_initial_times)
    _check_chronology_ref(stages)
    dates, validation = _get_dates(stages)
    return dates, validation, _build_stages(sim_ref, stages, verbose = verbose; kwargs...)

end


text:
NonlinearEigenproblems.contour_beyn(nep::NEP,display::Bool;params...)=contour_beyn(ComplexF64,nep,display;params...)

# the same as NonlinearEigenproblems.contour_beyn but compatible with ProgressMeter, notice extra display argument to distinguish from main method
function NonlinearEigenproblems.contour_beyn(::Type{T},
            nep::NEP,
            display::Bool;
            tol::Real=sqrt(eps(real(T))), # Note tol is quite high for this method
            σ::Number=zero(complex(T)),
            displaylevel::Integer=0,
            linsolvercreator::Function=backslash_linsolvercreator,
            neigs::Integer=2, # Number of wanted eigvals
            k::Integer=neigs+1, # Columns in matrix to integrate
            radius::Union{Real,Tuple,Array}=1, # integration radius
            quad_method::Symbol=:ptrapz, # which method to run. :quadg, :quadg_parallel, :quadgk, :ptrapz
            N::Integer=1000,  # Nof quadrature nodes
            # errmeasure::Function = NonlinearEigenproblems.default_errmeasure(nep::NEP),
            sanity_check=true,
            rank_drop_tol=tol # Used in sanity checking
            )where{T<:Number}

    # Geometry
    length(radius)==1 ? radius=(radius,radius) : nothing
    g(t) = complex(radius[1]*cos(t),radius[2]*sin(t)) # ellipse
    gp(t) = complex(-radius[1]*sin(t),radius[2]*cos(t)) # derivative

    n=size(nep,1);

    if (k>n)
        error("Cannot compute more eigenvalues than the size of the NEP with contour_beyn() k=",k," n=",n);

    end
    if (k<=0)
        error("k must be positive, k=",k,
              neigs==typemax(Int) ? ". The kwarg k must be set if you use neigs=typemax" : ".")
    end

    Random.seed!(10); # Reproducability
    Vh=Array{T,2}(randn(real(T),n,k)) # randn only works for real

    function local_linsolve(λ::TT,V::Matrix{TT}) where {TT<:Number}
        # @ifd(print("."))
        local M0inv::LinSolver = linsolvercreator(nep,λ+σ);
        # This requires that lin_solve can handle rectangular
        # matrices as the RHS
        return lin_solve(M0inv,V);
    end

    # Constructing integrands
    Tv(λ) = local_linsolve(T(λ),Vh)
    f(t) = Tv(g(t))*gp(t)
    # @ifd(print("Computing integrals"))

    pg = Progress(N; dt=.1, desc="Contour integration...")
    local A0,A1
    if (quad_method == :quadg_parallel)
        # @ifd(print(" using quadg_parallel"))
        error("disabled");
    elseif (quad_method == :quadg)
        # @ifd(print(" using quadg"))
        error("disabled");
    elseif (quad_method == :ptrapz)
        # @ifd(print(" using ptrapz"))
        (A0,A1)=ptrapz(f,g,0,2*pi,N,pg);
    elseif (quad_method == :ptrapz_parallel)
        # @ifd(print(" using ptrapz_parallel"))
        channel = RemoteChannel(()->Channel{Bool}(pg.n+1),1)
        @sync begin
            @async begin
                (A0,A1)=ptrapz_parallel(f,g,0,2*pi,N,channel);
                put!(channel,false)
            end
            @async while take!(channel)
                next!(pg)
            end
        end
    elseif (quad_method == :quadgk)
        error("disabled");
    else
        error("Unknown quadrature method:"*String(quad_method));
    end

    # Don't forget scaling
    A0[:,:] = A0 ./(2im*pi);
    A1[:,:] = A1 ./(2im*pi);

    # @ifd(print("Computing SVD prepare for eigenvalue extraction "))
    V,S,W = svd(A0)
    p = count( S/S[1] .> rank_drop_tol);
    # @ifd(println(" p=",p));

    V0 = V[:,1:p]
    W0 = W[:,1:p]
    B = (copy(V0')*A1*W0) * Diagonal(1 ./ S[1:p])

    # Extract eigenval and eigvec approximations according to
    # step 6 on page 3849 in the reference
    # @ifd(println("Computing eigenvalues "))
    λ,VB=eigen(B)
    λ[:] = λ .+ σ
    # @ifd(println("Computing eigenvectors "))
    V = V0 * VB;
    for i = 1:p
        normalize!(V[:,i]);
    end

    if (!sanity_check)
        sorted_index = sortperm(map(x->abs(σ-x), λ));
        inside_bool = (real(λ[sorted_index].-σ)/radius[1]).^2 + (imag(λ[sorted_index].-σ)/radius[2]).^2 .≤ 1
        inside_perm = sortperm(.!inside_bool)
        return (λ[sorted_index[inside_perm]],V[:,sorted_index[inside_perm]])
    end

    # Compute all the errors
    errmeasures=zeros(real(T),p);
    # for i = 1:p
        # errmeasures[i]=errmeasure(λ[i],V[:,i]);
    # end

    good_index=findall(errmeasures .< tol);

    # Index vector for sorted to distance from σ
    sorted_good_index=
       good_index[sortperm(map(x->abs(σ-x), λ[good_index]))];

    # check that eigenvalues are inside contour, move them to the end if they are outside
   inside_bool = (real(λ[sorted_good_index].-σ)/radius[1]).^2 + (imag(λ[sorted_good_index].-σ)/radius[2]).^2 .≤ 1
   if any(.!inside_bool)
       if neigs ≤ sum(inside_bool)
           @warn "found $(sum(.!inside_bool)) evals outside contour, $p inside. all $neigs returned evals inside, but possibly inaccurate. try increasing N, decreasing tol, or changing radius"
       else
           @warn "found $(sum(.!inside_bool)) evals outside contour, $p inside. last $(neigs-sum(inside_bool)) returned evals outside contour. possible inaccuracy. try increasing N, decreasing tol, or changing radius"
       end
   end
   sorted_good_inside_perm = sortperm(.!inside_bool)

    # Remove all eigpairs not sufficiently accurate
    # and potentially remove eigenvalues if more than neigs.
    local Vgood,λgood
    if( size(sorted_good_index,1) > neigs)
        # @ifd(println("Removing unwanted eigvals: neigs=",neigs,"<",size(sorted_good_index,1),"=found_eigvals"))
        Vgood=V[:,sorted_good_index[sorted_good_inside_perm][1:neigs]];
        λgood=λ[sorted_good_index[sorted_good_inside_perm][1:neigs]];
    else
        Vgood=V[:,sorted_good_index[sorted_good_inside_perm]];
        λgood=λ[sorted_good_index[sorted_good_inside_perm]];
    end

    if (p==k)
       @warn "Rank-drop not detected, your eigvals may be correct, but the algorithm cannot verify. Try to increase k. This warning can be disabled with `sanity_check=false`." S
    end

    if (size(λgood,1)<neigs  && neigs < typemax(Int))
       @warn "We found fewer eigvals than requested. Try increasing domain, or decreasing `tol`. This warning can be disabled with `sanity_check=false`." S
    end

    return (λgood,Vgood)
end

# Trapezoidal rule for a periodic function f
function ptrapz(f,g,a,b,N,pg)
    h = (b-a)/N
    t = range(a, stop = b-h, length = N)
    S0 = zero(f(t[1])); S1 = zero(S0)
    for i = 1:N
        temp = f(t[i])
        S0 += temp
        S1 += temp*g(t[i])
        next!(pg)
    end
    return h*S0, h*S1;
end

function ptrapz_parallel(f,g,a,b,N,channel)
    h = (b-a)/N
    t = range(a, stop = b-h, length = N)
    S = @distributed (+) for i = 1:N
        temp = f(t[i])
        put!(channel,true)
        [temp,temp*g(t[i])]
    end
    return h*S[1], h*S[2];
end


text:
function dtd2uftau(t, fft_u, fft_v, fft_du, fft_dv,
                      fft_du2, fft_dv2,
                      A1, A2, matr, conjmatr,
                      sigma, llambda, Ntaumm)

    sigma = 1

    u     = ifft(exp(1im*t*A1) * fft_u * matr)
    v     = ifft(exp(1im*t*A1) * fft_v * matr)
    dtu   = ifft(exp(1im*t*A1) * (1im*A1) * fft_u) * matr
    dtv   = ifft(exp(1im*t*A1) * (1im*A1) * fft_v) * matr

    du    = ifft(exp(1im*t*A1) * fft_du * matr)
    dv    = ifft(exp(1im*t*A1) * fft_dv * matr)
    dtdu  = ifft(exp(1im*t*A1) * (1im*A1) * fft_du) * matr
    dtdv  = ifft(exp(1im*t*A1) * (1im*A1) * fft_dv) * matr

    du2   = ifft(exp(1im*t*A1) * fft_du2 * matr)
    dv2   = ifft(exp(1im*t*A1) * fft_dv2 * matr)
    dtdu2 = ifft(exp(1im*t*A1) * (1im*A1) * fft_du2) * matr
    dtdv2 = ifft(exp(1im*t*A1) * (1im*A1) * fft_dv2) * matr

    z     = (u+conjugate(v))/2
    dtz   = (dtu+conjugate(dtv))/2
    dz    = (du+conjugate(dv))/2
    dtdz  = (dtdu+conjugate(dtdv))/2
    dz2   = (du2+conjugate(dv2))/2
    dtdz2 = (dtdu2+conjugate(dtdv2))/2

    fz1 = (  2*dtdz2*conjugate(z)*dz
           + 2*dz2*conjugate(dtz)*dz
           + 2*dz2*conjugate(z)*dtdz
           + 2*dtz*conjugate(dz2)*dz
           + 2*z*conjugate(dtdz2)*dz
           + 2*z*conjugate(dz2)*dtdz
           + 2*dtz*dz2*conjugate(dz)
           + 2*z*dtdz2*conjugate(dz)
           + 2*z*dz2*conjugate(dtdz) )

    champu1 = -1im*llambda*A2*exp(-1im*t*A1)*fft(conjmatr*fz1)
    champv1 = -1im*llambda*A2*exp(-1im*t*A1)*fft(conjmatr*conjugate(fz1))

    fz1 = (  2*dz2*conjugate(z)*dz
           + 2*z*conjugate(dz2)*dz
           + 2*z*dz2*conjugate(dz))

    champu2 = -1im*llambda*A2*exp(-1im*t*A1)*(-1im*A1)*fft(conjmatr*fz1)
    champv2 = -1im*llambda*A2*exp(-1im*t*A1)*(-1im*A1)*fft(conjmatr*conjugate(fz1))

    champu=champu1+champu2
    champv=champv1+champv2

    champu,champv

end


text:
module EGM96

    if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
        include("../deps/deps.jl")
    else
        error("EGM96 not properly installed. Please run Pkg.build(\"EGM96\")")
    end

    export undulation
    include("f477.jl")
end


text:
### Utilities specifically targeted at finding functions that can
### trigger a GC event
@doc """
`gcnames = findgc(dirname=".")` parses the `*.dot` files in the specified
directory, assembles the combined callgraph, and then finds the names of all
functions which might, directly or indirectly, trigger a garbage-collection
event by calling `jl_gc_collect`.
""" ->
function findgc(dirname=".")
    cgs = parsedots(dirname)
    calls, calledby = combine(cgs...)
    callers = findcallers(calledby, "jl_gc_collect")
    union!(callers, findcallers(calledby, "jl_apply"))
    uniquenames([callers...])
end

@doc """
`highlight(srcfilename, gcnames)` displays the specified C source file in the REPL,
coloring in red all function calls which might trigger garbage collection.
In some cases, this can be helpful in determining whether one might be missing a GC root.
""" ->
function highlightgc(srcfilename, gcnames::Set)
    lines = open(srcfilename) do file
        readlines(file)
    end
    modline = IOBuffer()
    for line in lines
        indx = 1
        while true
            m = match(ralphanum, line, indx)
            if m == nothing
                print(modline, line[indx:end])
                break
            end
            print(modline, line[indx:m.offset-1])
            if m.match in gcnames
                print_with_color(:red, modline, m.match)
            else
                print(modline, m.match)
            end
            indx = m.offset + length(m.match)
        end
        print(takebuf_string(modline))
    end
end

@doc """
```
emacs_highlighting(filename, funcnames, face="font-lock-warning-face")
```
Creates an emacs highlighting file from a list of function names.
""" ->
function emacs_highlighting(filename, funcnames, face="font-lock-warning-face")
    open(filename, "w") do file
        print(file, "(defvar CallGraphFuncs\n  '((\"\\\\<\\\\(")
        first = true
        for fname in funcnames
            if !first
                print(file, "\\\\|")
            end
            print(file, fname)
            first = false
        end
        println(file, "\\\\)\\\\>\" . '$face)))")
        println(file, "\n(font-lock-add-keywords 'c-mode CallGraphFuncs)")
    end
end


text:
module C_block_diag

include(dirname(@__FILE__)*"\\Moments\\Moments.jl")
using .Moments

export  split_expo,
        get_all_partitions,
        get_all_blocks,
        get_coef,
        get_γγᶥδδᶥζζᶥηηᶥ,
        get_ℜℑααᶥββᶥᴿ,
        get_real_blocks

## The complex blocks
"""Splits a vector of exponents into the register components ααᶥββᶥ →  (α,αᶥ,β,βᶥ)"""
split_expo(ααᶥββᶥ,d) =     (ααᶥββᶥ[1:d[1]],
                            ααᶥββᶥ[d[1]+1:2*d[1]],
                            ααᶥββᶥ[1+2*d[1]:2*d[1]+d[2]],
                            ααᶥββᶥ[1+2*d[1]+d[2]:end])

#@test split_expo([1,1,1,2,2,2,3,3,3,4,4,4],(3,3)) == ([1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4])

"""Partitions Iᵗ such that M( Iᵗ , Iᵗ ) would be block diagonal"""
get_r(v) = v[1] - v[2]
get_s(v) = v[3] - v[4]
function get_all_partitions(d,t)
    Iᵗ = Moments.make_mon_expo(d,t[1])
    Iᵗ_temp = map(x ->  sum.(split_expo(x,d)),Iᵗ)

    r_list = map(v -> get_r(v),Iᵗ_temp)
    s_list = map(v -> get_s(v),Iᵗ_temp)

    Iᵗ_dict = Dict()
    Iᵗ_perm = []
    for r ∈ -t[1]:t[1], s ∈ -t[1]:t[1]
        T = Iᵗ[(r_list .== r) .& (s_list .== s) ]
        isempty(T) ? continue : Iᵗ_dict[r,s] = T
    end
    return Iᵗ_dict
end

"""Gets the principal block of the moment matrix after block diagonalization"""
get_all_blocks(Iᵗd) = Dict(zip(keys(Iᵗd),[Iᵗd[k] .+ reshape(Iᵗd[k],1,:) for k in keys(Iᵗd)]))

## The real analogue

"""Returns all  γγᶥδδᶥ,ζζᶥηηᶥ ∈ (ℕᵈ)⁴ s.t.
ααᶥββᶥ =  γγᶥδδᶥ,ζζᶥηηᶥ
∀  ααᶥββᶥ ∈ Moment matrixₜ
"""
function get_γγᶥδδᶥζζᶥηηᶥ(d,t)
    M = Moments.make_mon_expo(d,t)
    M_vec = M[:,1]
    lok = Moments.make_mon_expo(d,2*t[1]) # Define variables in the moment matrix.
    nara = Dict()
    for k in lok
        fd = hcat(map(x -> [x[1],x[2]],findall(M .== [k]))...)
        nara[k] = [M_vec[fd[1,:]], M_vec[fd[2,:]]]
    end
    return nara
end

"""The coefficients function"""
pf(x) = prod(factorial.(x))
function get_coef(γγᶥδδᶥ,ηηᶥζζᶥ,d)
    γ,γᶥ,δ,δᶥ = split_expo(γγᶥδδᶥ,d)
    η,ηᶥ,ζ,ζᶥ = split_expo(ηηᶥζζᶥ,d)
    a = ((-1)^sum(δᶥ + ζᶥ) * (im)^sum(δ+δᶥ+ζ+ζᶥ))* prod(pf.([γ+δ γᶥ+δᶥ η+ζ ηᶥ+ζᶥ]))
    return a/prod(pf.([γ,γᶥ,δ,δᶥ,η,ηᶥ,ζ,ζᶥ]))
end

"""f"""
function get_ℜℑααᶥββᶥᴿ(d,ααᶥββᶥ,γdict,K)
    Tar = γdict[ααᶥββᶥ]
    Res_vec = []
    Res_coef = []
    for k ∈ 1:K
        if k < length(Tar[1])
            γ_,η_ = Tar[1][k],Tar[2][k]
            push!(Res_coef,get_coef(γ_,η_,d))

            γ,γᶥ,δ,δᶥ = split_expo(γ_,d)
            η,ηᶥ,ζ,ζᶥ = split_expo(η_,d)
            push!(Res_vec,vcat(γ+γᶥ,δ+δᶥ,η+ηᶥ,ζ+ζᶥ))
        else
            push!(Res_coef,0.0)
            push!(Res_vec,Int.(zeros(length(ααᶥββᶥ))))
        end
    end
    return  Res_vec,Res_coef
end

"""Get the exponents and coefficients in matrix form"""
function get_real_block(d,B,γdic,K)
    Res_vec = Dict()
    Res_coef = Dict()
    for i in 1:size(B)[1], j in i:size(B)[1]
            Res_vec[i,j], Res_coef[i,j] = get_ℜℑααᶥββᶥᴿ(d,B[i,j],γdic,K)
            Res_vec[j,i]  = Res_vec[i,j]
            Res_coef[j,i] = Res_coef[i,j]
    end

    Expo_mat_dict = Dict()
    Coef_mat_dict = Dict()
    for k in 1:K
        Expo_mat_dict[k] = [Res_vec[i,j][k] for i = 1:size(B)[1], j = 1:size(B)[1]]
        Coef_mat_dict[k] = [Res_coef[i,j][k] for i = 1:size(B)[1], j = 1:size(B)[1]]
    end
    return clean_real_blocks(Coef_mat_dict, Expo_mat_dict)
end

function clean_real_blocks(Coef_mat_dict, Expo_mat_dict)
    for k in keys(Coef_mat_dict)
        if iszero(Coef_mat_dict[k]) ## Get rid of zeros
            delete!(Coef_mat_dict,k)
            delete!(Expo_mat_dict,k)
            continue
        end

        if isreal(Coef_mat_dict[k]) ## Reduce to real
            Coef_mat_dict[k] = real(Coef_mat_dict[k])
        # elseif imag(Coef_mat_dict[k]) == Coef_mat_dict[k]
        #     Coef_mat_dict[k] = imag(Coef_mat_dict[k])
        end
    end
    return Coef_mat_dict, Expo_mat_dict
end

function get_real_blocks(d,t)
    Iᵗ_dict = C_block_diag.get_all_partitions(d,t)
    block_dict = C_block_diag.get_all_blocks(Iᵗ_dict)
    γγᶥδδᶥζζᶥηηᶥ_dict = C_block_diag.get_γγᶥδδᶥζζᶥηηᶥ(d,t)
    K = maximum([length(x[1]) for x in values(γγᶥδδᶥζζᶥηηᶥ_dict)])
    Coef_mat_dict = Dict()
    Expo_mat_dict = Dict()
    for k in keys(block_dict)
        Coef_mat_dict[k], Expo_mat_dict[k] = C_block_diag.get_real_block(d,block_dict[k],γγᶥδδᶥζζᶥηηᶥ_dict,K)
    end
    return Coef_mat_dict, Expo_mat_dict
end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
include("segbase.jl")
include("plot.jl")


text:
# WHAT IF caps[i] = 0
module DA
    include("matchings.jl")

    export deferred_acceptance, generate_random_prefs
end


text:
function find_species_in_compartment(list_of_species,compartment_string)

    # initialize -
    index_array = Int64[]

    # main -
    counter = 1
    for species_symbol in list_of_species

        if (occursin(compartment_string,species_symbol) == true)
            push!(index_array,counter)
        end

        # update the counter -
        counter = counter + 1
    end

    # return -
    return index_array
end

function reconstruct_reaction_string_list(cobra_dictionary)

    # initialize -
    reaction_string_buffer = String[]

    # get the stoichiometric array -
    stoichiometric_matrix = Matrix(cobra_dictionary["S"])
    list_of_reaction_tags = cobra_dictionary["rxns"]
    list_of_species = cobra_dictionary["mets"]
    list_of_reversible_reactions = cobra_dictionary["rev"]

    # what is the size?
    (number_of_species,number_of_reactions) = size(stoichiometric_matrix)
    for reaction_index = 1:number_of_reactions

        # initialize empty buffer -
        buffer = ""

        # get the reaction tag -
        reaction_tag_string = list_of_reaction_tags[reaction_index]

        # add the tag to the buffer -
        buffer *= "$(reaction_index),$(reaction_tag_string),"

        # find the reactants -
        idx_reactants = findall(stoichiometric_matrix[:,reaction_index].<0.0)
        if (isempty(idx_reactants) == true)
            buffer *= "[],"
        else

            # how many species do we have?
            number_of_species = length(idx_reactants)
            counter = 1
            for index in idx_reactants

                # get the metabolite -
                metabolite_string = list_of_species[index]
                stcoeff = stoichiometric_matrix[index,reaction_index]

                if (stcoeff != -1.0)
                    # add data to the buffer -
                    buffer *= "$(abs(stcoeff))*$(metabolite_string)"
                else
                    # add data to the buffer -
                    buffer *= "$(metabolite_string)"
                end

                # do we have more?
                if (counter < number_of_species)
                    buffer *= "+"
                else
                    buffer *= ","
                end

                counter = counter + 1
            end
        end

        # find the products -
        idx_products = findall(stoichiometric_matrix[:,reaction_index].>0.0)
        if (isempty(idx_products) == true)
            buffer *= "[],"
        else

            # how many species do we have?
            number_of_species = length(idx_products)
            counter = 1
            for index in idx_products

                # get the metabolite -
                metabolite_string = list_of_species[index]
                stcoeff = stoichiometric_matrix[index,reaction_index]

                if (stcoeff != 1.0)
                    # add data to the buffer -
                    buffer *= "$(stcoeff)*$(metabolite_string)"
                else
                    # add data to the buffer -
                    buffer *= "$(metabolite_string)"
                end

                # do we have more?
                if (counter < number_of_species)
                    buffer *= "+"
                else
                    buffer *= ","
                end

                counter = counter + 1
            end
        end

        # is this reaction reversible?
        rev_value = list_of_reversible_reactions[reaction_index]
        if (rev_value == 1.0)
            buffer *= "-inf,inf"
        else
            buffer *= "0,inf"
        end

        # add buffer to list of strings -
        push!(reaction_string_buffer,buffer)
    end

    # return reaction_string_buffer
    return reaction_string_buffer
end


text:
### RandomReal

@mkapprule RandomReal

@doap RandomReal() = return rand()

### Random

@sjdoc Random """
    Random(Integer)

return 0 or 1 with equal probability.

    Random(Integer, [lo,hi])

return a random integer between `lo` and `hi`, inclusive.

    Random(Real)

return a random floating-point number between 0 and 1.

    Random(Complex)

return a random complex floating-point number with parts between 0 and 1.
"""

@mkapprule Random

@doap Random() = rand()

@doap function Random(sym::SJSym)
    if sym == :Integer
        rand(0:1)
    elseif sym == :Real
        rand()
    elseif sym == :Complex
        complex(rand(),rand())
    else
        mx
    end
end

@doap function Random(sym::SJSym, lst::ListT)
    if sym == :Integer
        rand(UnitRange(margs(lst)...))
    end
end

### RandomInteger

@sjdoc RandomInteger """
    Random([nmin,nmax])

return a random integer between `nmin` and `nmax`.

    RandomInteger(nmax)

return a random integer between 0 and `nmax`.

    RandomInteger()

return 0 or 1 with equal probablility.

    RandomInteger(range, n)

return an array of `n` random integers
"""

@mkapprule RandomInteger

@doap RandomInteger() = rand(0:1)
@doap RandomInteger(n::Integer) = rand(0:n)

@doap RandomInteger(x::Mxpr{:List}) = _randintrange(margs(x)...)

_randintrange(n::Integer,m::Integer) = rand(n:m)

@doap function RandomInteger(x::Mxpr{:List},n::Integer)
    lims = x[1]:x[2]
    nargs = newargs(n)
    for i=1:n  nargs[i] = rand(lims) end
    mxpr(:List,nargs)
end