lib/model/iofv_module_subs.pm

   1 # {{{ Original Documentation
   2
   3 ###############################################################
   4 #                            INTRODUCTION                     #
   5 ###############################################################
   6
   7 # Individual objective function values can be a useful model
   8 # selection diagnostic. Their use in the selection of covariates
   9 # for inclusion in a non-linear mixed effects model is descibed
  10 # in: Sadray, Jonsson and Karlsson, PharmRes, 16(8) 1999,
  11 # pp. 1260-1265
  12
  13 # This script calculates the individual difference in objective
  14 # function value between two models ("basic" and "full").  Optionally
  15 # it also creates a file with the necessary fortran sub-routine.
  16
  17 # To use the script, you need perl installed on your system. Perl
  18 # can be obtained free of charge from www.perl.com. Perl, and this
  19 # script, can be used on both UNIX machines and PCs.
  20
  21 # Questions, comments and bug reports should be sent to
  22 # niclas.jonsson@biof.uu.se.
  23
  24 # (c) Copyright 1999 Niclas Jonsson and Mats Karlsson
  25
  26
  27 ################################################################
  28 #                  TO CREATE THE CONTR SUBROUTINE              #
  29 ################################################################
  30
  31 # To create the CONTR subroutine that is necessary to extract
  32 # the individual objective function values, do the following:
  33
  34 #  1. On the command line, type:
  35 #     perl ofv1 -p
  36
  37 #     This creates a file called iofvcont.f.
  38
  39 #  2. Open the iofvcont.f file for editing and change the
  40 #     parameter statement on line two to reflect the number of
  41 #     observations you have NONMEM compiled for. (The NONMEM
  42 #     default is 50.)
  43
  44 #  3. Open the NM-TRAN model file for which you want to obtain
  45 #     individual objective function values and add the following line
  46 #     before the $SUBROUTINE line:
  47 #     $CONTR DATA=(ID)
  48
  49 #  4. Change the $SUBROUTINE line to include the iofvcont.f file, e.g.:
  50 #     $SUBROUTINE ADVAN4 TRANS2 CONTR=iofvcont.f
  51
  52 #  5. When the model is run, a file called fort80 (or something similar)
  53 #     is created. It containes two columns. The first is the ID numbers
  54 #     and the second the individual objectives.
  55 #     There are (no of individuals)*(number of function evaluations) lines
  56 #     the file. Only the last (no of individuals) lines are of interest.
  57
  58 #########################################################################
  59 # TO COMPUTE THE INDIVIDUAL DIFFERENCES IN THE OBJECTIVE FUNCTION VALUE #
  60 #                          BETWEEN TWO MODELS                           #
  61 #########################################################################
  62
  63 # BACKGROUND
  64
  65 # The script assumes that the ID and associated individual objective
  66 # function value are stored in two files.  These files will
  67 # automatically be the output from the "iofvcont.f" file although
  68 # the user has to rename the output file after each run.
  69
  70 # This is what you have to do in order to get a list of individual
  71 # objective value differences between two models:
  72
  73 # 1. Create the iofvcont.f file (see above):
  74
  75 # 2. For each model you run, a file called "fort80" (or similar,
  76 #    depending on your compiler) will be created.  It will usually be
  77 #    very long and only the last number for each individual will be
  78 #    used, but that is taken care of by the script.  What you need to
  79 #    do is to rename this file after every run (to assure that it is
  80 #    not written over).  If basic and final models are runs 5 and 6,
  81 #    the following code can serve as example:
  82
  83 #    nmfe run5.mod run5.lst
  84 #    mv ftn80 iotab5
  85 #    nmfe run6.mod run6.lst
  86 #    mv ftn80 iotab6
  87
  88 # 3. When you have run the basic and full model run the perl script by
  89 #    typing:
  90 #    perl ofv1
  91
  92 #    Thereafter follows instructions. Names of output files will be given.
  93
  94 # }}}
  95
  96 # {{{ new
  97
  98 start new
  99 {
 100
 101   my $base_model = $this -> {'base_model'};
 102
 103   if( $base_model -> is_option_set( record => 'subroutine',
 104                                     name => 'CONTR' ) ){
 105
 106     'debug' -> die( message => 'CONTR in $SUBROUTINE is already set, iofv cannot be computed' );
 107
 108   }
 109
 110   $base_model -> set_option( record_name => 'subroutine',
 111                              option_name => 'CONTR',
 112                              option_value => 'iofvcont.f' );
 113
 114   $base_model -> set_records( type => 'contr',
 115                               record_strings => ['DATA=(ID)'] );
 116
 117
 118 }
 119 end new
 120
 121 # }}}
 122
 123 # {{{ post_process
 124
 125 start post_process
 126 {
 127
 128   unless( defined $self -> {'nm_version'} ){
 129     $self -> {'nm_version'} = 'default';
 130   }
 131
 132   my ($sizes_dir,$junk) = split(/,/ , $PsN::config -> {'nm_versions'} -> { $self -> {'nm_version'} } );
 133
 134   my $base_model = $self -> {'base_model'};
 135
 136   # Figure out if we have an sdtab and what number it has
 137   my ( $sd_ref, $junk ) = $base_model -> problems -> [0] ->
 138       _option_val_pos( name        => 'FILE',
 139                        record_name => 'table',
 140                        exact_match => 0 );
 141
 142   my $sdno = '1';
 143
 144   if( defined $sd_ref ) {
 145     foreach my $tabname ( @{$sd_ref} ) {
 146       if( $tabname =~ /[sd,pa]tab(\d+)/ ) {
 147         $sdno= $1;
 148       }
 149     }
 150   }
 151
 152   open(CONTR,">iofvcont.f");
 153   print CONTR << "EOF";
 154       subroutine contr (icall,cnt,ier1,ier2)
 155       INCLUDE "$sizes_dir/SIZES"
 156 C     parameter (no=50)
 157       common /rocm1/ y(no),data(no,3),nobs
 158       integer nobs, un
 159       double precision cnt,y
 160       OPEN(80,FILE=\'iotab$sdno\')
 161       if (icall.le.1) return
 162       call ncontr (cnt,ier1,ier2,l2r)
 163 C     individual obj. funct. value for indiv. jj = cnt
 164       write(80,10) data(1,1),cnt
 165    10 FORMAT(1E12.4E2,1E12.4E2)
 166       return
 167       end
 168
 169 EOF
 170
 171     close CONTR;
 172
 173 }
 174 end post_process
 175
 176 # }}}
 177
 178 # {{{ post_run_process
 179 start post_run_process
 180 {
 181
 182   my $base_model = $self -> {'base_model'};
 183
 184   # Figure out if we have an sdtab and what number it has
 185   my ( $sd_ref, $junk ) = $base_model -> problems -> [0] ->
 186       _option_val_pos( name        => 'FILE',
 187                        record_name => 'table',
 188                        exact_match => 0 );
 189
 190   my $sdno = '1';
 191
 192   if( defined $sd_ref ) {
 193     foreach my $tabname ( @{$sd_ref} ) {
 194       if( $tabname =~ /[sd,pa]tab(\d+)/ ) {
 195         $sdno= $1;
 196       }
 197     }
 198   }
 199
 200   my $data = $base_model -> datas -> [0];
 201
 202   my $ids = $data -> column_to_array( column => $data -> idcolumn -1 );
 203
 204   my ($first_id, $last_id) = ($ids -> [0],$ids -> [$#{$ids}]);
 205
 206   if( -e "iotab$sdno" ){
 207     open( IOTAB, "<iotab$sdno" );
 208     my @iotab = <IOTAB>;
 209     my @values;
 210     my @starting_points;
 211
 212     my $previous_id = $last_id;
 213
 214     for( my $i = 0;$i < scalar @iotab; $i++ ){
 215       my @line = split( ' ',$iotab[$i] );
 216
 217       push( @values, $line[1] );
 218
 219       if( $previous_id == $last_id and $line[0] == $first_id ){
 220         push( @starting_points, $i );
 221       }
 222
 223       $previous_id = $line[0];
 224     }
 225
 226     if( defined $base_model -> covariance ){
 227
 228       @iotab = @iotab[$starting_points[$#starting_points-1] .. $starting_points[$#starting_points]-1];
 229
 230     } else {
 231
 232       @iotab = @iotab[$starting_points[$#starting_points-1] .. $starting_points[$#starting_points]-1];
 233
 234     }
 235     close(IOTAB);
 236
 237     open IOTAB, ">iotab$sdno";
 238     print( IOTAB @iotab );
 239     close IOTAB;
 240
 241     my @eo;
 242     if ( defined $base_model -> extra_output() ) {
 243       @eo = @{$base_model -> extra_output()};
 244     }
 245
 246     push( @eo, "iotab$sdno" );
 247     $base_model -> extra_output( \@eo );
 248
 249   } else {
 250     'debug' -> die( message => "Unable to open iotab: iotab$sdno ." );
 251   }
 252 }
 253 end post_run_process
 254 # }}}
 255
 256 print("Enter the name of the file with the OFVs for the basic model:\n");
 257 $basic = <STDIN>;
 258
 259 print("Enter the name of the file with the OFVs for the full model:\n");
 260 $full = <STDIN>;
 261
 262 print("Enter the number of subjects:\n");
 263 $ind = <STDIN> ;
 264
 265 chomp($basic);
 266 chomp($full);
 267
 268 $out1 = "deltofv1.res";
 269
 270 ## Open the files with OFVs and extract the relevant values
 271 open (FILE1, $basic) || die("Couldn't open $basic for reading!\n");
 272 while(<FILE1>) {
 273   chomp;
 274   ($id,$basic_mof) = split(' ',$_);
 275   push @basic, $basic_mof;
 276   push @id, $id;
 277 }
 278 @basic = splice(@basic,-$ind);
 279 @ids   = splice(@id,-$ind);
 280 close FILE1;
 281
 282 open (FILE2, $full) || die("Couldn't open $full for reading!\n");
 283 while(<FILE2>) {
 284   chomp;
 285   ($junk,$full_mof) = split(' ',$_);
 286   push @full, $full_mof;
 287 }
 288 @full = splice(@full,-$ind);
 289 close FILE2;
 290
 291 ## Compute the differences
 292
 293 foreach $i (0..@basic-1) {
 294   push @diff, $basic[$i] - $full[$i];
 295 }
 296
 297 ## Create the arrays to be printed
 298 foreach $diff (@diff) {
 299   push @unsorted, sprintf "%12.4e%12.4e", $ids[$index],$diff[$index];
 300   $index++;
 301 }
 302
 303 @sorted =
 304   map  {$_-> [0]}
 305   sort {$a -> [1] <=> $b ->[1]}
 306   map  {[$_, (split ' ')[1]]}
 307 @unsorted;
 308
 309 open (OUTFILE1, ">$out1") ||
 310   die ("Cannot open output file $out1 for writing!\n");
 311
 312 print OUTFILE1 ("ID DELTAOFV SORTID SORTDELTAOFV\n");
 313 $index = 0;
 314 foreach $diff (@diff) {
 315   printf OUTFILE1 "%s    %s\n", $unsorted[$index],$sorted[$index];
 316   $index++;
 317 }
 318 close OUTFILE1;
 319 print ("Output are in file ",$out1,"\n");
 320