function RunningBlob( action, param )
%RunningBlob can be started from matlab by simply typing 'RunningBlob'
%The program is an adaptation of a java Program by Hartmut S. Loos.
%It will graph a function in the lower axis, and then in the two upper
%accessing the activity in the weight matrixes is plotted.
%The weight plots are plotted using the matlab function 'pcolor' which
%creates a psuedocolor plot of the values in the matrix.
%All variables are descibed at the beginning of the 'addGui' function.
%Andrew C. Chase
%March 21, 2002

        if nargin < 1,
                param = '';
                action = 'start';
        elseif nargin < 2,
                param = '';
        end;

        switch( action )
                case 'start'
                        addGui;
		case 'runCB'
			getValues;		
			compute;
		case 'quit'
			fig = gcf;
			quit_reply = questdlg('Really Quit?');
			if strcmp(quit_reply, 'Yes')
				close;
			end;

	end;

function addGui

	%set up data structure.
	scale = 5;
	gSize = 2 * scale + 1; %matrix size variable.
	hSize = 2 * gSize; %matrix size variable.
	
	delta = 0.7; %delta and end_time determine how long to continue
		     %looping in the 'compute' function.
	end_time = 20.0;
	
	loops = 5; %loops without updating the display.
	deltaX = 0; %size of one 'neuron'.	
	lPos = 0.0;    % self inhibition time scale - positiv dyn
	lNeg = 0.0;    % self inhibition time scale - negativ dyn
	alpha = 0.0;   % strength of the self-inhibition
	sigmas = 0.0;  % sigmas in the Gauss mask
	beta = 0.0;    % inhibitory base
	spike = 1.0;   % random spike pulse

	hIn = zeros(hSize, hSize);       % membrane potential
	selfInh = zeros(hSize, hSize);   % delayed self inhibition
	dyn = zeros(hSize, hSize);       % dynamic between hIn and selfInh
	hOut = zeros(hSize, hSize);      % cell activity
	hTransE = zeros(hSize, hSize);   % excitatory transfered signal
	hTrans = zeros(hSize, hSize);    % transfered signal: hTransE - hTransI

	ge = zeros(gSize, 1);          % excitatory connectivity kernel, 1 dim.
	mCol = zeros(gSize, 1);

	mainStruct = struct('loops', loops, 'deltaX', deltaX, 'lPos', lPos, ...
			'lNeg', lNeg, 'alpha', alpha, 'sigmas', sigmas, ...
			'beta', beta, 'spike', spike, 'hIn', hIn, ...
			'selfInh', selfInh, 'dyn', dyn, 'hOut', hOut, ...
			'hTransE', hTransE, 'hTrans', hTrans, 'ge', ge, ...
			'gSize', gSize, 'mCol', mCol, 'delta', delta, ...
			'end_time', end_time, 'hSize', hSize, 'scale', ...
			scale);
	
        height = 700; %Figure dimensions. If changed, the rest of the uicontrols
        width = 800;  %will probably need to be changed as well.

	%%%%%%%%%%%%%%%%%%Main Figure %%%%%%%%%%%%%%%%%%%%%%
        mainWindow = figure(    'ColorMap', [0.7 0.7 0.7], ...
                                'Position',[50 200 width height], ...
                                'Tag','mainWindow', ...
                                'Resize','off', ...
                                'Name','RunningBlob', ...
                                'MenuBar','none', ...
                                'ToolBar','none', ...
                                'UserData', mainStruct);

        %Set the title
        uicontrol( 'Parent', mainWindow, ...
                'Position', [(width/2 - 70) (height - 25) 140 15], ...
                'Style', 'text', ...
                'HorizontalAlignment', 'center', ...
                'String', 'RunningBlob');

        bottom = 5;
        inHeight = 20;
        %%%%%%%%%%%%%Add input fields and buttons%%%%%%%%%%%%%
        uicontrol('Parent', mainWindow, ...
                'Position', [40 bottom 50 inHeight], ...
                'Style', 'text',...
                'String', 'Beta=', ...
                'HorizontalAlignment', 'right');

        uicontrol('Parent', mainWindow, ...
                'Position', [90 bottom 50 inHeight], ...
                'Style', 'edit', ...
                'String', beta, ...
                'Tag', 'betaT');

        uicontrol('Parent', mainWindow, ...
                'Position', [160 bottom 50 inHeight], ...
                'Style', 'text',...
                'String', 'Sigmas=', ...
                'HorizontalAlignment', 'right');

        uicontrol('Parent', mainWindow, ...
                'Position', [200 bottom 50 inHeight], ...
               'Style', 'edit', ...
                'String', sigmas, ...
                'Tag', 'sigmasT');

        uicontrol('Parent', mainWindow, ...
                'Position', [260 bottom 50 inHeight], ...
                'Style', 'text',...
                'String', 'Alpha=', ...
                'HorizontalAlignment', 'right');

        uicontrol('Parent', mainWindow, ...
                'Position', [320 bottom 50 inHeight], ...
               'Style', 'edit', ...
                'String', alpha, ...
                'Tag', 'alphaT');

        uicontrol('Parent', mainWindow, ...
                'Position', [380 bottom 50 inHeight], ...
                'Style', 'text',...
                'String', 'Lambda- =', ...
                'HorizontalAlignment', 'right');

        uicontrol('Parent', mainWindow, ...
                'Position', [440 bottom 50 inHeight], ...
                'Style', 'edit', ...
                'String', lNeg, ...
                'Tag', 'lambdaNT');

        uicontrol('Parent', mainWindow, ...
                'Position', [500 bottom 50 inHeight], ...
                'Style', 'text',...
                'String', 'Lambda+ =', ...
                'HorizontalAlignment', 'right');

        uicontrol('Parent', mainWindow, ...
                'Position', [560 bottom 50 inHeight], ...
                'Style', 'edit', ...
                'String', lPos, ...
                'Tag', 'lambdaPT');

        uicontrol('Parent', mainWindow, ...
                'Position', [620 bottom 30 inHeight], ...
                'Style', 'pushbutton', ...
                'Tag', 'run', ...
		'Interruptible', 'on', ...
                'String', 'run', ...
                'Callback', 'RunningBlob runCB');

	uicontrol('Parent', mainWindow, ...
		'Style', 'pushbutton', ...
		'Position', [660 bottom 30 inHeight], ...
		'String', 'Quit', ...
		'Interruptible', 'off', ...
		'BusyAction', 'cancel', ...
		'Callback', 'RunningBlob quit');


	%%%%%%%%%%%%%Axes for the three graphs%%%%%%%%%%%%%%%%
	graph = axes('Parent', mainWindow, 'Tag', 'graph', ...
			'NextPlot', 'add', 'XGrid', 'on', ...
			'YGrid', 'on', 'Position', ...
			[ .05 .07 .9 .43 ] , 'Visible', 'on');
	inGraph = axes('Position', [ .01 .51 .49 .49], 'NextPlot', ...
			'add', ...
                        'YGrid', 'off', 'XGrid', 'off', ...
                        'Visible', 'on', 'YAxisLocation', 'left', ...
                        'TickDir', 'out', 'Tag', 'inGraph', ...
                        'XAxisLocation', 'top');


	outGraph = axes('Position', [ .51 .51 .49 .49], 'NextPlot', ...
                        'add', 'XAxisLocation', ...
                        'Top', 'YAxisLocation', 'right', ...
                        'TickDir', 'out', 'Tag', 'outGraph', ...
                        'Visible', 'on');

%retrieve values from the textfields.
function getValues
        main = findobj('Tag', 'mainWindow');
        struct = get(main, 'UserData');
        h = findobj('Tag', 'betaT');
        struct.beta = str2num(get(h, 'String'));
        h = findobj('Tag', 'lambdaPT');
        struct.lPos = str2num(get(h, 'String'));
        h = findobj('Tag', 'lambdaNT');
        struct.lNeg = str2num(get(h, 'String'));
        h = findobj('Tag', 'alphaT');
        struct.alpha = str2num(get(h, 'String'));
        h = findobj('Tag', 'sigmasT');
        struct.sigmas= str2num(get(h, 'String'));
        set(main, 'UserData', struct);

function y = geFunction(x)
	y = (exp(-x*x/2)) / (sqrt(2*pi));

function y = step(x)
	if x > 0,
		y = 1;
	else
		y = 0;
	end;

function c = mult_col_mat(a, b)
	sa = size(a);
	sa = sa(1);
	sb = size(b);
	sb = sb(1);
	sio = floor(sa/2);
	simax = sa + sb;
	c = zeros(sb, sb);
	d = zeros(simax, sb);

	d(sio:sio+sb -1, :) = b;
	
	for i = 1: sb,
		for j = 1: sb,
			val = 0;
			for m =1: sa,
				val = a(m) * d(i + m, j);
			end;
			c(i,j) = val;
		end;
	end;
	
function c = mult_vec_mat(a, b)
        sa = size(a);
        sa = sa(1);
        sb = size(b);
        sb = sb(1);
        sjo = floor(sa/2);
        simax = sa + sb;

	c = zeros(sb, sb);
	d = zeros(sb, simax);
	
	d(:,sjo: sjo+sb -1) = b;
		
        for i = 1: sb,
                for j = 1: sb,
                        val = 0;
                        for m =1: sa,
                                val = a(m) * d(i, j + m);
                        end;    
                        c(i,j) = val;  
                end;
        end;

%Initialize matrix, graph function.
function initM
        h = findobj('Tag', 'mainWindow');
        struct = get(h, 'UserData');
        %in the sample code this value never changes.
        d = 0.0;
        ds = 1/(1 + exp(d));

	struct.hIn = struct.hIn .* 0; 
	struct.selfInh = struct.selfInh .* 0;

	struct.deltaX = (2.0 * struct.sigmas)/struct.gSize;

	for k = 1: struct.gSize,
		struct.ge(k) = geFunction((k - struct.gSize/2) ... 
			* struct.deltaX);
	end;
	
	%find the 'graph' axes and set them as current.
	graph = findobj('Tag', 'graph');
	axes(graph);
	plot(struct.ge);
	%set scaling of axis
	axis([ 1 11 0 1 ]);
	drawnow;
	set(h, 'UserData', struct);	

%After calling 'initM' start looping through weight updates and plot them.
function compute
        h = findobj('Tag', 'mainWindow');
        struct = get(h, 'UserData');

	initM;
	
	for t =1 : struct.delta: struct.end_time,
		drawBlobs(struct.hOut, 2);
		for k = 1 : struct.loops,
			struct.hIn(floor(rand * struct.hSize + 1), ...
				floor(rand * struct.hSize + 1))  = ... 
                        struct.hIn(floor(rand * struct.hSize + 1), ...
                                floor(rand * struct.hSize + 1)) + struct.spike;
	
			hOut_sum = 0;
			
			for i = 1: struct.hSize,
				for j = 1: struct.hSize,
					struct.hOut(i,j) = ...
					step(struct.hIn(i,j));
					hOut_sum = hOut_sum + struct.hOut(i,j);
				end; 
			end;
			struct.hTransE = mult_col_mat(struct.ge, ...
				mult_vec_mat(struct.ge, struct.hOut));
			
	
			struct.hTrans = (struct.hTransE - ...
				struct.beta * hOut_sum) * ...
				struct.deltaX * struct.deltaX;
			struct.dyn = struct.hIn - struct.selfInh;

			for i = 1 : struct.hSize,
				for j = 1: struct.hSize,
					if struct.dyn(i,j) > 0,
						struct.dyn(i,j) = ... 
						struct.dyn(i,j)* ...
						struct.lPos;
					else
						struct.dyn(i,j) = ...
						struct.dyn(i,j) * ...
						struct.lNeg;
					end;
				end;
			end;

			struct.selfInh = struct.selfInh + ...
				struct.dyn .* struct.delta;

			%main IDE
			struct.hIn = struct.hIn + (-struct.hIn + ...
				struct.hTrans - struct.alpha .* ...
				struct.selfInh) .* struct.delta;
		end;
		drawBlobs(struct.hIn, 1);
	end;
	drawBlobs(struct.hOut, 2);
	
	set(h,'UserData', struct);

%draw the weight matrixes to either the top left or top right axis.
%if 'side' is 1 then draw to the top left, otherwise to the top right.
function drawBlobs(array_to_plot, side);
	h = findobj('Tag', 'mainWindow');
	struct = get(h, 'UserData');
	if side == 1
		axes_to_place = findobj('Tag', 'inGraph');
	else
		axes_to_place = findobj('Tag', 'outGraph');
	end;
	axes(axes_to_place);
	colormap(gray(20));
	pcolor(array_to_plot);
	axis ij;
	axis square;
	drawnow;

		
