## Software: ## (c) Vicenç Torra ## Current version: 20230316 ## ## 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, indices, etc. ## ## ## This software is about the optimal transport problem for non-additive measures. ## Use the following reference for this paper: ## V. Torra, The transport problem for non-additive measures, ## European Journal of Operational Research, 2023 ## https://doi.org/10.1016/j.ejor.2023.03.016 ## ## References on aggregation functions and non-additive measures: ## V. Torra, Y. Narukawa (2007) Modeling Decisions, Springer. ## V. Torra, Y. Narukawa, M. Sugeno (eds.) (2013) ## Non-additive measures: theory and applications, Springer. ## On the (max, \oplus)-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 ## ## Main function to compute the solution of the optimal cost ## fmtOptimalTransport (tCost, fmXT, nX, fmYT, nY) given cost table and transforms of fuzzy measures, return OT ## On costs: ## fmtBuildCostTable (cost, nX, nY) from a cost function for sets, build a cost table ## fmtCostEx10 (A,set2IntA,nX,B,set2IntB,nY) According Example 10 (Torra, EJOR, 2023) ## fmtCostPr21p (A,set2IntA,nX,B,set2IntB,nY) ## fmtCostPr21 (costXY) According Prop. 21 (Torra, EJOR, 2023) ## fmtCostMeasure (fmX, nX, fmY, nY, cost, costEmpty) According Def 23 (Torra, EJOR, 2023) ## Wasserstein distance: ## fmtWassersteinDistance (tCost, fmXT, nX, fmYT, nY): ## Functions to manipulate the ## fmtFromAssign2lol (assgn, fmX, fmY): ## fmtFromAssign2TabLatex (assgn, nX, fmX, fmXMxT, nY, fmY, fmYMxT): ## exec(open("prog.vectors.matrices.web.txt").read()) ## exec(open("prog.choquet.web.txt").read()) from scipy.optimize import linprog def fmtFromIJ2array (i,j,nColumns,nRows): #print("[i,j]="+str(i)+","+str(j)+","+str(i*nColumns + j)) return(i*nColumns + j) ## Function: ## Return the solution of the optimal transport problem ## Reference: corresponds to the solution of Def. 16 (Torra, EJOR 2023) ## Parameters: ## tCost: table of costs 2^X x 2^Y -> R ## fmXT: transform of a fuzzy measure on X ## nX: number of elements in X ## fmYT: transform of a fuzzy measure on Y ## nY: number of elements in Y ## Example: ## A set of examples is provided in the associated file ## ## NOTE: the function "fmtBuildCostTable" builds the table tCost from a function ## that given two sets computes its cost. def fmtOptimalTransport (tCost, fmXT, nX, fmYT, nY): sizeFmX = len(fmXT) sizeFmY = len(fmYT) nXY = sizeFmY * sizeFmX assign = [0]*nXY equationsLHS = [] equationsRHS = [] # Rows, given B, add all A: mySetB = firstSet(nY) # B set2IntB = 0 mySetB = nextSet(mySetB, nY) # Set empty, no need to add set2IntB = set2IntB+1 while set2IntB < len(fmYT): newRow = [0]*nXY mySetA = firstSet(nX) set2IntA = 0 while set2IntA < len(fmXT): # A newRow[fmtFromIJ2array(set2IntA, set2IntB, sizeFmX, sizeFmY)]=1 mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 equationsLHS.append(newRow) equationsRHS.append(fmYT[set2IntB]) # print("EQ:"+str(newRow)+":"+str(fmYT[set2IntB])) mySetB = nextSet(mySetB, nY) # Set empty, no need to add set2IntB = set2IntB+1 # # Rows, given A, add all B: mySetA = firstSet(nX) # A set2IntA = 0 mySetA = nextSet(mySetA, nX) # Set empty, no need to add set2IntA = set2IntA+1 while set2IntA < len(fmXT): # newRow = [0]*nXY mySetB = firstSet(nY) set2IntB = 0 while set2IntB < len(fmYT): # B newRow[fmtFromIJ2array(set2IntA, set2IntB, sizeFmX, sizeFmY)]=1 mySetB = nextSet(mySetB, nY) set2IntB = set2IntB+1 equationsLHS.append(newRow) equationsRHS.append(fmXT[set2IntA]) # print("EQ:"+str(newRow)+":"+str(fmXT[set2IntA])) mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 # print("PROBLEM") # print(equationsLHS) # print(equationsRHS) ### DEFINE PROBLEM opt = linprog(c=tCost, # A_ub=lhs_ineq, b_ub=rhs_ineq, A_eq=equationsLHS, b_eq=equationsRHS, bounds=None, method="revised simplex") ### Return the solution of the optimization problem return(opt) ## Function: ## from a function cost: 2^X x 2^Y -> R, build a table tc: 2^X x 2^Y -> R ## (returns a table that can be used by fmtOptimalTransport) ## Parameters: ## cost: the cost function. The cost function has 6 parameters ## nX: number of elements in reference set X ## nY: number of elements in reference set Y ## Example: ## fmtTCostEx10 = fmtBuildCostTable (fmtCostEx10, 3, 3) ## fmtTCostPr21 = fmtBuildCostTable (fmtCostPr21(fmtCostPr21XY), 3, 3) ## Note: the set function needs 6 parameters: ## the set A, integer representation of set A, number of elements in X ## the set B, integer representation of set B, number of elements in Y ## cost (A,set2IntA,nX,B,set2IntB,nY) -- A \subset X, B \subset Y def fmtBuildCostTable (cost, nX, nY): sizeFmX = 2**nX ## len(fmX) sizeFmY = 2**nY ## len(fmY) nXY = sizeFmY * sizeFmX ## Define cost function myCost = [0]*nXY # Rows, given A, add all B: mySetA = firstSet(nX) # A set2IntA = 0 while set2IntA < sizeFmX: ## len(fmX): mySetB = firstSet(nY) set2IntB = 0 while set2IntB < sizeFmY: ## len(fmY): # B myCost[fmtFromIJ2array(set2IntA, set2IntB, sizeFmX, sizeFmY)]=cost(mySetA,set2IntA,nX,mySetB,set2IntB,nY) mySetB = nextSet(mySetB, nY) set2IntB = set2IntB+1 mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 return myCost ## DEFINITIONS OF COST FUNCTIONS ## ## Cost functions for sets: They are functions of this form: ## cost (setA,set2IntA,nX,setB,set2IntB,nY) -> Value ## Parameters: ## setA,set2IntA,nX: the set A, integer representation of set A, number of elements in X ## setB,set2IntB,nY: the set B, integer representation of set B, number of elements in Y ## with A \subset X, B \subset Y ## Function: ## Given two sets A and B return a cost, ## based on Example 10 (Torra, EJOR, 2023) ## c(A,B)=1 for A \neq B ## c(A,B)=0 for A = B ## Example: ## fmtCostEx10([0,0,1],1,3,[0,0,1],1,3) ## fmtCostEx10([0,0,0],0,3,[0,0,0],0,3) ## fmtCostEx10([0,0,0],0,3,[0,0,1],1,3) def fmtCostEx10 (A,set2IntA,nX,B,set2IntB,nY): if A == B: cost = 0 else: cost = 1 return(cost) ## Function: ## Given two sets A and B return a cost, ## variation of Prop 21 (Torra, EJOR, 2023) ## when c(x,y) on singletons is ## c(x,y)=1 for x \neq y ## c(x,y)=0 for x = y ## with c_a(\emptyset, B) = c_a(A, \emptyset) = \infty ## c_a({x}, B) = c_a(A,{y})=\infty |A|, |B| \neq 1 ## c_a(A,B)=0 |A|, |B| > 1 ## Example: ## fmtCostPr21p([0,0,1],1,3,[0,0,1],1,3) ## fmtCostPr21p([0,0,0],0,3,[0,0,0],0,3) ## fmtCostPr21p([0,0,0],0,3,[0,0,1],1,3) def fmtCostPr21p (A,set2IntA,nX,B,set2IntB,nY): if nElemsSet(A, nX)==0 or nElemsSet(B, nY)==0: cost = 10000 elif nElemsSet(A, nX)>1 and nElemsSet(B, nY)==1: cost = 10000 elif nElemsSet(B, nY)>1 and nElemsSet(A, nX)==1: cost = 10000 elif A == B: cost = 0 elif nElemsSet(A,nX)==1 and nElemsSet(B,nY)==1: cost = 1 else: cost = 0 return(cost) ## Function: ## Given a cost for elements c: X x Y -> R ## the function builds a function cost: 2^X x 2^Y -> R ## according to Proposition 21. ## Example: ## fmtCostPr21XY = [[0,5,5],[5,0,5],[5,5,0]] ## fmtCostPr21(fmtCostPr21XY)([0,0,1],1,3,[0,0,1],1,3) ## fmtCostPr21(fmtCostPr21XY)([0,0,1],1,3,[0,0,0],0,3) ## fmtCostPr21(fmtCostPr21XY)([0,1,0],2,3,[0,0,1],1,3) def fmtCostPr21 (costXY): def myCostFM (A,set2IntA,nX,B,set2IntB,nY): if nElemsSet(A, nX)==1 and nElemsSet(B, nY)==1: i = (setSet2Index(A,nX))[0] j = (setSet2Index(B,nY))[0] # print("i,j="+str(i)+str(j)) cost = costXY[i][j] elif nElemsSet(A, nX)!=1 and nElemsSet(B, nY)==1: cost = 10000 elif nElemsSet(B, nY)!=1 and nElemsSet(A, nX)==1: cost = 10000 else: cost = 0 return(cost) return myCostFM ## Function: ## Returns a cost table built, according to Def. 23 in (Torra, EJOR, 2023), from ## a table of a cost function (on XxY) and ## two measures fmX and fmY ## Parameters: ## fmX: fuzzy measure, nX: number of elements of X (measure or transform) ## fmY: fuzzy measure, nY: number of elements of Y (measure or transform) ## cost: cost function on X x Y (table of two dimensions) ## Note: as this function only uses measure on singletons, measure or transform can be used ## Return: ## myCost: table ## Other variables: ## mXvSingl, mXiSingl: measure on singletons sorted increasing, and respective indices ## mYvSingl, mYiSingl: measure on singletons sorted increasing, and respective indices ## iX, iY: indices ## Example: ## ecsb = fmtCostMeasure (fmXp, 3, fmY, 3, [[0,1,2],[1,0,1],[2,1,0]], 1) def fmtCostMeasure (fmX, nX, fmY, nY, cost, costEmpty): mXvSingl, mXiSingl = sortWithOneListIncreasing( list(map(lambda i: fmX[setIndex2Int([i],nX)], range(0,nX))), list(range(0,nX))) mYvSingl, mYiSingl = sortWithOneListIncreasing( list(map(lambda i: fmY[setIndex2Int([i],nY)], range(0,nY))), list(range(0,nX))) sizeFmX = len(fmX) sizeFmY = len(fmY) nXY = sizeFmY * sizeFmX ## Define cost function myCost = [1]*nXY # Rows, given A, add all B: mySetA = firstSet(nX) # A set2IntA = 0 while set2IntA < len(fmX): mySetB = firstSet(nY) set2IntB = 0 while set2IntB < len(fmY): # B if set2IntB ==0 or set2IntA == 0: myCost[fmtFromIJ2array(set2IntA, set2IntB, sizeFmX, sizeFmY)] = costEmpty if set2IntB !=0 and set2IntA !=0: xSetIn = selectFv(range(0,nX), mySetA, lambda x:x==1) xSetMu = list(map(lambda i: fmX[setIndex2Int([i],nX)], xSetIn)) iMinX = argMin(xSetMu, xSetIn) # x_i s.t. arg min_{x_i \in A} mu(x_i) iX = mXiSingl.index(iMinX) # position of x_i in the ordered list ySetIn = selectFv(range(0,nY), mySetB, lambda x:x==1) ySetMu = list(map(lambda i: fmY[setIndex2Int([i],nY)], ySetIn)) iMinY = argMin(ySetMu, ySetIn) iY = mYiSingl.index(iMinY) # position of y_i in the ordered list nMinusInverseX = nX - (iX+1) nMinusInverseY = nY - (iY+1) # print((nMinusInverseX,nMinusInverseY,iMinX,iMinY,cost[iMinX][iMinY])) valueInCost = (1/2**nMinusInverseX)*(1/2**nMinusInverseY)*cost[iMinX][iMinY] ## cost(mySetA,set2IntA,nX,mySetB,set2IntB,nY) myCost[fmtFromIJ2array(set2IntA, set2IntB, sizeFmX, sizeFmY)]=valueInCost ## cost(mySetA,set2IntA,nX,mySetB,set2IntB,nY) mySetB = nextSet(mySetB, nY) set2IntB = set2IntB+1 mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 # return mXvSingl, mXiSingl, mYvSingl, mYiSingl, iX, iY, myCost return myCost ## Function: ## Wasserstein distance between two transforms of two fuzzy measures ## Parameters: ## tCost: table with costs ## fmXT: transform of the measure ## nX: num elements of set X ## fmYT: transform of the measure ## nY: num elements of set Y ## Example: ## wDisMxEx10 = fmtWassersteinDistance(fmtTCostEx10,fmXMxT,3,fmYMxT,3) ## wDisMxPr21 = fmtWassersteinDistance(fmtTCostPr21,fmXMxT,3,fmYMxT,3) def fmtWassersteinDistance (tCost, fmXT, nX, fmYT, nY): sol = fmtOptimalTransport(tCost,fmXT, nX, fmYT, nY) return sol['fun'] ## Function: ## transform an assignment to list of lists (a matrix) ## Parameters: ## assgn: the assignment ## fmX, fmY: the measures (only length of list is used) ## Example: ## fmtFromAssign2lol(solMxPr21['x'],fmX,fmY) def fmtFromAssign2lol (assgn, fmX, fmY): """ Function: From an assignment to lists of lists, elements in the same order. Example: fmtFromAssign2lol(solMxPr21['x'],fmX,fmY) """ shAssgn = list(map(lambda x:0 if x<0.0001 else x, assgn)) lol = [] for i in range(0,len(fmX)): lol.append(shAssgn[i*len(fmX):(i+1)*len(fmX):]) return(lol) ## Function: ## transform an assignment to a table in latex ## Parameters: ## assgn: the assignment ## nX, fmX, fmXMxT: dimension of X, fuzzy measure, and its transform ## nY, fmY, fmYMxT: dimension of Y, fuzzy measure, and its transform ## Example: ## fmtFromAssign2TabLatex (sol_TX_Y['x'], 3, fmX, fmXMxT, 3, fmY, fmYMxT) def fmtFromAssign2TabLatex (assgn, nX, fmX, fmXMxT, nY, fmY, fmYMxT): sizeFmX = len(fmX) sizeFmY = len(fmY) lol = fmtFromAssign2lol(assgn, fmX, fmY) string4columns = "l|l|l||"+("l"*(len(fmX))) outString = "\\begin{tabular}{"+string4columns+"} \n" ## Row Set outString = outString + "$\\nu(B)$ & $\\tau_{\\nu}$ & $set$ " mySetA = firstSet(nX) set2IntA = 0 while set2IntA < len(fmX): # A outString = outString + " & " mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 outString = outString + "\\\\ \n" # # Rows, given A, add all B: mySetB = firstSet(nY) # B set2IntB = 0 while set2IntB < len(fmY): outString = outString + str(fmY[set2IntB])+" & " + str(fmYMxT[set2IntB])+" & "+str(setSet2Index (mySetB,nY)) mySetA = firstSet(nX) set2IntA = 0 while set2IntA < len(fmX): # A outString = outString + "& "+str(assgn[fmtFromIJ2array(set2IntA, set2IntB, sizeFmX, sizeFmY)]) mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 outString = outString + "\\\\ \n" mySetA = nextSet(mySetB, nY) set2IntB = set2IntB+1 outString = outString + "\\hline \n" outString = outString + "\\hline \n" ## Row Set outString = outString + "set & & " mySetA = firstSet(nX) set2IntA = 0 while set2IntA < len(fmX): # B outString = outString + " & " + str(setSet2Index (mySetA,nX)) mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 outString = outString + "\\\\ \n" ## Row transform(mu_X)(A) outString = outString + "$\\tau_{\\mu}$ & & " mySetA = firstSet(nX) set2IntA = 0 while set2IntA < len(fmX): # B outString = outString + " & " + str(fmXMxT[set2IntA]) mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 outString = outString + "\\\\ \n" ## Row mu_X(A) outString = outString + "$\\mu(A)$ & & " mySetA = firstSet(nX) set2IntA = 0 while set2IntA < len(fmX): # B outString = outString + " & " + str(fmX[set2IntA]) mySetA = nextSet(mySetA, nX) set2IntA = set2IntA+1 outString = outString + "\\\\ \n\\end{tabular}" return(outString) ## latexTable = fmtFromAssign2TabLatex (sol1['x'], 3, fmX, fmXMxT, 3, fmY, fmYMxT) ## print(latexTable)