text:
# Copyright (C) 2014 Andreas M. Weller <andreas.m.weller@gmail.com>
#
# read a bedtools output file from
#
# bedtools coverage -abam Q2PL2_H01_N.bam -b TSB_148_gene_panel_HP_amplicons.bed -d > test_coverage.csv
#
# and find bases with coverage or strand_ratio below the threshold

# general modules
import pandas as pd
import sys
import subprocess
import os
import matplotlib.pylab as plt
import numpy as np
import pprint
import logging
import re

import seaborn as sns
sns.set(font="serif") # prevents error message about fonts

# personal modules from the same folder
import plot_exon_coverage_v2 as plotting
import plot_exon_coverage_all_samples as all_sample_plotting
import unite_coverage_files as UniteCoverage
import CoverageCheckClasses
from CoverageCheckConfig import *

script_folder = os.path.dirname(os.path.realpath(__file__)) + "/"

DISTNAME = 'coveragecheck'
MAINTAINER = 'Andreas Weller'
MAINTAINER_EMAIL = 'andreas.m.weller@gmail.com'
LICENSE = 'BSD (3-clause)'
DOWNLOAD_URL = 'https://github.com/aweller/CoverageCheck/'
VERSION = '0.1'

#####################################################################################################################

def parse_exons_into_dataframe_and_dict(exon_filename):
    """
    Parse a list of all human exons into 2 datastructures:
    1. exons (pandas df)
    2. exons_per_gene (dictionary of gene:[start, stop, exon_no])
    """

    #exon_filename = "/home/andreas/bioinfo/core/general/data/HumanExons_Ensembl_v65_merged.csv"
    #exon_filename = "/home/andreas/bioinfo/core/general/data/HumanExons_Ensembl_v75_all_genes_merged.csv"

    header = ["chrom", "exon_start", "exon_stop", "gene", "exon_no", "strand"]
    exons = pd.read_csv(exon_filename, sep="\t", names=header)
    exons["gene_upper"] = exons.gene.str.upper()
    exons = exons.sort(columns = ["gene", "exon_start", "exon_stop"])

    global exons_per_gene
    exons_per_gene = {}
    for _, row in exons.iterrows():
        gene = row["gene"].upper()
        start, stop = int(row["exon_start"]), int(row["exon_stop"])
        exon_no = row["exon_no"]

        if not exons_per_gene.get(gene):
            exons_per_gene[gene] = []

        exons_per_gene[gene].append((start, stop, exon_no))

    return exons, exons_per_gene

def parse_coverage_file_into_dataframe(coverage_file):
    """
    Parses the coverage file, transforms into a per-base matrix and returns a pandas DF
    """

    header = ["chr", "start", "stop", "amplicon", "na", "strand",  "amplicon_pos", "dp"]
    rawdf = pd.read_csv(coverage_file, sep="\t", names=header)
    rawdf["pos"] = rawdf.start +rawdf.amplicon_pos
    rawdf["chrompos"] = rawdf.apply(lambda x : "\t".join([str(x["chr"]), str(x["pos"]) ]), axis = 1 )
    rawdf["gene"] = rawdf.apply(lambda x : x["amplicon"].split("_")[0], axis = 1 )

    # df: per base

    df = dict(chrom = rawdf.chr.groupby(rawdf.chrompos).min(),
              pos = rawdf.pos.groupby(rawdf.chrompos).min(),
              gene = rawdf.gene.groupby(rawdf.chrompos).min(),
            start = rawdf.start.groupby(rawdf.chrompos).min(),
            stop = rawdf.stop.groupby(rawdf.chrompos).min(),
              minus_dp = rawdf[rawdf.strand == "-"].dp.groupby(rawdf.chrompos).max(),
              plus_dp = rawdf[rawdf.strand == "+"].dp.groupby(rawdf.chrompos).max(),
              dp = rawdf.dp.groupby(rawdf.chrompos).sum(),
    )
    df = pd.DataFrame(df).reset_index()

    def get_strand_ratios(row):
        minus = row["minus_dp"]
        plus = row["plus_dp"]

        if plus + minus == 0:
            return 0
        if plus == 0:
            return minus
        elif minus == 0:
            return plus
        else:
            larger = max([plus, minus])
            smaller = min([plus, minus])

            return abs(larger/float(smaller))

    df["strand_ratio"] = df.apply(get_strand_ratios, axis =1 )
    df = df.sort(columns = ["chrom", "pos"])

    return df

def parse_amplicons_into_df(bed, exons):
    """
    Parse the amplicons in the bedfile into a DF
    Get rid of amplicons that are too far away from any exon (as defined in the exon file)
    """

    header = ["chrom", "start", "stop", "amplicon", "na", "strand"]

    raw_amplicons = pd.read_csv(bed, sep="\t", names=header)
    raw_amplicons["gene"] = raw_amplicons.apply(lambda x: x["amplicon"].split("_")[0].upper(), axis=1)
    raw_amplicons["altgene"] = raw_amplicons.apply(lambda x: x["amplicon"].split("_")[1].upper(), axis=1)

    def crosscheck_gene_chromosome(row):
        """ Check if this row is close enough to an exon """
        offset = 5000

        gene = row["gene"]
        altgene = row["altgene"]
        chrom = row["chrom"]
        start = row["start"]
        stop = row["stop"]

        gexons = exons[exons.gene == gene]

        if len(gexons) > 0:
            altgexons = exons[exons.gene == altgene]
            if len(gexons) < len(altgexons):
                gexons = altgexons

            chrom_ok = chrom == gexons.chrom.unique()[0]
            start_ok = start > gexons.exon_start.min() - offset
            stop_ok = stop < gexons.exon_start.max() + offset

            return all([chrom_ok, start_ok, stop_ok])
        else:
            #print "gexons empty for", gene
            return False

    raw_amplicons["correct_chrom"] = raw_amplicons.apply(crosscheck_gene_chromosome, axis=1)
    amplicons = raw_amplicons[raw_amplicons.correct_chrom]
    amplicons = amplicons.drop_duplicates(cols=["start", "stop"]) # as in the  'separate strand' beds each exon is present twice

    return amplicons

#####################################################################################################################

def find_bad_positions(coverage_matrix, target_folder = None, trait = None, samplename = None,
                       trait_cutoff = None, whitelist = None):

    """
    Walk through all bases and find contigous regions of bases that fail the coverage/strandbias cutoff
    Create an output file of these regions
    If a base an expected variant (recorded in a ExpectedVariants instance ('whitelist'), record the trait to the instance)
    Return the updated whitelist
    """

    region = None
    last_gene = None
    last_pos = 0
    last_gene = None

    total_bases = 0
    bad_bases = 0

    sample = samplename
    bad_output_name = target_folder + sample + "_failed_regions_%s_cutoff_%s.csv" % (trait, trait_cutoff)
    bad_output = open(bad_output_name, "w")

    if trait == "strandbias":
        output_header = ["gene", "chrom", "start", "stop", "mean_strand_bias", "size"]
    elif trait == "coverage":
        output_header = ["gene", "chrom", "start", "stop", "mean_coverage", "size"]

    bad_output.write("\t".join(output_header) + "\n")

    for index, pandas_dict in coverage_matrix.iterrows():
        chrompos, chrom, dp, gene, minus_dp, plus_dp, pos, start, stop, strand_ratio = pandas_dict
        row = [chrom, start, stop, gene, ".", "NA", (pos-start), strand_ratio]
        row = "\t".join([str(x) for x in row])

        start = int(start)
        stop = int(stop)
        pos = int(pos)

        if whitelist:
            if whitelist.dict.get(chrompos):
                whitelist.add_coverage(chrompos, dp)
                whitelist.add_strand_ratio(chrompos, strand_ratio)

        ##########################################################

        pass_check = False

        if trait == "coverage":
            if dp > trait_cutoff:
                pass_check = True

        elif trait == "strandbias":
            if strand_ratio < trait_cutoff:
                pass_check = True

        total_bases += 1
        if not pass_check:
            bad_bases += 1

        ##########################################################

        if gene != last_gene or last_pos != (pos - 1):
            if region:
                result = region.print_output()
                bad_output.write(result + "\n")
                region = None

        elif not pass_check:
            if region:
                region.add_row(row)
            else:
                region = CoverageCheckClasses.BadRegion(row)
        else:
            if region:
                result = region.print_output()
                bad_output.write(result + "\n")
                region = None

        last_chrom = chrom
        last_pos = pos
        last_gene = gene

    bad_output.close()

    ######################################################################

    good_bases = total_bases - bad_bases
    pass_percent = 100 * round( good_bases / float(total_bases), 3)

    if trait == "strandbias":
        sampleinfo.add_strandbias(sample, [pass_percent, good_bases, total_bases])
        logging.info( "%s percent (%s/%s) of positions in %s have a strand bias below the threshold (%s:1)."
                     % (pass_percent, good_bases, total_bases, sample, trait_cutoff) )

    elif trait == "coverage":
        sampleinfo.add_coverage(sample, [pass_percent, good_bases, total_bases])
        logging.info( "%s percent (%s/%s) of positions in %s have at least the minimum coverage of %sX."
                     % (pass_percent, good_bases, total_bases, sample, trait_cutoff) )

    return whitelist

#####################################################################################################################

def run_bedtools_coverage(bam, output, bed = None):
    """ Run bedtools coverage to get the per-base coverage in the target bam. """

    bedtools_cmd = "bedtools coverage -s -d -abam %s -b %s > %s" % (bam, bed, output)
    logging.debug( bedtools_cmd )
    output_code = subprocess.call(bedtools_cmd, shell=True)

    if output_code != 0:
        logging.critical( "Bedtools coverage run error." )
        logging.critical( bedtools_cmd )
        logging.critical( "Sorry, aborting..." )
        sys.exit()

def run_bedtools_intersect(bed):
    """ Run bedtools intersect to reduce the exon file to the exons with an amplicon overlap. """

    output = bed.replace(".bed", "_covered_exon_locations.bed")

    bedtools_cmd = "intersectBed -u -a %s/input/%s -b %s > %s" % (script_folder, human_exon_file, bed, output)
    logging.debug( bedtools_cmd )
    output_code = subprocess.call(bedtools_cmd, shell=True)

    if output_code != 0:
        logging.critical( "Bedtools intersect run error." )
        logging.critical( bedtools_cmd )
        logging.critical( "Sorry, aborting..." )
        sys.exit()

    return output

def check_bed_filetype(filename):
    """
    Check if the bedfile is actually a bed and convert it to one if its an Illumina Manifest file
    """

    rows = open(filename).readlines()

    if not "Header" in rows[0]:

        if all([len(row.split("\t")) == 6 for row in rows]):
            logging.info( "Correct bed file detected: "+ filename )
            return filename

        else:
            logging.critical( "Error in input bed file:", filename )
            logging.critical( "Not all rows contain 6 fields as is expected." )
            logging.critical( "Sorry, aborting..." )
            sys.exit()

    else: # this is an Illumina manifest file

        logging.error( "Input is not a bed file, but an Illumina manifest file: " + filename )

        expected_bed = filename[:-4] + "_plusminus.bed"

        if os.path.exists(expected_bed):
            logging.info( "Switching to existing bed file: " + expected_bed )
            return expected_bed

        else:
            logging.info( "Converting manifest to bed." )

            import manifest2bed as m2b
            all_out, plusminus = m2b.convert_manifest(filename)

            logging.info( "Switching to newly created bed file: " + plusminus )
            return plusminus

def fix_gene_names_in_bedfile(bed, gene_alias_filename):
    """
    Creates a new bedfile that has all genenames changed according to the gene_aliases file
    Return the name of the new bedfile
    """

    fixed_bed_name = bed.replace(".bed", "_fixed_genenames.bed")

    logging.info( "Gene Name Aliases detected: " + gene_alias_filename )

    if os.path.exists(fixed_bed_name):
        logging.info( "Switching to existing fixed bed file: " + fixed_bed_name )
        return fixed_bed_name

    else:
        logging.info( "Fixing gene names in bedfile." )

        gene_aliases = {}
        with open(gene_alias_filename) as handle:
            for row in handle:
                row = row.strip()
                old, new = re.split("[\t ]", row) # split the row regardsless of tab or space
                gene_aliases[old] = new

        out = open(fixed_bed_name, "w")
        with open(bed) as handle:
            for row in handle:
                f = row.strip().split("\t")
                amplicon = f[3]
                gene = amplicon.split("_")[0]
                new_gene = gene_aliases.get(gene, gene)

                new_amplicon = amplicon.replace(gene, new_gene)

                f[3] = new_amplicon
                result = "\t".join(f)

                out.write(result + "\n")

                if amplicon != new_amplicon:
                    logging.debug("Replaced %s with %s" % (amplicon, new_amplicon))
        out.close()
        return fixed_bed_name

def remove_empty_files_from_folder(folder):
    for filename in os.listdir(folder):
        try:
            if os.path.getsize(folder +"/"+ filename) == 0:
                os.remove(folder +"/"+ filename)
        except:
            pass

#####################################################################################################################
#####################################################################################################################

def run(bed, target_folder, min_dp, max_strand_ratio, whitelist_filename=None, gene_alias_filename=None, target_bams_filename=None):

    if whitelist_filename == "None":
        whitelist_filename = None
    if gene_alias_filename == "None":
        gene_alias_filename = None
    if target_bams_filename == "None":
        target_bams_filename = None

    remove_empty_files_from_folder(target_folder) # remove empty files that might have been left over from previous runs

    if " " in target_folder:
        logging.critical( "The third-party tools used by CoverageCheck don't accept spaces in folder names." )
        logging.critical( "Please replace the spaces in your bam folder with underscores.")
        logging.critical( "Sorry, aborting..." )

    ##############################################################################################
    # configure logging to both sys.stdout and a file

    logging_filename = "CoverageCheck_log.txt"
    print "Log messages printed to %s" % (logging_filename)

    # set up logging to file
    logging.basicConfig(level=logging.INFO,
                  format='%(asctime)s %(name)-8s %(levelname)-8s %(message)s',
                  datefmt='%d-%m-%y %H:%M',
                  filename=logging_filename)

    # set up logging to console
    console = logging.StreamHandler()
    console.setLevel(logging.INFO)
    # set a format which is simpler for console use
    formatter = logging.Formatter('%(name)-12s: %(levelname)-8s %(message)s')
    # add the handler to the root logger
    logging.getLogger('').addHandler(console)

    logging.info( "-" * 100 )
    logging.info( "Minimum accepted coverage per base: %sX" % (min_dp) )
    logging.info( "Maximum accepted coverage ratio between strands: %s:1" % (max_strand_ratio) )

    #####################################################
    # setup variables and report information

    bed = check_bed_filetype(bed)
    if gene_alias_filename:
        bed = fix_gene_names_in_bedfile(bed, gene_alias_filename)

    # parse all exons and all proper amplicons into DataFrames
    exon_filename = run_bedtools_intersect(bed)
    exons, exons_per_gene = parse_exons_into_dataframe_and_dict(exon_filename)

    amplicons = parse_amplicons_into_df(bed, exons)

    # setup the class that collects stats on each sample
    global sampleinfo
    sampleinfo = CoverageCheckClasses.SampleInfo()

    # parse list of expected variants
    if whitelist_filename:
        variant_no = len([x for x in open(whitelist_filename).readlines() if x[0] != "#"])
        logging.info( "Found %s expected variants to check in %s" % (variant_no, whitelist_filename) )

        if variant_no == 0:
            logging.critical( "ERROR: expected variant file doesn't contain variants." )
            sys.exit()
    else:
        logging.info( "No expected variants specified." )

    # parse list of allowed bams
    target_folder = target_folder + "/" # just to be on the safe side
    bams = [x for x in os.listdir(target_folder) if x.endswith(".bam")]

    logging.info( "Found %s bams in the target folder %s" % (len(bams), target_folder))
    if target_bams_filename:
        target_bams = [x.strip() for x in open(target_folder + target_bams_filename).readlines()]
        target_bams = [x+".bam" for x in target_bams if not x.endswith(".bam")]
        logging.info( "Parsed list of %s allowed bams in %s" % (len(target_bams), target_bams_filename) )
        bams = [x for x in bams if x in target_bams]
        logging.info( "%s bams left to process in the target folder." % len(bams) )

    logging.info( "-" * 100 )

    ##############################################################################################
    # process each samples separately

    for bam in bams:

        samplename = bam.split(".")[0]
        sample_output_folder = "%s/%s_results/" % (target_folder, samplename)

        sampleinfo.add_sample(samplename)

        if not os.path.exists(sample_output_folder):
            os.makedirs(sample_output_folder)

        ##############################################################################################
        # create the coverage file (if necessary) and parse it into a pandas DF

        bedtools_output = target_folder + samplename + "_coverage.csv"

        if not os.path.exists(bedtools_output):
            logging.info( "Running bedtools coverage for: " + bam )
            run_bedtools_coverage(target_folder + bam, output = bedtools_output, bed=bed)

        coverage_matrix = parse_coverage_file_into_dataframe(bedtools_output)

        ##############################################################################################
        # Initialize the class instance thats collects information on each expected variant

        whitelist = None
        if whitelist_filename:
            whitelist = CoverageCheckClasses.ExpectedVariants(whitelist_filename, samplename=samplename, folder=sample_output_folder,
                                        dp_cutoff = min_dp, strandbias = max_strand_ratio, exon_dict = exons_per_gene)

        ##############################################################################################
        # run

        if output_undercovered_regions:
            whitelist = find_bad_positions(coverage_matrix, target_folder = sample_output_folder, trait = "coverage",
                               samplename = samplename, trait_cutoff = min_dp, whitelist=whitelist)

        if output_strandbiased_regions:
            whitelist = find_bad_positions(coverage_matrix, target_folder = sample_output_folder, trait = "strandbias",
                               samplename = samplename, trait_cutoff = max_strand_ratio, whitelist=whitelist)

        if output_undercovered_regions or output_strandbiased_regions:
            if whitelist_filename:
                whitelist.print_output()

        ##############################################################################################
        # plot

        if output_coverage_plots:
            sample_df = plotting.create_all_coverage_plots(bedtools_output, exons=exons, exons_per_gene = exons_per_gene,
                                        target_folder = sample_output_folder, whitelist=whitelist, amplicons=amplicons)
        logging.info(  "-" * 100 )

    ##############################################################################################
    # create summaries across all samples
    # this has to happen after the individual bam treatment because it assumes that coverage files for each bam are present

    sampleinfo.print_output()

    all_sample_filename = "all_samples.csv"
    UniteCoverage.unite(all_sample_filename, target_folder = target_folder)

    byte_size = os.path.getsize(target_folder+all_sample_filename)
    mb_size = byte_size/1.049e+6

    if output_all_sample_summary:
        if mb_size > 2000:
            logging.error( "The united coverage size of all samples is %sMB. Skipping summary plots." % (round(mb_size,2)) )
        else:
            all_sample_plotting.plot_exon_coverage(all_sample_filename, exons=exons, exons_per_gene = exons_per_gene,
                                                   target_folder = target_folder, whitelist=whitelist)

##################################################################################

if __name__ == "__main__":

    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument("-b", "--bed", help="region file in bed or Illumina Manifest format",
                        required=True)

    parser.add_argument("-x", "--expected_variants", help="expected variants in vcf or HGMD BioMart export format")
    parser.add_argument("-a", "--alias", help="gene name aliases, format: [old new]")
    parser.add_argument("-w", "--whitelist", help="list of bams to analyse")

    parser.add_argument("-c", "--min_coverage",  type=int, help="minimum coverage (Default: 50X)", default=default_min_coverage)
    parser.add_argument("-s", "--max_strandratio", type=float, help="maximum strand ratio (Default: 5)", default=default_max_strandbias)

    args = parser.parse_args()
    pprint.pprint(args)

    ##################################################################################

    bed = args.bed
    whitelist_filename = args.expected_variants
    gene_alias_filename = args.alias
    target_bams_filename = args.whitelist
    min_dp = args.min_coverage
    max_strand_ratio = args.max_strandratio

    if "/" in bed:
        target_folder = "/".join(bed.split("/")[:-1]) + "/"
    else:
        target_folder = "./"

    run(bed,
        target_folder,
        min_dp,
        max_strand_ratio,
        whitelist_filename=whitelist_filename,
        gene_alias_filename=gene_alias_filename,
        target_bams_filename=target_bams_filename)

max_stars_repo_name: aweller/CoverageCheck

max_issues_repo_name: aweller/CoverageCheck

max_forks_repo_name: aweller/CoverageCheck


text:
import os
import time
from datetime import datetime
import picamera
from picamera import PiCamera
import cv2
import numpy as np

def contrast_stretch(im):
    in_min = np.percentile(im, 5)
    in_max = np.percentile(im, 95)
    out_min = 0.0
    out_max = 255

    out = im - in_min
    out *= ((out_min - out_max) / (in_min - in_max))
    out += in_min

    return out

# Settings
countdown = 1  # Interval between image captures (seconds)
cam_width = 1280  # Cam sensor width
cam_height = 480  # Cam sensor height

scale_ratio = 0.5

# Camera resolution height must be dividable by 16, and width by 32
cam_width = int((cam_width + 31) / 32) * 32
cam_height = int((cam_height + 15) / 16) * 16
print("Used camera resolution: " + str(cam_width) + " x " + str(cam_height))

# Buffer for captured image settings
img_width = int(cam_width * scale_ratio)
img_height = int(cam_height * scale_ratio)
capture = np.zeros((img_height, img_width, 4), dtype=np.uint8)
print("Scaled image resolution: " + str(img_width) + " x " + str(img_height))

# Initialize camera
camera = PiCamera(stereo_mode='side-by-side', stereo_decimate=False)
camera.resolution = (cam_width, cam_height)
camera.framerate = 20
camera.hflip = False

counter = 0
t2 = datetime.now()
print("Starting photo sequence")
for frame in camera.capture_continuous(capture, format="bgra", use_video_port=True, resize=(img_width, img_height)):
    t1 = datetime.now()
    cntdwn_timer = countdown - int((t1 - t2).total_seconds())

    # If cowntdown is zero -> record next image
    if cntdwn_timer == -1:
        counter += 1

        # Save raw image
        cv2.imwrite('raw' + str(counter) + '.png', frame)

        # Split image to left and right
        imgLeft = frame[0:240, 0:320]
        imgRight = frame[0:240, 320:640]
        leftName = str(counter) + 'nir.png'
        rightName = str(counter) + 'r.png')
        cv2.imwrite(leftName, imgLeft)
        cv2.imwrite(rightName, imgRight)

        # Split left & right into components
        b, g, nir, a = cv2.split(imgLeft)
        b1, g1, r, a1 = cv2.split(imgRight)

        # Calculate NDVI
        bottom = (nir.astype(float) + r.astype(float))
        bottom[bottom == 0] = 0.01
        ndvi = (nir.astype(float) - r) / bottom
        ndvi = contrast_stretch(ndvi)
        ndvi = ndvi.astype(np.uint8)
        ndvi = cv2.applyColorMap(ndvi, cv2.COLORMAP_JET)
        cv2.imwrite('ndvi' + str(counter) + '.png', ndvi)

        # Calculate Chlorophyll Index
        b1 = g1.astype(float)
        b1[b1 == 0] = 0.01
        ci = nir.astype(float) / b1
        ci = contrast_stretch(ci)
        ci = ci.astype(np.uint8)
        ci = cv2.applyColorMap(ci, cv2.COLORMAP_SUMMER)
        cv2.imwrite('ci' + str(counter) + '.png', ci)

        t2 = datetime.now()
        time.sleep(1)
        cntdwn_timer = 0
        next

        key = cv2.waitKey(1)

        # Press q to quit
        if (key == ord("q")):
            break

    print("Photo sequence finished")

ext: py

max_stars_repo_name: etnsch/Senior-Project

max_issues_repo_name: etnsch/Senior-Project

max_forks_repo_name: etnsch/Senior-Project

reason: import numpy


text:
''' Assignment 3 - More Pandas
    This assignment requires more individual learning then the last one did - you are encouraged to check out the [pandas documentation](http://pandas.pydata.org/pandas-docs/stable/) to find functions or methods you might not have used yet, or ask questions on [Stack Overflow](http://stackoverflow.com/) and tag them as pandas and python related. And of course, the discussion forums are open for interaction with your peers and the course staff.

    Question 1 (20%)
    Load the energy data from the file `Energy Indicators.xls`, which is a list of indicators of [energy supply and renewable electricity production](Energy%20Indicators.xls) from the [United Nations](http://unstats.un.org/unsd/environment/excel_file_tables/2013/Energy%20Indicators.xls) for the year 2013, and should be put into a DataFrame with the variable name of **energy**.

    Keep in mind that this is an Excel file, and not a comma separated values file. Also, make sure to exclude the footer and header information from the datafile. The first two columns are unneccessary, so you should get rid of them, and you should change the column labels so that the columns are:

    `['Country', 'Energy Supply', 'Energy Supply per Capita', '% Renewable']`

    Convert `Energy Supply` to gigajoules (there are 1,000,000 gigajoules in a petajoule). For all countries which have missing data (e.g. data with "...") make sure this is reflected as `np.NaN` values.

    Rename the following list of countries (for use in later questions):

    "Republic of Korea": "South Korea",
    "United States of America": "United States",
    "United Kingdom of Great Britain and Northern Ireland": "United Kingdom",
    "China, Hong Kong Special Administrative Region": "Hong Kong"```

    There are also several countries with numbers and/or parenthesis in their name. Be sure to remove these,

    e.g.

    'Bolivia (Plurinational State of)'` should be `'Bolivia'`,
    `'Switzerland17'` should be `'Switzerland'`.

    Next, load the GDP data from the file `world_bank.csv`, which is a csv containing countries' GDP from 1960 to 2015 from [World Bank](http://data.worldbank.org/indicator/NY.GDP.MKTP.CD). Call this DataFrame **GDP**.

    Make sure to skip the header, and rename the following list of countries:

    "Korea, Rep.": "South Korea",
    "Iran, Islamic Rep.": "Iran",
    "Hong Kong SAR, China": "Hong Kong"```

    Finally, load the [Sciamgo Journal and Country Rank data for Energy Engineering and Power Technology](http://www.scimagojr.com/countryrank.php?category=2102) from the file `scimagojr-3.xlsx`, which ranks countries based on their journal contributions in the aforementioned area. Call this DataFrame **ScimEn**.

    Join the three datasets: GDP, Energy, and ScimEn into a new dataset (using the intersection of country names). Use only the last 10 years (2006-2015) of GDP data and only the top 15 countries by Scimagojr 'Rank' (Rank 1 through 15).

    The index of this DataFrame should be the name of the country, and the columns should be ['Rank', 'Documents', 'Citable documents', 'Citations', 'Self-citations',
    'Citations per document', 'H index', 'Energy Supply',
    'Energy Supply per Capita', '% Renewable', '2006', '2007', '2008',
    '2009', '2010', '2011', '2012', '2013', '2014', '2015'].

    *This function should return a DataFrame with 20 columns and 15 entries.*
'''
import pandas as pd
import numpy as np

energy = pd.read_excel('Energy Indicators.xls', skiprows=17)

# remove header and footer
energy = energy[:227]

# drop useless columns
energy = energy.drop(energy.columns[[0, 1]], axis=1)

# rename columns
energy.columns = ['Country', 'Energy Supply',
                  'Energy Supply per Capita', '% Renewable']

# replace ... with NaN
energy.replace('...', np.nan, inplace=True)

# convert Energy Supply units
energy['Energy Supply'] = energy['Energy Supply'] * 1000000

# remove all numbers from country names and parenthesis

def remove_numbers(country):
    return ''.join(filter(lambda x: not x.isdigit(), country))

def remove_parenthesis(country):
    start = country.find('(')
    if start > -1:
        return country[:start-1]
    return country

energy['Country'] = energy['Country'].apply(remove_numbers)
energy['Country'] = energy['Country'].apply(remove_parenthesis)

# replace country names
country_replace = {'Republic of Korea': 'South Korea',
                   'United States of America': 'United States',
                   'United Kingdom of Great Britain and Northern Ireland': 'United Kingdom',
                   'China, Hong Kong Special Administrative Region': 'Hong Kong',
                   }
# energy.replace(country_replace, inplace = True)
energy.replace({'Country': country_replace}, inplace=True)

GDP = pd.read_csv('world_bank.csv', skiprows=4)

# rename Country Name to Country so that merge will work later on.
GDP.rename(columns={'Country Name': 'Country'}, inplace=True)

country_replace = {'Korea, Rep.': 'South Korea',
                   'Iran, Islamic Rep.': 'Iran',
                   'Hong Kong SAR, China': 'Hong Kong'
                   }

GDP.replace({'Country': country_replace}, inplace=True)

ScimEn = pd.read_excel('scimagojr-3.xlsx')

def answer_one():
    global energy, GDP, ScimEn

    # merge all three dataframes
    combined_df = pd.merge(
        pd.merge(energy, GDP, on='Country'), ScimEn, on='Country')

    # set index to country name
    combined_df.set_index('Country', inplace=True)

    # keep specific columns
    combined_df = combined_df[['Rank', 'Documents', 'Citable documents', 'Citations',
                               'Self-citations', 'Citations per document', 'H index',
                               'Energy Supply', 'Energy Supply per Capita', '% Renewable',
                               '2006', '2007', '2008', '2009', '2010', '2011', '2012', '2013',
                               '2014', '2015']]

    # only include Ranking 15 or higher
    combined_df = combined_df[combined_df['Rank'] <= 15]

    # organize the Ranking
    combined_df = combined_df.sort_values(['Rank'], ascending=[1])

    return combined_df

''' Question 2 (6.6%)
    The previous question joined three datasets then reduced this to just the top 15 entries. When you joined the datasets, but before you reduced this to the top 15 items, how many entries did you lose?

    *This function should return a single number.*
'''

def answer_two():
    global energy, GDP, ScimEn

    un = pd.merge(pd.merge(energy, GDP, on='Country', how='outer'),
                  ScimEn, on='Country', how='outer')
    int = pd.merge(pd.merge(energy, GDP, on='Country'), ScimEn, on='Country')

    return (len(un) - len(int))

''' Answer the following questions in the context of only the top 15 countries by Scimagojr Rank (aka     the DataFrame returned by `answer_one()`)
'''

''' Question 3 (6.6%)
    What is the average GDP over the last 10 years for each country? (exclude missing values from this calculation.)

    *This function should return a Series named `avgGDP` with 15 countries and their average GDP sorted in descending order.*
'''

def answer_three():
    Top15 = answer_one()
    avgGDP = Top15[['2006', '2007', '2008', '2009', '2010', '2011', '2012',
                    '2013', '2014', '2015']].mean(axis=1).sort_values(ascending=False)

    return avgGDP

''' Question 4 (6.6%)
    By how much had the GDP changed over the 10 year span for the country with the 6th largest average GDP?

    *This function should return a single number.*
'''

def answer_four():
    Top15 = answer_one()
    avgGDP = answer_three().reset_index()

    return np.float64(Top15.filter(like=avgGDP.iloc[5]['Country'], axis=0)['2015'] - Top15.filter(like=avgGDP.iloc[5]['Country'], axis=0)['2006'])

''' Question 5 (6.6%)
    What is the mean `Energy Supply per Capita`?

    *This function should return a single number.*
'''

def answer_five():
    Top15 = answer_one()

    return Top15['Energy Supply per Capita'].mean()

''' Question 6 (6.6%)
    What country has the maximum % Renewable and what is the percentage?

    *This function should return a tuple with the name of the country and the percentage.*
'''

def answer_six():
    Top15 = answer_one()
    country = Top15.sort_values(by='% Renewable', ascending=False).iloc[0]

    return (country.name, country['% Renewable'])

''' Question 7 (6.6%)
    Create a new column that is the ratio of Self-Citations to Total Citations.
    What is the maximum value for this new column, and what country has the highest ratio?

    *This function should return a tuple with the name of the country and the ratio.*
'''

def answer_seven():
    Top15 = answer_one()
    Top15['Citation ratio'] = (Top15['Self-citations'] / Top15['Citations'])
    country = Top15.sort_values(by='Citation ratio', ascending=False).iloc[0]

    return (country.name, country['Citation ratio'])

''' Question 8 (6.6%)

    Create a column that estimates the population using Energy Supply and Energy Supply per capita.
    What is the third most populous country according to this estimate?

    *This function should return a single string value.*
'''

def answer_eight():
    Top15 = answer_one()
    Top15['Population'] = (Top15['Energy Supply'] /
                           Top15['Energy Supply per Capita'])

    return Top15.sort_values(by='Population', ascending=False).iloc[2].name

''' Question 9 (6.6%)
    Create a column that estimates the number of citable documents per person.
    What is the correlation between the number of citable documents per capita and the energy supply per capita? Use the `.corr()` method, (Pearson's correlation).

    *This function should return a single number.*
'''

def answer_nine():
    Top15 = answer_one()
    Top15['PopEst'] = Top15['Energy Supply'] / \
        Top15['Energy Supply per Capita']
    Top15['Citable docs per Capita'] = Top15['Citable documents'] / Top15['PopEst']

    return Top15[['Energy Supply per Capita', 'Citable docs per Capita']].corr().loc['Energy Supply per Capita', 'Citable docs per Capita']

''' Question 10 (6.6%)
    Create a new column with a 1 if the country's % Renewable value is at or above the median for all countries in the top 15, and a 0 if the country's % Renewable value is below the median.

    *This function should return a series named `HighRenew` whose index is the country name sorted in ascending order of rank.*
'''

def answer_ten():
    Top15 = answer_one()

    # T / F for % renewable over median or not
    Top15['HighRenew'] = Top15['% Renewable'] >= Top15['% Renewable'].median()
    Top15['HighRenew'] = Top15['HighRenew'].apply(lambda x: 1 if x else 0)

    # sorted by Rank
    Top15.sort_values(by='Rank', inplace=True)

    return Top15['HighRenew']

''' Question 11 (6.6%)
    Use the following dictionary to group the Countries by Continent, then create a dateframe that displays the sample size (the number of countries in each continent bin), and the sum, mean, and std deviation for the estimated population of each country.

    ContinentDict  = {'China':'Asia',
                  'United States':'North America',
                  'Japan':'Asia',
                  'United Kingdom':'Europe',
                  'Russian Federation':'Europe',
                  'Canada':'North America',
                  'Germany':'Europe',
                  'India':'Asia',
                  'France':'Europe',
                  'South Korea':'Asia',
                  'Italy':'Europe',
                  'Spain':'Europe',
                  'Iran':'Asia',
                  'Australia':'Australia',
                  'Brazil':'South America'}

*This function should return a DataFrame with index named Continent `['Asia', 'Australia', 'Europe', 'North America', 'South America']` and columns `['size', 'sum', 'mean', 'std']`*
'''

def answer_eleven():
    Top15 = answer_one()
    ContinentDict = {'China': 'Asia',
                     'United States': 'North America',
                     'Japan': 'Asia',
                     'United Kingdom': 'Europe',
                     'Russian Federation': 'Europe',
                     'Canada': 'North America',
                     'Germany': 'Europe',
                     'India': 'Asia',
                     'France': 'Europe',
                     'South Korea': 'Asia',
                     'Italy': 'Europe',
                     'Spain': 'Europe',
                     'Iran': 'Asia',
                     'Australia': 'Australia',
                     'Brazil': 'South America'
                     }

    groups = pd.DataFrame(columns=['size', 'sum', 'mean', 'std'])
    Top15['PopEst'] = Top15['Energy Supply'] / \
        Top15['Energy Supply per Capita']

    for continent, count in Top15.groupby(ContinentDict):
        #         print(count['PopEst'])
        #         print('*'*40)
        groups.loc[continent] = [len(count), count['PopEst'].sum(
        ), count['PopEst'].mean(), count['PopEst'].std()]

    return groups

''' Question 12 (6.6%)
    Cut % Renewable into 5 bins. Group Top15 by the Continent, as well as these new % Renewable bins. How many countries are in each of these groups?

    *This function should return a __Series__ with a MultiIndex of `Continent`, then the bins for `% Renewable`. Do not include groups with no countries.*
'''

def answer_twelve():
    Top15 = answer_one()
    ContinentDict = {'China': 'Asia',
                     'United States': 'North America',
                     'Japan': 'Asia',
                     'United Kingdom': 'Europe',
                     'Russian Federation': 'Europe',
                     'Canada': 'North America',
                     'Germany': 'Europe',
                     'India': 'Asia',
                     'France': 'Europe',
                     'South Korea': 'Asia',
                     'Italy': 'Europe',
                     'Spain': 'Europe',
                     'Iran': 'Asia',
                     'Australia': 'Australia',
                     'Brazil': 'South America'
                     }

    Top15['Continent'] = Top15.index.map(lambda c: ContinentDict[c])

    # https://pandas.pydata.org/pandas-docs/stable/generated/pandas.cut.html
    Top15['% Renewable bins'] = pd.cut(Top15['% Renewable'], 5)

    return Top15.groupby(['Continent', '% Renewable bins']).size()

''' Question 13 (6.6%)
    Convert the Population Estimate series to a string with thousands separator (using commas). Do not round the results.

    e.g. 317615384.61538464 -> 317,615,384.61538464

    *This function should return a Series `PopEst` whose index is the country name and whose values are the population estimate string.*
'''

def answer_thirteen():
    Top15 = answer_one()
    Top15['PopEst'] = Top15['Energy Supply'] / \
        Top15['Energy Supply per Capita']

    # https://stackoverflow.com/questions/22617/format-numbers-to-strings-in-python
    return Top15['PopEst'].apply(lambda str: '{0:,}'.format(str))

print('-'*40)
print(answer_one())
print('-'*40)
print(answer_two())
print('-'*40)
print(answer_three())
print('-'*40)
print(answer_four())
print('-'*40)
print(answer_five())
print('-'*40)
print(answer_six())
print('-'*40)
print(answer_seven())
print('-'*40)
print(answer_eight())
print('-'*40)
print(answer_nine())
print('-'*40)
print(answer_ten())
print('-'*40)
print(answer_eleven())
print('-'*40)
print(answer_twelve())
print('-'*40)
print(answer_thirteen())


text:
"""Utility functions for differential equation problems."""

from typing import Optional, Sequence, Union

import numpy as np

from probnum import problems, randprocs
from probnum.diffeq.odefilter import approx_strategies, information_operators
from probnum.typing import FloatArgType

__all__ = ["ivp_to_regression_problem"]

# The ODE information operator is not optional, because in order to create it
# one needs to know the order of the algorithm that is desired (i.e. num_prior_derivatives).
# Since this is a weird input for the function, it seems safer to just require
# the full operator.
def ivp_to_regression_problem(
    ivp: problems.InitialValueProblem,
    locations: Union[Sequence, np.ndarray],
    ode_information_operator: information_operators.InformationOperator,
    approx_strategy: Optional[approx_strategies.ApproximationStrategy] = None,
    ode_measurement_variance: Optional[FloatArgType] = 0.0,
    exclude_initial_condition=False,
):
    """Transform an initial value problem into a regression problem.

    Parameters
    ----------
    ivp
        Initial value problem to be transformed.
    locations
        Locations of the time-grid-points.
    ode_information_operator
        ODE information operator to use.
    approx_strategy
        Approximation strategy to use. Optional. Default is `EK1()`.
    ode_measurement_variance
        Artificial ODE measurement noise. Optional. Default is 0.0.
    exclude_initial_condition
        Whether to exclude the initial condition from the regression problem.
        Optional. Default is False, in which case the returned measurement model list
        consist of [`initcond_mm`, `ode_mm`, ..., `ode_mm`].

    Returns
    -------
    problems.TimeSeriesRegressionProblem
        Time-series regression problem.
    """

    # Construct data and solution
    N = len(locations)
    data = np.zeros((N, ivp.dimension))
    if ivp.solution is not None:
        solution = np.stack([ivp.solution(t) for t in locations])
    else:
        solution = None

    ode_information_operator.incorporate_ode(ivp)

    # Construct measurement models
    measmod_initial_condition, measmod_ode = _construct_measurement_models(
        ivp, ode_information_operator, approx_strategy, ode_measurement_variance
    )

    if exclude_initial_condition:
        measmod_list = [measmod_ode] * N
    else:
        measmod_list = [measmod_initial_condition] + [measmod_ode] * (N - 1)

    # Return regression problem
    return problems.TimeSeriesRegressionProblem(
        locations=locations,
        observations=data,
        measurement_models=measmod_list,
        solution=solution,
    )

def _construct_measurement_models(
    ivp, ode_information_operator, approx_strategy, ode_measurement_variance
):
    """Construct measurement models for the IVP."""

    transition_matrix = np.eye(
        ode_information_operator.output_dim, ode_information_operator.input_dim
    )
    shift_vector = -ivp.y0

    if ode_measurement_variance > 0.0:
        measmod_y0, measmod_ode = _construct_measurement_models_gaussian_likelihood(
            ode_information_operator,
            shift_vector,
            transition_matrix,
            approx_strategy,
            ode_measurement_variance,
        )
    else:
        measmod_y0, measmod_ode = _construct_measurement_models_dirac_likelihood(
            ode_information_operator,
            shift_vector,
            transition_matrix,
            approx_strategy,
        )
    return measmod_y0, measmod_ode

def _construct_measurement_models_gaussian_likelihood(
    ode_information_operator,
    shift_vector,
    transition_matrix,
    approx_strategy,
    ode_measurement_variance,
):
    """Construct measurement models for the IVP with Gaussian likelihoods."""

    def diff(t):
        return ode_measurement_variance * np.eye(ode_information_operator.output_dim)

    def diff_cholesky(t):
        return np.sqrt(ode_measurement_variance) * np.eye(
            ode_information_operator.output_dim
        )

    measmod_initial_condition = randprocs.markov.discrete.LTIGaussian(
        state_trans_mat=transition_matrix,
        shift_vec=shift_vector,
        proc_noise_cov_mat=diff(None),
        proc_noise_cov_cholesky=diff_cholesky(None),
    )
    if approx_strategy is not None:
        ode_information_operator = approx_strategy(ode_information_operator)
    measmod_ode = ode_information_operator.as_transition(
        measurement_cov_fun=diff, measurement_cov_cholesky_fun=diff_cholesky
    )

    return measmod_initial_condition, measmod_ode

def _construct_measurement_models_dirac_likelihood(
    ode_information_operator, shift_vector, transition_matrix, approx_strategy
):
    """Construct measurement models for the IVP with Dirac likelihoods."""
    measmod_initial_condition = randprocs.markov.discrete.LTIGaussian.from_linop(
        state_trans_mat=transition_matrix, shift_vec=shift_vector
    )
    if approx_strategy is not None:
        ode_information_operator = approx_strategy(ode_information_operator)
    measmod_ode = ode_information_operator.as_transition()

    return measmod_initial_condition, measmod_ode

max_issues_repo_name: simeoncarstens/probnum

max_forks_repo_name: JonathanWenger/probnum


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Sep 20 10:46:25 2021

@author: gregz
"""

from astropy.io import fits
from astropy.modeling.models import Gaussian2D
from astropy.modeling.fitting import LevMarLSQFitter
import numpy as np
import glob
import os.path as op
from input_utils import setup_logging

log = setup_logging('adr')
basedir = '/work/03946/hetdex/maverick/LRS2/STANDARDS'

filenames = glob.glob(op.join(basedir, 'multi_2021*uv.fits'))

N = 11
M = len(filenames)
k=0
XC, YC, WC = (np.zeros((M, 2*N)), np.zeros((M, 2*N)), np.zeros((2*N,)))

fitter = LevMarLSQFitter()

G = Gaussian2D()

for filename in filenames:
    log.info('Working on %s' % filename)
    j = 0
    for name in ['uv', 'orange']:
        f = fits.open(filename.replace('uv', name))
        wave = f[6].data[0]
        x = f[5].data[:, 0]
        y = f[5].data[:, 1]
        skysub = f[2].data
        chunks = [np.nanmedian(xo, axis=1) for xo in np.array_split(skysub, N,
                                                                    axis=1)]
        wc = [np.nanmedian(xo) for xo in np.array_split(wave, N)]
        xc = np.zeros((N,))
        yc = np.zeros((N,))
        wc = np.array(wc)
        for i, chunk in enumerate(chunks):
            ind = np.nanargmax(chunk)
            d = np.sqrt((x-x[ind])**2 + (y-y[ind])**2)
            xc[i] = x[ind]
            yc[i] = y[ind]
            G.amplitude.value = np.nanmax(chunk)
            G.x_mean.value = xc[i]
            G.y_mean.value = yc[i]
            sel = (d < 3.0) * np.isfinite(chunk)
            fit = fitter(G, x[sel], y[sel], chunk[sel])
            xc[i] = fit.x_mean.value * 1.
            yc[i] = fit.y_mean.value * 1.
        XC[k, j:j+N] = xc
        YC[k, j:j+N] = yc
        WC[j:j+N] = wc
        j += N

    XC[k, WC<4650.] += 0.20
    YC[k, WC<4650.] -= 0.24
    sel = np.isfinite(XC[k]) * np.isfinite(YC[k])
    p0 = np.polyfit(WC[sel], XC[k][sel], 3)
    p1 = np.polyfit(WC[sel], YC[k][sel], 3)
    xc = np.polyval(p0, 5500)
    yc = np.polyval(p1, 5500)
    XC[k] = XC[k] - xc
    YC[k] = YC[k] - yc
    if np.sqrt(xc**2 + yc**2)>2.5:
        XC[k] = np.nan
        YC[k] = np.nan
    log.info('Centroid: %0.2f %0.2f' % (xc, yc))
    k += 1

fits.HDUList([fits.PrimaryHDU(XC), fits.ImageHDU(YC), fits.ImageHDU(WC)]).writeto('test.fits', overwrite=True)

max_stars_repo_name: grzeimann/Panacea

max_issues_repo_name: grzeimann/Panacea

max_forks_repo_name: grzeimann/Panacea

reason: import numpy,from astropy


text:
import matplotlib.pyplot as plt
import numpy as np

plot_name = "botbowl-1"
exp_id_nopf = '7f556c38-7853-11eb-b93b-faffc23fefdb.dat'
exp_id_pf = '129055f2-7846-11eb-a07b-faffc23fefdb.dat'
window_size = 2

plot_name = "botbowl-3"
exp_id_nopf = '48ed35ea-79b7-11eb-ba58-19faed4b7487.dat'
exp_id_pf = '4cd8c832-794d-11eb-b2ed-19faed4b7487.dat'
window_size = 20

#plot_name = "botbowl-5"
#exp_id_nopf = 'f1e2354e-79d8-11eb-bb13-19faed4b7487.dat'
#exp_id_pf = '70310ce8-7a80-11eb-8061-19faed4b7487.dat'
#window_size = 100

logfile_nopf = f"logs/logs-pf/{plot_name}/nopf/{exp_id_nopf}"
logfile_pf = f"logs/logs-pf/{plot_name}/pf/{exp_id_pf}"
fig, axs = plt.subplots(1, 3, figsize=(4 * 3, 5))

axs[0].ticklabel_format(axis="x", style="sci", scilimits=(0, 0))
axs[0].set_title('Reward')
axs[1].ticklabel_format(axis="x", style="sci", scilimits=(0, 0))
axs[1].set_title('TD/Episode')
#axs[1].set_ylim(bottom=0.0)
#axs[1].set_xlim(left=0)
axs[2].ticklabel_format(axis="x", style="sci", scilimits=(0, 0))
axs[2].set_title('Win rate')
#axs[2].set_yticks(np.arange(0, 1.001, step=0.1))
#axs[2].set_xlim(left=0)

def plot(log_steps, log_mean_reward, log_td_rate, log_win_rate, label):
    axs[0].plot(log_steps, log_mean_reward, label=label)
    axs[1].plot(log_steps, log_td_rate, label=label)
    axs[2].plot(log_steps, log_win_rate, label=label)

def load_data(logfile, window_size):
    data = []
    with open(logfile, 'r') as file:
        for line in file:
            data.append([float(v) for v in line.strip().split(", ")])
    cols = []
    for i in range(len(data[0])):
        col = []
        vals = []
        for v in np.array(data)[:, i]:
            vals.append(v)
            if len(vals) == window_size:
                col.append(np.average(vals))
                vals.clear()
        cols.append(col)
    return np.array(cols).transpose()

data = load_data(logfile_nopf, window_size=window_size)
plot(log_steps=data[:, 2], log_mean_reward=data[:, 6], log_td_rate=data[:, 4], log_win_rate=data[:, 3], label="No Pathfinding")
data = load_data(logfile_pf, window_size=window_size)
plot(log_steps=data[:, 2], log_mean_reward=data[:, 6], log_td_rate=data[:, 4], log_win_rate=data[:, 3], label="Pathfinding")

# axs[0].legend(['Baseline', 'Pathfinding'])
# axs[1].legend(['Baseline', 'Pathfinding'])
axs[2].legend(['Baseline', 'Pathfinding'])

fig.tight_layout()
plotfile = f"plots/{plot_name}.png"
fig.savefig(plotfile)
plt.close('all')

max_stars_repo_name: hadarshavit/botbowl

max_issues_repo_name: ernestvmo/botbowl

max_forks_repo_name: ernestvmo/botbowl


text:
import _init_pathes
import logging
import datetime
import os.path as osp
import numpy as np
from configs.lung_cfgs import cfg as lung_cfg
from .lung_segmentor import LungSegmentor
from configs.cfgs import cfg as risk_cfg
from .risk_predictor import RiskPredictor
from data_process.data_processor import get_series_uids
from data_process.data_processor import load_series_volume
from data_process.data_processor import interpolate_volume
from data_process.data_processor import get_series_dcm_nums
from data_process.data_processor import get_convex_bbox_frm_3dmask
from data_process.data_processor import resize, hu2gray
import matplotlib.pyplot as plt

class CAD_Pipeline(object):
    def __init__(self):
        logging.basicConfig(level=logging.DEBUG,
                            format='%(asctime)s %(filename)s[line:%(lineno)d] %(levelname)s %(message)s',
                            datefmt='%a, %d %b %Y %H:%M:%S',
                            filename='nodule_detector.log',
                            filemode='w')
        console = logging.StreamHandler()
        console.setLevel(logging.INFO)
        formatter = logging.Formatter('%(name)-12s: %(levelname)-8s %(message)s')
        console.setFormatter(formatter)
        logging.getLogger('').addHandler(console)

    def initialize(self, debug_mode = True):
        '''
        Before pulmonary lesion detection, you have to call this function to initialize the networks.
        '''
        self.debug_mode = debug_mode
        self.is_ready = False
        #init lung segementation network
        self.lung_segmentor = LungSegmentor(lung_cfg)
        if not self.lung_segmentor.initialize():
            self._record_error_info(self.lung_segmentor.last_error)
            self.is_ready = False
            return False

        self.risk_predictor = RiskPredictor(risk_cfg)
        if not self.risk_predictor.initialize():
            self._record_error_info(self.risk_predictor.last_error)
            self.is_ready = False
            return False

        if osp.exists('../checkpoints/feature_minv.npy'):
            print('found minv file for normalization')
            self.patient_infominv = np.load('../checkpoints/feature_minv.npy', allow_pickle=True)

        else:
            self.patient_infominv = np.min(self.patients[:, 1:48], axis=0)
        self.patient_infominv[1] = 1

        if osp.exists('../checkpoints/feature_maxv.npy'):
            print('found maxv file for normalization')
            self.patient_infomaxv = np.load('../checkpoints/feature_maxv.npy', allow_pickle=True)

        else:
            self.patient_infomaxv = np.max(self.patients[:, 1:48], axis=0)
        self.patient_infomaxv[1] = 100

        self.patient_infominv = np.delete(self.patient_infominv, 2, axis=0)
        self.patient_infomaxv = np.delete(self.patient_infomaxv, 2, axis=0)
        self.is_ready = True
        self._record_debug_info('Networks are ready!')
        return self.is_ready

    def predict_risk_for(self, patient_root):
        '''
        Detection on scan of dcm files
        '''
        assert self.is_ready==True, 'Networks not ready!'
        if not osp.exists(patient_root):
            self._record_error_info('No patient path found at ' + patient_root)
            return False
        treatment_file = osp.join(patient_root, 'treatments.npy')
        clinical_file = osp.join(patient_root, 'clinical_features.npy')

        if not osp.exists(treatment_file) or not osp.exists(clinical_file):
            self._record_error_info('treatments.npy and clinical_features.npy required, but no file was found at ' + patient_root)
            return False
        try:
            treatments = np.load(treatment_file, allow_pickle=True).astype(np.float32)
            treatments = treatments.reshape((1, -1))
            clinicals = np.load(clinical_file, allow_pickle=True)
            int_clinicals = np.int32(clinicals.copy())
            keep = np.where(int_clinicals>0)

            clinicals = (clinicals-self.patient_infominv)/(self.patient_infomaxv-self.patient_infominv)
            clinicals[clinicals>1.0] = 1.0
            clinicals[clinicals<0.0] = 0.0
            clinicals = clinicals.reshape(1, -1)
            #processing ct scan if it exists
            feed_volume = self._try_to_load_ct(osp.join(patient_root, 'ct'))
            if feed_volume is None:
                feed_volume = np.zeros([1]+risk_cfg.im_feedsize+ [1], dtype=np.float32)

            feature_significane_coff, reg_pred = self.risk_predictor.predict(treatments,feed_volume, clinicals)
            feature_significane_coff = feature_significane_coff[0,...]

            reg_pred = reg_pred[0,...]
            if feature_significane_coff is None:
                return False

            np.save(osp.join(patient_root, 'predictions.npy'), reg_pred)

            self._generate_figures(patient_root, keep, feature_significane_coff, reg_pred)
            return True
        except Exception as e:
            self._record_error_info(str(e))

            return False

    def _try_to_load_ct(self, ct_root):
        if not osp.exists(ct_root):
            return None
        try:
            series_uids = get_series_uids(ct_root)
            if series_uids is None or len(series_uids)==0:
                self._record_error_info('No CT scan found at ' + ct_root)
                return None

            feed_volume = None
            for series_uid in series_uids:
                n_dcmfiles = get_series_dcm_nums(ct_root, series_uid)
                if n_dcmfiles < 60:
                    self._record_debug_info('CT scan ignored due to the number of dcm files<50 of series:' + series_uid)
                    continue

                self._record_debug_info('Start loading series:' + series_uid)

                volume, origin, spacing, slice_spacing, im_numbers = load_series_volume(ct_root, series_uid)
                if volume is None:
                    self._record_error_info('Load volume data failed.')
                    return None
                new_spacing = [0.65, 0.65, 1.5]
                new_volume = interpolate_volume(volume.copy(), spacing, new_spacing)
                nD,nH,nW = new_volume.shape

                self._record_debug_info('Lung segmentation ...')
                lung_mask_vol = self.lung_segmentor.segment_lung(new_volume)
                if lung_mask_vol is None or np.count_nonzero(lung_mask_vol)<(nD*nH*nW*0.01):
                    return None

                bbox = get_convex_bbox_frm_3dmask(lung_mask_vol, [2,2,2])

                feed_volume = resize(hu2gray(new_volume[bbox[2]:bbox[5], bbox[1]:bbox[4], bbox[0]:bbox[3]],WL=-500, WW=1200), risk_cfg.im_feedsize)
                feed_volume = feed_volume.reshape((1,)+feed_volume.shape + (1,))
                break
            return feed_volume
        except Exception as e:
            self._record_error_info(str(e))
            return None

    def _generate_figures(self, save_path, keep, coff, reg_pred):
        font={'family':'Arial',
              'style':'normal',
              'weight':'light',
              'color':'black',
              'size':13
              }

        plt.figure()
        #colors = ['#7281a7', '#af97ba', '#ef626c', '#ef626c']
        if reg_pred[-1]>0.5:
            color = 'red'
        else:
            color = 'green'

        pred_hitday = np.argmax(reg_pred)+1
        plt.plot([x for x in range(1,33)], reg_pred, color=color, label='Predicted probability distribution', linestyle='-', linewidth=3)
        plt.plot([pred_hitday, pred_hitday], [0, reg_pred[pred_hitday-1]], 'k--', lw=1, label='Predicted recovery day: ' + str(pred_hitday))

        plt.xlim([0, 33])
        plt.ylim([0, 0.5])
        plt.xlabel('Days', fontdict=font)
        plt.ylabel('Recovery probability', fontdict=font)
        plt.title('', fontdict=font)

        plt.legend(loc='upper right', prop={'size':14, 'family':'arial'})

        plt.savefig(osp.join(save_path, 'recovery_probability.png'), dpi=300)

        #saving top-5 features

        plt.figure()

        feature_names = np.array(['Gender', 'Age', 'Fever', 'Cough', 'Soreness', 'Weakness',\
                 'Headache', 'Nausea/Vomiting', 'Diarrhea', 'Expectoration',\
                 'Chest Congestion', 'Poor Appetite' ,'Poor Spirits',\
                 'Diabetes', 'Kidney Disease', 'Malignant Tumor',
                 'Chronic Hepatitis B' ,'Hypothyroidism', 'ARDS',
                 'SK' ,'IL-6' ,'WBC' ,'LAV',
                 'NAV', 'HG' ,'RBC',
                 'PC' ,'HS-CRP', 'CRP',
                 'PCT' ,'PT', 'TT',
                 'APTT' ,'FB' ,'DD', 'BGLU',
                 'ALT' ,'AST', 'TP',
                 'AM' ,'TB' ,'DB', 'UN',
                 'CRT', 'LDH', 'CK',
                 'CK-MB' ,'α-HBDH', 'CT'])

        colors = ['darkviolet', 'darkorange','aqua', 'yellowgreen']
        color_idx =[1, 1, 2, 2, 2, 2,\
                    2, 2, 2, 2,\
                    2, 2 ,2,\
                    0, 0, 0,
                    0 ,0, 0,
                    0 ,3 ,3 ,3,
                    3, 3 ,3,
                    3 ,3, 3,
                    3 ,3, 3,
                    3 ,3 ,3, 3,
                    3 ,3, 3,
                    3 ,3 ,3, 3,
                    3, 3, 3,
                    3 ,3,
                    3]
        name_colors = []
        for i in color_idx:
            name_colors.append(colors[i])
        name_colors = np.array(name_colors)

        coff = coff[keep]
        feature_names = feature_names[keep]
        name_colors = name_colors[keep]

        top10_ind = np.argsort(-np.array(coff))[0:8]

        top10_coffe = coff[top10_ind]

        top10_names = np.array(feature_names)[top10_ind]
        top10_colors = name_colors[top10_ind]

        plt.rcdefaults()
        fig, ax = plt.subplots()

        y_pos = np.arange(len(top10_names))
        ax.barh(y_pos, top10_coffe, align='center', color=top10_colors, edgecolor='black')
        ax.set_yticks(y_pos)
        ax.set_xlim(xmax=1.0, xmin=0)
        ax.set_yticklabels(top10_names, font)
        ax.invert_yaxis()  # labels read top-to-bottom
        ax.set_xlabel('Average impact on model output', font)

        plt.tight_layout()
        plt.savefig(osp.join(save_path, 'top5.png'), dpi=300)

    def _record_error_info(self, error_info):
        logging.error(error_info)
        if self.debug_mode:
            print('!!!!!!!!!!!!!!!!!!!!!!!!!!!')
            print (error_info)

    def _record_debug_info(self, debug_info):
        if self.debug_mode:
            logging.info(debug_info)
            print (debug_info)


text:
"""
Overview
--------

general info about this module

Classes and Inheritance Structure
----------------------------------------------
.. inheritance-diagram::

Summary
---------
.. autosummary::
   list of the module you want

Module API
----------
"""

from __future__ import absolute_import, division, print_function

from builtins import (bytes, str, open, super, range,
                      zip, round, input, int, pow, object, map, zip)

__author__ = "Andrea Tramacere"

# Standard library
# eg copy
# absolute import rg:from copy import deepcopy

# Dependencies
# eg numpy
# absolute import eg: import numpy as np

# Project
# relative import eg: from .mod import f
import yaml
import json

import threading
import time
import subprocess, os

class MyThread(threading.Thread):
    def run(self):
        dir_path = os.path.dirname(os.path.realpath(__file__))
        cmd = ["python", dir_path + '/'+'mock_worker.py N=%d job_id=%s session_id=%s scratch_dir=%s'%(self.N,self.job_id,self.session_id,self.scratch_dir)]
        print("command:", " ".join(cmd))
        os.spawnl(os.P_NOWAIT, cmd )

        # #mimic call_back
        # job_status={}
        #
        # job_status['status']='submitted'
        # job_status['fraction']=0.0
        # f_path=self.scratch_dir+'/'+'status.yml'
        # for i in range(self.N):
        #     print('i',i)
        #     time.sleep(1)
        #     job_status['fraction']=float(i)/(self.N-1)
        #     with open(f_path, 'w') as outfile:
        #         yaml.dump(job_status, outfile, default_flow_style=False)
        # #
        # job_status['status'] = 'done'
        # with open(f_path, 'w') as outfile:
        #     yaml.dump(job_status, outfile, default_flow_style=False)

        return

def spawn(N,job_id,session_id,scratch_dir):
    dir_path = os.path.dirname(os.path.realpath(__file__))
    cmd = ["python", dir_path + '/' + 'mock_worker.py',"%s %s %s %d & " % (job_id, session_id, scratch_dir,N) ]
    print("command:", " ".join(cmd))

    os.system(" ".join(cmd))

def mock_query(par_dic,session_id,job_id,scratch_dir):

    job_status = par_dic['job_status']
    session_id = par_dic['session_id']
    products=''

    if job_status == 'new':
        print('New Job --> id,session,dir',job_id, session_id, scratch_dir)

        job_status = submit_new_job(par_dic,session_id,job_id,scratch_dir)

    if job_status == 'submitted' or job_status == 'unacessible':
        print('Job Check --> id,session,dir', job_id, session_id, scratch_dir)
        job_status = mock_chek_job_status(job_id=job_id, session_id=session_id, scratch_dir=scratch_dir)

    if job_status == 'done':
        print('Job Done --> id,session,dir', job_id, session_id, scratch_dir)
        job_status = mock_chek_job_status(job_id=job_id, session_id=session_id, scratch_dir=scratch_dir)
        products = 'HELLO WORLD'

    return build_response(job_id,job_status,products=products)

def build_response(job_id,job_status,products='',exit_status=''):
    out_dict = {}
    out_dict['products'] = ''
    out_dict['exit_status'] = ''
    out_dict['job_id'] =job_id
    out_dict['job_status'] = job_status['status']
    out_dict['job_fraction'] = job_status['fraction']
    return out_dict

def submit_new_job(par_dic,job_id,session_id,scratch_dir):
    spawn(20, job_id, session_id, scratch_dir)
    job_status = {}
    job_status['status'] = 'submitted'
    job_status['fraction'] = ''
    f_path=scratch_dir+'/'+'query.yaml'
    with open(f_path, 'w') as outfile:
        my_json_str=json.dumps(par_dic,encoding='utf-8')
        #if isinstance(my_json_str, str):
        outfile.write(u'%s'%my_json_str)

    return job_status

def mock_chek_job_status(job_id,session_id,scratch_dir):
    f_path = scratch_dir + '/' + 'status.yml'
    print ('f_path',f_path)
    try:
        with open(f_path, 'r') as outfile:
            job_status=json.load(outfile,encoding='utf-8')
        print ('-->',job_status)
    except:
        job_status = {}
        job_status['status'] = 'unacessible'
        job_status['fraction'] = ''

    return job_status


text:
from matplotlib.widgets import RectangleSelector
from matplotlib.widgets import LassoSelector
import matplotlib.colors as clr
from matplotlib.path import Path
from matplotlib.patches import Rectangle, Ellipse,Polygon
import numpy as np
from matplotlib.ticker import MaxNLocator

class Highlighter(object):

    def __init__(self, ax_umap, ax_year, ax_mesh, dfpapers, colors='c', method='lasso', add_text=False,index_list=None,include_nlp=False):
        self.ax_umap = ax_umap
        self.ax_year = ax_year
        self.ax_mesh = ax_mesh
        self.canvas_umap = ax_umap.figure.canvas
        self.canvas_year = ax_year.figure.canvas
        self.canvas_mesh = ax_mesh.figure.canvas
        self.x, self.y = dfpapers['umap-x'], dfpapers['umap-y']
        self.add_text = add_text
        self.df_papers = dfpapers
        self.df_selected = None
        #self.selected_text = None
        self.full_index = dfpapers.index.values
        highlight_color = colors
        highlight_size = 20
        self.object_index = None
        self.mask = np.zeros(self.x.shape[0], dtype=bool)
        self.scaled_rgb = (65 / 256., 91 / 256., 169 / 256.)
        self.selection_method = method
        self.x0 = None
        self.y0 = None
        self.x1 = None
        self.y1 = None
        self.entities_dict = {}
        self._highlight = ax_umap.scatter(
           [],[], marker="o", lw=0, s=highlight_size, c=highlight_color
        )
        self.include_nlp = include_nlp
        if include_nlp:
            self.nlp = spacy.load('en_ner_bionlp13cg_md')
            self.nlp.add_pipe(self.nlp.create_pipe('sentencizer'))
        #if self.selected_text != None:
        #    self.selected_text.remove()

        if self.add_text:
            self.selected_text = self.ax_umap.text(
                0.03,
                0.97,
                "%3.2f (%i events)" % (0, 0),
                ha="left",
                va="top",
                transform=self.ax_umap.transAxes,
                color="k",
                fontname="Helvetica",
            )

        self.name_to_selector = {"lasso": LassoSelector}
        selector = self.name_to_selector[method.lower()]

        onselect_dict = {
            LassoSelector: self._onselect_lasso,
        }

        self.lasso = selector(self.ax_umap, onselect_dict[selector])

    # def remove_text(self):

    def disconnect(self):
        self.canvas_umap.mpl_disconnect(self.lasso)

    def _draw_points(self):
        self.object_index = self.full_index[self.mask]
        self.nevents = len(self.mask)
        self.nselected = np.sum(self.mask)
        if self.add_text:
            self.selected_text.remove()

        if self.add_text:
            self.selected_text = self.ax_umap.text(
                0.03,
                0.97,
                "%3.2f (%i events)"
                % (self.nselected * 100. / self.nevents, self.nselected),
                ha="left",
                va="top",
                transform=self.ax_umap.transAxes,
                color="k",
                fontname="Helvetica",
            )
        xy = np.column_stack([self.x[self.mask], self.y[self.mask]])
        self._highlight.set_offsets(xy)
        self.canvas_umap.draw()

    def _onselect_lasso(self, verts):
        self.verts = np.array(verts)
        p = Path(self.verts)
        pix = np.array([self.x, self.y]).T
        self.mask = p.contains_points(pix, radius=0)
        self._draw_points()
        self.df_selected = self.df_papers[self.mask]
        print("%s selected" % self.df_selected.shape[0])
        if self.df_selected.shape[0] == 0:
            self.entities_dict = {}

        self.ax_mesh.clear()
        self.ax_year.clear()

        mterms = []
        for idx in self.df_selected.index:
            yp = self.df_selected.loc[idx]['pubdate']
            mt = self.df_selected.loc[idx]['mesh_terms']
            t = self.df_selected.loc[idx]['title']
            #sdoc = self.nlp(self.df_selected.loc[idx]['abstract'])
            #sen = list(sdoc.sents)
            mterms.extend(sorted([mti.split(":")[-1] for mti in mt.split(";")]))
            print("(%s) %s" % (yp,t))
            doci = self.df_selected.loc[idx]['title'] + self.df_selected.loc[idx]['abstract']
            if self.include_nlp:
                sdoc = nlp(doci)
                #for si in sdoc.ents:
                self.entities_dict[idx] = sdoc.ents
            #print(" Mesh:",", ".join(mterms))
#            print(" ".join(item.split()[:100]))
        self.df_selected.to_csv("tmp_selected.csv")

        self.mesh_counts = Counter(mterms)
        self.mesh_counts_top = self.mesh_counts.most_common(50)
        dfcounts = pd.DataFrame.from_dict(dict(self.mesh_counts_top), orient='index')
        dfcounts.sort_values(0,inplace=True)
        if dfcounts.shape[0]!=0:
            dfcounts.plot(kind='barh',ax=self.ax_mesh,fontsize=6,legend=False)

        xmin = np.min(self.df_selected['pubdate'])
        xmax = np.max(self.df_selected['pubdate'])
        self.ax_year.hist(self.df_selected['pubdate'],bins=np.arange(xmin,xmax+1)-0.5)
        self.ax_year.set_xlim([xmin,xmax])
        self.ax_year.xaxis.set_major_locator(MaxNLocator(integer=True))
        for label in self.ax_year.get_xticklabels():
            label.set_rotation(45)
            label.set_fontsize(6)
            label.set_ha('right')

        self.canvas_year.draw()
        self.canvas_mesh.draw()

max_stars_repo_name: bgriffen/pubmedpandas

max_issues_repo_name: bgriffen/pubmedpandas

max_forks_repo_name: bgriffen/pubmedpandas


text:
import numpy as np

import torch
import torchvision.utils as vutils

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
#import pandas as pd
import seaborn as sns
sns.set()
sns.set_style('whitegrid')
sns.set_palette('colorblind')

def convert_npimage_torchimage(image):
    return 255*torch.transpose(torch.transpose(torch.from_numpy(image), 0, 2), 1, 2)

def get_scatter_plot(data, labels=None, n_classes=None, num_samples=1000, xlim=None, ylim=None):
    '''
    data   : 2d points, batch_size x data_dim (numpy array)
    labels : labels, batch_size (numpy array)
    '''
    batch_size, data_dim = data.shape
    num_samples = min(num_samples, batch_size)
    if labels is None:
        labels = np.zeros(batch_size, dtype=np.int)
    if n_classes is None:
        n_classes = len(np.unique(labels))

    # sub-samples
    if num_samples != batch_size:
        indices = np.random.permutation(batch_size)
        data = data[indices[:num_samples]]
        labels = labels[indices[:num_samples]]

    # init config
    palette = sns.color_palette(n_colors=n_classes)
    palette = [palette[i] for i in np.unique(labels)]

    # plot
    fig, ax = plt.subplots(figsize=(5, 5))
    data = {'x': data[:, 0],
            'y': data[:, 1],
            'class': labels}
    sns.scatterplot(x='x', y='y', hue='class', data=data, palette=palette)

    # set config
    if xlim is not None:
        plt.xlim((-xlim, xlim))
    if ylim is not None:
        plt.ylim((-ylim, ylim))

    # draw to canvas
    fig.canvas.draw()  # draw the canvas, cache the renderer
    image = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))

    # close figure
    plt.close()
    return image

def get_data_for_quiver_plot(val, num):
    _x = np.linspace(-val, val, num)
    _y = np.linspace(-val, val, num)
    _u, _v = np.meshgrid(_x, _y)
    _vis_data = np.stack([_u.reshape(num**2), _v.reshape(num**2)], axis=1)
    vis_data = torch.from_numpy(_vis_data).float()
    return vis_data, _x, _y

def get_quiver_plot(vec, x_pos, y_pos, xlim=None, ylim=None, scale=None):
    '''
    vec   : 2d points, batch_size x data_dim (numpy array)
    pos   : 2d points, batch_size x data_dim (numpy array)
    '''
    grid_size = x_pos.shape[0]
    batch_size = vec.shape[0]
    assert batch_size == grid_size**2
    assert y_pos.shape[0] == grid_size

    # get x, y, u, v
    X = x_pos #np.arange(-10, 10, 1)
    Y = y_pos #np.arange(-10, 10, 1)
    #U, V = np.meshgrid(X, Y)
    U = vec[:, 0].reshape(grid_size, grid_size)
    V = vec[:, 1].reshape(grid_size, grid_size)

    # plot
    fig, ax = plt.subplots(figsize=(5, 5))
    q = ax.quiver(X, Y, U, V, pivot='mid', scale=scale)
    #ax.quiverkey(q, X=0.3, Y=1.1, U=10,
    #                     label='Quiver key, length = 10', labelpos='E')

    # set config
    if xlim is not None:
        plt.xlim((-xlim, xlim))
    if ylim is not None:
        plt.ylim((-ylim, ylim))

    # tight
    plt.tight_layout()

    # draw to canvas
    fig.canvas.draw()  # draw the canvas, cache the renderer
    image = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))

    # close figure
    plt.close()
    return image

def get_data_for_heatmap(val=4, num=256):
    _x = np.linspace(-val, val, num)
    _y = np.linspace(-val, val, num)
    _u, _v = np.meshgrid(_x, _y)
    _data = np.stack([_u.reshape(num**2), _v.reshape(num**2)], axis=1)
    return _data, _x, _y

def energy_to_unnormalized_prob(energy):
    prob = torch.exp(-energy) # unnormalized prob
    return prob

def get_prob_from_energy_func_for_vis(energy_func, val=4, num=256):
    # get grid
    _z, _, _ = get_data_for_heatmap(val=val, num=num)
    z = torch.from_numpy(_z).float()

    # run energy func
    energy = energy_func(z)
    prob = energy_to_unnormalized_prob(energy)

    # convert to numpy array
    _prob = prob.cpu().float().numpy()
    _prob = _prob.reshape(num, num)
    return _prob

def get_imshow_plot(prob, val=4, use_grid=True):
    # plot
    fig, ax = plt.subplots(figsize=(5, 5))
    im = ax.imshow(prob, cmap='jet', extent=[-val, val, -val, val])
    ax.grid(False)
    if use_grid:
        plt.xticks(np.arange(-val, val+1, step=1))
        plt.yticks(np.arange(-val, val+1, step=1))
    else:
        plt.xticks([])
        plt.yticks([])

    # tight
    plt.tight_layout()

    # draw to canvas
    fig.canvas.draw()  # draw the canvas, cache the renderer
    image = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))

    # close figure
    plt.close()
    return image

def get_1d_histogram_plot(data, val=4, num=256, use_grid=True):
    xmin = 0
    xmax = val

    # get data
    x = data

    # get histogram
    hist, xedges = np.histogram(x, range=[xmin, xmax], bins=num)

    # plot heatmap
    fig, ax = plt.subplots(figsize=(5, 5))
    im = ax.bar(xedges[:-1], hist, width=0.5)#, color='#0504aa',alpha=0.7)

    ax.grid(False)
    if use_grid:
        plt.xticks(np.arange(0, val+1, step=1))
    else:
        plt.xticks([])

    # tight
    plt.tight_layout()

    # draw to canvas
    fig.canvas.draw()  # draw the canvas, cache the renderer
    image = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))

    # close figure
    plt.close()
    return image

def get_2d_histogram_plot(data, val=4, num=256, use_grid=True):
    xmin = -val
    xmax = val
    ymin = -val
    ymax = val

    # get data
    x = data[:, 0]
    y = data[:, 1]

    # get histogram
    heatmap, xedges, yedges = np.histogram2d(x, y, range=[[xmin, xmax], [ymin, ymax]], bins=num)
    extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

    # plot heatmap
    fig, ax = plt.subplots(figsize=(5, 5))
    im = ax.imshow(heatmap.T, extent=extent, cmap='jet')
    ax.grid(False)
    if use_grid:
        plt.xticks(np.arange(-val, val+1, step=1))
        plt.yticks(np.arange(-val, val+1, step=1))
    else:
        plt.xticks([])
        plt.yticks([])

    # tight
    plt.tight_layout()

    # draw to canvas
    fig.canvas.draw()  # draw the canvas, cache the renderer
    image = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))

    # close figure
    plt.close()
    return image

def get_grid_image(input, batch_size, nchannels, nheight, nwidth=None, ncol=8, pad_value=0, do_square=True):
    '''
    input : b x c x h x w (where h = w)
    '''
    if batch_size > ncol**2 and do_square:
        input = input[:ncol**2, :, :, :]
        batch_size = ncol**2
    nwidth = nwidth if nwidth is not None else nheight
    input = input.detach()
    output = input.view(batch_size, nchannels, nheight, nwidth).clone().cpu()
    output = vutils.make_grid(output, nrow=ncol, normalize=True, scale_each=True, pad_value=pad_value)
    #output = vutils.make_grid(output, normalize=False, scale_each=False)
    return output

#def get_canvas(fig):
#    fig.canvas.draw()  # draw the canvas, cache the renderer
#    image = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
#    image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))
#    return image
#    # close figure
#    plt.close()
#
#def get_contour_with_batch_size(model, batch_size=128, vmin=-10.0, vmax=10.0, title=None):
#    model.eval()
#    matplotlib.rcParams['xtick.direction'] = 'out'
#    matplotlib.rcParams['ytick.direction'] = 'out'
#    matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
#
#    # tmp
#    weight = next(model.parameters())
#
#    # gen grid⋅
#    delta = 0.1
#    xv, yv = torch.meshgrid([torch.arange(vmin, vmax, delta), torch.arange(vmin, vmax, delta)])
#    h = yv.size(0)
#    w = xv.size(0)
#    yv = yv.contiguous().view(-1)
#    xv = xv.contiguous().view(-1)
#    input = torch.cat([xv.unsqueeze(1), yv.unsqueeze(1)], dim=1).to(weight.device)
#
#    # forward
#    prob = model.prob(input, batch_size=batch_size)
#
#    # convert torch variable to numpy array
#    xv = xv.cpu().numpy().reshape(h, w)
#    yv = yv.cpu().numpy().reshape(h, w)
#    zv = prob.detach().cpu().numpy().reshape(h, w)
#
#    # plot and save⋅
#    fig = plt.figure()
#    CS1 = plt.contourf(xv, yv, zv)
#    CS2 = plt.contour(xv, yv, zv, alpha=.7, colors='k')
#    plt.clabel(CS2, inline=1, fontsize=10, colors='k')
#    #plt.title('Simplest default with labels')
#    if title is not None:
#        plt.title(title)
#    #plt.savefig(filename)
#    #plt.close()
#    image = get_canvas(fig)
#    plt.close()
#
#    return image
#
##def get_contour_with_data(model, data, vmin=-10.0, vmax=10.0, title=None):
##    model.eval()
##    matplotlib.rcParams['xtick.direction'] = 'out'
##    matplotlib.rcParams['ytick.direction'] = 'out'
##    matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
##
##    # gen grid⋅
##    delta = 0.1
##    xv, yv = torch.meshgrid([torch.arange(vmin, vmax, delta), torch.arange(vmin, vmax, delta)])
##    h = yv.size(0)
##    w = xv.size(0)
##    yv = yv.contiguous().view(-1)
##    xv = xv.contiguous().view(-1)
##    input = torch.cat([xv.unsqueeze(1), yv.unsqueeze(1)], dim=1).to(data.device)
##
##    # forward
##    prob = model.prob(input, data)
##
##    # convert torch variable to numpy array
##    xv = xv.cpu().numpy().reshape(h, w)
##    yv = yv.cpu().numpy().reshape(h, w)
##    zv = prob.detach().cpu().numpy().reshape(h, w)
##
##    # plot and save⋅
##    fig = plt.figure()
##    CS1 = plt.contourf(xv, yv, zv)
##    CS2 = plt.contour(xv, yv, zv, alpha=.7, colors='k')
##    plt.clabel(CS2, inline=1, fontsize=10, colors='k')
##    #plt.title('Simplest default with labels')
##    if title is not None:
##        plt.title(title)
##    #plt.savefig(filename)
##    #plt.close()
##    image = get_canvas(fig)
##    plt.close()
##
##    return image
#
#def get_contour_with_z(model, z, vmin=-10.0, vmax=10.0, title=None):
#    model.eval()
#    matplotlib.rcParams['xtick.direction'] = 'out'
#    matplotlib.rcParams['ytick.direction'] = 'out'
#    matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
#
#    # gen grid⋅
#    delta = 0.1
#    xv, yv = torch.meshgrid([torch.arange(vmin, vmax, delta), torch.arange(vmin, vmax, delta)])
#    h = yv.size(0)
#    w = xv.size(0)
#    yv = yv.contiguous().view(-1)
#    xv = xv.contiguous().view(-1)
#    input = torch.cat([xv.unsqueeze(1), yv.unsqueeze(1)], dim=1).to(z.device)
#
#    # forward
#    prob = model.prob(input, z=z)
#
#    # convert torch variable to numpy array
#    xv = xv.cpu().numpy().reshape(h, w)
#    yv = yv.cpu().numpy().reshape(h, w)
#    zv = prob.detach().cpu().numpy().reshape(h, w)
#
#    # plot and save⋅
#    fig = plt.figure()
#    CS1 = plt.contourf(xv, yv, zv)
#    CS2 = plt.contour(xv, yv, zv, alpha=.7, colors='k')
#    plt.clabel(CS2, inline=1, fontsize=10, colors='k')
#    #plt.title('Simplest default with labels')
#    if title is not None:
#        plt.title(title)
#    #plt.savefig(filename)
#    #plt.close()
#    image = get_canvas(fig)
#    plt.close()
#
#    return image


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright 1999-2018 Alibaba Group Holding Ltd.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import numpy as np

from .... import operands
from ..utils import infer_dtype
from .core import TensorBinOp, TensorConstant
from .utils import arithmetic_operand

@arithmetic_operand(sparse_mode='always_false')
class TensorFloorDiv(operands.FloorDiv, TensorBinOp):
    @classmethod
    def constant_cls(cls):
        return TensorFDivConstant

@arithmetic_operand(sparse_mode='binary_or_const')
class TensorFDivConstant(operands.FDivConstant, TensorConstant):
    @classmethod
    def _is_sparse(cls, x1, x2):
        if hasattr(x1, 'issparse') and x1.issparse():
            if x2 != 0:
                return True
            else:
                raise ZeroDivisionError('float division by zero')
        return False

@infer_dtype(np.floor_divide)
def floordiv(x1, x2, out=None, where=None, **kwargs):
    """
    Return the largest integer smaller or equal to the division of the inputs.
    It is equivalent to the Python ``//`` operator and pairs with the
    Python ``%`` (`remainder`), function so that ``b = a % b + b * (a // b)``
    up to roundoff.

    Parameters
    ----------
    x1 : array_like
        Numerator.
    x2 : array_like
        Denominator.
    out : Tensor, None, or tuple of Tensor and None, optional
        A location into which the result is stored. If provided, it must have
        a shape that the inputs broadcast to. If not provided or `None`,
        a freshly-allocated array is returned. A tuple (possible only as a
        keyword argument) must have length equal to the number of outputs.
    where : array_like, optional
        Values of True indicate to calculate the ufunc at that position, values
        of False indicate to leave the value in the output alone.
    **kwargs

    Returns
    -------
    y : Tensor
        y = floor(`x1`/`x2`)

    See Also
    --------
    remainder : Remainder complementary to floor_divide.
    divmod : Simultaneous floor division and remainder.
    divide : Standard division.
    floor : Round a number to the nearest integer toward minus infinity.
    ceil : Round a number to the nearest integer toward infinity.

    Examples
    --------
    >>> import mars.tensor as mt

    >>> mt.floor_divide(7,3).execute()
    2
    >>> mt.floor_divide([1., 2., 3., 4.], 2.5).execute()
    array([ 0.,  0.,  1.,  1.])
    """
    op = TensorFloorDiv(**kwargs)
    return op(x1, x2, out=out, where=where)

@infer_dtype(np.floor_divide, reverse=True)
def rfloordiv(x1, x2, **kwargs):
    op = TensorFloorDiv(**kwargs)
    return op.rcall(x1, x2)

max_stars_repo_name: ChenQuan/mars

max_issues_repo_name: ChenQuan/mars

max_forks_repo_name: ChenQuan/mars


text:
"""
Comments on notation found in the appendix.

Notation
--------
  P :
    A set of players, where P_i is one player.
  h :
    A node (i.e history), h is defined by all information of the current
    situation, including private knowledge known only to one player.
  A(h) :
    Denotes the actions available at a node
  P(h) :
    Denotes either the chance or player it is to act at that node.
  I :
    Imperfect information is represented by information sets (infosets) for
    each player P_i. For any infoset belonging to a player P_i, all nodes h to
    h' in the infoset I are indestinguishable to the player P_i. Moreover,
    every non-terminal node h in belongs to exactly one infoset for each
    player.
  sigma (lowercase) :
    A strategy (i.e a policy). Here sigma(I) is a probability vector over the
    actions for acting player P_i in infoset I.
"""
import numpy as np

def monte_carlo_cfr_with_pruning(t):
    """Conduct External-Sampling Monte Carlo CFR with pruning."""
    for player in players:
        pass

def calculate_strategy(R, I):
    """Caluclates the strategy based on regrets."""
    pass

def update_strategy(h, P):
    """Update the average strategy of P_i"""
    pass

def traverse_monte_carlo_cfr(h, P):
    """"Update the regrets for P_i"""
    pass

def traverse_monte_carlo_cfr_with_pruning(h, P):
    """MCCFR with pruning for very negative regrets."""
    pass

max_stars_repo_name: ZgodaRafal/pluribus-poker-AI

max_issues_repo_name: ZgodaRafal/pluribus-poker-AI

max_forks_repo_name: ZgodaRafal/pluribus-poker-AI


text:
import numpy as np
from PIL import Image
import colorsys

"""
Methods implementing the COLOR metric.
"""

def generate_bins() -> list:
    """
    Generate color bins.
    :return: List of bins
    """
    h_bins = [(x / 10.0, (x + 1) / 10.0) for x in range(0, 10)]
    h_bins[-1] = (h_bins[-1][0], 1.1)
    s_bins = [(0.0, 0.333), (0.333, 0.666), (0.666, 1.1)]
    l_bins = [(0.0, 0.333), (0.333, 0.666), (0.666, 1.1)]
    bins = []
    for h_bin in h_bins:
        for s_bin in s_bins:
            for l_bin in l_bins:
                bins.append((h_bin, s_bin, l_bin))
    return bins

def to_hsl(img: Image.Image) -> np.ndarray:
    """
    Converts an image to HSL format and return as a numpy array.
    :param img: Image to convert
    :return: Converted image.
    """
    img_rgb = img.convert("RGB")
    img_rgb_np = np.array(img_rgb) / 255.0
    rgb_hls = np.vectorize(colorsys.rgb_to_hls)
    np_h, np_l, np_s = rgb_hls(r=img_rgb_np[:, :, 0], g=img_rgb_np[:, :, 1], b=img_rgb_np[:, :, 2])
    img_hsl_np = np.concatenate([np.expand_dims(np_h, 2), np.expand_dims(np_s, 2), np.expand_dims(np_l, 2)], axis=2)
    return img_hsl_np

def compute_histogram(img_hsl_np: np.ndarray, bins: list) -> np.ndarray:
    """
    Computes HSL histogram for a given image in HSL color space.
    :param img_hsl_np: Image in HSL color space.
    :param bins: Bin specification.
    :return: Histogram
    """
    freqs = []
    for h_bound, s_bound, l_bound in bins:
        count = np.logical_and(
            np.logical_and(
                np.logical_and(img_hsl_np[:, :, 0] >= h_bound[0], img_hsl_np[:, :, 0] < h_bound[1]),
                np.logical_and(img_hsl_np[:, :, 1] >= s_bound[0], img_hsl_np[:, :, 1] < s_bound[1])),
            np.logical_and(img_hsl_np[:, :, 2] >= l_bound[0], img_hsl_np[:, :, 2] < l_bound[1])).sum()
        freqs.append(count)
    return np.array(freqs)

def hsl_hist_l1(pred: Image, gt:Image, bins:list):
    """
    Compute COLOR metric between prediction and ground truth using the specified bins.
    :param pred: Predicted texture
    :param gt: Ground truth texture
    :param bins: Bin specification returned from generate_bins
    :return: Metric value.
    """
    pred_hist = compute_histogram(to_hsl(pred), bins)
    gt_hist = compute_histogram(to_hsl(gt), bins)
    dif = gt_hist - pred_hist
    score = np.sum(np.abs(dif)) / (2.0 * pred.width * pred.height)
    return score