Vektorisasi loop molekular dinamik di MATLAB

Pada banyak kasus tidak selamanya dokumentasi resmi itu bisa dipercaya. Atau secara umum tidak semua buku yang kita baca di bangku kuliah itu bisa dipercaya. Contohnya waktu adanya sindrom HIV AIDS tahun 80-an. Semua orang menduga bahwa penyakit ini disebabkan oleh hubungan sesama jenis dan ini menjadi asumsi standar dalam prosedur  pemeriksaan dokter pada dekade tersebut. Akibatnya apa, banyak orang yang sebenarnya tertular penyakit tersebut, akan tetapi karena  semua dokter pada saat itu rujukannya pada buku yang menyatakan HIV itu hanya tertular melalui hubungan sesama jenis, tidak bisa mendapatkan penanganan yang sesuai. Nanti baru seorang Magic Johnson yang notabene orang berduit yang mengeluhkan gejala tersebut, barulah dokter-dokter berfikir di luar kotak. Padahal dengan diagnosis singkatpun seharusnya bisa disimpulakan: bagaimana mungkin seorang pebasket homo? Ajaibnya Magic Johnson mampu membayar untuk tes lebih lanjut.

Hal ini berlaku pula dalam dunia pemrograman. Dalam dokumentasi resmi MATLAB yang saya ambil di link http://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html disebutkan bahwa agar program MATLAB bisa berjalan dengan lebih cepat maka salah satu cara yang bisa digunakan adalah menggunakan pre-alocating array dan menghindari penggunaan looping:
http://stackoverflow.com/questions/12074373/matlab-vectorization-how-to-avoid-this-for-loop
http://stackoverflow.com/questions/16214891/how-to-avoid-a-for-loop-in-matlab-when-performing-a-operation-on-each-row-in-a
Akan tetapi baru-baru ini saya mencoba mengetes dugaan tersebut dalam simulasi molekular dinamik. Dan ternyata hasilnya tidak sesuai yang ada pada buku. Rupanya pengembangan MATLAB akhir-akhir ini cenderung melakukan optimisasi pada looping ketimbang prosedur standar (yakni operasi matriks) yang merupakan brand dari MATLAB.

Simulasi molekular dinamik sendiri saya translate dari sumber di
http://www.personal.psu.edu/auk183/MolDynamics/Molecular%20Dynamics%20Simulations.html
Yang aslinya dikembangkan dengan bahasa java. Setelah melakukan debugging yang susah payah, akhirnya saya berhasil mengembangkan versi MATLAB nya: menggunakan teknik vektorisasi yang menjadi brand MATLAB ketimbang menggunakan looping yang katanya  mengurangi performa. Berikut merupakan kode hasil vektorisasi tersebut

	% Belajar molekular dinamik
	function moldi1()
	clc;
	clear all;
	N = 500 ;
	f = figure('units', 'pix', 'menubar', 'none', 'position', [200, 200, 600 ,500]);
	uicontrol('units', 'pix', 'style', 'pushbutton','callback',@stop ,'position', ...
	[60, 10, 100, 30], 'string', 'STOP');
	axis([0 10 0 10]);
	axis manual;
	boxWidth = 100;
	boxHeight = 100;
	dt = 0.025;
	gravity = 0.00;
	wallStiffness = 50.0;
	forceCutoff = 3.0;
	Temperature = 0.5 ;
	rectangle('position', [0,0, boxWidth, boxHeight], 'curvature', [1 1], ...
	'facecolor' , 'r');
	x = zeros(N,1);
	y = zeros(N,1);
	vx = zeros(N,1);
	vy = zeros(N,1);
	ax = zeros(N,1);
	ay = zeros(N,1);
	dtOver2 = 0.5 * dt;
	dtSquaredOver2 = 0.5 * dt * dt;
	forceCutoff2 = forceCutoff^2 ;
	boxWidthMinusHalf = boxWidth - 0.5;
	boxHeightMinusHalf = boxHeight - 0.5;
	volumePerMolecule = (boxWidth * boxHeight) /N;
	averageSpacing = sqrt( volumePerMolecule);
	rows = ceil(boxHeight / averageSpacing);
	columns = ceil(boxWidth / averageSpacing);
	averageSpacing = min(boxHeight/rows , boxWidth/columns );
	if averageSpacing < 1.0 ,
	averageSpacing = 1.0;
	end
	xx = averageSpacing /2.0 ;
	yy = xx ;
	for i=1:N,
	x(i) = xx;
	y(i) = yy;
	vx(i) = rand(1,1) - 0.5;
	vy(i) = rand(1,1) - 0.5;
	ax(i) = 0.0;
	ay(i) = 0.0;
	x = x + averageSpacing ;
	if x > (boxWidth - 0.5),
	x = averageSpacing/2.0;
	y = y+ averageSpacing ;
	if y > (boxHeight - 0.5),
	break;
	end
	end
	end
	nv = 2*N;
	ng = nv - 4;
	ek = sum(vx.^2) + sum(vy.^2);
	vs = sqrt(1.0 * ek *nv /(ng*Temperature));
	vx = vx ./ vs ;
	vy = vy ./ vs;
	iter = 0 ;
	run = 1;
	time= 0;
	while(run),
	tic;
	iter = iter + 1;
	a = (vx(:,end) .* dt) + (ax(:,end) .* dtSquaredOver2);
	b = (vy(:,end) .* dt) + (ay(:,end) .* dtSquaredOver2);
	x = x + a;
	y = y + b;
	plot(x, y,'b.','MarkerSize',12);
	drawnow;
	% =======================================================================
	ax(:) = 0 ;
	index1 =x < 0.5;
	ax(index1) = wallStiffness *(0.5 - x(index1));
	index2 = x > boxWidthMinusHalf;
	ax(index2) = wallStiffness * (boxWidthMinusHalf - x(index2)) ;
	ay(:) = gravity;
	index3 = y < 0.5;
	ay(index3) = (wallStiffness * (0.5 -y(index3)));
	index4 = y > boxHeightMinusHalf;
	ay(index4) = (wallStiffness *(boxHeightMinusHalf) - y(index4));
	for i=1:N,
	dx=zeros(i-1,1);
	dy=zeros(i-1,1);
	rinv = zeros(i-1,1);
	attract = zeros(i-1,1);
	repel = zeros(i-1, 1);
	fOverR = zeros(i-1,1);
	fx = zeros(i-1,1);
	fy = zeros(i-1,1);
	dx(:,1)= x(i) - x(1:i-1,1);
	dx2 = dx.^2 ;
	dy(:,1) = y(i) - y(1:i-1,1) ;
	dy2= dy.^2;
	r = dx2 + dy2;
	index = ((dx2 < forceCutoff2) &(dy2 < forceCutoff2) &...
	(r < forceCutoff2)) ;
	rinv(:) = 1.0 ./ r;
	attract(:) = rinv.^3;
	repel(:) = attract.^2;
	fOverR(:) = 24.0.*((1.3.*repel) - attract).* rinv ;
	fx(index) = fOverR(index).* dx(index);
	fy(index) = fOverR(index).* dy(index);
	ax(1:i-1) = ax(1:i-1)-fx(:);
	ay(1:i-1) = ay(1:i-1)-fy(:);
	ax(i) = ax(i)+ sum(fx(:));
	ay(i) = ay(i)+ sum(fy(:));
	end
	c = ax(:).* dtOver2;
	d = ay(:).* dtOver2;
	vx(:) = vx(:) + c ;
	vy(:) = vy(:) + d ;
	time = time + toc;
	pause(0.1);
	end
	% =====================================================================
	function stop(varargin)
	disp('STOP looping');
	run = 0;
	disp(['rata-rata waktu ' num2str(time/iter)]);
	close(f);
	end
	end

view raw moldi1.m hosted with ❤ by GitHub

Adapun kode hasil translate secara 'mentah' adalah sebagai berikut

	function moldi()
	clc;
	clear all;
	N = 500 ;
	f = figure('units', 'pix', 'menubar', 'none', 'position', [200, 200, 600 ,500]);
	uicontrol('units', 'pix', 'style', 'pushbutton','callback',@stop ,'position', ...
	[60, 10, 100, 30], 'string', 'STOP');
	axis([0 10 0 10]);
	axis manual;
	boxWidth = 100;
	boxHeight = 100;
	dt = 0.025;
	gravity = 0.00;
	wallStiffness = 50.0;
	forceCutoff = 3.0;
	Temperature = 0.5 ;
	rectangle('position', [0,0, boxWidth, boxHeight], 'curvature', [1 1], ...
	'facecolor' , 'r');
	x = zeros(N,1);
	y = zeros(N,1);
	vx = zeros(N,1);
	vy = zeros(N,1);
	ax = zeros(N,1);
	ay = zeros(N,1);
	dtOver2 = 0.5 * dt;
	dtSquaredOver2 = 0.5 * dt * dt;
	forceCutoff2 = forceCutoff^2 ;
	boxWidthMinusHalf = boxWidth - 0.5;
	boxHeightMinusHalf = boxHeight - 0.5;
	volumePerMolecule = (boxWidth * boxHeight) /N;
	averageSpacing = sqrt( volumePerMolecule);
	rows = ceil(boxHeight / averageSpacing);
	columns = ceil(boxWidth / averageSpacing);
	averageSpacing = min(boxHeight/rows , boxWidth/columns );
	if averageSpacing < 1.0 ,
	averageSpacing = 1.0;
	end
	xx = averageSpacing /2.0 ;
	yy = xx ;
	for i=1:N,
	x(i) = xx;
	y(i) = yy;
	vx(i) = rand(1,1) - 0.5;
	vy(i) = rand(1,1) - 0.5;
	ax(i) = 0.0;
	ay(i) = 0.0;
	x = x + averageSpacing ;
	if x > (boxWidth - 0.5),
	x = averageSpacing/2.0;
	y = y+ averageSpacing ;
	if y > (boxHeight - 0.5),
	break;
	end
	end
	end
	nv = 2*N;
	ng = nv - 4;
	ek = sum(vx.^2) + sum(vy.^2);
	vs = sqrt(1.0 * ek *nv /(ng*Temperature));
	vx = vx ./ vs ;
	vy = vy ./ vs;
	iter = 0 ;
	run = 1;
	time= 0;
	while(run),
	tic;
	iter = iter + 1;
	a = (vx(:,end) .* dt) + (ax(:,end) .* dtSquaredOver2);
	b = (vy(:,end) .* dt) + (ay(:,end) .* dtSquaredOver2);
	x = x + a;
	y = y + b;
	plot(x, y,'b.','MarkerSize',12);
	drawnow;
	% =====================================================================
	for i=1:N,
	if x(i) < 0.5
	ax(i) = wallStiffness *(0.5 - x(i));
	else
	if x(i) > boxWidthMinusHalf,
	ax(i) = wallStiffness *(boxWidthMinusHalf - x(i));
	else
	ax(i) = 0.0;
	end
	end
	if y(i) < 0.5
	ay(i) = wallStiffness*(0.5 - y(i))+ gravity ;
	else
	if y(i) > boxHeightMinusHalf,
	ay(i) = (wallStiffness *(boxHeightMinusHalf - y(i))) + gravity;
	else
	ay(i) = gravity ;
	end
	end
	end
	for i =1:N
	for j=1:i-1
	dx = x(i) - x(j);
	dx2 = dx^2;
	if dx2 < forceCutoff2
	dy = y(i)- y(j);
	dy2 = dy^2 ;
	if dy2 < forceCutoff2
	rSquared = dx2 + dy2 ;
	if rSquared < forceCutoff2
	rSquaredInv = 1.0 / rSquared ;
	attract = rSquaredInv^3;
	repel = attract^2;
	fOverR = 24.0 *((1.3*repel) - attract ) * rSquaredInv ;
	fx = fOverR *dx ;
	fy = fOverR *dy ;
	ax(i) = ax(i) + fx ;
	ay(i) = ay(i) + fy ;
	ax(j) = ax(j) - fx ;
	ay(j) = ay(j) - fy ;
	end
	end
	end
	end
	end
	% =======================================================================
	c = ax(:).* dtOver2;
	d = ay(:).* dtOver2;
	vx(:) = vx(:) + c ;
	vy(:) = vy(:) + d ;
	time = time + toc;
	pause(0.1);
	end
	function stop(varargin)
	disp('STOP looping');
	run = 0;
	disp(['rata-rata waktu ' num2str(time/iter)]);
	close(f);
	end
	end

view raw moldi.m hosted with ❤ by GitHub

Setelah dieksekusi terlihat bahwa hasil penggunaan loop bukannya memperlambat performa, malah mempercepat performa dari program MATLAB sendiri. Rata-rata waktu eksekusi looping berada di sekitar 0.03 detik, sementara hasil vektorisasi ternyata membutuhkan waktu 3 kali lebih lama yani berada di sekitar 0.1 detik.