# create an image file of a sparse matrix
# Roman Pearce, CECM/SFU


# M = max dimension of image
# invert = white entries on black background
# gamma = gamma correction to accent non-zero

MatrixImage := proc(A::Matrix, M::posint:=500, {invert::truefalse:=false, gamma::truefalse:=true})
local m,n,X,Y,x,y,s,d,L,img,i,j,TIMER;
	m,n := op(1,A);
	if M >= max(n,m) then
		userinfo(2,'SparseMatrixImage',nprintf("%a x %a matrix -> %a x %a image", m, n, m, n));
		TIMER := time():
		L := LinearAlgebra:-Modular:-IndexList(A);
		userinfo(3,'SparseMatrixImage',nprintf("indices list: %.3f sec", time()-TIMER));
		TIMER := time():
		img := ImageTools:-Create(m,n,1,'background'=`if`(invert, 'black', 'white'));
		d := `if`(invert, 1, 0);
		for i to m do
			for j in L[i] do
				img[i,j] := d;
			end do;
		end do;
		userinfo(3,'SparseMatrixImage',nprintf("construct image: %.3f sec", 1, 1, time()-TIMER));
	else
		s := iquo(max(n,m), M, 'r') + `if`(r=0,0,1);
		X := iquo(m,s,'r') + `if`(r=0,0,1);
		Y := iquo(n,s,'r') + `if`(r=0,0,1);
		d := `if`(invert, 1.0/s^2, -1.0/s^2);
		userinfo(2,'SparseMatrixImage',nprintf("%a x %a matrix -> %a x %a image", m, n, X, Y));
		TIMER := time():
		L := LinearAlgebra:-Modular:-IndexList(A);
		userinfo(3,'SparseMatrixImage',nprintf("indices list: %.3f sec", time()-TIMER));
		TIMER := time():
		img := ImageTools:-Create(X,Y,1,'background'=`if`(invert, 'black', 'white'));
		i := 1;
		for x to X do
			to s while i <= m do
				for j in L[i] do
					y := iquo(j-1,s)+1;
					img[x,y] := img[x,y] + d;
				end do;
				i := i+1;
			end do;
		end do;
		userinfo(3,'SparseMatrixImage',nprintf("%a x %a sampling: %.3f sec", s, s, time()-TIMER));
		if gamma then
			TIMER := time():
			ImageTools:-Gamma(img, `if`(invert, 1/s, s), 'inplace'=true);
			userinfo(3,'SparseMatrixImage',nprintf("gamma correction: %.3f sec", time()-TIMER));
		end if;
	end if;
	img;
end proc:

infolevel[SparseMatrixImage] := 3:




# read MatrixMarket files
# sources of matrices:
#
# University of Florida Sparse Matrix Collection (MATLAB)
# http://www.cise.ufl.edu/research/sparse/matrices/
#
# NIST MatrixMarket
# http://math.nist.gov/MatrixMarket/

MatrixMarket := proc(fname::string, {output::name:=matrix, rhs::string:=""})
local spec, dtype, ex, dims, A, B, f, i, j, x;
	try:
	spec := readline(fname);
	if spec=0 then error "empty file" end if;

	# MatrixMarket specification
	spec := StringTools:-Words(spec);
	if spec[1] <> "MatrixMarket" then error "file not in MatrixMarket format" end if;

	# values read + default [Re,Im]
	if spec[4]="pattern" then
		dtype := "";
		ex := 1, 0;
	elif spec[4]="integer" then
		dtype := "%d";
		ex := 0;
	elif spec[4]="real" then
		dtype := "%f";
		ex := 0;
	elif spec[4]="complex" then
		dtype := "%f %f";
		ex := NULL;
	end if;

	# parse comments
	f := readline(fname);
	try
		while f <> 0 and f[1]="%" do
			userinfo(2, 'MatrixMarket', `NoName`, nprintf(f[2..-1]));
			f := readline(fname);
		end do;
	catch:
		while f <> 0 and f[1]="%" do
			f := readline(fname);
		end do;
	end try:
	dims := sscanf(f, "%d %d %d");

	# construct system
	if output='equations' then
		if spec[3]="coordinate" then
			userinfo(2, 'MatrixMarket', nprintf("equations: %d x %d with %d nonzero", dims[1], dims[2], dims[3]));
			dtype := cat("%d %d ", dtype);	# add coordinate data
			f := readline(fname);
			if member("symmetric", spec) then
				while f <> 0 do
					f := [op(sscanf(f, dtype)), ex];
					A[f[1]][f[2]] := f[3] + f[4]*I;
					A[f[2]][f[1]] := f[3] + f[4]*I;
					f := readline(fname);
				end do;
			else
				while f <> 0 do
					f := [op(sscanf(f, dtype)), ex];
					A[f[1]][f[2]] := f[3] + f[4]*I;
					f := readline(fname);
				end do;
			end if;
			A := [seq(`if`(assigned(A[i]), add(x[j]*A[i][j], j=map(op, [indices(A[i])])), 0), i=1..dims[1])];
			if rhs <> "" then
				B := procname(rhs, ':-output'='equations');
				A := A-B;
			end if;
		elif spec[3]="array" then
			userinfo(2, 'MatrixMarket', nprintf("vector: %d elements", dims[1]));
			f := readline(fname);
			for i from 1 while f <> 0 do
				f := [op(sscanf(f, dtype)), ex];
				A[i] := f[1] + f[2]*I;
				f := readline(fname);
			end do;
			A := [seq(A[j], j=1..i-1)];
			if rhs <> "" then
				B := procname(rhs, ':-output'='equations');
				A := expand(A*x[1]-B);
			end if;
		else
			error "unknown format", spec[3];
		end if;
	elif output='matrix' then
		if spec[3]="coordinate" then
			userinfo(2, 'MatrixMarket', nprintf("matrix: %d x %d with %d nonzero", dims[1], dims[2], dims[3]));
			dtype := cat("%d %d ", dtype);	# add coordinate data
			A := rtable(1..dims[1], 1..dims[2], 'subtype'='Matrix', 'storage'='sparse');
			f := readline(fname);
			if member("symmetric", spec) then
				while f <> 0 do
					f := [op(sscanf(f, dtype)), ex];
					A[f[1], f[2]] := f[3] + f[4]*I;
					A[f[2], f[1]] := f[3] + f[4]*I;
					f := readline(fname);
				end do;
			else
				while f <> 0 do
					f := [op(sscanf(f, dtype)), ex];
					A[f[1], f[2]] := f[3] + f[4]*I;
					f := readline(fname);
				end do;
			end if;
		elif spec[3]="array" then
			userinfo(2, 'MatrixMarket', nprintf("vector: %d elements", dims[1]));
			A := rtable(1..dims[1], 'subtype'='Vector[column]', 'storage'='rectangular');
			f := readline(fname);
			for i from 1 while f <> 0 do
				f := [op(sscanf(f, dtype)), ex];
				A[i] := f[1] + f[2]*I;
				f := readline(fname);
			end do;
		else			
			error "unknown format", spec[3];
		end if;
		if rhs <> "" then
			B := procname(rhs, ':-output'=':-matrix');
			A := A,B;
		end if;
	else
		error "unknown output type", output;
	end if;
	catch:
		try
		fclose(fname);
		catch:
		end try:
		error;
	end try:
	A;
end proc:


# print file comments
infolevel[MatrixMarket] := 2:



# construct matrix from equations
GenerateMatrix := proc(sys::{list,set}, {order::name:='default'})
local ord, cols, vars, eqns, A, c, m, n, i, j, TIMER;
	TIMER := time():
	ord := StringTools:-LowerCase(sprintf(order));

	if ord="default" then
		vars := [op(indets(sys))];
		eqns := [seq([[coeffs(i,indets(i),'m')],[m]], i=sys)];

		# map vars to column index
		cols := proc(v) 0 end proc:  forget(cols);
		for i to nops(vars) do cols(vars[i]) := i end do; cols(1) := i;

	elif ord="sparse" then
		vars := [op(indets(sys))];
		eqns := [seq([[coeffs(i,indets(i),'m')],[m]], i=sys)];
		eqns := map(attributes, sort([seq(setattribute(Float(nops(i[1]),0),i), i=eqns)], `<`));

		# map vars to column index
		cols := proc(v) 0 end proc:  forget(cols);
		for i to nops(vars) do cols(vars[i]) := i end do; cols(1) := i;

	elif ord="markowitz" then
		# sort vars by column weight
		cols := `+`(seq(op(indets(i)), i=sys));
		vars := indets(cols);
		cols := [coeffs(cols,vars,'vars')];
		vars := [vars];
		vars := map(attributes, sort([seq(setattribute(Float(cols[i],0),vars[i]), i=1..nops(cols))], `<`));

		userinfo(3,'GenerateMatrix',nprintf("sort columns: %.3f sec", time()-TIMER));
		TIMER := time():

		# map vars to column index
		cols := proc(v) 0 end proc:  forget(cols);
		for i to nops(vars) do cols(vars[i]) := i end do; cols(1) := i;

		# sort equations by lowest column
		eqns := [seq([[coeffs(i,indets(i),'m')],[m]], i=sys)];
		eqns := map(attributes, sort([seq(setattribute(Float(min(op(map(cols,i[2]))),0), i), i=eqns)], `<`));

	else
		error "unknown order %1", order;
	end if;

	userinfo(3,'GenerateMatrix',nprintf("sort equations: %.3f sec", time()-TIMER));
	TIMER := time():

	# build matrix
	A := Matrix(nops(eqns), nops(vars) + 1, 'storage'='sparse');
	for i to nops(eqns) do
		c, m := op(eqns[i]);
		n := nops(c);
		for j to n do
			A[i,cols(m[j])] := c[j]
		end do;
	end do;
	userinfo(3,'GenerateMatrix',nprintf("%a x %a matrix: %.3f sec", nops(eqns), nops(vars)+1, time()-TIMER));
	A;
end proc:	
	
infolevel[GenerateMatrix] := 3: