text:
module JuDGE

using MathOptInterface
using Printf
using Distributed
@everywhere using JuMP

include("tree.jl")
include("macros.jl")
include("convergence.jl")
include("risk.jl")
include("deteq.jl")

mutable struct Bounds
	UB::Float64
	LB::Float64
end

struct Column
	node::AbstractTree
	coeffs::Dict{Symbol,Any}
	obj::Float64
	var::VariableRef
	solution::Array{Float64,1}
end

struct JuDGEModel
	tree::AbstractTree
	master_problem::JuMP.Model
	sub_problems::Dict{AbstractTree,JuMP.Model}
	bounds::Bounds
	discount_factor::Float64
	master_solver#::Union{DataType,MathOptInterface.OptimizerWithAttributes,}
	probabilities::Dict{AbstractTree,Float64}
	risk::Union{Risk,Array{Risk,1}}
	sideconstraints::Union{Nothing,Function}
	log::Array{ConvergenceState,1}
	ext::Dict{Symbol,Any}
end

include("utilities.jl")

function JuDGEModel(tree::T where T <: AbstractTree, probabilities::Dict{AbstractTree,Float64}, master_problem::JuMP.Model, sub_problems::Dict{AbstractTree,JuMP.Model}, solver, bounds::Bounds, discount_factor::Float64, risk::Union{Risk,Array{Risk,1}}, sideconstraints::Union{Function,Nothing},ext::Dict{Symbol,Any})
	new_ext=Dict{Symbol,Any}()
	new_ext[:branches]=copy(ext[:branches])
	new_ext[:optimizer_settings]=deepcopy(ext[:optimizer_settings])
	JuDGEModel(tree,master_problem,sub_problems,bounds,discount_factor,solver,probabilities,risk,sideconstraints,Array{ConvergenceState,1}(),new_ext)
end

"""
	JuDGEModel(tree::AbstractTree,
               probabilities,
               sub_problem_builder::Function,
               solver;
               discount_factor=1.0,
               risk=RiskNeutral(),
               sideconstraints=nothing,
               parallel=false,
               sp_solver=nothing,
               check=true,
               perfect_foresight=false)

Define a JuDGE model.

### Required arguments
`tree` is a reference to a scenario tree

`probabilities` is either a function, which returns a dictionary of the probabilities
of all nodes in a tree, or simply the dictionary itself

`sub_problem_builder` is a function mapping a node to a JuMP model for each subproblems

`solver` is a reference to the optimizer used for the master problem (with appropriate settings);
 this can also be a tuple containing two optimizers (one for solving the relaxation, and one for
 solving the binary model)

### Optional arguments
`discount_factor` is a number between 0 and 1 defining a constant discount factor along each arc
in the scenario tree

`risk` can be either a `Risk` object, or a vector of such objects.

`sideconstraints` is a function which specifies side constraints in the master problem, see
[Tutorial 9: Side-constraints](@ref) for further details

`parallel` is a boolean, setting whether the sub-problems will be formulated in parallel

`sp_solver` if formulating the sub-problems in parallel, some solvers have issues; in these case
the sub-problem solver can be set using this argument

`check` is a boolean, which can be set to `false` to disable the validation of the JuDGE model.

`perfect_foresight` is a boolean; this is an experimental feature, which creates an array of
JuDGE models, one for each leaf node. This will enable users to easily compute the EVPI for
the stochastic program. Also can be used for regret-based risk implementations.

### Examples
	judge = JuDGEModel(tree, ConditionallyUniformProbabilities, sub_problems,
                                    Gurobi.Optimizer)
	judge = JuDGEModel(tree, probabilities, sub_problems, CPLEX.Optimizer,
                                    discount_factor=0.9, risk=(0.5,0.1)))
"""
function JuDGEModel(tree::T where T <: AbstractTree, probabilities, sub_problem_builder::Function, solver; discount_factor::Float64=1.0, risk::Union{Risk,Array{Risk,1}}=RiskNeutral(), sideconstraints::Union{Function,Nothing}=nothing, parallel::Bool=false, sp_solver=nothing, check::Bool=true, perfect_foresight::Bool=false)
	println("")
	if !perfect_foresight
		println("Establishing JuDGE model for tree: " * string(tree))
	else
		println("Establishing perfect foresight models for tree: " * string(tree))
	end

	if typeof(probabilities) <: Function
		probabilities = probabilities(tree)
	end
	if typeof(probabilities)!=Dict{AbstractTree,Float64}
		error("\'probabilities\' needs to be a dictionary mapping AbstractTree to Float64\nor a function that generates such a dictionary")
	end

	nodes=collect(tree)
	if parallel
		sps=pmap(sub_problem_builder,nodes)
		i=1
		sub_problems=Dict()
		for n in nodes
			sub_problems[n]=sps[i]
			i+=1
		end
	else
		sub_problems = Dict(i => sub_problem_builder(i) for i in nodes)
	end

	if sp_solver!=nothing
		for n in nodes
			set_optimizer(sub_problems[n],sp_solver)
		end
	end

	if check
		print("Checking sub-problem format...")
		check_specification_is_legal(sub_problems)
		println("Passed")
	else
		println("Skipping checks of sub-problem format")
	end
	scale_objectives(tree,sub_problems,discount_factor)

	if !perfect_foresight
		print("Building master problem...")
		master_problem = build_master(sub_problems, tree, probabilities, solver, discount_factor, risk, sideconstraints)
		println("Complete")
		ext=Dict{Symbol,Any}()
		ext[:branches]=Branch[]
		ext[:optimizer_settings]=Dict{Symbol,Any}()
		return JuDGEModel(tree,master_problem,sub_problems,Bounds(Inf,-Inf),discount_factor,solver,probabilities,risk,sideconstraints,Array{ConvergenceState,1}(),ext)
	else
		scenarios=Dict{AbstractTree,JuDGEModel}()
		pr=Dict{AbstractTree,Float64}()
		print("Building master problems...")
		scen_trees=get_scenarios(tree)
		for t in scen_trees
			leaf=getID(get_leafnodes(t)[1])
			sps = Dict(i => sub_problems[getID(i)] for i in collect(t))
			probs = Dict(i => 1.0 for i in collect(t))
			master_problem = build_master(sps, t, probs, solver, discount_factor, risk, sideconstraints)
			ext=Dict{Symbol,Any}()
			ext[:branches]=Branch[]
			ext[:optimizer_settings]=Dict{Symbol,Any}()
			scenarios[leaf]=JuDGEModel(t,master_problem,sps,Bounds(Inf,-Inf),discount_factor,solver,probabilities,risk,sideconstraints,Array{ConvergenceState,1}(),ext)
			pr[leaf]=probabilities[leaf]
		end
		println("Complete")
		return scenarios, pr
	end
end

include("branchandprice.jl")
include("master.jl")

function scale_objectives(tree::T where T <: AbstractTree,sub_problems,discount_factor::Float64)
	if discount_factor <= 0.0 || discount_factor > 1.0
		error("discount_factor must be greater than 0.0 and less than or equal to 1.0")
	end

	for (node,sp) in sub_problems
		if !haskey(sp.ext, :capitalcosts)
			sp.ext[:capitalcosts]=Dict()
			sp.ext[:capitalcosts][:constant]=0
		end
		if !haskey(sp.ext, :ongoingcosts)
			sp.ext[:ongoingcosts]=Dict()
			sp.ext[:ongoingcosts][:constant]=0
		end

		sp.ext[:objective]=@variable(sp, obj)
		sp.ext[:objective_con]=@constraint(sp, sp.ext[:objective]-objective_function(sp) == 0)

		@objective(sp, Min, 0.0)
		set_normalized_coefficient(sp.ext[:objective_con],sp.ext[:objective],1.0/(discount_factor^depth(node)))
	end
end

function Base.map(f, hello::JuMP.Containers.DenseAxisArray)
	JuMP.Containers.DenseAxisArray(map(f, hello.data), deepcopy(hello.axes),deepcopy(hello.lookup))
end

function Base.map(f, hello::JuMP.Containers.SparseAxisArray)
	JuMP.Containers.SparseAxisArray(Dict(i => f(hello[i]) for i in keys(hello.data)))
end

function add_variable_as_column(master, column)
	for constr in master.ext[:convexcombination][column.node]
		set_normalized_coefficient(constr, column.var, 1.0)
	end

	for (var,val) in column.coeffs
		if typeof(val)==Float64
			set_normalized_coefficient(master.ext[:coverconstraint][column.node][var], column.var, val)
		else
			for (i,val2) in val
				set_normalized_coefficient(master.ext[:coverconstraint][column.node][var][i], column.var, val2)
			end
		end
	end

	for node in collect(column.node)
		if typeof(node)==Leaf
			set_normalized_coefficient(master.ext[:scenprofit_con][node], column.var, column.obj)
		end
	end
end

# constraint that's used to recreate a discrete subproblem variable in the master problem
function add_mixed_cover(master, sp, column)
	if !(column.node in keys(master.ext[:discrete_con]))
		master.ext[:discrete_var][column.node]=Dict{Int64,VariableRef}()
		master.ext[:discrete_con][column.node]=Dict{Int64,ConstraintRef}()
		for i in 1:length(sp.ext[:discrete])
			master.ext[:discrete_var][column.node][i]=@variable(master)
			master.ext[:discrete_con][column.node][i]=@constraint(master,0.0==master.ext[:discrete_var][column.node][i])
		end
	end

	for i in 1:length(sp.ext[:discrete])
		set_normalized_coefficient(master.ext[:discrete_con][column.node][i],column.var,column.solution[i])
	end
end

function add_column(master::JuMP.Model, sub_problem::JuMP.Model, node::AbstractTree; branches=nothing)
	coeffs=Dict{Symbol,Any}()

	sol=[0.0]

	if sub_problem.ext[:form]==:mixed
		sol=JuMP.value.(sub_problem.ext[:discrete])
	end

	for (name,variable) in sub_problem.ext[:expansions]
		if typeof(variable) <: AbstractArray
			coeffs[name]=Dict{Any,Float64}()
			vals=JuMP.value.(variable)
			for i in eachindex(variable)
				if sub_problem.ext[:options][name][4]!=:Con
					rval=round(vals[i])
					if rval!=0.0
						coeffs[name][i]=rval
					end
				else
					coeffs[name][i]=vals[i]
				end
			end
		else
			if sub_problem.ext[:options][name][4]!=:Con
				rval=round(JuMP.value(variable))
				if rval!=0.0
					coeffs[name]=rval
				end
			else
				coeffs[name]=JuMP.value(variable)
			end
		end
	end

	UB = haskey(node.ext,:col_max) ? node.ext[:col_max] : 1.0

	column=Column(node,coeffs,JuMP.value(sub_problem.ext[:objective]),
					JuMP.add_variable(master, JuMP.build_variable(error, UnitIntervalInformation(UB=UB))),sol)

	if branches!=nothing
		for b in branches
			if b.filter!=nothing && b.filter(column)==:ban
				set_upper_bound(column.var,0.0)
				println("New column banned on current branch")
				break
			end
		end
	end

	add_variable_as_column(master, column)
	push!(master.ext[:columns][node],column)
	if sub_problem.ext[:form]==:mixed
		add_mixed_cover(master, sub_problem, column)
	end
	column
end

"""
	solve(judge::JuDGEModel;
	      termination::Termination=Termination(),
	      max_no_int::Int=2147483647,
	      blocks::Union{Nothing,Array{Array{AbstractTree,1},1}}=nothing,
	      warm_starts::Bool=false,
	      optimizer_attributes::Union{Nothing,Function}=nothing,
	      mp_callback::Union{Nothing,Function}=nothing,
	      prune::Float64=Inf,
	      heuristic::Union{Nothing,Function}=nothing,
	      verbose::Int=2)

Solve a JuDGEModel `judge` without branch-and-price.

### Required Arguments
`judge` is the JuDGE model that we wish to solve.

### Optional Arguments

`termination` is a `Termination` object containing all the stopping conditions.

`max_no_int` is the maximum number of iterations yielding a fractional solution before a MIP solve is
performed on the master. By default, the MIP solves will not occur until the relaxed bound gap is less
than the `relgap` / `absgap` stopping conditions. To override this, set `max_no_int` to the negative
of the number desired value.

`blocks` specifies the groups of nodes to solve in each iteration (these groups can be generated using
`JuDGE.get_groups()``, or created manually), after all nodes have been solved, a full pricing iteration
is used to compute an updated lower bound. See `advanced.jl` for more details.

`warm_starts` boolean specifing whether to use warm starts for subproblems and binary solves of master
problem.

`optimizer_attributes` can be set to a specific function that dynamically changes optimizer attributes
for the subproblems; this should only be used by people who have examined the `advanced.jl` example.

`mp_callback` is a user-defined function that specifies termination conditions for MIP solves of the
master problem. See examples/advanced.jl.

`prune` is used to stop the algorithm before convergence, if a known upper bound for the problem is
specified.

`heuristic` is a user-defined function that typically would perform an improvement heuristic when a
new incumbent is found.

`verbose` if 0, all output from solve will be suppressed, if 1, the subproblem solve process will be
suppressed. Default is 2.

### Examples
    JuDGE.solve(jmodel, termination=Termination(rlx_abstol=10^-6))
	JuDGE.solve(jmodel, termination=Termination(rlx_abstol=10^-6), max_no_int=-5)
"""
function solve(judge::JuDGEModel;
   termination::Termination=Termination(),
   max_no_int::Int=2147483647,
   blocks::Union{Nothing,Array{Array{AbstractTree,1},1}}=nothing,
   warm_starts::Bool=false,
   optimizer_attributes::Union{Nothing,Function}=nothing,
   mp_callback::Union{Nothing,Function}=nothing,
   prune::Float64=Inf,
   heuristic::Union{Nothing,Function}=nothing,
   verbose::Int=2)

	current = InitialConvergenceState()
	empty!(judge.log)
	push!(judge.log,current)
	if verbose>0
		print("")
		print("")
		println("Relaxed ObjVal  |   Upper Bound   Lower Bound  |  Absolute Diff   Relative Diff  |  Fractional  |      Time     Iter")
	end
	# set up times for use in convergence
	initial_time = time()
	obj = Inf
	heuristicobj=Inf
	nodes=collect(judge.tree)
	objduals=Dict{AbstractTree,Float64}()
	redcosts=Dict{AbstractTree,Float64}()
	for node in nodes
		objduals[node]=-1
		redcosts[node]=-1
	end
	no_int_count=0
	if optimizer_attributes!=nothing
		optimizer_attributes(judge,false,true)
	end

	max_char=length(nodes[end].name)+length(string(length(nodes)))
	function get_whitespace(name::String,number::Int64)
		blank="  "
		spaces=max_char-length(name)-length(string(number))
		for i in 1:spaces
			blank*=" "
		end
		blank
	end
	if blocks==nothing
		blocks=[nodes]
	end

	block=-1

	if judge.master_problem.ext[:mip]
		remove_binary(judge)
	end
	set_banned_variables!(judge)

	remove_branch_constraints!(judge)
	b_con=add_branch_constraints!(judge)

	optimize!(judge.master_problem)
	status=termination_status(judge.master_problem)
	if status==MOI.NUMERICAL_ERROR
		println("\nMaster problem returned a MOI.NUMERICAL_ERROR")
		return
	end

	while true
		if block<=0
			nodes2 = nodes
		else
			nodes2 = blocks[block]
		end
		# perform the main iterations
		for i in 1:length(nodes2)
			node=nodes2[i]
			sp=judge.sub_problems[node]
			updateduals(judge.master_problem, sp, node, status, current.iter)
			if verbose==2
				overprint("Solving subproblem for node "*node.name*get_whitespace(node.name,i)*string(i)*"/"*string(length(nodes2)))
			end

			optimize!(sp)

			if termination_status(sp)!=MathOptInterface.OPTIMAL && termination_status(sp)!=MathOptInterface.INTERRUPTED && termination_status(sp)!=MathOptInterface.TIME_LIMIT
				@warn("Solve for subproblem "*node.name*" exited with status "*string(termination_status(sp)))
				println("Subproblem is infeasible or unbounded")
				for (sp2,con) in b_con
					delete(sp2,con)
				end
				return :sp_infeasible
			end
			if warm_starts
				vars=all_variables(sp)
				set_start_value.(vars, JuMP.value.(vars))
			end
		end
		if verbose==2
			overprint("")
		end
		frac=0
		if status!=MathOptInterface.INFEASIBLE_OR_UNBOUNDED && status!=MathOptInterface.INFEASIBLE && status!=MathOptInterface.DUAL_INFEASIBLE
			for node in nodes2
				objduals[node]=objective_bound(judge.sub_problems[node])
				redcosts[node]=objective_value(judge.sub_problems[node])
			end
			if status!=MathOptInterface.OPTIMAL
				@warn("Master problem did not solve to optimality: "*string(status))
			elseif block==0
				getlowerbound(judge,objduals)
			end
		elseif judge.bounds.LB>-Inf
			println("\nMaster problem is infeasible or unbounded")
			for (sp,con) in b_con
				delete(sp,con)
			end
			judge.bounds.LB=Inf
			judge.bounds.UB=Inf
			return
		end

		num_var=num_variables(judge.master_problem)
		for node in nodes2
			if redcosts[node]<-10^-10 || status==MathOptInterface.INFEASIBLE_OR_UNBOUNDED || status==MathOptInterface.INFEASIBLE || status==MathOptInterface.DUAL_INFEASIBLE
				column = add_column(judge.master_problem, judge.sub_problems[node], node)
				if warm_starts
					set_start_value(column.var,0.0)
				end
			end
		end

		optimize!(judge.master_problem)
		status=termination_status(judge.master_problem)
		if status==MOI.NUMERICAL_ERROR
			println("\nMaster problem returned a MOI.NUMERICAL_ERROR")
			for (sp,con) in b_con
				delete(sp,con)
			end
			return
		elseif status!=MathOptInterface.INFEASIBLE_OR_UNBOUNDED && status!=MathOptInterface.INFEASIBLE && status!=MathOptInterface.DUAL_INFEASIBLE
			frac=fractionalcount(judge,termination.inttol)
			obj=objective_value(judge.master_problem)
			if frac==0
				if heuristic!=nothing && heuristicobj>obj
					if heuristic(judge)<0.0 && termination.allow_frac ∉ [:first_fractional,:no_binary_solve]
						solve_master_binary(judge,initial_time,termination,warm_starts,nothing,verbose)
					else
						optimize!(judge.master_problem)
					end
					obj=objective_value(judge.master_problem)
					heuristicobj=obj
				end
				judge.bounds.UB=obj

				no_int_count=0
			else
				no_int_count+=1
			end
		end
		current = ConvergenceState(obj, judge.bounds.UB, judge.bounds.LB, time() - initial_time, current.iter + 1, frac)
		if verbose>0
			display(current)
		end
		push!(judge.log,current)
		if prune<judge.bounds.LB
			if verbose>0
				println("\nDominated by incumbent.")
			end
			for (sp,con) in b_con
				delete(sp,con)
			end
			return
		elseif has_converged(termination, current)
			if frac>0
				solve_master_binary(judge,initial_time,termination,warm_starts,nothing,verbose)
			end
			if heuristic!=nothing
				if heuristic(judge)<0.0 && termination.allow_frac ∉ [:first_fractional,:no_binary_solve]
					solve_master_binary(judge,initial_time,termination,warm_starts,nothing,verbose)
					heuristicobj=judge.bounds.UB
				else
					optimize!(judge.master_problem)
				end
			end
			if verbose>0
				println("\nConvergence criteria met.")
			end
			for (sp,con) in b_con
				delete(sp,con)
			end
			return
		elseif termination.allow_frac==:first_fractional && frac>0
			if verbose>0
				println("\nFractional solution found.")
			end
			for (sp,con) in b_con
				delete(sp,con)
			end
			return
		elseif ((max_no_int>0 && no_int_count >= max_no_int && current.num_frac > 0) && (current.rlx_abs<termination.abstol || current.rlx_rel<termination.reltol)) || (max_no_int<0 && no_int_count>=-max_no_int && current.num_frac > 0)
			current=solve_master_binary(judge,initial_time,termination,warm_starts,mp_callback,verbose)
			if heuristic!=nothing && heuristicobj>judge.bounds.UB
				if heuristic(judge)<0.0
					current=solve_master_binary(judge,initial_time,termination,warm_starts,mp_callback,verbose)
				else
					optimize!(judge.master_problem)
				end
				heuristicobj=judge.bounds.UB
			end

			if has_converged(termination, current)
				if verbose>0
					println("\nConvergence criteria met.")
				end
				for (sp,con) in b_con
					delete(sp,con)
				end
				return
			elseif termination.allow_frac==:binary_solve
				remove_binary(judge)
				optimize!(judge.master_problem)
			end
			no_int_count=0
		end

		if optimizer_attributes==nothing
			if block==0 && num_var==num_variables(judge.master_problem)
				solve_master_binary(judge,initial_time,termination,warm_starts,nothing,verbose)
				if verbose>0
					println("\nStalled: exiting.")
				end
				for (sp,con) in b_con
					delete(sp,con)
				end
				return
			end
		elseif optimizer_attributes(judge,num_var==num_variables(judge.master_problem),false)
			if verbose>0
				println("\nStalled: exiting.")
			end
			for (sp,con) in b_con
				delete(sp,con)
			end
			return
		end
		if length(blocks)==1
			block=0
		else
			block = (block+1) % (length(blocks)+1)
		end
	end
	for (sp,con) in b_con
		delete(sp,con)
	end
end

function getlowerbound(judge::JuDGEModel,objduals::Dict{AbstractTree,Float64})
	lb = objective_value(judge.master_problem)
	for (i,objdual) in objduals
		if :sum_max in keys(i.ext)
			lb+=objdual*i.ext[:sum_max]
		else
			lb+=objdual
		end
	end
	if lb>judge.bounds.LB
		judge.bounds.LB=lb
	end
end

function fractionalcount(jmodel::JuDGEModel,inttol::Float64)
	count=0

	for (node,sp) in jmodel.sub_problems
		if sp.ext[:form]==:binarycolumns
			for col in jmodel.master_problem.ext[:columns][node]
				val=JuMP.value(col.var)
				count += min(val-floor(val),ceil(val)-val) > inttol ? 1 : 0
			end
			return count
		end
	end

	for node in collect(jmodel.tree)
		for x in keys(jmodel.master_problem.ext[:expansions][node])
			if jmodel.sub_problems[jmodel.tree].ext[:options][x][4]!=:Con
				var = jmodel.master_problem.ext[:expansions][node][x]
				if typeof(var) <: AbstractArray
					for key in eachindex(var)
						val=JuMP.value(var[key])
						count += min(val-floor(val),ceil(val)-val) > inttol ? 1 : 0
					end
				else
					val=JuMP.value(var)
					count += min(val-floor(val),ceil(val)-val) > inttol ? 1 : 0
				end
			end
		end
		if jmodel.sub_problems[node].ext[:form]==:mixed
			for (i,var) in jmodel.master_problem.ext[:discrete_var][node]
				val=JuMP.value(var)
				count += min(val-floor(val),ceil(val)-val) > inttol ? 1 : 0
			end
		end
	end

	count
end

function updateduals(master, sub_problem, node, status, iter)
	if status!=MathOptInterface.INFEASIBLE_OR_UNBOUNDED && status!=MathOptInterface.INFEASIBLE && status!=MathOptInterface.DUAL_INFEASIBLE
		for (name,var) in sub_problem.ext[:expansions]
			if typeof(var) <: AbstractArray
				for i in eachindex(var)
					set_objective_coefficient(sub_problem, var[i], -dual(master.ext[:coverconstraint][node][name][i]))
				end
			else
				set_objective_coefficient(sub_problem, var, -dual(master.ext[:coverconstraint][node][name]))
			end
		end
		nodes=collect(node)
		total=0.0
		for n in nodes
			if typeof(n)==Leaf
				total-=dual(master.ext[:scenprofit_con][n])
			end
		end
		set_objective_coefficient(sub_problem,sub_problem.ext[:objective],total)

		cc_sum=0.0
		for constr in master.ext[:convexcombination][node]
			set=typeof(constraint_object(constr).set)
			if set <: MathOptInterface.LessThan
				cc_sum+=dual(constr)
			elseif set <: MathOptInterface.GreaterThan
				cc_sum-=dual(constr)
			elseif set <: MathOptInterface.EqualTo
				cc_sum+=dual(constr)
			end
		end

		set_objective_function(sub_problem,
					objective_function(sub_problem)-objective_function(sub_problem).constant
					                               -cc_sum)
	else
		if iter%2==0
			oc=-10.0^10
		else
			oc=10.0^10
		end
		for (name,var) in sub_problem.ext[:expansions]
			if typeof(var) <: AbstractArray
				for i in eachindex(var)
					set_objective_coefficient(sub_problem, var[i], oc)
				end
			else
				set_objective_coefficient(sub_problem, var, oc)
			end
		end
		set_objective_coefficient(sub_problem,sub_problem.ext[:objective], 1.0)
	end
end

function get_objval(jmodel::JuDGEModel; risk=jmodel.risk)
	scenario_objs=Dict{Leaf,Float64}()

	for (leaf, var) in jmodel.master_problem.ext[:scenprofit_var]
		scenario_objs[leaf]=JuMP.value(var)
	end

	compute_objval(scenario_objs, jmodel.probabilities, risk)
end

"""
	resolve_subproblems(judge::JuDGEModel)

Once a JuDGE model has converged, it is necessary to re-solve the subproblems to find the optimal decisions within each node.

### Required Arguments
`jmodel` is the JuDGE model that we wish to solve.

### Examples
    resolve_subproblems(judge)
"""
function resolve_subproblems(jmodel::JuDGEModel)
	fix_expansions(jmodel)
	resolve_fixed(jmodel)
end

function fix_expansions(jmodel::JuDGEModel)
	if termination_status(jmodel.master_problem) != MathOptInterface.OPTIMAL && termination_status(jmodel.master_problem) != MathOptInterface.INTERRUPTED && termination_status(jmodel.master_problem) != MathOptInterface.LOCALLY_SOLVED  && termination_status(jmodel.master_problem) != MathOptInterface.INTERRUPTED
		error("You need to first solve the decomposed model.")
	end

	for node in collect(jmodel.tree)
		sp=jmodel.sub_problems[node]
		set_objective_function(sp, objective_function(sp)-objective_function(sp).constant)
		set_objective_coefficient(sp, sp.ext[:objective], 1.0)
		for (name,var) in jmodel.master_problem.ext[:expansions][node]
			var2=sp.ext[:expansions][name]
			if typeof(var) <: AbstractArray
				for i in eachindex(var)
					value=0.0
					if sp.ext[:options][name][1]==:state
						prev = node.parent
						if prev==nothing
							var3=jmodel.master_problem.ext[:expansions][node][name][i]
							value=JuMP.value(var3)-sp.ext[:options][name][7]
						else
							var3=jmodel.master_problem.ext[:expansions][node][name][i]
							var4=jmodel.master_problem.ext[:expansions][prev][name][i]
							value=JuMP.value(var3)-JuMP.value(var4)
						end
					else
						con_obj=constraint_object(jmodel.master_problem.ext[:coverconstraint][node][name][i])
						for prev in history(node)
							var3=jmodel.master_problem.ext[:expansions][prev][name][i]
							if var3 in keys(con_obj.func.terms)
								value+=JuMP.value(var3)*-con_obj.func.terms[var3]
							end
						end
					end
					if sp.ext[:options][name][1]==:shutdown
						JuMP.set_lower_bound(var2[i],value)
					elseif sp.ext[:options][name][1]==:expansion
						JuMP.set_upper_bound(var2[i],value)
					elseif sp.ext[:options][name][1]==:enforced || sp.ext[:options][name][1]==:state
						JuMP.fix(var2[i],value,force=true)
					end
					set_objective_coefficient(sp, var2[i], 0.0)
				end
			elseif isa(var,VariableRef)
				value=0.0
				if sp.ext[:options][name][1]==:state
					prev = node.parent
					if prev==nothing
						var3=jmodel.master_problem.ext[:expansions][node][name]
						value=JuMP.value(var3)-sp.ext[:options][name][7]
					else
						var3=jmodel.master_problem.ext[:expansions][node][name]
						var4=jmodel.master_problem.ext[:expansions][prev][name]
						value=JuMP.value(var3)-JuMP.value(var4)
					end
				else
					con_obj=constraint_object(jmodel.master_problem.ext[:coverconstraint][node][name])
					for prev in history(node)
						var3=jmodel.master_problem.ext[:expansions][prev][name]
						if var3 in keys(con_obj.func.terms)
							value+=JuMP.value(var3)*-con_obj.func.terms[var3]
						end
					end
				end
				if sp.ext[:options][name][1]==:shutdown
					JuMP.set_lower_bound(var2,value)
				elseif sp.ext[:options][name][1]==:expansion
					JuMP.set_upper_bound(var2,value)
				elseif sp.ext[:options][name][1]==:enforced || sp.ext[:options][name][1]==:state
					JuMP.fix(var2,value,force=true)
				end
				set_objective_coefficient(sp, var2, 0.0)
			end
		end
	end
end

function resolve_fixed(jmodel::JuDGEModel)
	for n in collect(jmodel.tree)
		JuMP.optimize!(jmodel.sub_problems[n])
	end

	scenario_objs=Dict{Leaf,Float64}()
	for (leaf,con) = jmodel.master_problem.ext[:scenprofit_con]
		obj=0.0
		for node in history(leaf)
			obj+=objective_value(jmodel.sub_problems[node])
			for key in keys(jmodel.master_problem.ext[:expansions][node])
				var = jmodel.master_problem.ext[:expansions][node][key]
				if isa(var,VariableRef)
					obj+=JuMP.value(var)*normalized_coefficient(con, var)
				elseif typeof(var) <: AbstractArray
					for v in eachindex(var)
						obj+=JuMP.value(var[v])*normalized_coefficient(con,var[v])
					end
				end
			end
		end
		scenario_objs[leaf]=obj;
	end

	compute_objval(scenario_objs, jmodel.probabilities, jmodel.risk)
end

include("model_verification.jl")

# pretty printing
function Base.show(io::IO, ::MIME"text/plain", judge::JuDGEModel)
	print(io, "JuDGE Model with:\n")
	println(io, "  Tree: ", judge.tree)
	print(io, "  Expansion variables: ")
	keys = judge.sub_problems[judge.tree].ext[:expansions]
	for (i,ii) in keys
		print(io, "$(i) ")
	end
end

function Base.show(io::IO, judge::JuDGEModel)
	print(io, "JuDGE Model with:\n")
	println(io, "  Tree: ", judge.tree)
	print(io, "  Expansion variables: ")
	keys = judge.sub_problems[judge.tree].ext[:expansions]
	for (i,ii) in keys
		print(io, "$(i) ")
	end
end

include("output.jl")

export @expansion, @shutdown, @enforced, @state, @expansionconstraint, @capitalcosts, @ongoingcosts, JuDGEModel, Risk, RiskNeutral, Leaf, Tree, AbstractTree, narytree, ConditionallyUniformProbabilities, UniformLeafProbabilities, get_node, tree_from_leaves, tree_from_nodes, tree_from_file, DetEqModel, resolve_subproblems, Termination

end


text:
using DashBootstrapComponents

items = [
    dbc_dropdownmenuitem("First"),
    dbc_dropdownmenuitem(divider = true),
    dbc_dropdownmenuitem("Second"),
];

dropdown = dbc_row(
    [
        dbc_col(
            dbc_dropdownmenu(
                label = "Dropdown (default)",
                children = items,
                direction = "down",
            ),
            width = "auto",
        ),
        dbc_col(
            dbc_dropdownmenu(label = "Dropstart", children = items, direction = "start"),
            width = "auto",
        ),
        dbc_col(
            dbc_dropdownmenu(label = "Dropend", children = items, direction = "end"),
            width = "auto",
        ),
        dbc_col(
            dbc_dropdownmenu(label = "Dropup", children = items, direction = "up"),
            width = "auto",
        ),
    ],
    justify = "between",
);

ext: jl

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
    simplicial_subdivision(k, E)

Generate generic rules for generating vertices when splitting an E-dimensional
simplex with a size reducing factor of k.

Note: Subsimplices containing any of the original vertices is not considered.

This function comes in two versions. Both versions returns information about
all the subsimplices of the splitted original simplex. They differ in what
information they return about the vertices furnishing these simplices.

The SimplicialSubdivisionMultiple function returns an array of vertices that
containONLY THE NEW VERTICES. This is to prevent duplication when considering
triangulations, in which simplices might share vertices.

SimplicialSubdivisionSingle(k, E) returns an array of all the vertices
furnishing the subidivded simplex, including the vertices of the original
simplex. Use this function when you want to split a single simplex.

SimplexSplittingOutput the decomposition part of a SimplexSplitting output.
This is a (k^E * (E+1) x k) array of integers. (so, ignore the orientations
	for now)
"""
function simplicial_subdivision(k::Int, E::Int)

    # Get simplex splitting rules, ignore orientations for now.
    splits::Array{Int, 2} = simplex_split(k, E, orientations = false)

    # Retrieve integer labels from the tensor decomposition. This is a a
    #  way of uniquely identifying the rows of the 'splits' array. Each integer
    # corresponds to a row in the splits array (but are not necessarily in order).
    integer_labels::Vector{Int} = (splits .- 1) * ((E+1) .^ collect(0:(k-1)))

    # Find the unique labels and sort them
    uniquelabels = sort(unique(integer_labels))

    # How many unique rows are there?
    n_rows = size(uniquelabels, 1)

    J = indexin(integer_labels, uniquelabels)

    # IndicesOfUnqiueIntegersinOriginal
    I = Vector{Int}(undef, length(uniquelabels))
    for i = 1:length(uniquelabels)
        I[i] = something(findfirst(isequal(uniquelabels[i]), integer_labels), 0)
    end
    V::Array{Int, 2} = repeat(uniquelabels, 1, E + 1)

    tmp = repeat((((E+1)^k - 1) ./ E) * transpose(collect(0:E)), n_rows, 1)

    tmp = heaviside0(-abs.(V - tmp)) .* repeat(collect(1:n_rows), 1, E + 1)

	#@show tmp
    #@show findall(tmp)
    #tmp = tmp[findall(tmp)]
    tmp = tmp[tmp .> 0]

    aux = round.(Int, tmp)
    Caux = round.(Int, complementary(tmp, n_rows))
    Ipp = [aux; Caux]

    Ip = zeros(Int, size(I))

    # for any i\in Ipp, Ip(i) is the number j\in {1,...,size(Ipp,1)} such that
    # Ipp(j)=i
    Ip[Ipp] = transpose(collect(1:size(I, 1)))

    # Indices of the non-repeated vertices that are not part of the original
    # simplex. These are the truly new vertices generated by the splitting.
    new_vertices = splits[I, :][Caux, :]

    # Subarrays of dimension E+1 x 1 vertically concatenated. Each subarray
    # contains the indices in the rows of NewVert corresponding to the vertices
    # furnishing the corresponding simplex (analogous to the output of
    # SimplexSplitting, check notes).
    subtriangulation_inds = Ip[J]

    subtriangulation_inds = copy(transpose(reshape(subtriangulation_inds, E + 1, k^E)))

    return new_vertices, subtriangulation_inds
end


text:
function cart2sphere(x,y,z)
r=sqrt(x^2+y^2+z^2);
phi=asin(z/r);

# ϕ = atan(z, hypot(y, x))
# if abs(ϕ) == 90
#     λ = zero(ϕ)
# else
#     λ = atan(y, x)
# end

lam=0;
if abs(abs(phi)-pi/2)>1.e-14
  lam=atan2(y,x); # TODO: check translation with Oswald
  if lam<0.0
    lam=lam+2*pi;
  end
end
return (lam,phi,r)
end


text:
# RolloutSimulator
# maintained by @zsunberg

"""
A fast simulator that just returns the reward

The simulation will be terminated when either
1) a terminal state is reached (as determined by `isterminal()` or
2) the discount factor is as small as `eps` or
3) max_steps have been executed

Keyword Arguments:
    - eps
    - max_steps

Usage (optional arguments in brackets):
    ro = RolloutSimulator()
    history = simulate(ro, pomdp, policy, [updater [, init_belief [, init_state]]])
"""
struct RolloutSimulator{RNG<:AbstractRNG} <: Simulator
    rng::RNG

    # optional: if these are null, they will be ignored
    max_steps::Nullable{Int}
    eps::Nullable{Float64}

    # DEPRECATED: remove in v0.2.8 or higher
    initial_state::Nullable{Any}
end

# These are the only safe constructors to use
RolloutSimulator(rng::AbstractRNG, d::Int=typemax(Int)) = RolloutSimulator(alphas, rng, Nullable{Int}(d), Nullable{Float64}(), Nullable{Any}())
function RolloutSimulator(;rng=MersenneTwister(rand(UInt32)),
                           initial_state=Nullable{Any}(),
                           eps=Nullable{Float64}(),
                           max_steps=Nullable{Int}())
    if !isnull(initial_state)
        warn("The initial_state argument for RolloutSimulator is deprecated. The initial state should be specified as the last argument to simulate(...).")
    end
    return RolloutSimulator{typeof(rng)}(rng, max_steps, eps, initial_state)
end

# COMPATIBILITY CONSTRUCTOR: DO NOT USE!
# Once version 0.2.7 is registered, start having this throw a warning in master
function RolloutSimulator(rng::AbstractRNG, is::Nullable{Any}, eps::Nullable{Float64}, ms::Nullable{Int})
    return RolloutSimulator(rng, ms, eps, is)
end

# @POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy) begin
#     @req updater(::typeof(policy))
#     bu = updater(policy)
#     @subreq simulate(sim, pomdp, policy, bu)
# end
#
# @POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, bu::Updater) begin
#     @req initial_state_distribution(::typeof(pomdp))
#     dist = initial_state_distribution(pomdp)
#     @subreq simulate(sim, pomdp, policy, bu, dist)
# end

function simulate(sim::RolloutSimulator, pomdp::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}, policy::Policy, bu::Updater=updater(policy))
    dist = initial_state_distribution(pomdp)
    return simulate(sim, pomdp, policy, bu, dist)
end

# @POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief) begin
#     @req rand(::typeof(sim.rng), ::typeof(initial_belief))
#     @subreq simulate(sim, pomdp, policy, updater, initial_belief, s)
# end

function simulate{S}(sim::RolloutSimulator, pomdp::Union{POMDP{S},IPOMDP{S},RPOMDP{S},RIPOMDP{S}}, policy::Policy, updater::Updater, initial_belief)

    if !isnull(sim.initial_state)
        s = convert(S, get(sim.initial_state))::S
    else
        s = rand(sim.rng, initial_belief)::S
    end

    return simulate(sim, pomdp, policy, updater, initial_belief, s)
end

# @POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief, s) begin
#     P = typeof(pomdp)
#     S = state_type(P)
#     A = action_type(P)
#     O = obs_type(P)
#     @req initialize_belief(::typeof(updater), ::typeof(initial_belief))
#     @req isterminal(::P, ::S)
#     @req discount(::P)
#     @req generate_sor(::P, ::S, ::A, ::typeof(sim.rng))
#     b = initialize_belief(updater, initial_belief)
#     @req action(::typeof(policy), ::typeof(b))
#     @req update(::typeof(updater), ::typeof(b), ::A, ::O)
# end

function minutil_gen(prob::Union{RPOMDP,RIPOMDP}, b::Vector{Float64}, a, alphavecs::Vector{Vector{Float64}})
    nz = n_observations(prob)
    ns = n_states(prob)
    nα = length(alphavecs)
    aind = action_index(prob, a)
    plower, pupper = dynamics(prob)
    m = Model(solver = ClpSolver())
    @variable(m, u[1:nz])
    @variable(m, p[1:ns, 1:nz, 1:ns])
    @objective(m, Min, sum(u))
    for zind = 1:nz, αind = 1:nα
        @constraint(m, u[zind] >= sum(b[sind] * p[:,zind,sind]' * alphavecs[αind] for sind = 1:ns))
    end
    @constraint(m, p .<= pupper[:,:,:,aind])
    @constraint(m, p .>= plower[:,:,:,aind])
    for sind = 1:ns
        @constraint(m, sum(p[:,:,sind]) == 1)
    end
    JuMP.solve(m)
    getvalue(u), getvalue(p)
end

function generate_sor_worst(prob::Union{RPOMDP,RIPOMDP}, b, s, a, rng::AbstractRNG, alphavecs)
    u, p = minutil_gen(prob, b, a, alphavecs)
    sind = state_index(prob, s)
    denom = sum(p[:,:,sind])
    tarray = [sum(p[spi,:,sind]) for spi in 1:n_states(prob)] ./ denom
    tdist = SparseCat(states(prob), tarray)
    # @show tdist
    sp = rand(rng, tdist)
    oarray = [sum(p[:,zind,sind]) for zind in 1:n_observations(prob)] ./ denom
    odist = SparseCat(observations(prob), oarray)
    # @show odist
    o = rand(rng, odist)
    r = reward(prob, b, s, a, sp)
    sp, o, r
end

function simulate(sim::RolloutSimulator, pomdp::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}, policy::Policy, updater::Updater, initial_belief, s; verbose = false)
    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))
    disc = 1.0
    # disc = discount(pomdp)
    r_total = 0.0
    b = initialize_belief(updater, initial_belief)
    step = 1
    bcorrect = 0
    while disc > eps && !isterminal(pomdp, s) && step <= max_steps # TODO also check for terminal observation
        a = action(policy, b)
        # println("State: $s    Belief: $(b.b)    Action: $a")
        if s == b.state_list[findmax(b.b)[2]]
            bcorrect += 1
        end
        sp, o, r = generate_sor(pomdp, b.b, s, a, sim.rng)
        verbose && (@show b.b)
        verbose && (@show a)
        verbose && (@show o)
        verbose && (@show r)
        r_total += disc*r
        s = sp
        bp = update(updater, b, a, o)
        b = bp
        disc *= discount(pomdp)
        pcorrect_inner = bcorrect / step
        verbose && (@show pcorrect_inner)
        step += 1
    end
    pcorrect = bcorrect/step
    return r_total, pcorrect
end

function simulate_worst(sim::RolloutSimulator, pomdp::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}, policy::Policy, updater::Updater, alphas::Vector{Vector{Float64}}; verbose = false)
    initial_belief = initial_state_distribution(pomdp)
    s = rand(sim.rng, initial_belief)
    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))
    disc = 1.0
    r_total = 0.0
    b = initialize_belief(updater, initial_belief)
    step = 1
    bcorrect = 0
    while disc > eps && !isterminal(pomdp, s) && step <= max_steps # TODO also check for terminal observation
        a = action(policy, b)
        # println("State: $s    Belief: $(b.b)    Action: $a")
        if s == b.state_list[findmax(b.b)[2]]
            bcorrect += 1
        end
        sp, o, r = generate_sor_worst(pomdp, b.b, s, a, sim.rng, alphas)
        verbose && (@show b.b)
        verbose && (@show a)
        verbose && (@show o)
        verbose && (@show r)
        r_total += disc*r
        s = sp
        bp = update(updater, b, a, o)
        b = bp
        disc *= discount(pomdp)
        step += 1
    end
    pcorrect = bcorrect/step
    return r_total, pcorrect
end

function simulate_worst(sim::RolloutSimulator,
        pomdp::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}, policy::Policy,
        updater::Updater, initial_belief, s,
        alphas::Vector{Vector{Float64}}; verbose = false)
    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))
    disc = 1.0
    r_total = 0.0
    b = initialize_belief(updater, initial_belief)
    step = 1
    bcorrect = 0
    while disc > eps && !isterminal(pomdp, s) && step <= max_steps # TODO also check for terminal observation
        a = action(policy, b)
        if s == b.state_list[findmax(b.b)[2]]
            bcorrect += 1
        end
        sp, o, r = generate_sor_worst(pomdp, b.b, s, a, sim.rng, alphas)
        verbose && (@show b.b)
        verbose && (@show a)
        verbose && (@show o)
        verbose && (@show r)
        r_total += disc*r
        s = sp
        bp = update(updater, b, a, o)
        b = bp
        disc *= discount(pomdp)
        pcorrect_inner = bcorrect / step
        verbose && (@show pcorrect_inner)
        step += 1
    end
    pcorrect = bcorrect/step
    return r_total, pcorrect
end

# worst-case simulator that tracks hellinger distance
function simulate_worst(sim::RolloutSimulator,
        pomdp::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}, policy::Policy,
        updater::Updater, reference_updater::Updater,
        initial_belief, s,
        alphas::Vector{Vector{Float64}}; verbose = false)
    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))
    disc = 1.0
    r_total = 0.0
    b = initialize_belief(updater, initial_belief)
    step = 1
    bcorrect = 0
    while disc > eps && !isterminal(pomdp, s) && step <= max_steps # TODO also check for terminal observation
        a = action(policy, b)
        if s == b.state_list[findmax(b.b)[2]]
            bcorrect += 1
        end
        sp, o, r = generate_sor_worst(pomdp, b.b, s, a, sim.rng, alphas)
        verbose && (@show b.b)
        verbose && (@show a)
        verbose && (@show o)
        verbose && (@show r)
        r_total += disc*r
        s = sp
        bp = update(updater, b, a, o)
        bp_reference = update(reference_updater, b, a, o)
        bdistance = hellinger(bp.b, bp_reference.b)
        verbose && (@show bdistance)
        bdistance_total += disc * bdistance
        pcorrect_inner = bcorrect / step
        verbose && (@show pcorrect_inner)
        b = bp
        disc *= discount(pomdp)
        step += 1
    end
    pcorrect = bcorrect/step
    return r_total, pcorrect
end

function simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy)
    istate=get(sim.initial_state, initial_state(mdp, sim.rng))
    simulate(sim, mdp, policy, istate)
end

function simulate{S}(sim::RolloutSimulator, mdp::Union{MDP{S}, POMDP{S}}, policy::Policy, initial_state::S)

    eps = get(sim.eps, 0.0)
    max_steps = get(sim.max_steps, typemax(Int))

    s = initial_state

    disc = 1.0
    r_total = 0.0
    step = 1

    while disc > eps && !isterminal(mdp, s) && step <= max_steps
        a = action(policy, s)

        sp, r = generate_sr(mdp, s, a, sim.rng)

        r_total += disc*r

        s = sp

        disc *= discount(mdp)
        step += 1
    end

    return r_total
end


text:
using Revise
using DifferentialEquations
using DynamicalSystems
using Basins
using Plots
using Printf

# Equations of motion:
function forced_pendulum(u, p, t)
    @inbounds begin
    d = p[1]; F = p[2]; omega = p[3]
    du1 = u[2]
    du2 = -d*u[2] - sin(u[1])+ F*cos(omega*t)
    return SVector{2}(du1, du2)
    end
end

# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
    if integrator.u[1] < 0
        integrator.u[1] += 2*π
    else
        integrator.u[1] -= 2*π
    end
end

condition(u,t,integrator) = (integrator.u[1] < -π  || integrator.u[1] > π)

cb = DiscreteCallback(condition,affect!)

#d, F ,w
F = 1.66
ω = 1.
d=0.2
p=[d, F, ω]
#p=[0.15, 0.2, 0.1]
df = ODEProblem(forced_pendulum,rand(2),(0.0,20.0), p)
integ_df  = init(df, alg=AutoTsit5(Rosenbrock23()); reltol=1e-9, save_everystep=false, callback=cb)

xres=200
yres=200

# range for forced pend
xg = range(-pi,pi,length=xres)
yg = range(-2.,4.,length=yres)

# compute basin
@time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)

# Basin entropy
@show Sb,Sbb = basin_entropy(bsn; eps_x=20, eps_y=20)

# Wada merge Haussdorff distances
@time max_dist,min_dist = detect_wada_merge_method(xg, yg, bsn)
epsilon = xg[2]-xg[1]
@show dmax = max_dist/epsilon
@show dmin = min_dist/epsilon

# Wada grid
W = detect_wada_grid_method(integ_df, bsn; max_iter=8)
@show W[:,end]

# Uncertainty exponent for these parameter and grid
bd = box_counting_dim(xg, yg, bsn)
α = 2 - bd

D = uncertainty_exponent(bsn, integ_df)
@show 2-D

println("---------------")
println("---------------")
println("Basin Report: ")
println("---------------")
println("---------------")

@printf("Basin entropy %.2f \n", Sb)
@printf("Boundary Basin Entropy: %.2f\n", Sbb)
@printf("Uncertainty exponent: α= %.2f\n", α )
@printf("Box counting dim: bd= %.2f\n", bd)
@printf("Uncertainty dim estimator: d = %.2f\n", 2-D[1])
@printf("Number of basins: %d\n", bsn.Na)
@printf("Merge Method: Max fattening parameter: %.2f\n", dmax)
@printf("Wada Grid Method: W_Na = %.2f\n ", W[end,end] )

plot(xg, yg, bsn.basin', seriestype=:heatmap)


text:
import Pkg
println(pwd())
Pkg.activate("../.")
# Pkg.add(url="https://github.com/angus-lewis/DiscretisedFluidQueues")
# Pkg.rm("DiscretisedFluidQueues")
Pkg.develop(url=pwd()*"/../../DiscretisedFluidQueues.jl")
Pkg.instantiate(verbose=true)
using DiscretisedFluidQueues, Plots
# plotly()
gr()


text:
load("utils/namespace.jl")

show(scanconsts_let(:(let a=5
    local k::Int
    local l=3
    n = 11
    m::Int = 3

    const x1=4
    const local y2=3
    local const z3=2
    f4(x)=x^2
    function f5(x,y)
        x*y
    end
    abstract U6
    abstract V7<:U6
    type S8; end
    type T9<:S8; end
    type P10{T9}; end
    type Q11{T9} <: P10{T9}; end
    typealias PInt12 P{Int}
    const h13 = x->(const x2=x^2; x2)
    begin
        const c14=4
    end

    let d=34
        const e=6
    end
    g = x->(const x2=x^2; x2)
    for i=1:5
        const j=3
    end

    global const gl=4
    const global lg=5

end)))

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
abstract type ButcherTableauType end

struct BE_1_0_1 <: ButcherTableauType end
struct SDIRK_2_1_2 <: ButcherTableauType end
struct TRBDF2_3_3_2 <: ButcherTableauType end

"""
Butcher tableau
"""
struct ButcherTableau{T <: ButcherTableauType}
  s::Int # stages
  p::Int # embedded order
  q::Int # order
  a::Matrix # A_ij
  b::Vector # b_j
  c::Vector # c_i
  d::Vector # d_j (embedded)
end

"""
ButcherTableau constructor
"""
function ButcherTableau(::BE_1_0_1)
  s = 1
  p = 0
  q = 1
  a = reshape([1.0],1,1)
  b = [1.0]
  c = [1.0]
  d = [0.0]
  ButcherTableau{BE_1_0_1}(s,p,q,a,b,c,d)
end

function ButcherTableau(type::SDIRK_2_1_2)
s = 2
p = 1
q = 2
a = [1.0 0.0; -1.0 1.0]
b = [0.5, 0.5]
c = [1.0, 0.0]
d = [1.0, 0.0]
ButcherTableau{SDIRK_2_1_2}(s,p,q,a,b,c,d)
end

function ButcherTableau(type::TRBDF2_3_3_2)
  s = 3
  p = 3
  q = 2
  aux = 2.0-√2.0
  a = [0.0 0.0 0.0; aux/2 aux/2 0.0; √2/4 √2/4 aux/2]
  b = [√2/4, √2/4, aux/2]
  c = [0.0, aux, 1.0]
  d = [(1.0-(√2/4))/3, ((3*√2)/4+1.0)/3, aux/6]
  ButcherTableau{TRBDF2_3_3_2}(s,p,q,a,b,c,d)
end

function ButcherTableau(type::Symbol)
  eval(:(ButcherTableau($type())))
end

"""
Runge-Kutta ODE solver
"""
struct RungeKutta <: ODESolver
  nls::NonlinearSolver
  dt::Float64
  bt::ButcherTableau
  function RungeKutta(nls,dt,type::Symbol)
    bt = ButcherTableau(type)
    new(nls,dt,bt)
  end
end

function solve_step!(uf::AbstractVector,
  solver::RungeKutta,
  op::ODEOperator,
  u0::AbstractVector,
  t0::Real,
  cache)

  # Unpack variables
  dt = solver.dt
  s = solver.bt.s
  a = solver.bt.a
  b = solver.bt.b
  c = solver.bt.c
  d = solver.bt.d

  # Create cache if not there
  if cache === nothing
    ode_cache = allocate_cache(op)
    vi = similar(u0)
    fi = [similar(u0)]
    nl_cache = nothing
  else
    ode_cache, vi, fi, nl_cache = cache
  end

  # Create RKNL operator
  nlop = RungeKuttaNonlinearOperator(op,t0,dt,u0,ode_cache,vi,fi,0,a)

  # Compute intermediate stages
  for i in 1:s

    # allocate space to store the RHS at i
    if (length(fi) < i)
      push!(fi,similar(u0))
    end

    # Skip stage solve if a_ii=0 => u_i=u_0, f_i = f_0
    if(a[i,i]==0)
      @assert c[i] == 0
      ti = t0
      update!(nlop,ti,fi,i)
      fi[i] = get_fi(u0,nlop,nl_cache)
    else
      # solve at stage i
      ti = t0 + c[i]*dt
      ode_cache = update_cache!(ode_cache,op,ti)
      update!(nlop,ti,fi,i)
      nl_cache = solve!(uf,solver.nls,nlop,nl_cache)
      fi[i] = get_fi(uf,nlop,nl_cache)
    end

  end

  # update
  uf = u0
  for i in 1:s
    uf = uf + dt*b[i]*fi[i]
  end

  cache = (ode_cache, vi, fi, nl_cache)

  tf = t0 + dt
  return (uf,tf,cache)

end

"""
Nonlinear operator that represents the Runge-Kutta nonlinear operator at a
given time step and stage, i.e., A(t,u_i,(u_i-u_n)/dt)
"""
mutable struct RungeKuttaNonlinearOperator <: NonlinearOperator
  odeop::ODEOperator
  ti::Float64
  dt::Float64
  u0::AbstractVector
  ode_cache
  vi::AbstractVector
  fi::AbstractVector
  i::Int
  a::Matrix
end

function residual!(b::AbstractVector,op::RungeKuttaNonlinearOperator,x::AbstractVector)
  # A(t,ui,∂ui/∂t) = ∂ui/∂t - a_ii * f(ui,ti) - ∑_{j<i} a_ij * f(uj,tj) = 0
  # b = [1/a_ii * ∂u/∂t - f(ui,ti)]
  # Res_ij = - a_ij/a_ii * f(uj,ti)
  # b + ∑_{j<i} Res_ij = 0
  @assert (abs(op.a[op.i,op.i]) > 0.0)
  ui = x
  vi = op.vi
  vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt)
  residual!(b,op.odeop,op.ti,(ui,vi),op.ode_cache)
  for j in 1:op.i-1
    b .= b - op.a[op.i,j]/op.a[op.i,op.i] * op.fi[j]
  end
  b
end

function jacobian!(A::AbstractMatrix,op::RungeKuttaNonlinearOperator,x::AbstractVector)
  @assert (abs(op.a[op.i,op.i]) > 0.0)
  ui = x
  vi = op.vi
  vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt)
  z = zero(eltype(A))
  fillstored!(A,z)
  jacobians!(A,op.odeop,op.ti,(ui,vi),(1/(op.a[op.i,op.i]*op.dt)),op.ode_cache)
end

function allocate_residual(op::RungeKuttaNonlinearOperator,x::AbstractVector)
  allocate_residual(op.odeop,x,op.ode_cache)
end

function allocate_jacobian(op::RungeKuttaNonlinearOperator,x::AbstractVector)
  allocate_jacobian(op.odeop,x,op.ode_cache)
end

function zero_initial_guess(op::RungeKuttaNonlinearOperator)
  x0 = similar(op.u0)
  fill!(x0,zero(eltype(x0)))
  x0
end

function get_fi(x::AbstractVector, op::RungeKuttaNonlinearOperator, cache::Nothing)
  ui = x
  vi = op.vi
  if(op.a[op.i,op.i]==0.0)
    vi=zero(x)
  else
    vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt)
  end
  b=similar(x)
  residual!(b,op.odeop,op.ti,(ui,vi),op.ode_cache)
  (vi-b) # store fi for future stages
end
function get_fi(x::AbstractVector, op::RungeKuttaNonlinearOperator, cache)
  ui = x
  vi = op.vi
  if(op.a[op.i,op.i]==0.0)
    vi=zero(x)
  else
    vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt)
  end
  residual!(cache.b,op.odeop,op.ti,(ui,vi),op.ode_cache)
  (vi-cache.b) # store fi for future stages
end

function update!(op::RungeKuttaNonlinearOperator,ti::Float64,fi::AbstractVector,i::Int)
  op.ti = ti
  op.fi = fi
  op.i = i
end


text:
const BosonStateType = Vector{Int}

struct BosonHilbertSpace{LatticeType<:AbstractSiteNetwork,IndexType<:AbstractUniqueVector{BosonStateType}} <: HilbertSpace{BosonStateType}
    lattice::LatticeType
    indexer::IndexType
end

BosonHilbertSpace(lattice) =
    BosonHilbertSpace(lattice, UniqueVector{BosonStateType}())

statetype(::BosonHilbertSpace) = BosonStateType

get_charge(::BosonHilbertSpace, site_state::Integer) = site_state
get_total_charge(::BosonHilbertSpace, state::BosonStateType) =
    sum(state)
get_total_charge(hs::BosonHilbertSpace, stateidx::Int) = get_total_charge(hs, hs.indexer[stateidx])

struct BosonParameters
    J::Float64
    U::Float64
    θ::Float64

    function BosonParameters(;
                             J::Real=0.0,
                             U::Real=0.0,
                             θ::Real=0.0)
        new(J, U, θ)
    end
end

boson_hamiltonian(; kwargs...) =
    boson_hamiltonian(BosonParameters(; kwargs...))

function boson_hamiltonian(p::BosonParameters)
    return function apply_hamiltonian(f, hs::BosonHilbertSpace, s_j::Integer)
        state = hs.indexer[s_j]

        # Hopping
        neighborsη(hs.lattice) do x::Int, x_r::Int, η::Rational{Int}
            # FIXME: what to do with η ??
            # FIXME: Also, we are assuming 1D.
            # XXX FIXME: type stability; currently this always returns a Complex number!

            # Rightward hop
            if state[x] > 0
                other = copy(state)
                phase = exp(im * p.θ * other[x_r])
                factor = √other[x]
                other[x] -= 1
                other[x_r] += 1
                factor *= √other[x_r]
                s_i = findfirst!(isequal(other), hs.indexer)
                f(s_i, -p.J * factor * phase)
            end
            # Leftward hop
            if state[x_r] > 0
                other = copy(state)
                factor = √other[x_r]
                other[x_r] -= 1
                other[x] += 1
                factor *= √other[x]
                phase = exp(-im * p.θ * other[x_r])
                s_i = findfirst!(isequal(other), hs.indexer)
                f(s_i, -p.J * factor * phase)
            end
        end

        # Bose-Hubbard U
        if !iszero(p.U)
            cnt_U = 0
            for x in 1:length(hs.lattice)
                cnt_U += state[x] * (state[x] - 1)
            end
            f(s_j, 0.5 * p.U * cnt_U)
        end

        nothing
    end
end

function apply_n(f, hs::BosonHilbertSpace, j::Integer, x1::Integer)
    state = hs.indexer[j]
    f(j, state[x1])
    nothing
end
export apply_n

function seed_state!(hs::BosonHilbertSpace, N::Integer)
    state = zeros(Int, length(hs.lattice))
    state[1] = N # XXX FIXME
    findfirst!(isequal(state), hs.indexer)
    hs
end

function reflectionη(hs::BosonHilbertSpace, j::Integer)
    state = hs.indexer[j]
    i = findfirst!(isequal(reverse(state)), hs.indexer)
    return i, 0//1
end


text:
int_rules_8_9 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.9*Product*logarithm*function =#
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') + ~(b') * ~x)) ^ ~p, ~x) => ((~a + ~b * ~x) * (~c * ProductLog(~a + ~b * ~x)) ^ ~p) / (~b * (~p + 1)) + (~p / (~c * (~p + 1))) * Antiderivative((~c * ProductLog(~a + ~b * ~x)) ^ (~p + 1) / (1 + ProductLog(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && LtQ(~p, -1)
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') + ~(b') * ~x)) ^ ~(p'), ~x) => ((~a + ~b * ~x) * (~c * ProductLog(~a + ~b * ~x)) ^ ~p) / ~b - ~p * Antiderivative((~c * ProductLog(~a + ~b * ~x)) ^ ~p / (1 + ProductLog(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && Not(LtQ(~p, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(c') * ProductLog(~a + ~(b') * ~x)) ^ ~(p'), ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand((~c * ProductLog(~x)) ^ ~p, ((~b * ~e - ~a * ~f) + ~f * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~p], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~n)) ^ ~(p'), ~x) => ~x * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p - ~n * ~p * Antiderivative((~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p / (1 + ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~n, ~p], ~x) && (EqQ(~n * (~p - 1), -1) || IntegerQ(~p - 1 / 2) && EqQ(~n * (~p - 1 / 2), -1))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~n)) ^ ~(p'), ~x) => (~x * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p) / (~n * ~p + 1) + ((~n * ~p) / (~c * (~n * ~p + 1))) * Antiderivative((~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p + 1) / (1 + ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~n], ~x) && (IntegerQ(~p) && EqQ(~n * (~p + 1), -1) || IntegerQ(~p - 1 / 2) && EqQ(~n * (~p + 1 / 2), -1))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~c * ProductLog(~a * (~x) ^ -(~n))) ^ ~p / (~x) ^ 2, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~c, ~p], ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p'), ~x) => ((~x) ^ (~m + 1) * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p) / (~m + 1) - ((~n * ~p) / (~m + 1)) * Antiderivative(((~x) ^ ~m * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p) / (1 + ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~m, ~n, ~p], ~x) && (NeQ(~m, -1) && (IntegerQ(~p - 1 / 2) && IGtQ(2 * Simplify(~p + (~m + 1) / ~n), 0) || Not(IntegerQ(~p - 1 / 2)) && IGtQ(Simplify(~p + (~m + 1) / ~n) + 1, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p'), ~x) => ((~x) ^ (~m + 1) * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p) / (~m + ~n * ~p + 1) + ((~n * ~p) / (~c * (~m + ~n * ~p + 1))) * Antiderivative(((~x) ^ ~m * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p + 1)) / (1 + ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~m, ~n, ~p], ~x) && (EqQ(~m, -1) || (IntegerQ(~p - 1 / 2) && ILtQ(Simplify(~p + (~m + 1) / ~n) - 1 / 2, 0) || Not(IntegerQ(~p - 1 / 2)) && ILtQ(Simplify(~p + (~m + 1) / ~n), 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * ~x)) ^ ~(p'), ~x) => Antiderivative(((~x) ^ ~m * (~c * ProductLog(~a * ~x)) ^ ~p) / (1 + ProductLog(~a * ~x)), ~x) + (1 / ~c) * Antiderivative(((~x) ^ ~m * (~c * ProductLog(~a * ~x)) ^ (~p + 1)) / (1 + ProductLog(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~m], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~c * ProductLog(~a * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~c, ~p], ~x) && (ILtQ(~n, 0) && (IntegerQ(~m) && NeQ(~m, -1)))
    @apply_utils Antiderivative(1 / (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x)), ~x) => (~a + ~b * ~x) / (~b * ~d * ProductLog(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(ProductLog(~(a') + ~(b') * ~x) / (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x)), ~x) => ~d * ~x - Antiderivative(1 / (~d + ~d * ProductLog(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') + ~(b') * ~x)) ^ ~p / (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x)), ~x) => (~c * (~a + ~b * ~x) * (~c * ProductLog(~a + ~b * ~x)) ^ (~p - 1)) / (~b * ~d) - ~c * ~p * Antiderivative((~c * ProductLog(~a + ~b * ~x)) ^ (~p - 1) / (~d + ~d * ProductLog(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && GtQ(~p, 0)
    @apply_utils Antiderivative(1 / (ProductLog(~(a') + ~(b') * ~x) * (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x))), ~x) => ExpIntegralEi(ProductLog(~a + ~b * ~x)) / (~b * ~d) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / (sqrt(~(c') * ProductLog(~(a') + ~(b') * ~x)) * (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x))), ~x) => (Rt(pi * ~c, 2) * erfi(sqrt(~c * ProductLog(~a + ~b * ~x)) / Rt(~c, 2))) / (~b * ~c * ~d) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && PosQ(~c)
    @apply_utils Antiderivative(1 / (sqrt(~(c') * ProductLog(~(a') + ~(b') * ~x)) * (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x))), ~x) => (Rt(-pi * ~c, 2) * erf(sqrt(~c * ProductLog(~a + ~b * ~x)) / Rt(-(~c), 2))) / (~b * ~c * ~d) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NegQ(~c)
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') + ~(b') * ~x)) ^ ~p / (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x)), ~x) => ((~a + ~b * ~x) * (~c * ProductLog(~a + ~b * ~x)) ^ ~p) / (~b * ~d * (~p + 1)) - (1 / (~c * (~p + 1))) * Antiderivative((~c * ProductLog(~a + ~b * ~x)) ^ (~p + 1) / (~d + ~d * ProductLog(~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && LtQ(~p, -1)
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') + ~(b') * ~x)) ^ ~(p') / (~d + ~(d') * ProductLog(~(a') + ~(b') * ~x)), ~x) => (Gamma(~p + 1, -(ProductLog(~a + ~b * ~x))) * (~c * ProductLog(~a + ~b * ~x)) ^ ~p) / (~b * ~d * (-(ProductLog(~a + ~b * ~x))) ^ ~p) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') / (~d + ~(d') * ProductLog(~a + ~(b') * ~x)), ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand(1 / (~d + ~d * ProductLog(~x)), ((~b * ~e - ~a * ~f) + ~f * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(((~(e') + ~(f') * ~x) ^ ~(m') * (~(c') * ProductLog(~a + ~(b') * ~x)) ^ ~(p')) / (~d + ~(d') * ProductLog(~a + ~(b') * ~x)), ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand((~c * ProductLog(~x)) ^ ~p / (~d + ~d * ProductLog(~x)), ((~b * ~e - ~a * ~f) + ~f * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(1 / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~n)), ~x) => -(Subst(Antiderivative(1 / ((~x) ^ 2 * (~d + ~d * ProductLog(~a * (~x) ^ -(~n)))), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~d], ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p') / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (~c * ~x * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p - 1)) / ~d <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && EqQ(~n * (~p - 1), -1)
    @apply_utils Antiderivative(ProductLog(~(a') * (~x) ^ ~(n')) ^ ~(p') / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => ((~a) ^ ~p * ExpIntegralEi(-(~p) * ProductLog(~a * (~x) ^ ~n))) / (~d * ~n) <-- FreeQ([~a, ~d], ~x) && (IntegerQ(~p) && EqQ(~n * ~p, -1))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~p / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (Rt(pi * ~c * ~n, 2) / (~d * ~n * (~a) ^ (1 / ~n) * (~c) ^ (1 / ~n))) * erfi(sqrt(~c * ProductLog(~a * (~x) ^ ~n)) / Rt(~c * ~n, 2)) <-- FreeQ([~a, ~c, ~d], ~x) && (IntegerQ(1 / ~n) && (EqQ(~p, 1 / 2 - 1 / ~n) && PosQ(~c * ~n)))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~p / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (Rt(-pi * ~c * ~n, 2) / (~d * ~n * (~a) ^ (1 / ~n) * (~c) ^ (1 / ~n))) * erf(sqrt(~c * ProductLog(~a * (~x) ^ ~n)) / Rt(-(~c) * ~n, 2)) <-- FreeQ([~a, ~c, ~d], ~x) && (IntegerQ(1 / ~n) && (EqQ(~p, 1 / 2 - 1 / ~n) && NegQ(~c * ~n)))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p') / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (~c * ~x * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p - 1)) / ~d - ~c * (~n * (~p - 1) + 1) * Antiderivative((~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p - 1) / (~d + ~d * ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~d], ~x) && (GtQ(~n, 0) && GtQ(~n * (~p - 1) + 1, 0))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p') / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (~x * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p) / (~d * (~n * ~p + 1)) - (1 / (~c * (~n * ~p + 1))) * Antiderivative((~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p + 1) / (~d + ~d * ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~d], ~x) && (GtQ(~n, 0) && LtQ(~n * ~p + 1, 0))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~n)) ^ ~(p') / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~n)), ~x) => -(Subst(Antiderivative((~c * ProductLog(~a * (~x) ^ -(~n))) ^ ~p / ((~x) ^ 2 * (~d + ~d * ProductLog(~a * (~x) ^ -(~n)))), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && ILtQ(~n, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') / (~d + ~(d') * ProductLog(~(a') * ~x)), ~x) => (~x) ^ (~m + 1) / (~d * (~m + 1) * ProductLog(~a * ~x)) - (~m / (~m + 1)) * Antiderivative((~x) ^ ~m / (ProductLog(~a * ~x) * (~d + ~d * ProductLog(~a * ~x))), ~x) <-- FreeQ([~a, ~d], ~x) && GtQ(~m, 0)
    @apply_utils Antiderivative(1 / (~x * (~d + ~(d') * ProductLog(~(a') * ~x))), ~x) => log(ProductLog(~a * ~x)) / ~d <-- FreeQ([~a, ~d], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') / (~d + ~(d') * ProductLog(~(a') * ~x)), ~x) => (~x) ^ (~m + 1) / (~d * (~m + 1)) - Antiderivative(((~x) ^ ~m * ProductLog(~a * ~x)) / (~d + ~d * ProductLog(~a * ~x)), ~x) <-- FreeQ([~a, ~d], ~x) && LtQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~(m') / (~d + ~(d') * ProductLog(~(a') * ~x)), ~x) => ((~x) ^ ~m * Gamma(~m + 1, -((~m + 1)) * ProductLog(~a * ~x))) / (~a * ~d * (~m + 1) * ℯ ^ (~m * ProductLog(~a * ~x)) * (-((~m + 1)) * ProductLog(~a * ~x)) ^ ~m) <-- FreeQ([~a, ~d, ~m], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(1 / (~x * (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n')))), ~x) => log(ProductLog(~a * (~x) ^ ~n)) / (~d * ~n) <-- FreeQ([~a, ~d, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~n)), ~x) => -(Subst(Antiderivative(1 / ((~x) ^ (~m + 2) * (~d + ~d * ProductLog(~a * (~x) ^ -(~n)))), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~d], ~x) && (IntegerQ(~m) && (ILtQ(~n, 0) && NeQ(~m, -1)))
    @apply_utils Antiderivative((~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p') / (~x * (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n')))), ~x) => (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p / (~d * ~n * ~p) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p')) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (~c * (~x) ^ (~m + 1) * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p - 1)) / (~d * (~m + 1)) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (NeQ(~m, -1) && EqQ(~m + ~n * (~p - 1), -1))
    @apply_utils Antiderivative(((~x) ^ ~(m') * ProductLog(~(a') * (~x) ^ ~(n')) ^ ~(p')) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => ((~a) ^ ~p * ExpIntegralEi(-(~p) * ProductLog(~a * (~x) ^ ~n))) / (~d * ~n) <-- FreeQ([~a, ~d, ~m, ~n], ~x) && (IntegerQ(~p) && EqQ(~m + ~n * ~p, -1))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~p) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => ((~a) ^ (~p - 1 / 2) * (~c) ^ (~p - 1 / 2) * Rt((pi * ~c) / (~p - 1 / 2), 2) * erf(sqrt(~c * ProductLog(~a * (~x) ^ ~n)) / Rt(~c / (~p - 1 / 2), 2))) / (~d * ~n) <-- FreeQ([~a, ~c, ~d, ~m, ~n], ~x) && (NeQ(~m, -1) && (IntegerQ(~p - 1 / 2) && (EqQ(~m + ~n * (~p - 1 / 2), -1) && PosQ(~c / (~p - 1 / 2)))))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~p) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => ((~a) ^ (~p - 1 / 2) * (~c) ^ (~p - 1 / 2) * Rt((-pi * ~c) / (~p - 1 / 2), 2) * erfi(sqrt(~c * ProductLog(~a * (~x) ^ ~n)) / Rt(-(~c) / (~p - 1 / 2), 2))) / (~d * ~n) <-- FreeQ([~a, ~c, ~d, ~m, ~n], ~x) && (NeQ(~m, -1) && (IntegerQ(~p - 1 / 2) && (EqQ(~m + ~n * (~p - 1 / 2), -1) && NegQ(~c / (~p - 1 / 2)))))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p')) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => (~c * (~x) ^ (~m + 1) * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p - 1)) / (~d * (~m + 1)) - ((~c * (~m + ~n * (~p - 1) + 1)) / (~m + 1)) * Antiderivative(((~x) ^ ~m * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p - 1)) / (~d + ~d * ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (NeQ(~m, -1) && GtQ(Simplify(~p + (~m + 1) / ~n), 1))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p')) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => ((~x) ^ (~m + 1) * (~c * ProductLog(~a * (~x) ^ ~n)) ^ ~p) / (~d * (~m + ~n * ~p + 1)) - ((~m + 1) / (~c * (~m + ~n * ~p + 1))) * Antiderivative(((~x) ^ ~m * (~c * ProductLog(~a * (~x) ^ ~n)) ^ (~p + 1)) / (~d + ~d * ProductLog(~a * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (NeQ(~m, -1) && LtQ(Simplify(~p + (~m + 1) / ~n), 0))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * ~x)) ^ ~(p')) / (~d + ~(d') * ProductLog(~(a') * ~x)), ~x) => ((~x) ^ ~m * Gamma(~m + ~p + 1, -((~m + 1)) * ProductLog(~a * ~x)) * (~c * ProductLog(~a * ~x)) ^ ~p) / (~a * ~d * (~m + 1) * ℯ ^ (~m * ProductLog(~a * ~x)) * (-((~m + 1)) * ProductLog(~a * ~x)) ^ (~m + ~p)) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~(c') * ProductLog(~(a') * (~x) ^ ~(n'))) ^ ~(p')) / (~d + ~(d') * ProductLog(~(a') * (~x) ^ ~(n'))), ~x) => -(Subst(Antiderivative((~c * ProductLog(~a * (~x) ^ -(~n))) ^ ~p / ((~x) ^ (~m + 2) * (~d + ~d * ProductLog(~a * (~x) ^ -(~n)))), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (NeQ(~m, -1) && (IntegerQ(~m) && LtQ(~n, 0)))
    @apply_utils Antiderivative(~u, ~x) => Subst(Antiderivative(SimplifyIntegrand((~x + 1) * ℯ ^ ~x * SubstFor(ProductLog(~x), ~u, ~x), ~x), ~x), ~x, ProductLog(~x)) <-- FunctionOfQ(ProductLog(~x), ~u, ~x)
end

max_stars_repo_path: src/rules/8 Special functions/8.9 Product logarithm function.jl


text:
@testset "mean" begin
    n = 9
    @testset "Basic" begin
        test_rrule(mean, randn(n))
    end
    @testset "with dims kwargs" begin
        test_rrule(mean, randn(n); fkwargs=(;dims=1))
        test_rrule(mean, randn(n,4); fkwargs=(;dims=2))
    end
end


text:
using Test
using PlotAdaptiveGrid1D

@testset "adapted_grid" begin
  @testset "sin(1/x) on [0,1]" begin
    f(x) = sin(1/x)
    xs = adapted_grid(f, (0.0, 1.0))
    @test xs == [7.629546109209478e-5, 0.00015259092218418955, 0.00022888638327628433, 0.0003051818443683791, 0.0003814773054604739, 0.00045777276655256866, 0.0005340682276446635, 0.0006103636887367582, 0.0006866591498288529, 0.0007629546109209478, 0.0008392500720130426, 0.0009155455331051373, 0.0009918409941972323, 0.001068136455289327, 0.0011444319163814217, 0.0012207273774735164, 0.0012970228385656112, 0.0013733182996577059, 0.0014496137607498006, 0.0015259092218418955, 0.0016022046829339905, 0.0016785001440260852, 0.00175479560511818, 0.0018310910662102746, 0.0019073865273023696, 0.0019836819883944645, 0.0020599774494865592, 0.002136272910578654, 0.0022125683716707487, 0.0022888638327628434, 0.002365159293854938, 0.002441454754947033, 0.0025177502160391276, 0.0025940456771312223, 0.002670341138223317, 0.0027466365993154118, 0.0028229320604075065, 0.002899227521499601, 0.002975522982591696, 0.003051818443683791, 0.0031281139047758862, 0.003204409365867981, 0.0032807048269600757, 0.0033570002880521704, 0.003433295749144265, 0.00350959121023636, 0.0035858866713284546, 0.0036621821324205493, 0.003738477593512644, 0.003814773054604739, 0.0038910685156968343, 0.003967363976788929, 0.004043659437881024, 0.0041199548989731185, 0.004196250360065213, 0.004272545821157308, 0.004348841282249403, 0.004425136743341497, 0.004501432204433592, 0.004577727665525687, 0.0046540231266177816, 0.004730318587709876, 0.004806614048801971, 0.004882909509894066, 0.0049592049709861604, 0.005035500432078255, 0.00511179589317035, 0.005188091354262445, 0.005264386815354539, 0.005340682276446634, 0.005416977737538729, 0.0054932731986308235, 0.005569568659722918, 0.005645864120815013, 0.005722159581907108, 0.005798455042999202, 0.005874750504091297, 0.005951045965183392, 0.006027341426275487, 0.006103636887367582, 0.006179932348459677, 0.0062562278095517725, 0.006332523270643867, 0.006408818731735962, 0.006485114192828057, 0.006561409653920151, 0.006637705115012246, 0.006714000576104341, 0.0067902960371964355, 0.00686659149828853, 0.006942886959380625, 0.00701918242047272, 0.007095477881564814, 0.007171773342656909, 0.007248068803749004, 0.007324364264841099, 0.007400659725933193, 0.007476955187025288, 0.007553250648117384, 0.007629546109209478, 0.007705841570301573, 0.007782137031393669, 0.007858432492485763, 0.007934727953577858, 0.008011023414669954, 0.008087318875762048, 0.008163614336854141, 0.008239909797946237, 0.008316205259038333, 0.008392500720130426, 0.00846879618122252, 0.008545091642314616, 0.008621387103406711, 0.008697682564498805, 0.0087739780255909, 0.008850273486682995, 0.00892656894777509, 0.009002864408867184, 0.009079159869959278, 0.009155455331051374, 0.00923175079214347, 0.009308046253235563, 0.009384341714327657, 0.009460637175419753, 0.009536932636511848, 0.009613228097603942, 0.009689523558696036, 0.009765819019788131, 0.01047128072407395, 0.011176742428359768, 0.011882204132645586, 0.012587665836931404, 0.013293127541217222, 0.01399858924550304, 0.014704050949788858, 0.015409512654074675, 0.01611497435836049, 0.01682043606264631, 0.017525897766932128, 0.018231359471217946, 0.018936821175503764, 0.019642282879789582, 0.0203477445840754, 0.021053206288361218, 0.021758667992647036, 0.02246412969693285, 0.02316959140121867, 0.023875053105504487, 0.024580514809790305, 0.025285976514076124, 0.025991438218361942, 0.02669689992264776, 0.027402361626933578, 0.028107823331219396, 0.028813285035505214, 0.029518746739791032, 0.03022420844407685, 0.03092967014836267, 0.03163513185264849, 0.032340593556934305, 0.03304605526122012, 0.03375151696550594, 0.03445697866979176, 0.03516244037407758, 0.03586790207836339, 0.036573363782649214, 0.037278825486935036, 0.03798428719122085, 0.038689748895506665, 0.039395210599792486, 0.04010067230407831, 0.04080613400836412, 0.04151159571264994, 0.04221705741693576, 0.04292251912122158, 0.043627980825507395, 0.04433344252979321, 0.04503890423407903, 0.04574436593836485, 0.04644982764265067, 0.04715528934693648, 0.047860751051222304, 0.048566212755508126, 0.04927167445979394, 0.049977136164079755, 0.05068259786836558, 0.0513880595726514, 0.05209352127693721, 0.05279898298122303, 0.05350444468550885, 0.05420990638979467, 0.054915368094080486, 0.0556208297983663, 0.05632629150265212, 0.057031753206937943, 0.05773721491122376, 0.05844267661550957, 0.059148138319795394, 0.059853600024081216, 0.06055906172836703, 0.061264523432652845, 0.06196998513693867, 0.06267544684122449, 0.0633809085455103, 0.06408637024979613, 0.06479183195408195, 0.06549729365836776, 0.06620275536265358, 0.06690821706693939, 0.0676136787712252, 0.06831914047551102, 0.06902460217979685, 0.06973006388408268, 0.07043552558836849, 0.0711409872926543, 0.07184644899694012, 0.07255191070122594, 0.07325737240551175, 0.07396283410979757, 0.0746682958140834, 0.07537375751836922, 0.07607921922265504, 0.07678468092694085, 0.07749014263122667, 0.0789010660397983, 0.08031198944836994, 0.08172291285694158, 0.08313383626551321, 0.08383929796979903, 0.08454475967408484, 0.08525022137837066, 0.08595568308265648, 0.08666114478694231, 0.08736660649122813, 0.08807206819551394, 0.08877752989979976, 0.08948299160408557, 0.09018845330837139, 0.0908939150126572, 0.09159937671694303, 0.09301030012551467, 0.0944212235340863, 0.09583214694265793, 0.09724307035122956, 0.09865399375980119, 0.10006491716837283, 0.10176072948516193, 0.10345654180195102, 0.10515235411874012, 0.10684816643552922, 0.1085439787523183, 0.1102397910691074, 0.1119356033858965, 0.11363141570268559, 0.11532722801947469, 0.11702304033626379, 0.11871885265305288, 0.12041466496984198, 0.12211047728663108, 0.12380628960342016, 0.12550210192020927, 0.12719791423699836, 0.12889372655378745, 0.13058953887057656, 0.13228535118736565, 0.13398116350415473, 0.13737278813773293, 0.14076441277131113, 0.1441560374048893, 0.1475476620384675, 0.1509392866720457, 0.15433091130562387, 0.16006606223189163, 0.16580121315815935, 0.1686687886212932, 0.17153636408442707, 0.17440393954756095, 0.17727151501069482, 0.1801390904738287, 0.18300666593696255, 0.1858742414000964, 0.18874181686323027, 0.19447696778949802, 0.20021211871576575, 0.20644210013042646, 0.21267208154508715, 0.21890206295974787, 0.2251320443744086, 0.23136202578906928, 0.23759200720373, 0.24382198861839072, 0.2500519700330514, 0.2730259372838398, 0.29599990453462816, 0.3083151650025936, 0.320630425470559, 0.33294568593852436, 0.3452609464064898, 0.35952938338950047, 0.3737978203725111, 0.3880662573555217, 0.4023346943385323, 0.42769490022477685, 0.45305510611102146, 0.47529372132562075, 0.49753233654022, 0.5217783965977698, 0.5460244566553195, 0.5704146177040659, 0.5948047787528121, 0.6487376626244042, 0.7000534407448434, 0.7534099435113505, 0.7983247536668452, 0.8531606459863641, 0.9047080021287085, 0.9872908640160599]
  end

  @testset "1/x on [0, 1]" begin
    f(x) = 1/x
    xs = adapted_grid(f, (0.0, 1.0))
    @test xs == [7.629546109209478e-5, 0.00015259092218418955, 0.00022888638327628433, 0.0003051818443683791, 0.0003814773054604739, 0.00045777276655256866, 0.0005340682276446635, 0.0006103636887367582, 0.0006866591498288529, 0.0007629546109209478, 0.0008392500720130426, 0.0009155455331051373, 0.001068136455289327, 0.0012207273774735164, 0.0013733182996577059, 0.0015259092218418955, 0.0016785001440260852, 0.0018310910662102746, 0.002136272910578654, 0.002441454754947033, 0.0027466365993154118, 0.003051818443683791, 0.0033570002880521704, 0.0036621821324205493, 0.004272545821157308, 0.004882909509894066, 0.0054932731986308235, 0.006103636887367582, 0.006714000576104341, 0.007324364264841099, 0.007934727953577858, 0.008545091642314616, 0.009155455331051374, 0.009765819019788131, 0.015409512654074675, 0.021053206288361218, 0.02669689992264776, 0.032340593556934305, 0.043627980825507395, 0.054915368094080486, 0.07749014263122667, 0.10006491716837283, 0.12719791423699836, 0.15433091130562387, 0.17727151501069482, 0.20021211871576575, 0.2500519700330514, 0.29599990453462816, 0.3452609464064898, 0.4023346943385323, 0.45305510611102146, 0.49753233654022, 0.5460244566553195, 0.5948047787528121, 0.6487376626244042, 0.7000534407448434, 0.7534099435113505, 0.7983247536668452, 0.8531606459863641, 0.9047080021287085, 0.9872908640160599]
  end

  @testset "abs(x) on [-1, 1]" begin
    xs = adapted_grid(abs, (-1.0, 1.0))
    @test xs == [-0.9804683619604238, -0.7998701656632544, -0.6913381773887523, -0.5995757625684686, -0.49989605993389713, -0.4080001909307437, -0.3094781071870203, -0.19533061132293544, -0.1446101995504463, -0.09388978777795712, -0.07165117256335783, -0.049412557348758546, -0.027173942134159263, -0.004935326919559979, 0.007187703109214885, 0.019310733137989748, 0.03143376316676461, 0.04355679319553948, 0.0678028532530892, 0.09204891331063893, 0.14082923540813155, 0.1896095575056242, 0.29747532524880854, 0.4001068814896869, 0.5068198870227011, 0.5966495073336903, 0.7063212919727282, 0.809416004257417, 0.9745817280321198]
  end
end


text:
# Doc from https://github.com/spglib/spglib/blob/d1cb3bd/src/spglib.h#L424-L439
"""
    get_ir_reciprocal_mesh(cell::Cell, mesh, is_shift=falses(3); is_time_reversal=true, symprec=1e-5)

Return k-points mesh and k-point map to the irreducible k-points.

Irreducible reciprocal grid points are searched from uniform
mesh grid points specified by `mesh` and `is_shift`.
`mesh` stores three integers. Reciprocal primitive vectors
are divided by the number stored in `mesh` with (0,0,0) point
centering. The centering can be shifted only half of one mesh
by setting `1` or `true` for each `is_shift` element. If `0` or `false` is set for
`is_shift`, it means there is no shift. This limitation of
shifting enables the irreducible k-point search significantly
faster when the mesh is very dense.

The reducible uniform grid points are returned in reduced
coordinates as `grid_address`. A map between reducible and
irreducible points are returned as `grid_mapping_table` as in the indices of
`grid_address`. The number of the irreducible k-points are
returned as the return value. The time reversal symmetry is
imposed by setting `is_time_reversal`.

!!! compat "Version 0.2"
    The returned mapping table is indexed starting at `1`, not `0` as in Python or C.
"""
function get_ir_reciprocal_mesh(
    cell::Cell,
    mesh,
    is_shift = falses(3);
    is_time_reversal = true,
    symprec = 1e-5,
)
    # Reference: https://github.com/unkcpz/LibSymspg.jl/blob/e912dd3/src/ir-mesh-api.jl#L1-L32
    if !(length(mesh) == length(is_shift) == 3)
        throw(DimensionMismatch("`grid` & `is_shift` must be both length-three vectors!"))
    end
    @assert all(isone(x) || iszero(x) for x in is_shift)
    # Prepare for input
    @unpack lattice, positions, types = _expand_cell(cell)
    mesh = Base.cconvert(Vector{Cint}, mesh)
    is_shift = Base.cconvert(Vector{Cint}, is_shift)
    is_time_reversal = Base.cconvert(Cint, is_time_reversal)
    num_atom = Base.cconvert(Cint, length(types))
    # Prepare for output
    npoints = prod(mesh)
    grid_address = Matrix{Cint}(undef, 3, npoints)  # Julia stores multi-dimensional data in column-major, not row-major (C-style) in memory.
    grid_mapping_table = Vector{Cint}(undef, npoints)
    nir = ccall(
        (:spg_get_ir_reciprocal_mesh, libsymspg),
        Cint,
        (
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Cint},
            Cint,
            Ptr{Cdouble},
            Ptr{Cdouble},
            Ptr{Cint},
            Cint,
            Cdouble,
        ),
        grid_address,
        grid_mapping_table,
        mesh,
        is_shift,
        is_time_reversal,
        lattice,
        positions,
        types,
        num_atom,
        symprec,
    )
    if nir <= 0
        throw(SpglibError("Something wrong happens when finding mesh!"))
    end
    grid_mapping_table .+= 1  # See https://github.com/singularitti/Spglib.jl/issues/56
    return nir, grid_mapping_table, grid_address
end

function get_stabilized_reciprocal_mesh(
    rotations,
    mesh,
    is_shift = falses(3);
    qpoints = [[0, 0, 0]],
    is_time_reversal = true,
)
    if !(length(mesh) == length(is_shift) == 3)
        throw(DimensionMismatch("`grid` & `is_shift` must be both length-three vectors!"))
    end
    @assert all(isone(x) || iszero(x) for x in is_shift)
    @assert size(qpoints, 2) == 3
    npoints = prod(mesh)
    grid_address = Matrix{Cint}(undef, npoints, 3)
    mapping = Vector{Cint}(undef, npoints)
    nir = ccall(
        (:spg_get_stabilized_reciprocal_mesh, libsymspg),
        Cint,
        (
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Cint},
            Cint,
            Cint,
            Ptr{Cint},
            Cint,
            Ptr{Cint},
        ),
        grid_address,
        mapping,
        mesh,
        is_shift,
        is_time_reversal,
        length(rotations),
        rotations,
        length(qpoints),
        qpoints,
    )
    if nir <= 0
        throw(SpglibError("Something wrong happens when finding mesh!"))
    end
    mapping .+= 1  # See https://github.com/singularitti/Spglib.jl/issues/56
    return nir, mapping, grid_address
end


text:
@testset "Join" begin

#@test TypedTables.testjoin(1, 1) == true
#@test TypedTables.testjoin(1, 2) == false
#@test TypedTables.dojoin(1, 1) == 1

#@test join(@row(A::Int64=1, B::Float64=2.0), @row(A::Int64=1, B::Float64=2.0)) === Nullable{Row{@index(A::Int64, B::Float64), Tuple{Int64, Float64}}}(@row(A::Int64=1, B::Float64=2.0))
#@test isnull(join(@row(A::Int64=1, B::Float64=2.0), @row(A::Int64=2, B::Float64=4.0)))
#@test join(@row(A::Int64=1, B::Float64=2.0), @row(B::Float64=2.0, C::Bool=true)) === Nullable{Row{@index(B::Float64, A::Int64, C::Bool), Tuple{Float64, Int64, Bool}}}(@row(B::Float64=2.0, A::Int64=1, C::Bool=true))

@test join(@Table(A=[1,2,3], B=[2.0,4.0,6.0]), @Table(B=[2.0,4.0,5.0,7.0], C=[true,false,true,false])) == @Table(B=[2.0,4.0], A=[1,2], C=[true,false])

end


text:
include("../utils.jl")

@testset "General group tests" begin
    @testset "Not implemented operation" begin
        G = GroupManifold(NotImplementedManifold(), NotImplementedOperation())
        @test repr(G) == "GroupManifold(NotImplementedManifold(), NotImplementedOperation())"
        x = [1.0, 2.0]
        v = [2.0, 3.0]
        eg = Identity(G)
        @test repr(eg) === "Identity($(G))"

        @test_throws Exception Identity(G, Val(true))
        @test_throws ErrorException Identity(G, Val(false))

        @test is_decorator_manifold(G) === Val(true)

        @test Manifolds.is_decorator_group(G) === Val(true)
        @test Manifolds.is_decorator_group(NotImplementedManifold()) === Val(false)
        @test Manifolds.is_decorator_group(G, Val(true)) === Val(false)
        @test Manifolds.is_decorator_group(G, Val(false)) === Val(false)

        @test base_group(G) === G

        if VERSION ≥ v"1.3"
            @test NotImplementedOperation(NotImplementedManifold()) === G
            @test (NotImplementedOperation())(NotImplementedManifold()) === G
        end

        @test_throws ErrorException copyto!(x, eg)

        @test_throws ErrorException inv!(G, x, x)
        @test_throws ErrorException inv!(G, x, eg)
        @test_throws ErrorException inv(G, x)

        @test_throws ErrorException identity!(G, x, x)
        @test_throws ErrorException identity(G, x)

        @test_throws ErrorException compose(G, x, x)
        @test_throws ErrorException compose(G, x, eg)
        @test_throws ErrorException compose!(G, x, eg, x)
        @test_throws ErrorException compose!(G, x, x, eg)
        @test_throws ErrorException compose!(G, x, x, x)
        @test_throws ErrorException compose!(G, x, eg, eg)

        @test_throws ErrorException translate(G, x, x)
        @test_throws ErrorException translate(G, x, x, LeftAction())
        @test_throws ErrorException translate(G, x, x, RightAction())
        @test_throws ErrorException translate!(G, x, x, x)
        @test_throws ErrorException translate!(G, x, x, x, LeftAction())
        @test_throws ErrorException translate!(G, x, x, x, RightAction())

        @test_throws ErrorException inverse_translate(G, x, x)
        @test_throws ErrorException inverse_translate(G, x, x, LeftAction())
        @test_throws ErrorException inverse_translate(G, x, x, RightAction())
        @test_throws ErrorException inverse_translate!(G, x, x, x)
        @test_throws ErrorException inverse_translate!(G, x, x, x, LeftAction())
        @test_throws ErrorException inverse_translate!(G, x, x, x, RightAction())

        @test_throws ErrorException translate_diff(G, x, x, v)
        @test_throws ErrorException translate_diff(G, x, x, v, LeftAction())
        @test_throws ErrorException translate_diff(G, x, x, v, RightAction())
        @test_throws ErrorException translate_diff!(G, v, x, x, v)
        @test_throws ErrorException translate_diff!(G, v, x, x, v, LeftAction())
        @test_throws ErrorException translate_diff!(G, v, x, x, v, RightAction())

        @test_throws ErrorException inverse_translate_diff(G, x, x, v)
        @test_throws ErrorException inverse_translate_diff(G, x, x, v, LeftAction())
        @test_throws ErrorException inverse_translate_diff(G, x, x, v, RightAction())
        @test_throws ErrorException inverse_translate_diff!(G, v, x, x, v)
        @test_throws ErrorException inverse_translate_diff!(G, v, x, x, v, LeftAction())
        @test_throws ErrorException inverse_translate_diff!(G, v, x, x, v, RightAction())
    end

    @testset "Action direction" begin
        @test switch_direction(LeftAction()) == RightAction()
        @test switch_direction(RightAction()) == LeftAction()
    end

    @testset "Addition operation" begin
        G = GroupManifold(NotImplementedManifold(), Manifolds.AdditionOperation())
        test_group(G, [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])

        x = [1.0, 2.0]
        ge = Identity(G)
        @test zero(ge) === ge
        @test number_eltype(ge) == Bool
        @test copyto!(ge, ge) === ge
        y = similar(x)
        copyto!(y, ge)
        @test y ≈ zero(x)
        @test ge - x == -x
        @test x - ge === x
        @test ge - ge === ge
        @test ge + x ≈ x
        @test x + ge ≈ x
        @test ge + ge === ge
        @test -ge === ge
        @test +ge === ge
        @test ge * 1 === ge
        @test 1 * ge === ge
        @test ge * ge === ge
        @test ge(x) ≈ zero(x)
        @test inv(G, x) ≈ -x
        @test inv(G, ge) === ge
        @test identity(G, x) ≈ zero(x)
        @test identity(G, ge) === ge
        y = similar(x)
        identity!(G, y, x)
        @test y ≈ zero(x)
        @test compose(G, x, x) ≈ x + x
        @test compose(G, x, ge) ≈ x
        @test compose(G, ge, x) ≈ x
        @test compose(G, ge, ge) == ge
        compose!(G, y, x, x)
        @test y ≈ x + x
        compose!(G, y, x, ge)
        @test y ≈ x
        compose!(G, y, ge, x)
        @test y ≈ x

        y = identity(G, x)
        @test isapprox(y, ge; atol=1e-10)
        @test isapprox(ge, y; atol=1e-10)
        @test isapprox(ge, ge)
    end

    @testset "Multiplication operation" begin
        G = GroupManifold(NotImplementedManifold(), Manifolds.MultiplicationOperation())
        test_group(G, [[1.0 2.0; 3.0 4.0], [2.0 3.0; 4.0 5.0], [3.0 4.0; 5.0 6.0]])

        x = [1.0 2.0; 2.0 3.0]
        ge = Identity(G)
        @test number_eltype(ge) == Bool
        @test copyto!(ge, ge) === ge
        y = similar(x)
        copyto!(y, ge)
        @test y ≈ one(x)
        @test one(ge) === ge
        @test transpose(ge) === ge
        @test det(ge) == 1
        @test ge * x ≈ x
        @test x * ge ≈ x
        @test ge * ge === ge
        @test inv(ge) === ge
        @test *(ge) === ge

        @test x / ge ≈ x
        @test ge \ x ≈ x
        @test ge / ge === ge
        @test ge \ ge === ge
        @test ge / x ≈ inv(x)
        @test x \ ge ≈ inv(x)
        y = similar(x)
        @test LinearAlgebra.mul!(y, x, ge) === y
        @test y ≈ x
        y = similar(x)
        @test LinearAlgebra.mul!(y, ge, x) === y
        @test y ≈ x
        y = similar(x)
        @test LinearAlgebra.mul!(y, ge, ge) === y
        @test y ≈ one(y)

        @test ge(x) ≈ one(x)
        @test inv(G, x) ≈ inv(x)
        @test inv(G, ge) === ge
        @test identity(G, x) ≈ one(x)
        @test identity(G, ge) === ge
        y = similar(x)
        identity!(G, y, x)
        @test y ≈ one(x)
        @test_throws ErrorException identity!(G, [0.0], ge)
        @test compose(G, x, x) ≈ x * x
        @test compose(G, x, ge) ≈ x
        @test compose(G, ge, x) ≈ x
        @test compose(G, ge, ge) == ge
        compose!(G, y, x, x)
        @test y ≈ x * x
        compose!(G, y, x, ge)
        @test y ≈ x
        compose!(G, y, ge, x)
        @test y ≈ x

        @testset "identity optimization" begin
            x2 = copy(x)
            identity!(G, x2, x)
            x3 = copy(x)
            invoke(identity!, Tuple{AbstractGroupManifold{Manifolds.MultiplicationOperation}, Any, AbstractMatrix}, G, x3, x)
            @test isapprox(G, x2, x3)
        end
    end
end

struct NotImplementedAction <: AbstractGroupAction{LeftAction}
end

@testset "General group action tests" begin
    @testset "Not implemented operations" begin
        A = NotImplementedAction()
        x = [1.0, 2.0]
        a = [1.0, 2.0]
        v = [1.0, 2.0]

        @test_throws ErrorException base_group(A)
        @test_throws ErrorException g_manifold(A)
        @test_throws ErrorException apply(A, a, x)
        @test_throws ErrorException apply!(A, x, a, x)
        @test_throws ErrorException inverse_apply(A, a, x)
        @test_throws ErrorException inverse_apply!(A, x, a, x)
        @test_throws ErrorException apply_diff(A, a, x, v)
        @test_throws ErrorException apply_diff!(A, v, x, a, v)
        @test_throws ErrorException inverse_apply_diff(A, a, x, v)
        @test_throws ErrorException inverse_apply_diff!(A, v, x, a, v)
        @test_throws ErrorException compose(A, a, a)
        @test_throws ErrorException compose!(A, a, a, a)
        @test_throws ErrorException optimal_alignment(A, x, x)
        @test_throws ErrorException optimal_alignment!(A, a, x, x)
    end
end


text:
struct Bus
    id::String
    generators::Vector{Generator}

    # Constructor
    function Bus(id::String)
        return new(id, Vector{Generator}())
    end
end

# Never use it directly, only use add function in Network!!!
function add_new_generator!(bus::Bus, generator_id::String,
                            type::GeneratorType, p_min::Float64, p_max::Float64,
                            start_cost::Float64, prop_cost::Float64, dmo::Dates.Second, dp::Dates.Second
                            )::Generator
    generator::Generator = Generator(generator_id, bus.id,
                                    type, p_min, p_max, start_cost, prop_cost, dmo, dp)
    push!(bus.generators, generator)
    return generator
end

function get_generators(bus::Bus)::Vector{Generator}
    return bus.generators
end

################
## INFO / LOG ##
################

function get_info(bus::Bus)::String
    info::String =
        string(bus.id)
    return info
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using Test
import Base: return_types

using Compat

import Distributed: RemoteChannel

using Neovim
import Neovim: get_buffers, set_line, get_line
import Neovim: vim_eval, command, get_var, set_var
import Neovim: on_notify, on_request
nvim, proc = Neovim.nvim_spawn()

# Test buffer
buf = get_buffers(nvim)[1]
@assert isa(buf, Buffer)
set_line(buf, 0, "some text")
text = get_line(buf, 0)
@assert text == "some text"

# test high-level buffer interface
buf = current_buffer(nvim)
@assert isa(buf, Buffer)
@assert buf[:] == ["some text"]
@assert buf[1] == "some text"

buf[:] = ["alpha", "beta", "gamma"]
@assert buf[:] == ["alpha", "beta", "gamma"]
@assert buf[1] == "alpha"
@assert buf[1:3] == ["alpha", "beta", "gamma"]
@assert buf[1:1] == ["alpha"]
@assert buf[end - 1:end] == ["beta", "gamma"]
@assert buf[1 + end - 1:end] == ["gamma"]
@assert buf[end - 1 + 1:1 + end - 1] == ["gamma"]
@assert buf[end - 1:2] == ["beta"]
@assert buf[end + 1 - 1] == "gamma"
@assert buf[1 + end - 1] == "gamma"
@assert buf[2:1] == []
@assert buf[1:0] == []

@test_throws BoundsError buf[0]
@test_throws BoundsError buf[-1]
@test_throws BoundsError buf[end + 1]
@test_throws BoundsError buf[-1:1]
@test_throws BoundsError buf[0:1]
@test_throws BoundsError buf[-1:0]
@test_throws BoundsError buf[end + 1:1]
@test_throws BoundsError buf[1:end + 1]
@test_throws BoundsError buf[1:1 + end]
@test_throws BoundsError buf[1 + end:1]
@test_throws BoundsError buf[end - 1:end + 1]
@test_throws BoundsError buf[end:end + 1]
@test_throws BoundsError buf[end + 1:end + 1]

@assert (buf[2] = "beta-ish") == "beta-ish"
@assert buf[:] == ["alpha", "beta-ish", "gamma"]
@assert (buf[end] = "gamma-ish") == "gamma-ish"
@assert buf[:] == ["alpha", "beta-ish", "gamma-ish"]
@assert (buf[2:3] = ["b", "c", "d"]) == ["b", "c", "d"]
@assert buf[:] == ["alpha", "b", "c", "d"]
@assert (buf[end - 1:end - 2] = ["boom"]) == ["boom"]
@assert buf[:] == ["alpha", "b", "boom", "c", "d"]
@assert (buf[end - 1:1 + end - 2] = ".") == "."
@assert buf[:] == ["alpha", "b", "boom", ".", "d"]

@test_throws BoundsError buf[0] = "twas"
@test_throws BoundsError buf[-1] = "brillyg"
@test_throws BoundsError buf[end + 1] = "and"
@test_throws BoundsError buf[-1:1] = "the"
@test_throws BoundsError buf[0:1] = "slythy"
@test_throws BoundsError buf[-1:0] = "toves"
@test_throws BoundsError buf[end + 1:1] = "did"
@test_throws BoundsError buf[1:end + 1] = "gyre"
@test_throws BoundsError buf[1:1 + end] = "and"
@test_throws BoundsError buf[1 + end:1] = "gymble"
@test_throws BoundsError buf[end - 1:end + 1] = "in"
@test_throws BoundsError buf[end:end + 1] = "the"
@test_throws BoundsError buf[end + 1:end + 1] = "wabe"

deleteat!(buf, 2)
@assert buf[:] == ["alpha", "boom", ".", "d"]
deleteat!(buf, 3:4)
@assert buf[:] == ["alpha", "boom"]

push!(buf, "the end")
@assert buf[:] == ["alpha", "boom", "the end"]
pushfirst!(buf, "stuff")
@assert buf[:] == ["stuff", "alpha", "boom", "the end"]
@assert length(buf) == 4

# test high-level cursor api
win = current_window(nvim)
cursor!(win, 1, 1)
@assert cursor(win) == (1, 1)
cursor!(win, 4, 2)
@assert cursor(win) == (4, 2)

# test eval
@assert vim_eval(nvim, "2+2") == 4

# test command and async behavior
command(nvim, "let g:test = []")
@sync for i in 1:5
    @async for j in 1:5
        command(nvim, "call add(g:test, [$i,$j])")
        if rand() > 0.5; sleep(0.001) end
    end
end

res = get_var(nvim, "test")
@assert length(res) == 25
for i in 1:5
    @assert all(indexin([[i,j] for j in 1:5], res) .> 0)
end

# test events
struct TestHandler
    r::RemoteChannel
end
on_notify(h::TestHandler, c, name, args) = put!(h.r, (name, args))
function on_request(h::TestHandler, c, serial, name, args)
    if name == "do_stuff"
        reply_result(c, serial, 10 * args[1] + args[2])
    end
end

# notification
ref = RemoteChannel()
nvim, proc = nvim_spawn(TestHandler(ref))
command(nvim, "call rpcnotify($(nvim.channel_id), 'mymethod', 10, 20)")
@assert take!(ref) == ("mymethod", Any[10, 20])

# request
@assert vim_eval(nvim, "100+rpcrequest($(nvim.channel_id), 'do_stuff', 2, 3)") == 123

# type stability of generated functions
# TODO(smolck): This doesn't pass now
# @assert return_types(Neovim.get_buffers, (NvimClient,)) == [Vector{Buffer}]

@assert return_types(Neovim.command, (NvimClient, String)) == [Nothing]
@assert return_types(Neovim.get_current_line, (NvimClient,)) == [String]
@assert return_types(Neovim.is_valid, (Tabpage,)) == [Bool]

# TODO(smolck): These don't pass now either
# @assert return_types(Neovim.get_height, (Window,)) == [Int]
# @assert return_types(Neovim.get_mark, (Buffer, String)) == [Tuple{Int,Int}]

# as ByteString isn't concrete anyway, this doesn't give that much really
# @assert return_types(Neovim.get_line_slice, (Buffer,Int,Int,Bool,Bool)) == [Vector{TypeVar(:_,None,ByteString)}]

# test host
hostdir = dirname(dirname(@__FILE__))
plugdir = joinpath(dirname(@__FILE__), "hosttest")

# fake initialization for :UpdateRemotePlugins
vimdir = mktempdir(cleanup=false)
nvimrc = joinpath(vimdir, "nvimrc")
open(f -> nothing, nvimrc,"w")
ENV["MYVIMRC"] = nvimrc
ENV["NEOVIM_JL_DEBUG"] = "templog"
println(vimdir)
ENV["NVIM_RPLUGIN_MANIFEST"] = vimdir * "/rplugin.vim"
ENV["NVIM_LOG_FILE"] = vimdir * "/.nvimlog"
rtp = "set rtp+=$hostdir,$plugdir"
juliap = "let g:julia_host_prog = '$(joinpath(Sys.BINDIR, "julia"))'"
cmd = `nvim -u $nvimrc -i NONE --cmd $rtp --cmd $juliap -c UpdateRemotePlugins -c q`
run(cmd)
run(`cat templog`)

try
    local ref = RemoteChannel()
    n, p = nvim_spawn(TestHandler(ref), cmd=`nvim -u $nvimrc -i NONE --cmd $rtp --cmd $juliap --embed --headless`)

    @assert vim_eval(n, "TestFun('a',3)") == "TestFun got a, 3"

    command(n, "call AsyncFun($(n.channel_id), {'alfa':1, 'omega':'theend'})")
    @assert take!(ref) == ("AsyncReply", Any[Dict("alfa" => 1, "omega" => "theend")])

    command(n, "new")
    b = current_buffer(n)
    b[1:5] = "line"
    command(n, "2,3JLCommand text")
    @assert b[:] == ["line","line","line","line","line","text"]
    @assert get_var(n, "therange") == [2,3]

# TODO: more specific test? (somewhere in there should be "KABOOM!")
    @test_throws ErrorException command(n, "Explode")

    set_var(n, "zinged", 5)
    command(n, "do User zing")
    @assert get_var(n, "zinged") == 6
    command(n, "do User yoink")
# shouldn't change
    @assert get_var(n, "zinged") == 6

    @assert vim_eval(n, "Global()") == 3
    @assert vim_eval(n, "Local()") == 37

finally
# for debugging tests:
    run(`cat templog`)
end


text:
function MIMMPSolver(A,b,C,d,VariableDef;Solver = "CB", XLB=[],XUB=[], rEPS = 1e-6,aEPS = 1e-6,TimeLim = 2e+10)
	VariableDef = Int.(VariableDef)
	if Solver == "CB"
		GUB, OptX = CBMIMMPSolver(A,b,C,d,VariableDef,XLB=XLB,XUB=XUB, rEPS = rEPS,aEPS = aEPS,TimeLim = TimeLim)
	else
		NODE,VARIABLE,ENHANCED_UB,MU,ENHANCED_CUT = 4,	3,	1,	1,1;
		GUB, OptX = DBMIMMPSolver(A,b,C,d,VariableDef,XLB=XLB,XUB=XUB,MU = MU,EnhancedUB = ENHANCED_UB,EnhancedCut = ENHANCED_CUT ,SelectVariableMode = VARIABLE,SelectNodeMode = NODE,TimeLim=TimeLim,rEPS = rEPS, aEPS = aEPS)
	end
	return GUB, OptX;
end