using StatsBase, Random, NSGAII

# include("timeseries.jl")

import Base: copyto!, isequal, iterate, getindex, length, ndims, size

# Abstract operator type.
abstract type AbstractOperator end

# Abstract unary operator type.
abstract type AbstractUnaryOperator                 <: AbstractOperator end
# Abstract binary operator type.
abstract type AbstractBinaryOperator                <: AbstractOperator end

# Abstract unary modal operator type.
abstract type AbstractUnaryModalOperator            <: AbstractUnaryOperator end

# Abstract binary boolean operator type.
abstract type AbstractBinaryBooleanOperator         <: AbstractBinaryOperator end

# Abstract unary interval operator type.
abstract type AbstractUnaryIntervalOperator         <: AbstractUnaryModalOperator end

# Abstract existential (unary) interval operator type.
abstract type AbstractExistentialIntervalOperator   <: AbstractUnaryIntervalOperator end
# Abstract universal (unary) interval operator type.
abstract type AbstractUniversalIntervalOperator     <: AbstractUnaryIntervalOperator end

# Concrete existential (unary) interval operator types.
struct ExistentialGlobally          <: AbstractExistentialIntervalOperator end
struct ExistentialLater             <: AbstractExistentialIntervalOperator end
struct ExistentialAfter             <: AbstractExistentialIntervalOperator end
struct ExistentialOverlaps          <: AbstractExistentialIntervalOperator end
struct ExistentialEnds              <: AbstractExistentialIntervalOperator end
struct ExistentialDuring            <: AbstractExistentialIntervalOperator end
struct ExistentialBegins            <: AbstractExistentialIntervalOperator end
struct ExistentialInverseLater      <: AbstractExistentialIntervalOperator end
struct ExistentialInverseAfter      <: AbstractExistentialIntervalOperator end
struct ExistentialInverseOverlaps   <: AbstractExistentialIntervalOperator end
struct ExistentialInverseEnds       <: AbstractExistentialIntervalOperator end
struct ExistentialInverseDuring     <: AbstractExistentialIntervalOperator end
struct ExistentialInverseBegins     <: AbstractExistentialIntervalOperator end
# Concrete universal (unary) interval operator types.
struct UniversalGlobally            <: AbstractUniversalIntervalOperator end
struct UniversalLater               <: AbstractUniversalIntervalOperator end
struct UniversalAfter               <: AbstractUniversalIntervalOperator end
struct UniversalOverlaps            <: AbstractUniversalIntervalOperator end
struct UniversalEnds                <: AbstractUniversalIntervalOperator end
struct UniversalDuring              <: AbstractUniversalIntervalOperator end
struct UniversalBegins              <: AbstractUniversalIntervalOperator end
struct UniversalInverseLater        <: AbstractUniversalIntervalOperator end
struct UniversalInverseAfter        <: AbstractUniversalIntervalOperator end
struct UniversalInverseOverlaps     <: AbstractUniversalIntervalOperator end
struct UniversalInverseEnds         <: AbstractUniversalIntervalOperator end
struct UniversalInverseDuring       <: AbstractUniversalIntervalOperator end
struct UniversalInverseBegins       <: AbstractUniversalIntervalOperator end

# abstract type AbstractBound end

# struct Top      <: AbstractBound end
# struct Bottom   <: AbstractBound end

# Concrete (binary) conjunction type.
struct Conjunction  <: AbstractBinaryBooleanOperator end
# Concrete (binary) disjunction type.
struct Disjunction  <: AbstractBinaryBooleanOperator end
# Concrete (binary) implication type.
struct Implication  <: AbstractBinaryBooleanOperator end

# Constant existential singletons.
const exG = ExistentialGlobally()
const exL = ExistentialLater()
const exA = ExistentialAfter()
const exO = ExistentialOverlaps()
const exE = ExistentialEnds()
const exD = ExistentialDuring()
const exB = ExistentialBegins()
const exL̄ = ExistentialInverseLater()
const exĀ = ExistentialInverseAfter()
const exŌ = ExistentialInverseOverlaps()
const exĒ = ExistentialInverseEnds()
const exD̄ = ExistentialInverseDuring()
const exB̄ = ExistentialInverseBegins()
# Constant universal singletons.
const univG = UniversalGlobally()
const univL = UniversalLater()
const univA = UniversalAfter()
const univO = UniversalOverlaps()
const univE = UniversalEnds()
const univD = UniversalDuring()
const univB = UniversalBegins()
const univL̄ = UniversalInverseLater()
const univĀ = UniversalInverseAfter()
const univŌ = UniversalInverseOverlaps()
const univĒ = UniversalInverseEnds()
const univD̄ = UniversalInverseDuring()
const univB̄ = UniversalInverseBegins()

# const top = Top()
# const bottom = Bottom()

order_relations = [≤, <, ==, >, ≥]
modal_relations = [exL,exA,exO,exE,exD,exB,exL̄,exĀ,exŌ,exĒ,exD̄,exB̄,univL,univA,univO,univE,univD,univB,univL̄,univĀ,univŌ,univĒ,univD̄,univB̄]
granularity_algebra_values = 0.1
algebra_values = [x for x ∈ 0.0:granularity_algebra_values:1.0]
minmd = 1 # max modal depth
maxmd = 2 # max modal depth
minr = 4
maxr = 6
minh = 1 # min horizon
maxh = 6 # max horizon


show(io::IO, ::ExistentialGlobally)        = print(io, "⟨G⟩")
show(io::IO, ::ExistentialLater)           = print(io, "⟨L⟩")
show(io::IO, ::ExistentialAfter)           = print(io, "⟨A⟩")
show(io::IO, ::ExistentialOverlaps)        = print(io, "⟨O⟩")
show(io::IO, ::ExistentialEnds)            = print(io, "⟨E⟩")
show(io::IO, ::ExistentialDuring)          = print(io, "⟨D⟩")
show(io::IO, ::ExistentialBegins)          = print(io, "⟨B⟩")
show(io::IO, ::ExistentialInverseLater)    = print(io, "⟨L̄⟩")
show(io::IO, ::ExistentialInverseAfter)    = print(io, "⟨Ā⟩")
show(io::IO, ::ExistentialInverseOverlaps) = print(io, "⟨Ō⟩")
show(io::IO, ::ExistentialInverseEnds)     = print(io, "⟨Ē⟩")
show(io::IO, ::ExistentialInverseDuring)   = print(io, "⟨D̄⟩")
show(io::IO, ::ExistentialInverseBegins)   = print(io, "⟨B⟩")
show(io::IO, ::UniversalGlobally)          = print(io, "[G]")
show(io::IO, ::UniversalLater)             = print(io, "[L]")
show(io::IO, ::UniversalAfter)             = print(io, "[A]")
show(io::IO, ::UniversalOverlaps)          = print(io, "[O]")
show(io::IO, ::UniversalEnds)              = print(io, "[E]")
show(io::IO, ::UniversalDuring)            = print(io, "[D]")
show(io::IO, ::UniversalBegins)            = print(io, "[B]")
show(io::IO, ::UniversalInverseLater)      = print(io, "[L̄]")
show(io::IO, ::UniversalInverseAfter)      = print(io, "[Ā]")
show(io::IO, ::UniversalInverseOverlaps)   = print(io, "[Ō]")
show(io::IO, ::UniversalInverseEnds)       = print(io, "[Ē]")
show(io::IO, ::UniversalInverseDuring)     = print(io, "[D̄]")
show(io::IO, ::UniversalInverseBegins)     = print(io, "[B̄]")

# Base.show(io::IO, ::Top)    = print("⊤")
# Base.show(io::IO, ::Bottom) = print("⊥")

# Constant conjunction singleton.
const conjunction   = Conjunction()
# Constant disjunction singleton.
const disjunction   = Disjunction()
# Constant implication singleton.
const implication   = Implication()

show(io::IO, ::Conjunction)    = print(io, " ∧ ")
show(io::IO, ::Disjunction)    = print(io, " ∨ ")
show(io::IO, ::Implication)    = print(io, " → ")

abstract type AbstractProposition end

# p = Proposition(1,"≤",10)
mutable struct Proposition <: AbstractProposition
    attribute::Int64
    relation::Function
    constant::Float64
end

show(io::IO, p::Proposition) = print(io, "A", p.attribute, " ", p.relation, " ", p.constant)

abstract type AbstractHeytingAlgebraValue end

mutable struct HeytingAlgebraValue <: AbstractHeytingAlgebraValue
    value::Float64
end

function show(io::IO, v::HeytingAlgebraValue)
    if v.value == 0.0
        print(io, "⊥")
    elseif v.value == 1.0
        print(io, "⊤")
    else
        print(io, v.value)
    end
end

mutable struct Node{T}
    data::T
    parent::Union{Node{AbstractProposition}, Node{AbstractOperator}, Node{AbstractHeytingAlgebraValue}}
    left::Union{Node{AbstractProposition}, Node{AbstractOperator}, Node{AbstractHeytingAlgebraValue}}
    right::Union{Node{AbstractProposition}, Node{AbstractOperator}, Node{AbstractHeytingAlgebraValue}}

    # Root constructor.
    Node{T}(data) where T = new{T}(data)
    # Child node constructor.
    # Node{T}(data, parent) where T = new{T}(data, parent)
end

parent!(n::Node, p::Node)       = n.parent  = p
leftchild!(n::Node, l::Node)    = n.left    = l
rightchild!(n::Node, r::Node)   = n.right   = r

abstract type AbstractRule end

mutable struct ClassificationRule <: AbstractRule
    antecedent::Node
    consequent::String
end

antecedent(ρ::ClassificationRule) = ρ.antecedent
consequent(ρ::ClassificationRule) = ρ.consequent

function size(ρ::ClassificationRule)
    return _size(antecedent(ρ))
end

function _size(n::Node)
    szl = isdefined(n, :left) ? _size(n.left) : 0
    szr = isdefined(n, :right) ? _size(n.right) : 0
    if typeof(n.data) <: Union{AbstractProposition, AbstractHeytingAlgebraValue}
        return 1 + szl + szr
    else
        return 2 + szl + szr
    end
end

function upwards_modal_depth(n::Node)
    md = isdefined(n, :parent) ? upwards_modal_depth(n.parent) : 0
    md += typeof(n.data) <: AbstractUnaryIntervalOperator ? 1 : 0
    return md
end

function is_leaf(n::Node)
    return (!isdefined(n, :left) && !isdefined(n, :right) ? true : false)
end

function is_internal_node(n::Node)
    return (!is_leaf(n) ? true : false)
end

show(io::IO, ρ::ClassificationRule) = print_rule(io, ρ)

mutable struct ClassificationRules
    rules::Vector{ClassificationRule}
    modalities::Vector{AbstractUnaryIntervalOperator}
    horizon::Int # TODO works only for specific functions; for example, in the future, add α and β for sigmoid-like functions

    ClassificationRules(rules;
        modalities::Vector{AbstractUnaryIntervalOperator}=modal_relations,
        horizon::Int=4) = new(rules, modalities, horizon)
end

rules(Γ::ClassificationRules) = Γ.rules
modalities(Γ::ClassificationRules) = Γ.modalities
horizon(Γ::ClassificationRules) = Γ.horizon
# numclasses(crules::ClassificationRules) = crules.numclasses
maxnumrules(Γ::ClassificationRules) = Γ.maxnumrules
getindex(Γ::ClassificationRules, inds...) = getindex(rules(Γ), inds...)
classes(Γ::ClassificationRules) = collect(Set([consequent(ρ) for ρ in Γ]))
complexity(Γ::ClassificationRules) = reduce(+, [size(ρ) for ρ in Γ])
mean_size(Γ::ClassificationRules) = reduce(+,[size(ρ) for ρ in Γ])/length(Γ)

function show(io::IO, Γ::ClassificationRules)
    println(io, "Horizon=", horizon(Γ))
    print_rules(io, Γ)
end

function show(io::IO, n::Node)
    _print_antecedent(io, n)
    println(io)
end

# Defined, otherwise it gives errors
length(crules::ClassificationRules) = length(rules(crules))
iterate(crules::ClassificationRules) = iterate(rules(crules))
iterate(crules::ClassificationRules, i::Int) = iterate(rules(crules), i)
firstindex(Γ::ClassificationRules) = 1
lastindex(Γ::ClassificationRules) = length(Γ)
# ndims(crules::Type{ClassificationRules}) = ndims(rules(crules)) # TODO necessary?
ndims(::Type{ClassificationRules}) = 1
size(crules::ClassificationRules) = (length(crules),)
# copyto!(crules::ClassificationRules, b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(identity),Tuple{Array{Any,1}}}) = copyto!(rules(crules), b)
copyto!(crules::ClassificationRules, b::Base.Broadcast.Broadcasted) = copyto!(rules(crules), b)

# p = Proposition(1,"≤",10)

# x = Node{AbstractOperator}(exL)
# y = Node{AbstractProposition}(p)

# leftchild!(x,y)

# println(x)


# TODO
#   * la proposizione deve essere calcolata sulla base del dominio dell'attributo i-esimo scelto
# function _rand_proposition(numprops::Int, relations::Array{String,1}, domain::Union{Array{Int,1},Array{Real,1}})
#     Proposition(rand(1:numprops), relations[rand(1:length(relations))],domain[rand(1:length(domain))])
# end

function _rand_proposition(numattributes::Int64, relations::Array{Function, 1}, domains::Array{Array{Float64, N} where N, 1})
    i = rand(1:numattributes)
    return Proposition(i, relations[rand(1:length(relations))], domains[i][rand(1:length(domains[i]))])
end

function _rand_proposition(numattributes::Int64, relations::Array{Function, 1}, domains::Dict{Symbol, Vector{Number}})
    i = rand(1:numattributes)
    return Proposition(i, relations[rand(1:length(relations))], domains[Symbol("A" * string(i))][rand(1:length(domains[Symbol("A" * string(i))]))])
end

function rand_rule(;
        modalities::Vector{AbstractUnaryIntervalOperator},
        minmd::Int64,
        maxmd::Int64,
        isroot::Bool)
    # TODO: mettere a parametro la probabilità di ritornare una proposizione
    antecedent = _rand_antecedent(; modalities=modalities, minmd=minmd, maxmd=maxmd, isroot=true)
    consequent = _rand_consequent(; classes=train.unique_classes) # edu
    return ClassificationRule(antecedent, consequent)
end

function _rand_antecedent(;
        modalities::Vector{AbstractUnaryIntervalOperator},
        minmd,
        maxmd,
        isroot)
    if maxmd == 0 || (rand() <= 0.5 && minmd <= 0)
        case = rand(1:100)
        # proposition with 90% probability
        if case ≤ 90
            return Node{AbstractProposition}(_rand_proposition(length(attributes(train)[1]),order_relations,train.domains[1])) # TODO works only for temporal frames, i.e. [1]
        # heyting algebra value with 10% probability
        else
            return Node{AbstractHeytingAlgebraValue}(HeytingAlgebraValue(algebra_values[rand(1:length(algebra_values))]))
        end
    else
        root = isroot ? Node{AbstractOperator}([exG, univG][rand(1:2)]) : Node{AbstractOperator}(modalities[rand(1:length(modalities))])
        # root = Node{AbstractOperator}(modalities[rand(1:length(modalities))])
        left = _rand_antecedent(; modalities=modalities, minmd=minmd-1, maxmd=maxmd-1, isroot=false)
        right = _rand_antecedent(; modalities=modalities, minmd=minmd-1, maxmd=maxmd-1, isroot=false)
        leftchild!(root, left)
        rightchild!(root, right)
        parent!(left, root)
        parent!(right, root)
        return root
    end
end

function _rand_consequent(; classes)
    return classes[rand(1:length(classes))]
end

# rules = rand_rules(; numrules=3)
# TODO
# * aggiungere a parametro la profondità modale, modalità (e proposizioni?) -- fatto
# * aggiungere la possibilità di avere solo esistenziali o universali; se chiamo con rand_rules(modalities=[exL]) da errore, ad esempio
# * togliere maxnumrules perché lo vado a catturare con un vincolo quando definisco il problema?
function rand_rules(;
        modalities::Vector{AbstractUnaryIntervalOperator}=modal_relations)
    # @assert [exL,univL] ⊆ modalities "both ⟨L⟩ and [L] must be in the set of modalities"
    # @assert numrules ≥ numclasses "there must be at least $numclasses rules; currently there are $numrules rules"
    # @assert maxnumrules ≥ numrules "there must be at most $maxnumrules rules; currently there are $numrules rules (hint: set both `numrules` and `maxnumrules` keyword arguments)"
    #rules = [rand_rule(modalities=modalities, minmd=minmd, maxmd=maxmd, isroot=true) for i ∈ 1:numrules]
    numrules = rand(minr:maxr)
    rules = ClassificationRule[]
    h = rand(minh:maxh)
    for _ ∈ 1:numrules
        push!(rules, rand_rule(modalities=modalities, minmd=minmd, maxmd=rand(minmd:maxmd), isroot=true))
    end
    return ClassificationRules(rules; modalities=modalities, horizon=h)
end

function modal_depth(rule::ClassificationRule)
    return _modal_depth(antecedent(rule))
end

function _modal_depth(node::Node)
    if typeof(node.data) <: AbstractProposition
        return 0
    elseif typeof(node.data) <: AbstractUnaryIntervalOperator # Union{AbstractUniversalIntervalOperator,AbstractExistentialIntervalOperator}
        return 1 + max(_modal_depth(node.left),_modal_depth(node.right))
    end
    return 0
end

function collect_subformulas(rule::ClassificationRule)
    subformulas = Node[]
    _collect_subformulas(antecedent(rule), subformulas)
    return collect(Set(subformulas))
end

function _collect_subformulas(node::Node, subformulas)
    push!(subformulas, node)
    if isdefined(node, :left) && isdefined(node, :right)
        _collect_subformulas(node.left, subformulas)
        _collect_subformulas(node.right, subformulas)
    end
end

function collect_sorted_subformulas(rule::ClassificationRule)
    return sort!(collect_subformulas(rule), by=x -> _size(x), alg=Base.Sort.QuickSort)
end

function collect_nodes(rule::ClassificationRule)
    nodes = Node[]
    _collect_nodes(antecedent(rule), nodes)
    return nodes
end

function _collect_nodes(node::Node, nodes)
    if isdefined(node, :left)
        _collect_nodes(node.left, nodes)
    end
    push!(nodes, node)
    if isdefined(node, :right)
        _collect_nodes(node.right, nodes)
    end
end

function print_rule(io::IO, rule::ClassificationRule)
    _print_antecedent(io, antecedent(rule))
    print(io, " ⟹  ") # TODO
    _print_consequent(io, consequent(rule))
    println(io)
    return
end

function _print_antecedent(io::IO, node::Node)
    if typeof(node.data) <: AbstractUniversalIntervalOperator
        print(io, node.data)
        print(io, "(")
        _print_antecedent(io, node.left)
        print(io, implication)
        _print_antecedent(io, node.right)
        print(io, ")")
        return
    elseif typeof(node.data) <: AbstractExistentialIntervalOperator
        print(io, node.data)
        print(io, "(")
        _print_antecedent(io, node.left)
        print(io, conjunction)
        _print_antecedent(io, node.right)
        print(io, ")")
        return
    elseif typeof(node.data) <: AbstractProposition
        print(io, node.data)
        return
    elseif typeof(node.data) <: AbstractHeytingAlgebraValue
        print(io, node.data)
        return
    else
        error("you should not be here")
    end
end

function _print_consequent(io::IO, class::String)
    print(io, class)
end

function print_rules(io::IO, crules::ClassificationRules)
    for i ∈ 1:length(crules)
        print_rule(io, crules[i])
    end
end

# rules1 = rand_rules(; numrules=2)
# rules2 = rand_rules(; numrules=3)
# rules_crossover!(rules1, rules2)
function rules_crossover!(parentA::ClassificationRules, parentB::ClassificationRules, childA::ClassificationRules, childB::ClassificationRules)


    # Compute random indices i ∈ {1,…,length(parentA)} and j ∈ {1,…,length(parentB)}.
    i = rand(1:length(parentA))
    j = rand(1:length(parentB))
    # Deep copy the parent genes (i.e., classification rules).
    childA = deepcopy(parentA)
    childB = deepcopy(parentB)
    # Put the j-th rule of parentB in the i-th position of childA.
    rules(childA)[i] = rules(parentB)[j]
    # Put the i-th rule of parentA in the j-th position of childB.
    rules(childB)[j] = rules(parentA)[i]



    """
    i = rand(1:length(parentA)-1)
    j = rand(1:length(parentB)-1)
    # parentA[1:i] and parentB[j+1:end]
    rulesA = [deepcopy(parentA[1:i]); deepcopy(parentB[j+1:length(parentB)])]
    # parentA[i+1:end] and parentB[1:j]
    rulesB = [deepcopy(parentA[i+1:length(parentA)]); deepcopy(parentB[1:j])]
    childA = ClassificationRules(rulesA)
    childB = ClassificationRules(rulesB)
    """

    return
end

# TODO:
# SUM = reduce(+,[_size(x) for x ∈ sub])
# W = [1 - _size(x)/SUM for x ∈ sub]
# P = [x/(length(sub)-1) for x ∈ W]

# levels: 40% 40% 20%
#   1) rules level:
#       a) add rule, and
#       b) remove rule
#   2) rule level that replaces a subformula with another subformula
#   4) class level that changes the class
function rules_mutation!(crules::ClassificationRules)
    level = rand(1:100)

    # rules level
    if level ≤ 60
        innerlevel = rand(1:2)
        # add rule
        if innerlevel == 1
            # TODO
            #   * cancellare questo controllo perché lo farò tramite un vincolo
            #if length(crules) < maxr
            push!(rules(crules), rand_rule(; modalities=modalities(crules), minmd=minmd, maxmd=rand(minmd:maxmd), isroot=true))
            return
            #end
                # numrules(crules) = numrules(crules) + 1
            #end
        # remove rule
        elseif innerlevel == 2
            if length(crules) > minr
                deleteat!(rules(crules), rand(1:length(crules)))
                return
            end
        end
    # rule level
    elseif 60 < level && level ≤ 90
        rule = rules(crules)[rand(1:length(crules))]

        """
        subformulas = collect_sorted_subformulas(rule)
        # sum of sizes
        s = reduce(+,[_size(sf) for sf ∈ subformulas])
        # probabilities: the leaves have higher probability than any internal node, and the root has the lowest probability
        w = [1 - _size(sf)/s for sf ∈ subformulas]./(length(subformulas)-1)
        # sample a subformula with the given probabilities
        sf = sample(subformulas, Weights(w))
        """

        nodes = collect_nodes(rule)
        s = reduce(+,[_size(node) for node ∈ nodes])
        w = [1 - _size(node)/s for node ∈ nodes]./(length(nodes)-1)
        sf = StatsBase.sample(nodes, Weights(w))
        #sf = nodes[rand(1:length(nodes))]

        # upwards modal depth of the selected subformula
        umd = upwards_modal_depth(sf)
        # upwards modal depth is the previously computed one if the subformula is a leaf;
        # otherwise, it is minus 1 because we need to replace an internal node which has a modality by definition
        umd = (is_leaf(sf) ? umd : umd-1)
        # new subtree
        subtree = isdefined(sf, :parent) ?
            _rand_antecedent(modalities=modalities(crules),minmd=minmd-umd, maxmd=maxmd-umd, isroot=false) :
            _rand_antecedent(modalities=modalities(crules),minmd=minmd, maxmd=maxmd, isroot=true)
        # set parent node, if exists
        if isdefined(sf, :parent)
            subtree.parent = sf.parent
        end
    # class level
    else
        i = rand(1:length(crules))
        C = train.unique_classes
        # c = rand(1:numclasses(crules))
        c = C[rand(1:length(C))]

        # TODO: farebbe loop infinito..
        while consequent(rules(crules)[i]) == c
            c = C[rand(1:length(C))]
        end

        rules(crules)[i].consequent = c
        return
    end
end

zz(crules::ClassificationRules) = convert(Float64, length(crules)), mean_size(crules)

# julia> result = NSGAII.nsga_max(10, 20, zz, rand_rules, fCV=CV, fcross=rules_crossover!, fmut=rules_mutation!)
# julia> result = NSGAII.nsga_max(10, 50, zz, rand_rules, fCV=CV, fcross=rules_crossover!, fmut=rules_mutation!)
# julia> result = NSGAII.nsga(50, 1000, zz, rand_rules, fcross=rules_crossover!, fmut=rules_mutation!)

# end

Base.:(==)(x::Proposition, y::Proposition) = x.attribute == y.attribute && x.relation == y.relation && x.constant == y.constant
Base.:(==)(x::HeytingAlgebraValue, y::HeytingAlgebraValue) = x.value == y.value
Base.:(==)(x::AbstractUnaryIntervalOperator, y::AbstractUnaryIntervalOperator) = typeof(x) == typeof(y)
Base.:(==)(x::Node, y::Node) = _isequal(x, y)
Base.:(==)(r::ClassificationRule, s::ClassificationRule) = _isequal(antecedent(r), antecedent(s)) && consequent(r) == consequent(s)
function _isequal(x::Node, y::Node)
    if isdefined(x, :left) && isdefined(x, :right) && isdefined(y, :left) && isdefined(y, :right)
        return x.data == y.data && _isequal(x.left, y.left) && _isequal(x.right, y.right)
    else
        return x.data == y.data
    end
end

function Base.:(==)(R::ClassificationRules, S::ClassificationRules)
    if length(rules(R)) == length(rules(S))
        res = true
        SS = Set(rules(S))
        for i ∈ 1:length(rules(R))
            (!(R[i] ∈ S)) && (res = false)
        end
        # (Set(rules(R)) != Set(rules(S))) && (res = false)
        (classes(R) != classes(S)) && (res = false)
        return res
    else
        return false
    end
end

hash(p::Proposition) = hash(p.attribute) + hash(p.relation) + hash(p.constant)
hash(v::HeytingAlgebraValue) = hash(v.value)
function hash(x::Node)
    if isdefined(x, :left) && isdefined(x, :right)
        return hash(x.data) + hash(x.left) + hash(x.right)
    else
        return hash(x.data)
    end
end
hash(r::ClassificationRule) = hash(antecedent(r)) + hash(consequent(r))
hash(R::ClassificationRules) = reduce(+, [hash(r) for r ∈ rules(R)])

# vuole la proposizione che ha un attributo e la serie di quel attributo
function f̃(p::Proposition, vals::Array{Union{Missing, Float64},1})
    return sum(map(x -> p.relation(x, p.constant), vals))/length(vals)
    #return length([val for val ∈ vals[p.attribute] if p.relation(val, p.constant)])/length(vals[p.attribute])
end

function f̃(p::Proposition, vals::Vector{Float64})
    return sum(map(x -> p.relation(x, p.constant), vals))/length(vals)
    #return length([val for val ∈ vals[p.attribute] if p.relation(val, p.constant)])/length(vals[p.attribute])
end

function ==̃ₕ(x::Int64, y::Int64; horizon::Int64=4)
    # |x-y| ≥ h
    if abs(x-y) ≥ horizon
        return float(0)
    # |x-y| < h
    else
        return (horizon - abs(x-y))/horizon
    end
end

function <̃ₕ(x::Int64, y::Int64; horizon::Int64=4)
    # y-x > h
    if y-x > horizon
        return float(1)
    # y ≦ x
    elseif y ≤ x
        return float(0)
    # y-x ≦ h
    elseif y-x ≤ horizon
        return (y-x)/horizon
    end
end

# ℋ = ⟨H, ∩, ∪, ↪, 0, 1⟩
∩(args...) = min(args...)
∪(args...) = max(args...)
(↪)(a::Float64, b::Float64) = a ≤ b ? float(1) : b

# w₁ = (x, y) and w₂ = (w, z)
function R̃ₓ(rel::AbstractUnaryIntervalOperator, w₁::NTuple{2, Int64}, w₂::NTuple{2, Int64}; horizon::Int64=4)
    x, y, w, z = w₁[1], w₁[2], w₂[1], w₂[2]
    # relations
    if typeof(rel) <: Union{ExistentialLater, UniversalLater}
        return <̃ₕ(y, w; horizon)
    elseif typeof(rel) <: Union{ExistentialAfter, UniversalAfter}
        return ==̃ₕ(y, w; horizon)
    elseif typeof(rel) <: Union{ExistentialOverlaps, UniversalOverlaps}
        return ∩(<̃ₕ(x, w; horizon), <̃ₕ(w, y; horizon), <̃ₕ(y, z; horizon))
    elseif typeof(rel) <: Union{ExistentialEnds, UniversalEnds}
        return ∩(<̃ₕ(x, w; horizon), ==̃ₕ(y, z; horizon))
    elseif typeof(rel) <: Union{ExistentialDuring, UniversalDuring}
        return ∩(<̃ₕ(x, w; horizon), <̃ₕ(z, y; horizon))
    elseif typeof(rel) <: Union{ExistentialBegins, UniversalBegins}
        return ∩(==̃ₕ(x, w; horizon), <̃ₕ(z, y; horizon))
    # inverse relations
    elseif typeof(rel) <: Union{ExistentialInverseLater, UniversalInverseLater}
        return <̃ₕ(z, x; horizon)
    elseif typeof(rel) <: Union{ExistentialInverseAfter, UniversalInverseAfter}
        return ==̃ₕ(z, x; horizon)
    elseif typeof(rel) <: Union{ExistentialInverseOverlaps, UniversalInverseOverlaps}
        return ∩(<̃ₕ(w, x; horizon), <̃ₕ(x, z; horizon), <̃ₕ(z, y; horizon))
    elseif typeof(rel) <: Union{ExistentialInverseEnds, UniversalInverseEnds}
        return ∩(<̃ₕ(w, x; horizon), ==̃ₕ(y, z; horizon))
    elseif typeof(rel) <: Union{ExistentialInverseDuring, UniversalInverseDuring}
        return ∩(<̃ₕ(w, x; horizon), <̃ₕ(y, z; horizon))
    elseif typeof(rel) <: Union{ExistentialInverseBegins, UniversalInverseBegins}
        return ∩(==̃ₕ(x, w; horizon), <̃ₕ(y, z; horizon))
    end
end

hash(x::Tuple{Node, Tuple{Int64, Int64}}) = hash(x[1]) + hash(x[2][1]) + hash(x[2][2])

global H = Dict{Node, Dict{ModalInstance, Dict{NTuple{2,Int}, Float64}}}()

function mc(I::ModalInstance, φ::ClassificationRule; horizon::Int=4, memo::Bool=false)
    if !memo
        return _mc(I,φ; horizon)
    else 
        return _mc(I,φ,H; horizon)
    end
end

function _mc(I::ModalInstance, φ::ClassificationRule, H::Dict{Node, Dict{ModalInstance, Dict{NTuple{2,Int}, Float64}}}; horizon::Int=4)

    if haskey(H, antecedent(φ)) && haskey(H[antecedent(φ)], I)
        return H[antecedent(φ)][I][(1,2)]
    end

    sub = collect_sorted_subformulas(φ)
    # L = Dict{Tuple{Union{Node{AbstractOperator}, Node{AbstractProposition}, Node{AbstractHeytingAlgebraValue}},NTuple{2,Int64}},Float64}()
    N = length(I, 1)  # TODO works only on the first (1) frame, i.e., time series
    intervals = [(x,y) for x in 1:N for y in 1:N if x < y]
    for ψ ∈ sub

        if haskey(H, ψ) && haskey(H[ψ],I)
            continue
        end

        if !haskey(H, ψ)
            H[ψ] = Dict{ModalInstance, Dict{NTuple{2,Int}, Float64}}()
            H[ψ][I] = Dict{NTuple{2,Int}, Float64}()
        else
            if !haskey(H[ψ], I)
                H[ψ][I] = Dict{NTuple{2,Int}, Float64}()
            end
        end

        if typeof(ψ.data) <: AbstractProposition
            if length(H[ψ][I]) == 0
                for (x, y) ∈ intervals 
                    # if !haskey(H[ψ][I], (x, y))
                    H[ψ][I][(x,y)] = f̃(ψ.data, I[1][ψ.data.attribute][x:y]) # TODO works only for time series
                    # end
                end
            end

        elseif typeof(ψ.data) <: AbstractHeytingAlgebraValue
            if length(H[ψ][I]) == 0
                for (x, y) ∈ intervals
                    H[ψ][I][(x,y)] = ψ.data.value
                end
            end

        elseif typeof(ψ.data) == ExistentialGlobally
            S = Float64[]
            for op in [exL,exA,exO,exE,exD,exB,exL̄,exĀ,exŌ,exĒ,exD̄,exB̄]
                for (x, y) ∈ intervals
                    s = 0.0
                    for (w, z) ∈ intervals
                        # if !haskey(H[ψ][I], (w, z))
                        s = ∪(s, ∩(R̃ₓ(op, (x, y), (w, z)), ∩(H[ψ.left][I][(w,z)], H[ψ.right][I][(w,z)])))
                        # end
                    end
                    push!(S, s)
                end
            end
            H[ψ][I][(1,2)] = ∪(S...)
            return H[ψ][I][(1,2)]

        elseif typeof(ψ.data) == UniversalGlobally
            S = Float64[]
            for op in [univL,univA,univO,univE,univD,univB,univL̄,univĀ,univŌ,univĒ,univD̄,univB̄]
                for (x, y) ∈ intervals
                    s = 1.0
                    for (w, z) ∈ intervals
                        # if !haskey(H[ψ][I], (w, z))
                        s = ∩(s, ↪(R̃ₓ(op, (x, y), (w, z)), ↪(H[ψ.left][I][(w,z)], H[ψ.right][I][(w,z)])))
                        # end
                    end
                    push!(S, s)
                end
            end
            H[ψ][I][(1,2)] = ∩(S...)
            return H[ψ][I][(1,2)]

        elseif typeof(ψ.data) <: AbstractExistentialIntervalOperator
            if length(H[ψ][I]) == 0
                for (x, y) ∈ intervals
                    # if !haskey(H[ψ][I], (x,y))
                    s = 0.0
                    for (w, z) ∈ intervals
                        # if !haskey(H[ψ][I], (w, z))
                        s = ∪(s, ∩(R̃ₓ(ψ.data, (x, y), (w, z)), ∩(H[ψ.left][I][(w,z)], H[ψ.right][I][(w,z)])))
                        # end
                    end
                    H[ψ][I][(x,y)] = s
                    # end
                end
            end

        elseif typeof(ψ.data) <: AbstractUniversalIntervalOperator
            if length(H[ψ][I]) == 0
                for (x, y) ∈ intervals
                    # if !haskey(H[ψ][I], (x,y))
                    s = 1.0
                    for (w, z) ∈ intervals
                        # if !haskey(H[ψ][I], (w, z))
                        s = ∩(s, ↪(R̃ₓ(ψ.data, (x, y), (w, z)), ↪(H[ψ.left][I][(w,z)], H[ψ.right][I][(w,z)])))
                        # end
                    end
                    H[ψ][I][(x,y)] = s
                    # end
                end
            end
        end
    end
    # return 0.0
end

function _mc(I::ModalInstance, φ::ClassificationRule; horizon::Int=4)
    sub = collect_sorted_subformulas(φ)
    L = Dict{Tuple{Union{Node{AbstractOperator}, Node{AbstractProposition}, Node{AbstractHeytingAlgebraValue}},NTuple{2,Int64}},Float64}()
    N = length(I, 1)  # TODO works only on the first (1) frame, i.e., time series
    intervals = [(x,y) for x in 1:N for y in 1:N if x < y]
    for ψ ∈ sub
        if typeof(ψ.data) <: AbstractProposition
            for (x, y) ∈ intervals # setdiff(𝕀, [(-1, 0)])
                # L[(ψ, (x, y))] = ((x, y) != (-1, 0) ? f̃(ψ.data, ts[ψ.data.attribute][x:y]) : 0.0) # TODO is (-1, 0) needed for a proposition? Perhaps not
                L[(ψ, (x ,y))] = f̃(ψ.data, I[1][ψ.data.attribute][x:y]) # TODO works only for time series
            end
        elseif typeof(ψ.data) <: AbstractHeytingAlgebraValue
            for (x, y) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                L[(ψ, (x, y))] = ψ.data.value
            end
        elseif typeof(ψ.data) == ExistentialGlobally
            S = Float64[]
            for op in [exL,exA,exO,exE,exD,exB,exL̄,exĀ,exŌ,exĒ,exD̄,exB̄]
                for (x, y) ∈ intervals
                    s = 0.0
                    for (w, z) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                        s = ∪(s, ∩(R̃ₓ(op, (x, y), (w, z)), ∩(L[(ψ.left, (w, z))], L[(ψ.right, (w, z))])))
                    end
                    push!(S, s)
                end
            end
            return ∪(S...)
        elseif typeof(ψ.data) == UniversalGlobally
            S = Float64[]
            for op in [univL,univA,univO,univE,univD,univB,univL̄,univĀ,univŌ,univĒ,univD̄,univB̄]
                for (x, y) ∈ intervals
                    s = 1.0
                    for (w, z) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                        s = ∩(s, ↪(R̃ₓ(op, (x, y), (w, z)), ↪(L[(ψ.left, (w, z))], L[(ψ.right, (w, z))])))
                    end
                    push!(S, s)
                end
            end
            return ∩(S...)
        elseif typeof(ψ.data) <: AbstractExistentialIntervalOperator
            for (x, y) ∈ intervals # (typeof(ψ.data) === ExistentialLater ? 𝕀 : setdiff(𝕀,[(-1, 0)]))
                s = 0.0
                for (w, z) ∈ intervals # setdiff(𝕀,[(-1, 0)])
                    s = ∪(s, ∩(R̃ₓ(ψ.data, (x, y), (w, z)), ∩(L[(ψ.left, (w, z))], L[(ψ.right, (w, z))])))
                end
                L[(ψ, (x, y))] = s
            end
        elseif typeof(ψ.data) <: AbstractUniversalIntervalOperator
            for (x, y) ∈ intervals # (typeof(ψ.data) === UniversalLater ? 𝕀 : setdiff(𝕀,[(-1,0)]))
                s = 1.0
                for (w, z) ∈ intervals # setdiff(𝕀,[(-1,0)])
                    s = ∩(s, ↪(R̃ₓ(ψ.data, (x, y), (w, z)), ↪(L[(ψ.left, (w, z))], L[(ψ.right, (w, z))])))
                end
                L[(ψ, (x, y))] = s
            end
        end
    end
end

# modal_relations = [exL,exA,exO,exE,exD,exB,exL̄,exĀ,exŌ,exĒ,exD̄,exB̄,univL,univA,univO,univE,univD,univB,univL̄,univĀ,univŌ,univĒ,univD̄,univB̄]

function mmmc(ℐ::ClassificationDataset, Γ::ClassificationRules)
    accuracy = 0.0

    # φ ↦ M ↦ (x, y) ↦ ℝ
    # H = Dict{Node, Dict{ModalInstance, Dict{NTuple{2,Int}, Float64}}}()

    for i in 1:length(ℐ.instances)
        I = ℐ.instances[i]
        class = ℐ.classes[i]
        L = Dict{ClassificationRule, Float64}()

        for ρ in Γ
            L[ρ] = mc(I, ρ; horizon=horizon(Γ), memo=true)
        end

        maxv = argmax(L)
        maxv_rules = [k for (k,v) ∈ L if v == L[maxv]]
        pred = consequent(maxv_rules[rand(1:length(maxv_rules))])

        if pred == class
            accuracy = accuracy + 1.0
        end
    end
    return accuracy/length(ℐ)
end


# function mmmc(𝓣::DataFrame, Γ::ClassificationRules)
#     accuracy = 0.0
#     for T ∈ eachrow(𝓣)
#         L = Dict{ClassificationRule, Float64}()
#         for ρ ∈ Γ
#             L[ρ] = mc(values(T), ρ)
#         end

#         # The key (i.e., rule) with the maximum mapped value.
#         maxv = argmax(L)
#         # The keys (i.e., rules) with the maximum mapped value.
#         maxv_rules = [k for (k,v) ∈ L if v == L[maxv]]
#         # Tights (i.e., more than one rule with the maximum mapped value) are handled randomly.
#         pred = consequent(maxv_rules[rand(1:length(maxv_rules))])

#         if pred == class(T)
#             accuracy = accuracy + 1.0
#         end
#     end
#     return accuracy/nrow(𝓣)
# end

# z(crules::ClassificationRules) = mmmc(instances(td), crules), -1 * mean_size(crules)

function z(crules::ClassificationRules)
    # sizes = [size(rule) for rule ∈ rules(crules)]
    # # max_size = max(sizes...)
    # max_size = 2^(maxdepth+1)-1
    round(mmmc(train, crules), digits=4), round(1/complexity(crules), digits=4) # 1 - sum(sizes)/(max_size*maxr) #1 - sum(sizes)/(max_size*100) # (mean_size(crules)/max_size) * (length(classes(td))/length(rules(crules))) # 1.0 - sum(sizes)/(max_size*12)
end

"""
function CV(crules::ClassificationRules)
    # return length(crules) ≥ numclasses(crules) ? 0 : length(crules)
    return length(classes(td)) - length(classes(crules)) # 1.0 #length(classes(crules)) == length(classes(td)) ? 0 : 1 #length(classes(td)) - length(classes(crules))
end
"""
function CV(crules::ClassificationRules)
    # return length(crules) ≥ numclasses(crules) ? 0 : length(crules)
    # return length(classes(crules)) == length(classes(td)) ? 0.0 : length(classes(td)) - length(classes(crules))
    return abs(length(train.unique_classes) - length(classes(crules)))
end

function init()
    # numrules = rand(minr:maxr)
    # maxmd = rand(2:maxdepth)
    # rand_rules(minmd=2, maxmd=maxdepth)
    rand_rules()
end

function hypervolume_indicator(y::Vector{Tuple{Float64,Float64}}; refpoint=(0,0))
    sort!(y, by=x->x[1], rev=true, alg=QuickSort)

    hi = 0.0
    # (index, tuple)
    for (i, t) in enumerate(y)
        # if it is the first, then
        if i == 1
            hi += abs(t[1] - refpoint[1]) * abs(t[2] - refpoint[2])
        else
            s = y[i-1]
            hi += abs(t[1] - refpoint[1]) * abs(t[2] - s[2])
        end
    end

    return round(hi, digits=4)
end


# init() = rand_rules(numrules=4, maxmd=3)

@time M = nsga_max(6,5,z,init,fcross=rules_crossover!,fmut=rules_mutation!,pmut=0.3)