diffNODE.001_01

#!/usr/bin/perl -w

use strict;
use FindBin qw ($Bin);
use lib $Bin;
use FileHandle;
#Frame was inlined in this script
#use Tools::Frame;
use Data::Dumper;
use POSIX qw (floor);
use Getopt::Long;

my %opts = (
            'spnorms' => 'VORTICITY,DIVERGENCE,TEMPERATURE,KINETIC ENERGY',
            'gpnorms' => '',
            'norm-max-diff' => .050000,
            'jo-max-diff'   => .030000,
           );

sub frame
{
  my ($t, $width, $border) = @_;
  my $len = length ($t);

  $border = 1 unless (defined ($border));

  my $S = ' ';

  my ($C, $V, $H) = ('*', '|', '-');

  ($C, $V, $H) = (' ', ' ', ' ')
    unless ($border);

  my $df = 3;
  $width ||= $len + 2 * $df;

  my $line1 = $C . ($H x ($width-2)) . $C;
  my $line2 = $V . ($S x ($width-2)) . $V;


  my $TEXT = '';

  $TEXT .= "$line1\n";
  for (1 .. ($df-1)/2)
    {
      $TEXT .= "$line2\n";
    }

  die ("Cannot frame text: `$t'\n")
    if ($width - 2 * $df <= 0);
  
  while ($t)
    {
      my $s = substr ($t, 0, $width - 2 * $df, '');

      my $i = 0;
      while (length ($s) < $width - 2 * $df)
        {
          if ($i % 2)
            {
              $s = " $s";
            }
          else
            {
              $s = "$s ";
            }
          $i++;
        }
      my $linet = $V . ($S x ($df-1)) . $s .  ($S x ($df-1)) . $V;
      $TEXT .= "$linet\n";
    }

  for (1 .. ($df-1)/2)
    {
      $TEXT .= "$line2\n";
    }
  $TEXT .= "$line1\n";
}


sub xave
{
  my $f = shift;
  my $fh = 'FileHandle'->new ("<$f");

  my @gpregs;

  if ($opts{gpnorms})
    {
      if ((scalar (@{$opts{gpnorms}}) == 1) && ($opts{gpnorms}[0] eq '*'))
        {
          @gpregs = (qr/^\s*GPNORM\s+\b(\S+)\b/);
        }
      else
        {
          @gpregs = map { qr/^\s*GPNORM\s+\b(\Q$_\E)\b/ } @{$opts{gpnorms}};
        }
    }

  my @x;
  MAIN: while (defined (my $line = <$fh>))
    {
      AGAIN:

      if ($line =~ m/^\s*GPNORM\s+/o)
        {
          for my $gpreg (@gpregs)
            {
              if ($line =~ $gpreg)
                {
                  my $F = $1;
                  $line = <$fh>;
                  next MAIN  unless ($line =~ m/AVE\s+0/o);
                  for ($line)
                    {
                      s/^\s*AVE\s*//o; 
                      s/\s+/\n/go; 
                    }
                  push @x, map { [ $F, $_ ] } split (m/\n/o, $line);
                  next MAIN;
                }
            }
        }

      if ($line =~ s/^\s*SPECTRAL\s+NORMS\s+-\s+//o)
        {
          AGAIN_SPNORMS:


          goto AGAIN
            unless (($line = <$fh>) =~ s/^\s+LEV\s+//o);

          my %index;
          %index = ();
          for my $spnorm (@{$opts{spnorms}})
            {
              my $index = index ($line, $spnorm);
              $index{$spnorm} = $index 
                if ($index >= 0);
            }

          my @spnormk = sort { $index{$a} <=> $index{$b} } 
                        grep { defined $index{$_} } 
                        @{$opts{spnorms}};

          goto AGAIN
            unless (($line = <$fh>) =~ s/^\s+AVE\s+//o);

          my @spnormv = split (m/\s+/o, $line);

          while (@spnormk)
            {
              my $spnormk = shift (@spnormk);
              my $spnormv = shift (@spnormv);
              die ("$spnormk, $spnormv\n")
                unless (defined ($spnormk) && defined ($spnormv));
              push @x, [ $spnormk, $spnormv ];
            }

          goto AGAIN_SPNORMS;

        }
    }

  return @x;
}

sub xobstype
{
  my $f = shift;
  my $fh = 'FileHandle'->new ("<$f");

  # keep final value for obs number & JO

  my @x;
  while (defined (my $line = <$fh>))
    {
      next unless ($line =~ m/ObsType\s+(\d+)\s+Total:\s*(\d+)\s+(\S+)\s+/o);
      @{$x[$1]}{qw (number JO)} = ($2, $3);
    }

  return @x;
}

sub xjog
{
  my $f = shift;
  my $fh = 'FileHandle'->new ("<$f");

  # keep final value for obs number & JO

  my %x;
  while (defined (my $line = <$fh>))
    {
      next unless ($line =~ m/Jo Global\s*:\s*(\d+)\s+(\S+)/o);
      @x{qw (number JO)} = ($1, $2);
    }

  return \%x;
}

sub center
{
  my ($s, $n) = @_;
  my $i = 0;
  while (length ($s) < $n) 
    {
      $s = $i % 2 ? " $s" : "$s ";
      $i++;
    }
  return $s;
}



&GetOptions 
  ('spnorms=s'       => \$opts{'spnorms'},
   'gpnorms=s'       => \$opts{'gpnorms'},
   'norm-max-diff=s' => \$opts{'norm-max-diff'},
   'jo-max-diff=s'   => \$opts{'jo-max-diff'},
  );

$opts{'spnorms'} = [ split (m/,/o, $opts{'spnorms'}) ];
$opts{'gpnorms'} = [ split (m/,/o, $opts{'gpnorms'}) ];

my ($f1, $f2) = @ARGV;

die ("Usage: $0 NODE.001_01 NODE.001_01.ref\n")
  unless ($f1 && $f2);

my @fx1 = &xave ($f1);
my @fx2 = &xave ($f2);

print &frame ("NORMS DIFFERENCES", 121);

my @x = ([]);
my %diff;
my $zero = 0;
my $numb = 0;

my $tag1 = "NORMDIFF";
my $tag2 = "NORMSTAT";

my $nout = 0;

while (defined (my $fx1 = shift (@fx1)) && defined (my $fx2 = shift (@fx2)))
  {
    my ($f1, $x1) = @$fx1;
    my ($f2, $x2) = @$fx2;

    die ("Field mismatch $f1 != $f2\n")
      unless ($f1 eq $f2);

    chomp ($x1); chomp ($x2);
    if (($x1 !~ m/^\s*$/o) && ($x2 !~ m/^\s*$/o))
      {
        my $dx = $x1 - $x2;
        my $dr = ($x1+$x2 > 0) ? 2*$dx/($x1+$x2) : 0.;

        my $sdx = sprintf ('%17.9e', $dx);
        my $sdr = sprintf ('%17.9e', $dr);

        $dx = $sdx; $dx = $dx + 0.;
        $dr = $sdr; $dr = $dr + 0.;

        push @{$x[-1]},
          sprintf (" $tag1 | %20s | %17.9e  |  %17.9e  |  %17s  |  %17s %s\n", &center ($f1, 20), $x1, $x2, $sdx, $sdr, 
                   $dr > $opts{'norm-max-diff'} ? '*' : '');

        $nout++ 
         if ($dr > $opts{'norm-max-diff'});

        if (abs ($dr) > 0)
          {
            my $n = &floor ((log (abs ($dr)) / log (10)));
            $diff{$n}++;

          }
        else
          {
            $zero++;
          }

        $numb++;
      }
    else
      {
        push @x, [];
      }
  }

printf " $tag1 |                      |%-19s | %-19s | %-19s | %-19s\n", 
        &center ("NORM(REF)", 19), &center ("NORM(EXP)", 19), 
        &center ("NORM(REF)-NORM(EXP)", 19), &center ("(NORM(REF)-NORM(EXP))", 19);

printf " $tag1 |                      |%-19s | %-19s | %-19s | %-19s\n", 
        '', '', '', &center ("/NORM(REF)", 19);

for (my $i = 0; $i <= $#x; $i++)
  {
    last unless (@{$x[$i]});
    print @{$x[$i]};
  }

print "\n";


my $diff_cumul = 0;
my $perc_cumul = 0;
for my $n1 (sort { $a <=> $b } keys (%diff))
  {
    my $n2 = $n1 + 1;
    my $diff = $diff{$n1};
    my $perc = 100 * $diff / $numb;
    $diff_cumul += $diff;
    $perc_cumul += $perc;
    printf (" $tag2 | 1e%+2.2d .. 1e%+2.2d | %3d / %3d | %3d / %3d | %6.2f %%, %6.2f %%\n", $n1, $n2, $diff, $numb, $diff_cumul, $numb, $perc, $perc_cumul);
  }

if ($nout)
  {
    print "\n";
    my $text = sprintf ("WARNING : SOME NORMS DIFFERENCES ARE OUTSIDE ALLOWED LIMIT OF %6.2f %%\n", 100 * $opts{'norm-max-diff'});
    print $text x 5;
  }

print "\n";

my @ot1 = &xobstype ($f1);
my @ot2 = &xobstype ($f2);

my $not1 = scalar (@ot1);
my $not2 = scalar (@ot2);

my $not = $not1 > $not2 ? $not1 : $not2;

goto END
  unless ($not > 0);

my $obs_fmtd  = " %10d |  %13.7e ";
my $obs_fmtds = " %10s |  %13s ";
my $obs_fmtp  = " %20d |      %9.4f %% ";
my $obs_fmtps = " %20s |  %15s ";

printf (" OBS_DIFF | %6s | ", 'Type');
printf ($obs_fmtds, 'NOBS(REF)', 'JO(REF)');
print " | ";
printf ($obs_fmtds, 'NOBS(EXP)', 'JO(EXP)');
print " | ";
printf ($obs_fmtps, 'NOBS(EXP)-NOBS(REF)', 'JO(EXP)-JO(EXP)');
printf ("\n");

my $pot = sub 
  {
    my ($ot, $obs_fmt) = @_;
    $obs_fmt ||= $obs_fmtd;
    (my $blank = sprintf ($obs_fmt, 0, 0)) =~ s/\S/ /go;
    if ($ot)
      {
        printf ($obs_fmt, $ot->{number}, $ot->{JO});
      }
    else
      {
        print $blank;
      }
  };

my $dot = sub
  {
    my ($ot1, $ot2) = @_;

    my $dot = $ot1 && $ot2 ? { number => abs ($ot1->{number} - $ot2->{number}), 
                               JO => 100 * 2 * abs ($ot1->{JO} - $ot2->{JO}) / abs ($ot1->{JO} + $ot2->{JO}), } : undef;
    return $dot;
  };

for my $i (1 .. $not-1)
  {
    my $ot1 = $ot1[$i];
    my $ot2 = $ot2[$i];
    my $dot12 = $dot->($ot1, $ot2);

    printf (" OBS_DIFF | %6d | ", $i);

    $pot->($ot1);
    print " | ";
    $pot->($ot2);
    print " | ";
    $pot->($dot12, $obs_fmtp);
    print "\n";


  }

my $jog1 = &xjog ($f1);
my $jog2 = &xjog ($f2);

if ($jog1 || $jog2)
  {
    my $dot12 = $dot->($jog1, $jog2);
    printf (" OBS_DIFF | GLOBAL | ");
    $pot->($jog1);
    print " | ";
    $pot->($jog2);
    print " | ";
    $pot->($dot12, $obs_fmtp);
    print "\n";

    
    if ($dot12->{JO} > 100 * $opts{'jo-max-diff'})
      {
        my $text = sprintf ("WARNING : GLOBAL JO DIFFERENCE IS OUTSIDE ALLOWED LIMIT OF %12.6f %%\n", 100 * $opts{'jo-max-diff'});
        print "\n";
        print $text x 5;
      }

  }


print "\n";

END:



exit $nout;