text:
#!/usr/bin/env python
PROGRAM_NAME = "cg-get-site-bw-avg.py"
PROGRAM_DESCRIPTION = """
CloudGenix Get-Site-BW Daily Averages Script
---------------------------------------
Writes the past 7-days bandwidth average dailies and the 7-day average for each SPOKE site.
Output is placed in a CSV file

optional arguments:
  -h, --help            show this help message and exit
  --token "MYTOKEN", -t "MYTOKEN"
                        specify an authtoken to use for CloudGenix
                        authentication
  --authtokenfile "MYTOKENFILE.TXT", -f "MYTOKENFILE.TXT"
                        a file containing the authtoken
  --csvfile csvfile, -c csvfile
                        the CSV Filename to write the BW averages to

"""

####Library Imports
from cloudgenix import API, jd
import os
import sys
import argparse
from csv import reader

def parse_arguments():
    CLIARGS = {}
    parser = argparse.ArgumentParser(
        prog=PROGRAM_NAME,
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description=PROGRAM_DESCRIPTION
            )
    parser.add_argument('--token', '-t', metavar='"MYTOKEN"', type=str,
                    help='specify an authtoken to use for CloudGenix authentication')
    parser.add_argument('--authtokenfile', '-f', metavar='"MYTOKENFILE.TXT"', type=str,
                    help='a file containing the authtoken')
    parser.add_argument('--csvfile', '-c', metavar='csvfile', type=str,
                    help='the CSV Filename to write the BW averages to', required=True)
    args = parser.parse_args()
    CLIARGS.update(vars(args))
    return CLIARGS

def authenticate(CLIARGS):
    print("AUTHENTICATING...")
    user_email = None
    user_password = None

    sdk = API()
    ##First attempt to use an AuthTOKEN if defined
    if CLIARGS['token']:                    #Check if AuthToken is in the CLI ARG
        CLOUDGENIX_AUTH_TOKEN = CLIARGS['token']
        print("    ","Authenticating using Auth-Token in from CLI ARGS")
    elif CLIARGS['authtokenfile']:          #Next: Check if an AuthToken file is used
        tokenfile = open(CLIARGS['authtokenfile'])
        CLOUDGENIX_AUTH_TOKEN = tokenfile.read().strip()
        print("    ","Authenticating using Auth-token from file",CLIARGS['authtokenfile'])
    elif "X_AUTH_TOKEN" in os.environ:              #Next: Check if an AuthToken is defined in the OS as X_AUTH_TOKEN
        CLOUDGENIX_AUTH_TOKEN = os.environ.get('X_AUTH_TOKEN')
        print("    ","Authenticating using environment variable X_AUTH_TOKEN")
    elif "AUTH_TOKEN" in os.environ:                #Next: Check if an AuthToken is defined in the OS as AUTH_TOKEN
        CLOUDGENIX_AUTH_TOKEN = os.environ.get('AUTH_TOKEN')
        print("    ","Authenticating using environment variable AUTH_TOKEN")
    else:                                           #Next: If we are not using an AUTH TOKEN, set it to NULL
        CLOUDGENIX_AUTH_TOKEN = None
        print("    ","Authenticating using interactive login")
    ##ATTEMPT AUTHENTICATION
    if CLOUDGENIX_AUTH_TOKEN:
        sdk.interactive.use_token(CLOUDGENIX_AUTH_TOKEN)
        if sdk.tenant_id is None:
            print("    ","ERROR: AUTH_TOKEN login failure, please check token.")
            sys.exit()
    else:
        while sdk.tenant_id is None:
            sdk.interactive.login(user_email, user_password)
            # clear after one failed login, force relogin.
            if not sdk.tenant_id:
                user_email = None
                user_password = None
    print("    ","SUCCESS: Authentication Complete")
    return sdk

def logout(sdk):
    print("Logging out")
    sdk.get.logout()

##########MAIN FUNCTION#############
def go(sdk, CLIARGS):
    ####CODE GOES BELOW HERE#########
    csvfilename = CLIARGS['csvfile']
    resp = sdk.get.tenants()
    if resp.cgx_status:
        tenant_name = resp.cgx_content.get("name", None)
        print("======== TENANT NAME",tenant_name,"========")
    else:
        logout()
        print("ERROR: API Call failure when enumerating TENANT Name! Exiting!")
        print(resp.cgx_status)
        sys.exit((vars(resp)))

    result = sdk.get.sites()
    if not result.cgx_status:
        sys.exit("Error Getting Site list")
    site_list = result.cgx_content.get("items")

    parsed_sites = []
    for site in site_list:
        if site['element_cluster_role'] == "SPOKE":
            parsed_sites.append([site['name'], site['id']])

    csv_out_array = [['SiteName', 'SiteID', '24Hr-Today', '24Hr-Yesterday',  '24Hr-3-days-ago', '24Hr-4-days-ago', '24Hr-5-days-ago', '24Hr-6-days-ago', '24Hr-7-days-ago',  '7-Day-AVG' ]]
    for site_data in parsed_sites:
        site_name = site_data[0]
        site_id = site_data[1]
        sum = 0
        csv_row = []
        csv_row.append(site_name) #1
        csv_row.append(site_id) #2
        for days_back in range(0,7): #7 days worth of data
            (start_time, end_time) = cgx_generate_timestamps_days(days_interval=1, offset_days=days_back)
            day_bw = cgx_get_bw_consumption(sdk, start_time, end_time, site_id=site_id)
            csv_row.append( day_bw ) # 3-10
            sum += day_bw
        sum = ( sum / 7 ) #7-day average
        csv_row.append( round(sum,2) ) # 11
        csv_out_array.append(csv_row) ## add the sites row

    result = False
    result = write_2d_list_to_csv(csvfilename, csv_out_array)
    if result:
        print("Wrote CSV File",csvfilename,"Successfully")
    else:
        print("Failure Writing CSV File")
    ####CODE GOES ABOVE HERE#########

#/----------------------
#| cgx_generate_timestamps_days - Generates CGX start and end Timestamps used in events, alarms, reporting for specified interval in days (default 1 day or 24 hours)
def cgx_generate_timestamps_days(days_interval=1, offset_days=0):
    from datetime import timedelta
    from datetime import datetime
    now = (datetime.utcnow() - timedelta(days=(offset_days)))
    end_time = now.isoformat()
    start_time = (now - timedelta(days=(days_interval))).isoformat()
    return (str(start_time)+"Z", str(end_time)+"Z") ## return (start_time, end_time) I.E. call with "(start_time, end_time) = cgx_generate_timestamps(days_interval=1) for the last 24 hours"
#\----------------------

#/----------------------
#| validate_2d_array - Validates an array is 2-dimensional for use in CSV Export
def validate_2d_array(test_list):
    import numpy as np
    np_list = np.array(test_list)
    if len(np_list.shape) == 2:
        return True
    return False
#\----------------------

#/----------------------
#| write_2d_list_to_csv - Writes a 2-Dimensional list to a CSV file
def write_2d_list_to_csv(csv_file, list_2d, write_mode="w"):
    import csv
    try:
        file = open(csv_file, write_mode)
        with file:
            write = csv.writer(file)
            write.writerows(list_2d)
            return True
        return False
    except:
        return False

#\----------------------

#/----------------------
#| cgx_average_series - takes a metrics series structure (input['metrics']['series']['data']['datapoints'][**list**]['value']) and averages it to the decimal places (default:2)
#|                  A typical call would look like:
#|                      metrics = cgx_get_bw_consumption(sdk,start,end) ## Calls sdk.post.monitor_metrics and returns Metrics with Series in them
#|                      for series in metrics.get("metrics",[{}])[0].get("series",[]):
#|                          average = cgx_average_series(series)
def cgx_average_series(metrics_series_structure, decimal_places=2):
    count = 0
    sum = 0
    for datapoints in metrics_series_structure.get("data",[{}])[0].get("datapoints",[{}]):
        if (datapoints.get("value",None) is not None):
            count += 1
            sum += datapoints.get("value",0)
    if count != 0:
        if decimal_places == 0:
            return round((sum/count))
        return round((sum/count),decimal_places)
    return 0
#\----------------------

#/----------------------
#| cgx_get_bw_consumption - Gets bandwidth consumption average for a give time period. If no site_id is given, Aggregate BW consumption for the tenant is provided.
def cgx_get_bw_consumption(sdk, start_time, end_time, site_id=None):
    true = True
    false = False
    post_request = {"start_time":start_time, "end_time":end_time, "interval":"5min","metrics":[{"name":"BandwidthUsage","statistics":["average"],"unit":"Mbps"}],"view":{},"filter":{"site":[]}}
    if site_id: post_request['filter']['site'].append(str(site_id))
    result = sdk.post.monitor_metrics(post_request)
    metrics = result.cgx_content
    series = metrics.get("metrics",[{}])[0].get("series",[])[0]
    return(cgx_average_series(series))
#\----------------------

if __name__ == "__main__":
    ###Get the CLI Arguments
    CLIARGS = parse_arguments()

    ###Authenticate
    SDK = authenticate(CLIARGS)

    ###Run Code
    go(SDK, CLIARGS)

    ###Exit Program
    logout(SDK)

ext: py

max_stars_repo_name: karl-sc/cg-get-site-bw-avg

max_issues_repo_name: karl-sc/cg-get-site-bw-avg

max_forks_repo_name: karl-sc/cg-get-site-bw-avg

reason: import numpy


text:
#!/usr/bin/env python

#python 3 compatibility
from __future__ import print_function

#stdlib imports
import sys
import collections
import datetime
import time
import os.path

#third party imports
import h5py
from scipy.io import netcdf
import numpy as np
from .multiple import MultiGrid
from .shake import ShakeGrid
from .grid2d import Grid2D
from .dataset import DataSetException
from .geodict import GeoDict

class MultiHazardGrid(MultiGrid):
    def __init__(self,layers,geodict,origin,header,metadata=None):
        """Construct a MultiHazardGrid object.

        :param layers:
          OrderedDict containing earthquake-induced hazard data layers (keys are 'pga', etc., values are 2D arrays of data).
        :param geodict:
          GeoDict object specifying the spatial extent,resolution and shape of the data.
        :param origin:
          Dictionary with elements:
            - id String of event ID (i.e., 'us2015abcd')
            - source String containing originating network ('us')
            - time Python datetime of origin time, UTC.
            - lat Float event latitude
            - lon Float event longitude
            - depth Float event depth
            - magnitude Datetime object representing event origin time.
        :param header:
          Dictionary with elements:
            - type Type of multi-layer earthquake induced hazard ('shakemap','gfe')
            - version Integer product version (1)
            - process_time Python datetime indicating when data was created (UTC)
            - code_version String version of code that created this file (i.e.,'4.0')
            - originator String representing network that created the hazard grid.
            - product_id String representing the ID of the product (may be different from origin ID)
            - map_status String, one of RELEASED, ??
            - event_type String, one of ['ACTUAL','SCENARIO']
        :param metadata:
          Dictionary of dictionaries containing other metadata users wish to preserve.
        :returns: A MultiHazardGrid object.
        """
        self._layers = collections.OrderedDict()
        self._geodict = geodict
        for (layerkey,layerdata) in layers.items():
            self._layers[layerkey] = Grid2D(data=layerdata,geodict=geodict)

        self.setHeader(header)
        self.setOrigin(origin)
        self.setMetadata(metadata)

    @classmethod
    def getFileGeoDict(cls,filename):
        """Get the spatial extent, resolution, and shape of grid inside MultiHazardGrid format HDF file.

        :param filename:
           File name of HDF file.
        :returns:
          GeoDict specifying spatial extent, resolution, and shape of grid inside HDF file.
        """
        f = h5py.File(filename, "r")
        geodict = {}
        xvar = f['x'][:]
        yvar = f['y'][:]
        geodict['xmin'] = xvar[0]
        geodict['xmax'] = xvar[-1]
        geodict['ymin'] = yvar[0]
        geodict['ymax'] = yvar[-1]
        geodict['ny'] = len(yvar)
        geodict['nx'] = len(xvar)
        geodict['dx'] = xvar[1]-xvar[0]
        geodict['dy'] = yvar[1]-yvar[0]
        f.close()

        return GeoDict(geodict)

    def _saveDict(self,group,mydict):
        """
        Recursively save dictionaries into groups in an HDF file.
        :param group:
          HDF group object to contain a given dictionary of data in HDF file.
        :param mydict:
          Dictionary of values to save in group.  Dictionary can contain objects of the following types:
            - str,unicode,int,float,long,list,tuple,np.ndarray,dict,datetime.datetime,collections.OrderedDict
        """
        ALLOWED = [str,int,float,
                   type(None),
                   list,tuple,np.ndarray,
                   dict,datetime.datetime,
                   collections.OrderedDict]
        if sys.version_info.major == 2:
            ALLOWED.append(unicode)
            ALLOWED.append(long)

        for (key,value) in mydict.items():
            tvalue = type(value)
            if tvalue not in ALLOWED:
                raise DataSetException('Unsupported metadata value type "%s"' % tvalue)
            if not isinstance(value,dict):
                if isinstance(value,datetime.datetime):
                    value = time.mktime(value.timetuple())
                group.attrs[key] = value
            else:
                subgroup = group.create_group(key)
                self._saveDict(subgroup,value)

    @classmethod
    def _loadDict(cls,group):
        """Recursively load dictionaries from groups in an HDF file.

        :param group:
          HDF5 group object.
        :returns:
          Dictionary of metadata (possibly containing other dictionaries).
        """
        tdict = {}
        for (key,value) in group.attrs.items(): #attrs are NOT subgroups
            if key.find('time') > -1:
                value = value = datetime.datetime.utcfromtimestamp(value)
            tdict[key] = value
        for (key,value) in group.items(): #these are going to be the subgroups
            tdict[key] = cls._loadDict(value)
        return tdict

    def save(self,filename):
        """
        Save MultiHazardGrid object to HDF file.
        Georeferencing information will be saved as datasets "x" and "y".  Layers will be saved as
        datasets named by layer keys.  Dictionaries contained in "origin", and "header" will be saved in
        groups of those same names.  Dictionaries contained in the "metadata" dictionary will be contained
        in a series of recursive groups under a group called "metadata".

        :param filename:
          Output desired filename (HDF format).
        """
        f = h5py.File(filename, "w")
        #Add in some attributes that will help make this GMT friendly...
        f.attrs['Conventions'] = 'COARDS, CF-1.5'
        f.attrs['title'] = 'filename'
        f.attrs['history'] = 'Created with python MultiHazardGrid.save(%s)' % filename
        f.attrs['GMT_version'] = 'NA'

        #create two top-level groups that should always be present
        header = f.create_group('header')
        self._saveDict(header,self._header)

        origin = f.create_group('origin')
        self._saveDict(origin,self._origin)

        #write out any other metadata, creating groups recursively as needed
        metadata = f.create_group('metadata')
        self._saveDict(metadata,self._metadata)

        xvar = np.linspace(self._geodict.xmin,self._geodict.xmax,self._geodict.nx)
        yvar = np.linspace(self._geodict.ymin,self._geodict.ymax,self._geodict.ny)
        x = f.create_dataset('x',data=xvar,compression='gzip',shape=xvar.shape,dtype=str(xvar.dtype))
        x.attrs['CLASS'] = 'DIMENSION_SCALE'
        x.attrs['NAME'] = 'x'
        x.attrs['_Netcdf4Dimid'] = 0 #no idea what this is
        x.attrs['long_name'] = 'x'
        x.attrs['actual_range'] = np.array((xvar[0],xvar[-1]))

        y = f.create_dataset('y',data=yvar,compression='gzip',shape=yvar.shape,dtype=str(yvar.dtype))
        y.attrs['CLASS'] = 'DIMENSION_SCALE'
        y.attrs['NAME'] = 'y'
        y.attrs['_Netcdf4Dimid'] = 1 #no idea what this is
        y.attrs['long_name'] = 'y'
        y.attrs['actual_range'] = np.array((yvar[0],yvar[-1]))

        for (layerkey,layer) in self._layers.items():
            dset = f.create_dataset(layerkey,data=layer.getData(),compression='gzip')
            dset.attrs['long_name'] = layerkey
            dset.attrs['actual_range'] = np.array((np.nanmin(layer._data),np.nanmax(layer._data)))

        f.close()

    @classmethod
    def load(cls,filename):
        """
        Load data from an HDF file into a MultiHazardGrid object.

        :param filename:
          HDF file containing data and metadata for ShakeMap or Secondary Hazards data.
        :returns:
          MultiHazardGrid object.
        """
        f = h5py.File(filename, "r")
        REQUIRED_GROUPS = ['origin','header']
        REQUIRED_DATASETS = ['x','y']
        for group in REQUIRED_GROUPS:
            if group not in f.keys():
                raise DataSetException('Missing required metadata group "%s"' % group)
        for dset in REQUIRED_DATASETS:
            if dset not in f.keys():
                raise DataSetException('Missing required data set "%s"' % dset)

        header = {}
        for (key,value) in f['header'].attrs.items():
            if key.find('time') > -1:
                value = datetime.datetime.utcfromtimestamp(value)
            header[key] = value

        origin = {}
        for (key,value) in f['origin'].attrs.items():
            if key.find('time') > -1:
                value = datetime.datetime.utcfromtimestamp(value)
            origin[key] = value

        if 'metadata' in f.keys():
            metadata = cls._loadDict(f['metadata'])

        geodict = {}
        xvar = f['x'][:]
        yvar = f['y'][:]
        geodict['xmin'] = xvar[0]
        geodict['xmax'] = xvar[-1]
        geodict['ymin'] = yvar[0]
        geodict['ymax'] = yvar[-1]
        geodict['ny'] = len(yvar)
        geodict['nx'] = len(xvar)
        geodict['dx'] = xvar[1]-xvar[0]
        geodict['dy'] = yvar[1]-yvar[0]
        gd = GeoDict(geodict)
        layers = collections.OrderedDict()
        dictDict = {}
        for key in f.keys():
            keytype = str(type(f[key]))
            if keytype.find('Dataset') > -1:
                if key in REQUIRED_DATASETS:
                    continue
                layers[key] = f[key][:]

        f.close()
        return cls(layers,gd,origin,header,metadata=metadata)

    def setHeader(self,header):
        """
        Set the header dictionary.

        :param header:
          Dictionary with elements:
            - type Type of multi-layer earthquake induced hazard ('shakemap','gfe')
            - version Integer product version (1)
            - process_time Python datetime indicating when data was created.
            - code_version String version of code that created this file (i.e.,'4.0')
            - originator String representing network that created the hazard grid.
            - product_id String representing the ID of the product (may be different from origin ID)
            - map_status String, one of RELEASED, ??
            - event_type String, one of ['ACTUAL','SCENARIO']
        """
        self._header = header.copy() #validate later

    def setOrigin(self,origin):
        """
        Set the origin dictionary.

        Dictionary with elements:
            - id String of event ID (i.e., 'us2015abcd')
            - source String containing originating network ('us')
            - time Float event magnitude
            - lat Float event latitude
            - lon Float event longitude
            - depth Float event depth
            - magnitude Datetime object representing event origin time.
        """
        self._origin = origin.copy() #validate later

    def setMetadata(self,metadata):
        """
        Set the metadata dictionary.

        :param metadata:
          Dictionary of dictionaries of metadata.  Each dictionary can contain any of the following types:
          str,unicode,int,float,long,list,tuple,np.ndarray,dict,datetime.datetime,collections.OrderedDict.
        """
        self._metadata = metadata.copy()

    def getHeader(self):
        """
        Return the header dictionary.

        :returns:
          Header dictionary (see setHeader()).
        """
        return self._header.copy()

    def getOrigin(self):
        """
        Return the origin dictionary.

        :returns:
          Origin dictionary (see setOrigin()).
        """
        return self._origin.copy()

    def getMetadata(self):
        """
        Return the dictionary of arbitrary metadata dictionaries.

        :returns:
          A dictionary of dictionaries containing arbitrary metadata.
        """
        return self._metadata.copy()

if __name__ == '__main__':
    shakefile = sys.argv[1]
    t1 = datetime.datetime.now()
    sgrid = ShakeGrid.load(shakefile)
    t2 = datetime.datetime.now()
    origin = {}
    origin['id'] = sgrid._eventDict['event_id']
    origin['source'] = sgrid._eventDict['event_network']
    origin['time'] = sgrid._eventDict['event_timestamp']
    origin['lat'] = sgrid._eventDict['lat']
    origin['lon'] = sgrid._eventDict['lon']
    origin['depth'] = sgrid._eventDict['depth']
    origin['magnitude'] = sgrid._eventDict['magnitude']

    header = {}
    header['type'] = 'shakemap'
    header['version'] = sgrid._shakeDict['shakemap_version']
    header['process_time'] = sgrid._shakeDict['process_timestamp']
    header['code_version'] = sgrid._shakeDict['code_version']
    header['originator'] = sgrid._shakeDict['shakemap_originator']
    header['product_id'] = sgrid._shakeDict['shakemap_id']
    header['map_status'] = sgrid._shakeDict['map_status']
    header['event_type'] = sgrid._shakeDict['shakemap_event_type']

    layers = collections.OrderedDict()
    for (layername,layerdata) in sgrid.getData().items():
        layers[layername] = layerdata.getData()

    tdict = {'name':'fred','family':{'wife':'wilma','daughter':'pebbles'}}
    mgrid = MultiHazardGrid(layers,sgrid.getGeoDict(),origin,header,metadata={'flintstones':tdict})
    mgrid.save('test.hdf')
    t3 = datetime.datetime.now()
    mgrid2 = MultiHazardGrid.load('test.hdf')
    t4 = datetime.datetime.now()
    xmlmb = os.path.getsize(shakefile)/float(1e6)
    hdfmb = os.path.getsize('test.hdf')/float(1e6)
    xmltime = (t2-t1).seconds + (t2-t1).microseconds/float(1e6)
    hdftime = (t4-t3).seconds + (t4-t3).microseconds/float(1e6)
    print('Input XML file size: %.2f MB (loading time %.3f seconds)' % (xmlmb,xmltime))
    print('Output HDF file size: %.2f MB (loading time %.3f seconds)' % (hdfmb,hdftime))
    os.remove('test.hdf')

max_stars_repo_name: mhearne-usgs/MapIO

max_issues_repo_name: mhearne-usgs/MapIO

max_forks_repo_name: mhearne-usgs/MapIO

reason: import numpy,from scipy


text:
#! /usr/bin/env python3
from popoff.potentials import BuckinghamPotential, buckingham_potentials
from popoff.atom_types import AtomType
from popoff.potential_parameters import BuckinghamParameter
import pytest
import numpy as np

def test_assert_buckingham_potential(buckingham_potential, pot_params):
    assert buckingham_potential.labels == ['Li','O']
    assert buckingham_potential.atype_index == [1,3]
    assert type(buckingham_potential.a) == type(pot_params[0])  #only compares the type, not the content.
    assert type(buckingham_potential.rho) == type(pot_params[1])  #only compares the type, not the content.
    assert type(buckingham_potential.c) == type(pot_params[2])  #only compares the type, not the content.

@pytest.mark.parametrize( 'labels', [(1),(1.0),(True), ([1,'a']), (['Li','O','a']) ])
def test_typeerror_for_labels_in_buckingham_potential(labels, buckingham_potential):
    with pytest.raises(TypeError):
        BuckinghamPotential(labels, buckingham_potential.atype_index, buckingham_potential.a, buckingham_potential.rho, buckingham_potential.c)

@pytest.mark.parametrize( 'atype_index', [(1),(1.0),(True), ([1,'a']), (['Li','O','a']) ])
def test_typeerror_for_atype_index_in_buckingham_potential(atype_index, buckingham_potential):
    with pytest.raises(TypeError):
        BuckinghamPotential(buckingham_potential.labels, atype_index, buckingham_potential.a, buckingham_potential.rho, buckingham_potential.c)

@pytest.mark.parametrize( 'a', [('test'),(1.0),(True),([1,2,3]),
                                (1),(AtomType(1, 'Li', 'Li', 1.555, 1.0))])
def test_typeerror_for_a_in_buckingham_potential(a, buckingham_potential):
    with pytest.raises(TypeError):
        BuckinghamPotential(buckingham_potential.labels, buckingham_potential.atype_index, a, buckingham_potential.rho, buckingham_potential.c)

@pytest.mark.parametrize( 'rho', [('test'),(1.0),(True),([1,2,3]),
                                  (1),(AtomType(1, 'Li', 'Li', 1.555, 1.0))])
def test_typeerror_for_rho_in_buckingham_potential(rho, buckingham_potential):
    with pytest.raises(TypeError):
        BuckinghamPotential(buckingham_potential.labels, buckingham_potential.atype_index, buckingham_potential.a, rho, buckingham_potential.c)

@pytest.mark.parametrize( 'c', [('test'),(1.0),(True),([1,2,3]),
                                (1),(AtomType(1, 'Li', 'Li', 1.555, 1.0))])
def test_typeerror_for_c_in_buckingham_potential(c, buckingham_potential):
    with pytest.raises(TypeError):
        BuckinghamPotential(buckingham_potential.labels, buckingham_potential.atype_index, buckingham_potential.a, buckingham_potential.rho, c)

def test_potential_string_in_buckingham_potential(buckingham_potential):
    assert buckingham_potential.potential_string() == 'pair_coeff 1 3 1.0000 0.1000 0.0000'

def test_assert_potentials_in_buckingham_potentials(potentials, atom_types, pot_params):
    assert potentials.keys() == {'Li-O', 'Ni-O', 'O-O'}
    assert potentials['Li-O'] == [10.00, 0.1, 0.0]
    assert potentials['Ni-O'] == [100.00, 0.1, 0.0]
    assert potentials['O-O'] == [1000.00, 1.0, 1.0]

def test_assert_atom_types_0_in_buckingham_potentials(atom_types):
    assert atom_types[0].atom_type_index == 1
    assert atom_types[0].label == 'Li'
    assert atom_types[0].element_type == 'Li'
    assert atom_types[0].mass == 1.555
    assert atom_types[0].charge == 1.0
    assert atom_types[0].formal_charge == 1.0
    assert atom_types[0].core_shell == None

def test_assert_atom_types_1_in_buckingham_potentials(atom_types):
    assert atom_types[1].atom_type_index == 2
    assert atom_types[1].label == 'Ni'
    assert atom_types[1].element_type == 'Ni'
    assert atom_types[1].mass == 56.111
    assert atom_types[1].charge == 3.0
    assert atom_types[1].formal_charge == 3.0
    assert atom_types[1].core_shell == None

def test_assert_atom_types_2_in_buckingham_potentials(atom_types):
    assert atom_types[2].atom_type_index == 3
    assert atom_types[2].label == 'O core'
    assert atom_types[2].element_type == 'O'
    assert atom_types[2].mass == 14.001
    assert atom_types[2].charge == -0.00
    assert atom_types[2].formal_charge == -0.00
    assert atom_types[2].core_shell == 'core'

def test_assert_atom_types_3_in_buckingham_potentials(atom_types):
    assert atom_types[3].atom_type_index == 3
    assert atom_types[3].label == 'O shell'
    assert atom_types[3].element_type == 'O'
    assert atom_types[3].mass == 1.998
    assert atom_types[3].charge == -2.0
    assert atom_types[3].formal_charge == -2.0
    assert atom_types[3].core_shell == 'shell'

def test_assert_pot_params_0_in_buckingham_potentials(pot_params):
    assert pot_params[0].label_string == 'Li_O_a'
    assert pot_params[0].param_type == 'a'
    assert pot_params[0].value == 1

def test_assert_pot_params_1_in_buckingham_potentials(pot_params):
    assert pot_params[1].label_string == 'Li_O_rho'
    assert pot_params[1].param_type == 'rho'
    assert pot_params[1].value == 0.1

def test_assert_pot_params_2_in_buckingham_potentials(pot_params):
    assert pot_params[2].label_string == 'Li_O_c'
    assert pot_params[2].param_type == 'c'
    assert pot_params[2].value == 0.0

@pytest.mark.parametrize( 'potentials', [(1),(1.0),(True), (['test', 1])])
def test_typeerror_for_potentials_in_buckingham_potentials(potentials, atom_types, pot_params):
    with pytest.raises(TypeError):
        buckingham_potentials(potentials, atom_types, pot_params)

@pytest.mark.parametrize( 'atom_types', [(1),(1.0),(True), (['test',AtomType(1, 'Li', 'Li', 1.555, 1.0)])])
def test_typeerror_for_atom_types_in_buckingham_potentials(potentials, atom_types, pot_params):
    with pytest.raises(TypeError):
        buckingham_potentials(potentials, atom_types, pot_params)

@pytest.mark.parametrize( 'pot_params', [(1),(1.0),(True), ([AtomType(1, 'Li', 'Li', 1.555, 1.0),AtomType(1, 'Li', 'Li', 1.555, 1.0)])])
def test_typeerror_for_pot_params_in_buckingham_potentials(pot_params, potentials, atom_types):
    with pytest.raises(TypeError):
        buckingham_potentials(potentials, atom_types, pot_params)

def test_output_in_buckingham_potentials(atom_types, pot_params):
    output = buckingham_potentials({'Li-O':[10.00, 0.1, 0.0]}, atom_types, pot_params)
    assert type(output) == list
    assert output[0].labels == ['Li', 'O']
    assert output[0].atype_index == [1, 3]
    assert type(output[0].a) == BuckinghamParameter
    assert type(output[0].rho) == BuckinghamParameter
    assert type(output[0].c) == BuckinghamParameter

max_stars_repo_name: pzarabadip/PopOff

max_issues_repo_name: pzarabadip/PopOff

max_forks_repo_name: pzarabadip/PopOff


text:
'''
Train VAE model on data created using an image crop batcher
final model saved into tf_vae/vae.json
'''

import os
os.environ["CUDA_VISIBLE_DEVICES"]="0" # can just override for multi-gpu systems

import tensorflow as tf
import random
import numpy as np
np.set_printoptions(precision=4, edgeitems=6, linewidth=100, suppress=True)
import argparse
import cv2
import matplotlib.pyplot as plt
import time
from vae import ConvVAE, reset_graph

import sys
sys.path.append('../tools')

from rgb_crop_batcher import rgb_crop_batcher

if __name__ == "__main__":
    """
    Main function for executing the .py script.
    Command:
        -p path/<filename>.npy
    """
    # construct the argument parser and parse the arguments
    ap = argparse.ArgumentParser()
    ap.add_argument("-n", "--npy", type=str,
                    help="path to numpy archive")
    args = vars(ap.parse_args())

    crops = np.load(args["npy"])

    # parameters
    z_size=16
    batch_size, crop_size, _, _ = crops.shape

    print("batch_size {}".format(batch_size))
    print("crop_size {}".format(crop_size))

    cv2.imshow("crop",crops[0])
    cv2.waitKey(0)

    reset_graph()
    test_vae = ConvVAE(z_size=z_size,
                       batch_size=batch_size,
                       is_training=False,
                       reuse=False,
                       gpu_mode=True)

    # show reconstruction example
    test_vae.load_json("../../models/0/vae_{}.json".format(180000))
    z = test_vae.encode(crops)
    print(z.shape)

    rec = test_vae.decode(z)
    print(rec.shape)

    np.save("../../output/z_{}.npy".format(batch_size), z)
    np.save("../../output/rec_{}.npy".format(batch_size), rec)

    #for img_idx, img in enumerate(test_batch):
    vis = np.concatenate((crops[0], rec[0]), axis=1)
    cv2.imshow("org vs. rec",cv2.resize(vis,(0,0),fx=4.0,fy=4.0))

    key = cv2.waitKey(0) & 0xFF

max_stars_repo_name: markpp/autoencoders

max_issues_repo_name: markpp/autoencoders

max_forks_repo_name: markpp/autoencoders


text:
import pandas as pd
import numpy as np

class NgramCounter():
    """A class to generate tables of ngram tokens and types from a list of sentences."""

    unk_sign = '<UNK>'
    sent_pad_signs = ['<s>','</s>']

    def __init__(self, sents:[], vocab:[], n:int=3):
        self.sents = sents # Expected to be normalized
        self.vocab = vocab # Can be extracted from another corpus
        self.n = n
        self.widx = [f'w{i}' for i in range(self.n)] # Used for cols and index names

    def generate(self):

        # Convert sentence list to dataframe
        self.S = pd.DataFrame(dict(sent_str=self.sents))

        # Pad sentences
        pad = (self.sent_pad_signs[0] + ' ') *  (self.n - 1)
        self.I = (pad + self.S.sent_str + ' ' + self.sent_pad_signs[1])\
            .str.split(expand=True).stack().to_frame('w0')

        # Set index names
        self.I.index.names = ['sent_num', 'token_num']

        # Remove OOV terms
        self.I.loc[~self.I.w0.isin(self.vocab + self.sent_pad_signs), 'w0'] = self.unk_sign

        # Get sentence lengths (these will include pads)
        self.S['len'] = self.I.groupby('sent_num').w0.count()

        # Add w columns
        for i in range(1, self.n):
            self.I[f'w{i}'] = self.I[f"w{i-1}"].shift(-1)

        # Generate ngrams
        self.NG = []
        for i in range(self.n):
            self.NG.append(self.I.iloc[:, :i+1].copy())

        # Remove spurious rows
        self.NG[i] = self.NG[i].dropna()

        # Generate raw ngram counts and MLEs
        self.LM = []
        for i in range(self.n):
            self.LM.append(self.NG[i].value_counts().to_frame('n'))
            self.LM[i]['mle'] = self.LM[i].n / self.LM[i].n.sum()
            self.LM[i] = self.LM[i].sort_index()

        # Hack to remove single value tuple from unigram table ...
        self.LM[0].index = [i[0] for i in self.LM[0].index]
        self.LM[0].index.name = 'w0'

class NgramLanguageModel():
    """A class to create ngram language models."""

    # Set the Lidstone Smoothing value; LaPlace = 1
    k:float = .5

    def __init__(self, ngc:NgramCounter):
        self.S = ngc.S
        self.LM = ngc.LM
        self.NG = ngc.NG
        self.n = ngc.n
        self.widx = ngc.widx

    def apply_smoothing(self):
        """Applies simple smoothing to ngram type counts to estimate the models."""

        # Z1 and Z2 will hold info about unseen ngrams
        self.Z1 = [None for _ in range(self.n)] # Unseen N grams, but seen N-1 grams
        self.Z2 = [None for _ in range(self.n)] # Unseen N-1 grams too

        # The base vocab size (same as number of unigram types)
        V = len(self.LM[0]) # Inlcides <s> and </s>

        # The number of ngram types
        B = [V**(i+1) for i in range(self.n)]

        # Handle unigram case (no need for smoothing)
        self.LM[0]['p'] = self.LM[0].n / self.LM[0].n.sum()
        self.LM[0]['log_p'] = np.log2(self.LM[0].p)

        # Handle higher order ngrams
        for i in range(1, self.n):

            # MLE
            self.LM[i]['mle2'] = self.LM[i].n / self.LM[i-1].n

            # Emplopy smoothing formula
            self.LM[i]['p'] = (self.LM[i].n + self.k) / (self.LM[i-1].n + B[i-1] * self.k)
            self.LM[i]['log_p'] = np.log2(self.LM[i].p)

            # Unseen N grams, but seen N-1 grams
            self.Z1[i] = np.log2(self.k / (self.LM[i-1].n + B[i-1] * self.k))

            # Unsess N-1 grams too
            self.Z2[i] = np.log2(self.k / B[i-1] * self.k)

            # Tidy up the index
            self.LM[i].sort_index(inplace=True)

    def predict(self, test:NgramCounter):
        """Predicts test sentences with estimated models."""
        self.T = test
        self.PP = []
        p_key = 'log_p'
        for i in range(self.n):
            ng = i + 1
            if i == 0:
                self.T.S[f'ng_{ng}_ll'] = self.T.NG[0]\
                    .join(self.LM[0].log_p, on=self.widx[:ng])\
                    .groupby('sent_num').log_p.sum()
            else:
                self.T.S[f'ng_{ng}_ll'] = self.T.NG[i]\
                    .join(self.LM[i][p_key], on=self.widx[:ng])\
                    .fillna(self.Z1[i]).fillna(self.Z2[i])\
                    .groupby('sent_num')[p_key].sum()

            self.T.S[f'pp{ng}'] = 2**( -self.T.S[f'ng_{ng}_ll'] / self.T.S['len'])

    def generate_text(self, n_sents=20):
        """Generate texts using Shannon Game method."""

        LM = self.LM # For convenience
        i = self.n - 1

        # Start with beginning sentence marker
        words = ['<s>' for _ in range(i)]

        # Sentence counter
        sent_count = 0

        # Generate sentences until we've reached our limit
        while sent_count < n_sents:

            # Get ngram context
            ng = tuple(words[-i:])

            # Get next word
            words.append(LM[i].loc[ng].sample(weights='mle').index.values[0])

            # Terminate when end-of-sentence marker found
            if words[-1] == '</s>':
                sent_count += 1
                if sent_count < n_sents:
                    words.append('<s>')

        # Create text from words
        text = ' '.join(words)

        sents = pd.DataFrame(dict(sent_str=text.split('<s> <s>')))
        sents['len'] = sents.sent_str.str.len()
        sents = sents[sents.len > 0]
        sents.sent_str = sents.sent_str.str.replace('<s> ', '')
        sents.sent_str = sents.sent_str.str.replace(' </s>', '')
        sents.sent_str = sents.sent_str.str.strip()
        sents.sent_str = sents.sent_str.str.replace(r" s ", "'s ", regex=True)

        _ = [print(f"{str(x+1).zfill(2)}. {sent}.\n".upper()) for x, sent in enumerate(sents.sent_str)]

        self.generated_sents = sents


text:
import numpy as np
from pyannote.core import SlidingWindowFeature
import random
import torch

class SpecAugmentor(object):
    """2 way spectrogram augmentation :
        1) Frequency masking
        2) Time masking
    """

    def _spec_augment(mel_spectrogram, frequency_masking_para=27,
                      time_masking_para=100, frequency_mask_num=1, time_mask_num=1, scheduler=None,
                      epoch=None, max_epoch=None):
        """Spectrogram augmentation function from https://arxiv.org/abs/1904.08779
        3 steps processing :
        i) Frequency masking : depending on parameter frequency_masking_para and frequency_mask_num
        ii) Time masking : depending on parameter time_masking_para and time_mask_num

        Parameters
        ----------
        mel_spectrogram :       `SlidingWindowFeature` or (n_samples, n_features ) `numpy.ndarray`
            Features.
        frequency_masking_para: frequency mask parameter F
        time_masking_para:      time mask parameter T
        frequency_mask_num:     number of frequency masks m_F
        time_mask_num:          number of time masks m_T

        Returns
        -------
        masked_melspec : `SlidingWindowFeature` or (n_samples, n_features ) `numpy.ndarray`
            augmented spectrogram
        """
        v = mel_spectrogram.shape[0]
        tau = mel_spectrogram.shape[1]

        if scheduler is not None and max_epoch is not None and epoch is not None:
            frequency_masking_para = int(frequency_masking_para * min(epoch, max_epoch)/max_epoch)
            time_masking_para = int(time_masking_para * min(epoch, max_epoch)/max_epoch)

        time_masking_para = int(min(time_masking_para, 0.2*tau))
        augmented_mel_spectrogram = mel_spectrogram.copy()

        # 1) Frequency masking
        for i in range(frequency_mask_num):
            f = np.random.uniform(low=0.0, high=frequency_masking_para)
            f = int(f)
            f0 = random.randint(0, v - f)
            augmented_mel_spectrogram[f0:f0 + f, :] = 0

        # 2) Time masking
        for i in range(time_mask_num):
            t = np.random.uniform(low=0.0, high=time_masking_para)
            t = int(t)
            t0 = random.randint(0, tau - t)
            augmented_mel_spectrogram[:, t0:t0 + t] = 0

        return augmented_mel_spectrogram

    def __call__(self, features, frequency_masking_para,
                 time_masking_para, nb_frequency_masks, nb_time_masks, scheduler=None,
                 epoch=None, max_epoch=None,
                 sliding_window=None):
        """Apply Google specAugment

        Parameters
        ----------
        features : `SlidingWindowFeature` or (n_samples, n_features ) `numpy.ndarray`
            Features.
        sliding_window : `SlidingWindow`, optional
            Not used.

        Returns
        -------
        normalized : `SlidingWindowFeature` or (n_samples, n_features ) `numpy.ndarray`
            Standardized features
        """

        if isinstance(features, SlidingWindowFeature):
            spec = features.data
        else:
            spec = features

        spec = SpecAugmentor._spec_augment(mel_spectrogram=spec,
                                           frequency_masking_para=frequency_masking_para,
                                           time_masking_para=time_masking_para,
                                           frequency_mask_num=nb_frequency_masks,
                                           time_mask_num=nb_time_masks,
                                           scheduler=scheduler,
                                           epoch=epoch,
                                           max_epoch=max_epoch)

        if isinstance(features, SlidingWindowFeature):
            return SlidingWindowFeature(spec, features.sliding_window)
        else:
            return spec

max_stars_repo_name: nhamilakis/pyannote-audio

max_issues_repo_name: nhamilakis/pyannote-audio

max_forks_repo_name: nhamilakis/pyannote-audio


text:
"""
"""
from __future__ import absolute_import, division, print_function

from unittest import TestCase
from astropy.tests.helper import pytest
import os
import shutil

from astropy.table import Table

try:
    import h5py
    HAS_H5PY = True
except ImportError:
    HAS_H5PY = False

from astropy.config.paths import _find_home

from . import helper_functions
from ..ptcl_table_cache_log_entry import PtclTableCacheLogEntry

# Determine whether the machine is mine
# This will be used to select tests whose
# returned values depend on the configuration
# of my personal cache directory files
aph_home = '/Users/aphearin'
detected_home = _find_home()
if aph_home == detected_home:
    APH_MACHINE = True
else:
    APH_MACHINE = False

__all__ = ('TestPtclTableCacheLogEntry',)

class TestPtclTableCacheLogEntry(TestCase):
    """ Class providing unit testing for `~halotools.sim_manager.PtclTableCacheLogEntry`.
    """
    hard_coded_log_attrs = ['simname', 'version_name', 'redshift', 'fname']

    def setUp(self):
        """ Pre-load various arrays into memory for use by all tests.
        """
        self.dummy_cache_baseloc = helper_functions.dummy_cache_baseloc

        try:
            shutil.rmtree(self.dummy_cache_baseloc)
        except:
            pass
        os.makedirs(self.dummy_cache_baseloc)

        self.simnames = ('bolshoi', 'consuelo',
            'bolshoi', 'bolshoi', 'multidark')
        self.version_names = ('v0', 'v1', 'v2', 'v3', 'v4')
        self.redshifts = (1.2, -0.1, 1.339, 1.3, 100.)

        self.basenames = ('non_existent.hdf5', 'existent.file',
            'existent.hdf5', 'existent.hdf5', 'good.hdf5')
        self.fnames = tuple(os.path.join(self.dummy_cache_baseloc, name)
            for name in self.basenames)

        self.table1 = Table({'x': [1, 2, 3]})

        self.good_table = Table(
            {'ptcl_id': [1, 2, 3],
            'x': [1, 2, 3],
            'y': [1, 2, 3],
            'z': [1, 2, 3],
            'vx': [1, 2, 3],
            'vy': [1, 2, 3],
            'vz': [1, 2, 3],
             })

    def get_scenario_kwargs(self, num_scenario):
        return ({'simname': self.simnames[num_scenario],
            'version_name': self.version_names[num_scenario],
            'redshift': self.redshifts[num_scenario],
            'fname': self.fnames[num_scenario]})

    @pytest.mark.skipif('not HAS_H5PY')
    def test_instantiation(self):
        """ We can instantiate the log entry with a complete set of metadata
        """

        for i in range(len(self.simnames)):
            constructor_kwargs = self.get_scenario_kwargs(i)
            log_entry = PtclTableCacheLogEntry(**constructor_kwargs)
            assert set(log_entry.log_attributes) == set(self.hard_coded_log_attrs)

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario0(self):
        num_scenario = 0
        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))
        assert log_entry.safe_for_cache is False
        assert "The input filename does not exist." in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario1(self):
        num_scenario = 1

        with open(self.fnames[num_scenario], 'w') as f:
            f.write('abc')
        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))
        assert log_entry.safe_for_cache is False
        assert "The input filename does not exist." not in log_entry._cache_safety_message
        assert "The input file must have '.hdf5' extension" in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario2(self):
        num_scenario = 2

        with open(self.fnames[num_scenario], 'w') as f:
            f.write('abc')
        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))
        assert log_entry.safe_for_cache is False
        assert "The input filename does not exist." not in log_entry._cache_safety_message
        assert "The input file must have '.hdf5' extension" not in log_entry._cache_safety_message
        assert "access the hdf5 metadata raised an exception." in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario2a(self):
        num_scenario = 2

        try:
            os.remove(self.fnames[num_scenario])
        except OSError:
            pass
        self.table1.write(self.fnames[num_scenario], path='data')

        f = h5py.File(self.fnames[num_scenario])
        k = list(f.attrs.keys())
        f.close()

        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))
        assert log_entry.safe_for_cache is False
        assert "access the hdf5 metadata raised an exception." not in log_entry._cache_safety_message
        assert "missing the following metadata" in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario2b(self):
        num_scenario = 2

        try:
            os.remove(self.fnames[num_scenario])
        except OSError:
            pass
        self.table1.write(self.fnames[num_scenario], path='data')

        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))

        f = h5py.File(self.fnames[num_scenario])
        for attr in self.hard_coded_log_attrs:
            f.attrs[attr] = getattr(log_entry, attr)
        f.close()

        assert log_entry.safe_for_cache is False
        assert "``particle_mass``" in log_entry._cache_safety_message

        f = h5py.File(self.fnames[num_scenario])
        f.attrs['Lbox'] = 100.
        f.attrs['particle_mass'] = 1.e8
        f.close()
        _ = log_entry.safe_for_cache
        assert "``particle_mass``" not in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario2c(self):
        num_scenario = 2

        try:
            os.remove(self.fnames[num_scenario])
        except OSError:
            pass
        self.table1.write(self.fnames[num_scenario], path='data')

        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))

        f = h5py.File(self.fnames[num_scenario])
        for attr in self.hard_coded_log_attrs:
            if attr != 'redshift':
                f.attrs[attr] = getattr(log_entry, attr)
            else:
                f.attrs[attr] = 0.4
        f.attrs['Lbox'] = 100.
        f.attrs['particle_mass'] = 1.e8
        f.close()

        assert log_entry.safe_for_cache is False
        assert "does not match" in log_entry._cache_safety_message

        f = h5py.File(self.fnames[num_scenario])
        f.attrs['redshift'] = 1.3390001
        f.close()
        assert log_entry.safe_for_cache is False
        assert "does not match" not in log_entry._cache_safety_message

        f = h5py.File(self.fnames[num_scenario])
        f.attrs['redshift'] = '1.3390001'
        f.close()
        assert log_entry.safe_for_cache is False
        assert "does not match" not in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario3(self):
        num_scenario = 3

        try:
            os.remove(self.fnames[num_scenario])
        except:
            pass
        self.table1.write(self.fnames[num_scenario], path='data')

        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))

        f = h5py.File(self.fnames[num_scenario])
        for attr in self.hard_coded_log_attrs:
            if attr != 'redshift':
                f.attrs[attr] = getattr(log_entry, attr)
            else:
                f.attrs[attr] = 0.4
        f.attrs['Lbox'] = 100.
        f.attrs['particle_mass'] = 1.e8
        f.close()

        assert log_entry.safe_for_cache is False
        substr = "must at a minimum have the following columns"
        assert substr in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_scenario4(self):
        num_scenario = 4

        try:
            os.remove(self.fnames[num_scenario])
        except:
            pass
        self.good_table.write(self.fnames[num_scenario], path='data')

        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))

        f = h5py.File(self.fnames[num_scenario])
        for attr in self.hard_coded_log_attrs:
            f.attrs[attr] = getattr(log_entry, attr)
        f.attrs['Lbox'] = 2.
        f.attrs['particle_mass'] = 1.e8
        f.close()

        assert log_entry.safe_for_cache is False
        substr = "must be bounded by [0, Lbox]."
        assert substr in log_entry._cache_safety_message

    @pytest.mark.skipif('not HAS_H5PY')
    def test_passing_scenario(self):
        num_scenario = 4

        try:
            os.remove(self.fnames[num_scenario])
        except:
            pass
        self.good_table.write(self.fnames[num_scenario], path='data')

        log_entry = PtclTableCacheLogEntry(**self.get_scenario_kwargs(num_scenario))

        f = h5py.File(self.fnames[num_scenario])
        for attr in self.hard_coded_log_attrs:
            f.attrs[attr] = getattr(log_entry, attr)
        f.attrs['Lbox'] = 100.
        f.attrs['particle_mass'] = 1.e8
        f.close()

        assert log_entry.safe_for_cache is True, log_entry._cache_safety_message
        substr = "The particle catalog is safe to add to the cache log."
        assert substr in log_entry._cache_safety_message

    def tearDown(self):
        try:
            shutil.rmtree(self.dummy_cache_baseloc)
        except:
            pass

reason: from astropy


text:
import numpy as np
import twopoint
import fitsio as fio
import healpy as hp
from numpy.lib.recfunctions import append_fields, rename_fields
from .stage import PipelineStage, NOFZ_NAMES
import subprocess
import os
import warnings

class nofz(PipelineStage):
    name = "nofz"
    inputs = {}
    outputs = {
        "weight"        : "weight.npy"          ,
        "nz_source"     : "nz_source_zbin.npy"  ,
        "nz_lens"       : "nz_lens_zbin.npy"    ,
        "nz_source_txt" : "source.nz"           ,
        "nz_lens_txt"   : "lens.nz"             ,
        "gold_idx"      : "gold_idx.npy"        ,
        "lens_idx"      : "lens_idx.npy"        ,
        "ran_idx"       : "ran_idx.npy"         ,
        "randoms"       : "randoms_zbin.npy"    ,
        "cov_ini"       : "cov.ini"             ,
        "2pt"           : "2pt.fits"            ,
        "metadata"      : "metadata.yaml"
    }

    def __init__(self, param_file):
        """
        Produces n(z)s from input catalogs.
        """
        super(nofz,self).__init__(param_file)

        # Load data
        self.load_data()
        if 'pzbin_col' in self.gold.dtype.names:
            print 'ignoring any specified bins, since a bin column has been supplied in gold file'

        # Construct new weight and cache - move to catalog
        if 'weight' in self.pz.dtype.names:
            self.weight = np.sqrt(self.pz['weight'] * self.shape['weight'])
        else:
            self.weight = self.shape['weight']
        filename = self.output_path("weight")
        np.save(filename, np.vstack((self.gold['objid'], self.weight)).T)
        # deal with photo-z weights for lenses later...
        # Setup binning
        if self.params['pdf_type']!='pdf':
            self.z       = (np.linspace(0.,4.,401)[1:]+np.linspace(0.,4.,401)[:-1])/2.+1e-4 # 1e-4 for buzzard redshift files
            self.dz      = self.z[1]-self.z[0]
            self.binlow  = self.z-self.dz/2.
            self.binhigh = self.z+self.dz/2.
        else:
            self.binlow  = np.array(self.params['pdf_z'][:-1])
            self.binhigh = np.array(self.params['pdf_z'][1:])
            self.z       = (self.binlow+self.binhigh)/2.
            self.dz      = self.binlow[1]-self.binlow[0]

        if hasattr(self.params['zbins'], "__len__"):
            self.tomobins = len(self.params['zbins']) - 1
            self.binedges = self.params['zbins']
        else:
            self.tomobins = self.params['zbins']
            self.binedges = self.find_bin_edges(self.pz['pzbin'][self.mask], self.tomobins, w = self.shape['weight'][self.mask])

        if self.params['lensfile'] != 'None':
            if hasattr(self.params['lens_zbins'], "__len__"):
                self.lens_tomobins = len(self.params['lens_zbins']) - 1
                self.lens_binedges = self.params['lens_zbins']
            else:
                self.lens_tomobins = self.params['lens_zbins']
                self.lens_binedges = self.find_bin_edges(self.lens_pz['pzbin'], self.lens_tomobins, w = self.lens['weight'])

        return

    def run(self):

        # Calculate source n(z)s and write to file
        if self.params['has_sheared']:
            pzbin = [self.pz['pzbin'],self.pz_1p['pzbin'],self.pz_1m['pzbin'],self.pz_2p['pzbin'],self.pz_2m['pzbin']]
        else:
            pzbin = self.pz['pzbin']

        if self.params['pdf_type']!='pdf':
            zbin, self.nofz = self.build_nofz_bins(
                               self.tomobins,
                               self.binedges,
                               pzbin,
                               self.pz_nofz['pzstack'],
                               self.params['pdf_type'],
                               self.weight,
                               shape=True)
        else:
            pdfs = np.zeros((len(self.pz),len(self.z)))
            for i in range(len(self.z)):
                pdfs[:,i] = self.pz['pzstack'+str(i)]
            zbin, self.nofz = self.build_nofz_bins(
                               self.tomobins,
                               self.binedges,
                               pzbin,
                               pdfs,
                               self.params['pdf_type'],
                               self.weight,
                               shape=True)

        self.get_sigma_e(zbin,self.tomobins,self.shape)
        self.get_neff(zbin,self.tomobins,self.shape)

        if self.params['has_sheared']:
            zbin,zbin_1p,zbin_1m,zbin_2p,zbin_2m=zbin
            np.save(self.output_path("nz_source"), np.vstack((self.gold['objid'], zbin,zbin_1p,zbin_1m,zbin_2p,zbin_2m)).T)
        else:
            np.save(self.output_path("nz_source"), np.vstack((self.gold['objid'], zbin)).T)

        # Calculate lens n(z)s and write to file
        if self.params['lensfile'] != 'None':
            lens_zbin, self.lens_nofz = self.build_nofz_bins(
                                         self.lens_tomobins,
                                         self.lens_binedges,
                                         self.lens_pz['pzbin'],
                                         self.lens_pz['pzstack'],
                                         self.params['lens_pdf_type'],
                                         self.lens['weight'])
            np.save(self.output_path("nz_lens")  , np.vstack((self.lens['objid'], lens_zbin)).T)

            ran_binning = np.digitize(self.randoms['ranbincol'], self.lens_binedges, right=True) - 1
            np.save(self.output_path("randoms"), ran_binning)

            self.get_lens_neff(lens_zbin,self.lens_tomobins,self.lens)

    def write(self):
        """
        Write lens and source n(z)s to fits file for tomographic and non-tomographic cases.
        """

        nz_source = twopoint.NumberDensity(
                     NOFZ_NAMES[0],
                     self.binlow,
                     self.z,
                     self.binhigh,
                     [self.nofz[i,:] for i in range(self.tomobins)])

        nz_source.ngal      = self.neff
        nz_source.sigma_e   = self.sigma_e
        nz_source.area      = self.area
        kernels             = [nz_source]
        np.savetxt(self.output_path("nz_source_txt"), np.vstack((self.binlow, self.nofz)).T)

        if self.params['lensfile'] != 'None':
            nz_lens      = twopoint.NumberDensity(
                            NOFZ_NAMES[1],
                            self.binlow,
                            self.z,
                            self.binhigh,
                            [self.lens_nofz[i,:] for i in range(self.lens_tomobins)])
            nz_lens.ngal = self.lens_neff
            nz_lens.area = self.area
            kernels.append(nz_lens)
            np.savetxt(self.output_path("nz_lens_txt"), np.vstack((self.binlow, self.lens_nofz)).T)

        data             = twopoint.TwoPointFile([], kernels, None, None)
        data.to_fits(self.output_path("2pt"), clobber=True)

        self.write_metadata()

    def write_metadata(self):
        import yaml
        data = {
            "neff": self.neff,
            "neffc": self.neffc,
            "tomobins": self.tomobins,
            "sigma_e": self.sigma_e,
            "sigma_ec": self.sigma_ec,
            "mean_e1":self.mean_e1,
            "mean_e2":self.mean_e2,
            "area": self.area,
            "repository_version:": find_git_hash(),
        }
        if 'pzbin_col' in self.gold.dtype.names:
            data["source_bins"] = "gold_file_bins"
        else:
            if type(self.binedges) is list:
                data.update({ "source_bins" : self.binedges })
            else:
                data.update({ "source_bins" : self.binedges.tolist() })

        if self.params['lensfile'] != 'None':
            data.update({ "lens_neff" : self.lens_neff,
                          "lens_tomobins" : self.lens_tomobins,
                          "lens_bins" : self.lens_binedges })
        print data
        filename = self.output_path('metadata')
        open(filename, 'w').write(yaml.dump(data))

    def load_array(self,d,file):

        if self.params['has_sheared'] & (file=='shapefile'):
            d['flags_1p'] = 'flags_select_1p'
            d['flags_1m'] = 'flags_select_1m'
            d['flags_2p'] = 'flags_select_2p'
            d['flags_2m'] = 'flags_select_2m'

        if self.params['pdf_type']=='pdf':
            keys = [key for key in d.keys() if (d[key] is not None)&(key is not 'pzstack')]
        else:
            keys = [key for key in d.keys() if (d[key] is not None)]

        if 'objid' in keys:
            dtypes = [('objid','i8')]
        else:
            raise ValueError('missing object id in '+file)
        dtypes += [(key,'f8') for key in keys if (key is not 'objid')]
        if self.params['pdf_type']=='pdf':
            dtypes  += [('pzstack_'+str(i),'f8') for i in range(len(self.params['pdf_z']))]

        fits = fio.FITS(self.params[file])[-1]
        array = fits.read(columns=[d[key] for key in keys])

        array = rename_fields(array,{v: k for k, v in d.iteritems()})

        if ('weight' not in array.dtype.names) & (file=='shapefile'):
            array = append_fields(array, 'weight', np.ones(len(array)), usemask=False)

        if self.params['pdf_type']=='pdf':
            for i in range(len(self.params['pdf_z'])):
                array['pzstack'+str(i)]=fits.read(columns=d['pzstack']+str(i))

        if np.any(np.diff(array['objid']) < 1):
            raise ValueError('misordered or duplicate ids in '+file)

        return array

    def load_data(self):
        """
        Load data files.
        """

        import time
        import importlib
        col = importlib.import_module('.'+self.params['dict_file'],'pipeline')

        # def lowercase_array(arr):
        #     old_names = arr.dtype.names
        #     new_names = [name.lower() for name in old_names]
        #     renames   = dict(zip(old_names, new_names))
        #     return rename_fields(arr, renames)
        t0 = time.time()

        self.gold      = self.load_array(col.gold_dict, 'goldfile')
        print 'Done goldfile',time.time()-t0,self.gold.dtype.names
        self.shape     = self.load_array(col.shape_dict, 'shapefile')
        print 'Done shapefile',time.time()-t0,self.shape.dtype.names
        self.pz        = self.load_array(col.pz_bin_dict, 'photozfile')
        print 'Done pzfile',time.time()-t0,self.pz.dtype.names
        if self.params['has_sheared']:
            self.pz_1p = self.load_array(col.pz_bin_dict, 'photozfile_1p')
            print 'Done pz1pfile',time.time()-t0,self.pz_1p.dtype.names
            self.pz_1m = self.load_array(col.pz_bin_dict, 'photozfile_1m')
            print 'Done pz1mfile',time.time()-t0,self.pz_1m.dtype.names
            self.pz_2p = self.load_array(col.pz_bin_dict, 'photozfile_2p')
            print 'Done pz2pfile',time.time()-t0,self.pz_2p.dtype.names
            self.pz_2m = self.load_array(col.pz_bin_dict, 'photozfile_2m')
            print 'Done pz2mfile',time.time()-t0,self.pz_2m.dtype.names
        self.pz_nofz   = self.load_array(col.pz_stack_dict, 'photozfile_nz')
        print 'Done pznofzfile',time.time()-t0,self.pz_nofz.dtype.names
        if self.params['lensfile'] != 'None':
            self.lens      = self.load_array(col.lens_dict, 'lensfile')
            print 'Done lensfile',time.time()-t0,self.lens.dtype.names
            self.lens_pz   = self.load_array(col.lens_pz_dict, 'lensfile')
            print 'Done lens_pzfile',time.time()-t0,self.lens_pz.dtype.names

        if 'm1' not in self.shape.dtype.names:
            self.shape = append_fields(self.shape, 'm1', self.shape['m2'], usemask=False)
        if 'm2' not in self.shape.dtype.names:
            self.shape = append_fields(self.shape, 'm2', self.shape['m1'], usemask=False)
        if self.params['oneplusm']==False:
            print 'converting m to 1+m'
            self.shape['m1'] = np.copy(self.shape['m1'])+1.
            self.shape['m2'] = np.copy(self.shape['m2'])+1.
        if 'c1' in self.shape.dtype.names:
            self.shape['e1'] -= self.shape['c1']
            self.shape['e2'] -= self.shape['c2']
            self.shape['c1'] = None
            self.shape['c2'] = None
        if self.params['flip_e2']==True:
            print 'flipping e2'
            self.shape['e2']*=-1
        if 'pzbin' not in self.lens_pz.dtype.names:
            self.lens_pz = append_fields(self.lens_pz, 'pzbin', self.lens_pz['pzstack'], usemask=False)
        if 'pzstack' not in self.lens_pz.dtype.names:
            self.lens_pz = append_fields(self.lens_pz, 'pzstack', self.lens_pz['pzbin'], usemask=False)

        if not ((len(self.gold)==len(self.shape))
            & (len(self.gold)==len(self.pz))
            & (len(self.gold)==len(self.pz_nofz))):
            raise ValueError('shape, gold, or photoz length mismatch')
        if self.params['has_sheared']:
            if not ((len(self.gold)==len(self.pz_1p))
                & (len(self.gold)==len(self.pz_1m))
                & (len(self.gold)==len(self.pz_2p))
                & (len(self.gold)==len(self.pz_2m))):
                raise ValueError('shape, gold, or photoz length mismatch')
        if self.params['lensfile'] != 'None':
            if (len(self.lens)!=len(self.lens_pz)):
                raise ValueError('lens and lens_pz length mismatch')

        if self.params['lensfile'] != 'None':
            keys = [key for key in col.ran_dict.keys() if (col.ran_dict[key] is not None)]
            fits = fio.FITS(self.params['randomfile'])[-1]

            dtypes=[(key,'f8') for key in keys]
            self.randoms = np.empty(fits.read_header()['NAXIS2'], dtype = dtypes)
            for key in keys:
                self.randoms[key]=fits.read(columns=[col.ran_dict[key]])

        if self.params['test_run']==True:
            idx = np.random.choice(np.arange(len(self.gold)),100000,replace=False)
            np.save(self.output_path("gold_idx"), idx)
            self.gold    = self.gold[idx]
            self.shape   = self.shape[idx]
            self.pz      = self.pz[idx]
            self.pz_nofz = self.pz_nofz[idx]
            if self.params['has_sheared']:
                self.pz_1p   = self.pz_1p[idx]
                self.pz_1m   = self.pz_1m[idx]
                self.pz_2p   = self.pz_2p[idx]
                self.pz_2m   = self.pz_2m[idx]
            if self.params['lensfile'] != 'None':
                idx = np.random.choice(np.arange(len(self.lens)),100000,replace=False)
                np.save(self.output_path("lens_idx"), idx)
                self.lens    = self.lens[idx]
                self.lens_pz = self.lens_pz[idx]
                idx = np.random.choice(np.arange(len(self.randoms)),100000,replace=False)
                np.save(self.output_path("ran_idx"), idx)
                self.randoms = self.randoms[idx]

        if 'pzbin_col' in self.gold.dtype.names:
            mask = (self.gold['pzbin_col'] >= 0)
        else:
            mask = (self.pz['pzbin'] > self.params['zlims'][0]) & (self.pz['pzbin'] <= self.params['zlims'][1])
            if self.params['has_sheared']:
                mask_1p = (self.pz_1p['pzbin'] > self.params['zlims'][0]) & (self.pz_1p['pzbin'] <= self.params['zlims'][1])
                mask_1m = (self.pz_1m['pzbin'] > self.params['zlims'][0]) & (self.pz_1m['pzbin'] <= self.params['zlims'][1])
                mask_2p = (self.pz_2p['pzbin'] > self.params['zlims'][0]) & (self.pz_2p['pzbin'] <= self.params['zlims'][1])
                mask_2m = (self.pz_2m['pzbin'] > self.params['zlims'][0]) & (self.pz_2m['pzbin'] <= self.params['zlims'][1])

        print 'ngal',np.sum(mask),np.sum(mask_1p),np.sum(mask_1m),np.sum(mask_2p),np.sum(mask_2m)
        if 'flags' in self.shape.dtype.names:
            mask = mask & (self.shape['flags']==0)
            if self.params['has_sheared']:
                mask_1p = mask_1p & (self.shape['flags_1p']==0)
                mask_1m = mask_1m & (self.shape['flags_1m']==0)
                mask_2p = mask_2p & (self.shape['flags_2p']==0)
                mask_2m = mask_2m & (self.shape['flags_2m']==0)
        print 'ngal',np.sum(mask),np.sum(mask_1p),np.sum(mask_1m),np.sum(mask_2p),np.sum(mask_2m)
        if 'flags' in self.pz.dtype.names:
            mask = mask & (self.pz['flags']==0)
            if self.params['has_sheared']:
                mask_1p = mask_1p & (self.pz_1p['flags']==0)
                mask_1m = mask_1m & (self.pz_1m['flags']==0)
                mask_2p = mask_2p & (self.pz_2p['flags']==0)
                mask_2m = mask_2m & (self.pz_2m['flags']==0)

        print 'hardcoded spt region cut'
        mask = mask & (self.shape['dec']<-35)
        if self.params['has_sheared']:
            mask_1p = mask_1p & (self.shape['dec']<-35)
            mask_1m = mask_1m & (self.shape['dec']<-35)
            mask_2p = mask_2p & (self.shape['dec']<-35)
            mask_2m = mask_2m & (self.shape['dec']<-35)

        np.save('radec.npy',np.vstack((self.shape['ra'],self.shape['dec'])).T[mask])

        print np.sum(mask)
        if 'footprintfile' in self.params.keys():
            print 'cutting catalog to footprintfile'
            footmask = np.in1d(hp.ang2pix(4096, np.pi/2.-np.radians(self.shape['dec']),np.radians(self.shape['ra']), nest=False),
                               fio.FITS(self.params['footprintfile'])[-1].read()['HPIX'],assume_unique=False)
            mask = mask & footmask
            if self.params['has_sheared']:
                mask_1p = mask_1p & footmask
                mask_1m = mask_1m & footmask
                mask_2p = mask_2p & footmask
                mask_2m = mask_2m & footmask

        print 'ngal final',np.sum(mask),np.sum(mask & mask_1p & mask_1m & mask_2p & mask_2m)

        if self.params['has_sheared']:
            full_mask     = mask | mask_1p | mask_1m | mask_2p | mask_2m
            self.pz_1p    = self.pz_1p[full_mask]
            self.pz_1m    = self.pz_1m[full_mask]
            self.pz_2p    = self.pz_2p[full_mask]
            self.pz_2m    = self.pz_2m[full_mask]
            self.mask     = mask[full_mask]
            self.mask_1p  = mask_1p[full_mask]
            self.mask_1m  = mask_1m[full_mask]
            self.mask_2p  = mask_2p[full_mask]
            self.mask_2m  = mask_2m[full_mask]
        else:
            full_mask  = mask
            self.mask  = mask[full_mask]
        self.gold      = self.gold[full_mask]
        self.shape     = self.shape[full_mask]
        self.pz        = self.pz[full_mask]
        self.pz_nofz   = self.pz_nofz[full_mask]

        return

    def build_nofz_bins(self, zbins, edge, bin_col, stack_col, pdf_type, weight,shape=False):
        """
        Build an n(z), non-tomographic [:,0] and tomographic [:,1:].
        """

        if shape&(self.params['has_sheared']):
            if 'pzbin_col' in self.gold.dtype.names:
                xbins = self.gold['pzbin_col']
            else:
                xbins0=[]
                for x in bin_col:
                    xbins0.append(np.digitize(x, edge, right=True) - 1)
                xbins = xbins0[0]
        else:
            if 'pzbin_col' in self.gold.dtype.names:
                xbins0 = self.gold['pzbin_col']
            else:
                xbins0 = np.digitize(bin_col, edge, right=True) - 1
            xbins=xbins0

        # Stack n(z)
        nofz  = np.zeros((zbins, len(self.z)))

        # MC Sample of pdf or redmagic (if redmagic, takes random draw from gaussian of width 0.01)
        if (pdf_type == 'sample') | (pdf_type == 'rm'):
            if pdf_type == 'rm':
                stack_col = np.random.normal(stack_col, self.lens_pz['pzerr']*np.ones(len(stack_col)))
            for i in range(zbins):
                mask        =  (xbins == i)
                if shape:
                    mask = mask&self.mask
                    if self.params['has_sheared']:
                        mask_1p = (xbins0[1] == i)&self.mask_1p
                        mask_1m = (xbins0[2] == i)&self.mask_1m
                        mask_2p = (xbins0[3] == i)&self.mask_2p
                        mask_2m = (xbins0[4] == i)&self.mask_2m
                        m1   = np.mean(self.shape['m1'][mask&self.mask])
                        m2   = np.mean(self.shape['m2'][mask&self.mask])
                        m1   += (np.mean(self.shape['e1'][mask_1p&self.mask_1p]) - np.mean(self.shape['e1'][mask_1m&self.mask_1m])) / (2.*self.params['dg'])
                        m2   += (np.mean(self.shape['e2'][mask_2p&self.mask_2p]) - np.mean(self.shape['e2'][mask_2m&self.mask_2m])) / (2.*self.params['dg'])
                        m1   = m1*np.ones(len(mask))
                        m2   = m2*np.ones(len(mask))
                    else:
                        m1 = self.shape['m1']
                        m2 = self.shape['m2']
                    weight *= (m1+m2)/2.
                nofz[i,:],b =  np.histogram(stack_col[mask], bins=np.append(self.binlow, self.binhigh[-1]), weights=weight[mask])
                nofz[i,:]   /= np.sum(nofz[i,:]) * self.dz

        # Stacking pdfs
        elif pdf_type == 'pdf':
            for i in xrange(zbins):
                mask      =  (xbins == i)
                if shape:
                    mask = mask&self.mask
                nofz[i,:] =  np.sum((stack_col[mask].T * weight[mask]).T, axis=0)
                nofz[i,:] /= np.sum(nofz[i,:]) * self.dz

        return xbins0, nofz

    def find_bin_edges(self, x, nbins, w=None):
        """
        For an array x, returns the boundaries of nbins equal (possibly weighted by w) bins.
        From github.com/matroxel/destest.
        """

        if w is None:
            xs = np.sort(x)
            r  = np.linspace(0., 1., nbins + 1.) * (len(x) - 1)
            return xs[r.astype(int)]

        fail = False
        ww   = np.sum(w) / nbins
        i    = np.argsort(x)
        k    = np.linspace(0.,1., nbins + 1.) * (len(x) - 1)
        k    = k.astype(int)
        r    = np.zeros((nbins + 1))
        ist  = 0
        for j in xrange(1,nbins):
            if k[j]  < r[j-1]:
                print 'Random weight approx. failed - attempting brute force approach'
                fail = True
                break

            w0 = np.sum(w[i[ist:k[j]]])
            if w0 <= ww:
                for l in xrange(k[j], len(x)):
                    w0 += w[i[l]]
                    if w0 > ww:
                        r[j] = x[i[l]]
                        ist  = l
                        break
            else:
                for l in xrange(k[j], 0, -1):
                    w0 -= w[i[l]]
                    if w0 < ww:
                        r[j] = x[i[l]]
                        ist  = l
                        break

        if fail:
            ist = np.zeros((nbins+1))
            ist[0]=0
            for j in xrange(1, nbins):
                wsum = 0.
                for k in xrange(ist[j-1].astype(int), len(x)):
                    wsum += w[i[k]]
                    if wsum > ww:
                        r[j]   = x[i[k-1]]
                        ist[j] = k
                        break

        r[0]  = x[i[0]]
        r[-1] = x[i[-1]]

        return r

    def get_neff(self, zbin, tomobins, cat):
        """
        Calculate neff for catalog.
        """

        if not hasattr(self,'area'):
            self.get_area()

        self.neff = []
        self.neffc = []
        for i in range(self.tomobins):
            if self.params['has_sheared']:
                mask = (zbin[0] == i)
                mask_1p = (zbin[1] == i)
                mask_1m = (zbin[2] == i)
                mask_2p = (zbin[3] == i)
                mask_2m = (zbin[4] == i)

                m1   = np.mean(cat['m1'][mask&self.mask])
                m2   = np.mean(cat['m2'][mask&self.mask])
                m1   += (np.mean(cat['e1'][mask_1p&self.mask_1p]) - np.mean(cat['e1'][mask_1m&self.mask_1m])) / (2.*self.params['dg'])
                m2   += (np.mean(cat['e2'][mask_2p&self.mask_2p]) - np.mean(cat['e2'][mask_2m&self.mask_2m])) / (2.*self.params['dg'])
                m1   = m1*np.ones(len(mask))
                m2   = m2*np.ones(len(mask))
            else:
                mask = (zbin == i)
                m1 = cat['m1']
                m2 = cat['m2']

            mask = mask & self.mask
            print i,np.sum(mask)

            if self.params['has_sheared']:
                e1  = cat['e1'][mask]
                e2  = cat['e2'][mask]
                w   = cat['weight'][mask]
                s   = (m1[mask]+m2[mask])/2.
                var = cat['cov00'][mask] + cat['cov11'][mask]
                var[var>2] = 2.
            else:
                e1  = cat['e1'][mask]
                e2  = cat['e2'][mask]
                w   = cat['weight'][mask]
                s = (m1[mask]+m2[mask])/2.
                snvar = 0.24#np.sqrt((cat['e1'][mask][cat['snr'][mask]>100].var()+cat['e2'][mask][cat['snr'][mask]>100].var())/2.)
                var = 1./w - snvar**2
                var[var < 0.] = 0.
                w[w > snvar**-2] = snvar**-2

            a    = np.sum(w)**2
            b    = np.sum(w**2)
            c    = self.area * 60. * 60.

            self.neff.append( a/b/c )
            self.neffc.append( ((self.sigma_ec[i]**2 * np.sum(w * s)**2) / np.sum(w**2 * (s**2 * self.sigma_ec[i]**2 + var/2.))) / self.area / 60**2 )

        return

    def get_lens_neff(self, zbin, tomobins, cat):
        """
        Calculate neff for catalog.
        """

        if not hasattr(self,'area'):
            self.get_area()

        self.lens_neff = []
        for i in range(tomobins):
          mask = (zbin == i)
          a    = np.sum(cat['weight'][mask])**2
          b    = np.sum(cat['weight'][mask]**2)
          c    = self.area * 60. * 60.

          self.lens_neff.append(np.asscalar( a/b/c ))

        return

    def get_sigma_e(self, zbin, tomobins, cat):
        """
        Calculate sigma_e for shape catalog.
        """

        if not hasattr(self,'area'):
            self.get_area()

        self.sigma_e = []
        self.sigma_ec = []
        self.mean_e1 = []
        self.mean_e2 = []
        for i in range(tomobins):
            if self.params['has_sheared']:
                mask = (zbin[0] == i)
                mask_1p = (zbin[1] == i)
                mask_1m = (zbin[2] == i)
                mask_2p = (zbin[3] == i)
                mask_2m = (zbin[4] == i)

                m1   = np.mean(cat['m1'][mask&self.mask])
                m2   = np.mean(cat['m2'][mask&self.mask])
                print 'rg',m1,m2
                m1   += (np.mean(cat['e1'][mask_1p&self.mask_1p]) - np.mean(cat['e1'][mask_1m&self.mask_1m])) / (2.*self.params['dg'])
                m2   += (np.mean(cat['e2'][mask_2p&self.mask_2p]) - np.mean(cat['e2'][mask_2m&self.mask_2m])) / (2.*self.params['dg'])
                print 'rs',(np.mean(cat['e1'][mask_1p&self.mask_1p]) - np.mean(cat['e1'][mask_1m&self.mask_1m])) / (2.*self.params['dg']),(np.mean(cat['e2'][mask_2p&self.mask_2p]) - np.mean(cat['e2'][mask_2m&self.mask_2m])) / (2.*self.params['dg'])
                m1   = m1*np.ones(len(mask))
                m2   = m2*np.ones(len(mask))
                print 'sn',i,np.sum(self.mask_1p),np.sum(self.mask_1m),np.sum(self.mask_2p),np.sum(self.mask_2m)
                print 'sn',i,np.sum(mask_1p&self.mask_1p),np.sum(mask_1m&self.mask_1m),np.sum(mask_2p&self.mask_2p),np.sum(mask_2m&self.mask_2m)
            else:
                mask = (zbin == i)
                m1 = cat['m1']
                m2 = cat['m2']

            mask = mask & self.mask

            if self.params['has_sheared']:
                e1  = cat['e1'][mask]
                e2  = cat['e2'][mask]
                w   = np.ones(len(cat[mask]))#cat['weight'][mask]
                s   = (m1[mask]+m2[mask])/2.
                var = cat['cov00'][mask] + cat['cov11'][mask]
                var[var>2] = 2.
            else:
                e1  = cat['e1'][mask]
                e2  = cat['e2'][mask]
                w   = cat['weight'][mask]
                s = (m1[mask]+m2[mask])/2.
                snvar = 0.24#np.sqrt((cat['e1'][mask][cat['snr'][mask]>100].var()+cat['e2'][mask][cat['snr'][mask]>100].var())/2.)
                print i,'snvar',snvar
                var = 1./w - snvar**2
                var[var < 0.] = 0.
                w[w > snvar**-2] = snvar**-2
                print 'var',var.min(),var.max()

            self.mean_e1.append(np.asscalar(np.average(cat['e1'][mask],weights=cat['weight'][mask])))
            self.mean_e2.append(np.asscalar(np.average(cat['e2'][mask],weights=cat['weight'][mask])))
            a1 = np.sum(w**2 * (e1-self.mean_e1[i])**2)
            a2 = np.sum(w**2 * (e2-self.mean_e2[i])**2)
            b  = np.sum(w**2)
            c  = np.sum(w * s)
            print len(w),w
            d  = np.sum(w)
            print i,a1,a2,b,c,d

            self.sigma_e.append( np.sqrt( (a1/c**2 + a2/c**2) * (d**2/b) / 2. ) )
            self.sigma_ec.append( np.sqrt( np.sum(w**2 * (e1**2 + e2**2 - var)) / (2.*np.sum(w**2 * s**2)) ) )

        return

    def get_area(self):

        if hasattr(self,'area'):
            return

        if self.params['area']=='None':

            import healpy as hp

            pix=hp.ang2pix(4096, np.pi/2.-np.radians(self.shape['dec']),np.radians(self.shape['ra']), nest=True)
            area=hp.nside2pixarea(4096)*(180./np.pi)**2
            mask=np.bincount(pix)>0
            self.area=np.sum(mask)*area
            self.area=float(self.area)
            print self.area

        else:

            self.area = self.params['area']

        return

def find_git_hash():
    try:
        dirname = os.path.dirname(os.path.abspath(__file__))
        head = subprocess.check_output("cd {0}; git show-ref HEADS".format(dirname), shell=True)
    except subprocess.CalledProcessError:
        head = "UNKNOWN"
        warnings.warn("Unable to find git repository commit ID in {}".format(dirname))
    return head.split()[0]

reason: import numpy,from numpy


text:
import pygad as pg
import pygad.plotting
import matplotlib.pyplot as plt
import numpy as np
import glob
from multiprocessing import Pool
from scipy import stats
import utils

filename = __file__

def cycle_r_max(args):
    halo = args[0]
    definition = args[1]

    path = '/ptmp/mpa/naab/REFINED/%s/SF_X/4x-2phase/out/snap_%s_4x_???' % (halo, halo)
    max = int(sorted(glob.glob(path))[-1][-3:])
    s, h, g = pg.prepare_zoom('/ptmp/mpa/naab/REFINED/%s/SF_X/4x-2phase/out/snap_%s_4x_%s' % (halo, halo, max), gas_trace='/u/mihac/data/%s/4x-2phase/gastrace_%s' % (halo, definition), star_form=None)

    cycle_r_max = s.gas['cycle_r_max'] #[s.gas['num_recycled'] > -1]
    particle_mass = s.gas['mass_at_ejection'] #[s.gas['num_recycled'] > -1]
    mask = np.max(s.gas['cycle_r_max'], axis=1) > '100 kpc'
    cycle_r_max_lim = cycle_r_max[mask]
    ID_mask = pg.IDMask(s.gas['ID'][mask])

    mass_in_bins, edges2, count2 = stats.binned_statistic(cycle_r_max.flatten(),
                                     particle_mass.flatten(), statistic='sum', bins=np.logspace(-1.1, 4, 90))
    mass_in_bins /= 1e10

    fig, ax = plt.subplots(4, figsize=(16,12))
    ax[0].set_xlabel("$r_{max}$ [kpc]")
    ax[0].set_ylabel("mass [$10^{10}\ M_{\odot}$]")
    ax[0].set_xscale('log')
    ax[0].set_xlim((1e-2, 1e4))
    ax[0].set_yscale('log')
    ax[0].grid(True)
    ax[0].bar(edges2[:-1], mass_in_bins, width=np.diff(edges2))

    ax[1].set_xlabel("$r_{max}$ [kpc]")
    ax[1].set_ylabel("count")
    ax[1].grid(True)
    ax[1].set_yscale('log')
    ax[1].hist(cycle_r_max.flatten(), bins=np.linspace(0, 300, 301))

    ax[2].set_title("same as above but logscale")
    ax[2].set_xlabel("$r_{max}$ [kpc]")
    ax[2].set_ylabel("count")
    ax[2].grid(True)
    ax[2].set_yscale('log')
    ax[2].set_xscale('log')
    ax[2].hist(cycle_r_max.flatten(), bins=np.logspace(-1.1, 4, 90), log=True)

    ax[3].set_title("same as above but just particles with max(cycle_r_max) > 100)")
    ax[3].set_xlabel("$r_{max}$ [kpc]")
    ax[3].set_ylabel("count")
    ax[3].grid(True)
    ax[3].set_yscale('log')
    ax[3].set_xscale('log')
    ax[3].hist(cycle_r_max_lim.flatten(), bins=np.logspace(-1.1, 4, 90), log=True)

    fig.tight_layout()

    plt.savefig(filename.split("/")[-1][:-3] + '_' + halo + '_' + definition + ".png", bbox_inches='tight')
    plt.clf()
    plt.close()
    #for snap in sorted(glob.glob(path))[70::20]:
    #    s, h, g = pg.prepare_zoom(snap, gas_trace=None, star_form=None)
    #    fig, ax = plt.subplots(1, 3)
    #    fig.suptitle(s.cosmic_time())
    #    pltargs = dict(clim=(10**3.5, 10**6.5))
    #    pg.plotting.image(s.gas, extent='400 kpc', ax=ax[0], **pltargs)
    #    pg.plotting.image(s.gas[ID_mask], extent='400 kpc', ax=ax[1], **pltargs)
    #    pg.plotting.image(s.gas[~ID_mask], extent='400 kpc', ax=ax[2], **pltargs)
    #    plt.tight_layout()
    #    print '%s/%s/%s_%s.png' % (halo, definition, definition, snap.split('/')[-1])
    #    plt.savefig('%s/%s/%s_%s.png' % (halo, definition, definition, snap.split('/')[-1]))
    #    plt.clf()
    #    plt.close()

p = Pool(4)
p.map(plot, utils.combinations)


text:
#################################################
### THIS FILE WAS AUTOGENERATED! DO NOT EDIT! ###
#################################################
# file to edit: dev_nb/Sentence_Weight.ipynb
import requests
from lexrank import LexRank
from lexrank.mappings.stopwords import STOPWORDS
import spacy
from spacy.tokens import Doc, Span
import numpy as np

urls= ['http://www.gutenberg.org/cache/epub/29444/pg29444.txt',
       'http://www.gutenberg.org/cache/epub/31428/pg31428.txt',
       'http://www.gutenberg.org/cache/epub/4908/pg4908.txt',
       'http://www.gutenberg.org/cache/epub/33504/pg33504.txt',
       'http://www.gutenberg.org/cache/epub/39017/pg39017.txt',
       'http://www.gutenberg.org/cache/epub/37157/pg37157.txt',
       'http://www.gutenberg.org/cache/epub/15207/pg15207.txt',
       'http://www.gutenberg.org/cache/epub/10773/pg10773.txt',
       'http://www.gutenberg.org/cache/epub/31624/pg31624.txt',
       'http://www.gutenberg.org/cache/epub/5192/pg5192.txt',
       'http://www.gutenberg.org/cache/epub/40030/pg40030.txt',
       'http://www.gutenberg.org/cache/epub/50992/pg50992.txt',
       'http://www.gutenberg.org/cache/epub/50880/pg50880.txt',
       'http://www.gutenberg.org/cache/epub/36525/pg36525.txt',
       'http://www.gutenberg.org/cache/epub/47167/pg47167.txt',
       'http://www.gutenberg.org/cache/epub/47464/pg47464.txt',
       'http://www.gutenberg.org/cache/epub/33397/pg33397.txt',
       'http://www.gutenberg.org/cache/epub/29782/pg29782.txt',
       'http://www.gutenberg.org/cache/epub/32857/pg32857.txt',
       'http://www.gutenberg.org/cache/epub/26262/pg26262.txt']
docs = [requests.get(u).text for u in urls]

lxr = LexRank(docs, stopwords=STOPWORDS['en'])
nlp =  spacy.load('en_core_web_md')

def tfidf(span):return sum([tok.prob for tok in span])*1.0/len(span)
Span.set_extension("weight", getter=tfidf, force=True)

def ranktfidf(doc): return np.argsort(np.array([sent._.weight for sent in doc.sents]))[::-1]
Doc.set_extension("ranktfidf", getter=ranktfidf, force=True)

max_stars_repo_name: visualbuffer/mcq

max_issues_repo_name: visualbuffer/mcq

max_forks_repo_name: visualbuffer/mcq


text:
from functools import namedtuple
import random
import numpy as np

from game2048.expectimax import board_to_move
from keras.models import Sequential, Model, Input
from keras.layers import Flatten, Dense, Dropout, Activation, concatenate
from keras.layers import Conv2D, MaxPooling2D, BatchNormalization
from keras.optimizers import RMSprop

from env_edit_reward import state2tensor, Env2048

Guide = namedtuple('Guides', ('state', 'action', 'next_state', 'reward'))
model_path = '/home/srtp/2048-api-master/zit_dqn.h5'

class Guides:
    def __init__(self, capacity):
        self.capacity = capacity
        self.memory = []
        self.position = 0

    def push(self, *args):
        if len(self.memory) < self.capacity:
            self.memory.append(None)
        self.memory[self.position] = Guide(*args)
        self.position = (self.position + 1) % self.capacity

    def sample(self, batch_size):
        return random.sample(self.memory, batch_size)

    def ready(self, batch_size):
        return len(self.memory) >= batch_size

    def __len__(self):
        return len(self.memory)

class DQN(object):
    def __init__(
            self,
            learning_rate=0.01,
            reward_decay=0.9,
            e_greedy=0.9,
            replace_target_iter=50,
            memory_size=2000,
            batch_size=128,
            e_greedy_increment=None,
            output_graph=False,
    ):

        self.lr = learning_rate
        self.gamma = reward_decay
        self.epsilon_max = e_greedy
        # 动作会有随机， 1-e_greedy的概率随机决策
        self.replace_target_iter = replace_target_iter

        self.batch_size = batch_size
        # 随着训练降低epsilon
        self.epsilon_increment = e_greedy_increment
        self.epsilon = 0 if e_greedy_increment is not None else self.epsilon_max
        self.learn_step_counter = 1
        self.memory = Guides(memory_size)
        self._build_net()

    def _build_net(self):
        # 构建网络
        # 待修改
        # ==================evualte_network=========
        filters = 128
        eval_inputs = Input((4, 4, 16))
        conv1 = Conv2D(filters, (1, 2), padding='same', kernel_initializer='he_uniform')(eval_inputs)
        conv2 = Conv2D(filters, (2, 1), padding='same', kernel_initializer='he_uniform')(eval_inputs)
        conv3 = Conv2D(filters, (1, 3), padding='same', kernel_initializer='he_uniform')(eval_inputs)
        conv4 = Conv2D(filters, (3, 1), padding='same', kernel_initializer='he_uniform')(eval_inputs)
        conv5 = Conv2D(filters, (4, 4), padding='same', kernel_initializer='he_uniform')(eval_inputs)

        hidden = concatenate([Flatten()(x) for x in [conv1, conv2, conv3, conv4, conv5]])
        eval_hidden = BatchNormalization()(hidden)
        eval_hidden = Activation('relu')(eval_hidden)

        for width in [512, 128]:
            eval_hidden = Dense(width, kernel_initializer='he_uniform')(eval_hidden)
            eval_hidden = BatchNormalization()(eval_hidden)
            eval_hidden = Activation('relu')(eval_hidden)

        self.q_eval = Dense(4)(eval_hidden)

        # ==================target network=============

        target_inputs = Input(shape=(4, 4, 16))
        conv21 = Conv2D(filters, (1, 2), padding='same', kernel_initializer='he_uniform')(target_inputs)
        conv22 = Conv2D(filters, (2, 1), padding='same', kernel_initializer='he_uniform')(target_inputs)
        conv23 = Conv2D(filters, (1, 3), padding='same', kernel_initializer='he_uniform')(target_inputs)
        conv24 = Conv2D(filters, (3, 1), padding='same', kernel_initializer='he_uniform')(target_inputs)
        conv25 = Conv2D(filters, (4, 4), padding='same', kernel_initializer='he_uniform')(target_inputs)

        hidden2 = concatenate([Flatten()(x) for x in [conv21, conv22, conv23, conv24, conv25]])
        target_hidden = BatchNormalization()(hidden2)
        target_hidden = Activation('relu')(target_hidden)

        for width in [512, 128]:
            target_hidden = Dense(width, kernel_initializer='he_uniform')(target_hidden)
            target_hidden = BatchNormalization()(target_hidden)
            target_hidden = Activation('relu')(target_hidden)

        self.q_next = Dense(4)(target_hidden)

        self.model1 = Model(target_inputs, self.q_next)
        self.model2 = Model(eval_inputs, self.q_eval)
        rmsprop = RMSprop(lr=self.lr)
        self.model1.compile(loss='mean_squared_error', optimizer=rmsprop, metrics=['accuracy'])
        self.model2.compile(loss='mean_squared_error', optimizer=rmsprop, metrics=['accuracy'])

        self.model1.summary()
        self.model2.summary()

    def target_replace_op(self):
        model_params = self.model2.get_weights()
        self.model1.set_weights(model_params)
        print('+++++++=======================+++++++++++++++')
        print("params has changed")
        print('+++++++=======================+++++++++++++++')

    def choose_action(self, state, greedy=False):
        # 决策
        if greedy or np.random.random() < self.epsilon:
            board = np.expand_dims(state, axis=0)
            actions_value = self.model1.predict(board)
            action = np.argmax(actions_value)

            print('=====================================')
            print('action_value', actions_value)
            print('-------------------------------------')
            print(action)
            print('=====================================')

        else:
            # action = np.random.randint(0, 4)
            gameboard = np.zeros((4, 4))
            for i in range(4):
                for j in range(4):
                    num = np.argmax(state[i][j])
                    if num == 0:
                        gameboard[i][j] = 0
                    else:
                        gameboard[i][j] = 2 ** num

            action = board_to_move(gameboard)

        return action

    def learn(self):
        # 模型训练
        if self.memory.ready(self.batch_size):
            guides = self.memory.sample(self.batch_size)

            learn_state = []
            learn_q_target = []

            for guide in guides:
                q_eval = self.model2.predict(np.expand_dims(guide.state, axis=0))
                if guide.next_state is None:
                    q_next = 0
                else:
                    q_next = self.model1.predict(np.expand_dims(guide.next_state, axis=0))

                q_target = q_eval.copy()

                """
                print(q_next)
                print(q_target)
                """
                learn_state.append(guide.state)
                reward_now = guide.reward
                action = guide.action

                q_target[0, action] = reward_now + self.gamma * np.max(q_next)
                learn_q_target.append(q_target.reshape(4))

            self.model2.fit(np.array(learn_state), np.array(learn_q_target), epochs=10)
            print('traning')
            self.model2.save(model_path)

            self.epsilon = self.epsilon + self.epsilon_increment if self.epsilon < self.epsilon_max else self.epsilon_max
            self.learn_step_counter += 1

    def training(self, num_episodes):
        last50 = []
        env = Env2048(score_to_win=1024)

        for i_episode in range(num_episodes):
            env.reset()
            state = env.to_tensor()
            for i_step in range(100000):
                action = self.choose_action(state)
                _, reward, done, _ = env.step(action)

                if not done:
                    next_state = env.to_tensor()
                else:
                    next_state = None

                self.memory.push(state, action, next_state, reward)

                # Move to the next state
                state = next_state

                if done:
                    print(env.get_state())
                    break

                # print(self.memory.memory[i_step])

            self.learn()

            if i_episode % self.replace_target_iter == 0:
                self.target_replace_op()

            last50.append(env.get_return())
            if len(last50) > 50:
                last50.pop(0)
            self.average_ = np.average(last50)

            print('Episode %3d, steps %4d, score %5d, avg %5.0f, eps %.4f' % (
            self.learn_step_counter, i_step, env.get_return(), self.average_, self.epsilon))

            self.learn()

if __name__ == '__main__':
    n_episodes = 2000
    player = DQN(e_greedy_increment=0.001)

    player.training(n_episodes)


text:
import numpy as np
import torch
import torchvision
import torchvision.transforms as transforms
from torch.utils.data.sampler import Sampler

def get_mnist(batch_size, train, path, augmentation=False, std=0.0, shuffle=True, adversarial=False, subset=1000,
              workers=0):
    classes = np.arange(0, 10)

    if augmentation:
        transform = transforms.Compose([
            transforms.ToTensor(),
            transforms.Lambda(lambda x: x + torch.Tensor(x.size()).normal_(mean=0.0, std=std)),  # add gaussian noise
            transforms.Lambda(lambda x: torch.clamp(x, 0.0, 1.0))
        ])
    else:
        transform = transforms.Compose([
            transforms.ToTensor(),
            transforms.Lambda(lambda x: torch.clamp(x, 0.0, 1.0))
        ])

    dataset = torchvision.datasets.MNIST(root=path,
                                         train=train,
                                         download=True,
                                         transform=transform)

    if adversarial:
        np.random.seed(123)  # load always the same random subset
        indices = np.random.choice(np.arange(dataset.__len__()), subset)

        subset_sampler = SubsetSampler(indices)

        dataloader = torch.utils.data.DataLoader(dataset,
                                                 batch_size=batch_size,
                                                 shuffle=False,
                                                 sampler=subset_sampler,
                                                 num_workers=workers
                                                 )

    else:
        dataloader = torch.utils.data.DataLoader(dataset,
                                                 batch_size=batch_size,
                                                 shuffle=shuffle,
                                                 num_workers=workers
                                                 )
    return dataloader, dataset, classes

def get_cifar10(batch_size, train, path, augmentation=False, std=0.0, shuffle=True, adversarial=False, subset=1000,
                workers=0):
    classes = ['plane', 'car', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship', 'truck']

    if augmentation:
        transform = transforms.Compose([
            transforms.ToTensor(),
            transforms.Lambda(lambda x: x + torch.Tensor(x.size()).normal_(0.0, std)),  # add gaussian noise
            transforms.Lambda(lambda x: torch.clamp(x, 0.0, 1.0))
        ])
    else:
        transform = transforms.Compose([
            transforms.ToTensor(),
            transforms.Lambda(lambda x: torch.clamp(x, 0.0, 1.0))
        ])

    dataset = torchvision.datasets.CIFAR10(root=path,
                                           train=train,
                                           download=True,
                                           transform=transform)

    if adversarial:
        np.random.seed(123)  # load always the same random subset
        indices = np.random.choice(np.arange(dataset.__len__()), subset)

        subset_sampler = SubsetSampler(indices)
        dataloader = torch.utils.data.DataLoader(dataset,
                                                 batch_size=batch_size,
                                                 shuffle=False,
                                                 sampler=subset_sampler,
                                                 num_workers=workers
                                                 )
    else:
        dataloader = torch.utils.data.DataLoader(dataset,
                                                 batch_size=batch_size,
                                                 shuffle=shuffle,
                                                 num_workers=workers
                                                 )
    return dataloader, dataset, classes

class SubsetSampler(Sampler):
    def __init__(self, indices):
        self.indices = indices
        self.shuffle = False

    def __iter__(self):
        if self.shuffle:
            np.random.shuffle(self.indices)
        return (self.indices[i] for i in range(len(self.indices)))

    def __len__(self):
        return len(self.indices)

    def set_shuffle(self, shuffle):
        self.shuffle = shuffle