Prepares RAPID streamflow for GSSHA simulation
498567
def prepare_rapid_streamflow(self, path_to_rapid_qout, connection_list_file):
        
        ihg_filename = '{0}.ihg'.format(self.project_manager.name)
        with tmp_chdir(self.project_manager.project_directory):
            # write out IHG file
            time_index_range = []
            with RAPIDDataset(path_to_rapid_qout, out_tzinfo=self.tz) as qout_nc:

                time_index_range = qout_nc.get_time_index_range(date_search_start=self.simulation_start,
                                                                date_search_end=self.simulation_end)

                if len(time_index_range) > 0:
                    time_array = qout_nc.get_time_array(return_datetime=True,
                                                        time_index_array=time_index_range)

                    # GSSHA STARTS INGESTING STREAMFLOW AT SECOND TIME STEP
                    if self.simulation_start is not None:
                        if self.simulation_start == time_array[0]:
                            log.warning("First timestep of streamflow skipped "
                                     "in order for GSSHA to capture the streamflow.")
                            time_index_range = time_index_range[1:]
                            time_array = time_array[1:]

                if len(time_index_range) > 0:
                    start_datetime = time_array[0]

                    if self.simulation_start is None:
                       self._update_simulation_start(start_datetime)

                    if self.simulation_end is None:
                        self.simulation_end = time_array[-1]

                    qout_nc.write_flows_to_gssha_time_series_ihg(ihg_filename,
                                                                 connection_list_file,
                                                                 date_search_start=start_datetime,
                                                                 date_search_end=self.simulation_end,
                                                                 )
                else:
                    log.warning("No streamflow values found in time range ...")

            if len(time_index_range) > 0:
                # update cards
                self._update_simulation_start_cards()

                self._update_card("END_TIME", self.simulation_end.strftime("%Y %m %d %H %M"))
                self._update_card("CHAN_POINT_INPUT", ihg_filename, True)

                # update duration
                self.set_simulation_duration(self.simulation_end-self.simulation_start)

                # UPDATE GMT CARD
                self._update_gmt()
            else:
                # cleanup
                os.remove(ihg_filename)
                self.project_manager.deleteCard('CHAN_POINT_INPUT', self.db_session)
64201
def run(self, subdirectory=None):
        
        with tmp_chdir(self.gssha_directory):
            if self.hotstart_minimal_mode:
                # remove all optional output cards
                for gssha_optional_output_card in self.GSSHA_OPTIONAL_OUTPUT_CARDS:
                    self._delete_card(gssha_optional_output_card)
                # make sure running in SUPER_QUIET mode
                self._update_card('SUPER_QUIET', '')
                if subdirectory is None:
                    # give execute folder name
                    subdirectory = "minimal_hotstart_run_{0}to{1}" \
                                   .format(self.event_manager.simulation_start.strftime("%Y%m%d%H%M"),
                                           self.event_manager.simulation_end.strftime("%Y%m%d%H%M"))
            else:
                # give execute folder name
                subdirectory = "run_{0}to{1}".format(self.event_manager.simulation_start.strftime("%Y%m%d%H%M"),
                                                     self.event_manager.simulation_end.strftime("%Y%m%d%H%M"))

            # ensure unique folder naming conventions and add to exisitng event manager
            prj_evt_manager = self.project_manager.projectFileEventManager
            prj_event = prj_evt_manager.add_event(name=subdirectory,
                                                  subfolder=subdirectory,
                                                  session=self.db_session)
            eventyml_path = self.project_manager.getCard('#GSSHAPY_EVENT_YML') \
                                                .value.strip("'").strip('"')
            prj_evt_manager.write(session=self.db_session,
                                  directory=self.gssha_directory,
                                  name=os.path.basename(eventyml_path))
            # ensure event manager not propagated to child event
            self.project_manager.deleteCard('#GSSHAPY_EVENT_YML',
                                            db_session=self.db_session)
            self.db_session.delete(self.project_manager.projectFileEventManager)
            self.db_session.commit()

            # make working directory
            working_directory = os.path.join(self.gssha_directory, prj_event.subfolder)
            try:
                os.mkdir(working_directory)
            except OSError:
                pass

            # move simulation generated files to working directory
            # PRECIP_FILE, HMET_NETCDF, HMET_ASCII, CHAN_POINT_INPUT
            # TODO: Move HMET_ASCII files
            for sim_card in self.simulation_modified_input_cards:
                if sim_card != 'MAPPING_TABLE':
                    self._update_card_file_location(sim_card, working_directory)

            mapping_table_card = self.project_manager.getCard('MAPPING_TABLE')
            if mapping_table_card:
                # read in mapping table
                map_table_object = self.project_manager.readInputFile('MAPPING_TABLE',
                                                                      self.gssha_directory,
                                                                      self.db_session,
                                                                      readIndexMaps=False)

                # connect index maps to main gssha directory
                for indexMap in map_table_object.indexMaps:
                    indexMap.filename = os.path.join("..", os.path.basename(indexMap.filename))

                # write copy of mapping table to working directory
                map_table_filename = os.path.basename(mapping_table_card.value.strip("'").strip('"'))
                map_table_object.write(session=self.db_session,
                                       directory=working_directory,
                                       name=map_table_filename,
                                       writeIndexMaps=False)

            # connect to other output files in main gssha directory
            for gssha_card in self.project_manager.projectCards:
                if gssha_card.name not in self.GSSHA_REQUIRED_OUTPUT_PATH_CARDS + \
                                            self.GSSHA_OPTIONAL_OUTPUT_PATH_CARDS + \
                                            tuple(self.simulation_modified_input_cards):
                    if gssha_card.value:
                        updated_value = gssha_card.value.strip('"').strip("'")
                        if updated_value:
                            if gssha_card.name == "READ_CHAN_HOTSTART":
                                # there are two required files
                                # the .dht and .qht
                                if os.path.exists(updated_value + '.dht') \
                                        and os.path.exists(updated_value + '.qht'):
                                    updated_path = os.path.join("..", os.path.basename(updated_value))
                                    gssha_card.value = '"{0}"'.format(updated_path)
                            elif os.path.exists(updated_value):
                                updated_path = os.path.join("..", os.path.basename(updated_value))
                                gssha_card.value = '"{0}"'.format(updated_path)
                            elif gssha_card.name == '#INDEXGRID_GUID':
                                path_split = updated_value.split()
                                updated_path = os.path.basename(path_split[0].strip('"').strip("'"))
                                if os.path.exists(updated_path):
                                    new_path = os.path.join("..", os.path.basename(updated_path))
                                    try:
                                        # Get WMS ID for Index Map as part of value
                                        gssha_card.value = '"{0}" "{1}"'.format(new_path, path_split[1])
                                    except:
                                        # Like normal if the ID isn't there
                                        gssha_card.value = '"{0}"'.format(new_path)
                                else:
                                    log.warning("{0} {1} not found in project directory ...".format("#INDEXGRID_GUID", updated_path))

            # make sure project path is ""
            self._update_card("PROJECT_PATH", "", True)

            # WRITE OUT UPDATED GSSHA PROJECT FILE
            self.project_manager.write(session=self.db_session,
                                       directory=working_directory,
                                       name=self.project_manager.name)

            with tmp_chdir(working_directory):
                # RUN SIMULATION
                if self.gssha_executable and find_executable(self.gssha_executable) is not None:
                    log.info("Running GSSHA simulation ...")

                    try:
                        run_gssha_command = [self.gssha_executable,
                                             os.path.join(working_directory, self.project_filename)]
                        # run GSSHA
                        out = subprocess.check_output(run_gssha_command)

                        # write out GSSHA output
                        log_file_path = os.path.join(working_directory, 'simulation.log')
                        with open(log_file_path, mode='w') as logfile:
                            logfile.write(out.decode('utf-8'))
                            # log to other logger if debug mode on
                            if log.isEnabledFor(logging.DEBUG):
                                for line in out.split(b'\n'):
                                    log.debug(line.decode('utf-8'))

                    except subprocess.CalledProcessError as ex:
                        log.error("{0}: {1}".format(ex.returncode, ex.output))

                else:
                    missing_exe_error = ("GSSHA executable not found. "
                                         "Skipping GSSHA simulation run ...")
                    log.error(missing_exe_error)
                    raise ValueError(missing_exe_error)

            return working_directory
991339
def _prep_callable_bed(in_file, work_dir, stats, data):
    
    out_file = os.path.join(work_dir, "%s-merge.bed.gz" % utils.splitext_plus(os.path.basename(in_file))[0])
    gsort = config_utils.get_program("gsort", data)
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(data, out_file) as tx_out_file:
            fai_file = ref.fasta_idx(dd.get_ref_file(data))
            cmd = ("{gsort} {in_file} {fai_file} | bedtools merge -i - -d {stats[merge_size]} | "
                   "bgzip -c > {tx_out_file}")
            do.run(cmd.format(**locals()), "Prepare SV callable BED regions")
    return vcfutils.bgzip_and_index(out_file, data["config"])
182442
def _createGsshaPyObjects(self, eventChunk):
        
        ## TODO: Add Support for RADAR file format type values

        # Create GSSHAPY PrecipEvent
        event = PrecipEvent(description=eventChunk['description'],
                            nrGag=eventChunk['nrgag'],
                            nrPds=eventChunk['nrpds'])

        # Associate PrecipEvent with PrecipFile
        event.precipFile = self

        gages = []
        for coord in eventChunk['coords']:
            # Create GSSHAPY PrecipGage object
            gage = PrecipGage(description=coord['description'],
                              x=coord['x'],
                              y=coord['y'])

            # Associate PrecipGage with PrecipEvent
            gage.event = event

            # Append to gages list for association with PrecipValues
            gages.append(gage)

        for valLine in eventChunk['valLines']:
            for index, value in enumerate(valLine['values']):
                # Create GSSHAPY PrecipValue object
                val = PrecipValue(valueType=valLine['type'],
                                  dateTime=valLine['dateTime'],
                                  value=value)

                # Associate PrecipValue with PrecipEvent and PrecipGage
                val.event = event
                val.gage = gages[index]
706514
def run_forecast(self):

        
        # ----------------------------------------------------------------------
        # LSM to GSSHA
        # ----------------------------------------------------------------------
        self.prepare_hmet()
        self.prepare_gag()

        # ----------------------------------------------------------------------
        # RAPID to GSSHA
        # ----------------------------------------------------------------------
        self.rapid_to_gssha()

        # ----------------------------------------------------------------------
        # HOTSTART
        # ----------------------------------------------------------------------
        self.hotstart()

        # ----------------------------------------------------------------------
        # Run GSSHA
        # ----------------------------------------------------------------------
        return self.run()
1094509
def _pooling_output_shape(input_shape, pool_size=(2, 2),
                          strides=None, padding='VALID'):
  """"""
  dims = (1,) + pool_size + (1,)  # NHWC
  spatial_strides = strides or (1,) * len(pool_size)
  strides = (1,) + spatial_strides + (1,)
  pads = padtype_to_pads(input_shape, dims, strides, padding)
  operand_padded = onp.add(input_shape, onp.add(*zip(*pads)))
  t = onp.floor_divide(onp.subtract(operand_padded, dims), strides) + 1
  return tuple(t)
722648
def _apply_vqsr(in_file, ref_file, recal_file, tranch_file,
                sensitivity_cutoff, filter_type, data):
    
    base, ext = utils.splitext_plus(in_file)
    out_file = "{base}-{filter}filter{ext}".format(base=base, ext=ext,
                                                   filter=filter_type)
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            broad_runner = broad.runner_from_config(data["config"])
            gatk_type = broad_runner.gatk_type()
            if gatk_type == "gatk4":
                params = ["-T", "ApplyVQSR",
                          "--variant", in_file,
                          "--output", tx_out_file,
                          "--recal-file", recal_file,
                          "--tranches-file", tranch_file]
            else:
                params = ["-T", "ApplyRecalibration",
                          "--input", in_file,
                          "--out", tx_out_file,
                          "--recal_file", recal_file,
                          "--tranches_file", tranch_file]
            params += ["-R", ref_file,
                       "--mode", filter_type]
            resources = config_utils.get_resources("gatk_apply_recalibration", data["config"])
            opts = resources.get("options", [])
            if not opts:
                if gatk_type == "gatk4":
                    opts += ["--truth-sensitivity-filter-level", sensitivity_cutoff]
                else:
                    opts += ["--ts_filter_level", sensitivity_cutoff]
            params += opts
            broad_runner.run_gatk(params)
    return out_file
1003084
def createfastq(self):
        """"""
        # Initialise samplecount
        samplecount = 0
        # If the fastq destination folder is not provided, make the default value of :path/:miseqfoldername
        self.fastqdestination = self.fastqdestination if self.fastqdestination else self.path + self.miseqfoldername
        # Make the path
        make_path(self.fastqdestination)
        # Initialise variables for storing index information
        index = ''
        indexlength = int()
        # bcl2fastq requires an older version of the sample sheet, this recreates the required version
        # Create the new sample sheet
        with open('{}/SampleSheet_modified.csv'.format(self.fastqdestination), "w") as modifiedsamplesheet:
            # Write the required headings to the file
            modifiedsamplesheet.write(
                "FCID,Lane,SampleID,SampleRef,Index,Description,Control,Recipe,Operator,SampleProject\n")
            for strain in self.samples:
                # Create a combined index of index1-index2
                try:
                    strain.run.modifiedindex = '{}-{}'.format(strain.run.index, strain.run.index2)
                    indexlength = 16
                    index = 'I8,I8'
                except KeyError:
                    strain.run.modifiedindex = strain.run.index
                    indexlength = 6
                    index = 'I6'
                # The list of items to print to each line of the modified sample sheet
                printlist = [self.flowcell, '1', strain.name, str(strain.run.SampleNumber), strain.run.modifiedindex,
                             strain.run.Description, 'N', 'NA',
                             strain.run.InvestigatorName, self.projectname]
                modifiedsamplesheet.write('{}\n'.format(",".join(printlist)))
                samplecount += 1
        # Set :forward/reverse length to :header.forward/reverse length if the argument is not provided, or it's 'full',
        # otherwise  use the supplied argument
        self.forwardlength = self.metadata.header.forwardlength if self.forwardlength.lower()\
            == 'full' else self.forwardlength
        # Set :reverselength to :header.reverselength
        self.reverselength = self.metadata.header.reverselength if self.reverselength.lower() \
            == 'full' else self.reverselength
        # As the number of cycles required is the number of forward reads + the index(8) + the second index(8)
        # Also set the basemask variable as required
        if self.reverselength != '0':
            self.readsneeded = int(self.forwardlength) + int(self.reverselength) + indexlength
            basemask = "Y{}n*,{},Y{}n*".format(self.forwardlength, index, self.reverselength)
            nohup = "nohup make -j 16 > nohup.out"
        else:
            #  + 1
            self.readsneeded = int(self.forwardlength) + indexlength
            basemask = "Y{}n*,{},n*".format(self.forwardlength, index)
            nohup = "nohup make -j 16 r1 > nohup.out"
        # Handle plurality appropriately
        samples = 'samples' if samplecount > 1 else 'sample'
        number = 'are' if samplecount > 1 else 'is'
        printtime('There {} {} {} in this run. '
                  'Running fastq creating module with the following parameters:\n'
                  'MiSeqPath: {},\n'
                  'MiSeqFolder: {},\n'
                  'Fastq destination: {},\n'
                  'SampleSheet: {}'
                  .format(number, samplecount, samples, self.miseqpath, self.miseqfolder,
                          self.fastqdestination, '{}/SampleSheet_modified.csv'.format(self.fastqdestination)),
                  self.start)
        # Count the number of completed cycles in the run of interest
        cycles = glob('{}Data/Intensities/BaseCalls/L001/C*'.format(self.miseqfolder))
        while len(cycles) < self.readsneeded:
            printtime('Currently at {} cycles. Waiting until the MiSeq reaches cycle {}'.format(len(cycles),
                      self.readsneeded), self.start)
            sleep(1800)
            cycles = glob('{}Data/Intensities/BaseCalls/L001/C*'.format(self.miseqfolder))
        # configureBClToFastq requires :self.miseqfolder//Data/Intensities/BaseCalls/config.xml in order to work
        # When you download runs from BaseSpace, this file is not provided. There is an empty config.xml file that
        # can be populated with run-specific values and moved to the appropriate folder
        if not os.path.isfile('{}Data/Intensities/BaseCalls/config.xml'.format(self.miseqfolder)):
            self.configfilepopulator()
        # Define the bcl2fastq system call
        bclcall = "configureBclToFastq.pl --input-dir {}Data/Intensities/BaseCalls " \
                  "--output-dir {} --force --sample-sheet {}/SampleSheet_modified.csv " \
                  "--mismatches 1 --no-eamss --fastq-cluster-count 0 --compression none --use-bases-mask {}"\
            .format(self.miseqfolder, self.fastqdestination, self.fastqdestination, basemask)
        # Define the nohup system call
        nohupcall = "cd {} && {}".format(self.fastqdestination, nohup)
        # fnull = open(os.devnull, 'wb')
        if not os.path.isdir("{}/Project_{}".format(self.fastqdestination, self.projectname)):
            # Call configureBclToFastq.pl
            printtime('Running bcl2fastq', self.start)
            # Run the commands
            threadlock = threading.Lock()
            outstr = ''
            outerr = ''
            out, err = run_subprocess(bclcall)
            outstr += out
            outerr += out
            out, err = run_subprocess(nohupcall)
            outstr += out
            outerr += out
            # call(bclcall, shell=True, stdout=fnull, stderr=fnull)
            # call(nohupcall, shell=True, stdout=fnull, stderr=fnull)
            threadlock.acquire()
            write_to_logfile(bclcall, bclcall, self.logfile)
            write_to_logfile(nohupcall, nohupcall, self.logfile)
            write_to_logfile(outstr, outerr, self.logfile)
            threadlock.release()
        # Populate the metadata
        for sample in self.metadata.samples:
            sample.commands = GenObject()
            sample.commands.nohup = nohupcall
            sample.commands.bcl = bclcall
            sample.run.forwardlength = self.forwardlength
            sample.run.reverselength = self.reverselength
        # Copy the fastq files to a central folder so they can be processed
        self.fastqmover()
103302
function(T, SA, fs, n, k, isbwt) {
        var C, B, RA;
        var i, j, b, c, m, p, q, name, pidx = 0, newfs;
        var c0, c1;
        var flags = 0;

        // allocate temporary storage [CSA]
        if (k <= 256) {
            C = Util.makeS32Buffer(k);
            if (k <= fs) { B = SA.subarray(n + fs - k); flags = 1; }
            else { B = Util.makeS32Buffer(k); flags = 3; }
        } else if (k <= fs) {
            C = SA.subarray(n + fs - k);
            if (k <= (fs - k)) { B = SA.subarray(n + fs - k * 2); flags = 0; }
            else if (k <= 1024) { B = Util.makeS32Buffer(k); flags = 2; }
            else { B = C; flags = 8; }
        } else {
            C = B = Util.makeS32Buffer(k);
            flags = 4 | 8;
        }

        /* stage 1: reduce the problem by at least 1/2
           sort all the LMS-substrings */
        getCounts(T, C, n, k);
        getBuckets(C, B, k, true); /* find ends of buckets */
        for (i = 0; i < n; i++) { SA[i] = 0; }
        b = -1; i = n - 1; j = n; m = 0; c0 = T[n - 1];
        do { c1 = c0; } while ((--i >= 0) && ((c0 = T[i]) >= c1));
        for (; i >= 0 ;) {
            do { c1 = c0; } while ((--i >= 0) && ((c0 = T[i]) <= c1));
            if ( i >= 0 ) {
                if ( b >= 0 ) { SA[b] = j; }
                b = --B[c1];
                j = i;
                ++m;
                do { c1 = c0; } while ((--i >= 0) && ((c0 = T[i]) >= c1));
            }
        }

        if (m > 1) {
            LMSsort(T, SA, C, B, n, k);
            name = LMSpostproc(T, SA, n, m);
        } else if (m === 1) {
            SA[b] = j + 1;
            name = 1;
        } else {
            name = 0;
        }

        /* stage 2: solve the reduced problem
           recurse if names are not yet unique */
        if(name < m) {
            if((flags & 4) !== 0) { C = null; B = null; }
            if((flags & 2) !== 0) { B = null; }
            newfs = (n + fs) - (m * 2);
            if((flags & (1 | 4 | 8)) === 0) {
                if((k + name) <= newfs) { newfs -= k; }
                else { flags |= 8; }
            }
            ASSERT( (n >>> 1) <= (newfs + m) );
            for (i = m + (n >>> 1) - 1, j = m * 2 + newfs - 1; m <= i; i--) {
                if(SA[i] !== 0) { SA[j--] = SA[i] - 1; }
            }
            RA = SA.subarray(m + newfs);
            SA_IS(RA, SA, newfs, m, name, false);
            RA = null;

            i = n - 1; j = m * 2 - 1; c0 = T[n - 1];
            do { c1 = c0; } while ((--i >= 0) && ((c0 = T[i]) >= c1));
            for (; i >= 0 ;) {
                do { c1 = c0; } while ((--i >= 0) && ((c0 = T[i]) <= c1));
                if ( i >= 0 ) {
                    SA[j--] = i + 1;
                    do { c1 = c0; } while ((--i >= 0) && ((c0 = T[i]) >= c1));
                }
            }

            for (i = 0; i < m; i++) { SA[i] = SA[m + SA[i]]; }
            if((flags & 4) !== 0) { C = B = Util.makeS32Buffer(k); }
            if((flags & 2) !== 0) { B = Util.makeS32Buffer(k); }
        }

        /* stage 3: induce the result for the original problem */
        if((flags & 8) !== 0) { getCounts(T, C, n, k); }
        /* put all left-most S characters into their buckets */
        if (m > 1) {
            getBuckets(C, B, k, true); /* find ends of buckets */
            i = m - 1; j = n; p = SA[m - 1]; c1 = T[p];
            do {
                q = B[c0 = c1];
                while (q < j) { SA[--j] = 0; }
                do {
                    SA[--j] = p;
                    if(--i < 0) { break; }
                    p = SA[i];
                } while((c1 = T[p]) === c0);
            } while (i >= 0 );
            while ( j > 0 ) { SA[--j] = 0; }
        }
        if (!isbwt) { induceSA(T, SA, C, B, n, k); }
        else { pidx = computeBWT(T, SA, C, B, n, k); }
        C = null; B = null;
        return pidx;
    }
294726
def prepare(args):
    
    from jcvi.formats.base import write_file

    p = OptionParser(prepare.__doc__ + FastqNamings)
    p.add_option("-K", default=45, type="int",
                 help="K-mer size [default: %default]")
    p.add_option("--assemble_1st_rank_only", default=False, action="store_true",
                 help="Assemble the first rank only, other libs asm_flags=2 [default: %default]")
    p.add_option("--scaffold",
                 help="Only perform scaffolding [default: %default]")
    p.add_option("--gapclose",
                 help="Only perform gap closure [default: %default]")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    fnames = args
    K = opts.K
    for x in fnames:
        assert op.exists(x), "File `{0}` not found.".format(x)

    a1st = opts.assemble_1st_rank_only

    cfgfile = "soap.config"
    gc_cfgfile = "soap.gc.config"
    fw = open(cfgfile, "w")
    fw_gc = open(gc_cfgfile, "w")

    libs = get_libs(fnames)
    rank = 0
    singletons = []
    max_rd_len = max(readlen([f]) for f in fnames)

    block = "max_rd_len={0}\n".format(max_rd_len)
    for stream in (sys.stderr, fw, fw_gc):
        print(block, file=stream)

    # Collect singletons first
    singletons = []
    for lib, fs in libs:
        if lib.size == 0:
            singletons += fs
            continue

    for lib, fs in libs:
        size = lib.size
        if size == 0:
            continue

        rank += 1
        block = "[LIB]\n"
        block += "avg_ins={0}\n".format(size)
        f = fs[0]
        block += "reverse_seq={0}\n".format(lib.reverse_seq)
        asm_flags = 2 if (rank > 1 and a1st) else lib.asm_flags
        block += "asm_flags={0}\n".format(asm_flags)
        block += "rank={0}\n".format(rank)
        if lib.reverse_seq:
            pair_num_cutoff = 3
            block += "pair_num_cutoff={0}\n".format(pair_num_cutoff)
        block += "map_len=35\n"

        for f in fs:
            if ".1." in f:
                tag = "q1"
            elif ".2." in f:
                tag = "q2"
            block += "{0}={1}\n".format(tag, f)

        if rank == 1:
            for s in singletons:
                tag = "q" if is_fastq(s) else "f"
                block += tag + "={0}\n".format(s)

        print(block, file=sys.stderr)
        print(block, file=fw)

        if asm_flags > 2:
            print(block, file=fw_gc)

    runfile = "run.sh"
    scaffold = opts.scaffold
    bb = 63 if K <= 63 else 127
    binary = "SOAPdenovo-{0}mer".format(bb)
    header = SOAPHEADER.format(opts.cpus, K, binary)
    if opts.gapclose:
        gapclose = opts.gapclose
        outfile = gapclose.rsplit(".", 1)[0] + ".closed.fasta"
        template = header + GCRUNG.format(gapclose, outfile)
    else:
        template = header + (SCFRUN % scaffold if scaffold else SOAPRUN)

    write_file(runfile, template)
    fw.close()
    fw_gc.close()
345429
def write_prepareref_fasta_file(outfile, gene_coords, genes_need_upstream, genes_non_upstream, upstream_before=100, upstream_after=100):
    ''''''
    tmp_dict = {}
    fasta_in = os.path.join(data_dir, 'NC_000962.3.fa.gz')
    pyfastaq.tasks.file_to_dict(fasta_in, tmp_dict)
    ref_seq = tmp_dict['NC_000962.3']

    with open(outfile, 'w') as f:
        for gene in genes_non_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start:end+1])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[end:start+1])
                gene_fa.revcomp()

            print(gene_fa, file=f)

        for gene in genes_need_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_before:start + upstream_after])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_after + 1:start + upstream_before + 1])
                gene_fa.revcomp()

            gene_fa.id += '_upstream'
            print(gene_fa, file=f)
766300
def re_pipe(FlowRate, Diam, Nu):
    """"""
    #Checking input validity
    ut.check_range([FlowRate, ">0", "Flow rate"], [Diam, ">0", "Diameter"],
                   [Nu, ">0", "Nu"])
    return (4 * FlowRate) / (np.pi * Diam * Nu)
267366
def precompute(self, cache_dir=None, swath_usage=0, **kwargs):
        """"""
        if kwargs.get('mask') is not None:
            LOG.warning("'mask' parameter has no affect during EWA "
                        "resampling")

        del kwargs
        source_geo_def = self.source_geo_def
        target_geo_def = self.target_geo_def

        if cache_dir:
            LOG.warning("'cache_dir' is not used by EWA resampling")

        # Satpy/PyResample don't support dynamic grids out of the box yet
        lons, lats = source_geo_def.get_lonlats()
        if isinstance(lons, xr.DataArray):
            # get dask arrays
            lons = lons.data
            lats = lats.data
        # we are remapping to a static unchanging grid/area with all of
        # its parameters specified
        chunks = (2,) + lons.chunks
        res = da.map_blocks(self._call_ll2cr, lons, lats,
                            target_geo_def, swath_usage,
                            dtype=lons.dtype, chunks=chunks, new_axis=[0])
        cols = res[0]
        rows = res[1]

        # save the dask arrays in the class instance cache
        # the on-disk cache will store the numpy arrays
        self.cache = {
            "rows": rows,
            "cols": cols,
        }

        return None
89416
@SuppressWarnings("deprecation")
  public static void writeFlowDetails(JsonGenerator jsonGenerator, Flow aFlow,
      DetailLevel selectedSerialization, Predicate<String> includeFilter)
      throws JsonGenerationException, IOException {
    jsonGenerator.writeStartObject();
    // serialize the FlowKey object
    filteredWrite("flowKey", includeFilter, aFlow.getFlowKey(), jsonGenerator);

    // serialize individual members of this class
    filteredWrite("flowName", includeFilter, aFlow.getFlowName(), jsonGenerator);
    filteredWrite("userName", includeFilter, aFlow.getUserName(), jsonGenerator);
    filteredWrite("jobCount", includeFilter, aFlow.getJobCount(), jsonGenerator);
    filteredWrite("totalMaps", includeFilter, aFlow.getTotalMaps(), jsonGenerator);
    filteredWrite("totalReduces", includeFilter, aFlow.getTotalReduces(), jsonGenerator);
    filteredWrite("mapFileBytesRead", includeFilter, aFlow.getMapFileBytesRead(), jsonGenerator);
    filteredWrite("mapFileBytesWritten", includeFilter, aFlow.getMapFileBytesWritten(), jsonGenerator);
    filteredWrite("reduceFileBytesRead", includeFilter, aFlow.getReduceFileBytesRead(), jsonGenerator);
    filteredWrite("hdfsBytesRead", includeFilter, aFlow.getHdfsBytesRead(), jsonGenerator);
    filteredWrite("hdfsBytesWritten", includeFilter, aFlow.getHdfsBytesWritten(), jsonGenerator);
    filteredWrite("mapSlotMillis", includeFilter, aFlow.getMapSlotMillis(), jsonGenerator);
    filteredWrite("reduceSlotMillis", includeFilter, aFlow.getReduceSlotMillis(), jsonGenerator);
    filteredWrite("megabyteMillis", includeFilter, aFlow.getMegabyteMillis(), jsonGenerator);
    filteredWrite("cost", includeFilter, aFlow.getCost(), jsonGenerator);
    filteredWrite("reduceShuffleBytes", includeFilter, aFlow.getReduceShuffleBytes(), jsonGenerator);
    filteredWrite("duration", includeFilter, aFlow.getDuration(), jsonGenerator);
    filteredWrite("wallClockTime", includeFilter, aFlow.getWallClockTime(), jsonGenerator);
    filteredWrite("cluster", includeFilter, aFlow.getCluster(), jsonGenerator);
    filteredWrite("appId", includeFilter, aFlow.getAppId(), jsonGenerator);
    filteredWrite("runId", includeFilter, aFlow.getRunId(), jsonGenerator);
    filteredWrite("version", includeFilter, aFlow.getVersion(), jsonGenerator);
    filteredWrite("hadoopVersion", includeFilter, aFlow.getHadoopVersion(), jsonGenerator);
    if (selectedSerialization == SerializationContext.DetailLevel.EVERYTHING) {
      filteredWrite("submitTime", includeFilter, aFlow.getSubmitTime(), jsonGenerator);
      filteredWrite("launchTime", includeFilter, aFlow.getLaunchTime(), jsonGenerator);
      filteredWrite("finishTime", includeFilter, aFlow.getFinishTime(), jsonGenerator);
    }
    filteredWrite(Constants.HRAVEN_QUEUE, includeFilter, aFlow.getQueue(), jsonGenerator);
    filteredWrite("counters", includeFilter, aFlow.getCounters(), jsonGenerator);
    filteredWrite("mapCounters", includeFilter, aFlow.getMapCounters(), jsonGenerator);
    filteredWrite("reduceCounters", includeFilter, aFlow.getReduceCounters(), jsonGenerator);

    // if flag, include job details
    if ((selectedSerialization == SerializationContext.DetailLevel.FLOW_SUMMARY_STATS_WITH_JOB_STATS)
        || (selectedSerialization == SerializationContext.DetailLevel.EVERYTHING)) {
      jsonGenerator.writeFieldName("jobs");
      jsonGenerator.writeObject(aFlow.getJobs());
    }
    jsonGenerator.writeEndObject();

  }
375986
def execute(self, nc_file_list, index_list, in_weight_table,
                out_nc, grid_type, mp_lock):

        """"""
        if not os.path.exists(out_nc):
            raise Exception("Outfile has not been created. "
                            "You need to run: generateOutputInflowFile "
                            "function ...")

        if len(nc_file_list) != len(index_list):
            raise Exception("ERROR: Number of runoff files not equal to "
                            "number of indices ...")

        demo_file_list = nc_file_list[0]
        if not isinstance(nc_file_list[0], list):
            demo_file_list = [demo_file_list]

        self.data_validation(demo_file_list[0])
        self.read_in_weight_table(in_weight_table)

        conversion_factor = self.get_conversion_factor(demo_file_list[0],
                                                       len(demo_file_list))

        # get indices of subset of data
        lon_ind_all = [int(i) for i in self.dict_list[self.header_wt[2]]]
        lat_ind_all = [int(j) for j in self.dict_list[self.header_wt[3]]]

        # Obtain a subset of  runoff data based on the indices in the
        # weight table
        min_lon_ind_all = min(lon_ind_all)
        max_lon_ind_all = max(lon_ind_all)
        min_lat_ind_all = min(lat_ind_all)
        max_lat_ind_all = max(lat_ind_all)
        lon_slice = slice(min_lon_ind_all, max_lon_ind_all + 1)
        lat_slice = slice(min_lat_ind_all, max_lat_ind_all + 1)
        index_new = []

        # combine inflow data
        for nc_file_array_index, nc_file_array in enumerate(nc_file_list):

            index = index_list[nc_file_array_index]

            if not isinstance(nc_file_array, list):
                nc_file_array = [nc_file_array]

            data_subset_all = None
            for nc_file in nc_file_array:
                # Validate the netcdf dataset
                self.data_validation(nc_file)

                # Read the netcdf dataset
                data_in_nc = Dataset(nc_file)

                # Calculate water inflows
                runoff_dimension_size = \
                    len(data_in_nc.variables[self.runoff_vars[0]].dimensions)
                if runoff_dimension_size == 2:
                    # obtain subset of surface and subsurface runoff
                    data_subset_runoff = \
                        data_in_nc.variables[self.runoff_vars[0]][
                            lat_slice, lon_slice]
                    for var_name in self.runoff_vars[1:]:
                        data_subset_runoff += \
                            data_in_nc.variables[var_name][
                                lat_slice, lon_slice]

                    # get runoff dims
                    len_time_subset = 1
                    len_lat_subset = data_subset_runoff.shape[0]
                    len_lon_subset = data_subset_runoff.shape[1]

                    # reshape the runoff
                    data_subset_runoff = data_subset_runoff.reshape(
                        len_lat_subset * len_lon_subset)

                elif runoff_dimension_size == 3:
                    # obtain subset of surface and subsurface runoff
                    data_subset_runoff = \
                        data_in_nc.variables[self.runoff_vars[0]][
                            :, lat_slice, lon_slice]
                    for var_name in self.runoff_vars[1:]:
                        data_subset_runoff += \
                            data_in_nc.variables[var_name][
                                :, lat_slice, lon_slice]

                    # get runoff dims
                    len_time_subset = data_subset_runoff.shape[0]
                    len_lat_subset = data_subset_runoff.shape[1]
                    len_lon_subset = data_subset_runoff.shape[2]
                    # reshape the runoff
                    data_subset_runoff = \
                        data_subset_runoff.reshape(
                            len_time_subset,
                            (len_lat_subset * len_lon_subset))

                data_in_nc.close()

                if not index_new:
                    # compute new indices based on the data_subset_surface
                    for r in range(0, self.count):
                        ind_lat_orig = lat_ind_all[r]
                        ind_lon_orig = lon_ind_all[r]
                        index_new.append(
                            (ind_lat_orig - min_lat_ind_all) * len_lon_subset
                            + (ind_lon_orig - min_lon_ind_all))

                # obtain a new subset of data
                if runoff_dimension_size == 2:
                    data_subset_new = data_subset_runoff[index_new]
                elif runoff_dimension_size == 3:
                    data_subset_new = data_subset_runoff[:, index_new]

                # FILTER DATA
                try:
                    # set masked values to zero
                    data_subset_new = data_subset_new.filled(fill_value=0)
                except AttributeError:
                    pass
                # set negative values to zero
                data_subset_new[data_subset_new < 0] = 0

                # combine data
                if data_subset_all is None:
                    data_subset_all = data_subset_new
                else:
                    data_subset_all = np.add(data_subset_all, data_subset_new)

            if runoff_dimension_size == 3 and len_time_subset > 1:
                inflow_data = np.zeros((len_time_subset, self.size_stream_id))
            else:
                inflow_data = np.zeros(self.size_stream_id)

            pointer = 0
            for stream_index in xrange(self.size_stream_id):
                npoints = int(self.dict_list[self.header_wt[4]][pointer])
                # Check if all npoints points correspond to the same streamID
                if len(set(self.dict_list[self.header_wt[0]][
                           pointer: (pointer + npoints)])) != 1:
                    print("ROW INDEX {0}".format(pointer))
                    print("COMID {0}".format(
                        self.dict_list[self.header_wt[0]][pointer]))
                    raise Exception(self.error_messages[2])

                area_sqm_npoints = \
                    np.array([float(k) for k in
                              self.dict_list[self.header_wt[1]][
                              pointer: (pointer + npoints)]])

                # assume data is incremental
                if runoff_dimension_size == 3:
                    data_goal = data_subset_all[:, pointer:(pointer + npoints)]
                else:
                    data_goal = data_subset_all[pointer:(pointer + npoints)]

                if grid_type == 't255':
                    # A) ERA Interim Low Res (T255) - data is cumulative
                    # from time 3/6/9/12
                    # (time zero not included, so assumed to be zero)
                    ro_first_half = \
                        np.concatenate([data_goal[0:1, ],
                                        np.subtract(data_goal[1:4, ],
                                                    data_goal[0:3, ])])
                    # from time 15/18/21/24
                    # (time restarts at time 12, assumed to be zero)
                    ro_second_half = \
                        np.concatenate([data_goal[4:5, ],
                                        np.subtract(data_goal[5:, ],
                                                    data_goal[4:7, ])])
                    ro_stream = \
                        np.multiply(
                            np.concatenate([ro_first_half, ro_second_half]),
                            area_sqm_npoints)

                else:
                    ro_stream = data_goal * area_sqm_npoints * \
                                conversion_factor

                # filter nan
                ro_stream[np.isnan(ro_stream)] = 0

                if ro_stream.any():
                    if runoff_dimension_size == 3 and len_time_subset > 1:
                        inflow_data[:, stream_index] = ro_stream.sum(axis=1)
                    else:
                        inflow_data[stream_index] = ro_stream.sum()

                pointer += npoints

            # only one process is allowed to write at a time to netcdf file
            mp_lock.acquire()
            data_out_nc = Dataset(out_nc, "a", format="NETCDF3_CLASSIC")
            if runoff_dimension_size == 3 and len_time_subset > 1:
                data_out_nc.variables['m3_riv'][
                    index*len_time_subset:(index+1)*len_time_subset, :] = \
                    inflow_data
            else:
                data_out_nc.variables['m3_riv'][index] = inflow_data
            data_out_nc.close()
            mp_lock.release()
1003255
private void prepareThreadManager(StreamConfiguration sc) {
    assert(!flacStreamIsOpen);
    if(MAX_THREADED_FRAMES > 0) {
      threadManager = new BlockThreadManager(this);
      threadedFrames = new Frame[MAX_THREADED_FRAMES];
      for(int i = 0; i < MAX_THREADED_FRAMES; i++) {
        threadedFrames[i]  = new Frame(this.streamConfig);
        threadManager.addFrameThread(threadedFrames[i]);
      }
    }
    else {
      threadManager = null;
    }
  }
1141350
def main():
    '''
    
    '''
    args = get_args()
    try:
        utils.make_output_dir(args.outdir)
        utils.init_logs(args)
        args.format = nanoplotter.check_valid_format(args.format)
        settings = vars(args)
        settings["path"] = path.join(args.outdir, args.prefix)
        sources = {
            "fastq": args.fastq,
            "bam": args.bam,
            "cram": args.cram,
            "fastq_rich": args.fastq_rich,
            "fastq_minimal": args.fastq_minimal,
            "summary": args.summary,
            "fasta": args.fasta,
            "ubam": args.ubam,
        }

        if args.pickle:
            datadf = pickle.load(open(args.pickle, 'rb'))
        else:
            datadf = get_input(
                source=[n for n, s in sources.items() if s][0],
                files=[f for f in sources.values() if f][0],
                threads=args.threads,
                readtype=args.readtype,
                combine="simple",
                barcoded=args.barcoded)
        if args.store:
            pickle.dump(
                obj=datadf,
                file=open(settings["path"] + "NanoPlot-data.pickle", 'wb'))
        if args.raw:
            datadf.to_csv(settings["path"] + "NanoPlot-data.tsv.gz",
                          sep="\t",
                          index=False,
                          compression="gzip")

        settings["statsfile"] = [make_stats(datadf, settings, suffix="")]
        datadf, settings = filter_and_transform_data(datadf, settings)
        if settings["filtered"]:  # Bool set when filter was applied in filter_and_transform_data()
            settings["statsfile"].append(make_stats(datadf, settings, suffix="_post_filtering"))

        if args.barcoded:
            barcodes = list(datadf["barcode"].unique())
            plots = []
            for barc in barcodes:
                logging.info("Processing {}".format(barc))
                settings["path"] = path.join(args.outdir, args.prefix + barc + "_")
                dfbarc = datadf[datadf["barcode"] == barc]
                if len(dfbarc) > 5:
                    settings["title"] = barc
                    plots.extend(
                        make_plots(dfbarc, settings)
                    )
                else:
                    sys.stderr.write("Found barcode {} less than 5x, ignoring...\n".format(barc))
                    logging.info("Found barcode {} less than 5 times, ignoring".format(barc))
        else:
            plots = make_plots(datadf, settings)
        make_report(plots, settings)
        logging.info("Finished!")
    except Exception as e:
        logging.error(e, exc_info=True)
        print("\n\n\nIf you read this then NanoPlot has crashed :-(")
        print("Please try updating NanoPlot and see if that helps...\n")
        print("If not, please report this issue at https://github.com/wdecoster/NanoPlot/issues")
        print("If you could include the log file that would be really helpful.")
        print("Thanks!\n\n\n")
        raise