text:
#=
using PyCall

@pyimport pickle

f = open(raw"D:\projets\PlanHeat\tmp\instances\not_dag_2\flows.pickle", "r")
flows=pickle.load(f)
close(f)

include("local_optimizer.jl")

flows = Dict{Tuple, Float32}((0, 1, 0) => 1,
             (1, 2, 0) => 1,
             (2, 3, 0) => 1,
             (1, 4, 0) => 1,
             (5, 6, 0) => 1,
             (6, 4, 0) => 1,
             (6, 7, 0) => 1,
             (7, 8, 0) => 1,
             (8, 10, 0) => 1,
             (10, 9, 0) => 1,
             (9, 10, 0) => 1,
             (9, 10, 1) => 1,
             (10, 3, 0) => 1)

flow_dag = get_flow_dag(flows)

edges = Set{Tuple}()
for (u, ns) in flow_dag
    for (v, ks) in ns
        for k in ks
            e = (u, v, k)
            push!(edges, e)
        end
    end
end
for e in edges
    re = (e[2], e[1], e[3])
    if re in edges
        println(e)
    end
end
start_node= "27307035"

k = find_cycle(flow_dag)

i_flows = copy(flows)

clean_flow_two_cycles(flows)

break_cycles(flows)

flow_dag = get_flow_dag(flows)

get_nodes(flow_dag)
descendants(flow_dag, "S_0")
descendants(flow_dag, start_node)

costs = Dict(e => 1.0 for e in keys(flows))
G = get_adj_dict(keys(flows))
dijkstra_multisource(G, Set(["S_0"]), start_node, costs)

optimization_graph = get_adj_dict(keys(flows))
=#

#=   ======================================================================== =#
include("FCFP.jl")
import FCFP.FCFPInstance
include("DSSP.jl")
using DSSP
using PyCall
using JSON
using DataStructures

@pyimport pickle

p = raw"D:\projets\PlanHeat\tmp\master_plugin\debug_local_optimizer2/"
f = open(p * "flows.pickle", "r")
flows=pickle.load(f)
close(f)

f = open(p * "costs.pickle", "r")
costs=pickle.load(f)
close(f)

f = open(p * "heat_demand.json", "r")
head_demand=JSON.parse(read(f, String))
close(f)

f = open(p * "production.json", "r")
productions = JSON.parse(read(f, String))
close(f)

include("local_optimizer.jl")
import LocalOptimizer:optimize_locally, MultiDiGraph, get_nodes, get_flow_dag, has_no_parent, add_edge!, get_edges
import LocalOptimizer:descendants, get_sub_graph, break_cycles, clean_flow_two_cycles, validate_move, LocalOptimizerInstance, get_flow_arteries, remove_node!, cascade_union, add_zero_weight_edges_from_nodes, dijkstra_multisource, get_edge_like_path, reroute_artery, lowest_common_ancestor
import LocalOptimizer:check_solution, find_cycle, VertexType, KeyType, EdgeType, MultiEdgeType, AdjencyView
import LocalOptimizer:FlowDictType, replace_arteries, WeightType, out_edges, in_edges, has_one_parent

optimization_graph = MultiDiGraph(keys(flows))
V = get_nodes(optimization_graph)
old_buildings = Set()
capacities = Dict()
fcfp = FCFPInstance(0.0, V, costs, head_demand, productions, capacities, old_buildings)
logger = println

using StatProfilerHTML
Profile.clear()
@profile new_flows = optimize_locally(flows, fcfp, logger)
statprofilehtml()

flows2 = FlowDictType(flows)
loi = LocalOptimizerInstance(flows, fcfp, logger)
initial_cost = sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)
check_solution(loi, loi.flows)

start = time()

loi = LocalOptimizerInstance(flows, fcfp, logger)
old_flows = copy(loi.flows)
initial_cost = sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)
clean_flow_two_cycles(loi.working_flows)
break_cycles(loi, loi.working_flows)
validate_move(loi)
while replace_arteries(loi)
    a = check_solution(loi, loi.flows)
    if a
        old_flows = copy(loi.flows)
    else
        println("now")
        break
    end
    gain = initial_cost - sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)
    loi.logger("Cost gain: $(gain)")
end

check_solution(loi, old_flows)
check_solution(loi, loi.flows)

f = open(p * "previous_correct_flows.pickle", "w")
pickle.dump(old_flows, f)
close(f)

################################################################################
################### Understanding while check is false #########################
################################################################################

flows = loi.flows
flow_dag = get_flow_dag(loi.optimization_graph, loi.flows)
c = find_cycle(flow_dag)
tol = 1e-5
nodes = get_nodes(flow_dag)
for n in nodes
    out_flow = sum(Array{WeightType}([flows[e] for e in out_edges(flow_dag, n)]))
    in_flow = sum(Array{WeightType}([flows[e] for e in in_edges(flow_dag, n)]))
    if n in loi.fcfp.building_nodes
        continue
    elseif n in loi.fcfp.production_nodes
        if out_flow > (loi.fcfp.production[n] + tol)
            println("here")
            break
        end
    else
        if abs(out_flow - in_flow) > tol
            println("or here")
            println(n)
            break
        end
    end
end

n= "intersect_S_2"
abs(sum([loi.flows[e] for e in out_edges(flow_dag, n)]) - sum([loi.flows[e] for e in in_edges(flow_dag, n)]))
v = "6341231519"
minimum(values(loi.fcfp.heat_demand))

ns = [u for u in get_nodes(flow_dag)]
v in ns
loi.flows[(v, n , 0)]

out_flow = sum(Array{WeightType}([flows[e] for e in out_edges(flow_dag, n)]))
in_flow = sum(Array{WeightType}([flows[e] for e in in_edges(flow_dag, n)]))

previous_flow = loi.flows
f = open(p * "flows_sources_move_pr.pickle", "w")
pickle.dump(previous_flow, f)
close(f)

################################################################################
################### Finding the artery that is involved ########################
################################################################################

loi.flows = copy(old_flows)
loi.working_flows = copy(old_flows)

loi.roots
[loi.fcfp.production[r] for r in loi.roots]

flow_dag = get_flow_dag(loi.optimization_graph, loi.flows)
arteries = get_flow_arteries(loi, flow_dag)
roots_descendants = Array{Set{VertexType}, 1}([descendants(flow_dag, r, false) for r in loi.roots])
roots_marginal_capacity = [loi.fcfp.production[r] - sum(get(loi.fcfp.heat_demand, n, 0.0)
                                                        for n in roots_descendants[r_index])
                                                                for (r_index, r) in enumerate(loi.roots)]
a = []
edge_like_path = []
while !isempty(arteries)
    a = dequeue!(arteries)
    # determine sub tree
    end_node = last(a)[2]
    sub_tree_nodes = descendants(flow_dag, end_node, false)
    if isempty(sub_tree_nodes)
       print("empty here")
    end
    sub_tree_demand = sum(get(loi.fcfp.heat_demand, n, 0.0) for n in sub_tree_nodes)
    # determine parent sub tree
    start_node = first(a)[1]
    current_roots_indexes = Set{Int64}()
    other_roots_indexes = Set{Int64}()
    for (d_index, d) in enumerate(roots_descendants)
       if (start_node in d)
           push!(current_roots_indexes, d_index)
       elseif (roots_marginal_capacity[d_index] >= sub_tree_demand)
           push!(other_roots_indexes, d_index)
       end
    end
    concerned_sub_trees_indexes = union(current_roots_indexes, other_roots_indexes)
    parent_dag_trees_nodes = cascade_union([roots_descendants[d_index] for d_index in concerned_sub_trees_indexes])
    parent_dag_trees_nodes = setdiff(parent_dag_trees_nodes, union(sub_tree_nodes,
                                                                  Set(e[2] for e in a)))
    if isempty(parent_dag_trees_nodes)
       print("empty here2")
    end
    # Add 0-weight edges in the sub-graph
    modified_costs = add_zero_weight_edges_from_nodes(sub_tree_nodes, loi.costs)
    path, cost = dijkstra_multisource(loi.optimization_graph,
                                     parent_dag_trees_nodes,
                                     end_node,
                                     modified_costs)
    edge_like_path = get_edge_like_path(path, loi.optimization_graph, loi.fcfp.fixed_costs)
    if edge_like_path != a
       print("hey")
       break
    end
end
edge_like_path
a
[loi.flows[e] for e in a]
[has_one_parent(flow_dag, e[2]) for e in a]
flow_dag = get_flow_dag(loi.optimization_graph, loi.working_flows)

potential_cycle_creation = false
if any((n != end_node) && (n in sub_tree_nodes) for n in path)
    potential_cycle_creation = true
end
new_start = first(path)
any(new_start in roots_descendants[d_index] for d_index in current_roots_indexes)
if any(new_start in roots_descendants[d_index] for d_index in current_roots_indexes)
    loi.logger("Moving flow inside current dag")
    reroute_artery(loi, a, edge_like_path, flow_dag, potential_cycle_creation)
else
    loi.logger("Moving flow to another dag")
    new_roots_candidates = Set{VertexType}([loi.roots[r] for r in other_roots_indexes if new_start in roots_descendants[r]])
    current_roots = Set{VertexType}([loi.roots[r] for r in current_roots_indexes])
    move_artery_to_other_dag(loi, a, edge_like_path, flow_dag, current_roots, new_roots_candidates, potential_cycle_creation)
end

###### debug : move_artery_to_other_dag ################################
new_a = edge_like_path
artery_flow = loi.working_flows[first(a)]
for e in a
    loi.working_flows[e] = zero(WeightType)
end
# Find path from old root to old start and remove flow
old_start = first(a)[1]
residual_artery_flow = artery_flow
r = first(current_roots)
for r in current_roots  # Removing artery flow taking into account multiple roots
    root_flow = sum(loi.working_flows[e] for e in out_edges(flow_dag, r))
    flow_to_remove = min(root_flow, residual_artery_flow)
    (back_path, _) = dijkstra_multisource(flow_dag, Set{VertexType}([r]), old_start, loi.costs)
    edge_back_path = get_edge_like_path(back_path, flow_dag, loi.fcfp.fixed_costs)
    for e in edge_back_path
        loi.working_flows[e] -= flow_to_remove
    end
    residual_artery_flow -= flow_to_remove
    if residual_artery_flow <= 0.0
        break
    end
end
# Find path from new root to new start and add flow
new_start = first(new_a)[1]
(down_path, _) = dijkstra_multisource(flow_dag, new_roots_candidates, new_start, loi.costs)
edge_down_path = get_edge_like_path(down_path, flow_dag, loi.fcfp.fixed_costs)
e = edge_down_path[2]
for e in edge_down_path
    loi.working_flows[e] += artery_flow
end
n = "intersect_S_2"

e2 = edge_down_path[2]
v1 = Float64(loi.flows[e2]) + Float64(artery_flow)
e1 = edge_down_path[1]
e3 = ("intersect_S_2", "intersect_5321182.0", 0)
v2 = Float64(loi.flows[e1]) + Float64(artery_flow) - Float64(loi.flows[e3])
v1 - v2

sum([loi.working_flows[e] for e in out_edges(flow_dag, n)]) - (sum([loi.working_flows[e] for e in in_edges(flow_dag, n)]))
[e for e in in_edges(flow_dag, n)]
# Set new artery flow
for e in new_a
    loi.working_flows[e] += artery_flow  # Adding flow because their may be already a flow on it
end
# Clean flows
clean_flow_two_cycles(loi.working_flows)
if potential_cycle_creation
    break_cycles(loi, loi.working_flows, Set{VertexType}([new_start]))
end

validate_move(loi)

new_a = edge_like_path
artery_flow = loi.working_flows[first(a)]
for e in a
    loi.working_flows[e] = zero(WeightType)
end
new_start = first(new_a)[1]
start_node = first(a)[1]
if new_start != start_node
    lca = lowest_common_ancestor(flow_dag, start_node, new_start)
    if lca === nothing
        throw(ErrorException("Unable to find common ancestor"))
    end
    # 1: Browse back the flow to not show flow in this direction
    if lca != start_node
        back_path, _ = dijkstra_multisource(flow_dag, Set{VertexType}([lca]), start_node, loi.costs)
        edge_back_path = get_edge_like_path(back_path, flow_dag, loi.fcfp.fixed_costs)
        for e in edge_back_path
            if artery_flow <= loi.working_flows[e]
                loi.working_flows[e] -= artery_flow
            else
                potential_cycle_creation = true
                loi.working_flows[(e[2], e[1], e[3])] = artery_flow - loi.working_flows[e]
                loi.working_flows[e] = zero(WeightType)
            end
        end
    end
    # 2: Browse down the flow to the start of the new artery
    if lca != new_start
        (down_path, _) = dijkstra_multisource(flow_dag, Set{VertexType}([lca]), new_start, loi.costs)
        edge_down_path = get_edge_like_path(down_path, flow_dag, loi.fcfp.fixed_costs)
        for e in edge_down_path
            loi.working_flows[e] += artery_flow
        end
    end
end
for e in new_a
    loi.working_flows[e] += artery_flow  # Adding flow because their may be already a flow on it
end
clean_flow_two_cycles(loi.working_flows)
if potential_cycle_creation
    break_cycles(loi, loi.working_flows, Set{VertexType}([new_start]))
end

gain = initial_cost - sum(loi.fcfp.fixed_costs[e] for (e, f) in loi.flows if f > 0)

optimization_time = time() - start
loi.logger("Arteries flows rerouting postprocess: gain=$(gain),  time: $(optimization_time)")
if check_solution(loi, loi.flows)
    return loi.flows
else
    loi.logger("Local optimization failed")
    return loi.init_flows
end

replace_arteries(loi)

import LocalOptimizer:replace_arteries, check_solution
check_solution(loi, loi.working_flows)

flow_dag = get_flow_dag(loi.optimization_graph, loi.flows)
arteries = get_flow_arteries(loi, flow_dag)
roots_descendants = Array{Set{VertexType}, 1}([descendants(flow_dag, r, false) for r in loi.roots])
roots_marginal_capacity = [loi.fcfp.production[r] - sum(get(loi.fcfp.heat_demand, n, 0.0)
                                                        for n in roots_descendants[r_index])
                           for (r_index, r) in enumerate(loi.roots)]

new_start = nothing
path = nothing
edge_like_path = nothing
a = nothing
potential_cycle_creation = nothing
using DataStructures

while !isempty(arteries)
    a = dequeue!(arteries)
    # determine sub tree
    end_node = last(a)[2]
    sub_tree_nodes = descendants(flow_dag, end_node, false)
    sub_tree_demand = sum(get(loi.fcfp.heat_demand, n, 0.0) for n in sub_tree_nodes)
    # determine parent sub tree
    start_node = first(a)[1]
    current_roots_indexes = Set{Int64}()
    other_roots_indexes = Set{Int64}()
    for (d_index, d) in enumerate(roots_descendants)
        if (start_node in d)
            push!(current_roots_indexes, d_index)
        elseif (roots_marginal_capacity[d_index] >= sub_tree_demand)
            push!(other_roots_indexes, d_index)
        end
    end
    concerned_sub_trees_indexes = union(current_roots_indexes, other_roots_indexes)
    parent_dag_trees_nodes = cascade_union([roots_descendants[d_index] for d_index in concerned_sub_trees_indexes])
    parent_dag_trees_nodes = setdiff(parent_dag_trees_nodes, union(sub_tree_nodes,
                                                                   Set(e[2] for e in a)))
    # Add 0-weight edges in the sub-graph
    modified_costs = add_zero_weight_edges_from_nodes(sub_tree_nodes, loi.costs)
    path, cost = dijkstra_multisource(loi.optimization_graph,
                                      parent_dag_trees_nodes,
                                      end_node,
                                      modified_costs)
    edge_like_path = get_edge_like_path(path, loi.optimization_graph, loi.fcfp.fixed_costs)
    if edge_like_path != a
        potential_cycle_creation = false
        if any((n != end_node) && (n in sub_tree_nodes) for n in path)
            potential_cycle_creation = true
        end
        new_start = first(path)
        if any(new_start in roots_descendants[d_index] for d_index in current_roots_indexes)
            loi.logger("Moving flow inside current dag")
            reroute_artery(loi, a, edge_like_path, flow_dag, potential_cycle_creation)
            break
        else
            loi.logger("Moving flow to another dag")
            println("NOWWW")
            break
        end
    end
end

edge_like_path

a

check_solution(loi, loi.working_flows)
validate_move(loi)

import LocalOptimizer: move_artery_to_other_dag, WeightType, VertexType, MultiEdgeType
other_roots_indexes = Set{Int64}(1)
new_roots_candidates = Set{VertexType}([loi.roots[r] for r in other_roots_indexes if new_start in roots_descendants[r]])
loi.working_flows = copy(loi.flows)

new_a = edge_like_path
artery_flow = loi.working_flows[first(a)]
for e in a
    loi.working_flows[e] = zero(WeightType)
end
new_start = first(new_a)[1]
# Find path from root to new start
(down_path, _) = dijkstra_multisource(flow_dag, new_roots_candidates, new_start, loi.costs)
edge_down_path = get_edge_like_path(down_path, flow_dag, loi.fcfp.fixed_costs)
for e in edge_down_path
    loi.working_flows[e] += artery_flow
end
# Set new artery flow
for e in new_a
    loi.working_flows[e] += artery_flow  # Adding flow because their may be already a flow on it
end
# Clean flows
clean_flow_two_cycles(loi.working_flows)
if potential_cycle_creation
    break_cycles(loi, loi.working_flows)
end

check_solution(loi, loi.working_flows)

flows = loi.working_flows
# Check that flow graph is a DAG
flow_dag = get_flow_dag(loi.optimization_graph, flows)
c = find_cycle(flow_dag)
if c !== nothing
    return false
end
# Check flow conservation
tol = 1e-5
nodes = get_nodes(flow_dag)
for n in nodes
    out_flow = sum(Array{WeightType}([flows[e] for e in out_edges(flow_dag, n)]))
    in_flow = sum(Array{WeightType}([flows[e] for e in in_edges(flow_dag, n)]))
    if n in loi.fcfp.building_nodes
        continue
    elseif n in loi.fcfp.production_nodes
        if out_flow > (loi.fcfp.production[n] + tol)
            println("$(n) : $(out_flow) > $((loi.fcfp.production[n] + tol))")
        end
    else
        if abs(out_flow - in_flow) > tol
            println("$(n) : $(out_flow) != $(in_flow)")
        end
    end
end
return true

import LocalOptimizer: out_edges, in_edges

previous_flow = loi.flows
f = open(raw"D:\projets\PlanHeat\tmp\instances\not_dag_3\flows_sources_move_pr.pickle", "w")
pickle.dump(previous_flow, f)
close(f)

f = open(raw"D:\projets\PlanHeat\tmp\instances\not_dag_3\flows_sources_move.pickle", "w")
pickle.dump(loi.working_flows, f)
close(f)
=#

[has_one_parent(flow_dag, e[2]) for e in current_artery]

a = dequeue!(arteries)

arteries = PriorityQueue{Array{MultiEdgeType, 1}, WeightType}()
explored_nodes = Set{VertexType}()
r = loi.roots[3]
nodes_stack = Array{VertexType, 1}([r])
while !isempty(nodes_stack)
    n = pop!(nodes_stack)
    if n in explored_nodes
        continue
    end
    push!(explored_nodes, n)
    out_flow_edges = collect(out_edges(flow_dag, n))
    println(n)
    if (length(out_flow_edges) != 1) || (! has_one_parent(flow_dag, n))
        for e in out_flow_edges
            push!(nodes_stack, e[2])
        end
        continue
    end
    println("here")
    println(n)
    break
end

n = "intersect_S_0_2"
out_flow_edges = collect(out_edges(flow_dag, n))
# A new artery is found, explore it in depth
current_artery = [first(in_edges(flow_dag, n)), out_flow_edges[1]]
sub_nodes_stack = [out_flow_edges[1][2]]
while !isempty(sub_nodes_stack)
    sub_node = pop!(sub_nodes_stack)
    push!(explored_nodes, sub_node)
    out_flow_edges = collect(out_edges(flow_dag, sub_node))
    if (length(out_flow_edges) != 1) || (! has_one_parent(flow_dag, n))
        # End of the artery
        artery_cost = - sum(loi.fcfp.fixed_costs[e] for e in current_artery)
        if current_artery in loi.treated_arteries
            artery_cost += 1e5
        end
        arteries[current_artery] = artery_cost
        for e in out_flow_edges
            push!(nodes_stack, e[2])
        end
        break
    else
        push!(current_artery, out_flow_edges[1])
        push!(sub_nodes_stack, out_flow_edges[1][2])
    end
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text:
module ANIMAL_ROUTER

using HTTP
# modified Animal struct to associate with specific user
mutable struct Animal
    id::Int
    userId::Base.UUID
    type::String
    name::String
end

# use a plain `Dict` as a "data store"
const ANIMALS = Dict{Int, Animal}()
const NEXT_ID = Ref(0)
function getNextId()
    id = NEXT_ID[]
    NEXT_ID[] += 1
    return id
end

# **simplified** "service" functions
function createAnimal(req::HTTP.Request, animal)
    animal.id = getNextId()
    ANIMALS[animal.id] = animal
    return animal
end

function getAnimal(req::HTTP.Request)
    animalId = HTTP.URIs.splitpath(req.target)[5] # /api/zoo/v1/animals/10, get 10
    return ANIMALS[animalId]
end

function updateAnimal(req::HTTP.Request, animal)
    ANIMALS[animal.id] = animal
    return animal
end

function deleteAnimal(req::HTTP.Request)
    animalId = HTTP.URIs.splitpath(req.target)[5] # /api/zoo/v1/animals/10, get 10
    delete!(ANIMALS, animal.id)
    return ""
end

end  # module ANIMAL_ROUTER

ext: jl


text:
module BuildTreeCoded

    using PhylogeneticTrees, JuMP, Gurobi
    @static if VERSION < v"0.7.0-DEV.2005"
        using Base.Test
    else
        using Test
    end

    println("========================================")
    println("Testing tree-building with default model")
    println("========================================")

    @testset "Solving a tree problem" begin

    	pd = PhylogeneticTrees.PopulationData(
    		"testdata/f3.Europe6.csv",
    		"testdata/f3.Europe6-covariance.csv"
    		)

        bt = PhylogeneticTrees.BinaryTree(3)
        tp = PhylogeneticTrees.TreeProblem(pd, bt,
            solver = GurobiSolver(OutputFlag = 0),
            binaryencoding = true)
        @time solve(tp.model)
        # 5.566190 seconds (10.90 M allocations: 532.489 MiB, 5.85% gc time)

        leaves = PhylogeneticTrees.getnodes(bt, bt.depth)

        # test integrality
        for a in 1:pd.npop, u in leaves
            @test isapprox(getvalue(tp.assign[a,u,1]), 1) || isapprox(getvalue(tp.assign[a,u,1]), 0)
        end
        for (a,b,(u,v)) in keys(tp.countedge)
            val = round(getvalue(tp.countedge[a,b,(u,v),1,1]),digits=3)
            @test isapprox(val,1) ||
                  isapprox(val,0)
        end

        xval = [leaves[findfirst(round.(getvalue(tp.assign[a,:,1])) .> 0)] for a in 1:pd.npop]
        @test isapprox(round(getvalue(tp.f3formula[1,1])), 478)
        @test isapprox(round(getvalue(tp.f3formula[1,2])), 33)
        @test isapprox(round(getvalue(tp.f3formula[2,2])), 247)
        for a in 1:pd.npop, b in a:pd.npop, u in leaves, v in leaves
            if round(getvalue(tp.assign[a,u,1]) + getvalue(tp.assign[b,v,1])) == 2
                for (u1,v1) in bt.edges
                    val = round(getvalue(tp.countedge[a,b,(u1,v1),1,1]), digits=3)
                    if in((u1,v1), intersect(bt.pathedges[u,tp.outgroupnode], bt.pathedges[v,tp.outgroupnode]))
                        @test isapprox(val,1)
                    else
                        @test isapprox(val,0)
                    end
                end
            end
        end

        @test isapprox(getobjectivevalue(tp.model), 267.7057665969791)

    end

end


text:
using Constructions
using Test

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


text:
module Server
using Rock.RPC
# TODO: Clean this up
import Rock.RPC: BaseMessage, Command , Response , NonCriticalCommand, ExternalNonCriticalCommand, InternalNonCriticalCommand, CriticalCommand , InternalCriticalCommand, ExternalCriticalCommand, ConsensusCommand , ConsensusResponse , Fate, Decided, Pending, Forgotten, Slot, PrepareArgs , PrepareReply , AcceptArgs , AcceptReply , DecideArgs , DecideReply , Ping , PingReply , DeadLeader , AddPeer , DeletePeer , SwapPeer , NoOp , ExtraInfo, WithExtraInfo, SlotNum, Hostname, Port

using LockUtils
using SQLite
using Logging
@Logging.configure(level=INFO)
############################
##    Base/Abstract types ##
############################
typealias ServerId Int64
typealias  LogIndex Int64
typealias Term Int64

# Our abstract notion of t
####################################
##   t t Components    ##
##   Not "integrated" to abstract ##
##   FS out                       ##
####################################
type Context
    port::Port
    peers::Array{Tuple{ASCIIString, Int64}, 1}
    user_defined_function::Function
    db_loc::AbstractString
    function Context(port::Port, peers::Array{Tuple{ASCIIString, Int64},1}, user_defined_function::Function, db_loc::ASCIIString)
        d = SQLite.DB(db_loc)
        SQLite.transaction(d) do
            try
                SQLite.query(d, "create table parameters (K TEXT Unique PRIMARY KEY, V BLOB)")
            catch SQLite.SQLiteException
            end
            try
                SQLite.query(d, "create table slots (slot INTEGER PRIMARY KEY AUTOINCREMENT, n_p INT64, n_a INT64, fate INT, v_a BLOB)")
            catch SQLite.SQLiteException
            end
            try
                SQLite.query(d, "create table peers (host TEXT, port Int, unique (host, port))")
                for (h, p)= peers
                    SQLite.query(d, "insert or ignore into peers values (?, ?)", [h, p])
                end
            catch SQLite.SQLiteException
            end

        end
        new(port, peers, user_defined_function, db_loc)

    end
end
function db(context::Context)
    SQLite.DB(context.db_loc)
end
type t
    self::Int64
    completed::Array{Int64, 1}
    forgot::Int64
    context::Context
    lock::Lock
    dblock::ReentrantLock
    rpc_task::Task
    log_task::Task
    function t(context::Context, self::Int64)

        d  = db(context)
        v = SQLite.query(d, "select V from parameters where (K = 'forgot')").data[1]
        forgot = 0
        if length(v) > 0 && !isnull(v[1])
            forgot = v[1].value
        end
        completed = fill(forgot,length(context.peers))
        v = SQLite.query(d, "select V from parameters where (K = 'completed')").data[1]
        if length(v) > 0 && !isnull(v[1])
            completed[self] = v[1].value
        end

        new(self, completed, forgot, context, Lock(), ReentrantLock())
    end
end

function start(px::t)
    @info "Starting Node on $(port(px))"
    px.log_task = @async LogCrawler(px)
    px.rpc_task = @async RPCServer(px)
    wait(px.rpc_task)
    # @async Heartbeat(px)
end

function port(px::t)#::Int64
    px.context.port
end
function peers(px::t)#::Array{Tuple{Hostname, Port},1}
    px.context.peers
    # d = db(px)
    # data = SQLite.query(d, "select host, port from peers").data
    # collect(zip(map(v -> v.value, data[1]), map(v -> v.value, data[2])))
end
function addpeer(px::t, host, port)#::Array{Tuple{Hostname, Port},1}
    d = db(px)
    data = SQLite.query(d,  "insert or ignore into peers values (?, ?)", [host, port])
end
function deletepeer(px::t, host, port)#::Array{Tuple{Hostname, Port},1}
    d = db(px)
    data = SQLite.query(d, "delete from peers where (host = ?, port = ?)", [host, port])
end
function user_defined_function(px::t, d::SQLite.DB, v::ExternalCriticalCommand)#:: Function
    px.context.user_defined_function(px, d, v)
end
function user_defined_function(px::t, d::SQLite.DB, v::ExternalNonCriticalCommand)#:: Function
    px.context.user_defined_function(px, d, v)
end
function user_defined_function(px::t)
    px.context.user_defined_function
end
function self(px::t)#::Int64
    px.self
end
function completed(px::t)#::Array{Int64, 1}
    px.completed
end
function forgot(px::t)##::Int64
    px.forgot
end
function context(px::t)#::Context
    px.context
end
function lock(px::t)#::Lock
    px.lock
end
function db(px::t)
    db(context(px))
end
function syncSlot(d::SQLite.DB, slot::Slot)
    SQLite.query(d, "insert or replace into slots values (?, ?, ?, ?, ?)",[slot.num, slot.n_p, slot.n_a, slot.fate, slot.v_a  ])
end

function updateCompleted(px::t, completed)
    px.completed = max(px.completed, completed) # Note this is vectorized
    px.forgot = min(px.completed...)
end
function syncCompleted(px::t, d::SQLite.DB)
    SQLite.query(d, "insert or replace into parameters values ('completed', ?)",[completed(px)[self(px)]])
end
function syncForgot(px::t, d::SQLite.DB)
    # Update the forgot index first because we want to prevent the
    # case where we delete a file, then crash, and think our
    # forgot index is too high
    blockingTransaction(px,d) do
        SQLite.query(d, "insert or replace into parameters values ('forgot', ?)",[forgot(px)])
        SQLite.query(d, "delete from slots where  (slot <= ?)",[forgot(px)])
        syncCompleted(px, d)
    end
end

############################################
#          Handlers                        #
############################################
function Handle(px::t, args::PrepareArgs)
    @debug "PREPARE"
    d = db(px)
    blockingTransaction(px, d) do
        s = unsafe_getSlot!(px,d,args.slot_num)
        if args.n > s.n_p
            persistently(d, s) do
                s.n_p = args.n
            end
            WithExtraInfo(PrepareReply(args.slot_num, true, s.n_p, s.n_a, s.v_a), ExtraInfo(px))
        else
            WithExtraInfo(PrepareReply(args.slot_num, false, s.n_p, s.n_a, s.v_a), ExtraInfo(px))
        end
    end
end
function Handle(px::t, args::AcceptArgs)
    d = db(px)
    blockingTransaction(px,d) do
        s = unsafe_getSlot!(px,d, args.slot_num)
        if args.n >= s.n_p
            persistently(d,s) do
                s.n_p = args.n
                s.n_a = args.n
                s.v_a = Nullable(args.v)
            end
            WithExtraInfo(AcceptReply(args.slot_num, true, s.n_p, args.v), ExtraInfo(px))
        else
            WithExtraInfo(AcceptReply(args.slot_num, false, 0, NoOp()), ExtraInfo(px))
        end
    end

end

function Handle(px::t, args::DecideArgs)
    d = db(px)
    blockingTransaction(px,d) do
        s = unsafe_getSlot!(px,d, args.slot_num)
        if s.fate != Decided
            persistently(d, s) do
                s.v_a = Nullable(args.v)
                s.fate = Decided
            end
            # @async signal(s.condition)
        end
        WithExtraInfo(DecideReply(args.slot_num, true), ExtraInfo(px))
    end
end
function Handle(px::t, args::Ping)
    @debug args.msg
    PingReply()
end
function Handle(px::t, args::PingReply)
    @warn "THIS SHOULD NOT HAPPEN EVER"
end
function Handle(px::t,args::NoOp)
end
#=
Piggybacking Commands
=#

function ExtraInfo(px::t)
    ExtraInfo(completed(px))
end
function safe_getExtraInfo(px::t)
    holding(lock(px), "getExtraInfo") do
        ExtraInfo(px)
    end
end
@inline function getExtraInfo(v::WithExtraInfo)
    v.info
end

function Handle(px::t, command::WithExtraInfo)
    try
        @async Forgettor(px, command.info.completed)
    catch err
        @error err
    end
    Handle(px, command.command)
end

function Handle(px::t, args::CriticalCommand)
    try
        d = db(px)
        while true
            t = holding(lock(px), "Handle Command") do
                @debug "GOT LOCK 1"

                blockingTransaction(px,d)  do
                @debug "GOT LOCK 2"
                    SQLite.query(d, "insert into slots values (null, ?, ?, ?, ?)", [0, 0,  Pending,Nullable{CriticalCommand}()])
                    # Only safe because on one connection
                    s = SQLite.query(d, "select last_insert_rowid()").data[1][1].value
                    @async Proposer(px, args, s, true)
                end
            end
                @debug "WAITING"
            if consume(t)

                @debug "FINISHED"
                return true
            end
            @debug "RETRYING"
        end
    catch err
        @debug err
        return false
    end

end
function Handle(px::t, args::ExternalNonCriticalCommand)
    user_defined_function(px,db(px), args)
end

##########################################
#         Utility                       ##
##########################################

function blockingTransaction(f::Function, px::t, d::SQLite.DB)
    # while true
    #     try
    #         @debug "Trying transaction"
    Base.lock(px.dblock)
    v = SQLite.transaction(d) do
        f()
    end
    Base.unlock(px.dblock)
    return v
    #     catch err
    #         if isa(err, SQLite.SQLiteException) && err.msg == "database is locked"
    #             @debug "Transaction Retrying"
    #             wait(Timer(0.01))
    #             continue
    #         else
    #             rethrow()
    #         end
    #     end
    # end
end

@inline function notNull(f::Function, m_v, def)
    if !isnull(m_v)
        f(m_v.value)
    else
        def
    end

end
@inline function notNull(f::Function, m_v)
    if !isnull(m_v)
        f(m_v.value)
    end
end

function message(px::t, m::BaseMessage, chan, i::Int64)
    host, port = peers(px)[i]
    conn = connect(host, port)
    serialize(conn, m)
    val = deserialize(conn)
    close(conn)
    v =  if isa(val,WithExtraInfo)
        @async Forgettor(px, val.info.completed)
        val.command
    else
        val
    end
    notNull(chan) do ch
        put!(ch, v)
    end
    v
end
function broadcast(px::t, m::BaseMessage, chan)
    @async for i = 1:length(peers(px))
        message(px, m, Nullable(chan), i)
    end
end
function broadcast(px::t, m::BaseMessage)
    @async for i = 1:length(peers(px))
        message(px, m, Nullable(),i)
    end
end

function unsafe_getSlot!(px::t,d::SQLite.DB, s::SlotNum)
    @debug "unsafe_getSlot! $(forgot(px)) $s"
    if s <= forgot(px)
        Slot(0,0,0,Nullable(), Forgotten)
    else
        obj = SQLite.query(d, "select slot, n_p,n_a, v_a, fate from slots where (slot=?) limit 1", [s]).data
        if length(obj[1]) > 0
            @debug "Got Slot $(Slot(obj[1][1].value, obj[2][1].value,obj[3][1].value, obj[4][1].value, obj[5][1].value))"

            Slot(obj[1][1].value, obj[2][1].value,obj[3][1].value, obj[4][1].value, obj[5][1].value)
        else
            v = Slot(s, 0,0,Nullable(), Pending)
            SQLite.query(d, "insert into slots values (?, ?, ?, ?, ?)", [s, 0, 0,  Pending, Nullable{CriticalCommand}()])
            v
        end
    end
end

function getSlot!(px::t, s::SlotNum)
    d = db(px)
    blockingTransaction(px,d) do
        unsafe_getSlot!(px,d, s)
    end
end
function safe_getSlot!(px::t, s::SlotNum)
    holding(lock(px), "safe_getSlot! $(s)") do
        getSlot!(px, s)
    end
end

function getSlot(px::t, s::SlotNum)
    d = db(px)
    if s <= forgot(px)
        Nullable(Slot(0,0,0,Nullable(), Forgotten))
    else
        obj = SQLite.query(d, "select slot, n_p,n_a, v_a, fate from slots where (slot=?) limit 1", [s]).data
        if length(obj[1]) > 0
            Nullable(Slot(obj[1][1].value, obj[2][1].value,obj[3][1].value, obj[4][1].value,(obj[5][1].value)))
        else
            Nullable()
        end
    end

end
function safe_getSlot(px::t, s::SlotNum)
    holding(lock(px), "safe_getSlot $(s)") do
        getSlot(px, s)
    end
end

function persistently(f::Function,  d::SQLite.DB, slot::Slot)
    f()
    syncSlot(d, slot)
end

######################
##     Tasklets     ##
######################

#=
#              Proposer functionality
=#
function AcceptPhase(px::t,slot_num::SlotNum, n::Int64, v::CriticalCommand, majority::Int64)
    accept = WithExtraInfo(AcceptArgs(slot_num, n, v),  safe_getExtraInfo(px))
    chan = Base.Channel{AcceptReply}(5)
    broadcast(px, accept, chan)
    positives = 0
    negatives = 0
    while true
        r = take!(chan)

        if r.Ok
            positives += 1
            if positives >= majority
                decide = WithExtraInfo(DecideArgs(slot_num, v), safe_getExtraInfo(px))
                broadcast(px, decide)
                return true
            end

        else
            negatives +=1
            if negatives >= majority
                return false
            end
        end

    end
end

function Proposer(px::t, v::CriticalCommand, slot_num::SlotNum, ret = false)
    epoch = -1
    majority = div(length(peers(px)) +2, 2)
    has_produced = false
    while true
        pending = holding(lock(px), "Proposer While Loop $slot_num") do
            @debug "GETTING FATE"
            getSlot!(px, slot_num).fate == Pending
        end
        if !pending
            break
        else
            epoch += 1
            n = length(peers(px))*epoch + self(px) + 1
            chan = Base.Channel{PrepareReply}(length(peers(px)))# We don't want it to block ever
            prep = WithExtraInfo(PrepareArgs(slot_num, n), safe_getExtraInfo(px))
            @info "Proposing on slot $slot_num on $(self(px)) with value $v"
            broadcast(px, prep, chan)
            positives = 0
            negatives = 0
            n_a = 0
            v_a = v
            loop = true
            while loop
                r = take!(chan)
                if r.Ok
                    positives += 1
                    if r.n_a > n_a
                        notNull(r.v_a) do v_new
                            if ret && v_a != v_new && !has_produced
                                has_produced = true
                                produce(false)
                            end
                            v_a = v_new
                        end
                    end
                    if positives >= majority
                        a = AcceptPhase(px, slot_num, n, v_a, majority)
                        if a
                            if !has_produced
                                produce(true)
                            end
                            loop = false
                        end
                    end
                else

                    negatives += 1
                end
                if negatives >= majority
                    loop = false
                end
            end

        end
    end

end

function RPCServer(px::t)
    server = listen(port(px))
    @info "RPC Server Starting on port $(port(px)), node $(self(px))"
    while true
        conn = accept(server)
        @info "Incoming Connection on node $(self(px))"
        @async begin
            while true
                try
                    args = deserialize(conn)::BaseMessage
                    @info "RPCServer $(self(px)) Args: $args"
                    try
                        r = Handle(px, args)
                        @info "RPCServer $(self(px)) Response: $r"
                        serialize(conn,r)
                        flush(conn)
                    catch err
                        @debug "RPCServer $(self(px)) Error Handling:\n    $args\n     $err"
                        close(conn)
                        break
                    end
                catch err
                    @debug "RPCServer $(self(px)) Connection Failed With $err"
                    close(conn)
                    break
                    # finally
                    #     close(conn)
                end
            end
        end
    end
end
function LogCrawler(px::t)
    @debug "Log Crawler Starting"
    while true
        changed = false
        while true
            next = safe_getSlot(px, completed(px)[self(px)]+1)
            if isnull(next)
                wait(Timer(0.01))
                break
            end
            next = next.value
            times = Int64(floor((rand()*10+5.0)/0.01)) # Wait at least 5 seconds worth, but randomize
            count = 0
            fired = false
            while true
                if count > times && !fired
                    #Timeout
                    fired = true
                    @warn "Timeout in Log Crawler, starting Proposer"
                    @async Proposer(px, NoOp(), next.num) # TODO: Can we guarantee that command has a No_Op more idomatically?
                    wait(Timer(0.01))
                end
                count += 1
                d = db(px)
                finished  = blockingTransaction(px,d) do
                    result = SQLite.query(d, "select v_a from slots where fate = ? and slot = ? limit 1", [Decided, next.num]).data
                    if length(result[1])> 0 && !isnull(result[1][1])
                        command = result[1][1].value.value # must not be null, otherwise consensus broken somewhere else
                        completed(px)[self(px)] += 1
                        syncCompleted(px,d)
                        @info "LogCrawler completed Command $command on $(self(px))"
                        if typeof(command) <: ExternalCriticalCommand
                            @debug "running user_defined_function function"
                            user_defined_function(px,d, command)
                        elseif typeof(command) <: InternalCriticalCommand
                        else
                        end
                        true
                    else
                        false
                    end
                end
                # @debug "Completed a log entry? $finished"
                if finished
                    break
                else
                    wait(Timer(0.01))
                end
            end
        end
        wait(Timer(1))
    end
end
function Forgettor(px::t, completed::Array{Int64,1})
    @debug "Forgetting up to $completed"
    holding(lock(px), "Forgettor") do
        updateCompleted(px, completed) # updates completed and forgotten fields and deletes elements from log
        d = db(px)
        syncForgot(px,d)

    end
end

function Heartbeat(px::t)
    while true
        @debug "Begin Heartbeat"
        arr = holding(lock(px)) do
            arr = fill(true, length(peers(px)))
            m = Ping("TEST")
            @async for i = 1:length(peers(px))
                host, port = peers(px)[i]
                try
                    conn = connect(host, port)
                    try
                        serialize(conn, m)
                        deserialize(conn)
                        # arr[i] = false
                    catch err
                        # Wait, then Retry!
                        wait(Timer(5))
                        conn2 = connect(host, port)
                        try
                            serialize(conn2, m)
                            flush(conn2)
                            deserialize(conn2)
                            # arr[i] = false
                        catch err
                            @debug err
                        finally
                            close(conn2)
                        end
                    finally
                        close(conn)
                    end
                catch err
                    @debug err
                end
            end
            arr
        end
        wait(Timer(10))
        for (i,v) = enumerate(arr)
            if v
                @warn "$i dead"
            end
        end
        @debug "End Heartbeat"

    end
end
end

max_stars_repo_name: JeremyRubin/Rock


text:
using DataFrames, CSV, CodecLz4

IN_CSV_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_CSV_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)
OUT_CSV_STREAM(p::String) = LZ4CompressorStream(open(p, "w"))
OUT_CSV_STREAM(F::Function, p::String) = open(F, LZ4CompressorStream, p, "w")

function mwe(N)
  local dfcsv::DataFrame
  #df = DataFrame(x1 = rand((string).(1:N),N), x2 = rand((string).(1:N),N))
  f = "test.csv"

  #sending the df to a CSV first breaks the code
  #open(LZ4CompressorStream, f, "w") do s
  #  CSV.write(s, df)
  #end
  #df |> CSV.write |> LZ4CompressorStream(open(f, "w"))

  #open(LZ4DecompressorStream, "$(pwd())\\data\\COMP-A.csv.lz4") do s
  #  dfcsv = CSV.read(s, threaded=true) |> DataFrame
  #end
  #df |> CSV.write(f)

  #open("$(pwd())\\data\\COMP-A.csv") do s
  #  dfcsv = CSV.read(s, threaded=false) |> DataFrame
  #end

  #dfcsv[1:75_000,8:9]  |> CSV.write(f)
  open(f) do s
    dfcsv = CSV.read(s, threaded=false) |> DataFrame
  end
  @info "unthreaded read seems to work fine"
  open(f) do s
    dfcsv = CSV.read(s, threaded=true) |> DataFrame
  end
  @info "Never gets here because the threaded read fails"

  #dfcsv = CSV.read("test.csv", copycols=false, threaded=true) |> DataFrame
  println("size dfcsv: ", size(dfcsv))

end

function mwe2()
  local dfcsv::DataFrame
  f = "test.csv"

  open(f) do s
    dfcsv = CSV.read(s, threaded=false) |> DataFrame
  end
  @info "unthreaded read seems to work fine"
  open(f) do s
    dfcsv = CSV.read(s, threaded=true) |> DataFrame
  end
  @info "Never gets here because the threaded read fails"

  println("size dfcsv: ", size(dfcsv))

end

mwe2()

@inline function fixcsv(df::DataFrame)::DataFrame
  #=for c::Int ∈ 1:size(df,2)
    df[!,c] = Array(df[!,c])
  end=#
  return DataFrame(view(df, :, :))
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using OceanWaveSpectralFitting
using Documenter

DocMeta.setdocmeta!(OceanWaveSpectralFitting, :DocTestSetup, :(using OceanWaveSpectralFitting); recursive=true)

makedocs(;
    modules=[OceanWaveSpectralFitting],
    authors="Jake Grainger <j.p.grainger2@outlook.com> and contributors",
    repo="https://github.com/JakeGrainger/OceanWaveSpectralFitting.jl/blob/{commit}{path}#{line}",
    sitename="OceanWaveSpectralFitting.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JakeGrainger.github.io/OceanWaveSpectralFitting.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Models" => [
            "univariate.md",
            "multivariate.md"
        ]
    ],
)

deploydocs(;
    repo="github.com/JakeGrainger/OceanWaveSpectralFitting.jl",
    devbranch="main",
)


text:
# Deprecated functionality

import Base: @deprecate_binding, @deprecate_moved


text:
src = "../src"
if !in(src, LOAD_PATH)
    push!(LOAD_PATH, src)
end
using Tel

#=
ev = readcsv("../data/event_type.csv", ASCIIString)
re = readcsv("../data/resource_type.csv", ASCIIString)
se = readcsv("../data/severity_type.csv", ASCIIString)
reformat(ev)
reformat(re)
reformat(se)
=#

log = readcsv("../data/log_feature.csv", ASCIIString)
Tel.reformat3(log)


text:
module ClangCompiler

const __DLEXT = Base.BinaryPlatforms.platform_dlext()
const __ARTIFACT_BINDIR = Sys.iswindows() ? "bin" : "lib"

include("jllshim.jl")
using .JLLShim

using Clang_jll
using libclangex_jll

# need manually add the output path below to `DYLD_LIBRARY_PATH`/`LD_LIBRARY_PATH`/`PATH`
const libclangex_include = joinpath(libclangex_jll.artifact_dir, "include") |> normpath

const julia_include_dir = joinpath(Sys.BINDIR, "..", "include", "julia") |> normpath

using LLVM
using LLVM.Interop: call_function
export call_function
import LLVM: dispose, name, value

import Base: dump, string

const CLANG_BIN = joinpath(Clang_jll.artifact_dir, "bin", "clang")
const CLANG_INC = joinpath(Clang_jll.artifact_dir, "lib", "clang", string(LLVM.version()), "include")

include("../lib/LibClangEx.jl")
using .LibClangEx

include("platform/JLLEnvs.jl")
using .JLLEnvs

function get_compiler_args(;is_cxx=true, version=v"7.1.0")
    args = JLLEnvs.get_default_args(is_cxx, version)
    push!(args, "-isystem"*CLANG_INC)
    is_cxx && push!(args, "-nostdinc++", "-nostdlib++")
    push!(args, "-nostdinc", "-nostdlib")
    pushfirst!(args, CLANG_BIN)  # Argv0
    return args
end
export get_compiler_args

# low-level
include("clang/utils.jl")
include("clang/core/core.jl")
include("clang/api/api.jl")
include("clang/ast.jl")
include("clang/basic.jl")
include("clang/codegen.jl")
include("clang/frontend.jl")
include("clang/lex.jl")
include("clang/parse.jl")
include("clang/qualtype.jl")
include("clang/sema.jl")
export AbstractClangType, AbstractQualType, AbstractBuiltinType

# high-level
include("types.jl")
export jlty_to_clty, clty_to_jlty, jlty_to_llvmty

include("parse.jl")
include("lookup.jl")
export DeclFinder, get_decl

include("template.jl")
export specialize

include("compiler.jl")
export AbstractCompiler, AbstractIRGenerator
export IRGenerator, IncrementalIRGenerator
export CXCompiler
export get_instance, get_context
export get_module, get_jit, get_dylib, get_codegen
export compile, dispose
export link_process_symbols
export get_modules, get_current_module
export incremental_parse

include("llvm.jl")
export lookup_function, link, link_crt
export get_buffer

include("utils.jl")
export jlty2llvmty

# boot
const BOOT_COMPILER_REF = Ref{CXCompiler}()

function __init__()
@static if haskey(ENV, "CLANGCOMPILER_ENABLE_BOOT")
    llvm_include_dir = joinpath(LLVM_full_jll.artifact_dir, "include") |> normpath
    @assert isdir(llvm_include_dir) "failed to find LLVM include dir."

    boot_include_dir = joinpath(@__DIR__, "..", "boot") |> normpath
    boot_src = joinpath(boot_include_dir, "boot.cpp")

    args = get_compiler_args(; version=v"7.1.0")
    push!(args, "-std=c++14")
    Sys.isapple() && push!(args, "-stdlib=libc++")
    push!(args, "-I$llvm_include_dir")
    push!(args, "-I$libclangex_include")
    push!(args, "-I$julia_include_dir")
    push!(args, "-I$boot_include_dir")

    jit = LLJIT(;tm=JITTargetMachine())
    irgen = IRGenerator(boot_src, args)
    BOOT_COMPILER_REF[] = CXCompiler(irgen, jit)

    link_process_symbols(BOOT_COMPILER_REF[])
    compile(BOOT_COMPILER_REF[])

    include(joinpath(@__DIR__, "boot.jl"))
end # @static
end

end


text:
function []=PlotFaceGrid(Grid,Subs,Trans,fig)
figure(fig)
clf(fig)
hold on;
dS=2/Subs;
ksi=zeros(2,1);
for iF=1:size(Grid.Faces,2)
  ksi(1)=-1;
  for is=1:Subs
    ksi(2)=-1;
    for js=1:Subs
      [J,D,X1]=Trans(ksi(1),ksi(2),Grid.Faces(iF),Grid);
      [J,D,X2]=Trans(ksi(1)+dS,ksi(2),Grid.Faces(iF),Grid);
      [J,D,X3]=Trans(ksi(1)+dS,ksi(2)+dS,Grid.Faces(iF),Grid);
      [J,D,X4]=Trans(ksi(1),ksi(2)+dS,Grid.Faces(iF),Grid);

      x=[X1(1) X2(1) X3(1) X4(1) X1(1)];
      y=[X1(2) X2(2) X3(2) X4(2) X1(2)];
      if size(X1,3)==2
        fill(x,y,1);
      else
        z=[X1(3) X2(3) X3(3) X4(3) X1(3)];
        fill3(x,y,z,1);
      end
      ksi(2)=ksi(2)+dS;
    end
    ksi(1)=ksi(1)+dS;
  end
end
hold off
end


text:
function compile(model::Model)
    observe = mapcatexpr(s -> :(observations[$(quot(s))][observationid] = $s), model.observed)
    declareobserved = mapcatexpr(var -> :($var = nothing), model.observed)
    declareconds = mapcatexpr(cond -> :($(condname(cond.first)) = $(cond.second)), model.conditionings)

    body = postwalk(model.body) do x
        if @capture(x, var_[idx_] ~ dist_)
            varbody(model, dist, var, brackets(idx)) >> footer(observe)
        elseif @capture(x, var_ ~ dist_)
            varbody(model, dist, var, identity) >> footer(observe)
        elseif @capture(x, return ret_)
            observe >> :(return traceroot)
        else x
        end
    end

    r = quote
        function $(gensym())(traceroot::Trace, observations::Dict{Symbol,Array}, observationid::Int)
            $(model.params) = $(model.arguments)
            $declareconds
            $declareobserved

            loglikelihood = 1.0
            previous = traceroot
            current = traceroot

            $body
            $observe
            return traceroot
        end
    end
    return eval(r)
end


text:
using MathLink
using MathLinkExtras

using MajoranaReps

using Test

@test W"a"*bx(1) == W"a"*bx(1)
@test W"a"*W"b"*bx(1) == W`a b`*bx(1)

@test W"a"*bx(1)*Ket() == W"a"*bx(1)*Ket()
@test W"a"*W"b"*bx(1)*Ket() == W`a b`*bx(1)*Ket()

####At the end also test the usual logic
include("TestMajoranas.jl")


text:
#=

 Unit tests from C-XSC version 2.5.4
 converted into portable ITL format by Oliver Heimlich.

 Copyright 1990-2000 Institut fuer Angewandte Mathematik,
                     Universitaet Karlsruhe, Germany
 Copyright 2000-2014 Wiss. Rechnen/Softwaretechnologie
                     Universitaet Wuppertal, Germany
 Copyright 2015-2016 Oliver Heimlich (oheim@posteo.de)

 This library is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public
 License as published by the Free Software Foundation; either
 version 2.1 of the License, or (at your option) any later version.

 This library is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public
 License along with this library; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

# Tests A+B, B+A, A-B, B-A, -A, +A
@testset "cxsc.intervaladdsub" begin
    @test Interval(10.0, 20.0) + Interval(13.0, 17.0) == Interval(23.0, 37.0)
    @test Interval(13.0, 17.0) + Interval(10.0, 20.0) == Interval(23.0, 37.0)
    @test Interval(10.0, 20.0) - Interval(13.0, 16.0) == Interval(-6.0, 7.0)
    @test Interval(13.0, 16.0) - Interval(10.0, 20.0) == Interval(-7.0, 6.0)
    @test -(Interval(10.0, 20.0)) == Interval(-20.0, -10.0)
    @test +Interval(10.0, 20.0) == Interval(10.0, 20.0)
end

# Tests A*B, B*A, A/B, B/A
@testset "cxsc.intervalmuldiv" begin
    @test Interval(1.0, 2.0) * Interval(3.0, 4.0) == Interval(3.0, 8.0)
    @test Interval(-1.0, 2.0) * Interval(3.0, 4.0) == Interval(-4.0, 8.0)
    @test Interval(-2.0, 1.0) * Interval(3.0, 4.0) == Interval(-8.0, 4.0)
    @test Interval(-2.0, -1.0) * Interval(3.0, 4.0) == Interval(-8.0, -3.0)
    @test Interval(1.0, 2.0) * Interval(-3.0, 4.0) == Interval(-6.0, 8.0)
    @test Interval(-1.0, 2.0) * Interval(-3.0, 4.0) == Interval(-6.0, 8.0)
    @test Interval(-2.0, 1.0) * Interval(-3.0, 4.0) == Interval(-8.0, 6.0)
    @test Interval(-2.0, -1.0) * Interval(-3.0, 4.0) == Interval(-8.0, 6.0)
    @test Interval(1.0, 2.0) * Interval(-4.0, 3.0) == Interval(-8.0, 6.0)
    @test Interval(-1.0, 2.0) * Interval(-4.0, 3.0) == Interval(-8.0, 6.0)
    @test Interval(-2.0, 1.0) * Interval(-4.0, 3.0) == Interval(-6.0, 8.0)
    @test Interval(-2.0, -1.0) * Interval(-4.0, 3.0) == Interval(-6.0, 8.0)
    @test Interval(1.0, 2.0) * Interval(-4.0, -3.0) == Interval(-8.0, -3.0)
    @test Interval(-1.0, 2.0) * Interval(-4.0, -3.0) == Interval(-8.0, 4.0)
    @test Interval(-2.0, -1.0) * Interval(-4.0, -3.0) == Interval(3.0, 8.0)
    @test Interval(1.0, 2.0) / Interval(4.0, 8.0) == Interval(0.125, 0.5)
    @test Interval(-1.0, 2.0) / Interval(4.0, 8.0) == Interval(-0.25, 0.5)
    @test Interval(-2.0, 1.0) / Interval(4.0, 8.0) == Interval(-0.5, 0.25)
    @test Interval(-2.0, -1.0) / Interval(4.0, 8.0) == Interval(-0.5, -0.125)
    @test Interval(1.0, 2.0) / Interval(-4.0, 8.0) == entireinterval(Float64)
    @test Interval(-1.0, 2.0) / Interval(-4.0, 8.0) == entireinterval(Float64)
    @test Interval(-2.0, 1.0) / Interval(-4.0, 8.0) == entireinterval(Float64)
    @test Interval(-2.0, -1.0) / Interval(-4.0, 8.0) == entireinterval(Float64)
    @test Interval(1.0, 2.0) / Interval(-8.0, 4.0) == entireinterval(Float64)
    @test Interval(-1.0, 2.0) / Interval(-8.0, 4.0) == entireinterval(Float64)
    @test Interval(-2.0, 1.0) / Interval(-8.0, 4.0) == entireinterval(Float64)
    @test Interval(-2.0, -1.0) / Interval(-8.0, 4.0) == entireinterval(Float64)
    @test Interval(1.0, 2.0) / Interval(-8.0, -4.0) == Interval(-0.5, -0.125)
    @test Interval(-1.0, 2.0) / Interval(-8.0, -4.0) == Interval(-0.5, 0.25)
    @test Interval(-2.0, 1.0) / Interval(-8.0, -4.0) == Interval(-0.25, 0.5)
    @test Interval(-2.0, -1.0) / Interval(-8.0, -4.0) == Interval(0.125, 0.5)
end

# Tests A|B, B|A, A&B, B&A
@testset "cxsc.intervalsetops" begin
    @test hull(Interval(-2.0, 2.0), Interval(-4.0, -3.0)) == Interval(-4.0, 2.0)
    @test ∪(Interval(-2.0, 2.0), Interval(-4.0, -3.0)) == Interval(-4.0, 2.0)
    @test hull(Interval(-2.0, 2.0), Interval(-4.0, -1.0)) == Interval(-4.0, 2.0)
    @test ∪(Interval(-2.0, 2.0), Interval(-4.0, -1.0)) == Interval(-4.0, 2.0)
    @test hull(Interval(-2.0, 2.0), Interval(-4.0, 4.0)) == Interval(-4.0, 4.0)
    @test ∪(Interval(-2.0, 2.0), Interval(-4.0, 4.0)) == Interval(-4.0, 4.0)
    @test hull(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == Interval(-2.0, 2.0)
    @test ∪(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == Interval(-2.0, 2.0)
    @test hull(Interval(-2.0, 2.0), Interval(1.0, 4.0)) == Interval(-2.0, 4.0)
    @test ∪(Interval(-2.0, 2.0), Interval(1.0, 4.0)) == Interval(-2.0, 4.0)
    @test hull(Interval(-2.0, 2.0), Interval(3.0, 4.0)) == Interval(-2.0, 4.0)
    @test ∪(Interval(-2.0, 2.0), Interval(3.0, 4.0)) == Interval(-2.0, 4.0)
    @test hull(Interval(-4.0, -3.0), Interval(-2.0, 2.0)) == Interval(-4.0, 2.0)
    @test ∪(Interval(-4.0, -3.0), Interval(-2.0, 2.0)) == Interval(-4.0, 2.0)
    @test hull(Interval(-4.0, -1.0), Interval(-2.0, 2.0)) == Interval(-4.0, 2.0)
    @test ∪(Interval(-4.0, -1.0), Interval(-2.0, 2.0)) == Interval(-4.0, 2.0)
    @test hull(Interval(-4.0, 4.0), Interval(-2.0, 2.0)) == Interval(-4.0, 4.0)
    @test ∪(Interval(-4.0, 4.0), Interval(-2.0, 2.0)) == Interval(-4.0, 4.0)
    @test hull(Interval(-1.0, 1.0), Interval(-2.0, 2.0)) == Interval(-2.0, 2.0)
    @test ∪(Interval(-1.0, 1.0), Interval(-2.0, 2.0)) == Interval(-2.0, 2.0)
    @test hull(Interval(1.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 4.0)
    @test ∪(Interval(1.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 4.0)
    @test hull(Interval(3.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 4.0)
    @test ∪(Interval(3.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 4.0)
    @test intersect(Interval(-2.0, 2.0), Interval(-4.0, -3.0)) == ∅
    @test ∩(Interval(-2.0, 2.0), Interval(-4.0, -3.0)) == ∅
    @test intersect(Interval(-2.0, 2.0), Interval(-4.0, -1.0)) == Interval(-2.0, -1.0)
    @test ∩(Interval(-2.0, 2.0), Interval(-4.0, -1.0)) == Interval(-2.0, -1.0)
    @test intersect(Interval(-2.0, 2.0), Interval(-4.0, 4.0)) == Interval(-2.0, 2.0)
    @test ∩(Interval(-2.0, 2.0), Interval(-4.0, 4.0)) == Interval(-2.0, 2.0)
    @test intersect(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == Interval(-1.0, 1.0)
    @test ∩(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == Interval(-1.0, 1.0)
    @test intersect(Interval(-2.0, 2.0), Interval(1.0, 4.0)) == Interval(1.0, 2.0)
    @test ∩(Interval(-2.0, 2.0), Interval(1.0, 4.0)) == Interval(1.0, 2.0)
    @test intersect(Interval(-2.0, 2.0), Interval(3.0, 4.0)) == ∅
    @test ∩(Interval(-2.0, 2.0), Interval(3.0, 4.0)) == ∅
    @test intersect(Interval(-4.0, -3.0), Interval(-2.0, 2.0)) == ∅
    @test ∩(Interval(-4.0, -3.0), Interval(-2.0, 2.0)) == ∅
    @test intersect(Interval(-4.0, -1.0), Interval(-2.0, 2.0)) == Interval(-2.0, -1.0)
    @test ∩(Interval(-4.0, -1.0), Interval(-2.0, 2.0)) == Interval(-2.0, -1.0)
    @test intersect(Interval(-4.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 2.0)
    @test ∩(Interval(-4.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 2.0)
    @test intersect(Interval(-1.0, 1.0), Interval(-2.0, 2.0)) == Interval(-1.0, 1.0)
    @test ∩(Interval(-1.0, 1.0), Interval(-2.0, 2.0)) == Interval(-1.0, 1.0)
    @test intersect(Interval(1.0, 4.0), Interval(-2.0, 2.0)) == Interval(1.0, 2.0)
    @test ∩(Interval(1.0, 4.0), Interval(-2.0, 2.0)) == Interval(1.0, 2.0)
    @test intersect(Interval(3.0, 4.0), Interval(-2.0, 2.0)) == ∅
    @test ∩(Interval(3.0, 4.0), Interval(-2.0, 2.0)) == ∅
end

# Tests A|B, B|A, A&B, B&A, B is scalar-type
@testset "cxsc.intervalmixsetops" begin
    @test hull(Interval(-2.0, 2.0), Interval(-4.0, -4.0)) == Interval(-4.0, 2.0)
    @test ∪(Interval(-2.0, 2.0), Interval(-4.0, -4.0)) == Interval(-4.0, 2.0)
    @test hull(Interval(-2.0, 2.0), Interval(1.0, 1.0)) == Interval(-2.0, 2.0)
    @test ∪(Interval(-2.0, 2.0), Interval(1.0, 1.0)) == Interval(-2.0, 2.0)
    @test hull(Interval(-2.0, 2.0), Interval(4.0, 4.0)) == Interval(-2.0, 4.0)
    @test ∪(Interval(-2.0, 2.0), Interval(4.0, 4.0)) == Interval(-2.0, 4.0)
    @test hull(Interval(-4.0, -4.0), Interval(-2.0, 2.0)) == Interval(-4.0, 2.0)
    @test ∪(Interval(-4.0, -4.0), Interval(-2.0, 2.0)) == Interval(-4.0, 2.0)
    @test hull(Interval(1.0, 1.0), Interval(-2.0, 2.0)) == Interval(-2.0, 2.0)
    @test ∪(Interval(1.0, 1.0), Interval(-2.0, 2.0)) == Interval(-2.0, 2.0)
    @test hull(Interval(4.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 4.0)
    @test ∪(Interval(4.0, 4.0), Interval(-2.0, 2.0)) == Interval(-2.0, 4.0)
    @test intersect(Interval(-2.0, 2.0), Interval(-4.0, -4.0)) == ∅
    @test ∩(Interval(-2.0, 2.0), Interval(-4.0, -4.0)) == ∅
    @test intersect(Interval(-2.0, 2.0), Interval(1.0, 1.0)) == Interval(1.0, 1.0)
    @test ∩(Interval(-2.0, 2.0), Interval(1.0, 1.0)) == Interval(1.0, 1.0)
    @test intersect(Interval(-2.0, 2.0), Interval(4.0, 4.0)) == ∅
    @test ∩(Interval(-2.0, 2.0), Interval(4.0, 4.0)) == ∅
    @test intersect(Interval(-4.0, -4.0), Interval(-2.0, 2.0)) == ∅
    @test ∩(Interval(-4.0, -4.0), Interval(-2.0, 2.0)) == ∅
    @test intersect(Interval(1.0, 1.0), Interval(-2.0, 2.0)) == Interval(1.0, 1.0)
    @test ∩(Interval(1.0, 1.0), Interval(-2.0, 2.0)) == Interval(1.0, 1.0)
    @test intersect(Interval(4.0, 4.0), Interval(-2.0, 2.0)) == ∅
    @test ∩(Interval(4.0, 4.0), Interval(-2.0, 2.0)) == ∅
end

# Tests A|B, B|A, A and B are scalar-type
@testset "cxsc.scalarmixsetops" begin
    @test hull(Interval(-2.0, -2.0), Interval(-4.0, -4.0)) == Interval(-4.0, -2.0)
    @test ∪(Interval(-2.0, -2.0), Interval(-4.0, -4.0)) == Interval(-4.0, -2.0)
    @test hull(Interval(-2.0, -2.0), Interval(-2.0, -2.0)) == Interval(-2.0, -2.0)
    @test ∪(Interval(-2.0, -2.0), Interval(-2.0, -2.0)) == Interval(-2.0, -2.0)
    @test hull(Interval(-2.0, -2.0), Interval(2.0, 2.0)) == Interval(-2.0, 2.0)
    @test ∪(Interval(-2.0, -2.0), Interval(2.0, 2.0)) == Interval(-2.0, 2.0)
    @test hull(Interval(-4.0, -4.0), Interval(-2.0, -2.0)) == Interval(-4.0, -2.0)
    @test ∪(Interval(-4.0, -4.0), Interval(-2.0, -2.0)) == Interval(-4.0, -2.0)
    @test hull(Interval(-2.0, -2.0), Interval(-2.0, -2.0)) == Interval(-2.0, -2.0)
    @test ∪(Interval(-2.0, -2.0), Interval(-2.0, -2.0)) == Interval(-2.0, -2.0)
    @test hull(Interval(2.0, 2.0), Interval(-2.0, -2.0)) == Interval(-2.0, 2.0)
    @test ∪(Interval(2.0, 2.0), Interval(-2.0, -2.0)) == Interval(-2.0, 2.0)
end

# Tests A<B, A>B, A<=B, A>=B, A==B
@testset "cxsc.intervalsetcompops" begin
    @test ⪽(Interval(-1.0, 2.0), Interval(-1.0, 2.0)) == false
    @test isinterior(Interval(-1.0, 2.0), Interval(-1.0, 2.0)) == false
    @test ⪽(Interval(-2.0, 1.0), Interval(-3.0, 2.0))
    @test isinterior(Interval(-2.0, 1.0), Interval(-3.0, 2.0))
    @test ⪽(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(-1.0, 2.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-1.0, 2.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(-2.0, 1.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-2.0, 1.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(-2.0, 3.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-2.0, 3.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(-3.0, 2.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-3.0, 2.0)) == false
    @test ⪽(Interval(-1.0, 2.0), Interval(-1.0, 2.0)) == false
    @test isinterior(Interval(-1.0, 2.0), Interval(-1.0, 2.0)) == false
    @test ⪽(Interval(-3.0, 2.0), Interval(-2.0, 1.0)) == false
    @test isinterior(Interval(-3.0, 2.0), Interval(-2.0, 1.0)) == false
    @test ⪽(Interval(-1.0, 1.0), Interval(-2.0, 2.0))
    @test isinterior(Interval(-1.0, 1.0), Interval(-2.0, 2.0))
    @test ⪽(Interval(-1.0, 2.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(-1.0, 2.0), Interval(-2.0, 2.0)) == false
    @test ⪽(Interval(-2.0, 1.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(-2.0, 1.0), Interval(-2.0, 2.0)) == false
    @test ⪽(Interval(-2.0, 3.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(-2.0, 3.0), Interval(-2.0, 2.0)) == false
    @test ⪽(Interval(-3.0, 2.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(-3.0, 2.0), Interval(-2.0, 2.0)) == false
    @test ⊆(Interval(-1.0, 2.0), Interval(-1.0, 2.0))
    @test ⊆(Interval(-2.0, 1.0), Interval(-3.0, 2.0))
    @test ⊆(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(-1.0, 2.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(-2.0, 1.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(-2.0, 3.0))
    @test ⊆(Interval(-2.0, 2.0), Interval(-3.0, 2.0))
    @test ⊆(Interval(-3.0, 2.0), Interval(-2.0, 1.0)) == false
    @test ⊆(Interval(-1.0, 1.0), Interval(-2.0, 2.0))
    @test ⊆(Interval(-1.0, 2.0), Interval(-2.0, 2.0))
    @test ⊆(Interval(-2.0, 1.0), Interval(-2.0, 2.0))
    @test ⊆(Interval(-2.0, 3.0), Interval(-2.0, 2.0)) == false
    @test ⊆(Interval(-3.0, 2.0), Interval(-2.0, 2.0)) == false
    @test ==(Interval(-1.0, 2.0), Interval(-1.0, 2.0))
    @test ==(Interval(-2.0, 1.0), Interval(-3.0, 2.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(-1.0, 1.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(-1.0, 2.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(-2.0, 1.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(-2.0, 3.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(-3.0, 2.0)) == false
end

# Tests A<B, A>B, A<=B, A>=B, A==B, B<A, B>A, B<=A, B>=A, B==A, where B is scalar
@testset "cxsc.intervalscalarsetcompops" begin
    @test ⪽(Interval(-1.0, 2.0), Interval(-2.0, -2.0)) == false
    @test isinterior(Interval(-1.0, 2.0), Interval(-2.0, -2.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(-2.0, -2.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(-2.0, -2.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(0.0, 0.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(0.0, 0.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(2.0, 2.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(2.0, 2.0)) == false
    @test ⪽(Interval(-2.0, 2.0), Interval(3.0, 3.0)) == false
    @test isinterior(Interval(-2.0, 2.0), Interval(3.0, 3.0)) == false
    @test ⪽(Interval(-1.0, -1.0), Interval(1.0, 1.0)) == false
    @test isinterior(Interval(-1.0, -1.0), Interval(1.0, 1.0)) == false
    @test ⪽(Interval(-1.0, -1.0), Interval(-1.0, -1.0)) == false
    @test isinterior(Interval(-1.0, -1.0), Interval(-1.0, -1.0)) == false
    @test ⪽(Interval(-2.0, -2.0), Interval(-1.0, 2.0)) == false
    @test isinterior(Interval(-2.0, -2.0), Interval(-1.0, 2.0)) == false
    @test ⪽(Interval(-2.0, -2.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(-2.0, -2.0), Interval(-2.0, 2.0)) == false
    @test ⪽(Interval(0.0, 0.0), Interval(-2.0, 2.0))
    @test isinterior(Interval(0.0, 0.0), Interval(-2.0, 2.0))
    @test ⪽(Interval(2.0, 2.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(2.0, 2.0), Interval(-2.0, 2.0)) == false
    @test ⪽(Interval(3.0, 3.0), Interval(-2.0, 2.0)) == false
    @test isinterior(Interval(3.0, 3.0), Interval(-2.0, 2.0)) == false
    @test ⪽(Interval(1.0, 1.0), Interval(-1.0, -1.0)) == false
    @test isinterior(Interval(1.0, 1.0), Interval(-1.0, -1.0)) == false
    @test ⪽(Interval(-1.0, -1.0), Interval(-1.0, -1.0)) == false
    @test isinterior(Interval(-1.0, -1.0), Interval(-1.0, -1.0)) == false
    @test ⊆(Interval(-1.0, 2.0), Interval(-2.0, -2.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(-2.0, -2.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(0.0, 0.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(2.0, 2.0)) == false
    @test ⊆(Interval(-2.0, 2.0), Interval(3.0, 3.0)) == false
    @test ⊆(Interval(-1.0, -1.0), Interval(1.0, 1.0)) == false
    @test ⊆(Interval(-1.0, -1.0), Interval(-1.0, -1.0))
    @test ⊆(Interval(-2.0, -2.0), Interval(-1.0, 2.0)) == false
    @test ⊆(Interval(-2.0, -2.0), Interval(-2.0, 2.0))
    @test ⊆(Interval(0.0, 0.0), Interval(-2.0, 2.0))
    @test ⊆(Interval(2.0, 2.0), Interval(-2.0, 2.0))
    @test ⊆(Interval(3.0, 3.0), Interval(-2.0, 2.0)) == false
    @test ⊆(Interval(1.0, 1.0), Interval(-1.0, -1.0)) == false
    @test ⊆(Interval(-1.0, -1.0), Interval(-1.0, -1.0))
    @test ==(Interval(-1.0, 2.0), Interval(-2.0, -2.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(-2.0, -2.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(0.0, 0.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(2.0, 2.0)) == false
    @test ==(Interval(-2.0, 2.0), Interval(3.0, 3.0)) == false
    @test ==(Interval(-1.0, -1.0), Interval(1.0, 1.0)) == false
    @test ==(Interval(-1.0, -1.0), Interval(-1.0, -1.0))
end

@testset "cxsc.intervalstdfunc" begin
    @test Interval(11.0, 11.0) ^2 == Interval(121.0, 121.0)
    @test Interval(11.0, 11.0) ^(2//1) == Interval(121.0, 121.0)
    @test Interval(11.0, 11.0) ^(2.0) == Interval(121.0, 121.0)
    @test Interval(0.0, 0.0) ^2 == Interval(0.0, 0.0)
    @test Interval(0.0, 0.0) ^(2//1) == Interval(0.0, 0.0)
    @test Interval(0.0, 0.0) ^(2.0) == Interval(0.0, 0.0)
    @test Interval(-9.0, -9.0) ^2 == Interval(81.0, 81.0)
    @test Interval(-9.0, -9.0) ^(2//1) == Interval(81.0, 81.0)
    @test Interval(-9.0, -9.0) ^(2.0) == Interval(81.0, 81.0)
    @test sqrt(Interval(121.0, 121.0)) == Interval(11.0, 11.0)
    @test (Interval(121.0, 121.0))^(1/2) == Interval(11.0, 11.0)
    @test (Interval(121.0, 121.0))^(0.5) == Interval(11.0, 11.0)
    @test (Interval(121.0, 121.0))^(1//2) == Interval(11.0, 11.0)
    @test sqrt(Interval(0.0, 0.0)) == Interval(0.0, 0.0)
    @test (Interval(0.0, 0.0))^(1/2) == Interval(0.0, 0.0)
    @test (Interval(0.0, 0.0))^(0.5) == Interval(0.0, 0.0)
    @test (Interval(0.0, 0.0))^(1//2) == Interval(0.0, 0.0)
    @test sqrt(Interval(81.0, 81.0)) == Interval(9.0, 9.0)
    @test (Interval(81.0, 81.0))^(1/2) == Interval(9.0, 9.0)
    @test (Interval(81.0, 81.0))^(0.5) == Interval(9.0, 9.0)
    @test (Interval(81.0, 81.0))^(1//2) == Interval(9.0, 9.0)
    @test Interval(2.0, 2.0) ^ Interval(2.0, 2.0) == Interval(4.0, 4.0)
    @test Interval(4.0, 4.0) ^ Interval(5.0, 5.0) == Interval(1024.0, 1024.0)
    # Negativ geht noch nicht
    @test Interval(2.0, 2.0) ^ Interval(3.0, 3.0) == Interval(8.0, 8.0)
end
# FactCheck.exitstatus()


text:
#holds prediction information
mutable struct VARPrediction
  dates::Vector{T} where T<:Union{Date, Int, Float64}
  X::Matrix{Float64}
  Y::Matrix{Float64}
  YHat::Matrix{Float64}
end

#holds fundamental information on the VAR
mutable struct VAR
  ϕ₀::Vector{Float64}
  ϕ₁::Matrix{Float64}
  compactϕ₀::Vector{Float64}
  compactϕ₁::Matrix{Float64}
  μ::Vector{Float64}
  L::Matrix{Float64}
  structures::Vector{T}  where T <: Structure
end

#holds the VAR and any estimation information
mutable struct VAREstimate
  var::VAR
  coeffErrors::Vector{T} where T <: CoefficientErrors
  N::Vector{T} where T <: RegressionN
  SSE::Vector{T} where T <: RegressionSS
  SSR::Vector{T} where T <: RegressionSS
  SST::Vector{T} where T <: RegressionSS
end

function lagDF!(DF::DataFrame, target::Symbol, lags::Int = 1; name= Symbol("$(target)L$lags"))::DataFrame
  DF[:,name] = similar(DF[:,target])
  DF[:,name] .= missing
  DF[(lags+1):end,name] = DF[1:(end-lags),target]

  return DF
end

function VAR(ϕ₀::Vector{Float64}, ϕ₁::Matrix{Float64}, L::Matrix{Float64},
    structures::Vector{T}  where T <: Structure)::VAR

  #calculate the mean
  K::Int = length(ϕ₀)
  μ::Vector{Float64} = ((Matrix{Float64}(I,K,K) .- ϕ₁)\Matrix{Float64}(I,K,K))*ϕ₀

  compactRows::BitArray = structures .== nothing

  compactϕ₀::Vector{Float64} = ϕ₀[compactRows]
  compactϕ₁::Matrix{Float64} = ϕ₁[compactRows,compactRows]

  return VAR(ϕ₀, ϕ₁, compactϕ₀, compactϕ₁, μ,  L, structures)
end

#not sure how well this works with restrictions, if at all
function propagateImpulse(var::VAR, Z::Vector{Float64}, lags::Int=20)::Matrix{Float64}
  compactShocks::Vector{Float64} = var.L * Z
  K::Int = length(var.ϕ₀)
  propagate::Matrix{Float64} = zeros(K, lags + 1)

  #expand the shock vector
  compactRows::BitArray = var.structures .== nothing
  propagate[compactRows,1] .= compactShocks

  #this accounts for the demeaning

  for i::Int ∈ 1:lags
    propagate[:, i+1] = var.ϕ₁*propagate[:,i]
  end

  return propagate
end

  #get the mean

#make period or multi-period ahead perdictions
function VARPrediction(var::VAR, dates::Vector{D}, X::Matrix{Float64}, Y::Matrix{Float64};
    τ::Int = 1)::VARPrediction where D <: Union{Date, Int, Float64}

  K::Int = length(var.ϕ₀)
  T::Int = size(Y,2)

  #housekeeping with lengths
  Y=Y[:,(τ+1):end]
  dates = dates[(τ+1):end]
  YHat::Matrix{Float64} = Matrix{Float64}(K, T-τ)
  ε::Matrix{Float64} = similar(YHat)
  predict::Vector{Float64} = Vector{Float64}(K)

  for i::Int ∈ 1:(T-τ)
    predict .= var.ϕ₀ .+ var.ϕ₁ * X[:,i] #first period ahead is based on data

    for j::Int ∈ 2:τ #subsequent periods based on previous prediction
      predict .= var.ϕ₀ .+ var.ϕ₁ * predict
    end
    YHat[:,i] .= predict
    ε[:,i] = Y[:, i] .- predict
  end

  return VARPrediction(dates, X, Y, YHat)
end

function varΣ(var::VAR)::Matrix{Float64}
  K::Int = length(var.compactϕ₀)

  #maybe a more efficient way to do this?
  bigMat::Matrix{Float64} = (eye(K^2) .- kron(var.compactϕ₁, var.compactϕ₁))\eye(K^2)
  Γ₀::Matrix{Float64} = reshape(bigMat*vec(var.L*var.L'),K,:)

  return Γ₀
end

#this is the low-level VAR estimator
function VAREstimate(DF::DataFrame, LHS::Vector{Symbol}, RHS::Vector{Symbol};
      errors::Function = getHomoskedΣ!, K::Int = length(LHS),
      structures::Vector{Structure}  = Vector{Structure}(fill(nothing, K)),
      ϕ::Matrix{Float64}=Matrix{Float64}(undef, K,K+1), #these are available for pre-allocaiton purposes
      σ::Vector{Float64} = Vector{Float64}(undef, K))

  #build up a couple of indices and helper variables
  K::Int = length(LHS)
  RHSIndex::Dict = Dict(RHS[i]=>i for i::Int ∈ 1:K)
  LHSIndex::Dict = Dict(LHS[i]=>i for i::Int ∈ 1:K)

  coeffErrors::Vector{CoefficientErrors} = Vector{Structure}(fill(nothing, K))
  N::Vector{RegressionN} = Vector{RegressionN}(fill(nothing, K))
  ε::Vector{Vector{Float64}} = Vector{Vector{Float64}}()

  #store the ANOVA information
  SSE::Vector{RegressionSS} = Vector{RegressionSS}(fill(nothing, K))
  SSR::Vector{RegressionSS} = similar(SSE)
  SST::Vector{RegressionSS} = similar(SSE)

  XExpr::CTExpr = Meta.parse(vec2String(RHS, " + "))

  #check for structure. If its not a structure (either specified or implied as a lag),
  #run a regression and estimate the coefficients
  for k::Int ∈ 1:K
    if structures[k] ≠ nothing
      ϕ[k,:] .= structures[k]
      σ[k] = 0.0
      coeffErrors[k], N[k], ε[k], SSE[k], SSR[k], SST[k]= fill(nothing, 6)
    elseif haskey(RHSIndex,LHS[k]) #by defintion, the lag must be restricted
      ϕ[k,:] .= 0.0
      ϕ[k,RHSIndex[LHS[k]]+1] = 1.0 #lookup the location and restict the outcome to an identity
      structures[k] = Vector{Float64}(ϕ[k,:])
      σ[k] = 0.0
      coeffErrors[k], N[k], ε[k], SSE[k], SSR[k], SST[k] = fill(nothing, 6)
    else
      structures[k] = nothing #probably redundant
      reg::CTLM = CTLM(DF,  XExpr, LHS[k], XNames=[:Intercept; RHS], YName=LHS[k])
      ϕ[k,:] .= reg.β
      push!(ε, reg.ε)
      σ[k] = std(ε[end])
      coeffErrors[k] = diag(getHomoskedΣ!(reg)).^0.5
      N[k] = reg.N
      SSE[k] = sum(reg.ε.^2)
      SSR[k] = sum((reg.Y .- reg.ε .- mean(reg.Y)).^2)
      SST[k] = sum((reg.Y .- mean(reg.Y)).^2)
    end
  end

  L::Matrix{Float64} = (cholesky(cov(hcat(ε...)))).L #(splat and calculate the cholesky root)

  var::VAR = VAR(Vector{Float64}(ϕ[:,1]), ϕ[:,2:end], L, structures)

  #note this only works for the linearly independent components

  return   VAREstimate(var, coeffErrors, N, SSE, SSR, SST)
end

max_stars_repo_path: Fall 17 Ivo Project 2.2/histtaxes/CTCode/CTVAR.jl