text:
using GeoStatsDevTools
using GeoStatsBase
using DataFrames: DataFrame
using Plots; gr(size=(600,400))
using VisualRegressionTests
using Test, Pkg

# list of maintainers
maintainers = ["juliohm"]

# environment settings
istravis = "TRAVIS" ∈ keys(ENV)
ismaintainer = "USER" ∈ keys(ENV) && ENV["USER"] ∈ maintainers
datadir = joinpath(@__DIR__,"data")

if ismaintainer
  Pkg.add("Gtk")
  using Gtk
end

# load data sets
data1D    = readgeotable(joinpath(datadir,"data1D.tsv"), delim='\t', coordnames=[:x])
data2D    = readgeotable(joinpath(datadir,"data2D.tsv"), delim='\t', coordnames=[:x,:y])
data3D    = readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t')
missdata  = readgeotable(joinpath(datadir,"missing.tsv"), delim='\t', coordnames=[:x,:y])

# list of tests
testfiles = [
  "distances.jl",
  "distributions.jl",
  "spatialdata.jl",
  "domains.jl",
  "paths.jl",
  "neighborhoods.jl",
  "mappers.jl",
  "problems.jl",
  "utils.jl"
]

@testset "GeoStatsDevTools.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
# MathOptInterface NLP

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    return model.nlp_data = nlp_data
end

# Keep loading of NLP constraints apart to load all NLP model all in once
# inside Knitro.
function load_nlp_constraints(model::Optimizer)
    num_nlp_constraints = length(model.nlp_data.constraint_bounds)

    # We need to load the NLP constraints inside Knitro.
    if num_nlp_constraints > 0
        num_cons = KN_add_cons(model.inner, num_nlp_constraints)

        for (ib, pair) in enumerate(model.nlp_data.constraint_bounds)
            if pair.upper == pair.lower
                KN_set_con_eqbnd(model.inner, num_cons[ib], pair.upper)
            else
                ub = check_value(pair.upper)
                lb = check_value(pair.lower)
                KN_set_con_upbnd(model.inner, num_cons[ib], ub)
                KN_set_con_lobnd(model.inner, num_cons[ib], lb)
            end
        end
        # Add constraint to index.
        model.nlp_index_cons = num_cons
    end
    return
end

function MOI.supports(::Optimizer, ::MOI.NLPBlockDualStart)
    return true
end

# TODO: FIXME
function MOI.set(model::Optimizer, ::MOI.NLPBlockDualStart, values)
    # @assert length(values) == length(model.nlp_index_cons)
    KN_set_con_dual_init_values(model.inner, Cint[0], values)
    return
end

ext: jl


text:
export Epilogue
module Epilogue

using CUDA
using GemmKernels
using GemmKernels.Tiling
using KernelAbstractions.Extras: @unroll

# ----------------
# Default epilogue
# ----------------

struct Default end

@inline function (ep::Default)(d, shmem_d, transform, ::Type{conf}) where {conf <: GemmKernels.Config}
    # Constants
    block_i = (blockIdx().x - 1) * conf.block_shape.M
    block_j = (blockIdx().y - 1) * conf.block_shape.N

    warpId = (threadIdx().x - 1) ÷ 32 + 1
    laneId = (threadIdx().x - 1) % 32 + 1

    gemm_sz = Tile(conf.matmul_shape)
    block_tile = Tile(conf.block_shape)

    # Cooperatively store a block_shape.M x block_shape.N tile of D from shared to global memory within one threadblock
    @unroll for warp_tile = parallellise(block_tile.MN, Tile(conf.mem_cd_warp), warpId, conf.warps_per_block)
        @unroll for thread_tile = parallellise(warp_tile, Tile(conf.mem_cd_thread), laneId, 32)
            x = Layout.load(conf.shared_d_layout, shmem_d, thread_tile)
            x = transform(x, thread_tile)
            Layout.store!(conf.global_d_layout, d, x, translate_base(thread_tile, (M = block_i, N = block_j)))
        end
    end
end

# -------------
# Bias epilogue
# -------------

struct Bias{B}
    bias_pointer::B
end

@inline function apply_bias(x, bias_pointer::CuPtr{Float32}, thread_tile)
    dev_ptr = reinterpret(Core.LLVMPtr{Float32, AS.Global}, bias_pointer)

    # Load bias value for this column
    col = thread_tile.index.N + 1
    b = unsafe_load(dev_ptr, col)

    return ntuple(k -> VecElement{Float32}(x[k].value + b), Val(4))
end

@inline function (ep::Bias{B})(d, shmem_d, transform, ::Type{conf}) where {B, conf <: GemmKernels.Config}
    # Constants
    block_i = (blockIdx().x - 1) * conf.block_shape.M
    block_j = (blockIdx().y - 1) * conf.block_shape.N

    warpId = (threadIdx().x - 1) ÷ 32 + 1
    laneId = (threadIdx().x - 1) % 32 + 1

    gemm_sz = Tile(conf.matmul_shape)
    block_tile = Tile(conf.block_shape)

    # Cooperatively store a block_shape.M x block_shape.N tile of D from shared to global memory within one threadblock
    @unroll for warp_tile = parallellise(block_tile.MN, Tile(conf.mem_cd_warp), warpId, conf.warps_per_block)
        @unroll for thread_tile = parallellise(warp_tile, Tile(conf.mem_cd_thread), laneId, 32)
            x = Layout.load(conf.shared_d_layout, shmem_d, thread_tile)
            x = apply_bias(x, ep.bias_pointer, translate_base(thread_tile, (M = block_i, N = block_j)))
            x = transform(x, thread_tile)
            Layout.store!(conf.global_d_layout, d, x, translate_base(thread_tile, (M = block_i, N = block_j)))
        end
    end
end

end


text:
"""
    poisson_time(a, b, u)

Obtaining waiting time for inhomogeneous Poisson Process
with rate of the form λ(t) = (a + b*t)^+, `a`,`b` ∈ R, `u` uniform random variable
"""
function poisson_time(a, b, u)
    if b > 0
        if a < 0
            return sqrt(-log(u)*2.0/b) - a/b
        else # a[i]>0
            return sqrt((a/b)^2 - log(u)*2.0/b) - a/b
        end
    elseif b == 0
        if a > 0
            return -log(u)/a
        else # a[i] <= 0
            return Inf
        end
    else # b[i] < 0
        if a <= 0
            return Inf
        elseif -log(u) <= -a^2/b + a^2/(2*b)
            return -sqrt((a/b)^2 - log(u)*2.0/b) - a/b
        else
            return Inf
        end
    end
end

"""
    poisson_time((a, b, c), u)

Obtaining waiting time for inhomogeneous Poisson Process
with rate of the form λ(t) = c + (a + b*t)^+,
where `c`> 0 ,`a, b` ∈ R, `u` uniform random variable
"""
function poisson_time((a,b,c)::NTuple{3}, u)
    if b > 0
        if a < 0
            if -c*a/b + log(u) < 0.0
                return sqrt(-2*b*log(u) + c^2 + 2*a*c)/b - (a+c)/b
            else
                return -log(u)/c
            end
        else # a >0
            return sqrt(-log(u)*2.0*b + (a+c)^2)/b - (a+c)/b
        end
    elseif b == 0
        if a > 0
            return -log(u)/(a + c)
        else # a <= 0
            return -log(u)/c
        end
    else # b < 0
        if a <= 0.0
            return -log(u)/c
        elseif  - c*a/b - a^2/(2*b)  + log(u) > 0.0
            return +sqrt((a+c)^2 - 2.0*log(u)*b)/b - (a+c)/b
        else
            return (-log(u)+ a^2/(2*b))/c
        end
    end
end

"""
    poisson_time(a[, u])

Obtaining waiting time for homogeneous Poisson Process
with rate of the form λ(t) = a, `a` ≥ 0, `u` uniform random variable
"""
function poisson_time(a::Number, u::Number)
    -log(u)/a
end

function poisson_time(rng, a)
    randexp(rng)/a
end
function poisson_time(a)
    randexp()/a
end

poisson_time((a, b)::Tuple, u=randn()) = poisson_time(a, b, u)


text:
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]

###############################
print("Testing ACQRS_sub: ")
nx = 4
i = nx+1
A,C,Q,R,S = ACQRS_sub(y,nx,i)
if sum(A.^2) > 1e-6
	println("OK")
end


text:
@pure length{T<:StaticArray}(a::Union{T,Type{T}}) = prod(size(a))

@pure function size{T<:StaticArray}(a::Union{T,Type{T}}, d::Integer)
    s = size(a)
    return (d <= length(s) ? s[d] : 1)
end

# This seems to confuse Julia a bit in certain circumstances (specifically for trailing 1's)
function Base.isassigned(a::StaticArray, i::Int...)
    ii = sub2ind(size(a), i...)
    1 <= ii <= length(a) ? true : false
end

Base.linearindexing{T<:StaticArray}(::Union{T,Type{T}}) = Base.LinearFast()

# Default type search for similar_type
@pure similar_type{SA<:StaticArray}(::Union{SA,Type{SA}}) = SA
@pure function similar_type{SA<:StaticArray,N,T}(::Union{SA,Type{SA}}, ::Type{T}, sizes::NTuple{N,Int})
    similar_type(similar_type(SA, T), sizes)
end
@pure function similar_type{SA<:StaticArray,T}(::Union{SA,Type{SA}}, ::Type{T}, size::Int)
    similar_type(similar_type(SA, T), size)
end
@generated function similar_type{SA<:StaticArray,T}(::Union{SA,Type{SA}}, ::Type{T})
    # This function has a strange error (on tests) regarding double-inference, if it is marked @pure
    if T == eltype(SA)
        return SA
    end

    primary_type = (SA.name.primary)
    eltype_param_number = super_eltype_param(primary_type)
    if isnull(eltype_param_number)
        if ndims(SA) == 1
            return SVector{length(SA), T}
        elseif ndims(SA) == 2
            sizes = size(SA)
            return SMatrix{sizes[1], sizes[2], T, prod(sizes)}
        else
            sizes = size(SA)
            return SArray{sizes, T, length(sizes), prod(sizes)}
        end
    end

    T_out = primary_type
    for i = 1:length(T_out.parameters)
        if i == get(eltype_param_number)
            T_out = T_out{T}
        else
            T_out = T_out{SA.parameters[i]}
        end
    end
    return T_out
end

@pure function super_eltype_param(T)
    T_super = supertype(T)
    if T_super == Any
        error("Unknown error")
    end

    if T_super.name.name == :StaticArray
        if isa(T_super.parameters[1], TypeVar)
            tv = T_super.parameters[1]
            for i = 1:length(T.parameters)
                if tv == T.parameters[i]
                    return Nullable{Int}(i)
                end
            end
            error("Unknown error")
        end
        return Nullable{Int}()
    else
        param_number = super_eltype_param(T_super)
        if isnull(param_number)
            return Nullable{Int}()
        else
            tv = T_super.parameters[get(param_number)]
            for i = 1:length(T.parameters)
                if tv == T.parameters[i]
                    return Nullable{Int}(i)
                end
            end
            error("Unknown error")
        end
    end
end

# Some fallbacks

@pure similar_type{SA<:StaticArray}(::Union{SA,Type{SA}}, size::Int) = SVector{size, eltype(SA)}
@pure similar_type{SA<:StaticArray}(::Union{SA,Type{SA}}, sizes::Tuple{Int}) = SVector{sizes[1], eltype(SA)}

@pure similar_type{SA<:StaticArray}(::Union{SA,Type{SA}}, sizes::Tuple{Int,Int}) = SMatrix{sizes[1], sizes[2], eltype(SA), sizes[1]*sizes[2]}

@pure similar_type{SA<:StaticArray,N}(::Union{SA,Type{SA}}, sizes::Tuple{Vararg{Int,N}}) = SArray{sizes, eltype(SA), N, prod(sizes)}

# Some specializations for the mutable case
@pure similar_type{MA<:Union{MVector,MMatrix,MArray}}(::Union{MA,Type{MA}}, size::Int) = MVector{size, eltype(MA)}
@pure similar_type{MA<:Union{MVector,MMatrix,MArray}}(::Union{MA,Type{MA}}, sizes::Tuple{Int}) = MVector{sizes[1], eltype(MA)}

@pure similar_type{MA<:Union{MVector,MMatrix,MArray}}(::Union{MA,Type{MA}}, sizes::Tuple{Int,Int}) = MMatrix{sizes[1], sizes[2], eltype(MA), sizes[1]*sizes[2]}

@pure similar_type{MA<:Union{MVector,MMatrix,MArray},N}(::Union{MA,Type{MA}}, sizes::Tuple{Vararg{Int,N}}) = MArray{sizes, eltype(MA), N, prod(sizes)}

# And also similar() returning mutable StaticArrays
@inline similar{SV <: StaticVector}(::SV) = MVector{length(SV),eltype(SV)}()
@inline similar{SV <: StaticVector, T}(::SV, ::Type{T}) = MVector{length(SV),T}()
@inline similar{SA <: StaticArray}(::SA, sizes::Tuple{Int}) = MVector{sizes[1], eltype(SA)}()
@inline similar{SA <: StaticArray}(::SA, size::Int) = MVector{size, eltype(SA)}()
@inline similar{T}(::StaticArray, ::Type{T}, sizes::Tuple{Int}) = MVector{sizes[1],T}()
@inline similar{T}(::StaticArray, ::Type{T}, size::Int) = MVector{size,T}()

@inline similar{SM <: StaticMatrix}(m::SM) = MMatrix{size(SM,1),size(SM,2),eltype(SM),length(SM)}()
@inline similar{SM <: StaticMatrix, T}(::SM, ::Type{T}) = MMatrix{size(SM,1),size(SM,2),T,length(SM)}()
@inline similar{SA <: StaticArray}(::SA, sizes::Tuple{Int,Int}) = MMatrix{sizes[1], sizes[2], eltype(SA), sizes[1]*sizes[2]}()
@inline similar(a::StaticArray, T::Type, sizes::Tuple{Int,Int}) = MMatrix{sizes[1], sizes[2], T, sizes[1]*sizes[2]}()

@inline similar{SA <: StaticArray}(m::SA) = MArray{size(SA),eltype(SA),ndims(SA),length(SA)}()
@inline similar{SA <: StaticArray,T}(m::SA, ::Type{T}) = MArray{size(SA),T,ndims(SA),length(SA)}()
@inline similar{SA <: StaticArray,N}(m::SA, sizes::NTuple{N, Int}) = MArray{sizes,eltype(SA),N,prod(sizes)}()
@inline similar{SA <: StaticArray,N,T}(m::SA, ::Type{T}, sizes::NTuple{N, Int}) = MArray{sizes,T,N,prod(sizes)}()

# This is used in Base.LinAlg quite a lot, and it impacts type stability
# since some functions like expm() branch on a check for Hermitian or Symmetric
# TODO much more work on type stability. Base functions are using similar() with
# size, which poses difficulties!!
@generated function Base.full{T,SM<:StaticMatrix}(sym::Symmetric{T,SM})
    exprs_up = [i <= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(m[$(sub2ind(size(SM), j, i))]) for i = 1:size(SM,1), j=1:size(SM,2)]
    exprs_lo = [i >= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(m[$(sub2ind(size(SM), j, i))]) for i = 1:size(SM,1), j=1:size(SM,2)]

    return quote
        $(Expr(:meta, :inline))
        m = sym.data
        if sym.uplo == 'U'
            @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_up...))))
        else
            @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_lo...))))
        end
    end
end

@generated function Base.full{T,SM<:StaticMatrix}(sym::Hermitian{T,SM})
    exprs_up = [i <= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(conj(m[$(sub2ind(size(SM), j, i))])) for i = 1:size(SM,1), j=1:size(SM,2)]
    exprs_lo = [i >= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(conj(m[$(sub2ind(size(SM), j, i))])) for i = 1:size(SM,1), j=1:size(SM,2)]

    return quote
        $(Expr(:meta, :inline))
        m = sym.data
        if sym.uplo == 'U'
            @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_up...))))
        else
            @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_lo...))))
        end
    end
end

# Reshape uses types to specify size (and also conveniently, the output type)
@generated function reshape{SA<:StaticArray}(a::StaticArray, ::Type{SA})
    if !(SA <: StaticVector) && length(a) != length(SA)
        error("Static array of size $(size(a)) cannot be reshaped to size $(size(SA))")
    end

    return quote
        $(Expr(:meta, :inline))
        return SA(tuple(a))
    end
end

function reshape{SA<:StaticArray}(a::AbstractArray, ::Type{SA})
    if !(SA <: StaticVector) && length(a) != length(SA)
        error("Static array of size $(size(a)) cannot be reshaped to size $(size(SA))")
    end

    return SA((a...))
end


text:
using HeadTurn
using Test

@testset "HeadTurn.jl" begin
    # Write your tests here.
end


text:
using PyCall
using PyPlot
using DelimitedFiles
using StatsBase
using ArchGDAL
const AG = ArchGDAL
using Dates
using DataFrames

const ICE_SAMPLE_NUM = 100
const CNFL_SAMPLE_NUM = 100

# Getting filenames for files where we may have locally saved datasets of one solar year's worth of PV system generation data
function get_mc_data_filename(num_samples, cnfl)
    provider = length(cnfl) == 0 ? "ALL" : (cnfl[1] == true ? "CNFL" : "ICE")
    num_samples = (length(cnfl) == 1) ? (cnfl[1] == true ? min(num_samples, CNFL_SAMPLE_NUM) : min(num_samples, ICE_SAMPLE_NUM)) : num_samples
    mc_filename = string("data/monte_carlo_data/", provider, "_", string(num_samples), ".txt")
    return (num_samples, mc_filename)
end

# Getting filenames for files where we may have locally saved the coordinates corresponding to the above PV system generation datasets
function get_mc_coords_filename(num_samples, cnfl)
    provider = length(cnfl) == 0 ? "ALL" : (cnfl[1] == true ? "CNFL" : "ICE")
    num_samples = (length(cnfl) == 1) ? (cnfl[1] == true ? min(num_samples, CNFL_SAMPLE_NUM) : min(num_samples, ICE_SAMPLE_NUM)) : num_samples
    mc_filename = string("data/monte_carlo_data/", provider, "_", string(num_samples), "_COORDS.txt")
    return mc_filename
end

# Tries to get those preexisting files, and averages all of the PV output timeseries with a strict equal weighting, in order to get an “average PV system output”. If the files aren’t there, it will call the next function, to make the NSRDB+SAM calls
function averaged_monte_carlo_solar_output(;num_samples=100, cnfl=[])
    # Get the individual data points for various locations across the country
    num_samples, mc_filename = get_mc_data_filename(num_samples, cnfl)
    if isfile(mc_filename)
        mc_pv_output = readdlm(mc_filename, '\t', Float64, '\n')
    else
        println("have to generate new data")
        mc_pv_output, mc_locations = monte_carlo_solar_output(num_samples, cnfl)
        open(mc_filename, "w") do io
           writedlm(io, mc_pv_output)
        end
        coords_filename = get_mc_coords_filename(num_samples, cnfl)
        open(coords_filename, "w") do io
           writedlm(io, mc_locations)
        end
    end

    # Average all of them, to return a normalized data set for a hypothetical "solar year"
    return sum(mc_pv_output, dims=2) ./ num_samples
end

# Calling to NSRDB+SAM to obtain Monte Carlo dataset of solar years
function monte_carlo_solar_output(num_samples, cnfl; use_cached=false)
    # Has an option to try and read from a cache, for testing
    if (use_cached)
        _, mc_filename = get_mc_data_filename(num_samples, cnfl)
        mc_coords_filename = get_mc_coords_filename(num_samples, cnfl)
        if isfile(mc_filename) && isfile(mc_coords_filename)
            mc_pv_output = readdlm(mc_filename, '\t', Float64, '\n')
            mc_coords = readdlm(mc_filename, '\t', Float64, '\n')
            return (mc_pv_output, mc_coords)
        end
    end

    # Loads GIS files for population density (to weight the location sampling), the CNFL area (for CNFL vs ICE split), the protected area (to exclude), the overall map of Costa Rica (to make sure that generated points are within the nation's borders)
    pop_density_filename = "data/distritos/distritos2008crtm05.shp"
    cnfl_gis_filename = "data/area_CNFL"
    area_protegidas_gis_filename = "data/Areaprotegidas"
    cr_map_filename = "data/CRI_adm/CRI_adm0.shp"
    ice_gis_filename = "data/fEmpresaDistribuidoras_2020"
    pv_outputs_array = Array{Float64,1}()

    coords = []

    AG.registerdrivers() do
        # Get data corresponding to the location of CNFL-serviced districts, for later location sampling
        AG.read(cnfl_gis_filename) do cnfl_gis_dataset
            cnfl_gis = AG.getlayer(cnfl_gis_dataset, 0)

            # Get data corresponding to the map of ICE service territory, for later location sampling
            AG.read(ice_gis_filename) do ice_gis_dataset
                ice_gis = AG.getlayer(ice_gis_dataset, 0)

                # Get data corresponding to the location of protected areas, to exclude from all location sampling
                AG.read(area_protegidas_gis_filename) do area_protegidas_gis_dataset
                    area_protegidas_gis = AG.getlayer(area_protegidas_gis_dataset, 0)

                    # The above 2 GIS data files are in a format which does not line up with latitude and longitude.
                    # We will run a transform to project it into lat-long
                    AG.importPROJ4("proj +proj=longlat") do target

                        # Get population density information
                        # This can either be derived from cnfl_gis, or we need a nation-wide district file
                        district_centroids = []
                        district_populations = []

                        if isassigned(cnfl, 1) && cnfl[1]
                            num_features = AG.nfeature(cnfl_gis)
                            for i in 1:(num_features)
                                ArchGDAL.getfeature(cnfl_gis, i - 1) do feature
                                    male_population = AG.asint(feature, AG.getfieldindex(feature, "POB_2000_M"))
                                    female_population = AG.asint(feature, AG.getfieldindex(feature, "POB_2000_H"))
                                    total_population = male_population + female_population
                                    push!(district_populations, total_population)

                                    geom = AG.getgeomfield(feature, 0)
                                    source = AG.getspatialref(geom)
                                    AG.createcoordtrans(source, target) do transform
                                        AG.transform!(geom, transform)
                                        centroid = AG.centroid(geom)
                                        push!(district_centroids, centroid)
                                    end
                                end
                            end
                        else
                            AG.read(pop_density_filename) do districts_gis_dataset
                                districts_gis = AG.getlayer(districts_gis_dataset, 0)
                                println(districts_gis)
                                num_features = AG.nfeature(districts_gis)
                                for i in 1:(num_features)
                                    ArchGDAL.getfeature(districts_gis, i - 1) do feature
                                        male_population = AG.asint(feature, AG.getfieldindex(feature, "POB_2000_M"))
                                        female_population = AG.asint(feature, AG.getfieldindex(feature, "POB_2000_H"))
                                        total_population = male_population + female_population
                                        push!(district_populations, total_population)

                                        geom = AG.getgeomfield(feature, 0)
                                        source = AG.getspatialref(geom)
                                        AG.createcoordtrans(source, target) do transform
                                            AG.transform!(geom, transform)
                                            centroid = AG.centroid(geom)
                                            push!(district_centroids, centroid)
                                        end
                                    end
                                end
                            end
                        end

                        # Convert population density into a set of weights for an Empirical CDF
                        weights = FrequencyWeights(collect(Iterators.flatten(district_populations)))

                        # Run a weighted sampling to obtain the coordinates to sample NSRDB from

                        # Obtain sample PV output for each location
                        while length(coords) < num_samples

                            possible_coords = sample(district_centroids, weights)

                            # Determine whether these coordinates have been already added, or are part of excluded parts of CR
                            if !in(possible_coords, coords)
                                can_add = false
                                ag_coords = possible_coords

                                # Is this point even on land, or is it over the ocean?
                                AG.read(cr_map_filename) do cr_map_dataset
                                    cr_map_gis = AG.getlayer(cr_map_dataset, 0)
                                    num_features = AG.nfeature(cr_map_gis)
                                    for i in 1:(num_features)
                                        ArchGDAL.getfeature(cr_map_gis, i - 1) do feature
                                            geom = AG.getgeomfield(feature, 0)
                                            if AG.contains(geom, ag_coords)
                                                can_add = true
                                            end
                                        end
                                    end
                                end

                                if !can_add
                                    continue
                                end

                                # Add a conditional check for whether the location is within CNFL or ICE service area
                                # Only do this when [cnfl] is defined: otherwise, assume we want national data
                                # If this is CNFL, we're using the CNFL districts so it's guaranteed to be correct
                                if isassigned(cnfl, 1) && !cnfl[1]
                                    can_add = false
                                    ArchGDAL.getfeature(ice_gis, 4) do feature
                                        geom = AG.getgeomfield(feature, 0)
                                        source = AG.getspatialref(geom)
                                        AG.createcoordtrans(source, target) do transform
                                            AG.transform!(geom, transform)
                                            if AG.contains(geom, ag_coords)
                                                can_add = true
                                                println(string(possible_coords, " is within the ICE area"))
                                            end
                                        end
                                    end

                                    if !can_add
                                        continue
                                    end
                                end

                                # Check that the generated point isn't in a nationally protected area/national forest
                                num_features = AG.nfeature(area_protegidas_gis)
                                for i in 1:(num_features)
                                    ArchGDAL.getfeature(area_protegidas_gis, i - 1) do feature
                                        geom = AG.getgeomfield(feature, 0)
                                        source = AG.getspatialref(geom)
                                        AG.createcoordtrans(source, target) do transform
                                            AG.transform!(geom, transform)
                                            if AG.contains(geom, ag_coords)
                                                println(string("but ", possible_coords, " is within a protected area so we can't use it"))
                                                can_add = false # Not using this point because it is in an uninhabited location
                                            end
                                        end
                                    end
                                end

                                if !can_add
                                    continue
                                end

                                # Make the call out to NSRDB + SAM to obtain pv_output for that location
                                pv_output = -1
                                try
                                    lat, lon = AG.gety(possible_coords, 0), AG.getx(possible_coords, 0)
                                    pv_output = convert(Array{Float64,1},get_nsrdb_sam_pv_output(lat=lat, lon=lon))
                                catch e
                                    println("but the NSRDB call failed, so we won't use it")
                                    can_add = false # Not using this point because NSRDB doesn't like it
                                end

                                if can_add
                                    push!(coords, possible_coords)
                                    println(length(coords))
                                    if length(pv_outputs_array) == 0
                                        append!(pv_outputs_array, pv_output)
                                    else
                                        pv_outputs_array = hcat(pv_outputs_array, pv_output)
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    # Return the outputs of SAM, along with a DataFrame of all of the coordinates for future reference
    coords = convert(Array{Float64,2}, DataFrame(Latitude = [AG.gety(coord, 0) for coord in coords], Longitude = [AG.getx(coord, 0) for coord in coords]))
    return (pv_outputs_array, coords)

end

# A wrapper to either call get_nsrdb_sam_df, or if we want to skip all this and just return a hardcoded value
function get_nsrdb_sam_pv_output(;lat=9.817934, lon=-84.070552, tz=-6, year=2010, pipeline=true)
    if pipeline == true
        nsrdb_sam_df = get_nsrdb_sam_df(lat, lon, tz, year)
        pv_output = values(nsrdb_sam_df["Generation"])
    else
        pv_output=readdlm("data/pv_output.txt", '\t', Float64, '\n')
    end
    return pv_output
end

# Obtain the NSRDB request URL, and invoke the python function which will actually call to the NREL services
function get_nsrdb_sam_df(lat, lon, tz, year)
    request_url = get_nsrdb_request_url(lat, lon, year);
    py"""
    import sys
    sys.path.insert(0, ".")
    sys.path.insert(0, "./functions")
    """
    call_nsrdb_and_ssc = pyimport("nsrdb_python")["call_nsrdb_and_ssc"];
    nsrdb_sam_df = call_nsrdb_and_ssc(request_url, lat, lon, tz);
    return nsrdb_sam_df;
    end

# Construct the URL to use to get the NSRDB data (this requires an external file with credentials)
function get_nsrdb_request_url(lat,lon,year)
    # Declare all variables as strings. Spaces must be replaced with "+", i.e., change "John Smith" to "John+Smith".

    # Set the attributes to extract (e.g., dhi, ghi, etc.), separated by commas.
    attributes = "ghi,dhi,dni,wind_speed,air_temperature,solar_zenith_angle"

    # Set leap year to true or false. True will return leap day data if present, false will not.
    leap_year = "false"

    # Set time interval in minutes, i.e., "30" is half hour intervals. Valid intervals are 30 & 60.
    interval = "60"

    # Specify Coordinated Universal Time (UTC), "true" will use UTC, "false" will use the local time zone of the data.
    # NOTE: In order to use the NSRDB data in SAM, you must specify UTC as "false". SAM requires the data to be in the
    # local time zone.
    utc = "false"

    # Have a file named "nrel-credentials.txt", with line-separated
    # - Full+Name
    # - Email address
    # - Your+Affiliation
    # - API Key
    credential_file = open("nrel-credentials.txt")
    nrel_credentials = readlines(credential_file)

    # Your full name, use "+" instead of spaces.
    your_name = nrel_credentials[1]

    # Your email address
    your_email = nrel_credentials[2]

    # Your affiliation
    your_affiliation = nrel_credentials[3]

    # You must have an NSRDB api key
    api_key = nrel_credentials[4]
    close(credential_file)

    # Please join our mailing list so we can keep you up-to-date on new developments.
    mailing_list = "false";

    # Your reason for using the NSRDB.
    reason_for_use = "data+analysis"

    # Declare url string
    url = "http://developer.nrel.gov/api/solar/nsrdb_psm3_download.csv?wkt=POINT($(lon)%20$(lat))&names=$(year)&leap_day=$(leap_year)&interval=$(interval)&utc=$(utc)&full_name=$(your_name)&email=$(your_email)&affiliation=$(your_affiliation)&mailing_list=$(mailing_list)&reason=$(reason_for_use)&api_key=$(api_key)&attributes=$(attributes)";

    return url;
end

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

# Unused function to plot PySAM output
function plot_pysam_output(df, i=5030, j=40)
    plt.style.use("seaborn")
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax2 = ax.twinx()
    ax.plot(df[[:GHI]][i:i+j,:][:1])
    ax.plot(df[[:DNI]][i:i+j,:][:1])
    ax.plot(df[[:DHI]][i:i+j,:][:1])
    ax.plot(df[[Symbol("Solar Zenith Angle")]][i:i+j,:][:1])
    ax2.plot(df[[:Generation]][i:i+j,:][:1], color="y")
    ax.grid()
    ax.set_ylabel("W/m2")
    ax2.set_ylabel("kW")
    ax.legend([:GHI, :DNI, :DHI, Symbol("Solar Zenith Angle")], loc="upper left")
    ax2.legend([:Generation], loc="upper right")
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
import Colors
using Colors: RGBA, N0f8

const RayColor = RGBA{N0f8}

raycolor(rgba...) = RayColor(map(Base.Fix1(reinterpret, N0f8)∘Cuchar, rgba)...)

"raylib defined Light Gray"
const LIGHTGRAY = raycolor( 200, 200, 200, 255 )
"raylib defined Gray"
const GRAY      = raycolor( 130, 130, 130, 255 )
"raylib defined Dark Gray"
const DARKGRAY  = raycolor( 80, 80, 80, 255 )
"raylib defined Yellow"
const YELLOW    = raycolor( 253, 249, 0, 255 )
"raylib defined Gold"
const GOLD      = raycolor( 255, 203, 0, 255 )
"raylib defined Orange"
const ORANGE    = raycolor( 255, 161, 0, 255 )
"raylib defined Pink"
const PINK      = raycolor( 255, 109, 194, 255 )
"raylib defined Red"
const RED       = raycolor( 230, 41, 55, 255 )
"raylib defined Maroon"
const MAROON    = raycolor( 190, 33, 55, 255 )
"raylib defined Green"
const GREEN     = raycolor( 0, 228, 48, 255 )
"raylib defined Lime"
const LIME      = raycolor( 0, 158, 47, 255 )
"raylib defined Dark Green"
const DARKGREEN = raycolor( 0, 117, 44, 255 )
"raylib defined Sky Blue"
const SKYBLUE   = raycolor( 102, 191, 255, 255 )
"raylib defined Blue"
const BLUE      = raycolor( 0, 121, 241, 255 )
"raylib defined Dark Blue"
const DARKBLUE  = raycolor( 0, 82, 172, 255 )
"raylib defined Purple"
const PURPLE    = raycolor( 200, 122, 255, 255 )
"raylib defined Violet"
const VIOLET    = raycolor( 135, 60, 190, 255 )
"raylib defined Dark Purple"
const DARKPURPLE = raycolor( 112, 31, 126, 255 )
"raylib defined Beige"
const BEIGE     = raycolor( 211, 176, 131, 255 )
"raylib defined Brown"
const BROWN     = raycolor( 127, 106, 79, 255 )
"raylib defined Dark Brown"
const DARKBROWN = raycolor( 76, 63, 47, 255 )

"raylib defined White"
const WHITE     = raycolor( 255, 255, 255, 255 )
"raylib defined Black"
const BLACK     = raycolor( 0, 0, 0, 255 )
"raylib defined Transparent"
const BLANK     = raycolor( 0, 0, 0, 0 )
"raylib defined Magenta"
const MAGENTA   = raycolor( 255, 0, 255, 255 )
"raylib defined Ray White"
const RAYWHITE  = raycolor( 245, 245, 245, 255 )


text:
#
# Manages organisms and their world
# Events:
#   created     {ManagerData}                Fires after manager creates it's data types
#   afterbackup {ManagerData}                Fires after backup has created
#   killorganism{ManagerData, Organism}      Fires if organism has die
#   bornorganism{ManagerData, Creature}      Fire if organism has borned
#   mutations   {ManagerData, Creature, Int} Fires if mutations have occures
#   clone       {ManagerData, UInt, UInt}    Fires if after clonning
#   yield       {ManagerData}                Fires if yield() is called
#   yieldto     {ManagerData}                Fires if yieldto() is called
#   iteration   {ManagerData, Float64}       Is fired on every iteration (see IPS)
#   updateenergy{ManagerData}                Is fired on every organisms energy update
#   killorganism{ManagerData, Organism}      Is called if organism has been killed
#   dotrequest  {ManagerData}                Is called if step (dot) request is occured
#   stepyield   {ManagerData}                Is called if one step related yield was called
#   backup      {ManagerData}                Fired after backup has created
#   ips         {ManagerData}                Fired after ips parameter has updated
# TODO: gather all events
#
#
# @singleton
# @author DeadbraiN
#
# TODO: describe included files and how to split Manager on files.
# TODO:   - Manager may call every function from any file
# TODO:   - File may call every function except other files
# TODO:   - If file calls Manager function it hould add Manager. before calling
# TODO: describe that manager is a mediator between all other objects
# TODO: like world, terminal and so on.
# TODO: add remote functions for changing period and probs
# TODO: add command line parameter for creating default config file
# TODO: add create method. It should return ManagerData() type
# TODO: describe frozen organisms conception
#
module Manager
  #
  # Macroses
  #
  import Config.@if_debug
  import Config.@if_test
  import Config.@if_profile
  #
  # Dependencies
  #
  import Creature
  import Mutator
  import World
  import Helper
  import Event
  import Server
  import Client
  import Connection
  import CommandLine
  import RpcApi
  import FastApi
  import Config
  import ManagerTypes
  import Dots

  export create
  export run
  export destroy
  #
  # This is how we collect Manager module from it's parts(files)
  #
  include("ManagerOrganism.jl")
  include("ManagerRpc.jl")
  include("ManagerBackup.jl")
  include("ManagerParams.jl")
  include("ManagerPlugins.jl")
  #
  # Creates manager related data instance. It will be passed to all
  # manager methods
  # @param cfg Configuration object for manager
  #
  function create(cfg::Config.ConfigData = Config.create())
    local man::ManagerTypes.ManagerData
    local obs::Event.Observer = Event.create()
    local tasks::Array{ManagerTypes.OrganismTask, 1} = Array{ManagerTypes.OrganismTask, 1}(cfg.worldMaxOrgs)
    local i::Int
    local task::Task = current_task()
    local org::Creature.Organism
    #
    # We create temporary organisms to prevent keeping empty slots in array. All
    # these organisms will be removed later, during application working.
    #
    for i = 1:cfg.worldMaxOrgs
      org       = Creature.create(cfg, i)
      org.alive = false
      tasks[i]  = ManagerTypes.OrganismTask(UInt(0), task, org)
    end

    man = ManagerTypes.ManagerData(
      cfg,                                                                      # cfg
      World.create(cfg.worldWidth, cfg.worldHeight),                            # world
      Dict{Int, Creature.Organism}(),                                           # positions
      Dict{UInt, Creature.Organism}(),                                          # organisms
      tasks,                                                                    # tasks
      CommandLine.create(),                                                     # params
      UInt(1),                                                                  # organismId
      UInt(0),                                                                  # totalOrganisms
      CommandLine.has(CommandLine.create(), ARG_QUIET),                         # quiet
      function() end,                                                           # dotCallback
      function() end,                                                           # moveCallback
      current_task(),                                                           # task
      0,                                                                        # maxEnergy
      0.0,                                                                      # ips
      Int[1:cfg.worldMaxOrgs...],                                               # killed
      Dict{String, ManagerTypes.Plugin}(),                                      # plugins
      obs                                                                       # obs
    )
    man.cons = _createConnections(man)
    initPlugins(man)

    man
  end
  #
  # Runs Manager instance, one world, server an so on... Blocking
  # function.
  # @param man Manager data type
  # @param recover true if we have to recover from last backup
  # @return {Bool} run status
  #
  function run(man::ManagerTypes.ManagerData, recover::Bool = false)
    #
    # must be started from 1! because we use % (div reminder) for checking
    # if some config related to iterations occures. In case of 0 value
    # all configs will be run at the beginning.
    #
    local counter   ::Int = 1
    local stamp     ::Float64 = time()
    local istamp    ::Float64 = stamp
    local bstamp    ::Float64 = stamp
    local ystamp    ::Float64 = stamp
    local cons      ::ManagerTypes.Connections = man.cons
    local tasks     ::Array{ManagerTypes.OrganismTask, 1} = man.tasks
    local cfg       ::Config.ConfigData = man.cfg
    local backups   ::Int = 0
    local needYield ::Bool = false
    @if_profile local i::Int = 0

    try
      #
      # This server is listening for all other managers and remote
      # terminal. It runs obtained commands and send answers back.
      # In other words, it works like RPC runner... Fast server is
      # listening for "fast" clients and works in "fast" mode.
      #
      Server.run(man.cons.server)
      # TODO: possibly, we don't need to run this server due to performance issue
      Server.run(man.cons.fastServer)
      #
      # If user set up some amount of organisms they will be created
      # in this call. If we are in recover mode, then this step should
      # be skipped.
      #
      if !recover
        if cfg.worldStartEnergyDots + cfg.worldStartEnergyInDot > 0
          setRandomEnergy(man, cfg.worldStartEnergyDots, cfg.worldStartEnergyInDot)
        end
        if cfg.orgStartAmount > 0 createOrganisms(man) end
      end
      #
      # This is main infinite loop. It manages input connections
      # and organism's tasks switching.
      #
      while true
        #
        # We have to wait while all clients are ready for streaming. This
        # is because the error in serializer. See issue for details:
        # https://github.com/JuliaLang/julia/issues/16746
        #
        if cons.streamInit yield(); Event.fire(man.obs, "yield", man); continue end
        #
        # This is global time stamp in seconds
        #
        stamp = time()
        #
        # This call runs all organism related tasks one by one
        #
        counter = _updateOrganisms(man, counter, needYield)
        #
        # Checks if total amount of energy in a world is less then
        # minimum, according to the configuration.
        #
        if counter % cfg.worldEnergyCheckPeriod === 0 && cfg.worldEnergyCheckPeriod > 0
          _updateWorldEnergy(man)
        end
        #
        # We have to update IPS (Iterations Per Second) every second
        #
        istamp = _updateIps(man, stamp, istamp)
        #
        # This call switches between all non blocking asynchronous
        # functions (see @async macro). For example, it handles all
        # input connections for current server. But we don't need to
        # call yield() every time, because it eats CPU cicles. We
        # have to wait some period and call yield() to check if
        # sockets have data.
        # TODO: this code may be optimized. We already call yield()
        # TODO: in _onDot() handler. Second idea, that we have to use
        # TODO: yieldto(), because we know all network related tasks
        # TODO: created by @async() macro.
        ystamp, needYield = _updateTasks(man, stamp, ystamp, needYield)
        #
        # Here we make auto-backup of application if there is a time
        #
        bstamp, backups = _updateBackup(man, stamp, bstamp, backups)
        #
        # Because of issue: https://github.com/JuliaLang/julia/issues/19013
        # we have to rerun our application after several eval calls. 2400
        # obtained experimentally. This code should be removed if this
        # issue will be fixed.
        #
        if cfg.orgEvals > 2400 && backups > 5 return false end
        #
        # Means that one iteration in a system has occured
        #
        Event.fire(man.obs, "iteration", man, stamp)
        #
        # This code is used for profiling of jevo. returning true means,
        # that the process will be stopped and second run will not occures.
        #
        @if_profile if (i += 1) > cfg.modeProfilePeriod return true end
        #
        # This line is for special testing mode, which is called "iterational".
        # In this mode we may run Manager iteration by iteration and measure
        # different parameters like amount of organisms or their energy. It's
        # easier to write unit tests with such mode...
        #
        @if_test produce(counter)
      end
    catch e
      Helper.error("Manager.run(): $e")
      @if_debug showerror(STDOUT, e, catch_backtrace())
      return false
    end
    #
    # true means, that everything is okay, false - something went wrong
    #
    true
  end
  #
  # Destroys the manager and all realted components like servers,
  # clients, etc...
  # @param man Manager data object
  #
  function destroy(man::ManagerTypes.ManagerData)
    World.destroy(man.world)
    Server.stop(man.cons.server)
    Server.stop(man.cons.fastServer)
    Client.stop(man.cons.left)
    Client.stop(man.cons.right)
    Client.stop(man.cons.up)
    Client.stop(man.cons.down)

    man.positions    = Dict{Int, Creature.Organism}()
    man.organisms    = Dict{UInt, Creature.Organism}()
    man.tasks        = []
    man.killed       = []
    man.params       = Dict{String, String}()
    man.dotCallback  = Helper.emptyFn
    man.moveCallback = Helper.emptyFn
    man.cons.frozen  = Dict{UInt, Creature.Organism}()
  end
  #
  # Exit means go to main task and never back in current task again. but
  # we have to check if stopped task is current one or not. Only current
  # task requires stopping, because it running right now.
  # @param man Manager data type
  # @param task Task we have to stop
  #
  function stopTask(man::ManagerTypes.ManagerData, task::Task)
    local curTask::Task = current_task()
    if curTask !== man.task && curTask === task yieldto(man.task) end
  end

  #
  # Updates total world's energy. It shouldn't be less then minimum
  # got from config: worldEnergyCheckPercent. We calculate this percent
  # as a ration between whole world (100%) and all energy points together.
  # @param man Manager data type
  # TODO: this method is very slow!!! should be optimized!
  function _updateWorldEnergy(man::ManagerTypes.ManagerData)
    local plane::Array{UInt16, 2} = man.world.data
    local total::Int = man.world.width * man.world.height
    local energy::Int = 0
    local positions::Dict{Int, Creature.Organism} = man.positions
    local pos::Helper.Point = Helper.Point(0, 0)
    local x::Int
    local y::Int

    for x in 1:size(plane)[2]
      for y in 1:size(plane)[1]
        pos.x = x
        pos.y = y
        if !Manager._isFree(man, pos) energy += 1 end
      end
    end
    #
    # Total amount of energy is less then in config
    #
    if div(Float64(energy * 100), total) <= man.cfg.worldEnergyCheckPercent
      setRandomEnergy(man, man.cfg.worldStartEnergyDots, man.cfg.worldStartEnergyInDot)
    end
  end
  #
  # Generates unique id by world position. This macro is
  # private insode Manager module
  # @param man Manager data type
  # @param {Helper.Point} pos Unique World position
  # @return {Int} Unique pos id
  #
  function _getPosId(man::ManagerTypes.ManagerData, pos::Helper.Point) pos.y * man.world.width + pos.x end
  #
  # Checks if specified position in a world is free. Other organism
  # or an energy block may be there at the moment.
  # @param pos Position we need to check
  # @return {Bool} true - free point, false - filled point
  # TODO: this method is very slow!!!
  function _isFree(man::ManagerTypes.ManagerData, pos::Helper.Point)
    !haskey(man.positions, _getPosId(man, pos)) && World.getEnergy(man.world, pos) === Dots.INDEX_EMPTY
  end
  #
  # Updates IPS (Iterations Per second) counter and stores it in config
  # @param man Manager data type
  # @param stamp Current UNIX time stamp value
  # @param istamp IPS last UNIX time stamp value
  # @return {Float64, Int} new ips and current UNIX time stamp
  #
  function _updateIps(man::ManagerTypes.ManagerData, stamp::Float64, istamp::Float64)
    local ts::Float64 = stamp - istamp
    if ts < man.cfg.worldIpsPeriod return istamp end
    local sock::Base.TCPSocket

    man.ips = man.cfg.codeRuns / ManagerTypes.orgAmount(man) / ts
    Event.fire(man.obs, "ips", man, stamp, man.cfg.codeRuns)
    man.cfg.codeRuns = 0
    @inbounds for sock in man.cons.fastServer.socks
      if Helper.isopen(sock)
        Server.request(sock, UInt8(FastApi.API_FLOAT64), man.ips)
        Event.fire(man.obs, "request", man)
      end
    end

    stamp
  end
  #
  # Checks if it's a time to make application backup. It also checks if
  # system reaches specified amount of backups for app reset.
  # @param man Manager data type
  # @param stamp Current UNIX timestamp
  # @param bstamp Backup last UNIX time stamp value
  # @param backups Amount of backups, created from the beginning
  # @return {(Float64, Int)} Updated time stamp and amount of backups
  #
  function _updateBackup(man::ManagerTypes.ManagerData, stamp::Float64, bstamp::Float64, backups::Int)
    if stamp - bstamp >= man.cfg.backupPeriod
      if ManagerTypes.orgAmount(man) > 0
        backups += 1
        backup(man)
      end
      return stamp, backups
    end

    bstamp, backups
  end
  # TODO: describe yield() call logic
  # Checks if active servers have bytes to read. It means, that we have to call
  # yield() for this reading. yield() function will be called in organisms
  # loop in ManagerOrganism.jl file.
  # @param man Manager data type
  # @param stamp Current UNIX time stamp
  # @param ystamp yield last UNIX time stamp
  # @param needYield Flag if we need for yield() call
  # @return {(Float64, Bool)}
  #
  function _updateTasks(man::ManagerTypes.ManagerData, stamp::Float64, ystamp::Float64, needYield::Bool)
    if stamp - ystamp >= man.cfg.conYieldPeriod
      yield()
      Event.fire(man.obs, "yield", man)
      # TODO: potential problem here. this list of sockets may be expanded
      # TODO: for example in many managers mode
      return stamp, (length(man.cons.server.socks) > 0 || man.cons.streamInit)
    end

    ystamp, needYield
  end
  #
  # Updates organism's task by creating new task with passed parameters
  # @param man Instamce of ManagerData type
  # @param task OrganismTask instance, where we have to update the task
  #
  function _updateOrgTask(man::ManagerTypes.ManagerData, task::ManagerTypes.OrganismTask)
    task.task = Task(() -> Creature.born(task.organism, man.cfg, man.task))
  end
  # # TODO: do i need this?
  # # Checks id data in sockets available for reading
  # # @param man Manager data type
  # # @return {Bool}
  # #
  # function _dataAvailable(man::ManagerTypes.ManagerData)
  #   local cons::ManagerTypes.Connections = man.cons
  #   local sock::Base.TCPSocket
  #
  #   @inbounds for sock in cons.server.socks
  #     if nb_available(sock) > 0 return true end
  #   end
  #   if cons.streamInit
  #     @inbounds for sock in cons.fastServer.socks
  #       if nb_available(sock) > 0 return true end
  #     end
  #   end
  #
  #   false
  # end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
function do_main_load_loop!(cur_package::Module, cur_cabinet::FileCabinet; is_sorted::Bool=false, verbose::Bool=false)

  verbose && println("\nmain load loop:\n")

  skip_undef = true

  loaded_files_count = 1

  while !iszero(loaded_files_count)

    verbose && println("-----")

    loaded_files_count = 0

    is_sorted || sort!(
      cur_cabinet.file_infos,
      by = ( cur_info -> cur_info.time )
    )

    delete_indices = Array{Integer}(0)

    for (cur_index, cur_info) in enumerate(cur_cabinet.file_infos)

      has_undef_var = (
        ( cur_info.undef != nothing ) &&
        !isdefined(cur_package, cur_info.undef)
      )

      has_undef_var && continue

      verbose && print(cur_info.name * " ")

      did_load = attempt_file_load!(cur_package, cur_info)

      did_load && ( loaded_files_count += 1 )

      verbose && println( did_load ? "✓" : "X" )

      isempty(cur_info.unloaded_shards) &&
        push!(delete_indices, cur_index)

      ( loaded_files_count >= 3 ) && break

    end

    append!(
      cur_cabinet.load_order,
      map(cur_info -> cur_info.name, cur_cabinet.file_infos[delete_indices])
    )

    deleteat!(cur_cabinet.file_infos, delete_indices)

    iszero(length(cur_cabinet.file_infos)) && break

    if skip_undef
      iszero(loaded_files_count) || continue

      loaded_files_count = 1

      foreach(
        cur_info -> cur_info.undef = nothing,
        cur_cabinet.file_infos
      )
    else
      iszero(loaded_files_count) && continue
    end

    skip_undef = !skip_undef

  end

end


text:
using HKQM
using Rotations
using Interpolations
using WGLMakie # Or some other Makie backend

## Helper functions to initialize the system

# Generates initial state that can be optimized in decent time
function initial_state(V::ScalarOperator)
    q = QuantumState(V.elementgrid, exp.(-V.vals))
    normalize!(q)
    return q
end

function benzene_like_potential(ceg; α=0.5u"bohr^-2", depth=1u"hartree")
    # Generate benzene carbon locations
    xyz = [139., 0, 0] .* u"pm"
    rot = RotZ(π/3) # Rotation matrix
    benz = [xyz]
    for i in 1:5
        push!(benz, rot*benz[end])
    end

    # Generate potential
    r = position_operator(ceg)
    V = sum(benz) do x
        q = r + x
         return -depth * exp(-α*(q⋅q))
    end
    return V
end

## Plotting tools

# Change from elment grid (6D) to regural grid (3D)
function get3d(ψ)
    s = size(ψ)
    a = s[1]*s[end]
    q = Array{Float64}(undef,a,a,a)
    for I in 1:s[4], J in 1:s[5], K in 1:s[6]
        q[(1:s[1]).+(I-1)*s[1], (1:s[2]).+(J-1)*s[2], (1:s[3]).+(K-1)*s[3]] =  ψ[:,:,:,I,J,K]
    end
    return q
end

function plot_potential(op::ScalarOperator; title="Potential Energy [Hartree]", levels=15)
    X = vec(xgrid(op.elementgrid)) |> collect
    X =  X.* u"bohr"
    x = uconvert.(u"Å", X)
    z = length(x)/2 |> Int
    ρ = get3d(op.vals)
    f = Figure()
    Axis(f[1,1], xlabel="x [Å]", ylabel="y [Å]", title=title )
    co = contourf!(ustrip(x), ustrip(x), ρ[:,:,z]; levels=levels, colormap=:plasma)
    Colorbar(f[1, 2], co)
    return f
end

function plot_density(ψ::QuantumState; title="Electron Density", levels=15)
    function get3d_density(ψ::QuantumState)
        ρ = conj!(ψ).*ψ
        return get3d(real.(ρ))
    end
    X = vec(xgrid(ψ.elementgrid)) |> collect
    X =  X.* u"bohr"
    x = uconvert.(u"Å", X)
    z = length(x)/2 |> Int
    ρ = get3d_density(ψ)
    f = Figure()
    Axis(f[1,1], xlabel="x [Å]", ylabel="y [Å]", title=title )
    co = contourf!(ustrip(x), ustrip(x), ρ[:,:,z]; levels=levels, colormap=:plasma)
    Colorbar(f[1, 2], co)
    return f
end

function plot_current(ψ::QuantumState, j; z=0, n=32, mode_3d=false, ng=30, arrow_size=0.06, title="title")
    X = vec(xgrid(ψ.elementgrid)) |> collect
    x = ustrip.( uconvert.(u"Å", X.*u"bohr") )
    u = interpolate((x,x,x), get3d(j[1]), Gridded(Linear()))
    v = interpolate((x,x,x), get3d(j[2]), Gridded(Linear()))
    w = interpolate((x,x,x), get3d(j[3]), Gridded(Linear()))

    f(rx,ry,rz) = Point3(
        u(rx,ry,rz),
        v(rx,ry,rz),
        w(rx,ry,rz)
    )
    f(rx,ry) = Point2(
        u(rx,ry,z),
        v(rx,ry,z)
    )
    xmin = minimum(x)
    xmax = maximum(x)
    r = range(xmin; stop=xmax, length=n)
    if mode_3d
        ngrid = min(ng, 15)
        fig = streamplot(f, r, r, r; colormap = :viridis, arrow_size = arrow_size, gridsize = (ngrid, ngrid))
        return fig
    else
        fig = Figure()
        Axis(fig[1,1], xlabel="x [Å]", ylabel="y [Å]", title=title )
        sp = streamplot!(f, r, r; colormap = :magma, arrow_size = arrow_size, gridsize = (ng, ng), title=title)
        return fig
    end
end

## No magnetic field calculation

# 2^3 elements with 64^3 Gauss points in element
ceg = CubicElementGrid(6u"Å", 2, 64)

# Well depth is 1 hartree. Well width is α, values 0.5-1.0 are ok.
V = benzene_like_potential(ceg; α=0.7u"bohr^-2", depth=3u"hartree")

H = HamiltonOperator(V)

ψ = initial_state(V)

# One orbital -> two electrons
sd = SlaterDeterminant(ψ)

@info "Solving Restricted Hartree-Fock for two electrons."
@info "Doing a single Helmholtz equation. Progress bar shows how long it takes per iteration."
helmholtz_equation!(sd, H; showprogress=true)

@info "Doing 12 more Helmholtz iterations to get good convergence."
@info "Energy is printed, so that progress can be followed."
for i in 1:12
    helmholtz_equation!(sd, H)
end

 J = coulomb_operator(sd)
@info "Final total energy is $(bracket(H+0.5J,sd))"

## Magnetic field calculation

@info "Adding magnetic field to system."

B = [0., 0., 100.0].*u"T"
A = vector_potential(ceg, B...)
@info "Magnetic field is set to Bx = $(B[1]), By = $(B[2]) and Bz = $(B[3])"

Hm = HamiltonOperatorMagneticField(V,A)

@info "Doing a single Helmholtz equation in magnetic field. This will take longer due to complex numbers."
sdm = helmholtz_equation(sd, Hm; showprogress=true)

@info "Doing 5 more iterations. See energy for convergence"
for i in 1:5
    global sdm = helmholtz_equation(sdm, Hm)
end

Jm = coulomb_operator(sdm)
@info "Final energy in magnetic field is $(real(bracket(Hm+0.5Jm, sdm)))"

ψm = sdm.orbitals[1]

# Magnetic current
j = magnetic_current(ψm, Hm);

# Para magnetic current
j_para =  magnetic_current(ψm, H);

# Dia magnetic current
j_dia = map((x,y)-> x.-y, j,j_para);

## Plotting results

fpot = plot_potential(H.V)

fpsi = plot_density(ψ)

# 2d plots of magnetic currents
f = plot_current(ψm, j; title="Magnetic Current in Plane")
fp = plot_current(ψm, j_para; title="Para Magnetic Current in Plane")
fd = plot_current(ψm, j_dia; title="Dia Magnetic Current in Plane")

# 3d plot of paramagnetic current
fp = plot_current(ψm, j; mode_3d=true, ng=11)


text:
const DiffOpList = Tuple{Vararg{AbstractDifferentialOp}}

"""
    RecombineMatrix{T} <: AbstractMatrix{T}

Matrix for transformation from coefficients of the recombined basis, to the
corresponding B-spline basis coefficients.

# Extended help

The transformation matrix ``M`` is defined by
```math
ϕ_j = M_{ij} b_i,
```
where ``b_j(x)`` and ``ϕ_i(x)`` are elements of the B-spline and recombined
bases, respectively.

This matrix allows to pass from known coefficients ``u_j`` in the recombined
basis ``ϕ_j``, to the respective coefficients ``v_i`` in the B-spline basis
``b_i``:
```math
\\bm{v} = \\mathbf{M} \\bm{u}.
```

Note that the matrix is not square: it has dimensions ``N × M``, where ``N``
is the length of the B-spline basis, and ``M = N - δ`` is that of the
recombined basis (see [`RecombinedBSplineBasis`](@ref) for details).

Due to the local support of B-splines, basis recombination can be performed
by combining just a small set of B-splines near the boundaries (as discussed in
[`RecombinedBSplineBasis`](@ref)).
This leads to a recombination matrix which is almost a diagonal of ones, plus
a few extra super- and sub-diagonal elements in the upper left and lower right
corners, respectively.
The matrix is stored in a memory-efficient way that also allows fast access to
its elements.

Efficient implementations of matrix-vector products (using the `*` operator or
`LinearAlgebra.mul!`) and of left division of vectors (using `\\` or
`LinearAlgebra.ldiv!`) are included.
These two operations can be used to transform between coefficients in the
original and recombined bases.

Note that, since the recombined basis forms a subspace of the original basis
(which is related to the rectangular shape of the matrix),
it is generally not possible to obtain recombined coefficients from original
coefficients, unless the latter already satisfy the constraints encoded in the
recombined basis.
The left division operation will throw a [`NoUniqueSolutionError`](@ref) if that
is not the case.

---

    RecombineMatrix(ops::Tuple{Vararg{AbstractDifferentialOp}},
                    B::BSplineBasis, [T])

Construct recombination matrix describing a B-spline basis recombination.

The default element type `T` is generally `Float64`, except for specific
differential operators which yield a matrix of zeroes and ones, for which `Bool`
is the default.

See the [`RecombinedBSplineBasis`](@ref) constructor for details on the
`ops` argument.
"""
struct RecombineMatrix{
        T, DiffOps <: DiffOpList,
        n, n1,   # nc = n1 - n
        Corner <: SMatrix{n1,n,T},
    } <: AbstractMatrix{T}
    ops :: DiffOps  # list of differential operators for BCs
    M :: Int      # length of recombined basis (= N - 2nc)
    N :: Int      # length of B-spline basis
    ul :: Corner  # upper-left corner of matrix, size (n1, n)
    lr :: Corner  # lower-right corner of matrix, size (n1, n)
    allowed_nonzeros_per_column :: NTuple{2, Int}

    function RecombineMatrix(
            ops::DiffOpList, N::Integer,
            ul::SMatrix{n1,n}, lr::SMatrix{n1,n};
            dropped_bsplines::NTuple{2, Integer},
        ) where {n1,n}
        if n1 < n
            throw(ArgumentError("matrices must have dimensions (m, n) with m ≥ n"))
        end

        nc = length(ops)
        @assert nc == n1 - n
        M = N - 2nc
        T = eltype(ul)

        # Locality condition: make sure corner arrays are as close to banded as possible.
        allowed_nonzeros_per_column = (nc + 1) .- dropped_bsplines

        _check_locality(ul, :upper_left, allowed_nonzeros_per_column[1])
        _check_locality(lr, :lower_right, allowed_nonzeros_per_column[2])

        Corner = typeof(ul)
        Ops = typeof(ops)

        new{T, Ops, n, n1, Corner}(ops, M, N, ul, lr, allowed_nonzeros_per_column)
    end
end

function _check_locality(A, corner, allowed_nonzeros_per_column)
    m, n = size(A)
    nc = m - n  # number of constraints
    is = axes(A, 1)
    js = axes(A, 2)

    @inbounds if corner == :upper_left
        dropped_bsplines = nc + 1 - allowed_nonzeros_per_column
        for j ∈ js
            istart = j + dropped_bsplines
            iend = istart + (allowed_nonzeros_per_column - 1)
            iallowed = istart:iend
            for i ∈ is
                @assert i ∈ iallowed || iszero(A[i, j])
            end
        end
    elseif corner == :lower_right
        # allowed_nonzeros_per_column = nc - dropped_bsplines + 1
        for j ∈ js
            iallowed = j:(j + allowed_nonzeros_per_column - 1)
            for i ∈ is
                @assert i ∈ iallowed || iszero(A[i, j])
            end
        end
    end

    nothing
end

# Default element type of recombination matrix.
# In some specific cases we can use Bool...
_default_eltype(::BoundaryCondition) = Float64
_default_eltype(::DiffOpList) = Float64
_default_eltype(::Derivative{0}) = Bool  # Dirichlet BCs
_default_eltype(::Derivative{1}) = Bool  # Neumann BCs

# Case (D(0), D(1), D(2), ...)
_default_eltype(::Derivative{0}, ::Derivative{1}, ::Vararg{Derivative}) = Bool

RecombineMatrix(ops, B::BSplineBasis) =
    RecombineMatrix(ops, B, _default_eltype(ops))

RecombineMatrix(op::AbstractDifferentialOp, B::BSplineBasis, args...) =
    RecombineMatrix((op,), B, args...)

RecombineMatrix(r::DerivativeUnitRange, B::BSplineBasis, args...) =
    RecombineMatrix(Tuple(r), B, args...)

# Specialisation for Dirichlet BCs.
function RecombineMatrix(op::Tuple{Derivative{0}}, B::BSplineBasis,
                         ::Type{T}) where {T}
    _check_bspline_order(op, B)
    N = length(B)
    ul = SMatrix{1,0,T}()
    lr = copy(ul)
    RecombineMatrix(op, N, ul, lr; dropped_bsplines = (1, 1))
end

# Specialisation for Neumann BCs.
function RecombineMatrix(op::Tuple{Derivative{1}}, B::BSplineBasis,
                         ::Type{T}) where {T}
    _check_bspline_order(op, B)
    N = length(B)
    ul = SMatrix{2,1,T}(1, 1)
    lr = copy(ul)
    RecombineMatrix(op, N, ul, lr; dropped_bsplines = (0, 0))
end

# Generalisation to higher orders and to more general differential operators
# (this includes Robin BCs, for instance).
#
# If more than one operator is passed, they must be of the form
#
#    (D(0), D(1), ..., D(m - 2), D(m - 1), [D(n)]),
#
# where the last D(n) is optional, and satisfies n ≥ m + 1.
# In this case, the first `m` B-splines are dropped, and the next `n - m + 1`
# B-splines are recombined into `q = n - m` functions.
#
# That last operator is allowed to be a linear combination of `Derivative`s.
# In that case, `n` is the maximum degree of the operator.
function RecombineMatrix(ops::DiffOpList, B::BSplineBasis, ::Type{T}) where {T}
    _check_bspline_order(ops, B)
    h = order(B) ÷ 2
    # Identify operators corresponding to natural splines
    if ops === _natural_ops(Val(h))
        return RecombineMatrix(Natural(), B, T)
    end
    op = last(ops)
    n = max_order(op)
    m = _bsplines_to_drop(ops...)
    _make_matrix(Val(n + 1), Val(m), ops, B, T)
end

# Generalised natural boundary conditions.
#
# This is equivalent to:
#
#       ops = (Derivative(2), Derivative(3), …, Derivative(k ÷ 2))
#
function RecombineMatrix(::Natural, B::BSplineBasis, ::Type{T}) where {T}
    k = order(B)
    isodd(k) && throw(ArgumentError(
        "`Natural` boundary condition only supported for even-order splines (got k = $k)"
    ))

    h = k ÷ 2
    xleft, xright = boundaries(B)
    N = length(B)

    # rhs = [h, 0, 0, …, 0] (the `h` at the beginning is kind of arbitrary)
    rhs = SVector(ntuple(i -> i == 1 ? T(h) : zero(T), Val(h + 1))...)

    # Left boundary
    Tl = let x = xleft
        js = 1:(h + 1)  # indices of left B-splines

        # Evaluate B-spline derivatives at boundary.
        bs = _natural_eval_derivatives(B, x, js, Val(h), T)

        # Construct linear systems for determining recombination matrix
        # coefficients.
        A = _natural_system_matrix(bs, 1)
        T1 = A \ rhs  # coefficients associated to recombined function ϕ₁

        A = _natural_system_matrix(bs, 2)
        T2 = A \ rhs  # coefficients associated to recombined function ϕ₂

        hcat(_remove_near_zeros(T1), _remove_near_zeros(T2))
    end

    Tr = let x = xright
        js = (N - h):N
        bs = _natural_eval_derivatives(B, x, js, Val(h), T)

        A = _natural_system_matrix(bs, 1)
        T1 = A \ rhs

        A = _natural_system_matrix(bs, 2)
        T2 = A \ rhs

        hcat(_remove_near_zeros(T1), _remove_near_zeros(T2))
    end

    ops = _natural_ops(Val(h))

    RecombineMatrix(ops, N, Tl, Tr; dropped_bsplines = (0, 0))
end

function _remove_near_zeros(A::SArray; rtol = 100 * eps(eltype(A)))
    v = maximum(abs, A)
    ϵ = rtol * v
    typeof(A)((abs(x) < ϵ ? zero(x) : x) for x ∈ A)
end

@inline function _natural_ops(::Val{h}) where {h}
    @assert h ≥ 2
    (_natural_ops(Val(h - 1))..., Derivative(h))
end
@inline _natural_ops(::Val{1}) = ()

# Case h = 1: return empty (0 × 2) matrix.
_natural_eval_derivatives(B, x, js, ::Val{1}, ::Type{T}) where {T} =
    zero(SMatrix{0, 2, T})

# Evaluate derivatives 2:h of B-splines 1:(h + 1) at the boundaries.
@generated function _natural_eval_derivatives(
        B, x, js, ::Val{h}, ::Type{T},
    ) where {h, T}
    @assert h ≥ 2
    ex = quote
        @assert length(js) == $h + 1
        M = zero(MMatrix{$h - 1, $h + 1, $T})
    end
    for i ∈ 1:(h - 1)
        jlast = Symbol(:jlast_, i)
        bs = Symbol(:bs_, i)
        ileft = i == 1 ? :(nothing) : Symbol(:jlast_, i - 1)
        ex = quote
            $ex
            $jlast, $bs = evaluate_all(B, x, Derivative($i + 1), $T; ileft = $ileft)
            @assert length($bs) == order(B)
            for n ∈ axes(M, 2)
                j = js[n]
                δj = $jlast + 1 - j
                if δj ∈ eachindex($bs)
                    @inbounds M[$i, n] = $bs[δj]
                end
            end
        end
    end
    quote
        $ex
        SMatrix(M)
    end
end

# On the left boundary, `i` is the index of the resulting recombined basis
# function ϕᵢ.
function _natural_system_matrix(bs::SMatrix{hm, hp}, i) where {hm, hp}
    h = hm + 1
    @assert hp == h + 1
    M = similar(bs, Size(hp, hp))
    fill!(M, 0)
    M[1, :] .= 1  # arbitrary condition
    M[2:h, :] .= bs
    @assert i ∈ (1, 2)
    # This is a locality condition: we want the matrix to be kind of banded.
    if i == 1
        M[hp, hp] = 1
    elseif i == 2
        M[hp, 1] = 1
    end
    SMatrix(M)
end

# Case of mixed BCs (D(0), D(1), ..., D(n - 1)).
# B-splines are dropped, and new functions are not created.
# This is a generalisation of the Dirichlet case.
function _make_matrix(::Val{n}, ::Val{n}, ops, B, ::Type{T}) where {n,T}
    @assert ops === Tuple(Derivative(0:(n - 1)))
    Nc = length(ops)
    @assert Nc >= 2  # case Nc = 1 is treated by different functions
    N = length(B)
    ul = SMatrix{Nc,0,T}()
    lr = copy(ul)
    RecombineMatrix(ops, N, ul, lr; dropped_bsplines = (n, n))
end

_droplast(a) = ()
_droplast(a, etc...) = (a, _droplast(etc...)...)

# Case of single BC, or mixed BCs where the last one is "different" from the
# others.
function _make_matrix(::Val{n1}, ::Val{m}, ops, B, ::Type{T}) where {n1,m,T}
    No = length(ops)

    # Check that operators look like:
    #
    #   (D{0}, D{1}, …, D{No - 2}, some other operator of order ≥ No - 1)
    #
    @assert _droplast(ops...) === Tuple(Derivative(0:(No - 2)))
    @assert max_order(last(ops)) ≥ No - 1

    n = n1 - 1
    q = n - m
    @assert q >= 1

    Nc = length(ops)  # this is the number of constraints
    a, b = boundaries(B)
    Ca = zero(MMatrix{q + Nc, q, T})
    Cb = copy(Ca)

    N = length(B)
    op = last(ops)

    let x = a
        # Indices of B-splines to recombine.
        is = ntuple(d -> m + d, Val(q + 1))

        # Coefficients of odd derivatives must change sign, since d/dn = -d/dx
        # on the left boundary.
        opn = dot(op, LeftNormal())

        # Evaluate n-th derivatives of bⱼ at the boundary.
        # TODO replace with analytical formula?
        bs = evaluate.(B, is, x, opn)
        for l = 1:q
            b0, b1 = bs[l], bs[l + 1]
            @assert !(b0 ≈ b1)
            r = 2 / (b0 - b1)  # normalisation factor
            Ca[l + Nc - 1, l] = -b1 * r
            Ca[l + Nc, l] = b0 * r
        end
    end

    let x = b
        M = N - m  # index of last undropped B-spline
        is = ntuple(d -> M - d + 1, Val(q + 1))
        opn = dot(op, RightNormal())
        @assert opn === op  # the right normal is in the x direction
        bs = evaluate.(B, is, x, opn)
        for l = 1:q
            b0, b1 = bs[q - l + 2], bs[q - l + 1]
            @assert !(b0 ≈ b1)
            r = 2 / (b0 - b1)
            Cb[l, l] = -b1 * r
            Cb[l + 1, l] = b0 * r
        end
    end

    ul = SMatrix(Ca)
    lr = SMatrix(Cb)

    RecombineMatrix(ops, N, ul, lr; dropped_bsplines = (m, m))
end

# Verify that the B-spline order is compatible with the given differential
# operators.
function _check_bspline_order(ops::Tuple, B::BSplineBasis)
    n = max_order(ops...)
    k = order(B)
    if n >= k
        throw(ArgumentError(
            "cannot resolve operators $ops with B-splines of order $k"))
    end
    nothing
end

# Single BC: no B-splines are dropped.
# (Except if op = Derivative{0}, but that case is treated separately.)
_bsplines_to_drop(op::AbstractDifferentialOp) = 0

_unval(::Val{n}) where {n} = n
_bsplines_to_drop(ops...) = _bsplines_to_drop(Val(0), ops...) |> _unval

function _bsplines_to_drop(::Val{n}, a, b, etc...) where {n}
    if a !== Derivative(n)
        throw(ArgumentError("unsupported combination of boundary conditions"))
    end
    _bsplines_to_drop(Val(n + 1), b, etc...)
end

# This is the case of the last operator, which is allowed to be something
# different than Derivative{n}.
function _bsplines_to_drop(::Val{n}, op) where {n}
    if op === Derivative(n)
        return Val(n + 1)
    end
    if max_order(op) ≤ n
        throw(ArgumentError("order of last operator is too low"))
    end
    Val(n)
end

Base.size(A::RecombineMatrix) = (A.N, A.M)

constraints(A::RecombineMatrix) = (A.ops, A.ops)
constraints(::UniformScaling) = ((), ())

num_constraints(A::RecombineMatrix) = map(length, constraints(A))
num_constraints(::UniformScaling) = (0, 0)  # case of non-recombined bases

# Returns (left, right) tuple with maximum BC order on each boundary.
_max_order(A::RecombineMatrix) = map(ops -> max_order(ops...), constraints(A))

# Returns number of basis functions that are recombined from the original basis
# near each boundary.
# For instance, Neumann BCs have a single recombined function, ϕ_1 = b_1 + b_2;
# while mixed Dirichlet + Neumann have none, ϕ_1 = b_3.
num_recombined(A::RecombineMatrix) =
    map((m, c) -> m + 1 - c, _max_order(A), num_constraints(A))
num_recombined(::UniformScaling) = (0, 0)

# j: index in recombined basis
# M: length of recombined basis
@inline function which_recombine_block(A::RecombineMatrix, j)
    @boundscheck checkbounds(A, :, j)
    M = A.M
    nl, nr = num_recombined(A)  # left/right number of recombined bases
    j <= nl && return 1
    j > M - nr && return 3
    2
end

"""
    nzrows(A::RecombineMatrix, col::Integer) -> UnitRange{Int}

Returns the range of row indices `i` such that `A[i, col]` is non-zero.
"""
@propagate_inbounds function nzrows(
        A::RecombineMatrix, j::Integer;
        block = nothing,  # to avoid recomputing it if already known
    )
    blk = something(block, which_recombine_block(A, j))
    j += num_constraints(A)[1]
    if blk == 1
        # We take advantage of the locality condition imposed when constructing
        # the recombination matrix.
        δ = A.allowed_nonzeros_per_column[1]
        (j - δ + 1):j
    elseif blk == 2
        j:j  # shifted diagonal of ones
    else
        δ = A.allowed_nonzeros_per_column[2]
        j:(j + δ - 1)
    end :: UnitRange{Int}
end

# Pretty-printing, adapted from BandedMatrices.jl code.
function Base.replace_in_print_matrix(
        A::RecombineMatrix, i::Integer, j::Integer, s::AbstractString)
    iszero(A[i, j]) ? Base.replace_with_centered_mark(s) : s
end

@inline function Base.getindex(
        A::RecombineMatrix, i::Integer, j::Integer;
        block = nothing,
    )
    @boundscheck checkbounds(A, i, j)
    T = eltype(A)
    @inbounds blk = something(block, which_recombine_block(A, j))

    cl, cr = num_constraints(A)  # left/right constraints

    if blk == 2
        return T(i == j + cl)  # δ_{i, j+c}
    end

    if blk == 1
        C = A.ul
        if i <= size(C, 1)
            return @inbounds C[i, j] :: T
        end
        return zero(T)
    end

    # The lower-right corner starts at column h + 1.
    M = size(A, 2)
    nl, nr = num_recombined(A)
    h = M - nr

    @assert j > h
    C = A.lr
    ii = i - h - cr
    if ii ∈ axes(C, 1)
        return @inbounds C[ii, j - h] :: T
    end

    zero(T)
end

# Efficient implementation of matrix-vector multiplication.
# Generally much faster than using a regular sparse array.
function LinearAlgebra.mul!(y::AbstractVector, A::RecombineMatrix,
                            x::AbstractVector)
    checkdims_mul(y, A, x)
    N, M = size(A)

    nl, nr = num_recombined(A)
    cl, cr = num_constraints(A)
    n1 = nl + cl
    h = M - nr

    @inbounds y[1:n1] = A.ul * @view x[1:nl]

    for i = (nl + 1):h
        @inbounds y[i + cl] = x[i]
    end

    @inbounds y[(N - nr - cr + 1):N] = A.lr * @view x[(h + 1):M]

    y
end

# Five-argument mul!.
# Note that, since I have this, I could remove the 3-argument mul!  (which is
# equal to the 5-argument one with α = 1 and β = 0), but it would be a bit more
# inefficient to compute.
function LinearAlgebra.mul!(y::AbstractVector, A::RecombineMatrix,
                            x::AbstractVector, α::Number, β::Number)
    checkdims_mul(y, A, x)
    N, M = size(A)

    nl, nr = num_recombined(A)
    cl, cr = num_constraints(A)
    n1 = nl + cl
    h = M - nr

    @inbounds y[1:n1] = α * A.ul * view(x, 1:nl) + β * SVector{n1}(view(y, 1:n1))

    for i = (nl + 1):h
        @inbounds y[i + cl] = α * x[i] + β * y[i + cl]
    end

    js = (N - nr - cr + 1):N
    @inbounds y[js] =
        α * A.lr * view(x, (h + 1):M) + β * SVector{n1}(view(y, js))

    y
end

function LinearAlgebra.ldiv!(x::AbstractVector, A::RecombineMatrix,
                             y::AbstractVector)
    checkdims_mul(y, A, x)
    N, M = size(A)

    nl, nr = num_recombined(A)
    cl, cr = num_constraints(A)
    n1 = nl + cl
    h = M - nr

    @inbounds x[1:nl] = _ldiv_unique_solution(A.ul, view(y, 1:n1))

    for i = (nl + 1):h
        @inbounds x[i] = y[i + cl]
    end

    js = (N - nr - cr + 1):N
    @inbounds x[(h + 1):M] = _ldiv_unique_solution(A.lr, view(y, js))

    x
end

function LinearAlgebra.:\(A::RecombineMatrix, y::AbstractVector)
    x = similar(y, size(A, 2))
    ldiv!(x, A, y)
end

"""
    NoUniqueSolutionError <: Exception

Exception thrown when solving linear system using [`RecombineMatrix`](@ref),
when the system has no unique solution.
"""
struct NoUniqueSolutionError <: Exception end

function Base.showerror(io::IO, ::NoUniqueSolutionError)
    print(io,
          "overdetermined system has no unique solution.",
          " This means that the given function expanded in the parent basis",
          " has no exact representation (i.e. cannot be expanded) in the recombined basis."
         )
end

# Solve A \ y for overdetermined system with rectangular matrix A, assuming that
# the system has exactly one solution. Throws error if that's not the case.
function _ldiv_unique_solution(A::SMatrix, y::AbstractVector)
    N, M = size(A)
    @assert length(y) == N
    @assert N > M "system must be overdetermined"

    # Solve the first non-trivial M equations, skipping rows of zeroes.
    n = 1
    if M != 0
        @inbounds while n ≤ N && all(iszero.(A[n, :]))
            n += 1
        end
    end
    m = n + M - 1
    @assert m <= N
    @inbounds Am = SMatrix{M,M}(view(A, n:m, :))
    @inbounds ym = SVector{M}(view(y, n:m))
    x = Am \ ym

    # Check that equations (m+1):N are satisfied.
    @inbounds for i = (m + 1):N
        res = zero(eltype(x))
        for j in axes(A, 2)
            res += A[i, j] * x[j]
        end
        res ≈ y[i] || throw(NoUniqueSolutionError())
    end

    x
end

function checkdims_mul(y, A, x)
    M, N = size(A)
    if length(y) != M
        throw(DimensionMismatch("first dimension of A, $M, " *
                                "does not match length of y, $(length(y))"))
    end
    if length(x) != N
        throw(DimensionMismatch("second dimension of A, $N, " *
                                "does not match length of x, $(length(x))"))
    end
    nothing
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
using Test
using ForwardDiff2: D, DI, ForwardDiff2
using StaticArrays
using OrdinaryDiffEq, ForwardDiff
using LinearAlgebra
using ModelingToolkit: @variables, Variable

@testset "Derivative, Gradient, Jacobian, and Hessian" begin
    # Derivative
    @test D(x->1)(1) * true === false # zero
    @test D(sin)(1) * true === cos(1.0)
    @test D(float)(1) * 1 === 1
    @test D(AbstractFloat)(1) * 1 === 1
    @test D(x->[x, x^2])(3) * 1 == [1, 6]
    @test DI(DI(sin))(1) === -sin(1)
    @test D(D(sin))(1) * 1 === -sin(1)

    # Gradient
    @test D(x->1)([1,2,3]) * I == zeros(3)' # zero
    @test D(x->1)([1,2,3]) * zeros(3, 4) == zeros(4)'
    @test D(x->1)([1,2,3]) * zeros(3, 2) == zeros(2)'
    @test D(x->1)([1,2,3]) * ones(3) === 0
    @test D(sum)([1,2,3]) * I == ones(3)'
    @test D(sum)([1,2,3]) * ones(3) === 3.0
    @test D(sum)([1,2,3]) * ones(Int, 3) === 3
    @test D(sum)([1,2,3]) * ones(Int, 3, 2) == fill(3, 2)'
    @test D(sum)([1,2,3]) * ones(Int, 3, 4) == fill(3, 4)'

    # Jacobian
    zero_J = D(x->[1,1,1])([1,2,3]) * I # zero
    @test eltype(zero_J) === Int
    @test zero_J == zeros(Int, 3, 3)
    @test D(x->@SVector([x[1]^x[2], x[3]^3, x[3]*x[2]*x[1]]))(@SVector[1,2,3.]) * I === @SMatrix [2.0 0 0; 0 0 27; 6 3 2]
    dcumsum = D(cumsum)
    j2 = dcumsum(@SVector([1,2,3])) * I
    tj2 = @MMatrix [1 0 0; 1 1 0; 1 1 1]
    @test j2 == tj2
    @test typeof(j2) == typeof(tj2)
    @test dcumsum([1,2,3]) * I == [1 0 0; 1 1 0; 1 1 1]
    jvp1 = dcumsum([1,2,3]) * [1, 2, 3]
    @test jvp1 isa Vector{Int}
    @test jvp1 == [1, 3, 6]
    jvp2 = dcumsum([1,2,3]) * @SVector [1, 2, 3]
    @test jvp2 isa MVector{3,Int}
    @test jvp2 == [1, 3, 6]
    jmp3 = dcumsum([1,2,3]) * @SArray [1 2; 3 4; 5 6]
    @test jmp3 isa MMatrix{3,2,Int}
    @test jmp3 == [1 2; 4 6; 9 12]
    @test_throws ArgumentError dcumsum([1,2,3]) * [1, 2]
    @test D(x->@SVector([x[1], x[2]]))(@SVector([1,2,3])) * I === @SMatrix [1 0 0; 0 1 0]
    @test D(x->@SVector([x[1], x[2]]))(@SVector([1,2,3])) * @SVector([1, 2, 3]) === @SVector [1, 2]

    # Hessian
    dh1 = D(DI(x->x[1]^x[2] + x[3]^3 + x[3]*x[2]*x[1]))
    dh2 = D(D(x->x[1]^x[2] + x[3]^3 + x[3]*x[2]*x[1]))
    for dh in (dh1, dh2)
        @test dh(@SVector[1,2,3]) * I === @SMatrix [2 4 2; 4 0 1; 2 1 18.0]
        @test dh([1,2,3]) * I == [2 4 2; 4 0 1; 2 1 18.0]
        @test dh([1,2,3]) * @SVector([1, 2, 3]) == @SVector [16, 7, 58.0]
    end
end

@testset "Inference" begin
    simple_D = @inferred D(x->exp(x) + x^x + cos(x) + tan(x) + 2^x)
    @inferred simple_D(1) * 1
    @inferred simple_D(1) * true

    complicated_D = D(x->exp(x) + x^x + cos(x) + tan(x) + 2^x + log(cos(x)) + sec(pi*x) - angle(x) + one(x) / log1p(sin(x)))
    # broken due to `Core._apply`
    @test_broken @inferred complicated_D(1) * 1
end

@testset "Differentiate an ODE solver" begin
    # TODO: drop derivative API
    OrdinaryDiffEq.constvalue(t::Type{<:ForwardDiff2.Dual}) = ForwardDiff2.value(t)
    function lorenz(du, u, p, t)
        du[1] = 10.0 * (u[2] - u[1])
        du[2] = u[1] * (28.0 - u[3]) - u[2]
        du[3] = u[1] * u[2] - (8 / 3) * u[3]
        return nothing
    end
    function odesolver(u0, value, tf)
        tspan = (0.0, tf)
        prob = ODEProblem(lorenz, u0, tspan)
        sol = solve(prob, Tsit5(), internalnorm=(u, t) -> value(norm(u)))
        return sol[end]
    end

    u0 = [1.0, 0.0, 0.0]
    for tf in [1, 10, 100.0]
        fd2_D = D(x->odesolver(x, ForwardDiff2.value, tf))(u0)
        fd2 = fd2_D * I
        fd2_jvp = fd2_D * [1, 2, 3.0]
        fd = ForwardDiff.jacobian(x->odesolver(x, ForwardDiff.value, tf), u0)
        @test fd2 ≈ fd rtol=1e-12 # Wow, chaos!
        @test fd2_jvp ≈ fd * [1, 2, 3.0] rtol=1e-12 # Wow, chaos!
    end
end

@testset "No perturbation confusion" begin
    @test D(x -> x * D(y -> x + y)(1))(1) * 1 === 1
end

@testset "Multivariable function with closure" begin
    x = [3, 5.0]
    p = [0.1, 0.2]
    g(x, p) = [sin(x[1])*p[2], tan(p[1]*x[1])]

    num_fd2_D = ForwardDiff2.D(p->(ForwardDiff2.D(x->g(x, p))(x) * I))(p)
    num_fd2 = num_fd2_D * I
    num_fd2_jvp = ForwardDiff2.D(p->(ForwardDiff2.D(x->g(x, p))(x) * I))(p) * [1, 2.0]
    num_fd = ForwardDiff.jacobian(p->ForwardDiff.jacobian(x->g(x, p), x), [0.1, 0.2])
    @test num_fd2 == num_fd
    @test num_fd2_jvp == num_fd * [1, 2.0]

    @variables x[1:2] p[1:2] v[1:2]
    g(x, p) = [sin(x[1])*p[2], tan(p[1]*x[1])]
    @test_nowarn ForwardDiff2.D(p->(ForwardDiff2.D(x->g(x, p))(x) * I))(p) * I
    @test_nowarn ForwardDiff2.D(p->(ForwardDiff2.D(x->g(x, p))(x) * I))(p) * v
    @test_nowarn ForwardDiff2.D(p->(ForwardDiff2.D(x->g(x, p))(x) * I))(p) * [1, 2]
    @test_nowarn ForwardDiff2.D(p->(ForwardDiff2.D(x->g(x, p))(x) * I))(p) * [1, 2.0]

    function fun(vars, theta)
        vars[1] + exp(vars[2])
    end
    @test all(iszero, DI(vars_arg -> DI(theta_arg -> fun(vars_arg, theta_arg))(ones(5)))(ones(6)))
end


text:
#Script to load various SI values which are needed to convert back from atomic rydberg units
#23/07/2017

finestrucconst = 1/137 #Fine structure constant
atomenergy = 27.2 #Energy units in eV
Ry = 1
rpceV = 6.58211889e-16 #Reduced Planck constant in eV s
rpcJ = 1.054571596e-34 #Reduced Planck constant in J s
lightc = 2.99792458e8 #Speed of light in m/s
vacpmtvty = 8.854187817e-12 #Vacuum permitivity in F/m
boltzconst = 1.3807e-23 #Boltzmann's constant in J/K
bohrrad = 0.5291772083e-10 #Bohr radius
eleccharge = 1.602176462e-19 #electron charge
espin = 0.5 #electron spin

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
# ------------------------------------------------------------------
# Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    iqsim(trainimg::AbstractArray{T,N}, tilesize::Dims{N},
          simsize::Dims{N}=size(trainimg);
          overlap::NTuple{N,Float64}=ntuple(i->1/6,N),
          soft::AbstractVector=[], hard::HardData=HardData(), tol::Real=.1,
          cut::Symbol=:boykov, path::Symbol=:raster, nreal::Integer=1,
          threads::Integer=cpucores(), gpu::Bool=false,
          debug::Bool=false, showprogress::Bool=false)

Performs image quilting simulation as described in Hoffimann et al. 2017.

## Parameters

### Required

* `trainimg` is any 3D array (add ghost dimension for 2D)
* `tilesize` is the tile size (or pattern size)

### Optional

* `simsize` is the size of the simulation grid (default to training image size)
* `overlap` is the percentage of overlap (default to 1/6 of tile size)
* `soft` is a vector of `(data,dataTI)` pairs (default to none)
* `hard` is an instance of `HardData` (default to none)
* `tol` is the initial relaxation tolerance in (0,1] (default to .1)
* `cut` is the cut algorithm (`:dijkstra` or `:boykov`)
* `path` is the simulation path (`:raster`, `:dilation` or `:random`)
* `nreal` is the number of realizations (default to 1)
* `threads` is the number of threads for the FFT (default to all CPU cores)
* `gpu` informs whether to use the GPU or the CPU (default to false)
* `debug` informs whether to export or not the boundary cuts and voxel reuse
* `showprogress` informs whether to show or not estimated time duration

The main output `reals` consists of a list of 3D realizations that can be indexed with
`reals[1], reals[2], ..., reals[nreal]`. If `debug=true`, additional output is generated:

```julia
reals, cuts, voxs = iqsim(..., debug=true)
```

`cuts[i]` is the boundary cut for `reals[i]` and `voxs[i]` is the associated voxel reuse.
"""
function iqsim(trainimg::AbstractArray{T,N}, tilesize::Dims{N},
               simsize::Dims{N}=size(trainimg);
               overlap::NTuple{N,Float64}=ntuple(i->1/6,N),
               soft::AbstractVector=[], hard::HardData=HardData(), tol::Real=.1,
               cut::Symbol=:boykov, path::Symbol=:raster, nreal::Integer=1,
               threads::Integer=cpucores(), gpu::Bool=false,
               debug::Bool=false, showprogress::Bool=false) where {T,N}

  # number of threads in FFTW
  set_num_threads(threads)

  # sanity checks
  @assert ndims(trainimg) == 3 "image is not 3D (add ghost dimension for 2D)"
  @assert all(0 .< tilesize .≤ size(trainimg)) "invalid tile size"
  @assert all(simsize .≥ tilesize) "invalid grid size"
  @assert all(0 .< overlap .< 1) "overlaps must be in range (0,1)"
  @assert 0 < tol ≤ 1 "tolerance must be in range (0,1]"
  @assert cut ∈ [:dijkstra,:boykov] "invalid cut algorithm"
  @assert path ∈ [:raster,:dilation,:random] "invalid simulation path"
  @assert nreal > 0 "invalid number of realizations"

  # soft data checks
  if !isempty(soft)
    for (aux, auxTI) in soft
      @assert ndims(aux) == 3 "soft data is not 3D (add ghost dimension for 2D)"
      @assert all(size(aux) .≥ simsize) "soft data size < grid size"
      @assert size(auxTI) == size(trainimg) "auxiliary TI must have the same size as TI"
    end
  end

  # hard data checks
  if !isempty(hard)
    coordinates = [coord[i] for i in 1:N, coord in coords(hard)]
    @assert all(maximum(coordinates, dims=2) .≤ simsize) "hard data coordinates outside of grid"
    @assert all(minimum(coordinates, dims=2) .> 0) "hard data coordinates must be positive indices"
  end

  # calculate the overlap size from given percentage
  ovlsize = @. ceil(Int, overlap*tilesize)

  # spacing in raster path
  spacing = @. tilesize - ovlsize

  # calculate the number of tiles from grid size
  ntiles = @. ceil(Int, simsize / max(spacing, 1))

  # simulation grid dimensions
  padsize = @. ntiles*(tilesize - ovlsize) + ovlsize

  # training image dimensions
  TIsize = size(trainimg)

  # total overlap volume in simulation grid
  ovlvol = prod(padsize) - prod(@. padsize - (ntiles - 1)*ovlsize)

  # warn in case of 1-voxel overlaps
  if any((tilesize .>  1) .& (ovlsize .== 1))
    @warn "Overlaps with only 1 voxel, check tilesize/overlap configuration"
  end

  # always work with floating point
  TI = Float64.(trainimg)

  # inactive voxels in the training image
  NaNTI = isnan.(TI); TI[NaNTI] .= 0

  # disable tiles in the training image if they contain inactive voxels
  disabled = Vector{CartesianIndex{N}}()
  for ind in findall(NaNTI)
    start  = @. max(ind.I - tilesize + 1, 1)
    finish = @. min(ind.I, TIsize - tilesize + 1)
    tile   = CartesianIndices(ntuple(i -> start[i]:finish[i], N))
    append!(disabled, tile)
  end

  # keep track of hard data and inactive voxels
  datum = Vector{CartesianIndex{N}}()
  skipped = Set{CartesianIndex{N}}()
  if !isempty(hard)
    # hard data in grid format
    hardgrid = zeros(padsize)
    preset = falses(padsize)
    activated = trues(padsize)
    for coord in coords(hard)
      if isnan(hard[coord])
        activated[coord] = false
      else
        hardgrid[coord] = hard[coord]
        preset[coord] = true
      end
    end

    # deactivate voxels beyond true grid size
    ax = axes(activated)
    for d=1:N
      slice = ntuple(i -> i == d ? (simsize[d]+1:padsize[d]) : ax[i], N)
      activated[CartesianIndices(slice)] .= false
    end

    # grid must contain active voxels
    any_activated = any(activated[CartesianIndices(simsize)])
    @assert any_activated "simulation grid has no active voxel"

    # determine tiles that should be skipped and tiles with data
    for tileind in CartesianIndices(ntiles)
      # tile corners are given by start and finish
      start  = @. (tileind.I - 1)*spacing + 1
      finish = @. start + tilesize - 1
      tile   = CartesianIndices(ntuple(i -> start[i]:finish[i], N))

      if all(.!activated[tile])
        push!(skipped, tileind)
      else
        if any(preset[tile])
          push!(datum, tileind)
        end
      end
    end
  end

  # preprocess soft data
  softgrid = Vector{Array{Float64,N}}()
  softTI   = Vector{Array{Float64,N}}()
  if !isempty(soft)
    for (aux, auxTI) in soft
      auxpad = padsize .- min.(padsize, size(aux))

      AUX = padarray(aux, Pad(:symmetric, ntuple(i->0,N), auxpad))
      AUX = parent(AUX)
      AUX[isnan.(AUX)] .= 0

      AUXTI = copy(auxTI)
      AUXTI[NaNTI] .= 0

      # always work with floating point
      AUX   = Float64.(AUX)
      AUXTI = Float64.(AUXTI)

      push!(softgrid, AUX)
      push!(softTI, AUXTI)
    end
  end

  # overwrite path option if data is available
  !isempty(datum) && (path = :datum)

  # select cut algorithm
  boundary_cut = cut == :dijkstra ? dijkstra_cut : boykov_kolmogorov_cut

  # main output is a vector of 3D grids
  realizations = Vector{Array{Float64,N}}()

  # for each realization we have:
  boundarycuts = Vector{Array{Float64,N}}()
  voxelreuse   = Vector{Float64}()

  # show progress and estimated time duration
  showprogress && (progress = Progress(nreal))

  # preallocate memory for distance calculations
  distance = Array{Float64}(undef, TIsize .- tilesize .+ 1)

  # preallocate memory for cut mask
  mask = Array{Bool}(undef, tilesize)

  for real=1:nreal
    # allocate memory for current simulation
    simgrid = zeros(padsize)
    debug && (cutgrid = zeros(padsize))

    # keep track of pasted tiles
    pasted = Set{CartesianIndex{N}}()

    # construct simulation path
    simpath = genpath(ntiles, path, datum)

    # loop simulation grid tile by tile
    for ind in simpath
      tileind = lin2cart(ntiles, ind)

      # skip tile if all voxels are inactive
      tileind ∈ skipped && continue

      # if not skipped, proceed and paste tile
      push!(pasted, tileind)

      # tile corners are given by start and finish
      start  = @. (tileind.I - 1)*spacing + 1
      finish = @. start + tilesize - 1
      tile   = CartesianIndices(ntuple(i -> start[i]:finish[i], N))

      # current simulation dataevent
      simdev = view(simgrid, tile)

      # compute overlap distance
      distance .= 0
      for d=1:N
        # Cartesian index of previous and next tiles along dimension
        prev = CartesianIndex(ntuple(i -> i == d ? (tileind[d]-1) : tileind[i], N))
        next = CartesianIndex(ntuple(i -> i == d ? (tileind[d]+1) : tileind[i], N))

        # compute overlap distance with previous tile
        if ovlsize[d] > 1 && prev ∈ pasted
          oslice = ntuple(i -> i == d ? (1:ovlsize[d]) : (1:tilesize[i]), N)
          ovl = view(simdev, CartesianIndices(oslice))

          D = convdist(TI, ovl)

          ax = axes(D)
          dslice = ntuple(i -> i == d ? (1:TIsize[d]-tilesize[d]+1) : ax[i], N)
          distance .+= view(D, CartesianIndices(dslice))
        end

        # compute overlap distance with next tile
        if ovlsize[d] > 1 && next ∈ pasted
          oslice = ntuple(i -> i == d ? (spacing[d]+1:tilesize[d]) : (1:tilesize[i]), N)
          ovl = view(simdev, CartesianIndices(oslice))

          D = convdist(TI, ovl)

          ax = axes(D)
          dslice = ntuple(i -> i == d ? (spacing[d]+1:TIsize[d]-ovlsize[d]+1) : ax[i], N)
          distance .+= view(D, CartesianIndices(dslice))
        end
      end

      # disable dataevents that contain inactive voxels
      distance[disabled] .= Inf

      # compute hard and soft distances
      auxdistances = Vector{Array{Float64,N}}()
      if !isempty(hard) && any(preset[tile])
        harddev = view(hardgrid, tile)
        D = convdist(TI, harddev, weights=preset[tile])

        # disable dataevents that contain inactive voxels
        D[disabled] .= Inf

        # swap overlap and hard distances
        push!(auxdistances, distance)
        distance = D
      end
      for n=1:length(softTI)
        softdev = view(softgrid[n], tile)
        D = convdist(softTI[n], softdev)

        # disable dataevents that contain inactive voxels
        D[disabled] .= Inf

        push!(auxdistances, D)
      end

      # current pattern database
      patterndb = isempty(auxdistances) ? findall(vec(distance .≤ (1+tol)minimum(distance))) :
                                          relaxation(distance, auxdistances, tol)

      # pattern probability
      patternprobs = tau_model(patterndb, distance, auxdistances)

      # pick a pattern at random from the database
      rind = sample(patterndb, weights(patternprobs))
      start  = lin2cart(size(distance), rind)
      finish = @. start.I + tilesize - 1
      rtile  = CartesianIndices(ntuple(i -> start[i]:finish[i], N))

      # selected training image dataevent
      TIdev = view(TI, rtile)

      # boundary cut mask
      mask .= false
      for d=1:N
        # Cartesian index of previous and next tiles along dimension
        prev = CartesianIndex(ntuple(i -> i == d ? (tileind[d]-1) : tileind[i], N))
        next = CartesianIndex(ntuple(i -> i == d ? (tileind[d]+1) : tileind[i], N))

        # compute mask with previous tile
        if ovlsize[d] > 1 && prev ∈ pasted
          oslice = ntuple(i -> i == d ? (1:ovlsize[d]) : (1:tilesize[i]), N)
          inds = CartesianIndices(oslice)
          A = view(simdev, inds); B = view(TIdev, inds)
          mask[inds] .|= boundary_cut(A, B, d)
        end

        # compute mask with next tile
        if ovlsize[d] > 1 && next ∈ pasted
          oslice = ntuple(i -> i == d ? (spacing[d]+1:tilesize[d]) : (1:tilesize[i]), N)
          inds = CartesianIndices(oslice)
          A = view(simdev, inds); B = view(TIdev, inds)
          mask[inds] .|= reverse(boundary_cut(reverse(A, dims=d), reverse(B, dims=d), d), dims=d)
        end
      end

      # paste quilted pattern from training image
      simdev[.!mask] = TIdev[.!mask]

      # save boundary cut
      debug && (cutgrid[tile] = mask)
    end

    # save voxel reuse
    debug && push!(voxelreuse, sum(cutgrid)/ovlvol)

    # hard data and shape correction
    if !isempty(hard)
      simgrid[preset] = hardgrid[preset]
      simgrid[.!activated] .= NaN
      debug && (cutgrid[.!activated] .= NaN)
    end

    # throw away voxels that are outside of the grid
    simgrid = view(simgrid, CartesianIndices(simsize))
    debug && (cutgrid = view(cutgrid, CartesianIndices(simsize)))

    # save and continue
    push!(realizations, simgrid)
    debug && push!(boundarycuts, cutgrid)

    # update progress bar
    showprogress && next!(progress)
  end

  debug ? (realizations, boundarycuts, voxelreuse) : realizations
end