1;

function [A,B] = ggllm (X, l, iter=500, alfa=0.01, prec=0.01, link="normal",plotting=0)
#
#	usage:	[A,B] = ggllm (X, l, iter, alfa, prec, link, plot)
#
# The function implements the Generalized^2 Linear^2 Model, a statistical
# estimator by decomposing data matrix X into a simpler representation:
#
# 			X=f(g(A)h(B)),
#
# where A and B are two low-rank matrices;
# f, g, and h are link functions.
#
# -------------------------------------------------------------------------------
#
#	X	- matrix to be analyzed
#	l	- the reduced set of features (the 2nd dimension of A)
#	iter	- number of iterations (default: 500); if 0 max time: 4 mins
#	alpha	- learning rate (default: 0.01)
#	prec	- precision, how close is the prediction to the measured data 
#		  (default: 0.01)
#	link	- link function, depending on the noise on the data matrix X
#	plot	- if you don't want any message to shown, and don't want to
#		  save your avg error results
#
# The selectable link functions are:
# ----------------------------------
# 	- "normal"	= Gaussian noise (default setting)
#	- "poisson"	= Poisson noise 
#	- "binomial"	= Bernoulli noise
#	- "gamma"	= Gamma noise
#	- "inverse"	= Inverse Gaussian noise
#
# -------------------------------------------------------------------------------
#
#Example of usage:
#
#	[a,b] = ggllm (hilb(5), 3, 0);
#
#Results:
#
#converged in 4962 iteration in 8 seconds.
#
#Avarage of errors: 0.002016
#
#error:
#ans =
#
#   0.0195930  -0.0105157  -0.0137456  -0.0133061  -0.0122210
#  -0.0105165   0.0100448   0.0070766   0.0038599   0.0016363
#  -0.0137457   0.0070761   0.0098048   0.0095502   0.0087065
#  -0.0133045   0.0038585   0.0095490   0.0114862   0.0119817
#  -0.0122211   0.0016358   0.0087064   0.0119828   0.0134309
#
#
#Data matrix
#X =
#
#   1.00000   0.50000   0.33333   0.25000   0.20000
#   0.50000   0.33333   0.25000   0.20000   0.16667
#   0.33333   0.25000   0.20000   0.16667   0.14286
#   0.25000   0.20000   0.16667   0.14286   0.12500
#   0.20000   0.16667   0.14286   0.12500   0.11111
#
#Prediction
#ans =
#
#   0.980407   0.510516   0.347079   0.263306   0.212221
#   0.510516   0.323289   0.242923   0.196140   0.165030
#   0.347079   0.242924   0.190195   0.157116   0.134151
#   0.263304   0.196141   0.157118   0.131371   0.113018
#   0.212221   0.165031   0.134151   0.113017   0.097680
#
#a =
#
#   0.13930   0.96704   0.16070
#   0.26716   0.45382   0.21438
#   0.26000   0.28850   0.19838
#   0.23740   0.20861   0.17746
#   0.21486   0.16212   0.15882
#
#b =
#
#   0.13925   0.26718   0.26004   0.23742   0.21489
#   0.96707   0.45380   0.28848   0.20859   0.16210
#   0.16066   0.21440   0.19839   0.17748   0.15884
#
# -------------------------------------------------------------------------------

#-----------------------------
#-Implemented by András Fülöp-
#-----------------------------

MAXTIME=240; #If the number of iterations is set to zero, there's a maximum computation time.
				 #By default it's 4 minutes, you can change it here.  
	
#If there's nothing to decompose
	if (nargin<2)
		printf("\n\tYou have to give more arguments!\n\tFor more help type: help ggllm!\n\n");
		return;
	endif

#Time to start
	start=time();	
	i = 0;

#setting the link function
	if(strcmp(link,"normal")) 
		f=g=h=@normal_function;
	else
		if(strcmp(link,"poisson"))
			f=g=h=@Poisson_function;
		else
			if(strcmp(link,"binomial"))
				f=g=h=@binomial_function;
			else
				if(strcmp(link,"gamma"))
					f=g=h=@gamma_function;
				else
					if(strcmp(link,"inverse"))
							f=g=h=@inverse_function;
					endif
				endif
			endif
		endif
	endif

#generating A and B, and precision matrix
	[n_,m_]=size(X);
	A_old=A_new=A=rande(n_,l).*1e-5;
	B_old=B_new=B=rande(l,m_).*1e-5;
	precmat=ones(n_,m_)*prec;	

#Selecting the type of computing limit
	if(iter==0)
#GGD algorithm			
		while(((sum(sum(abs(X .- f( g(A_old) * h(B_old))).-precmat)))>0)&&((time()-start)<MAXTIME))

			A_new=A_old.+alfa.*(X.-f(g(A_old)*h(B_old)))*h(B_old)';	#computing new value
			A_old=A_new;															#saving data
		
			B_new=B_old.+alfa.*g(A_old)'*(X.-f(g(A_old)*h(B_old)));	#computing new value
			B_old=B_new;															#saving data
			
			if(plotting)
				if(mod(i,50)==0)
					savedata(i,(mean(mean(X .- f( g(A_old) * h(B_old))))));	#save data for plotting
				endif
			endif
			++i;
		endwhile
	else
#GGD algorithm
		while(((sum(sum(abs(X .- f( g(A_old) * h(B_old))).-precmat)))>0)&&(i<iter))
			
			A_new=A_old.+alfa.*(X.-f(g(A_old)*h(B_old)))*h(B_old)';	#computing new value
			A_old=A_new;															#saving data
		
			B_new=B_old.+alfa.*g(A_old)'*(X.-f(g(A_old)*h(B_old)));	#computing new value
			B_old=B_new;															#saving data

			if(plotting)
				if(mod(i,50)==0)
					savedata(i,(mean(mean(X .- f( g(A_old) * h(B_old))))));	#save data for plotting
				endif
			endif
			++i;
		endwhile	
	endif

if(plotting)
	savedata(i,(mean(mean(X .- f( g(A_old) * h(B_old))))));	#save data for plotting
endif

	if(not(plotting))
#printing outputs	
		printf("\nconverged in %d iteration in %d seconds.\n\nAvarage of errors: %2.6f\n", i, time()-start, mean(mean(X .- f( g(A_old) * h(B_old)))));
			A=A_new;
			B=B_new;

		printf("\nerror:\n");
			X.-f(g(A)*h(B))

		printf("\nData matrix\n");
			X

		printf("\nPrediction\n");
			f(g(A)*h(B))
	endif

endfunction	



#save data for plotting
function savedata(i,err)
	fid = fopen("adat","a");
	fprintf(fid,"%d,%2.5f;\n",i,err);
	fclose(fid);
endfunction



#link functions
function y = normal_function (x)
	y=x;
endfunction

function y = Poisson_function (x)
	y=log(x);
endfunction

function y = binomial_function (x)
	y=log(x./1-x);
endfunction

function y = gamma_function (x)
	y=x.^(-1);
endfunction

function y = inverse_function (x)
	y=x.^(-2);
endfunction



function ggllmtest(number=5)
#		
#	usage: ggllmtest(number);
#
#Test function, it computes G^2L^2M, with a 5x5 size Hilbert matrix,
#'number' times; write the avg errors to a data file called 'adat',
#and then read from, to plot the results.
#
#

#-----------------------------
#-Implemented by András Fülöp-
#-----------------------------

	delete("adat");

	disp("Start");	
	exectime=time();

	for i=1:number
		ggllm(hilb(5),3,0,0.01,0.01,"normal",1);
	endfor

	disp(time()-exectime);

	disp("Loading...");
	load -ascii adat;

	disp("Plotting...");
	plot(adat(:,1),adat(:,2));

endfunction