def METHOD_NAME(args):
    data_np = pickle.load(open(args.data, "rb"))
    data = torch.from_numpy(data_np)
    if args.max_images is not None:
        data = data[: args.max_images]
    use_cuda = torch.cuda.is_available()
    print(f"Use cuda {use_cuda}")
    device = torch.device("cuda" if use_cuda else "cpu")
    data = data.to(device)
    N, D = data.shape
    for i in args.anchors:
        assert i < N
    assert len(args.anchors) >= 2
    n2 = (data * data).sum(-1, keepdim=True)
    B = args.batch_size
    ndist = torch.empty(data.shape[0], args.max_neighbors, device=device)
    neighbors = torch.empty(
        data.shape[0], args.max_neighbors, dtype=torch.long, device=device
    )
    for i in range(0, data.shape[0], B):
        # (a-b)^2 = a^2 + b^2 - 2ab
        print(f"Working on images {i}-{i+B}")
        batch_dist = n2[i : i + B] + n2.t() - 2 * torch.mm(data[i : i + B], data.t())
        ndist[i : i + B], neighbors[i : i + B] = batch_dist.topk(
            args.max_neighbors, dim=-1, largest=False
        )
    assert ndist.min() >= -1e-3, ndist.min()
    # convert d^2 to d
    ndist = ndist.clamp(min=0).pow(0.5)
    if args.avg_neighbors:
        total_neighbors = int(N * args.avg_neighbors)
        max_dist = ndist.view(-1).topk(total_neighbors, largest=False)[0][-1]
    else:
        max_dist = None
    print(
        f"Max dist between neighbors: {max_dist}  (to enforce average of {args.avg_neighbors} neighbors)"
    )
    if max_dist is not None:
        max_dist = max_dist.to("cpu")
    neighbors = neighbors.to("cpu")
    ndist = ndist.to("cpu")
    edges = []
    for i in range(neighbors.shape[0]):
        for j in range(neighbors.shape[1]):
            if max_dist is None or ndist[i, j] < max_dist:
                edges.append((int(i), int(neighbors[i, j]), float(ndist[i, j])))
    graph = Graph(edges)
    full_path = []
    for i in range(len(args.anchors) - 1):
        src, dest = args.anchors[i], args.anchors[i + 1]
        path, total_distance = graph.find_path(src, dest)
        dd = data[path].cpu().numpy()
        dists = ((dd[1:, :] - dd[0:-1, :]) ** 2).sum(axis=1) ** 0.5
        if path is not None:
            if full_path and full_path[-1] == path[0]:
                full_path.extend(path[1:])
            else:
                full_path.extend(path)
        print()
        if path is not None:
            print("Path:")
            for id in path:
                print(id)
            print()
            print(f"Total distance: {total_distance}")
            print()
            print("Neighbor distance:")
            for d in dists:
                print(d)
            print()
            print("Euclidean distance: {}".format(((dd[0] - dd[-1]) ** 2).sum() ** 0.5))
        else:
            print("Could not find path!")
    if not os.path.exists(os.path.dirname(args.o)):
        os.makedirs(os.path.dirname(args.o))
    print(args.o)
    np.savetxt(args.o, full_path, fmt="%d")
    print(args.out_z)
    np.savetxt(args.out_z, data_np[full_path])