text:
# -*- coding: utf-8 -*-
# -*- coding: utf-8 -*-
"""
Authors: Tim Hessels
         UNESCO-IHE 2016
Contact: t.hessels@unesco-ihe.org
Repository: https://github.com/wateraccounting/wa
Module: Collect/MOD16
"""

# import general python modules
import os
import numpy as np
import pandas as pd
import gdal
import urllib
import urllib2
import re
import glob
from joblib import Parallel, delayed
from bs4 import BeautifulSoup
import datetime
import math

# Water Accounting modules
import wa
import wa.General.raster_conversions as RC
import wa.General.data_conversions as DC

def DownloadData(Dir, Startdate, Enddate, latlim, lonlim, timestep, Waitbar, cores, hdf_library, remove_hdf):
    """
    This function downloads MOD13 16-daily data

    Keyword arguments:
    Dir -- 'C:/file/to/path/'
    Startdate -- 'yyyy-mm-dd'
    Enddate -- 'yyyy-mm-dd'
    latlim -- [ymin, ymax] (values must be between -90 and 90)
    lonlim -- [xmin, xmax] (values must be between -180 and 180)
    cores -- The number of cores used to run the routine. It can be 'False'
             to avoid using parallel computing routines.
    Waitbar -- 1 (Default) will print a waitbar
    """

    # Check start and end date and otherwise set the date
    if not Startdate:
        Startdate = pd.Timestamp('2000-01-01')
    if not Enddate:
        Enddate = pd.Timestamp('2014-12-31')

    # Make an array of the days of which the ET is taken
    if timestep == 'monthly':
        Dates = pd.date_range(Startdate,Enddate,freq = 'M')
        TIMESTEP = 'Monthly'
    elif timestep == '8-daily':
        Dates = Make_TimeStamps(Startdate,Enddate)
        TIMESTEP = '8_Daily'

    # Create Waitbar
    if Waitbar == 1:
        import wa.Functions.Start.WaitbarConsole as WaitbarConsole
        total_amount = len(Dates)
        amount = 0
        WaitbarConsole.printWaitBar(amount, total_amount, prefix = 'Progress:', suffix = 'Complete', length = 50)

    # Make directory for the MODIS ET data
    output_folder=os.path.join(Dir,'Evaporation','MOD16', TIMESTEP)
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    TilesVertical, TilesHorizontal = wa.Collect.MOD15.DataAccess.Get_tiles_from_txt(output_folder, hdf_library, latlim, lonlim)

    # Pass variables to parallel function and run
    args = [output_folder, TilesVertical, TilesHorizontal,latlim, lonlim, timestep, hdf_library]
    if not cores:
        for Date in Dates:
            RetrieveData(Date, args)
            if Waitbar == 1:
                amount += 1
                WaitbarConsole.printWaitBar(amount, total_amount, prefix = 'Progress:', suffix = 'Complete', length = 50)
        results = True
    else:
        results = Parallel(n_jobs=cores)(delayed(RetrieveData)(Date, args)
                                         for Date in Dates)

    if remove_hdf == 1:
        # Remove all .hdf files
        os.chdir(output_folder)
        files = glob.glob("*.hdf")
        for f in files:
            os.remove(os.path.join(output_folder, f))

        # Remove all .txt files
        files = glob.glob("*.txt")
        for f in files:
            os.remove(os.path.join(output_folder, f))

	return results

def RetrieveData(Date, args):
    """
    This function retrieves MOD16 ET data for a given date from the
    ftp://ftp.ntsg.umt.edu/ server.

    Keyword arguments:
    Date -- 'yyyy-mm-dd'
    args -- A list of parameters defined in the DownloadData function.
    """
    # Argument
    [output_folder, TilesVertical, TilesHorizontal,latlim, lonlim, timestep, hdf_library] = args

    # Collect the data from the MODIS webpage and returns the data and lat and long in meters of those tiles
    try:
        Collect_data(TilesHorizontal,TilesVertical,Date,output_folder, timestep, hdf_library)
    except:
        print "Was not able to download the file"

    # Define the output name of the collect data function
    name_collect = os.path.join(output_folder, 'Merged.tif')

    # Reproject the MODIS product to epsg_to
    epsg_to ='4326'
    name_reprojected = RC.reproject_MODIS(name_collect, epsg_to)

    # Clip the data to the users extend
    data, geo = RC.clip_data(name_reprojected, latlim, lonlim)

    if timestep == 'monthly':
         ETfileName = os.path.join(output_folder, 'ET_MOD16A2_mm-month-1_monthly_'+Date.strftime('%Y')+'.' + Date.strftime('%m')+'.01.tif')
    elif timestep == '8-daily':
         ETfileName = os.path.join(output_folder, 'ET_MOD16A2_mm-8days-1_8-daily_'+Date.strftime('%Y') + '.' + Date.strftime('%m') + '.' + Date.strftime('%d') + '.tif')

    DC.Save_as_tiff(name=ETfileName, data=data, geo=geo, projection='WGS84')

    # remove the side products
    os.remove(os.path.join(output_folder, name_collect))
    os.remove(os.path.join(output_folder, name_reprojected))

    return()

def Collect_data(TilesHorizontal,TilesVertical,Date,output_folder,timestep, hdf_library):
    '''
    This function downloads all the needed MODIS tiles from ftp.ntsg.umt.edu/pub/MODIS/NTSG_Products/MOD16/MOD16A2_MONTHLY.MERRA_GMAO_1kmALB/ as a hdf file.

    Keywords arguments:
    TilesHorizontal -- [TileMin,TileMax] max and min horizontal tile number
    TilesVertical -- [TileMin,TileMax] max and min vertical tile number
    Date -- 'yyyy-mm-dd'
    output_folder -- 'C:/file/to/path/'
    '''

    # Make a new tile for the data
    sizeX=int(TilesHorizontal[1]-TilesHorizontal[0]+1)*1200
    sizeY=int(TilesVertical[1]-TilesVertical[0]+1)*1200
    DataTot=np.ones((sizeY,sizeX))* -9999

    # Make a new tile for the lat and long info
    LatMet=np.zeros((sizeY))
    LongMet=np.zeros((sizeX))

    # Create the Lat and Long of the MODIS tile in meters
    for Vertical in range(int(TilesVertical[0]),int(TilesVertical[1])+1):
        Distance=926.625 # resolution of a MODIS pixel in meter
        countY=int((TilesVertical[1]-TilesVertical[0]+1)-(Vertical-TilesVertical[0]))
        LatMet[int((countY-1)*1200):int((countY)*1200)]=np.linspace(((8-Vertical)*1200+0.5)*Distance,((8-Vertical)*1200+1199.5)*Distance,1200)

        for Horizontal in range(int(TilesHorizontal[0]),int(TilesHorizontal[1])+1):
            countX=int(Horizontal-TilesHorizontal[0]+1)
            LongMet[int((countX-1)*1200):int((countX)*1200)]=np.linspace(((Horizontal-18)*1200+0.5)*Distance,((Horizontal-18)*1200+1199.5)*Distance,1200)

            # Set the download to zero again
            downloaded = 0

            if hdf_library is not None:
                os.chdir(hdf_library)
                hdf_name = glob.glob("MOD16A2.A%s%03s.h%02dv%02d.*" %(Date.strftime('%Y'), Date.strftime('%j'), Horizontal, Vertical))

                if len(hdf_name) == 1:
                    hdf_file = os.path.join(hdf_library, hdf_name[0])

                    if os.path.exists(hdf_file):
                        downloaded = 1
                        data = Open_mod16_data(hdf_file)
                        countYdata=(TilesVertical[1]-TilesVertical[0]+2)-countY
                        DataTot[(int(countYdata)-1)*1200:int(countYdata)*1200,(int(countX)-1)*1200:int(countX)*1200]=data*0.1

            while downloaded == 0:

                # Download the MODIS FPAR data
                if timestep == 'monthly':
                    url = 'http://files.ntsg.umt.edu/data/NTSG_Products/MOD16/MOD16A2_MONTHLY.MERRA_GMAO_1kmALB/Y%s/M%02s/' %(Date.strftime('%Y'), Date.strftime('%m'))

                if timestep == '8-daily':
                    url = 'http://files.ntsg.umt.edu/data/NTSG_Products/MOD16/MOD16A2.105_MERRAGMAO/Y%s/D%03s/' %(Date.strftime('%Y'), Date.strftime('%j'))

                # Get files on FTP server
                f = urllib2.urlopen(url)

                # Sum all the files on the server
                soup = BeautifulSoup(f, "lxml")

                try:
                    for i in soup.findAll('a', attrs = {'href': re.compile('(?i)(hdf)$')}):

                        # Find the file with the wanted tile number
                        nameHDF=str(i)
                        HDF_name = nameHDF.split('>')[-2][:-3]
                        Hfile=HDF_name[18:20]
                        Vfile=HDF_name[21:23]
                        if int(Vfile) is int(Vertical) and int(Hfile) is int(Horizontal):

                            HTTP_name = url + HDF_name
                            output_name = os.path.join(output_folder, HDF_name)
                            if not os.path.isfile(output_name):

                                while downloaded == 0:

                                    urllib.urlretrieve(HTTP_name,output_name)

                                    statinfo = os.stat(output_name)
                                    # Say that download was succesfull
                                    if int(statinfo.st_size) > 1000:
                                       downloaded = 1

                            # Open .hdf only band with ET and collect all tiles to one array
                            data = Open_mod16_data(output_name)
                            countYdata=(TilesVertical[1]-TilesVertical[0]+2)-countY
                            DataTot[(int(countYdata)-1)*1200:int(countYdata)*1200,(int(countX)-1)*1200:int(countX)*1200]=data*0.1
                            DataTot[DataTot>3000]=-9999
                            downloaded = 1
                            del data

                except:
                    proj='PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
                    data=np.ones((1200, 1200)) * (-9999)
                    countYdata=(TilesVertical[1] - TilesVertical[0] + 2) - countY
                    DataTot[(countYdata - 1) * 1200:countYdata * 1200,(countX - 1) * 1200:countX * 1200] = data * 0.1
                    downloaded = 1

	 # Make geotiff file
    name2 = os.path.join(output_folder, 'Merged.tif')
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(name2, DataTot.shape[1], DataTot.shape[0], 1, gdal.GDT_Float32, ['COMPRESS=LZW'])
    try:
        dst_ds.SetProjection(proj)
    except:
        proj='PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
        x1 = (TilesHorizontal[0] - 18) * 1200 * Distance
        x4 = (TilesVertical[0] - 9) * 1200 * -1 * Distance
        geo = [x1, Distance, 0.0, x4, 0.0, -Distance]
        geo_t = tuple(geo)
        dst_ds.SetProjection(proj)

    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo_t)
    dst_ds.GetRasterBand(1).WriteArray(DataTot)
    dst_ds = None

    return(DataTot,LatMet,LongMet)

def Open_mod16_data(output_name):

    dataset=gdal.Open(output_name)
    sdsdict=dataset.GetMetadata('SUBDATASETS')
    sdslist =[sdsdict[k] for k in sdsdict.keys() if '_1_NAME' in k]
    sds=[]

    sds=gdal.Open(sdslist[0])
    data = np.float_(sds.ReadAsArray())
    sds = None
    return(data)

def Make_TimeStamps(Startdate,Enddate):
    '''
    This function determines all time steps of which the FPAR must be downloaded
    The time stamps are 8 daily.

    Keywords arguments:
    Startdate -- 'yyyy-mm-dd'
    Enddate -- 'yyyy-mm-dd'
    '''

    # Define the DOY and year of the start day
    DOY = datetime.datetime.strptime(Startdate,'%Y-%m-%d').timetuple().tm_yday
    Year = datetime.datetime.strptime(Startdate,'%Y-%m-%d').timetuple().tm_year

    # Define the year of the end day
    YearEnd = datetime.datetime.strptime(Enddate,'%Y-%m-%d').timetuple().tm_year

    # Change the DOY of the start day into a DOY of MODIS day (8-daily) and create new startdate
    DOYstart = int(math.floor(DOY / 8.0) * 8) + 1
    DOYstart = str('%s-%s' %(DOYstart, Year))
    Day = datetime.datetime.strptime(DOYstart, '%j-%Y')
    Month = '%02d' % Day.month
    Day = '%02d' % Day.day
    Startdate = (str(Year) + '-' + str(Month) + '-' + str(Day))

    # Create the start and end data for the whole year
    YearStartDate = pd.date_range(Startdate, Enddate, freq = 'AS')
    YearEndDate = pd.date_range(Startdate, Enddate, freq = 'A')

    # Define the amount of years that are involved
    AmountOfYear = YearEnd - Year

    # If the startday is not in the same year as the enddate
    if AmountOfYear > 0:
        for i in range(0, AmountOfYear+1):
            if i is 0:
                Startdate1 = Startdate
                Enddate1 = YearEndDate[0]
                Dates = pd.date_range(Startdate1, Enddate1, freq = '8D')
            if i is AmountOfYear:
                Startdate1 = YearStartDate[-1]
                Enddate1 = Enddate
                Dates1 = pd.date_range(Startdate1, Enddate1, freq = '8D')
                Dates = Dates.union(Dates1)
            if (i is not AmountOfYear and i is not 0):
                Startdate1 = YearStartDate[i-AmountOfYear-1]
                Enddate1 = YearEndDate[i]
                Dates1 = pd.date_range(Startdate1, Enddate1, freq = '8D')
                Dates = Dates.union(Dates1)

    # If the startday is in the same year as the enddate
    if AmountOfYear is 0:
        Dates = pd.date_range(Startdate, Enddate, freq = '8D')

    return(Dates)

max_forks_repo_name: wateraccounting/wa


text:
# ******************************************************************************
"""
Checks if the two passed-in faces are the same person.

This module implements the face recognition based on the technique
'One-Shot Learning for Face Recognition'. This technique is useful for
recognizing a person from a small reference dataset.

Private Functions:
    . _parse                    parses the script arguments,
    . _preprocess_image         returns an image array from an image path,
    . _findCosineSimilarity     returns the cosine similarity between two images,
    . _findEuclideanDistance    eeturns the euclidean distance between two images,

Public Class:
    .  FaceMatcher              a class to recognize a face against a referencial,

Public Methods:
    . compare                   compares the target face from the source face,

@namespace      -
@author         <author_name>
@since          0.0.0
@version        0.0.0
@licence        MIT. Copyright (c) 2020 Mobilabs <contact@mobilabs.fr>
"""
# ******************************************************************************
import argparse

import keras
from keras.applications.imagenet_utils import preprocess_input
import numpy as np
from models.vggface import VGGFace

WEIGHTS = './weights/vgg_face_weights.h5'
# Cosine similarity threshold:
EPSILON = 0.40

# -- Private Functions ---------------------------------------------------------

def _parse():
    """Parses the script arguments.

    ### Parameters:
        param1 ():          none.

    ### Returns:
        (str):              returns the path of the source.
        (str):              returns the path of the target.

    ### Raises:
        none
    """
    parser = argparse.ArgumentParser()
    parser.add_argument('-s', '--source', type=str, help='the path of the image to recognize')
    parser.add_argument('-t', '--target', type=str, help='the path of the reference image')
    args = parser.parse_args()
    return args.source, args.target

def _preprocess_image(image_path):
    """Returns an image array from an image path.

    ### Parameters:
        param1 (str):       the image path.

    ### Returns:
        (array):            returns an image array.

    ### Raises:
        none
    """
    img = keras.preprocessing.image.load_img(image_path, target_size=(224, 224))
    img = keras.preprocessing.image.img_to_array(img)
    img = np.expand_dims(img, axis=0)
    img = preprocess_input(img)
    return img

def _findCosineSimilarity(vector_source, vector_target):
    """Returns the cosine similarity between the source and the target.

    Parameters:
        param1 (array):     the vectorized source face image.
        param2 (array):     the vectorized target face image.

    Returns:
        (num):              returns the cosine similarity.

    Raises:
        -
    """
    a = np.matmul(np.transpose(vector_source), vector_target)
    b = np.sum(np.multiply(vector_source, vector_source))
    c = np.sum(np.multiply(vector_target, vector_target))
    return 1 - (a / (np.sqrt(b) * np.sqrt(c)))

def _findEuclideanDistance(vector_source, vector_target):
    """Returns the euclidean distance between the source and the target.

    Parameters:
        param1 (array):     the vectorized source face image.
        param2 (array):     the vectorized target face image.

    Returns:
        (num):              returns the euclidian distance.

    Raises:
        -

    """
    euclidean_distance = vector_source - vector_target
    euclidean_distance = np.sum(np.multiply(euclidean_distance, euclidean_distance))
    euclidean_distance = np.sqrt(euclidean_distance)
    return euclidean_distance

# -- Public --------------------------------------------------------------------

class FaceMatcher:
    """A class to recognize a face against a referencial.

    ### Attributes:
        net (obj):          the VGGFace neural network.

    ### Methods:
        compare(source, target):
            Compares the target face from the source face.

    ### Raises:
        none
    """

    def __init__(self):
        """Creates and initializes the VGGFace neural network."""
        self.net = VGGFace(include_top=False)
        self.net.load_weights(WEIGHTS)

    def compare(self, source, target):
        """Compares the target face from the source face.

        ### Parameters:
            param1 (str):   the image path of the source image.
            param2 (str):   the image path of the target image.

        ### Returns:
            (num):          the euclidean distance,
            (num):          the cosine similarity,
            (bool):         true if the two face images show the same person,

        ### Raises:
            -
        """
        vector_source = self.net.predict(_preprocess_image(source))[0, :]
        vector_target = self.net.predict(_preprocess_image(target))[0, :]

        cosine_similarity = _findCosineSimilarity(vector_source, vector_target)
        euclidean_distance = _findEuclideanDistance(vector_source, vector_target)

        if cosine_similarity < EPSILON:
            match = True
        else:
            match = False

        return euclidean_distance, cosine_similarity, match

if __name__ == '__main__':
    matcher = FaceMatcher()
    source, target = _parse()
    distance, similarity, match = matcher.compare(source, target)
    if match is True:
        m = 'Yes!'
    else:
        m = 'No!'

    print('euclidean distance: ' + str(distance)
          + ', cosine similarity: ' + str(similarity)
          + ', match: ' + m)

# -- o ---

ext: py

reason: import numpy


text:
from .cost import Cost

from .layers.particle import Particle
from .layers.particle import ParticleInput
from .regularization.particle_regularize_l2 import ParticleRegularizeL2
from .regularization.particle_regularize_l2plus import ParticleRegularizeL2Plus
from .regularization.particle_regularize_orthogonal import ParticleRegularizeOrthogonal

import numpy as np
import json

class ParticleConvNetwork(object):

    def __init__(self, particle_input=None, cost="mse", regularizer=None):
        self.layers = []
        self.cost_name = cost
        self.cost_function = Cost.get(cost)
        self.cost_d_function = Cost.get_d(cost)
        self.lock_built = False
        self.regularizer = regularizer
        self.particle_input = particle_input

    def append(self, layer):
        """
        Appends a layer to the network

        :param layer:
        :return:
        """
        self.layers.append(layer)

    def build(self):
        """
        Handle networks layer dimensions checks, other possible initializations

        Release build lock

        :return:
        """
        self.lock_built = True

    def predict_single(self, data_X):
        """
        Same as predict, but only one sample
        """
        return self.predict(data_X.reshape((1, len(data_X))))

    def predict(self, data_X):
        """
        Pass given input through network to compute the output prediction

        :param data_X:
        :return:
        """
        a, r = self.particle_input.feed_forward(data_X)
        for layer in self.layers:
            a, r = layer.feed_forward(a, r)
        return a

    def feed_to_layer(self, data_X, end_layer=0):
        """
        Feed data forward until given end layer. Return the resulting activation

        :param data_X: input data
        :param end_layer: the index of the ending layer
        :return: resulting activation at end layer
        """
        if len(self.layers) <= end_layer < 0:
            return None

        a, r = self.particle_input.feed_forward(data_X)
        for l, layer in enumerate(self.layers):
            a, r = layer.feed_forward(a, r)
            if l == end_layer:
                return a

        return None

    def cost(self, data_X, data_Y):
        """
        Compute the cost for all input data corresponding to expected output

        :param data_X:
        :param data_Y:
        :return:
        """
        c = self.cost_function(data_Y, self.predict(data_X))

        if self.regularizer is not None:
            c += self.regularizer.cost(self.particle_input, self.layers)

        return c

    def cost_gradient_thread(self, data_XYt):
        """
        Wrapper for multithreaded call
        :param data_XY:
        :return:
        """
        return self.cost_gradient(data_XYt[0], data_XYt[1], thread_scale=data_XYt[2])

    def cost_gradient(self, data_X, data_Y, thread_scale=1):
        """
        Computes the gradient of the cost with respect to each weight and bias in the network

        :param data_X:
        :param data_Y:
        :return:
        """

        # Output gradients
        dc_db = []
        dc_dq = []
        dc_dr = [np.zeros_like(self.particle_input.r)]
        dc_drb = []
        for l, layer in enumerate(self.layers):
            dc_db.append(np.zeros(layer.b.shape))
            dc_dq.append(np.zeros(layer.q.shape))
            dc_dr.append(np.zeros(layer.r.shape))
            dc_drb.append(np.zeros(layer.rb.shape))

        sigma_Z = []
        A_scaled, _ = self.particle_input.feed_forward(data_X)
        A = [A_scaled]  # Note: A has one more element than sigma_Z
        prev_layer_rr = self.particle_input.get_rxyz()
        for l, layer in enumerate(self.layers):
            z = layer.compute_z(A[l], prev_layer_rr)
            a = layer.compute_a(z)
            A.append(a)
            sigma_Z.append(layer.compute_da(z))
            prev_layer_rr = layer.get_rxyz()

        delta_L = self.cost_d_function(data_Y, A[-1], sigma_Z[-1])

        # IMPORTANT:
        # For threaded calls, we need to divide the cost gradient by the number threads to account for the mean being
        # taken in the cost function. When data is split, the mean is off by a factor of the number of threads.
        if thread_scale > 1:
            delta_L /= thread_scale

        # For each piece of data
        for di, data in enumerate(data_X):
            dc_db[-1] += delta_L[di]

        l = -1
        layer = self.layers[l]
        prev_layer = self.particle_input if -(l-1) > len(self.layers) else self.layers[l-1]

        Al = A[l-1]
        Al_trans = Al.transpose()
        trans_delta_L = delta_L.transpose()
        trans_sigma_Z = []
        for sz in sigma_Z:
            trans_sigma_Z.append(np.asarray(sz).transpose())

        next_delta = np.zeros((len(prev_layer.r), len(data_X)))

        # Interaction gradient
        for j in range(layer.output_size):
            trans_delta_L_j = trans_delta_L[j]
            trans_sigma_Z_l = trans_sigma_Z[l-1] if -(l-1) <= len(self.layers) \
                else np.ones((prev_layer.output_size, len(data_X)))

            qj = layer.q[j]
            for i in range(prev_layer.output_size):
                delta_r = prev_layer.r[i] - layer.rb[j]
                r = np.sqrt(np.sum(delta_r * delta_r, axis=1))
                exp_dij = layer.potential(r)
                potential = np.sum(qj * exp_dij)
                next_delta[i] += potential * trans_delta_L_j * trans_sigma_Z_l[i]

                # Charge gradient
                sum_atj = np.sum(Al_trans[i] * trans_delta_L_j)
                dc_dq[l][j] += sum_atj * exp_dij

                # Position gradient
                tmp = (-qj * layer.d_potential(r) / r).reshape((-1, 1)) * delta_r
                dc_drb[l][j] += sum_atj * tmp
                dc_dr[l-1][i] -= np.sum(sum_atj * tmp, axis=0)

        l = -1
        while -l < len(self.layers):
            l -= 1
            # Gradient computation
            layer = self.layers[l]
            prev_layer = self.particle_input if -(l-1) > len(self.layers) else self.layers[l-1]

            Al = A[l-1]
            Al_trans = Al.transpose()

            this_delta = next_delta
            next_delta = np.zeros((prev_layer.output_size, len(data_X)))
            trans_sigma_Z_l = trans_sigma_Z[l-1] if -(l-1) <= len(self.layers) \
                else np.ones((prev_layer.output_size, len(data_X)))

            # Bias gradient
            trans_delta = this_delta.transpose()
            for di, data in enumerate(data_X):
                dc_db[l] += trans_delta[di]

            # Interaction gradient
            for j in range(layer.output_size):
                this_delta_j = this_delta[j]

                qj = layer.q[j]
                for i in range(prev_layer.output_size):
                    delta_r = prev_layer.r[i] - layer.rb[j]
                    r = np.sqrt(np.sum(delta_r * delta_r, axis=1))
                    exp_dij = layer.potential(r)
                    potential = np.sum(qj * exp_dij)
                    next_delta[i] += potential * this_delta_j * trans_sigma_Z_l[i]

                    # Charge gradient
                    sum_atj = np.sum(Al_trans[i] * this_delta_j)
                    dc_dq[l][j] += sum_atj * exp_dij

                    # Position gradient
                    tmp = (-qj * layer.d_potential(r) / r).reshape((-1, 1)) * delta_r
                    dc_drb[l][j] += sum_atj * tmp
                    dc_dr[l - 1][i] -= np.sum(sum_atj * tmp, axis=0)

        return dc_db, dc_dq, dc_drb, dc_dr

    def fit(self, data_X, data_Y, optimizer):
        """
        Run the optimizer for specified number of epochs

        :param data_X:
        :param data_Y:
        :return:
        """

        return optimizer.optimize(self, data_X, data_Y)

max_stars_repo_name: awlange/brainsparks

max_issues_repo_name: awlange/brainsparks

max_forks_repo_name: awlange/brainsparks


text:
import torch
import cv2
import dlib
from collections import defaultdict
import imutils
from imutils import face_utils
from imutils.face_utils import shape_to_np
import numpy as np
from .imgutils import is_tensor_and_convert, rgb_to_bgr, img_2_gray_scale
import time
# from openface import openface

torch.manual_seed(0)

def filter_landmark_by_index(selected_landmarks, all_landmarks):

    filtered_landmarks = defaultdict(lambda: {"idxs": None})

    for k, v in all_landmarks.items():

        filtered_idxs = [j for j in v["idxs"] if j in selected_landmarks]

        if filtered_idxs != []:
            filtered_landmarks[k]["idxs"] = filtered_idxs
    return filtered_landmarks

def filter_landmark_by_group(selected_landmarks, all_landmarks):

    return {k: v for k, v in all_landmarks.items() if k in selected_landmarks}

def process_landmarks(landmarks, selected_landmarks=-1):

    # TODO: This should be a parameter in case landmarks change
    all_landmarks = {
        "mouth": {"idxs": list(range(48, 68))},
        "nose": {"idxs": list(range(27, 36))},
        "head": {"idxs": list(range(17))},
        "left_eye": {"idxs": [36, 37, 38, 39, 40, 41]},
        "left_eyebrow": {"idxs": [17, 18, 19, 20, 21]},
        "right_eye": {"idxs": [42, 43, 44, 45, 46, 47]},
        "right_eyebrow": {"idxs": [22, 23, 24, 25, 26]},
    }

    raw_num_landmarks = sum(map(lambda x: len(x["idxs"]), all_landmarks.values()))
    processed_num_landmarks = int(raw_num_landmarks)

    # Filter landmarks. TODO: This could be a method itself.
    if selected_landmarks == -1:
        filtered_landmarks = all_landmarks.copy()

    elif isinstance(selected_landmarks, list):
        if all(map(lambda x: isinstance(x, str), selected_landmarks)):
            filtered_landmarks = filter_landmark_by_group(
                selected_landmarks, all_landmarks
            )

        elif all(map(lambda x: isinstance(x, int), selected_landmarks)):
            filtered_landmarks = filter_landmark_by_index(
                selected_landmarks, all_landmarks
            )
        else:
            raise TypeError(
                "Selected landmarks must be all str, all int or an empty list"
            )
    else:
        raise TypeError("selected_landmarks can be either -1 or a list")

    processed_landmarks = defaultdict(lambda: {"idxs": None, "values": None})

    for k, v in filtered_landmarks.items():
        processed_landmarks[k] = dict(zip(v["idxs"], landmarks[v["idxs"]]))

    # Create new landmarks based on original ones. TODO: This could be a method itself.
    if "left_eye" in processed_landmarks.keys():
        processed_num_landmarks += 1
        left_eye = np.array([*processed_landmarks["left_eye"].values()])
        processed_landmarks["left_eye_center"] = {
            processed_num_landmarks - 1: np.mean(left_eye, axis=0).round().astype(int)
        }

    if "right_eye" in processed_landmarks.keys():
        processed_num_landmarks += 1
        right_eye = np.array([*processed_landmarks["right_eye"].values()])
        processed_landmarks["right_eye_center"] = {
            processed_num_landmarks - 1: np.mean(right_eye, axis=0).round().astype(int)
        }

    processed_num_landmarks += 1

    if "head" in processed_landmarks.keys():
        processed_num_landmarks += 1
        chin = np.array([*processed_landmarks["head"].values()])[6:11, :]

        processed_landmarks["chin"] = {
            processed_num_landmarks - 1: np.mean(chin, axis=0).round().astype(int)
        }

    processed_num_landmarks += 1

    return processed_landmarks

def draw_landmarks(face, img_, gray_img_, selected_landmarks, landmark_radius=1):
    """
    img must be in BGR format
    """

    img = img_.copy()
    gray_img = gray_img_.copy()

    for group_landmarks in face["landmarks"].values():
        for (x, y) in group_landmarks.values():
            cv2.circle(img, (x, y), landmark_radius, (0, 0, 255), -1)

    # print(processed_landmarks)
    return img

def get_sorted_dict_values_by_key(dictionary, reverse=False):

    return [dictionary[key] for key in sorted(dictionary.keys(), reverse=reverse)]

def processed_landmark_dict_to_array(processed_landmarks):

    flattened_landmark_dict = {
        k: v for i in processed_landmarks.values() for k, v in i.items()
    }

    flattened_sorted_landmark_array = np.vstack(
        get_sorted_dict_values_by_key(flattened_landmark_dict)
    )

    return flattened_sorted_landmark_array

class DNNDetector:

    def __init__(self,
                 model_file="models/landmarks/res10_300x300_ssd_iter_140000.caffemodel",
                 config_file="models/landmarks/deploy.prototxt.txt"):

        self.model_file = model_file
        self.config_file = config_file
        self.net = cv2.dnn.readNetFromCaffe(config_file, model_file)

    def detect(self,img):
        '''
        Only returns one bbox for the most probable face
        '''

        h, w = img.shape[:2]
        blob = cv2.dnn.blobFromImage(cv2.resize(img, (300, 300)), 1.0,
        (300, 300), (104.0, 117.0, 123.0))
        self.net.setInput(blob)
        faces = self.net.forward()
        #to draw faces on image
        max_confidence = 0
        box = None
        for i in range(faces.shape[2]):
                confidence = faces[0, 0, i, 2]
                if confidence > max_confidence:
                    box = faces[0, 0, i, 3:7] * np.array([w, h, w, h])
                    (x, y, x1, y1) = box.astype("int")
                    max_confidence = max(max_confidence,confidence)

        rect = dlib.rectangle(x, y, x1, y1)

        return rect,max_confidence

class ProcessImage:
    def __init__(self, predictor, detector, shape_predictor):

        self.predictor = predictor
        self.detector = detector
        self.shape_predictor = shape_predictor
        self.time_face_detector = 0
        self.time_landmark_detector = 0

    def __detector_wrapper(self,detector,img,gray_img):
        if isinstance(self.detector,DNNDetector):
            rect,prob = self.detector.detect(img)
        else: #Default to HOG face detector
            rect = self.detector(gray_img, 1)[0]
        return rect
    def preprocess_image(self, img_, width=800, is_bgr=False):

        # Width sets a fix img width.

        img = is_tensor_and_convert(img_)

        if not is_bgr:
            img = rgb_to_bgr(img)

        if width is not None:
            img = imutils.resize(img, width=width)

        gray_img = img_2_gray_scale(img)

        return img, gray_img

    def crop_and_align(self, img_, gray_img_, method="face_chip", size=128):

        img = img_.copy()
        gray_img = gray_img_.copy()

        # TODO: Confirm when to use gray img

        tic_face_detector = time.time()

        rect = self.__detector_wrapper(self.detector,img,gray_img)

        toc_face_detector = time.time()
        self.time_face_detector += toc_face_detector-tic_face_detector

        faces = dlib.full_object_detections()

        tic_landmarks = time.time()

        faces.append(self.predictor(gray_img, rect))

        toc_landmarks = time.time()
        self.time_landmark_detector += toc_landmarks-tic_landmarks

        if len(faces) == 0:
            raise ValueError("No face detected")
        if len(faces) > 1:
            raise ValueError("only supports one face")

        if method == "face_chip":
            face_aligned = dlib.get_face_chip(img, faces[0], size=size)

            # Recalculate landmarks instead of transforming
            face_aligned_gray = img_2_gray_scale(face_aligned)
            #rects = self.detector(face_aligned_gray, 1)
            rect2 = self.__detector_wrapper(self.detector,face_aligned,face_aligned_gray)
            landmarks_aligned = self.predictor(face_aligned_gray, rect2)
            landmarks_aligned = shape_to_np(landmarks_aligned)

        elif method == "openface":

            '''
            face_aligner = openface.AlignDlib(self.shape_predictor)

            # Use openface to calculate and perform the face alignment
            face_aligned = face_aligner.align(
                size,
                img,
                rects[0],
                landmarkIndices=openface.AlignDlib.OUTER_EYES_AND_NOSE,
            )

            # Recalculate landmarks instead of transforming
            face_aligned_gray = img_2_gray_scale(face_aligned)
            rects = self.detector(face_aligned_gray, 1)
            landmarks_aligned = self.predictor(face_aligned_gray, rects[0])
            landmarks_aligned = shape_to_np(landmarks_aligned)
            '''

            pass
        elif method == "face_aligner":
            fa = FaceAligner(self.predictor, desiredFaceWidth=size)

            # extract the ROI of the *original* face, then align the face
            # using facial landmarks
            (x, y, w, h) = imutils.face_utils.rect_to_bb(rect)
            face_aligned, landmarks_aligned = fa.align(img, gray_img, rect)

        elif method is None:
            face_aligned = img.copy()
            landmarks_aligned = self.predictor(gray_img, rect)
            landmarks_aligned = shape_to_np(landmarks_aligned)
        else:
            raise ValueError(f"method: {method} not supported")

        face_aligned_gray = img_2_gray_scale(face_aligned)

        return rect,face_aligned, face_aligned_gray, landmarks_aligned

    def process(
        self,
        img,
        width=800,
        is_bgr=False,
        method="face_chip",
        size=128,
        selected_landmarks=-1,
    ):

        img, gray_img = self.preprocess_image(img, width, is_bgr)
        rect, img, gray_img, landmarks_aligned = self.crop_and_align(
            img, gray_img, method, size
        )

        processed_landmarks = process_landmarks(landmarks_aligned, selected_landmarks)

        # Left as dict b/c we might want to store more info about the face
        face = {"landmarks": processed_landmarks}

        return rect,img, gray_img, face

FACIAL_LANDMARKS_IDXS = face_utils.FACIAL_LANDMARKS_IDXS
FACIAL_LANDMARKS_68_IDXS = face_utils.FACIAL_LANDMARKS_68_IDXS
FACIAL_LANDMARKS_5_IDXS = face_utils.FACIAL_LANDMARKS_5_IDXS

class FaceAligner:
    def __init__(
        self,
        predictor,
        desiredLeftEye=(0.35, 0.35),
        desiredFaceWidth=256,
        desiredFaceHeight=None,
    ):
        # store the facial landmark predictor, desired output left
        # eye position, and desired output face width + height
        self.predictor = predictor
        self.desiredLeftEye = desiredLeftEye
        self.desiredFaceWidth = desiredFaceWidth
        self.desiredFaceHeight = desiredFaceHeight

        # if the desired face height is None, set it to be the
        # desired face width (normal behavior)
        if self.desiredFaceHeight is None:
            self.desiredFaceHeight = self.desiredFaceWidth

    def align(self, image, gray, rect):
        # convert the landmark (x, y)-coordinates to a NumPy array
        shape = self.predictor(gray, rect)
        shape = shape_to_np(shape)

        # simple hack ;)
        if len(shape) == 68:
            # extract the left and right eye (x, y)-coordinates
            (lStart, lEnd) = FACIAL_LANDMARKS_68_IDXS["left_eye"]
            (rStart, rEnd) = FACIAL_LANDMARKS_68_IDXS["right_eye"]
        else:
            (lStart, lEnd) = FACIAL_LANDMARKS_5_IDXS["left_eye"]
            (rStart, rEnd) = FACIAL_LANDMARKS_5_IDXS["right_eye"]

        leftEyePts = shape[lStart:lEnd]
        rightEyePts = shape[rStart:rEnd]

        # compute the center of mass for each eye
        leftEyeCenter = leftEyePts.mean(axis=0).astype("int")
        rightEyeCenter = rightEyePts.mean(axis=0).astype("int")

        # compute the angle between the eye centroids
        dY = rightEyeCenter[1] - leftEyeCenter[1]
        dX = rightEyeCenter[0] - leftEyeCenter[0]
        angle = np.degrees(np.arctan2(dY, dX)) - 180

        # compute the desired right eye x-coordinate based on the
        # desired x-coordinate of the left eye
        desiredRightEyeX = 1.0 - self.desiredLeftEye[0]

        # determine the scale of the new resulting image by taking
        # the ratio of the distance between eyes in the *current*
        # image to the ratio of distance between eyes in the
        # *desired* image
        dist = np.sqrt((dX ** 2) + (dY ** 2))
        desiredDist = desiredRightEyeX - self.desiredLeftEye[0]
        desiredDist *= self.desiredFaceWidth
        scale = desiredDist / dist

        # compute center (x, y)-coordinates (i.e., the median point)
        # between the two eyes in the input image
        eyesCenter = (
            (leftEyeCenter[0] + rightEyeCenter[0]) // 2,
            (leftEyeCenter[1] + rightEyeCenter[1]) // 2,
        )

        # grab the rotation matrix for rotating and scaling the face
        M = cv2.getRotationMatrix2D(
            (int(eyesCenter[0]), int(eyesCenter[1])), angle, scale
        )

        # update the translation component of the matrix
        tX = self.desiredFaceWidth * 0.5
        tY = self.desiredFaceHeight * self.desiredLeftEye[1]
        M[0, 2] += tX - eyesCenter[0]
        M[1, 2] += tY - eyesCenter[1]

        # apply the affine transformation
        (w, h) = (self.desiredFaceWidth, self.desiredFaceHeight)
        aligned_face = cv2.warpAffine(image, M, (w, h), flags=cv2.INTER_CUBIC)

        aligned_landmarks = np.squeeze(cv2.transform(np.expand_dims(shape, 1), M), 1)

        # TODO: aligned bounding box

        # return the aligned face
        return aligned_face, aligned_landmarks


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

import uuid
from itertools import chain

import typing
import pyspark

from pyspark.rddsampler import RDDSamplerBase

from fate_arch.abc import CTableABC
from fate_arch.common import log, hdfs_utils, hive_utils
from fate_arch.common.profile import computing_profile
from fate_arch.computing.spark._materialize import materialize, unmaterialize
from scipy.stats import hypergeom
from fate_arch.computing._type import ComputingEngine

LOGGER = log.getLogger()

class Table(CTableABC):
    def __init__(self, rdd):
        self._rdd: pyspark.RDD = rdd
        self._engine = ComputingEngine.SPARK

    @property
    def engine(self):
        return self._engine

    def __getstate__(self):
        pass

    def __del__(self):
        try:
            unmaterialize(self._rdd)
            del self._rdd
        except BaseException:
            return

    @computing_profile
    def save(self, address, partitions, schema, **kwargs):
        from fate_arch.common.address import HDFSAddress

        if isinstance(address, HDFSAddress):
            self._rdd.map(lambda x: hdfs_utils.serialize(x[0], x[1])).repartition(
                partitions
            ).saveAsTextFile(f"{address.name_node}/{address.path}")
            schema.update(self.schema)
            return

        from fate_arch.common.address import HiveAddress, LinkisHiveAddress

        if isinstance(address, (HiveAddress, LinkisHiveAddress)):
            # df = (
            #     self._rdd.map(lambda x: hive_utils.to_row(x[0], x[1]))
            #     .repartition(partitions)
            #     .toDF()
            # )
            LOGGER.debug(f"partitions: {partitions}")
            _repartition = self._rdd.map(lambda x: hive_utils.to_row(x[0], x[1])).repartition(partitions)
            _repartition.toDF().write.saveAsTable(f"{address.database}.{address.name}")
            schema.update(self.schema)
            return

        from fate_arch.common.address import LocalFSAddress

        if isinstance(address, LocalFSAddress):
            self._rdd.map(lambda x: hdfs_utils.serialize(x[0], x[1])).repartition(
                partitions
            ).saveAsTextFile(address.path)
            schema.update(self.schema)
            return

        raise NotImplementedError(
            f"address type {type(address)} not supported with spark backend"
        )

    @property
    def partitions(self):
        return self._rdd.getNumPartitions()

    @computing_profile
    def map(self, func, **kwargs):
        return from_rdd(_map(self._rdd, func))

    @computing_profile
    def mapValues(self, func, **kwargs):
        return from_rdd(_map_value(self._rdd, func))

    @computing_profile
    def mapPartitions(
        self, func, use_previous_behavior=True, preserves_partitioning=False, **kwargs
    ):
        if use_previous_behavior is True:
            LOGGER.warning(
                f"please use `applyPartitions` instead of `mapPartitions` "
                f"if the previous behavior was expected. "
                f"The previous behavior will not work in future"
            )
            return self.applyPartitions(func)
        return from_rdd(
            self._rdd.mapPartitions(func, preservesPartitioning=preserves_partitioning)
        )

    @computing_profile
    def mapReducePartitions(self, mapper, reducer, **kwargs):
        return from_rdd(self._rdd.mapPartitions(mapper).reduceByKey(reducer))

    @computing_profile
    def applyPartitions(self, func, **kwargs):
        return from_rdd(_map_partitions(self._rdd, func))

    @computing_profile
    def glom(self, **kwargs):
        return from_rdd(_glom(self._rdd))

    @computing_profile
    def sample(
        self,
        *,
        fraction: typing.Optional[float] = None,
        num: typing.Optional[int] = None,
        seed=None,
    ):
        if fraction is not None:
            return from_rdd(
                self._rdd.sample(fraction=fraction, withReplacement=False, seed=seed)
            )

        if num is not None:
            return from_rdd(_exactly_sample(self._rdd, num, seed=seed))

        raise ValueError(
            f"exactly one of `fraction` or `num` required, fraction={fraction}, num={num}"
        )

    @computing_profile
    def filter(self, func, **kwargs):
        return from_rdd(_filter(self._rdd, func))

    @computing_profile
    def flatMap(self, func, **kwargs):
        return from_rdd(_flat_map(self._rdd, func))

    @computing_profile
    def reduce(self, func, **kwargs):
        return self._rdd.values().reduce(func)

    @computing_profile
    def collect(self, **kwargs):
        #         return iter(self._rdd.collect())
        return self._rdd.toLocalIterator()

    @computing_profile
    def take(self, n=1, **kwargs):
        _value = self._rdd.take(n)
        if kwargs.get("filter", False):
            self._rdd = self._rdd.filter(lambda xy: xy not in [_xy for _xy in _value])
        return _value

    @computing_profile
    def first(self, **kwargs):
        return self.take(1)[0]

    @computing_profile
    def count(self, **kwargs):
        return self._rdd.count()

    @computing_profile
    def join(self, other: "Table", func=None, **kwargs):
        return from_rdd(_join(self._rdd, other._rdd, func=func))

    @computing_profile
    def subtractByKey(self, other: "Table", **kwargs):
        return from_rdd(_subtract_by_key(self._rdd, other._rdd))

    @computing_profile
    def union(self, other: "Table", func=None, **kwargs):
        return from_rdd(_union(self._rdd, other._rdd, func))

def from_hdfs(paths: str, partitions, in_serialized=True, id_delimiter=None):
    # noinspection PyPackageRequirements
    from pyspark import SparkContext

    sc = SparkContext.getOrCreate()
    fun = hdfs_utils.deserialize if in_serialized else lambda x: (x.partition(id_delimiter)[0],
                                                                  x.partition(id_delimiter)[2])
    rdd = materialize(
        sc.textFile(paths, partitions)
        .map(fun)
        .repartition(partitions)
    )
    return Table(rdd=rdd)

def from_localfs(paths: str, partitions, in_serialized=True, id_delimiter=None):
    # noinspection PyPackageRequirements
    from pyspark import SparkContext

    sc = SparkContext.getOrCreate()
    fun = hdfs_utils.deserialize if in_serialized else lambda x: (x.partition(id_delimiter)[0],
                                                                  x.partition(id_delimiter)[2])
    rdd = materialize(
        sc.textFile(paths, partitions)
        .map(fun)
        .repartition(partitions)
    )
    return Table(rdd=rdd)

def from_hive(tb_name, db_name, partitions):
    from pyspark.sql import SparkSession

    session = SparkSession.builder.enableHiveSupport().getOrCreate()
    rdd = materialize(
        session.sql(f"select * from {db_name}.{tb_name}")
        .rdd.map(hive_utils.from_row)
        .repartition(partitions)
    )
    return Table(rdd=rdd)

def from_rdd(rdd):
    rdd = materialize(rdd)
    return Table(rdd=rdd)

def _fail_on_stopiteration(fn):
    # noinspection PyPackageRequirements
    from pyspark import util

    return util.fail_on_stopiteration(fn)

def _map(rdd, func):
    def _fn(x):
        return func(x[0], x[1])

    def _func(_, iterator):
        return map(_fail_on_stopiteration(_fn), iterator)

    return rdd.mapPartitionsWithIndex(_func, preservesPartitioning=False)

def _map_value(rdd, func):
    def _fn(x):
        return x[0], func(x[1])

    def _func(_, iterator):
        return map(_fail_on_stopiteration(_fn), iterator)

    return rdd.mapPartitionsWithIndex(_func, preservesPartitioning=True)

def _map_partitions(rdd, func):
    def _func(_, iterator):
        return [(str(uuid.uuid1()), func(iterator))]

    return rdd.mapPartitionsWithIndex(_func, preservesPartitioning=False)

def _join(rdd, other, func=None):
    num_partitions = max(rdd.getNumPartitions(), other.getNumPartitions())
    rtn_rdd = rdd.join(other, numPartitions=num_partitions)
    if func is not None:
        rtn_rdd = _map_value(rtn_rdd, lambda x: func(x[0], x[1]))
    return rtn_rdd

def _glom(rdd):
    def _func(_, iterator):
        yield list(iterator)

    return rdd.mapPartitionsWithIndex(_func)

def _exactly_sample(rdd, num: int, seed: int):
    split_size = rdd.mapPartitionsWithIndex(
        lambda s, it: [(s, sum(1 for _ in it))]
    ).collectAsMap()
    total = sum(split_size.values())

    if num > total:
        raise ValueError(f"not enough data to sample, own {total} but required {num}")
    # random the size of each split
    sampled_size = {}
    for split, size in split_size.items():
        sampled_size[split] = hypergeom.rvs(M=total, n=size, N=num)
        total = total - size
        num = num - sampled_size[split]

    return rdd.mapPartitionsWithIndex(
        _ReservoirSample(split_sample_size=sampled_size, seed=seed).func,
        preservesPartitioning=True,
    )

class _ReservoirSample(RDDSamplerBase):
    def __init__(self, split_sample_size, seed):
        RDDSamplerBase.__init__(self, False, seed)
        self._split_sample_size = split_sample_size
        self._counter = 0
        self._sample = []

    def func(self, split, iterator):
        self.initRandomGenerator(split)
        size = self._split_sample_size[split]
        for obj in iterator:
            self._counter += 1
            if len(self._sample) < size:
                self._sample.append(obj)
                continue

            randint = self._random.randint(1, self._counter)
            if randint <= size:
                self._sample[randint - 1] = obj

        return self._sample

def _filter(rdd, func):
    def _fn(x):
        return func(x[0], x[1])

    def _func(_, iterator):
        return filter(_fail_on_stopiteration(_fn), iterator)

    return rdd.mapPartitionsWithIndex(_func, preservesPartitioning=True)

def _subtract_by_key(rdd, other):
    return rdd.subtractByKey(other, rdd.getNumPartitions())

def _union(rdd, other, func):
    num_partition = max(rdd.getNumPartitions(), other.getNumPartitions())
    if func is None:
        return rdd.union(other).coalesce(num_partition)
    else:

        def _func(pair):
            iter1, iter2 = pair
            val1 = list(iter1)
            val2 = list(iter2)
            if not val1:
                return val2[0]
            if not val2:
                return val1[0]
            return func(val1[0], val2[0])

        return _map_value(rdd.cogroup(other, num_partition), _func)

def _flat_map(rdd, func):
    def _fn(x):
        return func(x[0], x[1])

    def _func(_, iterator):
        return chain.from_iterable(map(_fail_on_stopiteration(_fn), iterator))

    return rdd.mapPartitionsWithIndex(_func, preservesPartitioning=False)

reason: from scipy


text:
import numpy as np
import struct
import random

# pip install Vose-Alias-Method
from vose_sampler import VoseAlias

# needs to match the value in (afl's) random.h
ALIAS_MAX = 1<<30

def weights2alias(w):
    # np.array(dtype=float64) -> dict of {outcome:proportion}
    dist = dict(enumerate(w))

    # Convert a distribution (dict of {outcome:proportion}) into an alias table (alias, prob)
    va = VoseAlias(dist)
    length = len(va.table_prob_list)

    alias = np.arange(length, dtype=np.uint32)
    for k,v in va.table_alias.items():
        alias[k] = v

    prob = np.zeros(length, dtype=np.uint32)
    for k,v in va.table_prob.items():
        prob[k] = v * ALIAS_MAX

    return alias, prob

def write_alias(alias, prob, path):
    with open(path, 'wb') as f:
        f.write(struct.pack('@I', len(alias)))
        f.write(alias.tobytes())
        f.write(prob.tobytes())

def read_alias(path):
    def read_fmt(fmt, f):
        return struct.unpack(fmt, f.read(struct.calcsize(fmt)))
    with open(path, 'rb') as f:
        length = read_fmt('@I', f)[0]
        print(length)
        alias = np.frombuffer(f.read(length * 4), dtype=np.uint32, count=length)
        prob = np.frombuffer(f.read(length * 4), dtype=np.uint32, count=length)
    return alias, prob

# This class can be used to read and sample from the saved alias tables in
# out/queue/.state/offset_weights/
class AliasTable(object):

    def __init__(self, alias, prob):
        super().__init__()
        self.table_prob_list = list(range(prob.size))
        self.table_alias = alias
        self.table_prob = prob

    # this routine is from Vose-Alias-Method, with ALIAS_MAX modification
    def alias_generation(self):
        """ Return a random outcome from the distribution. """
        # Determine which column of table_prob to inspect
        col = random.choice(self.table_prob_list)

        # Determine which outcome to pick in that column
        # if self.table_prob[col] >= random.uniform(0,1):
        if self.table_prob[col] >= random.randint(0, ALIAS_MAX-1):
            return col
        else:
            return self.table_alias[col]

    # this routine is from Vose-Alias-Method
    def sample_n(self, size):
        """ Return a sample of size n from the distribution."""
        # Ensure a non-negative integer as been specified
        n = int(size)
        if n <= 0:
            raise ValueError("Please enter a non-negative integer for the number of samples desired: %d" % n)

        return [self.alias_generation() for i in range(n)]

max_stars_repo_name: sandialabs/wafl

max_issues_repo_name: sandialabs/wafl

max_forks_repo_name: sandialabs/wafl


text:
"""
Script to restore subject information from the HANDS 2017 training set. It is known that each sample
in the dataset belongs to one out of 5 subjects. It can be assumed that each subject has a slightly
different hand shape. By running a simple clustering algorithm on the bone lengths the mapping from
subjects to sample indices can be restored.
"""

import os
import torch
import scipy.cluster.vq as scikmeans

from data_utils import datasets

dataset_name = 'HANDS17_DPREN_all'

dataset = datasets.PairedPoseDataset(dataset_name, use_preset=True)
all_labels = dataset[:].labels

wrist_to_tmcp_lengths = torch.norm(all_labels[:, 0] - all_labels[:, 1], dim=1).reshape(-1, 1)

whitened_lengths = scikmeans.whiten(wrist_to_tmcp_lengths)
means, _ = scikmeans.kmeans(whitened_lengths, 5, iter=10)
mapping, _ = scikmeans.vq(whitened_lengths, means)

output_file_path = os.path.join('data', 'clusterings', dataset_name + '_subjects.pt')
torch.save(torch.from_numpy(mapping), output_file_path)

reason: import scipy


text:
import config
import datetime
import json
import numpy as np
import pandas as pd
import requests
import sqlite3
import time
import yaml

# Special treatment for some claims

# LBRY Inc channels
inc = set(["f3da2196b5151570d980b34d311ee0973225a68e",
           "70b8a88fc6e5ce9e4d6e8721536688484ecd79f4",
           "3fda836a92faaceedfe398225fb9b2ee2ed1f01a",
           "e48d2b50501159034f68d53321f67b8aa5b1d771",
           "e8fed337dc4ee260f4bcfa6d24ae1e4dd75c2fb3"])

# LBRY Social channels
ls = set(["5bd299a92e7b31865d2bb3e2313402edaca41a94",
          "f8d6eccd887c9cebd36b1d42aa349279b7f5c3ed",
          "e11e2fc3056137948d2cc83fb5ca2ce9b57025ec",
          "1ba5acff747615510cf3f6089f54d5de669ad94f",
          "4506db7fb52d3ec5d3a024c870bf86fc35f7b6a3",
          "36b7bd81c1f975878da8cfe2960ed819a1c85bb5",
          "e5f33f22ef656cb1595140409850a04d60aa474b",
          "631ca9fce459f1116ae5317486c7f4af69554742",
          "4caa1f92fb477caed1ce07cb7762a2249050a59c",
          "56e86eb938c0b93beccde0fbaaead65755139a10",
          "60ea26a907f25bcbbc8215007eef2bf0fb846f5c",
          "d0174cf90b6ec4e26ee2fc013714b0803dec5dd1",
          "3849a35ae6122e0b7a035c2ba66e97b9e4ab9efa" ])

# Given mature tag by us
manual_mature = set(["f24ab6f03d96aada87d4e14b2dac4aa1cee8d787",
                 "fd4b56c7216c2f96db4b751af68aa2789c327d48"])

# Grey list (quietly disable link)
grey_list = set(["ca8cfeb5b6660a0b8874593058178b7ce6af5fed",
              "6c1119f18fd7a15fc7535fcb9eec3aa22af66b6b",
              "3097b755d3b8731e6103cc8752cb1b6c79da3b85",
              "11c2f6bb38f69a25dea3d0fbef67e2e3a83a1263",
              "7acf8b2fcd212afa2877afe289309a20642880c4",
              "b01a44af8b71c0c2001a78303f319ca960d341cf",
              "bc89d67d9f4d0124c347fd2c4a04e1696e8ba8b1",
              "14fcd92ad24c1f1bc50f6cbc1e972df79387d05c",
              "977cd1c90eefe4c9831f5c93b2359202733a9c2e" ])

# DMCA'd channels + rewards scammers (do not appear)
black_list = set([ "98c39de1c681139e43131e4b32c2a21272eef06e",
                "9ced2a722e91f28e9d3aea9423d34e08fb11e3f4",
                "d5557f4c61d6725f1a51141bbee43cdd2576e415",
                "35100b76e32aeb2764d334186249fa1b90d6cd74",
                "f2fe17fb1c62c22f8319c38d0018726928454112",
                "17db8343914760ba509ed1f8c8e34dcc588614b7",
                "06a31b83cd38723527861a1ca5349b0187f92193",
                "9b7a749276c69f39a2d2d76ca4353c0d8f75217d",
                "b1fa196661570de64ff92d031116a2985af6034c",
                "4e5e34d0ab3cae6f379dad75afadb0c1f683d30f",
                "86612188eea0bda3efc6d550a7ad9c96079facff",
                "00aa9655c127cccb2602d069e1982e08e9f96636",
                "4f2dba9827ae28a974fbc78f1b12e67b8e0a32c9",
                "c133c44e9c6ee71177f571646d5b0000489e419f",
                "eeb3c6452b240a9f6a17c06887547be54a90a4b9",
                "f625ef83a3f34cac61b6b3bdef42be664fd827da",
                "ed77d38da413377b8b3ee752675662369b7e0a49",
                "481c95bd9865dc17770c277ae50f0cc306dfa8af",
                "3c5aa133095f97bb44f13de7c85a2a4dd5b4fcbe",
                "bd6abead1787fa94722bd7d064f847de76de5655",
                "6114b2ce20b55c40506d4bd3f7d8f917b1c37a75",
                "0c65674e28f2be555570c5a3be0c3ce2eda359d1",
                "3395d03f379888ffa789f1fa45d6619c2037e3de",
                "cd31c9ddea4ac4574df50a1f84ee86aa17910ea2",
                "9d48c8ab0ad53c392d4d6052daf5f8a8e6b5a185",
                "51fbdb73893c1b04a7d4c4465ffcd1138abc9e93",
                "5183307ce562dad27367bdf94cdafde38756dca7",
                "56dca125e775b2fe607d3d8d6c29e7ecfa3cbd96",
                "a58926cb716c954bdab0187b455a63a2c592310e",
                "aa83130864bf22c66934c1af36182c91219233aa",
                "f3c1fda9bf1f54710b62ffe4b14be6990288d9ff",
                "6291b3b53dde4160ce89067281300585bdf51905",
                "eeef31480a14684a95898ecd3bcf3a5569e41a28",
                "8b8b3c8cd3e8364c37067b80bd5a20c09a0a0094",
                "725189cd101ff372edbce1c05ef04346864d3254",
                "35100b76e32aeb2764d334186249fa1b90d6cd74",
                "47beabb163e02e10f99838ffc10ebc57f3f13938",
                "e0bb55d4d6aec9886858df8f1289974e673309c7",
                "242734793097302d33b6a316c9db8d17b4beb18e",
                "71d3256c267ccc875df366258b9eff4766d6cb57",
                "dee09cad16900936d6af97154a6510a09587ad42",
                "357ce885e22f2a7bd426ac36224722d64fc90ce6",
                "c3ab2407e295cd267ced06d1fad2ed09b8d5643e",
                "37b96ce8ae7a5564174111573105ee7efe4cd2fc" ])

def subscriber_counts(num=200, preview=False):
    """
    Get subscriber counts for all channels.
    """
    now = time.time()

    # Open previous JSON
    f = open("json/subscriber_counts.json")
    old = json.load(f)
    f.close()

    # Get auth token
    f= open("secrets.yaml")
    auth_token = yaml.load(f, Loader=yaml.SafeLoader)["auth_token"]
    f.close()

    # Create a dict from the old JSON, where the claim_id can return
    # the subscribers and the rank
    old_dict = {}
    for i in range(len(old["ranks"])):
        old_dict[old["claim_ids"][i]] = (old["subscribers"][i], old["ranks"][i])

    # Open claims.db
    conn = sqlite3.connect(config.claims_db_file)
    c = conn.cursor()
    query = "select claim_name, claim_id, claim_hash from claim where claim_type = 2;"
    vanity_names = []
    claim_ids = []
    subscribers = []

    # Iterate over query results
    i = 0
    for row in c.execute(query):
        vanity_names.append(row[0])
        claim_ids.append(row[1])
        i = i + 1

    vanity_names = np.array(vanity_names)
    claim_ids = np.array(claim_ids)

    # Now get number of claims in each channel
    query = \
"""
select c2.claim_id claim_ids, count(*) num_claims
    from claim c1 inner join claim c2 on c2.claim_hash = c1.channel_hash
    group by c2.claim_hash
    having num_claims > 0;
"""

    claims_with_content = {}
    k = 0
    for row in c.execute(query):
        claims_with_content[row[0]] = None
        k += 1
        print("Getting channels with content...found {k} so far.".format(k=k))

    start = time.time()

    include = np.zeros(len(claim_ids), dtype=bool)
    for i in range(len(claim_ids)):
        include[i] = (claim_ids[i] in claims_with_content) and \
                            claim_ids[i] not in black_list

    vanity_names = vanity_names[include]
    claim_ids = claim_ids[include]

    k = 0
    while True:
        """
        Go in batches of 100 with a pause in between
        """
        time.sleep(3.0)

        # Cover a certain range of channels
        start = 100*k
        end = 100*(k+1)
        final = end >= len(claim_ids)
        if final:
            end = len(claim_ids)

        # Attempt the request until it succeeds
        while True:

            # Prepare the request to the LBRY API
            url = "https://api.lbry.com/subscription/sub_count?auth_token=" +\
                        auth_token + "&claim_id="
            for i in range(start, end):
                url += claim_ids[i] + ","
            url = url[0:-1] # No final comma

            f = open("url.txt", "w")
            f.write(url)
            f.close()

            try:
                # Do the request
                result = requests.get(url)
                result = result.json()
                break
            except:
                time.sleep(3.0)
                pass

        # Get sub counts from the result and put them in the subscribers list
        for x in result["data"]:
            subscribers.append(x)
            i = len(subscribers)-1

        print("Processed {end} channels.".format(end=end))
        if final:
            break
        k += 1

    # Sort by number of subscribers then by vanity name.
    # Zip subs with name
    s_n = []
    indices = []
    for i in range(len(vanity_names)):
        s_n.append((subscribers[i], vanity_names[i]))
        indices.append(i)
    indices = sorted(indices, key=lambda x: (s_n[x][0], s_n[x][1]))[::-1]

    vanity_names = np.array(vanity_names)[indices]
    claim_ids = np.array(claim_ids)[indices]
    subscribers = np.array(subscribers)[indices]

    # Put the top 200 into the dict
    my_dict = {}
    my_dict["unix_time"] = now
    my_dict["human_time_utc"] = str(datetime.datetime.utcfromtimestamp(int(now))) + " UTC"
    my_dict["old_unix_time"] = old["unix_time"]
    my_dict["old_human_time_utc"] = old["human_time_utc"]
    my_dict["interval_days"] = np.round((my_dict["unix_time"]\
                                        - my_dict["old_unix_time"])/86400.0, 2)
    my_dict["ranks"] = []
    my_dict["vanity_names"] = []
    my_dict["claim_ids"] = []
    my_dict["subscribers"] = []
    my_dict["change"] = []
    my_dict["rank_change"] = []
    my_dict["is_nsfw"] = []
    my_dict["grey"] = []
    my_dict["ls"] = []
    my_dict["inc"] = []

    for i in range(num):
        my_dict["ranks"].append(i+1)
        my_dict["vanity_names"].append(vanity_names[i])
        my_dict["claim_ids"].append(claim_ids[i])
        my_dict["subscribers"].append(int(subscribers[i]))
        my_dict["is_nsfw"].append(False)

        # Compute subscribers change
        my_dict["change"].append(None)
        my_dict["rank_change"].append(None)
        try:
            my_dict["change"][-1] = int(subscribers[i]) - \
                                        old_dict[claim_ids[i]][0]
            my_dict["rank_change"][-1] = old_dict[claim_ids[i]][1] - \
                                            int(my_dict["ranks"][-1])
        except:
            pass

        # Mark some channels NSFW manually
        if my_dict["claim_ids"][-1] in manual_mature:
            my_dict["is_nsfw"][-1] = True
        else:
            # Do SQL queries to see if there's a mature tag
            query = "SELECT tag.tag FROM claim INNER JOIN tag ON tag.claim_hash = claim.claim_hash WHERE claim_id = '"
            query += claim_ids[i] + "';"

            for row in c.execute(query):
                if row[0].lower() == "mature":
                    my_dict["is_nsfw"][-1] = True

        # Grey list
        my_dict["grey"].append(my_dict["claim_ids"][-1] in grey_list)

        # LS list
        my_dict["ls"].append(my_dict["claim_ids"][-1] in ls)
        my_dict["inc"].append(my_dict["claim_ids"][-1] in inc)

    if preview:
        f = open("json/subscriber_counts_preview.csv", "w")
        # Create data frame and make CSV
        df = pd.DataFrame()
        df["ranks"] = my_dict["ranks"]
        df["vanity_names"] = my_dict["vanity_names"]
        df["claim_ids"] = my_dict["claim_ids"]
        df["is_nsfw"] = my_dict["is_nsfw"]
        df["grey"] = my_dict["grey"]
        df["ls"] = my_dict["ls"]
        df["inc"] = my_dict["inc"]
        df["followers"] = my_dict["subscribers"]
        df["change"] = my_dict["change"]
        df["rank_change"] = my_dict["rank_change"]
        df.to_csv("json/subscriber_counts_preview.csv", index=False)

    else:
        f = open("json/subscriber_counts.json", "w")
        import update_rss
        update_rss.update(my_dict["human_time_utc"])
        f.write(json.dumps(my_dict, indent=4))
        f.close()

    conn.close()

# Main loop
if __name__ == "__main__":

    # Needs an initial JSON file to bootstrap from
    hour = 3600.0
    day = 24*hour
    week = 7*day

    f = open("json/subscriber_counts.json")
    t = json.load(f)["unix_time"]
    f.close()

    # Update frequency
    interval = 0.5*week

    while True:
        gap = time.time() - t

        msg = "{d} days until next update.".format(d=(interval - gap)/day)
        print(msg + "        ", end="\r", flush=True)
        time.sleep(1.0 - time.time()%1)

        if gap >= interval:
            subscriber_counts()

            f = open("json/subscriber_counts.json")
            t = json.load(f)["unix_time"]
            f.close()


text:
import csv
import glob

import numpy as np

from source import utils

class MetadataService(object):

    @staticmethod
    def get_all_files():
        project_root = str(utils.get_project_root())
        return glob.glob(project_root + "/data/mesa/polysomnography/edfs/*edf")

    @staticmethod
    def get_metadata_dictionary():
        metadata_dictionary = {}
        ref_dict = dict()
        ref_dict['ahi'] = 'ahiu35'
        ref_dict['age'] = 'sleepage5c'
        ref_dict['gender'] = 'gender1'  # 0 Female, 1 Male
        ref_dict['tst'] = 'slpprdp5'
        ref_dict['tib'] = 'time_bed5'
        ref_dict['waso'] = 'waso5'
        ref_dict['slp_eff'] = 'slp_eff5'
        ref_dict['time_rem'] = 'timerem5'
        ref_dict['time_n1'] = 'timest15'
        ref_dict['time_n2'] = 'timest25'
        ref_dict['time_n34'] = 'timest345'

        with open('../mesa/mesa-sleep-dataset-0.3.csv') as csv_file:
            csv_reader = csv.reader(csv_file, delimiter=',')
            is_first_row = True
            for row in csv_reader:
                if is_first_row:
                    is_first_row = False
                    count = 0
                    for col in row:
                        if col in ref_dict.values():
                            keys = list(ref_dict.keys())
                            key = keys[list(ref_dict.values()).index(col)]
                            ref_dict[key] = count  # Replace with index number where data is
                        count = count + 1
                else:
                    subject_dict = {}
                    for key in ref_dict:
                        val = row[ref_dict[key]]
                        if len(val) > 0:
                            subject_dict[key] = float(val)
                        else:
                            subject_dict[key] = np.nan
                    metadata_dictionary[int(row[0])] = subject_dict

        return metadata_dictionary

    @staticmethod
    def print_table(subject_ids):
        all_files = MetadataService.get_all_files()
        metadata_dictionary = MetadataService.get_metadata_dictionary()
        ahis = []
        ages = []
        genders = []
        tst = []
        tib = []
        waso = []
        slp_eff = []
        time_rem = []
        time_nrem = []

        for subject_index in subject_ids:  # Get all metadata for subjects
            subject = int(all_files[subject_index][-8:-4])
            subject_dict = metadata_dictionary[subject]
            ahis.append(subject_dict['ahi'])
            ages.append(subject_dict['age'])
            genders.append(subject_dict['gender'])
            tst.append(subject_dict['tst'])
            tib.append(subject_dict['tib'])
            waso.append(subject_dict['waso'])
            slp_eff.append(subject_dict['slp_eff'])
            time_rem.append(subject_dict['time_rem'])
            time_nrem.append(
                float(subject_dict['time_n1']) + float(subject_dict['time_n2']) + float(subject_dict['time_n34']))

        ahis = np.array(ahis)
        ages = np.array(ages)
        genders = np.array(genders)
        tst = np.array(tst)
        tib = np.array(tib)
        waso = np.array(waso)
        slp_eff = np.array(slp_eff)

        print('N women: ' + str(len(genders) - np.count_nonzero(genders)))

        latex = True
        if latex:
            print(
                '\\begin{table} \\caption{Sleep summary statistics - MESA subcohort}  \\label{tab:mesa} \\small  \\begin{tabularx}{\\columnwidth}{X | X | X  }\\hline Parameter & Mean (SD) & Range \\\\ \\hline')
        else:
            print('Parameter, Mean (SD), Range')

        print(MetadataService.data_to_line('Age (years)', ages, latex))
        print(MetadataService.data_to_line('TST (minutes)', tst, latex))
        print(MetadataService.data_to_line('TIB (minutes)', tib, latex))
        print(MetadataService.data_to_line('WASO (minutes)', waso, latex))
        print(MetadataService.data_to_line('SE (\%)', slp_eff, latex))
        print(MetadataService.data_to_line('AHI', ahis, latex))

        if latex:
            print('\\end{tabularx} \\end{table}')

    @staticmethod
    def data_to_line(name, data, latex=False):
        if latex:
            return name + ' & ' + str(round(np.mean(data), 2)) + ' (' + str(round(np.std(data), 2)) + ') & ' \
                   + str(round(np.min(data), 2)) + '-' + str(round(np.max(data), 2)) + '\\\\'
        else:
            return name + ',' + str(round(np.mean(data), 2)) + ' (' + str(round(np.std(data), 2)) + '),' \
                   + str(round(np.min(data), 2)) + '-' + str(round(np.max(data), 2))


text:
import tensorflow.compat.v2 as tf
import numpy as np
from . import core

epsilon = tf.keras.backend.epsilon()

def _counts_to_probs(ref_counts, tau, alphabet_size, dtype=tf.float64):
    """Get transition probabilities from reference.

    Parameters
    ----------
    ref_counts : dtype
        A tensor of shape [A1, ..., An, alphabet_size+1] of transition counts from the reference.
        All stops in the reference should be set to 0 so that they may instead be modelled
        by the net function.
    tau : dtype
        A positive constant describing the error rate 1-exp(-tau). Error is added by a Jukes-Cantor model.
    alphabet_size : int
    dtype : dtype, default = tf.float64

    Returns
    -------
    transition_probabilities : dtype
        A tensor of shape [A1, ..., An, alphabet_size+1] normalized in the last dimension and with
        zero probability of stopping.
    """
    # ref_counts have epsilon added to them except at stop, so normalization sums to 1.
    norm = tf.linalg.normalize(ref_counts, ord=1, axis=-1)
    # norm is only 0 at stop (by addition of epsilon when loading ref_counts).
    shape = tf.convert_to_tensor(np.r_[np.ones(alphabet_size), 0], dtype=dtype)
    return ((1/alphabet_size)*shape + tf.exp(-tau)*(norm[0]-((1/alphabet_size)*shape)))

def _make_ref_ar_func(lag, alphabet_size, make_net_func, af_kwargs, dtype=tf.float64):
    """Make an autoregressive function that uses the reference.

    Parameters
    ----------
    lag : int
    alphabet_size : int
    make_net_func : function
        Takes lag, alphabet_size, af_kwargs, dtype and returns an ar_func.
    af_kwargs : dict
        Keyword arguments for particular ar_func. For example, filter_width.
    dtype : dtype

    Returns
    -------
    ar_func : function, default = None
        A function that takes a tensor of shape [A1, ..., An, lag, alphabet_size+1]
        of dtype and returns a tensor of shape [A1, ..., An, alphabet_size+1] of dtype
        of transition probabilities for each kmer. The autoregressive function.
        Set to 0 if None.
    params : list
        List of parameters as tensorflow variables.
    """
    net_weight_signed = tf.Variable(-np.log(100), dtype=dtype, name='net_weight_signed')
    tau_signed = tf.Variable(np.log(1/30), dtype=dtype, name='tau_signed')
    net_func, ar_func_params = make_net_func(lag, alphabet_size, **af_kwargs, dtype=dtype)

    def ar_func(kmer_seqs, ref_counts):
        nw = tf.math.exp(net_weight_signed)
        tau = tf.math.exp(tau_signed)
        return ((nw*net_func(kmer_seqs)
                 + _counts_to_probs(ref_counts, tau, alphabet_size, dtype=dtype))
                / (nw + 1))
    return ar_func, ([tau_signed, net_weight_signed] + ar_func_params)

def _bear_kmer_counts(kmer_seqs, kmer_total_counts, ref_counts,
                      condition_trans_counts=None, h=None, ar_func=None):
    """Get random variable of kmer transition counts conditioned on observed counts.

    Parameters
    ----------
    kmer_seqs : dtype
        A tensor of shape [A1, ..., An, lag, alphabet_size+1] of one-hot encoded kmers.
    kmer_total_counts : int
        A tensor of shape [A1, ..., An] of counts of each kmer.
    condition_trans_counts : int, default = None
        A tensor of shape [A1, ..., An, alphabet_size+1] of transition counts of each
        kmer to condition on. set to 0 if None.
    h : dtype, default = None
        A positive constant of the :math:`h` parameter from the BEAR model. Set to 1 if None.
    ar_func : function, default = None
        A function that takes a tensor of shape [A1, ..., An, lag, alphabet_size+1]
        of dtype and returns a tensor of shape [A1, ..., An, alphabet_size+1] of dtype
        of transition probabilities for each kmer. The autoregressive function.
        Set to 0 if None.

    Returns
    -------
    x : tensorflow probability distribution
        Distribution of kmer transition counts for a BEAR model.
    """
    dtype = kmer_seqs.dtype
    if condition_trans_counts is None:
        condition_trans_counts = tf.constant(0., dtype)
    if h is None or ar_func is None:
        h = tf.constant(1., dtype)

        def ar_func(x, y):
            return tf.constant(0., dtype)
    concentrations = ar_func(kmer_seqs, ref_counts) / h + condition_trans_counts + epsilon
    x = core.tfpDirichletMultinomialPerm(kmer_total_counts, concentrations, name='x')
    return x

def _ar_kmer_counts(kmer_seqs, kmer_total_counts, ref_counts, ar_func):
    """Get Random variable of kmer transition counts conditioned on observed counts.

    Parameters
    ----------
    kmer_seqs : dtype
        A tensor of shape [A1, ..., An, lag, alphabet_size+1] of one-hot encoded kmers.
    kmer_total_counts : int
        A tensor of shape [A1, ..., An] of counts of each kmer.
    ar_func : function, default = None
        A function that takes a tensor of shape [A1, ..., An, lag, alphabet_size+1]
        of dtype and another of shape [A1, ..., An, alphabet_size+1] and returns a
        tensor of shape [A1, ..., An, alphabet_size+1] of dtype of transition probabilities
        for each kmer. The autoregressive function.

    Returns
    -------
    x : tensorflow probability distribution
        Distribution of kmer transition counts for an AR model.
    """
    probs = ar_func(kmer_seqs, ref_counts) + epsilon
    x = core.tfpMultinomialPerm(kmer_total_counts, probs, name='x')
    return x

def _create_params(lag, alphabet_size, make_ar_func,
                   af_kwargs, dtype=tf.float64):
    """Define and get parameters of BEAR or AR distribution.

    Parameters
    ----------
    lag : int
    alphabet_size : int
    make_ar_func : function
        Takes lag, alphabet_size, af_kwargs, dtype and returns a ar_func.
    af_kwargs : dict
        Keyword arguments for particular ar_func. For example, filter_width.
    dtype : tensorflow dtype, default = tf.float64
        dtype for the ar_func and h_signed.

    Returns
    -------
    params: list
        List of parameters as tensorflow variables.
    h_signed : dtype
        :math:`\log(h)` where :math:`h` is the BEAR concentration parameter.
    ar_func : function
        The autoregressive function.
    """
    ar_func, ar_func_params = _make_ref_ar_func(lag, alphabet_size, make_ar_func, af_kwargs, dtype)
    h_signed = tf.Variable(0, dtype=dtype, name='h_signed')
    params = ([h_signed] + ar_func_params)
    return params, h_signed, ar_func

def change_scope_params(lag, alphabet_size, make_ar_func,
                        af_kwargs, params, dtype=tf.float64):
    """Redefine and get parameters of BEAR or AR distribution in given scope.

    Used to to get unmirrored variables after training on multiple GPUs in parallel
    or to unpack a list of params into h and the autoregressive function.

    Parameters
    ----------
    lag : int
    alphabet_size : int
    make_ar_func : function
        Takes lag, alphabet_size, af_kwargs, dtype and returns a ar_func.
    af_kwargs : dict
        Keyword arguments for particular ar_func. For example, filter_width.
    params: list
        List of parameters as tensorflow variables.
    dtype : dtype, default = tf.float64
        dtype for the ar_func and h.

    Returns
    -------
    params : list
        List of parameters as tensorflow variables.
    h_signed : dtype
        :math:`\log(h)` where :math:`h` is the BEAR parameter.
    ar_func : function
        The autoregressive function.
    """
    pos_in_params = 0
    h_nmir = tf.Variable(0, dtype=dtype, name='h_signed')
    h_nmir.assign(params[pos_in_params])
    pos_in_params += 1
    ar_func_nmir, ar_func_params_nmir = _make_ref_ar_func(lag, alphabet_size, make_ar_func, af_kwargs, dtype)
    for param in ar_func_params_nmir:
        param.assign(params[pos_in_params])
        pos_in_params += 1
    params_nmir = ([h_nmir] + ar_func_params_nmir)
    return params_nmir, h_nmir, ar_func_nmir

def _train_step(batch, num_kmers, h_signed, ar_func,
                params, acc_grads, train_ar):
    """Add gradient of unbiased estimate of loss to accumulated gradients.

    Parameters
    ----------
    batch : list of two tensors of the same dtype.
        The first element is a one hot encoding of kmers of shape
        [kmer_batch_size, lag, alphabet_size+1] and the second is the transition
        counts of size [kmer_batch_size, alphabet_size+1].
    num_kmers : int
        Total number of kmers seen in data. Used to normalize estimate of loss.
    h_signed : dtype
        :math:`\log(h)` where :math:`h` is the BEAR concentration parameter.
    ar_func : function
        A function that takes a tensor of shape [A1, ..., An, lag, alphabet_size+1]
        of dtype and returns a tensor of shape [A1, ..., An, alphabet_size+1] of dtype
        of transition probabilities for each kmer. The autoregressive function.
    params : list
        List of parameters as tensorflow variables.
    acc_grads : list
        List of accumulated gradients for parameters.
    train_ar : bool
        Whether to evaluate the likelihood using an AR (True) or BEAR (False) model.

    Returns
    -------
    loss : dtype
        An unbiased estimate of the log likelihood of the data.
    """
    with tf.GradientTape() as grad_tape:
        kmer_batch_size = tf.shape(batch[0])[0]
        kmer_seqs = batch[0]
        transition_counts = batch[1]
        kmer_total_counts = tf.math.reduce_sum(transition_counts, axis=-1)
        ref_counts = batch[2]

        if train_ar:
            post = _ar_kmer_counts(kmer_seqs, kmer_total_counts, ref_counts, ar_func)
        else:
            post = _bear_kmer_counts(kmer_seqs, kmer_total_counts, ref_counts,
                                     h=tf.math.exp(h_signed), ar_func=ar_func)
        log_likelihood = tf.reduce_sum(
            post.counts_log_prob(transition_counts))

        elbo = (num_kmers / kmer_batch_size) * log_likelihood
        loss = -elbo

    gradients = grad_tape.gradient(loss, params)
    for tv, grad in zip(acc_grads, gradients):
        if grad is not None:
            tv.assign_add(grad)
    return loss

def train(data, num_kmers, epochs, ds_loc, ds_loc_ref, alphabet, lag, make_ar_func, af_kwargs,
          learning_rate, optimizer_name, train_ar, acc_steps=1,
          params_restart=None, writer=None, dtype=tf.float64):
    """Train a BEAR or AR model based on reference transition counts using all available GPUs in parallel.

    Parameters
    ----------
    data : tensorflow data object
        Load sequence data using tools in dataloader.py. Minibatch before passing.
    num_kmers : int
        Total number of kmers seen in the data. Used to normalize estimate of loss.
    epochs : int
    ds_loc : int
        Column in count data that corresponds with the training data.
    ds_loc_ref : int
        Column in count data that corresponds with the reference data.
    alphabet : str
        One of 'dna', 'rna', 'prot'.
    lag : int
    make_ar_func : function
        Takes lag, alphabet_size, af_kwargs, dtype and returns an ar_func.  See ar_funcs submodule.
    af_kwargs : dict
        Keyword arguments for particular ar_func. For example, filter_width.
    learning_rate : float
    optimizer_name : str
        For example 'Adam'.
    train_ar : bool
        Whether to train an AR (True) or BEAR (False) model.
    writer : tensorboard writer object, default = None
    acc_steps : int, default = 1
        Number of steps to accumulate gradients over.
    params_restart : list of tensorflow variables, default = None
        Pass the parameter list from a previous run to restart training.
    dtype : dtype, default=tf.float64

    Returns
    -------
    params: list
        List of parameters as tensorflow variables.
    h_signed : dtype
        :math:`\log(h)` where :math:`h` is the BEAR concentration parameter.
    ar_func : function
        The autoregressive function.
    """
    alphabet_size = len(core.alphabets_tf[alphabet]) - 1

    strategy = tf.distribute.MirroredStrategy()
    with strategy.scope():
        # Define parameters in parallel GPU usage scope.
        if params_restart is None:
            params, h_signed, ar_func = _create_params(
                lag, alphabet_size, make_ar_func,
                af_kwargs, dtype=dtype)
        else:
            params, h_signed, ar_func = change_scope_params(
                lag, alphabet_size, make_ar_func,
                af_kwargs, params_restart, dtype=dtype)

        # Set up accumulated gradient variables.p;';
        acc_grads = [tf.Variable(tf.zeros_like(tv), trainable=False) for tv in params]
        for tv in acc_grads:
            tv.assign(tf.zeros_like(tv))

        # Define optimizer in parallel GPU usage scope.
        optimizer = getattr(tf.keras.optimizers, optimizer_name)(
                                    learning_rate=learning_rate)

    # One hot encode kmers and get appropriate column from training data.
    not_stop = (1-tf.eye(alphabet_size+1, dtype=dtype)[-1])

    def map_(kmers, counts):
        return (core.tf_one_hot(kmers, alphabet),
                tf.gather(counts, ds_loc, axis=1),
                (tf.gather(counts, ds_loc_ref, axis=1) + epsilon) * not_stop)
    data = data.map(
        map_, num_parallel_calls=tf.data.experimental.AUTOTUNE).prefetch(10)
    # Get data iterable for use with GPU parallelization.
    data_iter = iter(strategy.experimental_distribute_dataset(data))

    # Define functions for updating parameters and accumulating gradients
    # with GPU parallelization.
    def add_grads(params, acc_grads, optimizer):
        optimizer.apply_gradients(zip(acc_grads, params))

    @tf.function
    def dist_add_grads(params, acc_grads, optimizer):
        strategy.run(add_grads, args=(params, acc_grads, optimizer))

    @tf.function
    def dist_train_step(batch, num_kmers, h_signed, ar_func,
                        params, acc_grads, train_ar):
        losses = strategy.run(_train_step, args=(
            batch, num_kmers, h_signed, ar_func,
            params, acc_grads, train_ar))
        return strategy.reduce(tf.distribute.ReduceOp.SUM, losses,
                               axis=None)

    # Training loop
    loss = 0.
    step = 1
    for batch in data_iter:
        # Accumulate gradients and update loss.
        loss += dist_train_step(batch, num_kmers, h_signed, ar_func,
                                params, acc_grads, train_ar)

        if step % acc_steps == 0:
            # Record loss
            if writer is not None:
                with writer.as_default():
                    tf.summary.scalar('elbo', - loss / acc_steps, step=step)
            # Update gradients
            dist_add_grads(params, acc_grads, optimizer)
            # Reset accumulated gradients and cumulative loss to zero.
            for tv in acc_grads:
                tv.assign(tf.zeros_like(tv))
            loss = 0

        step += 1

    # Remove the distributed scope from the parameters.
    params, h_signed, ar_func = change_scope_params(
        lag, alphabet_size, make_ar_func,
        af_kwargs, params, dtype=dtype)
    return params, h_signed, ar_func

@tf.function
def evaluation(data, ds_loc_train, ds_loc_test, ds_loc_ref,
               alphabet, h, ar_func, van_reg, dtype=tf.float64):
    """Evaluate a trained BEAR, AR or vanilla BEAR model.

    Parameters
    ----------
    data : tensorflow data object
        Load sequence data using tools in dataloader.py. Minibatch before passing.
    ds_loc_train : int
        Column in count data that corresponds with the training data.
    ds_loc_test : int
        Column in count data that corresponds with the testing data.
    alphabet : str
        One of 'dna', 'rna', 'prot'.
    h : dtype
        The :math:`h` parameter from the BEAR model.
    ar_func : function
        A function that takes a tensor of shape [A1, ..., An, lag, alphabet_size+1]
        of dtype and another of shape [A1, ..., An, alphabet_size+1] and returns a
        tensor of shape [A1, ..., An, alphabet_size+1] of dtype of transition probabilities
        for each kmer. The autoregressive function.
    van_reg : float
        Prior on vanilla BEAR model (Dirichlet concentration parameter).
    dtype : dtype, default=tf.float64

    Returns
    -------
    log_likelihood_ear, log_likelihood_arm, log_likelihood_van : float
        Total log likelihood of the data with the model evaluated as a BEAR,
        AR or vanilla BEAR model.
    perplexity_ear, perplexity_arm, perplexity_van : float
        Perplexity of the data with the model evaluated as a BEAR,
        AR or vanilla BEAR model.
    accuracy_ear, accuracy_arm, accuracy_van : float
        Accuracy of the data with the model evaluated as a BEAR,
        AR or vanilla BEAR model. Ties in maximum model probability are resolved randomly.
    """
    alphabet_size = len(core.alphabets_tf[alphabet]) - 1

    # Stops are removed from the reference counts as they are unlikely to be representative of
    # stops in read data.
    not_stop = (1-tf.eye(alphabet_size+1, dtype=dtype)[-1])

    def map_(kmers, counts):
        return (core.tf_one_hot(kmers, alphabet),
                tf.gather(counts, ds_loc_train, axis=1),
                tf.gather(counts, ds_loc_test, axis=1),
                (tf.gather(counts, ds_loc_ref, axis=1) + epsilon) * not_stop)
    data_iter = iter(data.map(map_))

    log_likelihood_ear = tf.constant(0., dtype=dtype)
    log_likelihood_arm = tf.constant(0., dtype=dtype)
    log_likelihood_van = tf.constant(0., dtype=dtype)
    correct_ear = tf.constant(0., dtype=dtype)
    correct_arm = tf.constant(0., dtype=dtype)
    correct_van = tf.constant(0., dtype=dtype)
    total_len = tf.constant(0., dtype=dtype)

    for batch in data_iter:
        kmer_seqs = batch[0]
        transition_counts_train = batch[1]
        transition_counts_test = batch[2]
        kmer_total_counts_test = tf.math.reduce_sum(transition_counts_test, axis=-1)
        ref_counts = batch[3]
        # Get posteriors.
        post_ear = _bear_kmer_counts(kmer_seqs, kmer_total_counts_test, ref_counts,
                                     condition_trans_counts=transition_counts_train,
                                     h=h, ar_func=ar_func)
        post_arm = _ar_kmer_counts(kmer_seqs, kmer_total_counts_test, ref_counts, ar_func)
        post_van = _bear_kmer_counts(kmer_seqs, kmer_total_counts_test, ref_counts,
                                     condition_trans_counts=transition_counts_train + van_reg)
        # Get likelihoods.
        log_likelihood_ear += tf.reduce_sum(
            post_ear.counts_log_prob(transition_counts_test))
        log_likelihood_arm += tf.reduce_sum(
            post_arm.counts_log_prob(transition_counts_test))
        log_likelihood_van += tf.reduce_sum(
            post_van.counts_log_prob(transition_counts_test))
        # Get most likely transition and accuracy.
        ml_ear = post_ear.ml_output()
        ml_arm = post_arm.ml_output()
        ml_van = post_van.ml_output()
        oh_ml_ear = tf.cast(tf.math.equal(ml_ear[..., None],
                                          tf.range(alphabet_size+1, dtype=dtype)),
                            dtype=dtype)
        oh_ml_arm = tf.cast(tf.math.equal(ml_arm[..., None],
                                          tf.range(alphabet_size+1, dtype=dtype)),
                            dtype=dtype)
        oh_ml_van = tf.cast(tf.math.equal(ml_van[..., None],
                                          tf.range(alphabet_size+1, dtype=dtype)),
                            dtype=dtype)
        correct_ear += tf.math.reduce_sum(transition_counts_test*oh_ml_ear)
        correct_arm += tf.math.reduce_sum(transition_counts_test*oh_ml_arm)
        correct_van += tf.math.reduce_sum(transition_counts_test*oh_ml_van)
        # Sum total number of transitions.
        total_len += tf.math.reduce_sum(transition_counts_test)
    return (log_likelihood_ear, log_likelihood_arm, log_likelihood_van,
            tf.exp(-log_likelihood_ear/total_len),
            tf.exp(-log_likelihood_arm/total_len),
            tf.exp(-log_likelihood_van/total_len),
            correct_ear/total_len, correct_arm/total_len, correct_van/total_len)

max_stars_repo_name: AlanNawzadAmin/BEAR

max_issues_repo_name: AlanNawzadAmin/BEAR

max_forks_repo_name: AlanNawzadAmin/BEAR


text:
import glob
import math
import sqlite3
import sys
from itertools import product

import logzero
import pandas as pd
from logzero import logger
from statsmodels.tsa.seasonal import seasonal_decompose

sys.path.append("../source")
import queries

# Connect to logzero log file
log_path = "logs/"
log_file = "dashboard_app.log"
logzero.logfile(
    log_path + log_file,
    maxBytes=1e5,
    backupCount=1,
    disableStderrLogger=True,
)
logger.info(f"ts_tools logger initialized")

def get_db_connection(db_path, db_filename):
    conn = sqlite3.connect(db_path + db_filename)
    logger.info(f"Connection made: {conn}")
    return conn

def get_db_zipcodes(conn):
    cursor = conn.cursor()
    cursor.execute(queries.select_distinct_zips)
    zipcodes = cursor.fetchall()
    zipcodes = [z[0] for z in zipcodes]
    logger.info(f"Distinct zip codes: {zipcodes}")
    return zipcodes

def get_column_names(conn, table_name):
    cursor = conn.cursor()
    cursor.execute(queries.select_column_names, {"table_name": table_name})
    names = cursor.fetchall()
    names = [name[0] for name in names]
    logger.info(f"Column Names: {names}")
    return names

def get_locale_data(conn, zipcode):
    """retrieve locale data from geo_zipcodes table of the database"""
    cursor = conn.cursor()
    cursor.execute(queries.select_locale_data, {"zipcode": zipcode})
    query_data = cursor.fetchone()
    locale_data = [qd for qd in query_data]
    logger.info(f"Locale data: {locale_data}")
    return locale_data

def get_db_files(db_path="./"):
    """retrieve database file names from the file system"""
    db_files = [
        file.split("/")[-1]
        for file in glob.glob(db_path + "*.db")
        if file.split("/")[-1] != "geo_zipcodes.db"
    ]

    return tuple(sorted(db_files))

def get_irr_data(conn, zipcode):
    """retrieve irradiance data from database"""
    params = {"zipcode": zipcode}

    df = pd.read_sql(
        queries.select_nsr_rows,
        conn,
        params=params,
        index_col="date_time",
        parse_dates=["date_time"],
    )

    df.sort_index(axis=0, inplace=True)

    return df

def get_plots_layout(num_columns=1, num_items=1):
    """row, column dimension calculation"""
    return {"rows": (math.ceil(num_items / num_columns)), "columns": num_columns}

def get_data_decomps(df, period=12):
    """data decomposition"""
    decomps = {}
    cols = df.columns.tolist()

    for col in cols:
        decomps.update({col: seasonal_decompose(df[col], model="additive", period=period)})

    return decomps

def get_train_test(df, test_len_yrs=1):
    """produce a time-series friendly train/test split"""

    months = 12
    total_len = len(df)
    test_len = months * test_len_yrs
    train_len = total_len - test_len

    train = df.iloc[:train_len]  # [columns[forecast_on_idx]]
    test = df.iloc[train_len:]  # [columns[forecast_on_idx]]

    return train, val

reason: from statsmodels


text:
#!/usr/bin/env python

from __future__ import print_function

import numpy as np
import h5py
from matplotlib import pyplot as plt
from argparse import ArgumentParser

def read_crops(args):
    if args.verbose:
        print('reading from file: ', args.crop_file)
    if args.mode == 'npy':
        if args.verbose:
            print('mode: npy')
        data = np.load(args.crop_file)
        crops = data['crops']
        return crops, None
    elif args.mode == 'hdf5':
        if args.verbose:
            print('mode: hdf5')
        h = h5py.File(args.crop_file, 'r')
        # if present this will return the pids - otherwise None
        if 'pids' in h.keys():
            pids = h['pids'][:]
        else:
            pids = None
        return h['crops'], pids
    else:
        raise NotImplemented('mode unknown' + args.output_mode)

def main(args):
    print('loading crops from: ', args.crop_file)
    crops, pids = read_crops(args)
    print('crops: ', crops.shape)
    if pids is not None:
        print('read pids...')

    indices = range(len(crops))
    if args.index:
        if args.index < 0 or args.index >= len(crops):
            print('invalid index - valid are 0..', len(crops))
            sys.exit(1)
        indices = [args.index]
    if args.random:
        np.random.shuffle(indices)
    for iii in indices:
        print('example idx: ', iii)
        c = crops[iii,:]
        t = c.transpose((2,1,0))
        o = args.sub
        if o > 0:
            plt.imshow(t[o:-o,o:-o,:])
        else:
            plt.imshow(t)

        title = 'example idx: ' + str(iii)
        if pids is not None:
            title += 'pid: ' + str(pids[iii])
        plt.title(title)
        plt.show()

if __name__ == '__main__':
    parser = ArgumentParser(description=__doc__)
    parser.add_argument('--verbose', '-v', action='store_true')
    parser.add_argument('--random', '-r', action='store_true', help='random order')
    parser.add_argument('--index', '-i', type=int, help='only visualize this example')
    parser.add_argument('--sub', type=int, default=0,
                        help='only show a centered subcrop (e.g. 8,16,32)')
    parser.add_argument('--mode', '-m',
                        choices=('npy', 'hdf5'),
                        default='hdf5',
                        help="data format of the prepared datasets")
    parser.add_argument('crop_file')
    args = parser.parse_args()

    main(args)