text:
import torch
# import torch.nn as nn
import argparse
# import torch.optim as optim
import torch.backends.cudnn as cudnn
from torch.utils.data import DataLoader
from torch.autograd import Variable
import os
import time
import random
from utils.dataset import HyperDatasetValid, HyperDatasetTrain  # Clean Data set
from models.models import create_model
# from utils.util import AverageMeter, initialize_logger, save_checkpoint, record_loss, LossTrainCSS, Loss_valid
from utils.visualizer import Visualizer
from utils.util import initialize_logger
from utils.EvalMetrics import computeMRAE
# from collections import OrderedDict
import numpy as np
from skimage.measure import compare_ssim
from skimage.measure import compare_psnr

parser = argparse.ArgumentParser(description="SSRGAN")
parser.add_argument("--batchSize", type=int, default=1, help="input batch size")
parser.add_argument("--SEED", type=int, default=2020, help="random seed")
parser.add_argument("--end_epoch", type=int, default=200+1, help="number of epochs")
parser.add_argument("--init_lr", type=float, default=1e-4, help="initial learning rate")
parser.add_argument("--decay_power", type=float, default=1.5, help="decay power")
parser.add_argument("--trade_off", type=float, default=10, help="trade_off")
parser.add_argument("--max_iter", type=float, default=300000, help="max_iter")  # patch48:380x450/32x100-534375; patch96:82x450/32x100-113906
parser.add_argument("--outf", type=str, default="SSRGAN_Results", help='path log files')
# custom options
parser.add_argument('--name', type=str, default='RGB2HSI', help='name of the experiment. It decides where to store samples and models')
parser.add_argument('--gpu_ids', type=str, default='0', help='gpu ids: e.g. 0  0,1,2, 0,2. use -1 for CPU')
parser.add_argument('--isTrain', type=bool, default=True, help='isTrain')
parser.add_argument('--checkpoints_dir', type=str, default='./checkpoints', help='models are saved here')
parser.add_argument('--continue_train', action='store_true', help='continue training: load the latest model')
parser.add_argument('--load_pretrain', type=str, default='', help='load the pretrained model from the specified location')
parser.add_argument('--which_epoch', type=str, default='latest', help='which epoch to load? set to latest to use latest cached model')
parser.add_argument('--display_freq', type=int, default=5, help='frequency of showing training results on screen')
parser.add_argument('--print_freq', type=int, default=100, help='frequency of showing training results on console')
parser.add_argument('--verbose', action='store_true', default=False, help='toggles verbose')

parser.add_argument('--instance_feat', action='store_true', help='if specified, add encoded instance features as input')
parser.add_argument('--label_feat', action='store_true', help='if specified, add encoded label features as input')
parser.add_argument('--load_features', action='store_true', help='if specified, load precomputed feature maps')
parser.add_argument('--niter', type=int, default=100, help='# of iter at starting learning rate')
parser.add_argument('--niter_decay', type=int, default=100, help='# of iter to linearly decay learning rate to zero')
parser.add_argument('--beta1', type=float, default=0.5, help='momentum term of adam')
parser.add_argument('--lr', type=float, default=0.0002, help='initial learning rate for adam')
parser.add_argument('--save_latest_freq', type=int, default=200, help='frequency of saving the latest results')
parser.add_argument('--save_epoch_freq', type=int, default=1, help='frequency of saving checkpoints at the end of epochs')
parser.add_argument('--loadSize', type=int, default=1024, help='scale images to this size')
parser.add_argument('--fineSize', type=int, default=512, help='then crop to this size')
parser.add_argument('--no_flip', action='store_true', default=False, help='if specified, do not flip the images for data argumentation')
parser.add_argument('--label_nc', type=int, default=35, help='# of input label channels')
parser.add_argument('--input_nc', type=int, default=3, help='# of input image channels')
parser.add_argument('--output_nc', type=int, default=31, help='# of output HSI channels')

parser.add_argument('--no_instance', action='store_true', help='if specified, do *not* add instance map as input')
parser.add_argument('--feat_num', type=int, default=3, help='vector length for encoded features')
parser.add_argument('--n_downsample_E', type=int, default=4, help='# of downsampling layers in encoder')
parser.add_argument('--nef', type=int, default=16, help='# of encoder filters in the first conv layer')
parser.add_argument('--n_clusters', type=int, default=10, help='number of clusters for features')

parser.add_argument('--netG', type=str, default='global', help='selects model to use for netG')
parser.add_argument('--ngf', type=int, default=64, help='# of gen filters in first conv layer')
parser.add_argument('--n_downsample_global', type=int, default=4, help='number of downsampling layers in netG')
parser.add_argument('--n_blocks_global', type=int, default=9, help='number of residual blocks in the global generator network')
parser.add_argument('--n_blocks_local', type=int, default=3, help='number of residual blocks in the local enhancer network')
parser.add_argument('--n_local_enhancers', type=int, default=1, help='number of local enhancers to use')
parser.add_argument('--niter_fix_global', type=int, default=0, help='number of epochs that we only train the outmost local enhancer')
parser.add_argument('--norm', type=str, default='instance', help='instance normalization or batch normalization')

parser.add_argument('--num_D', type=int, default=2, help='number of discriminators to use')
parser.add_argument('--n_layers_D', type=int, default=3, help='only used if which_model_netD==n_layers')
parser.add_argument('--ndf', type=int, default=64, help='# of discrim filters in first conv layer')
parser.add_argument('--lambda_feat', type=float, default=10.0, help='weight for feature matching loss')
parser.add_argument('--no_ganFeat_loss', action='store_true', help='if specified, do *not* use discriminator feature matching loss')
parser.add_argument('--no_vgg_loss', action='store_true', help='if specified, do *not* use VGG feature matching loss')
parser.add_argument('--no_lsgan', action='store_true', help='do *not* use least square GAN, if false, use vanilla GAN')
parser.add_argument('--pool_size', type=int, default=0, help='the size of image buffer that stores previously generated images')

opt = parser.parse_args()

os.environ["CUDA_DEVICE_ORDER"] = 'PCI_BUS_ID'
# os.environ["CUDA_VISIBLE_DEVICES"] = opt.gpu_ids

torch.manual_seed(opt.SEED)
torch.cuda.manual_seed(opt.SEED)
np.random.seed(opt.SEED)
random.seed(opt.SEED)

def main():
    cudnn.benchmark = True

    # load dataset
    print("\nloading dataset ...")
    train_data = HyperDatasetTrain(mode='train', opt=opt)
    print("Train datasets: %d " % (len(train_data)))
    val_data = HyperDatasetValid(mode='valid', opt=opt)
    print("Validation set samples: ", len(val_data))

    train_loader = DataLoader(dataset=train_data, batch_size=opt.batchSize, shuffle=True, num_workers=2, pin_memory=True, drop_last=True)
    val_loader = DataLoader(dataset=val_data, batch_size=1,  shuffle=False, num_workers=2, pin_memory=True)

    # model
    print("\nbuilding models_baseline ...")
    model = create_model(opt)
    print('Parameters number is %.3f M' % (sum(param.numel() for param in model.parameters())/1e6))
    visualizer = Visualizer(opt)

    optimizer_G, optimizer_D = model.optimizer_G, model.optimizer_D

    # Parameters, Loss and Optimizer
    start_epoch = 1 if not opt.continue_train else opt.which_epoch
    # iteration = 0
    total_steps = (start_epoch-1) * len(train_data)
    max_PSNR, max_epoch = 0, 0

    # visualzation
    os.makedirs(opt.outf, exist_ok=True)
    log_dir = os.path.join(opt.checkpoints_dir, opt.name, 'train_val_metric.log')
    logger = initialize_logger(log_dir)

    iter_path = os.path.join(opt.checkpoints_dir, opt.name, 'iter.txt')

    # start epoch
    for epoch in range(start_epoch, opt.end_epoch):
        start_time = time.time()
        train_PSNR, train_SSIM, val_PSNR, val_SSIM, val_MRAE, max_PSNR, max_epoch = train(train_loader, val_loader, model,
                                                                                optimizer_G, optimizer_D, epoch,
                                                                                total_steps, start_time, visualizer,
                                                                                iter_path, max_PSNR, max_epoch)
        end_time = time.time()
        epoch_time = end_time - start_time
        logger.info(" Epoch [%02d], Time:%.9f, Train PSNR: %.4f SSIM: %.4f Val PSNR: %.4f SSIM: %.4f MRAE: %.4f"
                    % (epoch, epoch_time, train_PSNR, train_SSIM, val_PSNR, val_SSIM, val_MRAE))

# Training
def train(train_loader, val_loader, model, optimizer_G, optimizer_D, epoch, total_steps, start_time, visualizer, iter_path, max_PSNR, max_epoch):
    save_fake = True if opt.display_freq else False
    psnr_list = []
    ssim_list = []
    mrae_list = []
    for i, (rgb, hyper) in enumerate(train_loader):
        total_steps += 1
        # print("rgb shape {} hyper shape {}".format(rgb.shape, hyper.shape))
        # rgb shape torch.Size([1, 3, 482, 512]) hyper shape torch.Size([1, 31, 482, 512])
        # print("rgb max {} rgb min {} hyper max {} hyper min {}".format(rgb.max(), rgb.min(), hyper.max(), hyper.min()))
        # rgb max 1.0 rgb min -1.0 hyper max 1.0 hyper min -0.9982895851135254

        # Forward Pass
        losses, generated = model(Variable(rgb),
                                  Variable(hyper),
                                  infer=save_fake)

        # sum per device losses
        losses = [torch.mean(x) if not isinstance(x, int) else x for x in losses]
        loss_dict = dict(zip(model.loss_names, losses))
        # print("loss_dict: ", loss_dict)
        # loss_dict:  {'G_GAN': tensor(12.2738, device='cuda:0', grad_fn=<MeanBackward0>),
        # 'G_GAN_Feat': tensor(15.5414, device='cuda:0', grad_fn=<MeanBackward0>),
        # 'G_CSS': tensor(4.9444, device='cuda:0', grad_fn=<MeanBackward0>),
        # 'D_real': tensor(5.7572, device='cuda:0', grad_fn=<MeanBackward0>),
        # 'D_fake': tensor(7.2036, device='cuda:0', grad_fn=<MeanBackward0>)}

        # calculate final loss scalar
        loss_D = (loss_dict['D_fake'] + loss_dict['D_real']) * 0.5
        loss_G = loss_dict['G_GAN'] + loss_dict.get('G_GAN_Feat', 0) + loss_dict.get('G_CSS', 0)

        # Backward Pass #
        # update generator weights
        optimizer_G.zero_grad()
        loss_G.backward()
        optimizer_G.step()

        # update discriminator weights
        optimizer_D.zero_grad()
        loss_D.backward()
        optimizer_D.step()

        # Display results and errors #
        psnr, ssim, l1, mrae = metrics(generated.cuda(), hyper.cuda())
        psnr_list.append(psnr.item())
        ssim_list.append(ssim.item())
        mrae_list.append(mrae.item())
        # print out errors
        if i % opt.display_freq == 0:
            errors = {k: v.data.item() if not isinstance(v, int) else v for k, v in loss_dict.items()}
            t = (time.time() - start_time) / opt.display_freq
            visualizer.print_current_errors(epoch, i, errors, t)
            # visualizer.plot_current_errors(errors, total_steps)

        # display output images
        if save_fake and i % opt.display_freq == 0:
            visualizer.display_samples(rgb, hyper, generated, epoch, total_steps, i, mode=True)

        # # save latest model
        # if total_steps % opt.save_latest_freq == 0:
        #     print('saving the latest model (epoch %d, total_steps %d)' % (epoch, total_steps))
        #     model.save('latest')
        #     np.savetxt(iter_path, (epoch, i), delimiter=',', fmt='%d')

    avg_psnr = np.average(psnr_list)
    avg_ssim = np.average(ssim_list)
    avg_mrae = np.average(mrae_list)
    print('End of epoch %d / %d \t PSNR %.4f SSIM %.4f MARE %.4f Time Taken: %d sec' %
          (epoch, opt.end_epoch, avg_psnr, avg_ssim, avg_mrae, time.time() - start_time))

    print("--- Start Validation ---")
    val_PSNR, val_SSIM, val_MRAE = validate(val_loader, model, epoch, total_steps, visualizer)

    # save model for this epoch
    if epoch % opt.save_epoch_freq == 0 and val_PSNR > max_PSNR:
        max_PSNR = val_PSNR
        max_epoch = epoch
        # print('Saving the model at the end of epoch %d, iters %d' % (epoch, total_steps))
        print('Saving the model at the end of epoch %d' % (epoch))
        model.save(epoch)
        np.savetxt(iter_path, (epoch+1, 0), delimiter=',', fmt='%d')
    model.save('latest')

    print('The best epoch is %d when max_PSNR is %.4f' % (max_epoch, max_PSNR))

    # instead of only training the local enhancer, train the entire network after certain iterations
    if (opt.niter_fix_global != 0) and (epoch == opt.niter_fix_global):
        model.update_fixed_params()

    # linearly decay learning rate after certain iterations
    if epoch > opt.niter:
        model.update_learning_rate()

    return avg_psnr, avg_ssim, val_PSNR, val_SSIM, val_MRAE, max_PSNR, max_epoch

# Validate
def validate(val_loader, model, epoch, total_steps, visualizer):
    start_time = time.time()
    psnr_list = []
    ssim_list = []
    mrae_list = []
    for i, (rgb, hyper) in enumerate(val_loader):
        # print("rgb max {} rgb min {} hyper max {} hyper min {}".format(rgb.max(), rgb.min(), hyper.max(), hyper.min()))
        generated = model.inference(Variable(rgb), Variable(hyper))
        visualizer.display_samples(rgb, hyper, generated, epoch, total_steps, i, mode=False)
        psnr, ssim, l1, mrae = metrics(generated.cuda(), hyper.cuda())
        psnr_list.append(psnr.item())
        ssim_list.append(ssim.item())
        mrae_list.append(mrae.item())

    avg_psnr = np.average(psnr_list)
    avg_ssim = np.average(ssim_list)
    avg_mrae = np.average(mrae_list)
    print('Epoch %d \t Validation PSNR %.4f SSIM %.4f MRAE %.4f Time Taken: %d sec' %
          (epoch, avg_psnr, avg_ssim, avg_mrae, time.time() - start_time))

    return avg_psnr, avg_ssim, avg_mrae

def metrics(inputs, gts):
    # inputs = self.postprocess(inputs)
    # gts = self.postprocess(gts)
    psnr_value = []
    l1_value = torch.mean(torch.abs(inputs-gts))

    [b, n, w, h] = inputs.size()
    inputs = (inputs*255.0).int().float()/255.0
    gts = (gts*255.0).int().float()/255.0

    for i in range(inputs.size(0)):
        inputs_p = inputs[i, :, :, :].cpu().numpy().astype(np.float32).transpose(1, 2, 0)
        gts_p = gts[i, :, :, :].cpu().numpy().astype(np.float32).transpose(1, 2, 0)
        psnr_value.append(compare_psnr(inputs_p, gts_p, data_range=1))

    psnr_value = np.average(psnr_value)
    inputs = inputs.view(b*n, w, h).cpu().numpy().astype(np.float32).transpose(1, 2, 0)
    gts = gts.view(b*n, w, h).cpu().numpy().astype(np.float32).transpose(1, 2, 0)
    ssim_value = compare_ssim(inputs, gts, data_range=1, win_size=51, multichannel=True)
    # print("inputs value: ", inputs.max(), inputs.min(), "gts value: ", gts.max(), gts.min())
    norm_inputs = (inputs - inputs.min()) / (inputs.max() - inputs.min())
    norm_gts = (gts - gts.min()) / (gts.max() - gts.min())
    # print("norm_gts value: ", norm_inputs.max(), norm_inputs.min(), "norm_gts value: ", norm_gts.max(), norm_gts.min())
    # print("norm_gts value: ", norm_inputs, "norm_gts value: ", norm_gts)
    mrae = computeMRAE(norm_inputs, norm_gts)

    return psnr_value, ssim_value, l1_value, mrae

if __name__ == '__main__':
    print(torch.__version__)
    str_ids = opt.gpu_ids.split(',')
    opt.gpu_ids = []
    for str_id in str_ids:
        id = int(str_id)
        if id >= 0:
            opt.gpu_ids.append(id)
    # set gpu ids
    if len(opt.gpu_ids) > 0:
        torch.cuda.set_device(opt.gpu_ids[0])
    main()


text:
"""
Write a function, persistence, that takes in a positive parameter num and returns its multiplicative persistence,
which is the number of times you must multiply the digits in num until you reach a single digit.

For example:

 persistence(39) === 3 // because 3*9 = 27, 2*7 = 14, 1*4=4
                       // and 4 has only one digit
 persistence(999) === 4 // because 9*9*9 = 729, 7*2*9 = 126,
                        // 1*2*6 = 12, and finally 1*2 = 2
 persistence(4) === 0 // because 4 is already a one-digit number
"""
import numpy as np

count = 0

def persistence(number):
    digits_list = []
    if number < 0 or isinstance(number, str):
        print('Your input is incorrect!')
    elif number < 10:
        print("Persistance for your integer is: 0")
    else:
        for digit in str(number):
            digits_list.append(int(digit))
        result = ((np.prod(np.array(digits_list))).item())
        result_list = []
        for elem in str(result):
            result_list.append(int(elem))
        print(result_list)
        global count
        count += 1
        if len(result_list) != 1:
            persistence(result)
        else:
            print("Persistance for your integer is: ", count)

persistence(999)

ext: py

reason: import numpy


text:
'''
Copyright 2022 Airbus SAS

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.
'''
from climateeconomics.core.core_witness.climateeco_discipline import ClimateEcoDiscipline
from sos_trades_core.tools.post_processing.charts.chart_filter import ChartFilter
from sos_trades_core.tools.post_processing.charts.two_axes_instanciated_chart import InstanciatedSeries,\
    TwoAxesInstanciatedChart
from climateeconomics.core.core_forest.forest_v1 import Forest
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from sos_trades_core.tools.post_processing.plotly_native_charts.instantiated_plotly_native_chart import \
    InstantiatedPlotlyNativeChart

class ForestDiscipline(ClimateEcoDiscipline):
    ''' Forest discipline
    '''

    # ontology information
    _ontology_data = {
        'label': 'Forest',
        'type': '',
        'source': '',
        'validated': '',
        'validated_by': '',
        'last_modification_date': '',
        'category': '',
        'definition': '',
        'icon': 'fas fa-tree fa-fw',
        'version': '',
    }
    default_year_start = 2020
    default_year_end = 2050

    deforestation_limit = 1000
    initial_emissions = 3.21

    DESC_IN = {'year_start': ClimateEcoDiscipline.YEAR_START_DESC_IN,
               'year_end': ClimateEcoDiscipline.YEAR_END_DESC_IN,
               'time_step': ClimateEcoDiscipline.TIMESTEP_DESC_IN,
               Forest.DEFORESTATION_SURFACE: {'type': 'dataframe', 'unit': 'Mha',
                                                      'dataframe_descriptor': {'years': ('float', None, False),
                                                                               'deforested_surface': ('float', [0, 1e9], True)}, 'dataframe_edition_locked': False,
                                                      'visibility': ClimateEcoDiscipline.SHARED_VISIBILITY, 'namespace': 'ns_witness'},
               Forest.LIMIT_DEFORESTATION_SURFACE: {'type': 'float', 'unit': 'Mha', 'default': deforestation_limit,
                                                            'namespace': 'ns_forest', },
               Forest.INITIAL_CO2_EMISSIONS: {'type': 'float', 'unit': 'GtCO2', 'default': initial_emissions,
                                                      'namespace': 'ns_forest', },
               Forest.CO2_PER_HA: {'type': 'float', 'unit': 'kgCO2/ha/year', 'default': 4000, 'namespace': 'ns_forest'},
               Forest.REFORESTATION_COST_PER_HA: {'type': 'float', 'unit': '$/ha', 'default': 3800, 'namespace': 'ns_forest'},
               Forest.REFORESTATION_INVESTMENT: {'type': 'dataframe', 'unit': 'G$',
                                                         'dataframe_descriptor': {'years': ('float', None, False),
                                                                                'forest_investment': ('float', [0, 1e9], True)}, 'dataframe_edition_locked': False,
                                                         'visibility': ClimateEcoDiscipline.SHARED_VISIBILITY, 'namespace': 'ns_invest'},
               }

    DESC_OUT = {
        'CO2_emissions_detail_df': {
            'type': 'dataframe', 'unit': 'GtCO2', 'namespace': 'ns_forest'},
        Forest.FOREST_SURFACE_DF: {
            'type': 'dataframe', 'unit': 'Gha', 'visibility': ClimateEcoDiscipline.SHARED_VISIBILITY, 'namespace': 'ns_witness'},
        Forest.FOREST_DETAIL_SURFACE_DF: {
            'type': 'dataframe', 'unit': 'Gha'},
        Forest.CO2_EMITTED_FOREST_DF: {
            'type': 'dataframe', 'unit': 'GtCO2', 'visibility': ClimateEcoDiscipline.SHARED_VISIBILITY, 'namespace': 'ns_witness'},
    }

    FOREST_CHARTS = 'Forest chart'

    def init_execution(self):
        inputs = list(self.DESC_IN.keys())
        param = self.get_sosdisc_inputs(inputs, in_dict=True)

        self.forest_model = Forest(param)

    def run(self):

        #-- get inputs
        #         inputs = list(self.DESC_IN.keys())
        #         inp_dict = self.get_sosdisc_inputs(inputs, in_dict=True)

        #-- compute
        in_dict = self.get_sosdisc_inputs()
        self.forest_model.compute(in_dict)

        outputs_dict = {
            Forest.CO2_EMITTED_DETAIL_DF: self.forest_model.CO2_emitted_df,
            Forest.FOREST_DETAIL_SURFACE_DF: self.forest_model.forest_surface_df,
            Forest.FOREST_SURFACE_DF: self.forest_model.forest_surface_df[['years', 'forest_surface_evol']],
            Forest.CO2_EMITTED_FOREST_DF: self.forest_model.CO2_emitted_df[['years', 'emitted_CO2_evol_cumulative']],
        }

        #-- store outputs
        self.store_sos_outputs_values(outputs_dict)

    def compute_sos_jacobian(self):
        """
        Compute jacobian for each coupling variable
        """
        in_dict = self.get_sosdisc_inputs()
        self.forest_model.compute(in_dict)

        # gradient for deforestation rate
        d_deforestation_surface_d_deforestation_surface = self.forest_model.d_deforestation_surface_d_deforestation_surface()
        d_cum_deforestation_d_deforestation_surface = self.forest_model.d_cum(
            d_deforestation_surface_d_deforestation_surface)
        d_forest_surface_d_invest = self.forest_model.d_forestation_surface_d_invest()
        d_cun_forest_surface_d_invest = self.forest_model.d_cum(
            d_forest_surface_d_invest)

        # forest surface vs deforestation grad
        self.set_partial_derivative_for_other_types(
            (Forest.FOREST_SURFACE_DF, 'forest_surface_evol'), (
                Forest.DEFORESTATION_SURFACE, 'deforested_surface'),
            d_deforestation_surface_d_deforestation_surface)
#         self.set_partial_derivative_for_other_types(
#             (Forest.FOREST_SURFACE_DF,
#              'forest_surface_evol_cumulative'),
#             (Forest.DEFORESTATION_SURFACE, 'deforested_surface'),
#             d_cum_deforestation_d_deforestation_surface)

        # forest surface vs forest invest
        self.set_partial_derivative_for_other_types(
            (Forest.FOREST_SURFACE_DF, 'forest_surface_evol'), (
                Forest.REFORESTATION_INVESTMENT, 'forest_investment'),
            d_forest_surface_d_invest)
#         self.set_partial_derivative_for_other_types(
#             (Forest.FOREST_SURFACE_DF,
#              'forest_surface_evol_cumulative'),
#             (Forest.REFORESTATION_INVESTMENT, 'forest_investment'),
#             d_cun_forest_surface_d_invest)

        # d_CO2 d deforestation
        d_CO2_emitted_d_deforestation_surface = self.forest_model.d_CO2_emitted(
            d_deforestation_surface_d_deforestation_surface)
        d_cum_CO2_emitted_d_deforestation_surface = self.forest_model.d_cum(
            d_CO2_emitted_d_deforestation_surface)

        self.set_partial_derivative_for_other_types(
            (Forest.CO2_EMITTED_FOREST_DF, 'emitted_CO2_evol_cumulative'),
            (Forest.DEFORESTATION_SURFACE, 'deforested_surface'),
            d_cum_CO2_emitted_d_deforestation_surface)

        # d_CO2 d invest
        d_CO2_emitted_d_invest = self.forest_model.d_CO2_emitted(
            d_forest_surface_d_invest)
        d_cum_CO2_emitted_d_invest = self.forest_model.d_cum(
            d_CO2_emitted_d_invest)

        self.set_partial_derivative_for_other_types(
            (Forest.CO2_EMITTED_FOREST_DF, 'emitted_CO2_evol_cumulative'),
            (Forest.REFORESTATION_INVESTMENT, 'forest_investment'),
            d_cum_CO2_emitted_d_invest)

    def get_chart_filter_list(self):

        # For the outputs, making a graph for tco vs year for each range and for specific
        # value of ToT with a shift of five year between then

        chart_filters = []

        chart_list = [ForestDiscipline.FOREST_CHARTS]

        # First filter to deal with the view : program or actor
        chart_filters.append(ChartFilter(
            'Charts filter', chart_list, chart_list, 'charts'))

        return chart_filters

    def get_post_processing_list(self, chart_filters=None):
        '''
        For the outputs, making a graph for tco vs year for each range and for specific
        value of ToT with a shift of five year between then
        '''
        instanciated_charts = []

        # Overload default value with chart filter
        if chart_filters is not None:
            for chart_filter in chart_filters:
                if chart_filter.filter_key == 'charts':
                    chart_list = chart_filter.selected_values

        if ForestDiscipline.FOREST_CHARTS in chart_list:

            forest_surface_df = self.get_sosdisc_outputs(
                Forest.FOREST_DETAIL_SURFACE_DF)
            years = forest_surface_df['years'].values.tolist()
            # values are *1000 to convert from Gha to Mha
            surface_evol_by_year = forest_surface_df['forest_surface_evol'].values * 1000
            surface_evol_cum = forest_surface_df['forest_surface_evol_cumulative'].values * 1000
            deforested_surface_by_year = forest_surface_df['deforested_surface'].values * 1000
            deforested_surface_cum = forest_surface_df['deforested_surface_cumulative'].values * 1000
            forested_surface_by_year = forest_surface_df['forested_surface'].values * 1000
            forested_surface_cum = forest_surface_df['forested_surface_cumulative'].values * 1000

            # forest evolution year by year chart
            new_chart = TwoAxesInstanciatedChart('years', 'Forest surface evolution [Mha / year]',
                                                 chart_name='Forest surface evolution', stacked_bar=True)

            deforested_series = InstanciatedSeries(
                years, deforested_surface_by_year.tolist(), 'Deforestation', 'bar')
            forested_series = InstanciatedSeries(
                years, forested_surface_by_year.tolist(), 'Reforestation', 'bar')
            total_series = InstanciatedSeries(
                years, surface_evol_by_year.tolist(), 'Surface evolution', InstanciatedSeries.LINES_DISPLAY)

            new_chart.add_series(deforested_series)
            new_chart.add_series(total_series)
            new_chart.add_series(forested_series)

            instanciated_charts.append(new_chart)

            # forest cumulative evolution chart
            new_chart = TwoAxesInstanciatedChart('years', 'Cumulative forest surface evolution [Mha]',
                                                 chart_name='Cumulative forest surface evolution', stacked_bar=True)

            deforested_series = InstanciatedSeries(
                years, deforested_surface_cum.tolist(), 'Deforested surface', 'bar')
            forested_series = InstanciatedSeries(
                years, forested_surface_cum.tolist(), 'Forested surface', 'bar')
            total_series = InstanciatedSeries(
                years, surface_evol_cum.tolist(), 'Surface evolution', InstanciatedSeries.LINES_DISPLAY)
            new_chart.add_series(deforested_series)
            new_chart.add_series(total_series)
            new_chart.add_series(forested_series)

            instanciated_charts.append(new_chart)

            # CO2 graph

            CO2_emissions_df = self.get_sosdisc_outputs(
                'CO2_emissions_detail_df')
            CO2_emitted_year_by_year = CO2_emissions_df['emitted_CO2']
            CO2_captured_year_by_year = CO2_emissions_df['captured_CO2']
            CO2_total_year_by_year = CO2_emissions_df['emitted_CO2_evol']
            CO2_emitted_cum = CO2_emissions_df['emitted_CO2_cumulative']
            CO2_captured_cum = CO2_emissions_df['captured_CO2_cumulative']
            CO2_total_cum = CO2_emissions_df['emitted_CO2_evol_cumulative']

            new_chart = TwoAxesInstanciatedChart('years', 'CO2 emission & capture [GtCO2 / year]',
                                                 chart_name='Yearly forest delta CO2 emissions', stacked_bar=True)

            CO2_emitted_series = InstanciatedSeries(
                years, CO2_emitted_year_by_year.tolist(), 'CO2 emissions', InstanciatedSeries.BAR_DISPLAY)
            CO2_captured_series = InstanciatedSeries(
                years, CO2_captured_year_by_year.tolist(), 'CO2 capture', InstanciatedSeries.BAR_DISPLAY)
            CO2_total_series = InstanciatedSeries(
                years, CO2_total_year_by_year.tolist(), 'CO2 evolution', InstanciatedSeries.LINES_DISPLAY)

            new_chart.add_series(CO2_emitted_series)
            new_chart.add_series(CO2_total_series)
            new_chart.add_series(CO2_captured_series)
            instanciated_charts.append(new_chart)

            # in Gt
            new_chart = TwoAxesInstanciatedChart('years', 'CO2 emission & capture [GtCO2]',
                                                 chart_name='Forest CO2 emissions', stacked_bar=True)
            CO2_emitted_series = InstanciatedSeries(
                years, CO2_emitted_cum.tolist(), 'CO2 emissions', InstanciatedSeries.BAR_DISPLAY)
            CO2_captured_series = InstanciatedSeries(
                years, CO2_captured_cum.tolist(), 'CO2 capture', InstanciatedSeries.BAR_DISPLAY)
            CO2_total_series = InstanciatedSeries(
                years, CO2_total_cum.tolist(), 'CO2 evolution', InstanciatedSeries.LINES_DISPLAY, custom_data=['width'])

            new_chart.add_series(CO2_emitted_series)
            new_chart.add_series(CO2_total_series)
            new_chart.add_series(CO2_captured_series)
            instanciated_charts.append(new_chart)
        return instanciated_charts

max_stars_repo_name: os-climate/witness-core

max_issues_repo_name: os-climate/witness-core

max_forks_repo_name: os-climate/witness-core


text:
from typing import Dict

import numpy as np
import torch
from easypl.callbacks.mixers.base import MixBaseCallback

AVAILABLE_DOMENS = ['classification', 'detection', 'segmentation']

class Mixup(MixBaseCallback):
    """
    Callback for mixuping data operations

    Attributes
    ----------
    on_batch: bool
        If True generate samples from batch otherwise from dataset.

    p: float
        Mix probability.

    num_workers: int
        Number of workers for mixing operation.

    alpha: float
        Parameter of mixup operation.

    domen: str
        Name of task, in which will be mixed samples. Available: ["classification, segmentation"].
    """

    def __init__(
            self,
            on_batch: bool = True,
            alpha: float = 0.4,
            num_workers: int = 1,
            p: float = 0.5,
            domen: str = 'classification'
    ):
        super().__init__(on_batch=on_batch, samples_per=1, num_workers=num_workers, p=p)
        self.alpha = alpha
        if domen not in AVAILABLE_DOMENS:
            raise NotImplementedError(
                f'Domen {domen} is not supported in Mixup callback. Available domens: {AVAILABLE_DOMENS}')
        self.domen = domen

    def __mix_classificate(self, sample1: dict, sample2: dict, alpha: float) -> dict:
        return {
            key: sample1[key] * alpha + sample2[key] * (1 - alpha) if isinstance(sample1[key], torch.Tensor) else
            sample1[key] for key in sample1
        }

    def __mix_detection(self, sample1: dict, sample2: dict, alpha: float) -> dict:
        raise NotImplementedError

    def mix(
            self,
            sample1: Dict,
            sample2: Dict
    ) -> Dict:
        """
        Mixup method for two samples.

        Attributes
        ----------
        sample1: Dict
            Sample of batch, which will be sampled with sample from `sample2`.

        sample2: Dict
            Sample from batch or dataset.

        Returns
        -------
        Dict
            Mixed sample.
        """
        sample2 = {key: sample2[key][0] for key in sample2}
        alpha = np.random.beta(self.alpha, self.alpha)
        if self.domen == 'classification':
            return self.__mix_classificate(sample1, sample2, alpha)
        elif self.domen == 'detection':
            return self.__mix_detection(sample1, sample2, alpha)
        elif self.domen == 'segmentation':
            return self.__mix_classificate(sample1, sample2, alpha)
        else:
            return sample1


text:
#!/usr/bin/env python
"""Tests of the geometry package."""
from math import sqrt

from numpy import (
    all,
    allclose,
    arange,
    array,
    insert,
    isclose,
    mean,
    ones,
    sum,
    take,
)
from numpy.linalg import inv, norm
from numpy.random import choice, dirichlet
from numpy.testing import assert_allclose

from cogent3.maths.geometry import (
    aitchison_distance,
    alr,
    alr_inv,
    center_of_mass,
    center_of_mass_one_array,
    center_of_mass_two_array,
    clr,
    clr_inv,
    distance,
    multiplicative_replacement,
    sphere_points,
)
from cogent3.util.unit_test import TestCase, main

__author__ = "Sandra Smit"
__copyright__ = "Copyright 2007-2020, The Cogent Project"
__credits__ = ["Sandra Smit", "Rob Knight", "Helmut Simon"]
__license__ = "BSD-3"
__version__ = "2020.6.30a"
__maintainer__ = "Sandra Smit"
__email__ = "sandra.smit@colorado.edu"
__status__ = "Production"

class CenterOfMassTests(TestCase):
    """Tests for the center of mass functions"""

    def setUp(self):
        """setUp for all CenterOfMass tests"""
        self.simple = array([[1, 1, 1], [3, 1, 1], [2, 3, 2]])
        self.simple_list = [[1, 1, 1], [3, 1, 1], [2, 3, 2]]
        self.more_weight = array([[1, 1, 3], [3, 1, 3], [2, 3, 50]])
        self.square = array([[1, 1, 25], [3, 1, 25], [3, 3, 25], [1, 3, 25]])
        self.square_odd = array([[1, 1, 25], [3, 1, 4], [3, 3, 25], [1, 3, 4]])
        self.sec_weight = array([[1, 25, 1], [3, 25, 1], [3, 25, 3], [1, 25, 3]])

    def test_center_of_mass_one_array(self):
        """center_of_mass_one_array should behave correctly"""
        com1 = center_of_mass_one_array
        self.assertEqual(com1(self.simple), array([2, 2]))
        self.assertEqual(com1(self.simple_list), array([2, 2]))
        self.assertFloatEqual(com1(self.more_weight), array([2, 2.785714]))
        self.assertEqual(com1(self.square), array([2, 2]))
        self.assertEqual(com1(self.square_odd), array([2, 2]))
        self.assertEqual(com1(self.sec_weight, 1), array([2, 2]))

    def test_CoM_one_array_wrong(self):
        """center_of_mass_one_array should fail on wrong input"""
        com1 = center_of_mass_one_array
        self.assertRaises(TypeError, com1, self.simple, "a")  # weight_idx wrong
        self.assertRaises(IndexError, com1, self.simple, 100)  # w_idx out of range
        # shape[1] out of range
        self.assertRaises(IndexError, com1, [1, 2, 3], 2)

    def test_center_of_mass_two_array(self):
        """center_of_mass_two_array should behave correctly"""
        com2 = center_of_mass_two_array
        coor = take(self.square_odd, (0, 1), 1)
        weights = take(self.square_odd, (2,), 1)
        self.assertEqual(com2(coor, weights), array([2, 2]))
        weights = weights.ravel()
        self.assertEqual(com2(coor, weights), array([2, 2]))

    def test_CoM_two_array_wrong(self):
        """center_of_mass_two_array should fail on wrong input"""
        com2 = center_of_mass_two_array
        weights = [1, 2]
        self.assertRaises(TypeError, com2, self.simple, "a")  # weight_idx wrong
        self.assertRaises(ValueError, com2, self.simple, weights)  # not aligned

    def test_center_of_mass(self):
        """center_of_mass should make right choice between functional methods
        """
        com = center_of_mass
        com1 = center_of_mass_one_array
        com2 = center_of_mass_two_array
        self.assertEqual(com(self.simple), com1(self.simple))
        self.assertFloatEqual(com(self.more_weight), com1(self.more_weight))
        self.assertEqual(com(self.sec_weight, 1), com1(self.sec_weight, 1))
        coor = take(self.square_odd, (0, 1), 1)
        weights = take(self.square_odd, (2,), 1)
        self.assertEqual(com(coor, weights), com2(coor, weights))
        weights = weights.ravel()
        self.assertEqual(com(coor, weights), com2(coor, weights))

    def test_distance(self):
        """distance should return Euclidean distance correctly."""
        # for single dimension, should return difference
        a1 = array([3])
        a2 = array([-1])
        self.assertEqual(distance(a1, a2), 4)
        # for two dimensions, should work e.g. for 3, 4, 5 triangle
        a1 = array([0, 0])
        a2 = array([3, 4])
        self.assertEqual(distance(a1, a2), 5)
        # vector should be the same as itself for any dimensions
        a1 = array([1.3, 23, 5.4, 2.6, -1.2])
        self.assertEqual(distance(a1, a1), 0)
        # should match hand-calculated case for an array
        a1 = array([[1, -2], [3, 4]])
        a2 = array([[1, 0], [-1, 2.5]])
        self.assertEqual(distance(a1, a1), 0)
        self.assertEqual(distance(a2, a2), 0)
        self.assertEqual(distance(a1, a2), distance(a2, a1))
        self.assertFloatEqual(distance(a1, a2), sqrt(22.25))

    def test_sphere_points(self):
        """tests sphere points"""
        self.assertEqual(sphere_points(1), array([[1.0, 0.0, 0.0]]))

class TestAitchison(TestCase):
    def setUp(self):
        x = choice(20, size=10) + 0.1
        self.x = x
        a = arange(1, 7)
        self.a = a
        d = dirichlet(a, size=2)
        self.d = d

    def test_Aitchison_transforms(self):
        """Test that alr_inv of alr is in fact the inverse
        of alr. Ditto for clr_inv and clr. Then test that clr
        transforms into hyperplane x1 + ... + xn=0."""
        length = len(self.x)
        for col in range(-1, length):
            y = alr_inv(self.x, col)
            assert allclose(self.x, alr(y, col)), (
                "Failed alr inverse test for col = " + str(col) + "."
            )

        z = dirichlet(self.x)
        y = clr(z)
        assert allclose(z, clr_inv(y)), "Failed clr inverse test."
        assert allclose(sum(y), 0), "Failed clr hyperplane test."

    def test_Aitchison_distance(self):
        x = self.d[0]
        y = self.d[1]
        assert allclose(
            aitchison_distance(x, y), norm(clr(x) - clr(y))
        ), "Failed distance test."

    def test_multiplicative_replacement(self):
        x1 = dirichlet(self.a)
        y1 = insert(x1, 3, 0)
        u = multiplicative_replacement(y1)
        assert allclose(
            y1, u, atol=1e-2
        ), "Multiplicative replacement peturbation is too large."
        assert isclose(
            sum(u), 1
        ), "Multiplicative replacement does not yield a composition."

if __name__ == "__main__":
    main()

reason: from numpy


text:
"""Some special pupropse layers for SSD."""

import keras.backend as K
from keras.engine.topology import InputSpec
from keras.engine.topology import Layer
import numpy as np
import tensorflow as tf

class Normalize(Layer):
    """Normalization layer as described in ParseNet paper.

    # Arguments
        scale: Default feature scale.

    # Input shape
        4D tensor with shape:
        `(samples, channels, rows, cols)` if dim_ordering='th'
        or 4D tensor with shape:
        `(samples, rows, cols, channels)` if dim_ordering='tf'.

    # Output shape
        Same as input

    # References
        http://cs.unc.edu/~wliu/papers/parsenet.pdf

    #TODO
        Add possibility to have one scale for all features.
    """
    def __init__(self, scale, **kwargs):
        if K.image_data_format() == 'channels_last':
            self.axis = 3
        else:
            self.axis = 1
        self.scale = scale
        super(Normalize, self).__init__(**kwargs)

    def build(self, input_shape):
        self.input_spec = [InputSpec(shape=input_shape)]
        shape = (input_shape[self.axis],)
        init_gamma = self.scale * np.ones(shape)
        self.gamma = K.variable(init_gamma, name='{}_gamma'.format(self.name))
        self.trainable_weights = [self.gamma]

    def call(self, x, mask=None):
        output = K.l2_normalize(x, self.axis)
        output *= self.gamma
        return output

class PriorBox(Layer):
    """Generate the prior boxes of designated sizes and aspect ratios.

    # Arguments
        img_size: Size of the input image as tuple (w, h).
        min_size: Minimum box size in pixels.
        max_size: Maximum box size in pixels.
        aspect_ratios: List of aspect ratios of boxes.
        flip: Whether to consider reverse aspect ratios.
        variances: List of variances for x, y, w, h.
        clip: Whether to clip the prior's coordinates
            such that they are within [0, 1].

    # Input shape
        4D tensor with shape:
        `(samples, channels, rows, cols)` if dim_ordering='th'
        or 4D tensor with shape:
        `(samples, rows, cols, channels)` if dim_ordering='tf'.

    # Output shape
        3D tensor with shape:
        (samples, num_boxes, 8)

    # References
        https://arxiv.org/abs/1512.02325

    #TODO
        Add possibility not to have variances.
        Add Theano support
    """
    def __init__(self, img_size, min_size, max_size=None, aspect_ratios=None,
                 flip=True, variances=[0.1], clip=True, **kwargs):
        if K.image_data_format() == 'channels_last':
            self.waxis = 2
            self.haxis = 1
        else:
            self.waxis = 3
            self.haxis = 2
        self.img_size = img_size
        if min_size <= 0:
            raise Exception('min_size must be positive.')
        self.min_size = min_size
        self.max_size = max_size
        self.aspect_ratios = [1.0]
        if max_size:
            if max_size < min_size:
                raise Exception('max_size must be greater than min_size.')
            self.aspect_ratios.append(1.0)
        if aspect_ratios:
            for ar in aspect_ratios:
                if ar in self.aspect_ratios:
                    continue
                self.aspect_ratios.append(ar)
                if flip:
                    self.aspect_ratios.append(1.0 / ar)
        self.variances = np.array(variances)
        self.clip = True
        super(PriorBox, self).__init__(**kwargs)

    def compute_output_shape(self, input_shape):
        num_priors_ = len(self.aspect_ratios)
        layer_width = input_shape[self.waxis]
        layer_height = input_shape[self.haxis]
        num_boxes = num_priors_ * layer_width * layer_height
        return (input_shape[0], num_boxes, 8)

    def call(self, x, mask=None):
        if hasattr(x, '_keras_shape'):
            input_shape = x._keras_shape
        elif hasattr(K, 'int_shape'):
            input_shape = K.int_shape(x)
        layer_width = input_shape[self.waxis]
        layer_height = input_shape[self.haxis]
        img_width = self.img_size[0]
        img_height = self.img_size[1]
        # define prior boxes shapes
        box_widths = []
        box_heights = []
        for ar in self.aspect_ratios:
            if ar == 1 and len(box_widths) == 0:
                box_widths.append(self.min_size)
                box_heights.append(self.min_size)
            elif ar == 1 and len(box_widths) > 0:
                box_widths.append(np.sqrt(self.min_size * self.max_size))
                box_heights.append(np.sqrt(self.min_size * self.max_size))
            elif ar != 1:
                box_widths.append(self.min_size * np.sqrt(ar))
                box_heights.append(self.min_size / np.sqrt(ar))
        box_widths = 0.5 * np.array(box_widths)
        box_heights = 0.5 * np.array(box_heights)
        # define centers of prior boxes
        step_x = img_width / layer_width
        step_y = img_height / layer_height
        linx = np.linspace(0.5 * step_x, img_width - 0.5 * step_x,
                           layer_width)
        liny = np.linspace(0.5 * step_y, img_height - 0.5 * step_y,
                           layer_height)
        centers_x, centers_y = np.meshgrid(linx, liny)
        centers_x = centers_x.reshape(-1, 1)
        centers_y = centers_y.reshape(-1, 1)
        # define xmin, ymin, xmax, ymax of prior boxes
        num_priors_ = len(self.aspect_ratios)
        prior_boxes = np.concatenate((centers_x, centers_y), axis=1)
        prior_boxes = np.tile(prior_boxes, (1, 2 * num_priors_))
        prior_boxes[:, ::4] -= box_widths
        prior_boxes[:, 1::4] -= box_heights
        prior_boxes[:, 2::4] += box_widths
        prior_boxes[:, 3::4] += box_heights
        prior_boxes[:, ::2] /= img_width
        prior_boxes[:, 1::2] /= img_height
        prior_boxes = prior_boxes.reshape(-1, 4)
        if self.clip:
            prior_boxes = np.minimum(np.maximum(prior_boxes, 0.0), 1.0)
        # define variances
        num_boxes = len(prior_boxes)
        if len(self.variances) == 1:
            variances = np.ones((num_boxes, 4)) * self.variances[0]
        elif len(self.variances) == 4:
            variances = np.tile(self.variances, (num_boxes, 1))
        else:
            raise Exception('Must provide one or four variances.')
        prior_boxes = np.concatenate((prior_boxes, variances), axis=1)
        prior_boxes_tensor = K.expand_dims(K.variable(prior_boxes), 0)
        if K.backend() == 'tensorflow':
            pattern = [tf.shape(x)[0], 1, 1]
            prior_boxes_tensor = tf.tile(prior_boxes_tensor, pattern)
        elif K.backend() == 'theano':
            #TODO
            pass
        return prior_boxes_tensor


text:
import numpy as np
np.random.seed(0)
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns

def array_info(array: np.ndarray) -> None:
    print(f"ndim: {array.ndim}")
    print(f"shape: {array.shape}")
    print(f"size: {array.size}")
    print(f"dtype: {array.dtype}")
    print(f"values:\n{array}\n")

print("=1=")
mean = [-5, 5]
cov = [[1, 0], [0, 1]]
n_samples = 100
x1, x2 = np.random.multivariate_normal(mean, cov, n_samples).T

y = x1 + x2

print("=2=")
# plt.scatter(x1, x2, c=y, cmap="coolwarm")
# plt.xlabel("x_1")
# plt.ylabel("x_2")
# plt.colorbar()
# plt.show()

print("=3=")
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import FormatStrFormatter
from matplotlib.ticker import LinearLocator

fig = plt.figure(figsize=(8, 8))
ax = Axes3D(fig, auto_add_to_figure=False)
fig.add_axes(ax)

fig3 = ax.scatter3D(x1, x2, y, c=y, cmap='coolwarm')
ax.set_xlabel("x_1")
ax.set_ylabel("x_2")
fig.colorbar(fig3)
plt.show()


text:
#!/usr/local/bin/python
import os
import json
import re
import cPickle

import numpy as np
import utils as tech

import matplotlib.pyplot as plt

from matplotlib import rcParams

rcParams['text.usetex'] = True

READ = 'r'
WRITE = 'w'
base = os.path.join(os.getcwd(),'data')
sources = ['harrison','wikipedia','braunwald','twitter']
keywords = [word.replace(' ','_') for word in open('keywords').read().splitlines()]
CORPUS_FILENAME = 'corpus.json'

'''
    Data files structured as:

     	./data/source/keyword
'''

#---LOAD DATA
if not os.path.isfile(CORPUS_FILENAME):
	#More expressive than itertools.product, small loops --> no important speed or memory difference
	corpus = {}
	for source in sources:
		corpus[source] = {}
		for disease in keywords:
				path = os.path.join(base,source,disease)
				text = ' '.join(open(os.path.join(path,filename),READ).read() for filename in os.listdir(path)
										if not os.path.isdir(os.path.join(path,filename)))
				text = text.replace('.',' ').replace("\n"," ")
				text = re.sub(r"[^\x00-\x7F]","",text) #Regexp faster than iterating through string to remove non-ASCII
				corpus[source][disease]  = list(tech.cleanse(text))
				#Cleanse returns type set. Type set is not JSON serializable. Type list is.
	json.dump(corpus,open(CORPUS_FILENAME,WRITE))

else:
	corpus = json.load(open(CORPUS_FILENAME,READ))

#--- CALCULATE JACCARD SIMILARITY

source_rubric = [[source for source in sources]
						 for source in sources]

filenames = ['jaccard-similarity-%s'%disease for disease in keywords]
filenames += ['jaccard-similarities.json']

if not all([os.path.isfile(filename) for filename in filenames]):
	jaccard_matrices = {disease:np.zeros((len(sources),len(sources))) for disease in keywords}
	for disease in keywords:
		jaccard_matrices[disease] = np.array([[tech.jaccard(corpus[sources[i]][disease],corpus[sources[j]][disease])
											for i in xrange(len(sources))]
											for j in xrange(len(sources))])

		fig = plt.figure()
		ax = fig.add_subplot(111)
		cax = ax.imshow(jaccard_matrices[disease],interpolation='nearest',aspect='equal',vmin=0,vmax=1)

		ax.set_xticks(range(len(sources)))
		ax.set_yticks(range(len(sources)))

		ax.set_xticklabels(map(tech.format,sources))
		ax.set_yticklabels(map(tech.format,sources))
		cbar = plt.colorbar(cax)
		cbar.set_label(tech.format('Jaccard Similarity'))
		fig.tight_layout()
		plt.savefig('jaccard-similarity-%s-w-twitter'%disease)

	cPickle.dump(jaccard_matrices,open('jaccard-similarities.json',WRITE))

#--- BOOTSTRAPPING
lens = [len(corpus[source][disease]) for source in sources] #Does order matter?
amalgamated_corpus = ' '.join(' '.join(corpus[source][disease]) for disease in keywords for source in sources)
#N.B. Don't depucliated -- must preserve original word frequencies for resampling
jaccard_distributions = tech.resample(amalgamated_corpus,n_partitions=len(lens),partition_sizes=lens,repetitions=10000,
	monitor=True,save=True)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(jaccard_distributions,color='k')
tech.adjust_spines(ax)
ax.set_xlabel(tech.format('Jaccard Similarity'))
ax.set_ylabel(tech.format('No. of occurence'))
plt.tight_layout()
plt.savefig('distribution-jaccard-similarities-w-twitter.tiff')


text:
# -*- coding: utf-8 -*-

import numpy as _np
from scipy.special import hermite, genlaguerre
from scipy.interpolate import RectBivariateSpline

USE_CV2 = False

if USE_CV2:
    import cv2

USE_SKIMAGE = False
if USE_SKIMAGE:
    from skimage.restoration import unwrap_phase as _unwrap_phase
else:
    #used in PhaseUnwrap
    # own implementation currently slower, but seems a little more stable
    # with jumpy phases and of course removes dependency on the extra package
    from .unwrap import unwrap_phase as _unwrap_phase

from .units import deg
from .field import Field
from .subs import Inv_Squares

def BeamMix(Fin1, Fin2):
    """
    Fout = BeamMix(F1, F2)

    :ref:`Addition of the fields F1 and F2. <BeamMix>`

    Args::

        F1, F2: input fields

    Returns::

        Fout: output field (N x N square array of complex numbers).

    Example:

    :ref:`Two holes interferometer <Young>`

    """
    if Fin1.field.shape != Fin2.field.shape:
        raise ValueError('Field sizes do not match')
    Fout = Field.copy(Fin1)
    Fout.field += Fin2.field
    return Fout

def CircAperture(R, x_shift, y_shift, Fin):
    """
    Fout = CircAperture(R, x_shift, y_shift, Fin)

    :ref:`Propagates the field through a circular aperture. <CircAperture>`

    Args::

        R: radius of the aperture
        x_shift, y_shift: shift from the center
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    Example:

    :ref:`Diffraction from a circular aperture <circ_aperture>`

    """

    #from
    #https://stackoverflow.com/questions/44865023/
    # circular-masking-an-image-in-python-using-numpy-arrays
    Fout = Field.copy(Fin)

    Y, X = Fout.mgrid_cartesian
    Y = Y - y_shift
    X = X - x_shift

    dist_sq = X**2 + Y**2 #squared, no need for sqrt

    Fout.field[dist_sq > R**2] = 0.0
    return Fout

def CircScreen(R, x_shift, y_shift, Fin):
    """
    Fout = CircScreen(R, x_shift, y_shift, Fin)

    :ref:`Diffracts the field by a circular screen. <CircScreen>`

    Args::

        R: radius of the screen
        x_shift, y_shift: shift from the center
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    Example:

    :ref:`Spot of Poisson <Poisson>`

    """

    #from
    #https://stackoverflow.com/questions/44865023/
    # circular-masking-an-image-in-python-using-numpy-arrays
    Fout = Field.copy(Fin)

    Y, X = Fout.mgrid_cartesian
    Y = Y - y_shift
    X = X - x_shift
    dist_sq = X**2 + Y**2 #squared, no need for sqrt

    Fout.field[dist_sq <= R**2] = 0.0
    return Fout

def GaussAperture(w, x_shift, y_shift, T, Fin):
    """
    Fout = GaussAperture(w, x_shift, y_shift, T, Fin)

    :ref:`Inserts an aperture with a Gaussian shape in the field. <GaussAperture>`

        :math:`F_{out}(x,y)= \\sqrt{T}e^{ -\\frac{ x^{2}+y^{2} }{2w^{2}} } F_{in}(x,y)`

    Args::

        w: 1/e intensity width
        x_shift, y_shift: shift from center
        T: center intensity transmission
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """

    Fout = Field.copy(Fin)

    Y, X = Fout.mgrid_cartesian
    Y = Y - y_shift
    X = X - x_shift

    w2=w*w*2
    SqrtT=_np.sqrt(T)
    Fout.field*=SqrtT*_np.exp(-(X*X+Y*Y)/w2)
    return Fout

def SuperGaussAperture(w, x_shift, y_shift, T, n, Fin):
    """
    Fout = SuperGaussAperture(w, x_shift, y_shift, T, n, Fin)

    :ref:`Inserts an aperture with a Gaussian shape in the field. <GaussAperture>`

        :math:`F_{out}(x,y)= \\sqrt{T}e^{ -\\left [ \\frac{ x^{2}+y^{2} }{2w^{2}} \\right ]^n } F_{in}(x,y)`

    Args::

        w: 1/e intensity width
        x_shift, y_shift: shift from center
        T: center intensity transmission
        n: power of the super Gaussian
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """

    Fout = Field.copy(Fin)

    Y, X = Fout.mgrid_cartesian
    Y = Y - y_shift
    X = X - x_shift

    w2=w*w*2
    SqrtT=_np.sqrt(T)
    Fout.field*=SqrtT*_np.exp(-((X*X+Y*Y)/w2)**n)
    return Fout

def GaussScreen(w, x_shift, y_shift, T, Fin):
    """
    Fout = GaussScreen(w, x_shift, y_shift, T, Fin)

    :ref:`Inserts a screen with a Gaussian shape in the field. <GaussScreen>`

        :math:`F_{out}(x,y)= \\sqrt{1-(1-T)e^{ -\\frac{ x^{2}+y^{2} }{w^{2}} }} F_{in}(x,y)`

    Args::

        w: 1/e intensity width
        x_shift, y_shift: shift from center
        T: center intensity transmission
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Fout = Field.copy(Fin)

    Y, X = Fout.mgrid_cartesian
    Y = Y - y_shift
    X = X - x_shift

    w2=w*w
    Fout.field*=_np.sqrt(1-(1-T)*_np.exp(-(X*X+Y*Y)/w2))
    return Fout

def GaussHermite( n, m, A, w0, Fin):
    """
    Fout = GaussHermite(m, n, A, w0, Fin)

    :ref:`Substitutes a Hermite-Gauss mode (beam waist) in the field. <GaussHermite>`

        :math:`F_{m,n}(x,y,z=0) = A H_m\\left(\\dfrac{\\sqrt{2}x}{w_0}\\right)H_n\\left(\\dfrac{\\sqrt{2}y}{w_0}\\right)e^{-\\frac{x^2+y^2}{w_0^2}}`

    Args::

        m, n: mode indices
        A: Amplitude
        w0: Guaussian spot size parameter in the beam waist (1/e amplitude point)
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    Reference::

        A. Siegman, "Lasers", p. 642

    """

    Fout = Field.copy(Fin)

    Y, X = Fout.mgrid_cartesian
    #Y = Y - y_shift
    #X = X - x_shift

    sqrt2w0=_np.sqrt(2.0)/w0
    w02=w0*w0

    Fout.field  = A * hermite(m)(sqrt2w0*X)*hermite(n)(sqrt2w0*Y)*_np.exp(-(X*X+Y*Y)/w02)
    return Fout

def GaussLaguerre(p, l, A, w0, Fin):
    """
    Fout = GaussLaguerre(p, l, A, w0, Fin)

    :ref:`Substitutes a Laguerre-Gauss mode (beam waist) in the field. <GaussLaguerre>`

        :math:`F_{p,l}(x,y,z=0) = A \\left(\\frac{\\rho}{2}\\right)^{\\frac{|l|}{2} }L^p_l\\left(\\rho\\right)e^{-\\frac{\\rho}{2}}\\cos(l\\theta)`,

        with :math:`\\rho=\\frac{2(x^2+y^2)}{w_0^2}`

    Args::

        p, l: mode indices
        A: Amplitude
        w0: Guaussian spot size parameter in the beam waist (1/e amplitude point)
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    Reference::

        A. Siegman, "Lasers", p. 642

    """

    Fout = Field.copy(Fin)
    R, Phi = Fout.mgrid_polar
    w02=w0*w0
    la=abs(l)
    rho = 2*R*R/w02
    Fout.field = A * rho**(la/2) * genlaguerre(p,la)(rho) * _np.exp(-rho/2) * _np.cos(l*Phi)
    return Fout

def IntAttenuator(att, Fin):
    """
    Fout = IntAttenuator(att, Fin)

    :ref:`Attenuates the intensity of the field. <IntAttenuator>`

        :math:`F_{out}(x,y)=\\sqrt{att}F_{in}(x,y)`

    Args::

        att: intensity attenuation factor
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Efactor = _np.sqrt(att) #att. given as intensity
    Fout = Field.copy(Fin)
    Fout.field *= Efactor
    return Fout

def Intensity(flag, Fin):
    """
    I=Intensity(flag,Fin)

    :ref:`Calculates the intensity of the field. <Intensity>`

    :math:`I(x,y)=F_{in}(x,y).F_{in}(x,y)^*`

    Args::

        flag: 0= no normalization, 1=normalized to 1, 2=normalized to 255 (for bitmaps)
        Fin: input field

    Returns::

        I: intensity distribution (N x N square array of doubles)

    """
    I = _np.abs(Fin.field)**2
    if flag > 0:
        Imax = I.max()
        if Imax == 0.0:
            raise ValueError('Cannot normalize because of 0 beam power.')
        I = I/Imax
        if flag == 2:
            I = I*255
    return I

def Interpol(new_size, new_N, x_shift, y_shift, angle, magnif, Fin):
    """
    Fout = Interpol(NewSize, NewN, x_shift, y_shift, angle, magnif, Fin)

    :ref:`Interpolates the field to a new grid size, grid dimension. <Interpol>`

    Args::

        NewSize: the new grid size
        NewN: the new grid dimension
        x_shift, y_shift: shift of the field
        angle: rotation of the field in degrees
        magnif: magnification of the field amplitude
        Fin: input field

    Returns::

        Fout: output field (Nnew x Nnew square array of complex numbers).

    """

    Fout = Field.begin(new_size, Fin.lam, new_N)
    Fout.field[:,:] = 0.0

    legacy = True
    if legacy:
        Pi = 3.141592654 #compare Cpp results numerically
    else:
        Pi = _np.pi #more accurate, but slightly different results
    angle *= Pi/180.
    cc=_np.cos(angle)
    ss=_np.sin(angle)

    if legacy:
        #dx defined differently
        size_old = Fin.siz
        old_number = Fin.N
        dx_old = size_old/(old_number-1)
        on21 = int(old_number/2)
        Xold = dx_old * _np.arange(-on21, old_number-on21)
        Yold = dx_old * _np.arange(-on21, old_number-on21)
    else:
        Xold = Fin.xvalues
        Yold = Fin.yvalues

    if legacy:
        dx_new = new_size/(new_N-1) #TODO legacy, once again without -1 seems correct
        nn21 = int(new_N/2)
        X0 = dx_new * _np.arange(-nn21, new_N-nn21)
        Y0 = dx_new * _np.arange(-nn21, new_N-nn21)
        X0, Y0 = _np.meshgrid(X0, Y0)
    else:
        dx_new = Fout.dx
        Y0, X0 = Fout.mgrid_cartesian #note swapped order!
    X0 -= x_shift
    Y0 -= y_shift
    Xnew = (X0*cc + Y0*ss)/magnif
    Ynew = (X0*(-ss) + Y0* cc)/magnif

    xmin, xmax = Xold[0], Xold[-1]
    ymin, ymax = Yold[0], Yold[-1]
    #filter strictly inside (not <=) since edge pixels seem wrong in interp
    filtmask = ((Xnew > xmin) & (Xnew < xmax) &
                (Ynew > ymin) & (Ynew < ymax))
    # same goes for Cpp lightpipes, interpolating a 20x20 grid to a 20x20 grid
    # of same size will have 0s along the edges and only 18x18 useful pixels

    #instead of calling interp for all pixels, only call for those new pixels
    # who's coordinates (transformed to old) are inside old grid box
    Xmask = Xnew[filtmask] #flat list of X-values, not meshgrid anymore
    Ymask = Ynew[filtmask]

    use_scipy_interp = False
    if use_scipy_interp:
        ks = 1 #spline order: linear or higher
        interp_real = RectBivariateSpline(Xold, Yold, Fin.field.real,
                                          kx=ks, ky=ks)
        interp_imag = RectBivariateSpline(Xold, Yold, Fin.field.imag,
                                          kx=ks, ky=ks)

        out_real = interp_real(Xmask, Ymask, grid=False)
        out_imag = interp_imag(Xmask, Ymask, grid=False)
        out_comp = out_real + 1j* out_imag
        Fout.field[filtmask] = out_comp
    else:
        out_z = Inv_Squares(Xmask, Ymask, Fin.field, dx_old)
        Fout.field[filtmask] = out_z
    Fout.field /= magnif
    return Fout

def MultIntensity(Intens, Fin):
    """
    Fout = MultIntensity(Intens, Fin)

    :ref:`Multiplies the field with a given intensity distribution. <MultIntensity>`

    Args::

        Intens: N x N square array of real numbers or scalar
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    if not _np.isscalar(Intens):
        if Intens.shape != Fin.field.shape:
            raise ValueError('Intensity pattern shape does not match field size')
    Fout = Field.copy(Fin)
    Efield = _np.sqrt(Intens)
    Fout.field *= Efield
    return Fout

def MultPhase(Phi, Fin):
    """
    Fout = MultPhase(Phase, Fin)

    :ref:`Multiplies the field with a given phase distribution. <MultPhase>`

    Args::

        Phase: N x N square array of real numbers or scalar
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    if not _np.isscalar(Phi):
        if Phi.shape != Fin.field.shape:
            raise ValueError('Phase pattern shape does not match field size')
    Fout = Field.copy(Fin)
    Fout.field *= _np.exp(1j*Phi)
    return Fout

def Normal(Fin):
    """
    Fout = Normal(Fin)

    :ref:`Normalizes the field using beam power. <Normal>`

        :math:`F_{out}(x,y)= \\frac{F_{in}(x,y)}{\\sqrt{P}}`

        with

        :math:`P=\\int \\int F_{in}(x,y)^2 dx dy`

    Args::

        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Fabs = _np.abs(Fin.field)**2
    Fabs *= Fin.dx**2
    Ptot = Fabs.sum()
    if Ptot == 0.0:
        raise ValueError('Error in Normal(Fin): Zero beam power!')
    Fout = Field.copy(Fin)
    Fout.field *= _np.sqrt(1/Ptot)
    return Fout

def Phase(Fin, unwrap = False, units='rad', blank_eps=0):
    """
    Phi=Phase(Fin)

    :ref:`Calculates the phase of the field. <Phase>`

    Args::

        Fin: input field
        unwrap: Call PhaseUnwarp on the extracted Phase, Default is False
        units: 'opd': returned in [meters] of optical path length
                'lam': returned in multiples of lambda
                'rad': returned in multiples of 2pi phase jumps (default)
        blank_eps: [fraction] of max. Intensity at which to blank the phase
            and replace the value with numpy.nan (e.g. 1e-3==0.1%)
            Set to 0 or None to disable

    Returns::

        Phi: phase distribution (N x N square array of doubles)

    """
    _2pi = 2*_np.pi
    Phi = _np.angle(Fin.field)
    if unwrap:
        Phi = PhaseUnwrap(Phi)

    if units=='opd':
        Phi = Phi/_2pi*Fin.lam #a PtV of 2pi will yield e.g. 1*lam=1e-6=1um
    elif units=='lam':
        Phi = Phi/_2pi #a PtV of 2pi=6.28 will yield 1 (as in 1 lambda)
    elif units=='rad':
        pass #a PtV of 2pi will yield 6.28 as requested
    else:
        raise ValueError('Unknown value for option units={}'.format(units))

    if blank_eps:
        I = Intensity(0,Fin)
        Phi[I<blank_eps*I.max()] = _np.nan

    return Phi

def PhaseSpiral(Fin, m = 1):
    """
    Fout = PhaseSpiral(Fin, m = 1)
    Multiplies Fin with a spiral phase distribution.

    Args::

        required:
        Fin: input field

        optional:
        m = 1: order of the spiral.
    """
    Fout = Field.copy(Fin)
    R, Phi = Fout.mgrid_polar
    Fout.field *= _np.exp(1j * m * Phi)
    return Fout

def PhaseUnwrap(Phi):
    """
    PhiOut=PhaseUnwrap(PhiIn)

    :ref:`Unwraps (removes jumps of pi radians) the phase. <PhaseUnwrap>`

    Args::

        PhiIn: input phase distribution

    Returns::

        PhiOut: unwrapped phase distribution (N x N square array of doubles)

    """
    PhiU = _unwrap_phase(Phi)
    return PhiU

def Power(Fin):
    """
    P = Power(Fin)

    :ref:`Calculates the total power. <Power>`

    Args::

        Fin: input field

    Returns::

        P: output power (real number).

    """
    #TODO why does Normal() also sum dx**2 (==integral) while this does not??
    I = _np.abs(Fin.field)**2
    return I.sum()

def RandomIntensity(seed, noise, Fin):
    """
    Fout = RandomIntensity(seed, noise, Fin)

    :ref:`Adds random intensity to the field <RandomIntensity>`

    Args::

        seed: seed number for the random noise generator
        noise: level of the noise
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    #TODO implementation error in original LP: field error, not I error!
    # need to sqrt for that
    Fout = Field.copy(Fin)
    _np.random.seed(int(seed))
    N = Fout.N
    ranint = _np.random.rand(N, N)*noise
    Fout.field += ranint
    return Fout

def RandomPhase(seed, maxPhase, Fin):
    """
    Fout = RandomPhase(seed, maxPhase, Fin)

    :ref:`Adds random phase to the field <RandomPhase>`

    Args::

        seed: seed number for the random noise generator
        maxPhase: maximum phase in radians
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    #2020023 - ldo - tested similar result as Cpp version, although not
    # 1:1 since seed is different in numpy
    Fout = Field.copy(Fin)
    _np.random.seed(int(seed))
    N = Fout.N
    ranphase = (_np.random.rand(N, N)-0.5)*maxPhase
    Fout.field *= _np.exp(1j * ranphase)
    return Fout

def RectAperture(sx, sy, x_shift, y_shift, angle, Fin):
    """
    Fout = RectAperture(w, h, x_shift, y_shift, angle, Fin)

    :ref:`Propagates the field through a rectangular aperture. <RectAperture>`

    Args::

        w: width of the aperture
        h: height of the aperture
        x_shift, y_shift: shift from the center
        angle: rotation angle in degrees
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Fout = Field.copy(Fin)
    yy, xx = Fout.mgrid_cartesian
    yy = yy - y_shift
    xx = xx - x_shift
    if angle!=0.0:
        ang_rad = -1*angle*deg #-1 copied from Cpp convention
        cc = _np.cos(ang_rad)
        ss = _np.sin(ang_rad)
        xxr = cc * xx + ss * yy
        yyr = -ss * xx + cc * yy
        yy, xx = yyr, xxr
    matchx = _np.abs(xx) > sx/2
    matchy = _np.abs(yy) > sy/2
    Fout.field[matchx | matchy] = 0.0
    return Fout

def RectScreen(sx, sy, x_shift, y_shift, angle, Fin):
    """
    Fout = RectScreen(w, h, x_shift, y_shift, angle, Fin)

    :ref:`Diffracts the field by a rectangular screen. <RectScreen>`

    Args::

        w: width of the screen
        h: height of the screen
        x_shift, y_shift: shift from the center
        angle: rotation angle in degrees
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Fout = Field.copy(Fin)
    yy, xx = Fout.mgrid_cartesian
    yy = yy - y_shift
    xx = xx - x_shift
    if angle!=0.0:
        ang_rad = -1*angle*deg #-1 copied from Cpp convention
        cc = _np.cos(ang_rad)
        ss = _np.sin(ang_rad)
        xxr = cc * xx + ss * yy
        yyr = -ss * xx + cc * yy
        yy, xx = yyr, xxr
    matchx = _np.abs(xx) <= sx/2
    matchy = _np.abs(yy) <= sy/2
    Fout.field[matchx & matchy] = 0.0
    return Fout

def Strehl(Fin):
    """
    S = Strehl( Fin)

    :ref:`Calculates the Strehl value of the field <Strehl>`

    Args::

        Fin: input field

    Returns::

        S: Strehl value (real number).

    """
    normsq = _np.abs(Fin.field).sum()**2
    if normsq == 0.0:
        raise ValueError('Error in Strehl: Zero beam power')
    strehl = _np.real(Fin.field).sum()**2 + _np.imag(Fin.field).sum()**2
    strehl = strehl/normsq
    return strehl

def SubIntensity(Intens, Fin):
    """
    Fout = SubIntensity(Intens, Fin)

    :ref:`Substitutes  a given intensity distribution in the field with. <SubIntensity>`

    Args::

        Intens: N x N square array of real numbers >= 0
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Fout = Field.copy(Fin)
    Intens = _np.asarray(Intens)
    if Intens.shape != Fout.field.shape:
        raise ValueError('Intensity map has wrong shape')
    phi = _np.angle(Fout.field)
    Efield = _np.sqrt(Intens)
    Fout.field = Efield * _np.exp(1j * phi)
    return Fout

def SubPhase(Phi, Fin):
    """
    Fout = SubPhase(Phi, Fin)

    :ref:`Substitutes  a given phase distribution in the field with. <SubPhase>`

    Args::

        Phase: N x N square array of real numbers
        Fin: input field

    Returns::

        Fout: output field (N x N square array of complex numbers).

    """
    Fout = Field.copy(Fin)
    Phi = _np.asarray(Phi)
    if Phi.shape != Fin.field.shape:
        raise ValueError('Phase map has wrong shape')
    oldabs = _np.abs(Fout.field)
    Fout.field = oldabs * _np.exp(1j * Phi)
    return Fout

max_stars_repo_name: ldoyle/lightpipes

max_issues_repo_name: ldoyle/lightpipes

max_forks_repo_name: ldoyle/lightpipes

reason: import numpy,from scipy


text:
# -----------------------------------------------------------------------------
# Copyright (c) 2009-2016 Nicolas P. Rougier. All rights reserved.
# Distributed under the (new) BSD License.
# -----------------------------------------------------------------------------
import numpy as np
from glumpy import app, gl, gloo

vertex_1 = """
attribute vec2 position;
void main()
{
    gl_Position = vec4(position,0.0,1.0);
}
"""
fragment_1 = """
void main()
{
    gl_FragColor = vec4(10,10,10,1);
}
"""

vertex_2 = """
attribute vec2 position;
varying vec2 v_texcoord;
void main()
{
    gl_Position = vec4(0.85*position,0.0,1.0);
    v_texcoord = (position+1.0)/2.0;
}
"""
fragment_2 = """
uniform sampler2D texture;
varying vec2 v_texcoord;
void main()
{
    vec4 color = texture2D(texture, v_texcoord);
    gl_FragColor = color / vec4(20,20,20,1);
}
"""

window = app.Window(width=512, height=512)

@window.event
def on_draw(dt):
    window.clear()
    framebuffer.activate()
    quad_1.draw(gl.GL_TRIANGLE_STRIP)
    # The texture can be read as np.ndarray.
    out_texture = framebuffer.color[0].get()
    framebuffer.deactivate()
    quad_2.draw(gl.GL_TRIANGLE_STRIP)

texture = np.zeros((window.height,window.width,4),np.float32).view(gloo.TextureFloat2D)
framebuffer = gloo.FrameBuffer(color=[texture])
quad_1 = gloo.Program(vertex_1, fragment_1, count=4)
quad_1["position"] = (-1,-1), (-1,+1), (+1,-1), (+1,+1)
quad_2 = gloo.Program(vertex_2, fragment_2, count=4)
quad_2["position"] = (-1,-1), (-1,+1), (+1,-1), (+1,+1)
quad_2["texture"] = texture

app.run()


text:
import networkx as nx

def convert_graph_to_db_format(input_graph: nx.Graph, with_weights=False, cast_to_directed=False):
    """Converts a given graph into a DB format, which consists of two or three lists

        1. **Index list:** a list where the i-th position contains the index of the beginning of the list of adjacent nodes (in the second list).
        2. **Node list:** for each node, we list (in order) all the nodes which are adjacent to it.
        3. **Weight list:** if the weight parameter is True, includes the weights of the edges, corresponds to the nodes list

    **Assumptions:**

    The code has several preexisting assumptions:

        a) The nodes are labeled with numbers
        b) Those numbers are the sequence [0,...,num_of_nodes-1]
        c) If there are weights, they are floats
        d) If there are weights, they are initialized for all edges
        e) If there are weights, the weight key is 'weight'

    .. Note::
        The code behaves differently for directed and undirected graphs.
	For undirected graph, every edge is actually counted twice (p->q and q->p).

    Example::
		For the simple directed graph (0->1, 0->2,0->3,2->0,3->1,3->2):

		`Indices: [0, 3, 3, 4, 6]`
		`Neighbors: [1, 2, 3, 0, 1, 2]`

        Note that index[1] is the same as index[2]. That is because 1 has no neighbors, and so his neighbor list is of size 0, but we still need to have an index for the node on.

    For the same graph when it is undirected:
        `Indices: [0, 3, 5, 7, 10]`
        `Neighbors: [1, 2, 3, 0, 3, 0, 3, 0, 1, 2]`

        Note that the number of edges isn't doubled because in the directed version there is a bidirectional edge.

    :param graph: the nx.Graph object to convert
    :param with_weights: whether to create a weight list. Defaults to False.
    :param cast_to_directed: whether to cast the graph into a directed format

    :return: two or three lists: index,nodes, [weights]
    """

    if cast_to_directed:
        graph = input_graph.to_directed()
    else:
        graph = input_graph.copy()

    if graph.is_directed():
        # Color printing taken from https://www.geeksforgeeks.org/print-colors-python-terminal/
        print("\033[93m {}\033[00m".format('Note that the graph is processed as a directed graph'))

    indices = [0]  # The first neighbor list always starts at index 0
    neighbor_nodes = []

    nodes = [node for node in graph.nodes()]
    # print(nodes)
    nodes.sort()
    neighbors = [sorted([x for x in graph.neighbors(node)]) for node in nodes]

    # Create the indices and neighbor nodes lists
    for neighbor_list in neighbors:
        neighbor_list.sort()
        # print(neighbor_list)
        neighbor_nodes.extend(neighbor_list)
        indices.append(indices[-1] + len(neighbor_list))

    if with_weights:
        try:
            weights = [0] * len(neighbor_nodes)
            current_index = 0
            for node in nodes:
                for x in neighbors[node]:
                    w = graph[node][x]['weight']
                    weights[current_index] = w
                    current_index += 1

            return indices, neighbor_nodes, weights
        except KeyError:
            # Print in red
            print("\033[91m {}\033[00m".format('No weights defined, returning an empty list of weights'))
            print()
            return indices, neighbor_nodes, []

    return indices, neighbor_nodes

def convert_graph_to_db_dict(graph: nx.Graph, with_weights=False, cast_to_directed=False):
    """
    Encapsulates the convert_graph_to_db_format function by wrapping the returned lists with a dictionary.
    The dictionary has the keys 'indices','neighbors' , ['weights'].

    :param graph: the nx graph to convert
    :param with_weights: whether to return a weights list
    :param cast_to_directed: whether to process the graph as a directed one.
    :return: a dictionary with the specified keys and the lists as values.
    """
    ret_dict = {}
    if with_weights:
        i, n, w = convert_graph_to_db_format(graph, with_weights, cast_to_directed)
        ret_dict['indices'] = i
        ret_dict['neighbors'] = n
        ret_dict['weights'] = w
        ret_dict['with_weights'] = True
    else:
        i, n = convert_graph_to_db_format(graph, with_weights, cast_to_directed)
        ret_dict['indices'] = i
        ret_dict['neighbors'] = n
        ret_dict['with_weights'] = False

    ret_dict['directed'] = graph.is_directed()
    return ret_dict

max_stars_repo_name: Unknown-Data/QGCN

max_issues_repo_name: Unknown-Data/QGCN

max_forks_repo_name: Unknown-Data/QGCN

reason: import networkx


text:
# -*- coding: utf-8 -*-
##
##  generate_return_periods.py
##  spt_lsm_autorapid_process
##
##  Created by Alan D. Snow and Scott D. Christensen.
##  Copyright © 2015-2016 Alan D Snow and Scott D. Christensen. All rights reserved.
##  License: BSD-3 Clause

from datetime import datetime
import netCDF4 as nc
import numpy as np
from RAPIDpy.dataset import RAPIDDataset

def generate_return_periods(qout_file, return_period_file, storm_duration_days=7):
    """
    Generate return period from RAPID Qout file
    """

    #get ERA Interim Data Analyzed
    with RAPIDDataset(qout_file) as qout_nc_file:
        print "Setting up Return Periods File ..."
        return_period_nc = nc.Dataset(return_period_file, 'w')

        return_period_nc.createDimension('rivid', qout_nc_file.size_river_id)

        timeSeries_var = return_period_nc.createVariable('rivid', 'i4', ('rivid',))
        timeSeries_var.long_name = (
            'Unique NHDPlus COMID identifier for each river reach feature')

        max_flow_var = return_period_nc.createVariable('max_flow', 'f8', ('rivid',))
        return_period_20_var = return_period_nc.createVariable('return_period_20', 'f8', ('rivid',))
        return_period_10_var = return_period_nc.createVariable('return_period_10', 'f8', ('rivid',))
        return_period_2_var = return_period_nc.createVariable('return_period_2', 'f8', ('rivid',))

        lat_var = return_period_nc.createVariable('lat', 'f8', ('rivid',),
                                                  fill_value=-9999.0)
        lat_var.long_name = 'latitude'
        lat_var.standard_name = 'latitude'
        lat_var.units = 'degrees_north'
        lat_var.axis = 'Y'

        lon_var = return_period_nc.createVariable('lon', 'f8', ('rivid',),
                                                  fill_value=-9999.0)
        lon_var.long_name = 'longitude'
        lon_var.standard_name = 'longitude'
        lon_var.units = 'degrees_east'
        lon_var.axis = 'X'

        return_period_nc.variables['lat'][:] = qout_nc_file.qout_nc.variables['lat'][:]
        return_period_nc.variables['lon'][:] = qout_nc_file.qout_nc.variables['lon'][:]

        river_id_list = qout_nc_file.get_river_id_array()
        return_period_nc.variables['rivid'][:] = river_id_list

        print "Extracting Data and Generating Return Periods ..."
        time_array = qout_nc_file.get_time_array()
        num_years = int((datetime.utcfromtimestamp(time_array[-1])-datetime.utcfromtimestamp(time_array[0])).days/365.2425)
        time_steps_per_day = (24*3600)/float((datetime.utcfromtimestamp(time_array[1])-datetime.utcfromtimestamp(time_array[0])).seconds)
        step = max(1,int(time_steps_per_day * storm_duration_days))

        for comid_index, comid in enumerate(river_id_list):

            filtered_flow_data = qout_nc_file.get_daily_qout_index(comid_index,
                                                                   steps_per_group=step,
				                                   mode="max")
            sorted_flow_data = np.sort(filtered_flow_data)[:num_years:-1]

            rp_index_20 = int((num_years + 1)/20.0)
            rp_index_10 = int((num_years + 1)/10.0)
            rp_index_2 = int((num_years + 1)/2.0)

            max_flow_var[comid_index] = sorted_flow_data[0]
            return_period_20_var[comid_index] = sorted_flow_data[rp_index_20]
            return_period_10_var[comid_index] = sorted_flow_data[rp_index_10]
            return_period_2_var[comid_index] = sorted_flow_data[rp_index_2]

        return_period_nc.close()


text:
#!/usr/bin/env python

import numpy as np
import tensorflow as tf
from tensorflow.examples.tutorials.mnist import input_data

mnist = input_data.read_data_sets("data/", one_hot=True)
trainimg = mnist.train.images
trainlabel = mnist.train.labels
testimg = mnist.test.images
testlabel = mnist.test.labels

learning_rate = 0.1
training_epochs = 1000
batch_size = 256
display_step = 10

n_input = 784
n_output = 10

x = tf.placeholder("float", [None, n_input])
y = tf.placeholder("float", [None, n_output])

# Convolute layer
w1 = tf.Variable(tf.random_normal([3, 3, 1, 64]))
b1 = tf.Variable(tf.random_normal([64]))
input_r = tf.reshape(x, shape=[-1, 28, 28, 1])
conv1 = tf.nn.conv2d(input_r, w1, strides=[1, 1, 1, 1], padding="SAME")
conv2 = tf.nn.bias_add(conv1, b1)
conv3 = tf.nn.relu(conv2)
pool = tf.nn.max_pool(conv3,
                      ksize=[1, 2, 2, 1],
                      strides=[1, 2, 2, 1],
                      padding="SAME")

# Full connected layer
w2 = tf.Variable(tf.random_normal([14 * 14 * 64, n_output]))
b2 = tf.Variable(tf.random_normal([n_output]))
dense = tf.reshape(pool, [-1, w2.get_shape().as_list()[0]])
pred = tf.add(tf.matmul(dense, w2), b2)
cost = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(pred, y))
optm = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(cost)

with tf.Session() as sess:
    sess.run(tf.initialize_all_variables())

    for epoch in range(training_epochs):
        avg_cost = 0
        total_batch = int(mnist.train.num_examples / batch_size)

        for i in range(total_batch):
            batch_xs, batch_ys = mnist.train.next_batch(batch_size)
            sess.run(optm, feed_dict={x: batch_xs, y: batch_ys})
            avg_cost += sess.run(cost, feed_dict={x: batch_xs,
                                              y: batch_ys}) / total_batch

        if epoch % display_step == 0:
            print("Epoch: {}, cost: {}".format(epoch, avg_cost))


text:
"""
Author: Nikhil Devraj
Adapted from Kitchen2D: https://github.com/zi-w/Kitchen2D
"""

from Box2D import *
from Box2D.b2 import *
from env_constants import *

import pygame
import numpy as np

class guiWorld:
    def __init__(self, caption='PyBox2D Simulator', overclock=None):
        """
        Args:
            caption:  caption on the window
            overclock: number of frames to skip when showing graphical interface
        """
        self.screen = pygame.display.set_mode((SCREEN_WIDTH_PX, SCREEN_HEIGHT_PX))
        pygame.display.set_caption(caption)
        self.clock = pygame.time.Clock()
        self.overclock = overclock
        self.screen_origin = b2Vec2(SCREEN_WIDTH/2., TABLE_HEIGHT)
        self.colors = {
            'countertop': (50, 50, 50, 255),
            'gripper': (244, 170, 66, 255),
            'cup': (15, 0, 100, 0),
            'stirrer': (163, 209, 224, 255),
            'spoon': (73, 11, 61, 255),
            'obstacle': (123, 128, 120, 255),
            'sensor': (255, 255, 255, 20),
            #'sensor': (155, 155, 155, 255), # To visualize the sensor
            'default': (81, 81, 81, 255),
            'faucet': (175, 175, 175, 255),
            'water': (26, 130, 252),
            'block': (0, 99, 0),
            'coffee': (165, 42, 42),
            'cream': (225, 225, 225),
            'sugar': (255, 192, 203),
        }


text:
import os
import unittest
import numpy as np
import pandas as pd
from context import performance as pi

# Change working directory
# This enable running tests from repository root
if os.getcwd() != os.path.abspath(os.path.dirname(__file__)):
    os.chdir('tests/')

# Test results
class ResultsCalmar(unittest.TestCase):
    # Input data
    test_yearly = pd.Series([100, 120, 110, 100, 110, 130])
    test_monthly = pd.read_csv('test_data/ohlc_monthly.csv').drop(columns=['Adj Close'])
    test_daily = pd.read_csv('test_data/series_daily.csv')

    # Expected results
    # Made with excel
    results_yearly = 0.05387395206 / 0.1666666667
    results_monthly = 0.4075254214 / 0.3070333382
    results_daily = -0.7543788728 / 0.06131828655
    results_daily_all_days = -0.8691260753 / 0.06131828655

    def test_result_yearly(self):
        '''result for yearly data expected'''
        result = pi.calmar(self.test_yearly, frequency='Y')
        self.assertAlmostEqual(self.results_yearly, result, places=4)

    def test_result_monthly_fraction(self):
        '''result for monthly data expected'''
        result = pi.calmar(self.test_monthly, frequency='M')
        self.assertAlmostEqual(self.results_monthly, result, places=4)

    def test_result_daily_fraction(self):
        '''result for daily data expected'''
        result = pi.calmar(self.test_daily, frequency='D')
        self.assertAlmostEqual(self.results_daily, result, places=4)

    def test_result_daily_fraction_all_days(self):
        '''result for daily with all days option'''
        result = pi.calmar(self.test_daily, frequency='D', only_business=False)
        self.assertAlmostEqual(self.results_daily_all_days, result, places=4)

# Test inputs
class BadInputSortino(unittest.TestCase):
    # Input data
    test_yearly = pd.Series(range(100, 201, 20))
    test_yearly_list = range(100, 201, 20)
    test_no_close = pd.read_csv('test_data/ohlc_monthly.csv').drop(columns=['Close', 'Adj Close'])
    test_more_close = pd.read_csv('test_data/ohlc_monthly.csv')

    def test_input_data_not_df_or_series(self):
        '''data is not a valid type raise TypeError'''
        self.assertRaises(TypeError, pi.calmar, self.test_yearly_list)

    def test_input_df_no_close(self):
        '''input df has no close column'''
        self.assertRaises(IndexError, pi.calmar, self.test_no_close)

    def test_input_df_more_close(self):
        '''input df has more than 1 close column'''
        self.assertRaises(KeyError, pi.calmar, self.test_more_close)

    def test_frequency_parameter(self):
        '''invalid string for frequency raise ValueError'''
        self.assertRaises(ValueError, pi.calmar, self.test_yearly, frequency='mo')

    def test_business_parameter(self):
        '''invalid value for only_business raise TypeError'''
        self.assertRaises(TypeError, pi.calmar, self.test_yearly, only_business='Yes')

if __name__ == '__main__':
    unittest.main()

max_stars_repo_name: gcontarini/finance-toolkit

max_issues_repo_name: gcontarini/finance-toolkit

max_forks_repo_name: gcontarini/finance-toolkit


text:
import numpy as np
import scipy.signal

class QuadAltitude(object):
    """Discrete-time model of altitude of a quadcopter

    Parameters
    ----------
    dt : float, optional
        Discretization time (zero-order hold), by default 0.001
    m : float, optional
        Mass of the quadcopter, by default 0.028
    g : float, optional
        Gravitational acceleration, by default 9.807
    """

    def __init__(self, dt: float = 0.001,
                 m: float = 0.028, g: float = 9.807):

        assert isinstance(g, float), 'Gravity needs to be of type float'
        assert isinstance(m, float), 'Mass of quadcopter \
            needs to be of type float'

        self.g = g
        self.m = m
        self.dt = dt

        self.vel = 0
        self.pos = 0
        self.acc = 0

    def __call__(self, u: float) -> float:
        """Input port to the quadcopter

        Parameters
        ----------
        u : float
            Thrust of the rotors

        Returns
        -------
        float:
            Current altitude of quadcopter

        """

        acc = u * (1/self.m) - self.g

        self.vel += acc * self.dt
        self.pos += self.vel * self.dt
        self.acc = acc

        # Include ground
        if self.pos <= 0:
            self.vel = 0
            self.pos = 0

        return self.pos

    @property
    def states(self):
        """Return the state of the quadcopter
        """

        return self.pos, self.vel, self.acc

    @states.setter
    def states(self, value):
        """Sets the state of the quadcopter to value
        value = (pos, vel, acc). If a state is None, it is not changed.

        Parameters
        ----------
        value : tuple
            new state vector of quadcopter
        """

        self.pos = value[0] if isinstance(value[0], (float, int)) else self.pos
        self.vel = value[1] if isinstance(value[1], (float, int)) else self.vel
        self.acc = value[2] if isinstance(value[2], (float, int)) else self.acc

    def reset(self):
        """Reset position and velocity to zero
        """

        self.pos = 0.
        self.vel = 0.

class System(object):

    """Python class to generate a first-or second-order process
    for simulation and verification purposes

    Parameters
    ----------
    K : float, optional
        system gain, by default 1.0
    T : float, optional
        time constant, by default 1.0
    D : float, optional
        damping factor, by default None
    delta : float, optional
        discretization time in seconds, by default 0.001
    """

    def __init__(self, K: float = 1.0, T: float = 1.0,
                 D: float = None, delta: float = 0.001) -> None:

        assert isinstance(K, float)
        assert isinstance(T, float)
        assert delta > 0, "Sampling time has to be positive"

        # First-order process
        num = np.array([K]).reshape(-1)
        den = np.array([T, 1]).reshape(-1)

        if D is not None:
            assert isinstance(D, float),\
                "Damping factor D must be type float"
            # Second-order process
            den = np.array([T**2, 2 * D * T, 1]).reshape(-1)

        system = scipy.signal.tf2ss(num, den)

        system = scipy.signal.cont2discrete(system, delta)

        self.A = system[0]
        self.B = system[1]
        self.C = system[2]

        self.x = np.zeros((len(den) - 1, 1), dtype=np.float64)

    def __call__(self, u: float) -> float:
        """System response w.r.t. control signal u

        Parameters
        ----------
        u : float
            control signal u

        Returns
        -------
        float
            system response y
        """

        self.x = self.A.dot(self.x) + self.B.dot(u)

        return float(self.C.dot(self.x))

class RandomSystem(object):

    def __init__(self, states=1) -> None:
        """Create random discrete state space system

        Args:
            states (int, optional): number of states. Defaults to 1.

        Returns:
            None: [description]
        """

        import control

        system = control.drss(states=states)

        self.A = np.array(system.A)
        self.B = np.array(system.B)
        self.C = np.array(system.C)
        self.D = np.array(system.D)

        self.x = np.zeros((states, 1), dtype=np.float64)

    def __call__(self, u: float) -> float:

        self.x = self.A.dot(self.x) + self.B.dot(u)
        return float(self.C.dot(self.x) + self.D.dot(u))

max_stars_repo_name: onguntoglu/pyadrc

max_issues_repo_name: onguntoglu/pyadrc

max_forks_repo_name: onguntoglu/pyadrc

reason: import numpy,import scipy


text:
import sys
import math
import numpy as np

"""This file gathers constant values used for different robots such as the lengths of the leg parts, the angle conventions, initial positions, etc.
Also, to maintain compatibility with legacy code, the user can chose whether the input and/or outputs of the kinematic functions are in (rads or degs) and (mm or meters).
Values in this file should however always be written in meters and in rads.
"""

BIOLOID = "BIOLOID"
PHANTOMX = "PHANTOMX"
PHANTOMX_SIMULATION = "PHANTOMX_SIMULATION"
ARM_SIMULATION = "ARM_SIMULATION"
AX12 = "AX12"
MOTOR_TYPE = AX12
ROBOT_TYPE = PHANTOMX_SIMULATION
USING_SIMPLE_ROBOT = True
USE_RADS_INPUT = False
USE_RADS_OUTPUT = False
USE_MM_INPUT = False
USE_MM_OUTPUT = False

if ROBOT_TYPE == PHANTOMX:
    constL1 = 0.001 * 54.8
    constL2 = 0.001 * 65.3
    constL3 = 0.001 * 133
    theta2Correction = 16.0 * math.pi / 180.0
    theta2ExtraCorrection = 0
    theta3Correction = 43.76 * math.pi / 180.0  # -theta2Correction ??
    THETA3_MOTOR_SIGN = 1
    THETA2_MOTOR_SIGN = -1
    THETA1_MOTOR_SIGN = 1
    USE_RADS_INPUT = False
    USE_RADS_OUTPUT = False
    USE_MM_INPUT = False
    USE_MM_OUTPUT = False
    Z_DIRECTION = -1
    LEG_ANGLES = [
        -math.pi / 4,
        math.pi / 4,
        math.pi / 2,
        3 * math.pi / 4,
        -3 * math.pi / 4,
        -math.pi / 2,
    ]
elif ROBOT_TYPE == BIOLOID:
    constL1 = 0.001 * 51
    constL2 = 0.001 * 63.7
    constL3 = 0.001 * 93
    # Angle to match the theory with reality for theta 2 (measures of the triangle are 22.5, 60.7, 63.7). => Angle =  -20.69
    theta2Correction = -20.69 * math.pi / 180.0
    # Same goes for theta 3 : +90 - 20.69 - a. Where a = asin(8.2/93) = 5.06
    theta3Correction = (90 + theta2Correction - 5.06) * math.pi / 180.0
    THETA3_MOTOR_SIGN = -1
    THETA2_MOTOR_SIGN = 1
    THETA1_MOTOR_SIGN = 1
    USE_RADS_INPUT = False
    USE_RADS_OUTPUT = False
    USE_MM_INPUT = True
    USE_MM_OUTPUT = True
    Z_DIRECTION = -1
    LEG_ANGLES = [math.pi / 2, 0, 0, -math.pi / 2, math.pi, math.pi]
elif ROBOT_TYPE == PHANTOMX_SIMULATION:
    constL1 = 0.054
    constL2 = 0.0645
    constL3 = 0.155  # Approximate value (different from real robot and couldn't find value in URDF yet)
    theta2Correction = -16.0 * math.pi / 180.0  # Measured on real robot
    theta3Correction = (
        -43.76 * math.pi / 180.0 + theta2Correction
    )  # Measured on real robot
    THETA3_MOTOR_SIGN = -1
    THETA2_MOTOR_SIGN = 1
    THETA1_MOTOR_SIGN = 1
    USE_RADS_INPUT = True
    USE_RADS_OUTPUT = True
    USE_MM_INPUT = False
    USE_MM_OUTPUT = False
    Z_DIRECTION = 1
    LEG_ANGLES = [
        math.pi / 4,
        -math.pi / 4,
        -math.pi / 2,
        -3 * math.pi / 4,
        3 * math.pi / 4,
        math.pi / 2,
    ]

    LEG_ANGLES_2 = np.negative(LEG_ANGLES)

    LEG_CENTER_POS = [
        (0.1248, -0.06164, 0.001116),
        (0.1248, 0.06164, 0.001116),
        (0, 0.1034, 0.001116),
        (-0.1248, 0.06164, 0.001116),
        (-0.1248, -0.06164, 0.001116),
        (0, -0.1034, 0.001116),
    ]

elif ROBOT_TYPE == ARM_SIMULATION:
    constL1 = 0.085
    constL2 = 0.185
    constL3 = 0.250
    theta2Correction = 0
    theta3Correction = 0
    THETA3_MOTOR_SIGN = 1
    THETA2_MOTOR_SIGN = 1
    THETA1_MOTOR_SIGN = 1
    USE_RADS_INPUT = True
    USE_RADS_OUTPUT = True
    USE_MM_INPUT = False
    USE_MM_OUTPUT = False
    Z_DIRECTION = 1
else:
    print("ERROR: Unknwon ROBOT_TYPE '{}'".format(ROBOT_TYPE))
    sys.exit()


text:
import sys
import pandas as pd
import numpy as np
import itertools
from sklearn.preprocessing import RobustScaler
from xgboost import XGBClassifier # Assumes XGBoost v0.6
from evaluate_model import evaluate_model

dataset = sys.argv[1]

pipeline_components = [RobustScaler, XGBClassifier]
pipeline_parameters = {}

n_estimators_values = [10, 50, 100, 500]
learning_rate_values = [0.01, 0.1, 0.5, 1.0, 10.0, 50.0, 100.0]
gamma_values = np.arange(0., 0.51, 0.05)
max_depth_values = [1, 2, 3, 4, 5, 10, 20, 50, None]
subsample_values = np.arange(0.0, 1.01, 0.1)
random_state = [324089]

all_param_combinations = itertools.product(n_estimators_values, learning_rate_values, gamma_values, max_depth_values, subsample_values, random_state)
pipeline_parameters[XGBClassifier] = \
   [{'n_estimators': n_estimators, 'learning_rate': learning_rate, 'gamma': gamma, 'max_depth': max_depth, 'subsample': subsample, 'seed': random_state, 'nthread': 1}
     for (n_estimators, learning_rate, gamma, max_depth, subsample, random_state) in all_param_combinations]

evaluate_model(dataset, pipeline_components, pipeline_parameters)

max_stars_repo_name: lacava/sklearn-benchmarks

max_issues_repo_name: LeMGarouani/AMLBID

max_forks_repo_name: arunsinghyadav/sklearn-benchmarks


text:
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 1, 1000)

def P(R):
    res = []
    for r in R:
        if 0 <= r < 1:
            val = 2*r*np.log(1+r) + (1 - r)*np.log(1 - r*r) - np.log(1 + r*r)
            res.append(val)
        elif r == 1:
            val = np.log(2)
            res.append(val)
    return np.array(res)

y = P(x)
y_2 = 1 - y
r_equal = 0.9464586434189278
r_mean = 0.67369968572800811
fig, ax = plt.subplots()
ax.plot(x, y, label=r"$P(\tau_{ro}'' < \tau_{lc}'')$", color='#1f77b4')
ax.plot(x, y_2, label=r"$P(\tau_{ro}'' > \tau_{lc}'')$", color='#d62728')
ax.plot(r_equal, 1/2, marker='x', color='black')
ax.plot(r_mean, P([r_mean]), marker='x', color='black')
ax.plot(r_mean, 1 - P([r_mean]), marker='x', color='black')
plt.xlabel(r"$r$")
plt.ylabel(r"$P$")
plt.xlim(0, 1)
plt.ylim(0, 1)
x_ticks = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, r_mean, 0.7, 0.8, 0.9, r_equal, 1.]
x_ticks2 = ['0', '0.1', '0.2', '0.3', '0.4', '0.5', '0.6', r'$r_{mean}$',
            '0.7', '0.8', '0.9', r'$r_{equal}$', '1']
plt.xticks(x_ticks, x_ticks2, rotation=-40)
plt.yticks(np.arange(0, 1.1, step=0.1))
plt.legend()
plt.savefig("RatioProbability.png", dpi=1000,
            bbox_inches='tight', pad_inches=0.01)
plt.savefig("RatioProbability.pdf",
            bbox_inches='tight', pad_inches=0.01)
plt.show()

max_stars_repo_name: Horep/Barn-And-Pole-Paradox

max_issues_repo_name: Horep/Barn-And-Pole-Paradox

max_forks_repo_name: Horep/Barn-And-Pole-Paradox