! Module file: refinementtarget
!
! CNS module
! **********
!
! Authors: Axel Brunger and Paul Adams
!
! copyright Yale University
!
! Function:
!    Defines the crystallographic refinement targets 
!
! Requirements:
!    Needs to be called within xray

module {refinementtarget}

(
  &target="residual";    {string}
  &sig_sigacv=0.07;      {real}
  &mbins=10;             {real}
  &fobs=fobs;            {reciprocal space array}
  &sigma=sigma;          {reciprocal space array}
  &weight=weight;        {reciprocal space array}
  &iobs=iobs;            {reciprocal space array}
  &sigi=sigi;            {reciprocal space array}
  &test=test;            {reciprocal space array}
  &fcalc=fcalc;          {reciprocal space array}
  &fpart=fpart;          {reciprocal space array}
  &pa=pa;                {reciprocal space array}
  &pb=pb;                {reciprocal space array}
  &pc=pc;                {reciprocal space array}
  &pd=pd;                {reciprocal space array}
  &phase=phase;          {reciprocal space array}
  &fom=fom;              {reciprocal space array}
  &sel=all;              {selection}
  &sel_test=none;        {selection}
  &statistics=false;     {logical}
)

set message ? end
evaluate ($message_old=$result)
set echo ? end
evaluate ($echo_old=$result)
if ( $log_level = verbose ) then
  set echo=on message=normal end
else
  set echo=off message=off end
end if

checkversion 1.3

if ( &target = "mlf" ) then

  declare name=r_test   domain=reciprocal type=real end
  declare name=r_eobs   domain=reciprocal type=real end
  declare name=r_ecalc  domain=reciprocal type=real end
  declare name=r_sigmaa domain=reciprocal type=real end
  query name=r_sigmad domain=reciprocal end
  if ( $object_exist = false ) then
    declare name=r_sigmad domain=reciprocal type=real end
  else
    do (r_sigmad=0) ( all )
  end if
  query name=r_dd domain=reciprocal end
  if ( $object_exist = false ) then
    declare name=r_dd     domain=reciprocal type=real end
  else
    do (r_dd=0) ( all ) 
  end if

  do (r_test=0) ( all )
  do (r_test=1) ( &sel_test )

  show sum (1) (&sel_test)
  if ( $result > 0 ) then
    evaluate ($test_exist=true)
    evaluate ($r_ntest=$result)
  else
    evaluate ($test_exist=false)
    display Error: no reflections flagged for cross-validation
    display        this target requires cross-validation
    abort
  end if

  {- Calculate number of bins used for E and r_sigmaa calculation.
     Make sure there are at least 50 reflections in each bin.
     Issue warning if there end up being less than 50 per bin   -}
  evaluate ($r_mbins=&mbins)
  if ( $r_mbins < 0 ) then
    show sum (1) (&sel_test)
    evaluate ($r_mbins=int($result/50))
    if ( $r_mbins <= 0 ) then
      display Error: there are less than 50 test set reflections
      display        please check the test set and the resolution limits
      abort
    end if
    if ( $r_mbins < 8 ) then
      display Warning: there are less than 50 reflections per bin
      display          you may want to increase the size of the test set
    end if
    evaluate ($r_mbins=max(8,$r_mbins))
    evaluate ($r_mbins=min($r_mbins,50))
  end if
  mbins=$r_mbins  

  {- Check to see that there are not too many bins with no
     test set reflections - limit is 1/5th of total bins -}
  declare name=r_temp domain=reciprocal type=real end
  do (r_temp=sum(r_test)) ( &sel )
  show sum(1) ( &sel )
  evaluate ($ref_per_bin=$select/$r_mbins)
  show sum(1) ( &sel and r_temp=0 )
  evaluate ($empty_bins=int($select/$ref_per_bin))
  undeclare name=r_temp domain=reciprocal end
  evaluate ($empty_ratio=$empty_bins/$r_mbins)
  if ( $empty_ratio >= 0.2 ) then
    display  Warning: too many bins with no test set reflections
    display           Is test set distributed randomly?
  end if

  {- compute r_sigmaa and r_sigmad -}
  if ( $r_ntest < 400 ) then
    do (r_eobs=norm(amplitude(&fobs)))                          ( &sel )
    do (r_ecalc=norm(amplitude(&fcalc+&fpart)))                 ( &sel )
  else
    do (r_eobs=norm(amplitude(&fobs)))                          ( &sel_test )
    do (r_ecalc=norm(amplitude(&fcalc+&fpart)))                 ( &sel_test )
  end if
  do (r_sigmaa=sigacv[sigma=&sig_sigacv](r_eobs,r_ecalc))     ( &sel_test )
  do (r_sigmad=sqrt(save(amplitude(&fobs)^2) (1-r_sigmaa^2))) ( &sel_test ) 
  {- compute d -}
  do (r_dd= r_sigmaa * sqrt(save(amplitude(&fobs)^2)  / 
                            save(amplitude(&fcalc+&fpart)^2))) ( &sel_test )
  {- Cross-validation active: need to expand to all selected
     reflections and compute the normalized structure factor
     for all data. -}
  do (r_sigmaa=distribute(r_sigmaa)) ( &sel )
  do (r_sigmad=distribute(r_sigmad)) ( &sel )
  do (r_dd=distribute(r_dd))         ( &sel )
  do (r_eobs=norm(amplitude(&fobs)))                     ( &sel )
  do (r_ecalc=norm(amplitude(&fcalc+&fpart))  )          ( &sel )

  {- define target and derivative -}
  target = (mlf(amplitude(&fobs),&sigma,(&fcalc+&fpart),r_dd,r_sigmad))
  dtarget = (dmlf(amplitude(&fobs),&sigma,(&fcalc+&fpart),r_dd,r_sigmad))

  if ( &statistics = true ) then
    display
    display  sigmaA statistics                   sigmaA   sigmaD     Delta     <Eobs^2>  <Ecalc^2>         
    statistics
      (r_sigmaa) (r_sigmad) (r_dd) (save(r_eobs^2)) (save(r_ecalc^2)) 
      selection=( &sel )
      output=OUTPUT
    end
  end if

  undeclare name=r_eobs   domain=reciprocal end
  undeclare name=r_ecalc  domain=reciprocal end  
  undeclare name=r_sigmaa domain=reciprocal end
  undeclare name=r_test   domain=reciprocal end

elseif ( &target = "mli" ) then

  declare name=r_test   domain=reciprocal type=real end
  declare name=r_eobs   domain=reciprocal type=real end
  declare name=r_ecalc  domain=reciprocal type=real end
  declare name=r_sigmaa domain=reciprocal type=real end
  query name=r_sigmad domain=reciprocal end
  if ( $object_exist = false ) then
    declare name=r_sigmad domain=reciprocal type=real end
  else
    do (r_sigmad=0) ( all )
  end if
  query name=r_dd domain=reciprocal end
  if ( $object_exist = false ) then
    declare name=r_dd     domain=reciprocal type=real end
  else
    do (r_dd=0) ( all ) 
  end if

  do (r_test=0) ( all )
  do (r_test=1) ( &sel_test )

  show sum (1) (&sel_test)
  if ( $result > 0 ) then
    evaluate ($test_exist=true)
    evaluate ($r_ntest=$result)
  else
    evaluate ($test_exist=false)
    display Error: no reflections flagged for cross-validation
    display        this target requires cross-validation
    abort
  end if

  {- Calculate number of bins used for E and r_sigmaa calculation.
     Make sure there are at least 50 reflections in each bin.
     Issue warning if there end up being less than 50 per bin   -}
  evaluate ($r_mbins=&mbins)
  if ( $r_mbins < 0 ) then
    show sum (1) (&sel_test)
    evaluate ($r_mbins=int($result/50))
    if ( $r_mbins <= 0 ) then
      display Error: there are less than 50 test set reflections
      display        please check the test set and the resolution limits
      abort
    end if
    if ( $r_mbins < 8 ) then
      display Warning: there are less than 50 reflections per bin
      display          you may want to increase the size of the test set
    end if
    evaluate ($r_mbins=max(8,$r_mbins))
    evaluate ($r_mbins=min($r_mbins,50))
  end if
  mbins=$r_mbins  

  {- Check to see that there are not too many bins with no
     test set reflections - limit is 1/5th of total bins -}
  declare name=r_temp domain=reciprocal type=real end
  do (r_temp=sum(r_test)) ( &sel )
  show sum(1) ( &sel )
  evaluate ($ref_per_bin=$select/$r_mbins)
  show sum(1) ( &sel and r_temp=0 )
  evaluate ($empty_bins=int($select/$ref_per_bin))
  undeclare name=r_temp domain=reciprocal end
  evaluate ($empty_ratio=$empty_bins/$r_mbins)
  if ( $empty_ratio >= 0.2 ) then
    display  Warning: too many bins with no test set reflections
    display           Is test set distributed randomly?
  end if

  {- compute r_sigmaa and r_sigmad -}
  if ( $r_ntest < 400 ) then
    do (r_eobs=norm(amplitude(&fobs)))                          ( &sel )
    do (r_ecalc=norm(amplitude(&fcalc+&fpart)))                 ( &sel )
  else
    do (r_eobs=norm(amplitude(&fobs)))                          ( &sel_test )
    do (r_ecalc=norm(amplitude(&fcalc+&fpart)))                 ( &sel_test )
  end if
  do (r_sigmaa=sigacv[sigma=&sig_sigacv](r_eobs,r_ecalc))     ( &sel_test )
  do (r_sigmad=sqrt(save(&iobs) (1-r_sigmaa^2))) ( &sel_test ) 
  {- compute d -}
  do (r_dd= r_sigmaa * sqrt(save(&iobs)  / 
                            save(amplitude(&fcalc+&fpart)^2))) ( &sel_test )
  {- Cross-validation active: need to expand to all selected
     reflections and compute the normalized structure factor
    for all data. -}
  do (r_sigmaa=distribute(r_sigmaa)) ( &sel )
  do (r_sigmad=distribute(r_sigmad)) ( &sel )
  do (r_dd=distribute(r_dd))         ( &sel )
  do (r_eobs=norm(amplitude(&fobs)))                     ( &sel )
  do (r_ecalc=norm(amplitude(&fcalc+&fpart)))            ( &sel )

  {- define target and derivative -}
  target = (mli(&iobs,&sigi,(&fcalc+&fpart),r_dd,r_sigmad))
  dtarget = (dmli(&iobs,&sigi,(&fcalc+&fpart),r_dd,r_sigmad))

  if ( &statistics = true ) then
    display
    display  sigmaA statistics                   sigmaA   sigmaD     Delta     <Eobs^2>  <Ecalc^2>         
    statistics
      (r_sigmaa) (r_sigmad) (r_dd) (save(r_eobs^2)) (save(r_ecalc^2)) 
      selection=( &sel )
      output=OUTPUT
    end
  end if

  undeclare name=r_eobs   domain=reciprocal end
  undeclare name=r_ecalc  domain=reciprocal end  
  undeclare name=r_sigmaa domain=reciprocal end
  undeclare name=r_test   domain=reciprocal end

elseif ( &target = "mlhl" ) then

  declare name=r_test   domain=reciprocal type=real end
  declare name=r_eobs   domain=reciprocal type=real end
  declare name=r_ecalc  domain=reciprocal type=real end
  declare name=r_sigmaa domain=reciprocal type=real end
  query name=r_sigmad domain=reciprocal end
  if ( $object_exist = false ) then
    declare name=r_sigmad domain=reciprocal type=real end
  else
     do (r_sigmad=0) ( all )
  end if
  query name=r_dd domain=reciprocal end
  if ( $object_exist = false ) then
    declare name=r_dd     domain=reciprocal type=real end
  else
    do (r_dd=0) ( all ) 
  end if

  do (r_test=0) ( all )
  do (r_test=1) ( &sel_test )

  show sum (1) (&sel_test)
  if ( $result > 0 ) then
    evaluate ($test_exist=true)
    evaluate ($r_ntest=$result)
  else
    evaluate ($test_exist=false)
    display Error: no reflections flagged for cross-validation
    display        this target requires cross-validation
    abort
  end if

  {- Calculate number of bins used for E and r_sigmaa calculation.
     Make sure there are at least 50 reflections in each bin.
     Issue warning if there end up being less than 50 per bin   -}
  evaluate ($r_mbins=&mbins)
  if ( $r_mbins < 0 ) then
    show sum (1) (&sel_test)
    evaluate ($r_mbins=int($result/50))
    if ( $r_mbins <= 0 ) then
      display Error: there are less than 50 test set reflections
      display        please check the test set and the resolution limits
      abort
    end if
    if ( $r_mbins < 8 ) then
      display Warning: there are less than 50 reflections per bin
      display          you may want to increase the size of the test set
    end if
    evaluate ($r_mbins=max(8,$r_mbins))
    evaluate ($r_mbins=min($r_mbins,50))
  end if
  mbins=$r_mbins  

  {- Check to see that there are not too many bins with no
     test set reflections - limit is 1/5th of total bins -}
  declare name=r_temp domain=reciprocal type=real end
  do (r_temp=sum(r_test)) ( &sel )
  show sum(1) ( &sel )
  evaluate ($ref_per_bin=$select/$r_mbins)
  show sum(1) ( &sel and r_temp=0 )
  evaluate ($empty_bins=int($select/$ref_per_bin))
  undeclare name=r_temp domain=reciprocal end
  evaluate ($empty_ratio=$empty_bins/$r_mbins)
  if ( $empty_ratio >= 0.2 ) then
    display  Warning: too many bins with no test set reflections
    display           Is test set distributed randomly?
  end if

  {- compute r_sigmaa and r_sigmad -}
  if ( $r_ntest < 400 ) then
    do (r_eobs=norm(amplitude(&fobs)))                          ( &sel )
    do (r_ecalc=norm(amplitude(&fcalc+&fpart)))                 ( &sel )
  else
    do (r_eobs=norm(amplitude(&fobs)))                          ( &sel_test )
    do (r_ecalc=norm(amplitude(&fcalc+&fpart)))                 ( &sel_test )
  end if
  do (r_sigmaa=sigacv[sigma=&sig_sigacv](r_eobs,r_ecalc))     ( &sel_test )
  do (r_sigmad=sqrt(save(amplitude(&fobs)^2) (1-r_sigmaa^2))) ( &sel_test ) 
  {- compute d -}
  do (r_dd= r_sigmaa * sqrt(save(amplitude(&fobs)^2)  / 
                            save(amplitude(&fcalc+&fpart)^2))) ( &sel_test )
  {- Cross-validation active: need to expand to all selected
     reflections and compute the normalized structure factor
     for all data. -}
  do (r_sigmaa=distribute(r_sigmaa)) ( &sel )
  do (r_sigmad=distribute(r_sigmad)) ( &sel )
  do (r_dd=distribute(r_dd))         ( &sel )
  do (r_eobs=norm(amplitude(&fobs)))                     ( &sel )
  do (r_ecalc=norm(amplitude(&fcalc+&fpart)))            ( &sel )

  {- define target and derivative -}
  target = (mlhl[phistep=5](amplitude(&fobs),&sigma,
                            phase(&fcalc+&fpart),
                            (&fcalc+&fpart),
                            &pa,&pb,&pc,&pd,r_dd,r_sigmad))
  dtarget = (dmlhl[phistep=5](amplitude(&fobs),&sigma,
                              phase(&fcalc+&fpart),
                              (&fcalc+&fpart),
                              &pa,&pb,&pc,&pd,r_dd,r_sigmad))

  if ( &statistics = true ) then
    display
    display  sigmaA statistics                   sigmaA   sigmaD     Delta     <Eobs^2>  <Ecalc^2>         
    statistics
      (r_sigmaa) (r_sigmad) (r_dd) (save(r_eobs^2)) (save(r_ecalc^2)) 
      selection=( &sel )
      output=OUTPUT
    end
  end if

  undeclare name=r_eobs   domain=reciprocal end
  undeclare name=r_ecalc  domain=reciprocal end  
  undeclare name=r_sigmaa domain=reciprocal end
  undeclare name=r_test   domain=reciprocal end

elseif ( &target = "mixed" ) then

  target=(  resi(amplitude(&fobs),(&fcalc+&fpart),(1-&fom)) +
          vector(combine(amplitude(&fobs),phase(&phase)),
                (&fcalc+&fpart),&fom))
  dtarget=(  dresi(amplitude(&fobs),(&fcalc+&fpart),(1-&fom)) +
           dvector(combine(amplitude(&fobs),phase(&phase)),
                  (&fcalc+&fpart),&fom))

elseif ( &target = "residual" ) then

  target=(resi(amplitude(&fobs),(&fcalc+&fpart),&weight))
  dtarget=(dresi(amplitude(&fobs),(&fcalc+&fpart),&weight))

elseif ( &target = "vector" ) then

  target=(vector(combine(amplitude(&fobs),phase(&phase)),
                (&fcalc+&fpart),&weight))
  dtarget=(dvector(combine(amplitude(&fobs),phase(&phase)),
                  (&fcalc+&fpart),&weight))

elseif ( &target = "f1f1" ) then

  target=(f1f1(amplitude(&fobs),(&fcalc+&fpart)))
  dtarget=(df1f1(amplitude(&fobs),(&fcalc+&fpart)))

elseif ( &target = "f2f2" ) then

  target=(f2f2(amplitude(&fobs),(&fcalc+&fpart)))
  dtarget=(df2f2(amplitude(&fobs),(&fcalc+&fpart)))

elseif ( &target = "e1e1" ) then

  target=(e1e1(amplitude(&fobs),(&fcalc+&fpart)))
  dtarget=(de1e1(amplitude(&fobs),(&fcalc+&fpart)))

elseif ( &target = "e2e2" ) then

  target=(e2e2(amplitude(&fobs),(&fcalc+&fpart)))
  dtarget=(de2e2(amplitude(&fobs),(&fcalc+&fpart)))

else
  
  display Target-error: unknown target definition -> abort
  abort
  
end if

!!if ( &statistics = true ) then
!!  statistics
!!    (rvalue(&fobs,&fcalc+&fpart))
!!    selection=(&sel and &sel_test)
!!    output=OUTPUT
!!  end
!!  statistics
!!    (rvalue(&fobs,&fcalc+&fpart))
!!    selection=(&sel and not &sel_test)
!!    output=OUTPUT
!!  end
!!end if

{- define refinement monitor -}
{- the R-value -}
monitor=(rvalue[overall](&fobs,(&fcalc+&fpart)))

set message=$message_old echo=$echo_old end