text:
# -*- coding: utf-8 -*-
import re
import time
import os
import errno
import numpy as np
import zipfile
import subprocess
from pprint import pprint
import uuid
import copy
import json
import psutil

from installed_clients.WorkspaceClient import Workspace
from installed_clients.KBaseReportClient import KBaseReport
from installed_clients.AssemblyUtilClient import AssemblyUtil
from installed_clients.kb_quastClient import kb_quast
from installed_clients.ReadsUtilsClient import ReadsUtils
from installed_clients.baseclient import ServerError

def log(message, prefix_newline=False):
    """Logging function, provides a hook to suppress or redirect log messages."""
    print(('\n' if prefix_newline else '') + '{0:.2f}'.format(time.time()) + ': ' + str(message))

def _mkdir_p(path):
    """
    _mkdir_p: make directory for given path
    """
    if not path:
        return
    try:
        os.makedirs(path)
    except OSError as exc:
        if exc.errno == errno.EEXIST and os.path.isdir(path):
            pass
        else:
            raise

class SPAdesUtils:
    """
    Define the SPAdesUtils functions
    """
    SPADES_VERSION = '3.13.0'
    SPADES_BIN = '/opt/SPAdes-' + SPADES_VERSION + '-Linux/bin'

    DISABLE_SPADES_OUTPUT = False  # should be False in production

    # Basic options
    PARAM_IN_SINGLE_CELL = 'single_cell'  # --sc
    PARAM_IN_METAGENOME = 'metagenomic'  # --meta
    PARAM_IN_PLASMID = 'plasmid'  # --plasmid
    PARAM_IN_RNA = 'rna'  # --rna
    PARAM_IN_IONTORRENT = 'iontorrent'  # --iontorrent

    # Pipeline options
    PARAM_IN_ONLY_ERROR_CORR = 'only-error-correction'  # --only-error-correction
    PARAM_IN_ONLY_ASSEMBLER = 'only-assembler'  # --only-assembler
    PARAM_IN_CAREFUL = 'careful'  # --careful
    PARAM_IN_CONTINUE = 'continue'  # --continue
    PARAM_IN_DISABLE_GZIP = 'disable-gzip-output'  # --disable-gzip-output

    # Input parameters
    PARAM_IN_WS = 'workspace_name'
    PARAM_IN_CS_NAME = 'output_contigset_name'
    PARAM_IN_READS = 'reads_libraries'
    PARAM_IN_LONG_READS = 'long_reads_libraries'
    PARAM_IN_KMER_SIZES = 'kmer_sizes'
    PARAM_IN_SKIP_ERR_CORRECT = 'skip_error_correction'
    PARAM_IN_MIN_CONTIG_LENGTH = 'min_contig_length'
    PARAM_IN_DNA_SOURCE = 'dna_source'
    PARAM_IN_PIPELINE_OPTION = 'pipeline_options'
    ASSEMBLE_RESULTS_DIR = 'assemble_results'

    INVALID_WS_OBJ_NAME_RE = re.compile('[^\\w\\|._-]')
    INVALID_WS_NAME_RE = re.compile('[^\\w:._-]')

    THREADS_PER_CORE = 3
    MAX_THREADS = 64  # per email thread with Anton Korobeynikov
    MAX_THREADS_META = 128  # Increase threads for metagenomic assemblies
    MEMORY_OFFSET_GB = 1  # 1GB
    MIN_MEMORY_GB = 5
    MAX_MEMORY_GB_SPADES = 500
    MAX_MEMORY_GB_META_SPADES = 1000
    GB = 1000000000

    # private method definition
    def __init__(self, prj_dir, config):
        self.workspace_url = config['workspace-url']
        self.callback_url = config['SDK_CALLBACK_URL']
        self.token = config['KB_AUTH_TOKEN']
        if 'shock-url' in config:
            self.shock_url = config['shock-url']
        if 'handle-service-url' in config:
            self.handle_url = config['handle-service-url']

        self.ws_client = Workspace(self.workspace_url, token=self.token)
        self.ru = ReadsUtils(self.callback_url, token=self.token, service_ver='release')
        self.au = AssemblyUtil(self.callback_url, token=self.token, service_ver='release')
        self.kbr = KBaseReport(self.callback_url)
        self.kbq = kb_quast(self.callback_url)
        self.proj_dir = prj_dir

        self.spades_version = 'SPAdes-' + os.environ['SPADES_VERSION']

    def _get_kbreads_info(self, wsname, reads_refs):
        """
        _get_kbreads_info--from a set of given KBase reads refs, fetches the corresponding
        reads info with as interleaved fastq files and returns a list of reads data in the
        following structure:
        reads_data = {
                'fwd_file': path_to_fastq_file,
                'type': reads_type,  # ('interleaved', 'paired', or 'single')
                'seq_tech': sequencing_tech,
                'reads_ref': KBase object ref for downstream convenience,
                'reads_name': KBase object name for downstream convenience,
                'rev_file': path_to_fastq_file,  # only if paired end
        }
        """
        obj_ids = []
        for r in reads_refs:
            if r:
                obj_ids.append({'ref': r if '/' in r else (wsname + '/' + r)})

        if not obj_ids:
            return []

        ws_info = self.ws_client.get_object_info_new({'objects': obj_ids})
        reads_params = []

        reftoname = {}
        for wsi, oid in zip(ws_info, obj_ids):
            ref = oid['ref']
            reads_params.append(ref)
            obj_name = wsi[1]
            reftoname[ref] = wsi[7] + '/' + obj_name

        typeerr = ('Supported types: KBaseFile.SingleEndLibrary ' +
                   'KBaseFile.PairedEndLibrary ' +
                   'KBaseAssembly.SingleEndLibrary ' +
                   'KBaseAssembly.PairedEndLibrary')
        try:
            reads = self.ru.download_reads({
                        'read_libraries': reads_params,
                        'interleaved': 'false'
                        })['files']
        except ServerError as se:
            log('logging stacktrace from dynamic client error')
            log(se.data)
            if typeerr in se.message:
                prefix = se.message.split('.')[0]
                raise ValueError(
                    prefix + '. Only the types ' +
                    'KBaseAssembly.SingleEndLibrary ' +
                    'KBaseAssembly.PairedEndLibrary ' +
                    'KBaseFile.SingleEndLibrary ' +
                    'and KBaseFile.PairedEndLibrary are supported')
            else:
                raise

        # log('Downloaded reads data from KBase:\n' + pformat(reads))
        reads_data = []
        for ref in reads_refs:
            reads_name = reftoname[ref]
            f = reads[ref]['files']
            seq_tech = reads[ref]['sequencing_tech']
            rds_info = {
                'fwd_file': f['fwd'],
                'reads_ref': ref,
                'type': f['type'],
                'seq_tech': seq_tech,
                'reads_name': reads_name
            }
            if f.get('rev', None):
                rds_info['rev_file'] = f['rev']
            reads_data.append(rds_info)

        return reads_data

    def _generate_output_file_list(self, out_dir):
        """
        _generate_output_file_list: zip result files and generate file_links for report
        """
        log('start packing result files')

        output_files = list()

        output_directory = os.path.join(self.proj_dir, str(uuid.uuid4()))
        _mkdir_p(output_directory)
        spades_output = os.path.join(output_directory, 'spades_output.zip')
        self._zip_folder(out_dir, spades_output)

        output_files.append({'path': spades_output,
                             'name': os.path.basename(spades_output),
                             'label': os.path.basename(spades_output),
                             'description': 'Output file(s) generated by {}'.format(
                                 self.spades_version)})

        return output_files

    def _zip_folder(self, folder_path, output_path):
        """
        _zip_folder: Zip the contents of an entire folder (with that folder included
        in the archive). Empty subfolders could be included in the archive as well
        if the commented portion is used.
        """
        with zipfile.ZipFile(output_path, 'w',
                             zipfile.ZIP_DEFLATED,
                             allowZip64=True) as ziph:
            for root, folders, files in os.walk(folder_path):
                for f in files:
                    absolute_path = os.path.join(root, f)
                    relative_path = os.path.join(os.path.basename(root), f)
                    # print "Adding {} to archive.".format(absolute_path)
                    ziph.write(absolute_path, relative_path)

        print("{} created successfully.".format(output_path))
        # with zipfile.ZipFile(output_path, "r") as f:
        #    print 'Checking the zipped file......\n'
        #    for info in f.infolist():
        #        print info.filename, info.date_time, info.file_size, info.compress_size

    def _load_stats(self, input_file_name):
        log('Starting conversion of FASTA to KBaseGenomeAnnotations.Assembly')
        log('Building Object.')
        if not os.path.isfile(input_file_name):
            raise Exception('The input file name {0} is not a file!'.format(input_file_name))
        with open(input_file_name, 'r') as input_file_handle:
            contig_id = None
            sequence_len = 0
            fasta_dict = dict()
            first_header_found = False
            # Pattern for replacing white space
            pattern = re.compile(r'\s+')
            for current_line in input_file_handle:
                if (current_line[0] == '>'):
                    # found a header line
                    # Wrap up previous fasta sequence
                    if not first_header_found:
                        first_header_found = True
                    else:
                        fasta_dict[contig_id] = sequence_len
                        sequence_len = 0
                    fasta_header = current_line.replace('>', '').strip()
                    try:
                        contig_id = fasta_header.strip().split(' ', 1)[0]
                    except (IndexError, ValueError, KeyError):
                        contig_id = fasta_header.strip()
                else:
                    sequence_len += len(re.sub(pattern, '', current_line))
        # wrap up last fasta sequence
        if not first_header_found:
            raise Exception("There are no contigs in this file")
        else:
            fasta_dict[contig_id] = sequence_len
        return fasta_dict

    def _parse_single_reads(self, reads_type, reads_list):
        """
        _parse_single_reads: given the reads_type and a list of reads, return an object
        defining the type and a list of fastq files.
        """
        single_reads_fqs = []
        ret_obj = {}
        if reads_list and isinstance(reads_list, list):
            for rds in reads_list:
                single_reads_fqs.append(rds['fwd_file'])
        if single_reads_fqs:
            ret_obj = {
                "type": reads_type,
                "single reads": single_reads_fqs
            }

        return ret_obj

    def _parse_pair_reads(self, reads_type, reads_list):
        """
        _parse_pair_reads: given the reads_type and a list of reads, return an object
        defining the type and a list of fastq files.
        """
        right_reads_fqs = []
        left_reads_fqs = []
        ret_obj = {}
        if reads_list and isinstance(reads_list, list):
            for rds in reads_list:
                right_reads_fqs.append(rds['fwd_file'])
                if rds.get('rev_file', None):
                    left_reads_fqs.append(rds['rev_file'])
            orent = reads_list[0]['orientation']

        if right_reads_fqs:
            ret_obj["right reads"] = right_reads_fqs
            ret_obj["orientation"] = orent
            ret_obj["type"] = reads_type
        if left_reads_fqs:
            ret_obj["left reads"] = left_reads_fqs

        return ret_obj
    # end of private methods

    # public method definitions

    def check_spades_params(self, params):
        """
        check_spades_params: checks params passed to run_HybridSPAdes method and set default values
        """
        # log('Start validating run_HybridSPAdes parameters:\n{}'.format(
        # json.dumps(params, indent=1)))

        # check for mandatory parameters
        if params.get(self.PARAM_IN_WS, None) is None:
            raise ValueError('Parameter {} is mandatory!'.format(self.PARAM_IN_WS))
        if self.INVALID_WS_NAME_RE.search(params[self.PARAM_IN_WS]):
            raise ValueError('Invalid workspace name: {}.'.format(params[self.PARAM_IN_WS]))

        if params.get(self.PARAM_IN_CS_NAME, None) is None:
            raise ValueError('Parameter {} is mandatory!'.format(self.PARAM_IN_CS_NAME))
        if self.INVALID_WS_OBJ_NAME_RE.search(params[self.PARAM_IN_CS_NAME]):
            raise ValueError('Invalid workspace object name: {}.'.format(
                params[self.PARAM_IN_CS_NAME]))

        if params.get(self.PARAM_IN_READS, None) is None:
            raise ValueError('Parameter {} is mandatory!'.format(self.PARAM_IN_READS))
        if type(params[self.PARAM_IN_READS]) != list:
            raise ValueError('Input reads {} must be a list.'.format(self.PARAM_IN_READS))
        if len(params[self.PARAM_IN_READS]) == 0:
            raise ValueError('Input parameter {} should have at least one reads.'.format(
                             self.PARAM_IN_READS))

        if self.PARAM_IN_MIN_CONTIG_LENGTH in params:
            if not isinstance(params[self.PARAM_IN_MIN_CONTIG_LENGTH], int):
                raise ValueError('{} must be of type int.'.format(self.PARAM_IN_MIN_CONTIG_LENGTH))

        if not params.get(self.PARAM_IN_KMER_SIZES, None):
            params[self.PARAM_IN_KMER_SIZES] = [21, 33, 55]
        kmer_sstr = ",".join(str(num) for num in params[self.PARAM_IN_KMER_SIZES])
        params[self.PARAM_IN_KMER_SIZES] = kmer_sstr
        print("KMER_SIZES: " + kmer_sstr)

        if params.get(self.PARAM_IN_SKIP_ERR_CORRECT, None):
            print("SKIP ERR CORRECTION: " + str(params[self.PARAM_IN_SKIP_ERR_CORRECT]))

        # check for basic option parameters
        if params.get(self.PARAM_IN_DNA_SOURCE, None):
            dna_src = params[self.PARAM_IN_DNA_SOURCE]
            if dna_src not in [self.PARAM_IN_SINGLE_CELL,
                               self.PARAM_IN_METAGENOME,
                               self.PARAM_IN_PLASMID,
                               self.PARAM_IN_RNA,
                               self.PARAM_IN_IONTORRENT]:
                params[self.PARAM_IN_DNA_SOURCE] = None
        else:
            params[self.PARAM_IN_DNA_SOURCE] = None

        # a list of basic options0
        params['basic_options'] = ['-o', self.ASSEMBLE_RESULTS_DIR]
        dna_src = params.get(self.PARAM_IN_DNA_SOURCE)
        if dna_src == self.PARAM_IN_SINGLE_CELL:
            params['basic_options'].append('--sc')
        elif dna_src == self.PARAM_IN_METAGENOME:
            params['basic_options'].append('--meta')
        elif dna_src == self.PARAM_IN_PLASMID:
            params['basic_options'].append('--plasmid')
        elif dna_src == self.PARAM_IN_RNA:
            params['basic_options'].append('--rna')
        elif dna_src == self.PARAM_IN_IONTORRENT:
            params['basic_options'].append('--iontorrent')

        # processing pipeline option parameters
        if params.get(self.PARAM_IN_PIPELINE_OPTION, None):
            pipe_opts = params[self.PARAM_IN_PIPELINE_OPTION]
            opts = [self.PARAM_IN_ONLY_ERROR_CORR,
                    self.PARAM_IN_ONLY_ASSEMBLER,
                    self.PARAM_IN_CONTINUE,
                    self.PARAM_IN_DISABLE_GZIP,
                    self.PARAM_IN_CAREFUL]
            if any(elem in opts for elem in pipe_opts):
                pass
            else:
                params[self.PARAM_IN_PIPELINE_OPTION] = [self.PARAM_IN_CAREFUL]
        else:
            params[self.PARAM_IN_PIPELINE_OPTION] = [self.PARAM_IN_CAREFUL]

        if '--meta' in params['basic_options']:
            # you cannot specify --careful, --mismatch-correction
            # or --cov-cutoff in metagenomic mode!
            try:
                params[self.PARAM_IN_PIPELINE_OPTION].remove(self.PARAM_IN_CAREFUL)
                params[self.PARAM_IN_PIPELINE_OPTION].remove('mismatch-correction')
                params[self.PARAM_IN_PIPELINE_OPTION].remove('cov-cutoff')
            except ValueError:
                pass

        if params.get('create_report', None) is None:
            params['create_report'] = 0

        return params

    def generate_report(self, fa_file_name, params, out_dir, wsname):
        """
        Generating and saving report
        """
        log('Generating and saving report')

        fa_file_with_path = os.path.join(out_dir, fa_file_name)
        fasta_stats = self._load_stats(fa_file_with_path)
        lengths = [fasta_stats[contig_id] for contig_id in fasta_stats]

        assembly_ref = wsname + '/' + params[self.PARAM_IN_CS_NAME]

        report_text = ''
        report_text += 'SPAdes results saved to: ' + wsname + '/' + out_dir + '\n'
        report_text += 'Assembly saved to: ' + assembly_ref + '\n'
        report_text += 'Assembled into ' + str(len(lengths)) + ' contigs.\n'
        report_text += 'Avg Length: ' + str(sum(lengths) / float(len(lengths))) + ' bp.\n'

        # compute a simple contig length distribution
        bins = 10
        counts, edges = np.histogram(lengths, bins)
        report_text += 'Contig Length Distribution (# of contigs -- min to max ' + 'basepairs):\n'
        for c in range(bins):
            report_text += ('   ' + str(counts[c]) + '\t--\t' + str(edges[c]) + ' to ' +
                            str(edges[c + 1]) + ' bp\n')
        print('Running QUAST')
        quastret = self.kbq.run_QUAST(
            {'files': [{'path': fa_file_with_path, 'label': params[self.PARAM_IN_CS_NAME]}]})

        output_files = self._generate_output_file_list(out_dir)

        print('Saving report')
        report_output = self.kbr.create_extended_report(
            {'message': report_text,
             'objects_created': [{'ref': assembly_ref, 'description': 'Assembled contigs'}],
             'direct_html_link_index': 0,
             'file_links': output_files,
             'html_links': [{'shock_id': quastret['shock_id'],
                             'name': 'report.html',
                             'label': 'QUAST report'}
                            ],
             'report_object_name': 'kb_spades_report_' + str(uuid.uuid4()),
             'workspace_name': params[self.PARAM_IN_WS]})

        return report_output['name'], report_output['ref']

    def get_hybrid_reads_info(self, input_params):
        """
        get_hybrid_reads_info--from a list of ReadsParams structures fetches the corresponding
        reads info with the ReadsParams[lib_ref]
        returns None or a tuple of nine reads data each is a list of the following structure:
        {
                'fwd_file': path_to_fastq_file,
                'orientation': (default value is "fr" (forward-reverse) for paired-end libraries
                                "rf" (reverse-forward) for mate-pair libraries), None for others
                'lib_type': ("paired-end", "mate-pairs", "hq-mate-pairs", "single", "pacbio",
                              "nanopore", "sanger", "trusted-contigs", "untrusted-contigs"),
                'type': reads_type, # 'interleaved', 'paired', or 'single'
                'seq_tech': sequencing_tech,
                'reads_ref': KBase object ref for downstream convenience,
                'reads_name': KBase object name for downstream convenience,
                'rev_file': path_to_fastq_file  # only if paired end
        }
        OR:
        {
                'fwd_file': path_to_fastq_file,
                'long_reads_type': ("pacbio-ccs", "pacbio-clr", "nanopore", "sanger",
                                    "trusted-contigs", "untrusted-contigs"),
                'type': reads_type, # 'interleaved', 'paired', or 'single'
                'seq_tech': sequencing_tech,
                'reads_ref': KBase object ref for downstream convenience,
                'reads_name': KBase object name for downstream convenience
        }
        """
        rds_params = copy.deepcopy(input_params)
        if rds_params.get(self.PARAM_IN_READS, None) is None:
            return ()  # an empty tuple

        wsname = rds_params[self.PARAM_IN_WS]

        sgl_rds_data = []  # single
        pe_rds_data = []   # paired-end
        mp_rds_data = []   # mate-pairs
        pb_ccs_data = []   # pacbio-ccs
        pb_clr_data = []   # pacbio-clr
        np_rds_data = []   # nanopore
        sgr_rds_data = []  # sanger
        tr_ctg_data = []   # trusted-contigs
        ut_ctg_data = []   # untrusted-contigs

        # a list of Illumina or IonTorrent paired-end/high-quality mate-pairs/unpaired reads
        rds_refs = []

        rds_libs = rds_params[self.PARAM_IN_READS]
        for rds_lib in rds_libs:
            if rds_lib.get('lib_ref', None):
                rds_refs.append(rds_lib['lib_ref'])
        kb_rds_data = self._get_kbreads_info(wsname, rds_refs)

        for rds_lib in rds_libs:
            for kb_d in kb_rds_data:
                if 'lib_ref' in rds_lib and rds_lib['lib_ref'] == kb_d['reads_ref']:
                    if rds_lib['lib_type'] == 'single':  # single end reads grouped params
                        kb_d['orientation'] = None
                        kb_d['lib_type'] = 'single'
                        sgl_rds_data.append(kb_d)
                    elif rds_lib['lib_type'] == 'paired-end':  # pairedEnd reads grouped params
                        kb_d['orientation'] = ('fr' if rds_lib.get('orientation', None) is None
                                               else rds_lib['orientation'])
                        kb_d['lib_type'] = 'paired-end'
                        pe_rds_data.append(kb_d)
                    elif rds_lib['lib_type'] == 'mate-pairs':
                        # mate-pairs reads grouped params
                        kb_d['orientation'] = ('rf' if rds_lib.get('orientation', None) is None
                                               else rds_lib['orientation'])
                        kb_d['lib_type'] = 'mate-pairs'
                        mp_rds_data.append(kb_d)

        # a list of PacBio (CCS or CLR), Oxford Nanopore Sanger reads
        # and/or additional contigs
        long_rds_refs = []
        if rds_params.get(self.PARAM_IN_LONG_READS, None):
            long_rds_libs = rds_params[self.PARAM_IN_LONG_READS]
            for lrds_lib in long_rds_libs:
                if lrds_lib.get('long_reads_ref', None):
                    long_rds_refs.append(lrds_lib['long_reads_ref'])
            kb_lrds_data = self._get_kbreads_info(wsname, long_rds_refs)

            for lrds_lib in long_rds_libs:
                for kb_ld in kb_lrds_data:
                    if ('long_reads_ref' in lrds_lib and
                            lrds_lib['long_reads_ref'] == kb_ld['reads_ref']):
                        if lrds_lib['long_reads_type'] == 'pacbio-ccs':
                            kb_ld['long_reads_type'] = lrds_lib['long_reads_type']
                            pb_ccs_data.append(kb_ld)
                        elif lrds_lib['long_reads_type'] == 'pacbio-clr':
                            kb_ld['long_reads_type'] = lrds_lib['long_reads_type']
                            pb_clr_data.append(kb_ld)
                        elif lrds_lib['long_reads_type'] == 'nanopore':
                            kb_ld['long_reads_type'] = lrds_lib['long_reads_type']
                            np_rds_data.append(kb_ld)
                        elif lrds_lib['long_reads_type'] == 'sanger':
                            kb_ld['long_reads_type'] = lrds_lib['long_reads_type']
                            sgr_rds_data.append(kb_ld)
                        elif lrds_lib['long_reads_type'] == 'trusted-contigs':
                            kb_ld['long_reads_type'] = lrds_lib['long_reads_type']
                            tr_ctg_data.append(kb_ld)
                        elif lrds_lib['long_reads_type'] == 'untrusted-contigs':
                            kb_ld['long_reads_type'] = lrds_lib['long_reads_type']
                            ut_ctg_data.append(kb_ld)

        return (sgl_rds_data, pe_rds_data, mp_rds_data, pb_ccs_data, pb_clr_data, np_rds_data,
                sgr_rds_data, tr_ctg_data, ut_ctg_data)

    def construct_yaml_dataset_file(self, sgl_libs=None, pe_libs=None, mp_libs=None,
                                    pb_ccs=None, pb_clr=None, np_libs=None,
                                    sgr_libs=None, tr_ctgs=None, ut_ctgs=None):
        """
        construct_yaml_dataset_file: Specifying input data with YAML data set file (advanced)
        An alternative way to specify an input data set for SPAdes is to create a YAML
        data set file.
        By using a YAML file you can provide an unlimited number of paired-end, mate-pair
        and unpaired libraries. Basically, YAML data set file is a text file, in which input
        libraries are provided as a comma-separated list in square brackets. Each library is
        provided in braces as a comma-separated list of attributes.

        The following attributes are available:

            - orientation ("fr", "rf", "ff")
            - type ("paired-end", "mate-pairs", "hq-mate-pairs", "single", "pacbio", "nanopore",
                "sanger", "trusted-contigs", "untrusted-contigs")
            - interlaced reads (comma-separated list of files with interlaced reads)
            - left reads (comma-separated list of files with left reads)
            - right reads (comma-separated list of files with right reads)
            - single reads (comma-separated list of files with single reads or unpaired reads from
                paired library)
            - merged reads (comma-separated list of files with merged reads)

        To properly specify a library you should provide its type and at least one file with reads.
        For ONT, PacBio, Sanger and contig libraries you can provide only single reads. Orientation
        is an optional attribute. Its default value is "fr" (forward-reverse) for paired-end
        libraries and "rf" (reverse-forward) for mate-pair libraries.

        The value for each attribute is given after a colon. Comma-separated lists of files should
        be given in square brackets.
        For each file you should provide its full path in double quotes. Make sure that files with
        right reads are given in the same order as corresponding files with left reads.

        For example, if you have one paired-end library splitted into two pairs of files:
            lib_pe1_left_1.fastq
            lib_pe1_right_1.fastq
            lib_pe1_left_2.fastq
            lib_pe1_right_2.fastq

        one mate-pair library:
            lib_mp1_left.fastq
            lib_mp1_right.fastq

        and PacBio CCS and CLR reads:
            pacbio_ccs.fastq
            pacbio_clr.fastq

        YAML file should look like this:
        ------------------------------------------------
        [
            {
                orientation: "fr",
                type: "paired-end",
                right reads: [
                "/FULL_PATH_TO_DATASET/lib_pe1_right_1.fastq",
                "/FULL_PATH_TO_DATASET/lib_pe1_right_2.fastq"
                ],
                left reads: [
                "/FULL_PATH_TO_DATASET/lib_pe1_left_1.fastq",
                "/FULL_PATH_TO_DATASET/lib_pe1_left_2.fastq"
                ]
            },
            {
                orientation: "rf",
                type: "mate-pairs",
                right reads: [
                "/FULL_PATH_TO_DATASET/lib_mp1_right.fastq"
                ],
                left reads: [
                "/FULL_PATH_TO_DATASET/lib_mp1_left.fastq"
                ]
            },
            {
                type: "single",
                single reads: [
                "/FULL_PATH_TO_DATASET/pacbio_ccs.fastq"
                ]
            },
            {
                type: "pacbio",
                single reads: [
                "/FULL_PATH_TO_DATASET/pacbio_clr.fastq"
                ]
            }
        ]
        ------------------------------------------------

        Once you have created a YAML file save it with .yaml extension (e.g. as my_data_set.yaml)
        and run SPAdes using the --dataset option:
        e.g., <SPAdes_bin_dir>/spades.py --dataset <your YAML file> -o spades_output

        """
        # STEP 1: get the working folder housing the .yaml file and the SPAdes results
        if not os.path.exists(self.proj_dir):
            os.makedirs(self.proj_dir)
        yaml_file_path = os.path.join(self.proj_dir, 'input_data_set.yaml')

        # STEP 2: construct and save the 'input_data_set.yaml' file
        # generate the object array
        input_data_set = []

        if pe_libs:
            pair_libs = self._parse_pair_reads('paired-end', pe_libs)
            if pair_libs:
                input_data_set.append(pair_libs)

        if mp_libs:
            pair_libs = self._parse_pair_reads('mate-pairs', mp_libs)
            if pair_libs:
                input_data_set.append(pair_libs)

        # for reads_type = 'single'
        if sgl_libs:
            single_libs = self._parse_single_reads("single", sgl_libs)
            if single_libs:
                input_data_set.append(single_libs)

        # for long_reads_type = 'pacbio-ccs', treated as type of 'single'
        if pb_ccs:
            single_libs = self._parse_single_reads("single", pb_ccs)
            if single_libs:
                input_data_set.append(single_libs)

        # for long_reads_type = 'pacbio-clr'
        if pb_clr:
            single_libs = self._parse_single_reads("pacbio", pb_clr)
            if single_libs:
                input_data_set.append(single_libs)

        # for long_reads_type = 'nanopore'
        if np_libs:
            single_libs = self._parse_single_reads("nanopore", np_libs)
            if single_libs:
                input_data_set.append(single_libs)

        # for long_reads_type = 'sanger'
        if sgr_libs:
            single_libs = self._parse_single_reads("sanger", sgr_libs)
            if single_libs:
                input_data_set.append(single_libs)

        # for long_reads_type = 'trusted-contigs'
        if tr_ctgs:
            single_libs = self._parse_single_reads("trusted-contigs", tr_ctgs)
            if single_libs:
                input_data_set.append(single_libs)

        # for long_reads_type = 'untrusted-contigs'
        if ut_ctgs:
            single_libs = self._parse_single_reads("untrusted-contigs", ut_ctgs)
            if single_libs:
                input_data_set.append(single_libs)

        if input_data_set == []:
            print('Empty input data set!!')
            return ''

        pprint(input_data_set)
        try:
            with open(yaml_file_path, 'w') as yaml_file:
                json.dump(input_data_set, yaml_file)
        except IOError as ioerr:
            log('Creation of the {} file raised error:\n'.format(yaml_file_path))
            pprint(ioerr)
            return ''
        else:
            return yaml_file_path

    def run_assemble(self, yaml_file, kmer_sizes, dna_source=None,
                     basic_opts=None, pipeline_opts=['careful']):
        """
        run_assemble: run the SPAdes assemble with given input parameters/options
        """
        exit_code = 1
        if not os.path.isfile(yaml_file):
            log("The input data set yaml file DOES NOT exist at {}\n".format(yaml_file))
            return exit_code

        log("The input data set yaml file exists at {}\n".format(yaml_file))
        yf_dir, yf_nm = os.path.split(yaml_file)

        mem = (psutil.virtual_memory().available / self.GB - self.MEMORY_OFFSET_GB)
        if mem < self.MIN_MEMORY_GB:
            raise ValueError(
                'Only ' + str(psutil.virtual_memory().available) +
                ' bytes of memory are available. The SPAdes wrapper will' +
                ' not run without at least ' +
                str(self.MIN_MEMORY_GB + self.MEMORY_OFFSET_GB) +
                ' gigabytes available')

        if dna_source and dna_source == self.PARAM_IN_METAGENOME:
            max_mem = self.MAX_MEMORY_GB_META_SPADES
            max_threads = self.MAX_THREADS_META
        else:
            max_mem = self.MAX_MEMORY_GB_SPADES
            max_threads = self.MAX_THREADS

        threads = min(max_threads, psutil.cpu_count() * self.THREADS_PER_CORE)

        if mem > max_mem:
            mem = max_mem

        tmpdir = os.path.join(self.proj_dir, 'spades_tmp_dir')
        if not os.path.exists(tmpdir):
            os.makedirs(tmpdir)

        a_cmd = [os.path.join(self.SPADES_BIN, 'spades.py')]
        a_cmd += ['--threads', str(threads), '--memory', str(mem)]
        a_cmd += ['--tmp-dir', tmpdir]
        a_cmd += ['--dataset', yaml_file]

        if kmer_sizes is not None:
            a_cmd += ['-k ' + kmer_sizes]

        if basic_opts is None:
            basic_opts = ['-o', self.ASSEMBLE_RESULTS_DIR]
        if isinstance(basic_opts, list):
            a_cmd += basic_opts

        if pipeline_opts and isinstance(pipeline_opts, list):
            for p_opt in pipeline_opts:
                if p_opt == self.PARAM_IN_CAREFUL:
                    a_cmd += ['--careful']
                if p_opt == self.PARAM_IN_ONLY_ERROR_CORR:
                    a_cmd += ['--only-error-correction']
                if p_opt == self.PARAM_IN_ONLY_ASSEMBLER:
                    a_cmd += ['--only-assembler']
                if p_opt == self.PARAM_IN_CONTINUE:
                    a_cmd += ['--continue']
                if p_opt == self.PARAM_IN_DISABLE_GZIP:
                    a_cmd += ['--disable-gzip-output']

        # Last check of command options before the call
        if '--meta' in a_cmd:
            # you cannot specify --careful, --mismatch-correction
            # or --cov-cutoff in metagenomic mode!
            try:
                a_cmd.remove(self.PARAM_IN_CAREFUL)
                a_cmd.remove('mismatch-correction')
                a_cmd.remove('cov-cutoff')
            except ValueError:
                pass

        log("**************The HybridSPAdes assembling command is:\n{}".format(' '.join(a_cmd)))
        assemble_out_dir = os.path.join(self.proj_dir, self.ASSEMBLE_RESULTS_DIR)
        if not os.path.exists(assemble_out_dir):
            os.makedirs(assemble_out_dir)

        p = subprocess.Popen(a_cmd, cwd=yf_dir, shell=False)
        exit_code = p.wait()
        log('Return code: ' + str(exit_code))

        if p.returncode != 0:
            raise ValueError('Error running spades.py, return code: ' + str(p.returncode) + '\n')
        else:
            exit_code = p.returncode
        return exit_code

    def save_assembly(self, fa_file_path, wsname, a_name, min_ctg_length=0):
        """
        save_assembly: save the assembly to KBase workspace
        """
        if os.path.isfile(fa_file_path):
            log('Uploading FASTA file to Assembly...')
            if min_ctg_length > 0:
                self.au.save_assembly_from_fasta(
                            {'file': {'path': fa_file_path},
                             'workspace_name': wsname,
                             'assembly_name': a_name,
                             'min_contig_length': min_ctg_length})
            else:
                self.au.save_assembly_from_fasta(
                            {'file': {'path': fa_file_path},
                             'workspace_name': wsname,
                             'assembly_name': a_name})
        else:
            log("The resulting sequence file {} is not found.".format(fa_file_path))

    # end of public methods


text:
import logging
import joblib
from pathlib import Path

import numpy as np
from sklearn.metrics import classification_report, confusion_matrix
from sklearn.model_selection import train_test_split

from sklearn.decomposition import KernelPCA
from sklearn.datasets import load_iris
from sklearn.pipeline import Pipeline
from sklearn.ensemble import RandomForestClassifier
from sklearn import svm

logger = logging.getLogger(__name__)

def run_sklearn_pipeline(config):
    """Demo of a sklearn pipeline.

    read:
    - https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html
    - https://gist.github.com/amberjrivera/8c5c145516f5a2e894681e16a8095b5c
    """
    # load iris data
    data = load_iris()
    X_train, X_test, y_train, y_test = train_test_split(
        data.data, data.target, test_size=0.25, random_state=999
    )

    # import model
    model_name = config["model"]["name"]
    model_params = config["model"]["params"]
    # @note: compare the two cases
    # if you want to run reproducible experiments at scale,
    # you'd better not hard-code parameters like below in the source code.
    try:
        if model_name == "random_forest":
            # case - A
            clf = RandomForestClassifier(max_depth=2, random_state=999)
        elif model_name == "support_vector_machine":
            # case - B
            kernel = model_params["svm_kernel"]
            gamma = model_params["svm_gamma"]
            C = model_params["svm_C"]
            clf = svm.SVC(kernel=kernel, gamma=gamma, C=C)
        else:
            raise ValueError(f"Unsupported model_name was given {model_name}.")
    except Exception as e:
        logger.error(f"{e}")

    # create a pipeline
    pipeline = Pipeline([("kernel_pca", KernelPCA()), ("model", clf)])

    # train
    pipeline.fit(X_train, y_train)

    # evaluate
    y_pred = pipeline.predict(X_test)

    report = classification_report(y_test, y_pred)
    logger.info(f"Accuracy report:\n{report}")

    cm = confusion_matrix(y_test, y_pred)
    logger.info(f"Confusion Matrix:\n{cm}")

    save_path = Path(config["experiment"]["output_dirname"]) / f"{model_name}.joblib"
    save_sklearn_model(pipeline, save_path)
    logger.info(f"Saved sklearn pipeline model at {save_path}")

    # check that the saved model behaves the same as before
    checkpoint = load_sklearn_model(save_path)
    assert np.all(checkpoint.predict(X_test) == pipeline.predict(X_test))

def save_sklearn_model(model, save_path):
    try:
        joblib.dump(model, save_path)
    except Exception as e:
        logger.error(f"{e}")

def load_sklearn_model(load_path):
    try:
        model = joblib.load(load_path)
    except Exception as e:
        logger.error(f"{e}")
    return model

ext: py

reason: import numpy


text:
import argparse
import csv
import math
import os
import sys
import tensorflow as tf
from scipy.stats import chi2
import numpy as np
import tensorflow.contrib.slim as slim
import tensorflow.contrib.metrics as metrics

ARGS = None

def get_file_dims(input_file, header=True):
    num_columns = 0
    num_genes = 0
    with open(input_file) as fp:
        line_no = 0
        for line in fp:
            line_no += 1
            if (line_no == 1 and header == True):
                next
            else:
                new_num_columns = len(line.split("\t"))
                if (num_columns != 0 and num_columns != new_num_columns):
                    print("Error the number of columns isn't consistent")
                    print(line)
                    print("Num genes " + str(num_genes))
                    print("New columns " + str(new_num_columns))
                    print("Old columns " + str(num_columns))
                    exit(1)
                num_columns = new_num_columns
                num_genes += 1
    return num_genes, num_columns

def load_data(input_file, num_features):
    evidence = np.loadtxt(
        input_file,
        dtype="float64",
        skiprows=1,
        delimiter="\t",
        usecols=tuple([i for i in range(2, 2 + num_features)]))
    data = np.loadtxt(
        input_file, dtype="float64", skiprows=1, delimiter="\t", usecols=1)
    genes = np.loadtxt(
        input_file, dtype="str", skiprows=1, delimiter="\t", usecols=0)
    return genes, data, evidence

def mstep(new_a, loglik, scope, reuse=True):
    with tf.variable_scope(scope, reuse=reuse):
        w = tf.get_variable("w", dtype=tf.float64)
        a = tf.get_variable("a", dtype=tf.float64)
    a_update = tf.assign(
        a,
        tf.maximum(
            tf.convert_to_tensor(0.000001, tf.float64), tf.minimum(new_a, 1)))
    w_update = tf.train.AdamOptimizer(0.01).minimize(-1 * loglik)
    return a_update, w_update, w

def estep(x, y, n, m, scope, reuse=False):
    # ESTEP: Compute Posteriors
    with tf.variable_scope(scope, reuse=reuse):
        w = tf.get_variable(
            "w", shape=[m, 1], dtype=tf.float64)
        a = tf.get_variable(
            "a",
            dtype=tf.float64,
            trainable=False,
            initializer=tf.random_uniform([],
                                          minval=0,
                                          maxval=1,
                                          dtype=tf.float64))
        #tf.add_to_collection("NON_GRADIENT", a)

    # Compute Pr(zg=1| xg1..xgm, w1...wm)
    if m == 1:
        prob_zg1 = tf.sigmoid(
            tf.squeeze(tf.matmul(tf.ones([1, 1], dtype=tf.float64), w)))
    else:
        prob_zg1 = tf.sigmoid(tf.squeeze(tf.matmul(x, w)))
    # Compute p(yg, zg=1 | xg1..xgm, w1...wm)
    prob_pg1 = a * tf.pow(y, a - 1) * prob_zg1
    # Compute p(yg, zg=0 | xg1..xgm, w1...wm)
    prob_pg0 = 1 - prob_zg1
    # Compute p(yg| xg1..xgm, w1...wm)
    prob_pg = prob_pg1 + prob_pg0
    # Compute p(zg=1|yg, xg1..xgm, w1...wm)
    post_zg1 = tf.div(prob_pg1, prob_pg)
    # Compute total loglik
    loglik = tf.reduce_sum(tf.log(prob_pg))
    # update a
    new_a = tf.div(
        tf.reduce_sum(-1 * post_zg1), tf.reduce_sum(post_zg1 * tf.log(y)))
    return loglik, new_a, post_zg1

def inference(x, y, n, m, scope, reuse=False):
    loglik, new_a, post_zg1 = estep(x, y, n, m, scope, reuse)
    a_update, w_update, new_w = mstep(new_a, loglik, scope)
    with tf.control_dependencies([a_update, w_update]):
        new_loglik, blah2, new_post_zg1 = estep(x, y, n, m, scope, True)
    return new_loglik, a_update, w_update, new_w

def next_step(sess, ops, scope, i, feed_dict=None):
    output = sess.run(ops, feed_dict=feed_dict)
    model = {
        'name': scope,
        'i': i,
        'loglik': output[0],
        'a': output[1],
        'w': output[3].tolist()
    }
    return model

def get_model(sess, loglik_op, scope, i, feed_dict=None):
    loglik = sess.run(
        loglik_op, feed_dict=feed_dict)
    a = sess.run(get_var(scope, 'a'))
    w = sess.run(get_var(scope, 'w'))
    model = {'name': scope, 'i': i, 'loglik': loglik, 'a': a, 'w': w.tolist()}
    return model

def print_models(models):
    for model in models:
        print_model(model)

def print_model(model):
    print(str_model(model) + "\n")

def str_model(model):
    w_str = str(np.squeeze(model["w"])).replace('[', '').replace(']',
                                                                 '').replace(
                                                                     ',', '\t')
    return str(model["name"]) + "\t" + str(model["i"]) + "\t" + str(
        model['loglik']) + "\t" + w_str + "\t" + str(model["a"])

def get_loglik(x, y, n, m, scope, reuse=True):
    loglik, new_a, post_zg1 = estep(x, y, n, m, scope, reuse)
    return loglik

def get_var(scope, var):
    with tf.variable_scope(scope, reuse=True):
        return tf.get_variable(var, dtype=tf.float64)

def run_graph(input_file,
              evidence,
              data,
              num_genes,
              num_features,
              min_steps=5000,
              max_steps=10000,
              delta=0.001,
              last_n=1000):
    global_step = tf.Variable(0, trainable=False, name='global_step')
    x = tf.convert_to_tensor(evidence)
    y = tf.convert_to_tensor(data)
    h1 = inference(x, y, num_genes, num_features, input_file)

    with tf.Session() as sess:
        sess.run(tf.global_variables_initializer())
        h1_prev = get_model(
            sess, h1[0], input_file, 'INITIAL_0', feed_dict=None)
        h1_init = h1_prev
        i = 1
        delta_loglik = delta + 1
        h1_last_n_loglik = h1_prev["loglik"]
        while ((i <= min_steps) or (delta_loglik > delta and i <= max_steps)):
            h1_next = next_step(
                sess, h1, input_file, 'FINAL_' + str(i), feed_dict=None)
            if (i % last_n == 0):
                delta_loglik = h1_next["loglik"] - h1_last_n_loglik
                h1_last_n_loglik = h1_next["loglik"]
            i += 1

        h1_init_str = str_model(h1_init)
        h1_next_str = str_model(h1_next)
        print(h1_init_str + "\t" + h1_next_str)

def main(_):
    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
    seed = ARGS.seed
    input_file = ARGS.input_file
    tf.set_random_seed(seed)
    num_genes, num_columns = get_file_dims(input_file)
    num_features = num_columns - 2
    genes, data, evidence = load_data(input_file, num_features)
    run_graph(input_file, evidence, data, num_genes, num_features)

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--input_file",
        type=str,
        default="input.tsv",
        help='File to run pGENMi on')
    parser.add_argument(
        '--seed',
        type=int,
        default='10902',
        help='Set seed for reproducibility')
    ARGS, unparsed = parser.parse_known_args()
    tf.app.run(main=main, argv=[sys.argv[0]] + unparsed)

max_stars_repo_name: KnowEnG/pgenmi

max_issues_repo_name: KnowEnG/pgenmi

max_forks_repo_name: KnowEnG/pgenmi

reason: import numpy,from scipy


text:
# ============================================================================
# 付録 A 石油温水暖房機及び石油給湯温水暖房機
# ============================================================================

import numpy as np
import pyhees.section4_7_h as appendix_H

# ============================================================================
# A.2 エネルギー消費量
# ============================================================================

# ============================================================================
# A.2.1 灯油消費量
# ============================================================================

def calc_E_K_hs(Q_out_H_hs, e_rtd, hs_type, Theta_SW_hs, Theta_RW_hs, region, A_A, A_MR, A_OR, mode_MR, mode_OR,
                has_MR_hwh, has_OR_hwh):
    """温水暖房用熱源機の灯油消費量 (1)

    Args:
      Q_out_H_hs(ndarray): 熱源機暖房出力 (MJ/h)
      e_rtd(float): 当該給湯機の効率
      hs_type(str): 温水暖房用熱源機の種類
      Theta_SW_hs(ndarray): 温水暖房用熱源機の往き温水温度
      Theta_RW_hs(ndarray): 1時間平均の熱源機の戻り温水温度 (℃)
      region(int): 省エネルギー地域区分
      A_A(float): 床面積の合計 (m2)
      A_MR(float): 主たる居室の床面積 (m2)
      A_OR(float): その他の居室の床面積 (m2)
      mode_MR(str): 主たる居室の運転モード 'い', 'ろ', 'は'
      mode_OR(str): その他の居室の運転モード 'い', 'ろ', 'は'
      has_MR_hwh(bool): 温水暖房の放熱器を主たる居室に設置する場合はTrue
      has_OR_hwh(bool): 温水暖房の放熱器をその他の居室に設置する場合はTrue

    Returns:
      ndarray: 温水暖房用熱源機の灯油消費量

    """
    # 定格試験時の温水暖房用熱源機の筐体放熱損失 (4)
    Q_body_rtd = get_Q_body_rtd(hs_type)

    # 温水暖房熱源機の定格能力
    q_rtd_hs = calc_q_rtd_hs(region, A_A, A_MR, A_OR, mode_MR, mode_OR, has_MR_hwh, has_OR_hwh)

    # 1時間平均の温水暖房用熱源機の熱交換効率 (3)
    e_ex = calc_e_ex(e_rtd, hs_type, Theta_RW_hs, q_rtd_hs, Q_body_rtd)

    # 温水暖房用熱源機の筐体放熱損失 (2)
    Q_body = get_Q_body(hs_type, Theta_SW_hs)

    # 温水暖房用熱源機の灯油消費量 (1)
    E_K_hs = (Q_out_H_hs + Q_body) / e_ex
    E_K_hs[Q_out_H_hs == 0] = 0

    return E_K_hs

def get_Q_body(hs_type, Theta_SW_hs):
    """温水暖房用熱源機の筐体放熱損失 (2)

    Args:
      hs_type(str): 温水暖房用熱源機の種類
      Theta_SW_hs(ndarray): 温水暖房用熱源機の往き温水温度

    Returns:
      ndarray: 温水暖房用熱源機の筐体放熱損失

    """
    if hs_type in ['石油従来型暖房機', '石油従来型温水暖房機', '石油従来型給湯温水暖房機', '不明']:
        # (2a)
        return [234 * 3600 * 10 ** (-6)] * 24 * 365
    elif hs_type in ['石油潜熱回収型暖房機', '石油潜熱回収型温水暖房機', '石油潜熱回収型給湯温水暖房機']:
        # (2b)
        return (5.3928 * Theta_SW_hs - 71.903) * 3600 * 10 ** (-6)
    else:
        raise ValueError(hs_type)

def calc_e_ex(e_rtd, hs_type, Theta_RW_hs, q_rtd_hs, Q_body_rtd):
    """1時間平均の温水暖房用熱源機の熱交換効率 (3)

    Args:
      e_rtd(float): 当該給湯機の効率
      hs_type(str): 温水暖房用熱源機の種類
      Theta_RW_hs(ndarray): 1時間平均の熱源機の戻り温水温度 (℃)
      q_rtd_hs(float): 温水暖房用熱源機の定格能力 (W)
      Q_body_rtd(ndarray): 1時間当たりの定格試験時の温水暖房用熱源機の筐体放熱損失（MJ/h）

    Returns:
      ndarray: 1時間平均の温水暖房用熱源機の熱交換効率 (3)

    """
    # 1時間平均の定格効率を補正する係数
    f_rtd = get_f_rtd(hs_type, Theta_RW_hs)

    return e_rtd * f_rtd * (q_rtd_hs * 3600 * 10 ** (-6) + Q_body_rtd) / (q_rtd_hs * 3600 * 10 ** (-6))

def get_Q_body_rtd(hs_type):
    """定格試験時の温水暖房用熱源機の筐体放熱損失 (4)

    Args:
      hs_type(str): 温水暖房用熱源機の種類

    Returns:
      ndarray: 定格試験時の温水暖房用熱源機の筐体放熱損失 (4)

    """
    if hs_type in ['石油従来型暖房機', '石油従来型温水暖房機', '石油従来型給湯温水暖房機']:
        # (4a)
        return 234 * 3600 * 10 ** (-6) * np.ones(24*365)
    elif hs_type in ['石油潜熱回収型暖房機', '石油潜熱回収型温水暖房機', '石油潜熱回収型給湯温水暖房機']:

        # 温水暖房用熱源機の定格試験時の往き温水温度
        Theta_SW_hs_rtd = get_Theta_SW_hs_rtd()

        # (4b)
        return (5.3928 * Theta_SW_hs_rtd - 71.903) * 3600 * 10 ** (-6)
    else:
        raise ValueError(hs_type)

def get_f_rtd(hs_type, Theta_RW_hs):
    """1時間平均の定格効率を補正する係数 (5)

    Args:
      hs_type(str): 温水暖房用熱源機の種類
      Theta_RW_hs(ndarray): 1時間平均の熱源機の戻り温水温度 (℃)

    Returns:
      ndarray: 1時間平均の定格効率を補正する係数 (5)

    """
    f_rtd = np.zeros(24 * 365)

    if hs_type in ['石油従来型暖房機', '石油従来型温水暖房機', '石油従来型給湯温水暖房機']:
        # (5a)
        f_rtd[:] = 0.946
    elif hs_type in ['石油潜熱回収型暖房機', '石油潜熱回収型温水暖房機', '石油潜熱回収型給湯温水暖房機']:
        # (5b)

        Theta_RW_hs_rtd = get_Theta_RW_hs_rtd()

        # 条件1
        f1 = np.logical_and(Theta_RW_hs >= 46.5, Theta_RW_hs_rtd >= 46.5)
        f_rtd[f1] = 0.9768

        # 条件2
        f2 = np.logical_and(Theta_RW_hs < 46.5, Theta_RW_hs_rtd >= 46.5)
        f_rtd[f2] = (-0.0023 * Theta_RW_hs[f2] + 1.014) / 0.907 * 0.9768

        # 条件3
        f3 = np.logical_and(Theta_RW_hs >= 46.5, Theta_RW_hs_rtd < 46.5)
        f_rtd[f3] = 0.907 / (-0.0023 * Theta_RW_hs[f3] + 1.014) * 0.9768

        # 条件4
        f4 = np.logical_and(Theta_RW_hs < 46.5, Theta_RW_hs_rtd < 46.5)
        f_rtd[f4] = (-0.0023 * Theta_RW_hs[f4] + 1.014) / (-0.0023 * Theta_RW_hs_rtd + 1.014) * 0.9768

    else:
        raise ValueError(hs_type)

    return f_rtd

def get_e_rtd_default(hs_type):
    """定格効率(規定値)

    Args:
      hs_type(str): 温水暖房用熱源機の種類

    Returns:
      float: 定格効率(規定値)

    """
    if hs_type in ['石油従来型暖房機', '石油従来型温水暖房機', '石油従来型給湯温水暖房機']:
        return 0.82
    elif hs_type in ['石油潜熱回収型暖房機', '石油潜熱回収型温水暖房機', '石油潜熱回収型給湯温水暖房機']:
        return 0.91
    else:
        raise ValueError(hs_type)

def get_Theta_SW_hs_rtd():
    """温水暖房用熱源機の定格試験時の往き温水温度

    Args:

    Returns:
      float: 温水暖房用熱源機の定格試験時の往き温水温度

    """
    return 70

def get_Theta_RW_hs_rtd():
    """温水暖房熱源機の定格試験時の戻り温水温度

    Args:

    Returns:
      float: 温水暖房熱源機の定格試験時の戻り温水温度

    """
    return 50

def calc_q_rtd_hs(region, A_A, A_MR, A_OR, mode_MR, mode_OR, has_MR_hwh, has_OR_hwh):
    """温水暖房熱源機の定格能力

    Args:
      region(int): 省エネルギー地域区分
      A_A(float): 床面積の合計 (m2)
      A_MR(float): 主たる居室の床面積 (m2)
      A_OR(float): その他の居室の床面積 (m2)
      mode_MR(str): 主たる居室の運転モード 'い', 'ろ', 'は'
      mode_OR(str): その他の居室の運転モード 'い', 'ろ', 'は'
      has_MR_hwh(bool): 温水暖房の放熱器を主たる居室に設置する場合はTrue
      has_OR_hwh(bool): 温水暖房の放熱器をその他の居室に設置する場合はTrue

    Returns:
      float: 温水暖房熱源機の定格能力

    """
    # 付録Hに定める温水暖房用熱源機の最大能力 q_max_hs に等しい
    return appendix_H.calc_q_max_hs(region, A_A, A_MR, A_OR, mode_MR, mode_OR, has_MR_hwh, has_OR_hwh)

# ============================================================================
# A.2.2 消費電力量
# ============================================================================

def calc_E_E_hs(hs_type, r_WS_hs, E_K_hs):
    """温水暖房用熱源機の消費電力量 (6)

    Args:
      hs_type(str): 温水暖房用熱源機の種類
      r_WS_hs(ndarray): 1時間平均の温水暖房用熱源機の温水供給運転率
      E_K_hs(ndarray): 1時間当たりの温水暖房用熱源機の灯油消費量（MJ/h）

    Returns:
      ndarray: 温水暖房用熱源機の消費電力量

    """
    # 温水暖房用熱源機の送水ポンプの消費電力量
    E_E_hs_pump = calc_E_E_hs_pmp(hs_type, r_WS_hs)

    # 温水暖房用の熱源機の排気ファンの消費電力量
    E_E_hs_fan = calc_E_E_hs_fan(E_K_hs)

    return E_E_hs_pump + E_E_hs_fan

def calc_E_E_hs_pmp(hs_type, r_WS_hs):
    """温水暖房用熱源機の送水ポンプの消費電力量 (7)

    Args:
      hs_type(str): 温水暖房用熱源機の種類
      r_WS_hs(ndarray): 1時間平均の温水暖房用熱源機の温水供給運転率

    Returns:
      ndarray: 温水暖房用熱源機の送水ポンプの消費電力量

    """
    # 温水暖房用熱源機の送水ポンプの消費電力
    P_hs_pmp = get_P_hs_pmp(hs_type)

    return P_hs_pmp * r_WS_hs * 10 ** (-3)

def get_P_hs_pmp(hs_type):
    """温水暖房用熱源機の送水ポンプの消費電力

    Args:
      hs_type(str): 温水暖房用熱源機の種類

    Returns:
      float: 温水暖房用熱源機の送水ポンプの消費電力

    """
    if hs_type in ['石油従来型暖房機', '石油従来型温水暖房機', '石油従来型給湯温水暖房機']:
        return 90
    elif hs_type in ['石油潜熱回収型暖房機', '石油潜熱回収型温水暖房機', '石油潜熱回収型給湯温水暖房機']:
        return 70
    else:
        raise ValueError(hs_type)

def calc_E_E_hs_fan(E_K_hs):
    """温水暖房用の熱源機の排気ファンの消費電力量 (8)

    Args:
      E_K_hs(ndarray): 1時間当たりの温水暖房用熱源機の灯油消費量（MJ/h）

    Returns:
      ndarray: 温水暖房用の熱源機の排気ファンの消費電力量 (8)

    """
    # 排気ファンの効率
    gamma = get_gamma()

    return E_K_hs * gamma * 10 ** 3 / 3600

def get_gamma():
    """排気ファンの効率

    Args:

    Returns:
      float: 排気ファンの効率

    """
    return 0.003

# ============================================================================
# A.2.3 ガス消費量
# ============================================================================

def get_E_G_hs():
    """ガス消費量

    Args:

    Returns:
      ndarray: ガス消費量

    """
    return np.zeros(24 * 365)

# ============================================================================
# A.2.4 その他の燃料による一次エネルギー消費量
# ============================================================================

def get_E_M_hs():
    """その他の燃料による一次エネルギー消費量

    Args:

    Returns:
      ndarray: その他の燃料による一次エネルギー消費量

    """
    return np.zeros(24 * 365)

# ============================================================================
# A.3 最大暖房出力
# ============================================================================

def get_Q_max_H_hs(q_rtd_hs):
    """最大暖房出力 (9)

    Args:
      q_rtd_hs(float): 温水暖房用熱源機の定格能力 (W)

    Returns:
      ndarray: 最大暖房出力

    """
    return np.ones(24*365) * q_rtd_hs * 3600 / (10 ** 6)

max_stars_repo_name: jjj-design/pyhees

max_issues_repo_name: jjj-design/pyhees

max_forks_repo_name: jjj-design/pyhees


text:
import numpy as np

from rllab.envs.base import Env

from sandbox.rocky.tf import spaces

class MockPolicy(object):

    def __init__(self, action_space):
        self.action_space = action_space
        self.recurrent = False
        self.counter = None

    def reset(self, dones=None):
        if self.counter is None:
            self.counter = np.zeros(len(dones))
        for i, done in enumerate(dones):
            if done:
                self.counter[i] = 0

    def get_actions(self, observations):
        n_samples = len(observations)
        actions = [self.action_space.sample() for i in range(n_samples)]
        agent_infos = dict(
            prob=[[1,0] for _ in range(n_samples)],
            latent=[[self.counter[i]] for i in range(n_samples)],
            latent_info=dict(
                prob=[[0,1] for _ in range(n_samples)],
                latent=[[self.counter[i]] for i in range(n_samples)],
                update=[[self.counter[i] % 2 == 0] for i in range(n_samples)],
                latent_info=dict(
                    latent=[[self.counter[i] * 2]  for i in range(n_samples)]
                )
            )
        )
        self.counter += 1
        return actions, agent_infos

    @property
    def distribution(self):
        return MockDist()

class MockAlgo(object):
    def __init__(
            self,
            policy=None,
            env=None,
            baseline=None,
            discount=1.,
            batch_size=10,
            max_path_length=1000,
            gae_lambda=1.
        ):
        self.policy = policy
        self.env = env
        self.baseline = baseline
        self.discount = discount
        self.batch_size = batch_size
        self.max_path_length = max_path_length
        self.gae_lambda = gae_lambda
        self.center_adv = False
        self.positive_adv = False

class MockDist(object):

    def entropy(self, *args):
        return 0

class MockEnv(Env):
    def __init__(self, observation_space, action_space):
        self._observation_space = observation_space
        self._action_space = action_space
    def reset(self):
        return self.observation_space.sample()
    def step(self, action):
        assert self.action_space.contains(action)
        next_obs = self.observation_space.sample()
        return next_obs, 0, False, {}
    @property
    def observation_space(self):
        return self._observation_space
    @property
    def action_space(self):
        return self._action_space

class MockEnvSpec(object):
    def __init__(
            self,
            observation_space=spaces.Box(np.array([0]),np.array([1])),
            action_space=spaces.Discrete(2),
            num_envs=1):
        self.observation_space = observation_space
        self.action_space = action_space
        self.num_envs = num_envs

max_stars_repo_name: Kailiangdong/hgail

max_issues_repo_name: Kailiangdong/hgail

max_forks_repo_name: Kailiangdong/hgail


text:
import sys
import copy
import pprint
pp = pprint.PrettyPrinter(width=120)
from typing import Any, Dict, Optional
import inspect

import numpy as np
import qcelemental as qcel
from qcelemental.models import ResultInput

import qcengine as qcng
from qcengine.programs.util import PreservingDict
from qcengine.programs.psi4 import Psi4Harness
#from qcengine.programs.cfour.keywords import format_keywords, format_keyword

from .. import __version__
from .. import qcvars
#from ..driver.driver_helpers import print_variables
from ..exceptions import *
from ..util import provenance_stamp
from ..molecule import Molecule
from ..pdict import PreservingDict
#from .worker import psi4_subprocess
from .botanist import muster_inherited_options

#def run_psi4_old(name, molecule, options, **kwargs):
#
#    if kwargs['ptype'] not in ['energy', 'properties', 'gradient', 'hessian']:
#        raise ValidationError("""run_psi4: ptype not regonized: {}""".format(ptype))
#
#    prov = {}
#    prov['creator'] = 'QCDB'
#    prov['version'] = __version__
#    prov['routine'] = sys._getframe().f_code.co_name
#
#    jobrec = {}
#    jobrec['error'] = ''
#    jobrec['success'] = None
#    jobrec['return_output'] = True
#    jobrec['provenance'] = [prov]
#    jobrec['molecule'] = molecule.to_dict(np_out=False)
#    jobrec['method'] = name
#    jobrec['driver'] = kwargs['ptype']
#    jobrec['kwargs'] = kwargs
#    jobrec['options'] = copy.deepcopy(options)
#    jobrec['hooks'] = kwargs.get('hooks', {})
#
#    try:
#        jobrec['molecule']
#        jobrec['method']
#    except KeyError as err:
#        jobrec['error'] += repr(err) + 'Required fields missing from ({})'.format(jobrec.keys())
#        return jobrec
#
#
#    psi4rec = {}
#    psi4rec['json'] = {}
#
#    opts = jobrec['options']
#    # NOTE TODO very limited OPTIONS HANDSHAKE
#    muster_inherited_options(opts)
#
#    psi4rec['json']['schema_name'] = 'qcschema_input'
#    psi4rec['json']['schema_version'] = 1
#    omem = opts.scroll['QCDB'].pop('MEMORY')
#    psi4rec['json']['memory'] = omem.value
#    psi4rec['json']['molecule'] = qcel.molparse.to_schema(jobrec['molecule'], dtype=2)
#    psi4rec['json']['driver'] = jobrec['driver']
#    mtd = jobrec['method']
#    psi4rec['json']['model'] = {
#        'method': mtd[3:] if mtd.startswith('p4-') else mtd,
#        'basis': '(auto)',
#    }
#    psi4rec['json']['extras'] = {'wfn_qcvars_only': True}
#    psi4rec['json']['kwargs'] = jobrec['kwargs']
#    psi4rec['json']['return_output'] = True
#
#    popts = {}
#    for k, v in opts.scroll['QCDB'].items():
#        if v.disputed():
#            popts[k] = v.value
#
#    for k, v in opts.scroll['PSI4'].items():
#        if v.disputed():
#            popts[k] = v.value
#    psi4rec['json']['keywords'] = popts
#    if 'BASIS' in psi4rec['json']['keywords']:
#        psi4rec['json']['model']['basis'] =  psi4rec['json']['keywords']['BASIS']
#
#    psi4rec['command'] = ['psi4', '--json', '--nthread', '6']  # TODO
#
#    psi4_subprocess(psi4rec)  # updates psi4rec
#
#    psi4rec = psi4rec['json']  # TODO NOT how this should be done figure out 1-tier/2-tier
#
#    try:
#        pass
#    except KeyError as err:
#        raise KeyError(
#            'Required fields missing from ({})'.format(jobrec.keys())) from err
#
#    try:
#        psi4rec['raw_output']
#    except KeyError as err:
#        raise KeyError('Required fields missing from ({})'.format(
#            psi4rec.keys())) from err
#
#    if not psi4rec['success']:
#        raise RuntimeError(psi4rec['error'])
#
#    for fl in psi4rec.keys():
#        if fl.startswith('outfile_'):
#            jobrec[fl] = psi4rec[fl]
#
#    # Absorb results into qcdb data structures
#    progvars = PreservingDict(psi4rec['extras']['qcvars'])
#    for k in list(progvars.keys()):
#        if k in ['DETCI AVG DVEC NORM', 'MCSCF TOTAL ENERGY']:
#            progvars.pop(k)
#
#    qcvars.build_out(progvars)
#    calcinfo = qcvars.certify(progvars, plump=True, nat=len(jobrec['molecule']['mass']))
#
#    jobrec['raw_output'] = psi4rec['raw_output']
#    jobrec['qcvars'] = calcinfo
#
#    jobrec['provenance'].append(psi4rec['provenance'])
#
#    return jobrec

def run_psi4(name, molecule, options, **kwargs):

    resi = ResultInput(
        **{
            'driver': inspect.stack()[1][3],
            'extras': {
                'qcdb:options': copy.deepcopy(options),
            },
            'model': {
                'method': name,
                'basis': '(auto)',
            },
            'molecule': molecule.to_schema(dtype=2),
            'provenance': provenance_stamp(__name__),
        })

    jobrec = qcng.compute(resi, "qcdb-psi4", raise_error=True).dict()
    hold_qcvars = jobrec['extras'].pop('qcdb:qcvars')
    jobrec['qcvars'] = {key: qcel.Datum(**dval) for key, dval in hold_qcvars.items()}
    pp.pprint(jobrec)
    print(jobrec.keys())
    print(jobrec['success'])
    return jobrec

def _print_helper(label, dicary, do_print):
    if do_print:
        print(label + ' <<<')
        pp.pprint(dicary)
        print('>>>')

class QcdbPsi4Harness(Psi4Harness):

    def compute(self, input_model: 'ResultInput', config: 'JobConfig') -> 'Result':
        self.found(raise_error=True)

        verbose = 1
        _print_helper(f'[1] {self.name} RESULTINPUT PRE-PLANT', input_model.dict(), verbose >= 3)

        input_data = self.qcdb_build_input(input_model, config)
        input_model = ResultInput(**input_data)

        _print_helper(f'[2] {self.name} RESULTINPUT PRE-ENGINE', input_model.dict(), verbose >= 4)

        output_model = Psi4Harness.compute(self, input_model=input_model, config=config)

        _print_helper(f'[3] {self.name} RESULT POST-ENGINE', output_model.dict(), verbose >= 4)

        # ???
        if not output_model.success:
            return output_model

        _print_helper(f'[4a] {self.name} RESULT POST-HARVEST', output_model.dict(), verbose >= 5)

        output_model = self.qcdb_post_parse_output(input_model, output_model)

        _print_helper(f'[4] {self.name} RESULT POST-POST-HARVEST', output_model.dict(), verbose >= 2)

        return output_model

    def qcdb_build_input(self, input_model: 'ResultInput', config: 'JobConfig',
                         template: Optional[str] = None) -> Dict[str, Any]:
        input_data = input_model.dict()

        molrec = qcel.molparse.from_schema(input_model.molecule.dict())
        ropts = input_model.extras['qcdb:options']

        muster_inherited_options(ropts)
        mtd = input_data['model']['method']
        mtd =  mtd[3:] if mtd.startswith('p4-') else mtd
        input_data['model']['method'] = mtd

        # should we put this memory in the JobConfig object? I don't think the units agree
        omem = ropts.scroll['QCDB'].pop('MEMORY')
        #print(config.memory, '!!')
        #config.memory = omem.value #???
        #print(config.memory, '!!')

        input_data['extras'] = {'wfn_qcvars_only' : True}
        #input_data['kwargs'] = jobrec['kwargs']
        #input_data['return_output'] = True

        popts = {}
        for k, v in ropts.scroll['QCDB'].items():
            if v.disputed():
                popts[k] = v.value

        for k, v in ropts.scroll['PSI4'].items():
            if v.disputed():
                popts[k] = v.value
        input_data['keywords'] = popts

        if 'BASIS' in input_data['keywords']:
            input_data['model']['basis'] =  input_data['keywords']['BASIS']

        return input_data

    def qcdb_post_parse_output(self, input_model: 'ResultInput', output_model: 'Result') -> 'Result':

        dqcvars = PreservingDict(copy.deepcopy(output_model.extras['qcvars']))
        qcvars.build_out(dqcvars)
        calcinfo = qcvars.certify(dqcvars, plump=True, nat=len(output_model.molecule.symbols))
        output_model.extras['qcdb:qcvars'] = calcinfo

        return output_model

max_stars_repo_name: nuwandesilva/qcdb

max_issues_repo_name: nuwandesilva/qcdb

max_forks_repo_name: nuwandesilva/qcdb


text:
import os
import argparse
import matplotlib.pyplot as plt
from numpy.core.multiarray import empty
import pandas as pd
import numpy as np
import sys
from keras.models import load_model
from keras.backend import clear_session
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from keras.layers import *
from sklearn.preprocessing import MinMaxScaler

if len(sys.argv) < 5:
    print("Too few arguments")
    quit()
elif len(sys.argv) > 5:
    print("Too many arguments")
    quit()

# Read the arguments
for i in range (1, len(sys.argv)):
    if(sys.argv[i] == "-d"):
        dataset_name = sys.argv[i+1]
    elif(sys.argv[i] == "-n"):
        num_time_series = int(sys.argv[i+1])

csv_path = os.path.join(os.path.abspath(__file__), "../../../dir/")
csv_path = os.path.join(csv_path, dataset_name)

model_path = os.path.join(os.path.abspath(__file__), "../../../models/")

# Read the input file
df = pd.read_csv(csv_path, header=None, delimiter='\t')
file_ids = df.iloc[:, [0]].values
df = df.drop(df.columns[0], axis=1)
df = df.transpose()

# Training data: 80%, Test data: 20%
train_size = int(len(df) * 0.80)
test_size = len(df) - train_size

sc = MinMaxScaler(feature_range = (0, 1))

# For each time series get a prediction
for step in range (0, num_time_series):
    X_train = []
    y_train = []

    training_set = df.iloc[:train_size, [step]].values
    training_set_scaled = sc.fit_transform(training_set)

    for i in range(60, train_size):
        X_train.append(training_set_scaled[i-60:i, 0])
        y_train.append(training_set_scaled[i, 0])

    X_train, y_train = np.array(X_train), np.array(y_train)
    X_train = np.reshape(X_train, (X_train.shape[0], X_train.shape[1], 1))

    dataset_train = df.iloc[:train_size, [step]]
    dataset_test = df.iloc[train_size:, [step]]
    dataset_total = pd.concat((dataset_train, dataset_test), axis = 0)

    inputs = dataset_total[len(dataset_total) - len(dataset_test) - 60:].values
    inputs = inputs.reshape(-1,1)
    inputs = sc.transform(inputs)
    X_test = []
    for i in range(60, test_size+60):
        X_test.append(inputs[i-60:i, 0])

    X_test = np.array(X_test)
    X_test = np.reshape(X_test, (X_test.shape[0], X_test.shape[1], 1))

    # Load model for this time series
    model = load_model(os.path.join(model_path, "Forecast_model" + str(step) +".h5"))

    # Get the prediction
    predicted_stock_price = model.predict(X_test)
    predicted_stock_price = sc.inverse_transform(predicted_stock_price)

    # Visualising the results
    plt.plot(dataset_test.values, color = 'red', label = 'Real values')
    plt.plot(predicted_stock_price, color = 'blue', label = 'Predicted values')
    plt.title('Time Series Prediction')
    plt.xlabel('Time')
    plt.ylabel('Value')
    plt.legend()
    fig = plt.gcf()
    fig.set_size_inches(18.5, 10.5)
    plt.savefig(os.path.join(os.path.abspath(__file__), "../../../dir/exports/q1/var1/graph"+ str(step) +".png"))
    plt.clf()

    del model
    clear_session()

reason: import numpy,from numpy


text:
import numpy
from astropy import units as u
from astropy.constants import c, G, M_sun, R_sun, au
from matplotlib import pyplot as plt
from PyCom import integrated_flux, b_of_z

wavelength = 10**-6  # [m]
path = ''
list_b0 = []
list_flux = []
new_d0 = 1
base_flux = integrated_flux(wavelength=wavelength, z=(550 * au) / u.meter, d0=new_d0)
for i in range(546, 2200):
    z = (i * au) / u.meter
    new_flux = integrated_flux(wavelength=wavelength, z=z, d0=new_d0)
    b = b_of_z(z)
    list_b0.append(b)
    list_flux.append(new_flux / base_flux)

plt.rc('font',  family='serif', serif='Computer Modern Roman')
plt.rc('text', usetex=True)
size = 4.5
aspect_ratio = 1.5
plt.figure(figsize=(size, size / aspect_ratio))
ax = plt.gca()
ax.get_yaxis().set_tick_params(direction='out')
ax.get_xaxis().set_tick_params(direction='out')
ax.get_yaxis().set_tick_params(which='both', direction='out')
ax.get_xaxis().set_tick_params(which='both', direction='out')
ax.set_xlabel(r'$R$ $(R_{\odot})$')
ax.set_ylabel(r'Relative flux')
ax.set_xlim(1, 2)
ax.set_ylim(1, 2)
plt.plot(list_b0, list_flux, color='black')
ay2 = ax.twiny()
ay2.get_xaxis().get_major_formatter().set_useOffset(False)
ay2.get_xaxis().get_major_formatter().set_scientific(False)
ay2.get_xaxis().set_tick_params(direction='out')
ay2.get_xaxis().set_tick_params(direction='out')
ay2.get_xaxis().set_tick_params(which='both', direction='out')
ay2.get_xaxis().set_tick_params(which='both', direction='out')
ay2.set_xlim([546, 2184])
ay2.set_xlabel(r'Heliocentric distance $z$ (au)')
plt.savefig(path + 'figure_z.pdf', bbox_inches='tight')

max_stars_repo_name: hippke/communication

max_issues_repo_name: hippke/communication

max_forks_repo_name: hippke/communication

reason: import numpy,from astropy


text:
import torch
import torch.nn as nn
from torch.autograd import Variable
import numpy as np

def get_rnn(rnn_type, ninputs, nhid, nlayers, dropout=0, bidirectional=False):
    if rnn_type in ['LSTM', 'GRU']:
        return getattr(nn, rnn_type)(ninputs, nhid, nlayers, dropout=dropout,
                                     bidirectional=bidirectional)
    else:
        try:
            nonlinearity = {'RNN_TANH': 'tanh', 'RNN_RELU': 'relu'}[rnn_type]
        except KeyError:
            raise ValueError("An invalid option for `rnn_type` was supplied, "
                             "options are ['LSTM', 'GRU', 'RNN_TANH' or "
                             "'RNN_RELU']")
        return nn.RNN(ninputs, nhid, nlayers, nonlinearity=nonlinearity,
                      dropout=dropout, bidirectional=bidirectional)

class Encoder(nn.Module):
    """Encodes text x of style y into content vector z"""
    def __init__(self, nhid, embeddings, rnn_type='LSTM',
                 nlayers=1, dropout=0.5):
        super(Encoder, self).__init__()

        self.embeddings = embeddings
        self.drop = nn.Dropout(dropout)
        self.rnn = get_rnn(rnn_type, self.embeddings.embedding_dim,
                           nhid, nlayers, dropout)

        self.rnn_type = rnn_type
        self.nhid = nhid
        self.initial_state_dim = nhid * nlayers
        if rnn_type == 'LSTM':
            self.initial_state_dim *= 2

    def forward(self, style, data, seq_lens):
        """
        Here style, encoded style of a given text,
        serves as initial hidden state of encoder.
        """
        emb = self.drop(self.embeddings(data))
        output, _ = self.rnn(emb, style)

        _, batch_size, hidden_dim = output.size()
        arange = seq_lens.data.new(*seq_lens.size()).copy_(
                torch.arange(0, batch_size))

        last_states = output[seq_lens.data - 1, arange]
        return last_states

    def split_initial_hidden(self, initial_hidden):
        initial_hidden = initial_hidden.view(
        initial_hidden.size(0), -1, self.nhid)
        initial_hidden = initial_hidden.transpose(0, 1).contiguous()
        if self.rnn_type == 'LSTM':
            initial_hidden = initial_hidden.chunk(2)
        return initial_hidden

class Generator(nn.Module):
    """Generates text x based on content z and style y"""
    def __init__(self, nhid, embeddings, eos_id, gamma=0.001,
                 rnn_type='LSTM', nlayers=1, dropout=0.5):
        super(Generator, self).__init__()

        self.embeddings = embeddings
        self.drop = nn.Dropout(dropout)
        self.rnn = get_rnn(rnn_type, self.embeddings.embedding_dim,
                           nhid, nlayers, dropout)
        self.projection = nn.Linear(nhid, self.embeddings.num_embeddings)
        self.softmax = nn.Softmax()
        self.gamma = gamma
        self.rnn_type = rnn_type
        self.eos = eos_id
        self.nhid = nhid

        self.initial_state_dim = nhid * nlayers
        if rnn_type == 'LSTM':
            self.initial_state_dim *= 2

    def forward(self, init_hidden, data, teacher_forcing=True):
        hidden = init_hidden

        if teacher_forcing:
            emb = self.drop(self.embeddings(data))
            outputs, hidden = self.rnn(emb, hidden)
        else: # sampling
            outputs = []
            max_len, batch_size = data.size()
            eos = Variable(data.data.new([self.eos]))
            emb = self.embeddings(eos).expand(
                    [1, batch_size, self.embeddings.embedding_dim])

            for step in range(max_len):
                emb = self.drop(emb)
                output, hidden = self.rnn(emb, hidden)
                outputs.append(output)

                decoded = self.projection(
                        output.view([-1, output.size(2)]))

                emb = torch.matmul(
                        self.softmax(decoded/self.gamma), self.embeddings.weight)
                emb = emb.view([1, batch_size, self.embeddings.embedding_dim])

            outputs = torch.cat(outputs, 0)

        outputs = self.drop(outputs)
        decoded = self.projection(
            outputs.view(outputs.size(0) * outputs.size(1), outputs.size(2)))
        decoded = decoded.view(
            outputs.size(0), outputs.size(1), decoded.size(1))

        return decoded, outputs, hidden

    def split_initial_hidden(self, initial_hidden):
        initial_hidden = initial_hidden.view(
        initial_hidden.size(0), -1, self.nhid)
        initial_hidden = initial_hidden.transpose(0, 1).contiguous()
        if self.rnn_type == 'LSTM':
            initial_hidden = initial_hidden.chunk(2)
        return initial_hidden

class Discriminator(nn.Module):
    def __init__(self, dim, filter_sizes, n_filters,
            activation=nn.LeakyReLU(), dropout=0.5):
        super(Discriminator, self).__init__()

        self.activation = activation
        self.drop = nn.Dropout(dropout)
        self.projection = nn.Linear(len(filter_sizes)*n_filters, 1)

        convolutions = []
        for f_size in filter_sizes:
            convolutions.append(nn.Conv2d(1, n_filters, [f_size, dim], padding=(2, 0)))

        self.convolutions = nn.ModuleList(convolutions)

    def forward(self, data):
        data = data.transpose(0, 1).unsqueeze(1)
        data = self.drop(data)
        outputs = []
        for conv in self.convolutions:
            hid = self.activation(conv(data))
            # maxpool over time:
            # Note that torch.max over a dimension
            # returns tuple of tensors (max, argmax)
            pooled = torch.max(hid, 2)[0].contiguous()
            outputs.append(pooled.view([-1, conv.out_channels]))

        logits = self.projection(self.drop(torch.cat(outputs, dim=1)))
        return logits.squeeze(-1)

class Model(nn.Module):
    def __init__(self, ntokens, style_dim,
                 encoder_emb_dim, generator_emb_dim,
                 encoder_kwargs, generator_kwargs,
                 discriminator_kwargs, dropout=0.5,
                 lmb=1.0, tie_embeddings=True):
        super(Model, self).__init__()

        generator_embeddings = nn.Embedding(ntokens, generator_emb_dim)

        if tie_embeddings:
#            if encoder_kwargs["nhid"] != generator_kwargs["nhid"]:
            if encoder_emb_dim != generator_emb_dim:
                raise ValueError("When using tied flag hidden states' dimenstions of encoder and generator must be the same")
            encoder_embeddings = generator_embeddings
        else:
            encoder_embeddings = nn.Embedding(ntokens, encoder_emb_dim)

        self.lmb = lmb
        self.drop = nn.Dropout(dropout)
        self.style_encoder = nn.Embedding(2, style_dim)
        self.generator = Generator(embeddings=generator_embeddings, **generator_kwargs)
        self.encoder = Encoder(embeddings=encoder_embeddings, **encoder_kwargs)

        self.encoder_init_projection = nn.Linear(
                style_dim, self.encoder.initial_state_dim)
        self.generator_init_projection = nn.Linear(
                style_dim+encoder_kwargs["nhid"],
                self.generator.initial_state_dim)

        self.discriminator = nn.ModuleList([
                Discriminator(dim=generator_kwargs['nhid'],
                    **discriminator_kwargs)
                for i in range(2)
                ])

        self.rec_criterion = nn.CrossEntropyLoss(
                size_average=False, ignore_index=-1)
        self.adv_criterion = nn.BCEWithLogitsLoss()

    def get_style_encoding(self, style_id, data):
        style = data.data.new([style_id])
        style = self.style_encoder(Variable(style))
        style = style.expand([data.size(1),
            self.style_encoder.embedding_dim])

        return style

    def mask_invalid_hiddens(self, hiddens, seq_lens):
        arange = seq_lens.data.new(hiddens.size(0)).copy_(
                torch.arange(0, hiddens.size(0)))
        arange = Variable(arange)
        mask = (arange.view([-1, 1]) < seq_lens).float().unsqueeze(2)
        return hiddens * mask

    def transfer_style(self, data, in_style, out_style):
        assert data.size()[1:] == torch.Size([1])
        self.eval()

        in_style = self.get_style_encoding(in_style, data)
        out_style = self.get_style_encoding(out_style, data)

        encoder_init = self.encoder.split_initial_hidden(
                self.encoder_init_projection(in_style))

        emb = self.encoder.embeddings(data)
        output, _ = self.encoder.rnn(emb, encoder_init)
        content = output[-1]

        generator_init = self.generator_init_projection(
                torch.cat([out_style, content], 1))
        generator_init = self.generator.split_initial_hidden(
                generator_init)

        transferred = []
        hidden = generator_init
        eos = Variable(data.data.new([self.generator.eos]))
        emb = self.generator.embeddings(eos).expand(
                [1, 1, self.generator.embeddings.embedding_dim])

        for step in range(80):
            output, hidden = self.generator.rnn(emb, hidden)

            decoded = self.generator.projection(
                    output.view([-1, output.size(2)])).max(-1)[1]

            transferred.append(decoded.data[0])
            if transferred[-1] == self.generator.eos:
                break

            emb = self.generator.embeddings(decoded)
            emb = emb.view([1, 1, self.generator.embeddings.embedding_dim])

        return transferred

    def compute_losses(self, data, targets, seq_lens):
        assert (data[0].size(1) == data[1].size(1))

        batch_size = data[0].size(1)

        style = [self.get_style_encoding(0, data[0]),
                self.get_style_encoding(1, data[1])]

        labels = [
            Variable(data[0].data.new([0])).expand([batch_size]).float(),
            Variable(data[1].data.new([1])).expand([batch_size]).float()
            ]

        rec_loss = 0.0
        adv_loss = [0.0, 0.0]

        for (p, q) in [(0,1), (1,0)]:
            encoder_init = self.encoder_init_projection(self.drop(style[p]))
            encoder_init = self.encoder.split_initial_hidden(encoder_init)
            content = self.encoder(encoder_init, data[p], seq_lens[p])

            generator_init_true = self.generator_init_projection(
                    self.drop(torch.cat([style[p], content], 1)))
            generator_init_true = self.generator.split_initial_hidden(
                    generator_init_true)
            output, hiddens_true, _ = self.generator(
                    generator_init_true, data[p], teacher_forcing=True)

            generator_init_false = self.generator_init_projection(
                    self.drop(torch.cat([style[q], content], 1)))
            generator_init_false = self.generator.split_initial_hidden(
                    generator_init_false)
            _, hiddens_false, _ = self.generator(
                    generator_init_false, data[p], teacher_forcing=False)

            output_flat = output.view(-1, output.size(-1))
            targets_flat = targets[p].view(-1)

            rec_loss += self.rec_criterion(output_flat, targets_flat)

            hiddens_true = self.mask_invalid_hiddens(
                    hiddens_true, seq_lens[p])
            hiddens_false = self.mask_invalid_hiddens(
                    hiddens_false, seq_lens[p])

            if self.generator.rnn_type == "LSTM":
                init_true = generator_init_true[0]
                init_false = generator_init_false[0]
            else:
                init_true = generator_init_true
                init_false = generator_init_false

            hiddens_true = torch.cat([init_true, hiddens_true], 0)
            hiddens_false = torch.cat([init_false, hiddens_false], 0)

            adv_loss[p] += self.adv_criterion(
                    self.discriminator[p](hiddens_true),
                    labels[1])

            adv_loss[q] += self.adv_criterion(
                    self.discriminator[q](hiddens_false),
                    labels[0])

        return rec_loss, adv_loss[0], adv_loss[1]

    def eval_on(self, batch_iterator):
        # Turn on evaluation mode which disables dropout.
        self.eval()

        total_rec = 0.0
        total_adv0 = 0.0
        total_adv1 = 0.0
        total_words = 0.0
        for batch_no, (data, targets, seq_lens) in enumerate(batch_iterator):
            rec, adv0, adv1 = self.compute_losses(data, targets, seq_lens)
            total_words += sum(
                (t.view(-1) != self.rec_criterion.ignore_index).float().sum()
                for t in targets).data[0]

            total_rec += rec.data[0]
            total_adv0 += adv0.data[0] if type(adv0) is Variable else adv0
            total_adv1 += adv1.data[0] if type(adv1) is Variable else adv1

        total_rec /= total_words
        total_adv0 /= (batch_no + 1.0)
        total_adv1 /= (batch_no + 1.0)

        print (90 * '#')
        print ("End of epoch: pplx: %f, rec_loss: %f, ae_loss: %f, adv0: %f, adv1: %f" %
                (np.exp(total_rec), total_rec,
                    total_rec - self.lmb * (total_adv0 + total_adv1),
                total_adv0, total_adv1))
        print (90 * '#')

        return {'total_ae_loss':
                total_rec - self.lmb * (total_adv0 + total_adv1),
                'total_adv0' : total_adv0,
                'total_adv1' : total_adv1,
                }

    def train_on(self, batch_iterator, optimizer_step):
        # Turn on training mode which enables dropout.
        self.train()

        for batch_no, (data, targets, seq_lens) in enumerate(batch_iterator):
            rec, adv0, adv1 = self.compute_losses(data, targets, seq_lens)
            batch_words = sum(
                (t.view(-1) != self.rec_criterion.ignore_index).float().sum()
                for t in targets).data[0]

            rec /= batch_words
            optimizer_step(rec, adv0, adv1, batch_no)

            rec = rec.data[0]
            adv0 = adv0.data[0] if type(adv0) is Variable else adv0
            adv1 = adv1.data[0] if type(adv1) is Variable else adv1

            if batch_no % 250 == 0:
                print ("batch %d: pplx: %f, rec_loss: %f, ae_loss: %f, adv0: %f, adv1: %f" %
                        (batch_no, np.exp(rec), rec,
                            (rec - self.lmb * (adv0 + adv1)),
                            adv0, adv1))


text:
# coding=utf-8
# @author:
# - Meiling Chen, brookechen@sjtu.edu.cn
# select outliers with qp or mip model

from numpy import inf

def qp_gurobipy(x, y0, n, m, N):

    from gurobipy import Model, GRB, quicksum

    # set lamda
    lamda = 1

    # generate model
    model = Model('linear regression_qp')

    # variables
    r = [model.addVar(name='r({})'.format(i), vtype=GRB.CONTINUOUS) for i in range(n)]
    q = [model.addVar(name='q({})'.format(i), vtype=GRB.BINARY) for i in range(n)]
    b = [model.addVar(name='b({})'.format(j), lb=-inf, vtype=GRB.CONTINUOUS) for j in range(m)]

    # objective
    model.setObjective(quicksum(r[i] * q[i] + lamda * (r[i]*r[i] + q[i]*q[i]) for i in range(n)), GRB.MINIMIZE)

    # constraints
    for i in range(n):
        model.addConstr(r[i] >= y0[i] - quicksum(x[j, i] * b[j] for j in range(m)))
        model.addConstr(r[i] >= quicksum(x[j, i] * b[j] for j in range(m)) - y0[i])
    model.addConstr(quicksum(q[i] for i in range(n)) >= N)

    # optimize
    model.optimize()

    # get results
    r_sol = [r[i].x for i in range(n)]
    q_sol = [round(q[i].x) for i in range(n)]
    b_sol = [b[j].x for j in range(m)]

    return r_sol, q_sol, b_sol

def mip_gurobipy(x, y0, n, m, N):

    from gurobipy import Model, GRB, quicksum

    # set big M
    M = 100

    # generate model
    model = Model('linear regression_mip')

    # variables
    r = [model.addVar(name='r({})'.format(i), vtype=GRB.CONTINUOUS) for i in range(n)]
    t = [model.addVar(name='t({})'.format(i), vtype=GRB.CONTINUOUS) for i in range(n)]
    q = [model.addVar(name='q({})'.format(i), vtype=GRB.BINARY) for i in range(n)]
    b = [model.addVar(name='b({})'.format(j), lb=-inf, vtype=GRB.CONTINUOUS) for j in range(m)]

    # objective
    model.setObjective(quicksum(r[i] - t[i] for i in range(n)), GRB.MINIMIZE)

    # constraints
    for i in range(n):
        model.addConstr(r[i] >= y0[i] - quicksum(x[j, i] * b[j] for j in range(m)))
        model.addConstr(r[i] >= quicksum(x[j, i] * b[j] for j in range(m)) - y0[i])
        model.addConstr(t[i] <= M * (1 - q[i]))
        model.addConstr(t[i] <= r[i])
    model.addConstr(quicksum(q[i] for i in range(n)) >= N)

    # optimize
    model.optimize()

    # get results
    r_sol = [r[i].x for i in range(n)]
    q_sol = [round(q[i].x) for i in range(n)]
    b_sol = [b[j].x for j in range(m)]
    t_sol = [t[i].x for i in range(n)]

    return r_sol, q_sol, b_sol

def mip_mip(x, y0, n, m, N, solver_name):

    from mip import Model, xsum, minimize

    # set big M
    M = 100

    # generate model
    model = Model(name='linear regression_mip', solver_name=solver_name)
    model.verbose = 0

    # variables
    r = [model.add_var(name='r({})'.format(i), var_type='C') for i in range(n)]
    t = [model.add_var(name='t({})'.format(i), var_type='C') for i in range(n)]
    q = [model.add_var(name='q({})'.format(i), var_type='B') for i in range(n)]
    b = [model.add_var(name='b({})'.format(j), lb=-inf, var_type='C') for j in range(m)]

    # objective
    model.objective = minimize(xsum(r[i] - t[i] for i in range(n)))

    # constraints
    for i in range(n):
        model += r[i] >= y0[i] - xsum(x[j, i] * b[j] for j in range(m))
        model += r[i] >= xsum(x[j, i] * b[j] for j in range(m)) - y0[i]
        model += t[i] <= M * (1 - q[i])
        model += t[i] <= r[i]
    model += xsum(q[i] for i in range(n)) >= N

    # optimize
    model.optimize()

    # get results
    r_sol = [r[i].x for i in range(n)]
    q_sol = [round(q[i].x) for i in range(n)]
    b_sol = [b[j].x for j in range(m)]
    t_sol = [t[i].x for i in range(n)]

    return r_sol, q_sol, b_sol

max_stars_repo_name: brentian/artoopt

max_issues_repo_name: brentian/artoopt

max_forks_repo_name: brentian/artoopt

reason: from numpy


text:
import numpy as np
import pandas as pd
import h5py

from sklearn.feature_extraction.text import CountVectorizer
from scipy.sparse import save_npz, load_npz

from tools import pad_integers

# Importing the data
indir = 'C:/data/addm/'
corpus = pd.read_csv(indir + 'corpus_with_lemmas_clean.csv')
text = corpus.dx
target = corpus.aucaseyn

# Vectorizing the text
vec = CountVectorizer(ngram_range=(1, 2), binary=False)
X = vec.fit_transform(text)

# Saving the docterm matrix as npz
save_npz(indir + 'doctermat.npz', X)

# Making a binary unigram doctermat for fastText
univec = CountVectorizer(ngram_range=(1, 1), binary=True)
X = univec.fit_transform(text)
save_npz(indir + 'unigram_doctermat.npz', X)

# Converting the docterm matrix to bag-of-integers for fastText
indices = [np.nonzero(row)[1] for row in X]
max_length = np.max([len(doc) for doc in indices])
int_sents = np.array([pad_integers(doc, max_length) for doc in indices])

# Writing the integer sequences to HDF5
out = h5py.File(indir + 'unigram_ints.hdf5', mode='w')
out['sents'] = int_sents
out.close()

# Making a compressed bigram matrix for fastText
bivec = CountVectorizer(ngram_range=(1, 2), binary=True, min_df=5)
X = bivec.fit_transform(text)
save_npz(indir + 'bigram_doctermat.npz', X)

# Converting the docterm matrix to bag-of-integers for fastText
indices = [np.nonzero(row)[1] for row in X]
max_length = np.max([len(doc) for doc in indices])
int_sents = np.array([pad_integers(doc, max_length) for doc in indices])

# Writing the integer sequences to HDF5
out = h5py.File(indir + 'bigram_ints.hdf5', mode='w')
out['sents'] = int_sents
out.close()


text:
# Title: General Demo file
# Description: Demo testing on spiked covariance different online PSP algorithms
# Author: Victor Minden (vminden@flatironinstitute.org)and Andrea Giovannucci (agiovannucci@flatironinstitute.org)
# Notes: Adapted from code by Andrea Giovannucci and Cengiz Pehlevan
# Reference: (Cardot and Degras, 2015)
##############################
import numpy as np
import pylab as pl
import time

from online_psp.incremental_pca import IPCA
from online_psp.ccipca import CCIPCA
from online_psp.fast_similarity_matching import FSM
from online_psp.similarity_matching import SM

from online_psp.util import subspace_error, generate_samples

# Simulation parameters
K = 10
n_epoch = 5
err_its = 64

scale_data = True
init_ortho = True

D, N = 50, 10000
X, U, sigma2 = generate_samples(K, N, D, method='spiked_covariance', scale_data=scale_data)

# %% initialization
lambda_1 = 1e-8 * np.ones(shape=(K,))
Uhat0 = X[:, :K] / np.sqrt((X[:, :K] ** 2).sum(0))

if init_ortho:
    # Optionally orthogonalize the initial guess
    Uhat0, _ = np.linalg.qr(Uhat0)

# Build fitters of each type
ccipca = CCIPCA(K, D, Uhat0=Uhat0, sigma2_0=lambda_1,
                cython='auto')
lambda_1 *= 0
ipca = IPCA(K, D, Uhat0=Uhat0, sigma2_0=lambda_1)

scal = 100
lr = lambda t: 1/(t + 5)
fsm = FSM(K, D, W0=Uhat0.T / scal,
                Minv0=scal * np.eye(K),
                learning_rate=lr)
sm = SM(K, D, W0=Uhat0.T / scal,
                M0= np.eye(K) / scal,
                learning_rate=lr)

algorithms = {'IPCA': ipca, 'CCIPCA': ccipca, 'SM': sm, 'FSM': fsm}

times = {}
errs  = {}

for name, algo in algorithms.items():
    print('Starting algorithm %s' % name)
    err = []
    time_1 = time.time()
    for _ in range(n_epoch):
        for its, x in enumerate(X.T):
            algo.fit_next(x)
            if its % err_its == 0:
                err.append(subspace_error(algo.get_components(), U[:, :K]))
    time_2 = time.time() - time_1
    errs[name] = err
    times[name] = time_2

# DISPLAY RESULTS
keys = list(algorithms.keys())
keys.sort()

for name in keys:
    pl.loglog(errs[name])
    pl.ylabel('relative subspace error (pop.)')
    pl.xlabel('samples')

    print('Final subspace error ' + name + ':' + str(errs[name][-1]))

pl.legend(keys)
pl.show()

max_stars_repo_name: eWaterCycle/streamingDataAssimilation

max_issues_repo_name: eWaterCycle/streamingDataAssimilation

max_forks_repo_name: eWaterCycle/streamingDataAssimilation


text:
import multiprocessing as mp
import queue
import time
import numpy as np

def my_fn(x):
    time.sleep(np.random.rand())
    print(-x)
    return x

def worker(fn, jobs_q, results_q):

    for args in iter(jobs_q.get, None):

        res = fn(*args)
        results_q.put(res)

def apply_async_ll(fn, iterable, nprocs=2):

    jobs_q = mp.Queue(maxsize=nprocs)
    results_q = mp.Queue()

    procs = [mp.Process(target=worker, args=(fn, jobs_q, results_q)) for _ in range(nprocs)]
    for proc in procs:
        proc.start()

    it = iter(iterable)
    results = []

    consumed = True
    num_in, num_out = 0, 0
    while True:

        try:
            if consumed:
                x = next(it)
                consumed = False
            jobs_q.put(x, block=False)
        except queue.Full:
            res = results_q.get(block=True)
            results.append(res)
            num_out += 1
        except StopIteration:
            break
        else:
            num_in += 1
            consumed = True

    for _ in range(nprocs):
        jobs_q.put(None, block=True)

    # Need this bit since the process might not exit if there is still data left in the results_q. When objects in the
    # queue are large they are buffered. The process won't exit until the buffer is flushed
    # https://stackoverflow.com/questions/26025486/python-processes-not-joining
    while num_out < num_in:
        res = results_q.get(block=True)
        results.append(res)
        num_out += 1

    for proc in procs:
        proc.join()

    return results

if __name__ == '__main__':

    apply_async_ll(my_fn, ([x] for x in range(10)))


text:
import numpy as np, h5py
import pandas as pd
import argparse
if __name__=="__main__":
    psr = argparse.ArgumentParser()
    psr.add_argument('-i', dest='ipt', nargs='+', help='input h5 file')
    psr.add_argument('-o', dest='opt', help='output h5 file')
    psr.add_argument('-c', dest='channel', default=2, type=int, help='channel for test pmt')
    psr.add_argument('--ref', default=3, type=int, help='channel for refer pmt')
    args = psr.parse_args()
    testR, refR = 1, 1
    for f in args.ipt:
        with h5py.File(f, 'r') as ipt:
            QEinfo = pd.DataFrame(ipt['QEinfo'][:])
            QEinfo = QEinfo.set_index('ch')
            testR *= QEinfo.loc[args.channel]['ratio']
            refR *= QEinfo.loc[args.ref]['ratio']
    print(np.sqrt(testR/refR))

max_stars_repo_name: greatofdream/pmtTest

max_issues_repo_name: greatofdream/pmtTest

max_forks_repo_name: greatofdream/pmtTest


text:
# -*- coding: utf-8 -*-
"""
Created on Sun Dec  8 22:37:00 2019

@author: for_y
"""

import numpy as np
from AnnoDomini.hamilton_mc import HMC, describe
import sys
import os

sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))

def test_hmc():
    def norm_function(x):
        var = 1
        denom = (2*np.pi*var)**.5
        num = np.exp(-x**2/2)
        return num/denom

    start_point = 0
    chain,accepts_ratio = HMC(target_pdf = norm_function, burn_in=0, thinning=1,chain_len=100, q_init=[start_point],epsilon = 0.05)
    assert len(chain) == 100
    q = chain[:,0]
    assert max(q) < 10
    assert np.isnan(q).sum() == 0

    def neg_log_weibull(lam = 1, k = 0.5):
        def w(x):
            if x > 0:
                return -(np.log(k / lam) + (k-1) * np.log(x/lam) - (x/lam) ** k)
            else:
                return float('inf')
        return w

    start_point = 3
    func = neg_log_weibull(k = 1.5)
    chain,accepts_ratio = HMC(U = func, burn_in=0, thinning=1,chain_len=100, q_init=[start_point],epsilon = 0.02)
    assert len(chain) == 100
    q = chain[:,0]
    assert max(q) < 10
    assert min(q) > 0
    assert np.isnan(q).sum() == 0


text:
import argparse
import logging
import string
import codecs
import numpy as np
from tqdm import tqdm
from scipy.stats import pearsonr
from scipy.stats import spearmanr

logging.basicConfig(level=logging.DEBUG,
                    format='%(asctime)s - %(filename)s[line:%(lineno)d] - %(levelname)s: %(message)s')

parser = argparse.ArgumentParser()
parser.add_argument("--hyp_file", type=str, help="path to hypothesis file")
parser.add_argument("--ref_file", type=str, help="path to reference file")
parser.add_argument("--num_test", type=int, help="number of test cases")
parser.add_argument("--save_path",  type=str, help='path to save system level fm score')
args = parser.parse_args()

def calc_fm(hyp, ref):
    return max(1/hyp, 1/ref)/min(1/hyp, 1/ref)

if __name__=='__main__':

    logging.info("Reading hypothesis probs -------------------------------------------------------")
    with codecs.open(args.hyp_file, mode='r', encoding='utf-8') as rf:
        hyp_sent_ppl = rf.read()
    hyp_sent_ppl = eval(hyp_sent_ppl)
    assert len(hyp_sent_ppl) == args.num_test, 'number of hypotheses not equal to {}'.format(args.num_test)

    logging.info("Reading references probs -------------------------------------------------------")
    with codecs.open(args.ref_file, mode='r', encoding='utf-8') as rf:
        ref_sent_ppl = rf.read()
    ref_sent_ppl = eval(ref_sent_ppl)
    assert len(ref_sent_ppl) == args.num_test, 'number of references not equal to {}'.format(args.num_test)

    logging.info("computing FM scores -------------------------------------------------------")
    # calculate FM score
    full_scores = []
    for hyp, ref in zip(hyp_sent_ppl, ref_sent_ppl):
        score = calc_fm(hyp, ref)
        full_scores.append(score)

    with codecs.open(args.save_path, mode='w', encoding='utf-8') as wf:
        for score in full_scores:
            wf.write(str(score) + '\n')

    logging.info("Done writing FM scores to {}-------------------------------------------------------".format(args.save_path))

max_stars_repo_name: AliOsm/deep-amfm

max_issues_repo_name: AliOsm/deep-amfm

max_forks_repo_name: AliOsm/deep-amfm