lib/model/cwres_module_subs.pm

   1 # {{{ new
   2
   3 start new
   4 {
   5
   6   unless( defined $this -> {'nm_version'} ){
   7     $this -> {'nm_version'} = 'default';
   8   }
   9
  10   my $mirror_name = $this -> {'mirror_plots'} ? 'sim' : '';
  11
  12   my ($nmdir,$version) = split(/,/ , $PsN::config -> { 'nm_versions' } -> { $this -> {'nm_version'} } );
  13
  14   if( not defined $version ){
  15     'debug' -> die( message => "CWRES: No NONMEM version \"".$this -> {'nm_version'}."\" in \"$nmdir\" defined in psn.conf. Format should be: name=directory,version" );
  16   } else {
  17     unless( ($version == '5') or ($version == '6') ){
  18       'debug' -> die( message => "CWRES: unknown NONMEM version: $version" );
  19     }
  20   }
  21
  22   $this -> {'nm_version'} = $version;
  23
  24   # Problem is the modelfile problem we are modifing to compute CWRES.
  25
  26   my $prob = $this -> {'problem'};
  27
  28   # Get number of etas and eps;
  29   my $nthetas = $prob -> record_count( record_name => 'theta' );
  30   my $netas = $prob -> nomegas();
  31   my $neps = $prob -> nsigmas();
  32
  33   # Get current comres number
  34   my $comresno;
  35   my ( $crno_ref, $junk ) = $prob ->
  36       _option_val_pos( name        => 'COMRES',
  37                        record_name => 'abbreviated',
  38                        exact_match => 0 );
  39   if( defined $crno_ref ) {
  40     $comresno = $crno_ref -> [0];
  41   }
  42
  43   # Add $ABBREVIATED if necessary
  44   if ( defined $comresno ) {
  45     $prob -> remove_option( record_name  => 'abbreviated',
  46                             option_name  => 'COMRES' );
  47     $prob -> add_option( record_name  => 'abbreviated',
  48                          option_name  => 'COMRES',
  49                          option_value => ($netas+$neps+$comresno) );
  50   } else {
  51     $prob -> add_records( type => 'abbreviated',
  52                           record_strings => [ "COMRES=".($netas+$neps) ] );
  53   }
  54
  55   # get the table names. They are needed below and further down
  56   my @cwtab_names = @{$this -> cwtab_names};
  57
  58   # Figure out if we have an sdtab and what number it has
  59   my ( $sd_ref, $junk ) = $prob ->
  60       _option_val_pos( name        => 'FILE',
  61                        record_name => 'table',
  62                        exact_match => 0 );
  63   if( defined $sd_ref ) {
  64     foreach my $tabname ( @{$sd_ref} ) {
  65       if( $tabname =~ /[sd,pa]tab(\d+)/i ) {
  66         my $sdno = $1;
  67         if( $sdno eq '' ){
  68           $sdno = 1;
  69         }
  70         $this -> sdno($sdno);
  71         for( my $i = 0; $i <= $#cwtab_names; $i++ ) {
  72
  73           # This regular expression is probably quite unneccessary. It
  74           # matches evertyhing before the first 'dot' in a filename,
  75           # the dot, and the rest of the name(dots included). We can
  76           # then inject a number before the first dot. It also handles
  77           # no dots, the number will then be injected at the end of
  78           # the filename.
  79
  80           if( $cwtab_names[$i] =~ /([^\.]+)(\.{0,1})(.*)/ ) {
  81             $cwtab_names[$i] = $1.$sdno.$mirror_name.$2.$3;
  82           }
  83         }
  84         $this -> cwtab_names( \@cwtab_names);
  85         last;
  86       }
  87     }
  88   }
  89
  90
  91   # Figure out wheter we have and 'ADVAN' option. By not using
  92   # "exact_match" we can search for a prefix of the different ADVAN
  93   # options.
  94
  95   my ($advan,$junk) = $prob -> _option_val_pos( record_name => 'subroutine',
  96                                                 name => 'ADVAN',
  97                                                 exact_match => 0);
  98
  99   my $have_advan = scalar(@{$advan}) > 0;
 100
 101   if( $version == 5 ){
 102     if( $have_advan ){
 103       # infn.f will be written in "post_process"
 104     } else {
 105       my $code = $prob -> preds -> [0] -> verbatim_first;
 106       unless( defined $code ){
 107         $code = [];
 108         $prob -> preds -> [0] -> verbatim_first($code);
 109       }
 110       unshift(@{$code},
 111               # {{{ fortan code
 112               ('"      COMMON /ROCM6/ THETAF(40),OMEGAF(30,30),SIGMAF(30,30)',
 113                '"      COMMON /ROCM7/ SETH(40),SEOM(30,30),SESIG(30,30)',
 114                '"      COMMON /ROCM8/ OBJECT ',
 115                '"      DOUBLE PRECISION THETAF, OMEGAF, SIGMAF ',
 116                '"      DOUBLE PRECISION OBJECT ',
 117                '"      REAL SETH,SEOM,SESIG ',
 118                '"      INTEGER J,I ',
 119                '"      INTEGER MODE ',
 120                '"      INTEGER NTH,NETA,NEPS ',
 121                "\"      DATA NTH,NETA,NEPS/$nthetas,$netas,$neps/ ")
 122               # }}}
 123               );
 124
 125       # Abbrev code
 126       $code = $prob -> preds -> [0] -> code;
 127
 128       # {{{ fortran code
 129       push( @{$code},
 130             ('"      IF (ICALL.EQ.0) THEN',
 131              '"C     open files here, if necessary',
 132              '"         OPEN(50,FILE=\'cwtab'.$this -> sdno().$mirror_name.'.est\')') );
 133
 134       #for( my $i = 0; $i <= $#cwtab_names; $i++ ) {
 135         #push( @{$code}, '"         OPEN(5'.$i.',FILE=\''.$cwtab_names[$i].'\')' );
 136       #}
 137
 138       push( @{$code},
 139               # {{{ fortan code
 140             ('"      ENDIF',
 141              '"      IF (ICALL.EQ.3) THEN',
 142              '"         MODE=0',
 143              '"         CALL PASS(MODE)',
 144              '"         MODE=1',
 145              '"         WRITE (50,*) \'ETAS\'',
 146              '" 20      CALL PASS(MODE)',
 147              '"         IF (MODE.EQ.0) GO TO 30',
 148              '"         IF (NEWIND.NE.2) THEN',
 149              '"            CALL GETETA(ETA)',
 150              '"            WRITE (50,97) (ETA(I),I=1,NETA)',
 151              '"         ENDIF',
 152              '"         GO TO 20',
 153              '" 30      CONTINUE',
 154              '"         WRITE (50,*) \'THETAS\'',
 155              '"         WRITE (50,99) (THETAF(J),J=1,NTH)',
 156              '"         WRITE (50,*) \'OMEGAS\'',
 157              '"         DO 7000 I=1,NETA',
 158              '" 7000       WRITE (50,99) (OMEGAF(I,J),J=1,NETA)',
 159              '"            WRITE (50,*) \'SIGMAS\'',
 160              '"            DO 7999 I=1,NEPS',
 161              '" 7999          WRITE (50,99) (SIGMAF(I,J),J=1,NEPS)',
 162              '"            ENDIF',
 163              '" 99         FORMAT (20E15.7)',
 164              '" 98         FORMAT (2I8)',
 165              '" 97         FORMAT (10E15.7)')
 166
 167             # }}}
 168             );
 169
 170     }
 171   } elsif ( $version == 6 ) {
 172
 173     my $code;
 174
 175     if( $have_advan ){
 176       unless( $prob -> infns ){
 177         $prob -> add_records( type => 'infn',
 178                               record_strings => [] );
 179       }
 180
 181       $code = $prob -> infns -> [0] -> code;
 182
 183     } else {
 184
 185       $code = $prob -> preds -> [0] -> code;
 186
 187     }
 188
 189     push( @{$code},
 190           'IF (ICALL.EQ.3) THEN',
 191           '  OPEN(50,FILE=\'cwtab'.$this -> sdno().$mirror_name.'.est\')',
 192           '  WRITE (50,*) \'ETAS\'',
 193           '  DO WHILE(DATA)',
 194           '    IF (NEWIND.LE.1) WRITE (50,*) ETA',
 195           '  ENDDO',
 196           '  WRITE (50,*) \'THETAS\'',
 197           '  WRITE (50,*) THETA',
 198           '  WRITE (50,*) \'OMEGAS\'',
 199           '  WRITE (50,*) OMEGA(BLOCK)',
 200           '  WRITE (50,*) \'SIGMAS\'',
 201           '  WRITE (50,*) SIGMA(BLOCK)',
 202           'ENDIF' );
 203   }
 204
 205
 206   my $code_records;
 207   if( $have_advan ){
 208     # We have and ADVAN option in $SUBROUTINE, get $ERROR code
 209     $code_records = $prob -> errors;
 210
 211     # If we also use version 5, we must include "infn.f" in $SUBROUTINE
 212     if( $version == 5 ){
 213       $prob -> add_option( record_name => 'subroutine',
 214                            option_name => 'INFN',
 215                            option_value=> 'infn.f' );
 216     }
 217
 218   } else {
 219     # No ADVAN subroutine, we should modify $PRED code
 220     $code_records = $prob -> preds;
 221   }
 222
 223   # Get code array reference, so we can update the code inplace.
 224   my $code = $code_records -> [0] -> verbatim_last;
 225
 226   unless( defined $code ){
 227     $code = [];
 228     $code_records -> [0] -> verbatim_last($code);
 229   }
 230
 231   my $com = defined $comresno ? $comresno + 1 : 1;
 232
 233   my @table_row;
 234
 235   for( 1..$netas ){
 236     push( @{$code},"\"  COM($com)=G($_,1)" );
 237     push( @table_row, "COM($com)=G$_"."1");
 238     $com++;
 239   }
 240
 241   for( 1..$neps ){
 242     if( $have_advan ){
 243       push( @{$code},"\"  COM($com)=HH($_,1)" );
 244       push( @table_row, "COM($com)=H$_"."1" );
 245     } else {
 246       push( @{$code},"\"  COM($com)=H($_,1)" );
 247       push( @table_row, "COM($com)=H$_"."1" );
 248     }
 249     $com++;
 250   }
 251
 252   my ($mdv,$junk) = $prob -> _option_val_pos( record_name => 'input',
 253                                               name => 'MDV' );
 254
 255   if( defined $prob -> preds and scalar(@{$mdv}) == 0  ){
 256     $mdv = '';
 257   } else {
 258     $mdv = 'MDV';
 259   }
 260
 261   $prob -> add_records( type => 'table',
 262                         record_strings => ['ID ',
 263                                            join(' ',@table_row),
 264                                            "IPRED DV $mdv NOPRINT ".
 265                                            "ONEHEADER FILE=cwtab".$this -> sdno().'.deriv'] );
 266 #  $prob -> add_records( type => 'table',
 267 #                       record_strings => ['ID ',
 268 #                                          join(' ',@table_row),
 269 #                                          "IPRED DV $mdv NOPRINT ".
 270 #                                          "ONEHEADER FILE=cwtab".$this -> sdno().'.cwres'] );
 271 }
 272 end new
 273
 274 # }}}
 275
 276 # {{{ post_process
 277 start post_process
 278 {
 279
 280   my ($advan,$junk) = $self -> {'problem'} -> _option_val_pos( record_name => 'subroutine',
 281                                                                name => 'ADVAN',
 282                                                                exact_match => 0);
 283
 284   my $mirror_name = $self -> {'mirror_plots'} ? 'sim' : '';
 285
 286   if( $self -> {'nm_version'} == 5 and scalar(@{$advan}) > 0 ){
 287
 288     my $ntheta = $self -> {'problem'} -> record_count( record_name => 'theta' );
 289     my $neta = $self -> {'problem'} -> nomegas;
 290     my $neps = $self -> {'problem'} -> nsigmas;
 291
 292     open(INFN, ">infn.f");
 293
 294   print INFN << "EOF";
 295       SUBROUTINE INFN(ICALL,THETA,DATREC,INDXS,NEWIND)
 296       DIMENSION THETA(*),DATREC(*),INDXS(*)
 297       DOUBLE PRECISION THETA
 298       COMMON /ROCM6/ THETAF(40),OMEGAF(30,30),SIGMAF(30,30)
 299       COMMON /ROCM7/ SETH(40),SEOM(30,30),SESIG(30,30)
 300       COMMON /ROCM8/ OBJECT
 301       COMMON /ROCM9/ IERE,IERC
 302       DOUBLE PRECISION THETAF, OMEGAF, SIGMAF
 303       DOUBLE PRECISION OBJECT
 304       REAL SETH,SEOM,SESIG
 305       DOUBLE PRECISION ETA(10)
 306       INTEGER J,I
 307       INTEGER IERE,IERC
 308       INTEGER MODE
 309       INTEGER NTH,NETA,NEPS
 310 EOF
 311
 312     print INFN "      DATA NTH,NETA,NEPS/$ntheta,$neta,$neps/\n";
 313
 314     print INFN "      IF (ICALL.EQ.0) THEN\nC     open files here, if necessary\n";
 315     print INFN "      OPEN(50,FILE='cwtab".$self -> sdno().$mirror_name.".est')\n";
 316
 317 #    for( my $i = 0; $i <= $#cwtab_names; $i++ ) {
 318 #      print INFN '         OPEN(5'.$i.',FILE=\''.$cwtab_names[$i].'\')'."\n";
 319 #    }
 320
 321   print INFN << "EOF";
 322       ENDIF
 323       IF (ICALL.EQ.3) THEN
 324          MODE=0
 325          CALL PASS(MODE)
 326          MODE=1
 327          WRITE (50,*) 'ETAS'
 328  20      CALL PASS(MODE)
 329          IF (MODE.EQ.0) GO TO 30
 330          IF (NEWIND.NE.2) THEN
 331             CALL GETETA(ETA)
 332             WRITE (50,97) (ETA(I),I=1,NETA)
 333          ENDIF
 334          GO TO 20
 335  30      CONTINUE
 336          WRITE (50,*) 'THETAS'
 337          WRITE (50,99) (THETAF(J),J=1,NTH)
 338          WRITE (50,*) 'OMEGAS'
 339          DO 7000 I=1,NETA
 340  7000       WRITE (50,99) (OMEGAF(I,J),J=1,NETA)
 341             WRITE (50,*) 'SIGMAS'
 342             DO 7999 I=1,NEPS
 343  7999          WRITE (50,99) (SIGMAF(I,J),J=1,NEPS)
 344             ENDIF
 345  99         FORMAT (20E15.7)
 346  98         FORMAT (2I8)
 347  97         FORMAT (10E15.7)
 348             RETURN
 349             END
 350 EOF
 351
 352     close INFN;
 353   }
 354 }
 355 end post_process
 356
 357 # }}}