text:
#!/bin/julia

using CSV
using DelimitedFiles

function main(args)
    input = args[1]
    output = args[2]
    k = parse(Float64, args[3])
    equity = CSV.read(input)
    range = zeros(Int, 13, 13)
    for i in 1:13
        for j in 1:13
            range[i, j] = if equity[i, j] >= k 1 else 0 end
        end
    end
    writedlm(output, ["A" "K" "Q" "J" "T" "9" "8" "7" "6" "5" "4" "3" "2"], ",")
    open(output, "a") do io
        writedlm(io, range, ",")
    end
end

main(ARGS)

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
module TestLibPQ
using Test
using LibPQ
using StructDatabaseMapping
using Dates

DB_FILE = "test_db"

include("../includes/basic_test.jl")

host = get(ENV, "POSTGRES_HOST", "localhost")
port = get(ENV, "INPUT_POSTGRES_PORT", 5432)
db_name = get(ENV, "POSTGRES_DB", "sdm_test")
user = get(ENV, "POSTGRES_USER", "luciano")
password = get(ENV, "POSTGRES_PASSWORD", "")
conn_str = "host=$host port=$port user=$user dbname=$db_name"
conn_str = conn_str * (length(password) > 0 ? " password=$password" : "" )
function cleanup()
    conn = LibPQ.Connection(conn_str)
    execute(conn, "DROP TABLE book")
    execute(conn, "DROP TABLE author")

end

function test_create_tables()
    mapper = DBMapper(()->LibPQ.Connection(conn_str))

    register!(mapper, Author)
    register!(mapper, Book)
    configure_relation(mapper, Book, :author, on_delete=Cascade())
    @test (StructDatabaseMapping.create_table_query(mapper, Author)
    == "CREATE TABLE IF NOT EXISTS author (" *
       "age INTEGER NOT NULL, " *
       "country VARCHAR NOT NULL, " *
       "date TIMESTAMP NOT NULL, " *
       "id SERIAL PRIMARY KEY, " *
       "name VARCHAR NOT NULL)")
    @test (StructDatabaseMapping.create_table_query(mapper, Book)
     == "CREATE TABLE IF NOT EXISTS book (" *
         "author_id INTEGER NOT NULL, " *
         "data JSON NOT NULL, " *
         "id VARCHAR PRIMARY KEY, " *
         "title VARCHAR NOT NULL, " *
         "FOREIGN KEY(author_id) REFERENCES author(id) ON DELETE CASCADE ON UPDATE NO ACTION)")
end
function test()
    test_create_tables()
    test_postgres()
end
function test_postgres()
    _test_basic_functionalities(()->LibPQ.Connection(conn_str))
end

end

ext: jl

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
struct InstrumentData <: JAXTAMData
    mission    :: Mission
    instrument :: Symbol
    obsid      :: String
    events     :: DataFrames.DataFrame
    gtis       :: DataFrames.DataFrame
    start      :: Number
    stop       :: Number
    src_ctrate :: Union{Missing, Number}
    bkg_ctrate :: Union{Missing, Number}
    header     :: DataFrames.DataFrame
end

"""
    _read_fits_hdu(fits_file::FITS, hdu_id::String; cols="auto")

Reads the HDU `hdu_id` from the loaded FITS file `fits_file`, returns the HDU data

Cannot read `BitArray` type columns due to `FITSIO` limitations
"""
function _read_fits_hdu(fits_file::FITS, hdu_id::String; cols="auto")
    # Columns read individually, instead of using the FITSIO function to read
    # all the header, as if header contains non-supported column then that function
    # just fails. This avoids that problem and throws warning for non-supported columns
    if cols == "auto"
        fits_cols_events = Array{String,1}

        fits_cols_events = FITSIO.colnames(fits_file[hdu_id])
    else
        fits_cols_events = cols
    end

    fits_hdu_data = DataFrame()

    for col in fits_cols_events
        try
            fits_col_data = read(fits_file[hdu_id], col)

            if ndims(fits_col_data) == 1
                fits_hdu_data[Symbol(col)] = fits_col_data
            end
        catch error
            @warn "$col not supported by FITSIO, skipped - $error"
        end
    end

    return fits_hdu_data
end

"""
    _read_fits_event(fits_path::String)

Reads the standard columns for timing analysis ("TIME", "PI", "GTI") from a FITS file,
returns `InstrumentData` type filled with the relevant data
"""
function _read_fits_event(mission::Mission, fits_path::String)
    @info "Loading $fits_path"
    fits_file   = FITS(fits_path)

    fits_header = read_header(fits_file["EVENTS"])

    instrument_name = Symbol(fits_header["INSTRUME"])
    fits_telescope  = lowercase(fits_header["TELESCOP"])
    @assert fits_telescope == _mission_name(mission) "Mission names do not match"
    fits_telescope  = _mission_symbol_to_type(Symbol(lowercase(fits_telescope)))

    fits_events_df = _read_fits_hdu(fits_file, "EVENTS"; cols=String["TIME", "PI"])

    fits_gtis_df = _read_fits_hdu(fits_file, "GTI")

    fits_obsid = fits_header["OBS_ID"]
    fits_start = fits_header["TSTART"]
    fits_stop  = fits_header["TSTOP"]

    fits_header_df = DataFrame()

    for (i, key) in enumerate(keys(fits_header))
        key = Symbol(replace(key, '-', '_')) # Colnames with `-` don't behave well with DataFrames/Feather
        if typeof(fits_header[i]) == Nothing || fits_header[i] == ""
            fits_header_df[Symbol(key)] =  "empty" # Have to write something, or Feahter.jl errors saving ""
        else
            fits_header_df[Symbol(key)] = fits_header[i]
        end
    end

    # Hacky fix to not create DF with col type of missing
    src_rt = Array{Union{Float64,Missing},1}(undef, 1)
    bkg_rt = Array{Union{Float64,Missing},1}(undef, 1)

    total_gti_time    = sum(fits_gtis_df[:STOP] .- fits_gtis_df[:START])
    total_event_count = size(fits_events_df, 1)

    src_rt[1] = total_event_count/total_gti_time
    bkg_rt[1] = missing # TODO: Add in support for background count rates

    fits_header_df[:SRC_RT] = src_rt
    fits_header_df[:BKG_RT] = bkg_rt
    fits_header_df[:TELESCOP] = lowercase.(fits_header_df[:TELESCOP])

    close(fits_file)
    @info "Finished FITS loading"

    return InstrumentData(fits_telescope, instrument_name, fits_obsid,
        fits_events_df, fits_gtis_df, fits_start, fits_stop,
        src_rt[1], bkg_rt[1],
        fits_header_df)
end

"""
    read_cl_fits(mission::Symbol, obs_row::DataFrames.DataFrame)

Reads in FITS data for an observation, returns a `Dict{Symbol,InstrumentData}`, with the
symbol as the instrument name. So `instrument_data[:XTI]` works for NICER, and either
`instrument_data[:FPMA]` or `instrument_data[:FPMB]` work for NuSTAR
"""
function read_cl_fits(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}, file_path::String)
    if isfile(file_path)
        @info "Found: $file_path"
    elseif isfile(file_path[1:end-3]) # Check for files ending in .evt.gz too
        file_path = file_path[1:end-1]
        @info "Found: $file_path"
    else
        throw(SystemError("attempted to open:\n\t$file_path\n\t$(file_path[1:end-3])\n", 2, nothing))
    end

    instrument_data = _read_fits_event(mission, file_path)
    return Dict(instrument_data.instrument => instrument_data)
end

"""
    _save_cl_feather(feather_dir::String, instrument::Union{String,Symbol},
        fits_events_df::DataFrame, fits_gtis_df::DataFrame, fits_meta_df::DataFrame)

Due to Feather file restrictions, cannot save all the event and GTI data in one, so
they are split up into three files: `events`, `gtis`, and `meta`. The `meta` file contains
just the mission name, obsid, and observation start and stop times
"""
function _save_cl_feather(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}, instrument::Union{String,Symbol},
        feather_dir::String, fits_events_df::DataFrames.DataFrame, fits_gtis_df::DataFrames.DataFrame, fits_meta_df::DataFrames.DataFrame;
        log=true)
    @info "Saving '$instrument' cl files to $feather_dir"
    mkpath(feather_dir)

    path_events = joinpath(feather_dir, "$(instrument)_events.feather")
    path_gtis   = joinpath(feather_dir, "$(instrument)_gtis.feather")
    path_meta   = joinpath(feather_dir, "$(instrument)_meta.feather")

    Feather.write(path_events, fits_events_df)
    Feather.write(path_gtis, fits_gtis_df)
    Feather.write(path_meta, fits_meta_df)

    if log
        _log_add(mission, obs_row,
            Dict("data" =>
                Dict(:feather_cl =>
                    Dict(instrument =>
                        Dict(
                            :path_events => path_events,
                            :path_gtis   => path_gtis,
                            :path_meta   => path_meta,
                        )
                    )
                )
            )
        )
    end

    @info "Saved $instrument feather file"

    return path_events, path_gtis, path_meta
end

function _read_cl_feather(path_events::String, path_gtis::String, path_meta::String)
    data_events = Feather.read(path_events)
    data_gtis   = Feather.read(path_gtis)
    data_meta   = Feather.read(path_meta)

    meta_missn = _mission_symbol_to_type(Symbol(data_meta[1, :TELESCOP]))
    meta_inst  = data_meta[1, :INSTRUME]
    meta_obsid = data_meta[1, :OBS_ID]
    meta_start = data_meta[1, :TSTART]
    meta_stop  = data_meta[1, :TSTOP]
    meta_srcrt = data_meta[1, :SRC_RT]
    meta_bkgrt = data_meta[1, :BKG_RT]

    return InstrumentData(meta_missn, meta_inst, meta_obsid,
        data_events, data_gtis, meta_start, meta_stop,
        meta_srcrt, meta_bkgrt, data_meta)
end

"""
    read_cl(mission::Symbol, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}; overwrite=false)

Attempts to read saved (feather) data, if none is found then the `read_cl_fits` function is ran
and the data is saved with `_save_cl_feather` for future use
"""
function read_cl(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}; overwrite=false)
    obsid       = obs_row[:obsid]
    instruments = _mission_instruments(mission)

    cl_files_feather = _log_query(mission, obs_row, "data", :feather_cl)
    cl_files_raw     = _obs_files_cl(mission, obs_row)

    total_src_ctrate = 0.0
    instrument_data = Dict{Symbol,JAXTAM.InstrumentData}()
    for instrument in instruments
        if ismissing(cl_files_feather) || !haskey(cl_files_feather, instrument) || overwrite
            file_path = cl_files_raw[instrument]

            @info "Missing feather_cl for $instrument"
            @info "Processing $instrument cl fits"

            current_instrument = read_cl_fits(mission, obs_row, file_path)[instrument]

            if size(current_instrument.events, 1) == 0
                @warn "No events found in $instrument observation $obsid"
                throw(JAXTAMError("Unable to construct InstrumentData from empty DataFrame. Observation likely has no events", :read_cl))
            else
                cl_path = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM", _log_query_path(; category=:data, kind=:feather_cl))
                path_events, path_gtis, path_meta = _save_cl_feather(mission, obs_row, instrument, cl_path,
                    current_instrument.events, current_instrument.gtis, current_instrument.header)
            end

            total_src_ctrate += current_instrument.header[1, :SRC_RT]
            instrument_data[instrument] = current_instrument
        else
            @info "Loading feather cl files for '$instrument'"
            feather_paths = cl_files_feather[instrument]
            instrument_data[instrument] = _read_cl_feather(feather_paths[:path_events], feather_paths[:path_gtis], feather_paths[:path_meta])
        end
    end

    if ismissing(_log_query(mission, obs_row, "meta", :countrates, :raw))
        total_src_ctrate = total_src_ctrate/length(instruments)
        _log_add(mission, obs_row, Dict("meta" => Dict(:countrates => Dict(:raw => total_src_ctrate))))
    end

    return instrument_data
end

function read_cl(mission::Mission, master_df::DataFrames.DataFrame, obsid::String; overwrite=false)
    obs_row = master_query(master_df, :obsid, obsid)

    return read_cl(mission, obs_row; overwrite=overwrite)
end

function read_cl(mission::Mission, obsid::String; overwrite=false)
    return read_cl(mission, master(mission), obsid; overwrite=overwrite)
end


text:
using Blend2D, ColorTypes, Images, FileIO

blRuntimeInit()

s = BLSizeI(16, 7)

image = BLImageCore()
@show blImageInitAs(image, s.w, s.h, BL_FORMAT_PRGB32)

path = BLPathCore()
blPathInit(path)
blPathMoveTo(path, 0.0, 0.0)
blPathLineTo(path, Float64(s.w-1), Float64(s.h-1))

ctx = BLContextCore()
@show blContextInit(ctx)
@show blContextBegin(ctx, image, BLContextCreateInfo())
@show blContextSetFillStyleRgba32(ctx, 0xffffffff)
@show blContextFillAll(ctx)
@show blContextSetStrokeStyleRgba32(ctx, 0xff0000ff)
@show blContextStrokePathD(ctx, path)
@show blContextEnd(ctx)

data = Blend2D.BLImageData()
@show Blend2D.blImageGetData(image, data)
@show data

arr = unsafe_wrap(Array, convert(Ptr{RGBA{N0f8}}, data.pixelData), (data.size.w, data.size.h))
arr2 = PermutedDimsArray(arr, (2, 1))

save("test2.png", arr2)


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

module PowerModelsReliability

using Compat
using JuMP
using PowerModels
PMs = PowerModels

import Compat: @__MODULE__

using Compat.LinearAlgebra
using Compat.SparseArrays

include("core/variable.jl")
include("core/constraint_template.jl")
include("core/constraint.jl")
include("core/solution.jl")
include("core/base.jl")
include("core/objective.jl")
include("core/data.jl")
include("form/acp.jl")
include("form/dcp.jl")
include("prob/tfopf.jl")
include("prob/unittfopf.jl")
include("prob/scunittfopf.jl")
end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e3a36f34-f16e-11eb-0480-43a3e559f803
begin
    using Pkg
    Pkg.activate("..")
    using Revise,
        Javis, Animations, PlutoUI, Random, Colors, FileIO, JUtils, Distributions, StatsBase
    φ = Base.MathConstants.φ
    import Luxor
end

# ╔═╡ 3d896c63-b553-49a6-a901-c34688fedeeb
md"""
# Test Video
"""

# ╔═╡ d9a4d4e5-16ca-4d26-b18d-6be136ba6ee9
@bind parole PlutoUI.TextField(default = "Ciao")

# ╔═╡ 3bb08101-95c8-4dff-9516-67226c852392
md"""
# Utility
"""

# ╔═╡ 87ace88b-2fe6-4d01-8c01-f4b8a2c31384
function ground(args...)
    background("white")
    sethue("black")
end

# ╔═╡ 975c3747-7922-4619-85c9-62bcb0a7a1fe
let
    # max_frames = 50

    # test_video = Video(350, 350)
    # Background(1:max_frames, ground)
    # # Object(1:max_frames, (args...)-> circle(O, 3, :fill))
    # for i in 1:1
    # d = make_banner(pos=Point(0,0), sizex=300, sizey=200, frames=1:max_frames, rng=43,
    # 			words=split(parole))
    # moving_object(d, start_point=Point(0,0), final_point=Point(10,0), frames=25:50)
    # end
    # render(test_video, pathname="test.gif")
end

# ╔═╡ 383c4044-8e9b-42b4-b93b-cb7578ba5dc0
let
    n_frames = 400
    n_parts = 12
    height_text = -150
    height_1 = 0
    video_width = 800
    video_height = 400
    video = Video(video_width, video_height)
    starts = [1; [n_frames ÷ n_parts * i for i in 1:(n_parts - 1)]]
    stops = [n_frames ÷ n_parts * i for i in 1:n_parts]

    bg = Background(1:n_frames, ground)

    # ALL TEXTS
    texts = [
        make_text(
            body = "Have an Idea",
            text_frames = starts[1]:stops[2],
            action_frames = 1:(stops[2] - stops[1] ÷ 3),
            size = 50,
            pos = Point(0, height_text),
            halign = :center,
        ),
        make_text(
            body = "Gather Data",
            text_frames = starts[4]:stops[6],
            action_frames = 1:(stops[1] - stops[1] ÷ 3),
            size = 50,
            pos = Point(0, height_text),
            halign = :center,
        ),
    ]

    # ALL ARROWS
    arrows = [
        Object(
            starts[3]:stops[6],
            (args...) -> arrow(Point(-150, height_1), Point(150, height_1)),
            linewidth = 5,
        ),
        Object(
            starts[9]:n_frames,
            (args...) -> arrow(Point(-100, height_1), Point(100, height_1)),
            linewidth = 5,
        ),
    ]

    act!(arrows[1], Action(1:stops[1], appear(:scale)))

    # ROUND 1
    bulb = Object(
        stops[1]:stops[6],
        (args...) -> draw_lamp(O, 40, n_frames),
        Point(-video_width ÷ 2 + 100, height_1 - 40),
    )

    act!(bulb, Action(1:stops[1], appear(:fade)))

    # ROUND 2
    lens = make_lens(pos = Point(0, height_1 - 70), size = 20, frames = starts[3]:stops[6])

    # ROUND 3

    r_p = draw_random_point(100, 200)
    real_doc = make_doc(pos = r_p, sizex = 100, sizey = 110, frames = starts[5]:n_frames)

    moving_object(
        real_doc,
        start_point = r_p,
        final_point = Point(video_width ÷ 2 - 100, height_1),
        frames = 1:stops[1],
    )

    obs = [real_doc]
    for i in 1:9
        r_p = draw_random_point(100, 200)
        ob = make_doc(pos = r_p, sizex = 100, sizey = 110, frames = starts[5]:stops[6])
        push!(obs, ob)
        moving_object(
            ob,
            start_point = r_p,
            final_point = Point(video_width ÷ 2 - 100, height_1),
            frames = 1:stops[1],
        )
    end

    moving_object(
        real_doc,
        start_point = Point(video_width ÷ 2 - 100, height_1),
        final_point = Point(-video_width ÷ 2 + 150, height_1 +
                                                    # add correction because scale does not keep the center
                                                    110 ÷ 2),
        frames = starts[3]:stops[3],
    )

    act!(real_doc, Action(starts[3]:stops[3], anim_scale(1.5)))

    banner = make_banner(
        pos = Point(video_width ÷ 2 - 150, height_1),
        sizex = 200,
        frames = starts[10]:n_frames,
        rng = 43,
        # body="DESIGN"
    )

    act!(banner, Action(1:stops[1], appear(:fade)))

    render(video, pathname = "../output/process_animation.gif")
end

# ╔═╡ Cell order:
# ╠═e3a36f34-f16e-11eb-0480-43a3e559f803
# ╟─3d896c63-b553-49a6-a901-c34688fedeeb
# ╟─d9a4d4e5-16ca-4d26-b18d-6be136ba6ee9
# ╟─975c3747-7922-4619-85c9-62bcb0a7a1fe
# ╠═383c4044-8e9b-42b4-b93b-cb7578ba5dc0
# ╟─3bb08101-95c8-4dff-9516-67226c852392
# ╟─87ace88b-2fe6-4d01-8c01-f4b8a2c31384


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

using Markdown
using InteractiveUtils

# ╔═╡ 2ad11ff8-b8c8-11eb-1c46-f3babe900e58
using Luxor

# ╔═╡ de5764b5-7203-41d5-94e1-438bcd8c75fb
md"Based on [inconvergent/differential](https://github.com/inconvergent/differential-line/tree/d951ed459a99fd87713b8b5a2aee06685119c007)"

# ╔═╡ 328c272f-4d85-4bb5-928a-8d95883498ba
function draw_points(points, fix_first_last)
	n_points = length(points)
	for i in 1:length(points)
		# don't draw circular
		if fix_first_last & (i == n_points)
			continue
		end

		next_i = i < length(points) ? i + 1 : 1
		line(points[i], points[next_i], :stroke)
	end
end

# ╔═╡ 4b35feec-702d-402f-8015-5d77ea9a3b67
function attraction_step(points, farl, nearl, step, fix_first_last)
	n_points = length(points)
	d_points = [Point(0, 0) for _ in 1:n_points]
	# for every point, attract and reject
	for i in 1:n_points
		if fix_first_last & ((i == 1) | (i == n_points))
			continue
		end

		next_i = i < n_points ? i + 1 : 1
		prev_i = i > 1 ? i - 1 : n_points

		v = points[i]

		# for every other point, attract or reject
		for j in 1:n_points
			w = points[j]
			diff = v - w
			norm = sqrt(diff.x^2 + diff.y^2)
			if j == i
				continue
			end

			if j == next_i || j == prev_i
				# it is one of the neighbors
				# don't want to be too close
				if norm >= nearl
					# attract
					d_points[i] -= diff/norm*step
				end
			else
				# reject if not too far
				# closer points reject stronger
				if norm < farl
					# reject
					d_points[i] += diff*(farl/norm-1)*step
				end
			end
		end
	end

	for i in 1:n_points
		points[i] += d_points[i]
	end
end

# ╔═╡ 9d425a51-74a4-416e-a61a-39da78af1649
function add_points(points, nearl, limit, fix_first_last)
	n_points = length(points)

	# random points
	for i in 1:n_points
		if rand() > limit
			continue
		end

		# don't add point before first
		if fix_first_last & (i == 1)
			continue
		end

		# split between i-1 and i
		prev_i = i > 1 ? i - 1 : n_points

		v = points[i]
		w = points[prev_i]
		diff = v - w
		norm = sqrt(diff.x^2 + diff.y^2)
		if norm >= nearl
			# split edge if it is long enough
			midpoint = w + (diff/2)
			insert!(points, i, midpoint)
		end
	end
end

# ╔═╡ 3b64b5ea-d488-4092-ad65-f4f692df53c6
function start_circle(radius, n_init_points)
	angles = sort(rand(n_init_points)*2pi)

	[Point(radius*cos(a), radius*sin(a)) for a in angles]
end

# ╔═╡ 542a45f3-cdee-4724-bf77-e2178b0fe8fb
function start_sinusoid(size, n_init_points)
	xs = sort(rand(n_init_points)*2pi)

	[Point(size*(x-pi), size*sin(x)) for x in xs]
end

# ╔═╡ 4fd4ddd5-7a1c-454f-938b-328205aa63dd
function start_line(size, n_init_points)
	xs = sort(rand(n_init_points)*size)

	[Point(x, x + rand()/10) for x in xs]
end

# ╔═╡ 4e291e01-6601-4102-8038-6f082e0ffedd
@draw begin
	n_init_points = 40
	n_max_points = 10000
	n_steps = 4000
	draw_every = 30
	stepsize = 0.1
	farl = 40 # maximum distance to be considered 'far'
	nearl = 2 # maximum distance to be considered 'near'
	setline(0.1)

	result = []

	# points = start_circle(10, n_init_points)
	# points = start_sinusoid(20, n_init_points)
	points = start_line(20, n_init_points)
	fix_first_last = true # not circular: first and last point are fixed

	draw_points(points, fix_first_last)

	for i in 1:n_steps
		attraction_step(points, farl, nearl, stepsize, fix_first_last)
		if length(points) < n_max_points
			add_points(points, nearl, 0.001, fix_first_last)
		end
		if i % draw_every == 0
			draw_points(points, fix_first_last)
			append!(result, [points])
		end
	end
	draw_points(points, fix_first_last)
	append!(result, [points])
end 2000 2000

# ╔═╡ Cell order:
# ╠═de5764b5-7203-41d5-94e1-438bcd8c75fb
# ╠═2ad11ff8-b8c8-11eb-1c46-f3babe900e58
# ╠═328c272f-4d85-4bb5-928a-8d95883498ba
# ╠═4b35feec-702d-402f-8015-5d77ea9a3b67
# ╠═9d425a51-74a4-416e-a61a-39da78af1649
# ╠═3b64b5ea-d488-4092-ad65-f4f692df53c6
# ╠═542a45f3-cdee-4724-bf77-e2178b0fe8fb
# ╠═4fd4ddd5-7a1c-454f-938b-328205aa63dd
# ╠═4e291e01-6601-4102-8038-6f082e0ffedd


text:
### Domains

abstract type AbstractDomain{d} end

struct EmptyDomain{d} <: AbstractDomain{d}
    states::Vector{Symbol}
end

EmptyDomain{d}() where d = EmptyDomain{d}([Symbol("x$i") for i=1:d])

Base.show(io::IO, ed::EmptyDomain{d}) where d = print(io,join(fill("∅",d)))

struct CartesianDomain{d} <: AbstractDomain{d}
    states
    min::SVector{d,Float64}
    max::SVector{d,Float64}
end

CartesianDomain(states, min::Vector, max::Vector) = CartesianDomain{length(min)}(
    states,
    SVector(min...),
    SVector(max...)
)

function Base.show(io::IO, cd::CartesianDomain{d}) where d
    rgs = [((cd.min[i]), (cd.max[i])) for i=1:d]
    print(io, replace(join(rgs, "×"), "Inf"=>"∞"))
end

function CartesianDomain(min, max)
    d = length(min)
    states = [string("x", i) for i =1:d]
    return CartesianDomain(states, min, max)
end

struct DiscreteDomain{d} <: AbstractDomain{d}
    points::Vector{Point{d}}
end

struct ProductDomain{d, T,S} <: AbstractDomain{d}
    exo::T
    endo::S
end

function Base.show(io::IO, pd::ProductDomain{d}) where d
    print(io, join( [], "×"))
end

ProductDomain(d1, d2) = ProductDomain{ndims(d1)+ndims(d2), typeof(d1), typeof(d2)}(d1,d2)

function Base.show(io::IO, pd::ProductDomain) where d
    print(io, string(pd.exo, "⊗", pd.endo))
end

⊗(dom1::AbstractDomain{d1}, dom2::AbstractDomain) where d1 where d2 = ProductDomain(dom1, dom2)

struct PDomain{d, T<:Tuple} <: Grid{d}
    domains::T
end

function Base.show(io::IO, pd::PDomain) where d
    print(io, join( [pd.domains], "×"))
end

ndims(t::AbstractDomain{d}) where d = d

ndims(dom::CartesianDomain) = length(dom.min)

discretize(dom::DiscreteDomain{d}) where d = UnstructuredGrid(dom.points)

function discretize(dom::CartesianDomain; n=Union{Int, Vector{Int}})
    if typeof(n)<:Int
        nv = fill(n, ndims(dom))
    else
        nv = n
    end
    min = dom.min
    max = dom.max
    return CartesianGrid(min, max, n)
end


text:
module Memoize

export @memoize, clear_cache, get_cache

function make_key(::Function, args...)
    args
end

function make_key(::Function, x)
    x
end

# function memoize{T}(f, ::Type{T})
#     cache = Dict{Any, T}()
#     function g(args...)::T
#         key = make_key(args...)
#         get!(cache, key) do
#             f(args...)
#         end
#     end
#     g
# end

# const fib = memoize(x::Int -> begin
#     if x == 2
#         return 2
#     end
#     if x == 1
#         return 1
#     end
#     fib(x - 1) + fib(x - 2)
# end, Int)

# @code_warntype fib(3)
# @code_warntype fib.f(3)

type FunctionAst
    func_name::Symbol
    parameters::Vector{Any}
    return_type::Any
    body::Vector{Any}
end

function decompose_parameter(x::Symbol)
    x, :Any
end

function decompose_parameter(x::Expr)
    @assert x.head == :(::)
    name, variable_type = x.args
end

function parameter_names_and_types(ast::FunctionAst)
    if isempty(ast.parameters)
        return (), ()
    else
        parameter_names, parameter_types = zip(decompose_parameter.(ast.parameters)...)
        return parameter_names, parameter_types
    end
end

function FunctionAst(ast)
    @assert ast.head == :function
    header, body = ast.args
    # Has `return_type`.
    @assert header.head == :(::)
    header_token, return_type = header.args
    @assert header_token.head == :call
    func_name = header_token.args[1]
    parameters = header_token.args[2:end]
    @assert body.head == :block

    FunctionAst(func_name, parameters, return_type, body.args)
end

function compose(x::FunctionAst)
    body = Expr(:block, x.body...)
    header_token = Expr(:call, x.func_name, x.parameters...)
    header = Expr(:(::), header_token, x.return_type)
    full = Expr(:function, header, body)
end

# ast = (quote
#     function f(x::Int, y, z)::Int
#         t1()
#         t2()
#     end
# end).args[2]

# Meta.show_sexpr(ast)
# compose(FunctionAst(ast))
# eval(compose(FunctionAst(ast)))

# ast = (quote
#     function f()::Int
#         t1()
#         t2()
#     end
# end).args[2]

# Meta.show_sexpr(ast)
# compose(FunctionAst(ast))
# eval(compose(FunctionAst(ast)))

const CACHE = Dict{Tuple{Module, Symbol}, Dict}()

function clear_cache(m::Module, func_name::Symbol)
    for cache in values(CACHE[m, func_name])
        empty!(cache)
    end
end

clear_cache(f::Function) = clear_cache(typeof(f).name.module, Base.function_name(f))

function get_cache(m::Module, func_name::Symbol)
    CACHE[m, func_name]
end

get_cache(f::Function) = get_cache(typeof(f).name.module, Base.function_name(f))

macro memoize(f)
    func_ast = FunctionAst(f)

    return_type = esc(func_ast.return_type)
    func_name = esc(func_ast.func_name)
    parameters = esc.(func_ast.parameters)
    pn, pt = parameter_names_and_types(func_ast)
    parameter_names = esc.(pn)
    parameter_types = esc.(pt)

    cache = esc(:cache)
    key = esc(:key)
    func_cache = get!(Dict{Tuple, Dict}, CACHE, (current_module(), func_ast.func_name))

    func_ast.func_name = :unmemoized
    func_ast.parameters = []
    unmemoized = esc(compose(func_ast))

    quote
        let
            local pt = ($(parameter_types...),)
            local $cache = get!(Dict{Any, $return_type}, $func_cache, pt)
            global $func_name
            function $func_name($(parameters...))::$return_type
                $unmemoized
                $key = make_key($func_name, $(parameter_names...))
                get!($(esc(:unmemoized)), $cache, $key)
            end
        end
    end
end

end


text:
@testset "Linear Stretching" begin

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        #=
        Stretching an image consisting of a linear ramp should not change the image
        if the specified minimum and maximum values match the minimum and maximum of
        the image.
        =#
        img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
        minval = minimum(img)
        maxval = maximum(img)
        ret =  adjust_histogram(img, LinearStretching(minval = minval, maxval = maxval))
        if T <: Gray{Float32} || T <: Gray{Float64}
            @test all(map((i, r) -> isapprox(i, r), img, ret))
        else
            @test ret == img
        end

        # Verify that NaN is also handled correctly.
        if T <: Gray{Float32} || T <: Gray{Float64}
            img[10] = NaN
            ret = adjust_histogram(img, LinearStretching(minval = minval, maxval = maxval))
            @test isapprox(first(img), first(ret))
            @test isapprox(last(img), last(ret))
            @test isnan(ret[10])
        end

        # Verify that the smallest and largest values match the specified minval and maxval.
        img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
        minval = minimum(img)
        maxval = maximum(img)
        ret =  adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
        @test isapprox(0, first(ret))
        @test isapprox(1, last(ret))
        @test isapprox(0, minimum(ret[.!isnan.(ret)]))
        @test isapprox(1, maximum(ret[.!isnan.(ret)]))

        # Verify that the return type matches the input type.
        img = T.(testimage("mandril_gray"))
        ret = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
        @test eltype(ret) == eltype(img)
        @test isapprox(0, minimum(ret))
        @test isapprox(1, maximum(ret))

        ret = adjust_histogram(img, LinearStretching(minval = 0.2, maxval = 0.8))
        @test eltype(ret) == eltype(img)
        @test isapprox(0.2, minimum(ret))
        @test isapprox(0.8, maximum(ret))
    end

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
        #=
        Create a color image that spans a narrow graylevel range.  Then
        quantize the 256 bins down to 32 and determine how many bins have
        non-zero counts.
        =#
        imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = colorview(RGB,imgg,imgg,imgg)
        img = T.(img)
        _, counts_before = build_histogram(img,32, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Stretch the histogram. Then quantize the 256 bins down to 32 and
        verify that all 32 bins have non-zero counts. This will confirm that
        the dynamic range of the original image has been increased.
        =#
        ret = adjust_histogram(img, LinearStretching(minval = 0, maxval = 1))
        edges, counts_after = build_histogram(ret, 32, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 32
        @test eltype(img) == eltype(ret)
    end

end


text:
"""
   add_input(state, agent_idx, post_list, config)

Create new directed edges from other agents.

# Arguments
- `state`: A tuple of the current graph and agent_list
- `agent_idx`: Agent index
- `post_list`: List of all published posts in network
- `config`: `Config` object as provided by `Config`

See also: [`Config`](@ref), [`drop_input!`](@ref)
"""
function add_input!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    post_list::AbstractArray,
    config::Config
)

    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    input_candidates = Int64[]  # neighbors of neighbors

    for neighbor in inneighbors(graph, agent_idx)
        append!(
            input_candidates,
            setdiff(
                inneighbors(graph, neighbor),
                inneighbors(graph, agent_idx),
                agent_idx
            )
        )
    end

    if length(input_candidates) == 0

        input_queue = Array{Tuple{Int64,Int64}, 1}()
        not_neighbors = setdiff(
            [1:(agent_idx - 1); (agent_idx + 1):nv(graph)],
            inneighbors(graph, agent_idx)
        )  # all non-adjacent agents, self excluded

        for candidate in not_neighbors
            if (
                abs(this_agent.opinion - agent_list[candidate].opinion)
                < config.opinion_threshs.follow
            )
                push!(input_queue, (candidate, outdegree(graph, candidate)))
            end
        end

        if (
            length(input_queue) == 0
            && indegree(graph, agent_idx) <= config.agent_props.min_input_count
            && config.mechanics.dynamic_net
        )
                set_inactive!(state, agent_idx, post_list)
                this_agent.inactive_ticks = -1
                return state
        end

        input_queue = first.(sort(input_queue, by=last, rev=true))

    else
        shuffle!(input_candidates)
        input_queue = first.(
            sort(
                collect(countmap(input_candidates)),
                by=last,
                rev=true
            )
        )  #= order neighbors by frequency of occurence
              in input_candidates descending =#

    end

    if (length(input_queue) - config.network.new_follows) < 0
        new_input_count = length(input_queue)
    else
        new_input_count = config.network.new_follows
    end

    for _ in 1:new_input_count
        new_neighbor = popfirst!(input_queue)
        add_edge!(graph, new_neighbor, agent_idx)
        if (
            (
                abs(this_agent.opinion - agent_list[new_neighbor].opinion)
                < config.opinion_threshs.follow
            )
            && outdegree(graph, agent_idx) > outdegree(graph, new_neighbor)
            && rand() < 0.5
        )
            add_edge!(graph, agent_idx, new_neighbor)
        end
    end

    return state
end

# suppress output of include()
;


text:
using Apriori

itemset1 = [1,2,4]
itemset2 = [1,2,5]
itemset3 = [1,2,7]
itemset4 = [1,3,7]
itemset5 = [1,3,4]

itemsets = [itemset1, itemset2, itemset3, itemset4, itemset5]

test_merge4 = merge_vectors(itemsets)
test_merge5 = merge_vectors(test_merge4)
print(test_merge4)
print(test_merge5)
print(dummy_dataset_biased(10,10, ['a'], ['b'], 0.5))


text:
using VortexDistributions, Test, SafeTestsets

@safetestset "Point Vortex" begin include("point_vortex_test.jl") end
@safetestset "Single Vortex Accuracy" begin include("single_vort_acc_test.jl") end
@safetestset "Multiple Vortices" begin include("multivort_test.jl") end
@safetestset "Periodic Dipole" begin include("periodic_dipole_test.jl") end
@safetestset "Creation functions" begin include("creations_test.jl") end
@safetestset "Multiple Vortices 2" begin include("multivort_test2.jl") end
@safetestset "Unwrap" begin include("test_unwrap.jl") end
@safetestset "Unwrap2D" begin include("test_unwrap2D.jl") end
@safetestset "Phase Jumps" begin include("phase_jumps.jl") end
@safetestset "Periodic Dipole" begin include("periodic_dipole_test.jl") end
@safetestset "Scalar Vortex Arguments" begin include("vortex_types.jl") end
@safetestset "Single Vortex Sphere Field" begin include("single_vort_sphere.jl") end
@safetestset "Keep Vortices" begin include("keep_vortices_test.jl") end


text:
using ..SearchProblem
using Random

puzzle_actions = Dict{String,Tuple{Int,Int}}(
        "UP" => (-1, 0),
        "DOWN" => (1, 0),
        "LEFT" => (0, -1),
        "RIGHT" => (0, 1))

mutable struct City
    position::Tuple{Float32,Float32}
    num_places::Int
    np::Int
    side::Float32
    places_positions::Array{Tuple{Float32,Float32}}
    places_neighbors::Array{Array{Int}}
    costs::Array{Float32,2}
    function City(position::Tuple{Real,Real}, num_places::Integer, np::Integer, side::Real=1)
        city = new(position, num_places, np, side)
        city.places_positions = Array{Tuple{Int,Int}}(undef, num_places)
        city.places_neighbors = fill(Int[], num_places)
        city.costs = fill(Inf32, num_places, num_places)
        return city
    end
end

@inline euclid_distance(point1, point2) = sqrt(sum(@. (point1 - point2)^2))

function reset_city!(city::City, position::Tuple{Real,Real}, rng::MersenneTwister)
    @inline connected(place1, place2) = city.costs[place1, place2] < Inf
    @inline function connect!(place1, place2)
        if !connected(place1, place2)
            @inbounds distance = euclid_distance(city.places_positions[place1], city.places_positions[place2])
            @inbounds city.costs[place1, place2] = city.costs[place2, place1] = distance + (1 + 0.1 * rand(rng, Float32))
            push!(city.places_neighbors[place1], place2)
            push!(city.places_neighbors[place2], place1)
        end
        return nothing
    end

    city.position = position
    city.places_neighbors = [Int[] for _ in 1:city.num_places]
    fill!(city.costs, Inf32)

    @inbounds for place in 1:city.num_places
        city.places_positions[place] = place == 1 ? city.position : Tuple(city.position .+ rand(rng, Float32, 2) * city.side)
        city.costs[place, place] = 0
        place > 1 ? connect!(place, place - 1) : nothing
        place == city.num_places ? connect!(place, 1) : nothing
    end
    @inbounds for place in 1:city.num_places
        nearest = sort(1:city.num_places, by=p -> euclid_distance(city.places_positions[place], city.places_positions[p]))
        for p in @views nearest[2:city.np + 1]
            connect!(place, p)
        end
    end
end

mutable struct CityNavigation <: SearchProblem.AbstractSearchProblem{Tuple{Int, Int}, String}
    num_cities::Int
    num_places::Int
    nc::Int
    np::Int
    side::Float64
    city_side::Float64
    cities::Array{City,1}
    cities_neighbors::Array{Array{Int}}
    costs::Array{Float32,2}
    rng::MersenneTwister
    start_loc::Tuple{Int,Int}
    end_loc::Tuple{Int,Int}
    function CityNavigation(num_cities::Integer, num_places::Integer, nc::Integer, np::Integer, side::Real, city_side::Real)
        cnp = new(num_cities, num_places, nc, np, side, city_side)
        cnp.rng = MersenneTwister()
        cnp.cities = [City((0, 0), num_places, np, city_side) for _ in 1:num_cities]
        cnp.cities_neighbors = fill(Int[], num_cities)
        cnp.costs = fill(Inf32, num_cities, num_cities)
        cnp.start_loc = cnp.end_loc = (0, 0)
    return cnp
    end
end

function SearchProblem.reset!(cnp::CityNavigation)
    @inline connected(cityid1, cityid2) = cnp.costs[cityid1, cityid2] < Inf
    @inline function connect!(cityid1, cityid2)
        if !connected(cityid1, cityid2)
            distance = euclid_distance(cnp.cities[cityid1].position, cnp.cities[cityid2].position)
            @inbounds cnp.costs[cityid1, cityid2] = cnp.costs[cityid2, cityid1] = distance + 2
            push!(cnp.cities_neighbors[cityid1], cityid2)
            push!(cnp.cities_neighbors[cityid2], cityid1)
        end
        return nothing
    end

    cnp.cities_neighbors = [Int[] for _ in 1:cnp.num_cities]
    fill!(cnp.costs, Inf32)
    @inbounds for cityid in 1:cnp.num_cities
        reset_city!(cnp.cities[cityid], Tuple(rand(cnp.rng, Float32, 2) * cnp.side), cnp.rng)
        cnp.costs[cityid, cityid] = 0
        cityid > 1 ? connect!(cityid, cityid - 1) : nothing
        cityid == cnp.num_cities ? connect!(cityid, 1) : nothing
    end
    @inbounds for cityid in 1:cnp.num_cities
        nearest = sort(1:cnp.num_cities, by=cid -> euclid_distance(cnp.cities[cityid].position, cnp.cities[cid].position))
        for cid in @views nearest[2:cnp.nc + 1]
            connect!(cityid, cid)
        end
    end

    cnp.start_loc = cnp.end_loc = (1, 1)
    while cnp.start_loc == cnp.end_loc
        cnp.start_loc = rand(cnp.rng, 1:cnp.num_cities), rand(cnp.rng, 1:cnp.num_places)
        cnp.end_loc = rand(cnp.rng, 1:cnp.num_cities), rand(cnp.rng, 1:cnp.num_places)
    end
end

@inline SearchProblem.start_state(cnp::CityNavigation) = cnp.start_loc

@inline function SearchProblem.cost(cnp::CityNavigation, state::Tuple{Integer,Integer}, action::String, next_state::Tuple{Integer,Integer})::Float64
    cityid1, placeid1 = state
    cityid2, placeid2 = next_state
    return cityid1 == cityid2 ? cnp.cities[cityid1].costs[placeid1, placeid2] : cnp.costs[cityid1, cityid2]
end

@inline function SearchProblem.goal_test(cnp::CityNavigation, state::Tuple{Int, Int})::Bool
    return state == cnp.end_loc
end

@inline function SearchProblem.heuristic(cnp::CityNavigation, state::Tuple{Int, Int})::Float64
    cityid, placeid = state
    dest_cityid, dest_placeid = cnp.end_loc
    return euclid_distance(cnp.cities[cityid].places_positions[placeid],
                        cnp.cities[dest_cityid].places_positions[dest_placeid])
end

function SearchProblem.successors(cnp::CityNavigation, state::Tuple{Int, Int})
    cur_city_id, cur_place_id = state
    cur_city = cnp.cities[cur_city_id]
    succ = Set{Tuple{Tuple{Int,Int},String}}()
    for neighbor_place_id in cur_city.places_neighbors[cur_place_id]
        action = "goto ($cur_city_id, $neighbor_place_id)"
        next_state = (cur_city_id, neighbor_place_id)
        push!(succ, (next_state, action))
    end
    if cur_place_id == 1
        # can go to other cities as well
        for neighbor_city_id in cnp.cities_neighbors[cur_city_id]
            action = "goto city $neighbor_city_id"
            next_state = (neighbor_city_id, 1)
            push!(succ, (next_state, action))
        end
    end
    return succ
end


text:
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("Semantics are still enforced" ) do

    nan = fill(NaN, 3)
    inf = fill(Inf, 3)
    neg = fill(-Inf, 3)

    context("NaN math adds up") do
        @fact isnan(sum(nan)) --> true
    end

    context("Inf math adds up") do
        @fact sum(inf) --> Inf
    end

    context("NaN math adds up") do
        @fact sum(neg) --> -Inf
    end
end

facts("Types correctly constructed" ) do

    here     = HereSomewhere(NaN)
    consumed = Consumed(NaN)
    nonsense = Nonsense(NaN)
    shifted  = Shifted(NaN)

    context("HereSomewhere") do
        @fact isnan(here.val) --> true
        @fact_throws HereSomewhere(1.0)
    end

    context("Consumed") do
        @fact isnan(consumed.val) --> true
        @fact_throws Consumed(1.0)
    end

    context("Nonsense") do
        @fact isnan(nonsense.val) --> true
        @fact_throws Nonsense(1.0)
    end

    context("Shifted") do
        @fact isnan(shifted.val) --> true
        @fact_throws Shifted(1.0)
    end
end


text:
using Literate

function notebook_filter(str)
  re1 = r"(?<!`)``(?!`)"  # Two backquotes not preceded by nor followed by another
  re2 = r"\[\^(\d)\]:?" # e.g. "[^1]" or "[^1]:" while capturing the digit
  str |> x -> replace(x, re1 => "\$") |> x -> replace(x, re2 => s"$^\1$")
end

examples_list = ["symmetric_strip.jl", "resistive_square_patch.jl",
                 "cross_on_dielectric_substrate.jl",
                 "square_loop_absorber.jl", "band_pass_filter.jl",
                 "cpss_optimization.jl", "cpss2.jl"]

# Adds examples as subsections to the Examples page:
function replace_unknowns(str)
    findstr = """EditURL = "<unknown>"""
    replace(str, findstr => """EditURL = "https://github.com/simonp0420/PSSFSS.jl/tree/main/docs/literate""")
end

flist = ["manual.jl"]
for file in flist
    Literate.markdown(file, "../src", codefence=("```@repl manual" => "```"), credit=false,
    postprocess=replace_unknowns)
    Literate.notebook(file, "../notebooks", preprocess=notebook_filter, execute=false)
end

for (i,file) in enumerate(examples_list)
    fnpre = splitext(file)[1]
    Literate.markdown(file, ".", credit=false)
    #Literate.markdown(file, ".", codefence=("```@example $i" => "```"), credit=false)
    Literate.notebook(file, "../notebooks", preprocess=notebook_filter, execute=false)
end

function postinclude(str)
  str = replace_unknowns(str)
  for file in examples_list
    mdfile = splitext(file)[1] * ".md"
    str *= replace_unknowns(read(mdfile, String))
    #rm(mdfile)
  end
  str
end

Literate.markdown("examples.jl", "../src", postprocess=postinclude, credit=false)


text:
InputFolder = './Images/NucleiDAB/';
OutputFolder = './Results/Images/NucleiDAB2/';
Rescale = 1;
Dilate = 2;
Fill = 1;
Lbl = 1;

@iA = '*_C0000*.tif';		% Image filter
@iL = '*_C0001*.tif';		% Image annotations filter

@fxg_sLocMax [iA] > [iS];
params.GRad = 7;
params.LocalMaxBox = [5 5];
params.ThrLocalMax = 127;
params.Polarity = -1;
/endf

@fxgs_lPatchClassify [iA, iS, iL] > [C];
params.FeatType = 'HoG';
params.BoxSize = 7;
params.Sub = 5;
params.NumBins = 8;
params.ClassifierType = 'RF';
params.ClassifierFile = './Classifiers/ClassifierDABNuc_HoGRF.mat';
params.ExportAnnotations = './Classifiers/Annotation_DABNuc_HoG.tif';
/endf

/show iA > C;
/keep C > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
module BlockFieldArrays

using Gridap.Fields
using Gridap.Arrays
using Gridap.TensorValues
using BlockArrays
using FillArrays
#using Gridap.Fields: MaskedField, MaskedFieldArray

using Gridap.Fields: BlockFieldArrayCoo
using Gridap.Fields: BlockFieldArrayCooMap
using Gridap.Arrays: BlockArrayCooMap

using Test

p = Point(1.0,2.0)
np = 4
x = fill(p,np)
z = fill(p,0)

v1 = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
f1 = MockFieldArray(v1)
v2 = Float64[1,2,3,4]
f2 = MockFieldArray(v2)

axs = (append_ranges([axes(f1)[1],axes(f2)[1]]),)
b1 = BlockFieldArrayCoo(axs,[(1,)],f1)
b2 = BlockFieldArrayCoo(axs,[(2,)],f2)

b1p = evaluate(b1,p)
r = BlockArrayCoo(axs,[(1,)],[v1,])
test_array(b1p,r)

b2p = evaluate(b2,p)
r = BlockArrayCoo(axs,[(2,)],[v2,])
test_array(b2p,r)

b1x = evaluate(b1,x)
test_array(b1x,collect(b1x))

b2x = evaluate(b2,x)
test_array(b2x,collect(b2x))

#using BenchmarkTools
##blockids = [(2,)]
##@btime BlockFieldArrayCoo($axs,$blockids,$f2)
##cache = return_cache(b1,p)
##@btime evaluate!($cache,$b1,$p)
#cache = return_cache(b1,x)
#@btime evaluate!($cache,$b1,$x)

∇b2 = Broadcasting(∇)(b2)
∇b2x = evaluate(∇b2,x)
@test ∇b2x.blockids == b2x.blockids

b2t = transpose(b2)
@test b2t.blockids[1][2] == b2.blockids[1][1]
@test size(b2t) == (1,length(b2))

b2tp = evaluate(b2t,p)
test_array(b2tp,collect(b2tp))

b2tx = evaluate(b2t,x)
test_array(b2tx,collect(b2tx))

∇b2t = Broadcasting(∇)(b2t)
∇b2tx = evaluate(∇b2t,x)
@test ∇b2tx.blockids == b2tx.blockids

# Block of Blocks

axS = (append_ranges([axs[1],axs[1]]),)
b1P = BlockFieldArrayCoo(axS,[(1,)],b1)
b1M = BlockFieldArrayCoo(axS,[(2,)],b1)
b2P = BlockFieldArrayCoo(axS,[(1,)],b2)
b2M = BlockFieldArrayCoo(axS,[(2,)],b2)

b2Pp = evaluate(b2P,p)
@test is_zero_block(b2Pp,2)
@test is_nonzero_block(b2Pp,1)
@test is_zero_block(b2Pp[Block(1)],1)
@test is_nonzero_block(b2Pp[Block(1)],2)
test_array(b2Pp[Block(1)],b2p)

b2Px = evaluate(b2P,x)
@test is_zero_block(b2Px,1,2)
@test is_nonzero_block(b2Px,1,1)
@test is_zero_block(b2Px[Block(1,1)],1,1)
@test is_nonzero_block(b2Px[Block(1,1)],1,2)
test_array(b2Px[Block(1,1)],b2x)

∇b2P = Broadcasting(∇)(b2P)
∇b2Px = evaluate(∇b2P,x)
@test ∇b2Px.blockids == b2Px.blockids

b2Pt = transpose(b2P)
b2Ptp = evaluate(b2Pt,p)

# Global optimizations
ncells = 10
cell_p = Fill(p,ncells)
cell_x = Fill(x,ncells)

cell_f1 = Fill(f1,ncells)
cell_f2 = fill(f2,ncells)

cell_axs = Fill(axs,ncells)

cell_b1 = lazy_map(BlockFieldArrayCooMap((2,),[(1,)]),cell_axs,cell_f1)
cell_b2 = lazy_map(BlockFieldArrayCooMap((2,),[(2,)]),cell_axs,cell_f2)

cell_b1p = lazy_map(evaluate,cell_b1,cell_p)
cell_b2p = lazy_map(evaluate,cell_b2,cell_p)
@test isa(cell_b2p.maps.value,BlockArrayCooMap)

cell_b1x = lazy_map(evaluate,cell_b1,cell_x)
cell_b2x = lazy_map(evaluate,cell_b2,cell_x)
@test isa(cell_b2x.maps.value,BlockArrayCooMap)

cell_∇b2 = lazy_map(Broadcasting(∇),cell_b2)
@test isa(cell_∇b2.maps.value,BlockFieldArrayCooMap)

phi(x) = 2*x
cell_phi = Fill(GenericField(phi),ncells)
cell_b2ophi = lazy_map(Broadcasting(∘),cell_b2,cell_phi)
@test isa(cell_b2ophi.maps.value,BlockFieldArrayCooMap)

face_cell = [3,2,1,2]
face_b2 = lazy_map(Reindex(cell_b2),face_cell)
@test isa(face_b2.maps.value,BlockFieldArrayCooMap)

cell_b2t = lazy_map(transpose,cell_b2)
@test isa(cell_b2t.maps.value,BlockFieldArrayCooMap)

cell_b2tx = lazy_map(evaluate,cell_b2t,cell_x)
test_array(cell_b2tx,Fill(b2tx,ncells))

cell_g2 = lazy_map(Broadcasting(Operation(-)),cell_b2)
cell_g2x = lazy_map(evaluate,cell_g2,cell_x)
@test isa(cell_g2x.maps.value,BlockArrayCooMap)
test_array(cell_g2x,-collect(cell_b2x))

cell_f = Fill(ConstantField(2.0),ncells)
cell_g2 = lazy_map(Broadcasting(Operation(*)),cell_b2,cell_f)
cell_g2x = lazy_map(evaluate,cell_g2,cell_x)
@test isa(cell_g2x.maps.value,BlockArrayCooMap)
r = map(i->2.0*i,cell_b2x)
test_array(cell_g2x,r)

cell_f = Fill(ConstantField(VectorValue(1,1)),ncells)
cell_g1 = lazy_map(Broadcasting(Operation(⋅)),cell_b1,cell_f)
cell_g1x = lazy_map(evaluate,cell_g1,cell_x)
@test isa(cell_g1x.maps.value,BlockArrayCooMap)

cell_g = lazy_map(Broadcasting(Operation(-)),cell_g1,cell_g2)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)

cell_g = lazy_map(Broadcasting(Operation(*)),cell_g1,cell_b2t)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)

cell_g = lazy_map(Broadcasting(Operation(*)),cell_b2t,cell_g1)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)

ϕ = GenericField(x->2*x)
cell_ϕ = fill(ϕ,ncells)
cell_J = lazy_map(∇,cell_ϕ)
w = rand(np)
cell_w = Fill(w,ncells)

cell_i = lazy_map(integrate,cell_b2,cell_x,cell_w,cell_J)
@test isa(cell_i.maps.value,BlockArrayCooMap)

cell_g = lazy_map(Broadcasting(Operation(*)),cell_g1,cell_b2t)
cell_i = lazy_map(integrate,cell_g,cell_x,cell_w,cell_J)
@test isa(cell_i.maps.value,BlockArrayCooMap)

f(a) = 2*a
cell_g = lazy_map(Broadcasting(Operation(f)),cell_b1)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
test_array(cell_gx,map((a)->f.(a),cell_b1x))

cell_f = Fill(ConstantField(3),ncells)
cell_fx = lazy_map(evaluate,cell_f,cell_x)
f(a,b) = 2*a + a*b
cell_g = lazy_map(Broadcasting(Operation(f)),cell_b2,cell_f)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
test_array(cell_gx,map((a,b)->f.(a,b),cell_b2x,cell_fx))

f(b,a) = 2*a + a*b
cell_g = lazy_map(Broadcasting(Operation(f)),cell_f,cell_b2)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
r = map((a,b)->f.(a,b),cell_fx,cell_b2x)
test_array(cell_gx,r)

# Blocks of blokcs

cell_axsS = lazy_map(cell_axs) do axs
  (append_ranges([axs[1],axs[1]]),)
end

cell_b1L = lazy_map(BlockFieldArrayCooMap((2,),[(1,)]),cell_axsS,cell_b1)
cell_b1R = lazy_map(BlockFieldArrayCooMap((2,),[(2,)]),cell_axsS,cell_b1)
cell_b2L = lazy_map(BlockFieldArrayCooMap((2,),[(1,)]),cell_axsS,cell_b2)
cell_b2R = lazy_map(BlockFieldArrayCooMap((2,),[(2,)]),cell_axsS,cell_b2)

cell_b1Lx = lazy_map(evaluate,cell_b1L,cell_x)
cell_b1Rx = lazy_map(evaluate,cell_b1R,cell_x)
cell_b2Lx = lazy_map(evaluate,cell_b2L,cell_x)
cell_b2Rx = lazy_map(evaluate,cell_b2R,cell_x)
@test isa(cell_b1Lx.maps.value,BlockArrayCooMap)
@test isa(cell_b2Rx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_b1Lx[1],1,1) == false
@test is_zero_block(cell_b1Lx[1],1,2) == true
@test is_zero_block(cell_b1Lx[1][Block(1,1)],1,1) == false
@test is_zero_block(cell_b1Lx[1][Block(1,1)],1,2) == true
@test is_zero_block(cell_b2Lx[1],1,1) == false
@test is_zero_block(cell_b2Lx[1],1,2) == true
@test is_zero_block(cell_b2Lx[1][Block(1,1)],1,1) == true
@test is_zero_block(cell_b2Lx[1][Block(1,1)],1,2) == false

cell_∇b1R = lazy_map(Broadcasting(∇),cell_b1R)
cell_∇b2L = lazy_map(Broadcasting(∇),cell_b2L)
cell_∇b1Rx = lazy_map(evaluate,cell_∇b1R,cell_x)
cell_∇b2Lx = lazy_map(evaluate,cell_∇b2L,cell_x)
@test isa(cell_∇b1Rx.maps.value,BlockArrayCooMap)
@test isa(cell_∇b2Lx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_∇b1Rx[1],1,1) == true
@test is_zero_block(cell_∇b1Rx[1],1,2) == false
@test is_zero_block(cell_∇b1Rx[1][Block(1,2)],1,1) == false
@test is_zero_block(cell_∇b1Rx[1][Block(1,2)],1,2) == true
@test is_zero_block(cell_∇b2Lx[1],1,1) == false
@test is_zero_block(cell_∇b2Lx[1],1,2) == true
@test is_zero_block(cell_∇b2Lx[1][Block(1,1)],1,1) == true
@test is_zero_block(cell_∇b2Lx[1][Block(1,1)],1,2) == false

cell_b1Rt = lazy_map(transpose,cell_b1R)
cell_b2Lt = lazy_map(transpose,cell_b2L)
cell_b2Rt = lazy_map(transpose,cell_b2R)
cell_b1Rtx = lazy_map(evaluate,cell_b1Rt,cell_x)
cell_b2Ltx = lazy_map(evaluate,cell_b2Lt,cell_x)
cell_b2Rtx = lazy_map(evaluate,cell_b2Rt,cell_x)
@test isa(cell_b1Rtx.maps.value,BlockArrayCooMap)
@test isa(cell_b2Ltx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_b1Rtx[1],1,1,1) == true
@test is_zero_block(cell_b1Rtx[1],1,1,2) == false
@test is_zero_block(cell_b1Rtx[1][Block(1,1,2)],1,1,1) == false
@test is_zero_block(cell_b1Rtx[1][Block(1,1,2)],1,1,2) == true
@test is_zero_block(cell_b2Ltx[1],1,1,1) == false
@test is_zero_block(cell_b2Ltx[1],1,1,2) == true
@test is_zero_block(cell_b2Ltx[1][Block(1,1,1)],1,1,1) == true
@test is_zero_block(cell_b2Ltx[1][Block(1,1,1)],1,1,2) == false

cell_∇b1Rt = lazy_map(Broadcasting(∇),cell_b1Rt)
cell_∇b2Lt = lazy_map(Broadcasting(∇),cell_b2Lt)
cell_∇b1Rtx = lazy_map(evaluate,cell_∇b1Rt,cell_x)
cell_∇b2Ltx = lazy_map(evaluate,cell_∇b2Lt,cell_x)
@test isa(cell_∇b1Rtx.maps.value,BlockArrayCooMap)
@test isa(cell_∇b2Ltx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_∇b1Rtx[1],1,1,1) == true
@test is_zero_block(cell_∇b1Rtx[1],1,1,2) == false
@test is_zero_block(cell_∇b1Rtx[1][Block(1,1,2)],1,1,1) == false
@test is_zero_block(cell_∇b1Rtx[1][Block(1,1,2)],1,1,2) == true
@test is_zero_block(cell_∇b2Ltx[1],1,1,1) == false
@test is_zero_block(cell_∇b2Ltx[1],1,1,2) == true
@test is_zero_block(cell_∇b2Ltx[1][Block(1,1,1)],1,1,1) == true
@test is_zero_block(cell_∇b2Ltx[1][Block(1,1,1)],1,1,2) == false

f(a) = 3*a
cell_g = lazy_map(Broadcasting(Operation(f)),cell_b2Lt)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_gx[1],1,1,1) == false
@test is_zero_block(cell_gx[1],1,1,2) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,1) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,2) == false
r = map(i->3.0*i,cell_b2Ltx)
test_array(cell_gx,r)

cell_f = Fill(ConstantField(3),ncells)
cell_fx = lazy_map(evaluate,cell_f,cell_x)
cell_g = lazy_map(Broadcasting(Operation(*)),cell_b2Lt,cell_f)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_gx[1],1,1,1) == false
@test is_zero_block(cell_gx[1],1,1,2) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,1) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,2) == false
r = map(i->3.0*i,cell_b2Ltx)
test_array(cell_gx,r)

cell_f = Fill(ConstantField(3),ncells)
cell_fx = lazy_map(evaluate,cell_f,cell_x)
cell_g = lazy_map(Broadcasting(Operation(*)),cell_f,cell_b2Lt)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_gx[1],1,1,1) == false
@test is_zero_block(cell_gx[1],1,1,2) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,1) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,2) == false
r = map(i->3.0*i,cell_b2Ltx)
test_array(cell_gx,r)

cell_g = lazy_map(Broadcasting(Operation(-)),cell_b2Lt,cell_b2Rt)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_gx[1],1,1,1) == false
@test is_zero_block(cell_gx[1],1,1,2) == false
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,1) == true
@test is_zero_block(cell_gx[1][Block(1,1,1)],1,1,2) == false
@test is_zero_block(cell_gx[1][Block(1,1,2)],1,1,1) == true
@test is_zero_block(cell_gx[1][Block(1,1,2)],1,1,2) == false
r = map((i,j)->i-j,cell_b2Ltx,cell_b2Rtx)
test_array(cell_gx,r)

cell_g = lazy_map(Broadcasting(Operation(*)),cell_b2Lt,cell_b2R)
cell_gx = lazy_map(evaluate,cell_g,cell_x)
@test isa(cell_gx.maps.value,BlockArrayCooMap)
@test is_zero_block(cell_gx[1],1,1,1) == true
@test is_zero_block(cell_gx[1],1,1,2) == true
@test is_zero_block(cell_gx[1],1,2,1) == false
@test is_zero_block(cell_gx[1],1,2,2) == true
@test is_zero_block(cell_gx[1][Block(1,2,1)],1,1,1) == true
@test is_zero_block(cell_gx[1][Block(1,2,1)],1,1,2) == true
@test is_zero_block(cell_gx[1][Block(1,2,1)],1,2,1) == true
@test is_zero_block(cell_gx[1][Block(1,2,1)],1,2,2) == false

#b1 = BlockFieldArrayCoo(axs,[(1,)],f1)
#b2 = BlockFieldArrayCoo(axs,[(2,)],f2)

#c = Broadcasting(Operation(*))(b2,b2t)
#@show typeof(c)
#display(evaluate(c,p))
#kk

#v = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
#f = MockField.(v)
#
##fi = first(f)
#
#mi = MaskedField(fi,false)
#test_field(mi,p,fi(p))
#test_field(mi,x,fi(x))
#
#mi = MaskedField(fi,true)
#test_field(mi,p,0*fi(p))
#test_field(mi,x,0*fi(x))
#
#m = MaskedFieldArray(axes(f),f,false)
#test_field_array(m,p,evaluate(f,p))
#test_field_array(m,x,evaluate(f,x))
#
#m = MaskedFieldArray(axes(f),f,true)
#test_field_array(m,p,0*evaluate(f,p))
#test_field_array(m,x,0*evaluate(f,x))
#
#
#v1 = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
#f1 = MockField.(v1)
#v2 = Float64[1,2,3,4]
#f2 = MockField.(v2)
#
#axs = (append_ranges([axes(f1)[1],axes(f2)[1]]),)
#
#f11 = MaskedFieldArray(axes(f1),f1,false)
#f12 = MaskedFieldArray(axes(f2),f1,true)
#
#b1 = BlockArrayCoo(axs,[(1,)],[f11],[1,-1],[f12])
#
#display(b1)
#
#display(evaluate(b1,p))
#
#
#kk

#display(m)
##test_field_array(m,p,evaluate(f,p))
##test_field_array(m,x,evaluate(f,x))
#
#m = MaskedFieldArray(axes(f),f,true)
##test_field_array(m,p,0*evaluate(f,p))
##test_field_array(m,x,0*evaluate(f,x))

#v1 = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
#f1 = MockField.(v1)
#
#v2 = Float64[1,2,3]
#f2 = MockField.(v1)

end # module


text:
using Base.Test

include("transpose.jl")

@testset "empty string" begin
    @test transpose_strings([]) == []
end

@testset "two characters in a row" begin
    @test transpose_strings(["A1"]) == ["A", "1"]
end

@testset "two characters in a column" begin
    @test transpose_strings(["A", "1"]) == ["A1"]
end

@testset "simple" begin
    @test transpose_strings(["ABC", "123"]) == [
          "A1",
          "B2",
          "C3"
        ]
end

@testset "single line" begin
    @test transpose_strings(["Single line."]) == [
          "S",
          "i",
          "n",
          "g",
          "l",
          "e",
          " ",
          "l",
          "i",
          "n",
          "e",
          "."
        ]
end

@testset "first line longer than second line" begin
    @test transpose_strings([
          "The fourth line.",
          "The fifth line."
        ]) == [
          "TT",
          "hh",
          "ee",
          "  ",
          "ff",
          "oi",
          "uf",
          "rt",
          "th",
          "h ",
          " l",
          "li",
          "in",
          "ne",
          "e.",
          "."
        ]
end

@testset "second line longer than first line" begin
    @test transpose_strings([
          "The first line.",
          "The second line."
        ]) == [
          "TT",
          "hh",
          "ee",
          "  ",
          "fs",
          "ie",
          "rc",
          "so",
          "tn",
          " d",
          "l ",
          "il",
          "ni",
          "en",
          ".e",
          " ."
        ]
end

@testset "square" begin
    @test transpose_strings([
          "HEART",
          "EMBER",
          "ABUSE",
          "RESIN",
          "TREND"
        ]) == [
          "HEART",
          "EMBER",
          "ABUSE",
          "RESIN",
          "TREND"
        ]
end

@testset "rectangle" begin
    @test transpose_strings([
          "FRACTURE",
          "OUTLINED",
          "BLOOMING",
          "SEPTETTE"
        ]) == [
          "FOBS",
          "RULE",
          "ATOP",
          "CLOT",
          "TIME",
          "UNIT",
          "RENT",
          "EDGE"
        ]
end

@testset "triangle" begin
    @test transpose_strings([
          "T",
          "EE",
          "AAA",
          "SSSS",
          "EEEEE",
          "RRRRRR"
        ]) == [
          "TEASER",
          " EASER",
          "  ASER",
          "   SER",
          "    ER",
          "     R"
        ]
end

@testset "many lines" begin
    @test transpose_strings([
          "Chor. Two households, both alike in dignity,",
          "In fair Verona, where we lay our scene,",
          "From ancient grudge break to new mutiny,",
          "Where civil blood makes civil hands unclean.",
          "From forth the fatal loins of these two foes",
          "A pair of star-cross'd lovers take their life;",
          "Whose misadventur'd piteous overthrows",
          "Doth with their death bury their parents' strife.",
          "The fearful passage of their death-mark'd love,",
          "And the continuance of their parents' rage,",
          "Which, but their children's end, naught could remove,",
          "Is now the two hours' traffic of our stage;",
          "The which if you with patient ears attend,",
          "What here shall miss, our toil shall strive to mend."
        ]) == [
          "CIFWFAWDTAWITW",
          "hnrhr hohnhshh",
          "o oeopotedi ea",
          "rfmrmash  cn t",
          ".a e ie fthow ",
          " ia fr weh,whh",
          "Trnco miae  ie",
          "w ciroitr btcr",
          "oVivtfshfcuhhe",
          " eeih a uote  ",
          "hrnl sdtln  is",
          "oot ttvh tttfh",
          "un bhaeepihw a",
          "saglernianeoyl",
          "e,ro -trsui ol",
          "h uofcu sarhu ",
          "owddarrdan o m",
          "lhg to'egccuwi",
          "deemasdaeehris",
          "sr als t  ists",
          ",ebk 'phool'h,",
          "  reldi ffd   ",
          "bweso tb  rtpo",
          "oea ileutterau",
          "t kcnoorhhnatr",
          "hl isvuyee'fi ",
          " atv es iisfet",
          "ayoior trr ino",
          "l  lfsoh  ecti",
          "ion   vedpn  l",
          "kuehtteieadoe ",
          "erwaharrar,fas",
          "   nekt te  rh",
          "ismdsehphnnosa",
          "ncuse ra-tau l",
          " et  tormsural",
          "dniuthwea'g t ",
          "iennwesnr hsts",
          "g,ycoi tkrttet",
          "n ,l r s'a anr",
          "i  ef  'dgcgdi",
          "t  aol   eoe,v",
          "y  nei sl,u; e",
          ",  .sf to l   ",
          "     e rv d  t",
          "     ; ie    o",
          "       f, r   ",
          "       e  e  m",
          "       .  m  e",
          "          o  n",
          "          v  d",
          "          e  .",
          "          ,"]
end


text:
using Sudoku, Test
#### n = 2 ####
A = sudoku(Sudoku.puzz21)
@test sudoku_check(A)

A = sudoku(Sudoku.puzz22)
@test sudoku_check(A)

A = sudoku(Sudoku.puzz23)
@test sudoku_check(A)

sudoku_print(A)

#### n = 3 ####
A = sudoku(Sudoku.puzz1)
@test sudoku_check(A)

A = sudoku(Sudoku.puzz2)
@test sudoku_check(A)

A = sudoku(Sudoku.puzz3)
@test sudoku_check(A)

sudoku_print(A)

#### n = 4 ####
A = sudoku(Sudoku.puzz41)
@test sudoku_check(A)
sudoku_print(A)


text:
using DelimitedFiles
using FastGaussQuadrature
include("standardQuadrature.jl")

for n=2:2:200
	p,w = createPointStandards(n)
	m = [p w]
	display(n)
	name = string(n)*"gauss.txt"
	writedlm(name,m)
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
using Test
@time @testset "MeshCore" begin
    include("test_core.jl")
    include("attributes_tests.jl")
    include("shapedesc_tests.jl")
end


text:
include("../common/PWGrid_v01.jl")
include("../common/gen_lattice.jl")
include("gen_dr.jl")

function test_main(Ns; verbose=false)

    #LatVecs = 6.0*diagm(ones(3))
    #LatVecs = gen_lattice_hexagonal( 7.0, ca=1.5 )
    LatVecs = gen_lattice_fcc( 4.0 )
    #LatVecs = gen_lattice_bcc( 4.0 )
    println(LatVecs)

    pw = PWGrid( Ns, LatVecs )

    dr = gen_dr( pw.R, LatVecs )

    lenx = 0.5*sum(LatVecs[:,1])
    leny = 0.5*sum(LatVecs[:,2])
    lenz = 0.5*sum(LatVecs[:,3])
    MaxDist = sqrt( lenx^2 + leny^2 + lenz^2 )
    @printf("VecLen: (%18.10f,%18.10f,%18.10f)\n", lenx, leny, lenz)
    @printf("Max dist, dr: %18.10f , %18.10f\n", MaxDist, maximum(dr))
end

test_main([4,4,4], verbose=true)
#test_main([75,75,75])


text:
# Centroid
# ---------------------------------------------------------------------------------------------
# (https://en.wikipedia.org/wiki/Centroid#Of_a_polygon)
function centroid(poly::Polygon{N, 2, T}) where {N, T}
    a = zero(T) # Scalar
    c = SVector{2,T}(0,0)
    for i ∈ 1:N-1
        subarea = poly[i] × poly[i+1]
        c += subarea*(poly[i] + poly[i+1])
        a += subarea
    end
    subarea = poly[N] × poly[1]
    c += subarea*(poly[N] + poly[1])
    a += subarea
    return Point(c/(3a))
end

centroid(tri::Triangle2D) = Point((tri[1] + tri[2] + tri[3])/3)
centroid(tri::Triangle3D) = Point((tri[1] + tri[2] + tri[3])/3)

centroid(quad8::QuadraticQuadrilateral2D) = centroid(quad8, Val(3))
function centroid(quad8::QuadraticQuadrilateral{Dim, T}, ::Val{N}) where {Dim, T, N}
    # Gauss-Legendre quadrature over a quadrilateral is used.
    # Let Q(r,s) be the interpolation function for quad8,
    #            1  1                        N   N
    # A = ∬ dA = ∫  ∫ ‖∂Q/∂r × ∂Q/∂s‖ds dr = ∑   ∑ wᵢwⱼ‖∂Q/∂r(rᵢ,sⱼ) × ∂Q/∂s(rᵢ,sⱼ)‖
    #     S      0  0                       j=1 i=1
    #                  1  N   N
    # 𝗖 = (∬ 𝘅 dA)/A = -  ∑   ∑ 𝘅ᵢⱼwᵢwⱼ‖∂Q/∂r(rᵢ,sⱼ) × ∂Q/∂s(rᵢ,sⱼ)‖
    #      S           A j=1 i=1
    w, r = gauss_legendre_quadrature(T, Val(N))
    A = zero(T)
    𝗖 = @SVector zeros(T, Dim)
    for j ∈ 1:N, i ∈ 1:N
        J = 𝗝(quad8, r[i], r[j])
        weighted_val = w[i]*w[j]*norm(view(J, :, 1) × view(J, :, 2))
        𝗖 += weighted_val * quad8(r[i], r[j]).coord
        A += weighted_val
    end
    return Point(𝗖)/A
end

centroid(tri6::QuadraticTriangle2D) = centroid(tri6, Val(6))
function centroid(tri6::QuadraticTriangle{Dim, T}, ::Val{N}) where {Dim, T, N}
    # Gauss-Legendre quadrature over a triangle is used.
    # Let F(r,s) be the interpolation function for tri6,
    #            1 1-r                       N
    # A = ∬ dA = ∫  ∫ ‖∂F/∂r × ∂F/∂s‖ds dr = ∑ wᵢ‖∂F/∂r(rᵢ,sᵢ) × ∂F/∂s(rᵢ,sᵢ)‖
    #     S      0  0                       i=1
    #                  1  N
    # 𝗖 = (∬ 𝘅 dA)/A = -  ∑ 𝘅 wᵢ‖∂F/∂r(rᵢ,sᵢ) × ∂F/∂s(rᵢ,sᵢ)‖
    #      S           A i=1
    w, r, s = gauss_legendre_quadrature(tri6, Val(N))
    A = zero(T)
    𝗖 = @SVector zeros(T, Dim)
    for i ∈ 1:N
        J = 𝗝(tri6, r[i], s[i])
        weighted_val = w[i] * norm(view(J, :, 1) × view(J, :, 2))
        𝗖 += weighted_val * tri6(r[i], s[i])
        A += weighted_val
    end
    return Point(𝗖)/A
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
const USE_GPU = haskey(ENV, "USE_GPU") ? parse(Bool, ENV["USE_GPU"]) : true
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 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
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_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

@parallel function compute_timesteps!(dτVx::Data.Array, dτVy::Data.Array, dτVz::Data.Array, dτPt::Data.Array, Musτ::Data.Array, Vsc::Data.Number, Ptsc::Data.Number, min_dxyz2::Data.Number, β_n::Data.Number, max_nxyz)
    @all(dτVx) = 1.0/6.1/(1.0+β_n)/Vsc*min_dxyz2/@av_xi(Musτ)
    @all(dτVy) = 1.0/6.1/(1.0+β_n)/Vsc*min_dxyz2/@av_yi(Musτ)
    @all(dτVz) = 1.0/6.1/(1.0+β_n)/Vsc*min_dxyz2/@av_zi(Musτ)
    @all(dτPt) = 6.1*(1.0+β_n)/Ptsc/max_nxyz*@all(Musτ)
    return
end

@parallel function compute_P!(∇V::Data.Array, Pt::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, dτPt::Data.Array, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(∇V)  = @d_xa(Vx)/dx + @d_ya(Vy)/dy + @d_za(Vz)/dz
    @all(Pt)  = @all(Pt) - @all(dτPt)*@all(∇V)
    return
end

@parallel function compute_τ!(∇V::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Mus::Data.Array, β_n::Data.Number, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(τxx) = 2.0*@inn_yz(Mus)*(@d_xi(Vx)/dx  - 1.0/3.0*@inn_yz(∇V) + β_n*@inn_yz(∇V))
    @all(τyy) = 2.0*@inn_xz(Mus)*(@d_yi(Vy)/dy  - 1.0/3.0*@inn_xz(∇V) + β_n*@inn_xz(∇V))
    @all(τzz) = 2.0*@inn_xy(Mus)*(@d_zi(Vz)/dz  - 1.0/3.0*@inn_xy(∇V) + β_n*@inn_xy(∇V))
    @all(τxy) = 2.0*@av_xyi(Mus)*(0.5*(@d_yi(Vx)/dy + @d_xi(Vy)/dx))
    @all(τxz) = 2.0*@av_xzi(Mus)*(0.5*(@d_zi(Vx)/dz + @d_xi(Vz)/dx))
    @all(τyz) = 2.0*@av_yzi(Mus)*(0.5*(@d_zi(Vy)/dz + @d_yi(Vz)/dy))
    return
end

@parallel function compute_dV!(Rx::Data.Array, Ry::Data.Array, Rz::Data.Array, dVxdτ::Data.Array, dVydτ::Data.Array, dVzdτ::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, dampX::Data.Number, dampY::Data.Number, dampZ::Data.Number, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(Rx)    = @d_xa(τxx)/dx + @d_ya(τxy)/dy + @d_za(τxz)/dz - @d_xi(Pt)/dx
    @all(Ry)    = @d_ya(τyy)/dy + @d_xa(τxy)/dx + @d_za(τyz)/dz - @d_yi(Pt)/dy
    @all(Rz)    = @d_za(τzz)/dz + @d_xa(τxz)/dx + @d_ya(τyz)/dy - @d_zi(Pt)/dz
    @all(dVxdτ) = dampX*@all(dVxdτ) + @all(Rx)
    @all(dVydτ) = dampY*@all(dVydτ) + @all(Ry)
    @all(dVzdτ) = dampZ*@all(dVzdτ) + @all(Rz)
    return
end

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, dVxdτ::Data.Array, dVydτ::Data.Array, dVzdτ::Data.Array, dτVx::Data.Array, dτVy::Data.Array, dτVz::Data.Array)
    @inn(Vx) = @inn(Vx) + @all(dτVx)*@all(dVxdτ)
    @inn(Vy) = @inn(Vy) + @all(dτVy)*@all(dVydτ)
    @inn(Vz) = @inn(Vz) + @all(dτVz)*@all(dVzdτ)
    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
    μs0        = 1.0               # matrix viscosity
    μsi        = 1e-3              # inclusion viscosity
    εbg        = 1.0               # background strain-rate
    # Numerics
    iterMax    = 6e4               # maximum number of pseudo-transient iterations
    nout       = 500               # error checking frequency
    Vdmp       = 4.0               # damping paramter for the momentum equations
    β_n        = 3.0               # numerical compressibility
    Vsc        = 1.0               # relaxation paramter for the momentum equations pseudo-timesteps limiters
    Ptsc       = 5.0               # relaxation paramter for the pressure equation pseudo-timestep limiter
    ε          = 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    = (8, 4, 4)         # boundary width for comm/comp overlap
    # Derived numerics
    me, dims   = 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
    min_dxyz2  = min(dx,dy,dz)^2
    max_nxyz   = max(nx_g(),ny_g(),nz_g())
    dampX      = 1.0-Vdmp/nx_g()   # damping term for the x-momentum equation
    dampY      = 1.0-Vdmp/ny_g()   # damping term for the y-momentum equation
    dampZ      = 1.0-Vdmp/nz_g()   # damping term for the z-momentum equation
    # Array allocations
    Pt         = @zeros(nx  ,ny  ,nz  )
    dτPt       = @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)
    Rx         = @zeros(nx-1,ny-2,nz-2)
    Ry         = @zeros(nx-2,ny-1,nz-2)
    Rz         = @zeros(nx-2,ny-2,nz-1)
    dVxdτ      = @zeros(nx-1,ny-2,nz-2)
    dVydτ      = @zeros(nx-2,ny-1,nz-2)
    dVzdτ      = @zeros(nx-2,ny-2,nz-1)
    dτVx       = @zeros(nx-1,ny-2,nz-2)
    dτVy       = @zeros(nx-2,ny-1,nz-2)
    dτVz       = @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
    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
        if (me==0) ENV["GKSwstype"]="nul"; !ispath("../../figures") && mkdir("../../figures") 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)
        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)
        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_timesteps!(dτVx, dτVy, dτVz, dτPt, Musτ, Vsc, Ptsc, min_dxyz2, β_n, max_nxyz)
    err=2*ε; iter=0; err_evo1=[]; err_evo2=[]
    while err > ε && iter <= iterMax
        if (iter==11)  tic()  end
        @parallel compute_P!(∇V, Pt, Vx, Vy, Vz, dτPt, dx, dy, dz)
        @parallel compute_τ!(∇V, τxx, τyy, τzz, τxy, τxz, τyz, Vx, Vy, Vz, Mus, β_n, dx, dy, dz)
        @parallel compute_dV!(Rx, Ry, Rz, dVxdτ, dVydτ, dVzdτ, Pt, τxx, τyy, τzz, τxy, τxz, τyz, dampX, dampY, dampZ, dx, dy, dz)
        @hide_communication b_width begin # communication/computation overlap
            @parallel compute_V!(Vx, Vy, Vz, dVxdτ, dVydτ, dVzdτ, dτVx, dτVy, dτVz)
            @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
            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("Total steps = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_Rz=%1.3e, norm_∇V=%1.3e] \n", iter, err, norm_Rx, norm_Ry, norm_Rz, norm_∇V) end
        end
    end
    # Performance
    wtime    = toc()
    A_eff    = (4*2)/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/(iter-10)                         # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                          # Effective memory throughput [GB/s]
    if (me==0) @printf("Total steps = %d, err = %1.3e, time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", iter, err, wtime, round(T_eff, sigdigits=2)) end
    # Visualisation
    if do_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)
        if (me==0)
            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)
            plot(p1, p2, p4, p5)
            savefig("../../figures/Stokes_3D0_$(nx_g()).png")
        end
    end
    if me==0 && do_save
        !ispath("../../output") && mkdir("../../output")
        open("../../output/out_Stokes3D0.txt","a") do io
            println(io, "$(nx_g()) $(ny_g()) $(nz_g()) $(iter)")
        end
    end
    finalize_global_grid()
    return
end

Stokes3D()