#! /bin/tcsh -f
#
#	test of French-Wilson "truncate" procedure in different programs
#
# script requires:
# phenix
# CCP4 suite
# xdsconv
#

set CELL = ( 20 20 20 )
set grid = ( 256 256 256 )

# aim for a particular Matthews volume
set N    = `echo $CELL | awk '{print int($1*$2*$3/14/2.4)}'`
set B    = 20

# pick some random seeds
set seed = `awk 'BEGIN{srand();for(i=1;i<10;++i) print int(sqrt(rand()**2)*10000)}'`

echo "making a random pdb file of $N atoms with B=$B and cell: $CELL 90 90 90"
echo "$CELL $N $B $seed[1]" |\
awk '{a=$1;b=$2;c=$3;N=$4;B=$5;srand($6);\
	for(i=1;i<=N;++i){\
	    x=a*rand();\
	    y=b*rand();\
	    z=c*rand();\
	    print x,y,z,B}\
     }' |\
awk '{printf"ATOM %6d  %-3s %3s %1s%4d     %7.3f %7.3f %7.3f %5.2f%6.2f\n",\
          1, "OW1", "WAT", " ", 1, $1, $2, $3, 1, $4}' |\
cat >! temp.pdb
pdbset xyzin temp.pdb xyzout random.pdb << EOF > /dev/null
CELL $CELL 90 90 90
SPACE 1
EOF

# compute the structure factors
echo "calculating structure factors to 0.5A with phenix.fmodel"
rm -f fmodel.mtz
phenix.fmodel random.pdb k_sol=0 high_res=0.5 file_name=fmodel.mtz \
 algorithm=direct > /dev/null

cp fmodel.mtz fmodel_wilsonme.mtz

# add some noise using two sets so we can compute CC1/2
echo "adding some noise"
rm -f fmodel_noise.mtz
sftools << EOF >! sftools.log
read fmodel.mtz
calc SEED $seed[2]
calc F col noise1 = RAN_G 1.4142135623731 /
calc SEED $seed[3]
calc F col noise2 = RAN_G 1.4142135623731 /
calc J col I = col FMODEL col FMODEL *
calc J col Ihalf = col I 2 /
calc J col Inoisy1half = col Ihalf col noise1 +
calc J col Inoisy2half = col Ihalf col noise2 +
calc J col IMEAN = col Inoisy1half col Inoisy2half +
calc Q col SIGIMEAN = 1.0
write fmodel_noise.mtz 
exit
y
EOF
echo "truncate"
truncate hklin fmodel_noise.mtz hklout truncated_wrongscale.mtz << EOF >! truncate.log
resolution 0.5
EOF
# fix the wrong scale from truncate
cad hklin1 fmodel_noise.mtz hklin2 truncated_wrongscale.mtz hklout scaleme.mtz << EOF > /dev/null
labin file 1 E1=FMODEL E2=PHIFMODEL E3=SIGIMEAN
labin file 2 E1=F E2=SIGF
EOF
echo "scaleit"
scaleit hklin scaleme.mtz hklout scaled.mtz << EOF >! scaleit.log
refine aniso
labin FP=FMODEL SIGFP=SIGIMEAN FPH1=F SIGFPH1=SIGF
EOF
egrep "Sc_kraut|TOTALS" scaleit.log
set scale = `awk '$1=="Derivative" && ! /itle/{print $3}' scaleit.log | tail -1`
set B     = `awk '/equivalent iso/{print $NF}' scaleit.log | tail -1`
echo "scale= $scale B= $B"
cad hklin1 scaled.mtz hklout truncated_wilsonme.mtz << EOF > /dev/null
labin file 1 E1=F E2=SIGF
EOF



echo "ctruncate"
echo "" |\
ctruncate -mtzin fmodel_noise.mtz -mtzout ctruncated_wrongscale.mtz \
  -colin '/*/*/[IMEAN,SIGIMEAN]' >! ctruncate.log
# fix the wrong scale from ctruncate
cad hklin1 fmodel_noise.mtz hklin2 ctruncated_wrongscale.mtz hklout scaleme.mtz << EOF > /dev/null
labin file 1 E1=FMODEL E2=PHIFMODEL E3=SIGIMEAN
labin file 2 E1=F E2=SIGF
EOF
echo "scaleit"
scaleit hklin scaleme.mtz hklout scaled.mtz << EOF >! scaleit.log
refine aniso
labin FP=FMODEL SIGFP=SIGIMEAN FPH1=F SIGFPH1=SIGF
EOF
egrep "Sc_kraut|TOTALS" scaleit.log
set scale = `awk '$1=="Derivative" && ! /itle/{print $3}' scaleit.log | tail -1`
set B     = `awk '/equivalent iso/{print $NF}' scaleit.log | tail -1`
echo "scale= $scale B= $B"
cad hklin1 scaled.mtz hklout ctruncated_wilsonme.mtz << EOF > /dev/null
labin file 1 E1=F E2=SIGF
EOF




echo "xdsconv"
cat << EOF >! forxds.ahkl
!FORMAT=XDS_ASCII    MERGE=TRUE     FRIEDEL'S_LAW=FALSE
!DATE=28-Mar-2013
!GENERATED BY: XSCALE (VERSION  July 4, 2012)          
!OUTPUT_FILE=forxds.ahkl                                       
!COMPRISES THE FOLLOWING SCALED INPUT FILES:
! ISET=   1 INPUT_FILE=XDS_ASCII.HKL                                     
! ISET=   1 X-RAY_WAVELENGTH=   0 (<0 if unknown)
!
!SPACE_GROUP_NUMBER=  1
!UNIT_CELL_CONSTANTS=     20.00     20.00     20.00  90.000  90.000  90.000
!NUMBER_OF_ITEMS_IN_EACH_DATA_RECORD=   5
!ITEM_H=1
!ITEM_K=2
!ITEM_L=3
!ITEM_IOBS=4
!ITEM_SIGMA(IOBS)=5
!END_OF_HEADER
EOF
cad hklin1 fmodel_noise.mtz hklout dumpme.mtz << EOF > /dev/null
labin file 1 E1=IMEAN E2=SIGIMEAN
#resolution over_all 1.0
EOF
echo "FORMAT '(3i6,2g40.20)'\nnref -1" | mtzdump hklin dumpme.mtz |\
awk 'substr($0,1,18)==sprintf("%6d%6d%6d",$1,$2,$3){\
   printf("%6d%6d%6d%11.3E%11.3E\n",$1,$2,$3,$4,$5)}' |\
cat >> forxds.ahkl
echo '\!END_OF_DATA' >> forxds.ahkl

cat << EOF >! XDSCONV.INP
OUTPUT_FILE=temp.hkl CCP4_F
INPUT_FILE=forxds.ahkl
WILSON_STATISTICS=TRUE
EOF
xdsconv >! xdsconv.log
awk '{gsub("FP","F"); print}' F2MTZ.INP | f2mtz HKLOUT xdsconv.mtz > /dev/null
cad hklin1 xdsconv.mtz hklout xdsconv_wilsonme.mtz << EOF > /dev/null
labin file 1 E1=F E2=SIGF
EOF



echo "making Wilson plots..."
foreach file ( fmodel truncated ctruncated xdsconv )
echo "$file"
echo "nref -1\nlreso\nFORMAT '(3i6,3g35.20)'" |\
 mtzdump hklin ${file}_wilsonme.mtz |\
awk 'substr($0,1,18)==sprintf("%6d%6d%6d",$1,$2,$3) && $5+0 !~ -999{print $4/4,$5,$6}' |\
awk -v bs=0.005 '{bin=sprintf("%.0f",$1/bs);sum1[bin]+=$2*$2;sum2[bin]+=$3;++count[bin]}\
     END{for(bin in sum1) print bin*bs,sum1[bin]/count[bin],sum2[bin]/count[bin],count[bin]}' |\
sort -g |\
tee ${file}_wilson.txt | wc
end


echo "plotting..."
gnuplot -persist << EOF
set log y
set title "Wilson plot"
set xlabel "(sin(theta)/lambda)^2"
set ylabel "average F^2"
set grid
plot [:0.5] [1e-5:1e5] 'fmodel_wilson.txt' ti "true F" with l,\
  'xdsconv_wilson.txt' ti "xdsconv" with l,\
  'truncated_wilson.txt' ti "truncate" with l, \
  'ctruncated_wilson.txt' ti "ctruncate" lt 6
set terminal png
set output "truncated_wilsons.png"
replot
set terminal x11
EOF






display truncated_wilsons.png