def optimal_myopic(domain, domain_mu_prior, domain_cm_prior, xObs, yObs,\
                   xres, yres, ysdbounds, lenscale, sigvar, noisevar2):
    out = {}
    bounds = [domain[0], domain[-1]]
    xcandidates = linspace(bounds[0], bounds[1], xres)
    out['x'] = xcandidates
    out['ev'] = zeros_like(xcandidates)
    ysdcandidates = linspace(ysdbounds[0], ysdbounds[1], yres)
    pysdcandidates = norm.pdf(ysdcandidates)

    for ix0, x0 in enumerate(xcandidates):
        x0 = a([x0])
        xObsPlusX0 = append(xObs, x0)
        xcmpri0 = jbgp.K_se(x0, x0, lenscale, sigvar)
        ymu0 = jbgp.conditioned_mu(x0, xObs, yObs, lenscale, sigvar, noisevar2)
        ycm0 = jbgp.conditioned_covmat(x0, atleast_2d(xcmpri0), xObs, lenscale, sigvar, noisevar2)
        ysd0 = diag(ycm0)
        ycands = a([ymu0 + (ysd0 * n) for n in ysdcandidates])
        for iy0, y0 in enumerate(ycands):
            yObsPlusY0 = append(yObs, y0)
            py0 = pysdcandidates[iy0]
            mu0 = jbgp.conditioned_mu(domain, xobsPlusX0, yObsPlusY0, lenscale, sigvar, noisevar2)
            evmax = mu0.max()
            out['ev'][ix0] += evmax * py0

    return out

 def getNormal(self, u, v):
     if self.bump is None:
         return a([0,0,0])
     else:
         u, v = map(int,(u*(self.bump.size[0]-1),v*(self.bump.size[1]-1)))
         n = 2*a(self.bump.getpixel((u,v)), float)/255 - a([1.,1.,1])
         return Vec3(*n).normalized()

 def getColor(self, u, v):
     if self.tex is None:
         return a([1.,1.,1.])
     else:
         u, v = map(int,(u*(self.tex.size[0]-1),v*(self.tex.size[1]-1)))
         c = a(self.tex.getpixel((u,v)), float)
         return c/255

def draw_simplex_3d_euclidean_bounds(simplex, obj_nr=0, L=[1., 1.], points=[]):
    '''2d->2d simplex lower bound surface.
    Draws a single objective below 2D simplex'''
    from matplotlib import pyplot as plt
    from matplotlib import cm

    t = a([simplex[0][:-1], simplex[1][:-1], simplex[2][:-1]])
    y = a([simplex[0][-1]['obj'], simplex[1][-1]['obj'], simplex[2][-1]['obj']])   # (obj1, obj2) for A, B, C

    X1 =  np.arange(-0.1, 2.1, 0.05)
    X2 =  np.arange(-0.1, 2.1, 0.05)
    X1, X2 = np.meshgrid(X1, X2)
    fig = plt.figure()
    ax = fig.gca(projection='3d')

    LB = lower_bound_surface(X1, X2, t, y, L)

    ax.plot_surface(X1, X2, LB[:,:,obj_nr], linewidth=0, rstride=1, cstride=1, cmap=cm.coolwarm)
    ax.plot_wireframe(np.hstack((t[:,0], t[0,0])), np.hstack(([t[:,1], t[0,1]])), np.hstack((y[:,obj_nr],y[0,obj_nr])), zorder=1000)  ## Line surface

    # points = [simplex[-1]['mins_ABC'][1]]
    for p in points:
        ax.plot([p[0]], [p[1]], [p[2]], 'go')

    # points = [[1.5892857095626787, 1.5892857194077434, 1.186929245703222]]
    # for p in points:
    #     ax.plot([p[0]], [p[1]], [p[2]], 'ro')
    plt.show()

def emep(domainbounds, xObs, yObs, lenscale, sigvar, noisevar2, xres, yres, ysdbounds):
    """expected value of the max expected value of the posterior"""
    out = {}
    xcandidates = linspace(domainbounds[0], domainbounds[1], xres)
    out['x'] = xcandidates
    out['emep'] = zeros_like(xcandidates)
    ySDcandidates = linspace(ysdbounds[0], ysdbounds[1], yres)
    pdfySDcandidates = norm.pdf(ySDcandidates)
    cdfySDcandidates = norm.cdf(ySDcandidates)
    for ix0, xx0 in enumerate(xcandidates):
        # print 'x: ' + str(ix0)
        x0 = a([xx0])
        # get ymu and ysd for x0 so know what points to consider for generating prob-weighted ev of max of posterior
        ymu0 = jbgp.conditioned_mu(x0, xObs, yObs, lenscale, sigvar, noisevar2)
        xcmpri0 = jbgp.K_se(x0, x0, lenscale, sigvar)  # get covmat for xSam
        ycm0 = jbgp.conditioned_covmat(x0, atleast_2d(xcmpri0), xObs, lenscale, sigvar, noisevar2)
        ysd0 = diag(ycm0)
        # y-vals to consider with probs pysdcandidates
        ycands = a([ymu0 + (ysd0 * d) for d in ySDcandidates])
        xObsPlusX0 = append(xObs, x0)  # add considered point to xObs locations
        # run simulations of what happens with certain y-vals
        mep = zeros_like(cdfySDcandidates)

        for iy0, y0 in enumerate(ycands):
            # print 'y: ' + str(iy0)
            yObsPlusY0 = append(yObs, y0)
            py0 = pdfySDcandidates[iy0]
            mu0 = jbgp.conditioned_mu(domain, xObsPlusX0, yObsPlusY0, lenscale, sigvar, noisevar2)
            mep[iy0] = mu0.max()

        out['emep'][ix0] = trapz(maxevpost, cdfySDcandidates)
    return out

    def test_compute_expectation(self):
        '''
        Intermediate results expected

        Kernel matrix:
        1.0 .135335283
        .135335283 1.0

        K^-1:
        1.018657360 -0.137860282
        -0.137860282  1.018657360

        Kernel vector:
        .324652467
        .882496903

        k' * Kinv:
        0.209048353 0.854205285
        '''
        expected = predict_f(
            x=a([
                [-1.0, 1.0],
                [1.0, 1.0],
            ]),
            fmap=a([
                [1.0],
                [3.0],
            ]),
            xnew=a([0.5, 1.0]),
            kernelfunc=default_kernel
        )
        self.assertAlmostEqual(expected, 2.771664208)

def analyse_beam(data):
    beam = generate_beam_from_json(data)

    beam.calculations()

    sections = [(section.start, section.end) for section in beam.sections]
    distance = a([np.linspace(section.start, section.end) for section in beam.sections])

    return {
        'sections': sections,
        'distance': [distance[i].tolist() for i in range(beam.num_sections)],
        'shear_force': {
            'equations': beam.shear_force_eqs,
            'values': [beam.shear_force_values[i].tolist() for i in range(beam.shear_force_values.shape[0])],
            'plot': {
                'x': distance.reshape(-1).tolist(),
                'y': beam.shear_force_values.reshape(-1).tolist(),
                'backup': {
                    'x': a([np.linspace(s.start, s.end) for s in beam.sections]).reshape(-1).tolist(),
                    'y': a([np.linspace(beam.shear_force_values[i], beam.shear_force_values[i]) for i in
                            range(beam.num_sections)]).reshape(-1).tolist()
                }
            }
        },
        'bending_moment': {
            'equations': beam.bending_moment_eqs,
            'values': [beam.bending_moment_values[i].tolist() for i in range(beam.bending_moment_values.shape[0])],
            'plot': {
                'x': distance.reshape(-1).tolist(),
                'y': beam.bending_moment_values.reshape(-1).tolist()
            }
        }
    }

 def test_select_point_to_exploit(self):
     # We attempt to force exploitation by covering most of the input
     # space and expecting that the maximization algorithm will choose
     # the point between the highest outputs, given a symmetric output function.
     next_point = acquire(
         x=a([
             [-0.75],
             [-0.25],
             [0.25],
             [0.75],
         ]),
         # I got these fmap and Cmap values from running our optimizer
         # on the input data with comparisons [1, 0], [1, 3], [2, 0], [2, 3].
         fmap=a([
             [0.08950024],
             [0.21423927],
             [0.21423927],
             [0.08950024],
         ]),
         Cmap=a([
             [0.15672336, -0.07836168, -0.07836168, 0.0],
             [-0.07836168, 0.15672336, 0.0, -0.07836168],
             [-0.07836168, 0.0, 0.15672336, -0.07836168],
             [0.0, -0.07836168, -0.07836168, 0.15672336],
         ]),
         bounds=a([
             [-1.0, 1.0],
         ]),
         kernelfunc=default_kernel
     )
     self.assertTrue(next_point[0] > -.25)
     self.assertTrue(next_point[0] < .25)

def show_lower_pareto_bound(simplexes):
    from matplotlib import pyplot as plt
    '''For 2D -> 2D problem show the lower pareto bound.'''
    # Warning:  unfinished, untested.
    def lower_bound_for_interval(t, y, dist=None, L=[1.,1.], verts=[0,1]):
        '''Returns lower L bound line in (f1,f2) space for an interval in
        multidimensional feasible region.'''
        def lower_bound_cracks(x1, x2, y1, y2, dist):
            '''Computes lower L bound crack points for the given interval.'''
            if dist is None:
                dist = enorm(x2-x1)
            t1 = (y1[0]-y2[0])/(2.*L[0]) + dist/2.
            t2 = (y1[1]-y2[1])/(2.*L[1]) + dist/2.
            if t1 >= t2:
                p1 = [y1[0] - L[0]*t2, y1[1] - L[1]*t2]
                p2 = [y1[0] - L[0]*t1, y2[1] - dist*L[1] + L[1]*t1]
            else:
                p2 = [y2[0] - dist*L[0] + L[0]*t2, y1[1] - L[1]*t2]
                p1 = [y1[0] - L[0]*t1, y1[1] - L[1]*t1]
            return [p1, p2]
        p1, p2 = lower_bound_cracks(t[verts[0]], t[verts[1]], y[verts[0]], y[verts[1]], dist)
        return a([y[verts[0]], p1, p2, y[verts[1]]])

    # For each simplex we have dimensions+1 vertex, so what does this mean.
    # Patobulinimas: We could check if any of them are dominated and keep the
    # least dominated.
    for simplex in simplexes:
        t = a([simplex[0][:-1], simplex[1][:-1], simplex[2][:-1]])
        y = a([simplex[0][-1]['obj'], simplex[1][-1]['obj'], simplex[2][-1]['obj']])   # (obj1, obj2) for A, B, C
        lb = lower_bound_for_interval(t, y)
        plt.plot(lb[:,0], lb[:,1])

    # Draw longest (or nondominated) edge lower bound and mark these vertexes as stars.
    plt.show()

 def test_get_distinct_x_from_comparisons(self):
     comp = a([
         [[0.0], [1.0]],
         [[2.0], [3.0]]
     ])
     x = get_distinct_x(comp)
     self.assertAlmostEqual(x, a([[0.0], [1.0], [2.0], [3.0]]))

 def test_skip_repetitions_within_comparison(self):
     comp = a([
         [[0.0], [1.0]],
         [[2.0], [2.0]]
     ])
     x = get_distinct_x(comp)
     self.assertAlmostEqual(x, a([[0.0], [1.0], [2.0]]))

 def test_run_optimization(self):
     f, _ = newton_rhapson(
         x=a([
             [0.0, 1.0],
             [1.0, 1.0],
             [-1.0, 0.0],
             [2.0, 2.0],
             [1.5, -1.0]
         ]),
         f0=a([
             [0.0],
             [0.0],
             [0.0],
             [0.0],
             [0.0],
         ]),
         comparisons=a([
             [3, 1],
             [0, 1],
             [2, 1],
             [4, 0],
             [2, 4],
         ]),
         kernelfunc=default_kernel,
         Hfunc=compute_H,
         gfunc=compute_g,
         sigma=2,
         maxiter=20,
     )
     self.assertTrue(f[3][0] > f[1][0])
     self.assertTrue(f[0][0] > f[1][0])
     self.assertTrue(f[2][0] > f[1][0])
     self.assertTrue(f[4][0] > f[0][0])
     self.assertTrue(f[2][0] > f[4][0])

def draw_bounds(x1, x2, L=[1.,1.]):
    '''Draws 3 plots describing bounds for 1D->2D problems'''
    from matplotlib import pyplot as plt
    import matplotlib.ticker as plticker
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18,6))

    ## Normalize x1, x2 for 2D variable space
    x2 = [enorm(a(x2[:-1]) - a(x1[:-1])), x2[-1]]
    x1 = [0, x1[-1]]

    ax1 = draw_objective_bounds_and_hat_epsilon(x1, x2, ax=ax1, L=L)
    ax2 = draw_tolerance_change(x1, x2, ax=ax2, L=L)
    ax3 = draw_objective_bounds(x1, x2, ax=ax3, L=L)

    ax1.xaxis.set_major_locator(plticker.MultipleLocator(base=0.2))
    ax1.yaxis.set_major_locator(plticker.MultipleLocator(base=0.1))
    # ax1.axis([0,1.,0,1.2])

    ax2.xaxis.set_major_locator(plticker.MultipleLocator(base=0.2))
    ax2.yaxis.set_major_locator(plticker.MultipleLocator(base=0.1))

    ax3.xaxis.set_major_locator(plticker.MultipleLocator(base=0.2))
    ax3.yaxis.set_major_locator(plticker.MultipleLocator(base=0.1))
    # ax3.axis([0,1.2,0,1.2])

    plt.show()

 def test_skip_repetitions_across_comparisons(self):
     comp = a([
         [[1.0], [2.0]],
         [[2.0], [3.0]]
     ])
     x = get_distinct_x(comp)
     self.assertAlmostEqual(x, a([[1.0], [2.0], [3.0]]))

 def test_get_distinct_x_from_2_dimensional_input_data(self):
     comp = a([
         [[1.0, 2.0], [2.0, 2.0]],
         [[2.0, 2.0], [3.0, 3.0]]
     ])
     x = get_distinct_x(comp)
     self.assertAlmostEqual(x, a([[1.0, 2.0], [2.0, 2.0], [3.0, 3.0]]))

 def test_select_point_to_explore_in_two_dimensions(self):
     '''
     This test case just runs code on a multi-dimensional input dataset.
     It doesn't test any conditions, but could presumably do so in the future.
     '''
     acquire(
         x=a([
             [-0.5, -0.5],
             [-0.5, 0.5],
         ]),
         # Fake fmap and Cmap values from another set of x's
         fmap=a([
             [0.03254087],
             [-0.03254087],
         ]),
         Cmap=a([
             [0.07894662, -0.07894662],
             [-0.07894662, 0.07894662],
         ]),
         bounds=a([
             [-1.0, 1.0],
             [-1.0, 1.0],
         ]),
         kernelfunc=default_kernel
     )

 def setUp(self):
     self.default_f = a([
         [6.0],
         [1.0],
         [2.0],
     ])
     self.default_comparison = a([0, 2], dtype=np.int)
     self.default_sigma = 2.0

 def test_compute_kernel_vector(self):
     x = a([
         [0.0, 1.0],
         [1.0, 1.0]
     ])
     xnew = a([0.0, 0.0])
     k = kernel_vector(default_kernel, x, xnew)
     self.assertAlmostEqual(k, a([[.60653066], [.367879441]]))

 def testNorm1(self):
     '''test Norm1.'''
     example = ( ( a( [1,1,1] ), 1 ) ,
                 ( a( [[1,3,4],[2,-2,1]] ), a( [4,2] )) 
                )
     for x,y in example:
         rslt = Band.norm1(x)
         npt.assert_allclose(rslt, y)

 def testBlockSub2CenterCarWithoutBand(self):
     '''test BlockSub2CenterCarWithoutBand'''
     example = ( ( 0, a([-1.9,-1.9,-1.9])), 
                 ( 1, a([-1.7,-1.9,-1.9])) )
     for ind,right in example:
         sub = self.bnd.BlockInd2SubWithoutBand(ind)
         rslt = self.bnd.BlockSub2CenterCarWithoutBand(sub)
         npt.assert_array_equal(rslt, right)