text:
"""
    AbstractTreeNode

Abstract type of tree nodes
"""

abstract type AbstractTreeNode end

"""
    TreeNode

Tree node stores information of stage problem

    - `id`: index of node
    - `stage_builder`: adds variables and constraints to model with input (tree::Tree, subtree::SubTree, node::SubTreeNode)
    - `parent`: index of parent node
    - `children`: indices of child nodes
    - `stage`: current stage
    - `ξ`: current scenario
    - `p`: probability of node
"""

mutable struct TreeNode <: AbstractTreeNode
    id::Int
    stage_builder::Union{Nothing,Function}
    parent::Int
    children::Vector{Tuple{Int, Float64}}
    stage::Int
    ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}
    p::Float64

    function TreeNode(id::Int, parent::Int, stage::Int, ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}, p::Float64)
        tn = new()
        tn.id = id
        tn.stage_builder = nothing
        tn.parent = parent
        tn.children = Vector{Tuple{Int, Float64}}()
        tn.stage = stage
        tn.ξ = ξ
        tn.p = p
        return tn
    end
end

function TreeNode(ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}})
    return TreeNode(1, 0, 1, ξ, 1.0)
end

get_id(node::TreeNode) = node.id
get_parent(node::TreeNode) = node.parent
get_children(node::TreeNode) = node.children
get_stage(node::TreeNode) = node.stage
get_scenario(node::TreeNode) = node.ξ
get_probability(node::TreeNode) = node.p
function set_stage_builder!(node::TreeNode, func::Function)
    node.stage_builder = func
end

"""
    AbstractTree

Abstract type of tree
"""

abstract type AbstractTree end

"""
    Tree

Tree keeps information of tree nodes.

    - `nodes`: list of nodes
"""

mutable struct Tree <: AbstractTree
    nodes::Dict{Int,TreeNode}
end

Tree() = Tree(Dict{Int,TreeNode}())

Tree(ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}) = Tree(Dict{Int,TreeNode}(1 => TreeNode(ξ)))

"""
    add_node!

adds abstract tree node to abstract tree

# Arguments
    - `tree`: Abstract tree
    - `node`: Abstract tree node
"""

function add_node!(tree::AbstractTree, node::AbstractTreeNode)
    @assert !haskey(tree.nodes, get_id(node))
    id = get_id(node)
    tree.nodes[id] = node
end

"""
    add_child!

creates child nodee and adds to tree

# Arguments
    - `tree`: Tree
    - `pt`: Parent ID
    - 'ξ': Dictionary of random variable name as symbols and scenarios
    - `prob`: probability of transitioning from the parent to this node
"""

function add_child!(tree::Tree, pt::Int, ξ::Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}, prob::Float64)::Int
    #   adds child node to tree.nodes[pt]
    @assert haskey(tree.nodes, pt)                                  # check if pt is valid
    stage = get_stage(tree, pt) + 1                                 # get new stage value
    p = get_probability(tree, pt) * prob                            # get new node probability
    id = length(tree.nodes) + 1                                     # get node id
    add_node!(tree, TreeNode(id, pt, stage, ξ, p ))                 # create node and add to tree
    push!(get_children(tree, pt), (id, prob))                      # push id to parent node children
    return id
end

get_children(tree, id) = get_children(tree.nodes[id])
get_parent(tree, id) = get_parent(tree.nodes[id])
get_stage(tree, id) = get_stage(tree.nodes[id])
get_scenario(tree, id) = get_scenario(tree.nodes[id])
get_probability(tree, id) = get_probability(tree.nodes[id])
function set_stage_builder!(tree, id, func::Function)
    set_stage_builder!(tree.nodes[id], func)
end

"""
    SubTreeNode

Contains TreeNode and other information that are used for dual decomposition.
    - `treenode`: an instance of TreeNode
    - `weight`: multiplied to the objective
    - `in`: dictionary of incoming variables from the previous stage
    - `out`: dictionary of outgoing variablees to the subsequent stage
    - `obj`: objective expression
"""

mutable struct SubTreeNode <: AbstractTreeNode
    treenode::TreeNode
    weight::Float64
    in::Dict{Symbol, Union{JuMP.VariableRef, <:AbstractArray{JuMP.VariableRef}}}    # incoming variables
    out::Dict{Symbol, Union{JuMP.VariableRef, <:AbstractArray{JuMP.VariableRef}}}   # outgoing variables
    obj::Union{Float64, JuMP.AbstractJuMPScalar}
    function SubTreeNode(treenode::TreeNode, weight::Float64)
        stn = new()
        stn.treenode = treenode
        stn.weight = weight
        stn.in = Dict{Symbol, Union{JuMP.VariableRef, <:AbstractArray{JuMP.VariableRef}}}()
        stn.out = Dict{Symbol, Union{JuMP.VariableRef, <:AbstractArray{JuMP.VariableRef}}}()
        stn.obj = 0.0
        return stn
    end
end

get_id(node::SubTreeNode) = get_id(node.treenode)
get_parent(node::SubTreeNode) = get_parent(node.treenode)
get_children(node::SubTreeNode) = get_children(node.treenode)
get_stage(node::SubTreeNode) = get_stage(node.treenode)
get_scenario(node::SubTreeNode) = get_scenario(node.treenode)
get_probability(node::SubTreeNode) = get_probability(node.treenode)

function set_input_variable!(nd::SubTreeNode, symb::Symbol, var::Union{JuMP.VariableRef, <:AbstractArray{JuMP.VariableRef}})
    nd.in[symb] = var
end

function set_output_variable!(nd::SubTreeNode, symb::Symbol, var::Union{JuMP.VariableRef, <:AbstractArray{JuMP.VariableRef}})
    nd.out[symb] = var
end

function set_stage_objective(nd::SubTreeNode, obj::Union{Float64, JuMP.AbstractJuMPScalar})
    nd.obj = obj
end

"""
    SubTree

Used for creeating and keeeping subproblems

    - `block_id`: ID of the subtrees
    - `nodes`: dictionary of SubTreeNodes
    - `model`: JuMP model
    - `parent`: ID of parent if exists (default is 0)
    - `root`: ID of the root node of the subtree (default is 1)
"""

mutable struct SubTree <: AbstractTree
    block_id::Int
    nodes::Dict{Int,SubTreeNode}     # list of nodes
    model::JuMP.AbstractModel
    parent::Int
    root::Int
    function SubTree(block_id::Int)
        return new(block_id, Dict{Int,SubTreeNode}(), JuMP.Model(), 0, 1)
    end
end

"""
    create_subtree!

creates subtree from vector of nodes and adds coupling variables

# Arguments
    - `tree`: Tree
    - `block_id`: ID of subtree
    - 'coupling_variables': list of coupling variables to be modified
    - `nodes`: vector of nodes
"""

function create_subtree!(tree::Tree, block_id::Int, coupling_variables::Vector{CouplingVariableRef}, nodes::Vector{Tuple{TreeNode,Float64}})::SubTree
    subtree = SubTree(block_id)
    # add nodes to subtree
    for (node, weight) in nodes
        subnode = SubTreeNode(node, weight)
        add_node!(subtree, subnode)
    end

    obj = 0
    for (id, subnode) in subtree.nodes
        subnode.treenode.stage_builder(tree, subtree, subnode) # make macro for changing variable names and constraint names to include node id
        obj += subnode.weight * subnode.obj
    end
    set_objective!(subtree, MOI.MIN_SENSE, obj)

    #
    for (id, subnode) in subtree.nodes
        couple_common_variables!(coupling_variables, block_id, subnode)
        parent = get_parent(subnode)
        if parent!=0 && haskey(subtree.nodes, parent)
            add_links!(subtree, id, parent)
        elseif parent!=0 # assuming 1st stage node is 1
            subtree.parent = parent
            subtree.root = id
            couple_incoming_variables!(coupling_variables, block_id, subnode)
        end
    end
    return subtree
end

function set_objective!(tree::SubTree, sense::MOI.OptimizationSense, obj::Union{Float64, JuMP.AbstractJuMPScalar})
    JuMP.set_objective(tree.model, sense, obj)
end

"""
    add_links!

creates linking constraints for nodes within the subtree by connecting the incoming and outgoing variables

# Arguments
    - `tree`: SubTree
    - `id`: ID of node
    - 'pt': ID of parent node
"""

function add_links!(tree::SubTree, id::Int, pt::Int)
    node = tree.nodes[id]
    parent = tree.nodes[pt]
    for (symb, var1) in node.in
        var2 = parent.out[symb]
        @constraint(tree.model, var1 .== var2)
    end
end

"""
    couple_common_variables!

couple outgoing variables

# Arguments
    - `coupling_variables`: vector of coupling variables
    - `block_id`: ID of block
    - `node`: SubTreeNode
"""

function couple_common_variables!(coupling_variables::Vector{CouplingVariableRef}, block_id::Int, node::SubTreeNode)
    label = get_id(node)
    for (symb, var) in node.out
        couple_variables!(coupling_variables, block_id, label, symb, var)
    end
end

"""
    couple_incoming_variables!

couple incoming variables with the root node of the subtree

# Arguments
    - `coupling_variables`: vector of coupling variables
    - `block_id`: ID of block
    - `child`: SubTreeNode
"""

function couple_incoming_variables!(coupling_variables::Vector{CouplingVariableRef}, block_id::Int, child::SubTreeNode)
    label = get_parent(child)
    for (symb, var) in child.in
        couple_variables!(coupling_variables, block_id, label, symb, var)
    end
end

"""
    couple_variables!

adds variables to coupling_variables

# Arguments
    - `coupling_variables`: vector of coupling variables
    - `block_id`: ID of block
    - `label`: ID of node
    - `symb`: symbol for variable
    - `var`: variable or vector of variables
"""

function couple_variables!(coupling_variables::Vector{CouplingVariableRef}, block_id::Int, label::Int, symb::Symbol,
        var::JuMP.VariableRef)
    push!(coupling_variables, CouplingVariableRef(block_id, [label, symb], var))
end

function couple_variables!(coupling_variables::Vector{CouplingVariableRef}, block_id::Int, label::Int, symb::Symbol,
        var::Array{JuMP.VariableRef})
    for (index, value) in pairs(var)
        push!(coupling_variables, CouplingVariableRef(block_id, [label, symb, Tuple(index)], value))
    end
end

function couple_variables!(coupling_variables::Vector{CouplingVariableRef}, block_id::Int, label::Int, symb::Symbol,
        var::JuMP.Containers.DenseAxisArray{JuMP.VariableRef})
    for (index, value) in pairs(var.data)
        push!(coupling_variables, CouplingVariableRef(block_id, [label, symb, keys(var)[index].I], value))
    end
end

function couple_variables!(coupling_variables::Vector{CouplingVariableRef}, block_id::Int, label::Int, symb::Symbol,
        var::JuMP.Containers.SparseAxisArray{JuMP.VariableRef})
    for (index, value) in var.data
        push!(coupling_variables, CouplingVariableRef(block_id, [label, symb, index], value))
    end
end

"""
    decomposition_not

outputs the entire tree

# Arguments
    - `tree`: Tree
"""

function decomposition_not(tree::Tree)::Vector{Vector{Tuple{TreeNode,Float64}}}
    nodes = Vector{Tuple{TreeNode,Float64}}()
    for (id, node) in tree.nodes
        push!(nodes,(node, get_probability(node)))
    end
    return [nodes]
end

"""
    decomposition_scenario

outputs the scenario decomposition at each leaf nodes

# Arguments
    - `tree`: Tree
"""

function decomposition_scenario(tree::Tree)::Vector{Vector{Tuple{TreeNode,Float64}}}
    node_cluster = Vector{Vector{Tuple{TreeNode,Float64}}}()
    for (id, node) in tree.nodes
        if check_leaf(node)
            nodes = Vector{Tuple{TreeNode,Float64}}()
            prob = get_probability(node)
            current = node
            while true
                push!(nodes,(current, prob))
                pt = get_parent(current)
                current = tree.nodes[pt]
                if check_root(current)
                    push!(nodes,(current, prob))
                    break
                end
            end
            push!(node_cluster, nodes)
        end
    end
    return node_cluster
end

"""
    decomposition_temporal

outputs the temporal decomposition at each nodes

# Arguments
    - `tree`: Tree
"""
#need fix

function decomposition_temporal(tree::Tree)::Vector{Vector{Tuple{TreeNode,Float64}}}
    node_cluster = Vector{Vector{Tuple{TreeNode,Float64}}}()
    for (id, node) in tree.nodes
        push!(node_cluster,[(node, get_probability(node))])
    end
    return node_cluster
end

function check_leaf(node::TreeNode)::Bool
    if length(get_children(node)) == 0
        return true
    else
        return false
    end
end

function check_root(node::TreeNode)::Bool
    if get_parent(node) == 0
        return true
    else
        return false
    end
end


text:
visitors
rooms
hoursNumber

type Individual
	meetings
	indstrialistSchedule
	researcherSchedule
	roomSchedule
	f
	g
end

#PRE Constraints respected for the appointment
#PRE: the meeting was not appointed
function appointment!(individual::Individual, industrialist::Int, researcher::Int, room::Int, timeSlot::Int)
	individual.meetings[industrialist][researcher] = [room,timeSlot]
	individual.indstrialistSchedule[industrialist][timeSlot] = [researcher,room]
	individual.researcherSchedule[researcher][timeSlot] = [industrialist,room]
	individual.roomSchedule[room][timeSlot] = [industrialist,researcher]
	individual.f = individual.f + visitors[industrialist][researcher]
	if timeSlot > 1
		industrialistRoom = individual.indstrialistSchedule[industrialist][timeSlot-1][1]
		researcherRoom = researcher.ResearcherSchedule[researcher][timeSlot-1][1]
		if industrialistRoom > 0
			individual.g = individual.g + rooms[industrialistRoom][room]
		end
		if researcherRoom > 0
			individual.g = individual.g + rooms[researcherRoom][room]
		end
	end
	if timeSlot < hoursNumber
		industrialistRoom = individual.indstrialistSchedule[industrialist][timeSlot+1][1]
		researcherRoom = researcher.ResearcherSchedule[researcher][timeSlot+1][1]
		if industrialistRoom > 0
			individual.g = individual.g + rooms[industrialistRoom][room]
		end
		if researcherRoom > 0
			individual.g = individual.g + rooms[researcherRoom][room]
		end
	end
end

function remove(individual::Individual, industrialist::Int, researcher::Int, room::Int, timeSlot::Int)
	individual.f = individual.f - visitors[industrialist][researcher]
	if timeSlot > 1
		industrialistRoom = individual.indstrialistSchedule[industrialist][timeSlot-1][1]
		researcherRoom = researcher.ResearcherSchedule[researcher][timeSlot-1][1]
		if industrialistRoom > 0
			individual.g = individual.g - rooms[industrialistRoom][room]
		end
		if researcherRoom > 0
			individual.g = individual.g - rooms[researcherRoom][room]
		end
	end
	if timeSlot < hoursNumber
		industrialistRoom = individual.indstrialistSchedule[industrialist][timeSlot+1][1]
		researcherRoom = researcher.ResearcherSchedule[researcher][timeSlot+1][1]
		if industrialistRoom > 0
			individual.g = individual.g - rooms[industrialistRoom][room]
		end
		if researcherRoom > 0
			individual.g = individual.g - rooms[researcherRoom][room]
		end
	end
	individual.meetings[industrialist][researcher] = [0,0]
	individual.indstrialistSchedule[industrialist][timeSlot] = [0,0]
	individual.researcherSchedule[researcher][timeSlot] =[0,0]
	individual.roomSchedule[room][timeSlot] = [0,0]
end

function crossoverMeetingScheduleIntersection(father::Individual, mother::Individual)
	choice=ones(::Int, hoursNumber)
	for h in 1:hoursNumber
		choice[h]=rand(1:2)
	end
	child = Individual(father)
	for i in 1:visitors[2]
		for r in 1:visitors[3]
			hour = father.meetings[i][r][2]
			if hour == mother.meetings[i][r][2]
				if choice[hour] == 1
					appointment!(child, i, r, father.meetings[i][r][1], hour)
				else
					appointment!(child, i, r, mother.meetings[i][r][1], hour)
				end
			end
		end
	end
	grasp!(child)
	return child
end

function crossoverindustrialistScheduleIntersection(father::Individual, mother::Individual)
	child = Individual(father)
	choice=ones(::Int, hoursNumber)
	for h in 1:hoursNumber
		choice[h]=rand(1:2)
	end
	for i in 1:visitors[2]
		for h in 1:hoursNumber
			room = father.indstrialistSchedule[i][h][2]
			if room == mother.indstrialistSchedule[i][h][2]
				if choice[h] == 1
					r = father.indstrialistSchedule[i][h][1]
					if child.meetings[i][r][1]==0
						appointment!(child, i, r, room, h)
					end
				else
					r = mother.indstrialistSchedule[i][h][1]
					if child.meetings[i][r][1]==0
						appointment!(child, i, r, room, h)
					end
				end
			else

			end
		end
	end

	graps!(child)
	return child
end

function crossoverResearcherScheduleIntersection(father::Individual, mother::Individual)
	child = Individual(father)
	choice=ones(::Int, hoursNumber)
	for h in 1:hoursNumber
		choice[h]=rand(1:2)
	end
	for r in 1:visitors[3]
		for h in 1:hoursNumber
			room = father.researcherSchedule[r][h][2]
			if room == mother.researcherSchedule[r][h][2]
				if choice[h] == 1
					i = father.researcherSchedule[r][h][1]
					if child.meetings[i][r][1]==0
						appointment!(child, i, r, room, h)
					end
				else
					i = mother.researcherSchedule[r][h][1]
					if child.meetings[i][r][1]==0
						appointment!(child, i, r, room, h)
					end
				end
			else

			end
		end
	end

	graps!(child)
	return child
end

function crossoverRandomScheduleUnion(father::Individual, mother::Individual)
	child = Individual(father)

	for h in 1:hoursNumber
		for s in 1:rooms[2]
			i1 = father.roomSchedule[s][h][1]
			r1 = father.roomSchedule[s][h][2]
			i2 = mother.roomSchedule[s][h][1]
			r2 = mother.roomSchedule[s][h][2]
			if i1 > 0 && i2 == 0 && child.meetings[i1,r1] == 0
				appointment!(child, i1, r1, s, h)
			elseif i1 == 0 && i2 > 0 && child.meetings[i2,r2] == 0
				appointment!(child, i2, r2, s, h)
			elseif i1 > 0 && i2 > 0
				if child.meetings[i2,r2] == 0
					if child.meetings[i1,r1] == 0
						choice = rand(1:2)
						if choice == 1
							appointment!(child, i1, r1, s, h)
						else
							appointment!(child, i2, r2, s, h)
						end # choice == 1
					else
							appointment!(child, i2, r2, s, h)
					end # child.meetings[i1,father.roomSchedule[s][h][2]] == 0
				elseif child.meetings[i1,r1] == 0
					appointment!(child, i1, r1, s, h)
				end # child.meetings[i2,mother.roomSchedule[s][h][2]] == 0
			end	# i1 > 0 && i2 == 0 && child.meetings[i1,father.roomSchedule[s][h][2]] == 0
		end # s in 1:rooms[2]
	end # h in 1:hoursNumber
end # crossoverRandomScheduleUnion function

	graps!(child)
	return child
end

# #return the distance
# #beforeAfter: int=-1 if we look the change with the precedent hour, +1 if we look the change with the succeed hour
# #PRE: not room1=room2=0
# function exchangeI(beforeAfter, visitors, rooms, hoursNumber, timeSlot, room1, room2, indexi)
# 	if beforeAfter+timeSlot > 0 && beforeAfter+timeSlot < hoursNumber
# 		indexRoom = IndstrialistSchedule[indexi][timeSlot+beforeAfter][1]
# 	end
# 	if 	indexRoom > 0
# 		if room1 == 0
# 			return rooms[1][indexRoom][room2]
# 		elseif room2=0
# 			return -rooms[1][indexRoom][room2]
# 		else
# 			return -rooms[1][indexRoom][room1] + rooms[1][indexRoom][room2]
# 		end
# 	else
# 		return 0
# 	end
# else
# 	return 0
# end
# end
#
# #return the distance
# #beforeAfter: int=-1 if we look the change with the precedent hour, +1 if we look the change with the succeed hour
# function exchangeR(beforeAfter, visitors, rooms, hoursNumber, timeSlot, room1, room2, indexr)
# 	if beforeAfter + timeSlot > 0 && beforeAfter + timeSlot < hoursNumber
# 		indexRoom = ResearcherSchedule[indexr][timeSlot+beforeAfter][1]
# 	end
# 	if 	indexRoom > 0
# 		if room1 == 0
# 			return rooms[1][indexRoom][room2]
# 		elseif room2 = 0
# 			return -rooms[1][indexRoom][room2]
# 		else
# 			return -rooms[1][indexRoom][room1] + rooms[1][indexRoom][room2]
# 		end
# 	else
# 		return 0
# 	end
# else
# 	return 0
# end
# end
#
# function roomExchange!(visitors, rooms, timeSlot, room1, room2)
# 	#pre-requis: les rendez vous fixé aux lieux room1 et room2 eu temps timeSlot existent
# 	RDV = RoomSchedule[room1][timeSlot]
# 	indexi1 = RDV[1]
# 	indexr1 = RDV[2]
# 	RDV = RoomSchedule[room2][timeSlot]
# 	indexi2 = RDV[1]
# 	indexr2 = RDV[2]
#
# 	x[indexi1][indexr1][room1][timeSlot] = 0
# 	x[indexi2][indexr2][room2][timeSlot] = 0
# 	x[indexi1][indexr1][room2][timeSlot]=1
# 	x[indexi2][indexr2][room1][timeSlot]=1
#
# 	secondObjectiv = g
# 	secondObjectiv=secondObjectiv + exchangeI(-1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexi1)
# 	secondObjectiv=secondObjectiv + exchangeI(-1,visitors,rooms,hoursNumber,timeSlot,room2,room1,indexi2)
# 	secondObjectiv=secondObjectiv + exchangeR(-1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexr1)
# 	secondObjectiv=secondObjectiv + exchangeR(-1,visitors,rooms,hoursNumber,timeSlot,room2,room1,indexr2)
# 	secondObjectiv=secondObjectiv + exchangeI(1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexi1)
# 	secondObjectiv=secondObjectiv + exchangeI(1,visitors,rooms,hoursNumber,timeSlot,room2,room1,indexi2)
# 	secondObjectiv=secondObjectiv + exchangeR(1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexr1)
# 	secondObjectiv=secondObjectiv + exchangeR(1,visitors,rooms,hoursNumber,timeSlot,room2,room1,indexr2)
# 	g = secondObjectiv
# end
#
# function swapMeeting!(visitors, rooms, hoursNumber, room1, room2, timeSlot) #Cross two meetings
# 	#pre-requis: les rendez vous fixé aux lieux room1 et room2 eu temps timeSlot existent
# 	RDV=RoomSchedule[room1][timeSlot]
# 	indexi1 = RDV[1]
# 	indexr1 = RDV[2]
# 	RDV=RoomSchedule[room2][timeSlot]
# 	indexi2 = RDV[1]
# 	indexr2 = RDV[2]
# 	secondObjectiv=g
#
# 	if meetings[indexi1][indexr2][1] == 0 && meetings[indexi2][indexr1][1] == 0
# 		f = f - visitors[1][indexi1][indexr1] - visitors[1][indexi2][indexr2] + visitors[1][indexi1][indexr2] + visitors[1][indexi2][indexr1]
# 		g1 = secondObjectiv+exchangeI(-1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexi1)
# 		g1 = g1 + exchangeI(1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexi1)
# 		g1 = secondObjectiv+exchangeI(-1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexi2)
# 		g1 = g1 + exchangeI(1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexi2)
# 		g2 = secondObjectiv+exchangeR(-1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexr1)
# 		g2 = g2 + exchangeR(1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexr1)
# 		g2 = secondObjectiv+exchangeR(-1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexr2)
# 		g2 = g2 + exchangeR(1,visitors,rooms,hoursNumber,timeSlot,room1,room2,indexr2)
# 		if g1 < g2
# 			x[indexi1][indexr1][room1][timeSlot]=0
# 			x[indexi2][indexr1][room1][timeSlot]=1
# 			x[indexi1][indexr2][room2][timeSlot]=1
# 			x[indexi2][indexr2][room2][timeSlot]=0
# 			g = g1
# 		else
# 			x[indexi1][indexr1][room1][timeSlot]=0
# 			x[indexi2][indexr2][room2][timeSlot]=0
# 			x[indexi1][indexr2][room1][timeSlot]=1
# 			x[indexi2][indexr1][room2][timeSlot]=1
# 			g = g2
# 		end
# 	end
# end
#
# function delayMeeting(visitors, rooms, hoursNumber, room, timeSlot1, timeSlot2)
# RDV = RoomSchedule[room][timeSlot1]
# indexi = RDV[1]
# indexr = RDV[2]
# if IndstrialistSchedule[indexi][timeSlot2] == 0 && ResearcherSchedule[indexr][timeSlot2] == 0
# 	available = false
# 	s = 1
# 	while !available && s <= rooms[2]
# 		if RoomSchedule[s][timeSlot2][1] == 0
# 			available = true
# 			indexs = s
# 		end
# 		s = s+1
# 	end
# 	if available
#
# 		x[indexi][indexr][room][timeSlot1]=0
# 		x[indexi][indexr][indexs][timeSlot2]=1
#
# 		secondObjectiv = g
# 		secondObjectiv = secondObjectiv + exchangeI(-1, visitors, rooms, hoursNumber, timeSlot1, room, 0, indexi)
# 		secondObjectiv = secondObjectiv + exchangeI(-1, visitors, rooms, hoursNumber, timeSlot2, 0, indexs, indexi)
# 		secondObjectiv = secondObjectiv + exchangeR(-1, visitors, rooms, hoursNumber, timeSlot1, room, 0, indexr)
# 		secondObjectiv = secondObjectiv + exchangeR(-1, visitors, rooms, hoursNumber, timeSlot2, 0, indexs, indexr)
# 		secondObjectiv = secondObjectiv + exchangeI(1, visitors, rooms, hoursNumber, timeSlot1, room, 0, indexi)
# 		secondObjectiv = secondObjectiv + exchangeI(1, visitors, rooms, hoursNumber, timeSlot2, 0, indexs, indexi)
# 		secondObjectiv = secondObjectiv + exchangeR(1, visitors, rooms, hoursNumber, timeSlot1, room, 0, indexr)
# 		secondObjectiv = secondObjectiv + exchangeR(1, visitors, rooms, hoursNumber, timeSlot2, 0, indexs, indexr)
# 		g = secondObjectiv
# 	end
# end

ext: jl


text:
using Game2048
using Game2048Core
using Game2048Core: value, count0
using Dates

mutable struct MctsTreeNode
    parent::Union{Nothing,MctsTreeNode}
    nvisit::Int32
    totreward::Float32
    state::Bitboard
    children_prior::Vector{Float32} # probability of children
    children::Vector{MctsTreeNode}
    children_are_moves::Bool # Its children are moves instead of add tile
    policy::Function

    function MctsTreeNode(; state, policy, children_are_moves, parent)
        new(parent, 0, 0, state, [], [], children_are_moves, policy)
    end
end

function Base.show(io::IO, tree::MctsTreeNode)
    show(io, tree.nvisit)
end

function _display(tree::MctsTreeNode, indent)
    if length(indent) >= 4
        return
    end
    display(tree.state)
    for (prior, child) in zip(tree.children_prior, tree.children)
        status = string((totreward = child.totreward, prior = prior, nvisit = child.nvisit))
        display(indent * status)
        _display(child, indent * "    ")
    end
end

function Base.display(tree::MctsTreeNode)
    _display(tree, "")
end

# do I need to bubble up the reward? Since we can choose the action with the most visits.
# propate the reward (bubble up)
function bubble_up_reward!(child_node::MctsTreeNode)
    if isnothing(child_node.parent) || child_node.parent.children_are_moves
        if all(child_node.children_prior .== 0)
            # if no more moves
            child_node.nvisit += 1
            child_node.totreward = value(child_node.state)
        else
            # did not receive a value so I am updating it
            # it's children are all addtiles, so you don't get to choose but can only average them
            maxreward = mapreduce(node -> node.totreward, max, child_node.children)
            totvisits = mapreduce(node -> node.nvisit, +, child_node.children)

            # this function should update the reward and visit count
            child_node.totreward = maxreward
            child_node.nvisit = totvisits
        end

        if isnothing(child_node.parent)
            return
        else
            bubble_up_reward!(child_node.parent)
        end
    else
        # did not receive a value so I am updating it
        # its children are all `add_tile`s, so you don't get to choose but can only average them
        rewards = map(node -> node.totreward, child_node.children)
        visits = map(node -> node.nvisit, child_node.children)

        # this function should update the reward and visit count
        child_node.totreward = sum(rewards .* visits) / sum(visits)
        child_node.nvisit = sum(visits)

        # upto here
        bubble_up_reward!(child_node.parent)
    end
end

"""
Select a move randomly according to policy distribution

Return:

    the same board and left move if no move is possible
"""
function select_move(bboard::Bitboard, probs::Vector{Float32})
    # mostly copied from play_one_move_with_policy
    # select a direction to play
    tot_prob = sum(probs)
    while !(tot_prob ≈ 0.0)
        for i in 1:4
            if rand() < probs[i] / tot_prob
                # select this move
                newboard = move(bboard, DIRS[i])
                if newboard == bboard # if nothing changed then the move is not possible
                    # this move will never be selected again
                    probs[i] = zero(Float32)
                else
                    return newboard, DIRS[i]
                end
            end
            tot_prob -= probs[i]
        end

        # if no move got selected. this could just be unlucky; just keep going
        tot_prob = sum(probs)
    end
    return bboard, left
end

const C = 2 # exploration constant

"""
    Internal function do a single simulation given a bitboard
"""
function _rollout(state::Bitboard, policy)
    # a rollout will just play the game at random
    # select a children to do rollout
    end_state = play_game_with_policy(policy, state)

    return value(end_state)
end

"""
    rollout(tree::MctsTreeNode)

Do one rolll out of the game using the policy
"""
function rollout(tree::MctsTreeNode)
    if !tree.children_are_moves
        tmp_state = add_tile(tree.state)
    end
    _rollout(tree.state, tree.policy)
end

# use this to decide which child node to do a rollout
function select_child_node(tree::MctsTreeNode)
    if tree.children_are_moves
        if all(tree.children_prior .== 0)
            # if no move is posible return itself as it's not possible go any further
            return tree
        end
        best_child_val = 0.0
        best_child = tree.children[1]
        for (prior_prob, child) in zip(tree.children_prior, tree.children)
            # Note: due to expand_node!, each child must have had a visit and rollout already
            child_val = tree.totreward + prior_prob / child.nvisit
            if best_child_val < child_val
                best_child_val = child_val
                best_child = child
            end
        end
        return best_child
    else
        best_child_val = 0.0
        if length(tree.children) == 0
            println("meh")
        end
        best_child = tree.children[1]

        for (prior_prob, child) in zip(tree.children_prior, tree.children)
            # Note: due to expand_node!, each child must have had a visit and rollout already

            # for notes that's all about add tile; we can't choose by the reward hence the exploration
            # is purely guide by probabilities
            child_val = prior_prob / child.nvisit
            if best_child_val < child_val
                best_child_val = child_val
                best_child = child
            end
        end
        return best_child
    end
end

function expand_node!(tree::MctsTreeNode)
    @assert length(tree.children) == 0
    if tree.children_are_moves
        tree.children = Vector{MctsTreeNode}(undef, 4)
        tree.children_prior = tree.policy(tree.state)

        rescale_prior = false

        # expand the node
        for (i, dir) in enumerate(DIRS)
            tmp_state = move(tree.state, dir)
            # if the move has no effect, skip it
            if tmp_state == tree.state
                # set it's prior to 0
                tree.children_prior[i] = 0

                rescale_prior = true # rescale the priors so they sum to 1
            end

            tmp_child = MctsTreeNode(;
                state = tmp_state,
                policy = tree.policy,
                children_are_moves = false,
                parent = tree)

            tree.children[Int(dir)+1] = tmp_child
            tmp_child.totreward += rollout(tmp_child)
            # could've set1 it to 1 as there is not reason why the child just created would be anything
            # other than 0
            tmp_child.nvisit += 1
        end

        if rescale_prior
            if all(tree.children_prior .== 0)
                # this means the game has ended as no move is possible
                return tree
            end
            tree.children_prior .= tree.children_prior ./ sum(tree.children_prior)
        end

        return tree
    else
        # find all the empty spots
        num_of_possible_pos = count0(tree.state)

        tree.children_prior = vcat(
            fill(1 / num_of_possible_pos * 0.9, num_of_possible_pos),
            fill(1 / num_of_possible_pos * 0.1, num_of_possible_pos),
        )

        # make sure it sums to 1
        tree.children_prior .= tree.children_prior ./ sum(tree.children_prior)

        resize!(tree.children, num_of_possible_pos * 2)

        for i in 1:num_of_possible_pos
            # add the possible of a 2-tile being generated
            tmp_state = add_tile(tree.state, i, 1)
            tmp_child = MctsTreeNode(;
                state = tmp_state,
                policy = tree.policy,
                children_are_moves = true,
                parent = tree)

            tree.children[i] = tmp_child

            tmp_child.nvisit += 1
            tmp_child.totreward += rollout(tmp_child)

            # add the possible of a 4-tile being generated
            tmp_state = add_tile(tree.state, i, 2)

            tmp_child = MctsTreeNode(;
                state = tmp_state,
                policy = tree.policy,
                children_are_moves = true,
                parent = tree)

            tree.children[i+num_of_possible_pos] = tmp_child

            tmp_child.nvisit += 1
            tmp_child.totreward += rollout(tree.children[i])
        end
        return tree
    end
end

# every `expand_node!` does one roll out of its children
# so every node has had at least one visit
function select_node_to_expand(tree)
    # traverse the tree and find a node to expand
    selected_child = select_child_node(tree)
    # println("this is the next child:")
    # display(selected_child.state)

    # if it does not have children then visit is just 1
    if (selected_child.nvisit == 1) || all(==(0), selected_child.children_prior)
        # println("got here")
        # display(selected_child)
        return selected_child
    else
        # if the child node has been visited do it recursively
        select_node_to_expand(selected_child)
    end
end

function inittree(policy::Function)
    state = initbboard()
    # function monte_carlo_tree_search(bboard, policy)
    # probs = policy(state)

    # create the base node
    tree = MctsTreeNode(; state, policy, children_are_moves = true, parent = nothing)

    # always expand the base node
    # when you expand, the children is populated
    # and one rollout for each child is done
    # when I expand a node should I add to its visit count?
    expand_node!(tree)
end

function displaytop(tree)
    for (prior, child) in zip(tree.children_prior, tree.children)
        status = string((totreward = child.totreward, prior = prior, nvisit = child.nvisit))
        display(status)
    end
end

function mcts!(tree)
    ## now select a child node
    child_node = select_node_to_expand(tree)
    # expand it; all its children have now done a roll out
    if length(child_node.children) == 0
        expand_node!(child_node)
    end
    bubble_up_reward!(child_node)
    tree
end

function mcts!(tree, n)
    for _ in 1:n
        mcts!(tree)
    end
    tree
end

function maxdepth(tree, atdepth = 0)
    if length(tree.children) == 0
        return atdepth
    else
        return atdepth + 1 + maximum(maxdepth, tree.children)
    end
end

function play_game_via_mcts_w_timer(policy; ms = 10, verbose = false)
    tree = inittree(policy)

    while true
        # select a move by playing mcts
        t = Dates.now()
        while (Dates.now() - t).value < ms # within one second
            mcts!(tree)
        end

        # select the best move
        children_and_priors = zip(tree.children, tree.children_prior)
        best_child_idx = argmax(map(((c, prior),) -> c.totreward * (prior != 0), children_and_priors))

        dir = DIRS[best_child_idx]

        pre_add_tile_state = move(tree.state, dir)
        if count0(pre_add_tile_state) == 0
            if verbose
                println("game ended")
                display(pre_add_tile_state)
            end
            return pre_add_tile_state
        end
        grand_child_state = add_tile(pre_add_tile_state)

        grand_children = tree.children[best_child_idx].children

        gc_idx = findfirst(child -> child.state == grand_child_state, grand_children)

        # if isnothing(gc_idx)
        #     @infiltrate
        # end

        # look the grandchild
        tree = tree.children[best_child_idx].children[gc_idx]
        tree.parent = nothing

        if verbose
            display(grand_child_state)
            println("the maxdepth is now $(maxdepth(tree))")
        end
    end
end

function play_game_via_mcts_w_tries(policy; n = 10, verbose = false)
    tree = inittree(policy)

    while true
        # select a move by playing mcts
        mcts!(tree, n)

        # select the best move
        children_and_priors = zip(tree.children, tree.children_prior)
        best_child_idx = argmax(map(((c, prior),) -> c.totreward * (prior != 0), children_and_priors))

        dir = DIRS[best_child_idx]

        pre_add_tile_state = move(tree.state, dir)
        if count0(pre_add_tile_state) == 0
            if verbose
                println("game ended")
                display(pre_add_tile_state)
            end
            return pre_add_tile_state
        end
        grand_child_state = add_tile(pre_add_tile_state)

        grand_children = tree.children[best_child_idx].children

        gc_idx = findfirst(child -> child.state == grand_child_state, grand_children)

        # if isnothing(gc_idx)
        #     @infiltrate
        # end

        # look the grandchild
        tree = tree.children[best_child_idx].children[gc_idx]
        tree.parent = nothing

        if verbose
            display(grand_child_state)
            println("the maxdepth is now $(maxdepth(tree))")
        end
    end
end


text:
# ---
# title: 932. Beautiful Array
# id: problem932
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Divide and Conquer
# link: <https://leetcode.com/problems/beautiful-array/description/>
# hidden: true
# ---
#
# For some fixed `N`, an array `A` is _beautiful_ if it is a permutation of the
# integers `1, 2, ..., N`, such that:
#
# For every `i < j`, there is **no**  `k` with `i < k < j` such that `A[k] * 2 =
# A[i] + A[j]`.
#
# Given `N`, return **any** beautiful array `A`.  (It is guaranteed that one
# exists.)
#
#
#
# **Example 1:**
#
#
#
#     Input: 4
#     Output: [2,1,4,3]
#
#
# **Example 2:**
#
#
#
#     Input: 5
#     Output: [3,1,2,5,4]
#
#
#
# **Note:**
#
#   * `1 <= N <= 1000`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
"""
    gradient(f, p::AbstractParticles)

Calculate the gradient of `f` in `p`. This corresponds to a smoothed finite-difference approximation where the smoothing kernel is given by the distribution of `p`.
Return mean and std.
"""
function gradient(f,p::MonteCarloMeasurements.AbstractParticles)
    r = 2(p\f(p))
    mean(r), std(r)
end
function gradient(f,p::Union{Integer, AbstractFloat})
    p = p  ± 0.000001
    r = 2(p\f(p))
    mean(r)
end

function gradient(f::Function,p::MonteCarloMeasurements.MvParticles)
    r = (p-mean(p))\(f(p) - f(mean(p)))
end

"""
    jacobian(f::Function, p::Vector{<:AbstractParticles})

Calculate the Jacobian of `f` in `p`. This corresponds to a smoothed finite-difference approximation where the smoothing kernel is given by the distribution of `p`.
"""
function jacobian(f::Function,p::MonteCarloMeasurements.MvParticles)
    r = (p-mean(p))\(f(p) - f(mean(p)))
end


text:
#=
Goal : compare neuron-like input with fake connection input
type, make sure they match perfectly
=#

push!(LOAD_PATH, abspath(@__DIR__,".."))

using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark) #; plotlyjs();
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using BenchmarkTools
using ProgressMeter
import FileIO: save # save raster png

using InvertedIndices

using Random; Random.seed!(0)

##

# define E population of conductance based EIF neurons
const dt = 0.1E-3
const Ttot = 1.0
const Ne = 5

myτe = 20E-3 # seconds
τrefr= 15E-3 # refractoriness
vth_e = 20.   # mV
vthexp = -50.0 # actual threshold for spike-generation
eifslope = 2.0
v_rest_e = -70.0
v_rev_e = 0.0
v_leak_e = v_rest_e
v_reset_e = v_rest_e
Cap = 300.0 #capacitance mF

myτe_ker = 20E-3 # synaptic kernel
myτe_ker_plus = 8E-3
myτe_ker_minus = 3E-3

# synaptic_kernel = S.SKExp(myτe_ker)
synaptic_kernel = S.SKExpDiff(myτe_ker_plus,myτe_ker_minus)

nt_e,ps_e = let sgen = S.SpikeGenEIF(vthexp,eifslope)
  nt = S.NTLIFConductance(synaptic_kernel,sgen,myτe,Cap,
    vth_e,v_reset_e,v_rest_e,τrefr,v_rev_e)
  ps = S.PSLIFConductance(nt,Ne)
  (nt,ps)
end

## define input trains
const in_rate = 30.0
const w_in_vec = rand(Uniform(80.0,200.0),Ne)
onetrain(rate,Ttot) = sort(rand(Uniform(0,Ttot),round(Integer,rate*Ttot)))
trains = [ onetrain(in_rate,Ttot) for _ in 1:Ne ]
spkgen = S.SGTrains(trains)

nt_in_type1 = S.NTInputConductance(spkgen,synaptic_kernel,v_rev_e)
ps_in_type1 = S.PSInputConductance(nt_in_type1,Ne)

nt_in_type2 = S.NTInputConductance(spkgen,synaptic_kernel,v_rev_e)
ps_in_type2 = S.PSInputPoissonConductance(nt_in_type2,w_in_vec)

## define connection for type 1 input
conn_e_in = let wmat=sparse(Diagonal(w_in_vec))
  S.ConnGeneralIF2(wmat)
end

## connect !
pop_et1 = S.Population(ps_e,(conn_e_in,ps_in_type1))
pop_in1 = S.Population(ps_in_type1,(S.FakeConnection(),ps_in_type1))
pop_et2 = S.Population(ps_e,(S.FakeConnection(),ps_in_type2))

ntw1 = S.RecurrentNetwork(dt,pop_et1,pop_in1)
ntw2 = S.RecurrentNetwork(dt,pop_et2)

## initial conditions
initial_states = rand(Uniform(v_reset_e,-49.0),Ne)

##
const krec = 1
rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot)
rec_spikes = S.RecSpikes(ps_e,50.0,Ttot)

## Run
times = (0:dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_state_e,rec_spikes])
S.reset!.([ps_e,ps_in_type1,ps_in_type2])
S.reset!(conn_e_in)
# initial conditions
ps_e.state_now .= initial_states

testinput1=Vector{Float64}[]
@time begin
  @showprogress 1.0 "network simulation " for (k,t) in enumerate(times)
    rec_state_e(t,k,ntw1)
    rec_spikes(t,k,ntw1)
    push!(testinput1,ps_e.input)
    S.dynamics_step!(t,ntw1)
  end
end
S.add_fake_spikes!(vth_e,rec_state_e,rec_spikes)

testinput1=hcat(testinput1...)

## plot the state for one neuron

_ = let neu = 4,
  plt=plot()
  plot!(plt,rec_state_e.times,rec_state_e.state_now[neu,:];linewidth=2,leg=false,
    xlims=(0,1))
end

spkdic=S.get_spiketimes_dictionary(rec_spikes)

## now try the other one

rec_state_e2 = S.RecStateNow(ps_e,krec,dt,Ttot)
rec_spikes2 = S.RecSpikes(ps_e,50.0,Ttot)

## Run
# clean up
S.reset!.([rec_state_e2,rec_spikes2])
S.reset!.([ps_e,ps_in_type1,ps_in_type2])
# initial conditions
ps_e.state_now .= initial_states

testinput2=Vector{Float64}[]

@time begin
  @showprogress 1.0 "network simulation " for (k,t) in enumerate(times)
    rec_state_e2(t,k,ntw2)
    rec_spikes2(t,k,ntw2)
    push!(testinput2,ps_e.input)
    S.dynamics_step!(t,ntw2)
  end
end
S.add_fake_spikes!(vth_e,rec_state_e2,rec_spikes2)
testinput2=hcat(testinput2...)
## plot the state for one neuron

_ = let neu = 4,
  plt=plot()
  plot!(plt,rec_state_e2.times,rec_state_e2.state_now[neu,:];linewidth=2,leg=false)
  plot!(plt,rec_state_e.times,rec_state_e.state_now[neu,:];linewidth=2,leg=false,
    linestyle=:dash)
end

#spkdic=S.get_spiketimes_dictionary(rec_spikes2)


text:
export load_policy

using Pkg.Artifacts

function load_policy(s::String, T = AbstractPolicy)
    if isfile(s)
        RLCore.load(s, T)
    elseif isdir(s)
        load_policy(joinpath(s, "policy.bson"), T)
    else
        dir = ensure_artifact_installed(s, find_artifacts_toml(@__DIR__))
        return load_policy(dir, T)
    end
end


text:
export AdvancedHMC

using Random:
    AbstractRNG,
    GLOBAL_RNG
using AdvancedHMC
using AdvancedHMC:
    Hamiltonian,
    AbstractProposal,
    AbstractIntegrator,
    PhasePoint,
    resize,
    sample_init,
    step
import AdvancedHMC: sample

function update_metric(h::Hamiltonian, q::AbstractVector)
    return resize(h, q) # Ensure h.metric has the same dim as q.
end

function make_phasepoint(rng::AbstractRNG, h::Hamiltonian, q::AbstractVector)
    h, t = sample_init(rng, h, q)
    return t.z
end

make_phasepoint(h, q) = make_phasepoint(GLOBAL_RNG, h, q)

position(z::PhasePoint) = z.θ

momentum(z::PhasePoint) = z.r

step_size(i::AbstractIntegrator) = i.ϵ

function sample(rng::AbstractRNG, h::Hamiltonian, τ::AbstractProposal, z::PhasePoint)
    t = step(rng, h, τ, z)
    return t.z, t.stat
end

sample(h, τ, z) = sample(GLOBAL_RNG, h, τ, z)


text:
function toCSP(data::BPPData)
    weights = [ i for i in unique(data.weights) ]
    demands = [ count(==(i), data.weights) for i in unique(data.weights) ]
    return CSPData(data.name, data.capacity, weights, demands, data.lb, data.ub)
end

function toBPP(data::CSPData)
    weights = []
    for i in 1:length(data.weights)
        weights = vcat(weights, repeat([ data.weights[i] ], data.demands[i]))
    end
    return BPPData(data.name, data.capacity, weights, data.lb, data.ub)
end

function loadBounds(name::String)
    file_name = joinpath(data_path, "Bounds.zip")
    file = ZipFile.Reader(file_name)
    values = split(read(file.files[1], String))

    index = findfirst(isequal(name), values)
    if index !== nothing
        return parse(Int64, values[index + 1]), parse(Int64, values[index + 2])
    else
        return 0, typemax(Int64)
    end
end


text:
module Accuracy

using JLD2

using PanGraph
using PanGraph.Graphs
using Hungarian
using Statistics

# ------------------------------------------------------------------------
# cost computation

δ(q, r; L=100) = mod(q-r, L)
distance(qry, ref; L=100) = min.(δ.(qry, ref'; L=L), δ.(ref', qry; L=L))

function cutoff(position, len, low; match=nothing)
    measurable = copy(position)

    @label loop
    D = mod.(measurable.- measurable', len)
    i = findfirst(0 .< D .< low)
    i !== nothing || @goto endloop

        if match === nothing
            deleteat!(measurable, i.I[2])
        else
            m1 = minimum(distance(measurable[i.I[1]], match; L=len))
            m2 = minimum(distance(measurable[i.I[2]], match; L=len))
            if m1 < m2
                deleteat!(measurable, i.I[2])
            else
                deleteat!(measurable, i.I[1])
            end
        end

    @goto  loop
    @label endloop

    return measurable
end

function nearestbreaks(graph)
    return [
        let
            L = length(sequence(known))
            guess = graph.guess.sequence[name]

            x = cutoff(guess.position,L,100)
            y = cutoff(known.position,L,100; match=x)
            d = distance(x,y; L=L)

            _, cost = hungarian(d)
            cost / min(length(x),length(y))
        end for (name,known) in graph.known.sequence
    ]
end

struct Partition
    dx :: Int64
    b1 :: Graphs.Block
    b2 :: Graphs.Block
end

function breakpoints(known, guess, name)
    L = length(sequence(known))
    #  x1->     x2->     x3->
    #  | node 1 | node 2 |
    bps = [
        (x,j,p.node[i].block) for (j,p) in enumerate([known,guess])
        for (i,x) in enumerate(p.position[1:end-1])
    ]
    sort!(bps, by=first)

    # get the start block
    blocks = Array{Graphs.Block,1}(undef, 2)
    i, j = (bps[1][2] == 1) ? (1, 2) : (2, 1)

    blocks[i] = bps[1][3]
    for k in length(bps):-1:1
        if bps[k][2] == j
            blocks[j] = bps[k][3]
            break
        end
    end

    # edge case: if length(bps) == 2 (the whole genome for both)
    if length(bps) == 2
        return [Partition(L, known.node[1].block, guess.node[1].block)]
    end

    # edge case, if i = 1 and i = 2 correspond to equal break points
    part = Partition[]; i = 2
    if bps[i][1] == bps[i-1][1]
        blocks[bps[i][2]] = bps[i][3]
        i += 1
    end

    while i <= length(bps)
        dx = bps[i][1] - bps[i-1][1]
        push!(part, Partition(dx, blocks[1], blocks[2]))
        if (i <= length(bps)-1) && (bps[i][1] == bps[i+1][1])
            blocks[bps[i+0][2]] = bps[i+0][3]
            blocks[bps[i+1][2]] = bps[i+1][3]
            i += 2
        else
            blocks[bps[i][2]] = bps[i][3]
            i += 1
        end
    end
    # edge case: have to close the circular genome manually
    dx = δ(bps[1][1], bps[end][1]; L=L)
    push!(part, Partition(dx, blocks[1], blocks[2]))

    return part
end

function tiling(breaks)
    map = Dict{Tuple{Graphs.Block,Graphs.Block}, Int64}()
    for bp in breaks
        for p in bp
            key = (p.b1, p.b2)
            if key in keys(map)
                map[key] += p.dx
            else
                map[key] = p.dx
            end
        end
    end
    return map
end

function mutualentropy(graph)
    self  = (breakpoints(known, known, name) for (name, known) in graph.known.sequence)
    cross = (breakpoints(known, graph.guess.sequence[name], name) for (name, known) in graph.known.sequence)

    function entropy(tile)
        len = sum(values(tile))
        return -sum(v/len * log(v/len) for v in values(tile))
    end

    return entropy(tiling(cross)) - entropy(tiling(self))
end

function compare(path)
	graph = (
		known = open(unmarshal,path.known),
		guess = open(unmarshal,path.guess),
	)

    μ = [Graphs.Blocks.diversity(b) for b in values(graph.guess.block)]
    l = [Graphs.Blocks.length(b) for b in values(graph.guess.block)]
    n = mean( length(p.node) for p in values(graph.known.sequence) )
    return (
        filter((l)->l≤1000, nearestbreaks(graph)),
        mutualentropy(graph),
        sum(μ.*l) ./ sum(l),
        n,
        graph
    )
end

# ------------------------------------------------------------------------
# main point of entry

function unpack(message)
    entry = split(message,';')

    return(
        hgt = parse(Float64, entry[1]),
        snp = parse(Float64, entry[2]),
        nit = parse(Int64,   entry[3]),
        known = entry[4],
        guess = entry[5]
    )
end

function usage()
	println("usage: julia script/assay-alignment.jl <input-fifo> <output.jld2>")
	exit(2)
end

if abspath(PROGRAM_FILE) == @__FILE__
	length(ARGS) == 2 || usage()
    pipe = ARGS[1]
    data = ARGS[2]

    jldopen(data, "w") do database; open(pipe) do io
        for msg in eachline(io)
            println("Recieved: ", msg)
            param = unpack(msg)
            group = JLD2.Group(database, "$(param.hgt)/$(param.snp)/$(param.nit)")
            try
                costs, tiles, dists, nblks, input = compare((known=param.known, guess=param.guess))
                group["input"] = nothing #input # NOTE: increases the stored data massively
                group["costs"] = costs
                group["tiles"] = tiles
                group["nblks"] = nblks
                group["dists"] = dists

                rm(param.known); rm(param.guess)
            catch
                println("PROBLEM: ", param.known, " ", param.guess)
                continue # skip the message
            finally
            end
        end
    end; end
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using Distributed
addprocs(16)
@everywhere include("../src/get_sens.jl")
si = LinRange(0.,1.0,10)
s = zeros(4,10)
s[2,:] = si
#s[3,:] = si
get_sens(s)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using HDF5
import ProgressMeter: Progress
using JLD2, FileIO

mutable struct NER
    worddict::Dict
    chardict::Dict
    tagset
    model
end

function NER()
    words = h5read(wordembeds_file, "key")
    worddict = Dict(words[i] => i for i=1:length(words))
    chardict = Dict("UNKNOWN" => 1)
    for word in words
        for c in Vector{Char}(word)
            get!(chardict, string(c), length(chardict)+1)
            get!(chardict, string(uppercase(c)), length(chardict)+1)
        end
    end
    NER(worddict, chardict, BIOES(), nothing)
end

function encode_word(ner::NER, words::Vector{String})
    worddict = ner.worddict
    unkword = worddict["UNKNOWN"]
    ids = map(words) do w
        w = lowercase(w)
        # w = replace(word, r"[0-9]", '0')
        get(worddict, w, unkword)
    end
    Var(ids)
end

function encode_char(ner::NER, words::Vector{String})
    chardict = ner.chardict
    unkchar = chardict["UNKNOWN"]
    batchdims = Int[]
    ids = Int[]
    for w in words
        # w = replace(word, r"[0-9]", '0')
        chars = Vector{Char}(w)
        push!(batchdims, length(chars))
        for c in chars
            push!(ids, get(chardict,string(c),unkchar))
        end
    end
    Var(ids, batchdims)
end

function readdata!(ner::NER, path::String)
    data = NTuple{3,Var}[]
    words, tags = String[], String[]
    lines = open(readlines, path)
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line) || i == length(lines)
            isempty(words) && continue
            w = encode_word(ner, words)
            c = encode_char(ner, words)
            t = Var(encode(ner.tagset,tags))
            push!(data, (w,c,t))
            empty!(words)
            empty!(tags)
        else
            items = split(line, "\t")
            push!(words, String(items[1]))
            push!(tags, String(items[2]))
        end
    end
    data
end

function train(ner::NER, traindata::Vector, testdata::Vector)
    info("#Training examples:\t$(length(traindata))")
    info("#Testing examples:\t$(length(testdata))")
    info("#Words:\t$(length(ner.worddict))")
    info("#Chars:\t$(length(ner.chardict))")
    info("#Tags:\t$(length(ner.tagset.tag2id))")
    testdata = batch(testdata, 100)

    wordembeds = embeddings(h5read(wordembeds_file,"value"))
    charembeds = embeddings(Float32, length(ner.chardict), 20, init_w=Normal(0,0.05))
    ner.model = setup_model(wordembeds, charembeds, length(ner.tagset.tag2id))
    batchsize = 10
    opt = SGD()
    for epoch = 1:50
        println("epoch:\t$epoch")
        opt.rate = 0.0005 * batchsize / sqrt(batchsize) / (1 + 0.05*(epoch-1))
        println("learning rate: $(opt.rate)")

        shuffle!(traindata)
        batches = batch(traindata, batchsize)
        prog = Progress(length(batches))
        loss = 0.0
        for i in 1:length(batches)
            w, c, t = batches[i]
            h = ner.model(w, c)
            y = softmax_crossentropy(t, h)
            loss += sum(y.data)
            params = gradient!(y)
            foreach(opt, params)
            ProgressMeter.next!(prog)
        end
        loss /= length(batches)
        println("Loss:\t$loss")

        # test
        println("Testing...")
        pred = Int[]
        gold = Int[]
        for (w,c,t) in testdata
            y = ner.model(w, c)
            append!(pred, vec(argmax(y.data,1)))
            append!(gold, t.data)
        end
        length(pred) == length(gold) || throw("Length mismatch: $(length(pred)), $(length(gold))")

        ranges_p = decode(ner.tagset, pred)
        ranges_g = decode(ner.tagset, gold)
        fscore(ranges_g, ranges_p)
        println()
    end
end

function fscore{T}(golds::Vector{T}, preds::Vector{T})
    set = intersect(Set(golds), Set(preds))
    count = length(set)
    prec = round(count/length(preds), 5)
    recall = round(count/length(golds), 5)
    fval = round(2*recall*prec/(recall+prec), 5)
    println("Prec:\t$prec")
    println("Recall:\t$recall")
    println("Fscore:\t$fval")
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
# x[i] posicion X del sensor bluetooth i
# y[i] posicion Y del sensor bluetooth i
# r[i] distacia calculada del sensor bluetooth i

function pos = triangulacion(x, y, r)
  N = length(x); # Numero de sensores bluetooth

  D = -0.00680102923817849*r.^3 - 1.04905123190747*r.^2 - 59.2087843354658*r - 1106.35595941215;

  # Vectores con dos veces la diferencia entre una coordenada y el resto
  h_len = int((N - 1) * N / 2);
  hx    = zeros(1, h_len);    # Para el tamaño del vector, es igual a:
  hy    = zeros(1, h_len);    # sumatoria de i = 0 a i = N - 1 de i
  C     = zeros(1, h_len);

  n = 0;
  for j = 1 : N-1
    for i = 1+j : N

      n = n + 1;  # Índice del vector hx, hy
      hx[n] = 2 .* (x[j] - x[i]);
      hy[n] = 2 .* (y[j] - y[i]);

      C[n] = D[i].^2 - D[j].^2 + x[j].^2 - x[i].^2 + y[j].^2 - y[i].^2;

    end
  end

  C   = C';
  H   = [hx' hy'];
  pos = (inv(H' * H) * H') * C # Posición x e y del objeto

  # Pintando figuras
  # TODO: Esto en mi ordenador no va porque PyPlot para Julia
  # usa una versión más moderna de Qt que la que hay instalada por paquetes
  # y eso produce una excepción
  #
  #using PyPlot
  #figure(1);
  #hold on;
  #stem(x, y, '-b');
  #stem(m[1], m[2], '-r', 's');
  #hold off;
end

max_stars_repo_name: Prometheus-ETSIIT/locaviewer

max_issues_repo_name: Prometheus-ETSIIT/locaviewer

max_forks_repo_name: Prometheus-ETSIIT/locaviewer


text:
# # Creating slabs with ASE
#
# ASE is short for the
# [atomistic simulation environment](https://wiki.fysik.dtu.dk/ase/index.html),
# a Python package to simplify the process of setting up, running and
# analysing results from atomistic simulations across different programs.
# Extremely powerful in this respect are the routines this code provides
# for setting up complicated systems (including surface-adsorption scenarios,
# defects, nanotubes, etc).
# See also the [ASE installation instructions](https://wiki.fysik.dtu.dk/ase/install.html).
#
# This example shows how to use ASE to setup a particular gallium arsenide surface
# and run the resulting calculation in DFTK.
# If you are less interested in having access to the full playground of options in DFTK,
# but more interested in performing analysis in ASE itself,
# have a look at [asedftk](https://github.com/mfherbst/asedftk).
# This package provides an ASE-compatible calculator class based on DFTK,
# such that one may write the usual Python scripts against ASE,
# but the calculations are still run in DFTK.
#
# The particular example we consider the (1, 1, 0) GaAs surface separated by vacuum
# with the setup slightly adapted from [^RCW2001].
#
# [^RCW2001]:
#     D. Raczkowski, A. Canning, and L. W. Wang
#     *Thomas-Fermi charge mixing for obtaining self-consistency in density functional calculations*
#     Phys. Rev. B **64**, 121101(R).

# Parameters of the calculation. Since this surface is far from easy to converge,
# we made the problem simpler by choosing a smaller `Ecut` and smaller values
# for `n_GaAs` and `n_vacuum`.
# More interesting settings are `Ecut = 15` and `n_GaAs = n_vacuum = 20`.
#
miller = (1, 1, 0)   # Surface Miller indices
n_GaAs = 2           # Number of GaAs layers
n_vacuum = 4         # Number of vacuum layers
Ecut = 5             # Hartree
kgrid = (4, 4, 1);   # Monkhorst-Pack mesh

# Use ASE to build the structure:
using PyCall

ase_build = pyimport("ase.build")
a = 5.6537  # GaAs lattice parameter in Ångström (because ASE uses Å as length unit)
gaas = ase_build.bulk("GaAs", "zincblende", a=a)
surface = ase_build.surface(gaas, miller, n_GaAs, 0, periodic=true);

# Get the amount of vacuum in Ångström we need to add
d_vacuum = maximum(maximum, surface.cell) / n_GaAs * n_vacuum
surface = ase_build.surface(gaas, miller, n_GaAs, d_vacuum, periodic=true);

# Write an image of the surface and embed it as a nice illustration:
pyimport("ase.io").write("surface.png", surface * (3, 3, 1),
                         rotation="-90x, 30y, -75z")

#md # ```@raw html
#md # <img src="../surface.png" width=500 height=500 />
#md # ```
#nb # <img src="surface.png" width=500 height=500 />

# Use the `load_atoms` and `load_lattice` functions
# to convert to DFTK datastructures.
# These two functions not only support importing ASE atoms into DFTK,
# but a few more third-party datastructures as well.
# Typically the imported `atoms` use a bare Coulomb potential,
# such that appropriate pseudopotentials need to be attached in a post-step:

using DFTK

atoms = load_atoms(surface)
atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional="pbe")) => position
         for (el, position) in atoms]
lattice = load_lattice(surface);

# We model this surface with (quite large a) temperature of 0.01 Hartree
# to ease convergence. Try lowering the SCF convergence tolerance (`tol`
# or the `temperature` to see the full challenge of this system.
model = model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe],
                  temperature=0.001, smearing=DFTK.Smearing.Gaussian())
basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)

scfres = self_consistent_field(basis, tol=1e-4, mixing=KerkerMixing());
#-
scfres.energies


text:
using InteractiveUtils
using Test

using Trixi

@testset "TrixiDebug.jl" begin

versioninfo(verbose=true)

# Find version of CPUSummary and display additional information, see
# https://discourse.julialang.org/t/how-to-find-out-the-version-of-a-package-from-its-module/37755/15
using CPUSummary
using TOML
CPUSummary_project = TOML.parsefile(joinpath(pkgdir(CPUSummary), "Project.toml"))
CPUSummary_version = VersionNumber(CPUSummary_project["version"])
if CPUSummary_version == v"0.1.14"
  @info "CPUSummary" CPUSummary_version isdefined(CPUSummary, :safe_topology_load!) CPUSummary.USE_HWLOC
else
  @info "CPUSummary" CPUSummary_version isdefined(CPUSummary, :safe_topology_load!)
end

@testset "3D DGMulti failures" begin
  @time @testset "elixir_euler_weakform.jl (Hexahedral elements)" begin
    # trixi_include(
    #   joinpath(examples_dir(), "dgmulti_3d", "elixir_euler_weakform.jl"),
    #   element_type=Hex()
    # )
    @info "$(@__FILE__), l. $(@__LINE__): New setup"

    dg = DGMulti(polydeg = 3, element_type = Hex(),
                 surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
                 volume_integral = VolumeIntegralWeakForm())
    @info "$(@__FILE__), l. $(@__LINE__): Constructed dg"

    equations = CompressibleEulerEquations3D(1.4)
    initial_condition = initial_condition_convergence_test
    source_terms = source_terms_convergence_test
    @info "$(@__FILE__), l. $(@__LINE__): Constructed equations"

    # example where we tag two separate boundary segments of the mesh
    top_boundary(x, tol=50*eps()) = abs(x[2] - 1) < tol
    rest_of_boundary(x, tol=50*eps()) = !top_boundary(x, tol)
    is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)

    mesh = DGMultiMesh(dg, cells_per_dimension=(4, 4, 4), is_on_boundary=is_on_boundary)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed mesh"

    boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
    boundary_conditions = (; :top => boundary_condition_convergence_test,
                            :rest => boundary_condition_convergence_test)

    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                        source_terms = source_terms,
                                        boundary_conditions = boundary_conditions)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed semi"

    tspan = (0.0, 0.1)
    ode = semidiscretize(semi, tspan)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed ode"

    u_ode = copy(ode.u0)
    du_ode = zero(u_ode)
    Trixi.rhs!(du_ode, u_ode, semi, first(tspan))
    @info "$(@__FILE__), l. $(@__LINE__): Finished calling rhs!" sum(du_ode)
  end

  @time @testset "elixir_euler_curved.jl (Hex elements, GaussSBP, flux differencing)" begin
    # trixi_include(
    #   joinpath(examples_dir(), "dgmulti_3d", "elixir_euler_curved.jl"),
    #   approximation_type=GaussSBP()
    # )
    @info "$(@__FILE__), l. $(@__LINE__): New setup"

    dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type=GaussSBP(),
                surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
                volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
    @info "$(@__FILE__), l. $(@__LINE__): Constructed dg"

    equations = CompressibleEulerEquations3D(1.4)
    initial_condition = initial_condition_convergence_test
    source_terms = source_terms_convergence_test
    @info "$(@__FILE__), l. $(@__LINE__): Constructed equations"

    # example where we tag two separate boundary segments of the mesh
    top_boundary(x, tol=50*eps()) = abs(x[2] - 1) < tol
    rest_of_boundary(x, tol=50*eps()) = !top_boundary(x, tol)
    is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)

    function mapping(xi, eta, zeta)
      x = xi   + 0.1 * sin(pi * xi) * sin(pi * eta)
      y = eta  + 0.1 * sin(pi * xi) * sin(pi * eta)
      z = zeta + 0.1 * sin(pi * xi) * sin(pi * eta)
      return SVector(x, y, z)
    end
    cells_per_dimension = (4, 4, 4)
    mesh = DGMultiMesh(dg, cells_per_dimension, mapping, is_on_boundary=is_on_boundary)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed mesh"

    boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
    boundary_conditions = (; :top => boundary_condition_convergence_test,
                            :rest => boundary_condition_convergence_test)

    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                        source_terms = source_terms,
                                        boundary_conditions = boundary_conditions)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed semi"

    tspan = (0.0, 0.1)
    ode = semidiscretize(semi, tspan)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed ode"

    u_ode = copy(ode.u0)
    du_ode = zero(u_ode)
    Trixi.rhs!(du_ode, u_ode, semi, first(tspan))
    @info "$(@__FILE__), l. $(@__LINE__): Finished calling rhs!" sum(du_ode)
  end

  @time @testset "elixir_euler_weakform_periodic.jl (Hexahedral elements)" begin
    # trixi_include(
    #   joinpath(examples_dir(), "dgmulti_3d", "elixir_euler_weakform_periodic.jl"),
    #   element_type=Hex()
    # )
    @info "$(@__FILE__), l. $(@__LINE__): New setup"

    dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(),
                surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
                volume_integral = VolumeIntegralWeakForm())
    @info "$(@__FILE__), l. $(@__LINE__): Constructed dg"

    equations = CompressibleEulerEquations3D(1.4)
    initial_condition = initial_condition_convergence_test
    source_terms = source_terms_convergence_test
    @info "$(@__FILE__), l. $(@__LINE__): Constructed equations"

    mesh = DGMultiMesh(dg, cells_per_dimension=(4, 4, 4), periodicity=true)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed mesh"

    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                        source_terms = source_terms)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed semi"

    tspan = (0.0, 0.1)
    ode = semidiscretize(semi, tspan)
    @info "$(@__FILE__), l. $(@__LINE__): Constructed ode"

    u_ode = copy(ode.u0)
    du_ode = zero(u_ode)
    Trixi.rhs!(du_ode, u_ode, semi, first(tspan))
    @info "$(@__FILE__), l. $(@__LINE__): Finished calling rhs!" sum(du_ode)
  end
end

end # "TrixiDebug.jl"


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 6inch)

plot(dataset("MASS", "mammals"), x="Body", y="Brain",
     label=1, Scale.x_log10, Scale.y_log10, Geom.point, Geom.label)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Function to include a table in a file.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export include_pt_in_file

"""
    include_pt_in_file(filename::AbstractString, mark::AbstractString, args...; kwargs...)

Include a table in the file `filename` using the `mark`.

This function will print a table using the arguments `args` and keywords
`kwargs` in the function `pretty_table` (**the IO must not be passed to `args`
here**). Then, it will search inside the file `filename` for the following
section:

    <PrettyTables mark>
    ...
    </PrettyTables>

and will **replace everything between the marks** with the printed table. If the
closing tag is in a separate line, then all characters before it will be kept.
This is important to add comment tags.

If the user wants to also remove the opening and ending tags, then pass the
keyword `remove_tags = true`.

The keyword `tag_append` can be used to pass a string that can be used to add a
text after the opening tag. This is important for HTML where the comments have
openning and closing tags. Thus, if `tag_append = " -->"`, then the following
can be used to add a table into HTML files:

    <!-- <PrettyTables mark> -->
    ...
    <!-- </PrettyTables> -->

By default, this function will copy the original file to `filename_backup`. If
this is not desired, then pass the keyword `backup_file = false` to the
function.
"""
function include_pt_in_file(
    filename::AbstractString,
    mark::AbstractString,
    args...;
    backup_file = true,
    remove_tags = false,
    tag_append::String = "",
    kwargs...
)
    orig = ""

    open(filename, "r") do f
        orig = read(f, String)
    end

    GC.gc()

    # First, print the table into a string.
    io = IOBuffer()
    pretty_table(io, args...; kwargs...)
    str = String(take!(io))

    # Write the output to a temporary file.
    path, io = mktemp()

    if !remove_tags
        r = Regex("(?<=<PrettyTables $mark>$tag_append)(?:.|\n)*?(?=.*</PrettyTables>)")
        write(io, replace(orig, r => "\n$str"))
    else
        r = Regex("<PrettyTables $mark>$tag_append(?:.|\n)*?</PrettyTables>")
        write(io, replace(orig, r => "$str"))
    end
    close(io)

    # Backup the original file if required.
    backup_file && mv(filename, filename * "_backup"; force = true)

    # Copy the temporary file to `filename`.
    #
    # If we user `mv`, then we get some problems related to `libuv` in
    # Windows. This seems related to this issue:
    #
    #   https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278
    cp(path, filename; force = true)

    return nothing
end


text:
import Dance
import DanceWebpack

include("utils.jl")

Dance.start_project("demo")
cd("demo")
project_settings()
add_project_toml()
@test_logs (:info, "Webpack files added to `static` dir in project.\nPlease cd into `static` dir and add dependencies by installing NodeJS and running `npm install`") DanceWebpack.setup()

cd("static")
for line in readlines("package.json")
    for item in [
        ["authors", "[\"Chris \\\"Yoh\\\" Meyers <48554065+yoh-meyers@users.noreply.github.com>\"]"],
        ["name", "Demo_1"],
        ["version", "1.0.1"]
    ]
        if occursin(item[1], line)
            value::String = strip(
                strip(
                    rstrip(split(line, ':')[2], ',')
                ),
                '"'
            )
            @test value==item[2]
        end
    end
end

delete_project()


text:
module Spikes

using MAT, PyPlot, HDF5,  Images, StatsBase, PyCall, HypothesisTests, Distributions

@pyimport matplotlib.collections as C

include("types.jl")
include("stats.jl")
include("rates.jl")
include("tuning.jl")
#include("dimensionreduce.jl")
#include("jpsth.jl")
include("neuron_dropping.jl")
include("information.jl")
include("classifiers.jl")
include("kinematics.jl")

end


text:
using Test
using Random

using StatsModelsTools
using StatsModels

@testset "Test Polynomial Formulas" begin

    data = (
            x=float.(1:10),
            y=float.(1:10).^3,
           )

    f = @formula(y ~ poly(x, 3))
    sch = apply_schema(f, schema(data))
    @inferred modelcols(sch.rhs, data)

    @test poly(data.x, 3) == data.x.^3
    @test modelcols(sch.rhs, data) == [data.x data.x.^2 data.x.^3]

end

@testset "Eval RMS Formulas" begin

    @test_throws Exception rcs([1,2,3], 0)
    @test_throws Exception rcs([1,2,3], 1)

    data = (
            x=(float.(1:10)),
            y=float.(1:10),
            z=collect(1:10),
           )

    xs = collect(0:10)
    @test rms_generate_knots(float.(0:100), 3) == [10, 50, 90]
    @test rms_generate_knots((0:100), 3) == [10, 50, 90]
    @test rms_generate_knots(vcat(float.(0:10), [missing]), 3) == [1, 5, 9]
    @test rms_generate_knots(xs, 3) == [1, 5, 9]
    @test rms_generate_knots(reverse(xs) |> collect, 3) == [1, 5, 9]
    @test rms_generate_knots(shuffle(xs) |> collect, 3) == [1, 5, 9]

    @test rcs(xs, [1, 5, 9]) == [0.0  0.0;
                                 1.0  0.0;
                                 2.0  0.015625;
                                 3.0  0.125;
                                 4.0  0.421875;
                                 5.0  1.0;
                                 6.0  1.921875;
                                 7.0  3.125;
                                 8.0  4.515625;
                                 9.0  6.0;
                                 10.0 7.5]

    xs_with_missing = [xs; missing]
    @inferred rms_generate_knots(xs, 3)
    @inferred rms_generate_knots(xs_with_missing, 3)
    @inferred rcs(xs_with_missing, [1, 5, 9])
    @test rms_generate_knots(xs_with_missing, 3) == [1, 5, 9]

    knots = [1, 5, 9]
    result = rcs(xs, knots);
    result_with_missing = rcs(xs_with_missing, knots)
    @test all(result_with_missing .=== [result; [missing missing]])

    f = @formula(y ~ rcs(x, 3))
    sch = apply_schema(f, schema(data))
    @inferred modelcols(sch.rhs, data)
    @test coefnames(sch)[2] == ["rcs(x, 3)[1]", "rcs(x, 3)[2]"]

end


text:
module GatekeepHelperFlagSummitGem

using ..Ahorn, Maple

# Thanks to Communal Helper for most of this implementation

@mapdef Entity "GatekeepHelper/FlagSummitGem" FlagSummitGem(
    x::Integer,
    y::Integer,
    index::Integer=0,
    sprite::String="",
    flag::String="",
    particleColor::String="",
)

const placements = Ahorn.PlacementDict(
    "Flag Summit Gem (Gatekeep Helper)" => Ahorn.EntityPlacement(
        FlagSummitGem,
    ),
)

const sprites = ["collectables/summitgems/$i/gem00" for i in 0:7]

function getSprite(index)
    if index > length(sprites)
        return sprites[end]
    end
    return sprites[index]
end

# positive numbers only
function getClampedIndex(entity::FlagSummitGem)
    index = Int(get(entity.data, "index", 0))
    if index < 0
        index = 0
    end
    entity.data["index"] = index
    return index
end

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

    return Ahorn.getSpriteRectangle(getSprite(getClampedIndex(entity) + 1), x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagSummitGem) =
    Ahorn.drawSprite(ctx, getSprite(getClampedIndex(entity) + 1), 0, 0)
end


text:
logscale(params::AbstractVector{<:Real}) = /(extrema(params)...) < 2e-2 && minimum(params) > eps()

logscale(params) = false

@recipe function plot(ho::Hyperoptimizer)
    N = length(ho.params)
    layout --> N
    for i = 1:N
        params = getindex.(ho.history, i)
        perm = sortperm(params)
        yguide --> "Function value"
        seriestype --> :scatter
        @series begin
            xguide --> ho.params[i]
            subplot --> i
            label --> "Sampled points"
            legend --> false
            logscale(params) && (xscale --> :log10)
            logscale(ho.results) && (yscale --> :log10)
            params[perm], ho.results[perm]
        end

    end
end

# robust_acq(model, x::Vector{<:Real}) = robust_acq(model, [[x] for x in x]')
robust_acq(model, x::AbstractVector{<:Real}) = robust_acq(model, reshape(x,1,:))
tovec(x::Number) = [x]
tovec(x::AbstractVector) = x

function robust_acq(s, x=s.model.x)
    acqfunc = BayesianOptimization.acquisitionfunction(BayesianOptimization.ExpectedImprovement(maximum(s.model.y)), s.model)
    map(1:size(x,2)) do i
        try
            return acqfunc(tovec(x[:,i]))
        catch ex
            @warn("BayesianOptimization failed, error: ", ex); return -Inf
        end
    end
end

@recipe function plot(s::GPSampler)
    model = s.model
    ndims = size(model.x,1)
    ms, var = GaussianProcesses.predict_f(model, model.x)
    sig = sqrt.(var)
    layout := (ndims == 1 ? 2 : ndims)
    seriestype := :scatter
    xvals = from_logspace([model.x[i,:] for i in 1:ndims], s.logdims)
    xrange = [extrema(s.candidates[i]) for i in 1:ndims]
    for i = 1:ndims
        logscale([model.y[:]; Int(s.sense)*ms]) && (yscale --> :log10)
        @series begin
            subplot := i
            label := "Model function"
            yerror := sqrt.(var)
            xvals[i], Int(s.sense)*ms
        end
        @series begin
            subplot := i
            label := "Obs"
            xlims := xrange[i]
            markershape := :cross
            xvals[i], Int(s.sense)*model.y[:]
        end
    end
    if ndims == 1
        @series begin
            xa = LinRange(extrema(s.candidates[1])...,100)
            aqy = robust_acq(s, xa)
            title -> "Expected improvement"
            subplot := 2
            lab="Aq function"
            xa, aqy
        end
    end
end

# function extend_dim(x::AbstractMatrix,d,n=200)
#     new_dim = LinRange(extrema(x[d,:])..., n)
#     new_x = repeat(x,1,n)
#     new_x[d,:] = repeat(new_dim, 1, size(x,2))
#     new_x
# end
#
# function extend_dim(x::AbstractVector{<:AbstractVector},d,n=200)
#     new_dim = LinRange(extrema(x[d])..., n)
#     new_x = repeat(x,1,n)
#     new_x[d,:] = repeat(new_dim, 1, size(x,2))
#     new_x
# end


text:
# # Getting Started
#
# Here we document how to do some basic stuff, including learning and inference in kernel parameters, with Stheno.jl. To do this, we that makes use of a variety of packages from the Julia ecosystem. In particular, we'll make use of
# - [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl) to perform Bayesian inference in our model parameters,
# - [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) for point-estimates of our model parameters,
# - [ParameterHandling.jl](https://github.com/invenia/ParameterHandling.jl) to make it easy to work with our model's parameters, and to ensure that it plays nicely with Optim and AdvancedHMC,
# - [Zygote.jl](https://github.com/FluxML/Zygote.jl/) to compute gradients.

# This guide assumes that you know roughly what's going on conceptually with GPs. If you're
# new to Gaussian processes, an excellent introduction can be found in either
# [this video lecture](http://videolectures.net/gpip06_mackay_gpb/) or
# [this one](https://www.youtube.com/watch?v=92-98SYOdlY).

# ## Exact Inference in a GP in 2 Minutes
#
# This is only a slightly more interesting version of the first example on the README.
# It's slightly more interesting in that we give the kernels some learnable parameters.

# Import the packages we'll need for this bit of the demo.
using AbstractGPs
using LinearAlgebra
using Stheno
using Plots

# Short length-scale and small variance.
l1 = 0.4
s1 = 0.2

# Long length-scale and larger variance.
l2 = 5.0
s2 = 1.0

# Specify a GaussianProcessProbabilisticProgramme object, which is itself a GP
# built from other GPs.
f = @gppp let
    f1 = s1 * stretch(GP(Matern52Kernel()), 1 / l1)
    f2 = s2 * stretch(GP(SEKernel()), 1 / l2)
    f3 = f1 + f2
end;

# Generate a sample from f3, one of the processes in f, at some random input locations.
# Add some iid observation noise, with zero-mean and variance 0.02.
const x = GPPPInput(:f3, collect(range(-5.0, 5.0; length=100)));
σ²_n = 0.02;
fx = f(x, σ²_n);
const y = rand(fx);

# Compute the log marginal likelihood of this observation, just because we can.
logpdf(fx, y)

# `fx` should be thought of as "`f` at `x`", and is just as a multivariate Normal
# distribution, with zero mean and covariance matrix
heatmap(cov(f, x) + σ²_n * I)

# As such samples can be drawn from it, and the log probability any particular value under
# it can be computed, in the same way that you would an `MvNormal` from
# [Distributions.jl](https://github.com/JuliaStats/Distributions.jl).

# We can visualise `x` and `y` with [Plots.jl](https://github.com/JuliaPlots/Plots.jl)

plt = plot();
scatter!(plt, x.x, y; color=:red, label="");
display(plt)

# It's straightforward to compute the posterior over `f`:
f_posterior = posterior(fx, y);
# `f_posterior` is another GP, the posterior over `f` given noisy observations `y` at inputs `x`.

# The plotting recipes from AbstractGPs.jl can be utilised to easily print whichever
# component of the GPPP is needed.
x_plot = range(-7.0, 7.0; length=1000);
xp = GPPPInput(:f3, x_plot);
plot!(
    plt, x_plot, f_posterior(xp);
    ribbon_scale=3, label="", color=:blue, fillalpha=0.2, linewidth=2,
)
plot!(
    plt, x_plot, rand(f_posterior(xp, 1e-9), 10);
    samples=10, markersize=1, alpha=0.3, label="", color=:blue,
);
plt

# So you've built a simple GP probabilistic programme, performed inference in it, and
# looked at the posterior.
# We've only looked at one component of it though -- we could look at others.
# Consider `f2`:
xp2 = GPPPInput(:f2, x_plot);
plot!(
    plt, x_plot, f_posterior(xp2);
    ribbon_scale=3, label="", color=:red, fillalpha=0.2, linewidth=2,
);
plot!(
    plt, x_plot, rand(f_posterior(xp2, 1e-9), 10);
    alpha=0.3, label="", color=:red,
);
plt

# ## Fit a GP with NelderMead in 2 Minutes

# Stheno.jl is slightly unusual in that it declines to provide a `fit` or `train` function. Why is this? In short, because there's really no need -- the ecosystem now contains everything that is needed to easily do this yourself. By declining to insist on an interface, Stheno.jl is able to interact with a wide array of tools, that you can use in whichever way you please.

# Optim requires that you provide an objective function with a single `Vector{<:Real}` parameter for most of its optimisers.
# We'll use ParameterHandling.jl to build one of these in a way that doesn't involve manually writing code to convert between a structured, human-readable, representation of our parameters (in a `NamedTuple`) and a `Vector{Float64}`.

# First, we'll put the model from before into a function:
function build_model(θ::NamedTuple)
    return @gppp let
        f1 = θ.s1 * stretch(GP(SEKernel()), 1 / θ.l1)
        f2 = θ.s2 * stretch(GP(SEKernel()), 1 / θ.l2)
        f3 = f1 + f2
    end
end

# We've assumed that the parameters will be provided as a `NamedTuple`, so let's build one and check that the model can be constructed:
using ParameterHandling

θ = (
    ## Short length-scale and small variance.
    l1 = positive(0.4),
    s1 = positive(0.2),

    ## Long length-scale and larger variance.
    l2 = positive(5.0),
    s2 = positive(1.0),

    ## Observation noise variance -- we'll be learning this as well.
    s_noise = positive(0.1),
)

# We've used `ParameterHandling.jl`s `positive` constraint to ensure that all of the
# parameters remain positive during optimisation.
# Note that there's no magic here, and `Optim` knows nothing about `positive`.
# Rather, `ParameterHandling` knows how to make sure that `Optim` will optimise the log of
# the parameters which we want to be positive.

# We can make this happen with the following:
using ParameterHandling
using ParameterHandling: value, flatten

θ_flat_init, unflatten = flatten(θ);

# Concrete types used for clarity only.
unpack = value ∘ unflatten;

# We can now easily define a function which accepts the flat form of the parameters, and
# return the negative log marginal likelihood (nlml) of the parameters θ:
function nlml(θ::NamedTuple)
    f = build_model(θ)
    return -logpdf(f(x, θ.s_noise + 1e-6), y)
end

# We can use any gradient-free optimisation technique from
# [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) to find the parameters whose
# negative log marginal likelihood is locally minimal:
using Optim
results = Optim.optimize(
    nlml ∘ unpack,
    θ_flat_init + randn(length(θ_flat_init)),
    NelderMead(),
)
θ_opt = unpack(results.minimizer);
# Note that we just added some noise to the initial values to make the optimisation more
# interesting.

# We can now use this to construct the posterior GP and look at the posterior in comparison
# to the true posterior with the known hyperparameters
f_opt = build_model(θ_opt);
f_posterior_opt = posterior(f_opt(x, θ_opt.s_noise), y);
ms_opt = marginals(f_posterior_opt(xp));
plot!(
    plt, x_plot, mean.(ms_opt);
    ribbon=3std.(ms_opt), label="", color=:green, fillalpha=0.2, linewidth=2,
);
plot!(
    plt, x_plot, rand(f_posterior_opt(xp, 1e-9), 10);
    alpha=0.3, label="", color=:green,
);
plt

# Of course, the exact posterior has not been recovered because the exact hyperparameters
# cannot be expected to be recovered given a finite amount of data over a finite width
# window.

# ## Fit a GP with BFGS in 2 minutes

# The BFGS algorithm is generally the preferred choice when optimising the hyperparameters
# of fairly simple GPs. It requires access to the gradient of our `nlml` function, which
# can be straightforwardly obtained via reverse-mode algorithmic differentiation, which is
# provided by [Zygote.jl](https://github.com/FluxML/Zygote.jl):

using Zygote: gradient

# This will probably take a while to get going as Zygote needs to compile.
results = Optim.optimize(
    nlml ∘ unpack,
    θ->gradient(nlml ∘ unpack, θ)[1],
    θ_flat_init + 0.1 * randn(length(θ_flat_init)),
    BFGS(),
    Optim.Options(
        show_trace=true,
    );
    inplace=false,
)
θ_bfgs = unpack(results.minimizer);

# Once more visualising the results:
f_bfgs = build_model(θ_bfgs);
f_posterior_bfgs = posterior(f_bfgs(x, θ_bfgs.s_noise), y);
ms_bfgs = marginals(f_posterior_bfgs(xp));
plot!(
    plt, x_plot, mean.(ms_bfgs);
    ribbon=3std.(ms_bfgs), label="", color=:orange, fillalpha=0.2, linewidth=2,
);
plot!(
    plt, x_plot, rand(f_posterior_bfgs(xp, 1e-9), 10);
    alpha=0.3, label="", color=:orange,
);
plt

# Notice that the two optimisers produce (almost) indistinguishable results.

# ## Inference with NUTS in 2 minutes

# [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl/) provides a
# state-of-the-art implementation of the No-U-Turns sampler, which we can use to perform
# approximate Bayesian inference in the hyperparameters of the GP.
# This is slightly longer than the previous examples, but it's all set up associated with
# AdvancedHMC, which is literally a copy-paste from that package's README:
using AdvancedHMC, Zygote

# Define the log marginal joint density function and its gradient
ℓπ(θ_flat) = -nlml(unpack(θ_flat)) - 0.5 * sum(abs2, θ_flat)
function ∂ℓπ∂θ(θ_flat)
    lml, back = Zygote.pullback(ℓπ, θ_flat)
    ∂θ_flat = first(back(1.0))
    return lml, ∂θ_flat
end

# Sampling parameter settings
n_samples, n_adapts = 500, 20

# Perturb the initialisation a little bit.
θ0_flat = θ_flat_init + 0.1 * randn(length(θ_flat_init))

# Define metric space, Hamiltonian, sampling method and adaptor
metric = DiagEuclideanMetric(5)
h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
int = Leapfrog(find_good_eps(h, θ0_flat))
prop = NUTS{MultinomialTS, GeneralisedNoUTurn}(int)
adaptor = StanHMCAdaptor(n_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, int.ϵ))

# Perform inference.
samples, stats = sample(h, prop, θ0_flat, n_samples, adaptor, n_adapts; progress=true)

# Inspect posterior distribution over hyperparameters.
hypers = map(unpack, samples);
h_l1 = histogram(getindex.(hypers, :l1); label="l1");
h_l2 = histogram(getindex.(hypers, :l2); label="l2");
h_s1 = histogram(getindex.(hypers, :s1); label="s1");
h_s2 = histogram(getindex.(hypers, :s2); label="s2");
plot(h_l1, h_l2, h_s1, h_s2; layout=(2, 2))

# As expected, the sampler converges to the posterior distribution quickly.
# One could combine this code with that from the previous sections to make predictions under
# the posterior over the hyperparameters.

# ## Conclusion

# So you now know how to fit GPs using Stheno.jl, and to investigate their posterior
# distributions. It's also straightforward to utilise Stheno.jl inside probabilistic
# programming frameworks like Soss.jl and Turing.jl (see examples folder).


text:
using Parameters
import Base: length

using SparseArrays

abstract type AbstractPhase end
export AbstractPhase

abstract type IdealPhase <: AbstractPhase end
export IdealPhase

struct EmptyPhase <: AbstractPhase end
export EmptyPhase

include("Calculators/Ratevec.jl")
include("Calculators/Thermovec.jl")
include("Reaction.jl")

@with_kw struct IdealGas{W<:Tuple,W2,W3} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{ElementaryReaction,1}
    spcdict::Dict{String,Int64}
    stoichmatrix::W2
    Nrp::Array{Float64,1}
    rxnarray::Array{Int64,2}
    veckinetics::W
    veckineticsinds::Array{Int64,1}
    vecthermo::NASAvec
    otherreactions::Array{ElementaryReaction,1}
    electronchange::W3
    reversibility::Array{Bool,1}
    diffusionlimited::Bool = false
end

function IdealGas(species,reactions; name="",diffusionlimited=false)
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,reversible=rxn.reversible,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    therm = getvecthermo(species)
    rxnarray = getreactionindices(species,rxns)
    M,Nrp = getstoichmatrix(rxnarray,species)
    echangevec = getfield.(rxns,:electronchange)
    if all(echangevec .== 0)
        electronchange = nothing
    else
        electronchange = convert(echangevec,Array{Float64,1})
    end
    reversibility = getfield.(rxns,:reversible)
    return IdealGas(species=species,reactions=rxns,name=name,
        spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
        veckineticsinds=posinds, vecthermo=therm, otherreactions=otherrxns, electronchange=electronchange,
        reversibility=reversibility,diffusionlimited=diffusionlimited,)
end
export IdealGas

@with_kw struct IdealDiluteSolution{W<:Tuple,W2,W3} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{ElementaryReaction,1}
    solvent::Solvent
    stoichmatrix::W2
    Nrp::Array{Float64,1}
    rxnarray::Array{Int64,2}
    veckinetics::W
    veckineticsinds::Array{Int64,1}
    vecthermo::NASAvec
    otherreactions::Array{ElementaryReaction,1}
    electronchange::W3
    spcdict::Dict{String,Int64}
    reversibility::Array{Bool,1}
    diffusionlimited::Bool = true
end

function IdealDiluteSolution(species,reactions,solvent; name="",diffusionlimited=true)
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,reversible=rxn.reversible,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    therm = getvecthermo(species)
    rxnarray = getreactionindices(species,rxns)
    M,Nrp = getstoichmatrix(rxnarray,species)
    echangevec = getfield.(rxns,:electronchange)
    if all(echangevec .== 0)
        electronchange = nothing
    else
        electronchange = convert(echangevec,Array{Float64,1})
    end
    reversibility = getfield.(rxns,:reversible)

    return IdealDiluteSolution(species=species,reactions=rxns,solvent=solvent,name=name,
        spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
        veckineticsinds=posinds,vecthermo=therm,otherreactions=otherrxns,electronchange=electronchange,
        reversibility=reversibility,diffusionlimited=diffusionlimited)
end
export IdealDiluteSolution

@with_kw struct IdealSurface{W<:Tuple,W2,W3,W4,W5} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{ElementaryReaction,1}
    stoichmatrix::W2
    Nrp::W5
    rxnarray::Array{Int64,2}
    veckinetics::W
    veckineticsinds::Array{Int64,1}
    vecthermo::W4
    otherreactions::Array{ElementaryReaction,1}
    electronchange::W3
    spcdict::Dict{String,Int64}
    reversibility::Array{Bool,1}
    sitedensity::Float64
    diffusionlimited::Bool = false
end
function IdealSurface(species,reactions,sitedensity;name="",diffusionlimited=false)
    @assert diffusionlimited==false "diffusionlimited=true not supported for IdealSurface"
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,electronchange=rxn.electronchange,reversible=rxn.reversible,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    therm = getvecthermo(species)
    rxnarray = getreactionindices(species,rxns)
    M,Nrp = getstoichmatrix(rxnarray,species)
    echangevec = getfield.(rxns,:electronchange).*F
    if all(echangevec .== 0)
        electronchange = nothing
    else
        electronchange = convert(typeof(Nrp),echangevec)
    end
    reversibility = getfield.(rxns,:reversible)
    return IdealSurface(species=species,reactions=rxns,name=name,
        spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
        veckineticsinds=posinds,vecthermo=therm,otherreactions=otherrxns,electronchange=electronchange,
        reversibility=reversibility,sitedensity=sitedensity,diffusionlimited=diffusionlimited)
end
export IdealSurface

"""
construct the stochiometric matrix for the reactions and the reaction molecule # change
"""
function getstoichmatrix(rxnarray,spcs)
    M = spzeros(size(rxnarray)[2],length(spcs))
    Nrp = zeros(size(rxnarray)[2])
    for i in 1:size(rxnarray)[2]
        n = 0.0
        for (j,ind) in enumerate(rxnarray[:,i])
            if ind == 0
                continue
            elseif j > 3
                M[i,ind] -= 1.0
                n += 1.0
            else
                M[i,ind] += 1.0
                n -= 1.0
            end
        end
        Nrp[i] = n
    end
    return M,Nrp
end

"""
Split the reactions into groups with the same kinetics
"""
function splitreactionsbykinetics(rxns)
    tps = []
    rxnlists = []
    for (i,rxn) in enumerate(rxns)
        typ = getkineticstype(rxn.kinetics)
        tind = findfirst(x->x==typ,tps)
        if tind == nothing
            push!(rxnlists,[])
            push!(tps,typ)
            push!(rxnlists[end],rxn)
        else
            push!(rxnlists[tind],rxn)
        end
    end
    return (tps,rxnlists)
end
export splitreactionsbykinetics

"""
create vectorized kinetics calculators for the reactions
"""
function getveckinetics(rxns)
    tps,rxnlists = splitreactionsbykinetics(rxns)
    posinds = Array{Int64,1}()
    fs = []
    otherrxns = Array{ElementaryReaction,1}()
    vecinds = Array{Int64,1}()
    otherrxninds = Array{Int64,1}()
    for (i,tp) in enumerate(tps)
        if typeof(tp)<:Tuple
            typ = split(tp[1],".")[end] #this split needs done because in pyrms the type names come as ReactionMechanismSimulator.Arrhenius instead of Arrhenius in RMS proper
        else
            typ = split(tp,".")[end]
        end
        rinds = [findfirst(isequal(rxn),rxns) for rxn in rxnlists[i]]
        fcn = Symbol(typ * "vec")
        if !(fcn in allowedfcnlist) || occursin("Troe",typ) #no vectorized kinetics or for the moment stuff with efficiencies
            append!(otherrxns,rxnlists[i])
            append!(otherrxninds,rinds)
        else
            append!(vecinds,rinds)
            fcn = eval(fcn)
            x = fcn([x.kinetics for x in rxnlists[i]])
            push!(fs,x)
            if posinds == Array{Int64,1}()
                push!(posinds,length(rinds))
            else
                push!(posinds,length(rinds)+posinds[end])
            end
        end
    end
    vectuple = tuple(fs...)
    return (vectuple,vecinds,otherrxns,otherrxninds,posinds)
end
export getveckinetics

"""
create vectorized thermo calculator for species
"""
function getvecthermo(spcs)
    thermo = [sp.thermo for sp in spcs]
    if isa(thermo[1],NASA)
        typeassert.(thermo,NASA)
        return NASAvec([sp.thermo for sp in spcs])
    elseif isa(thermo[1],ConstantG)
        typeassert.(thermo,ConstantG)
        return ConstantGvec([th.G for th in thermo],thermo[1].T)
    else
        t = typeof(thermo[1])
        error("Thermo type $t unsupported!")
    end
end
export getvecthermo

function getC0(ph::X,T) where {X<:Union{IdealDiluteSolution,IdealGas}}
    return 1.0e5/(R*T)
end

function getC0(ph::X,T) where {X<:IdealSurface}
    return ph.sitedensity
end

export getC0

length(p::T) where {T<:AbstractPhase} = 1
export length

iterate(p::T) where {T<:AbstractPhase} = p
export iterate

Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
export broadcastable

function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
    arr = zeros(Int64,(6,length(rxns)))
    names = [spc.name for spc in spcs]
    for (i,rxn) in enumerate(rxns)
        inds = [findfirst(isequal(spc),spcs) for spc in rxn.reactants]
        for (j,spc) in enumerate(rxn.reactants)
            ind = findfirst(isequal(spc),spcs)
            arr[j,i] = ind
            rxn.reactantinds[j] = ind
        end
        for (j,spc) in enumerate(rxn.products)
            ind = findfirst(isequal(spc),spcs)
            arr[j+3,i] = ind
            rxn.productinds[j] = ind
        end
        if hasproperty(rxn.kinetics,:efficiencies) && length(rxn.kinetics.nameefficiencies) > 0
            while length(rxn.kinetics.efficiencies) > 0
                pop!(rxn.kinetics.efficiencies)
            end
            for (key,val) in rxn.kinetics.nameefficiencies
                ind = findfirst(isequal(key),names)
                rxn.kinetics.efficiencies[ind] = val
            end
        end
    end
    return arr
end
export getreactionindices


text:
module Cars

using SearchLight
import Base: @kwdef

export Car

@kwdef mutable struct Car <: AbstractModel
  id::DbId        = DbId()
  make::String    = ""
  max_speed::Int  = 220
end

end


text:
struct Nuclear <: Tabular end
url(::Nuclear) = "https://github.com/JuliaStats/RDatasets.jl/raw/master/data/boot/nuclear.csv.gz"
checksum(::Nuclear) = "0283f0f95a2df8e66ea40035483b6dcecb11fd856823f190ea943d774d0b9a44"
prep(::Nuclear) = path -> preprocess(extract(path), Nuclear())
target(::Nuclear) = 1
headers(::Nuclear) = 1
size(::Nuclear) = (32, 0, 0)
problem(::Nuclear) = Regression
message(::Nuclear) = "Nuclear Power Station Construction Data"

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
# ---
# title: 264. Ugly Number II
# id: problem264
# author: Indigo
# date: 2021-06-05
# difficulty: Medium
# categories: Math, Dynamic Programming, Heap
# link: <https://leetcode.com/problems/ugly-number-ii/description/>
# hidden: true
# ---
#
# Write a program to find the `n`-th ugly number.
#
# Ugly numbers are **positive numbers** whose prime factors only include `2, 3,
# 5`.
#
# **Example:**
#
#
#
#     Input: n = 10
#     Output: 12
#     Explanation:1, 2, 3, 4, 5, 6, 8, 9, 10, 12 is the sequence of the first 10 ugly numbers.
#
# **Note:**
#
#   1. `1` is typically treated as an ugly number.
#   2. `n` **does not exceed 1690**.
#
#
## @lc code=start
using LeetCode

function nth_ugly_number(n::Int)
    dp = Vector{Int}(undef, n)
    p2 = p3 = p5 = dp[1] = 1
    for i in 2:n
        nums = dp[p2] * 2, dp[p3] * 3, dp[p5] * 5
        dp[i] = minimum(nums)
        p2 += dp[i] == nums[1]
        p3 += dp[i] == nums[2]
        p5 += dp[i] == nums[3]
    end
    return dp[n]
end
## @lc code=end