add chap6 & minor changes

fig/plottingDistriInOptimizedArrangement-gen.py

@@ -0,0 +1,161 @@
+import math
+import numpy as np
+import numpy.matlib
+import numpy.linalg
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import MultipleLocator
+from matplotlib import cm
+import scipy
+import scipy.integrate as integrate
+import scipy.optimize as optimize
+from scipy import special
+from scipy import signal
+
+sigmaF=0.10895354*(2**0.5)
+sigmaA=0.078188270041128
+# sigmaA2=sigmaA*(2**0.5)
+gamma=0.5488594815581366
+# thres=0.103468491232405
+
+OX=lambda x:(1.6971354708*x**2 + -3.3520704577*x + 1.6552008479)
+INTOX=lambda x:(1.6971354708*x**3/3 + -3.3520704577*x**2/2 + 1.6552008479*x)
+OXint=integrate.quad(OX,0,1)[0]
+X=lambda x:(OX(x)/OXint)
+INTX=lambda x:(INTOX(x)/OXint)
+
+def R(x,sigma=sigmaA):
+    return (special.erf(x/(2**0.5*sigma))+1)/2
+
+def omR(x,sigma=sigmaA):
+    return (1-special.erf(x/(2**0.5*sigma)))/2
+
+def rangeR(l,r,sigma=sigmaA):
+    return 1-R(l,sigma)-omR(r,sigma)
+
+def P(x,sigma=sigmaA):
+    return np.exp(-(x**2)/(2*sigma**2))/((2*math.pi)**0.5*sigma)
+
+def C(x):
+    return integrate.quad((lambda s:(X(s)*P(x-s)/rangeR(0-s,1-s))),0,1)[0]
+
+cP=256
+rP=1/cP
+
+def P2(x):
+    return rangeR(x-rP/2,x+rP/2)
+
+def getvec_x(x,tss): # prob,probx,sum
+    allprob=rangeR(0-x,1-x)
+    once_x=np.arange(rP/2,1,rP,dtype=np.float64)
+    once_h=P2(once_x-x)/allprob
+#     print(sum(once_h))
+    assert(abs(sum(once_h)-1)<1e-8)
+    cur_h=once_h.copy()
+    ret_prob=0
+    ret_sum=0
+#     once_h2=signal.convolve(once_h,[0.5,0.5])
+    for i in range(len(tss)):
+        ncP=cP*(i+1)
+        nrP=1/ncP
+        cur_x=np.arange(nrP/2,1,nrP,dtype=np.float64)
+        assert(len(cur_h)==len(cur_x))
+        t=tss[i]
+        ret_sum+=np.sum(cur_h)
+        d=np.searchsorted(cur_x,t-nrP/2,side='right')
+        cur_h[:d]=0
+        if d<len(cur_h):
+            cur_h[d]*=((cur_x[d]+nrP/2)-t)/nrP
+        if i<len(tss)-1:
+            cur_h=signal.convolve(cur_h,once_h)
+            cur_h=(np.append([0],cur_h)+np.append(cur_h,[0]))/2
+#             cur_h=cur_h[0::2]+(np.append(cur_h[1::2],[0])+np.append([0],cur_h[1::2]))/2
+#             cur_h=(curh[0::2]+np.append([0],curh[0::2]))/2+(curh[1::2]+6*np.append([0],curh[1::2])+np.append([0,0],curh[1::2]))/8
+#             cur_h=(np.append(signal.convolve(cur_h[0::2],once_h[0::2]),[0])+np.append([0],signal.convolve(cur_h[1::2],once_h[1::2])))*2
+    ret_prob=np.sum(cur_h)
+    return np.array([ret_prob,ret_prob*x,ret_sum])
+
+# mp=dict([])
+
+def getvec(tss):
+#     print(tss[2])
+#     if tuple(tss) in mp:
+#         return mp[tuple(tss)]
+    ret=integrate.quad_vec((lambda x:getvec_x(x,tss)*X(x)),0,1,epsabs=1e-8,epsrel=1e-8)[0]
+#     mp[tuple(tss)]=ret
+    return ret
+
+def totprob(tss):
+    r=getvec(tss)
+    return r[0]
+
+def ex(tss):
+    r=getvec(tss)
+    return r[1]/r[0]
+
+def cost(tss):
+    r=getvec(tss)
+    return r[2]
+
+mp=dict([])
+
+lastrt=0
+
+def getvec2(tss,show=False):
+    global lastrt
+    if (tuple(tss) in mp) and not show:
+        return mp[tuple(tss)]
+    tss=list(tss)
+    if totprob(tss+[0])<0.1+1e-9:
+        mp[tuple(tss)]=getvec(tss+[0])
+        return mp[tuple(tss)]
+    tmp=scipy.optimize.root_scalar(lambda x:totprob(tss+[x])-0.1, bracket=[0,1], method='brentq',xtol=1e-9)
+    rt=tmp.root
+    if show:
+        lastrt=rt
+        print(tss+[rt])
+    mp[tuple(tss)]=getvec(tss+[rt])
+    return mp[tuple(tss)]
+
+def totprob2(tss):
+    r=getvec2(tss)
+    return r[0]
+
+def ex2(tss):
+    r=getvec2(tss)
+    return r[1]/r[0]
+
+def cost2(tss):
+    r=getvec2(tss)
+    return r[2]
+
+plt.rcParams['figure.figsize'] = (8.0, 4.0)
+plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
+plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
+plt.rcParams['text.usetex']=False
+
+ax2=plt.subplot(1,1,1)
+ax=ax2.twinx()
+xs=np.arange(0,1,0.01)
+xx=[0.11394000787380809, 0.3616143710651146, 0.5336339315602933]
+ys=[getvec_x(cx,xx)[0] for cx in xs]
+print(sum(xs[i]*ys[i]*X(xs[i]) for i in range(len(xs)))/sum(ys[i]*X(xs[i]) for i in range(len(xs))))
+ax.plot(xs,ys,color='blue',label='$P(x)$的图像',linestyle='--')
+# ax.vlines(ex(xx),-10,10,colors='orange',linestyles='dotted',linewidths=1.5)
+# at=scipy.optimize.root_scalar(lambda x:INTX(x)-0.9, bracket=[0,1], method='brentq',xtol=1e-9).root
+# mnv=scipy.integrate.quad(lambda x:X(x)*x,at,1)[0]/0.1
+# ax.vlines(mnv,-10,10,colors='red',linestyles='dotted',linewidths=1.5)
+ax.set_ylim(0-0.05,1+0.05)
+ax.set_xlim(0,1)
+ax.set_ylabel('通过概率')
+ax.legend(loc=(0.8,0.3))
+ax2.plot(xs,X(xs),color='red',label='$\\mathcal{X}_{A_1}(x)$的图像')
+ax2.plot(xs,[X(cx)*getvec_x(cx,xx)[0] for cx in xs],color='orange',label='$P(x)\\cdot\\mathcal{X}_{A_1}(x)$的图像')
+ax2.set_ylabel('分布密度')
+ax2.set_xlabel('期望得分')
+ax2.legend(loc=(0.01,0.3))
+# getvec2([0],show=True)
+# xx=[0]+[lastrt]
+# ys=[getvec_x(cx,xx)[0] for cx in xs]
+# print(sum(xs[i]*ys[i]*X(xs[i]) for i in range(len(xs)))/sum(ys[i]*X(xs[i]) for i in range(len(xs))))
+# plt.plot(xs,ys)
+plt.savefig(fname="plottingDistriInOptimizedArrangement.pdf",format="pdf",bbox_inches='tight',pad_inches=0.05)

fig/plottingH1H2H3-gen.py

@@ -0,0 +1,116 @@
+costs=[1.9854540490170365,1.8897541219398266,1.783125960258785,1.697661177237779,1.5833520920472837,1.499519431563351,1.3999052841943496,1.299999325270233]
+x1s=[0.11394000787380809,0.1552115117282269,0.20341885780262836,0.23489489080820136,0.2850241829465864,0.3243740272124554,0.38111643418482016,0.4486791582817367,]
+x2s=[0.3616143710651146,0.35727851169949304,0.3588230683776886,0.3796779340021746,0.4059465485621914,0.43086125275598436,0.4589843750741278,0.49046097762175916,]
+x3s=[0.5336339315602933,0.533633931821313,0.5336339304450538,0.5336339101399719,0.5336330310173099,0.5336205377222956,0.5333756468004085,0.5297378279288261,]
+
+import math
+import numpy as np
+import numpy.matlib
+import numpy.linalg
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import MultipleLocator
+from matplotlib import cm
+import scipy
+import scipy.integrate as integrate
+import scipy.optimize as optimize
+from scipy import special
+from scipy import signal
+
+sigmaF=0.10895354*(2**0.5)
+sigmaA=0.078188270041128
+# sigmaA2=sigmaA*(2**0.5)
+gamma=0.5488594815581366
+# thres=0.103468491232405
+
+OX=lambda x:(1.6971354708*x**2 + -3.3520704577*x + 1.6552008479)
+OXint=integrate.quad(OX,0,1)[0]
+X=lambda x:(OX(x)/OXint)
+
+def R(x,sigma=sigmaA):
+    return (special.erf(x/(2**0.5*sigma))+1)/2
+
+def omR(x,sigma=sigmaA):
+    return (1-special.erf(x/(2**0.5*sigma)))/2
+
+def rangeR(l,r,sigma=sigmaA):
+    return 1-R(l,sigma)-omR(r,sigma)
+
+def P(x,sigma=sigmaA):
+    return np.exp(-(x**2)/(2*sigma**2))/((2*math.pi)**0.5*sigma)
+
+def C(x):
+    return integrate.quad((lambda s:(X(s)*P(x-s)/rangeR(0-s,1-s))),0,1)[0]
+
+cP=256
+rP=1/cP
+
+def P2(x):
+    return rangeR(x-rP/2,x+rP/2)
+
+def getvec_x(x,tss):
+    allprob=rangeR(0-x,1-x)
+    once_x=np.arange(rP/2,1,rP,dtype=np.float64)
+    once_h=P2(once_x-x)/allprob
+    assert(abs(sum(once_h)-1)<1e-8)
+    cur_h=once_h.copy()
+    ret=[]
+    for i in range(len(tss)):
+        ncP=cP*(i+1)
+        nrP=1/ncP
+        cur_x=np.arange(nrP/2,1,nrP,dtype=np.float64)
+        assert(len(cur_h)==len(cur_x))
+        t=tss[i]
+        d=np.searchsorted(cur_x,t-nrP/2,side='right')
+        cur_h[:d]=0
+        ret.append(sum(cur_h))
+        if d<len(cur_h):
+            cur_h[d]*=((cur_x[d]+nrP/2)-t)/nrP
+        if i<len(tss)-1:
+            cur_h=signal.convolve(cur_h,once_h)
+            cur_h=(np.append([0],cur_h)+np.append(cur_h,[0]))/2
+    return np.array(ret)
+
+# mp=dict([])
+
+def getvec(tss):
+#     print(tss[2])
+#     if tuple(tss) in mp:
+#         return mp[tuple(tss)]
+    ret=integrate.quad_vec((lambda x:getvec_x(x,tss)*X(x)),0,1,epsabs=1e-8,epsrel=1e-8)[0]
+#     mp[tuple(tss)]=ret
+    return ret
+
+plt.rcParams['figure.figsize'] = (9.0, 9.0)
+plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
+plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
+plt.rcParams['text.usetex']=False
+
+plta=plt.subplot(2,1,1)
+pltb=plt.subplot(2,1,2)
+
+plta.plot(costs,x1s,'s-',label='第1轮分数线',color='darkblue')#,markerfacecolor='none')
+plta.plot(costs,x2s,'s-',label='第2轮分数线',color='blue')#,markerfacecolor='none')
+plta.plot(costs,x3s,'s-',label='第3轮分数线',color='cornflowerblue')#,markerfacecolor='none')
+tmp=[getvec([x1s[i],x2s[i],x3s[i]]) for i in range(len(costs))]
+tmp1=[f[0] for f in tmp]
+tmp2=[f[1] for f in tmp]
+tmp3=[f[2] for f in tmp]
+pltb.plot(costs,tmp1,'^-',label='第1轮选拔人数',color='darkgreen')#,markerfacecolor='none')
+pltb.plot(costs,tmp2,'^-',label='第2轮选拔人数',color='forestgreen')#,markerfacecolor='none')
+pltb.plot(costs,tmp3,'^-',label='第3轮选拔人数',color='limegreen')#,markerfacecolor='none')
+# plt.vlines(costs,[min(x1s[i],tmp3[i]) for i in range(len(x1s))],[max(x3s[i],tmp1[i]) for i in range(len(x1s))],linestyles='dotted',colors='orange',linewidths=1.5)
+# plt.legend(aa,['x1','x2','x3'],loc=0)
+# plt.legend(bb,['r1','r2','r3'],loc=0)
+plta.vlines(costs,x1s,x3s,linestyles='dotted',colors='orange',linewidths=1.5)
+plta.legend(loc=0)
+plta.grid(linestyle='--')
+plta.set_xlabel('开销')
+plta.set_ylabel('分数')
+plta.set_ylim((0,None))
+pltb.vlines(costs,tmp3,tmp1,linestyles='dotted',colors='orange',linewidths=1.5)
+pltb.legend(loc=0)
+pltb.grid(linestyle='--')
+pltb.set_xlabel('开销')
+pltb.set_ylabel('人数（占初始人数的比例）')
+pltb.set_ylim((0,None))
+plt.savefig(fname="plottingH1H2H3.pdf",format="pdf",bbox_inches='tight',pad_inches=0.05)