def crb_compare(state0, samples0, state1, samples1, crb0=None, crb1=None,<EOL>zlayer=None, xlayer=None):
s0 = state0<EOL>s1 = state1<EOL>h0 = np.array(samples0)<EOL>h1 = np.array(samples1)<EOL>slicez = zlayer or s0.image.shape[<NUM_LIT:0>]//<NUM_LIT:2><EOL>slicex = xlayer or s0.image.shape[<NUM_LIT:2>]//<NUM_LIT:2><EOL>slicer1 = np.s_[slicez,s0.pad:-s0.pad,s0.pad:-s0.pad]<EOL>slicer2 = np.s_[s0.pad:-s0.pad,s0.pad:-s0.pad,slicex]<EOL>center = (slicez, s0.image.shape[<NUM_LIT:1>]//<NUM_LIT:2>, slicex)<EOL>mu0 = h0.mean(axis=<NUM_LIT:0>)<EOL>mu1 = h1.mean(axis=<NUM_LIT:0>)<EOL>std0 = h0.std(axis=<NUM_LIT:0>)<EOL>std1 = h1.std(axis=<NUM_LIT:0>)<EOL>mask0 = (s0.state[s0.b_typ]==<NUM_LIT:1.>) & (<EOL>analyze.trim_box(s0, mu0[s0.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)))<EOL>mask1 = (s1.state[s1.b_typ]==<NUM_LIT:1.>) & (<EOL>analyze.trim_box(s1, mu1[s1.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)))<EOL>active0 = np.arange(s0.N)[mask0]<EOL>active1 = np.arange(s1.N)[mask1]<EOL>pos0 = mu0[s0.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)[active0]<EOL>pos1 = mu1[s1.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)[active1]<EOL>rad0 = mu0[s0.b_rad][active0]<EOL>rad1 = mu1[s1.b_rad][active1]<EOL>link = analyze.nearest(pos0, pos1)<EOL>dpos = pos0 - pos1[link]<EOL>drad = rad0 - rad1[link]<EOL>drift = dpos.mean(axis=<NUM_LIT:0>)<EOL>log.info('<STR_LIT>'.format(drift))<EOL>dpos -= drift<EOL>fig = pl.figure(figsize=(<NUM_LIT>,<NUM_LIT:10>))<EOL>gs0 = ImageGrid(fig, rect=[<NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>], nrows_ncols=(<NUM_LIT:2>,<NUM_LIT:3>), axes_pad=<NUM_LIT:0.1>)<EOL>lbl(gs0[<NUM_LIT:0>], '<STR_LIT:A>')<EOL>for i,slicer in enumerate([slicer1, slicer2]):<EOL><INDENT>ax_real = gs0[<NUM_LIT:3>*i+<NUM_LIT:0>]<EOL>ax_fake = gs0[<NUM_LIT:3>*i+<NUM_LIT:1>]<EOL>ax_diff = gs0[<NUM_LIT:3>*i+<NUM_LIT:2>]<EOL>diff0 = s0.get_model_image() - s0.image<EOL>diff1 = s1.get_model_image() - s1.image<EOL>a = (s0.image - s1.image)<EOL>b = (s0.get_model_image() - s1.get_model_image())<EOL>c = (diff0 - diff1)<EOL>ptp = <NUM_LIT>*max([np.abs(a).max(), np.abs(b).max(), np.abs(c).max()])<EOL>cmap = pl.cm.RdBu_r<EOL>ax_real.imshow(a[slicer], cmap=cmap, vmin=-ptp, vmax=ptp)<EOL>ax_real.set_xticks([])<EOL>ax_real.set_yticks([])<EOL>ax_fake.imshow(b[slicer], cmap=cmap, vmin=-ptp, vmax=ptp)<EOL>ax_fake.set_xticks([])<EOL>ax_fake.set_yticks([])<EOL>ax_diff.imshow(c[slicer], cmap=cmap, vmin=-ptp, vmax=ptp)<EOL>ax_diff.set_xticks([])<EOL>ax_diff.set_yticks([])<EOL>if i == <NUM_LIT:0>:<EOL><INDENT>ax_real.set_title(r"<STR_LIT>", fontsize=<NUM_LIT>)<EOL>ax_fake.set_title(r"<STR_LIT>", fontsize=<NUM_LIT>)<EOL>ax_diff.set_title(r"<STR_LIT>", fontsize=<NUM_LIT>)<EOL>ax_real.set_ylabel('<STR_LIT>')<EOL><DEDENT>else:<EOL><INDENT>ax_real.set_ylabel('<STR_LIT>')<EOL><DEDENT><DEDENT>gs1 = GridSpec(<NUM_LIT:1>,<NUM_LIT:3>, left=<NUM_LIT>, bottom=<NUM_LIT>, right=<NUM_LIT>, top=<NUM_LIT>,<EOL>wspace=<NUM_LIT>, hspace=<NUM_LIT>)<EOL>spos0 = std0[s0.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)[active0]<EOL>spos1 = std1[s1.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)[active1]<EOL>srad0 = std0[s0.b_rad][active0]<EOL>srad1 = std1[s1.b_rad][active1]<EOL>def hist(ax, vals, bins, *args, **kwargs):<EOL><INDENT>y,x = np.histogram(vals, bins=bins)<EOL>x = (x[<NUM_LIT:1>:] + x[:-<NUM_LIT:1>])/<NUM_LIT:2><EOL>y /= len(vals)<EOL>ax.plot(x,y, *args, **kwargs)<EOL><DEDENT>def pp(ind, tarr, tsim, tcrb, var='<STR_LIT:x>'):<EOL><INDENT>bins = <NUM_LIT:10>**np.linspace(-<NUM_LIT:3>, <NUM_LIT:0.0>, <NUM_LIT:30>)<EOL>bin2 = <NUM_LIT:10>**np.linspace(-<NUM_LIT:3>, <NUM_LIT:0.0>, <NUM_LIT:100>)<EOL>bins = np.linspace(<NUM_LIT:0.0>, <NUM_LIT>, <NUM_LIT:30>)<EOL>bin2 = np.linspace(<NUM_LIT:0.0>, <NUM_LIT>, <NUM_LIT:100>)<EOL>xlim = (<NUM_LIT:0.0>, <NUM_LIT>)<EOL>ylim = (<NUM_LIT>, <NUM_LIT:30>)<EOL>ticks = ticker.FuncFormatter(lambda x, pos: '<STR_LIT>'.format(np.log10(x)))<EOL>scaler = lambda x: x <EOL>ax_crb = pl.subplot(gs1[<NUM_LIT:0>,ind])<EOL>ax_crb.hist(scaler(np.abs(tarr)), bins=bins,<EOL>normed=True, alpha=<NUM_LIT>, histtype='<STR_LIT>', lw=<NUM_LIT:1>)<EOL>ax_crb.hist(scaler(np.abs(tcrb)).ravel(), bins=bin2,<EOL>normed=True, alpha=<NUM_LIT:1.0>, histtype='<STR_LIT>', ls='<STR_LIT>', lw=<NUM_LIT>, color='<STR_LIT:k>')<EOL>ax_crb.hist(scaler(np.abs(tsim).ravel()), bins=bin2,<EOL>normed=True, alpha=<NUM_LIT:1.0>, histtype='<STR_LIT>', lw=<NUM_LIT:3>)<EOL>ax_crb.set_xlabel(r"<STR_LIT>" % (var,var), fontsize=<NUM_LIT>)<EOL>ax_crb.set_xlim(xlim)<EOL>ax_crb.grid(b=False, which='<STR_LIT>', axis='<STR_LIT>')<EOL>if ind == <NUM_LIT:0>:<EOL><INDENT>lbl(ax_crb, '<STR_LIT:B>')<EOL>ax_crb.set_ylabel(r"<STR_LIT>")<EOL><DEDENT>else:<EOL><INDENT>ax_crb.set_yticks([])<EOL><DEDENT>ax_crb.locator_params(axis='<STR_LIT:x>', nbins=<NUM_LIT:3>)<EOL><DEDENT>f,g = <NUM_LIT>, <NUM_LIT><EOL>sim = f*sim_crb_diff(spos0[:,<NUM_LIT:1>], spos1[:,<NUM_LIT:1>][link])<EOL>crb = g*sim_crb_diff(crb0[<NUM_LIT:0>][:,<NUM_LIT:1>][active0], crb1[<NUM_LIT:0>][:,<NUM_LIT:1>][active1][link])<EOL>pp(<NUM_LIT:0>, dpos[:,<NUM_LIT:1>], sim, crb, '<STR_LIT:x>')<EOL>sim = f*sim_crb_diff(spos0[:,<NUM_LIT:0>], spos1[:,<NUM_LIT:0>][link])<EOL>crb = g*sim_crb_diff(crb0[<NUM_LIT:0>][:,<NUM_LIT:0>][active0], crb1[<NUM_LIT:0>][:,<NUM_LIT:0>][active1][link])<EOL>pp(<NUM_LIT:1>, dpos[:,<NUM_LIT:0>], sim, crb, '<STR_LIT:z>')<EOL>sim = f*sim_crb_diff(srad0, srad1[link])<EOL>crb = g*sim_crb_diff(crb0[<NUM_LIT:1>][active0], crb1[<NUM_LIT:1>][active1][link])<EOL>pp(<NUM_LIT:2>, drad, sim, crb, '<STR_LIT:a>')<EOL>gs2 = GridSpec(<NUM_LIT:2>,<NUM_LIT:2>, left=<NUM_LIT>, bottom=<NUM_LIT>, right=<NUM_LIT>, top=<NUM_LIT>,<EOL>wspace=<NUM_LIT>, hspace=<NUM_LIT>)<EOL>ax_hist = pl.subplot(gs2[<NUM_LIT:0>,<NUM_LIT:0>])<EOL>ax_hist.hist(std0[s0.b_pos], bins=np.logspace(-<NUM_LIT>, <NUM_LIT:0>, <NUM_LIT:50>), alpha=<NUM_LIT>, label='<STR_LIT>', histtype='<STR_LIT>')<EOL>ax_hist.hist(std0[s0.b_rad], bins=np.logspace(-<NUM_LIT>, <NUM_LIT:0>, <NUM_LIT:50>), alpha=<NUM_LIT>, label='<STR_LIT>', histtype='<STR_LIT>')<EOL>ax_hist.set_xlim((<NUM_LIT:10>**-<NUM_LIT>, <NUM_LIT:1>))<EOL>ax_hist.semilogx()<EOL>ax_hist.set_xlabel(r"<STR_LIT>")<EOL>ax_hist.set_ylabel(r"<STR_LIT>")<EOL>ax_hist.legend(loc='<STR_LIT>')<EOL>lbl(ax_hist, '<STR_LIT:C>')<EOL>imdiff = ((s0.get_model_image() - s0.image)/s0._sigma_field)[s0.image_mask==<NUM_LIT:1.>].ravel()<EOL>mu = imdiff.mean()<EOL>x = np.linspace(-<NUM_LIT:5>,<NUM_LIT:5>,<NUM_LIT>)<EOL>ax_diff = pl.subplot(gs2[<NUM_LIT:0>,<NUM_LIT:1>])<EOL>ax_diff.plot(x, <NUM_LIT:1.0>/np.sqrt(<NUM_LIT:2>*np.pi) * np.exp(-(x-mu)**<NUM_LIT:2> / <NUM_LIT:2>), '<STR_LIT:->', alpha=<NUM_LIT>, color='<STR_LIT:k>', lw=<NUM_LIT:2>)<EOL>ax_diff.hist(imdiff, bins=<NUM_LIT:1000>, histtype='<STR_LIT>', alpha=<NUM_LIT>, normed=True)<EOL>ax_diff.semilogy()<EOL>ax_diff.set_ylabel(r"<STR_LIT>")<EOL>ax_diff.set_xlabel(r"<STR_LIT>")<EOL>ax_diff.locator_params(axis='<STR_LIT:x>', nbins=<NUM_LIT:5>)<EOL>ax_diff.grid(b=False, which='<STR_LIT>', axis='<STR_LIT:y>')<EOL>ax_diff.set_xlim(-<NUM_LIT:5>, <NUM_LIT:5>)<EOL>ax_diff.set_ylim(<NUM_LIT>, <NUM_LIT>)<EOL>lbl(ax_diff, '<STR_LIT:D>')<EOL>pos = mu0[s0.b_pos].reshape(-<NUM_LIT:1>,<NUM_LIT:3>)<EOL>rad = mu0[s0.b_rad]<EOL>mask = analyze.trim_box(s0, pos)<EOL>pos = pos[mask]<EOL>rad = rad[mask]<EOL>gx, gy = analyze.gofr(pos, rad, mu0[s0.b_zscale][<NUM_LIT:0>], resolution=<NUM_LIT>,mask_start=<NUM_LIT:0.5>)<EOL>mask = gx < <NUM_LIT:5><EOL>gx = gx[mask]<EOL>gy = gy[mask]<EOL>ax_gofr = pl.subplot(gs2[<NUM_LIT:1>,<NUM_LIT:0>])<EOL>ax_gofr.plot(gx, gy, '<STR_LIT:->', lw=<NUM_LIT:1>)<EOL>ax_gofr.set_xlabel(r"<STR_LIT>")<EOL>ax_gofr.set_ylabel(r"<STR_LIT>")<EOL>ax_gofr.locator_params(axis='<STR_LIT>', nbins=<NUM_LIT:5>)<EOL>lbl(ax_gofr, '<STR_LIT:E>')<EOL>gx, gy = analyze.gofr(pos, rad, mu0[s0.b_zscale][<NUM_LIT:0>], method='<STR_LIT>')<EOL>mask = gx < <NUM_LIT:5><EOL>gx = gx[mask]<EOL>gy = gy[mask]<EOL>gy[gy <= <NUM_LIT:0.>] = gy[gy><NUM_LIT:0>].min()<EOL>ax_gofrs = pl.subplot(gs2[<NUM_LIT:1>,<NUM_LIT:1>])<EOL>ax_gofrs.plot(gx, gy, '<STR_LIT:->', lw=<NUM_LIT:1>)<EOL>ax_gofrs.set_xlabel(r"<STR_LIT>")<EOL>ax_gofrs.set_ylabel(r"<STR_LIT>")<EOL>ax_gofrs.locator_params(axis='<STR_LIT>', nbins=<NUM_LIT:5>)<EOL>ax_gofrs.grid(b=False, which='<STR_LIT>', axis='<STR_LIT:y>')<EOL>lbl(ax_gofrs, '<STR_LIT:F>')<EOL>ylim = ax_gofrs.get_ylim()<EOL>ax_gofrs.set_ylim(gy.min(), ylim[<NUM_LIT:1>])<EOL>
To run, do:

s,h = pickle...
s1,h1 = pickle...
    i.e. /media/scratch/bamf/vacancy/vacancy_zoom-1.tif_t002.tif-featured-v2.pkl
    i.e. /media/scratch/bamf/frozen-particles/0.tif-featured-full.pkl
crb0 = diag_crb_particles(s); crb1 = diag_crb_particles(s1)
crb_compare(s,h[-25:],s1,h1[-25:], crb0, crb1)