#   This procedure computes a position for a new rectangle with given 
#   width and height that such that the weight of the covered area is 
#   minimised.
#   The call syntax is
#   \begin{quote}
#     |best_rectangle_position| \word{boundary} \word{old rectangles} 
#     \word{width} \word{height} \word{visible-var}\regopt
#   \end{quote}
#   and the return value is the pair $(x,y)$ of coordinates for the 
#   new value. The \word{boundary} is a rectangle in the form of a 
#   list
#   \begin{quote}
#     \word{min-x} \word{min-y} \word{max-x} \word{max-y} 
#   \end{quote}
#   within which the new rectangle should be positioned. The 
#   \word{old rectangles} is a list of lists with the structure
#   \begin{quote}
#     \word{min-x} \word{min-y} \word{max-x} \word{max-y} 
#     \word{weight}
#   \end{quote}
#   where the first four elements are integer coordinates of the 
#   sides of the rectangle, and the \word{weight} is a double that will 
#   be multiplied by the area of this rectangle. The order of the 
#   rectangles is assumed to be from bottom to top, so that later 
#   rectangles cover earlier ones. The \word{width} and \word{height} 
#   are of course the sides of the new rectangle to position.
#   
#   The \word{visible-var} is, if specified, 
#   the name of a variable in the calling context that will be set to 
#   the list of indices of those \word{old rectangles} which are not 
#   completely covered by other rectangles. A |-1| in this list means 
#   that some part of the window is not covered by any rectangle.
#   
#   From a strictly mathematical prespective, the problem considered 
#   here is to find a global minimum of ``covered area counted with 
#   weight'' seen as a function of position. This function is 
#   piecewise bilinear and thus possible to optimise through the 
#   general method of checking all critical points, all singular 
#   points, and all boundary points, but in view of (i) the 
#   bilinearity and (ii) the rather large number of singular points 
#   (typically, every boundary of a rectangle will be a set of 
#   singular points), it seems unlikely that there will be any 
#   minimum in a critical point that is not attained also in some 
#   boundary point. (Possibly this can even be proved, but I don't 
#   really have the patience to work it out right now.) Hence the 
#   procedure follows the simpler algorithm of locating all 
#   coordinate values where there may be a boundary and finds the 
#   minimum among these; maybe it isn't completely optimal, but I 
#   suspect it will be good enough.

proc best_rectangle_position {boundary oldL width height {visiblevar ""}} {
#   The first step is to find the $x$- and $y$-coordinates of points 
#   where there may be a boundary. This is (surprisingly?) simple.
   set xL [list [lindex $boundary 0] [lindex $boundary 2]]
   set yL [list [lindex $boundary 1] [lindex $boundary 3]]
   foreach rect $oldL {
      foreach {x1 y1 x2 y2 w} $rect break
      lappend xL $x1 $x2; lappend yL $y1 $y2
   }
   set xL [lsort -integer -unique $xL]
   set yL [lsort -integer -unique $yL]
#   The next step is to construct a description of the existing 
#   tiling. Since all boundaries are present in |xL| and |yL|, the 
#   cartesian product of these two contain not only all corners of 
#   the rectangles, but also all corners of the visible parts of the 
#   rectangles. This means |xL| and |yL| induce a partition of the 
#   container rectangle into smaller rectangles such that each of 
#   these are visible as part of precisely one of the rectangles from 
#   |oldL|. This means that the weights of different parts of the 
#   plane can be conveniently stored in a matrix.
#   
#   It is also convenient to precompute lists of the interval lengths 
#   in the $x$- and $y$-axis partitions.
   set dxL [list]; set x0 [lindex $xL 0]
   foreach x1 [lrange $xL 1 end] {
      lappend dxL [expr {$x1-$x0}]
      set x0 $x1
   }
   set dyL [list]; set y0 [lindex $yL 0]
   foreach y1 [lrange $yL 1 end] {
      lappend dyL [expr {$y1-$y0}]
      set y0 $y1
   }
#   In order to be consistent about the order of indices, let it be 
#   so that |lindex $weightM |$i$| |$j$ returns the weight associated 
#   with points $(x,y)$ such that |lindex $xL |$i$\({}\leqslant x 
#   \leqslant{}\)|lindex $xL |$i$|+1| and |lindex $yL |$j$\({}
#   \leqslant x \leqslant{}\)|lindex $yL |$j$|+1|.
   set weightM [lrepeat [llength $dxL] [lrepeat [llength $dyL] 0.0]]
   set visM [lrepeat [llength $dxL] [lrepeat [llength $dyL] -1]]
#   Now go through the list of rectangles and update the weight for 
#   those elements which are covered (here is thus where we rely on 
#   the stacking order). A trick that is used (and will be reused 
#   below) is to first compile lists of the index elements which are 
#   of interest here, and then |foreach| over these lists, rather 
#   than to use a |for|.
   set rcount 0
   foreach old_rect $oldL {
      foreach {x1 y1 x2 y2 weight} $old_rect break
      set iL [list]; set i 0
      foreach x $xL {
	 if {$x >= $x1 && $x < $x2} then {lappend iL $i}
	 incr i
      }
      set jL [list]; set j 0
      foreach y $yL {
	 if {$y >= $y1 && $y < $y2} then {lappend jL $j}
	 incr j
      }
      foreach i $iL {
	 foreach j $jL {
	    lset weightM $i $j $weight
	    lset visM $i $j $rcount
	 }
      }
      incr rcount
   }
#   As an optimisation, a second matrix giving the weighted areas for 
#   each little rectangle is computed as well. These will be needed 
#   comparatively often.
   set wgtaM [list]
   foreach column $weightM dx $dxL {
      set c2 [list]
      foreach w $column dy $dyL {
	 lappend c2 [expr {$w * $dx * $dy}]
      }
      lappend wgtaM $c2
   }
#   Now we get to the loop over all possible placements of the 
#   rectangle and the evaluation of which is best. An important part 
#   of this work is carried out by |interval_positions|, which is 
#   first used to find all possible positions on each axis. The parts 
#   of computing weighted areas can then be handled by two nested 
#   loops.
#   
#   In case several positions would cover the same minimal weighted 
#   area (which is not entirely unlikely if the rectangles being 
#   covered are large in comparison to the rectangle to position), 
#   then the area of the smallest rectangle with one corner at a 
#   boundary corner and containing the postitioned rectangle is used 
#   as a tie-breaker. This has the effect of choosing more peripheral 
#   positions whenever possible. The so far best weighted area is 
#   kept in |best|, whereas the associated peripherial area is kept 
#   in |best2|.
   set ipL [interval_positions $xL $width]
   set jpL [interval_positions $yL $height]
   set best infinity; set best2 infinity
   set midx [expr {
     ([lindex $boundary 2] - $width + [lindex $boundary 0]) / 2
   }]
   set midy [expr {
     ([lindex $boundary 3] - $height + [lindex $boundary 1]) / 2
   }]
   foreach ip $ipL {
      foreach jp $jpL {
	 set area 0
	 foreach i [lindex $ip 2] {
	    foreach j [lindex $jp 2] {
	       set area [expr {$area + [lindex $wgtaM $i $j]}]
	    }
	    foreach {j dy} [lindex $jp 3] {
	       set area [expr {$area + 
		 [lindex $weightM $i $j] * [lindex $dxL $i] * $dy
	       }]
	    }
	 }
	 foreach {i dx} [lindex $ip 3] {
	    foreach j [lindex $jp 2] {
	       set area [expr {$area + 
		 [lindex $weightM $i $j] * $dx * [lindex $dyL $j]
	       }]
	    }
	    foreach {j dy} [lindex $jp 3] {
	       set area [expr {$area + 
		 [lindex $weightM $i $j] * $dx * $dy
	       }]
	    }
	 }
	 if {$area <= $best} then {
	    set area2 [expr {
	      ( [lindex $ip 0]<$midx ? 
		[lindex $ip 1] - [lindex $boundary 0] :
		[lindex $boundary 2] - [lindex $ip 0] ) *
	      ( [lindex $jp 0]<$midy ? 
		[lindex $jp 1] - [lindex $boundary 1] :
		[lindex $boundary 3] - [lindex $jp 0] )
	    }]
	    if {$area < $best || $area==$best && $area2<$best2} then {
	       set ip0 $ip; set jp0 $jp
	       set best $area
	       set best2 $area2
	    }
	 }
      }
   }
#   Finally, if visibility information is requested then the 
#   positioned rectangle is imposed here and the result is returned.
   if {$visiblevar ne ""} then {
      foreach i [lindex $ip0 2] {
	 foreach j [lindex $jp0 2] {
	    lset visM $i $j [llength $oldL]
	 }
      }
      uplevel 1 [list ::set $visiblevar [
	 lsort -unique -integer [eval [list concat] $visM]
      ]]
   }
   return [list [lindex $ip0 0] [lindex $jp0 0]]
}

#   The |interval_positions| procedure is a helper for 
#   |best_rectangle_position| that determines all the intervals of a 
#   given length that have one endpoint in a given partition. This is 
#   the one-dimensional analogue of placing a rectangle with one 
#   corner in one of the given points, and is applied separately to 
#   the $x$- and $y$-coordinates.
#   
#   The call syntax is
#   \begin{quote}
#     |interval_positions| \word{partition} \word{length}
#   \end{quote}
#   where \word{partition} is a (sorted) list of points (numbers) 
#   in the partition and \word{length} is the length of the interval 
#   which should be positioned. The return value is the list of all 
#   possible positions. A position is a list
#   \begin{quote}
#     \word{lower endpoint} \word{upper endpoint}
#     \word{fully covered} \word{partially covered}
#   \end{quote}
#   where \word{lower endpoint} and \word{upper endpoint} are the 
#   endpoints of the interval at this position. 
#   \word{fully covered} is a list of indices (of intervals 
#   in the partition, i.e., indices into \word{partition}) for those 
#   intervals which are completely covered in this position, and 
#   \word{partially covered} is a list with the structure
#   \begin{quote}
#     \begin{regblock}[\regstar]\word{index} \word{width}\end{regblock}
#   \end{quote}
#   where each \word{index} is an index of a partially covered 
#   interval and the \word{width} is the width of the covered part of 
#   that interval. There are typically zero or one such pair in this 
#   list.
#   
#   The main loop simply goes through all the points of the 
#   partition, first looking backwards from and point and then 
#   forwards. In between, the index counter |n| of the interval in 
#   the partition that contains the current end of the positioned 
#   interval is incremented. Positions without a partial overlap 
#   interval are only added when looking forward.

proc interval_positions {partL length} {
   set min [lindex $partL 0]
   set max [lindex $partL end]
   set res [list]
   set n -1
   foreach p $partL {
      if {$p - $length >= $min} then {
	 set fullL [list]
	 set c_prev 0
	 for {set m $n} {[
	    set covered [expr {$p - [lindex $partL $m]}]
	 ] < $length} {incr m -1} {
	    lappend fullL $m
	    set c_prev $covered
	 }
	 if {$covered > $length} then {
	    lappend res [list [expr {$p-$length}] $p $fullL\
	      [list $m [expr {$length - $c_prev}]]]
	 }
      }
      incr n
      if {$p + $length <= $max} then {
	 set fullL [list]
	 set covered 0
	 for {set m $n}\
	   {[set q [lindex $partL [expr {$m+1}]]] < $p + $length}\
	   {incr m} {
	    lappend fullL $m
	    set covered [expr {$q - $p}]
	 }
	 lappend res [list $p [expr {$p+$length}] {} {}]
	 if {$q == $p + $length} then {
	    lappend fullL $m
	 } else {
	    lset res end 3 [list $m [expr {$length - $covered}]]
	 }
	 lset res end 2 $fullL
      }
   }
   set res
}