text:
##################################################
# Copyright (c) Xuanyi Dong [GitHub D-X-Y], 2019 #
######################################################################################
# One-Shot Neural Architecture Search via Self-Evaluated Template Network, ICCV 2019 #
######################################################################################
import os, sys, time, glob, random, argparse
import numpy as np
from copy import deepcopy
import torch
import torch.nn as nn
from pathlib import Path
lib_dir = (Path(__file__).parent / '..' / 'lib').resolve()
if str(lib_dir) not in sys.path: sys.path.insert(0, str(lib_dir))
from config_utils import load_config, dict2config, configure2str
from datasets     import get_datasets, get_nas_search_loaders
from procedures   import prepare_seed, prepare_logger, save_checkpoint, copy_checkpoint, get_optim_scheduler
from procedures.transfer import get_search_methods
from utils        import get_model_infos, obtain_accuracy
from log_utils    import AverageMeter, time_string, convert_secs2time, write_results
from models       import get_cell_based_tiny_net, get_search_spaces, load_net_from_checkpoint, FeatureMatching, CellStructure as Structure
from nas_201_api  import NASBench201API as API

def get_lr(optimizer):
    for param_group in optimizer.param_groups:
        return param_group['lr']

def get_best_arch(xloader, network, n_samples):
  # setn evaluation
  with torch.no_grad():
    network.eval()
    archs, valid_accs = network.module.return_topK(n_samples), []
    #print ('obtain the top-{:} architectures'.format(n_samples))
    loader_iter = iter(xloader)
    for i, sampled_arch in enumerate(archs):
      network.module.set_cal_mode('dynamic', sampled_arch)
      try:
        inputs, targets = next(loader_iter)
      except:
        loader_iter = iter(xloader)
        inputs, targets = next(loader_iter)

      _, logits = network(inputs)
      val_top1, val_top5 = obtain_accuracy(logits.cpu().data, targets.data, topk=(1, 5))

      valid_accs.append( val_top1.item() )
      #print ('--- {:}/{:} : {:} : {:}'.format(i, len(archs), sampled_arch, val_top1))

    best_idx = np.argmax(valid_accs)
    best_arch, best_valid_acc = archs[best_idx], valid_accs[best_idx]
    return best_arch, best_valid_acc

def main(args):
  assert torch.cuda.is_available(), 'CUDA is not available.'
  torch.backends.cudnn.enabled   = True
  torch.backends.cudnn.benchmark = False
  torch.backends.cudnn.deterministic = True
  torch.set_num_threads( args.workers )
  prepare_seed(args.rand_seed)
  logger = prepare_logger(args)
  total_time = time.time()

  train_data, valid_data, xshape, class_num = get_datasets(args.dataset, args.data_path, args.cutout_length)
  config = load_config(args.config_path, {'class_num': class_num, 'xshape': xshape}, logger)
  search_loader, _, valid_loader = get_nas_search_loaders(train_data, valid_data, args.dataset, 'configs/nas-benchmark/', \
                                        config.batch_size if not hasattr(config, "test_batch_size") else (config.batch_size, config.test_batch_size), args.workers)
  logger.log('||||||| {:10s} ||||||| Search-Loader-Num={:}, Valid-Loader-Num={:}, batch size={:}'.format(args.dataset, len(search_loader), len(valid_loader), config.batch_size))
  logger.log('||||||| {:10s} ||||||| Config={:}'.format(args.dataset, config))

  search_space = get_search_spaces('cell', args.search_space_name)
  if args.fixed_genotype:
    fixed_genotype = Structure.str2structure( args.fixed_genotype )
    model_config = dict2config({'name': args.nas_name, 'C': args.channel, 'N': args.num_cells,
                                'max_nodes': args.max_nodes, 'num_classes': class_num,
                                'space'    : search_space,
                                'affine'   : False, 'track_running_stats': bool(args.track_running_stats),
                                "fixed_genotype"    :   fixed_genotype, "search_position" : args.pos}, None)
  else:
    model_config = dict2config({'name': args.nas_name, 'C': args.channel, 'N': args.num_cells,
                                'max_nodes': args.max_nodes, 'num_classes': class_num,
                                'space'    : search_space,
                                'affine'   : False, 'track_running_stats': bool(args.track_running_stats)}, None)
  logger.log('search space : {:}'.format(search_space))
  logger.log('model-config : {:}'.format(model_config))

  search_model = get_cell_based_tiny_net(model_config)
  w_optimizer, w_scheduler, criterion = get_optim_scheduler(search_model.get_weights(), config)
  a_optimizer = torch.optim.Adam(search_model.get_alphas(), lr=args.arch_learning_rate, betas=(0.5, 0.999), weight_decay=args.arch_weight_decay)
  logger.log('w-optimizer : {:}'.format(w_optimizer))
  logger.log('a-optimizer : {:}'.format(a_optimizer))
  logger.log('w-scheduler : {:}'.format(w_scheduler))
  logger.log('criterion   : {:}'.format(criterion))
  flop, param  = get_model_infos(search_model, xshape)
  logger.log('{:}'.format(search_model))
  logger.log('FLOP = {:.2f} M, Params = {:.2f} MB'.format(flop, param))
  logger.log('search-space : {:}'.format(search_space))
  if args.search_space_name != "nas-bench-201" or args.num_cells != 5:
    api = None
  else:
    api = API(args.arch_nas_dataset)
    logger.log('{:} create API = {:} done'.format(time_string(), api))
  last_info, model_base_path, model_best_path = logger.path('info'), logger.path('model'), logger.path('best')

  network, criterion = torch.nn.DataParallel(search_model).cuda(), criterion.cuda()

  if last_info.exists() and not args.overwrite: # automatically resume from previous checkpoint
    logger.log("=> loading checkpoint of the last-info '{:}' start".format(last_info))
    last_info   = torch.load(last_info)
    start_epoch = last_info['epoch']
    checkpoint  = torch.load(last_info['last_checkpoint'])
    genotypes   = checkpoint['genotypes']
    arch_params = checkpoint['arch_params']
    search_losses = checkpoint['search_losses']
    valid_losses = checkpoint['valid_losses']
    search_arch_losses = checkpoint['search_arch_losses']
    search_model.load_state_dict( checkpoint['search_model'] )
    w_optimizer.load_state_dict ( checkpoint['w_optimizer'] )
    w_scheduler.load_state_dict ( checkpoint['w_scheduler'] )
    a_optimizer.load_state_dict ( checkpoint['a_optimizer'] )
    matching_layers.load_state_dict( checkpoint['matching_layers'] )
    logger.log("=> loading checkpoint of the last-info '{:}' start with {:}-th epoch.".format(last_info, start_epoch))
  else:
    logger.log("=> do not find the last-info file : {:}".format(last_info))
    start_epoch, genotypes = 0, {-1: search_model.genotype()}
    search_losses, search_arch_losses = {}, {}
    valid_losses, valid_acc1s, valid_acc5s = {}, {'best': -1}, {}
    arch_params = {}
  # start training
  (search_w_func, search_a_func), valid_func = get_search_methods(args.nas_name, 20)
  # search w training
  start_time, search_w_time, epoch_time, total_epoch = time.time(), AverageMeter(), AverageMeter(), config.epochs + config.warmup
  for epoch in range(start_epoch, total_epoch):
      w_scheduler.update(epoch, 0.0)
      need_time = 'Time Left: {:}'.format( convert_secs2time(epoch_time.val * (total_epoch-epoch), True) )
      epoch_str = '{:03d}-{:03d}'.format(epoch, total_epoch)
      if args.nas_name == "GDAS":
          search_model.set_tau( args.tau_max - (args.tau_max-args.tau_min) * epoch / (total_epoch-1) )
      logger.log('\n[Search the {:}-th epoch] {:}, LR={:}'.format(epoch_str, need_time, min(w_scheduler.get_lr())))
      search_w_loss, search_w_top1, search_w_top5 \
            = search_w_func(search_loader, network, criterion, w_scheduler, w_optimizer, epoch_str, args.print_freq, logger)
      search_w_time.update(time.time() - start_time)
      logger.log('[{:}] search [base] : loss={:.2f}, accuracy@1={:.2f}%, accuracy@5={:.2f}%, time-cost={:.1f} s'.format(epoch_str, search_w_loss, search_w_top1, search_w_top5, search_w_time.sum))
      search_losses[epoch] = search_w_loss
      # measure elapsed time
      epoch_time.update(time.time() - start_time)
      start_time = time.time()
  # save checkpoint
  logger.log('<<<--->>> Supernet Train Complete.')
  supernet_save_path = logger.path('log') / 'seed-{:}-supernet.pth'.format(logger.seed)
  save_path = save_checkpoint(
             {'epoch' : epoch + 1,
              'args'  : deepcopy(args),
              'shared_cnn'  : search_model.state_dict(),
              'w_optimizer' : w_optimizer.state_dict(),
              'w_scheduler' : w_scheduler.state_dict(),
              'search_losses' : search_losses},
              supernet_save_path, logger)
  # search a training
  valid_time = AverageMeter()
  start_time, search_a_time, epoch_time, total_epoch = time.time(), AverageMeter(), AverageMeter(), 250 + config.warmup #config.epochs + config.warmup
  for epoch in range(start_epoch, total_epoch):
      need_time = 'Time Left: {:}'.format( convert_secs2time(epoch_time.val * (total_epoch-epoch), True) )
      epoch_str = '{:03d}-{:03d}'.format(epoch, total_epoch)
      if args.nas_name == "GDAS":
          search_model.set_tau( args.tau_max - (args.tau_max-args.tau_min) * epoch / (total_epoch-1) )
      logger.log('\n[Search the {:}-th epoch] {:}, LR={:}'.format(epoch_str, need_time, min(w_scheduler.get_lr())))
      search_a_loss, search_a_top1, search_a_top5 \
            = search_a_func(search_loader, network, criterion, a_optimizer, epoch_str, args.print_freq, logger)
      search_a_time.update(time.time() - start_time)
      logger.log('[{:}] search [arch] : loss={:.2f}, accuracy@1={:.2f}%, accuracy@5={:.2f}%'.format(epoch_str, search_a_loss, search_a_top1, search_a_top5))
      # validation
      valid_start_time = time.time()
      if args.nas_name == "SETN":
          genotype, _ = get_best_arch(valid_loader, network, args.select_num)
          network.module.set_cal_mode('dynamic', genotype)
      else:
          genotype = search_model.genotype()
      valid_a_loss , valid_a_top1 , valid_a_top5  = valid_func(valid_loader, network, criterion)
      valid_time.update(time.time() - valid_start_time)
      logger.log('[{:}] evaluate : loss={:.2f}, accuracy@1={:.2f}%, accuracy@5={:.2f}%, time-cost={:1f} s'.format(epoch_str, valid_a_loss, valid_a_top1, valid_a_top5, valid_time.sum))
      # check the best accuracy
      search_arch_losses[epoch] = search_a_loss
      valid_losses[epoch] = valid_a_loss
      valid_acc1s[epoch] = valid_a_top1
      valid_acc5s[epoch] = valid_a_top5
      genotypes[epoch] = genotype
      with torch.no_grad():
          arch_param = nn.functional.softmax(search_model.arch_parameters, dim=-1).cpu().numpy()
      arch_params[epoch] = arch_param
      logger.log('<<<--->>> The {:}-th epoch : {:}'.format(epoch_str, genotypes[epoch]))
      if valid_a_top1 > valid_acc1s['best']:
          valid_acc1s['best'] = valid_a_top1
          genotypes['best']   = genotypes[epoch]
          arch_params['best'] = arch_param
          find_best = True
      else: find_best = False
      # save checkpoint
      save_path = save_checkpoint({'epoch' : epoch + 1,
                  'args'  : deepcopy(args),
                  'search_model': search_model.state_dict(),
                  'w_optimizer' : w_optimizer.state_dict(),
                  'a_optimizer' : a_optimizer.state_dict(),
                  'w_scheduler' : w_scheduler.state_dict(),
                  'arch_params' : arch_params,
                  'genotypes'   : deepcopy(genotypes),
                  'fixed_genotype' : args.fixed_genotype,
                  "search_position" : args.pos,
                  "search_losses" : deepcopy(search_losses),
                  "search_arch_losses" : deepcopy(search_arch_losses),
                  "valid_losses" : deepcopy(valid_losses),
                  "valid_acc1s" : deepcopy(valid_acc1s),
                  "valid_acc5s" : deepcopy(valid_acc5s),
                  },
                  model_base_path, logger)
      last_info = save_checkpoint({
              'epoch': epoch + 1,
              'args' : deepcopy(args),
              'last_checkpoint': save_path,
              }, logger.path('info'), logger)
      if find_best:
          logger.log('<<<--->>> The {:}-th epoch : find the highest validation accuracy : {:.2f}%.'.format(epoch_str, valid_a_top1))
          copy_checkpoint(model_base_path, model_best_path, logger)
      logger.log('arch-parameters :\n{:}'.format( arch_param ))
      if api is not None: logger.log('{:}'.format(api.query_by_arch( genotype )))
      # measure elapsed time
      epoch_time.update(time.time() - start_time)
      start_time = time.time()

  logger.log('\n' + '-'*100)
  # check the performance from the architecture dataset
  logger.log('{:} : run {:} epochs, cost w-{:.1f} + a-{:.1f} s, last-geno is {:}.'.format(args.nas_name, total_epoch, search_w_time.sum, search_a_time.sum, genotypes[total_epoch-1]))
  if api is not None: logger.log('{:}'.format( api.query_by_arch(genotypes[total_epoch-1]) ))
  logger.log('The best-geno is {:} with Valid Acc {:}.'.format(genotypes['best'], valid_acc1s['best']))
  if api is not None: logger.log('{:}'.format( api.query_by_arch(genotypes['best']) ))
  logger.log("[Time cost] total: {:}, search w: {:}, search a: {:}, valid: {:}".format(convert_secs2time(time.time() - total_time, True), convert_secs2time(search_w_time.sum, True), convert_secs2time(search_a_time.sum, True), convert_secs2time(valid_time.sum, True)))
  logger.close()

if __name__ == '__main__':
  parser = argparse.ArgumentParser("Part Re-Search2")
  parser.add_argument('--exp_name',           type=str,  default="",     help='Experiment name')
  parser.add_argument('--overwrite',          type=bool, default=False,  help='Overwrite the existing results')
  parser.add_argument("--nas_name",           type=str,  default="GDAS", help="NAS algorithm to use")
  parser.add_argument("--fixed_genotype",     type=str,   help="Part cell search architecture")
  parser.add_argument("--pos",                type=int,   help="Part cell search stage: [0,1,2]")
  # data
  parser.add_argument('--data_path',          type=str,   default=os.environ['TORCH_HOME'] + "/cifar.python", help='Path to dataset')
  parser.add_argument('--dataset',            type=str,   default='cifar100', choices=['cifar10', 'cifar100', 'ImageNet16-120'], help='Choose between Cifar10/100 and ImageNet-16.')
  parser.add_argument('--cutout_length',      type=int,   default=-1,      help='The cutout length, negative means not use.')
  # channels and number-of-cells
  parser.add_argument('--search_space_name',  type=str,   default="nas-bench-201", help='The search space name.')
  parser.add_argument('--max_nodes',          type=int,   default=4, help='The maximum number of nodes.')
  parser.add_argument('--channel',            type=int,   default=16, help='The number of channels.')
  parser.add_argument('--num_cells',          type=int,   default=2, help='The number of cells in one stage.')
  parser.add_argument('--track_running_stats',type=int,   default=0, choices=[0,1],help='Whether use track_running_stats or not in the BN layer.')
  parser.add_argument('--config_path',        type=str,   default="configs/research/nasnet-search-E250.config", help='The path of the configuration.')
  parser.add_argument('--model_config',       type=str,   help='The path of the model configuration. When this arg is set, it will cover max_nodes / channels / num_cells.')
  # architecture leraning rate
  parser.add_argument('--arch_learning_rate', type=float, default=3e-4, help='learning rate for arch encoding')
  parser.add_argument('--arch_weight_decay',  type=float, default=1e-3, help='weight decay for arch encoding')
  # GDAS
  parser.add_argument('--tau_min',            type=float, default=0.1,  help='The minimum tau for Gumbel')
  parser.add_argument('--tau_max',            type=float, default=10,   help='The maximum tau for Gumbel')
  # SETN
  parser.add_argument("--select_num",         type=int,   default=100,  help="The number of architectures to be sampled for evaluation")
  # log
  parser.add_argument('--workers',            type=int,   default=8,    help='number of data loading workers')
  parser.add_argument('--save_dir',           type=str,   default="./output/nasnet-search2",     help='Folder to save checkpoints and log.')
  parser.add_argument('--arch_nas_dataset',   type=str,   default=os.environ['TORCH_HOME'] + "/NAS-Bench-201-v1_1-096897.pth", help='The path to load the architecture dataset (tiny-nas-benchmark).')
  parser.add_argument('--print_freq',         type=int,   default=100, help='print frequency (default: 100)')
  parser.add_argument('--rand_seed',          type=int,   default=-1, help='manual seed')
  args = parser.parse_args()
  if args.rand_seed is None or args.rand_seed < 0: args.rand_seed = random.randint(1, 100000)
  if args.exp_name != "":
      args.save_dir = "./output/{}-n{}/nasnet-search2/{}".format(args.dataset, args.num_cells, args.exp_name)
  main(args)

max_stars_repo_name: yongsu-bloo/AutoDL-Projects

max_issues_repo_name: yongsu-bloo/AutoDL-Projects

max_forks_repo_name: yongsu-bloo/AutoDL-Projects


text:
"""
Cardinality
- The values of a categorical variable are selected from a group of categories, also called labels.
- For example, in the variable gender the categories or labels are male and female,
    whereas in the variable city the labels can be London, Manchester, Brighton and so on.

Different categorical variables contain different number of labels or categories.
    The variable gender contains only 2 labels,
    but a variable like city or postcode, can contain a huge number of different labels.

The number of different labels within a categorical variable is known as cardinality.
    A high number of labels within a variable is known as high cardinality.

Are multiple labels in a categorical variable a problem?
High cardinality may pose the following problems:
- Variables with too many labels tend to dominate over those with only a few labels,
    particularly in Tree based algorithms.

- A big number of labels within a variable may introduce noise with little,
    if any, information, therefore making machine learning models prone to over-fit.

- Some of the labels may only be present in the training data set, but not in the test set,
    therefore machine learning algorithms may over-fit to the training set.

- Contrarily, some labels may appear only in the test set, therefore leaving the machine learning
    algorithms unable to perform a calculation over the new (unseen) observation.

- In particular, tree methods can be biased towards variables with lots of labels (variables with high cardinality).
    Thus, their performance may be affected by high cardinality.

Below, I will show the effect of high cardinality of variables on the performance of different machine learning
algorithms, and how a quick fix to reduce the number of labels, without any sort of data insight,
already helps to boost performance.

In this example:
We will:
- Learn how to quantify cardinality
- See examples of high and low cardinality variables
- Understand the effect of cardinality when preparing train and test sets
- Visualise the effect of cardinality on Machine Learning Model performance

We will use the Titanic dataset.
"""
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import AdaBoostClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split

# let's load the titanic dataset
data = pd.read_csv('dataset/titanic.csv')

data.head()
"""
   pclass  survived                                             name     sex  \
0       1         1                    Allen, Miss. Elisabeth Walton  female
1       1         1                   Allison, Master. Hudson Trevor    male
2       1         0                     Allison, Miss. Helen Loraine  female
3       1         0             Allison, Mr. Hudson Joshua Creighton    male
4       1         0  Allison, Mrs. Hudson J C (Bessie Waldo Daniels)  female

       age  sibsp  parch  ticket      fare cabin embarked boat   body  \
0  29.0000      0      0   24160  211.3375    B5        S    2    NaN
1   0.9167      1      2  113781  151.5500   C22        S   11    NaN
2   2.0000      1      2  113781  151.5500   C22        S  NaN    NaN
3  30.0000      1      2  113781  151.5500   C22        S  NaN  135.0
4  25.0000      1      2  113781  151.5500   C22        S  NaN    NaN

                         home.dest
0                     St Louis, MO
1  Montreal, PQ / Chesterville, ON
2  Montreal, PQ / Chesterville, ON
3  Montreal, PQ / Chesterville, ON
4  Montreal, PQ / Chesterville, ON """

"""
The categorical variables in this dataset are Name, Sex, Ticket, Cabin and Embarked.

Note that Ticket and Cabin contain both letters and numbers, so they could be treated as Mixed Variables.
For this example, I will treat them as categorical.
 """

# let's inspect the cardinality, this is the number
# of different labels, for the different categorical variables

print('Number of categories in the variable Name: {}'.format(len(data.name.unique())))

print('Number of categories in the variable Gender: {}'.format(len(data.sex.unique())))

print('Number of categories in the variable Ticket: {}'.format(len(data.ticket.unique())))

print('Number of categories in the variable Cabin: {}'.format(len(data.cabin.unique())))

print('Number of categories in the variable Embarked: {}'.format(len(data.embarked.unique())))

"""
Number of categories in the variable Name: 1307
Number of categories in the variable Gender: 2
Number of categories in the variable Ticket: 929
Number of categories in the variable Cabin: 182
Number of categories in the variable Embarked: 4 """

print('Total number of passengers in the Titanic: {}'.format(len(data)))
# Total number of passengers in the Titanic: 1309

"""
- While the variable Sex contains only 2 categories and
- Embarked 4 (low cardinality),
- the variables Ticket, Name and Cabin, as expected, contain a huge number of different labels (high cardinality).

To demonstrate the effect of high cardinality in train and test sets and machine learning performance,
- I will work with the variable Cabin.
- I will create a new variable with reduced cardinality.
"""

# let's explore the values / categories of Cabin

# we know from the previous outcome that there are 148
# different cabins, therefore the variable
# is highly cardinal

data.cabin.unique()
""" array(['B5', 'C22', 'E12', 'D7', 'A36', 'C101', nan, 'C62', 'B35', 'A23',
       'B58', 'D15', 'C6', 'D35', 'C148', 'C97', 'B49', 'C99', 'C52', 'T',
       'A31', 'C7', 'C103', 'D22', 'E33', 'A21', 'B10', 'B4', 'E40',
       'B38', 'E24', 'B51', 'B96', 'C46', 'E31', 'E8', 'B61', 'B77', 'A9',
       'C89', 'A14', 'E58', 'E49', 'E52', 'E45', 'B22', 'B26', 'C85',
       'E17', 'B71', 'B20', 'A34', 'C86', 'A16', 'A20', 'A18', 'C54',
       'C45', 'D20', 'A29', 'C95', 'E25', 'C111', 'C23', 'E36', 'D34',
       'D40', 'B39', 'B41', 'B102', 'C123', 'E63', 'C130', 'B86', 'C92',
       'A5', 'C51', 'B42', 'C91', 'C125', 'D10', 'B82', 'E50', 'D33',
       'C83', 'B94', 'D49', 'D45', 'B69', 'B11', 'E46', 'C39', 'B18',
       'D11', 'C93', 'B28', 'C49', 'B52', 'E60', 'C132', 'B37', 'D21',
       'D19', 'C124', 'D17', 'B101', 'D28', 'D6', 'D9', 'B80', 'C106',
       'B79', 'C47', 'D30', 'C90', 'E38', 'C78', 'C30', 'C118', 'D36',
       'D48', 'D47', 'C105', 'B36', 'B30', 'D43', 'B24', 'C2', 'C65',
       'B73', 'C104', 'C110', 'C50', 'B3', 'A24', 'A32', 'A11', 'A10',
       'B57', 'C28', 'E44', 'A26', 'A6', 'A7', 'C31', 'A19', 'B45', 'E34',
       'B78', 'B50', 'C87', 'C116', 'C55', 'D50', 'E68', 'E67', 'C126',
       'C68', 'C70', 'C53', 'B19', 'D46', 'D37', 'D26', 'C32', 'C80',
       'C82', 'C128', 'E39', 'D', 'F4', 'D56', 'F33', 'E101', 'E77', 'F2',
       'D38', 'F', 'E121', 'E10', 'G6', 'F38'], dtype=object) """

"""
- Let's now reduce the cardinality of the variable.

How?
- instead of using the entire cabin value, I will capture only the first letter.
- Rationale: the first letter indicates the deck on which the cabin was located, and is therefore an indication of
    both social class status and proximity to the surface of the Titanic.
    Both are known to improve the probability of survival.
 """
# let's capture the first letter of Cabin
data['Cabin_reduced'] = data['cabin'].astype(str).str[0]

data[['cabin', 'Cabin_reduced']].head()
""" cabin	Cabin_reduced
0	B5	B
1	C22	C
2	C22	C
3	C22	C
4	C22	C """

print('Number of categories in the variable Cabin: {}'.format(len(data.cabin.unique())))

print('Number of categories in the variable Cabin reduced: {}'.format(len(data.Cabin_reduced.unique())))

"""
Number of categories in the variable Cabin: 182
Number of categories in the variable Cabin reduced: 9

- We reduced the number of different labels from 182 to 9.
 """

# let's separate into training and testing set
# in order to build machine learning models

use_cols = ['cabin', 'Cabin_reduced', 'sex']

# this functions comes from scikit-learn
X_train, X_test, y_train, y_test = train_test_split(
    data[use_cols],
    data['survived'],
    test_size=0.3,
    random_state=0)

X_train.shape, X_test.shape
# ((916, 3), (393, 3))

"""
High cardinality leads to uneven distribution of categories in train and test sets
When a variable is highly cardinal, often some categories land only on the training set, or only on the testing set.
- If present only in the training set, they may lead to over-fitting.
- If present only on the testing set, the machine learning algorithm will not know how to handle them,
    as it has not seen them during training. """

# Let's find out labels present only in the training set

unique_to_train_set = [x for x in X_train.cabin.unique() if x not in X_test.cabin.unique()]

len(unique_to_train_set)
# 113

# There are 113 Cabins only present in the training set, and not in the testing set.

# Let's find out labels present only in the test set

unique_to_test_set = [x for x in X_test.cabin.unique() if x not in X_train.cabin.unique()]

len(unique_to_test_set)
# 36

"""
- Variables with high cardinality tend to have values (i.e., categories) present in the training set,
    that are not present in the test set, and vice versa.
- This will bring problems at the time of training (due to over-fitting) and scoring of new data
    (how should the model deal with unseen categories?).

This problem is almost overcome by reducing the cardinality of the variable. See below.
 """
# Let's find out labels present only in the training set
# for Cabin with reduced cardinality

unique_to_train_set = [
    x for x in X_train['Cabin_reduced'].unique()
    if x not in X_test['Cabin_reduced'].unique()
]

len(unique_to_train_set)
# 1

# Let's find out labels present only in the test set
# for Cabin with reduced cardinality

unique_to_test_set = [
    x for x in X_test['Cabin_reduced'].unique()
    if x not in X_train['Cabin_reduced'].unique()
]

len(unique_to_test_set)
# 0
"""
- Observe how by reducing the cardinality there is now only 1 label in the training set that is not present in the test set.
- And no label in the test set that is not contained in the training set as well.

Effect of cardinality on Machine Learning Model Performance
- In order to evaluate the effect of categorical variables in machine learning models,
- I will quickly replace the categories by numbers. See below.
 """
# Let's re-map Cabin into numbers so we can use it to train ML models

# I will replace each cabin by a number
# to quickly demonstrate the effect of
# labels on machine learning algorithms

##############
# Note: this is neither the only nor the best
# way to encode categorical variables into numbers
# there is more on these techniques in the section
# "Encoding categorical variales"
##############

cabin_dict = {k: i for i, k in enumerate(X_train.cabin.unique(), 0)}
cabin_dict
""" {nan: 0,
 'E36': 1,
 'C68': 2,
 'E24': 3,
 'C22': 4,
 'D38': 5,
 'B50': 6,
 'A24': 7,
 'C111': 8,
 'F': 9,
 'C6': 10,
 'C87': 11,
 'E8': 12,
 'B45': 13,
 'C93': 14,
 'D28': 15,
 'D36': 16,
 'C125': 17,
 'B35': 18,
 'T': 19,
 'B73': 20,
 'B57': 21,
 'A26': 22,
 'A18': 23,
 'B96': 24,
 'G6': 25,
 'C78': 26,
 'C101': 27,
 'D9': 28,
 'D33': 29,
 'C128': 30,
 'E50': 31,
 'B26': 32,
 'B69': 33,
 'E121': 34,
 'C123': 35,
 'B94': 36,
 'A34': 37,
 'D': 38,
 'C39': 39,
 'D43': 40,
 'E31': 41,
 'B5': 42,
 'D17': 43,
 'F33': 44,
 'E44': 45,
 'D7': 46,
 'A21': 47,
 'D34': 48,
 'A29': 49,
 'D35': 50,
 'A11': 51,
 'B51': 52,
 'D46': 53,
 'E60': 54,
 'C30': 55,
 'D26': 56,
 'E68': 57,
 'A9': 58,
 'B71': 59,
 'D37': 60,
 'F2': 61,
 'C55': 62,
 'C89': 63,
 'C124': 64,
 'C23': 65,
 'C126': 66,
 'E49': 67,
 'E46': 68,
 'D19': 69,
 'B58': 70,
 'C82': 71,
 'B52': 72,
 'C92': 73,
 'E45': 74,
 'C65': 75,
 'E25': 76,
 'B3': 77,
 'D40': 78,
 'C91': 79,
 'B102': 80,
 'B61': 81,
 'A20': 82,
 'B36': 83,
 'C7': 84,
 'B77': 85,
 'D20': 86,
 'C148': 87,
 'C105': 88,
 'E38': 89,
 'B86': 90,
 'C132': 91,
 'C86': 92,
 'A14': 93,
 'C54': 94,
 'A5': 95,
 'B49': 96,
 'B28': 97,
 'B24': 98,
 'C2': 99,
 'F4': 100,
 'A6': 101,
 'C83': 102,
 'B42': 103,
 'A36': 104,
 'C52': 105,
 'D56': 106,
 'C116': 107,
 'B19': 108,
 'E77': 109,
 'E101': 110,
 'B18': 111,
 'C95': 112,
 'D15': 113,
 'E33': 114,
 'B30': 115,
 'D21': 116,
 'E10': 117,
 'C130': 118,
 'D6': 119,
 'C51': 120,
 'D30': 121,
 'E67': 122,
 'C110': 123,
 'C103': 124,
 'C90': 125,
 'C118': 126,
 'C97': 127,
 'D47': 128,
 'E34': 129,
 'B4': 130,
 'D50': 131,
 'C62': 132,
 'E17': 133,
 'B41': 134,
 'C49': 135,
 'C85': 136,
 'B20': 137,
 'C28': 138,
 'E63': 139,
 'C99': 140,
 'D49': 141,
 'A10': 142,
 'A16': 143,
 'B37': 144,
 'C80': 145,
 'B78': 146}

 """

# replace the labels in Cabin, using the dic created above
X_train.loc[:, 'Cabin_mapped'] = X_train.loc[:, 'cabin'].map(cabin_dict)
X_test.loc[:, 'Cabin_mapped'] = X_test.loc[:, 'cabin'].map(cabin_dict)

X_train[['Cabin_mapped', 'cabin']].head(10)
"""
Cabin_mapped	cabin
501	0	NaN
588	0	NaN
402	0	NaN
1193	0	NaN
686	0	NaN
971	0	NaN
117	1	E36
540	0	NaN
294	2	C68
261	3	E24 """

# We see how NaN takes the value 0 in the new variable, E36 takes the value 1, C68 takes the value 2, and so on.
# Now I will replace the letters in the reduced cabin variable with the same procedure

# create replace dictionary
cabin_dict = {k: i for i, k in enumerate(X_train['Cabin_reduced'].unique(), 0)}

# replace labels by numbers with dictionary
X_train.loc[:, 'Cabin_reduced'] = X_train.loc[:, 'Cabin_reduced'].map(cabin_dict)
X_test.loc[:, 'Cabin_reduced'] = X_test.loc[:, 'Cabin_reduced'].map(cabin_dict)

X_train[['Cabin_reduced', 'cabin']].head(20)
"""
    Cabin_reduced	cabin
501	0	NaN
588	0	NaN
402	0	NaN
1193	0	NaN
686	0	NaN
971	0	NaN
117	1	E36
540	0	NaN
294	2	C68
261	1	E24
587	0	NaN
489	0	NaN
2	2	C22
405	0	NaN
1284	0	NaN
338	0	NaN
356	0	NaN
985	0	NaN
182	0	NaN
1027	0	NaN """

# We see now that E36 and E24 take the same number, 1, because we are capturing only the letter.
# They both start with E.

# re-map the categorical variable Sex into numbers

X_train.loc[:, 'sex'] = X_train.loc[:, 'sex'].map({'male': 0, 'female': 1})
X_test.loc[:, 'sex'] = X_test.loc[:, 'sex'].map({'male': 0, 'female': 1})

X_train.sex.head()
"""
501     1
588     1
402     1
1193    0
686     1
Name: sex, dtype: int64 """

# check if there are missing values in these variables

X_train[['Cabin_mapped', 'Cabin_reduced', 'sex']].isnull().sum()
"""
Cabin_mapped     0
Cabin_reduced    0
sex              0
dtype: int64 """

X_test[['Cabin_mapped', 'Cabin_reduced', 'sex']].isnull().sum()
"""
Cabin_mapped     41
Cabin_reduced     0
sex               0
dtype: int64 """

# In the test set, there are now 30 missing values for the highly cardinal variable.
# These were introduced when encoding the categories into numbers.
#  How?
# Many categories exist only in the test set. Thus, when we created our encoding dictionary using only the train set,
#   we did not generate a number to replace those labels present only in the test set.
#   As a consequence, they were encoded as NaN. We will see in future examples how to tackle this problem.
#   For now, I will fill those missing values with -1.

# let's check the number of different categories in the encoded variables
len(X_train.Cabin_mapped.unique()), len(X_train.Cabin_reduced.unique())
# (147, 9)

""" From the above we note immediately that from the original 148 cabins in the dataset,
    only 121 are present in the training set. We also see how we reduced the number of different categories to
    just 9 in our previous step.

Let's go ahead and evaluate the effect of labels in machine learning algorithms.
 """

# Random Forests
# model built on data with high cardinality for cabin

# call the model
rf = RandomForestClassifier(n_estimators=200, random_state=39)

# train the model
rf.fit(X_train[['Cabin_mapped', 'sex']], y_train)

# make predictions on train and test set
pred_train = rf.predict_proba(X_train[['Cabin_mapped', 'sex']])
pred_test = rf.predict_proba(X_test[['Cabin_mapped', 'sex']].fillna(0))

print('Train set')
print('Random Forests roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Random Forests roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Random Forests roc-auc: 0.853790650048556
Test set
Random Forests roc-auc: 0.7691361097284443

- We observe that the performance of the Random Forests on the training set is quite superior to its performance in
 the test set. This indicates that the model is over-fitting, which means that it does a great job at predicting
 the outcome on the dataset it was trained on, but it lacks the power to generalise the prediction to unseen data.
 """
# model built on data with low cardinality for cabin

# call the model
rf = RandomForestClassifier(n_estimators=200, random_state=39)

# train the model
rf.fit(X_train[['Cabin_reduced', 'sex']], y_train)

# make predictions on train and test set
pred_train = rf.predict_proba(X_train[['Cabin_reduced', 'sex']])
pred_test = rf.predict_proba(X_test[['Cabin_reduced', 'sex']])

print('Train set')
print('Random Forests roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Random Forests roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Random Forests roc-auc: 0.8163420365403872
Test set
Random Forests roc-auc: 0.8017670482827277

- We can see now that the Random Forests no longer over-fit to the training set.
-  In addition, the model is much better at generalising the predictions
        (compare the roc-auc of this model on the test set vs the roc-auc of the model above also in the test set: 0.81 vs 0.80).

I would like to point out, that likely we can overcome the effect of high cardinality by adjusting the hyper-parameters
of the random forests. That goes beyond the scope. Here, I want to show you that given a same model,
with identical hyper-parameters, high cardinality may cause the model to over-fit.
 """

# AdaBoost
# model build on data with plenty of categories in Cabin

# call the model
ada = AdaBoostClassifier(n_estimators=200, random_state=44)

# train the model
ada.fit(X_train[['Cabin_mapped', 'sex']], y_train)

# make predictions on train and test set
pred_train = ada.predict_proba(X_train[['Cabin_mapped', 'sex']])
pred_test = ada.predict_proba(X_test[['Cabin_mapped', 'sex']].fillna(0))

print('Train set')
print('Adaboost roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Adaboost roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Adaboost roc-auc: 0.8296861713101102
Test set
Adaboost roc-auc: 0.7604391350035948 """

# model build on data with fewer categories in Cabin Variable

# call the model
ada = AdaBoostClassifier(n_estimators=200, random_state=44)

# train the model
ada.fit(X_train[['Cabin_reduced', 'sex']], y_train)

# make predictions on train and test set
pred_train = ada.predict_proba(X_train[['Cabin_reduced', 'sex']])
pred_test = ada.predict_proba(X_test[['Cabin_reduced', 'sex']].fillna(0))

print('Train set')
print('Adaboost roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Adaboost roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Adaboost roc-auc: 0.8161256723642566
Test set
Adaboost roc-auc: 0.8001078480172557

- Similarly, the Adaboost model trained on the variable with high cardinality is overfit to the train set.
- Whereas the Adaboost trained on the low cardinal variable is not overfitting and therefore does a better job in
     generalising the predictions.

In addition, building an AdaBoost on a model with less categories in Cabin, is a) simpler and b) should a different category in the test set appear, by taking just the front letter of cabin, the ML model will know how to handle it because it was seen during training.
 """

# Logistic Regression
# model build on data with plenty of categories in Cabin variable

# call the model
logit = LogisticRegression(random_state=44, solver='lbfgs')

# train the model
logit.fit(X_train[['Cabin_mapped', 'sex']], y_train)

# make predictions on train and test set
pred_train = logit.predict_proba(X_train[['Cabin_mapped', 'sex']])
pred_test = logit.predict_proba(X_test[['Cabin_mapped', 'sex']].fillna(0))

print('Train set')
print('Logistic regression roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Logistic regression roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Logistic regression roc-auc: 0.8133909298124677
Test set
Logistic regression roc-auc: 0.7750815773463858 """

# model build on data with fewer categories in Cabin Variable

# call the model
logit = LogisticRegression(random_state=44, solver='lbfgs')

# train the model
logit.fit(X_train[['Cabin_reduced', 'sex']], y_train)

# make predictions on train and test set
pred_train = logit.predict_proba(X_train[['Cabin_reduced', 'sex']])
pred_test = logit.predict_proba(X_test[['Cabin_reduced', 'sex']].fillna(0))

print('Train set')
print('Logistic regression roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Logistic regression roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Logistic regression roc-auc: 0.8123468468695123
Test set
Logistic regression roc-auc: 0.8008268347989602

We can draw the same conclusion for Logistic Regression:
- reducing the cardinality improves the performance and generalisation of the algorithm.
 """

# Gradient Boosted Classifier
# model build on data with plenty of categories in Cabin variable

# call the model
gbc = GradientBoostingClassifier(n_estimators=300, random_state=44)

# train the model
gbc.fit(X_train[['Cabin_mapped', 'sex']], y_train)

# make predictions on train and test set
pred_train = gbc.predict_proba(X_train[['Cabin_mapped', 'sex']])
pred_test = gbc.predict_proba(X_test[['Cabin_mapped', 'sex']].fillna(0))

print('Train set')
print('Gradient Boosted Trees roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Gradient Boosted Trees roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))
"""
Train set
Gradient Boosted Trees roc-auc: 0.862631390919749
Test set
Gradient Boosted Trees roc-auc: 0.7733117637298823 """

# model build on data with plenty of categories in Cabin variable

# call the model
gbc = GradientBoostingClassifier(n_estimators=300, random_state=44)

# train the model
gbc.fit(X_train[['Cabin_reduced', 'sex']], y_train)

# make predictions on train and test set
pred_train = gbc.predict_proba(X_train[['Cabin_reduced', 'sex']])
pred_test = gbc.predict_proba(X_test[['Cabin_reduced', 'sex']].fillna(0))

print('Train set')
print('Gradient Boosted Trees roc-auc: {}'.format(roc_auc_score(y_train, pred_train[:,1])))
print('Test set')
print('Gradient Boosted Trees roc-auc: {}'.format(roc_auc_score(y_test, pred_test[:,1])))

"""
Train set
Gradient Boosted Trees roc-auc: 0.816719415917359
Test set
Gradient Boosted Trees roc-auc: 0.8015181682429069

- Gradient Boosted trees are indeed over-fitting to the training set in those cases where the variable Cabin has
    a lot of labels. This was expected as tree methods tend to be biased to variables with plenty of categories.

That is all for this example. I hope you enjoyed the information, and see you in the next one. """

ext: py

max_stars_repo_name: Akshaykumarcp/ML-Feature-Engineering

max_issues_repo_name: Akshaykumarcp/ML-Feature-Engineering

max_forks_repo_name: Akshaykumarcp/ML-Feature-Engineering

reason: import numpy


text:
from sklearn.metrics import jaccard_score
from sklearn.metrics import f1_score
from sklearn.metrics import accuracy_score

from sklearn.linear_model import PassiveAggressiveClassifier
import numpy as np

def _passiveaggressiveclassifier(*, train, test, x_predict=None, metrics, C=1.0, fit_intercept=True, max_iter=1000,
tol=0.001, early_stopping=False, validation_fraction=0.1, n_iter_no_change=5, shuffle=True,
verbose=0, loss='hinge', n_jobs=None, random_state=None, warm_start=False, class_weight=None, average=False):
    """For for info visit :
        https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.PassiveAggressiveClassifier.html#sklearn.linear_model.PassiveAggressiveClassifier
    """

    model = PassiveAggressiveClassifier(C=C, fit_intercept=fit_intercept, max_iter=max_iter, tol=tol, early_stopping=early_stopping,
validation_fraction=validation_fraction, n_iter_no_change=n_iter_no_change, shuffle=shuffle, verbose=verbose, loss=loss,
 n_jobs=n_jobs, random_state=random_state, warm_start=warm_start, class_weight=class_weight, average=average)
    model.fit(train[0], train[1])
    model_name = 'PassiveAggressiveClassifier'
    y_hat = model.predict(test[0])

    if metrics == 'f1_score':
        accuracy = f1_score(test[1], y_hat)
    if metrics == 'jaccard_score':
        accuracy = jaccard_score(test[1], y_hat)
    if metrics == 'accuracy_score':
        accuracy = accuracy_score(test[1], y_hat)

    if x_predict is None:
        return (model_name, accuracy, None)

    y_predict = model.predict(x_predict)
    return (model_name, accuracy, y_predict)

max_stars_repo_name: kritikaparmar-programmer/DataScientist

max_issues_repo_name: kritikaparmar-programmer/DataScientist

max_forks_repo_name: kritikaparmar-programmer/DataScientist


text:
import os.path as op
import re
import sys
from toolshed import reader, nopen
from collections import defaultdict

import numpy as np
from itertools import cycle
import pylab as pl
import seaborn

colors = cycle(seaborn.color_palette('Set1', 8))

BASES = False

def name(bam):
    return op.basename(bam).rsplit(".", 1)[0].split("-")[0]
    return op.basename(bam).rsplit(".", 1)[0].split("-")[0] + ("-trim" if "trim" in bam else "")
    return op.dirname(bam)

def count_bases(cigar, patt=re.compile("\d+[IM]")):
    return sum(int(p[:-1]) for p in patt.findall(cigar))

def count_on_off(bam, flags, pad, bases=BASES):

    #if not "last" in bam: flags = flags + " -f 0x2"

    on_count = [0] * 256
    off_count = [0] * 256
    rcounts = defaultdict(int)

    posn = re.compile(".*?_(.+):(\d+)-(\d+)")

    # chr1b:3001315:+__chr1b:3001467:-        99      chr1    3001316 60 100M
    print >>sys.stderr, bam
    for toks in reader("|samtools view {flags} {bam}".format(**locals()),
            header=False):
        if toks[0][0] == "@": continue
        if toks[0].endswith(("_R1", "_R2")):
            toks[0] = toks[0][:-3]

        cigar = toks[5]

        qual = int(toks[4])
        if toks[0].endswith("_BAD"):
            off_count[qual] += (count_bases(cigar) if bases else 1)
            continue

        rname, flag = toks[0], int(toks[1])
        rcounts[rname] += 1
        # @2_chr4:72040172-72040686_R1
        if rcounts[rname] > 2:
            raise Exception("%s:%s", (bam, rname))

        try:
            chrom, start, end = re.match(posn, rname).groups(0)
        except:
            print >>sys.stderr, rname
            raise

        pos = int(start if flag & 0x40 else end)
        on = chrom == toks[2] and abs(pos - int(toks[3])) <= pad
        if on:
            on_count[qual] += (count_bases(cigar) if bases else 1)
        else:
            off_count[qual] += (count_bases(cigar) if bases else 1)

    on_count = np.cumsum(on_count[::-1])[::-1]
    off_count = np.cumsum(off_count[::-1])[::-1]
    return off_count, on_count

FLAGS="-F%i" % (0x4 | 0x100 | 0x200)
def main(bams, reads=None, flags=FLAGS, pad=12002):
    if not reads.isdigit():
        reads = 2 * float(nopen("|bioawk -c fastx 'END { print NR }' %s" % reads).next())
    else:
        reads = 2.0 * int(reads)

    if any('trim' in b for b in bams):
        assert all('trim' in b for b in bams), [b for b in bams if not 'trim'
                in b]
    counts = {}
    names, pts = [], []

    denom = float(reads)
    if BASES: denom *= 100.
    for bam in bams:
        counts[bam] = count_on_off(bam, flags, pad)

    out = sys.stdout
    print >>out, "qual\tmethod\toff\ton"

    for qual in range(0, 256):
        for b in bams:
            print >>out, "{qual}\t{bam}\t{off}\t{on}".format(
            qual=qual, bam=name(b),
            off=counts[b][0][qual] / denom,
            on=counts[b][1][qual] / denom)
    print >>sys.stderr, "wrote", out.name

    for bam in bams:

        # multiply by 100 bases per read.
        color = next(colors)

        symbol = 'o' if len(set(counts[bam][0])) < 3 else '.'
        nmax = max(i for i, v in enumerate(counts[bam][1]) if v > 0)
        p, = pl.plot(counts[bam][0][1:nmax + 1] / denom,
                counts[bam][1][1:nmax + 1] / denom,
                color=color,
                mec=color,
                mfc=color,
                marker=symbol, label=name(bam))
        pts.append(p)
        names.append(name(bam))
    pl.xlabel('off target')
    pl.ylabel('on target')
    pl.legend(pts, names, loc='lower right')
    pl.xlim(xmin=0)
    pl.ylim(ymin=0)
    print >>sys.stderr, reads

    pl.show()

if __name__ == "__main__":
    import argparse
    p = argparse.ArgumentParser()
    p.add_argument("--reads", help="reads file", required=True)
    p.add_argument("bams", nargs="+")

    a = p.parse_args()

    main(a.bams, reads=a.reads)


text:
##########################################################################
# NSAp - Copyright (C) CEA, 2017
# Distributed under the terms of the CeCILL-B license, as published by
# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
# for details.
##########################################################################

"""
Common tools to filter tractograms.
"""

# System import
from __future__ import print_function
import os

# Third party import
import nibabel
from nibabel.streamlines.tractogram import Tractogram
import numpy
from dipy.viz.colormap import line_colors
from dipy.viz import fvtk
from dipy.core.gradients import gradient_table
from dipy.tracking.distances import approx_polygon_track
from dipy.tracking.utils import length
import dipy.core.optimize as opt
import dipy.tracking.life as dpilife
import matplotlib.pyplot as plt
import matplotlib

def life(dwifile, bvecsfile, bvalsfile, tractogramfile, outdir,
         display_tracks=False, verbose=0):
    """ Linear fascicle evaluation (LiFE)
    Evaluating the results of tractography algorithms is one of the biggest
    challenges for diffusion MRI. One proposal for evaluation of tractography
    results is to use a forward model that predicts the signal from each of a
    set of streamlines, and then fit a linear model to these simultaneous
    prediction.

    Parameters
    ----------
    dwifile: str
        the path to the diffusion dataset.
    bvecsfile: str
        the path to the diffusion b-vectors.
    bvalsfile: str
        the path to the diffusion b-values.
    tractogramfile: str
        the path to the tractogram.
    outdir: str
        the destination folder.
    display_tracks: bool, default False
        if True render the tracks.
    verbose: int, default 0
        the verbosity level.

    Returns
    -------
    life_weights_file: str
        a file containing the fiber track weights.
    life_weights_snap: str
        a snap with the distrubution of weights.
    spatial_error_file: str
        the model root mean square error.
    tracks_snap: str
        a snap with the tracks.
    """
    # Load diffusion data and tractogram
    bvecs = numpy.loadtxt(bvecsfile)
    bvals = numpy.loadtxt(bvalsfile)
    gtab = gradient_table(bvals, bvecs)
    im = nibabel.load(dwifile)
    data = im.get_data()
    trk = nibabel.streamlines.load(tractogramfile)
    if verbose > 0:
        print("[info] Diffusion shape: {0}".format(data.shape))
        print("[info] Number of tracks: {0}".format(len(trk.streamlines)))

    # Express the tractogram in voxel coordiantes
    inv_affine = numpy.linalg.inv(trk.affine)
    trk = [
        numpy.dot(
            numpy.concatenate(
                (
                    streamline,
                    numpy.ones(
                        (len(streamline), 1)
                    )
                ), axis=1
            ),
            inv_affine) for streamline in trk.streamlines]
    trk = [track[..., :3] for track in trk]

    # Create a viewer
    tracks_snap = None
    if display_tracks:
        nb_tracks = len(trk)
        if nb_tracks < 5000:
            downsampling = 1
        else:
            downsampling = nb_tracks // 5000
        tracks_snap = os.path.join(outdir, "tracks.png")
        streamlines_actor = fvtk.line(trk[::downsampling],
                                      line_colors(trk[::downsampling]))
        vol_actor = fvtk.slicer(data[..., 0])

        vol_actor.display(data.shape[0] // 2, None, None)
        ren = fvtk.ren()
        fvtk.add(ren, streamlines_actor)
        fvtk.add(ren, vol_actor)
        fvtk.record(ren, n_frames=1, out_path=tracks_snap, size=(800, 800))
        if verbose > 1:
            fvtk.show(ren)

    # Fit the Life model and save the associated weights
    fiber_model = dpilife.FiberModel(gtab)
    fiber_fit = fiber_model.fit(data, trk, affine=numpy.eye(4))
    life_weights = fiber_fit.beta
    life_weights_file = os.path.join(outdir, "life_weights.txt")
    numpy.savetxt(life_weights_file, life_weights)
    life_weights_snap = os.path.join(outdir, "life_weights.png")
    fig, ax = plt.subplots(1)
    ax.hist(life_weights, bins=100, histtype="step")
    ax.set_xlabel("Fiber weights")
    ax.set_ylabel("# Fibers")
    fig.savefig(life_weights_snap)

    # Invert the model and predict back either the data that was used to fit
    # the model: compute the prediction error of the diffusion-weighted data
    # and calculate the root of the mean squared error.
    model_predict = fiber_fit.predict()
    model_error = model_predict - fiber_fit.data
    model_rmse = numpy.sqrt(numpy.mean(model_error[:, 10:] ** 2, -1))
    data_rmse = numpy.ones(data.shape[:3]) * numpy.nan
    data_rmse[fiber_fit.vox_coords[:, 0],
              fiber_fit.vox_coords[:, 1],
              fiber_fit.vox_coords[:, 2]] = model_rmse
    model_error_file = os.path.join(outdir, "model_rmse.nii.gz")
    error_im = nibabel.Nifti1Image(data_rmse, im.affine)
    nibabel.save(error_im, model_error_file)

    # As a baseline against which we can compare, we assume that the weight
    # for each streamline is equal to zero. This produces the naive prediction
    # of the mean of the signal in each voxel.
    life_weights_baseline = numpy.zeros(life_weights.shape[0])
    pred_weighted = numpy.reshape(
        opt.spdot(fiber_fit.life_matrix, life_weights_baseline),
        (fiber_fit.vox_coords.shape[0], numpy.sum(~gtab.b0s_mask)))
    mean_pred = numpy.empty(
        (fiber_fit.vox_coords.shape[0], gtab.bvals.shape[0]))
    S0 = fiber_fit.b0_signal

    # Since the fitting is done in the demeaned S/S0 domain, we need to add
    # back the mean and then multiply by S0 in every voxels.
    mean_pred[..., gtab.b0s_mask] = S0[:, None]
    mean_pred[..., ~gtab.b0s_mask] = (
            (pred_weighted + fiber_fit.mean_signal[:, None]) * S0[:, None])
    mean_error = mean_pred - fiber_fit.data
    mean_rmse = numpy.sqrt(numpy.mean(mean_error ** 2, -1))
    data_rmse = numpy.ones(data.shape[:3]) * numpy.nan
    data_rmse[fiber_fit.vox_coords[:, 0],
              fiber_fit.vox_coords[:, 1],
              fiber_fit.vox_coords[:, 2]] = mean_rmse
    mean_error_file = os.path.join(outdir, "mean_rmse.nii.gz")
    error_im = nibabel.Nifti1Image(data_rmse, im.affine)
    nibabel.save(error_im, mean_error_file)

    # Compute the improvment array
    data_rmse = numpy.ones(data.shape[:3]) * numpy.nan
    data_rmse[fiber_fit.vox_coords[:, 0],
              fiber_fit.vox_coords[:, 1],
              fiber_fit.vox_coords[:, 2]] = mean_rmse - model_rmse
    improvment_error_file = os.path.join(outdir, "improvment_rmse.nii.gz")
    error_im = nibabel.Nifti1Image(data_rmse, im.affine)
    nibabel.save(error_im, improvment_error_file)

    return (life_weights_file, life_weights_snap, model_error_file,
            mean_error_file, improvment_error_file, tracks_snap)

def lossy_compression_of_tractogram(tractogramfile, outdir, rate=0.392,
                                    search_optimal_rate=False,
                                    weightsfile=None, weights_thr=0.,
                                    max_search_dist=2.2, verbose=0):
    """ Reduce the number of points of the track by keeping
    intact the start and endpoints of the track and trying to remove
    as many points as possible without distorting much the shape of
    the track, ie. more points in curvy regions and less points in less curvy
    regions.

    Parameters
    ----------
    tractogramfile: str
        the path to the tractogram.
    outdir: str
        the destination folder.
    rate: float, default 0.392
        the compression rate, ie. smoothing parameter (<0.392 smoother,
        >0.392 rougher).
    search_optimal_rate: bool, default False
        determine the optimal compression rate.
    weightsfile: str, default None
        use these weights to remove unsignificant streamlines.
    weights_thr: float, default 0.
        the threshold used to identify unsignificant streamlines.
    max_search_dist: float, default 2.2
        the maximum distance between the initial and downsampled streamlines
        allowed during the best rate search.
    verbose: int, default 0
        the verbosity level.

    Returns
    -------
    compressed_tractogramfile: str
        the compressed tractogram.
    nb_points_file: str
        the compression result compared to the original sampling.
    """
    # Load the tractogram
    trk = nibabel.streamlines.load(tractogramfile)
    if verbose > 0:
        print("[info] Number of tracks: {0}".format(len(trk.streamlines)))

    # Keep only significant streamlines
    tracks = trk.streamlines
    if weightsfile is not None:
        weights = numpy.loadtxt(weightsfile)
        keep_indices = numpy.where(weights > weights_thr)[0]
        tracks = list(numpy.array(tracks)[keep_indices])
        weights = weights[numpy.where(keep_indices)[0]]
        if verbose > 0:
            print("[info] Number of significant tracks: {0}".format(
                len(tracks)))
    else:
        weights = None

    # Compress tractogram
    # > dynamic compression rate
    if search_optimal_rate:
        rate = "dynamic"
        ref_lengths = list(length(tracks))
        rates = numpy.linspace(1, 0, 21)
        opt_lengths = numpy.zeros((len(rates), len(tracks)))
        for idx, optrate in enumerate(rates):
            if verbose > 0:
                print("[info] Grid search at rate '{0}'.".format(optrate))
            decimated_tracks = [
                approx_polygon_track(t, optrate) for t in tracks]
            opt_lengths[idx] = list(length(decimated_tracks))
            opt_lengths[idx] -= ref_lengths
        opt_lengths = numpy.abs(opt_lengths)
        if verbose > 2:
            print("[debug] Optimal lengths: {0}".format(opt_lengths))
        opt_lengths[numpy.where(opt_lengths > max_search_dist)] = 0
        opt_rate_indices = numpy.argmax(opt_lengths, axis=0)
        if verbose > 2:
            print("[debug] Optimal rate indices: {0}".format(opt_rate_indices))
        tracks = [approx_polygon_track(t, rates[i])
                  for t, i in zip(tracks, opt_rate_indices)]
    # > static compression rate
    else:
        tracks = [approx_polygon_track(t, rate) for t in tracks]
    compressed_tractogramfile = os.path.join(
        outdir, "compressed_tractogram.trk")
    compressed_trk = Tractogram(
        streamlines=tracks,
        affine_to_rasmm=trk.affine)
    nibabel.streamlines.save(compressed_trk, compressed_tractogramfile)

    # Summary graph
    n_pts_initial = [len(t) for t in trk.streamlines]
    n_pts_compressed = [len(t) for t in compressed_trk.streamlines]
    nb_points_file = os.path.join(outdir, "nb_points.png")
    fig, ax = plt.subplots(1)
    ax.hist(n_pts_initial, color="r", histtype="step", label="initial")
    ax.hist(n_pts_compressed, color="b", histtype="step",
            label="compressed ({0})".format(rate))
    ax.set_xlabel("Number of points")
    ax.set_ylabel("Count")
    plt.legend()
    plt.savefig(nb_points_file)

    return compressed_tractogramfile, nb_points_file

max_stars_repo_name: neurospin/pyconnectome

max_issues_repo_name: neurospin/pyconnectome

max_forks_repo_name: neurospin/pyconnectome


text:
''' this is a file which is intended to store the non class functions pertaining to the map '''

import graphics as g
import numpy as np
import random
import math

# creates 'num' amount of points with centers in ranges
# xrange_ and yrange_, which are tuples denoting the min and max
# x/y coordinates for each of the points
def generate_random_points(num, xrange_, yrange_):
	return np.array([ g.Point(random.randint(xrange_[0], xrange_[1]), random.randint(yrange_[0], yrange_[1])) for k in range(num) ])

# returns the straight line distance between two points
# assumes the points are of type Point from the graphics lib
def sl_distance(pt1, pt2):
	return math.sqrt(((pt2.x - pt1.x) ** 2) + ((pt2.y - pt1.y) ** 2))

# the same as above but doesn't use a sqrt to save computation time
# this makes sense in a more preformant language
def fast_sl_distance(pt1, pt2):
	return ((pt2.x - pt1.x) ** 2) + ((pt2.y - pt1.y) ** 2)

# gets a subarray of arr that starts at index a and goes until index b
# wrapping around if need be
# assumes both arrays are numpy arrays
def get_circular_range(arr, a, b):
	if a > b:
		return np.concatenate((arr[a:], arr[:b]))
	else:
		# this is done so the original arr is not mutated
		# it can lead to frustration if it is
		return np.array(list(arr[a:b]))

# takes in a numpy array and two points, reverses their order in the path,
# and places them back in the path
def swap_and_reverse(arr, a, b):
	# if a and b are the same, there is no swapping to be done
	# so just return
	if a == b:
		return arr
	# get the array by value not reference, we don't want to mutate the list
	_arr = np.array(list(arr))
	# get the subarray and reverse it
	subarray = get_circular_range(_arr, a, b)
	subarray = np.flip(subarray, 0)
	# now replace the indicies between a and b with the reversed subarray
	# different logic depending on whether or not a is smaller than b
	# basically the range differs
	if a < b:
		for k in range(abs(b-a)):
			_arr[(a + k) % len(_arr)] = subarray[k]
	else:
		for k in range(len(arr) - abs(b-a)):
			_arr[(a + k) % len(_arr)] = subarray[k]
	# return the array
	return _arr

# generates the random points and prints them to a .csv file
# also generates the number of points printed to the file
# assumes filename does not have an extension
# also generates a file in the same directory
def generate_random_points_to_file(num, xrange_, yrange_, filename):
	# get the random points
	pts = generate_random_points(num, xrange_, yrange_)
	# then print them to the file
	file = open(filename + ".csv", 'w')
	file.write("{}\n".format(num))
	for point in pts:
		file.write("{},{},\n".format(point.x, point.y))
	file.close()

# grabs points from a text file
# to be used to test algorithms on static inputs
def get_points_from_txt(filename):
	file = open(filename, 'r')
	num = int(file.readline())
	# initialize the array we wish to return
	arr = np.array([g.Point(0,0) for k in range(num)])
	count = 0
	# for each of the lines, parse them and store them in the array
	for line in file:
		arr[count] = g.Point(line.split(',')[0], line.split(',')[1])
		count += 1
	# close the file and return the array
	file.close()
	return arr


text:
import numpy as np
import pandas as pd
from typing import Tuple, List, Callable, Type

from plotly.subplots import make_subplots

from IMLearn import BaseModule
from IMLearn.desent_methods import GradientDescent, FixedLR, ExponentialLR
from IMLearn.desent_methods.modules import L1, L2
from IMLearn.learners.classifiers.logistic_regression import LogisticRegression
from IMLearn.metrics import misclassification_error
from IMLearn.model_selection import cross_validate
from IMLearn.utils import split_train_test
from sklearn.metrics import roc_curve, auc

import plotly.graph_objects as go

def plot_descent_path(module: Type[BaseModule],
                      descent_path: np.ndarray,
                      title: str = "",
                      xrange=(-1.5, 1.5),
                      yrange=(-1.5, 1.5)) -> go.Figure:
    """
    Plot the descent path of the gradient descent algorithm

    Parameters:
    -----------
    module: Type[BaseModule]
        Module type for which descent path is plotted

    descent_path: np.ndarray of shape (n_iterations, 2)
        Set of locations if 2D parameter space being the regularization path

    title: str, default=""
        Setting details to add to plot title

    xrange: Tuple[float, float], default=(-1.5, 1.5)
        Plot's x-axis range

    yrange: Tuple[float, float], default=(-1.5, 1.5)
        Plot's x-axis range

    Return:
    -------
    fig: go.Figure
        Plotly figure showing module's value in a grid of [xrange]x[yrange] over which regularization path is shown

    Example:
    --------
    fig = plot_descent_path(IMLearn.desent_methods.modules.L1, np.ndarray([[1,1],[0,0]]))
    fig.show()
    """

    def predict_(w):
        return np.array([module(weights=wi).compute_output() for wi in w])

    from utils import decision_surface
    return go.Figure([decision_surface(predict_, xrange=xrange, yrange=yrange, density=70, showscale=False),
                      go.Scatter(x=descent_path[:, 0], y=descent_path[:, 1], mode="markers+lines",
                                 marker_color="black")],
                     layout=go.Layout(xaxis=dict(range=xrange),
                                      yaxis=dict(range=yrange),
                                      title=f"GD Descent Path {title}"))

def get_gd_state_recorder_callback() -> Tuple[Callable[[], None], List[np.ndarray], List[np.ndarray]]:
    """
    Callback generator for the GradientDescent class, recording the objective's value and parameters at each iteration

    Return:
    -------
    callback: Callable[[], None]
        Callback function to be passed to the GradientDescent class, recoding the objective's value and parameters
        at each iteration of the algorithm

    values: List[np.ndarray]
        Recorded objective values

    weights: List[np.ndarray]
        Recorded parameters
    """
    values, weights = [], []

    def callback(gradient_descent, args):
        values.append(args[1])
        weights.append(args[0])

    return callback, values, weights

def compare_fixed_learning_rates(init: np.ndarray = np.array([np.sqrt(2), np.e / 3]),
                                 etas: Tuple[float] = (1, .1, .01, .001)):
    lowest_loss = []
    L1_convergence = go.Figure()
    L2_convergence = go.Figure()
    for eta in etas:
        for model, name in ((L1, "L1"), (L2, "L2")):
            cur_LR = FixedLR(eta)
            callback, values, weights = get_gd_state_recorder_callback()

            GD_solver = GradientDescent(cur_LR, callback=callback)
            GD_solver.fit(model(init.copy()), None, None)
            title = f"{name} module as function of iteration number (eta = {eta})".format(name=name, eta=eta)
            plot_descent_path(model, np.array(weights), title)  # .show()
            if eta == 0.01:
                plot_descent_path(model, np.array(weights), title).show()

            fig = go.Scatter(x=np.arange(GD_solver.max_iter_), y=values, mode='lines',
                             name=f"learning rate {eta}".format(name=name, eta=eta))
            if name == "L1":
                L1_convergence.add_trace(fig)
            else:
                L2_convergence.add_trace(fig)
            lowest_loss.append(np.min(values))
    L1_convergence.update_layout(title="L1 convergence rate as function of iteration number")
    L1_convergence.update_xaxes(title_text="Iteration Number")
    L1_convergence.update_yaxes(title_text=f"{name} Norm".format(name=name))
    L1_convergence.show()

    L2_convergence.update_layout(title="L2 convergence rate as function of iteration number")
    L2_convergence.update_xaxes(title_text="Iteration Number")
    L2_convergence.update_yaxes(title_text=f"{name} Norm".format(name=name))
    L2_convergence.show()

    print(lowest_loss)
    lowest_loss = np.array(lowest_loss)
    lowest_loss.shape = (len(etas), 2)
    for i in range(len(lowest_loss)):
        print("eta = ", etas[i], "L1: ", lowest_loss[i][0], "  L2: ", lowest_loss[i][1])

def compare_exponential_decay_rates(init: np.ndarray = np.array([np.sqrt(2), np.e / 3]),
                                    eta: float = .1,
                                    gammas: Tuple[float] = (.9, .95, .99, 1)):
    # Optimize the L1 objective using different decay-rate values of the exponentially decaying learning rate
    fig = make_subplots(rows=2, cols=2, subplot_titles=gammas, horizontal_spacing=0.01, vertical_spacing=0.05)
    figs = []
    for i, gamma in enumerate(gammas):
        cur_LR = ExponentialLR(eta, gamma)
        convergences = []
        callback, values, weights = get_gd_state_recorder_callback()
        GD_solver = GradientDescent(cur_LR, callback=callback)

        GD_solver.fit(L1(init.copy()), None, None)

        fig.add_trace(go.Scatter(x=np.arange(GD_solver.max_iter_), y=values, name=gamma),
                      row=(i // 2) + 1, col=(i % 2) + 1)
        print(f"min L1 with decay rate of {gamma} is {np.min(values)}")

    # Plot algorithm's convergence for the different values of gamma
    fig.update_layout(title="L1 Convergence rate as a function of iteration number")
    # fig.show()

    # Plot descent path for gamma=0.95
    L1_module = L1(init.copy())
    lr = ExponentialLR(eta, 0.95)
    callback, values, weights = get_gd_state_recorder_callback()
    GradientDescent(learning_rate=lr, callback=callback).fit(L1_module, None, None)
    fig = plot_descent_path(L1, np.array(weights), title="L1 norm with exponential learning rate gamma=0.95")
    fig.show()

def load_data(path: str = "../datasets/SAheart.data", train_portion: float = .8) -> \
        Tuple[pd.DataFrame, pd.Series, pd.DataFrame, pd.Series]:
    """
    Load South-Africa Heart Disease dataset and randomly split into a train- and test portion

    Parameters:
    -----------
    path: str, default= "../datasets/SAheart.data"
        Path to dataset

    train_portion: float, default=0.8
        Portion of dataset to use as a training set

    Return:
    -------
    train_X : DataFrame of shape (ceil(train_proportion * n_samples), n_features)
        Design matrix of train set

    train_y : Series of shape (ceil(train_proportion * n_samples), )
        Responses of training samples

    test_X : DataFrame of shape (floor((1-train_proportion) * n_samples), n_features)
        Design matrix of test set

    test_y : Series of shape (floor((1-train_proportion) * n_samples), )
        Responses of test samples
    """
    df = pd.read_csv(path)
    df.famhist = (df.famhist == 'Present').astype(int)
    return split_train_test(df.drop(['chd', 'row.names'], axis=1), df.chd, train_portion)

def fit_logistic_regression():
    # Load and split SA Heard Disease dataset
    X_train, y_train, X_test, y_test = load_data()
    X_train, y_train, X_test, y_test = (X_train - X_train.mean()) / X_train.std(), y_train, (
            X_test - X_test.mean()) / X_test.std(), y_test

    X_train, y_train = X_train.to_numpy(), y_train.to_numpy()
    X_test, y_test = X_test.to_numpy(), y_test.to_numpy()

    # # Plotting convergence rate of logistic regression over SA heart disease data
    callback, values, weights = get_gd_state_recorder_callback()
    gd = GradientDescent(callback=callback, learning_rate=FixedLR(1e-4), max_iter=20000)
    lg = LogisticRegression(solver=gd, include_intercept=True).fit(X_train, y_train)
    fpr, tpr, thresholds = roc_curve(y_train, lg.predict_proba(X_train))
    alphas = np.arange(0, 1, 0.01)

    fig = go.Figure(data=[go.Scatter(x=[0, 1], y=[0, 1], mode="lines", line=dict(color="black", dash='dash'),
                                     name="Random Class Assignment", showlegend=False),
                          go.Scatter(x=fpr, y=tpr, mode='markers+lines', text=alphas, name="", showlegend=False,
                                     marker_size=5,
                                     hovertemplate="Threshold:%{text:.3f}
FPR: %{x:.3f}
TPR: %{y:.3f}")],
                    layout=go.Layout(title=rf"$\text{{ROC Curve Of Fitted Model - AUC}}={auc(fpr, tpr):.6f}$",
                                     xaxis=dict(title=r"$\text{False Positive Rate (FPR)}$"),
                                     yaxis=dict(title=r"$\text{True Positive Rate (TPR)}$")))
    # fig.show()
    best_alpha_index = np.argmax(tpr - fpr)
    best_alpha_value = thresholds[best_alpha_index]
    print(f"Best alpha: {best_alpha_value}")
    error = misclassification_error(y_test, lg.predict_proba(X_test) >= best_alpha_value)
    print("Test error for alpha: ", error)

    # # Fitting l1- and l2-regularized logistic regression models, using cross-validation to specify values
    # of regularization parameter
    # alpha = 0.05
    # modules = ["l1", "l2"]
    # lambdas = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1]
    #
    # callback, values, weights = get_gd_state_recorder_callback()
    # for module in modules:
    #     cv_scores = []
    #     for lam in lambdas:
    #         lg = LogisticRegression(include_intercept=True, penalty=module, alpha=alpha, lam=lam,
    #                                 solver=GradientDescent(learning_rate=FixedLR(1e-4), max_iter=20000,
    #                                                        callback=callback))
    #         train_score, validation_score = cross_validate(lg, X_train, y_train, misclassification_error)
    #         cv_scores.append(validation_score)
    #
    #     best_lambda = lambdas[np.argmin(np.array(cv_scores))]
    #     print(f"best lambda for {module} penalty: {best_lambda}")
    #     callback, values, weights = get_gd_state_recorder_callback()
    #     lg = LogisticRegression(penalty=module, alpha=alpha, lam=best_lambda,
    #                             solver=GradientDescent(learning_rate=FixedLR(1e-4), max_iter=20000,
    #                                                    callback=callback)).fit(X_train, y_train)
    #     error = misclassification_error(y_test, lg.predict(X_test))
    #     print(f"Test error for {module} penalty and {best_lambda} is: {error}")
    lambdas = np.array([0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1])
    for i in range(1, 3):
        callback, values, weights = get_gd_state_recorder_callback()
        validation_scores = np.zeros(lambdas.shape[0])
        train_scores = np.zeros(lambdas.shape[0])
        for j, lam in enumerate(lambdas[0:5]):
            model = LogisticRegression(penalty=f"l{i}", lam=lam, solver=GradientDescent(learning_rate=FixedLR(1e-4), max_iter=20000,
                                                       callback=callback))
            train_score, validation_score = cross_validate(model,
                                                           X_train,
                                                           y_train,
                                                           misclassification_error)
            validation_scores[j] = validation_score
            train_scores[j] = train_score
        best_lambda = lambdas[np.argmin(validation_scores)]
        model = LogisticRegression(penalty=f"l{i}", lam=best_lambda, solver=GradientDescent(learning_rate=FixedLR(1e-4), max_iter=20000,
                                                       callback=callback))
        model.fit(X_train, y_train)
        loss = model.loss(X_test, y_test)
        print(f"Best lambda for model with regularization L{i} is: {loss} with lambda {best_lambda}")

if __name__ == '__main__':
    np.random.seed(0)
    # compare_fixed_learning_rates()
    # compare_exponential_decay_rates()
    fit_logistic_regression()


text:
from settings import *

def loadXMLAttributes(d, attributes):
    for attr in attributes.keys():
        sattr = str(attr)
        val = attributes[attr].value
        try:
            sval = str(val)
            if sval.isdigit():
                d[sattr] = int(sval)
            else:
                d[sattr] = sval
        except:
            d[sattr] = val
    return d

def getNoneTextChildNodes(xmlNode):
    """
    """

    return filter(lambda x:x.nodeType != x.TEXT_NODE, xmlNode.childNodes)

def testLU(fileName):
    a = LexicalUnit()

    try:
        if not a.loadXML(FRAMENET_PATH+"/"+LU_DIR_ENV+"/"+fileName):
            print >> sys.stderr, 'loading:', fileName, 'failed'
    except:
        print >> sys.stderr, 'loading:', fileName, 'failed'
    print a.getLexemes()
    return a

def initialize():
     """
     """

     fn = FrameNet()
     fn.initialize()

     pass

from lexical_unit import LexicalUnit
from frame import Frame
from frame_relation import FrameRelation

class FrameNet(dict):
    """
    This the master class for FrameNet
    """

    def __init__(self):
        """
        Constructor, loads the frames and the frame relations.
        """
        pickledFramePath = FRAMENET_PATH + '/' + PICKLED_FRAME_FILE
        pickledFrameRelationPath = FRAMENET_PATH  + '/' + PICKLED_FRAME_RELATIONS_FILE
        pickledLUPath = FRAMENET_PATH + '/' + PICKLED_LU_FILE

        try:
            print >> sys.stderr, 'Loading the frames ...',
            frames = cPickle.load(open(pickledFramePath))
            self['frames'] = frames
            print >> sys.stderr, 'done'

            print >> sys.stderr, 'Loading the frame relations ...',
            frRelations = cPickle.load(open(pickledFrameRelationPath))
            self['frameRelations'] = frRelations
            print >> sys.stderr, 'done'

            print >> sys.stderr, 'Loading the lexical units data ...',
            luData = cPickle.load(open(pickledLUPath))
            print >> sys.stderr, 'done'
            self['luIndex'] = luData

        except:
            print >> sys.stderr, 'Framenet not initialized, doing it now'
            self.initialize()
            print >> sys.stderr, 'Loading the frames ...',
            frames = cPickle.load(open(pickledFramePath))
            self['frames'] = frames
            print >> sys.stderr, 'done'

            print >> sys.stderr, 'Loading the frame relations ...',
            frRelations = cPickle.load(open(pickledFrameRelationPath))
            self['frameRelations'] = frRelations
            print >> sys.stderr, 'done'

            print >> sys.stderr, 'Loading the lexical units data ...',
            luData = cPickle.load(open(pickledLUPath))
            print >> sys.stderr, 'done'
            self['luIndex'] = luData

        self['luCache'] = {}
        self._generateFrameIndex()
        pass

    def _generatePickledFrames(self):
        """
        Initialises all the frames
        """

        framePath = FRAMENET_PATH+ '/' + FRAME_DIR_ENV + '/' + FRAME_FILE
        pickledFramePath = FRAMENET_PATH + '/' + PICKLED_FRAME_FILE

        print >> sys.stderr, 'Loading xml for frames ...',
        doc = xml.dom.minidom.parse(framePath)
        frameNodes = getNoneTextChildNodes(doc.childNodes[1])
        print >> sys.stderr, 'done',

        frames = {}
        print >> sys.stderr, 'parsing each frame ...',
        for fn in frameNodes:
            f = Frame()
            f.loadXMLNode(fn)
            frames[f['ID']] = f
        print >> sys.stderr, 'done',

        print >> sys.stderr, 'saving the frames ...',
        cPickle.dump(frames, open(pickledFramePath, 'w'), cPickle.HIGHEST_PROTOCOL)
        print >> sys.stderr, 'done'

        pass

    def _generatePickledFrameRelations(self):
        """
        Initialises all the frames
        """

        frRelationPath = FRAMENET_PATH + '/' + FRAME_DIR_ENV + '/' + FRAME_RELATION_FILE
        pickledFrameRelationPath = FRAMENET_PATH + '/' + PICKLED_FRAME_RELATIONS_FILE

        print >> sys.stderr, 'Loading xml for frame relations ...',
        frRelations = FrameRelation()
        frRelations.loadXML(frRelationPath)
        print >> sys.stderr, 'done',

        print >> sys.stderr, 'saving the frame relationships ...',
        cPickle.dump(frRelations, open(pickledFrameRelationPath, 'w'), cPickle.HIGHEST_PROTOCOL)
        print >> sys.stderr, 'done'

        pass

    def _generatePickledLexicalUnitsIndex(self):
        """
        Initialises all the frames
        """

        baseDir = FRAMENET_PATH + '/' + LU_DIR_ENV
        pickledLUPath = FRAMENET_PATH + '/' + PICKLED_LU_FILE

        lexicalUnitsIndexByName = {}
        for _f in os.listdir(baseDir):
            if _f.lower().startswith('lu') and _f.lower().endswith('.xml'):
                print >> sys.stderr, 'Loading:', _f, ": ",
                lu = LexicalUnit()
                lu.loadXML(baseDir + '/' + _f)
                print >> sys.stderr, lu['name'], '...',
                print >> sys.stderr, 'done'

                if lexicalUnitsIndexByName.has_key(lu['name']):
                    lexicalUnitsIndexByName[lu['name']][lu['ID']] = 1
                else:
                    lexicalUnitsIndexByName[lu['name']] = {lu['ID']:1}

        print >> sys.stderr, 'Saving the pickled lu files ...',
        cPickle.dump(lexicalUnitsIndexByName, open(pickledLUPath, 'w'), cPickle.HIGHEST_PROTOCOL)
        print >> sys.stderr, 'done'
        pass

    def _generateFrameIndex(self):
        """
        Generate an index from frame name to Frame objects
        """

        self['frameIndex'] = {}

        for _id in self['frames'].keys():
            f = self['frames'][_id]
            name = f['name']
            if self['frameIndex'].has_key(name):
                print >> sys.stderr, 'Error, multiple frame name:', name, 'found'
                sys.exit()

            self['frameIndex'][name] = f

        pass

    def initialize(self):
        """
        This function should be called only once before the whole thing can be used.
        """
        self._generatePickledFrames()
        self._generatePickledFrameRelations()
        self._generatePickledLexicalUnitsIndex()

        return True

    def lookupLexicalUnit(self, headWord, pos):
        """
        This function will look up a given word by its pos. The word must be already
        lemmatised and in lower case. The pos must also be in lower case and it can be
        one of the following:
        (1) v  -- for verb
        (2) n  -- for noun
        (3) a  -- for adjective
        (4) adv -- for adverb
        (5) prep -- for preposition
        (6) num -- for numbers
        (7) intj -- for interjections

        This function will return a dictionary of lexical units which match the (headWord, pos)
        pair. The keys to the dictionary will be the IDs of the lexical units, and the values of
        the dictionary will be LexicalUnit objects.
        """

        pickledLUPath = FRAMENET_PATH + '/' + LU_DIR_ENV

        w = headWord + '.' + pos

        if self['luCache'].has_key(w):
            return self['luCache'][w]

        if not self['luIndex'].has_key(w):
            return {}

        objects = {}
        for _id in self['luIndex'][w].keys():
            inputFile = pickledLUPath + '/lu' + str(_id) + '.xml'
            lu = LexicalUnit()
            lu.loadXML(inputFile)
            objects[lu['ID']] = lu

        self['luCache'][w] = objects
        return objects

    def getFrameGraph(self, frameText):
        import networkx as nx
        frame = self.lookupFrame(frameText)
        g = nx.DiGraph()
        frame_ids = frame.keys()
        for relations in self['frameRelations']['relation-types'].values():
            relationName = relations['name']
            superFrameName = relations['superFrameName']
            childFrameName = relations['subFrameName']
            rels = relations['frame-relations']
            print "Loading", relationName, superFrameName, childFrameName
            for frame_id in frame_ids:
                if True:
                #for k,rel in rels.items():
                    for k2, rel2 in rels.items():
                        if rel2['subFrameName'] == frame['name']:
                            print rel2

                        if rel2['ID'] == frame_id:
                            print "MATCH ON ID"
                            print rel2['subFrameName']
                        if int(rel2['subId']) == frame_id:
                            print rel2
                        elif int(rel2['supId']) == frame_id:
                            print rel2

    def lookupFrame(self, frame):
        """
        This function takes a string as input, the string should be the name
        of the Frame that you want to lookup, and its case sensitive. If not
        found, None will be returned.
        """

        if self['frameIndex'].has_key(frame):
            return self['frameIndex'][frame]

        return None

#----------------------------------------------------------------------------

max_stars_repo_name: dasmith/FrameNet-python

reason: import networkx