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