2014-10-16 91 views
4

我有一个问题可以归结为寻找一种将三角矩阵映射到跳过对角线的向量的方法。在跳过对角线的向量上映射上三角矩阵

基本上我需要使用Gecode库

// implied constraints 
for (int k=0, i=0; i<n-1; i++) 
    for (int j=i+1; j<n; j++, k++) 
    rel(*this, d[k], IRT_GQ, (j-i)*(j-i+1)/2); 

向该MiniZinc(功能性)代码

constraint 
    forall (i in 1..m-1 , j in i+1..m) 
     ((differences[?]) >= (floor(int2float((j-i)*(j-i+1))/int2float(2)))); 

把这种C代码和我需要在differences[?]弄清楚的索引。

MiniZinc是一种功能/数学语言,没有适当的循环。 因此,我必须映射那些触及所有且仅接触上三角矩阵的单元格的索引i和j,跳过其对角线,将k指向从0到任意值的单元格。

如果这是一个普通的三角矩阵(它不是),解决like this会做

index = x + (y+1)*y/2 

我处理的矩阵是一个正方形n*n矩阵指数从0到n-1,但为n*m矩阵提供更通用的解决方案将会很好。

下面是完整的Minizinc代码

% modified version of the file found at https://github.com/MiniZinc/minizinc-benchmarks/blob/master/golomb/golomb.mzn 

include "alldifferent.mzn"; 

int: m; 
int: n = m*m; 
array[1..m] of var 0..n: mark; 
array[int] of var 0..n: differences = [mark[j] - mark[i] | i in 1..m, j in i+1..m]; 

constraint mark[1] = 0; 

constraint forall (i in 1..m-1) (mark[i] < mark[i+1]); 

% this version of the constraint works 
constraint forall (i in 1..m-1 , j in i+1..m) 
    ((mark[j] - mark[i]) >= (floor(int2float((j-i)*(j-i+1))/int2float(2)))); 

%this version does not 
%constraint forall (i in 1..m-1, j in i+1..m) 
% ((differences[(i-1) + ((j-2)*(j-1)) div 2]) >= (floor(int2float((j-i)*(j-i+1))/int2float(2)))); 

constraint alldifferent(differences); 

constraint differences[1] < differences[(m*(m-1)) div 2]; 

solve :: int_search(mark, input_order, indomain, complete) minimize mark[m]; 

output ["golomb ", show(mark), "\n"]; 

感谢。

回答

4

要小心。您从该链接找到的公式包括index = x + (y+1)*y/2,包括对角线条目,并且对于下三角矩阵,我收集的是而不是您想要的。您正在寻找的确切公式是index = x + ((y-1)y)/2 (查看:https://math.stackexchange.com/questions/646117/how-to-find-a-function-mapping-matrix-indices)。

再次注意,我给你的这个公式假设你的指数:x,y,是零基。您的MiniZinc代码使用的索引为i,j 从1开始(1 < = i < = m),1 < = j < = m))。对于从1开始的指数,公式为T(i,j) = i + ((j-2)(j-1))/2。所以,你的代码应该是这样的:

constraint 
    forall (i in 1..m-1 , j in i+1..m) 
     ((distances[(i + ((j-2)*(j-1)) div 2]) >= ... 

注意(j-2)(j-1)永远是2的倍数,所以我们可以只使用整数除法与除数2(不用担心彩车转换为/)。


上面假定您使用的是矩形m*m矩阵。
推广到一个M*N矩形矩阵,一个公式可以是:

general formula

其中0 < = I <男,0 < = j的< N [如果再次,需要你的索引从1开始,将i替换为i-1,将j替换为上式中的j-1]。这涉及上三角矩阵的所有单元格以及当N> M时出现的正方形的“额外的块”。也就是说,它触及所有单元格(i,j),使得对于0而言,0 < j < = I <男,0 < = j的< N.

extra block on the side


全码:

% original: https://github.com/MiniZinc/minizinc-benchmarks/blob/master/golomb/golomb.mzn 

include "alldifferent.mzn"; 

int: m; 
int: n = m*m; 
array[1..m] of var 0..n: mark; 
array[1..(m*(m-1)) div 2] of var 0..n: differences; 

constraint mark[1] = 0; 
constraint forall (i in 1..m-1) (mark[i] < mark[i+1]); 
constraint alldifferent(differences); 
constraint forall (i,j in 1..m where j > i) 
    (differences[i + ((j-1)*(j-2)) div 2] = mark[j] - mark[i]); 
constraint forall (i,j in 1..m where j > i) 
    (differences[i + ((j-1)*(j-2)) div 2] >= (floor(int2float((j-i)*(j-i+1))/int2float(2)))); 
constraint differences[1] < differences[(m*(m-1)) div 2]; 

solve :: int_search(mark, input_order, indomain, complete) 
    minimize mark[m]; 

output ["golomb ", show(mark), "\n"]; 

下三角版本(以前面的代码和交换i和j,其中necessar y):

% original: https://github.com/MiniZinc/minizinc-benchmarks/blob/master/golomb/golomb.mzn 

include "alldifferent.mzn"; 

int: m; 
int: n = m*m; 
array[1..m] of var 0..n: mark; 
array[1..(m*(m-1)) div 2] of var 0..n: differences; 

constraint mark[1] = 0; 
constraint forall (i in 1..m-1) (mark[i] < mark[i+1]); 
constraint alldifferent(differences); 
constraint forall (i,j in 1..m where i > j) 
    (differences[j + ((i-1)*(i-2)) div 2] = mark[i] - mark[j]); 
constraint forall (i,j in 1..m where i > j) 
    (differences[j + ((i-1)*(i-2)) div 2] >= (floor(int2float((i-j)*(i-j+1))/int2float(2)))); 
constraint differences[1] < differences[(m*(m-1)) div 2]; 

solve :: int_search(mark, input_order, indomain, complete) 
    minimize mark[m]; 

output ["golomb ", show(mark), "\n"]; 
+0

我对这个答案仍然很感兴趣。我将再次使用这个非常公式来进行另一个项目。你可以看看代码吗?谢谢。 – Agostino 2015-01-12 21:29:33

+1

对不起,长时间回复。看看我的编辑。关键点:用T(i,j)= i +(j-1)*(j-2))/ 2来映射{(i,j)|。 1 <= i j(下三角形),但是我的代码使用元组j> i(上三角形)。无论哪种方式将工作,但你问上部三角问题,所以我给了代码。如果你想要下三角的代码,随时问。 – 2015-01-13 21:53:00

+1

删除列表理解。 var 0..n:differences = [mark [j] - mark [i] |的array [int]而不是var [int] i in 1..m,j in i + 1..m];''你只需要var 0..n的差阵[1 ..(m *(m-1))div 2]'。 – 2015-01-13 23:59:49