Commit 5cfac261 authored by Fabrizio Detassis's avatar Fabrizio Detassis
Browse files

WIP Cplex / SMT implementation

parent 116efe62
# Input Parameters:
# - n_points -> number of examples.
# - n_classes -> number of classes
import numpy as np
from docplex.mp.model import Model as CPModel
def didi_r(x, y, pfeat):
"""
Compute the disparate impact discrimination index of a given dataset.
"""
n_points = len(y)
tot = .0
for ix_feat in pfeat:
vals = np.unique(x[:, ix_feat])
for v in vals:
mask = (1 * (x[:, ix_feat] == v)).reshape(-1, 1)
yp = mask * y
# print(yp)
Np = np.sum(mask)
if Np > 0:
tmp = (1.0 / Np) * np.sum(yp) - \
(1.0 / n_points) * np.sum(y)
tot += np.abs(tmp)
return tot
n_points = 10
pfeat = [1]
x1 = np.hstack([np.random.rand(n_points, 1),
np.full(shape=(n_points, 1), fill_value=0)])
x2 = np.hstack([np.random.rand(n_points, 1),
np.full(shape=(n_points, 1), fill_value=1)])
x = np.vstack([x1, x2])
y1 = np.random.normal(loc=0.0, scale=1.0, size=(n_points, 1))
y2 = np.random.normal(loc=3.0, scale=3.0, size=(n_points, 1))
y = np.vstack([y1, y2])
y = (y - y.min()) / (y.max() - y.min())
didi_tr = didi_r(x, y, pfeat)
constraint_value = .6 * didi_tr
print("DIDI train: %.3f -> CT value: %.3f" % (didi_tr, constraint_value))
# Build a model
mod = CPModel('Fairness problem')
# Variable declaration.
n_points = len(y)
idx_var = [i for i in range(n_points)]
z = mod.continuous_var_list(keys=idx_var, lb=0.0, ub=1.0, name='z')
# Fairness constraint: instead of adding a penalization term in the objective function - as done by
# Phebe et al - I impose the fairness term to stay below a certain threshold.
constraint = .0
abs_val = mod.continuous_var_list(keys=len(np.unique(x[:, pfeat])), name='abs_val')
# Add fairness constraint.
for ix_feat in pfeat:
vals = np.unique(x[:, ix_feat])
for i, v in enumerate(vals):
mask = 1 * (x[:, ix_feat] == v).reshape(-1, 1)
yp = mask * y
Np = np.sum(mask)
# print("Np", Np)
# print("yp", yp)
if Np > 0:
tmp = (1.0 / n_points) * mod.sum(z) -\
(1.0 / Np) * mod.sum([yp[j][0] * z[j] for j in idx_var])
# Linearization of the absolute value.
mod.add_constraint(abs_val[i] >= tmp)
mod.add_constraint(abs_val[i] >= -tmp)
constraint += mod.sum(abs_val)
mod.add_constraint(constraint <= constraint_value, ctname='fairness_cnst')
# Objective Function.
y_loss = (1.0 / n_points) * mod.sum([(y[i][0] - z[i]) * (y[i][0] - z[i]) for i in idx_var])
# y_loss = (1.0 / n_points) * mod.sum([mod.abs(y[i][0] - z[i]) for i in idx_var])
mod.minimize(y_loss)
# Solve the problem
sol = mod.solve()
# print(mod.get_constraint_by_index(0))
print("Objective fct value: %.2f" % sol.get_value(y_loss))
print("Fairness ct: ", sol.get_value(constraint))
if sol:
sat = mod.get_solve_status()
print("Status: " + str(sat))
y_new = np.array([
sol.get_value(z[i]) for i in range(n_points)
], dtype=float).reshape(-1, 1)
print()
print("Final data: ")
print("(x1, x2) \t y \t y_new")
for i in range(len(x)):
print("(%.2f, %d) \t %.2f \t %.2f" % (x[i, 0], x[i, 1], y[i], y_new[i]))
print("Original data: didi = %.2f" % didi_r(x, y, pfeat))
print("Tot: %.2f \t x2 = 0: %.2f \t x2 = 1: %.2f" % (np.mean(y),
np.mean(y[x[:, 1] == 0]),
np.mean(y[x[:, 1] == 1])))
print("New data: didi = %.2f" % didi_r(x, y_new, pfeat))
print("Tot: %.2f \t x2 = 0: %.2f \t x2 = 1: %.2f" % (np.mean(y_new),
np.mean(y_new[x[:, 1] == 0]),
np.mean(y_new[x[:, 1] == 1])))
else:
print("No solution found")
import numpy as np
from abc import ABC, abstractmethod
import utils
class AbstractConstraint(ABC):
......@@ -9,7 +10,7 @@ class AbstractConstraint(ABC):
pass
class InequalityGlobalConstraint(AbstractConstraint):
class InequalityRegGlobalConstraint(AbstractConstraint):
def __init__(self, name, sense, rhs):
self.name = name
......@@ -43,8 +44,90 @@ class InequalityGlobalConstraint(AbstractConstraint):
c = np.apply_along_axis(self._checker, axis=0, arr=y)
return np.all(c)
def __str__(self):
return "InequalityRegGlobalConstraint: y " + str(self.sense) + " " + str(self.rhs)
class InequalityClsGlobalConstraint(AbstractConstraint):
def __init__(self, name, sense, rhs):
self.name = name
self.sense = sense
self.rhs = rhs
if sense == '<=':
self._checker = self.leq
elif sense == '>=':
self._checker = self.geq
elif sense == '>':
self._checker = self.gt
elif sense == '<':
self._checker = self.lt
else:
raise ValueError("Sense " + str(sense) + " not understood!")
def leq(self, value):
return value <= self.rhs
def lt(self, value):
return value < self.rhs
def geq(self, value):
return value >= self.rhs
def gt(self, value):
return value > self.rhs
class FairnessConstraint(AbstractConstraint):
# @todo
def is_satisfied(self, x, y):
pass
classes = np.argmax(y, axis=1)
c = np.apply_along_axis(self._checker, axis=0, arr=y)
return np.all(c)
def __str__(self):
return "InequalityClsGlobalConstraint: y " + str(self.sense) + " " + str(self.rhs)
class BalanceConstraint(AbstractConstraint):
def __init__(self, name, value):
self.name = name
self.value = value
def is_satisfied(self, x, y):
lab, cnts = np.unique(y, return_counts=True)
n_points = len(y)
return np.max(cnts) <= int(np.ceil(n_points / len(lab))) + self.value
def __str__(self):
return "BalanceConstraint with value " + str(self.value)
class FairnessRegConstraint(AbstractConstraint):
def __init__(self, name, pfeat, value):
self.name = name
self.pfeat = pfeat
self.value = value
def is_satisfied(self, x, y):
didi = utils.didi_r(x, y, self.pfeat)
return didi <= self.value
def __str__(self):
return "FairnessConstraint: didi <= " + str(self.value)
class FairnessClsConstraint(AbstractConstraint):
def __init__(self, name, pfeat, value):
self.name = name
self.pfeat = pfeat
self.value = value
def is_satisfied(self, x, y):
didi = utils.didi_c(x, y, self.pfeat)
return didi <= self.value
def __str__(self):
return "FairnessConstraint: didi <= " + str(self.value)
......@@ -98,7 +98,7 @@ f_LAP = GrpcFunction(host, port, lap_uri)
# Finally, we can apply mean square error to data:
weights = f_MSE.apply(example_data)
# weights = f_MSE.apply(example_data)
# print("Weights:", weights)
......
......@@ -27,7 +27,7 @@ class MovingTarget(ABC):
M = self.initialize_ext(d)
self.add_constraints(M, C, d)
self.add_constraints(M, C, d[0], d[1])
# self.set_loss_function(M, L, d)
# d_data = self.get_pysmt_data(y_k)
......@@ -49,7 +49,7 @@ class MovingTarget(ABC):
return y_k
@abstractmethod
def add_constraints(self, M, C, y_s):
def add_constraints(self, M, C, x_s, y_s):
"""Convert collection of constraint expressions."""
@abstractmethod
......
......@@ -4,7 +4,7 @@ from moving_target_abc import MovingTarget
from pysmt.shortcuts import Symbol, LE, GE, LT, GT, Int, And, Equals, Plus, Minus, Div, Times, Max, is_sat, get_model, \
Real, Solver, Ite
from pysmt.typing import REAL
from constraint import InequalityGlobalConstraint
from constraint import InequalityRegGlobalConstraint
class SMTModel:
......@@ -97,7 +97,7 @@ class MovingTargetRegSMT(MovingTarget):
ctype, cvar, cval = c
# Store the constraint in the class.
cstr = InequalityGlobalConstraint('ct', ctype, cval)
cstr = InequalityRegGlobalConstraint('ct', ctype, cval)
self.constraints.append(cstr)
# Translate the constraint to SMT language.
......
import numpy as np
def didi_r(x, y, pfeat):
"""
Compute the disparate impact discrimination index of a given dataset.
"""
y = y.reshape(-1, 1)
n_points = len(y)
tot = .0
for ix_feat in pfeat:
vals = np.unique(x[:, ix_feat])
for v in vals:
mask = (1 * (x[:, ix_feat] == v)).reshape(-1, 1)
yp = mask * y
# print(yp)
Np = np.sum(mask)
if Np > 0:
tmp = (1.0 / Np) * np.sum(yp) - \
(1.0 / n_points) * np.sum(y)
tot += np.abs(tmp)
return tot
def didi_c(x, y, pfeat):
"""
Compute the disparate impact discrimination index of a given dataset.
"""
y = y.reshape(-1, 1)
n_points = len(y)
n_classes = len(np.unique(y))
tot = .0
for ix_feat in pfeat:
xvals = np.unique(x[:, ix_feat])
for xv in xvals:
mask = (1 * (x[:, ix_feat] == xv)).reshape(-1, 1)
# print(yp)
for c in range(n_classes):
Np = np.sum(mask)
if Np > 0:
tmp = (1.0 / Np) * np.sum(mask * (y == c)) - \
(1.0 / n_points) * np.sum(y == c)
tot += np.abs(tmp)
return tot
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment