NatComm2020/In-vivo/Calibration and longitudinal assessment/LongitudinalAssessment.py

# -*- coding: utf-8 -*-
"""
Created on Sun Apr 28 15:50:58 2019

@author: aemdlabs
"""
import deepdish as dd
import matplotlib.pyplot as plt
from glob import glob
import os
import numpy as np
import matplotlib.colors as colors
import matplotlib as mpl
import csv
from matplotlib.colors import LinearSegmentedColormap


MCSMapI={'SE1':'Ch03',
         'SE2':'Ch05',
         'SE3':'Ch01',
         'SE4':'Ch02',
         'SE5':'Ch22',
         'SE6':'Ch06',
         'SE7':'Ch16',
         'SE8':'Ch37',
         'SE9':'Ch20',
         'SE10':'Ch10',
         'SE11':'Ch24',
         'SE12':'Ch08',
         'SE13':'Ch14',
         'SE14':'Ch04',
         'SE15':'Ch18',
         'SE16':'Ch33',
         'SE17':'Ch34',
         'SE18':'Ch60',
         'SE19':'Ch38',
         'SE20':'Ch64',
         'SE21':'Ch40',
         'SE22':'Ch56',
         'SE23':'Ch42',
         'SE24':'Ch70',
         'SE25':'Ch66',
         'SE26':'Ch65',
         'SE27':'Ch68',
         'SE28':'Ch67',
         'SE29':'Ch55',
         'SE30':'Ch62',
         'SE31':'Ch58',
         'SE32':'Ch69',
         'ME1':'Ch57',
         'ME2':'Ch61',
         'ME3':'Ch53',
         'ME4':'Ch63',
         'ME5':'Ch52',
         'ME6':'Ch41',
         'ME7':'Ch49',
         'ME8':'Ch51',
         'ME9':'Ch46',
         'ME10':'Ch45',
         'ME11':'Ch44',
         'ME12':'Ch39',
         'ME13':'Ch54',
         'ME14':'Ch43',
         'ME15':'Ch50',
         'ME16':'Ch47',
         'ME17':'Ch32',
         'ME18':'Ch27',
         'ME19':'Ch30',
         'ME20':'Ch29',
         'ME21':'Ch28',
         'ME22':'Ch25',
         'ME23':'Ch26',
         'ME24':'Ch07',
         'ME25':'Ch21',
         'ME26':'Ch11',
         'ME27':'Ch17',
         'ME28':'Ch15',
         'ME29':'Ch13',
         'ME30':'Ch31',
         'ME31':'Ch19',
         'ME32':'Ch09'}

                          #Col, Row  
MCSMapFacingDown={'Ch58':(0,1),
                  'Ch57':(0,2),
                  'Ch56':(0,3),
                  'Ch55':(0,4),
                  'Ch54':(0,5),
                  'Ch53':(0,6),
                  'Ch52':(0,7),
                  'Ch51':(0,8),
                  'Ch50':(0,9),
                  'Ch49':(0,10),
                  'Ch60':(1,0),
                  'Ch61':(1,1),
                  'Ch62':(1,2),
                  'Ch63':(1,3),
                  'Ch64':(1,4),
                  'Ch65':(1,5),
                  'Ch43':(1,6),
                  'Ch44':(1,7),
                  'Ch45':(1,8),
                  'Ch46':(1,9),
                  'Ch47':(1,10),
                  'Ch70':(2,0),
                  'Ch69':(2,1),
                  'Ch68':(2,2),
                  'Ch67':(2,3),
                  'Ch66':(2,4),
                  'Ch42':(2,5),
                  'Ch41':(2,6),
                  'Ch40':(2,7),
                  'Ch39':(2,8),
                  'Ch38':(2,9),
                  'Ch37':(2,10),
                  'Ch01':(3,0),
                  'Ch02':(3,1),
                  'Ch03':(3,2),
                  'Ch04':(3,3),
                  'Ch05':(3,4),
                  'Ch06':(3,5),
                  'Ch30':(3,6),
                  'Ch31':(3,7),
                  'Ch32':(3,8),
                  'Ch33':(3,9),
                  'Ch34':(3,10),
                  'Ch11':(4,0),
                  'Ch10':(4,1),
                  'Ch09':(4,2),
                  'Ch08':(4,3),
                  'Ch07':(4,4),
                  'Ch29':(4,5),
                  'Ch28':(4,6),
                  'Ch27':(4,7),
                  'Ch26':(4,8),
                  'Ch25':(4,9),
                  'Ch24':(4,10),
                  'Ch12':None,
                  'Ch59':None,
                  'Ch13':(5,1),
                  'Ch14':(5,2),
                  'Ch15':(5,3),
                  'Ch16':(5,4),
                  'Ch17':(5,5),
                  'Ch18':(5,6),
                  'Ch19':(5,7),
                  'Ch20':(5,8),
                  'Ch21':(5,9),
                  'Ch22':(5,10)}
def MeanStd(Data):
    Arr = np.zeros([len(Data.keys()),len(Data[Data.keys()[0]])])
    for iT,TrtName in enumerate(Data.keys()):
        Arr[iT,:] = Data[TrtName]
    
    return np.mean(Arr,0), np.std(Arr,0)

def MeanStdCNP(Data):
    MeanCNP = np.array([])
    StdCNP = np.array([])
    for iT,TrtName in enumerate(Data.keys()):
        Mean = np.mean(Data[TrtName])
        Std = np.std(Data[TrtName])
        MeanCNP = np.append(MeanCNP, Mean)
        StdCNP = np.append(StdCNP, Std)
    return MeanCNP, StdCNP


plt.close('all')

#%% DC charact plotting

plt.close('all')

Files = [
         'IV-charact/IV-Day1-2/20190731_033840 Current vs Vgs_5sStep.txt',
         'IV-charact/IV-Day1-2/20190801_031608 Current vs Vgs5s-afterRec6-highRange.txt',
         'IV-charact/IV-Day2/20190802_115002 Current vs Vgs-24h2-PostRec7.txt',
         'IV-charact/IV-Day4/20190803_031630 Current vs Vgs-PostRec12.txt',
         'IV-charact/IV-Day6/20190805_113528 Current vs Vgs-PreRec1.txt',
         'IV-charact/IV-Day12/20190811_090313 Current vs Vgs-PostRec1.txt',
         'IV-charact/IV-Day26/20190825_080628 Current vs Vgs.txt']

color = ['r','orange','g','b','purple','m','cyan']
cmap = mpl.cm.rainbow

legend= ['Day 1','Day 2','Day 3','Day 5','Day 7','Day 12','Day 28']
time = [1,2,3,4,6,12,26]
CNPDict = {}
GMDict = {}
GMforSN = {}
counter = -1
for iFile, File in enumerate(Files):
    Fin = open(File)
    
    plt.figure(1)    
    reader = csv.reader(Fin, delimiter='\t')
    lenVgs = 1000
    Vgs = np.zeros(lenVgs)
    Ids = np.zeros([lenVgs,8])
    GMforSN[iFile] = {}
    
    reader = csv.reader(Fin, delimiter='\t')
    for il, e in enumerate(reader):
        if il == 0:
            continue
        print(e)
        
        vg = e[0]
        vg = vg.replace(',','.')
        Vgs[il-1] = float(vg)
        for iids, ids in enumerate(e[1:]):
            ids = ids.replace(',','.')
            Ids[il-1,iids] = float(ids)
    
        
    CNP = np.array([])
    GM = np.array([])
    IdsTrim = {}
    for i in range(8):
        if i == 2:
            continue
        IdsTrim[i] = np.trim_zeros(Ids[:,i])
        if i == 0:
            Vgs = np.append(Vgs[0:2],np.trim_zeros(Vgs[2:]))
            plt.plot(Vgs,IdsTrim[i],color = cmap(iFile / float(len(Files))), label = legend[iFile])
            plt.legend()
        else:
            Vgs = np.append(Vgs[0:2],np.trim_zeros(Vgs[2:]))
            plt.plot(Vgs,IdsTrim[i],color =cmap(iFile / float(len(Files))))
        VgsFit = np.linspace(Vgs[0],Vgs[-1],100)
        fitVal=np.polyval(np.polyfit(Vgs,IdsTrim[i],9),VgsFit)
        
        cnp = VgsFit[np.where(fitVal==np.min(fitVal))]
        CNP = np.append(CNP,cnp)
    
        gm = np.max(abs(np.polyval(np.polyder(np.polyfit(Vgs,IdsTrim[i],9)),np.linspace(cnp-0.18, cnp,100))))*1e-2 #mS/V
        GMforSN[iFile][i] = gm
        GM = np.append(GM,gm)
    CNPDict[iFile] = CNP
    GMDict[iFile] = GM

#### Calculate Gm
    
    if iFile not in [0,2,4,5,6]:
        continue
    counter += 1

    
    plt.figure(2)
    MeanIds, StdIds = MeanStd(IdsTrim)
    plt.plot(Vgs,MeanIds, color =color[counter], label = legend[iFile])
    plt.fill_between(Vgs,MeanIds+StdIds,MeanIds-StdIds, color =color[counter], alpha=0.3)

plt.legend(bbox_to_anchor=(0.0, 0.36, 0.7, 0.0),fontsize=13)
plt.xlabel('V$_{gs}$ (V)',fontsize=15)
plt.ylabel('I$_{ds}$ ($\mu$A)',fontsize=15)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)

MeanCNP, StdCNP = MeanStdCNP(CNPDict)
CNPArray = np.zeros([len(Files),7])
GMArray = np.zeros([len(Files),7])
w = 0.1
width = lambda p, w: 10**(np.log10(p)+w/2.)-10**(np.log10(p)-w/2.)
for ikey,key in enumerate(CNPDict.keys()):

    CNPArray[ikey,:] = CNPDict[key].T
    GMArray[ikey,:] = GMDict[key].T

plt.figure()
plt.boxplot(CNPArray.T, positions = time, widths =  width(time,w))
plt.xlabel('time (days)')
plt.ylabel('CNP (V)')
width = lambda p, w: 10**(np.log10(p)+w/2.)-10**(np.log10(p)-w/2.)

plt.figure()
plt.boxplot(GMArray.T, positions = time, widths = width(time,w))
plt.xlabel('time(h)')
plt.ylabel('Gm (mS/V)')
plt.xscale('log')



#%% Noise Plotting
dates = [
         'Day1-2',
         'Day2',
         'Day4',
         'Day6',
         'Day12',
         'Day26'
         ]

cmap = mpl.cm.rainbow


DCch = ['ME5','ME7','ME29','ME31','SE5','SE7','SE29','SE31']

Directory = 'PSD-charact/{}/'.format(dates[0])
ImportNoise = dd.io.load(Directory+'Power200Hz5')
noise = np.ones((len(dates),len(ImportNoise)))*(-10000)

GmDCNoise = np.zeros((len(dates),len(DCch)))
GMDC = np.zeros((len(dates),len(DCch)))

days = [1,2,3,6,11,29]
NoiseMin = 2.5e-10
GmMin = 1.5e-9
DayNum = days

figN,axN =plt.subplots()

#axN = axS.twinx()
for iday, day in enumerate(dates):
    Directory = 'PSD-charact/{}/'.format(day)
    ImportNoise = dd.io.load(Directory+'Power200Hz5')
    
    print(day)
    for ikey, key in enumerate(ImportNoise.keys()):
        print(iday,key,ikey)
        
        if key in DCch:
            if DCch.index(key) == 2:
                continue
            
            GmDCNoise[iday,DCch.index(key)] = 1e6*np.sqrt(ImportNoise[key]*200*np.log(10/1)/np.sqrt(2))/(GMforSN[iday][DCch.index(key)]*1e-4)
            GMDC[iday,DCch.index(key)] = (GMforSN[iday][DCch.index(key)])

        Aparam = np.sqrt(ImportNoise[key]*200)
        if  Aparam >= NoiseMin:
            noise[iday,ikey] = Aparam
        else:
            noise[iday,ikey] = -100000
            print('{} damaged -Trt{}'.format(iday,key))
        
        
        
        figN
        w = 0.04
        width = lambda p, w: 10**(np.log10(p)+w/2.)-10**(np.log10(p)-w/2.)
        axN.semilogy(np.random.normal(iday+1, w, size=1)-6*w, Aparam, 'k*',markersize=2)
    axN.set_xlabel('Time(days)')
    axN.set_ylabel('Estimate Irms(A) @1Hz')


    
figN

w = 0.035
width = lambda p, w: 10**(np.log10(p)+w/2.)-10**(np.log10(p)-w/2.)
PropS = {}
PropS['markeredgecolor'] = 'r'

N = axN.boxplot(noise.T,widths = 4*w)#,positions = np.array(days)+1.7*width(days,w),widths=width(days,w))
axN.set_xticklabels(days,)
#axN.set_xscale("log")
axN.set_xlabel('time (days)')
axN.set_ylabel('Estimate I$_{ds-rms}$ (A) @1Hz')


plt.legend()

for element in ['boxes', 'whiskers', 'fliers', 'means', 'medians', 'caps']:
    plt.setp(N[element], color='k',markersize=2)

for patch in N['boxes']:
    patch.set(color='k') 




fig2,ax2 = plt.subplots()

w = 0.1
width = lambda p, w: 10**(np.log10(p)+w/2.)-10**(np.log10(p)-w/2.)

bp = ax2.boxplot(GmDCNoise.T)#,positions = days, widths=width(days,w))

ax2.set_xticklabels(days,)
for element in ['boxes', 'whiskers', 'fliers', 'means', 'medians', 'caps']:
    plt.setp(bp[element], color='k')

for patch in bp['boxes']:
    patch.set(color='k') 

plt.xlabel('time (days)')
plt.ylabel('V$_{gs-rms}$ ($\mu$V)')