%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Source coding with Huffman codes" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Alexander Knaub
% date: 04.03.2016   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This function creates a huffman code for a given source with specified 
% ocurrence probabilities. The result is a tree graph annotated at leaves and branches.


function [tikz_link, pdf_link] = create_code(...
	up_bit,...
	extend_en,... 
    annotate_joins,...
	src_type,...  
	n_symbols,... 
	prob_file,... 
	filename)


	% define formats
	landscape = '-S930,350';
	portrait  = '-S640,480';

	% select format
	output_format = portrait;

	% Create an invisible figure.
	fh = figure(1);
	set(fh, 'visible', 'off');

	% define style
	axes_font_size = 12;
	text_font_size = 14;
	set(0, 'defaultlinelinewidth', 4);
	set(0, 'defaultaxesfontsize', axes_font_size);
	set(0, 'defaultaxesfontname', 'Helvetica');
	set(0, 'defaulttextfontsize', text_font_size);
	set(0, 'defaulttextfontname', 'Helvetica');



	p1 = [.1 .3 .3 .2 .1];
	p2 = [.1 .7 .2];

	[p11 p12] = meshgrid(p1, p2);
	p = p11 .* p12;
	p = [0.11 0.08 0.08 0.08 0.07 0.07 0.07 0.07,... 
		 0.06 0.05 0.05 0.05 0.04 0.04 0.04 0.04];
	
	
	if(strcmp(src_type, 'hard_coded1'))
		fprintf('Hard coded source\n');
		probs = [0.11 0.08 0.08 0.08 0.07 0.07 0.07 0.07,... 
			0.06 0.05 0.05 0.05 0.04 0.04 0.04 0.04];
		probs=reshape(probs, [], 1);
		
	elseif(strcmp(src_type, 'random'))
		if(n_symbols > 36 || n_symbols < 2)
			disp('Bad number of symbols!');
			error_text = ['\n\nBad number of symbols entered!\n',... 
			'Please provide an integer value between 2 and 36 for the number of symbols!\n'];
			error(error_text) ;
		end;
		probs = abs(rand(n_symbols,1));
		probs = probs ./ sum(probs);
		% quantize to 2 decimal digits:
		probs = round(probs*100)/100;
		probs(probs==0) = 0.01;
		quant_error = 1 - sum(probs);
		% re-normalize by adding the total quanization error 
		% to the first max probable symbol:
		index = find(probs==max(probs));
		probs(index(1)) += quant_error;
		
	elseif(strcmp(src_type, 'uploaded'))
		if(strcmp(prob_file, 'undefined'))
			disp('Undefined Source!');
			error_text = ['\n\nPlease specify a source text file!\n',... 
			'The file has to contain at least two symbols, one symbol per line.\n',...
			'Each line has to be formatted like <symbol_name> <symbol_probability>\n\n'];
			error(error_text) ;
		end;
		fprintf('Uploaded source with file %s\n', prob_file);
		[names probs] = parse_probs(prob_file);
	end;
	
    
   	if abs(sum(probs)-1) > 1e-3 
		error(sprintf('Bad distribution! Sum is not 100%% (found %0.2f%%).', sum(probs)*100));
	end;
	
	if(up_bit ~= 1 && up_bit ~= 0)
		error('The selected bit labelling is not yet implemented!');
	end;
	
    names = [1:length(probs)];
	
	% Perform the code creation algorithm:
	tic;
    [codewords, joined_symbols] = make_huffcode(probs, up_bit);
    fprintf('Tree created in %0.2fs\n', toc);
	
	% Display the results:
	tic;
    cw_length = zeros(length(codewords), 1);
    k = 1;
	p_zeros = 0;
	p_ones = 0;
	
    for cw = codewords
        cw_length(k) = length(cw{1});
        p_ones += sum(cw{1}) * probs(k);
        p_zeros += sum(abs(1-cw{1})) * probs(k);
        k = k+1;
    end;
	mean_cw_len = sum(cw_length.*probs);
	p_ones=p_ones/mean_cw_len;
	p_zeros=p_zeros/mean_cw_len;
	entropy = sum(probs .* log2(1./probs));
	
	fprintf('Avg. %0.2f ones and %0.2f zeros per codeword\n', p_ones, p_zeros);
	
	
	% show the tree:
	symbols_row = '\\node{%d}; & \\node{%0.2f}; & \\node{%s};& \\node[dspnodefull, tree line](sym%d){};\\\\ ';
	extend_parallel = ['\draw[tree line] (sym%left) -- (sym%left-|sym%prev) coordinate[at end](sym%left);',...
                        char(10),...
                        '\draw[tree line] (sym%right) -- (sym%right-|sym%prev) coordinate[at end](sym%right);'];
    join_start = '\draw[tree line] (sym%left) -- (sym%left-|{$(sym%right)+(\tree_step,0)$}) coordinate[at end](sym%left);';
    collision_resolution = '\draw[tree line] (sym%d) -- ($(sym%d)+(\tree_step,0)$) coordinate[at end](sym%d);';
    symbols_connection = ['\draw[tree line] (sym%left) coordinate (last_node) -- (last_node-|{$(sym%right)+(\tree_step,0)$}) node[bit assignment]{%leftb} ',... 
                            '-- ($(sym%right)+(\tree_step,0)$) node[midway,dspnodefull](sym%new){} node[bit assignment]{%rightb} -- (sym%right);'];
                            
    new_symbol_draw = '\node[probability] at (sym%new) {%prob}; \node[new symbol name] at (sym%new){\small(%new)};';
	symbols_matrix = '';
	join_tree = '';
	tikz_template = fopen('tikz_template.tex');
    nlines = fskipl (tikz_template, Inf)
    frewind(tikz_template)
    tikz_code = ''
    for line=0:nlines
        cur_line = fgets(tikz_template);
        tikz_code = sprintf('%s%s', tikz_code, cur_line);
    end;
    fclose(tikz_template);
    
	if(extend_en~=0)
		% extended plot requested... 
		% Tree definitely gets maximally wide wide 
		% => leave enough space for the table...
		space_for_tree = length(probs) % in plot units
		space_for_table = 3 + max(cw_length) + 2 % characters
		space_for_table *= (text_font_size)*.8 % in pixels,  0.8 for font aspect ratio
		% 640 - 68 = 592 Pixels for the plot area going from 0 to x + space_for_tree in plot inits
		% (x+space_for_tree) is 592px = space_for_table + y*space_for_tree
		x_unit_px = (592 - space_for_table)/space_for_tree;
		x_offset = space_for_table / x_unit_px + 1;
		cw_separator_x = 3.1*text_font_size/x_unit_px
		probs_position_x = 10 / x_unit_px
	else
		% make the tree maximally compact
		x_offset = 4;
		% separator line between probability values and codewords
		cw_separator_x = 1.45;
		probs_position_x = 0.5;
    end;
	
	
	for (join = 0:floor(joined_symbols(end,1)/2))
        if(join == 0)
			% prepare the axes
            axis 'ij';

            set(gca,'YTick', 1:length(probs));
            set(gca,'YTickLabel',names(joined_symbols(1:length(probs),1)));
            
			if(extend_en ~= 0)
				xtick_pos = round((x_offset-1+[0:length(probs)])*100)/100;
				fprintf('%s\n', mat2str(xtick_pos));
				set(gca,'XTick', xtick_pos);
				xtick_labels = [1:length(probs)+1];
				fprintf('%s\n', mat2str(xtick_labels));
				set(gca,'XTickLabel',xtick_labels);
	            set(gca,'XLabel', 'Iteration #');
			else
				set(gca,'XTick', []);
			end;


            set(gca,'YLabel', 'Symbol #');
            box off; grid off;
            hold on;
			title(sprintf(...
				'N=%d; H(X)=%0.2fbit/Symb; avg. CW-length m=%0.2fbit/CW; H_0(X)=%0.2fbit/Symb; P[0]=%0.02f; P[1]=%0.02f',... 
				length(probs), entropy, mean_cw_len, log2(length(probs)), p_zeros, p_ones));
				
            tikz_code = strrep(tikz_code, '*n_sym*', sprintf('%d', length(probs)));
            tikz_code = strrep(tikz_code, '*entropy*', sprintf('%0.2f', entropy));
            tikz_code = strrep(tikz_code, '*cw_len*', sprintf('%0.2f', mean_cw_len));
            tikz_code = strrep(tikz_code, '*max_entropy*', sprintf('%0.2f', log2(length(probs))));
            tikz_code = strrep(tikz_code, '*p0*', sprintf('%0.2f', p_zeros));
            tikz_code = strrep(tikz_code, '*p1*', sprintf('%0.2f', p_ones));
			symbol_position = zeros(length(joined_symbols(:,1)),2);
			for k=1:length(probs)
                text(probs_position_x, k, sprintf('%0.2f', joined_symbols(k,4)));

                symbol_position(joined_symbols(k,1),:) = [x_offset-1 k];
                symbols_matrix = sprintf('%s\n%s', symbols_matrix,... 
                        sprintf(symbols_row, ...
                                joined_symbols(k,1), ...
                                joined_symbols(k,4),... 
                                num2str(codewords{joined_symbols(k,1)}'), ...
                                joined_symbols(k,1)));
            end;
			
			text(0.1, 0, 'P[x_i]');
			text(cw_separator_x*1.1, 0, 'CW_i');
			% separator line between probability values and codewords
			plot([cw_separator_x cw_separator_x], [.5 length(probs)+.5], '-k');
			tikz_code = strrep(tikz_code, '%symbols_matrix%', symbols_matrix);
%			fprintf(tikz_code)
			
		else
            new_symbol = joined_symbols(length(probs)+join,1);
            more_prob = joined_symbols(length(probs)+join,2);
            less_prob = joined_symbols(length(probs)+join,3);
			
			if symbol_position(more_prob,1) < symbol_position(less_prob,1)
                left_sym = more_prob;
                right_sym = less_prob;
                left_bit = up_bit;
                right_bit = 1-up_bit;
            else
                left_sym = less_prob;
                right_sym = more_prob;
                left_bit = 1-up_bit;
                right_bit = up_bit;
            end;


			if(extend_en ~= 0)
				new_position =... 
					[x_offset,... 
						0.5*(symbol_position(more_prob,2)... 
						+ symbol_position(less_prob,2))];
				x_offset +=1
				
                tikz_current = strrep(extend_parallel, '%left', sprintf('%d', left_sym));
                tikz_current = strrep(tikz_current, '%right', sprintf('%d', right_sym));
                tikz_current = strrep(tikz_current, '%prev', sprintf('%d', new_symbol-1));
                join_tree = sprintf('%s\n%s', join_tree, sprintf('%s', tikz_current));
				
			else
				new_position =... 
					[1+max(symbol_position(more_prob,1),... 
						   symbol_position(less_prob,1)),... 
						0.5*(symbol_position(more_prob,2)... 
						+ symbol_position(less_prob,2))];
						
                
                tikz_current = strrep(join_start, '%left', sprintf('%d', left_sym));
                tikz_current = strrep(tikz_current, '%right', sprintf('%d', right_sym));
                join_tree = sprintf('%s\n%s', join_tree, tikz_current);
				% resolve possible collisions:
				collision_resolved = 0;
				while(collision_resolved == 0)
					if(sum(symbol_position(1:new_symbol-1,1)==new_position(1)) == 0)
						% This symbol is alone!
						collision_resolved = 1;
						continue;
					end;
					% If we got to this point, a joined symbol has been at this position
					% is there a collision? => check vertical coordinates
%					fprintf('Possible collision of symbol %d(%0.02f) @ (%0.1f; %0.2f)\n',... 
%							 	new_symbol, joined_symbols(new_symbol, 4),... 
%								new_position(1), new_position(2));
					y1 = symbol_position(more_prob, 2);
					y2 = symbol_position(less_prob, 2);
%					fprintf('\tJoining coords %0.2f and %0.2f\n', y1, y2);
					collision_symbols = joined_symbols(symbol_position(1:new_symbol-1,1)==new_position(1),:);
					collision_detected = 0;
					for k = 1:length(collision_symbols(:,1))
						more_prob_old = symbol_position(collision_symbols(k,2),2);
						less_prob_old = symbol_position(collision_symbols(k,3),2);
%						fprintf('Checking symbol %d (%0.2f) %0.2f <-> %0.2f\n',... 
%									collision_symbols(k,1), collision_symbols(k,4),... 
%									more_prob_old, less_prob_old);
						if(sum(max(y1, y2) > min(more_prob_old, less_prob_old)) ~=0 )
%							fprintf('...Collision detected!\n');
							new_position(1)+=1;
							tikz_current = strrep(collision_resolution, '%d', sprintf('%d', more_prob));
							join_tree = sprintf('%s\n%s', join_tree, tikz_current);
							tikz_current = strrep(collision_resolution, '%d', sprintf('%d', less_prob));
							join_tree = sprintf('%s\n%s', join_tree, tikz_current);
							collision_detected = 1;
							break;
						else
%							fprintf('... No collision!\n');
						end;					
					end;
					collision_resolved = 1 - collision_detected;
				end;
			end;

			
			symbol_position(new_symbol,:) = new_position;
			

			
			% connect the symbols:
            tikz_current = strrep(symbols_connection, '%leftb', sprintf('%d', left_bit));
            tikz_current = strrep(tikz_current, '%left', sprintf('%d', left_sym));
            tikz_current = strrep(tikz_current, '%rightb', sprintf('%d', right_bit));
            tikz_current = strrep(tikz_current, '%right', sprintf('%d', right_sym));
            tikz_current = strrep(tikz_current, '%new', sprintf('%d', new_symbol));
			join_tree = sprintf('%s\n%s', join_tree, tikz_current);
			
			tikz_current = strrep(new_symbol_draw, '%new', sprintf('%d', new_symbol));
			tikz_current = strrep(tikz_current, '%prob', sprintf('%0.2f', joined_symbols(new_symbol, 4)));
%  			sprintf('%s\n%s', tikz_current, sprintf('\\node[probability] at (sym%d) {%0.2f}; ', new_symbol, joined_symbols(new_symbol, 4)));
            join_tree = sprintf('%s\n%s', join_tree, tikz_current);
            
            % horizontal lines up to this point
            plot([symbol_position(more_prob, 1) symbol_position(new_symbol,1)],...
                 [symbol_position(more_prob, 2) symbol_position(more_prob, 2)], 'k');
            plot([symbol_position(less_prob, 1) symbol_position(new_symbol,1)],...
                 [symbol_position(less_prob, 2) symbol_position(less_prob,2)], 'k');

            % vertical line:
		    plot([symbol_position(new_symbol,1) symbol_position(new_symbol,1)],...
	              [symbol_position(more_prob,2) symbol_position(less_prob,2)] , 'k');
		plot(new_position(1), new_position(2), 'k', 'markersize', 12);
			text(new_position(1)+0.1, new_position(2)-0.2,... 
				 sprintf('%0.2f', joined_symbols(length(probs)+join,4)),... 
				 'color', 'red');
		end;
	end; % end displaying the tree
%	text(0.5,0.5, 'Test!');
	new_symbol = joined_symbols(end,1);
	
	tikz_current = '\draw[tree line] (sym%d) -- ($(sym%d)+(\tree_step,0)$);';
	tikz_current = strrep(tikz_current, '%d', sprintf('%d', new_symbol));
	join_tree = sprintf('%s\n%s', join_tree, tikz_current);
	
    plot([symbol_position(new_symbol, 1) symbol_position(new_symbol,1)+1],...
         [symbol_position(new_symbol, 2) symbol_position(new_symbol, 2)], 'k');
	axis([0, symbol_position(new_symbol, 1)+1, -.5, length(probs)+1]);	

	
	% fill in the codewords:	
    for(k=1:length(codewords))
        text(cw_separator_x+.1, k, sprintf('%c', num2str(codewords{joined_symbols(k,1)}')));
    end;
	fprintf('Plot created in %0.2fs\n', toc);
	tic;
	% Create an image file. This png image is displayed on the webpage. Its size 
	% must be 640x480.
	print(filename, '-dpng', output_format);
	fprintf('Image saved in %0.2fs\n', toc);
	
	tikz_code = strrep(tikz_code, '%join_tree%', join_tree);
	
	basename = strrep(filename, '.png', '');
    tikz_file = fopen(sprintf('%s.tex', basename), 'w');
    fwrite(tikz_file, tikz_code);
    fclose(tikz_file)
    
    tikz_link = ['<a href=./',sprintf('%s.pdf', basename),'> Hi-Res PDF image</a>'];
    pdf_link = ['<a href=./',sprintf('%s.tex', basename),'> TikZ code of the tree</a>'];
    
    
%      cmd = sprintf('pdflatex -output-directory=./tmp %s.tex; pdf2ps -eps %s.pdf; convert -density 200 %s.eps %s.tikz.png', filename, filename, filename, filename);
    cmd = 'pdflatex -output-directory=./tmp %s.tex; pdf2ps %s.pdf %s.ps; ps2eps -q %s.ps; convert -density 200 -resize 640x480 %s.eps %s.tikz.png ';
    cmd = strrep(cmd, '%s', basename); 
%                  filename,... % TEX
%                  filename,... % PDF
%                  filename,... % PS from PDF 
%                  filename,... % PS 
%                  filename,... % EPS 
%                  filename... % PNG 
%                  );
    fprintf('about to execute the command\n%s', cmd);
    [status, output] = system(cmd);
    copyfile(sprintf('%s.tikz.png', basename), filename);
    
    
end;

function [codewords, joined_symbols] = make_huffcode(probs, top_bit)
    
   
    joined_symbols = zeros(2*length(probs)-1, 4);
    n_joins = 0;
    bottom_bit = abs(top_bit - 1);
    
    
    symbols_probs = [[1:length(probs)]' probs];
    symbol_position = zeros(length(probs)*2-1,2);
    
    
    
    %% build the tree:
    while (size(symbols_probs, 1) > 1)
        
        % save backup
        prev_symprobs = symbols_probs;
        
        % re-sort the symbols:
        symbols_probs = huffman_sort(symbols_probs, n_joins);
        
        % save two symbols jont at this iteration
        if(n_joins == 0)
            joined_symbols(1:length(probs),1) = symbols_probs(1:end,1);
            joined_symbols(1:length(probs),2:3) = NaN*ones(length(probs),2); 
            joined_symbols(1:length(probs),4) = symbols_probs(1:end,2);
                
            joined_symbols(length(probs)+1,1:3) =... 
                [length(probs)+1, ...
                 symbols_probs(end-1,1),...  
                 symbols_probs(end,1)];
        else
            joined_symbols(length(probs) + n_joins+1,1:3) =... 
                [joined_symbols(length(probs) + n_joins,1)+1,... 
                 symbols_probs(end-1,1),...  
                 symbols_probs(end,1)];
        end;
        
        % replace the two least-probable symbols by the new one
        if(length(symbols_probs) ~= 2)
            symbols_probs = [symbols_probs(1:end-2,:);... 
                             [joined_symbols(length(probs) + n_joins+1,1),... 
                              symbols_probs(end-1,2) + symbols_probs(end,2)]];
        else
            symbols_probs = [joined_symbols(length(probs) + n_joins+1,1),... 
                              symbols_probs(1,2) + symbols_probs(2,2)];
        end;
        joined_symbols(length(probs) + n_joins+1, 4) = symbols_probs(end,2);
        n_joins += 1;
    end;
    
    % create the codewords:
    codewords = cell(joined_symbols(length(probs) + n_joins,1), 1)';
    
    root_index = joined_symbols(end,1);
    codewords = fill_up(joined_symbols, codewords, root_index, top_bit, bottom_bit);
    codewords = fill_down(joined_symbols, codewords, root_index, top_bit, bottom_bit);
    
    codewords = codewords(1:(joined_symbols(end,1)+1)/2);
    
   
    
end;

function sorted = huffman_sort(symbols_probabilities, presorted)
    
    if(presorted == 0)
        % use octave's sort algorithm:
        [s i]  = sort(symbols_probabilities, 1, 'descend');
        sorted = [symbols_probabilities(i(:,2),1), s(:,2)];
    else
        % only the last element has changed, find the topmost position for it:
        to_insert = symbols_probabilities(end,:);
        k = 1;

        while(symbols_probabilities(k,2) > to_insert(2))
            k += 1;
        end;
        % now, k points to the position where to insert the new element:
        if(k==1)
            % insert at top:
            sorted = [to_insert; symbols_probabilities(1:end-1, :)];
        else
            % insert after k leading elements:
            sorted = [symbols_probabilities(1:k-1,:);...
                      to_insert;...
                      symbols_probabilities(k:end-1, :)];
        end;
    end;
end;


function codewords = fill_up(joined_symbols, codewords, root_index, up_bit, down_bit)
    index = joined_symbols(joined_symbols(:,1) == root_index, 2);
%    fprintf('Up (%d) => %d\n', root_index, index)
    codewords{index} = [codewords{root_index} up_bit];
    if(index > (joined_symbols(end,1)+1)/2)
        codewords = fill_up(joined_symbols, codewords, index, up_bit, down_bit);
        codewords = fill_down(joined_symbols, codewords, index, up_bit, down_bit);
    end
end;

function codewords = fill_down(joined_symbols, codewords, root_index, up_bit, down_bit)
    index = joined_symbols(joined_symbols(:,1) == root_index, 3);
%    fprintf('Down (%d) => %d \n', root_index, index)
    codewords{index} = [codewords{root_index} down_bit];
    if(index > (joined_symbols(end,1)+1)/2)
        codewords = fill_up(joined_symbols, codewords, index, up_bit, down_bit);
        codewords = fill_down(joined_symbols, codewords, index, up_bit, down_bit);
    end
end;


function [names probs] = parse_probs(prob_file)
	
	[values] = textread(prob_file, '%f');
	%disp(strings);
	disp(values);
	names=['Sym1' 'Sym2']';
	probs=values;
end;
