def run(cosmology, zi=0, Mi=1e12, z=False, com=True, mah=True,
        filename=None, verbose=None, retcosmo=None):
    """ Run commah code on halo of mass 'Mi' at redshift 'zi' with
        accretion and profile history at higher redshifts 'z'
        This is based on Correa et al. (2015a,b,c)

    Parameters
    ----------
    cosmology : str or dict
        Can be named cosmology, default WMAP7 (aka DRAGONS), or
        DRAGONS, WMAP1, WMAP3, WMAP5, WMAP7, WMAP9, Planck13, Planck15
        or dictionary similar in format to:
        {'N_nu': 0,'Y_He': 0.24, 'h': 0.702, 'n': 0.963,'omega_M_0': 0.275,
         'omega_b_0': 0.0458,'omega_lambda_0': 0.725,'omega_n_0': 0.0,
         'sigma_8': 0.816, 't_0': 13.76, 'tau': 0.088,'z_reion': 10.6}

    zi : float / numpy array, optional
        Redshift at which halo has mass 'Mi'. If float then all
        halo masses 'Mi' are assumed to be at this redshift.
        If array but Mi is float, then this halo mass is used across
        all starting redshifts. If both Mi and zi are arrays then they
        have to be the same size for one - to - one correspondence between
        halo mass and the redshift at which it has that mass. Default is 0.
    Mi : float / numpy array, optional
        Halo mass 'Mi' at a redshift 'zi'. If float then all redshifts 'zi'
        are solved for this halo mass. If array but zi is float, then this
        redshift is applied to all halo masses. If both Mi and zi are
        arrays then they have to be the same size for one - to - one
        correspondence between halo mass and the redshift at which it
        has that mass. Default is 1e12 Msol.
    z : float / numpy array, optional
        Redshift to solve commah code at. Must have zi<z else these steps
        are skipped. Default is False, meaning commah is solved at z=zi

    com : bool, optional
        If true then solve for concentration-mass,
        default is True.
    mah : bool, optional
        If true then solve for accretion rate and halo mass history,
        default is True.
    filename : bool / str, optional
        If str is passed this is used as a filename for output of commah
    verbose : bool, optional
        If true then give comments, default is None.
    retcosmo : bool, optional
        Return cosmological parameters used as a dict if retcosmo = True,
        default is None.

    Returns
    -------
    dataset : structured dataset
        dataset contains structured columns of size
        (size(Mi) > size(z)) by size(z)

        If mah = True and com = False then columns are
        ('zi',float),('Mi',float),('z',float),('dMdt',float),('Mz',float)
        where 'zi' is the starting redshift, 'Mi' is halo mass at zi
        'z' is output redshift (NB z>zi), 'dMdt' is accretion rate [Msol/yr]
        and 'Mz' is the halo mass at 'z' for a halo which was 'Mi' massive
        at starting redshift 'zi'

        If mah = False and com = True then columns are
        ('zi',float),('Mi',float),('z',float),('c',float),('sig',float),('nu',float),('zf',float)
        where 'zi' is the starting redshift, 'Mi' is halo mass at zi
        'z' is output redshift (NB z>zi), 'c' is NFW concentration of halo
        at the redshift 'z', 'sig' is the mass variance 'sigma',
        'nu' is the dimensionless fluctuation for halo mass 'Mi' at 'zi',
        'zf' is the formation redshift for a halo of mass 'Mi' at redshift 'zi'

        If mah = True and com = True then columns are:
        ('zi',float),('Mi',float),('z',float),('dMdt',float),('Mz',float),
        ('c',float),('sig',float),('nu',float),('zf',float)

    file : structured dataset with name 'filename' if passed

    Raises
    ------
    Output -1
        If com = False and mah = False as user has to select something.
    Output -1
        If 'zi' and 'Mi' are arrays of unequal size. Impossible to match
        corresponding masses and redshifts of output.

    Examples
    --------
    Examples should be written in doctest format, and should illustrate how
    to use the function.

    >>> import examples
    >>> examples.runcommands() #  A series of ways to query structured dataset
    >>> examples.plotcommands() #  Examples to plot data

    """

    # Check user choices...
    if not com and not mah:
        print("User has to choose com=True and / or mah=True ")
        return(-1)

    # Convert arrays / lists to np.array
    # and inflate redshift / mass axis
    # to match each other for later loop
    results = _checkinput(zi, Mi, z=z, verbose=verbose)

    # Return if results is -1
    if(results == -1):
        return(-1)
    # If not, unpack the returned iterable
    else:
        zi, Mi, z, lenz, lenm, lenzout = results
    # At this point we will have lenm objects to iterate over

    # Get the cosmological parameters for the given cosmology
    cosmo = getcosmo(cosmology)

    # Create  output file if desired
    if filename:
        print("Output to file %r" % (filename))
        fout = open(filename, 'wb')

    # Create the structured dataset
    try:
        if mah and com:
            if verbose:
                print("Output requested is zi, Mi, z, dMdt, Mz, c, sig, nu, "
                      "zf")
            if filename:
                fout.write(_getcosmoheader(cosmo)+'\n')
                fout.write("# Initial z - Initial Halo  - Output z - "
                           " Accretion -  Final Halo  - concentration - "
                           "    Mass   -    Peak    -  Formation z "+'\n')
                fout.write("#           -     mass      -          -"
                           "    rate    -     mass     -               - "
                           " Variance  -   Height   -              "+'\n')
                fout.write("#           -    (M200)     -          - "
                           "  (dM/dt)  -    (M200)    -               - "
                           "  (sigma)  -    (nu)    -              "+'\n')
                fout.write("#           -    [Msol]     -          - "
                           " [Msol/yr] -    [Msol]    -               - "
                           "           -            -              "+'\n')
            dataset = np.zeros((lenm, lenzout), dtype=[('zi', float),
                               ('Mi', float), ('z', float), ('dMdt', float),
                               ('Mz', float), ('c', float), ('sig', float),
                               ('nu', float), ('zf', float)])
        elif mah:
            if verbose:
                print("Output requested is zi, Mi, z, dMdt, Mz")
            if filename:
                fout.write(_getcosmoheader(cosmo)+'\n')
                fout.write("# Initial z - Initial Halo  - Output z -"
                           "   Accretion - Final Halo "+'\n')
                fout.write("#           -     mass      -          -"
                           "     rate    -   mass     "+'\n')
                fout.write("#           -    (M200)     -          -"
                           "    (dm/dt)  -  (M200)    "+'\n')
                fout.write("#           -    [Msol]     -          -"
                           "   [Msol/yr] -  [Msol]    "+'\n')
            dataset = np.zeros((lenm, lenzout), dtype=[('zi', float),
                               ('Mi', float), ('z', float),
                               ('dMdt', float), ('Mz', float)])
        else:
            if verbose:
                print("Output requested is zi, Mi, z, c, sig, nu, zf")
            if filename:
                fout.write(_getcosmoheader(cosmo)+'\n')
                fout.write("# Initial z - Initial Halo  - Output z - "
                           " concentration - "
                           "  Mass    -    Peak    -  Formation z "+'\n')
                fout.write("#           -     mass      -          -"
                           "                -"
                           " Variance  -   Height   -              "+'\n')
                fout.write("#           -   (M200)      -          - "
                           "               - "
                           " (sigma)  -    (nu)    -              "+'\n')
                fout.write("#           -   [Msol]      -          - "
                           "               - "
                           "          -            -            "+'\n')
            dataset = np.zeros((lenm, lenzout), dtype=[('zi', float),
                               ('Mi', float), ('z', float), ('c', float),
                               ('sig', float), ('nu', float), ('zf', float)])

        # Now loop over the combination of initial redshift and halo mamss
        for i_ind, (zval, Mval) in enumerate(_izip(zi, Mi)):
            if verbose:
                print("Output Halo of Mass Mi=%s at zi=%s" % (Mval, zval))
            # For a given halo mass Mi at redshift zi need to know
            # output redshifts 'z'
            # Check that all requested redshifts are greater than
            # input redshift, except if z is False, in which case
            # only solve z at zi, i.e. remove a loop
            if z is False:
                ztemp = np.array(zval, ndmin=1, dtype=float)
            else:
                ztemp = np.array(z[z >= zval], dtype=float)

            # Loop over the output redshifts
            if ztemp.size:
                # Return accretion rates and halo mass progenitors at
                # redshifts 'z' for object of mass Mi at zi
                dMdt, Mz = MAH(ztemp, zval, Mval, **cosmo)
                if mah and com:
                    # More expensive to return concentrations
                    c, sig, nu, zf = COM(ztemp, Mz, **cosmo)
                    # Save all arrays
                    for j_ind, j_val in enumerate(ztemp):
                        dataset[i_ind, j_ind] =\
                            (zval, Mval, ztemp[j_ind], dMdt[j_ind], Mz[j_ind],
                             c[j_ind], sig[j_ind], nu[j_ind], zf[j_ind])
                        if filename:
                            fout.write(
                                "{}, {}, {}, {}, {}, {}, {}, {}, {} \n".format(
                                    zval, Mval, ztemp[j_ind], dMdt[j_ind],
                                    Mz[j_ind], c[j_ind], sig[j_ind], nu[j_ind],
                                    zf[j_ind]))
                elif mah:
                    # Save only MAH arrays
                    for j_ind, j_val in enumerate(ztemp):
                        dataset[i_ind, j_ind] =\
                            (zval, Mval, ztemp[j_ind], dMdt[j_ind], Mz[j_ind])
                        if filename:
                            fout.write("{}, {}, {}, {}, {} \n".format(
                                zval, Mval, ztemp[j_ind], dMdt[j_ind],
                                Mz[j_ind]))
                else:
                    # Output only COM arrays
                    c, sig, nu, zf = COM(ztemp, Mz, **cosmo)
                    # For any halo mass Mi at redshift zi
                    # solve for c, sig, nu and zf
                    for j_ind, j_val in enumerate(ztemp):
                        dataset[i_ind, j_ind] =\
                            (zval, Mval, ztemp[j_ind], c[j_ind], sig[j_ind],
                             nu[j_ind], zf[j_ind])
                        if filename:
                            fout.write("{}, {}, {}, {}, {}, {}, {} \n".format(
                                zval, Mval, ztemp[j_ind], c[j_ind], sig[j_ind],
                                nu[j_ind], zf[j_ind]))

    # Make sure to close the file if it was opened
    finally:
        fout.close() if filename else None

    if retcosmo:
        return(dataset, cosmo)
    else:
        return(dataset)
Run commah code on halo of mass 'Mi' at redshift 'zi' with
        accretion and profile history at higher redshifts 'z'
        This is based on Correa et al. (2015a,b,c)

    Parameters
    ----------
    cosmology : str or dict
        Can be named cosmology, default WMAP7 (aka DRAGONS), or
        DRAGONS, WMAP1, WMAP3, WMAP5, WMAP7, WMAP9, Planck13, Planck15
        or dictionary similar in format to:
        {'N_nu': 0,'Y_He': 0.24, 'h': 0.702, 'n': 0.963,'omega_M_0': 0.275,
         'omega_b_0': 0.0458,'omega_lambda_0': 0.725,'omega_n_0': 0.0,
         'sigma_8': 0.816, 't_0': 13.76, 'tau': 0.088,'z_reion': 10.6}

    zi : float / numpy array, optional
        Redshift at which halo has mass 'Mi'. If float then all
        halo masses 'Mi' are assumed to be at this redshift.
        If array but Mi is float, then this halo mass is used across
        all starting redshifts. If both Mi and zi are arrays then they
        have to be the same size for one - to - one correspondence between
        halo mass and the redshift at which it has that mass. Default is 0.
    Mi : float / numpy array, optional
        Halo mass 'Mi' at a redshift 'zi'. If float then all redshifts 'zi'
        are solved for this halo mass. If array but zi is float, then this
        redshift is applied to all halo masses. If both Mi and zi are
        arrays then they have to be the same size for one - to - one
        correspondence between halo mass and the redshift at which it
        has that mass. Default is 1e12 Msol.
    z : float / numpy array, optional
        Redshift to solve commah code at. Must have zi<z else these steps
        are skipped. Default is False, meaning commah is solved at z=zi

    com : bool, optional
        If true then solve for concentration-mass,
        default is True.
    mah : bool, optional
        If true then solve for accretion rate and halo mass history,
        default is True.
    filename : bool / str, optional
        If str is passed this is used as a filename for output of commah
    verbose : bool, optional
        If true then give comments, default is None.
    retcosmo : bool, optional
        Return cosmological parameters used as a dict if retcosmo = True,
        default is None.

    Returns
    -------
    dataset : structured dataset
        dataset contains structured columns of size
        (size(Mi) > size(z)) by size(z)

        If mah = True and com = False then columns are
        ('zi',float),('Mi',float),('z',float),('dMdt',float),('Mz',float)
        where 'zi' is the starting redshift, 'Mi' is halo mass at zi
        'z' is output redshift (NB z>zi), 'dMdt' is accretion rate [Msol/yr]
        and 'Mz' is the halo mass at 'z' for a halo which was 'Mi' massive
        at starting redshift 'zi'

        If mah = False and com = True then columns are
        ('zi',float),('Mi',float),('z',float),('c',float),('sig',float),('nu',float),('zf',float)
        where 'zi' is the starting redshift, 'Mi' is halo mass at zi
        'z' is output redshift (NB z>zi), 'c' is NFW concentration of halo
        at the redshift 'z', 'sig' is the mass variance 'sigma',
        'nu' is the dimensionless fluctuation for halo mass 'Mi' at 'zi',
        'zf' is the formation redshift for a halo of mass 'Mi' at redshift 'zi'

        If mah = True and com = True then columns are:
        ('zi',float),('Mi',float),('z',float),('dMdt',float),('Mz',float),
        ('c',float),('sig',float),('nu',float),('zf',float)

    file : structured dataset with name 'filename' if passed

    Raises
    ------
    Output -1
        If com = False and mah = False as user has to select something.
    Output -1
        If 'zi' and 'Mi' are arrays of unequal size. Impossible to match
        corresponding masses and redshifts of output.

    Examples
    --------
    Examples should be written in doctest format, and should illustrate how
    to use the function.

    >>> import examples
    >>> examples.runcommands() #  A series of ways to query structured dataset
    >>> examples.plotcommands() #  Examples to plot data