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

class dmfull(process):
  name="density modification with Fourier recycling"
  short_name="density modification"
  supported_progs=["parrot", "shelxe"]
  supported_procs=["dm", "phcomb", "ref"]
  supported_params = {}
  supported_params['target'] = par( desc='Experiment/refinement target', typ=str, cap=True, share=True )
  supported_params['dmcyc'] = par( desc='Number of density modification iterations', typ=int )
  supported_params['biascyc'] = par( desc='Number of iterations used for bias estimation', typ=int )
  supported_params['beta'] = par( desc='Bias parameter beta', typ=float, share=True )
  supported_params['solvent_content'] = par( desc='(Expected) solvent fraction of crystal', typ=float, share=True )
  supported_params['no_bias_est'] = par( desc='Skip bias estimation', typ=bool )
  supported_params['ncs_det'] = par( desc='Determine NCS from heavy atoms or partial model (Parrot only)', typ=bool )
  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 )
  supported_params['optimize_solvent'] = par( desc='Try to optimize solvent content fraction', typ=bool )
  supported_params['incr_radius_threshold_fom'] = par( desc='Filtering radius is increased if FOM after cycle 5 is below this value (Parrot only)', typ=float )


  def TreatInOutPar(self, set_all_par=False):
    # defaults if no programs and no subprocesses are specified
    if not self.GetProg(supported=True) and not self.GetProcess(supported=True):
      if self.GetVirtPar('target') in ('SAD','SIRAS'):
        self.AddProcess('dm')
        self.AddProcess('phcomb')
      else:
        self.AddProg('parrot')
    if self.GetProg(supported=True):
      # any complete DM with recycling within a single program stuff shall go here
      if not self.IsInputtedParam('no_bias_est'):
        self.SetParam('no_bias_est', True)
    else:
      # parameters specific for our multivariate and/or bias reduced density modification recycling
      # many of these parameters are our DM specific, so here we also set their defaults
      # now first assign the child processes to specific class attributes for convenience
      # and set defaults if an incomplete subprocesses list is specified:
      self.dm = self.GetOrAddProcess('dm')
      self.ph = self.GetProcess('phcomb')
      if not self.ph:
        self.ph = self.GetProcess('ref')
      if not self.ph:
        self.ph = self.AddProcess('phcomb')
      if self.ph.nick=='ref':
        self.ph.SetVirtPar('phcomb')
      # check and set some defaults
      if self.GetParam('dmcyc') is None:
        self.SetParam('dmcyc',25)
        if self.parent_process and self.parent_process.nick=='handdet':
          self.SetParam('dmcyc',10)
      # bias estimation parameters
      if self.GetParam('biascyc') is None:
        self.SetParam('biascyc',4)
      if self.GetParam('no_bias_est') is None and (self.GetParam('beta') \
         or (self.parent_process and self.parent_process.nick=='handdet')):
        self.SetParam('no_bias_est',True)
      # disable ncs detection from heavy atoms by default if it would likely take too long
      if not self.IsInputtedParam('ncs_det_ha') and self.GetParam('ncs_det') is not False:
        self.SetParam('ncs_det_ha', True)
        # real time decision making - disabled due to problems with how to present in the GUI...
        #if self.inp.Get('model',typ=('substr','partial+substr'),has_atomtypes=True,filetype='pdb'):
        #  # determine the actual number of atoms in the pdb
        #  pdbcur=self.AddProg('pdbcur', propagate_out=False)
        #  for at in self.inp.Get('model',typ=('substr','partial+substr'),has_atomtypes=True).GetAtomTypes():
        #    pdbcur.SetKey('delatom', "{0}[{0}]:*".format(at))
        #  pdbcur.Run()
        #  if pdbcur.GetStat('atoms_deleted')>35:
        #    self.SetParam('ncs_det_ha', False)
        #  self.programs.remove(pdbcur)
        # this should be a reasonable solution in the case of substructure detection...
        if self.inp.Get(has_num_atoms=True) and self.inp.Get(has_num_atoms=True).exp_num_atoms>30:
          self.SetParam('ncs_det_ha', False)
        if self.GetParam('incr_radius_threshold_fom') is None:
          self.SetParam('incr_radius_threshold_fom', 0.4)
      self.ph.TreatInOutPar(set_all_par)
      self.dm.TreatInOutPar(set_all_par)
      # this will have to be enabled for free lunch or other extension schemes
      #self.dm.GetProg().keep_nodata=True
    process.TreatInOutPar(self,set_all_par)


  def RunBody(self,*args,**kwargs):
    # if program is to be run directly (eg "original" parrot, shelxe)
    if any(p for p in self.programs if p.nick in self.supported_progs):
      process.RunBody(self)
    # dm program and phase comb program combination
    else:
      # make sure there is container with input data mtz file in self.inp 
      # (copy it from program input if needed) - ie the dmfull input data mtz is saved here and used later
      if self.dm.GetProg().inp.Get('fsigf',typ='average',col='f'):
        self.inp.AddCopy( self.dm.GetProg().inp.Get('fsigf',typ='average',col='f'), propagate=False )
      else:
        common.Error('(Average) F/SIGF is not available for input to {0}'.format(self.name))
      self.inp.SetCopy( self.dm.GetProg().inp.GetAll('mapcoef',stored_order=True), propagate=False )
      self.CreateWorkingSetMTZ(new_seed=True)
      if self.GetVirtPar('optimize_solvent'):
        dmf_copy=self.GetProcessCopy()
      self.RunComb()
      if self.GetVirtPar('optimize_solvent'):
        self.OptimizeSolvent(dmf_copy)
      #self.Clean(mtz_in_saved)
      #if dm_be is not None:
      #  self.Clean(dm_be.mtzinitwork)

  def OptimizeSolvent(self, dmf_copy):
    #tries to optimize solvent by FOM based multiple jobs with different solvent content
    #requires a dmfull object on input which is assumed to be run already (thus having self.fom*)
    # first make sure all is available and consistent
    if not self.dm.inp.Get('sequence'):
      common.Warning('Solvent optimization will be skipped as sequence not known.')
      return
    if self.dm.inp.Get(has_monomers_asym=True):
      num_mon = self.dm.inp.Get(has_monomers_asym=True).monomers_asym
    else:
      matthews=dmf_copy.GetProcess(self.dm.nick).GetOrAddProcess('matthews')
      matthews.out.never_propagate=True
      matthews.Run()
      num_mon = matthews.GetProg().GetStat('monomers_asym')
      if self.dm.GetParam('solvent_content') and \
        abs(matthews.GetProg().GetStat('solvent_content'),self.dm.GetParam('solvent_content'))>0.1:
          common.Warning('Solvent optimization will be skipped as inputted solvent content differs too much from the estimation.')
          return
    # now try the optimization
    num_mon_start=num_mon
    fom_best=sum(self.ph.fom_all[5:10])
    dmf_best=self
    num_mon_change=-1
    if num_mon<=1:
      num_mon_change=1
    while num_mon_change!=0:
      num_mon = num_mon+num_mon_change
      dmf_new = dmf_copy.GetProcessCopy()
      dmf_new.Info("\nTesting solvent content corresponding to {0} monomers.".format(num_mon))
      dmf_new.out.never_propagate=True
      dmf_new.SetRunDir(dmf_new.rundir+'_'+str(num_mon))
      dmf_new.ph=dmf_new.GetProcess(self.ph.nick)
      dmf_new.dm=dmf_new.GetProcess(self.dm.nick)
      dmf_new.sft=dmf_new.GetProcess('createfree').GetProg('sftools')
      dmf_new.sft.seed=self.sft.seed
      dmf_new.sft.SetRunDir(os.path.join(dmf_new.rundir,'betafromfree'))
      dmf_new.dm.inp.Get('sequence').monomers_asym=num_mon
      dmf_new.inp.Set(dmf_new.dm.inp.Get('sequence'))
      matthews=dmf_new.dm.GetOrAddProcess('matthews')
      try:
        matthews.Run()
      except (matthews.GetProg().ProgramRunError,):
        self.Info("Matthews coeficients for {0} monomers too unlikely.".format(num_mon))
        fom=0.0
      else:
        dmf_new.RunComb()
        fom=sum(dmf_new.ph.fom_all[5:10])
      if abs(fom-fom_best)<0.02:
        self.Info("Very weak discrimination - the solvent content estimation may be unreliable.")
      if fom>fom_best:
        fom_best=fom
        dmf_best=dmf_new
        if num_mon<=1:
          num_mon_change=0
      elif num_mon_change<0 and dmf_best is self:
        num_mon = num_mon_start
        num_mon_change=1
      else:
        num_mon_change=0
    # evaluate and copy the best results
    if dmf_best is self:
      self.Info("More likely solvent content not found.")
    else:
      self.dm.inp.Set( dmf_best.dm.inp.Get('sequence') )
      self.dm.out.Add( dmf_best.dm.inp.Get('sequence') )
      self.dm.out.Add( dmf_best.out.Get('mapcoef',typ='densmod') )
      self.ph.out.Add( dmf_best.out.Get('mapcoef',typ='combined') )
      self.Info("Solvent content of {0} corresponding to {1} monomers will be used.".format(
        self.dm.inp.Get('sequence').solvent_content, self.dm.inp.Get('sequence').monomers_asym))

  def RunPostprocess(self,restore=True,*args,**kwargs):
    if self.GetProcess('createfree'):
      self.processes.remove(self.createfree)
    self.PrintActualLogGraph(finished=True)
    # a workaround for multicomb replacing the residue type (for example IOD always becomes I)
    # this can (and is better to) be removed after this is fixed in multicomb
    if self.GetProcesses(supported=True) and self.ph.GetProg().nick=='multicomb' and \
       self.ph.out.Get('model',typ='substr',filetype='pdb'):
      fixsub=self.ph.AddProcess('fixsubstrpdb',propagate_inp=False)
      fixsub.inp.Add(self.ph.out.Get('model',typ='substr',filetype='pdb'))
      fixsub.AddProg('pdbcur')
      fixsub.GetProg().outfilename['pdb'] = "fixed_"+self.ph.nick+".pdb"
      fixsub.Run()
      self.ph.processes.remove(fixsub)
    process.RunPostprocess(self,restore,*args,**kwargs)

  def PrintActualLogGraph(self, finished=False, beta=False):
    if self.opened_loggraph:
      self.GetLogGraphHandle().seek(0,0)
      self.LGInfo(self.GetCCP4Header())
      self.LGInfo(self.GetLogGraphPrograms())
      if not self.GetVirtPar('no_bias_est'):
        self.LGInfo('\n$TEXT: Beta estimation result: $$ Beta estimation result $$')
        self.LGInfo('Beta bias parameter estimated as {0}\n$$'.format(self.ph.GetParam('beta')))
      if self.GetProcesses(supported=True) and self.ph.fom_all:
        self.LGInfo('\n $TABLE : Fom per cycle:')
        self.LGInfo('$GRAPHS    : FOM per cycle :A:1,2,3:\n$$')
        self.LGInfo('Cycle FOM $$ $$')
        for i,n in enumerate(zip(self.ph.fom_all)):
          self.LGInfo('{0:3}  {1:8.3f}'.format(i+1,n[0]))
        self.LGInfo('$$\n')
        if finished:
          self.LGInfo('\n$TEXT:Result: $$ Final result $$')
          self.LGInfo('The final FOM is {0}\n$$'.format(self.ph.fom_all[-1]))
    if self.rv_report is not None:
      if beta:
        self.rv_report.Text("<BR><i><small>Beta (bias) parameter estimation result:</small></i>")
        self.rv_report.Text("&emsp;Beta parameter estimated as {0}".format(self.ph.GetParam('beta')), flush=True)
      elif self.GetProcesses(supported=True) and self.ph.fom_all:
        if not hasattr(self,'rv_plot'):
          self.rv_plot = self.rv_report.Plot( 'FOM vs cycle', "Cycle", "Figure of merit")
          self.rv_fom = self.rv_plot.PlotLine( ["x"], ["FOM"] )
        #else:
        self.rv_fom.Data( len(self.ph.fom_all), self.ph.fom_all[-1], flush=True )
        if finished:
          self.rv_report.Text("<i><small>Result:</small></i>") # title
          self.rv_report.Text('&emsp;'+'The final FOM is {0}'.format(self.ph.fom_all[-1]))
          self.rv_report.Text("<BR><i><small>Output:</small></i>")
          self.rv_files = self.rv_report.DataFile(self.out.Get('model').GetFileName(), "xyz", "<small>DM improved map</small>")
          self.rv_files.DataFile(self.out.Get('mapcoef',typ='combined').GetFileName(), "hkl:map", flush=True)

  def RunComb(self):
    # contin is used so that the process can be repeated if needed
    contin=1
    num_failed_bias=0
    while contin:
      if hasattr(self,'rv_plot'):
        self.rv_plot.parent.parent.Remove()
        del self.rv_plot
      # estimate density modification bias
      extra_corr=False
      if not self.GetParam('no_bias_est') and self.GetParam('biascyc')>0:
        if num_failed_bias>3:
          self.Info('Beta estimation seems unstable and will be skipped, with beta set to 1.\n')
          self.ph.SetParam('beta',1.0)
          self.SetParam('no_bias_est', True)
        else:
          # possibly restarting due to detected NCS estimation issues
          if self.RunWithGivenBeta(bias_est=True):
            continue
          self.EstBias()
          # if beta is low then we should estimate the correction separately - ie more cycles needed...
          if self.ph.GetParam('beta')<0.7 or self.GetParam('optimize_solvent'):
            extra_corr=True
      # run density modification + check bias estimation, apply correction (if bias is estimated)
      if self.RunWithGivenBeta(extra_corr=extra_corr):
        # if bias estimation problem then let's reseed and repeat bias estimation
        self.CreateWorkingSetMTZ()
        num_failed_bias+=1
      else:
        contin=0

  def ReInitPhDM(self):
    self.dm.inp.SetCopy(self.inp.fsigf)
    self.dm.inp.SetCopy(self.inp.mapcoef)
    self.dm.inp.SetCopy(self.inp.model)
    # this will have to be disabled for free lunch or other extension schemes!
    # done here as a lot of extra nodata can hinder the dm performance.
    self.dm.GetProg().keep_nodata=False
    self.ph.inp.SetCopy(self.inp.fsigf)
    self.ph.inp.SetCopy(self.inp.model)
    self.ph.inp.SetCopy(self.inp.mapcoef)

  def RunWithGivenBeta(self, bias_est=False, extra_corr=False, no_corr=False, verbose=1):
    """Simple DM+phase combination recycling"""
    self.ph.R_all, self.ph.corr_all, self.ph.fom_all, self.dm.contr_inv_all = [], [], [], []
    beta_saved = 0
    self.ReInitPhDM()
    if bias_est or extra_corr:
      ncyc = self.GetParam('biascyc')
      if extra_corr:
        beta_saved=self.ph.GetParam('beta')
      self.ph.SetParam('beta',1.0)
    else:
      ncyc = self.GetParam('dmcyc')
    if verbose and ncyc:
      if bias_est:
        self.Info("Starting {0} bias estimation cycles\n".format(ncyc))
      elif extra_corr:
        self.Info("Starting {0} bias correction cycles\n".format(ncyc))
      else:
        #refi="and refinement" if mb_mod else ""
        self.Info("Starting {0} density modification cycles\n".format(ncyc))
    self.ph.out.Set(self.ph.inp.model, propagate=False)
    # the recycling starts here
    for self.cyc in range(1, ncyc+1):
      if verbose:
        self.Info(" Cycle {0}".format(self.cyc))
      self.out.ClearAll()
      self.SetNCSParams()  # parrot only
      self.dm.GetProg().runname=self.dm.GetProg().name+'_cyc'+str(self.cyc)
      self.dm.Run()
      self.ph.inp.AddCopy(self.dm.out.Get('mapcoef',typ='densmod'))
      if self.ph.out.model: # take the latest model if available (not available for MLHL)
        self.ph.inp.Set(self.ph.out.model)
      self.ph.GetProg().runname=self.ph.GetProg().name+'_cyc'+str(self.cyc)
      self.ph.Run()
      if self.ph.GetProg().nick=='refmac': 
        self.ph.inp.Clear('datafile')
        if self.cyc>2 and self.cyc<ncyc:
          self.ph.inp.Set(self.ph.out.Get('datafile',typ='dluz'))
      if not self.dm.GetProg('parrot'):
        self.dm.inp.Clear('mapcoef')
      self.dm.inp.AddCopy(self.ph.out.Get('mapcoef',typ='weighted'))
      # substr. model can be used for ncs detection parrot (otherwise it is ignored)
      if self.out.Get('model',typ='substr'):
        self.dm.inp.Set(self.out.Get('model',typ='substr'))
      if bias_est and self.cyc>0:
        self.dm.GetProg().keep_nodata=True
        self.MakeFreeDMInpAbsent()
      self.GetActualStats(verbose, update_loggraph=not bias_est and not extra_corr)
      if self.GetNCSParams():  # parrot only
        self.ReInitPhDM()
        return 1
      # automatically increase the filtering radius for weak phases (only for parrot as of now)
      if self.dm.GetProg('parrot') and not self.dm.GetProg().GetKey('solvent-mask-filter-radius') \
         and self.cyc==5 and self.ph.fom<self.GetParam('incr_radius_threshold_fom'):
        self.dm.GetProg().SetKey( 'solvent-mask-filter-radius', self.dm.GetProg().GetStat('radius_auto')*1.5 )
      if not bias_est and not self.GetParam('no_bias_est') and not no_corr and \
         self.cyc==self.GetParam('biascyc') and self.EstBias(correction=True,extra_correction=beta_saved):
        self.ReInitPhDM()
        return 1
    if extra_corr:
      if self.RunWithGivenBeta(no_corr=True):
        self.ReInitPhDM()
        return 1
    self.ReInitPhDM()
    return 0

  def GetActualStats(self,verbose,update_loggraph):
    self.ph.fom = self.ph.GetProg().GetStat('fom')
    self.ph.fom_all.append(self.ph.fom)
    if verbose:
      self.Info("  Overall MEAN FOM is {0}".format(self.ph.fom))
    if self.ph.nick=='ref':
      self.ph.R = self.ph.GetProg().GetStat('rfact')
      self.ph.R_all.append(self.ph.R)
      if verbose:
        self.Info("  Overall Rcomb-factor is {0}".format(self.ph.R))
    else:
      self.ph.corr = self.ph.GetProg().GetStat('correl')
      self.ph.corr_all.append(self.ph.corr)
      if verbose and (self.ph.nick=='ref' or self.ph.GetProg().nick!='refmac'):
        self.Info("  Overall Correlation is {0}".format(self.ph.corr))
      if self.dm.GetProg('solomon'):
        self.dm.contr_inv=self.dm.GetProg('solomon').GetStat('contrast_inv')
        self.dm.contr_inv_all.append(self.dm.contr_inv)
        if verbose:
          self.Info("  Inverse contrast is {0}".format(self.dm.contr_inv))
    if verbose and update_loggraph:
      self.PrintActualLogGraph()
      if self.parent_process and self.parent_process.nick=='comb_phdmmb' and self.ph.nick=='ref':
        self.parent_process.UpdateResults('ph')

  def SetNCSParams(self):
    if not self.dm.GetProg('parrot') or self.GetParam('ncs_det') is False:
      self.dm.SetParam('ncs_det', False)
      return
    # we do not want to run ncs determination every cycle
    if (self.cyc-2)%10==0 or self.cyc==4:
      self.dm.SetParam('ncs_det', self.GetParam('ncs_det'))
      self.dm.GetProg('parrot').SetKey('ncs-operator',False,keep_previous=False)
    else:
      self.dm.SetParam('ncs_det', False)
    if self.cyc==1:
      self.dm.GetProg('parrot').SetKey('ncs-operator',False,keep_previous=False)

  def GetNCSParams(self):
    #if not self.dm.GetProg('parrot') or self.GetParam('ncs_det') is False:
    if not self.dm.GetProg('parrot'):
      return 0
    #if self.dm.GetParam('ncs_det') is not False:
    ncs_opers=self.dm.GetProg('parrot').GetStat('ncs_operator')
    self.dm.GetProg('parrot').SetKey('ncs-operator',False,keep_previous=False)
    for ncs in ncs_opers:
      self.dm.GetProg('parrot').SetKey('ncs-operator',ncs)
    if ncs_opers and self.dm.GetParam('ncs_det') is not False:
      self.Info('  {0} NCS operators detected by {1}.'.format(len(ncs_opers), self.dm.GetProg('parrot').name))
    # if NCS detection seems to go wrong disable it; this assumes NCS detection is run every (self.cyc-2)%10==0 cycle
    # also this assumes the phases should improve (ie starting phases are from phasing)!
    if self.GetParam('ncs_det') is None and \
       self.ph.nick=='phcomb' and not self.GetParam('optimize_solvent') and \
       (self.cyc-3)%10==0 and self.ph.fom_all[self.cyc-1]+0.003<self.ph.fom_all[self.cyc-3]:
      self.Info("Disabling automatic NCS detection and restarting DM.\n")
      self.SetParam('ncs_det', False)
      self.dm.GetProg('parrot').SetKey('ncs-operator',False,keep_previous=False)
      return 1
    return 0

  def EstBias(self, correction=False, extra_correction=False): 
    """Estimate DM bias from performed bias estimation jobs"""
    # run sftools to calculate the correlations
    self.sft.ClearAnyParams()
    self.sft.runname=self.sft.name+'_correl'
    if not correction and not extra_correction:
      self.sft.runname+='_betaest'
    if extra_correction:
      self.ph.SetParam('beta',extra_correction)
    #self.sft.MergeMTZ(*self.sft.inp.GetAll())
    exclude_sft=self.sft.out.Get('exclude',col='free')
    mapc_dm_out=self.dm.out.Get('mapcoef',typ='densmod',filetype='mtz',col='f')
    mapc_obs_in=self.inp.Get('fsigf',typ='average',filetype='mtz',col='f')
    self.sft.SetKey( 'read', (exclude_sft.GetFile('mtz').name, 'col', "\"{0}\"".format(exclude_sft.GetLabel('free'))) )
    self.sft.SetKey( 'read', (mapc_dm_out.GetFile('mtz').name, 'col', "\"{0}\"".format(mapc_dm_out.GetLabel('f'))) )
    self.sft.SetKey( 'read', (mapc_obs_in.GetFile('mtz').name, 'col', "\"{0}\"".format(mapc_obs_in.GetLabel('f'))) )
    self.sft.SetKey( 'select', ('col', exclude_sft.GetLabel('free'), '=1') )
    self.sft.SetKey( 'checkhkl' )
    self.sft.SetKey( 'correl', ('col', mapc_dm_out.GetLabel('f')+' '+mapc_obs_in.GetLabel('f')) )
    self.sft.SetKey( 'select', 'invert' )
    self.sft.SetKey( 'checkhkl' )
    self.sft.SetKey( 'correl', ('col', mapc_dm_out.GetLabel('f')+' '+mapc_obs_in.GetLabel('f')) )
    self.sft.SetKey( 'quit\nY' )
    self.sft.Run(restore=False, clear_out=False)
    # get the correlations from the log
    with open(self.sft.GetLogFileName(),"r") as f:
      corr=re.findall("OVERALL STATISTICS\s+\S+\s+(\S+)", f.read())
      if len(corr)!=2:
        common.Error( "Unexpected number of matches in {0}{1}: {2}".format(
          self.sft.runname,self.sft.log_suffix,len(corr)) )
      corr=(float(corr[0]), float(corr[1]))
    # now calculate the estimation or correction to the estimation
    if not correction:
      self.ph.SetParam('beta',round(corr[1]/corr[0],3))
    else:
      perc=corr[1]/corr[0]-1.0
      if (perc>0.15 or perc<-0.15) and (corr[0]>0.3 or perc>0.3 or perc<-0.3):
        self.Info('Too unrepresentative free set. Will need to generate new free set and start from beginning again.\n')
        return 1
      else:
        self.ph.SetParam('beta', round(self.ph.GetParam('beta')*corr[0]/corr[1],3))
    if self.ph.GetParam('beta')<0.2:
      self.ph.SetParam('beta',0.2)
    elif self.ph.GetParam('beta')>1.0:
      self.ph.SetParam('beta',1.0)
    if correction:
      self.Info("Beta corrected to {0}\n".format(self.ph.GetParam('beta')))
    else:
      self.Info("Beta estimated as {0}\n".format(self.ph.GetParam('beta')))
      self.PrintActualLogGraph(beta=True)
    return 0


  def CreateWorkingSetMTZ(self,new_seed=False):
    """Create MTZ containing only working set (for bias estimation) reflections"""
    self.createfree=self.GetOrAddProcess('createfree')
    self.createfree.SetVirtPar('freetype','freebias')
    self.createfree.inp.Set(self.inp.Get('fsigf',typ='average',filetype='mtz'))
    self.createfree.SetRunDir('betafromfree')
    self.sft=self.createfree.GetOrAddProg('sftools')
    self.sft.SetRunDir()
    # this makes sure that nodata are purged
    self.sft.MergeMTZ( self.createfree.inp.Get(filetype='mtz'), force_merge=True )
    self.sft.ClearAnyParams()
    self.sft.runname=self.sft.name+'_createfree'
    if new_seed:
      self.sft.seed=0
    else:
      self.sft.seed+=1
    self.sft.SetKey( 'calc', ('seed', self.sft.seed), keep_previous=False )
    self.createfree.Run()

  def MakeFreeDMInpAbsent(self):
    exclude_sft=self.sft.out.Get('exclude',col='free')
    self.sft.ClearAnyParams()
    self.sft.inp.ClearAll()
    self.sft.runname=self.sft.name+'_workonly'
    # catch all objects needed for input of the DM program
    self.dm.GetProg().BackupAnyPars()
    for o in self.dm.GetProg().CatchMTZObjects():
      self.sft.inp.Add(o)
    # merge and treat parameters and input
    sft_rundir=self.sft.rundir
    self.sft.SetRunDir(os.path.join(sft_rundir,'workingonly'))
    self.sft.MergeMTZ( *self.sft.inp.GetAll(), force_merge=True )
    for obj in self.sft.inp.GetAll(filetype='mtz'):
      excl=self.sft.out.AddCopy(obj,propagate=False)
      excl.AddCustomTag('workonly')
      self.sft.out.SetFileToChild( excl, 'workonly.mtz', 'mtz' )
    self.sft.SetKey( 'read', obj.GetFileName('mtz') )
    self.sft.SetKey( 'read', (exclude_sft.GetFileName('mtz'), 'col', "\"{0}\"".format(exclude_sft.GetLabel('free'))) )
    for obj in self.sft.inp.GetAll(filetype='mtz'):
      for c,l in obj.GetAllLabels():
        self.sft.SetKey( 'absent', ('col', l, 'if col', exclude_sft.GetLabel('free'), '= 0') )
    self.sft.SetKey( 'write', (excl.GetFileName('mtz'), '\nY') )
    self.sft.SetKey( 'quit\nY' )
    self.sft.Run(restore=False, clear_out=False)
    self.sft.SetRunDir(os.path.join(sft_rundir))
    self.dm.inp.Set(self.sft.out.fsigf)
    self.dm.inp.Set(self.sft.out.mapcoef)