This problem can be solved with straight-forward dynamic programming. Faced with building this desired shape, consider any solution. Either there's a solution that includes a long strip that covers the whole length of the bottom (then the top becomes another subproblem):

             ________                         ________
            |        |____                   |        |____    
      ______|             |            ______|             |   
     |                    |           |                    |   
 ____|                    |  =>   ____|    subproblem      | 
|                         |      |_________________________|   
|                         |      |     a single strip      |   
|_________________________|      |_________________________|
[download]

             ________                         ________
            |        |____                   |        |____    
      ______|             |            ______|             |   
     |                    |           |      |             |   
 ____|                    |  =>   ____|      |             | 
|                         |      |           | subproblem  |   
|                         |      | subproblem|             |   
|_________________________|      |___________|_____________|
[download]

             ________                            ________
            |        |____                      |________|____    
      ______|             |               ______|             |   
     |                    |  NEVER!      |______|_____________|   
 ____|                    |    =>    ____|____________________| 
|                         |         |___________|             |   
|                         |         |           |_____________|   
|_________________________|         |___________|_____________|
[download]

use List::Util qw(min sum);
use strict;
use Memoize;

memoize 'solve';

my @thicknesses = (4, 8, 16, 24, 32, 47, 64);
my @goal = ( [10,20], [15,28], [20,48], [10,40] );

sub argmin(&@) {
    my $i = undef;
    my $min = undef;
    my $block = shift;
    for (@_) {
         my $this = $block->($_);
         if ( not defined $min or $min > $this ) {
             $min = $this;
             $i = $_;
         }
    }
    return wantarray ? ($i, $min) : $i;
}

my @decomp;
my $INFINITY = 1e10;

sub solve {
    my ($left, $right, $height) = @_;

    ## we could either we can split into two pyramids..

    my ($split_index, $split_value) = (undef, $INFINITY);
    if ($left != $right) {
        ($split_index, $split_value) =
            argmin { solve($left, $_, $height)
                      + solve($_+1, $right, $height) } $left .. $right
+-1;
    }


    ## ... OR we could try to place a strip across the whole length

    my $min_height = min map { $_->[1] - $height } @goal[$left .. $rig
+ht];
    my @strips = grep { $_ <= $min_height } @thicknesses;

    ## (degenerate case, nothing left to do!)
    return 0 if $min_height == 0 and $left == $right;

    my ($strip_size, $strip_value) = (undef, $INFINITY);
    if (@strips) {
        ($strip_size, $strip_value) =
            argmin { 1 + solve($left, $right, $height+$_) } @strips;
    }

    ## which one of these is best?

    if ($strip_value < $split_value) {
        $decomp[$left][$right][$height] =
            [ strip => sum( map { $goal[$_][0] } $left .. $right ),
                       $strip_size ];
        return $strip_value;
    } else {
        $decomp[$left][$right][$height] = [ split => $split_index ];
        return $split_value;
    }
}   

printf "you can solve this with %d strips\n",
    solve(0, $#goal, 0);

reconstruct(0, $#goal, 0);

sub reconstruct {
    my ($left, $right, $height) = @_;
    my $how = $decomp[$left][$right][$height];

    return if not ref $how;
    if ($how->[0] eq "strip") {
        print "strip of length $how->[1], height $how->[2]\n";
        reconstruct($left, $right, $height + $how->[2]);
    } else {
        my $split = $how->[1];
        reconstruct($left, $split, $height);
        reconstruct($split+1, $right, $height);
    }
}
[download]

you can solve this with 5 strips
strip of length 55, height 16
strip of length 25, height 4
strip of length 15, height 8
strip of length 20, height 32
strip of length 10, height 24
[download]

Again, this only is correct (well, I hope it's correct) when you can build the shape exactly. It's nontrivial, but not too difficult, to modify this algorithm to do approximate matching, once you specify what you want to optimize. But I think I've already spent enough time on this for one afternoon ;)

You may be able to convince yourself, but I'm convinced that you're wrong about that decomposition being possible. Let me try to convince you.

Consider the case where you have blocks with width/thickness pairs of 4,3, 1,5 and 3,2. They can be arranged into the following pattern to form a 5,8 square:

 _______ _
|       | |
|       | |
|_______| |
| |     | |
| |_____|_|
| |       |
| |       |
|_|_______|
[download]

(I've actually drawn each horizontal space as 2 characters...)

I assert that that example can't be broken down as you claimed it could. A brute force proof is fairly straightforward. You have to put some block in the top left corner. Let's examine each possibility.

• Top left is a 4,3 block. Then going clockwise looking at width and depth it is easy to show you need a 1,5, then a 4,3, then a 1,5 leaving the 3,2 hole in the middle. (This is the diagram above.)
• Top left is a 1,5 block. Going counter-clockwise it is easy to see that you need a 4,3 then 1,5 then 4,3 leaving a 3,2 hole in the middle. (This is the mirror image of the diagram above.)
• Top left is a 3,2 block. Going clockwise we need it 2 wider, so we have to add 2 1,5s. To add 3 to the depth below them you need a 4,3, and then to the left of that we have to add 1 more width so we need a 1,5 that unfortunately intersects our original 3,2 block. This case is therefore impossible.
So you see that by brute force there are 2 ways to solve this problem, and neither decomposes as you'd hoped.

Thanks Blockhead,
I will take a look at this and see what I can do.

#!/usr/bin/perl

# First line of DATA is desired total and the tolerance.

my $Finish = <DATA>;
my ($Wanted, $Tolerance) = ($Finish =~ /(-?\d+)\s+(\d+)/);

# Other lines of DATA are the available elements.

my @Items = (<DATA>); chomp @Items;

#---------------------------------------------------------------------
+-------

my $Stack = ();

sub build
{
    my $index = shift;
    my $total = shift;

    while ($index < @Items)
    {
        my $num = $Items[$index];
        push(@Stack, $num);

        if ($num + $total > ($Wanted+$Tolerance))
        {
            next;
        }
        elsif ($num + $total >= ($Wanted-$Tolerance) and
               $num + $total <= ($Wanted+$Tolerance))
        {
            print "( @Stack ) = ", $num+$total, "\n";
        }
        else
        {
            build($index+1, $num + $total);
        }
    }
    continue
    {
        pop(@Stack);
        $index++;
    }
}

#---------------------------------------------------------------------
+-------

build(0, 0);

__DATA__
30 5
-5
5
10
15
20
10
25
25
[download]

I am trying to arrive at a solution of building irregular 2d "pyrimids" or shapes.

To follow up blockhead's advice about dynamic programming, you might also want to look at greedy algorithms.

A greedy algorithm breaks a problem down into recursive subproblems just like dynamic programming does, but it picks the best solution at each level before moving on to the subproblem. It's like making change for 47c. A quarter is the largest coin that fits in a 47c 'hole'. That leaves 22c. A dime is the largest coin that fits in a 22c 'hole'. That leaves 12c. By the same process, we get another dime and two pennies.

Greedy algorithms perform well on average, but they do not guarantee an optimal solution every time. If we tried the above experiment with coins worth 35c, 30c, and 15c, the greedy algorithm would tell us that a single 35c coin is the best match. But in fact, a 30-15 combination would get us closer to 47c.

If you want to guarantee an optimal solution, you need to use dynamic programmming. If you're willing to tolerate imperfect solutions, or you know that your pieces are partitioned well enough to solve any problem, you can use greedy algorithms.

Thanks mstone, Looking at the discreet thicknesses: 0.4, 0.8, 1.6, 2.4, 3.2, 4.7, 6.4, All of these (save 4.7) are divisible by 0.4, so a greedy algorythm may work just fine.

Yes some example code would be helpful. From the outside looking in, it sounds like you're trying to solve a 2d packing problem.

Or trying to win at tetris...

Well, I am at the brainstorming stage right now, I was hoping that someone clever has seen a module or some code that does something similar, or has a clever suggestion.
With that said my plans are to:


• Parse the input to seperate the groups (apply some rounding rules to find the 0 thicknesses)
• Starting with the base (first level), build that level by iterating over the discreet thicknesses, keeping the closest match.
• Subtract the thickness of the first level from each of that groups thicknesses and re-apply the rounding rule to further breack down into groups.
• Repeat this process until all thicknesses are found.


The problem with this is that it does not determing the best overall thicknesses, the latter thicknesses would be further limited by the early choices, so I would have to further determine possiblities by including a minimum agreement thickness and then iterate over all possibilities assigning a score to each iteration and getting a solution based off of scores.

This is very brute force, but it is my first thought.

I guess I don't understand the correlation between your discreet thicknesses and the pairs entered. How does

1,2
1.5,3
[download]

correlate to the descreet thicknesses?

0.4, 0.8, 1.6, 2.4, 3.2, 4.7, 6.4
[download]

Also, why is the pair 3,0 different? Is it a delimiter?

To clear it up, the first number is the width that you want, the second number is the thickness that you want. You can't always get what you want, but it should be close (that is where the discreet thicknesses of the building blocks comes in).

Lets say that you wanted a width/thickness of 5,0.2, well there is no constraint on width, but there is a 0.4 minimum thickness for the building blocks, so you would have to decide if you would span the 5 width (this could be cm,mm...) with 0 thickness (rounded down) or 0.4 thickness(rounded up).

The pair 3,0 is different because it has 0 thickness.
So if you were builing a 2 humped shape and had a 0 thickness in the middle, you could then look at it as 2 seperate shapes, but if it was a small thickness, say 0.4 in the middle, then it would be 1 complete shape.

I hope that this clears things up.

I too am pretty well lost. I think you've attempted to simplify your description of the problem too much, as it's fairly obvious to me, there are many unstated rules that must be followed to get from the input to the output that is desired.

How do you determine where a given block can start on a given row, how do you know when a row is finished? Etc.

-Scott

Thank you Scott,

I did try to just get to the basics of what I need and left out some important info. Here is the totality of what I am trying to accomplish:

I am a Medical Physicist (deals with treating cancer with radiation). When we irradiate the entire body (this is done to supress the immune system before a transplant, i.e bone marrow), we want the entire body to receive uniform dose. If you look at the body from the side, and consider its verying thicknesses (neck, shoulder, hips, knees,...), it becomes clear that we need to attenuate some of the radiation at different areas. For example, at the midline of the patient, the radiation has to travel though more tissue at the shoulders then it does at the neck, so we place a lead filter of a certain thickness for only the radiation that will go though the neck so that the dose at the midline of the patient will be the same. Ideally we want to build up these filters, if you just place the filters next to eachother, then there is a chance that radiation can escape unfiltered between the filters leading to hotspots.

Now, when filling out the info the width/thickness of the filters would start from the top of the patient and would continue in series towards the patients feet.

i.e. 3,2
5,4

Would indicate that the first filter would be 3cm wide and 2mm thick and the second filter would be 5cm wide (starting after the first filter) and 4mm thick.

This would be built using an 8cm piece of 1.6mm lead, an 8cm piece of 0.4mm lead, a 5cm piece of 1.6mm lead, and a 5cm piece of 0.4mm lead and would look like the following:



             ____________________
            |____________________| 0.4mm
            |                    |
 ___________|____________________| 1.6mm
|________________________________| 0.4mm
|                                |
|________________________________| 1.6mm
[download]

Let me know if I can clear anything up on this.

Thanks again everyone!

I am not familiar with Prolog, but I will do some reading on it. Hopefully it will be something new that I can add to my arsenal.

Well I only saw those metal masks you put where the beam hits, where you shape the hole to fit the shape of the tumor, this is totally different from what you are doing I believe.

I imagine if you look down on a patient lying on his back, if the radiation is coming from above, you basically would want something like a topographical map (like a 3D contour map). See this map to get an idea of what I mean. Mountain peaks are where it gets a few layers thick. Surely this would be optimum? You would use a single thickness lead sheet and cut it to make each layer of the map. Then you also wouldn't have vertical seams that could maybe burn the patient? There would be no packing problem then either.

That said, the (brass?) masks I saw did indeed seem to have been created by a computerized cutter, as the holes were made of many fine vertical cuts, perhaps <1mm in thickness and maybe more resolution in the y axis.

So my questions are:

1. Is this really the best solution? I mean you have to type in all those numbers, wouldn't it be better if you could go directly from some scanned data to a lead solution?

2. The other question is, this is all seems to be 2D, i.e. looking from the side at a cross section. But we're really talking about 3D, that is, an arrangement of many slabs all over the patient, more slabs in some places than others.

3. I feel pretty uncomfortable about how these threads tend to die quickly but you might need more help. I honestly don't think it's appropriate to be asking such important questions on this kind of a board when misunderstanding of a problem could jepoardize someone's life (though the basic ideas of packing and dynamic programming are of course good things to know). For example in the beginning I thought it was a student with a homework problem!
Maybe it's a problem with perlmonks that things are so evanescent (I still don't know how to view the list of nodes that showed on the homepage 2 days before). In particular I was wondering about whether vertical seams are dangerous, since it made me worry about a completely unattenuated beam burning a line into the patient if for example the pattern has a completely vertical seam, then how it is manufactured, i.e. could it come apart there, etc. That is how I came up with the idea of cutting each layer of a contour map with a jigsaw so there would be no seams and it could closely match the results of a CT scan, allowing much higher resolution than if you just type in the numbers (possibly making mistakes) yourself.

Every once in a while heavy particle radiation makes it back in the news, but it has lost favor with many hospitals because of the size of the accelerators. Most accelerators that are currently being used emit photons and/or electrons. The heavy particle accelerators emit either protons or neutrons. An electron is ~ 1000 times lighter than protons or neutrons, so they are much easier to accelerate. neutrons due not have charge, so the way that neutron beams are created is by colliding protons into the (although we are dealing with subatomic particles and quantum physics rule this domain, general mechanics of equal mass collisions are still a good approximate). Getting back to the size issue, an electron (linear) accelerator can fit in a room easily, but a proton accelerator is the size of a small building.

The metal masks that you speak of are called blocks and are usually made of cerrobend.


Now to answer your questions:

1. The setup and treatment of total body irradiation is more complicated than taking a CT and building the filters that are required. BUT on a theoretic basis, yes you are correct, a continuous distribution would be best. BUT we do not live in a theoretical world. We have to look at complicity and cost v.s. gains. For TBIs, we are not treating cancer, we are trying to suppress the immune system to prepare the body for a transplant. It is not crittical to get absolute even dose over the entire body +/- 10% works fine.


2. I think that #1 answers this, we are looking to correct gross changes. Also the beam enters the patient laterally (through the side), so imaging that you could see that it does become a 2d problem.


3. Yes, these thread can die quickly, but then again there are threads that seem to never end. It all depends on the intrest level, and the area of the question. As far as burning the patient, we take every precaution to ensure that the patient is not mistreated. Currently this involve painstaking measures of creating the filters 4-6 hours per patient. We also take extreme care in the setup so that there are not any light leaks. That being said, we are changing the way that we make the blocks from a column type to a pyrimid type, and I am trying to write some software (as a tool) to quickly figure out the best way to do this. It is easy to do by hand, but can be time consumming. The software will only be a tool, and every thing that we do will be checked and double checked.



Thank you everyone for your ideas and comments!

Thank you very much for your detailed comments and very interesting project.

Yes, the site I visited was as you discuss, a medium sized building which cost $300 million to build, very expensive. Apparently using heavy particles (actually not just a proton but a whole carbon atom I believe) works well for deep targets but you know all this. It was very interesting though it was wierd to see even small children following the route which requires you to climb over and under the beam tube. Though the staff did not appreciate my asking (quietly to the side) if you would die from lethal x-ray dose etc. if you are in the donut shaped corridor holding the beam when it is turned on.. apparently it is accelerated to 80% of the speed of light.

About the dose entering from the side it does seem (superficially to me sorry and I know nothing of course of diffraction in the body but..) there would be a shadow from bones or at least more radiation absorbed on the side of the body the radiation comes from? That is, I thought it is absorbed within a foot of the surface. Oh well sorry to bug you with silly questions, the packing and pyramid shaping ideas are quite interesting!

I have not used it in large problems but just for fun there is Damian Conway's module series Quantum::Superpositions which provides the "any" and "all" words. It seems that some problems can be more easily stated in those terms, but really the best seems to me to try and make lots of blocks by hand and then translate some of the procedures you develop as you go along into bits of the program. Also possibly dynamic programming (blokhead mentioned this above it seems..) such as that used to find alignments of genetic sequences, where you solve many sub-problems first based on a set of scoring rules, then find a path through the solutions, to find the most efficient solution.

The above approach is different from the way I was thinking about it, but the part I should stress is that depending on the scoring rules you create, a dynamic programming approach can create different outcomes as to what is the most efficient. Obviously, you might think, but it can be subtle. For example there is the Smith-Watermann method in genomics. You can end up with more sensitivity, or more preference for a certain type of solution, depending on whether you want to penalize certain kinds of mutations, gaps, or end features. In these cases it seems best to experiment with variations, and also consider statistical/practical issues (for example statistically where beam is coming from or what shapes are more dangerous, or what shapes are easier to build).

If you are aligning two genetic sequences with dynamic programming, you end up drawing a matrix with one sequence (of letters A,T,G or C) for each axis. You start in the upper left corner and go right, then go to the next line when done. Basically you use the scoring rules to decide what score to write into a given cell, based on its coordinates and the contents of nearby (i.e. the previous) cell. When done filling in the whole matrix, you then go to the lower right hand corner and using another rule basically draw arrows towards the upper left corner, seeking the highest scoring valid cell. Basically this is supposed to give you the best sequence alignment and can be done by hand. So dynamic programming would suggest a similar kind of a procedure for your problem.

Also it appears to me that we really want to solve this for the whole length of the body, in other words not just a single pyramid but a number of pyramids running from feet to head. So continuity will be important as a way to reduce cost (well unless you have only short pieces available) and to reduce seams. However it will become more computationally complex to solve for the whole body, and very hard to test a solution by eyeball to see if it makes sense. I wonder if I understand this correctly. If so, the seam between two pyramid units should be considered as a penalty and you have to have a piece that bridges the two.

My own guess is that I would penalize cases in which vertical seams are created at all (since a single piece makes it easier to manufacture), and highly penalize a double seam (since a lot of radiation could leak through to sensitive tissue), if I am interpreting the pyramid to mean the beam is coming from the top of the pyramid, i.e. the base of the pyramid is actually against the patient's body, sideways, with the beam coming horizontally into the body from the top of the pyramid in through the base of the pyramid and then into the patient. In which case a vertical seam is bad, which is why you want to stagger or "build up" the pyramid like bricks.

In this case the drawing tilly provided might not be selected, depending on the rules, because it has lots of vertical seams, both dangerous radiation wise as well as perhaps using too many strips (we are ignoring how many thick pieces you have available).

Also you need to penalize or disallow solutions that are impossible to build (if they are?) such as a pyramid with a hole in it. To me it seems that you should be inputting not the widths and thicknesses of metal pieces, but the desired beam density at each coordinate, and a way to determine the beam density you get from a given solution, if that is what you are trying to determine.

Anyway, the most important thing to do is to determine how to score potential solutions.

Anyway, if you need something to get the juices flowing, here are some pages to look at: algorithm design and Design Patterns in Dynamic Programming. I just googled for "how to" and "dynamic programming" but also noticed that there are lots of kinds of "heuristic programming" methods. Also dynamic programming can take a long time for large data sets, maybe this is not a problem of yours.

Finally, and this is not really my field, but there is a whole field of bin packing, set covering and partioning problems. One famous problem is the "knapsack problem". There are different algorithms that go by names like linear programming, simulated annealing, continuous relaxation and others. I'm not qualified to give you an overview of this sort of thing and my guess is that you don't have enough data to be worrying so much about the search algorithm, you should pay most attention to the definition of the problem.

By the way I did find (google for radiation and "dynamic programming") a reference to computation of radiation doses PDF and as HTML from google. These people talk about something like what I hypothesized, using the CT scan to compute doses, but say it is too big a problem computationally and so use a heuristic called a Markov decision problem, or MDP. Anyway, google says that variants of dynamic programming are indeed used widely for radiation planning over time, though this is possibly a different problem from yours.

Hope this thread helps you in working out this problem.