import config import getData import loadData import fitData import numpy import segmentation import plotData import os import datetime import sys import json 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') print('Running {} for {}'.format(stage,config.getCode(r,setup))) 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])