Run quality control in QC environment on a single sample.

    Enables peddy integration with CWL runs.
492311
def run_qc(_, data, out_dir):
    
    if cwlutils.is_cwl_run(data):
        qc_data = run_peddy([data], out_dir)
        if tz.get_in(["summary", "qc", "peddy"], qc_data):
            return tz.get_in(["summary", "qc", "peddy"], qc_data)
610464
def fixpoint_runner(self):
        """"""
        q = self.cfg.nodes

        while q != []:
            x_i = constraint_table[q[0]]  # x_i = q[0].old_constraint
            self.analysis.fixpointmethod(q[0])  # y = F_i(x_1, ..., x_n);
            y = constraint_table[q[0]]  # y = q[0].new_constraint

            if y != x_i:
                for node in self.analysis.dep(q[0]):  # for (v in dep(v_i))
                    q.append(node)  # q.append(v):
                constraint_table[q[0]] = y  # q[0].old_constraint = q[0].new_constraint # x_i = y
            q = q[1:]
930969
def check(context: click.Context, family: str):
    """"""
    analysis_obj = context.obj['store'].analyses(family=family).first()
    if analysis_obj is None:
        LOG.error('no analysis found')
        context.abort()

    config_path = Path(analysis_obj.config_path)
    if not config_path.exists():
        LOG.error(f"analysis config not found: {config_path}")
        context.abort()
    config_raw = ruamel.yaml.safe_load(config_path.open())
    config_data = files.parse_config(config_raw)

    sampleinfo_raw = ruamel.yaml.safe_load(Path(config_data['sampleinfo_path']).open())
    sampleinfo_data = files.parse_sampleinfo(sampleinfo_raw)

    qcmetrics_path = Path(sampleinfo_data['qcmetrics_path'])
    if not qcmetrics_path.exists():
        LOG.error(f"qc metrics not found: {str(qcmetrics_path)}")
        context.abort()
    qcmetrics_raw = ruamel.yaml.safe_load(qcmetrics_path.open())
    qcmetrics_data = files.parse_qcmetrics(qcmetrics_raw)

    samples = {
        'sample': [],
        'type': [],
        'ped': [],
        'chanjo': [],
        'peddy': [],
        'plink': [],
        'duplicates': [],
    }
    for sample_data in config_data['samples']:
        LOG.debug(f"{sample_data['id']}: parse analysis config")
        samples['sample'].append(sample_data['id'])
        samples['type'].append(sample_data['type'])

    for sample_data in sampleinfo_data['samples']:
        LOG.debug(f"{sample_data['id']}: parse sample info")
        samples['ped'].append(sample_data['sex'])

        with Path(sample_data['chanjo_sexcheck']).open() as chanjo_handle:
            sexcheck_data = files.parse_chanjo_sexcheck(chanjo_handle)

        predicted_sex = sexcheck_data['predicted_sex']
        xy_ratio = sexcheck_data['y_coverage'] / sexcheck_data['x_coverage']
        samples['chanjo'].append(f"{predicted_sex} ({xy_ratio:.3f})")

    for sample_data in qcmetrics_data['samples']:
        LOG.debug(f"{sample_data['id']}: parse qc metrics")
        samples['plink'].append(sample_data['plink_sex'])
        duplicates_percent = sample_data['duplicates'] * 100
        samples['duplicates'].append(f"{duplicates_percent:.3f}%")

    peddy_path = Path(sampleinfo_data['peddy']['sex_check'])
    if peddy_path.exists():
        with peddy_path.open() as sexcheck_handle:
            peddy_data = files.parse_peddy_sexcheck(sexcheck_handle)

        for sample_id in samples['sample']:
            LOG.debug(f"{sample_id}: parse peddy")
            predicted_sex = peddy_data[sample_id]['predicted_sex']
            het_ratio = peddy_data[sample_id]['het_ratio']
            samples['peddy'].append(f"{predicted_sex} ({het_ratio})")
    else:
        LOG.warning(f"missing peddy output: {peddy_path}")

    print(tabulate(samples, headers='keys', tablefmt='psql'))
798771
def run_analysis(self, argv):
        """"""
        args = self._parser.parse_args(argv)

        # Read the input maps
        ccube_dirty = HpxMap.create_from_fits(args.ccube_dirty, hdu='SKYMAP')
        bexpcube_dirty = HpxMap.create_from_fits(args.bexpcube_dirty, hdu='HPXEXPOSURES')
        ccube_clean = HpxMap.create_from_fits(args.ccube_clean, hdu='SKYMAP')
        bexpcube_clean = HpxMap.create_from_fits(args.bexpcube_clean, hdu='HPXEXPOSURES')

        # Decide what HEALPix order to work at
        if args.hpx_order:
            hpx_order = args.hpx_order
        else:
            hpx_order = ccube_dirty.hpx.order

        # Cast all the input maps to match ccube_clean
        cube_dict = ResidualCR._match_cubes(ccube_clean, ccube_dirty,
                                            bexpcube_clean, bexpcube_dirty, hpx_order)

        # Intenstiy maps
        intensity_clean = ResidualCR._compute_intensity(cube_dict['ccube_clean'],
                                                        cube_dict['bexpcube_clean'])
        intensity_dirty = ResidualCR._compute_intensity(cube_dict['ccube_dirty'],
                                                        cube_dict['bexpcube_dirty'])
        # Mean & ratio of intensity maps
        intensity_mean = ResidualCR._compute_mean(intensity_dirty,
                                                  intensity_clean)
        intensity_ratio = ResidualCR._compute_ratio(intensity_dirty,
                                                    intensity_clean)
        # Selecting the bright pixels for Aeff correction and to mask when filling output map
        bright_pixel_select = ResidualCR._make_bright_pixel_mask(intensity_mean,
                                                                 args.select_factor)
        bright_pixel_mask = ResidualCR._make_bright_pixel_mask(intensity_mean,
                                                               args.mask_factor)
        # Compute thte Aeff corrections using the brightest pixels
        aeff_corrections = ResidualCR._get_aeff_corrections(intensity_ratio,
                                                            bright_pixel_select)
        # Apply the Aeff corrections and get the intensity residual
        corrected_dirty = ResidualCR._apply_aeff_corrections(intensity_dirty,
                                                             aeff_corrections)
        corrected_ratio = ResidualCR._compute_ratio(corrected_dirty,
                                                    intensity_clean)
        intensity_resid = ResidualCR._compute_diff(corrected_dirty,
                                                   intensity_clean)
        # Replace the masked pixels with the map mean to avoid features associates with sources
        filled_resid = ResidualCR._fill_masked_intensity_resid(intensity_resid,
                                                               bright_pixel_mask)
        # Smooth the map
        smooth_resid = ResidualCR._smooth_hpx_map(filled_resid,
                                                  args.sigma)
        # Convert to a differential map
        out_model = ResidualCR._intergral_to_differential(smooth_resid)

        # Make the ENERGIES HDU
        out_energies = ccube_dirty.hpx.make_energies_hdu()

        # Write the maps
        cubes = dict(SKYMAP=out_model)
        fits_utils.write_maps(None, cubes,
                              args.outfile, energy_hdu=out_energies)

        if args.full_output:
            # Some diagnostics
            check = ResidualCR._differential_to_integral(out_model)
            check_resid = ResidualCR._compute_diff(smooth_resid, check)
            counts_resid =\
                ResidualCR._compute_counts_from_intensity(intensity_resid,
                                                          cube_dict['bexpcube_dirty'])
            pred_counts\
                = ResidualCR._compute_counts_from_model(out_model,
                                                        cube_dict['bexpcube_dirty'])
            pred_resid = ResidualCR._compute_diff(pred_counts, counts_resid)

            out_ebounds = ccube_dirty.hpx.make_energy_bounds_hdu()
            cubes = dict(INTENSITY_CLEAN=intensity_clean,
                         INTENSITY_DIRTY=intensity_dirty,
                         INTENSITY_RATIO=intensity_ratio,
                         CORRECTED_DIRTY=corrected_dirty,
                         CORRECTED_RATIO=corrected_ratio,
                         INTENSITY_RESID=intensity_resid,
                         PIXEL_SELECT=bright_pixel_select,
                         PIXEL_MASK=bright_pixel_mask,
                         FILLED_RESID=filled_resid,
                         SMOOTH_RESID=smooth_resid,
                         CHECK=check,
                         CHECK_RESID=check_resid,
                         COUNTS_RESID=counts_resid,
                         PRED_COUNTS=pred_counts,
                         PRED_RESID=pred_resid)

            fits_utils.write_maps(None, cubes,
                                  args.outfile.replace('.fits', '_full.fits'),
                                  energy_hdu=out_ebounds)
424881
def run(self):
        """"""
        if self.target_eff_nsamples and self.checkpoint_interval is None:
            raise ValueError("A checkpoint interval must be set if "
                             "targetting an effective number of samples")
        # get the starting number of samples:
        # "nsamples" keeps track of the number of samples we've obtained (if
        # target_eff_nsamples is not None, this is the effective number of
        # samples; otherwise, this is the total number of samples).
        # _lastclear is the number of iterations that the file already
        # contains (either due to sampler burn-in, or a previous checkpoint)
        if self.new_checkpoint:
            self._lastclear = 0
        else:
            with self.io(self.checkpoint_file, "r") as fp:
                self._lastclear = fp.niterations
        if self.target_eff_nsamples is not None:
            target_nsamples = self.target_eff_nsamples
            with self.io(self.checkpoint_file, "r") as fp:
                nsamples = fp.effective_nsamples
        elif self.target_niterations is not None:
            # the number of samples is the number of iterations times the
            # number of walkers
            target_nsamples = self.nwalkers * self.target_niterations
            nsamples = self._lastclear * self.nwalkers
        else:
            raise ValueError("must set either target_eff_nsamples or "
                             "target_niterations; see set_target")
        self._itercounter = 0
        # figure out the interval to use
        iterinterval = self.checkpoint_interval
        if iterinterval is None:
            iterinterval = self.target_niterations
        # run sampler until we have the desired number of samples
        while nsamples < target_nsamples:
            # adjust the interval if we would go past the number of iterations
            if self.target_niterations is not None and (
                    self.niterations + iterinterval > self.target_niterations):
                iterinterval = self.target_niterations - self.niterations
            # run sampler and set initial values to None so that sampler
            # picks up from where it left off next call
            logging.info("Running sampler for {} to {} iterations".format(
                self.niterations, self.niterations + iterinterval))
            # run the underlying sampler for the desired interval
            self.run_mcmc(iterinterval)
            # update the itercounter
            self._itercounter = self._itercounter + iterinterval
            # dump the current results
            self.checkpoint()
            # update nsamples for next loop
            if self.target_eff_nsamples is not None:
                nsamples = self.effective_nsamples
                logging.info("Have {} effective samples post burn in".format(
                    nsamples))
            else:
                nsamples += iterinterval * self.nwalkers
462415
def _verify_app_source_dir_impl(src_dir, temp_dir, mode, enforce=True):
    
    _lint(os.path.join(src_dir, "dxapp.json"), mode)

    # Check that the entry point file parses as the type it is going to
    # be interpreted as. The extension is irrelevant.
    manifest = json.load(open(os.path.join(src_dir, "dxapp.json")))
    if "runSpec" in manifest:
        if "interpreter" not in manifest['runSpec']:
            raise dxpy.app_builder.AppBuilderException('runSpec.interpreter field was not present')

        if "release" not in manifest['runSpec'] or "distribution" not in manifest['runSpec']:
            warn_message = 'runSpec.distribution or runSpec.release was not present. These fields '
            warn_message += 'will be required in a future version of the API. Recommended value '
            warn_message += 'for distribution is \"Ubuntu\" and release - \"14.04\".'
            logger.warn(warn_message)

        if manifest['runSpec']['interpreter'] in ["python2.7", "bash"]:
            if "file" in manifest['runSpec']:
                entry_point_file = os.path.abspath(os.path.join(src_dir, manifest['runSpec']['file']))
                try:
                    _check_file_syntax(entry_point_file, temp_dir, override_lang=manifest['runSpec']['interpreter'], enforce=enforce)
                except IOError as e:
                    raise dxpy.app_builder.AppBuilderException(
                        'Could not open runSpec.file=%r. The problem was: %s' % (entry_point_file, e))
                except DXSyntaxError:
                    raise dxpy.app_builder.AppBuilderException('Entry point file %s has syntax errors, see above for details. Rerun with --no-check-syntax to proceed anyway.' % (entry_point_file,))
            elif "code" in manifest['runSpec']:
                try:
                    _check_syntax(manifest['runSpec']['code'], manifest['runSpec']['interpreter'], temp_dir, enforce=enforce)
                except DXSyntaxError:
                    raise dxpy.app_builder.AppBuilderException('Code in runSpec.code has syntax errors, see above for details. Rerun with --no-check-syntax to proceed anyway.')

            cluster_bootstrap_scripts = _get_all_cluster_bootstrap_script_names(manifest)
            for filename in cluster_bootstrap_scripts:
                abs_filename = os.path.abspath(os.path.join(src_dir, filename))
                try:
                    _check_file_syntax(abs_filename,
                                       temp_dir,
                                       override_lang=manifest['runSpec']['interpreter'],
                                       enforce=enforce)
                except IOError as e:
                    raise dxpy.app_builder.AppBuilderException(
                        'Could not open cluster bootstrap script %r. The problem was: %s' % (abs_filename, e))
                except DXSyntaxError:
                    raise dxpy.app_builder.AppBuilderException('Code in cluster bootstrapScript %r has syntax errors, see above for details. Rerun with --no-check-syntax to proceed anyway.' % filename)

        if 'execDepends' in manifest['runSpec']:
            if not isinstance(manifest['runSpec']['execDepends'], list):
                raise dxpy.app_builder.AppBuilderException('Expected runSpec.execDepends to be an array. Rerun with --no-check-syntax to proceed anyway.')
            if not all(isinstance(dep, dict) for dep in manifest['runSpec']['execDepends']):
                raise dxpy.app_builder.AppBuilderException('Expected runSpec.execDepends to be an array of hashes. Rerun with --no-check-syntax to proceed anyway.')
            if any(dep.get('package_manager', 'apt') != 'apt' for dep in manifest['runSpec']['execDepends']):
                if not isinstance(manifest.get('access'), dict) or 'network' not in manifest['access']:
                    msg = '\n'.join(['runSpec.execDepends specifies non-APT dependencies, but no network access spec is given.',
                    'Add {"access": {"network": ["*"]}} to allow dependencies to install.',
                    'See https://wiki.dnanexus.com/Developer-Tutorials/Request-Additional-App-Resources#Network-Access.',
                    'Rerun with --no-check-syntax to proceed anyway.'])
                    raise dxpy.app_builder.AppBuilderException(msg)

    if 'authorizedUsers' in manifest:
        if not isinstance(manifest['authorizedUsers'], list) or isinstance(manifest['authorizedUsers'], basestring):
            raise dxpy.app_builder.AppBuilderException('Expected authorizedUsers to be a list of strings')
        for thing in manifest['authorizedUsers']:
            if thing != 'PUBLIC' and (not isinstance(thing, basestring) or not re.match("^(org-|user-)", thing)):
                raise dxpy.app_builder.AppBuilderException('authorizedUsers field contains an entry which is not either the string "PUBLIC" or a user or org ID')

    if "pricingPolicy" in manifest:
        error_message = "\"pricingPolicy\" at the top level is not accepted. It must be specified "
        error_message += "under the \"regionalOptions\" field in all enabled regions of the app"
        raise dxpy.app_builder.AppBuilderException(error_message)

    # Check all other files that are going to be in the resources tree.
    # For these we detect the language based on the filename extension.
    # Obviously this check can have false positives, since the app can
    # execute (or not execute!) all these files in whatever way it
    # wishes, e.g. it could use Python != 2.7 or some non-bash shell.
    # Consequently errors here are non-fatal.
    files_with_problems = []
    for dirpath, dirnames, filenames in os.walk(os.path.abspath(os.path.join(src_dir, "resources"))):
        for filename in filenames:
            # On Mac OS, the resource fork for "FILE.EXT" gets tarred up
            # as a file named "._FILE.EXT". To a naive check this
            # appears to be a file of the same extension. Therefore, we
            # exclude these from syntax checking since they are likely
            # to not parse as whatever language they appear to be.
            if not filename.startswith("._"):
                try:
                    _check_file_syntax(os.path.join(dirpath, filename), temp_dir, enforce=True)
                except IOError as e:
                    raise dxpy.app_builder.AppBuilderException(
                        'Could not open file in resources directory %r. The problem was: %s' %
                        (os.path.join(dirpath, filename), e)
                    )
                except DXSyntaxError:
                    # Suppresses errors from _check_file_syntax so we
                    # only print a nice error message
                    files_with_problems.append(os.path.join(dirpath, filename))

    if files_with_problems:
        # Make a message of the form:
        #    "/path/to/my/app.py"
        # OR "/path/to/my/app.py and 3 other files"
        files_str = files_with_problems[0] if len(files_with_problems) == 1 else (files_with_problems[0] + " and " + str(len(files_with_problems) - 1) + " other file" + ("s" if len(files_with_problems) > 2 else ""))
        logging.warn('%s contained syntax errors, see above for details' % (files_str,))
553579
def sample(config_file, profile, solver_def, biases, couplings, random_problem,
           num_reads, verbose):
    """"""

    # TODO: de-dup wrt ping

    def echo(s, maxlen=100):
        click.echo(s if verbose else strtrunc(s, maxlen))

    try:
        client = Client.from_config(
            config_file=config_file, profile=profile, solver=solver_def)
    except Exception as e:
        click.echo("Invalid configuration: {}".format(e))
        return 1
    if config_file:
        echo("Using configuration file: {}".format(config_file))
    if profile:
        echo("Using profile: {}".format(profile))
    echo("Using endpoint: {}".format(client.endpoint))

    try:
        solver = client.get_solver()
    except SolverAuthenticationError:
        click.echo("Authentication error. Check credentials in your configuration file.")
        return 1
    except (InvalidAPIResponseError, UnsupportedSolverError):
        click.echo("Invalid or unexpected API response.")
        return 2
    except SolverNotFoundError:
        click.echo("Solver with the specified features does not exist.")
        return 3

    echo("Using solver: {}".format(solver.id))

    if random_problem:
        linear, quadratic = generate_random_ising_problem(solver)
    else:
        try:
            linear = ast.literal_eval(biases) if biases else []
        except Exception as e:
            click.echo("Invalid biases: {}".format(e))
        try:
            quadratic = ast.literal_eval(couplings) if couplings else {}
        except Exception as e:
            click.echo("Invalid couplings: {}".format(e))

    echo("Using qubit biases: {!r}".format(linear))
    echo("Using qubit couplings: {!r}".format(quadratic))
    echo("Number of samples: {}".format(num_reads))

    try:
        result = solver.sample_ising(linear, quadratic, num_reads=num_reads).result()
    except Exception as e:
        click.echo(e)
        return 4

    if verbose:
        click.echo("Result: {!r}".format(result))

    echo("Samples: {!r}".format(result['samples']))
    echo("Occurrences: {!r}".format(result['occurrences']))
    echo("Energies: {!r}".format(result['energies']))
207893
def build_spec(user, repo, sha=None, prov=None, extraMetadata=[]):
    """"""
    loader = grlc.utils.getLoader(user, repo, sha=sha, prov=prov)

    files = loader.fetchFiles()
    raw_repo_uri = loader.getRawRepoUri()

    # Fetch all .rq files
    items = []

    allowed_ext = ["rq", "sparql", "json", "tpf"]
    for c in files:
        glogger.debug('>>>>>>>>>>>>>>>>>>>>>>>>>c_name: {}'.format(c['name']))
        extension = c['name'].split('.')[-1]
        if extension in allowed_ext:
            call_name = c['name'].split('.')[0]

            # Retrieve extra metadata from the query decorators
            query_text = loader.getTextFor(c)

            item = None
            if extension == "json":
                query_text = json.loads(query_text)

            if extension in ["rq", "sparql", "json"]:
                glogger.debug("===================================================================")
                glogger.debug("Processing SPARQL query: {}".format(c['name']))
                glogger.debug("===================================================================")
                item = process_sparql_query_text(query_text, loader, call_name, extraMetadata)
            elif "tpf" == extension:
                glogger.debug("===================================================================")
                glogger.debug("Processing TPF query: {}".format(c['name']))
                glogger.debug("===================================================================")
                item = process_tpf_query_text(query_text, raw_repo_uri, call_name, extraMetadata)
            else:
                glogger.info("Ignoring unsupported source call name: {}".format(c['name']))

            if item:
                items.append(item)
    return items
1355096
def analyze_multiple(
            benchmark_runs,
            analysis_sets = [],
            # Singleton arguments
            analysis_directory = None,
            remove_existing_analysis_directory = True,
            quick_plots = False,
            use_multiprocessing = True,
            verbose = True,
            compile_pdf = True,
            limit_to_complete_presence = True,
            all_by_all_comparisons = False,
    ):
        ''''''
        if remove_existing_analysis_directory and os.path.isdir(analysis_directory):
            shutil.rmtree(analysis_directory)

        unique_ajps = BenchmarkRun.get_unique_ajps( benchmark_runs )

        if limit_to_complete_presence:
            BenchmarkRun.get_common_subset( benchmark_runs, verbose = not use_multiprocessing )

        unique_ajps = BenchmarkRun.get_unique_ajps( benchmark_runs )

        ###  Process each benchmark run object individually
        if use_multiprocessing:
            pool = mp.Pool()
        singleton_chapters = []
        calculated_brs = []
        def save_latex_report(t):
            br, unique_name, latex_report = t
            latex_report.set_title_page( title = unique_name )
            singleton_chapters.append( latex_report )
            calculated_brs.append( br )
        for br in benchmark_runs:
            for analysis_set in analysis_sets:
                unique_name = br.get_definitive_name(unique_ajps, join_character = '\n')
                filepath_unique_name = br.get_definitive_name(unique_ajps, join_character = '-')
                subdir = os.path.join(analysis_directory, os.path.join('analysis_sets', os.path.join(analysis_set, filepath_unique_name) ) )
                if use_multiprocessing:
                    pool.apply_async( _full_analysis_mp_alias, ( br, analysis_set, subdir, unique_name, False, quick_plots ), callback = save_latex_report )
                else:
                    print 'Individual report saving in:', subdir
                    save_latex_report( _full_analysis_mp_alias( br, analysis_set, subdir, unique_name, True, quick_plots ) )
        if use_multiprocessing:
            pool.close()
            pool.join()
        benchmark_runs = calculated_brs

        ### Pointwise all-by-all comparison
        comparison_chapters = []
        if all_by_all_comparisons:
            if use_multiprocessing:
                pool = mp.Pool()
            def save_latex_report(t):
                latex_report = t
                comparison_chapters.append( latex_report )
            comparisons_subdir = os.path.join(analysis_directory, 'comparison_analysis_sets')
            for analysis_set in analysis_sets:
                analysis_set_subdir = os.path.join(comparisons_subdir, analysis_set)
                for i, br_i in enumerate(benchmark_runs):
                    for j, br_j in enumerate(benchmark_runs):
                        if i > j:
                            if use_multiprocessing:
                                br_i_copy = copy.deepcopy( br_i )
                                br_j_copy = copy.deepcopy( br_j )
                                pool.apply_async( _compare_mp_alias, (br_i_copy, br_j_copy, analysis_set, analysis_set_subdir, unique_ajps, False), callback = save_latex_report )
                            else:
                                save_latex_report( _compare_mp_alias(br_i, br_j, analysis_set, analysis_set_subdir, unique_ajps, True) )
            if use_multiprocessing:
                pool.close()
                pool.join()

        intro_report = lr.LatexReport()
        intro_report.set_title_page('All data comparison')
        # All data series comparison
        # Get joined stats comparison dataframe
        stats_df = BenchmarkRun.get_stats_comparison_dataframe(
            benchmark_runs, unique_ajps,
            output_csv = os.path.join(analysis_directory, 'analysis_metrics.csv'),
        )
        intro_report.add_section_page( title = 'Case comparison tables' )
        intro_report.content.extend( BenchmarkRun.make_case_description_tables( stats_df ) )

        intro_report.add_section_page('All data plots')
        subplot_directory = os.path.join(analysis_directory, 'subplots')
        if not os.path.isdir( subplot_directory ):
            os.makedirs(subplot_directory)
        runtime_df = benchmark_runs[0]._get_dataframe_columns( ['RunTime'] )
        runtime_df.columns = [ benchmark_runs[0].get_definitive_name(unique_ajps, join_character = '\n', prepend_label = False) ]
        for br in benchmark_runs[1:]:
            inner_runtime_df = br._get_dataframe_columns( ['RunTime'] )
            inner_runtime_df.columns = [ br.get_definitive_name(unique_ajps, join_character = '\n', prepend_label = False) ]
            runtime_df = runtime_df.merge(
                inner_runtime_df,
                left_index = True,
                right_index = True,
            )

        intro_report.add_plot(
            general_matplotlib.plot_box(
                runtime_df,
                output_directory = subplot_directory,
                plot_title = 'Prediction Run Times',
                output_name = 'runtimes',
                fig_height = 6.7,
                fig_width = 10,
                ylabel = 'Run time (minutes)',
                xlabel = 'Prediction Set',
                verbose = verbose,
                xtick_fontsize = 4,
                log_y = True,
                label_n = False,
                rotation_angle = 45,
            ),
            plot_title = 'Run times'
        )

        # Report concatenation
        main_latex_report = lr.LatexReport()
        main_latex_report.set_title_page('$\Delta\Delta G$ Report')
        main_latex_report.add_chapter(intro_report)
        for chapter in comparison_chapters:
            main_latex_report.add_chapter(chapter)
        for chapter in singleton_chapters:
            main_latex_report.add_chapter(chapter)
        main_latex_report.generate_pdf_report(
            os.path.join( analysis_directory, 'report.pdf' ),
            verbose = verbose,
            compile_pdf = compile_pdf,
        )
        print os.path.join( analysis_directory, 'report.pdf' )
1236314
def _check(sample, data):
    """"""
    if dd.get_chip_method(sample).lower() == "atac":
        return [sample]
    if dd.get_phenotype(sample) == "input":
        return None
    for origin in data:
        if dd.get_batch(sample) in (dd.get_batches(origin[0]) or []) and dd.get_phenotype(origin[0]) == "input":
            sample["work_bam_input"] = origin[0].get("work_bam")
            return [sample]
    return [sample]
159065
def _correct_qualimap_genome_results(samples):
    
    for s in samples:
        if verify_file(s.qualimap_genome_results_fpath):
            correction_is_needed = False
            with open(s.qualimap_genome_results_fpath, 'r') as f:
                content = f.readlines()
                metrics_started = False
                for line in content:
                    if ">> Reference" in line:
                        metrics_started = True
                    if metrics_started:
                        if line.find(',') != -1:
                            correction_is_needed = True
                            break
            if correction_is_needed:
                with open(s.qualimap_genome_results_fpath, 'w') as f:
                    metrics_started = False
                    for line in content:
                        if ">> Reference" in line:
                            metrics_started = True
                        if metrics_started:
                            if line.find(',') != -1:
                                line = line.replace(',', '')
                        f.write(line)
833806
def _single_multichannel_mp_run(X, Phi, bound, selection_rule, stop_crit,
                                max_iter, verbose=False, pad=0,
                                random_state=None, memory=Memory(None)):
    
    rng = check_random_state(random_state)

    # padding as v stak
    pad = int(pad)
    n_channels = X.shape[0]
    X = np.hstack((np.zeros((n_channels, pad)), X,
                   np.zeros((n_channels, pad))))
    n_samples = X.shape[1]
    n_projs = Phi.doth(X).shape[1]
    err_mse = {}

    # Initialisation
    residual = np.hstack((X.copy(), np.zeros((n_channels,
                                              max(Phi.sizes) / 2))))

    s_rep = np.zeros((n_channels, n_projs))
    X_est = np.zeros((n_channels, n_samples))
    # Main algorithm
    coeffs = np.zeros((n_channels, n_projs))

    it_number = 0
    current_lambda = 1
    for c_idx in range(n_channels):
        err_mse[c_idx] = []
        err_mse[c_idx].append(linalg.norm(residual[c_idx, :]))

    # Decomposition loop: stopping criteria is either SNR or iteration number
    while (current_lambda > bound) & (it_number < max_iter):

        # pick a shift at random : in each size
        rndshifts = {}
        for c_idx in range(n_channels):
            rndshifts[c_idx] = []
        for s_idx, L in enumerate(Phi.sizes):
            shift = rng.randint(low=0, high=L / 4)
            for c_idx in range(n_channels):
                coeffs[c_idx, s_idx * n_samples:(s_idx + 1) * n_samples] = \
                    mdct(residual[c_idx, shift:shift + n_samples], L).ravel()
                rndshifts[c_idx].append(shift)

        # Multichannel mode : we combine projections
        combined = selection_rule(coeffs ** 2)

        # Select a new element
        idx = np.argmax(np.abs(combined))
        # find scale and frequency bin of selected atom
        s_idx = idx // n_samples
        L = Phi.sizes[s_idx]
        F = n_samples // (L // 2)
        frame = (idx - (s_idx * n_samples)) % F
        freq_bin = ((idx - (s_idx * n_samples))) // F

        mdct_wf = memory.cache(mdct_waveform)

        # Update coefficients and residual
        current_lambda_array = np.zeros(n_channels)
        for c_idx in range(n_channels):
            s_rep[c_idx, idx] += coeffs[c_idx, idx]

            # Only one method now : local update via a cached waveform
            pos = (frame * L / 2) - L / 4 + rndshifts[c_idx][s_idx]
            residual[c_idx, pos:pos + L] -= coeffs[c_idx, idx] * \
                mdct_wf(L, freq_bin)

            # also add it to the reconstruction
            X_est[c_idx, pos:pos + L] += coeffs[c_idx, idx] * \
                mdct_wf(L, freq_bin)

            # error computation (err_mse)
            err_mse[c_idx].append(linalg.norm(residual[c_idx, :]))

            current_lambda_array[c_idx] = np.sqrt(
                1. - err_mse[c_idx][-1] / err_mse[c_idx][-2])

        current_lambda = stop_crit(current_lambda_array)

        if verbose:
            print("Iteration %d : Current lambda of %1.4f" % (
                  it_number, current_lambda))
        it_number += 1

    return X_est[:, pad: -pad], err_mse
48350
def composite_qc(df_orig, size=(16, 12)):
    
    df = df_orig.rename(columns={"hli_calc_age_sample_taken": "Age",
                       "hli_calc_gender": "Gender",
                       "eth7_max": "Ethnicity",
                       "MeanCoverage": "Mean coverage",
                       "Chemistry": "Sequencing chemistry",
                       "Release Client": "Cohort",

                      })

    fig = plt.figure(1, size)
    ax1 = plt.subplot2grid((2, 7), (0, 0), rowspan=1, colspan=2)
    ax2 = plt.subplot2grid((2, 7), (0, 2), rowspan=1, colspan=2)
    ax3 = plt.subplot2grid((2, 7), (0, 4), rowspan=1, colspan=3)
    ax4 = plt.subplot2grid((2, 7), (1, 0), rowspan=1, colspan=2)
    ax5 = plt.subplot2grid((2, 7), (1, 2), rowspan=1, colspan=2)
    ax6 = plt.subplot2grid((2, 7), (1, 4), rowspan=1, colspan=3)

    sns.distplot(df["Age"].dropna(), kde=False, ax=ax1)
    sns.countplot(x="Gender", data=df, ax=ax2)
    sns.countplot(x="Ethnicity", data=df, ax=ax3,
                    order = df['Ethnicity'].value_counts().index)
    sns.distplot(df["Mean coverage"].dropna(), kde=False, ax=ax4)
    ax4.set_xlim(0, 100)
    sns.countplot(x="Sequencing chemistry", data=df, ax=ax5)
    sns.countplot(x="Cohort", data=df, ax=ax6,
                    order = df['Cohort'].value_counts().index)
    # Anonymize the cohorts
    cohorts = ax6.get_xticklabels()
    newCohorts = []
    for i, c in enumerate(cohorts):
        if c.get_text() == "Spector":
            c = "TwinsUK"
        elif c.get_text() != "Health Nucleus":
            c = "C{}".format(i + 1)
        newCohorts.append(c)
    ax6.set_xticklabels(newCohorts)

    for ax in (ax6,):
        ax.set_xticklabels(ax.get_xticklabels(), ha="right", rotation=30)

    for ax in (ax1, ax2, ax3, ax4, ax5, ax6):
        ax.set_title(ax.get_xlabel())
        ax.set_xlabel("")

    plt.tight_layout()

    root = fig.add_axes((0, 0, 1, 1))
    labels = ((.02, .96, "A"),
              (.3, .96, "B"),
              (.6, .96, "C"),
              (.02, .52, "D"),
              (.3, .52, "E"),
              (.6, .52, "F"))
    panel_labels(root, labels)
    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()
245123
def _pair_samples_with_pipelines(run_info_yaml, config):
    
    samples = config_utils.load_config(run_info_yaml)
    if isinstance(samples, dict):
        resources = samples.pop("resources")
        samples = samples["details"]
    else:
        resources = {}
    ready_samples = []
    for sample in samples:
        if "files" in sample:
            del sample["files"]
        # add any resources to this item to recalculate global configuration
        usample = copy.deepcopy(sample)
        usample.pop("algorithm", None)
        if "resources" not in usample:
            usample["resources"] = {}
        for prog, pkvs in resources.items():
            if prog not in usample["resources"]:
                usample["resources"][prog] = {}
            if pkvs is not None:
                for key, val in pkvs.items():
                    usample["resources"][prog][key] = val
        config = config_utils.update_w_custom(config, usample)
        sample["resources"] = {}
        ready_samples.append(sample)
    paired = [(x, _get_pipeline(x)) for x in ready_samples]
    d = defaultdict(list)
    for x in paired:
        d[x[1]].append([x[0]])
    return d, config
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))
851238
def setup_interval_coinc_inj(workflow, hdfbank, full_data_trig_files, inj_trig_files,
              stat_files, background_file, veto_file, veto_name, out_dir, tags=None):
    
    if tags is None:
        tags = []
    make_analysis_dir(out_dir)
    logging.info('Setting up coincidence for injection')

    if len(hdfbank) > 1:
        raise ValueError('This coincidence method only supports a '
                         'pregenerated template bank')
    hdfbank = hdfbank[0]

    if len(workflow.ifos) > 2:
        raise ValueError('This coincidence method only supports two-ifo searches')

    # Wall time knob and memory knob
    factor = int(workflow.cp.get_opt_tags('workflow-coincidence', 'parallelization-factor', tags))

    ffiles = {}
    ifiles = {}
    for ifo, ffi in zip(*full_data_trig_files.categorize_by_attr('ifo')):
        ffiles[ifo] = ffi[0]
    ifos, files = inj_trig_files.categorize_by_attr('ifo')  # ifos list is used later
    for ifo, ifi in zip(ifos, files):
        ifiles[ifo] = ifi[0]
    ifo0, ifo1 = ifos[0], ifos[1]
    combo = [(FileList([ifiles[ifo0], ifiles[ifo1]]), "injinj"),
             (FileList([ifiles[ifo0], ffiles[ifo1]]), "injfull"),
             (FileList([ifiles[ifo1], ffiles[ifo0]]), "fullinj"),
            ]
    bg_files = {'injinj':[], 'injfull':[], 'fullinj':[]}

    for trig_files, ctag in combo:
        findcoinc_exe = PyCBCFindCoincExecutable(workflow.cp, 'coinc',
                                              ifos=workflow.ifos,
                                              tags=tags + [ctag], out_dir=out_dir)
        for i in range(factor):
            group_str = '%s/%s' % (i, factor)
            coinc_node = findcoinc_exe.create_node(trig_files, hdfbank,
                                                   stat_files,
                                                   veto_file, veto_name,
                                                   group_str,
                                                   tags=[str(i)])
            bg_files[ctag] += coinc_node.output_files
            workflow.add_node(coinc_node)

    return setup_statmap_inj(workflow, bg_files, background_file, hdfbank, out_dir, tags=tags)
458703
def variant(samples):
    
    checkpoints = _variant_checkpoints(samples)
    if checkpoints["align"]:
        align_wf = _alignment(checkpoints)
        alignin = [["files"], ["analysis"],
                   ["config", "algorithm", "align_split_size"],
                   ["reference", "fasta", "base"],
                   ["rgnames", "pl"], ["rgnames", "sample"], ["rgnames", "pu"],
                   ["rgnames", "lane"], ["rgnames", "rg"], ["rgnames", "lb"],
                   ["reference", "aligner", "indexes"],
                   ["config", "algorithm", "aligner"],
                   ["config", "algorithm", "trim_reads"],
                   ["config", "algorithm", "adapters"],
                   ["config", "algorithm", "bam_clean"],
                   ["config", "algorithm", "variant_regions"],
                   ["config", "algorithm", "mark_duplicates"]]
        if checkpoints["hla"]:
            alignin.append(["config", "algorithm", "hlacaller"])
        if checkpoints["umi"]:
            alignin.append(["config", "algorithm", "umi_type"])
        align = [s("alignment_to_rec", "multi-combined", alignin,
                   [cwlout("alignment_rec", "record")],
                   "bcbio-vc",
                   disk={"files": 1.5}, cores=1, no_files=True),
                 w("alignment", "multi-parallel", align_wf,
                   [["align_split"], ["process_alignment_rec"],
                    ["work_bam"], ["config", "algorithm", "quality_format"]])]
    else:
        align = [s("organize_noalign", "multi-parallel",
                   ["files"],
                   [cwlout(["align_bam"], ["File", "null"], [".bai"]),
                    cwlout(["work_bam_plus", "disc"], ["File", "null"]),
                    cwlout(["work_bam_plus", "sr"], ["File", "null"]),
                    cwlout(["hla", "fastq"], ["File", "null"])],
                   "bcbio-vc", cores=1)]
    align_out = [["rgnames", "sample"], ["align_bam"]]
    pp_align, pp_align_out = _postprocess_alignment(checkpoints)
    if checkpoints["umi"]:
        align_out += [["umi_bam"]]
    vc, vc_out = _variant_vc(checkpoints)
    sv, sv_out = _variant_sv(checkpoints)
    hla, hla_out = _variant_hla(checkpoints)
    qc, qc_out = _qc_workflow(checkpoints)
    steps = align + pp_align + hla + vc + sv + qc
    final_outputs = align_out + pp_align_out + vc_out + hla_out + sv_out + qc_out
    return steps, final_outputs