% Simulates irreversible edge tile binding
% B     B
%  B   B
%   B B
%    C
edge_tiles = 10000;
corner_tiles = 500;
simulations = 10;
corner_assembly_sizes = zeros(corner_tiles*2,simulations);
for i=1:simulations

  i
  corner_tile_left_sizes = zeros(1,corner_tiles);
  number_of_left_additions = zeros(1,corner_tiles);
  corner_tile_right_sizes = zeros(1,corner_tiles);
  number_of_right_additions = zeros(1,corner_tiles);

  edge_assemblies = ones(1,edge_tiles);
  number_of_edge_assemblies = edge_tiles;
  edge_bond = zeros(1,edge_tiles);

  percentage_completion = 0.95;
  circular_assemblies = 0;

  while number_of_edge_assemblies > ...
	((1-percentage_completion) * edge_tiles)
    % Find a random possible bond
    corner_bond_ratio = corner_tiles / ...
	(2 * corner_tiles + number_of_edge_assemblies * 2);
    corner_bond_ratios(number_of_edge_assemblies) = ...
	corner_bond_ratio * 2;

    bond = random ('Uniform', 0, 1);
    
    % And a random possible other bond:
    lucky_other_assembly = random('Discrete Uniform', ...
				  number_of_edge_assemblies);
    
    if bond < corner_bond_ratio
      % Bond is from a left corner tile: pick the assembly
      lucky_assembly = ceil(bond / corner_bond_ratio * corner_tiles);
      
      corner_tile_left_sizes(lucky_assembly) = ...
	  corner_tile_left_sizes(lucky_assembly) + ...
	  edge_assemblies(lucky_other_assembly);
      number_of_left_additions(lucky_assembly) = ...
	  number_of_left_additions(lucky_assembly)+1;
      edge_assemblies(lucky_other_assembly) = ...
	  edge_assemblies(number_of_edge_assemblies);
      number_of_edge_assemblies = number_of_edge_assemblies - 1;
      edge_bond(number_of_edge_assemblies) = 1;
    elseif bond < 2 * corner_bond_ratio
      % Bond is from a right corner tile
      lucky_assembly = ceil((bond - corner_bond_ratio)/corner_bond_ratio ...
			    * corner_tiles);
      
      corner_tile_right_sizes(lucky_assembly) = ...
	  corner_tile_right_sizes(lucky_assembly) + ...
	  edge_assemblies(lucky_other_assembly);
      number_of_right_additions(lucky_assembly) = ...
	  number_of_right_additions(lucky_assembly)+1;
      edge_assemblies(lucky_other_assembly) = ...
	  edge_assemblies(number_of_edge_assemblies);
      number_of_edge_assemblies = number_of_edge_assemblies - 1;
      edge_bond(number_of_edge_assemblies) = 1;
    else
      % Bond is between two edge assemblies
      lucky_assembly = ceil((bond - 2 * corner_bond_ratio) / ...
			    (1 - 2 * corner_bond_ratio) * ...
			    number_of_edge_assemblies);
      if lucky_assembly == lucky_other_assembly 
	'Circular assembly; skipping';
	circular_assemblies = circular_assemblies + 1;
      else
	edge_assemblies(lucky_assembly) = edge_assemblies(lucky_assembly) ...
	    + edge_assemblies(lucky_other_assembly);
	edge_assemblies(lucky_other_assembly) = ...
	    edge_assemblies(number_of_edge_assemblies);
	number_of_edge_assemblies = number_of_edge_assemblies - 1;
      end;
    end;
  end;
  
  corner_assembly_sizes(:,i) = [corner_tile_left_sizes, ...
		    corner_tile_right_sizes(1:corner_tiles)]';
end;
  
% Display the current assemblies
  
max_corner_size = max(max(corner_assembly_sizes));
min_corner_size = min(min(corner_assembly_sizes));
bins = 12

size(histc(corner_assembly_sizes(:,1),linspace(min_corner_size, ...
						  max_corner_size,bins)))

for i=1:simulations
  corner_assembly_sizes_hist(:,i) = histc(corner_assembly_sizes(:,i),linspace(min_corner_size, ...
						  max_corner_size,bins));
end;

corner_assembly_sizes_mean = mean(corner_assembly_sizes_hist,2);
corner_assembly_sizes_std = std(corner_assembly_sizes_hist,0,2);

figure(5);
axis([-5 125 0 50]);
hold off;
bar(linspace(min_corner_size,max_corner_size,bins),corner_assembly_sizes_mean,'y');
hold on;
errorbar(linspace(min_corner_size,max_corner_size,bins), ...
	 corner_assembly_sizes_mean,corner_assembly_sizes_std,'.');
hold off;
xlabel('Corner Assembly Length, in Number of Tiles','FontSize',14);
ylabel('Number of Assemblies','FontSize',14);
title(['Size Distribution, 95% Completed Irreversible Aggregation'],'FontSize',14);
saveas(gcf, 'irreversible_aggregate_incomplete_length_distribution_100_sims', 'fig');


figure(6);
% Arm histogram
arm_threshold = 5;
long_enough_arms = zeros(simulations,5);
i=1;
while (i < corner_tiles)
  for j=1:simulations
    number_of_long_arms = sum(corner_assembly_sizes(i:i+3,j) > arm_threshold);
    long_enough_arms(j,:) = long_enough_arms(j,:) + ((0:4) <= ...
						     number_of_long_arms);
  end;
  i=i+4;
end;

long_enough_arms_mean = mean(long_enough_arms,1);
long_enough_arms_std = std(long_enough_arms,0,1);


hold off;
bar(0:4,long_enough_arms_mean,'y');
hold on;
GRAPH = errorbar(0:4,long_enough_arms_mean,long_enough_arms_std,'.');

xlabel('Number of Attached Edges of Length >5','FontSize',14);
ylabel('Number of Assemblies','FontSize',14);
title(['X-Seeded Shapes,95% Completed Unseeded' ...
       ' Irreversible Aggregation'],'FontSize',14);
saveas(gcf, 'irreversible_aggregation_incomplete_long_enough_arms_100_sims', 'fig');
hold off;