numpy.r_

Here are the examples of the python api numpy.r_ taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

200 Examples 7

Example 51

Project: skl-groups
Source File: features.py
View license
    def __add__(self, oth):
        if isinstance(oth, Features):
            meta = {k: np.r_[self.meta[k], oth.meta[k]]
                    for k in self.meta if k in oth.meta}
            oth_features = oth.features
        elif isinstance(oth, list):
            meta = {}
            oth_features = np.empty(len(oth), object)
            oth_features[:] = oth
        else:
            return NotImplemented

        return Features(np.r_[self.features, oth_features],
                        stack=False, copy=True, **meta)

Example 52

Project: rayopt
Source File: utils.py
View license
@public
def gl_roots(n):
    """Gauss Lobatto roots and weights for [-1, 1]
    with -1 first and 1 last
    """
    leg = orthogonal.legendre(n - 1)
    x = np.r_[-1, leg.deriv().roots, 1]
    w = 2/(n*(n - 1)*leg(x)**2)
    return x, w

Example 53

Project: rayopt
Source File: utils.py
View license
@public
def gl_roots(n):
    """Gauss Lobatto roots and weights for [-1, 1]
    with -1 first and 1 last
    """
    leg = orthogonal.legendre(n - 1)
    x = np.r_[-1, leg.deriv().roots, 1]
    w = 2/(n*(n - 1)*leg(x)**2)
    return x, w

Example 54

Project: rapprentice
Source File: retiming.py
View license
def retime2(positions, vel_limits_j):
    good_inds = np.r_[True,(abs(positions[1:] - positions[:-1]) >= 1e-6).any(axis=1)]
    positions = positions[good_inds]
    
    move_nj = positions[1:] - positions[:-1]
    dur_n = (np.abs(move_nj) / vel_limits_j[None,:]).max(axis=1) # time according to velocity limit
    times = np.cumsum(np.r_[0,dur_n])
    return times, good_inds    

Example 55

Project: rapprentice
Source File: retiming.py
View license
def make_traj_with_limits(positions, vel_limits_j, acc_limits_j):
    times, inds = retime(positions, vel_limits_j, acc_limits_j)
    positions = positions[inds]

    deg = min(3, len(positions) - 1)
    s = len(positions)*.001**2
    (tck, _) = si.splprep(positions.T, s=s, u=times, k=deg)
    smooth_positions = np.r_[si.splev(times,tck,der=0)].T
    smooth_velocities = np.r_[si.splev(times,tck,der=1)].T    
    return smooth_positions, smooth_velocities, times    

Example 56

Project: Haystack
Source File: haystack_motifs_CORE.py
View license
def smooth(x,window_len=200,window='hanning'):
    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"


    s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
    if window == 'flat': #moving average
        w=np.ones(window_len,'d')
    else:
        w=eval('np.'+window+'(window_len)')
    y=np.convolve(w/w.sum(),s,mode='valid')
    return y[window_len/2:-window_len/2+1]

Example 57

Project: Haystack
Source File: haystack_motifs_CORE.py
View license
def smooth(x,window_len=200,window='hanning'):
    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"


    s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
    if window == 'flat': #moving average
        w=np.ones(window_len,'d')
    else:
        w=eval('np.'+window+'(window_len)')
    y=np.convolve(w/w.sum(),s,mode='valid')
    return y[window_len/2:-window_len/2+1]

Example 58

Project: mne-python
Source File: search_light.py
View license
def _sl_init_pred(y_pred, X):
    """Aux. function to _SearchLight to initialize y_pred."""
    n_sample, n_iter = X.shape[0], X.shape[-1]
    if y_pred.ndim > 1:
        # for estimator that generate multidimensional y_pred,
        # e.g. clf.predict_proba()
        y_pred = np.zeros(np.r_[n_sample, n_iter, y_pred.shape[1:]],
                          y_pred.dtype)
    else:
        # for estimator that generate unidimensional y_pred,
        # e.g. clf.predict()
        y_pred = np.zeros((n_sample, n_iter), y_pred.dtype)
    return y_pred

Example 59

Project: mne-python
Source File: search_light.py
View license
def _sl_init_pred(y_pred, X):
    """Aux. function to _SearchLight to initialize y_pred."""
    n_sample, n_iter = X.shape[0], X.shape[-1]
    if y_pred.ndim > 1:
        # for estimator that generate multidimensional y_pred,
        # e.g. clf.predict_proba()
        y_pred = np.zeros(np.r_[n_sample, n_iter, y_pred.shape[1:]],
                          y_pred.dtype)
    else:
        # for estimator that generate unidimensional y_pred,
        # e.g. clf.predict()
        y_pred = np.zeros((n_sample, n_iter), y_pred.dtype)
    return y_pred

Example 60

Project: mne-python
Source File: _stockwell.py
View license
def _precompute_st_windows(n_samp, start_f, stop_f, sfreq, width):
    """Precompute stockwell gausian windows (in the freq domain)."""
    tw = fftpack.fftfreq(n_samp, 1. / sfreq) / n_samp
    tw = np.r_[tw[:1], tw[1:][::-1]]

    k = width  # 1 for classical stowckwell transform
    f_range = np.arange(start_f, stop_f, 1)
    windows = np.empty((len(f_range), len(tw)), dtype=np.complex)
    for i_f, f in enumerate(f_range):
        if f == 0.:
            window = np.ones(len(tw))
        else:
            window = ((f / (np.sqrt(2. * np.pi) * k)) *
                      np.exp(-0.5 * (1. / k ** 2.) * (f ** 2.) * tw ** 2.))
        window /= window.sum()  # normalisation
        windows[i_f] = fftpack.fft(window)
    return windows

Example 61

Project: mne-python
Source File: _stockwell.py
View license
def _precompute_st_windows(n_samp, start_f, stop_f, sfreq, width):
    """Precompute stockwell gausian windows (in the freq domain)."""
    tw = fftpack.fftfreq(n_samp, 1. / sfreq) / n_samp
    tw = np.r_[tw[:1], tw[1:][::-1]]

    k = width  # 1 for classical stowckwell transform
    f_range = np.arange(start_f, stop_f, 1)
    windows = np.empty((len(f_range), len(tw)), dtype=np.complex)
    for i_f, f in enumerate(f_range):
        if f == 0.:
            window = np.ones(len(tw))
        else:
            window = ((f / (np.sqrt(2. * np.pi) * k)) *
                      np.exp(-0.5 * (1. / k ** 2.) * (f ** 2.) * tw ** 2.))
        window /= window.sum()  # normalisation
        windows[i_f] = fftpack.fft(window)
    return windows

Example 62

Project: apsg
Source File: plotting.py
View license
    def format_coord(self, x, y):
        a, b = self.Ti.dot(np.r_[x, y]-self.C)
        c = 1 - a - b
        if a<0 or b<0 or c<0:
            return ''
        else:
            return 'P:{:0.2f} G:{:0.2f} R:{:0.2f}'.format(a,b,c)

Example 63

Project: deep-learning-from-scratch
Source File: util.py
View license
def smooth_curve(x):
    """損失関数のグラフを滑らかにするために用いる

    参考:http://glowingpython.blogspot.jp/2012/02/convolution-with-numpy.html
    """
    window_len = 11
    s = np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
    w = np.kaiser(window_len, 2)
    y = np.convolve(w/w.sum(), s, mode='valid')
    return y[5:len(y)-5]

Example 64

Project: numdifftools
Source File: test_fornberg.py
View license
    @staticmethod
    def test_weights():
        x = np.r_[-1, 0, 1]
        xbar = 0
        k = 1
        weights = fornberg_weights(x, xbar, k)
        np.testing.assert_allclose(weights, [-.5, 0, .5])

Example 65

Project: numdifftools
Source File: test_numdifftools.py
View license
    @staticmethod
    def test_directional_diff():
        v = np.r_[1, -1]
        v = v / np.linalg.norm(v)
        x0 = [2, 3]
        directional_diff = np.dot(nd.Gradient(rosen)(x0), v)
        assert_allclose(directional_diff, 743.87633380824832)
        dd, _info = nd.directionaldiff(rosen, x0, v, full_output=True)
        assert_allclose(dd, 743.87633380824832)

Example 66

Project: popupcad
Source File: spline_functions.py
View license
def control_from_fit_points(fp,p,plot = True):
    import scipy.linalg
    
    r = len(fp)-1
    middle = numpy.r_[1:r]
    knots = numpy.r_[[0]*(p+1),middle,[r]*(p+1)]
    m = len(knots)-1
    n = m-p-1
    w = [1]*(n+1)
    
    domain = numpy.r_[0:r+1]
    y = calc_spline(knots,n,w,domain)
    yi = scipy.linalg.pinv(y)
    
    fp = numpy.array(fp)
    
    cp2 = yi.dot(fp)
    if plot:
        plt.plot(domain,y)
    
    return cp2,knots

Example 67

Project: popupcad
Source File: vertex.py
View license
    def getpos3D(self):
        return tuple(numpy.r_[self.getpos(),0].tolist())

Example 68

Project: attention-lvcsr
Source File: speed_test_conv.py
View license
def exec_multilayer_conv_nnet_old(conv_mode, ss, bsize, imshp, kshps, nkerns,
        unroll_batch=0, unroll_kern=0, img=T.dmatrix(), validate=True,
        conv_op_py=False, do_print=True, repeat=1,
        unroll_patch=False, unroll_patch_size=False, verbose=0):

        # build actual input images
        imgval = global_rng.rand(bsize, imshp[0], imshp[1], imshp[2])

        a = T.dmatrix()
        kerns = [a for i in nkerns]
        inputs4 = dmatrix4()
        kerns4 = dmatrix4()

        # for each layer
        ntot = 0
        tctot = 0
        tpytot = 0

        for kshp, kern, nkern, n_layer in zip(kshps, kerns, nkerns,
                                              xrange(len(nkerns))):
            if do_print:
                print('************* layer %i ***************' % n_layer)
                
                print(conv_mode, ss, n_layer, kshp, nkern)

            # actual values
            w = global_rng.random_sample(N.r_[nkern, imshp[0], kshp])
            w_flip = flip(w, kshp).reshape(w.shape)

            # manual implementation
            # check first stage
            padimg = imgval
            if conv_mode == 'full':
                padimg_shp = N.array(imshp[1:]) + 2*(N.array(kshp) - N.array([1, 1]))
                padimg = N.zeros(N.r_[bsize, imshp[0], padimg_shp])
                padimg[:, :, kshp[0]-1:-kshp[0]+1,
                             kshp[1]-1:-kshp[1]+1] = imgval

            outshp = N.hstack((nkern, ConvOp.getOutputShape(imshp[1:], kshp, ss, conv_mode)))

            time1 = time.time()
            outval = N.zeros(N.r_[bsize, outshp])
            if validate:
                # causes an atexit problem
                from scipy.signal.sigtools import _convolve2d
                from scipy.signal.signaltools import  _valfrommode, _bvalfromboundary
                val = _valfrommode(conv_mode)
                bval = _bvalfromboundary('fill')
                for b in xrange(bsize):  # loop over batches
                    for n in xrange(nkern):  # loop over filters
                        for i in xrange(imshp[0]):  # loop over input feature maps
                            outval[b, n, ...] +=  _convolve2d(\
                                imgval[b, i, ...], w_flip[n, i, ...], 1, val, bval, 0)[0::ss[0], 0::ss[1]]
                ntot += time.time() - time1

            # ConvOp
            if unroll_patch and not unroll_patch_size:
                conv_op = ConvOp(dx=ss[0], dy=ss[1], output_mode=conv_mode,
                                 unroll_patch=unroll_patch, verbose=verbose)(inputs4, kerns4)
            else:
                conv_op = ConvOp(imshp, kshp, nkern, bsize, ss[0], ss[1], conv_mode,
                                 unroll_batch=unroll_batch, unroll_kern=unroll_kern, unroll_patch=unroll_patch, verbose=verbose)(inputs4, kerns4)
            l1shp = N.hstack((nkern,
                            ConvOp.getOutputShape(imshp[1:], kshp, ss, conv_mode)))
            propup2 = function([inputs4, kerns4], conv_op)
            propup3 = function([inputs4, kerns4], conv_op, mode=Mode(linker="py"))

            time1 = time.time()
            for i in xrange(repeat):
                hidval2_ = propup2(imgval, w_flip)
            hidval2 = hidval2_  # [:,:,0::ss[0],0::ss[1]]
            tctot += time.time() - time1

            if conv_op_py:
                time1 = time.time()
                for i in xrange(repeat):
                    hidval3_ = propup3(imgval, w_flip)
                hidval3 = hidval3_  # [:,:,0::ss[0],0::ss[1]]
                tpytot += time.time() - time1
                assert (N.abs(hidval2-hidval3) < 1e-5).all()
            else:
                tpytot += 0

            if validate:
                temp = N.abs(outval - hidval2)
                assert (temp < 1e-5).all()
            if validate and conv_op_py:
                temp = N.abs(outval - hidval3)
                assert (temp < 1e-5).all()

            imshp = tuple(outshp)
            imgval = outval.reshape(bsize, outshp[0], outshp[1], outshp[2])

        return tctot, tpytot, ntot

Example 69

Project: attention-lvcsr
Source File: speed_test_conv.py
View license
def exec_multilayer_conv_nnet(conv_mode, ss, bsize, imshp, kshps, nkerns,
        unroll_batch=0, unroll_kern=0, img=T.dmatrix(),
        do_print=True, repeat=1,
        unroll_patch=False, unroll_patch_size=False, verbose=0):

        # build actual input images
        imgval = global_rng.rand(bsize, imshp[0], imshp[1], imshp[2])

        a = T.dmatrix()
        kerns = [a for i in nkerns]
        inputs4 = dmatrix4()
        kerns4 = dmatrix4()

        # for each layer
        ntot = 0
        tctot = 0
        tpytot = 0

        for kshp, kern, nkern, n_layer in zip(kshps, kerns, nkerns, xrange(len(nkerns))):
            if do_print:
                print('************* layer %i ***************' % n_layer)
                
                print(conv_mode, ss, n_layer, kshp, nkern)

            # actual values
            w = global_rng.random_sample(N.r_[nkern, imshp[0], kshp])
            w_flip = flip(w, kshp).reshape(w.shape)

            outshp = N.hstack((nkern, ConvOp.getOutputShape(imshp[1:], kshp, ss, conv_mode)))

            time1 = time.time()
            outval = N.zeros(N.r_[bsize, outshp])

            # ConvOp
            if unroll_patch and not unroll_patch_size:
                conv_op = ConvOp(dx=ss[0], dy=ss[1], output_mode=conv_mode,
                                 unroll_patch=unroll_patch, verbose=verbose)(inputs4, kerns4)
            else:
                conv_op = ConvOp(imshp, kshp, nkern, bsize, ss[0], ss[1], conv_mode,
                                 unroll_batch=unroll_batch, unroll_kern=unroll_kern, unroll_patch=unroll_patch, verbose=verbose)(inputs4, kerns4)
            l1shp = N.hstack((nkern,
                            ConvOp.getOutputShape(imshp[1:], kshp, ss, conv_mode)))
            propup2 = function([inputs4, kerns4], conv_op)

            time1 = time.time()
            for i in xrange(repeat):
                hidval2_ = propup2(imgval, w_flip)
            hidval2 = hidval2_  # [:,:,0::ss[0],0::ss[1]]
            tctot += time.time() - time1

            imshp = tuple(outshp)
            imgval = outval.reshape(bsize, outshp[0], outshp[1], outshp[2])

        return tctot, tpytot, ntot

Example 70

Project: lfd
Source File: move_rope.py
View license
def create_augmented_traj(robot, pick_pos, drop_pos, pick_R, drop_R, move_height, pos_displacement_per_step=.02):
    pos_traj = np.array([pick_pos + np.r_[0,0,move_height], 
                         pick_pos, 
                         pick_pos + np.r_[0,0,move_height], 
                         drop_pos + np.r_[0,0,move_height], 
                         drop_pos, 
                         drop_pos + np.r_[0,0,move_height]])
    R_traj = np.array([pick_R, pick_R, pick_R, drop_R, drop_R, drop_R])
    ee_traj = np.empty((len(pos_traj), 4, 4))
    ee_traj[:] = np.eye(4)
    ee_traj[:,:3,3] = pos_traj
    ee_traj[:,:3,:3] = R_traj
    open_finger_traj = np.array([False, False, False, False, True, False])
    close_finger_traj = np.array([False, True, False, False, False, False])
    open_finger_value = sim_util.get_binary_gripper_angle(True)
    closed_finger_value = sim_util.get_binary_gripper_angle(False)
    finger_traj = np.array([open_finger_value, open_finger_value, closed_finger_value, closed_finger_value, closed_finger_value, open_finger_value])[:,None]
    
    lr = 'r' # use right arm/gripper
    aug_traj = AugmentedTrajectory(lr2ee_traj={lr: ee_traj}, lr2finger_traj={lr: finger_traj}, lr2open_finger_traj={lr: open_finger_traj}, lr2close_finger_traj={lr: close_finger_traj})
    
    # resample augmented trajectory according to the position displacement
    pos_dists = np.linalg.norm(np.diff(pos_traj, axis=0), axis=1)
    new_times = np.r_[0, np.cumsum(pos_dists/pos_displacement_per_step)]
    timesteps_rs = np.interp(np.arange(new_times[-1]+1), new_times, np.arange(len(new_times)))
    aug_traj_rs = aug_traj.get_resampled_traj(timesteps_rs)
    
    # do motion planning for aug_traj_rs
    manip_name = {"l":"leftarm", "r":"rightarm"}[lr]
    ee_link_name = "%s_gripper_tool_frame"%lr
    ee_link = robot.GetLink(ee_link_name)
    dof_vals = robot.GetManipulator(manip_name).GetArmDOFValues()
    init_traj = np.tile(dof_vals, (aug_traj_rs.n_steps,1))
    arm_traj, _, _ = planning.plan_follow_traj(robot, manip_name, ee_link, aug_traj_rs.lr2ee_traj['r'], init_traj, no_collision_cost_first=True)
    aug_traj_rs.lr2arm_traj[lr] = arm_traj
    
    return aug_traj_rs

Example 71

Project: lfd
Source File: move_rope.py
View license
def main():
    # define simulation objects
    table_height = 0.77
    cyl_radius = 0.025
    cyl_height = 0.3
    cyl_pos0 = np.r_[.7, -.15, table_height+cyl_height/2]
    cyl_pos1 = np.r_[.7,  .15, table_height+cyl_height/2]
    cyl_pos2 = np.r_[.4, -.15, table_height+cyl_height/2]
    rope_poss = np.array([[.2, -.2, table_height+0.006], 
                          [.8, -.2, table_height+0.006], 
                          [.8,  .2, table_height+0.006], 
                          [.2,  .2, table_height+0.006]])
    
    sim_objs = []
    sim_objs.append(XmlSimulationObject("robots/pr2-beta-static.zae", dynamic=False))
    sim_objs.append(BoxSimulationObject("table", [1, 0, table_height-.1], [.85, .85, .1], dynamic=False))
    cyl_sim_objs = create_cylinder_grid(cyl_pos0, cyl_pos1, cyl_pos2, cyl_radius, cyl_height)
    sim_objs.extend(cyl_sim_objs)
    rope_sim_obj = create_rope(rope_poss)
    sim_objs.append(rope_sim_obj)
    
    # initialize simulation world and environment
    sim = DynamicSimulationRobotWorld()
    sim.add_objects(sim_objs)
    sim.create_viewer()
    
    sim.robot.SetDOFValues([0.25], [sim.robot.GetJoint('torso_lift_joint').GetJointIndex()])
    sim_util.reset_arms_to_side(sim)
    
    color_cylinders(cyl_sim_objs)
    
    env = environment.LfdEnvironment(sim, sim)
    
    # define augmented trajectory
    pick_pos = rope_poss[0] + .1 * (rope_poss[1] - rope_poss[0])
    drop_pos = rope_poss[3] + .1 * (rope_poss[2] - rope_poss[3]) + np.r_[0, .2, 0]
    pick_R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
    drop_R = np.array([[0, 1, 0], [0, 0, -1], [-1, 0, 0]])
    move_height = .2
    aug_traj = create_augmented_traj(sim.robot, pick_pos, drop_pos, pick_R, drop_R, move_height)
    
    env.execute_augmented_trajectory(aug_traj)

Example 72

Project: lfd
Source File: simulation.py
View license
    def execute_trajectory(self, full_traj, step_viewer=1, interactive=False,
                           max_cart_vel_trans_traj=.05, sim_callback=None):
        # TODO: incorporate other parts of sim_full_traj_maybesim
        if sim_callback is None:
            sim_callback = lambda i: self.step()
        
        traj, dof_inds = full_traj
        
    #     # clip finger joint angles to the binary gripper angles if necessary
    #     for lr in 'lr':
    #         joint_ind = self.robot.GetJoint("%s_gripper_l_finger_joint"%lr).GetDOFIndex()
    #         if joint_ind in dof_inds:
    #             ind = dof_inds.index(joint_ind)
    #             traj[:,ind] = np.minimum(traj[:,ind], get_binary_gripper_angle(True))
    #             traj[:,ind] = np.maximum(traj[:,ind], get_binary_gripper_angle(False))
        
        # in simulation mode, we must make sure to gradually move to the new starting position
        self.robot.SetActiveDOFs(dof_inds)
        curr_vals = self.robot.GetActiveDOFValues()
        transition_traj = np.r_[[curr_vals], [traj[0]]]
        sim_util.unwrap_in_place(transition_traj, dof_inds=dof_inds)
        transition_traj = ropesim.retime_traj(self.robot, dof_inds, transition_traj,
                                              max_cart_vel=max_cart_vel_trans_traj)
        animate_traj.animate_traj(transition_traj, self.robot, restore=False, pause=interactive,
                                  callback=sim_callback, step_viewer=step_viewer if self.viewer else 0)
        
        traj[0] = transition_traj[-1]
        sim_util.unwrap_in_place(traj, dof_inds=dof_inds)
        traj = ropesim.retime_traj(self.robot, dof_inds, traj)  # make the trajectory slow enough for the simulation
    
        animate_traj.animate_traj(traj, self.robot, restore=False, pause=interactive,
                                  callback=sim_callback, step_viewer=step_viewer if self.viewer else 0)
        if self.viewer and step_viewer:
            self.viewer.Step()
        return True

Example 73

Project: lfd
Source File: sim_util.py
View license
def mirror_arm_joints(x):
    "mirror image of joints (r->l or l->r)"
    return np.r_[-x[0],x[1],-x[2],x[3],-x[4],x[5],-x[6]]

Example 74

Project: lfd
Source File: math_utils.py
View license
def remove_duplicate_rows(mat):
    diffs = mat[1:] - mat[:-1]
    return mat[np.r_[True,(abs(diffs) >= 1e-5).any(axis=1)]]

Example 75

Project: lfd
Source File: plotting_plt.py
View license
def plot_tps_registration_3d(x_nd, y_md, f, res, x_color, y_color, xwarped_color):
    # set interactive
    plt.ion()
     
    fig = plt.figure('3d plot')
    fig.clear()

    ax = fig.add_subplot(121, projection='3d')
    ax.set_aspect('equal')

    ax.scatter(x_nd[:,0], x_nd[:,1], x_nd[:,2], c=x_color, edgecolors=x_color, marker=',', s=5)

    # manually set axes limits at a cube's bounding box since matplotlib doesn't correctly set equal axis in 3D
    xwarped_nd = f.transform_points(x_nd)
    max_pts = np.r_[x_nd, y_md, xwarped_nd].max(axis=0)
    min_pts = np.r_[x_nd, y_md, xwarped_nd].min(axis=0)
    max_range = (max_pts - min_pts).max()
    center = 0.5*(max_pts + min_pts)
    ax.set_xlim(center[0] - 0.5*max_range, center[0] + 0.5*max_range)
    ax.set_ylim(center[1] - 0.5*max_range, center[1] + 0.5*max_range)
    ax.set_zlim(center[2] - 0.5*max_range, center[2] + 0.5*max_range)

    grid_means = .5 * (x_nd.max(axis=0) + x_nd.min(axis=0))
    grid_mins = grid_means - (x_nd.max(axis=0) - x_nd.min(axis=0))
    grid_maxs = grid_means + (x_nd.max(axis=0) - x_nd.min(axis=0))
    plot_warped_grid_3d(lambda xyz: xyz, grid_mins, grid_maxs, xres=res[0], yres=res[1], zres=res[2], draw=False)

    ax = fig.add_subplot(122, projection='3d')
    ax.set_aspect('equal')

    ax.scatter(y_md[:,0], y_md[:,1], y_md[:,2], c=y_color, marker='+', s=50)
    xwarped_nd = f.transform_points(x_nd)
    ax.scatter(xwarped_nd[:,0], xwarped_nd[:,1], xwarped_nd[:,2], edgecolors=xwarped_color, facecolors='none', marker='o', s=50)

    ax.set_xlim(center[0] - 0.5*max_range, center[0] + 0.5*max_range)
    ax.set_ylim(center[1] - 0.5*max_range, center[1] + 0.5*max_range)
    ax.set_zlim(center[2] - 0.5*max_range, center[2] + 0.5*max_range)

    plot_warped_grid_3d(f.transform_points, grid_mins, grid_maxs, xres=res[0], yres=res[1], zres=res[2], draw=False)
    
    plt.draw()

Example 76

Project: lfd
Source File: retiming.py
View license
def remove_duplicate_rows(mat):
    diffs = mat[1:] - mat[:-1]
    return mat[np.r_[True,(abs(diffs) >= 1e-6).any(axis=1)]]

Example 77

Project: lfd
Source File: rope_initialization.py
View license
def remove_duplicate_rows(mat):
    diffs = mat[1:] - mat[:-1]
    return mat[np.r_[True,(abs(diffs) >= 1e-5).any(axis=1)]]

Example 78

Project: lfd
Source File: tps_registration.py
View license
def main():
    import argparse, h5py, os
    import matplotlib.pyplot as plt
    from rapprentice import clouds, plotting_plt
    import registration
    import time
    
    parser = argparse.ArgumentParser()
    parser.add_argument("input_file", type=str)
    parser.add_argument("--output_folder", type=str, default="")
    parser.add_argument("--i_start", type=int, default=0)
    parser.add_argument("--i_end", type=int, default=-1)
    regtype_choices = ['rpm', 'rpm-cheap', 'rpm-bij', 'rpm-bij-cheap']
    parser.add_argument("--regtypes", type=str, nargs='*', choices=regtype_choices, default=regtype_choices)
    parser.add_argument("--plot_color", type=int, default=1)
    parser.add_argument("--proj", type=int, default=1, help="project 3d visualization into 2d")
    parser.add_argument("--visual_prior", type=int, default=1)
    parser.add_argument("--plotting", type=int, default=1)

    args = parser.parse_args()
    
    def plot_cb_gen(output_prefix, args, x_color, y_color):
        def plot_cb(x_nd, y_md, corr_nm, f, iteration):
            if args.plot_color:
                plotting_plt.plot_tps_registration(x_nd, y_md, f, x_color = x_color, y_color = y_color, proj_2d=args.proj)
            else:
                plotting_plt.plot_tps_registration(x_nd, y_md, f, proj_2d=args.proj)
            # save plot to file
            if output_prefix is not None:
                plt.savefig(output_prefix + "_iter" + str(iteration) + '.png')
        return plot_cb

    def plot_cb_bij_gen(output_prefix, args, x_color, y_color):
        def plot_cb_bij(x_nd, y_md, xtarg_nd, corr_nm, wt_n, f):
            if args.plot_color:
                plotting_plt.plot_tps_registration(x_nd, y_md, f, res = (.3, .3, .12), x_color = x_color, y_color = y_color, proj_2d=args.proj)
            else:
                plotting_plt.plot_tps_registration(x_nd, y_md, f, res = (.4, .3, .12), proj_2d=args.proj)
            # save plot to file
            if output_prefix is not None:
                plt.savefig(output_prefix + "_iter" + str(iteration) + '.png')
        return plot_cb_bij

    # preprocess and downsample clouds
    DS_SIZE = 0.025
    infile = h5py.File(args.input_file)
    source_clouds = {}
    target_clouds = {}
    for i in range(args.i_start, len(infile) if args.i_end==-1 else args.i_end):
        source_cloud = clouds.downsample(infile[str(i)]['source_cloud'][()], DS_SIZE)
        source_clouds[i] = source_cloud
        target_clouds[i] = []
        for (cloud_key, target_cloud) in infile[str(i)]['target_clouds'].iteritems():
            target_cloud = clouds.downsample(target_cloud[()], DS_SIZE)
            target_clouds[i].append(target_cloud)
    infile.close()
    
    tps_costs = []
    tps_reg_costs = []
    for regtype in args.regtypes:
        start_time = time.time()
        costs = []
        reg_costs = []
        for i in range(args.i_start, len(source_clouds) if args.i_end==-1 else args.i_end):
            source_cloud = source_clouds[i]
            for target_cloud in target_clouds[i]:
                if args.visual_prior:
                    vis_cost_xy = ab_cost(source_cloud, target_cloud)
                else:
                    vis_cost_xy = None
                if regtype == 'rpm':
                    f, corr_nm = tps_rpm(source_cloud[:,:-3], target_cloud[:,:-3],
                                         vis_cost_xy = vis_cost_xy,
                                         plotting=args.plotting, plot_cb = plot_cb_gen(os.path.join(args.output_folder, str(i) + "_" + cloud_key + "_rpm") if args.output_folder else None,
                                                                                       args,
                                                                                       source_cloud[:,-3:],
                                                                                       target_cloud[:,-3:]))
                elif regtype == 'rpm-cheap':
                    f, corr_nm = tps_rpm(source_cloud[:,:-3], target_cloud[:,:-3],
                                         vis_cost_xy = vis_cost_xy, n_iter = N_ITER_CHEAP, em_iter = EM_ITER_CHEAP, 
                                         plotting=args.plotting, plot_cb = plot_cb_gen(os.path.join(args.output_folder, str(i) + "_" + cloud_key + "_rpm_cheap") if args.output_folder else None,
                                                                                       args,
                                                                                       source_cloud[:,-3:],
                                                                                       target_cloud[:,-3:]))
                elif regtype == 'rpm-bij':
                    x_nd = source_cloud[:,:3]
                    y_md = target_cloud[:,:3]
                    scaled_x_nd, _ = registration.unit_boxify(x_nd)
                    scaled_y_md, _ = registration.unit_boxify(y_md)
                    f,g = registration.tps_rpm_bij(scaled_x_nd, scaled_y_md, rot_reg=np.r_[1e-4, 1e-4, 1e-1], n_iter=50, reg_init=10, reg_final=.1, outlierfrac=1e-2, vis_cost_xy=vis_cost_xy,
                                                   plotting=args.plotting, plot_cb=plot_cb_bij_gen(os.path.join(args.output_folder, str(i) + "_" + cloud_key + "_rpm_bij") if args.output_folder else None,
                                                                                                   args,
                                                                                                   source_cloud[:,-3:],
                                                                                                   target_cloud[:,-3:]))
                elif regtype == 'rpm-bij-cheap':
                    x_nd = source_cloud[:,:3]
                    y_md = target_cloud[:,:3]
                    scaled_x_nd, _ = registration.unit_boxify(x_nd)
                    scaled_y_md, _ = registration.unit_boxify(y_md)
                    f,g = registration.tps_rpm_bij(scaled_x_nd, scaled_y_md, rot_reg=np.r_[1e-4, 1e-4, 1e-1], n_iter=10, outlierfrac=1e-2, vis_cost_xy=vis_cost_xy, # Note registration_cost_cheap in rope_qlearn has a different rot_reg and outlierfrac
                                                   plotting=args.plotting, plot_cb=plot_cb_bij_gen(os.path.join(args.output_folder, str(i) + "_" + cloud_key + "_rpm_bij_cheap") if args.output_folder else None,
                                                                                                   args,
                                                                                                   source_cloud[:,-3:],
                                                                                                   target_cloud[:,-3:]))
                costs.append(f._cost)
                reg_costs.append(registration.tps_reg_cost(f))
        tps_costs.append(costs)
        tps_reg_costs.append(reg_costs)
        print regtype, "time elapsed", time.time() - start_time

    np.set_printoptions(suppress=True)
    
    print ""
    print "tps_costs"
    print args.regtypes
    print np.array(tps_costs).T
    
    print ""
    print "tps_reg_costs"
    print args.regtypes
    print np.array(tps_reg_costs).T

Example 79

Project: lfd
Source File: planning.py
View license
def plan_follow_trajs(robot, manip_name, ee_link_names, ee_trajs, old_traj, 
                     no_collision_cost_first=False, use_collision_cost=True, start_fixed=False, joint_vel_limits=None,
                     beta_pos=settings.BETA_POS, beta_rot=settings.BETA_ROT, gamma=settings.GAMMA):
    orig_dof_inds = robot.GetActiveDOFIndices()
    orig_dof_vals = robot.GetDOFValues()
    
    n_steps = len(ee_trajs[0])
    dof_inds = sim_util.dof_inds_from_name(robot, manip_name)
    assert old_traj.shape[0] == n_steps
    assert old_traj.shape[1] == len(dof_inds)
    assert len(ee_link_names) == len(ee_trajs)
        
    if no_collision_cost_first:
        init_traj, _, _ = plan_follow_trajs(robot, manip_name, ee_link_names, ee_trajs, old_traj, 
                                        no_collision_cost_first=False, use_collision_cost=False, start_fixed=start_fixed, joint_vel_limits=joint_vel_limits,
                                        beta_pos = beta_pos, beta_rot = beta_rot, gamma = gamma)
    else:
        init_traj = old_traj.copy()

    if start_fixed:
        init_traj = np.r_[robot.GetDOFValues(dof_inds)[None,:], init_traj[1:]]
        sim_util.unwrap_in_place(init_traj, dof_inds)
        init_traj += robot.GetDOFValues(dof_inds) - init_traj[0,:]

    request = {
        "basic_info" : {
            "n_steps" : n_steps,
            "manip" : manip_name,
            "start_fixed" : start_fixed
        },
        "costs" : [
        {
            "type" : "joint_vel",
            "params": {"coeffs" : [gamma/(n_steps-1)]}
        },            
        ],
        "constraints" : [
        ],
        "init_info" : {
            "type":"given_traj",
            "data":[x.tolist() for x in init_traj]
        }
    }
    
    if use_collision_cost:
        request["costs"].append(
            {
                "type" : "collision",
                "params" : {
                  "continuous" : True,
                  "coeffs" : [1000],  # penalty coefficients. list of length one is automatically expanded to a list of length n_timesteps
                  "dist_pen" : [0.025]  # robot-obstacle distance that penalty kicks in. expands to length n_timesteps
                }
            })
    
    if joint_vel_limits is not None:
        request["constraints"].append(
             {
                "type" : "joint_vel_limits",
                "params": {"vals" : joint_vel_limits,
                           "first_step" : 0,
                           "last_step" : n_steps-1
                           }
              })

    for (ee_link_name, ee_traj) in zip(ee_link_names, ee_trajs):
        poses = [openravepy.poseFromMatrix(hmat) for hmat in ee_traj]
        for (i_step,pose) in enumerate(poses):
            if start_fixed and i_step == 0:
                continue
            request["costs"].append(
                {"type":"pose",
                 "params":{
                    "xyz":pose[4:7].tolist(),
                    "wxyz":pose[0:4].tolist(),
                    "link":ee_link_name,
                    "timestep":i_step,
                    "pos_coeffs":[np.sqrt(beta_pos/n_steps)]*3,
                    "rot_coeffs":[np.sqrt(beta_rot/n_steps)]*3
                 }
                })

    s = json.dumps(request)
    with openravepy.RobotStateSaver(robot):
        orig_dof_vals
        with util.suppress_stdout():
            prob = trajoptpy.ConstructProblem(s, robot.GetEnv()) # create object that stores optimization problem
            result = trajoptpy.OptimizeProblem(prob) # do optimization
    traj = result.GetTraj()
    
    pose_costs = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts() if cost_type == "pose"])
    pose_err = []
    with openravepy.RobotStateSaver(robot):
        for (ee_link_name, ee_traj) in zip(ee_link_names, ee_trajs):
            ee_link = robot.GetLink(ee_link_name)
            for (i_step, hmat) in enumerate(ee_traj):
                if start_fixed and i_step == 0:
                    continue
                robot.SetDOFValues(traj[i_step], dof_inds)
                new_hmat = ee_link.GetTransform()
                pose_err.append(openravepy.poseFromMatrix(mu.invertHmat(hmat).dot(new_hmat)))
    pose_err = np.asarray(pose_err)
    pose_costs2 = (beta_rot/n_steps) * np.square(pose_err[:,1:4]).sum() + (beta_pos/n_steps) * np.square(pose_err[:,4:7]).sum()

    joint_vel_cost = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts() if cost_type == "joint_vel"])
    joint_vel_err = np.diff(traj, axis=0)
    joint_vel_cost2 = (gamma/(n_steps-1)) * np.square(joint_vel_err).sum()
    sim_util.unwrap_in_place(traj, dof_inds)
    joint_vel_err = np.diff(traj, axis=0)
    
    collision_costs = [cost_val for (cost_type, cost_val) in result.GetCosts() if "collision" in cost_type]
    if len(collision_costs) > 0:
        collision_err = np.asarray(collision_costs)
        collision_costs = np.sum(collision_costs)

    obj_value = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts()])
    
    print "{:>15} | {:>10} | {:>10}".format("", "trajopt", "computed")
    print "{:>15} | {:>10}".format("COSTS", "-"*23)
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("joint_vel", joint_vel_cost, joint_vel_cost2)
    if np.isscalar(collision_costs):
        print "{:>15} | {:>10,.4} | {:>10}".format("collision(s)", collision_costs, "-")
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("pose(s)", pose_costs, pose_costs2)
    print "{:>15} | {:>10,.4} | {:>10}".format("total_obj", obj_value, "-")
    print ""

    print "{:>15} | {:>10} | {:>10}".format("", "abs min", "abs max")
    print "{:>15} | {:>10}".format("ERRORS", "-"*23)
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("joint_vel (deg)", np.rad2deg(np.abs(joint_vel_err).min()), np.rad2deg(np.abs(joint_vel_err).max()))
    if np.isscalar(collision_costs):
        print "{:>15} | {:>10,.4} | {:>10,.4}".format("collision(s)", np.abs(-collision_err).min(), np.abs(-collision_err).max())
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("rot pose(s)", np.abs(pose_err[:,1:4]).min(), np.abs(pose_err[:,1:4]).max())
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("trans pose(s)", np.abs(pose_err[:,4:7]).min(), np.abs(pose_err[:,4:7]).max())
    print ""

    # make sure this function doesn't change state of the robot
    assert not np.any(orig_dof_inds - robot.GetActiveDOFIndices())
    assert not np.any(orig_dof_vals - robot.GetDOFValues())
    
    return traj, obj_value, pose_costs

Example 80

Project: lfd
Source File: planning.py
View license
def plan_follow_finger_pts_trajs(robot, manip_name, flr2finger_link_names, flr2finger_rel_pts, flr2finger_pts_trajs, old_traj, 
                                no_collision_cost_first=False, use_collision_cost=True, start_fixed=False, joint_vel_limits=None,
                                beta_pos=settings.BETA_POS, gamma=settings.GAMMA):
    orig_dof_inds = robot.GetActiveDOFIndices()
    orig_dof_vals = robot.GetDOFValues()
    
    n_steps = old_traj.shape[0]
    dof_inds = sim_util.dof_inds_from_name(robot, manip_name)
    assert old_traj.shape[1] == len(dof_inds)
    for flr2finger_pts_traj in flr2finger_pts_trajs:
        for finger_pts_traj in flr2finger_pts_traj.values():
            assert len(finger_pts_traj)== n_steps
    assert len(flr2finger_link_names) == len(flr2finger_pts_trajs)
    
    if no_collision_cost_first:
        init_traj, _ = plan_follow_finger_pts_trajs(robot, manip_name, flr2finger_link_names, flr2finger_rel_pts, flr2finger_pts_trajs, old_traj, 
                                                   no_collision_cost_first=False, use_collision_cost=False, start_fixed=start_fixed, joint_vel_limits=joint_vel_limits,
                                                   beta_pos = beta_pos, gamma = gamma)
    else:
        init_traj = old_traj.copy()

    if start_fixed:
        init_traj = np.r_[robot.GetDOFValues(dof_inds)[None,:], init_traj[1:]]
        sim_util.unwrap_in_place(init_traj, dof_inds)
        init_traj += robot.GetDOFValues(dof_inds) - init_traj[0,:]

    request = {
        "basic_info" : {
            "n_steps" : n_steps,
            "manip" : manip_name,
            "start_fixed" : start_fixed
        },
        "costs" : [
        {
            "type" : "joint_vel",
            "params": {"coeffs" : [gamma/(n_steps-1)]}
        },            
        ],
        "constraints" : [
        ],
        "init_info" : {
            "type":"given_traj",
            "data":[x.tolist() for x in init_traj]
        }
    }
    
    if use_collision_cost:
        request["costs"].append(
            {
                "type" : "collision",
                "params" : {
                  "continuous" : True,
                  "coeffs" : [1000],  # penalty coefficients. list of length one is automatically expanded to a list of length n_timesteps
                  "dist_pen" : [0.025]  # robot-obstacle distance that penalty kicks in. expands to length n_timesteps
                }
            })
    
    if joint_vel_limits is not None:
        request["constraints"].append(
             {
                "type" : "joint_vel_limits",
                "params": {"vals" : joint_vel_limits,
                           "first_step" : 0,
                           "last_step" : n_steps-1
                           }
              })

    for (flr2finger_link_name, flr2finger_pts_traj) in zip(flr2finger_link_names, flr2finger_pts_trajs):
        for finger_lr, finger_link_name in flr2finger_link_name.items():
            finger_rel_pts = flr2finger_rel_pts[finger_lr]
            finger_pts_traj = flr2finger_pts_traj[finger_lr]
            for (i_step, finger_pts) in enumerate(finger_pts_traj):
                if start_fixed and i_step == 0:
                    continue
                request["costs"].append(
                    {"type":"rel_pts",
                     "params":{
                        "xyzs":finger_pts.tolist(),
                        "rel_xyzs":finger_rel_pts.tolist(),
                        "link":finger_link_name,
                        "timestep":i_step,
                        "pos_coeffs":[np.sqrt(beta_pos/n_steps)]*4, # there is a coefficient for each of the 4 points
                     }
                    })

    s = json.dumps(request)
    with openravepy.RobotStateSaver(robot):
        with util.suppress_stdout():
            prob = trajoptpy.ConstructProblem(s, robot.GetEnv()) # create object that stores optimization problem
            result = trajoptpy.OptimizeProblem(prob) # do optimization

    traj = result.GetTraj() 


    rel_pts_costs = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts() if cost_type == "rel_pts"])
    rel_pts_err = []
    with openravepy.RobotStateSaver(robot):
        for (flr2finger_link_name, flr2finger_pts_traj) in zip(flr2finger_link_names, flr2finger_pts_trajs):
            for finger_lr, finger_link_name in flr2finger_link_name.items():
                finger_link = robot.GetLink(finger_link_name)
                finger_rel_pts = flr2finger_rel_pts[finger_lr]
                finger_pts_traj = flr2finger_pts_traj[finger_lr]
                for (i_step, finger_pts) in enumerate(finger_pts_traj):
                    if start_fixed and i_step == 0:
                        continue
                    robot.SetDOFValues(traj[i_step], dof_inds)
                    new_hmat = finger_link.GetTransform()
                    rel_pts_err.append(finger_pts - (new_hmat[:3,3][None,:] + finger_rel_pts.dot(new_hmat[:3,:3].T)))
    rel_pts_err = np.concatenate(rel_pts_err, axis=0)
    rel_pts_costs2 = (beta_pos/n_steps) * np.square(rel_pts_err).sum() # TODO don't square n_steps

    joint_vel_cost = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts() if cost_type == "joint_vel"])
    joint_vel_err = np.diff(traj, axis=0)
    joint_vel_cost2 = (gamma/(n_steps-1)) * np.square(joint_vel_err).sum()
    sim_util.unwrap_in_place(traj, dof_inds)
    joint_vel_err = np.diff(traj, axis=0)

    collision_costs = [cost_val for (cost_type, cost_val) in result.GetCosts() if "collision" in cost_type]
    if len(collision_costs) > 0:
        collision_err = np.asarray(collision_costs)
        collision_costs = np.sum(collision_costs)

    obj_value = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts()])
    
    print "{:>15} | {:>10} | {:>10}".format("", "trajopt", "computed")
    print "{:>15} | {:>10}".format("COSTS", "-"*23)
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("joint_vel", joint_vel_cost, joint_vel_cost2)
    if np.isscalar(collision_costs):
        print "{:>15} | {:>10,.4} | {:>10}".format("collision(s)", collision_costs, "-")
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("rel_pts(s)", rel_pts_costs, rel_pts_costs2)
    print "{:>15} | {:>10,.4} | {:>10}".format("total_obj", obj_value, "-")
    print ""

    print "{:>15} | {:>10} | {:>10}".format("", "abs min", "abs max")
    print "{:>15} | {:>10}".format("ERRORS", "-"*23)
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("joint_vel (deg)", np.rad2deg(np.abs(joint_vel_err).min()), np.rad2deg(np.abs(joint_vel_err).max()))
    if np.isscalar(collision_costs):
        print "{:>15} | {:>10,.4} | {:>10,.4}".format("collision(s)", np.abs(-collision_err).min(), np.abs(-collision_err).max())
    print "{:>15} | {:>10,.4} | {:>10,.4}".format("rel_pts(s)", np.abs(rel_pts_err).min(), np.abs(rel_pts_err).max())
    print ""

    # make sure this function doesn't change state of the robot
    assert not np.any(orig_dof_inds - robot.GetActiveDOFIndices())
    assert not np.any(orig_dof_vals - robot.GetDOFValues())
    
    obj_value = np.sum([cost_val for (cost_type, cost_val) in result.GetCosts()])
    return traj, obj_value, rel_pts_costs

Example 81

Project: scikit-image
Source File: draw3d.py
View license
def ellipsoid(a, b, c, spacing=(1., 1., 1.), levelset=False):
    """
    Generates ellipsoid with semimajor axes aligned with grid dimensions
    on grid with specified `spacing`.

    Parameters
    ----------
    a : float
        Length of semimajor axis aligned with x-axis.
    b : float
        Length of semimajor axis aligned with y-axis.
    c : float
        Length of semimajor axis aligned with z-axis.
    spacing : tuple of floats, length 3
        Spacing in (x, y, z) spatial dimensions.
    levelset : bool
        If True, returns the level set for this ellipsoid (signed level
        set about zero, with positive denoting interior) as np.float64.
        False returns a binarized version of said level set.

    Returns
    -------
    ellip : (N, M, P) array
        Ellipsoid centered in a correctly sized array for given `spacing`.
        Boolean dtype unless `levelset=True`, in which case a float array is
        returned with the level set above 0.0 representing the ellipsoid.

    """
    if (a <= 0) or (b <= 0) or (c <= 0):
        raise ValueError('Parameters a, b, and c must all be > 0')

    offset = np.r_[1, 1, 1] * np.r_[spacing]

    # Calculate limits, and ensure output volume is odd & symmetric
    low = np.ceil((- np.r_[a, b, c] - offset))
    high = np.floor((np.r_[a, b, c] + offset + 1))

    for dim in range(3):
        if (high[dim] - low[dim]) % 2 == 0:
            low[dim] -= 1
        num = np.arange(low[dim], high[dim], spacing[dim])
        if 0 not in num:
            low[dim] -= np.max(num[num < 0])

    # Generate (anisotropic) spatial grid
    x, y, z = np.mgrid[low[0]:high[0]:spacing[0],
                       low[1]:high[1]:spacing[1],
                       low[2]:high[2]:spacing[2]]

    if not levelset:
        arr = ((x / float(a)) ** 2 +
               (y / float(b)) ** 2 +
               (z / float(c)) ** 2) <= 1
    else:
        arr = ((x / float(a)) ** 2 +
               (y / float(b)) ** 2 +
               (z / float(c)) ** 2) - 1

    return arr

Example 82

Project: scikit-learn
Source File: test_iforest.py
View license
def test_iforest_performance():
    """Test Isolation Forest performs well"""

    # Generate train/test data
    rng = check_random_state(2)
    X = 0.3 * rng.randn(120, 2)
    X_train = np.r_[X + 2, X - 2]
    X_train = X[:100]

    # Generate some abnormal novel observations
    X_outliers = rng.uniform(low=-4, high=4, size=(20, 2))
    X_test = np.r_[X[100:], X_outliers]
    y_test = np.array([0] * 20 + [1] * 20)

    # fit the model
    clf = IsolationForest(max_samples=100, random_state=rng).fit(X_train)

    # predict scores (the lower, the more normal)
    y_pred = - clf.decision_function(X_test)

    # check that there is at most 6 errors (false positive or false negative)
    assert_greater(roc_auc_score(y_test, y_pred), 0.98)

Example 83

Project: scikit-learn
Source File: ranking.py
View license
def precision_recall_curve(y_true, probas_pred, pos_label=None,
                           sample_weight=None):
    """Compute precision-recall pairs for different probability thresholds

    Note: this implementation is restricted to the binary classification task.

    The precision is the ratio ``tp / (tp + fp)`` where ``tp`` is the number of
    true positives and ``fp`` the number of false positives. The precision is
    intuitively the ability of the classifier not to label as positive a sample
    that is negative.

    The recall is the ratio ``tp / (tp + fn)`` where ``tp`` is the number of
    true positives and ``fn`` the number of false negatives. The recall is
    intuitively the ability of the classifier to find all the positive samples.

    The last precision and recall values are 1. and 0. respectively and do not
    have a corresponding threshold.  This ensures that the graph starts on the
    x axis.

    Read more in the :ref:`User Guide <precision_recall_f_measure_metrics>`.

    Parameters
    ----------
    y_true : array, shape = [n_samples]
        True targets of binary classification in range {-1, 1} or {0, 1}.

    probas_pred : array, shape = [n_samples]
        Estimated probabilities or decision function.

    pos_label : int or str, default=None
        The label of the positive class

    sample_weight : array-like of shape = [n_samples], optional
        Sample weights.

    Returns
    -------
    precision : array, shape = [n_thresholds + 1]
        Precision values such that element i is the precision of
        predictions with score >= thresholds[i] and the last element is 1.

    recall : array, shape = [n_thresholds + 1]
        Decreasing recall values such that element i is the recall of
        predictions with score >= thresholds[i] and the last element is 0.

    thresholds : array, shape = [n_thresholds <= len(np.unique(probas_pred))]
        Increasing thresholds on the decision function used to compute
        precision and recall.

    Examples
    --------
    >>> import numpy as np
    >>> from sklearn.metrics import precision_recall_curve
    >>> y_true = np.array([0, 0, 1, 1])
    >>> y_scores = np.array([0.1, 0.4, 0.35, 0.8])
    >>> precision, recall, thresholds = precision_recall_curve(
    ...     y_true, y_scores)
    >>> precision  # doctest: +ELLIPSIS
    array([ 0.66...,  0.5       ,  1.        ,  1.        ])
    >>> recall
    array([ 1. ,  0.5,  0.5,  0. ])
    >>> thresholds
    array([ 0.35,  0.4 ,  0.8 ])

    """
    fps, tps, thresholds = _binary_clf_curve(y_true, probas_pred,
                                             pos_label=pos_label,
                                             sample_weight=sample_weight)

    precision = tps / (tps + fps)
    recall = tps / tps[-1]

    # stop when full recall attained
    # and reverse the outputs so recall is decreasing
    last_ind = tps.searchsorted(tps[-1])
    sl = slice(last_ind, None, -1)
    return np.r_[precision[sl], 1], np.r_[recall[sl], 0], thresholds[sl]

Example 84

Project: scikit-learn
Source File: test_lof.py
View license
def test_lof_performance():
    # Generate train/test data
    rng = check_random_state(2)
    X = 0.3 * rng.randn(120, 2)
    X_train = np.r_[X + 2, X - 2]
    X_train = X[:100]

    # Generate some abnormal novel observations
    X_outliers = rng.uniform(low=-4, high=4, size=(20, 2))
    X_test = np.r_[X[100:], X_outliers]
    y_test = np.array([0] * 20 + [1] * 20)

    # fit the model
    clf = neighbors.LocalOutlierFactor().fit(X_train)

    # predict scores (the lower, the more normal)
    y_pred = -clf._decision_function(X_test)

    # check that roc_auc is good
    assert_greater(roc_auc_score(y_test, y_pred), .99)

Example 85

Project: scikit-learn
Source File: test_svm.py
View license
def test_oneclass_decision_function():
    # Test OneClassSVM decision function
    clf = svm.OneClassSVM()
    rnd = check_random_state(2)

    # Generate train data
    X = 0.3 * rnd.randn(100, 2)
    X_train = np.r_[X + 2, X - 2]

    # Generate some regular novel observations
    X = 0.3 * rnd.randn(20, 2)
    X_test = np.r_[X + 2, X - 2]
    # Generate some abnormal novel observations
    X_outliers = rnd.uniform(low=-4, high=4, size=(20, 2))

    # fit the model
    clf = svm.OneClassSVM(nu=0.1, kernel="rbf", gamma=0.1)
    clf.fit(X_train)

    # predict things
    y_pred_test = clf.predict(X_test)
    assert_greater(np.mean(y_pred_test == 1), .9)
    y_pred_outliers = clf.predict(X_outliers)
    assert_greater(np.mean(y_pred_outliers == -1), .9)
    dec_func_test = clf.decision_function(X_test)
    assert_array_equal((dec_func_test > 0).ravel(), y_pred_test == 1)
    dec_func_outliers = clf.decision_function(X_outliers)
    assert_array_equal((dec_func_outliers > 0).ravel(), y_pred_outliers == 1)

Example 86

Project: scipy
Source File: _bsplines.py
View license
def make_interp_spline(x, y, k=3, t=None, bc_type=None, axis=0,
                       check_finite=True):
    """Compute the (coefficients of) interpolating B-spline.

    Parameters
    ----------
    x : array_like, shape (n,)
        Abscissas.
    y : array_like, shape (n, ...)
        Ordinates.
    k : int, optional
        B-spline degree. Default is cubic, k=3.
    t : array_like, shape (nt + k + 1,), optional.
        Knots.
        The number of knots needs to agree with the number of datapoints and
        the number of derivatives at the edges. Specifically, ``nt - n`` must
        equal ``len(deriv_l) + len(deriv_r)``.
    bc_type : 2-tuple or None
        Boundary conditions.
        Default is None, which means choosing the boundary conditions
        automatically. Otherwise, it must be a length-two tuple where the first
        element sets the boundary conditions at ``x[0]`` and the second
        element sets the boundary conditions at ``x[-1]``. Each of these must
        be an iterable of pairs ``(order, value)`` which gives the values of
        derivatives of specified orders at the given edge of the interpolation
        interval.
    axis : int, optional
        Interpolation axis. Default is 0.
    check_finite : bool, optional
        Whether to check that the input arrays contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.
        Default is True.

    Returns
    -------
    b : a BSpline object of the degree ``k`` and with knots ``t``.

    Examples
    --------

    Use cubic interpolation on Chebyshev nodes:

    >>> def cheb_nodes(N):
    ...     jj = 2.*np.arange(N) + 1
    ...     x = np.cos(np.pi * jj / 2 / N)[::-1]
    ...     return x

    >>> x = cheb_nodes(20)
    >>> y = np.sqrt(1 - x**2)

    >>> from scipy.interpolate import BSpline, make_interp_spline
    >>> b = make_interp_spline(x, y)
    >>> np.allclose(b(x), y)
    True

    Note that the default is a cubic spline with a not-a-knot boundary condition

    >>> b.k
    3

    Here we use a 'natural' spline, with zero 2nd derivatives at edges:

    >>> l, r = [(2, 0)], [(2, 0)]
    >>> b_n = make_interp_spline(x, y, bc_type=(l, r))
    >>> np.allclose(b_n(x), y)
    True
    >>> x0, x1 = x[0], x[-1]
    >>> np.allclose([b_n(x0, 2), b_n(x1, 2)], [0, 0])
    True

    Interpolation of parametric curves is also supported. As an example, we
    compute a discretization of a snail curve in polar coordinates

    >>> phi = np.linspace(0, 2.*np.pi, 40)
    >>> r = 0.3 + np.cos(phi)
    >>> x, y = r*np.cos(phi), r*np.sin(phi)  # convert to Cartesian coordinates

    Build an interpolating curve, parameterizing it by the angle

    >>> from scipy.interpolate import make_interp_spline
    >>> spl = make_interp_spline(phi, np.c_[x, y])

    Evaluate the interpolant on a finer grid (note that we transpose the result
    to unpack it into a pair of x- and y-arrays)

    >>> phi_new = np.linspace(0, 2.*np.pi, 100)
    >>> x_new, y_new = spl(phi_new).T

    Plot the result

    >>> import matplotlib.pyplot as plt
    >>> plt.plot(x, y, 'o')
    >>> plt.plot(x_new, y_new, '-')
    >>> plt.show()

    See Also
    --------
    BSpline : base class representing the B-spline objects
    CubicSpline : a cubic spline in the polynomial basis
    make_lsq_spline : a similar factory function for spline fitting
    UnivariateSpline : a wrapper over FITPACK spline fitting routines
    splrep : a wrapper over FITPACK spline fitting routines

    """
    if bc_type is None:
        bc_type = (None, None)
    deriv_l, deriv_r = bc_type

    # special-case k=0 right away
    if k == 0:
        if any(_ is not None for _ in (t, deriv_l, deriv_r)):
            raise ValueError("Too much info for k=0: t and bc_type can only "
                             "be None.")
        x = _as_float_array(x, check_finite)
        t = np.r_[x, x[-1]]
        c = np.asarray(y)
        c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
        return BSpline.construct_fast(t, c, k, axis=axis)

    # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
    if k == 1 and t is None:
        if not (deriv_l is None and deriv_r is None):
            raise ValueError("Too much info for k=1: bc_type can only be None.")
        x = _as_float_array(x, check_finite)
        t = np.r_[x[0], x, x[-1]]
        c = np.asarray(y)
        c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
        return BSpline.construct_fast(t, c, k, axis=axis)

    # come up with a sensible knot vector, if needed
    if t is None:
        if deriv_l is None and deriv_r is None:
            if k == 2:
                # OK, it's a bit ad hoc: Greville sites + omit
                # 2nd and 2nd-to-last points, a la not-a-knot
                t = (x[1:] + x[:-1]) / 2.
                t = np.r_[(x[0],)*(k+1),
                           t[1:-1],
                           (x[-1],)*(k+1)]
            else:
                t = _not_a_knot(x, k)
        else:
            t = _augknt(x, k)

    x = _as_float_array(x, check_finite)
    y = _as_float_array(y, check_finite)
    t = _as_float_array(t, check_finite)
    k = int(k)

    axis = axis % y.ndim
    y = np.rollaxis(y, axis)    # now internally interp axis is zero

    if x.ndim != 1 or np.any(x[1:] <= x[:-1]):
        raise ValueError("Expect x to be a 1-D sorted array_like.")
    if k < 0:
        raise ValueError("Expect non-negative k.")
    if t.ndim != 1 or np.any(t[1:] < t[:-1]):
        raise ValueError("Expect t to be a 1-D sorted array_like.")
    if x.size != y.shape[0]:
        raise ValueError('x and y are incompatible.')
    if t.size < x.size + k + 1:
        raise ValueError('Got %d knots, need at least %d.' %
                         (t.size, x.size + k + 1))
    if (x[0] < t[k]) or (x[-1] > t[-k]):
        raise ValueError('Out of bounds w/ x = %s.' % x)

    # Here : deriv_l, r = [(nu, value), ...]
    if deriv_l is not None:
        deriv_l_ords, deriv_l_vals = zip(*deriv_l)
    else:
        deriv_l_ords, deriv_l_vals = [], []
    deriv_l_ords, deriv_l_vals = np.atleast_1d(deriv_l_ords, deriv_l_vals)
    nleft = deriv_l_ords.shape[0]

    if deriv_r is not None:
        deriv_r_ords, deriv_r_vals = zip(*deriv_r)
    else:
        deriv_r_ords, deriv_r_vals = [], []
    deriv_r_ords, deriv_r_vals = np.atleast_1d(deriv_r_ords, deriv_r_vals)
    nright = deriv_r_ords.shape[0]

    # have `n` conditions for `nt` coefficients; need nt-n derivatives
    n = x.size
    nt = t.size - k - 1

    if nt - n != nleft + nright:
        raise ValueError("number of derivatives at boundaries.")

    # set up the LHS: the collocation matrix + derivatives at boundaries
    kl = ku = k
    ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
    _bspl._colloc(x, t, k, ab, offset=nleft)
    if nleft > 0:
        _bspl._handle_lhs_derivatives(t, k, x[0], ab, kl, ku, deriv_l_ords)
    if nright > 0:
        _bspl._handle_lhs_derivatives(t, k, x[-1], ab, kl, ku, deriv_r_ords,
                                offset=nt-nright)

    # set up the RHS: values to interpolate (+ derivative values, if any)
    extradim = prod(y.shape[1:])
    rhs = np.empty((nt, extradim), dtype=y.dtype)
    if nleft > 0:
        rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
    rhs[nleft:nt - nright] = y.reshape(-1, extradim)
    if nright > 0:
        rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)

    # solve Ab @ x = rhs; this is the relevant part of linalg.solve_banded
    if check_finite:
        ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
    gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
    lu, piv, c, info = gbsv(kl, ku, ab, rhs,
            overwrite_ab=True, overwrite_b=True)

    if info > 0:
        raise LinAlgError("Collocation matix is singular.")
    elif info < 0:
        raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)

    c = np.ascontiguousarray(c.reshape((nt,) + y.shape[1:]))
    return BSpline.construct_fast(t, c, k, axis=axis)

Example 87

Project: scipy
Source File: _discrete_distns.py
View license
    def _entropy(self, n, p):
        k = np.r_[0:n + 1]
        vals = self._pmf(k, n, p)
        return np.sum(entr(vals), axis=0)

Example 88

Project: scipy
Source File: _discrete_distns.py
View license
    def _entropy(self, M, n, N):
        k = np.r_[N - (M - n):min(n, N) + 1]
        vals = self.pmf(k, M, n, N)
        return np.sum(entr(vals), axis=0)

Example 89

Project: sfepy
Source File: lcbc_operators.py
View license
    def finalize(self):
        """
        Call this after all LCBCs of the variable have been added.

        Initializes the global column indices and DOF counts.
        """
        keys = self.variables.adi.var_names
        self.n_master = {}.fromkeys(keys, 0)
        self.n_slave = {}.fromkeys(keys, 0)
        self.n_new = {}.fromkeys(keys, 0)
        ics = {}
        for ii, op in enumerate(self):
            self.n_master[op.var_names[0]] += op.n_mdof

            if op.var_names[1] is not None:
                self.n_slave[op.var_names[1]] += op.n_sdof

            self.n_new[op.var_names[0]] += op.n_new_dof
            ics.setdefault(op.var_names[0], []).append((ii, op.n_new_dof))

        self.ics = {}
        for key, val in six.iteritems(ics):
            iis, ics = zip(*val)
            self.ics[key] = (iis, nm.cumsum(nm.r_[0, ics]))

        self.n_free = {}
        self.n_active = {}
        n_dof = self.variables.adi.n_dof
        for key in six.iterkeys(self.n_master):
            self.n_free[key] = n_dof[key] - self.n_master[key]
            self.n_active[key] = self.n_free[key] + self.n_new[key]

        def _dict_to_di(name, dd):
            di = DofInfo(name)
            for key in keys:
                val = dd[key]
                di.append_raw(key, val)
            return di

        self.lcdi = _dict_to_di('lcbc_active_state_dof_info', self.n_active)
        self.fdi = _dict_to_di('free_dof_info', self.n_free)
        self.ndi = _dict_to_di('new_dof_info', self.n_new)

Example 90

Project: sfepy
Source File: refine.py
View license
def refine_2_3(mesh_in):
    """
    Refines mesh out of triangles by cutting cutting each edge in half
    and making 4 new finer triangles out of one coarser one.
    """
    cmesh = mesh_in.cmesh

    # Unique edge centres.
    e_centres = cmesh.get_centroids(cmesh.dim - 1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres]

    o1 = mesh_in.n_nod

    cc = cmesh.get_conn(cmesh.dim, cmesh.dim - 1)

    conn = mesh_in.get_conn('2_3')
    n_el = conn.shape[0]

    e_nodes = cc.indices.reshape((n_el, 3)) + o1

    c = nm.c_[conn, e_nodes].T

    new_conn = nm.vstack([c[0], c[3], c[5],
                          c[3], c[4], c[5],
                          c[1], c[4], c[3],
                          c[2], c[5], c[4]]).T
    new_conn = new_conn.reshape((4 * n_el, 3))

    new_mat_id = cmesh.cell_groups.repeat(4)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
                          [new_mat_id], mesh_in.descs )

    return mesh

Example 91

Project: sfepy
Source File: refine.py
View license
def refine_2_4(mesh_in):
    """
    Refines mesh out of quadrilaterals by cutting cutting each edge in
    half and making 4 new finer quadrilaterals out of one coarser one.
    """
    cmesh = mesh_in.cmesh

    # Unique edge centres.
    e_centres = cmesh.get_centroids(cmesh.dim - 1)

    # Unique element centres.
    centres = cmesh.get_centroids(cmesh.dim)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]

    cc = cmesh.get_conn(cmesh.dim, cmesh.dim - 1)

    conn = mesh_in.get_conn('2_4')
    n_el = conn.shape[0]

    e_nodes = cc.indices.reshape((n_el, 4)) + o1
    nodes = nm.arange(n_el) + o2

    c = nm.c_[conn, e_nodes, nodes].T

    new_conn = nm.vstack([c[0], c[4], c[8], c[7],
                          c[1], c[5], c[8], c[4],
                          c[2], c[6], c[8], c[5],
                          c[3], c[7], c[8], c[6]]).T
    new_conn = new_conn.reshape((4 * n_el, 4))

    new_mat_id = cmesh.cell_groups.repeat(4)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
                          [new_mat_id], mesh_in.descs )

    return mesh

Example 92

Project: sfepy
Source File: refine.py
View license
def refine_3_4(mesh_in):
    """
    Refines tetrahedra by cutting each edge in half and making 8 new
    finer tetrahedra out of one coarser one. Old nodal coordinates come
    first in `coors`, then the new ones. The new tetrahedra are similar
    to the old one, no degeneration is supposed to occur as at most 3
    congruence classes of tetrahedra appear, even when re-applied
    iteratively (provided that `conns` are not modified between two
    applications - ordering of vertices in tetrahedra matters not only
    for positivity of volumes).

    References:

    - Juergen Bey: Simplicial grid refinement: on Freudenthal s algorithm and 
      the optimal number of congruence classes, Numer.Math. 85 (2000), 
      no. 1, 1--29, or
    - Juergen Bey: Tetrahedral grid refinement, Computing 55 (1995), 
      no. 4, 355--378, or
      http://citeseer.ist.psu.edu/bey95tetrahedral.html
    """
    cmesh = mesh_in.cmesh

    # Unique edge centres.
    e_centres = cmesh.get_centroids(cmesh.dim - 2)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres]

    o1 = mesh_in.n_nod

    cc = cmesh.get_conn(cmesh.dim, cmesh.dim - 2)

    conn = mesh_in.get_conn('3_4')
    n_el = conn.shape[0]

    e_nodes = cc.indices.reshape((n_el, 6)) + o1

    c = nm.c_[conn, e_nodes].T

    new_conn = nm.vstack([c[0], c[4], c[6], c[7],
                          c[4], c[1], c[5], c[8],
                          c[6], c[5], c[2], c[9],
                          c[7], c[8], c[9], c[3],
                          c[4], c[6], c[7], c[8],
                          c[4], c[6], c[8], c[5],
                          c[6], c[7], c[8], c[9],
                          c[6], c[5], c[9], c[8]]).T
    new_conn = new_conn.reshape((8 * n_el, 4))

    new_mat_id = cmesh.cell_groups.repeat(8)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
                          [new_mat_id], mesh_in.descs )

    return mesh

Example 93

Project: sfepy
Source File: refine.py
View license
def refine_3_8(mesh_in):
    """
    Refines hexahedral mesh by cutting cutting each edge in half and
    making 8 new finer hexahedrons out of one coarser one.
    """
    cmesh = mesh_in.cmesh

    # Unique edge centres.
    e_centres = cmesh.get_centroids(cmesh.dim - 2)

    # Unique face centres.
    f_centres = cmesh.get_centroids(cmesh.dim - 1)

    # Unique element centres.
    centres = cmesh.get_centroids(cmesh.dim)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, f_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]
    o3 = o2 + f_centres.shape[0]

    ecc = cmesh.get_conn(cmesh.dim, cmesh.dim - 2)
    fcc = cmesh.get_conn(cmesh.dim, cmesh.dim - 1)

    conn = mesh_in.get_conn('3_8')
    n_el = conn.shape[0]

    st = nm.vstack

    e_nodes = ecc.indices.reshape((n_el, 12)) + o1
    f_nodes = fcc.indices.reshape((n_el, 6)) + o2
    nodes = nm.arange(n_el) + o3

    c = nm.c_[conn, e_nodes, f_nodes, nodes].T

    new_conn = st([c[0], c[8], c[20], c[11], c[16], c[22], c[26], c[21],
                   c[1], c[9], c[20], c[8], c[17], c[24], c[26], c[22],
                   c[2], c[10], c[20], c[9], c[18], c[25], c[26], c[24],
                   c[3], c[11], c[20], c[10], c[19], c[21], c[26], c[25],
                   c[4], c[15], c[23], c[12], c[16], c[21], c[26], c[22],
                   c[5], c[12], c[23], c[13], c[17], c[22], c[26], c[24],
                   c[6], c[13], c[23], c[14], c[18], c[24], c[26], c[25],
                   c[7], c[14], c[23], c[15], c[19], c[25], c[26], c[21]]).T
    new_conn = new_conn.reshape((8 * n_el, 8))

    new_mat_id = cmesh.cell_groups.repeat(8)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, [new_conn],
                          [new_mat_id], mesh_in.descs )

    return mesh

Example 94

Project: sfepy
Source File: test_poly_spaces.py
View license
    def test_continuity(self):
        ok = True
        orders = {'2_3' : 3, '2_4' : 3, '3_4' : 4, '3_8' : 3}

        bads = []
        bad_families = set()
        for (geom, poly_space_base, qp_weights, mesh, ir, ic,
             field, ps, rrc, rcell, crc, ccell, vec,
             edofs, fdofs) in _gen_common_data(orders, self.gels, self.report):

            if poly_space_base == 'lagrange':
                rbf = ps.eval_base(rrc)
                cbf = ps.eval_base(crc)

            else:
                rbf = ps.eval_base(rrc, ori=field.ori[:1])
                cbf = ps.eval_base(crc, ori=field.ori[1:])

            dofs = nm.r_[edofs, fdofs]

            res = nm.zeros((2, dofs.shape[0]), dtype=nm.int32)
            res[0, :] = dofs
            for ii, ip in enumerate(dofs):
                vec.fill(0.0)
                vec[ip] = 1.0

                evec = vec[field.econn]

                rvals = nm.dot(rbf, evec[rcell])
                cvals = nm.dot(cbf, evec[ccell])

                _ok = nm.allclose(rvals, cvals, atol=1e-14, rtol=0.0)
                res[1, ii] = _ok
                if not _ok:
                    bads.append([geom, poly_space_base, ir, ic, ip])
                    bad_families.add((geom, poly_space_base))

                ok = ok and _ok

            self.report('results (dofs, status: 1 ok, 0 failure):\n%s' % res)

        if not ok:
            self.report('continuity errors:\n', bads)
            self.report('%d in total!' % len(bads))
            self.report('continuity errors occurred in these spaces:\n',
                        bad_families)

        return ok

Example 95

Project: sfepy
Source File: test_poly_spaces.py
View license
    def test_gradients(self):
        from sfepy.discrete.fem.mappings import VolumeMapping

        ok = True
        orders = {'2_3' : 3, '2_4' : 3, '3_4' : 4, '3_8' : 3}

        bads = []
        bad_families = set()
        for (geom, poly_space_base, qp_weights, mesh, ir, ic,
             field, ps, rrc, rcell, crc, ccell, vec,
             edofs, fdofs) in _gen_common_data(orders, self.gels, self.report):
            gel = self.gels[geom]
            conn = mesh.get_conn(gel.name)

            geo_ps = field.gel.poly_space
            rmapping = VolumeMapping(mesh.coors, conn[rcell:rcell+1],
                                     poly_space=geo_ps)
            rori = field.ori[:1] if field.ori is not None else None
            rvg = rmapping.get_mapping(rrc, qp_weights,
                                       poly_space=ps, ori=rori)
            rbfg = rvg.bfg

            cmapping = VolumeMapping(mesh.coors, conn[ccell:ccell+1],
                                     poly_space=geo_ps)
            cori = field.ori[1:] if field.ori is not None else None
            cvg = cmapping.get_mapping(crc, qp_weights,
                                       poly_space=ps, ori=cori)
            cbfg = cvg.bfg

            dofs = nm.r_[edofs, fdofs]

            res = nm.zeros((2, dofs.shape[0]), dtype=nm.int32)
            res[0, :] = dofs
            for ii, ip in enumerate(dofs):
                vec.fill(0.0)
                vec[ip] = 1.0

                evec = vec[field.econn]

                rvals = nm.dot(rbfg, evec[rcell])[0]
                cvals = nm.dot(cbfg, evec[ccell])[0]

                okx = nm.allclose(rvals[:, 0], cvals[:, 0],
                                  atol=1e-12, rtol=0.0)
                if gel.dim == 2:
                    oky = nm.allclose(rvals[:, 1], -cvals[:, 1],
                                      atol=1e-12, rtol=0.0)
                    _ok = okx and oky

                else:
                    oky = nm.allclose(rvals[:, 1], cvals[:, 1],
                                      atol=1e-12, rtol=0.0)
                    okz = nm.allclose(rvals[:, 2], -cvals[:, 2],
                                      atol=1e-12, rtol=0.0)
                    _ok = okx and oky and okz

                res[1, ii] = _ok
                if not _ok:
                    bads.append([geom, poly_space_base, ir, ic, ip])
                    bad_families.add((geom, poly_space_base))

                ok = ok and _ok

            self.report('results (dofs, status: 1 ok, 0 failure):\n%s' % res)

        if not ok:
            self.report('gradient continuity errors:\n', bads)
            self.report('%d in total!' % len(bads))
            self.report('gradient continuity errors occurred in these'
                        ' spaces:\n', bad_families)

        return ok

Example 96

Project: simpeg
Source File: FDEM.py
View license
def MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0):
    """
    Analytical solution for a dipole in a whole-space.

    Equation 2.57 of Ward and Hohmann

    TODOs:
        - set it up to instead take a mesh & survey
        - add E-fields
        - handle multiple frequencies
        - add divide by zero safety


    .. plot::

        from SimPEG import EM
        import matplotlib.pyplot as plt
        from scipy.constants import mu_0
        freqs = np.logspace(-2,5,100)
        Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace([0,100,0], [0,0,0], 1e-2, freqs, moment=1, orientation='Z')
        plt.loglog(freqs, np.abs(Bz.real)/mu_0, 'b')
        plt.loglog(freqs, np.abs(Bz.imag)/mu_0, 'r')
        plt.legend(('real','imag'))
        plt.show()


    """

    if not isinstance(orientation, str):
        if np.allclose(orientation, np.r_[1., 0., 0.]):
            orientation = 'X'
        elif np.allclose(orientation, np.r_[0., 1., 0.]):
            orientation = 'Y'
        elif np.allclose(orientation, np.r_[0., 0., 1.]):
            orientation = 'Z'
        else:
            raise NotImplementedError('arbitrary orientations not implemented')

    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)

    dx = XYZ[:,0]-srcLoc[0]
    dy = XYZ[:,1]-srcLoc[1]
    dz = XYZ[:,2]-srcLoc[2]

    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
    k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
    kr = k*r

    front = moment / (4.*pi * r**3.) * np.exp(-1j*kr)
    mid   = -kr**2. + 3.*1j*kr + 3.

    if orientation.upper() == 'X':
        Hx = front*( (dx/r)**2. * mid + (kr**2. - 1j*kr - 1.) )
        Hy = front*( (dx*dy/r**2.) * mid )
        Hz = front*( (dx*dz/r**2.) * mid )

    elif orientation.upper() == 'Y':
        Hx = front*( (dy*dx/r**2.) * mid )
        Hy = front*( (dy/r)**2. * mid + (kr**2. - 1j*kr - 1.) )
        Hz = front*( (dy*dz/r**2.) * mid )

    elif orientation.upper() == 'Z':
        Hx = front*( (dx*dz/r**2.) * mid )
        Hy = front*( (dy*dz/r**2.) * mid )
        Hz = front*( (dz/r)**2. * mid + (kr**2. - 1j*kr - 1.) )

    Bx = mu*Hx
    By = mu*Hy
    Bz = mu*Hz

    if Bx.ndim is 1:
        Bx = Utils.mkvc(Bx,2)

    if By.ndim is 1:
        By = Utils.mkvc(By,2)

    if Bz.ndim is 1:
        Bz = Utils.mkvc(Bz,2)

    return Bx, By, Bz

Example 97

Project: simpeg
Source File: BoundaryUtils.py
View license
def getxBCyBC_CC(mesh, alpha, beta, gamma):

    """
    This is a subfunction generating mixed-boundary condition:

    .. math::

        \\nabla \cdot \\vec{j} = -\\nabla \cdot \\vec{j}_s = q

        \\rho \\vec{j} = -\\nabla \phi

        \\alpha \phi + \\beta \\frac{\partial \phi}{\partial r} = \gamma
        \quad \\vec{r} \subset \partial \Omega

        xBC = f_1(\\alpha, \\beta, \\gamma)

        yBC = f(\\alpha, \\beta, \\gamma)


    Computes xBC and yBC for cell-centered discretizations

    """

    if mesh.dim == 1:  # 1D
        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
            raise Exception("Lenght of list, alpha should be 2")
        fCCxm, fCCxp = mesh.cellBoundaryInd
        nBC = fCCxm.sum()+fCCxp.sum()
        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]

        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]

        # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
        h_xm, h_xp = mesh.hx[0], mesh.hx[-1]

        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)

        xBC_xm = 0.5*a_xm
        xBC_xp = 0.5*a_xp/b_xp
        yBC_xm = 0.5*(1.-b_xm)
        yBC_xp = 0.5*(1.-1./b_xp)

        xBC = np.r_[xBC_xm, xBC_xp]
        yBC = np.r_[yBC_xm, yBC_xp]

    elif mesh.dim == 2:  # 2D
        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
            raise Exception("Lenght of list, alpha should be 4")

        fxm, fxp, fym, fyp = mesh.faceBoundaryInd
        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()

        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]

        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]

        h_xm = mesh.hx[0]*np.ones_like(alpha_xm)
        h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
        h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
        h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)

        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)

        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)

        xBC_xm = 0.5*a_xm
        xBC_xp = 0.5*a_xp/b_xp
        yBC_xm = 0.5*(1.-b_xm)
        yBC_xp = 0.5*(1.-1./b_xp)
        xBC_ym = 0.5*a_ym
        xBC_yp = 0.5*a_yp/b_yp
        yBC_ym = 0.5*(1.-b_ym)
        yBC_yp = 0.5*(1.-1./b_yp)

        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
                                      np.arange(mesh.nFx)[fxp]])
        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym],
                                      np.arange(mesh.nFy)[fyp]])

        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]

        xBC = np.r_[xBC_x, xBC_y]
        yBC = np.r_[yBC_x, yBC_y]

    elif mesh.dim == 3:  # 3D
        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
            raise Exception("Lenght of list, alpha should be 6")
        # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
        fxm, fxp, fym, fyp, fzm, fzp = mesh.faceBoundaryInd
        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()

        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
        alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4]
        alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5]

        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
        # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]

        h_xm = mesh.hx[0]*np.ones_like(alpha_xm)
        h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
        h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
        h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
        h_zm = mesh.hz[0]*np.ones_like(alpha_zm)
        h_zp = mesh.hz[-1]*np.ones_like(alpha_zp)

        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)

        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)

        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)

        xBC_xm = 0.5*a_xm
        xBC_xp = 0.5*a_xp/b_xp
        yBC_xm = 0.5*(1.-b_xm)
        yBC_xp = 0.5*(1.-1./b_xp)
        xBC_ym = 0.5*a_ym
        xBC_yp = 0.5*a_yp/b_yp
        yBC_ym = 0.5*(1.-b_ym)
        yBC_yp = 0.5*(1.-1./b_yp)
        xBC_zm = 0.5*a_zm
        xBC_zp = 0.5*a_zp/b_zp
        yBC_zm = 0.5*(1.-b_zm)
        yBC_zp = 0.5*(1.-1./b_zp)

        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
                                      np.arange(mesh.nFx)[fxp]])
        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym],
                                      np.arange(mesh.nFy)[fyp]])
        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm],
                                      np.arange(mesh.nFz)[fzp]])

        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]

        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]

        xBC = np.r_[xBC_x, xBC_y, xBC_z]
        yBC = np.r_[yBC_x, yBC_y, yBC_z]

    return xBC, yBC

Example 98

View license
def run(plotIt=True):
    """
        EM: Schenkel and Morrison Casing Model
        ======================================

        Here we create and run a FDEM forward simulation to calculate the
        vertical current inside a steel-cased. The model is based on the
        Schenkel and Morrison Casing Model, and the results are used in a 2016
        SEG abstract by Yang et al.

        .. code-block:: text

            Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on
            potential field measurements using downhole current sources:
            Geophysical prospecting, 38, 663-686.


        The model consists of:

        - Air: Conductivity 1e-8 S/m, above z = 0
        - Background: conductivity 1e-2 S/m, below z = 0
        - Casing: conductivity 1e6 S/m
            - 300m long
            - radius of 0.1m
            - thickness of 6e-3m

        Inside the casing, we take the same conductivity as the background.

        We are using an EM code to simulate DC, so we use frequency low enough
        that the skin depth inside the casing is longer than the casing length
        (f = 1e-6 Hz). The plot produced is of the current inside the casing.

        These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
        resistivity modeling of steel casing for reservoir monitoring using
        equivalent resistor network. The solver used to produce these results
        and achieve the CPU time of ~30s is Mumps, which was installed using
        pymatsolver_

        .. _pymatsolver: https://github.com/rowanc1/pymatsolver

        This example is on figshare:
        https://dx.doi.org/10.6084/m9.figshare.3126961.v1

        If you would use this example for a code comparison, or build upon it,
        a citation would be much appreciated!

    """

    if plotIt:
        import matplotlib.pylab as plt

    # ------------------ MODEL ------------------
    sigmaair = 1e-8  # air
    sigmaback = 1e-2  # background
    sigmacasing = 1e6  # casing
    sigmainside = sigmaback  # inside the casing

    casing_t = 0.006  # 1cm thickness
    casing_l = 300  # length of the casing

    casing_r = 0.1
    casing_a = casing_r - casing_t/2.  # inner radius
    casing_b = casing_r + casing_t/2.  # outer radius
    casing_z = np.r_[-casing_l, 0.]

    # ------------------ SURVEY PARAMETERS ------------------
    freqs = np.r_[1e-6]  # [1e-1, 1, 5] # frequencies
    dsz = -300  # down-hole z source location
    src_loc = np.r_[0., 0., dsz]
    inf_loc = np.r_[0., 0., 1e4]

    print('Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs])

    # ------------------ MESH ------------------
    # fine cells near well bore
    csx1, csx2 = 2e-3, 60.
    pfx1, pfx2 = 1.3, 1.3
    ncx1 = np.ceil(casing_b/csx1+2)

    # pad nicely to second cell size
    npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
    hx1a = Utils.meshTensor([(csx1, ncx1)])
    hx1b = Utils.meshTensor([(csx1, npadx1, pfx1)])
    dx1 = sum(hx1a)+sum(hx1b)
    dx1 = np.floor(dx1/csx2)
    hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)

    # second chunk of mesh
    dx2 = 300.  # uniform mesh out to here
    ncx2 = np.ceil((dx2 - dx1)/csx2)
    npadx2 = 45
    hx2a = Utils.meshTensor([(csx2, ncx2)])
    hx2b = Utils.meshTensor([(csx2, npadx2, pfx2)])
    hx = np.hstack([hx1a, hx1b, hx2a, hx2b])

    # z-direction
    csz = 0.05
    nza = 10
    # cell size, number of core cells, number of padding cells in the
    # x-direction
    ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68
    # vector of cell widths in the z-direction
    hz = Utils.meshTensor([(csz, npadzd, -1.3), (csz, ncz),
                          (csz, npadzu, 1.3)])

    # Mesh
    mesh = Mesh.CylMesh([hx, 1., hz], [0., 0., -np.sum(hz[:npadzu+ncz-nza])])

    print('Mesh Extent xmax: {0:f},: zmin: {1:f}, zmax: {2:f}'.format(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()))
    print('Number of cells', mesh.nC)

    if plotIt is True:
        fig, ax = plt.subplots(1, 1, figsize=(6, 4))
        ax.set_title('Simulation Mesh')
        mesh.plotGrid(ax=ax)
        plt.show()

    # Put the model on the mesh
    sigWholespace = sigmaback*np.ones((mesh.nC))

    sigBack = sigWholespace.copy()
    sigBack[mesh.gridCC[:, 2] > 0.] = sigmaair

    sigCasing = sigBack.copy()
    iCasingZ = ((mesh.gridCC[:, 2] <= casing_z[1]) &
                (mesh.gridCC[:, 2] >= casing_z[0]))
    iCasingX = ((mesh.gridCC[:, 0] >= casing_a) &
                (mesh.gridCC[:, 0] <= casing_b))
    iCasing = iCasingX & iCasingZ
    sigCasing[iCasing] = sigmacasing

    if plotIt is True:
        # plotting parameters
        xlim = np.r_[0., 0.2]
        zlim = np.r_[-350., 10.]
        clim_sig = np.r_[-8, 6]

        # plot models
        fig, ax = plt.subplots(1, 1, figsize=(4, 4))

        f = plt.colorbar(mesh.plotImage(np.log10(sigCasing), ax=ax)[0], ax=ax)
        ax.grid(which='both')
        ax.set_title('Log_10 (Sigma)')
        ax.set_xlim(xlim)
        ax.set_ylim(zlim)
        f.set_clim(clim_sig)

        plt.show()

    # -------------- Sources --------------------
    # Define Custom Current Sources

    # surface source
    sg_x = np.zeros(mesh.vnF[0], dtype=complex)
    sg_y = np.zeros(mesh.vnF[1], dtype=complex)
    sg_z = np.zeros(mesh.vnF[2], dtype=complex)

    nza = 2  # put the wire two cells above the surface
    ncin = 2

    # vertically directed wire
    # hook it up to casing at the surface
    sgv_indx = ((mesh.gridFz[:, 0] > casing_a) &
                (mesh.gridFz[:, 0] < casing_a + csx1))
    sgv_indz = ((mesh.gridFz[:, 2] <= +csz*nza) &
                (mesh.gridFz[:, 2] >= -csz*2))
    sgv_ind = sgv_indx & sgv_indz
    sg_z[sgv_ind] = -1.

    # horizontally directed wire
    sgh_indx = ((mesh.gridFx[:, 0] > casing_a) &
                (mesh.gridFx[:, 0] <= inf_loc[2]))
    sgh_indz = ((mesh.gridFx[:, 2] > csz*(nza-0.5)) &
                (mesh.gridFx[:, 2] < csz*(nza+0.5)))
    sgh_ind = sgh_indx & sgh_indz
    sg_x[sgh_ind] = -1.

    # hook it up to casing at the surface
    sgv2_indx = ((mesh.gridFz[:, 0] >= mesh.gridFx[sgh_ind, 0].max()) &
                 (mesh.gridFz[:, 0] <= inf_loc[2]*1.2))
    sgv2_indz = ((mesh.gridFz[:, 2] <= +csz*nza) &
                 (mesh.gridFz[:, 2] >= -csz*2))
    sgv2_ind = sgv2_indx & sgv2_indz
    sg_z[sgv2_ind] = 1.

    # assemble the source
    sg = np.hstack([sg_x, sg_y, sg_z])
    sg_p = [FDEM.Src.RawVec_e([], _, sg/mesh.area) for _ in freqs]

    # downhole source
    dg_x = np.zeros(mesh.vnF[0], dtype=complex)
    dg_y = np.zeros(mesh.vnF[1], dtype=complex)
    dg_z = np.zeros(mesh.vnF[2], dtype=complex)

    # vertically directed wire
    dgv_indx = (mesh.gridFz[:, 0] < csx1)  # go through the center of the well
    dgv_indz = ((mesh.gridFz[:, 2] <= +csz*nza) &
                (mesh.gridFz[:, 2] > dsz + csz/2.))
    dgv_ind = dgv_indx & dgv_indz
    dg_z[dgv_ind] = -1.

    # couple to the casing downhole
    dgh_indx = mesh.gridFx[:, 0] < casing_a + csx1
    dgh_indz = (mesh.gridFx[:, 2] < dsz + csz)  & (mesh.gridFx[:, 2] >= dsz)
    dgh_ind = dgh_indx & dgh_indz
    dg_x[dgh_ind] = 1.

    # horizontal part at surface
    dgh2_indx = mesh.gridFx[:, 0] <= inf_loc[2]*1.2
    dgh2_indz = sgh_indz.copy()
    dgh2_ind = dgh2_indx & dgh2_indz
    dg_x[dgh2_ind] = -1.

    # vertical part at surface
    dgv2_ind = sgv2_ind.copy()
    dg_z[dgv2_ind] = 1.

    # assemble the source
    dg = np.hstack([dg_x, dg_y, dg_z])
    dg_p = [FDEM.Src.RawVec_e([], _, dg/mesh.area) for _ in freqs]

    # ------------ Problem and Survey ---------------
    survey = FDEM.Survey(sg_p + dg_p)
    mapping = [('sigma', Maps.IdentityMap(mesh))]
    problem = FDEM.Problem3D_h(mesh, mapping=mapping, Solver=solver)
    problem.pair(survey)

    # ------------- Solve ---------------------------
    t0 = time.time()
    fieldsCasing = problem.fields(sigCasing)
    print('Time to solve 2 sources', time.time() - t0)

    # Plot current

    # current density
    jn0 = fieldsCasing[dg_p, 'j']
    jn1 = fieldsCasing[sg_p, 'j']

    # current
    in0 = [mesh.area*fieldsCasing[dg_p, 'j'][:, i] for i in range(len(freqs))]
    in1 = [mesh.area*fieldsCasing[sg_p, 'j'][:, i] for i in range(len(freqs))]

    in0 = np.vstack(in0).T
    in1 = np.vstack(in1).T

    # integrate to get z-current inside casing
    inds_inx = ((mesh.gridFz[:, 0] >= casing_a) &
                (mesh.gridFz[:, 0] <= casing_b))
    inds_inz = (mesh.gridFz[:, 2] >= dsz ) & (mesh.gridFz[:, 2] <= 0)
    inds_fz = inds_inx & inds_inz

    indsx = [False]*mesh.nFx
    inds = list(indsx) + list(inds_fz)

    in0_in = in0[np.r_[inds]]
    in1_in = in1[np.r_[inds]]
    z_in = mesh.gridFz[inds_fz, 2]

    in0_in = in0_in.reshape([in0_in.shape[0]//3, 3])
    in1_in = in1_in.reshape([in1_in.shape[0]//3, 3])
    z_in = z_in.reshape([z_in.shape[0]//3, 3])

    I0 = in0_in.sum(1).real
    I1 = in1_in.sum(1).real
    z_in = z_in[:, 0]

    if plotIt is True:
        fig, ax = plt.subplots(1, 2, figsize=(12, 4))

        ax[0].plot(z_in, np.absolute(I0), z_in, np.absolute(I1))
        ax[0].legend(['top casing', 'bottom casing'], loc='best')
        ax[0].set_title('Magnitude of Vertical Current in Casing')

        ax[1].semilogy(z_in, np.absolute(I0), z_in, np.absolute(I1))
        ax[1].legend(['top casing', 'bottom casing'], loc='best')
        ax[1].set_title('Magnitude of Vertical Current in Casing')
        ax[1].set_ylim([1e-2, 1.])

        plt.show()

Example 99

Project: simpeg
Source File: DiffOperators.py
View license
    def getBCProjWF(self, BC, discretization='CC'):
        """

        The weak form boundary condition projection matrices.

        Examples::
            # Neumann in all directions
            BC = 'neumann'

            # 3D, Dirichlet in y Neumann else
            BC = ['neumann', 'dirichlet', 'neumann']

            # 3D, Neumann in x on bottom of domain, Dirichlet else
            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
        """

        if discretization is not 'CC':
            raise NotImplementedError('Boundary conditions only implemented'
                                      'for CC discretization.')

        if(type(BC) is str):
            BC = [BC for _ in self.vnC]  # Repeat the str self.dim times
        elif(type(BC) is list):
            assert len(BC) == self.dim, 'BC list must be the size of your mesh'
        else:
            raise Exception("BC must be a str or a list.")

        for i, bc_i in enumerate(BC):
            BC[i] = checkBC(bc_i)

        def projDirichlet(n, bc):
            bc = checkBC(bc)
            ij = ([0, n], [0, 1])
            vals = [0, 0]
            if(bc[0] == 'dirichlet'):
                vals[0] = -1
            if(bc[1] == 'dirichlet'):
                vals[1] = 1
            return sp.csr_matrix((vals, ij), shape=(n+1, 2))

        def projNeumannIn(n, bc):
            bc = checkBC(bc)
            P = sp.identity(n+1).tocsr()
            if(bc[0] == 'neumann'):
                P = P[1:, :]
            if(bc[1] == 'neumann'):
                P = P[:-1, :]
            return P

        def projNeumannOut(n, bc):
            bc = checkBC(bc)
            ij   = ([0, 1], [0, n])
            vals = [0,0]
            if(bc[0] == 'neumann'):
                vals[0] = 1
            if(bc[1] == 'neumann'):
                vals[1] = 1
            return sp.csr_matrix((vals, ij), shape=(2, n+1))

        n = self.vnC
        indF = self.faceBoundaryInd
        if(self.dim == 1):
            Pbc = projDirichlet(n[0], BC[0])
            indF = indF[0] | indF[1]
            Pbc = Pbc*sdiag(self.area[indF])

            Pin = projNeumannIn(n[0], BC[0])

            Pout = projNeumannOut(n[0], BC[0])

        elif(self.dim == 2):
            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
            Pbc = Pbc*sdiag(self.area[indF])

            P1 = sp.kron(speye(n[1]), projNeumannIn(n[0], BC[0]))
            P2 = sp.kron(projNeumannIn(n[1], BC[1]), speye(n[0]))
            Pin = sp.block_diag((P1, P2), format="csr")

            P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
            P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
            Pout = sp.block_diag((P1, P2), format="csr")

        elif(self.dim == 3):
            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] |
                          indF[5])]
            Pbc = Pbc*sdiag(self.area[indF])

            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
            P2 = kron3(speye(n[2]), projNeumannIn(n[1], BC[1]), speye(n[0]))
            P3 = kron3(projNeumannIn(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pin = sp.block_diag((P1, P2, P3), format="csr")

            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannOut(n[0], BC[0]))
            P2 = kron3(speye(n[2]), projNeumannOut(n[1], BC[1]), speye(n[0]))
            P3 = kron3(projNeumannOut(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pout = sp.block_diag((P1, P2, P3), format="csr")

        return Pbc, Pin, Pout

Example 100

Project: simpeg
Source File: DiffOperators.py
View license
    def getBCProjWF_simple(self, discretization='CC'):
        """
        The weak form boundary condition projection matrices
        when mixed boundary condition is used
        """

        if discretization is not 'CC':
            raise NotImplementedError('Boundary conditions only implemented'
                                      'for CC discretization.')

        def projBC(n):
            ij = ([0, n], [0, 1])
            vals = [0, 0]
            vals[0] = 1
            vals[1] = 1
            return sp.csr_matrix((vals, ij), shape=(n+1, 2))

        def projDirichlet(n, bc):
            bc = checkBC(bc)
            ij = ([0, n], [0, 1])
            vals = [0, 0]
            if(bc[0] == 'dirichlet'):
                vals[0] = -1
            if(bc[1] == 'dirichlet'):
                vals[1] = 1
            return sp.csr_matrix((vals, ij), shape=(n+1, 2))

        BC = [['dirichlet', 'dirichlet'], ['dirichlet', 'dirichlet'],
              ['dirichlet', 'dirichlet']]
        n = self.vnC
        indF = self.faceBoundaryInd

        if(self.dim == 1):
            Pbc = projDirichlet(n[0], BC[0])
            B = projBC(n[0])
            indF = indF[0] | indF[1]
            Pbc = Pbc*sdiag(self.area[indF])

        elif(self.dim == 2):
            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
            B1 = sp.kron(speye(n[1]), projBC(n[0]))
            B2 = sp.kron(projBC(n[1]), speye(n[0]))
            B = sp.block_diag((B1, B2), format="csr")
            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
            Pbc = Pbc*sdiag(self.area[indF])

        elif(self.dim == 3):
            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
            B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0]))
            B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0]))
            B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0]))
            B = sp.block_diag((B1, B2, B3), format="csr")
            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
            Pbc = Pbc*sdiag(self.area[indF])

        return Pbc, B.T