dunno.

idmsweep.py

@@ -37,11 +37,15 @@ def linfty(re, n=25, tmax=100, dt=0.1, timeplot=False):
 def idmgo(t, a, b, s, v, fmt=None):
     idm = IDMLin(t=t, a=a, b=b, s0=s)
     res = arange(100)/100.
-    idm.go(v, res)
+    mr = ones(size(res))
+    for i, re in enumerate(res):
+        idm.go(v, re)
+        mr[i] = idm.maxresponse()
+
     if fmt is None:
-        return res, idm.os()
+        return res, mr
     else:
-        return res, idm.os(), fmt
+        return res, mr, fmt
 
 def sweep(fn, rngs, axis=None, separate=True, logx=False):
     if logx:
@@ -77,20 +81,21 @@ if __name__ == "__main__":
     ### What happens when you vary a?
     ### ( enhance )
 
-    avals = linspace(0.5, 0.7, 20)
+    avals = linspace(0.6, 1.5, 10)
+    bvals = linspace(0.6, 1.5, 10)
     dots = ["-"] * len(avals)
 
-    arng = {"name": "a",   "color": "", "dots": dots, "default": 1.0, "values": avals }
-    trng = {"name": "T",   "default": 1.5}
-    brng = {"name": "b",   "default": 3.0}
-    srng = {"name": "s_0", "default": 2.0}
-    vrng = {"name": "v_0", "default": 30.}
+    arng = {"name": "a",   "color": "", "dots": dots, "default": 0.5, "values": avals }
+    trng = {"name": "T",   "default": 1.2}
+    brng = {"name": "b",   "color": "", "dots": dots, "default": 0.5, "values": bvals }
+    srng = {"name": "s_0", "default": 1.5}
+    vrng = {"name": "v_0", "default": 16.1}
 
     rngs = [trng, arng, brng, srng, vrng]
 
-    sweep(idmgo, rngs, axis=[0,0.1,-.05,.45])
-    go = linfty(0.5, n=50, tmax=500)
-    sweep(go, rngs)
+    sweep(idmgo, rngs)#, axis=[0,0.1,-.05,.45])
+    # go = linfty(0.5, n=50, tmax=500)
+    # sweep(go, rngs)
 
     '''
     ### What happens when you vary IDM parameters?

sstability.py

@@ -5,11 +5,12 @@ from scipy.optimize import fmin
 from idmlin import IDMLin
 
 
-def maxStringResponse(idm, n, tau = 0):
-    gammainv = humanstring(idm, n, tau)
+def maxStringResponse(fn, idm, n, **kwargs):
+    gammainv = fn(idm, n, **kwargs)
     wmax = fmin(gammainv, 10)
+    gmax = -gammainv(wmax)
     # print "wmax =",  wmax
-    return -gammainv(wmax)
+    return gmax < 1e-2 and -1 or gmax
 
 def human(idm, tau=0):
     f, g, h = idm.f, idm.g, idm.h
@@ -21,30 +22,30 @@ def human(idm, tau=0):
         return num/den
     return gamma
 
-def humanstring(idm, n, tau = 0):
+def humanstring(idm, n, tau=0):
     gamma = human(idm, tau)
     def dbgammainv(w):
         db = 20 * log10(absolute(gamma(w)))
         return -n*db
     return dbgammainv
 
-def followedrobot(idm):#, ff, fb, gf, gb, hr, taur=0, tau=0):
-    tau = 0
+def followedrobotstring(idm, n, tau=0, alpha=0.5):#, ff, fb, gf, gb, hr, taur=0, tau=0):
     taur = 0
-    ff = idm.f/2
-    fb = idm.f/2
-    gf = idm.g/2
-    gb = idm.g/2
+    fb = idm.f*alpha
+    ff = idm.f*(1-alpha)
+    gb = idm.g*alpha
+    gf = idm.g*(1-alpha)
     hr = idm.h
 
     gamma = human(idm, tau)
+    stringinv = humanstring(idm, n, tau)
     def dbgammainv(w):
         s = w * 1j
         num = ff + gf * s
         den = ff + fb + (gf + gb + hr) * s + exp(taur * s) * s * s - (fb + gb * s) * gamma(w)
         dbnum = 20 * log10(absolute(num))
         dbden = 20 * log10(absolute(den))
-        return (dbden-dbnum)
+        return (dbden-dbnum)+stringinv(w)
     return dbgammainv
 
 def teststability(idm, N):
@@ -61,49 +62,70 @@ idm = IDMLin(t=t, a=a, b=b, s0=s)
 
 nrange = logspace(1, 2.5, 20)
 rerange = logspace(-2, -0.2, 20)
+rerange = linspace(0.1, 0.4, 20)
 
-mgrid = zeros((len(nrange), len(rerange)))
-mgrid2 = zeros((len(nrange), len(rerange)))
+mgrid_h = zeros((len(nrange), len(rerange)))
+mgrid_r25 = zeros((len(nrange), len(rerange)))
+mgrid_r50 = zeros((len(nrange), len(rerange)))
+mgrid_r75 = zeros((len(nrange), len(rerange)))
 
 for i, n in enumerate(nrange):
     for j, re in enumerate(rerange):
         idm.go(v, re)
-        mgrid[i,j] = teststability(idm, n)
-        mgrid2[i,j] = maxStringResponse(idm, n)
+        mgrid_h[i,j] = maxStringResponse(humanstring, idm, n)
+        mgrid_r25[i,j] = maxStringResponse(followedrobotstring, idm, n, alpha=0.25)
+        mgrid_r50[i,j] = maxStringResponse(followedrobotstring, idm, n, alpha=0.5)
+        mgrid_r75[i,j] = maxStringResponse(followedrobotstring, idm, n, alpha=0.75)
 
-fig, axarr = plt.subplots(1,3, sharex=True, sharey=True)
+fig, axarr = plt.subplots(2,3, sharex=True, sharey=True)
 
-im = axarr[2].imshow(mgrid, cmap=plt.cm.Reds, 
+im = axarr[0,0].imshow(mgrid_h, cmap=plt.cm.bwr,
                 interpolation='none', 
                 aspect='auto',
-		vmin=0, vmax=1,
+		vmin=-1, vmax=1,
                 )
-im = axarr[1].imshow(mgrid2, cmap=plt.cm.Reds, 
+im = axarr[0,1].imshow(mgrid_r50, cmap=plt.cm.bwr, 
                 interpolation='none', 
                 aspect='auto',
-		vmin=0, vmax=1,
+		vmin=-1, vmax=1,
                 )
-im = axarr[0].imshow((mgrid2-mgrid)*100, cmap=plt.cm.Reds, 
+im = axarr[0,2].imshow(mgrid_r50-mgrid_h, cmap=plt.cm.bwr, 
                 interpolation='none', 
                 aspect='auto',
-		vmin=0, vmax=1,
+		vmin=-1, vmax=1,
+                )
+im = axarr[1,0].imshow(mgrid_r25, cmap=plt.cm.bwr, 
+                interpolation='none', 
+                aspect='auto',
+		vmin=-1, vmax=1,
+                )
+im = axarr[1,2].imshow(mgrid_r75, cmap=plt.cm.bwr, 
+                interpolation='none', 
+                aspect='auto',
+		vmin=-1, vmax=1,
                 )
 
-axarr[0].set_yticks(range(len(nrange)))
-axarr[1].set_yticks(range(len(nrange)))
-axarr[1].set_xticks(range(len(rerange)))
+axarr[0,0].set_title("Human cars only")
+axarr[0,1].set_title("Human cars followed by BCM50 robot")
+axarr[0,2].set_title("Difference")
+axarr[1,0].set_title("Human cars followed by BCM25 robot")
+axarr[1,2].set_title("Human cars followed by BCM75 robot")
 
-axarr[0].set_yticklabels(nrange)
-axarr[1].set_yticklabels(nrange)
-axarr[1].set_xticklabels(around(rerange, 2))
+axarr[0,0].set_yticks(range(0,len(nrange),4))
+axarr[0,0].set_xticks(range(0,len(rerange),4))
 
-axarr[0].set_ylabel("N")
-axarr[1].set_ylabel("N")
-axarr[1].set_xlabel("R_e")
+axarr[0,0].set_yticklabels(around(nrange[::4], 0))
+axarr[0,0].set_xticklabels(around(rerange, 2)[::4])
+
+axarr[0,0].set_ylabel("N")
+axarr[0,0].set_xlabel("R_e")
 
 fig.subplots_adjust(right=0.8)
 cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
 fig.colorbar(im, cax=cbar_ax)
 
+print nrange
+print rerange
+
 plt.show()