text:
export updateSettings!,updateSettingsMod!
export Splitdef,Rulepath,Leaf,Node,Tree,BoostedTree,SplittingCriterion,CVOptions,BaggedTree
export ModelSettings,copySettingsToCurrentType

#struct pdaMod end #this is a legacy type. It is only defined here, because we have some 'old' unused functions in the code which have a singature with this type (and the signature needs an update....)

#Splitting Criterions use multiple dispatch
	abstract type SplittingCriterion  end
	struct DifferenceSplit      <: SplittingCriterion end
    struct msePointwiseSplit      <: SplittingCriterion end
    struct sseSplit <:SplittingCriterion end
    struct mseSplit <:SplittingCriterion end
    struct NormalDevianceDifferenceToMeanFitWEIGHTEDSplit <:SplittingCriterion end
	struct PoissonDevianceSplit      <: SplittingCriterion end
	struct GammaDevianceSplit      <: SplittingCriterion end
	struct MaxValueSplit      <: SplittingCriterion end
	struct MaxAbsValueSplit      <: SplittingCriterion end
	struct MaxMinusValueSplit      <: SplittingCriterion end
	#struct RankOptSplit         <: SplittingCriterion end
	#struct ROptMinRLostPctSplit	   <: SplittingCriterion end
	#struct ROptMinRLostSplit	   <: SplittingCriterion end

	PoissonOrGamma = Union{PoissonDevianceSplit,GammaDevianceSplit}

#sortby options for categorical splits
	abstract type SortBy end
	struct SortByMean 	<: SortBy end

function canbeint(a)
	try convert(Int,parse(Float64,a))
	   return true
	catch
		return false
	end
end

function canbefloat(a)
	try Float64(a)
	   return true
	catch
		return false
	end
end

mutable struct Chart
	sheet::String
	chartDict::Dict{AbstractString,Dict{AbstractString,Any}}
	location::String
end

mutable struct ExcelSheet
	name::String #I think this name is currently not used anywhere
	data::DataFrame
end

mutable struct ExcelData
	sheets::Array{ExcelSheet,1}
    charts::Array{Chart,1}
    function ExcelData()
        return new(Array{ExcelSheet}(undef,0),Array{Chart}(undef,0))
    end
    function ExcelData(a,b)
        return  new(a,b)
    end
end

#errorStats
struct ErrorStats
	mse::Float64
	mae::Float64
	mrse::Float64
	mrae::Float64
	pearsoncorrelation::Float64
	function ErrorStats()
		new(NaN,NaN,NaN,NaN,NaN)
	end
	function ErrorStats(a::Float64,b::Float64,c::Float64,d::Float64,e::Float64)
		new(a,b,c,d,e)
	end
end

#subsets which are considered by splits
abstract type DTSubsets end

#custom loop over subsets (Gray Code) ONLY UNTIL the middle - USING SYMMETRY of the list as we consider binary splits
	struct MyGrayCodeSubsetsHALF <: DTSubsets
		size_of_set::Int
		stopat::Int #(1 << it.size_of_set)-2
	end

	eltype(it::MyGrayCodeSubsetsHALF) = Int #Array{eltype(1),1}
	length(it::MyGrayCodeSubsetsHALF) = it.stopat+1 #1 << it.size_of_set
	start(it::MyGrayCodeSubsetsHALF) = zero(1) #state = simply count through the subsets
	done(it::MyGrayCodeSubsetsHALF,state) = state>it.stopat  #if we have equality here, then the full set will be shifted / this is not needed for a binary split, hence we stop beforehand

	function next(it::MyGrayCodeSubsetsHALF, state)
	#returns the element which needs to be flipped
	  bit_to_flip=1+trailing_ones(state) #in the function print_graycode bit_to_flip is defined without the plus 1
	  state+=one(1)
	  bit_to_flip,state
	end
    @inline iterate(x::MyGrayCodeSubsetsHALF) = next(x, start(x))
    @inline iterate(x::MyGrayCodeSubsetsHALF, i) = done(x, i) ? nothing : next(x, i)

	bitflip_graycode_subsetsHALF(xs) = MyGrayCodeSubsetsHALF(length(xs),(1 << (length(xs)-1))-2)

#normal subset list for numeric splitting points ("increasing" subset)
	struct MyIncreasingSubsets <: DTSubsets
		size_of_set::Int
	end

	eltype(it::MyIncreasingSubsets) = Array{eltype(1),1}
	length(it::MyIncreasingSubsets)=it.size_of_set
	start(it::MyIncreasingSubsets)=zero(1) #state = simply count through the subsets
	done(it::MyIncreasingSubsets,state)=state>=it.size_of_set-1 #we do not need to flip/change the last element (as this would resutl in an empty left child)
    @inline iterate(x::MyIncreasingSubsets) = next(x, start(x))
    @inline iterate(x::MyIncreasingSubsets, i) = done(x, i) ? nothing : next(x, i)

	function next(it::MyIncreasingSubsets,state)
		state+=one(1)
		return copy(state),state
	end
	increasing_subsets(xs)=MyIncreasingSubsets(length(xs))

#in constrast to types, immutables cannot (or should not?) be changed over time
#however immutables should increase performance (todo: need to check this or read more about it)
mutable struct DTMTable
	key::Vector{String}
	trnidx::Vector{Int}
	validx::Vector{Int}
	numerator::Vector{Float64}
	denominator::Vector{Float64}
	weight::Vector{Float64}
	features::DataFrame
	candMatWOMaxValues::Vector{Vector{Float64}}
	mappings::Vector{Vector{String}}
	function DTMTable(key,trnidx,validx,numerator,denominator,weight,features,cands,maps)
		@assert length(validx)+length(trnidx)==length(key)
		@assert length(key)==length(numerator)
		@assert length(key)==length(denominator)
		@assert length(key)==length(weight)
		@assert length(key)==size(features,1)
		@assert length(cands)+length(maps)==size(features,2)
		@assert length(findall((in)(validx),trnidx))==0 "DTM: trnidx and validx are not mutually exclusive"
		return new(key,trnidx,validx,numerator,denominator,weight,features,cands,maps)
	end
end
==(x::DTMTable,y::DTMTable)= (x.key==y.key) && (x.trnidx==y.trnidx) && (x.validx==y.validx) && (x.numerator==y.numerator) && (x.denominator==y.denominator)  && (x.weight==y.weight)  && (x.features==y.features) && (x.candMatWOMaxValues==y.candMatWOMaxValues) && (x.mappings==y.mappings)

Base.size(d::DTMTable) = size(d.features)
hash(x::DTMTable)=hash(x,0)
function hash(x::DTMTable,h)
    h=hash(x.key,h)
    h=hash(x.trnidx,h)
    h=hash(x.validx,h)
    h=hash(x.numerator,h)
    h=hash(x.denominator,h)
    h=hash(x.weight,h)
    h=hash(x.features,h)
    h=hash(x.candMatWOMaxValues,h)
    h=hash(x.mappings,h)
    return h
end

function getindex_internal(r::DTMTable,idx;compact_features=false) where T <: Number
    #this is currently not working for boolean indices (BitVector)
@assert !compact_features
key=r.key[idx]
numerator=r.numerator[idx]
denominator=r.denominator[idx]
weight=r.weight[idx]
features=r.features[idx,:]

#redefine trn and validx
#version with intermediate arrays
#found=findall((in)(idx),trnidx)
#integersNew=trnidx[found]
#trnidxNEW=findall((in)(integersNew),idx)

#version with fewer intermediate arrays
found=findall((in)(idx),r.trnidx)
found=r.trnidx[found]
trnidxNEW=findall((in)(found),idx)

found=findall((in)(idx),r.validx)
found=r.validx[found]
validxNEW=findall((in)(found),idx)

(length(validxNEW)!=0)||@warn("DTM: Validation index of the DTMTable has length zero! You may need to redefine validx. The function resample_trnvalidx! might be helpful.")
(length(trnidxNEW)!=0)||@warn("DTM: Training index of the DTMTable has length zero! You may need to redefine trnidx. The function resample_trnvalidx! might be helpful.")

#I think this is generally not desirable (e.g. when a dtmtable is created based on a large table, then models are fit to subsets of the data: in that case we want the pools to match perfectly!)
if compact_features
    @assert false
    @warn("for this to work one would need to update mappings and candMatWOMaxValues too!")
    tmp_number_of_num_features=sum(map(x->eltype(r.features[x].pool),1:size(r.features,2)).!=String)
    #"compact" features, it is possible that we could do this in a better/faster way
    for i=tmp_number_of_num_features+1:size(features,2)
        if eltype(features[i])!=UInt8
            features[i]=DecisionTrees.PooledArraysDTM.PooledArray(Array(features[i]))
        end
    end
end

dt=DTMTable(key,trnidxNEW,validxNEW,numerator,denominator,weight,features,deepcopy(r.candMatWOMaxValues),deepcopy(r.mappings))
return dt
end

function Base.getindex(r::DTMTable,idx::UnitRange{T};compact_features=false) where T <: Number
    return getindex_internal(r,idx,compact_features=compact_features)
end

function Base.getindex(r::DTMTable,idx::Vector{T};compact_features=false) where T <: Number
    return getindex_internal(r,idx,compact_features=compact_features)
end

#todo check the difference between immutable and type!! are we using immutables correctly here?
struct Splitdef{T<:Unsigned}
	featid::Int
	featid_new_positive::Int
    featurename::Symbol
    subset::Array{T,1}
    splitvalue::Float64 #Depends on the criterion function used (e.g. difference)
    weightl::Float64
    weightr::Float64
end
==(x::Splitdef,y::Splitdef)= (x.featid==y.featid) && (x.featurename==y.featurename) && (x.featid_new_positive==y.featid_new_positive) && (x.splitvalue==y.splitvalue) && (x.subset==y.subset)  && (x.weightl==y.weightl)  && (x.weightr==y.weightr)
#Splitdef(x::Tuple{Int,Float64,Array{UInt8,1},Float64,Float64,Float64})= Splitdef(x[1],x[2],x[3],x[4],x[5],x[6])
hash(x::Splitdef)=hash(x,UInt(9))
function hash(x::Splitdef,h::UInt)
  h=hash(x.featid,h)
  h=hash(x.featid_new_positive,h)
  h=hash(x.featurename,h)
  h=hash(x.subset,h)
  h=hash(x.splitvalue,h)
  h=hash(x.weightl,h)
  h=hash(x.weightr,h)
  return h
end

struct Rulepath{T<:Unsigned}
    featid::Int
    subset::Vector{T}
    isLeftChild::Bool
end

Rulepath{T}() where {T<:Unsigned} = Rulepath(0,T[],false)

==(x::Rulepath,y::Rulepath)= (x.featid==y.featid) && (x.isLeftChild==y.isLeftChild) && (x.subset==y.subset)
hash(x::Rulepath)=hash(x,UInt(9))
function hash(x::Rulepath,h::UInt)
  h=hash(x.featid,h)
  h=hash(x.isLeftChild,h)
  for i in x.subset
    h = hash(i, h)
  end
  return h
end

mutable struct Leaf{T<:Unsigned} #leafs are mutable, since we determine the number of the leaf a posteriori
    #NOTE: the idx array can have an unexpected meaning/size in the case where we perform subsampling!
	rowcount::Int #number of observations of the training data set that fell into thta leaf
    mean::Float64
    fitted::Float64 #this is usually identical to the mean, but can be different for some instances
  	size::Float64 #  this is equal to sum(weight)
    depth::Int
    rule_path::Array{Rulepath{T},1}
	sumnumerator::Float64
	sumdenominator::Float64
    id::Int #a number 1<=id<=#number of leaves in the tree
end
function Leaf{T}() where T
        rp=Vector{Rulepath{T}}(undef,0)
        return Leaf{T}(0,0.0,0.0,0.0,0,rp,0.0,0.0,0)
    end
function Leaf(rowc,me,fi,sz,depth,rp::Vector{Rulepath{T}},sumn,sumd,id) where T
    return Leaf{T}(rowc,me,fi,sz,depth,rp,sumn,sumd,id)
end
==(x::Leaf,y::Leaf)= (x.rowcount==y.rowcount) && (x.mean==y.mean) && (x.fitted==y.fitted)  && (x.size==y.size)  && (x.depth==y.depth)  && (x.rule_path==y.rule_path)  && (x.sumnumerator==y.sumnumerator)  && (x.sumdenominator==y.sumdenominator) && (x.id==y.id)
hash(x::Leaf)=hash(x,UInt(9))
function hash(x::Leaf,h::UInt)
  h=hash(x.rowcount,h)
  h=hash(x.mean,h)
  h=hash(x.fitted,h)
  h=hash(x.size,h)
  h=hash(x.depth,h)
  h=hash(x.rule_path,h)
  h=hash(x.sumnumerator,h)
  h=hash(x.sumdenominator,h)
  h=hash(x.id,h)
  return h
end

struct Node{T<:Unsigned}
    #todo, add depth!
    featid::Int
    featid_new_positive::Int
    #we have redundance in the tree structure, the rulepath contains the same information as the members of the node; todo - remove it!
    subset::Array{T,1}
    left::Union{Leaf{T},Node{T}} #,Node{UInt16}}
    right::Union{Leaf{T},Node{T}} #,Node{UInt16}}
    rule_path::Array{Rulepath{T},1}
end

mutable struct ModelSettings
	model_type::String #1
	minWeight::Float64 #2
	randomw::Float64 #3
	crit::SplittingCriterion # fn version of 4
	maxSplittingPoints::Int #5
	iterations::Int #6
	learningRate::Float64 #7 moderationfactor
	nScores::Int #8
	adaptiveLearningRate::Float64 #9
	prem_buffer::Int #11
	startAtMean::Bool #12
	writeTree::Bool #13
	number_of_num_features::Int #14
	spawnsmaller::Bool #26
    boolRankOptimization::Bool #35
	boolRandomizeOnlySplitAtTopNode::Bool #37
	subsampling_prop::Float64 #38 #NOTE Friedman suggests to sample without replacement see Friedman 2002 Stochastic Gradient Boosting,” Computational Statistics and Data Analysis 38(4):367-378
	subsampling_features_prop::Float64 #39
	version::String #40
	preppedJLDFileExists::Bool #41
	catSortByThreshold::Int #42
	catSortBy::SortBy #43
	scorebandsstartingpoints::Array{Int,1} #44
	showTimeUsedByEachIteration::Bool #45
	smoothEstimates::String #46
	deriveFitPerScoreFromObservedRatios::Bool #deriveFitPerScoreFromObservedRatios::Bool (if this is true then the fit per score will be based on the observed ratio per score (instead of the fitted ratio per score)) . Note that this does not influence the structure of the tree, but only the fitted ratio per score (and scoreband)
    roptForcedPremIncr::Bool #47
    writeSasCode::Bool #49
	writeIterationMatrix::Bool #50
	writeResult::Bool #51
	writeStatistics::Bool #52
	boolCreateZipFile::Bool #53
	writeCsharpCode::Bool
	writeVbaCode::Bool
	nDepthToStartParallelization::Int
	baggingWeightTreesError::String
	cBB_niterBoosting::Int
	cBB_niterBagging::Int
	fixedinds::Array{Int,1}
	boolNumeratorStats::Bool
	statsByVariables::Array{Int,1}
	statsRandomByVariable::Int
	saveJLDFile::Bool #if this is set to false no *.jld2 file is saved #todo we should rename this to ensure it is clear that this refers to the prepped data only!
	saveResultAsJLDFile::Bool #this refers to the model output/result
	print_details::Bool #whether julia should print details to the log or not (note to developers: warnings/errors and certain important infos should always be printed!)
	seed::Int #optional entropy
	graphvizexecutable::String #e.g. C:\Program Files (x86)\Graphviz2.38\ if this
	showProgressBar_time::Bool # bar form iterators package (or some other package)
	prroduceEstAndLeafMatrices::Bool
	write_dot_graph::Bool
	calculateGini::Bool #whether or not to calculate the gini (which requires sorting of the data (which takes time!))
	calculatePoissonError::Bool
    performanceMeasure::String
	fitForStatsAndCharts::String
	ignoreZeroDenominatorValues::Bool

	#the following are treated specially
	df_name_vector::Array{String,1}
    number_of_char_features::Int #this is calculated by Julia from number_of_num_features and ncolsdfIndata (which is defined by the data)

	chosen_apply_tree_fn::String #defined below
	moderationvector::Array{Float64,1} #50 onwards

  function ModelSettings()
	#initialize default settings
	model_type="build_tree" #1
	minWeight=-0.1 #2
	randomw=0.0 #3
	crit=DifferenceSplit()	 #4
	maxSplittingPoints=250 #5 #arbitrary choice. Preliminary testing has shown no performance difference for different values (say 10 or 200) for this setting
	iterations=2 #6
	learningRate=0.1 #7
	nScores=1000 #8
	adaptiveLearningRate=1.0 #9
	prem_buffer=0 #11
	startAtMean=true #12
	writeTree=true #13
	number_of_num_features=-1 #14
	spawnsmaller=true #26
	boolRankOptimization=false #35
	boolRandomizeOnlySplitAtTopNode=true #37
	subsampling_prop=1.0 #38
	subsampling_features_prop=1.0 #39
	version="not_initialized" #40
	preppedJLDFileExists=false #41
	catSortByThreshold=7 #42
	catSortBy=SortByMean() #43
	scorebandsstartingpoints=[parse(Int,x) for x in split("1,100,200,300,400,500,600,700,800,900",',')] #scorebandsstartingpoints=parse(Int,split("1,100,200,300,400,500,600,700,800,900",',')) #44
	showTimeUsedByEachIteration=false #true #45
	smoothEstimates="1" #46 this is a string for now (as there could be different smoothing methods specified by this string)
	deriveFitPerScoreFromObservedRatios=true
	roptForcedPremIncr=false #47
    writeSasCode=false #49
	writeIterationMatrix=false#50
	writeResult=false #true #51
	writeStatistics=true #52
	boolCreateZipFile=false #true #53
	writeCsharpCode=false #true
	writeVbaCode=false
	nDepthToStartParallelization=-1
	baggingWeightTreesError="mae" # ::String
	#combined bagging&boosting model
	cBB_niterBoosting=0
	cBB_niterBagging=0
	fixedinds=Array{Int}(undef,0)
	boolNumeratorStats=false
	statsByVariables=Int[]
	statsRandomByVariable=5
	saveJLDFile=true
	saveResultAsJLDFile=false
	print_details=true
	seed=90210
	graphvizexecutable="" #C:\\Program Files (x86)\\Graphviz2.38\\bin\\dot.exe""
	showProgressBar_time=true
	prroduceEstAndLeafMatrices=false
	write_dot_graph=false
	calculateGini=false
	calculatePoissonError=false
    performanceMeasure="Lift Val"
	fitForStatsAndCharts="rawRelativities" #rawRelativities,unsmoothedPerScore,smoothedPerScore
	ignoreZeroDenominatorValues=false

	#the following are treated specially
	df_name_vector=Array{String}(undef,0) #
    number_of_char_features=-1 #
	chosen_apply_tree_fn="apply_tree_by_leaf" # #apply_tree_by_row does not seem to work for (certain?) boosting models
	moderationvector=[0.1] #

	return new(model_type,minWeight,randomw,crit,maxSplittingPoints,iterations,learningRate,nScores,adaptiveLearningRate,prem_buffer,startAtMean,writeTree,number_of_num_features,spawnsmaller,boolRankOptimization,boolRandomizeOnlySplitAtTopNode,subsampling_prop,subsampling_features_prop,version,preppedJLDFileExists,catSortByThreshold,catSortBy,scorebandsstartingpoints,showTimeUsedByEachIteration,smoothEstimates,deriveFitPerScoreFromObservedRatios,roptForcedPremIncr,writeSasCode,writeIterationMatrix,writeResult,writeStatistics,boolCreateZipFile,writeCsharpCode,writeVbaCode,nDepthToStartParallelization,baggingWeightTreesError,cBB_niterBoosting,cBB_niterBagging,fixedinds,boolNumeratorStats,statsByVariables,statsRandomByVariable,saveJLDFile,saveResultAsJLDFile,print_details,seed,graphvizexecutable,showProgressBar_time,prroduceEstAndLeafMatrices,write_dot_graph,calculateGini,calculatePoissonError,performanceMeasure,fitForStatsAndCharts,ignoreZeroDenominatorValues
	 ,df_name_vector,number_of_char_features,chosen_apply_tree_fn,moderationvector)
  end  # ModelSettings()

end #type definition

function edit_statsByVariables(col,stringValue,s::ModelSettings)
if col=="statsbyvariables"&&length(stringValue)>0
	if in(',',stringValue)||in(';',stringValue)
		println("statsbyvariables=$(stringValue)")
		error("DTM: statsbyvariables must be space separated not comma separated!")
	end
	try
	chosenvars=[lowercase(x) for x in split(stringValue,' ')]
	integerindices=findall(in(chosenvars), [lowercase(x) for x in s.df_name_vector])
	if length(integerindices)!=length(chosenvars)
		@warn("DMT: Some variables of statsbyvariables were not found in the data")
		@show s.df_name_vector
		@show chosenvars
		@show stringValue
		@show integerindices
		@show [lowercase(x) for x in s.df_name_vector]
		error("DTM: Abort")
		else
		integerindices2=[x<=s.number_of_num_features ? x : -1*(x-s.number_of_num_features) for x in integerindices]
		return string(integerindices2)[2:end-1] #remove leading and trailing '[' and ']'
	end
	catch er
		@show er
		error("DTM: Unable to parse statsbyvariables option")
	end
end
return stringValue
end

function copySettingsToCurrentType(oldSetting)
	s=ModelSettings()
	newfields=fieldnames(typeof(s))
	oldfields=fieldnames(typeof(oldSetting))

	for field in fieldnames(typeof(oldSetting))
		if in(field,newfields)
			v=getfield(oldSetting,field)
			setfield!(s,field,v)
		end
    end

	#consider new and dropped fields
		newfs=setdiff(newfields,oldfields)
		droppedfs=setdiff(oldfields,newfields)
	@info "DTM: Copied settings to new type."
	if length(newfs)>0
		@info "DTM: New type has the following new fields:"
		for f in newfs
			v=getfield(s,f)
			@show f,v
		end
	end
	if length(droppedfs)>0
		@info "DTM: Old type had the following fields which were dropped:"
		for f in droppedfs
			v=getfield(oldSetting,f)
			@show f,v
		end
	end
    return s
end

"""
updateSettings!(s::ModelSettings;args...)
use this function to update the model settings
try for instance:

updateSettings!(sett,minWeight=-0.3,iterations=30)
"""
function updateSettings!(s::ModelSettings;args...)
	updateSettingsMod!(s;args...)
end

"""
updateSettingsMod!(s::ModelSettings;args...)
use this function to update the model settings
try for instance:

updateSettingsMod!(sett,minWeight=-0.3,iterations=30)
"""
function updateSettingsMod!(s::ModelSettings;args...)
	seen=[]
	for (smember,v) in args
		if in(smember,seen)
			error("DTM: You provided multiple values for the field $(smember)")
		end
		push!(seen,smember)

		oldValue=getfield(s,smember)
		if typeof(v)!=typeof(oldValue)
			valConverted=convertFromString(oldValue,string(v)) #tbd maybe string(v) could/should be replaced with v here
			setfield!(s,smember,valConverted)
		else
			setfield!(s,smember,v)
		end
	end
	res=checkIfSettingsAreValid(s)
	@assert res==nothing "Some settings are invalid! Abort."
end

myParseAsFloat(val::T) where {T<:AbstractString} = parse(Float64,val)
myParseAsFloat(val::T) where {T<:Number} = float(val)
convertFromString(oldvalue::T,val) where {T <: Any}=convert(T,val) #generic method catchall
convertFromString(oldvalue::T,val) where {T <: AbstractString}=convert(T,string(val))

#convertFromString(oldvalue::T,val::U) where {T <: AbstractFloat,U<:AbstractString}=parse(Float64,val)
convertFromString(oldvalue::T,val) where {T <: AbstractFloat}=myParseAsFloat(val)
convertFromString(oldvalue::T,val) where {T <: Integer}=convert(Integer,myParseAsFloat(val))
#convertFromString(oldvalue::T,val::U) where {T <: Integer,U<:AbstractString}=convert(Integer,parse(Float64,val))
#convertFromString(oldvalue::T,val::U) where {T <: Unsigned,U<:AbstractString}=uint(parse(Float64,val))
convertFromString(oldvalue::T,val) where {T <: Unsigned}=uint(myParseAsFloat(val))
convertFromString(oldvalue::T,val) where {T <: Bool}=(lowercase(val)=="t"||lowercase(val)=="true"||val=="1")
function convertFromString(oldvalue::T,val) where {T <: SortBy}
	@assert lowercase(val)=="sortbymean"||lowercase(val)=="mean"
	return SortByMean()  #currently only mean is possible here
end
function convertFromString(oldvalue::T,critfnstr) where {T <: SplittingCriterion}
	local allowedinputstrings=[lowercase(x) for x in ["difference" "mse" "msePointwise" "poisson" "gamma" "maxvalue" "maxabsvalue" "maxminusvalue"]]
	critfnstr=lowercase(critfnstr)
	if (critfnstr=="difference"||critfnstr=="default")
		crit=DifferenceSplit()
	elseif (critfnstr=="msepointwise") #||critfnstr=="normaldeviancepointwise" ||critfnstr=="normalpointwise" ||critfnstr=="gaussianpointwise" || critfnstr=="gaussiandeviancepointwise" || critfnstr=="rsspointwise")
		crit=msePointwiseSplit()
    elseif (critfnstr=="mse")
        crit=mseSplit()
    elseif (critfnstr=="sse")
        crit=sseSplit()
	elseif (critfnstr=="poissondeviance"||critfnstr=="poisson")
		crit=PoissonDevianceSplit()
	elseif (critfnstr=="gammadeviance"||critfnstr=="gamma")
		crit=GammaDevianceSplit()
	elseif critfnstr=="maxvalue" #pick the split which leads to the maximum average value (in either the right or left child)
		crit=MaxValueSplit()
	elseif critfnstr=="maxabsvalue" #pick the split which leads to the absolute maximum average value (in either the right or left child)
		crit=MaxAbsValueSplit()
	elseif critfnstr=="maxminusvalue" #pick the split which leads to the maximum negative value (in either the right or left child)
		crit=MaxMinusValueSplit()
	#elseif critfnstr=="roptminrlostpct"
	#	crit=ROptMinRLostPctSplit()
	#elseif critfnstr=="roptminrlost"
	#	crit=ROptMinRLostSplit()
	else
		@info "Invalid splitting criterion: $(critfnstr). Currently only the following strings are allowed as input:"
		@info string(allowedinputstrings)
		error("DTM: Abort due to invalid splitting criterion (see above).")
	end
	return crit
end

function convertFromString(oldvalue::Array{Int,1},val)
local arr
    try
		arr=Int[parse(Int,x) for x in split(val,',')]
	catch er
		@show er
		error("DTM: Unable to parse input string as integer vector: input= $(val)")
	end
return arr
end

function convertFromString(oldvalue::Array{Float64,1},val)
local arr
    try
		arr=Float64[parse(Float64,x) for x in split(val,',')]
	catch er
		@show er
		error("DTM: Unable to parse input string as integer vector: input= $(val)")
	end
return arr
end

function checkForMandatorySettings(headerLowercase)
	@assert in("number_of_num_features",headerLowercase) global_number_of_num_f_warning_mandatory_field
	return nothing
end

function checkIfSettingsAreValid(s::ModelSettings)
	  #@assert in(method,globalConstAllowableMethodsForDefineCandidates)
	  @assert in(s.fitForStatsAndCharts,globalValidfitForStatsAndCharts)
      @assert s.statsRandomByVariable<typemax(UInt8) #more than 255 random groups is really not meaningful; we do work with a UInt8 list later on
	  @assert in(s.model_type,["boosted_tree" "build_tree" "bagged_tree" "bagged_boosted_tree"])
      @assert (s.randomw>=0) & (s.randomw<=1)
	  #Note, we could easily support more splitting Points (i.e. UInt), but it is likely that the user has not intended to exceed typemax(UInt16), thus we catch this here
	  max_splitting_points_num_INTERNAL=Int(typemax(UInt16))-1
      @assert s.maxSplittingPoints<max_splitting_points_num_INTERNAL "Maximum number of splitting points is currently limited to $(max_splitting_points_num_INTERNAL)"
	  @assert s.maxSplittingPoints>0
	  if (s.model_type=="bagged_boosted_tree")
		@assert min(cBB_niterBagging,cBB_niterBoosting)>0
		@assert abs(subsampling_prop)<1.0 #without subsampling bagged boosting probably does not make too much sense
		@assert subsampling_features_prop<1.0 #without sampling features bagged boosting probably does not make too much sense
		@assert s.iterations==0 "Please set iterations to 0 for Combined Bagging&Boosting Models" #I want to avoid confusion where the use thinkgs that this parameter is of relevance.
	  else
		@assert s.iterations>0
	  end
	  @assert in(s.baggingWeightTreesError,["mse" "pearsoncorrelation" "mrae" "mrse" "mae" "uniform"])
	  #if length(s.moderationvector)=1
	  if in(s.model_type,["build_tree"])
		s.moderationvector=[s.learningRate]
	  end
	  if length(s.moderationvector)==1
		if s.moderationvector[1]!=s.learningRate && s.model_type=="boosted_tree"
			@info "Replacing moderationvector[1] with the value of learningRate=$(s.learningRate)"
			s.moderationvector[1]=s.learningRate
		end
	  end
	  if s.model_type=="boosted_tree"
		if !((length(s.moderationvector)==1)&&(s.moderationvector[1]==s.learningRate))
			@show s.moderationvector
			@show s.learningRate
			error("DTM: Invalid learningRate or moderationvector settings")
	  	end
		end
      @assert (s.learningRate>=0 && s.learningRate<=1)
	  #end
      @assert (s.adaptiveLearningRate<=1.0 && s.adaptiveLearningRate>=0.0)
      #@assert s.subsampling_prop<=1.0 #positive value => without replacement
	  @assert s.subsampling_prop>=-1.0 #negative value => with replacement, smaller than -1 is not possible wherease larger than 1 is possible
	  if s.subsampling_prop<1.0
		@warn("Subsampling of the data is enabled. Ensure the value of subsampling_prop=$(s.subsampling_prop) has a meaningful value relative to minWeight=$(s.minWeight)")
	  end
	  #if (s.model_type!="bagged_tree")&&(abs(s.subsampling_prop)!=1.0);@warn("Subsampling with Boosting methods is experimental!");end
	  #todo/tbd need to fix this! see goodnessOfFit
		if (s.model_type=="bagged_tree")&&(abs(s.subsampling_prop)!=1.0);@warn("Bagging is currently EXPERIMENTAL!") ;end
      #@assert (abs(s.subsampling_prop)==1.0)||(s.model_type=="bagged_tree") "ERROR: (bk) Sumbsampling is currently only available for Bagging, as it will require the usage of apply tree after each iteration (since not all estimates are available after the construction of the tree)! Thus (for now) you need to set subsampling_prop to 1.0."
      @assert (abs(s.subsampling_features_prop>0) && abs(s.subsampling_features_prop<=1))
	  @assert s.scorebandsstartingpoints[1]==1
	  @assert issorted(s.scorebandsstartingpoints)
	  @assert s.scorebandsstartingpoints[end]<s.nScores
	  @assert s.nScores>8 "Please set nScores to a value greater than 8" #Less than 8 scores makes it hard/impossible for the moving average
	  @assert in(s.smoothEstimates,["0" "1"]) #this is currently a string since there could (in the future) be multiple smoothing methods
      #check graphvizexecutable

      if s.write_dot_graph
        gcmd = `$(s.graphvizexecutable) "-V"`
        try
            tmptxt = mktemp()[1]
            tmperr = mktemp()[1]
            pip = pipeline(gcmd, stdout=tmptxt, stderr=tmperr)
            gres = run(pip)
            if gres.exitcode != 0
                @warn("Please review the graphvizexecutable path in the settings")
                @show s.write_dot_graph
                @show s.graphvizexecutable
                @show gcmd
                @show gres.exitcode
            end
        catch
            @warn("Please review the graphvizexecutable path in the settings")
            println("Unable to issue command:\r\n $(gcmd)")
            @show s.write_dot_graph
            @show s.graphvizexecutable
        end
    end # if s.write_dot_graph

	if s.roptForcedPremIncr
		#@assert (s.crit==ROptMinRLostPctSplit()||s.crit==ROptMinRLostSplit())
	end
	@assert !s.boolRankOptimization "this is currently not working in this version of the Code"
	#@assert s.ishift==0 #is zero for normal models (historically this was 1 for rkoptmodels (labels_orig and so on....)
	#@show s.df_name_vector
	#@show s.number_of_char_features,s.number_of_num_features
	if length(s.df_name_vector)!=(s.number_of_char_features+s.number_of_num_features)
		@show length(s.df_name_vector),s.number_of_char_features,s.number_of_num_features
		@assert length(s.df_name_vector)==(s.number_of_char_features+s.number_of_num_features)
	end
    @assert in(s.performanceMeasure,global_statsperiter_header)
    if (size(s.moderationvector,1)!=s.iterations && size(s.moderationvector,1)!=1)
		@assert false "DTM: This should not have happend: size(moderationvector) is neither 1 nor #iterations"
    end
    @assert s.startAtMean "DTM: This is currently not working: please set startAtMean=true)"
    #@assert !s.BoolAdaptiveLearningRate "DTM: this is not yet updated to reflect numerator/denominator objectives!!"
    @assert ((s.adaptiveLearningRate>=1.0)||(s.adaptiveLearningRate<=0.0)) "DTM: (adaptiveLearningRate is not in [0,1]). This is not yet updated to reflect numerator/denominator objectives!!"

     return nothing
end

abstract type DTModel end
mutable struct Tree <: DTModel
	rootnode::Union{Leaf,Node{UInt8},Node{UInt16}}
	intVarsUsed::Array{Array{Int,1},1}
	candMatWOMaxValues::Array{Array{Float64,1},1}
	charMappings::Array{Array{String,1},1}
	inds_considered::Array{Int,1}
	settings::ModelSettings
	variableImp1Dim::Array{Int,1}
	variableImp2Dim::Array{Int,1}
	modelstats::DataFrame # #overall model statistics and performance (trn and val), this needs to be a matrix (Any) with a header row NOTE: this is a feature of multiple model types! (tree, boosting, bagging) and we need consistency because of the function run_model_multirow_settings
	exceldata::ExcelData
	featurepools::Array{Union{Array{Float64,1},Array{String,1}},1}
	function Tree(rootnode::Union{Leaf,Node{UInt8},Node{UInt16}},intVarsUsed::Array{Array{Int,1},1},candMatWOMaxValues::Array{Array{Float64,1},1},charMappings::Array{Array{String,1},1},inds_considered::Array{Int,1},settings::ModelSettings,exceldata::ExcelData,fp::Array{Union{Array{Float64,1},Array{String,1}},1})
		nf=settings.number_of_char_features+settings.number_of_num_features
		@assert length(settings.df_name_vector)==nf
		variableImp1Dim=zeros(Int,nf)
		variableImp2Dim=zeros(Int,div(nf*(nf-1),2))
		ms = DataFrame()

		return new(rootnode,intVarsUsed,candMatWOMaxValues,charMappings,inds_considered,settings,variableImp1Dim,variableImp2Dim,ms,exceldata,fp)
	end
	function Tree()
        rootnode=Leaf{UInt8}()
        intVarsUsed=Int[]
        candMatWOMaxValues=[Float64[]]
        charMappings=[String[]]
        inds_considered=Int[]
        settings=ModelSettings()
        variableImp1Dim=Int[]
        variableImp2Dim=Int[]
        ms=DataFrame()
        exceldata=ExcelData()
        fp=[Float64[]]
		return new(rootnode,intVarsUsed,candMatWOMaxValues,charMappings,inds_considered,settings,variableImp1Dim,variableImp2Dim,ms,exceldata,fp)
	end
end

#this type was only created, such that the @distributed results can be created more easilty (with a vcat operation)
#there would have been other workarounds which would have worked too
mutable struct TreeWithErrorStats
	tree::Tree
	err::ErrorStats
end

abstract type Ensemble <: DTModel end

struct EmtpyDTModel <:DTModel
    modelstats::DataFrame #this might not be needed
    function EmtpyDTModel()
        return new(DataFrame())
    end
end

#This is currently only used if multi run fails (as a return value of an 'empty' model)
struct EmptyEnsemble <: Ensemble
    modelstats::DataFrame #this might not be needed
    function EmptyEnsemble()
        return new(DataFrame())
    end
end

struct BoostedTree <: Ensemble
    trees::Vector{Union{Leaf,Node{UInt8},Node{UInt16}}} #todo tbd make this to trees!

	settings::ModelSettings
	intVarsUsed::Array{Array{Array{Int,1},1},1}
	candMatWOMaxValues::Array{Array{Float64,1},1}
	charMappings::Array{Array{String,1},1}
	inds_considered::Array{Array{Int,1},1}

    moderationvector::Array{Float64,1}
    scores::Array{Int,1}
    rawrelativities::Array{Float64,1}
    maxRawRelativityPerScoreSorted::Array{Float64,1}
    meanobserved::Float64
    startAtMean::Bool
	ScoreToSmoothedEstimate::Array{Float64,1}
	rawObservedRatioPerScore::Array{Float64,1}
	iterationmatrix::Array{Float64,2}
	modelstats::DataFrame #overall model statistics and performance (trn and val), this needs to be a matrix (Any) with a header row NOTE: this is a feature of multiple model types! (tree, boosting, bagging) and we need consistency because of the function run_model_multirow_settings
	exceldata::ExcelData
	trnidx_one_zero_full::Array{UInt,1}
	featurepools::Array{Union{Array{Float64,1},Array{String,1}},1}
end

struct BaggedTree <: Ensemble
    trees::Vector{Node} #todo tbd make this to trees! (? really)
	weightPerTree::Array{Float64,1}
	oobagErr::Array{ErrorStats,1}

	settings::ModelSettings
	intVarsUsed::Array{Array{Array{Int,1},1},1}
	candMatWOMaxValues::Array{Array{Float64,1},1}
	charMappings::Array{Array{String,1},1}
	inds_considered::Array{Array{Int,1},1}

    scores::Array{Int,1}
    rawrelativities::Array{Float64,1}
    maxRawRelativityPerScoreSorted::Array{Float64,1}
    meanobserved::Float64
    startAtMean::Bool
	ScoreToSmoothedEstimate::Array{Float64,1}
	rawObservedRatioPerScore::Array{Float64,1}
	iterationmatrix::Array{Float64,2}
	modelstats::DataFrame
	trnidx_one_zero_full::Array{UInt,1}
	featurepools::Array{Union{Array{Float64,1},Array{String,1}},1}
end

"""
cvo=CVOptions(folds, training_proportion, use_all_data)

folds is an integer. If folds<0 then we consider n disjoint training sets

    use_all_data::Bool #if set to false, the CVOptions will never use the 'original' validation rows in the data

    training_proportion::Float64 #size of trn data. The value must be >0 and <1

More precisely the evaluation is as follows

if cvo.training_proportion>0

\t   training_proportion > 0 provided. Performing random sampling without replacement 'cvo.folds' samples)

else

\t if cvo.folds<0

\t DTM: Performing 'DISJOINT training sets' k-fold cross validation 'k=cvo.folds'

\t else

\t\t #we have folds>=0
\t\t Performing k-fold cross validation (disjoint validation sets, 'k=cvo.folds')

\t end

end
"""
mutable struct CVOptions
	folds::Int #number of cv models
	#if folds<0 then we consider n disjoint training sets
	training_proportion::Float64 #size of trn data
	use_all_data::Bool #if false, the CVOptions will never use the 'original' validation rows in the data
	function CVOptions()
		x=new(10,0.7,true)
		@assert check_cvoptions(x)
		return x
	end
	function CVOptions(a::T,b::N,c::Bool) where {T<:Integer,N<:Number}
		x=new(Int(a),float(b),c)
		@assert check_cvoptions(x)
		return x
	end
end

function check_cvoptions(cvo::CVOptions)
	# folds::UInt16 #number of cv models
	# disjoint::Bool #whether we split the data into 'folds' disjoing subsets, if this is false then random sampling is performed
	# training_proportion::Float64 #size of trn data
	# use_all_data::Bool #if false, the CVOptions will never use the 'original' validation rows in the data
	#@assert cvo.folds>0
	# training_proportion == 1 could be possible but is a somewhat degenerate case
	@assert (cvo.training_proportion>=0 && cvo.training_proportion<1)
	return true
end

function Base.convert(::Type{Splitdef{UInt16}},x::Splitdef{UInt8})
    return Splitdef{UInt16}(x.featid,x.featid_new_positive,x.featurename,Vector{UInt16}(x.subset),x.splitvalue,x.weightl,x.weightr)
end


text:
"""
Pseudopotential correction energy. TODO discuss the need for this.
"""
struct PspCorrection end
(::PspCorrection)(basis) = TermPspCorrection(basis)

struct TermPspCorrection <: Term
    energy::Real  # precomputed energy
end
function TermPspCorrection(basis::PlaneWaveBasis)
    model = basis.model
    if model.n_dim != 3 && any(attype isa ElementPsp for (attype, _) in model.atoms)
        error("The use of pseudopotentials is only sensible for 3D systems.")
    end
    TermPspCorrection(energy_psp_correction(model.lattice, model.atoms))
end

function ene_ops(term::TermPspCorrection, basis::PlaneWaveBasis, ψ, occ; kwargs...)
    (E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
end

"""
Compute the correction term for properly modelling the interaction of the pseudopotential
core with the compensating background charge induced by the `Ewald` term.
"""
function energy_psp_correction(lattice, atoms)
    T = eltype(lattice)

    # Early return for cases without atoms or psp atoms
    any(attype isa ElementPsp for (attype, _) in atoms) || return T(0)

    # Total number of explicitly treated (i.e. valence) electrons
    n_electrons = sum(length(positions) * n_elec_valence(attype)
                      for (attype, positions) in atoms)

    correction_per_cell = sum(
        length(positions) * eval_psp_energy_correction(T, attype.psp, n_electrons)
        for (attype, positions) in atoms
        if attype isa ElementPsp
    )

    correction_per_cell / compute_unit_cell_volume(lattice)
end

ext: jl


text:
  """
  Module containing a gallery of examples of nonlinear eigenvalue problems.\\
  Look at the function 'nep_gallery()' for further instructions.
  """
module Gallery
    using NonlinearEigenproblems.Serialization
    using ..NEPCore
    using ..NEPTypes
    using Random
    using LinearAlgebra
    using SparseArrays
    using SpecialFunctions

    export nep_gallery

    function __init__()
        # Add the search-path to the extra galleries - so that unloaded modules can be easily loaded from top level
        this_path = string(@__DIR__, "/gallery_extra")
        if !(this_path in LOAD_PATH)
            push!(LOAD_PATH, this_path)
        end
    end

    # Add the search-path to the extra galleries - temporarily so so needed files can be included
    push!(LOAD_PATH, string(@__DIR__, "/gallery_extra"))
    include("gallery_extra/basic_random_examples.jl")
    include("gallery_extra/gallery_examples.jl")
    include("gallery_extra/distributed_example.jl")
    include("gallery_extra/periodic_dde.jl")
    include("gallery_extra/NLEVP_native.jl")
    include("gallery_extra/bem_hardcoded/bem_hardcoded.jl");
    include("gallery_extra/load_dtn_dimer.jl");

  """
     nep=nep_gallery(name)
     nep=nep_gallery(name,params)
     nep=nep_gallery(name,params;kwargs)

Collection of nonlinear eigenvalue problems.
Returns a NEP object from a gallery of examples of nonlinear eigenvalue problems.
The parameter `name` decides which NEP.

# Supported problems:
The following list describes the NEP with a certain `name` and the associated parameters (`params`) and keyword arguments (`kwargs`), if any.\\

* `dep0`: A random [`DEP`](@ref) with one delay tau = 1, generated with a pseudorandom number generator. \\
    *params:* size (default = 5)

* `dep0_sparse`: A random [`DEP`](@ref) with sparse matrices and one delay tau = 1,  generated with a pseudorandom number generator.\\
    *Params:* Two optional params determining the size (default = 5) and the fill (default = 0.25)

* `dep0_tridiag`: A random [`DEP`](@ref) with sparse tridiagonal matrices and one delay tau = 1, generated with a pseudorandom number generator.\\
    *Params:* size (default = 100)

* `dep_symm_double`: A [`DEP`](@ref) with double eigenvalues and sparse symmetric matrices and one delay tau = 1.\\
    *Params:* size (default = 100)\\
    *Reference:* Example from [H. Voss and M. M. Betcke, Restarting iterative projection methods for Hermitian nonlinear eigenvalue problems with minmax property, Numer. Math., 2017](https://doi.org/10.1007/s00211-016-0804-3)

* `dep_double`: A [`DEP`](@ref) with a double non-semisimple eigenvalue in λ=3πi.\\
    *Reference:* Example from [E. Jarlebring, Convergence factors of Newton methods for nonlinear eigenvalue problems, LAA, 2012](https://doi.org/10.1016/j.laa.2010.08.045)

* `dep1`: A [`DEP`](@ref) with one eigenvalue equal to one.

* `pep0`: A random [`PEP`](@ref), generated with a pseudorandom number generator\\
    *Params:* size (default = 200)

* `pep0_sym`: A random symmetric [`PEP`](@ref), generated with a pseudorandom number generator\\
    *Params:* size (default = 200)

* `pep0_sparse`: A random [`PEP`](@ref) with sparse matrices, generated with a pseudorandom number generator\\
    *Params:* Two optional `params` determining the size (default = 200) and the fill (default = 0.03)

* `real_quadratic`: A quadratic [`PEP`](@ref) with real eigenvalues.\\
    *Solutions:*     Four smallest eigenvalues of the problem:\\
    ∘ `-2051.741417993845`\\
    ∘ ` -182.101627437811`\\
    ∘ `  -39.344930222838`\\
    ∘ `   -4.039879577113`\\

* `dep_distributed`: A NEP corresponding to a distributed delay time-delay system\\
    *Reference:* Example in [E. Jarlebring and W. Michiels and K. Meerbergen, The infinite Arnoldi method and an application to time-delay systems with distributed delays, Delay Systems - Methods, Applications and New Trends, 2012](https://doi.org/10.1007/978-3-642-25221-1_17)\\
    *Solutions:*  Some correct eigenvalues:\\
    ∘ `-0.400236388049641 + 0.970633098237807i`\\
    ∘ `-0.400236388049641 - 0.970633098237807i`\\
    ∘ `2.726146249832675 +  0.000000000000000i`\\
    ∘ `-1.955643591177653 + 3.364550574688863i`\\
    ∘ `-1.955643591177653 - 3.364550574688863i`\\
    ∘ `4.493937056300693 +  0.000000000000000i`\\
    ∘ `-1.631513006819252 + 4.555484848248613i`\\
    ∘ `-1.631513006819252 - 4.555484848248613i`\\
    ∘ `-1.677320660400946 + 7.496870451838560i`\\
    ∘ `-1.677320660400946 - 7.496870451838560i`

* `qdep0`:  A quadratic delay eigenvalue problem.\\
   *Reference:* [S. W. Gaaf and E. Jarlebring, The infinite Bi-Lanczos method for nonlinear eigenvalue problems, SIAM J. Sci. Comput., 2017](https://doi.org/10.1137/16M1084195)

* `qdep1`:  A quadratic delay eigenvalue problem.\\
   *Reference:* [E. Jarlebring and W. Michiels and K. Meerbergen, A linear eigenvalue algorithm for the  nonlinear eigenvalue problem, Numer. Math., 2011](https://doi.org/10.1007/s00211-012-0453-0)

* `qep_fixed_eig`: A quadratic eigenvalue problem with chosen eigenvalues.\\
   *Params:* Two optional `params` determining the size (default = 5) and a vector containing the eigenvalues (default = uniform [-1,1])

* `neuron0`: A [`DEP`](@ref) that stems from the modling of a coupled neuron.\\
   *Reference:* [L. P. Shayer and S. A. Campbell, Stability, bifurcation and multistability in a system of two coupled neurons with multiple time delays, SIAM J. Applied Mathematics, 2000](https://doi.org/10.1137/S0036139998344015). It is also a benchmark example in DDE-BIFTOOL.

* `schrodinger_movebc`: A NEP stemming from the discretization of a Schrödinger equation as described in the NEP-PACK online tutorial.  The nonlinearity contains ``sinh()``, ``cosh()`` and ``sqrt()``.\\
   *Params:* The optional parameters are size of discretization `n`  and domain and potential description `L0`,`L1`,`α` and `V0`.

* `beam`: A [`DEP`](@ref) modelling a beam with delayed stabilizing feedback described. The A1-term is rank one.\\
   *Params:*  size of the matrix (defalut = 100)\\
   *Reference:* [R. Van Beeumen, E. Jarlebring, and W. Michiels,  A rank-exploiting infinite Arnoldi algorithm for nonlinear eigenvalue problems, 2016.](https://doi.org/10.1002/nla.2043)

* `sine`: A NEP formed by the sum of a polynomial and a sine-function. The sine-term has a rank-one matrix coefficient.\\
   *Reference:*  [R. Van Beeumen, E. Jarlebring, and W. Michiels,  A rank-exploiting infinite Arnoldi algorithm for nonlinear eigenvalue problems, 2016.](https://doi.org/10.1002/nla.2043)

* `bem_fichera`:  Represents a boundary element discretization of Helmholtz equation for a domain consisting of the unit cube, except one removed corner (Fichera corner). The mesh is hardcoded.
    *Params:* The parameter `N` determines the size of the problem (default `N = 5`).\\
    *Reference:* The model stems from the model in these papers:  [A boundary element method for solving PDE eigenvalue problems, M. Steinlechner, bachelor thesis, ETH Zürich, 2010](http://sma.epfl.ch/~anchpcommon/students/steinlechner.pdf) and  [Effenberger and Kressner, Chebyshev interpolation for nonlinear eigenvalue problems, BIT Numerical Mathematics, December 2012, Volume 52, Issue 4, pp 933–951](https://doi.org/10.1007/s10543-012-0381-5)

* `dtn_dimer`: NEP with quotients of Bessel functions stemming from the modeling of resonances \\
    *Params:*  This NEP takes two parameters: `data_dir::String` and `l::Int`. The `data_dir` specifies the directory of the dowloaded FEM-matrices (available here [https://umu.app.box.com/s/b52yux3z9rcl8y0l7la22k0vi062cvu5](https://umu.app.box.com/s/b52yux3z9rcl8y0l7la22k0vi062cvu5)). The integer `l` specifies the number of DtN-terms: 2l+1.\\
    *Reference:* [J. Araujo-Cabarcas, C. Engström and E. Jarlebring, Efficient resonance computations for Helmholtz problems based on a Dirichlet-to-Neumann map, J. Comput. Appl. Math., 330:177-192, 2018](http://arxiv.org/pdf/1606.09547))

The MATLAB-package described in
[T. Betcke, N. J. Higham, V. Mehrmann, Ch. Schröder, F. Tisseur, NLEVP: A Collection of Nonlinear Eigenvalue Problems, ACM Transactions on Mathematical Software 39(2), January 2011](https://doi.org/10.1145/2427023.2427024)
provides a number of benchmark problems for NEPs. These are available in NEP-PACK in two different ways.
We have native implementations of some problems (referred to as `nlevp_native_`) and the separate `GalleryNLEVP`.
The native implementation is preferred since the `GalleryNLEVP`  interfaces with MATLAB and is therefore considerably slower.

* `nlevp_native_gun`:  The benchmark problem from the NLEVP-collection called "gun", represented in the native NEP-PACK format.
    [B.-S. Liao, Z. Bai, L.-Q. Lee, and K. Ko. Nonlinear Rayleigh-Ritz iterative method for solving large scale
    nonlinear eigenvalue problems.  Taiwan. Journal of Mathematics, 14(3):869–883, 2010](https://doi.org/10.11650/twjm/1500405872)

* `nlevp_native_cd_player`:  The benchmark problem from the NLEVP-collection called "cd\\_player", represented in the native NEP-PACK format.
    [Y. Chahlaoui, and P. M. Van Dooren, Benchmark examples for model reduction of linear time-
    invariant dynamical systems. In Dimension Reduction of Large-Scale Systems, P. Benner, V. Mehrmann,
    and D. C. Sorensen, Eds. Lecture Notes in Computational Science and Engineering Series, vol. 45.
    Springer-Verlag, Berlin, 380–392, 2005.](https://doi.org/10.1007/3-540-27909-1_24)
    and
    P. M. R. Wortelboer, M. Steinbuch, and  O. H. Bosgra, Closed-loop balanced reduction with
    application to a compact disc mechanism. In Selected Topics in Identification, Modeling and Control.
    Vol. 9. Delft University Press, 47–58, 1996.
    and
    [W. Draijer, M. Steinbuch, and  O. H. Bosgra, Adaptive control of the radial servo system of a
    compact disc player. Automatica 28, 3, 455–462. 1992.](https://doi.org/10.1016/0005-1098(92)90171-B)

* `nlevp_native_fiber`: The benchmark problem from the NLEVP-collection called "fiber", represented in the native NEP-PACK format.
    One of terms in this problem is approximated by interpolation, and may not always coincide with the benchmark.
    [L. Kaufman, Eigenvalue problems in fiber optic design. SIAM J. Matrix Anal. Appl. 28, 1, 105–117, 2006.](https://doi.org/10.1137/S0895479803432708)
    and
    [X. Huang, Z. Bai, and Y. Su, Nonlinear rank-one modification of the symmetric eigenvalue problem. J. Comput. Math. 28, 2, 218–234, 2010](https://doi.org/10.4208/jcm.2009.10-m1002)

* `nlevp_native_hadeler`: The benchmark problem from the NLEVP-collection called "hadeler", represented in the native NEP-PACK format. The problem is of the form ``M(λ)=(e^λ-1)B+A_0+A_2λ^2``. \\
    [Hadeler K.  P.  1967.  Mehrparametrige  und  nichtlineare  Eigenwertaufgaben. Arch.  Rational  Mech. Anal. 27, 4, 306–328.](https://doi.org/10.1007/BF00281717)

* `nlevp_native_pdde_stability`: The benchmark problem from the NLEVP-collection called "pdde_stability", represented in the native NEP-PACK format.
    This problem is a quadratic eigenvalue with arbitrary given size `n`. See
    [E. Jarlebring, The Spectrum of Delay-Differential Equations:
    Numerical Methods, Stability and Perturbation, PhD thesis,
    TU Braunschweig, Institut Computational Mathematics, Germany, 2008](https://nbn-resolving.org/urn:nbn:de:gbv:084-22041) and
    [H. Fassbender, N. Mackey, D. S. Mackey and C. Schroeder, Structured
    Polynomial Eigenproblems Related to Time-Delay Systems, ETNA, 2008, vol 31, pp 306-330](http://etna.mcs.kent.edu/vol.31.2008/pp306-330.dir/)

* `nlevp_native_loaded_string`:  The benchmark problem from the NLEVP-collection called "pdde_stability", represented in the native NEP-PACK format.
    The parameters are (n,kappa,m) where n is the size, and the NEP is a SPMF with rational terms and the coefficient
    matrices are rank one modifications of Toeplitz matrices.\\
    [S. I. Solov"ev. Preconditioned iterative methods for a class of nonlinear eigenvalue problems. Linear Algebra Appl., 415 (2006), pp.210-229.](https://doi.org/10.1016/j.laa.2005.03.034)

# Example
```julia-repl
julia> nep=nep_gallery("dep0",100);
julia> norm(compute_Mlincomb(nep,1.0+1.0im,ones(size(nep,1))))
57.498446538064954
```

# See also the following galleries:
* GalleryNLEVP
* GalleryWaveguide
"""
    nep_gallery(name::String, params...; kwargs...)=nep_gallery(NEP, name, params...; kwargs...)
    function nep_gallery(::Type{T}, name::String, params...; kwargs...) where T<:NEP
        global Gallery
        if !haskey(GALLERY, name)
            error("$name not supported")
        end

        return GALLERY[name](params...; kwargs...)
    end

    const GALLERY = Dict(
        "dep0" => dep0,
        "dep1" => dep1,
        "dep0_sparse" => dep0_sparse,
        "dep0_tridiag" => dep0_tridiag,
        "dep_symm_double" => dep_symm_double,
        "dep_double" => dep_double,
        "pep0" => pep0,
        "pep0_sym" => pep0_sym,
        "pep0_sparse" => pep0_sparse,
        "real_quadratic" => real_quadratic,
        "dep_distributed" => gallery_dep_distributed,
        "qdep0" => qdep0,
        "qdep1" => qdep1,
        "qep_fixed_eig" => qep_fixed_eig,
        "periodicdde" => (params...; kwargs...) -> periodic_dde_gallery(PeriodicDDE_NEP, params...; kwargs...),
        "schrodinger_movebc" => schrodinger_movebc,
        "neuron0" => neuron0,
        "nlevp_native_gun" => nlevp_native_gun,
        "nlevp_native_hadeler" => nlevp_native_hadeler,
        "nlevp_native_cd_player" => nlevp_native_cd_player,
        "nlevp_native_fiber" => nlevp_native_fiber,
        "nlevp_native_pdde_stability" => nlevp_native_pdde_stability,
        "nlevp_native_loaded_string" => nlevp_native_loaded_string,
        "beam" => beam,
        "sine" => sine_nep,
        "bem_fichera" => bem_fichera,
        "dtn_dimer" => load_dtn_dimer
    )

end # end of Gallery-module


text:
using Documenter, GenericSchur, LinearAlgebra

makedocs(
    modules = [GenericSchur],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == true),
    sitename = "GenericSchur.jl",
    pages = ["Overview" => "index.md",
             "Library" => "library.md"
             ]
)

# or maybe just the pkg site?
deploydocs(repo = "github.com/RalphAS/GenericSchur.jl.git")


text:
module XLNet

using Flux
using Requires
using Requires: @init
using BSON
using Transformers

# Write your package code here.

end


text:
include("handler/index_page.jl")
include("handler/state.jl")
include("handler/misc.jl")
include("handler/callback_context.jl")
include("handler/make_handler.jl")


text:
# using CSV
# using Dates
# using Plots
# using JLD

## Main simulation function

function simulateODEs(model_def, simul_def)

    # First make a folder where to save the results
    cudi = pwd(); #@__DIR__;
    if !isdir(string(cudi,"\\Results"))
        mkdir(string(cudi,"\\Results"))
    end
    if !isdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
        mkdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
    end

    # Check which of the 2 types of simul_def have been given and make a choice on how to proceed
    if haskey(simul_def, "switchT")
        simul_def = checkStructSimul(model_def, simul_def)
    else
        simul_def = extractSimulCSV(model_def, simul_def)
    end

    # Make empty structure where the simulations will be
    simuls = Dict();

    # Get ODEs function
    ODEsFun = Symbol(string(join(model_def["NameF"]),"_SolveAll"));

    # loop over the different experiments to simulate
    for i in 1:simul_def["Nexp"]
        # Define simulation parameters
        global ts = collect(0.0:round(simul_def["finalTime"][i]))'; # 24 hour experiment
        global sp = [convert(Int,v) for v in (round.(simul_def["switchT"][i])')]; # 8 step experiment for now
        global ivss = simul_def["y0"][i]; # Experimental values for Y0 (IPTG, aTc, RFP, GFP)

        if model_def["Y0eqs"] != [] # ON inputs
            global pre = simul_def["preInd"][i];
        else
            global pre = [];
        end
        global samps = convert.(Int, round.(simul_def["tsamps"][i]));
        global theta = simul_def["theta"];

        # Inputs case
        if model_def["nInp"] != 0
            r = convert.(Int, 1:model_def["nInp"]:(length(simul_def["uInd"][i])));
            global inputs = zeros(convert.(Int,length(simul_def["uInd"][i])));
            for j in 1:convert.(Int,length(simul_def["uInd"][i])/model_def["nInp"])
                for k in 0:(model_def["nInp"]-1)
                    inputs[r[j]+k] = simul_def["uInd"][i][j,(k+1)];
                end
            end
        else
            r=0
        inputs=[];
        end

        simul = @eval $ODEsFun(ts, theta, sp, inputs, ivss, samps, pre);

        simuls[string("Exp_", i)] = simul;
    end

    simul_def["savepath"] = string(cudi, "\\Results\\", model_def["NameF"],"_",today());
    simul_def["savename"] = string(model_def["NameF"],"_",today(), "_SimulationResults_",simul_def["flag"],".jld");

    if simul_def["plot"] == true
        plotSimsODE(simuls,model_def,simul_def)
        println("")
        println("----------------------------------------- PLOTS -----------------------------------------")
        println("Simulation PLOTS are saved in the directory: ")
        println(string("                 ", simul_def["savepath"]))
        println(string("Under the name PlotSimulation_Exp(i)_",simul_def["flag"],".png"))
        println("--------------------------------------------------------------------------------------")
        println("")
    end

    save(string(simul_def["savepath"], "\\", simul_def["savename"]), "simuls", simuls, "model_def", model_def, "simul_def", simul_def)

    println("")
    println("----------------------------------------- RESULTS -----------------------------------------")
    println("Simulation results are saved in the directory: ")
    println(string("                 ", simul_def["savepath"]))
    println(string("Under the name ",simul_def["savename"]))
    println("--------------------------------------------------------------------------------------")
    println("")

    return simuls, model_def, simul_def

end

## Simulation Data structure Dict
function defSimulStruct()

    # The structure of things user should give
    simul_def = Dict()

    simul_def["Nexp"] = []; # Integer indicating the number of experiments to be simulated
    simul_def["finalTime"] = []; # -> Vector of final times for each simulation (initial time will allways be asumed as 0, so please consider that)
    simul_def["switchT"] = []; # Array with the switching times of the inducer in the simulation (time 0 and final time need to be considered)
    simul_def["y0"] = []; # Array of Y0s for the simulations for each experiment. If you are computing the steady state this vector might not be used, however you still need to introduce it
    simul_def["preInd"] = []; # Level of the inducer in the ON. This entry can be empty if no Steady State equations are given and only the Y0 given by the user is considered.
    simul_def["uInd"] = []; # Array with the values for the inducer for each experiment at each step
    simul_def["theta"] = []; # Vector/Matrix with the parameter samples or directory and file location of CSV file with them
    simul_def["tsamps"] = []; # Array of Sampling times vectors
    simul_def["plot"] = []; # Bollean or yes/no string to save the resulting simulations in the results directory (false will be considered as default)
    simul_def["flag"] = []; # String to attach a unique flag to the generated files so it is not overwritten. If empty, nothing will be added.

    return(simul_def)

end

# Check Simulations function
# Function to check the content of the fields the user has given
function checkStructSimul(model_def, simul_def)

    # Check taht all the dictionary entries are correct
    entries = ["Nexp", "finalTime", "switchT", "y0", "preInd", "uInd", "theta", "tsamps", "plot", "flag"]
    if symdiff(entries,keys(simul_def))!=[] && symdiff(entries,keys(simul_def)) != ["savepath", "savename"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the entries of the dictionary, there is something wrong...")
        println(symdiff(entries,keys(simul_def)))
        return
    end

    # Check one to be sure there is no empty entry
    nee = ["Nexp", "finalTime", "switchT", "y0", "uInd", "theta", "tsamps"]; # No empty entries
    for i in 1:length(nee)
        if simul_def[nee[i]] == []
            println("-------------------------- Process STOPPED!!! --------------------------")
            println(string("Plese, check the field ", nee[i], ", this cannot be empty"))
            return
        end
    end

    # Checks to see if the content of the fields is the expected (also consier if single value entries have been given as a number or a vector)
    if ((typeof(simul_def["Nexp"])!=Array{Int,1}) || length(simul_def["Nexp"])!=1) && (typeof(simul_def["Nexp"])!=Int)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field Nexp! This should be an integer. ")
        return
    elseif (typeof(simul_def["finalTime"]) != Array{Int,1}) && (typeof(simul_def["finalTime"]) != Array{Float64,1}) && (typeof(simul_def["finalTime"]) != Array{Float32,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field finalTime! This should be a vector of final times. ")
        return
    elseif (typeof(simul_def["switchT"]) != Array{Array{Int,1},1}) && (typeof(simul_def["switchT"]) != Array{Array{Float64,1},1}) && (typeof(simul_def["switchT"]) != Array{Array{Float32,1},1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field switchT! This should be an array of arrays of numbers. ")
        return
    elseif (typeof(simul_def["y0"]) != Array{Array{Int,1},1}) && (typeof(simul_def["y0"]) != Array{Array{Float64,1},1}) && (typeof(simul_def["y0"]) != Array{Array{Float32,1},1}) &&
        (typeof(simul_def["y0"]) != Array{Array{Int,N} where N,1}) && (typeof(simul_def["y0"]) != Array{Array{Float64,N} where N,1}) && (typeof(simul_def["y0"]) != Array{Array{Float32,N} where N,1}) &&
        (typeof(simul_def["y0"]) != Array{Array{Int,2},1}) && (typeof(simul_def["y0"]) != Array{Array{Float64,2},1}) && (typeof(simul_def["y0"]) != Array{Array{Float32,2},1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field y0! This should be an array of arrays of numbers. ")
        return
    elseif (typeof(simul_def["preInd"]) != Array{Array{Int,1},1}) && (typeof(simul_def["preInd"]) != Array{Array{Float64,1},1}) && (typeof(simul_def["preInd"]) != Array{Array{Float32,1},1}) && simul_def["preInd"] != []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field preInd! This should be an array of arrays of numbers. ")
        return
    # elseif (typeof(simul_def["uInd"]) != Array{Array{Int,1},1}) && (typeof(simul_def["uInd"]) != Array{Array{Float64,1},1}) && (typeof(simul_def["uInd"]) != Array{Array{Float32,1},1})
    #     println("-------------------------- Process STOPPED!!! --------------------------")
    #     println("Please, check the field uInd! This should be an array of arrays of numbers. ")
    #     return
    elseif (typeof(simul_def["theta"]) != Array{Float64,1}) && (typeof(simul_def["theta"]) != Array{Float32,1}) &&
        (typeof(simul_def["theta"]) != Array{Int,1}) && (typeof(simul_def["theta"]) != Array{Int,1}) &&
        (typeof(simul_def["theta"]) != Array{Float64,2}) && (typeof(simul_def["theta"]) != Array{Float32,2}) &&
        (typeof(simul_def["theta"]) != Array{String,1}) && (typeof(simul_def["theta"]) != String)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field theta! This should be a vector, matrix or string. ")
        return
    elseif (typeof(simul_def["tsamps"]) != Array{Array{Int,1},1}) && (typeof(simul_def["tsamps"]) != Array{Array{Float64,1},1}) && (typeof(simul_def["tsamps"]) != Array{Array{Float32,1},1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field tsamps! This should be an array of arrays of numbers. ")
        return
    elseif (simul_def["plot"] != [])
        try
            if (typeof(simul_def["plot"]) != Bool) && (typeof(simul_def["plot"]) != Array{Bool,1}) &&
            (typeof(simul_def["plot"]) != String) && (typeof(simul_def["plot"]) != Array{String,1})
                println("-------------------------- Process STOPPED!!! --------------------------")
                println("Please, check the field plot! This should be true, fals, \"Yes\", \"yes\", \"No\", \"no\" or an empty array ")
                return
            end
        catch
             if (typeof(simul_def["plot"]) != Array{Bool,1}) &&
            (typeof(simul_def["plot"]) != Array{String,1})
                println("-------------------------- Process STOPPED!!! --------------------------")
                println("Please, check the field plot! This should be true, fals, \"Yes\", \"yes\", \"No\", \"no\" or an empty array ")
                return
            end
        end
    elseif (typeof(simul_def["flag"]) != Array{String,1}) && (typeof(simul_def["flag"]) != String) && ((simul_def["flag"]) != [])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field flag! This should be a vector, matrix or string. ")
    return

    end

    # Extract necessary elements to ease generalisation
    if typeof(simul_def["Nexp"]) == Array{Int,1}
        simul_def["Nexp"] = simul_def["Nexp"][1];
    end

    if simul_def["plot"]==[true] || simul_def["plot"]==["Yes"] || simul_def["plot"]==["yes"] || simul_def["plot"]=="Yes" || simul_def["plot"]=="yes"
        simul_def["plot"]=true
    elseif simul_def["plot"]==[false] || simul_def["plot"]==["No"] || simul_def["plot"]==["no"] || simul_def["plot"]=="No" || simul_def["plot"]=="no" || simul_def["plot"]==[]
        simul_def["plot"]=false
    end

    if (typeof(simul_def["theta"]) == Array{String,1})
        simul_def["theta"] = simul_def["theta"][1];
    end

    if typeof(simul_def["flag"]) == Array{String,1}
        simul_def["flag"] = simul_def["flag"][1];
    elseif typeof(simul_def["flag"]) == String
        simul_def["flag"] = simul_def["flag"];
    elseif ((simul_def["flag"]) == [])
        simul_def["flag"] = "";
    end

    # Check that all the contents make sense
    if (simul_def["Nexp"] != length(simul_def["finalTime"])) || (simul_def["Nexp"] != length(simul_def["switchT"])) ||
        (simul_def["Nexp"] != length(simul_def["y0"])) || (simul_def["Nexp"] != length(simul_def["uInd"])) ||
        (simul_def["Nexp"] != length(simul_def["tsamps"]))
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, some field does not match the number of experiments selected. Check that the number of")
        println("entries for finalTime, switchT, y0, uInd or tsamps matches the number of experiments in Nexp.")
        return
    end

    if (model_def["Y0eqs"] != []) && simul_def["preInd"]==[]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("You specified computation of Y0 as steady-state but you have not introduced an inducer value for it!")
        return
    end

    if (typeof(simul_def["theta"]) == Array{Float64,1}) || (typeof(simul_def["theta"]) == Array{Float32,1})
        if length(simul_def["theta"]) != model_def["nPar"]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Number of parameters introduced does not match the specified")
            return
        end
    elseif (typeof(simul_def["theta"]) == Array{Float64,2}) || (typeof(simul_def["theta"]) == Array{Float32,2})
        if size(simul_def["theta"])[1] != model_def["nPar"] && size(simul_def["theta"])[2] != model_def["nPar"]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Number of parameters introduced does not match the specified")
            return
        end
    elseif (typeof(simul_def["theta"]) == String)
        if simul_def["theta"][end-3:end] != ".csv"
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, add the .csv termination to the theta file! And no spaces after!")
            return
        end
        if isfile(simul_def["theta"])
            simul_def["theta"] = Matrix(CSV.read(simul_def["theta"]));
        else
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Sorry, but the file path you introduced does not exist!")
            return
        end
    end

    for i in 1:simul_def["Nexp"]
        if simul_def["tsamps"][i][end] > simul_def["finalTime"][i]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, check finalTime. You have selected a sampling point past it.")
            return
        end
        if model_def["nInp"] != 0
            if length(simul_def["uInd"][i]) != (length(simul_def["switchT"][i])-1)
                println("-------------------------- Process STOPPED!!! --------------------------")
                println("Please, check uInd and switchT. Number of steps does not match the number of values for the inputs.")
                return
            end
        end

        if length(simul_def["y0"][i])/model_def["nStat"] == 1
            if length(simul_def["y0"][i]) != model_def["nStat"]
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but Y0 does not match the number or states in experiment ", i))
                return
            end
        else
            if size(simul_def["y0"][i])[1] != model_def["nStat"] && size(simul_def["y0"][i])[2] != model_def["nStat"]
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but Y0 does not match the number or states in experiment ", i))
                return
            end
            if length(simul_def["y0"][i])/model_def["nStat"] != length(simul_def["theta"])/model_def["nPar"]
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but Y0 does not match the dimensions of theta "))
                return
            end
            if size(simul_def["y0"][i])[1] == size(simul_def["y0"][i])[2]
                println("-------------------------- WARNING --------------------------")
                println(string("Sorry, but the number of rows and columns of the y0 matrix in experiment ",i," is the same, so the checks on "))
                println("correct orientation will not work. Please make sure that the dimensions follow: ")
                println("y0[samples, states]")
            end
        end
    end

    # Check warning in case the matrix introduced is symmetric
    if length(size(simul_def["theta"])) == 2
        if size(simul_def["theta"])[1] == size(simul_def["theta"])[2]
            println("-------------------------- WARNING --------------------------")
            println("Sorry, but the number of rows and columns of the theta matrix is the same, so the checks on ")
            println("correct orientation will not work. Please make sure that the dimensions follow: ")
            println("theta[samples, parameters]")
        end
    end

    # Check that no step is no smaller than 2 unit of time
    for i in 1:simul_def["Nexp"]
        for j in 1:length(simul_def["uInd"][i])
            if (simul_def["switchT"][i][j+1]-simul_def["switchT"][i][j])<=4
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but 2 of the steps in experiment ", i, " are too close. This package cannot "))
                println("handle this for now. We are working on it!")
                return
            end
        end
    end

    # Check that the inputs for the experiment are given as collumns
    try
        if model_def["nInp"]>1
            for i in 1:simul_def["Nexp"]
                if size(simul_def["uInd"][i])[2] != model_def["nInp"];
                    simul_def["uInd"][i] = simul_def["uInd"][i]';
                end
            end
        end
    catch
        println("-------------------------- Process STOPPED!!! --------------------------")
        println(string("Sorry, but there is some issue with the contents of the field uInd."))
        return
    end

    return(simul_def)
end

## CSV introduction description
function fileStructInfo()

    println("ALL THE FILES MUST BE CSV FILES!!!")
    println("")
    println("---------------------------------------------------------------------------------------------------------")
    println("")
    println("The observables file entry should have the following structure: ")
    println("Column 1: Sampling Times for the simulations")
    println("Column 2 + Number of states: Y0 value for each state (in order) in each column. ")
    println("           Value might need to be repeated across all the column, but only the first row will be considered")
    println("")
    println("---------------------------------------------------------------------------------------------------------")
    println("")
    println("The event inputs file entry should have the following structure: ")
    println("Column 1: Each Switching time (final time not considered)")
    println("Column 2: Final time (can be repeated for as many entrances as column 1 has)")
    println("Column 3 + number of inducers: Repeated column for the value of the inducers in the ON. ")
    println("              If no pre-inducer will be used in the simulations, any number could be used")
    println("Column 4 + number of inducers: Value of each inducer at each switching time.")
    println("")
    println("This is an Example of an experiment with 3 steps and 2 inducers: ")
    println(" |Switching|  |FinalTime|  |IPTGpre|  |aTcPre|  |  IPTG  |  |  aTc  |")
    println("      0           1350         1          0          1          0    ")
    println("     400          1350         1          0         0.3         50   ")
    println("     700          1350         1          0         0.7         35   ")
    println("     950          1350         1          0         0.1        100   ")
    println("")
    println("---------------------------------------------------------------------------------------------------------")
    println("")
    println("For the main directory, remember that in Julia you have to type \\ twice!")

end

## CSV introduction structure
function defSimulStructFiles()

    fileStructInfo()

    # The structure of things user should give
    simul_def = Dict()

    simul_def["ObservablesFile"] = []; # File containing information about the states Y0 and sampling times
    simul_def["EventInputsFile"] = []; # -> File containing all the information about the stimuli
    simul_def["theta"] = []; # Vector/Matrix with the parameter samples or directory and file location of CSV file with them
    simul_def["MainDir"] = []; # Directory path where all files are located (to avoid repetition). This can be empty.
    simul_def["plot"] = [];
    simul_def["flag"] = [];

    return(simul_def)

end

## Extract data from CSV
# Function to extract contents of CSV files if user sellects that
function extractSimulCSV(model_def, simul_def)

    # Check taht all the dictionary entries are correct
    entries = ["ObservablesFile", "EventInputsFile", "theta", "MainDir", "plot", "flag"]
    if symdiff(entries,keys(simul_def))!=[]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the entries of the dictionary, there is something wrong...")
        println(symdiff(entries,keys(simul_def)))
        return
    end

    # Checks on the different strings introduced
    if simul_def["MainDir"] == [] || simul_def["MainDir"] == ""
        simul_def["MainDir"] = "";
    elseif typeof(simul_def["MainDir"])==Array{String,1}
        simul_def["MainDir"] = simul_def["MainDir"][1]
    elseif typeof(simul_def["MainDir"]) == String
        nothing
    else
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there is something wrong in the type of the field MainDir")
        return
    end

    if typeof(simul_def["ObservablesFile"]) != Array{String,1}
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there is something wrong in the type of the field ObservablesFile. This should be a Array{String,1}")
        return
    end

    if length(simul_def["ObservablesFile"]) != length(simul_def["EventInputsFile"])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there number of entries for ObservablesFile and EventInputsFile should be the same")
        return
    end

    for i in 1:length(simul_def["ObservablesFile"])
        if simul_def["ObservablesFile"][i][end-3:end] != ".csv"
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, add the .csv termination to the ObservablesFile file! And no spaces after!")
            return
        end
    end

    if typeof(simul_def["EventInputsFile"]) != Array{String,1}
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there is something wrong in the type of the field EventInputsFile. This should be a Array{String,1}")
        return
    end

    for i in 1:length(simul_def["EventInputsFile"])
        if simul_def["EventInputsFile"][i][end-3:end] != ".csv"
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, add the .csv termination to the EventInputsFile file! And no spaces after!")
            return
        end
    end

    simstu = defSimulStruct();

    simstu["theta"] = simul_def["theta"];

    tmpInp = Array{Any}(undef, 1,length(simul_def["EventInputsFile"]))
    for i in 1:length(simul_def["EventInputsFile"])
        if simul_def["MainDir"] != [] && simul_def["MainDir"] != ""
            tmpInp[i] = Matrix(CSV.read(string(simul_def["MainDir"], "\\\\",simul_def["EventInputsFile"][i])));
        else
            tmpInp[i] = Matrix(CSV.read(simul_def["EventInputsFile"][i]));
        end
    end

    tmpObs = Array{Any}(undef, 1,length(simul_def["ObservablesFile"]))
    for i in 1:length(simul_def["ObservablesFile"])
        if simul_def["MainDir"] != [] && simul_def["MainDir"] != ""
            tmpObs[i] = Matrix(CSV.read(string(simul_def["MainDir"], "\\\\",simul_def["ObservablesFile"][i])));
        else
            tmpObs[i] = Matrix(CSV.read(simul_def["ObservablesFile"][i]));
        end
    end

    try
        simstu["Nexp"] = length(simul_def["EventInputsFile"]);
        simstu["finalTime"] = [tmpInp[i][1,2] for i in 1:length(simul_def["EventInputsFile"])];
        simstu["switchT"] = [vcat(tmpInp[i][:,1], tmpInp[i][1,2]) for i in 1:length(simul_def["EventInputsFile"])];
        simstu["y0"] = [tmpObs[i][1,2:(1+model_def["nStat"])] for i in 1:length(simul_def["ObservablesFile"])];
        simstu["preInd"] = [tmpInp[i][1,3:(2+model_def["nInp"])] for i in 1:length(simul_def["EventInputsFile"])];
        simstu["uInd"] = [tmpInp[i][:,(3+model_def["nInp"]):(3+model_def["nInp"]+(model_def["nInp"]-1))] for i in 1:length(simul_def["EventInputsFile"])];
        simstu["theta"] = simul_def["theta"];
        simstu["tsamps"] = [tmpObs[i][:,1] for i in 1:length(simul_def["ObservablesFile"])];
        simstu["plot"] = simul_def["plot"];
        simstu["flag"] = simul_def["flag"];
    catch
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but there seems to be some issue with the internal structure of the files you have given...")
        return
    end

    simstu = checkStructSimul(model_def, simstu);

    return(simstu)

end

## Plot simulations function
function plotSimsODE(simuls,model_def,simul_def)
    cudi = pwd(); # @__DIR__;
    for i in 1:simul_def["Nexp"]
        # Title for each subplot
        tit = "";
        yl1 = "";
        for k in 1:model_def["nStat"]
            tit = hcat(tit, string(model_def["stName"][k]));
            yl1 = hcat(yl1, "y");
        end
        tit = tit[:,2:end];
        yl1 = yl1[:,2:end];

        titu = "";
        yl2 = "";
        if model_def["nInp"] != 0
            for k in 1:model_def["nInp"]
                titu = hcat(titu, string(model_def["inpName"][k]));
                yl2 = hcat(yl2, "u");
            end
            titu = titu[:,2:end];
            yl2 = yl2[:,2:end];
        end

        tuu = hcat(tit, titu);
        yuu = hcat(yl1, yl2);

        # Elements for the inducers
        tu = Array{Array{Any,1}}(undef, model_def["nStat"]+model_def["nInp"])
        [tu[k] = [] for k in 1:model_def["nStat"]]
        su = Array{Array{Any,1}}(undef, model_def["nStat"]+model_def["nInp"])
        [su[k] = [] for k in 1:model_def["nStat"]]
        for k in 1:model_def["nInp"]
            tu[k+model_def["nStat"]] = round.(simul_def["switchT"][i]);
            su[k+model_def["nStat"]] = vcat(simul_def["uInd"][i][:,(k)], simul_def["uInd"][i][end,(k)])
        end

        t = round.(simul_def["tsamps"][i]);
        s = simuls[string("Exp_",i)];

        pl = plot(t, s[:,:,1], layout=model_def["nStat"]+model_def["nInp"], label = "", title = tit, xlabel = "time", ylabel = yuu,
        color = "blue", size = [2000, 1200])

        if model_def["nStat"]>1
            for p in 2:size(simuls[string("Exp_",i)])[3]
                pl = plot!(t,s[:,:,p], layout=model_def["nStat"], label = "", color = "blue");
            end
        end

        pl = plot!(tu, su, layout=model_def["nStat"]+model_def["nInp"], label = "", linetype=:step, title = tuu);

        savefig(string(cudi, "\\Results\\", model_def["NameF"],"_",today(), "\\PlotSimulation_Exp", i,"_", simul_def["flag"], ".png"));

    end

end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


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

using Markdown
using InteractiveUtils

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

# ╔═╡ da1d65a0-ec42-11ea-0141-334c9eeeb035
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add(["PlutoUI", "Images", "ImageMagick"])
	using PlutoUI
	using Images
end

# ╔═╡ 5ef51c3a-70a7-11eb-2023-31113399a57f
html"""

Section 1.2

 Intro to Abstractions

Lecture Video

"""

# ╔═╡ 60ae819a-70a7-11eb-31d4-750c7f5dc6ca
PlutoUI.TableOfContents(aside=true)

# ╔═╡ 792c6a62-ec41-11ea-01f3-73e7eee23cc7
md"""
#### Intializing packages

_When running this notebook for the first time, this could take up to 15 minutes. Hang in there!_
"""

# ╔═╡ ef1bfa16-70ea-11eb-189c-a54db292cd6f
md"
## Introduction

The goal of this section is to introduce you to the notion of abstraction. You can think of abstraction as an opposite to specialization. We will illustrate this by looking at the following example.

### What is _one_?

Before we get lost talking about the foundations of number theory, I will present you with a few examples that represent one to me.
"

# ╔═╡ 6fcac482-70ee-11eb-0b80-ff41c708053b
md"Each of the items in this list is a specific, or **_specialized_** representation of _one_:
1. as an integer
1. as a float
1. as a string
1. as a rational number
1. as a cute picture
1. as an 2x2 identity matrix
1. as a singular dog

Of course, these are just a few examples of _one_. People have been representing _one_ for ages in different langauges, scripts, artistic expression, etc.

The difference between these ones to me is clear. In fact, I just articulated it to you. Now, let's turn to how a computer sees _one_ differently based on what I type.
"

# ╔═╡ 9ebc079a-70f0-11eb-07d9-f9e80f3f4584
md"So to a computer, all of these are different types."

# ╔═╡ 15f7f90a-70f0-11eb-0d41-63677e4023f4
md"### What is a collection of _one_s?

Now, I want to make a collection of ones for some reason. Below is a way for you to experiment building this collection with different _one_s. As you do this experiment, I want you to look at what stays in the same in the Julia output, and what doesn't."

# ╔═╡ f6886d90-70ed-11eb-07c4-471ee267e7c1
md"""
Before we even look at the output, I am amazed that this code even ran. Are you telling me that the computer doesn't care which _one_ I am using in my array?

Yes! That's exactly what abstraction is. By stepping back, we can now think and operate at a level that doesn't care about which _specific_ one I am using. This is what we mean by **abstraction is the opposite of specialization**.

The information that Julia gives back is quite informative. Here is an example of the first line of the output for a few different types:

```
array = 3x4 Array{Int64, 2}
array = 3x4 Array{Float64, 2}
array = 3x4 Array{Rational{Int64}, 2}
```

Notice that for all of these, we have the same `3x4 Array{***, 2}`.
"""

# ╔═╡ 3c1a3cf8-70f8-11eb-3c18-375207f321eb
md"""
## First Taste of Abstraction

Now, I want to do something to a collection of ones, that doesn't care about which one I'm using. So I'm going to write a function that takes in my collection, and add a corgi whereever I desire.
"""

# ╔═╡ 19f4ddb0-ec44-11ea-20b9-5d97fb2b1cf4
function insert(new, A, i, j)
	B = copy(A)
	B[i,j] = new
	return B
end

# ╔═╡ 424f5f10-ec44-11ea-076d-f3cba4435e0c
begin
	A = fill(0, 3, 4)
	md"""
	$(@bind i Slider(1: size(A,1), show_value=true))
	$(@bind j Slider(1: size(A,2), show_value=true))
	"""
end

# ╔═╡ 71ac08ea-7145-11eb-237d-5506adfb9533
begin
	one_number_array = fill(1,3,4)
	insert(5, one_number_array, i, j)
end

# ╔═╡ ee43d808-70fa-11eb-0cc6-337279f41494
md"This is still amazing. I wrote one function that just cares about how to insert an object into an array, without knowing anything about what's inside, and it worked for two completely different arrays, _collections of ones of **any kind**_."

# ╔═╡ 263a8a0a-70ee-11eb-236d-c941ba63dff3
md"
## Conclusion
The key idea here is that a computer language should allow you to do operations that make sense. Often times, an operation can make sense for many different objects. So we can abstract away the specifics of the object in our implementation. It should let you step back from there.
"

# ╔═╡ 52461588-ea1a-4e7d-aec2-3de388d31656
md"""
## Appendix
"""

# ╔═╡ 1a2a9000-ec43-11ea-3f39-8312ea286a92
begin
	oneimage = load(download("https://gallery.yopriceville.com/var/albums/Free-Clipart-Pictures/Decorative-Numbers/Cute_Number_One_PNG_Clipart_Image.png?m=1437447301"))
	corgi = load(download("https://i.barkpost.com/wp-content/uploads/2015/01/corgi2.jpg?q=70&fit=crop&crop=entropy&w=808&h=500"))
	nothing
end

# ╔═╡ 0504ac94-70ee-11eb-1c4e-977d9e7d35c9
one = [
	1,
	1.0,
	"one",
	1//1,
	oneimage,
	[1 0; 0 1],
	corgi,
]

# ╔═╡ 0b1668ba-ec42-11ea-3e50-ed97c5b17ced
computer_ones = typeof.(one)

# ╔═╡ b2239b96-70ef-11eb-0b85-21ecab25dc9f
begin
	one_keys = ["1", "1.0", "one", "1//1", "Cute One", "2x2 Identity", "One Corgi"]
	selections = one_keys .=> one_keys
	lookup_element = Dict(one_keys .=> one)
	md"$(@bind element_key Select(selections))"
end

# ╔═╡ 4251f668-70aa-11eb-3d89-35f8d53b7d9b
# your chosen one
element = lookup_element[element_key]

# ╔═╡ f1568d10-ec41-11ea-3dd2-a9cb273ce5b8
#its type
typeof(element)

# ╔═╡ ab02d850-ec41-11ea-10b2-a1b600b12658
# a 3x4 array of this one.
array = fill(element,3,4)

# ╔═╡ 5363a400-ec44-11ea-284e-d13a8872551c
begin
	one_image_array = fill(oneimage,3,4)
	insert(corgi, one_image_array, i, j)
end

# ╔═╡ Cell order:
# ╟─5ef51c3a-70a7-11eb-2023-31113399a57f
# ╟─60ae819a-70a7-11eb-31d4-750c7f5dc6ca
# ╟─792c6a62-ec41-11ea-01f3-73e7eee23cc7
# ╟─ef1bfa16-70ea-11eb-189c-a54db292cd6f
# ╠═0504ac94-70ee-11eb-1c4e-977d9e7d35c9
# ╟─6fcac482-70ee-11eb-0b80-ff41c708053b
# ╠═0b1668ba-ec42-11ea-3e50-ed97c5b17ced
# ╟─9ebc079a-70f0-11eb-07d9-f9e80f3f4584
# ╟─15f7f90a-70f0-11eb-0d41-63677e4023f4
# ╠═b2239b96-70ef-11eb-0b85-21ecab25dc9f
# ╠═4251f668-70aa-11eb-3d89-35f8d53b7d9b
# ╠═f1568d10-ec41-11ea-3dd2-a9cb273ce5b8
# ╠═ab02d850-ec41-11ea-10b2-a1b600b12658
# ╟─f6886d90-70ed-11eb-07c4-471ee267e7c1
# ╟─3c1a3cf8-70f8-11eb-3c18-375207f321eb
# ╠═19f4ddb0-ec44-11ea-20b9-5d97fb2b1cf4
# ╠═424f5f10-ec44-11ea-076d-f3cba4435e0c
# ╠═5363a400-ec44-11ea-284e-d13a8872551c
# ╠═71ac08ea-7145-11eb-237d-5506adfb9533
# ╟─ee43d808-70fa-11eb-0cc6-337279f41494
# ╟─263a8a0a-70ee-11eb-236d-c941ba63dff3
# ╟─52461588-ea1a-4e7d-aec2-3de388d31656
# ╠═da1d65a0-ec42-11ea-0141-334c9eeeb035
# ╠═1a2a9000-ec43-11ea-3f39-8312ea286a92


text:
// wrong type of main

void main(int argc)
{
   return;
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
#=
This file contains a few algorithms that are used to compute convex
hulls in two dimensions.

Currently included are:

* Graham Scan Algorithm
* Monotone Chain Algorithm

Algorithms that should be added are:

* Quickhull Algorithm
* Chan Algorithm
* ?
=#
module CHull2d

# Load dependencies
using FixedSizeArrays

# Base methods
import Base: *, in, intersect, isfinite, isless, isequal, length

# Include file of all types
include("CHullTypes.jl")
include("utils.jl")

#
# Include all algorithms
#
include("intersect.jl")
include("algorithms.jl")

const ALGDICT = Dict{Symbol, Function}()
ALGDICT[:MonotoneChain] = _monotonechain
ALGDICT[:GrahamScan] = _grahamscan
# ALGDICT[:QuickHull] = _quickhull

"""
Creates the convex hull of a set of points

* points : The set of points that you would like to take the convex hull of
* algorithm : Which convex hull algorithm to use. Currently implemented are
              [:GrahamScan, :MonotoneChain]
* _at : Whether or not to use the Akl Toussaint pruning
"""
function convexhull{T}(points::Vector{Point{2, T}}; algorithm::Symbol=:MonotoneChain, _at::Bool=true)
    # TODO: Should make sure that there are at least 3 not co-linear points
    # First prune points
    points = _at ? _akltoussaint(points) : points

    # Apply algorithm
    ep = ALGDICT[algorithm](points)

    return ConvexHull(ep)
end

export ConvexHull, LineSegment, Quadrant, Point, # Types
       convexhull, wrappoints, pt_out_tri,
       evaluatex, evaluatey, ccw, cw, orientation

end


text:
import ..TransferOperators:
	AbstractTransferOperator

import LinearAlgebra:
	norm

"""
    invariantmeasure(to::AbstractTransferOperator; N::Int = 200)

Compute the invariant probability distribution from a transfer operator.

## Computing an invariant probability distribution
The distribution is taken as a left eigenvector of the transfer matrix,
obtained by repeated application of the transfer operator on a randomly
initialised distribution until the probability distribution converges.

## Keyword arguments
- `N`: the maximum number of iterations.

## Examples

```julia
# Create some random points and make sure they are invariant
pts = [rand(3) for i = 1:30]
invariant_pts = invariantize(pts)

# Compute the transfer operator for a triangulation of the points
to = transferoperator_triangulation_approx(invariant_pts)

# Compute the invariant measure over the simplices of the triangulation
ivm = invariantmeasure(to)
```
"""
function invariantmeasure(to::AbstractTransferOperator;
			N::Int = 200,
            tolerance::Float64 = 1e-8,
            delta::Float64 = 1e-8)
    #=
    # Start with a random distribution `Ρ` (big rho). Normalise it so that it
    # sums to 1 and forms a true probability distribution over the simplices.
    =#
    Ρ = rand(Float64, 1, size(to.transfermatrix, 1))
    Ρ = Ρ ./ sum(Ρ, dims = 2)

    #=
    # Start estimating the invariant distribution. We could either do this by
    # finding the left-eigenvector of M, or by repeated application of M on Ρ
    # until the distribution converges. Here, we use the latter approach,
    # meaning that we iterate until Ρ doesn't change substantially between
    # iterations.
    =#
    distribution = Ρ * to.transfermatrix

    distance = norm(distribution - Ρ) / norm(Ρ)

    check = floor(Int, 1 / delta)
    check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N))
    check_pts = check_pts[check_pts .> 0]
    num_checkpts = size(check_pts, 1)
    check_pts_counter = 1

    counter = 1
    while counter <= N && distance >= tolerance
        counter += 1
        Ρ = distribution

        # Apply the Markov matrix to the current state of the distribution
        distribution = Ρ * to.transfermatrix

        if (check_pts_counter <= num_checkpts &&
           counter == check_pts[check_pts_counter])

            check_pts_counter += 1
            colsum_distribution = sum(distribution, dims = 2)[1]
            if abs(colsum_distribution - 1) > delta
                distribution = distribution ./ colsum_distribution
            end
        end

        distance = norm(distribution - Ρ) / norm(Ρ)
    end
    distribution = dropdims(distribution, dims = 1)

    # Do the last normalisation and check
    colsum_distribution = sum(distribution)

    if abs(colsum_distribution - 1) > delta
        distribution = distribution ./ colsum_distribution
    end

    # Find partition elements with strictly positive measure.
    δ = tolerance/size(to.transfermatrix, 1)
    simplex_inds_nonzero = findall(distribution .> δ)

    # Extract the elements of the invariant measure corresponding to these indices
    return InvariantDistribution(distribution,simplex_inds_nonzero)
end

export invariantmeasure


text:
# debug why discontinuous at border.

include("../examples/helpers/visualization.jl")

PyPlot.close("all")

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

Random.seed!(25)

fig_num = 1

radius = 0.3 #3.0
δ = 1e-5

#xq = [-2.58, 9.13]
#xq =  [-2.038690328209032, 9.387816600304607]
#xq = X_parts[2][3]

#xq =  [-2.468; -6.576]
#xq =  [-2.331; -6.462]

#xq =  [-2.491; -6.688]
#xq =  [-2.375; -6.574]

#xq = reverse(xq)

#xq =  [-0.766; -5.782] # dark.
#xq = [-0.775; -5.671] # brighter.

# linee across boundary.
xq_st = [0.16; -4.441]
xq_fin = [1.053; -4.513]
t_xq = LinRange(0, 1.0, 100)
dir_xq = xq_fin-xq_st
Xq = collect( xq_st + t_xq[i] .* dir_xq for i = 1:length(t_xq) )

# individual models..
q_all = xx->collect( PatchMixtureKriging.queryinner(xx, X_set[m], θ, η.c_set[m]) for m = 1:length(X_set) )
q_all_evals_tuple = q_all.(Xq)

s_select_a = 3
q_evals_a = collect(q_all_evals_tuple[n][s_select_a][1] for n = 1:length(t_xq))

s_select_b = 4
q_evals_b = collect(q_all_evals_tuple[n][s_select_b][1] for n = 1:length(t_xq))

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(q_evals_a, label = "$(s_select_a)")
PyPlot.plot(q_evals_b, "--", label = "$(s_select_b)")

PyPlot.title("queries from GPs")
PyPlot.legend()

 # I am here. figure out why the weighted query has a dip.

# the weighted models.
q = xx->PatchMixtureKriging.querymixtureGP(xx, η, root, levels, radius, δ, θ, σ², weight_θ;
debug_flag = true)

q_Xq_tuple = q.(Xq)

q_evals = collect( q_Xq_tuple[n][1] for n = 1:length(t_xq))

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(q_evals)

PyPlot.title("weighted query vs. t index")
PyPlot.legend()

## point evaluate.
xq = [0.536; 4.494]
Yq1, Vq1 =  q(xq)

println("xq = ", xq)
println("Yq1 = ", Yq1)
println("Vq1 = ", Vq1)
println()

#### debug.

debug_vars_set = collect( q_Xq_tuple[n][end] for n = 1:length(t_xq))

N_regions = length(X_parts)

Ws_tilde, Us, Ps, Rs = parsedebugvarsset(N_regions, debug_vars_set, 1.0)

# diagnose weights.
Ws_tilde_display = collect( collect( Ws_tilde[n][i] for n = 1:length(t_xq) ) for i = 1:N_regions )

PyPlot.figure(fig_num)
fig_num += 1

for i = 1:N_regions
    PyPlot.plot(Ws_tilde_display[i], label = "region $(i)")
end

PyPlot.title("weights vs t index")
PyPlot.legend()

# diagnose u.
Us_display = collect( collect( Us[n][i] for n = 1:length(t_xq) ) for i = 1:N_regions )

PyPlot.figure(fig_num)
fig_num += 1

for i = 1:N_regions
    PyPlot.plot(Us_display[i], label = "region $(i)")
end

PyPlot.title("u vs t index")
PyPlot.legend()

println("q_all(Xq[86]) = ")
display(q_all(Xq[86]))

println("Us[86] = ")
display(Us[86])

println("Ps[86] = ")
display(Ps[86])

println("Rs[86] = ")
display(Rs[86])

println()

# figure out why blending is not available for one side.
# next, try adaptive RKHS.

@assert 1==2

#### visualize.
println("Global scope for debug")
p = xq

# find the region for p.
p_region_ind = PatchMixtureKriging.findpartition(p, root, levels)
println("p_region_ind = ", p_region_ind)
println()

# get all hyperplanes.
hps = PatchMixtureKriging.fetchhyperplanes(root)
@assert length(hps) == length(X_parts) - 1 # sanity check.

# radius = 0.2
# δ = 1e-5
#region_inds, ts, zs, hps_keep_flags = PatchMixtureKriging.findneighbourpartitions(p, radius, root, levels, hps, p_region_ind; δ = δ)
region_inds, ts, zs, hps_keep_flags = findneighbourpartitions(p, radius, root, levels, hps, p_region_ind; δ = δ)

# debug.
hps_kept = hps[hps_keep_flags]
hp = hps_kept[1]
u = hp.v
c = hp.c

z_kept = zs[hps_keep_flags]
t_kept = ts[hps_keep_flags]
#z = z_kept[1]

dists = collect( norm(z_kept[i]-p) for i = 1:length(z_kept) )
@assert norm( dists - abs.(t_kept) ) < 1e-10 # sanity check.

# visualize on fresh plot.
fig_num, ax = visualize2Dpartition(X_set, y_set, t_set, fig_num, "levels = $(levels)")

PyPlot.scatter([p[1];], [p[2];], marker = "x", s = 600.0, label = "p")

for i = 1:length(z_kept)
    PyPlot.scatter([z_kept[i][1];], [z_kept[i][2];], marker = "d", s = 600.0, label = "z = $(i)")
end

PyPlot.axis([p[1]-radius; p[1]+radius; p[2]-radius; p[2]+radius])
PyPlot.axis("scaled")

PyPlot.legend()

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
@testset "719.find-k-th-smallest-pair-distance.jl" begin
    @test smallest_distance_pair([1, 3, 1], 1) == 0
    @test smallest_distance_pair([1, 5, 3, 1], 3) == 2
end


text:
using Statistics
N=100000;
val=zeros(N)
for i=1:N
    gamma=-2.1 # minimum activation probability
    ActMax = 1
    agent_act =
        (
            (ActMax .^ (gamma + 1) - ActMin .^ (gamma + 1)) .* rand(N) .+
            ActMin .^ (gamma + 1)
        ) .^ (1.0 / (gamma + 1))
    val[i]=mean(agent_act)+sqrt(mean(agent_act .* agent_act))
end
# mmmean=
# s=1.5104924655094145

# s1=-0.019221243887172085

println(mean(val))