text:
using Statistics
using JLD2
using Printf
using Oceananigans
using Oceananigans.Units

using Oceananigans.Fields: interpolate
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis
using Oceananigans.Architectures: arch_array
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using CUDA: @allowscalar
using Oceananigans.Operators: Δzᵃᵃᶜ

include("cyclic_interpolate_utils.jl")

#####
##### Grid
#####

arch = GPU()
reference_density = 1035

latitude = (-84.375, 84.375)
Δφ = latitude[2] - latitude[1]

# 2.8125 degree resolution for the fluxes data
Nx_fluxes = 128
Ny_fluxes = 60

# 0.25 degree resolution
Nx = 1440
Ny = 675
Nz = 36

output_prefix = "annual_cycle_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_fine"

#####
##### Load forcing files roughly from CORE2 paper
##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html)
#####

using DataDeps

path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/lat_lon_bathymetry_and_fluxes/"

dh = DataDep("near_global_lat_lon",
    "Forcing data for global latitude longitude simulation",
    [path * "bathymetry_lat_lon_128x60_FP32.bin",
     path * "sea_surface_temperature_25_128x60x12.jld2",
     path * "tau_x_128x60x12.jld2",
     path * "tau_y_128x60x12.jld2"]
)

DataDeps.register(dh)

datadep"near_global_lat_lon"

#####
##### Load forcing files roughly from CORE2 paper
##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html)
#####

# filename = [:sea_surface_temperature_25_128x60x12, :tau_x_128x60x12, :tau_y_128x60x12]

# for name in filename
#     datadep_path = @datadep_str "near_global_lat_lon/" * string(name) * ".jld2"
#     file = Symbol(:file_, name)
#     @eval $file = jldopen($datadep_path)
# end

# bathymetry_data = Array{Float32}(undef, Nx_fluxes*Ny_fluxes)
# bathymetry_path = @datadep_str "near_global_lat_lon/bathymetry_lat_lon_128x60_FP32.bin"
# read!(bathymetry_path, bathymetry_data)

# bathymetry_data = bswap.(bathymetry_data) |> Array{Float64}
# bathymetry_data = reshape(bathymetry_data, Nx_fluxes, Ny_fluxes)
# bathymetry      = interpolate_fluxes(bathymetry_data, Nx_fluxes, Ny_fluxes, Nx, Ny)

# τˣ_data = zeros(Nx_fluxes, Ny_fluxes, Nmonths)
# τʸ_data = zeros(Nx_fluxes, Ny_fluxes, Nmonths)
# T★_data = zeros(Nx_fluxes, Ny_fluxes, Nmonths)

# τˣ = zeros(Nx, Ny, Nmonths)
# τʸ = zeros(Nx, Ny, Nmonths)
# T★ = zeros(Nx, Ny, Nmonths)

# for month in 1:Nmonths
#     τˣ_data[:, :, month] = file_tau_x_128x60x12["tau_x/$month"] ./ reference_density
#     τʸ_data[:, :, month] = file_tau_y_128x60x12["tau_y/$month"] ./ reference_density
#     T★_data[:, :, month] = file_sea_surface_temperature_25_128x60x12["sst25/$month"]

#     τˣ[:, :, month] = interpolate_fluxes(τˣ_data[:, :, month], Nx_fluxes, Ny_fluxes, Nx, Ny)
#     τʸ[:, :, month] = interpolate_fluxes(τʸ_data[:, :, month], Nx_fluxes, Ny_fluxes, Nx, Ny)
#     T★[:, :, month] = interpolate_fluxes(T★_data[:, :, month], Nx_fluxes, Ny_fluxes, Nx, Ny)
# end

# bathymetry = arch_array(arch, bathymetry)
# τˣ = arch_array(arch, τˣ)
# τʸ = arch_array(arch, τʸ)
# target_sea_surface_temperature = T★ = arch_array(arch, T★)

# H = 3600.0
# # H = - minimum(bathymetry)

# # Uncomment for a flat bottom:
# # bathymetry = - H .* (bathymetry .< -10)

# # stretched grid
# k_center = collect(1:Nz)
# Δz_center = @. 11.6 * 1.104^(Nz - k_center)

# const Lz = sum(Δz_center)

# z_faces = vcat([-Lz], -Lz .+ cumsum(Δz_center))
# z_faces[Nz+1] = 0

# # A spherical domain
# @show underlying_grid = LatitudeLongitudeGrid(arch,
#                                               size = (Nx, Ny, Nz),
#                                               longitude = (-180, 180),
#                                               latitude = latitude,
#                                               halo = (3, 3, 3),
#                                               z = z_faces,
#                                               precompute_metrics = true)

# grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))

# #####
# ##### Physics and model setup
# #####

# νh = 1e+4
# νz = 1e+1
# κh = 1e+3
# κz = 1e-4

# background_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz,
#                                                                        time_discretization = VerticallyImplicit())

# convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0)

# #####
# ##### Boundary conditions / constant-in-time surface forcing
# #####

# Δz_top    = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.grid)
# Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.grid)

# @inline function surface_temperature_relaxation(i, j, grid, clock, fields, p)
#     time = clock.time

#     n₁ = current_time_index(time)
#     n₂ = next_time_index(time)

#     @inbounds begin
#         T★₁ = p.T★[i, j, n₁]
#         T★₂ = p.T★[i, j, n₂]
#         T_surface = fields.T[i, j, grid.Nz]
#     end

#     T★ = cyclic_interpolate(T★₁, T★₂, time)

#     return p.λ * (T_surface - T★)
# end

# T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation,
#                                                 discrete_form = true,
#                                                 parameters = (λ = Δz_top/3days, T★ = target_sea_surface_temperature))

# @inline function wind_stress(i, j, grid, clock, fields, τ)
#     time = clock.time
#     n₁ = current_time_index(time)
#     n₂ = next_time_index(time)

#     @inbounds begin
#         τ₁ = τ[i, j, n₁]
#         τ₂ = τ[i, j, n₂]
#     end

#     return cyclic_interpolate(τ₁, τ₂, time)
# end

# u_wind_stress_bc = FluxBoundaryCondition(wind_stress, discrete_form = true, parameters = τˣ)
# v_wind_stress_bc = FluxBoundaryCondition(wind_stress, discrete_form = true, parameters = τʸ)

# @inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
# @inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]

# # Linear bottom drag:
# μ = Δz_bottom / 10days

# u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ)
# v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ)

# u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc)
# v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc)
# T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc)

# free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver)

# buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-4), constant_salinity=true)

# model = HydrostaticFreeSurfaceModel(grid = grid,
#                                     free_surface = free_surface,
#                                     momentum_advection = VectorInvariant(),
#                                     tracer_advection = WENO5(grid = underlying_grid),
#                                     coriolis = HydrostaticSphericalCoriolis(),
#                                     boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs),
#                                     buoyancy = buoyancy,
#                                     tracers = (:T, ),
#                                     closure = (background_diffusivity, convective_adjustment))

# @show model.free_surface.implicit_step_solver.matrix_iterative_solver.preconditioner_method

# #####
# ##### Initial condition:
# #####

# u, v, w = model.velocities
# η = model.free_surface.η
# T = model.tracers.T
# T .= -1

# #####
# ##### Simulation setup
# #####

# # Time-scale for gravity wave propagation across the smallest grid cell
# g = model.free_surface.gravitational_acceleration
# gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

# minimum_Δx = abs(grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ[1:grid.Ny])) * deg2rad(grid.Δλᶜᵃᵃ))
# minimum_Δy = abs(grid.radius * deg2rad(grid.Δφᵃᶜᵃ))
# wave_propagation_time_scale = min(minimum_Δx, minimum_Δy) / gravity_wave_speed

# if model.free_surface isa ExplicitFreeSurface
#     Δt = 10seconds
# else
#     Δt = 3minutes
# end

# simulation = Simulation(model, Δt = Δt, stop_time = 30years)

# start_time = [time_ns()]

# function progress(sim)
#     wall_time = (time_ns() - start_time[1]) * 1e-9

#     η = model.free_surface.η
#     u = model.velocities.u
#     @info @sprintf("Time: % 12s, iteration: %d, max(|η|): %.2e m, max(|u|): %.2e ms⁻¹, wall time: %s",
#                     prettytime(sim.model.clock.time),
#                     sim.model.clock.iteration,
#                     maximum(abs, η), maximum(abs, u),
#                     prettytime(wall_time))

#     start_time[1] = time_ns()

#     return nothing
# end

# simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

# u, v, w = model.velocities
# T = model.tracers.T
# η = model.free_surface.η

# output_fields = (; u, v, T, η)
# save_interval = 5days

# simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, (; u, v, T, η),
#                                                               schedule = TimeInterval(save_interval),
#                                                               prefix = output_prefix * "_surface",
#                                                               field_slicer = FieldSlicer(k=grid.Nz),
#                                                               force = true)

# simulation.output_writers[:bottom_fields] = JLD2OutputWriter(model, (; u, v, T),
#                                                              schedule = TimeInterval(save_interval),
#                                                              prefix = output_prefix * "_bottom",
#                                                              field_slicer = FieldSlicer(k=1),
#                                                              force = true)

# simulation.output_writers[:checkpointer] = Checkpointer(model,
#                                                         schedule = TimeInterval(1year),
#                                                         prefix = output_prefix * "_checkpoint",
#                                                         cleanup = true,
#                                                         force = true)

# # Let's goo!
# @info "Running with Δt = $(prettytime(simulation.Δt))"

# run!(simulation)

# @info """

#     Simulation took $(prettytime(simulation.run_wall_time))
#     Background diffusivity: $background_diffusivity
#     Minimum wave propagation time scale: $(prettytime(wave_propagation_time_scale))
#     Free surface: $(typeof(model.free_surface).name.wrapper)
#     Time step: $(prettytime(Δt))
# """

####
#### Visualize solution
####

using GLMakie

surface_file = jldopen(output_prefix * "_surface.jld2")
bottom_file = jldopen(output_prefix * "_bottom.jld2")

iterations = parse.(Int, keys(surface_file["timeseries/t"]))
iterations = iterations[15:end]

iter = Node(0)

ηi(iter) = surface_file["timeseries/η/" * string(iter)][:, :, 1]
ui(iter) = surface_file["timeseries/u/" * string(iter)][:, :, 1]
vi(iter) = surface_file["timeseries/v/" * string(iter)][:, :, 1]
Ti(iter) = surface_file["timeseries/T/" * string(iter)][:, :, 1]
ti(iter) = string(surface_file["timeseries/t/" * string(iter)] / day)

ubi(iter) = bottom_file["timeseries/u/" * string(iter)][:, :, 1]
vbi(iter) = bottom_file["timeseries/v/" * string(iter)][:, :, 1]

η = @lift ηi($iter)
u = @lift ui($iter)
v = @lift vi($iter)
T = @lift Ti($iter)

ub = @lift ubi($iter)
vb = @lift vbi($iter)

max_η = 4
min_η = - max_η
max_u = 0.2
min_u = - max_u
max_T = 32
min_T = 0

fig = Figure(resolution = (1200, 900))

ax = Axis(fig[1, 1], title="Free surface displacement (m)")
hm = GLMakie.heatmap!(ax, η, colorrange=(min_η, max_η), colormap=:balance)
cb = Colorbar(fig[1, 2], hm)

ax = Axis(fig[2, 1], title="Sea surface temperature (ᵒC)")
hm = GLMakie.heatmap!(ax, T, colorrange=(min_T, max_T), colormap=:thermal)
cb = Colorbar(fig[2, 2], hm)

ax = Axis(fig[1, 3], title="East-west surface velocity (m s⁻¹)")
hm = GLMakie.heatmap!(ax, u, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[1, 4], hm)

ax = Axis(fig[2, 3], title="North-south surface velocity (m s⁻¹)")
hm = GLMakie.heatmap!(ax, v, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[2, 4], hm)

ax = Axis(fig[3, 1], title="East-west bottom velocity (m s⁻¹)")
hm = GLMakie.heatmap!(ax, ub, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[3, 2], hm)

ax = Axis(fig[3, 3], title="North-south bottom velocity (m s⁻¹)")
hm = GLMakie.heatmap!(ax, vb, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[3, 4], hm)

title_str = @lift "Earth day = " * ti($iter)
ax_t = fig[0, :] = Label(fig, title_str)

GLMakie.record(fig, output_prefix * ".mp4", iterations, framerate=8) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end

display(fig)

close(surface_file)


text:
using DataStructures
using StatsBase

# Problem 1
nWashers = 5
nTrials = 1000000
success = 0
for i in 1:nTrials
    rn = rand(1:nWashers,nWashers)
    if get(countmap(rn),1,0) >= 4
        success+=1
    end
end
success/nTrials

# Problem 2
wait1 = 5
wait2 = 7
nTrials = 1000000
success = 0
for i in 1:nTrials
    rn = rand(2)
    if (0 <= (rn[2]-rn[1])*30 <= wait1) ||
       (0 <= (rn[1]-rn[2])*30 <= wait2)
        success+=1
    end
end
success/nTrials

# Problem 3
# nCars=3;  0.6666666666644826
# nCars=4;  0.6667395
# nCars=5;  0.6663580000064863
# nCars=6;  0.6665606666653675
# nCars=7;  0.6665282857209788
# nCars=8;  0.66674625
# nCars=9;  0.6667986666646687
# nCars=10; 0.6665119999960215
# nCars=11; 0.6665550909138054
# nCars=12; 0.6666493333327062
# nCars=20; 0.6666518999954775
# nCars=30; 0.6665245333273143

nCars = 7
nTrials = 1000000
success=0
for i in 1:nTrials
    rn = rand(nCars)
    rn = sort(rn)
    cn = zeros(nCars)
    cn[1]=2
    cn[nCars]=nCars-1
    for j in 2:nCars-1
        if rn[j+1]-rn[j] >= rn[j]-rn[j-1]
            cn[j]=j-1
        else
            cn[j]=j+1
        end
    end
    cnTotal = 0
    j = 1
    while j <= nCars-1
        if (cn[j]==j+1 && cn[j+1]==j)
            j+=2
            cnTotal+=2
        else
            j+=1
        end
    end
    success+=cnTotal/nCars
end
success/nTrials

ext: jl

max_stars_repo_path: julia/Digital Dice.jl

max_stars_repo_name: candychiu/randomsnippet

max_issues_repo_path: julia/Digital Dice.jl

max_issues_repo_name: candychiu/randomsnippet

max_forks_repo_path: julia/Digital Dice.jl

max_forks_repo_name: candychiu/randomsnippet


text:
#module ParallelMatrix

export SparseBinMatrix, ParallelLogic, ParallelSBM, balanced_parallelsbm
export sort_hilbert
export AtA_mul_B!

using Compat

type SparseBinMatrix
  m::Int64
  n::Int64
  mrange::UnitRange{Int32}
  nrange::UnitRange{Int32}
  ## data
  rows::Vector{Int32}
  cols::Vector{Int32}
end

function SparseBinMatrix(m, n, rows::Vector{Int32}, cols::Vector{Int32})
  length(rows) == length(cols) || throw(DimensionMismatch("length(rows) must equal length(cols)"))
  return SparseBinMatrix(m, n, minimum(rows):maximum(rows), minimum(cols):maximum(cols), rows, cols)
end
SparseBinMatrix(rows::Vector{Int32}, cols::Vector{Int32}) = SparseBinMatrix(maximum(rows), maximum(cols), rows, cols)
SparseBinMatrix(rows::Vector{Int64}, cols::Vector{Int64}) = SparseBinMatrix(convert(Vector{Int32}, rows), convert(Vector{Int32}, cols))

import Base.getindex
function getindex(sbm::SparseBinMatrix, rows::Vector{Bool})
  length(rows) == sbm.m || throw(DimensionMismatch("length(rows) must equal size(sbm,1)"))
  idx = rows[sbm.rows]
  colidx = sbm.cols[idx]
  rsum   = cumsum(rows)
  rowidx = rsum[ sbm.rows[idx] ]
  out = SparseBinMatrix(sum(rows), sbm.n, one(Int32):convert(Int32,sum(rows)), sbm.nrange, rowidx, colidx)
  return out
end

function getindex(sbm::SparseBinMatrix, rows::Vector{Bool}, cols::Colon)
  return getindex(sbm, rows)
end

function getindex(sbm::SparseBinMatrix, rows::Vector{Bool}, cols::UnitRange{Int64})
  return getindex(sbm, rows)
end

import Base.ndims
ndims(sbm::SparseBinMatrix) = 2

type ParallelLogic
  ## for parallel compute
  mblocks::Vector{UnitRange{Int32}}
  nblocks::Vector{UnitRange{Int32}}
  mblock_order::Vector{Int64}
  nblock_order::Vector{Int64}

  ## local vectors
  localm::Vector{Float64}
  localn::Vector{Float64}

  ## semaphores
  tmp::SharedVector{Float64}
  sems::Vector{SharedArray{UInt32,1}}

  ## constructor for non-shared elements
  ParallelLogic(mblocks, nblocks, mblock_order, nblock_order, localm, localn) = new(mblocks, nblocks, mblock_order, nblock_order, localm, localn)
  ParallelLogic(mblocks, nblocks, mblock_order, nblock_order, localm, localn, tmp, sems) = new(mblocks, nblocks, mblock_order, nblock_order, localm, localn, tmp, sems)
end

nonshared(A::ParallelLogic) = ParallelLogic(A.mblocks, A.nblocks, A.mblock_order, A.nblock_order, A.localm, A.localn)

type ParallelSBM
  m::Int64
  n::Int64
  pids::Vector{Int64}
  sbms::Vector{RemoteRef}  ## SparseBinMatrix
  logic::Vector{RemoteRef} ## ParallelLogic

  numblocks::Int                ## number of blocks, each has semaphore
  tmp::SharedVector{Float64}    ## for storing middle vector in A'A
  sh1::SharedVector{Float64}    ## length(sh1) = A.n
  sh2::SharedVector{Float64}    ## length(sh2) = A.n
  sems::Vector{SharedArray{UInt32,1}} ## semaphores

  ## constructors
  ParallelSBM(m, n, pids, sbms, logic, numblocks) = new(m, n, pids, sbms, logic, numblocks)
  ParallelSBM(m, n, pids, sbms, logic, numblocks, tmp, sh1, sh2, sems) = new(m, n, pids, sbms, logic, numblocks, tmp, sh1, sh2, sems)
end

nonshared(A::ParallelSBM) = ParallelSBM(A.m, A.n, A.pids, A.sbms, A.logic, A.numblocks)

function make_sems(numblocks::Int, pids::Vector{Int})
  sems = SharedArray{UInt32,1}[SharedArray(UInt32, 16, pids=pids) for i=1:numblocks]
  for sem in sems
    sem_init(sem)
  end
  return sems
end

function block_order(counts, blocks)
  bcounts = counts[blocks]
  return blocks[sortperm(bcounts, rev=true)]
end

function ParallelSBM(rows::Vector{Int}, cols::Vector{Int}, pids::Vector{Int}, weights=ones(length(pids)), m=convert(Int32, maximum(rows)), n=convert(Int32, maximum(cols)), numblocks=length(pids)*2)
  return ParallelSBM(convert(Vector{Int32}, rows), convert(Vector{Int32}, cols), pids, weights=weights, m=m, n=n, numblocks=numblocks)
end

function ParallelSBM(rows::Vector{Int32}, cols::Vector{Int32}, pids::Vector{Int}=Int[]; weights=ones(length(pids)), m=maximum(rows), n=maximum(cols), numblocks=length(pids)*2 )
  length(rows) == length(cols) || throw(DimensionMismatch("length(rows) must equal length(cols)"))

  shtmp = SharedArray(Float64, convert(Int, m), pids=pids)
  sh1   = SharedArray(Float64, convert(Int, n), pids=pids)
  sh2   = SharedArray(Float64, convert(Int, n), pids=pids)
  sems  = make_sems(numblocks, pids)
  ps = ParallelSBM(m, n, pids, RemoteRef[], RemoteRef[], numblocks, shtmp, sh1, sh2, sems)
  ranges  = make_lengths(length(rows), weights)
  mblocks = make_blocks(m, convert(Int32, numblocks) )
  nblocks = make_blocks(n, convert(Int32, numblocks) )
  mblock_grid = zeros(Int, numblocks, length(pids))
  nblock_grid = zeros(Int, numblocks, length(pids))
  for i in 1:length(pids)
    sbm = SparseBinMatrix(m, n, rows[ranges[i]], cols[ranges[i]])
    sbm_ref = @spawnat pids[i] fetch(sbm)
    push!(ps.sbms, sbm_ref)
    mblock_grid[:,i] = [! isempty(intersect(sbm.mrange, i)) for i in mblocks]
    nblock_grid[:,i] = [! isempty(intersect(sbm.nrange, i)) for i in nblocks]
  end
  mblock_counts = vec(sum(mblock_grid, 2))
  nblock_counts = vec(sum(nblock_grid, 2))
  for i in 1:length(pids)
    mb = block_order(mblock_counts, find(mblock_grid[:,i]))
    nb = block_order(nblock_counts, find(nblock_grid[:,i]))
    pl_ref = @spawnat pids[i] ParallelLogic(mblocks, nblocks, mb, nb, zeros(m), zeros(n), ps.tmp, ps.sems )
    push!(ps.logic, pl_ref)
  end
  return ps
end

## copies ParallelSBM to new pids
function copyto(F::ParallelSBM, pids::Vector{Int})
  length(pids) == length(F.pids) || throw(DimensionMismatch("length(pids)=$(length(pids)) must equal length(F.pids)=$(length(F.pids))"))

  shtmp = SharedArray(Float64, size(F, 1), pids=pids)
  sh1   = SharedArray(Float64, size(F, 2), pids=pids)
  sh2   = SharedArray(Float64, size(F, 2), pids=pids)
  sems  = make_sems(F.numblocks, pids)
  ps    = ParallelSBM(F.m, F.n, pids, RemoteRef[], RemoteRef[], F.numblocks, shtmp, sh1, sh2, sems)

  for i in 1:length(F.sbms)
    push!(ps.sbms, @spawnat(pids[i], fetch(F.sbms[i])) )

    l = fetch(@spawnat F.pids[i] nonshared(fetch(F.logic[i])))
    logic_ref = @spawnat pids[i] ParallelLogic(l.mblocks, l.nblocks, l.mblock_order, l.nblock_order, zeros(F.m), zeros(F.n), ps.tmp, ps.sems)
    push!(ps.logic, logic_ref)
  end
  return ps
end

sem_init(x::SharedArray)    = ccall(:sem_init, Cint, (Ptr{Void}, Cint, Cuint), x, 1, one(UInt32))
sem_wait(x::SharedArray)    = ccall(:sem_wait, Cint, (Ptr{Void},), x)
sem_trywait(x::SharedArray) = ccall(:sem_trywait, Cint, (Ptr{Void},), x)
sem_post(x::SharedArray)    = ccall(:sem_post, Cint, (Ptr{Void},), x)

gmean(x) = prod(x) ^ (1 / length(x))
pretty(x) = "[" * join([@sprintf("%.3f", i) for i in x], ", ") * "]"

function balanced_parallelsbm(rows::Vector{Int32}, cols::Vector{Int32}, pids::Vector{Int}; numblocks=length(pids)*2, niter=4, ntotal=30, keeplast=4, verbose=false)
  weights = ones(length(pids))
  y = SharedArray(Float64, convert(Int, maximum(rows)) )
  x = SharedArray(Float64, convert(Int, maximum(cols)) )
  local psbm
  for i = 1:niter
    psbm   = ParallelSBM(rows, cols, pids, numblocks=numblocks, weights=weights)
    times  = A_mul_B!_time(y, psbm, x, ntotal)
    ctimes = vec(mean(times[:,end-keeplast:end], 2))
    meantime  = gmean(ctimes)
    weights .*= (meantime ./ ctimes) .^ (1/(1+0.2*i))
    weights   = weights ./ sum(weights)
    verbose && println("$i. ctimes  = ", pretty(ctimes) )
    verbose && println("$i. weights = ", pretty(weights) )
  end
  return psbm
end

function make_blocks(n::Int32, nblocks::Int32)
  bsize   = 8 * ceil(n / nblocks / 8)
  if (bsize-1) * nblocks > n
    bsize = ceil(n / nblocks)
  end
  return  UnitRange{Int32}[convert(Int32, 1+(i-1)*bsize) : convert(Int32, min(n, i*bsize))
                           for i in 1:nblocks]
end

function make_lengths(total::Int, weights)
    wnorm  = weights ./ sum(weights)
    n      = convert(Vector{Int64}, round(wnorm * total))
    excess = sum(n) - total
    i = 0
    while excess != 0
        n[i+1] -= sign(excess)
        excess  = sum(n) - total
        i = (i + 1) % length(n)
    end
    ranges = UnitRange{Int}[]
    k = 1
    for len in n
      push!(ranges, k:(k+len-1))
      k += len
    end
    return ranges
end

## waits 100x to see that x[n] == value
function busywait(x::SharedArray{Int,1}, n::Int, value::Int, ntimes=100)
  for i=1:ntimes
    if x[n] != value
      return false
    end
  end
  return true
end

import Base.size
import Base.isempty
size(X::SparseBinMatrix) = (X.m, X.n)
size(X::SparseBinMatrix, d::Int) = d==1 ? X.m : X.n
isempty(X::SparseBinMatrix)      = X.m == 0 || X.n == 0

size(X::ParallelSBM) = (X.m, X.n)
size(X::ParallelSBM, d::Int) = d==1 ? X.m : X.n
isempty(X::ParallelSBM)      = X.m == 0 || X.n == 0

import Base.A_mul_B!
import Base.At_mul_B!
import Base.At_mul_B
import Base.*

## multiplication: y = A * x
*{Tx}(A::SparseBinMatrix, x::AbstractArray{Tx,1}) = (y = zeros(Tx, A.m); A_mul_B!(y, A, x); y)
*{Tx}(A::ParallelSBM,     x::AbstractArray{Tx,1}) = (y = zeros(Tx, A.m); A_mul_B!(y, A, x); y)

function A_mul_B!{Tx}(y::AbstractArray{Tx,1}, A::SparseBinMatrix, x::AbstractArray{Tx,1})
    A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
    A.m == length(y) || throw(DimensionMismatch("A.m=$(A.m) must equal length(y)=$(length(y))"))
    fill!(y, zero(Tx) )
    rows = A.rows
    cols = A.cols
    @inbounds for i = 1:length(rows)
        y[rows[i]] += x[cols[i]]
    end
    return
end

## multiplication: y = A' * x
At_mul_B{Tx}(A::SparseBinMatrix, x::AbstractArray{Tx,1}) = (y = zeros(Tx, A.n); At_mul_B!(y, A, x); y)

function At_mul_B!{Tx}(y::AbstractArray{Tx,1}, A::SparseBinMatrix, x::AbstractArray{Tx,1})
    A.n == length(y) || throw(DimensionMismatch("A.n=$(A.n) must equal length(y)=$(length(y))"))
    A.m == length(x) || throw(DimensionMismatch("A.m=$(A.m) must equal length(x)=$(length(x))"))
    fill!(y, zero(Tx) )
    rows = A.rows
    cols = A.cols
    @inbounds for i = 1:length(rows)
        y[cols[i]] += x[rows[i]]
    end
    return
end

function A_mul_B!{Tx}(y::SharedArray{Tx,1}, A::ParallelSBM, x::SharedArray{Tx,1})
  A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
  A.m == length(y) || throw(DimensionMismatch("A.m=$(A.m) must equal length(y)=$(length(y))"))
  y[1:end] = zero(Tx)
  np = length(A.pids)
  @sync begin
    for p in 1:length(A.pids)
      pid = A.pids[p]
      if pid != myid() || np == 1
        @async begin
          remotecall_wait(pid, partmul_ref, y, A.sbms[p], A.logic[p], x)
        end
      end
    end
  end
  ## done
end

function At_mul_B!{Tx}(y::SharedArray{Tx,1}, A::ParallelSBM, x::SharedArray{Tx,1})
  A.m == length(x) || throw(DimensionMismatch("A.m=$(A.m) must equal length(x)=$(length(x))"))
  A.n == length(y) || throw(DimensionMismatch("A.n=$(A.n) must equal length(y)=$(length(y))"))
  y[1:end] = zero(Tx)
  np = length(A.pids)
  @sync begin
    for p in 1:length(A.pids)
      pid = A.pids[p]
      if pid != myid() || np == 1
        @async begin
          remotecall_wait(pid, partmul_t_ref, y, A.sbms[p], A.logic[p], x)
        end
      end
    end
  end
end

import Base.*

function *{Tx}(A::ParallelSBM, x::AbstractArray{Tx,1})
  A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
  A.sh1[1:end] = x
  A_mul_B!(A.tmp, A, A.sh1)
  y = zeros(Tx, length(A.tmp))
  y[1:end] = A.tmp
  return y
end

function At_mul_B{Tx}(A::ParallelSBM, x::AbstractArray{Tx,1})
  A.m == length(x) || throw(DimensionMismatch("A.m=$(A.m) must equal length(x)=$(length(x))"))
  A.tmp[1:end] = x
  At_mul_B!(A.sh1, A, A.tmp)
  y = zeros(Tx, A.n)
  y[1:end] = A.sh1
  return y
end

function prod_copy!{Tx}(y::SharedArray{Tx,1}, v::Tx, x::SharedArray{Tx,1})
  length(y) == length(x) || throw(DimensionMismatch("length(y)=$(length(y)) must equal length(x)=$(length(x))"))
  @inbounds @simd for i = 1:length(y)
    y[i] = v * x[i]
  end
  return nothing
end

## computes A'A*x + lambda*x
function AtA_mul_B!{Tx}(y::SharedArray{Tx,1}, A::ParallelSBM, x::SharedArray{Tx,1}, lambda::Float64)
  A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
  A.n == length(y) || throw(DimensionMismatch("A.n=$(A.n) must equal length(y)=$(length(y))"))
  tmp = A.tmp
  tmp[1:end] = zero(Tx)
  np = length(A.pids)
  ## doing tmp = A * x
  @sync begin
    for p in 1:length(A.pids)
      pid = A.pids[p]
      if pid != myid() || np == 1
        @async begin
          remotecall_wait(pid, partmul_ref, tmp, A.sbms[p], A.logic[p], x)
        end
      end
    end
    prod_copy!(y, lambda, x)
  end
  ## doing y += A' * tmp
  @sync begin
    for p in 1:length(A.pids)
      pid = A.pids[p]
      if pid != myid() || np == 1
        @async begin
          remotecall_wait(pid, partmul_t_ref, y, A.sbms[p], A.logic[p], tmp)
        end
      end
    end
  end
  return nothing
end

function A_mul_B!_time{Tx}(y::SharedArray{Tx,1}, A::ParallelSBM, x::SharedArray{Tx,1}, ntimes::Int)
  A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
  A.m == length(y) || throw(DimensionMismatch("A.m=$(A.m) must equal length(y)=$(length(y))"))
  ptime = zeros(length(A.pids), ntimes)
  for i = 1:ntimes
    y[1:end] = zero(Tx)
    ## clearing warmup results
    @sync begin
      for p in 1:length(A.pids)
        pid = A.pids[p]
        if pid != myid() || length(A.pids) == 1
          @async begin
            ptime[p, i] += remotecall_fetch(pid, partmul_time, y, A.sbms[p], A.logic[p], x)
          end
        end
      end
    end
  end
  return ptime
end

function copy!{Tx}(to::AbstractArray{Tx,1}, from::AbstractArray{Tx,1}, range)
  @inbounds @simd for i in range
    to[i] = from[i]
  end
  return nothing
end

function add!{Tx}(to::AbstractArray{Tx,1}, from::AbstractArray{Tx,1}, range)
  @inbounds @simd for i in range
    to[i] += from[i]
  end
  return nothing
end

function rangefill!{Tx}(x::AbstractArray{Tx,1}, v::Tx, range)
  @inbounds @simd for i in range
    x[i] = v
  end
  return nothing
end

function partmul_time{Tx}(y::SharedArray{Tx,1}, Aref::RemoteRef, logicref::RemoteRef, x::SharedArray{Tx,1})
  A     = fetch(Aref)::SparseBinMatrix
  logic = fetch(logicref)::ParallelLogic
  tic();
  partmul(y, A, logic, x)
  return toq()
end

function partmul_ref{Tx}(y::SharedArray{Tx,1}, Aref::RemoteRef, logicref::RemoteRef, x::SharedArray{Tx,1})
  A     = fetch(Aref)::SparseBinMatrix
  logic = fetch(logicref)::ParallelLogic
  partmul(y, A, logic, x)
end

function partmul_t_ref{Tx}(y::SharedArray{Tx,1}, Aref::RemoteRef, logicref::RemoteRef, x::SharedArray{Tx,1})
  A     = fetch(Aref)::SparseBinMatrix
  logic = fetch(logicref)::ParallelLogic
  partmul_t(y, A, logic, x)
end

## assumes sizes are correct
function partmul{Tx}(y::SharedArray{Tx,1}, A::SparseBinMatrix, logic::ParallelLogic, x::SharedArray{Tx,1})
  ylocal = logic.localm
  xlocal = logic.localn
  rangefill!(ylocal, zero(Tx), A.mrange)
  copy!(xlocal, x, A.nrange)
  ## standard y = A * x
  rows = A.rows
  cols = A.cols
  @inbounds for i = 1:length(rows)
      ylocal[rows[i]] += xlocal[cols[i]]
  end
  ## adding the result to shared array
  addshared!(y, ylocal, logic.sems, logic.mblocks, logic.mblock_order, A.mrange)
  return nothing
end

## assumes sizes are correct
function partmul_t{Tx}(y::SharedArray{Tx,1}, A::SparseBinMatrix, logic::ParallelLogic, x::SharedArray{Tx,1})
  ylocal = logic.localn
  xlocal = logic.localm
  rangefill!(ylocal, zero(Tx), A.nrange)
  copy!(xlocal, x, A.mrange)
  ## standard y = A' * x
  rows = A.cols
  cols = A.rows
  @inbounds for i = 1:length(rows)
      ylocal[rows[i]] += xlocal[cols[i]]
  end
  ## adding the result to shared array
  addshared!(y, ylocal, logic.sems, logic.nblocks, logic.nblock_order, A.nrange)
  return nothing
end

## does y += x
function addshared!{Tx}(y::SharedArray{Tx,1}, x::AbstractArray{Tx,1}, sems, ranges, order, yrange)
  blocks = copy(order)
  nblocks = length(blocks)
  pid = myid()
  while true
    i = findfirst(blocks)
    i <= 0 && return nothing ## done
    for j = i:nblocks
      block = blocks[j]
      block == 0 && continue
      ## try to get a lock
      if sem_trywait( sems[block] ) < 0
        ## didn't get lock
        continue
      end
      ## copying result to shared array
      add!(y, x, intersect(ranges[blocks[j]], yrange) )
      sem_post( sems[block] )
      blocks[j] = 0
      break
    end
  end
  return nothing
end

######## parallel operations on Frefs
function solve_cg2(Frefs::Vector{RemoteRef}, rhs::Matrix{Float64}, lambda::Float64; tol=1e-6, maxiter=size(rhs,1))
  beta = zeros(size(rhs,1), size(rhs,2))
  D    = size(rhs,2)
  i    = 1
  # function to produce the next work item from the queue.
  # in this case it's just an index.
  nextidx() = (idx=i; i+=1; idx)
  @sync begin
    for ref in Frefs
      @async begin
        while true
          idx = nextidx()
          idx > D && break
          beta[:,idx] = remotecall_fetch(ref.where, cg_AtA_ref, ref, rhs[:,idx], lambda, tol, maxiter)
        end
      end
    end
  end
  return beta
end

function A_mul_B_ref(Fref::RemoteRef, x::AbstractVector{Float64})
  F = fetch(Fref)
  return F * x
end

function At_mul_B_ref(Fref::RemoteRef, x::AbstractVector{Float64})
  F = fetch(Fref)
  return At_mul_B(F, x)
end

## computes y = F * x in parallel (along columns of x)
function Frefs_mul_B(Frefs::Vector{RemoteRef}, x::Matrix{Float64})
  m, n = fetch( @spawnat Frefs[1].where size(fetch(Frefs[1])) )
  n == size(x, 1) || throw(DimensionMismatch("Frefs.n=$(n) must equal length(x)=$(length(x))"))
  y = zeros(m, size(x,2))
  D = size(x, 2)
  i = 1
  nextidx() = (idx=i; i+=1; idx)
  @sync begin
    for ref in Frefs
      @async begin
        while true
          idx = nextidx()
          idx > D && break
          y[:,idx] = remotecall_fetch(ref.where, A_mul_B_ref, ref, x[:,idx])
        end
      end
    end
  end
  return y
end

## computes y = F' * x in parallel (along columns of x)
function Frefs_t_mul_B(Frefs::Vector{RemoteRef}, x::Matrix{Float64})
  m, n = fetch( @spawnat Frefs[1].where size(fetch(Frefs[1])) )
  m == size(x, 1) || throw(DimensionMismatch("Frefs.m=$(m) must equal length(x)=$(length(x))"))
  y = zeros(n, size(x,2))
  D = size(x, 2)
  i = 1
  nextidx() = (idx=i; i+=1; idx)
  @sync begin
    for ref in Frefs
      @async begin
        while true
          idx = nextidx()
          idx > D && break
          y[:,idx] = remotecall_fetch(ref.where, At_mul_B_ref, ref, x[:,idx])
        end
      end
    end
  end
  return y
end

######## Hilbert ordering ########

## convert (x,y) to d
function xy2d(n, x, y)
    local rx::Bool, ry::Bool
    d = 0::Int64
    s = div(n, 2)
    while s > 0
        rx = (x & s) > 0;
        ry = (y & s) > 0;
        d += s * s * ((3 * rx) $ ry);
        x, y = rot(s, x, y, rx, ry);
        s = div(s, 2)
    end
    return d
end

function rot(n, x, y, rx::Bool, ry::Bool)
    if ry == false
        if rx == true
            x = n - 1 - x
            y = n - 1 - y
        end
        return (y, x)
    end
    return (x, y)
end

function sort_hilbert(rows, cols)
  maxrc = max(maximum(rows), maximum(cols))
  n = 2 ^ round(Int, ceil(log2( maxrc )))
  h = zeros(Int, length(rows))
  @inbounds @simd for i = 1:length(h)
    h[i] = xy2d(n, rows[i] - 1, cols[i] - 1)
  end
  hsorted = sortperm(h)
  return rows[hsorted], cols[hsorted]
end


text:
# https://projecteuler.net/problem=1
#If we list all the natural numbers below 10 that are multiples of 3 or 5,
# we get 3, 5, 6 and 9. The sum of these multiples is 23.
#
#Find the sum of all the multiples of 3 or 5 below 1000.

elements = Set()

function add_multiples(number, elements)
    i = number
    while i < 1000
        push!(elements, i)
        i += number
    end
end

add_multiples(3, elements)
add_multiples(5, elements)

println(sum(elements))

max_stars_repo_name: the-gigi/project-euler

max_issues_repo_name: the-gigi/project-euler

max_forks_repo_name: the-gigi/project-euler


text:
function runall(sin = true)
    T = 200.0
    clock = 200.0
    L = 7
    if sin == true
        α = 0.7 #sin
        u = - 3π
        v = + 3π
        X = SinSDE(α)
    else
        K = 2000
        r = 0.1
        β = 0.1
        u = -log(50)/β     # end points in the lamperti transform
        v= -log(1000)/β
        X = LogGrowthSDE(r, K, β)   #end points in the lamperti tranform
    end
    XX = zigzagsampler(X, T, L, u, v, clock)
    burning = 10.0    #burning
    f = clock - 1.0; n = 30
    db = (f-burning)/n
    b =  burning:db:f
    if sin == true
        p = plotmixing(XX, b, T, L, u, v)
        hline!(p, [n*π for n in -5:2:5])
    else
        p = plotmixing(XX, b, T, L, u, v, x -> exp(- x*β))
        hline!(p, [K])
    end
    display(p)
    return XX
end

x = runall(false)
error("STOP HERE")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
# Weights and points for Gauss-Legendre quadrature for a variety of shapes

# The weights and points for Gauss-Legendre quadrature on the interval F[0,1]
# ∑wᵢ= 1, rᵢ∈ F[0, 1]
#
# N that have entries in this function: N = F[ 1, 2, 3, 4, 5, 10, 15, 20].
function gauss_legendre_quadrature(::Type{F}, ::Val{1}) where {F <: AbstractFloat}
    weights = @SVector F[1.0]
    r = @SVector F[0.5]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{2}) where {F <: AbstractFloat}
    weights = @SVector F[0.5, 0.5]
    r = @SVector F[0.21132486540518713, 0.7886751345948129]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{3}) where {F <: AbstractFloat}
    weights = @SVector F[0.2777777777777776,
                        0.4444444444444444,
                        0.2777777777777776]
    r = @SVector F[0.1127016653792583,
                  0.5,
                  0.8872983346207417]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{4}) where {F <: AbstractFloat}
    weights = @SVector F[
                  0.17392742256872684,
                  0.3260725774312732,
                  0.3260725774312732,
                  0.17392742256872684
               ]
    r = @SVector F[
          0.06943184420297371,
          0.33000947820757187,
          0.6699905217924281,
          0.9305681557970262
         ]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{5}) where {F <: AbstractFloat}
    weights = @SVector F[
                  0.1184634425280945,
                  0.23931433524968326,
                  0.28444444444444444,
                  0.23931433524968326,
                  0.1184634425280945
               ]
    r = @SVector F[
          0.04691007703066802,
          0.23076534494715845,
          0.5,
          0.7692346550528415,
          0.9530899229693319,
         ]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{10}) where {F <: AbstractFloat}
    weights = @SVector F[
                  0.03333567215434387,
                  0.07472567457529025,
                  0.1095431812579911,
                  0.1346333596549981,
                  0.14776211235737646,
                  0.14776211235737646,
                  0.1346333596549981,
                  0.1095431812579911,
                  0.07472567457529025,
                  0.03333567215434387
               ]
    r = @SVector F[
          0.013046735741414184,
          0.06746831665550773,
          0.16029521585048778,
          0.2833023029353764,
          0.4255628305091844,
          0.5744371694908156,
          0.7166976970646236,
          0.8397047841495122,
          0.9325316833444923,
          0.9869532642585859
         ]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{15}) where {F <: AbstractFloat}
    weights = @SVector F[
                    0.015376620998058315,
                    0.03518302374405407,
                    0.053579610233586,
                    0.06978533896307715,
                    0.083134602908497,
                    0.0930805000077811,
                    0.0992157426635558,
                    0.10128912096278064,
                    0.0992157426635558,
                    0.0930805000077811,
                    0.083134602908497,
                    0.06978533896307715,
                    0.053579610233586,
                    0.03518302374405407,
                    0.015376620998058315
                ]
    r = @SVector F[
          0.006003740989757311,
          0.031363303799647024,
          0.0758967082947864,
          0.13779113431991497,
          0.21451391369573058,
          0.3029243264612183,
          0.39940295300128276,
          0.5,
          0.6005970469987172,
          0.6970756735387817,
          0.7854860863042694,
          0.862208865680085,
          0.9241032917052137,
          0.968636696200353,
          0.9939962590102427
         ]
    return weights, r
end

function gauss_legendre_quadrature(::Type{F}, ::Val{20}) where {F <: AbstractFloat}
    weights = @SVector F[
                  0.00880700356957606,
                  0.02030071490019347,
                  0.03133602416705453,
                  0.04163837078835238,
                  0.0509650599086202,
                  0.0590972659807592,
                  0.0658443192245883,
                  0.07104805465919105,
                  0.07458649323630186,
                  0.07637669356536295,
                  0.07637669356536295,
                  0.07458649323630186,
                  0.07104805465919105,
                  0.0658443192245883,
                  0.0590972659807592,
                  0.0509650599086202,
                  0.04163837078835238,
                  0.03133602416705453,
                  0.02030071490019347,
                  0.00880700356957606
               ]
    r = @SVector F[
          0.0034357004074525577,
          0.018014036361043095,
          0.04388278587433703,
          0.08044151408889061,
          0.1268340467699246,
          0.1819731596367425,
          0.24456649902458644,
          0.3131469556422902,
          0.38610707442917747,
          0.46173673943325133,
          0.5382632605667487,
          0.6138929255708225,
          0.6868530443577098,
          0.7554335009754136,
          0.8180268403632576,
          0.8731659532300754,
          0.9195584859111094,
          0.956117214125663,
          0.981985963638957,
          0.9965642995925474
         ]
    return weights, r
end

# # The weights and points for Gauss-Legendre quadrature on the parametric unit triangle
# # ∑wᵢ= 1/2, rᵢ∈ F[0, 1], sᵢ∈ F[0, 1], rᵢ + sᵢ ≤ 1
# # N that have entries in this function: N = F[3, 4, 6, 12, 27, 48, 79]
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{1}) where {F <: AbstractFloat}
#     # P1. 0 negative weights, 0 points outside of the triangle
#     w = @SVector F[0.500000000000000]
#
#     r = @SVector F[0.333333333333333]
#
#     s = @SVector F[0.333333333333333]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{3}) where {F <: AbstractFloat}
#     # P2. 0 negative weights, 0 points outside of the triangle
#     w = @SVector F[0.166666666666667,
#                   0.166666666666667,
#                   0.166666666666667]
#
#     r = @SVector F[0.166666666666667,
#                   0.166666666666667,
#                   0.666666666666667]
#
#     s = @SVector F[0.166666666666667,
#                   0.666666666666667,
#                   0.166666666666667]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{4}) where {F <: AbstractFloat}
#     # P3. 1 negative weights, 0 points outside of the triangle
#     w = @SVectorF[-0.281250000000000,
#                   0.260416666666667,
#                   0.260416666666667,
#                   0.260416666666667]
#
#     r = @SVector F[0.333333333333333,
#                   0.200000000000000,
#                   0.200000000000000,
#                   0.600000000000000]
#
#     s = @SVectorF[ 0.333333333333333,
#                   0.200000000000000,
#                   0.600000000000000,
#                   0.200000000000000]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{6}) where {F <: AbstractFloat}
#     # P4. 0 negative weights, 0 points outside of the triangle
#     w = @SVector F[0.111690794839005,
#                   0.111690794839005,
#                   0.111690794839005,
#                   0.054975871827661,
#                   0.054975871827661,
#                   0.054975871827661]
#     r = @SVector F[0.445948490915965,
#                   0.445948490915965,
#                   0.108103018168070,
#                   0.091576213509771,
#                   0.091576213509771,
#                   0.816847572980459]
#
#     s = @SVector F[0.445948490915965,
#                   0.108103018168070,
#                   0.445948490915965,
#                   0.091576213509771,
#                   0.816847572980459,
#                   0.091576213509771]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{12}) where {F <: AbstractFloat}
#     # P6. 0 negative weights, 0 points outside of the triangle
#     w = @SVector F[0.058393137863189,
#           0.058393137863189,
#           0.058393137863189,
#           0.025422453185104,
#           0.025422453185104,
#           0.025422453185104,
#           0.041425537809187,
#           0.041425537809187,
#           0.041425537809187,
#           0.041425537809187,
#           0.041425537809187,
#           0.041425537809187]
#
#     r = @SVector F[0.249286745170910,
#           0.249286745170910,
#           0.501426509658179,
#           0.063089014491502,
#           0.063089014491502,
#           0.873821971016996,
#           0.310352451033785,
#           0.636502499121399,
#           0.053145049844816,
#           0.310352451033785,
#           0.636502499121399,
#           0.053145049844816]
#
#     s = @SVector F[0.249286745170910,
#           0.501426509658179,
#           0.249286745170910,
#           0.063089014491502,
#           0.873821971016996,
#           0.063089014491502,
#           0.636502499121399,
#           0.053145049844816,
#           0.310352451033785,
#           0.053145049844816,
#           0.310352451033785,
#           0.636502499121399]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{27}) where {F <: AbstractFloat}
#     # P11. 0 negative weights, 3 points outside of the triangle
#     w = @SVector F[0.00046350316448,
#           0.00046350316448,
#           0.00046350316448,
#           0.03857476745740,
#           0.03857476745740,
#           0.03857476745740,
#           0.02966148869040,
#           0.02966148869040,
#           0.02966148869040,
#           0.01809227025170,
#           0.01809227025170,
#           0.01809227025170,
#           0.00682986550135,
#           0.00682986550135,
#           0.00682986550135,
#           0.02616855598110,
#           0.02616855598110,
#           0.02616855598110,
#           0.02616855598110,
#           0.02616855598110,
#           0.02616855598110,
#           0.01035382981955,
#           0.01035382981955,
#           0.01035382981955,
#           0.01035382981955,
#           0.01035382981955,
#           0.01035382981955]
#
#     r = @SVector F[+0.5346110482710,
#           -0.0692220965415,
#           +0.5346110482710,
#           +0.3989693029660,
#           +0.2020613940680,
#           +0.3989693029660,
#           +0.2033099004310,
#           +0.5933801991370,
#           +0.2033099004310,
#           +0.1193509122830,
#           +0.7612981754350,
#           +0.1193509122830,
#           +0.0323649481113,
#           +0.9352701037770,
#           +0.0323649481113,
#           +0.5932012134280,
#           +0.0501781383105,
#           +0.3566206482610,
#           +0.0501781383105,
#           +0.3566206482610,
#           +0.5932012134280,
#           +0.8074890031600,
#           +0.0210220165362,
#           +0.1714889803040,
#           +0.0210220165362,
#           +0.1714889803040,
#           +0.8074890031600]
#
#     s = @SVector F[+0.5346110482710,
#           +0.5346110482710,
#           -0.0692220965415,
#           +0.3989693029660,
#           +0.3989693029660,
#           +0.2020613940680,
#           +0.2033099004310,
#           +0.2033099004310,
#           +0.5933801991370,
#           +0.1193509122830,
#           +0.1193509122830,
#           +0.7612981754350,
#           +0.0323649481113,
#           +0.0323649481113,
#           +0.9352701037770,
#           +0.3566206482610,
#           +0.5932012134280,
#           +0.0501781383105,
#           +0.3566206482610,
#           +0.5932012134280,
#           +0.0501781383105,
#           +0.1714889803040,
#           +0.8074890031600,
#           +0.0210220165362,
#           +0.1714889803040,
#           +0.8074890031600,
#           +0.0210220165362]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{48}) where {F <: AbstractFloat}
#     # P15. 0 negative weights, 9 points outside of the triangle
#     w = @SVector F[0.000958437821425,
#           0.000958437821425,
#           0.000958437821425,
#           0.022124513635550,
#           0.022124513635550,
#           0.022124513635550,
#           0.025593274359450,
#           0.025593274359450,
#           0.025593274359450,
#           0.011843867935350,
#           0.011843867935350,
#           0.011843867935350,
#           0.006644887845000,
#           0.006644887845000,
#           0.006644887845000,
#           0.002374458304095,
#           0.002374458304095,
#           0.002374458304095,
#           0.019275036299800,
#           0.019275036299800,
#           0.019275036299800,
#           0.019275036299800,
#           0.019275036299800,
#           0.019275036299800,
#           0.013607907160300,
#           0.013607907160300,
#           0.013607907160300,
#           0.013607907160300,
#           0.013607907160300,
#           0.013607907160300,
#           0.001091038683400,
#           0.001091038683400,
#           0.001091038683400,
#           0.001091038683400,
#           0.001091038683400,
#           0.001091038683400,
#           0.010752659923850,
#           0.010752659923850,
#           0.010752659923850,
#           0.010752659923850,
#           0.010752659923850,
#           0.010752659923850,
#           0.003836971315525,
#           0.003836971315525,
#           0.003836971315525,
#           0.003836971315525,
#           0.003836971315525,
#           0.003836971315525]
#
#     r = @SVector F[+0.5069729168580,
#           -0.0139458337165,
#           +0.5069729168580,
#           +0.4314063542830,
#           +0.1371872914340,
#           +0.4314063542830,
#           +0.2776936448470,
#           +0.4446127103060,
#           +0.2776936448470,
#           +0.1264648910410,
#           +0.7470702179170,
#           +0.1264648910410,
#           +0.0708083859747,
#           +0.8583832280510,
#           +0.0708083859747,
#           +0.0189651702411,
#           +0.9620696595180,
#           +0.0189651702411,
#           +0.6049544668930,
#           +0.1337341619670,
#           +0.2613113711400,
#           +0.1337341619670,
#           +0.2613113711400,
#           +0.6049544668930,
#           +0.5755865555130,
#           +0.0363666773969,
#           +0.3880467670900,
#           +0.0363666773969,
#           +0.3880467670900,
#           +0.5755865555130,
#           +0.7244626630770,
#           -0.0101748831266,
#           +0.2857122200500,
#           -0.0101748831266,
#           +0.2857122200500,
#           +0.7244626630770,
#           +0.7475564660520,
#           +0.0368438698759,
#           +0.2155996640720,
#           +0.0368438698759,
#           +0.2155996640720,
#           +0.7475564660520,
#           +0.8839645740920,
#           +0.0124598093312,
#           +0.1035756165760,
#           +0.0124598093312,
#           +0.1035756165760,
#           +0.8839645740920]
#
#     s = @SVector F[+0.5069729168580,
#           +0.5069729168580,
#           -0.0139458337165,
#           +0.4314063542830,
#           +0.4314063542830,
#           +0.1371872914340,
#           +0.2776936448470,
#           +0.2776936448470,
#           +0.4446127103060,
#           +0.1264648910410,
#           +0.1264648910410,
#           +0.7470702179170,
#           +0.0708083859747,
#           +0.0708083859747,
#           +0.8583832280510,
#           +0.0189651702411,
#           +0.0189651702411,
#           +0.9620696595180,
#           +0.2613113711400,
#           +0.6049544668930,
#           +0.1337341619670,
#           +0.2613113711400,
#           +0.6049544668930,
#           +0.1337341619670,
#           +0.3880467670900,
#           +0.5755865555130,
#           +0.0363666773969,
#           +0.3880467670900,
#           +0.5755865555130,
#           +0.0363666773969,
#           +0.2857122200500,
#           +0.7244626630770,
#           -0.0101748831266,
#           +0.2857122200500,
#           +0.7244626630770,
#           -0.0101748831266,
#           +0.2155996640720,
#           +0.7475564660520,
#           +0.0368438698759,
#           +0.2155996640720,
#           +0.7475564660520,
#           +0.0368438698759,
#           +0.1035756165760,
#           +0.8839645740920,
#           +0.0124598093312,
#           +0.1035756165760,
#           +0.8839645740920,
#           +0.0124598093312]
#     return w, r, s
# end
#
# function gauss_legendre_quadrature(tri::Triangle6_2D{F}, ::Val{79}) where {F <: AbstractFloat}
#     # P20. 3 negative weights, 9 points outside of the triangle
#     w = @SVector F[+0.016528527770800,
#           +0.000433509592831,
#           +0.000433509592831,
#           +0.000433509592831,
#           +0.005830026358200,
#           +0.005830026358200,
#           +0.005830026358200,
#           +0.011438468178200,
#           +0.011438468178200,
#           +0.011438468178200,
#           +0.015224491336950,
#           +0.015224491336950,
#           +0.015224491336950,
#           +0.015312445862700,
#           +0.015312445862700,
#           +0.015312445862700,
#           +0.012184028838400,
#           +0.012184028838400,
#           +0.012184028838400,
#           +0.007998716016000,
#           +0.007998716016000,
#           +0.007998716016000,
#           +0.003849150907800,
#           +0.003849150907800,
#           +0.003849150907800,
#           -0.000316030248744,
#           -0.000316030248744,
#           -0.000316030248744,
#           +0.000875567150595,
#           +0.000875567150595,
#           +0.000875567150595,
#           +0.008232919594800,
#           +0.008232919594800,
#           +0.008232919594800,
#           +0.008232919594800,
#           +0.008232919594800,
#           +0.008232919594800,
#           +0.002419516770245,
#           +0.002419516770245,
#           +0.002419516770245,
#           +0.002419516770245,
#           +0.002419516770245,
#           +0.002419516770245,
#           +0.012902453267350,
#           +0.012902453267350,
#           +0.012902453267350,
#           +0.012902453267350,
#           +0.012902453267350,
#           +0.012902453267350,
#           +0.004235545527220,
#           +0.004235545527220,
#           +0.004235545527220,
#           +0.004235545527220,
#           +0.004235545527220,
#           +0.004235545527220,
#           +0.009177457053150,
#           +0.009177457053150,
#           +0.009177457053150,
#           +0.009177457053150,
#           +0.009177457053150,
#           +0.009177457053150,
#           +0.000352202338954,
#           +0.000352202338954,
#           +0.000352202338954,
#           +0.000352202338954,
#           +0.000352202338954,
#           +0.000352202338954,
#           +0.005056342463750,
#           +0.005056342463750,
#           +0.005056342463750,
#           +0.005056342463750,
#           +0.005056342463750,
#           +0.005056342463750,
#           +0.001786954692975,
#           +0.001786954692975,
#           +0.001786954692975,
#           +0.001786954692975,
#           +0.001786954692975,
#           +0.001786954692975]
#
#     r = @SVector F[+0.3333333333333,
#           +0.5009504643520,
#           -0.0019009287044,
#           +0.5009504643520,
#           +0.4882129579350,
#           +0.0235740841305,
#           +0.4882129579350,
#           +0.4551366869500,
#           +0.0897266360994,
#           +0.4551366869500,
#           +0.4019962593180,
#           +0.1960074813630,
#           +0.4019962593180,
#           +0.2558929097590,
#           +0.4882141804810,
#           +0.2558929097590,
#           +0.1764882559950,
#           +0.6470234880100,
#           +0.1764882559950,
#           +0.1041708553370,
#           +0.7916582893260,
#           +0.1041708553370,
#           +0.0530689638409,
#           +0.8938620723180,
#           +0.0530689638409,
#           +0.0416187151960,
#           +0.9167625696080,
#           +0.0416187151960,
#           +0.0115819214068,
#           +0.9768361571860,
#           +0.0115819214068,
#           +0.6064026461060,
#           +0.0487415836648,
#           +0.3448557702290,
#           +0.0487415836648,
#           +0.3448557702290,
#           +0.6064026461060,
#           +0.6158426144570,
#           +0.0063141159486,
#           +0.3778432695950,
#           +0.0063141159486,
#           +0.3778432695950,
#           +0.6158426144570,
#           +0.5590480003900,
#           +0.1343165205470,
#           +0.3066354790620,
#           +0.1343165205470,
#           +0.3066354790620,
#           +0.5590480003900,
#           +0.7366067432630,
#           +0.0139738939624,
#           +0.2494193627750,
#           +0.0139738939624,
#           +0.2494193627750,
#           +0.7366067432630,
#           +0.7116751422870,
#           +0.0755491329098,
#           +0.2127757248030,
#           +0.0755491329098,
#           +0.2127757248030,
#           +0.7116751422870,
#           +0.8614027171550,
#           -0.0083681532082,
#           +0.1469654360530,
#           -0.0083681532082,
#           +0.1469654360530,
#           +0.8614027171550,
#           +0.8355869579120,
#           +0.0266860632587,
#           +0.1377269788290,
#           +0.0266860632587,
#           +0.1377269788290,
#           +0.8355869579120,
#           +0.9297561715570,
#           +0.0105477192941,
#           +0.0596961091490,
#           +0.0105477192941,
#           +0.0596961091490,
#           +0.9297561715570]
#
#     s = @SVector F[ +0.3333333333333,
#            +0.5009504643520,
#            +0.5009504643520,
#            -0.0019009287044,
#            +0.4882129579350,
#            +0.4882129579350,
#            +0.0235740841305,
#            +0.4551366869500,
#            +0.4551366869500,
#            +0.0897266360994,
#            +0.4019962593180,
#            +0.4019962593180,
#            +0.1960074813630,
#            +0.2558929097590,
#            +0.2558929097590,
#            +0.4882141804810,
#            +0.1764882559950,
#            +0.1764882559950,
#            +0.6470234880100,
#            +0.1041708553370,
#            +0.1041708553370,
#            +0.7916582893260,
#            +0.0530689638409,
#            +0.0530689638409,
#            +0.8938620723180,
#            +0.0416187151960,
#            +0.0416187151960,
#            +0.9167625696080,
#            +0.0115819214068,
#            +0.0115819214068,
#            +0.9768361571860,
#            +0.3448557702290,
#            +0.6064026461060,
#            +0.0487415836648,
#            +0.3448557702290,
#            +0.6064026461060,
#            +0.0487415836648,
#            +0.3778432695950,
#            +0.6158426144570,
#            +0.0063141159486,
#            +0.3778432695950,
#            +0.6158426144570,
#            +0.0063141159486,
#            +0.3066354790620,
#            +0.5590480003900,
#            +0.1343165205470,
#            +0.3066354790620,
#            +0.5590480003900,
#            +0.1343165205470,
#            +0.2494193627750,
#            +0.7366067432630,
#            +0.0139738939624,
#            +0.2494193627750,
#            +0.7366067432630,
#            +0.0139738939624,
#            +0.2127757248030,
#            +0.7116751422870,
#            +0.0755491329098,
#            +0.2127757248030,
#            +0.7116751422870,
#            +0.0755491329098,
#            +0.1469654360530,
#            +0.8614027171550,
#            -0.0083681532082,
#            +0.1469654360530,
#            +0.8614027171550,
#            -0.0083681532082,
#            +0.1377269788290,
#            +0.8355869579120,
#            +0.0266860632587,
#            +0.1377269788290,
#            +0.8355869579120,
#            +0.0266860632587,
#            +0.0596961091490,
#            +0.9297561715570,
#            +0.0105477192941,
#            +0.0596961091490,
#            +0.9297561715570,
#            +0.0105477192941]
#     return w, r, s
# end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# Just a simple lattice for examples and test.
# We note that there is package Lattices.jl but quite out-dated and poor-maintained. (does not support cubic lattices.)
export OBC, Cubic, graph, locations, Lattice, edgeline

using Reexport
@reexport using LightGraphs, GeometryBasics
abstract type Lattice end
abstract type AbstractBC end
abstract type OBC <: AbstractBC end

"""
    A Lattice has two components.
    G: a Graph including vertices and edges.
    L: Locations as Point for each vertices.
"""
struct Cubic{bc<:AbstractBC} <: Lattice
    graph::AbstractGraph
    locations::Array
end

graph(l::Lattice) = l.graph
locations(l::Lattice) = l.locations

function Cubic{OBC}(Lx::Int, Ly::Int, Lz::Int)
    dim = [Lx; Ly; Lz]
    N = prod(dim)
    D = 3
    g = SimpleGraph(N)
    locations = zeros(Float64, N, 3)
    f(i, j, k) = i + (j - 1) * Lx + (k - 1) * Lx * Ly

    function tadd!(g, pos, dir)
        if pos[dir] != dim[dir]
            tpos = [i == dir ? 1 : 0 for i = 1:D]
            add_edge!(g, f(pos...), f((pos + tpos)...))
        end
    end

    for i = 1:Lx
        for j = 1:Ly
            for k = 1:Lz
                for dir = 1:D
                    tadd!(g, [i; j; k], dir)
                end
                locations[f(i, j, k), :] = [i; j; k]
            end
        end
    end
    return Cubic{OBC}(g, locations)
end

"""
    Function edgeline
    create a Line from an edge
"""
edgeline(l2D::Array, edge) = Line(Point(l2D[src(edge), :]...), Point(l2D[dst(edge), :]...))

# """
#     obtain raw data
# """
# function raw(lattice::Lattice)

# end

max_stars_repo_name: qiyang-ustc/FTensorKit

max_issues_repo_name: qiyang-ustc/FTensorKit

max_forks_repo_name: qiyang-ustc/FTensorKit


text:
"
    Abstract environment for an infinite state
"
abstract type AbstractInfEnv <: Cache end;

function leftenv(envs::AbstractInfEnv,pos::Int,state)
    check_recalculate!(envs,state);
    envs.lw[pos,:]
end

function rightenv(envs::AbstractInfEnv,pos::Int,state)
    check_recalculate!(envs,state);
    envs.rw[pos,:]
end

function leftenv(envs::AbstractInfEnv,row::Int,col::Int,state)
    check_recalculate!(envs,state);
    envs.lw[row,col]
end

function rightenv(envs::AbstractInfEnv,row::Int,col::Int,state)
    check_recalculate!(envs,state);
    envs.rw[row,col]
end

leftenv(envs,pos::CartesianIndex,state) = leftenv(envs,Tuple(pos)...,state)
rightenv(envs,pos::CartesianIndex,state) = rightenv(envs,Tuple(pos)...,state)

function check_recalculate!(envs,state)
    if !(envs.dependency === state)
        #acquire the lock
        lock(envs) do
            if !(envs.dependency === state)
                recalculate!(envs,state);
            end
        end
    end

    return envs;
end

Base.lock(fun::Function,env::AbstractInfEnv) = lock(fun,env.lock)
Base.lock(env::AbstractInfEnv) = lock(env.lock);
Base.unlock(env::AbstractInfEnv) = unlock(env.lock);


text:
using Test

m = page_model()
include("../src/components/N2Oemissions.jl")

add_comp!(m, n2oemissions)

set_param!(m, :n2oemissions, :e0_baselineN2Oemissions, readpagedata(m,"data/e0_baselineN2Oemissions.csv"))
set_param!(m, :n2oemissions, :er_N2Oemissionsgrowth, readpagedata(m, "data/er_N2Oemissionsgrowth.csv"))

##running Model
run(m)

# Generated data
emissions= m[:n2oemissions,  :e_regionalN2Oemissions]
# Recorded data
emissions_compare=readpagedata(m, "test/validationdata/e_regionalN2Oemissions.csv")

@test emissions ≈ emissions_compare rtol=1e-3


text:
module IntegralUtils

using LinearAlgebra
using Decimals

export parse_function, parse_num, Φ, transform, divergence, round_float

"""
    ∂(f::Function, var::Symbol, P₀::Array{T, 1}; Δ::Real = 1e-3)::Union{Real, Array{Real, 1}} where T <: Real

Compute a partial derivative of `f` with respect to `var` at a point `P₀`.

# Examples
```
julia> ∂((x, y, z) -> cos(x)sin(y)exp(z), :y, [π, 0, log(2)])
-1.9999999999996667

julia> ∂((x, y, z) -> cos(x)sin(y)exp(z), :y, [π, 0, log(2)]; Δ = 1e-9)
-2.0

julia> ∂((x, y, z) -> [x, y, z], :z, [1, 1, 1])
3-element Array{Float64,1}:
 0.0
 0.0
 0.9999999999177334

julia> ∂((u, v) -> [cos(u), sin(u), v], :u, [0, 1])
3-element Array{Float64,1}:
 -5.000444502911705e-7
  0.9999999999998334
  0.0
```
"""
function ∂(f::Function, var::Symbol, P₀::Array{T1, 1}; Δ::Real = 1e-6)::Union{Real, Array{T2, 1} where T2 <: Real} where T1 <: Real
    if length(P₀) == 3
        limits = Dict(sym=>vec for (sym, vec) in
                zip((:x, :y, :z), [Δ*Diagonal(ones((3, 3)))[3i-2:3i] for i in 1:3]))
    elseif length(P₀) == 2
        limits = Dict(sym=>vec for (sym, vec) in
                zip((:u, :v), [Δ*Diagonal(ones((2, 2)))[2i-1:2i] for i in 1:2]))
    end
    return (f((P₀ + limits[var])...) - f(P₀...)) / Δ
end

"""
    divergence(F::Function, P₀::Array{T, 1})::Real where T <: Real

Calculate the divergence of a vector field `F` at a point `P₀`.

# Examples
```
julia> divergence((x, y, z) -> [z * cos(x), z * sin(y), z - cos(x+y)], [π, π, 1])
4.998225722196992e-7

julia> divergence((x, y, z) -> [z * cos(x), z * sin(y), z - cos(x+y)], [0, 0, 1])
1.9999994998731165
```
"""
function divergence(F::Function, P₀::Array{T, 1})::Real where T <: Real
   return sum(hcat([∂(F, sym, P₀) for sym in (:x, :y, :z)]...)' .* Diagonal(ones(3, 3)))
end

"""
    normal(r::Function, P₀::Array{T, 1})::Array{Real, 1} where T <: Real

Compute a normal vector to the surface parametrized by `r` at point `P₀`.

# Examples
```
julia> normal((u, v) -> [cos(u)cos(v), sin(u)cos(v), sin(v)], [π/2, 0])
3-element Array{Real,1}:
 -5.000444502910872e-7
  0.9999999999174001
  5.000444502499502e-7

julia> normal((u, v) -> [cos(u), sin(u), v], [π/2, 1])
3-element Array{Real,1}:
 -5.000444502500335e-7
 0.9999999998353001
 0.0
```
"""
function normal(r::Function, P₀::Array{T, 1})::Array{Real, 1} where T <: Real
    rᵤ = ∂(r, :u, P₀)
    rᵥ = ∂(r, :v, P₀)
    return cross(rᵤ, rᵥ)
end

"""
    transform(F::Function, r::Function, P₀::Array{T, 1})::Real where T <: Real

Convert ``F(P₀)⋅dS`` to ``F(r(P₀))dudv``.

# Examples
```
julia> transform((x, y, z) -> [x, y, z], (u, v) -> [cos(u)cos(v), sin(u)cos(v), sin(v)], [π/2, 0])
0.9999999999174001

julia> transform((x, y, z) -> [-y, x, 0], (u, v) -> [(cos(u) + 2)cos(v), (cos(u) + 2)sin(v), sin(u)], [0, 0])
-4.500400052619784e-6
```
"""
function transform(F::Function, r::Function, P₀::Array{T, 1})::Real where T <: Real
   return dot(F(r(P₀...)...), normal(r, P₀))
end

function try_eval(f::Function, args::Float64...)::Float64
    try
        return f(args...)
    catch e
        return 0.0
    end
end

"""
    create_weights(X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function, ξ::Int64, υ::Int64, ζ::Int64)::Array{Float64, 3}

Build a 3d matrix according to composite Simpson's rule for triple integrals.

# Examples
```
julia> create_weights((0, 1), (x, y) -> 0, (x, y) -> 1 - y - x, x -> 0, x -> 1 - x, 2, 2, 2)
3×3×3 Array{Float64,3}:
[:, :, 1] =
 1.0  1.0  0.0
 2.0  2.0  0.0
 0.0  0.0  0.0

[:, :, 2] =
 4.0  4.0  0.0
 8.0  8.0  0.0
 0.0  0.0  0.0

[:, :, 3] =
 1.0  1.0  0.0
 2.0  2.0  0.0
 0.0  0.0  0.0

 julia> create_weights((-1, 1), (x, y) -> sqrt(2-y^2-x^2), (x, y) -> 1, x -> 0, x -> sqrt(1-x^2), 4, 4, 2)
 5×5×3 Array{Float64,3}:
 [:, :, 1] =
  0.0  -1.11847  -0.828427  -1.11847  0.0
  0.0  -4.22673  -3.13553   -4.22673  0.0
  0.0  -1.73205  -1.2915    -1.73205  0.0
  0.0  -2.11231  -1.59166   -2.11231  0.0
  0.0   0.0       0.0        0.0      0.0

 [:, :, 2] =
  0.0   -4.4739    -3.31371   -4.4739   0.0
  0.0  -16.9069   -12.5421   -16.9069   0.0
  0.0   -6.9282    -5.16601   -6.9282   0.0
  0.0   -8.44925   -6.36665   -8.44925  0.0
  0.0    0.0        0.0        0.0      0.0

 [:, :, 3] =
  0.0  -1.11847  -0.828427  -1.11847  0.0
  0.0  -4.22673  -3.13553   -4.22673  0.0
  0.0  -1.73205  -1.2915    -1.73205  0.0
  0.0  -2.11231  -1.59166   -2.11231  0.0
  0.0   0.0       0.0        0.0      0.0
```
"""
function create_weights(X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function,
    ξ::Int64, υ::Int64, ζ::Int64)::Array{Float64, 3}
    Xᵢ = LinRange(X... , ξ + 1)
    areas = reshape([try_eval((a, b) -> η(a, b) - ρ(a, b), x, y) * try_eval(a -> ψ(a) - ϕ(a), x) for x in Xᵢ for y in LinRange(ϕ(x), ψ(x), υ + 1)], (1, υ + 1, ξ + 1))
    areas[areas .< 0] .= 0.0
    return vcat([1], 4 * ones(ξ - 2) - repeat([0, 2], (ξ - 2) ÷ 2), [4, 1])  .*
            vcat([1], 4 * ones(υ - 2) - repeat([0, 2], (υ - 2) ÷ 2), [4, 1])' .*
            reshape(vcat([1], 4 * ones(ζ - 2) - repeat([0, 2], (ζ - 2) ÷ 2), [4, 1]), (1, 1, ζ + 1)) .* areas
end

"""
    create_weights(U::Tuple{Real, Real}, ϕ::Function, ψ::Function, μ::Int64, ν::Int64)::Array{Float64, 2}

Create a 2d matrix of weights according to composite simpson rule for double integrals.

# Examples
```
julia> create_weights((0, 1), x -> 0, x -> x, 6, 4)
7×5 Array{Float64,2}:
 0.0        0.0      0.0       0.0      0.0
 0.666667   2.66667  1.33333   2.66667  0.666667
 0.666667   2.66667  1.33333   2.66667  0.666667
 2.0        8.0      4.0       8.0      2.0
 1.33333    5.33333  2.66667   5.33333  1.33333
 3.33333   13.3333   6.66667  13.3333   3.33333
 1.0        4.0      2.0       4.0      1.0

julia> create_weights((0, π/4), x -> sin(x), x -> cos(x), 8, 4)
9×5 Array{Float64,2}:
 1.0           4.0          2.0           4.0          1.0
 3.58867      14.3547       7.17734      14.3547       3.58867
 1.57139       6.28556      3.14278       6.28556      1.57139
 2.66662      10.6665       5.33325      10.6665       2.66662
 1.08239       4.32957      2.16478       4.32957      1.08239
 1.6421        6.56839      3.2842        6.56839      1.6421
 0.551799      2.2072       1.1036        2.2072       0.551799
 0.554469      2.21787      1.10894       2.21787      0.554469
 1.11022e-16   4.44089e-16  2.22045e-16   4.44089e-16  1.11022e-16
```
"""
function create_weights(U::Tuple{Real, Real}, ϕ::Function, ψ::Function, μ::Int64, ν::Int64)::Array{Float64, 2}
    Uᵢ = LinRange(U..., μ + 1)
    lengths = ψ.(Uᵢ) - ϕ.(Uᵢ)
    return vcat([1], 4 * ones(μ - 2) - repeat([0, 2], (μ - 2) ÷ 2), [4, 1]) .*
            vcat([1], 4 * ones(ν - 2) - repeat([0, 2], (ν - 2) ÷ 2), [4, 1])' .* lengths
end

"""
    split_region(X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function, ξ::Int64, υ::Int64, ζ::Int64)::Array{Array{Float64, 1}, 3}

Divide a 3d region Int64o points on space.

# Examples
```
julia> split_region((0, 1), (x, y) -> 0, (x, y) -> 1, x -> 0, x -> 1, 4, 2, 2)
5×3×3 Array{Array{Float64,1},3}:
[:, :, 1] =
 [0.0, 0.0, 0.0]  [0.0, 0.5, 1.0]   [0.25, 0.0, 0.5]
 [0.0, 0.0, 0.5]  [0.0, 1.0, 0.0]   [0.25, 0.0, 1.0]
 [0.0, 0.0, 1.0]  [0.0, 1.0, 0.5]   [0.25, 0.5, 0.0]
 [0.0, 0.5, 0.0]  [0.0, 1.0, 1.0]   [0.25, 0.5, 0.5]
 [0.0, 0.5, 0.5]  [0.25, 0.0, 0.0]  [0.25, 0.5, 1.0]

[:, :, 2] =
 [0.25, 1.0, 0.0]  [0.5, 0.0, 1.0]  [0.5, 1.0, 0.5]
 [0.25, 1.0, 0.5]  [0.5, 0.5, 0.0]  [0.5, 1.0, 1.0]
 [0.25, 1.0, 1.0]  [0.5, 0.5, 0.5]  [0.75, 0.0, 0.0]
 [0.5, 0.0, 0.0]   [0.5, 0.5, 1.0]  [0.75, 0.0, 0.5]
 [0.5, 0.0, 0.5]   [0.5, 1.0, 0.0]  [0.75, 0.0, 1.0]

[:, :, 3] =
 [0.75, 0.5, 0.0]  [0.75, 1.0, 1.0]  [1.0, 0.5, 0.5]
 [0.75, 0.5, 0.5]  [1.0, 0.0, 0.0]   [1.0, 0.5, 1.0]
 [0.75, 0.5, 1.0]  [1.0, 0.0, 0.5]   [1.0, 1.0, 0.0]
 [0.75, 1.0, 0.0]  [1.0, 0.0, 1.0]   [1.0, 1.0, 0.5]
 [0.75, 1.0, 0.5]  [1.0, 0.5, 0.0]   [1.0, 1.0, 1.0]

julia> split_region((-1, 1), (x, y) -> 0, (x, y) -> sqrt(2-x^2-y^2), x -> -sqrt(1-x^2), x -> sqrt(1-x^2), 4, 4, 2)
5×5×3 Array{Array{Float64,1},3}:
[:, :, 1] =
 [-1.0, 0.0, 0.0]  [-1.0, 0.0, 1.0]  [-1.0, 0.0, 0.5]  [-0.5, -0.866025, 0.0]    [-0.5, -0.433013, 1.25]
 [-1.0, 0.0, 0.5]  [-1.0, 0.0, 0.0]  [-1.0, 0.0, 1.0]  [-0.5, -0.866025, 0.5]    [-0.5, 0.0, 0.0]
 [-1.0, 0.0, 1.0]  [-1.0, 0.0, 0.5]  [-1.0, 0.0, 0.0]  [-0.5, -0.866025, 1.0]    [-0.5, 0.0, 0.661438]
 [-1.0, 0.0, 0.0]  [-1.0, 0.0, 1.0]  [-1.0, 0.0, 0.5]  [-0.5, -0.433013, 0.0]    [-0.5, 0.0, 1.32288]
 [-1.0, 0.0, 0.5]  [-1.0, 0.0, 0.0]  [-1.0, 0.0, 1.0]  [-0.5, -0.433013, 0.625]  [-0.5, 0.433013, 0.0]

[:, :, 2] =
 [-0.5, 0.433013, 0.625]  [0.0, -1.0, 0.0]       [0.0, -0.5, 1.32288]  [0.0, 0.5, 0.661438]  [0.5, -0.866025, 0.0]
 [-0.5, 0.433013, 1.25]   [0.0, -1.0, 0.5]       [0.0, 0.0, 0.0]       [0.0, 0.5, 1.32288]   [0.5, -0.866025, 0.5]
 [-0.5, 0.866025, 0.0]    [0.0, -1.0, 1.0]       [0.0, 0.0, 0.707107]  [0.0, 1.0, 0.0]       [0.5, -0.866025, 1.0]
 [-0.5, 0.866025, 0.5]    [0.0, -0.5, 0.0]       [0.0, 0.0, 1.41421]   [0.0, 1.0, 0.5]       [0.5, -0.433013, 0.0]
 [-0.5, 0.866025, 1.0]    [0.0, -0.5, 0.661438]  [0.0, 0.5, 0.0]       [0.0, 1.0, 1.0]       [0.5, -0.433013, 0.625]

[:, :, 3] =
 [0.5, -0.433013, 1.25]  [0.5, 0.433013, 0.625]  [1.0, 0.0, 0.0]  [1.0, 0.0, 1.0]  [1.0, 0.0, 0.5]
 [0.5, 0.0, 0.0]         [0.5, 0.433013, 1.25]   [1.0, 0.0, 0.5]  [1.0, 0.0, 0.0]  [1.0, 0.0, 1.0]
 [0.5, 0.0, 0.661438]    [0.5, 0.866025, 0.0]    [1.0, 0.0, 1.0]  [1.0, 0.0, 0.5]  [1.0, 0.0, 0.0]
 [0.5, 0.0, 1.32288]     [0.5, 0.866025, 0.5]    [1.0, 0.0, 0.0]  [1.0, 0.0, 1.0]  [1.0, 0.0, 0.5]
 [0.5, 0.433013, 0.0]    [0.5, 0.866025, 1.0]    [1.0, 0.0, 0.5]  [1.0, 0.0, 0.0]  [1.0, 0.0, 1.0]
```
"""
function split_region(X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function,
    ξ::Int64, υ::Int64, ζ::Int64)::Array{Array{Float64, 1}, 3}
    Xᵢ = LinRange(X..., ξ + 1)
    return reshape([[x, y, z] for x in Xᵢ for y in LinRange(ϕ(x), ψ(x), υ + 1) for z in LinRange(try_eval(ρ, x, y), try_eval(η, x, y), ζ + 1)], :, υ + 1, ξ + 1)
end

"""
    split_region(Uᵢ::Tuple{Real, Real}, ϕ::Function, ψ::Function, μ::Int64, ν::Int64)::Array{Array{Float64, 1}, 2}

Divide a 2d region Int64o points on plane.

# Examples
```
julia> split_region((0, 1), x -> 0, x -> x, 6, 4)
7×5 Array{Array{Float64,1},2}:
 [0.0, 0.0]             [0.166667, 0.0833333]  [0.333333, 0.333333]  [0.666667, 0.166667]  [0.833333, 0.625]
 [0.0, 0.0]             [0.166667, 0.125]      [0.5, 0.0]            [0.666667, 0.333333]  [0.833333, 0.833333]
 [0.0, 0.0]             [0.166667, 0.166667]   [0.5, 0.125]          [0.666667, 0.5]       [1.0, 0.0]
 [0.0, 0.0]             [0.333333, 0.0]        [0.5, 0.25]           [0.666667, 0.666667]  [1.0, 0.25]
 [0.0, 0.0]             [0.333333, 0.0833333]  [0.5, 0.375]          [0.833333, 0.0]       [1.0, 0.5]
 [0.166667, 0.0]        [0.333333, 0.166667]   [0.5, 0.5]            [0.833333, 0.208333]  [1.0, 0.75]
 [0.166667, 0.0416667]  [0.333333, 0.25]       [0.666667, 0.0]       [0.833333, 0.416667]  [1.0, 1.0]

julia> split_region((0, 2), x -> -sqrt(4-x^2), x -> sqrt(4-x^2), 6, 6)
7×7 Array{Array{Float64,1},2}:
 [0.0, -2.0]       [0.333333, -1.97203]   [0.666667, -1.88562]   [1.0, -1.73205]  [1.33333, -1.49071]   [1.66667, -1.10554]   [2.0, 0.0]
 [0.0, -1.33333]   [0.333333, -1.31468]   [0.666667, -1.25708]   [1.0, -1.1547]   [1.33333, -0.993808]  [1.66667, -0.737028]  [2.0, 0.0]
 [0.0, -0.666667]  [0.333333, -0.657342]  [0.666667, -0.628539]  [1.0, -0.57735]  [1.33333, -0.496904]  [1.66667, -0.368514]  [2.0, 0.0]
 [0.0, 0.0]        [0.333333, 0.0]        [0.666667, 0.0]        [1.0, 0.0]       [1.33333, 0.0]        [1.66667, 0.0]        [2.0, 0.0]
 [0.0, 0.666667]   [0.333333, 0.657342]   [0.666667, 0.628539]   [1.0, 0.57735]   [1.33333, 0.496904]   [1.66667, 0.368514]   [2.0, 0.0]
 [0.0, 1.33333]    [0.333333, 1.31468]    [0.666667, 1.25708]    [1.0, 1.1547]    [1.33333, 0.993808]   [1.66667, 0.737028]   [2.0, 0.0]
 [0.0, 2.0]        [0.333333, 1.97203]    [0.666667, 1.88562]    [1.0, 1.73205]   [1.33333, 1.49071]    [1.66667, 1.10554]    [2.0, 0.0]
```
"""
function split_region(U::Tuple{Real, Real}, ϕ::Function, ψ::Function, μ::Int64, ν::Int64)::Array{Array{Float64, 1}, 2}
    Uᵢ = LinRange(U..., μ + 1)
    return reshape([[u, v] for u in Uᵢ for v in LinRange(ϕ(u), ψ(u), ν + 1)], μ + 1, :)
end

"""
    coeff(A::Tuple{Real, Real}, n::Int64)::Float64

Compute a coefficient for composite Simpson's rule with `n` midpoints and interval `A`.

# Examples
```
julia> coeff((0, 1), 2)
0.16666666666666666

julia> coeff((0, 4), 4)
0.3333333333333333
```
"""
coeff(A::Tuple{Real, Real}, n::Real)::Float64 = (A[2] - A[1]) / (3 * n)

"""
    ∯(F::Function, X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function, ξ::Int64, υ::Int64, ζ::Int64)::Float64

Determine flux of vector field `F` through a 3d region using Gauss-Ostrogradski theorem and Simpson's rule for ``ξ⋅υ⋅ζ`` nodes.

# Examples
```
julia> ∯((x, y, z) -> [x, y, z], (0, 1), (x, y) -> 0, (x, y) -> 1, x -> 0, x -> 1, 2, 2, 2)
3.0000000000163776

julia> ∯((x, y, z) -> [x^2, y^2, z^2], (-1, 1), (x, y) -> 0, (x, y) -> x^2+y^2, x -> -1, x -> 1, 6, 4, 2)
2.5488967009810803
```
"""
function ∯(F::Function, X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function, ξ::Int64, υ::Int64, ζ::Int64)::Float64

    weights = create_weights(X, ρ, η, ϕ, ψ, ξ, υ, ζ)
    points = split_region(X, ρ, η, ϕ, ψ, ξ, υ, ζ)
    return sum(divergence.(F, points) .* weights) * prod([coeff(interval, steps)
            for (interval, steps) in zip((X, (0, 1), (0, 1)), (ξ, υ, ζ))])
end

"""
    ∯(F::Function, r::Function, U::Tuple{Real, Real}, ϕ::Function, ψ::Function, μ::Int64, ν::Int64)::Float64

Determine flux of vector field `F` through a 2d region using Jacobian and Simpson's rule for ``μ⋅ν`` nodes.

# Examples
```
julia> ∯((x, y, z) -> [x, y, z], (u, v) -> [cos(u), sin(u), v], (0, 2π), u -> 0, u -> u, 2, 2)
19.739208804825942

julia> ∯((x, y, z) -> [x^2, y^2, z^2], (u, v) -> [u, v, 4-u-v], (0, 4), u -> 0, u -> 4-u, 8, 10)
64.00000000410644
```
"""
function ∯(F::Function, r::Function, U::Tuple{Real, Real}, ϕ::Function, ψ::Function, μ::Int64, ν::Int64)::Float64
    points = split_region(U, ϕ, ψ, μ, ν)
    weights = create_weights(U, ϕ, ψ, μ, ν)
    return sum(transform.(F, r, points) .* weights') * prod(coeff(interval, steps)
            for (interval, steps) in zip((U, (0, 1)), (μ, ν)))
end

"""
    draw_points(Xᵢ::LinRange{Float64}, Yᵢ::LinRange{Float64}, Zᵢ::LinRange{Float64}, N::Int64)::Array{Array{Float64, 1}, 1}

Choose randomly 3d points from a cuboid.

# Examples
```
julia> draw_points(LinRange(0, 1, 5), LinRange(0, 1, 5), LinRange(0, 1, 5), 11)
11-element Array{Array{Float64,1},1}:
 [0.75, 0.0, 0.5]
 [1.0, 0.25, 1.0]
 [1.0, 0.0, 1.0]
 [1.0, 0.5, 0.25]
 [0.5, 0.25, 1.0]
 [0.5, 0.0, 0.0]
 [0.5, 0.5, 0.75]
 [1.0, 0.75, 1.0]
 [0.5, 1.0, 0.0]
 [0.5, 1.0, 1.0]
 [0.75, 0.75, 1.0]

julia> draw_points(LinRange(-1, 1, 5), LinRange(0, π, 19), LinRange(2, 2, 2), 3)
3-element Array{Array{Float64,1},1}:
 [-0.5, 2.0943951023931953, 2.0]
 [0.0, 2.96705972839036, 2.0]
 [0.0, 2.443460952792061, 2.0]
```
"""
function draw_points(Xᵢ::LinRange{Float64}, Yᵢ::LinRange{Float64}, Zᵢ::LinRange{Float64}, N::Int64)::Array{Array{Float64, 1}, 1}
    return [[rand(Xᵢ), rand(Yᵢ), rand(Zᵢ)] for i in 1:N]
end

"""
    draw_points(Xᵢ::LinRange{Float64}, Yᵢ::LinRange{Float64}, N::Int64)::Array{Array{Float64, 1}, 1}

Choose randomly 2d points from a rectangle.

# Examples
```
julia> draw_points(LinRange(0, 1, 5), LinRange(0, 1, 5), 11)
11-element Array{Array{Float64,1},1}:
 [0.5, 1.0]
 [0.75, 0.75]
 [0.5, 0.75]
 [0.0, 0.25]
 [1.0, 1.0]
 [0.75, 0.25]
 [0.25, 0.25]
 [0.5, 0.75]
 [1.0, 0.5]
 [0.5, 1.0]
 [0.0, 0.25]

julia> draw_points(LinRange(-1, 1, 5), LinRange(2, 2, 19), 3)
3-element Array{Array{Float64,1},1}:
 [0.5, 2.0]
 [-1.0, 2.0]
 [-0.5, 2.0]
```
"""
function draw_points(Xᵢ::LinRange{Float64}, Yᵢ::LinRange{Float64}, N::Int64)::Array{Array{Float64, 1}, 1}
    return [[rand(Xᵢ), rand(Yᵢ)] for i in 1:N]
end

"""
    surface_part(F::Function, r::Function, S::Array{LinRange{Float64}, 1}, N::Int64, ϕ::Function, ψ::Function)::Float64

Compute flux over one of parts of a rectangle using Monte Carlo method.

# Examples
```
julia> surface_part((x, y, z) -> [x, y, z], (u, v) -> [u, v, 0], [LinRange(0, 1, 100), LinRange(0, 1, 100)], 1000, u -> 0, u -> u)
0.0

julia> surface_part((x, y, z) -> [x^2, y^2, z^2], (u, v) -> [u, v, u-v], [LinRange(0, 1, 100), LinRange(-1, 1, 100)], 100000, u -> -u, u -> u)
0.33893245180194237
```
"""
function surface_part(F::Function, r::Function, S::Array{LinRange{Float64}, 1}, N::Int64, ϕ::Function, ψ::Function)::Float64
    Uᵢ, Vᵢ = S
    points = draw_points(Uᵢ, Vᵢ, N)
    validate(P::Array{Float64, 1}) = ϕ(P[1]) <= P[2] <= ψ(P[1]) ? transform(F, r, P) : 0.0
    return sum(validate.(points)) * (max(Uᵢ...) - min(Uᵢ...)) * (max(Vᵢ...) - min(Vᵢ...)) / N
end

"""
    surface_part(F::Function, S::Array{LinRange{Float64}, 1}, N::Int64, ρ::Function, η::Function, ϕ::Function, ψ::Function)::Float64

Compute flux over one of parts of a cuboid using Monte Carlo method.

# Examples
```
julia> surface_part((x, y, z) -> [x, y, z], [LinRange(0, 1, 100), LinRange(0, 1, 100), LinRange(0, 1, 100)], 1000, (x, y) -> 0, (x, y) -> 1, x -> 0, x -> x)
1.4940000000104536

julia> surface_part((x, y, z) -> [x^2, y^2, z^2], [LinRange(0, 1, 100), LinRange(0, 1, 100), LinRange(-1, 1, 100)], 100000, (x, y) -> 0, (x, y) -> x+y, x -> -x, x -> x)
1.301682061621842
```
"""
function surface_part(F::Function, S::Array{LinRange{Float64}, 1}, N::Int64, ρ::Function, η::Function, ϕ::Function, ψ::Function)::Float64
    Xᵢ, Yᵢ, Zᵢ = S
    points = draw_points(Xᵢ, Yᵢ, Zᵢ, N)
    validate(P::Array{Float64, 1}) = ϕ(P[1]) <= P[2] <= ψ(P[1]) && ρ(P[1:2]...) <= P[3] <= η(P[1:2]...) ? divergence(F, P) : 0.0
    return sum(validate.(points)) * (max(Xᵢ...) - min(Xᵢ...)) * (max(Yᵢ...) - min(Yᵢ...)) * (max(Zᵢ...) - min(Zᵢ...)) / N
end

"""
    ∯(F::Function, r::Function, U::Tuple{Real, Real}, ϕ::Function, ψ::Function, N::Int64, parts::Float64)::Float64

Determine flux of vector field `F` through a 3d region using Gauss-Ostrogradski theorem and Monte Carlo method for `N` random points.

# Examples
```
julia> ∯((x, y, z) -> [x, y, z], (0, 1), (x, y) -> 0, (x, y) -> 1, x -> 0, x -> 1, 10000, 16.0)
2.999400020649781

julia> ∯((x, y, z) -> [x^2, y^2, z^2], (-1, 1), (x, y) -> 0, (x, y) -> x^2+y^2, x -> -1, x -> 1, 100, 256.0)
2.4227719867790514
```
"""
function ∯(F::Function, r::Function, U::Tuple{Real, Real}, ϕ::Function, ψ::Function, N::Int64, parts::Float64)::Float64
    x_part = Int64(parts / 2)
    y_part = Int64(parts / x_part)
    Uᵢ = LinRange(U..., x_part*N)
    Vᵢ = LinRange(min(ϕ.(Uᵢ)...), max(ψ.(Uᵢ)...), y_part*N)
    S = [[Uᵢ[N*i+1:(i+1)N], Vᵢ[N*j+1:(j+1)N]] for i in 0:(x_part-1), j in 0:(y_part-1)]
    return sum(surface_part.(F, r, S, N, ϕ, ψ))
end

"""
    ∯(F::Function, X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function, N::Int64, parts::Float64)::Float64

Determine flux of vector field `F` through a 2d region using Jacobian and Monte Carlo method for `N` random points.

# Examples
```
julia> ∯((x, y, z) -> [x, y, z], (u, v) -> [cos(u), sin(u), v], (0, 2π), u -> 0, u -> u, 10000, 16.0)
19.776707787586332

julia> ∯((x, y, z) -> [x^2, y^2, z^2], (u, v) -> [u, v, 4-u-v], (0, 4), u -> 0, u -> 4-u, 100000, 32.0)
63.9916073475161
```
"""
function ∯(F::Function, X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function, N::Int64, parts::Float64)::Float64
    x_part = y_part = Int64(parts / 4)
    z_part = Int64(parts / 2x_part)
    Xᵢ = LinRange(X..., x_part*N)
    Yᵢ = LinRange(min(ϕ.(Xᵢ)...), max(ψ.(Xᵢ)...), y_part*N)
    unzip(f, P) = try_eval(f, P...)
    Zᵢ = LinRange(min(unzip.(ρ, zip(Xᵢ, Yᵢ))...), max(unzip.(η, zip(Xᵢ, Yᵢ))...), z_part*N)
    S = [[Xᵢ[N*i+1:(i+1)N], Yᵢ[N*j+1:(j+1)N], Zᵢ[N*k+1:(k+1)N]] for i in 0:(x_part-1), j in 0:(y_part-1), k in 0:(z_part-1)]
    return sum(surface_part.(F, S, N, ρ, η, ϕ, ψ))
end

"""
    round_float(α::Float64, ϵ::Float64)::Union{Float64, Int64}

Round `α` to the same Real of digits as `ϵ` or nearest integer with `ϵ` accuracy.

# Examples
```
julia> round_float(π, 1e-12)
3.14159265359

julia> round_float(exp(1), 1e-3)
2.718

julia> round_float(1.999, 1e-3)
2
```
"""
round_float(α::Float64, ϵ::Float64)::Union{Int64, Float64} = abs(α - round(α)) < ϵ ? Int64(round(α)) : round(α, digits = abs(Decimal(ϵ).q))

"""
    arguments(args::Symbol...)

Parse arguments' symbols to a string.

# Examples
```
julia> arguments(:u, :v)
"(u, v)"

julia> arguments(:x, :y, :z)
"(x, y, z)"
```
"""
arguments(args::Symbol...)::String = '(' * String(args[1]) * reduce(*, ", " .* String.(args[2:end])) * ')'

"""
    parse_function(name::Symbol, body::String, args::Symbol...)::Function

Parse a body of a function and arguments' and function name's symbols to the actual function.

# Examples
```
julia> parse_function(:f, "cos(x)sin(y)+exp(z)", :x, :y, :z)
f (generic function with 1 method)

julia> parse_function(:g, "u^2-cos(exp(1/v))", :u, :v)
g (generic function with 1 method)
```
"""
parse_function(body::String, args::Symbol...)::Function = eval(Meta.parse(arguments(args...) * "->" * body))

"""
    parse_function(name::Symbol, body::String, args::Symbol)::Function

Parse a body of a function and argument and function name's symbols to the actual function.

# Examples
```
julia> parse_function(:f, "2cos(x)", :x)
f (generic function with 1 method)

julia> parse_function(:g, "u^2-7u", :u)
g (generic function with 1 method)
```
"""
parse_function(body::String, args::Symbol)::Function = eval(Meta.parse(string(args) * "->" * body))

"""
    parse_num(value::String)::Float64

Parse a Real or an expression from string to float.

# Examples
```
julia> parse_num("2pi")
6.283185307179586

julia> parse_num("2pi - 6/3 + 7.2")
11.483185307179586
```
"""
parse_num(value::String)::Float64 = eval(Meta.parse(value))

"""
    Φ(F::Function, X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function; N::Int64 = 100, technique::String = "Simpson")::Union{Int64, Float64}

Determine a flux of a vector field `F` through a surface.

# Examples
```
julia> Φ((x, y, z) -> [x, y, z], (0, 1), (x, y) -> 0, (x, y) -> 1, x -> 0, x -> 1)
3

julia> Φ((x, y, z) -> [-y, x, 0], (-1, 1), (x, y) -> 0, (x, y) -> sqrt(2-y^2-x^2), x -> -sqrt(1-x^2), x -> sqrt(1-x^2); technique = "Monte Carlo", ϵ = 0.1)
0
```
"""
function Φ(F::Function, X::Tuple{Real, Real}, ρ::Function, η::Function, ϕ::Function, ψ::Function;
     N::Int64 = 100, technique::String = "Simpson")::Union{Int64, Float64}
    if technique in ("Simpson", "Monte Carlo")
        surface_integral(n) = technique == "Simpson" ? ∯(F, X, ρ, η, ϕ, ψ, 2n, 2n, 2n) : ∯(F, X, ρ, η, ϕ, ψ, 100n, 4.0n)
    else
        throw(ArgumentError("Invalid technique name."))
    end
    return abs(surface_integral(N))
end

"""
    Φ(F::Function, r::Function, U::Tuple{Real, Real}, ϕ::Function, ψ::Function; N::Int64 = 100, technique::String = "Simpson")::Union{Int64, Float64}

Determine a flux of a vector field `F` through a surface parametrized by `r`.

# Examples
```
julia> Φ((x, y, z) -> (x^2, y^2, z^2), (u, v) -> [u, v, 4-u-v], (0, 4), u -> 0, u -> 4-u)
64

julia> Φ((x, y, z) -> [x-y, y-z, z-x], (u, v) -> [(cos(u)+2)cos(v), (cos(u)+2)sin(v), sin(u)], (0, 2π), u -> 0, u -> 2π; technique = "Monte Carlo", ϵ = 0.01)
-118.25
```
"""
function Φ(F::Function, r::Function, U::Tuple{Real, Real}, ϕ::Function, ψ::Function;
     N::Int64 = 100, technique::String = "Simpson")::Union{Int64, Float64}
    if technique in ("Simpson", "Monte Carlo")
        surface_integral(n::Int64) = technique == "Simpson" ? ∯(F, r, U, ϕ, ψ, 2n, 2n) : ∯(F, r, U, ϕ, ψ, 100n, 4.0n)
    else
        throw(ArgumentError("Invalid technique name."))
    end
    return abs(surface_integral(N))
end

end


text:
@defcomp EquityWeighting begin
    region = Index()

    # Basic information
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")

    # Impacts across all gases
    pop_population = Parameter(index=[time, region], unit="million person")

    #Total and Per-Capita Abatement and Adaptation Costs
    tct_percap_totalcosts_total = Parameter(index=[time, region], unit="\$/person")
    act_adaptationcosts_total = Parameter(index=[time, region], unit="\$million")
    act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")

    # Consumption
    cons_percap_consumption = Parameter(index=[time, region], unit="\$/person") # Called "CONS_PER_CAP"
    cons_percap_consumption_0 = Parameter(index=[region], unit="\$/person")
    cons_percap_aftercosts = Parameter(index=[time, region], unit="\$/person")

    # Calculation of weighted costs
    emuc_utilityconvexity = Parameter(unit="none", default=1.1666666667)

    wtct_percap_weightedcosts = Variable(index=[time, region], unit="\$/person")
    eact_percap_weightedadaptationcosts = Variable(index=[time, region], unit="\$/person")
    wact_percap_partiallyweighted = Variable(index=[time, region], unit="\$/person")
    wact_partiallyweighted = Variable(index=[time, region], unit="\$million")

    # Amount of equity weighting variable (0, (0, 1), or 1)
    equity_proportion = Parameter(unit="fraction", default=1.0)

    pct_percap_partiallyweighted = Variable(index=[time, region], unit="\$/person")
    pct_partiallyweighted = Variable(index=[time, region], unit="\$million")
    pct_g_partiallyweighted_global = Variable(index=[time], unit="\$million")

    # Discount rates
    ptp_timepreference = Parameter(unit="%/year", default=1.0333333333) # <0.1,1, 2>
    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year")
    popgrw_populationgrowth = Parameter(index=[time, region], unit="%/year")

    dr_discountrate = Variable(index=[time, region], unit="%/year")
    yp_yearsperiod = Variable(index=[time], unit="year") # defined differently from yagg
    dfc_consumptiondiscountrate = Variable(index=[time, region], unit="1/year")

    df_utilitydiscountrate = Variable(index=[time], unit="fraction")

    # Discounted costs
    pcdt_partiallyweighted_discounted = Variable(index=[time, region], unit="\$million")
    pcdt_g_partiallyweighted_discountedglobal = Variable(index=[time], unit="\$million")

    pcdat_partiallyweighted_discountedaggregated = Variable(index=[time, region], unit="\$million")
    tpc_totalaggregatedcosts = Variable(unit="\$million")

    wacdt_partiallyweighted_discounted = Variable(index=[time, region], unit="\$million")

    # Equity weighted impact totals
    rcons_percap_dis = Parameter(index=[time, region], unit="\$/person")

    wit_equityweightedimpact = Variable(index=[time, region], unit="\$million")
    widt_equityweightedimpact_discounted = Variable(index=[time, region], unit="\$million")

    yagg_periodspan = Parameter(index=[time], unit="year")

    addt_equityweightedimpact_discountedaggregated = Variable(index=[time, region], unit="\$million")
    addt_gt_equityweightedimpact_discountedglobal = Variable(unit="\$million")

    civvalue_civilizationvalue = Parameter(unit="\$million", default=5.3e10) # Called "CIV_VALUE"
    td_totaldiscountedimpacts = Variable(unit="\$million")

    aact_equityweightedadaptation_discountedaggregated = Variable(index=[time, region], unit="\$million")
    tac_totaladaptationcosts = Variable(unit="\$million")

    # Final result: total effect of climate change
    te_totaleffect = Variable(unit="\$million")

    function run_timestep(p, v, d, tt)
        if is_first(tt)
            v.tpc_totalaggregatedcosts = 0
            v.addt_gt_equityweightedimpact_discountedglobal = 0
            v.tac_totaladaptationcosts = 0
            v.te_totaleffect = 0
        end

        v.df_utilitydiscountrate[tt] = (1 + p.ptp_timepreference / 100)^(-(p.y_year[tt] - p.y_year_0))

        for rr in d.region

            ## Gas Costs Accounting
            # Weighted costs (Page 23 of Hope 2009)
            v.wtct_percap_weightedcosts[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_consumption[tt, rr]^(1 - p.emuc_utilityconvexity) - (p.cons_percap_consumption[tt, rr] - p.tct_percap_totalcosts_total[tt, rr])^(1 - p.emuc_utilityconvexity))

            # Add these into consumption
            v.eact_percap_weightedadaptationcosts[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_consumption[tt, rr]^(1 - p.emuc_utilityconvexity) - (p.cons_percap_consumption[tt, rr] - p.act_percap_adaptationcosts[tt, rr])^(1 - p.emuc_utilityconvexity))

            # Do partial weighting
            if p.equity_proportion == 0
                v.pct_percap_partiallyweighted[tt, rr] = p.tct_percap_totalcosts_total[tt, rr]
                v.wact_percap_partiallyweighted[tt, rr] = p.act_percap_adaptationcosts[tt, rr]
            else
                v.pct_percap_partiallyweighted[tt, rr] = (1 - p.equity_proportion) * p.tct_percap_totalcosts_total[tt, rr] + p.equity_proportion * v.wtct_percap_weightedcosts[tt, rr]
                v.wact_percap_partiallyweighted[tt, rr] = (1 - p.equity_proportion) * p.act_percap_adaptationcosts[tt, rr] + p.equity_proportion * v.eact_percap_weightedadaptationcosts[tt, rr]
            end

            v.pct_partiallyweighted[tt, rr] = v.pct_percap_partiallyweighted[tt, rr] * p.pop_population[tt, rr]
            v.wact_partiallyweighted[tt, rr] = v.wact_percap_partiallyweighted[tt, rr] * p.pop_population[tt, rr]

            # Discount rate calculations
            v.dr_discountrate[tt, rr] = p.ptp_timepreference + p.emuc_utilityconvexity * (p.grw_gdpgrowthrate[tt, rr] - p.popgrw_populationgrowth[tt, rr])
            if is_first(tt)
                v.yp_yearsperiod[1] = p.y_year[1] - p.y_year_0
            else
                v.yp_yearsperiod[tt] = p.y_year[tt] - p.y_year[tt-1]
            end

            if is_first(tt)
                v.dfc_consumptiondiscountrate[1, rr] = (1 + v.dr_discountrate[1, rr] / 100)^(-v.yp_yearsperiod[1])
            else
                v.dfc_consumptiondiscountrate[tt, rr] = v.dfc_consumptiondiscountrate[tt - 1, rr] * (1 + v.dr_discountrate[tt, rr] / 100)^(-v.yp_yearsperiod[tt])
            end

            # Discounted costs
            if p.equity_proportion == 0
                v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
                v.wacdt_partiallyweighted_discounted[tt, rr] = p.act_adaptationcosts_total[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
            else
                v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.df_utilitydiscountrate[tt]
                v.wacdt_partiallyweighted_discounted[tt, rr] = v.wact_partiallyweighted[tt, rr] * v.df_utilitydiscountrate[tt]
            end

            v.pcdat_partiallyweighted_discountedaggregated[tt, rr] = v.pcdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]

            ## Equity weighted impacts (end of page 28, Hope 2009)
            v.wit_equityweightedimpact[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_aftercosts[tt, rr]^(1 - p.emuc_utilityconvexity) - p.rcons_percap_dis[tt, rr]^(1 - p.emuc_utilityconvexity)) * p.pop_population[tt, rr]

            v.widt_equityweightedimpact_discounted[tt, rr] = v.wit_equityweightedimpact[tt, rr] * v.df_utilitydiscountrate[tt]

            v.addt_equityweightedimpact_discountedaggregated[tt, rr] = v.widt_equityweightedimpact_discounted[tt, rr] * p.yagg_periodspan[tt]
            v.aact_equityweightedadaptation_discountedaggregated[tt, rr] = v.wacdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]
        end

        v.pct_g_partiallyweighted_global[tt] = sum(v.pct_partiallyweighted[tt, :])
        v.pcdt_g_partiallyweighted_discountedglobal[tt] = sum(v.pcdt_partiallyweighted_discounted[tt, :])
        v.tpc_totalaggregatedcosts = v.tpc_totalaggregatedcosts + sum(v.pcdat_partiallyweighted_discountedaggregated[tt, :])

        v.addt_gt_equityweightedimpact_discountedglobal = v.addt_gt_equityweightedimpact_discountedglobal + sum(v.addt_equityweightedimpact_discountedaggregated[tt, :])

        v.tac_totaladaptationcosts = v.tac_totaladaptationcosts + sum(v.aact_equityweightedadaptation_discountedaggregated[tt, :])

        v.td_totaldiscountedimpacts = min(v.addt_gt_equityweightedimpact_discountedglobal, p.civvalue_civilizationvalue)

        # Total effect of climate change
        v.te_totaleffect = min(v.td_totaldiscountedimpacts + v.tpc_totalaggregatedcosts + v.tac_totaladaptationcosts, p.civvalue_civilizationvalue)
    end
end


text:
module Test09

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp
using Test

@testset "Test09" begin
    tokenizer = Automa.compile(
        re"a"      => :(emit(:a, ts:te)),
        re"a*b"    => :(emit(:ab, ts:te)),
    )
    ctx = Automa.CodeGenContext()

    @eval function tokenize(data)
        $(Automa.generate_init_code(ctx, tokenizer))
        p_end = p_eof = sizeof(data)
        tokens = Tuple{Symbol,String}[]
        emit(kind, range) = push!(tokens, (kind, data[range]))
        while p ≤ p_eof && cs > 0
            $(Automa.generate_exec_code(ctx, tokenizer))
        end
        if cs < 0
            error()
        end
        return tokens
    end

    @test tokenize("") == []
    @test tokenize("a") == [(:a, "a")]
    @test tokenize("b") == [(:ab, "b")]
    @test tokenize("aa") == [(:a, "a"), (:a, "a")]
    @test tokenize("ab") == [(:ab, "ab")]
    @test tokenize("aaa") == [(:a, "a"), (:a, "a"), (:a, "a")]
    @test tokenize("aab") == [(:ab, "aab")]
    @test tokenize("abaabba") == [(:ab, "ab"), (:ab, "aab"), (:ab, "b"), (:a, "a")]
    @test_throws ErrorException tokenize("c")
    @test_throws ErrorException tokenize("ac")
    @test_throws ErrorException tokenize("abc")
    @test_throws ErrorException tokenize("acb")
end

end


text:
module TempleEye

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Temple Eye (Small)" => Ahorn.EntityPlacement(
        Maple.TempleEye
    )
)

function Ahorn.selection(entity::Maple.TempleEye)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 6, y - 6, 12, 12)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TempleEye, room::Maple.Room)
    x, y = Ahorn.position(entity)
    tx, ty = floor(Int, x / 8) + 1, floor(Int, y / 8) + 1

    layer = get(room.fgTiles.data, (ty, tx), '0') == '0' ? "bg" : "fg"

    Ahorn.drawSprite(ctx, "scenery/temple/eye/$(layer)_eye.png", 0, 0)
    Ahorn.drawSprite(ctx, "scenery/temple/eye/$(layer)_lid00.png", 0, 0)
    Ahorn.drawSprite(ctx, "scenery/temple/eye/$(layer)_pupil.png", 0, 0)
end

end