Matlab: VECTORISATION

[SnakePlinsky]

Tramite search ho visto che non c'è nessuna discussione sulla "vettorizzazzione", quindi provo ad iniziarne una io (anche perchè avrei bisogno di suggerimenti).

La "vectorization" è un modo per rendere più veloci i codici in Matlab, evitando di utilizzare cicli e rimpiazzandoli con operazioni direttamente sui vettori.

Alcuni link per capire di quello di cui stò parlando:
http://www.eng.cam.ac.uk/help/tpl/progr ... ricks.html (breve guida)
http://www.mathworks.com/support/tech-n ... 1109.shtml ( Mathworks Code Vectorization Guide)

Alcuni esempi (tratti dal primo link, in rosso la soluzione "vectorized")
# Vectorise the following, where data is 200x400

for i = 1:200
   for j = 1:400
     if data(i,j) > 0
           data(i,j) = 0;
     end
   end
end

data(data > 0) = 0;

for i = 1:100
    for j = 1:100
      r(i,j) = sqrt(i^2+j^2);
    end
  end

[i,j]=meshgrid(1:100,1:100);
r = sqrt(i.^2+j.^2);

n=1000;
x(1)=1;
for j=1:n-1,
   x(j+1) = x(j) + n - j;
end

n=1000;
x(1)=1;
j=1:n-1;
x(j+1) = n - j;
cumsum(x);


Per chi ha avuto la pazienza di seguirmi fin qui, propongo il seguente problema:

Un primo esempio:
Ho due vettori A e B, entrambi Nx1.

Voglio ottenere un vettore dummy, che abbia valore 1 quando il valore nella cella A(i,1) > B(i,1), -1 quando il valore nella cella A(i,1) < B(i,1) e 0 quando il valore nella cella A(i,1) = B(i,1).

Modo classico:

for i=1:N
        if A(i,1) > B(i,1)
        dummy(i,1)=1;
        elseif A(i,1) < B(i,1)
        dummy(i,1)=-1;
        elseif A(i,1) > B(i,1)
        dummy(i,1)=0;
        end
end

Tramite vectorisation:

dummy(A>B)=1;
dummy(A dummy(A==B)=1;

Il riparmio in tempo di calcolo è evidente.

Fin qui ci siamo, ora però mi serve una matrice che mi sommi il numero di volte che la dummy varia:

es dummy' = [ 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 -1 -1 -1 1 1 1 1]

abbiamo in totale 4 variazioni (evidenziate in rosso).

La domanda è : come faccio ad ottenere nel caso generale il numero di variazioni senza ricorrere acicli for e if? Ogni tipo di aiuto è gradito.

[SnakePlinsky]

"SnakePlinsky":
.....

Fin qui ci siamo, ora però mi serve una matrice che mi sommi il numero di volte che la dummy varia:

es dummy' = [ 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 -1 -1 -1 1 1 1 1]

abbiamo in totale 4 variazioni (evidenziate in rosso).

La domanda è : come faccio ad ottenere nel caso generale il numero di variazioni senza ricorrere acicli for e if? Ogni tipo di aiuto è gradito.

Detto fatto :

a=1:N-1;
p(a+1)=dummy(a+1)-dummy(a);
p(p~=0)=1;
Numero_di_variazioni=sum(p)

Se qualcuno trova un modo alternativo o più breve deve ovviamente postarlo

[amel3]

Si potrebbe anche fare:

sum(dummy(2:end)~=dummy(1:end-1))

direi che funziona (o almeno credo...), comunque è praticamente equivalente a quello che hai scritto tu.

[SnakePlinsky]

"Sergio":Se ho capito cosa vuoi fare...

octave:1> n = 10;
octave:2> A = rand(n,1);
octave:3> B = rand(n,1);
octave:4> cfr = zeros(n, 1);
octave:5> cfrg = cfr; cfrg(A>B)=1
cfrg =

   0
   1
   1
   1
   0
   0
   0
   1
   0
   1

octave:6> sum(cfrg != cfr)
ans =  5
octave:7> cfrl = cfrg; cfrl(A<B)=-1
cfrl =

  -1
   1
   1
   1
  -1
  -1
  -1
   1
  -1
   1

octave:8> sum(cfrl != cfrg)
ans =  5

Ciao Sergio, ho visto ora il tuo messaggio. Si intendevo proprio un modo per vedere quante volte una dummy (immaginala come una serie temporale) varia, ovvero il numero di mutamenti nel tempo.
Io ho risolto nel modo che puoi vedere... ( per il mio codice basta e avanza la soluzione che ho adottato) ; se riesci a trovare un modo più parsimonioso tanto di cappello.

Comunque la discussione è aperta al problema della "vectorisation(*)" più in generale, visto che non c'era l'argomento: chi vuole postare suoi problemi (codici da velocizzare) faccia pure.

Un saluto.

[SnakePlinsky]

"amel":Si potrebbe anche fare:

sum(dummy(2:end)~=dummy(1:end-1))

direi che funziona (o almeno credo...), comunque è praticamente equivalente a quello che hai scritto tu.

No anzi, è molto più elegante del mio.


Visto che accettate la sfida ve ne posto uno veramente difficile:

Abbiamo una matrice RL di apiezza Nxn

Vogliamo ottenere un vettore riga 1xN in cui in ogni cella in colonna b abbiamo il valore : sum(RL(:,b))/log(1-(mean(RL(:,b)(RL(:,b)
con un ciclo for

for b=1:n
        P_f(1,b)=sum(RL(:,b))/log(1-(mean(RL(:,b)(RL(:,b)<quantile(RL(:,b),alpha_VAR)))));
end

Ammetto che non sò se sia risolvibile, in quanto utilizzo già all'interno dell'espressione un formula "vectorized", per ora quindi nel mio codice uso il ciclo postato.

A chi riesce a eliminarmi il ciclo for gli offro da bere (vi avverto, non sforzatevi troppo, non sò nemmeno se la cosa sia possibile, io ci ho rinunciato)

(*) per chi volesse sapere cosa diavolo significa quell'espressione: http://www.smartfolio.com/theory/detail ... ices/c/#95 (STARR Ratio)

[amel3]

Ma se tu fai così non è la stessa cosa?

P_f(1,1:n)=sum(RL(:,1:n))./log(1-(mean(RL(:,1:n)(RL(:,1:n)<quantile(RL(:,1:n),alpha_VAR)))));

No, eh?

[SnakePlinsky]

"amel":Ma se tu fai così non è la stessa cosa?

P_f(1,1:n)=sum(RL(:,1:n))./log(1-(mean(RL(:,1:n)(RL(:,1:n)<quantile(RL(:,1:n),alpha_VAR)))));

No, eh?

Ci avevo provato, ma mi da questo error message:

??? Error: ()-indexing must appear last in an index expression.

Più in generale ho notato che quando cerco di "vettorizzare" con sub-matrici tramite indexing mi da lo stesso error

PEr esempio se voglio rendere pari a 1 le celle nella colonna 1 della matrice magic 5x5 a seconda che i valori delle corrispondenti celle sia maggiore di 10 nella colonna 2

A=magic(5)

A(:,1)(A(:,2)>10)=1

??? Error: ()-indexing must appear last in an index expression.

mi da lo stesso errore.

[SnakePlinsky]

Devo moltiplicare ogni colonna j di una matriceA mxn per il valore j di un vettore b 1xn

con ciclo for

for j = 1:n
     X(:,j)=A(:,j).*b(1,j);
end

Soluzione vectorised:

X= A.*repmat(b,m,1)

Alla prossima

[_luca.barletta]

Un esercizio:
sono dati un vettore di valori x[1:n], un vettore in cui sono specificate le occorrenze f[1:n],
bisogna costruire un nuovo vettore in cui ogni elemento di x, x(i) per i=1...n, compaia esattamente f(i) volte.
Ad es. $x={1,2,3}$ e $f={2,4,1}$, allora il risultato è $r={1,1,2,2,2,2,3}$.
Provate a farlo con il minor numero di cicli in matlab o, ancora meglio, senza cicli.

[SnakePlinsky]

"luca.barletta":Un esercizio:
sono dati un vettore di valori x[1:n], un vettore in cui sono specificate le occorrenze f[1:n],
bisogna costruire un nuovo vettore in cui ogni elemento di x, x(i) per i=1...n, compaia esattamente f(i) volte.
Ad es. $x={1,2,3}$ e $f={2,4,1}$, allora il risultato è $r={1,1,2,2,2,2,3}$.
Provate a farlo con il minor numero di cicli in matlab o, ancora meglio, senza cicli.

Coi vettori come da te definiti, il meglio che mi è venuto in mente, ha un 1 ciclo for

r=repmat(x(1,1),1,f(1,1));
for a=2:length(x)
    r=cat(2,r,(repmat(x(1,a),1,f(1,a))));
end

ma io sono un principiante, penso si possa fare di meglio

[_luca.barletta]

"SnakePlinsky":
ma io sono un principiante, penso si possa fare di meglio

sarà la fortuna del principiante, ma pare che lavori un pelino meglio del codice che avevo in mente
Anch'io avevo pensato ad un ciclo for, e non credo si possa fare di meglio (sarò felice di essere smentito).