Adding regularization to compartment fits

pythonScripts/config.py

@@ -104,8 +104,9 @@ def decode(code,xconfig):
    return {getPatientField(xconfig):fid,getVisitField(xconfig):vid}
 
 #standardized file names
-def getPattern(ftype,code,nclass=0,ir=0,ic=0,fitType='global',qaName='fits',timepoint=0,iseg=0,qLambda=0):
+def getPattern(ftype,code,nclass=0,ir=0,ic=0,fitType='global',qaName='fits',timepoint=0,iseg=0,qLambda=0,qLambdaC=0):
    qLambdaString='{:.3f}'.format(qLambda).replace('.','p')
+   qLambdaCString='{:.3f}'.format(qLambdaC).replace('.','p')
    replacePatterns={
       '_code_':code,
       '_nclass_':nclass,
@@ -115,7 +116,9 @@ def getPattern(ftype,code,nclass=0,ir=0,ic=0,fitType='global',qaName='fits',time
       '_qaName_':qaName,
       '_timepoint_':timepoint,
       '_iseg_':iseg,
-      '_qLambda_':qLambdaString}
+      '_qLambda_':qLambdaString,
+      '_qLambdaC_':qLambdaCString}
+
    map={'CT':'_code__CT.nrrd',
       'SPECT':'_code__Volume_timepoint_.nrrd',
       'Dummy':'_code__Dummy.mcsv',
@@ -129,8 +132,8 @@ def getPattern(ftype,code,nclass=0,ir=0,ic=0,fitType='global',qaName='fits',time
       'pixelFitPNG':'_code___nclass___ir__Pixel__fitType__centers_ic_.png',
       'fitIVF':'_code___nclass___ir___qLambda__fitParIVF.txt',
       'plotIVF':'_code___nclass___ir___qLambda__plotIVF__qaName_.png',
-      'fitCompartment':'_code___nclass___ir___qLambda__fitCompartment__iseg___qaName_.txt',
-      'plotCompartment':'_code___nclass___ir___qLambda__plotCompartment__iseg___qaName_.png'}
+      'fitCompartment':'_code___nclass___ir___qLambda__fitCompartment__iseg___qaName___qLambdaC_.txt',
+      'plotCompartment':'_code___nclass___ir___qLambda__plotCompartment__iseg___qaName___qLambdaC_.png'}
    w=map[ftype]
    for x in replacePatterns:
       w=re.sub(x,'{}'.format(replacePatterns[x]),w)

pythonScripts/fitData.py

@@ -97,21 +97,21 @@ def generateIVFFinite():
 def generateCModel():
     #generate approx candidate
     ig=initialValueGenerator()
-    ig.add('loggaus',[-3,1],[0,numpy.inf])
+    ig.add('loggaus',[-3,1],[0,0.1])
     ig.add('flat',[0,1],[0,1])
-    ig.add('loggaus',[-7,2],[0,numpy.inf])
-    ig.add('gaus',[10,5],[0,numpy.inf])
+    ig.add('loggaus',[-7,2],[0,0.1])
+    ig.add('gaus',[10,5],[0,30])
     return ig
 
 def generateCModelFinite():
     #generate approx candidate
     ig=initialValueGenerator()
     #k1
-    ig.add('loggaus',[-3,1],[0,1])
+    ig.add('loggaus',[-3,1],[0,0.1])
     #BVF
     ig.add('flat',[0,1],[0,1])
     #k2
-    ig.add('loggaus',[-7,2],[0,1])
+    ig.add('loggaus',[-7,2],[0,0.1])
     #dt
     ig.add('gaus',[10,2],[5,50])
     return ig
@@ -286,7 +286,7 @@ def fitCompartment(ivfFit, t, qf ,nfit=10,nbatch=20, useJac=False):
    return samples
 
 
-def fitCompartmentGlobal(ivfFit, t, qf ,nfit=10,nbatch=20, useJac=False):
+def fitCompartmentGlobal(ivfFit, t, qf ,nfit=10,nbatch=20, useJac=False,qLambda=0):
     #nclass=setup['nclass'][0]
     #t,dt=loadData.loadTime(r,setup)
     #centers=loadData.loadCenters(r,setup,ir)
@@ -317,11 +317,15 @@ def fitCompartmentGlobal(ivfFit, t, qf ,nfit=10,nbatch=20, useJac=False):
         w=numpy.ones(t.shape)
         fc1=functools.partial(fitModel.fCompartment,ivfPar)
         funScalar=functools.partial(fitModel.fDiffScalar,fc1,t,qf,w)  
+        funScalarRegularized=functools.partial(fitModel.fDiffScalarRegularized,funScalar,fitModel.fCRegBVF,qLambda)
         jac=functools.partial(fitModel.jacDep,ivfPar,t)
         jacScalar=functools.partial(fitModel.jacScalar,fc1,t,qf,w,jac)
-        minSetup=dict(method='L-BFGS-B',jac=jacScalar)
+        jacScalarRegularized=functools.partial(fitModel.jacScalarRegularized,jacScalar,fitModel.jacDepRegBVF,qLambda)
+
+        #minSetup=dict(method='L-BFGS-B',jac=jacScalar)
+        minSetup=dict(method='L-BFGS-B',jac=jacScalarRegularized)
         
-        result=scipy.optimize.dual_annealing(func=funScalar,bounds=boundsScalar,minimizer_kwargs=minSetup)
+        result=scipy.optimize.dual_annealing(func=funScalarRegularized,bounds=boundsScalar,minimizer_kwargs=minSetup)
 
         print(f'[{j}/{nfit}]')
         
@@ -329,7 +333,7 @@ def fitCompartmentGlobal(ivfFit, t, qf ,nfit=10,nbatch=20, useJac=False):
         qx=result.x
         qx[0]/=scale
         qx[1]/=scale
-        print(qx)
+        print('{} {} {}'.format(scale,result.fun/scale,qx))
         samples[1:nC+1,j]=qx
         samples[(1+nC):,j]=ivfPar
         

pythonScripts/fitModel.py

@@ -321,3 +321,14 @@ def jacDep(ivfPar,t,par):
 
    return jac
 
+def fCRegBVF(par):
+   return par[1]
+
+def jacDepRegBVF(par):
+   jac=numpy.zeros(par.shape[0])
+   jac[1]=1
+   return jac
+
+  
+
+

pythonScripts/loadData.py

@@ -184,7 +184,7 @@ def readIVF(r,xsetup,ir=0,qLambda=0):
    samples=fw[1:,:]
    return m,samples
 
-def saveSamples(r,xsetup,samples,m,name,iseg=0,ir=0):
+def saveSamples(r,xsetup,samples,m,name,iseg=0,ir=0,qLambda=0,qLambdaC=0):
    nclass=xsetup['nclass'][0]
    code=config.getCode(r,xsetup)
    fm=numpy.zeros(samples.shape[1])
@@ -192,15 +192,17 @@ def saveSamples(r,xsetup,samples,m,name,iseg=0,ir=0):
    for i in range(n):
       fm[i]=m[i]+1
    fw=numpy.vstack((fm,samples))
-   f=getData.getLocalPath(r,xsetup,config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg))
+   pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
+   f=getData.getLocalPath(r,xsetup,pat)
    print(f'Saving samples to {f}')
    numpy.savetxt(f,fw,delimiter=',')
 
-def readSamples(r,xsetup,name,iseg=0,ir=0):
+def readSamples(r,xsetup,name,iseg=0,ir=0,qLambda=0,qLambdaC=0):
    m=[]
    nclass=xsetup['nclass'][0]
    code=config.getCode(r,xsetup)
-   f=getData.getLocalPath(r,xsetup,config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg))
+   pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
+   f=getData.getLocalPath(r,xsetup,pat)
    print(f'Reading from  {f}')
    fw=numpy.loadtxt(f,delimiter=',')
    samples=fw[1:,:]
@@ -211,17 +213,19 @@ def readSamples(r,xsetup,name,iseg=0,ir=0):
       m.append(fm[i]-1)
    return m,samples
 
-def saveTAC(r,xsetup,tac,name,iseg=0,ir=0):
+def saveTAC(r,xsetup,tac,name,iseg=0,ir=0,qLambda=0,qLambdaC=0):
    nclass=xsetup['nclass'][0]
    code=config.getCode(r,xsetup)
-   f=getData.getLocalPath(r,xsetup,config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg))
+   pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
+   f=getData.getLocalPath(r,xsetup,pat)
    print(f'Saving samples to {f}')
    numpy.savetxt(f,tac,delimiter=',')
 
-def readTAC(r,xsetup,name,iseg=0,ir=0):
+def readTAC(r,xsetup,name,iseg=0,ir=0,qLambda=0,qLambdaC=0):
    nclass=xsetup['nclass'][0]
    code=config.getCode(r,xsetup)
-   f=getData.getLocalPath(r,xsetup,config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg))
+   pat=config.getPattern('fitCompartment',code=code,nclass=nclass,ir=ir,qaName=name,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC)
+   f=getData.getLocalPath(r,xsetup,pat)
    print(f'Reading from {f}')
    return numpy.loadtxt(f,delimiter=',')
 

pythonScripts/workflow.py

@@ -6,6 +6,7 @@ import numpy
 import segmentation
 import plotData
 import os
+import datetime
 
 def getRows(setup,returnFB=False):
     qFilter=config.getFilter(setup)
@@ -16,84 +17,90 @@ def getRows(setup,returnFB=False):
         return fb,db,rows
     return rows
 
-def listRequiredFiles(stage,r,setup):
-    code=config.getCode(r,setup)
-    nclass=setup['nclass'][0]
-    nr=setup['nr']
-    nt=20
-    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']=[]
-        qLambda=setup['qLambda']
-        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)
+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['segmentation']=[segmentation.getSegmentationFileName(r,setup)]
-        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
+      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']=[]
-        qLambda=setup['qLambda']
-        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)]
-        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':
-        names={}
-        names['center']=[]
-        names['fitIVF']=[]
-        names['fitCompartment']=[]
-        nseg=setup['nseg']
-        qLambda=setup['qLambda']
-        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):
-                rel=[config.getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=qLambda) for qn in sNames]
-                names['fitCompartment'].extend(rel)
-        names['segmentation']=[segmentation.getSegmentationFileName(r,setup)]
-        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
+   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 {}
+   return {}
 
 def listCreatedFiles(stage,r,setup):
     code=config.getCode(r,setup)
     nclass=setup['nclass'][0]
-    qLambda=setup['qLambda']
+    qLambda=getRobust(setup,'qLambda')
+    qLambdaC=getRobust(setup,'qLambdaC')
     nr=setup['nr']
-    try:
-       nseg=setup['nseg']
-    except KeyError:
-      nseg=0
+    nseg=getRobust(setup,'nseg')
     names={}
                                 
     if stage=='setCenters':
@@ -126,7 +133,7 @@ def listCreatedFiles(stage,r,setup):
         sNames=['kmeansFit','localFit','kmeansTAC','localTAC']
         for ir in range(nr):
             for iseg in range(nseg):
-                rel=[config.getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=qLambda) for qn in sNames]
+                rel=[config.getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) for qn in sNames]
                 names[xc].extend(rel)
     if stage=='plotCompartment':
         xc='plotCompartment'
@@ -134,16 +141,16 @@ def listCreatedFiles(stage,r,setup):
         sNames=['realizations','diff']
         for ir in range(nr):
             for iseg in range(nseg):
-                rel=[config.getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=qLambda) for qn in sNames]
+                rel=[config.getPattern(xc,code=code,nclass=nclass,ir=ir,qaName=qn,iseg=iseg,qLambda=qLambda,qLambdaC=qLambdaC) for qn in sNames]
                 names[xc].extend(rel)
     
     return []
     
 
-def getRequiredFiles(stage,r,setup,fb,names=None):
+def getRequiredFiles(stage,r,setup,fb,names=None,db=None):
     #fb,r=getRow(setup,True)
     if not names:
-        names=listRequiredFiles(stage,r,setup)
+        names=listRequiredFiles(stage,r,setup,db=db)
     for f in names:
         _copyFromServer=getData.copyFromServer
         if f=='segmentation':
@@ -151,10 +158,10 @@ def getRequiredFiles(stage,r,setup,fb,names=None):
         _copyFromServer(fb,r,setup,names[f])
     return fb,r
 
-def checkRequiredFiles(stage,r,setup,names=None,fb=None,doPrint=False):
+def checkRequiredFiles(stage,r,setup,names=None,fb=None,doPrint=False,db=None):
     ok=True
     if not names:
-        names=listRequiredFiles(stage,r,setup)
+        names=listRequiredFiles(stage,r,setup,db=db)
     for f in names:
         nm=names[f]
         for x in nm:
@@ -231,16 +238,14 @@ def workflow(r,setup,stage,fb=None,db=None):
     setIVF=False
     plotIVF=False
     setC=True
-    try:
-        qLambda=setup['qLambda']
-    except KeyError:
-        qLambda=0
+    qLambda=getRobust(setup,'qLambda')
+    qLambdaC=getRobust(setup,'qLambdaC')
                                    
     if stage=='setCenters':
-        names=listRequiredFiles(stage,r,setup)
+        names=listRequiredFiles(stage,r,setup,db=db)
         if fb:
-            getRequiredFiles(stage,r,setup,fb,names=names)
-        if not checkRequiredFiles(stage,r,setup,names=names,fb=fb,doPrint=True):
+            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)
@@ -248,7 +253,7 @@ def workflow(r,setup,stage,fb=None,db=None):
     if stage=='fitIVF':
         #get required files
         #stage='fitIVF'
-        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True):
+        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,db=db):
             return
         
         loadData.saveIVF(r,setup,nfit=30,qLambda=qLambda)
@@ -256,10 +261,10 @@ def workflow(r,setup,stage,fb=None,db=None):
     
     if stage=='plotIVF':
         ir=0
-        names=listRequiredFiles(stage,r,setup)
+        names=listRequiredFiles(stage,r,setup,db=db)
         if fb:
-            getRequiredFiles(stage,r,setup,fb,names=names)
-        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,names=names):
+            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')
@@ -295,10 +300,10 @@ def workflow(r,setup,stage,fb=None,db=None):
 
     if stage=='fitCompartment':
         ir=0
-        names=listRequiredFiles(stage,r,setup)
+        names=listRequiredFiles(stage,r,setup,db=db)
         if fb:
-            getRequiredFiles(stage,r,setup,fb,names=names)
-        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,names=names):
+            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
@@ -342,22 +347,22 @@ def workflow(r,setup,stage,fb=None,db=None):
                 s+=1
             qCenter/=s
             qData/=s
-            samplesC=fitData.fitCompartmentGlobal(ivfFit,t,qCenter,useJac=True,nfit=20)
-            samplesC1=fitData.fitCompartmentGlobal(ivfFit,t,qData,nfit=20,useJac=True)
+            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)
-            loadData.saveSamples(r,setup,samplesC1,[-1],'localFit',iseg=x,ir=ir,qLambda=qLambda)
-            loadData.saveTAC(r,setup,qCenter,'kmeansTAC',iseg=x,ir=ir,qLambda=qLambda)
-            loadData.saveTAC(r,setup,qData,'localTAC',iseg=x,ir=ir,qLambda=qLambda)
+            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)
+        names=listRequiredFiles(stage,r,setup,db=db)
         if fb:
-            getRequiredFiles(stage,r,setup,fb,names=names)
-        if not checkRequiredFiles(stage,r,setup,fb=fb,doPrint=True,names=names):
+            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'
@@ -369,10 +374,10 @@ def workflow(r,setup,stage,fb=None,db=None):
         setup['nseg']=len(segmentIds)
         t,dt=loadData.loadTime(r,setup)
         for iseg in segmentIds:
-            m,samplesC=loadData.readSamples(r,setup,'kmeansFit',ir=ir,iseg=iseg,qLambda=qLambda)
-            m1,samplesC1=loadData.readSamples(r,setup,'localFit',ir=ir,iseg=iseg,qLambda=qLambda)
-            qCenter=loadData.readTAC(r,setup,'kmeansTAC',ir=ir,iseg=iseg,qLambda=qLambda)
-            qData=loadData.readTAC(r,setup,'localTAC',ir=ir,iseg=iseg,qLambda=qLambda)
+            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,:]
@@ -391,13 +396,14 @@ def workflow(r,setup,stage,fb=None,db=None):
             row['mean']=k1
             row['std']=stdK1
             row['regionId']=iseg
-            row['fitPlot']=config.getPattern(tag,code=code,ir=0,nclass=nclass,qaName='realizations',iseg=iseg,qLambda=qLambda)
-            row['diffPlot']=config.getPattern(tag,code=code,ir=0,nclass=nclass,qaName='diff',iseg=iseg,qLambda=qLambda)
+            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
             row['qLambda']=qLambda
+            row['qLambdaC']=qLambdaC
             if db:
                 db.modifyRows('insert',setup['project'],'lists','Summary',[row,row1])
             
@@ -412,6 +418,8 @@ def workflow(r,setup,stage,fb=None,db=None):
 
 
     uploadCreatedFiles(stage,fb,r,setup)
+   #setup could be modified
+    return setup
 
 def main(setupFile):
    with open(setupFile,'r') as f: