OR-Tools CP-SAT实现射线法遇除数错误及布尔变量切换问题

我正尝试用OR-Tools CP-SAT实现射线法（对应pnpoly算法），用于判断点是否在多边形内部。先完成了2顶点+1个测试点的实现，但扩展到多顶点场景时，虽设置除数diff_3的定义域不含0，却触发The Divisior cannot span accross zero in constraint错误；同时不清楚如何在CP-SAT中切换布尔变量。参考Laurent的方案实现多顶点版本后仍存在问题，特寻求技术帮助。

原pnpoly算法代码 #

int pnpoly(int nvert, double *vertx, double *verty, double testx, double testy)
{
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
     (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
 }

2顶点场景实现代码 #

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

# 2 points
vert_x1 = model.NewIntVar(1, 1, 'x1')
vert_y1 = model.NewIntVar(4, 4, 'y1')

vert_x2 = model.NewIntVar(6, 6, 'x2')
vert_y2 = model.NewIntVar(0, 0, 'y2')

# 1 test point
testx = model.NewIntVar(2, 2, 'xc')
testy = model.NewIntVar(2, 2, 'yc')

result = model.NewBoolVar("result")

# booleans for y-axis conditions
b1 = model.NewBoolVar('b1')
b2 = model.NewBoolVar('b2')
b3 = model.NewBoolVar('b3')

#booleans for x-axis conditions
c1 = model.NewBoolVar('c1')

# booleans for final decision
d1 = model.NewBoolVar('d1')
d2 = model.NewBoolVar('d2')

#temp_variable for x_conditions
diff_1 = model.NewIntVar(-100, 100, 'diff_1')
diff_2 = model.NewIntVar(-100, 100, 'diff_2')
diff_3 = model.NewIntVarFromDomain(cp_model.Domain.FromIntervals([[-100, -1]]), 'diff_3')
add_1 = model.NewIntVar(-100, 100, 'add_1')
div_1 = model.NewIntVar(-100, 100, 'div_1')
prod_1 = model.NewIntVar(-1000, 1000, 'prod_1')

"""
(verty[i] > testy) && (verty[j] > testy)
"""

#check conditions for y_axis
model.Add(vert_y1 > testy).OnlyEnforceIf(b1)
model.Add(vert_y1 <= testy).OnlyEnforceIf(b1.Not())

model.Add(vert_y2 > testy).OnlyEnforceIf(b2)
model.Add(vert_y2 <= testy).OnlyEnforceIf(b2.Not())

model.Add(b1 != b2).OnlyEnforceIf(b3)
model.Add(b1 == b2).OnlyEnforceIf(b3.Not())

"""
(testx < (vertx[j] - vertx[i]) * (testy - vert[i]) / (verty[j] - verty[i])  + vertx[i] )
"""
#check conditions for x_axis
model.Add(diff_1 == vert_x2 - vert_x1)
model.Add(diff_2 == testy - vert_y1)
model.Add(diff_3 == vert_y2 - vert_y1)
model.AddMultiplicationEquality(prod_1, diff_1, diff_2)
model.AddDivisionEquality(div_1, prod_1, diff_3 )
model.Add(add_1 == div_1 + vert_x1)

model.Add(testx < add_1).OnlyEnforceIf(c1)
model.Add(testx >= add_1).OnlyEnforceIf(c1.Not())

"""
check if both conditions are True)
"""
model.Add(c1 == 1).OnlyEnforceIf(d1)
model.Add(b3 == 1).OnlyEnforceIf(d2)

"""toggle result
(I dont know how to implement toggle)
"""
model.Add(d1 != d2).OnlyEnforceIf(result.Not())

solver = cp_model.CpSolver()
solver.parameters.log_search_progress = True
status = solver.Solve(model)

if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
    print("b1 = ", solver.Value(b1))
    print("b2 = ", solver.Value(b2))
    print("b3 = ", solver.Value(b3))
    print("c1 = ", solver.Value(c1))
    print("d1 = ", solver.Value(d1))
    print("d2 = ", solver.Value(d2))
    print("result : ", solver.Value(result))
    print("")
    print("diff_1", solver.Value(diff_1))
    print("diff_2", solver.Value(diff_2))
    print("diff_3", solver.Value(diff_3))
    print("prod_1 = ", solver.Value(prod_1))
    print("div_1", solver.Value(div_1))
    print("addd_1", solver.Value(add_1))

多顶点扩展问题及代码 #

扩展到多顶点场景时，注释掉除法相关逻辑后可正常输出，说明其他约束无问题，但保留除法逻辑会触发The Divisior cannot span accross zero in constraint错误。

diff_3_arr的定义域设置：

diff_3_arr = [model.NewIntVar(-1000, 1000, 'diff_3_arr_{i}') for i in range(num_points)]

多顶点扩展代码：

denum_pos = [model.NewIntVar(1, 100, 'denum_pos_{i}') for i in range(num_points)]
denum_neg = [model.NewIntVar(-100, -1, 'denum_neg_{i}') for i in range(num_points)]
denum_is_positive = [model.NewBoolVar('denum_is_positive_{i}') for i in range(num_points)]
target_pos = [model.NewIntVar(1, 100, 'target_pos_{i}') for i in range(num_points)]
target_neg = [model.NewIntVar(-100, -1, 'target_neg_{i}') for i in range(num_points)]


for i in range(0, num_points):
    j = i - 1 if i != 0 else num_points - 1
    model.Add(diff_1_arr[i] == xp[j] - xp[i])
    model.Add(diff_2_arr[i] == yc - yp[i])
    model.Add(diff_3_arr[i] == yp[j] - yp[i])
    model.AddMultiplicationEquality(prod_1_arr[i], diff_1_arr[i], diff_2_arr[i])

    #perform division
    model.Add(diff_3_arr[i] > 0).OnlyEnforceIf(denum_is_positive[i])
    model.Add(diff_3_arr[i] < 0).OnlyEnforceIf(denum_is_positive[i].Not())
    
    model.Add(denum_pos[i] == diff_3_arr[i]).OnlyEnforceIf(denum_is_positive[i])
    model.Add(denum_pos[i] == 1).OnlyEnforceIf(denum_is_positive[i].Not())
    model.Add(denum_neg[i] == 1).OnlyEnforceIf(denum_is_positive[i])
    model.Add(denum_neg[i] == diff_3_arr[i]).OnlyEnforceIf(denum_is_positive[i].Not())
    
    model.AddDivisionEquality(target_pos[i], prod_1_arr[i], denum_pos[i])
    model.AddDivisionEquality(target_neg[i], prod_1_arr[i], denum_neg[i])
    
    model.Add(div_1_arr[i] == target_pos[i]).OnlyEnforceIf(denum_is_positive[i])
    model.Add(div_1_arr[i] == target_neg[i]).OnlyEnforceIf(denum_is_positive[i].Not())
    model.Add(add_1_arr[i] == div_1_arr[i] + xp[i])

内容的提问来源于stack exchange，提问作者Ken Adams