{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "remoteSourcesURL https://git0.fmf.uni-lj.si/studen/nixSuite/raw/master/remoteResources/resources.json\n",
      "{'labkeyInterface': {'url': 'https://git0.fmf.uni-lj.si/studen/labkeyInterface/archive/master.zip', 'branch': 'master', 'modules': []}, 'irAEMM': {'url': 'https://git0.fmf.uni-lj.si/studen/iraemm/archive/master.zip', 'branch': 'master', 'modules': ['iraemmBrowser']}, 'SlicerLabkeyExtension': {'url': 'https://git0.fmf.uni-lj.si/studen/SlicerLabkeyExtension/archive/SlicerExtensionIndex.zip', 'branch': 'SlicerExtensionIndex', 'modules': ['labkeyBrowser']}, 'limfomiPET': {'url': 'https://git0.fmf.uni-lj.si/studen/limfomiPET/archive/master.zip', 'branch': 'master', 'modules': ['imageBrowser']}, 'parseConfig': {'url': 'https://git0.fmf.uni-lj.si/studen/parseConfig/archive/master.zip', 'branch': 'master', 'modules': []}, 'orthancInterface': {'url': 'https://git0.fmf.uni-lj.si/studen/orthancInterface/archive/master.zip', 'branch': 'master', 'modules': []}}\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import os\n",
    "import SimpleITK\n",
    "import numpy\n",
    "import matplotlib.pyplot\n",
    "import subprocess\n",
    "import json\n",
    "import re\n",
    "import config\n",
    "import getData\n",
    "\n",
    "sys.path.append(os.path.join(os.path.expanduser('~'),'software','src','nixSuite','wrapper'))\n",
    "import nixWrapper\n",
    "nixWrapper.loadLibrary('labkeyInterface')\n",
    "import labkeyInterface\n",
    "import labkeyFileBrowser\n",
    "import labkeyDatabaseBrowser\n",
    "\n",
    "def calculateCenters(db,setup,nclass,nr):\n",
    "    #path of scripts is ../scripts\n",
    "    #one up\n",
    "    baseDir=os.path.dirname(os.getcwd())\n",
    "    #scriptPath\n",
    "    scriptPath=os.path.join(baseDir,'scripts','generateCenters.sh')\n",
    "    rows=getData.getPatients(db,setup)\n",
    "    for r in rows:\n",
    "        #download \n",
    "        code=config.getCode(r,setup)\n",
    "        for nc in nclass:\n",
    "            print('centers: {}/{},{}'.format(code,nc,nr))\n",
    "            print(subprocess.run([scriptPath,code,str(nc),str(nr)], \\\n",
    "                             check=True, stdout=subprocess.PIPE).stdout)\n",
    "            \n",
    "def doAnalysis(db,setup,nclass,nr,mode):\n",
    "    baseDir=os.path.dirname(os.getcwd())#one up\n",
    "    #this is if IVF is inferred together with k1 fits\n",
    "    #runScript='doAnalysis.sh'\n",
    "    #this is if IVF is taken from previous cluster fits\n",
    "    #runScript='doAnalysisIVF.sh'\n",
    "    \n",
    "    \n",
    "    if mode=='global':\n",
    "        runScript='doAnalysis.sh'\n",
    "        analysisType=''\n",
    "    if mode=='IVF':\n",
    "        runScript='doAnalysisIVF.sh'\n",
    "        analysisType='IVF_'\n",
    "    \n",
    "    try:\n",
    "        print('Setting runScript to {}'.format(runScript))\n",
    "    except NameError:\n",
    "        print('Mode can be one of (global,IVF))')\n",
    "        return\n",
    "        \n",
    "    #either doAnalysis or doAnalysisIVF.sh\n",
    "    scriptPath=os.path.join(baseDir,'scripts',runScript)\n",
    "\n",
    "    rows=getData.getPatients(db,setup)\n",
    "    for r in rows:\n",
    "        dataPath=config.getLocalDir(r,setup)\n",
    "        code=config.getCode(r,setup)\n",
    "        for nc in nclass:\n",
    "            for rId in numpy.arange(nr):\n",
    "                print('{} [{} {}/{}]'.format(code,nc,rId+1,nr))\n",
    "                #this is a duplicate of path generated in fitCenters.m\n",
    "                fName=os.path.join(dataPath,'{}_{}_{}_{}fitParFinal.txt'.format(code,nc,rId+1,analysisType))\n",
    "                if os.path.isfile(fName):\n",
    "                    print('Skipping; {} available.'.format(fName))\n",
    "                    continue\n",
    "                subprocess.run([scriptPath,code,str(nc),str(rId+1)],\\\n",
    "                         check=True, stdout=subprocess.PIPE)\n",
    "\n",
    "def doPixelAnalysis(db,setup,sigma2,mode='IVF'):\n",
    "    \n",
    "    baseDir=os.path.dirname(os.getcwd())#one up\n",
    "                \n",
    "    #in global mode, IVF parameters are inferred together with fits to classes\n",
    "    #this is essentially repeat of the above, except that classes are taken as\n",
    "    #time-response curves from pixels in the sigma2-defined neighborhood of\n",
    "    #target pixels\n",
    "    if mode=='global':\n",
    "        runScript='doPixelAnalysis.sh'\n",
    "        analysisType=''\n",
    "    #in IVF mode, the parameters of input function are taken from the cluster fit \n",
    "    #(doAnalysis above, with mode=general). The rest is the same as for global mode\n",
    "    if mode=='IVF':\n",
    "        runScript='doPixelIVFAnalysis.sh'\n",
    "        analysisType='IVF_'\n",
    "   \n",
    "    #either doPixelAnalysis or doPixelAnalysisIVF.sh\n",
    "    scriptPath=os.path.join(baseDir,'scripts',runScript)\n",
    "\n",
    "\n",
    "    try:\n",
    "        print('Setting runScript to {}'.format(runScript))\n",
    "    except NameError:\n",
    "        print('Mode can be one of (global,IVF))')\n",
    "        return\n",
    "    \n",
    "    rows=getData.getPatients(db,setup)\n",
    "    for r in rows:\n",
    "        dataPath=config.getLocalDir(r,setup)\n",
    "        code=config.getCode(r,setup)\n",
    "        sName=os.path.join(dataPath,'{}_Segmentation.txt'.format(code))\n",
    "        sFile=os.path.join(dataPath,sName)\n",
    "        x=numpy.loadtxt(sFile)\n",
    "        nc=x.shape[0]\n",
    "        for s2 in sigma2:\n",
    "            sigmaCode='{:.2f}'.format(s2)\n",
    "            sigmaCode=re.sub('\\.','p',sigmaCode)\n",
    "            fName=os.path.join(dataPath,'{}_{}_{}_Pixel{}_fitParFinal.txt'.format(code,nc,sigmaCode,analysisType))\n",
    "            if os.path.isfile(fName):\n",
    "                print('Skipping; {} available.'.format(fName))\n",
    "                continue\n",
    "            subprocess.run([scriptPath,code,str(s2)], \\\n",
    "                             check=True, stdout=subprocess.PIPE)\n",
    "            \n",
    "            \n",
    "def getIWeights(r,setup,nclass,realizationId,ic):\n",
    "    locDir=config.getLocalDir(r,setup)\n",
    "    code=config.getCode(r,setup)\n",
    "    fname='{}_{}_{}_center{}_centerWeigth.nrrd'.\\\n",
    "            format(code,nclass,realizationId+1,ic+1)\n",
    "    uFile=os.path.join(locDir,fname)\n",
    "            \n",
    "    imU=SimpleITK.ReadImage(uFile)\n",
    "    return SimpleITK.GetArrayFromImage(imU)\n",
    "\n",
    "def getGaussianWeight(nU,pt,sigma2,na):\n",
    "    #gaussian weighted summation of surrounding pixels\n",
    "    \n",
    "    #find point closest to the target point\n",
    "    idx=[int(x) for x in pt]\n",
    "    #running offset        \n",
    "    fidx=numpy.zeros(3)\n",
    "    #half of the neighborhood\n",
    "    na2=0.5*(na-1)\n",
    "    \n",
    "    fs=0\n",
    "    fWeight=0\n",
    "    for i in numpy.arange(na):\n",
    "        fidx[0]=idx[0]+(i-na2)\n",
    "        for j in numpy.arange(na):\n",
    "            fidx[1]=idx[1]+(j-na2)\n",
    "            for k in numpy.arange(na):\n",
    "                fidx[2]=idx[2]+(k-na2)\n",
    "                fidx=[int(x) for x in fidx]\n",
    "                fd=numpy.array(fidx)-numpy.array(pt)\n",
    "                dist2=numpy.sum(fd*fd)\n",
    "                fw=numpy.exp(-0.5*dist2/sigma2);\n",
    "                fs+=fw\n",
    "                fWeight+=fw*nU[tuple(fidx)]\n",
    "                #print('\\t{}/{}: {}/{:.2f} {:.2g} {:.3g} {:.2g}'.format(idx,fidx,numpy.sqrt(dist2),fw,nU[tuple(fidx)],fs,fWeight))\n",
    "    fWeight/=fs\n",
    "    return fWeight\n",
    "            \n",
    "\n",
    "\n",
    "#gets weights by class for a particular realization and sigma2 averaging\n",
    "def getWeights(db,r,setup,nclass,realizationId,sigma2,na):\n",
    "    #for w1, classes are in 0 to nclass-1 space\n",
    "    #na is the size of the neighborhood\n",
    "    idFilter=config.getIdFilter(r,setup)\n",
    "    visitFilter=config.getVisitFilter(r,setup)\n",
    "    code=config.getCode(r,setup)\n",
    "    rows=getData.getSegmentation(db,setup,[idFilter,visitFilter])\n",
    "    pts={r['regionId']:[float(x) for x in [r['x'],r['y'],r['z']]] for r in rows}\n",
    "    \n",
    "    w={region:numpy.zeros(nclass) for region in pts}\n",
    "    na2=0.5*(na-1)\n",
    "    \n",
    "    for c in numpy.arange(nclass):\n",
    "        nU=getIWeights(r,setup,nclass,realizationId,c)\n",
    "        \n",
    "        for region in w:\n",
    "            #print(pts[region])\n",
    "            #print('{} {}/{} {}'.format(code,c,nclass,region))\n",
    "            #gaussian weighted summation of surrounding pixels\n",
    "            w[region][c]=getGaussianWeight(nU,pts[region],sigma2,na)\n",
    "            \n",
    "            \n",
    "    return w\n",
    "\n",
    "\n",
    "#gets fitPar for a particular realization in [0..nr-1] range\n",
    "def getFitPar(r,setup,nclass,realizationId,mode=''):\n",
    "    #fitGoodnes A tau alpha delay [k1 BVF k2 delay]xNcenters\n",
    "    allowedModes=['','IVF','Pixel','PixelIVF']\n",
    "    if mode not in allowedModes:\n",
    "        print('Mode should be one of {}'.format(allowedModes))\n",
    "        return []\n",
    "    \n",
    "    if mode=='PixelIVF':\n",
    "        #4th parameter is sigma, not realizationId\n",
    "        rCode='{:.2f}'.format(realizationId)\n",
    "        rCode=re.sub('\\.','p',rCode)\n",
    "    else:\n",
    "        #add one to match matlab 1..N array indexing\n",
    "        rCode='{}'.format(realizationId+1)\n",
    "    \n",
    "    if len(mode)>0:\n",
    "        mode=mode+'_'\n",
    "    \n",
    "    code=config.getCode(r,setup)\n",
    "    fname='{}_{}_{}_{}fitParFinal.txt'.format(code,nclass,rCode,mode)\n",
    "    locDir=config.getLocalDir(r,setup)\n",
    "    uFile=os.path.join(locDir,fname)\n",
    "    return numpy.genfromtxt(uFile,delimiter='\\t')\n",
    "\n",
    "def getK1(fitPar,iclass):\n",
    "    #fitGoodnes A tau alpha delay [k1 BVF k2 delay]xNcenters\n",
    "    #0 to nclass-1 space\n",
    "    return fitPar[4*iclass+5]\n",
    "\n",
    "def calculateK1(w,fitPar):\n",
    "    #calculateK1 for region weights\n",
    "    #return the k1 that belongs to the \n",
    "    #maximum class in region (M) and \n",
    "    #a weighted average (W)\n",
    "    k1={region:{'M':0,'W':0} for region in w}\n",
    "    for region in w:\n",
    "        cmax=numpy.argmax(w[region])\n",
    "        k1[region]['M']=getK1(fitPar,cmax)\n",
    "        fs=0\n",
    "        for c in numpy.arange(len(w[region])):\n",
    "            fs+=w[region][c]*getK1(fitPar,c)\n",
    "        k1[region]['W']=fs\n",
    "    return k1\n",
    "\n",
    "def getPatientValuesByNclass(db,r,setup,nclass,nrealizations,sigma2,na):\n",
    "    #summary script\n",
    "    #db is for database; needs segmentations\n",
    "    #r,setup identify patient\n",
    "    #nclass and nrealizations select strategy\n",
    "    #sigma2 is for combining output from adjacent pixels\n",
    "    #na is neighborhood size where smoothing/combination is done\n",
    "    k1={}\n",
    "    for rId in numpy.arange(nrealizations):\n",
    "        w=getWeights(db,r,setup,nclass,rId,sigma2,na)\n",
    "        fitPar=getFitPar(r,setup,nclass,rId,'IVF')\n",
    "        qk1=calculateK1(w,fitPar)\n",
    "        for region in w:\n",
    "            for type in qk1[region]:\n",
    "                try:\n",
    "                    k1[region][type].append(qk1[region][type])\n",
    "                except KeyError:\n",
    "                    k1={region:{type:[] for type in qk1[region]} for region in w}\n",
    "                    print(type)\n",
    "                    k1[region][type].append(qk1[region][type])\n",
    "        print('[{}] {}/{}'.format(nclass,rId+1,nrealizations))\n",
    "    return k1   \n",
    "\n",
    "def getSummaryPatientValuesByNclass(db,r,setup,nclass,nrealizations,sigma2,na):\n",
    "    #summary script, same arguments as above\n",
    "    #also return deviation over realization\n",
    "    k1=getPatientValuesByNclass(db,r,setup,nclass,nrealizations,sigma2,na)\n",
    "    avgType=['M','W']\n",
    "    summaryK1={type:{region:{\n",
    "        'mean':numpy.mean(k1[region][type]), \n",
    "        'std':numpy.std(k1[region][type]), \n",
    "        'median':numpy.median(k1[region][type])} for region in k1}\n",
    "               for type in avgType}\n",
    "    \n",
    "    return summaryK1\n",
    "\n",
    "def fullSummary(db,setup,classes,nr,sigma2,na):\n",
    "    rows=getData.getPatients(db,setup)\n",
    "    return \\\n",
    "        {config.getCode(r,setup):\\\n",
    "         {c:getSummaryPatientValuesByNclass(db,r,setup,c,nr,sigma2,na) for c in classes} for r in rows}\n",
    "          \n",
    "def storeSummary(db,setup,summary,sigma2,na):\n",
    "    #dsM=db.selectRows(project,'study','Imaging',[])\n",
    "    for rCode in summary:\n",
    "        r=config.decode(rCode,setup)\n",
    "        idFilter=config.getIdFilter(r,setup)\n",
    "        visitFilter=config.getVisitFilter(r,setup)\n",
    "        for c in summary[rCode]:\n",
    "            cFilter={'variable':'nclass','value':str(c),'oper':'eq'}\n",
    "            for t in summary[rCode][c]:\n",
    "                tFilter={'variable':'option','value':t,'oper':'eq'}\n",
    "                for rId in summary[rCode][c][t]:\n",
    "                    rFilter={'variable':'regionId','value':str(rId),'oper':'eq'}\n",
    "                    rows=getData.getSummary(db,setup,[idFilter,visitFilter,cFilter,tFilter,rFilter])\n",
    "                    if len(rows)>0:\n",
    "                        qrow=rows[0]\n",
    "                        mode='update'\n",
    "                    else:\n",
    "                        qrow={qr:r[qr] for qr in r}\n",
    "                        qrow['nclass']=c\n",
    "                        qrow['option']=t\n",
    "                        qrow['regionId']=rId\n",
    "                        seqNum=config.getTargetSeqNum(r,setup)\n",
    "                        qrow['SequenceNum']=100+seqNum+c+0.001*rId\n",
    "                        if t=='M':\n",
    "                            qrow['SequenceNum']+=0.0001\n",
    "                        mode='insert'\n",
    "                    for v in summary[rCode][c][t][rId]:\n",
    "                        qrow[v]=summary[rCode][c][t][rId][v]\n",
    "                    qrow['sigma2']=sigma2\n",
    "                    qrow['na']=na\n",
    "                    getData.updateSummary(db,setup,mode,[qrow])\n",
    "                    \n",
    "def summaryPixelIVF(db,setup,sigma2):\n",
    "    #for second type of analysis (pixel based regions)\n",
    "    rows=getData.getPatients(db,setup)\n",
    "    return \\\n",
    "        {config.getCode(r,setup):\\\n",
    "         {s2:getPixelIVF(db,r,setup,s2) for s2 in sigma2} for r in rows}\n",
    "\n",
    "    \n",
    "def storeIVF(db,setup,summary):\n",
    "    for rCode in summary:\n",
    "        r=config.decode(rCode,setup)\n",
    "        idFilter=config.getIdFilter(r,setup)\n",
    "        visitFilter=config.getVisitFilter(r,setup)\n",
    "        for s2 in summary[rCode]:\n",
    "            sigmaFilter={'variable':'sigma2','value':str(s2),'oper':'eq'}\n",
    "            nr=len(summary[rCode][s2])\n",
    "            for rId in summary[rCode][s2]:\n",
    "                rFilter={'variable':'regionId','value':str(rId),'oper':'eq'}\n",
    "                typeFilter={'variable':'option','value':'D','oper':'eq'}\n",
    "                rows=getData.getSummary(db,setup,[idFilter,visitFilter,sigmaFilter,rFilter,typeFilter])\n",
    "                if len(rows)>0:\n",
    "                    qrow=rows[0]\n",
    "                    mode='update'\n",
    "                else:\n",
    "                    qrow={qr:r[qr] for qr in r}\n",
    "                    qrow['sigma2']=s2\n",
    "                    qrow['regionId']=rId\n",
    "                    seqNum=config.getTargetSeqNum(r,setup)\n",
    "                    qrow['SequenceNum']=140+seqNum+0.01*rId+0.001*s2\n",
    "                    mode='insert'\n",
    "                qrow['mean']=summary[rCode][s2][rId]\n",
    "                qrow['na']=7\n",
    "                qrow['nclass']=nr\n",
    "                qrow['option']='D'\n",
    "                getData.updateSummary(db,setup,mode,[qrow])\n",
    "                \n",
    "def getPixelIVF(db,r,setup,sigma2):\n",
    "    idFilter=config.getIdFilter(r,setup)\n",
    "    visitFilter=config.getVisitFilter(r,setup)\n",
    "    rows=getData.getSegmentation(db,setup,[idFilter,visitFilter])\n",
    "    nclassIVF=len(rows)\n",
    "    fitPar=getFitPar(r,setup,nclassIVF,sigma2,'PixelIVF')\n",
    "    k1={r['regionId']:getK1(fitPar,r['regionId']) for r in rows}\n",
    "    return k1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#setup\n",
    "fsetup='../template/cardiacSPECT.json'\n",
    "with open(fsetup,'r') as f:\n",
    "    setup=json.load(f)\n",
    "\n",
    "net=labkeyInterface.labkeyInterface()\n",
    "fconfig=os.path.join(os.path.expanduser('~'),'.labkey',setup['network'])\n",
    "net.init(fconfig)\n",
    "#net.getCSRF()\n",
    "db=labkeyDatabaseBrowser.labkeyDB(net)\n",
    "fb=labkeyFileBrowser.labkeyFileBrowser(net)\n",
    "    \n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#doAnalysis\n",
    "\n",
    "#0.) convert to NRRD \n",
    "# ~/software/src/SlicerLabkeyExtenstion/slicerScripts/runSlicer.sh convertToNRRD.py ../template/cardiacSPECT.json)\n",
    "\n",
    "#1.) downloadFiles\n",
    "\n",
    "#getData.downloadPatientFiles(db,fb,setup)\n",
    "\n",
    "#2.) generate (c-cluster) centers\n",
    "\n",
    "#setup the analysis (class sizes, number of realizations)\n",
    "nclass=[10,20,30]\n",
    "#nclass=[10]\n",
    "nr=20\n",
    "#nr=1\n",
    "\n",
    "#takes multiple hours for nr=20, nclass=[10,20,30]\n",
    "#runs matlab - maybe the first time licensing will be checked and error will be returned\n",
    "#just run once\n",
    "#calculateCenters(db,setup,nclass,nr)\n",
    "\n",
    "#3.) now do analysis - global gets the average IVF fit (see point one of the slides for both branches of analysis)\n",
    "#this again takes multiple hours for four images\n",
    "#from what I see this could be shortend (since only global parameters are sought)\n",
    "#doAnalysis(db,setup,nclass,nr,'global')\n",
    "\n",
    "#4.) do cluster analysis (bullets two and three under ML analysis)\n",
    "#doAnalysis(db,setup,nclass,nr,'IVF')\n",
    "\n",
    "#5.) sort out the segmentation points -> segmentation.ipynb and generate Segmentation.txt files read by pixelAnalysis\n",
    "\n",
    "#6.) run doPixelAnalysisIVF to calculate parameters from locations\n",
    "sigma2=[0.1,1,4]\n",
    "#doPixelAnalysis(db,setup,sigma2)\n",
    "\n",
    "#7.) store summary to database\n",
    "na=3\n",
    "s2=0.1\n",
    "#summary=fullSummary(db,setup,nclass,nr,s2,na)\n",
    "#storeSummary(db,setup,summary,s2,na)\n",
    "\n",
    "#summaryIVF=summaryPixelIVF(db,setup,sigma2)\n",
    "#storeIVF(db,setup,summaryIVF)\n",
    "\n",
    "#8.) Plot regions through chart wizard on labkey\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}