#!/usr/bin/python
import os,sys,re,shutil,copy
from process import process, crvapi
import common
par=common.parameter

class comb_phdmmb(process):
  name="combined iterative model building with density modification with phased refinement"
  short_name="combined model building"
  supported_procs = ["dmfull", "mb"]
  supported_params={}
  supported_params['target'] = par( desc='Experiment/refinement target', typ=str, cap=True, share=True )
  supported_params['minbigcyc'] = par( desc='Minimal number of model building cycles', typ=int )
  supported_params['maxbigcyc'] = par( desc='Maximal number of model building cycles', typ=int )
  supported_params['low_res_par'] = par( desc='Use low resolution parameters in refinement and model building.', typ=bool )
  supported_params['refcyc_finish'] = par( desc='Number of refinement cycles in "finishing mode"', typ=int )
  supported_params['ncs_det'] = par( desc='Determine NCS from heavy atoms (Parrot only)', typ=bool, share=True )
  supported_params['ncs_det_ha'] = par( desc='Determine NCS from heavy atoms (Parrot only)', typ=bool, share=True )
  supported_params['ncs_det_mr'] = par( desc='Determine NCS from partial model (Parrot only)', typ=bool, share=True )
  supported_params['no_bias_est'] = par( desc='Do not estimate bias in density modification', typ=bool, share=True )
  supported_params['skip_initial_build'] = par( desc='Skip model bulding in the first "big" cycle', typ=bool )
  supported_params['start_dm'] = par( desc='Include initial phase improvement without model building', typ=bool )
  supported_params['rebuild_only'] = par( desc='Do not include cycles with building "from scratch"', typ=bool )
  supported_params['always_exclude_free'] = par( desc='Exclude inputted free refl. during the whole process (default: in last mb cycle only)', typ=bool )
  supported_params['optimize_solvent'] = par( desc='Disable solvent content optimization', typ=bool )
  supported_params['R_threshold'] = par( desc='Rcomb-factor threshold for switching to finishing mode', typ=float )
  supported_params['fom_threshold'] = par( desc='FOM threshold for switching to finishing mode', typ=float )
  supported_params['max_atoms_from_anom_map'] = par( desc='max. number of atoms to be picked from anomalous diff. map', typ=int )


  def TreatInOutPar(self, set_all_par=False):
    # assign the child processes to specific class attributes for convenience
    if self.GetVirtPar('target') != 'SAD':
      common.Error('Sorry, the "combined" algorithm is only available for SAD as of now ({0} not supported).'.format(
        self.GetVirtPar('target')))
    self.dmbr=self.GetOrAddProcess('dmfull')
    self.ph = self.dmbr.GetOrAddProcess('ref')
    self.dm = self.dmbr.GetOrAddProcess('dm')
    self.mb = self.GetOrAddProcess('mb')
    # check and set some defaults
    if not self.GetParam('minbigcyc'):
      self.SetParam('minbigcyc',5)
    self.minbigcyc_init = self.GetParam('minbigcyc')
    if not self.GetParam('maxbigcyc'):
      self.SetParam('maxbigcyc',50)
    if not self.dmbr.GetParam('dmcyc'):
      self.dmbr.SetParam('dmcyc',8)
    if self.GetVirtPar('low_res_par') is not None and self.ph.GetVirtPar('low_res_par') is None:
      self.ph.SetVirtPar('low_res_par', self.GetVirtPar('low_res_par'))
    self.dmbr.SetRunDir()
    self.dmbr.TreatInOutPar()
    if not self.IsInputtedParam('ncs_det_ha'):
      self.SetParam('ncs_det_ha', self.dmbr.GetParam('ncs_det_ha'))
    self.mb.SetRunDir()
    self.mb.TreatInOutPar()
    #if os.environ['CCP4_LIB'] and os.path.isfile(os.path.join(os.environ['CCP4_LIB'],'data','reference_structures','reference-1ajr.pdb')):
    #  self.mb.GetProg().SetKey('pdbin-ref',os.path.join(os.environ['CCP4_LIB'],'data','reference_structures','reference-1ajr.pdb'))
    #  self.mb.GetProg().SetKey('mtzin-ref',os.path.join(os.environ['CCP4_LIB'],'data','reference_structures','reference-1ajr.mtz'))
    #if not self.mb.GetProg().GetKey('resolution'):
    #  self.mb.GetProg().SetKey('resolution',2.0)
    if not self.GetParam('refcyc_finish'):
      if self.ph.GetVirtPar('low_res_par'):
        self.SetParam('refcyc_finish',15)
      else:
        self.SetParam('refcyc_finish',10)
    if not self.mb.GetProg().GetKey('sequence-reliability'):
      self.mb.GetProg().SetKey('sequence-reliability',0.99)
    if not self.ph.inp.Get('model', typ=('substr','partial+substr'), has_atomtypes=True):
      common.Error('No substructure model or atomtype(s) inputted for {0}.'.format(self.name))
    self.dm_cyc_init=self.dmbr.GetParam('dmcyc')
    self.finish_mode = False
    process.TreatInOutPar(self,set_all_par)
    if self.GetParam('R_threshold'):
      self.R_tresh=self.GetParam('R_threshold')
    if self.GetParam('fom_threshold'):
      self.fom_tresh_max=self.GetParam('fom_threshold')
    if not self.GetParam('always_exclude_free'):
      self.SetParam('always_exclude_free', False)
    if self.GetParam('max_atoms_from_anom_map') in (False,0):
      self.find_new = False

  def RunPreprocess(self,*args,**kwargs):
    process.RunPreprocess(self,*args,**kwargs)
    # remove exclude set as it will be only used for the last cycle by default
    if not self.GetVirtPar('always_exclude_free') and self.inp.Get('exclude',typ='freeR'):
      self.dmbr.inp.Delete(self.inp.Get('exclude',typ='freeR'))
      self.mb.inp.Delete(self.inp.Get('exclude',typ='freeR'))
    # prepare mlhl-comb object
    self.ref_mlhl=self.AddProcess('ref')
    self.ref_mlhl.SetParam('target', 'MLHL')
    self.ref_mlhl.SetParam('cycles', 15)
    self.ref_mlhl.TreatInOutPar()
    # various other preparations
    self.SetParam('minbigcyc', self.GetParam('minbigcyc')*2)
    self.dmbr.separate_models=self.dmbr.AddProcess('sepsubstrprot',propagate_inp=False)
    self.ref_mlhl.separate_models=self.ref_mlhl.AddProcess('sepsubstrprot',propagate_inp=False)
    self.dmbr.logfilehandle=self.logfilehandle
    self.best_pdb=os.path.join(self.rundir,"best.pdb")
    self.best_mtz=os.path.join(self.rundir,"best.mtz")

  def UpdateResults(self, update):
    # defining variables just for the convenience
    mb_stat=self.mb.GetProg(supported=True).GetStat
    ph_stat=self.ph.GetProg(supported=True).GetStat
    if update=='mb':
      self.mb_res_all.append( (self.mb.cyc, mb_stat('res_built'), mb_stat('frag_built'), 
                              mb_stat('compl_chain'), mb_stat('compl_res')) )
    elif update=='ph':
      self.ref_res_all.append( (ph_stat('rfact'), ph_stat('fom')) )
    self.PrintActualLogGraph(update)

  def RunBody(self,*args,**kwargs):
    # prepare phases for initial model building
    if self.GetVirtPar('start_dm'):
      phcmb=self.dmbr.AddProcess('phcomb')
      phcmb.SetParam('target',self.ph.GetParam('target'))
      dmbr_backup=self.dmbr.BackupAnyPars()
      self.dmbr.SetParam('dmcyc',25)
      self.dmbr.Run()
      self.dmbr.RestoreAnyPars(dmbr_backup)
      self.ph.out.Set(phcmb.out.mapcoef)
      self.ph.out.Set(phcmb.out.model)
      self.dmbr.processes.remove(phcmb)
    else:
      self.ph.out.Set(self.dmbr.inp.mapcoef)
      self.ph.out.Set(self.dmbr.inp.model)
    # separate substructure/partial if needed
    obj,N=self.ph.GetProg().inp.GetTargetObjects(self.ph.GetParam('target'),modtyp='pdb')
    if obj['mod'].GetType()=='partial+substr':
      self.dmbr.separate_models.inp.Set(obj['mod'])
      self.dmbr.separate_models.Run()
    self.ph.SetParam('cycles', 1)
    # go with model building!
    for self.mb.cyc in range(1, self.GetParam('maxbigcyc')+1):
      self.out.ClearAll()
      # build model
      if not self.GetVirtPar("skip_initial_build") or self.mb.cyc>1:
        self.Info("\nStarting model building cycle {0}".format(self.mb.cyc))
        self.mb.inp.Set(self.dmbr.out.mapcoef)
        self.mb.GetProg().runname=self.mb.GetProg().name+'_cyc'+str(self.mb.cyc)
        self.mb.Run()
        self.UpdateResults('mb')
      else:
        self.mb.out.Set(self.dmbr.out.Get('model',typ='partial'))
      # setup and run density modification recycling (with bias reduction at specified mb cycles)
      self.dmbr.inp.Set([self.dmbr.out.Get('model',typ='substr',filetype='pdb')]+self.mb.out.model)
      self.dmbr.inp.Set(self.inp.fsigf)
      self.dmbr.inp.Set(self.ph.out.Get('mapcoef',typ=('best','combined')))
      # DM specific parameters
      if self.dmbr.nick == 'dmfull':
        if not self.GetParam('no_bias_est'):
          self.dmbr.SetParam( 'no_bias_est', bool((self.mb.cyc-1)%self.bias_reest_cyc) )
        # ncs parameters will be ignored for solomon
        if self.GetParam('ncs_det') is not False:
          self.dmbr.SetParam( 'ncs_det', self.GetParam('ncs_det') if not (self.mb.cyc-2)%self.bias_reest_cyc else False)
      self.dmbr.Run()
      # update results in case of mlhl-comb ref (dmfull updates results from itself on the fly)
      if self.dmbr.nick == 'ref':
        self.UpdateResults('ph')
      # simple (mainly FOM based) decision making
      self.decide()
      if self.GetParam('minbigcyc')<=0:
        break
      # try some simple extra actions before new bias estimation if we are not successful yet
      if self.mb.cyc%self.bias_reest_cyc==0:
        self.try_something()
    # this file cleanup stuff should be properly organized. at the moment just ad hoc like this.
    os.remove(self.best_mtz)
    os.remove(self.best_pdb)
    #self.Clean(self.dm_be.mtzinitwork)
    self.WriteXMLToContinue()
    # rough FOM based result estimation and printing
    self.FOMresult()
    self.UpdateResults('end')

  def WriteXMLToContinue(self):
    continue_copy = self.from_xml(self.init_xmlelem, dummy=True)
    if continue_copy.GetVirtPar('minbigcyc')<10:
      continue_copy.SetVirtPar('minbigcyc',10)
    continue_copy.inp.Add( self.ph.out.Get('mapcoef',typ=('best','combined')), propagate=False )
    continue_copy.inp.Add( self.ph.out.Get('model',typ='partial+substr'), propagate=False )
    if self.GetCrankParent():
      crank_p=copy.copy(self.GetCrankParent())
    else:
      crank_p=process.from_name('crank',None)
    crank_p.processes=[continue_copy,]
    common.WriteXML(crank_p, os.path.join(self.rundir, self.nick+'_continue.xml'))


  def __init__(self,*args,**kwargs):
    super(comb_phdmmb,self).__init__(*args,**kwargs)
    self.fom_tresh_min = 0.38
    self.fom_tresh_min_init = self.fom_tresh_min
    self.fom_tresh_max=0.60
    self.R_tresh=0.445
    self.bias_reest_cyc=3
    #self.find_rem_heavy=find_remove_refine_heavy(ph)
    self.mb_correlcyc=2
    self.try_dm=0
    self.R_best=1.0
    self.cyc_best=0
    self.mb_res_all, self.ref_res_all = [], []
    self.find_new=1
    #if debug:
    #  self.debug=debug
    #  self.ph.debug=debug
    #  self.dm.debug=debug
    #  self.mb.debug=debug


  def decide(self):
    """FOM and Rcomb based decision making"""
    # for convenience
    ph=self.ph
    dmbr=self.dmbr
    mb=self.mb
    if self.dmbr.nick == 'dmfull':
      # increase weight of refinement as we are getting close
      if ph.R<self.R_tresh+0.07 and ph.fom>=self.fom_tresh_max-0.17 and not self.finish_mode:
        dmbr.SetParam('dmcyc', self.dm_cyc_init-1  if self.dm_cyc_init>1  else 1)
        ph.SetParam('cycles', 2)
      if ph.fom>=self.fom_tresh_max-0.07 and not self.finish_mode:
        dmbr.SetParam('dmcyc', self.dm_cyc_init-2  if self.dm_cyc_init>2  else 1)
        if ph.R<self.R_tresh+0.05:
          ph.SetParam('cycles', 3)
          self.mb_correlcyc=2
      # switch to finishing stage
      if ph.fom>=self.fom_tresh_max and not self.finish_mode:
        if ph.R<self.R_tresh+0.05:
          dmbr.SetParam('dmcyc', self.dm_cyc_init-3  if self.dm_cyc_init>3  else 1)
        if ph.R<self.R_tresh:
          ph.SetParam('cycles', self.GetParam('refcyc_finish'))
          ph.SetParam('beta',1.0)
          dmbr.SetParam('biascyc', 0)
          self.mb_correlcyc=3
          # the number of finishing cycles depends on the number of cycles it took to get to the finishing stage
          self.SetParam( 'minbigcyc', min(self.minbigcyc_init+mb.cyc-1, 2*self.minbigcyc_init) )
          self.finish_mode = 1
          self.Info("Switching to finishing mode ({0} cycles).".format(self.GetParam('minbigcyc')))
    if self.fom_tresh_min<0.5:
      self.fom_tresh_min+=0.005
    self.Info("FOM after {0}. model building cycle is {1}".format(mb.cyc,ph.fom))
    # decrease number of min cycles if we are doing very well or very badly
    if (ph.fom<self.fom_tresh_min or self.finish_mode) and not (ph.R<self.R_best-0.01 and self.GetParam('minbigcyc')<3):
      self.SetParam('minbigcyc', self.GetParam('minbigcyc')-1)
    cyc2go = min(self.GetParam('minbigcyc'), self.GetParam('maxbigcyc')-mb.cyc)
    # "anneal" from time to time
    if mb.cyc%(self.bias_reest_cyc*4)==0 and \
       cyc2go>self.mb_correlcyc+(cyc2go!=(self.GetParam('maxbigcyc')-mb.cyc)):
      self.try_dm=1
      self.fom_tresh_min=self.fom_tresh_min_init+0.05
      #self.SetParam('minbigcyc', self.GetParam('minbigcyc')+1)
    # take the best model to the last cycle
    if self.finish_mode and cyc2go<=1 and self.R_best<ph.R-0.01:
      self.Info("Taking the results from cycle {0} into the last cycle.".format(self.cyc_best))
      shutil.copy(self.best_mtz, dmbr.out.Get('mapcoef',typ='combined').GetFileName('mtz'))
      shutil.copy(self.best_pdb, dmbr.out.Get('model',typ='partial+substr').GetFileName('pdb'))
    # separate substructure/partial if available
    if dmbr.out.Get('model',typ='partial+substr'):
      dmbr.separate_models.inp.Set(dmbr.out.Get('model',typ='partial+substr'))
      dmbr.separate_models.Run()
    # if only substr. was refined then copy it otherwise it might get rewritten by partial+susbtr in the next cycles
    elif dmbr.out.Get('model',typ='substr'):
      newfile=os.path.join(ph.rundir,'substr.pdb')
      shutil.copy(dmbr.out.Get('model',typ='substr').GetFileName('pdb'),newfile)
      dmbr.out.SetFileToChild(dmbr.out.Get('model',typ='substr'),newfile,'pdb')
    # switch between correlation building and standard building
    if mb.GetProg('buccaneer').inp.Get('model',typ=('partial','partial+substr')):
      mb.GetProg('buccaneer').inp.model.remove(mb.GetProg('buccaneer').inp.Get('model',typ=('partial','partial+substr')))
    mb.SetVirtPar('from_weighted_map',False)
    if self.GetParam('rebuild_only') or   ( mb.cyc>1 and \
       (mb.cyc%self.mb_correlcyc!=0 or (cyc2go<self.mb_correlcyc and self.finish_mode) or ph.R<self.R_best-0.01) ):
      if self.out.Get('model',typ='partial'):
        mb.GetProg('buccaneer').inp.Add( self.out.Get('model',typ='partial') )
      if ph.GetVirtPar('low_res_par') and self.finish_mode and mb.cyc%self.mb_correlcyc>1:
        mb.SetVirtPar('from_weighted_map')
    # if we have the best model then save it and take it as result at the end
    # usually the last model is the best but time to time this helps to keep a bit more complete model
    if ph.R<self.R_best:
      self.R_best=ph.R
      self.cyc_best=mb.cyc
      #if self.finish_mode:
      shutil.copy(ph.out.Get('mapcoef',typ='combined').GetFileName('mtz'), self.best_mtz)
      shutil.copy(ph.out.Get('model',typ='partial+substr').GetFileName('pdb'), self.best_pdb)
    # use MLHL-comb (or return to dmfull from MLHL-comb)
    #or (mb.cyc%self.bias_reest_cyc==2 and (self.finish_mode or self.mb.cyc>self.bias_reest_cyc) and mb.cyc%(self.bias_reest_cyc*4)!=(self.bias_reest_cyc*4-1)):
    if self.finish_mode and cyc2go<=1 and self.GetParam('maxbigcyc')>1:
      if self.inp.Get('exclude',typ='freeR') and not self.GetParam('always_exclude_free'):
        self.ref_mlhl.inp.Set(self.inp.Get('exclude',typ='freeR'))
        #self.mb.inp.Set(self.inp.Get('exclude',typ='freeR'))
      self.ref_mlhl.out.Add(dmbr.out.Get('model',typ='substr'))
      self.ref_mlhl.out.Add(dmbr.out.Get('model',typ='partial'))
      self.ref_mlhl.out.Set(dmbr.out.Get('mapcoef',typ='combined'))
      self.ph = self.ref_mlhl
      self.dmbr = self.ref_mlhl
    elif self.dmbr.nick != 'dmfull':
      self.dmbr = self.GetProcess('dmfull')
      self.ph = self.dmbr.GetProcess('ref')
      self.dmbr.out.Add(self.ref_mlhl.out.Get('model',typ='substr'))
      self.dmbr.out.Add(self.ref_mlhl.out.Get('model',typ='partial'))
      self.ph.out.Set(self.ref_mlhl.out.Get('mapcoef',typ=('combined','best')))


  def FOMresult(self):
    """rough FOM based result estimation and its printing"""
    if self.ph.fom<self.fom_tresh_min:
      self.result_str="Wrong substructure or very weak phases suspected: the model building seems unsuccessful."
    elif self.finish_mode:
      self.result_str="Majority of model should be successfully built."
    else:
      self.result_str="A partial (or no) model is probably correctly built."
    self.Info('\n'+self.result_str)

  def PrintActualLogGraph(self, update):
    if self.opened_loggraph:
      self.GetLogGraphHandle().seek(0,0)
      self.LGInfo(self.GetCCP4Header())
      self.LGInfo(self.GetLogGraphPrograms())
      if self.ref_res_all:
        self.LGInfo('\n $TABLE : Rcomb-factor and FOM per cycle:')
        self.LGInfo('$GRAPHS    : Rcomb-factor and FOM per cycle :A:1,2,3:\n$$')
        self.LGInfo('Cycle  Rcomb-factor   FOM  $$ $$')
        for cyc,(R,fom) in enumerate(self.ref_res_all):
          self.LGInfo('{0:3}  {1:8.3f}  {2:8.3f}'.format(cyc+1,R,fom))
        self.LGInfo('$$\n')
      if self.mb_res_all:
        self.LGInfo('\n $TABLE : Built per building cycle:')
        self.LGInfo('$GRAPHS    : Residues and fragments per cycle :A:1,2,3:\n$$')
        self.LGInfo('Cycle Residues Fragments  $$ $$')
        for cyc,res,fr,c1,c2 in self.mb_res_all:
          self.LGInfo('{0:3}  {1:8}  {2:8}'.format(cyc,res,fr))
        self.LGInfo('$$\n')
      if update=='end':
        self.LGInfo('\n$TEXT:Result: $$ Final result $$')
        self.LGInfo("{0}\n$$\n".format(self.result_str))
    if self.rv_report is not None:
      if self.mb_res_all and update=='mb':
        if not hasattr(self,'rv_plotmb'):
          self.rv_plotmb = self.rv_report.Plot( 'Residues and fragments vs build cycle', "Cycle", "Number of built residues and fragments", block="Model building statistics" )
          self.rv_res = self.rv_plotmb.PlotLine(["x"],["residues"])
          self.rv_fr = self.rv_plotmb.PlotLine(["x"],["fragments"])
        cyc,res,fr,c1,c2 = self.mb_res_all[-1]
        self.rv_res.Data( cyc, res, int ),  self.rv_fr.Data( cyc, fr, int, flush=True )
      if self.ref_res_all and update=='ph':
        if not hasattr(self,'rv_plotref'):
          self.rv_plotref = self.rv_plotmb.parent.parent.Plot( 'Rcomb and FOM vs refinement cycle', "Cycle", "Rcomb-factor and FOM", block="Refinement statistics" )
          self.rv_fom = self.rv_plotref.PlotLine(["x"],["FOM"])
          self.rv_R = self.rv_plotref.PlotLine(["x"],["Rcomb"])
          # this displays the actual pdb+mtz updated cycle by cycle
          # the files are continuously rewritten which might cause issues - we'll need to test!
          out_pdb = self.out.Get('model',typ=('partial','partial+substr')).GetFileName()
          out_mtz = self.out.Get('mapcoef',typ='combined').GetFileName()
          self.rv_report.Text("<BR><i><small>Output (updating after each refinement cycle):</small></i>")
          self.rv_files = self.rv_report.DataFile(out_pdb, "xyz", "<small>Built model and map</small>")
          self.rv_files.DataFile(out_mtz, "hkl:map", flush=True)
        R,fom = self.ref_res_all[-1]
        self.rv_fom.Data( len(self.ref_res_all), fom, int ),  self.rv_R.Data( len(self.ref_res_all), R, int, flush=True )
      if update=='end':
        self.rv_report.Text("<i><small>Result:</small></i>")
        self.rv_report.Text('&emsp;'+self.result_str, flush=True)
        #self.rv_report.Text("<BR><i><small>Output:</small></i>")
        #self.rv_files = self.rv_report.DataFile(out_pdb, "xyz", "<small>Built model and map</small>")
        #self.rv_files.DataFile(out_mtz, "hkl:map", flush=True)

  def try_something(self):
    # try searching for new atoms and removing low occupancy ones 
    # restricted to only once in the building due to too many related fake atoms that freeze NCS detection by Parrot
    if self.find_new and not self.finish_mode and self.mb.cyc<=self.bias_reest_cyc:
      # find new atoms
      atomsfrommap=self.AddProcess('atomsfrommap')
      atomsfrommap.SetParam('rms_threshold', 6.0)
      if self.GetParam('max_atoms_from_anom_map'):
        atomsfrommap.SetParam('max_new_atoms', self.GetParam('max_atoms_from_anom_map'))
      atomsfrommap.inp.Add(self.ph.GetProg().out.Get('mapcoef',typ='anom-diff'))
      atomsfrommap.inp.Add(self.dmbr.out.Get('model',typ='substr',filetype='pdb'))
      atomsfrommap.Run()
    # remove low occ atoms
    pdbcur=self.AddProg('pdbcur')
    pdbcur.runname=pdbcur.name+'_removelowocc'
    pdbcur.inp.Add(self.out.Get('model',typ='substr',filetype='pdb'))
    pdbcur.SetKey('cutocc',0.02)
    pdbcur.Run()
    self.ph.GetProg().out.Add(pdbcur.out.model[-1])
    self.programs.remove(pdbcur)
    #if not self.finish_mode and self.mb.cyc==2*self.bias_reest_cyc:
    #  # try increasing the DM filtering radius
    #  if hasattr(self.dmbr,'dm') and self.dmbr.dm.GetProg().nick=='parrot':
    #    self.dmbr.dm.GetProg().SetKey( 'solvent-mask-filter-radius', self.dmbr.dm.GetProg().GetStat('radius_auto')*2 )
    # try including standard DM if DM2 might got stuck
    if self.try_dm:
      self.Info("Phasing-DM annealing will follow.")
      # first run 10 cycles of current model refinement without phase combination
      ph_backup=self.ph.BackupAnyPars()
      self.ph.SetParam('cycles',10)
      self.ph.SetParam('phcomb',False)
      self.ph.SetParam('beta',1.0)
      self.ph.GetProg().runname=self.ph.GetProg().nick+'_ref'
      self.ph.Run()
      self.ph.RestoreAnyPars(ph_backup)
      # now run standard DM
      phcomb=self.dmbr.AddProcess('phcomb')
      phcomb.SetParam('target',self.ph.GetParam('target'))
      phcomb.SetParam('cycles',0)
      phcomb.AddProg(self.ph.GetProg().nick)
      dmbr_backup=self.dmbr.BackupAnyPars()
      self.dmbr.SetParam('dmcyc',15)
      self.dmbr.SetParam('biascyc',4)
      self.dmbr.SetParam('no_bias_est',False)
      # try solvent optimization - unless the user inputted solvent content or number of monomers
      mon_obj=self.dmbr.inp.Get(has_monomers_asym=True)
      solv_obj=self.dmbr.inp.Get(has_solvent_content=True)
      if self.mb.cyc<=5*self.bias_reest_cyc and self.GetParam('optimize_solvent') is not False and \
         not self.dmbr.IsInputtedParam('solvent_content') and \
         (not mon_obj or not mon_obj.IsInputtedAttr('monomers_asym')) and \
         (not solv_obj or not solv_obj.IsInputtedAttr('solvent_contest')):
        self.dmbr.SetParam('optimize_solvent')
      #self.dmbr.SetParam('ncs_det', None)
      self.dmbr.Run()
      self.ph.out.Set(phcomb.out.mapcoef,propagate=False)
      self.dmbr.processes.remove(phcomb)
      self.dmbr.RestoreAnyPars(dmbr_backup)
      self.try_dm=0