2014-10-16 136 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