<jupyter_start><jupyter_text>Now do a cube<jupyter_code>results = []
uncertainties = []

# myerr = np.ones_like(mydata)*0.01*(onp.nanmedian(np.abs(mydata)))
seps_in = np.linspace(90,500,20)
calib = np.median(calibrators,axis=0)

for j, sep_in in enumerate(tqdm(seps_in)):
#     img = preprocess_like(tgt_cube[:,:,j],b)
    
    truth = [sep_in,90,25]
    binary = sim_binary(cal,*truth,43.1)
    img_sim = preprocess_like(binary,b)
    a = c.copy()
    
    a.extract_KPD_single_frame(
    img_sim, pscale, wl, recenter=True,method='LDFT1')

    a.kpi.name = "simulation" #  # labels the data
#     calib  = pca.inverse_transform(pca.transform(a.KPDT[0]))[0]
    mydata = a.KPDT[0][0]-calib
    
    p0 = truth
    p0[1] = np.mod(p0[1],360.)

    # mfit = a.binary_model_fit(p0)
    mfit = ben_binary_model_fit(p0,kpo=a)
    
    cvis_b = xara.core.cvis_binary(
        a.kpi.UVC[:,0], a.kpi.UVC[:,1], wl, p2) # binary
    ker_theo = a.kpi.KPM.dot(np.angle(cvis_b))

    if myerr is not None:
        chi2 = np.sum(((mydata - ker_theo)/(this_error))**2) / a.kpi.nbkp
    else:
        chi2 = np.sum(((mydata - ker_theo))**2) / a.kpi.nbkp

    added_error = add_error()
    this_error = np.sqrt(myerr**2+added_error**2)

    p1 = mfit[0] # the best fit parameter vector (sep, P.A., contrast)
    p1[1] = np.mod(p1[1],360.)
    results.append(p1)

    def negative_log_posterior(params):
        temp = binary_model(params,a)
        model = np.array(a.kpi.KPM).dot(np.angle(temp))
        error = mydata-model
        error /= (this_error)
        chi2 = np.sum((error)**2)
        return chi2

    def approx_covariance_matrix(params):
        # evaluate the covariance matrix of the approximate normal
        return np.linalg.inv(hessian(negative_log_posterior)(params))

    covar = approx_covariance_matrix(1.0*np.array(p1))
    uncertainties.append(np.sqrt(np.diag(covar)))

seps_out = np.array([result[0] for result in results])
thetas_out = np.array([result[1] for result in results])
cons_out = np.array([result[2] for result in results])

dseps_out = np.array([uncertainty[0] for uncertainty in uncertainties])
dthetas_out = np.array([uncertainty[1] for uncertainty in uncertainties])
dcons_out = np.array([uncertainty[2] for uncertainty in uncertainties])
pscale
fig, (ax1,ax2,ax3) = plt.subplots(1,3,figsize=(16.0,4.0))
inds = np.arange(len(seps_in))

ax1.plot(seps_in,seps_out-seps_in,'.')
ax1.axhline(0,color='k',linestyle='--',alpha=0.5)
ax1.errorbar(seps_in,seps_out-seps_in,yerr=dseps_out,ls='none',color=colours[0])

# ax1.plot(seps_in,seps_in,'--k',alpha=0.5)

ax2.plot(seps_in,thetas_out,'.')
ax2.axhline(truth[1],color='k',linestyle='--',alpha=0.5)
ax2.errorbar(seps_in,thetas_out,yerr=dthetas_out,ls='none',color=colours[0])

ax3.plot(seps_in,cons_out,'.')
ax3.errorbar(seps_in,dcons_out,yerr=dcons_out,ls='none',color=colours[0])

ax3.axhline(truth[2],color='k',linestyle='--',alpha=0.5)<jupyter_output><empty_output><jupyter_text>Now try Peter's<jupyter_code>myerr = np.ones_like(mydata)*0.01*(onp.nanmedian(np.abs(mydata)))
myerr = np.std(calibrators,axis=0)
calib = np.median(calibrators,axis=0)
chi2s = []

for do_cal in [False,True]:
    results = []
    uncertainties = []
    ssteps, thsteps = [], []
    # plt.plot(0,0,'*')
    for i in range(0,20):
        step=np.array((-0.17,0.28))*i + np.array((3.0,2.0))
    #     print(step)
        sep_step = 43.1 * (np.hypot(*step))
        th_step = 180./np.pi*(np.angle(step[1]-1.j*step[0]))
    #     plt.plot(step[0],step[1],'.')

        ssteps.append(sep_step)
        thsteps.append(th_step)
    ssteps = np.array(ssteps)
    thsteps = np.mod(np.array(thsteps),360)
    # plt.axis('square')

    for j in tqdm(range(20)):
        img = preprocess_like(tgt_cube[:,:,j],b)

        a = c.copy()

        a.extract_KPD_single_frame(
        img, pscale, wl, recenter=True,method='LDFT1')

        a.kpi.name = "simulation" #  # labels the data
        if do_cal:
            mydata = a.KPDT[0][0]-calib#b.KPDT[0][0]
        else:
            mydata = a.KPDT[0][0]

    #     print("\ncomputing colinearity map...")
    #     gsize = 100 # gsize x gsize grid
    #     gstep = 5 # grid step in mas
    #     xx, yy = np.meshgrid(
    #             np.arange(gsize) - gsize/2, np.arange(gsize) - gsize/2)
    #     azim = -np.arctan2(xx, yy) * 180.0 / np.pi
    #     dist = np.hypot(xx, yy) * gstep

    # #     mmap = kpo1.kpd_binary_match_map(100, 10, mydata/myerr, norm=True)
    #     mmap = a.kpd_binary_match_map(gsize, gstep, mydata, norm=True)
    #     x0, y0 = np.argmax(mmap) % gsize, np.argmax(mmap) // gsize
    #     print("max colinearity found for sep = %.2f mas and ang = %.2f deg" % (
    #             dist[y0, x0], azim[y0, x0]))

    #     p0 =  [dist[y0, x0], azim[y0, x0], mmap.max()] # good starting point
        p0 = [ssteps[j], thsteps[j],20]
        p0[1] = np.mod(p0[1],360.)
    #     p0[0] = 150+j*10

        # mfit = a.binary_model_fit(p0)
        mfit = ben_binary_model_fit(p0,kpo=a)

        p1 = mfit[0] # the best fit parameter vector (sep, P.A., contrast)
        p1[1] = np.mod(p1[1],360.)
        p1[0] = np.abs(p1[0])
        results.append(p1)

        cvis_b = xara.core.cvis_binary(a.kpi.UVC[:,0], a.kpi.UVC[:,1], wl, p1) # binary
        ker_theo = a.kpi.KPM.dot(np.angle(cvis_b))

        def get_chi2(addederror):
            return np.sum(((mydata - np.array(ker_theo))/np.sqrt(myerr**2+addederror**2))**2) / (1.0*a.kpi.nbkp)
        chi2 = get_chi2(0)
        chi2s.append(chi2)

        if chi2 > 1.0:
            from jax import grad

            def add_error():
                def get_objective(addederror):
                    return np.abs(get_chi2(addederror)-1.)
                jac = grad(get_objective)
                return minimize(get_objective,0.0025,method='BFGS')['x']

            added_error = add_error()
            this_error = np.sqrt(myerr**2+added_error**2)
        else:
            this_error = myerr 
        
        def negative_log_posterior(params):
            temp = binary_model(params,a)
            model = np.array(a.kpi.KPM).dot(np.angle(temp))
            error = mydata-model
            error /= (this_error)
            chi2 = np.sum((error)**2)

            return chi2

        def grad_negative_log_posterior(params):
            return -1*jacrev(negative_log_posterior)(params)

        def approx_covariance_matrix(params):
            # evaluate the covariance matrix of the approximate normal
            return np.linalg.inv(hessian(negative_log_posterior)(params))

        covar = approx_covariance_matrix(1.0*np.array(p1))
        uncertainties.append(np.sqrt(np.diag(covar)))
    print('chi2s:',(np.array(chi2s)))
    
    if do_cal:
        seps_out_nocal = np.array([result[0] for result in results])
        thetas_out_nocal = np.array([result[1] for result in results])
        cons_out_nocal = np.array([result[2] for result in results])
        
        dseps_out_nocal = np.array([uncertainty[0] for uncertainty in uncertainties])
        dthetas_out_nocal = np.array([uncertainty[1] for uncertainty in uncertainties])
        dcons_out_nocal = np.array([uncertainty[2] for uncertainty in uncertainties])

    else:
        seps_out = np.array([result[0] for result in results])
        thetas_out = np.array([result[1] for result in results])
        cons_out = np.array([result[2] for result in results])
        dseps_out = np.array([uncertainty[0] for uncertainty in uncertainties])
        dthetas_out = np.array([uncertainty[1] for uncertainty in uncertainties])
        dcons_out = np.array([uncertainty[2] for uncertainty in uncertainties])


fig, (ax1,ax2,ax3) = plt.subplots(3,1,figsize=(6,9),sharex=True)
inds = np.arange(len(ssteps))
ax1.plot(ssteps,'.',label='Injected')
ax1.plot(seps_out,'.',label='Calibrated',color=colours[1])
ax1.errorbar(inds,seps_out,yerr=dseps_out,ls='none',color=colours[1])
ax1.plot(seps_out_nocal,'.',label='Uncalibrated')
ax1.errorbar(inds,seps_out_nocal,yerr=dseps_out_nocal,ls='none',color=colours[2])
ax1.set_ylabel('Separation (mas)')
# ax1.plot(seps_in,seps_in,'--k',alpha=0.5)
ax1.legend(fontsize=16,frameon=False)
ax1.set_xticks([])

ax2.plot(thsteps,'.')
# ax2.plot(thetas_out,'.')
# ax2.plot(thetas_out_nocal,'.')
ax2.plot(thetas_out,'.',label='Calibrated',color=colours[1])
ax2.errorbar(inds,thetas_out,yerr=dthetas_out,ls='none',color=colours[1])
ax2.plot(thetas_out_nocal,'.',label='Uncalibrated')
ax2.errorbar(inds,thetas_out_nocal,yerr=dthetas_out_nocal,ls='none',color=colours[2])

ax2.set_ylabel('Position Angle (deg)')
ax2.set_xticks([])

# ax2.axhline(truth[1],color='k',linestyle='--',alpha=0.5)
# ax3.plot(cons_out,'.',color=colours[1])
# ax3.plot(cons_out_nocal,'.',color=colours[2])
ax3.plot(cons_out,'.',label='Calibrated',color=colours[1])
ax3.errorbar(inds,cons_out,yerr=dcons_out,ls='none',color=colours[1])
ax3.plot(cons_out_nocal,'.',label='Uncalibrated')
ax3.errorbar(inds,cons_out_nocal,yerr=dcons_out_nocal,ls='none',color=colours[2])
ax3.axhline(20)
ax3.set_ylabel('Contrast')
ax3.set_xlabel('Epoch')
ax2.set_xticks(np.arange(20)[::5])

plt.subplots_adjust(wspace=0, hspace=0)

# ax3.axhline(truth[2],color='k',linestyle='--',alpha=0.5)
# plt.savefig('hare_hounds.png',bbox_inches='tight')
# paramlimits=[p1[0]*0.9,p1[0]*1.1,0,360,15,30.]

# def binary_model(params,kpo):
#     u = kpo.kpi.UVC[:,0]
#     v = kpo.kpi.UVC[:,1]
#     wl = kpo.CWAVEL
#     detpa = 0
#     return(xara.core.cvis_binary_jax(u,v,wl, params, detpa))

# @jit
# def kp_loglikelihood(params,kpo):
#     temp = binary_model(params,kpo)
#     model = np.array(kpo.kpi.KPM).dot(np.angle(temp))
#     error = mydata-model
#     error /= (myerr)
#     chi2 = -np.sum((error)**2)
#     return chi2 

# def negative_log_posterior(params):
#     temp = binary_model(params,a)
#     model = np.array(a.kpi.KPM).dot(np.angle(temp))
#     error = mydata-model
#     error /= (myerr)
#     chi2 = np.sum((error)**2)

#     return chi2

# def grad_negative_log_posterior(params):
#     return -1*jacrev(negative_log_posterior)(params)

# @jit
# def approx_covariance_matrix(params):
#     # evaluate the covariance matrix of the approximate normal
#     return np.linalg.inv(hessian(negative_log_posterior)(params))


# def lnprior(params):
#     if paramlimits[0] < params[0] < paramlimits[1] and paramlimits[2] < params[1] < paramlimits[3] and paramlimits[4] < params[2] < paramlimits[5]:
#         return -np.log(params[0]) -np.log(params[2])
#     return -np.inf

# def lnprob(params,kpo):
#     return lnprior(params) + kp_loglikelihood(params,kpo)

# ndim=3
# nwalkers=100
# plot=False
# burnin=100
# nsteps=1000

# import time

# np.sqrt(np.diag(approx_covariance_matrix(1.0*np.array(p1))))
# import emcee

# ivar = np.array(p1)  # initial parameters for model-fit

# ball = np.array([ivar + 0.01*ivar*onp.random.rand(ndim) for i in range(nwalkers)]) # initialise walkers in a ball

# print('Running emcee now!')

# t0 = time.time()

# sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[a])

# # burn in
# pos,prob,state = sampler.run_mcmc(ball, burnin)
# sampler.reset()

# t1 = time.time()

# print('Burnt in! Took %.3f seconds' %(t1-t0))

# # restart
# sampler.run_mcmc(pos,nsteps)

# tf = time.time()

# print('Time elapsed = %.3f s' %(tf-t0))

# seps = sampler.flatchain[:,0]
# ths = sampler.flatchain[:,1]

# meansep = np.mean(seps)
# dsep = np.std(seps)

# meanth = np.mean(ths)
# dth = np.std(ths)

# cs = sampler.flatchain[:,2]
# bestcon = np.mean(cs)
# conerr = np.std(cs)

# print('Separation %.3f pm %.3f mas' % (meansep,dsep))
# print('Position angle %.3f pm %.3f deg' % (meanth,dth))
# print('Contrast at',wl,'um %.3f pm %.3f' % (bestcon,conerr))

# import corner
# # Plot it.
# figure = corner.corner(sampler.flatchain, labels=[r"$\rho$", r"$\theta$", r"$c$",],
#                        quantiles=[0.16, 0.5, 0.84],
#                        show_titles=True, title_kwargs={"fontsize": 12},truths=p0)

p0<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17199099174243057, 'fraction_numerical': 0.03336391692384686, 'mean_word_length': 3.5109104589917233, 'pattern_counts': {'":': 1, '<': 15, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '17456', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5086, 'n_tokens_neox': 4887, 'n_words': 717, 'path': 'notebooks/nicmos_hare_and_hounds.ipynb', 'repo_name': 'benjaminpope/poppy', 'size': '11989'}