%C = (1=>.5, 0=>.5);
%{$S{1}} = (1=>.1, 0=>.9);
%{$S{0}} = (1=>.5, 0=>.5);
%{$R{1}} = (1=>.8, 0=>.2);
%{$R{0}} = (1=>.2, 0=>.8);
%{$W{"11"}} = (1=>.99, 0=>.01);
%{$W{"10"}} = (1=>.9, 0=>.1);
%{$W{"01"}} = (1=>.9, 0=>.1);
%{$W{"00"}} = (1=>0, 0=>1);

$dist = "P(C|SR)";
foreach $s (0..1) {
  foreach $r (0..1) {
    $cond = "$s$r";
    foreach $c (0..1) {
      $val = "$c";
      $joint = $C{$c}*$S{$c}{$s}*$R{$c}{$r};
      $Distr{$dist}{$cond}{$val}=$joint;
      $Tot{$dist}{$cond}+=$joint;
    }
  }
}

$dist = "P(S|CRW)";
foreach $c (0..1) {
  foreach $r (0..1) {
    foreach $w (0..1) {
      $cond = "$c$r$w";
      foreach $s (0..1) {
        $val = "$s";
        $joint = $C{$c}*$S{$c}{$s}*$R{$c}{$r}*$W{"$s$r"}{$w};
        $Distr{$dist}{$cond}{$val}=$joint;
        $Tot{$dist}{$cond}+=$joint;
      }
    }
  }
}

$dist = "P(R|CSW)";
foreach $c (0..1) {
  foreach $s (0..1) {
    foreach $w (0..1) {
      $cond = "$c$s$w";
      foreach $r (0..1) {
        $val = "$r";
        $joint = $C{$c}*$S{$c}{$s}*$R{$c}{$r}*$W{"$s$r"}{$w};
        $Distr{$dist}{$cond}{$val}=$joint;
        $Tot{$dist}{$cond}+=$joint;
      }
    }
  }
}

$dist = "P(W|SR)";
foreach $s (0..1) {
  foreach $r (0..1) {
    $cond = "$s$r";
    foreach $w (0..1) {
      $val = "$w";
      $joint = $W{"$s$r"}{$w};
      $Distr{$dist}{$cond}{$val}=$joint;
      $Tot{$dist}{$cond}+=$joint;
    }
  }
}

# normalize
foreach $dist (keys %Distr){
  foreach $cond (keys %{$Distr{$dist}}){
    foreach $val (0..1){
      $Distr{$dist}{$cond}{$val}/=$Tot{$dist}{$cond};
    }
  }
}

print "\n";
foreach $dist (keys %Distr){
  print "$dist:\n";
  foreach $cond (keys %{$Distr{$dist}}){
    print "$cond :";
    foreach $val (0..1){
      printf( " %0.12f",$Distr{$dist}{$cond}{$val} );
    }
    print "\n";
  }
  print "\n";
}

$s=1;
$w=1;
foreach $c (0..1){
  foreach $r (0..1){
    $Distr{"P(r|sw)"}{"$s$w"}{$r}+=$C{$c}*$S{$c}{$s}*$R{$c}{$r}*$W{"$s$r"}{$w};
    $Tot{"P(r|sw)"}{"$s$w"}      +=$C{$c}*$S{$c}{$s}*$R{$c}{$r}*$W{"$s$r"}{$w};
  }
}
printf "FJD result = %0.12f\n",$Distr{"P(r|sw)"}{"$s$w"}{1}/$Tot{"P(r|sw)"}{"$s$w"};


foreach $c (0..1){
  foreach $s (0..1){
    foreach $r (0..1){
      foreach $w (0..1){
        printf " $c$s$r$w ";
      }
    }
  }
}
print "\n";
foreach $c (0..1){
  foreach $s (0..1){
    foreach $r (0..1){
      foreach $w (0..1){
        printf "%0.3f ",$C{$c}*$S{$c}{$s}*$R{$c}{$r}*$W{"$s$r"}{$w};
      }
    }
  }
}
print "\n";


$N=0;
$D=0;
$c=(((rand)<($C{1}))?1:0);
$s=1;
$r=(((rand)<($R{$c}{1}))?1:0);
$w=1;
#$w=(rand(1.0)>$W{"$s$r"}{0});
print "$c $s $r $w\n";
foreach $t (0..100000){
  $c=(((rand)<($Distr{"P(C|SR)"}{"$s$r"}{1}))?1:0);
  $r=(((rand)<($Distr{"P(R|CSW)"}{"$c$s$w"}{1}))?1:0);
  $N+=((($r==0)?0:1));
  $D+=(1);
}
printf "MCMC result = %0.12f\n",$N/$D;