text:
import csv
import glob
import os
import numpy as np

import torch
import tqdm
from torch.nn.utils import clip_grad_norm_

def train_one_epoch(model, optimizer, train_loader, model_func, explained_model, explainer, attr_func, epoch_obj_cnt,
                    epoch_tp_obj_cnt, epoch_fp_obj_cnt, pred_score_file_name, pred_score_field_name, cur_epoch,
                    lr_scheduler, accumulated_iter, optim_cfg, rank, tbar, total_it_each_epoch, dataloader_iter,
                    cls_names, dataset_name, attr_loss, gt_infos, score_thresh, aggre_method, attr_sign, tb_log=None,
                    leave_pbar=False, box_selection="tp/fp", xc_goal="higher"):
    if total_it_each_epoch == len(train_loader):
        dataloader_iter = iter(train_loader)

    if rank == 0:
        pbar = tqdm.tqdm(total=total_it_each_epoch, leave=leave_pbar, desc='train', dynamic_ncols=True)

    # set to False in case we want verify something without altering the model
    enable_training = True
    constant_lr = False
    for cur_it in range(total_it_each_epoch):
        print("\nThe {}th batch\n".format(cur_it))
        try:
            batch = next(dataloader_iter)
        except StopIteration:
            dataloader_iter = iter(train_loader)
            batch = next(dataloader_iter)
            print('new iters')

        if not constant_lr:
            lr_scheduler.step(accumulated_iter)

        try:
            cur_lr = float(optimizer.lr)
        except:
            cur_lr = optimizer.param_groups[0]['lr']

        try:
            beta1, beta2 = float(optimizer.betas[0]), float(optimizer.betas[1])
        except:
            beta1, beta2 = optimizer.param_groups[0]['betas'][0], optimizer.param_groups[0]['betas'][1]

        if tb_log is not None:
            tb_log.add_scalar('meta_data/learning_rate', cur_lr, accumulated_iter)
        print("\nIn tools/train_utils/train_utils_new_loss.py, type(batch): {}".format(type(batch)))

        model.train()
        optimizer.zero_grad()

        # # Loading parameters from the training model to the explained model
        # explained_model.load_state_dict(model.state_dict(), strict=False)

        loss, tb_dict, disp_dict = model_func(model, batch)
        # print("\nin tools/train_utils/train_utils_new_loss.py, train_one_epoch, loss.shape {}".format(loss.shape)
        # print("loss data type is {}\n".format(type(loss)))
        if box_selection == "tp/fp":
            xc, far_attr, pap, fp_xc, fp_far_attr, fp_pap = attr_func(
                explained_model, explainer, batch, dataset_name, cls_names, cur_it=cur_it, gt_infos=gt_infos,
                score_thresh=score_thresh, object_cnt=epoch_obj_cnt, tp_object_cnt=epoch_tp_obj_cnt,
                fp_object_cnt=epoch_fp_obj_cnt, box_selection=box_selection, pred_score_file_name=pred_score_file_name,
                cur_epoch=cur_epoch, pred_score_field_name=pred_score_field_name, aggre_method=aggre_method,
                attr_sign=attr_sign, train_mode=True
            )
        else:
            xc, far_attr, pap = attr_func(
                explained_model, explainer, batch, dataset_name, cls_names, cur_it=cur_it, gt_infos=gt_infos,
                score_thresh=score_thresh, object_cnt=epoch_obj_cnt, tp_object_cnt=epoch_tp_obj_cnt,
                fp_object_cnt=epoch_fp_obj_cnt, box_selection=box_selection, pred_score_file_name=pred_score_file_name,
                cur_epoch=cur_epoch, pred_score_field_name=pred_score_field_name, aggre_method=aggre_method,
                attr_sign=attr_sign, train_mode=True
            )

        # need a scaling ratio to normalize the losses to the previous setting (top 3 per frame, batch_size = 2, taking
        # the average of frames in a batch)
        # TODO: handle the case when batch_box_cnt is zero
        skip_attr_loss = False
        # upper limit for the attr losses
        upper_limit = 0.3
        lambda_xc = 0.2
        lambda_pap = 0.0005

        if cur_epoch > 19 and cur_epoch < 60:
            lambda_pap = 0.00025
            lambda_xc = 0.1

        if xc is not None:
            valid_xc_flag = ~torch.isnan(xc)  # indicating where the xc values are not NaN
            zero_tensor = torch.zeros(xc.size(), dtype=xc.dtype).cuda()
            batch_box_cnt = torch.sum(valid_xc_flag).cuda()
            scaling_ratio = batch["batch_size"] * 3 / batch_box_cnt if batch_box_cnt != 0 else 0
            xc_val_raw = torch.sum(torch.where(valid_xc_flag, xc, zero_tensor)) / batch["batch_size"] * scaling_ratio
            xc_val_record = xc_val_raw.item()
            pap_val = torch.sum(torch.where(valid_xc_flag, pap, zero_tensor)) / batch["batch_size"] * scaling_ratio
            #far_attr_val = torch.sum(torch.where(valid_xc_flag, far_attr, zero_tensor)) / batch["batch_size"] * scaling_ratio
            if xc_goal == "higher":
                print('\nincreasing xc\n')
                three_val = 3 * torch.ones(xc_val_raw.size(), dtype=xc_val_raw.dtype).cuda()
                xc_val = torch.sub(three_val, xc_val_raw)  # xc_val = 3 - xc_val_raw, 3 since there are 3 boxes per frame
                # xc_val = 3 - xc_val_raw
            else:
                print('\nreducing xc\n')
                xc_val = xc_val_raw

        if (xc is None) or (batch_box_cnt == 0): # handles the case where we have no tp boxes
            skip_attr_loss = True
            xc_loss = 0
            pap_loss = 0
            far_attr_loss = 0
        elif xc is not None:
            # epsilon = 0.000001  # to avoid division by very small number
            # #
            # xc_val = xc_val if xc_val > epsilon else epsilon / xc_val.item() * xc_val
            # xc_val = torch.maximum(epsilon, xc_val)

            # lambda_tensor = lambda_ * torch.ones(xc_val.size(), dtype=xc_val.dtype).cuda()
            xc_loss_raw = lambda_xc * xc_val
            pap_loss_raw = lambda_pap * pap_val
            #far_attr_loss_raw = 0.01 * far_attr_val

            # to put an upper limit on the attr losses
            xc_loss = xc_loss_raw if xc_loss_raw < upper_limit else upper_limit / xc_loss_raw.item() * xc_loss_raw
            pap_loss = pap_loss_raw if pap_loss_raw < upper_limit else upper_limit / pap_loss_raw.item() * pap_loss_raw
            #far_attr_loss = far_attr_loss_raw if far_attr_loss_raw < upper_limit else upper_limit / far_attr_loss_raw.item() * far_attr_loss_raw
            print("\nxc_loss.requires_grad: {}".format(xc_loss.requires_grad))
            print("\npap_loss.requires_grad: {}".format(pap_loss.requires_grad))
            #print("\nfar_attr_loss.requires_grad: {}".format(far_attr_loss.requires_grad))
            print("\nloss.requires_grad: {}".format(loss.requires_grad))
        if box_selection == 'tp/fp':
            # assuming we always have at least 3 fp predictions in each frame
            print("\ncomputing fp_xc_loss\n")
            valid_fp_xc_flag = ~torch.isnan(fp_xc)
            fp_zero_tensor = torch.zeros(fp_xc.size(), dtype=fp_xc.dtype).cuda()
            fp_xc_loss_raw = lambda_xc * torch.sum(torch.where(valid_fp_xc_flag, fp_xc, fp_zero_tensor)) / batch["batch_size"]
            fp_pap_loss_raw = lambda_pap * torch.sum(torch.where(valid_fp_xc_flag, fp_pap, fp_zero_tensor)) / batch["batch_size"]
            fp_xc_loss = fp_xc_loss_raw if fp_xc_loss_raw < upper_limit else upper_limit / fp_xc_loss_raw.item() * fp_xc_loss_raw
            fp_pap_loss = fp_pap_loss_raw if fp_pap_loss_raw < upper_limit else upper_limit / fp_pap_loss_raw.item() * fp_pap_loss_raw
            # TODO: a proper far_attr loss for fp
        if attr_loss == 'xc' or attr_loss == 'XC':
            if not skip_attr_loss:  # only applicable in the tp and tp/fp mode
                loss += xc_loss
            if box_selection == 'tp/fp':
                print("\nadding fp_xc_loss\n")
                loss += fp_xc_loss
        elif attr_loss == 'pap' or attr_loss == 'PAP':
            if not skip_attr_loss:  # only applicable in the tp and tp/fp mode
                loss += pap_loss
            if box_selection == 'tp/fp':
                print("\nadding fp_xc_loss\n")
                loss += fp_pap_loss
        #elif attr_loss == 'far_attr' or attr_loss == 'FAR_ATTR':
            #if not skip_attr_loss:  # only applicable in the tp and tp/fp mode
                #loss += far_attr_loss
            # TODO: a proper far_attr loss for fp
        # print("loss: {}".format(loss))
        if enable_training:
            loss.backward()
            clip_grad_norm_(model.parameters(), optim_cfg.GRAD_NORM_CLIP)
            optimizer.step()

        accumulated_iter += 1
        disp_dict.update({'loss': loss.item(), 'lr': cur_lr})

        # log to console and tensorboard
        if rank == 0:
            pbar.update()
            pbar.set_postfix(dict(total_it=accumulated_iter))
            tbar.set_postfix(disp_dict)
            tbar.refresh()

            if tb_log is not None:
                tb_log.add_scalar('train/loss', loss, accumulated_iter)
                if xc is not None:
                    tb_log.add_scalar('train/xc', xc_val_record, accumulated_iter)
                    tb_log.add_scalar('train/xc_loss', xc_loss, accumulated_iter)
                    #tb_log.add_scalar('train/far_attr_loss', far_attr_loss, accumulated_iter)
                    tb_log.add_scalar('train/pap_loss', pap_loss, accumulated_iter)
                    if box_selection == 'tp/fp':
                        tb_log.add_scalar('train/fp_xc_loss', fp_xc_loss, accumulated_iter)
                        tb_log.add_scalar('train/fp_pap_loss', fp_pap_loss, accumulated_iter)
                    #tb_log.add_scalar('train/far_attr', far_attr_val, accumulated_iter)
                    tb_log.add_scalar('train/pap', pap_val, accumulated_iter)
                tb_log.add_scalar('meta_data/learning_rate', cur_lr, accumulated_iter)
                tb_log.add_scalar('meta_data/beta1', beta1, accumulated_iter)
                tb_log.add_scalar('meta_data/beta2', beta1, accumulated_iter)
                for key, val in tb_dict.items():
                    tb_log.add_scalar('train/' + key, val, accumulated_iter)
    if rank == 0:
        pbar.close()
    return accumulated_iter

def train_one_epoch_plain(model, optimizer, train_loader, model_func, explained_model, explainer, attr_func, epoch_obj_cnt,
                    epoch_tp_obj_cnt, epoch_fp_obj_cnt, pred_score_file_name, pred_score_field_name, cur_epoch,
                    lr_scheduler, accumulated_iter, optim_cfg, rank, tbar, total_it_each_epoch, dataloader_iter,
                    cls_names, dataset_name, attr_loss, gt_infos, score_thresh, aggre_method, attr_sign, tb_log=None,
                    leave_pbar=False, box_selection="tp/fp", xc_goal="higher"):
    if total_it_each_epoch == len(train_loader):
        dataloader_iter = iter(train_loader)

    if rank == 0:
        pbar = tqdm.tqdm(total=total_it_each_epoch, leave=leave_pbar, desc='train', dynamic_ncols=True)

    # set to False in case we want verify something without altering the model
    enable_training = True
    constant_lr = False
    for cur_it in range(total_it_each_epoch):
        print("\nThe {}th batch\n".format(cur_it))
        try:
            batch = next(dataloader_iter)
        except StopIteration:
            dataloader_iter = iter(train_loader)
            batch = next(dataloader_iter)
            print('new iters')

        if not constant_lr:
            lr_scheduler.step(accumulated_iter)

        try:
            cur_lr = float(optimizer.lr)
        except:
            cur_lr = optimizer.param_groups[0]['lr']

        try:
            beta1, beta2 = float(optimizer.betas[0]), float(optimizer.betas[1])
        except:
            beta1, beta2 = optimizer.param_groups[0]['betas'][0], optimizer.param_groups[0]['betas'][1]

        if tb_log is not None:
            tb_log.add_scalar('meta_data/learning_rate', cur_lr, accumulated_iter)
        print("\nIn tools/train_utils/train_utils_new_loss.py, type(batch): {}".format(type(batch)))

        model.train()
        optimizer.zero_grad()

        # # Loading parameters from the training model to the explained model
        # explained_model.load_state_dict(model.state_dict(), strict=False)

        loss, tb_dict, disp_dict = model_func(model, batch)

        # print("loss: {}".format(loss))
        if enable_training:
            loss.backward()
            clip_grad_norm_(model.parameters(), optim_cfg.GRAD_NORM_CLIP)
            optimizer.step()

        accumulated_iter += 1
        disp_dict.update({'loss': loss.item(), 'lr': cur_lr})

        # log to console and tensorboard
        if rank == 0:
            pbar.update()
            pbar.set_postfix(dict(total_it=accumulated_iter))
            tbar.set_postfix(disp_dict)
            tbar.refresh()

            if tb_log is not None:
                tb_log.add_scalar('train/loss', loss, accumulated_iter)
                tb_log.add_scalar('meta_data/learning_rate', cur_lr, accumulated_iter)
                tb_log.add_scalar('meta_data/beta1', beta1, accumulated_iter)
                tb_log.add_scalar('meta_data/beta2', beta1, accumulated_iter)
                for key, val in tb_dict.items():
                    tb_log.add_scalar('train/' + key, val, accumulated_iter)
    if rank == 0:
        pbar.close()
    return accumulated_iter

def train_model(model, optimizer, train_loader, logger, model_func, explained_model, explainer, attr_func,
                lr_scheduler, optim_cfg, start_epoch, total_epochs, start_iter, rank, tb_log, ckpt_save_dir, gt_infos,
                score_thresh, output_dir, aggre_method, attr_sign,
                train_sampler=None, lr_warmup_scheduler=None, ckpt_save_interval=1,
                max_ckpt_save_num=50, merge_all_iters_to_one_epoch=False, cls_names=['Car', 'Pedestrian', 'Cyclist'],
                dataset_name="KittiDataset", attr_loss="XC", box_selection="tp/fp", xc_goal="higher"):
    # setup logging files for training-related boxes info
    obj_cnt_file_name = output_dir / 'interested_obj_cnt.csv'
    pred_score_file_name = output_dir / 'interested_pred_scores.csv'
    obj_cnt_field_name = ['epoch', 'tp_cnt', 'fp_cnt', 'tp_car_cnt', 'tp_pede_cnt', 'tp_cyc_cnt', 'fp_car_cnt',
                          'fp_pede_cnt', 'fp_cyc_cnt', 'car_cnt', 'pede_cnt', 'cyc_cnt']
    pred_score_field_name = ['epoch', 'batch', 'tp/fp', 'pred_label', 'pred_score', 'xc', 'far_attr', 'pap']
    if box_selection != "tp/fp" and box_selection != "tp":
        obj_cnt_field_name = ['epoch', 'car_cnt', 'pede_cnt', 'cyc_cnt']
        pred_score_field_name = ['epoch', 'batch', 'pred_label', 'pred_score', 'xc', 'far_attr', 'pap']
    with open(obj_cnt_file_name, 'w', newline='') as csvfile:
        writer = csv.DictWriter(csvfile, delimiter=',', fieldnames=obj_cnt_field_name)
        writer.writeheader()
    with open(pred_score_file_name, 'w', newline='') as csvfile:
        writer = csv.DictWriter(csvfile, delimiter=',', fieldnames=pred_score_field_name)
        writer.writeheader()
    accumulated_iter = start_iter
    if not (attr_loss == "XC" or attr_loss == "xc" or attr_loss == "PAP" or attr_loss == "pap" or
            attr_loss == "far_attr" or attr_loss == "FAR_ATTR" or attr_loss == "None" or
            attr_loss == "none_but_calc" or attr_loss == "None_but_calc"):
        raise NotImplementedError
    with tqdm.trange(start_epoch, total_epochs, desc='epochs', dynamic_ncols=True, leave=(rank == 0)) as tbar:
        total_it_each_epoch = len(train_loader)
        if merge_all_iters_to_one_epoch:
            assert hasattr(train_loader.dataset, 'merge_all_iters_to_one_epoch')
            train_loader.dataset.merge_all_iters_to_one_epoch(merge=True, epochs=total_epochs)
            total_it_each_epoch = len(train_loader) // max(total_epochs, 1)

        dataloader_iter = iter(train_loader)
        epoch_obj_cnt = {}
        epoch_tp_obj_cnt = {}
        epoch_fp_obj_cnt = {}
        for i in range(len(cls_names)):
            epoch_obj_cnt[i] = 0
            epoch_tp_obj_cnt[i] = 0
            epoch_fp_obj_cnt[i] = 0
        for cur_epoch in tbar:
            if train_sampler is not None:
                train_sampler.set_epoch(cur_epoch)
            # train one epoch
            if lr_warmup_scheduler is not None and cur_epoch < optim_cfg.WARMUP_EPOCH:
                cur_scheduler = lr_warmup_scheduler
            else:
                cur_scheduler = lr_scheduler

            if attr_loss == "None" or attr_loss == "none":
                # when no attr_loss is applied
                accumulated_iter = train_one_epoch_plain(
                    model, optimizer, train_loader, model_func, explained_model, explainer, attr_func, epoch_obj_cnt,
                    epoch_tp_obj_cnt, epoch_fp_obj_cnt,
                    cur_epoch=cur_epoch,
                    lr_scheduler=cur_scheduler,
                    accumulated_iter=accumulated_iter, optim_cfg=optim_cfg,
                    rank=rank, tbar=tbar, tb_log=tb_log,
                    leave_pbar=(cur_epoch + 1 == total_epochs),
                    total_it_each_epoch=total_it_each_epoch,
                    dataloader_iter=dataloader_iter,
                    cls_names=cls_names, dataset_name=dataset_name, attr_loss=attr_loss, gt_infos=gt_infos,
                    score_thresh=score_thresh, box_selection=box_selection, pred_score_file_name=pred_score_file_name,
                    pred_score_field_name=pred_score_field_name, aggre_method=aggre_method, attr_sign=attr_sign
                )
            else:
                accumulated_iter = train_one_epoch(
                    model, optimizer, train_loader, model_func, explained_model, explainer, attr_func, epoch_obj_cnt,
                    epoch_tp_obj_cnt, epoch_fp_obj_cnt,
                    cur_epoch=cur_epoch,
                    lr_scheduler=cur_scheduler,
                    accumulated_iter=accumulated_iter, optim_cfg=optim_cfg,
                    rank=rank, tbar=tbar, tb_log=tb_log,
                    leave_pbar=(cur_epoch + 1 == total_epochs),
                    total_it_each_epoch=total_it_each_epoch,
                    dataloader_iter=dataloader_iter,
                    cls_names=cls_names, dataset_name=dataset_name, attr_loss=attr_loss, gt_infos=gt_infos,
                    score_thresh=score_thresh, box_selection=box_selection, pred_score_file_name=pred_score_file_name,
                    pred_score_field_name=pred_score_field_name, aggre_method=aggre_method, attr_sign=attr_sign,
                    xc_goal=xc_goal
                )
            # log object count
            with open(obj_cnt_file_name, 'a', newline='') as csvfile:
                writer = csv.DictWriter(csvfile, delimiter=',', fieldnames=obj_cnt_field_name)
                data_dict = {}
                data_dict["epoch"] = cur_epoch
                # TODO: fill the data dict
                if box_selection == "tp/fp" or box_selection == "tp":
                    data_dict["tp_car_cnt"], data_dict["fp_car_cnt"] = epoch_tp_obj_cnt[0], epoch_fp_obj_cnt[0]
                    data_dict["tp_pede_cnt"], data_dict["fp_pede_cnt"] = epoch_tp_obj_cnt[1], epoch_fp_obj_cnt[1]
                    data_dict["tp_cyc_cnt"], data_dict["fp_cyc_cnt"] = epoch_tp_obj_cnt[2], epoch_fp_obj_cnt[2]
                    data_dict["tp_cnt"], data_dict["fp_cnt"] = sum(epoch_tp_obj_cnt.values()), sum(
                        epoch_fp_obj_cnt.values())
                data_dict["car_cnt"], data_dict["pede_cnt"], data_dict["cyc_cnt"] = epoch_obj_cnt[0], epoch_obj_cnt[1], \
                                                                                epoch_obj_cnt[2]
                writer.writerow(data_dict)
            logger.info("{}:{} {}:{} {}:{} after training for {} epochs".format(
                cls_names[0], epoch_obj_cnt[0], cls_names[1], epoch_obj_cnt[1], cls_names[2], epoch_obj_cnt[2],
                cur_epoch))
            if box_selection == "tp/fp":
                logger.info("tp {}:{} {}:{} {}:{} after training for {} epochs".format(
                    cls_names[0], epoch_tp_obj_cnt[0], cls_names[1], epoch_tp_obj_cnt[1], cls_names[2], epoch_tp_obj_cnt[2],
                    cur_epoch))
                logger.info("fp {}:{} {}:{} {}:{} after training for {} epochs".format(
                    cls_names[0], epoch_fp_obj_cnt[0], cls_names[1], epoch_fp_obj_cnt[1], cls_names[2], epoch_fp_obj_cnt[2],
                    cur_epoch))
            # save trained model
            trained_epoch = cur_epoch + 1
            if trained_epoch % ckpt_save_interval == 0 and rank == 0:

                ckpt_list = glob.glob(str(ckpt_save_dir / 'checkpoint_epoch_*.pth'))
                ckpt_list.sort(key=os.path.getmtime)

                if ckpt_list.__len__() >= max_ckpt_save_num:
                    for cur_file_idx in range(0, len(ckpt_list) - max_ckpt_save_num + 1):
                        os.remove(ckpt_list[cur_file_idx])

                ckpt_name = ckpt_save_dir / ('checkpoint_epoch_%d' % trained_epoch)
                save_checkpoint(
                    checkpoint_state(model, optimizer, trained_epoch, accumulated_iter), filename=ckpt_name,
                )

def model_state_to_cpu(model_state):
    model_state_cpu = type(model_state)()  # ordered dict
    for key, val in model_state.items():
        model_state_cpu[key] = val.cpu()
    return model_state_cpu

def checkpoint_state(model=None, optimizer=None, epoch=None, it=None):
    optim_state = optimizer.state_dict() if optimizer is not None else None
    if model is not None:
        if isinstance(model, torch.nn.parallel.DistributedDataParallel):
            model_state = model_state_to_cpu(model.module.state_dict())
        else:
            model_state = model.state_dict()
    else:
        model_state = None

    try:
        import pcdet
        version = 'pcdet+' + pcdet.__version__
    except:
        version = 'none'

    return {'epoch': epoch, 'it': it, 'model_state': model_state, 'optimizer_state': optim_state, 'version': version}

def save_checkpoint(state, filename='checkpoint'):
    if False and 'optimizer_state' in state:
        optimizer_state = state['optimizer_state']
        state.pop('optimizer_state', None)
        optimizer_filename = '{}_optim.pth'.format(filename)
        torch.save({'optimizer_state': optimizer_state}, optimizer_filename)

    filename = '{}.pth'.format(filename)
    torch.save(state, filename)


text:
"""
=================================
Travelling Salesman Problem (TSP)
=================================

Implementation of approximate algorithms
for solving and approximating the TSP problem.

Categories of algorithms which are implemented:

- Christofides (provides a 3/2-approximation of TSP)
- Greedy
- Simulated Annealing (SA)
- Threshold Accepting (TA)
- Asadpour Asymmetric Traveling Salesman Algorithm

The Travelling Salesman Problem tries to find, given the weight
(distance) between all points where a salesman has to visit, the
route so that:

- The total distance (cost) which the salesman travels is minimized.
- The salesman returns to the starting point.
- Note that for a complete graph, the salesman visits each point once.

The function `travelling_salesman_problem` allows for incomplete
graphs by finding all-pairs shortest paths, effectively converting
the problem to a complete graph problem. It calls one of the
approximate methods on that problem and then converts the result
back to the original graph using the previously found shortest paths.

TSP is an NP-hard problem in combinatorial optimization,
important in operations research and theoretical computer science.

http://en.wikipedia.org/wiki/Travelling_salesman_problem
"""
import math

import networkx as nx

from networkx.utils import py_random_state, not_implemented_for, pairwise

__all__ = [
    "traveling_salesman_problem",
    "christofides",
    "asadpour_atsp",
    "greedy_tsp",
    "simulated_annealing_tsp",
    "threshold_accepting_tsp",
]

def swap_two_nodes(soln, seed):
    """Swap two nodes in `soln` to give a neighbor solution.

    Parameters
    ----------
    soln : list of nodes
        Current cycle of nodes

    seed : integer, random_state, or None (default)
        Indicator of random number generation state.
        See :ref:`Randomness<randomness>`.

    Returns
    -------
    list
        The solution after move is applied. (A neighbor solution.)

    Notes
    -----
        This function assumes that the incoming list `soln` is a cycle
        (that the first and last element are the same) and also that
        we don't want any move to change the first node in the list
        (and thus not the last node either).

        The input list is changed as well as returned. Make a copy if needed.

    See Also
    --------
        move_one_node
    """
    a, b = seed.sample(range(1, len(soln) - 1), k=2)
    soln[a], soln[b] = soln[b], soln[a]
    return soln

def move_one_node(soln, seed):
    """Move one node to another position to give a neighbor solution.

    The node to move and the position to move to are chosen randomly.
    The first and last nodes are left untouched as soln must be a cycle
    starting at that node.

    Parameters
    ----------
    soln : list of nodes
        Current cycle of nodes

    seed : integer, random_state, or None (default)
        Indicator of random number generation state.
        See :ref:`Randomness<randomness>`.

    Returns
    -------
    list
        The solution after move is applied. (A neighbor solution.)

    Notes
    -----
        This function assumes that the incoming list `soln` is a cycle
        (that the first and last element are the same) and also that
        we don't want any move to change the first node in the list
        (and thus not the last node either).

        The input list is changed as well as returned. Make a copy if needed.

    See Also
    --------
        swap_two_nodes
    """
    a, b = seed.sample(range(1, len(soln) - 1), k=2)
    soln.insert(b, soln.pop(a))
    return soln

@not_implemented_for("directed")
def christofides(G, weight="weight", tree=None):
    """Approximate a solution of the traveling salesman problem

    Compute a 3/2-approximation of the traveling salesman problem
    in a complete undirected graph using Christofides [1]_ algorithm.

    Parameters
    ----------
    G : Graph
        `G` should be a complete weighted undirected graph.
        The distance between all pairs of nodes should be included.

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    tree : NetworkX graph or None (default: None)
        A minimum spanning tree of G. Or, if None, the minimum spanning
        tree is computed using :func:`networkx.minimum_spanning_tree`

    Returns
    -------
    list
        List of nodes in `G` along a cycle with a 3/2-approximation of
        the minimal Hamiltonian cycle.

    References
    ----------
    .. [1] Christofides, Nicos. "Worst-case analysis of a new heuristic for
       the travelling salesman problem." No. RR-388. Carnegie-Mellon Univ
       Pittsburgh Pa Management Sciences Research Group, 1976.
    """
    # Remove selfloops if necessary
    loop_nodes = nx.nodes_with_selfloops(G)
    try:
        node = next(loop_nodes)
    except StopIteration:
        pass
    else:
        G = G.copy()
        G.remove_edge(node, node)
        G.remove_edges_from((n, n) for n in loop_nodes)
    # Check that G is a complete graph
    N = len(G) - 1
    # This check ignores selfloops which is what we want here.
    if any(len(nbrdict) != N for n, nbrdict in G.adj.items()):
        raise nx.NetworkXError("G must be a complete graph.")

    if tree is None:
        tree = nx.minimum_spanning_tree(G, weight=weight)
    L = G.copy()
    L.remove_nodes_from([v for v, degree in tree.degree if not (degree % 2)])
    MG = nx.MultiGraph()
    MG.add_edges_from(tree.edges)
    edges = nx.min_weight_matching(L, maxcardinality=True, weight=weight)
    MG.add_edges_from(edges)
    return _shortcutting(nx.eulerian_circuit(MG))

def _shortcutting(circuit):
    """Remove duplicate nodes in the path"""
    nodes = []
    for u, v in circuit:
        if v in nodes:
            continue
        if not nodes:
            nodes.append(u)
        nodes.append(v)
    nodes.append(nodes[0])
    return nodes

def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, method=None):
    """Find the shortest path in `G` connecting specified nodes

    This function allows approximate solution to the traveling salesman
    problem on networks that are not complete graphs and/or where the
    salesman does not need to visit all nodes.

    This function proceeds in two steps. First, it creates a complete
    graph using the all-pairs shortest_paths between nodes in `nodes`.
    Edge weights in the new graph are the lengths of the paths
    between each pair of nodes in the original graph.
    Second, an algorithm (default: `christofides` for undirected and
    `asadpour_atsp` for directed) is used to approximate the minimal Hamiltonian
    cycle on this new graph. The available algorithms are:

     - christofides
     - greedy_tsp
     - simulated_annealing_tsp
     - threshold_accepting_tsp
     - asadpour_atsp

    Once the Hamiltonian Cycle is found, this function post-processes to
    accommodate the structure of the original graph. If `cycle` is ``False``,
    the biggest weight edge is removed to make a Hamiltonian path.
    Then each edge on the new complete graph used for that analysis is
    replaced by the shortest_path between those nodes on the original graph.

    Parameters
    ----------
    G : NetworkX graph
        A possibly weighted graph

    nodes : collection of nodes (default=G.nodes)
        collection (list, set, etc.) of nodes to visit

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    cycle : bool (default: True)
        Indicates whether a cycle should be returned, or a path.
        Note: the cycle is the approximate minimal cycle.
        The path simply removes the biggest edge in that cycle.

    method : function (default: None)
        A function that returns a cycle on all nodes and approximates
        the solution to the traveling salesman problem on a complete
        graph. The returned cycle is then used to find a corresponding
        solution on `G`. `method` should be callable; take inputs
        `G`, and `weight`; and return a list of nodes along the cycle.

        Provided options include :func:`christofides`, :func:`greedy_tsp`,
        :func:`simulated_annealing_tsp` and :func:`threshold_accepting_tsp`.

        If `method is None`: use :func:`christofides` for undirected `G` and
        :func:`threshold_accepting_tsp` for directed `G`.

        To specify parameters for these provided functions, construct lambda
        functions that state the specific value. `method` must have 2 inputs.
        (See examples).

    Returns
    -------
    list
        List of nodes in `G` along a path with an approximation of the minimal
        path through `nodes`.

    Raises
    ------
    NetworkXError
        If `G` is a directed graph it has to be strongly connected or the
        complete version cannot be generated.

    Examples
    --------
    >>> tsp = nx.approximation.traveling_salesman_problem
    >>> G = nx.cycle_graph(9)
    >>> G[4][5]["weight"] = 5  # all other weights are 1
    >>> tsp(G, nodes=[3, 6])
    [3, 2, 1, 0, 8, 7, 6, 7, 8, 0, 1, 2, 3]
    >>> path = tsp(G, cycle=False)
    >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
    True

    Build (curry) your own function to provide parameter values to the methods.

    >>> SA_tsp = nx.approximation.simulated_annealing_tsp
    >>> method = lambda G, wt: SA_tsp(G, "greedy", weight=wt, temp=500)
    >>> path = tsp(G, cycle=False, method=method)
    >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
    True

    """
    if method is None:
        if G.is_directed():
            method = asadpour_atsp
        else:
            method = christofides
    if nodes is None:
        nodes = list(G.nodes)

    dist = {}
    path = {}
    for n, (d, p) in nx.all_pairs_dijkstra(G, weight=weight):
        dist[n] = d
        path[n] = p

    if G.is_directed():
        # If the graph is not strongly connected, raise an exception
        if not nx.is_strongly_connected(G):
            raise nx.NetworkXError("G is not strongly connected")
        GG = nx.DiGraph()
    else:
        GG = nx.Graph()
    for u in nodes:
        for v in nodes:
            if u == v:
                continue
            GG.add_edge(u, v, weight=dist[u][v])
    best_GG = method(GG, weight)

    if not cycle:
        # find and remove the biggest edge
        (u, v) = max(pairwise(best_GG), key=lambda x: dist[x[0]][x[1]])
        pos = best_GG.index(u) + 1
        while best_GG[pos] != v:
            pos = best_GG[pos:].index(u) + 1
        best_GG = best_GG[pos:-1] + best_GG[:pos]

    best_path = []
    for u, v in pairwise(best_GG):
        best_path.extend(path[u][v][:-1])
    best_path.append(v)
    return best_path

@not_implemented_for("undirected")
def asadpour_atsp(G, weight="weight", seed=None, source=None):
    """
    Returns an approximate solution to the traveling salesman problem.

    This approximate solution is one of the best known approximations for the
    asymmetric traveling salesman problem developed by Asadpour et al,
    [1]_. The algorithm first solves the Held-Karp relaxation to find a lower
    bound for the weight of the cycle. Next, it constructs an exponential
    distribution of undirected spanning trees where the probability of an
    edge being in the tree corresponds to the weight of that edge using a
    maximum entropy rounding scheme. Next we sample that distribution
    $2 \\lceil \\ln n \\rceil$ times and save the minimum sampled tree once the
    direction of the arcs is added back to the edges. Finally, we augment
    then short circuit that graph to find the approximate tour for the
    salesman.

    Parameters
    ----------
    G : nx.DiGraph
        The graph should be a complete weighted directed graph. The
        distance between all paris of nodes should be included and the triangle
        inequality should hold. That is, the direct edge between any two nodes
        should be the path of least cost.

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    seed : integer, random_state, or None (default)
        Indicator of random number generation state.
        See :ref:`Randomness<randomness>`.

    source : node label (default=`None`)
        If given, return the cycle starting and ending at the given node.

    Returns
    -------
    cycle : list of nodes
        Returns the cycle (list of nodes) that a salesman can follow to minimize
        the total weight of the trip.

    Raises
    ------
    NetworkXError
        If `G` is not complete or has less than two nodes, the algorithm raises
        an exception.

    NetworkXError
        If 'source` is not `None` and is not a node in `G`, the algorithm raises
        an exception.

    NetworkXNotImplemented
        If `G` is an undirected graph.

    References
    ----------
    .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
       An o(log n/log log n)-approximation algorithm for the asymmetric
       traveling salesman problem, Operations research, 65 (2017),
       pp. 1043–1061

    Examples
    --------
    >>> import networkx as nx
    >>> import networkx.algorithms.approximation as approx
    >>> G = nx.complete_graph(3, create_using=nx.DiGraph)
    >>> nx.set_edge_attributes(G, {(0, 1): 2, (1, 2): 2, (2, 0): 2, (0, 2): 1, (2, 1): 1, (1, 0): 1}, "weight")
    >>> tour = approx.asadpour_atsp(G,source=0)
    >>> tour
    [0, 2, 1, 0]
    """
    from math import log as ln
    from math import exp
    from math import ceil

    # Check that G is a complete graph
    N = len(G) - 1
    if N < 2:
        raise nx.NetworkXError("G must have at least two nodes")
    # This check ignores selfloops which is what we want here.
    if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
        raise nx.NetworkXError("G is not a complete DiGraph")
    # Check that the source vertex, if given, is in the graph
    if source is not None and source not in G.nodes:
        raise nx.NetworkXError("Given source node not in G.")

    opt_hk, z_star = held_karp_ascent(G, weight)

    # Test to see if the ascent method found an integer solution or a fractional
    # solution. If it is integral then z_star is a nx.Graph, otherwise it is
    # a dict
    if not isinstance(z_star, dict):
        # Here we are using the shortcutting method to go from the list of edges
        # returned from eularian_circuit to a list of nodes
        return _shortcutting(nx.eulerian_circuit(z_star, source=source))

    # Create the undirected support of z_star
    z_support = nx.MultiGraph()
    for u, v in z_star:
        if (u, v) not in z_support.edges:
            edge_weight = min(G[u][v][weight], G[v][u][weight])
            z_support.add_edge(u, v, **{weight: edge_weight})

    # Create the exponential distribution of spanning trees
    gamma = spanning_tree_distribution(z_support, z_star)

    # Write the lambda values to the edges of z_support
    z_support = nx.Graph(z_support)
    lambda_dict = {(u, v): exp(gamma[(u, v)]) for u, v in z_support.edges()}
    nx.set_edge_attributes(z_support, lambda_dict, "lambda_key")
    del gamma, lambda_dict

    # Sample 2 * ceil( ln(n) ) spanning trees and record the minimum one
    minimum_sampled_tree = None
    minimum_sampled_tree_weight = math.inf
    for _ in range(2 * ceil(ln(G.number_of_nodes()))):
        sampled_tree = sample_spanning_tree(z_support, "lambda_key", seed)
        sampled_tree_weight = sampled_tree.size(weight)
        if sampled_tree_weight < minimum_sampled_tree_weight:
            minimum_sampled_tree = sampled_tree.copy()
            minimum_sampled_tree_weight = sampled_tree_weight

    # Orient the edges in that tree to keep the cost of the tree the same.
    t_star = nx.MultiDiGraph()
    for u, v, d in minimum_sampled_tree.edges(data=weight):
        if d == G[u][v][weight]:
            t_star.add_edge(u, v, **{weight: d})
        else:
            t_star.add_edge(v, u, **{weight: d})

    # Find the node demands needed to neutralize the flow of t_star in G
    node_demands = {n: t_star.out_degree(n) - t_star.in_degree(n) for n in t_star}
    nx.set_node_attributes(G, node_demands, "demand")

    # Find the min_cost_flow
    flow_dict = nx.min_cost_flow(G, "demand")

    # Build the flow into t_star
    for source, values in flow_dict.items():
        for target in values:
            if (source, target) not in t_star.edges and values[target] > 0:
                # IF values[target] > 0 we have to add that many edges
                for _ in range(values[target]):
                    t_star.add_edge(source, target)

    # Return the shortcut eulerian circuit
    circuit = nx.eulerian_circuit(t_star, source=source)
    return _shortcutting(circuit)

def held_karp_ascent(G, weight="weight"):
    """
    Minimizes the Held-Karp relaxation of the TSP for `G`

    Solves the Held-Karp relaxation of the input complete digraph and scales
    the output solution for use in the Asadpour [1]_ ASTP algorithm.

    The Held-Karp relaxation defines the lower bound for solutions to the
    ATSP, although it does return a fractional solution. This is used in the
    Asadpour algorithm as an initial solution which is later rounded to a
    integral tree within the spanning tree polytopes. This function solves
    the relaxation with the branch and bound method in [2]_.

    Parameters
    ----------
    G : nx.DiGraph
        The graph should be a complete weighted directed graph.
        The distance between all paris of nodes should be included.

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    Returns
    -------
    OPT : float
        The cost for the optimal solution to the Held-Karp relaxation
    z : dict or nx.Graph
        A symmetrized and scaled version of the optimal solution to the
        Held-Karp relaxation for use in the Asadpour algorithm.

        If an integral solution is found, then that is an optimal solution for
        the ATSP problem and that is returned instead.

    References
    ----------
    .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
       An o(log n/log log n)-approximation algorithm for the asymmetric
       traveling salesman problem, Operations research, 65 (2017),
       pp. 1043–1061

    .. [2] M. Held, R. M. Karp, The traveling-salesman problem and minimum
           spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
           pp.1138-1162
    """
    import numpy as np
    import scipy.optimize as optimize

    def k_pi():
        """
        Find the set of minimum 1-Arborescences for G at point pi.

        Returns
        -------
        Set
            The set of minimum 1-Arborescences
        """
        # Create a copy of G without vertex 1.
        G_1 = G.copy()
        minimum_1_arborescences = set()
        minimum_1_arborescence_weight = math.inf

        # node is node '1' in the Held and Karp paper
        n = next(G.__iter__())
        G_1.remove_node(n)

        # Iterate over the spanning arborescences of the graph until we know
        # that we have found the minimum 1-arborescences. My proposed strategy
        # is to find the most extensive root to connect to from 'node 1' and
        # the least expensive one. We then iterate over arborescences until
        # the cost of the basic arborescence is the cost of the minimum one
        # plus the difference between the most and least expensive roots,
        # that way the cost of connecting 'node 1' will by definition not by
        # minimum
        min_root = {"node": None, weight: math.inf}
        max_root = {"node": None, weight: -math.inf}
        for u, v, d in G.edges(n, data=True):
            if d[weight] < min_root[weight]:
                min_root = {"node": v, weight: d[weight]}
            if d[weight] > max_root[weight]:
                max_root = {"node": v, weight: d[weight]}

        min_in_edge = min(G.in_edges(n, data=True), key=lambda x: x[2][weight])
        min_root[weight] = min_root[weight] + min_in_edge[2][weight]
        max_root[weight] = max_root[weight] + min_in_edge[2][weight]

        min_arb_weight = math.inf
        for arb in nx.ArborescenceIterator(G_1):
            arb_weight = arb.size(weight)
            if min_arb_weight == math.inf:
                min_arb_weight = arb_weight
            elif arb_weight > min_arb_weight + max_root[weight] - min_root[weight]:
                break
            # We have to pick the root node of the arborescence for the out
            # edge of the first vertex as that is the only node without an
            # edge directed into it.
            for N, deg in arb.in_degree:
                if deg == 0:
                    # root found
                    arb.add_edge(n, N, **{weight: G[n][N][weight]})
                    arb_weight += G[n][N][weight]
                    break

            # We can pick the minimum weight in-edge for the vertex with
            # a cycle. If there are multiple edges with the same, minimum
            # weight, We need to add all of them.
            #
            # Delete the edge (N, v) so that we cannot pick it.
            edge_data = G[N][n]
            G.remove_edge(N, n)
            min_weight = min(G.in_edges(n, data=weight), key=lambda x: x[2])[2]
            min_edges = [
                (u, v, d) for u, v, d in G.in_edges(n, data=weight) if d == min_weight
            ]
            for u, v, d in min_edges:
                new_arb = arb.copy()
                new_arb.add_edge(u, v, **{weight: d})
                new_arb_weight = arb_weight + d
                # Check to see the weight of the arborescence, if it is a
                # new minimum, clear all of the old potential minimum
                # 1-arborescences and add this is the only one. If its
                # weight is above the known minimum, do not add it.
                if new_arb_weight < minimum_1_arborescence_weight:
                    minimum_1_arborescences.clear()
                    minimum_1_arborescence_weight = new_arb_weight
                # We have a 1-arborescence, add it to the set
                if new_arb_weight == minimum_1_arborescence_weight:
                    minimum_1_arborescences.add(new_arb)
            G.add_edge(N, n, **edge_data)

        return minimum_1_arborescences

    def direction_of_ascent():
        """
        Find the direction of ascent at point pi.

        See [1]_ for more information.

        Returns
        -------
        dict
            A mapping from the nodes of the graph which represents the direction
            of ascent.

        References
        ----------
        .. [1] M. Held, R. M. Karp, The traveling-salesman problem and minimum
           spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
           pp.1138-1162
        """
        # 1. Set d equal to the zero n-vector.
        d = {}
        for n in G:
            d[n] = 0
        del n
        # 2. Find a 1-Aborescence T^k such that k is in K(pi, d).
        minimum_1_arborescences = k_pi()
        while True:
            # Reduce K(pi) to K(pi, d)
            # Find the arborescence in K(pi) which increases the lest in
            # direction d
            min_k_d_weight = math.inf
            min_k_d = None
            for arborescence in minimum_1_arborescences:
                weighted_cost = 0
                for n, deg in arborescence.degree:
                    weighted_cost += d[n] * (deg - 2)
                if weighted_cost < min_k_d_weight:
                    min_k_d_weight = weighted_cost
                    min_k_d = arborescence

            # 3. If sum of d_i * v_{i, k} is greater than zero, terminate
            if min_k_d_weight > 0:
                return d, min_k_d
            # 4. d_i = d_i + v_{i, k}
            for n, deg in min_k_d.degree:
                d[n] += deg - 2
            # Check that we do not need to terminate because the direction
            # of ascent does not exist. This is done with linear
            # programming.
            c = np.full(len(minimum_1_arborescences), -1, dtype=int)
            a_eq = np.empty((len(G) + 1, len(minimum_1_arborescences)), dtype=int)
            b_eq = np.zeros(len(G) + 1, dtype=int)
            b_eq[len(G)] = 1
            arb_count = 0
            for arborescence in minimum_1_arborescences:
                n_count = len(G) - 1
                for n, deg in arborescence.degree:
                    a_eq[n_count][arb_count] = deg - 2
                    n_count -= 1
                a_eq[len(G)][arb_count] = 1
                arb_count += 1
            program_result = optimize.linprog(c, A_eq=a_eq, b_eq=b_eq)
            bool_result = program_result.x >= 0
            if program_result.status == 0 and np.sum(bool_result) == len(
                minimum_1_arborescences
            ):
                # There is no direction of ascent
                return None, minimum_1_arborescences

            # 5. GO TO 2

    def find_epsilon(k, d):
        """
        Given the direction of ascent at pi, find the maximum distance we can go
        in that direction.

        Parameters
        ----------
        k_xy : set
            The set of 1-arborescences which have the minimum rate of increase
            in the direction of ascent

        d : dict
            The direction of ascent

        Returns
        -------
        float
            The distance we can travel in direction `d`
        """
        min_epsilon = math.inf
        for e_u, e_v, e_w in G.edges(data=weight):
            if (e_u, e_v) in k.edges:
                continue
            # Now, I have found a condition which MUST be true for the edges to
            # be a valid substitute. The edge in the graph which is the
            # substitute is the one with the same terminated end. This can be
            # checked rather simply.
            #
            # Find the edge within k which is the substitute. Because k is a
            # 1-arborescence, we know that they is only one such edges
            # leading into every vertex.
            if len(k.in_edges(e_v, data=weight)) > 1:
                raise Exception
            sub_u, sub_v, sub_w = next(k.in_edges(e_v, data=weight).__iter__())
            k.add_edge(e_u, e_v, **{weight: e_w})
            k.remove_edge(sub_u, sub_v)
            if (
                max(d for n, d in k.in_degree()) <= 1
                and len(G) == k.number_of_edges()
                and nx.is_weakly_connected(k)
            ):
                # Ascent method calculation
                if d[sub_u] == d[e_u] or sub_w == e_w:
                    # Revert to the original graph
                    k.remove_edge(e_u, e_v)
                    k.add_edge(sub_u, sub_v, **{weight: sub_w})
                    continue
                epsilon = (sub_w - e_w) / (d[e_u] - d[sub_u])
                if 0 < epsilon < min_epsilon:
                    min_epsilon = epsilon
            # Revert to the original graph
            k.remove_edge(e_u, e_v)
            k.add_edge(sub_u, sub_v, **{weight: sub_w})

        return min_epsilon

    # I have to know that the elements in pi correspond to the correct elements
    # in the direction of ascent, even if the node labels are not integers.
    # Thus, I will use dictionaries to made that mapping.
    pi_dict = {}
    for n in G:
        pi_dict[n] = 0
    del n
    original_edge_weights = {}
    for u, v, d in G.edges(data=True):
        original_edge_weights[(u, v)] = d[weight]
    dir_ascent, k_d = direction_of_ascent()
    while dir_ascent is not None:
        max_distance = find_epsilon(k_d, dir_ascent)
        for n, v in dir_ascent.items():
            pi_dict[n] += max_distance * v
        for u, v, d in G.edges(data=True):
            d[weight] = original_edge_weights[(u, v)] + pi_dict[u]
        dir_ascent, k_d = direction_of_ascent()
    # k_d is no longer an individual 1-arborescence but rather a set of
    # minimal 1-arborescences at the maximum point of the polytope and should
    # be reflected as such
    k_max = k_d

    # Search for a cycle within k_max. If a cycle exists, return it as the
    # solution
    for k in k_max:
        if len([n for n in k if k.degree(n) == 2]) == G.order():
            # Tour found
            return k.size(weight), k

    # Write the original edge weights back to G and every member of k_max at
    # the maximum point. Also average the number of times that edge appears in
    # the set of minimal 1-arborescences.
    x_star = {}
    size_k_max = len(k_max)
    for u, v, d in G.edges(data=True):
        edge_count = 0
        d[weight] = original_edge_weights[(u, v)]
        for k in k_max:
            if (u, v) in k.edges():
                edge_count += 1
                k[u][v][weight] = original_edge_weights[(u, v)]
        x_star[(u, v)] = edge_count / size_k_max
    # Now symmetrize the edges in x_star and scale them according to (5) in
    # reference [1]
    z_star = {}
    scale_factor = (G.order() - 1) / G.order()
    for u, v in x_star.keys():
        frequency = x_star[(u, v)] + x_star[(v, u)]
        if frequency > 0:
            z_star[(u, v)] = scale_factor * frequency
    del x_star
    # Return the optimal weight and the z dict
    return next(k_max.__iter__()).size(weight), z_star

def total_spanning_tree_weight(G, weight=None):
    """
    Apply Kirchhoff's Tree Matrix Theorem a graph in order to find the total
    weight of all spanning trees.

    The theorem states that the determinant of any cofactor of the Laplacian
    matrix of a graph is the number of spanning trees in the graph. For a
    weighted Laplacian matrix, it is the sum across all spanning trees of the
    multiplicative weight of each tree. That is, the weight of each tree is the
    product of its edge weights.

    Parameters
    ----------
    G : NetworkX Graph
        The graph to use Kirchhoff's theorem on.

    weight : string or None
        The key for the edge attribute holding the edge weight. If `None`, then
        each edge is assumed to have a weight of 1.

    Returns
    -------
    float
        The sum of the total multiplicative weight for all spanning trees in the
        graph.
    """
    import numpy as np

    G_laplacian = nx.laplacian_matrix(G, weight=weight).toarray()
    # Determinant ignoring first row and column
    return abs(np.linalg.det(G_laplacian[1:, 1:]))

def spanning_tree_distribution(G, z):
    """
    Find the asadpour exponential distribution of spanning trees.

    Solves the Maximum Entropy Convex Program in the Asadpour algorithm [1]_
    using the approach in section 7 to build an exponential distribution of
    undirected spanning trees.

    This algorithm ensures that the probability of any edge in a spanning
    tree is proportional to the sum of the probabilities of the tress
    containing that edge over the sum of the probabilities of all spanning
    trees of the graph.

    Parameters
    ----------
    G : nx.MultiGraph
        The undirected support graph for the Held Karp relaxation

    z : dict
        The output of `held_karp_ascent()`, a scaled version of the Held-Karp
        solution.

    Returns
    -------
    gamma : dict
        The probability distribution which approximately preserves the marginal
        probabilities of `z`.
    """
    from math import exp
    from math import log as ln

    def q(e):
        """
        The value of q(e) is described in the Asadpour paper is "the
        probability that edge e will be included in a spanning tree T that is
        chosen with probability proportional to exp(gamma(T))" which
        basically means that it is the total probability of the edge appearing
        across the whole distribution.

        Parameters
        ----------
        e : tuple
            The `(u, v)` tuple describing the edge we are interested in

        Returns
        -------
        float
            The probability that a spanning tree chosen according to the
            current values of gamma will include edge `e`.
        """
        # Create the laplacian matrices
        for u, v, d in G.edges(data=True):
            d[lambda_key] = exp(gamma[(u, v)])
        G_Kirchhoff = total_spanning_tree_weight(G, lambda_key)
        G_e = nx.contracted_edge(G, e, self_loops=False)
        G_e_Kirchhoff = total_spanning_tree_weight(G_e, lambda_key)

        # Multiply by the weight of the contracted edge since it is not included
        # in the total weight of the contracted graph.
        return exp(gamma[(e[0], e[1])]) * G_e_Kirchhoff / G_Kirchhoff

    # initialize gamma to the zero dict
    gamma = {}
    for u, v, _ in G.edges:
        gamma[(u, v)] = 0

    # set epsilon
    EPSILON = 0.2

    # pick an edge attribute name that is unlikely to be in the graph
    lambda_key = "spanning_tree_distribution's secret attribute name for lambda"

    while True:
        # We need to know that know that no values of q_e are greater than
        # (1 + epsilon) * z_e, however changing one gamma value can increase the
        # value of a different q_e, so we have to complete the for loop without
        # changing anything for the condition to be meet
        in_range_count = 0
        # Search for an edge with q_e > (1 + epsilon) * z_e
        for u, v in gamma:
            e = (u, v)
            q_e = q(e)
            z_e = z[e]
            if q_e > (1 + EPSILON) * z_e:
                delta = ln(
                    (q_e * (1 - (1 + EPSILON / 2) * z_e))
                    / ((1 - q_e) * (1 + EPSILON / 2) * z_e)
                )
                gamma[e] -= delta
                # Check that delta had the desired effect
                new_q_e = q(e)
                desired_q_e = (1 + EPSILON / 2) * z_e
                if round(new_q_e, 8) != round(desired_q_e, 8):
                    raise nx.NetworkXError(
                        f"Unable to modify probability for edge ({u}, {v})"
                    )
            else:
                in_range_count += 1
        # Check if the for loop terminated without changing any gamma
        if in_range_count == len(gamma):
            break

    # Remove the new edge attributes
    for _, _, d in G.edges(data=True):
        if lambda_key in d:
            del d[lambda_key]

    return gamma

@py_random_state(2)
def sample_spanning_tree(G, lambda_key, seed=None):
    """
    Sample a spanning tree using the edges weights of the graph.

    The edge weights are multiplicative, so the probability of each tree is
    proportional to the product of edge weights.

    The algorithm itself uses algorithm A8 in [1]_ .

    We 'shuffle' the edges in the graph, and then probabilistically
    determine weather to add the edge conditioned on all of the previous
    edges which where added to the tree. Probabilities are calculated using
    Kirchhoff's Matrix Tree Theorem and a weighted Laplacian matrix.

    At each iteration, we contract the edges we have decided to include in the
    sampled tree and delete those which we have decided not to include.

    Parameters
    ----------
    G : nx.Graph
        An undirected version of the original graph.

    lambda_key : string
        The edge key for the edge attribute holding edge weight.

    seed : integer, random_state, or None (default)
        Indicator of random number generation state.
        See :ref:`Randomness<randomness>`.

    Returns
    -------
    nx.Graph
        A spanning tree using the distribution defined by `gamma`.

    References
    ----------
    .. [1] V. Kulkarni, Generating random combinatorial objects, Journal of
       algorithms, 11 (1990), pp. 185–207
    """

    def find_node(merged_nodes, n):
        """
        We can think of clusters of contracted nodes as having one
        representative in the graph. Each node which is not in merged_nodes
        is still its own representative. Since a representative can be later
        contracted, we need to recursively search though the dict to find
        the final representative, but once we know it we can use path
        compression to speed up the access of the representative for next time.

        Parameters
        ----------
        merged_nodes : dict
            The dict storing the mapping from node to representative
        n
            The node whose representative we seek

        Returns
        -------
        The representative of the `n`
        """
        if n not in merged_nodes:
            return n
        else:
            rep = find_node(merged_nodes, merged_nodes[n])
            merged_nodes[n] = rep
            return rep

    def prepare_graph():
        """
        For the graph `G`, remove all edges not in the set `V` and then
        contract all edges in the set `U`.

        Returns
        -------
        A copy of `G` which has had all edges not in `V` removed and all edges
        in `U` contracted.
        """

        # The result is a MultiGraph version of G so that parallel edges are
        # allowed during edge contraction
        result = nx.MultiGraph(incoming_graph_data=G)

        # Remove all edges not in V
        edges_to_remove = set(result.edges()).difference(V)
        result.remove_edges_from(edges_to_remove)

        # Contract all edges in U
        #
        # Imagine that you have two edges to contract and they share an
        # endpoint like this:
        #                        [0] ----- [1] ----- [2]
        # If we contract (0, 1) first, the contraction function will always
        # delete the second node it is passed so the resulting graph would be
        #                             [0] ----- [2]
        # and edge (1, 2) no longer exists but (0, 2) would need to be contracted
        # in its place now. That is why I use the below dict as a merge-find
        # data structure with path compression to track how the nodes are merged.
        merged_nodes = {}

        for u, v in U:
            u_rep = find_node(merged_nodes, u)
            v_rep = find_node(merged_nodes, v)
            # We cannot contract a node with itself
            if u_rep == v_rep:
                continue
            nx.contracted_nodes(result, u_rep, v_rep, self_loops=False, copy=False)
            merged_nodes[v_rep] = u_rep

        return merged_nodes, result

    U = set()
    V = set(G.edges())
    shuffled_edges = list(G.edges())
    seed.shuffle(shuffled_edges)

    for u, v in shuffled_edges:
        node_map, prepared_G = prepare_graph()
        G_total_tree_weight = total_spanning_tree_weight(prepared_G, lambda_key)
        # Add the edge to U so that we can compute the total tree weight
        # assuming we include that edge
        U.add((u, v))
        # Now, if (u, v) cannot exist in G because it is fully contracted out
        # of existence, then it by definition cannot influence G_e's Kirchhoff
        # value. But, we also cannot pick it.
        _, prepared_G_e = prepare_graph()
        rep_edge = (find_node(node_map, u), find_node(node_map, v))
        # Check to see if the 'representative edge' for the current edge is
        # in prepared_G. If so, then we can pick it.
        if rep_edge in prepared_G.edges:
            G_e_total_tree_weight = total_spanning_tree_weight(prepared_G_e, lambda_key)
        else:
            G_e_total_tree_weight = 0.0
        z = seed.uniform(0.0, 1.0)
        # This will be useful if I move this random spanning tree method to
        # the boarder NetworkX library
        e_weight = G[u][v][lambda_key] if lambda_key is not None else 1
        if z > e_weight * G_e_total_tree_weight / G_total_tree_weight:
            # Remove the edge from U since we did not decide to include it in
            # the sampled spanning tree. Also remove the edge from V because if
            # we did not decide to include it we must reject it.
            U.remove((u, v))
            V.remove((u, v))
        # If we decide to keep an edge, it may complete the spanning tree.
        elif len(U) == G.number_of_nodes() - 1:
            spanning_tree = nx.Graph()
            spanning_tree.add_edges_from(U)
            return spanning_tree

def greedy_tsp(G, weight="weight", source=None):
    """Return a low cost cycle starting at `source` and its cost.

    This approximates a solution to the traveling salesman problem.
    It finds a cycle of all the nodes that a salesman can visit in order
    to visit many nodes while minimizing total distance.
    It uses a simple greedy algorithm.
    In essence, this function returns a large cycle given a source point
    for which the total cost of the cycle is minimized.

    Parameters
    ----------
    G : Graph
        The Graph should be a complete weighted undirected graph.
        The distance between all pairs of nodes should be included.

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    source : node, optional (default: first node in list(G))
        Starting node.  If None, defaults to ``next(iter(G))``

    Returns
    -------
    cycle : list of nodes
        Returns the cycle (list of nodes) that a salesman
        can follow to minimize total weight of the trip.

    Raises
    ------
    NetworkXError
        If `G` is not complete, the algorithm raises an exception.

    Examples
    --------
    >>> from networkx.algorithms import approximation as approx
    >>> G = nx.DiGraph()
    >>> G.add_weighted_edges_from({
    ...     ("A", "B", 3), ("A", "C", 17), ("A", "D", 14), ("B", "A", 3),
    ...     ("B", "C", 12), ("B", "D", 16), ("C", "A", 13),("C", "B", 12),
    ...     ("C", "D", 4), ("D", "A", 14), ("D", "B", 15), ("D", "C", 2)
    ... })
    >>> cycle = approx.greedy_tsp(G, source="D")
    >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
    >>> cycle
    ['D', 'C', 'B', 'A', 'D']
    >>> cost
    31

    Notes
    -----
    This implementation of a greedy algorithm is based on the following:

    - The algorithm adds a node to the solution at every iteration.
    - The algorithm selects a node not already in the cycle whose connection
      to the previous node adds the least cost to the cycle.

    A greedy algorithm does not always give the best solution.
    However, it can construct a first feasible solution which can
    be passed as a parameter to an iterative improvement algorithm such
    as Simulated Annealing, or Threshold Accepting.

    Time complexity: It has a running time $O(|V|^2)$
    """
    # Check that G is a complete graph
    N = len(G) - 1
    # This check ignores selfloops which is what we want here.
    if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
        raise nx.NetworkXError("G must be a complete graph.")

    if source is None:
        source = nx.utils.arbitrary_element(G)

    if G.number_of_nodes() == 2:
        neighbor = next(G.neighbors(source))
        return [source, neighbor, source]

    nodeset = set(G)
    nodeset.remove(source)
    cycle = [source]
    next_node = source
    while nodeset:
        nbrdict = G[next_node]
        next_node = min(nodeset, key=lambda n: nbrdict[n].get(weight, 1))
        cycle.append(next_node)
        nodeset.remove(next_node)
    cycle.append(cycle[0])
    return cycle

@py_random_state(9)
def simulated_annealing_tsp(
    G,
    init_cycle,
    weight="weight",
    source=None,
    temp=100,
    move="1-1",
    max_iterations=10,
    N_inner=100,
    alpha=0.01,
    seed=None,
):
    """Returns an approximate solution to the traveling salesman problem.

    This function uses simulated annealing to approximate the minimal cost
    cycle through the nodes. Starting from a suboptimal solution, simulated
    annealing perturbs that solution, occasionally accepting changes that make
    the solution worse to escape from a locally optimal solution. The chance
    of accepting such changes decreases over the iterations to encourage
    an optimal result.  In summary, the function returns a cycle starting
    at `source` for which the total cost is minimized. It also returns the cost.

    The chance of accepting a proposed change is related to a parameter called
    the temperature (annealing has a physical analogue of steel hardening
    as it cools). As the temperature is reduced, the chance of moves that
    increase cost goes down.

    Parameters
    ----------
    G : Graph
        `G` should be a complete weighted undirected graph.
        The distance between all pairs of nodes should be included.

    init_cycle : list of all nodes or "greedy"
        The initial solution (a cycle through all nodes returning to the start).
        This argument has no default to make you think about it.
        If "greedy", use `greedy_tsp(G, weight)`.
        Other common starting cycles are `list(G) + [next(iter(G))]` or the final
        result of `simulated_annealing_tsp` when doing `threshold_accepting_tsp`.

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    source : node, optional (default: first node in list(G))
        Starting node.  If None, defaults to ``next(iter(G))``

    temp : int, optional (default=100)
        The algorithm's temperature parameter. It represents the initial
        value of temperature

    move : "1-1" or "1-0" or function, optional (default="1-1")
        Indicator of what move to use when finding new trial solutions.
        Strings indicate two special built-in moves:

        - "1-1": 1-1 exchange which transposes the position
          of two elements of the current solution.
          The function called is :func:`swap_two_nodes`.
          For example if we apply 1-1 exchange in the solution
          ``A = [3, 2, 1, 4, 3]``
          we can get the following by the transposition of 1 and 4 elements:
          ``A' = [3, 2, 4, 1, 3]``
        - "1-0": 1-0 exchange which moves an node in the solution
          to a new position.
          The function called is :func:`move_one_node`.
          For example if we apply 1-0 exchange in the solution
          ``A = [3, 2, 1, 4, 3]``
          we can transfer the fourth element to the second position:
          ``A' = [3, 4, 2, 1, 3]``

        You may provide your own functions to enact a move from
        one solution to a neighbor solution. The function must take
        the solution as input along with a `seed` input to control
        random number generation (see the `seed` input here).
        Your function should maintain the solution as a cycle with
        equal first and last node and all others appearing once.
        Your function should return the new solution.

    max_iterations : int, optional (default=10)
        Declared done when this number of consecutive iterations of
        the outer loop occurs without any change in the best cost solution.

    N_inner : int, optional (default=100)
        The number of iterations of the inner loop.

    alpha : float between (0, 1), optional (default=0.01)
        Percentage of temperature decrease in each iteration
        of outer loop

    seed : integer, random_state, or None (default)
        Indicator of random number generation state.
        See :ref:`Randomness<randomness>`.

    Returns
    -------
    cycle : list of nodes
        Returns the cycle (list of nodes) that a salesman
        can follow to minimize total weight of the trip.

    Raises
    ------
    NetworkXError
        If `G` is not complete the algorithm raises an exception.

    Examples
    --------
    >>> from networkx.algorithms import approximation as approx
    >>> G = nx.DiGraph()
    >>> G.add_weighted_edges_from({
    ...     ("A", "B", 3), ("A", "C", 17), ("A", "D", 14), ("B", "A", 3),
    ...     ("B", "C", 12), ("B", "D", 16), ("C", "A", 13),("C", "B", 12),
    ...     ("C", "D", 4), ("D", "A", 14), ("D", "B", 15), ("D", "C", 2)
    ... })
    >>> cycle = approx.simulated_annealing_tsp(G, "greedy", source="D")
    >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
    >>> cycle
    ['D', 'C', 'B', 'A', 'D']
    >>> cost
    31
    >>> incycle = ["D", "B", "A", "C", "D"]
    >>> cycle = approx.simulated_annealing_tsp(G, incycle, source="D")
    >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
    >>> cycle
    ['D', 'C', 'B', 'A', 'D']
    >>> cost
    31

    Notes
    -----
    Simulated Annealing is a metaheuristic local search algorithm.
    The main characteristic of this algorithm is that it accepts
    even solutions which lead to the increase of the cost in order
    to escape from low quality local optimal solutions.

    This algorithm needs an initial solution. If not provided, it is
    constructed by a simple greedy algorithm. At every iteration, the
    algorithm selects thoughtfully a neighbor solution.
    Consider $c(x)$ cost of current solution and $c(x')$ cost of a
    neighbor solution.
    If $c(x') - c(x) <= 0$ then the neighbor solution becomes the current
    solution for the next iteration. Otherwise, the algorithm accepts
    the neighbor solution with probability $p = exp - ([c(x') - c(x)] / temp)$.
    Otherwise the current solution is retained.

    `temp` is a parameter of the algorithm and represents temperature.

    Time complexity:
    For $N_i$ iterations of the inner loop and $N_o$ iterations of the
    outer loop, this algorithm has running time $O(N_i * N_o * |V|)$.

    For more information and how the algorithm is inspired see:
    http://en.wikipedia.org/wiki/Simulated_annealing
    """
    if move == "1-1":
        move = swap_two_nodes
    elif move == "1-0":
        move = move_one_node
    if init_cycle == "greedy":
        # Construct an initial solution using a greedy algorithm.
        cycle = greedy_tsp(G, weight=weight, source=source)
        if G.number_of_nodes() == 2:
            return cycle

    else:
        cycle = list(init_cycle)
        if source is None:
            source = cycle[0]
        elif source != cycle[0]:
            raise nx.NetworkXError("source must be first node in init_cycle")
        if cycle[0] != cycle[-1]:
            raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")

        if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
            raise nx.NetworkXError("init_cycle should be a cycle over all nodes in G.")

        # Check that G is a complete graph
        N = len(G) - 1
        # This check ignores selfloops which is what we want here.
        if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
            raise nx.NetworkXError("G must be a complete graph.")

        if G.number_of_nodes() == 2:
            neighbor = next(G.neighbors(source))
            return [source, neighbor, source]

    # Find the cost of initial solution
    cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))

    count = 0
    best_cycle = cycle.copy()
    best_cost = cost
    while count <= max_iterations and temp > 0:
        count += 1
        for i in range(N_inner):
            adj_sol = move(cycle, seed)
            adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
            delta = adj_cost - cost
            if delta <= 0:
                # Set current solution the adjacent solution.
                cycle = adj_sol
                cost = adj_cost

                if cost < best_cost:
                    count = 0
                    best_cycle = cycle.copy()
                    best_cost = cost
            else:
                # Accept even a worse solution with probability p.
                p = math.exp(-delta / temp)
                if p >= seed.random():
                    cycle = adj_sol
                    cost = adj_cost
        temp -= temp * alpha

    return best_cycle

@py_random_state(9)
def threshold_accepting_tsp(
    G,
    init_cycle,
    weight="weight",
    source=None,
    threshold=1,
    move="1-1",
    max_iterations=10,
    N_inner=100,
    alpha=0.1,
    seed=None,
):
    """Returns an approximate solution to the traveling salesman problem.

    This function uses threshold accepting methods to approximate the minimal cost
    cycle through the nodes. Starting from a suboptimal solution, threshold
    accepting methods perturb that solution, accepting any changes that make
    the solution no worse than increasing by a threshold amount. Improvements
    in cost are accepted, but so are changes leading to small increases in cost.
    This allows the solution to leave suboptimal local minima in solution space.
    The threshold is decreased slowly as iterations proceed helping to ensure
    an optimum. In summary, the function returns a cycle starting at `source`
    for which the total cost is minimized.

    Parameters
    ----------
    G : Graph
        `G` should be a complete weighted undirected graph.
        The distance between all pairs of nodes should be included.

    init_cycle : list or "greedy"
        The initial solution (a cycle through all nodes returning to the start).
        This argument has no default to make you think about it.
        If "greedy", use `greedy_tsp(G, weight)`.
        Other common starting cycles are `list(G) + [next(iter(G))]` or the final
        result of `simulated_annealing_tsp` when doing `threshold_accepting_tsp`.

    weight : string, optional (default="weight")
        Edge data key corresponding to the edge weight.
        If any edge does not have this attribute the weight is set to 1.

    source : node, optional (default: first node in list(G))
        Starting node.  If None, defaults to ``next(iter(G))``

    threshold : int, optional (default=1)
        The algorithm's threshold parameter. It represents the initial
        threshold's value

    move : "1-1" or "1-0" or function, optional (default="1-1")
        Indicator of what move to use when finding new trial solutions.
        Strings indicate two special built-in moves:

        - "1-1": 1-1 exchange which transposes the position
          of two elements of the current solution.
          The function called is :func:`swap_two_nodes`.
          For example if we apply 1-1 exchange in the solution
          ``A = [3, 2, 1, 4, 3]``
          we can get the following by the transposition of 1 and 4 elements:
          ``A' = [3, 2, 4, 1, 3]``
        - "1-0": 1-0 exchange which moves an node in the solution
          to a new position.
          The function called is :func:`move_one_node`.
          For example if we apply 1-0 exchange in the solution
          ``A = [3, 2, 1, 4, 3]``
          we can transfer the fourth element to the second position:
          ``A' = [3, 4, 2, 1, 3]``

        You may provide your own functions to enact a move from
        one solution to a neighbor solution. The function must take
        the solution as input along with a `seed` input to control
        random number generation (see the `seed` input here).
        Your function should maintain the solution as a cycle with
        equal first and last node and all others appearing once.
        Your function should return the new solution.

    max_iterations : int, optional (default=10)
        Declared done when this number of consecutive iterations of
        the outer loop occurs without any change in the best cost solution.

    N_inner : int, optional (default=100)
        The number of iterations of the inner loop.

    alpha : float between (0, 1), optional (default=0.1)
        Percentage of threshold decrease when there is at
        least one acceptance of a neighbor solution.
        If no inner loop moves are accepted the threshold remains unchanged.

    seed : integer, random_state, or None (default)
        Indicator of random number generation state.
        See :ref:`Randomness<randomness>`.

    Returns
    -------
    cycle : list of nodes
        Returns the cycle (list of nodes) that a salesman
        can follow to minimize total weight of the trip.

    Raises
    ------
    NetworkXError
        If `G` is not complete the algorithm raises an exception.

    Examples
    --------
    >>> from networkx.algorithms import approximation as approx
    >>> G = nx.DiGraph()
    >>> G.add_weighted_edges_from({
    ...     ("A", "B", 3), ("A", "C", 17), ("A", "D", 14), ("B", "A", 3),
    ...     ("B", "C", 12), ("B", "D", 16), ("C", "A", 13),("C", "B", 12),
    ...     ("C", "D", 4), ("D", "A", 14), ("D", "B", 15), ("D", "C", 2)
    ... })
    >>> cycle = approx.threshold_accepting_tsp(G, "greedy", source="D")
    >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
    >>> cycle
    ['D', 'C', 'B', 'A', 'D']
    >>> cost
    31
    >>> incycle = ["D", "B", "A", "C", "D"]
    >>> cycle = approx.threshold_accepting_tsp(G, incycle, source="D")
    >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
    >>> cycle
    ['D', 'C', 'B', 'A', 'D']
    >>> cost
    31

    Notes
    -----
    Threshold Accepting is a metaheuristic local search algorithm.
    The main characteristic of this algorithm is that it accepts
    even solutions which lead to the increase of the cost in order
    to escape from low quality local optimal solutions.

    This algorithm needs an initial solution. This solution can be
    constructed by a simple greedy algorithm. At every iteration, it
    selects thoughtfully a neighbor solution.
    Consider $c(x)$ cost of current solution and $c(x')$ cost of
    neighbor solution.
    If $c(x') - c(x) <= threshold$ then the neighbor solution becomes the current
    solution for the next iteration, where the threshold is named threshold.

    In comparison to the Simulated Annealing algorithm, the Threshold
    Accepting algorithm does not accept very low quality solutions
    (due to the presence of the threshold value). In the case of
    Simulated Annealing, even a very low quality solution can
    be accepted with probability $p$.

    Time complexity:
    It has a running time $O(m * n * |V|)$ where $m$ and $n$ are the number
    of times the outer and inner loop run respectively.

    For more information and how algorithm is inspired see:
    https://doi.org/10.1016/0021-9991(90)90201-B

    See Also
    --------
    simulated_annealing_tsp

    """
    if move == "1-1":
        move = swap_two_nodes
    elif move == "1-0":
        move = move_one_node
    if init_cycle == "greedy":
        # Construct an initial solution using a greedy algorithm.
        cycle = greedy_tsp(G, weight=weight, source=source)
        if G.number_of_nodes() == 2:
            return cycle

    else:
        cycle = list(init_cycle)
        if source is None:
            source = cycle[0]
        elif source != cycle[0]:
            raise nx.NetworkXError("source must be first node in init_cycle")
        if cycle[0] != cycle[-1]:
            raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")

        if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
            raise nx.NetworkXError("init_cycle is not all and only nodes.")

        # Check that G is a complete graph
        N = len(G) - 1
        # This check ignores selfloops which is what we want here.
        if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
            raise nx.NetworkXError("G must be a complete graph.")

        if G.number_of_nodes() == 2:
            neighbor = list(G.neighbors(source))[0]
            return [source, neighbor, source]

    # Find the cost of initial solution
    cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))

    count = 0
    best_cycle = cycle.copy()
    best_cost = cost
    while count <= max_iterations:
        count += 1
        accepted = False
        for i in range(N_inner):
            adj_sol = move(cycle, seed)
            adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
            delta = adj_cost - cost
            if delta <= threshold:
                accepted = True

                # Set current solution the adjacent solution.
                cycle = adj_sol
                cost = adj_cost

                if cost < best_cost:
                    count = 0
                    best_cycle = cycle.copy()
                    best_cost = cost
        if accepted:
            threshold -= threshold * alpha

    return best_cycle

ext: py

max_stars_repo_name: amcandio/networkx

max_issues_repo_name: amcandio/networkx

max_forks_repo_name: amcandio/networkx

reason: import numpy,import scipy,import networkx,from networkx


text:
#!/usr/bin/env python3
# Copyright (c) Facebook, Inc. and its affiliates.

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import logging

import numpy as np
from ..processor.commons import AbstractProcessor
from ..processor.matrix import get_result

class LinearRegressionCoefficients:
    """
    Internal class to be used by LinearRegression
    You should not call this class directly.
    """

    def __init__(self, n, xtx, xty, yty, x_s2, df_corrected=True):
        if xtx.shape[0] != xtx.shape[1]:
            raise ValueError("xtx should be square matrix!")
        if xtx.shape[0] != xty.shape[0]:
            raise ValueError("xtx and xty have unmatched dimension")
        if any(s2 < 0.0 for s2 in x_s2):
            raise ValueError("x_s2 should not contain negative elements")
        self.n = n
        self.xtx = xtx
        self.xty = xty
        self.yty = yty
        self.x_s2 = x_s2
        self.df_corrected = df_corrected
        self.k = self.xtx.shape[0]
        self.beta = None
        self.sigma_sq = None
        self.omega = None

    def estimate_beta(self):
        self.omega = self.xtx.copy()
        self.omega[np.diag_indices(self.k)] = np.diag(self.omega) - self.x_s2

        # check whether omega is positive-definite or not
        # reference on threshold https://fburl.com/7y4wyaev
        eigenvalues = np.linalg.eigvals(self.omega)
        if eigenvalues.max() < 0 or (
            eigenvalues.min() < eigenvalues.max() * self.k * np.finfo(np.float_).eps
        ):
            logging.warning("omega is not positive definite")
            self.beta = np.linalg.lstsq(self.omega, self.xty, rcond=None)[0]
        else:
            self.beta = np.linalg.solve(self.omega, self.xty)

    def estimate_residual_var(self):
        if self.yty is None:
            self.sigma_sq = None
            return
        df = self.n
        if self.df_corrected:
            # I think it should be df -= self.k
            # but use df -= self.k - 1 for now to be aligned with the R package
            df -= self.k - 1
        term1 = (
            self.yty
            - 2 * self.beta.T.dot(self.xty)
            + self.beta.T.dot(self.xtx).dot(self.beta)
        ) * self.n / df
        term2 = (self.beta ** 2 * self.x_s2).sum()
        sigma_sq = term1 - term2
        # check whether variance is positive or not
        if sigma_sq < 0.0:
            logging.warning("negative residual variance!")
        self.sigma_sq = sigma_sq

    def estimate_all(self):
        self.estimate_beta()
        self.estimate_residual_var()
        return self.beta, self.sigma_sq, self.omega

class LinearRegressionVariance:
    """
    Internal class to be used by LinearRegressionFit
    You should not call this class directly.
    """

    def __init__(
        self,
        n,
        xtx,
        xty,
        yty,
        sigma_sq,
        omega,
        x_s2,
        n_replications=500,
        random_state=None,
    ):
        self.n = n
        self.xtx = xtx
        self.xty = xty
        self.yty = yty
        self.sigma_sq = sigma_sq
        self.omega = omega
        self.x_s2 = x_s2
        self.n_replications = n_replications
        self.random_state = random_state
        self.k = self.xtx.shape[0]

    def estimate_vcov_xx_xx(self, k, l, j, m):
        return (
            self.omega[k, l] * self.x_s2[j] * (j == m)
            + self.omega[k, m] * self.x_s2[j] * (j == l)
            + self.omega[j, l] * self.x_s2[k] * (k == m)
            + self.omega[j, m] * self.x_s2[k] * (k == l)
            + self.x_s2[k] * (k == l) * self.x_s2[j] * (j == m)
            + self.x_s2[k] * (k == m) * self.x_s2[j] * (j == l)
        ) / self.n

    def estimate_vcov_xy_xy(self, k, j):
        return (
            self.sigma_sq * self.omega[k, j]
            + self.x_s2[k] * (k == j) * self.yty
        ) / self.n

    def estimate_vcov_xy_xx(self, k, j, m):
        return (
            self.x_s2[k] * (k == m) * self.xty[j]
            + self.x_s2[k] * (k == j) * self.xty[m]
        ) / self.n

    def simulate_distribution(self):
        index = np.triu_indices(self.k)
        t = np.concatenate((self.xtx[index], self.xty), axis=0)

        d1, d2 = self.k * (self.k + 1) // 2, self.k
        v_t = np.zeros((d1 + d2, d1 + d2))
        for i1 in range(d1 + d2):
            for i2 in range(i1 + 1):
                if i1 < d1 and i2 < d1:
                    k, j = index[0][i1], index[1][i1]
                    l, m = index[0][i2], index[1][i2]
                    v_t[i1, i2] = self.estimate_vcov_xx_xx(k, l, j, m)
                elif i2 < d1 <= i1:
                    k = i1 - d1
                    j, m = index[0][i2], index[1][i2]
                    v_t[i1, i2] = self.estimate_vcov_xy_xx(k, j, m)
                else:
                    k = i1 - d1
                    j = i2 - d1
                    v_t[i1, i2] = self.estimate_vcov_xy_xy(k, j)
                v_t[i2, i1] = v_t[i1, i2]

        return t, v_t

    def transform_vector_to_matrix(self, sample):
        k = self.k
        d = k * (k + 1) // 2
        xty = sample[d:]
        xtx = np.zeros((k, k))
        xtx[np.triu_indices(k)] = sample[:d]
        i, j = np.triu_indices(k, 1)
        xtx[(j, i)] = xtx[(i, j)]
        return xtx, xty

    def simulate_beta_vcov(self):
        t, v_t = self.simulate_distribution()
        np.random.seed(self.random_state)
        t_samples = np.random.multivariate_normal(
            t, v_t, self.n_replications, check_valid="ignore"
        )

        simu_beta_list = []
        for i in range(self.n_replications):
            simu_xtx, simu_xty = self.transform_vector_to_matrix(t_samples[i, :])
            simu_beta, _, _ = LinearRegressionCoefficients(
                self.n, simu_xtx, simu_xty, None, self.x_s2
            ).estimate_all()
            simu_beta_list.append(simu_beta)

        beta_vcov = np.cov(np.array(simu_beta_list), rowvar=False)
        return beta_vcov

class LinearRegression:
    def __init__(
        self,
        x_columns,
        y_column,
        x_s2,
        fit_intercept=True,
        df_corrected=True,
        n_replications=500,
        random_state=None,
    ):
        """
        LinearRegression fits a linear regression where some variables
        in the training data are subject to naturally occurring measurement
        error, and implements a computationally efficient algorithm for
        statistical inference in large datasets.

        :param x_columns: list of strings.
        Specify the column names for the features in the data.
        :param y_column: string
        Specify the column name for the response in the data.
        :param x_s2: list of doubles
        Specify the variance of noise that was added to each of the features.
        It’s length should be equal to that of x_columns. If no noise is
        added to a feature, use 0.0 in the corresponding place.
        :param fit_intercept: bool, optional, default to True
        Whether to include the intercept into the model. If set to False,
        the design matrix will only consist of the features specified.
        If set to True, a column of 1 will be added to the design matrix
        besides features specified.
        :param df_corrected: bool, optional, default to True
        Whether to adjust the degree of freedom when estimating the error
        variance. If set to False, the degree of freedom is n, where n is
        the sample size. If set to True, the degree of freedom is n - p,
        where p is the number of columns in the design matrix.
        :param n_replications: int, optional, default to 500
        The number of simulation replicates. It should be a positive integer.
        :param random_state: int or None, optional, default to None
        Determines random number generation for dataset creation. Pass an
        int for reproducible output across multiple function calls.
        """
        self.x_columns = x_columns
        self.y_column = y_column
        self.x_s2 = ([0.0] if fit_intercept else []) + x_s2
        self.fit_intercept = fit_intercept
        self.df_corrected = df_corrected
        self.n_replications = n_replications
        self.random_state = random_state

        self.beta = None
        self.sigma_sq = None
        self.beta_vcov = None
        self.beta_standarderror = None

    def _preprocess_data(self, data: AbstractProcessor):
        logging.info("data is an instance of %s", type(data))
        x, y = data.prepare_xy(self.x_columns, self.y_column, self.fit_intercept)
        z = get_result({
            "xtx": x.cross(x),
            "xty": x.cross(y),
            "yty": y.cross(y),
        }, data.run_query)
        n = z["sample_size"]
        xtx = z["xtx"]
        xty = z["xty"][:, 0]
        yty = z["yty"][0, 0]
        return n, xtx, xty, yty

    def fit(self, data):
        """
        fit model
        :return: a tuple; the 1st element is the estimated beta,
        the 2nd element is the estimated sigma squared, the 3rd element is
        the estimated variance-covariance matrix of the estimated beta.
        """
        n, xtx, xty, yty = self._preprocess_data(data)
        self.beta, self.sigma_sq, omega = LinearRegressionCoefficients(
            n, xtx, xty, yty, self.x_s2, df_corrected=self.df_corrected
        ).estimate_all()

        self.beta_vcov = LinearRegressionVariance(
            n,
            xtx,
            xty,
            yty,
            self.sigma_sq,
            omega,
            self.x_s2,
            n_replications=self.n_replications,
            random_state=self.random_state,
        ).simulate_beta_vcov()

        self.beta_standarderror = np.sqrt(np.diag(self.beta_vcov))
        return self

max_stars_repo_name: facebookresearch/svinfer

max_issues_repo_name: facebookresearch/svinfer

max_forks_repo_name: facebookresearch/svinfer

reason: import numpy


text:
import pandas as pd
import warnings
from featuregenerator import Featureengineering
from utils import get_cat_cont_feats,set_seed,split_dataset,cat_dim
from dataset import ClassificationDataset
from model import ClassificationEmbdNN
from engine import train_network
import torch
import numpy as np
from predict import predict_pipe

warnings.filterwarnings('ignore')

def run():
    train=pd.read_csv('../input/train.csv')
    test=pd.read_csv('../input/test.csv')

    traintest=pd.concat([train,test],ignore_index=True)
    traintest=traintest.drop('id',axis=1)

    feateng=Featureengineering(traintest,cyc_feats=['month','day'],bin_feats=['bin_0', 'bin_1', 'bin_2', 'bin_3', 'bin_4'],
                                ord_feats=['ord_0', 'ord_1', 'ord_2', 'ord_3', 'ord_4', 'ord_5'],nom_feats=['nom_0', 'nom_1',
                                'nom_2', 'nom_3', 'nom_4', 'nom_5', 'nom_6', 'nom_7', 'nom_8', 'nom_9'])

    continuous=['month_sin', 'month_cos','day_sin', 'day_cos']

    dataset=feateng.engineer()

    train=dataset[dataset['target']!=-1]
    test=dataset[dataset['target']==-1]

    traindata=train.drop(['target'],axis=1)
    target=train['target']

    testdata=test.drop(['target'],axis=1)

    print(traindata.shape,testdata.shape)

    cats,conts=get_cat_cont_feats(train,continuous)

    set_seed(2020)

    train_dataset = ClassificationDataset(traindata,
                                      targets=target,
                                      cat_cols=cats,
                                      cont_cols=conts)
    test_dataset = ClassificationDataset(testdata,
                                     cat_cols=cats,
                                     cont_cols=conts,
                                     is_train=False)

    train_loader, valid_loader = split_dataset(train_dataset,
                                           valid_size=0.2,batch_size=128)

    cat_dims=cat_dim(train,cats)

    model = ClassificationEmbdNN(emb_dims=cat_dims,
                             no_of_cont=len(continuous))

    loss_func = torch.nn.BCELoss()
    optimizer = torch.optim.Adam(model.parameters(), lr=0.1)

    train_losses, valid_losses = train_network(model=model,
                                           train_loader=train_loader,
                                           valid_loader=valid_loader,
                                           loss_func=loss_func,
                                           optimizer=optimizer,
                                           n_epochs=3,
                                           saved_model='../output/simple_nn.pt')

    return model,test_dataset

if __name__=='__main__':
    model,test=run()
    inference=True

    if inference:
        predictions=predict_pipe(model,test)
        sumbission=pd.read_csv('../input/sample_submission.csv')
        nn_predictions_df = pd.DataFrame({'id': sumbission['id'], 'target': predictions.squeeze()})
        nn_predictions_df.to_csv('../output/submission.csv',index=None)