text:
abstract type abstract_ls_info end

include("frac_boundary.jl")
include("move.jl")
include("filter_ls.jl")
include("agg_ls.jl")
include("stable_ls.jl")
include("kkt_ls.jl")

function scale_direction(dir::Class_point, step_size::Float64)
    new_dir = deepcopy(dir)
    new_dir.x *= step_size
    new_dir.y *= step_size
    new_dir.s *= step_size
    new_dir.mu *= step_size
    new_dir.primal_scale *= step_size

    return new_dir
end

function Blank_ls_info()
    this = Class_stable_ls(0.0,0.0, 0, 0.0, 0.0)
    return this
end

function basic_checks(iter::Class_iterate, f_it::Class_iterate, dir::Class_point)
  primal_accurate = LinearAlgebra.norm(eval_jac_prod(f_it, dir.x) - eval_jac_prod(iter, dir.x),Inf) < LinearAlgebra.norm(get_primal_res(iter),Inf) #* 10.0
  dual_accurate = LinearAlgebra.norm(eval_jac_T_prod(f_it, dir.y) - eval_jac_T_prod(iter, dir.y),Inf) < LinearAlgebra.norm(eval_grad_lag(iter,get_mu(iter)),Inf) #* 10.0
  #comp_accurate = -50.0 * get_mu(iter) < minimum(comp_predicted(iter,dir,1.0)) && maximum(comp_predicted(iter,dir,1.0)) < 50.0 * get_mu(iter)
  comp_accurate = LinearAlgebra.norm(comp_predicted(iter,dir,1.0),Inf) < 50.0 * get_mu(iter)

  return comp_accurate #&& primal_accurate && dual_accurate
  #return true
end

function simple_ls(iter::Class_iterate, f_it::Class_iterate, dir::Class_point, accept_type::Symbol, filter::Array{Class_filter,1},  pars::Class_parameters, min_step_size::Float64, timer::class_advanced_timer)
    start_advanced_timer(timer, "SIMPLE_LS")

    # compute fraction to boundary
    lb_s = lb_s_predict(iter, dir, pars)
    step_size_P = simple_max_step(iter.point.s, dir.s, lb_s)
    step_size_D = NaN

    #step_size_P = min(step_size_P,simple_max_step(iter.point.y, dir.y, 0.2 * iter.point.y))

    if accept_type == :accept_stable
      accept_obj = Class_stable_ls(iter, dir, pars)
    elseif accept_type == :accept_aggressive
      accept_obj = Class_agg_ls(iter, dir, pars)
    elseif accept_type == :accept_filter
      accept_obj = Class_filter_ls(iter, dir, pars)
    elseif accept_type == :accept_kkt
      accept_obj = Class_kkt_ls(iter, dir, pars)
    elseif accept_type == :accept_comp
      accept_obj = Class_comp_ls(iter, dir, pars)
    else
      error("acceptance function not defined")
    end

    accept_obj.num_steps = 0

    if !accept_obj.do_ls
        pause_advanced_timer(timer, "SIMPLE_LS")

        return :predict_red_non_negative, iter, accept_obj
    end

    if pars.output_level >= 5
      println("")
      println(pd("α_P"), pd("α_D"), pd("merit_diff"), pd("comp_diff"), pd("phi_diff"), pd("kkt_diff"), pd("dx"), pd("dy"), pd("ds"), pd("status"))
    end

    for i = 1:pars.ls.num_backtracks
      status = :none

      if step_size_P >= min_step_size
        start_advanced_timer(timer,"SIMPLE_LS/move/primal")
        candidate, move_status = move_primal(iter, dir, step_size_P, pars, timer)
        pause_advanced_timer(timer,"SIMPLE_LS/move/primal")

        accept_obj.step_size_P = step_size_P
        accept_obj.num_steps = i

        if move_status == :success
            start_advanced_timer(timer,"SIMPLE_LS/move/dual_bounds")
            lb, ub = dual_bounds(candidate, candidate.point.y, dir.y, pars.ls.comp_feas)
            lb_y_new = lb_y(iter,dir,pars) #min(dir.y .* (dir.s + LinearAlgebra.norm(dir.x,Inf)^2) ./ candidate.point.s, pars.ls.fraction_to_boundary * candidate.point.y)
            ub = min(ub, simple_max_step(candidate.point.y, dir.y, lb_y_new))
            pause_advanced_timer(timer,"SIMPLE_LS/move/dual_bounds")

            if lb >= ub
              move_status = :dual_infeasible
            end
            if !pars.ls.move_primal_seperate_to_dual && !(lb <= step_size_P && step_size_P <= ub)
              move_status = :dual_infeasible
            end
        end

        if move_status == :success
            no_nan = update_J!(candidate, timer, pars) && update_grad!(candidate, timer, pars)
            if no_nan
              start_advanced_timer(timer,"SIMPLE_LS/move/dual")
              candidate, move_status, step_size_D = move_dual(candidate, dir, step_size_P, lb, ub, pars, timer)
              pause_advanced_timer(timer,"SIMPLE_LS/move/dual")
              accept_obj.step_size_D = step_size_D
            else
              move_status = :NaN_ERR
            end
        end

        if move_status != :success
            suggested_step_size_P = step_size_P * pars.ls.backtracking_factor
        end

        if move_status == :success
          start_advanced_timer(timer,"SIMPLE_LS/accept?")
          no_nan = update_obj!(candidate, timer, pars)
          if no_nan
            status, suggested_step_size_P = accept_func!(accept_obj, iter, candidate, dir, step_size_P, filter, pars, timer)
          else
            status = :NaN_ERR
          end
          pause_advanced_timer(timer,"SIMPLE_LS/accept?")
        else
          status = move_status
        end

        @assert(is_updated_correction(iter.cache))

        # SOME DEBUGGING CODE
        if pars.debug_mode >= 1
          obj = copy(get_fval(iter))
          update_obj!(iter, timer, pars)
          #@show get_fval(iter), obj, get_fval(candidate)
          if(obj != get_fval(iter))
              error("this shouldn't happen!!!")
          end
        end

        # OUTPUT INFORMATION ABOUT LINE SEARCH
        if pars.output_level >= 5
          diff = eval_merit_function(candidate, pars) - eval_merit_function(iter, pars)
          comp_diff = LinearAlgebra.norm(comp(candidate),Inf) - LinearAlgebra.norm(comp(iter), Inf)
          mu = iter.point.mu
          phi_diff = eval_phi(candidate, mu) - eval_phi(iter, mu)

          dx = step_size_P * LinearAlgebra.norm(dir.x,2);
          dy = LinearAlgebra.norm(candidate.point.y - iter.point.y,2); ds = norm(candidate.point.s - iter.point.s,2);
          kkt_diff = norm(eval_grad_lag(candidate, get_mu(candidate)),Inf) / LinearAlgebra.norm(eval_grad_lag(iter, get_mu(iter)),Inf)
          println(rd(step_size_P), rd(step_size_D), rd(diff), rd(comp_diff), rd(phi_diff), rd(kkt_diff), rd(dx), rd(dy), rd(ds), pd(status))
        end

        if status == :success
          pause_advanced_timer(timer,"SIMPLE_LS")
          start_advanced_timer(timer,"SIMPLE_LS_SUCCESS")
          pause_advanced_timer(timer,"SIMPLE_LS_SUCCESS")
          if pars.output_level >= 5
              println("--end ls--")
          end

          return :success, candidate, accept_obj
        elseif status == :predict_red_non_negative
          pause_advanced_timer(timer, "SIMPLE_LS")

          if pars.output_level >= 5
              println("--end ls--")
          end
          return status, iter, accept_obj
        end

        #if status == :s_bound
        #  step_size_P *= pars.ls.backtracking_factor # min(pars.ls.backtracking_factor * step_size_P, 1.0 / maximum( (get_s(iter) - get_s(candidate)) ./ get_s(iter) ))
        #else
        #step_size_P *= pars.ls.backtracking_factor
        #end
        step_size_P = suggested_step_size_P
      else
        if pars.output_level >= 5
          println(rd(step_size_P), pd("N/A"), pd("N/A"), pd(:min_α))
        end

        if pars.output_level >= 5
            println("--end ls--")
        end

        pause_advanced_timer(timer, "SIMPLE_LS")
        return :min_α, iter, accept_obj
      end
    end

    if pars.output_level >= 5
        println("--end ls--")
    end

    pause_advanced_timer(timer, "SIMPLE_LS")
    return :max_ls_it, iter, accept_obj
end

function eigenvector_ls(iter::Class_iterate, dir::Class_point, pars::Class_parameters)
    step_size_P = 1.0

    best_candidate = iter;
    intial_val = eval_merit_function(iter, pars)
    best_val = intial_val

    max_it = 10;

    counter_i = 1;
    for i = 1:max_it
      counter_i = i
      candidate_pos, is_feas, step_size_D = move(iter, dir, step_size_P, pars, timer)
      candidate_neg, is_feas, step_size_D = move(iter, dir, -step_size_P, pars, timer)

      better = false

      new_val = eval_merit_function(candidate_pos, pars)
      if new_val < best_val
        best_candidate = candidate_pos
        best_val = new_val
        better = true
      end

      new_val = eval_merit_function(candidate_neg, pars)
      if new_val < best_val
        best_candidate = candidate_neg
        best_val = new_val
        better = true
      end

      if !better
        break
      end

      step_size_P *= 6.0
    end

    @show step_size_P, LinearAlgebra.norm(dir.x), best_val - intial_val

    if counter_i == max_it
      my_warn("max it reached for eig search")
    end

    return :success, best_candidate
end


text:
using Flux, Images
using StatsBase: sample, shuffle

const PATH = joinpath(@__DIR__, "../../res/data/output_train/")
const BP = joinpath.(string(PATH, "bp/"), readdir(string(PATH, "bp/")))
const WP = joinpath.(string(PATH, "wp/"), readdir(string(PATH, "wp/")))
const BR = joinpath.(string(PATH, "br/"), readdir(string(PATH, "br/")))
const WR = joinpath.(string(PATH, "wr/"), readdir(string(PATH, "wr/")))
const BN = joinpath.(string(PATH, "bn/"), readdir(string(PATH, "bn/")))
const WN = joinpath.(string(PATH, "wn/"), readdir(string(PATH, "wn/")))
const BB = joinpath.(string(PATH, "bb/"), readdir(string(PATH, "bb/")))
const WB = joinpath.(string(PATH, "wb/"), readdir(string(PATH, "wb/")))
const BQ = joinpath.(string(PATH, "bq/"), readdir(string(PATH, "bq/")))
const WQ = joinpath.(string(PATH, "wq/"), readdir(string(PATH, "wq/")))
const BK = joinpath.(string(PATH, "bq/"), readdir(string(PATH, "bq/")))
const WK = joinpath.(string(PATH, "wq/"), readdir(string(PATH, "wq/")))
const EMPTY = joinpath.(string(PATH, "empty/"), readdir(string(PATH, "empty/")))

function load_batch(n = 130, nsize = (227,227); path = PATH)
  imgs_paths = shuffle(vcat(
    sample(BP, Int(n/13)),
    sample(WP, Int(n/13)),
    sample(BR, Int(n/13)),
    sample(WR, Int(n/13)),
    sample(BN, Int(n/13)),
    sample(WN, Int(n/13)),
    sample(BB, Int(n/13)),
    sample(WB, Int(n/13)),
    sample(BQ, Int(n/13)),
    sample(WQ, Int(n/13)),
    sample(BK, Int(n/13)),
    sample(WK, Int(n/13)),
    sample(EMPTY, Int(n/13))
  ))
  labels = map(x -> occursin("bp", x) ? 0 :
                    occursin("wp", x) ? 1 :
                    occursin("br", x) ? 2 :
                    occursin("wr", x) ? 3 :
                    occursin("bn", x) ? 4 :
                    occursin("wn", x) ? 5 :
                    occursin("bb", x) ? 6 :
                    occursin("wb", x) ? 7 :
                    occursin("bq", x) ? 8 :
                    occursin("wq", x) ? 9 :
                    occursin("bk", x) ? 10 :
                    occursin("wk", x) ? 11 : 12, imgs_paths)
  labels = Flux.onehotbatch(labels, 0:12)
  imgs = Images.load.(imgs_paths)
  # imgs = map(img -> Images.imresize(img, nsize...), imgs)
  imgs = map(img -> permutedims(channelview(img), (3,2,1)), imgs)
  imgs = cat(imgs..., dims = 4)
  Float32.(imgs), labels
end

ext: jl

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
mutable struct Param <: AbstractParam
    ## number of latent classes
    G::Int

    ## number of empty latent classes
    G0::Int

    ## multinomial model
    xi::Matrix{Float64} # G x (p+1) matrix

    ### longitudinal model
    ## fixed effects
    beta::Matrix{Float64} # G x nfix

    # ## random effects
    mu::Matrix{Float64} # G x nrand

    ## random effects standard errors
    S::Matrix{Float64} ##Vector{Float64} ## q

    # ## random effects correlation matrix Cholesky factor
    L::Array{Float64,3} ##Matrix{Float64} # nrand x nrand

    # ## measurement error standard deviation
    sigma::Float64

    ### survival model
    bl_shape::Vector{Float64} # G
    bl_scale::Vector{Float64} # G
    delta::Matrix{Float64} # G x r

    b::Matrix{Float64}
    cind::Vector{Int}
end

mutable struct ParamSpline <: AbstractParam
    ## number of latent classes
    G::Int

    ## number of empty latent classes
    G0::Int

    ## multinomial model
    xi::Matrix{Float64} # G x (p+1) matrix

    ### longitudinal model
    ## fixed effects
    beta::Matrix{Float64} # G x nfix

    # ## random effects
    mu::Matrix{Float64} # G x nrand

    ## random effects standard errors
    S::Matrix{Float64} ##Vector{Float64} ## q

    # ## random effects correlation matrix Cholesky factor
    L::Array{Float64,3} ##Matrix{Float64} # nrand x nrand

    # ## measurement error standard deviation
    sigma::Float64

    ### survival model
    w::Matrix{Float64} # G x (knots + 4)
    delta::Matrix{Float64} # G x r

    b::Matrix{Float64}
    cind::Vector{Int}
end

Param(N, xi, beta, mu, S, L, sigma, shape, scale, delta) = Param(size(beta,1), 0, xi, beta, mu, S, L, sigma,
                                                                 shape, scale, delta, zeros(eltype(mu), N, size(mu,2)),
                                                                 sample(1:size(mu,1), N))

ParamSpline(N, xi, beta, mu, S, L, sigma, w, delta) = Param(size(beta)[1], 0, xi, beta, mu, S, L, sigma,
                                                            w, delta, zeros(eltype(mu), N, size(mu,2)),
                                                            sample(1:size(mu,1), N))

function Param(data::Data, G::Int; set_zero=false)
    N, p = size(data.Xp)

    ## number of hazard predictors
    r = size(data.Xe,2)

    ## number of fixed effects in LMM model
    l = size(data.Xl,3)

    ## number of random effects in LMM model
    q = size(data.Z,3)

    if set_zero
        Param(G, 0, zeros(Float64, G, p), zeros(Float64, G, l), zeros(Float64, G, q), zeros(Float64, G, q),
              zeros(Float64, G, q, q), ##diagm(0 => ones(Float64, q)),
              0.0, zeros(Float64, G), zeros(Float64, G), zeros(Float64, G, r),
              zeros(Float64, N, q), zeros(Float64, N))
    else
        param = Param(G, 0, zeros(Float64,G,p), randn(G,l), randn(G,q), ones(Float64, G, q),
                      zeros(Float64, G, q, q), ##diagm(0 => ones(Float64, q)),
                      1.0, rand(Uniform(0.5, 1.5), G), rand(G) .+ 0.05, rand(G, r) .- 0.5,
                      randn(N, q), sample(1:G, N))

        for g in 1:G
            param.L[g,:,:] = diagm(0 => ones(Float64, q))
        end
        param
    end
end

function ParamSpline(data::Data, G::Int; set_zero=false)
    N, p = size(data.Xp)

    ## number of hazard predictors
    r = size(data.Xe,2)

    ## number of fixed effects in LMM model
    l = size(data.Xl,3)

    ## number of random effects in LMM model
    q = size(data.Z,3)

    k = size(data.haz_basis,2)

    if set_zero
        ParamSpline(G, 0, zeros(Float64, G, p), zeros(Float64, G, l), zeros(Float64, G, q), zeros(Float64, G, q),
              zeros(Float64, G, q, q), ##diagm(0 => ones(Float64, q)),
              0.0, zeros(Float64, G, k), zeros(Float64, G, r),
              zeros(Float64, N, q), zeros(Float64, N))
    else
        param = ParamSpline(G, 0, zeros(Float64, G, p), randn(G, l), randn(G, q), ones(Float64, G, q),
                            zeros(Float64, G, q, q), ##diagm(0 => ones(Float64, q)),
                            1.0, rand(G, k) .- 1.0, rand(G, r) .- 0.5,
                            randn(N, q), sample(1:G, N))
        for g in 1:G
            param.L[g,:,:] = diagm(0 => ones(Float64, q))
        end
        param
    end
end

function slice_param(par::Tpar, sel) where Tpar <: AbstractParam
    new_par = deepcopy(par)
    new_par.b = par.b[sel,:]
    new_par.cind = par.cind[sel]
    new_par
end


text:
# Prototype for test_parser.jl

using Revise
using TerminalPager
using Test
using OilData
using Dates
using ProgressMeter

# May need to run runtests.jl first to load required functions
# Test, learn, understand here then move to test_parser.jl

# function add_day(start::Dates.DateTime, durations::Vector)
#     ret = fill(start, size(durations))
#     #ret = add_day

#     return d
# end

function learn_add_day()

end

function learn_rsm()
    nt = run_read_rsm()
    df_body = nt.body
    df_meta = nt.meta

    colname = find_column_name(df_meta, "WOPR", "ADA-1762", "")
    @test "WOPR_3" == colname
    days = df_body[!,"TIME"]
    @test days isa Vector{Float64}
    wopr = df_body[!, colname]
    @test wopr isa Vector{Float64}

    # Create vector of datetime from start date and days
    path = joinpath(data_dir(), "test_start_date.PRT")
    start_datetime = find_prt_start_date(path)
    timevec = add_day.(start_datetime, days)
end

function learn_schedule()
    path = joinpath(data_dir(), "test_start_date.PRT")
    start_datetime = find_prt_start_date(path)
    @test DateTime(2020, 1, 18) == start_datetime

    path = joinpath(data_dir(), "test_schedule.ixf")
    end_datetime = find_schedule_end_date(path)
    @test DateTime(2021,2,4) == end_datetime

    #@show delta = end_datetime - start_datetime
    day_count = convert(Day, end_datetime - start_datetime)
    @test Day(383) == day_count
    @test 383 == day_count.value
    @test day_count.value isa Int64

    day_count = convert(Day, DateTime(2021,2,5) - DateTime(2021,2,4))
    @test 1 == day_count.value
end

"""

Keep only `====` sections

```
======================================================================================================================================+
SECTION  The simulation has reached 05-Apr-2020 78 d.  E 3708s = 1h01m48s | C [3691s,3741s] | M [1.4G,1.6G,15.0G] | P [2.2G,2.8G,23.5G]
======================================================================================================================================-
```
"""
function sanitize_prt()
    path = joinpath(data_dir(), "Pad-6A_15_35.PRT")
    open(path, "r") do in
        path = joinpath(data_dir(), "test_current_date.PRT")
        open(path, "w") do out
            for line in eachline(in)
                if !startswith(line, "=============")
                    continue
                end

                println(out, line)
                println(out, readline(in))
                println(out, readline(in))
                println(out)
                #println(out)
            end
        end
    end
end

function learn_progress_meter()
    path = joinpath(data_dir(), "test_start_date.PRT")
    start_datetime = find_prt_start_date(path)

    path = joinpath(data_dir(), "test_schedule.ixf")
    end_datetime = find_schedule_end_date(path)

    day_count = convert(Day, end_datetime - start_datetime)
    progress = Progress(day_count.value, dt=0.5, desc="Running...", color=:green)
    done = Threads.Atomic{Bool}(false)
    path = joinpath(data_dir(), "progress.prt")
    rm(path, force=true)

    writefile = function()
        in  = open(joinpath(data_dir(), "test_current_date.PRT"), "r")
        out = open(path, "w")
        try
            for line in eachline(in)
                sleep(1/1600) # 1600 lines
                for i in 1:100
                     # print gargage to force write to disk
                    println(out, "$i .....................................................................")
                end
                println(out, line)
                #flush(out)
            end
        finally
            close(in)
            close(out)
        end

        Threads.atomic_or!(done, true)
    end

    readfile = function()
        while !isfile(path)
            sleep(0.1)
        end

        pos = 0
        while !done[]
            open(path, "r") do io
                try
                    datetime = find_prt_current_date(io, pos)
                    pos = position(io)
                    days = convert(Day, datetime - start_datetime)
                    update!(progress, days.value)
                catch e
                    if e isa ErrorException
                        # ignore
                    end
                end
            end
            sleep(0.5)
        end

        finish!(progress)
    end

    read_task = Threads.@spawn readfile()
    write_task = Threads.@spawn writefile()

    wait(write_task)
    wait(read_task)
end

function learn_read_write()
    done = Threads.Atomic{Bool}(false)

    path = joinpath(data_dir(), "progress.prt")
    rm(path, force=true)

    writefile = function()
        in  = open(joinpath(data_dir(), "test_current_date.PRT"), "r")
        out = open(path, "w")
        try
            for line in eachline(in)
                sleep(10/1600) # 1600 lines
                for i in 1:100
                     # print gargage to force write to disk
                    println(out, "$i .....................................................................")
                end
                println(out, line)
                #flush(out)
            end
        finally
            close(in)
            close(out)
        end

        Threads.atomic_or!(done, true)
    end

    readfile = function()
        while !isfile(path)
            sleep(0.1)
        end

        pos = 0
        while !done[]
            open(path, "r") do io
                skip(io, pos)

                for line in eachline(io)
                    if !startswith(line, "SECTION")
                        continue
                    end
                    if length(line) < 48
                        continue
                    end
                    println(line)
                end

                pos = position(io)
            end
        end
    end

    read_task = Threads.@spawn readfile()
    write_task = Threads.@spawn writefile()

    wait(write_task)
    wait(read_task)
end

#nt = run_read_rsm()
learn_add_day()
learn_rsm()
learn_schedule()
#sanitize_prt()
#learn_read_write()
learn_progress_meter()

nothing

max_stars_repo_name: ykyang/OilData

max_issues_repo_name: ykyang/OilData

max_forks_repo_name: ykyang/OilData


text:
import QuantumOps
import ElectronicStructure
import QiskitQuantumInfo

include(joinpath(pkgdir(QuantumOps), "test", "runtests.jl"))
include(joinpath(pkgdir(ElectronicStructure), "test", "runtests.jl"))
include(joinpath(pkgdir(QiskitQuantumInfo), "test", "runtests.jl"))

max_stars_repo_name: jlapeyre/qiskit-alt

max_issues_repo_name: jlapeyre/qiskit-alt

max_forks_repo_name: Qiskit-Extensions/qiskit-alt


text:
include("../intcode.jl")

tape = set_memory!(Tape(parse_data(readline(open("input.txt")))), 100000)

@async run_program!(tape)
output = []
while true
    ch = take!(tape.outch)
    if ch == CodeFinished()
        break
    end
    push!(output, ch)
end

instr = reshape(output[1:end], 3, div(length(output[1:end]), 3))

println(sum(instr[3, :].==2))

screen = transpose(reshape(instr[3, :], 41, 24))
for i=1:24
    for j=1:41
        val = screen[i,j]
        if val == 1 ch = '|' elseif val == 2 ch = '#' elseif val == 3 ch = '=' elseif val == 4 ch= 'O' else ch = ' ' end
        print(ch)
    end
    println()
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
function xbacktest(x)
  2x
end


text:
@generated function findsteps!(I::Array{Tuple{Float64,Float64},1}, f::Function, xl::Float64, xr::Float64,
            fl::Int, fr::Int, fmin::Int, fmax::Int, args...)
    quote
        if fl == fr || fr <= fmin || fl > fmax
            return
        end

        if fr - fl == 1
            push!(I, (xl,xr))  # found an interval (xl,xr)
            return
        end

        # subdivide the interval
        xmid = (xl + xr)/2
        fmid = f(xmid,args...)::Int

        # recursive search
        findsteps!(I, f, xl, xmid, fl, fmid, fmin, fmax, args...)
        findsteps!(I, f, xmid, xr, fmid, fr, fmin, fmax, args...)
    end
end

doc"""
`findsteps!(I,f,xl,xr,[fmin,fmax];args=())`

like `findsteps` except that the intervals $(x_l,x_r)$ are appended to `I`.
"""
function findsteps!(I::Array{Tuple{Float64,Float64},1}, f::Function,
        xl::Float64, xr::Float64, fmin::Int, fmax::Int; args=())
    findsteps!(I, f, xl, xr, f(xl,args...)::Int, f(xr,args...)::Int, fmin, fmax, args...)
end

function findsteps!(I::Array{Tuple{Float64,Float64},1}, f::Function, xl::Float64, xr::Float64; args=())
    fl = f(xl,args...)::Int
    fr = f(xr,args...)::Int
    findsteps!(I, f, xl, xr, fl, fr, fl, fr, args...)
end

doc"""
`findsteps(f,xl,xr,[fmin,fmax];args=())`

takes a function `f(x,args...)`, assumed to be a
monotonically increasing integer function of `x`.

optionally one can specify the range of `f` to search
with `fmin`, `fmax`

performs a binary search and returns an ordered list of
intervals $[(x_{l,1},x_{r,1}),(x_{l,2},x_{r,2}),\dots]$
within the specificed range `(xl,xr)` which contains exactly one step.

e.g.
```julia
f = (x,c) -> round(Int,x^2+c)
findsteps(f,0.0,2.0,args=1.0)
```
"""
function findsteps(f::Function, xl::Float64, xr::Float64; args=())
    I = Array{Tuple{Float64,Float64}}(0)
    findsteps!(I, f, xl, xr, args=args)
    I
end

function findsteps(f::Function, xl::Float64, xr::Float64, fmin::Int, fmax::Int; args=())
    I = Array{Tuple{Float64,Float64}}(0)
    findsteps!(I, f, xl, xr, fmin, fmax, args=args)
    I
end

# function findsteps!(I::Array{Tuple{Float64,Float64},1}, f::Function, xl::Float64, xr::Float64,
#             fl::Int, fr::Int, fmin::Int, fmax::Int)
#     if fl == fr || fr <= fmin || fl > fmax
#         return
#     end

#     if fr - fl == 1
#         push!(I, (xl,xr))  # found an interval (xl,xr)
#         return
#     end

#     # subdivide the interval
#     xmid = (xl + xr)/2
#     fmid = f(xmid)::Int

#     # recursive search
#     findsteps!(I, f, xl, xmid, fl, fmid, fmin, fmax)
#     findsteps!(I, f, xmid, xr, fmid, fr, fmin, fmax)
# end

# function findsteps!(I::Array{Tuple{Float64,Float64},1}, f::Function,
#         xl::Float64, xr::Float64, fmin::Int, fmax::Int)
#     findsteps!(I, f, xl, xr, f(xl)::Int, f(xr)::Int, fmin, fmax)
# end

# function findsteps!(I::Array{Tuple{Float64,Float64},1}, f::Function, xl::Float64, xr::Float64)
#     fl = f(xl)::Int
#     fr = f(xr)::Int
#     findsteps!(I, f, xl, xr, fl, fr, fl, fr)
# end

# function findsteps(f::Function, xl::Float64, xr::Float64, fmin::Int, fmax::Int)
#     I = Array{Tuple{Float64,Float64}}(0)
#     findsteps!(I, f, xl, xr, fmin, fmax)
#     I
# end

# doc"""
# `findsteps(f,xl,xr)`

# takes a function $f(x)$, assumed to be a
# monotonically increasing integer function of $x$.

# performs a binary search and returns an ordered list of
# intervals $[(x_{l,1},x_{r,1}),(x_{l,2},x_{r,2}),\dots]$
# within the specificed range `(xl,xr)` which contains exactly one step.

# """
# function findsteps(f::Function, xl::Float64, xr::Float64)
#     I = Array{Tuple{Float64,Float64}}(0)
#     findsteps!(I, f, xl, xr)
#     I
# end

max_stars_repo_name: johncsnyder/Element


text:
using ITensors, Test

import ITensors: In, Out, Neither

@testset "QN Index" begin
  @testset "hasqns function" begin
    i = Index(4, "i")
    @test hasqns(i) == false
    j = Index(QN(0) => 1, QN(1) => 1)
    @test hasqns(j) == true
  end

  @testset "Array of QN Constructor" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    @test hasqns(i)
    @test dim(i) == 3
    @test hastags(i, "i")
  end

  @testset "Vararg Constructor" begin
    i = Index(QN(0) => 1, QN(1) => 2; tags="i")
    @test hasqns(i)
    @test dim(i) == 3
    @test hastags(i, "i")
    @test dir(i) == Out
    @test dir(i => 2) == Out

    j = Index(QN(0) => 1, QN(1) => 2; tags="j", dir=In)
    @test hasqns(j)
    @test dim(j) == 3
    @test hastags(j, "j")
    @test dir(j) == In
    @test dir(j => 2) == In
  end

  @testset "flux and qn" begin
    i = dag(Index([QN(0) => 2, QN(1) => 2], "i"))

    @test flux(i => 1) == QN(0)
    @test flux(i => 2) == QN(0)
    @test flux(i => 3) == QN(-1)
    @test flux(i => 4) == QN(-1)
    @test flux(i => Block(1)) == QN(0)
    @test flux(i => Block(2)) == QN(-1)

    @test qn(i => 1) == QN(0)
    @test qn(i => 2) == QN(0)
    @test qn(i => 3) == QN(1)
    @test qn(i => 4) == QN(1)
    @test qn(i => Block(1)) == QN(0)
    @test qn(i => Block(2)) == QN(1)
  end

  @testset "directsum" begin
    i = Index([QN(0) => 1, QN(1) => 2], "i")
    j = Index([QN(2) => 3, QN(3) => 4], "j")
    ij = ITensors.directsum(i, j; tags="test")
    @test dim(ij) == dim(i) + dim(j)
    @test hastags(ij, "test")
    @test flux(ij, Block(1)) == QN(0)
    @test flux(ij, Block(2)) == QN(1)
    @test flux(ij, Block(3)) == QN(2)
    @test flux(ij, Block(4)) == QN(3)
    @test dim(ij, Block(1)) == 1
    @test dim(ij, Block(2)) == 2
    @test dim(ij, Block(3)) == 3
    @test dim(ij, Block(4)) == 4
    @test_throws ErrorException ITensors.directsum(i, dag(j))
  end
end

nothing


text:
using Ramnet.Utils: stack
using Ramnet.Models: train!, predict, predict_response
using Ramnet.Models: MultiDiscriminatorClassifier, BleachingDiscriminator

@testset "MultiDiscriminatorClassifier" begin
    all_active   = ones(Bool, 9)
    all_inactive = zeros(Bool, 9)

    X = stack(all_active, all_inactive)
    y = String["On", "Off"]

    ## Training with single inputs
    model_1 = MultiDiscriminatorClassifier{String}(9, 3)

    train!(model_1, all_active, "On")
    train!(model_1, all_inactive, "Off")

    @test predict(model_1, all_active) == "On"
    @test predict_response(model_1, all_active) == Dict("On" => 3, "Off" => 0)
    @test predict(model_1, all_inactive) == "Off"
    @test predict_response(model_1, all_inactive) == Dict("On" => 0, "Off" => 3)

    ## Training with multiple inputs
    model_2 = MultiDiscriminatorClassifier{String}(9, 3)

    train!(model_2, X, y)

    @test predict(model_2, X) == ["On", "Off"]
    @test predict_response(model_2, X) == Dict("On" => [3,0], "Off" => [0,3])

    ## Dimensions must match when training with multiple inputs
    model_3 = MultiDiscriminatorClassifier{String}(9, 3)

    @test_throws DimensionMismatch train!(model_3, X, ["On"])

    ## Parameter `width` can be infered from data
    @test_throws DomainError MultiDiscriminatorClassifier{string}(0)
    @test_throws DomainError MultiDiscriminatorClassifier{String}(3; seed=-1)

    model_4 = MultiDiscriminatorClassifier{String}(3)

    train!(model_4, X, y)

    ## Classifiers with nodes that support bleaching
    model_b = MultiDiscriminatorClassifier{String,BleachingDiscriminator}(9, 3)

    train!(model_b, X, y)

    @test predict(model_b, X) == ["On", "Off"]
    @test predict_response(model_b, X) == Dict("On" => [3,0], "Off" => [0,3])
    @test predict_response(model_b, X; b=1) == Dict("On" => [0, 0], "Off" => [0, 0])

    train!(model_b, all_active, "On")

    @test predict_response(model_b, X; b=1) ==  Dict("On" => [3, 0], "Off" => [0, 0])

    ## Linear search bleaching
    model_b = MultiDiscriminatorClassifier{String,BleachingDiscriminator}(9, 3)

    train!(model_b, X, y)
    # Simulating ambiguity. "Off" target is also trained with the all_active input
    train!(model_b, stack(all_active, all_active), ["On", "Off"])

    @test predict_response(model_b, all_active; b=:linear) == Dict("On" => 3, "Off" => 0)
    @test predict_response(model_b, all_inactive; b=:linear) == Dict("On" => 0, "Off" => 3)
    # Test with multiple inputs (last to in one go)
    @test predict_response(model_b, X; b=:linear) == Dict("On" => [3,0], "Off" => [0,3])
    @test predict(model_b, X; b=:linear) == ["On", "Off"]

    # Future tests:
    # - predict over model trained on a single pattern
    # - predict over a model that hasn't been trained at all
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
global test_suite_polyhedra

for N in [Float64, Rational{Int}, Float32]
    # random polytopes
    if test_suite_polyhedra
        rand(HPolytope)
    else
        @test_throws AssertionError rand(HPolytope)
    end
    rand(VPolytope)

    # -----
    # H-rep
    # -----

    # constructor from matrix and vector
    A = [N(1) N(2); N(-1) N(1)]
    b = [N(1), N(2)]
    p = HPolytope(A, b)
    c = p.constraints
    @test c isa Vector{LinearConstraint{N}}
    @test c[1].a == N[1, 2] && c[1].b == N(1)
    @test c[2].a == N[-1, 1] && c[2].b == N(2)
    @test_throws AssertionError HPolytope(N[1 0; 0 1], N[1, 1];
                                         check_boundedness=true)

    # convert back to matrix and vector
    A2, b2 = tosimplehrep(p)
    @test A == A2 && b == b2

    # 2D polytope
    p = HPolytope{N}()
    c1 = LinearConstraint(N[2, 2], N(12))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(0))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p, c3)
    addconstraint!(p, c1)
    addconstraint!(p, c4)
    addconstraint!(p, c2)

    # support vector
    d = N[1, 0]
    @test σ(d, p) == N[4, 2]
    d = N[0, 1]
    @test σ(d, p) == N[2, 4]
    d = N[-1, 0]
    @test σ(d, p) == N[-1, 1]
    d = N[0, -1]
    @test σ(d, p) == N[0, 0]

    # support vector of polytope with no constraints
    @test_throws ErrorException σ(N[0], HPolytope{N}())

    # boundedness
    @test isbounded(p) && isbounded(p, false)
    p2 = HPolytope{N}()
    @test isbounded(p2) && !isbounded(p2, false)

    # membership
    @test N[5 / 4, 7 / 4] ∈ p
    @test N[4, 1] ∉ p

    # singleton list (only available with Polyhedra library)
    if test_suite_polyhedra
        @test length(singleton_list(p)) == 4
    else
        @test_throws AssertionError singleton_list(p)
    end

    if test_suite_polyhedra
        # conversion to and from Polyhedra's VRep data structure
        cl = constraints_list(HPolytope(polyhedron(p)))
        @test length(p.constraints) == length(cl)
    end

    # vertices_list of "universal polytope" (strictly speaking: illegal input)
    @test vertices_list(HPolytope{N}()) == Vector{Vector{N}}()

    if test_suite_polyhedra
        # convert hyperrectangle to a HPolytope
        H = Hyperrectangle(N[1, 1], N[2, 2])
        P = convert(HPolytope, H)
        vlist = vertices_list(P)
        @test ispermutation(vlist, [N[3, 3], N[3, -1], N[-1, -1], N[-1, 3]])
        # check boundedness after conversion
        HPolytope(constraints_list(H); check_boundedness=true)

        # isempty
        @test !isempty(p)
        @test !isempty(HPolytope{N}())

        # H-representaion of an empty v-polytope
        @test tohrep(VPolytope{N}()) == EmptySet{N}()
    end

    # remove redundant constraints
    P = HPolytope([HalfSpace(N[1, 0], N(1)),
                   HalfSpace(N[0, 1], N(1)),
                   HalfSpace(N[-1, -0], N(1)),
                   HalfSpace(N[-0, -1], N(1)),
                   HalfSpace(N[1, 0], N(2))]) # redundant

    Pred = remove_redundant_constraints(P)
    @test length(Pred.constraints) == 4
    @test length(P.constraints) == 5

    # test in-place removal of redundancies
    remove_redundant_constraints!(P)
    @test length(P.constraints) == 4

    # translation
    P2 = translate(P, N[1, 2])
    @test P2 isa HPolytope && ispermutation(constraints_list(P2),
        [HalfSpace(N[1, 0], N(2)), HalfSpace(N[0, 1], N(3)),
         HalfSpace(N[-1, -0], N(0)), HalfSpace(N[-0, -1], N(-1))])

    # subset
    H = BallInf(N[0, 0], N(1))
    P = convert(HPolytope, H)
    @test BallInf(N[0, 0], N(1)) ⊆ P
    @test !(BallInf(N[0, 0], N(1.01)) ⊆ P)

    # =====================
    # Concrete linear map
    # =====================
    LM = linear_map(N[2 3; 1 2], P) # invertible matrix
    @test LM isa HPolytope
    if test_suite_polyhedra
        LM = linear_map(N[2 3; 0 0], P)  # non-invertible matrix
        @test LM isa VPolygon
    end

    M = N[2 1; 0 1]
    L1 = linear_map(M, P, use_inv=true)  # calculates inv(M) explicitly
    L2 = linear_map(M, P, use_inv=false) # uses transpose(M) \ c.a for each constraint c of P
    L3 = linear_map(M, P, cond_tol=1e3)  # set a custom tolerance for the condition number (invertibility check)
    # needs M * an_element(Li) to be stable for rational, see #1105
    p = convert(Vector{N}, an_element(P))
    @assert p ∈ P
    @test all([M * p ∈ Li for Li in [L1, L2, L3]])

    # do not check for invertibility => use the vertices
    L4 = linear_map(M, P, check_invertibility=false)
    @test L4 isa VPolygon

    # linear map for mixed types
    M = [2 1; 0 1] # Int's
    LM = linear_map(M, P)
    @test LM isa HPolytope{N}

    # -----
    # V-rep
    # -----

    # constructor from a VPolygon
    polygon = VPolygon([N[0, 0], N[1, 0], N[0, 1]])
    p = VPolytope(polygon)
    @test vertices_list(polygon) == vertices_list(p)

    # dim
    @test dim(p) == 2

    # support vector
    d = N[1, 0]
    @test σ(d, p) == N[1, 0]

    # boundedness
    @test isbounded(p)

    # isempty
    @test !isempty(p)
    @test isempty(VPolytope{N}())

    # vertices_list function
    @test vertices_list(p) == p.vertices

    if test_suite_polyhedra
        V = VPolytope(polyhedron(p))

        # conversion to and from Polyhedra's VRep data structure
        vl = vertices_list(V)
        @test length(p.vertices) == length(vl) && vl ⊆ p.vertices

        # list of constraints of a VPolytope; calculates tohrep
        @test ispermutation(constraints_list(V), constraints_list(tohrep(V)))

        # convert empty VPolytope to a polyhedron
        Vempty = VPolytope()
        @test_throws ErrorException polyhedron(Vempty) # needs to pass the (relative) dim
        Pe = polyhedron(Vempty, relative_dimension=2)
    end

    # membership
    @test N[.49, .49] ∈ p
    @test N[.51, .51] ∉ p

    # translation
    @test translate(p, N[1, 2]) == VPolytope([N[1, 2], N[2, 2], N[1, 3]])

    # copy (see #1002)
    p, q = [N(1)], [N(2)]
    P = VPolytope([p, q])
    Pcopy = copy(P)
    p[1] = N(5)
    # test that Pcopy is independent of P ( = deepcopy)
    @test Pcopy.vertices[1] == [N(1)]
end

# default Float64 constructors
@test HPolytope() isa HPolytope{Float64}
@test VPolytope() isa VPolytope{Float64}

# Polyhedra tests that only work with Float64
if test_suite_polyhedra
    for N in [Float64]
        # -----
        # H-rep
        # -----
        # support function/vector
        d = N[1, 0]
        p_unbounded = HPolytope([LinearConstraint(N[-1, 0], N(0))])
        @test_throws ErrorException σ(d, p_unbounded)
        @test_throws ErrorException ρ(d, p_unbounded)
        p_infeasible = HPolytope([LinearConstraint(N[1], N(0)),
                                  LinearConstraint(N[-1], N(-1))])
        @test_throws ErrorException σ(N[1], p_infeasible)
        @test_throws ErrorException ρ(N[1], p_infeasible)

        # empty intersection
        A = [N(0) N(-1); N(-1) N(0); N(1) N(1)]
        b = N[0, 0, 1]
        p1 = HPolytope(A, b)
        A = [N(0) N(1); N(1) N(0); N(-1) N(-1)]
        b = N[2, 2, -2]
        p2 = HPolytope(A, b)
        @test intersection(p1, p2) isa EmptySet{N}
        @test intersection(p1, p2; backend=Polyhedra.default_library(2, N)) isa EmptySet{N}
        @test intersection(p1, p2; backend=CDDLib.Library()) isa EmptySet{N}

        # intersection with half-space
        hs = HalfSpace(N[2, -2], N(-1))
        c1 = constraints_list(intersection(p1, hs))
        c2 = constraints_list(intersection(hs, p1))
        @test length(c1) == 3 && ispermutation(c1, c2)

        # Cartesian product
        A = [N(1) N(-1)]'
        b = N[1, 0]
        p1 = HPolytope(A, b)
        p2 = copy(p1)
        cp = cartesian_product(p1, p2)
        cl = constraints_list(cp)
        @test length(cl) == 4

        # vertices_list
        vl = vertices_list(p1)
        @test ispermutation(vl, [N[0], N[1]])

        # checking for emptiness
        P = HPolytope([LinearConstraint(N[1, 0], N(0))])      # x <= 0
        @test !isempty(P)
        addconstraint!(P, LinearConstraint(N[-1, 0], N(-1)))  # x >= 1
        @test isempty(P)

        # checking for empty intersection (also test symmetric calls)
        P = convert(HPolytope, BallInf(zeros(N, 2), N(1)))
        Q = convert(HPolytope, BallInf(ones(N, 2), N(1)))
        R = convert(HPolytope, BallInf(3*ones(N, 2), N(1)))
        res, w = isdisjoint(P, Q, true)
        @test !isdisjoint(P, Q) && !res && w ∈ P && w ∈ Q
        @test !isdisjoint(P, Q; algorithm="sufficient")
        res, w = isdisjoint(Q, P, true)
        @test !isdisjoint(Q, P) && !res && w ∈ P && w ∈ Q
        res, w = isdisjoint(Q, R, true)
        @test !isdisjoint(Q, R) && !res && w ∈ Q && w ∈ R
        res, w = isdisjoint(R, Q, true)
        @test !isdisjoint(R, Q) && !res && w ∈ Q && w ∈ R
        res, w = isdisjoint(P, R, true)
        @test isdisjoint(P, R) && res && w == N[]
        res, w = isdisjoint(R, P, true)
        @test isdisjoint(R, P) && res && w == N[]
        res, w = isdisjoint(P, R, true; algorithm="sufficient")
        @test isdisjoint(P, R; algorithm="sufficient") && res && w == N[]

        # test that one can pass a sparse vector as the direction (see #1011)
        P = HPolytope([HalfSpace(N[1, 0], N(1)),
                       HalfSpace(N[0, 1], N(1)),
                       HalfSpace(N[-1, -1], N(-1))])
        @test an_element(P) ∈ P

        # tovrep from HPolytope
        A = N[0 1; 1 0; -1 -1]
        b = N[0.25, 0.25, -0]
        P = HPolytope(A, b)
        @test tovrep(P) isa VPolytope{N}
        @test tohrep(P) isa HPolytope{N}  # no-op

        # -----
        # V-rep
        # -----

        # intersection
        p1 = VPolytope(vertices_list(BallInf(N[0, 0], N(1))))
        p2 = VPolytope(vertices_list(BallInf(N[2, 2], N(1))))
        cap = intersection(p1, p2)
        @test vertices_list(cap) ≈ [N[1, 1]]
        # other polytopic sets
        p3 = VPolygon(vertices_list(p2))
        cap = intersection(p1, p3)
        @test vertices_list(cap) ≈ [N[1, 1]]
        p4 = BallInf(N[2, 2], N(1))
        cap = intersection(p1, p4)
        vlist = vertices_list(cap)  # contains duplicates (see #1405)
        @test all(v -> v == N[1, 1], vlist)

        # convex hull
        v1 = N[1, 0]
        v2 = N[0, 1]
        v3 = N[-1, 0]
        v4 = N[0, -1]
        v5 = N[0, 0]
        p1 = VPolytope([v1, v2, v5])
        p2 = VPolytope([v3, v4])
        ch = convex_hull(p1, p2)
        vl = vertices_list(ch)
        @test ispermutation(vl, [v1, v2, v3, v4])

        # Cartesian product
        p1 = VPolytope([N[0, 0], N[1, 1]])
        p2 = VPolytope([N[2]])
        cp = cartesian_product(p1, p2)
        vl = vertices_list(cp)
        @test ispermutation(vl, [N[0, 0, 2], N[1, 1, 2]])

        # tohrep from VPolytope
        P = VPolytope([v1, v2, v3, v4, v5])
        @test tohrep(P) isa HPolytope{N}
        @test tovrep(P) isa VPolytope{N}  # no-op

        # subset (see #974)
        H = BallInf(N[0, 0], N(1))
        P = convert(VPolytope, H)
        @test BallInf(N[0, 0], N(1)) ⊆ P
        @test !(BallInf(N[0, 0], N(1.01)) ⊆ P)

        # concrete minkowski sum
        B = convert(VPolytope, BallInf(N[0, 0, 0], N(1)))
        X = minkowski_sum(B, B)
        twoB = 2.0*B
        @test X ⊆ twoB && twoB ⊆ X

        # same but specifying a custom polyhedral computations backend (CDDLib)
        X = minkowski_sum(B, B, backend=CDDLib.Library())
        twoB = 2.0*B
        @test X ⊆ twoB && twoB ⊆ X

        P1 = VPolytope([N[0, 0, 0], N[0, 1, 0]])
        P2 = VPolytope([N[0, 0, 0], N[1, 0, 0]])
        Q = minkowski_sum(P1, P2)
        @test ispermutation(vertices_list(Q), [N[0, 0, 0], N[0, 1, 0], N[1, 0, 0], N[1, 1, 0]])

        # -----------------
        # mixed H-rep/V-rep
        # -----------------

        # intersection
        p1 = convert(HPolytope, BallInf(N[0, 0], N(1)))
        p2 = convert(VPolytope, BallInf(N[1, 1], N(1)))
        cap1 = intersection(p1, p2)
        cap2 = intersection(p2, p1)
        @test ispermutation([round.(v) for v in vertices_list(cap1)],
                            [N[1, 1], N[0, 1], N[1, 0], N[0, 0]])
        @test ispermutation([round.(v) for v in vertices_list(cap1)],
                            [round.(v) for v in vertices_list(cap2)])
    end
end


text:
module TestRecursionLimitError

using SplittablesBase
using SplittablesBase.Testing: RecursionLimitError
using Test

struct BrokenSplittable end

Base.iterate(::BrokenSplittable, n = 1) = n > 3 ? nothing : (n, n + 1)
Base.eltype(::Type{BrokenSplittable}) = Int
Base.length(::BrokenSplittable) = 3

# Invalid definition:
SplittablesBase.halve(::BrokenSplittable) = (BrokenSplittable(), ())

@testset begin
    # Make sure it's OK as an iterable:
    @test collect(BrokenSplittable()) == 1:3

    err = try
        SplittablesBase.Testing.recursive_vcat(BrokenSplittable())
        nothing
    catch err
        err
    end
    @test err isa RecursionLimitError
    msg = sprint(showerror, err)
    @debug "Testing an error from `recursive_vcat(BrokenSplittable())`" Text(msg)
    @test occursin("RecursionLimitError", msg)
    @test occursin("BrokenSplittable", msg)
end

@testset "smoke tests" begin
    @test sprint(showerror, RecursionLimitError(identity, ())) isa AbstractString
    @test sprint(showerror, RecursionLimitError(identity, (1,))) isa AbstractString
    @test sprint(showerror, RecursionLimitError(identity, (), (a = 1,))) isa AbstractString
    @test sprint(showerror, RecursionLimitError(identity, (1,), (a = 1,))) isa
          AbstractString
end

end  # module


text:
module Penopt
using LinearAlgebra: norm
using Optim

export penopt

function penopt(f, c, x_0, Λ=vcat([0], [2 .^k for k in 0:10]); p=2)
    x = x_0
    fx=0.0
    for λ ∈ Λ
        xopt = optimize(arg -> f(arg) + λ*posnorm(c(arg)).^p, x
            , NelderMead())
            # , BFGS())
            # , ConjugateGradient(), Optim.Options())
        x=xopt.minimizer
        fx=xopt.minimum
        # println("f(",x,")=", fx, ", λ=", λ);
    end
    x,fx
end

function posnorm(cx)
    cxp = filter(arg -> arg > 0, cx)
    norm(cxp)
end

end


text:
# 2D FWI on Overthrust model with L-BFGS using NLopt library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: December 2017
#

using Pkg; Pkg.activate("JUDI")
using Statistics, Random, LinearAlgebra
using JUDI.TimeModeling, HDF5, NLopt, SeisIO

# Load starting model
n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./ model0.m)
vmin = ones(Float32, model0.n) .* 1.3f0
vmax = ones(Float32, model0.n) .* 6.5f0

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data
block = segy_read("../../data/overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0)    # 8 Hz wavelet
q = judiVector(src_geometry, wavelet)

############################### FWI ###########################################

# optimization parameters
batchsize = 16
count = 0

# NLopt objective function
println("No.  ", "fval         ", "norm(gradient)")
function f!(x,grad)

    # Update model
    model0.m = convert(Array{Float32, 2}, reshape(x, model0.n))

    # Seclect batch and calculate gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, gradient = fwi_objective(model0, q[i], d_obs[i])

    # Reset gradient in water column to zero
    gradient = reshape(gradient, model0.n); gradient[:, 1:21] .= 0f0
    grad[1:end] = vec(gradient)

    global count; count += 1
    println(count, "    ", fval, "    ", norm(grad))
    return convert(Float64, fval)
end

# Optimization parameters
opt = Opt(:LD_LBFGS, prod(model0.n))
lower_bounds!(opt, mmin); upper_bounds!(opt, mmax)
min_objective!(opt, f!)
maxeval!(opt, 10)
(minf, minx, ret) = optimize(opt, vec(model0.m))


text:
mul(::Any, x, y) = x * y
mul(C, A, B, α, β) = A * B * α + C * β


text:
using Mimi
using Random

page_years = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
page_year_0 = 2015

function getpageindexfromyear(year)
    i = findfirst(isequal(year), page_years)
    if i == 0
        error("Invalid PAGE year: $year.")
    end
    return i
end

function getperiodlength(year)      # same calculations made for yagg_periodspan in the model
    i = getpageindexfromyear(year)

    if year==page_years[1]
        start_year = page_year_0
    else
        start_year = page_years[i - 1]
    end

    if year == page_years[end]
        last_year = page_years[end]
    else
        last_year = page_years[i + 1]
    end

    return (last_year - start_year) / 2
end

@defcomp PAGE_marginal_emissions begin
    er_CO2emissionsgrowth = Variable(index=[time,region], unit = "%")
    marginal_emissions_growth = Parameter(index=[time,region], unit = "%", default = zeros(10,8))
    function run_timestep(p, v, d, t)
        if is_first(t)
            v.er_CO2emissionsgrowth[:, :] = p.marginal_emissions_growth[:, :]
        end
    end
end

"""
compute_scc(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing)
Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiPAGE2009 model.
If no model is provided, the default model from MimiPAGE2009.get_model() is used.
Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model.
If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
The emission pulse can be set manually and otherwise defaults to 100000.
"""

# changed default pulse_size from 100,000 to 10,000 Mtonne/year :
function compute_scc(m::Model = getpage(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing,
                    pulse_size::Union{Float64, Nothing} = 2000.)
    return compute_scc_mm(m; year=year, eta=eta, prtp=prtp, pulse_size=pulse_size)[:scc]
end

"""
compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing)
Returns a NamedTuple (scc=scc, scc_disaggregated=scc_disaggregated, mm=mm) of the social cost of carbon, its disaggregation by region and year, and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiPAGE2009 model.
If no model is provided, the default model from MimiPAGE2009.get_model() is used.
Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model.
If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
The emission pulse can be set manually and otherwise defaults to 100000.
"""
function compute_scc_mm(m::Model = getpage(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing,
                            pulse_size::Union{Float64, Nothing} = 2000.)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    eta == nothing ? nothing : update_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : update_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size = pulse_size)   # Returns a marginal model that has already been run
    if use_annual
        scc = mm[:EquityWeighting, :td_totaldiscountedimpacts_ann] # CHANGED: made to new annual variable
        scc_disaggregated = mm[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated_ann] # CHANGED: made to new annual variable
    else
        scc = mm[:EquityWeighting, :td_totaldiscountedimpacts]
        scc_disaggregated = mm[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated]
    end

    return (scc = scc, scc_disaggregated = scc_disaggregated, mm = mm)
end

"""
get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing)
Returns a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
If no Model m is provided, the default model from MimiPAGE2009.get_model() is used as the base model.
Note that the returned MarginalModel has already been run.
"""
function get_marginal_model(m::Model = getpage(); year::Union{Int, Nothing} = nothing, pulse_size::Union{Float64, Nothing} = nothing)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $page_years.") : nothing
    pulse_size === nothing ? error("Must specify an emission pulse size. Try `get_marginal_model(m, year=2020, pulse_size = 10000.)`.") : nothing

    mm = create_marginal_model(m, pulse_size)

    add_comp!(mm.marginal, PAGE_marginal_emissions, :marginal_emissions; before = :co2emissions)
    connect_param!(mm.marginal, :co2emissions=>:er_CO2emissionsgrowth, :marginal_emissions=>:er_CO2emissionsgrowth)
    connect_param!(mm.marginal, :AbatementCostsCO2=>:er_emissionsgrowth, :marginal_emissions=>:er_CO2emissionsgrowth)

    i = getpageindexfromyear(year)

    randomdraw = rand(1:1000000000000) # to prevent different variability effects between the marginal models

    # Base model
    Random.seed!(randomdraw);
    run(mm.base)
    base_glob0_emissions = mm.base[:CO2Cycle, :e0_globalCO2emissions]
    er_co2_a = mm.base[:co2emissions, :er_CO2emissionsgrowth][i, :]
    e_co2_g = mm.base[:co2emissions, :e_globalCO2emissions]

    # Calculate pulse
    ER_SCC = 100 * -1 * pulse_size / (base_glob0_emissions * getperiodlength(year))
    pulse = er_co2_a - ER_SCC * (er_co2_a/100) * (base_glob0_emissions / e_co2_g[i])
    marginal_emissions_growth = copy(mm.base[:co2emissions, :er_CO2emissionsgrowth])
    marginal_emissions_growth[i, :] = pulse
    # println(string(copy(mm.base[i,:][:co2emissions, :er_CO2emissionsgrowth])))
    # println(string(pulse))

    # Marginal emissions model
    update_param!(mm.marginal, :marginal_emissions_growth, marginal_emissions_growth)
    Random.seed!(randomdraw);
    run(mm.marginal)
    # print(string(" ... it runs with pulse size ", pulse_size," ... "))

    # print(string("Original global CO2 emissions are: ", mm.base[:co2emissions, :e_globalCO2emissions][i, :]))
    # print(string("Pulsed global CO2 emissions are: ", mm.marginal[:co2emissions, :e_globalCO2emissions][i, :]))

    return mm
end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
module Utils
using OffsetArrays

export @something_nothing
export @aoc_str
export split_blocks
export slurp
export per_split
export per_line
export per_split_parse
export per_line_parse
export find_input
export zerobased
export curry
export curry2nd
export currylast
export assertequal
export boundaryconditions

using Chain
using Underscores

function _something_nothing_impl(thing)
    quote
        local evaluated = $(esc(thing))
        if isa(evaluated, Some)
            evaluated.value
        else
            evaluated
        end
    end
end

function _something_nothing_impl(thing, rest...)
    quote
        local evalued = $(esc(thing))
        if isnothing(evalued)
            $(_something_nothing_impl(rest...))
        else
            something(evalued)
        end
    end
end

macro something_nothing(things...)
    _something_nothing_impl(things...)
end

function find_input(julia_file, lookfor=nothing)
    if isnothing(lookfor)
        num = escape_string(match(r"(\d+)\.jl$", julia_file)[1])
        lookfor = Regex("\\Q$num\\E[^\\/]*\\.txt\$")
    end
    dir, _, files = @chain julia_file begin
        dirname
        walkdir
        first
    end
    @chain files begin
        filter(fname -> occursin(lookfor, fname), _)
        Iterators.only
        joinpath(dir, _)
    end
end

slurp(filename::AbstractString) = read(filename, String)

function split_blocks(f, buffer::Base.IO)
    re = String[]
    block = String[]
    for line in eachline(buffer)
        if f(line)
            push!(re, join(block, '\n'))
            block = String[]
        else
            push!(block, line)
        end
    end
    push!(re, join(block, '\n'))
    re
end

split_blocks(str::Base.String, buffer::Base.IO) = @_ split_blocks(==(str), buffer)

split_blocks(buffer::Base.IO) = split_blocks(isempty, buffer)

per_split(str::AbstractString, sep=isspace, keepempty=true) = @_ str |> split(__, sep, keepempty=keepempty)

per_split(f::Function, str::AbstractString, sep=isspace, keepempty=true) = @_ str |> split(__, sep, keepempty=keepempty) |> map(f, __)

const NEWLINE = r"\r?\n"

per_line(str::AbstractString, keepempty=true) = per_split(str, NEWLINE, keepempty)

per_line(f::Function, str::AbstractString, keepempty=true) = per_split(f, str, NEWLINE, keepempty)

function per_split_parse(str::AbstractString, sep, ::Type{T} = Int) where {T <: Number}
    per_split(str, sep, false) do line
        parse(T, line)
    end
end

function per_line_parse(str::AbstractString, ::Type{T} = Int) where {T <: Number}
    per_line(str, false) do line
        parse(T, line)
    end
end

zerobased(array::AbstractArray{T,1}) where {T} = OffsetArray(array, OffsetArrays.Origin(0))
zerobased(array::AbstractArray{T,2}) where {T} = OffsetArray(array, OffsetArrays.Origin(0, 0))
zerobased(array::AbstractArray{T,3}) where {T} = OffsetArray(array, OffsetArrays.Origin(0, 0, 0))

curry(f::Function, x)::Function = Base.Fix1(f, x)
curry2nd(f::Function, x)::Function = Base.Fix2(f, x)
currylast(f::Function, x)::Function = (xs...) -> f(xs..., x)

function assertequal(result, compareto=nothing, _ = Union{})
    if !isnothing(compareto)
        if result == compareto
            println("Correct   result: $(result)")
        else
            println("INCORRECT result: $(result) should be $(compareto)")
        end
    else
        println("Unchecked result: $(result)")
    end
end

function assertequal(result::AbstractString)
    assertequal("\"$(escape_string(result))\"", nothing)
end

function assertequal(result::AbstractString, compareto::AbstractString)
    assertequal("\"$(escape_string(result))\"", "\"$(escape_string(compareto))\"", Union{})
end

function boundaryconditions(vector::AbstractArray{T,1}, fillvalue::Union{Nothing,T}=nothing) where T
    if fillvalue === nothing
        [vector[begin]; vector; vector[end]]
    else
        [fillvalue; vector; fillvalue]
    end
end

function boundaryconditions(matrix::AbstractArray{T,2}, fillvalue::Union{Nothing,T}=nothing) where T
    sub = mapslices(curry2nd(boundaryconditions, fillvalue), matrix, dims=2)
    if fillvalue === nothing
        head = permutedims(sub[firstindex(sub, 1),:])
        tail = permutedims(sub[lastindex(sub, 1),:])
        cat(head, sub, tail; dims=1)
    else
        filler = permutedims(fill(fillvalue, size(sub, 2)))
        cat(filler, sub, filler; dims=1)
    end
end

boundaryconditions(matrix::AbstractArray{T,N}, fillvalue::U) where {N,T,U} = boundaryconditions(matrix, convert(T, fillvalue))

end

max_stars_repo_name: CmdQ/AoC


text:
module Diffusion

using OceanTurb

# Just one field: "c"
@solution c

Base.@kwdef struct Parameters{T} <: AbstractParameters
    K :: T = 0.0
    W :: T = 0.0
    μ :: T = 0.0
end

struct Model{P, TS, G, T, S} <: AbstractModel{TS, G, T}
    clock       :: Clock{T}
    grid        :: G
    timestepper :: TS
    solution    :: S
    bcs         :: BoundaryConditions
    parameters  :: P
end

function Model(; N=10, L=1.0, K=0.1, W=0.0, μ=0.0,
          grid = UniformGrid(N, L),
    parameters = Parameters(K, W, μ),
       stepper = :ForwardEuler,
           bcs = BoundaryConditions(ZeroFluxBoundaryConditions())
    )

    solution = Solution(CellField(grid))

      K = (c=Kc,)
      W = (c=Wc,)
      L = (c=Lc,)

    eqn = Equation(K=K, M=W, L=L)
    lhs = build_lhs(solution)

    timestepper = Timestepper(stepper, eqn, solution, lhs)

    return Model(Clock(), grid, timestepper, solution, bcs, parameters)
end

@inline Kc(m, i) = m.parameters.K
@inline Wc(m, i) = m.parameters.W
@inline Lc(m, i) = m.parameters.μ

end # module


text:
#========================================================================================#
#	Ingolstadt
#
# This include file defines the structure of Activities and Sessions.
#
# Author: Niall Palfreyman, 7/12/2021
#========================================================================================#
#-----------------------------------------------------------------------------------------
# Concrete types:

#-----------------------------------------------------------------------------------------
"""
	Activity

A single activity to be performed and responded to by the learner.

An Activity is the basic unit of learning in Ingolstadt. It contains ...

# Fields
* `affordance`	: Text affording an activity to the learner.
* `hint`		: Text suggesting to the learner how to respond to the affordance.
* `success`		: Boolean function that defines the criteria for a successful learner response.

# Notes
* None

# Examples
```julia
julia> act = Activity( "What is 2+3?", "Add the numbers together", x -> x==5)
```
"""
struct Activity
	affordance :: String				# Affordance text
	hint :: String						# A hint to the learner, if she requires one
	success :: Function					# The criterion for a successful response
end

#-----------------------------------------------------------------------------------------
"""
	Session

Current state of a laboratory of learning activities.

A Session maintains the current state of a learning laboratory experience. It contains ...

# Fields
* `lnr_file :: String`				: Registration file of the current learner
* `lab_path :: String`				: Path of the file containing the current laboratory
* `lab_num :: Int`					: Number of the current laboratory
* `is_pluto :: Bool`				: Is the current laboratory a Pluto workshop?
* `current_act :: Int`				: Number of the current activity in the laboratory
* `activities :: Vector{Activity}`	: Complete list of activities in this laboratory

# Notes
* None

# Examples
```julia
julia> sesh = Session()
```
"""
mutable struct Session
	lnr_file :: String					# Learner registration file
	lab_path :: String					# Path of current labfile
	lab_num :: Int						# Number of current lab
	is_pluto :: Bool					# Is this a Pluto session?
	current_act :: Int					# Number of current activity
	activities :: Vector{Activity}		# Set of activities in this session

	# One and only constructor of Sessions:
	Session() = new( "", "", 1, false, 1, Vector{Activity}())
end

#-----------------------------------------------------------------------------------------
# Methods:

#-----------------------------------------------------------------------------------------
"""
pose( act::Activity)

Present the affordance of the given Activity to the learner.
"""
function pose( act::Activity)
	print( act.affordance)
end

#-----------------------------------------------------------------------------------------
"""
evaluate( act::Activity, response)

Apply the success criterion of the activity to the given response from a learner.
"""
function evaluate( act::Activity, response)
	success = true
	try
		success = act.success(response)
	catch e
		# Response is badly wrong:
		success = false
	end

	if success
		# Success criterion is fulfilled:
		congratulate()
	else
		# Success criterion is not fulfilled:
		commiserate( act)
	end

	success
end

#-----------------------------------------------------------------------------------------
"""
congratulate()

Provide uplifting feedback to a successful learner response.
"""
function congratulate()
	congrats = [
		"Well done - great work!",
		"Great job - well done!",
		"You're doing great - keep it up!",
		"Good job!",
		"Good work!",
		"Nicely done!",
		"Yay! Well done!",
	]

	println()
	print( rand(congrats), " :) ")
end

#-----------------------------------------------------------------------------------------
"""
commiserate( act::Activity)

Provide supportive feedback to an unsuccessful learner response.
"""
function commiserate( act::Activity)
	commiseration = [
		"Not quite right, I'm afraid",
		"No, sorry - not quite right",
		"Oh, hard luck!",
		"No, not quite",
		"Sorry - that's not it",
		"Think again",
		"Nope, sorry",
	]

	println()
	print( rand(commiseration), " - you could ask for a hint(). ")
end

#-----------------------------------------------------------------------------------------
"""
labfile( path::String, labnum::Int)

Construct a valid labfile name from the given path and laboratory number.
"""
function labfile( path, labnum)
	joinpath( path, "INLab" * lpad("$labnum",3,'0')[1:3] * ".jl")
end


text:
mutable struct Agent
    opinion::Float64
    inclin_interact::Float64
    perceiv_publ_opinion::Float64
    check_regularity::Float64
    active::Bool
    inactive_ticks::Int16
    feed::Array{Tweet, 1}
    liked_Tweets::Array{Tweet, 1}
    retweeted_Tweets::Array{Tweet, 1}
    function Agent(opinion, inclin_interact, check_regularity)
        # check if opinion value is valid
        if opinion < -1 || opinion > 1
            error("invalid opinion value")
        end
        # check if value for inclination to interact is valid
        if inclin_interact < 0
            error("invalid value for inclination to interact")
        end
        new(
            opinion,
            inclin_interact,
            opinion,
            check_regularity,
            true,
            0,
            Array{Tweet, 1}(undef, 0),
            Array{Tweet, 1}(undef, 0),
            Array{Tweet, 1}(undef, 0)
        )
    end
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : false
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 64 - 1
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 64 - 1
const nz          = haskey(ENV, "NZ"         ) ? parse(Int , ENV["NZ"]         ) : 64 - 1
###
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using ImplicitGlobalGrid, Plots, Printf, Statistics, LinearAlgebra, MAT
import MPI
# Global reductions
mean_g(A)    = (mean_l = mean(A); MPI.Allreduce(mean_l, MPI.SUM, MPI.COMM_WORLD)/MPI.Comm_size(MPI.COMM_WORLD))
norm_g(A)    = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))
maximum_g(A) = (max_l  = maximum(A); MPI.Allreduce(max_l,  MPI.MAX, MPI.COMM_WORLD))
minimum_g(A) = (min_l  = minimum(A); MPI.Allreduce(min_l,  MPI.MIN, MPI.COMM_WORLD))
# CPU functions
@views av_xza(A) = (A[1:end-1,:,1:end-1] .+ A[1:end-1,:,2:end] .+ A[2:end,:,1:end-1] .+ A[2:end,:,2:end]).*0.25
@views av_zi(A)  = (A[2:end-1,2:end-1,2:end-2] .+ A[2:end-1,2:end-1,3:end-1]).*0.5
@views av_za(A)  = (A[:,:,1:end-1] .+ A[:,:,2:end]).*0.5
@views inn(A)    =  A[2:end-1,2:end-1,2:end-1]

@parallel function smooth!(A2::Data.Array, A::Data.Array, fact::Data.Number)
    @inn(A2)   = @inn(A) + 1.0/6.1/fact*(@d2_xi(A) + @d2_yi(A) + @d2_zi(A))
    return
end

@parallel function compute_maxloc!(Musτ2::Data.Array, Musτ::Data.Array)
    @inn(Musτ2) = @maxloc(Musτ)
    return
end

macro Mu_eff() esc(:(1.0/(1.0/@all(Musτ) + 1.0/(G*dt)))) end
@parallel function compute_iter_params!(dτ_Rho::Data.Array, Gdτ::Data.Array, Musτ::Data.Array, Vpdτ::Data.Number, G::Data.Number, dt::Data.Number, Re::Data.Number, r::Data.Number, max_lxyz::Data.Number)
    @all(dτ_Rho) = Vpdτ*max_lxyz/Re/@Mu_eff()
    @all(Gdτ)    = Vpdτ^2/@all(dτ_Rho)/(r+2)
    return
end

@parallel_indices (ix,iy,iz) function assign_τ!(τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array,  τxz::Data.Array,  τyz::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τzz_o::Data.Array, τxy_o::Data.Array, τxz_o::Data.Array, τyz_o::Data.Array)
    if (ix<=size(τxx,1) && iy<=size(τxx,2) && iz<=size(τxx,3))  τxx_o[ix,iy,iz] = τxx[ix,iy,iz]  end
    if (ix<=size(τyy,1) && iy<=size(τyy,2) && iz<=size(τyy,3))  τyy_o[ix,iy,iz] = τyy[ix,iy,iz]  end
    if (ix<=size(τzz,1) && iy<=size(τzz,2) && iz<=size(τzz,3))  τzz_o[ix,iy,iz] = τzz[ix,iy,iz]  end
    if (ix<=size(τxy,1) && iy<=size(τxy,2) && iz<=size(τxy,3))  τxy_o[ix,iy,iz] = τxy[ix,iy,iz]  end
    if (ix<=size(τxz,1) && iy<=size(τxz,2) && iz<=size(τxz,3))  τxz_o[ix,iy,iz] = τxz[ix,iy,iz]  end
    if (ix<=size(τyz,1) && iy<=size(τyz,2) && iz<=size(τyz,3))  τyz_o[ix,iy,iz] = τyz[ix,iy,iz]  end
    return
end

macro inn_yz_Gdτ(ix,iy,iz) esc(:( Gdτ[$ix  ,$iy+1,$iz+1] )) end
macro inn_xz_Gdτ(ix,iy,iz) esc(:( Gdτ[$ix+1,$iy  ,$iz+1] )) end
macro inn_xy_Gdτ(ix,iy,iz) esc(:( Gdτ[$ix+1,$iy+1,$iz  ] )) end
macro inn_yz_Mus(ix,iy,iz) esc(:( Mus[$ix  ,$iy+1,$iz+1] )) end
macro inn_xz_Mus(ix,iy,iz) esc(:( Mus[$ix+1,$iy  ,$iz+1] )) end
macro inn_xy_Mus(ix,iy,iz) esc(:( Mus[$ix+1,$iy+1,$iz  ] )) end
macro av_xyi_Gdτ(ix,iy,iz) esc(:( (Gdτ[$ix  ,$iy  ,$iz+1] + Gdτ[$ix+1,$iy  ,$iz+1] + Gdτ[$ix  ,$iy+1,$iz+1] + Gdτ[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_xzi_Gdτ(ix,iy,iz) esc(:( (Gdτ[$ix  ,$iy+1,$iz  ] + Gdτ[$ix+1,$iy+1,$iz  ] + Gdτ[$ix  ,$iy+1,$iz+1] + Gdτ[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_yzi_Gdτ(ix,iy,iz) esc(:( (Gdτ[$ix+1,$iy  ,$iz  ] + Gdτ[$ix+1,$iy+1,$iz  ] + Gdτ[$ix+1,$iy  ,$iz+1] + Gdτ[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_xyi_Mus(ix,iy,iz) esc(:( (Mus[$ix  ,$iy  ,$iz+1] + Mus[$ix+1,$iy  ,$iz+1] + Mus[$ix  ,$iy+1,$iz+1] + Mus[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_xzi_Mus(ix,iy,iz) esc(:( (Mus[$ix  ,$iy+1,$iz  ] + Mus[$ix+1,$iy+1,$iz  ] + Mus[$ix  ,$iy+1,$iz+1] + Mus[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro av_yzi_Mus(ix,iy,iz) esc(:( (Mus[$ix+1,$iy  ,$iz  ] + Mus[$ix+1,$iy+1,$iz  ] + Mus[$ix+1,$iy  ,$iz+1] + Mus[$ix+1,$iy+1,$iz+1])*0.25 )) end
macro inn_yz_Gr(ix,iy,iz) esc(:( @inn_yz_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro inn_xz_Gr(ix,iy,iz) esc(:( @inn_xz_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro inn_xy_Gr(ix,iy,iz) esc(:( @inn_xy_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro av_xyi_Gr(ix,iy,iz) esc(:( @av_xyi_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro av_xzi_Gr(ix,iy,iz) esc(:( @av_xzi_Gdτ($ix,$iy,$iz)/(G*dt) )) end
macro av_yzi_Gr(ix,iy,iz) esc(:( @av_yzi_Gdτ($ix,$iy,$iz)/(G*dt) )) end
@parallel_indices (ix,iy,iz) function compute_Pt_τ!(Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τzz_o::Data.Array, τxy_o::Data.Array, τxz_o::Data.Array, τyz_o::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Mus::Data.Array, Gdτ::Data.Array, r::Data.Number, G::Data.Number, dt::Data.Number, _dx::Data.Number, _dy::Data.Number, _dz::Data.Number)
    if (ix<=size(Pt,1)  && iy<=size(Pt,2)  && iz<=size(Pt,3))    Pt[ix,iy,iz] = Pt[ix,iy,iz] - r*Gdτ[ix,iy,iz]*(_dx*(Vx[ix+1,iy,iz] - Vx[ix,iy,iz]) + _dy*(Vy[ix,iy+1,iz] - Vy[ix,iy,iz]) + _dz*(Vz[ix,iy,iz+1] - Vz[ix,iy,iz]))  end
    if (ix<=size(τxx,1) && iy<=size(τxx,2) && iz<=size(τxx,3))  τxx[ix,iy,iz] = ( τxx[ix,iy,iz] + τxx_o[ix,iy,iz]*@inn_yz_Gr(ix,iy,iz) + 2.0*@inn_yz_Gdτ(ix,iy,iz)*(_dx*(Vx[ix+1,iy+1,iz+1] - Vx[ix,iy+1,iz+1])))/(1.0 + @inn_yz_Gdτ(ix,iy,iz)/@inn_yz_Mus(ix,iy,iz) + @inn_yz_Gr(ix,iy,iz))  end
    if (ix<=size(τyy,1) && iy<=size(τyy,2) && iz<=size(τyy,3))  τyy[ix,iy,iz] = ( τyy[ix,iy,iz] + τyy_o[ix,iy,iz]*@inn_xz_Gr(ix,iy,iz) + 2.0*@inn_xz_Gdτ(ix,iy,iz)*(_dy*(Vy[ix+1,iy+1,iz+1] - Vy[ix+1,iy,iz+1])))/(1.0 + @inn_xz_Gdτ(ix,iy,iz)/@inn_xz_Mus(ix,iy,iz) + @inn_xz_Gr(ix,iy,iz))  end
    if (ix<=size(τzz,1) && iy<=size(τzz,2) && iz<=size(τzz,3))  τzz[ix,iy,iz] = ( τzz[ix,iy,iz] + τzz_o[ix,iy,iz]*@inn_xy_Gr(ix,iy,iz) + 2.0*@inn_xy_Gdτ(ix,iy,iz)*(_dz*(Vz[ix+1,iy+1,iz+1] - Vz[ix+1,iy+1,iz])))/(1.0 + @inn_xy_Gdτ(ix,iy,iz)/@inn_xy_Mus(ix,iy,iz) + @inn_xy_Gr(ix,iy,iz))  end
    if (ix<=size(τxy,1) && iy<=size(τxy,2) && iz<=size(τxy,3))  τxy[ix,iy,iz] = ( τxy[ix,iy,iz] + τxy_o[ix,iy,iz]*@av_xyi_Gr(ix,iy,iz) + 2.0*@av_xyi_Gdτ(ix,iy,iz)*(0.5*(_dy*(Vx[ix+1,iy+1,iz+1] - Vx[ix+1,iy,iz+1]) + _dx*(Vy[ix+1,iy+1,iz+1] - Vy[ix,iy+1,iz+1])))) / (1.0 + @av_xyi_Gdτ(ix,iy,iz)/@av_xyi_Mus(ix,iy,iz) + @av_xyi_Gr(ix,iy,iz))  end
    if (ix<=size(τxz,1) && iy<=size(τxz,2) && iz<=size(τxz,3))  τxz[ix,iy,iz] = ( τxz[ix,iy,iz] + τxz_o[ix,iy,iz]*@av_xzi_Gr(ix,iy,iz) + 2.0*@av_xzi_Gdτ(ix,iy,iz)*(0.5*(_dz*(Vx[ix+1,iy+1,iz+1] - Vx[ix+1,iy+1,iz]) + _dx*(Vz[ix+1,iy+1,iz+1] - Vz[ix,iy+1,iz+1])))) / (1.0 + @av_xzi_Gdτ(ix,iy,iz)/@av_xzi_Mus(ix,iy,iz) + @av_xzi_Gr(ix,iy,iz))  end
    if (ix<=size(τyz,1) && iy<=size(τyz,2) && iz<=size(τyz,3))  τyz[ix,iy,iz] = ( τyz[ix,iy,iz] + τyz_o[ix,iy,iz]*@av_yzi_Gr(ix,iy,iz) + 2.0*@av_yzi_Gdτ(ix,iy,iz)*(0.5*(_dz*(Vy[ix+1,iy+1,iz+1] - Vy[ix+1,iy+1,iz]) + _dy*(Vz[ix+1,iy+1,iz+1] - Vz[ix+1,iy,iz+1])))) / (1.0 + @av_yzi_Gdτ(ix,iy,iz)/@av_yzi_Mus(ix,iy,iz) + @av_yzi_Gr(ix,iy,iz))  end
    return
end

macro av_xi_dτ_Rho(ix,iy,iz)  esc(:( (dτ_Rho[$ix  ,$iy+1,$iz+1] + dτ_Rho[$ix+1,$iy+1,$iz+1])*0.5 )) end
macro av_yi_dτ_Rho(ix,iy,iz)  esc(:( (dτ_Rho[$ix+1,$iy  ,$iz+1] + dτ_Rho[$ix+1,$iy+1,$iz+1])*0.5 )) end
macro av_zi_dτ_Rho(ix,iy,iz)  esc(:( (dτ_Rho[$ix+1,$iy+1,$iz  ] + dτ_Rho[$ix+1,$iy+1,$iz+1])*0.5 )) end
@parallel_indices (ix,iy,iz) function compute_V!(Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, dτ_Rho::Data.Array, _dx::Data.Number, _dy::Data.Number, _dz::Data.Number, nx_1, nx_2, ny_1, ny_2, nz_1, nz_2)
    if (ix<=nx_1 && iy<=ny_2 && iz<=nz_2)  Vx[ix+1,iy+1,iz+1] = Vx[ix+1,iy+1,iz+1] + (_dx*(τxx[ix+1,iy,iz] - τxx[ix,iy,iz]) + _dy*(τxy[ix,iy+1,iz] - τxy[ix,iy,iz]) + _dz*(τxz[ix,iy,iz+1] - τxz[ix,iy,iz]) - _dx*(Pt[ix+1,iy+1,iz+1] - Pt[ix,iy+1,iz+1]))*@av_xi_dτ_Rho(ix,iy,iz)  end
    if (ix<=nx_2 && iy<=ny_1 && iz<=nz_2)  Vy[ix+1,iy+1,iz+1] = Vy[ix+1,iy+1,iz+1] + (_dy*(τyy[ix,iy+1,iz] - τyy[ix,iy,iz]) + _dx*(τxy[ix+1,iy,iz] - τxy[ix,iy,iz]) + _dz*(τyz[ix,iy,iz+1] - τyz[ix,iy,iz]) - _dy*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy,iz+1]))*@av_yi_dτ_Rho(ix,iy,iz)  end
    if (ix<=nx_2 && iy<=ny_2 && iz<=nz_1)  Vz[ix+1,iy+1,iz+1] = Vz[ix+1,iy+1,iz+1] + (_dz*(τzz[ix,iy,iz+1] - τzz[ix,iy,iz]) + _dx*(τxz[ix+1,iy,iz] - τxz[ix,iy,iz]) + _dy*(τyz[ix,iy+1,iz] - τyz[ix,iy,iz]) - _dz*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy+1,iz]))*@av_zi_dτ_Rho(ix,iy,iz)  end
    return
end

@parallel_indices (ix,iy,iz) function compute_Res!(∇V::Data.Array, Rx::Data.Array, Ry::Data.Array, Rz::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, _dx::Data.Number, _dy::Data.Number, _dz::Data.Number)
    if (ix<=size(∇V,1) && iy<=size(∇V,2) && iz<=size(∇V,3))  ∇V[ix,iy,iz] = _dx*(Vx[ix+1,iy,iz] - Vx[ix,iy,iz]) + _dy*(Vy[ix,iy+1,iz] - Vy[ix,iy,iz]) + _dz*(Vz[ix,iy,iz+1] - Vz[ix,iy,iz]) end
    if (ix<=size(Rx,1) && iy<=size(Rx,2) && iz<=size(Rx,3))  Rx[ix,iy,iz] = _dx*(τxx[ix+1,iy,iz] - τxx[ix,iy,iz]) + _dy*(τxy[ix,iy+1,iz] - τxy[ix,iy,iz]) + _dz*(τxz[ix,iy,iz+1] - τxz[ix,iy,iz]) - _dx*(Pt[ix+1,iy+1,iz+1] - Pt[ix,iy+1,iz+1]) end
    if (ix<=size(Ry,1) && iy<=size(Ry,2) && iz<=size(Ry,3))  Ry[ix,iy,iz] = _dy*(τyy[ix,iy+1,iz] - τyy[ix,iy,iz]) + _dx*(τxy[ix+1,iy,iz] - τxy[ix,iy,iz]) + _dz*(τyz[ix,iy,iz+1] - τyz[ix,iy,iz]) - _dy*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy,iz+1]) end
    if (ix<=size(Rz,1) && iy<=size(Rz,2) && iz<=size(Rz,3))  Rz[ix,iy,iz] = _dz*(τzz[ix,iy,iz+1] - τzz[ix,iy,iz]) + _dx*(τxz[ix+1,iy,iz] - τxz[ix,iy,iz]) + _dy*(τyz[ix,iy+1,iz] - τyz[ix,iy,iz]) - _dz*(Pt[ix+1,iy+1,iz+1] - Pt[ix+1,iy+1,iz]) end
    return
end

@parallel_indices (iy,iz) function bc_x!(A::Data.Array)
    A[  1, iy,  iz] = A[    2,   iy,   iz]
    A[end, iy,  iz] = A[end-1,   iy,   iz]
    return
end

@parallel_indices (ix,iz) function bc_y!(A::Data.Array)
    A[ ix,  1,  iz] = A[   ix,    2,   iz]
    A[ ix,end,  iz] = A[   ix,end-1,   iz]
    return
end

@parallel_indices (ix,iy) function bc_z!(A::Data.Array)
    A[ ix,  iy,  1] = A[   ix,   iy,    2]
    A[ ix,  iy,end] = A[   ix,   iy,end-1]
    return
end

@views function Stokes3D_()
    # Physics
    lx, ly, lz = 10.0, 10.0, 10.0  # domain extends
    ξ          = 1.0               # Maxwell relaxation time
    μs0        = 1.0               # matrix viscosity
    μsi        = 1e-3              # inclusion viscosity
    G          = 1.0               # elastic shear modulus
    εbg        = 1.0               # background strain-rate
    dt         = μs0/(G*ξ)
    # Numerics
    nt         = 1#5                 # number of time steps
    iterMax    = 200#2e5               # maximum number of pseudo-transient iterations
    nout       = 2000               # error checking frequency
    Re         = 6.0π              # Reynolds number
    r          = 1.0               # Bulk to shear elastic modulus ratio
    CFL        = 0.8/sqrt(3)       # CFL number
    ε          = 1e-8              # nonlinear absolute tolerence
    # nx, ny, nz = 127, 127, 127     # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    b_width    = (16, 8, 4)         # boundary width for comm/comp overlap
    # Derived numerics
    me, dims, nprocs = init_global_grid(nx, ny, nz) # MPI initialisation
    @static if USE_GPU select_device() end    # select one GPU per MPI local rank (if >1 GPU per node)
    dx, dy, dz = lx/nx_g(), ly/ny_g(), lz/nz_g()      # cell sizes
    max_lxyz   = max(lx,ly,lz)
    Vpdτ       = min(dx,dy,dz)*CFL
    _dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz
    xc, yc, zc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny), LinRange(dz/2, lz-dz/2, nz)
    # Array allocations
    Pt         = @zeros(nx  ,ny  ,nz  )
    dτ_Rho     = @zeros(nx  ,ny  ,nz  )
    Gdτ        = @zeros(nx  ,ny  ,nz  )
    ∇V         = @zeros(nx  ,ny  ,nz  )
    Vy         = @zeros(nx  ,ny+1,nz  )
    τxx        = @zeros(nx  ,ny-2,nz-2)
    τyy        = @zeros(nx-2,ny  ,nz-2)
    τzz        = @zeros(nx-2,ny-2,nz  )
    τxy        = @zeros(nx-1,ny-1,nz-2)
    τxz        = @zeros(nx-1,ny-2,nz-1)
    τyz        = @zeros(nx-2,ny-1,nz-1)
    τxx_o      = @zeros(nx  ,ny-2,nz-2)
    τyy_o      = @zeros(nx-2,ny  ,nz-2)
    τzz_o      = @zeros(nx-2,ny-2,nz  )
    τxy_o      = @zeros(nx-1,ny-1,nz-2)
    τxz_o      = @zeros(nx-1,ny-2,nz-1)
    τyz_o      = @zeros(nx-2,ny-1,nz-1)
    Rx         = @zeros(nx-1,ny-2,nz-2)
    Ry         = @zeros(nx-2,ny-1,nz-2)
    Rz         = @zeros(nx-2,ny-2,nz-1)
    Mus2       = @zeros(nx  ,ny  ,nz  )
    Musτ       = @zeros(nx  ,ny  ,nz  )
    # Initial conditions
    Rad2       =  zeros(nx  ,ny  ,nz  )
    Vx         =  zeros(nx+1,ny  ,nz  )
    Vz         =  zeros(nx  ,ny  ,nz+1)
    Rad2      .= [(x_g(ix,dx,Rad2) +0.5*dx -0.5*lx)^2 + (y_g(iy,dy,Rad2) +0.5*dy -0.5*ly)^2 + (z_g(iz,dz,Rad2) +0.5*dz -0.5*lz)^2 for ix=1:size(Rad2,1), iy=1:size(Rad2,2), iz=1:size(Rad2,3)]
    Vx         = Data.Array( -εbg.*[(x_g(ix,dx,Vx) +0.5*dx -0.5*lx) for ix=1:size(Vx,1), iy=1:size(Vx,2), iz=1:size(Vx,3)] )
    Vz         = Data.Array(  εbg.*[(z_g(iz,dz,Vz) +0.5*dz -0.5*lz) for ix=1:size(Vz,1), iy=1:size(Vz,2), iz=1:size(Vz,3)] )
    Mus        = μs0*ones(nx,ny,nz)
    Mus[Rad2.<1.0] .= μsi
    Mus        = Data.Array(Mus)
    Mus2      .= Mus
    for ism=1:10#15
        @hide_communication b_width begin # communication/computation overlap
            @parallel smooth!(Mus2, Mus, 1.0)
            Mus, Mus2 = Mus2, Mus
            update_halo!(Mus)
        end
    end
    Musτ      .= Mus
    @hide_communication b_width begin # communication/computation overlap
        @parallel compute_maxloc!(Musτ, Mus)
        @parallel (1:size(Musτ,2), 1:size(Musτ,3)) bc_x!(Musτ)
        @parallel (1:size(Musτ,1), 1:size(Musτ,3)) bc_y!(Musτ)
        @parallel (1:size(Musτ,1), 1:size(Musτ,2)) bc_z!(Musτ)
        update_halo!(Musτ)
    end
    nx_1, nx_2, ny_1, ny_2, nz_1, nz_2 = nx-1, nx-2, ny-1, ny-2, nz-1, nz-2
    len_Rx_g   = ((nx-2-1)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    len_Ry_g   = ((nx-2-2)*dims[1]+2)*((ny-2-1)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    len_Rz_g   = ((nx-2-2)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-1)*dims[3]+2)
    len_∇V_g   = ((nx-2  )*dims[1]+2)*((ny-2  )*dims[2]+2)*((nz-2  )*dims[3]+2)
    # Preparation of visualisation
    if do_viz || do_save_viz
        if (me==0) ENV["GKSwstype"]="nul"; if do_viz !ispath("../../figures") && mkdir("../../figures") end; end
        nx_v, ny_v, nz_v = (nx-2)*dims[1], (ny-2)*dims[2], (nz-2)*dims[3]
        if (nx_v*ny_v*nz_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        Pt_v   = zeros(nx_v, ny_v, nz_v) # global array for visu
        Vz_v   = zeros(nx_v, ny_v, nz_v)
        Rz_v   = zeros(nx_v, ny_v, nz_v)
        Mus_v  = zeros(nx_v, ny_v, nz_v)
        τxz_v  = zeros(nx_v, ny_v, nz_v)
        Pt_inn = zeros(nx-2, ny-2, nz-2) # no halo local array for visu
        Vz_inn = zeros(nx-2, ny-2, nz-2)
        Rz_inn = zeros(nx-2, ny-2, nz-2)
        Mus_inn= zeros(nx-2, ny-2, nz-2)
        τxz_inn= zeros(nx-2, ny-2, nz-2)
        y_sl2, y_sl = Int(ceil((ny_g()-2)/2)), Int(ceil(ny_g()/2))
        Xi_g, Zi_g  = dx+dx/2:dx:(lx-dx-dx/2), dz+dz/2:dz:(lz-dz-dz/2) # inner points only
    end
    # Time loop
    @parallel compute_iter_params!(dτ_Rho, Gdτ, Musτ, Vpdτ, G, dt, Re, r, max_lxyz)
    t=0.0; ittot=0; evo_t=[]; evo_τzz=[]
    for it = 1:nt
        err=2*ε; iter=0; err_evo1=[]; err_evo2=[]
        @parallel assign_τ!(τxx, τyy, τzz, τxy, τxz, τyz, τxx_o, τyy_o, τzz_o, τxy_o, τxz_o, τyz_o)
        # Pseudo-transient iteration
        while err > ε && iter <= iterMax
            if (it==1 && iter==11)  tic()  end
            @parallel compute_Pt_τ!(Pt, τxx, τyy, τzz, τxy, τxz, τyz, τxx_o, τyy_o, τzz_o, τxy_o, τxz_o, τyz_o, Vx, Vy, Vz, Mus, Gdτ, r, G, dt, _dx, _dy, _dz)
            @hide_communication b_width begin # communication/computation overlap
                @parallel compute_V!(Vx, Vy, Vz, Pt, τxx, τyy, τzz, τxy, τxz, τyz, dτ_Rho, _dx, _dy, _dz, nx_1, nx_2, ny_1, ny_2, nz_1, nz_2)
                @parallel (1:size(Vy,2), 1:size(Vy,3)) bc_x!(Vy)
                @parallel (1:size(Vz,2), 1:size(Vz,3)) bc_x!(Vz)
                @parallel (1:size(Vx,1), 1:size(Vx,3)) bc_y!(Vx)
                @parallel (1:size(Vz,1), 1:size(Vz,3)) bc_y!(Vz)
                @parallel (1:size(Vx,1), 1:size(Vx,2)) bc_z!(Vx)
                @parallel (1:size(Vy,1), 1:size(Vy,2)) bc_z!(Vy)
                update_halo!(Vx, Vy, Vz)
            end
            iter += 1
            if iter % nout == 0
                @parallel compute_Res!(∇V, Rx, Ry, Rz, Vx, Vy, Vz, Pt, τxx, τyy, τzz, τxy, τxz, τyz, _dx, _dy, _dz)
                Vmin, Vmax = minimum_g(Vx), maximum_g(Vx)
                Pmin, Pmax = minimum_g(Pt), maximum_g(Pt)
                norm_Rx    = norm_g(Rx)/(Pmax-Pmin)*lx/sqrt(len_Rx_g)
                norm_Ry    = norm_g(Ry)/(Pmax-Pmin)*lx/sqrt(len_Ry_g)
                norm_Rz    = norm_g(Rz)/(Pmax-Pmin)*lx/sqrt(len_Rz_g)
                norm_∇V    = norm_g(∇V)/(Vmax-Vmin)*lx/sqrt(len_∇V_g)
                # norm_Rx = norm_g(Rx)/len_Rx_g; norm_Ry = norm_g(Ry)/len_Ry_g; norm_Rz = norm_g(Rz)/len_Rz_g; norm_∇V = norm_g(∇V)/len_∇V_g
                err = maximum([norm_Rx, norm_Ry, norm_Rz, norm_∇V])
                if isnan(err) error("NaN") end
                push!(err_evo1,maximum([norm_Rx, norm_Ry, norm_Rz, norm_∇V])); push!(err_evo2,iter)
                if (me==0) @printf("Step = %d, iter = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_Rz=%1.3e, norm_∇V=%1.3e] \n", it, iter, err, norm_Rx, norm_Ry, norm_Rz, norm_∇V) end
            end
        end
        ittot += iter; t += dt
        # push!(evo_t, t); push!(evo_τzz, maximum(τzz))
    end
    # Performance
    wtime    = toc()
    A_eff    = (10*2 + 6*1 + 3)/1e9*nx*ny*nz*sizeof(Data.Number) # Effective main memory access per iteration [GB] (Lower bound of required memory access: Te has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    wtime_it = wtime/(ittot-10)                           # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                             # Effective memory throughput [GB/s]
    if (me==0) @printf("Total iters = %d (%d steps), time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", ittot, nt, wtime, round(T_eff, sigdigits=3)) end
    # Visualisation
    if do_viz || do_save_viz
        Pt_inn .= inn(Pt);   gather!(Pt_inn, Pt_v)
        Vz_inn .= av_zi(Vz); gather!(Vz_inn, Vz_v)
        Rz_inn .= av_za(Rz); gather!(Rz_inn, Rz_v)
        Mus_inn.= inn(Mus);  gather!(Mus_inn, Mus_v)
        τxz_inn.= av_xza(τxz); gather!(τxz_inn, τxz_v)
        if me==0 && do_viz
            p1 = heatmap(Xi_g, Zi_g, Pt_v[:,y_sl,:]', aspect_ratio=1, xlims=(Xi_g[1],Xi_g[end]), zlims=(Zi_g[1],Zi_g[end]), c=:viridis, title="Pressure")
            p2 = heatmap(Xi_g, Zi_g, Vz_v[:,y_sl,:]', aspect_ratio=1, xlims=(Xi_g[1],Xi_g[end]), zlims=(Zi_g[1],Zi_g[end]), c=:viridis, title="Vz")
            p4 = heatmap(Xi_g, Zi_g, log10.(abs.(Rz_v[:,y_sl2,:]')), aspect_ratio=1,  xlims=(Xi_g[1],Xi_g[end]), zlims=(Zi_g[1],Zi_g[end]), c=:viridis, title="log10(Rz)")
            #p5 = plot(err_evo2,err_evo1, legend=false, xlabel="# iterations", ylabel="log10(error)", linewidth=2, markershape=:circle, markersize=3, labels="max(error)", yaxis=:log10)
            p3 = plot(evo_t, evo_τzz, legend=false, xlabel="time", ylabel="max(τzz)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
               #plot!(evo_t, 2.0.*εbg.*μs0.*(1.0.-exp.(.-evo_t.*G./μs0)), linewidth=2.0) # analytical solution
            plot(p1, p2, p4, p3)
            savefig("../../figures/Stokes_3D_ve_perf_xpu_$(nx_g()).png")
        end
    end
    if me==0 && do_save
        !ispath("../../output") && mkdir("../../output")
        open("../../output/out_Stokes3D_ve_xpu_perf.txt","a") do io
            println(io, "$(nprocs) $(nx_g()) $(ny_g()) $(nz_g()) $(ittot) $(wtime) $(A_eff) $(wtime_it) $(T_eff)")
        end
    end
    if me==0 && do_save_viz
        !ispath("../../out_visu") && mkdir("../../out_visu")
        matwrite("../../out_visu/Stokes_3D_ve_perf_xpu.mat", Dict("Pt_3D"=> Pt_v, "Mus_3D"=> Mus_v, "Txz_3D"=> τxz_v, "Vz_3D"=> Vz_v, "dx_3D"=> dx, "dy_3D"=> dy, "dz_3D"=> dz); compress = true)
    end
    finalize_global_grid()
    return xc, yc, zc, Pt
end

if use_return
    xc, yc, zc, P = Stokes3D_();
else
    Stokes3D = begin Stokes3D_(); return; end
end


text:
struct ImGuiViewportDataGlfw
    Window::Ptr{Cvoid}
    WindowOwned::Bool
    IgnoreWindowPosEventFrame::Cint
    IgnoreWindowSizeEventFrame::Cint
end

function Base.getproperty(x::Ptr{ImGuiViewportDataGlfw}, f::Symbol)
    f === :Window && return Ptr{Ptr{Cvoid}}(x + fieldoffset(ImGuiViewportDataGlfw,1))
    f === :WindowOwned && return Ptr{Bool}(x + fieldoffset(ImGuiViewportDataGlfw,2))
    f === :IgnoreWindowPosEventFrame && return Ptr{Cint}(x + fieldoffset(ImGuiViewportDataGlfw,3))
    f === :IgnoreWindowSizeEventFrame && return Ptr{Cint}(x + fieldoffset(ImGuiViewportDataGlfw,4))
    return getfield(x, f)
end

function Base.setproperty!(x::Ptr{ImGuiViewportDataGlfw}, f::Symbol, v)
    unsafe_store!(getproperty(x, f), v)
end

function ImGui_ImplGlfw_CreateWindow(viewport::Ptr{ImGuiViewport})
    io::Ptr{ImGuiIO} = igGetIO()
    ctx = unsafe_pointer_to_objref(unsafe_load(io.UserData))

    data::Ptr{ImGuiViewportDataGlfw} = Libc.malloc(sizeof(ImGuiViewportDataGlfw))
    data.Window = C_NULL
    data.WindowOwned = true
    data.IgnoreWindowPosEventFrame = -1
    data.IgnoreWindowSizeEventFrame = -1
    viewport.PlatformUserData = data

    GLFW.WindowHint(GLFW.VISIBLE, false)
    GLFW.WindowHint(GLFW.FOCUSED, false)
    GLFW.WindowHint(GLFW.FOCUS_ON_SHOW, false)
    GLFW.WindowHint(GLFW.DECORATED, !(unsafe_load(viewport.Flags) & ImGuiViewportFlags_NoDecoration == ImGuiViewportFlags_NoDecoration))
    GLFW.WindowHint(GLFW.FLOATING, unsafe_load(viewport.Flags) & ImGuiViewportFlags_TopMost == ImGuiViewportFlags_TopMost)

    share_window = GLFW.Window(ctx.Window)
    wx, wy = trunc(Cint, unsafe_load(viewport.Size.x)), trunc(Cint, unsafe_load(viewport.Size.y))
    data_Window = GLFW.CreateWindow(wx, wy, "No Title Yet", GLFW.Monitor(C_NULL), share_window)
    viewport.PlatformHandle = data_Window.handle
    if Sys.iswindows()
        viewport.PlatformHandleRaw = ccall((:glfwGetWin32Window, GLFW.libglfw), Ptr{Cvoid}, (GLFW.Window,), data_Window)
    end
    data.Window = data_Window.handle

    GLFW.SetWindowPos(data_Window, trunc(Cint, unsafe_load(viewport.Pos.x)), trunc(Cint, unsafe_load(viewport.Pos.y)))

    # install callbacks
    GLFW.SetMouseButtonCallback(data_Window, ImGui_ImplGlfw_MouseButtonCallback)
    GLFW.SetScrollCallback(data_Window, ImGui_ImplGlfw_ScrollCallback)
    GLFW.SetKeyCallback(data_Window, ImGui_ImplGlfw_KeyCallback)
    GLFW.SetCharCallback(data_Window, ImGui_ImplGlfw_CharCallback)
    GLFW.SetWindowCloseCallback(data_Window, ImGui_ImplGlfw_WindowCloseCallback)
    GLFW.SetWindowPosCallback(data_Window, ImGui_ImplGlfw_WindowPosCallback)
    GLFW.SetWindowSizeCallback(data_Window, ImGui_ImplGlfw_WindowSizeCallback)

    GLFW.MakeContextCurrent(data_Window)
    GLFW.SwapInterval(0)

    return nothing
end

function ImGui_ImplGlfw_DestroyWindow(viewport::Ptr{ImGuiViewport})
    io::Ptr{ImGuiIO} = igGetIO()
    ctx = unsafe_pointer_to_objref(unsafe_load(io.UserData))

    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    if data != C_NULL
        if unsafe_load(data.WindowOwned)
            # release any keys that were pressed in the window being destroyed and are still held down,
            # because we will not receive any release events after window is destroyed.
            win = unsafe_load(data.Window)
            for i = 0:length(ctx.KeyOwnerWindows)-1
                if ctx.KeyOwnerWindows[i+1].handle == win
                    ImGui_ImplGlfw_KeyCallback(GLFW.Window(win), i, 0, GLFW.RELEASE, 0)
                end
            end
            GLFW.DestroyWindow(GLFW.Window(win))
        end
        data.Window = C_NULL
        Libc.free(data)
    end

    viewport.PlatformUserData = C_NULL
    viewport.PlatformHandle = C_NULL

    return nothing
end

function ImGui_ImplGlfw_ShowWindow(viewport::Ptr{ImGuiViewport})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)

    if Sys.iswindows()
        # TODO: impl
    end

    GLFW.ShowWindow(GLFW.Window(unsafe_load(data.Window)))

    return nothing
end

function ImGui_ImplGlfw_GetWindowPos(viewport::Ptr{ImGuiViewport})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    x, y = GLFW.GetWindowPos(GLFW.Window(unsafe_load(data.Window)))
    return ImVec2(Cfloat(x), Cfloat(y))
end

function ImGui_ImplGlfw_SetWindowPos(viewport::Ptr{ImGuiViewport}, pos::ImVec2)
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    data.IgnoreWindowPosEventFrame = igGetFrameCount()
    x, y = trunc(Cint, pos.x), trunc(Cint, pos.y)
    GLFW.SetWindowPos(GLFW.Window(unsafe_load(data.Window)), x, y)
    return nothing
end

function ImGui_ImplGlfw_GetWindowSize(viewport::Ptr{ImGuiViewport})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    w, h = GLFW.GetWindowSize(GLFW.Window(unsafe_load(data.Window)))
    return ImVec2(Cfloat(w), Cfloat(h))
end

function ImGui_ImplGlfw_SetWindowSize(viewport::Ptr{ImGuiViewport}, size::ImVec2)
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    data.IgnoreWindowSizeEventFrame = igGetFrameCount()
    x, y = trunc(Cint, size.x), trunc(Cint, size.y)
    GLFW.SetWindowSize(GLFW.Window(unsafe_load(data.Window)), x, y)
    return nothing
end

function ImGui_ImplGlfw_SetWindowTitle(viewport::Ptr{ImGuiViewport}, title::Ptr{Cchar})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    glfwSetWindowTitle(GLFW.Window(unsafe_load(data.Window)), title)
    return nothing
end

function ImGui_ImplGlfw_SetWindowFocus(viewport::Ptr{ImGuiViewport})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    GLFW.FocusWindow(GLFW.Window(unsafe_load(data.Window)))
    return nothing
end

function ImGui_ImplGlfw_GetWindowFocus(viewport::Ptr{ImGuiViewport})::Bool
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    return GLFW.GetWindowAttrib(GLFW.Window(unsafe_load(data.Window)), GLFW.FOCUSED) != 0
end

function ImGui_ImplGlfw_GetWindowMinimized(viewport::Ptr{ImGuiViewport})::Bool
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    return GLFW.GetWindowAttrib(GLFW.Window(unsafe_load(data.Window)), GLFW.ICONIFIED) != 0
end

function ImGui_ImplGlfw_SetWindowAlpha(viewport::Ptr{ImGuiViewport}, alpha::Cfloat)
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    GLFW.SetWindowOpacity(GLFW.Window(unsafe_load(data.Window)), alpha)
    return nothing
end

function ImGui_ImplGlfw_RenderWindow(viewport::Ptr{ImGuiViewport}, userdata::Ptr{Cvoid})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    # FIXME: pass through userdata(add a new field in ImGuiViewportDataGlfw)
    io::Ptr{ImGuiIO} = igGetIO()
    ctx = unsafe_pointer_to_objref(unsafe_load(io.UserData))
    GLFW.MakeContextCurrent(GLFW.Window(unsafe_load(data.Window)))
    return nothing
end

function ImGui_ImplGlfw_SwapBuffers(viewport::Ptr{ImGuiViewport}, userdata::Ptr{Cvoid})
    data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
    # FIXME: pass through userdata(add a new field in ImGuiViewportDataGlfw)
    io::Ptr{ImGuiIO} = igGetIO()
    ctx = unsafe_pointer_to_objref(unsafe_load(io.UserData))
    GLFW.MakeContextCurrent(GLFW.Window(unsafe_load(data.Window)))
    GLFW.SwapBuffers(GLFW.Window(unsafe_load(data.Window)))
    return nothing
end

function ImGui_ImplGlfw_InitPlatformInterface(ctx::Context)
    # register platform interface (will be coupled with a renderer interface)
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    platform_io.Platform_CreateWindow = @cfunction(ImGui_ImplGlfw_CreateWindow, Cvoid, (Ptr{ImGuiViewport},))
    platform_io.Platform_DestroyWindow = @cfunction(ImGui_ImplGlfw_DestroyWindow, Cvoid, (Ptr{ImGuiViewport},))
    platform_io.Platform_ShowWindow = @cfunction(ImGui_ImplGlfw_ShowWindow, Cvoid, (Ptr{ImGuiViewport},))
    platform_io.Platform_SetWindowPos = @cfunction(ImGui_ImplGlfw_SetWindowPos, Cvoid, (Ptr{ImGuiViewport}, ImVec2))
    platform_io.Platform_GetWindowPos = @cfunction(ImGui_ImplGlfw_GetWindowPos, ImVec2, (Ptr{ImGuiViewport},))
    platform_io.Platform_SetWindowSize = @cfunction(ImGui_ImplGlfw_SetWindowSize, Cvoid, (Ptr{ImGuiViewport}, ImVec2))
    platform_io.Platform_GetWindowSize = @cfunction(ImGui_ImplGlfw_GetWindowSize, ImVec2, (Ptr{ImGuiViewport},))
    platform_io.Platform_SetWindowFocus = @cfunction(ImGui_ImplGlfw_SetWindowFocus, Cvoid, (Ptr{ImGuiViewport},))
    platform_io.Platform_GetWindowFocus = @cfunction(ImGui_ImplGlfw_GetWindowFocus, Bool, (Ptr{ImGuiViewport},))
    platform_io.Platform_GetWindowMinimized = @cfunction(ImGui_ImplGlfw_GetWindowMinimized, Bool, (Ptr{ImGuiViewport},))
    platform_io.Platform_SetWindowTitle = @cfunction(ImGui_ImplGlfw_SetWindowTitle, Cvoid, (Ptr{ImGuiViewport}, Ptr{Cchar}))
    platform_io.Platform_RenderWindow = @cfunction(ImGui_ImplGlfw_RenderWindow, Cvoid, (Ptr{ImGuiViewport}, Ptr{Cvoid}))
    platform_io.Platform_SwapBuffers = @cfunction(ImGui_ImplGlfw_SwapBuffers, Cvoid, (Ptr{ImGuiViewport}, Ptr{Cvoid}))
    platform_io.Platform_SetWindowAlpha = @cfunction(ImGui_ImplGlfw_SetWindowAlpha, Cvoid, (Ptr{ImGuiViewport}, Cfloat))
    # TODO: support Vulkan
    # platform_io.Platform_CreateVkSurface = ImGui_ImplGlfw_CreateVkSurface;
    # TODO: support IME
    # platform_io.Platform_SetImeInputPos = ImGui_ImplWin32_SetImeInputPos;

    # register main window handle (which is owned by the main application, not by us)
    # this is mostly for simplicity and consistency, so that our code (e.g. mouse handling etc.) can use same logic for main and secondary viewports.
    main_viewport::Ptr{ImGuiViewport} = igGetMainViewport()
    data::Ptr{ImGuiViewportDataGlfw} = Libc.malloc(sizeof(ImGuiViewportDataGlfw))
    data.Window = ctx.Window.handle
    data.WindowOwned = false
    data.IgnoreWindowPosEventFrame = -1
    data.IgnoreWindowSizeEventFrame = -1
    main_viewport.PlatformUserData = data
    main_viewport.PlatformHandle = Ptr{Cvoid}(ctx.Window.handle)

    return nothing
end

function ImGui_ImplGlfw_ShutdownPlatformInterface(ctx::Context)
    main_viewport::Ptr{ImGuiViewport} = igGetMainViewport()
    Libc.free(unsafe_load(main_viewport.PlatformUserData))
    main_viewport.PlatformUserData = C_NULL
    return nothing
end

function ImGui_ImplOpenGL3_RenderWindow(viewport::Ptr{ImGuiViewport}, gl_ctx_ptr::Ptr{Cvoid})::Cvoid
    if !(unsafe_load(viewport.Flags) & ImGuiViewportFlags_NoRendererClear == ImGuiViewportFlags_NoRendererClear)
        glClearColor(0.0f0, 0.0f0, 0.0f0, 1.0f0)
        glClear(GL_COLOR_BUFFER_BIT)
    end
    # TODO: submit this to upstream
    dpi_scale = unsafe_load(viewport.DpiScale)
    draw_data = unsafe_load(viewport.DrawData)
    if unsafe_load(draw_data.FramebufferScale.x) != dpi_scale
        draw_data.FramebufferScale.x = dpi_scale
        draw_data.FramebufferScale.y = dpi_scale
    end
    ImGui_ImplOpenGL3_RenderDrawData(unsafe_pointer_to_objref(gl_ctx_ptr), draw_data)
    return nothing
end

function ImGui_ImplOpenGL3_InitPlatformInterface()
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    platform_io.Renderer_RenderWindow = @cfunction(ImGui_ImplOpenGL3_RenderWindow, Cvoid, (Ptr{ImGuiViewport}, Ptr{Cvoid}))
    return true
end

ImGui_ImplOpenGL3_ShutdownPlatformInterface() = igDestroyPlatformWindows()