NatComm2020/In-vivo/Calibration and longitudinal assessment/Day1-2/Calibration scripts/CalcUd0.py

# -*- coding: utf-8 -*-
"""
Created on Sat Oct 13 18:50:52 2018

@author: aemdlabs
"""

from PhyREC.NeoInterface import NeoSegment, NeoSignal#, ReadMCSFile
import PhyREC.SignalAnalysis as Ran
import PhyREC.PlotWaves as Rplt
import quantities as pq
import matplotlib.pyplot as plt
import numpy as np
import neo
import deepdish as dd
import csv
from datetime import datetime
import PyGFET.DataClass as Gfetdat


def ReadMCSFile(McsFile, OutSeg=None, SigNamePrefix=''):
    import McsPy.McsData as McsData

    Dat = McsData.RawData(McsFile)
    Rec = Dat.recordings[0]
    NSamps = Rec.duration

    if OutSeg is None:
        OutSeg = NeoSegment()

    for AnaStrn, AnaStr in Rec.analog_streams.iteritems():
        if len(AnaStr.channel_infos) == 1:
            continue

        for Chn, Chinfo in AnaStr.channel_infos.iteritems():
            print('Analog Stream ', Chinfo.label, Chinfo.sampling_frequency)
            ChName = str(SigNamePrefix + Chinfo.label)
            print(ChName)

            Fs = Chinfo.sampling_frequency
            Var, Unit = AnaStr.get_channel_in_range(Chn, 0, NSamps)
            sig = neo.AnalogSignal(pq.Quantity(Var, Chinfo.info['Unit']),
                                   t_start=0*pq.s,
                                   sampling_rate=Fs.magnitude*pq.Hz,
                                   name=ChName)

            OutSeg.AddSignal(sig)
    return OutSeg

def ReadLogFile(File):
    Fin = open(File)
    
    reader = csv.reader(Fin, delimiter='\t')
    
    LogVals = {}
    ValsPos = {}
    for il, e in enumerate(reader):
        if il == 0:
            for ih, v in enumerate(e):
                ValsPos[ih] = v
                LogVals[v] = []
        else:
            for ih, v in enumerate(e):
                par = ValsPos[ih]
                if (par=='Vgs') or (par=='Vds') or (par=='Vref'):
                    LogVals[par].append(float(v.replace(',','.')))
                elif par == 'Date/Time':
                    LogVals[par].append(datetime.strptime(v, '%d.%m.%Y %H:%M:%S'))
                else:
                    LogVals[par].append(v)
    
    deltas = np.array(LogVals['Date/Time'])[:]-LogVals['Date/Time'][0]
    LogVals['Time'] = []
    for d in deltas:
        LogVals['Time'].append(d.total_seconds())
        
    Fin.close()
    
    return LogVals

def GetSwitchTimes(Sig, Thres=-1e-4, Plot=True):
    s = Sig.GetSignal(None)
    ds = np.abs(np.diff(np.array(s), axis=0))

    if Plot:
        plt.figure()
        plt.plot(s.times, s)
        plt.plot(s.times[1:], ds)

    ds = Sig.duplicate_with_new_array(signal=ds)
    Times = Ran.threshold_detection(ds,
                                    threshold=Thres,
                                    RelaxTime=5*pq.s)
    return Times

def CalcVgeffAC(Sig, TACchar, Vgs, vgs):
    i = list(np.around(Vgs,3)).index(vgs)
    gm = TACchar[i]*pq.A/pq.V
    return neo.AnalogSignal(Sig/gm,
                            units='V',
                            t_start=Sig.t_start,
                            sampling_rate=Sig.sampling_rate,
                            name=str(Sig.name),
                            file_origin=Sig.file_origin)
    
    
if __name__ == '__main__':
    #### Test autosweep


    Path = '../../Day0/IV charact/'
    PBSInFileM = Path + '2019-07-30T15-39-02B12784O18-T3-ACDC-PreImplant-NoMemb.h5' ############
    PBSInFileS = Path + '2019-07-30T15-39-02B12784O18-T3-ACDC-PreImplant-NoMemb_2.h5' ##########
    PBSLogFile = Path + 'B12784O18-T3-ACDC-PBS-PreImplant-20s-NoMemb.txt'
    StartCycle = 0

    Directory = '../IV charact/'
    File = '20190731_033840 Current vs Vgs_5sStep.txt'
    Fin = open(Directory+File)
     

    LogVals = ReadLogFile(PBSLogFile)
    delta = np.mean([t-LogVals['Time'][it] for it, t in enumerate(LogVals['Time'][1:])])*pq.s
    Delay = delta * StartCycle

    DCch = ('ME5', 'ME7', 'ME29', 'ME31', 'SE5', 'SE7', 'SE29', 'SE31')

    
    TrigChannel = 'ME29'
    TrigThres = 1e-5
    Vgs = np.array(LogVals['Vgs']) 
    VgsInt = np.linspace(Vgs[0],Vgs[-1],len(Vgs)*10)
    Vds = [LogVals['Vds'][0]]
    
    
    ivgain1 = 12e3*pq.V  
    ivgain2 = 101
    ACgain = 10*1/pq.V
    DCgain = 1*pq.V
    Fsig = 10
    StabTime = 2*pq.s
    GuardTime = 1*pq.s
    BW = 100
    ivgainDC = 118.8*pq.V 
    ivgainAC = 1188*pq.V 


    Rec = ReadMCSFile(PBSInFileM,
                      OutSeg=None,
                      SigNamePrefix='M')

    Rec = ReadMCSFile(PBSInFileS,
                      OutSeg=Rec,
                      SigNamePrefix='S')


# %% In vivo Ids
    plt.close('all')  
    ChKo = {}
    reader = csv.reader(Fin, delimiter='\t')
    lenVgs = 1000
    Vgs = np.zeros(lenVgs)
    Ids = np.zeros([lenVgs,8])
    
    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)
    
        
    color = ['m','r','orange','y','g','c','b','violet']
    plt.figure()
    Ud0 = []
    for i in range(8):
        Vgs = np.append(Vgs[0:2],np.trim_zeros(Vgs[2:]))
        ids = np.trim_zeros(Ids[:,i]) 
        plt.plot(Vgs,ids,color = color[i])
        if np.any(ids>=5.0):
            
            index = np.where(ids == np.amin(ids))
            print(ids, np.amin(ids), index)
            Ud0.append(Vgs[index[0][0]])
    Ud0meanInVivo = np.mean(np.array(Ud0))
#%%
    
    SwTimes = GetSwitchTimes(Sig=Rec.GetSignal(TrigChannel),
                             Thres=TrigThres,
                             Plot=True)
    SlotsDC = []
    SlotsAC = []

    for sig in Rec.Signals():
        sig.__class__ = NeoSignal

        if sig.name  in DCch:
            if sig.name.startswith('M'):
                col = 'r'
            else:
                col = 'g'
            SlotsDC.append(Rplt.WaveSlot(sig,
                                         Position=0,
                                         Color=col,
                                         Alpha=0.5))
        if sig.name.startswith('M'):
            col = 'k'
        else:
            col = 'b'
        SlotsAC.append(Rplt.WaveSlot(sig,
                                   Position=1,
                                   Color=col,
                                   Alpha=0.5))
            
    Splots = Rplt.PlotSlots(SlotsDC)
    Splots.PlotChannels(Time=None,
                        Units='mV')
    Splots.PlotEvents(SwTimes,
                      color='k')
    Splots.PlotEvents(LogVals['Time']*pq.s+SwTimes[0]+Delay,
                      color='r')
    
    SwTimes = LogVals['Time']*pq.s+SwTimes[0]+Delay
    figId, Axiv = plt.subplots()
    fig, Axt = plt.subplots()
    Ids = {}
    for sl in SlotsDC:
        Ids[sl.name] = []
        for isw, (t, vg) in enumerate(zip(SwTimes, Vgs)):
            if isw == len(SwTimes)-1:
                ts = sl.Signal.t_stop
            else:
                ts = SwTimes[isw+1]
            
            TWind = (t+StabTime, ts-GuardTime)
            s = sl.GetSignal(TWind, Units='V')
            Axt.plot(s)
            vio = np.mean(s).magnitude
            ids = (vio*101-(-vg+Vds))/12e3
            if ids <= 7e-6:
                ChKo[sl.name] = True
            else:
                ChKo[sl.name] = False
                
            Ids[sl.name].append(ids)
        if ChKo[sl.name] == False:
            Axiv.plot(Vgs,Ids[sl.name], label=sl.name)
   
    DevACvals = {}
    Ud0mean = 0
    dev = Ids
    Ud0 = np.zeros(len(dev))
    for ich, (ch, ids) in enumerate(dev.iteritems()):
        
        if ChKo[ch] == False:
            print(ich)
            DevACvals[ch] = { 'Ids': np.array(ids[:-1]),
                              'Vgs': Vgs[:-1],
                              'Vds': Vds}
            DataDC = Gfetdat.DataCharDC(DevACvals[ch])
            Ud0[ich] = DataDC.GetUd0()
            print(Ud0[ich])
            Ud0mean=Ud0mean+Ud0[ich]
            Axiv.plot(VgsInt, DataDC.GetIds(VgsInt),'--')
        Ud0meanPBS = Ud0mean/(ich+1)




#%%  Calc GM
    GM = {}
    CalFile = {}
    figgm, axgm = plt.subplots()
    for sig in Rec.Signals():
        
        GM[sig.name] = []
        fig, (AxPsd, Axt) = plt.subplots(2,1)        
        for isw, (t, vg) in enumerate(zip(SwTimes, Vgs)):
            if isw == len(SwTimes)-1:
                ts = sl.Signal.t_stop
            else:
                ts = SwTimes[isw+1]
            TWind = (t+StabTime, ts-GuardTime)
            s = sig.GetSignal(TWind, Units ='V')
            
            
            if s.name in DCch:
                s = (s*ivgain2-(-vg+Vds)*pq.V)/ivgain1
            else:
                s = s/(ivgain1*ACgain/ivgain2)
    
        
            Axt.plot(s.times, s, label=vg, alpha=0.5)     

            PSD = Ran.PlotPSD((s,),
                              Time = TWind,
                              Ax=AxPsd,
                              FMin=1,
                              Label=str(vg),
                              scaling='spectrum')
            
            psd = PSD[sig.name]['psd']
            Fpsd = PSD[sig.name]['ff']
           
            indicesPeak = np.where( ((Fpsd >= Fsig-4) & (Fpsd<=Fsig+4)))   

            IDSpeak = np.sqrt(psd[np.argmax(psd[indicesPeak])+indicesPeak[0][0]]+
                             psd[np.argmax(psd[indicesPeak])+indicesPeak[0][0]+1]+
                             psd[np.argmax(psd[indicesPeak])+indicesPeak[0][0]-1])
            
            gm = IDSpeak*1000/0.707

            GM[sig.name]=np.append(GM[sig.name],gm)
         
        if np.any(GM[sig.name]>= 2e-5) :
            ChKo[sig.name] = False
            color = 'k'
        else:
            ChKo[sig.name] = True      
            color = 'r'
        if sig.name in DCch:
            color = 'g'
        axgm.plot(Vgs, GM[sig.name], label=sig.name, color = color)   

        plt.close(fig)
        
    axgm.legend()
CalFile = {'GM': GM,
           'Ud0PBS':Ud0meanPBS,
           'Ud0InVivo':Ud0meanInVivo,
           'DeltaUd0':Ud0meanPBS-Ud0meanInVivo,
           'Vgs': Vgs-(Ud0meanPBS-Ud0meanInVivo)}
CalFile['ChKo'] = ChKo



dd.io.save('../Cal Data/'+'Cal-Gm'+File[:15]+'.h5', CalFile)