def estimate_from_ssr(histograms, readout_povm, channel_ops, settings):
        """
        Estimate a density matrix from single shot histograms obtained by measuring bitstrings in
        the Z-eigenbasis after application of given channel operators.

        :param numpy.ndarray histograms: The single shot histograms, `shape=(n_channels, dim)`.
        :param DiagognalPOVM readout_povm: The POVM corresponding to the readout plus classifier.
        :param list channel_ops: The tomography measurement channels as `qutip.Qobj`'s.
        :param TomographySettings settings: The solver and estimation settings.
        :return: The generated StateTomography object.
        :rtype: StateTomography
        """
        nqc = len(channel_ops[0].dims[0])
        pauli_basis = grove.tomography.operator_utils.PAULI_BASIS ** nqc
        pi_basis = readout_povm.pi_basis

        if not histograms.shape[1] == pi_basis.dim:  # pragma no coverage
            raise ValueError("Currently tomography is only implemented for two-level systems.")

        # prepare the log-likelihood function parameters, see documentation
        n_kj = np.asarray(histograms)
        c_jk_m = _prepare_c_jk_m(readout_povm, pauli_basis, channel_ops)
        rho_m = cvxpy.Variable(pauli_basis.dim)

        p_jk = c_jk_m * rho_m
        obj = -n_kj.ravel() * cvxpy.log(p_jk)

        p_jk_mat = cvxpy.reshape(p_jk, pi_basis.dim, len(channel_ops))  # cvxpy has col-major order

        # Default constraints:
        # MLE must describe valid probability distribution
        # i.e., for each k, p_jk must sum to one and be element-wise non-negative:
        # 1. \sum_j p_jk == 1  for all k
        # 2. p_jk >= 0         for all j, k
        # where p_jk = \sum_m c_jk_m rho_m
        constraints = [
            p_jk >= 0,
            np.matrix(np.ones((1, pi_basis.dim))) * p_jk_mat == 1,
        ]

        rho_m_real_imag = sum((rm * o_ut.to_realimag(Pm)
                               for (rm, Pm) in ut.izip(rho_m, pauli_basis.ops)), 0)

        if POSITIVE in settings.constraints:
            if tomography._SDP_SOLVER.is_functional():
                constraints.append(rho_m_real_imag >> 0)
            else:  # pragma no coverage
                _log.warning("No convex solver capable of semi-definite problems installed.\n"
                             "Dropping the positivity constraint on the density matrix.")

        if UNIT_TRACE in settings.constraints:
            # this assumes that the first element of the Pauli basis is always proportional to
            # the identity
            constraints.append(rho_m[0, 0] == 1. / pauli_basis.ops[0].tr().real)

        prob = cvxpy.Problem(cvxpy.Minimize(obj), constraints)

        _log.info("Starting convex solver")
        prob.solve(solver=tomography.SOLVER, **settings.solver_kwargs)
        if prob.status != cvxpy.OPTIMAL:  # pragma no coverage
            _log.warning("Problem did not converge to optimal solution. "
                         "Solver settings: {}".format(settings.solver_kwargs))

        return StateTomography(np.array(rho_m.value).ravel(), pauli_basis, settings)
Estimate a density matrix from single shot histograms obtained by measuring bitstrings in
        the Z-eigenbasis after application of given channel operators.

        :param numpy.ndarray histograms: The single shot histograms, `shape=(n_channels, dim)`.
        :param DiagognalPOVM readout_povm: The POVM corresponding to the readout plus classifier.
        :param list channel_ops: The tomography measurement channels as `qutip.Qobj`'s.
        :param TomographySettings settings: The solver and estimation settings.
        :return: The generated StateTomography object.
        :rtype: StateTomography