use application 'polytope'; # Find the number of lattice points for the 3-leaf tree. use vars '\$s', '\$t', '\$B', '\$m', '\$sm', '\$ineq'; use vars '\$fix', '\$i', '\$j', '\$k', '\$l','\$eq', '\$a', '\$point', '\$counter', '@points'; # Define the polytope associated to the Kimura 3-parameter model on the 3-leaf tree. \$s=new Polytope(POINTS=><<"."); 1 1 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 . # Apply the vertex lattice normalization to \$s. \$t=vertex_lattice_normalization(\$s); # Basis transformation matrix taking vertices of \$s to vertices of \$t. \$B=new Matrix(<<"."); 1 -3 -3 -2 -2 -1 -1 -3 -1 0 0 3/4 3/2 1/2 1 1/4 1/2 3/2 1/2 0 0 3/4 1/2 1/2 1 1/4 1/2 1/2 0 0 0 3/4 1/2 1/2 0 1/4 1/2 1/2 1/2 -1/2 0 3/4 1/2 1/2 0 1/4 -1/2 1/2 0 1/2 0 3/2 3/4 1 1/2 1/2 1/4 0 0 0 0 1/2 3/4 1 1/2 1/2 1/4 1 1/2 0 0 1/2 3/4 0 1/2 1/2 1/4 1 0 1/2 0 1/2 3/4 0 1/2 -1/2 1/4 1 1/2 -1/2 0 3/4 3/4 1/2 1/2 1/4 1/4 3/2 1/2 0 0 3/4 3/4 1/2 1/2 1/4 1/4 1/2 0 0 0 3/4 3/4 1/2 1/2 1/4 1/4 1/2 1/2 -1/2 0 3/4 3/4 1/2 1/2 1/4 1/4 1/2 0 1/2 . # Matrix \$m for multiplying the vertices of \$s by 3. \$m=new Matrix(<<"."); 1 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 3 . # Define the polytope \$sm as the convex hull of the vertices of \$s multiplied by 3. \$sm=new Polytope(VERTICES=>(\$s->VERTICES)*\$m); # \$ineq contains facet inequalities for \$sm. \$ineq=(\$sm->FACETS); # Loop that finds the number of lattice points for the polytope 3*\$s with the first 4 coordinates fixed. \$counter=0; # The 2nd to 5th coordinate of \$sm can take integer values # from 0 to 3, which sum up to 3. for \$i (0...3){ for \$j (0...(3-\$i)){ for \$k (0...(3-\$i-\$j)){ \$l=3-\$i-\$j-\$k; # Matrix \$fix for equalities that fix 2nd to 5th coordinate. \$fix=new Matrix(<<"."); -\$i 1 0 0 0 0 0 0 0 0 0 0 0 -\$j 0 1 0 0 0 0 0 0 0 0 0 0 -\$k 0 0 1 0 0 0 0 0 0 0 0 0 -\$l 0 0 0 1 0 0 0 0 0 0 0 0 . # Define a new set of equalities by taking the affine hull of \$sm and \$fix. \$eq=(\$sm->AFFINE_HULL)/\$fix; # Define a new polytope \$a using the new set of equalities \$eq and the inequalities of \$sm. \$a=new Polytope(INEQUALITIES=>\$ineq,EQUATIONS=>\$eq); # Apply the basis transformation matrix to \$a and count lattice points. \$a=new Polytope(POINTS=>(\$a->VERTICES)*\$B); \$point=(\$a->N_LATTICE_POINTS); # Remember all the results in the array \$points. \$points[\$counter]=\$point; \$counter++; } } } #Find the number of lattice points for 4-leaf tree. use vars '\$p', '\$r', '\$A', '\$m', '\$pm', '\$ineq'; use vars '\$i', '\$j', '\$k', '\$l', '\$w', '\$x', '\$y', '\$z', '\$e', '\$f', '\$g', '\$h'; use vars '\$fix', '\$eq', '\$a', '\$point', '\$sum1', '\$sum2'; use vars '\$counter1', '\$counter2', '@points1', '@points2'; # Define the polytope of the Kimura 3-parameter model on the 4-leaf tree. \$p=new Polytope(POINTS=><<"."); 1 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 1 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 1 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 1 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 . # Apply vertex lattice normalization to \$p. \$r=vertex_lattice_normalization(\$p); # Basis transformation matrix taking vertices of \$p to vertices of \$r. \$A=new Matrix(<<"."); 1 3 -1 0 3 1 -1 1 3 -3 -1 1 -1 -1 -1 -2 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 0 0 0 0 0 0 0 1/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -1 1 1 0 0 0 1/2 3/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -1 1 0 0 1/2 1/2 0 1/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -1 1 0 -1 1/2 1/2 1/2 3/4 0 -3/4 1/4 0 0 -1/4 0 0 -3/4 0 1/4 -1/4 0 0 0 1/4 0 -3/4 1/4 0 -1 -1/4 1 0 -3/4 1 1/4 -1/4 0 0 1/2 3/4 0 -3/4 1/4 0 -1 -1/4 0 0 -3/4 1 1/4 -1/4 1/2 1/2 0 1/4 0 -3/4 1/4 0 -1 -1/4 0 -1 -3/4 1 1/4 -1/4 1/2 1/2 1/2 3/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -3/4 3/2 1/4 -1/4 1/2 1/4 1/2 1/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -3/4 1/2 1/4 -1/4 1/2 1/4 0 1/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -3/4 1/2 1/4 -1/4 0 1/4 1/2 3/4 0 -3/4 1/4 0 -3/4 -1/4 1/4 -1/4 -3/4 1/2 1/4 -1/4 0 1/4 0 3/4 0 0 0 -3/4 -3/4 0 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 1/4 1/4 1/4 0 -1 1 1/4 -3/4 0 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 1/4 1/4 3/4 0 -1 0 1/4 -3/4 0 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 3/4 1/4 3/4 0 -1 0 1/4 -3/4 -1 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 -1/4 1/4 1/4 0 -3/4 1/4 3/4 -3/4 -1/4 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 1/2 1/4 1 0 -3/4 1/4 -1/4 -3/4 -1/4 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 1/2 1/4 1/2 0 -3/4 1/4 -1/4 -3/4 -1/4 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 0 1/4 1/2 0 -3/4 1/4 -1/4 -3/4 -1/4 1/4 -1/4 -3/4 3/4 1/4 -1/4 1/4 0 1/4 0 . # Matrix \$m for multiplying the vertices of \$p by 3. \$m=new Matrix(<<"."); 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 . # Define the polytope \$pm as the convex hull of the vertices of \$p multiplied by 3. \$pm=new Polytope(VERTICES=>(\$p->VERTICES)*\$m); # \$ineq contains facet inequalities for \$pm. \$ineq=(\$pm->FACETS); # Loop that finds the number of lattice points for polytope # 3*\$p with 14th to 21th coordinate fixed (corresponding to 2 # incident leaf edges) and with 14th to 17th coordinate fixed # (corresponding to 1 leaf edge). \$counter1=0; # The 14th to 21th coordinate of 3*\$s can take integer values # from 0 to 3, coordinates 14 to 17 sum up to 3 and # coordinates 18 to 21 sum up to 3. for \$i (0...3){ for \$j (0...(3-\$i)){ for \$k (0...(3-\$i-\$j)){ \$sum1=0; \$counter2=0; for \$w (0...3){ for \$x (0...(3-\$w)){ for \$y (0...(3-\$w-\$x)){ # Loop that fixes 2nd to 5th coordinate corresponding to an # additional leaf edge to make computations faster. We find # the number of lattice points separately for each possible # integer fixation and then sum over all the values. \$sum2=0; for \$e (0...3){ for \$f (0...(3-\$e)){ for \$g (0...(3-\$e-\$f)){ \$l=3-\$i-\$j-\$k; \$z=3-\$w-\$x-\$y; \$h=3-\$e-\$f-\$g; # Matrix for equalities that fix coordinates under # consideration. \$fix=new Matrix(<<"."); -\$i 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 -\$j 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 -\$k 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 -\$l 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 -\$w 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 -\$x 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -\$y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -\$z 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -\$e 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -\$f 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -\$g 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -\$h 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . # Define a new set of equalities by taking the affine hull of # \$pm and \$fix. \$eq=(\$pm->AFFINE_HULL)/\$fix; # Define a new polytope \$a using the new set of equalities \$eq and # the inequalities of \$pm. \$a=new Polytope(INEQUALITIES=>\$ineq,EQUATIONS=>\$eq); # Apply the basis transformation matrix to \$a and count # lattice points. \$a=new Polytope(POINTS=>(\$a->VERTICES)*\$A); \$point=\$a->N_LATTICE_POINTS; # Sum the number of lattice points over all possible integer # values corresponding to the additional leaf edge. \$sum2=\$sum2+\$point; } } } # Array @points2 holds the number of lattice points when 14th # to 21st coordinate are fixed. \$points2[(\$counter1*20)+\$counter2]=\$sum2; \$sum1=\$sum1+\$sum2; \$counter2++; } } } # Array @points1 holds the number of lattice points when 14th # to 17th coordinate are fixed. \$points1[\$counter1]=\$sum1; \$counter1++; } } } use vars '\$sum'; # Finding the total number of lattice points by using the # formulas (there are 20 possibilities for integer values # corresponding to one edge). \$sum=0; for \$i(0...19){ for \$j(0...19){ \$sum=\$sum+(\$points2[(\$i*20)+\$j]*\$points[\$i]*\$points[\$j]); } } print "Snowflake has \$sum lattice points. \n \n"; \$sum=0; for \$i(0...19){ \$sum=\$sum+(\$points1[\$i]*\$points1[\$i]); } print "3-caterpillar has \$sum lattice points. \n";