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

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

D := {x : x[i] ϵ C[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

Обозначение и «подготовительная» часть кода связаны с мое приложение.


set cartesian-product stata
person Frank    schedule 28.01.2015    source источник

Ответы (2)


arrow_upward
2
arrow_downward

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

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, имеют разную длину.

person BeingQuisitive    schedule 29.01.2015
comment
Спасибо! Я принимаю ваше элегантное решение, но также публикую свое собственное, так как подозреваю, что оно немного быстрее. - person Frank; 03.02.2015

arrow_upward
0
arrow_downward

Каждый столбец результата должен быть mm_expanded только дважды. Должно быть быстрее построить столбцы непосредственно таким образом (а не за 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, поэтому я не подтвердил свои подозрения по поводу прироста эффективности.

person Frank    schedule 02.02.2015