%prog simulation inversion.txt translocation.txt maps.txt multimaps.txt

    Plot ALLMAPS accuracy across a range of simulated datasets.
767623
def simulation(args):
    
    p = OptionParser(simulation.__doc__)
    opts, args, iopts = p.set_image_options(args, dpi=300)

    if len(args) != 4:
        sys.exit(not p.print_help())

    dataA, dataB, dataC, dataD = args
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])
    A = fig.add_axes([.12, .62, .35, .35])
    B = fig.add_axes([.62, .62, .35, .35])
    C = fig.add_axes([.12, .12, .35, .35])
    D = fig.add_axes([.62, .12, .35, .35])
    dataA = import_data(dataA)
    dataB = import_data(dataB)
    dataC = import_data(dataC)
    dataD = import_data(dataD)
    subplot(A, dataA, "Inversion error rate", "Accuracy", xlim=.5)
    subplot(B, dataB, "Translocation error rate", "Accuracy", xlim=.5,
                      legend=("intra-chromosomal", "inter-chromosomal",
                              "75\% intra + 25\% inter"))
    subplot(C, dataC, "Number of input maps", "Accuracy", xcast=int)
    subplot(D, dataD, "Number of input maps", "Accuracy", xcast=int)

    labels = ((.03, .97, "A"), (.53, .97, "B"),
              (.03, .47, "C"), (.53, .47, "D"))
    panel_labels(root, labels)

    normalize_axes(root)
    image_name = "simulation." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1226302
def run_multisample_qualimap(output_dir, work_dir, samples, targqc_full_report):
    
    plots_dirpath = join(output_dir, 'plots')
    individual_report_fpaths = [s.qualimap_html_fpath for s in samples]
    if isdir(plots_dirpath) and not any(
            not can_reuse(join(plots_dirpath, f), individual_report_fpaths)
            for f in listdir(plots_dirpath) if not f.startswith('.')):
        debug('Qualimap miltisample plots exist - ' + plots_dirpath + ', reusing...')
    else:
        # Qualimap2 run for multi-sample plots
        if len([s.qualimap_html_fpath for s in samples if s.qualimap_html_fpath]) > 0:
            if find_executable() is not None:  # and get_qualimap_type(find_executable()) == 'full':
                qualimap_output_dir = join(work_dir, 'qualimap_multi_bamqc')

                _correct_qualimap_genome_results(samples)
                _correct_qualimap_insert_size_histogram(samples)

                safe_mkdir(qualimap_output_dir)
                rows = []
                for sample in samples:
                    if sample.qualimap_html_fpath:
                        rows += [[sample.name, sample.qualimap_html_fpath]]

                data_fpath = write_tsv_rows(([], rows), join(qualimap_output_dir, 'qualimap_results_by_sample.tsv'))
                qualimap_plots_dirpath = join(qualimap_output_dir, 'images_multisampleBamQcReport')
                cmdline = find_executable() + ' multi-bamqc --data {data_fpath} -outdir {qualimap_output_dir}'.format(**locals())
                run(cmdline, env_vars=dict(DISPLAY=None),
                    checks=[lambda _1, _2: verify_dir(qualimap_output_dir)], reuse=cfg.reuse_intermediate)

                if not verify_dir(qualimap_plots_dirpath):
                    warn('Warning: Qualimap for multi-sample analysis failed to finish. TargQC will not contain plots.')
                    return None
                else:
                    if exists(plots_dirpath):
                        shutil.rmtree(plots_dirpath)
                    shutil.move(qualimap_plots_dirpath, plots_dirpath)
            else:
                warn('Warning: Qualimap for multi-sample analysis was not found. TargQC will not contain plots.')
                return None

    targqc_full_report.plots = []
    for plot_fpath in listdir(plots_dirpath):
        plot_fpath = join(plots_dirpath, plot_fpath)
        if verify_file(plot_fpath) and plot_fpath.endswith('.png'):
            targqc_full_report.plots.append(relpath(plot_fpath, output_dir))
332271
def _map_arguments(self, input_dict):
        
        data = input_dict.get('data')
        comp = input_dict.get('comp')
        library = input_dict.get('library')
        models = input_dict.get('models')
        hpx_order = input_dict.get('hpx_order_fitting')
        dry_run = input_dict.get('dry_run', False)

        self._set_link('init-model', InitModel,
                       comp=comp, data=data,
                       library=library,
                       models=models,
                       hpx_order=hpx_order,
                       dry_run=dry_run)

        self._set_link('assemble-model', AssembleModel_SG,
                       comp=comp, data=data,
                       models=models)
459956
def main(): # mini/unit test
    
    options = optionsParser().parse_args()
    params = getSelectionParams(options)
    if options.list or options.details:
        specifiedMaps = filterMapNames(
            options.mapname,
            records = filterMapAttrs(**params),
            excludeRegex=options.exclude,
            closestMatch=options.best
        )
        if specifiedMaps:
            for v in specifiedMaps:
                if options.details: v.display()
                else:               print(v)
            print("Found %d maps that match given criteria."%(len(specifiedMaps)))
        else:
            print("No matching maps found.")
    else:
        try:
            specifiedMaps = selectMap(
                options.mapname,
                excludeName =options.exclude,
                closestMatch=options.best,
                **params)
        except Exception as e:
            specifiedMaps = []
            print("No matching maps found: %s"%e)
        if not isinstance(specifiedMaps, list):
            specifiedMaps = [specifiedMaps]
        for m in specifiedMaps:
            if options.path:    print(m.path)
            else:               print(m.name)
706485
def plot_accuracy(data, output_dir_path='.', output_filename='accuracy.png',
                  width=10, height=8):
    
    output_path = os.path.join(output_dir_path, output_filename)

    max_val_data = get_epoch_max_val_acc(data)
    max_val_label = round(max_val_data['acc'].values[0], 4)

    # max_val_epoch = max_val_data['epoch'].values[0]
    max_epoch_data = data[data['epoch'] == data['epoch'].max()]

    plot = ggplot(data, aes('epoch', 'acc', color='factor(data)')) + \
        geom_line(size=1, show_legend=False) + \
        geom_vline(aes(xintercept='epoch', color='data'),
                   data=max_val_data, alpha=0.5, show_legend=False) + \
        geom_label(aes('epoch', 'acc'), data=max_val_data,
                   label=max_val_label, nudge_y=-0.02, va='top', label_size=0,
                   show_legend=False) + \
        geom_text(aes('epoch', 'acc', label='data'), data=max_epoch_data,
                  nudge_x=2, ha='center', show_legend=False) + \
        geom_point(aes('epoch', 'acc'), data=max_val_data,
                   show_legend=False) + \
        labs(y='Accuracy', x='Epochs') + \
        theme_bw(base_family='Arial', base_size=15) + \
        scale_color_manual(['#ef8a62', '#67a9cf', "#f7f7f7"])

    plot.save(output_path, width=width, height=height)
672819
def _map_arguments(self, args):
        
        comp_file = args.get('comp', None)
        datafile = args.get('data', None)
        do_ltsum = args.get('do_ltsum', False)

        NAME_FACTORY.update_base_dict(datafile)
        outdir_base = os.path.join(NAME_FACTORY.base_dict['basedir'], 'counts_cubes')
        num_files = args.get('nfiles', 96)

        self.comp_dict = yaml.safe_load(open(comp_file))
        coordsys = self.comp_dict.pop('coordsys')

        for key_e, comp_e in sorted(self.comp_dict.items()):

            if 'mktimefilters' in comp_e:
                mktimelist = comp_e['mktimefilters']
            else:
                mktimelist = ['none']

            if 'evtclasses' in comp_e:
                evtclasslist_vals = comp_e['evtclasses']
            else:
                evtclasslist_vals = [NAME_FACTORY.base_dict['evclass']]

            for mktimekey in mktimelist:
                zcut = "zmax%i" % comp_e['zmax']
                kwargs_mktime = dict(zcut=zcut,
                                     ebin=key_e,
                                     psftype='ALL',
                                     coordsys=coordsys,
                                     mktime=mktimekey)

                if do_ltsum:
                    ltsum_listfile = 'ltsumlist_%s_%s' % (key_e, mktimekey)
                    ltsum_outfile = 'ltsum_%s_%s' % (key_e, mktimekey)
                    linkname = 'ltsum_%s_%s' % (key_e, mktimekey)
                    self._set_link(likname, Gtlink_ltsum,
                                   infile1=ltsum_listfile,
                                   infile2=None,
                                   outfile=ltsum_outfile,
                                   logfile=os.path.join(outdir_base, "%s.log" % linkname))

                for evtclassval in evtclasslist_vals:
                    for psf_type in sorted(comp_e['psf_types'].keys()):
                        fullkey = "%s_%s_%s_%s"%(key_e, mktimekey, evtclassval, psf_type)
                        linkname = 'coadd_%s' % (fullkey)
                        kwargs_bin = kwargs_mktime.copy()
                        kwargs_bin['psftype'] = psf_type
                        kwargs_bin['evclass'] = evtclassval
                        ccube_name =\
                            os.path.basename(NAME_FACTORY.ccube(**kwargs_bin))
                        outputfile = os.path.join(outdir_base, ccube_name)
                        args = _make_input_file_list(outputfile, num_files)
                        
                        self._set_link(linkname,
                                       Link_FermipyCoadd,
                                       args=args,
                                       output=outputfile,
                                       logfile=os.path.join(outdir_base, "%s.log" % linkname))
345955
def plot(args):
    
    from jcvi.graphics.base import savefig

    p = OptionParser(plot.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="8x7", format="png")

    if len(args) != 3:
        sys.exit(not p.print_help())

    workdir, sample_key, chrs = args
    chrs = chrs.split(",")
    hmm = CopyNumberHMM(workdir=workdir)
    hmm.plot(sample_key, chrs=chrs)

    image_name = sample_key + "_cn." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
333814
def simulate(args):
    """"""
    if args.fasta:
        name = None
        seq = ""
        reads = dict()
        with open(args.fasta) as in_handle:
            for line in in_handle:
                if line.startswith(">"):
                    if name:
                        reads.update(_generate_reads(seq, name))
                    seq = ""
                    name = line[1:-1]
                else:
                    seq += line.strip()

        reads.update(_generate_reads(seq, name))
    _write_reads(reads, args.out)
1050372
def _print_map_dict(self, argkey, filename, append):
        """"""
        result = []
        skeys = list(sorted(self.curargs[argkey].keys()))
        for key in skeys:
            result.append("'{}' => {}".format(key, self.curargs[argkey][key]))
        self._redirect_output('\n'.join(result), filename, append, msg.info)
439103
def thconfig(txtfiles, cavename='cave'):
    """"""
    fmts_model = ['lox', '3d', 'dxf', 'kml', 'plt', 'vrml']
    outfilename = 'thconfig'

    with open(outfilename, 'w') as outfile:
        for txtfilename in txtfiles:
            print >> outfile, 'source "%s"' % convert_filename(txtfilename)
        print >> outfile

        for fmt in fmts_model:
            print >> outfile, 'export model -output "%s"' % os.path.join('models', '%s.%s' % (cavename, fmt))

        print >> outfile, 'export map -projection plan -format esri -output "%s"' % os.path.join('models', 'gis')

        print 'Wrote Therion project file %s .' % outfilename
467051
def cli_action(args):
    
    problem = read_param_file(args.paramfile)
    param_values = sample(problem, seed=args.seed)
    np.savetxt(args.output, param_values, delimiter=args.delimiter,
               fmt='%.' + str(args.precision) + 'e')
933185
def main(args=None):

    # parse command-line options
    parser = argparse.ArgumentParser(prog='rect_wpoly_for_mos')
    # required arguments
    parser.add_argument("input_list",
                        help="TXT file with list JSON files derived from "
                             "longslit data")
    parser.add_argument("--fitted_bound_param", required=True,
                        help="Input JSON with fitted boundary parameters",
                        type=argparse.FileType('rt'))
    parser.add_argument("--out_MOSlibrary", required=True,
                        help="Output JSON file with results",
                        type=lambda x: arg_file_is_new(parser, x))
    # optional arguments
    parser.add_argument("--debugplot",
                        help="Integer indicating plotting & debugging options"
                             " (default=0)",
                        default=0, type=int,
                        choices=DEBUGPLOT_CODES)
    parser.add_argument("--echo",
                        help="Display full command line",
                        action="store_true")
    args = parser.parse_args(args)

    if args.echo:
        print('\033[1m\033[31m% ' + ' '.join(sys.argv) + '\033[0m\n')

    # ---

    # Read input TXT file with list of JSON files
    list_json_files = list_fileinfo_from_txt(args.input_list)
    nfiles = len(list_json_files)
    if abs(args.debugplot) >= 10:
        print('>>> Number of input JSON files:', nfiles)
        for item in list_json_files:
            print(item)
    if nfiles < 2:
        raise ValueError("Insufficient number of input JSON files")

    # read fitted boundary parameters and check that all the longslit JSON
    # files have been computed using the same fitted boundary parameters
    refined_boundary_model = RefinedBoundaryModelParam._datatype_load(
        args.fitted_bound_param.name)
    for ifile in range(nfiles):
        coef_rect_wpoly = RectWaveCoeff._datatype_load(
            list_json_files[ifile].filename)
        uuid_tmp = coef_rect_wpoly.meta_info['origin']['bound_param']
        if uuid_tmp[4:] != refined_boundary_model.uuid:
            print('Expected uuid:', refined_boundary_model.uuid)
            print('uuid for ifile #' + str(ifile + 1) + ": " + uuid_tmp)
            raise ValueError("Fitted boundary parameter uuid's do not match")

    # check consistency of grism, filter, DTU configuration and list of
    # valid slitlets
    coef_rect_wpoly_first_longslit = RectWaveCoeff._datatype_load(
        list_json_files[0].filename)
    filter_name = coef_rect_wpoly_first_longslit.tags['filter']
    grism_name = coef_rect_wpoly_first_longslit.tags['grism']
    dtu_conf = DtuConfiguration.define_from_dictionary(
        coef_rect_wpoly_first_longslit.meta_info['dtu_configuration']
    )
    list_valid_islitlets = list(range(1, EMIR_NBARS + 1))
    for idel in coef_rect_wpoly_first_longslit.missing_slitlets:
        list_valid_islitlets.remove(idel)
    for ifile in range(1, nfiles):
        coef_rect_wpoly = RectWaveCoeff._datatype_load(
            list_json_files[ifile].filename)
        filter_tmp = coef_rect_wpoly.tags['filter']
        if filter_name != filter_tmp:
            print(filter_name)
            print(filter_tmp)
            raise ValueError("Unexpected different filter found")
        grism_tmp = coef_rect_wpoly.tags['grism']
        if grism_name != grism_tmp:
            print(grism_name)
            print(grism_tmp)
            raise ValueError("Unexpected different grism found")
        coef_rect_wpoly = RectWaveCoeff._datatype_load(
            list_json_files[ifile].filename)
        dtu_conf_tmp = DtuConfiguration.define_from_dictionary(
            coef_rect_wpoly.meta_info['dtu_configuration']
        )
        if dtu_conf != dtu_conf_tmp:
            print(dtu_conf)
            print(dtu_conf_tmp)
            raise ValueError("Unexpected different DTU configurations found")
        list_valid_islitlets_tmp = list(range(1, EMIR_NBARS + 1))
        for idel in coef_rect_wpoly.missing_slitlets:
            list_valid_islitlets_tmp.remove(idel)
        if list_valid_islitlets != list_valid_islitlets_tmp:
            print(list_valid_islitlets)
            print(list_valid_islitlets_tmp)
            raise ValueError("Unexpected different list of valid slitlets")

    # check consistency of horizontal bounding box limits (bb_nc1_orig and
    # bb_nc2_orig) and ymargin_bb, and store the values for each slitlet
    dict_bb_param = {}
    print("Checking horizontal bounding box limits and ymargin_bb:")
    for islitlet in list_valid_islitlets:
        islitlet_progress(islitlet, EMIR_NBARS)
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        dict_bb_param[cslitlet] = {}
        for par in ['bb_nc1_orig', 'bb_nc2_orig', 'ymargin_bb']:
            value_initial = \
                coef_rect_wpoly_first_longslit.contents[islitlet - 1][par]
            for ifile in range(1, nfiles):
                coef_rect_wpoly = RectWaveCoeff._datatype_load(
                    list_json_files[ifile].filename)
                value_tmp = coef_rect_wpoly.contents[islitlet - 1][par]
                if value_initial != value_tmp:
                    print(islitlet, value_initial, value_tmp)
                    print(value_tmp)
                    raise ValueError("Unexpected different " + par)
                dict_bb_param[cslitlet][par] = value_initial
    print('OK!')

    # ---

    # Read and store all the longslit data
    list_coef_rect_wpoly = []
    for ifile in range(nfiles):
        coef_rect_wpoly = RectWaveCoeff._datatype_load(
            list_json_files[ifile].filename)
        list_coef_rect_wpoly.append(coef_rect_wpoly)

    # ---

    # Initialize structure to save results into an ouptut JSON file
    outdict = {}
    outdict['refined_boundary_model'] = refined_boundary_model.__getstate__()
    outdict['instrument'] = 'EMIR'
    outdict['meta_info'] = {}
    outdict['meta_info']['creation_date'] = datetime.now().isoformat()
    outdict['meta_info']['description'] = \
        'rectification and wavelength calibration polynomial coefficients ' \
        'as a function of csu_bar_slit_center for MOS'
    outdict['meta_info']['recipe_name'] = 'undefined'
    outdict['meta_info']['origin'] = {}
    outdict['meta_info']['origin']['wpoly_longslits'] = {}
    for ifile in range(nfiles):
        cdum = 'longslit_' + str(ifile + 1).zfill(3) + '_uuid'
        outdict['meta_info']['origin']['wpoly_longslits'][cdum] = \
            list_coef_rect_wpoly[ifile].uuid
    outdict['tags'] = {}
    outdict['tags']['grism'] = grism_name
    outdict['tags']['filter'] = filter_name
    outdict['dtu_configuration'] = dtu_conf.outdict()
    outdict['uuid'] = str(uuid4())
    outdict['contents'] = {}

    # include bb_nc1_orig, bb_nc2_orig and ymargin_bb for each slitlet
    # (note that the values of bb_ns1_orig and bb_ns2_orig cannot be
    # computed at this stage because they depend on csu_bar_slit_center)
    for islitlet in list_valid_islitlets:
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        outdict['contents'][cslitlet] = dict_bb_param[cslitlet]

    # check that order for rectification transformations is the same for all
    # the slitlets and longslit configurations
    order_check_list = []
    for ifile in range(nfiles):
        tmpdict = list_coef_rect_wpoly[ifile].contents
        for islitlet in list_valid_islitlets:
            ttd_order = tmpdict[islitlet - 1]['ttd_order']
            if ttd_order is not None:
                order_check_list.append(ttd_order)
            ttd_order_modeled = \
                tmpdict[islitlet - 1]['ttd_order_longslit_model']
            order_check_list.append(ttd_order_modeled)
    # remove duplicates in list
    order_no_duplicates = list(set(order_check_list))
    if len(order_no_duplicates) != 1:
        print('order_no_duplicates:', order_no_duplicates)
        raise ValueError('tdd_order is not constant!')
    ttd_order = int(order_no_duplicates[0])
    ncoef_rect = ncoef_fmap(ttd_order)
    if abs(args.debugplot) >= 10:
        print('>>> ttd_order........:', ttd_order)
        print('>>> ncoef_rect.......:', ncoef_rect)

    # check that polynomial degree in frontiers and spectrails are the same
    poldeg_check_list = []
    for ifile in range(nfiles):
        tmpdict = list_coef_rect_wpoly[ifile].contents
        for islitlet in list_valid_islitlets:
            tmppoly = tmpdict[islitlet - 1]['frontier']['poly_coef_lower']
            poldeg_check_list.append(len(tmppoly) - 1)
            tmppoly = tmpdict[islitlet - 1]['frontier']['poly_coef_upper']
            poldeg_check_list.append(len(tmppoly) - 1)
            tmppoly = tmpdict[islitlet - 1]['spectrail']['poly_coef_lower']
            poldeg_check_list.append(len(tmppoly) - 1)
            tmppoly = tmpdict[islitlet - 1]['spectrail']['poly_coef_middle']
            poldeg_check_list.append(len(tmppoly) - 1)
            tmppoly = tmpdict[islitlet - 1]['spectrail']['poly_coef_upper']
            poldeg_check_list.append(len(tmppoly) - 1)
    # remove duplicates in list
    poldeg_no_duplicates = list(set(poldeg_check_list))
    if len(poldeg_no_duplicates) != 1:
        print('poldeg_no_duplicates:', poldeg_no_duplicates)
        raise ValueError('poldeg is not constant in frontiers and '
                         'spectrails!')
    poldeg_spectrails = int(poldeg_no_duplicates[0])
    if abs(args.debugplot) >= 10:
        print('>>> poldeg spectrails:', poldeg_spectrails)

    # check that polynomial degree of wavelength calibration is the same for
    # all the slitlets
    poldeg_check_list = []
    for ifile in range(nfiles):
        tmpdict = list_coef_rect_wpoly[ifile].contents
        for islitlet in list_valid_islitlets:
            tmppoly = tmpdict[islitlet - 1]['wpoly_coeff']
            poldeg_check_list.append(len(tmppoly) - 1)
            tmppoly = tmpdict[islitlet - 1]['wpoly_coeff_longslit_model']
            poldeg_check_list.append(len(tmppoly) - 1)
    # remove duplicates in list
    poldeg_no_duplicates = list(set(poldeg_check_list))
    if len(poldeg_no_duplicates) != 1:
        print('poldeg_no_duplicates:', poldeg_no_duplicates)
        raise ValueError('poldeg is not constant in wavelength calibration '
                         'polynomials!')
    poldeg_wavecal = int(poldeg_no_duplicates[0])
    if abs(args.debugplot) >= 10:
        print('>>> poldeg wavecal...:', poldeg_wavecal)

    # ---

    # csu_bar_slit_center values for each slitlet
    print("CSU_bar_slit_center values:")
    for islitlet in list_valid_islitlets:
        islitlet_progress(islitlet, EMIR_NBARS)
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        list_csu_bar_slit_center = []
        for ifile in range(nfiles):
            tmpdict = list_coef_rect_wpoly[ifile].contents[islitlet - 1]
            csu_bar_slit_center = tmpdict['csu_bar_slit_center']
            list_csu_bar_slit_center.append(csu_bar_slit_center)
        # check that list_csu_bar_slit_center is properly sorted
        if not np.all(list_csu_bar_slit_center[:-1] <=
                  list_csu_bar_slit_center[1:]):
            print('cslitlet: ', cslitlet)
            print('list_csu_bar_slit_center: ', list_csu_bar_slit_center)
            raise ValueError('Unsorted list_csu_bar_slit_center')
        outdict['contents'][cslitlet]['list_csu_bar_slit_center'] = \
            list_csu_bar_slit_center
    print('OK!')

    # ---

    # rectification polynomial coefficients

    # note: when aij and bij have not been computed, we use the modeled
    # version aij_longslit_model and bij_longslit_model
    print("Rectification polynomial coefficients:")
    for islitlet in list_valid_islitlets:
        islitlet_progress(islitlet, EMIR_NBARS)
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        outdict['contents'][cslitlet]['ttd_order'] = ttd_order
        outdict['contents'][cslitlet]['ncoef_rect'] = ncoef_rect
        for keycoef in ['ttd_aij', 'ttd_bij', 'tti_aij', 'tti_bij']:
            for icoef in range(ncoef_rect):
                ccoef = str(icoef).zfill(2)
                list_cij = []
                for ifile in range(nfiles):
                    tmpdict = \
                        list_coef_rect_wpoly[ifile].contents[islitlet - 1]
                    cij = tmpdict[keycoef]
                    if cij is not None:
                        list_cij.append(cij[icoef])
                    else:
                        cij_modeled = tmpdict[keycoef + '_longslit_model']
                        if cij_modeled is None:
                            raise ValueError("Unexpected cij_modeled=None!")
                        else:
                            list_cij.append(cij_modeled[icoef])
                        if abs(args.debugplot) >= 10:
                            print("Warning: using " + keycoef +
                                  "_longslit_model for " + cslitlet +
                                  " in file " +
                                  list_json_files[ifile].filename)
                cdum = 'list_' + keycoef + '_' + ccoef
                outdict['contents'][cslitlet][cdum] = list_cij
    print('OK!')

    # ---

    # wavelength calibration polynomial coefficients

    # note: when wpoly_coeff have not been computed, we use the
    # wpoly_coeff_longslit_model
    print("Wavelength calibration polynomial coefficients:")
    for islitlet in list_valid_islitlets:
        islitlet_progress(islitlet, EMIR_NBARS)
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        outdict['contents'][cslitlet]['wpoly_degree'] = poldeg_wavecal
        for icoef in range(poldeg_wavecal + 1):
            ccoef = str(icoef).zfill(2)
            list_cij = []
            for ifile in range(nfiles):
                tmpdict = list_coef_rect_wpoly[ifile].contents[islitlet - 1]
                cij = tmpdict['wpoly_coeff']
                if cij is not None:
                    list_cij.append(cij[icoef])
                else:
                    cij_modeled = tmpdict['wpoly_coeff_longslit_model']
                    if cij_modeled is None:
                        raise ValueError("Unexpected cij_modeled=None!")
                    else:
                        list_cij.append(cij_modeled[icoef])
                    if abs(args.debugplot) >= 10:
                        print("Warning: using wpoly_coeff_longslit_model" +
                              " for " + cslitlet +
                              " in file " +
                              list_json_files[ifile].filename)
            outdict['contents'][cslitlet]['list_wpoly_coeff_' + ccoef] = \
                list_cij
    print('OK!')

    # ---

    # OBSOLETE
    # Save resulting JSON structure
    '''
    
    '''

    # --

    # Create object of type MasterRectWave with library of coefficients
    # for rectification and wavelength calibration
    master_rectwv = MasterRectWave(instrument='EMIR')
    master_rectwv.quality_control = numina.types.qc.QC.GOOD
    master_rectwv.tags['grism'] = grism_name
    master_rectwv.tags['filter'] = filter_name
    master_rectwv.meta_info['dtu_configuration'] = outdict['dtu_configuration']
    master_rectwv.meta_info['refined_boundary_model'] = {
        'parmodel': refined_boundary_model.meta_info['parmodel']
    }
    master_rectwv.meta_info['refined_boundary_model'].update(
        outdict['refined_boundary_model']['contents']
    )
    master_rectwv.total_slitlets = EMIR_NBARS
    master_rectwv.meta_info['origin'] = {
        'bound_param': 'uuid' + refined_boundary_model.uuid,
        'longslit_frames': ['uuid:' + list_coef_rect_wpoly[ifile].uuid
                            for ifile in range(nfiles)]
    }
    for i in range(EMIR_NBARS):
        islitlet = i + 1
        dumdict = {'islitlet': islitlet}
        cslitlet = 'slitlet' + str(islitlet).zfill(2)
        if cslitlet in outdict['contents']:
            dumdict.update(outdict['contents'][cslitlet])
        else:
            dumdict.update({
                'bb_nc1_orig': 0,
                'bb_nc2_orig': 0,
                'ymargin_bb': 0,
                'list_csu_bar_slit_center': [],
                'ttd_order': 0,
                'ncoef_rect': 0,
                'wpolydegree': 0
            })
            master_rectwv.missing_slitlets.append(islitlet)
        master_rectwv.contents.append(dumdict)
    master_rectwv.writeto(args.out_MOSlibrary.name)
    print('>>> Saving file ' + args.out_MOSlibrary.name)
80705
def from_mapping(self, *mapping, **kwargs):
        
        mappings = []
        if len(mapping) == 1:
            if hasattr(mapping[0], 'items'):
                mappings.append(mapping[0].items())
            else:
                mappings.append(mapping[0])
        elif len(mapping) > 1:
            raise TypeError('expected at most 1 positional argument, got %d' % len(mapping))

        mappings.append(kwargs.items())
        for mapping in mappings:
            for key, value in mapping:
                if key.isupper():
                    self[key] = value
        return True
114191
def histogram(args):
    
    from jcvi.graphics.histogram import histogram_multiple

    p = OptionParser(histogram.__doc__)
    p.add_option("--bins", dest="bins", default=40, type="int",
            help="number of bins to plot in the histogram [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(not p.print_help())

    gff_files = args
    # metrics = ("Exon_Length", "Intron_Length", "Gene_Length", "Exon_Count")
    colors = ("red", "green", "blue", "black")
    vmaxes = (1000, 1000, 4000, 20)
    xlabels = ("bp", "bp", "bp", "number")
    for metric, color, vmax, xlabel in zip(metrics, colors, vmaxes, xlabels):
        logging.debug("Parsing files in `{0}`..".format(metric))
        numberfiles = [op.join(metric, op.basename(x).split(".")[0] + ".txt") \
                        for x in gff_files]

        histogram_multiple(numberfiles, 0, vmax, xlabel, metric,
                       bins=opts.bins, facet=True, fill=color,
                       prefix=metric + ".")
57761
def plot_results(fout_png, goea_results, *args, **kws):
    """"""
    if "{NS}" not in fout_png:
        plt_goea_results(fout_png, goea_results, *args, **kws)
    else:
        # Plot separately by NS: BP, MF, CC
        ns2goea_results = cx.defaultdict(list)
        for rec in goea_results:
            ns2goea_results[rec.NS].append(rec)
        for ns_name, ns_res in ns2goea_results.items():
            png = fout_png.format(NS=ns_name)
            plt_goea_results(png, ns_res, *args, **kws)
1188031
public static void main(String[] args) throws FileParsingException {
    if (args.length == 0)
      System.err.println("Supply arguments that are paths to pepxml files");
    List<Path> paths = Arrays.stream(args).map(s -> Paths.get(s)).collect(Collectors.toList());
    if (paths.stream().anyMatch(p -> !Files.exists(p)))
      System.err.println("Not all input files exist");

    for (Path p : paths) {
      System.out.println("\nParsing: " + p.toString());
      MsmsPipelineAnalysis pepxml = PepXmlParser.parse(p);
      int sum = pepxml.getMsmsRunSummary().stream()
          .flatMapToInt(msmsRunSummary -> msmsRunSummary.getSpectrumQuery().stream()
              .flatMapToInt(spectrumQuery -> spectrumQuery.getSearchResult().stream()
                  .mapToInt(value -> value.getSearchHit().size()))
          ).sum();
      System.out.printf("Done, found %d PSMs\n", sum);
    }
  }
1123553
def compare2(args):
    
    p = OptionParser(compare2.__doc__)
    p.add_option('--maxinsert', default=300, type="int",
                 help="Maximum number of repeats")
    add_simulate_options(p)
    opts, args, iopts = p.set_image_options(args, figsize="10x5")

    if len(args) != 0:
        sys.exit(not p.print_help())

    depth = opts.depth
    readlen = opts.readlen
    distance = opts.distance
    max_insert = opts.maxinsert
    fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1,
                                   figsize=(iopts.w, iopts.h))
    plt.tight_layout(pad=2)

    # ax1: lobSTR vs TREDPARSE with haploid model
    lobstr_results = parse_results("lobstr_results_homo.txt")
    tredparse_results = parse_results("tredparse_results_homo.txt")
    title = SIMULATED_HAPLOID + \
            r" ($D=%s\times, L=%dbp, V=%dbp$)" % (depth, readlen, distance)
    plot_compare(ax1, title, tredparse_results, lobstr_results,
                 max_insert=max_insert)

    # ax2: lobSTR vs TREDPARSE with diploid model
    lobstr_results = parse_results("lobstr_results_het.txt", exclude=20)
    tredparse_results = parse_results("tredparse_results_het.txt", exclude=20)
    title = SIMULATED_DIPLOID + \
            r" ($D=%s\times, L=%dbp, V=%dbp$)" % (depth, readlen, distance)
    plot_compare(ax2, title, tredparse_results, lobstr_results,
                 max_insert=max_insert)

    for ax in (ax1, ax2):
        ax.set_xlim(0, max_insert)
        ax.set_ylim(0, max_insert)

    root = fig.add_axes([0, 0, 1, 1])
    pad = .03
    panel_labels(root, ((pad / 2, 1 - pad, "A"), (1 / 2., 1 - pad, "B")))
    normalize_axes(root)

    image_name = "tredparse." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)