abs and square supported in CP SAT? sum of boolean comparisons?

[David Olmo Pérez]

Hi

I'm trying to move a problem that previously was solved by the pywrapcp (original CP solver) to CP SAT solver as I have noticed the huge performance improvement in other problems

With pywrapcp I could perform Absolute (abs) and Square (numpy.square, or 2**2) operation over solver.IntVars:

balance_horas = {}
for i in ids.keys():
    balance_horas[i] = abs(solver.Sum([var_decision[i,hor] for hor in cols_requisitos]) - jornadas[i] + bolsa_horas[i])

suma_total_bolsa_horas = solver.IntVar(0,10000,"suma bolsa horas")
solver.Add(suma_total_bolsa_horas == solver.Sum([balance_horas[i] for i in ids.keys()]))

In addition, with the original CP solver I could sum a list of boolean comparison on IntVars and then use it on a constraint:

 solver.Add(solver.Sum([var_decision[(i, turno)] > 0 for i in ids.keys()]) >= temp_personal_necesario)

I've tried to do the same with CP SAT on model.NewEnumeratedIntVar() but it seems like it is not supported. I've needed to convert the enumerated vars to Boolean Vars

Is there anyway to do any of these 3 points with CP SAT?

[Laurent Perron]

Abs: no quick path right now

x = model.NewIntVar5, 10, "x")

y = model.NewIntVar(-10, 5, "")

model.Add(x + y == 0)

IntVar abs_x = model.NewIntVar(0, 10, "");

model.AddMaxEquality(abs_x, [x, y])

Square:

IntVar square_x = model.NewIntVar(0, 100, "");

model.AddProdEquality(square_x, [x, x])

Sum(x[i] >= 0):

see https://github.com/google/or-tools/blob/master/ortools/sat/doc/channeling.md#implementing-an-if-then-else-expression

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

--
You received this message because you are subscribed to the Google Groups "or-tools-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to or-tools-discu...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[David Olmo Pérez]

Thanks Laurent, I understood Abs and Square operations

Regarding the last point (if-then-else expression) do you know if there is any SAT example using this method? I have looked over nearly all SAT examples in github but couldn't find one. I'm asking this because it is being difficult for me to understand from that link and would like to see another example

[Laurent Perron]

x = array of int vars.

non_zeroes = []

for v in x:

  b = model.NewBoolVar('%s in not zero' % x)

  model.Add(x > 0).OnlyEnforceIf(b)

  model.Add(x <= 0).OnlyEnforceIf(b.Not())

  non_zeroes.append(b)

model,.Add(sum_total_bolsa_horas == sum(non_zeroes))

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[David Olmo Pérez]

Thanks Laurent, I understood it now. I just applied this if-else strategy successfully

Regarding the Square, I'm doing something very similar to the code below, which should work fine:

z = [model.NewBoolVar("") for i in range(5)]

x = model.NewIntVar(0, 5, "x")

y = model.NewIntVar(0, 5, "y")

x_constant = 5

y_constant = 3

model.Add( x == sum([z[i] for i in range(3)]) + x_constant)

model.Add( y == sum([z[i] for i in range(3,5)]) + y_constant)

x_square = model.NewIntVar(0,1000,"square x")

y_square = model.NewIntVar(0,1000,"square y")

model.AddProdEquality(x_square, [x,x])

model.AddProdEquality(y_square, [y,y])

all_squared = model.NewIntVar(0,100000,"x and y squared")

model.Add(all_squared == sum([x_square,y_square]))

model.Minimize(all_squared)

solver = cp_model.CpSolver()

solver.Solve(model)

solver.ObjectiveValue()

But I'm getting this error:

WARNING: Logging before InitGoogleLogging() is written to STDERR

F1001 02:20:27.917237 16036 integer_expr.h:507] Not supported

*** Check failure stack trace: ***

Process finished with exit code -1073740791 (0xC0000409)

Is there something I'm missing?

Square:

IntVar square_x = model.NewIntVar(0, 100, "");

model.AddProdEquality(square_x, [x, x])

[Laurent Perron]

This is fixed in the latest version (6.9, out for python3 on pypi):

from ortools.sat.python import cp_model

model = cp_model.CpModel()

z = [model.NewBoolVar("") for i in range(5)]

x = model.NewIntVar(0, 5, "x")

y = model.NewIntVar(0, 5, "y")

x_constant = 5

y_constant = 3

model.Add( x == sum([z[i] for i in range(3)]) + x_constant)

model.Add( y == sum([z[i] for i in range(3,5)]) + y_constant)

x_square = model.NewIntVar(0, 1000, "square x")

y_square = model.NewIntVar(0, 1000, "square y")

model.AddProdEquality(x_square, [x, x])

model.AddProdEquality(y_square, [y, y])

all_squared = model.NewIntVar(0, 100000, "x and y squared")

model.Add(all_squared == sum([x_square, y_square]))

model.Minimize(all_squared)

solver = cp_model.CpSolver()

solver.Solve(model)

print (solver.ResponseStats())

outputs:

CpSolverResponse:

status: OPTIMAL

objective: 34.000000

best_bound: 34.000000

booleans: 5

conflicts: 0

branches: 2

propagations: 3

integer_propagations: 9

walltime: 0.015469

usertime: 0.007020

deterministic_time: 0.000001

Please note that you do not need all_squared, just do 

  model.Minimize(x_square + y_square)

BTW, the bounds on X look suspicious, but this is just a toy example.

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[David Olmo Pérez]

Ok Laurent, got it. Thanks again

One last question: with pywrapcp I could use solver.ScalProd (scalar product between 2 vectors, possibly a vector of decision variables and a constant vector of the same size). How can this be done in CP-SAT?

[David Olmo Pérez]

Ok nevermind, I just realized that I can replicate ScalProd behavior easily using a simple list comprehension:

model.Add(sum([var_decision[i,turno]*repartidores_con_moto[i] for i in ids.keys()]) <= 18)

Where "var_decision" are decision variables and "repartidores_con_moto" are constant vector (both of the same size) and i is an index

In case there is a more elegant/efficient way of doing this please let me know

[David Olmo Pérez]

Laurent, congratulations for your package. I'm amazed of the performance increase I've had when switching my problem to SAT CP solver. It has had a 1000% time reduction or more (difficult to say, as now the wall time is 0.01 sec, even adding complexities that I didn't have before)  

[Laurent Perron]

This is why the original copy solver is in maintenance mode :-)

[Etics Etics]

Hi,

sorry but I don't understand how can I do ABS in C#.

I try this:

IntegerExpression diff = var1 - var2;

IntegerExpression negDiff = -diff;

model.Add(diff  + negDiff  == 0);

model.AddMaxEquality(abs, new IntVar[] { diff, negDiff  });

but diff and negDiff are not IntVar.

How can I cast IntegerExpression to IntVar?

Thank you.

Bye.

[Laurent Perron]

y = abs(x)

Create neg_x another intvar with the opposite domain of x.

add (x  + neg_x == 0)

AddMaxEquality(y, new IntVar[] { x, neg_x});

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[Etics Etics]

Hi, 

the problem is that I need to get Abs of a difference and the difference of two IntVar is a IntegerExpression.

How can I get an IntVar from a difference?

And how can I get the opposite domain of the difference?

Thank you.

Bye.

[Laurent Perron]

z = abs(x - y)

intvar z1, z2

z1 == x - y

z2 == y - x

z == max(z1, z2)

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[Laurent Perron]

with domain 

D(z1) = [min(x) - max(y), max(x) - min(y)]

D(z2) = [min(y) - max(x), max(y) - min(x)]

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[David Olmo Pérez]

Laurent, 

I'm coming again since I found some problems with the ProdEquality you explained above and as it seems a similar erro with intExpr as Etics Etics.

Consider this toy example (this is happening as well in my resource allocation problem and it is exactly the same error):

from ortools.sat.python import cp_model

model = cp_model.CpModel()

multipliers = [5,5,5]
constants = [1,2,3]
num_elements = len(multipliers)
boolvars = [model.NewBoolVar("var {}".format(i)) for i in range(num_elements)]
squared = {}
minus = {}

for i in range(num_elements):
    minus[i] = model.NewIntVar(-10000,10000,"minus {}".format(i))
    model.Add(minus[i] == sum([boolvars[i]*multipliers[i] for i in range(num_elements)]) - constants[i])
    squared[i] = model.NewIntVar(-10000,100000,"squared {}".format(i))
    model.AddProdEquality(squared[i],[minus[i] ,minus[i] ])

model.Minimize(sum([squared[i] for i in range(num_elements)]))

solver = cp_model.CpSolver()
status = solver.StatusName(solver.Solve(model))

This will prompt the following error:

WARNING: Logging before InitGoogleLogging() is

written to STDERR
F1011 20:33:40.517401 13700 integer_expr.h:514] Not supported

*** Check failure stack trace: ***
Process finished with exit code -1073740791 (0xC0000409)

What is happening here?

[Laurent Perron]

I have not implemented product on integer variables spanning over 0.

Do you really need it?

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[David Olmo Pérez]

Oh I see

I have been working on employee scheduling and resource allocations problems. This is a very common operations that I had implemented in every script I did with pywrapcp. It also is used by Hakank in most of his scripts (which I used to learn)

I guess that I would not really need it if there is any way to re-formulate the problem to avoid doing this. Do you know if there is any strategy to overcome this? I haven't found any in SAT

[Laurent Perron]

What do you want to model?

--Laurent

[Laurent Perron]

A bit slower, create y = abs(x)

and then square = y * y.

Laurent Perron | Operations Research | lpe...@google.com | (33) 1 42 68 53 00

[David Olmo Pérez]

Oh this worked beautifully! Thanks Laurent

I'm copying the same example with the reformulation just in case anyone in the future wants to check how to do this. I solved it by doing the abs using your previous strategy:

z = abs(x - y)

intvar z1, z2
z1 == x - y
z2 == y - x
z == max(z1, z2)

The code:

from ortools.sat.python import cp_model

model = cp_model.CpModel()

multipliers = [5,5,5]
constants = [1,2,3]
num_elements = len(multipliers)
boolvars = [model.NewBoolVar("var {}".format(i)) for i in range(num_elements)]
squared = {}

minus_1 = {}
minus_2 = {}
abs_minus = {}
for i in range(num_elements):
    minus_1[i] = model.NewIntVar(-10000,10000,"minus {}".format(i))
    minus_2[i] = model.NewIntVar(-10000, 10000, "minus {}".format(i))
    model.Add(minus_1[i] == sum([boolvars[i]*multipliers[i] for i in range(num_elements)]) - constants[i])
    model.Add(minus_2[i] == constants[i] - sum([boolvars[i] * multipliers[i] for i in range(num_elements)]))

    abs_minus[i] = model.NewIntVar(0,10000," abs minus {}".format(i))
    model.Add(abs_minus[i] == max(minus_1[i],minus_2[i]))

    squared[i] = model.NewIntVar(-10000,100000,"squared {}".format(i))

    model.AddProdEquality(squared[i],[abs_minus[i] ,abs_minus[i] ])

model.Minimize(sum([squared[i] for i in range(num_elements)]))

solver = cp_model.CpSolver()
status = solver.StatusName(solver.Solve(model))

[David Olmo Pérez]

Ok, in this case of using Max in order to calculate the absolute vlaue,  using the python max operator will give unstable results and infeasibility when working with negative numbers in the following line:

model.Add(abs_minus[i] == max(minus_1[i],minus_2[i]))

just for the record and for other people, in my previous script it should be used model.AddMaxEquality instead of max python operator:

from ortools.sat.python import cp_model
model = cp_model.CpModel()

multipliers = [5,5,5]

constants = [-1,2,3]

num_elements = len(multipliers)
boolvars = [model.NewBoolVar("var {}".format(i)) for i in range(num_elements)]
squared = {}
minus_1 = {}
minus_2 = {}
abs_minus = {}
for i in range(num_elements):

    minus_1[i] = model.NewIntVar(-10000000000,10000000000,"minus {}".format(i))
    minus_2[i] = model.NewIntVar(-10000000000, 10000000000, "minus {}".format(i))

    model.Add(minus_1[i] == sum([boolvars[i]*multipliers[i] for i in range(num_elements)]) - constants[i])
    model.Add(minus_2[i] == constants[i] - sum([boolvars[i] * multipliers[i] for i in range(num_elements)]))

    abs_minus[i] = model.NewIntVar(0,10000000000," abs minus {}".format(i))
    model.AddMaxEquality(abs_minus[i],[minus_1[i],minus_2[i]])


    squared[i] = model.NewIntVar(0,10000000000,"squared {}".format(i))

    model.AddProdEquality(squared[i],[abs_minus[i] ,abs_minus[i] ])

# model.Minimize(sum([squared[i] for i in range(num_elements)]))

solver = cp_model.CpSolver()
status = solver.StatusName(solver.Solve(model))

print(status)