Plots and analysis for chapter 3: QCA characterization

In [1]:
# imports
import math
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
import qca
import coma
import configure_matplotlib
print(qca.__version__)
print(coma.__version__)
8.2.0-2.g60b16a3.dirty
unknown
In [2]:
# increases the resolution of the pngs displayed inline, 
# so they are big enough even for a small figsize
%matplotlib inline
c = %config InlineBackend.rc
c['savefig.dpi'] = 150
%config InlineBackend.figure_format='png'
%config InlineBackend.rc = c
%matplotlib inline

configure_matplotlib.configure_matplotlib()
In [3]:
savefig = False

Preliminary: Polarization at 45 degrees at larger distances

I had observed, for example in notebook 4 or the paper draft, that for diagonal wires ($\theta = 45^{\circ}$), a qualitative change in the polarization seems to occur as I go from shorter to longer cell-cell distances. At short distances, cell polarizations alternate, as expected. But at long distances the polarization of the second cell in a three-cell wire, as an example, flips and then all active cells' polarizations have the same sign. This is in disagreement with the derived $J$ at $45^{\circ}$ which should be always negative, yielding always alternating cell polarizations. Here we have a look at this again.

In [4]:
def Ps_diagonal():
    V1s,doas,theta,PD = [40,200],[3.0,6.0],45,1.0
    
    rs = []
    for V1 in V1s:
        for doa in doas:
            s = qca.QcaBond()
            s.T = 1
            s.q = 0.5
            s.l = qca.AngleWire(3,V1,doa,theta,PD)
            s.init()
            s.run()
            Ps = s.results['P']
            rs.append([V1,doa] + Ps)
    print(''.join([s.ljust(p) for s,p in 
                   zip(['V_1','doa','P_1','P_2','P_3'],[5,5,10,10,10])]))
    ss = ['{:3d}{:5.1f}{:10.5f}{:+10.5f}{:10.5f}'.format(*r) for r  in rs]
    for s in ss:
        print(s)
In [5]:
Ps_diagonal()
V_1  doa  P_1       P_2       P_3       
 40  3.0  -0.55058  +0.28950  -0.15535
 40  6.0  -0.01787  -0.00024  -0.00006
200  3.0  -0.99601  +0.99231  -0.98913
200  6.0  -0.09325  +0.00588  -0.00066

The table above prints the cell polarizations of a diagonal three-cell wire with an input polarization $P_D = 1$ at two different cell-cell distances and for two different $V_1$. At $V_1 = 40$ and $d/a = 3.0$ the cell polarizations alternate, as expected. But going up to $d/a = 6.0$ the second cell now has a negative polarization as well, which is confusing. However, for a larger $V_1$ the expected behaviour is recovered, even for $d/a = 6$.

Here is my (new) explanation of what is going on: We know that generally the response of a cell to the driver cell is stronger (gain, non-linear) than the response to another active cell (no gain, linear). Therefore, when the polarization of the adjacent (active) cell is very small, as for the second cell at $V_1 = 40$ and $d/a = 6.0$, then the influence of the driver cell, even though it is farther away, may outweigh the response to the adjacent cell. In this example, the driver cell polarization is $P_D = 1$, $J$ is negative and thus $P_2 < 0$. For larger $V_1$ the adjacent cell's polarization is larger and thus dominates over the influence of the driver cell.

Characterizing a three-cell finite-temperature wire

In [6]:
def text(p,t,loc='upper right',dx=0.035,dy=0.035):
    if loc=='upper right': x,y,v,h = 1-dx,1-dy,'top','right'
    elif loc=='upper left': x,y,v,h = dx,1-dy,'top','left'
    elif loc=='lower right': x,y,v,h = 1-dx,dy,'bottom','right'
    elif loc=='lower left': x,y,v,h = dx,dy,'bottom','left'
    p.text(x,y,t,verticalalignment=v,horizontalalignment=h,transform=p.transAxes)
In [7]:
# This is taken from my QCA paper notebook

from scipy.optimize import curve_fit

class CharacterizeAWire(object):
    def __init__(self):
        self.V1 = 40
        self.doa = 2.2
        self.T = 1
        self.N_plot = 100
        self._PP = []
        self._POverD = []
        self._POverDLog = []
        self._POverTheta = []
        self.legend_py = 0.3
        self.legend_px = 0.3
    
    def PP(self):
        if len(self._PP) == 0:
            self._PP = [(0.0,self.calculatePP(0.0)),
                        (0.5,self.calculatePP(0.5))]
        return self._PP
    
    def POverD(self):
        if len(self._POverD) == 0:
            self._POverD = [(0.0,self.calculatePOverD(0.0)),
                            (0.5,self.calculatePOverD(0.5))]
        return self._POverD

    def POverDLog(self):
        if len(self._POverDLog) == 0:
            self._POverDLog = [(0.0,self.calculatePOverDLog(0.0)),
                               (0.5,self.calculatePOverDLog(0.5))]
        return self._POverDLog

    def POverTheta(self):
        if len(self._POverTheta) == 0:
            self._POverTheta = [(0.0,self.calculatePOverTheta(0.0)),
                                (0.5,self.calculatePOverTheta(0.5))]
        return self._POverTheta
            
    def calculatePP(self, q): 
        r = []
        s = qca.QcaBond()
        for P in np.linspace(-1,1,self.N_plot):
            c = [P]
            s.q = q
            s.T = self.T
            s.t = 1
            s.l = qca.Wire(3,self.V1,self.doa-1,P)
            s.init()
            s.run()
            c.extend(s.results['P'])
            r.append(c)
        return np.array(r)

    def calculatePOverD(self, q):
        r = []
        s = qca.QcaBond()
        for boa in np.linspace(0.1,5,self.N_plot):
            c = [boa+1]
            s.q = q
            s.T = self.T
            s.t = 1
            s.l = qca.Wire(3,self.V1,boa,1.0)
            s.init()
            s.run()
            c.extend(s.results['P'])
            r.append(c)
        return np.array(r)
    
    def calculatePOverDLog(self, q):
        r = []
        s = qca.QcaBond()
        for boa in np.logspace(-1,2,self.N_plot):
            c = [boa+1]
            s.q = q
            s.T = self.T
            s.t = 1
            s.l = qca.Wire(3,self.V1,boa,1.0)
            s.init()
            s.run()
            c.extend(s.results['P'])
            r.append(c)
        return np.array(r)
    
    def calculatePOverTheta(self, q):
        r = []
        s = qca.QcaBond()
        for theta in np.linspace(0,180,self.N_plot):
            c = [theta]
            s.q = q
            s.T = self.T
            s.t = 1
            s.l = qca.AngleWire(3,self.V1,self.doa,theta,1.0)
            s.init()
            s.run()
            c.extend(s.results['P'])
            r.append(c)
        return np.array(r)
    
    def plotPP(self):
        a = self.PP()
        cols = 2
        rows = 1
        fig,ps = plt.subplots(rows,cols,figsize=(2.5*cols,2.5*rows))
        for p,(q,d) in zip(ps,a):
            p.plot(d[:,0],d[:,1],label='$P_1(P_D)$')
            p.plot(d[:,1],d[:,2],label='$P_2(P_1)$')
            p.plot(d[:,2],d[:,3],label='$P_3(P_2)$')
            p.set_xlabel('polarization')
            p.set_ylabel('polarization')
            p.set_xlim((-1,1))
            p.set_ylim((-1,1))
            text(p,'$q = {}$'.format(q),loc='lower right')
        ps[1].legend(loc='upper left',
                     labelspacing=self.legend_py,
                     handletextpad=self.legend_px)
        fig.tight_layout()
        plt.close()
        return fig
    
    def plotPOverD_old(self):
        a = self.POverD()
        cols = 2
        rows = 1
        fig,ps = plt.subplots(rows,cols,figsize=(2.5*cols,2.5*rows))
        for p,(q,d) in zip(ps,a):
            for i in range(3):
                p.plot(d[:,0],d[:,i+1],label='$P_{}$'.format(i+1))
            p.set_xlabel('$d/a$')
            p.set_ylabel('polarization')
            text(p,'$q = {}$'.format(q),loc='lower right')
        ps[1].legend(loc='upper right',
             labelspacing=self.legend_py,
             handletextpad=self.legend_px)
        fig.tight_layout()
        plt.close()
        return fig

    def plotPOverD(self):
        a_ = [self.POverD(),self.POverDLog()]
        cols = 2
        rows = 2
        fig,ps = plt.subplots(rows,cols,figsize=(2.5*cols,2.5*rows))
        for p,(q,d) in zip(ps[0],a_[0]):
            for i in range(3):
                p.plot(d[:,0],d[:,i+1],label='$P_{}$'.format(i+1))
            p.set_xlabel('$d/a$')
            p.set_ylabel('polarization')
            text(p,'$q = {}$'.format(q),loc='lower right',dy=0.07)
        
        for p,(q,d) in zip(ps[1],a_[1]):
            for i in range(3):
                p.plot(d[:,0],d[:,i+1],label='$P_{}$'.format(i+1))
            p.set_xlabel('$d/a$')
            p.set_ylabel('polarization')
            p.set_xscale('log')
            p.set_yscale('log')
            p.set_yticks(np.logspace(-10,0,6))
            p.set_xlim((1,100))
            p.set_ylim((1E-10,1))
            text(p,'$q = {}$'.format(q),loc='lower right',dx=0.16)

        text(ps[0,0],'(a)',loc='upper left',dy=0.05)
        text(ps[0,1],'(b)',loc='lower left')
        text(ps[1,0],'(c)',loc='lower left')
        text(ps[1,1],'(d)',loc='lower left')
        ps[0,1].legend(loc='upper right',
             labelspacing=self.legend_py,
             handletextpad=self.legend_px)
        fig.tight_layout()
        plt.close()
        return fig

    def plotPOverDLog(self):
        a = self.POverDLog()
        cols = 2
        rows = 1
        fig,ps = plt.subplots(rows,cols,figsize=(2.5*cols,2.5*rows))
        for p,(q,d) in zip(ps,a):
            for i in range(3):
                p.plot(d[:,0],d[:,i+1],label='$P_{}$'.format(i+1))
            p.set_xlabel('$d/a$')
            p.set_ylabel('polarization')
            p.set_xscale('log')
            p.set_yscale('log')
            p.set_xlim((1,100))
            p.set_ylim((1E-10,1))
            text(p,'$q = {}$'.format(q),loc='lower right')
        ps[1].legend(loc='upper right',
             labelspacing=self.legend_py,
             handletextpad=self.legend_px)
        fig.tight_layout()
        plt.close()
        return fig

    def fitPOverDLog(self):
        f = lambda x,a,b: a*x**b
        r = []
        a = self.POverDLog()
        for q,d in a:
            #print('q = {}'.format(q))
            #e = d[ (d[:,0] > 2) & (d[:,0] < 3) ] # also very interesting
            row = []
            e = d[d[:,0] > 10]
            for i in range(3):
                #print('P_{}'.format(i+1))
                ps,cov = curve_fit(f, e[:,0], e[:,i+1])
                es = np.sqrt(cov.diagonal())
                row.append(zip(ps,es))
                #print('  a = {:.4f} +- {:.4f}\n  b = {:.4f} +- {:.4f}'.format(ps[0], es[0], ps[1], es[1]))
            r.append((q,row))
        return r
    
    def plotPOverTheta(self):
        a = self.POverTheta()
        rows = 1
        cols = 2
        fig,ps = plt.subplots(rows,cols,figsize=(2.5*cols,2.5*rows))
        for p,(q,d) in zip(ps,a):
            for i in range(3):
                p.plot(d[:,0], d[:,i+1],label='$P_{}$'.format(i+1))
            p.set_xlabel('$\\theta$')
            p.set_ylabel('polarization')
            p.set_xticks([0,45,90,135,180])
            text(p,'$q = {}$'.format(q),loc='lower right')
        ps[1].legend(loc='upper right',
             labelspacing=self.legend_py,
             handletextpad=self.legend_px)
        text(ps[0],'(a)',dx=0.05,loc='lower left')
        text(ps[1],'(b)',loc='lower left')
        fig.tight_layout()
        plt.close()
        return fig
In [8]:
e = coma.Experiment('../experiments/experiment.000045/')
rs = e.retrieve_results(
    (
         ('P_D','parameters/layout/P'),
         ('P','results/P'),
    ),
    (
         ('model','parameters/model'),
         ('q','parameters/q'),
         ('T','parameters/T'),
         ('boa','parameters/layout/boa'),
    )
)
#rs
In [9]:
def PP_figure(rs):
    cols = 2
    rows = 2
    ls = ['(a)','(b)','(c)','(d)']
    ls.reverse()
    fig,ps = plt.subplots(rows,cols,figsize=(2.5*cols,2.5*rows))
    for ps_,T in zip(ps,[1,1E-9]):
        for p,r in zip(ps_,rs['T',T]):
            d = np.row_stack((-r.table,r.table))
            d = d[ d[:,0].argsort() ]
            p.plot(d[:,0],d[:,1],label='$P_1(P_D)$')
            p.plot(d[:,1],d[:,2],label='$P_2(P_1)$')
            p.plot(d[:,2],d[:,3],label='$P_3(P_2)$')
            p.set_xlabel('polarization')
            p.set_ylabel('polarization')
            p.set_xlim((-1,1))
            p.set_ylim((-1,1))
            q = r.parameters.q
            T = r.parameters.T
            if T < 0.01: T = 0.0
            l = ls.pop()
            text(p,l,loc='upper left')
            dy = 0.38 if l=='(b)' else 0.035
            text(p,'$T {} {:d}$\n$q = {}$'
                 .format('=' if T!=0 else '\sim',int(T),q),
                 dy=dy,
                 loc='lower right')
    ps[0,1].legend(loc='lower right',
                   labelspacing=0.3,
                   handletextpad=0.3)
    fig.tight_layout()
    plt.close()
    return fig
In [10]:
w = CharacterizeAWire()
w.N_plot = 100
# w.doa = 3.0 # PP plot is not very meaningful for doa = 3.0
# w.T = 2 # again: PP plot not good -- only very small gain; not fully polarized
f1 = PP_figure(rs)
f2 = w.plotPOverD()
fit = w.fitPOverDLog()
f3 = w.plotPOverTheta()

if savefig:
    f1.savefig('../plots/chapter03/three_cells_PP.pdf')
    f2.savefig('../plots/chapter03/three_cells_P_over_d.pdf')
    f3.savefig('../plots/chapter03/three_cells_P_over_theta.pdf')

display(f1)
display(f2)
display(f3)
In [11]:
# Find the b/a values corresponding to maximum / large values of P_3
def calculate_P3_over_boa(boas, q):
    P3s = []
    for boa in boas:
        s = qca.QcaBond()
        s.T = 1
        s.t = 1
        s.q = q
        s.l = qca.Wire(3,40,boa,1)
        s.init()
        s.run()
        P3s.append(s.results['P'][2])
    m = np.column_stack((boas,np.array(P3s)))
    return m

m1 = calculate_P3_over_boa(np.linspace(1.8,2.2,100),0)
m1_max = m1[ m1[:,1].argmax() ]
print('q = 0, boa_max = {}, P_3max = {}'.format(m1_max[0], m1_max[1]))

m2 = calculate_P3_over_boa(np.linspace(0.2,1.5,100),0.5)
m2_ = m2[ m2[:,1]>0.9 ]
print('q = 1/2, P_3 > 0.9 in the range boa = {} ... {}'.format(m2_[0,0],m2_[-1,0]))
q = 0, boa_max = 1.92121212121, P_3max = 0.0551161973124
q = 1/2, P_3 > 0.9 in the range boa = 0.331313131313 ... 1.26363636364

Workable parameters

In [12]:
e1 = coma.Experiment('../experiments/experiment.000046/')
rs1 = e1.retrieve_results(
    (
        ('V1', 'parameters/V1'),
        ('P_D', 'results/P_D'),
        ('P_D_2', 'parameters/model/parameters/layout/P'),
        ('P', 'parameters/model/results/P')
    ),
    (
        ('model', 'parameters/type'),
        ('boa', 'parameters/boa')
    ))
e2 = coma.Experiment('../experiments/experiment.000047/')
rs2 = e2.retrieve_results(
    (
         ('N', 'parameters/N'),
         ('V1_crit','results/V1_crit'),
    ),
    (
         ('model','parameters/model_type'),
         ('boa','parameters/boa'),
    )
)
#rs1,rs2
In [13]:
# Error of Ising versus bond at boa = 1.2
rs2['boa',1.2][1].table[:,1] - rs2['boa',1.2][0].table[0:5,1]
Out[13]:
array([ 2.77491474,  1.85757179,  0.797192  , -0.24793053, -1.25681705])
In [14]:
# Error of Ising versus bond at boa = 2.0
rs2['boa',2.0][1].table[:,1] - rs2['boa',2.0][0].table[0:5,1]
Out[14]:
array([ 0.89064274,  0.41003895, -0.02860737, -0.46153221, -0.89445705])
In [15]:
def plotSelfConsistentPOverV1(p, r):
    d = r.table
    d = d[ d[:,0].argsort() ] # Sort by first column
    p.plot(d[:,0],d[:,1],label='$P_D$')
    for i in range(3):
        p.plot(d[:,0],d[:,3+i],label='$P_{}$'.format(i+1))
    p.legend(loc='lower right',
             labelspacing=0.3,
             handletextpad=0.3)
    p.set_xlabel('$V_1$')
    p.set_ylabel('polarization')
    #p.set_xlim((0,60))
    p.set_xlim((0,250))

def plotV1critOverN(p,rs):
    d1 = rs['model','QcaIsing'][0].table
    d2 = rs['model','QcaBond'][0].table
    p.plot(d1[:,0],d1[:,1],'+-',label='Ising')
    p.plot(d2[:,0],d2[:,1],'+', color='red', markersize=8,label='Bond')
    p.set_xlabel('number of active cells')
    p.set_ylabel('$V_{1crit}$')
    p.set_xlim((0,13))
    #p.set_ylim((0,50))
    #p.set_ylim((50,200))
    p.legend(loc='lower right')
    
def createMeanFieldFigure():
    rows = 1
    cols = 2
    fig,ps = plt.subplots(rows,cols,figsize=(cols*2.5,rows*2.5))
    
    boa = 2.0
    #boa = 1.2
    
    p = ps[0]
    r = rs1['model','QcaBond']['boa',boa][0]
    plotSelfConsistentPOverV1(p, r)
    text(p,'(a)',loc='upper left')
    
    p = ps[1]
    plotV1critOverN(p,rs2['boa',boa])
    text(p,'(b)',loc='upper left')

    fig.tight_layout()
    plt.close()
    return fig
In [16]:
f = createMeanFieldFigure()
if savefig:
    f.savefig('../plots/chapter03/critical_V1.pdf')
display(f)
In [17]:
r = rs1['boa',2.0]['model','QcaBond'][0]
d = r.table
P_sc = d[ d[:,0] == 200 ][0][1]
print('Self-consistent polarization at V1 = 200: P_D = {:.2f}'.format(P_sc))
Self-consistent polarization at V1 = 200: P_D = 0.97
In [18]:
d1 = rs2['model','QcaBond']['boa',1.2][0].table
d2 = rs2['model','QcaIsing']['boa',1.2][0].table
print('Some critical V1 for boa = 1.2')
print(' 1 active cell:  V_1crit_bond = {:.2f}, V_1crit_ising = {:.2f}'.format(d1[0,1],d2[0,1]))
print(' 5 active cells: V_1crit_bond = {:.2f}, V_1crit_ising = {:.2f}'.format(d1[4,1],d2[4,1]))
print('12 active cells: V_1crit_bond = n.a.   V_1crit_ising = {:.2f}'.format(d2[11,1]))
Some critical V1 for boa = 1.2
 1 active cell:  V_1crit_bond = 19.02, V_1crit_ising = 16.24
 5 active cells: V_1crit_bond = 30.07, V_1crit_ising = 31.33
12 active cells: V_1crit_bond = n.a.   V_1crit_ising = 45.20
In [19]:
# A couple of V1_crit values for different T, boa, for referencing in the thesis

e = coma.Experiment('../experiments/experiment.000050/')
rs = e.retrieve_results(
    (
         ('T', 'parameters/T'),
         ('boa','parameters/boa'),
         ('V1_crit','results/V1_crit'),
    ))
r = rs[0]
print(r.table_columns)
print(r.table)
('T', 'boa', 'V1_crit')
[[   1.            1.2          14.71359739]
 [   1.            1.5          24.1082572 ]
 [   1.            2.           55.09766588]
 [   1.            2.5         117.72682409]
 [   1.            3.          230.23960428]
 [   1.            3.5         416.87026186]
 [   1.            4.          708.92288656]
 [   2.            1.2          25.99443636]
 [   2.            1.5          44.41186018]
 [   2.            2.          107.89733229]
 [   2.            2.5         234.29422197]
 [   2.            3.          459.87285814]
 [   2.            3.5         833.4011754 ]]

Cell polarization along the wire

In [20]:
e = coma.Experiment('../experiments/experiment.000049/')
rs = e.retrieve_results(
    (
         ('P','results/P'),
    ),
    (
         ('model','parameters/model'),
         ('T','parameters/T'),
         ('N','parameters/layout/N'),
         ('V1','parameters/layout/V1'),
         ('boa','parameters/layout/boa'),
         ('PD','parameters/layout/P'),
    )
)
#rs
In [21]:
# adapted from notebook 7


# fitting
from scipy.optimize import curve_fit

f = lambda x,a,b: a*b**x

def fit_to_exponential(d):
    xs = np.arange(1,len(d)+1,1)
    #xs = np.arange(0,len(d),1) # alternate definition of \chi_e  and \chi
    ys = d
    ps,cov = curve_fit(f, xs, ys)
    es = np.sqrt(cov.diagonal())
    return np.concatenate((ps,es))

def fit_rs(rs,printit=False):
    rows = []
    for r in rs:
        if r.parameters.N < 3:
            continue
        ps = fit_to_exponential(r.table[0])
        rows.append(np.concatenate((np.array([r.parameters.N]), ps)))
    d = np.vstack(rows)
    
    n = len(d[:,1])
    
    a_min = min(d[:,1])
    a_max = max(d[:,1])
    a_avg = sum(d[:,1])/n
    e_a_avg = sum(d[:,3])/n
    
    b_min = min(d[:,2])
    b_max = max(d[:,2])
    b_avg = sum(d[:,2])/n
    e_b_avg = sum(d[:,4])/n
    
    if printit:
        print('Fitting to f(x) = a b^x')
        print('-----------------------')
        print('Parameters: ' + str(r.parameters))
        print('Average a: {:.8g}'.format(a_avg))
        print('Average error for a: {:.8g}'.format(e_a_avg))
        print('Average b: {:.8g}'.format(b_avg))
        print('Average error for b: {:.8g}'.format(e_b_avg))
        print('\n')
    
    return d,((a_min,a_max,a_avg,e_a_avg),(b_min,b_max,b_avg,e_b_avg))

def compute_chis(rs):
    chis = []
    for r in rs:
        if r.parameters.N < 2:
            continue
        d = np.concatenate(([1.0],r.table[0]))
        chis_ = [d[i+1]/d[i] for i in range(len(d)-1)]
        chis.append([r.parameters.N] + chis_)
    return chis

# plots
def plot_P_for_different_wires(p,rs):
    params = [
        # V1,boa
        (200,1.2),
        (200,2.0,)
    ]
    cs = ['b','g','r','c','m','y'] # 'k'
    Ns = [12,10,8,6,4,2]
    for V1,boa in params:
        for N,c in zip(Ns,cs):
            r = rs['V1',V1]['boa',boa]['N',N][0]
            d = r.table
            params = r.parameters
            p.plot(range(1,params.N+1),d[0],'-',
                   label='{}-cell wire'.format(params.N),color=c)
    p.set_xlim((0,13))
    p.set_ylim((0,1.05))
    p.set_xlabel('cell number')
    p.set_ylabel('polarization')
    text(p,'(a)',loc='lower left')
    p.text(8.5,0.9,'$d/a = 2.2$')
    p.text(8.5,0.4,'$d/a = 3.0$')
    #p.legend(loc='upper left', bbox_to_anchor=(1,1))

def plot_fit_parameters_over_N(p,rs):
    params = [
        # V1,boa
        (200,1.2),
        (200,2.0)
    ]
    for V1,boa in params:
        d,v = fit_rs(rs['V1',V1]['boa',boa])
        p.plot(d[:,0],d[:,1],'-',label='$\\chi_e(d/a={})$'.format(boa+1))
        p.plot(d[:,0],d[:,2],'-',label='$\\chi(d/a={})$'.format(boa+1))
    p.set_xlim((2,13))
    p.set_ylim((0.7,1.01))
    p.set_xlabel('number of cells in wire')
    p.set_ylabel('fit parameters')
    text(p,'(b)',loc='lower left')
    p.legend(loc='lower right',                    
             labelspacing=0.3,
             handletextpad=0.3)

def polarization_of_wires_figure(P_D=1.0):
    n_r,n_c = 1,2
    fig,ps = plt.subplots(n_r,n_c,figsize=(2.5*n_c,2.5*n_r))
    plot_P_for_different_wires(ps[0],rs['V1',200]['PD',P_D])
    plot_fit_parameters_over_N(ps[1],rs['V1',200]['PD',P_D])
    fig.tight_layout()
    plt.close()
    return fig

def polarization_response_figure(P_D=1.0):
    fig,ps = plt.subplots(1,2,figsize=(5,2.5))
    boas = [2.0,3.0]
    for p,boa in zip(ps,boas):
        chis = compute_chis(rs['V1',200]['boa',boa]['PD',P_D])
        for chis_ in chis:
            p.plot(np.arange(1,len(chis_[1:])+1,1),chis_[1:])
        p.set_xlabel('cell number')
        p.set_ylabel('$P_i / P_{i-1}$')
        p.set_xlim((0,13))
        #p.set_ylim((0,1))
        text(p,'$d/a = {}$'.format(boa+1),loc='lower right')
    fig.tight_layout()
    plt.close()
    return fig
In [22]:
f_ = polarization_of_wires_figure()
if savefig:
    f_.savefig('../plots/chapter03/wire_polarization_old__.pdf')
display(f_)
In [23]:
def number_of_cells(P,a,b):
    return math.log(P/a)/math.log(b)

def print_wire_polarization_estimates(P_D=1.0):
    boas = [1.2,2.0]
    V1 = 200
    for boa in boas:
        d,((_,_,a,_),(_,_,b,_)) = fit_rs(rs['V1',V1]['boa',boa]['PD',P_D],printit=True)
        print('b/a = {}'.format(boa))
        print('Number of cells before polarization drops below 0.9: {:.2f}'
              .format(number_of_cells(0.9,a,b)))
        print('Number of cells before polarization drops below 0.1: {:.2f}'
              .format(number_of_cells(0.1,a,b)))
        print('\n\n')
In [24]:
print_wire_polarization_estimates(P_D=1.0)
Fitting to f(x) = a b^x
-----------------------
Parameters: (model=QcaIsing,T=2.0,N=12,V1=200,boa=1.2,PD=1.0)
Average a: 0.99911477
Average error for a: 5.11369e-07
Average b: 0.99999702
Average error for b: 1.4143317e-07


b/a = 1.2
Number of cells before polarization drops below 0.9: 35040.81
Number of cells before polarization drops below 0.1: 771988.54



Fitting to f(x) = a b^x
-----------------------
Parameters: (model=QcaIsing,T=2.0,N=12,V1=200,boa=2.0,PD=1.0)
Average a: 0.97367401
Average error for a: 0.0022100875
Average b: 0.88404441
Average error for b: 0.00069302451


b/a = 2.0
Number of cells before polarization drops below 0.9: 0.64
Number of cells before polarization drops below 0.1: 18.47



In [25]:
polarization_response_figure()
Out[25]:

This is a new graph and goes slightly beyond the analysis in notebook 7. I don't think I want to inlcude it in the thesis, but it is still a nice plot: it shows how the polarization response really looks like inside the wire. The response is constant over the majority of cells, but falls off over the last two cells. It's also clearly visible that the response of the first cell---to the driver cell---is very different. At larger distances $d/a = 4$, the behaviour is qualitatively different (note that here the actual polarizations are almost zero).

It is clear that the right most cells' response falls off. An argment why the response grows for the first couple of cells: Maybe the cells to the left still see the driver cell and hence have a slightly higher polarization --- this would make the response with respect to them smaller.

Maybe the most significant finding of this graph is that the polarization response falls off towards both edges (left and right; but less to the left). I should definitely mention that in the thesis.

In [26]:
# Just to check: two graphs for different input polarizations
#polarization_of_wires_figure(P_D=0.5)
polarization_of_wires_figure(P_D=0.1)
Out[26]:
In [27]:
#polarization_response_figure(P_D=0.5)
polarization_response_figure(P_D=0.1)
Out[27]:

The bottom line of the above two graphs is that the our simple physical model actually works very well. When changing $P_D$, it is only $\chi_e$ that changes. $\chi$ remains virtually unchanged. For small $P_D$, $\chi_e$ seems to drop off faster with the number of cells per wire, for $d/a = 2.2$. Due to the gain in the response to the driver cell, even an input polarization of only $P_D = 0.5$ still yields very good results, especially at $d/a = 2.2$.

In [28]:
e1 = coma.Experiment('../experiments/experiment.000051/')
rs1 = e1.retrieve_results(
    (
         ('boa','parameters/layout/boa'),
         ('P','results/P'),
    ),
    (
         ('model','parameters/model'),
         ('T','parameters/T'),
         ('N','parameters/layout/N'),
         ('V1','parameters/layout/V1'),
         ('PD','parameters/layout/P'),
    )
)
e2 = coma.Experiment('../experiments/experiment.000052/')
rs2 = e2.retrieve_results(
    (
         ('N','parameters/model/parameters/layout/N'),
         ('boa','results/boa'),
    ),
    (
         ('T','parameters/model/parameters/T'),
         ('V1','parameters/model/parameters/layout/V1'),
         ('PD','parameters/model/parameters/layout/P')
    )
)

#rs1,rs2
In [29]:
def plot_output_polarization_over_cell_cell_distance(p,rs):
    Ns = [2,4,6,8,10,12]
    for N in Ns:
        r = rs['N',N][0]
        d = r.table
        p.plot(d[:,0]+1,d[:,-1],label='{} cells'.format(r.parameters.N))
    p.set_ylim((0,1.02))
    p.set_ylabel('output polarization')
    p.set_xlabel('$d/a$')
    p.legend(loc='upper right',
             labelspacing=0.3,
             handletextpad=0.3)
    
def plot_threshold_boa(p,rs):
    for r in rs:
        d = r.table
        p.plot(d[:,0],d[:,1]+1,label='$V_1 = {}$'.format(r.parameters.V1))
    p.set_ylim((2.0,3.4))
    p.set_xlim((0,13))
    p.set_xlabel('number of cells in wire')
    p.set_ylabel('$d_{th}/a$')
    p.legend(loc='upper right',
             labelspacing=0.3,
             handletextpad=0.3)

def polarization_cell_cell_distance_figure(P_D=1.0):
    fig,ps = plt.subplots(1,2,figsize=(5,2.5))
    
    p = ps[0]
    plot_output_polarization_over_cell_cell_distance(p,rs1['V1',200]['T',2.0]['PD',P_D])
    text(p,'(a)',loc='lower left')
    
    p = ps[1]
    plot_threshold_boa(p,rs2['PD',P_D])
    text(p,'(b)',loc='lower left')
    
    fig.tight_layout()
    plt.close()
    return fig
In [30]:
f = polarization_cell_cell_distance_figure()
if savefig:
    f.savefig('../plots/chapter03/cell_cell_distance_threshold_old__.pdf')
display(f)
In [31]:
# Check how it changes for small polarizations
polarization_cell_cell_distance_figure(P_D=0.1)
/Users/burkhard/anaconda/lib/python2.7/site-packages/matplotlib/axes.py:4747: UserWarning: No labeled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labeled objects found. "
Out[31]:

Even though here the region where we have perfect output polarization is much smaller (and at unphysical small cell-cell distances), it remains true that over a broad range $d/a \sim 1 \ldots 2.5$ the output polarization is the same independent of the number of cells in the wire. And this is the same range as for $P_D = 1.0$ above. This is consistent with the graph at $P_D = 0.1$ above (cell polarization along wire, and also polarization response along wire), where we had seen that $\chi$ is independent of $P_D$. Because $\chi \sim 1$ for the chosen parameters at small cell-cell distances, the output polarization for all wires here is virtually the same. In the region $d/a > 2.5$ we have $\chi < 1$ and that's why the polarization drops rapidly to zero, and more rapidly for longer wires.

The Silicon system

In [22]:
# adapted from notebook 5a

import scipy.optimize
from scipy.constants import Boltzmann,elementary_charge

def fit(f, xs, ys):
    ps,cov = scipy.optimize.curve_fit(f,xs,ys)
    es = np.sqrt(cov.diagonal())
    labels = ['a','b','c','d','e','f','g']
    print('Fit results: ')
    for l,p,e in zip(labels, ps, es):
        print('   {} = {:.4g} +- {:4g}'.format(l,p,e))
    return (lambda x: f(x,*ps),) + tuple(ps)

def plot_silicon():
    # data from Pitters et al. The Journal of Chemical Physics, 
    # 134, 064712 (2011), DOI:10.1063/1.3514896
    # Columns: distance in Angstrom, t_ij in eV, V_ij in eV
    table1 = [[3.84,0.1538,0.3208],
              [7.68,0.0439,0.1754],
              [9.93,0.0303,0.1275],
              [10.86,0.0254,0.1133],
              [11.52,0.0223,0.1046],
              [11.74,0.0213,0.1019],
              [12.55,0.0179,0.0927]]
    
    # fit
    d = np.array(table1)
    xs = np.linspace(1,30,50)
    d_ = d[ d[:,0]>8 ]
    f_t,a_t,b_t = fit(lambda x,a,b: a*np.exp(b*x), d_[:,0],d_[:,1])
    d_ = d
    f_V,b_V = fit(lambda x,b: b/x, d_[:,0],d_[:,2])
    
    # plot
    fig,ps = plt.subplots(1,3,figsize=(7.5,2.5))

    p = ps[0]
    p.plot(d[:,0],d[:,1],'o-',label='data')
    p.plot(xs,f_t(xs),
           label='$f_t \sim {:.2f} e^{{ {:.2f} r}}$'.format(a_t,b_t))
    p.set_yscale('log')
    p.set_xlabel('distance $(\AA)$')
    p.set_ylabel('$t_{ij} (\mathrm{eV})$')
    p.legend(loc='upper right')
    
    p = ps[1]
    p.plot(d[:,0],d[:,2],'o-',label='data')
    p.plot(xs,f_V(xs),
           label='$f_V \sim {:.2f} r^{{-1}}$'.format(b_V))
    p.set_xlabel('distance $(\AA)$')
    p.set_ylabel('$V_{ij} (\mathrm{eV})$')
    p.legend(loc='upper right')
    
    p = ps[2]
    p.plot(xs,f_V(xs)/f_t(xs),
           label='$f_V/f_t$')
    p.set_xlabel('distance $(\AA)$')
    p.set_ylabel('$V_{ij}/t_{ij}$')
    p.legend(loc='upper right')
    
    fig.tight_layout()
    plt.close()
    return fig,f_t,f_V

def print_silicon_values(f_t,f_V,x=32,T=1):
    T_over_t = T
    T_in_Kelvin = T_over_t * f_t(x) * elementary_charge / Boltzmann
    print('x = {}\nT/t = {:.2f}\nV1/t = {:.2f} \nT = {:.2f} K\n'
          't = {:.2f} meV\nV1 = {:.2f} meV\n'
          .format(x,T_over_t,f_V(x)/f_t(x), T_in_Kelvin,f_t(x)*1000.0,f_V(x)*1000.0))
In [23]:
fig,f_t,f_V = plot_silicon()
display(fig)
#print_silicon_values(f_t,f_V,36,1)
print_silicon_values(f_t,f_V,32,1)
print_silicon_values(f_t,f_V,22,0.1)
Fit results: 
   a = 0.2177 +- 0.00810395
   b = -0.1982 +- 0.0033719
Fit results: 
   a = 1.243 +- 0.0188579
x = 32
T/t = 1.00
V1/t = 101.45 
T = 4.44 K
t = 0.38 meV
V1 = 38.85 meV

x = 22
T/t = 0.10
V1/t = 20.33 
T = 3.23 K
t = 2.78 meV
V1 = 56.50 meV

Majority gate

In [34]:
e = coma.Experiment('../experiments/experiment.000040/')
rs = e.retrieve_results(
    (
         ('I1','parameters/layout/I1'),
         ('I3','parameters/layout/I3'),
         ('P','results/P'),
    ),
    (
         ('T','parameters/T'),
         ('N_lead','parameters/layout/N_lead'),
         ('V1','parameters/layout/V1'),
         ('doa','parameters/layout/doa'),
         ('I2','parameters/layout/I2'),
    )
)
#rs
In [35]:
def create_majority_gate_plot(rs):
    n_r,n_c = 2,2
    fig,ps = plt.subplots(n_r,n_c,figsize=(5,4.1),sharex=True,sharey=True)
    
    V = np.arange(-1.0,1.01,0.1)
    doas = [3.0,2.2]
    Is = [-1,1]
    col = lambda d,c: d[:,c].reshape(30,30)
    ls = ['(a)','(b)','(c)','(d)']
    ls.reverse()
    cs = []
    
    for doa,ps_ in zip(doas,ps):
        for I,p in zip(Is,ps_):
            d = rs['I2',I]['doa',doa][0].table[:900]
            cs1 = p.contourf(col(d,0), col(d,1), col(d,10),V)
            cs2 = p.contour(col(d,0), col(d,1), col(d,10),[0.0],
                            colors='black',linewidths=2)
            cs.append((cs1,cs2))
            p.grid(True, color='#666666')
            text(p, '$d/a = {}$'.format(doa), loc='lower right')
            text(p, ls.pop(), loc='lower left')
    ps[0,0].set_title('$I_1$ AND $I_3$')
    ps[0,1].set_title('$I_1$ OR $I_3$')
    ps[0,0].set_ylabel('$I_3$')
    ps[1,0].set_ylabel('$I_3$')
    ps[1,0].set_xlabel('$I_1$')
    ps[1,1].set_xlabel('$I_1$')

    fig.tight_layout()
    
    for cs_,ps_ in zip(cs,ps):
        cbar = fig.colorbar(cs_[0],ax=ps_.tolist(),ticks=[-1,-0.5,0,0.5,1])
        cbar.add_lines(cs_[1])
    plt.close()
    return fig
In [36]:
fig = create_majority_gate_plot(rs)
# if savefig:
#     fig.savefig('../plots/chapter03/majority_gate_contour.pdf')
display(fig)

This was my original majority gate contour plot, as in notebook 8. It then occured to me that it is more meaningful to plot not against the driver cell polarizations, but the cells right next to it. The cells are numbered:

          3
          2
      5 4 1 8 9
          6
          7

Therefore, let's plot $P_9 \doteq O$ against $P_3 \doteq I_1^{\prime}$ and $P_7 \doteq I_3^{\prime}$; we leave I_2 fixed. This new graph is below.

In [37]:
# r = rs['doa',3.0]['I2',-1][0]
# r.table[:900,(0,1,4,6,8,10)][:10]
# r.table[:900,(4,8,10)][:10]
#
# cells are numbered:
#           3
#           2
#       5 4 1 8 9
#           6
#           7
In [38]:
e = coma.Experiment('../experiments/experiment.000054/')
rs = e.retrieve_results(
    (
         ('I1','parameters/layout/I1'),
         ('I3','parameters/layout/I3'),
         ('P','results/P'),
    ),
    (
         ('T','parameters/T'),
         ('N_lead','parameters/layout/N_lead'),
         ('V1','parameters/layout/V1'),
         ('doa','parameters/layout/doa'),
         ('I2','parameters/layout/I2'),
    )
)
#rs
In [39]:
def create_majority_gate_plot_prime(rs):
    n_r,n_c = 2,2
    fig,ps = plt.subplots(n_r,n_c,figsize=(5,4.1),sharex=True,sharey=True)
    
    V = np.arange(-1.0,1.01,0.1)
    doas = [2.2,2.6]
    Is = [-1,1]
    col = lambda d,c: d[:,c].reshape(30,30)
    ls = ['(a)','(b)','(c)','(d)']
    ls.reverse()
    cs = []
    
    for doa,ps_ in zip(doas,ps):
        for I,p in zip(Is,ps_):
            d = rs['I2',I]['doa',doa][0].table[:900]
            cs1 = p.contourf(col(d,4), col(d,8), col(d,10),V)
            cs2 = p.contour(col(d,4), col(d,8), col(d,10),[0.0],
                            colors='black',linewidths=2)
            cs.append((cs1,cs2))
            p.grid(True, color='#666666')
            text(p, '$d/a = {}$'.format(doa), loc='lower right')
            text(p, ls.pop(), loc='lower left')
    ps[0,0].set_title('$I_1$ AND $I_3$')
    ps[0,1].set_title('$I_1$ OR $I_3$')
    ps[0,0].set_ylabel('$I_3^{\prime}$')
    ps[1,0].set_ylabel('$I_3^{\prime}$')
    ps[1,0].set_xlabel('$I_1^{\prime}$')
    ps[1,1].set_xlabel('$I_1^{\prime}$')

    fig.tight_layout()
    
    for cs_,ps_ in zip(cs,ps):
        cbar = fig.colorbar(cs_[0],ax=ps_.tolist(),ticks=[-1,-0.5,0,0.5,1])
        cbar.add_lines(cs_[1])
    plt.close()
    return fig
In [40]:
fig = create_majority_gate_plot_prime(rs)
if savefig:
    fig.savefig('../plots/chapter03/majority_gate_contour.pdf')
display(fig)

Why is the ouput polarization asymmetric for larger d/a? I guess it is because the polarization from cell to cell decays so fast so that the diagonal interactions become more important. E.g. if $P_2 \sim 1$ and $P_6 \sim 1$, then they will tend to set $P_8 \sim -1$ via their diagonal interaction. If then polarizations decay fast, so that e.g. $|P_1| < |P_2|$ and additionally, as an exmaple, $I_2 = -1$, then this will shift the point at which the output polarization switches from positive to negative.

More on cell polarizations / workable parameters

In [41]:
e1 = coma.Experiment('../experiments/experiment.000057/')
rs1 = e1.retrieve_results(
    (
         ('V1','parameters/layout/V1'),
         ('P','results/P'),
    ),
    (
         ('model','parameters/model'),
         ('T','parameters/T'),
         ('N','parameters/layout/N'),
         ('boa','parameters/layout/boa'),
         ('PD','parameters/layout/P'),
    )
)
e2 = coma.Experiment('../experiments/experiment.000058/')
rs2 = e2.retrieve_results(
    (
         ('boa','parameters/layout/boa'),
         ('P','results/P'),
    ),
    (
         ('model','parameters/model'),
         ('T','parameters/T'),
         ('N','parameters/layout/N'),
         ('V1','parameters/layout/V1'),
         ('PD','parameters/layout/P'),
    )
)
#rs1,rs2
In [42]:
def compute_chis(d):
    return np.column_stack([d[:,0]] + [d[:,i+2] / d[:,i+1] for i in range(d.shape[1]-2)])

def plot_chis_over_V1(p,rs):
    models = ['QcaIsing','QcaBond','QcaFixedCharge']
    ls = {'QcaIsing': 'Ising','QcaBond':'Bond','QcaFixedCharge':'Fixed'}
    cs = [[4],[1],[0]]
    # Ts = [0.1,0.01]
    # T = 0.01 actually gives a slightly higher minimum V_1 than T = 0.1
    Ts = [2.0,1.0,0.1] 
    boa = 1.2
    for T in Ts:
        cols = ['b','g','r']
        for m,c,col in zip(models,cs,cols):
            chis = compute_chis(rs['T',T]['model',m]['boa',boa][0].table)
            chis = chis[ chis[:,0].argsort()]
            for c_ in c:
                p.plot(chis[:,0],chis[:,1+c_],color=col,label=ls[m])
    p.set_xlabel('$V_1$')
    p.set_ylabel('$\chi$')
    p.set_xlim((0,140))
    p.set_ylim((0.9,1.015))
    p.text(74,0.97,'$T = 2$')
    p.text(46,0.983,'$T = 1$')
    p.text(10,1.002,'$T = 0.1$')
    
    handles,labels = p.get_legend_handles_labels()
    p.legend(handles[:3],labels[:3],loc='lower right')

def plot_chis_over_doa(p,rs):
    models = ['QcaFixedCharge']
    cs = [[0]]
    ls = {'QcaIsing': 'Ising','QcaBond':'Bond','QcaFixedCharge':'Fixed'}
    # Ts = [0.1]
    # V1s = [20,40,60,100,200]
    # at T = 0.1 and V1 = 20: chi < 1 for all doa
    Ts = [1.0]
    V1s = [40,60,100,200]
    for T in Ts:
        for V1 in V1s:
            for m,c in zip(models,cs):
                chis = compute_chis(rs['V1',V1]['model',m]['T',T][0].table)
                chis = chis[ chis[:,0].argsort()]
                for c_ in c:
                    p.plot(chis[:,0]+1,chis[:,1+c_],label=' $V_1 = {}$'.format(V1))
    p.set_xlabel('$d/a$')
    p.set_ylabel('$\chi$')
    p.set_ylim((0.9,1.015))
    p.legend(loc='lower right',
             labelspacing=0.15,
             handletextpad=-0.2)

def create_chi_figure(rs1,rs2):
    fig,ps = plt.subplots(1,2,figsize=(5,2.5))
    
    p = ps[0]
    plot_chis_over_V1(p,rs1)
    text(p,'(a)',loc='upper right')
    
    p = ps[1]
    plot_chis_over_doa(p,rs2)
    text(p,'(b)',loc='upper right')
    
    fig.tight_layout()
    plt.close()
    return fig
In [43]:
fig = create_chi_figure(rs1,rs2)
if savefig:
    fig.savefig('../plots/chapter03/chis.pdf')
display(fig)