PDA

Visualizza la versione completa : [c++ vs matlab] ottimizzare codice c++ per configgere matlab


mamo139
16-11-2011, 02:52
Ciao a tutti!
spero che l'argomento possa essere di interesse a qualcuno visto che per la mia testolina da programmatore della domenica non è banale :zizi:

allora sto scrivendo un software che attraverso simulazione di montecarlo trovi il prezzo di un certo derivato esotico (precisamente una opzione magrabe ma non credo vi sia di interesse)

l'ho fatto sia in matlab che in c++, sicuro che il c++ avrebbe stracciato matlab di tanto... e invece mi ritrovo con un codice matlab piu veloce :confused:

eccovi le due funzioni incriminate

c++


//using normal distribution object from Technical Report 1
typedef std::tr1::ranlux64_base_01 Myeng;
typedef std::tr1::normal_distribution<double> ndistr;

double Margrabe_option_MC_price(double r, double t, double s1, double s2, double q1, double q2,
double std1, double std2, double corr, double iterations){
//variables and objects initialization
register long x;
double price = 0; //the final price will be stored here
double payoff;
double z1, z2; //to store generated random numbers
double matr21, matr22;
double sigma1, sigma2;
double part1, part2, sqrtt;

clock_t start, stop;

Myeng eng;
ndistr stdnorm(0, 1); //creating the standard normal distribution generator

//seeding the random number generator
eng.seed((unsigned long) time(NULL));

//calculating sigmas
matr21 = corr * std2;
matr22 = sqrt(1-corr*corr) * std2;

//precalculating fixed part of drive formula
part1 = exp((r - q1 - 0.5*std1*std1)*t);
part2 = exp((r - q2 - 0.5*std2*std2)*t);
sqrtt = sqrt(t);

for(x=0; x < iterations ;x++){
//generating random numbers
z1 = stdnorm(eng);
z2 = stdnorm(eng);

//creating sigmas
sigma1 = z1 * std1;
sigma2 = matr21 * z1 + matr22 * z2;

//calculating payoff
payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
price += (payoff > 0) ? payoff : 0;
}
price /= iterations;
price *= exp((-r)*t);

return price;
}


matlab


function [p,ci] = ExchangeMC(V0,U0,sigmaV,sigmaU,rho,T,r,NRepl)
eps1 = randn(1,NRepl);
eps2 = rho*eps1 + sqrt(1-rho^2)*randn(1,NRepl);
VT = V0*exp((r - 0.5*sigmaV^2)*T + sigmaV*sqrt(T)*eps1);
UT = U0*exp((r - 0.5*sigmaU^2)*T + sigmaU*sqrt(T)*eps2);
DiscPayoff = exp(-r*T)*max(VT-UT, 0);
[p,s,ci] = normfit(DiscPayoff);


ebbene matlab per fare 100.000 generazioni ci mette 24 centesimi di secondo, contro i 38 del codice c++.

Considerando che ho sempre pensato il c/c++ essere i linguaggi di alto livello piu performanti che esistono (correggetemi se sbaglio) devo aver sbagliato io qualcosa.

come rendereste voi il codice c piu efficiente?

da ignorante mi è venuto da pensare che il problema potrebbe essere nella funziona di generazione di numeri casuali di una normale standard!!
ho controllato e matlab lo risolve con questi calcoli
eps1 = randn(1,NRepl);
eps2 = rho*eps1 + sqrt(1-rho^2)*randn(1,NRepl);
in un elegante 8 centesimi

mentre la funzione che ho trovato io nelle librerie del boost che da quel che ho capito riscreano il pacchetto aggiuntivo Technical Report 1
for(x=0; x < iterations ;x++){
//generating random numbers
z1 = stdnorm(eng);
z2 = stdnorm(eng);
}
se ne esce con un goffo 26 centesimi di secondo

come posso migliorare questo aspetto?
inoltre avete altre osservazioni sui codici??

ci tengo particolarmente a rendere il c++ piu performante... :)

grazie
Marco

MItaly
16-11-2011, 09:35
Con che compilatore/opzioni lo stai compilando? Visto che c'è molto calcolo in virgola mobile potrebbe beneficiare delle SSE...
Provato a cercare i colli di bottiglia con un profiler?

Inoltre, sei sicuro della misurazione del tempo? Una singola prova "a freddo" del codice non è significativa (cache fredda, branch predictor ancora non "istruito", ...) , dovresti ripetere il codice un tot di volte (ottenendo un tempo più grande e quindi meno affetto da errori di misura) e vedere quanto ci mette in tutto, eventualmente ripetendo la prova è confrontando poi media e deviazione standard con quella che viene da MatLab.

mamo139
16-11-2011, 12:13
Originariamente inviato da MItaly
Con che compilatore/opzioni lo stai compilando? Visto che c'è molto calcolo in virgola mobile potrebbe beneficiare delle SSE...
Provato a cercare i colli di bottiglia con un profiler?

Inoltre, sei sicuro della misurazione del tempo? Una singola prova "a freddo" del codice non è significativa (cache fredda, branch predictor ancora non "istruito", ...) , dovresti ripetere il codice un tot di volte (ottenendo un tempo più grande e quindi meno affetto da errori di misura) e vedere quanto ci mette in tutto, eventualmente ripetendo la prova è confrontando poi media e deviazione standard con quella che viene da MatLab.

Sto compilando con vc++ 2010 express.

Purtroppo non sono molto esperto e devo chiederti qualche delucidazione!
Non trovo nulla riguardante SSE, sai per caso come attivare questa feature in vc++?

Ehm non ho mai usato un profiler in vita mia... consigli?

il tempo lo rilevo con la funzione clock()


start = clock();
for(x=0;x<trials;x++){
price = Margrabe_option_MC_price(0.02, 0.5, 20, 15, 0, 0,
0.3, 0.2, 0.5, 1000000);
}
stop = clock();

printf("Elapsed time: %f\nAverage elapsed time per MC: %f\n\n", (double)(stop-start)/CLOCKS_PER_SEC, (double)(stop-start)/CLOCKS_PER_SEC/trials);


va bene come misuratore?

grazie

MItaly
16-11-2011, 12:55
Partiamo dalle basi: hai compilato in modalità Release o Debug (c'è una casella in alto)? Perché in debug sono disattivate le ottimizzazioni...

mamo139
16-11-2011, 13:25
Originariamente inviato da MItaly
Partiamo dalle basi: hai compilato in modalità Release o Debug (c'è una casella in alto)? Perché in debug sono disattivate le ottimizzazioni...

Buona idea :D
Ho compilato in modalità Release!

MItaly
16-11-2011, 14:17
Ok; nelle proprietà del progetto (per quanto riguarda la configurazione Release) cerca e attiva "Ottimizzazione globale", "Ometti puntatore al frame", "Ottimizza per velocità". Da qualche parte ci sarà "set di istruzioni aggiuntivi", abilita MMX, SSE e SSE2.

mamo139
16-11-2011, 14:25
Visto che non so come usare un profiler (appena ho tempo imparerò) ho spezzettato i punti del codice per vedere dove matlab va meglio di c++.

ho separato i calcoli dalla generazione dei numberi casuali e questi sono i risulstati:

matlab


Elapsed time is 0.066459 seconds.
Elapsed time is 0.161765 seconds.


c++


Elapsed time: 0.170000
Elapsed time: 0.149000


Quindi in c++ sono gia molto piu veloce nella parte di calcoli (considera che matlab lavora gia in multithreading e il mio pc ha due core, mentre il mio codice c++ deve essere ancora modificato per far andare piu calcoli in parallelo

direi che il problema è il generatore di numeri casuali!!

eccoti il codice con cui ho registrato i tempi:


inline double Margrabe_option_MC_price(double r, double t, double s1, double s2, double q1, double q2,
double std1, double std2, double corr, long iterations){
clock_t start, stop;
start = clock();

//variables and objects initialization
register long x,y;
double price = 0; //the final price will be stored here
double payoff;
double *z1; //to store generated random numbers
double matr21, matr22;
double sigma1, sigma2;
double part1, part2, sqrtt;

//memory allocation
z1 = new double[iterations*2];

//seeding the random number generator
#ifdef RANDOM_MODE_TR1
Myeng eng;
ndistr stdnorm(0, 1); //creating the standard normal distribution generator
eng.seed((unsigned long) time(NULL));
#endif
#ifdef RANDOM_MODE_Mersenne_twister
init_genrand((int)time(NULL));
#endif
//calculating sigmas
matr21 = corr * std2;
matr22 = sqrt(1-corr*corr) * std2;

//precalculating fixed part of drive formula
part1 = exp((r - q1 - 0.5*std1*std1)*t);
part2 = exp((r - q2 - 0.5*std2*std2)*t);
sqrtt = sqrt(t);

stop = clock();
printf("Elapsed time: %f\n", (double)(stop-start)/CLOCKS_PER_SEC);
start = clock();

for(x=0; x < iterations ;x++){
y=x+x;

//generating random numbers
#ifdef RANDOM_MODE_TR1
z1[y] = stdnorm(eng);
z1[y+1] = stdnorm(eng);
#endif
#ifdef RANDOM_MODE_Mersenne_twister
z1[y] = normsinv(genrand_real3());
z1[y+1] = normsinv(genrand_real3());
#endif
}


stop = clock();
printf("Elapsed time: %f\n", (double)(stop-start)/CLOCKS_PER_SEC);
start = clock();

for(x=0; x < iterations ;x++){
y=x+x;

//creating sigmas
sigma1 = z1[y] * std1;
sigma2 = matr21 * z1[y] + matr22 * z1[y+1];

//calculating payoff
payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
if(payoff > 0) price += payoff;
}
price /= iterations;
price *= exp((-r)*t);

delete [] z1;

stop = clock();
printf("Elapsed time: %f\n", (double)(stop-start)/CLOCKS_PER_SEC);

return price;
}

mamo139
16-11-2011, 14:25
Originariamente inviato da MItaly
Ok; nelle proprietà del progetto (per quanto riguarda la configurazione Release) cerca e attiva "Ottimizzazione globale", "Ometti puntatore al frame", "Ottimizza per velocità". Da qualche parte ci sarà "set di istruzioni aggiuntivi", abilita MMX, SSE e SSE2.

ok ora provo grazie :)

MacApp
17-11-2011, 01:51
Originariamente inviato da mamo139
l'ho fatto sia in matlab che in c++, sicuro che il c++ avrebbe stracciato matlab di tanto... e invece mi ritrovo con un codice matlab piu veloce :confused:

Non escludere a priori che magari matlab stesso sia stato scritto in C++ ;-)

mamo139
20-11-2011, 16:38
Originariamente inviato da MacApp
Non escludere a priori che magari matlab stesso sia stato scritto in C++ ;-)

in realtà sono abbastanza sicuro che matlab sia scritto in c o c++
ma proprio per questo un programma in c puro dovrebbe vincere :D

alla fine ho risolto il problema e battuto matlab ottimizzato per multithreading (utilizzato su un dual core) con un codice C non ancora ottimizzato :)
semplicemente ho ricercato lo stesso algoritmo usato da matlab e l'ho implementato! in questo caso era tutta una questione di algoritmo di generazione dei numeri random piuttosto che un problema di ottimizzazione del codice!

ecco l'algoritmo utilizzato per chiunque fosse interessato
http://www.mathworks.it/company/newsletters/news_notes/clevescorner/spring01_cleve.html

avrò provato una decina di algoritmi trovati su dei paper prima di arrivare a questo e nessuno andava così bene! quindi se vi interessa l'argomento consiglio la lettura! :D

Loading