text:
#This script illustate how to set up constraints and project a 2D model onto an intersection
# of an intersection and a Minkowski sum of intersections of sets
# Bas Peters, 2019

using Distributed
@everywhere using SetIntersectionProjection
using MAT
using LinearAlgebra
ENV["MPLBACKEND"]="qt5agg"
using PyPlot

mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

#PARSDMM options:
options=PARSDMM_options()
options.FL=Float32
#options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 500

options.feas_tol     = 0.001
options.obj_tol      = 0.001
options.evol_rel_tol = 0.00001
set_zero_subnormals(true)
BLAS.set_num_threads(2)

#select working precision
if options.FL==Float64
  TF = Float64
  TI = Int64
elseif options.FL==Float32
  TF = Float32
  TI = Int32
end

#load image to project
file = matopen(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Data/compass_velocity.mat"))
m    = read(file, "Data");close(file)
m    = m[1:341,200:600]
m    = permutedims(m,[2,1])

#set up computational grid (25 and 6 m are the original distances between grid points)
comp_grid = compgrid((TF(25.0), TF(6.0)),(size(m,1), size(m,2)))
m         = convert(Vector{TF},vec(m))

#define axis limits and colorbar limits
xmax = comp_grid.d[1]*comp_grid.n[1]
zmax = comp_grid.d[2]*comp_grid.n[2]
vmi  = 1500.0
vma  = 4500.0

#####################################################################################################
# approximate monotonically increasing velocity with depth + TV decomposition #######################

#constraints for total model (sum of two compontents)
constraint_sum = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
m_min     = 1500.0
m_max     = 4500.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_sum, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#constraints for component 1 of the model
constraint_c1 = Vector{SetIntersectionProjection.set_definitions}()

#bounds component 1:
m_min     = 1500f0.*ones(TF,comp_grid.n)
m_max     = 4500f0.*ones(TF,comp_grid.n)
m_max[:,1:35] .= 1500f0
m_min     = vec(m_min)
m_max     = vec(m_max)
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c1, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#slope constraints component 1(vertical)
m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c1, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#constraints for component 2 of the model
constraint_c2 = Vector{SetIntersectionProjection.set_definitions}()

#bounds component 2:
m_min     = -1500.0
m_max     = 1500.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c2, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#TV constraint component 2
m_min     = 0
(TV,dummy1,dummy2,dummy3) = get_TD_operator(comp_grid,"TV",options.FL)
m_max = 0.15f0*norm(TV*vec(m),1)
set_type  = "l1"
TD_OP     = "TV"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c2, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#get all projectors for sub-problems and linear operators
(P_sub_sum,TD_OP_sum,set_Prop_sum) = setup_constraints(constraint_sum,comp_grid,options.FL)
(P_sub_c1,TD_OP_c1,set_Prop_c1)    = setup_constraints(constraint_c1,comp_grid,options.FL)
(P_sub_c2,TD_OP_c2,set_Prop_c2)    = setup_constraints(constraint_c2,comp_grid,options.FL)

#pre-computations
options.parallel  = false
options.Minkowski = true
(TD_OP,set_Prop,AtA,l,y) = PARSDMM_precompute_distribute_Minkowski(TD_OP_c1,TD_OP_c2,TD_OP_sum,set_Prop_c1,set_Prop_c2,set_Prop_sum,comp_grid,options)

#put all projectors for sub-problems in a vector
P_sub = deepcopy(P_sub_c1);
append!(P_sub,P_sub_c2)
append!(P_sub,P_sub_sum)

options.zero_ini_guess=true

#project
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

#plot; output is component 1 (m1), component 2 (m2) stacked in a vector
m1 = x[1:length(m)]
m2 = x[1+length(m):end]
ms = m1+m2; #total model output

figure(figsize=(10.5,1.75));
imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]);colorbar();title("Original")
savefig("original_model.png",bbox_inches="tight")
figure(figsize=(10.5,1.75));imshow(permutedims(reshape(m1,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]);colorbar(); title("Projection (component 1, positive vertical slope)")
savefig("projected_model_component1.png",bbox_inches="tight")
figure(figsize=(10.5,1.75));imshow(permutedims(reshape(m2,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=-1500,vmax=1500,extent=[0,  xmax, zmax, 0]);colorbar(); title("Projection (component 2, low TV)")
savefig("projected_model_component2.png",bbox_inches="tight")
figure(figsize=(10.5,1.75));imshow(permutedims(reshape(ms,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]);colorbar(); title("Projection (sum of components)")
savefig("projected_model_component_sum.png",bbox_inches="tight")

figure();
subplot(3, 3, 3);semilogy(log_PARSDMM.r_pri)          ;title("r primal")
subplot(3, 3, 4);semilogy(log_PARSDMM.r_dual)         ;title("r dual")
subplot(3, 3, 1);semilogy(log_PARSDMM.obj)            ;title(L"$ \frac{1}{2} || \mathbf{m}-\mathbf{x} ||_2^2 $")
subplot(3, 3, 2);semilogy(log_PARSDMM.set_feasibility);title("TD feasibility violation")
subplot(3, 3, 5);plot(log_PARSDMM.cg_it)              ;title("nr. of CG iterations")
subplot(3, 3, 6);semilogy(log_PARSDMM.cg_relres)      ;title("CG rel. res.")
subplot(3, 3, 7);semilogy(log_PARSDMM.rho)            ;title("rho")
subplot(3, 3, 8);plot(log_PARSDMM.gamma)              ;title("gamma")
subplot(3, 3, 9);semilogy(log_PARSDMM.evol_x)         ;title("x evolution")
tight_layout()
#tight_layout(pad=0.0, w_pad=0.0, h_pad=1.0)
savefig("PARSDMM_logs.png",bbox_inches="tight")

#################################################################################################################################
# approximate monotonically increasing velocity with depth + monotonically decreasing velocity with depth #######################

#constraints for total model
constraint_sum = Vector{SetIntersectionProjection.set_definitions}()

#bounds:
m_min     = 1500.0
m_max     = 4500.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_sum, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#constraints for component 1 of the model
constraint_c1 = Vector{SetIntersectionProjection.set_definitions}()

#bounds component 1:
m_min     = 1500.0f0*ones(TF,prod(comp_grid.n))
m_max     = 4500.0f0.*ones(TF,comp_grid.n)
m_max[:,1:35] .= 1500.0f0
m_max     = vec(m_max)
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c1, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#slope constraints component 1(vertical, positive)
m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c1, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#constraints for component 2 of the model
constraint_c2 = Vector{SetIntersectionProjection.set_definitions}()

#bounds component 2:
m_min     = -3000.0
m_max     = 0.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c2, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#slope constraints component 1(vertical, negative)
m_min     = -1e6
m_max     = 0.0
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint_c2, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

(P_sub_sum,TD_OP_sum,set_Prop_sum) = setup_constraints(constraint_sum,comp_grid,options.FL)
(P_sub_c1,TD_OP_c1,set_Prop_c1) = setup_constraints(constraint_c1,comp_grid,options.FL)
(P_sub_c2,TD_OP_c2,set_Prop_c2) = setup_constraints(constraint_c2,comp_grid,options.FL)

options.parallel  = false
options.Minkowski = true
(TD_OP,set_Prop,AtA,l,y) = PARSDMM_precompute_distribute_Minkowski(TD_OP_c1,TD_OP_c2,TD_OP_sum,set_Prop_c1,set_Prop_c2,set_Prop_sum,comp_grid,options)

P_sub = deepcopy(P_sub_c1);
append!(P_sub,P_sub_c2)
append!(P_sub,P_sub_sum)

options.zero_ini_guess=true

@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
@time (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

#plot
m1 = x[1:length(m)]
m2 = x[1+length(m):end]
ms = m1+m2;

figure(figsize=(10.5,1.75));imshow(permutedims(reshape(m,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]);colorbar();title("Original")
savefig("original_model.png",bbox_inches="tight")
figure(figsize=(10.5,1.75));imshow(permutedims(reshape(m1,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]);colorbar(); title("Projection (component 1, positive vertical slope)")
savefig("projected_model_component1.png",bbox_inches="tight")
figure(figsize=(10.5,1.75));imshow(permutedims(reshape(m2,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=-1500,vmax=1500,extent=[0,  xmax, zmax, 0]);colorbar(); title("Projection (component 2, negative vertical slope)")
savefig("projected_model_component2.png",bbox_inches="tight")
figure(figsize=(10.5,1.75));imshow(permutedims(reshape(ms,(comp_grid.n[1],comp_grid.n[2])),[2,1]),cmap="jet",vmin=vmi,vmax=vma,extent=[0,  xmax, zmax, 0]);colorbar(); title("Projection (sum of components)")
savefig("projected_model_component_sum.png",bbox_inches="tight")

figure();
subplot(3, 3, 3);semilogy(log_PARSDMM.r_pri)          ;title("r primal")
subplot(3, 3, 4);semilogy(log_PARSDMM.r_dual)         ;title("r dual")
subplot(3, 3, 1);semilogy(log_PARSDMM.obj)            ;title(L"$ \frac{1}{2} || \mathbf{m}-\mathbf{x} ||_2^2 $")
subplot(3, 3, 2);semilogy(log_PARSDMM.set_feasibility);title("TD feasibility violation")
subplot(3, 3, 5);plot(log_PARSDMM.cg_it)              ;title("nr. of CG iterations")
subplot(3, 3, 6);semilogy(log_PARSDMM.cg_relres)      ;title("CG rel. res.")
subplot(3, 3, 7);semilogy(log_PARSDMM.rho)            ;title("rho")
subplot(3, 3, 8);plot(log_PARSDMM.gamma)              ;title("gamma")
subplot(3, 3, 9);semilogy(log_PARSDMM.evol_x)         ;title("x evolution")
tight_layout()
#tight_layout(pad=0.0, w_pad=0.0, h_pad=1.0)
savefig("PARSDMM_logs.png",bbox_inches="tight")

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Notes:
#   - import raw data, initial cleaning
#   - subset and compute data for figure
#   - loop over months and add to DataFrame
#   - export data
#   - import data and plot figures (different script)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# timespan
years = string.(2001:2021)
months = lpad.(1:12, 2, '0')

# initialize dataframes
df_fig1 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[],
                    VALUE_IN_EUROS=Float64[], QUANTITY_IN_KG=Float64[], PRICE_INDEX=Float64[])

df_fig2 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[], PRODUCT_NC_digits=String[],
                    VALUE_IN_EUROS=Float64[], QUANTITY_IN_KG=Float64[], PRICE_INDEX=Float64[])

df_fig3 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[],
                    VALUE_SHARE=Float64[], QUANTITY_SHARE=Float64[])

df_tab1 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], YEAR=String[], PRODUCT_NC_digits=String[],
                    VALUE_SHARE_PARTNER=Float64[], QUANTITY_SHARE_PARTNER=Float64[], VALUE_SHARE_TOTAL=Float64[], QUANTITY_SHARE_TOTAL=Float64[])

df_tab1 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PRODUCT_NC_digits=String[], PERIOD=Int64[],
                    VALUE_IN_EUROS=Float64[], QUANTITY_IN_KG=Float64[])

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

for i in years
    for j in months

        #import and clean data
        df = initial_cleaning(i, j)

        # create EU as partner
        df = append_EU(df, ctrys_EU27)

        # append data for figures/tables
        append!(df_fig1, data_fig1(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"]))
        append!(df_fig2, data_fig2(df, ["BE"], ["GB"], 2))
        append!(df_fig3, data_fig3(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"]))
        append!(df_tab1, tab1(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"], 2))

        println(" ✓ Data for $j/$i has been successfully added. \n")

    end
end

# export locally
#CSV.write(dir_io * "clean/" * "df_fig1" * ".csv", df_fig1)
#CSV.write(dir_io * "clean/" * "df_fig2" * ".csv", df_fig2)
#CSV.write(dir_io * "clean/" * "df_fig3" * ".csv", df_fig3)
CSV.write(dir_io * "clean/" * "df_tab1" * ".csv", df_tab1)

# export Dropbox
CSV.write(dir_dropbox * "results/" * "df_fig1" * ".csv", df_fig1)
CSV.write(dir_dropbox * "results/" * "df_fig2" * ".csv", df_fig2)
CSV.write(dir_dropbox * "results/" * "df_fig3" * ".csv", df_fig3)
CSV.write(dir_dropbox * "results/" * "df_tab1" * ".csv", df_tab1)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

ext: jl


text:
# SPDX-License-Identifier: MIT

@testset "logaddexp" begin
    @test (@inferred MarkovModels.Semirings.logaddexp(2.0, 3.0)) ≈ log(exp(2.0) + exp(3.0))
    @test (@inferred MarkovModels.Semirings.logaddexp(10002.0, 10003.0)) ≈
        10000 + MarkovModels.Semirings.logaddexp(2.0, 3.0)
end

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
abstract type Metric end

struct MetricAccumulator{T <: Metric}
  val :: Float64
  cnt :: Int
end

function merge{T <: Metric}(ms :: Vector{MetricAccumulator{T}})
  val = sum(m.val for m in ms)
  cnt = sum(m.cnt for m in ms)
  return MetricAccumulator{T}(val, cnt)
end

function Base.show{T <: Metric}(io :: IO, m :: MetricAccumulator{T})
  print(io, @sprintf("%s = %8.4f", name(T), m.val / m.cnt))
end

type Accuracy <: Metric
end

"""
`evaluate` takes a real-valued prediction vector and a +1/-1 integer vector of labels.
It returns a metric accumulator.
"""
function evaluate(:: Type{Accuracy}, pred :: Vector{Float64}, y :: Vector{Int})
  cnt = length(y)
  acc = 0.0
  for (a, b) in zip(pred, y)
    acc += (a*b > 0)
  end
  return MetricAccumulator{Accuracy}(acc, cnt)
end

function name(:: Type{Accuracy})
  return "ACC"
end


text:
@graphs begin
	generic_move = """
	Move according to the latest statement. Select army at random if several can perform such a move.
	"""

	settlement_cond = """
	An army can move into an emtpy settlement of a nation at war.
	"""
	settlement_cond_2 = """
	An army can move into an empty settlement of a nation at war, without increasing the distance to its *target* (the *target* may change).
	"""
	settlement_move_unit = """
	Move 1 unit into an empty settlement of nation at war.

	Army priority: Random

	Unit priority:
	1. Regular
	2. Elite
	3. Random
	"""

	merge_cond = """
	Two armies, where at least one is not *mobile*, can merge.
	And, merging the armies would increase the number of *mobile* armies.
	Or, the merged army would have higher *value* then either of the two currently have.
	"""

	merge_move = """
	Merge two armies.

	Priority:
	1. Merge decrease greatest distance to *target* (*target* may change)
	2. Highest *value* of resulting army
	3. Moves the army furthest from its *target*
	4. Least number of units left behind after move
	5. Region where armies merge contains a stronghold
	6. Random
	"""

	move_target_cond = """
	A *mobile* army can move towards, or attack, its *target*.
	"""

	move_target = """
	Select a *mobile* army and move towards or attack *target*.

	Priority:
	1. Army is adjacent to its *target*
	2. Army's *target* is in a nation at war
	3. Move/attack doesn't activate a nation
	4. Move/attack doesn't change a nation to "at war"
	5. Army whose *target* is highest in the priority list in the definition of *target*
	6. Army with hightest *value*
	7. Move/attack does not block another *mobile* army's shortest route to their *target*
	8. Destination region contains the Fellowship
	9. Random
	"""

	basic_move_cond = """
	A Shadow army is on the board.
	"""

	basic_move = """
	Move an army.

	Priority:
	1. Move doesn't change a nation to "at war"
	2. Merge two armies to create the highest army *value* possible
	3. Army's *target* is adjacent to a passive Shadow army
	4. Movement ends adjacent to another Shadow army
	5. Decreases distance to *target* (*target* may change)
	6. Army with hightest *value*
	7. Random
	"""

	################################################################################
	@node movement_attack_besiege = Start() -> mv_1

	@node mv_1 = BinaryCondition("""
								 A *mobile* army is adjacent to *target* not under siege.
								 """) -> [n_true = mv_1_yes, n_false = mv_1_return]
	@node mv_1_return = ReturnFromGraph() -> []
	@node mv_1_yes = UseActiveDie() -> mv_1_action
	@node mv_1_action = PerformAction("""
									  Attack with army adjacent to *target* not under siege.

									  Priority:
									  1. Army whose *target* is in a nation at war
									  2. Army whose attack would not put a nation at war.
									  3. Army whose *target* is in an active nation
									  4. Highest *value* army
									  5. Random
									  """) -> mv_1_end
	@node mv_1_end = End() -> []

	################################################################################
	@node movement_attack_corr = Start() -> mv_2
	@node mv_2 = BinaryCondition("""
								 There are dice in the hunt pool.
								 And, no army is in the Fellowship's region
								 And, the Fellowship do not reach mordor with the current progress.
								 And, an army can move into the Fellowship's region without increasing the distance to its *target* (the *target* may change).
								 """) -> [n_true = mv_2_yes, n_false = mv_3]
	@node mv_2_yes = UseActiveDie() -> mv_2_action
	@node mv_2_action = PerformAction(generic_move) -> mv_2_army_die_to_move
	@node mv_2_army_die_to_move = CheckActiveDie('A') -> [n_true=mv_2_movement_remains, n_false=mv_2_end]
	@node mv_2_movement_remains = BinaryCondition("""
												  The Army Die have one move remaining.
												  """) -> [n_true = mv_2, n_false = mv_2_end]
	@node mv_2_end = End() -> []

	########################################
	@node mv_3 = BinaryCondition(settlement_cond) -> [n_true = mv_3_1, n_false = mv_4]
	@node mv_3_1 = BinaryCondition(settlement_cond_2) -> [n_true = mv_3_yes, n_false = mv_3_no]
	@node mv_3_yes = UseActiveDie() -> mv_3_yes_action
	@node mv_3_yes_action = PerformAction(generic_move) -> mv_3_movement_remains
	@node mv_3_no = UseActiveDie() -> mv_3_no_action
	@node mv_3_no_action = PerformAction(settlement_move_unit) -> mv_3_army_die_to_move
	@node mv_3_army_die_to_move = CheckActiveDie('A') -> [n_true=mv_3_movement_remains, n_false=mv_3_end]
	@node mv_3_movement_remains = BinaryCondition("""
												  The Army Die have one move remaining.
												  """) -> [n_true = mv_3, n_false = mv_3_end]
	@node mv_3_end = End() -> []

	########################################
	@node mv_4 = BinaryCondition(merge_cond) -> [n_true = mv_4_yes, n_false = mv_5]
	@node mv_4_yes = UseActiveDie() -> mv_4_action
	@node mv_4_action = PerformAction(merge_move) -> mv_4_army_die_to_move
	@node mv_4_army_die_to_move = CheckActiveDie('A') -> [n_true=mv_4_movement_remains, n_false=mv_4_end]
	@node mv_4_movement_remains = BinaryCondition("""
												  The Army Die have one move remaining.
												  """ ) -> [n_true = mv_4, n_false = mv_4_end]
	@node mv_4_end = End() -> []

	################################################################################
	@node movement_attack_basic = Start() -> mv_5
	@node mv_5 = BinaryCondition(move_target_cond) -> [n_true = mv_5_yes, n_false = mv_6]
	@node mv_5_yes = UseActiveDie() -> mv_5_action
	@node mv_5_action = PerformAction(move_target) -> mv_5_army_die_to_move
	@node mv_5_army_die_to_move = CheckActiveDie('A') -> [n_true=mv_5_movement_remains, n_false=mv_5_end]
	@node mv_5_movement_remains = BinaryCondition("""
												  The Army Die have one move remaining.
												  """) -> [n_true = mv_5, n_false = mv_5_end]
	@node mv_5_end = End() -> []

	########################################
	@node mv_6 = BinaryCondition(basic_move_cond) -> [n_true = mv_6_yes, n_false = mv_return]
	@node mv_return = ReturnFromGraph() -> []
	@node mv_6_yes = UseActiveDie() -> mv_6_action
	@node mv_6_action = PerformAction(basic_move) -> mv_6_army_die_to_move
	@node mv_6_army_die_to_move = CheckActiveDie('A') -> [n_true=mv_6_movement_remains, n_false=mv_6_end]
	@node mv_6_movement_remains = BinaryCondition("""
												  The Army Die have one move remaining.
												  """) -> [n_true = mv_6, n_false = mv_6_end]
	@node mv_6_end = End() -> []

	################################################################################
	@node movement_attack_card = Start() -> mv_c_1

	@node mv_c_1 = BinaryCondition(settlement_cond) -> [n_true = mv_c_1_1, n_false = mv_c_2]
	@node mv_c_1_1 = BinaryCondition(settlement_cond_2) -> [n_true = mv_c_1_army, n_false = mv_c_1_unit]
	@node mv_c_1_army = PerformAction(generic_move) -> mv_c_1_army_end
	@node mv_c_1_army_end = End() -> []

	@node mv_c_1_unit = PerformAction(settlement_move_unit) -> mv_c_1_unit_end
	@node mv_c_1_unit_end = End() -> []

	@node mv_c_2 = BinaryCondition(merge_cond) -> [n_true = mv_c_2_move, n_false = mv_c_3]
	@node mv_c_2_move = PerformAction(merge_move) -> mv_c_2_end
	@node mv_c_2_end = End() -> []

	@node mv_c_3 = BinaryCondition(move_target_cond) -> [n_true = mv_c_3_move, n_false = mv_c_4_move]
	@node mv_c_3_move = PerformAction(move_target) -> mv_c_3_end
	@node mv_c_3_end = End() -> []

	@node mv_c_4_move = PerformAction(basic_move) -> mv_c_4_end
	@node mv_c_4_end = End() -> []

end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
using SmolkaFormat
using Test

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


text:
struct LockandKeyLookup{ A <: Function, B <: Function, C <: Function, K, T }
    key::K
    tumbler::T
    key_lookup_fn::A
    pin_lookup_fn::B
    emission_fn::C
    key_length::Int
    tumblers::Vector{ Int }
end

#Sizes of the return collection cannot be known at run time! Base.SizeUnknown() is the secret sauce to make true Julia Generators!
Base.length(::LockandKeyLookup) = Base.SizeUnknown()
Base.IteratorSize(::LockandKeyLookup) = Base.SizeUnknown()

"""
    LockandKeyLookup(   key, tumbler,
                        key_lookup_fn, pin_lookup_fn,
                        emitter_fn = ( k, t ) -> key_lookup_fn( k ) == pin_lookup_fn( t ) )

Instantiates a `LockandKeyLookup` iterable object.
The `key` is the iterator to find a matching condition for across the `tumbler` iterators or internally called "pins".

the `_lookup_fn` functions are the functions used to `select` the data from each iteration for matching.
Ie: `key_lookup_fn`     = row -> row[!, [ :id, :time, :address ] ]
    `tumbler_lookup_fn` = row -> row[!, [ :row_id, :t, :Address ] ]

The `emitter_fn` is the function used to determine if there is infact a match, by default it asseses if the lookup functions equate.

"""
function LockandKeyLookup(  key, tumbler,
                            key_lookup_fn, pin_lookup_fn,
                            emitter_fn = ( k, t ) -> k == t)
    return LockandKeyLookup(    key, tumbler,
                                key_lookup_fn, pin_lookup_fn,
                                emitter_fn,
                                first( size( key ) ), length.(tumbler) )
end

"""
    get_smallest_pin( pin_values, not_nothing )::Int

An internal function for the tumbler portion of the iterator. It scans an array of available tumbler iterators - aka Pins - and finds the smallest index.

Note: In cases where there is a tie, the leftmost tied pin will be selected.
In cases where all pins are of type `Nothing` this function returns an integer value of 0.

"""
function get_smallest_pin( pin_values, not_nothing )::Int
    is_something    = sum( not_nothing )
    if is_something == 0
        return 0
    elseif is_something == 1
        return findfirst( not_nothing )
    elseif is_something > 1
        return findall( not_nothing )[ argmin( pin_values[ not_nothing ] ) ]
    end
end

"""
    unpack_2tuples( x )

This function unpacks collections of 2-tuples into 2 arrays in one pass.
Note: Nothings are handled by expanding into `(nothing, nothing)`

"""
unpack_2tuples( x ) = collect.(zip( [ isnothing( i ) ? [ i, i ] : i  for i in x ]... ) )

"""
    Base.iterate( lkl::LockandKeyLookup, state::T = ( iterate( lkl.key ), iterate.(lkl.tumbler) ) ) where T <: Any

This would be the main iterator for a given `LockandKeyLookup` instance.

"""
function Base.iterate( lkl::LockandKeyLookup, state::T = ( iterate( lkl.key ), iterate.(lkl.tumbler) ) ) where T <: Any
    # Get Latest Key:Tumbler State
    ( key_value, key_state ), tumbler_values_and_states = state
    tumbler_values_and_states = collect( tumbler_values_and_states )
    tumbler_values_and_states = convert(Vector{ Union{ eltype(tumbler_values_and_states), Nothing } }, tumbler_values_and_states)
    #( last( key_state ) > lkl.key_length ) && return nothing
    tumbler_values, tumbler_states = unpack_2tuples( tumbler_values_and_states )
    tumbler_values = convert(Vector{ Union{ eltype(tumbler_values), Nothing } }, tumbler_values)
    tumbler_states = convert(Vector{ Union{ eltype(tumbler_states), Nothing } }, tumbler_states)
    # Get Lowest Pin
    get_key               = lkl.key_lookup_fn( key_value )
    tumbler_values        = collect(tumbler_values)
    notnothing            = .!isnothing.( tumbler_values )
    get_pins              = map( x -> isnothing(x) ? x : lkl.pin_lookup_fn( x ), tumbler_values )
    smallest_pin          = get_smallest_pin( get_pins, notnothing )
    (smallest_pin == 0) && return nothing
    #pre-allocate output
    key_idx, tumbler_idx, pin_idx  = last(key_state), smallest_pin, last( tumbler_states[smallest_pin] )
    # If match( Key, Lowest Pin )
    if lkl.emission_fn( get_key, get_pins[ smallest_pin ] )
        tumbler_idx = smallest_pin
        pin_idx     = last( tumbler_states[ smallest_pin ] )
        values_and_states = Base.iterate( lkl.tumbler[ smallest_pin ], tumbler_states[ smallest_pin ] )
        tumbler_values[ smallest_pin ], tumbler_states[ smallest_pin ] = isnothing( values_and_states ) ? [nothing,nothing] : values_and_states
        #println("111 Kick tumbler match")
    else
        while ( !lkl.emission_fn( get_key, get_pins[ smallest_pin ] ) )
            if get_key < get_pins[ smallest_pin ]
                #println("222 Kick key")
                while ( get_key < get_pins[ smallest_pin ] )
                    key_value_and_state = Base.iterate( lkl.key, key_state )
                    isnothing( key_value_and_state ) && return nothing
                    ( key_value, key_state ) = key_value_and_state
                    get_key = lkl.key_lookup_fn( key_value )
                end
                key_idx = last( key_state )
            elseif get_key > get_pins[smallest_pin]
                #println("333 Kick tumbler")
                while ( get_key > get_pins[ smallest_pin ] )
                    values_and_states = Base.iterate( lkl.tumbler[ smallest_pin ], tumbler_states[ smallest_pin ] )
                    tumbler_values[ smallest_pin ], tumbler_states[ smallest_pin ] = isnothing( values_and_states ) ? [nothing,nothing] : values_and_states
                    notnothing            = .!isnothing.( tumbler_values )
                    get_pins              = map( x -> isnothing(x) ? x : lkl.pin_lookup_fn( x ), tumbler_values )
                    smallest_pin          = get_smallest_pin( get_pins, notnothing )
                    ( smallest_pin == 0 ) && return nothing
                end
                key_idx = last( key_state )
            end
            #Update the tumbler if we have a match
            tumbler_idx = smallest_pin
            pin_idx     = last( tumbler_states[ smallest_pin ] )
            if lkl.emission_fn( get_key, get_pins[ smallest_pin ] )
                values_and_states = Base.iterate( lkl.tumbler[ smallest_pin ], tumbler_states[ smallest_pin ] )
                tumbler_values[ smallest_pin ], tumbler_states[ smallest_pin ] = isnothing( values_and_states ) ? [nothing,nothing] : values_and_states
            end
        end
    end
    return  ( key_idx => ( tumbler_idx, pin_idx ) ), ( ( key_value, key_state ) , zip( tumbler_values, tumbler_states ) )
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Reset MDP.
"""
function reset!(mdp::ASTMDP)
    reset!(mdp.sim)
    reset!(mdp.reward.heuristic)
    return nothing
end

"""
Return set of available actions as sampleable object. # TODO: not type-stable.
"""
function actions(mdp::ASTMDP{<:State, SampleAction})
    env = environment(mdp.sim)
    return mdp.flatten ? Flat(env) : env
end

actions(::ASTMDP{<:State, SeedAction}) = UInt32

"""
Return observation of environment.
"""
observe(mdp::ASTMDP{ObservableState, <:Action}) = Float32.(observe(mdp.sim))

"""
Step simulation and return log probability of transition.
"""
function evaluate!(mdp::ASTMDP{<:State, A}, a::A) where A <: SampleAction
    env = environment(mdp.sim)
    val = a.sample
    logp = logprob(env, val, mdp.reward.marginalize)
    step!(mdp.sim, val)
    return logp
end

function evaluate!(mdp::ASTMDP{<:State, A}, a::A) where A <: SeedAction
    copy!(RNG_TEMP[], mdp.rng)
    Random.seed!(mdp.rng, a.seed)
    logp = step!(mdp.rng, mdp.sim)
    copy!(mdp.rng, RNG_TEMP[])
    return logp
end

"""
Apply raw action to environment and return reward.
"""
act!(mdp::ASTMDP{<:State, <:Action}, action) = act!(mdp, convert_a(mdp, action))

"""
Apply converted action to environment and return reward.
"""
function act!(mdp::ASTMDP{<:State, A}, a::A) where A <: Action
    # rewards (partial application)
    event = mdp.episodic ? isterminal(mdp.sim) && isevent(mdp.sim) : isevent(mdp.sim)
    bonus = event ? mdp.reward.event_bonus : 0.0
    heur = mdp.reward.heuristic(mdp)
    rew = custom_reward(mdp.sim, a)

    # stepping and likelihood evaluation
    logp = evaluate!(mdp, a)

    # final reward calculation
    return mdp.reward.reward_function(logp, bonus, heur(mdp)) + rew(mdp.sim)
end

"""
Return Boolean indicating termination status.
"""
terminated(mdp::ASTMDP) = isterminal(mdp.sim) || !mdp.episodic && isevent(mdp.sim)

# Connection between AdaStress internal functions and CommonRLInterface

CommonRLInterface.reset!(mdp::AbstractASTMDP) = reset!(mdp)

CommonRLInterface.actions(mdp::AbstractASTMDP) = actions(mdp)

CommonRLInterface.observe(mdp::AbstractASTMDP) = observe(mdp)

CommonRLInterface.act!(mdp::AbstractASTMDP, action) = act!(mdp, action)

CommonRLInterface.terminated(mdp::AbstractASTMDP) = terminated(mdp)


text:
"""
    initialize_sbm_gwf_model(config::Config)

Initial part of the sbm_gwf model concept. The model contains:
    - the vertical SBM concept
    - the following lateral components:
        - 1-D kinematic wave for river flow
        - 1-D kinematic wave for overland flow
        - unconfined aquifer with groundwater flow in four directions (adjacent cells)

The unconfined aquifer contains a recharge, river and a drain (optional) boundary.

The initial part reads the input settings and data as defined in the Config object.
Will return a Model that is ready to run.
"""
function initialize_sbm_gwf_model(config::Config)

    # unpack the paths to the NetCDF files
    tomldir = dirname(config)
    static_path = joinpath(tomldir, config.input.path_static)
    cyclic_path = joinpath(tomldir, config.input.path_static)
    dynamic_path = joinpath(tomldir, config.input.path_forcing)
    instate_path = joinpath(tomldir, config.state.path_input)
    output_path = joinpath(tomldir, config.output.path)

    Δt = Second(config.timestepsecs)
    sizeinmetres = get(config.model, "sizeinmetres", false)
    reinit = get(config.model, "reinit", true)
    do_snow = get(config.model, "snow", false)
    do_reservoirs = get(config.model, "reservoirs", false)
    do_lakes = get(config.model, "lakes", false)
    do_drains = get(config.model, "drains", false)
    do_constanthead = get(config.model, "constanthead", false)

    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
    kw_land_tstep = get(config.model, "kw_land_tstep", 0)

    nc = NCDataset(static_path)
    dims = dimnames(nc[param(config, "input.subcatchment")])

    # There is no need to permute the dimensions of the data, since the active indices are
    # correctly calculated in both ways.
    # The dimension order only needs to be known for interpreting the LDD directions
    # and creating the coordinate maps.
    dims_xy = dims[2] in ("y", "lat")

    subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
    # indices based on catchment
    inds, rev_inds = active_indices(subcatch_2d, missing)
    n = length(inds)
    modelsize_2d = size(subcatch_2d)

    river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
    river = river_2d[inds]
    riverwidth_2d =
        ncread(nc, param(config, "input.lateral.river.width"); type = Float64, fill = 0)
    riverwidth = riverwidth_2d[inds]
    riverlength_2d =
        ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
    riverlength = riverlength_2d[inds]

    # read x, y coordinates and calculate cell length [m]
    y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
    x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
    if dims_xy
        y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
    else
        y = repeat(y_nc, outer = (1, length(x_nc)))[inds]
    end
    cellength = abs(mean(diff(x_nc)))

    xl = fill(mv, n)
    yl = fill(mv, n)
    riverfrac = fill(mv, n)

    for i = 1:n
        xl[i] = sizeinmetres ? cellength : lattometres(y[i])[1] * cellength
        yl[i] = sizeinmetres ? cellength : lattometres(y[i])[2] * cellength
        riverfrac[i] =
            river[i] ? min((riverlength[i] * riverwidth[i]) / (xl[i] * yl[i]), 1.0) : 0.0
    end

    # initialize vertical SBM concept
    sbm = initialize_sbm(nc, config, riverfrac, xl, yl, inds)

    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
    nriv = length(inds_riv)

    # reservoirs
    pits = zeros(Bool, modelsize_2d)
    if do_reservoirs
        reservoirs, resindex, reservoir, pits = initialize_simple_reservoir(config, nc, inds_riv, nriv, pits)
    else
        reservoir = ()
    end

    # lakes
    if do_lakes
        lakes, lakeindex, lake, pits = initialize_natural_lake(config, nc, inds_riv, nriv, pits)
    else
        lake = ()
    end

    # overland flow (kinematic wave)
    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float64)
    clamp!(βₗ, 0.00001, Inf)
    ldd_2d = ncread(nc, param(config, "input.ldd"); allow_missing = true)
    ldd = ldd_2d[inds]
    dl = fill(mv, n)
    dw = fill(mv, n)
    sw = fill(mv, n)

    for i = 1:n
        dl[i] = detdrainlength(ldd[i], xl[i], yl[i])
        dw[i] = (xl[i] * yl[i])/dl[i]
        sw[i] = river[i] ? max(dw[i] - riverwidth[i], 0.0) : dw[i]
    end

    n_land = ncread(
        nc,
        param(config, "input.lateral.land.n", nothing);
        sel = inds,
        defaults = 0.072,
        type = Float64,
    )

    olf = SurfaceFlow(
        sl = βₗ,
        n = n_land,
        dl = dl,
        Δt = tosecond(Δt),
        its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
        width = sw,
        wb_pit = pits[inds],
    )

    pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
    graph = flowgraph(ldd, inds, pcr_dir)

    # river flow (kinematic wave)
    riverslope = ncread(
        nc,
        param(config, "input.lateral.river.slope");
        sel = inds_riv,
        type = Float64,
    )
    clamp!(riverslope, 0.00001, Inf)
    riverlength = riverlength_2d[inds_riv]
    riverwidth = riverwidth_2d[inds_riv]
    n_river = ncread(
        nc,
        param(config, "input.lateral.river.n", nothing);
        sel = inds_riv,
        defaults = 0.036,
        type = Float64,
    )
    ldd_riv = ldd_2d[inds_riv]
    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)

    # the indices of the river cells in the land(+river) cell vector
    index_river = filter(i -> !isequal(river[i], 0), 1:n)
    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)

    rf = SurfaceFlow(
        sl = riverslope,
        n = n_river,
        dl = riverlength,
        Δt = tosecond(Δt),
        its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) : kw_river_tstep,
        width = riverwidth,
        reservoir_index = do_reservoirs ? resindex : fill(0, nriv),
        lake_index = do_lakes ? lakeindex : fill(0, nriv),
        reservoir = do_reservoirs ? reservoirs : nothing,
        lake = do_lakes ? lakes : nothing,
        rivercells = river,
    )

    # unconfined aquifer
    if do_constanthead
        constanthead = ncread(
            nc,
            param(config, "input.lateral.subsurface.constant_head", nothing);
            sel = inds,
            type = Float64,
            fill = mv,
        )
        index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n)
        constant_head = ConstantHead(constanthead[index_constanthead],
                                    index_constanthead,)

    else
        constant_head = ConstantHead[]
    end

    conductivity = ncread(
        nc,
        param(config, "input.lateral.subsurface.conductivity", nothing);
        sel = inds,
        type = Float64,
    )
    specific_yield = ncread(
        nc,
        param(config, "input.lateral.subsurface.specific_yield", nothing);
        sel = inds,
        type = Float64,
    )

    connectivity = Connectivity(inds, rev_inds, sbm.xl, sbm.yl)
    initial_head = sbm.altitude .- 0.10 # cold state for groundwater head
    initial_head[index_river] = sbm.altitude[index_river]

    if do_constanthead
        initial_head[constant_head.index] = constant_head.head
    end

    aquifer = UnconfinedAquifer(
        initial_head,
        conductivity,
        sbm.altitude,
        sbm.altitude .- sbm.soilthickness ./ 1000.0,
        xl .* yl,
        specific_yield,
        fill(0.0, connectivity.nconnection),  # conductance
    )

    # reset zi and satwaterdepth with groundwater head from unconfined aquifer
    sbm.zi .= (sbm.altitude .- initial_head) .* 1000.0
    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* sbm.θₑ

    # river boundary of unconfined aquifer
    infiltration_conductance = ncread(
        nc,
        param(config, "input.lateral.subsurface.infiltration_conductance", nothing);
        sel = inds_riv,
        type = Float64,
    )
    exfiltration_conductance = ncread(
        nc,
        param(config, "input.lateral.subsurface.exfiltration_conductance", nothing);
        sel = inds_riv,
        type = Float64,
    )
    river_bottom = ncread(
        nc,
        param(config, "input.lateral.subsurface.river_bottom", nothing);
        sel = inds_riv,
        type = Float64,
    )

    river_flux = fill(mv,nriv)
    river_stage = fill(mv,nriv)
    river = River(
        river_stage,
        infiltration_conductance,
        exfiltration_conductance,
        river_bottom,
        river_flux,
        index_river,
    )

    # drain boundary of unconfined aquifer (optional)
    if do_drains
        drain = ncread(
            nc,
            param(config, "input.lateral.subsurface.drain", nothing);
            sel = inds,
            type = Bool,
        )
        drain_elevation = ncread(
            nc,
            param(config, "input.lateral.subsurface.drain_elevation", nothing);
            sel = inds,
            type = Float64,
            fill = mv,
        )
        drain_conductance = ncread(
            nc,
            param(config, "input.lateral.subsurface.drain_conductance", nothing);
            sel = inds,
            type = Float64,
            fill = mv,
        )
        index_drain = filter(i -> !isequal(drain[i], 0), 1:n)
        drain_flux = fill(mv,length(index_drain))
        drain = Drainage(
            drain_elevation[index_drain],
            drain_conductance[index_drain],
            drain_flux,
            index_drain,
        )
    else
        drain = Drainage[]
    end

    # recharge boundary of unconfined aquifer
    r = fill(mv,n)
    recharge = Recharge(r, fill(0.0, n), collect(1:n))

    gwf = GroundwaterFlow(
        aquifer,
        connectivity,
        constant_head,
        AquiferBoundaryCondition[recharge, river, drain],
    )

    state_ncnames = ncnames(config.state)

    reader = prepare_reader(dynamic_path, cyclic_path, config)

    modelmap = (vertical = sbm,
                lateral = (subsurface = (flow = gwf, recharge = gwf.boundaries[1], river = gwf.boundaries[2], drain = gwf.boundaries[3]),
                land = olf,
                river = rf))
    indices_reverse = (
        land = rev_inds,
        river = rev_inds_riv,
        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
        lake = isempty(lake) ? nothing : lake.reverse_indices,
    )
    writer = prepare_writer(
        config,
        reader,
        output_path,
        modelmap,
        state_ncnames,
        indices_reverse,
        x_nc,
        y_nc,
        dims_xy,
        nc,
        maxlayers = sbm.maxlayers,
    )

    # for each domain save the directed acyclic graph, the traversion order,
    # and the indices that map it back to the two dimensional grid
    # for reservoirs and lakes this is information is also available (except the graph)
    # from the initialization functions
    land = (
        graph = graph,
        order = topological_sort_by_dfs(graph),
        indices = inds,
        reverse_indices = rev_inds,
    )
    river = (
        graph = graph_riv,
        order = topological_sort_by_dfs(graph_riv),
        indices = inds_riv,
        reverse_indices = rev_inds_riv,
    )

    model = Model(
        config,
        (; land, river, reservoir, lake, index_river, frac_toriver),
        (subsurface = (flow = gwf, recharge = gwf.boundaries[1], river = gwf.boundaries[2], drain = gwf.boundaries[3]),
        land = olf,
        river = rf),
        sbm,
        Clock(config.starttime, 1, Δt),
        reader,
        writer,
    )

    # read and set states in model object if reinit=false
    if reinit  == false
        set_states(
            instate_path,
            model,
            state_ncnames,
            type = Float64
        )
    end

    # make sure the forcing is already loaded
    # it's fine to run twice, and may help catching errors earlier
    update_forcing!(model)
    if haskey(config.input, "cyclic")
        update_cyclic!(model)
    end

    return model
end

"update the sbm_gwf model for a single timestep"
function update_sbm_gwf(model)
    @unpack lateral, vertical, network, clock, config = model

    inds_riv = network.index_river
    kinwave_it = get(config.model, "kin_wave_iteration", false)

    update_forcing!(model)
    if haskey(config.input, "cyclic")
        update_cyclic!(model)
    end

    # extract water levels h_av [m] from the land and river domains
    # this is used to limit open water evaporation
    vertical.waterlevel_land .= lateral.land.h_av .* 1000.0
    vertical.waterlevel_river[inds_riv] .= lateral.river.h_av .* 1000.0

    # vertical sbm concept is updated until snow state, after that (optional)
    # snow transport is possible
    update_until_snow(vertical, config)

    # lateral snow transport
    if get(config.model, "masswasting", false)
        snowflux_frac =
            min.(0.5, lateral.land.sl ./ 5.67) .* min.(1.0, vertical.snow ./ 10000.0)
        maxflux = snowflux_frac .* vertical.snow
        vertical.snow .= accucapacityflux(network.land, vertical.snow, maxflux)
        vertical.snowwater .= accucapacityflux(network.land, vertical.snowwater, vertical.snowwater .* snowflux_frac)
    end

    # update vertical sbm concept until recharge [mm]
    update_until_recharge(vertical, config)

    # set river stage (groundwater) to average h from kinematic wave
    lateral.subsurface.river.stage .= lateral.river.h_av .+ lateral.subsurface.river.bottom

    # determine stable time step for groundwater flow
    Δt_gw = stable_timestep(lateral.subsurface.flow.aquifer) # time step in day (Float64)
    Δt_sbm = (vertical.Δt / 86400.0) # vertical.Δt is in seconds (Float64)
    if Δt_gw < Δt_sbm
        @warn("stable time step Δt $Δt for groundwater flow is smaller than sbm Δt $Δt_sbm")
    end

    Q = zeros(vertical.n)
    # exchange of recharge between vertical sbm concept and groundwater flow domain
    lateral.subsurface.recharge.rate .= vertical.recharge ./ 1000.0
    # update groundwater domain
    update(lateral.subsurface.flow, Q, Δt_sbm)

    # determine excess head [m] (exfiltwater) in groundwater domain (head > surface) and reset head
    exfiltwater = lateral.subsurface.flow.aquifer.head .- min.(lateral.subsurface.flow.aquifer.head, lateral.subsurface.flow.aquifer.top)
    lateral.subsurface.flow.aquifer.head .= min.(lateral.subsurface.flow.aquifer.head, lateral.subsurface.flow.aquifer.top)

    # update vertical sbm concept (runoff, ustorelayerdepth and satwaterdepth)
    update_after_lateralflow(
        vertical,
        (vertical.altitude .- lateral.subsurface.flow.aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
        exfiltwater .* 1000.0,
    )

    # update overland flow based on vertical runoff [mm] from vertical sbm concept
    lateral.land.qlat .=
        (vertical.runoff .* vertical.xl .* vertical.yl .* 0.001) ./ lateral.land.Δt ./
        lateral.land.dl
    update(
        lateral.land,
        network.land,
        frac_toriver = network.frac_toriver,
        river = lateral.river.rivercells,
        do_iter = kinwave_it,
    )

    # update river domain with net runoff from vertical sbm concept, overland flow
    # and groundwater flow to the river cells
    net_runoff_river =
        (vertical.net_runoff_river[inds_riv] .* vertical.xl[inds_riv] .*
        vertical.yl[inds_riv] .* 0.001) ./ vertical.Δt

    # flux from groundwater domain to river (Q to river from drains and groundwater)
    flux_gw = zeros(vertical.n)
    flux_gw[lateral.subsurface.river.index] = -lateral.subsurface.river.flux
    flux_gw[lateral.subsurface.drain.index] = flux_gw[lateral.subsurface.drain.index] - lateral.subsurface.drain.flux

    lateral.river.qlat .=
        (
            flux_gw[inds_riv] ./ lateral.river.Δt .+
            lateral.land.to_river[inds_riv] .+ net_runoff_river
        ) ./ lateral.river.dl
    update(lateral.river, network.river, do_iter = kinwave_it, doy = dayofyear(clock.time))

    write_output(model)

    # update the clock
    advance!(clock)

    return model
end


text:
"""
    simulate(N, models::OrderedDict{Symbol,LinearMixedModel}, saveresults::Function, β, σ, θ)
`N` replications of data are simulated from model `sim` in the ordered dictionary
`models` using the parameter values `β`, `σ` and `θ`.  All the models are fit to
this response, after which `saveresults` is called with arguments `models` and `i`,
the current index in the range `1:N`.   A simple `saveresults` callback would insert
the current value of the objective (i.e. -2*loglikelihood(m)) in a
`length(models) × N` array.
"""
function simulate(N, models::OrderedDict{Symbol,LinearMixedModel}, sim::Symbol,
                  saveresults::Function, β, σ, θ)
    simmod = models[sim]
    for i in 1 : N
        y = simulate!(simmod, β = β, σ = σ, θ = θ)
        for m in values(models)
            refit!(m, y)
        end
        saveresults(models, i)
    end
end

const dat = makedata(50, 20, DataFrame(C = [-0.5, 0.5]));

const mods = OrderedDict{Symbol, LinearMixedModel}(
    :max1 => lmm(y ~ 1 + C + (1 + C | S) + (1 + C | I), dat),
    :max0 => lmm(y ~ 1 + (1 + C | S) + (1 + C | I), dat),
    :nrc1 => lmm(y ~ 1 + C + (1 | S) + (0 + C | S) + (1 | I) + (0 + C | I), dat),
    :nrc0 => lmm(y ~ 1 + (1 | S) + (0 + C | S) + (1 | I) + (0 + C | I), dat),
    :zis1 => lmm(y ~ 1 + C + (1 | S) + (0 + C | S) + (1 | I), dat),
    :zis0 => lmm(y ~ 1 + (1 | S) + (0 + C | S) + (1 | I), dat),
#    :zss1 => lmm(y ~ 1 + C + (1 | S) + (1 | I) + (0 + C | I), dat),  # currently fails in refit
#    :zss0 => lmm(y ~ 1 + (1 | S) + (1 | I) + (0 + C | I), dat),  # currently fails in refit
    :rin1 => lmm(y ~ 1 + C + (1 | S) + (1 | I), dat),
    :rin0 => lmm(y ~ 1 + (1 | S) + (1 | I), dat)
)


text:
module PIAssistant
    # Developed by:
    # PhD. Arick Castillo Landero
    # PhD. Arturo Jimenez Gutierrez
    # PhD. Kelvyn B. Sánchez Sánchez
    # Instituto Tecnológico de Celaya/TecNM - México.

    # Function to call the GUI
    export mainPI
    export PIAssistantGUI

    # Include the main file .fl
    include("index.jl")
    include("main.jl")
end


text:
module zypher
  include("display.jl")
  export display_cpa, display_dpa, display_ta, cpa, buildTemplates, ta, get_pois
  export SPNdata, loadSPNtrace, loadTA_TRAINtrace, loadTA_TESTtrace
end

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
module Scotch

# Define before execution, for instance in your /~.juliarc.jl add the following line
# const SCOTCH_LIBSCOTCH = "/path/to/scotch/lib/libscotch.dylib"
const libscotch = Main.SCOTCH_LIBSCOTCH

mutable struct Graph
end

mutable struct Strat
end

mutable struct Ordering
end

mutable struct File # libc's File
end

function version()
    ver = Cint[0]
    rel = Cint[0]
    pat = Cint[0]
    ccall((:SCOTCH_version, libscotch), Void,
            (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
            pointer(ver), pointer(rel), pointer(pat))
    return (ver[1], rel[1], pat[1])
end

function graphAlloc()
    graph = ccall((:SCOTCH_graphAlloc, libscotch), Ptr{Graph}, ())
    if graph == C_NULL
        error("Error in graphAlloc")
    end
    return graph
end

function graphInit(graph::Ptr{Graph})
    err = ccall((:SCOTCH_graphInit, libscotch), Cint,
                (Ptr{Graph},),
                graph)
    err == 0 || error("Error in graphInit")
    return
end

function graphLoad(graph::Ptr{Graph}, filename, baseval::Integer, flagval::Integer)
    m = "r"
    file = ccall((:fopen, "libc"), Ptr{File},
             (Ptr{Char}, Ptr{Char}),
             pointer(filename), pointer(m))
    file != C_NULL || error("Error in graphInit, fopen")
    err = ccall((:SCOTCH_graphLoad, libscotch), Cint,
                (Ptr{Graph}, Ptr{File}, Cint, Cint),
                graph, file, Cint(baseval), Cint(flagval))
    err == 0 || error("Error in graphInit")
    return
end

function graphBuild(graph::Ptr{Graph}, A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti <: Integer}
    A.n == A.m || error("Matrix should be square")
    # Remove self loops
    vertnbr = A.n
    verttab = Ptr{Cint}(Libc.malloc(sizeof(Cint) * (vertnbr+1)))
    edgetab = Ptr{Cint}(Libc.malloc(sizeof(Cint) * length(A.rowval)))
    unsafe_store!(verttab, 1)
    k = 1
    edgenbr = 0
    for i = 1:vertnbr
        for j = A.colptr[i]:(A.colptr[i+1]-1)
            if A.rowval[j] != i
                unsafe_store!(edgetab, A.rowval[j], k)
                k += 1
                edgenbr += 1
            end
        end
        unsafe_store!(verttab, k, i+1)
    end
    # FIXME: resize edgetab to k
    err = ccall((:SCOTCH_graphBuild, libscotch), Cint,
          (Ptr{Graph}, Cint, Cint,    Ptr{Cint},  Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint,    Ptr{Cint},  Ptr{Cint}),
          graph,       1,    vertnbr, verttab,    C_NULL,    C_NULL,    C_NULL,    edgenbr, edgetab,    C_NULL)
    err == 0 || error("Error in graphBuild")
    return
end

function graphCheck(graph::Ptr{Graph})
    err = ccall((:SCOTCH_graphCheck, libscotch), Cint,
                (Ptr{Graph},),
                graph)
    err == 0 || error("Error in graphCheck")
end

function graphSize(graph::Ptr{Graph})
    nv = Cint[0]
    ne = Cint[0]
    ccall((:SCOTCH_graphSize, libscotch), Void,
               (Ptr{Graph}, Ptr{Cint}, Ptr{Cint}),
               graph, pointer(nv), pointer(ne))
    return (nv[1], ne[1])
end

function graphFree(graph::Ptr{Graph})
    ccall((:SCOTCH_graphFree, libscotch), Void,
          (Ptr{Graph},),
          graph)
end

function graphExit(graph::Ptr{Graph})
    ccall((:SCOTCH_graphExit, libscotch), Void,
          (Ptr{Graph},),
          graph)
end

function graphStat(graph::Ptr{Graph})
    A1 = Cint[0,0,0,0,0,0,0,0]
    A2 = Cdouble[0,0,0,0,0,0]
    @show A1, A2
    ccall((:SCOTCH_graphStat, libscotch), Void,
          (Ptr{Graph},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cint}, Ptr{Cint},            Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}),
          graph,
          pointer(A1, 1), pointer(A1, 2), pointer(A1, 3), pointer(A2, 1), pointer(A2, 2),
          pointer(A1, 4), pointer(A1, 5),                 pointer(A2, 3), pointer(A2, 4),
          pointer(A1, 6), pointer(A1, 7), pointer(A1, 8), pointer(A2, 5), pointer(A2, 6))
    @show A1, A2
    return (A1[1], A1[2], A1[3], A2[1], A2[2],
            A1[4], A1[5],        A2[3], A2[4],
            A1[6], A1[7], A1[8], A2[5], A2[6])
end

# Strategy

function stratAlloc()
    strat = ccall((:SCOTCH_stratAlloc, libscotch), Ptr{Strat}, ())
    strat != C_NULL || error("Error in stratAlloc")
    return strat
end

function stratInit(strat::Ptr{Strat})
    err = ccall((:SCOTCH_stratInit, libscotch), Cint,
                (Ptr{Strat},),
                strat)
    err == 0 || error("Error in stratInit")
    return
end

function stratGraphOrder(strat::Ptr{Strat}, st::String)
    err = ccall((:SCOTCH_stratGraphOrder, libscotch), Cint,
                (Ptr{Strat}, Ptr{Char}),
                strat, pointer(st))
    err == 0 || error("Error in stratGraphOrder")
    return
end

function stratExit(strat::Ptr{Strat})
    ccall((:SCOTCH_stratExit, libscotch), Void,
          (Ptr{Strat},),
          strat)
end

# Ordering

function graphOrder(graph::Ptr{Graph}, strat::Ptr{Strat})
    (nv, ne) = graphSize(graph)
    permtab = Array{Cint, 1}(nv)
    peritab = Array{Cint, 1}(nv)
    cblkptr = Array{Cint, 1}(1)
    rangtab = Array{Cint, 1}(nv + 1)
    treetab = Array{Cint, 1}(nv)
    err = ccall((:SCOTCH_graphOrder, libscotch), Cint,
            (Ptr{Graph}, Ptr{Strat}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
            graph,       strat,      pointer(permtab), pointer(peritab), pointer(cblkptr), pointer(rangtab), pointer(treetab))
    err == 0 || error("Error in graphOrder")
    return (permtab, peritab, cblkptr, rangtab, treetab)
end

function orderAlloc()
    optr = ccall((:SCOTCH_orderAlloc, libscotch), Ptr{Ordering}, ())
    optr != C_NULL || error("Error in orderAlloc")
    return optr
end

function graphOrderInit(graph::Ptr{Graph}, ordering::Ptr{Ordering}, permtab::Array{Cint,1}, peritab::Array{Cint,1}, cblkptr::Array{Cint,1}, rangtab::Array{Cint,1}, treetab::Array{Cint,1})
    err = ccall((:SCOTCH_graphOrderInit, libscotch), Cint,
                (Ptr{Graph}, Ptr{Ordering}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                graph, ordering, pointer(permtab), pointer(peritab), pointer(cblkptr), pointer(rangtab), pointer(treetab))
    err == 0 || error("Error in graphOrderInit")
    return
end

function graphOrderComputeList(graph::Ptr{Graph}, ordering::Ptr{Ordering}, listtab::Array{Cint,1}, strat::Ptr{Strat})
    err = ccall((:SCOTCH_graphOrderComputeList, libscotch), Cint,
                (Ptr{Graph}, Ptr{Ordering}, Cint, Ptr{Cint}, Ptr{Strat}),
                graph, ordering, Cint(length(listtab)), pointer(listtab), strat)
    err == 0 || error("Error in graphOrderComputeList")
    return
end

function graphOrderCompute(graph::Ptr{Graph}, ordering::Ptr{Ordering}, strat::Ptr{Strat})
    err = ccall((:SCOTCH_graphOrderCompute, libscotch), Cint,
                (Ptr{Graph}, Ptr{Ordering}, Ptr{Strat}),
                graph, ordering, strat)
    err == 0 || error("Error in graphOrderCompute")
    return
end

function graphOrderCheck(graph::Ptr{Graph}, ordering::Ptr{Ordering})
    err = ccall((:SCOTCH_graphOrderCheck, libscotch), Cint,
                (Ptr{Graph}, Ptr{Ordering}),
                graph, ordering)
    err == 0 || error("Error in graphOrderCheck")
end

function graphOrderExit(graph::Ptr{Graph}, ordering::Ptr{Ordering})
    ccall((:SCOTCH_graphOrderExit, libscotch), Void,
                (Ptr{Graph}, Ptr{Ordering}),
                graph, ordering)
end

end # module Scotch


text:
include("setup.jl")

@testset "JET.jl" begin
    @testset "abstractinterpret" begin
        @testset "inferenceerrorreport.jl" begin
            include("abstractinterpret/test_inferenceerrorreport.jl")
        end

        @testset "typeinfer.jl" begin
            include("abstractinterpret/test_typeinfer.jl")
        end
    end

    @testset "toplevel" begin
        @testset "virtualprocess.jl" begin
            include("toplevel/test_virtualprocess.jl")
        end
    end

    @testset "ui" begin
        # tests with Windows-paths is just an hell
        @static Sys.iswindows() || @testset "print.jl" begin
            include("ui/test_print.jl")
        end

        @testset "vscode.jl" begin
            include("ui/test_vscode.jl")
        end
    end

    @testset "misc" begin
        include("test_misc.jl")
    end

    @testset "Test.jl integration" begin
        include("test_Test.jl")
    end

    @testset "JETInterface" begin
        include("test_JETInterface.jl")
    end

    @testset "analyzers" begin
        @testset "JETAnalyzer" begin
            include("analyzers/test_jetanalyzer.jl")
        end

        @testset "OptAnalyzer" begin
            include("analyzers/test_optanalyzer.jl")
        end
    end

    # run simple performance benchmark
    # and also run OptAnalyzer on JET
    @testset "performance" begin
        include("performance.jl")
    end

    @testset "self check !!!" begin
        target_modules = (JET,)
        JETAnalyzerT   = typeof(JET.JETAnalyzer())
        OptAnalyzerT   = typeof(JET.OptAnalyzer())
        InferenceState = Core.Compiler.InferenceState

        # error analysis
        # ==============

        # JETAnalyzer
        test_call(JET.analyze_frame!, (JETAnalyzerT, InferenceState);
                  target_modules)
        # OptAnalyzer
        test_call(JET.analyze_frame!, (OptAnalyzerT, InferenceState);
                  target_modules)
        # top-level
        test_call(JET.virtual_process, (String, String, JETAnalyzerT, JET.ToplevelConfig);
                  target_modules)
        # entries
        test_call(JET.report_file, (String,);
                  target_modules)
        test_call(JET.report_package, (Union{String,Module,Nothing},);
                  target_modules)

        # optimization analysis
        # =====================

        function function_filter(@nospecialize ft)
            if ft === typeof(JET.widenconst) ||
               ft === typeof(JET.print) ||
               ft === typeof(Base.CoreLogging.handle_message) ||
               ft == Type{<:JET.InferenceErrorReport} # the constructor used in `restore_cached_report` is very dynamic
                return false
            end
            return true
        end
        # JETAnalyzer
        test_opt(JET.analyze_frame!, (JETAnalyzerT, InferenceState);
                 target_modules,
                 function_filter)
        # OptAnalyzer
        test_opt(JET.analyze_frame!, (OptAnalyzerT, InferenceState);
                 target_modules,
                 function_filter)
    end
end


text:
using LinearAlgebra: dot, norm, diagm
# using SymPy: symbols

@testset "Test length in a hexagonal lattice" begin
    g = MetricTensor(1, 1, 2, 90, 90, 120)  # Primitive hexagonal
    @test g ≈ MetricTensor([1, 0, 0], [-1 // 2, sqrt(3) / 2, 0], [0, 0, 2])
    @test g ≈ MetricTensor(Lattice([
        1.0 -1/2 0.0
        0.0 sqrt(3)/2 0.0
        0.0 0.0 2.0
    ]))
    a = [1, 2, 1]
    @test dot(a, g, a) ≈ 7
    @test norm([1, 2, 1], g)^2 ≈ 7
end

@testset "Test distance between atoms in a hexagonal lattice" begin
    g = MetricTensor(1, 1, 2, 90, 90, 120)  # Primitive hexagonal
    a = [1, 1, 1]
    b = [1 // 3, 1 // 3, 1 // 2]
    @test distance(a, g, b)^2 == 13 / 9
end

@testset "Test direction cosine in a tetragonal lattice" begin
    g = MetricTensor(2, 2, 3, 90, 90, 90)  # Primitive tetragonal
    @test g == MetricTensor(Lattice(diagm([2, 2, 3])))
    @test g == MetricTensor([2, 0, 0], [0, 2, 0], [0, 0, 3])
    a = [1, 2, 1]
    b = [0, 0, 1]
    @test directioncosine(a, g, b)^2 ≈ 9 / 29
end

# @testset "Symbolic calculation" begin
#     a, c = symbols("a, c", positive = true)
#     @test MetricTensor(a, a, c, 90, 90, 120) ==
#           MetricTensor([a^2 -0.5*a^2 0; -0.5*a^2 a^2 0; 0 0 c^2])  # Primitive hexagonal
# end

# From https://ssd.phys.strath.ac.uk/wp-content/uploads/Crystallographic_maths.pdf
@testset "Find angle θ between [100] and [111] directions in GaN" begin
    a, c = 3.19u"angstrom", 5.19u"angstrom"
    g = MetricTensor(a, a, c, 90, 90, 120)
    @test directioncosine([1, 0, 0], g, [1, 1, 1]) == a / 2 / sqrt(a^2 + c^2)
end


text:
include("./lib/metadata.jl")
include("./lib/pccf.jl")
include("./1naive_update.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")
#include("./sim1.jl")
#include("./exp1.jl")

using BayesChpOrig
using Metadata
using Roc
#using Sim1
#using Exp1
using SignalGenerator

# test oise adder
sig = rand(100)
addNoiseRec!(sig, [10, 40, 70, 90], 50, width = 2, maxlen=10)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
module Unwrap

include("interface.jl")
include("unwrap_1d.jl")
include("unwrap_common.jl")
include("unwrap_2d.jl")
include("unwrap_3d.jl")

export unwrap, unwrap!

end # module

max_forks_repo_name: platawiec/Unwrap


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures for the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export TextFormat, Highlighter

"""
    TextFormat

# Fields

* `up_right_corner`: Character in the up right corner.
* `up_left_corner`: Character in the up left corner.
* `bottom_left_corner`: Character in the bottom left corner.
* `bottom_right_corner`: Character in the bottom right corner.
* `up_intersection`: Character in the intersection of lines in the up part.
* `left_intersection`: Character in the intersection of lines in the left part.
* `right_intersection`: Character in the intersection of lines in the right
                        part.
* `middle_intersection`: Character in the intersection of lines in the middle of
                         the table.
* `bottom_intersection`: Character in the intersection of the lines in the
                         bottom part.
* `column`: Character in a vertical line inside the table.
* `left_border`: Character used as the left border.
* `right_border`: Character used as the right border.
* `row`: Character in a horizontal line inside the table.
* `hlines`: Horizontal lines that must be drawn by default.
* `vlines`: Vertical lines that must be drawn by default.

# Pre-defined formats

The following pre-defined formats are available: `unicode` (**default**),
`mysql`, `compact`, `markdown`, `simple`, `ascii_rounded`, and `ascii_dots`.

"""
@with_kw struct TextFormat
    up_right_corner::Char                 = '┐'
    up_left_corner::Char                  = '┌'
    bottom_left_corner::Char              = '└'
    bottom_right_corner::Char             = '┘'
    up_intersection::Char                 = '┬'
    left_intersection::Char               = '├'
    right_intersection::Char              = '┤'
    middle_intersection::Char             = '┼'
    bottom_intersection::Char             = '┴'
    column::Char                          = '│'
    row::Char                             = '─'
    hlines::Vector{Symbol}                = [:begin,:header,:end]
    vlines::Union{Symbol,Vector{Symbol}}  = :all
end

"""
    Highlighter

Defines the default highlighter of a table when using the text backend.

# Fields

* `f`: Function with the signature `f(data,i,j)` in which should return `true`
       if the element `(i,j)` in `data` must be highlighter, or `false`
       otherwise.
* `fd`: Function with the signature `f(h,data,i,j)` in which `h` is the
        highlighter. This function must return the `Crayon` to be applied to the
        cell that must be highlighted.
* `crayon`: The `Crayon` to be applied to the highlighted cell if the default
            `fd` is used.

# Remarks

This structure can be constructed using three helpers:

    Highlighter(f::Function; kwargs...)

where it will construct a `Crayon` using the keywords in `kwargs` and apply it
to the highlighted cell,

    Highlighter(f::Function, crayon::Crayon)

where it will apply the `crayon` to the highlighted cell, and

    Highlighter(f::Function, fd::Function)

where it will apply the `Crayon` returned by the function `fd` to the
highlighted cell.

"""
@with_kw struct Highlighter
    f::Function
    fd::Function = (h,data,i,j)->h.crayon

    # Private
    crayon::Crayon = Crayon()
end

Highlighter(f; kwargs...) = Highlighter(f = f, crayon = Crayon(;kwargs...))
Highlighter(f, crayon::Crayon) = Highlighter(f = f, crayon = crayon)
Highlighter(f::Function, fd::Function) = Highlighter(f = f, fd = fd)

"""
    Screen

Store the information of the screen and the current cursor position. Notice that
this is not the real cursor position with respect to the screen, but with
respect to the point in which the table is printed.

# Fields

* `size`: Screen size.
* `row`: Current row.
* `col`: Current column.
* `has_color`: Indicates if the screen has color support.

"""
@with_kw mutable struct Screen
    size::Tuple{Int,Int} = (-1,-1)
    row::Int             = 1
    col::Int             = 0
    has_color::Bool      = false
end


text:
module Colab

export ipynb_to_colab
include("functions.jl")

end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

function BoxCalculation(Vertices)
	Minx=minimum(Vertices[1,:])
	Maxx=maximum(Vertices[1,:])
	Miny=minimum(Vertices[2,:])
	Maxy=maximum(Vertices[2,:])
	dx=Maxx-Minx
	dy=Maxy-Miny
	Box=dx*dy
	if size(Vertices,1)==3
		Minz=minimum(Vertices[3,:])
		Maxz=maximum(Vertices[3,:])
		dz=Maxz-Minz
		Box=Box*dz
	end
	return Box
end

@testset "circle" begin
	@test BoxCalculation(Lar.circle()()[1])==4
	@test BoxCalculation(Lar.circle(2., 2*pi)()[1])==16
	@test BoxCalculation(Lar.circle(3, 2*pi)()[1])==36
	@test BoxCalculation(Lar.circle(5, pi/2)()[1])==25
	#@test size(Lar.circle(3,2*pi)(60)[1],2)==60
	@test length(Lar.circle(3,2*pi)([60])[2])==60
end

@testset "helix" begin
	@test BoxCalculation(Lar.helix()()[1])==8
	@test BoxCalculation(Lar.helix(2, 2, 2)()[1])==64
	@test BoxCalculation(Lar.helix(1, 2, 5)()[1])==40
	@test BoxCalculation(Lar.helix(1, 1, 3)()[1])==12
	#=
	1, 1, 2, is contained
	in a box that has volume (radius*2)^2*pitch*nturns
	=#
	@test size(Lar.helix(5, 7, 9)()[1],2)==325
	@test length(Lar.helix(5, 7, 9)()[2])==324
end

@testset "disk" begin
	@test BoxCalculation(Lar.disk(1., 2*pi)([36, 1])[1])==4
	@test BoxCalculation(Lar.disk(2., 2*pi)()[1])==16
	@test BoxCalculation(Lar.disk(2., pi)()[1])==8
	@test BoxCalculation(Lar.disk(2., pi/2)()[1])==4
	@test size(Lar.disk(10, pi/7)()[1],2)==75
	@test length(Lar.disk(10, pi/7)()[2])==108
end

@testset "helicoid" begin
	@test BoxCalculation(Lar.helicoid()()[1])==8
	@test BoxCalculation(Lar.helicoid(2,1,2,2)()[1])==64
	@test BoxCalculation(Lar.helicoid(1,.5,2,5)()[1])==40
	@test BoxCalculation(Lar.helicoid(1,.3,1,3)()[1])==12
	#=
	1, 1, 2, is contained in a
	box that has volume (radius*2)^2*pitch*nturns
	=#
	@test size(Lar.helicoid(1.3,.7,1,3)()[1],2)==327
	@test length(Lar.helicoid(1.3,.7,1,3)()[2])==432
end

@testset "ring" begin
	@test BoxCalculation(Lar.ring(1.,3.,2*pi)()[1])==36
	@test BoxCalculation(Lar.ring(1.,2,pi)()[1])==8
	@test BoxCalculation(Lar.ring(1.,5,pi/2)()[1])==25
	#(radius*2)^2
	@test size(Lar.ring(1,5,pi/2)()[1],2)==74
	@test length(Lar.ring(5,10,pi/6)()[2])==72
end

@testset "cylinder" begin
	@test BoxCalculation(Lar.cylinder(1,5,2*pi)()[1])==20
	@test BoxCalculation(Lar.cylinder(2,2,pi)()[1])==16
	@test BoxCalculation(Lar.cylinder(1,4,pi)()[1])==8
	#((radius*2)^2)*height
	@test size(Lar.cylinder(3.4,20,pi/7)()[1],2)==74
	@test length(Lar.cylinder(3.4,20,pi/7)()[2])==72
end

@testset "sphere" begin
	@test BoxCalculation(Lar.sphere(2,pi,2*pi)()[1])==64
	@test BoxCalculation(Lar.sphere(6,pi,pi)()[1])==864
	@test BoxCalculation(Lar.sphere(4,pi,2*pi)()[1])==8^3
	@test size(Lar.sphere(2.5,pi/3,pi/5)()[1],2)==703
	@test length(Lar.sphere(2.5,pi/3,pi/5)()[2])==1296
end

@testset "toroidal" begin
	@test BoxCalculation(Lar.toroidal(1,3,2*pi,2*pi)()[1])==128
	@test BoxCalculation(Lar.toroidal(2,3,2*pi,2*pi)()[1])==400
	#(((R+r)*2)^2)*(r*2)
	@test size(Lar.toroidal(1.3,4.6,pi/4,pi/7)()[1],2)==925
	@test length(Lar.toroidal(1.3,4.6,pi/4,pi/7)()[2])==1728
end

@testset "crown" begin
	@test BoxCalculation(Lar.crown(1., 3., 2*pi)()[1])==128
	@test BoxCalculation(Lar.crown(2., 3., 2*pi)()[1])==400
	#(((R+r)*2)^2)*(r*2)
	@test size(Lar.crown(1.5, 5.6, pi/8)()[1],2)==481
	@test length(Lar.crown(1.5, 5.6, pi/8)()[2])==864
end

@testset "ball" begin
	@test BoxCalculation(Lar.ball(1, pi, 2*pi)()[1])==8
	@test BoxCalculation(Lar.ball(6, pi, pi)()[1])==864
	@test BoxCalculation(Lar.ball(3, pi, 2*pi)()[1])==6^3
	@test size(Lar.ball(2.6, pi/5, pi/9)()[1],2)==2813
	@test length(Lar.ball(2.6, pi/5, pi/9)()[2])==2592
end

@testset "rod" begin
	@test BoxCalculation(Lar.rod(1,5,2*pi)()[1])==20
	@test BoxCalculation(Lar.rod(2,2,pi)()[1])==16
	@test BoxCalculation(Lar.rod(1,4,pi)()[1])==8
	#((radius*2)^2)*height
	@test size(Lar.rod(3.7,8.9,pi/9)()[1],2)==74
	@test length(Lar.rod(3.7,8.9,pi/9)()[2])==1
end

@testset "hollowCyl" begin
	@test BoxCalculation(Lar.hollowCyl(0,1.,5,2*pi)()[1])==20
	@test BoxCalculation(Lar.hollowCyl(1,2.,4,pi)()[1])==32
	@test BoxCalculation(Lar.hollowCyl(1,4,3,pi/2)()[1])==48
	#((radius*2)^2)*height
	@test size(Lar.hollowCyl(3,4.,7.8,pi/5)()[1],2)==148
	@test length(Lar.hollowCyl(3,4.,7.8,pi/5)()[2])==36
end

@testset "hollowBall" begin
	@test BoxCalculation(Lar.hollowBall(1,2,pi,2*pi)([36,36,1])[1]) ==64
	@test BoxCalculation(Lar.hollowBall(1,6,pi,pi)([36,36,1])[1])== 864
	@test BoxCalculation(Lar.hollowBall(2,4,pi,2*pi)([36,36,1])[1])==8^3
	#(radius*2)^3
	@test size(Lar.hollowBall(1.5,6.7,pi/3,2*pi/3)()[1],2)==3700
	@test length(Lar.hollowBall(1.5,6.7,pi/3,2*pi/3)()[2])==2592
end

@testset "torus" begin
	@test BoxCalculation(Lar.torus(1.,2.,.5,2*pi,2*pi)()[1])==147
	@test BoxCalculation(Lar.torus(2,3,.5,2*pi,2*pi)()[1])==605.0
	#(((R+r)*2)^2)*(r*2)
	@test size(Lar.torus(5.2,7,.5,pi/3,pi/4)()[1],2)==4625
	@test length(Lar.torus(5.2,7,.5,pi/3,pi/4)()[2])==3456
end


text:
"""
Returns the a total loss after trying a 'Direct multiple shooting' on ODE data, predictions on the whole ODE data and an array of predictions from the each of the groups (smaller intervals).
In Direct Multiple Shooting, the Neural Network divides the interval into smaller intervals and solves for them separately.
The default group size is 5 implying the whole dataset would be divided in groups of 5 and the Neural Network will solve on them individually.
The default continuity term is 100 implying any losses arising from the non-continuity of 2 different groups will be scaled by 100.

```julia
multiple_shoot(p,ode_data,tsteps,prob,loss_function,grp_size,continuity_strength=100)
```
Arguments:
- `p`: The parameters of the Neural Network to be trained.
- `ode_data`: Original Data to be modelled.
- `tsteps`: Timesteps on which ode_data was calculated.
- `prob`: ODE problem that the Neural Network attempts to solve.
- `loss_function`: Any arbitrary function to calculate loss.
- `grp_size`: The group size achieved after splitting the ode_data into equal sizes.
- `continuity_term`: Multiplying factor to ensure continuity of predictions throughout different groups.

!!!note
The parameter 'continuity_term' should be a relatively big number to enforce a large penalty whenever the last point of any group doesn't coincide with the first point of next group.
"""
function multiple_shoot(
    p::AbstractArray,
    ode_data::AbstractArray,
    tsteps,
    prob::ODEProblem,
    loss_function::Function,
    solver::DiffEqBase.AbstractODEAlgorithm,
    grp_size::Integer,
    continuity_term::Real=100
)
	datasize = length(ode_data[1,:])

	if grp_size < 1 || grp_size > datasize
        throw(DomainError(grp_size, "grp_size can't be < 1 or > number of data points"))
    end

	tot_loss = 0

	if(grp_size == datasize)
		grp_predictions = [solve(remake(prob, p=p, tspan=(tsteps[1],tsteps[datasize]), u0=ode_data[:,1]), solver, saveat=tsteps)]
		tot_loss = loss_function(ode_data, grp_predictions[1][:,1:grp_size])
		return tot_loss, grp_predictions[1], grp_predictions
	end

	if(grp_size == 1)
		# store individual single shooting predictions for each group
		grp_predictions = [solve(remake(prob, p=p, tspan=(tsteps[i],tsteps[i+1]), u0=ode_data[:,i]), solver, saveat=tsteps[i:i+1]) for i in 1:datasize-1]

		# calculate multiple shooting loss from the single shoots done in above step
		for i in 1:datasize-1
		tot_loss += loss_function(ode_data[:,i:i+1], grp_predictions[i][:, 1:grp_size]) + (continuity_term * sum(abs,grp_predictions[i][:,2] - ode_data[:,i+1]))
		end

		# single shooting predictions from ode_data[:,1] (= u0)
		pred = solve(remake(prob, p=p, tspan=(tsteps[1],tsteps[datasize]), u0=ode_data[:,1]), solver, saveat=tsteps)
		return tot_loss, pred, grp_predictions
	end

	# multiple shooting predictions
	grp_predictions = [solve(remake(prob, p=p, tspan=(tsteps[i],tsteps[i+grp_size-1]), u0=ode_data[:,i]), solver, saveat=tsteps[i:i+grp_size-1]) for i in 1:grp_size-1:datasize-grp_size]

	# calculate multiple shooting loss
	for i in 1:grp_size-1:datasize-grp_size
		# The term `trunc(Integer,(i-1)/(grp_size-1)+1)` goes from 1, 2, ... , N where N is the total number of groups that can be formed from `ode_data` (In other words, N = trunc(Integer, (datasize-1)/(grp_size-1)))
		tot_loss += loss_function(ode_data[:,i:i+grp_size-1], grp_predictions[trunc(Integer,(i-1)/(grp_size-1)+1)][:, 1:grp_size]) + (continuity_term * sum(abs,grp_predictions[trunc(Integer,(i-1)/(grp_size-1)+1)][:,grp_size] - ode_data[:,i+grp_size-1]))
	end

	# single shooting prediction
	pred = solve(remake(prob, p=p, tspan=(tsteps[1],tsteps[datasize]), u0=ode_data[:,1]), solver, saveat=tsteps)
	return tot_loss, pred, grp_predictions
end


text:
# module for logging robot's trajectory and observation
# logged data is used for graph based slam simulation

module TrajObsrvLogger
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function animate_per_time(time_interval, world, map, robot)
    draw(world)
    annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
    draw!(map)
    draw!(robot)
  end

  function main(time_interval=3.0; is_test=false)
    w = World(-5.0, 5.0, -5.0, 5.0)

    # create true map
    m = Map()
    landmark_positions = [[-4.0, 2.0],[2.0, -3.0],[3.0, 3.0],[0.0, 4.0], [1.0, 1.0], [-3.0, -1.0]]
    for (idx, pos) in enumerate(landmark_positions)
      add_object(m, Object(pos[1], pos[2], id=idx))
    end

    s = PsiSensor(m, dist_noise_rate=0.1, dir_noise=pi/90,
                  dist_bias_rate_stddev=0.1, dir_bias_stddev=pi/90)

    init_pose = [0.0, -3.0, 0.0]
    a = LoggerAgent(0.2, 5.0/180*pi, time_interval=time_interval,
                    init_pose=init_pose, is_test=is_test)
    r = DifferentialWheeledRobot(init_pose, 0.2, "black",
                                 a, time_interval,
                                 noise_per_meter=5,
                                 noise_std=pi/30,
                                 bias_rate_stds=[0.1, 0.1],
                                 sensor=s)

    if is_test == false
      anime = @animate for t in 0:time_interval:180
        animate_per_time(t, w, m, r)
      end
      path = "src/slam/graph_based_slam/traj_obsrv_edge_input_log.gif"
      gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], path))
      close_log_file(a)
    else
      for t in 0:time_interval:90
        animate_per_time(t, w, m, r)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# ---
# title: 172. Factorial Trailing Zeroes
# id: problem172
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/factorial-trailing-zeroes/description/>
# hidden: true
# ---
#
# Given an integer `n`, return _the number of trailing zeroes in`n!`_.
#
# **Follow up:** Could you write a solution that works in logarithmic time
# complexity?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3
#     Output: 0
#     Explanation:  3! = 6, no trailing zero.
#
#
# **Example 2:**
#
#
#
#     Input: n = 5
#     Output: 1
#     Explanation:  5! = 120, one trailing zero.
#
#
# **Example 3:**
#
#
#
#     Input: n = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 104`
#
#
## @lc code=start
using LeetCode

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


text:
module XConv

# Dependencies
using Random
using LinearAlgebra, IterTools
using CUDA, GPUArrays
using NNlib
using ChainRulesCore

import NNlib: colmajor

export grad_ev, initXConv

# Setup valuesa
const DEFAULT = "TrueGrad"
const EV = "EVGrad"
const _params = Dict(:p_size => 2^4, :mode => DEFAULT)

# Probing
include("gemm.jl")
include("probe.jl")

# EV setup function
function initXConv(probe_size::Integer, mode::String)
    _params[:mode] = mode
    _params[:p_size] = probe_size
end

# rrule
function ChainRulesCore.rrule(::typeof(NNlib.conv), X, w, cdim::DenseConvDims; kw...)
    Y = conv(X, w, cdim)
    back = _params[:mode] == EV ? Δconv_ev(X, w, cdim; kw...) : Δconv_std(X, w, cdim; kw...)
    return Y, back
end

function Δconv_std(x, w, cdim; kw...)
    function back(Δ)
        Δ = colmajor(Δ)
        return (
            NoTangent(),
            @thunk(∇conv_data(Δ, w, cdim, kw...)),
            @thunk(∇conv_filter(x, Δ, cdim, kw...)),
            NoTangent(),
        )
    end
    return back
end

function Δconv_ev(X, w, cdim; kw...)
    seed, eX = probe_X(X)
    function back(Δ)
        Δ = colmajor(Δ)
        return (
            NoTangent(),
            @thunk(∇conv_data(Δ, w, cdim, kw...)),
            @thunk(grad_ev(seed, eX, Δ, w, NNlib.stride(cdim))),
            NoTangent(),
        )
    end
    return back
end

end

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv