text:
import argparse
import copy
import json
import logging
import os

import numpy as np
import torch
import transformers
from data import get_data
from evaluate import EVAL_FUNC_MAP
from model import MODEL_CLASS_MAP
from tqdm import tqdm

logging.basicConfig(
    format="%(asctime)s - %(levelname)s - %(name)s -   %(message)s",
    datefmt="%m/%d/%Y %H:%M:%S",
    level=logging.INFO,
)

logger = logging.getLogger(__name__)

try:
    import torch_xla.core.xla_model as xm
    import torch_xla.distributed.data_parallel as xla_dp
except ImportError as e:
    logger.error(f"Failed to import XLA. {e}")

def set_seed(seed):
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = True

def cal_running_avg_loss(loss, running_avg_loss, decay=0.99):
    if running_avg_loss == 0:
        return loss
    running_avg_loss = running_avg_loss * decay + (1 - decay) * loss
    return running_avg_loss

def _run_epoch(model, loader, device=None, context=None, **kwargs):
    config = kwargs["config"]
    is_train = kwargs["is_train"]

    avg_loss = 0
    results = []
    batch_num = len(loader)

    if is_train:
        model.train()
        if config.use_tpu:
            optimizer = context.getattr_or(
                "optimizer",
                lambda: model.get_optimizer(),
            )
            scheduler = context.getattr_or(
                "scheduler",
                lambda: model.get_scheduler(batch_num, optimizer),
            )
        else:
            optimizer = kwargs["optimizer"]
            scheduler = kwargs["scheduler"]
    else:
        model.eval()

    is_master = True
    if config.use_tpu:
        is_master = xm.is_master_ordinal()

    pbar = tqdm(enumerate(loader), total=batch_num, disable=not is_master, dynamic_ncols=True)
    for i, inputs in pbar:

        if not config.use_tpu:
            for k, v in inputs.items():
                if isinstance(v, torch.Tensor):
                    inputs[k] = v.to(device)

        outputs = model(inputs)
        loss = outputs.loss.mean()
        avg_loss = cal_running_avg_loss(loss.item(), avg_loss)
        loss /= config.gradient_accumulation_steps

        if is_train:
            loss.backward()
            if i % config.gradient_accumulation_steps == 0 or i == batch_num - 1:

                if config.max_grad_norm > 0:
                    torch.nn.utils.clip_grad_norm_(model.parameters(), config.max_grad_norm)

                if config.use_tpu:
                    xm.optimizer_step(optimizer)
                else:
                    optimizer.step()
                optimizer.zero_grad()

                if scheduler is not None:
                    scheduler.step()
        else:
            result = (model.module if hasattr(model, "module") else model).eval_step(inputs, outputs)
            results.extend(result)

        if is_master:
            pbar.set_description(
                f"epoch: {kwargs['epoch'] + 1}, {('train' if is_train else 'valid')} loss: {min(100, round(avg_loss, 4))}"
            )

    return {
        "loss": avg_loss,
        "result": results,
    }

def run_epoch(**kwargs):
    model = kwargs.pop("model")
    if kwargs["config"].use_tpu:
        results = model(_run_epoch, **kwargs)
    else:
        results = _run_epoch(model, **kwargs)

    if isinstance(results, list):
        loss = sum([result["loss"] for result in results]) / len(results)
        result = []
        for res in results:
            result.extend(res["result"])
        results = {"loss": loss, "result": result}

    return results

def run(parser):
    # NOTE Remove redundant bigbird logs
    transformers.logging.get_logger("transformers.models.big_bird.modeling_big_bird").setLevel(logging.ERROR)

    args, _ = parser.parse_known_args()

    model = MODEL_CLASS_MAP.get(args.task, None)
    if model is None:
        raise Exception(f"Invalid model task {args.task}")

    parser = model.add_args(parser)
    config = parser.parse_args()

    set_seed(config.seed)

    model = model(config)

    logger.info(f"configuration: {str(config)}")

    if config.use_tpu:
        devices = xm.get_xla_supported_devices()
        model_dp = xla_dp.DataParallel(model, device_ids=devices)
    else:
        if torch.cuda.is_available():
            gpu_count = torch.cuda.device_count()
            logger.info(f"{gpu_count} GPU device detected")
            devices = ["cuda:{}".format(i) for i in range(gpu_count)]
            model_dp = torch.nn.DataParallel(model, device_ids=devices)
            model.to(devices[0])
        else:
            devices = ["cpu"]
            model_dp = model

    config.world_size = len(devices)
    if config.do_train:
        train_loader = get_data(config, tokenizer=model.tokenizer)
    valid_loader = get_data(config, tokenizer=model.tokenizer, is_train=False)

    optimizer = None
    scheduler = None
    if not config.use_tpu and config.do_train:
        optimizer = model.get_optimizer()
        scheduler = model.get_scheduler(len(train_loader), optimizer)

    params = {
        "config": config,
        "model": model_dp,
        "optimizer": optimizer,
        "scheduler": scheduler,
    }
    if not config.use_tpu:
        params["device"] = devices[0]

    def do_eval(epoch):
        with torch.no_grad():
            results = run_epoch(loader=valid_loader, epoch=epoch, is_train=False, **params)["result"]
            results = EVAL_FUNC_MAP[config.task](
                config=config,
                model=model,
                loader=valid_loader,
                tokenizer=model.tokenizer,
                results=results,
            )

        logger.info("Eval results.")
        for k, v in results["results"].items():
            logger.info(f"{k} : {v}")

        return results["best_score"]

    if config.do_train:
        best_score = 0
        for epoch in range(config.num_train_epochs):
            run_epoch(loader=train_loader, epoch=epoch, is_train=True, **params)

            score = 0
            if config.do_eval_during_train:
                score = do_eval(epoch)

            if score >= best_score:
                best_score = score
                output_dir = os.path.join(config.output_dir, config.task, config.dataset, f"{epoch}-{best_score}-ckpt")
                copy.deepcopy(
                    model_dp.module
                    if hasattr(model_dp, "module")
                    else model_dp._models[0]
                    if hasattr(model_dp, "_models")
                    else model_dp
                ).cpu().save_pretrained(output_dir)
                with open(os.path.join(output_dir, "finetune_config.json"), "w") as save_config:
                    json.dump(vars(config), save_config, sort_keys=True, indent=4)
                logger.info(f"Checkpoint {output_dir} saved.")

    if config.do_eval:
        do_eval(-1)

def main(parser):
    args, _ = parser.parse_known_args()

    if not os.path.exists(args.cache_dir):
        os.makedirs(args.cache_dir)

    output_dir = os.path.join(args.output_dir, args.task, args.dataset)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    run(parser)

if __name__ == "__main__":
    parser = argparse.ArgumentParser()

    parser.add_argument("--task", default=None, type=str, required=True)
    parser.add_argument("--dataset", default=None, type=str, required=True)
    parser.add_argument("--cache_dir", default="cache", type=str)
    parser.add_argument("--output_dir", default="output", type=str)
    parser.add_argument("--do_train", action="store_true")
    parser.add_argument("--do_eval_during_train", action="store_true")
    parser.add_argument("--do_eval", action="store_true")
    parser.add_argument("--use_tpu", action="store_true")
    parser.add_argument("--threads", default=4, type=int)
    parser.add_argument("--seed", default=42, type=int)

    main(parser)

max_stars_repo_name: sackoh/KoBigBird

max_issues_repo_name: sackoh/KoBigBird

max_forks_repo_name: sackoh/KoBigBird


text:
import numpy
import random
random.seed(1234)   # ONLY FOR TESTING, COMMENT OUT FOR BETTER RANDOMNESS
import PyPluMA

EPS=1e-8

class ExponentialPlugin:
   def input(self, filename):
      filestuff = open(filename, 'r')
      firstline = filestuff.readline()
      self.bacteria = firstline.split(',')
      if (self.bacteria.count('\"\"') != 0):
         self.bacteria.remove('\"\"')
      self.n = len(self.bacteria)

      self.ADJ = numpy.matrix(numpy.zeros([self.n, self.n]))
      self.P = numpy.matrix(numpy.zeros([self.n, 1]))

      i = 0
      for line in filestuff:
         contents = line.split(',')
         for j in range(self.n):
            if (i != j):
              self.ADJ[i,j] = float(contents[j+1])
         i += 1

      for i in range(self.n):
         self.bacteria[i] = self.bacteria[i].strip()
         self.P[i,0] = 1.0/self.n

      self.EPS = random.random()*0.000001
      self.MU = (numpy.amax(self.ADJ) - numpy.amin(self.ADJ)) / 2.0 + self.EPS

   def run(self):
      oldP = numpy.zeros([self.n, 1])
      iteration = 0
      while (numpy.amax(self.P - oldP) > self.EPS):
         oldP = self.P.copy()
         tmp = numpy.exp(1.0/self.MU)*numpy.transpose(self.ADJ)*self.P
         self.P = tmp / numpy.linalg.norm(tmp)
         iteration += 1
      self.k = numpy.transpose(self.ADJ)*self.P
      PyPluMA.log("*** Converged after "+str(iteration)+" iterations ***")

   def output(self, filename):
     UG = []
     for i in range(self.n):
      UG.append((abs(self.k[i,0]), self.bacteria[i]))
     UG.sort()
     UG.reverse()
     # Formatted for Cytoscape, borrowed from ATria
     outfile = open(filename, 'w')
     outfile.write("Name\tCentrality\tRank\n")
     centvals = numpy.zeros([len(UG)])
     for i in range(len(UG)):
       PyPluMA.log(str((UG[i][1], UG[i][0])))
       bac = UG[i][1]
       rank = len(UG)-i
       if (i != 0 and abs(UG[i][0]-UG[i-1][0]) < EPS):
          rank = oldrank
       oldrank = rank
       if (bac[0] == '\"'):
          bac = bac[1:len(bac)-1]
       #if (UG[i][0] != UG[len(UG)-1][0]):
       outfile.write(bac+"\t"+str(abs(UG[i][0]))+"\t"+str(rank)+"\n")
       #else:
       #  outfile.write(bac+"\t"+str(abs(UG[i][0]))+"\t"+"0\n")
       centvals[i] = abs(UG[i][0])

     PyPluMA.log("Wrote file: "+filename)
     PyPluMA.log("Min centrality: "+str(numpy.min(centvals)))
     PyPluMA.log("Max centrality: "+str(numpy.max(centvals)))
     mymean = numpy.mean(centvals)
     stddev = numpy.std(centvals)
     PyPluMA.log("Standard Deviation: "+str(stddev))
     PyPluMA.log("Two STDs back: "+str(mymean - 2*stddev))
     PyPluMA.log("One STD back: "+str(mymean - stddev))
     PyPluMA.log("One STD forward: "+str(mymean + stddev))
     PyPluMA.log("Two STDs forward: "+str(mymean + 2*stddev))

ext: py

reason: import numpy


text:
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8
#
# Capriqorn --- CAlculation of P(R) and I(Q) Of macRomolcules in solutioN
#
# Copyright (c) Juergen Koefinger, Klaus Reuter, and contributors.
# See the file AUTHORS.rst for the full list of contributors.
#
# Released under the GNU Public Licence, v2 or any higher version, see the file LICENSE.txt.

"""Capriqorn deltaH filter.
"""

from __future__ import division
from past.utils import old_div

import numpy as np
import copy

from cadishi import base
from cadishi import util

from ...lib import rdf
from ...lib import selection
from ...lib import formFactor as ff
from ...lib import pyntensities as pynt

def _make_2d(dict_in):
    """
    Convert a dictionary of 1d NumPy arrays into a
    2d NumPy array.  Return a tuple containing the
    2d array and the sorted list of keys.
    """
    # ---
    assert ('radii' in list(dict_in.keys()))
    # ---
    np_2d = np.zeros((len(dict_in['radii']), len(list(dict_in.keys()))))
    np_2d[:, 0] = (dict_in['radii'])[:]
    np_2d_keys = sorted(dict_in.keys())
    np_2d_keys.remove('radii')
    idx = 1
    for key in np_2d_keys:
        np_2d[:, idx] = (dict_in[key])[:]
        idx += 1
    # ---
    return (copy.deepcopy(np_2d), copy.deepcopy(np_2d_keys))

def _separate_ghosts_from_histogram(H_in):
    """
    Create two disjunct dictionaries from H_in, one (H) containing the real
    species histograms, another (Hx) containing the X particle histograms.
    """
    H = {}
    Hx = {}
    for key in H_in:
        if (key == 'radii'):
            H[key] = H_in[key].copy()
            Hx[key] = H_in[key].copy()
            continue
        if ('X' in key):
            Hx[key] = H_in[key].copy()
        else:
            H[key] = H_in[key].copy()
    return (H, Hx)

# ---
# Implementation mostly follows <lsz/deltaH.py>.
# ---
class DeltaH(base.Filter):
    """A filter that computes the intensity."""
    _depends = ["Solvent"]
    _conflicts = []

    def __init__(self,
                 form_factor_file='atomsf.dat',
                 nq=300,
                 dq=0.01,
                 source=-1,
                 debug=False,
                 verbose=False):
        self.form_factor_file = form_factor_file
        self.nq = nq
        self.dq = dq
        # --- Note: The following two parameters are obtained from the pipeline log!
        self.geometry = None
        self.x_particle_method = None
        # ---
        self.src = source
        self.debug = debug
        self.verb = verbose
        # ---
        self._depends.extend(super(base.Filter, self)._depends)
        self._conflicts.extend(super(base.Filter, self)._conflicts)

    def get_meta(self):
        """
        Return information on the present filter,
        ready to be added to a frame object's list of
        pipeline meta information.
        """
        meta = {}
        label = 'DeltaH'
        param = {'form_factor_file': self.form_factor_file,
                 'nq': self.nq,
                 'dq': self.dq,
                 # --- we keep the following two entries for log purposes
                 'geometry': self.geometry,
                 'x_particle_method': self.x_particle_method}
        meta[label] = param
        return meta

    def __iter__(self):
        return self

    def __next__(self):
        for hs in next(self.src):
            # hs was chosen to mean 'histogram set' ...
            if hs is not None:
                assert isinstance(hs, base.Container)

                # --- obtain information from the pipeline log
                dr = hs.query_meta('histograms/histogram/dr')
                assert (dr is not None)
                # ---
                self.geometry = hs.get_geometry()
                assert (self.geometry is not None)
                geometry_param = hs.query_meta(self.geometry)
                assert (geometry_param is not None)
                if (self.geometry == 'Sphere' or self.geometry == 'Cuboid' or self.geometry == 'Ellipsoid'):
                    V = geometry_param['volume']
                    assert (V is not None)
                    # VAvg is neede below for MultiReference, where V=1 for histograms
                    # but VAvg is used for densities
                    VAvg = V
                else:
                    # --- for non-sphere geometries, we calculate V using x particles below
                    V = None
                # ---
                virtual_param = hs.query_meta('VirtualParticles')
                if (virtual_param is not None):
                    self.x_particle_method = virtual_param['method']
                    xrho = virtual_param['x_density']
                # ---
                rdf_header = list((hs.get_data(base.loc_solv_match + '/g_scaled')).keys())
                rdf_header.remove('radii')
                rdf_elements = util.get_elements(rdf_header)
                n_solv = len(rdf_elements)
                hs_full = selection.get_full(hs)
                n_part_elements = util.get_elements(list(hs_full.get_data(base.loc_nr_particles).keys()))
                histo_elements = util.get_elements(list(hs_full.get_data(base.loc_histograms).keys()))
                # --- print '###', n_part_elements, histo_elements
                assert (n_part_elements == histo_elements)
                n_prot = len(histo_elements)
                if (self.debug):
                    print(" rdf_header", rdf_header)
                    print(" rdf_elements", rdf_elements)
                    print(" n_solv =", n_solv)
                    print(" n_prot=", n_prot)

                histo = copy.deepcopy(hs_full.get_data(base.loc_histograms))
                nr = len(hs_full.get_data(base.loc_histograms + '/radii'))
                for key in list(histo.keys()):
                    if (key == 'radii'):
                        continue
                    histo[key] *= 2.0

                # --- NOTE : in the following, we're sorting out the virtual particles
                H, Hx = _separate_ghosts_from_histogram(histo)

                # --- consistency checks of the radii
                assert (len(hs.get_data(base.loc_solv_match + '/g_scaled/radii')) >= nr)
                radii_1 = hs.get_data(base.loc_solv_match + '/g_scaled/radii')[0:nr]
                radii_2 = hs_full.get_data(base.loc_histograms + '/radii')[0:nr]
                eps = 1.e-9
                assert np.all(np.fabs(radii_1 - radii_2) < eps)
                #
                rdf_expanded = {}
                for key in list((hs.get_data(base.loc_solv_match + '/g_scaled')).keys()):
                    # Note: g_scaled dictionary includes the radii array
                    rdf_expanded[key] = np.zeros(nr)
                    (rdf_expanded[key])[0:nr] = (hs.get_data(base.loc_solv_match + '/g_scaled/' + key))[0:nr]
                rdf_org = copy.deepcopy(rdf_expanded)
                # ---
                if (self.debug):
                    hs.put_data(base.loc_delta_h + '/debug/rdf_org', rdf_org)
                rho = hs.get_data(base.loc_solv_match + '/rho')
                # ---
                for key in list(rdf_expanded.keys()):
                    # zero densities as done in the reference code <deltaH.py>
                    if (key == 'radii'):
                        continue
                    (rdf_expanded[key])[0] = 0.0

                if (self.geometry == 'Sphere') and (self.x_particle_method is None):
                    # print "### delta_H :: SPHERE BRANCH ###"
                    R = geometry_param['radius']
                    dHs = copy.deepcopy(H)
                    for key in list(dHs.keys()):
                        if (key == 'radii'):
                            continue
                        (dHs[key])[:] = 0.0
                    dHp = copy.deepcopy(dHs)
                    dH = copy.deepcopy(dHs)
                    # ---
                    assert hs.contains_key(base.loc_len_histograms)
                    # --- construct h matrix
                    x_list = hs.get_data(base.loc_len_histograms + '/radii')
                    r_list = H['radii']
                    # --- implement h using a dictionary
                    h_dict = {}
                    for key in hs.get_keys(base.loc_len_histograms):
                        h_dict[key] = np.zeros(len(r_list))
                    h_dict['radii'][:] = r_list[:]
                    # ---
                    pSphere = np.zeros((r_list.shape[0], 2))
                    pSphere[:, 0] = r_list[:]
                    # ---
                    for j, rv in enumerate(r_list):
                        pSphere[j, 1] = rdf.P(rv, R)
                    if (self.debug):
                        hs.put_data(base.loc_delta_h + '/debug/pSphere', pSphere)
                    # ---
                    dx = dr  # DANGER! However, these two should be the same anyway!
                    # ---
                    getSr = rdf.SrFast
                    print()
                    for j, rv in enumerate(r_list):
                        if self.verb and (j % 1000 == 0):
                            print(" DeltaH:bin ", j)
                        Sr = getSr(x_list, R, rv)
                        for key in hs.get_keys(base.loc_len_histograms, skip_keys=['radii']):
                            (h_dict[key])[j] += (hs.get_data(base.loc_len_histograms + '/' + key)[:] * Sr[:]).sum(axis=0)

                    for key in hs.get_keys(base.loc_len_histograms, skip_keys=['radii']):
                        h_dict[key] *= 2.0 * dx / V
                        # --- patch rho[] such that it can be used in the block below
                        if key not in rho:
                            rho[key] = 0.0

                    for key in sorted(H.keys()):
                        if (key == 'radii'):
                            continue
                        # --- patch rdf_expanded such that it can be used below
                        if key not in rdf_expanded:
                            rdf_expanded[key] = np.zeros(nr)
                            # ---
                        (el1, el2) = key.split(',')
                        if (el1 == el2):
                            fac = 1.0
                        else:
                            fac = 2.0
                        # ---
                        (dHs[key])[:] = fac * ((rdf_expanded[key])[:] - 2.0) \
                            * rho[el1] * rho[el2] * V * V * pSphere[:, 1] * dr  # Hx[:,ixx]
                        # ---
                        if (fac == 1):
                            (dHp[key])[:] = (H[key])[:] - V * rho[el1] * (h_dict[el1])[:]  # *Hx[:,ix1]
                        else:
                            (dHp[key])[:] = (H[key])[:] - V * (rho[el2] * (h_dict[el1])[:] +
                                                               rho[el1] * (h_dict[el2])[:])
                        # ---
                        (dH[key])[:] = (dHp[key])[:] - (dHs[key])[:]
                else:
                    # --- virtual particle method branch
                    # print "delta_h "
                    nx = np.mean(hs_full.get_data(base.loc_nr_particles + '/X'), dtype=np.float64)
                    if not (self.geometry == 'Sphere' or self.geometry == 'Cuboid' or self.geometry == 'Ellipsoid'):
                        V = old_div(nx, xrho)
                    VAvg = V
                    if (self.geometry == 'MultiReferenceStructure'):
                        V = 1.
                        # print " V =", V, "VAvg =", VAvg
                    else:
                        for key in Hx:
                            if key in 'radii':
                                continue
                            pair = key.split(',')
                            assert (len(pair) == 2)
                            if (pair[0] == 'X') and (pair[1] == 'X'):
                                Hx[key] /= (Hx[key]).sum()
                            else:
                                Hx[key] /= nx
                    # ---
                    dHs = copy.deepcopy(H)
                    for key in list(dHs.keys()):
                        if (key == 'radii'):
                            continue
                        (dHs[key])[:] = 0.0
                    dHp = copy.deepcopy(dHs)
                    dH = copy.deepcopy(dHs)
                    # ---
                    for key in sorted(H.keys()):
                        if (key == 'radii'):
                            continue
                        # --- patch rdf_expanded such that it can be used below
                        if key not in rdf_expanded:
                            rdf_expanded[key] = np.zeros(nr)
                        # ---
                        (el1, el2) = key.split(',')
                        # --- patch rho such that it can be used below
                        for el in [el1, el2]:
                            if el not in rho:
                                rho[el] = 0.0
                        # ---
                        if (el1 == el2):
                            fac = 1.0
                        else:
                            fac = 2.0
                        # ---
                        (dHs[key])[:] = fac * ((rdf_expanded[key])[:] - 2.0) \
                            * rho[el1] * rho[el2] * V * V * Hx['X,X']
                        # ---
                        if (fac == 1):
                            (dHp[key])[:] = (H[key])[:] - V * rho[el1] * (Hx[el1 + ',X'])[:]
                        else:
                            (dHp[key])[:] = (H[key])[:] - V * (rho[el2] * (Hx[el1 + ',X'])[:] +
                                                               rho[el1] * (Hx[el2 + ',X'])[:])
                        # ---
                        (dH[key])[:] = (dHp[key])[:] - (dHs[key])[:]
                # ---
                if (self.debug):
                    hs.put_data(base.loc_delta_h + '/debug/rdf_expanded', rdf_expanded)
                    hs.put_data(base.loc_delta_h + '/debug/dHs', dHs)
                    hs.put_data(base.loc_delta_h + '/debug/dHp', dHp)
                    hs.put_data(base.loc_delta_h + '/debug/dH', dH)
                    hs.put_data(base.loc_delta_h + '/debug/H', H)
                    hs.put_data(base.loc_delta_h + '/debug/rho', rho)
                    hs.put_data(base.loc_delta_h + '/debug/h_dict', h_dict)

                # ---
                hs.put_data(base.loc_intensity + '/dH', dH)

                # --- calculate average particle numbers
                nr_part_avg = {}
                for key in hs_full.get_keys(base.loc_nr_particles):
                    nr_part = hs_full.get_data(base.loc_nr_particles + '/' + key)
                    nr_part_avg[key] = nr_part.sum(axis=0) / float(nr_part.shape[0])
                # --- calculate difference
                nr_part_diff = {}
                for key in nr_part_avg:
                    if key in rho:
                        nr_part_diff[key] = nr_part_avg[key] - rho[key] * VAvg
                    else:
                        nr_part_diff[key] = nr_part_avg[key]

                hs.put_data(base.loc_intensity + '/dN', nr_part_diff)

                # --- INTENSITY CALCULATION ---

                # --- read atom form factors from file
                ff_dict = ff.readAtomSFParam(self.form_factor_file)

                # --- prepare 2D input array for library routine
                (dH_2d, dH_2d_keys) = _make_2d(dH)
                if (self.debug):
                    hs.put_data(base.loc_delta_h + '/debug/dH_2d', dH_2d)
                # ---
                _pIntInter, dInter = pynt.intensitiesFFFaster(self.nq, self.dq,
                                                              dH_2d, dH_2d_keys, ff_dict, 1)
                # # --- prepare input
                # nr_part_keys = sorted(nr_part_diff.keys())
                # nr_part_vals = []
                # for key in nr_part_keys:
                #     nr_part_vals.append(nr_part_diff[key])
                # # ---
                # _pIntIntra, dIntra = pynt.intensitiesFFIntraAtom(self.nq, self.dq, \
                #                             nr_part_vals, nr_part_keys, ff_dict)
                # # ---
                # hs.put_data(base.loc_intensity+'/dI_inter', dInter)
                # hs.put_data(base.loc_intensity+'/dI_intra', dIntra)
                dI = dInter.copy()
                # TODO: check, the following line was commented out in <deltaH.py>
                # dI[:,1] += dIntra[:,1]
                hs.put_data(base.loc_intensity + '/dI', dI)

                # --- prepare 2D input array for library routine
                #     (cannot use _make_2d() routine here)
                n_row = len(rdf_org['radii'])
                n_col = len(list(dH.keys()))
                rdf_org_2d = np.zeros((n_row, n_col))
                rdf_org_2d[:, 0] = (rdf_org['radii'])[:]
                rdf_org_2d_keys = sorted(rdf_org.keys())
                rdf_org_2d_keys.remove('radii')
                for key in rdf_org_2d_keys:
                    if key in dH_2d_keys:
                        idx = dH_2d_keys.index(key) + 1
                        rdf_org_2d[:, idx] = (rdf_org[key])[:]
                # ---
                rho_list = []
                rho_keys = sorted(rho.keys())
                for key in ['X', 'X1', 'X2']:
                    if key in rho_keys:
                        rho_keys.remove(key)
                for key in rho_keys:
                    rho_list.append(rho[key])
                # ---
                if (self.debug):
                    hs.put_data(base.loc_delta_h + '/debug/rdf_org_2d', rdf_org_2d)
                # ---
                bulkH = pynt.getBulkIntegrand(rdf_org_2d, rho_list)
                bulkH_keys = dH_2d_keys
                # print "bulkH_keys", bulkH_keys
                _pIntInter, dIntensityInter = pynt.intensitiesFFFaster(self.nq,
                                                                       self.dq, bulkH, bulkH_keys, ff_dict, dr)
                _pIntIntra, dIntensityIntra = pynt.intensitiesFFIntraAtom(self.nq,
                                                                          self.dq, rho_list, rho_keys, ff_dict)
                dIntensity = dIntensityInter.copy()
                dIntensity[:, 1] += dIntensityIntra[:, 1]
                # ---
                hs.put_data(base.loc_intensity + '/I_solv_inter', dIntensityInter)
                hs.put_data(base.loc_intensity + '/I_solv_intra', dIntensityIntra)
                hs.put_data(base.loc_intensity + '/I_solv', dIntensity)

                # TODO: Think about stripping unnecessary information
                # (histograms, length histograms) from hs at this point.

                hs.put_meta(self.get_meta())
                if self.verb:
                    print("DeltaH.next() :", hs.i)
                yield hs
            else:
                yield None

max_stars_repo_name: bio-phys/capriqorn

max_issues_repo_name: bio-phys/capriqorn

max_forks_repo_name: bio-phys/capriqorn


text:
import numpy as np
import pandas as pd
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from PIL import Image
import PIL.ImageOps
x,y=fetch_openml("mnist_784",version=1,return_X_y=True)
X_train,X_test,Y_train,Y_test=train_test_split(x,y,random_state=9,train_size=7500,test_size=2500)
X_trainscaled=X_train/255.0
X_testS=X_test/255.0
clf=LogisticRegression(solver="saga",multi_class="multinomial").fit(X_trainscaled,Y_train)
def get_predict(image):
    im_pil = Image.open(image)
    image_bw = im_pil.convert('L')
    image_bw_resized = image_bw.resize((28,28), Image.ANTIALIAS)
    pixel_filter = 20
    min_pixel = np.percentile(image_bw_resized, pixel_filter)
    image_bw_resized_inverted_scaled = np.clip(image_bw_resized-min_pixel, 0, 255)
    max_pixel = np.max(image_bw_resized)
    image_bw_resized_inverted_scaled = np.asarray(image_bw_resized_inverted_scaled)/max_pixel
    test_sample = np.array(image_bw_resized_inverted_scaled).reshape(1,784)
    test_pred = clf.predict(test_sample)
    return test_pred[0]

max_stars_repo_name: venkasam/jin

max_issues_repo_name: venkasam/jin

max_forks_repo_name: venkasam/jin


text:
import random
import networkx as nx

def mutation(graph, del_orig=True):

    edges = list(graph.edges)
    nonedges = list(nx.non_edges(graph))

    chosen_edge = random.choice(edges)
    chosen_nonedge = random.choice([x for x in nonedges if chosen_edge[0] == x[0]])

    if del_orig:
        # delete chosen edge
        graph.remove_edge(*chosen_edge)
        # add new edge
    graph.add_edge(*chosen_nonedge)
    return graph

def mutation_simple(graph):
    nonedge = random.choice(list(nx.non_edges(graph)))
    graph.add_edge(*nonedge)
    return nx.minimum_spanning_tree(graph)

max_stars_repo_name: sz-piotr/alhe-rail-network

max_issues_repo_name: sz-piotr/alhe-rail-network

max_forks_repo_name: sz-piotr/alhe-rail-network

reason: import networkx


text:
""" Tests for analyse effect of sensing dot on fitting of tunnel barrier."""

import unittest
import numpy as np
from qtt.algorithms.chargesensor import DataLinearizer

class TestOneDot(unittest.TestCase):

    @staticmethod
    def test_data_linearizer():
        x = np.arange(0, 10, .1)
        y = x + .05 * x**2
        dl = DataLinearizer(x, y)
        _ = dl.forward([1])
        _ = dl.forward_curve([1])

max_stars_repo_name: codecrap/qtt

max_issues_repo_name: codecrap/qtt

max_forks_repo_name: codecrap/qtt


text:
#coding=utf-8
from model import *
import redis
#{shop_id:[user_pay]}
import datetime
import pandas as pd
from pandas import read_csv
from pandas import DataFrame, Series
import statsmodels.api as sm
from statsmodels.tsa.arima_model import ARIMA
from sklearn.metrics import mean_squared_error

from matplotlib import pyplot as plt
import statsmodels.api as sm
import numpy as np
import math

def connect_redis():
    client =redis.Redis(host='127.0.0.1',port=6379,db=0)
    return client

def distinct_weekday_weekend(str_left,str_right,date_format):
    weekday = []
    weekend = []
    start = datetime.datetime.strptime(str_left, date_format)
    while  start <= datetime.datetime.strptime(str_right, date_format):
        if start.weekday() in range(0, 5):
            weekday.append(start.strftime(date_format))
        else:
            weekend.append(start.strftime(date_format))
        start += datetime.timedelta(days = 1)
    return weekday, weekend

def mysql_to_redis(client):
    #count user by day
    shops_sales = {}
    for shop_id in range(1,2001):
        print shop_id
        #group by sql too show
        #shop_query = UserPay.select(fn.DATE(UserPay.pay_time).alias('day'),fn.Count(UserPay.user_id).alias('user_per_shop')).where(UserPay.shop_id== str(shop_id)).group_by(SQL('day')).order_by(SQL('day asc'))
        shop_query = UserPay.select(UserPay.pay_time).where(UserPay.shop_id == str(shop_id))
        #shops_sales[shop_id] = []
        day_count = {}
        for shop in shop_query:
            day_count[shop.pay_time.strftime('%Y-%m-%d')] = day_count.setdefault(shop.pay_time.strftime('%Y-%m-%d'), 0) + 1

        shops_sales[shop_id] = day_count
    #write to redis
    for shop_id in shops_sales.keys():
        for day in shops_sales[shop_id].keys():
            client.hset(shop_id, day, shops_sales[shop_id][day])
    #put day count into array

def sorted_dict_values(adict):
    keys = adict.keys()
    keys.sort()
    return [(k, adict[k]) for k in keys]

def get_pq_with_aic(data):
    import warnings
    with warnings.catch_warnings():
    # catch a hessian inversion and convergence failure warning
        warnings.simplefilter("ignore")
        return sm.tsa.arma_order_select_ic(np.array(data),max_ar=6,max_ma=4,ic='aic')['aic_min_order']

def get_pq_with_test(data):
    pmax=6
    qmax=4
    bic_matrix=[] #bic矩阵
    for p in range(pmax+1):
    	tmp=[]
    for q in range(qmax+1):
        try: #存在部分报错，所以用try来跳过报错。
            tmp.append(ARIMA(data,(p,1,q)).fit().aic)
        except:
            tmp.append(None)
    bic_matrix.append(tmp)
    bic_matrix=pd.DataFrame(bic_matrix) #从中可找出最小值
    p,q=bic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小位置。
    return p, q

def arma_train(data):
    from statsmodels.tsa.stattools import adfuller as ADF
    p, q = get_pq_with_aic(data)
    model = sm.tsa.ARMA(data,(p, q)).fit()
    return model

def sarima_train(data):
    model = sm.tsa.statespace.SARIMAX(data, trend='n', order=(0,1,0), seasonal_order=(1,1,1,7)).fit()
    return model

def sarima_predict(data, start, end,  model):
    predicts = model.predict(start, end, dynamic = True)
    #print predicts, len(predicts)
    #print data, len(data)
    #plt.plot(range(14), predicts)
    #plt.show()
    #plot_data(data, predicts)
    return map(int, predicts)

def arima_train(data):
    '''
        train with algorithm arima
    '''
    #train unique shop with history data
    #ADF test
    print '====================='
    print len(data), data
    from statsmodels.tsa.stattools import adfuller as ADF
    print ADF(pd.Series(np.array(data)))[1]
    p, q = get_pq_with_test(data)
    #p, q = get_pq_with_aic(data)
    print 'ARIMA p值和q值为: %s、%s'%(p,q) #0,1
    #p, q = 0, 0

    model=ARIMA(data,(p,1,q)).fit() #建立ARIMA模型
    return model

def arma_predict(data, start, model):
    predicts = model.predict(start, start + len(data) - 2, dynamic = False)
    print predicts, len(predicts)
    print data, len(data)
    plot_data(data, predicts)

def arima_train_broadcast(data):
    X = data
    size = int(len(X) * 0.66)
    train, test = X[0:size], X[size:len(X)]
    history = [x for x in train]
    predictions = list()
    for t in range(len(test)):
        model = ARIMA(history, order=(5,1,0))
        model_fit = model.fit(disp=0)
        output = model_fit.forecast()
        yhat = output[0][0]
        predictions.append(yhat)
        obs = test[t]
	history.append(yhat)
	#print('predicted=%f, expected=%f' % (yhat, obs))
    train.extend(predictions)
    plot_data(X, train)

def plot_data(origin, predict):
    df = pd.DataFrame()
    print len(origin), len(predict)
    error = mean_squared_error(origin, predict)
    print('Test MSE: %.3f' % error)
    df['original'] = origin
    df['predict'] = predict
    df.plot()
    plt.show()

def arima_predict(data,start, model):
    #predicts = []
    #for t in range(len(data)):
    #    predicts.append(model.forecast()[0])
    predicts = model.predict(start, start + len(data) - 2, typ='levels')
    print predicts, len(predicts)
    print data, len(data)
    plot_data(data, predicts)

if __name__ == '__main__':
    #get day from mysql orm to redis
    #print distinct_weekday_weekend('2015-07-01', '2016-10-31', '%Y-%m-%d')
    #connect redis
    client = connect_redis()
    #mysql_to_redis()
    res = []
    for x in range(1,2001):
        print "shop id is %d" % x
        day_users = client.hgetall(str(x))
        users = [ int(day_user[1]) for day_user in sorted_dict_values(day_users) ]
        #arima_train_broadcast(users)
        #split_pos = int(math.floor(len(users) * 0.66))
        split_pos = int(math.floor(len(users)))
        try:
            model = sarima_train(users[:split_pos])
            shop_res = sarima_predict([], split_pos, split_pos + 13, model)
            head = str(x) + ','
            res.append(head + ','.join(map(str,shop_res)))
        except ValueError:
            continue
    with open('result.csv', 'w') as wf:
        wf.write('\n'.join(res))

max_stars_repo_name: jiaorenyu/koubei

max_issues_repo_name: jiaorenyu/koubei

max_forks_repo_name: jiaorenyu/koubei

reason: import numpy,import statsmodels,from statsmodels


text:
from keras.applications.vgg16 import VGG16, preprocess_input, decode_predictions
from keras.preprocessing import image
import keras.backend as K
import numpy as np
import cv2
import sys

model = VGG16(weights=None)
model.load_weights('/home/f410/ZY1802122/vgg16_weights_tf_dim_ordering_tf_kernels.h5')
#img_path = sys.argv[1]
model.summary()
img_path='/home/f410/ZY1802122/relitu/009.png'
img = image.load_img(img_path, target_size=(224, 224))
x = image.img_to_array(img)
x = np.expand_dims(x, axis=0)
x = preprocess_input(x)

preds = model.predict(x)
class_idx = np.argmax(preds[0])
class_output = model.output[:, class_idx]
last_conv_layer = model.get_layer("block5_conv3")

grads = K.gradients(class_output, last_conv_layer.output)[0]
pooled_grads = K.mean(grads, axis=(0, 1, 2))
iterate = K.function([model.input], [pooled_grads, last_conv_layer.output[0]])
pooled_grads_value, conv_layer_output_value = iterate([x])
for i in range(512):
    conv_layer_output_value[:, :, i] *= pooled_grads_value[i]

heatmap = np.mean(conv_layer_output_value, axis=-1)
heatmap = np.maximum(heatmap, 0)
heatmap /= np.max(heatmap)
img = cv2.imread(img_path)
heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))
heatmap = np.uint8(255 * heatmap)
heatmap = cv2.applyColorMap(heatmap, cv2.COLORMAP_JET)
superimposed_img = cv2.addWeighted(img, 0.6, heatmap, 0.4, 0)
cv2.imshow("Original", img)
cv2.imshow("GradCam", superimposed_img)
cv2.waitKey(0)


text:
#!/usr/bin/env python

__author__ = "C. Clayton Violand"
__copyright__ = "Copyright 2017"

## Preps ImageNet input data (1000) for system by subtracting specified mean (image or pixel).
## Adds labels as first column.
## Writes result to file as csv.
##

from PIL import Image
import cv2
import numpy as np
import glob
import os
import sys

def load_image(infile):
    image = Image.open(infile)
    image.load()
    data = np.asarray(image, dtype="float")

    return data

def main(option):
    print "Taking " + option

    READ_DIR = 'helper/imagenet_tools/imagenet_1000/'
    LABEL_FILE = 'helper/imagenet_tools/imagenet_1000_labels.csv'

    if option == "mean_pixel":
        WRITE_FILE = 'helper/imagenet_tools/mean_pixel_reduced/imagenet_pxl_norm_1000.csv'
    else:
        MEAN_IMAGE = 'helper/imagenet_tools/imagenet_mean.npy'
        WRITE_FILE = 'helper/imagenet_tools/mean_image_reduced/imagenet_img_norm_1000.csv'

    # Load the mean image
    if option == "mean_pixel":
        pass
    else:
        mean_image = np.load(MEAN_IMAGE)
        mean_image = mean_image.astype(float)

    # Load the labels
    labels = np.genfromtxt(LABEL_FILE, dtype=int, delimiter=',')
    labels = labels[:,0]

    index = 0
    result = np.empty((0, 150529))
    for i in sorted(glob.glob(READ_DIR + '*.JPEG')):
        print i

        image = cv2.imread(i)
        image = np.moveaxis(image, 2, 0)
        image = image.astype(float)

        # Deal with grayscale images.
        if len(image.shape) == 2:
            image = np.stack((image, ) * 3)

        if option == "mean_pixel":
            image[0,:,:] = np.subtract(image[0,:,:], 103.939)
            image[1,:,:] = np.subtract(image[1,:,:], 116.779)
            image[2,:,:] = np.subtract(image[2,:,:], 123.68)
        else:
            image -= mean_image

        flattened_image = image.reshape(1, image.shape[0] * image.shape[1] * image.shape[2])
        flattened_image = np.insert(flattened_image, 0, labels[index]).reshape(1, 150529)
        result = np.append(result, flattened_image, axis=0)

        index += 1

    results = np.append(labels, result)
    np.savetxt(WRITE_FILE, result, fmt='%i', delimiter=',')

if __name__ == "__main__":
    try:
        option = sys.argv[1].lower()
    except IndexError:
        print "Usage: python helper/preprocess_imagenet.py <mean_image> | <mean_pixel>"
        sys.exit(1)

    if option != "mean_image" and option != "mean_pixel":
        print "Usage: python helper/preprocess_imagenet.py <mean_image> | <mean_pixel>"
        sys.exit(1)

    main(option)


text:
#!/usr/bin/env python
# coding: utf-8
#
#    Project: Azimuthal integration
#             https://github.com/silx-kit/pyFAI
#
#    Copyright (C) 2015-2018 European Synchrotron Radiation Facility, Grenoble, France
#
#    Principal author:       Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

"test suite for polarization corrections"

__author__ = "Jérôme Kieffer"
__contact__ = "Jerome.Kieffer@ESRF.eu"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
__date__ = "16/10/2020"

import unittest
import numpy
import logging
logger = logging.getLogger(__name__)
from ..azimuthalIntegrator import AzimuthalIntegrator

class TestPolarization(unittest.TestCase):

    def setUp(self):
        unittest.TestCase.setUp(self)
        self.shape = (13, 13)
        Y, X = numpy.ogrid[-6:7, -6:7]
        self.rotY = numpy.radians(30.0 * Y)
        self.rotX = numpy.radians(30.0 * X)
        self.tth = numpy.sqrt(self.rotY ** 2 + self.rotX ** 2)
        chi = numpy.arctan2(self.rotY, self.rotX)
        self.ai = AzimuthalIntegrator(dist=1, pixel1=0.1, pixel2=0.1)
        self.ai._cached_array["2th_center"] = self.tth
        self.ai._cached_array["chi_center"] = chi
        self.epsilon = 1e-15

    def tearDown(self):
        unittest.TestCase.tearDown(self)
        self.shape = self.rotY = self.rotX = self.tth = self.ai = None

    def testNoPol(self):
        "without polarization correction should be 1"
        self.assertTrue(abs(self.ai.polarization(factor=None) - numpy.ones(self.shape)).max() == 0, "without polarization correction should be 1")

    def testCircularPol(self):
        "Circular polarization should decay in (1+(cos2θ)^2)/2"
        pol = ((1.0 + numpy.cos(self.tth) ** 2) / 2.0).astype("float32")
        # print([abs(self.ai.polarization(factor=0, axis_offset=i) - pol).max() for i in range(6)])
        self.assertTrue(abs(self.ai.polarization(factor=0) - pol).max() == 0, "with circular polarization correction is independent of chi")
        self.assertTrue(abs(self.ai.polarization(factor=0, axis_offset=1) - pol).max() == 0, "with circular polarization correction is independent of chi, 1")
        self.assertTrue(abs(self.ai.polarization(factor=0, axis_offset=2) - pol).max() == 0, "with circular polarization correction is independent of chi, 2")
        self.assertTrue(abs(self.ai.polarization(factor=0, axis_offset=3) - pol).max() == 0, "with circular polarization correction is independent of chi, 3")

    def testHorizPol(self):
        "horizontal polarization should decay in (cos2θ)**2 in horizontal plane and no correction in vertical one"
        self.assertTrue(abs(self.ai.polarization(factor=1)[:, 6] - numpy.ones(13)).max() == 0, "No correction in the vertical plane")
        self.assertTrue(abs(self.ai.polarization(factor=1)[6] - numpy.cos(self.rotX) ** 2).max() < self.epsilon, "cos(2th)^2 like in the horizontal plane")

    def testVertPol(self):
        "Vertical polarization should decay in (cos2θ)**2 in vertical plane and no correction in horizontal one"
        self.assertTrue(abs(self.ai.polarization(factor=-1)[6] - numpy.ones(13)).max() == 0, "No correction in the horizontal plane")
        self.assertTrue(abs(self.ai.polarization(factor=-1)[:, 6] - (numpy.cos((2 * self.rotX)) + 1) / 2).max() < self.epsilon, "cos(2th)^2 like in the verical plane")

    def testoffsetPol(self):
        "test for the rotation of the polarization axis"
        self.assertTrue(abs(self.ai.polarization(factor=1, axis_offset=numpy.pi / 2)[6] - numpy.ones(13)).max() == 0, "No correction in the horizontal plane")
        self.assertTrue(abs(self.ai.polarization(factor=1, axis_offset=numpy.pi / 2)[:, 6] - (numpy.cos((2 * self.rotX)) + 1) / 2).max() < self.epsilon, "cos(2th)^2 like in the verical plane")

def suite():
    loader = unittest.defaultTestLoader.loadTestsFromTestCase
    testsuite = unittest.TestSuite()
    testsuite.addTest(loader(TestPolarization))
    return testsuite

if __name__ == '__main__':
    runner = unittest.TextTestRunner()
    runner.run(suite())

max_stars_repo_name: tacaswell/pyFAI

max_issues_repo_name: tacaswell/pyFAI

max_forks_repo_name: tacaswell/pyFAI


text:
import os.path
from collections import OrderedDict

import h5py

from .glm_matlab import glmnet_cv_best_result, convert_sklearn_kfold_to_glmnet_foldid
from scipy.stats import pearsonr
# from .fitting_util import num_fold_global
num_fold_global = 5
import numpy as np

def compute_all_r(y_true, y_predict, foldid):
    assert np.array_equal(np.unique(foldid), np.arange(num_fold_global) + 1)
    # print(y_true.shape, y_predict.shape, foldid.shape)
    assert y_true.shape == y_predict.shape == foldid.shape == (foldid.size,)
    r_all = []
    for fold in np.arange(num_fold_global):
        mask_this = (foldid == fold + 1)
        r_all.append(pearsonr(y_true[mask_this],
                              y_predict[mask_this])[0])
    r_all = np.array(r_all)
    return r_all

def one_train_loop_cv(X, y, cv_seed, training_tool_params):
    assert training_tool_params.keys() >= {'reg', 'family'}
    family = training_tool_params['family']
    assert family in {'gaussian', 'poisson', 'softplus'}
    reg = training_tool_params['reg']
    if reg == 'L1':
        alpha_value = 1.0
        standardize = False
    elif reg == 'L2':
        alpha_value = 0.0
        standardize = False
    # add these to try more variations on regularization.
    elif reg == 'L1_ns':
        alpha_value = 1.0
        standardize = True
    elif reg == 'L2_ns':
        alpha_value = 0.0
        standardize = True
    else:
        raise NotImplementedError('wrong value {}'.format(reg))
    foldid_to_use = convert_sklearn_kfold_to_glmnet_foldid(len(X), num_fold_global, seed=cv_seed)
    # then pass in everything.
    # I think standardize = False makes sense, as scale in PCAed data carry information

    y_to_return = glmnet_cv_best_result(X, y, standardize=standardize,
                                        family=family,
                                        foldid=foldid_to_use, alpha=alpha_value)

    # # then let's compute mean r2.
    # if family == 'poisson':
    #     y_to_return = np.exp(y_this_predicted_debug)
    # elif family == 'gaussian':
    #     y_to_return = y_this_predicted_debug
    # else:
    #     raise NotImplementedError('no such family')
    r_all = compute_all_r(y.ravel(), y_to_return, foldid_to_use)
    return y_to_return, r_all

reason: import numpy,from scipy


text:
from math import sqrt
from numpy import linspace
from scipy import optimize
from scipy.linalg import norm

def fit(p, *n):
    if n == ():
        n = 1000
    n_reparam = 1000
    x_data = []
    y_data = []
    z_data = []
    d = []
    y_reparam = []
    z_reparam = []
    u_reparam = []
    u_coeffs = []
    u_scaled = linspace(0, 1, n_reparam)
    ut = linspace(0, 1, n)
    chord = []
    du = ut[1] - ut[0]
    parameter = []

    for i in range(len(p)):
        pt = p[i]
        x_data.append(pt[0])
        y_data.append(pt[1])
        z_data.append(pt[2])

    for i in range(1, len(p)):
        chord.append(sqrt((z_data[i]-z_data[i-1])**2+(y_data[i]-y_data[i-1])**2))

    chord_arc = sum(chord)

    for i in range(len(chord)):
        parameter.append(sum(chord[:i])/chord_arc)
    parameter.append(1)

    y_coeffs, y_cov = optimize.curve_fit(y_func, parameter, y_data)
    z_coeffs, z_cov = optimize.curve_fit(z_func, parameter, z_data)

    for i in range(n):
        par = ut[i]

        yen = y_func(par + du, *y_coeffs)
        ye = y_func(par, *y_coeffs)
        zen = z_func(par + du, *z_coeffs)
        ze = z_func(par, *z_coeffs)

        d.append(sqrt((yen-ye)**2+(zen-ze)**2))
    perimeter = sum(d)

    if y_data[0] < y_data[-1]:
        direction = 1
    else:
        direction = -1

    for i in range(n_reparam):
        y_reparam.append(y_func(u_scaled[i], *y_coeffs))
        z_reparam.append(z_func(u_scaled[i], *z_coeffs))

    for i in range(1, n_reparam):
        u_reparam.append(sqrt((z_reparam[i]-z_reparam[i-1])**2+(y_reparam[i]-y_reparam[i-1])**2))

    u_chord_arc = sum(u_reparam)

    for i in range(len(u_reparam)):
        u_coeffs.append(sum(u_reparam[:i]))
    u_coeffs.append(u_chord_arc)

    u_coeffs, u_cov = optimize.curve_fit(u_func, u_coeffs, u_scaled)

    return perimeter, x_data[0], [y_coeffs, z_coeffs, u_coeffs, direction]

def y_func(x, a, b, c, d):
    return a * x**3 + b * x**2 + c * x + d

def y_func_diff(x, a, b, c, d):
    return 3 * a * x**2 + 2 * b * x + c

def z_func(x, a, b, c, d, e):
    return a * x**4 + b * x**3 + c * x**2 + d * x + e

def z_func_diff(x, a, b, c, d, e):
    return 4 * a * x**3 + 3 * b * x**2 + 2 * c * x + d

def u_func(x, a, b, c, d, e, f):
    return a * x**5 + b * x**4 + c * x**3 + d * x**2 + e * x + f

def surf(uv_list, params, offset):
    u = uv_list[0]
    v = uv_list[1]
    y_coeffs = params[0]
    z_coeffs = params[1]
    u_coeffs = params[2]
    direction = params[3]
    x = []
    y = []
    z = []
    nx = []
    ny = []
    nz = []

    for i in range(len(u)):
        par = u_func(u[i], *u_coeffs)
        x.append(v[i] + offset)
        y.append(y_func(par, *y_coeffs))
        z.append(z_func(par, *z_coeffs))

        x_diff = 0
        y_diff = y_func_diff(par, *y_coeffs)
        z_diff = z_func_diff(par, *z_coeffs)

        scale = -direction*norm([x_diff, y_diff, z_diff])

        nx.append(0)
        ny.append(z_diff/scale)
        nz.append(-y_diff/scale)

    return [x, y, z, nx, ny, nz]

reason: from numpy,from scipy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 13 16:35:12 2017

@author: smussend
"""
# An exercise in using pandas to analyze data, using the 2016 Maryland state employees salary data, via Baltimore Sun. Data link: https://github.com/baltimore-sun-data/maryland-salaries. Interactive database link: http://data.baltimoresun.com/salaries/state/cy2016/.

# Do your best to work with pandas as much as possible here, as it's the standard for working with tabular data. The DataCamp lesson due next week goes into detail on pandas. Keep in mind that may be multiple ways to answer the question; find one that works!

# Four resources you may find useful, in addition to DataCamp:
# Pandas cheat sheet: https://github.com/pandas-dev/pandas/blob/master/doc/cheatsheet/Pandas_Cheat_Sheet.pdf
# Pandas documentation: https://pandas.pydata.org/pandas-docs/stable/basics.html
# Pandas tutorials: https://pandas.pydata.org/pandas-docs/stable/tutorials.html
# And, of course, Google

# Important note about the data: though this data was published by the Baltimore Sun, it's not totally clean.  Typically, before beginning our analysis, we'd walk through a series of steps to prepare the data, searching for flaws that could compromise our analysis.  And as we analyze the data and uncover other problems, we'd update our cleaning scripts. Because the purpose of this exercise is to get a little more familiar with the basics of pandas -- and not an exercise in data cleaning -- we're going to work with it even though we haven't done a robust cleaning yet.  We'll do that in future exercises.

# First, import the numpy and pandas packages.

# Use pandas to read in the data and store the dataframe as a variable you can call later. Note, reading in the data may fail if you haven't set the working directory.

# Question 1: Can you produce a list of all the column names? column_names = df.columns.values

# Question 2: What was the highest YTD_gross_earnings for an employee in 2016? Hint: you can use one short bit of code to get the answers to question 2 through 5. df['ytd_gross_earnings'].max()

# Question 3: What was the lowest YTD_gross_earnings for an employee in 2016? df['ytd_gross_earnings'].min()

# Question 4: What was the average (mean) salary (YTD_gross_earnings) in 2016? df['ytd_gross_earnings'].mean()

# Question 5: What was the median (50th percentile) salary (YTD_gross_earnings) in 2016? 35,000

# Question 6: How much did the state spend on employee salaries (YTD_gross_earnings) in 2016? 5672042000

# Question 7: Can you produce a list of all the organizations represented in the data set? After you generate it, closely review the list.  What problem do you see that might hinder your ability to ask questions about some organizations?  Think about how you might correct the problem you've identified (no need to write code at this point, just think about it).

# Question 8: Which 10 employees got the most money (YTD_gross_earnings) from a single organization in 2016, and how much did each of them make? What do they have in common?

# Question 9: Who was the highest paid non-University of Maryland employee?

# Question 10: What was the mean salary (YTD_gross_earnings) for University of Maryland employees? How does that compare to the mean for St. Mary’s College of Maryland employees?

# Question 11: What other questions can you ask of this dataset?

import pandas as pd
import numpy as np

df = pd.read_csv('../data/state-cy-2016.csv')

column_names = df.columns.values

# df['ytd_gross_earnings'].sum()

organizations = pd.DataFrame(df.organization.unique())
organizations

df.sort_values('ytd_gross_earnings',

max_stars_repo_name: kvasudevan/pythonclass-

max_issues_repo_name: kvasudevan/pythonclass-

max_forks_repo_name: kvasudevan/pythonclass-


text:
import numpy as np
import os
import matplotlib.pyplot
import time

time0 = time.time()

# # Append mapparam to pchain
# cmd='cp ../Step1_ParametricFit/Gumbel/pchain.dat .'
# os.system(cmd)
# cmd='cp ../Step1_ParametricFit/Gumbel/mapparam.dat .'
# os.system(cmd)

cmd = "./transpose_file.x mapparam.dat >> pchain.dat"
os.system(cmd)

# Run model inference code in a non-inference, propagation-only regime
# Parameters have to be very carefully selected to match them with weibull_inference code
METHOD="gausmarg"
DATANOISE=0.0
cmd = 'echo "0" > rndInd_file.dat'
os.system(cmd)
cmd = 'echo "1" >> rndInd_file.dat'
os.system(cmd)

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

std_ModErr, std_PusFor, std_ModErrMAP = np.loadtxt('pvars.txt', usecols = (0,1,2), unpack=True)
p, p_MAP = np.loadtxt('pmeans.txt', usecols = (0,1,), unpack=True)

std_p0 = np.sqrt(std_ModErr[0]+std_PusFor[0])
std_p1 = np.sqrt(std_ModErr[1]+std_PusFor[1])

amid = p[0]
bmid = p[1]
adel = 3.0*std_p0
bdel = 3*std_p1

print "amid = ", amid
print "adel = ", adel
print "bmid = ", bmid
print "bdel = ", bdel

cmd="awk '{print ($1-amid)/adel,$2/adel,($3-bmid)/bdel,$4/bdel,$5/bdel}' amid=" + str(amid) + " adel=" + str(adel) + " bmid=" + str(bmid) + " bdel=" + str(bdel) +  " pchain.dat > pchain_sc.dat"
os.system(cmd)

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

# cmd="awk '{print ($1-amid)/adel,$2/adel,($3-bmid)/bdel,$4/bdel,$5/bdel}' amid=1.0 adel=1.0 bmid=1.0 bdel=0.4 pchain.dat > pchain_sc.dat"
# os.system(cmd)

cmd = "$UQBINDIR/model_inf -z -x Depth.dat -y Depth.dat -f pc -p pchain_sc.dat -l" + str(METHOD) + " -e " + str(DATANOISE) + " -m0 -o 3 -s pct -d 2 -t Depth.dat -r rndInd_file.dat"
os.system(cmd)

# Postprocess
cmd = "paste Depth.dat fmeans.dat fvars.dat > full_propagation.dat"
os.system(cmd)

cmd = "cp full_propagation.dat ../Step4_Postprocessing/full_propagation.dat"
os.system(cmd)

# Clean the rest
cmd = "rm -rf *.dat *.txt"
os.system(cmd)

time1 = time.time()

print "Time elapsed = ", time1 - time0, " sec"

max_stars_repo_name: ORNL-Fusion/xolotl

max_issues_repo_name: ORNL-Fusion/xolotl

max_forks_repo_name: ORNL-Fusion/xolotl


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

"""Tests for the NumPy ops."""

import itertools
from typing import Callable
from absl.testing import parameterized

import numpy as np
from scipy import special

import tensorflow as tf

from dedal import alignment
from dedal import smith_waterman_np as npy_ops

# For test purposes.
def _sw_general(sim_mat,
                gap_func):
  """Computes soft Smith-Waterman with general gap function.

  Args:
    sim_mat: a np.ndarray<float>[len1, len2] containing
     the substitution values for pairs of sequences.
    gap_func: a function of the form gap_func(k), where k is an integer.

  Returns:
    Smith-Waterman value (float).
  """
  len_1, len_2 = sim_mat.shape

  match = np.zeros((len_1 + 1, len_2 + 1))

  for i in range(1, len_1 + 1):
    for j in range(1, len_2 + 1):
      tmp = match[i - 1, j - 1] + sim_mat[i-1, j-1]
      delete = max([match[i - k, j] - gap_func(k) for k in range(1, i + 1)])
      insert = max([match[i, j - k] - gap_func(k) for k in range(1, j + 1)])
      match[i, j] = max(tmp, delete, insert, 0)

  return max(match[1:, 1:].ravel())

def scores_brute_force(weights):
  len_1, len_2, _ = weights.shape
  ret = []
  for alignment_mat in npy_ops.alignment_matrices(len_1, len_2):
    ret.append(np.vdot(alignment_mat, weights))
  return np.array(ret)

def max_brute_force(weights):
  scores = scores_brute_force(weights)
  return np.max(scores)

def lse_brute_force(weights, temperature=1.0):
  scores = scores_brute_force(weights)
  return temperature * special.logsumexp(scores / temperature)

class SmithWatermanNumpyTest(parameterized.TestCase, tf.test.TestCase):

  def setUp(self):
    super().setUp()

    single_sub = 6*tf.eye(5) - 5*tf.ones((5, 5))
    second_sub = tf.tensor_scatter_nd_update(single_sub,
                                             indices=[[0, 0], [1, 1]],
                                             updates=[-5, -5])
    third_sub = tf.tensor_scatter_nd_update(single_sub,
                                            indices=[[0, 0], [2, 2]],
                                            updates=[-5, -5])
    fourth_sub = tf.tensor_scatter_nd_update(single_sub,
                                             indices=[[0, 0], [4, 4]],
                                             updates=[-5, -5])
    self._toy_sub = tf.stack([single_sub, second_sub, third_sub, fourth_sub])
    self._toy_gap_open = 0.03 * tf.ones((self._toy_sub.shape[0],))
    self._toy_gap_extend = 0.02 * tf.ones((self._toy_sub.shape[0],))
    self._weights = alignment.weights_from_sim_mat(self._toy_sub,
                                                   self._toy_gap_open,
                                                   self._toy_gap_extend)

  def test_max_against_brute_force(self):
    maxes = [max_brute_force(self._weights[i])
             for i in range(self._weights.shape[0])]
    maxes2 = npy_ops.soft_sw_affine(self._toy_sub.numpy(),
                                    self._toy_gap_open.numpy(),
                                    self._toy_gap_extend.numpy(),
                                    temperature=0)
    self.assertAllClose(maxes, maxes2)

  def test_max_against_sw_general(self):
    def make_gap_func(gap_open, gap_extend):
      def gap_func(k):
        return gap_open + gap_extend * (k-1)
      return gap_func

    toy_sub = self._toy_sub.numpy()
    toy_gap_open = self._toy_gap_open.numpy()
    toy_gap_extend = self._toy_gap_extend.numpy()

    maxes = [_sw_general(toy_sub[i],
                         make_gap_func(toy_gap_open[i], toy_gap_extend[i]))
             for i in range(self._weights.shape[0])]
    maxes2 = npy_ops.soft_sw_affine(toy_sub, toy_gap_open, toy_gap_extend,
                                    temperature=0)
    self.assertAllClose(maxes, maxes2)

  def test_lse_against_brute_force(self):
    for temp in (1e-3, 1.0, 1e3):
      softmaxes = np.array([lse_brute_force(self._weights[i],
                                            temperature=temp)
                            for i in range(self._weights.shape[0])])
      softmaxes2 = npy_ops.soft_sw_affine(self._toy_sub.numpy(),
                                          self._toy_gap_open.numpy(),
                                          self._toy_gap_extend.numpy(),
                                          temperature=temp)
      self.assertAllClose(softmaxes, softmaxes2)

  @parameterized.parameters(itertools.product([0, 0.1, 1, 10], range(4)))
  def test_gradients_against_finite_difference(self, temp=1.0, j=0):
    def num_grad(f, x, eps=1e-4):
      shape = x.shape
      x_flat = x.ravel()
      grad_flat = np.zeros_like(x_flat)
      for j in range(len(x_flat)):
        e = np.zeros_like(x_flat)
        e[j] = 1
        call1 = f((x_flat + eps * e).reshape(shape))
        call2 = f((x_flat - eps * e).reshape(shape))
        grad_flat[j] = (call1 - call2) / (2 * eps)
      return grad_flat.reshape(shape)

    toy_sub = self._toy_sub.numpy()
    toy_gap_open = self._toy_gap_open.numpy()
    toy_gap_extend = self._toy_gap_extend.numpy()

    _, g_sim_mat, g_gap_open, g_gap_extend = (
        npy_ops._soft_sw_affine(sim_mat=toy_sub[j], gap_open=toy_gap_open[j],
                                gap_extend=toy_gap_extend[j], temperature=temp,
                                ret_grads=True))
    # Gradient w.r.t. sim_mat.
    def f_sim_mat(x):
      return npy_ops._soft_sw_affine(sim_mat=x, gap_open=toy_gap_open[j],
                                     gap_extend=toy_gap_extend[j],
                                     temperature=temp, ret_grads=False)

    g_sim_mat_num = num_grad(f_sim_mat, toy_sub[j].astype(np.float64))
    self.assertAllClose(g_sim_mat, g_sim_mat_num)

    # Gradient w.r.t. gap_open.
    def f_gap_open(x):
      return npy_ops._soft_sw_affine(sim_mat=toy_sub[j], gap_open=x,
                                     gap_extend=toy_gap_extend[j],
                                     temperature=temp, ret_grads=False)

    g_gap_open_num = num_grad(f_gap_open, toy_gap_open[j].astype(np.float64))
    self.assertAllClose(g_gap_open, g_gap_open_num)

    # Gradient w.r.t. gap_extend.
    def f_gap_ext(x):
      return npy_ops._soft_sw_affine(sim_mat=toy_sub[j],
                                     gap_open=toy_gap_open[j], gap_extend=x,
                                     temperature=temp, ret_grads=False)

    g_gap_extend_num = num_grad(f_gap_ext, toy_gap_extend[j].astype(np.float64))
    self.assertAllClose(g_gap_extend, g_gap_extend_num)

  def test_gradient_shape(self):
    toy_sub = self._toy_sub.numpy()
    toy_gap_open = self._toy_gap_open.numpy()
    toy_gap_extend = self._toy_gap_extend.numpy()

    values, g_sim_mat, g_gap_open, g_gap_extend = npy_ops.soft_sw_affine(
        sim_mat=toy_sub,
        gap_open=toy_gap_open,
        gap_extend=toy_gap_extend,
        ret_grads=True)
    self.assertEqual(values.shape[0], toy_sub.shape[0])
    self.assertAllEqual(g_sim_mat.shape, toy_sub.shape)
    self.assertAllEqual(g_gap_open.shape, toy_gap_open.shape)
    self.assertAllEqual(g_gap_extend.shape, toy_gap_extend.shape)

if __name__ == '__main__':
  tf.test.main()

max_stars_repo_name: ParikhKadam/google-research

max_issues_repo_name: ParikhKadam/google-research

max_forks_repo_name: admariner/google-research