def read_double_out(self, out_file):
        with open(out_file, 'rb') as f:
            out_arr = np.fromstring(f.read()[13:], dtype='>d').reshape(
                (-1, self.nt + self.nr + 1)).T
            out_arr = np.asarray(out_arr, dtype=np.float64).T

        # Sort by the absolute value of the stoichiometry
        sort_inds= np.abs(out_arr[:, :self.nr]).sum(1).argsort()
        out_arr = out_arr[sort_inds]

        unbounded = out_arr[np.isclose(out_arr[:,-1], 0.)]
        bounded = out_arr[~np.isclose(out_arr[:,-1], 0.)]

        if bounded.size:  # Test if its empty
            bounded /= np.atleast_2d(bounded[:,-1]).T

        unbounded_df = pd.DataFrame(
            unbounded[:, :self.nr], 
            columns=[r.id for r in self.model.reactions],
            index=['UEV{}'.format(i) 
                   for i in range(1, unbounded.shape[0] + 1)])

        bounded_df = pd.DataFrame(
            bounded[:, :self.nr], 
            columns=[r.id for r in self.model.reactions],
            index=('BEV{}'.format(i) 
                   for i in range(1, bounded.shape[0] + 1)))

        return unbounded_df.append(bounded_df)