декартово произведение в Mata

Question

Чтобы построить набор векторов, мне нужно взять декартово произведение множеств С [1] .. C [D],

D: = {х: х [г] ε С [I], I = 1..d}

Пример: Если *C[1]=(5,6,7)';*C[2]=(3,5,6)';*C[3]=(1,3,5)', то некоторые элементы D являются (5,3,1), (5,3,3) ...

Я хотел бы знать: Каков наилучший способ взять декартово произведение в Мата, в общем? Я нашел неуклюжий подход для d = 3, как показано ниже.

---

Подробный пример. Этот код должен иллюстрировать то, что я пробовал, и желаемый результат. Функция mm_expand поступает от ssc install moremata.

mata 

// prep 

lo = (5,3,1)' 
hi = (7,6,5)' 
all = uniqrows((lo\hi)) 

n_cols = length(lo) 
n_vals = length(all) 

c_list = J(1,n_cols,NULL) 
c_lens = J(1,n_cols,0) 

for (i=1;i<=n_cols;i++){ 
    c_list[i] = &(select(all,all :>= lo[i] :& all :<= hi[i])) 
    c_lens[i] = length(*c_list[i]) 
} 

// question: How should I take this Cartesian product? 

grid_box = 
mm_expand(*c_list[1],c_lens[2]*c_lens[3],0,1), 
mm_expand(mm_expand(*c_list[2],c_lens[1],0,1),c_lens[3],0,0), 
mm_expand(*c_list[3],c_lens[1]*c_lens[2],0,0) 

// (just fyi) my next step 

is_decr = ! rowsum(grid_box[,1..(n_cols-1)]-grid_box[,2..n_cols] :< 0) 
select(grid_box,is_decr) 

end 

нотация и «подготовительную» часть кода связаны с my application.

Answer 1

Самый простой способ состоит в использовании рекурсии:

real matrix cart_prod(pointer vector c_list ,| real scalar curr_i){ 
    if(curr_i==.) curr_i=1 
    myret = (*c_list[curr_i]) 
    if (curr_i<length(c_list)){ 
     ret = cart_prod(c_list, curr_i+1) 
     myret = mm_expand(myret,rows(ret),1,1), mm_expand(ret, rows(myret),1,0) 
    } 
    return(myret) 
} 
cart_prod(c_list) 

Это будет работать, даже если векторы указывают на из c_list имеют разную длину.

Answer 2

Каждый столбец результата должен быть только mm_expand ed дважды. Быстро строить столбцы напрямую (а не в шагах d-k для k-го столбца, как в ответе @ BeingQuisitive).

function prod(x,| real scalar need_int){ 
    y = exp(sum(log(x))) 
    if (need_int==0) return(y) 
    return(round(y)) 
} 

function cartem(pointer vector vlist){ 

    // input: vlist should point to column vectors 

    d = length(vlist) 
    lens = J(1,d,.) 
    for (i=1;i<=d;i++){ 
     lens[i] = length (*vlist[i]) 
    } 
    tot_len = prod(lens) 
    out = J(tot_len,d,.) 

    out[,1] = mm_expand(*vlist[1],round(tot_len/lens[1]),0,1) 
    out[,d] = mm_expand(*vlist[d],round(tot_len/lens[d]),0,0) 

    if (d == 2) return(out) 

    for (i=2;i<=d-1;i++){ 
     out[,i] = mm_expand(
     mm_expand(*vlist[i],prod(lens[1..i-1]),0,0) 
     ,prod(lens[i+1..d]),0,1) 
    } 

    return(out) 
} 

Вот пример, основанный на моем прецеденте. Это показывает, что приведенный выше код производит тот же самый (желаемый) результат как @ BeingQuisitive решения которого:

c_list = (&(7\10),&(5\6\7),&(3\5\6),&(1\3\5)) 

cartem(c_list) == cart_prod(c_list) 
//^it's true 

Я не очень хорошо разбираюсь в бенчмаркинг инструментах STATA, так что я не подтвердил мои подозрения о повышении эффективности.