%
%    This file computes Heegaard Floer homology starting form the data in data.m
%    It also needs the files setprod.m, loop.m and maslov.m
%
%    Ciprian Manolescu - May 18, 2005
%
%    bug fixed (about negative genus) - September 15, 2005

[Alpha, Gamma, Beta] = data;

% computes genus
g = length (Alpha);
if (length(Gamma) ~= g) 
	error ('Alpha and Gamma should have the same length.')
end

% computes number of vertices on each beta curve
for i=1:g
	N(i) = length(Beta(i).curve);
end

% changes Beta.point to the absolute values
% makes a new entry Beta.sign with the signs

for i=1:g
	Beta(i).sign = sign (Beta(i).point);
	Beta(i).point = abs (Beta(i).point);
	if any(Beta(i).point == 1)
		error ('The first point on an alpha curve should be on a gamma curve.')
	end
end

% forms a matrix V with the intersections of alpha and beta curves
% first column tells you which alpha curve, second which point on that 
% alpha curve, third which beta curve, fourth which point on the beta
% fifth column is made of 1's - code for an alpha-beta intersection 
% later we'll add beta-gamma and alpha-gamma intersections 

% also writes alpha and beta curves in terms of vertices
% the results are cell arrays Av and Bv
% so far they show 0 when there is an intersection with gamma
% we also make a cell array Aw by eliminating the zeros in Av i.e.
% Aw encodes only the intersections with the beta curves

V = [];
Av = [];
Bv = [];
Aw = [];
Cv = [];

k=1;
for i=1:g
	jj=1;
	for j=2:Alpha(i)
		V(k,1)=i;
		V(k, 2)=j;
		for l=1:g
			for s=1:N(l)
			 if ((Beta(l).curve(s)==i) & (Beta(l).point(s)==j)) 
					V(k, 3)=l;
					V(k, 4)=s;
					V(k, 5)=1;
					Bv{l}(s)=k;
					Av{i}(j)=k; 
					Aw{i}(jj)=k;
					jj = jj+1;
					continue
			  end 
			end
		end 	
		k=k+1;
	end
end


% generates all possible permutations of the beta curves
PBv = perms (Bv);

% finds all generators

G = [];
UV = [];
for i = 1:factorial(g)
	for j = 1:g
        	UV{j} = intersect (Aw{j}, PBv{i,j});
	end
	G=[G; setprod(UV)];
end

% number of generators for the Heegaard Floer complex
[numg, whocares] = size(G);

% now we're adding intersections of alpha and gamma curves to our V
% first column tells you which alpha curve, second which point on that
% alpha curve, third which gamma curve, fourth which point on the gamma.
% fifth column is made of 2's - code for an alpha-gamma intersection

% we also fill in Av, and form a new cell array Cv with the gamma curves

for i=1:g
	V(k, 1)=i;
	V(k, 2)=i;
	V(k, 3)=i;
	V(k, 4)=i;
	V(k, 5)=2;
	Av{i}(1)=k;
	Cv{i}(1)=k;
	k=k+1;
end 

% now we're adding intersections of beta and gamma curves to our V
% first column tells you which beta curve, second which point on that
% beta curve, third which gamma curve, fourth which point on the gamma.
% fifth column is made of 3's - code for an beta-gamma intersection
                                                                                
% we also fill in Bv and Cv

for i=1:g
        for j=2:Gamma(i)
                V(k,3)=i;
                V(k,4)=j;
                for l=1:g
                        for s=1:N(l)
                         if ((Beta(l).curve(s)==-i) & (Beta(l).point(s)==j))
                                        V(k, 1)=l;
                                        V(k, 2)=s;
                                        V(k, 5)=3;
                                        Bv{l}(s)=k;
                                        Cv{i}(j)=k;
                                        continue
                          end
                        end
                end
                k=k+1;
        end
end


% number of 0-cells (vertices) in our complex
[numv, whocares] = size(V);

disp(' ')
disp(['This diagram has ' num2str(numv) ' vertices.'])
answer = input ('Do you want to see them? (y/n) ', 's');
if (answer=='y')
	disp (' ')
	disp ('Each vertex can be identified with a point on exactly two curves.')
	disp ('For example, a2(5) b3(7) means that the vertex is the fifth point on the second') 
	disp ('alpha curve and the seventh point on the third beta curve.')
	disp ('Press any key to continue.')
	disp (' ')
	pause
        for i=1:numv
		if V(i, 5) == 1
			disp ([num2str(i), '  a', num2str(V(i,1)),'(', num2str(V(i,2)), ') b', num2str(V(i,3)), '(', num2str(V(i,4)), ')'])		
		elseif V(i,5) == 2
			disp ([num2str(i), '  a', num2str(V(i,1)),'(', num2str(V(i,2)), ') c', num2str(V(i,3)), '(', num2str(V(i,4)), ')'])
		else
			disp ([num2str(i), '  b', num2str(V(i,1)),'(', num2str(V(i,2)), ') c', num2str(V(i,3)), '(', num2str(V(i,4)), ')'])
		end
	end
	disp (' ')
end

answer='y';
while (answer=='y')
	answer = input ('Do you want to see a curve? (y/n) ', 's');
	if (answer ~= 'y') break, end
	type=input ('Type of curve (a/b/c) : ', 's');
	numc = input (['Number of curve (an integer between 1 and ' num2str(g) '): ']); 
	disp ('Here are the vertices on the curve, in order:')
	if (type=='a')
		disp (Av{numc})
	elseif (type=='b')
		disp (Bv{numc})
	elseif (type=='c')
		disp (Cv{numc})
	end				
end


% reorder the g elements inside each generator, so that the first is a
% point on the first beta curve, the second on the second, etc.
% the result is the matrix Gbeta

Gbeta = [];
                                                                                         
for i=1:numg
        for j =1:g
                whichbeta = V(G(i, j), 3);
                Gbeta(i, whichbeta) = G (i, j);
        end
end
                                                                                         
% the number of generators on each beta curve
                                                                                         
for i = 1:g
        Sizeofbeta(i)=length(Beta(i).curve);
end

% encode all edges into a matrix with three columns
% the first shows the type of curve on which it lies:
% 1 for alpha, 2 for beta, 3 for gamma
% the second shows the number of the respective a/b/c curve
% the third shows the number of its starting point on the curve

% EE is a cell array saying which edges appear on each a/b/c curve

EE = [];
Edge = [];

j=1;
for i = 1:g
	for k =1:(Alpha(i))
		Edge (j, 1) = 1;
		Edge (j, 2) = i;
		Edge (j, 3) = k;
		EE {1, i} (k) = j;
		j = j +1;
	end
end

for i = 1:g
        for k =1:(Sizeofbeta(i))
                Edge (j, 1) = 2;
                Edge (j, 2) = i;
                Edge (j, 3) = k;
		EE {2, i} (k) = j;
                j = j +1;
        end
end

for i = 1:g
        for k =1:(Gamma(i))
                Edge (j, 1) = 3;
                Edge (j, 2) = i;
                Edge (j, 3) = k;
		EE {3, i} (k) = j;
                j = j +1;
        end
end

% number of 1-cells (edges) in our complex
[nume, whocares] = size(Edge);


% find the domains (2-cells) in the diagram
% we encode them in a cell array Dom that indicates
% the edges on the boundary, counterclockwise 
% a minus sign means the orientation on the edge is the opposite one

% we run through all the edges as long as it's necessary
% in the process we form a matrix ED (edges in terms of domains)
% the first column tells the domain on its left
% the second the domain on its right
% the third is 0 at the beginning, 2 at the end (the number of times
% it has appeared as part of a domain)

Dom = [];
ED = zeros (nume, 3);
j = 0;



for i = 1:nume
for u = 1:2
if ((ED (i, u))  == 0)
	ED (i, 3) = ED (i, 3) + 1;
	% construct a new domain to the left or to the right, according to whether u = 1 or 2
	j = j+1;
	if u==1
		edg = i;
		ED (i, 1) = j;
	else
		edg = -i;
		ED (i, 2) = j;
	end
	k=1;
	Dom{j}(1) = edg;
	aedg = abs (edg);
	semn = sign (edg);
	while 1
		if ((Edge(aedg, 1)) == 1) & (semn == 1)
			% what to do for an alpha curve in the right direction !!!!!!!!!!!!!!
			e = Edge (aedg, 2);
			if (Edge(aedg,3) == (Alpha(e)))
				% go up on the gamma curve
				edg = EE {3, e} (1);	
				ED (edg, 1) = j;			
			else
				% go up on the beta curve
				vertex = Av{e}(Edge(aedg,3) +1);
				whichbeta = V (vertex, 3);
				whichpoint = V (vertex, 4);
				if (Beta(whichbeta).sign(whichpoint)) == 1
					edg = EE {2, whichbeta} (whichpoint);	
					ED (edg, 1) = j;
				else				
					if (whichpoint == 1)
						edg = -EE {2, whichbeta} (Sizeofbeta(whichbeta));	
					else
						edg = -EE {2, whichbeta} (whichpoint -1);
					end
					ED (-edg, 2) = j;
				end
			end
		elseif ((Edge(aedg, 1)) == 1) & (semn == -1)
			% what to do for an alpha curve in the wrong direction !!!!!!!!!!!!!!!
                        e = Edge (aedg, 2);
                        if (Edge(aedg,3) == 1)
                                % go down on the gamma curve
                                edg = -EE {3, e} (Gamma(e));
                                ED (-edg, 2) = j;
                        else
                                % go down on the beta curve
                                vertex = Av{e}(Edge(aedg,3));
                                whichbeta = V (vertex, 3);
                                whichpoint = V (vertex, 4);
                                if (Beta(whichbeta).sign(whichpoint)) == (-1)
                                        edg = EE {2, whichbeta} (whichpoint);
                                        ED (edg, 1) = j;
                                else
                                        if (whichpoint == 1)
                                                edg = -EE {2, whichbeta} (Sizeofbeta(whichbeta));
                                        else
                                                edg = -EE {2, whichbeta} (whichpoint -1);
                                        end
                                        ED (-edg, 2) = j;
                                end
			end

		elseif ((Edge(aedg, 1)) == 2) & (semn == 1)
			% what to do for a beta curve in the right direction !!!!!!!!!!!!!!!!!
			e = Edge (aedg, 2);
			vertex = Edge (aedg, 3);
			if vertex == (Sizeofbeta(e))			
				vertex = 1;
			else 
				vertex = vertex + 1;
			end
			
			whichcurve = Beta(e).curve(vertex);
			if (sign (whichcurve)) == 1
				% we continue on an alpha curve
				whichpoint = Beta(e).point(vertex);
				if (Beta(e).sign(vertex)) == 1
					% we go left on the alpha curve
					vertex = Bv{e}(vertex);
					vertex = V(vertex, 2) - 1;
					if (vertex == 0) 
						vertex = Alpha(whichcurve);
					end
					edg = -EE {1, whichcurve} (vertex);
					ED (-edg, 2) = j;
				else
					% we go right on the alpha curve
					vertex = Bv{e}(vertex);
					vertex = V (vertex, 2);
					edg = EE {1, whichcurve} (vertex);
					ED (edg, 1) = j;
				end

			else
				% we continue on a gamma curve
				whichcurve = - whichcurve;
				whichpoint = Beta(e).point(vertex);
                                if (Beta(e).sign(vertex)) == -1
                                        % we go clockwise on the gamma curve
                                        vertex = Bv{e}(vertex);
                                        vertex = V(vertex, 4) - 1;
                                        if (vertex == 0)
                                                vertex = Gamma(whichcurve);
                                        end
                                        edg = -EE {3, whichcurve} (vertex);
                                        ED (-edg, 2) = j;
                                else
                                        % we go counterclockwise on the gamma curve
                                        vertex = Bv{e}(vertex);
                                        vertex = V (vertex, 4);
                                        edg = EE {3, whichcurve} (vertex);
                                        ED (edg, 1) = j;
                                end				
			end


		elseif ((Edge(aedg, 1)) == 2) & (semn == -1)
			% what to do for a beta curve in the wrong direction !!!!!!!!!!!!!!!!!!!
			e = Edge (aedg, 2);
                        vertex = Edge (aedg, 3);

                        whichcurve = Beta(e).curve(vertex);
			if (sign (whichcurve)) == 1
                                % we continue on an alpha curve
                                whichpoint = Beta(e).point(vertex);
                                if (Beta(e).sign(vertex)) == -1
                                        % we go left on the alpha curve
                                        vertex = Bv{e}(vertex);
                                        vertex = V(vertex, 2) - 1;
                                        if (vertex == 0)
                                                vertex = Alpha(whichcurve);
                                        end
                                        edg = -EE {1, whichcurve} (vertex);
                                        ED (-edg, 2) = j;
                                else
                                        % we go right on the alpha curve
                                        vertex = Bv{e}(vertex);
                                        vertex = V (vertex, 2);
                                        edg = EE {1, whichcurve} (vertex);
                                        ED (edg, 1) = j;
                                end

                        else
                                % we continue on a gamma curve
                                whichcurve = - whichcurve;
                                whichpoint = Beta(e).point(vertex);
                                if (Beta(e).sign(vertex)) == 1
                                        % we go clockwise on the gamma curve
                                        vertex = Bv{e}(vertex);
                                        vertex = V(vertex, 4) - 1;
                                        if (vertex == 0)
                                                vertex = Gamma(whichcurve);
                                        end
                                        edg = -EE {3, whichcurve} (vertex);
                                        ED (-edg, 2) = j;
                                else
                                        % we go counterclockwise on the gamma curve
                                        vertex = Bv{e}(vertex);
                                        vertex = V (vertex, 4);
                                        edg = EE {3, whichcurve} (vertex);
                                        ED (edg, 1) = j;
                                end
                        end


		elseif ((Edge(aedg, 1)) == 3) & (semn == 1)
			 % what to do for a gamma curve in the right direction !!!!!!!!!!!!!!!!!
			e = Edge (aedg, 2);
			if (Edge(aedg,3) == (Gamma(e)))
                                % go left on the alpha curve
                                edg = -EE {1, e} (Alpha(e));
                                ED (-edg, 2) = j;
                        else
                                % go inward on the beta curve
                                vertex = Cv{e}(Edge(aedg,3) +1);
                                whichbeta = V (vertex, 1);
                                whichpoint = V (vertex, 2);
                                if (Beta(whichbeta).sign(whichpoint)) == (-1)
                                        edg = EE {2, whichbeta} (whichpoint);
                                        ED (edg, 1) = j;
                                else
                                        if (whichpoint == 1)
                                                edg = -EE {2, whichbeta} (Sizeofbeta(whichbeta));
                                        else
                                                edg = -EE {2, whichbeta} (whichpoint -1);
                                        end
                                        ED (-edg, 2) = j;
                                end
                        end



		elseif ((Edge(aedg, 1)) == 3) & (semn == -1)
			 % what to do for a gamma curve in the wrong direction !!!!!!!!!!!!!!!!!!!
			e = Edge (aedg, 2);
                        if (Edge(aedg,3) == 1)
                                % go right on the alpha curve
                                edg = EE{1, e}(1);
                                ED (edg, 1) = j;
                        else
                                % go down on the beta curve
                                vertex = Cv{e}(Edge(aedg,3));
                                whichbeta = V (vertex, 1);
                                whichpoint = V (vertex, 2);
                                if (Beta(whichbeta).sign(whichpoint)) == 1
                                        edg = EE {2, whichbeta} (whichpoint);
                                        ED (edg, 1) = j;
                                else
                                        if (whichpoint == 1)
                                                edg = -EE {2, whichbeta} (Sizeofbeta(whichbeta));
                                        else
                                                edg = -EE {2, whichbeta} (whichpoint -1);
                                        end
                                        ED (-edg, 2) = j;
                                end
                        end


		end

		aedg = abs (edg);
		semn = sign (edg);
		if (edg == (Dom{j}(1))), break, end
		ED (aedg, 3) = ED (aedg, 3) + 1;
		k = k+1;
		Dom{j}(k) = edg;
	end		
end
end
end

% number of domains
numd = j;

disp (' ')
disp(['This diagram has ' num2str(numd) ' two-cells.'])


% display the domains in terms of vertices
% encode them into a cell array Dov similar to Dom
% Typedom tells the number of vertices for each domian

Dov = [];
Typedom = [];
for i=1:numd
	Typedom(i) = length(Dom{i});
	for j=1:Typedom(i)
		edg = Dom{i}(j);
		aedg = abs(edg);
		semn = sign(edg);		
		if semn == 1
			if (Edge(aedg,1) == 1)
				% on an alpha curve			
				Dov{i}(j) = Av{Edge(aedg,2)}(Edge(aedg,3));
			elseif (Edge(aedg,1) == 2)			
				% on a beta curve
				Dov{i}(j) = Bv{Edge(aedg,2)}(Edge(aedg,3));
			else
				% on a gamma curve
				Dov{i}(j) = Cv{Edge(aedg,2)}(Edge(aedg,3));
			end
		else
			if (Edge(aedg,1) == 1)
                                % on an alpha curve
				if (Edge(aedg,3) == Alpha(Edge(aedg,2))) 
					vertex = 1;
				else
					vertex = Edge(aedg,3) + 1;
				end				
                                Dov{i}(j) = Av{Edge(aedg,2)}(vertex);
                        elseif (Edge(aedg,1) == 2)
                                % on a beta curve
				if (Edge(aedg,3) == Sizeofbeta(Edge(aedg,2)))
                                        vertex = 1;
                                else
                                        vertex = Edge(aedg,3) + 1;
                                end
                                Dov{i}(j) = Bv{Edge(aedg,2)}(vertex);
                        else
                                % on a gamma curve
				if (Edge(aedg,3) == Gamma(Edge(aedg,2)))
                                        vertex = 1;
                                else
                                        vertex = Edge(aedg,3) + 1;
				end
                                Dov{i}(j) = Cv{Edge(aedg,2)}(vertex);
			end
		end
	end	
end



answer = input ('Do you want to see them? (y/n) ', 's');
if (answer=='y')
	disp (' ')
	disp ('Each two-cell is shown in terms of its vertices:')
	disp (' ')
        for i=1:numd
		disp (['Cell ' num2str(i) ' : '])
		disp (Dov{i})
	end
end



% set the basepoint
basept = input ('In which cell do you want to put the basepoint? ');

% the matrices BA and BG contain the intersection numbers
% between beta(i) and alpha(j), and beta(i) and gamma(j), respectively

BA = zeros (g, g);
BG = zeros (g, g);
for i=1:g
        for j=1:Sizeofbeta(i)
                if (Beta(i).curve(j)) > 0
                         ind = Beta(i).curve(j);
                         BA(i, ind) = BA(i, ind) + Beta(i).sign(j);
                else
                        ind = -Beta(i).curve(j);
                        BG(i, ind) = BG(i, ind) + Beta(i).sign(j);
                end
        end
end




if (det(BA) == 0)
        error ('This is not a rational homology sphere.')
end


% we decompose the generators in G into spin^c structures
% the vector Spstr has Spstr(i) = the spinc structure of G(i)

Spstr = zeros (numg, 1);
index = 1;
spc = 1;
jj=1;

while 1
        Spstr (index) = jj;
        for ii = index : numg
                if ((Spstr(ii)) == 0)
                        pone = index;
                        ptwo = ii;
                        loop;
                        if (spinc == 1)
                                Spstr (ii) = jj;
                        end
                end
        end
        ii = index;
        while ii < (numg+1)
                if ((Spstr(ii)) == 0), break, end
                ii=ii+1;
        end
        if (ii == (numg + 1)) break, end
        index = ii;
        jj = jj+1;
end

% number of spin^c structures
nums = jj;

% count the elements in each spin^c structure

Spincount = zeros (nums, 1);
for i = 1:numg
        j = Spstr (i);
        Spincount (j) = Spincount (j) + 1;
end
Spincount = Spincount' ;


disp (' ')
disp (['This diagram has ' num2str(nums) ' spinc structures.'])
disp ('Number of elements in each spinc structure: ')
disp (Spincount)


% make a matrix NewG by adding a first column with numbers 1:numg, and one with
% the corresponding spin^c structures

New = [1:numg]';
NewG = [New, G, Spstr];

disp(['This diagram has ' num2str(numg) ' generators.'])
answer = input ('Do you want to see them? (y/n) ', 's');
if (answer=='y')
	disp (' ')
	disp (['Each generator consists of ', num2str(g), ' vertices. After an index, we show']) 
	disp ('these vertices and the corresponding spinc structure.') 
	disp ('Press any key to continue.')
	pause
	disp (' ')
        disp (NewG)
end


while 1

whichspin = input ('Which spinc structure are you interested in? ');
disp (' ')
disp ('Here are the generators and their relative Maslov gradings: ')
disp (' ')

% the vector Ourspin contains all the elemnts in that spinc structure
Ourspin = [];
for i = 1:numg
        if (Spstr(i) == whichspin)
                Ourspin = [Ourspin; i];
        end
end

% compute the Maslov indices with respect to Ourspin(1)
% encode the result in Ourmaslov

Ourmaslov = [];

ss = Spincount(whichspin);
ptwo = Ourspin (1);
for ii = 1:ss
        pone = Ourspin(ii);
	power = 0;
        maslov;
        Ourmaslov (ii) = round(masl);
end

minm = min(Ourmaslov);
maxm = max(Ourmaslov);
diff = maxm - minm + 1;

% put the elements in order of their gradings
% encode the result into a cell array Result

for i=1:diff
        Result{i} = [];
end
for i = 1:ss
        j = Ourmaslov (i) - minm + 1;
        Result{j} = [Result{j}, Ourspin(i)];
end

% display the result
for i=1:diff
        disp (['Grading ', num2str(diff -i), ':  ', num2str(Result{diff + 1 - i}) ])
end


% guess the differentials

disp (' ')
disp ('For each domain of index 1, a potential holomorphic disk in the symmetric')
disp ('product has a preimage in Sym^{g-1}(\Sigma) \times Sigma which consists of t')
disp ('trivial disks and a connected surface of genus g with b boundary components.')
disp (' ')
disp ('There are positive domains of index 1 between the following pairs: ')
disp (' ')

% look at domains with power of U = 0 first
for ii=1:(diff-1)
	for jj = 1:length(Result{diff + 1 - ii})
		pone = Result{diff + 1 - ii}(jj);
		for kk = 1:length(Result{diff - ii})
			ptwo = Result{diff - ii}(kk);
			power = 0;		
			maslov;
			if all(MD > (-1))
				% added in program: Sept. 14, 2005
				% the quantity genus can be negative when the domain is disconnected, and more than one piece is nontrivial
				% in that case we have a positive domain of index 1, but it cannot support any differentials, so we do not display it
				if (genus > (-1))
					disp ([ '(', num2str(pone), ',' ,num2str(ptwo), ') t=', num2str(triv), ' g=', num2str(genus), ' b=', num2str(cycles)])
				end
			end 			
		end
	end
end

% look at domains with positive powers of U
for hh = 1:(ceil((diff-1)/2))
	for ii=(diff+1-(2*hh)):-1:1
        for jj = 1:length(Result{ii})
                pone = Result{ii}(jj);
                for kk = 1:length(Result{ii - 1 + (2*hh)})
                        ptwo = Result{ii-1+(2*hh)}(kk);
                        power = hh;
                        maslov;
			if (masl ~= 1)
				error ('Something is terribly wrong.');
			end
			if all(MD > (-1)) 
			% added in program: Sept. 14, 2005
                        % the quantity genus can be negative when the domain is disconnected, and more than one piece is nontrivial
                        % in that case we have a positive domain of index 1, but it cannot support any differentials, so we do not display it
                                if (genus > (-1))
	                                disp ([ '(', num2str(pone), ',(U^' ,num2str(power), ')', num2str(ptwo), ') t=', num2str(triv), ' g=', num2str(genus), ' b=', num2str(cycles)]);
				end
			end
                end
        end
	end
end





% look more closely at differentials
answer='y';
while (answer=='y')
	disp (' ')
        answer = input ('Do you want to look at a domain? (y/n) ', 's');
        if (answer ~= 'y') break, end
        pone=input ('First generator: ');
	power = input ('Power of U for the second generator: ');
	ptwo = input ('Second generator: ');
	maslov;
	disp (' ')
	disp ('This domain consist of the following two-cells, with the given multiplicities:')
	disp (' ')
	disp(MDred)
end


disp (' ');
answer = input ('Do you want to go to a new spinc structure? (y/n) ', 's');
if (answer ~= 'y') break, end
end