dynamicSPECT/pythonScripts/workflow.py

import config
import getData
import loadData
import fitData
import numpy
import segmentation
import plotData
import os
import datetime

def getRows(setup,returnFB=False):
    qFilter=config.getFilter(setup)
    db,fb=getData.connectDB(setup['network'])
    rows=getData.getPatients(db,setup,qFilter)
    #print(rows)
    if returnFB:
        return fb,db,rows
    return rows

def getRobust(setup,var,default=0):
   #get from dict, return default if not set
   try:
      return setup[var]
   except KeyError:
      return default

def listRequiredFiles(stage,r,setup,db=None):
   code=config.getCode(r,setup)
   nclass=setup['nclass'][0]
   nr=setup['nr']
   nt=20
   qLambda=getRobust(setup,'qLambda')
   qLambdaC=getRobust(setup,'qLambdaC')

   if stage=='setCenters':
      names={x:[config.getPattern(x,code)] for x in ['CT','Dummy']}
      names['SPECT']=[config.getPattern('SPECT',code=code,timepoint=i) for i in range(nt)]
      return names
   if stage=='fitIVF':
      names={x:[config.getPattern(x,code)] for x in ['Dummy']}
      names['center']=[]
      for ir in range(nr):
         rel=[config.getPattern('center',code=code,nclass=nclass,ir=ir,ic=i) for i in range(nclass)]
         names['center'].extend(rel)
      return names
   if stage=='plotIVF':
      names={x:[config.getPattern(x,code)] for x in ['Dummy']}
      names['center']=[]
      names['fitIVF']=[]
      for ir in range(nr):
         rel=[config.getPattern('center',code=code,nclass=nclass,ir=ir,ic=i) for i in range(nclass)]
         names['center'].extend(rel)
         names['fitIVF'].append(config.getPattern('fitIVF',code=code,nclass=nclass,ir=ir,qLambda=qLambda))
         #names['center'].append(config.getPattern('centerMap',code=code,nclass=nclass,ir=ir))
         rel=[config.getPattern('centerNRRD',code=code,nclass=nclass,ir=ir,qaName=x) for x in ['CT','SPECT']]
         names['center'].extend(rel)
            
      names.update({x:[config.getPattern(x,code)] for x in ['CT']})
      names['SPECT']=[config.getPattern('SPECT',code=code,timepoint=i) for i in range(nt)]
      return names
    
   if stage=='fitCompartment':
      names={}
      names['center']=[]
      names['fitIVF']=[]
      for ir in range(nr):
         rel=[config.getPattern('center',code=code,nclass=nclass,ir=ir,ic=i) for i in range(nclass)]
         names['center'].extend(rel)
         names['fitIVF'].append(config.getPattern('fitIVF',code=code,nclass=nclass,ir=ir,qLambda=qLambda))
         names['center'].append(config.getPattern('centerMap',code=code,nclass=nclass,ir=ir))
      names['segmentation']=[segmentation.getSegmentationFileName(r,setup,db=db)]
      return names
    
   if stage=='plotCompartment':
      sNames=['kmeansFit','localFit','kmeansTAC','localTAC']
      names={}
      names['center']=[]
      names['fitIVF']=[]
      names['fitCompartment']=[]
      nseg=setup['nseg']
      for ir in range(nr):
         rel=[config.getPattern('center',code=code,nclass=nclass,ir=ir,ic=i) for i in range(nclass)]
         names['center'].extend(rel)
         names['fitIVF'].append(config.getPattern('fitIVF',code=code,nclass=nclass,ir=ir,qLambda=qLambda))
         names['center'].append(config.getPattern('centerMap',code=code,nclass=nclass,ir=ir))
         for iseg in range(nseg):
            xc='fitCompartment'
            rel=[config.getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) for qn in sNames]
            names['fitCompartment'].extend(rel)
      names['segmentation']=[segmentation.getSegmentationFileName(r,setup,db=db)]
      names.update({x:[config.getPattern(x,code)] for x in ['CT','Dummy']})
      names['SPECT']=[config.getPattern('SPECT',code=code,timepoint=i) for i in range(nt)]
      return names
        
   return {}

def listCreatedFiles(stage,r,setup):
   code=config.getCode(r,setup)
   nclass=setup['nclass'][0]
   qLambda=getRobust(setup,'qLambda')
   qLambdaC=getRobust(setup,'qLambdaC')
   nr=setup['nr']
   nseg=getRobust(setup,'nseg')
   print(f'List created files nseg={nseg}')
   names={}
   getPattern=config.getPattern
                                
   if stage=='setCenters':
      names['center']=[]
      for ir in range(nr):
         rel=[getPattern('center',code=code,nclass=nclass,ir=ir,ic=i) for i in range(nclass)]
         names['center'].extend(rel)
         #rel=[config.getPattern('centerWeight',code=code,nclass=nclass,ir=ir,ic=i) for i in range(nclass)]
         #names['center'].extend(rel)
         names['center'].append(getPattern('centerMap',code=code,nclass=nclass,ir=ir))
         rel=[getPattern('centerNRRD',code=code,nclass=nclass,ir=ir,qaName=x) for x in ['CT','SPECT']]
         names['center'].extend(rel)
            
      return names
   if stage=='fitIVF':
      names['fitIVF']=[]
      for ir in range(nr):
         names['fitIVF'].append(getPattern('fitIVF',code=code,nclass=nclass,ir=ir,qLambda=qLambda))
      return names


   if stage=='plotIVF':
      names['plotIVF']=[]
      sNames=['fits','diff','generatedIVF','centerIVFSPECT','centerIVFCT']
      for ir in range(nr):
         x=[getPattern('plotIVF',code=code,nclass=nclass,ir=ir,qaName=y,qLambda=qLambda) for y in sNames]
         names['plotIVF'].extend(x)
      return names


   if stage=='fitCompartment':
      xc='fitCompartment'
      names[xc]=[]
      sNames=['kmeansFit','localFit','kmeansTAC','localTAC']
      for ir in range(nr):
         for iseg in range(nseg):
            _q=qLambda
            _qC=qLambdaC
            rel=[getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=_q,qLambdaC=_qC) for qn in sNames]
            names[xc].extend(rel)
      return names

   if stage=='plotCompartment':
      xc=stage
      names[xc]=[]
      sNames=['realizations','diff']
      for ir in range(nr):
         for s in range(nseg):
            iseg=s+1
            _q=qLambda
            _qC=qLambdaC
            rel=[getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=_q,qLambdaC=_qC) for qn in sNames]
            names[xc].extend(rel)
      return names
    
   return []
    

def getRequiredFiles(stage,r,setup,fb,names=None,db=None):
    #fb,r=getRow(setup,True)
    if not names:
        names=listRequiredFiles(stage,r,setup,db=db)
    for f in names:
        _copyFromServer=getData.copyFromServer
        if f=='segmentation':
            _copyFromServer=segmentation.copyFromServer
        _copyFromServer(fb,r,setup,names[f])
    return fb,r

def checkRequiredFiles(stage,r,setup,names=None,fb=None,doPrint=False,db=None):
    ok=True
    if not names:
        names=listRequiredFiles(stage,r,setup,db=db)
    for f in names:
        nm=names[f]
        for x in nm:
            avail=os.path.isfile(getData.getLocalPath(r,setup,x))
            if not avail:
                print(f'Missing {x}')
                if fb:
                    _getURL=getData.getURL
                    if f=='segmentation':
                        _getURL=segmentation.getURL
                    availRemote=fb.entryExists(_getURL(fb,r,setup,x))
                    print(f'Available remote: {availRemote}')
                ok=False
            if doPrint:
                print(f'[{avail}] {x}')
    return ok

def uploadCreatedFiles(stage,fb,r,setup,names=None):
   if not names:
      names=listCreatedFiles(stage,r,setup)
   for f in names:
      _copyToServer=getData.copyToServer
      _getURL=getData.getURL
      if f=='segmentation':
         _copyToServer=segmentation.copyToServer
         _getURL=segmentation.getURL
      _copyToServer(fb,r,setup,names[f])
      for x in names[f]:
         print('[{}] Uploaded {}'.format(fb.entryExists(_getURL(fb,r,setup,x)),x))


#this is a poor fit for workflow, but no better logical unit was found, so here it is
def makeMap(segs,kClass,tac):
    map={}
    vals=[(kClass[i],tac[:,i]) for i in range(len(kClass))]
    
    for (i,v) in zip(segs,vals):
        try:
            map[i].append(v)
        except KeyError:
            map[i]=[v]
    return map

#def getDataAtPixels(data,loc) replaced by loadData.getTACAtPixels(data,loc)

def updateDatabase(r,setup,stage,fb=None,db=None,categories=[]):
   #set database entry
   try:
      qLam=setup['qLambda']
   except KeyError:
      qLam=0

   nclass=setup['nclass'][0]
   code=config.getCode(r,setup) 
   if stage=='plotIVF':
      m,samples=loadData.readIVF(r,setup,qLambda=qLam)
      chi2=samples[0,:]
      threshold=numpy.median(chi2)
      fit=fitData.getFit(samples,threshold)
      row={x:r[x] for x in ['PatientId','visitCode']}
      row['nclass']=nclass
      row['mean']=fit.mu[0]
      row['std']=fit.cov[0,0]
      row['qLambda']=qLam

      fNames={x:config.getPattern('plotIVF',code=code,ir=0,nclass=nclass,qaName=x,qLambda=qLam) for x in categories}
      row.update(fNames)
      if db:
         db.modifyRows('insert',setup['project'],'lists','SummaryIVF',[row])


def workflow(r,setup,stage,fb=None,db=None):
    setCenters=False
    setIVF=False
    plotIVF=False
    setC=True
    qLambda=getRobust(setup,'qLambda')
    qLambdaC=getRobust(setup,'qLambdaC')
                                   
    if stage=='setCenters':
        names=listRequiredFiles(stage,r,setup,db=db)
        if fb:
            getRequiredFiles(stage,r,setup,fb,names=names,db=db)
        if not checkRequiredFiles(stage,r,setup,names=names,fb=fb,doPrint=True,db=db):
            return
        
        loadData.saveCenters(r,setup)

    if stage=='fitIVF':
        #get required files
        #stage='fitIVF'
        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,db=db):
            return
        
        loadData.saveIVF(r,setup,nfit=30,qLambda=qLambda)

    
    if stage=='plotIVF':
        ir=0
        names=listRequiredFiles(stage,r,setup,db=db)
        if fb:
            getRequiredFiles(stage,r,setup,fb,names=names,db=db)
        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,names=names,db=db):
            return
        
        print('Loading files to memory')
        m,samples=loadData.readIVF(r,setup,qLambda=qLambda)

        data=loadData.loadData(r,setup)
        ct=loadData.loadCT(r,setup)
        centerMapSPECT,centerMapCT=loadData.loadCenterMapNRRD(r,setup,ir=ir)

        t,dt=loadData.loadTime(r,setup)
        centers=loadData.loadCenters(r,setup,ir=ir)
        ivf=centers[m]
        chi2=samples[0,:]
        threshold=numpy.median(chi2)
        code=config.getCode(r,setup)
        ir=0
          
        categories=['fits','diff','generatedIVF','centerIVFSPECT','centerIVFCT']

        fNames={x:config.getPattern('plotIVF',code=code,ir=0,nclass=nclass,qaName=x,qLambda=qLambda) for x in categories}
        files={x:getData.getLocalPath(r,setup,fNames[x]) for x in fNames}
                
        
        plotData.plotIVF(t,ivf,samples,threshold,file0=files['fits'],file1=files['diff'])
        
        plotData.plotIVFRealizations(t,ivf,samples,threshold,file=files['generatedIVF'])
        
        #temporarily blocking center generation
        plotData.plotIVFCenter(centerMapSPECT,centerMapCT,m,data,ct,file0=files['centerIVFSPECT'],
                               file1=files['centerIVFCT'])

        updateDatabase(r,setup,stage,db=db,fb=fb,categories=categories)

    if stage=='fitCompartment':
        ir=0
        names=listRequiredFiles(stage,r,setup,db=db)
        if fb:
            getRequiredFiles(stage,r,setup,fb,names=names,db=db)
        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,names=names,db=db):
            return
        
        #load class classification
        u=loadData.loadCenterMap(r,setup)
        print(u.shape)
        
        #load segmentation
        seg=segmentation.getNRRDImage(r,setup,names)
        loc=numpy.nonzero(seg)
        vClass=[int(x) for x in u[loc]]
        segments=[int(x) for x in seg[loc]]
        print(segments)
        data=loadData.loadData(r,setup)
        tac=loadData.getTACAtPixels(data,loc)
        
        segMap=makeMap(segments,vClass,tac)
        #for x in segMap:
        #    print('{} {}'.format(x,segMap[x]))
        #return
                  
        m1,samples=loadData.readIVF(r,setup,qLambda=qLambda)
        chi2=samples[0,:]
        threshold=numpy.median(chi2)
        ivfFit=fitData.getFit(samples,threshold)
        t,dt=loadData.loadTime(r,setup)
        centers=loadData.loadCenters(r,setup)
        #save segmentation pixels
        setup['nseg']=len(segMap.keys())
        for x in segMap:
            mArray=segMap[x]
            qCenter=numpy.zeros(t.shape[0])
            qData=numpy.zeros(t.shape[0])
            s=0
            #average over contributions for each segmentation included in map
            kCenters=[]
            for m in mArray:
                #m is a tuple of classId and tac
                kCenters.append(m[0])
                qCenter+=centers[m[0]]
                qData+=m[1]
                s+=1
            qCenter/=s
            qData/=s
            samplesC=fitData.fitCompartmentGlobal(ivfFit,t,qCenter,useJac=True,nfit=20,qLambda=qLambdaC)
            samplesC1=fitData.fitCompartmentGlobal(ivfFit,t,qData,nfit=20,useJac=True,qLambda=qLambdaC)
            
            loadData.saveSamples(r,setup,samplesC,kCenters,'kmeansFit',iseg=x,ir=ir,qLambda=qLambda,qLambdaC=qLambdaC)
            loadData.saveSamples(r,setup,samplesC1,[-1],'localFit',iseg=x,ir=ir,qLambda=qLambda,qLambdaC=qLambdaC)
            loadData.saveTAC(r,setup,qCenter,'kmeansTAC',iseg=x,ir=ir,qLambda=qLambda,qLambdaC=qLambdaC)
            loadData.saveTAC(r,setup,qData,'localTAC',iseg=x,ir=ir,qLambda=qLambda,qLambdaC=qLambdaC)

            
            
    if stage=='plotCompartment':
        ir=0
        names=listRequiredFiles(stage,r,setup,db=db)
        if fb:
            getRequiredFiles(stage,r,setup,fb,names=names,db=db)
        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,names=names,db=db):
            return
        
        tag='plotCompartment'
        seg=segmentation.getNRRDImage(r,setup,names)
        loc=numpy.nonzero(seg)
        segmentIds=list(set([int(x) for x in seg[loc]]))
        nclass=setup['nclass'][0]
        code=config.getCode(r,setup)
        setup['nseg']=len(segmentIds)
        print('Setting setup[nseg]={}'.format(setup['nseg']))
        t,dt=loadData.loadTime(r,setup)
        for iseg in segmentIds:
            m,samplesC=loadData.readSamples(r,setup,'kmeansFit',ir=ir,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
            m1,samplesC1=loadData.readSamples(r,setup,'localFit',ir=ir,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
            qCenter=loadData.readTAC(r,setup,'kmeansTAC',ir=ir,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
            qData=loadData.readTAC(r,setup,'localTAC',ir=ir,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
            chi2C=samplesC[0,:]
            threshold=numpy.median(chi2C)
            chi2C1=samplesC1[0,:]
            threshold1=numpy.median(chi2C1)
            fit=fitData.getFit(samplesC,threshold)
            fit1=fitData.getFit(samplesC1,threshold1)
            k1=fit.mu[0]
            stdK1=fit.cov[0,0]
            k11=fit1.mu[0]
            stdK11=fit1.cov[0,0]
            #update database with entries
            row={x:r[x] for x in ['PatientId','visitCode']}
            row['Date']=datetime.datetime.now().isoformat()
            row['nclass']=nclass
            row['option']='kmeansFit'
            row['mean']=k1
            row['std']=stdK1
            row['regionId']=iseg
            row['qLambda']=qLambda
            row['qLambdaC']=qLambdaC
            row['fitPlot']=config.getPattern(tag,code=code,ir=0,nclass=nclass,qaName='realizations',\
               iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
            row['diffPlot']=config.getPattern(tag,code=code,ir=0,nclass=nclass,qaName='diff',\
               iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
            row1={x:row[x] for x in row}
            row1['option']='localFit'
            row1['mean']=k11
            row1['std']=stdK11
            if db:
                db.modifyRows('insert',setup['project'],'lists','Summary',[row,row1])
            
                            
            evalArray=[(samplesC,qCenter,'blue'),
                       (samplesC1,qData,'orange')]
            
            file0=getData.getLocalPath(r,setup,row['fitPlot'])
            file1=getData.getLocalPath(r,setup,row['diffPlot'])
            plotData.plotSamples(t,evalArray,file0=file0,file1=file1)


    nseg=getRobust(setup,'nseg')
    print(f'At upload nseg={nseg}')
    uploadCreatedFiles(stage,fb,r,setup)
   #setup could be modified
    return setup

def main(setupFile):
   with open(setupFile,'r') as f:
      setup=json.load(f)
   db,fb,rows=getRows(setup,returnFB=True)
   for r in rows:
      workflow(r,setup,setup['stage'],fb=fb,db=db)

if __name__=="__main__":
   main(sys.argv[1])