#!/usr/bin/perl -w
if(!@ARGV)
{ use File::Basename;
  print "
Usage:  ${\basename($0)}  <required arguments>  <options>

Required arguments:

-hklin   filename.mtz   - Input MTZ file (output from a refinement).
-pdbin   filename.pdb   - Input PDB file.
-- or --
-cifin   filename.mmcif - Input mmCIF file.

NOTE that the PDB or mmCIF file MUST have been output in the SAME
refinement run that produced the MTZ file, or incorrect results will be
obtained!

Optional arguments:

-atomsf  <filename>     - Atomic scattering factor file (default = get
                          experiment type from the PDB/mmCIF header).
-flabel  <string>       - Column label for Fobs (by default tries F, FP,
                          FOSC, F_xxx, FP_xxx where xxx = any string).
-grid    <integer>      - Grid sampling (default = 4.5; must be >= 4).
-byatom                 - Use by-atom averages (default is by-residue).
-select  <string>       - Residue selection (see program documentation
                          for details, except enclosing quotes are not
                          needed here unless the string contains
                          spaces).
-bfac    <string>       - Statistic to place in Biso column of output
                          PDB file: R, RG, SRG, CCS, CCP, ZCCP, ZD-,
                          ZD+, ZD, ZO or NONE (default = ZO).
-occu    <string>       - Statistic to place in occupancy column of
                          output PDB file (options as for -bfac except
                          default = ZD-).
-null    <string>       - String to indicate null entries in residue
                          table (max. 4 chars; default = 'n/a').
-usealt                 - Flag to indicate atoms are to be grouped
                          according to the alternate location indicator.
-useh                   - Flag to indicate H atoms are to be included.
-log     <filename>     - Log file (default is standard output).
-output  <filename>     - Output data file (default = 'edstats.out').
-hisout  <filename>     - Output density histogram plot (default =
                          none).
-ppdout  <filename>     - Output P-P difference plot (default = none).
-qqdout  <filename>     - Output Q-Q difference plot (default = none).
-mapout1 <filename>     - Output rescaled rho_obs map (default = none).
-mapout2 <filename>     - Output rescaled rho_diff map (default = none).
-pdbout  <filename>     - Output PDB file with atom scores in occupancy
                          and Biso columns (default = none).
-cifout  <filename>     - Output mmCIF file with atom Z scores in
                          occupancy & Biso columns (default = none).
-tmpdir  <dir>          - Directory for temporary files
                          (default = '/tmp/\$USER/').
-exec    <path>         - Path for executables (default = none).
-test    <integer>      - Test level (default = none).
-noerror                - Use input MTZ file after error in MTZFIX (NOT
                          RECOMMENDED!).
-debug                  - Debug (don't delete temp dir).

The columns and column labels in the output file are:

 1. RT:   Residue type (3-letter code).
 2. CI:   2 characters: chain ID (reset to underscore if blank) and
          alternate location indicator (if -usealt is specified).
 3. RN:   Residue name (including insertion code if present).

Then for the main-chain atoms in the residue (including CB):

 4. BAm:   Weighted average Biso.
 5. NPm:   No of statistically independent grid points covered by atoms.
 6. Rm:    Real-space R factor (RSR).
 7. RGm:   Real-space RG factor (RSRG).
 8. SRGm:  Standard uncertainty of RSRG.
 9. CCSm:  Real-space sample correlation coefficient (RSCC).
10. CCPm:  Real-space 'population' correlation coefficient (RSPCC).
11. ZCCPm: Z-score of real-space correlation coefficient.
12. ZOm:   Real-space Zobs score (RSZO).
13. ZDm:   Real-space Zdiff score (RSZD) i.e. max(-RSZD-,RSZD+).
14. ZD-m:  Real-space Zdiff score for negative differences (RSZD-).
15. ZD+m:  Real-space Zdiff score for positive differences (RSZD+).

Columns 16-27 (with labels 'BAs' etc.) contain the same information as
columns 4-15 above (i.e. add 12) but for the side-chain atoms
(excluding CB) if present.  Columns 28-39 (with labels 'BAa' etc.
contain the corresponding all-atom (main + side-chain) statistics.

40. MN:    Model number.
41. CP:    PDB chain ID.
42. NR:    PDB residue no.

";
  exit;
}

# Default settings.
$opt_grid = 4.5;
$opt_output = 'edstats.out';
$opt_usealt = $opt_exec = $opt_noerror = '';

# GetOptions zaps the argument list, so make a copy.
@args = @ARGV;

use Getopt::Long;

GetOptions('hklin=s','pdbin=s','xyzin=s','cifin=s','atomsf=s',
  'flabel=s','grid=i','byatom','select=s','bfac=s','occu=s','null=s',
  'usealt','useh','log=s','output=s','hisout:s','ppdout:s','qqdout:s',
  'mapout1:s','mapout2:s','pdbout:s','xyzout:s','cifout:s','tmpdir=s',
  'exec=s','test:i','noerror','debug') or
  die "$/ERROR in processing command line arguments,";

# If debugging echo command line for reassurance.
if($opt_debug)
{ $command = $0;
  $command .= /^$|[\s\#]/?" '$_'":" $_" for @args;
  print "\nCommand line: $command\n\n";
}

# Check we have CCP4 and all input files.
$ENV{CCP4} or die "$/ERROR: you must set up the CCP4 environment,";

$opt_hklin or die "$/ERROR: MTZ file not specified,";

if(!-e $opt_hklin)
{ -e "$opt_hklin.mtz" or die
    "$/ERROR: input MTZ file $opt_hklin not found,";
  $opt_hklin .= '.mtz';
}

$n=0;
++$n if $opt_pdbin;
$n+=2 if $opt_xyzin;
$n+=4 if $opt_cifin;

if($n && $n<=2)
{ $opt_pdbin = $opt_xyzin if $n==2;
  if(!-e $opt_pdbin)
  { -e "$opt_pdbin.pdb" or die
      "$/ERROR: input PDB file $opt_pdbin not found,";
    $opt_pdbin .= '.pdb';
  }

  if(!$opt_atomsf)
  { open IN,$opt_pdbin or die "$/ERROR opening $opt_pdbin for read,";
    while(<IN>)
    { if(/^EXPDTA +(.+)/)
      { $exp = $1;
	last;
      }
    }

    close IN;
  }

  $in="PDBIN $opt_pdbin";
}
elsif($n==4)
{ if(!-e $opt_cifin)
  { -e "$opt_cifin.mmcif" or die
      "$/ERROR: input mmCIF file $opt_cifin not found,";
    $opt_cifin .= '.mmcif';
  }

  if(!$opt_atomsf)
  { open IN,$opt_cifin or die "$/ERROR opening $opt_cifin for read,";
    while(<IN>)
    { if(/^_exptl.method +['"](.+)['"]/)
      { $exp = $1;
	last;
      }
    }

    close IN;
  }

  $opt_cifin="./$opt_cifin" if $opt_cifin!~/\// && $opt_cifin=~/\.cif$/;
  $in = "CIFIN $opt_cifin";
}
else {die "$/ERROR: One of PDB or mmCIF file must be specified,"}

$opt_pdbout = 'edstats.pdb' if defined $opt_pdbout && !$opt_pdbout;
$opt_xyzout = 'edstats.pdb' if defined $opt_xyzout && !$opt_xyzout;
$opt_cifout = 'edstats.mmcif' if defined $opt_cifout && !$opt_cifout;

$n=0;
++$n if $opt_pdbout;
$n+=2 if $opt_xyzout;
$n+=4 if $opt_cifout;

if(!$n) {$out=''}
elsif($n<=2)
{ $opt_pdbin or die "$/ERROR: PDB input file required,";
  $opt_pdbout = $opt_xyzout if $n==2;
  $opt_pdbout=~/\./ or $opt_pdbout .= '.pdb';
  $out=" PDBOUT $opt_pdbout";
}
elsif($n==4)
{ $opt_cifin or die "$/ERROR: mmCIF input file required,";
  $opt_cifout="./$opt_cifout" if $opt_cifout!~/\// &&
    $opt_cifout=~/\.cif$/;
  $opt_cifout=~/\./ or $opt_cifout .= '.mmcif';
  $out=" CIFOUT $opt_cifout";
}
else {die "$/ERROR: Only one of PDBOUT or CIFOUT may be specified,"}

if(!$opt_atomsf)
{ if(!$exp) {warn "$/WARNING: Experiment type is not defined, assuming",
    " 'X-RAY DIFFRACTION'$/"}
  elsif($exp=~/^NEUTRON/i) {$opt_atomsf='atomsf_neutron.lib'}
  elsif($exp=~/^ELECTRON/i) {$opt_atomsf='atomsf_electron.lib'}
  elsif($exp!~/^X-RAY/i) {die "$/ERROR: Invalid experiment type '$exp',"}
}

# $opt_exec will be the full path to the EDSTATS executable (filename =
# 'edstats'); $path is the directory part for MTZINFO & MTZFIX only.
if(!($path = $opt_exec)) {$opt_exec = 'edstats'}
elsif(-d $path) {$opt_exec =~ s=/?$=/edstats=}
else {$path =~ s=[^/]+$==}

# Make sure executables can be found somewhere.
use Env '@PATH';

$mtzinfo='mtzinfo';
grep(-x "$_/mtzinfo",@PATH) || -x ($mtzinfo="${path}mtzinfo") or
  die "$/ERROR: mtzinfo executable '${path}mtzinfo' not found in path,";

$mtzfix='mtzfix';
grep(-x "$_/mtzfix",@PATH) || -x ($mtzfix="${path}mtzfix") or
  die "$/ERROR: mtzfix executable '${path}mtzfix' not found in path,";

grep(-x "$_/$opt_exec",@PATH) || -x $opt_exec or
  die "$/ERROR: edstats executable '$opt_exec' not found in path,";

# Options for MTZFIX.
$opt_flabel = $opt_flabel?" FLABEL $opt_flabel":'';

# Options for FFT.
$opt_grid=~/^\d+(\.\d*)?$/ && $opt_grid>=4 or
  die "$/ERROR: grid sampling must be >= 4,";

# Options for EDSTATS.
$opt_byatom = $opt_byatom?',MAIN=atom,SIDE=atom':'';
$opt_select = $opt_select?",SELE=\\'$opt_select\\'":'';
$opt_bfac = $opt_bfac?",BFAC=$opt_bfac":'';
$opt_occu = $opt_occu?",OCCU=$opt_occu":'';

$opt_hisout = 'edstats.his' if defined $opt_hisout && !$opt_hisout;
$opt_hisout = $opt_hisout?" HISOUT $opt_hisout":'';

$opt_ppdout = 'edstats.ppd' if defined $opt_ppdout && !$opt_ppdout;
$opt_ppdout = $opt_ppdout?" PPDOUT $opt_ppdout":'';

$opt_qqdout = 'edstats.qqd' if defined $opt_qqdout && !$opt_qqdout;
$opt_qqdout = $opt_qqdout?" QQDOUT $opt_qqdout":'';

$opt_mapout1 = 'edstats-fo.map' if defined $opt_mapout1 &&
  !$opt_mapout1;
$opt_mapout1 = $opt_mapout1?" MAPOUT1 $opt_mapout1":'';

$opt_mapout2 = 'edstats-df.map' if defined $opt_mapout2 &&
  !$opt_mapout2;
$opt_mapout2 = $opt_mapout2?" MAPOUT2 $opt_mapout2":'';

$opt_log = $opt_log?">$opt_log":'';

# Create temporary files directory if it doesn't exist.
if($opt_tmpdir)
{ if(-d $opt_tmpdir) {$keep = 1}
  else {mkdir $opt_tmpdir,0755 or die "$/ERROR in mkdir $opt_tmpdir,"}
}
else
{ $opt_tmpdir = "/tmp/$ENV{USER}";
  mkdir $opt_tmpdir,0755 or die "$/ERROR in mkdir $opt_tmpdir,"
    if !-e $opt_tmpdir;
  use File::Temp qw(tempdir);
  $opt_tmpdir = tempdir('edstats-XXXXXX',DIR=>$opt_tmpdir);
}

# Get column labels & resolution limits from the MTZ header.
open IN,"$mtzinfo $opt_hklin |";

while(<IN>)
{ print;
  if(/^LABELS/i) {@a = split}
  elsif(/^XDATA/i) {($rl,$rh) = (split)[7,8]}
}
close IN;

# Set up column labels for FFT.
if(grep(/2FOFCWT/,@a) && grep(/PH2FOFCWT/,@a) && grep(/FOFCWT/,@a) &&
  grep(/PHFOFCWT/,@a))
  {($fo,$po,$fd,$pd) = ('2FOFCWT','PH2FOFCWT','FOFCWT','PHFOFCWT')}

elsif(grep(/FWT/,@a) && grep(/PHWT/,@a) && grep(/DELFWT/,@a) &&
  grep(/PHDELWT/,@a))
  {($fo,$po,$fd,$pd) = ('FWT','PHWT','DELFWT','PHDELWT')}

else {die "$/ERROR - standard map coefficient labels not found,"}

# Fix up the map coefficients: FLABEL specifies the label for Fobs &
# sigma(Fobs) (defaults are F/SIGF or FOSC/SIGFOSC or FP/SIGFP or
# F_xxx/SIGF_xxx or FP_xxx/SIGFP_xxx).  Here, 'in.mtz' is the output
# reflexion file from the refinement program in MTZ format.
$cmd="rm -f $opt_tmpdir/mtzfix.mtz; $mtzfix$opt_flabel".
  " HKLIN $opt_hklin HKLOUT $opt_tmpdir/mtzfix.mtz".
  " >$opt_tmpdir/mtzfix.log";

print "$/$cmd$/" if $opt_debug;
!system $cmd or $opt_noerror or die "$/ERROR in mtzfix,";
  
# If no fix-up was needed, or error in MTZFIX, use the original file.
if(!-e "$opt_tmpdir/mtzfix.mtz")
{ if($opt_hklin!~/^\//)
  { use Cwd;
    $opt_hklin = &cwd.'/'.$opt_hklin;
  }

  print "symlink $opt_hklin,\"$opt_tmpdir/mtzfix.mtz\"$/" if $opt_debug;
  symlink $opt_hklin,"$opt_tmpdir/mtzfix.mtz";
}

# Compute the 2mFo-DFc map (mFo for centrics).  Note that EDSTATS needs
# only 1 asymmetric unit (but will also work with more).  Grid sampling
# must be at least 4.
$cmd="echo 'LABI F1=$fo PHI=$po\nXYZL asu\nGRID samp $opt_grid'".
  " |fft HKLIN $opt_tmpdir/mtzfix.mtz MAPOUT $opt_tmpdir/fo.map".
  " >$opt_tmpdir/fft-fo.log";

print "$/$cmd$/" if $opt_debug;
!system $cmd or die "$/ERROR in fft(1),";

# Compute the 2(mFo-DFc) map (mFo-DFc for centrics).
$cmd="echo 'LABI F1=$fd PHI=$pd\nXYZL asu\nGRID samp $opt_grid'".
  " |fft HKLIN $opt_tmpdir/mtzfix.mtz MAPOUT $opt_tmpdir/df.map".
  " >$opt_tmpdir/fft-df.log";

print "$/$cmd$/" if $opt_debug;
!system $cmd or die "$/ERROR in fft(2),";

# ATOMSF option for EDSTATS.
$opt_atomsf = $opt_atomsf?" ATOMSF $opt_atomsf":'';

# NULL option for EDSTATS.
$opt_null = $opt_null?",NULL=$opt_null":'';

# USEALT option for EDSTATS.
$opt_usealt = $opt_usealt?",USEA=t":'';

# USEH option for EDSTATS.
$opt_useh = $opt_useh?',USEH=t':'';

# Test option for EDSTATS.
$opt_test = -1 if defined $opt_test && !$opt_test;
$opt_test = $opt_test?",TEST=$opt_test":'';

# Main- & side-chain residue statistics.
$cmd="echo RESL=$rl,RESH=$rh$opt_byatom$opt_select$opt_bfac$opt_occu".
  "$opt_null$opt_usealt$opt_useh$opt_test |$opt_exec".
  " MAPIN1 $opt_tmpdir/fo MAPIN2 $opt_tmpdir/df $in$opt_atomsf".
  " OUT $opt_output$opt_hisout$opt_ppdout$opt_qqdout$opt_mapout1".
  "$opt_mapout2$out$opt_log";

print "$/$cmd$/" if $opt_debug;
!system $cmd or die "$/ERROR in $opt_exec,";

# Delete temp directory if not required.
$opt_debug || $keep or system "rm -fr $opt_tmpdir";
