text:
# Copyright 2018 The TensorFlow 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.
# ==============================================================================
"""Benchmarks for `tf.data.Dataset.batch()`."""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import time

import numpy as np

from tensorflow.python.client import session
from tensorflow.python.data.ops import dataset_ops
from tensorflow.python.framework import dtypes
from tensorflow.python.framework import sparse_tensor
from tensorflow.python.ops import array_ops
from tensorflow.python.platform import test

# TODO(b/119837791): Add eager benchmarks.
class BatchBenchmark(test.Benchmark):
  """Benchmarks for `tf.data.Dataset.batch()`."""

  def benchmarkBatchSparse(self):
    non_zeros_per_row_values = [0, 1, 5, 10, 100]
    batch_size_values = [1, 32, 64, 128, 1024]

    sparse_placeholder = array_ops.sparse_placeholder(dtype=dtypes.int64)
    batch_size_placeholder = array_ops.placeholder(dtype=dtypes.int64, shape=[])

    dataset = dataset_ops.Dataset.from_tensors(sparse_placeholder).repeat(
        ).batch(batch_size_placeholder)
    iterator = dataset_ops.make_initializable_iterator(dataset)
    next_element = iterator.get_next()

    for non_zeros_per_row in non_zeros_per_row_values:

      sparse_value = sparse_tensor.SparseTensorValue(
          indices=np.arange(non_zeros_per_row, dtype=np.int64)[:, np.newaxis],
          values=np.arange(non_zeros_per_row, dtype=np.int64),
          dense_shape=[1000])

      for batch_size in batch_size_values:

        with session.Session() as sess:
          sess.run(iterator.initializer, feed_dict={
              sparse_placeholder: sparse_value,
              batch_size_placeholder: batch_size})
          # Run five steps to warm up the session caches before taking the
          # first measurement.
          for _ in range(5):
            sess.run(next_element.indices.op)
          deltas = []
          for _ in range(100):
            start = time.time()
            for _ in range(100):
              sess.run(next_element.indices.op)
            end = time.time()
            deltas.append(end - start)

        median_wall_time = np.median(deltas) / 100.0

        print("Batch sparse dataset non-zeros per row: %d batch_size: %d "
              "wall time: %f"
              % (non_zeros_per_row, batch_size, median_wall_time))
        self.report_benchmark(
            iters=10000, wall_time=median_wall_time,
            name="batch_sparse_dataset_nnz_%d_batch_size_%d" % (
                non_zeros_per_row, batch_size))

if __name__ == "__main__":
  test.main()

max_stars_repo_name: owennewo/kfserving

max_issues_repo_name: owennewo/kfserving

max_forks_repo_name: owennewo/kfserving


text:
# -*- coding: utf-8 -*-
"""
@author: david
"""

import numpy as np
from src.clean_data import DataCleaner as dc
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler as ss
from src.generate_features import GenerateFeatures as gf

from models.model_longitud import ModelLongitud as MLongitud
from models.model_custom_classifier import ModelCustomClassifier as MCClf
from sklearn.linear_model import LogisticRegression as LR
from sklearn.ensemble import RandomForestClassifier as RF

from src.model_evaluation import ModelEvaluation as ME

def main(file, regenerate=True):
    clean_file = 'data/clean_data.csv'
    dataset = dc.loadCleanData(data=file, clean_data=clean_file, regenerate=regenerate)
    contrasenyes = dataset.values[:,0]
    seguretat = dataset.values[:,1].astype(np.int64)
    pipeLongitud = Pipeline([("model", MLongitud())])
    clf_args = {"clf1_max_iter": 1e6, "clf1_solver": "liblinear",\
                "clf1_C": 10, "clf1_penalty": "l2",\
                "clf2_n_estimators": 20, "clf2_max_depth": None,\
                "clf2_criterion": "gini", "clf2_max_features": "log2"}
    pipeCClf = Pipeline([("extraction", gf([gf.mitja_cadena, gf.ponderacio_cadena,\
                                        gf.teMinuscules, gf.teMajuscules, gf.teXifres, gf.teEspecials],\
                          ferRecompte=False)), ("scaling", ss()),\
                         ("model", MCClf(classifier1 = LR, classifier2 = RF, **clf_args))])
    evaluacio = {'kFold': 10, 'class_weighted': True,\
                 'plot': ['confusion', 'percentage', 'AUC', 'ROC'],\
                'score': ['all', 'accuracy', 'class accuracy']}
    ME.evaluate(pipeLongitud, contrasenyes, seguretat, name="Model sobre la longitud", **evaluacio)
    ME.evaluate(pipeCClf, contrasenyes, seguretat, name="Model propi sense la longitud", **evaluacio)

if __name__ == "__main__":
    main('data/data.csv', regenerate=False)

ext: py

reason: import numpy


text:
"""
Licensed Materials - Property of IBM
Restricted Materials of IBM
20190891
© Copyright IBM Corp. 2021 All Rights Reserved.
"""
import logging
import joblib
import numpy as np
import time

from sklearn.linear_model import SGDClassifier, SGDRegressor
from sklearn.utils.validation import check_is_fitted
from sklearn.exceptions import NotFittedError

from ibmfl.util import fl_metrics
from ibmfl.util import config
from ibmfl.model.sklearn_fl_model import SklearnFLModel
from ibmfl.model.model_update import ModelUpdate
from ibmfl.exceptions import LocalTrainingException, \
    ModelInitializationException, ModelException

logger = logging.getLogger(__name__)

class SklearnSGDFLModel(SklearnFLModel):
    """
    Wrapper class for sklearn.linear_model.SGDClassifier and
    sklearn.linear_model.SGDRegressor.
    """

    def __init__(self, model_name, model_spec, sklearn_model=None, **kwargs):
        """
        Create a `SklearnSGDFLModel` instance from a
        sklearn.linear_model.SGDClassifier or a
        sklearn.linear_model.SGDRegressor model.
        If sklearn_model is provided, it will use it; otherwise it will take
        the model_spec to create the model.

        :param model_name: A name specifying the type of model, e.g., \
        linear_SVM
        :type model_name: `str`
        :param model_spec: A dictionary contains model specification
        :type model_spec: `dict`
        :param sklearn_model: Compiled sklearn model
        :type sklearn_model: `sklearn.linear_model`
        :param kwargs: A dictionary contains other parameter settings on \
         to initialize a sklearn.linear_model model.
        :type kwargs: `dict`
        """
        super().__init__(model_name, model_spec,
                         sklearn_model=sklearn_model,
                         **kwargs)
        self.model_type = 'Sklearn-linear-SGD'
        if sklearn_model:
            if not issubclass(type(sklearn_model), (SGDClassifier,
                                                    SGDRegressor)):
                raise ValueError('Compiled sklearn model needs to be provided'
                                 '(sklearn.linear_model). '
                                 'Type provided: ' + str(type(sklearn_model)))

            self.model = sklearn_model

        if type(self.model) is SGDClassifier:
            self.is_classification = True
        elif type(self.model) is SGDRegressor:
            self.is_classification = False

        # checking if classes_ is provided for classification problem
        if isinstance(self.model, SGDClassifier) and \
                hasattr(self.model, "classes_"):
            # check the type
            if isinstance(self.model.classes_, np.ndarray) or \
                    self.model.classes_ is None:
                self.classes = self.model.classes_
            else:
                raise ValueError(
                    "Provided SGDClassifier model has wrong type of "
                    "`classes_` attribute. "
                    "`classes_` should be of type `numpy.ndarray`. "
                    "Instead it is of type " + str(type(self.model.classes_)))
        else:
            self.classes = None

    def fit_model(self, train_data, fit_params=None, **kwargs):
        """
        Fits current model with provided training data.

        :param train_data: Training data a tuple \
        given in the form (x_train, y_train).
        :type train_data: `np.ndarray`
        :param fit_params: (Optional) Dictionary with hyperparameters that \
        will be used to call sklearn.linear_model fit function. \
        Provided hyperparameter should only contains parameters that \
        match sklearn expected values, e.g., `learning_rate`, which provides \
        the learning rate schedule. \
        If no `learning_rate` or `max_iter` is provided, a default value will \
        be used ( `optimal` and `1`, respectively).
        :return: None
        """
        # Default values
        max_iter = 1
        warm_start = True

        # Extract x_train and y_train, by default,
        # label is stored in the last column
        x = train_data[0]
        y = train_data[1]

        hyperparams = fit_params.get('hyperparams', {}) or {} if fit_params else {}
        sample_weight = fit_params.get('sample_weight', None) if fit_params else {}

        local_hp = hyperparams.get('local', {}) or {}
        training_hp = local_hp.get('training', {}) or {}

        try:

            if 'max_iter' not in training_hp:
                training_hp['max_iter'] = max_iter
                logger.info('Using default max_iter: ' + str(max_iter))

            # set warm_start to True
            training_hp['warm_start'] = warm_start
            logger.info('Set warm_start as ' + str(warm_start))

            for key, val in training_hp.items():
                self.model.set_params(**{key: val})
        except Exception as e:
            logger.exception(str(e))
            raise LocalTrainingException(
                'Error occurred while setting up model parameters')

        try:
            # Compute the classes based on labels if classification problem
            if isinstance(self.model, SGDClassifier):
                if self.classes is None:
                    logger.warning(
                        "Obtaining class labels based on local dataset. "
                        "This may cause failures during aggregation "
                        "when parties have distinctive class labels. ")
                    self.classes = self.get_classes(labels=y)
                # sklearn `partial_fit` uses max_iter = 1,
                # manually start the local training cycles
                for iter in range(training_hp['max_iter']):
                    logger.info("Local training epoch " + str(iter+1) + ":")
                    self.model.partial_fit(x, y,
                                           classes=self.classes,
                                           sample_weight=sample_weight)
            elif isinstance(self.model, SGDRegressor):
                for iter in range(training_hp['max_iter']):
                    logger.info("Local training epoch " + str(iter + 1) + ":")
                    self.model.partial_fit(x, y, sample_weight=sample_weight)
            else:
                raise NotImplementedError
        except Exception as e:
            logger.info(str(e))
            raise LocalTrainingException(
                'Error occurred while performing model.fit'
            )

    def update_model(self, model_update):
        """
        Update sklearn model with provided model_update, where model_update
        should contains `coef_` and `intercept_` having the same dimension
        as expected by the sklearn.linear_model.
        `coef_` : np.ndarray, shape (1, n_features) if n_classes == 2
        else (n_classes, n_features)
        `intercept_` : np.ndarray, shape (1,) if n_classes == 2 else (n_classes,)

        :param model_update: `ModelUpdate` object that contains the coef_ and \
        the intercept vectors that will be used to update the model.
        :type model_update: `ModelUpdate`
        :return: None
        """
        if isinstance(model_update, ModelUpdate):
            weights = model_update.get('weights')

            if isinstance(self.model, SGDClassifier):
                coef = np.array(weights)[:, :-1]
                intercept = np.array(weights)[:, -1]
            elif isinstance(self.model, SGDRegressor):
                coef = np.array(weights)[:-1]
                intercept = np.array(weights)[-1].reshape(1,)
            else:
                raise LocalTrainingException(
                    "Expecting scitkit-learn model of "
                    "type either "
                    "sklearn.linear_model.SGDClassifier "
                    "or sklearn.linear_model.SGDRegressor."
                    "Instead provided model is of type "
                    + str(type(self.model)))
            try:
                self.model.coef_ = coef
                self.model.intercept_ = intercept
            except Exception as e:
                raise LocalTrainingException('Error occurred during '
                                             'updating the model weights. ' +
                                             str(e))
        else:
            raise LocalTrainingException('Provided model_update should be of '
                                         'type ModelUpdate. '
                                         'Instead they are: ' +
                                         str(type(model_update)))

    def get_model_update(self):
        """
        Generates a `ModelUpdate` object that will be sent to other entities.

        :return: ModelUpdate
        :rtype: `ModelUpdate`
        """

        coef = self.model.coef_
        intercept = self.model.intercept_
        if isinstance(self.model, SGDClassifier):
            n_classes = len(self.get_classes())
            # intercept is of shape (1,) if n_classes == 2 else (n_classes,)
            if n_classes == 2:
                intercept = np.reshape(intercept, [1, 1])
            else:
                intercept = np.reshape(intercept, [n_classes, 1])
            w = np.append(coef, intercept, axis=1)
        elif isinstance(self.model, SGDRegressor):
            w = np.append(coef, intercept)
        else:
            raise LocalTrainingException("Expecting scitkit-learn model of "
                                         "type either "
                                         "sklearn.linear_model.SGDClassifier "
                                         "or sklearn.linear_model.SGDRegressor."
                                         "Instead provided model is of type "
                                         + str(type(self.model)))

        return ModelUpdate(weights=w.tolist(),
                           coef=self.model.coef_,
                           intercept=self.model.intercept_)

    def evaluate(self, test_dataset, **kwargs):
        """
        Evaluates the model given testing data.
        :param test_dataset: Testing data, a tuple given in the form \
        (x_test, test) or a datagenerator of of type `keras.utils.Sequence`,
        `keras.preprocessing.image.ImageDataGenerator`
        :type test_dataset: `np.ndarray`

        :param kwargs: Dictionary of metrics available for the model
        :type kwargs: `dict`
        """

        if type(test_dataset) is tuple:
            x_test = test_dataset[0]
            y_test = test_dataset[1]

            return self.evaluate_model(x_test, y_test)

        else:
            raise ModelException("Invalid test dataset!")

    def evaluate_model(self, x, y, **kwargs):
        """
        Evaluates the model given test data x and the corresponding labels y.

        :param x: Samples with shape as expected by the model.
        :type x: `np.ndarray`
        :param y: Corresponding true labels to x
        :type y: `np.ndarray`
        :param kwargs: Optional sample weights accepted by model.score
        :return: score, mean accuracy on the given test data and labels
        :rtype: `dict`
        """
        acc = {}
        acc['score'] = self.model.score(x, y, **kwargs)

        y_pred = self.predict(x, **kwargs)

        if self.is_classification:
            acc['acc'] = acc['score']
            additional_metrics = fl_metrics.get_eval_metrics_for_classificaton(
                y, y_pred)
        else:
            additional_metrics = fl_metrics.get_eval_metrics_for_regression(
                y, y_pred)

        dict_metrics = {**acc, **additional_metrics}

        return dict_metrics

    def predict_proba(self, x):
        """
        Perform prediction for the given input.

        :param x: Samples with shape as expected by the model.
        :type x: `np.ndarray`
        :return: Array of predictions
        :rtype: `np.ndarray`
        """
        return self.model.predict_proba(x)

    def get_classes(self, labels=None):
        """
        Returns an array of shape (n_classes,). If self.classes is not None,
        return self.classes, else obtains the array based on provided labels.

        :param labels: Provided class labels to obtain the array of classes.
        :type labels: `numpy.ndarray`
        :return: An array of shape `(n_classes,)`.
        :rtype: `numpy.ndarray`
        """
        if isinstance(self.model, SGDRegressor):
            raise ValueError(
                "SGDRegressor is used for regression problems and "
                "does not have `classes_`.")
        if self.classes is not None:
            return self.classes
        elif hasattr(self.model, "classes_"):
            return self.model.classes_
        elif labels is not None:
            return np.unique(labels)
        else:
            raise NotFittedError(
                "The scikit-learn model has not been initialized with "
                "`classes_` attribute, "
                "please either manually specify `classes_` attribute as "
                "an array of shape (n_classes,) or "
                "provide labels to obtain the array of classes. ")

    def save_model(self, filename=None, path=None):
        """
        Save a sklearn model to file in the format specific
        to the framework requirement.

        :param filename: Name of the file where to store the model.
        :type filename: `str`
        :param path: Path of the folder where to store the model. If no path \
        is specified, the model will be stored in the default data location of \
        the library `DATA_PATH`.
        :type path: `str`
        :return: filename
        """
        if filename is None:
            file = self.model_name if self.model_name else self.model_type
            filename = '{}_{}.pickle'.format(file, time.time())

        full_path = super().get_model_absolute_path(filename)

        # check if the classes_ attribute exists for SGDClassifier
        # this attribute is required for prediction and scoring
        if isinstance(self.model, SGDClassifier) and \
                not hasattr(self.model, "classes_"):
            logger.warning(
                "The classification model to be saved has no `classes_` "
                "attribute and cannot be used for prediction!")

        with open(full_path, 'wb') as f:
            joblib.dump(self.model, f)
        logger.info('Model saved in path: %s.', full_path)
        return filename

    @staticmethod
    def load_model_from_spec(model_spec):
        """
        Loads model from provided model_spec, where model_spec is a `dict`
        that contains the following items: model_spec['model_definition']
        contains the model definition as
        type sklearn.linear_model.SGDClassifier
        or sklearn.linear_model.SGDRegressor.

        :param model_spec: Model specification contains \
        a compiled sklearn model.
        :param model_spec: `dict`
        :return: model
        :rtype: `sklearn.cluster`
        """
        model = None
        try:
            if 'model_definition' in model_spec:
                model_file = model_spec['model_definition']
                model_absolute_path = config.get_absolute_path(model_file)

                with open(model_absolute_path, 'rb') as f:
                    model = joblib.load(f)

                if not issubclass(type(model), (SGDClassifier, SGDRegressor)):
                    raise ValueError('Provided compiled model in model_spec '
                                     'should be of type sklearn.linear_model.'
                                     'Instead it is: ' + str(type(model)))
        except Exception as ex:
            raise ModelInitializationException('Model specification was '
                                               'badly formed. '+ str(ex))
        return model

    def is_fitted(self):
        """
        Return a boolean value indicating if the model is fitted or not.

        :return: res
        :rtype: `bool`
        """
        try:
            check_is_fitted(self.model, ("coef_", "intercept_"))
        except NotFittedError as ex:
            logger.warning(
                "Model has no attribute `coef_` and `intercept_`, "
                "and hence is not fitted yet.")
            return False
        return True

max_stars_repo_name: imec-int/federated-learning-lib

max_issues_repo_name: imec-int/federated-learning-lib

max_forks_repo_name: imec-int/federated-learning-lib


text:
"""
Preprocessing scripts.

Written by Ed Oughton.

Winter 2020

"""
import configparser
import glob
import gzip
import json
import math
import os
import shutil
# import geoio
import numpy as np
import pandas as pd
import geopandas
import requests

import urllib.request
from shapely.geometry import MultiPolygon, Polygon, mapping, box, shape
from shapely.ops import cascaded_union
from fiona.crs import from_epsg
import rasterio
from rasterio.mask import mask
from rasterstats import zonal_stats
from geopy.distance import distance
from pathos.multiprocessing import Pool, cpu_count
from tqdm import tqdm

from pytal.utils import load_config #,create_folder_lookup,map_roads,line_length

CONFIG = configparser.ConfigParser()
CONFIG.read(os.path.join(os.path.dirname(__file__), 'script_config.ini'))
BASE_PATH = CONFIG['file_locations']['base_path']

DATA_RAW = os.path.join(BASE_PATH, 'raw')
DATA_INTERMEDIATE = os.path.join(BASE_PATH, 'intermediate')
DATA_PROCESSED = os.path.join(BASE_PATH, 'processed')

def process_country_shapes():
    """
    Created a set of global country shapes. Adds the single national boundary for
    each country to each country folder.

    """
    path_processed = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')

    if not os.path.exists(path_processed):

        print('Working on global_countries.shp')
        path_raw = os.path.join(DATA_RAW, 'gadm36_levels_shp', 'gadm36_0.shp')
        countries = geopandas.read_file(path_raw)

        print('Excluding Antarctica')
        countries = countries.loc[~countries['NAME_0'].isin(['Antarctica'])]

        print('Excluding small shapes')
        countries['geometry'] = countries.apply(exclude_small_shapes,axis=1)

        print('Simplifying geometries')
        countries['geometry'] = countries.simplify(tolerance = 0.005, preserve_topology=True) \
            .buffer(0.01).simplify(tolerance = 0.005, preserve_topology=True)

        print('Adding ISO country codes')
        glob_info_path = os.path.join(BASE_PATH, 'global_information.csv')
        load_glob_info = pd.read_csv(glob_info_path, encoding = "ISO-8859-1")

        countries = countries.merge(load_glob_info,left_on='GID_0', right_on='ISO_3digit')

        print('Writing global_countries.shp to file')
        countries.to_file(path_processed, driver='ESRI Shapefile')

    else:
        countries = geopandas.read_file(path_processed)

    for name in countries.GID_0.unique():

        path = os.path.join(DATA_INTERMEDIATE, name)

        if not os.path.exists(path):
            print('Creating directory {}'.format(path))
            os.makedirs(path)

        shape_path = os.path.join(path, 'national_outline.shp')

        if not os.path.exists(shape_path):

            print('working on {}'.format(name))

            single_country = countries[countries.GID_0 == name]

            single_country.to_file(shape_path)

    return print('All country shapes complete')

def process_regions(level):
    """
    Function for processing subnational regions.

    """
    filename = 'global_regions_{}.shp'.format(level)
    path_processed = os.path.join(DATA_INTERMEDIATE, filename)

    # if not os.path.exists(path_processed):

    print('Working on global_countries.shp')
    filename = 'gadm36_{}.shp'.format(level)
    path_regions = os.path.join(DATA_RAW, 'gadm36_levels_shp', filename)
    regions = geopandas.read_file(path_regions)

    print('Excluding Antarctica')
    regions = regions.loc[~regions['NAME_0'].isin(['Antarctica'])]
    regions = regions.loc[~regions['GID_1'].isin(['RUS.12_1'])] # Chukot
    regions = regions.loc[~regions['GID_1'].isin(['FJI.3_1'])] #Northern

    print('Excluding small shapes')
    regions['geometry'] = regions.apply(exclude_small_shapes,axis=1)

    print('Simplifying geometries')
    regions['geometry'] = regions.simplify(tolerance = 0.005, preserve_topology=True) \
        .buffer(0.01).simplify(tolerance = 0.005, preserve_topology=True)

    print('Adding ISO country codes')
    glob_info_path = os.path.join(BASE_PATH, 'global_information.csv')
    load_glob_info = pd.read_csv(glob_info_path, encoding = "ISO-8859-1")

    regions = regions.merge(load_glob_info,left_on='GID_0',right_on='ISO_3digit')
    regions.rename(columns={'coordinates':'coordinate'}, inplace=True)

    regions.reset_index(drop=True,inplace=True)

    print('Writing global_regions.shp to file')
    regions.to_file(path_processed, driver='ESRI Shapefile')

    print('Completed processing of regional shapes level {}'.format(level))

    # else:
    #     regions = geopandas.read_file(path_processed)

    return regions

def create_full_regional_layer_2():
    """

    """
    filename = 'global_regions_2.shp'
    path_processed = os.path.join(DATA_INTERMEDIATE, filename)

    if not os.path.exists(path_processed):

        #load regions and process
        regions_1 = process_regions(1)
        regions_2 = process_regions(2)

        #select desired regions
        regions_1 = regions_1.loc[regions_1['GID_0'].isin(['ESH', 'LBY', 'LSO'])]
        regions_1 = regions_1.drop_duplicates(subset='GID_1', keep="first")
        regions_2 = regions_2.loc[regions_2['gid_region'] == 2]
        regions_2 = regions_2.drop_duplicates(subset='GID_2', keep="first")

        #concatinate
        regions = pd.concat([regions_1, regions_2], sort=True)

        regions.reset_index(drop=True,inplace=True)

        print('Writing global_regions.shp to file')
        regions.to_file(path_processed, driver='ESRI Shapefile')

    return 'done'

def assemble_global_regional_layer():
    """
    Assemble a single layer of sub-national regional shapes.

    Set the desired level to be used in the 'global_information.csv' folder.

    Write each sub-national region to an individual shape in each nation's
    regions folder.

    """
    filename = 'global_regions.shp'
    path_processed = os.path.join(DATA_INTERMEDIATE, filename)

    if not os.path.exists(path_processed):

        #load regions and process
        regions_1 = process_regions(1)
        regions_2 = process_regions(2)

        #select desired regions
        regions_1 = regions_1.loc[regions_1['gid_region'] == 1]
        regions_2 = regions_2.loc[regions_2['gid_region'] == 2]

        #concatinate
        regions = pd.concat([regions_1, regions_2], sort=True)

        regions.reset_index(drop=True,inplace=True)

        print('Writing global_regions.shp to file')
        regions.to_file(path_processed, driver='ESRI Shapefile')

        print('Completed processing of regional shapes level')

        for name in regions.GID_0.unique():

            single_country = regions[regions.GID_0 == name]

            #try for level two regions
            try:
                for name_region in single_country.GID_2.unique():

                    single_region = single_country[single_country.GID_2 == name_region]

                    shape_path = os.path.join(path_processed, '{}.shp'.format(name_region))
                    single_region.to_file(shape_path)
            except:
                pass

            #try for level one regions
            try:
                for name_region in single_country.GID_1.unique():
                    single_region = single_country[single_country.GID_1 == name_region]

                    shape_path = os.path.join(path_processed, '{}.shp'.format(name_region))
                    single_region.to_file(shape_path)
            except:
                pass

    return print('Completed processing of regional shapes')

def find_country_list(continent_list):
    """

    """
    path_processed = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')
    countries = geopandas.read_file(path_processed)

    subset = countries.loc[countries['continent'].isin(continent_list)]

    country_list = []
    country_regional_levels = []

    for name in subset.GID_0.unique():

        if not name in ['LBY', 'ESH']:
            continue

        country_list.append(name)

        if name in ['LBY', 'ESH'] :
            regional_level =  1
        else:
            regional_level = 2

        country_regional_levels.append({
            'country': name,
            'regional_level': regional_level,
        })

    return country_list, country_regional_levels

def process_lowest_regions(country_list, country_regional_levels):
    """
    Function for processing subnational regions.

    """
    regions = []
    for country in country_list:

        level = get_region_level(country, country_regional_levels)

        filename = 'regions_{}_{}.shp'.format(level, country)
        folder = os.path.join(DATA_INTERMEDIATE, country, 'regions_lowest')
        path_processed = os.path.join(folder, filename)

        if not os.path.exists(path_processed):

            if not os.path.exists(folder):
                os.mkdir(folder)

            print('Working on regions')
            filename = 'gadm36_{}.shp'.format(level)
            path_regions = os.path.join(DATA_RAW, 'gadm36_levels_shp', filename)
            regions = geopandas.read_file(path_regions)

            path_country = os.path.join(DATA_INTERMEDIATE, country, 'national_outline.shp')
            countries = geopandas.read_file(path_country)

            for name in countries.GID_0.unique():

                if not name == country:
                    continue

                print('Working on {}'.format(name))
                regions = regions[regions.GID_0 == name]

                # print('Excluding small shapes')
                # regions['geometry'] = regions.apply(exclude_small_shapes,axis=1)

                # print('Simplifying geometries')
                # regions['geometry'] = regions.simplify(tolerance = 0.005, preserve_topology=True) \
                #     .buffer(0.01).simplify(tolerance = 0.005, preserve_topology=True)
                try:
                    print('Writing global_regions.shp to file')
                    regions.to_file(path_processed, driver='ESRI Shapefile')
                except:
                    pass
        print('Completed processing of regional shapes level {}'.format(level))

    return print('complete')

def segment_lowest_layer(country_list, country_regional_levels):
    """

    """

    path_countries = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')

    if os.path.exists(path_countries):
        countries = geopandas.read_file(path_countries)
    else:
        print('Must generate global_countries.shp first' )

    for country in countries.GID_0.unique():

        print('Working on segmenting lowest layer in {}'.format(country))

        if not country in country_list:
            continue

        level = get_region_level(country, country_regional_levels)

        filename = 'regions_{}_{}.shp'.format(level, country)
        path_regions = os.path.join(DATA_INTERMEDIATE, country, 'regions_lowest', filename)

        try:
            regions = geopandas.read_file(path_regions)

            level_name = 'GID_{}'.format(level)

            for name_region in regions[level_name]:

                single_region = regions[regions[level_name] == name_region]

                folder = os.path.join(DATA_INTERMEDIATE, country, 'regions_lowest')
                shape_path = os.path.join(folder, '{}.shp'.format(name_region))
                single_region.to_file(shape_path)
        except:
            pass

    return print('Completed processing of regional shapes')

def process_settlement_layer(country_list):
    """
    Clip the settlement layer to each country boundary and place in each country folder.

    """
    path_settlements = os.path.join(DATA_RAW,'settlement_layer','ppp_2020_1km_Aggregated.tif')

    settlements = rasterio.open(path_settlements, 'r+')
    settlements.nodata = 255
    settlements.crs = {"init": "epsg:4326"}

    path_countries = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')

    if os.path.exists(path_countries):
        countries = geopandas.read_file(path_countries)
    else:
        print('Must generate global_countries.shp first' )

    num = 0
    for name in countries.GID_0.unique():

        if not name in country_list:
            continue

        print('working on {}'.format(name))

        path_country = os.path.join(DATA_INTERMEDIATE, name)
        shape_path = os.path.join(path_country, 'settlements.tif')

        # if os.path.exists(shape_path):
        #     continue

        single_country = countries[countries.GID_0 == name]

        bbox = single_country.envelope
        geo = geopandas.GeoDataFrame()
        # print(bbox)
        geo = geopandas.GeoDataFrame({'geometry': bbox})
        # print(geo)
        num += 1
        coords = [json.loads(geo.to_json())['features'][0]['geometry']]

        #chop on coords
        out_img, out_transform = mask(settlements, coords, crop=True)

        # Copy the metadata
        out_meta = settlements.meta.copy()

        out_meta.update({"driver": "GTiff",
                        "height": out_img.shape[1],
                        "width": out_img.shape[2],
                        "transform": out_transform,
                        "crs": 'epsg:4326'})

        with rasterio.open(shape_path, "w", **out_meta) as dest:
                dest.write(out_img)

    print('Completed processing of settlement layer')

def process_night_lights(country_list):
    """
    Clip the nightlights layer and place in each country folder.

    """
    path_processed = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')

    filename = 'F182013.v4c_web.stable_lights.avg_vis.tif'
    path_night_lights = os.path.join(DATA_RAW, 'nightlights', '2013', filename)

    countries = geopandas.read_file(path_processed)

    for name in countries.GID_0.unique():

        if not name in country_list:
            continue

        print('working on {}'.format(name))

        path_country = os.path.join(DATA_INTERMEDIATE, name)

        existing_files = glob.glob(os.path.join(path_country, 'night_lights*'))
        for existing_file in existing_files:
            if os.path.exists(existing_file):
                os.remove(existing_file)

        single_country = countries[countries.GID_0 == name]

        bbox = single_country.envelope

        geo = geopandas.GeoDataFrame()
        geo = geopandas.GeoDataFrame({'geometry': bbox}, crs=from_epsg('4326'))

        coords = [json.loads(geo.to_json())['features'][0]['geometry']]

        night_lights = rasterio.open(path_night_lights, "r+")
        night_lights.nodata = 0

        out_img, out_transform = mask(night_lights, coords, crop=True)

        out_meta = night_lights.meta.copy()

        out_meta.update({"driver": "GTiff",
                        "height": out_img.shape[1],
                        "width": out_img.shape[2],
                        "transform": out_transform,
                        "crs": 'epsg:4326'})

        shape_path = os.path.join(path_country,'night_lights.tif')

        with rasterio.open(shape_path, "w", **out_meta) as dest:
                dest.write(out_img)

    return print('Completed processing of night lights layer')

def get_regional_data(country_list, country_regional_levels):
    """
    Extract luminosity values.

    """
    path_processed = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')

    countries = geopandas.read_file(path_processed)

    for name in countries.GID_0.unique():

        if not name in country_list:
            continue

        single_country = countries[countries.GID_0 == name]

        print('working on {}'.format(name))

        path_country = os.path.join(DATA_INTERMEDIATE, name, 'regions_lowest')

        if os.path.exists(os.path.join(path_country, '..', 'luminosity.csv')):
            os.remove(os.path.join(path_country, '..', 'luminosity.csv'))

        path_night_lights = os.path.join(DATA_INTERMEDIATE, name, 'night_lights.tif')
        path_settlements = os.path.join(DATA_INTERMEDIATE, name, 'settlements.tif')

        path_regions = glob.glob(os.path.join(path_country, '*.shp'))

        level = get_region_level(name, country_regional_levels)

        results = []

        for path_region in path_regions:

            region = geopandas.read_file(path_region)
            # gid_level = 'GID_{}'.format(level)
            # print(region[gid_level].values[0])
            try:
                # print('get night light values')
                with rasterio.open(path_night_lights) as src:
                    # src.nodata = np.nan
                    affine = src.transform
                    array = src.read(1)

                    #set missing data (-999) to 0
                    # array[array == 0] = 0
                    array[array <= 0] = 0

                    # kwargs.update({'nodata': 0})
                    #get luminosity values
                    luminosity_median = [
                        d['median']
                        for d in zonal_stats(region, array, stats=['median'], affine=affine, nodata=0)
                    ][0]
                    luminosity_summation = [
                        d['sum']
                        for d in zonal_stats(region, array, stats=['sum'], affine=affine, nodata=0)
                    ][0]

                # print('#get settlement values')
                with rasterio.open(path_settlements) as src:
                    # src.nodata = np.nan
                    affine = src.transform
                    array = src.read(1)

                    # array[array == None] = 0
                    array[array <= 0] = 0

                    #get luminosity values
                    population_summation = [
                        d['sum']
                        for d in zonal_stats(region, array, stats=['sum'], affine=affine, nodata=0)
                    ][0]

                #get geographic area
                region.crs = "epsg:4326"
                region = region.to_crs("epsg:3857")

                area_km2 = region['geometry'].area[0] / 10**6

                if luminosity_median == None:
                    luminosity_median = 0
                if luminosity_summation == None:
                    luminosity_summation = 0

                gid_level = 'GID_{}'.format(level)

                results.append({
                    'GID_0': single_country.GID_0.values[0],
                    gid_level: region[gid_level].values[0],
                    'median_luminosity': luminosity_median,
                    'sum_luminosity': luminosity_summation,
                    'mean_luminosity_km2': luminosity_summation / area_km2,
                    'population': population_summation,
                    'area_km2': area_km2,
                    'population_km2': population_summation / area_km2,
                })

            except:
                pass

        results_df = pd.DataFrame(results)

        results_df.to_csv(os.path.join(path_country, '..', 'regional_data.csv'), index=False)

    return print('Completed night lights data querying')

def get_region_level(name, country_regional_levels):

    for country_regional_level in country_regional_levels:
        if country_regional_level['country'] == name:
            level = country_regional_level['regional_level']

    return level

def exclude_small_shapes(x,regionalized=False):
    """
    This function will remove the small shapes of multipolygons. Will reduce the size
        of the file.

    Arguments:
        *x* : a geometry feature (Polygon) to simplify. Countries which are very large will
        see larger (unhabitated) islands being removed.

    Optional Arguments:
        *regionalized*  : Default is **False**. Set to **True** will use lower threshold
        settings (default: **False**).

    Returns:
        *MultiPolygon* : a shapely geometry MultiPolygon without tiny shapes.

    """
    # if its a single polygon, just return the polygon geometry
    if x.geometry.geom_type == 'Polygon':
        return x.geometry

    # if its a multipolygon, we start trying to simplify and remove shapes if its too big.
    elif x.geometry.geom_type == 'MultiPolygon':

        if regionalized == False:
            area1 = 0.1
            area2 = 250

        elif regionalized == True:
            area1 = 0.01
            area2 = 50

        # dont remove shapes if total area is already very small
        if x.geometry.area < area1:
            return x.geometry
        # remove bigger shapes if country is really big

        if x['GID_0'] in ['CHL','IDN']:
            threshold = 0.01
        elif x['GID_0'] in ['RUS','GRL','CAN','USA']:
            if regionalized == True:
                threshold = 0.01
            else:
                threshold = 0.01

        elif x.geometry.area > area2:
            threshold = 0.1
        else:
            threshold = 0.001

        # save remaining polygons as new multipolygon for the specific country
        new_geom = []
        for y in x.geometry:
            if y.area > threshold:
                new_geom.append(y)

        return MultiPolygon(new_geom)

def clean_coverage(x):

    # if its a single polygon, just return the polygon geometry
    if x.geometry.geom_type == 'Polygon':
        if x.geometry.area > 1e7:
            return x.geometry

    # if its a multipolygon, we start trying to simplify and remove shapes if its too big.
    elif x.geometry.geom_type == 'MultiPolygon':

        threshold = 1e7
        # hole_threshold_size = 1e7

        # save remaining polygons as new multipolygon for the specific country
        new_geom = []
        for y in x.geometry:

            if y.area > threshold:
                # print(y.boundary)
                # holes =[]

                # for hole in y.interiors:
                #     if Polygon(hole).area > hole_threshold_size:
                #         holes.append(hole)
                # # print(Polygon(y.exterior.coords,[holes]))
                # poly = Polygon(y.exterior.coords,[holes])

                new_geom.append(y)

        return MultiPolygon(new_geom)

def process_coverage_shapes():
    """
    Load in coverage maps, process and export for each country.

    """
    print('Working on coverage_2g.shp')
    path_raw = os.path.join(DATA_RAW, 'coverage_maps', 'african_coverage_2g_tiles.shp')
    coverage = geopandas.read_file(path_raw)

    print('Setting crs for tiles')
    coverage.crs = 'epsg:3857'

    print('Converting tiles to epsg 4326 (WGS84)')
    coverage = coverage.to_crs('epsg:4326')

    print('Importing global_countries.shp to get each country name')
    path_countries = os.path.join(DATA_INTERMEDIATE,'global_countries.shp')
    countries = geopandas.read_file(path_countries)

    for name in countries.GID_0.unique():

        # if not name == 'MWI':
        #     continue

        print('Working on {}'.format(name))
        path = os.path.join(DATA_INTERMEDIATE, name, 'national_outline.shp')
        national_outline = geopandas.read_file(path)

        print('Intersecting coverage tiles with country outline')
        intersection = geopandas.overlay(coverage, national_outline, how='intersection')

        if len(intersection) > 0:
            print('Exporting country coverage shape')
            output_path = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2g_unprocessed.shp')
            intersection.to_file(output_path, driver='ESRI Shapefile')
        else:
            print('Nothing to write for {}'.format(name))
            continue

    for name in countries.GID_0.unique():

        # if not name == 'MWI':
        #     continue

        path = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2g_unprocessed.shp')

        if os.path.exists(path):

            print('Working on {}'.format(name))

            print('Loading coverage shape')
            coverage = geopandas.read_file(path)

            print('Dissolving polygons')
            coverage['dissolve'] = 1
            coverage = coverage.dissolve(by='dissolve', aggfunc='sum')

            coverage = coverage.to_crs('epsg:3857')

            print('Excluding small shapes')
            coverage['geometry'] = coverage.apply(clean_coverage,axis=1)

            print('Removing empty and null geometries')
            coverage = coverage[~(coverage.is_empty)]
            coverage = coverage[coverage['geometry'].notnull()]

            if len(coverage) > 0:

                print('Simplifying geometries')
                coverage['geometry'] = coverage.simplify(tolerance = 0.005, \
                    preserve_topology=True).buffer(0.01).simplify(tolerance = 0.005, \
                    preserve_topology=True)

                print('Exporting country coverage shape')
                output_path = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2g.shp')
                coverage.to_file(output_path, driver='ESRI Shapefile')

                path1 = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2_unprocessed.shp')
                path2 = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2_unprocessed.shx')
                path3 = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2_unprocessed')
                path4 = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2_unprocessed.dbf')
                path5 = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2_unprocessed.cpg')
                path6 = os.path.join(DATA_INTERMEDIATE, name, 'coverage_2_unprocessed.prj')

                if os.path.exists(path1):
                    os.remove(path1)
                if os.path.exists(path2):
                    os.remove(path2)
                if os.path.exists(path3):
                    os.remove(path3)
                if os.path.exists(path4):
                    os.remove(path4)
                if os.path.exists(path5):
                    os.remove(path5)
                if os.path.exists(path5):
                    os.remove(path6)

            else:
                print('Nothing to write for {}'.format(name))
                continue

    print('Processed coverage shapes')

def load_extents(glob_interator):
    """
    Check the extent of each DEM tile, save to dict for future reference.

    For London DEM, this needs to be in EPSG: 4326.

    """
    extents = {}
    for tile_path in glob_interator:
        dataset = rasterio.open(tile_path)
        extents[tuple(dataset.bounds)] = tile_path

    return extents

def get_tile_path_for_point(extents, x, y):
    for (left, bottom, right, top), path in extents.items():
        if x >= left and x <= right and y <= top and y >= bottom:
            return path
    raise ValueError("No tile includes x {} y {}".format(x, y))

def query_data(regions, filepath_nl, filepath_lc, filepath_pop):
    """
    Query raster layer for each shape in regions.

    """
    shapes = []
    csv_data = []

    for region in tqdm(regions):

        geom = shape(region['geometry'])

        population = get_population(geom, filepath_pop)

        pop_density_km2, area_km2 = get_density(geom, population, 'epsg:4326', 'epsg:3857')

        shapes.append({
            'type': region['type'],
            'geometry': mapping(geom),
            # 'id': region['id'],
            'properties': {
                'population': population,
                'pop_density_km2': pop_density_km2,
                'area_km2': area_km2,
                'geotype': define_geotype(pop_density_km2),
                'GID_2': region['properties']['GID_2'],
                'GID_3': region['properties']['GID_3'],
            }
        })

        csv_data.append({
            'population': population,
            'pop_density_km2': pop_density_km2,
            'area_km2': area_km2,
            'geotype': define_geotype(pop_density_km2),
            'GID_2': region['properties']['GID_2'],
            'GID_3': region['properties']['GID_3'],
        })

    return shapes, csv_data

def get_population(geom, filepath_pop):
    """
    Get sum of population within geom.

    """
    population = zonal_stats(geom, filepath_pop, stats="sum", nodata=0)[0]['sum']

    try:
        if population >= 0:
            return population
        else:
            return 0
    except:
        return 0

def aggregate(shapes_lower, regions_upper):
    """
    Using zonal_stats on large areas can lead to a runtime warning due to an overflow.

    This function takes lower (smaller) regions and aggregates them into larger areas.

    """
    idx = index.Index()
    [idx.insert(0, shape(region['geometry']).bounds, region) for region in regions_upper]

    population_data = []

    for shape_lower in shapes_lower:
        for n in idx.intersection((shape(shape_lower['geometry']).bounds), objects=True):
            centroid = shape(shape_lower['geometry']).centroid
            upper_region_shape = shape(n.object['geometry'])
            if upper_region_shape.contains(centroid):
                population_data.append({
                    'population': shape_lower['properties']['population'],
                    'GID_2': n.object['properties']['GID_2'],
                    'GID_3': shape_lower['properties']['GID_3'],
                    })

    shapes = []
    csv_data = []

    for region in tqdm(regions_upper):

        region_id = region['properties']['GID_2']

        geom = shape(region['geometry'])

        population = aggregate_population(population_data, region_id)

        pop_density_km2, area_km2 = get_density(geom, population, 'epsg:4326', 'epsg:3857')

        shapes.append({
            'type': region['type'],
            'geometry': mapping(geom),
            'properties': {
                'population': population,
                'pop_density_km2': pop_density_km2,
                'area_km2': area_km2,
                'geotype': define_geotype(pop_density_km2),
                'GID_2': region['properties']['GID_2'],
            }
        })

        csv_data.append({
            'population': population,
            'pop_density_km2': pop_density_km2,
            'area_km2': area_km2,
            'geotype': define_geotype(pop_density_km2),
            'GID_2': region['properties']['GID_2'],
        })

    return shapes, csv_data

def aggregate_population(population_data, region_id):
    """
    Sum the population from lower level regions to
    upper level regions.

    """
    population = 0

    for item in population_data:
        if item['GID_2'] == region_id:
            population += item['population']

    return population

if __name__ == '__main__':

    # country_list = ['UGA', 'ETH', 'BGD', 'PER', 'MWI', 'ZAF']
    # country_regional_levels = [
    #     {'country': 'UGA', 'regional_level': 3},
    #     {'country': 'ETH', 'regional_level': 3},
    #     {'country': 'BGD', 'regional_level': 3},
    #     {'country': 'PER', 'regional_level': 3},
    #     {'country': 'MWI', 'regional_level': 3},
    #     {'country': 'ZAF', 'regional_level': 3},
    #     ]

    # ###create 'global_countries.shp' if not already processed
    # ###create each 'national_outline.shp' if not already processed
    # process_country_shapes()

    # # ###create 'global_regions.shp' if not already processed
    # # ###create each subnational region .shp if not already processed
    # assemble_global_regional_layer()

    create_full_regional_layer_2()

    country_list, country_regional_levels = find_country_list(['Africa', 'South America'])

    # # print('Processing lowest regions')
    # process_lowest_regions(country_list, country_regional_levels)

    # segment_lowest_layer(country_list, country_regional_levels)

    # print('Processing settlement layer')
    # process_settlement_layer(country_list)

    # print('Processing night lights')
    # process_night_lights(country_list)

    get_regional_data(country_list, country_regional_levels)

    # # # process_coverage_shapes()

max_stars_repo_name: tomalrussell/pytal

max_issues_repo_name: tomalrussell/pytal

max_forks_repo_name: tomalrussell/pytal


text:
import numpy as np
import pandas as pd
import pytest

from nilearn.interfaces.fmriprep.load_confounds_scrub import (
    _optimize_scrub, _get_outlier_cols, _extract_outlier_regressors
)

from pandas.testing import assert_frame_equal

@pytest.mark.parametrize(
    "original_motion_outliers_index,expected_optimal",
    [
        (
            [3, 11, 13, 80, 97],
            np.array([0, 1, 2, 3, 11, 12, 13, 80, 97, 98, 99]),
        ),
        (
            [11, 13, 16, 44, 50, 80],  # middle volumes
            np.array([11, 12, 13, 14, 15, 16, 44, 50, 80]),
        ),
        ([4], np.array([0, 1, 2, 3, 4])),  # head volumes
        ([96], np.array([96, 97, 98, 99])),  # tail volumes
        ([5], np.array([5])),
    ],
)  # no optimisation needed
def test_optimize_scrub(original_motion_outliers_index, expected_optimal):
    """Check the segment removal is acting correctly."""
    # simulated labels with 100 time frames and remove any segment under
    # 5 volumes
    optimised_index = _optimize_scrub(original_motion_outliers_index, 100, 5)
    assert np.array_equal(optimised_index, expected_optimal)

def test_get_outlier_cols():
    """Check the non-steady state columns are deteched."""
    col_names = ["confound_regressor"]
    non_steady_state = [f"non_steady_state_outlier{i:02d}" for i in range(3)]
    col_names += non_steady_state
    col_names = pd.Index(col_names)
    outlier_cols, confounds_cols = _get_outlier_cols(col_names)
    assert confounds_cols == ["confound_regressor"]
    assert outlier_cols == non_steady_state

def test_extract_outlier_regressors():
    """Check outlier regressors of different types."""
    # Create a fake confound dataframe
    n_scans = 50
    fake_confounds = pd.DataFrame(
        np.random.rand(n_scans, 1), columns=["confound_regressor"]
    )

    # scrubbed volume one-hot, overlap with non-steady-state
    idx_scrubbed = [2, 4, 34, 44]
    scrub_vol = pd.DataFrame(
        np.eye(n_scans)[:, idx_scrubbed],
        columns=[f"motion_outlier{i:02d}" for i in range(len(idx_scrubbed))],
    )
    # First three volumes are non-steady-state
    non_steady_vol = pd.DataFrame(
        np.eye(n_scans)[:, :3],
        columns=[f"non_steady_state_outlier{i:02d}" for i in range(3)],
    )

    # non-steady only
    non_steady_conf = pd.concat([fake_confounds, non_steady_vol], axis=1)
    sample_mask, confounds, outliers = _extract_outlier_regressors(
        non_steady_conf
    )
    assert np.array_equal(sample_mask, np.arange(n_scans)[3:]) is True
    assert_frame_equal(outliers, non_steady_vol)
    assert_frame_equal(confounds, fake_confounds)

    # scrub only
    srub_conf = pd.concat([fake_confounds, scrub_vol], axis=1)
    make_mask = np.delete(np.arange(n_scans), idx_scrubbed)
    sample_mask, confounds, outliers = _extract_outlier_regressors(
        srub_conf
    )
    assert np.array_equal(sample_mask, make_mask) is True
    assert_frame_equal(outliers, scrub_vol)
    assert_frame_equal(confounds, fake_confounds)

    # scrub and non-steady state
    all_conf = pd.concat([fake_confounds, non_steady_vol, scrub_vol], axis=1)
    make_mask = np.delete(np.arange(n_scans), idx_scrubbed)[2:]
    make_outliers = pd.concat([non_steady_vol, scrub_vol], axis=1)
    make_outliers = make_outliers.reindex(
        sorted(make_outliers.columns), axis=1
    )
    make_outliers = make_outliers.drop(columns="non_steady_state_outlier02")

    sample_mask, confounds, outliers = _extract_outlier_regressors(
        all_conf
    )
    assert len(sample_mask) == 44
    assert np.array_equal(sample_mask, make_mask) is True
    assert_frame_equal(outliers, make_outliers)
    assert_frame_equal(confounds, fake_confounds)

max_stars_repo_name: ctw/nilearn

max_issues_repo_name: ctw/nilearn

max_forks_repo_name: ctw/nilearn


text:
from numpy import random
from numpy.random.mtrand import rand, randint
import random
import numpy as np

target = [1,2,3,4]

def fitness(populacao):

    points = 0
    for h in range(len(target)):

        for j in range(len(populacao)):
            if populacao[j] == target[h]:
                points+=1

    points = points / len(target)

    return points

def selection(populacao):
    selecionados = []

    for i in range(len(populacao)):
        n = fitness(populacao[i]) * 100

        x = 0
        while(x<n):
            selecionados.append(populacao[i])
            x+=1
    return selecionados

def crossover(selecionados):

    a = randint(0,len(selecionados))
    b = randint(0,len(selecionados))

    parentA = selecionados[a]
    parentB = selecionados[b]

    child = []
    n = int(int(len(parentA)) /2)
    for i in range(n):
        child.append(parentA[i])

    for i in range(n):
        child.append(parentB[n+i])

    return child

def mutation(child,pop):

    MutationRate = 0.1
    r = random.uniform(0,1)

    rand = randint(0,len(child))

    if r < MutationRate:
        print("Child: ",child)

        print('Mutação ocorrida')
        child[rand] = np.random.randint(0,len(pop))

    print('Melhor individuo: ', child)
    return child

reason: import numpy,from numpy


text:
"""
Tyson Reimer
University of Manitoba
November 15th, 2021
"""

import os
import numpy as np
import matplotlib.pyplot as plt

from umbms import get_proj_path
from umbms.loadsave import load_pickle

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

__DIR = os.path.join(get_proj_path(), 'output/g3/')

# Define colors for plotting
das_col = [0, 0, 0]
dmas_col = [80, 80, 80]
gd_col = [160, 160, 160]
das_col = [ii / 255 for ii in das_col]
dmas_col = [ii / 255 for ii in dmas_col]
gd_col = [ii / 255 for ii in gd_col]

# Define additional colors for plotting
das_col2 = [124, 0, 0]
dmas_col2 = [255, 0, 0]
gd_col2 = [231, 167, 167]
das_col2 = [ii / 255 for ii in das_col2]
dmas_col2 = [ii / 255 for ii in dmas_col2]
gd_col2 = [ii / 255 for ii in gd_col2]

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

def get_aucs():
    ####### 2022-03-01 stuff....

    das_tps_adi = (das_adi_sens) / 100
    das_fps_adi = (100 - das_adi_spec) / 100
    dmas_tps_adi = (dmas_adi_sens) / 100
    dmas_fps_adi = (100 - dmas_adi_spec) / 100
    gd_tps_adi = (gd_adi_sens) / 100
    gd_fps_adi = (100 - gd_adi_spec) / 100

    das_tps_fib = (das_fib_sens) / 100
    das_fps_fib = (100 - das_fib_spec) / 100
    dmas_tps_fib = (dmas_fib_sens) / 100
    dmas_fps_fib = (100 - dmas_fib_spec) / 100
    gd_tps_fib = (gd_fib_sens) / 100
    gd_fps_fib = (100 - gd_fib_spec) / 100

    ##########

    das_adi_sort = np.argsort(das_fps_adi)
    das_tps_adi = das_tps_adi[das_adi_sort]
    das_fps_adi = das_fps_adi[das_adi_sort]
    dmas_adi_sort = np.argsort(dmas_fps_adi)
    dmas_tps_adi = dmas_tps_adi[dmas_adi_sort]
    dmas_fps_adi = dmas_fps_adi[dmas_adi_sort]
    gd_adi_sort = np.argsort(gd_fps_adi)
    gd_tps_adi = gd_tps_adi[gd_adi_sort]
    gd_fps_adi = gd_fps_adi[gd_adi_sort]

    das_fib_sort = np.argsort(das_fps_fib)
    das_tps_fib = das_tps_fib[das_fib_sort]
    das_fps_fib = das_fps_fib[das_fib_sort]
    dmas_fib_sort = np.argsort(dmas_fps_fib)
    dmas_tps_fib = dmas_tps_fib[dmas_fib_sort]
    dmas_fps_fib = dmas_fps_fib[dmas_fib_sort]
    gd_fib_sort = np.argsort(gd_fps_fib)
    gd_tps_fib = gd_tps_fib[gd_fib_sort]
    gd_fps_fib = gd_fps_fib[gd_fib_sort]

    das_adi_auc = 0
    dmas_adi_auc = 0
    gd_adi_auc = 0

    das_fib_auc = 0
    dmas_fib_auc = 0
    gd_fib_auc = 0

    for ii in range(len(das_tps_adi) - 1):

        das_d_fps = das_fps_adi[ii + 1] - das_fps_adi[ii]
        dmas_d_fps = dmas_fps_adi[ii + 1] - dmas_fps_adi[ii]
        gd_d_fps = gd_fps_adi[ii + 1] - gd_fps_adi[ii]

        if das_d_fps != 0:
            das_adi_auc += das_d_fps * das_tps_adi[ii + 1]

        if dmas_d_fps != 0:
            dmas_adi_auc += dmas_d_fps * dmas_tps_adi[ii + 1]

        if gd_d_fps != 0:
            gd_adi_auc += gd_d_fps * gd_tps_adi[ii + 1]

        ###

        das_fib_d_fps = das_fps_fib[ii + 1] - das_fps_fib[ii]
        dmas_fib_d_fps = dmas_fps_fib[ii + 1] - dmas_fps_fib[ii]
        gd_fib_d_fps = gd_fps_fib[ii + 1] - gd_fps_fib[ii]

        if das_fib_d_fps != 0:
            das_fib_auc += das_fib_d_fps * das_tps_fib[ii + 1]

        if dmas_fib_d_fps != 0:
            dmas_fib_auc += dmas_fib_d_fps * dmas_tps_fib[ii + 1]

        if gd_fib_d_fps != 0:
            gd_fib_auc += gd_fib_d_fps * gd_tps_fib[ii + 1]

    return (das_adi_auc, dmas_adi_auc, gd_adi_auc,
            das_fib_auc, dmas_fib_auc, gd_fib_auc)

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

if __name__ == "__main__":

    # Load the sensitivities and specificities when adipose-only scans
    # are used for reference subtraction
    (_, das_adi_sens, dmas_adi_sens, gd_adi_sens) = \
        load_pickle(os.path.join(__DIR, 'adi_ref_sensitvities.pickle'))
    (_, das_adi_spec, dmas_adi_spec, gd_adi_spec) = \
        load_pickle(os.path.join(__DIR, 'adi_ref_specificities.pickle'))

    # Load the sensitivities and specificities when
    # adipose-fibroglandular (healthy) scans are used for
    # reference subtraction
    (_, das_fib_sens, dmas_fib_sens, gd_fib_sens) = \
        load_pickle(os.path.join(__DIR, 'fib_ref_sensitvities.pickle'))
    (_, das_fib_spec, dmas_fib_spec, gd_fib_spec) = \
        load_pickle(os.path.join(__DIR, 'fib_ref_specificities.pickle'))

    (das_adi_auc, dmas_adi_auc, gd_adi_auc,
     das_fib_auc, dmas_fib_auc, gd_fib_auc) = get_aucs()

    # Make the figure
    plt.figure(figsize=(10, 8))
    plt.rc('font', family='Times New Roman')
    plt.tick_params(labelsize=(20))
    plt.gca().set_aspect('equal', adjustable='box')

    # Plot the ROC curves when adipose-only references are used
    plt.plot(100 - das_adi_spec, das_adi_sens, c=das_col, linestyle='-',
             label="DAS, Adipose Reference (AUC: %.1f%%)"
                   % (100 * das_adi_auc))
    plt.plot(100 - dmas_adi_spec, dmas_adi_sens, c=dmas_col, linestyle='--',
             label="DMAS, Adipose Reference (AUC: %.1f%%)"
                   % (100 * dmas_adi_auc))
    plt.plot(100 - gd_adi_spec, gd_adi_sens, c=gd_col, linestyle='--',
             label="ORR, Adipose Reference (AUC: %.1f%%)"
                   % (100 * gd_adi_auc))

    # Plot the ROC curves when healthy references are used
    plt.plot(100 - das_fib_spec, das_fib_sens, c=das_col2, linestyle='-',
             label="DAS, Healthy Reference (AUC: %.1f%%)"
                   % (100 * das_fib_auc))
    plt.plot(100 - dmas_fib_spec, dmas_fib_sens, c=dmas_col2, linestyle='--',
             label="DMAS, Healthy Reference (AUC: %.1f%%)"
                   % (100 * dmas_fib_auc))
    plt.plot(100 - gd_fib_spec, gd_fib_sens, c=gd_col2, linestyle='--',
             label="ORR, Healthy Reference (AUC: %.1f%%)"
                   % (100 * gd_fib_auc))

    # Plot the ROC curve of a random classifier
    plt.plot(np.linspace(0, 100, 100), np.linspace(0, 100, 100),
             c=[0, 108 / 255, 255 / 255], linestyle='--',
             label='Random Classifier')

    # Make the legend, axes labels, etc
    plt.legend(fontsize=16, loc='center right')
    plt.xlim([0, 100])
    plt.ylim([0, 100])
    plt.xlabel("False Positive Rate (%)", fontsize=22)
    plt.ylabel("True Positive Rate (%)", fontsize=22)
    plt.tight_layout()
    plt.show()
    # Save the fig
    # plt.savefig(os.path.join(__DIR, 'rocs.png'), dpi=300,
    #             transparent=True)

max_stars_repo_name: TysonReimer/ORR-Algorithm

max_issues_repo_name: TysonReimer/ORR-Algorithm

max_forks_repo_name: TysonReimer/ORR-Algorithm


text:
#!/usr/bin/python
import socket
import sys
import time
import traceback

import cv2
import numpy

capture = None

def runClient(sock):
    global capture

    # Set those constants for easy access
    TCP_IP = 'MAXPC.local'

    # Grab the port number from the command line
    TCP_PORT = int(sys.argv[1])

    # Connect to the socket with the previous information
    print("Connecting to Socket")
    sock.connect((TCP_IP, TCP_PORT))

    # Enable instant reconnection and disable timeout system
    sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)

    # Spread the good news
    print("Connected!")

    # If the constructor is an integer, make it one.  Otherwise leave it as a string.
    try:
        vid = int(sys.argv[2])
    except:
        vid = sys.argv[2]

    # Create the VideoCapture object with the adjusted parameter
    capture = cv2.VideoCapture(vid)

    # C R O N C H  that image for minimal latency
    capture.set(3, 480)
    capture.set(4, 360)

    while True:

        # Grab a frame from the webcam
        ret, frame = capture.read()

        # C R O N C H   that image down to half size
        newX, newY = frame.shape[1] * .5, frame.shape[0] * .5
        frame = cv2.resize(frame, (int(newX), int(newY)))

        # C R O N C H   that image down to extremely compressed JPEG
        encode_param = [int(cv2.IMWRITE_JPEG_QUALITY), 7]
        result, imgencode = cv2.imencode('.jpg', frame, encode_param)

        # Encode that JPEG string into a NumPy array for compatibility on the other side
        data = numpy.array(imgencode)

        # Turn NumPy array into a string so we can ship er' on over the information superhighway
        stringData = data.tostring()

        print(sys.getsizeof(stringData))

        # Send the size of the data for efficient unpacking
        sock.send(str(len(stringData)).ljust(16).encode())

        # Might as well send the actual data while we're sending things
        sock.send(stringData)

        # Arbitrary OpenCV wait statement that I still don't understand the purpose of but when I take it out it
        # doesn't work so here it will stay.
        k = cv2.waitKey(5) & 0xFF
        if k == 27:
            break

    cv2.waitKey(10)
    cv2.destroyAllWindows()

# This horrific piece of garbage just verifies that even if it crashes it gets up and tries again.  Please don't throw
# things at me.

def start():
    global capture

    while True:
        # Create the socket object
        sock = socket.socket()

        try:
            # Start the client
            runClient(sock)
        except Exception as e:
            # If it dies, give it a second to rest and force it to try again.
            sock.close()
            capture.release()
            print(traceback.format_exc())
            time.sleep(1)
            start()

# Start the Client
start()

max_forks_repo_name: rc-bellergy/pxpi


text:
import numpy as np
import os
import sys

def plot_separate(file_name, loss_list, title_list):
	import matplotlib.pyplot as plt
	from matplotlib.backends.backend_agg import FigureCanvasAgg
	from matplotlib.figure import Figure
	rows, cols, size = 1,1,5
	font_size = 30; label_size = 25; line_width = 2.5
	fig = Figure(tight_layout=True,figsize=(8, 6)); ax = fig.subplots(rows,cols)
	ax.plot(loss_list[0],linewidth = line_width);ax.plot(loss_list[1], linewidth=line_width)
	ax.set_ylabel(title_list[0], fontsize = font_size);ax.set_xlabel('Epochs',fontsize = font_size);ax.legend(['train','valid'],fontsize = font_size)
	ax.tick_params(axis = 'x', labelsize = label_size); ax.tick_params(axis = 'y', labelsize = label_size);
	ax.set_xlim([0,len(loss_list[0])])
	canvas = FigureCanvasAgg(fig); canvas.print_figure(file_name, dpi=100)

model_root_dir = '/data/models/report_results/'
model_names = os.listdir(model_root_dir)
# nb_epochs = 100
nb_epochs = 80
#model_name = 'single-net-Unet-bone-efficientnetb2-pre-True-epoch-200-batch-3-lr-0.0005-dim-992-train-1100-rot-0-set-cell_cycle_1984-ext-True-loss-focal+dice-up-upsampling-filters-256'
for model_name in model_names:
	train_dir = os.path.join(model_root_dir, model_name, 'train_dir')
	if os.path.exists(train_dir):
		train_loss = 4*np.loadtxt(train_dir+'/train_loss.txt')
		val_loss = 4*np.loadtxt(train_dir+'/val_loss.txt')

		train_dice = np.loadtxt(train_dir+'/train_f1-score.txt')
		val_dice = np.loadtxt(train_dir+'/val_f1-score.txt')
		if nb_epochs < len(train_loss):
			train_loss = train_loss[:nb_epochs]
			val_loss = val_loss[:nb_epochs]
			train_dice = train_dice[:nb_epochs]
			val_dice = val_dice[:nb_epochs]
		file_name = train_dir + '/loss_dice_history.png'
		loss_list = [train_loss, val_loss]
		metric_list = [train_dice, val_dice]
		#title_list = ['Focal loss', 'Dice score']
		#plot_history(file_name, loss_list, metric_list, title_list)
		file_name = train_dir + '/loss_history_{}.png'.format(nb_epochs)
		plot_separate(file_name, loss_list, ['Focal loss'])
		file_name = train_dir + '/dice_history.png'
		plot_separate(file_name, metric_list, ['Dice score'])


text:
# Copyright (C) 2018-2021 Intel Corporation
# SPDX-License-Identifier: Apache-2.0

import numpy as np

from openvino.tools.mo.front.common.replacement import FrontReplacementOp
from openvino.tools.mo.graph.graph import Graph
from openvino.tools.mo.ops.const import Const
from openvino.tools.mo.ops.strided_slice import StridedSlice

class MXSliceToStridedSliceReplacer(FrontReplacementOp):
    op = 'MXSlice'
    enabled = True

    def replace_sub_graph(self, graph: Graph, match: dict):
        node = match['op']

        strided_slice_node = StridedSlice(graph, dict(name=node.id + '/strided_slice_',
                                                      shrink_axis_mask=np.zeros(len(node.crop_begin), dtype=np.int64),
                                                      new_axis_mask=np.zeros(len(node.crop_begin), dtype=np.int64),
                                                      ellipsis_mask=np.zeros(len(node.crop_begin), dtype=np.int64),
                                                      begin_mask=np.ones(len(node.crop_begin), dtype=np.int64),
                                                      end_mask=np.ones(len(node.crop_end), dtype=np.int64))).create_node()
        node.in_port(0).get_connection().set_destination(strided_slice_node.in_port(0))
        node.out_port(0).get_connection().set_source(strided_slice_node.out_port(0))

        crop_begin_node = Const(graph, dict(value=node.crop_begin,
                                            symbol_dict={'name': node.id + '/crop_begin_const'})).create_node()
        crop_end_node = Const(graph, dict(value=node.crop_end,
                                          symbol_dict={'name': node.id + '/crop_end_const'})).create_node()
        strided_slice_node.in_port(1).get_connection().set_source(crop_begin_node.out_port(0))
        strided_slice_node.in_port(2).get_connection().set_source(crop_end_node.out_port(0))

        if len(node.step) > 0:
            stride_node = Const(graph, dict(value=node.step,
                                            symbol_dict={'name': node.id + '/steps_const'})).create_node()
            strided_slice_node.in_port(3).get_connection().set_source(stride_node.out_port(0))

max_stars_repo_name: pazamelin/openvino

max_issues_repo_name: pazamelin/openvino

max_forks_repo_name: pazamelin/openvino


text:
import cv2
import numpy as np
from scipy.signal import convolve
import time
from commonly_used import *
import matplotlib.pyplot as plt
from time import time
import scipy.io as sio
import os
import xlrd
from multiprocessing import Pool, Process
from math import exp
from scipy.signal import convolve, convolve2d
import torch
import torch.backends.cudnn as cudnn
from Utils import Rec_Loss, Augmentation, GaussPSF

# 相机透镜模型 u是物距
def uvf(v, f):
	u = (v * f) / (v - f)
	return u

# F镜头焦距、A光圈、fd focus depth、gtd ground truth depth、N 相机的f-number、p 像素尺寸到实际尺寸的比值、s 传感单元与相机输出尺寸大小、p s is the quintessence of this formula
def getCoC(F, A, fd, gtd):
	# A=F/N
	c = abs((gtd - fd)) / gtd * F / (fd - F)
	# print(c,end=' ')
	# c=c/(p*s)
	c = A * c
	# print(c)
	return c

# 生产高斯卷积核
def getGaussianKernel(n, sigma):
	kernel = np.zeros((n, n), dtype=np.float)
	# print(kernel.shape,kernel.dtype)
	for i in range(-n // 2, n // 2):
		for j in range(-n // 2, n // 2):
			kernel[i + n // 2][j + n // 2] = (1 / (2 * np.pi * sigma * sigma)) * np.exp(
				(-1) * (i * i + j * j) / (2 * sigma))
	# print(i,j,kernel[i+n//2][j+n//2])
	kernel = kernel / np.sum(kernel)
	return kernel

def load_fd(root, f, C, c0):
	fd = np.ones(C)
	v = np.ones(C)
	print("读取传感器数据......", end=' ')
	if os.path.isfile(root + '/fd.npy') != True or read == True:
		print("从excel读取：")
		v=np.load(root+'/fd/v.npy')
		#table = xlrd.open_workbook(root + '/fd/v.xlsx').sheets()[0]  # 读取excel文件中的第一个sheet
		for i in range(C):
			#v[i] = float(table.cell(i, 0).value[1:-1])
			# print(v[i],end=' ')
			fd[i] = uvf(v[i] + c0, f)
		# print(fd[i])
		np.save(root + '/fd/fd.npy', fd)
		np.save(root + '/fd/v.npy', v)
	else:
		print("读取现有数据")
		fd = np.load(root + '/fd.npy')
	return fd

def read_img(img, d):
	dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
	dftshift = np.fft.fftshift(dft)
	rows, cols = img.shape
	crow, ccol = int(rows / 2), int(cols / 2)
	mask = np.zeros((rows, cols, 2), np.uint8)
	mask[crow - d:crow + d, ccol - d:ccol + d] = 1
	fshit = dftshift * mask
	ishit = np.fft.ifftshift(fshit)
	iImg = cv2.idft(ishit)
	iImg = cv2.magnitude(iImg[:, :, 0], iImg[:, :, 1])
	plt.subplot(121), plt.imshow(img, cmap='gray')
	plt.title('orignal'), plt.axis('off')
	plt.subplot(122), plt.imshow(iImg, cmap='gray')
	plt.title('inverse'), plt.axis('off')
	plt.show()

def loadPosition(root):
	# 此函数用于读入三张包含坐标信息的图并压入padding 1
	# 将三张图表读入
	tic = time()
	print("读入每个点的Calib XYZ")
	load_path = os.path.join(os.getcwd(), root + "/3D/X.mat")
	X_mat = np.expand_dims(np.array(sio.loadmat(load_path)["X"]), 2)
	load_path = os.path.join(os.getcwd(), root + "/3D/Y.mat")
	Y_mat = np.expand_dims(np.array(sio.loadmat(load_path)["Y"]), 2)
	load_path = os.path.join(os.getcwd(), root + "/3D/Z.mat")
	Z_mat = np.expand_dims(np.array(sio.loadmat(load_path)["Z"]), 2)
	# 将三张图标融合在一起，total形状为[row, col, 3]
	total = np.append(np.append(X_mat, Y_mat, 2), Z_mat, 2)
	toc = time()
	print("读入完成\t\t用时%f秒" % (toc - tic))
	print("--------------------------------")
	# 输出total的形状为[row, col, 3]
	return total

def calGTDepth(RootPath, TransMat):
	# 此函数用于计算Ground True Depth
	calib_Pos = loadPosition(RootPath)
	[row, col, _] = calib_Pos.shape
	# 设定Padding， 将Padding融合至calib_Pos中
	Padding = np.ones([row, col, 1])
	# 此时calib_Pos形状为[row, col, 4]
	calib_Pos = np.append(calib_Pos, Padding, 2)
	Cam_Pos = np.matmul(calib_Pos, TransMat)
	depth = np.abs(Cam_Pos[:, :, 2])
	print('true depth shape', depth.shape)
	depth = cv2.resize(depth, img_shape)
	# 获得其中所有的深度种类
	ini_depth_list = np.unique(depth)
	min_depth, max_depth = np.min(ini_depth_list), np.max(ini_depth_list)
	depth_list = np.linspace(min_depth, max_depth, 400)
	print("所有的深度在(%f, %f)之间 " % (depth_list[0], depth_list[-1]))
	print("--------------------------------")
	np.save(RootPath + '/3D/depth.npy', depth)
	np.save(RootPath + '/3D/depth_list.npy', depth_list)
	return depth, depth_list

def depth_ini(shape):
	# 此函数用于深度的初始化
	# 方形深度，从400变化至649，共250种深度
	print("检查是否存在已有深度表")
	check = os.path.isfile("depth_400_649.npy") and os.path.isfile("depth_list.npy")
	if not check:
		print("制作深度表")
		tic = time()
		num_depth = shape[0] // 2
		# 设定一个存储depth标签和数值的两列空间，便于后面切片时寻找位置
		depth_list = np.zeros(num_depth)
		depth = np.zeros(shape)
		row, col = shape[0], shape[1]
		for i in range(num_depth):
			depth[i:row - i, i:col - i].fill(400 + i)
			depth_list[i] = 400 + i
		toc = time()
		print("深度表制作完成\t用时 %f秒" % (toc - tic))
		print("--------------------------------")
		np.save("depth_400_649.npy", depth)
		np.save("depth_list.npy", depth_list)

	else:
		print("读入深度表")
		print("--------------------------------")
		depth = np.load("depth_400_649.npy")
		depth_list = np.load("depth_list.npy")
	return depth, depth_list

def PSF(img, depth, F, A, ksize, fd):
	# 检查输入图像通道数 目前只支持单通道
	# 检查图像和深度图的形状是否匹配
	simu_img = np.ones((img.shape))
	for i in range(0, img.shape[0]):
		print(i)
		for j in range(0, img.shape[1]):
			c = getCoC(F=F, A=A, fd=fd, gtd=depth[i][j])
			if c != 0:
				kernel = getGaussianKernel(ksize, 500 * c)
				simu_img[i][j] = point_convolve(img, kernel, i, j)
			else:
				simu_img[i][j] = img[i][j]

	simu_img = simu_img.astype("uint8")
	# cv2.imshow("test",simu_img)
	# cv2.imwrite("psf"+str(fd)+".png",simu_img)
	# edge=cv2.Laplacian(simu_img,cv2.CV_8UC1,ksize=7)
	# cv2.imshow("edge",edge)

	return simu_img

def PSF_silce(img, depth, depth_list, F, A, w, fd, dd, id):
	num_depth = depth_list.shape[0]
	img_list = np.zeros(shape=[num_depth, img.shape[0], img.shape[1]])
	simu_img = np.ones(img.shape)
	# 把每一种depth带来的卷积操作执行完成
	print(id, end=',')
	# print("针对每一种深度，执行卷积操作")
	tic = time()
	for i in range(num_depth):
		# 针对depth_list种的每一个深度，寻找一次c，并做一次卷积操作
		c = getCoC(F=F, A=A, fd=fd, gtd=depth_list[i])
		# ksize=int(c)//2*2+1
		# print('A:',A,' C:',c)
		ksize = int((w * c) // 2) * 2 + 1
		img_list[i, :, :] = cv2.GaussianBlur(img, (ksize, ksize), c / 2)
	# ksize=int((w*c)//2)*2+1
	# img_list[i, :, :] = cv2.GaussianBlur(img, (w, w),c/2)
	# 开始挑选点
	toc = time()
	# print("卷积操作完成\t用时%f秒"%(toc-tic))
	# print("--------------------------------")
	row, col = img.shape
	# # 由于形状的特殊性，这里采用fill的方式填充
	# for i in range(row//2):
	#     simu_img[i:row-i, i:row-i] = img_list[i, i:row-i, i:row-i]

	# 对于一般情况, 采用以下这种方式
	# print("按照每个点的深度，依次寻找对应卷积值")
	tic = time()
	for i in range(row):
		for j in range(col):
			# 首先是读入这个点的深度 depth[i, j]
			cur_depth = depth[i, j]
			# 寻找index，并修改此处图像点
			index = np.argwhere(depth_list >= cur_depth)[0][0]
			simu_img[i, j] = img_list[index, i, j]
	toc = time()
	# print("查询完成\t\t用时%f秒"%(toc-tic))
	# print("--------------------------------")
	simu_img = simu_img.astype("uint8")
	# print(root+"/simu"+str(id)+".bmp")
	if not os.path.exists(os.path.join(root, "temp2")):
		# 判断文件夹是否存在，不存在则创建文件夹
		os.mkdir(os.path.join(root, "temp2"))
	cv2.imwrite(root + "/temp2/simu" + str(A) + '_' + str(w) + '_' + str(dd) + ".bmp", simu_img)
	return simu_img

def generate_simu_img(clear_img, depth, depth_list, fd, F, rangeA, rangew, ranged, num_pool=12):
	num = 0
	p = Pool(num_pool)
	for ww in rangew:
		for aa in rangeA:
			for dd in ranged:
				# 生成3张图片
				# print(aa,ww,dd)
				# print(num)
				p.apply_async(PSF_silce, args=(clear_img, depth, depth_list, F, aa, ww, fd[dd], dd, num))
				num += 1
	p.close()
	p.join()

# 输入全为CPU下numpy格式的数据
def generate_simu_img_cuda(clear_img, depth, fds, F, rangeA, rangew, ranged, theta):

	# ························第2种方式························#
	R = clear_img[:, :, 2]
	G = clear_img[:, :, 1]
	B = clear_img[:, :, 0]
	H = cv2.cvtColor(clear_img, cv2.COLOR_BGR2GRAY)
	weights_hsv = theta * 1 / (R / 700 + G / 546.1 + B / 435.8) / 180.038470
	weights_hsv = cv2.normalize(weights_hsv, None, alpha=0.8, beta=1.0, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
	# weights_hsv=np.ones(clear_img.shape)
	weights_hsv[:]=1
	print(weights_hsv)
	# cv2.imshow('hsl_imagev',weights_hsv)
	# cv2.waitKey(0)
	weights_hsv_temp = torch.Tensor(weights_hsv).float().cuda()
	# weights_hsv=hsl_image

	# ························第3种方式························#
	# B435.8 0.777，G5460.975   R700 1.248
	weights_rgb = np.array([1, 1, 1])
	weights_rgb = torch.Tensor(weights_rgb).float().cuda()

	# 清晰图归一化处理
	clear_img = cv2.normalize(clear_img, None, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
	clear_img = np.transpose(clear_img, (2, 0, 1))
	# print(clear_img.shape)
	image = torch.torch.zeros(1, 3, clear_img.shape[1], clear_img.shape[2])
	gt_depth = torch.torch.zeros(1, clear_img.shape[1], clear_img.shape[2])
	print(gt_depth.shape)
	print(image.shape)
	print(depth.shape)
	print(clear_img.shape)
	gt_depth[0, :, :] = torch.torch.from_numpy(depth)

	weights_hsv = torch.torch.zeros(1, img_shape[1], img_shape[0])
	print("weights_hsv.shape",weights_hsv.shape)
	weights_hsv[0, :, :] = weights_hsv_temp

	image[0] = torch.torch.from_numpy(clear_img)
	focal_length = torch.Tensor([F]).float().cuda()
	for w in range(len(rangew)):
		for ape in range(len(rangeA)):
			print(f'ape={rangeA[ape]}', end=',')
			for fd in range(len(ranged)):
				focal_depth = torch.Tensor([fds[ranged[fd]]]).float().cuda()
				aperture = torch.Tensor([rangeA[ape]]).float().cuda()
				# print(focal_depth)

				focused = generator(image, gt_depth, focal_depth, aperture, 5, focal_length, weights_hsv, weights_rgb)

				focused1 = (focused[0, :, :].permute(1, 2, 0).cpu().detach().numpy() * 255).astype('uint8')
				#focused1=np.rot90(np.flip(focused1,axis=1),1)
				cv2.imshow('focused1', focused1)
				cv2.waitKey(20)
				if not os.path.exists(os.path.join(root, "temp2")):
					os.mkdir(os.path.join(root, "temp2"))
				cv2.imwrite(
					root + "/temp2/simu" + str(rangeA[ape]) + '_' + str(theta) + '_' + str(ranged[fd]) + ".bmp",
					focused1)
		# print(ape,fd)
	cv2.destroyAllWindows()
	print('')

def load_gt_img(ranged):
	# print(len(ranged))
	N = len(ranged)
	gt_img = np.zeros((img_shape[1], img_shape[0], N), dtype='uint8')
	for i in range(N):
		temp = cv2.imread(root + '/' + img_name + '/' + img_name + '_' + str(ranged[i] + 1) + '.bmp', 0)

		gt_img[:, :, i] = cv2.resize(temp, img_shape)
		#cv2.imshow('test',gt_img[:, :, i] )
		#cv2.waitKey(0)
	return gt_img

def comp_loss(gt_img, rangeA, rangew, ranged, num_pool=5):
	loss = np.zeros((len(rangeA), 1))
	# np.save(root+'/temp_loss.npy',loss)
	# print(loss.shape)
	results = []
	p = Pool(num_pool)
	for a in range(len(rangeA)):
		#for w in range(len(rangew)):
			# print('a=',a,'w=',w)
		results.append(p.apply_async(loss_all, args=(root, gt_img, rangeA[a],rangew , ranged)))
	p.close()
	p.join()
	# print(results)
	num = 0
	for a in range(len(rangeA)):
		#for w in range(len(rangew)):
		loss[a, 0] = results[num].get()
		num += 1
	# print(loss[a,w])
	#print(loss)
	#print(loss.shape)
	##print(type(loss))
	return loss

# loss 多线程调用
def loss_all(root, gt_img, a, w, ranged):
	loss0, loss1, loss2, loss3 = 0, 0, 0, 0
	# loss=np.load(root+'/temp_loss.npy')
	# print('a=',a,' w=',w)
	for dd in range(len(ranged)):
		#print('a=',a,' w=',w,'d=',ranged[dd],end=' ')
		simu_img = cv2.imread(root + "/temp2/simu" + str(a) + '_' + str(w) + '_' + str(ranged[dd]) + ".bmp", 0)
		true_blur_img = gt_img[:, :, dd]
		#cv2.imshow('true'+str(a),true_blur_img)
		#cv2.imshow('simu'+str(a),simu_img)
		#cv2.waitKey(0)
		# print(simu_img.shape,true_blur_img.shape)
		# print(true_blur_img)
		loss0 += L2_loss(simu_img, true_blur_img) / (img_shape[0] * img_shape[1]) * 50000
		loss1 += defocus_loss(simu_img, true_blur_img)
		loss2 += defocus_num_loss(simu_img, true_blur_img) * 4
		loss3 += calRecError(simu_img, true_blur_img) * 5000
	print(loss0, loss1, loss2, loss3)
	# loss=(loss0+loss1+loss2)/len(ranged)
	loss = (loss0 + loss1 + loss2 + loss3) / (len(ranged) * 4)
	# loss=(loss3)/(len(ranged))

	print('a=', a, ' w=', w, ' loss=', loss)
	# np.save(root+'/temp_loss.npy',loss)
	return loss

# LOSS 之间计算F范数
def L2_loss(img1, img2):
	return (np.linalg.norm(img1.astype('float').flatten() - img2.astype('float').flatten(), ord=2))

# 对图像求一次二阶导然后求F范数
def defocus_loss(img1, img2):
	edge1 = cv2.Laplacian(img1, cv2.CV_64F, ksize=5)
	edge1 = cv2.convertScaleAbs(edge1)
	edge2 = cv2.Laplacian(img2, cv2.CV_64F, ksize=5)
	edge2 = cv2.convertScaleAbs(edge2)
	loss = np.linalg.norm(edge1.astype('float').flatten() - edge2.astype('float').flatten(), ord=2) / (
				img1.shape[0] * img1.shape[1] / 10000)
	return loss

# 统计每个清晰度的个数 并求loss 后期有空加入格式输入检查 优化算法
def defocus_num_loss(img1, img2):
	tic = time()
	# 拉普拉斯 图像清晰度可以看作是二阶导
	edge1 = cv2.Laplacian(img1, cv2.CV_64F, ksize=5)
	edge2 = cv2.Laplacian(img2, cv2.CV_64F, ksize=5)

	# 把Laplace得到的数据也就是defocus数据分为1000项
	classes = 1000
	minab = np.min((np.min(edge1), np.min(edge2)))
	maxab = np.max((np.max(edge1), np.max(edge2)))

	defocus_num_list = np.linspace(minab, maxab, classes + 1)
	# print(defocus_num_list,defocus_num_list.shape)
	# 统计每个清晰度区间的个数 很慢 不知道有没有更快的操作
	num1 = np.zeros(classes)
	num2 = np.zeros(classes)
	for i in range(classes):
		num1[i] = ((edge1 >= defocus_num_list[i]) & (edge1 <= defocus_num_list[i + 1])).sum()
		num2[i] = ((edge2 >= defocus_num_list[i]) & (edge2 <= defocus_num_list[i + 1])).sum()
	# print(num1)
	# print(num2)
	toc = time()
	# 归一化 取消图像大小对loss影响
	loss = L2_loss(num1, num2) / (img1.shape[0] * img1.shape[1] / 10000)
	# print(loss,(toc-tic),'s')
	return loss

def create_gaussian(window_size, sigma):
	# 此函数用于经计算SSIM中的高斯核
	gauss = np.array([exp(-(x - window_size // 2) ** 2 / float(2 * sigma ** 2)) for x in range(window_size)])
	gauss = gauss / np.sum(gauss)
	window_1 = np.expand_dims(gauss, 1)
	window = np.matmul(window_1, window_1.transpose())
	return window

def calSSIM(img1, img2, window_size, sigma=1.5):
	# 此函数用于计算SSIM
	window = create_gaussian(window_size, sigma)
	# 分别求公式中的mu 和 sigma
	mu1 = convolve2d(img1, window, "same")
	mu2 = convolve2d(img2, window, "same")
	mu1_sq = np.power(mu1, 2)
	mu2_sq = np.power(mu2, 2)
	mu1_mu2 = mu1 * mu2

	sigma1_sq = convolve2d(img1 * img1, window, "same") - mu1_sq
	sigma2_sq = convolve2d(img2 * img2, window, "same") - mu2_sq
	sigma12 = convolve2d(img1 * img2, window, "same") - mu1_mu2

	C1 = 0.01 ** 2
	C2 = 0.03 ** 2
	ssim = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2)) / ((mu1_sq + mu2_sq + C1) * (sigma1_sq + sigma2_sq + C2))
	return ssim

def calRecError(img1, img2, window_size=11, sigma=1.5):
	alpha = 1
	img1 = img1.astype('float')
	img2 = img2.astype('float')
	ssim_loss = np.mean(1 - calSSIM(img1, img2, window_size, sigma=1.5)) / 2
	rec_loss = np.linalg.norm(img1.flatten() - img2.flatten(), ord=1)
	return alpha * ssim_loss + (1 - alpha) * rec_loss

roots = [
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.27true1',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.27true2',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.27true3',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.27true4',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.27true5',

	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.25boox1',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.25boox2',

	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.23mask1',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.23mask2',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.23mask3',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.23mask4',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.23mask5',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.27true6',

	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.07mask1',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.07mask2',
	'F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/9.07mask3'

]

#np_old_load=np.load

# img_name=root[-5:]
# print(img_name)
# img_shape = (500, 500)
# root='F:/OneDrive - stu.cqupt.edu.cn/code/SFF_code/data/7.24mask1'
# root='I:/heka/SFF/PSF/data/7.24mask1'
# img_name='mask1'
img_shape = ( 613,514)
#img_shape = ( 514,613)
#
if __name__ == '__main__':

	for i in range(len(roots)):
		global root, img_name
		root = roots[i]
		img_name = root[-5:]
		# ············cuda初始化············#
		cudnn.benchmark = True
		generator = GaussPSF(7, near=1e-5, far=1e7)
		generator = torch.nn.DataParallel(generator)
		generator = generator.cuda()
		torch.set_num_threads(1)

		# ···························································读入清晰图片 clear image···························································#
		print('读入清晰图片')
		tic = time()
		clear_img = cv2.imread(root + "/clear.tif")
		print(root)
		print(clear_img.shape)
		clear_img = cv2.resize(clear_img, img_shape)
		toc = time()
		print(clear_img.shape)
		print(f'读入完成    耗时{(toc - tic):.2f}s')
		print("--------------------------------")
		# ···························································读入深度 ground truth depth···························································#

		if os.path.isfile(root + '/3D/depth.npy') != True :
			depth, depth_list = calGTDepth(root, Trans_mat0907)
			print(depth.shape)
		else:
			print("读取现有数据")
			depth = np.load(root + '/3D/depth.npy',allow_pickle=True)
			depth_list = np.load(root + '/3D/depth_list.npy',allow_pickle=True)
			depth=cv2.resize(depth,img_shape)
			print("读取完毕!")
			print("--------------------------------")

		#depth=np.flip(depth,axis=1)
		#cv2.imshow('test', depth.astype('uint8'))
		#cv2.waitKey(0)
		#fig = show_3D_picture(depth)
		#depth=depth.T

		depth=np.flip(depth,axis=1)
		#depth=np.flip(depth,axis=0)
		#depth = np.rot90(depth, 2)

		#fig=show_3D_picture(depth[000:500,000:500])
		#depth = depth.T
		print('depthshape',depth.shape)
		#depth=np.rot90(depth,-1)
		depth = np.ascontiguousarray(depth)
		print('depthshape',depth.shape)

		cv2.imshow("test_depth",cv2.normalize(depth, None, alpha=1, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8UC1))
		cv2.imshow('clear',clear_img)
		cv2.waitKey(10)
		#cv2.destroyAllWindows()
		#fig=show_3D_picture(depth)

		#fig=show_3D_picture(depth)
		# depth2=np.flip(depth2,axis=0)
		# fig=show_3D_picture(depth2)
		# ···························································读入焦深 focus depth···························································#
		read = 1  # 1 强制重新读取 0 读入现有数据
		F = 50
		C = 200
		c0 = 24.5
		fd = load_fd(root, F, C, c0)
		print("读取完毕!")
		print(f"所有的focus depth在({fd[0]}, {fd[-1]})之间")
		print("--------------------------------")

		# ··························································暴力搜索A w范围···························································#
		F = 50  # 设置镜头参数
		# 设置求解参数的范围 A为光圈参数A/(p*s) w为卷积核大小 使用高斯函数时只能为单数 d为选择的照片序列
		rangeA = np.array(range(400, 1200, 50))
		rangew = np.array(range(4, 5, 1))
		# ranged=np.array(range(0,200,10))
		ranged = np.array(range(0, 200, 10))
		rangec = np.array(range(0, 10, 1))

		# 初始化loss矩阵
		loss = np.zeros((len(rangeA), len(rangec)))
		print(loss.shape)
		tic = time()

		loss_c = np.zeros(len(rangec))
		# 可以对每个w计算一次loss 也可以全部把图片求出来之后再计算loss 前一种方法实时性较好
		for cc in range(len(rangec)):
			# for i in range(len(rangew)):
			c0=23+cc*0.2
			#c0 = 24.2
			fds = load_fd(root, F, C, c0)
			print(f"c0为{c0},所有的focus depth在({fds[0]}, {fds[-1]})之间")
			# 生成模拟模糊图片
			#generate_simu_img(clear_img,depth,depth_list,fds,F,rangeA,rangew,ranged,num_pool=3)
			theta = 1 + cc
			print(f'theta:{theta}')
			#generate_simu_img_cuda(clear_img, depth, fds, F, rangeA, rangew, ranged, theta)
			# 载入真实模拟图片
			gt_img = load_gt_img(ranged)
			# 计算loss
			# temploss= comp_loss(gt_img,rangeA,rangew[i:i+1],ranged)
			temploss = comp_loss(gt_img, rangeA, theta, ranged, num_pool=5)
			# 通过第一种方法计算出来只是一条loss 需要拼接一次
			# loss [:,i]=temploss[:,0]
			# np.save(root+'/preloss.npy',loss)
			#print(type(temploss))
			#print(temploss.shape)
			loss[:, cc] = temploss[:,0]
			np.save(root + '/preloss.npy', loss)
			print('当前最小值为', np.min(loss[:, 0:cc + 1]), '位于', np.where(loss == np.min(loss[:, 0:cc + 1])))

		np.save(root + '/preloss.npy', loss)
		print(np.min(loss), np.where(loss == np.min(loss)))

reason: import numpy,import scipy,from scipy


text:
import pandas as pd
import lightgbm as lgb
from sklearn.metrics import roc_auc_score
import numpy as np

print("Loading resampled train data")
train_X = pd.read_csv("../input/AllData_v4_os.train")
train_X.pop("Unnamed: 0")

print("Loading resampled train labels")
train_y = pd.read_csv("../input/AllData_v4_os.label")
train_y = train_y.pop("TARGET")

print("Loading resampled validation data")
valid_X = pd.read_csv("../input/AllData_v4_os_valid.train")
valid_X.pop("Unnamed: 0")

print("Loading resampled validation labels")
valid_y = pd.read_csv("../input/AllData_v4_os_valid.label")
valid_y = valid_y.pop("TARGET")

print("Loading application test data")
test_X = pd.read_csv("../input/AllData_v4.test")

print("Replacing Inf values in valid_X")
valid_X = valid_X.replace([np.inf, -np.inf], np.nan)
print("Replacing NaN values in valid_X")
valid_X = valid_X.fillna(0)

print("Replacing Inf values in test_X")
test_X = test_X.replace([np.inf, -np.inf], np.nan)
print("Replacing NaN values in test_X")
test_X = test_X.fillna(0)

print("train_y shape: " + str(train_y.shape))
print("train_X shape: " + str(train_X.shape))
print("valid_y shape: " + str(valid_y.shape))
print("valid_X shape: " + str(valid_X.shape))
print("test_X shape: " + str(test_X.shape))

lgb_train = lgb.Dataset(train_X, train_y)
lgb_test = lgb.Dataset(valid_X)

params = {'task'                    :'train',
          'objective'               :'binary',
          'learning_rate'           :0.01,
          'num_leaves'              :10,
          'max_depth'               :3,
          'min_data_in_leaf'        :80,
          'min_sum_hessian_in_leaf' :0.001,
          'lambda_l1'               :0.2,
          'lambda_l2'               :0,
          'scale_pos_weight'        :1,
          'metric'                  :'auc',
          'verbose'                 :-1}

print("Model training started...")
model   = lgb.train(params=params,
                    train_set=lgb_train,
                    num_boost_round=10000,
                    verbose_eval=True)
print("Model training completed...")

model.save_model('GridSearch/AllData_v4_OS_LGBM_v2_Model')

print("Predicting validation set...")
valid_preds = model.predict(valid_X)
print("Validation set prediction completed...")

print("Predicting test set...")
test_preds = model.predict(test_X)
print("Test set prediction completed...")

auc = roc_auc_score(valid_y, valid_preds)
print("Validation AUC: " + str(auc))

valid_preds = pd.DataFrame(valid_preds)
valid_preds.to_csv("GridSearch/AllData_v4_OS_LGBM_ValidPreds_v2.csv", index=False)

sub = pd.read_csv('../input/sample_submission.csv')
sub['TARGET'] = test_preds
sub.to_csv('GridSearch/AllData_v4_OS_LGBM_Preds_v2.csv', index=False)