Python numpy.number() Examples

def set_input_map(self, input_map):
        """Set the input map for an interconnected I/O system.

        Parameters
        ----------
        input_map : 2D array
             Specify the matrix that will be used to multiply the vector of
             system inputs to obtain the vector of subsystem inputs.  These
             values are added to the inputs specified in the connection map.

        """
        # Figure out the number of internal inputs
        ninputs = sum(sys.ninputs for sys in self.syslist)

        # Make sure the input map is the right size
        if input_map.shape[0] != ninputs:
            ValueError("Input map is not the right shape")
        self.input_map = input_map
        self.ninputs = input_map.shape[1] 

def __truediv__(self, other):
        """Divide two LTI objects."""

        if isinstance(other, (int, float, complex, np.number)):
            return FRD(self.fresp * (1/other), self.omega,
                       smooth=(self.ifunc is not None))
        else:
            other = _convertToFRD(other, omega=self.omega)

        if (self.inputs > 1 or self.outputs > 1 or
            other.inputs > 1 or other.outputs > 1):
            raise NotImplementedError(
                "FRD.__truediv__ is currently only implemented for SISO "
                "systems.")

        return FRD(self.fresp/other.fresp, self.omega,
                   smooth=(self.ifunc is not None) and
                          (other.ifunc is not None))

    # TODO: Remove when transition to python3 complete 

def one_from_two_pdm(two_pdm, nelec):
    '''Return a 1-rdm, given a 2-rdm to contract.

    Args:
        two_pdm : ndarray
            A (spin-free) 2-particle reduced density matrix.
        nelec: int
            The number of electrons contributing to the RDMs.

    Returns:
        one_pdm : ndarray
            The (spin-free) 1-particle reduced density matrix.
    '''

    # Last two indices refer to middle two second quantized operators in the 2RDM
    one_pdm = numpy.einsum('ikjj->ik', two_pdm)
    one_pdm /= (numpy.sum(nelec)-1)
    return one_pdm 

def one_from_two_pdm(two_pdm, nelec):
    '''Return a 1-rdm, given a 2-rdm to contract.

    Args:
        two_pdm : ndarray
            A (spin-free) 2-particle reduced density matrix.
        nelec: int
            The number of electrons contributing to the RDMs.

    Returns:
        one_pdm : ndarray
            The (spin-free) 1-particle reduced density matrix.
    '''

    # Last two indices refer to middle two second quantized operators in the 2RDM
    one_pdm = numpy.einsum('ikjj->ik', two_pdm)
    one_pdm /= (numpy.sum(nelec)-1)
    return one_pdm 

def test_ticket_1539(self):
        dtypes = [x for x in np.typeDict.values()
                  if (issubclass(x, np.number)
                      and not issubclass(x, np.timedelta64))]
        a = np.array([], np.bool_)  # not x[0] because it is unordered
        failures = []

        for x in dtypes:
            b = a.astype(x)
            for y in dtypes:
                c = a.astype(y)
                try:
                    np.dot(b, c)
                except TypeError:
                    failures.append((x, y))
        if failures:
            raise AssertionError("Failures: %r" % failures) 

def revisitFilter(self,sInds,tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead was added in MJD form
        Returns:
            sInds - indices of stars still in observation list
        """

        tovisit = np.zeros(self.TargetList.nStars, dtype=bool)
        if len(sInds) > 0:
            tovisit[sInds] = ((self.starVisits[sInds] == min(self.starVisits[sInds])) \
                    & (self.starVisits[sInds] < self.nVisitsMax))#Checks that no star has exceeded the number of revisits and the indicies of all considered stars have minimum number of observations
            #The above condition should prevent revisits so long as all stars have not been observed
            if self.starRevisit.size != 0:
                dt_rev = np.abs(self.starRevisit[:,1]*u.day - tmpCurrentTimeNorm)
                ind_rev = [int(x) for x in self.starRevisit[dt_rev < self.dt_max,0] 
                        if x in sInds]
                tovisit[ind_rev] = (self.starVisits[ind_rev] < self.nVisitsMax)
            sInds = np.where(tovisit)[0]
        return sInds 

def evalfr(self, omega):
        """Evaluate a transfer function at a single angular frequency.

        self._evalfr(omega) returns the value of the frequency response
        at frequency omega.

        Note that a "normal" FRD only returns values for which there is an
        entry in the omega vector. An interpolating FRD can return
        intermediate values.

        """
        warn("FRD.evalfr(omega) will be deprecated in a future release "
             "of python-control; use sys.eval(omega) instead",
             PendingDeprecationWarning)         # pragma: no coverage
        return self._evalfr(omega)

    # Define the `eval` function to evaluate an FRD at a given (real)
    # frequency.  Note that we choose to use `eval` instead of `evalfr` to
    # avoid confusion with :func:`evalfr`, which takes a complex number as its
    # argument.  Similarly, we don't use `__call__` to avoid confusion between
    # G(s) for a transfer function and G(omega) for an FRD object. 

def set_states(self, states, prefix='x'):
        """Set the number/names of the system states.

        Parameters
        ----------
        states : int, list of str, or None
            Description of the system states.  This can be given as an integer
            count or as a list of strings that name the individual signals.
            If an integer count is specified, the names of the signal will be
            of the form `u[i]` (where the prefix `u` can be changed using the
            optional prefix parameter).
        prefix : string, optional
            If `states` is an integer, create the names of the states using
            the given prefix (default = 'x').  The names of the input will be
            of the form `prefix[i]`.

        """
        self.nstates, self.state_index = \
            self._process_signal_list(states, prefix=prefix) 

def set_outputs(self, outputs, prefix='y'):
        """Set the number/names of the system outputs.

        Parameters
        ----------
        outputs : int, list of str, or None
            Description of the system outputs.  This can be given as an integer
            count or as a list of strings that name the individual signals.
            If an integer count is specified, the names of the signal will be
            of the form `u[i]` (where the prefix `u` can be changed using the
            optional prefix parameter).
        prefix : string, optional
            If `outputs` is an integer, create the names of the states using
            the given prefix (default = 'y').  The names of the input will be
            of the form `prefix[i]`.

        """
        self.noutputs, self.output_index = \
            self._process_signal_list(outputs, prefix=prefix) 

def __rmul__(sys1, sys2):
        """Pre-multiply an input/output systems by a scalar/matrix"""
        if isinstance(sys2, (int, float, np.number)):
            # TODO: Scale the output
            raise NotImplemented("Scalar multiplication not yet implemented")
        elif isinstance(sys2, np.ndarray):
            # TODO: Post-multiply by a matrix
            raise NotImplemented("Matrix multiplication not yet implemented")
        elif isinstance(sys1, StateSpace) and isinstance(sys2, StateSpace):
            # Special case: maintain linear systems structure
            new_ss_sys = StateSpace.__rmul__(sys1, sys2)
            # TODO: set input and output names
            new_io_sys = LinearIOSystem(new_ss_sys)

            return new_io_sys
        elif not isinstance(sys2, InputOutputSystem):
            raise ValueError("Unknown I/O system object ", sys1)
        else:
            # Both systetms are InputOutputSystems => use __mul__
            return InputOutputSystem.__mul__(sys2, sys1) 

def issiso(sys, strict=False):
    """
    Check to see if a system is single input, single output

    Parameters
    ----------
    sys : LTI system
        System to be checked
    strict: bool (default = False)
        If strict is True, do not treat scalars as SISO
    """
    if isinstance(sys, (int, float, complex, np.number)) and not strict:
        return True
    elif not isinstance(sys, LTI):
        raise ValueError("Object is not an LTI system")

    # Done with the tricky stuff...
    return sys.issiso()

# Return the timebase (with conversion if unspecified) 

def kernel(self, h1e, eri, norb, nelec, ci0=None,
               tol=None, lindep=None, max_cycle=None, max_space=None,
               nroots=None, davidson_only=None, pspace_size=None,
               orbsym=None, wfnsym=None, ecore=0, **kwargs):
        if isinstance(nelec, (int, numpy.number)):
            nelecb = nelec//2
            neleca = nelec - nelecb
        else:
            neleca, nelecb = nelec
        link_indexa = cistring.gen_linkstr_index(range(norb), neleca)
        link_indexb = cistring.gen_linkstr_index(range(norb), nelecb)
        e, c = direct_spin1.kernel_ms1(self, h1e, eri, norb, nelec, ci0,
                                       (link_indexa,link_indexb),
                                       tol, lindep, max_cycle, max_space, nroots,
                                       davidson_only, pspace_size, ecore=ecore,
                                       **kwargs)
        self.eci, self.ci = e, c
        return e, c 

def timebase(sys, strict=True):
    """Return the timebase for an LTI system

    dt = timebase(sys)

    returns the timebase for a system 'sys'.  If the strict option is
    set to False, dt = True will be returned as 1.
    """
    # System needs to be either a constant or an LTI system
    if isinstance(sys, (int, float, complex, np.number)):
        return None
    elif not isinstance(sys, LTI):
        raise ValueError("Timebase not defined")

    # Return the sample time, with converstion to float if strict is false
    if (sys.dt == None):
        return None
    elif (strict):
        return float(sys.dt)

    return sys.dt

# Check to see if two timebases are equal 

def from_fcivec(ci0, norb, nelec, frozen=None):
    '''Extract CISD coefficients from FCI coefficients'''
    if not (frozen is None or frozen == 0):
        raise NotImplementedError

    if isinstance(nelec, (int, numpy.number)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    nocc = neleca
    nvir = norb - nocc
    t1addr, t1sign = t1strs(norb, nocc)

    c0 = ci0[0,0]
    c1 = ci0[0,t1addr] * t1sign
    c2 = numpy.einsum('i,j,ij->ij', t1sign, t1sign, ci0[t1addr[:,None],t1addr])
    c1 = c1.reshape(nocc,nvir)
    c2 = c2.reshape(nocc,nvir,nocc,nvir).transpose(0,2,1,3)
    return amplitudes_to_cisdvec(c0, c1, c2) 

def write_ci(fout, fcivec, norb, nelec, ncore=0):
    from pyscf import fci
    if isinstance(nelec, (int, numpy.number)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    fout.write('NELACTIVE, NDETS, NORBCORE, NORBACTIVE\n')
    fout.write(' %5d %5d %5d %5d\n' % (neleca+nelecb, fcivec.size, ncore, norb))
    fout.write('COEFFICIENT/ OCCUPIED ACTIVE SPIN ORBITALS\n')

    nb = fci.cistring.num_strings(norb, nelecb)
    stringsa = fci.cistring.gen_strings4orblist(range(norb), neleca)
    stringsb = fci.cistring.gen_strings4orblist(range(norb), nelecb)
    def str2orbidx(string, ncore):
        bstring = bin(string)
        return [i+1+ncore for i,s in enumerate(bstring[::-1]) if s == '1']

    addrs = numpy.argsort(abs(fcivec.ravel()))
    for iaddr in reversed(addrs):
        addra, addrb = divmod(iaddr, nb)
        idxa = ['%3d' % x for x in str2orbidx(stringsa[addra], ncore)]
        idxb = ['%3d' % (-x) for x in str2orbidx(stringsb[addrb], ncore)]
        #TODO:add a cuttoff and a counter for ndets
        fout.write('%18.10E %s %s\n' % (fcivec[addra,addrb], ' '.join(idxa), ' '.join(idxb))) 

def reorder_eri(eri, norb, orbsym):
    if orbsym is None:
        return [eri], numpy.arange(norb), numpy.zeros(norb,dtype=numpy.int32)
# map irrep IDs of Dooh or Coov to D2h, C2v
# see symm.basis.linearmole_symm_descent
    orbsym = numpy.asarray(orbsym) % 10
# irrep of (ij| pair
    trilirrep = (orbsym[:,None]^orbsym)[numpy.tril_indices(norb)]
# and the number of occurence for each irrep
    dimirrep = numpy.asarray(numpy.bincount(trilirrep), dtype=numpy.int32)
# we sort the irreps of (ij| pair, to group the pairs which have same irreps
# "order" is irrep-id-sorted index. The (ij| paired is ordered that the
# pair-id given by order[0] comes first in the sorted pair
# "rank" is a sorted "order". Given nth (ij| pair, it returns the place(rank)
# of the sorted pair
    old_eri_irrep = numpy.asarray(trilirrep, dtype=numpy.int32)
    rank_in_irrep = numpy.empty_like(old_eri_irrep)
    p0 = 0
    eri_irs = [numpy.zeros((0,0))] * TOTIRREPS
    for ir, nnorb in enumerate(dimirrep):
        idx = numpy.asarray(numpy.where(trilirrep == ir)[0], dtype=numpy.int32)
        rank_in_irrep[idx] = numpy.arange(nnorb, dtype=numpy.int32)
        eri_irs[ir] = lib.take_2d(eri, idx, idx)
        p0 += nnorb
    return eri_irs, rank_in_irrep, old_eri_irrep 

def __add__(self, other):
        """Add two LTI systems (parallel connection)."""

        # Check for a couple of special cases
        if isinstance(other, (int, float, complex, np.number)):
            # Just adding a scalar; put it in the D matrix
            A, B, C = self.A, self.B, self.C
            D = self.D + other
            dt = self.dt
        else:
            other = _convertToStateSpace(other)

            # Check to make sure the dimensions are OK
            if ((self.inputs != other.inputs) or
                    (self.outputs != other.outputs)):
                raise ValueError("Systems have different shapes.")

            # Figure out the sampling time to use
            if self.dt is None and other.dt is not None:
                dt = other.dt       # use dt from second argument
            elif (other.dt is None and self.dt is not None) or \
                    (timebaseEqual(self, other)):
                dt = self.dt        # use dt from first argument
            else:
                raise ValueError("Systems have different sampling times")

            # Concatenate the various arrays
            A = concatenate((
                concatenate((self.A, zeros((self.A.shape[0],
                                           other.A.shape[-1]))),axis=1),
                concatenate((zeros((other.A.shape[0], self.A.shape[-1])),
                                other.A),axis=1)
                            ),axis=0)
            B = concatenate((self.B, other.B), axis=0)
            C = concatenate((self.C, other.C), axis=1)
            D = self.D + other.D

        return StateSpace(A, B, C, D, dt)

    # Right addition - just switch the arguments 

def __mul__(self, other):
        """Multiply two LTI objects (serial connection)."""

        # Convert the second argument to a transfer function.
        if isinstance(other, (int, float, complex, np.number)):
            return FRD(self.fresp * other, self.omega,
                       smooth=(self.ifunc is not None))
        else:
            other = _convertToFRD(other, omega=self.omega)

        # Check that the input-output sizes are consistent.
        if self.inputs != other.outputs:
            raise ValueError(
                "H = G1*G2: input-output size mismatch: "
                "G1 has %i input(s), G2 has %i output(s)." %
                (self.inputs, other.outputs))

        inputs = other.inputs
        outputs = self.outputs
        fresp = empty((outputs, inputs, len(self.omega)),
                      dtype=self.fresp.dtype)
        for i in range(len(self.omega)):
            fresp[:, :, i] = dot(self.fresp[:, :, i], other.fresp[:, :, i])
        return FRD(fresp, self.omega,
                   smooth=(self.ifunc is not None) and
                          (other.ifunc is not None)) 

def __add__(sys1, sys2):
        """Add two input/output systems (parallel interconnection)"""
        # TODO: Allow addition of scalars and matrices
        if not isinstance(sys2, InputOutputSystem):
            raise ValueError("Unknown I/O system object ", sys2)
        elif isinstance(sys1, StateSpace) and isinstance(sys2, StateSpace):
            # Special case: maintain linear systems structure
            new_ss_sys = StateSpace.__add__(sys1, sys2)
            # TODO: set input and output names
            new_io_sys = LinearIOSystem(new_ss_sys)

            return new_io_sys

        # Make sure number of input and outputs match
        if sys1.ninputs != sys2.ninputs or sys1.noutputs != sys2.noutputs:
            raise ValueError("Can't add systems with different numbers of "
                             "inputs or outputs.")
        ninputs = sys1.ninputs
        noutputs = sys1.noutputs

        # Create a new system to handle the composition
        newsys = InterconnectedSystem((sys1, sys2))

        # Set up the input map
        newsys.set_input_map(np.concatenate(
            (np.eye(ninputs), np.eye(ninputs)), axis=0))
        # TODO: set up input names

        # Set up the output map
        newsys.set_output_map(np.concatenate(
            (np.eye(noutputs), np.eye(noutputs)), axis=1))
        # TODO: set up output names

        # Return the newly created system
        return newsys

    # TODO: add __radd__ to allow postaddition by scalars and matrices 

def __truediv__(self, other):
        """Divide two LTI objects."""

        if isinstance(other, (int, float, complex, np.number)):
            other = _convert_to_transfer_function(
                other, inputs=self.inputs,
                outputs=self.inputs)
        else:
            other = _convert_to_transfer_function(other)

        if (self.inputs > 1 or self.outputs > 1 or
                other.inputs > 1 or other.outputs > 1):
            raise NotImplementedError(
                "TransferFunction.__truediv__ is currently \
                implemented only for SISO systems.")

        # Figure out the sampling time to use
        if self.dt is None and other.dt is not None:
            dt = other.dt       # use dt from second argument
        elif (other.dt is None and self.dt is not None) or \
             (self.dt == other.dt):
            dt = self.dt        # use dt from first argument
        else:
            raise ValueError("Systems have different sampling times")

        num = polymul(self.num[0][0], other.den[0][0])
        den = polymul(self.den[0][0], other.num[0][0])

        return TransferFunction(num, den, dt)

    # TODO: Remove when transition to python3 complete 

def __neg__(sys):
        """Negate an input/output systems (rescale)"""
        if isinstance(sys, StateSpace):
            # Special case: maintain linear systems structure
            new_ss_sys = StateSpace.__neg__(sys)
            # TODO: set input and output names
            new_io_sys = LinearIOSystem(new_ss_sys)

            return new_io_sys
        if sys.ninputs is None or sys.noutputs is None:
            raise ValueError("Can't determine number of inputs or outputs")

        # Create a new system to hold the negation
        newsys = InterconnectedSystem((sys,), dt=sys.dt)

        # Set up the input map (identity)
        newsys.set_input_map(np.eye(sys.ninputs))
        # TODO: set up input names

        # Set up the output map (negate the output)
        newsys.set_output_map(-np.eye(sys.noutputs))
        # TODO: set up output names

        # Return the newly created system
        return newsys

    # Utility function to parse a list of signals 

def __rmul__(self, other):
        """Right multiply two LTI objects (serial connection)."""

        # Check for a couple of special cases
        if isinstance(other, (int, float, complex, np.number)):
            # Just multiplying by a scalar; change the input
            A, C = self.A, self.C
            B = self.B * other
            D = self.D * other
            return StateSpace(A, B, C, D, self.dt)

        # is lti, and convertible?
        if isinstance(other, LTI):
            return _convertToStateSpace(other) * self

        # try to treat this as a matrix
        try:
            X = _ssmatrix(other)
            C = np.dot(X, self.C)
            D = np.dot(X, self.D)
            return StateSpace(self.A, self.B, C, D, self.dt)

        except Exception as e:
            print(e)
            pass
        raise TypeError("can't interconnect systems")

    # TODO: __div__ and __rdiv__ are not written yet. 

def _finalize(self):
        log = logger.Logger(self.stdout, self.verbose)
        if log.verbose >= logger.NOTE and getattr(self.fcisolver, 'spin_square', None):
            ncore = self.ncore
            ncas = self.ncas
            mocas = (self.mo_coeff[0][:,ncore[0]:ncore[0]+ncas],
                     self.mo_coeff[1][:,ncore[1]:ncore[1]+ncas])
            mocore = (self.mo_coeff[0][:,:ncore[0]],
                      self.mo_coeff[1][:,:ncore[1]])
            ovlp_ao = self._scf.get_ovlp()
            ss_core = self._scf.spin_square(mocore, ovlp_ao)
            if isinstance(self.e_cas, (float, numpy.number)):
                ss = fci.spin_op.spin_square(self.ci, self.ncas, self.nelecas,
                                             mocas, ovlp_ao)
                log.note('UCASCI E = %.15g  E(CI) = %.15g  S^2 = %.7f',
                         self.e_tot, self.e_cas, ss[0]+ss_core[0])
            else:
                for i, e in enumerate(self.e_cas):
                    ss = fci.spin_op.spin_square(self.ci[i], self.ncas, self.nelecas,
                                                 mocas, ovlp_ao)
                    log.note('UCASCI root %d  E = %.15g  E(CI) = %.15g  S^2 = %.7f',
                             i, self.e_tot[i], e, ss[0]+ss_core[0])
        else:
            if isinstance(self.e_cas, (float, numpy.number)):
                log.note('UCASCI E = %.15g  E(CI) = %.15g', self.e_tot, self.e_cas)
            else:
                for i, e in enumerate(self.e_cas):
                    log.note('UCASCI root %d  E = %.15g  E(CI) = %.15g',
                             i, self.e_tot[i], e)
        #if self.verbose >= logger.INFO:
        #    self.analyze(mo_coeff, self.ci, verbose=self.verbose)
        return self 

def fix_spin_(self, shift=PENALTY, ss=None):
        r'''Use level shift to control FCI solver spin.

        .. math::

            (H + shift*S^2) |\Psi\rangle = E |\Psi\rangle

        Kwargs:
            shift : float
                Energy penalty for states which have wrong spin
            ss : number
                S^2 expection value == s*(s+1)
        '''
        fci.addons.fix_spin_(self.fcisolver, shift, ss)
        return self 

def tn_addrs_signs(norb, nelec, n_excite):
    '''Compute the FCI strings (address) for CIS n-excitation amplitudes and
    the signs of the coefficients when transferring the reference from physics
    vacuum to HF vacuum.
    '''
    if n_excite > nelec:
        print("Warning: Not enough occupied orbitals to excite.")
        return [0], [0]
    nocc = nelec

    hole_strs = cistring.gen_strings4orblist(range(nocc), nocc - n_excite)
    # For HF vacuum, hole operators are ordered from low-lying to high-lying
    # orbitals. It leads to the opposite string ordering.
    hole_strs = hole_strs[::-1]
    hole_sum = numpy.zeros(len(hole_strs), dtype=int)
    for i in range(nocc):
        hole_at_i = (hole_strs & (1<<i)) == 0
        hole_sum[hole_at_i] += i

    # The hole operators are listed from low-lying to high-lying orbitals
    # (from left to right).  For i-th (0-based) hole operator, the number of
    # orbitals which are higher than i determines the sign.  This number
    # equals to nocc-(i+1).  After removing the highest hole operator, nocc
    # becomes nocc-1, the sign for next hole operator j will be associated to
    # nocc-1-(j+1).  By iteratively calling this procedure, the overall sign
    # for annihilating three holes is (-1)**(3*nocc - 6 - sum i)
    sign = (-1) ** (n_excite * nocc - n_excite*(n_excite+1)//2 - hole_sum)

    particle_strs = cistring.gen_strings4orblist(range(nocc, norb), n_excite)
    strs = hole_strs[:,None] ^ particle_strs
    addrs = cistring.strs2addr(norb, nocc, strs.ravel())
    signs = numpy.vstack([sign] * len(particle_strs)).T.ravel()
    return addrs, signs 

def from_fcivec(ci0, norb, nelec, frozen=None):
    '''Extract CISD coefficients from FCI coefficients'''
    if not (frozen is None or frozen == 0):
        raise NotImplementedError

    if isinstance(nelec, (int, numpy.number)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec

    norba = norbb = norb
    nocca, noccb = neleca, nelecb
    nvira = norba - nocca
    nvirb = norbb - noccb
    t1addra, t1signa = cisd.tn_addrs_signs(norba, nocca, 1)
    t1addrb, t1signb = cisd.tn_addrs_signs(norbb, noccb, 1)

    na = cistring.num_strings(norba, nocca)
    nb = cistring.num_strings(norbb, noccb)
    ci0 = ci0.reshape(na,nb)
    c0 = ci0[0,0]
    c1a = (ci0[t1addra,0] * t1signa).reshape(nocca,nvira)
    c1b = (ci0[0,t1addrb] * t1signb).reshape(noccb,nvirb)

    c2ab = numpy.einsum('i,j,ij->ij', t1signa, t1signb, ci0[t1addra[:,None],t1addrb])
    c2ab = c2ab.reshape(nocca,nvira,noccb,nvirb).transpose(0,2,1,3)
    t2addra, t2signa = cisd.tn_addrs_signs(norba, nocca, 2)
    t2addrb, t2signb = cisd.tn_addrs_signs(norbb, noccb, 2)
    c2aa = (ci0[t2addra,0] * t2signa).reshape(nocca*(nocca-1)//2, nvira*(nvira-1)//2)
    c2aa = _unpack_4fold(c2aa, nocca, nvira)
    c2bb = (ci0[0,t2addrb] * t2signb).reshape(noccb*(noccb-1)//2, nvirb*(nvirb-1)//2)
    c2bb = _unpack_4fold(c2bb, noccb, nvirb)

    return amplitudes_to_cisdvec(c0, (c1a,c1b), (c2aa,c2ab,c2bb)) 

def write_head(fout, nmo, nelec, ms=0, orbsym=None):
    if not isinstance(nelec, (int, numpy.number)):
        ms = abs(nelec[0] - nelec[1])
        nelec = nelec[0] + nelec[1]
    fout.write(' &FCI NORB=%4d,NELEC=%2d,MS2=%d,\n' % (nmo, nelec, ms))
    if orbsym is not None and len(orbsym) > 0:
        fout.write('  ORBSYM=%s\n' % ','.join([str(x) for x in orbsym]))
    else:
        fout.write('  ORBSYM=%s\n' % ('1,' * nmo))
    fout.write('  ISYM=1,\n')
    fout.write(' &END\n') 

def rss(states=1, outputs=1, inputs=1):
    """
    Create a stable *continuous* random state space object.

    Parameters
    ----------
    states : integer
        Number of state variables
    inputs : integer
        Number of system inputs
    outputs : integer
        Number of system outputs

    Returns
    -------
    sys : StateSpace
        The randomly created linear system

    Raises
    ------
    ValueError
        if any input is not a positive integer

    See Also
    --------
    drss

    Notes
    -----
    If the number of states, inputs, or outputs is not specified, then the
    missing numbers are assumed to be 1.  The poles of the returned system
    will always have a negative real part.

    """

    return _rss_generate(states, inputs, outputs, 'c') 

def drss(states=1, outputs=1, inputs=1):
    """
    Create a stable *discrete* random state space object.

    Parameters
    ----------
    states : integer
        Number of state variables
    inputs : integer
        Number of system inputs
    outputs : integer
        Number of system outputs

    Returns
    -------
    sys : StateSpace
        The randomly created linear system

    Raises
    ------
    ValueError
        if any input is not a positive integer

    See Also
    --------
    rss

    Notes
    -----
    If the number of states, inputs, or outputs is not specified, then the
    missing numbers are assumed to be 1.  The poles of the returned system
    will always have a magnitude less than 1.

    """

    return _rss_generate(states, inputs, outputs, 'd') 

def evalfr(sys, x):
    """
    Evaluate the transfer function of an LTI system for a single complex
    number x.

    To evaluate at a frequency, enter x = omega*j, where omega is the
    frequency in radians

    Parameters
    ----------
    sys: StateSpace or TransferFunction
        Linear system
    x: scalar
        Complex number

    Returns
    -------
    fresp: ndarray

    See Also
    --------
    freqresp
    bode

    Notes
    -----
    This function is a wrapper for StateSpace.evalfr and
    TransferFunction.evalfr.

    Examples
    --------
    >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
    >>> evalfr(sys, 1j)
    array([[ 44.8-21.4j]])
    >>> # This is the transfer function matrix evaluated at s = i.

    .. todo:: Add example with MIMO system
    """
    if issiso(sys):
        return sys.horner(x)[0][0]
    return sys.horner(x)