Python Numeric.dot函数代码示例

 def rtz(self, pt):
     d = pt - self.p0
     z = dot(d, self.zn)
     d = d - z * self.zn
     r = vlen(d)
     theta = Numeric.arctan2(dot(d, self.v), dot(d, self.u))
     return Numeric.array((r, theta, z), "d")

 def __call__(self,input,error=False):
     X = self.scale_input(input)
     if not error:
         Y = dot(self.w,X)
     else:
         Y = dot(self.w,X), zeros(self.num_outputs)
     return self.scale_output(Y)

 def __init__(self, point1 = None, point2 = None, slab = None):
     # Huaicai 4/23/05: added some comments below to help understand the code.
     if slab:
         # convert from 2d (x, y) coordinates into its 3d world (x, y, 0) 
         #coordinates(the lower-left and upper-right corner). 
         #In another word, the 3d coordinates minus the z offset of the plane
         x = dot(A(point1), A(point2))
         # Get the vector from upper-right point to the lower-left point
         dx = subtract.reduce(x)
         # Get the upper-left and lower right corner points
         oc = x[1] + V(point2[0]*dot(dx,point2[0]), 
                       point2[1]*dot(dx,point2[1]))
         # Get the four 3d cooridinates on the bottom cookie-cutting plane
         sq1 = cat(x,oc) + slab.normal*dot(slab.point, slab.normal)
         # transfer the above 4 3d coordinates in parallel to get that on
         #the top plane, put them together
         sq1 = cat(sq1, sq1+slab.thickness*slab.normal)
         self.data = V(maximum.reduce(sq1), minimum.reduce(sq1))
     elif point2:
         # just 2 3d points
         self.data = V(maximum(point1, point2), minimum(point1, point2))
     elif point1:
         # list of points: could be 2d or 3d?  +/- 1.8 to make the bounding 
         #box enclose the vDw ball of an atom?
         self.data = V(maximum.reduce(point1) + BBOX_MARGIN, 
                       minimum.reduce(point1) - BBOX_MARGIN)
     else:
         # a null bbox
         self.data = None

    def _leftDown_preparation_for_dragging(self, objectUnderMouse, event):
        """
        Handle left down event. Preparation for rotation and/or selection
        This method is called inside of self.leftDown.
        @param event: The mouse left down event.
        @type  event: QMouseEvent instance
        @see: self.leftDown
        @see: self.leftDragRotation
        Overrides _superclass._leftDown_preparation_for_dragging
        """

        _superclass._leftDown_preparation_for_dragging(self, objectUnderMouse, event)
        self.o.SaveMouse(event)
        self.picking = True
        self.dragdist = 0.0
        # delta for constrained rotations.
        self.rotDelta = 0

        if self.rotateOption == "ROTATEDEFAULT":
            self.o.trackball.start(self.o.MousePos[0], self.o.MousePos[1])
        else:
            if self.rotateOption == "ROTATEX":
                ma = V(1, 0, 0)  # X Axis
                self.axis = "X"
            elif self.rotateOption == "ROTATEY":
                ma = V(0, 1, 0)  # Y Axis
                self.axis = "Y"
            elif self.rotateOption == "ROTATEZ":
                ma = V(0, 0, 1)  # Z Axis
                self.axis = "Z"
            elif self.rotateOption == "ROT_TRANS_ALONG_AXIS":
                # The method 'self._leftDown_preparation_for_dragging should
                # never be reached if self.rotateOption is 'ROT_TRANS_ALONG_AXIS'
                # If this code is reached, it indicates a bug. So fail gracefully
                # by calling self.leftADown()
                if debug_flags.atom_debug:
                    print_compact_stack(
                        "bug: _leftDown_preparation_for_dragging" " called for rotate option" "'ROT_TRANS_ALONG_AXIS'"
                    )

                self.leftADown(objectUnderMouse, event)
                return

            else:
                print "Move_Command: Error - unknown rotateOption value =", self.rotateOption
                return

            ma = norm(V(dot(ma, self.o.right), dot(ma, self.o.up)))
            # When in the front view, right = 1,0,0 and up = 0,1,0, so ma will
            # be computed as 0,0.This creates a special case problem when the
            # user wants to constrain rotation around the Z axis because Zmat
            # will be zero.  So we have to test for this case (ma = 0,0) and
            # fix ma to -1,0.  This was needed to fix bug 537.  Mark 050420
            if ma[0] == 0.0 and ma[1] == 0.0:
                ma = [-1.0, 0.0]
            self.Zmat = A([ma, [-ma[1], ma[0]]])

        self.leftDownType = "ROTATE"

        return

def drawLinearDimension(color,      # what color are we drawing this in
                        right, up,  # screen directions mapped to xyz coords
                        bpos,       # position of the handle for moving the text
                        p0, p1,     # positions of the ends of the dimension
                        text, highlighted=False):
    outOfScreen = cross(right, up)
    bdiff = bpos - 0.5 * (p0 + p1)
    csys = CylindricalCoordinates(p0, p1 - p0, bdiff, right)
    # This works OK until we want to keep the text right side up, then the
    # criterion for right-side-up-ness changes because we've changed 'up'.
    br, bt, bz = csys.rtz(bpos)
    e0 = csys.xyz((br + 0.5, 0, 0))
    e1 = csys.xyz((br + 0.5, 0, 1))
    drawline(color, p0, e0)
    drawline(color, p1, e1)
    v0 = csys.xyz((br, 0, 0))
    v1 = csys.xyz((br, 0, 1))
    if highlighted:
        drawline(color, v0, v1, width=THICKLINEWIDTH)
    else:
        drawline(color, v0, v1)
    # draw arrowheads at the ends
    a1, a2 = 0.25, 1.0 * csys.zinv
    arrow00 = csys.xyz((br + a1, 0, a2))
    arrow01 = csys.xyz((br - a1, 0, a2))
    drawline(color, v0, arrow00)
    drawline(color, v0, arrow01)
    arrow10 = csys.xyz((br + a1, 0, 1-a2))
    arrow11 = csys.xyz((br - a1, 0, 1-a2))
    drawline(color, v1, arrow10)
    drawline(color, v1, arrow11)
    # draw the text for the numerical measurement, make
    # sure it goes from left to right
    xflip = dot(csys.z, right) < 0
    # then make sure it's right side up
    theoreticalRight = (xflip and -csys.z) or csys.z
    theoreticalOutOfScreen = cross(theoreticalRight, bdiff)
    yflip = dot(theoreticalOutOfScreen, outOfScreen) < 0
    if debug_flags.atom_debug:
        print "DEBUG INFO FROM drawLinearDimension"
        print csys
        print theoreticalRight, theoreticalOutOfScreen
        print xflip, yflip
    if yflip:
        def fx(y):
            return br + 1.5 - y / (1. * HEIGHT)
    else:
        def fx(y):
            return br + 0.5 + y / (1. * HEIGHT)
    if xflip:
        def fz(x):
            return 0.9 - csys.zinv * x / (1. * WIDTH)
    else:
        def fz(x):
            return 0.1 + csys.zinv * x / (1. * WIDTH)
    def tfm(x, y, fx=fx, fz=fz):
        return csys.xyz((fx(y), 0, fz(x)))
    f3d = Font3D()
    f3d.drawString(text, tfm=tfm, color=color)

 def Q(self,state,action=None):
     """
     Compute Q(s,a) from W*s.
     """
     if action is None:
         return dot(self.w, state)
     else:
         return dot(self.w[action],state)

  def go(self, xn):
    if (self.x == None):
      self.x = ones(len(self.b)) * xn * 1.0
      self.y = ones(len(self.a)) * xn * 1.0 * sum(self.b) / (1+sum(self.a))

    self.x = concatenate([[xn], self.x[:-1]])
    yn = dot(self.b, self.x) - dot(self.a, self.y)
    self.y = concatenate([[yn], self.y[:-1]])
    return yn

 def constrainedPosition(self):
     a = self.atoms
     p0, p1, p2, p3 = a[0].posn(), a[1].posn(), a[2].posn(), a[3].posn()
     axis = norm(p2 - p1)
     midpoint = 0.5 * (p1 + p2)
     return _constrainHandleToAngle(self.center() + self.handle_offset,
                                    p0 - dot(p0 - midpoint, axis) * axis,
                                    midpoint,
                                    p3 - dot(p3 - midpoint, axis) * axis)

    def rotateAboutPoint(self):
        """
        Rotates the selected entities along the specified vector, about the
        specified pivot point (pivot point it the starting point of the
        drawn vector.
        """

        if len(self.mouseClickPoints) != self.mouseClickLimit:
            print_compact_stack("Rotate about point bug: mouseclick points != mouseclicklimit: ")
            return


        pivotPoint = self.mouseClickPoints[0]
        ref_vec_endPoint = self.mouseClickPoints[1]
        rot_vec_endPoint = self.mouseClickPoints[2]

        reference_vec = norm(ref_vec_endPoint - pivotPoint)

        lineVector = norm(rot_vec_endPoint - pivotPoint)


        #lineVector = endPoint - startPoint

        quat1 = Q(lineVector, reference_vec)

        #DEBUG Disabled temporarily . will not be used
        if dot(lineVector, reference_vec) < 0:
            theta = math.pi - quat1.angle
        else:
            theta = quat1.angle

        #TEST_DEBUG-- Works fine
        theta = quat1.angle

        rot_axis = cross(lineVector, reference_vec)


        if dot(lineVector, reference_vec) < 0:
            rot_axis = - rot_axis

        cross_prod_1 = norm(cross(reference_vec, rot_axis))
        cross_prod_2 = norm(cross(lineVector, rot_axis))

        if dot(cross_prod_1, cross_prod_2) < 0:
            quat2 = Q(rot_axis, theta)
        else:
            quat2 = Q(rot_axis, - theta)

        movables = self.graphicsMode.getMovablesForLeftDragging()
        self.assy.rotateSpecifiedMovables(
            quat2,
            movables = movables,
            commonCenter = pivotPoint)

        self.glpane.gl_update()
        return

def calc_torsion_angle(atom_list):
    """
    Calculates torsional angle defined by four atoms, A1-A2-A3-A4,
    Return torsional angle value between atoms A2 and A3.

    @param atom_list: list of four atoms describing the torsion bond
    @type atom_list: list

    @return: value of the torsional angle (float)
    """
    # Note: this appears to be very general and perhaps ought to be moved to a more
    # general place (someday), perhaps VQT.py or nearby. [bruce 080828 comment]

    from Numeric import dot
    from math import atan2, pi, sqrt
    from geometry.VQT import cross

    if len(atom_list) != 4:
        # The list has to have four members.
        return 0.0

    # Calculate pairwise distances
    v12 = atom_list[0].posn() - atom_list[1].posn()
    v43 = atom_list[3].posn() - atom_list[2].posn()
    v23 = atom_list[1].posn() - atom_list[2].posn()

    # p is a vector perpendicular to the plane defined by atoms 1,2,3
    # p is perpendicular to v23_v12 plane
    p = cross(v23, v12)

    # x is a vector perpendicular to the plane defined by atoms 2,3,4.
    # x is perpendicular to v23_v43 plane
    x = cross(v23, v43)

    # y is perpendicular to v23_x plane
    y = cross(v23, x)

    # Calculate lengths of the x, y vectors.
    u1 = dot(x, x)
    v1 = dot(y, y)

    if u1 < 0.0 or \
       v1 < 0.0:
        return 360.0

    u2 = dot(p, x) / sqrt(u1)
    v2 = dot(p, y) / sqrt(v1)

    if u2 != 0.0 and \
       v2 != 0.0:
        # calculate the angle
        return atan2(v2, u2) * (180.0 / pi)
    else:
        return 360.0

    def rotateAboutPoint(self):
        """
        Rotates the selected entities along the specified vector, about the
        specified pivot point (pivot point it the starting point of the
        drawn vector.
        """
        startPoint = self.mouseClickPoints[0]
        endPoint = self.mouseClickPoints[1]
        pivotAtom = self.graphicsMode.pivotAtom
        #initial assignment of reference_vec. The selected movables will be
        #rotated by the angle between this vector and the lineVector
        reference_vec = self.glpane.right
        if isinstance(pivotAtom, Atom) and not pivotAtom.molecule.isNullChunk():
            mol = pivotAtom.molecule
            reference_vec, node_junk = mol.getAxis_of_self_or_eligible_parent_node(
                atomAtVectorOrigin = pivotAtom)
            del node_junk
        else:
            reference_vec = self.glpane.right

        lineVector = endPoint - startPoint

        quat1 = Q(lineVector, reference_vec)

        #DEBUG Disabled temporarily . will not be used
        ##if dot(lineVector, reference_vec) < 0:
            ##theta = math.pi - quat1.angle
        ##else:
            ##theta = quat1.angle

        #TEST_DEBUG-- Works fine
        theta = quat1.angle

        rot_axis = cross(lineVector, reference_vec)

        if dot(lineVector, reference_vec) < 0:
            rot_axis = - rot_axis

        cross_prod_1 = norm(cross(reference_vec, rot_axis))
        cross_prod_2 = norm(cross(lineVector, rot_axis))

        if dot(cross_prod_1, cross_prod_2) < 0:
            quat2 = Q(rot_axis,  theta)
        else:
            quat2 = Q(rot_axis,  - theta)


        movables = self.graphicsMode.getMovablesForLeftDragging()
        self.assy.rotateSpecifiedMovables(
            quat2,
            movables = movables,
            commonCenter = startPoint)

        self.glpane.gl_update()

 def project_2d_noeyeball(self, pt):
     """
     Bruce: Project a point into our plane (ignoring eyeball). Warning: arbitrary origin!
        
     Huaicai 4/20/05: This is just to project pt into a 2d coordinate 
     system (self.right, self.up) on a plane through pt and parallel to the screen 
     plane. For perspective projection, (x, y) on this plane is different than that on the plane 
     through pov.
     """
     x, y = self.right, self.up
     return V(dot(pt, x), dot(pt, y))

  def kfupdate(self, dt, rs):
    self.ab  = array((rs.ax, -rs.ay, -rs.az))

    ph = self.phi
    th = self.theta
    P  = self.pmat

    A  = array(((-rs.q*cos(ph)*tan(th)+rs.r*sin(ph)*tan(th), (-rs.q*sin(ph)+rs.r*cos(ph))/(cos(th)*cos(th))),
                (rs.q*sin(ph)-rs.r*cos(ph)                 , 0)))

    dph = rs.p - rs.q*sin(ph)*tan(th) - rs.r*cos(ph)*tan(th)
    dth =      - rs.q*cos(ph)         - rs.r*sin(ph)
    dP  = dot(A, P) + dot(P, transpose(A)) + self.Q

    ph = ph + dph * dt
    th = th + dth * dt
    P  = P  + dP  * dt

    Cx = array((0               , cos(th)))
    Cy = array((-cos(th)*cos(ph), sin(th)*sin(ph)))
    Cz = array((cos(th)*sin(ph) , sin(th)*cos(ph)))
    C  = array((Cx, Cy, Cz))

    L = dot(dot(P, transpose(C)), inverse(self.R + dot(dot(C, P), transpose(C))))
    h = array((sin(th), -cos(th)*sin(ph), -cos(th)*cos(ph)))

    P  = dot(identity(2) - dot(L, C), P)
    ph = ph + dot(L[0], self.ab - h)
    th = th + dot(L[1], self.ab - h) 

    ph = ((ph+pi) % (2*pi)) - pi;
    th = ((th+pi) % (2*pi)) - pi;

    self.pmat  = P
    self.phi   = ph 
    self.theta = th

    psidot = rs.q * sin(ph) / cos(th) + rs.r * cos(ph) / cos(th);
    self.psi += psidot * dt;

    self.quat = eul2quat(ph,th,self.psi)
    # self.dcmb2e = quat2dcm(quatinv(self.quat))

    # self.ae = dot(self.dcmb2e, self.ab)
    self.ae = quatrotate(quatinv(self.quat), self.ab)
    self.ae[2] = -self.ae[2]-1
    self.ae[1] = -self.ae[1]

    self.ve += self.ae * dt * 9.81
    self.xe += self.ve * dt

    def setup_quat_center(self, atomList = None):
        """
        Setup the plane's quat using a list of atoms.

        If no atom list is supplied, the plane is centered in the glpane
        and parallel to the screen.

        @param atomList: A list of atoms.
        @type  atomList: list

        """
        if atomList:    
            self.atomPos = []
            for a in atomList:
                self.atomPos += [a.posn()]    
            planeNorm = self._getPlaneOrientation(self.atomPos)
            if dot(planeNorm, self.glpane.lineOfSight) < 0:
                planeNorm = -planeNorm                
            self.center = add.reduce(self.atomPos) / len(self.atomPos)
            self.quat   = Q(V(0.0, 0.0, 1.0), planeNorm) 
        else:       
            self.center = V(0.0, 0.0, 0.0)
            # Following makes sure that Plane edges are parallel to
            # the 3D workspace borders. Fixes bug 2448
            x, y ,z = self.glpane.right, self.glpane.up, self.glpane.out
            self.quat  = Q(x, y, z) 
            self.quat += Q(self.glpane.right, pi)

def selection_polyhedron(basepos, borderwidth = 1.8):
    """
    Given basepos (a Numeric array of 3d positions), compute and return (as a list of vertices, each pair being an edge to draw)
    a simple bounding polyhedron, convenient for designating the approximate extent of this set of points.
    (This is used on the array of atom and open bond positions in a chunk to designate a selected chunk.)
       Make every face farther out than it needs to be to enclose all points in basepos, by borderwidth (default 1.8).
    Negative values of borderwidth might sometimes work, but are likely to fail drastically if the polyhedron is too small.
       The orientation of the faces matches the coordinate axes used in basepos (in typical calls, these are chunk-relative).
    [#e Someday we might permit caller-specified orientation axes. The axes of inertia would be interesting to try.]
    """
    #bruce 060226/060605 made borderwidth an argument (non-default values untested); documented retval format
    #bruce 060119 split this out of shakedown_poly_evals_evecs_axis() in chunk.py.
    # Note, it has always had some bad bugs for certain cases, like long diagonal rods.

    if not len(basepos):
        return [] # a guess
    
    # find extrema in many directions
    xtab = dot(basepos, polyXmat)
    mins = minimum.reduce(xtab) - borderwidth
    maxs = maximum.reduce(xtab) + borderwidth

    polyhedron = makePolyList(cat(maxs,mins))
        # apparently polyhedron is just a long list of vertices [v0,v1,v2,v3,...] to be passed to GL_LINES
        # (via drawlinelist), which will divide them into pairs and draw lines v0-v1, v2-v3, etc.
        # Maybe we should generalize the format, so some callers could also draw faces.
        # [bruce 060226/060605 comment]
    return polyhedron