File contents

def makeHsqcPeakListFromShifts(argServer):
  
  from ccpnmr.analysis.UnitConverter   import ppm2pnt, pnt2ppm
  from ccpnmr.analysis.AssignmentBasic import assignResToDim
  from ccpnmr.analysis.MoleculeBasic   import areResonancesBound
  from ccpnmr.analysis.PeakBasic       import pickPeak
  from ccpnmr.analysis.Util            import getDataDimAssignmentTolerance
  
  hsqcNames = ['H[N]','H[C]']

  project = argServer.getProject()

  argServer.showInfo('Select shift list')

  shiftList = argServer.getShiftList()

  assignedOnly = argServer.askYesNo('Use assigned resonances only?')
  
  spectra = []
  for experiment in project.nmrExperiments:
    refExperiment = experiment.refExperiment
    if refExperiment and (refExperiment.name in hsqcNames):
      for spec in experiment.dataSources:
        spectra.append(spec)
  
  if not spectra:
    return
  
  if len(spectra) > 1:
    argServer.showInfo('Select HSQC spectrum to associate with new peak list')
    spectrum = argServer.getSpectrum(spectra)

  else:
    spectrum = spectra[0]

  if shiftList and spectrum:
  
    N = len(spectrum.dataDims)
    isotopes    = []
    fullRegion  = []
    for dataDim in spectrum.dataDims:
      dataDimRef = dataDim.dataDimRefs[0]
      expDimRef  = dataDimRef.expDimRef
      valueMin   = expDimRef.minAliasedFreq or pnt2ppm(1,dataDimRef)
      valueMax   = expDimRef.maxAliasedFreq or pnt2ppm(dataDim.numPoints,dataDimRef)

      fullRegion.append([valueMin,valueMax])
      isotopes.append( list(expDimRef.isotopeCodes) )
  
    dimShifts = []
    for i in range(N):
      dimShifts.append([])
    
    for shift in shiftList.measurements:
      ppm         = shift.value
      resonance   = shift.resonance
      isotopeCode = resonance.isotopeCode
      
      for i in range(N):
        if isotopeCode in isotopes[i]:
          if (ppm > fullRegion[i][0]) and (ppm < fullRegion[i][1]):
            if assignedOnly:
              if resonance.resonanceSet:
                dimShifts[i].append(shift)
                
            else:
              dimShifts[i].append(shift)


     
    def getShiftIntersections(shiftsList, intersections=None, shifts=None, dim=0):
      
      numDim = len(shiftsList)
      
      if shifts is None:
        shifts = []

      if intersections is None:
        intersections = []
      
      for shift in shiftsList[dim]:
        shifts0 = list(shifts) # copy
        if (not shifts0) or ( areResonancesBound(shift.resonance,shifts0[-1].resonance) ):

          shifts0.append(shift)
 
          if (dim+1) < len(shiftsList):
 
            getShiftIntersections(shiftsList, intersections, shifts0, dim+1)
 
          elif len(shifts0) == numDim: # last dimension

            intersections.append(shifts0)
      
      return intersections
      

    intersections = getShiftIntersections(dimShifts)

    tolerances = []
    for dataDim in spectrum.dataDims:
      tolerances.append( getDataDimAssignmentTolerance(dataDim) )

    i = 0
    I = len(intersections)
    if I:
      peakList = spectrum.newPeakList(isSimulated=True)

    print "Making %d peaks" % I
    
    for intersection in intersections:
      if i and i % 10 == 0:
        print "  Done %d of %d" % (i, I)
        
      position = []
      for j in range(N):
        ppm = intersection[j].value
        position.append(ppm)

      peak = pickPeak(peakList, position, unit=shiftList.unit)
      for j in range(N):
        assignResToDim(peak.peakDims[j],intersection[j].resonance)
      
      i += 1
 
    print "Done"