text:
### SIMULATION STEP ###

# Agent
function step!(model::Simulation, t::EventTime, agent::Commuter, orig_map_nodes_num::Int)
	if agent.pos > orig_map_nodes_num
		# An agent can interact maximum with 10 agents on TTC station.
		max_inter_num = 10
	else
		# An agent can interact maximum with 3 agents on a street.
		max_inter_num = 3
	end
	lala = 0
    # infections
    if agent.infected == false # if healthy -> look at infected agents already in the node
        inf_no = sum([model.agents[ag].infected for ag in model.nodes_agents[agent.pos]])
		 # all infections have their own separate realizations of probabilities.
		if inf_no >= 1
			p = 1 - (1-model.p0)^min(inf_no,max_inter_num)
			agent.infected = (p > rand())
			if agent.infected == true
				agent.infected_in = :street
				model.infected_agents_count_current += 1
				model.infected_agents_street_current += 1
			end
		end
    else # for agents that are infected -> possibly infect healthy agents standing at the intersection
		num_of_interactions = 0 # an agent can interact maximum with 3 agents
        for i in model.nodes_agents[agent.pos]
            if !model.agents[i].infected
                p = model.p0  # just from that one agent!
                model.agents[i].infected = (p > rand())
                if model.agents[i].infected == true
                    model.agents[i].infected_in = :street
					model.infected_agents_count_current += 1
					model.infected_agents_street_current += 1
                end
            end

			num_of_interactions += 1
			if num_of_interactions >= max_inter_num
				break
			end
        end
    end
    # after the infections from previous position are handled, we can move the agent
    move_agent!(agent, model, orig_map_nodes_num)
    # println("S1! ","AGENT: ", agent.id," POS: ", agent.pos, " DIRECTION: ", agent.direction," TIME: ", t[1], " HOME:WORK: ", agent.home_loc,":",agent.work_loc)
    # If the current position is not in agent`s path it means that this node is the last one and we need to generate a new path.
    if agent.pos ∉ keys(agent.path)
        if agent.direction == :home_to_work
            working_day = 1 #rand(Normal(8,0.5),1)[1]*60*60
            queue_key = tuple(t[1] + working_day,rand())
            enqueue!(model.events, agent, queue_key)
            # println("I am at work. I will be here for ", working_day/60/60, " hours.")
        else
            time_to_rest_at_home = 1 # rand(Normal(16,0.5),1)[1] * 60 * 60
            queue_key = tuple(t[1] + time_to_rest_at_home, rand())
            enqueue!(model.events, agent, queue_key)
            # println("I am at home. I will take ", time_to_rest_at_home/60/60, " hours of rest.")
        end
    else
        # If an agent is not inside a TTC car we add next movement for him/her.
        if (agent.pos <= orig_map_nodes_num)
            queue_key = tuple(t[1] + agent.path_distances[agent.pos], rand())
            enqueue!(model.events, agent, queue_key)

        # If an agent is inside a TTC infrastucture, but the next node is outside the infrastructure we also add next movement.
        # For example, underpass.
        elseif (agent.pos > orig_map_nodes_num) && (agent.path[agent.pos] <= orig_map_nodes_num)
            queue_key = tuple(t[1] + agent.path_distances[agent.pos],rand())
            enqueue!(model.events, agent, queue_key)

        end
        # If an agent is inside a TTC car we don`t care because the TTC car will manage the queue.
    end
end

# TTC car
function step!(model::Simulation, t::EventTime, wagon::Wagon, orig_map_nodes_num::Int)
    # If it is the very beginning we just take passengers. TTC car will move only during the second step.
    if t[1] > 0.0
        move_agent!(wagon, model, t)
    end

    # If the current node is the last one on this route we will remove this TTC car.
    if wagon.pos == collect(values(wagon.path))[end]
        leave_passengers!(model,t,wagon)
        update_passengers_pos!(model,wagon, orig_map_nodes_num)
        model.max_passengers_per_TTC_car[(wagon.id, wagon.line, t[1])] = wagon.max_passengers_ever
        #This wagon have arrived to the last stop on the route, so it will die.
        filter!(e->e≠wagon,model.wagons)
    else
        # A wagon continues its path to the last node on the route.
        leave_passengers!(model,t,wagon)
        take_passengers!(model,t,wagon)
        update_passengers_pos!(model,wagon,orig_map_nodes_num)
        count_passengers_and_update_max_passengers_ever!(model,wagon)
        queue_key = tuple(t[1] + wagon.path_distances[wagon.pos],rand())
        enqueue!(model.events, wagon, queue_key)
    end

    # if it is time to run the next TTC car on this route we will run one.
    if (wagon.next_wagon_requested == false) && (t[1] >= (wagon.time_start + wagon.run_every))
       add_new_wagon!(model, wagon, orig_map_nodes_num)
    end

end

# Artificial agent
function step!(model::Simulation, t::EventTime, agent::Artificial_Commuter_Stats, orig_map_nodes_num::Int)
	model.infected_agents_count[floor(Int,t[1])+1:end] .= model.infected_agents_count_current
	model.infected_agents_wagon[floor(Int,t[1])+1:end] .= model.infected_agents_wagon_current
	model.infected_agents_street[floor(Int,t[1])+1:end] .= model.infected_agents_street_current
	queue_key = tuple(t[1] + agent.enter_every_sec, rand())
	enqueue!(model.events, agent, queue_key)
end

### Movement of an agent ###
function move_agent!(agent::Commuter, model::Simulation, orig_map_nodes_num::Int)
    # If an agent is inside a TTC the pos for agent will be the same as for the TTC.
        if agent.wagon_id != -999
            pos = filter(wag->wag.id==agent.wagon_id,model.wagons)[1].pos
        else
            #If an agent reached the last node a new travel destination will be selected.
            if !(agent.pos in keys(agent.path))
                if agent.direction == :home_to_work
                    ind = sample(1:length(keys(agent.paths_work_home)))
                    agent.path = agent.paths_work_home[ind]
                    agent.path_distances = agent.distances_work_home[ind]
                    agent.direction = :work_to_home
                else
                    ind = sample(1:length(keys(agent.paths_home_work)))
                    agent.path = agent.paths_home_work[ind]
                    agent.path_distances = agent.distances_home_work[ind]
                    agent.direction = :home_to_work
                end
                agent.no_of_trips += 1
                #println("agent ", agent.id, " walking ", agent.direction, ", route: ", agent.path)
                pos = agent.pos

            else
                # If an agent will use TTC and not only an underpass. We should not assign a new node for him, because
                # It will be assigned thanks to a TTC car.
                if (agent.pos > orig_map_nodes_num) && (agent.path[agent.pos] > orig_map_nodes_num)
                    # If agent changes line (transfer), so we need to change his pos
                    if find_route_id_by_stop_id(agent.pos - orig_map_nodes_num) != find_route_id_by_stop_id(agent.path[agent.pos] - orig_map_nodes_num)
                        node2 = agent.path[agent.pos]
                    else
                        node2 = agent.pos
                    end
                else
                    node2 = agent.path[agent.pos]
                end
                pos = deepcopy(node2)
        end
    end

    # Remove an agent from previous location and add to a new one.
    delete!(model.nodes_agents[agent.pos], agent.id)
    push!(model.nodes_agents[pos], agent.id)
    agent.pos = pos
end

### TTC CARS DYNAMIC ###

# Movement of a TTC car
function move_agent!(wagon::Wagon, model::Simulation, t::EventTime)
    # A wagon moves from one node to another on the route and disappear at the end.
	pos = wagon.path[wagon.pos]
    # Delete this wagon from the previous position and add to a new one.
    delete!(model.nodes_wagons[wagon.pos], wagon.id)
    push!(model.nodes_wagons[pos], wagon.id)
    wagon.pos = pos
end

# Request new wagons
function add_new_wagon!(model::Simulation, wagon::Wagon, orig_map_nodes_num)
    # New wagon definition

    new_wagon = Wagon(id = model.recent_wagon + 1
        ,pos = collect(keys(wagon.path))[1]
        ,line = wagon.line
        ,type = wagon.type
        ,path = wagon.path
        ,path_distances = wagon.path_distances
        ,time_start = wagon.time_start + wagon.run_every
		,run_every = wagon.run_every
        ,max_load = wagon.max_load)
    # Adding to the simulation
    push!(model.wagons, new_wagon)
    model.recent_wagon = new_wagon.id

    # Our old wagon requsted a new one.
    wagon.next_wagon_requested = true

    # New wagon gather all people which are on the first station
    time_new_wagon_starts = tuple(new_wagon.time_start,rand())
    leave_passengers!(model,time_new_wagon_starts,new_wagon)
    take_passengers!(model,time_new_wagon_starts,new_wagon)
    update_passengers_pos!(model,new_wagon, orig_map_nodes_num)
    # We create a new event for the wagon for the second station
    queue_key = tuple(time_new_wagon_starts[1]+new_wagon.path_distances[new_wagon.pos],rand())
    enqueue!(model.events, new_wagon, queue_key)
end

# Take passengers waiting on the station
function take_passengers!(model::Simulation, t::EventTime, wagon::Wagon)
    # A wagon takes passengers only when it has free places.
    # A wagon checks all passengers who are on the station.
    # It take only those passengers for whom the next node on the route is the same as for this wagon.
    for pass in model.nodes_agents[wagon.pos]
        if length(wagon.passengers) < ceil(wagon.max_load*model.max_load_gov_restr)
            ag = filter(ag->ag.id==pass,model.agents)[1]
            if (wagon.path[wagon.pos] == ag.path[wagon.pos]) && (ag.wagon_id == -999)
				# Tell anyone inside that this agent is infected
				if ag.infected == true
					for copass in wagon.passengers
						if copass.infected == false
							push!(copass.met_infected_in_wagon,ag.id)
						end
					end
				end
                # Update of co-riders
                push!(wagon.passengers,ag)
                # For each passenger we update in which wagon he/she is.
                ag.wagon_id = wagon.id
            end
        end
    end
end

# Leave passengers who need to leave a TTC car on specific point in time.
function leave_passengers!(model::Simulation, t::EventTime, wagon::Wagon)
    # If it is the last station everyone will be asked to leave a TTC car.
    if wagon.pos == collect(values(wagon.path))[end]
        current_passengers = copy(wagon.passengers) # copy because if there is no copy we will iterate through and modify the same array, so the loop will stop after the first change.
        # We remove all passengers and schedule next events for them.
        for pass in current_passengers
			if !pass.infected
				inf_no = length(pass.met_infected_in_wagon)
				if inf_no >= 1
					p0_TTC = model.p0 * 19
					# An agent can interact maximum with 10 agents in the TTC.
					p = 1 - (1-p0_TTC)^min(inf_no,10)
					pass.infected = (min(1, p) > rand())
					if pass.infected == true
						pass.infected_in = :wagon
						model.infected_agents_count_current += 1
						model.infected_agents_wagon_current += 1
					end
					pass.met_infected_in_wagon = Set{Int}()
				end
			end
			# Delete passengers who left a TTC car and add them new events.
            delete!(model.nodes_agents[pass.pos], pass.id)
            filter!(e->e!=pass,wagon.passengers)
            pass.wagon_id = -999

            pass.pos = deepcopy(wagon.pos)

            push!(model.nodes_agents[pass.pos], pass.id)

            queue_key = tuple(t[1] + pass.path_distances[pass.pos],rand())
            enqueue!(model.events, pass, queue_key)
        end
    else
        # If it is not the last station we check for each agent if this station is one to leave a TTC car.
        # If a next station for TTC is not included in an agent path it means that the current station is the last one.
        # If it is so we schedule a new event.
        current_passengers = copy(wagon.passengers) # copy because if there is no copy we will iterate through and modify the same array, so the loop will stop after the first change.
        for pass in current_passengers
           if wagon.path[wagon.pos] ∉ keys(pass.path)
                # Caculation of infection p and infection assgiment right after exiting a TTC car.
                if !pass.infected
					inf_no = length(pass.met_infected_in_wagon)
                    if inf_no >= 1
						p0_TTC = model.p0 * 19
						# An agent can interact maximum with 10 agents in the TTC.
						p = 1 - (1-p0_TTC)^min(inf_no,10)
                        pass.infected = (min(1, p) > rand())
                        if pass.infected == true
                            pass.infected_in = :wagon
							model.infected_agents_count_current += 1
							model.infected_agents_wagon_current += 1
                        end
						pass.met_infected_in_wagon = Set{Int}()
                    end
                end

                # Delete passengers who left a TTC car and add them new events.
				delete!(model.nodes_agents[pass.pos], pass.id)
                filter!(e->e!=pass,wagon.passengers)
                pass.wagon_id = -999

                pass.pos = deepcopy(wagon.pos)

                push!(model.nodes_agents[pass.pos], pass.id)

                queue_key = tuple(t[1] + pass.path_distances[pass.pos], rand())
                enqueue!(model.events, pass, queue_key)
            end
        end

    end
end

# For agents who are inside a TTC car. The TTC car updates their positions on each station.
function update_passengers_pos!(model::Simulation, wagon::Wagon, orig_map_nodes_num)
    # We update positions of all passengers according to the TTC car they are inside.
    current_passengers = copy(wagon.passengers) # copy because if there is no copy we will iterate through and modify the same array, so the loop will stop after the first change.
    move_agent!.(current_passengers, Ref(model), orig_map_nodes_num)
end

# Count passengers inside a TTC car on each station. Thanks to this we can check if there was a max_load issue.
function count_passengers_and_update_max_passengers_ever!(model::Simulation, wagon::Wagon)
    wagon.max_passengers_ever = max(wagon.max_passengers_ever,length(wagon.passengers))
end

println("Simulation dynamic functions were read")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
julia> if 1
         println("true")
       end
ERROR: type: non-boolean (Int64) used in boolean context

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CrossEntropyMethod
using POMDPs
using POMDPPlayback
using POMDPModels
using POMDPPolicies
using Distributions
using Random
using Test

mdp = SimpleGridWorld(tprob = 1)

backup_policy = RandomPolicy(mdp)
Distributions.logpdf(mdp::SimpleGridWorld, h) = length(h)*log(0.25)
Distributions.logpdf(p::RandomPolicy, h) = length(h)*log(0.25)

d = Dict{Symbol, Tuple{Sampleable, Int64}}(:a => (Categorical(4), 30))
loss = mdp_loss(mdp, backup_policy)
d_opt = cross_entropy_method(loss, d, max_iter = 10, N=1000)
@test argmax(d_opt[:a][1].p) == 4

mdp = SimpleGridWorld(tprob = 1, rewards = Dict(GWPos(1,1)=>0., GWPos(10,10)=>1.))
loss = mdp_loss(mdp, backup_policy)
weight = mdp_weight(mdp)
h = add_categorical_entropy([1., 1., 1., 1.])
d_opt = cross_entropy_method(loss, d, max_iter = 10, N=100, weight_fn = weight, add_entropy = h)

@test argmax(d_opt[:a][1].p) == 4 || argmax(d_opt[:a][1].p) == 1
d_opt[:a][1].p


text:
# The default purifier tags for LPDOs
const default_purifier_tags = ts"Purifier"

# Locally purified density operator
# L = prime(X, !purifier_tags(X)) * dag(X)
struct LPDO{XT<:Union{MPS,MPO}}
  X::XT
  purifier_tags::TagSet
  function LPDO(M::Union{MPO,MPS}, purifier_tags)
    sites = siteinds(all, M)
    if any(s -> has_indpairs(s, 0 => 1), sites)
      error(
        "In `LPDO(X::MPO, purifier_tags)`, `MPO` `X` must not have pairs of primed and unprimed site indices, since the `LPDO` is interpreted as `prime(X; tags = !purifier_tags) * dag(X)`.",
      )
    end
    return new{typeof(M)}(M, TagSet(purifier_tags))
  end
end

# MPS defaults to having no purifier tag
LPDO(M::MPS) = LPDO(M, random_tags())

length(L::LPDO) = length(L.X)

copy(L::LPDO) = LPDO(copy(L.X), L.purifier_tags)

# TODO: define this (not defined for MPS/MPO yet)
#Base.lastindex(L::LPDO) = lastindex(L.X)

function getindex(L::LPDO, args...)
  return error(
    "`getindex(L::LPDO, args...)` is purposefully not implemented yet. For the LPDO `L = X X†`, you can get the jth tensor `X[j]` with `L.X[j]`.",
  )
end

function setindex!(L::LPDO, args...)
  return error(
    "`setindex!(L::LPDO, args...)` is purposefully not implemented yet. For the LPDO `L = X X†`, you can set the jth tensor `X[j]` with `L.X[j] = A`.",
  )
end

purifier_tags(L::LPDO) = L.purifier_tags

ket(L::LPDO) = prime(L.X, !purifier_tags(L))
ket(L::LPDO, j::Int) = prime(L.X[j], !purifier_tags(L))
bra(L::LPDO) = dag(L.X)
bra(L::LPDO, j::Int) = dag(L.X[j])

function Base.iterate(L::LPDO, state=1)
  if state > 2 * length(L)
    return nothing
  elseif isodd(state)
    T = bra(L, (state + 1) ÷ 2)
  else
    T = ket(L, state ÷ 2)
  end
  return T, state + 1
end

ket(L::LPDO{MPS}) = prime(L.X)
ket(L::LPDO{MPS}, j::Int) = prime(L.X[j])
bra(L::LPDO{MPS}) = dag(L.X)
bra(L::LPDO{MPS}, j::Int) = dag(L.X[j])

tr(L::LPDO) = inner(L.X, L.X)

logtr(L::LPDO) = loginner(L.X, L.X)

norm(L::LPDO{MPS}) = norm(L.X)

norm(L::LPDO{MPO}) = norm(MPO(L))

# TODO: define siteinds, firstsiteind, allsiteind, etc.

"""
    normalize!(ψ::MPS; localnorms! = [])
    normalize!(M::MPO; localnorms! = [])
    normalize!(L::LPDO; sqrt_localnorms! = [])

Normalize the MPS/MPO/LPDO and returns the log of the norm and a vector of the local norms of each site.

An MPS `|ψ⟩` is normalized by `√⟨ψ|ψ⟩`, and the resulting MPS will have the property `norm(ψ) ≈ 1`.

An MPO `M` is normalized by `tr(M)`, and the resulting MPO will have the property `tr(M) ≈ 1`.

An LPDO `L = X X†` is normalized by `tr(L) = tr(X X†)`, so each `X` is normalized by `√tr(L) = √tr(X X†)`. The resulting LPDO will have the property `tr(L) ≈ 1`.

Passing a vector `v` as the keyword arguments `localnorms!` (`sqrt_localnorms!`) will fill the vector with the (square root) of the normalization factor per site. For an MPS `ψ`, `prod(v) ≈ norm(ψ)`. For an MPO `M`, `prod(v) ≈ tr(M)`. For an LPDO `L`, `prod(v)^2 ≈ tr(L)`.
"""
function normalize!(M::MPO; plev=0 => 1, tags=ts"" => ts"", (localnorms!)=[])
  N = length(M)
  resize!(localnorms!, N)
  blob = tr(M[1]; plev=plev, tags=tags)
  localZ = norm(blob)
  blob /= localZ
  M[1] /= localZ
  localnorms![1] = localZ
  for j in 2:(N - 1)
    blob *= M[j]
    blob = tr(blob; plev=plev, tags=tags)
    localZ = norm(blob)
    blob /= localZ
    M[j] /= localZ
    localnorms![j] = localZ
  end
  blob *= M[N]
  localZ = tr(blob; plev=plev, tags=tags)
  M[N] /= localZ
  localnorms![N] = localZ
  return M
end

function normalize!(L::LPDO; (sqrt_localnorms!)=[], localnorm=1.0)
  N = length(L)
  resize!(sqrt_localnorms!, N)
  # TODO: replace with:
  #blob = noprime(ket(L, 1) * siteind(L, 1)) * bra(L, 1)
  blob = noprime(ket(L, 1), "Site") * bra(L, 1)
  localZ = norm(blob)
  blob /= localZ
  L.X[1] /= sqrt(localZ / localnorm)
  sqrt_localnorms![1] = sqrt(localZ / localnorm)
  for j in 2:(length(L) - 1)
    # TODO: replace with:
    # noprime(ket(L, j), siteind(L, j))
    blob = blob * noprime(ket(L, j), "Site")
    blob = blob * bra(L, j)
    localZ = norm(blob)
    blob /= localZ
    L.X[j] /= sqrt(localZ / localnorm)
    sqrt_localnorms![j] = sqrt(localZ / localnorm)
  end
  # TODO: replace with:
  # noprime(ket(L, N), siteind(L, N))
  blob = blob * noprime(ket(L, N), "Site")
  blob = blob * bra(L, N)
  localZ = norm(blob)
  L.X[N] /= sqrt(localZ / localnorm)
  sqrt_localnorms![N] = sqrt(localZ / localnorm)

  return L
end

function normalize!(ψ::MPS; (localnorms!)=[])
  normalize!(LPDO(ψ); (sqrt_localnorms!)=localnorms!)
  return ψ
end

"""
    MPO(L::LPDO)

Contract the purifier indices to get the MPO
`ρ = L.X L.X†`. This contraction is performed exactly,
in the future we will support approximate contraction.
"""
function MPO(L::LPDO{MPO}; kwargs...)
  X = L.X
  X′ = prime(X; tags=!purifier_tags(L))
  X
  return *(X′, dag(X); kwargs...)
end

# TODO: implement in terms of `outer(L.X', L.X; kwargs...)
# It could also maybe just call `MPO(::LPDO{MPO})` if it
# is generic enough.
ITensors.MPO(L::LPDO{MPS}; kwargs...) = MPO(L.X; kwargs...)

"""
    tr(L::LPDO, tag::String)

Trace `"Input"` or `"Output"` qubits.
"""
function tr(L::LPDO, tag::String)
  N = length(L)

  Φ = ITensor[]

  tmp = noprime(ket(L, 1); tags=tag) * bra(L, 1)
  Cdn = combiner(commonind(tmp, L.X[2]), commonind(tmp, L.X[2]'))
  push!(Φ, tmp * Cdn)

  for j in 2:(N - 1)
    tmp = noprime(ket(L, j); tags=tag) * bra(L, j)
    Cup = Cdn
    Cdn = combiner(commonind(tmp, L.X[j + 1]), commonind(tmp, L.X[j + 1]'))
    push!(Φ, tmp * Cup * Cdn)
  end
  tmp = noprime(ket(L, N); tags=tag) * bra(L, N)
  Cup = Cdn
  push!(Φ, tmp * Cup)
  return MPO(Φ)
end

function HDF5.write(parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, L::LPDO)
  g = create_group(parent, name)
  attributes(g)["type"] = String(Symbol(typeof(L)))
  return write(parent, "state/X", L.X)
end

function HDF5.read(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{LPDO{XT}}
) where {XT}
  g = open_group(parent, name)
  X = read(g, "X", XT)
  return LPDO(X, ts"Purifier")
end


text:
# Code originally by Kristoffer Carlsson (2016-ish)
# Code modified by Jesse Chan (2020)
struct SparseMatrixBSC{Tv, Ti <: Integer}
    R::Int                 # Block size in rows
    C::Int                 # Block size in columns
    m::Int                 # Number of rows
    n::Int                 # Number of columns
    colptr::Vector{Ti}     # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::Vector{Ti}     # Row values of blocks
    nzval::Array{Tv, 3}    # Nonzero values, one "matrix" per block, nzval[i, j, block]

    function SparseMatrixBSC(R::Integer, C::Integer,
                             m::Integer, n::Integer,
                             colptr::Vector{Ti}, rowval::Vector{Ti},  nzval::Array{Tv, 3}) where {Ti, Tv}
      m < 0 && throw(ArgumentError("number of rows (m) must be ≥ 0, got $m"))
      n < 0 && throw(ArgumentError("number of columns (n) must be ≥ 0, got $n"))

      R < 0 && throw(ArgumentError("block size y must be ≥ 0, got $x_block_size"))
      C < 0 && throw(ArgumentError("block size x must be ≥ 0, got $y_block_size"))

      m % R != 0 && throw(ArgumentError("row block size: $(R) must evenly divide number of rows: $m"))
      n % C != 0 && throw(ArgumentError("column block size: $(C) must evenly divide number of rows: $n"))
      new{Tv, Ti}(Int(R), Int(C), Int(m), Int(n), colptr, rowval, nzval)
  end
end

function SparseMatrixBSC(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Array{Tv, 3}) where {Tv, Ti<:Integer}
    R, C = size(nzval, 1), size(nzval, 2)
    SparseMatrixBSC(R, C, m, n, colptr, rowval, nzval)
end

###########################
# Conversions CSC <-> BSC #
###########################

nblocks(A::SparseMatrixBSC) = (length(A.colptr) - 1, A.n ÷ A.C)
nblocks(A::SparseMatrixBSC, i::Int) = nblocks(A)[i]
blocksize(A::SparseMatrixBSC) = A.R, A.C
blocksize(A::SparseMatrixBSC, i::Int) = blocksize(A)[i]

# Column i is in colptr[i]:(colptr[i+1]-1)
nzblockrange(A::SparseMatrixBSC, col::Integer) =  Int(A.colptr[col]):Int(A.colptr[col + 1] - 1) # returns range for blocks in column "col"

Base.size(A::SparseMatrixBSC) = (A.m, A.n)
Base.size(A::SparseMatrixBSC, i) = size(A)[i]
Base.eltype(A::SparseMatrixBSC) = eltype(A.nzval)
SparseArrays.nnz(A::SparseMatrixBSC) = length(A.nzval)

function SparseArrays.SparseMatrixCSC(A::SparseMatrixBSC{Tv, Ti}) where {Tv, Ti <: Integer}
    if blocksize(A) == (1,1)
        return SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, vec(A.nzval))
    end
    rowval = zeros(Ti, nnz(A))
    colptr = zeros(Ti, size(A, 2) + 1)
    nzval = zeros(Tv, length(A.nzval))

    count_row = 1
    count_col = 2
    colptr[1] = 1
    @inbounds for col in 1:nblocks(A, 2)
        blockrange = nzblockrange(A, col)
        n_blocks_col = length(blockrange)
        nnz_values_col = n_blocks_col * blocksize(A, 1)
        for j_blk in 1:blocksize(A, 2)
            # The new colptr is the previous one plus the number of nonzero elements in this column.
            colptr[count_col] = colptr[count_col - 1] + nnz_values_col
            count_col += 1
            for block in blockrange
                i_offset = (A.rowval[block] - 1) * blocksize(A, 1)
                for i_blk in 1:blocksize(A, 1)
                    nzval[count_row] = A.nzval[i_blk, j_blk, block]
                    rowval[count_row] = i_blk + i_offset
                    count_row += 1
                end
            end
        end
    end
    return SparseMatrixCSC(A.m, A.n, colptr, rowval, nzval)
end

function SparseMatrixBSC(A::SparseMatrixCSC{Tv, Ti}, R::Integer, C::Integer) where {Tv, Ti}
    if (R, C) == (1,1)
        return SparseMatrixBSC(1, 1, A.m, A.n, A.colptr, A.rowval, reshape(A.nzval, length(A.nzval), 1, 1))
    end

    A.m % R != 0 && throw(ArgumentError("row block size: $(R) must evenly divide number of rows: $(A.m)"))
    A.n % C != 0 && throw(ArgumentError("column block size: $(C) must evenly divide number of rows: $(A.n)"))

    Anzval = A.nzval
    Arowval = A.rowval
    Acolptr = A.colptr

    # Upper bound of number of nonzero blocks is nnz(A).
    rowval = zeros(Ti, nnz(A))
    colptr = zeros(Ti, A.n ÷ C + 1)

    n_colblocks = div(A.n, C)

    colptr[1] = 1
    row_counter = 1
    rows = Int[]

    # Strategy to compute rowval and colptr:
    # For each column block accumulate all the rowvalues in the CSC matrix for that block.
    # Convert these to what rowblock they represent
    # Each unique rowblock should now be entered in order into rowval.
    # Colptr for this column block is incremented by the number of unique rowblock values.
    @inbounds for colblock in 1:n_colblocks
        j_offset = (colblock - 1) * C
        row_block_counter = 1

        # Count the number of non zero values for the columns in this column block
        nzvals_block = Acolptr[j_offset + C + 1] - Acolptr[j_offset + 1]
        if nzvals_block == 0 # No nz values in this column block, exit early
            colptr[colblock + 1] = row_counter
            continue
        end

        # Accumulate rowvals for this block
        resize!(rows, nzvals_block)
        for j_blk in 1:C
            col = j_offset + j_blk
            nz_range = Acolptr[col]:Acolptr[col + 1] - 1
            for r in nz_range
                rows[row_block_counter] = Arowval[r]
                row_block_counter += 1
            end
        end

        # Convert from row values -> block rows
        @simd for i in 1:length(rows)
            rows[i] = (rows[i] - 1) ÷ R + 1
        end

        # Pick out the unique block rows and put them into rowval

        #######################
        sort!(rows) # <- A bit of a bottle enck, takes about 30% of tot time.
        #######################
        rowval[row_counter] = rows[1] # We have at least one value in rows so this is ok
        row_counter += 1
        for i in 2:length(rows)
            if rows[i] > rows[i-1]
                rowval[row_counter] = rows[i]
                row_counter += 1
            end
        end
        colptr[colblock + 1] = row_counter
    end

    # We now know the true number of non zero blocks so we reshape rowval
    # and allocate the exact space we need for nzval
    deleteat!(rowval, row_counter:length(rowval))
    nzval = zeros(Tv, R, C, length(rowval))

    @inbounds for colblock in 1:n_colblocks
        j_offset = (colblock - 1) * C
        for j_blk in 1:C
            current_block = colptr[colblock]
            col = j_offset + j_blk
            for r in Acolptr[col]:Acolptr[col + 1] - 1
                row = Arowval[r]
                # Looking for the correct block for this column
                while row > rowval[current_block] * R
                    current_block += 1
                end
                i_blk = row - (rowval[current_block] - 1) * R
                nzval[i_blk, j_blk, current_block] = Anzval[r]
            end
        end
    end

    SparseMatrixBSC(A.m, A.n, colptr, rowval, nzval)
end

# for fast CSC conversion: convert to SparseMatrixCSC nzval ordering
function getCSCordering(A::SparseMatrixBSC)

    # loop through all active columns
    R,C = blocksize(A)
    ntotalblocks = size(A.nzval,3)
    flattened_indices = reshape(1:R*C*ntotalblocks,R,C,ntotalblocks) # UnitRange avoids allocation
    CSC_permuted_indices = zeros(Int,R,C*ntotalblocks) # [col1, col2, ..., col_blocksize*nblocks]

    sk = 1
    for i = 1:length(A.colptr)-1
        # find i,j indices of each block in a col
        blocks_in_column = nzblockrange(A,i)
        for col = 1:C
            for block_id in blocks_in_column
                CSC_permuted_indices[:,sk] .= flattened_indices[:,col,block_id]
                sk += 1
            end
        end
    end

    return CSC_permuted_indices[:]
end

# fast conversion from BSC to CSC matrix
function SparseMatrixCSC!(B::SparseMatrixCSC,A::SparseMatrixBSC{Tv, Ti}, CSCpermuted_indices) where {Tv, Ti <: Integer}
    B.nzval .= A.nzval[:][CSCpermuted_indices] # need to add permutation, but this should be super fast
end

############
# Base interfaces #
############

function Base.show(io::IO, A::SparseMatrixBSC;
                   header::Bool=true, repr=false)
    print(io, A.m, "×", A.n, " sparse block matrix with ", size(A.nzval,3), " ",
        eltype(A), " blocks of size ", blocksize(A, 1),"×",blocksize(A, 2))
end

# convert from block index to CSC index
function getCSCindex(A::SparseMatrixBSC{Tv,Ti}, index::Block{2,Ti}) where {Tv,Ti}
    row,col  = index.n
    colrange = nzblockrange(A,col)
    rows     = A.rowval[colrange]
    index    = findfirst(rows.==row)
    if isnothing(index)
        return nothing
    else
        return colrange[index]
    end
end

function Base.getindex(A::SparseMatrixBSC{Tv,Ti}, blockindex::Block{2,Ti}) where {Tv,Ti}
    CSC_id = getCSCindex(A,blockindex)
    if isnothing(CSC_id)
        return zeros(blocksize(A)...)
    else
        return A.nzval[:,:,CSC_id]
    end
end

function Base.setindex!(A::SparseMatrixBSC{Tv,Ti}, val::Array{Tv,2}, blockindex::Block{2,Ti}) where {Tv,Ti}
    CSC_id = getCSCindex(A,blockindex)
    if isnothing(CSC_id)
        error("block index not in original sparsity pattern. adding new sparse blocks not currently supported")
    end
    if size(val)!=blocksize(A)
        error("assigned block is size ", size(val), ", block size is ", blocksize(A))
    end
    A.nzval[:,:,CSC_id] .= val
end

# ##########
# # LinAlg #
# ##########
#
# # We do dynamic dispatch here so that the size of the blocks are known at compile time
# function LinearAlgebra.mul!(b::Vector{Tv}, A::SparseMatrixBSC{Tv}, x::Vector{Tv}) where {Tv}
#     if blocksize(A) == (1,1)
#         mul!(b, SparseMatrixCSC(A), x)
#     else
#         _mul!(b, A, x, Val(A.C), Val(A.R))
#     end
# end
#
# function _mul!(b::Vector{Tv}, A::SparseMatrixBSC{Tv}, x::Vector{Tv},
#                              ::Val{C}, ::Val{R}) where {Tv, C, R}
#     fill!(b, 0.0)
#     n_cb, n_rb = nblocks(A)
#     for J in 1:n_cb
#         j_offset = (J - 1)  * blocksize(A, 2)
#         for r in nzblockrange(A, J)
#             @inbounds i_offset = (A.rowval[r] - 1) * blocksize(A, 1)
#             matvec_kernel!(b, A, x, r, i_offset, j_offset, Val(C), Val(R))
#         end
#     end
#     return b
# end
#
# # TODO: Possibly SIMD.jl and do a bunch of simd magic coolness for small mat * vec
# @inline function matvec_kernel!(b::Vector{T}, A::SparseMatrixBSC{T}, x::Vector{T}, r,
#                                 i_offset, j_offset, ::Val{C}, ::Val{R}) where {C, R, T}
#     @inbounds for j in 1:C
#         for i in 1:R
#             b[i_offset + i] += A.nzval[i, j, r] * x[j_offset + j]
#         end
#     end
# end
#
# function Base.:*(A::SparseMatrixBSC{Tv}, x::Vector{Tv}) where {Tv}
#     mul!(similar(x, A.m), A, x)
# end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using Bandits
using Distributions
using Base.Test

T = 100
sb = staticbandit([Normal(0,1), Normal(100, 1)])
@test sb.optimal_expected_return == 100
@test regret(sb, 2) == 0

sb = staticbandit([Normal(0,1), Normal(100, 1)])
greedy_policy = GreedyPolicy()
epsilon_greedy_policy = EpsilonGreedyPolicy(0.1)
thompson_sampling = ThompsonSampling()
ucb1 = UCB1()
explore_then_exploit = ExploreThenExploit(T / 4)

agent = BasicAgent(greedy_policy, sb, 0, [0.0, 1.0])
stats = simulate(sb, agent, T)
@test stats.regret ≈ 0 atol=0.01
@test stats.arm_counts[2] ≈ T atol=1.10

sb = staticbandit([Normal(0, 1), Normal(100, 1)])
agent_epsilon = BasicAgent(epsilon_greedy_policy, sb, 0, [0.0, 1.0])
stats = simulate(sb, agent_epsilon, T)
@test stats.regret > 0
@test stats.arm_counts[2] < T

sb = staticbandit([Normal(0, 1), Normal(100, 1)])
agent_bad_prior = BasicAgent(greedy_policy, sb, 0, [1.0, 0.0])
stats = simulate(sb, agent_bad_prior, T)
@test stats.regret > 0
@test stats.arm_counts[2] < T

sb = staticbandit([Normal(0, 100), Normal(0.01, 100)])
normal_agent = NormalAgent(thompson_sampling, sb)
stats = simulate(sb, normal_agent, T)
@test stats.regret > 0

sb = staticbandit([Bernoulli(0.5), Bernoulli(0.6)])
beta_agent = BetaBernoulliAgent([0.6, 0.5], ucb1, sb)
stats = simulate(sb, beta_agent, T)
@test stats.regret > 0
@test stats.arm_counts[1] < T

sb = staticbandit([Bernoulli(0.5), Bernoulli(0.6)])
beta_agent = BetaBernoulliAgent([0.6, 0.5], thompson_sampling, sb)
stats = simulate(sb, beta_agent, T)
@test stats.regret > 0
@test stats.arm_counts[1] < T

sb = staticbandit([Bernoulli(0.5), Bernoulli(0.6)])
beta_agent = BetaBernoulliAgent([0.6, 0.5], explore_then_exploit, sb)
stats = simulate(sb, beta_agent, T)
@test stats.regret > 0
@test stats.arm_counts[1] < T


text:
# ======================= READING ====================

# Resolve `GlobalRef` instances to literals.
resolve(gr::GlobalRef) = getproperty(gr.mod, gr.name)
resolve(gr) = gr

function Pattern(ex, mod=@__MODULE__, resolve_fun=false)
    if !istree(ex)
        return PatLiteral(ex)
    end

    if ex isa Expr
        ex = cleanast(ex)
    end
    head = exprhead(ex)
    op = operation(ex)
    args = arguments(ex)

    if istree(op)
        error("Cannot yet match on composite expressions as function symbols")
    end

    if resolve_fun && op isa Symbol
        op = resolve(GlobalRef(mod, f))
    end

    patargs = map(i -> Pattern(i, mod, resolve_fun), args)

    if head == :(::)
        v = patargs[1]
        t = patargs[2]
        ty = Union{}
        if t isa PatVar
            ty = getfield(mod, t.name)
        elseif t isa PatLiteral{<:Type}
            ty = t.val
        end
        return PatTypeAssertion(v, ty)
    elseif head == :(...)
        v = patargs[1]
        if v isa PatVar
            return PatSplatVar(v)
        end
    end

    PatTerm(head, op, patargs)
end

function Pattern(x::Symbol, mod=@__MODULE__, resolve_fun=false)
    PatVar(x)
end

function Pattern(x::QuoteNode, mod=@__MODULE__, resolve_fun=false)
    if x.value isa Symbol
        PatLiteral(x.value)
    else
        PatLiteral(x)
    end
end

function Pattern(p::Pattern, mod=@__MODULE__, resolve_fun=false)
    p
end

macro pat(ex)
    Pattern(ex, __module__, false)
end

macro pat(ex, resolve_fun::Bool)
    Pattern(ex, __module__, resolve_fun)
end


text:
using DrWatson, Test, BSON, JLD2

T = 1000
N = 50 # spatial extent
Δt = 0.05 # resolution of integration
every = 10
seed = 1111

simulation = @ntuple T N Δt every seed

function f(simulation)
    @test typeof(simulation.T) <: Real
    a = rand(10); b = [rand(10) for _ in 1:10]
    return @strdict a b simulation
end

for ending ∈ ("bson", "jld2")
    @test !isfile(savename(simulation, "bson"))
    sim = produce_or_load(simulation, f; suffix = ending)
    @test isfile(savename(simulation, ending))
    @test sim["simulation"].T == T
    sim = produce_or_load(simulation, f; suffix = ending)
    @test sim["simulation"].T == T
    rm(savename(simulation, ending))
    @test !isfile(savename(simulation, ending))
end


text:
using JuliaExcel
using Dates
using Test

@testset "JuliaExcel.jl" begin
   @test JuliaExcel.encode_for_xl(1) == "&1"
   @test JuliaExcel.encode_for_xl(1.0) == "#1.0"
   @test JuliaExcel.encode_for_xl(Int8(1)) == "S1"
   @test JuliaExcel.encode_for_xl(Int16(1)) == "S1"
   @test JuliaExcel.encode_for_xl(Int32(1)) == "&1"
   @test JuliaExcel.encode_for_xl(Int64(1)) == "&1"
   @test JuliaExcel.encode_for_xl(true) == "T"
   @test JuliaExcel.encode_for_xl(false) == "F"
   @test JuliaExcel.encode_for_xl("foo") == "£foo"
   @test JuliaExcel.encode_for_xl('x') == "£x"
   @test JuliaExcel.encode_for_xl(:x) == "£:x"
   @test JuliaExcel.encode_for_xl(nothing) == "E"
   @test JuliaExcel.encode_for_xl(missing) == "E"
   @test JuliaExcel.encode_for_xl(Inf) =="!2036"
   @test JuliaExcel.encode_for_xl(-Inf) =="!2036"
   @test JuliaExcel.encode_for_xl(NaN) =="!2042"
   @test JuliaExcel.encode_for_xl(Date("2021-11-8")) == "D44508"
   @test JuliaExcel.encode_for_xl(DateTime("2021-11-8T12:00:00")) == "D44508.5"
   @test JuliaExcel.encode_for_xl(Int64) == "£Int64"
   @test JuliaExcel.encode_for_xl(v"1.2.3") == "£1.2.3"
   @test JuliaExcel.encode_for_xl((1,2)) == "*1,2;2,2,;&1&2"
   @test JuliaExcel.encode_for_xl([1,2,3]) == "*1,3;2,2,2,;&1&2&3"
   @test JuliaExcel.encode_for_xl(Any[1,2,3.0]) == "*1,3;2,2,4,;&1&2#3.0"
   @test JuliaExcel.encode_for_xl([1,true,"x"]) == "*1,3;2,1,2,;&1T£x"
   @test JuliaExcel.encode_for_xl([1,[2,3]]) == "*1,2;2,14,;&1*1,2;2,2,;&2&3"
   @test JuliaExcel.encode_for_xl(Dict("a"=>1,"b"=>2)) == "^2;2,2,2,2,;£b&2£a&1"
end


text:
using OrdinaryDiffEq

@testset "Gravitational Functional Test" begin
    G = 1
    @testset "The well-known figure Eight choreography" begin

        m1 = MassBody(SVector(-0.995492, 0.0, 0.0), SVector(-0.347902, -0.53393, 0.0), 1.0)
        m2 = MassBody(SVector(0.995492, 0.0, 0.0), SVector(-0.347902, -0.53393, 0.0), 1.0)
        m3 = MassBody(SVector(0.0, 0.0, 0.0), SVector(0.695804, 1.067860, 0.0), 1.0)
        t1 = 0.0
        t2 = 2pi
        tspan = (t1, t2);
        system = GravitationalSystem([m1, m2, m3], G)
        simulation = NBodySimulation(system, tspan)
        sim_result = run_simulation(simulation)
        solution_simo_3 = sim_result.solution;
        ε = 0.1
        for j = 1:3, i = 1:3
            @test solution_simo_3[1][i,j] ≈ solution_simo_3[end][i,j] atol = ε
        end

        @testset "Analyzing simulation result" begin
            r1 = get_position(sim_result, tspan[1], 1)
            r1_by_index = get_position(sim_result, tspan[1], 1)
            v1 = get_velocity(sim_result, tspan[1], 1)
            masses = get_masses(system)
            for i = 1:3
                @test r1[i] == (-0.995492, 0.0, 0.0)[i]
                @test r1_by_index[i] == (-0.995492, 0.0, 0.0)[i]
                @test v1[i] == (-0.347902, -0.53393, 0.0)[i]
                @test masses[i] == 1.0
            end

            ε = 0.001
            e_kin = 1.218
            @test e_kin ≈ kinetic_energy(sim_result, t1) atol = ε
        end

        @testset "Using convertion into SecondOrderODEProblem" begin
            sim_result = run_simulation(simulation, DPRKN6())
            solution_simo_3_2nd = sim_result.solution;
            ε = 0.001
            for i = 1:3, j = 1:3
                @test solution_simo_3_2nd[1][9 + 3(i - 1) + j] ≈ solution_simo_3_2nd[end][9 + 3(i - 1) + j] atol = ε
            end
        end

        @testset "Using symplectic integrators" begin
            sim_result = run_simulation(simulation, VelocityVerlet(), dt=pi / 130)
            solution_simo_3_2nd = sim_result.solution;
            ε = 0.001
            for i = 1:3, j = 1:3
                @test solution_simo_3_2nd[1][9 + 3(i - 1) + j] ≈ solution_simo_3_2nd[end][9 + 3(i - 1) + j] atol = ε
            end

            sim_result = run_simulation(simulation, Yoshida6(), dt=pi / 12)
            solution_simo_3_2nd = sim_result.solution;
            ε = 0.001
            for i = 1:3, j = 1:3
                @test solution_simo_3_2nd[1][9 + 3(i - 1) + j] ≈ solution_simo_3_2nd[end][9 + 3(i - 1) + j] atol = ε
            end
        end
    end

    @testset "5-body choreography" begin
        m1 = MassBody(SVector(1.657666, 0.0, 0.0), SVector(0.0, -0.593786, 0.0), 1.0)
        m2 = MassBody(SVector(0.439775, -0.169717, 0.0), SVector(1.822785, 0.128248, 0.0), 1.0)
        m3 = MassBody(SVector(-1.268608, -0.267651, 0.0), SVector(1.271564, 0.168645, 0.0), 1.0)
        m4 = MassBody(SVector(-1.268608, 0.267651, 0.0), SVector(-1.271564, 0.168645, 0.0), 1.0)
        m5 = MassBody(SVector(0.439775, 0.169717, 0.0), SVector(-1.822785, 0.128248, 0.0), 1.0)
        tspan = (0.0, 2pi);
        system = GravitationalSystem([m1, m2, m3, m4, m5], G)
        simulation = NBodySimulation(system, tspan)
        sim_result = run_simulation(simulation, Tsit5(), abstol=1e-10, reltol=1e-10)
        solution_simo_5 = sim_result.solution;

        ε = 0.01
        for j = 1:5, i = 1:3
            @test solution_simo_5[1][i,j] ≈ solution_simo_5[end][i,j] atol = ε
        end
    end

    @testset "Constructing electorstatic potential parameters entity" begin
        default_potential = GravitationalParameters()
        @test 6.67408e-11 == default_potential.G

        io = IOBuffer()
        potential1 = GravitationalParameters()
        potential2 = GravitationalParameters(35.67)
        @test sprint(io -> show(io, potential1)) == "Gravitational:\n\tG:6.67408e-11\n"
        @test sprint(io -> show(io, potential2)) == "Gravitational:\n\tG:35.67\n"
    end
end


text:
function Export(Format::String,Args::AbstractVector{Any})
    if (Format=="Blender")
        Bl_Export(Args[1],Args[2],Args[3])
    end
end

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
# =======================================================================
# 2KP

using vOptSpecific

# ==========================================
# instance

n  = 3
p1 = [ 3,  9 , 7]
p2 = [16, 15 , 6]
w  = [ 2,  3 , 4]
c  = 8

# ==========================================

myInstance = set2UKP( n , p1 , p2 , w , c )
algorithm = UKP_Jorge2010()
z1, z2, r, x = vSolve( myInstance , algorithm )

for i = 1:length(z1)
    println("(" , z1[i], " | ", z2[i], ") : ", x[i,:])
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
using ValkyrieRobot
using ValkyrieRobot.BipedControlUtil
using RigidBodyTreeInspector
using Base.Test

@testset "side" begin
    @test -left == right
    @test -right == left

    sides = [rand(Side) for i = 1 : 10000];
    @test isapprox(count(side -> side == left, sides) / length(sides), 0.5; atol = 0.05)

    @test flipsign_if_right(2., left) == 2.
    @test flipsign_if_right(2., right) == -2.
end

@testset "load geometries" begin
    val = Valkyrie()
    geometries = parse_urdf(ValkyrieRobot.urdfpath(), val.mechanism; package_path = [ValkyrieRobot.packagepath()]);
    @test length(geometries) == 81
end


text:
"""
    struct HeaderSearch <: Any
Hold a pointer to a `clang::HeaderSearch` object.
"""
struct HeaderSearch
    ptr::CXHeaderSearch
end

Base.unsafe_convert(::Type{CXHeaderSearch}, x::HeaderSearch) = x.ptr
Base.cconvert(::Type{CXHeaderSearch}, x::HeaderSearch) = x


text:
using BenchmarkTools
using LoopVectorization
using LinearAlgebra
using Random

#Original functions (smooth_line and smooth_seq using @simd)
@inline function smooth_line(nrm1,j,i1,sl,rl,ih2,denom)
    @fastmath @inbounds @simd for i=i1:2:nrm1
            # @show i,j,denom,sl[i,j],rl[i,j],sl[i,j-1],sl[i,j+1],sl[i-1,j],sl[i+1,j]
            sl[i,j]=denom*(rl[i,j]+ih2*(sl[i,j-1]+sl[i-1,j]+sl[i+1,j]+sl[i,j+1]))
    end
end

@inline function smooth_seq(nrows,ncols,sl,rl,ih2)
    denom=1/(4ih2)
    nrm1=nrows-1
    #even (i+j)%2==0 from odds
    @inbounds for j=2:2:ncols-1
        smooth_line(nrm1,j,2,sl,rl,ih2,denom)
        smooth_line(nrm1,j+1,3,sl,rl,ih2,denom)
    end
    # #odds (i+j)%2==1 from evens
    @inbounds for j=2:2:ncols-1
        smooth_line(nrm1,j,3,sl,rl,ih2,denom)
        smooth_line(nrm1,j+1,2,sl,rl,ih2,denom)
    end

end

#RBArray is a Type acting like a 2D array with a data layout suitable to
# Red Black Gauss-Seidel algorithms
#  contiguous red (iseven(i+j)) values are stored in [i2,1,j]
#  contiguous odd (isodd(i+j)) values are stores in [i2,2,j]
struct RBArray
    data::Array{Float64,3}
end
#Ctor : should use an inner Ctor to ensure that nx is even
#Note that >>> means unsigned >>. Hence x >>> 1 returns div(x,2) is x is a positive integer
RBArray(nx,ny) = RBArray(Array{Float64,3}(undef,nx>>>1,2,ny))
#Helper functions to acess the the RB arrays with 2 indexes rb[i,j]
@inline Base.getindex(rb::RBArray,i,j) = rb.data[(i+1)>>>1,1+isodd(i+j),j]
@inline Base.setindex!(rb::RBArray,value,i,j) = rb.data[(i+1)>>>1,1+isodd(i+j),j]=value
Base.size(rb::RBArray) = (size(rb.data)[1]*2,size(rb.data)[3])
Base.fill!(rb::RBArray,val) = fill!(rb.data,val)
Base.eltype(b::RBArray) = Float64

#RBArray from Array
function RBArray(a::Array{Float64,2})
    nx,ny=size(a)
    @assert iseven(nx)
    nxs2=nx>>>1
    rb=RBArray(nx,ny)
    for i in 1:nx,j in 1:ny
        rb[i,j]=a[i,j]
    end
    for j in 1:2:ny
        #j odd
        for i2 in 0:nxs2-1
            rb.data[i2+1,1,j]=a[2i2+1,j] #i=2i2+1 is odd  => i+j is even (rb.data[:,1,:])
            rb.data[i2+1,2,j]=a[2i2+2,j] #i=2i2+2 is even => i+j is odd (rb.data[:,2,:])
        end
        #j+1 even
        for i2 in 0:nxs2-1 #i=2i2 is even
            rb.data[i2+1,2,j+1]=a[2i2+1,j+1] #i=2i2+1 is odd  => i+j is even (rb.data[:,2,:])
            rb.data[i2+1,1,j+1]=a[2i2+2,j+1] #i=2i2+2 is even => i+j is odd (rb.data[:,1,:])
        end
    end
    rb
end

#Array from RBArray
function Array{Float64,2}(rb::RBArray)
    nxs2,two,ny=size(rb.data)
    nx=2nxs2
    a=Array{Float64,2}(undef,nx,ny)
    #The following blocks works but is probably slow
    # for i in 1:nx,j in 1:ny
    #     a[i,j]=rb[i,j]
    # end
    for j in 1:2:ny
        #j odd
        for i2 in 0:nxs2-1 #i=2i2 is even
            a[2i2+1,j]=rb.data[i2+1,1,j]
            a[2i2+2,j]=rb.data[i2+1,2,j]
        end
        #j+1 even
        for i2 in 0:nxs2-1 #i=2i2 is even
            a[2i2+1,j+1]=rb.data[i2+1,2,j+1]
            a[2i2+2,j+1]=rb.data[i2+1,1,j+1]
        end
    end
    a
end

@inline function smooth_line_rb(nrm1s2,j,di,rbout,rbin,sl,rl,ih2,denom)
    @avx for i2 in 1:nrm1s2
        sl[i2+di,rbout,j]=denom*(rl[i2+di,rbout,j]+ih2*(sl[i2+di,rbin,j-1]+sl[i2+di,rbin,j+1]+sl[i2,rbin,j]+sl[i2+1,rbin,j]))
    end
end

# #Version with broadcast (slow but suitable to test CuArrays)
@inline function smooth_line_rb_brcst(nrm1s2,j,di,rbout,rbin,sl,rl,ih2,denom)
    ib=1+di
    ie=nrm1s2+di
    @inbounds @fastmath @views @. sl[ib:ie,rbout,j]=denom*(rl[ib:ie,rbout,j]+ih2*(sl[ib:ie,rbin,j-1]+sl[ib:ie,rbin,j+1]+sl[1:nrm1s2,rbin,j]+sl[2:nrm1s2+1,rbin,j]))
end

#Specialized smooth_seq method for RedBlack arrays.
#Note that non specialized smooth_seq method should produce the same result(but slowly)
@inline function smooth_seq(nrows,ncols,slrb::RBArray,rlrb::RBArray,ih2)
    denom=1/(4ih2)
    nrm1,ncm1=nrows-1,ncols-1
    nrm1s2=nrm1 >>> 1
    sl,rl=slrb.data,rlrb.data

    #evens from odds
    @inbounds for j in 2:2:ncm1
        smooth_line_rb(nrm1s2,j  ,0,1,2,sl,rl,ih2,denom)
        smooth_line_rb(nrm1s2,j+1,1,1,2,sl,rl,ih2,denom)
    end
     #odds from evens
     @inbounds for j in 2:2:ncm1
        smooth_line_rb(nrm1s2,j  ,1,2,1,sl,rl,ih2,denom)
        smooth_line_rb(nrm1s2,j+1,0,2,1,sl,rl,ih2,denom)
    end
end

#function that call the benchmarks
function benchmark_smooth_gs_rb(n)
    @assert iseven(n)
    # initialize two 2D arrays
    sl=rand(n+2,n+2)
    rl=rand(n+2,n+2)
    # the following blocks check if smooth_seq and smooth_seq_avx
    # return the same result
    sl_ref=deepcopy(sl)
    slrb=RBArray(sl)
    rlrb=RBArray(rl)
    (nrows,ncols)=size(sl)
    ih2=rand()
    smooth_seq(nrows,ncols,sl,rl,ih2)
    smooth_seq(nrows,ncols,slrb,rlrb,ih2)
    res=norm(sl-Array{Float64,2}(slrb))

    # #The actual timings
    # bs=@btime smooth_seq($nrows,$ncols,$slrb,$rlrb,$ih2)
    tsimd=@belapsed smooth_seq($nrows,$ncols,$sl,$rl,$ih2)
    trb=@belapsed smooth_seq($nrows,$ncols,$slrb,$rlrb,$ih2)
    @show n,res,tsimd,trb,tsimd/trb
    nothing
end

function interpolate_correct!(uf,uc)
    (nrows,ncols)=size(uc)
    # @show size(uc)
    ncm1=ncols-1
    nrm1=nrows-1

    @inbounds @fastmath @simd for j=2:ncm1
        fj=2j-1
        for i=2:nrm1
            fi=2i-1
             v=uc[i,j]
             uf[fi  ,fj  ]+=v
             uf[fi-1,fj  ]+=v
             uf[fi  ,fj-1]+=v
             uf[fi-1,fj-1]+=v
        end
    end

end

@inline function scatter!(df,v,i,fj)
    df[i-1,1,fj-1]+=v
    df[i,1,fj]+=v
    df[i-1,2,fj]+=v
    df[i,2,fj-1]+=v
end

function interpolate_correct!(uf::RBArray,uc::RBArray)
    (nrows,ncols)=size(uc)
    # @show size(uc)
    ncm1=ncols-1
    nrm1=nrows-1

    df,dc=uf.data,uc.data

    ncm1s2=ncm1>>>1
    nrm1s2=nrm1>>>1

    @inbounds @fastmath @simd for j2=1:ncm1s2
    # @avx @fastmath for j2=1:ncm1s2
            for i2=1:nrm1s2
            # @inbounds @simd for j2=1:ncm1>>>1
            #     @fastmath @simd for i2=1:nrm1>>>1
            v1,v2,v3,v4=dc[i2,1,2j2],dc[i2+1,2,2j2],dc[i2,2,2j2+1],dc[i2+1,1,2j2+1]

            df[2i2-1,1,4j2-2]+=v1
            df[2i2  ,1,4j2-2]+=v2

            df[2i2  ,1,4j2-1]+=v1
            df[2i2+1,1,4j2-1]+=v2

            df[2i2-1,1,4j2  ]+=v3
            df[2i2  ,1,4j2  ]+=v4

            df[2i2  ,1,4j2+1]+=v3
            df[2i2+1,1,4j2+1]+=v4

            df[2i2-1,2,4j2-1]+=v1
            df[2i2,  2,4j2-2]+=v1

            df[2i2  ,2,4j2-1]+=v2
            df[2i2+1,2,4j2-2]+=v2

            df[2i2-1,2,4j2+1]+=v3
            df[2i2,  2,4j2  ]+=v3

            df[2i2  ,2,4j2+1]+=v4
            df[2i2+1,2,4j2  ]+=v4

        end
    end

    # @inbounds @simd for j2=1:ncm1>>>1
    #     @fastmath @simd for i2=1:nrm1>>>1
    #         v1,v2,v3,v4=dc[i2,1,2j2],dc[i2+1,2,2j2],dc[i2,2,2j2+1],dc[i2+1,1,2j2+1]

    #         df[2i2-1,1,4j2-2]+=v1
    #         df[2i2  ,1,4j2-1]+=v1
    #         df[2i2-1,2,4j2-1]+=v1
    #         df[2i2,  2,4j2-2]+=v1

    #         df[2i2  ,1,4j2-2]+=v2
    #         df[2i2+1,1,4j2-1]+=v2
    #         df[2i2  ,2,4j2-1]+=v2
    #         df[2i2+1,2,4j2-2]+=v2

    #         df[2i2-1,1,4j2  ]+=v3
    #         df[2i2  ,1,4j2+1]+=v3
    #         df[2i2-1,2,4j2+1]+=v3
    #         df[2i2,  2,4j2  ]+=v3

    #         df[2i2  ,1,4j2  ]+=v4
    #         df[2i2+1,1,4j2+1]+=v4
    #         df[2i2  ,2,4j2+1]+=v4
    #         df[2i2+1,2,4j2  ]+=v4

    #         # scatter!(df,v1,2i2,4j2-1)
    #         # scatter!(df,v2,2i2+1,4j2-1)
    #         # scatter!(df,v3,2i2,4j2+1)
    #         # scatter!(df,v4,2i2+1,4j2+1)
    #     end
    # end

    # @inbounds @simd for j2=1:ncm1>>>1
    #     @fastmath @simd for i2=1:nrm1>>>1
    #         v1,v2,v3,v4=dc[i2,1,2j2],dc[i2+1,2,2j2],dc[i2,2,2j2+1],dc[i2+1,1,2j2+1]

    #         scatter!(df,v1,2i2,4j2-1)
    #         scatter!(df,v2,2i2+1,4j2-1)
    #         scatter!(df,v3,2i2,4j2+1)
    #         scatter!(df,v4,2i2+1,4j2+1)
    #     end
    # end

end

function interpolate_correct!(uf::RBArray,uc::RBArray)
    (nrows,ncols)=size(uc)
    # @show size(uc)
    ncm1=ncols-1
    nrm1=nrows-1

    df,dc=uf.data,uc.data

    ncm1s2=ncm1>>>1
    nrm1s2=nrm1>>>1

    @inbounds @fastmath for jcoarse2=1:ncm1s2
        for djcoarse=0:1
            jcoarse=2jcoarse2+djcoarse

            @simd for icoarse2=1:ncm1s2

                for dicoarse=0:1
                    icoarse=2icoarse2+dicoarse
                    v=dc[(icoarse+1)>>>1,1+isodd(jcoarse+icoarse),jcoarse]

                    # ifine=2icoarse-2
                    jfine=2jcoarse-2

                    # uf[ifine  ,jfine  ]+=v
                    # uf[ifine+1,jfine  ]+=v
                    # uf[ifine  ,jfine+1]+=v
                    # uf[ifine+1,jfine+1]+=v

                    df[icoarse-1,1,jfine]+=v
                    df[icoarse  ,2,jfine]+=v
                    df[icoarse-1,2,jfine+1]+=v
                    df[icoarse  ,1,jfine+1]+=v
                end
            end
        end
    end
end

#function that call the benchmarks
function benchmark_interpolate_rb(nf)
    @assert iseven(nf)
    # initialize two 2D arrays (fine and coarse)
    nc=nf>>>1
    uf=rand(nf+2,nf+2)
    uc=rand(nc+2,nc+2)
    # the following blocks check if smooth_seq and smooth_seq_avx
    # return the same result
    ufrb=RBArray(uf)
    ucrb=RBArray(uc)
    interpolate_correct!(uf,uc)
    interpolate_correct!(ufrb,ucrb)
    res=norm(uf-Array{Float64,2}(ufrb))
    @show res
    # #The actual timings
    tsimd=@belapsed interpolate_correct!($uf,$uc)
    trb=@belapsed interpolate_correct!($ufrb,$ucrb)
    @show "interp",nf,res,tsimd,trb,tsimd/trb
    nothing
end

@inline function fres(x,y,i,j,ih2,nx,ny)
    # @show x[i,j]+ih2*(y[i+1,j]+y[i,j+1]+y[i,j-1]+y[i-1,j]-eltype(x)(4)*y[i,j])
    # @show x[i,j]+ih2*(y[i+1,j]+y[i,j+1]+y[i,j-1]+y[i-1,j]-eltype(x)(4)*y[i,j])
    # x[i,j]+ih2*(y[i+1,j]+y[i,j+1])
    @inbounds x[i,j]+ih2*(y[i+1,j]+y[i,j+1]+y[i,j-1]+y[i-1,j]-eltype(x)(4)*y[i,j])
end

#rlc = rl coarse
function restrict_residual!(rlc,slf,rlf,ih2)

    (nrows,ncols) = size(rlc)
    (nx,ny) = (nrows-1,ncols-1)
    fill!(rlc,0.)

    @inline finer(i,j)=fres(rlf,slf,i,j,ih2,nx,ny)
    oneOverFour=eltype(rlc)(1/4)

    @simd for j=2:ncols-1
        fj=2j-2
        for i=2:nrows-1
            fi=2i-2
            # rlc[i,j]=oneOverFour*(finer(fi,fj)+finer(fi+1,fj))
            # rlc[i,j]=oneOverFour*(finer(fi+1,fj+1))
             @inbounds rlc[i,j]=oneOverFour*(finer(fi,fj)+finer(fi+1,fj)+finer(fi,fj+1)+finer(fi+1,fj+1))
        end
    end
end

# #rlc = rl coarse
# function restrict_residual!(rlc::RBArray,slf::RBArray,rlf::RBArray,ih2)

#     (nrows,ncols) = size(rlc)
#     (nx,ny) = (nrows-1,ncols-1)
#     fill!(rlc,0.)

#     @inline finer(i,j)=fres(rlf,slf,i,j,ih2,nx,ny)
#     oneOverFour=eltype(rlc)(1/4)
#     four=eltype(rlc)(4)
#     x,y=rlf.data,slf.data

#     nrm1s2 = (nrows>>>1)-1
#     ncm1s2 = (ncols>>>1)-1

#     @inbounds for j2=1:ncm1s2
#         # j=2j2
#         # fj=2j-2
#         @fastmath @simd for i2=1:nrm1s2
#             rlc.data[i2  ,1,2j2] =
#             0.25*(
#                 x[2i2-1,1,4j2-2  ]+2ih2*(y[2i2-1,2,4j2-2  ]+y[2i2-1,2,4j2-2+1]+y[2i2-1,2,4j2-2-1]+y[2i2  ,2,4j2-2  ]-4.0*y[2i2-1,1,4j2-2  ])+
#                 x[2i2+0,2,4j2-2  ]+2ih2*(y[2i2+0,1,4j2-2  ]+y[2i2  ,1,4j2-2+1]+y[2i2  ,1,4j2-2-1]+y[2i2-1,1,4j2-2  ]-4.0*y[2i2  ,2,4j2-2  ])+
#                 x[2i2-1,2,4j2-2+1]+2ih2*(y[2i2+0,1,4j2-2+1]+y[2i2  ,1,4j2-2+1]+y[2i2-1,1,4j2-2+2]+y[2i2-1,1,4j2-2+1]-4.0*y[2i2-1,2,4j2-2+1])+
#                 x[2i2+0,1,4j2-2+1]+2ih2*(y[2i2+0,2,4j2-2+1]+y[2i2  ,2,4j2-2+2]+y[2i2  ,2,4j2-2  ]+y[2i2-1,2,4j2-2+1]-4.0*y[2i2,1  ,4j2-2+1]))

#             # rlc.data[i2  ,1,j] = restrict_cell(x,y,2i2,2j-2,ih2)
#             # rlc.data[i2+1,2,j] = restrict_cell(x,y,2i2+1,2j-2,ih2)
#         end
#         # j=2j2+1
#         # fj=2j-2
#         # @fastmath  @simd for i2=1:nrm1s2
#         #     i=2i2
#         #     fi=2i-2
#         #     rlc.data[i2,2,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#         #     i=2i2+1
#         #     fi=2i-2
#         #     rlc.data[i2+1,1,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#         #     # @inbounds rlc[i,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#         # end

@inline fres_ee(x,y,i,fj,ih2) = x[i,1,fj]+ih2*(y[i  ,2,fj]+y[i,2,fj+1]+y[i,2,fj-1]+y[i+1,2,fj]-4.0*y[i,1,fj])
@inline fres_oe(x,y,i,fj,ih2) = x[i,2,fj]+ih2*(y[i  ,1,fj]+y[i,1,fj+1]+y[i,1,fj-1]+y[i-1,1,fj]-4.0*y[i,2,fj])
@inline fres_eo(x,y,i,fj,ih2) = x[i,2,fj]+ih2*(y[i+1,1,fj]+y[i,1,fj+1]+y[i,1,fj-1]+y[i  ,1,fj]-4.0*y[i,2,fj])
@inline fres_oo(x,y,i,fj,ih2) = x[i,1,fj]+ih2*(y[i  ,2,fj]+y[i,2,fj+1]+y[i,2,fj-1]+y[i-1,2,fj]-4.0*y[i,1,fj])

# #rlc = rl coarse
# function restrict_residual!(rlc::RBArray,slf::RBArray,rlf::RBArray,ih2)

#     (nrows,ncols) = size(rlc)
#     (nx,ny) = (nrows-1,ncols-1)
#     fill!(rlc,0.)

#     @inline finer(i,j)=fres(rlf,slf,i,j,ih2,nx,ny)
#     oneOverFour=eltype(rlc)(1/4)
#     four=eltype(rlc)(4)
#     x,y=rlf.data,slf.data

#     @simd for j=2:ncols-1
#         fj=2j-2
#         for i=2:nrows-1
#             fi=2i-2
#             @inbounds rlc[i,j]=oneOverFour*(
#                 fres_ee(x,y,i-1,fj,ih2,nx,ny)+
#                 fres_oe(x,y,i,fj,ih2,nx,ny)+
#                 fres_eo(x,y,i-1,fj+1,ih2,nx,ny)+
#                 fres_oo(x,y,i,fj+1,ih2,nx,ny)
#                 # fres_e(x,y,i,fj+1,ih2,nx,ny)
#                 # rlf[fi,fj]+ih2*(slf[fi+1,fj]+slf[fi,fj+1]+slf[fi,fj-1]+slf[fi-1,fj]-four*slf[fi,fj])+
#                 # rlf[fi+1,fj]+ih2*(slf[fi+2,fj]+slf[fi+1,fj+1]+slf[fi+1,fj-1]+slf[fi  ,fj]-four*slf[fi+1,fj])+
#                 # rlf[fi,fj+1]+ih2*(slf[fi+1,fj+1]+slf[fi,fj+2]+slf[fi,fj  ]+slf[fi-1,fj+1]-four*slf[fi,fj+1])+
#                 # rlf[fi+1,fj+1]+ih2*(slf[fi+2,fj+1]+slf[fi+1,fj+2]+slf[fi+1,fj  ]+slf[fi  ,fj+1]-four*slf[fi+1,fj+1])
#             )
#         end
#     end
# end

# @inline function restrict_cell_ifj(x,y,i,fj,ih2)
#     0.25*(
#         x[i-1,1,fj  ]+ih2*(y[i-1,2,fj  ]+y[i-1,2,fj+1]+y[i-1,2,fj-1]+y[i  ,2,fj  ]-4.0*y[i-1,1,fj  ])+
#         x[i  ,2,fj  ]+ih2*(y[i  ,1,fj  ]+y[i  ,1,fj+1]+y[i  ,1,fj-1]+y[i-1,1,fj  ]-4.0*y[i  ,2,fj  ])+
#         x[i-1,2,fj+1]+ih2*(y[i  ,1,fj+1]+y[i  ,1,fj+1]+y[i-1,1,fj+2]+y[i-1,1,fj+1]-4.0*y[i-1,2,fj+1])+
#         x[i  ,1,fj+1]+ih2*(y[i  ,2,fj+1]+y[i  ,2,fj+2]+y[i  ,2,fj  ]+y[i-1,2,fj+1]-4.0*y[i,1  ,fj+1]))
#     # fres_ee(x,y,i-1,fj,ih2)+
#     # fres_oe(x,y,i,fj,ih2)+
#     # fres_eo(x,y,i-1,fj+1,ih2)+
#     # fres_oo(x,y,i,fj+1,ih2))
# end

# @inline function restrict_cell(x,y,i,fj,ih2)
#     0.25*(
#     fres_ee(x,y,i-1,fj,ih2)+
#     fres_oe(x,y,i,fj,ih2)+
#     fres_eo(x,y,i-1,fj+1,ih2)+
#     fres_oo(x,y,i,fj+1,ih2))
# end

# #rlc = rl coarse
# function restrict_residual!(rlc::RBArray,slf::RBArray,rlf::RBArray,ih2)

#     (nrows,ncols) = size(rlc)
#     (nx,ny) = (nrows-1,ncols-1)
#     fill!(rlc,0.)

#     @inline finer(i,j)=fres(rlf,slf,i,j,ih2,nx,ny)
#     oneOverFour=eltype(rlc)(1/4)
#     four=eltype(rlc)(4)
#     x,y=rlf.data,slf.data

#     nrm1s2 = (nrows>>>1)-1
#     ncm1s2 = (ncols>>>1)-1

#     @inbounds for j2=1:ncm1s2
#         # j=2j2
#         # fj=2j-2
#         @fastmath @simd for i2=1:nrm1s2
#             rlc.data[i2  ,1,2j2] =
#             0.25*(
#                 x[2i2-1,1,4j2-2  ]+2ih2*(y[2i2-1,2,4j2-2  ]+y[2i2-1,2,4j2-2+1]+y[2i2-1,2,4j2-2-1]+y[2i2  ,2,4j2-2  ]-4.0*y[2i2-1,1,4j2-2  ])+
#                 x[2i2  ,2,4j2-2  ]+2ih2*(y[2i2  ,1,4j2-2  ]+y[2i2  ,1,4j2-2+1]+y[2i2  ,1,4j2-2-1]+y[2i2-1,1,4j2-2  ]-4.0*y[2i2  ,2,4j2-2  ])+
#                 x[2i2-1,2,4j2-2+1]+2ih2*(y[2i2  ,1,4j2-2+1]+y[2i2  ,1,4j2-2+1]+y[2i2-1,1,4j2-2+2]+y[2i2-1,1,4j2-2+1]-4.0*y[2i2-1,2,4j2-2+1])+
#                 x[2i2  ,1,4j2-2+1]+2ih2*(y[2i2  ,2,4j2-2+1]+y[2i2  ,2,4j2-2+2]+y[2i2  ,2,4j2-2  ]+y[2i2-1,2,4j2-2+1]-4.0*y[2i2,1  ,4j2-2+1]))

#             # rlc.data[i2  ,1,j] = restrict_cell(x,y,2i2,2j-2,ih2)
#             # rlc.data[i2+1,2,j] = restrict_cell(x,y,2i2+1,2j-2,ih2)
#         end
#         # j=2j2+1
#         # fj=2j-2
#         # @fastmath  @simd for i2=1:nrm1s2
#         #     i=2i2
#         #     fi=2i-2
#         #     rlc.data[i2,2,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#         #     i=2i2+1
#         #     fi=2i-2
#         #     rlc.data[i2+1,1,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#         #     # @inbounds rlc[i,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#         # end

#     end

#     # @simd for j=2:ncols-1
#     #     fj=2j-2

#     #     for i2=1:nrm1s2
#     #         i=2i2
#     #         fi=2i-2
#     #         @inbounds rlc[i,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#     #         i=2i2+1
#     #         fi=2i-2
#     #         @inbounds rlc[i,j]=restrict_cell_ifj(x,y,i,fj,ih2)
#     #     end

#     # end

#     # @simd for j2=1:ncm1s2
#     #     fj=j2-2
#     #     for i2=1:nrm1s2
#     #         fi=i2-2
#     #         @inbounds rlc[i2,j2]=restrict_cell_ifj(x,y,i2,fj,ih2)
#     #         fi=i2-1
#     #         @inbounds rlc[i2+1,j2]=restrict_cell_ifj(x,y,i2+1,fj,ih2)

#     #     end
#     #     fj=j2-1
#     #     for i2=1:nrm1s2
#     #         fi=i2-2
#     #         @inbounds rlc[i2,j2+1]=restrict_cell_ifj(x,y,i2,fj,ih2)
#     #         fi=i2-1
#     #         @inbounds rlc[i2+1,j2+1]=restrict_cell_ifj(x,y,i2+1,fj,ih2)

#     #     end
#     # end
# end

#function that call the benchmarks
function benchmark_restrict_rb(nf)
    @assert iseven(nf)
    # initialize two 2D arrays (fine and coarse)
    nc=nf>>>1
    Random.seed!(1234);

    slf=rand(nf+2,nf+2)
    rlf=rand(nf+2,nf+2)
    rlc=rand(nc+2,nc+2)
    ih2=rand()
    # the following blocks check if smooth_seq and smooth_seq_avx
    # return the same result
    slfrb=RBArray(slf)
    rlfrb=RBArray(rlf)
    rlcrb=RBArray(rlc)
    restrict_residual!(rlc,slf,rlf,ih2)
    restrict_residual!(rlcrb,slfrb,rlfrb,ih2)
    res=norm(rlc-Array{Float64,2}(rlcrb))
    @show res
    # #The actual timings
    tsimd=@belapsed restrict_residual!($rlc,$slf,$rlf,$ih2)
    trb=@belapsed restrict_residual!($rlcrb,$slfrb,$rlfrb,$ih2)
    @show "restrict",nf,res,tsimd,trb,tsimd/trb
    nothing
end

# for n in [32,1024]
#     benchmark_smooth_gs_rb(n)
# end

for n in [32,1024]
    benchmark_interpolate_rb(n)
end

# for n in [32,128]
#     benchmark_restrict_rb(n)
# end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 65beba5e-776d-40d5-96a6-d0a2f3a9a4ea
using Pkg; Pkg.activate("/Users/jj/JuliaProjects/LaserLab/")

# ╔═╡ 8a67d331-35ae-44d6-942c-a61db4e7afa4
begin
	using PlutoUI
	using CSV
	using DataFrames
	using Images
	using Colors
	using Plots
	using Printf
	using Interpolations
	using QuadGK
	using Markdown
	using InteractiveUtils
	using LsqFit
	using Statistics
	using StatsBase
	import Glob
end

# ╔═╡ cbac975d-cda5-43fc-bd22-010ce10e8f92

# ╔═╡ d1b57350-bd9d-11ec-21f4-3d7a8656a31c
import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W,
    A, N, mol, mmol, V, L, M

# ╔═╡ 53943a33-9cf3-4c9f-add5-278fdc3bc0fe
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ b829dbf6-486e-4bd1-914c-54fe6d389b9c
lfi = ingredients("../src/LaserLab.jl")

# ╔═╡ 7aafa82d-4fe3-4bdc-a14c-ec48a7178260
PlutoUI.TableOfContents(title="G2 Bold73 analysis", indent=true)

# ╔═╡ 9da70d55-28b9-4c12-ba6a-c560832a5194
md"""
# Read Data
"""

# ╔═╡ 6502fdee-cd60-4d32-ac77-c42c845b12ba
begin
	proot ="/Users/jj/JuliaProjects/LaserLab/labdata"
	sroot = "/Users/jj/JuliaProjects/LaserLab/data/G2Sl"
	psample = "G2_BOLD_073_A1"
	prun = "G2SL_20220407"
	path = joinpath(proot,psample, prun)
	pfiles ="*.dat"
	fls = joinpath(path,pfiles)
	nfiles = length(fls)
	sflts = string.(collect(2:11));
	sreps = string.(collect(0:10));
	snum = ["1", "2", "c"];
	stypes=["Dark", "Imag"];
md"
- Path for glob: $path
- files : $fls
- number of files: $nfiles
"
end

# ╔═╡ dd51775e-a254-49da-9062-c166ec4a0922
files = Glob.glob(pfiles, path);

# ╔═╡ 6fade49e-845d-4844-9bbc-748f8c958432
typeof(files)

# ╔═╡ 02e2aa4a-40ff-499b-ab31-1004023b5faa
names = [split(f,"/")[end] for f in files];

# ╔═╡ a2ce79cf-a447-4cab-8843-a2cd691b4004
md"""
# Select Filter and Rep
"""

# ╔═╡ 90bda847-8afe-4968-9a02-50ad7c7cc572
md""" Select Filter number (1-10 if 0 add all): $(@bind whichf Select(sflts))"""

# ╔═╡ 1cd87997-3ff3-40fa-98b7-ee9655848a76
md""" Select repetition number (1-10 if 0 first): $(@bind whichr Select(sreps))"""

# ╔═╡ 01aeef0d-12e1-4e32-974f-f5405a661d95
md"""
# Reconstruct spectra
"""

# ╔═╡ 8feb58fd-3cb7-4860-bc22-52bcee93d40a
scs = [string(c) for c in collect(5:10)]

# ╔═╡ 18a0121f-b0a5-46e5-8178-9bc800476b65
@bind zrepx CheckBox()

# ╔═╡ 3aa1788a-131a-402f-9493-32d353020219
fdf = lfi.LaserLab.load_df_from_csv(sroot, "G2Bold073A1Fluo.csv", lfi.LaserLab.enG)

# ╔═╡ 8740f867-109f-47c1-9e90-57ece3b83886

# ╔═╡ 8c52eac3-388f-49ea-bcc4-6d9a42085c52
begin
	plot(collect(1:10), fdf.r1, lw=2, label="rep5 S")
	plot!(collect(1:10), fdf.r2, lw=2, label="rep6 S")
	plot!(collect(1:10), fdf.r3, lw=2, label="rep7 S")
	plot!(collect(1:10), fdf.r4, lw=2, label="rep8 S")
	plot!(collect(1:10), fdf.r5, lw=2, label="rep9 S")
	plot!(collect(1:10), fdf.r6, lw=2, label="rep10 S")
end

# ╔═╡ cec73687-de04-46e6-a73e-a4d1d3beafc5
md"# Functions"

# ╔═╡ 6fdaa0d2-b6af-40ca-b8e2-5cb8176bbadb
function select_image(inames::Vector{String}, fnumber::String, repetition::String,
	                  itype::String="Imag", inum::String="1")

	names = [split(f,"/")[end] for f in inames]
	if repetition == "0"
		fhead = string("Filter_", fnumber, "_ExpoT_1000000ms_")
	else
		fhead = string("Filter_", fnumber, "_rep_",repetition,"_ExpoT_1000000ms_")
	end

	if itype == "Dark"
		ftail = string(itype, "_", inum, ".dat")
	elseif itype == "Imag" && inum=="1"
		ftail = string("Imag_1.dat")
	elseif itype == "Imag" && inum=="c"
		ftail = string("Imag_corrected.dat")
	else
		@info "not implemented"
	end

	ff = string(fhead, ftail)
	indx = findall([occursin(ff, name) for name in names])[1]
	inames[indx]
end

# ╔═╡ edabab28-72bb-43f2-b90b-7b6ba939f2eb
function dftomat(img::DataFrame)
	mimg = Array{Float64}(undef, size(img)...);
	for i in 1:size(img)[1]
		for j in 1:size(img)[2]
			if img[i,j] >= 0.0 && !isnan(img[i,j])
    			mimg[i, j] = img[i,j]
			else
				mimg[i, j] = 0.0
			end
		end
	end
	mimg
end

# ╔═╡ b8af0d87-f9e1-4ea7-9fa0-1c675233a8d3
function findpeaks(img::Matrix{Float64}; xrange::Int64=10)
	maxx, indxt = findmax(img)
	ixmin = max(1, indxt[1] - xrange)
	ixmax = min(size(img)[1], indxt[1] + xrange -1)
	iymin = max(1, indxt[2] - xrange)
	iymax = min(size(img)[2], indxt[2] + xrange -1)

	img2 = Array{Float64}(undef, ixmax-ixmin+1, iymax-iymin+1);

	stot = 0.0
	for (i, ix) in enumerate(ixmin:ixmax)
		for (j, iy) in enumerate(iymin:iymax)
			img2[i,j]=img[ix,iy]
			stot += img[ix,iy]
		end
	end
	img2, stot
end

# ╔═╡ 5dc25eb7-7841-4c04-aa5e-7d63acc5a312
function get_images(files::Vector{String}, whichf::String, whichr::String)
	ffimg = select_image(files, whichf, whichr, "Imag", "1");
	imgdf = DataFrame(CSV.File(ffimg, header=false,delim="\t"));
	imgm  = dftomat(imgdf)

	ffdrk = select_image(files, whichf, whichr, "Dark", "1");
	drkdf = DataFrame(CSV.File(ffdrk, header=false,delim="\t"));
	drkm  = dftomat(drkdf)

	ffdrk2 = select_image(files, whichf, whichr, "Dark", "2");
	drkdf2 = DataFrame(CSV.File(ffdrk2, header=false,delim="\t"));
	sdrk   = drkdf .- drkdf2
	sdrkm  = dftomat(sdrk)

	simgdf = imgdf .- drkdf
	simgm  = dftomat(simgdf)

	imgm, drkm, sdrkm, simgm
end

# ╔═╡ 21f71c85-9e22-4834-a7d7-76e3d1ab681b
imgm, drkm, sdrkm, simgm = get_images(files, whichf, whichr);

# ╔═╡ 4aab4565-280c-432c-8f6d-6774a4dbe8be
begin
	simgmn = simgm ./maximum(simgm);
	gimg = Gray.(simgmn);
	plot(gimg, layout=(1,1), titlefontsize=10)
end

# ╔═╡ 6d08fd5c-6757-43dd-a0e3-bb161e78cd26
pimgx, stot = findpeaks(simgm; xrange=30);

# ╔═╡ f56df91d-53ce-4ca7-a9c8-7740cc3bc430
maxx, indxt = findmax(pimgx);

# ╔═╡ 3444eb77-b397-4ac3-99a2-972abc97090d
begin
	pimgn = pimgx ./maximum(pimgx);
	gimg2  =Gray.(pimgn);
	plot(gimg2,layout=(1,1), titlefontsize=10)
end

# ╔═╡ 17eac4d7-ea71-4205-87d1-e239317c2d7e
typeof(simgm)

# ╔═╡ 64ed6a8c-1bca-4c19-83ff-d507ba4359a5
sdrkmn = sdrkm ./maximum(sdrkm);

# ╔═╡ 1313a166-29f3-4f5e-8e79-dbfab22a69d9
begin
	gdrk = Gray.(sdrkmn)
	plot(gdrk, layout=(1,1), titlefontsize=10)
end

# ╔═╡ e6816855-d168-4317-9ddc-8e6b13c1d0a1
pimgd, stotd = findpeaks(sdrkm; xrange=30);

# ╔═╡ a9d90521-a7e5-4cd2-92ca-45de8c03e514
maxd, indxd = findmax(pimgd)

# ╔═╡ 55a34a0a-ae88-4b80-b6e1-515da142faa8
thrz = maxd + sqrt(maxd)

# ╔═╡ 87369d14-7318-417c-bc33-6a1bc807c0ee
begin
	pimgdn = pimgd ./maximum(pimgd)
	gimgd  =Gray.(pimgdn)
	plot(gimgd, layout=(1,1), titlefontsize=10)
end

# ╔═╡ b776cef6-90d6-4897-bd6e-b81a3a58b416
md"""
- Filter = $whichf
- repetition = $whichr
- **signal**
- value of max = $maxx
- position of max = $indxt
- Total light in spot (DC subtracted) = $(round(stot, sigdigits=3))
- **DC**
- value of max = $maxd
- position of max = $indxd
- DC in spot  = $(round(stotd, sigdigits=3))
"""

# ╔═╡ f41cf5c7-d208-40f9-a4ee-5e992f73ca36
function recospec(files::Vector{String}; rep::String="5", xrange::Int64=30)
	STOT = []
	MAX  = []
	IX   = []
	IY   = []
	STOTd = []
	MAXd  = []
	IXd  = []
	IYd   = []
	for flt in string.(collect(2:11))
		imgm, drkm, sdrkm, simgm = get_images(files, flt, rep)
		pimgx, stot = findpeaks(simgm; xrange=xrange)
		maxx, indxt = findmax(pimgx)
		pimgd, stotd = findpeaks(sdrkm; xrange=xrange)
		maxd, indxd = findmax(pimgd)
		push!(STOT, stot)
		push!(MAX, maxx)
		push!(IX, indxt[1])
		push!(IY, indxt[2])
		push!(STOTd, stotd)
		push!(MAXd, maxd)
		push!(IXd, indxd[1])
		push!(IYd, indxd[2])
	end
	Dict("stot" =>STOT, "smax"=>MAX, "six"=>IX,"siy"=>IY,
	     "dtot" =>STOTd, "dmax"=>MAXd, "dix"=>IXd,"diy"=>IYd)
end

# ╔═╡ 2c71142a-34b2-4804-9af7-dc8147ab30a3
if zrepx
	repx = [recospec(files; rep=c, xrange=30) for c in scs]
	rp = [repx[i]["stot"] for i in 1:6]

	fluodf = DataFrame("r1" => repx[1]["stot"],
		               "r2" => repx[2]["stot"],
		               "r3" => repx[3]["stot"],
		               "r4" => repx[4]["stot"],
		               "r5" => repx[5]["stot"],
		               "r6" => repx[6]["stot"])

	dfname = joinpath(sroot, "G2Bold073A1Fluo.csv")
	CSV.write(dfname, fluodf)

	plot(collect(1:10), repx[1]["stot"], lw=2, label="rep5 S")
	plot!(collect(1:10), repx[1]["dtot"], label="rep5 D")
	plot!(collect(1:10), repx[2]["stot"], lw=2, label="rep6 S")
	plot!(collect(1:10), repx[3]["stot"], lw=2, label="rep7 S")
	plot!(collect(1:10), repx[4]["stot"], lw=2, label="rep8 S")
	plot!(collect(1:10), repx[5]["stot"], lw=2, label="rep9 S")
	plot!(collect(1:10), repx[6]["stot"], lw=2, label="rep10 S")
end

# ╔═╡ fca14add-ed54-449b-82a3-bffd03f88cdd
function plot_images(imgm, drkm, sdrkm, simgm)
	imghm = heatmap(imgm, fill_z=imgm, title=string("rep=", whichr, " filt=", whichf, " Img" ), titlefontsize=10)

	drkhm = heatmap(drkm, title=string("rep=", whichr, " filt=", whichf, " Dark" ), titlefontsize=10)

	sdrkhm = heatmap(sdrkm, title=string("rep=", whichr, " filt=", whichf, "Dark - dark2" ), titlefontsize=10)

	simghm = heatmap(simgm, title=string("rep=", whichr, " filt=", whichf, "Img - Dark" ))

	plot(imghm, drkhm, sdrkhm, simghm, layout=(2,2), titlefontsize=10)
end

# ╔═╡ dfc2478c-9834-427d-a976-38068b9e59fc
plot_images(imgm, drkm, sdrkm, simgm)

# ╔═╡ 4480a6a1-3299-43b9-b871-0b9f525577e5
function img_stats(img)
	xmean = describe(img).mean
	xmax = describe(img).max
	xmin = describe(img).min
	hxmean = lfi.LaserLab.hist1d(xmean, 50, minimum(xmean), maximum(xmean))
	hxmin = lfi.LaserLab.hist1d(xmin, 50, minimum(xmin), maximum(xmin))
	hxmax = lfi.LaserLab.hist1d(xmax, 50, minimum(xmax), maximum(xmax))
	(xmean = xmean, xmin = xmin, xmax = xmax,
	hxmean = hxmean, hxmin = hxmin, hxmax = hxmax)
end

# ╔═╡ 454cab5c-e715-46bd-a799-0e9e2e9a7f67
function img_title(label)
	string(label,"_Filter_", whichf, "_rep_", whichr, "_type_", whicht, "_number_", whichn)
end

# ╔═╡ 4c1c878c-fcfb-4069-a27a-5b1288896955
with_terminal() do
	for (i, ix) in enumerate(10:12)
		for (j, iy) in enumerate(20:22)
			#println("i = ", i, " j = ",j ," ix = ", ix, " iy = ", iy)
		end
	end
end

# ╔═╡ Cell order:
# ╠═8a67d331-35ae-44d6-942c-a61db4e7afa4
# ╠═cbac975d-cda5-43fc-bd22-010ce10e8f92
# ╠═65beba5e-776d-40d5-96a6-d0a2f3a9a4ea
# ╠═d1b57350-bd9d-11ec-21f4-3d7a8656a31c
# ╠═53943a33-9cf3-4c9f-add5-278fdc3bc0fe
# ╠═b829dbf6-486e-4bd1-914c-54fe6d389b9c
# ╠═7aafa82d-4fe3-4bdc-a14c-ec48a7178260
# ╠═9da70d55-28b9-4c12-ba6a-c560832a5194
# ╠═6502fdee-cd60-4d32-ac77-c42c845b12ba
# ╠═dd51775e-a254-49da-9062-c166ec4a0922
# ╠═6fade49e-845d-4844-9bbc-748f8c958432
# ╠═02e2aa4a-40ff-499b-ab31-1004023b5faa
# ╠═a2ce79cf-a447-4cab-8843-a2cd691b4004
# ╠═90bda847-8afe-4968-9a02-50ad7c7cc572
# ╠═1cd87997-3ff3-40fa-98b7-ee9655848a76
# ╠═21f71c85-9e22-4834-a7d7-76e3d1ab681b
# ╠═dfc2478c-9834-427d-a976-38068b9e59fc
# ╠═4aab4565-280c-432c-8f6d-6774a4dbe8be
# ╠═6d08fd5c-6757-43dd-a0e3-bb161e78cd26
# ╠═17eac4d7-ea71-4205-87d1-e239317c2d7e
# ╠═f56df91d-53ce-4ca7-a9c8-7740cc3bc430
# ╠═3444eb77-b397-4ac3-99a2-972abc97090d
# ╠═64ed6a8c-1bca-4c19-83ff-d507ba4359a5
# ╠═1313a166-29f3-4f5e-8e79-dbfab22a69d9
# ╠═e6816855-d168-4317-9ddc-8e6b13c1d0a1
# ╠═a9d90521-a7e5-4cd2-92ca-45de8c03e514
# ╠═55a34a0a-ae88-4b80-b6e1-515da142faa8
# ╠═87369d14-7318-417c-bc33-6a1bc807c0ee
# ╠═b776cef6-90d6-4897-bd6e-b81a3a58b416
# ╠═01aeef0d-12e1-4e32-974f-f5405a661d95
# ╠═8feb58fd-3cb7-4860-bc22-52bcee93d40a
# ╠═18a0121f-b0a5-46e5-8178-9bc800476b65
# ╠═2c71142a-34b2-4804-9af7-dc8147ab30a3
# ╠═3aa1788a-131a-402f-9493-32d353020219
# ╠═8740f867-109f-47c1-9e90-57ece3b83886
# ╠═8c52eac3-388f-49ea-bcc4-6d9a42085c52
# ╠═cec73687-de04-46e6-a73e-a4d1d3beafc5
# ╠═6fdaa0d2-b6af-40ca-b8e2-5cb8176bbadb
# ╠═edabab28-72bb-43f2-b90b-7b6ba939f2eb
# ╠═b8af0d87-f9e1-4ea7-9fa0-1c675233a8d3
# ╠═f41cf5c7-d208-40f9-a4ee-5e992f73ca36
# ╠═5dc25eb7-7841-4c04-aa5e-7d63acc5a312
# ╠═fca14add-ed54-449b-82a3-bffd03f88cdd
# ╠═4480a6a1-3299-43b9-b871-0b9f525577e5
# ╠═454cab5c-e715-46bd-a799-0e9e2e9a7f67
# ╠═4c1c878c-fcfb-4069-a27a-5b1288896955

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab


text:
include("settings.jl")

makedocs(sitename = sitename; settings...)


text:
function Dims_create!(nnodes::Integer, ndims::Integer, dims::MPIBuffertype{T}) where {T <: Integer}
    # int MPI_Dims_create(int nnodes, int ndims, int dims[])
    @mpichk ccall((:MPI_Dims_create, libmpi), Cint,
                  (Cint, Cint, Ptr{Cint}), nnodes, ndims, dims)
end

function Dims_create!(nnodes::Integer, dims::AbstractArray{T,N}) where {T<:Integer, N}
    cdims = Cint.(dims[:])
    ndims = length(cdims)
    Dims_create!(nnodes, ndims, cdims)
    dims[:] .= cdims
end

function Cart_create(comm_old::Comm, ndims::Integer, dims::MPIBuffertype{Cint}, periods::MPIBuffertype{Cint}, reorder)
    comm_cart = Comm()
    # int MPI_Cart_create(MPI_Comm comm_old, int ndims, const int dims[],
    #                     const int periods[], int reorder, MPI_Comm *comm_cart)
    @mpichk ccall((:MPI_Cart_create, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{MPI_Comm}),
                  comm_old, ndims, dims, periods, reorder, comm_cart)
    if comm_cart.val != MPI_COMM_NULL
        refcount_inc()
        finalizer(free, comm_cart)
    end
    comm_cart
end

function Cart_create(comm_old::Comm, dims::AbstractArray{T,N}, periods::Array{T,N}, reorder) where T <: Integer where N
    cdims    = Cint.(dims[:])
    cperiods = Cint.(periods[:])
    ndims    = length(cdims)
    Cart_create(comm_old, ndims, cdims, cperiods, reorder)
end

function Cart_coords!(comm::Comm, rank::Integer, maxdims::Integer, coords::MPIBuffertype{Cint})
    # int MPI_Cart_coords(MPI_Comm comm, int rank, int maxdims, int coords[])
    @mpichk ccall((:MPI_Cart_coords, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}),
                  comm, rank, maxdims, coords)
end

function Cart_coords(comm::Comm, maxdims::Integer)
    ccoords = Vector{Cint}(undef, maxdims)
    rank    = Comm_rank(comm)
    Cart_coords!(comm, rank, maxdims, ccoords)
    Int.(ccoords)
end

function Cart_shift(comm::Comm, direction::Integer, disp::Integer)
    rank_source = Ref{Cint}()
    rank_dest   = Ref{Cint}()
    # int MPI_Cart_shift(MPI_Comm comm, int direction, int disp,
    #                    int *rank_source, int *rank_dest)
    @mpichk ccall((:MPI_Cart_shift, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}, Ptr{Cint}),
                  comm, direction, disp, rank_source, rank_dest)
    Int(rank_source[]), Int(rank_dest[])
end


text:
# This file includes modified source code from https://github.com/odow/SDDP.jl
# as at 5b3ba3910f6347765708dd6e7058e2fcf7d13ae5

#  Copyright 2017, Oscar Dowson and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################

rhs   = [1, 2, 3]
duals = [0.5, 1.5, 2.5]
newrhs = [1.1, 1.2, 1.3]
bounds = [1.0, 10]

@testset "Single element" begin
    @testset "Binary variable" begin
        # Binary states should do exactly the same things as states in SDDP.jl
        m = SDDP.Subproblem()
        @binarystate(m, x, x0 == rhs[1], Bin)
        @test SDDP.nstates(m) == 1
        @binarystate(m, y >= 0, y0 == rhs[2])
        @test SDDP.nstates(m) == 2
        @binarystate(m, 1 <= z <= 10, z0 == rhs[3])
        @test SDDP.nstates(m) == 3
        v = [x, y, z]
        v0 = [x0, y0, z0]
        m.linconstrDuals = duals
        for i in 1:3
            s = SDDP.states(m)[i]
            @test s.variable == v[i]
            c = m.linconstr[s.constraint.idx]
            @test c.lb == c.ub == rhs[i]
            SDDP.setvalue!(s, newrhs[i])
            @test c.lb == c.ub == newrhs[i]
            @test c.terms == JuMP.AffExpr(v0[i])
            @test JuMP.getdual(s) == duals[i]
        end
    end

    @testset "Integer variable" begin
        m = SDDP.Subproblem()

        @binarystate(m, y <= 10, y0 == 10, Int)
        @test isa(y, JuMP.Variable)
        @test isa(y0, JuMP.Variable)
        @test SDDP.nstates(m) == 4
    #     # when moved from set sddip solver, test ub/lb to be 0/1
        @test isapprox(m.colUpper[y.col], 10.0)
        @test m.colLower[y.col] == -Inf

        @binarystate(m, 1 <= z <= 10, z0 == 3, Int)
        @test SDDP.nstates(m) == 4 + 4
        @test isapprox(m.colUpper[z.col], 10.0)
        @test isapprox(m.colLower[z.col], 1.0)

        @binarystate(m, bounds[1] <= w <= bounds[2], w0 == rhs[3], Int)
        @test SDDP.nstates(m) == 4 + 4 + 4
        @test isapprox(m.colUpper[w.col], 10.0)
        @test isapprox(m.colLower[w.col], 1.0)

        @binarystate(m, bounds[1] <= x <= bounds[2] + 4, x0 == rhs[3] - rhs[1], Int)
        @test SDDP.nstates(m) == 4 + 4 + 4 + 4
        @test isapprox(m.colUpper[x.col], 14.0)
        @test isapprox(m.colLower[x.col], 1.0)
    end

    @testset "Continuous variable" begin
        m = SDDP.Subproblem()
        @binarystate(m, x <= 10, x0 == 10, Cont)
        @test length(SDDP.states(m)) == 7
        @binarystate(m, y <= 10, y0 == 10, Cont, 0.01)
        @test SDDP.nstates(m) == 7 + 10
    end
end

@testset "Vector elements" begin
    @testset "Binary variable" begin
        m = SDDP.Subproblem()
        @binarystate(m, bounds[1] <= x[i=1:3] <= bounds[2], x0 == rhs[i], Bin)
        @test SDDP.nstates(m) == 3

        # Binary states should do exactly the same things as states in SDDP.jl
        v = [x[1], x[2], x[3]]
        v0 = [x0[1], x0[2], x0[3]]
        for i in 1:3
            s = SDDP.states(m)[i]
            @test s.variable == v[i]
            c = m.linconstr[s.constraint.idx]
            @test c.lb == c.ub == rhs[i]
            @test c.terms == JuMP.AffExpr(v0[i])
        end
    end

    @testset "Integer variable" begin
        m = SDDP.Subproblem()
        @binarystate(m, x[i=1:3] <= 10, x0 == rhs[i], Int)
        @test isa(x, Vector{JuMP.Variable})
        @test isa(x0, Vector{JuMP.Variable})
        @test SDDP.nstates(m) == 4 * 3
        m.solver = GLPKSolverLP()
        JuMP.solve(m)
        @test isapprox(JuMP.getvalue(x0), rhs)

        letters = [:a, :b]
        idxs = Dict(zip(letters, [1; 2]))
        @binarystate(m, y[i=1:3, j=letters] <= 10, y0 == ones(3, 2)[i, idxs[j]], Int)
        @test isa(y, JuMP.JuMPArray)
        @test isa(y0, JuMP.JuMPArray)
        @test SDDP.nstates(m) == 4 * 3 + 4 * 6
        JuMP.solve(m)
        @test isapprox(JuMP.getvalue(y0).innerArray, ones(3, 2))

    end

    @testset "Continuous variable" begin
        m = SDDP.Subproblem()
        @binarystate(m, 1 <= x[i=1:3] <= 10, x0 == rhs[i], Cont)
        @test SDDP.nstates(m) == 7 * 3
        @binarystate(m, y[i=1:3] <= 10, y0 == rhs[i], Cont, 0.01)
        @test length(SDDP.states(m)) == 7 * 3 + 10 * 3
        m.solver = GLPKSolverLP()
        JuMP.solve(m)
        @test isapprox(JuMP.getvalue(y0), rhs)
    end

    @testset "@binarystates" begin

        m = SDDP.Subproblem()
        @binarystates(m, begin
            x[i=1:3] <= bounds[2], x0 == rhs[i], Bin
            y <= 1, y0==0.5, Cont, 0.01
            z <= 5, z0==5, Int
        end)

        # 3 state variables for x, nstates(1*0.01)=7 for y, and nstates(5)=3 for z
        @test SDDP.nstates(m) == 3 + 7 + 3

        v = [x[1], x[2], x[3]]
        v0 = [x0[1], x0[2], x0[3]]
        for i in 1:3
            s = SDDP.states(m)[i]
            @test s.variable == v[i]
            c = m.linconstr[s.constraint.idx]
            @test c.lb == c.ub == rhs[i]
            @test c.terms == JuMP.AffExpr(v0[i])
        end
        for v in [y, y0, z, z0]
            @test isa(v, JuMP.Variable)
        end
        m.solver = GLPKSolverLP()
        JuMP.solve(m)
        @test isapprox(JuMP.getvalue(y0), 0.5)
        @test isapprox(JuMP.getvalue(z0), 5)
    end
end


text:
function print_colored(in::String,color1)
	if (VERSION.major==0&&VERSION.minor<=6)
		return print_with_color(color1,in)
	else
		return printstyled(in,color=color1)
	end
end

test_list = ["testRealNumbers.jl";"testComplexNumbers.jl";"testForwardDiff.jl";"testDual_.jl";"testHyperDualNumbers.jl";"testDates.jl"]

println("Running tests:\n")
for (current_test,i) in zip(test_list,1:length(test_list))
	println("------------------------------------------------------------")
    println("  * $(current_test) *")
    include(current_test)
	println("------------------------------------------------------------")
	if (i<length(test_list))
		println("")
	end
end


text:
module Seismic
    using Grid,Requires,Compat
    include("ReadWrite/ReadWrite.jl")
    include("Utils/Utils.jl")
    include("Processing/Processing.jl")
    include("Imaging/Imaging.jl")
    include("Modelling/Modelling.jl")
    include("Operators/Operators.jl")
    include("Solvers/Solvers.jl")
    include("Tools/Tools.jl")
    include("Wavelets/Wavelets.jl")
    include("Windows/Windows.jl")
    @require PyPlot include("Plotting/Plotting.jl")
end


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

const ALEInterface = Void
const ALEPtr = Ptr{ALEInterface}
const ALEState = Void
const ALEStatePtr = Ptr{ALEState}

ALE_new() = ccall((:ALE_new, libale_c), ALEPtr, ())
ALE_del(ale::ALEPtr) = ccall((:ALE_del, libale_c), Void, (ALEPtr,), ale)

function getString(ale::ALEPtr, key::String)
    res = ccall((:getString, libale_c), Ptr{Cchar}, (ALEPtr, Ptr{Cchar}),
        ale, key)
    bytestring(res)
end
getInt(ale::ALEPtr, key::String) = ccall((:getInt, libale_c), Cint,
    (ALEPtr, Ptr{Cchar}), ale, key)
getBool(ale::ALEPtr, key::String) = ccall((:getBool, libale_c), Cint,
    (ALEPtr, Ptr{Cchar}), ale, key) > 0
getFloat(ale::ALEPtr, key::String) = ccall((:getFloat, libale_c), Float32,
    (ALEPtr, Ptr{Cchar}), ale, key)

setString(ale::ALEPtr, key::String, value::String) =
    ccall((:setString, libale_c), Void, (ALEPtr, Ptr{Cchar}, Ptr{Cchar}),
        ale, key, value)
setInt(ale::ALEPtr, key::String, value::Cint) = ccall((:setInt, libale_c),
    Void, (ALEPtr, Ptr{Cchar}, Cint), ale, key, value)
setBool(ale::ALEPtr, key::String, value::Bool) = ccall((:setBool, libale_c),
    Void, (ALEPtr, Ptr{Cchar}, Cint), ale, key, value)
setFloat(ale::ALEPtr, key::String, value::Cfloat) =
    ccall((:setFloat, libale_c), Void, (ALEPtr, Ptr{Cchar}, Cfloat),
        ale, key, value)

loadROM(ale::ALEPtr, rom_file::String) = ccall((:loadROM, libale_c), Void,
    (ALEPtr, Ptr{Cchar}), ale, rom_file)

act(ale::ALEPtr, action::Cint) =
    ccall((:act, libale_c), Cint, (ALEPtr, Cint), ale, action)
game_over(ale::ALEPtr) =
    ccall((:game_over, libale_c), Cint, (ALEPtr,), ale) > 0
reset_game(ale::ALEPtr) = ccall((:reset_game, libale_c), Void, (ALEPtr,), ale)

function getLegalActionSet(ale::ALEPtr)
    actions = Array{Cint}(0)
    getLegalActionSet!(ale, actions)
    actions
end
function getLegalActionSet!(ale::ALEPtr, actions::Vector{Cint})
    resize!(actions, getLegalActionSize(ale))
    ccall((:getLegalActionSet, libale_c), Void, (ALEPtr, Ptr{Cint}),
        ale, actions)
end
getLegalActionSize(ale::ALEPtr) =
    ccall((:getLegalActionSize, libale_c), Cint, (ALEPtr,), ale)

function getMinimalActionSet(ale::ALEPtr)
    actions = Array{Cint}(0)
    getMinimalActionSet!(ale, actions)
    actions
end
function getMinimalActionSet!(ale::ALEPtr, actions::Vector{Cint})
    resize!(actions, getMinimalActionSize(ale))
    ccall((:getMinimalActionSet, libale_c), Void, (ALEPtr, Ptr{Cint}),
        ale, actions)
end
getMinimalActionSize(ale::ALEPtr) =
    ccall((:getMinimalActionSize, libale_c), Cint, (ALEPtr,), ale)

getFrameNumber(ale::ALEPtr) =
    ccall((:getFrameNumber, libale_c), Cint, (ALEPtr,), ale)
lives(ale::ALEPtr) = ccall((:lives, libale_c), Cint, (ALEPtr,), ale)
getEpisodeFrameNumber(ale::ALEPtr) =
    ccall((:getEpisodeFrameNumber, libale_c), Cint, (ALEPtr,), ale)

function getScreen(ale::ALEPtr)
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    screen_data = Array{Cuchar}(w*h) # row-major order
    getScreen!(ale, screen_data)
    screen_data
end
getScreen!(ale::ALEPtr, screen_data::Vector{Cuchar}) =
    ccall((:getScreen, libale_c), Void, (ALEPtr, Ptr{Cuchar}), ale, screen_data)
getRAM(ale::ALEPtr, ram::Vector{Cuchar}) =
    ccall((:getRAM, libale_c), Void, (ALEPtr, Ptr{Cuchar}), ale, ram)
getRAMSize(ale::ALEPtr) = ccall((:getRAMSize, libale_c), Cint, (ALEPtr,), ale)
getScreenWidth(ale::ALEPtr) =
    ccall((:getScreenWidth, libale_c), Cint, (ALEPtr,), ale)
getScreenHeight(ale::ALEPtr) =
    ccall((:getScreenHeight, libale_c), Cint, (ALEPtr,), ale)

getScreenRGB(ale::ALEPtr, output_buffer::Vector{Cuchar}) =
    ccall((:getScreenRGB, libale_c), Void, (ALEPtr, Ptr{Cuchar}),
        ale, output_buffer)
getScreenGrayscale(ale::ALEPtr, output_buffer::Vector{Cuchar}) =
    ccall((:getScreenGrayscale, libale_c), Void, (ALEPtr, Ptr{Cuchar}),
        ale, output_buffer)

saveState(ale::ALEPtr) = ccall((:saveState, libale_c), Void, (ALEPtr,), ale)
loadState(ale::ALEPtr) = ccall((:loadState, libale_c), Void, (ALEPtr,), ale)
cloneState(ale::ALEPtr) = ccall((:cloneState, libale_c), ALEStatePtr,
    (ALEPtr,), ale)
restoreState(ale::ALEPtr, state::ALEStatePtr) =
    ccall((:restoreState, libale_c), Void, (ALEPtr, ALEStatePtr), ale, state)
cloneSystemState(ale::ALEPtr) = ccall((:cloneSystemState, libale_c),
    ALEStatePtr, (ALEPtr,), ale)
restoreSystemState(ale::ALEPtr, state::ALEStatePtr) =
    ccall((:restoreSystemState, libale_c), Void, (ALEPtr, ALEStatePtr),
        ale, state)
deleteState(state::ALEStatePtr) = ccall((:deleteState, libale_c), Void,
    (ALEStatePtr,), state)
saveScreenPNG(ale::ALEPtr, filename::String) = ccall((:saveScreenPNG, libale_c),
    Void, (ALEPtr, Ptr{Cchar}), ale, filename)

function encodeState(state::ALEStatePtr)
    len = encodeStateLen(state)
    buf = Array{Cchar}(len)
    ccall((:encodeState, libale_c), Void, (ALEStatePtr, Ptr{Cchar}, Cint),
        state, buf, len)
    buf
end
encodeStateLen(state::ALEStatePtr) = ccall((:encodeStateLen, libale_c),
    Cint, (ALEStatePtr,), state)
decodeState(buf::Array{Cchar,1}) = ccall((:decodeState, libale_c),
    ALEStatePtr, (Ptr{Cchar}, Cint), buf, length(buf))

export
    # Types
    ALEInterface,
    ALEPtr,
    # Functions
    ALE_new,
    ALE_del,
    getInt,
    getBool,
    getFloat,
    setString,
    setInt,
    setBool,
    setFloat,
    loadROM,
    act,
    game_over,
    reset_game,
    getLegalActionSet,
    getLegalActionSet!,
    getLegalActionSize,
    getMinimalActionSet,
    getMinimalActionSet!,
    getMinimalActionSize,
    getFrameNumber,
    lives,
    getEpisodeFrameNumber,
    getScreen,
    getScreen!,
    getRAM,
    getRAMSize,
    getScreenWidth,
    getScreenHeight,
    getScreenRGB,
    getScreenGrayscale,
    saveState,
    loadState,
    cloneState,
    restoreState,
    cloneSystemState,
    restoreSystemState,
    deleteState,
    saveScreenPNG,
    encodeState,
    encodeStateLen,
    decodeState


text:
function run(s)
    # Your code here
    code = readlines(IOBuffer(s))
    seen = Set{Int}()
    curr = 1
    acc = 0
    while !(curr in seen)
        push!(seen, curr)
        line = code[curr]
        op = line[1:3]
        if op == "acc"
            acc += parse(Int, line[5:end])
            curr += 1
        elseif op == "jmp"
            curr += parse(Int, line[5:end])
        else
            curr += 1
        end
    end
    return acc
end

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

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

main()


text:
# ---
# title: 70. Climbing Stairs
# id: problem70
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/climbing-stairs/description/>
# hidden: true
# ---
#
# You are climbing a staircase. It takes `n` steps to reach the top.
#
# Each time you can either climb `1` or `2` steps. In how many distinct ways can
# you climb to the top?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 2
#     Output: 2
#     Explanation: There are two ways to climb to the top.
#     1. 1 step + 1 step
#     2. 2 steps
#
#
# **Example 2:**
#
#
#
#     Input: n = 3
#     Output: 3
#     Explanation: There are three ways to climb to the top.
#     1. 1 step + 1 step + 1 step
#     2. 1 step + 2 steps
#     3. 2 steps + 1 step
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 45`
#
#
## @lc code=start
using LeetCode

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


text:
# COMPARISON & EQUALITY
<(a::DoubleFloat64, b::DoubleFloat64) = a.hi + a.lo < b.hi + b.lo
<(a::DoubleFloat64, b::Float64) = a.hi < b || (a.hi == b) && a.lo < 0.0
<(a::Float64, b::DoubleFloat64) = a < b.hi || (a == b.hi) && b.lo > 0.0

Base.isless(a::DoubleFloat64, b::DoubleFloat64) = isless(a.hi + a.lo, b.hi + b.lo)

<=(a::DoubleFloat64, b::DoubleFloat64) = !(b < a)
<=(a::DoubleFloat64, b::Float64) = !(b < a)
<=(a::Float64, b::DoubleFloat64) = !(b < a)

==(a::DoubleFloat64, b::Float64) = a.hi == b && a.lo == 0.0
==(a::Float64, b::DoubleFloat64) = b == a
Base.iszero(a::DoubleFloat64) = a.hi == 0.0

if VERSION > v"0.7.0-DEV.1319"
    Base.isone(a::DoubleFloat64) = a.hi == 1.0 && a.lo == 0.0
else
    isone(a::DoubleFloat64) = a.hi == 1.0 && a.lo == 0.0
end

Base.abs(a::DoubleFloat64) = a.hi < 0.0 ? -a : a

Base.eps(::DoubleFloat64) = 4.93038065763132e-32 # 2^-104
Base.eps(::Type{DoubleFloat64}) = 4.93038065763132e-32 # 2^-104
Base.eps(::Type{FastDouble}) = 4.93038065763132e-32 # 2^-104
Base.eps(::Type{AccurateDouble}) = 4.93038065763132e-32 # 2^-104

Base.realmin(::Type{DoubleFloat64}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
Base.realmin(::Type{FastDouble}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
Base.realmin(::Type{AccurateDouble}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
Base.realmax(::Type{DoubleFloat64}) = DoubleFloat64(1.79769313486231570815e+308, 9.97920154767359795037e+291);
Base.realmax(::Type{FastDouble}) = DoubleFloat64(1.79769313486231570815e+308, 9.97920154767359795037e+291);
Base.realmax(::Type{AccurateDouble}) = DoubleFloat64(1.79769313486231570815e+308, 9.97920154767359795037e+291);

Base.isnan(a::DoubleFloat64) = isnan(a.hi) || isnan(a.lo)
Base.isinf(a::DoubleFloat64) = isinf(a.hi)
Base.isfinite(a::DoubleFloat64) = isfinite(a.hi)

#
# ROUNDING
#
@inline function Base.round(a::DoubleFloat64)
    hi = round(a.hi)
    lo = 0.0

    if hi == a.hi
        # High word is an integer already.  Round the low word.
        lo = round(a.lo)

        # Renormalize. This is needed if hi = some integer, lo = 1/2.
        hi, lo = quick_two_sum(hi, lo)
    else
        # High word is not an integer.
        if abs(hi-a.hi) == 0.5 && a.lo < 0.0
            # There is a tie in the high word, consult the low word to break the tie.
            hi -= 1.0
        end
    end

     DoubleFloat64(hi, lo)
end

@inline function Base.floor(a::DoubleFloat64)
    hi = floor(a.hi)
    lo = 0.0

    if hi == a.hi
        lo = floor(a.lo)
        hi, lo = quick_two_sum(hi, lo)
    end

    DoubleFloat64(hi, lo)
end

@inline function Base.ceil(a::DoubleFloat64)
    hi = ceil(a.hi)
    lo = 0.0

    if hi == a.hi
        lo = ceil(a.lo)
        hi, lo = quick_two_sum(hi, lo)
    end

    DoubleFloat64(hi, lo)
end

Base.trunc(a::DoubleFloat64) = a.hi ≥ 0.0 ? floor(a) : ceil(a)
Base.isinteger(x::DoubleFloat64) = iszero(x - trunc(x))

function Base.rand(rng::AbstractRNG, S::Type{DoubleFloat64{T}}) where T<:ComputeMode
    u = rand(rng, UInt64)
    f = Float64(u)
    uf = UInt64(f)
    ur = uf > u ? uf - u : u - uf
    DoubleFloat64{T}(5.421010862427522e-20 * f, 5.421010862427522e-20 * Float64(ur))
end
Base.rand(rng::AbstractRNG, ::Type{DoubleFloat64}) = rand(rng, FastDouble)
# Base.rand(rng, ::Type{FastDouble}) = rand(rng, FastDouble)
# Base.rand(rng, ::Type{AccurateDouble}) = rand(rng, AccurateDouble)
Base.rand(::Type{FastDouble}) = rand(Base.Random.GLOBAL_RNG, FastDouble)
Base.rand(::Type{AccurateDouble}) = rand(Base.Random.GLOBAL_RNG, AccurateDouble)
Base.rand(::Type{DoubleFloat64}) = rand(Base.Random.GLOBAL_RNG, FastDouble)

# function Base.rand(rng, T::Type{<:DoubleFloat64}, d1, dims::Vararg{Int, N}) where N
#     rands = Array{T}(d1, dims...)
#     for l in eachindex(rands)
#         rands[l] = rand(rng, T)
#     end
#     rands
# end
function Base.rand(rng::AbstractRNG, T::Type{<:DoubleFloat64}, dims::Vararg{Int, N}) where N
    rands = Array{T}(dims)
    for l in eachindex(rands)
        rands[l] = rand(rng, T)
    end
    rands
end
Base.rand(T::Type{<:DoubleFloat64}, dims::Vararg{Int, N}) where N = rand(Base.Random.GLOBAL_RNG, T, dims)

Base.rand(::Type{Complex{DoubleFloat64}}) = rand(Complex{FastDouble})
Base.rand(::Type{Complex{DoubleFloat64}}, dims::Vararg{Int, N}) where N = rand(Complex{FastDouble}, dims)
Base.rand(rng::AbstractRNG, ::Type{Complex{DoubleFloat64}}) = rand(rng, Complex{FastDouble})
Base.rand(rng::AbstractRNG, ::Type{Complex{DoubleFloat64}}, dims::Vararg{Int, N}) where N = rand(rng, Complex{FastDouble}, dims)

function Base.decompose(a::DoubleFloat64) ::Tuple{Int128, Int, Int}
      hi, lo = a.hi, a.lo
      num1, pow1, den1 = Base.decompose(hi)
      num2, pow2, den2 = Base.decompose(lo)

      num = Int128(num1)

      pdiff = pow1 - pow2
      shift = min(pdiff, 52)
      signed_num = den1 * (Int128(num) << shift) # den1 is +1/-1
      signed_num += den2 * (num2 >> (pdiff - shift)) # den2 is +1/-1

      num = abs(signed_num)
      den = signed_num ≥ 0 ? 1 : -1
      pow = pow1 - shift

      num, pow, den
end


text:
### Day 5 - Problem 1 & 2
## How About a Nice Game of Chess
## Author: Thanasis Georgiou

using Nettle;

# Read input
input = "reyedfim"

# Password
password1 = Array{Char}(8)
password2 = Array{Char}(8)
fill!(password1, ' ')  # 8 Spaces
fill!(password2, ' ')  # 8 Spaces

# Test a couple of hashes
i = 0
while length(find(c -> c == ' ', password1)) != 0 || length(find(c -> c == ' ', password2)) != 0
    i += 1

    # Calculate hash
    md5 = hexdigest("md5", "$(input)$i")
    if md5[1:5] == "00000"
        # First password (in order)
        md5Char = md5[6]
        spaces = find(c -> c == ' ', password1)
        if length(spaces) != 0
            password1[spaces[1]] = md5Char
        end

        # Second password (with index)
        md5Index = parse(Int64, md5[6]) + 1
        md5Char = md5[7]

        if md5Index <= 8 && md5Index >= 1 && password2[md5Index] == ' '
            password2[md5Index] = md5Char
        end
    end
end

println("Password 1: $(password1)")
println("Password 2: $(password2)")