## Software:
## (c) Vicenç Torra
## Current version: 20241118
##
## Software for non-additive measures and integrals
##    (these measures are also known as fuzzy measures, monotonic games,
##    and include belief functions, k-additive measures,
##    distorted probabilities) 
## Software includes Choquet and Sugeno integrals, 
##    measure identification from examples,
##    Shapley values, Upsilon values (here called tau),
##    and other indices (interactions), transforms, etc.
## 
## 
## This software is similar and extends previous software in R
##        http://www.mdai.cat/ifao/prog.ifao.models.r
## Reference on this module in python (functions about model learning):
##   E. Turkarslan, V. Torra (2022) Measure Identification for the
##        Choquet integral: a Python module, Int. J. of
##        Comp. Intel. Systems 15:89 (2022) 
##        https://doi.org/10.1007/s44196-022-00146-w
##
## References on aggregation functions: 
##     V. Torra, Y. Narukawa (2007) Modeling Decisions, Springer.
## References on the Upsilon values:
##     V. Torra (2024) $\Upsilon$-Values: Power Indices à La Orness for
##          Nonadditive Measures. IEEE Trans. Fuzzy Syst. 32:7 4099-4108
##          https://doi.org/10.1109/TFUZZ.2024.3392268
## References on the (max,+)-transform:
##     V. Torra, (Max,\oplus)-transforms and genetic algorithms for
##          fuzzy measure identification,
##          Fuzzy Sets and Systems 451 (2022) 253-265
##          https://doi.org/10.1016/j.fss.2022.09.008
## References on the WOWA: 
##     V. Torra, The Weighted OWA operator, Int. J. of Intel. Systems,
##          12 (1997) 153-166.
##          http://dx.doi.org/10.1002/(SICI)1098-111X(199702)12:2<153::AID-INT3>3.0.CO;2-P
##   Original interpolation: 
##     V. Torra, The WOWA operator and the interpolation function W*: Chen and
##          Otto's interpolation method revisited, Fuzzy Sets and Systems, 113:3
##          (2000) 389-396. http://dx.doi.org/10.1016/S0165-0114(98)00040-2
##   Discussion on intepolation functions
##     V. Torra, Z. Lv, On the WOWA operator and its interpolation function.
##          Int. J. Intell. Syst. 2:(1 (2009): 1039-105
##          http://dx.doi.org/10.1002/int.20371
##   An overview of the WOWA and related results: 
##     V. Torra, The WOWA Operator: A Review, in R. R. Yager, J. Kacprzyk,
##          G. Beliakov (eds.) Recent Developments in the Ordered Weighted
##          Averaging Operators, Springer 2011 17-28.
##          https://doi.org/10.1007/978-3-642-17910-5_2
##

import random
import numpy as np
import numpy.linalg as lina
import matplotlib.pyplot as plt
import statistics
import math
import sklearn ## -- DEPRECATED
## import scikit-learn
import time

## Necessary to solve the quadratic problem with constraints: 
##      $ pip install cvxopt
##      from cvxopt import matrix, solvers
##      import cvxopt.matrix as co.matrix

import cvxopt as opt

## exec(open("prog.vectors.matrices.web.txt").read())


## -----------------------------------------------------
## Aggregation functions: AM, WM, OWA, QWM, WOWA
## -----------------------------------------------------

def ciAM (data):
    return sum(data)/len(data)

def ciWM (data, w):
    """
    Function: Weighted mean
    Example: 
    ciWM([1,2,3,4],[0.1,0.2,0.3,0.4])
    ciWM([1,2,3,4],[0.25,0.25,0.25,0.25])
    """
    return sum(map(lambda x,y:x*y, data, w))

def ciHM (data, w):
    return 1/sum(map(lambda x,w:w/x, data,w))

def ciOWA (data, p):
    """
    Function: OWA operator 
    Example: 
    ciOWA([1,3,2,4],[0.1,0.2,0.3,0.4])      ##   2.0
    ciOWA([1,2,3,4],[0.25,0.25,0.25,0.25])  ##   2.5
    ciOWA([1,4,3,2],[0,1,0,0])              ##   3     -- 2nd largest
    """
    return sum(map(lambda x,y:x*y, sorted(data, reverse=True), p))

def ciQWM (data, w, f, fm1):
    """
    Function: Weighted mean
    Example: 
    ciQWM([1,2,3,4],[0.1,0.2,0.3,0.4], lambda x:x, lambda x:x)
    ciQWM([1,2,3,4],[0.25,0.25,0.25,0.25], lambda x:x, lambda x:x)
    ciQWM([1,2,3,4],[0.1,0.2,0.3,0.4], lambda x:x*x, lambda x:math.sqrt(x))
    ciQWM([1,2,3,4],[0.25,0.25,0.25,0.25], lambda x:x*x, lambda x:math.sqrt(x))
    """
    return fm1(sum(map(lambda x,y:f(x)*y, data, w)))

def ciPartialSum (v):
    """
    Function:
    Given a vector (x_1, x_2, ..., x_i, ..., x_n)
    returns the vector (x_1, x_1+x_2, ..., sum_{j<=i} x_j, ..., sum x_j)
    Input:   v: a vector (x_1, ..., x_n)
    Output:  a vector of partial sums
    Example: ciPartialSum([1,2,3,4])
    """
    acc = 0
    cpv = []
    for i in range(0,len(v)):
        acc = acc + v[i]
        cpv.append(acc)
    return(cpv)


def ciEvalLinearQ (w, x):
    """
    Function:
    #   A function for linear interpolation and evaluation at a given point
    #   interpolates points {(0,0)} cup {(i/n, sum_{j<=i} x_j)}_{i=1, ..., n}
    #   This function is used to compute the function w* in WOWA 
    Input:
    #   w: the vector (w_1, ..., w_n)
    #        in WOWA w_i should be positive and add to one
    #   x: a point where to evaluate the interpolation function
    #        in WOWA x in [0,1]
    Output: the evaluation of the interpolated polynomial at x
    Example: ciEvalLinearQ([0.1,0.2,0.3,0.4],0.2)
    """
    if (x<=0):
        return 0
    elif (x>=1):
        return 1
    else:
        n = len(w)
        accw = ciPartialSum(w)
        i = 0
        while (x>i/n):
            i=i+1
        i = i-1
        if i==0:
            return x*(accw[0])/(1/n)
        else:
            return (accw[i-1] + (x-i/n)*(accw[i]-accw[i-1])/(1/n))

def ciLinearQ (w):
    """
    Example: 
    #    ciLinearQ([0.1,0.2,0.3,0.4])(0.5)
    #    list(map(ciLinearQ([0.2,0.3,0.4,0.1]),[0.1,0.2,0.3,0.4,0.5]))
    """
    return lambda x: ciEvalLinearQ(w,x)



def ciWOWA (data, owaWeightsVect, wmWeightsVect):
    """
    Function: 
    #   The WOWA operator (with w* defined by means of a linear interpolation)
    Input:
    #   data: vector of data elements (a_1, ..., a_n)
    #   owaWeightsVect: vector of OWA weights (w_1, ..., w_n)
    #   wmWeightsVect: vector of weights (p_1, ..., p_n)
    Output:
    #   sum a_{sigma(i)}*omega_i
    #     (where sigma is a permutation such that a_{sigma(i-1)}>=a_{sigma(i)}
    #     omega_i=w*(sum_{j leq i}p_{sigma(j)}) - w*(sum_{j<i}p_{sigma(j)})
    #     w* interpolates {(i/n),sum_{j leq i} w_j}
    Notes:
    #    w*: is wStar in this function
    #    omega: is weights in this function
    Examples: 
    ciWOWA([4,3,2,1],[0.1,0.2,0.3,0.4],[0.25,0.25,0.25,0.25]) # 2.0
    ciWOWA([1,4,3,2],[0.1,0.2,0.3,0.4],[0.25,0.25,0.25,0.25]) # 2.0
    ciWOWA([4,3,2,1],[0.25,0.25,0.25,0.25],[0.1,0.2,0.3,0.4]) # 2.0
    ciWOWA([1,4,3,2],[0.25,0.25,0.25,0.25],[0.1,0.2,0.3,0.4]) # 2.6
    """
    wStar = ciLinearQ(owaWeightsVect)
    sortedData, sortedWeights = sortWithOneListDecreasing(data, wmWeightsVect)
    accWeights = list(map(wStar,ciPartialSum(sortedWeights)))
    weights = [accWeights[0]] # accWeights
    for i in range(1,len(accWeights)):
        weights.append(accWeights[i]-accWeights[i-1])
    return vectorProduct(sortedData,weights)


def ciWOWAq (data, wStar, wmWeightsVect): 
    """
    Function: 
    #   The WOWA operator using a fuzzy quantifier w*
    Input:
    #   data: vector of data elements (a_1, ..., a_n)
    #   wStar: Fuzzy quantifier, a function w*:[0,1]->[0,1]
    #   wmWeightsVect: vector of weights (p_1, ..., p_n)
    Output:
    #   sum a_{sigma(i)}*omega_i
    #     (where sigma is a permutation such that a_{sigma(i-1)}>=a_{sigma(i)}
    #     omega_i=w*(sum_{j leq i}p_{sigma(j)}) - w*(sum_{j<i}p_{sigma(j)})
    Notes:
    #    w*: is wStar in this function
    #    omega: is weights in this function
    Examples: 
    #    ciWOWAq([0.5,0.25,0.8,0.0,0.75],lambda x: x**(1/4),[0.1,0.1,0.6,0.1,0.1]) # 0.7526147305039615
    """
    #wStar = ciLinearQ(owaWeightsVect)
    sortedData, sortedWeights = sortWithOneListDecreasing(data, wmWeightsVect)
    accWeights = list(map(wStar,ciPartialSum(sortedWeights)))
    weights = [accWeights[0]] # accWeights
    for i in range(1,len(accWeights)):
        weights.append(accWeights[i]-accWeights[i-1])
    return vectorProduct(sortedData,weights)

## Example 6.15 @ V. Torra, Y. Narukawa (2007) Modeling decisions, Springer.
## Examples of weights and fuzzy quantifiers
ciP1 = [0.1,0.1,0.6,0.1,0.1]
ciP2 = [0.1,0.1,0.1,0.6,0.1]
ciQQ1 = lambda x: x**(1/4)
ciQQ3 = lambda x: x**4
ciQQ4 = lambda x: x

def ciQQ2(x):
    if x<0.2:
        return 0
    elif (x <= 4/5):
        return((5/3)*(x-1/5))
    return(1)

## Example 6.29 @ V. Torra, Y. Narukawa (2007) Modeling decisions, Springer. 
def ciQ629(x):
    return 1 if x>=0.75 else x/0.75

## This is equivalent to the OWA_q
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ1,[0.2,0.2,0.2,0.2,0.2]) # 0.6887180338897262
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ2,[0.2,0.2,0.2,0.2,0.2]) # 0.5000000000000001
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ3,[0.2,0.2,0.2,0.2,0.2]) # 0.14128000000000004
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ4,[0.2,0.2,0.2,0.2,0.2]) # 0.4600000000000001 
#
## Example 6.16 @ V. Torra, Y. Narukawa (2007) Modeling decisions, Springer. 
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ1,ciP1) # 0.7526147305039615
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ2,ciP1) # 0.5
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ1,ciP2) # 0.579140525818457
# ciWOWAq([0.5,0.25,0.8,0.0,0.75],ciQQ2,ciP2) # 0.7083333333333333 



## -----------------------------------------------------------------
## Choquet integral and Sugeno integral
## computation from fuzzy measure and computation from Mobius transform
## -----------------------------------------------------------------

def choquetIntegral (f, fm):
    """
    Function: compute the Choquet integral of the function f w.r.t fm
    Input: 
      f: function to be integrated
      fm: fuzzy measure
      *NOTE*: order of the elements in f, same in the measure
    Output: 
      the Choquet integral of f with respect to fm
    Examples:
    choquetIntegral([0,0.5,1,1],fmW)        #// Output 0.4666
    choquetIntegral([0.5,0.5,0.5,0.5],fmW)  #// Output 0.5
    choquetIntegral([1,0.5,0,0.5],fmW)      #// Output 0.8333
    choquetIntegral([1,0.5,0.5,0.5],fmW)    #// Output 0.8333
    choquetIntegral([0.5,1,1,0.5],fmW)      #// Output 0.8
    choquetIntegral([1,1,0.5,1],fmW)        #// Output 1.0
    """
    items = list(map(lambda i: 2**(len(f)-1-i), range(0,len(f)))) # 2^(n-1), .. 4,2,1
    (fSorted, itemsSorted) = sortWithOneListIncreasing (f.copy(), items)     # f_{s(i)}
    total = fSorted[0]
    fmOfSetItems = 2**(len(f))-1
    i = 1
    while (i<len(f)):
        fmOfSetItems = fmOfSetItems - itemsSorted[i-1]
        total = total + (fSorted[i]-fSorted[i-1])*fm[fmOfSetItems]
        i = i+1
    return(total)

def choquetIntegralFMasFunction (f, fm):
    """
    Function: compute the Choquet integral of the function f w.r.t fm
    Input: 
      f: function to be integrated
      fm: fuzzy measure (a function [0,0,1,0,1] => fm([0,0,1,0,1]))
      *NOTE*: order of the elements in f, same in the measure
    Output: 
      the Choquet integral of f with respect to fm
    Examples:
    choquetIntegralFMasFunction([0,0.5,1,1],lambda s:fmW[set2Int(s)])        #// Output 0.4666
    choquetIntegralFMasFunction([0.5,0.5,0.5,0.5],lambda s:fmW[set2Int(s)])  #// Output 0.5
    choquetIntegralFMasFunction([1,0.5,0,0.5],lambda s:fmW[set2Int(s)])      #// Output 0.8333
    choquetIntegralFMasFunction([1,0.5,0.5,0.5],lambda s:fmW[set2Int(s)])    #// Output 0.8333
    choquetIntegralFMasFunction([0.5,1,1,0.5],lambda s:fmW[set2Int(s)])      #// Output 0.8
    choquetIntegralFMasFunction([1,1,0.5,1],lambda s:fmW[set2Int(s)])        #// Output 1.0
    """
    items = list(range(0,len(f)))           #  X={0,...,N-1}
    (fSorted, itemsSorted) = sortWithOneListIncreasing (f.copy(), items)     # f_{s(i)}
    setItems = [1]*len(f)                   # A_{s(1)} = all elements
    # fmOfSetItems = set2Int(setItems)
    total = fSorted[0] * fm(setItems)
    i = 1
    while (i<len(f)):
        setItems[itemsSorted[i-1]]=0
        total = total + (fSorted[i]-fSorted[i-1])*fm(setItems)
        i = i+1
    return(total)


# product(coefCIMoebius(values), fromFM2Moebius(fmW,4))
def ciFromMoebius (values, moebius):
    """
    Function: Computes the CI from the Moebius transform
    Examples: 
    ciFromMoebius([0,0.5,1,1],fromFM2Moebius(fmW,4))       # Output 0.4666
    cifromMoebius([0.5,0.5,0.5,0.5],fromFM2Moebius(fmW,4)) # Output 0.5
    ciFromMoebius([1,0.5,0,0.5],fromFM2Moebius(fmW,4))     # Output 0.8333
    ciFromMoebius([1,0.5,0.5,0.5],fromFM2Moebius(fmW,4))   # Output 0.8333
    ciFromMoebius([0.5,1,1,0.5],fromFM2Moebius(fmW,4))     # Output 0.8
    ciFromMoebius([1,1,0.5,1],fromFM2Moebius(fmW,4))       # Output 1.0
    """
    return(vectorProduct(coefCIMoebius(values), moebius))

def ciSugenoIntegral (f, fm):
    """
    Function: compute the Sugeno integral of the function f w.r.t fm
    Input: 
      f: function to be integrated
      fm: fuzzy measure
      *NOTE*: order of the elements in f, same in the measure
    Output: 
      the Sugeno integral of f with respect to fm
    Examples:
    ciSugenoIntegral([0,0.5,1,1],fmW)        #// Output 0.5
    ciSugenoIntegral([0.5,0.5,0.5,0.5],fmW)  #// Output 0.5
    ciSugenoIntegral([1,0.5,0,0.5],fmW)      #// Output 0.666
    ciSugenoIntegral([1,0.5,0.5,0.5],fmW)    #// Output 0.666
    ciSugenoIntegral([0.5,1,1,0.5],fmW)      #// Output 0.666
    ciSugenoIntegral([1,1,0.5,1],fmW)        #// Output 1.0
    """
    items = list(range(0,len(f)))           #  X={0,...,N-1}
    (fSorted, itemsSorted) = sortWithOneListIncreasing (f.copy(), items)     # f_{s(i)}
    setItems = [1]*len(f)                   # A_{s(1)} = all elements
    fmOfSetItems = set2Int(setItems)
    total = min(fSorted[0], fm[fmOfSetItems])
    i = 1
    while (i<len(f)):
        setItems[itemsSorted[i-1]]=0   ## essa(i)*As(i)
        total = max(total, min(fSorted[i],fm[set2Int(setItems)]))
        i = i+1
    return(total)


def coefCIMoebius(values):
    """
    Function: Given a vector of values, it builds the vector so that 
    when multiplied by the Moebius transform it results into the Choquet integral. 
    Observe that the length of the vector is 2^len(values)
    Example: 
    coefCIMoebius([1,2,3,4])
    len(coefCIMoebius([1,2,3,4]))==16
    """
    n = len(values)
    positions = list(range(0,n))
    (sortedVal, sortedPos) = sortWithOneListIncreasing (values, positions)
    sortedLowHigh = [0]+sortedVal
    mySet = firstSet(n)
    set2IntS = 0
    aOfMySet = [0]*(2**n)
    while notLastSet(mySet):
        mySet = nextSet(mySet, n)
        set2IntS = set2IntS+1
        Asi = [1]*n
        for i in range(1,n+1):
            diff=sortedLowHigh[i]-sortedLowHigh[i-1]
            if (subseteq(mySet, Asi, n)):
                aOfMySet[set2IntS] = aOfMySet[set2IntS] + diff
            Asi[sortedPos[i-1]] = 0
    return(aOfMySet)

def allCoefCIMoebius (matrixDataWithOutput):
    """
    Function: Given a matrix (m pairs (input, output) = [x_1, dots, x_n, output])
    builds the coefficients for each input. So, the output is of size m * 2^n
    Example: 
    allCoefCIMoebius(data85)
    """
    n = len(matrixDataWithOutput[0])-1
    matrixCoef = list(map(lambda row: coefCIMoebius(row[:n:]), matrixDataWithOutput))
    vectorOut  = list(map(lambda row: row[n], matrixDataWithOutput))
    return(matrixCoef, vectorOut)

## ------------------------------------------------------------
## Determine a fuzzy measure from examples.
## Main function:
##      ciSolveMSE (matrixDataWithOutput, belief=False, kAdditive=False, shapleyValue=False)
##      ciFindMoebius (matrixDataWithOutput, belief=False, kAdditive=False, shapleyValue=False)
## Also:
##      ciFindOWAweights (matrixDataWithOutput)
##      ciFindMoebiusOWA (matrixDataWithOutput)
##      def ciFindWMweights (matrixDataWithOutput)    ## Return weights (additive weights)
## ------------------------------------------------------------


# Equations about the monotonicity of the Moebius transform.
# Inequality constraints
#   OBS! They should be of the form Gx <= h  (i.e., less than or equal to zero!)
def buildEquationsMoebius (n):
    """
    Function: Builds all inequality equations for Moebius when we have n input data
    Examples:
    buildEquationsMoebius(2)
    buildEquationsMoebius(3)
    """
    allEquations = []
    mySet = firstSet(n)
    set2IntS = 0
    while notLastSet(mySet):
        mySet = nextSet(mySet, n)
        set2IntS = set2IntS+1
        if nElemsSet(mySet,n)==1: 
            equation = [0.0]*(2**n)
            equation[set2IntS]=-1.0
            allEquations.append(equation)            
        else:  
            indices = [x[0] for x in zip(range(0,len(mySet)),mySet) if x[1]==1]
            nIndices = len(indices)
            EQs = [[0.0]*(2**n) for i in range(nIndices)] ## [[0.0]*(2**n)]*nIndices
            mySetB = firstSet(nIndices)
            set2IntB = 0
            while notLastSet(mySetB):
                mySetB = nextSet(mySetB, nIndices)
                set2IntB = set2IntB + 1
                for i in range(0,nIndices):
                    if mySetB[i]==1:
                        selectIndices = [x[0] for x in zip(indices, mySetB) if x[1]==1]
                        posMySetB = setIndex2Int(selectIndices,n)
                        EQs[i][posMySetB]=-1.0
            allEquations = allEquations + EQs
    allEquationsNoEmptySet = list(map(lambda row:row[1:], allEquations))
    return(allEquationsNoEmptySet, [0.0]*len(allEquations))


def buildOFMoebius (matrixDataWithOutput):
    """
    Function:
    Example: 
    buildOFMoebius (data85)
    """
    myMatrix, vectorOutm = allCoefCIMoebius (matrixDataWithOutput)
    matrixCoefA = list(map(lambda row:row[1:], myMatrix))  ## First column in matrixCoefA is zero
    Q = matrixProduct(transpose(matrixCoefA),matrixCoefA)
    Q = matrixProductConstant(Q,2)
    q = matrixProduct(transpose(matrixCoefA),list(map(lambda x:[x],vectorOutm)))
    q = matrixProductConstant(q, -2)
    q = list(map(lambda x:x[0],q))
    return(Q,q)

def buildEquationsMoebiusPositive (n):
    positiveEQ = [[0.0]*n for i in range(n)]
    for i in range(0,n):
        positiveEQ[i][i]=-1.0
    return (positiveEQ, [0.0]*n)

def buildEquationsKAdditive(n, kAdditive):
    """
    Function: Equations for k-additive for given k
    Moebius transform position corresponding to emptyset removed. 
    buildEquationsKAdditive(3, 1)
    buildEquationsKAdditive(3, 2)
    """
    alleqs=[]
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        if nElemsSet(mySetA, n)>kAdditive:
            thisEq = [0.0]*(2**n)
            thisEq[set2IntA]=-1.0
            alleqs.append(thisEq)
    nalleqs = list(map(lambda x:x[1:], alleqs))
    # print("additive:")
    # print(nalleqs)
    return(nalleqs, [0.0]*len(alleqs))

def buildEquationsShapley (n, shapley):
    """ Function: 
    Build equations for the Shapley values. 
    Position corresponding to emptyset removed. 
    Example: 
    buildEquationsShapley(4,[0.25,0.25,0.25,0.25])
    """
    eqShapley = [[0.0]*(2**n) for i in range(n)]
    # shapley = [0]*n
    # 
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        #
        s = nElemsSet(mySetA, n)
        # print("s"+str(s)+",s!="+str(sfactor))
        for i in range(0,n):
            if belongSet(i, mySetA, n):
                eqShapley[i][set2IntA]=(1/s)
                ## shapley[i]=shapley[i] + (1/s)*(moebius[set2IntA])
                # print(shapley)
    neqShapley = list(map(lambda x:x[1:], eqShapley))
    return((neqShapley, shapley))

# qp(P, q, G=None, h=None, A=None, b=None, solver=None, kktsolver=None, initvals=None, **kwargs)
#    Solves a quadratic program
#    
#        minimize    (1/2)*x'*P*x + q'*x
#        subject to  G*x <= h
#                    A*x = b.

# old name of this function: create_solution
def ciSolveMSE (matrixDataWithOutput, belief=False, kAdditive=False, shapleyValue=False): 
    """
    Function: 
    Return the solution of the optimization problem 
    (i.e., the Moebius transform of a fuzzy mesure that minimizes the error -- MSE)
    Example:
    data85_B_sol_prob = ciSolveMSE(data85, belief=False, kAdditive=2)
    """
    n = len(matrixDataWithOutput[0])-1
    Q,q = buildOFMoebius (matrixDataWithOutput)
    nVars = len(Q)
    mQ=opt.matrix(np.array(Q))
    mb=opt.matrix(q, (nVars, 1))
    G,h = buildEquationsMoebius (n)
    if belief:
        posG, posH = buildEquationsMoebiusPositive(nVars)
        G = G + posG
        h = h + posH
    nConstraints = len(G)
    mG=opt.matrix(np.array(G))
    mh=opt.matrix(h, (nConstraints, 1))
    addOne = [[1.0]*nVars]
    addOneLeft = [1.0]
    if shapleyValue!=False:
        shapleyEqLeft, shapleyEqRight = buildEquationsShapley (n, shapleyValue)
        addOne = shapleyEqLeft      ## If shapley, we do not need sum m = 1
        addOneLeft = shapleyEqRight ##    it is implicit as sum shapley = 1
        #print("nVars="+str(nVars))
        #print(addOne)
        #print(addOneLeft)
        #print("addOne="+str(len(addOne)))
        #print("addOneLeft="+str(len(addOneLeft)))
        #print(len(addOne[0]))
        #print(len(addOne[1]))
    if kAdditive!=False:
        zeroMEq, zeroMLeft = buildEquationsKAdditive(n, kAdditive)
        addOne = addOne + zeroMEq
        addOneLeft = addOneLeft + zeroMLeft
        # print("addOne=")
        # print(addOne)
    mEA = opt.matrix(transpose(addOne))
    mEb = opt.matrix(addOneLeft, (len(addOne), 1))
    #print("mEA="+str(mEA))
    #print("mEb="+str(mEb))
    sol = opt.solvers.qp(mQ, mb, G=mG, h=mh, A=mEA, b=mEb)
    return(sol)

def ciFindMoebiusCIwithMSE (matrixDataWithOutput, belief=False, kAdditive=False, shapleyValue=False):
    """
    Function: 
    Find the Moebius transform that minimizes MSE for a dataset. 
    Returns:(0) the Moebius transform, (1) the solution of the optimization problem 
    Example:
    data85_B_Mob = ciFindMoebiusCIwithMSE(data85, belief=False, kAdditive=2)
    """
    this_sol = ciSolveMSE (matrixDataWithOutput, belief=belief, kAdditive=kAdditive, shapleyValue=shapleyValue)
    return (([0]+list(this_sol['x']), this_sol))

def ciFindMoebius (matrixDataWithOutput, belief=False, kAdditive=False, shapleyValue=False):
    """
    Function: 
    Find the Moebius transform that minimizes MSE for a dataset. 
    Example:
    data85_B_Mob = ciFindMoebius(data85, belief=False, kAdditive=2)
    """
    this_sol = ciSolveMSE (matrixDataWithOutput, belief=belief, kAdditive=kAdditive, shapleyValue=shapleyValue)
    return ([0]+list(this_sol['x']))


## Version: 20230623
def ciModel_of (matrixDataWithOutput, moebius): ## (numInstancesClass, classesPerDevice, nCopies):
    """
    Function: value of the objective function 
    Output: 
       err:  sum_i (CI(x_i)-y_i)^2   -- computed in matrix form from Mobius
       erSq: sum_i (CI(x_i)-y_i)^2
       vErr: [(CI(x_1)-y_1), ..., (CI(x_n)-y_n)]
       prediction: [CI(x_1), ..., CI(x_n)]
    Example: 
       ciModel_of(data85, data85Moebius)
    """
    Q,q = buildOFMoebius (matrixDataWithOutput)
    n = len(matrixDataWithOutput[0])-1
    vecy = list(map(lambda row: row[n], matrixDataWithOutput))
    #print("partial results")
    #print("1:"+str(matrixProduct(transpose(vec2mat(vecy)),vec2mat(vecy))))
    #print("2:"+str(matrixProduct(transpose(vec2mat(q)), vec2mat(moebius[1:]))))
    #print("3:"+str(matrixProductConstant(matrixProduct(transpose(vec2mat(moebius[1:])),matrixProduct(Q,vec2mat(moebius[1:]))),0.5)))
    err = matrixSum(matrixProduct(transpose(vec2mat(vecy)),vec2mat(vecy)),
                    matrixSum(matrixProduct(transpose(vec2mat(q)), vec2mat(moebius[1:])),
                              matrixProductConstant(matrixProduct(transpose(vec2mat(moebius[1:])),
                                                                  matrixProduct(Q,vec2mat(moebius[1:]))),0.5)))
    prediction = list(map(lambda inps: ciFromMoebius (inps[:-1], moebius) , matrixDataWithOutput))
    vErr = list(map(lambda inps, out: ((ciFromMoebius (inps[:-1], moebius)) - out), matrixDataWithOutput, vecy))
    erSq = vectorProduct(vErr, vErr)
    return(err, erSq, vErr, prediction)

def ciFindOWAweights (matrixDataWithOutput):
    """ 
    Function: 
    Find the OWA weights that minimizes MSE for the OWA. 
    This can be used later to build a symmetric measure. 
    Example: 
    ciFindOWAweights (data85)
    """
    sortedMatrix = list(map(lambda r: sorted(r[:-1], reverse=True)+[r[-1]], matrixDataWithOutput))
    mobius = ciFindMoebius(sortedMatrix, belief=False, kAdditive=1)
    symmetricWeights = list(map(lambda i:mobius[2**i],range(len(matrixDataWithOutput[0])-1)))
    symmetricWeights.reverse()
    return symmetricWeights

def ciFindMoebiusOWA (matrixDataWithOutput):
    """ 
    Function: 
    Find the Moebius transform that minimizes MSE for the OWA. 
    This is a symmetric measure. 
    """
    owaWeights = ciFindOWAweights (matrixDataWithOutput)
    fm = ciBuildSymmetric (owaWeights)
    return fromFM2Moebius (fm, len(matrixDataWithOutput[0])-1)

def ciFindWMweights (matrixDataWithOutput):
    """ 
    Function: 
    Find the OWA weights that minimizes MSE for the OWA. 
    This can be used later to build a symmetric measure. 
    Example: 
    ciFindWMweights (data85)
    """
    mobius = ciFindMoebius(matrixDataWithOutput, belief=False, kAdditive=1)
    weights = list(map(lambda i:mobius[2**i],range(len(matrixDataWithOutput[0])-1)))
    weights.reverse()
    return weights


## -----------------------------------------------------
## Functions to work with fuzzy measures
## -----------------------------------------------------
## 
##   Functions to build a transform from a measure
##         and to build the measure from a transform
##   Function to normalize a fuzzy measure

def fromMoebius2FM (moebius, n):
    """
    Function: 
    It computes the non-additive measure associated to a Moebius transform
    Example: 
    fromMoebius2FM(fromFM2Moebius(fmW,4),4)
    fromMoebius2FM(fromFM2Moebius(fmM,4),4)
    """
    fm=[0]*(2**n)
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        # Consider B subseteq A, they are in positions 0..set2IntA
        mySetB = firstSet(n)
        set2IntB = 0
        while not(equalSets(mySetA, mySetB,n)):
            mySetB = nextSet(mySetB, n)
            set2IntB = set2IntB+1
            if subseteq(mySetB, mySetA, n):
                fm[set2IntA]=fm[set2IntA]+moebius[set2IntB]
    return(fm)


def fromFM2Moebius (fm, n):
    """
    Function: It computes the Moebius transform of a non-additive measure
    Example: 
    fromFM2Moebius(fmW,4)
    fromFM2Moebius(fmM,4)
    """
    moebius=[0]*(2**n)
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        # Consider B subseteq A, they are in positions 0..set2IntA
        mySetB = firstSet(n)
        set2IntB = 0
        while not(equalSets(mySetA, mySetB,n)):
            mySetB = nextSet(mySetB, n)
            set2IntB = set2IntB+1
            if subseteq(mySetB, mySetA, n):
                sign = (-1)**(abs(nElemsSet(mySetA,n) - nElemsSet(mySetB,n)))
                moebius[set2IntA]=moebius[set2IntA]+sign*fm[set2IntB]
    return(moebius)


## In java this function is called: computeMeasureMax
def fromMaxSum2FM (maxSum, n):
    """
    Function: This function transforms a (max,+)-transform into a FM
    Example: 
    fromMaxSum2FM([0,1,2,5,10,3,7,11],3)
    fromMaxSum2FM([0.0,1.0,2.0,5.0,10.0,3.0,7.0,3.0],3)
    fromMaxSum2FM([0.0,0.1,0.4,0.2,0.5,0.0,0.2,0.3],3)
    """
    measure = [0]*len(maxSum)
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        # Consider B subseteq A, they are in positions 0..set2IntA
        mySetB = firstSet(n)
        set2IntB=0
        thisMeasure = 0
        while set2IntB<(set2IntA-1):
            mySetB = nextSet(mySetB,n)
            set2IntB = set2IntB+1
            if subseteq(mySetB, mySetA, n):
                thisMeasure = max(thisMeasure, measure[set2IntB])
        measure[set2IntA]=thisMeasure + maxSum[set2IntA]
    return (measure)


def fromFM2MaxSum (fm, n):
    """
    Function: This function transforms a FM into a (max,+)-transform
    Example: 
    fromFM2MaxSum([0, 1, 2, 7, 10, 13, 17, 28],3)
    fromFM2MaxSum([0, 1.0, 2.0, 7.0, 10.0, 13.0, 17.0, 20.0],3)
    fromFM2MaxSum([0, 0.1, 0.4, 0.6000000000000001, 0.5, 0.5, 0.7, 1.0],3)
    """
    msTransf = [0]*len(fm)
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        # Consider all subsets of mySet with one element less, take the max of fm
        fmNoElement = fm[set2Int(mySetA)]
        mySetACopyWElem = mySetA.copy()
        maxFMofSubsets = 0
        for i in range(0,n):
            if belongSet(i, mySetA, n):
                mySetACopyWElem = removeElement2Set(i, mySetACopyWElem, n)
                maxFMofSubsets = max(maxFMofSubsets, fm[set2Int(mySetACopyWElem)])
                # print("i="+str(i)+", sfactor="+str(sfactor)+" diff="+str(fm[set2Int(mySetACopyWElem)] - fmNoElement))
                mySetACopyWElem = addElement2Set(i, mySetACopyWElem, n)
        msTransf[set2IntA]=fm[set2IntA]-maxFMofSubsets
    return (msTransf)



def normalizeFM (fm):
    """
    normalizeFM ([0.0,1.0,2.0,5.0,10.0,3.0,7.0,3.0])
    """
    lastPos = len(fm)-1
    nFM = fm.copy()
    for i in range(0,len(fm)):
        nFM[i]=fm[i]/fm[lastPos]
    return(nFM)


## -----------------------------------------
## --   Indices and evaluation methods   ---
## -----------------------------------------

def ciShapley (fm,n):
    """ Function: Computation of the Shapley value of the fuzzy measure/game
    Example: 
    ciShapley(fmM,4)
    ciShapley(fmO,4)
    ciShapley(fmW,4)
    """
    shapley = [0]*n
    nf = math.factorial(n)
    # 
    mySetA = firstSet(n)
    #set2IntA = 0
    while notLastSet(mySetA):
        #
        s = nElemsSet(mySetA, n)
        sfactor = math.factorial(s)*(math.factorial(n-s-1))/nf
        # print("s"+str(s)+",s!="+str(sfactor))
        fmNoElement = fm[set2Int(mySetA)]
        mySetACopyWElem = mySetA.copy()
        for i in range(0,n):
            if not(belongSet(i, mySetA, n)):
                mySetACopyWElem = addElement2Set(i, mySetACopyWElem, n)
                # print("index="+str(set2Int(mySetACopyWElem)))
                shapley[i]=shapley[i] + sfactor*(fm[set2Int(mySetACopyWElem)] - fmNoElement)
                # print("i="+str(i)+", sfactor="+str(sfactor)+" diff="+str(fm[set2Int(mySetACopyWElem)] - fmNoElement))
                mySetACopyWElem = removeElement2Set(i, mySetACopyWElem, n)
                #
                ## print(shapley)
        mySetA = nextSet(mySetA, n)
        #set2IntA = set2IntA+1
    return(shapley)

def ciShapleyMoebius (moebius, n):
    """ Function: Computation of the Shapley value of the fuzzy measure/game
    given in terms of its Moebius transform
    Example: 
    ciShapleyMoebius(fromFM2Moebius(fmZero,4),4)
    ciShapleyMoebius(fromFM2Moebius(fmOne,4),4)
    ciShapleyMoebius(fromFM2Moebius (fmM,4),4)
    ciShapleyMoebius(fromFM2Moebius (fmO,4),4)
    ciShapley(fmM,4)
    ciShapley(fmO,4)
    """
    shapley = [0]*n
    # 
    mySetA = firstSet(n)
    set2IntA = 0
    while notLastSet(mySetA):
        mySetA = nextSet(mySetA, n)
        set2IntA = set2IntA+1
        #
        s = nElemsSet(mySetA, n)
        # print("s"+str(s)+",s!="+str(sfactor))
        for i in range(0,n):
            if belongSet(i, mySetA, n):
                shapley[i]=shapley[i] + (1/s)*(moebius[set2IntA])
                # print(shapley)
    return(shapley)

def ciInteraction (fm,n,vecIndicesS):
    """ 
    Function: Computation of the interaction index of the 
    fuzzy measure/game (fm), the number of elements in the 
    reference set (n), and a given set S (vecIndicesS). 
    The set is described by a list of indices. 
    Example: 
    ciInteraction(fmM,4,[0])
    ciInteraction(fmM,4,[1])
    ciInteraction(fmM,4,[2])
    ciInteraction(fmM,4,[3])
    ciInteraction(fmM,4,[0,1])
    ciInteraction(fmO,4,[0])
    ciShapley(fmO,4)
    ciInteraction(fmO,4,[0,1])
    ciInteraction(fmO,4,[0,2])
    """
    vecIndicesNotS = [i for i in range(0,n) if i not in vecIndicesS]
    lNotS = len(vecIndicesNotS)
    lS    = len(vecIndicesS)
    # print("S="+str(vecIndicesS)+", nS="+str(vecIndicesNotS))
    nsfactor = math.factorial(n-lS+1)
    interaction = 0
    # 
    mySetTofNotS = firstSet(lNotS)
    # set2IntA = 0
    # for cada subconjunt de NOTS:
    # while notLastSet(mySetTofNotS):
    for ii in range(0,2**lNotS):
        lT = nElemsSet(mySetTofNotS, lNotS)
        ssfact = math.factorial(n-lT-lS)*math.factorial(lT)/nsfactor
        # print("ssFact="+str(ssfact))
        ss = 0
        mySetKofS = firstSet(lS)
        # set2IntA = 0
        # for cada subconjunt k de S:
        # while notLastSet(mySetKofS):
        for jj in range(0,2**lS):
            lK = nElemsSet(mySetKofS, lS)
            TcupK = setIndex2Set(setAllElems2ElemsInSet (mySetKofS, vecIndicesS)+
                                 setAllElems2ElemsInSet (mySetTofNotS, vecIndicesNotS),n)
            # print("TcupK="+str(TcupK)+",mu(TcupK)="+str(fm[set2Int(TcupK)]))
            ss = ss + (-1)**(lS-lK) * fm[set2Int(TcupK)]
            mySetKofS = nextSet(mySetKofS, lS)
            #set2IntA = set2IntA+1
        interaction = interaction + ssfact*ss
        mySetTofNotS = nextSet(mySetTofNotS, lNotS)
        #set2IntA = set2IntA+1
        # print(mySetTofNotS)
    return(interaction)

def ci2Interactions (fm,n):
    """
    Function: It computes Shapley (i.e., interaction index of singletons)
    and the interaction index of two indices. 
    Output: A matrix where position (i,i) represents the Shapley index 
    and position (i,j) for i neq j represents an interaction index. 
    Examples: 
    ciShapley(fmO,4)
    ciInteraction(fmO,4,[0,1])
    ciInteraction(fmO,4,[0,2])
    ciInteraction(fmO,4,[0,3])
    ciInteraction(fmM,4,[0,3])
    ci2Interactions (fmO,4)
    """
    shapley = ciShapley (fm,n)
    mat = []
    for i in range(0,n):
        matRow = []
        for j in range(0,n):
            if j<i:
                matRow.append(mat[j][i])
            elif j==i:
                matRow.append(shapley[i])
            else:
                matRow.append(ciInteraction(fm,n,[i,j]))
        mat.append(matRow)
    return(mat)

def ciTauI (fm, n, iIndex):
    """
    Function: tau index for iIndex in 1..n for the fuzzy measure fm
    """
    #return list(map(lambda aPerm: fm[setIndex2Int(aPerm[:i], n)] - fm[setIndex2Int(aPerm[:i-1], n)],
    diffs = list(map(lambda aPerm:
                     fm[setIndex2Int(aPerm[:iIndex], n)] - fm[setIndex2Int(aPerm[:iIndex-1], n)],
                     permutations(list(range(0,n)))))
    return (sum(diffs)/math.factorial(n))

def ciTau (fm, n):
    """ Function: Computation of the Tau index of a measure
    Example: 
    ciTau(fmW,4)
    ciTau(fmM,4)
    ciTau(fmO,4)
    """
    return list(map(lambda i: ciTauI(fm, n, i), range(1, n+1)))


def ciOrnessTau (fm, n):
    """
    Examples: 
    # ciOrnessTau(fmW,4)  ## 0.6295972222222225
    # ciOrnessTau(fmM,4)  ## 0.5
    # ciOrnessTau(fmO,4)  ## 0.6666333333333336
    """
    tau = list(map(lambda i:ciTauI(fm,n,i), range(1,n+1)))
    return (sum(map(lambda t, j: (n-j)*t, tau, range(1,n+1))) / (n-1))


## ----------------------------------
## Definition of fuzzy measures
## ----------------------------------

def ciDefDP (Q, w):
    """ 
    Function: From a distortion function / quantifier, and 
    a probability distribution, returns the measure as a function: 
    set -> Q(w(set))
    Example: 
    """
    return lambda mySet: Q(vectorProduct(mySet, w))

def ciBuildDP (Q,w):
    """
    Function: Builds Distorted Probability
    Input: distortion/quantifier and probability, 
    Ouput: an array for FM
    Examples: 
    ciBuildDP(lambda x:x, [0.1,0.2,0.3,0.4])
    ciBuildDP(lambda x:x*x, [0.1,0.2,0.3,0.4])
    """
    mySet = firstSet(len(w))
    allSetsM = [ciDefDP(Q, w)(mySet)]
    while notLastSet(mySet):
        mySet = nextSet(mySet,len(w))
        allSetsM.append(ciDefDP(Q, w)(mySet))
    return allSetsM

def ciBuildFM (fFM, n):
    """
    Function: Builds an array fuzzy measure from a function set->measure
    Input: 
    # function fFM: set -> mu(set), 
    # n: number of elements in the reference set    
    Ouput: an array for FM
    Examples: 
    ##          V, S,   MP,  C,  L, KD,  M, SD     (.SE 2022)
    mandats = [24, 107, 18, 24, 16, 19, 68, 73]
    fmRiksdagsval2022 = ciBuildFM (fmWinningCoalition(mandats,175),8)
    ciShapley(fmRiksdagsval2022, 8)
    ci2Interactions(fmRiksdagsval2022,8)
    """
    mySet = firstSet(n)
    allSetsM = [fFM(mySet)]
    while notLastSet(mySet):
        mySet = nextSet(mySet,n)
        allSetsM.append(fFM(mySet))
    return allSetsM


def ciBuildSymmetric (owaWeights):
    """
    Function: Given OWA weights, it builds a summetric fuzzy measure
    Example: 
    ciBuildSymmetric ([0.1,0.2,0.3,0.4])
    """
    n = len(owaWeights)
    accOwaWeights = [0]*n
    accOwaWeights[0] = owaWeights[0] 
    for i in range(1,n):
        accOwaWeights[i] = accOwaWeights[i-1] + owaWeights[i]
    mySet = firstSet(n)
    allSetsM = [0]
    while notLastSet(mySet):
        mySet = nextSet(mySet,n)
        ### allSetsM.append(fFM(mySet))
        allSetsM.append(accOwaWeights[sum(mySet)-1])
    return allSetsM




## ----------------------------------------------------------------------------------
## Examples of fuzzy measures:
## ----------------------------------------------------------------------------------

## 
## Grabisch measure: {M,P,L}  (Example 5.12 in (Torra & Narukawa,2007))
##              -,L,  P,   PL, M,   ML, MP, MPL
fmG = [0,0.3,0.45,0.9,0.45,0.9,0.5,1]
fmGMoebius = fromFM2Moebius (fmG, 3)
## Grabisch/Sugeno lambda measure: {M,P,L}
##                            (Example 5.27 in (Torra & Narukawa,2007))
fmL = [0,0.1842,0.2237,0.5061911,0.2237,0.5061911,0.5667687,1]
## \mu_{p,w}  -- 5 elements   (Example 5.66, Table 5.9 in (Torra & Narukawa,2007))
fmT = [0.0,0.04296875,0.1375,0.2375,0.06909722,0.1375, 0.3,0.5,
       0.18333333,0.3,0.61666666,0.7625,0.38333333,0.61666666,0.81666666,0.9,
       0.1,0.18333333,0.38333333,0.61666666,0.2375,0.38333333,0.70000000,0.81666666,
       0.5,0.7,0.8625,0.93090277,0.7625,0.8625,0.95703125,1.0]
## \mu_{WM}   -- 4 elements   (Table 6.4 in (Torra & Narukawa,2007)) // w=(0.05,0.15,0.3,0.5)
fmM = [0,0.05,0.15,0.2, 0.3, 0.35,0.45,0.5,0.5,0.55,0.65,0.7,0.8,0.85,0.95,1.0]
## \mu_{OWA}  -- 4 elements   (Table 6.4 in (Torra & Narukawa,2007))
fmO = [0.0,0.3333,0.3333,0.6666,0.3333,0.6666,0.6666,1.0,0.3333,0.6666,0.6666,1.0,0.6666,1.0,1.0,1.0]
## \mu_{WOWA} -- 4 elements   (Table 6.4 in (Torra & Narukawa,2007))
fmW = [0.0,0.0666,0.2,0.2666,0.4,0.4666,0.6,0.6666,0.6666,0.7333,0.8666,0.9333,1.0,1.0,1.0,1.0]

fmZero = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0]
fmOne  = [0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]

def fmOneDim (n):
    return [0.0] + [1.0]*(2**n-1)

def fmWinningCoalition (seats, minimumVotes):
    """
    mandats = [24, 107, 18, 24, 16, 19, 68, 73]
    fmWinningCoalition (mandats, 175)([0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0])
    ciBuildFM (fmWinningCoalition(mandats,175),len(mandats))
    """
    return (lambda mySet: 1.0 if ((vectorProduct(mySet, seats)) >= minimumVotes) else 0.0)

def fmUnanimitySet (unanimitySet):
    """
    Build a fuzzy measure from a set: an unanimity set. 
    fmUnanimitySet ([0,1,1,0,0])
    ciBuildFM (fmUnanimitySet ([0,1,1,0,0]),5)
    """
    return fmWinningCoalition (unanimitySet, sum(unanimitySet))


fmOneD5 = fmOneDim (5)
   
fm4 = [0.0,  0.1, 0.2, 0.2, 
       0.3, 0.31, 0.32,0.4,
       0.4, 0.41, 0.42,0.42,
       0.7, 0.70, 0.70,1.00]

## Measures from the following paper: 
## V. Torra, (Max, +)-transforms and genetic algorithms for fuzzy measure identification, Fuzzy sets and systems.
##
## Solution found for the problem dataLiEtal2013, when the model is the Choquet integral (using GA as described in the paper): 
## 
fmFSS2022CI = [0.0,0.13274336283185842,0.17699115044247787,0.336283185840708,0.19469026548672566,
               0.4247787610619469,0.415929203539823,0.48672566371681414,0.2743362831858407,0.5132743362831859,
               0.49557522123893805,0.6814159292035398,0.5398230088495575,0.7522123893805309,0.6548672566371682,1.0]
##
## Solution found for the problem dataLiEtal2013, when the model is the Sugeno integral (using GA as described in the paper): 
## 
fmFSS2022SI = [0.0,0.11688311688311688,0.06493506493506493,0.16883116883116883,0.03896103896103896,
               0.5064935064935064,0.19480519480519481,0.5064935064935064,0.3116883116883117,0.6883116883116883,
               0.6233766233766234,0.7142857142857143,0.4935064935064935,0.8831168831168831,0.8961038961038961,1.0]
##
## Solution found using the software in this package, when the model is the Choquet integral: 
## Execution: 
##    dataLi_sol = ciFindMoebiusCIwithMSE(dataLiEtal2013)
##    dataLi_Mob = dataLi_sol[0]                     # Select Mobius transform
##    dataLi_FM  = fromMoebius2FM (dataLi_Mob, 4)
## 
dataLi_FM = [0, 0.11902721152246414, 0.2192511179972149, 0.36944035435601363, 0.278254585089656, 0.4216047035265328,
             0.38662134767666456, 0.4859785348532103, 0.40336076705004414, 0.5130985615950798, 0.5008868475489031,
             0.6440459784023589, 0.5477764557651374, 0.6801695626428382, 0.5592789623500151, 1.0]




## Data for measure identification ("ML", "P", "M", "L", "G", "Subjective evaluation")
##                            (Table 8.5 in (Torra & Narukawa,2007))

data85 = [
    [     0.8  ,  0.9  ,  0.8 ,   0.1 ,   0.1, 0.7],
    [     0.7  ,  0.6  ,  0.9 ,   0.2 ,   0.3, 0.6],
    [     0.7  ,  0.7  ,  0.7 ,   0.2 ,   0.6, 0.6],
    [     0.6  ,  0.9  ,  0.9 ,   0.4 ,   0.4, 0.8],
    [     0.8  ,  0.6  ,  0.3 ,   0.9 ,   0.9, 0.8],
    [     0.2  ,  0.4  ,  0.2 ,   0.8 ,   0.1, 0.3],
    [     0.1  ,  0.2  ,  0.4 ,   0.1 ,   0.2, 0.1],
    [     0.3  ,  0.3  ,  0.3 ,   0.8 ,   0.3, 0.4],
    [     0.5  ,  0.2  ,  0.1 ,   0.2 ,   0.1, 0.3],
    [     0.8  ,  0.2  ,  0.2 ,   0.5 ,   0.1, 0.5]]

## Data for measure identification from: 
##   Li, C., Zeng-tai, G., Gang, D. (2013) Genetic Algorithm Optimization for Determining Fuzzy
##   Measures from Fuzzy Data, Journal of Applied Mathematics, Article ID 542153.
##   https://www.hindawi.com/journals/jam/2013/542153/
## 
## Also used in Torra (FSS, 2022), Table 1. 

dataLiEtal2013 = [
    [0, 	0,	0.3,	1   , 0.2],
    [0.3,	0.8,	0.5,	0.6 , 0.5],
    [0.7,	0.9,	1,	0.7 , 0.8],
    [0.1,	0,	0.4,	0.3 , 0.15],
    [0.5,	0.4,	0.5,	0.6 , 0.5],
    [1,   	0.3,	0.8,	0.8 ,  0.7],	
    [0.3,	1,	0,	1   , 0.5],	
    [0.9,	0.6,	0.5,	0.7 , 0.7],	
    [0.6,	0.8,	0.2,	0.4 , 0.5],	
    [0.4,	1,	0.5,	0   , 0.4],	
    [0.8,	0.4,	1,	0.3 , 0.6],	
    [1, 	0.8,	0.7,	0.5 , 0.75],	
    [0.2,	0.6,	1,	0.7 , 0.5],	
    [0, 	0.3,	0.2,	0.9 , 0.2],	
    [0.4 ,	0.2,	0.8,	0   , 0.3]]




## -----------------------------------------------------
## andness and orness
## -----------------------------------------------------

## python -m pip install cubature

## pip install --upgrade numba
## pip install --upgrade numpy
## pip install --upgrade cubature

import cubature

from cubature import cubature

from scipy import integrate

def ciOrness (agop, dim):
    return ciOrnessCubature (agop, dim)

def ciOrnessCubature (agop, dim):
    """
    ciOrness(ciAM,4)
    ciOrnessCubature(ciAM,4)
    ciOrness(lambda x: ciWM(x,[0.1,0.2,0.3,0.4]),4)
    """
    vol  = cubature(agop, dim, 1, [0]*dim, [1]*dim)
    # vMin = integrate.nquad(min, [[0,1]]*dim)
    # vMax = integrate.nquad(max, [[0,1]]*dim)
    vMin = ciVolumeMin(dim)
    vMax = ciVolumeMax(dim)
    return ((vol[0]-vMin) / (vMax-vMin))


def ciOrnessIntegrateNQUAD (agop, dim):
    """
    ciOrness(ciAM,4)
    ciOrness(lambda x: ciWM(x,[0.1,0.2,0.3,0.4]),4)
    """
    vol  = integrate.nquad(ciWrap4F(agop), [[0,1]]*dim)
    # vMin = integrate.nquad(min, [[0,1]]*dim)
    # vMax = integrate.nquad(max, [[0,1]]*dim)
    vMin = ciVolumeMin(dim)
    vMax = ciVolumeMax(dim)
    return ((vol[0]-vMin) / (vMax-vMin))

def ciAndness (agop, dim):
    return 1-ciOrnessCubature (agop, dim)

## Function that permits to have a function receiving a list instead of n parameters
def ciWrap4F(f):
    return lambda *args: f(list(args))

def ciVolume (agop, dim):
    return integrate.nquad(ciWrap4F(agop), [[0,1]]*dim)

def ciVolumeMin (dim):
    return 1/(dim+1)

def ciVolumeMax (dim):
    return dim/(dim+1)