Bayes 8: impossible

By thins point I'm done with re-telling the things I've read, and we're going into the things I've discovered in practice, and the solutions that I've come up with. They aren't always perfect, and often I have more than one solution with different pros and cons, and you might be able to find some better ones, but you can't deal with the real data without some sort of these solutions.

Let's look at another example fed into the same table as before:

evC,0
evA,0
evB,0

The event evC goes first to highlight an interesting peculiarity. When we run this input into the system, we get:

$ perl ex06_01run.pl tab06_01_01.txt in08_01_01.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evA, c=0.000000
Impossible events, division by -2.22045e-16
 at ex06_01run.pl line 131

If we feed the same events in a different order

evA,0
evB,0
evC,0

we get:

$ perl ex06_01run.pl tab06_01_01.txt in08_01_01a.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evA, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.00000,1.00000,0.66667,0.66667,
hyp2   ,1.00000,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.00000,1.00000,0.66667,0.66667,
hyp2   ,1.00000,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
Impossible events, division by 0
 at ex06_01run.pl line 131

What happened here? Why is there a division by 0, or even by the negative numbers? The negative number comes out of a rounding error. The data in the table has been rounded to only 5 digits after the decimal points to start with, and as the computations run, they introduce more and more rounding errors. Because of that the probability might go slightly out of the normal range, and go very slightly below 0 (as here, looks like it went negative by the lest-significant bit) or above 1. This very tiny negative value should have been 0 if everything had been perfect. If you really want, you could probably compare all the resulting probabilities for being in the proper range and bring them into it if they diverge slightly. Remember though that if your probabilities start to diverge substantially from the proper range, it means that you're getting some formulas terribly wrong, and you better take notice of that. Because of this I wouldn't recommend just blibdly forcing the values into the range, it's easier to let them be, and the small divergences won't matter much while the large divergences will become noticeable and alarming.

So why did it try to divide by 0? Let's run the computation once more with a verbose output:

$ perl ex06_01run.pl -v tab06_01_01.txt in08_01_01.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
P(E)=0.777779
H hyp1 *0.333330/0.222221
H hyp2 *0.000000/0.222221
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evA, c=0.000000
P(E)=1.000000
Impossible events, division by -2.22045e-16

The reason is that P(E)=1, based on the previous seen events the model thinks that this event must absolutely be true, and we're telling it that it's false. The division by 0 means that the events form a combination that is impossible from the standpoint of the model.

But why is it impossible? Let's look again at the training data:

         evA evB evC
4 * hyp1 1   1   1
2 * hyp1 1   0   0
3 * hyp2 0   0   1

No lines in it contain evC=0 and evA=0. The model knows only as much as the training data tells it. We're trying to feed the input (0, 0, 0) which says that none of the symptoms are present, and which wasn't present in the training data.

For this particular case we can add an item to the training data, representing the situation when there is nothing wrong with the subject of the expert system:

            evA evB evC
4    * hyp1 1   1   1
2    * hyp1 1   0   0
3    * hyp2 0   0   1
0.01 * ok   0   0   0

I've added it with a real small count 0.01, so that it would not affect the probabilities of other outcomes very much. How can a count be not a whole number? Remember, it's all relative, from the computation standpoint it's not really counts but weights. Having 9 "real" cases and 0.01 "ok" case is the same as having 900 "real" cases and 1 "ok" case.

This training data translates to the table:

!,,evA,evB,evC
hyp1,0.66593,1,0.66667,0.66667
hyp2,0.33296,0,0,1
ok,0.00111,0,0,0

And let's apply the same input data:

$ perl ex06_01run.pl -v tab08_01.txt in08_01_01.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66593,1.00000,0.66667,0.66667,
hyp2   ,0.33296,0.00000,0.00000,1.00000,
ok     ,0.00111,0.00000,0.00000,0.00000,
--- Applying event evC, c=0.000000
P(E)=0.776916
H hyp1 *0.333330/0.223084
H hyp2 *0.000000/0.223084
H ok *1.000000/0.223084
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.99502,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
ok     ,0.00498,0.00000,0.00000,0.00000,
--- Applying event evA, c=0.000000
P(E)=0.995024
H hyp1 *0.000000/0.004976
H hyp2 *1.000000/0.004976
H ok *1.000000/0.004976
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
ok     ,1.00000,0.00000,0.00000,0.00000,
--- Applying event evB, c=0.000000
P(E)=0.000000
H hyp1 *0.333330/1.000000
H hyp2 *1.000000/1.000000
H ok *1.000000/1.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
ok     ,1.00000,0.00000,0.00000,0.00000,

After the first event evC the hypothesis "ok" inches up a little bit but after the event evA it jumps all the way to 1, becoming the winning one. The former impossibility now became the strong indication for "ok". Note that the hypothesis hyp2 would have surged up in the same way on evA but its probability is already 0 by that time and no amount of multiplication can make it change from 0.

Now let's look at another input set:

evA,0
evB,1
evC,1

It should come as no surprise to you that it also causes a division by 0, as this combination wasn't present in the training data either. In this case it might be reasonable to just stop and say "we don't know, this doesn't match any known outcome". Or it might mean that more than one thing is wrong with the thing we're trying to diagnose, and the symptoms of these illnesses are interacting with each other. It might be still work trying to continue the diagnosis to try figuring out the leading hypotheses. Of course, if there are only two hypotheses available, saying that "both might be true" sounds less useful, but with a large number of hypotheses narrowing down to a small number of them is still useful.

The trick I've come up with for this situation is to cap the confidence in the events. Instead of saying "I'm absolutely sure at confidence 1" let's say "I'm almost absolutely sure at confidence 0.99". And the same thing for 0, replacing it with 0.01. The code I've shown before contains an option (-c) that implements this capping without having to change the input data.

$ perl ex06_01run.pl -v -c 0.01 tab06_01_01.txt in08_01_02.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evA, c=0.010000
P(E)=0.666670
H hyp1 *0.010000/0.336663
H hyp2 *0.990000/0.336663
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.01980,1.00000,0.66667,0.66667,
hyp2   ,0.98020,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.990000
P(E)=0.013202
H hyp1 *0.663337/0.022938
H hyp2 *0.010000/0.022938
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.57267,1.00000,0.66667,0.66667,
hyp2   ,0.42733,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.990000
P(E)=0.809113
H hyp1 *0.663337/0.802931
H hyp2 *0.990000/0.802931
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.47311,1.00000,0.66667,0.66667,
hyp2   ,0.52689,0.00000,0.00000,1.00000,

The first event ~evA drives the evidence strongly towards the hypothesis hyp2 but then the impossible event evB about evens things out, and evC moves the balance a little towards hyp2 again.

The value of the "unsurety cap" will vary the outcomes a little. If you use the smaller cap values, such as 0.00001, it will drive the balance a little sharper towards hyp2, a larger cap, such as 0.1, will produce the more even results. This effect gets more pronounced if you have hundreds of events to look at, and the larger cap values will muddle your results a lot and you won't be able to pick any hypothesis ever. In reality the values like 1e-8 are more suitable for the cap.

The same trick with capping would also work on all events being negative without the hypothesis "ok" but you probably wouldn't want to get a rather random diagnosis when there is nothing wrong with the patient. It's a fine line between being able to handle the cases that vary just a little from the known training data and giving the random diagnoses on the data that doesn't resemble anything known at all.

You might want to detect and count the cases when model goes though the "eye of a needle" of an almost-impossible event, and if it happens more that a certain small number of times, give up and admit that there is no good hypothesis that matches the data. How do you detect them? These are the cases when either the expected P(E) is very close to 0 and the event comes true, or the expected P(E) is very close to 1 and the event comes false. How close is "close"? It depends on your cap value. If you use a small enough cap, something within (cap*10) should be reasonable. Having a small cap value really pays in this case. Note that you still need the capping to make it work. You can't just go without capping, detect the division by 0 and drop the event that caused it. The division by 0 is caused not by one event but by a combination of multiple events and dropping one of them that happened to be the last one isn't right. The capping approach allows you to bring the probabilities back into the reasonable range when you see the last event in an "impossible" sequence.

It also happens that not all the impossible event combinations cause the division by 0. Look at this input:

evA,1
evB,0
evC,1

$ perl ex06_01run.pl -v tab06_01_01.txt in08_01_03.txt
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evA, c=1.000000
P(E)=0.666670
H hyp1 *1.000000/0.666670
H hyp2 *0.000000/0.666670
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.000000
P(E)=0.666670
H hyp1 *0.333330/0.333330
H hyp2 *1.000000/0.333330
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evC, c=1.000000
P(E)=0.666670
H hyp1 *0.666670/0.666670
H hyp2 *1.000000/0.666670
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,

It doesn't match any input in the training data, yet the model never tries to divide by 0 without any capping needed, and it's absolutely sure that the hypothesis hyp1 is correct. Why? As it happens, this input (1, 0, 1) represents a mix of inputs (1, 1, 1) and (1, 0, 0), both of which point to hyp1. When these training data are minced up into the table of probabilities, the information about the correlation of evB and evC that is showing in the cases for hyp1 becomes lost. It just thinks "if evA is true and evB might be true, and evC might be true, then it's hyp1". There are some ways to deal with it, and we'll look at them yet.

Bayes 7: training

Where does the table of probabilities come from? From the training data. What is the training data? It's the information about the previously known cases with the known correct diagnosis and all the known symptoms.

Let's look at what kind of training data (admittedly, pulled from the thin air) I used to create the table in the last example. To remind, the table was:

!,,evA,evB,evC
hyp1,0.66667,1,0.66667,0.66667
hyp2,0.33333,0,0,1

Basically, you start by writing down the known cases. Some of them might have happened in the same way more than once, so instead of writing down the same case multiple times, you just write the count of this kind of cases.

I've made up the following imaginary cases:

----     evA evB evC
4 * hyp1 1   1   1
2 * hyp1 1   0   0
3 * hyp2 0   0   1

9 cases total, with hypothesis hyp1 happeing with two different combinations of symptoms (4 and 2 cases), and the hypothesis hyp2 happening with one combination of symptoms (3 cases).

The prior probabilities P(H) of hypotheses get computed by dividing the number of cases with this hypothesis by the total number of cases:

P(hyp1) = (4+2)/(4+2+3) = 6/9 = 0.66667
P(hyp2) = 3/(4+2+3) = 3/9 = 0.33333

The conditional probabilities of the events P(E/H) are computed by dividing the number of cases where this event is true for this hypothesis by the total number of the cases for this hypothesis:

P(evA|hyp1) = (4+2)/(4+2) = 6/6 = 1
P(evB|hyp1) = 4/(4+2) = 4/6 = 0.66667
P(evC|hyp1) = 4/(4+2) = 4/6 = 0.66667

P(evA|hyp2) = 0/3 = 0
P(evB|hyp2) = 0/3 = 0
P(evC|hyp2) = 3/3 = 1

These numbers go into the table. Very simple. Of course, with the real-sized data sets of many thousands of cases you'd have to do a lot more calculations but it's fundamentally straightforward. At least for this basic approach. There are more complicated approaches, we'll get to them later.

Bayes 6: basic code example

After going through in words, I want to show you the code, with a simple example. The code is fairly independent of the actual knowledge that is contained in the probability tables, so I'll make the probability tables loadable.

The probability tables format will be like this:

# this is a simple table example
!,,evA,evB,evC
hyp1,0.66667,1,0.66667,0.66667
hyp2,0.33333,0,0,1

The comments start with a "#".

The data lines contain the comma-separated values (CSV). The lines, except the one that starts with a "!" describe the hypotheses. The first field is the hypothesis name, the second field is the prior probability P(H), and the rest of the fields are the conditional probabilities P(E|H) for every event and this hypothesis. The lines may contain a trailing comma, it will get ignored:

hyp_name,P(H),P(E1/H),P(E2/H),...

The line that starts with "!" contains the event names. If this line is absent, the events will be assigned the numeric names, starting from 1. The first two fields of the line are ignored, and the rest contain the names of the events, forming the table header, where each table column correcponds to an event and its conditional probabilities.

The same table can also be formatted with the nice alignment:

!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,

When my code runs, it will print its information about the tables in this aligned format. I've arbitrarily limited the width of the columns at 7, so any longer names will be truncated in the print-out. Good enough for the examples. The spaces around the hypothesis and event names will be ignored.

The other data item will be the input data describing a particular case. It consists of the CSV lines in format

event_name,confidence

The event name must match one of the names defined in the table, and the confidence is in the range [0..1], with 0 meaning that the event is definitely false, and 1 meaning that it's definitely true.

Now to the code. Eventually I'll make the code downloadable as files but for now I'll probably do the changes in it for some time, so I'l just copy-paste the current version here to reduce the confusion between the code and text. I'll talk more about how and why it works but first just the pasted code and some examples of its use.

#!/usr/bin/perl
#
# Running of a Bayes expert system.

use strict;
use Carp;

our @evname; # the event names
our %evhash; # hash of event names to indexes
our @hypname; # the hypothesis names
our @phyp; # the hypothesis probabilities
our @peh; # probabilities P(E/H) by hypothesis by event

# options
our $verbose = 0; # verbose printing during the computation
our $cap = 0; # cap on the confidence, factor accounting for fuzziness of training data,
  # limits the confidence to the range [$cap..1-$cap]

sub load_table($) # (filename)
{
  my $filename = shift;

  @evname = ();
  %evhash = ();
  @hypname = ();
  @phyp = ();
  @peh = ();

  my $nev = undef; # number of events minus 1

  confess "Failed to open '$filename': $!\n"
    unless open(INPUT, "<", $filename);
  while(<INPUT>) {
    chomp;
    s/,\s*$//; # remove the trailing comma if any
    if (/^\#/ || /^\s*$/) {
      # a comment line
    } elsif (/^\!/) {
      # row with event names
      @evname = split(/,/); # CSV format, the first 2 elements gets skipped
      shift @evname;
      shift @evname;
    } else {
      my @s = split(/,/); # CSV format for hypothesis
      push @hypname, shift @s;
      push @phyp, shift @s;

      if (defined $nev) {
        if ($nev != $#s) {
          close(INPUT);
          my $msg = sprintf("Wrong number of events, expected %d, got %d in line: %s\n",
            $nev+1, $#s+1, $_);
          confess $msg;
        }
      } else {
        $nev = $#s;
      }
      # the rest of fields are probabilities, just take the whole array by reference
      push @peh, \@s;
    }
  }
  close(INPUT);

  if ($#evname >= 0) {
    if ($#evname != $nev) {
      my $msg = sprintf("Wrong number of event names, %d events in the table, %d names\n",
        $nev+1, $#evname+1);
      confess $msg;
    }
  } else {
    for (my $i = 0; $i <= $nev; $i++) {
      push @evname, ($i+1)."";
    }
  }

  for (my $i = 0; $i <= $#evname; $i++) {
    $evname[$i] =~ s/^\s+//;
    $evname[$i] =~ s/\s+$//;
    $evhash{$evname[$i]} = $i;
  }
  for (my $i = 0; $i <= $#hypname; $i++) {
    $hypname[$i] =~ s/^\s+//;
    $hypname[$i] =~ s/\s+$//;
  }
}

sub print_table()
{
  my $w = 7; # width of every field
  my $tfmt = "%-${w}.${w}s,"; # format of one text field
  my $nw = $w-2; # width after the dot in the numeric fields field
  my $nfmt = "%-${w}.${nw}f,"; # format of one numeric field (does the better rounding)
  
  # the title row
  printf($tfmt, "!");
  printf($tfmt, "");
  foreach my $e (@evname) {
    printf($tfmt, $e);
  }
  print("\n");
  # the hypotheses
  for (my $i = 0; $i <= $#hypname; $i++) {
    printf($tfmt, $hypname[$i]);
    printf($nfmt, $phyp[$i]);
    foreach my $e (@{$peh[$i]}) {
      printf($nfmt, $e);
    }
    print("\n");
  }
}

# Apply one event
# evi - event index in the array
# conf - event confidence [0..1]
sub apply_event($$) # (evi, conf)
{
  my ($evi, $conf) = @_;

  # compute P(E)
  my $pe = 0;
  for (my $i = 0; $i <= $#hypname; $i++) {
    $pe += $phyp[$i] * $peh[$i][$evi];
  }

  if ($verbose) {
    printf("P(E)=%f\n", $pe);
  }
  
  #compute the divisor
  my $div = $pe*$conf + (1. - $pe)*(1. - $conf);
  if ($div < 1e-10) {
    confess (sprintf("Impossible events, division by %g\n", $div));
  }

  # update the hypothesis probabilities
  # (it works more efficiently if the array $peh is indexed by the event first
  # but for the simple example who cares)
  for (my $i = 0; $i <= $#hypname; $i++) {
    my $mul = $peh[$i][$evi]*$conf + (1. - $peh[$i][$evi])*(1. - $conf);
    if ($verbose) {
      printf("H %s *%f/%f\n", $hypname[$i], $mul, $div);
    }
    $phyp[$i] *= $mul/$div;
  }
}

# Apply an input file
sub apply_input($) # (filename)
{
  my $filename = shift;

  confess "Failed to open the input '$filename': $!\n"
    unless open(INPUT, "<", $filename);
  while(<INPUT>) {
    chomp;
    next if (/^\#/ || /^\s*$/); # a comment

    my @s = split(/,/);
    $s[0] =~ s/^\s+//;
    $s[0] =~ s/\s+$//;

    confess ("Unknown event name '" . $s[0] . "' in the input\n")
      unless exists $evhash{$s[0]};
    my $evi = $evhash{$s[0]};

    my $conf = $s[1] + 0.;
    if ($conf < $cap) {
      $conf = $cap;
    } elsif ($conf > 1.-$cap) {
      $conf = 1. - $cap;
    }
    printf("--- Applying event %s, c=%f\n", $s[0], $conf);
    &apply_event($evi, $conf);
    &print_table;
  }
  close(INPUT);
}

# main()
while ($ARGV[0] =~ /^-(.*)/) {
  if ($1 eq "v") {
    $verbose = 1;
  } elsif ($1 eq "c") {
    shift @ARGV;
    $cap = $ARGV[0]+0.;
  } else {
    confess "Unknown switch -$1";
  }
  shift @ARGV;
}
&load_table($ARGV[0]);
&print_table;
if ($#ARGV >= 1) {
  &apply_input($ARGV[1]);
}

How to run it: Suppose the table above is locate din the file tab06_01_01.txt, and the code is in ex06_01run.pl. The first the code can do is just load the table and print it back with the nice formatting:

$ perl ex06_01run.pl tab06_01_01.txt
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,

Next, it can apply the input data for computation. Suppose the input data is in in06_01_01_01.txt:

evA,1
evB,0
evC,0

Let's apply it:

$ perl ex06_01run.pl tab06_01_01.txt in06_01_01_01.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evA, c=1.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,

It prints the table with the new probabilities after applying every event, so that we can trace, how the probabilities change. The only thing that changes is really the second column, P(H) of the hypotheses. The conditional probabilities of the events don't change.

In this case the very first event decided the fate: it drove the probability all the way for hyp1, and the rest of the events just stayed consistent with that. If we reorder the events, the intermediate computations will change:

evB,0
evA,1
evC,0

$ perl ex06_01run.pl tab06_01_01.txt in06_01_01_04.txt
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.40000,1.00000,0.66667,0.66667,
hyp2   ,0.60000,0.00000,0.00000,1.00000,
--- Applying event evA, c=1.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,

If evB is applied first, it drives the probabilities toward hyp2, but then the next event brings it back to the full glory of hyp1. Which basically means that if evA is true, there is no way hyp2 could happen, so hyp1 wins big time.

The order in which the events are applied doesn't matter. Either way it will come to the same result. In reality it might not be the best result since the table in this form doesn't represent the cross-correlation between the events quite right. But for a given table the result is the same independently of the order in which the events are applied.

The Bayes formula is applied for one event at a time in the function apply_event(). The formula is exactly as described before.

P(H|C(E)) = P(H) * ( C(E)*(0+P(E|H)) + (1-C(E))*(1-P(E|H)) )
    / ( C(E)*(0+P(E)) + (1-C(E))*(1-P(E)) )

The only catch is that P(E) is not present anywhere in the table. It has to be calculated from the other available data. This is because the probabilities of the events change after the knowledge about the previous events has been taken into the account. It's the biggest effect of the cross-correlation of the events and it absolutely has to be taken into account or you'll be seeing the probability values around 15 or so. There are other subtler effects as well but more about that later.

The event probability is computed from adding up its weighed conditional probabilities:

P(E) = sum(P(H) * P(E/H))

Since the P(H) values change after every events and P(E/H) stay the same, the values of P(E) also change.

If used with the option -v, the program will print the verbose intermediate data as it preforms the computations. With the original event order it will print:

$ perl ex06_01run.pl -v tab06_01_01.txt in06_01_01_01.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evA, c=1.000000
P(E)=0.666670
H hyp1 *1.000000/0.666670
H hyp2 *0.000000/0.666670
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.000000
P(E)=0.666670
H hyp1 *0.333330/0.333330
H hyp2 *1.000000/0.333330
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
P(E)=0.666670
H hyp1 *0.333330/0.333330
H hyp2 *0.000000/0.333330
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,

With the different order:

$ perl ex06_01run.pl -v tab06_01_01.txt in06_01_01_04.txt 
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.66667,1.00000,0.66667,0.66667,
hyp2   ,0.33333,0.00000,0.00000,1.00000,
--- Applying event evB, c=0.000000
P(E)=0.444449
H hyp1 *0.333330/0.555551
H hyp2 *1.000000/0.555551
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,0.40000,1.00000,0.66667,0.66667,
hyp2   ,0.60000,0.00000,0.00000,1.00000,
--- Applying event evA, c=1.000000
P(E)=0.400001
H hyp1 *1.000000/0.400001
H hyp2 *0.000000/0.400001
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,
--- Applying event evC, c=0.000000
P(E)=0.666670
H hyp1 *0.333330/0.333330
H hyp2 *0.000000/0.333330
!      ,       ,evA    ,evB    ,evC    ,
hyp1   ,1.00000,1.00000,0.66667,0.66667,
hyp2   ,0.00000,0.00000,0.00000,1.00000,

As you can see, the probabilities for evA and evB were wildly different in the first two rounds. But either way they were applied, by the time the evC gets into the fray, it has the same probability in both cases.

Bayes 5: the big picture

It's not enough to just have the magic formulas, an expert system has to actually use them. So, how do you use them?

You start with the table of probabilities built from the training data. We'll look at the build later, for now let's just assume that the table is there.

Each event represents some symptom, so you look at the symptoms one by one, check if they are present, making the events true or false (or some confidence value in between), and adjust the probabilities of the hypotheses. At the end we pick the hypothesis that is most likely. It's usually a good idea to pick some absolute probability requirement as well. Basically, if you have 98 hypotheses at 1% probability and one hypothesis at 2% probability, it probably doesn't mean that the 2% hypothesis is the right one. It probably means that the system has no clue about the right answer, and the answer it gives should be "I don't know". You usually want the winning hypothesis to be at least somewhere in the vicinity of 70% probability. The optimal value might take some fine-tuning of your system. There is another important consideration to this point: are the hypotheses mutually exclusive or could multiple of them be true at the same time? I'll discuss this consideration in detail later, and for mow just leave it at that.

There are two ways to look for the symptoms. In one way you have some automated way to look for all the symptoms you know: in medicine that might be some panel of lab tests, in car diagnostics it might be some 107-point inspection, and of course as far as the computers or software are concerned, it's just the perfect application for the automated tests. The uses like the spam detection are just the perfect example: all the information you can get is right there in the e-mail, just go and extract it. If you have the results for all the symptoms, you just go and apply them all, one after another, and at the end you have the probabilities for all the hypotheses and can pick the best one.

In the other way the examination of a symptom requires some manual procedure, such as having a doctor poke the patient in some way, or ordering some test. There is some cost associated with performing the test. In some cases it might even be an an attempt at fixing the issue, mixing the symptoms and hypotheses a bit: such as "try this antibiotic" or "replace the alternator, did the problem go away"? When the symptoms are like this, there are two consequences: first, to save on the cost of diagnostics, you want to stop after examining as few symptoms as possible, and second, the order of examining the symptoms matters.

I've dealt in reality only with the situation when all the symptoms are available up-front, so I'll re-tell the solutions for the expensive symptoms I've read in that old British book (perhaps there are some better ones invented by now?).

A certain way to stop early is in the case if one hypothesis becomes certainly more probable than all the others. What does that mean, "certainly"? It means that if we look at all the symptoms we haven't examined yet and suppose that they go the worst way for this hypothesis, it would still beat all the other hypotheses, satisfying the requirements to accept it as the best one.

Another certain way to stop is if no hypothesis, even if all the yet unexamined symptoms go the best way for them, could satisfy the acceptance requirements. It's just the early admission of impending defeat without the useless thrashing beforehand.

The way you do it is by calculating the best and worst cases for each hypothesis after applying each symptom. You can also optimize it by remembering which hypotheses start failing the cut unconditionally, and skipping them in all the future checks. Obviously, if you have all the events available up-front, there is no point in stopping early, it's faster to just apply all the events rather than re-calculating the best and worst cases after every event.

How exacly do you compute the best and the worst cases? Like this:

• Start with the current probability of the hypothesis, and assign it as teh prior values Pbest(H) and Pworst(H)
• For all the un-examined symptoms:
• Check if P(E|H)/P(E) > 1 (i.e. the event has a positive correlation with this hypothesis)
• if this check is true, compute: Pbest(H) *= P(E|H)/P(E); Pworst(H) *= P(~E|H)/P(~E)
• else, compute Pbest(H): *= P(~E|H)/P(~E); Pworst(H) *= P(E|H)/P(E)

Now, what is the best next symptom to try out? The one that will have the strongest effect by affecting the probabilities of the hypotheses in the strongest way. I don't remember how that book defined the strength but here are some logical considerations, and it would probably work out in the same way as they did in that book.

We can use the difference between the probabilities as the measure of the effect of the event on a hypothesis:

abs(P(H|E) - P(H|~E))

Another possible way would be to divide the probabilities by each other but that is fraught with troubles. Let me show you why. Suppose we
define the effect as

max(P(H|E), P(H|~E)) / min(P(H|E), P(H|~E))

Max() and min() are needed to keep the value always above 1. But the minimum may always become 0, which would result in the infinities, which is not good. If we try to divide the other way around, min by max, we will produce the same results for 0/0.0001 and 0/0.9999, which is obviously also a bad metric. Thus the subtraction it is.

And then the obvious way to get the effect of an event on the whole set of hypotheses is to add up its effects on each individual hypotheses. Whichever event gets the highest metric is the most important event.

But there is also the question of the cost of checking the symptom. As everyone who performs the diagnostics knows, the reasonable approach is to try the cheapest things first. Which means that we should probably multiply that cumulative effect of an event by the cost of checking it, and use that as the metric.

The cost approach can be taken further: each hypothesis means a diagnosis that implies a prescription for a certain treatment. That treatment also has some cost. And again the reasonable thing is to hope for the diagnoses with the cheap treatments. So when computing the effect of an event on the whole set of hypotheses we can give the preference to the hypotheses with a cheaper treatment, for example by multiplying the effect on this hypothesis by 1/cost_of_treatment or log(1/cost_of_treatment) or some such.

What if a single test provides the knowledge about multiple symptoms? Then obviously we need to look for the choice not of the next symptom but of the whole next test. The effect of a test can be computed by adding up the effects of all the symptoms it resolves. Well, not necessarily exactly adding, this has issues with the correlated symptoms, but it's a reasonable first approach, and I'll talk more about this issue later.

Another interesting variation (that's not from the book, that's my own musings) is if the hypothesis are used to mean not just a diagnosis but a confirmed diagnosis, meaning "we digagnosed, treated, and confirmed that the treatment had cured the problem". In this case not just the symptoms as such but also the prescriptions for the treatment become the event (with the event being true if the treatment had worked). This nicely folds the knowledge about the alternative treatments into the same model. Obviously, this approach needs to have some stricter criteria of which treatments are worthy enough to try, the same as for the hypotheses in the more basic approach, to prevent the system from thrashing in vain when it's unable to make any useful conclusion.

This approach can be merged with the automated tests as well: We start by applying the results of all the automated tests as events, and then choose one of the possible treatments as the "next event to test" using the principles described above. If it fixes the issue, great, our job is done. Otherwise we can try the next treatment, and maybe re-run the automated tests, see if anything has changed, and take the changes into account before recommending the next treatment.

Bayes 4: the fuzzy logic

Before moving on to the other aspects, I want to wrap up with the basic formula, this time in its fuzzy-logic version. As I've already mentioned, sometimes it's hard to tell if the event is true or not. It would be nice to have some way to express the idea of the event being very likely or slighlty likely or slightly or very unlikely. It would be even nicer to represent our confidence in the event happening as a continuous range, with 1 meaning "the event is definitely true", 0 meaning "the event is definitely false", 0.5 meaning "we can't tell at all if the event is true or false", and the points in between representing the values in between. Let's use the letter C to denote this confidence, or maybe C(E) to show the confidence of which event we're talking about.

There is a way to generalize the Bayes formula for this confidence range. Way back when I've read it in that British book without explantions and marveled at how people come up with these formulas but in the more recent times I was able to re-create it from scratch from the logical consideration, so I'll show you, where does this marvelous formula come from (and a little later I'll show you where the bayes formula comes from too).

First, let's look at the boundary conditions. Kind of obviously, with C(E)=1 the formula must give the same result as the basic Bayes formula with E being true, with C(E)=0 it must give the same result as the basic formula for E being false, and for C(E)=0.5 it must leave the P(H) unchanged.

Let's look at the basic formulas again:

P(H|E) = P(H) * ( P(E|H) / P(E) )
P(H|~E) = P(H) * ( P(~E|H) / P(~E) )

If we substitute

P(~E|H) = 1 - P(E|H)
P(~E) = 1 - P(E)

then the formulas become:

P(H|E) = P(H) * ( P(E|H) / P(E) )
P(H|~E) = P(H) * ( (1-P(E|H)) / (1-P(E)) )

In both cases we multiply the same P(H) by some coefficient. The coefficient in both cases is computed by division of something done with P(E|H) by something done with P(E). For symmetry we can write in the first formula P(E|H) as 0+P(E|H) and P(E) as 0+P(E), then the formulas become even more like each other:

P(H|E)  = P(H) * ( (0+P(E|H)) / (0+P(E)) )
P(H|~E) = P(H) * ( (1-P(E|H)) / (1-P(E)) )

If we represent the points 0+P(E|H) and 1-P(E|H) graphically as a line section, we can see that in both cases they are located at the distance P(E|H), only one is to the right from 0, and another one is to the left of 1:

     0   P(E|H)         1-P(E|H)  1
.....|......|..............|......|.......
     ------>\--------------/<------

And the same works for P(E). This way it looks even more symmetrical.

The natural follow-up from here is that we can use the confidence value C(E) to split the sections [ 0+P(E|H), 1-P(E|H) ] and [ 0+P(E), 1-P(E) ] proportionally. For example, if C(E) = 0.75, it will look like this:

           P(0.75 E|H)
                |
     0   P(E|H) |       1-P(E|H)  1
.....|......|...|..........|......|.......
             0.25   0.75

The (0+P(E|H)) gets taken with the weight C(E) and (1-P(E|H)) gets taken with the weight 1-C(E) and they get added together. If C(E) is 0.75, it means that the split point will be located closer to (0+P(E|H)), dividing the original range at its quarter-point.

In the formula form this can be expressed as:

P(H|C(E)) = P(H) * ( C(E)*(0+P(E|H)) + (1-C(E))*(1-P(E|H)) )
 / ( C(E)*(0+P(E)) + (1-C(E))*(1-P(E)) )

Let's check that this formula conforms to the boundary conditions formulated above:

P(H|C(E)=1) = P(H) * ( 1*(0+P(E|H)) + (1-1)*(1-P(E|H)) )
 / ( 1*(0+P(E)) + (1-1)*(1-P(E)) )
 = P(H) * P(E|H) / P(E)

P(H|C(E)=0) = P(H) * ( 0*(0+P(E|H)) + (1-0)*(1-P(E|H)) )
 / ( 0*(0+P(E)) + (1-0)*(1-P(E)) )
 = P(H) * (1-P(E|H)) / (1-P(E))

P(H|C(E)=0.5) = P(H) * ( 0.5*(0+P(E|H)) + (1-0.5)*(1-P(E|H)) )
 / ( 0.5*(0+P(E)) + (1-0.5)*(1-P(E)) )
 = P(H) * ( 0.5*P(E|H)) + 0.5*(1-P(E|H)) )
 / ( 0.5*P(E) + 0.5*(1-P(E)) )
 = P(H) * ( 0.5*P(E|H)) + 0.5 - 0.5*P(E|H)) )
 / ( 0.5*P(E) + 0.5 - 0.5*P(E)) )
 = P(H) * 0.5 / 0.5
 = P(H)

They all match. So this is the right formula even though a bit complicated. We can make it a bit simpler by splitting it into two formulas, one representing the function for the proportional splitting, and another one computing the probability with the use of the first function:

Pc(P, C) = P*C + (1-P)*(1-C)
P(H|C(E)) = P(H) * Pc(P(E|H), C(E))/Pc(P(E), C(E))

Now you might complain "but I don't like the shape of this dependency I feel that the confidence should be represented in a more non-linear way in one way or another". Well, if it's too linear for you, you don't have to use it directly. You can define your own confidence function Cmy(E) and map its values into C(E) in any kind of non-linear form, and then use the resulting C(E) values for computation.

In fact, I like a variation of that trick myself. I like limiting the values of C(E) to something like 0.00001 and 0.99999 instead of 0 and 1. It helps the model recover from the situations that it would otherwise think impossible and provide a way to handle the problem of the overfitting. But more on that later.

Bayes 3: the magic formula

The next thing is the conditional probability. It's written as P(H|E) or P(E|H).

P(H|E) means the probability of the hypothesis H being true in case if we know that the event E is true.

P(E|H) means the probability of the event E being true in case if the hypothesis H turns out to be true.

And now finally the magic formula discovered by Bayes:

P(H|E) = P(H) * P(E|H) / P(E)

It tells us how the probability of a hypothesis gets adjusted when we learn some related data. I haven't watched "House, MD" closely enough to give an example from the medical area, so I'll give one about the car repair.

Suppose our hypothesis is "the battery is discharged". One way to test it is to turn on the headlights (with the engine not running) and see that they come on at full brightness, that would be our event.

Suppose our original ("prior") probability P(H) was 0.1. And suppose the probability of the headlights being able to come up was 0.7. You might think that no way would the headlights shine brightly with the discharged battery, but the headlights consume a lot less current than the starter, and the battery might be capable of providing the small amount of current without too much of a voltage drop, so the headlights might still look decently bright to the eye. So let's say that P(E|H), the probability of the headlights shining brightly with a discharged battery, is 0.05.

Where do all these numbers come from? Rights now I've pulled them out of my imagination. In a real expert system they would come from two sources: the training data and the computation of the previous events. Some people might say that the numbers might come from the expert opinions but spit them in the eye, getting the good probabilities without a good set of training data is pretty much impossible. We'll get to the details of the training data a little later, now let's continue with the example.

So we've tested the headlights, found them shining bright, and now can plug the numbers into the formula:

P(H|E) = 0.1 * 0.05 / 0.7 = 0.00715

As you can see, the probability of this hypothesis took a steep plunge. Why did it happen? Basically, the effect depends on the difference between how probable this effect is in general (P(E)) and how probable it is if this hypothesis is true (P(E|H)). The formula can be reformatted to emphasise this:

P(H|E) = P(H) * ( P(E|H) / P(E) )

If P(E|H) is greater than P(E) then the probability of the hypothesis gets pulled up. If it's less then the probability gets pulled down.

What if we try the experiment and the result is false? Obviously, we can't just use this formula, since E is not true. But then the event ~E will be true, and we should instead use the similar formula for ~E and also adjust the probability of the hypothesis:

P(H|~E) = P(H) * P(~E|H) / P(~E)

Where do P(~E|H) and P(~E) come from? Again, from the training data, and they are connected with the direct values in the obvious way:

P(~E|H) = 1 - P(E|H)
P(~E) = 1 - P(E)

A lot of examples of Bayesian logic I've found on the Internet don't adjust the probabilities if an event is found to be false, and this is very, very wrong. Basically, they treat the situation "the event is known to be false" as "the event is unknown". If you have a spam filter and know that a message doesn't contain the word "viagra", you shouldn't treat it the same way as if you don't know if this word is present. The result will be that you'll be driving the probabilities of the message being spam way high on seeing any of the suspicious words. In case if you wonder if ignoring the negative event is sheer stupidity, the answer is no, it isn't, there are reasons to why people do this, they are trying to compensate for the other deficiencies of the bayesian method. Only a rare spam message indeed would contain the whole set of suspicious words, so if you properly test each word and apply it as positive or negative, you end up with pretty low probability values. I'll talk more about it later, although I don't have the perfect solutions either.

Continuing the example, suppose the headlights didn't come on brightly, and we have the probabilities

P(~E|H) = 1 - 0.05 = 0.95
P(~E) = 1 - 0.7 = 0.3

Then we can substitute them into the formula:

P(H|~E) = 0.1 * 0.95 / 0.3 = 0.31667

The hypothesis of a dead battery took quite a boost!

Now let me show you something else. Suppose the prior probability P(H) was not 0.1 but 0.5. We put the data with this single adjustment into the same formula and get:

P(H|~E) = 0.5 * 0.95 / 0.3 = 1.58333

Shouldn't the probability be in the range between 0 and 1? It sure should. If you get a probability value outside of this range, it means that something is wrong with your input data. Either you've used some expert estimations of probabilities and these experts turned out to be not any good (that's why using a proper set of training data is important), or you forgot to adjust the probabilities for the previous computations you've done. Or maybe the rounding errors from the previous computations have accumulated to become greater than the actual values.

Bayes 2: hypotheses and events

A bayesian expert system deals with the probabilities of hypotheses. It tries to compute these probabilities based on the experimentation with the black box, and choose the most probable one (or maybe more than one).

The probability of a hypothesis H is commonly written as P(H). Naturally, it's a real number in the range between 0 ("definitely false") and 1 ("definitely true"). If you haven't dealt with the probabilities before, now you know this.

Each hypothesis H has an opposite hypothesis that I'll write as ~H (in the mathematical texts it's written as H with a line above it but it's inconvenient to type, so I chose a notation like the C language). If the hypothesis Hlupus is "the patient has lupus" then the opposite ~Hlupus will be "the patient doesn't have lupus". I'll also call ~H the "negative hypothesis", as opposed to the "positive" H.

The probabilities of positive and negative hypotheses are connected:

P(H) + P(~H) = 1

Which means that whatever happens, exactly one of H and ~H is true. We might not know for sure which one is true but if we think that H has the probability 0.7 (or in other words 70%) of being true, this means that it has the probability 0.3 of being false, while the opposite ~H must have the probability 0.3 of being true and 0.7 of being false. This is a very rigid connection, and we can use it for various deductions.

An event is not just an experiment but an experiment with a particular outcome. Not just "look at the patient's throat" but "I've looked at the patient's throat and found that it's sore". Just as the hypotheses, the events have their opposites too. If the event E is "I've looked at the patient's throat and found that it's sore", its opposite ~E will be "I've looked at the patient's throat and found that it's NOT sore". The experiment is the same but its outcome is opposite for the negative event. It's real important to note that the case "I haven't looked at the patient's throat yet" is different: until we look at the throat, neither E or ~E have happened. But once we look at the throat, we can tell that either E is true and ~E is false for this particular patient or the other way around, E is false and ~E is true. A typical mistake (one found in Wikipedia for example) is to mix up the cases when we haven't looked at the throat and when we looked at the throat and found that it's not sore. If you mix them up, your calculations will all be way wrong. Beware of it.

A clever reader might say "what if I looked at the patient's throat and can't say for sure whether it's sore or not, it's slightly reddish but not totally so"? Glad you asked, that's the department of the fuzzy logic, and I'll get to it later. But for now let's assume that we can tell whether it's one or the other. If we can't tell then we'll just treat it the same way as we haven't done the experiment at all (which is by the way consistent with how it's treated by the fuzzy logic).

We can't know the outcome of an experiment before we do the experiment, whether E or ~E will be true. But we can estimate the probability of E in advance, based on the other information we have collected so far. This advance estimation is called the probability of an event, P(E). For example, if we know that 30% of patients visiting a doctor have a sore throat, we can estimate P(Esorethroat)=0.3 based on our knowledge that this person came to see a doctor.

Just as with hypotheses, the probabilities of the complementary events are connected:

P(E) + P(~E) = 1

And it's really not any different from a hypothesis. It really is a hypothesis about the outcome of an experiment we haven't done yet, which then collapses to the true or false after the experiment is done and the result is observed.

Some experiments may have more than two possible outcomes. For example, suppose our black box under investigation contains some colored balls in it, and we can stick a hand in it and pull out a ball that might be of one of three colors: blue, yellow or green. This can be treated as three separate events:

Eblue: the ball is blue
Eyellow: the ball is yellow
Egreen: the ball is green

The total of all three probabilities will still be 1:

P(Eblue) + P(Eyellow) + P(Egreen) = 1

The connection between the probabilities of events means that these events are not independent. Observing any of these three events changes the probabilities of the other two events.

But the sum of probabilities in each pair will also be 1:

P(Eblue) + P(~Eblue) = 1
P(Eyellow) + P(~Eyellow) = 1
P(Egreen) + P(~Egreen) = 1

The connected events can be processed sequentially. We can predict the probabilities of the events that haven't been processed yet after each step of processing. We can achieve this in the real world by having an assistant pull out the ball out of the box and then asking him questions with the answers of "yes" or "no". Suppose we know for some reason that there is an equal probability for each color. Thus initially:

P(Eblue) = 1/3
P(~Eblue) = 2/3

P(Eyellow) = 1/3
P(~Eyellow) = 2/3

P(Egreen) = 1/3
P(~Egreen) = 2/3

Then we have the assistant pull out the ball and ask him: "Is the ball blue?" If the answer is "No" then the event probabilities change:

P(Eblue) = 0
P(~Eblue) = 1

P(Eyellow) = 1/2
P(~Eyellow) = 1/2

P(Egreen) = 1/2
P(~Egreen) = 1/2

The probabilities of Eblue and ~Eblue have collapsed to 0 and 1 but the remaining two outcomes are now equi-probable between two, not three.

Now we can ask the next question: "Is the ball yellow?". If the answer is again "No" then the probabilities change like this:

P(Eblue) = 0
P(~Eblue) = 1

P(Eyellow) = 0
P(~Eyellow) = 1

P(Egreen) = 1
P(~Egreen) = 0

At this point we also become quite sure that the ball must be green, and we can as well skip asking about it. Not in the real world though. In the real world something might have changed since we've made our original estimations and more colors became available to the box-builder. Or maybe the fourth outcome is extremely rare and our training data didn't contain any example of it. Thus it could happen that we ask "Is the ball green?" expecting to hear "Yes" and get a "No". Obviously, this would be a curve-ball that disturbs the model, but I'll show a possible way of dealing with it.

Bayes 1: basic terms

A couple of years ago, while at Google, I've done a "20% project" of a bayesian expert system. The expert system existed before me but I was able to improve it. While doing so, I've figured out a few things about the Bayes logic. Maybe they're commonplace knowledge but I haven't seen them on the Internet, and moreover, most of the descriptions I've seen on the Internet tend to get things somewhat wrong. Especially in Wikipedia. Obviously, I can't talk about the expert system itself, and after a couple of years I don't remember that much about it anyway. But I took notes about the general formulas I've figured out. I've recently found those notes, and want to share them.

This has really nothing to do with CEP, though I guess the CEP can be used to drive the expert systems. I just want to share this knowledge somewhere, and why not here?

Some 25 years or so ago I've read a British collection of articles on the artificial intelligence. One of these articles was on the expert systems based on the bayesian formulas, which quite impressed me. Then I've played a little more with this stuff in college. I vaguely remembered that the author of that article had developed the system called Microexpert but now what I can find about Microexpert doesn't seem to quite match it, so maybe that one was from another article.

When I started working on that recent expert system and went to refresh my knowledge, I've thought, by now Bayes is used all over the place, and the descriptions should be all over the Internet, right? Nope. Not only there isn't much, but also a lot of them apply the Bayes formula somewhat wrong. So I've had to restore this stuff from vague memories and logical reasoning.

First I'll describe what I've restored from my memory and then I'll go into what I've come up with while actually mucking with the realities of the expert systems. This might be well known and published in some obscure articles or maybe even in well-known books, but I don't know where (if you have pointers, you're welcome to comment). But it's not on the Internet and not in a popular form. I want to put it there.

Let's start with the basics. The basic description of the formula created by Bayes actually can be found all over the Internet. You can look at Wikipedia (in a crappy article), or say read a lengthy explanation at http://yudkowsky.net/rational/bayes or http://lesswrong.com/lw/2b0/bayes_theorem_illustrated_my_way/ . But I'll try to give my own short descriptions and then go further. Read on.

The probability theory talks about hypotheses and events. Both of them relate to some system being examined, some kind of a "black box" whose contents we're trying to guess based on observations. If we could just look inside the box, there would be no need for all this, but in reality we could only poke the box in certain ways and then have to deduce its contents to the best of our knowledge based on that limited experimentation. The contents of the box is usually not very random. We normally know a good deal about what could be in the black boxes we get, we just don't know what is in this particular black box. The goal of an expert system is to look at a black box, do some poking, and classify the contents of the black box to one of the known categories with a good certainty. The expert system example in that old British book was trying to diagnose the automobile breakages. A medical expert system would be trying to diagnose what sickness a patient has.

It's kind of intuitive that the hypotheses are the descriptions of what might be in the box, and thus we're trying to select which contents seems most likely based on out experimentations, what hypothesis is most probable.

The events are the experiments. We might be able to select the experiments to perform or might be limited to observing what the mechanism in the black box does on its own.

The hypotheses are usually denoted H or Hsomething. The events are usually denoted E or Esomething.

something old

Hey, it turns out that the patent for the dynamic schema modification that I worked on back in the Aleri days got actually issued, way back in 2009: http://www.google.com/patents/US20090037769

I've totally missed it and thought that it was still mired in the process. Looks like it took only 8 months from the final filing to the issue.

On simplicity

Today someone brought to my attention the quote from http://blog.rongarret.info/2009/02/css-and-meaning-of-life.html :

<<And then there is one element of my quality metric which seems to be at  the heart of the controversy: I believe that computers are meant to  serve people and not the other way around.  That means that if something
 inherently simple is difficult to do, it's wrong.

Not  everyone agrees with this.  A surprising number of people like things to be difficult.  I have never understood this, but it is a fact.  Some  people revel in the intellectual studliness of mastering, say, C++
template programming or x86 assembler.  I don't.  I think C++ is a pain in the ass.  Life is short, and I would rather spend my allotted time skiing or making love to my wife than worrying about whether or not I need to define a virtual destructor in my base class in order to avoid memory leaks.  And that means that when I want, say, for something to appear centered on the screen, I want to be able to do it by typing "center" rather than "margin:auto".  If I can't tell a computer to center something by saying "center", then as far as I'm concerned the computer is wrong, not me.

Like I said, not everyone buys in to this worldview.  I don't get it, but it is surprisingly common.  I think it might have something to do with job security.  I first encountered the complicated-is-beautiful mindset when I was working for NASA.  I had an idea for how to simplify spacecraft sequencing and make it less error prone.  I figured everyone would love it.  To the contrary, the establishment fought the idea tooth-and-nail (and still does to this day).  And this attitude is everywhere.  It's the reason the tax code is so fucked up.  It's the reason the law is so byzantine.  It's the reason that languages like C++ and Perl thrive.>>

I see a major issue with this kind of thinking. Now, let me get this straight: I don't normally do the web stuff, so i don't have a whole lot of opinion about CSS as such. For what opinion I have, I agree with Ron that CSS sucks and personally I like to see everything in one file without artificially splitting it. Nothing is worse than dealing with lots of small files with cross-dependencies. For what it's worth, I hate JavaScript and the web sites written with it, why do they have to be so slow and annoying?

But let's return to the subject, about the simplicity. I think there is more than one failure on Ron's part.

First of all, he seems to miss that everyone has a different definition of simplicity. And that applies not only to him but to many other people, and not only to the simplicity but to many other subjects. There are many people who like to quote the Principle of the Least Surprise. But then they go and propose the solutions that surprise me immensely.

My favorite example is a poll where people were asked, whether they think that they're better drivers than average, and more than 80% of them answered "yes". And then the article about that poll went on to talk about how people overestimate themselves. What the article failed to appreciate though is that everyone has a different definition of what a "good driver" is. Some people think that it means driving faster, some think that it means driving slower, some think that showing the turn signals is important, some don't know about the existence of the turn signals. Sure enough, most people strive to be good drivers by their own definition, and really become the good drivers by their own definition. The remaining 20% must have been the really lousy drivers.

So in my reckoning C++ and Perl are the simple languages. They allow for simple and easy writing and reading of the code. Yeah, C++ has some strange points, like why not make these virtual destructors mandatory if you define the virtual methods? And the early C++ really was a pain. But after the addition of templates it became awesome. It's not that you want to use the templates for everything, and you specifically don't want to use the stuff from the Alexandrescu book in any actual project. And they really need the explicit loops instead of the stupid functional recursion, and a better error diagnostics. But when they're used right, they're awesome. The modern C++ really combines the simplicity aspects of C and Perl. And Perl is an awesome simple language. It's driven by how to make things easy and convenient to use.

See, there are at least two definitions of simplicity. One definition thinks that simplicity equals logic and consistency, that everything that is logical and consistent is simple. The other definition thinks that things should be easy to use. The typical problems need to be solvable with minimal effort, and if that requires inconsistency, so be it, the consistency is of not much use by itself. If we look at the classic programming languages, C is convenient while Pascal is all-consistent (and oh so much pain to use). Perl is convenient, Python is consistent (and oh so much pain to use, have you ever tried to read the code in Python?). C++ is convenient, Java is consistent (and is just incredibly painful to use, it's like raising the sun with manual labor). Shell is convenient, PowerShell is consistent (and well, you've got my drift).

Second, the wrong kind of simplicity is bad. During my time at Google, I've actually heard with my own ears someone say: "The system has got so many features to it that it became difficult to modify, so we had to throw it away and rewrite from scratch without the extra features". This is I think a pinnacle of how things can go wrong. Yes, a system without many features is much easier to write. And probably keeping just 80% of the features would make is much easier to maintain. But there is a good reason for the features, they make the system easy to use. Nobody needs a system with 80% of the features because it becomes unusable. And well, if someone can't understand and maintain it, well, they're probably the wrong people to do the job.

That said, of course there is a limit: there are features that make the life of the users convenient and the features for the sake of features. The first kind is good, the second kind is bad, and telling which is which may sometimes be pretty hard. It might not even be that someone (er, marketing) had included the features for the sake of features, they might have been a good idea at some point and then got supplanted by the better ways. But things evolve, and sometimes there is time for some features to die. But not wholesale and not by the principle of simplicity of implementation or for some dumb consistency. The best and most simple-to-use features are often difficult to implement and full of quirks. But that's because the tasks they're solving are difficult and full of quirks, and the features make these tasks easy by covering the complexity and adopting to the quirks.

On performace of reading the unaligned data

I've made a template that was supposed to optimize the reading of the unaligned data:

template <typename T>
T getUnaligned(const T *ptr)
{      
    if ((intptr_t)ptr & (sizeof(T)-1)) {
        T value;
        memcpy(&value, ptr, sizeof(T));
        return value;
    } else {
        return *ptr;
    }
}

And I've been wondering if it really makes a difference. So I've made a test.The test compares 5 cases:

Just read the aligned memory directly;
Read the aligned memory through memcpy;
Read the aligned memory through memcpy;
Read the aligned memory through the macro;
Read the unaligned memory through the macro.

All the cases follow this general pattern:

UTESTCASE readAlignedDirect(Utest *utest)
{
    int64_t n = findRunCount();
    int64_t x = 0;

    double tstart = now();
    for (int64_t i = 0; i < n; i++) {
        x += apdata[i % ARRAYSZ];
    }
    double tend = now();
    printf("        [%lld] %lld iterations, %f seconds, %f iter per second\n", (long long) x, (long long)n, (tend-tstart), (double)n / (tend-tstart));
}


The results on an AMD64 machine with -O3 came out like this:

   case 1/6: warmup
        [0] 4294967296 iterations, 4.601818 seconds, 933319757.968194 iter per second
  case 2/6: readAlignedDirect
        [0] 4294967296 iterations, 4.579874 seconds, 937791527.916276 iter per second
  case 3/6: readAlignedMemcpy
        [0] 4294967296 iterations, 4.537628 seconds, 946522579.007429 iter per second
  case 4/6: readUnalignedMemcpy
        [0] 4294967296 iterations, 6.215574 seconds, 691000961.046305 iter per second
  case 5/6: readAlignedTemplate
        [0] 4294967296 iterations, 4.579680 seconds, 937831317.452678 iter per second
  case 6/6: readUnalignedTemplate
        [0] 4294967296 iterations, 5.052706 seconds, 850033049.741168 iter per second

Even more interestingly, the code generated for readAlignedDirect(), readAlignedTemplate() and readUnalignedTemplate() is exactly the same:

.L46:
    movq    %rax, %rdx
    addq    $1, %rax
    andl    $15, %edx
    addq    (%rcx,%rdx,8), %rbx
    cmpq    %rbp, %rax
    jne .L46

The code generated for readAlignedMemcpy() and readUnalignedMemcpy() is slightly different but very similar:

.L34:
    movq    %rax, %rdx
    addq    $1, %rax
    andl    $15, %edx
    salq    $3, %rdx
    addq    (%rcx), %rdx
    movq    (%rdx), %rdx
    addq    %rdx, %rbx
    cmpq    %rax, %rbp
    movq    %rdx, (%rsi)
    jg  .L34

Basically, the compiler is smart enough to know that the machine doesnt' care about alignment, and generates memcpy() as a plain memory read, and then for condition in the template it's also smart enough to know that both branches end up with the same code and eliminate the condition check.

Smart, huh? So now I'm not sure if keeping this template as it is makes sense, or should I just change it to a plain memcpy().

BTW,  I must say it took me a few attempts to make up the code that would not get optimized out by the compiler. I've found out that memcpy() is not suitable for reading the volatile data nowadays. It defines its second argument as (const void *), and so the volatility has to be cast away for passing the pointer to it, and then the compiler just happily takes the whole memory read out of the loop. And it was even smart enough to sometimes replace the whole loop with a multiplication instead of collecting a sum. That's where the modulo operator came handy to confuse this extreme smartness.

A separate interesting note is that each iteration of the loop executes in an about one-billionth of a second. Since this is a 3GHz CPU, this means that every iteration takes only about 3 CPU cycles.  That's 6 or 10 instructions executed in 3 CPU cycles. Such are the modern miracles.

ThreadedClient improvement of timeout handling

The work on updating the tests for the new ordered index brought a nice side effect: the handling of the timeouts in the Perl class X:ThreadedClient has improved. Now it returns not only the error message but also the data received up to that point. This helps a lot with diagnosing the errors in the automated tests of TQL: the error messages get returned in a clear way.

Ordered Index implemented in C++

The ordered index implemented in Perl has been pretty slow, so I've finally got around to implementing one in C++. It became much faster. Here is an excerpt from the performance test results (with an optimized build):

Table insert makeRowArray (single hashed idx, direct) 0.551577 s, 181298.48 per second.
  Excluding makeRowArray 0.311936 s, 320578.44 per second.
Table insert makeRowArray (single ordered int idx, direct) 0.598462 s, 167095.09 per second.
  Excluding makeRowArray 0.358821 s, 278690.37 per second.
Table insert makeRowArray (single ordered string idx, direct) 1.070565 s, 93408.64 per second.
  Excluding makeRowArray 0.830924 s, 120347.91 per second.
Table insert makeRowArray (single perl sorted idx, direct) 21.810374 s, 4584.97 per second.
  Excluding makeRowArray 21.570734 s, 4635.91 per second.
Table lookup (single hashed idx) 0.147418 s, 678342.02 per second.
Table lookup (single ordered int idx) 0.174585 s, 572785.77 per second.
Table lookup (single ordered string idx) 0.385963 s, 259092.38 per second.
Table lookup (single perl sorted idx) 6.840266 s, 14619.32 per second.

Here the "ordered" is the new ordered index in C++ (on an integer field and on a string field), and "perl sorted" is the ordering in Perl, on an integer field. Even though the ordered index is a little slower than the hashed index, it's about 40-70 times faster than the ordering implemented in Perl.

Naturally, ordering by a string field is slower than by an integer one, since it not only has to compare more bytes one-by-one but also honors the locale-specific collation order.

In Perl, the ordered index type is created very similarly to the hashed index type:

Triceps::IndexType->newOrdered(key => [ "a", "!b" ])

The single option "key" specifies the array with the names key fields of the index. The "!" in front of the field name says that the order by this field is descending, and otherwise the order is ascending. This is more compact and arguably easier-to-read format than the one used by the SimpleOrderedIndex.

As usual, the NULL values in the key fields are permitted, and are considered less than any non-NULL value.

The array fields may also be used in the ordered indexes.

The ordered index does return the list of key fields with getKey(), so it can be used in joins and can be found by key fields with TableType::findIndexPathForKeys() and friends, just like the hashed index.

The getKey() for an ordered index returns the list of plain field names, without the indications of descending order. Thus the finding of the index by key fields just works out of the box, unchanged.

To get the key fields as they were specified in the index creation, including the possible "!", use the new method:

@fields = $indexType->getKeyExpr();

The idea of this method is that the contents of the array returned by it depends on the index type and is an "expression" that can be used to build another instance of the same index type. For the hashed index it simply returns the same data as getKey(). For the ordered index it returns the list of keys with indications. For the indexes with Perl conditions it still returns nothing, though in the future might be used to store the condition.

In the C++ API the index type gets created with the constructor:

OrderedIndexType(NameSet *key = NULL);
static OrderedIndexType *make(NameSet *key = NULL);

The key can also be set after construction:

OrderedIndexType *setKey(NameSet *key);

Same as in Perl, the field names in the NameSet can be prefixed with a "!" to specify the descending order on that field.

The new method in the IndexType is:

virtual const NameSet *getKeyExpr() const;

And the new constant for the index id is IT_ORDERED, available in both C++ and Perl.

more of performance numbers, with optimization

I've realized that the previous published  performance numbers were produced by a build without optimization. So I've enabled the optimization with -O3 -fno-strict-aliasing, and the numbers improved, some of them more than twice (that's still an old Core 2 Duo 3GHz laptop):

Performance test, 100000 iterations, real time.
Empty Perl loop 0.005546 s, 18030711.03 per second.
Empty Perl function of 5 args 0.031840 s, 3140671.52 per second.
Empty Perl function of 10 args 0.027833 s, 3592828.57 per second.
Row creation from array and destruction 0.244782 s, 408526.42 per second.
Row creation from hash and destruction 0.342986 s, 291557.00 per second.
Rowop creation and destruction 0.133347 s, 749922.94 per second.
Calling a dummy label 0.047770 s, 2093352.57 per second.
Calling a chained dummy label 0.049562 s, 2017695.16 per second.
  Pure chained call 0.001791 s, 55827286.04 per second.
Calling a Perl label 0.330590 s, 302489.48 per second.
Row handle creation and destruction 0.140406 s, 712218.40 per second.
Repeated table insert (single hashed idx, direct) 0.121556 s, 822665.80 per second.
Repeated table insert (single hashed idx, direct & Perl construct) 0.337209 s, 296551.58 per second.
  RowHandle creation overhead in Perl 0.215653 s, 463707.00 per second.
Repeated table insert (single sorted idx, direct) 1.083628 s, 92282.62 per second.
Repeated table insert (single hashed idx, call) 0.153614 s, 650981.13 per second.
Table insert makeRowArray (single hashed idx, direct) 0.553364 s, 180713.02 per second.
  Excluding makeRowArray 0.308581 s, 324063.65 per second.
Table insert makeRowArray (double hashed idx, direct) 0.638617 s, 156588.31 per second.
  Excluding makeRowArray 0.393835 s, 253913.40 per second.
  Overhead of second index 0.085254 s, 1172969.41 per second.
Table insert makeRowArray (single sorted idx, direct) 22.355793 s, 4473.11 per second.
  Excluding makeRowArray 22.111011 s, 4522.63 per second.
Table lookup (single hashed idx) 0.142762 s, 700466.78 per second.
Table lookup (single sorted idx) 6.929484 s, 14431.09 per second.
Lookup join (single hashed idx) 2.944098 s, 33966.26 per second.
Nexus pass (1 row/flush) 0.398944 s, 250661.96 per second.
Nexus pass (10 rows/flush) 0.847021 s, 1180608.79 per row per second.
  Overhead of each row 0.049786 s, 2008583.51 per second.
  Overhead of flush 0.349157 s, 286403.84 per second.

I've also tried to run the numbers in a newer laptop with 2GHz Core i7 CPU, in a VM configured with 2CPUs. On the build without the optimization, the numbers came out very similar to the old laptop. On the build with optimization they came up to 50% better but not consistently so (perhaps, running in a VM added variability).

Triceps 2.0.1 released

This time I'm trying to update the docs following the code changes, and do the small releases.

The release 2.0.1 is here! It includes:

• Fixed the version information that was left incorrect (at 0.99).
• Used a more generic pattern in tests for Perl error messages that have changed in the more recent versions of Perl (per CPAN report #99268).
• Added the more convenient way to wrap the error reports in Perl, Triceps::nestfess() and Triceps::wrapfess().
• Added functions for the nicer printing of auto-generated code, Triceps::alignsrc() and Triceps::numalign().
• In the doc chapter on the templates, fixed the output of the examples: properly interleaved the inputs and outputs.

formatting the source code snippets

I've got to printing the snippets of the source code in the error messages when the said code fails to compile and realized that it could be printed better. So I wrote a couple of helper functions to print it better.

First of all, it could be indented better, to match the indenting of the rest of the error message.

Second (but connected), it could be left-aligned better:  the code snippets tend to have the extra spacing on the left.

And third, it could use the line numbers, to make the error location easier to find.

The magic is done with two functions, Triceps::Code::alignsrc() and Triceps::Code::numalign(). They work in the exact same way, have the same arguments, only numalign() adds the line numbers while alignsrc() doesn't.

Here is an example of use:

        confess "$myname: error in compilation of the generated function:\n  $@function text:\n"
        . Triceps::Code::numalign($gencode, "  ") . "\n";

It can produce an error message like this (with a deliberately introduced syntax error):

Triceps::Fields::makeTranslation: error in compilation of the generated function:
  syntax error at (eval 27) line 13, near "})
"
function text:
     2 sub { # (@rows)
     3   use strict;
     4   use Carp;
     5   confess "template internal error in Triceps::Fields::makeTranslation: result translation expected 1 row args, received " . ($#_+1)
     6     unless ($#_ == 0);
     7   # $result_rt comes at compile time from Triceps::Fields::makeTranslation
     8   return $result_rt->makeRowArray(
     9     $_[0]->get("one"),
    10     $_[0]->get("two"),
    11   );
    12 })
 at /home/babkin/w/triceps/trunk/perl/Triceps/blib/lib/Triceps/Fields.pm line 219
    Triceps::Fields::makeTranslation('rowTypes', 'ARRAY(0x2943cb0)', 'filterPairs', 'ARRAY(0x2943b30)', '_simulateCodeError', 1) called at t/Fields.t line 205
    eval {...} called at t/Fields.t line 204


The first argument of alignsrc() and numalign() is the code snippet string to be aligned. The following will be done to it:

• The empty lines at the front will be removed. numalign() is smart enough to take the removed lines into account and adjust the numbering. That's why the numbering in the error message shown above starts with 2. You can also get the number of the removed lines afterwards, from the global variable $Triceps::Code::align_removed_lines.
• The \n at the end of the snippet will be chomped. But only one, the rest of the empty lines at the end will be left alone.
• Then the "baseline" indenting of the code will be determined by looking at the first three and last two lines. The shortest non-empty indenting will be taken as the baseline. If some lines examined start with spaces and some start with tabs, the lines starting with tabs will be preferred as the baseline indenting.
• The baseline indenting will be removed from the front of all lines. If some lines in the middle of the code have a shorter indenting, they will be left unchanged.
• The tabs will be replaced by two spaces each. If you prefer a different replacement, you can specify it as the third argument of the function.
• In numalign() the line numbers will be prepended to the lines.
• The indenting from the second argument of the function will be prepended to the lines.

The second argument  contains the new indenting that allows to align the code nicely with the rest of the error message. Technically, it's optional, and will default to an empty string.

The third argument is optional and allows to provide an alternative replacement to the tab characters in the code. If it's an empty string, it will revert to the default two spaces "  ".  To keep the tabs unchanged, specify it as "\t".

error reporting in the templates

When writing the Triceps templates, it's always good to make them report any usage errors in the terms of the template (though the extra detail doesn't hurt either). That is, if a template builds a construction out of the lower-level primitives, and one of these primitives fail, the good approach is to not just pass through the error from the primitive but wrap it into a high-level explanation.

This is easy to do if the primitives report the errors by returning them directly, as Triceps did in the version 1. Just check for the error in the result, and if an error is found, add the high-level explanation and return it further.

It becomes more difficult when the errors are reported like the exceptions, which means in Perl by die() or confess(). The basic handling is easy, there is just no need to do anything to let the exception propagate up, but adding the extra information becomes difficult. First, you've got to explicitly check for these errors by catching them with eval() (which is more difficult than checking for the errors returned directly), and only then can you add the extra information and re-throw. And then there is this pesky problem of the stack traces: if the re-throw uses confess(), it will likely add a duplicate of at least a part of the stack trace that came with the underlying error, and if it uses die(), the stack trace might be incomplete since the native XS code includes the stack trace only to the nearest eval() to prevent the same problem when unrolling the stacks mixed between Perl and Triceps scheduling.

Because of this, some of the template error reporting got worse in Triceps 2.0.

Well, I've finally come up with the solution. The solution is not even limited to Triceps, it can be used with any kind of Perl programs. Here is a small example of how this solution is used, from Fields::makeTranslation():

    my $result_rt = Triceps::wrapfess
        "$myname: Invalid result row type specification:",
        sub { Triceps::RowType->new(@rowdef); };

The function Triceps::wrapfess() takes care of wrapping the confessions. It's very much like the try/catch, only it has the hardcoded catch logic that adds the extra error information and then re-throws the exception.

Its first argument is the error message that describes the high-level problem. This message will get prepended to the error report when the error propagates up (and the original error message will get a bit of extra indenting, to nest under that high-level explanation).

The second argument is the code that might throw an error, like the try-block. The result from that block gets passed through as the result of wrapfess().

The full error message might look like this:

Triceps::Fields::makeTranslation: Invalid result row type specification:
  Triceps::RowType::new: incorrect specification:
    duplicate field name 'f1' for fields 3 and 2
    duplicate field name 'f2' for fields 4 and 1
  Triceps::RowType::new: The specification was: {
    f2 => int32[]
    f1 => string
    f1 => string
    f2 => float64[]
  } at blib/lib/Triceps/Fields.pm line 209.
    Triceps::Fields::__ANON__ called at blib/lib/Triceps.pm line 192
    Triceps::wrapfess('Triceps::Fields::makeTranslation: Invalid result row type spe...', 'CODE(0x1c531e0)') called at blib/lib/Triceps/Fields.pm line 209
    Triceps::Fields::makeTranslation('rowTypes', 'ARRAY(0x1c533d8)', 'filterPairs', 'ARRAY(0x1c53468)') called at t/Fields.t line 186
    eval {...} called at t/Fields.t line 185

It contains both the high-level and the detailed description of the error, and the stack trace.

The stack trace doesn't get indented, no matter how many times the message gets wrapped. wrapfess() uses a slightly dirty trick for that: it assumes that the error messages are indented by the spaces while the stack trace from confess() is indented by a single tab character. So the extra spaces of indenting are added only to the lines that don't start with a tab.

Note also that even though wrapfess() uses eval(), there is no eval above it in the stack trace. That's the other part of the magic: since that eval is not meaningful, it gets cut from the stack trace, and wrapfess() also uses it to find its own place in the stack trace, the point from which a simple re-confession would dump the duplicate of the stack. So it cuts the eval and everything under it in the original stack trace, and then does its own confession, inserting the stack trace again. This works very well for the traces thrown by the XS code, which actually doesn't write anything below that eval; wrapfess() then adds the missing part of the stack.

Wrapfess() can do a bit more. Its first argument may be another code reference that generates the error message on the fly:

    my $result_rt = Triceps::wrapfess sub {
            "$myname: Invalid result row type specification:"
        },
        sub { Triceps::RowType->new(@rowdef); };

In this small example it's silly but if the error diagnostics is complicated and requires some complicated printing of the data structures, it will be called only if the error actually occurs, and the normal code path will avoid the extra overhead.

It gets even more flexible: the first argument of wrapfess() might also be a reference to a scalar variable that contains either a string or a code reference. I'm not sure yet if it will be very useful but it was easy to implement. The idea there is that it allows to write only one wrapfess() call and then change the error messages from inside the second argument, providing different error reports for its different sections. Something like this:

    my $eref;
    return Triceps::wrapfess \$eref,
        sub {
          $eref = "Bad argument foo";
          buildTemplateFromFoo();
          $eref = sub {
            my $bardump = $argbar->dump();
            $bardump =~ s/^/    /mg;
            return "Bad argument bar:\n  bar value is:\n$bardump";
           }
          buildTemplateFromBar();
          ...

       };

It might be too flexible, we'll see how it works.

Internally, wrapfess() uses the function Triceps::nestfess() to re-throw the error. Nestfess() can also be used directly, like this:

eval {
  buildTemplatePart();
};
if ($@) {

  Triceps::nestfess("High-level error message", $@);
}

The first argument is the high-level descriptive message to prepend, the second argument is the original error caught by eval. Nestfess() is responsible for all the smart processing of the indenting and stack traces, wrapfess() is really just a bit of syntactic sugar on top of it.

nested record and ORM

Some more thoughts on the subject. First of all, why should anyone care about the nested records?

Nowadays people are used to the relational model. But next to all of this sits Google with its protobufs. It has published the information about protobufs, and about some of its systems that use protobufs but I think there is not enough appreciation of how much the protobufs mean for Google. If you look at the published information, you'll realize that everything at Google, totally everything, is based on protobufs. They're used as the arguments of the RPC calls and as the returned values, the map-reduce is built on them, as well as all kinds of the no-SQL DBMS and a few of SQL DBMS. And there are good reasons for why they're used so much. Well, for the no-SQL databases it's a no-brainer: if all it can do is store the key-value pairs, at least the values must be capable of representing complexity. To understand the other reasons, let's look at the difference between the nested records and the Object Relation Managers.

Entirely coincidentally, I've encountered recently a reference to the article at http://blogs.tedneward.com/2006/06/26/The%2BVietnam%2BOf%2BComputer%2BScience.aspx, on the subject of ORM. It's an interesting viewpoint though the conclusions don't seem quite right to me, as well as a few of the pre-suppositions.

The idea that the object model has the concept of identity and the relational model doesn't, doesn't seem right. In reality the primary keys are used all over the real databases. And the primary keys are exactly that, the identity. Even more so if the primary key is not a part of data but an auto-generated unique value. Here we can go into the discussions about "but it's just a hack, not a part of the theory" but think about this: these auto-generated unique keys are the bread and butter of the data normalization, which certainly is the part of the theory. The point of normalization is to get rid of the duplicate data in the rows and move it to a separate table where it will exist as a single copy, with multiple short references to it from the original table. But so are the references between the objects, multiple objects may refer to the same subordinate object. And the reference to that object's address (its identity) is just the same thing as the references to that unique key in the subordinate table.

So, what ORM does is map one-to-one the objects and records, and maps the references between the objects by address to these references between the records by the unique id values. The joins are then the way to navigate the references between the objects.

The nested records in protobufs (or their equivalents) have a different meaning and are used differently. The whole point of a nested record is that it represents a self-contained unit of data. This data is all stored together because it's normally used together. It's obvious to see when the protobufs are used for RPC but it works the same when the protobufs are used in a database. Obviously, not everything in the world can be covered by this usage and obviously there are downsides. An obvious downside is that a nested record  in a protobuf can't be modified piecemeal: if you need to change it, you have to build a whole new record and put it back in place. And obviously the data is stored by value, not by reference, so the duplicate data has to be duplicated. But for the situations where the self-containment is real important, they work great. As in all the massively-parallel processing at Google. Each nested record has enough information in it to be processed separately from the others. Again, not everything in the world can be expressed like this, but what can, benefits tremendously from the nested records.

Now I think I can formulate, what is the difference between the ORM and the nested records: the ORMs connect the objects by reference while the nested records bunch up the data by value. Different uses, different trade-offs, different implementations. (Though for the in-memory implementations the boundary is more blurred).

And I think I've found the explanation to why they say that the Oracle's nested records are slower than the manual joins, even though they amount to joins, the the ORM-like way. Well, or course it could just be a limitation of the query optimizer: if the optimizer runs first and only then the nesting is expanded, it obviously would not be as optimal. But that should be easy to fix, just expand the nesting and then run the optimizer, so this is unlikely to be the problem. A more fundamental problem seems to be that the records are sent into the database with the data by value, and then they get broken up into the flat tables and stored with connections by reference. It successfully combines the data duplication of by-value and the cost of the extra joins from by-reference. Both kinds of shortcomings rolled together.

Coming back to that point that not everything in the world can be represented as the nested records in protobufs. This means that the references between the tables (or objects, if you prefer to think about them that way) are still needed, and the nested records still need to use them. Which means fitting them into the relational model. And that's where my whitepaper comes in: how to represent the nested records in the relational model. They can be logically dis-joined into the "underlying tables" and the relational operations can be applied to them in the usual way.

Which can also solve problem if the representation of the object inheritance and polymorphism in ORMs, described in the article I listed above. Well, representing the polymorphic data as nested records is fairly straightforward as long as the records contain the indication of subtypes.  And as I've shown in the whitepaper, mapping the nested records to the relational use is also straightforward.  Put them together and you have the solution.

more on nested records (and Oracle)

As people have pointed out to me, Oracle has the concept of "nested tables" and "variable arrays" (VARRAY). Since at least the version 9i, a very long time ago.

The nested tables have a difference in how they store the data: they break out the repeating sections into a separate table and automatically create the id fields to link them together. So the select just does an implicit join. It's like an embedded ORM (by the way, there is also a concept of references to rows, also apparently implemented in an ORM way). Strangely, they say that the nested tables are slower than the usual joins. Not sure, why.

The VARRAYs are more like the real nested structures, they store the data directly in the row, and they preserve the order of nested elements. However the maximal size of VARRAYs has to be declared in advance.

Obviously, the nested tables (the ORM version) would be much faster on updates, they can change the nested parts very easily without changing the rest of the records. And also faster on scanning of the non-repeating parts.

The version with the repeated parts stored directly in the records, such as protobufs or varrays should be faster on reading these repeated parts.

The syntax for accessing the nested parts is like this:

select
  upper.key,
  nested.*
from
  main_table upper,
  table(upper.nested) nested;

So you can select the nested fields as a single field of a record type or you can flatten it as shown above. And the joins within the same top-level records should be possible to express in the same syntax like this:

from
  main_table upper,
  table(upper.nested1) nested1,
  table(upper.nested2) nested2;

There also are the examples of deleting the nested parts, like:

delete from table (
  select nested from upper where key = 1234
) nested
where
  nested.field = 5678;

Or even by the whole value of the nested part.

And by the same example, you can insert a nested part into the record, or update one.

But the part that seems to be missing is the un-flattening. Can you select data from one table with nested data, re-format the nested parts, create different nested records and either return them from select or insert them into another table? I haven't found any examples of that. It looks like the nested records aren't the first-class citizen in Oracle either.