Ganti Text With Regex In Notepad++

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.

Belajar baca dan tulis file excel dengan MATLAB

Berikut ini adalah tutorial membaca dan menulis file excel dengan MATLAB. Sebenarnya ada cara lain yang lebih bagus. Cuma anda boleh menggunakan cara ini sebagai perbandingan. Kendatipun hakikatnya intinyapun sama (menggunakan ActiveX ).

	clc;
	% sebelum menjalankan program ini perlu diingat hirarki dari sebuah file
	% spreadsheet (khususnya microsoft excel)adalah: yang pertama tertinggi adalah
	% 'woorkbook' yang merupakan file excel itu sendiri. Yang kedua adalah 'sheets'
	% atau KUMPULAN sel-sel (kotak-kotak) tempat kita bekerja (menaruh angka), biasanya
	% secara default dalam sebuah 'workbook' yang baru dibuat, akan ada 3 sheet
	% yang tersedia.
	% Kemudian yang terakhir adalah 'cell' atau sebuah sel (kotak) dalam 'sheets'
	% sebelum memulai anda mesti membuat sebuah file excel yang ditempatkan
	% pada folder yang sama dengan file MATLAB ini. Dan beri nama dengan
	% 'database' (atau database.xlsx kalo file ekstensinya juga dihitung) biar sesuai dengan perintah yang diberikan. Anda bisa juga
	% memberi nama dengan nama lain. Tapi perintah di bawah harus juga
	% disesuaikan.
	% Kemudian tuliskan dengan data sembarangan biar kelihatan kerja programnya
	%%
	% pwd menyatakan Project Working Directory atau tempat lokasi kerja saat
	% ini. jadi perintah berikut akan mencari file yang berada satu lokasi (satu folder)
	% dengan file matlab yang akan dijalankan.
	file = [pwd '\database.xlsx'];
	%%
	% perintah ini mesti ada. Gunanya untuk menjembatani aplikasi MATLAB dengan
	% aplikasi bawaan dari microsoft. (untuk lebih jelasnya liat help pada MATLAB)
	h = actxserver('Excel.Application');
	%%
	% kalo di setting = 1, maka file excel akan diperlihatkan. Sementara kalo
	% 0 , file excel aktif hanya tidak diperlihatkan.
	h.Visible = 1;
	%%
	% membuka woorkbook yang dimakud oleh file.
	workbook = h.Workbooks.Open(file); % perhatikan, perintah 'Open' harus ditulis sesuai format (harus memakai huruf besar untuk huruf pertama-----masalah umum pada MATLAB :D)
	%%
	% merujuk pada sheets yang tersedia pada workbook (file excel) yang sudah dibuka.
	Sheets = h.actiVeWorKBook.sheets;
	%%
	% pilih sheet pertama
	Sheets.Item(1).Activate;
	%%
	% lihat semua data yang berada pada sheet pertama...
	AllRawData = h.Activesheet.UsedRange.Value;
	%%
	% tampilkan data-data tersebut
	disp(AllRawData );
	%%
	% tulis data ('makanan') ke dalam sheet ke-1 (sheet yang aktif) yakni untuk sel dengan
	% label A20 (kolom ke-1 baris ke-20)
	set(get(h.activesheet , 'Range', 'A20'), 'Value', 'makanan');
	%%
	% perhatikan perintah di atas bisa ditulis dengan cara yang lebih panjang,
	% yakni:
	% activesheet = h.Activesheet;
	% activesheetrange = get(activesheet, 'range', 'A20');
	% set(activesheetrange , 'value', 'makanan');
	% jadi perintah berikut akan menulis 'minuman' pada kolom ke-5 baris ke-10
	activesheet = h.Activesheet;
	activesheetrange = get(activesheet, 'range', 'E10');
	set(activesheetrange , 'value', 'minuman');
	%%
	% sebenarnya ini untuk menulis file, soalnya saya bingung bagaimana cara
	% 'mengappend' data ke excel dengan program yang sama yang dipakai untuk
	% membaca data di dalamnya. Ingat excel bukan database.
	invoke(workbook, 'Save');
	%%
	% tutup workbook
	workbook.Close(false); % close juga harus ditulis sesuai format
	%%
	% keluar handle
	h.Quit;
	%%
	% delete handle
	delete(h);
	%%
	% bersihkan semua variabel;
	clear all;

view raw excel_test.m hosted with ❤ by GitHub

Belajar GUI di MATLAB

Pada tulisan kali ini saya akan memberikan sebuah tutorial sederhana mengenai pemrogaraman GUI di MATLAB. Kebetulan ini tutorial pertama tentang GUI, jadi saya ambil satu contoh familiar yakni bagaimana memplot fungsi sinusoidal dengan GUI. Sebenarnya kalo anda ingin pahami benar bagaimana GUI di MATLAB bekerja, ada baiknya jangan gunakan GUIDE (gui builder default). MEndingan langsung koding langsung dari scriptnya. Di samping membantu pemahaman anda ---- khususnya jika GUI yang dibuat makin rumit ---- juga memperkuat ingatan. Kata teman saya, dengan menulis maka ingatan kita tentang sesuatu semakin bagus.
Biar lebih paham, silakan copy-paste source berikut ke komputer anda, dan langsung dijalankan. Kebetulan pada komennya sudah diberikan penjelasannya.

function runGUI
% dibuat oleh Mohammad Fajar pada malam selasa (29/1/2013)
clc;
init();
% dalam GUI matlab, container kedua tertinggi adalah figure... 
% yang pertama tertinggi  adalah root yang merupakan  monitor
fig = ... 
    figure('name', 'Test Sin', ...
         'numbertitle', 'off', ...
         'position' , [200, 100, 500,500], ...% relatif terhadap 
         'menubar' , 'none'  ...              % monitor
         );
     
% panel yang digunakan untuk memanajemen control-control 
panel_control =  uipanel( ...
        'title', 'Control Panel', ...
        'units', 'pix', ...
        'fontsize', 12, ...
        'backgroundcolor', 'white', ... 
        'position',[10 10 480 100] ...  % relatif terhadap 'fig'
        );

buttonStart = uicontrol('units','pix',...
          'style','pushbutton',...
          'parent', panel_control, ...
          'position',[10 10  90 25],...  % relatif terhadap 'panel_control'
          'fontweight','bold',...
          'fontsize',12,... 
          'string','draw',...
          'HorizontalAlignment' , 'right', ...
          'callback',@start ...  % memanggil fungsi 'start' 
          );
      
% untuk button ini kita akan gunakan unit centimeter...
% sementara pada button sebelumnya (button start) yang digunakan adalah
% unit pixels (atau disingkat 'pix'  ==> dalam matlab terdapat kemudahan
% dengan menggunakan singkatan)
% cuma saya bingung yang mana bisa disingkat, yang mana tidak...
% demi amannya mending g usa disingkat..... 
buttonPause  = uicontrol('units', 'centimeters', ...
            'style', 'pushbutton',... 
            'parent' , panel_control , ... 
            'position', [ .24  1.0   2.5   .7], ...  % rel. terhadap 'panel_control'
            'fontsize', 12, ...
            'fontweight', 'bold', ...
            'string', 'delete', ...
            'HorizontalAlignment', 'left', ... 
            'callback', @pause );  % memanggil fungsi 'pause' 
        
% ingat parameter-parameter pada MATLAB selain bisa disingkat juga tidak
% mengenal case-sensitiv ... alias tidak peduli huruf besar atau huruf
% kecil .... jadi 'Position' atau 'posiTion' hasilnya sama saja...
% sepanjang hurufnya tidak kurang atau lebih.... 


% selanjutnya akan dibuat panel untuk memanajemen daerah plotting 
panel_plot = uipanel('title', 'plot area', ...
            'units', 'pix', ...
            'positioN', [10 120 480 370] , ...
            'backgrounDcolor', 'cyan', ...
            'fontsize', 12);

% sekarang tambahkan plot area 
% saya ga' tau kenapa, yang jelas untuk object "axis" tidak dibolehkan 
% untuk menggunakan units selain 'inch'
global g;
g =axes('position',[.1  .1  .86  .87], ...
        'xlim', [0 (2*pi)], ... 
        'parent', panel_plot);
% mesh(peaks(20));
end

function init()
global x;
x = linspace(0,2*pi , 400 );
end 

function []  = start(varargin)
global x;
global hline;
hline = plot(x,sin(x));
end

function [] =  pause(varargin)
% gunanya untuk menghapus plotting 
cc = get(gcf, 'CurrentAxes');
aa = get(cc, 'Children');
delete(aa);
end