$nr_clusters = 71;
$nr_cell_types = 15;
$cell_sample_size = 500;
$max_comparison = 50;

open DATAFILE, 'Cell_stage_clique_umap_info.txt';
$i = 0;
while (<DATAFILE>) 
{
	chomp;
	@x = split '	',$_;
	if ($i==0){print "Nr columns = ",$#x+1,"\n";}
	else
	{
		for ($j=1;$j<=$#x+1;$j++) 
		{
			$datafile[$i][$j]=$x[$j-1];
		};
	};
	$i += 1;
};
close DATAFILE;
$nr_cells = $i-1;
print "Nr cells = ",$nr_cells,"\n";
$nr_ret_cells=0;
$nr_ips_cells=0;
for ($j=1;$j<=$nr_clusters;$j++)
{
	$clust_nr[$j][1]=0;
	$clust_nr[$j][2]=0;
};
for ($j=1;$j<=$nr_cell_types;$j++)
{
	$clique_nr[$j][1]=0;
	$clique_nr[$j][2]=0;
};
$nr_P0=0;
$nr_P5=0;
$nr_P9=0;
$nr_E13=0;
$nr_DD13=0;
$nr_DD21=0;
$nr_DD25=0;
$nr_DD29=0;


for ($i=1;$i<=$nr_cells;$i++)
{
	if ($datafile[$i][6] eq "Retinal")
	{
		$nr_ret_cells+=1;
		$source = 1;
	}
	elsif ($datafile[$i][6] eq "iPS")
	{
		$nr_ips_cells+=1;
		$source = 2;
	}
	else {print "Sample neither Ret nor IPS but ",$datafile[$i][6],"\n";sleep;};
	
	$found=0;
	$j=1;
	while ($j<=$nr_clusters && $found==0)
	{
		if ($datafile[$i][14] == $j-1){$clust_nr[$j][$source]+=1;$found=1;};
		$j+=1;
	};
	
	if ($datafile[$i][15] eq "Amacrine Cells"){$clique_nr[1][$source]+=1;}
	elsif ($datafile[$i][15] eq "Late RPCs"){$clique_nr[2][$source]+=1;}
	elsif ($datafile[$i][15] eq "Neurogenic"){$clique_nr[3][$source]+=1;}
	elsif ($datafile[$i][15] eq "PhotoPrecurs."){$clique_nr[4][$source]+=1;}
	elsif ($datafile[$i][15] eq "RNeurogenic"){$clique_nr[5][$source]+=1;}
	elsif ($datafile[$i][15] eq "Early RPCs"){$clique_nr[6][$source]+=1;}
	elsif ($datafile[$i][15] eq "NE"){$clique_nr[7][$source]+=1;}
	elsif ($datafile[$i][15] eq "RGCs"){$clique_nr[8][$source]+=1;}
	elsif ($datafile[$i][15] eq "Presumptive RPE"){$clique_nr[9][$source]+=1;}
	elsif ($datafile[$i][15] eq "Müller Glia"){$clique_nr[10][$source]+=1;}
	elsif ($datafile[$i][15] eq "Bipolar Cells"){$clique_nr[11][$source]+=1;}
	elsif ($datafile[$i][15] eq "Horizontal Cells"){$clique_nr[12][$source]+=1;}
	elsif ($datafile[$i][15] eq "Cones"){$clique_nr[13][$source]+=1;}
	elsif ($datafile[$i][15] eq "Rods"){$clique_nr[14][$source]+=1;}
	elsif ($datafile[$i][15] eq "NI"){$clique_nr[15][$source]+=1;}
	else {print "Add ",$datafile[$i][15],"\n";sleep;};
	
	if ($datafile[$i][5] eq "P0"){$nr_P0+=1;}
	elsif ($datafile[$i][5] eq "P5"){$nr_P5+=1;}
	elsif ($datafile[$i][5] eq "P9"){$nr_P9+=1;}
	elsif ($datafile[$i][5] eq "E13"){$nr_E13+=1;}
	elsif ($datafile[$i][5] eq "DD13"){$nr_DD13+=1;}
	elsif ($datafile[$i][5] eq "DD21"){$nr_DD21+=1;}
	elsif ($datafile[$i][5] eq "DD25"){$nr_DD25+=1;}
	elsif ($datafile[$i][5] eq "DD29"){$nr_DD29+=1;}
	else {print "Add ",$datafile[$i][5],"\n";sleep;};
};

print "Nr retinal cells = ",$nr_ret_cells,"\n";
print "Nr ips cells = ",$nr_ips_cells,"\n\n";
print "Nr ret cells in cluster ";
for ($i=1;$i<=$nr_clusters;$i++){print " ",$i,"-",$clust_nr[$i][1];};
print "\n";
print "Nr ips cells in cluster ";
for ($i=1;$i<=$nr_clusters;$i++){print " ",$i,"-",$clust_nr[$i][2];};
print "\n\n";
print "Nr ret cells in clique ";
for ($i=1;$i<=$nr_cell_types;$i++){print " ",$i,"-",$clique_nr[$i][1];};
print "\n";
print "Nr ips cells in clique ";
for ($i=1;$i<=$nr_cell_types;$i++){print " ",$i,"-",$clique_nr[$i][2];};
print "\n\n";
print "nr_P0 = ",$nr_P0," nr_P5 = ",$nr_P5," nr_P9 = ",$nr_P9," nr_E13 = ",$nr_E13,"\n";
print "nr_DD13 = ",$nr_DD13," nr_DD21 = ",$nr_DD21," nr_DD25 = ",$nr_DD25," nr_DD19 = ",$nr_DD29,"\n\n";

## Compute sampling rate for each experiment
if ($nr_P0 < $nr_DD13){$rate_P0 = 1; $rate_DD13 = $nr_P0/$nr_DD13;}
elsif ($nr_P0 == $nr_DD13){$rate_P0 = 1; $rate_DD13 = 1;}
elsif ($nr_P0 > $nr_DD13){$rate_P0 = $nr_DD13/$nr_P0; $rate_DD13 = 1;};

if ($nr_P5 < $nr_DD21){$rate_P5 = 1; $rate_DD21 = $nr_P5/$nr_DD21;}
elsif ($nr_P5 == $nr_DD21){$rate_P5 = 1; $rate_DD21 = 1;}
elsif ($nr_P5 > $nr_DD21){$rate_P5 = $nr_DD21/$nr_P5; $rate_DD21 = 1;};

if ($nr_P9 < $nr_DD25){$rate_P9 = 1; $rate_DD25 = $nr_P9/$nr_DD25;}
elsif ($nr_P9 == $nr_DD25){$rate_P9 = 1; $rate_DD25 = 1;}
elsif ($nr_P9 > $nr_DD25){$rate_P9 = $nr_DD25/$nr_P9; $rate_DD25 = 1;};

if ($nr_E13 < $nr_DD29){$rate_E13 = 1; $rate_DD29 = $nr_E13/$nr_DD29;}
elsif ($nr_E13 == $nr_DD29){$rate_E13 = 1; $rate_DD29 = 1;}
elsif ($nr_E13 > $nr_DD29){$rate_E13 = $nr_DD29/$nr_E13; $rate_DD29 = 1;};

open RET, '>>ret.txt';
open IPS, '>>ips.txt';
open SUM, '>>sum.txt';
for ($l=1;$l<=$max_comparison;$l++)
{
	$sum_average_nearest_dist_ret[0][$l]=0;
	$sum_average_nearest_dist_ips[0][$l]=0;
};
for ($i=1;$i<=$cell_sample_size;$i++)
{ ## Loop over reference cell pairs
print "Cell pair ",$i;
## Select RET and IPS cell randomly
$sel_ret_cell_nr = int(rand($nr_ret_cells))+1;
$j=0;
$counter=0;
while ($counter < $sel_ret_cell_nr)
{
		$j+=1;
		if ($datafile[$j][6] eq "Retinal"){$counter+=1;};
};
$sel_ret_cell_line = $j;
print ": ret cell ",$sel_ret_cell_nr," file line ",$sel_ret_cell_line; 

$sel_ips_cell_nr = int(rand($nr_ips_cells))+1;
$j=0;
$counter=0;
while ($counter < $sel_ips_cell_nr)
{
		$j+=1;
		if ($datafile[$j][6] eq "iPS"){$counter+=1;};
};
$sel_ips_cell_line = $j;
print " - ips cell ",$sel_ips_cell_nr," file line ",$sel_ips_cell_line,"\n"; 

## Select number-adjusted series of cells with which to measure distance, store distance
## RETINA
$retained_nr_cells=0;
for ($j=1;$j<=$nr_cells;$j++)
{
		if ($datafile[$j][6] eq "Retinal")
		{
			if ($datafile[$j][5] eq "P0")
			{
				$random = rand();
				if ($random < $rate_P0)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			}
			elsif ($datafile[$j][5] eq "P5")
			{
				$random = rand();
				if ($random < $rate_P5)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			}
			elsif ($datafile[$j][5] eq "P9")
			{
				$random = rand();
				if ($random < $rate_P9)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			}
			elsif ($datafile[$j][5] eq "E13")
			{
				$random = rand();
				if ($random < $rate_E13)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			};
		};
};
print "Nr of retained RET cells for comparison = ",$retained_nr_cells,"\n";
## Find the n nearest neighbors
$sum_nearest_distance_ret[$i][0]=0;
for ($l=1;$l<=$max_comparison;$l++)
{
	$present_shortest_dist = 1000;
	for ($j=1;$j<=$retained_nr_cells;$j++)
	{
		if ($line_nr[$j] != $sel_ret_cell_line)
		{
				if ($eucl_dist[$j] < $present_shortest_dist)
				{
					$present_shortest_dist = $eucl_dist[$j];
					$present_closest_cell = $j; 
				};
		};
	};
	$nearest_distance_ret[$i][$l] = $present_shortest_dist;
	$sum_nearest_distance_ret[$i][$l] = $sum_nearest_distance_ret[$i][$l-1]+$present_shortest_dist;
	$average_nearest_distance_ret[$i][$l] = $sum_nearest_distance_ret[$i][$l]/$l;
	$eucl_dist[$present_closest_cell] = 1000;
	
	$sum_average_nearest_dist_ret[$i][$l]=$sum_average_nearest_dist_ret[$i-1][$l]+$average_nearest_distance_ret[$i][$l];
	$average_average_nearest_dist_ret[$i][$l]=$sum_average_nearest_dist_ret[$i][$l]/$i;
};
print RET "Ref;",$i,";ND";
for ($j=1;$j<=$max_comparison;$j++)
{
	print RET ";",$nearest_distance_ret[$i][$j];	
};
print RET "\n";
print "SND:";
for ($j=1;$j<=$max_comparison;$j++)
{
	print " ",$sum_nearest_distance_ret[$i][$j];	
};
print "\n";
print "AND:";
for ($j=1;$j<=$max_comparison;$j++)
{
	print " ",$average_nearest_distance_ret[$i][$j];	
};
print "\n";
print "AAD:";
for ($j=1;$j<=$max_comparison;$j++)
{
	print " ",$average_average_nearest_dist_ret[$i][$j];	
};
print "\n";
print SUM "Ref;",$i,";RET";
for ($j=1;$j<=$max_comparison;$j++)
{
	print SUM ";",$average_average_nearest_dist_ret[$i][$j];	
};
print SUM "\n";



$retained_nr_cells=0;
for ($j=1;$j<=$nr_cells;$j++)
{
		if ($datafile[$j][6] eq "iPS")
		{
			if ($datafile[$j][5] eq "DD13")
			{
				$random = rand();
				if ($random < $rate_DD13)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			}
			elsif ($datafile[$j][5] eq "DD21")
			{
				$random = rand();
				if ($random < $rate_DD21)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			}
			elsif ($datafile[$j][5] eq "DD25")
			{
				$random = rand();
				if ($random < $rate_DD25)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			}
			elsif ($datafile[$j][5] eq "DD29")
			{
				$random = rand();
				if ($random < $rate_DD29)
				{
					$retained_nr_cells+=1;
					$line_nr[$retained_nr_cells]=$j;
					$eucl_dist[$retained_nr_cells]=sqrt(($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])*($datafile[$sel_ret_cell_line][16]-$datafile[$j][16])+($datafile[$sel_ret_cell_line][17]-$datafile[$j][17])*($datafile[$sel_ret_cell_line][17]-$datafile[$j][17]));
				};
			};
		};	
};
print "Nr of retained IPS cells for comparison = ",$retained_nr_cells,"\n";
## Find the n nearest neighbors
$sum_nearest_distance_ips[$i][0]=0;
for ($l=1;$l<=$max_comparison;$l++)
{
	$present_shortest_dist = 1000;
	for ($j=1;$j<=$retained_nr_cells;$j++)
	{
		if ($line_nr[$j] != $sel_ips_cell_line)
		{
				if ($eucl_dist[$j] < $present_shortest_dist)
				{
					$present_shortest_dist = $eucl_dist[$j];
					$present_closest_cell = $j; 
				};
		};
	};
	$nearest_distance_ips[$i][$l] = $present_shortest_dist;
	$sum_nearest_distance_ips[$i][$l] = $sum_nearest_distance_ips[$i][$l-1]+$present_shortest_dist;
	$average_nearest_distance_ips[$i][$l] = $sum_nearest_distance_ips[$i][$l]/$l;
	$eucl_dist[$present_closest_cell] = 1000;
	
	$sum_average_nearest_dist_ips[$i][$l]=$sum_average_nearest_dist_ips[$i-1][$l]+$average_nearest_distance_ips[$i][$l];
	$average_average_nearest_dist_ips[$i][$l]=$sum_average_nearest_dist_ips[$i][$l]/$i;
};
print IPS "Ref;",$i,";ND";
for ($j=1;$j<=$max_comparison;$j++)
{
	print IPS ";",$nearest_distance_ips[$i][$j];	
};
print IPS "\n";
print "SND:";
for ($j=1;$j<=$max_comparison;$j++)
{
	print " ",$sum_nearest_distance_ips[$i][$j];	
};
print "\n";
print "AND:";
for ($j=1;$j<=$max_comparison;$j++)
{
	print " ",$average_nearest_distance_ips[$i][$j];	
};
print "\n";
print SUM "Ref;",$i,";IPS";
for ($j=1;$j<=$max_comparison;$j++)
{
	print SUM ";",$average_average_nearest_dist_ips[$i][$j];	
};
print SUM "\n";

## IPS
## Compute average distance with growing number of neirest neighbors => write to file

## Update averages, print and write


}; ## Loop over reference cell pairs 

close RET;
close IPS;
close SUM;
open RESULTS, '>distr.txt';

print RESULTS "Cluster";
for ($i=1;$i<=$nr_clusters;$i++)
{
	print RESULTS ";",$i;  	
};
print RESULTS ";Total\n";
print RESULTS "Ret";
for ($i=1;$i<=$nr_clusters;$i++)
{
	print RESULTS ";",$clust_nr[$i][1];
};
print RESULTS ";",$nr_ret_cells,"\n";
print RESULTS "IPS";
for ($i=1;$i<=$nr_clusters;$i++)
{
	print RESULTS ";",$clust_nr[$i][2];
};
print RESULTS ";",$nr_ips_cells,"\n";
print RESULTS "PROP_IPS";
for ($i=1;$i<=$nr_clusters;$i++)
{
	#if (($clust_nr[$i][1]/$nr_ret_cells) > ($clust_nr[$i][2]/$nr_ips_cells))
	#{
		print RESULTS ";",($clust_nr[$i][2]/$nr_ips_cells)/(($clust_nr[$i][1]/$nr_ret_cells) + ($clust_nr[$i][2]/$nr_ips_cells));
	#}
	#else
	#{
	#	print RESULTS ";",($clust_nr[$i][1]/$nr_ret_cells)/(($clust_nr[$i][1]/$nr_ret_cells) + ($clust_nr[$i][2]/$nr_ips_cells));
	#}
};
print RESULTS "\n\n";

print RESULTS "Clique";
for ($i=1;$i<=$nr_cell_types;$i++)
{
	print RESULTS ";",$i;  	
};
print RESULTS ";Total\n";
print RESULTS "Ret";
for ($i=1;$i<=$nr_cell_types;$i++)
{
	print RESULTS ";",$clique_nr[$i][1];
};
print RESULTS ";",$nr_ret_cells,"\n";
print RESULTS "IPS";
for ($i=1;$i<=$nr_cell_types;$i++)
{
	print RESULTS ";",$clique_nr[$i][2];
};
print RESULTS ";",$nr_ips_cells,"\n";
print RESULTS "PROP_IPS";
for ($i=1;$i<=$nr_cell_types;$i++)
{
	#if (($clust_nr[$i][1]/$nr_ret_cells) > ($clust_nr[$i][2]/$nr_ips_cells))
	#{
		print RESULTS ";",($clique_nr[$i][2]/$nr_ips_cells)/(($clique_nr[$i][1]/$nr_ret_cells) + ($clique_nr[$i][2]/$nr_ips_cells));
	#}
	#else
	#{
	#	print RESULTS ";",($clust_nr[$i][1]/$nr_ret_cells)/(($clust_nr[$i][1]/$nr_ret_cells) + ($clust_nr[$i][2]/$nr_ips_cells));
	#}
};
print RESULTS "\n";



close RESULTS;