text:
#!/usr/bin/env python

"""
description: script to compare pwm hits
in grammar regions vs non-grammar regions

"""

import os
import sys
import h5py
import logging
import argparse

import networkx as nx
import numpy as np
import pandas as pd

from numpy.random import RandomState

from tronn.util.scripts import setup_run_logs
from tronn.util.utils import DataKeys

def parse_args():
    """parser
    """
    parser = argparse.ArgumentParser(
        description="annotate grammars with functions")

    parser.add_argument(
        "--grammar",
        help="grammar gml file")
    parser.add_argument(
        "--score_files", nargs="+",
        help="h5 file with pwm scores (raw sequence)")
    parser.add_argument(
        "--synergy_only", action="store_true",
        help="filter for synergy only")
    parser.add_argument(
        "--synergy_file",
        help="h5 file with synergy scores")
    parser.add_argument(
        "--compare_keys", nargs="+", default=["ATAC_SIGNALS.NORM", "H3K27ac_SIGNALS.NORM"],
        help="which score keys to compare")
    parser.add_argument(
        "--compare_bigwigs", nargs="+", default=[],
        help="other bigwigs for extracting signals for comparison")
    parser.add_argument(
        "--map_chain_file",
        help="liftover file")
    parser.add_argument(
        "--eqtl_bed_file",
        help="eqtl BED file to see if any eQTLs in grammar/null regions (rsid must be in name col)")
    parser.add_argument(
        "--eqtl_effects_file",
        help="table with effect sizes")
    parser.add_argument(
        "--n_null", type=int, default=100,
        help="number of null distributions to try")
    parser.add_argument(
        "-o", "--out_dir", dest="out_dir", type=str,
        default="./",
        help="out directory")
    parser.add_argument(
        "--prefix",
        help="file prefix")

    # parse args
    args = parser.parse_args()

    return args

def get_bed_from_metadata_list(examples, bed_file, interval_key="active", merge=True):
    """make a bed file from list of metadata info
    """
    with open(bed_file, "w") as fp:
        for region_metadata in examples:
            interval_types = region_metadata.split(";")
            interval_types = dict([
                interval_type.split("=")[0:2]
                for interval_type in interval_types])
            interval_string = interval_types[interval_key]

            chrom = interval_string.split(":")[0]
            start = interval_string.split(":")[1].split("-")[0]
            stop = interval_string.split("-")[1]
            fp.write("{}\t{}\t{}\n".format(chrom, start, stop))

    if merge:
        tmp_bed_file = "{}.tmp.bed".format(bed_file.split(".bed")[0])
        os.system("mv {} {}".format(bed_file, tmp_bed_file))
        os.system((
            "cat {} | "
            "sort -k1,1 -k2,2n | "
            "bedtools merge -i stdin | "
            "gzip -c > {}").format(
                tmp_bed_file, bed_file))
        os.system("rm {}".format(tmp_bed_file))

    return

def get_signal_from_bigwig(bed_file, bigwig_file, map_chain_file=None, liftover=True):
    """
    """
    # if liftover, first convert the bed file to correct coords
    out_dir = os.path.dirname(bed_file)
    tmp_bed = "{}/liftover.tmp.bed".format(out_dir)
    unmapped_file = "{}/unmapped.txt".format(out_dir)
    if liftover:
        liftover = "liftOver {} {} {} {}".format(
            bed_file, map_chain_file, tmp_bed, unmapped_file)
        os.system(liftover)
    else:
        tmp_bed = bed_file

    # give each region a name
    new_bed = "{}/signal.bed".format(out_dir)
    name_regions = (
        "cat {} | "
        "awk -F '\t' '{{ print $0\"\tregion\"NR }}' "
        "> {}").format(tmp_bed, new_bed)
    os.system(name_regions)

    # then bigWigAverageOverBed
    avg_signal_file = "{}/signal.avg.tmp.txt".format(out_dir)
    get_signal = "bigWigAverageOverBed {} {} {}".format(
        bigwig_file, new_bed, avg_signal_file)
    os.system(get_signal)

    # then read in data
    avg_signal = pd.read_csv(avg_signal_file, sep="\t", header=None)
    avg_signal = avg_signal.iloc[:,4].values
    #print "{}: median signal {}".format(bigwig_file, np.median(avg_signal))

    return avg_signal

def get_overlapping_variants(
        bed_file,
        variant_bed_file,
        lookup_table,
        variant_type="GTEx"):
    """get variants that are in BED regions as well as effect sizes
    """
    tmp_dir = os.path.dirname(bed_file)
    # get overlapping variants
    tmp_overlap_file = "{}/overlap.tmp.txt.gz".format(tmp_dir)
    overlap_cmd = (
        "bedtools intersect -wa -a {} -b {} | "
        "awk -F '\t' '{{ print $4 }}' | "
        "gzip -c > {}").format(
            variant_bed_file, bed_file, tmp_overlap_file)
    os.system(overlap_cmd)
    variant_ids = pd.read_csv(tmp_overlap_file, header=None).values[:,0]

    # get filtered set from lookup
    if variant_type == "GTEx":
        variant_results = lookup_table[lookup_table["variant_id"].isin(variant_ids)]
    else:
        variant_results = lookup_table[lookup_table["CausalSNP"].isin(variant_ids)]

    # get effect sizes
    if variant_type == "GTEx":
        col_key = "slope"
    else:
        col_key = "P VALUE"
    variant_effects = np.absolute(variant_results[col_key].values)
    #print "{}: {} variants, median effect {}".format(
    #    bed_file, variant_results.shape[0], np.median(variant_effects))

    return variant_effects

def build_score_matched_bins(target_scores, full_scores, num_increments=21):
    """build bins
    """
    # TODO ideally build intervals for both GC and scores...
    # then to build bins, do a nested loop?
    # for score interval, for gc interval:
    # collect all target that match <- track count as num_target
    # collect all that match that interval <- select num_target from that set
    # trick part is then you need to match the num target in that set
    # so then you want to do min/max scores and divide equally

    # build intervals
    sorted_scores = np.sort(target_scores)
    intervals = np.interp(
        np.linspace(0, len(sorted_scores), num=num_increments),
        range(len(sorted_scores)),
        sorted_scores)

    # go through bins and collect sequences in bins
    binned_indices = []
    for i in range(len(intervals)-1):
        range_min = intervals[i]
        range_max = intervals[i+1]
        full_indices = np.where(
            (full_scores >= range_min) &
            (full_scores < range_max))[0]
        binned_indices.append(full_indices)

    return binned_indices

def build_matched_null_set(
        target_scores,
        background_scores,
        target_gc_scores,
        background_gc_scores,
        num_increments=11,
        rand_seed=1):
    """given a set of target scores and all scores, give back a matched
    set
    """
    # build target score intervals
    sorted_scores = np.sort(target_scores)
    intervals = np.interp(
        np.linspace(0, len(sorted_scores), num=num_increments),
        range(len(sorted_scores)),
        sorted_scores)

    # build gc score intervals
    sorted_gc = np.sort(target_gc_scores)
    gc_intervals = np.interp(
        np.linspace(0, len(sorted_gc), num=num_increments),
        range(len(sorted_gc)),
        sorted_gc)

    # for each score interval
    target_counts_per_bin = []
    indices_per_bin = []
    for score_interval_idx in range(len(intervals)-1):
        score_min = intervals[score_interval_idx]
        score_max = intervals[score_interval_idx+1]
        #print score_min, score_max

        # get target examples in score range
        target_score_range_indices = np.where(
            (target_scores >= score_min) &
            (target_scores < score_max))[0]

        # get examples in score range
        score_range_indices = np.where(
            (background_scores >= score_min) &
            (background_scores < score_max))[0]
        #print background_scores.shape
        #print score_range_indices.shape

        # for each gc interval
        for gc_interval_idx in range(len(gc_intervals)-1):
            gc_min = gc_intervals[gc_interval_idx]
            gc_max = gc_intervals[gc_interval_idx+1]
            #print "  ", gc_min, gc_max

            # get examples in gc range
            target_gc_range_indices = np.where(
                (target_gc_scores >= gc_min) &
                (target_gc_scores < gc_max))[0]

            # get examples in gc range
            gc_range_indices = np.where(
                (background_gc_scores >= gc_min) &
                (background_gc_scores < gc_max))[0]

            # intersect to get those in score and GC range
            target_indices_in_bin = set(target_score_range_indices).intersection(
                set(target_gc_range_indices))
            target_counts_per_bin.append(len(target_indices_in_bin))
            #print len(target_indices_in_bin)

            background_indices_in_bin = set(score_range_indices).intersection(
                set(gc_range_indices))
            #print len(background_indices_in_bin)
            #background_indices_in_bin = background_indices_in_bin.difference(
            #    target_indices_in_bin)
            #print len(background_indices_in_bin)
            indices_per_bin.append(list(background_indices_in_bin))

    # for each bin, get scores from all scores
    rand_state = RandomState(rand_seed)
    keep_indices = []
    for bin_idx in range(len(target_counts_per_bin)):
        target_count = target_counts_per_bin[bin_idx]
        if target_count == 0:
            continue
        indices = indices_per_bin[bin_idx]
        selected = rand_state.choice(indices, size=target_count, replace=False)
        keep_indices.append(selected)
    keep_indices = np.concatenate(keep_indices)

    return keep_indices

def _get_data_from_h5_files(h5_files, key):
    """get data concatenated
    """
    data = []
    for file_idx in range(len(h5_files)):
        h5_file = h5_files[file_idx]
        with h5py.File(h5_file, "r") as hf:
            data.append(hf[key][:])
    data = np.concatenate(data)

    return data

def get_results(data_indices, metadata, args, prefix="grammar.instances"):
    """collect all results that cover this bed file
    """
    # make bed file
    bed_file = "{}/{}.bed.gz".format(args.out_dir, prefix)
    get_bed_from_metadata_list(
        metadata[data_indices],
        bed_file,
        interval_key="active",
        merge=True)

    # save results to dict
    results = {}

    # get results for compare keys (ATAC, H3K27ac)
    for key in args.compare_keys:
        signal = _get_data_from_h5_files(
            args.score_files, key)
        results[key] = signal[data_indices]

    # get grammar variants
    variants = get_overlapping_variants(
        bed_file, args.eqtl_bed_file, args.eqtl_lookup)
    results["variants"] = variants

    # for each bigwig, get signal
    for bigwig_file in args.compare_bigwigs:
        bigwig_signal = get_signal_from_bigwig(
            bed_file,
            bigwig_file,
            map_chain_file=args.map_chain_file,
            liftover=True)
        results[os.path.basename(bigwig_file)] = bigwig_signal

    return results

def get_sig(results, null_results):
    """sort nulls, figure out where results fit in
    """
    null_results_sorted = np.sort(null_results)

    idx = np.searchsorted(null_results, results)

    pval = 1 - float(idx) / len(null_results)

    return pval

def save_results(
        filename,
        results,
        null_results,
        pwm_names,
        save_str):
    """save out to file with a comment line at front
    """
    data = pd.DataFrame(data=results, columns=["signal"])
    data["variable"] = "target"

    for i in range(len(null_results)):
        null_data = pd.DataFrame(data=null_results[i], columns=["signal"])
        null_data["variable"] = pwm_names[i]
        data = pd.concat([data, null_data], axis=0)

    with open(filename, "w") as fp:
        fp.write("#{}\n".format(save_str))
        data.to_csv(fp, index=False, sep="\t")

    # just make the plot each time?
    plot_file = "{}.pdf".format(filename.split(".txt")[0])
    plot_cmd = "Rscript /datasets/inference.2019-03-12/dmim.shuffle/quick_plot.R {} {}".format(
        filename, plot_file)
    os.system(plot_cmd)

    return

def main():
    """run annotation
    """
    # set up args
    args = parse_args()
    # make sure out dir exists
    os.system("mkdir -p {}".format(args.out_dir))
    setup_run_logs(args, os.path.basename(sys.argv[0]).split(".py")[0])
    prefix = "{}/{}".format(args.out_dir, args.prefix)

    # read in grammar
    grammar = nx.read_gml(args.grammar)
    grammar_instances = grammar.graph["examples"].split(",")

    # to consider, subset those that are diff?
    if args.synergy_only:
        with h5py.File(args.synergy_file, "r") as hf:
            synergy_diff = hf["{}.0".format(DataKeys.SYNERGY_DIFF_SIG)][:]
            diff_indices = np.where(
                np.any(synergy_diff!=0, axis=1))[0]
            synergy_metadata = hf[DataKeys.SEQ_METADATA][:,0][diff_indices]
        synergy_metadata = synergy_metadata[np.isin(synergy_metadata, grammar_instances)]
        grammar_instances = synergy_metadata

    # get metadata and get grammar indices
    metadata = _get_data_from_h5_files(
        args.score_files,
        DataKeys.SEQ_METADATA)[:,0]
    grammar_indices = np.where(
        np.isin(metadata, grammar_instances))[0]

    # get pwm scores and GC content
    pwm_raw_scores = _get_data_from_h5_files(
        args.score_files,
        DataKeys.ORIG_SEQ_PWM_SCORES_SUM)
    gc_scores = _get_data_from_h5_files(
        args.score_files,
        DataKeys.GC_CONTENT)

    # load eqtl table
    args.eqtl_lookup = pd.read_csv(args.eqtl_effects_file, sep="\t")

    # get grammar results
    results = get_results(
        grammar_indices, metadata, args, prefix="grammar.instances")
    #for key in sorted(results.keys()):
    #    print key, results[key].shape, np.median(results[key])

    # for each pwm in grammar, get null distributions and run tests
    null_results_all = []
    plot_data = None
    for pwm in grammar.nodes:
        print "\nRunning {}\n".format(pwm)
        pwm_idx = int(grammar.nodes[pwm]["pwmidx"])

        # get target scores (GC and pwms)
        target_gc_scores = gc_scores[grammar_indices]
        target_pwm_scores = pwm_raw_scores[grammar_indices,0,pwm_idx]
        background_pwm_scores = pwm_raw_scores[:,0,pwm_idx]

        # set up null store
        nulls = {}
        for key in results.keys():
            nulls[key] = []

        # build n null sets
        for null_idx in range(args.n_null):

            # logging
            if null_idx % 10 == 0:
                print null_idx

            # get null set
            matched_background_indices = build_matched_null_set(
                target_pwm_scores,
                background_pwm_scores,
                target_gc_scores,
                gc_scores,
                rand_seed=null_idx)

            # get null results
            null_results = get_results(
                matched_background_indices,
                metadata,
                args,
                prefix="grammar.non_instances")
            #for key in sorted(null_results.keys()):
            #    print key, null_results[key].shape, np.median(null_results[key])
            #print ""

            # save out null results
            for key in sorted(null_results.keys()):
                nulls[key].append(null_results[key])

        # finally save out to overall list
        null_results_all.append(nulls)

        if False:
            # quick test save for plot
            test_idx = 10 # 10
            plot_background = background_signal[:,test_idx]
            plot_background_data = pd.DataFrame(data=plot_background, columns=["signal"])
            plot_background_data["variable"] = pwm

            if plot_data is None:
                plot_data = plot_background_data.copy()
            else:
                plot_data = pd.concat([plot_data, plot_background_data], axis=0)

            print np.mean(target_signal, axis=0)
            print np.mean(background_signal, axis=0)

            print ""

    # get medians, calc significance
    for key in sorted(results.keys()):

        file_prefix = "{}/results.{}".format(args.out_dir, key.replace(".", "_"))

        pval_thresh = 0.05
        print key

        grammar_results = results[key]
        grammar_median = np.median(grammar_results, axis=0)

        if len(grammar_median.shape) > 0:
            # split out
            for i in range(grammar_median.shape[0]):
                split_file_prefix = "{}.idx-{}".format(file_prefix, i)
                null_sigs = []
                null_median_medians = []
                save_nulls = []
                split_grammar_results = grammar_results[:,i]
                split_grammar_median = grammar_median[i]
                for null_idx in range(len(null_results_all)):
                    null_results = null_results_all[null_idx][key]
                    split_null_results = [
                        vals[:,i]
                        for vals in null_results]
                    split_null_medians = [
                        np.median(val, axis=0)
                        for val in split_null_results]
                    pval = get_sig(split_grammar_median, split_null_medians)
                    null_sigs.append(pval)
                    # get the median null dist
                    if len(split_null_medians) % 2 == 0:
                        split_null_medians = split_null_medians[1:]
                    null_median_median = np.median(split_null_medians)
                    null_median_medians.append(null_median_median)
                    save_idx = split_null_medians.index(null_median_median)
                    save_nulls.append(split_null_results[save_idx])

                #print split_grammar_median, null_median_medians, null_sigs
                filename = "{}.txt".format(split_file_prefix)
                save_str = "filename={};pwms={};grammar_avg={};null_avgs={},null_sigs={}".format(
                    filename,
                    grammar.nodes,
                    split_grammar_median,
                    null_median_medians,
                    null_sigs)
                print save_str
                save_results(
                    filename,
                    split_grammar_results,
                    save_nulls,
                    list(grammar.nodes),
                    save_str)

        else:
            null_sigs = []
            null_median_medians = []
            save_nulls = []
            for null_idx in range(len(null_results_all)):
                null_results = null_results_all[null_idx][key]
                null_medians = [
                    np.median(vals)
                    for vals in null_results]
                pval = get_sig(grammar_median, null_medians)
                null_sigs.append(pval)
                # get the median null dist
                if len(null_medians) % 2 == 0:
                    null_medians = null_medians[1:]
                null_median_median = np.median(null_medians)
                null_median_medians.append(null_median_median)
                save_idx = null_medians.index(null_median_median)
                save_nulls.append(null_results[save_idx])

            #print grammar_median, null_median_medians, null_sigs
            filename = "{}.txt".format(file_prefix)
            save_str = "filename={};pwms={};grammar_avg={};null_avgs={};null_sigs={}".format(
                filename,
                grammar.nodes,
                grammar_median,
                null_median_medians,
                null_sigs)
            print save_str
            save_results(
                filename,
                grammar_results,
                save_nulls,
                list(grammar.nodes),
                save_str)

    return None

if __name__ == "__main__":
    main()

max_stars_repo_name: kundajelab/tronn

max_issues_repo_name: kundajelab/tronn

max_forks_repo_name: kundajelab/tronn

reason: import numpy,from numpy,import networkx


text:
__copyright__ = "Copyright (c) 2021 Jina AI Limited. All rights reserved."
__license__ = "Apache-2.0"

import os

import numpy as np
from PIL import Image

from .. import ImageCropper

cur_dir = os.path.dirname(os.path.abspath(__file__))

def create_random_img_array(img_height, img_width):
    return np.random.randint(0, 256, (img_height, img_width, 3))

def test_crop():
    img_size = 217
    img_array = create_random_img_array(img_size, img_size)
    left = 2
    top = 17
    width = 30
    height = 20
    crafter = ImageCropper(top=top, left=left, width=width, height=height)
    crafted_docs = crafter.craft(np.stack([img_array, img_array]))

    assert len(crafted_docs) == 2
    for crafted_doc in crafted_docs:
        np.testing.assert_array_equal(
            crafted_doc['blob'], np.asarray(img_array[top:top + height, left:left + width, :]),
            'img_array: {}\ntest: {}\ncontrol:{}'.format(
                img_array.shape,
                crafted_doc['blob'].shape,
                np.asarray(img_array[left:left + width, top:top + height, :]).shape))
        crop = Image.fromarray(np.uint8(crafted_doc['blob']))
        crop_width, crop_height = crop.size
        assert crop_width == width
        assert crop_height == height
        assert crafted_doc['location'] == (top, left)

def test_crop_file_image():
    tmp_fn = os.path.join(cur_dir, 'imgs/cars.jpg')
    img = Image.open(tmp_fn).convert('RGB')
    img_array = np.array(img).astype('float32')
    crafter = ImageCropper(top=541, left=992, width=24, height=67)
    crafted_docs = crafter.craft(np.stack([img_array, img_array]))
    assert len(crafted_docs) == 2
    for crafted_doc in crafted_docs:
        assert crafted_doc['blob'].shape == (67, 24, 3)
        crop_real_img = Image.open(os.path.join(cur_dir, 'imgs/faster_rcnn/person-0.png'))
        crop_real_img_array = np.array(crop_real_img).astype('float32')
        np.testing.assert_array_almost_equal(crafted_doc['blob'], crop_real_img_array)

ext: py

max_stars_repo_name: strawberrypie/jina-hub

max_issues_repo_name: strawberrypie/jina-hub

max_forks_repo_name: strawberrypie/jina-hub

reason: import numpy


text:
#!/usr/bin/env python3
import numpy as np
np.set_printoptions(linewidth=200)

# https://www.youtube.com/watch?v=cmy7LBaWuZ8&t=74s

D = 8

N = D
A = np.random.rand(N,N)
B = np.random.rand(N,N)
#A = np.arange(1,N*N+1).reshape(N,N)
#B = np.ones((N,N))*2
#B = np.arange(1,17).reshape(4,4) + 10
C = A @ B

def reset(D):
  global acc, acache, wcache
  acc = np.zeros((D,D))
  acache = np.zeros((D,D))
  wcache = np.zeros((D,D))

def mxu(a):
  global acc, acache, wcache
  assert a.shape == (D,)
  acache = np.concatenate([a[:, None], acache[:, :-1]], axis=1)
  ret = np.copy(acc[0])
  acc = np.concatenate([acc[1:], np.zeros((1, D), acc.dtype)], axis=0)
  acc += acache * wcache
  """
  print("****")
  print(acache)
  print(acc)
  print(ret)
  """
  return ret

def apad(a):
  ret = np.zeros((a.shape[0], a.shape[1]+(2*a.shape[0])-1)) # length of the flowing dim can be theoretically infinite
  for i in range(0, a.shape[0]):                            # A.shape = (X, D), B.shape = (D, D) -> C.shape = (X, D)
    ret[i, i+a.shape[0]:i+a.shape[1]+a.shape[0]] = a[i]     # padding needed for a.shape[0] steps on each side
  return ret

def unapad(a):
  ret = np.zeros(((a.shape[0] - a.shape[1] +1), a.shape[1]))
  for i in range(0, ret.shape[1]):
    ret[:, i] = a[a.shape[1]-i-1:ret.shape[0]+a.shape[1]-i-1, i]
  return ret

print(A)
print(B)
print(C)
print("**************************")
print("Weight Stationary Systolic Array") # static weights inside systolic array units

reset(D)
AA = apad(A.T)
wcache = B.copy()
out = []
for n in range(3*N-1): # since we padded on both sides we have N zeros on the end of the array to finish accumulation
  r = mxu(AA[:, -(n+1)])
  out.append(r)
ret = unapad(np.array(out)[N:][::-1])
print(ret)
assert np.allclose(C, ret)

def reset_mxm(N, X):
  global acc, acache, bcache
  acc = np.zeros((N,N))
  acache = np.zeros((N,N))
  bcache = np.zeros((N,N))

def mxm(a_col, b_row):
  global acc, acache, bcache

  # length of the flowing dim can be theoretically infinite
  # A.shape = (D, X) , B.shape = (X, D) -> C.shape = (D, D)
  # usually more useful than fixed weight matrix since you can chain matmuls easier that way
  # but supporting both ways in hardware shouldnt be too hard (just a simple "if" to decide what to pass on)

  acache = np.concatenate([a_col[:, None], acache[:, :-1]], axis=1) # flow the cols of A
  bcache = np.concatenate([b_row[None, :], bcache[:-1, :]], axis=0) # flow the cols of B
  acc += acache * bcache                                            # FMAC

def apad_mxm(a):
  ret = np.zeros((a.shape[0], a.shape[1]+(2*a.shape[0])-1))
  for i in range(0, a.shape[0]):
    ret[i, i:i+a.shape[1]] = a[i]
  return ret

print("**************************")
print("Output Stationary Systolic Array") # static accumulators inside systolic array units

reset_mxm(D, D)
AA = apad_mxm(A).T
BB = apad_mxm(B.T)

for n in range(3*N-1):
  mxm(AA[n, :], BB[:, n]) # no need for unpad since acc has results

ret2 = acc
print(ret2)
assert np.allclose(C, ret2)

print("")
print("**************************")
print("")
print("**************************")
print("Weight Stationary Systolic Array - Dynamic Flow Size")

D = 8
X = 16 # Free dimension, matmul efficiency goes to 1 as X->inf (X / (X + 2D))

N = D
A = np.random.rand(X,N)
B = np.random.rand(N,N)
C = A @ B
print(C)

reset(D)
AA = apad(A.T)
wcache = B.copy()
out = []
for n in range(2*N+X-1): # since we padded on both sides we have N zeros on the end of the array to finish accumulation
  r = mxu(AA[:, -(n+1)])
  out.append(r)
ret3 = unapad(np.array(out)[N:][::-1])
print("**************************")
print(ret3)
assert np.allclose(C, ret3)

print("")
print("**************************")
print("")
print("**************************")
print("Output Stationary Systolic Array - Dynamic Flow Size")

D = 8
X = 16 # Free dimension, matmul efficiency goes to 1 as X->inf (X / (X + 2D))

N = D
A = np.random.rand(N,X)
B = np.random.rand(X,N)
C = A @ B
print(C)

reset_mxm(N,X)
AA = apad_mxm(A).T
BB = apad_mxm(B.T)
for n in range(2*N+X-1):
  mxm(AA[n, :], BB[:, n]) # no need for unpad since acc has results

ret3 = acc

print("**************************")
print(ret3)
assert np.allclose(C, ret3)

print("")
print("**************************")
print("")
print("**************************")
print("Output Stationary Systolic Array - Dynamic Flow Size, broadcasted dim loading")

def bmxm(a_col, b_row): # broadcasing the inputs across dims on every clk makes padding unnecessary
  global acc, acache, bcache

  # length of the flowing dim can be theoretically infinite
  # A.shape = (D, X) , B.shape = (X, D) -> C.shape = (D, D)
  # this impl just broadcasts an input across the dimension that way we can dump the padding
  # we can support all 3 ways but it might be a problem with regard to wire density
  acache = np.repeat(a_col[:, None], len(a_col), axis=1) # flow the cols of A
  bcache = np.repeat(b_row[None, :], len(b_row), axis=0) # flow the rows of B
  acc += acache * bcache # FMAC

D = 8
X = 16 # Free dimension

N = D
A = np.random.rand(N,X)
B = np.random.rand(X,N)
C = A @ B
print(C)

reset_mxm(N,X)
AA = A
BB = B
for x in range(X):
  bmxm(AA[:, x], BB[x, :]) # no need for unpad since acc has results
ret3 = acc

print("**************************")
print(ret3)
assert np.allclose(C, ret3)

max_stars_repo_name: EladDv/tinygrad

max_issues_repo_name: EladDv/tinygrad

max_forks_repo_name: EladDv/tinygrad


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import rospy
import tf
import tf2_ros
from std_msgs.msg import String
from std_srvs.srv import Trigger,TriggerResponse,SetBool,SetBoolResponse
from geometry_msgs.msg import Twist,TransformStamped
from vapaa_cruiser.msg import NavState
from vapaa_cruiser.srv import TriggerWithInfo,TriggerWithInfoResponse
import math
import numpy as np
import json
import sqlite3
import rospkg

class AutoNavigation():
    def __init__(self):
        self.updateRate = rospy.get_param("~updateRate",30)
        rospack = rospkg.RosPack()
        self.dbFile = rospy.get_param("~dbFile",rospack.get_path("vapaa_cruiser")+"/../../../web/vapaa_cruiser.db")
        self.useFakeGPS = rospy.get_param("~useFakeGPS", False)
        self.fakeInitLat = rospy.get_param("~fakeInitLat", 23)
        self.fakeInitLng = rospy.get_param("~fakeInitLng", 121)
        self.fakeInitAngle = rospy.get_param("~fakeInitAngle", 180)

        self.subCarState = rospy.Subscriber("car_state",String,self.UpdateState)
        self.pubNavState = rospy.Publisher("nav_state",NavState,queue_size=1)
        self.pubCarCmd = rospy.Publisher("car_cmd", Twist, queue_size=1)

        self.navState = "stop"
        self.navPause = False
        self.navLoop = False
        self.curPath = None
        self.targetIndex = 0
        self.drive = {
            "targetForward": 0, "targetTurn": 0,
            "curForward": 0, "curTurn": 0,
            "errForward": 0, "errTurn": 0,
            "errSumForward": 0, "errSumTurn": 0,
            "forwardP": 1, "forwardI": 0, "forwardD": 0,
            "turnP": 1, "turnI": 0, "turnD": 0,
            "adjustByUS": True
        }
        self.lat = -9999
        self.lng = -9999
        self.usDist = {"LF":-1,"F":-1,"RF":-1,"LB":-1,"B":-1,"RB":-1}
        self.Trans = {  #for odom to gps transform computation
            #y依angle旋轉後乘scale，加上offsetY即為lat
            #x依angle旋轉後乘scale，除以cos(lat)調整，再加上offsetX即為lng
            "maxNum": 100,
            "mappingArr": [],   #存原始資料
            "dxSum": 0,
            "dySum": 0,
            "dLatSum": 0,
            "dLngSum": 0,
            "angle": 0,
            "scale": 0,
            "offsetX": 0,
            "offsetY": 0,
            "inited": False
        }
        self.srvStartNav = rospy.Service("autoNavigation/start", TriggerWithInfo, self.ServiceStartNav)
        self.srvStopNav = rospy.Service("autoNavigation/stop", Trigger, self.ServiceStopNav)
        self.srvPauseNav = rospy.Service("autoNavigation/pause", SetBool, self.ServiceSetPause)
        self.srvSetLoop = rospy.Service("autoNavigation/setLoop", SetBool, self.ServiceSetLoop)

        self.tfBuffer = tf2_ros.Buffer()
        self.tfListener = tf2_ros.TransformListener(self.tfBuffer)

        #realsense向左看，對z軸旋轉-90度才是車體正前方
        #tf tree: /map -> /odom -> /camera_link -> /car
        broadcaster = tf2_ros.StaticTransformBroadcaster()
        t = TransformStamped()
        t.header.stamp = rospy.Time.now()
        t.header.frame_id = "camera_link"
        t.child_frame_id = "car"
        t.transform.translation.x = 0
        t.transform.translation.y = 0
        t.transform.translation.z = 0
        quat = tf.transformations.quaternion_from_euler(0,0,-math.pi*0.5)
        t.transform.rotation.x = quat[0]
        t.transform.rotation.y = quat[1]
        t.transform.rotation.z = quat[2]
        t.transform.rotation.w = quat[3]
        broadcaster.sendTransform(t)

    def UpdateState(self,msg):
        data = msg.data.split(",")
        self.lat = float(data[0])
        self.lng = float(data[1])
        self.usDist["FL"] = float(data[2])
        self.usDist["F"] = float(data[3])
        self.usDist["FR"] = float(data[4])
        self.usDist["BL"] = float(data[5])
        self.usDist["B"] = float(data[6])
        self.usDist["BR"] = float(data[7])

    def ServiceStartNav(self,request):
        info = json.loads(request.info)
        conn = sqlite3.connect(self.dbFile)
        c = conn.cursor()
        cmd = "SELECT * FROM NavigationPath WHERE id='%s';" % (info["id"])
        row = c.execute(cmd).fetchone()
        conn.close()
        if row is None:
            return TriggerWithInfoResponse(
                success=False,
                message="path not exist"
            )
        self.curPath = {
            "id": row[0],
            "userID": row[1],
            "path": json.loads(row[2])
        }
        self.targetIndex = 0
        self.navState = "start"
        self.navPause = True
        return TriggerWithInfoResponse(
            success=True,
            message=self.curPath["id"]
        )

    def ServiceStopNav(self,request):
        self.navState = "stop"
        self.curPath = None
        self.StopDrive()
        return TriggerResponse(
            success=True,
            message=""
        )

    def ServiceSetPause(self,request):
        self.navPause = request.data
        if not self.navPause and self.curPath is not None:
            target = self.curPath["path"]["ptArr"][self.targetIndex]
            print("continue go to target %d, lat: %f, lng: %f" % (self.targetIndex,target["lat"],target["lng"]))
        else:
            self.StopDrive()

        return SetBoolResponse(
            success=True,
            message=""
        )

    def ServiceSetLoop(self,request):
        self.navLoop = request.data
        return SetBoolResponse(
            success=True,
            message=""
        )

    def UpdateOdomToGPSTransform(self,x,y,lat,lng):
        #假設地球是正圓形，經度隨距離變化量隨緯度越高變化越快
        #將經度變化量乘以cos(緯度)消除緯度的影響
        #此時經度、緯度變化360度的距離 = 2pi*地球半徑
        r = 6378137 #地球半徑(m)
        self.Trans["scale"] = 180.0/(r*math.pi)

        dSize = len(self.Trans["mappingArr"])
        if dSize < 2:
            self.Trans["mappingArr"].append({
                "x": x, "y": y, "lat": lat, "lng": lng,
                "dx": 0, "dy": 0, "dLat":0, "dLng":0
            })
        else:  #compute difference
            newPt = {
                "x": x, "y":y, "lat":lat, "lng":lng
            }
            lastPt = self.Trans["mappingArr"][dSize-1]
            newPt["dx"] = newPt["x"] - lastPt["x"]
            newPt["dy"] = newPt["y"] - lastPt["y"]
            if newPt["dx"]*newPt["dx"]+newPt["dy"]*newPt["dy"] < 0.01:  #距離太近，不加入計算
                return
            newPt["dLat"] = (newPt["lat"] - lastPt["lat"])
            newPt["dLng"] = (newPt["lng"] - lastPt["lng"])
            newPt["dLng"] *= math.cos(lastPt["lat"]*math.pi/180.0)
            self.Trans["mappingArr"].append(newPt)
            self.Trans["dxSum"] += newPt["dx"]
            self.Trans["dySum"] += newPt["dy"]
            self.Trans["dLatSum"] += newPt["dLat"]
            self.Trans["dLngSum"] += newPt["dLng"]

        if len(self.Trans["mappingArr"]) > self.Trans["maxNum"]:    #moving window
            firstPt = self.Trans["mappingArr"][0]
            self.Trans["dxSum"] -= firstPt["dx"]
            self.Trans["dySum"] -= firstPt["dy"]
            self.Trans["dLatSum"] -= firstPt["dLat"]
            self.Trans["dLngSum"] -= firstPt["dLng"]
            self.Trans["mappingArr"].pop(0)
        dSize = len(self.Trans["mappingArr"])

        #compute angle between xy & latlng
        xyNorm = math.sqrt(self.Trans["dxSum"]*self.Trans["dxSum"]+self.Trans["dySum"]*self.Trans["dySum"])
        if xyNorm == 0:
            return
        latLngNorm = math.sqrt(self.Trans["dLatSum"]*self.Trans["dLatSum"]+self.Trans["dLngSum"]*self.Trans["dLngSum"])
        if latLngNorm == 0:
            return
        xyVec = [self.Trans["dxSum"]/xyNorm, self.Trans["dySum"]/xyNorm]
        latLngVec = [self.Trans["dLngSum"]/latLngNorm, self.Trans["dLatSum"]/latLngNorm]
        dot = xyVec[0]*latLngVec[0]+xyVec[1]*latLngVec[1]
        cross = xyVec[0]*latLngVec[1]-xyVec[1]*latLngVec[0]
        self.Trans["angle"] = math.atan2(cross,dot)

        #compute offset between xy & latlng
        cosAngle = dot
        sinAngle = cross
        offsetXSum = 0
        offsetYSum = 0
        for pt in self.Trans["mappingArr"]:
            x = self.Trans["scale"]*(cosAngle*pt["x"]-sinAngle*pt["y"])
            y = self.Trans["scale"]*(sinAngle*pt["x"]+cosAngle*pt["y"])
            cosLat = math.cos(pt["lat"]*math.pi/180.0)
            if cosLat == 0: #極地
                continue
            x = x/cosLat
            offsetXSum += (pt["lng"]-x)
            offsetYSum += (pt["lat"]-y)
        self.Trans["offsetX"] = offsetXSum/dSize
        self.Trans["offsetY"] = offsetYSum/dSize

        #print([self.Trans["scale"],self.Trans["angle"],self.Trans["offsetX"],self.Trans["offsetY"]])

    def XYToLatLng(self,x,y):
        cosAngle = math.cos(self.Trans["angle"])
        sinAngle = math.sin(self.Trans["angle"])
        lat = self.Trans["scale"]*(sinAngle*x+cosAngle*y)+self.Trans["offsetY"]
        lng = self.Trans["scale"]*(cosAngle*x-sinAngle*y)/math.cos(lat*math.pi/180.0)+self.Trans["offsetX"]
        return {"lat":lat, "lng":lng}

    def LatLngToXY(self,lat,lng):
        y = (lat-self.Trans["offsetY"])/self.Trans["scale"]
        x = (lng-self.Trans["offsetX"])/self.Trans["scale"]
        cosAngle = math.cos(-self.Trans["angle"])
        sinAngle = math.sin(-self.Trans["angle"])
        rotX = cosAngle*x-sinAngle*y
        rotY = sinAngle*x+cosAngle*y
        return {"x":rotX, "y":rotY}

    def GenFakeLatLng(self,x,y): #generate fake lat lng for testing
        offsetX = self.fakeInitLng
        offsetY = self.fakeInitLat
        angle = self.fakeInitAngle*math.pi/180.0
        scale = 180.0/(6378137*math.pi)
        #print([scale,angle,offsetX,offsetY])

        cosAngle = math.cos(angle)
        sinAngle = math.sin(angle)
        self.lat = scale*(sinAngle*x+cosAngle*y)+offsetY
        #直接assign給self.lng python會莫名其妙出現typeerror: can't multiply sequence by non-int of type 'float'
        lat = scale*(cosAngle*x-sinAngle*y)/math.cos(self.lat*math.pi/180.0)+offsetX
        self.lng = lat
        #print([x,y,self.lat,self.lng])

    def CheckGPSValid(self,lat,lng):
        lat = float(lat)
        lng = float(lng)
        if math.isnan(lat) or math.isnan(lng):
            return False
        if lat < -90 or lat > 90 or lng < -180 or lng > 180:
            return False
        return True

    def StopDrive(self):
        msg = Twist()
        msg.linear.x = 0
        msg.angular.z = 0
        self.pubCarCmd.publish(msg)

    def AutoDrive(self, curPose):
        if self.navState != "start":
            return
        if self.navPause:
            return
        if self.curPath is None:
            return
        if self.targetIndex < 0 or self.targetIndex >= len(self.curPath["path"]["ptArr"]):
            return

        distTolerance = 0.2
        turnScale = 1.0
        forwardScale = 0.3
        target = self.curPath["path"]["ptArr"][self.targetIndex]
        targetXY = self.LatLngToXY(target["lat"],target["lng"])
        diffXY = {
            "x": targetXY["x"]-curPose.translation.x,
            "y": targetXY["y"]-curPose.translation.y
        }
        #check if arrive target
        dist = math.sqrt(diffXY["x"]*diffXY["x"]+diffXY["y"]*diffXY["y"])
        if dist < distTolerance:
            self.targetIndex+=1
            if self.targetIndex >= len(self.curPath["path"]["ptArr"]):
                if self.navLoop:
                    self.targetIndex = 0
                else:
                    self.navState = "stop"
                    self.curPath = None
                    print("finish navigation")
                    return

            target = self.curPath["path"]["ptArr"][self.targetIndex]
            print("go to target %d, lat: %f, lng: %f" % (self.targetIndex,target["lat"],target["lng"]))

        #compute drive dirrection
        q = [curPose.rotation.x,curPose.rotation.y,curPose.rotation.z,curPose.rotation.w]
        curAngle = tf.transformations.euler_from_quaternion(q)[2]
        targetAngle = math.atan2(diffXY["y"],diffXY["x"])
        self.drive["targetTurn"] = -(targetAngle-curAngle)*turnScale
        self.drive["targetForward"] = dist*forwardScale
        #print([self.drive["targetForward"],self.drive["targetTurn"]])

        #依超音波測距調整command
        if self.drive["adjustByUS"]:
            turnForce = 0
            turnExpScale = 2e-3
            breakFactor = 1
            breakExpScale = 4e-3
            activeDist = 1000
            if self.drive["targetForward"] > 0:
                if self.usDist["FL"] > 0 and self.usDist["FL"] < activeDist:
                    turnForce += math.exp(-self.usDist["FL"]*turnExpScale)
                if self.usDist["FR"] > 0 and self.usDist["FR"] < activeDist:
                    turnForce -= math.exp(-self.usDist["FR"]*turnExpScale)
                if self.usDist["F"] > 0 and self.usDist["F"] < activeDist:
                    breakFactor = (1-math.exp(-self.usDist["F"]*breakExpScale))
            else:
                if self.usDist["BL"] > 0 and self.usDist["BL"] < activeDist:
                    turnForce -= math.exp(-self.usDist["BL"]*turnExpScale)
                if self.usDist["BR"] > 0 and self.usDist["BR"] < activeDist:
                    turnForce += math.exp(-self.usDist["BR"]*turnExpScale)
                if self.usDist["B"] > 0 and self.usDist["B"] < activeDist:
                    breakFactor = (1-math.exp(-self.usDist["B"]*breakExpScale))
            #print([self.drive["targetForward"],self.drive["targetTurn"],breakFactor,turnForce])
            self.drive["targetTurn"] += turnForce
            self.drive["targetForward"] *= breakFactor

        #limit target magnitude
        if self.drive["targetForward"] > 1:
            self.drive["targetForward"] = 1
        elif self.drive["targetForward"] < -1:
            self.drive["targetForward"] = -1
        if self.drive["targetTurn"] > 1:
            self.drive["targetTurn"] = 1
        elif self.drive["targetTurn"] < -1:
            self.drive["targetTurn"] = -1
        print([self.drive["targetForward"],self.drive["targetTurn"]])

        #update command by PID controller
        errForward = self.drive["targetForward"] - self.drive["curForward"]
        errTurn = self.drive["targetTurn"] - self.drive["curTurn"]
        errDiffForward = errForward - self.drive["errForward"]
        errDiffTurn = errTurn - self.drive["errTurn"]
        self.drive["errSumForward"] += errForward
        self.drive["errSumTurn"] += errTurn

        self.drive["curForward"] += errForward*self.drive["forwardP"] + self.drive["errSumForward"]*self.drive["forwardI"] + errDiffForward*self.drive["forwardD"]

        self.drive["curTurn"] += errTurn*self.drive["turnP"] + self.drive["errSumTurn"]*self.drive["turnI"] + errDiffTurn*self.drive["turnD"]

        self.drive["errForward"] = errForward
        self.drive["errTurn"] = errTurn

        #limit command
        if self.drive["curForward"] > 1:
            self.drive["curForward"] = 1
        elif self.drive["curForward"] < -1:
            self.drive["curForward"] = -1
        if self.drive["curTurn"] > 1:
            self.drive["curTurn"] = 1
        elif self.drive["curTurn"] < -1:
            self.drive["curTurn"] = -1
        #print([self.drive["curForward"],self.drive["curTurn"]])

        #publish command message
        msg = Twist()
        msg.linear.x = self.drive["curForward"]
        msg.angular.z = self.drive["curTurn"]
        self.pubCarCmd.publish(msg)

    def Run(self):
        rate = rospy.Rate(self.updateRate)
        start = rospy.Time.now()
        while not rospy.is_shutdown():
            try:
                #get car pose from odometry
                pose = self.tfBuffer.lookup_transform("map", "car", rospy.Time()).transform
                elapse = (rospy.Time.now()-start)
                t = elapse.secs+elapse.nsecs*1e-9
                #pose.translation.x = t
                #pose.translation.y = t
                if self.useFakeGPS:
                    self.GenFakeLatLng(pose.translation.x,pose.translation.y)
                if not self.CheckGPSValid(self.lat,self.lng):
                    print("invalid gps")
                    continue
                if not self.Trans["inited"]:
                    self.Trans["angle"] = 0
                    self.Trans["offsetX"] = self.lng
                    self.Trans["offsetY"] = self.lat
                    self.Trans["inited"] = True
                self.UpdateOdomToGPSTransform(pose.translation.x,pose.translation.y,self.lat,self.lng)

                self.AutoDrive(pose)

                #publish navigation state
                latlng = self.XYToLatLng(pose.translation.x,pose.translation.y)
                q = [pose.rotation.x,pose.rotation.y,pose.rotation.z,pose.rotation.w]
                angle = tf.transformations.euler_from_quaternion(q)
                #print([pose.translation.x,pose.translation.y,(angle[2]+self.Trans["angle"])*180/math.pi])
                navStateMsg = NavState()
                navStateMsg.lat = latlng["lat"]
                navStateMsg.lng = latlng["lng"]
                navStateMsg.angle = (angle[2]+self.Trans["angle"])*180/math.pi
                navStateMsg.state = self.navState
                navStateMsg.loop = self.navLoop
                navStateMsg.pause = self.navPause
                if self.curPath is not None:
                    navStateMsg.pathID = self.curPath["id"]
                    navStateMsg.targetIndex = self.targetIndex
                else:
                    navStateMsg.pathID = ""
                    navStateMsg.targetIndex = 0
                self.pubNavState.publish(navStateMsg)

            except (tf2_ros.LookupException, tf2_ros.ConnectivityException, tf2_ros.ExtrapolationException):
                #rospy.logwarn("auto_navigation lookup transform error")
                pass
            rate.sleep()

if __name__ == '__main__':
    rospy.init_node('auto_navigation_node')
    rospy.loginfo("auto_navigation_node started")
    an = AutoNavigation()
    an.Run()


text:
import os
import operator
import collections
import itertools

import numpy as np
from scipy.ndimage import gaussian_filter

from skimage import io, color, feature, morphology, filters, exposure, restoration, segmentation

import matplotlib.pyplot as plt

def main():
    pass

class ContactAngle(object):

    def __init__(self, img):
        """Initialize object with image.

        Args:
            img (np.ndarray, str): Image as array or /path/to/file
        """
        # if os.path.isfile(img) and isinstance(str, img):
        #     self.img = io.imread(img, as_grey=True)
        if False:
            pass
        elif isinstance(img, np.ndarray):
            if len(img.shape) > 2:
                self.img = img
                # self.img = color.rgb2gray(img)
            else:
                self.img = img.copy()
        else:
            raise ValueError("Must pass image (as numpy array) or /path/to/file.")

        self.img_processed = None
        self.labels = None
        self.measured_parameters = None

    def crop(self):
        """Crop image.

        Returns:
            None
        """
        shape_x, shape_y = self.img.shape[0], self.img.shape[1]
        x_lo, x_hi = int(shape_x * .2), int(shape_x * .8)
        y_lo, y_hi = int(shape_y * .2), int(shape_y * .8)
        img2 = self.img[x_lo: x_hi, y_lo: y_hi]
        self.img = img2.copy()

    def show(self, title='Processed'):
        """Plot image (matplotlib).

        args:
            title (str): title for image

        Returns:
            None
        """
        if self.img_processed is None:
            f, a = plt.subplots(figsize=(18, 6))
            a.imshow(self.img, cmap="Greys_r")
        else:
            f, a = plt.subplots(ncols=2, figsize=(18, 6))
            a[0].imshow(self.img, cmap="gray")
            a[0].set_title('Original')

            a[1].imshow(self.img, cmap="gray")
            # a[1].imshow(self.img_processed, cmap="gray")
            if self.labels is not None:
                indices = np.argwhere(self.labels > 0)
                a[1].scatter(indices[:, 1], indices[:, 0], c='r', s=1)
            if self.measured_parameters is not None:
                d = self.measured_parameters
                a[1].scatter(*d['p1'], c='b', s=50, edgecolor='w')
                a[1].scatter(*d['p2'], c='g', s=50, edgecolor='w')
                a[1].scatter(*d['p1_b'], c='b', s=50, edgecolor='w')
                a[1].scatter(*d['p2_b'], c='g', s=50, edgecolor='w')
                a[1].scatter(*d['midpt'], c='y', s=50, edgecolor='w')
                a[1].scatter(*d['top_pt'], c='y', s=50, edgecolor='w')
                a[1].plot([d['p1'][0], d['p2'][0]], [d['p1'][1], d['p2'][1]], c='w')
                a[1].plot([d['midpt'][0], d['top_pt'][0]], [d['midpt'][1], d['top_pt'][1]], c='w')
                text = 'Left angle: {:.3f}$^o$\nRight angle: {:.3f}$^o$\nBase: {:.3f}\nHeight: {:.3f}'\
                    .format(d['Left contact angle'], d['Right contact angle'], d['Width'], d['Height'])
                a[1].text(50, 200, text, bbox={'facecolor': 'white', 'alpha': 1, 'pad': 5}, fontdict={'size': 12})
            a[1].set_title(title)
            f.tight_layout()

    def get_input_image(self):
        """Decide which image should be processed (default or a 'workspace' image).

        Returns:
            np.ndarray: image to process.
        """
        input_img = self.img if self.img_processed is None else self.img_processed
        return input_img

    def canny_edge(self, sigma=1):
        """Canny edge detection on image.

        Args:
            sigma (float):  parameter for canny edge (from skimage)

        Returns:
            None
        """
        input_img = self.get_input_image()

        self.img_processed = feature.canny(input_img, sigma=sigma)

    def roberts_edge(self):
        """Perform roberts edge detection.

        Returns:
            None
        """
        input_img = self.get_input_image()

        self.img_processed = filters.roberts(input_img)

    def dilation(self, structure=np.ones((3, 3))):
        """Perform dilation on image.

        Args:
            structure (np.ndarray): structure for dilation.

        Returns:
            None
        """
        input_img = self.get_input_image()

        self.img_processed = morphology.dilation(input_img, structure)

    def closing(self, structure=np.ones((3, 3))):
        """Perform closing.

        Args:
            structure (np.ndarray):  structure for closing

        Returns:
            None
        """
        input_img = self.get_input_image()

        self.img_processed = morphology.closing(input_img, structure)

    def label_regions(self):
        """Find continously connected regions.

        Returns:
            None
        """
        if self.img_processed is None:
            raise RuntimeError("Labelling raw image isn't allowed.")
        input_img = self.img_processed
        labels = morphology.label(input_img)
        region_sizes = collections.Counter(labels.flatten())
        region_sizes = {key: value for key, value in region_sizes.items() if key > 0}
        largest_region = max(region_sizes.items(), key=operator.itemgetter(1))[0]
        labels[labels != largest_region] = 0
        labels[labels == largest_region] = 1
        self.labels = labels

    def furthest_two_points(self):
        """Find furthest two points in an array.  Points to consider must have values = 1.  All other elements
        should be 0.

        Returns:
            tuple: coordinates of two points, ordered by x-value
        """
        furthest = -1
        furthest_i, furthest_j = None, None
        contour = np.argwhere(self.labels >= 1)
        low_pct, high_pct = np.percentile(contour[:, 1], 5), np.percentile(contour[:, 1], 95)
        contour_low = contour[contour[:, 1] < low_pct]
        contour_high = contour[contour[:, 1] > high_pct]
        for x in itertools.product(contour_low, contour_high):
            i, j = x
            dist = np.linalg.norm(i - j)
            if dist > furthest:
                furthest = dist
                furthest_i = i
                furthest_j = j

        p1 = (furthest_i[1], furthest_i[0])
        p2 = (furthest_j[1], furthest_j[0])

        if p1[0] < p2[0]:
            return p1, p2
        else:
            return p2, p1

    def calc_midpt(self, p1, p2):
        """Calculate midpoint between two points.

        Args:
            p1 (tuple): (x, y) of first point
            p2 (tuple): (x, y) of second point

        Returns:
            tuple: (x, y) of midpoint
        """
        midpt = [(p1[0] + p2[0]) // 2, (p1[1] + p2[1]) //2]
        return midpt

    def equalize_hist(self):
        input_img = self.get_input_image()

        self.img_processed = exposure.equalize_adapthist(input_img)

    def erosion(self, structure=np.ones((3, 3))):
        input_img = self.get_input_image()

        self.img_processed = morphology.erosion(input_img, selem=structure)

    def adjust_gamma(self, gamma=1):
        input_img = self.get_input_image()

        self.img_processed = exposure.adjust_gamma(input_img, gamma=gamma)

    def adjust_histogram(self, gamma=1):
        input_img = self.get_input_image()

        self.img_processed = exposure.equalize_hist(input_img, nbins=100)

    def extract(self):
        """Extract contact angles, height, and width of droplet in image

        Returns:
            dict: measured parameters

            keys of most interest are:
                'Left contact angle', 'Right contact angle', 'Width', 'Height'

            other keys used to derive these parameters are:
                'p1', 'p1_b' -- used to calculate left contact angle
                'p2', 'p2_b' -- used to calculate right contact angle
                'midpt', 'top_pt' -- used to calculate contact angles, width, and height
        """
        p1, p2 = self.furthest_two_points()
        midpt = self.calc_midpt(p1, p2)

        p1_by = p1[1] - 20
        p1_bx = np.argmax(self.labels[p1_by, :])
        p1_b = (p1_bx, p1_by)

        p2_by = p2[1] - 20
        p2_bx = self.labels.shape[1] - np.argmax(self.labels[p2_by, :][::-1])
        p2_b = (p2_bx, p2_by)

        p1_angle = self.angle_between_lines(p1, midpt, p1_b)
        p2_angle = self.angle_between_lines(p2, midpt, p2_b)

        top_pt = self.find_90_deg_point(midpt, p2)
        height = self.calc_dist(top_pt, midpt)
        length = self.calc_dist(p1, p2)

        self.measured_parameters = {'p1': p1, 'p1_b': p1_b, 'Left contact angle': p1_angle,
                                    'p2': p2, 'p2_b': p2_b, 'Right contact angle': p2_angle,
                                    'midpt': midpt, 'top_pt': top_pt, 'Height': height, 'Width': length}

        return self.measured_parameters

    def angle_between_lines(self, p0, p1, p2):
        """Calculate angle between three points.

        Args:
            p0 (tuple): (x, y) of vertex
            p1 (tuple): (x, y) of non-vertex
            p2 (tuple): (x, y) of non-vertex

        Returns:
            float: angle (in degrees)
        """
        p0 = np.asarray(p0)
        p1 = np.asarray(p1)
        p2 = np.asarray(p2)
        p1_p0 = p1 - p0
        p2_p0 = p2 - p0

        cos_theta = np.dot(p1_p0, p2_p0) / (np.linalg.norm(p1_p0) * np.linalg.norm(p2_p0))
        return np.degrees(np.arccos(cos_theta))

    def find_90_deg_point(self, p0, p1):
        """Find point closest to 90 degrees from two other points.

        Args:
            p0 (tuple): (x, y) of vertex
            p1 (tuple): (x, y) of non-vertex

        Returns:
            tuple: (x, y) of point
        """
        best_angle = 1e6
        best_p = None
        values = np.argwhere(self.labels >= 1)
        for p2 in values:
            p2 = p2[::-1]
            angle = self.angle_between_lines(p0, p1, p2)
            if np.abs(angle - 90) < best_angle and p2[1] < p1[1]:
                best_angle = np.abs(angle - 90)
                best_p = p2
        return best_p

    def calc_dist(self, p1, p2):
        """Calculate distance between two points.

        Args:
            p1 (tuple): (x, y) of point 1
            p2 (tuple): (x, y) of point 2

        Returns:
            float: distance
        """
        return np.linalg.norm(np.asarray(p1) - np.asarray(p2), 2)

    def run_pipeline(self, sigma=1):
        """Run 'current' contact angle extraction pipeline.

        Args:
            sigma (float):  sigma value for canny edge detection

        Returns:
            dict: parameters of contact angle -- see extract() docs
        """
        self.crop()
        self.canny_edge(sigma=sigma)
        self.dilation(morphology.disk(3))
        self.label_regions()
        params = self.extract()
        self.show(title='Final')
        return params

    def _explore_func(self):
        """Function used for exploratory analysis.  There should be *NO* expectation that this function provides
        consistent functionality.

        Returns:
            None
        """
        # self.crop()
        # self.img_processed = self.img[:, :, 0]
        # self.img_processed = filters.rank.minimum(
        #     filters.rank.maximum(self.img, morphology.disk(5)), morphology.disk(5))
        # self.img_processed = filters.rank.enhance_contrast_percentile(self.img, morphology.disk(5), p0=.2, p1=.5)
        # self.img_processed = filters.rank.enhance_contrast(self.img, morphology.disk(3))
        # self.img_processed = filters.rank.mean_bilateral(self.img, np.ones((3, 3)))
        # self.img_processed = filters.rank.mean_bilateral(self.img, morphology.disk(5))
        # self.img_processed = filters.rank.mean_bilateral(self.img, morphology.disk(10))
        # self.img_processed = filters.rank.(self.img, morphology.disk(10))
        # self.img_processed = exposure.equalize_hist(self.img, nbins=128)
        # self.img_processed = filters.threshold_otsu(self.img)
        # self.img_processed = filters.rank.equalize(self.img, morphology.disk(50))

        # image = gaussian_filter(self.img, .1)
        # image = self.img.copy()

        # seed = np.copy(image)
        # seed[1:-1, 1:-1] = image.min()
        # mask = image

        # self.img_processed = morphology.reconstruction(seed, mask, method='dilation', selem=morphology.disk(3))

        self.canny_edge(sigma=5)
        self.dilation(np.ones((10, 10)))
        # self.dilation(np.ones((3, 3)))
        # self.show()
        # self.dilation(morphology.disk(140))
        # self.show()
        # self.dilation()
        f, a = plt.subplots(ncols=2, figsize=(18, 6))
        a[0].imshow(self.img, cmap='gray')
        a[1].imshow(self.img_processed, cmap='gray')
        # self.label_regions()
        # self.extract()
        # self.show()

if __name__ == '__main__':
    main()

max_stars_repo_name: hackingmaterials/duramat

max_issues_repo_name: hackingmaterials/duramat

max_forks_repo_name: hackingmaterials/duramat

reason: import numpy,from scipy


text:
from cvxpy import Variable, quad_form, Parameter, Problem, sum, Minimize
import numpy as np
from posdef import nearestPD

def MarkowitzOpt(mean, variance, covariance, interest_rate, min_return):
    n = mean.size + 1                   # Number of assets (number of stocks + interest rate)

    mu = mean.values                    # Mean returns of n assets
    temp = np.full(n, interest_rate)
    temp[:-1] = mu
    mu = temp

    counter = 0
    Sigma = np.zeros((n,n))                 # Covariance of n assets
    for i in np.arange(n-1):
        for j in np.arange(i, n-1):
            if i==j:
                Sigma[i,j] = variance[i]
            else:
                Sigma[i,j] = covariance[counter]
                Sigma[j,i] = Sigma[i,j]
                counter+=1
    Sigma = nearestPD(Sigma)                # Converting covariance to the nearest positive-definite matrix

    # Ensuring feasability of inequality contraint
    if mu.max() < min_return:
        min_return = interest_rate

    w = Variable(n)                         # Portfolio allocation vector
    ret = mu.T*	w
    risk = quad_form(w, Sigma)
    min_ret = Parameter(nonneg=True)
    min_ret.value = min_return
    prob = Problem(Minimize(risk),          # Restricting to long-only portfolio
                   [ret >= min_ret,
                   sum(w) == 1,
                   w >= 0])
    prob.solve()
    return w.value

reason: import numpy,from cvxpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import glob
from PIL import Image
from PSNN import model, layers, activation

# this get inputs arrays form images in directory
#
# @param directory str location of dir with images
# @param target np array of target (output) for images in this dir
# @return tuple of [0] -> images as inputs arrays [1] -> traget (output) array for same input
def getITArray(directory, target):
  imgs = glob.glob(directory + "/*.png")
  out = ([], [])
  for loc in imgs:
    img = Image.open(loc).convert('L')  # convert image to 8-bit grayscale
    WIDTH, HEIGHT = img.size
    data = list(img.getdata()) # convert image data to a list of integers
    out[0].append([data[offset:offset+WIDTH] for offset in range(0, WIDTH*HEIGHT, WIDTH)])
    out[1].append(target)

  return out

netmodel = model([
  layers.flatten(),
  layers.activation(activation.tanh()),
  layers.dense(30, activation.tanh()),
  layers.dense(20, activation.tanh()),
  layers.dense(10, activation.sigmoid())
])

netmodel.debug = True

netmodel.create(
  inputs=(28,28)
)

# get learn data from images
learn = {
  "inputs": [],
  "targets": []
}

# fill imputs and targets
for i in range(10):
  # one hot encoding
  tar = np.zeros((10))
  tar[i] = 1
  # get inputs arrays and targets arrays from images
  inputs, targets = getITArray("mnist/training/" + str(i), tar)
  learn['inputs'] += inputs
  learn['targets'] += targets

# convert lists to np array
learn['inputs'] = np.array(learn['inputs'])
learn['targets'] = np.array(learn['targets'])

# load images for test
test = np.array(getITArray("mnist/testing", None)[0])

# start learning
print("learning form " + str(len(learn['inputs'])) + " samples")

netmodel.fit(
  rate=0.1,
  inputs=learn['inputs'],
  targets=learn['targets'],
  epochs=100,
  type='backPropagation'
)

print("try predict")
print("\n")
for sample in test:
  # print image
  print('\n'.join([''.join(['{:4}'.format(item) for item in row])
        for row in sample]))
  print("\n")
  # predict with image as inputs
  predict = netmodel.predict(sample)
  # print prediction
  print("it is " + str(np.argmax(predict)) + " with score " + str(np.max(predict)))
  print("\n\n")


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
This script saves bid and ask data for specified ETFs to files for each day
during market open hours.

It assumes the computer is at US East Coast Time.

@author: mark
"""

import os
import sys
import time

import pandas as pd
import numpy as np

import alpaca_trade_api as ati
from datetime import datetime

import unittest
from get_quotes_alpaca_polygon import get_time_averages

s = pd.Series({0:5, 1:6, 2:7})
averages = get_time_averages(s, 1, 0 , 4)
correct_averages = [5, 6, 7, 7]
assert 4 == len(averages)
for a, b in zip(averages.values, correct_averages):
    assert a == b

s = pd.Series({0:5, 0.5:6, 2.5:7, 3.5:8})
averages = get_time_averages(s, 1, 0 , 4)
correct_averages = [5.5, 6, 6.5, 7.5]
assert 4 == len(averages)
for a, b in zip(averages.values, correct_averages):
    assert a == b

s = pd.Series({0:5, 0.1:6, 1.75:7, 2:8})
averages = get_time_averages(s, 1, 0 , 4)
correct_averages = [5.9, 6.25, 8, 8]
assert 4 == len(averages)
for a, b in zip(averages.values, correct_averages):
    assert a == b

max_stars_repo_name: goldmanm/bid-ask-visualization

max_issues_repo_name: goldmanm/bid-ask-visualization

max_forks_repo_name: goldmanm/bid-ask-visualization


text:
"""
Validate our RNN. Basically just runs a validation generator on
about the same number of videos as we have in our test set.
"""

"""
Validate our RNN. Basically just runs a validation generator on
about the same number of videos as we have in our test set.
"""
import tensorflow as tf
from keras.backend.tensorflow_backend import set_session
config = tf.ConfigProto()
config.gpu_options.per_process_gpu_memory_fraction = 0.9
config.gpu_options.visible_device_list = "1"
config.allow_soft_placement=True
config.log_device_placement=False

set_session(tf.Session(config=config))

if 'session' in locals() and session is not None:
    print('Close interactive session')
    session.close()

from keras.callbacks import TensorBoard, ModelCheckpoint, CSVLogger
from model_caps_iv3 import ResearchModels
from data import DataSet
import functools
import keras.metrics

import numpy as np

from sklearn.metrics import classification_report, confusion_matrix

def validate(data_type, model, seq_length=60, saved_model=None,
             class_limit=None, image_shape=None):
    batch_size = 16

    # Get the data and process it.
    if image_shape is None:
        data = DataSet(
            seq_length=seq_length,
            class_limit=class_limit
        )
    else:
        data = DataSet(
            seq_length=seq_length,
            class_limit=class_limit,
            image_shape=image_shape
        )

    val_generator = data.frame_generator(batch_size, 'test', data_type)

    # Get the model.
    rm = ResearchModels(len(data.classes), model, seq_length, saved_model)

    # Evaluate!
    results = rm.model.evaluate_generator(
        generator=val_generator,
        val_samples=15,
        use_multiprocessing=True,
        workers=1)

    print(results)
    print(rm.model.metrics_names)

        #val_generator = data.frame_generator(1,'test', data_type)
    # Get the model.
    rm = ResearchModels(len(data.classes), model, seq_length, saved_model)
    # Evaluate!
    scores=np.zeros([14])
    total=np.zeros([14])

    val_trues=[]
    val_preds=[]
    for X,y in data.gen_test('test', data_type):
        results = rm.model.predict(X)
        predicted=np.argmax(results, axis=-1)
        idx=np.where(np.array(y)==1)
        true_label=idx[1]
        print(true_label)

        total[true_label]=total[true_label]+1

        print(len(predicted))
        print(len(true_label))
        if predicted[0]==true_label[0]:
            scores[true_label]=scores[true_label]+1

    #val_preds = np.argmax(results, axis=-1)

    print('Confusion Matrix')
    tn, fp, fn, tp =confusion_matrix(true_label, predicted).ravel()
    print("\n *** Confusion matrix**")
    print(tp)
    print(tn)
    print(fp)
    print(fn)
    print('\n****************')

    print(classification_report(true_label, predicted))
    print(scores)
    print('\n****************')
    print(total)

def main():
    model = 'lrcn'  #capsule,lstm,conv_3d,lrcn,c3d
    #saved_model = '/home/abdkhan/covid/data/checkpoints_caps/lstm-features.015-0.658-iv3-caps.hdf5'
    #saved_model = '/home/abdkhan/covid/data/checkpoints_caps/lstm-features.028-0.616-iv3-caps.hdf5'
    #saved_model = '/home/abdkhan/covid/data/checkpoints_caps/conv_3d-images.015-0.523-iv3-caps.hdf5'
    #saved_model = '/home/abdkhan/covid/data/checkpoints_caps/lrcn-images.009-2.551-iv3-caps.hdf5'
    #saved_model = '/home/abdkhan/covid/data/checkpoints_caps/c3d-images.004-0.507-iv3-caps.hdf5'

    ################### Inverted#####################
    saved_model = '/home/abdkhan/covid/data/checkpoints_caps/capsule-features.004-0.582-iv3-caps.hdf5'
    saved_model = '/home/abdkhan/covid/data/checkpoints_caps/lrcn-images.001-3.222-iv3-caps.hdf5'

        # Chose images or features and image shape based on network.
    if model in ['conv_3d', 'c3d', 'lrcn']:
        data_type = 'images'
        image_shape = (80, 80, 3)
    elif model in ['capsule','lstm']:
        data_type = 'features'
        image_shape = None
    else:
        raise ValueError("Invalid model. See train.py for options.")

    validate(data_type, model, saved_model=saved_model,
             image_shape=image_shape, class_limit=2)

if __name__ == '__main__':
    main()


text:
#
# minerals.py -- Some tables for ISM abundances and depletion factors
# that are useful for calculating dust mass and dust-to-gas ratios
#
# 2016.01.22 - lia@space.mit.edu
##----------------------------------------------------------------

import numpy as np

amu   = {'H':1.008,'He':4.0026,'C':12.011,'N':14.007,'O':15.999,'Ne':20.1797,
         'Na':22.989,'Mg':24.305,'Al':26.981,'Si':28.085,'P':30.973,'S':32.06,
         'Cl':35.45,'Ar':39.948,'Ca':40.078,'Ti':47.867,'Cr':51.9961,'Mn':54.938,
         'Fe':55.845,'Co':58.933,'Ni':58.6934}
amu_g = 1.661e-24  # g
mp    = 1.673e-24  # g (proton mass)

wilms = {'H':12.0, 'He':10.99, 'C':8.38, 'N':7.88, 'O':8.69, 'Ne':7.94,
         'Na':6.16, 'Mg':7.40, 'Al':6.33, 'Si':7.27, 'P':5.42, 'S':7.09,
         'Cl':5.12, 'Ar':6.41, 'Ca':6.20, 'Ti':4.81, 'Cr':5.51, 'Mn':5.34,
         'Fe':7.43, 'Co':4.92, 'Ni':6.05}  # 12 + log A_z

# Fraction of elements still in gas form
wilms_1mbeta = {'H':1.0, 'He':1.0, 'C':0.5, 'N':1.0, 'O':0.6, 'Ne':1.0, 'Na':0.25,
                'Mg':0.2, 'Al':0.02, 'Si':0.1, 'P':0.6, 'S':0.6, 'Cl':0.5, 'Ar':1.0,
                'Ca':0.003, 'Ti':0.002, 'Cr':0.03, 'Mn':0.07, 'Fe':0.3, 'Co':0.05,
                'Ni':0.04}

class Mineral(object):
    """
    | Use a dictionary to define the composition.
    | e.g. Olivines of pure MgFe^{2+}SiO_4 composition would be
    | olivine_halfMg = Mineral( {'Mg':1.0, 'Fe':1.0, 'Si':1.0, 'O':4.0} )
    | self.composition : dictionary containing elements and their weights
    |
    | @property
    | self._weight_amu : amu weight of unit crystal
    | self.weight_g    : g weight of unit crystal
    """
    def __init__(self, comp):
        self.composition = comp

    @property
    def weight_amu(self):
        result = 0.0
        for atom in self.composition.keys():
            result += self.composition[atom] * amu[atom]
        return result

    @property
    def weight_g(self):
        return self.weight_amu * amu_g

def calc_mass_conversion(elem, mineral):
    """
    | Returns the number of atoms per gram of a particular mineral object
    | Useful for converting mass column to a number density column for an element
    """
    assert type(mineral) == Mineral
    assert type(elem) == str
    return mineral.composition[elem] / mineral.weight_g  # g^{-1}

def calc_element_column(NH, fmineral, atom, mineral, d2g=0.009):
    """
    Calculate the column density of an element for a particular NH value,
    assuming a dust-to-gas ratio (d2g) and the fraction of dust in that
    particular mineral species (fmineral)
    """
    dust_mass = NH * mp * d2g * fmineral  # g cm^{-2}
    print('Dust mass = %.3e g cm^-2' % (dust_mass))
    return calc_mass_conversion(atom, mineral) * dust_mass  # cm^{-2}

def get_ISM_abund(elem, abund_table=wilms):
    """
    Given an abundance table, calculate the number per H atom of a
    given element in any ISM form
    """
    assert type(elem) == str
    assert type(abund_table) == dict
    return np.power(10.0, abund_table[elem] - 12.0)  # number per H atom

def get_dust_abund(elem, abund_table=wilms, gas_ratio=wilms_1mbeta):
    """
    Given an abundance table (dict) and a table of gas ratios (dict),
    calculate the number per H atom of a given ISM element in *solid* form
    """
    assert type(elem) == str
    assert type(abund_table) == dict
    assert type(gas_ratio) == dict
    return get_ISM_abund(elem, abund_table) * (1.0 - gas_ratio[elem])  # number per H atom

max_stars_repo_name: eblur/newdust

max_issues_repo_name: eblur/newdust

max_forks_repo_name: eblur/newdust


text:
import numpy as np

from lib.transition import Transition
from lib.buffer_utils import BufferUtils

class Fuzz(Transition):
    """
    """

    def __init__(self, app):
        Transition.__init__(self, app)

    def __str__(self):
        return "Fuzz"

    def reset(self):
        self.buffer_size = BufferUtils.get_buffer_size()
        self.mask = np.tile(False, (self.buffer_size, 3))

        num_elements = np.ndarray.flatten(self.mask).size / 3
        np.random.seed()
        self.rand_index = np.arange(num_elements)
        np.random.shuffle(self.rand_index)

        self.last_idx = 0

    def get(self, start, end, progress):
        idx = int(progress * len(self.rand_index))
        for i in range(self.last_idx, idx):
            offset = self.rand_index[i] * 3
            self.mask.flat[offset] = True
            self.mask.flat[offset + 1] = True
            self.mask.flat[offset + 2] = True
        self.last_idx = idx

        start[self.mask] = 0.0
        end[np.invert(self.mask)] = 0.0

        return (start) + (end)

max_stars_repo_name: ikea-lisp-code/firemix

max_issues_repo_name: ikea-lisp-code/firemix

max_forks_repo_name: ikea-lisp-code/firemix


text:
import sys
import numpy as np

if __name__== "__main__":
    """
    Arguments: dimensions num-components num-points output-file-name.
    Output file format: let d=dimensions, c=num-components, n=num-points
      d c n
      c lines for each component
        weight m_1 ... m_d (mean) v_00 ... v_d0 v_01 ... v_dd (covariance)
      n lines describing the points
        x_1 ... x_d
    """

    #parse input
    d = int(sys.argv[1])
    c = int(sys.argv[2])
    n = int(sys.argv[3])
    seed = int(sys.argv[4])
    print("Dimensions:",d,", Num Components:",c,", Num Points",n,", Seed", seed)
    outName = sys.argv[5]
    print("Output file name:",outName);

    #ground truth
    np.random.seed(66)
    c_means = np.random.normal(size=[c, d]) * 10
    c_covariances = np.abs(np.random.normal(size=[c, d, d])) * 1 + (1 if n>0 else 5) * np.array([(3+np.abs(np.random.normal()))*np.diag([1]*d) for i in range(c)])
    c_weights = np.abs(np.random.normal(size=[c])) + 1

    #pertubate data
    np.random.seed(seed)
    c_means += np.random.normal(size=[c, d]) * 5
    c_covariances += 1 * np.random.normal(size=[c, 1, 1]) * np.eye(d).reshape((1, d, d))
    c_covariances = np.matmul(c_covariances, np.transpose(c_covariances, (0,2,1)))
    c_weights += (np.random.normal(size=[c]) * 0.01 + 0.01)
    c_weights /= np.sum(c_weights)

    result = np.zeros((n, d), dtype=np.float32)
    for i in range(n):
        comp = np.random.choice(np.array(range(c)), p=c_weights)
        result[i] = np.random.multivariate_normal(
            c_means[comp], c_covariances[comp]
        )

    # write result
    with open(outName, 'w') as f:
        f.write('%d %d %d\n'%(d,c,n))
        for ic in range(c):
            f.write('%f'%c_weights[ic]);
            for x in range(d):
                f.write(' %f'%c_means[ic,x])
            for x in range(d):
                for y in range(d):
                    f.write(' %f'%c_covariances[ic,x,y])
            f.write('\n')
        for i in range(n):
            for x in range(d):
                f.write('%f '%result[i,x])
            f.write('\n')

    print("Done")

max_stars_repo_name: chrismile/cuMat

max_issues_repo_name: chrismile/cuMat

max_forks_repo_name: chrismile/cuMat


text:
import json
import re
import pandas as pd
import numpy as np

import io

from logzero import logger
from flask import Response, Blueprint, render_template, request, url_for, jsonify, redirect, flash, abort

from base.utils.auth import jwt_required, admin_required, get_current_user
from base.forms import PairwiseIndelForm

from caendr.services.dataset_release import get_browser_tracks_path
from caendr.utils.constants import CHROM_NUMERIC
from caendr.utils.data import get_object_hash

from caendr.services.indel_primer import get_sv_strains, query_indels_and_mark_overlaps, create_new_indel_primer, get_indel_primer, fetch_ip_data, fetch_ip_result, get_all_indel_primers, get_user_indel_primers

# Tools blueprint
indel_primer_bp = Blueprint('indel_primer',
                            __name__,
                            template_folder='tools')

from caendr.utils.constants import CHROM_NUMERIC

@indel_primer_bp.route('/pairwise_indel_finder', methods=['GET'])
@jwt_required()
def indel_primer():
  form = PairwiseIndelForm(request.form)
  title = "Pairwise Indel Finder"
  strains = get_sv_strains()
  chroms = CHROM_NUMERIC.keys()
  fluid_container = True
  alt_parent_breadcrumb = {"title": "Tools", "url": url_for('tools.tools')}

  # TODO: REPLACE THESE STATIC VALUES
  divergent_track_summary_url = "//storage.googleapis.com/elegansvariation.org/browser_tracks/lee2020.divergent_regions_all.bed.gz"
  fasta_url = "//storage.googleapis.com/elegansvariation.org/browser_tracks/c_elegans.PRJNA13758.WS276.genomic.fa"

  return render_template('tools/indel_primer/submit.html', **locals())

@indel_primer_bp.route("/pairwise_indel_finder/all-results")
@admin_required()
def indel_primer_all_results():
  title = "All Indel Primer Results"
  alt_parent_breadcrumb = {"title": "Tools", "url": url_for('tools.tools')}
  user = get_current_user()
  items = get_all_indel_primers()
  return render_template('tools/indel_primer/list-all.html', **locals())

@indel_primer_bp.route("/pairwise_indel_finder/my-results")
@jwt_required()
def indel_primer_user_results():
  title = "My Indel Primer Results"
  alt_parent_breadcrumb = {"title": "Tools", "url": url_for('tools.tools')}
  user = get_current_user()
  items = get_user_indel_primers(user.name)
  return render_template('tools/indel_primer/list-user.html', **locals())

@indel_primer_bp.route("/pairwise_indel_finder/query_indels", methods=["POST"])
@jwt_required()
def pairwise_indel_finder_query():
  form = PairwiseIndelForm()
  if form.validate_on_submit():
    strain_1 = form.data['strain_1']
    strain_2 = form.data['strain_2']
    chrom = form.data['chromosome']
    start = form.data['start']
    stop = form.data['stop']
    results = query_indels_and_mark_overlaps(strain_1, strain_2, chrom, start, stop)
    return jsonify(results=results)
  return jsonify({"errors": form.errors})

@indel_primer_bp.route('/pairwise_indel_finder/submit', methods=["POST"])
@jwt_required()
def submit_indel_primer():
  user = get_current_user()
  data = request.get_json()
  data_hash = get_object_hash(data, length=32)

  props = {'site': data['site'],
            'strain_1': data['strain_1'],
            'strain_2': data['strain_2'],
            'size': data['size'],
            'data_hash': data_hash,
            'username': user.name}

  p = create_new_indel_primer(**props)
  return jsonify({'started': True,
                  'data_hash': data_hash,
                  'id': p.id })

# TODO: Move internals of this to a service function
@indel_primer_bp.route("/indel_primer/result/<id>")
@indel_primer_bp.route("/indel_primer/result/<id>/tsv/<filename>")
@jwt_required()
def pairwise_indel_query_results(id, filename = None):
    alt_parent_breadcrumb = {"title": "Tools", "url": url_for('tools.tools')}
    user = get_current_user()
    ip = get_indel_primer(id)

    if (not ip._exists) or (ip.username != user.name):
      flash('You do not have access to that report', 'danger')
      abort(401)

    data_hash = ip.data_hash

    title = "Indel Primer Results"
    data = fetch_ip_data(ip)
    result = fetch_ip_result(ip)
    ready = False

    if data is None:
        return abort(404, description="Indel primer report not found")

    data = json.loads(data.download_as_string().decode('utf-8'))
    logger.debug(data)
    # Get trait and set title
    title = f"Indel Primer Results {data['site']}"
    subtitle = f"{data['strain_1']} | {data['strain_2']}"

    # Set indel information
    size = data['size']
    chrom, indel_start, indel_stop = re.split(":|-", data['site'])
    indel_start, indel_stop = int(indel_start), int(indel_stop)

    if result:
        result = result.download_as_string().decode('utf-8')
        result = pd.read_csv(io.StringIO(result), sep="\t")

        # Check for no results
        empty = True if len(result) == 0 else False
        ready = True
        ip.status = 'COMPLETE'
        ip.empty = empty
        ip.save()
        if empty is False:
            # left primer
            result['left_primer_start'] = result.amplicon_region.apply(lambda x: x.split(":")[1].split("-")[0]).astype(int)
            result['left_primer_stop'] = result.apply(lambda x: len(x['primer_left']) + x['left_primer_start'], axis=1)

            # right primer
            result['right_primer_stop'] = result.amplicon_region.apply(lambda x: x.split(":")[1].split("-")[1]).astype(int)
            result['right_primer_start'] = result.apply(lambda x:  x['right_primer_stop'] - len(x['primer_right']), axis=1)

            # Output left and right melting temperatures.
            result[["left_melting_temp", "right_melting_temp"]] = result["melting_temperature"].str.split(",", expand = True)

            # REF Strain and ALT Strain
            ref_strain = result['0/0'].unique()[0]
            alt_strain = result['1/1'].unique()[0]

            # Extract chromosome and amplicon start/stop
            result[[None, "amp_start", "amp_stop"]] = result.amplicon_region.str.split(pat=":|-", expand=True)

            # Convert types
            result.amp_start = result.amp_start.astype(int)
            result.amp_stop = result.amp_stop.astype(int)

            result["N"] = np.arange(len(result)) + 1
            # Setup output table
            format_table = result[["N",
                                "CHROM",
                                "primer_left",
                                "left_primer_start",
                                "left_primer_stop",
                                "left_melting_temp",
                                "primer_right",
                                "right_primer_start",
                                "right_primer_stop",
                                "right_melting_temp",
                                "REF_product_size",
                                "ALT_product_size"]]

            # Format table column names
            COLUMN_NAMES = ["Primer Set",
                            "Chrom",
                            "Left Primer (LP)",
                            "LP Start",
                            "LP Stop",
                            "LP Melting Temp",
                            "Right Primer (RP)",
                            "RP Start",
                            "RP Stop",
                            "RP Melting Temp",
                            f"{ref_strain} (REF) amplicon size",
                            f"{alt_strain} (ALT) amplicon size"]

            format_table.columns = COLUMN_NAMES

            records = result.to_dict('records')

    if request.path.endswith("tsv"):
        # Return TSV of results
        return Response(format_table.to_csv(sep="\t"), mimetype="text/tab-separated-values")

    return render_template("tools/indel_primer/view.html", **locals())

max_stars_repo_name: AndersenLab/CAENDR

max_issues_repo_name: AndersenLab/CAENDR

max_forks_repo_name: AndersenLab/CAENDR


text:
from config import CIFAR10_config as config
from keras.models import load_model
import numpy as np
import pandas as pd

# Load the test set and model
print('Loading test data and trained model...')
X_test = np.load(config.DATASET + 'X_test_array.npy')
model = load_model(config.MODEL + 'best_model.h5')

# Predict classes
y_test = model.predict(X_test)

# Convert y_test to class names
classes = [
    'airplane', 'automobile', 'bird', 'cat', 'deer',
    'dog', 'frog', 'horse', 'ship', 'truck'
]
pred = []
for i in y_test.argmax(axis=1):
    pred.append(classes[i])
pred = np.array(pred)

# Save to csv in required format
df = pd.DataFrame({
    'id': np.arange(1,X_test.shape[0]+1),
    'label': pred
})
df.to_csv(config.OUTPUT + 'kaggle_predictions.csv', index=False)
print('Predictions saved!')


text:
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm # Displays a progress bar

import torch
from torch import nn
from torch import optim
import torch.nn.functional as F
from torchvision import datasets, transforms , utils
from torch.utils.data import Dataset, Subset, DataLoader, random_split
from torch.utils.data.sampler import SubsetRandomSampler

from dataset import EnableDataset

from tempfile import TemporaryFile
import pickle
from skimage import io, transform
import pdb

from networks import Network

########## SETTINGS  ########################

## for N way classifieres
numb_class= [5,2,2,2,2]
len_class = len(numb_class)

# numb_class = 6
SAVE_MODEL = False
num_epoch = 60

BATCH_SIZE = 32
LEARNING_RATE = 1e-4*0.8
WEIGHT_DECAY = 1e-4

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

def save_object(obj, filename):
    with open(filename, 'wb') as output:  # Overwrites any existing file.
        pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL)

# Load the dataset and train, val, test splits
print("Loading datasets...")

BIO_trains=[]
BIO_tests=[]
events = ['RC','RT','LT','LC']
# for i in range(len_class):
#     for event in events:
#         BIO_trains.append(EnableDataset(subject_list= ['156','185','186','188','189','190', '191', '192', '193', '194'],data_range=(1,46),window_size=500,label=i+1,event_label=event))
#         BIO_tests.append(EnableDataset(subject_list= ['156','185','186','188','189','190', '191', '192', '193', '194'],data_range=(46,50),window_size=500,label=i+1,event_label=event))
#         print("label_{}_event_{} loaded".format(i,event))

# for i in range(len_class):
#     for event in events:
#         BIO_trains.append(EnableDataset(subject_list= ['156'],data_range=(1,4),window_size=500,label=i+1,event_label=event))
#         BIO_tests.append(EnableDataset(subject_list= ['156'],data_range=(46,50),window_size=500,label=i+1,event_label=event))
#         print("label_{}_event_{} loaded".format(i,event))

# save_object(BIO_trains,'BIO_trains_reduced_size.pkl')
# save_object(BIO_tests,'BIO_tests_reduced_size.pkl')

with open('BIO_trains_reduced_size.pkl', 'rb') as input:
    BIO_trains = pickle.load(input)
with open('BIO_tests_reduced_size.pkl', 'rb') as input:
    BIO_tests = pickle.load(input)

## check the class distribution

def weight_classes(dataset):
    trainloader = DataLoader(dataset, shuffle=False,batch_size=BATCH_SIZE)
    classes = [0,0,0,0,0,0,0]
    for data, labels in trainloader:
        for x in range(labels.size()[0]):
            classes[labels[x]] +=1
    # print(classes)

    classes= classes[1:-1]

    ## with sample
    weights=[]
    sum_classes = np.sum(classes)
    for idx in classes:
        if idx != 0 :
            weights.append(sum_classes/idx)
        else:
            continue

    print(weights)
    weights = torch.FloatTensor(weights)

    return weights

weights_list=[]
j = 0
for i in range(len_class):
    for event in events:
        weights_list.append(weight_classes(BIO_trains[j]))
        j+=1

validation_split = .2
shuffle_dataset = True
random_seed= 42

# Creating data indices for training and validation splits:
train_samplers=[]
valid_samplers=[]
val_lens=[]
j = 0
for i in range(len_class):
    for event in events:
        dataset_size = len(BIO_trains[j])
        indices = list(range(dataset_size))
        split = int(np.floor(validation_split * dataset_size))
        if shuffle_dataset :
            np.random.seed(random_seed)
            np.random.shuffle(indices)
        train_indices, val_indices = indices[split:], indices[:split]
        val_len=len(val_indices)
        val_lens.append(val_len)

        # Creating PT data samplers and loaders:
        train_sampler = SubsetRandomSampler(train_indices)
        valid_sampler = SubsetRandomSampler(val_indices)
        train_samplers.append(train_sampler)
        valid_samplers.append(valid_sampler)
        j +=1

trainloaders=[]
valloaders=[]
testloaders=[]
j = 0
for i in range(len_class):
    for event in events:
        trainloaders.append(DataLoader(BIO_trains[j], shuffle=False,batch_size=BATCH_SIZE,sampler = train_samplers[j]))
        valloaders.append(DataLoader(BIO_trains[j], shuffle=False,batch_size=BATCH_SIZE,sampler = valid_samplers[j]))
        testloaders.append(DataLoader(BIO_tests[j], shuffle=False,batch_size=BATCH_SIZE))
        j +=1

device = "cuda" if torch.cuda.is_available() else "cpu" # Configure device
print('GPU USED?',torch.cuda.is_available())

################# MODELS#####################

models =[]
optimizers =[]
criterions =[]
j = 0
for i in range(len_class):
    for event in events:
            model = Network(output_dim=numb_class[i])
            model = model.to(device)
            model.eval()
            models.append(model)
            weights_list[j] = weights_list[j].to(device)
            criterions.append(nn.CrossEntropyLoss(weight=weights_list[j]))
            optimizers.append(optim.Adam(model.parameters(), lr=LEARNING_RATE, weight_decay=WEIGHT_DECAY))
            del model
            j +=1

#### MODE-SPECIFIC CLASSIFIER SCHEME
# 1 : LW, 2:RA, 3:RD, 4:SA, 5:SD, 6:Stand
# MODEL1: LW(1)-> LW(1), RA(2), RD(3), SA(4), SD(5), Stand(6)
# MODEL2: RA(2)-> LW(1), RA(2)
# MODEL3: RD(3)-> LW(1), RD(3)
# MODEL4: SA(4)-> LW(1), SA(4)
# MODEL5: SD(5)-> LW(1), SD(5)
# MODEL6: Stand(6)-> LW(1), SD(6) ## PAPER DOES NOT DO THIS

def train(model, criterion, optimizer, loader, valloader, num_epoch = 20,label_no=None, event=None,vallen=None): # Train the model
    loss_history=[]
    val_history=[]
    print("Start training...")
    model.train()

    for i in range(num_epoch):
        model.train()
        running_loss = []
        for batch, label in tqdm(loader):
            batch = batch.to(device)
            label = label -1 # indexing start from 1 (removing sitting conditon)
            if label_no>2:
                label = label/torch.LongTensor([label_no-1])
            label = label.to(device)
            optimizer.zero_grad()
            pred = model(batch)
            loss = criterion(pred, label)
            running_loss.append(loss.item())
            loss.backward()
            optimizer.step()
            loss_history.append(np.mean(running_loss))
        _,_,val_acc = evaluate(model, valloader,label_no,vallen)
        val_history.append(val_acc)
        print("Epoch {} loss:{} val_acc:{}".format(i+1,np.mean(running_loss),val_acc)) # Print the average loss for this epoch

    if SAVE_MODEL:
        torch.save(model, "model_class{}_{}.pth".format(label_no,event))
        print('model saved label:{} event:{}'.format(label_no,event))
    print("Done!")
    return loss_history, val_history

def evaluate(model, loader,label_no,vallen=None): # Evaluate accuracy on validation / test set
    model.eval()
    correct = 0
    with torch.no_grad():
        for batch, label in tqdm(loader):
            batch = batch.to(device)

            pred = model(batch)
            label = label-1
            if label_no>2:
                label = label/torch.LongTensor([label_no-1])
            label = label.to(device)
            correct += (torch.argmax(pred,dim=1)==label).sum().item()
    if len(loader.dataset) != 0 and vallen != None:
        acc = correct/vallen if vallen != 0 else 0
        print("Evaluation accuracy: {}".format(acc))
    elif vallen == None:
        acc = correct/len(loader.dataset) if len(loader.dataset) != 0 else 0
        print("Evaluation accuracy: {}".format(acc))
    else:
        acc =0
        print('empty dataset!!', label_no)
    return correct, len(loader.dataset), acc

loss_historys=[]
val_historys=[]
# j = 15
j = 0
for i in range(len_class):
    for event in events:

        # print(event)
        loss_history, val_history =train(models[j], criterions[j], optimizers[j], trainloaders[j], valloaders[j], num_epoch, label_no=i+1,event=event,vallen=val_lens[j])
        loss_historys.append(loss_history)
        val_historys.append(val_history)
        print('*****class{}_{}*********'.format(i+1,event))
        j +=1

print("Evaluate on test set")

corrs=[]
len_datas=[]
j =0
for i in range(len_class):
    for event in events:
        corr, len_data,_ =evaluate(models[j], testloaders[j],label_no=i+1 )
        corrs.append(corr)
        len_datas.append(len_data)
        print('Evlauation on class{}_{}*********'.format(i+1,event))
        j +=1

corr_total =0
len_data_total=0
j =0
for i in range(len_class):
    for event in events:
        corr_total = corrs[j] + corr_total
        len_data_total = len_data_total + len_datas[j]
        j +=1

acc_total = corr_total/len_data_total
print("Total Evaluation accuracy: {}".format(acc_total))

np.savetxt('loss_history.txt',loss_history)
np.savetxt('val_history.txt',val_history)

fig =plt.figure()
plt.plot(loss_history,label='train loss')
plt.xlabel('iteration')
plt.ylabel('loss')
plt.legend()
plt.show()
fig.savefig('train_loss.jpg')

fig =plt.figure()
plt.plot(val_history,label='vaidation accuracy')
plt.xlabel('epoch')
plt.ylabel('accuracy')
list_idx = [ i for i in range(num_epoch)]
plt.xticks(np.array(list_idx))
plt.legend()
plt.show()
fig.savefig('val_acc.jpg')