original_docstring: 
          Constructs a projector onto the non-gauge parameter space.
  
          This is useful for isolating the gauge degrees of freedom from the non-gauge
          degrees of freedom.
  
          Parameters
          ----------
          item_weights : dict, optional
              Dictionary of weighting factors for individual gates and spam operators.
              Keys can be gate, state preparation, POVM effect, spam labels, or the
              special strings "gates" or "spam" whic represent the entire set of gate
              or SPAM operators, respectively.  Values are floating point numbers.
              These weights define the metric used to compute the non-gauge space,
              *orthogonal* the gauge space, that is projected onto.
  
          non_gauge_mix_mx : numpy array, optional
              An array of shape (n_non_gauge_params,n_gauge_params) specifying how to
              mix the non-gauge degrees of freedom into the gauge degrees of
              freedom that are projected out by the returned object.  This argument
              essentially sets the off-diagonal block of the metric used for
              orthogonality in the "gauge + non-gauge" space.  It is for advanced
              usage and typically left as None (the default).
  
          Returns
          -------
          numpy array
              The projection operator as a N x N matrix, where N is the number
              of parameters (obtained via num_params()).  This projector acts on
              parameter-space, and has rank equal to the number of non-gauge
              degrees of freedom.
          
docstring: Constructs a projector onto the non-gauge parameter space.
  This is useful for isolating the gauge degrees of freedom from the non-gauge
  degrees of freedom.
  
  Parameters
  
  item_weights : dict, optional
  Dictionary of weighting factors for individual gates and spam operators.
  Keys can be gate, state preparation, POVM effect, spam labels, or the
  special strings "gates" or "spam" whic represent the entire set of gate
  or SPAM operators, respectively.  Values are floating point numbers.
  These weights define the metric used to compute the non-gauge space,
  orthogonal* the gauge space, that is projected onto.
  
  non_gauge_mix_mx : numpy array, optional
  An array of shape (n_non_gauge_params,n_gauge_params) specifying how to
  mix the non-gauge degrees of freedom into the gauge degrees of
  freedom that are projected out by the returned object.  This argument
  essentially sets the off-diagonal block of the metric used for
  orthogonality in the "gauge + non-gauge" space.  It is for advanced
  usage and typically left as None (the default).
  
  Returns
  
  numpy array
  The projection operator as a N x N matrix, where N is the number
  of parameters (obtained via num_params()).  This projector acts on
  parameter-space, and has rank equal to the number of non-gauge
  degrees of freedom.
docstring_tokens: Constructs a projector onto the non - gauge parameter space . This is useful for isolating the gauge degrees of freedom from the non - gauge degrees of freedom . Parameters item_weights : dict optional Dictionary of weighting factors for individual gates and spam operators . Keys can be gate state preparation POVM effect spam labels or the special strings " gates " or " spam " whic represent the entire set of gate or SPAM operators respectively . Values are floating point numbers . These weights define the metric used to compute the non - gauge space orthogonal * the gauge space that is projected onto . non_gauge_mix_mx : numpy array optional An array of shape ( n_non_gauge_params n_gauge_params ) specifying how to mix the non - gauge degrees of freedom into the gauge degrees of freedom that are projected out by the returned object . This argument essentially sets the off - diagonal block of the metric used for orthogonality in the " gauge + non - gauge " space . It is for advanced usage and typically left as None ( the default ) . Returns numpy array The projection operator as a N x N matrix where N is the number of parameters ( obtained via num_params () ) . This projector acts on parameter - space and has rank equal to the number of non - gauge degrees of freedom .
code: def nongauge_projector(self, item_weights=None, non_gauge_mix_mx=None):
          """
          Constructs a projector onto the non-gauge parameter space.
  
          This is useful for isolating the gauge degrees of freedom from the non-gauge
          degrees of freedom.
  
          Parameters
          ----------
          item_weights : dict, optional
              Dictionary of weighting factors for individual gates and spam operators.
              Keys can be gate, state preparation, POVM effect, spam labels, or the
              special strings "gates" or "spam" whic represent the entire set of gate
              or SPAM operators, respectively.  Values are floating point numbers.
              These weights define the metric used to compute the non-gauge space,
              *orthogonal* the gauge space, that is projected onto.
  
          non_gauge_mix_mx : numpy array, optional
              An array of shape (n_non_gauge_params,n_gauge_params) specifying how to
              mix the non-gauge degrees of freedom into the gauge degrees of
              freedom that are projected out by the returned object.  This argument
              essentially sets the off-diagonal block of the metric used for
              orthogonality in the "gauge + non-gauge" space.  It is for advanced
              usage and typically left as None (the default).
  
          Returns
          -------
          numpy array
              The projection operator as a N x N matrix, where N is the number
              of parameters (obtained via num_params()).  This projector acts on
              parameter-space, and has rank equal to the number of non-gauge
              degrees of freedom.
          """
  
          # We want to divide the Model-space H (a Hilbert space, 56-dimensional in the 1Q, 3-gate, 2-vec case)
          # into the direct sum of gauge and non-gauge spaces, and find projectors onto each
          # sub-space (the non-gauge space in particular).
          #
          # Within the Model H-space lies a gauge-manifold of maximum chi2 (16-dimensional in 1Q case),
          #  where the gauge-invariant chi2 function is constant.  At a single point (Model) P on this manifold,
          #  chosen by (somewhat arbitrarily) fixing the gauge, we can talk about the tangent space
          #  at that point.  This tangent space is spanned by some basis (~16 elements for the 1Q case),
          #  which associate with the infinitesimal gauge transform ?generators? on the full space.
          #  The subspace of H spanned by the derivatives w.r.t gauge transformations at P (a Model) spans
          #  the gauge space, and the complement of this (in H) is the non-gauge space.
          #
          #  An element of the gauge group can be written gg = exp(-K), where K is a n x n matrix.  If K is
          #   hermitian then gg is unitary, but this need not be the case.  A gauge transform acts on a
          #   gatset via Model => gg^-1 G gg, i.e. G => exp(K) G exp(-K).  We care about the action of
          #   infinitesimal gauge tranformations (b/c the *derivative* vectors span the tangent space),
          #   which act as:
          #    G => (I+K) G (I-K) = G + [K,G] + ignored(K^2), where [K,G] := KG-GK
          #
          # To produce a projector onto the gauge-space, we compute the *column* vectors
          #  dG_ij = [K_ij,G], where K_ij is the i,j-th matrix unit (56x1 in the 1Q case, 16 such column vectors)
          #  and then form a projector in the standard way.
          #  (to project onto |v1>, |v2>, etc., form P = sum_i |v_i><v_i|)
          #
          #Typically nOpParams < len(dG_ij) and linear system is overconstrained
          #   and no solution is expected.  If no solution exists, simply ignore
          #
          # So we form P = sum_ij dG_ij * transpose(dG_ij) (56x56 in 1Q case)
          #              = dG * transpose(dG)              where dG_ij form the *columns* of dG (56x16 in 1Q case)
          # But dG_ij are not orthonormal, so really we need a slight modification,
          #  otherwise P'P' != P' as must be the case for a projector:
          #
          # P' = dG * (transpose(dG) * dG)^-1 * transpose(dG) (see wikipedia on projectors)
          #
          #    or equivalently (I think)
          #
          # P' = pseudo-inv(P)*P
          #
          #  since the pseudo-inv is defined by P*pseudo-inv(P) = I and so P'P' == P'
          #  and P' is our gauge-projector!
  
          # Note: In the general case of parameterized gates (i.e. non-fully parameterized gates), there are fewer
          #   gate parameters than the size of dG_ij ( < 56 in 1Q case).  In this general case, we want to know
          #   what (if any) change in gate parameters produces the change dG_ij of the operation matrix elements.
          #   That is, we solve dG_ij = derivWRTParams * dParams_ij for dParams_ij, where derivWRTParams is
          #   the derivative of the operation matrix elements with respect to the gate parameters (derivWRTParams
          #   is 56x(nOpParams) and dParams_ij is (nOpParams)x1 in the 1Q case) -- then substitute dG_ij
          #   with dParams_ij above to arrive at a (nOpParams)x(nOpParams) projector (compatible with
          #   hessian computed by model).
          #
          #   Actually, we want to know if what changes gate parameters
          #   produce changes in the span of all the dG_ij, so really we want the intersection of the space
          #   defined by the columns of derivWRTParams (the "gate-parameter range" space) and the dG_ij vectors.
          #
          #   This intersection is determined by nullspace( derivWRTParams | dG ), where the pipe denotes
          #     concatenating the matrices together.  Then, if x is a basis vector of the nullspace
          #     x[0:nOpParams] is the basis vector for the intersection space within the gate parameter space,
          #     that is, the analogue of dParams_ij in the single-dG_ij introduction above.
          #
          #   Still, we just substitue these dParams_ij vectors (as many as the nullspace dimension) for dG_ij
          #   above to get the general case projector.
  
          nParams = self.Np
          dPG = self._buildup_dpg()
  
          #print("DB: shapes = ",dP.shape,dG.shape)
          nullsp = _mt.nullspace_qr(dPG)  # columns are nullspace basis vectors
          gen_dG = nullsp[0:nParams, :]  # take upper (gate-param-segment) of vectors for basis
          # of subspace intersection in gate-parameter space
          #Note: gen_dG == "generalized dG", and is (nParams)x(nullSpaceDim==gaugeSpaceDim), so P
          #  (below) is (nParams)x(nParams) as desired.
  
          #DEBUG
          #nElements = self.num_elements
          #for iRow in range(nElements):
          #    pNorm = _np.linalg.norm(dP[iRow])
          #    if pNorm < 1e-6:
          #        print "Row %d of dP is zero" % iRow
          #
          #print "------------------------------"
          #print "nParams = ",nParams
          #print "nElements = ",nElements
          #print " shape(M) = ",M.shape
          #print " shape(dG) = ",dG.shape
          #print " shape(dP) = ",dP.shape
          #print " shape(gen_dG) = ",gen_dG.shape
          #print " shape(nullsp) = ",nullsp.shape
          #print " rank(dG) = ",_np.linalg.matrix_rank(dG)
          #print " rank(dP) = ",_np.linalg.matrix_rank(dP)
          #print "------------------------------"
          #assert(_np.linalg.norm(_np.imag(gen_dG)) < 1e-9) #gen_dG is real
  
          if non_gauge_mix_mx is not None:
              msg = "You've set both non_gauge_mix_mx and item_weights, both of which"\
                  + " set the gauge metric... You probably don't want to do this."
              assert(item_weights is None), msg
  
              # BEGIN GAUGE MIX ----------------------------------------
              # nullspace of gen_dG^T (mx with gauge direction vecs as rows) gives non-gauge directions
              # columns are the non-gauge directions *orthogonal* to the gauge directions
              nonGaugeDirections = _mt.nullspace_qr(gen_dG.T)
  
              #for each column of gen_dG, which is a gauge direction in model parameter space,
              # we add some amount of non-gauge direction, given by coefficients of the
              # numNonGaugeParams non-gauge directions.
              gen_dG = gen_dG + _np.dot(nonGaugeDirections, non_gauge_mix_mx)  # add non-gauge direction components in
              # dims: (nParams,n_gauge_params) + (nParams,n_non_gauge_params) * (n_non_gauge_params,n_gauge_params)
              # non_gauge_mix_mx is a (n_non_gauge_params,n_gauge_params) matrix whose i-th column specifies the
              #  coefficents to multipy each of the non-gauge directions by before adding them to the i-th
              #  direction to project out (i.e. what were the pure gauge directions).
  
              #DEBUG
              #print "gen_dG shape = ",gen_dG.shape
              #print "NGD shape = ",nonGaugeDirections.shape
              #print "NGD rank = ",_np.linalg.matrix_rank(nonGaugeDirections, P_RANK_TOL)
              #END GAUGE MIX ----------------------------------------
  
          # Build final non-gauge projector by getting a mx of column vectors
          # orthogonal to the cols fo gen_dG:
          #     gen_dG^T * gen_ndG = 0 => nullspace(gen_dG^T)
          # or for a general metric W:
          #     gen_dG^T * W * gen_ndG = 0 => nullspace(gen_dG^T * W)
          # (This is instead of construction as I - gaugeSpaceProjector)
  
          if item_weights is not None:
              metric_diag = _np.ones(self.Np, 'd')
              opWeight = item_weights.get('gates', 1.0)
              spamWeight = item_weights.get('spam', 1.0)
              for lbl, gate in self.operations.items():
                  metric_diag[gate.gpindices] = item_weights.get(lbl, opWeight)
              for lbl, vec in _itertools.chain(iter(self.preps.items()),
                                               iter(self.effects.items())):
                  metric_diag[vec.gpindices] = item_weights.get(lbl, spamWeight)
              metric = _np.diag(metric_diag)
              #OLD: gen_ndG = _mt.nullspace(_np.dot(gen_dG.T,metric))
              gen_ndG = _mt.nullspace_qr(_np.dot(gen_dG.T, metric))
          else:
              #OLD: gen_ndG = _mt.nullspace(gen_dG.T) #cols are non-gauge directions
              gen_ndG = _mt.nullspace_qr(gen_dG.T)  # cols are non-gauge directions
              # print("DB: nullspace of gen_dG (shape = %s, rank=%d) = %s" \
              #       % (str(gen_dG.shape),_np.linalg.matrix_rank(gen_dG),str(gen_ndG.shape)))
  
          # ORIG WAY: use psuedo-inverse to normalize projector.  Ran into problems where
          #  default rcond == 1e-15 didn't work for 2-qubit case, but still more stable than inv method below
          P = _np.dot(gen_ndG, _np.transpose(gen_ndG))  # almost a projector, but cols of dG are not orthonormal
          Pp = _np.dot(_np.linalg.pinv(P, rcond=1e-7), P)  # make P into a true projector (onto gauge space)
  
          # ALT WAY: use inverse of dG^T*dG to normalize projector (see wikipedia on projectors, dG => A)
          #  This *should* give the same thing as above, but numerical differences indicate the pinv method
          #  is prefereable (so long as rcond=1e-7 is ok in general...)
          #  Check: P'*P' = (dG (dGT dG)^1 dGT)(dG (dGT dG)^-1 dGT) = (dG (dGT dG)^1 dGT) = P'
          #invGG = _np.linalg.inv(_np.dot(_np.transpose(gen_ndG), gen_ndG))
          #Pp_alt = _np.dot(gen_ndG, _np.dot(invGG, _np.transpose(gen_ndG))) # a true projector (onto gauge space)
          #print "Pp - Pp_alt norm diff = ", _np.linalg.norm(Pp_alt - Pp)
  
          #OLD: ret = _np.identity(nParams,'d') - Pp
          # Check ranks to make sure everything is consistent.  If either of these assertions fail,
          #  then pinv's rcond or some other numerical tolerances probably need adjustment.
          #print "Rank P = ",_np.linalg.matrix_rank(P)
          #print "Rank Pp = ",_np.linalg.matrix_rank(Pp, P_RANK_TOL)
          #print "Rank 1-Pp = ",_np.linalg.matrix_rank(_np.identity(nParams,'d') - Pp, P_RANK_TOL)
          #print " Evals(1-Pp) = \n","\n".join([ "%d: %g" % (i,ev) \
          #       for i,ev in enumerate(_np.sort(_np.linalg.eigvals(_np.identity(nParams,'d') - Pp))) ])
  
          try:
              rank_P = _np.linalg.matrix_rank(P, P_RANK_TOL)  # original un-normalized projector
  
              # Note: use P_RANK_TOL here even though projector is *un-normalized* since sometimes P will
              #  have eigenvalues 1e-17 and one or two 1e-11 that should still be "zero" but aren't when
              #  no tolerance is given.  Perhaps a more custom tolerance based on the singular values of P
              #  but different from numpy's default tolerance would be appropriate here.
  
              assert(rank_P == _np.linalg.matrix_rank(Pp, P_RANK_TOL))  # rank shouldn't change with normalization
              #assert( (nParams - rank_P) == _np.linalg.matrix_rank(ret, P_RANK_TOL) ) # dimension of orthogonal space
          except(_np.linalg.LinAlgError):
              _warnings.warn("Linear algebra error (probably a non-convergent"
                             "SVD) ignored during matric rank checks in "
                             "Model.nongauge_projector(...) ")
  
          return Pp
code_tokens: def nongauge_projector ( self , item_weights = None , non_gauge_mix_mx = None ) : nParams = self . Np dPG = self . _buildup_dpg ( ) nullsp = _mt . nullspace_qr ( dPG ) gen_dG = nullsp [ 0 : nParams , : ] if non_gauge_mix_mx is not None : msg = "You've set both non_gauge_mix_mx and item_weights, both of which" + " set the gauge metric... You probably don't want to do this." assert ( item_weights is None ) , msg nonGaugeDirections = _mt . nullspace_qr ( gen_dG . T ) gen_dG = gen_dG + _np . dot ( nonGaugeDirections , non_gauge_mix_mx ) if item_weights is not None : metric_diag = _np . ones ( self . Np , 'd' ) opWeight = item_weights . get ( 'gates' , 1.0 ) spamWeight = item_weights . get ( 'spam' , 1.0 ) for lbl , gate in self . operations . items ( ) : metric_diag [ gate . gpindices ] = item_weights . get ( lbl , opWeight ) for lbl , vec in _itertools . chain ( iter ( self . preps . items ( ) ) , iter ( self . effects . items ( ) ) ) : metric_diag [ vec . gpindices ] = item_weights . get ( lbl , spamWeight ) metric = _np . diag ( metric_diag ) gen_ndG = _mt . nullspace_qr ( _np . dot ( gen_dG . T , metric ) ) else : gen_ndG = _mt . nullspace_qr ( gen_dG . T ) P = _np . dot ( gen_ndG , _np . transpose ( gen_ndG ) ) Pp = _np . dot ( _np . linalg . pinv ( P , rcond = 1e-7 ) , P ) try : rank_P = _np . linalg . matrix_rank ( P , P_RANK_TOL ) assert ( rank_P == _np . linalg . matrix_rank ( Pp , P_RANK_TOL ) ) except ( _np . linalg . LinAlgError ) : _warnings . warn ( "Linear algebra error (probably a non-convergent" "SVD) ignored during matric rank checks in " "Model.nongauge_projector(...) " ) return Pp
short_docstring: Constructs a projector onto the non-gauge parameter space.
short_docstring_tokens: Constructs a projector onto the non - gauge parameter space .
comment: """
          Constructs a projector onto the non-gauge parameter space.
  
          This is useful for isolating the gauge degrees of freedom from the non-gauge
          degrees of freedom.
  
          Parameters
          ----------
          item_weights : dict, optional
              Dictionary of weighting factors for individual gates and spam operators.
              Keys can be gate, state preparation, POVM effect, spam labels, or the
              special strings "gates" or "spam" whic represent the entire set of gate
              or SPAM operators, respectively.  Values are floating point numbers.
              These weights define the metric used to compute the non-gauge space,
              *orthogonal* the gauge space, that is projected onto.
  
          non_gauge_mix_mx : numpy array, optional
              An array of shape (n_non_gauge_params,n_gauge_params) specifying how to
              mix the non-gauge degrees of freedom into the gauge degrees of
              freedom that are projected out by the returned object.  This argument
              essentially sets the off-diagonal block of the metric used for
              orthogonality in the "gauge + non-gauge" space.  It is for advanced
              usage and typically left as None (the default).
  
          Returns
          -------
          numpy array
              The projection operator as a N x N matrix, where N is the number
              of parameters (obtained via num_params()).  This projector acts on
              parameter-space, and has rank equal to the number of non-gauge
              degrees of freedom.
          """$SEP$# We want to divide the Model-space H (a Hilbert space, 56-dimensional in the 1Q, 3-gate, 2-vec case)$SEP$# into the direct sum of gauge and non-gauge spaces, and find projectors onto each$SEP$# sub-space (the non-gauge space in particular).$SEP$#$SEP$# Within the Model H-space lies a gauge-manifold of maximum chi2 (16-dimensional in 1Q case),$SEP$#  where the gauge-invariant chi2 function is constant.  At a single point (Model) P on this manifold,$SEP$#  chosen by (somewhat arbitrarily) fixing the gauge, we can talk about the tangent space$SEP$#  at that point.  This tangent space is spanned by some basis (~16 elements for the 1Q case),$SEP$#  which associate with the infinitesimal gauge transform ?generators? on the full space.$SEP$#  The subspace of H spanned by the derivatives w.r.t gauge transformations at P (a Model) spans$SEP$#  the gauge space, and the complement of this (in H) is the non-gauge space.$SEP$#$SEP$#  An element of the gauge group can be written gg = exp(-K), where K is a n x n matrix.  If K is$SEP$#   hermitian then gg is unitary, but this need not be the case.  A gauge transform acts on a$SEP$#   gatset via Model => gg^-1 G gg, i.e. G => exp(K) G exp(-K).  We care about the action of$SEP$#   infinitesimal gauge tranformations (b/c the *derivative* vectors span the tangent space),$SEP$#   which act as:$SEP$#    G => (I+K) G (I-K) = G + [K,G] + ignored(K^2), where [K,G] := KG-GK$SEP$#$SEP$# To produce a projector onto the gauge-space, we compute the *column* vectors$SEP$#  dG_ij = [K_ij,G], where K_ij is the i,j-th matrix unit (56x1 in the 1Q case, 16 such column vectors)$SEP$#  and then form a projector in the standard way.$SEP$#  (to project onto |v1>, |v2>, etc., form P = sum_i |v_i><v_i|)$SEP$#$SEP$#Typically nOpParams < len(dG_ij) and linear system is overconstrained$SEP$#   and no solution is expected.  If no solution exists, simply ignore$SEP$#$SEP$# So we form P = sum_ij dG_ij * transpose(dG_ij) (56x56 in 1Q case)$SEP$#              = dG * transpose(dG)              where dG_ij form the *columns* of dG (56x16 in 1Q case)$SEP$# But dG_ij are not orthonormal, so really we need a slight modification,$SEP$#  otherwise P'P' != P' as must be the case for a projector:$SEP$#$SEP$# P' = dG * (transpose(dG) * dG)^-1 * transpose(dG) (see wikipedia on projectors)$SEP$#$SEP$#    or equivalently (I think)$SEP$#$SEP$# P' = pseudo-inv(P)*P$SEP$#$SEP$#  since the pseudo-inv is defined by P*pseudo-inv(P) = I and so P'P' == P'$SEP$#  and P' is our gauge-projector!$SEP$# Note: In the general case of parameterized gates (i.e. non-fully parameterized gates), there are fewer$SEP$#   gate parameters than the size of dG_ij ( < 56 in 1Q case).  In this general case, we want to know$SEP$#   what (if any) change in gate parameters produces the change dG_ij of the operation matrix elements.$SEP$#   That is, we solve dG_ij = derivWRTParams * dParams_ij for dParams_ij, where derivWRTParams is$SEP$#   the derivative of the operation matrix elements with respect to the gate parameters (derivWRTParams$SEP$#   is 56x(nOpParams) and dParams_ij is (nOpParams)x1 in the 1Q case) -- then substitute dG_ij$SEP$#   with dParams_ij above to arrive at a (nOpParams)x(nOpParams) projector (compatible with$SEP$#   hessian computed by model).$SEP$#$SEP$#   Actually, we want to know if what changes gate parameters$SEP$#   produce changes in the span of all the dG_ij, so really we want the intersection of the space$SEP$#   defined by the columns of derivWRTParams (the "gate-parameter range" space) and the dG_ij vectors.$SEP$#$SEP$#   This intersection is determined by nullspace( derivWRTParams | dG ), where the pipe denotes$SEP$#     concatenating the matrices together.  Then, if x is a basis vector of the nullspace$SEP$#     x[0:nOpParams] is the basis vector for the intersection space within the gate parameter space,$SEP$#     that is, the analogue of dParams_ij in the single-dG_ij introduction above.$SEP$#$SEP$#   Still, we just substitue these dParams_ij vectors (as many as the nullspace dimension) for dG_ij$SEP$#   above to get the general case projector.$SEP$#print("DB: shapes = ",dP.shape,dG.shape)$SEP$# columns are nullspace basis vectors$SEP$# take upper (gate-param-segment) of vectors for basis$SEP$# of subspace intersection in gate-parameter space$SEP$#Note: gen_dG == "generalized dG", and is (nParams)x(nullSpaceDim==gaugeSpaceDim), so P$SEP$#  (below) is (nParams)x(nParams) as desired.$SEP$#DEBUG$SEP$#nElements = self.num_elements$SEP$#for iRow in range(nElements):$SEP$#    pNorm = _np.linalg.norm(dP[iRow])$SEP$#    if pNorm < 1e-6:$SEP$#        print "Row %d of dP is zero" % iRow$SEP$#$SEP$#print "------------------------------"$SEP$#print "nParams = ",nParams$SEP$#print "nElements = ",nElements$SEP$#print " shape(M) = ",M.shape$SEP$#print " shape(dG) = ",dG.shape$SEP$#print " shape(dP) = ",dP.shape$SEP$#print " shape(gen_dG) = ",gen_dG.shape$SEP$#print " shape(nullsp) = ",nullsp.shape$SEP$#print " rank(dG) = ",_np.linalg.matrix_rank(dG)$SEP$#print " rank(dP) = ",_np.linalg.matrix_rank(dP)$SEP$#print "------------------------------"$SEP$#assert(_np.linalg.norm(_np.imag(gen_dG)) < 1e-9) #gen_dG is real$SEP$# BEGIN GAUGE MIX ----------------------------------------$SEP$# nullspace of gen_dG^T (mx with gauge direction vecs as rows) gives non-gauge directions$SEP$# columns are the non-gauge directions *orthogonal* to the gauge directions$SEP$#for each column of gen_dG, which is a gauge direction in model parameter space,$SEP$# we add some amount of non-gauge direction, given by coefficients of the$SEP$# numNonGaugeParams non-gauge directions.$SEP$# add non-gauge direction components in$SEP$# dims: (nParams,n_gauge_params) + (nParams,n_non_gauge_params) * (n_non_gauge_params,n_gauge_params)$SEP$# non_gauge_mix_mx is a (n_non_gauge_params,n_gauge_params) matrix whose i-th column specifies the$SEP$#  coefficents to multipy each of the non-gauge directions by before adding them to the i-th$SEP$#  direction to project out (i.e. what were the pure gauge directions).$SEP$#DEBUG$SEP$#print "gen_dG shape = ",gen_dG.shape$SEP$#print "NGD shape = ",nonGaugeDirections.shape$SEP$#print "NGD rank = ",_np.linalg.matrix_rank(nonGaugeDirections, P_RANK_TOL)$SEP$#END GAUGE MIX ----------------------------------------$SEP$# Build final non-gauge projector by getting a mx of column vectors$SEP$# orthogonal to the cols fo gen_dG:$SEP$#     gen_dG^T * gen_ndG = 0 => nullspace(gen_dG^T)$SEP$# or for a general metric W:$SEP$#     gen_dG^T * W * gen_ndG = 0 => nullspace(gen_dG^T * W)$SEP$# (This is instead of construction as I - gaugeSpaceProjector)$SEP$#OLD: gen_ndG = _mt.nullspace(_np.dot(gen_dG.T,metric))$SEP$#OLD: gen_ndG = _mt.nullspace(gen_dG.T) #cols are non-gauge directions$SEP$# cols are non-gauge directions$SEP$# print("DB: nullspace of gen_dG (shape = %s, rank=%d) = %s" \$SEP$#       % (str(gen_dG.shape),_np.linalg.matrix_rank(gen_dG),str(gen_ndG.shape)))$SEP$# ORIG WAY: use psuedo-inverse to normalize projector.  Ran into problems where$SEP$#  default rcond == 1e-15 didn't work for 2-qubit case, but still more stable than inv method below$SEP$# almost a projector, but cols of dG are not orthonormal$SEP$# make P into a true projector (onto gauge space)$SEP$# ALT WAY: use inverse of dG^T*dG to normalize projector (see wikipedia on projectors, dG => A)$SEP$#  This *should* give the same thing as above, but numerical differences indicate the pinv method$SEP$#  is prefereable (so long as rcond=1e-7 is ok in general...)$SEP$#  Check: P'*P' = (dG (dGT dG)^1 dGT)(dG (dGT dG)^-1 dGT) = (dG (dGT dG)^1 dGT) = P'$SEP$#invGG = _np.linalg.inv(_np.dot(_np.transpose(gen_ndG), gen_ndG))$SEP$#Pp_alt = _np.dot(gen_ndG, _np.dot(invGG, _np.transpose(gen_ndG))) # a true projector (onto gauge space)$SEP$#print "Pp - Pp_alt norm diff = ", _np.linalg.norm(Pp_alt - Pp)$SEP$#OLD: ret = _np.identity(nParams,'d') - Pp$SEP$# Check ranks to make sure everything is consistent.  If either of these assertions fail,$SEP$#  then pinv's rcond or some other numerical tolerances probably need adjustment.$SEP$#print "Rank P = ",_np.linalg.matrix_rank(P)$SEP$#print "Rank Pp = ",_np.linalg.matrix_rank(Pp, P_RANK_TOL)$SEP$#print "Rank 1-Pp = ",_np.linalg.matrix_rank(_np.identity(nParams,'d') - Pp, P_RANK_TOL)$SEP$#print " Evals(1-Pp) = \n","\n".join([ "%d: %g" % (i,ev) \$SEP$#       for i,ev in enumerate(_np.sort(_np.linalg.eigvals(_np.identity(nParams,'d') - Pp))) ])$SEP$# original un-normalized projector$SEP$# Note: use P_RANK_TOL here even though projector is *un-normalized* since sometimes P will$SEP$#  have eigenvalues 1e-17 and one or two 1e-11 that should still be "zero" but aren't when$SEP$#  no tolerance is given.  Perhaps a more custom tolerance based on the singular values of P$SEP$#  but different from numpy's default tolerance would be appropriate here.$SEP$# rank shouldn't change with normalization$SEP$#assert( (nParams - rank_P) == _np.linalg.matrix_rank(ret, P_RANK_TOL) ) # dimension of orthogonal space