text:
from __future__ import print_function
from builtins import range
import random
import pandas as pd
import time
import math
import numpy as np
from copy import deepcopy
import itertools

class GeneticAlgorithm(object):
    """
    A python implementation of real-valued, genetic algorithm for solving optimization problems.

    Parameters
    ----------
    evaluate : function
        The objective function that has to be optimized. The first parameter of the objective function is a list of the trial values of the hyper-parameters in the order in which they are declared in the space variable. The objective function should always return a tuple with the metric/metrics for single/multi-objective optimization.

    space : tuple,
        A tuple of dict objects specifying the hyper-parameter space to search in.
        Each hyper-parameter should be a python dict object with the name of the hyper-parameter as the key.
        Value is also a dict object with one mandatory key among: 'uniform', 'int' and 'choice' for defining floating point, integer and choice variables respectively.
        Values for these keys should be a list defining the valid hyper-parameter search space (lower and upper bounds for 'int' and 'uniform', and all valid choices for 'choice').
        For uniform, a 'mutation' key is also required for which the value is [mean, standard deviation] for the gaussian distribution.
        Example:
                ({'alpha': {'uniform': [0.001, 1],
                            'mutation': [0, 1]}},
                {'layers': {'int': [1, 3]}},
                {'neurons': {'choice': range(0,200,20)}})

    fitness : tuple, optional (default = ('Max',)
        A tuple of string(s) for Maximizing (Max) or minimizing (Min) the objective function(s).

    pop_size : integer, optional (default = 50)
        Size of the population

    crossover_size : int, optional (default = 30)
        Number of individuals to select for crossover.

    mutation_size : int, optional (default = 20)
        Number of individuals to select for mutation.

    crossover_type : string, optional (default = "Blend")
        Type of crossover: SinglePoint, DoublePoint, Blend, Uniform

    mutation_prob : float, optional (default = 0.4)
        Probability of mutation.

    algorithm : int, optional (default=1)
        The algorithm to use for the search. Look at the 'search' method for a description of the various algorithms.

            - Algorithm 1:
                Initial population is instantiated.
                Roulette wheel selection is used for selecting individuals for crossover and mutation.
                The initial population, crossovered and mutated individuals form the pool of individuals from which the best
                n members are selected as the initial population for the next generation, where n is the size of population.

            - Algorithm 2:
                Same as algorithm 1 but when selecting individuals for next generation, n members are selected using Roulette wheel selection.

            - Algorithm 3:
                Same as algorithm 1 but when selecting individuals for next generation, best members from each of the three pools (initital population, crossover and mutation) are selected according to the input parameters in the search method.

            - Algorithm 4:
                Same as algorithm 1 but mutation population is selected from the crossover population and not from the parents directly.

    initial_population : list, optional (default=None)
        The initial population for the algorithm to start with. If not provided, initial population is randomly generated.

    """

    def __init__(self,
                evaluate,
                space,
                fitness=("Max", ),
                pop_size=50,
                crossover_size=30,
                mutation_size=20,
                crossover_type="Blend",
                fused_cutoff = 5,
                mutation_prob=0.6,
                algorithm=3,
                initial_population=None):

        self.chromosome_length = len(space)
        if self.chromosome_length < 1:
            print("Space variable not defined. Aborting.")
            exit(code=1)
        if self.chromosome_length <2 and crossover_type == "SinglePoint": raise Exception('Single point crossover not possible for chromosome length 1.')
        if self.chromosome_length <3 and crossover_type == "DoublePoint": raise Exception('Double point crossover not possible for chromosome length 2.')
        self.chromosome_type, self.bit_limits, self.mutation_params, self.var_names = [], [], [], []

        uni = 0
        for param_dict in space:
            for name in param_dict:
                self.var_names.append(name)
                var = param_dict[name]
                if 'uniform' in var:
                    self.chromosome_type.append('uniform')
                    self.bit_limits.append(var['uniform'])
                    uni += 1

                elif 'int' in var:
                    self.chromosome_type.append('int')
                    self.bit_limits.append(var['int'])

                elif 'choice' in var:
                    self.chromosome_type.append('choice')
                    self.bit_limits.append(var['choice'])

                if 'mutation' in var:
                    self.mutation_params.append(var['mutation'])
                else: self.mutation_params.append(None)

        self.evaluate = evaluate
        self.pop_size = pop_size
        self.mutation_prob = mutation_prob
        self.crossover_type = crossover_type
        self.fused_cutoff = fused_cutoff
        self.crossover_size = crossover_size
        self.mutation_size = mutation_size
        self.algo = algorithm
        self.initial_pop = initial_population
        self.fit_val, self.population, self.fitness_dict, self.global_cm_list = [], None, {}, None
        for i in fitness:
            if i.lower() == 'max': self.fit_val.append(1)
            else: self.fit_val.append(-1)
        # if there is no uniform type in the space parameter, use a pre-defined crossover and mutation list as follows:
        if uni == 0:
            if len(space) == 1: self.global_cm_list = self.bit_limits[0]
            else:
                gcl = []
                for i, t in zip(self.bit_limits, self.chromosome_type):
                    if t == 'int':
                        gcl.append(list(range(i[0], i[1]+1)))
                    else: gcl.append(i)
                self.global_cm_list = list(itertools.product(*gcl))

    def pop_generator(self, n):
        pop = []
        if self.initial_pop is not None:
            for i in self.initial_pop:
                pop.append(tuple(i))
        else:
            for x in range(n):
                while True:
                    new = self.chromosome_generator(x, n)
                    if new not in pop:
                        pop.append(new)
                        break
        return pop

    def chromosome_generator(self, x, n):
        chsome = []
        for i, j in zip(self.bit_limits, self.chromosome_type):
            if j == 'uniform':
                chsome.append(np.linspace(i[0], i[1], n)[x])
            elif j == 'int':
                chsome.append(random.randint(i[0], i[1]))
            elif j == 'choice':
                chsome.append(random.choice(i))
        return tuple(chsome)

    def SinglePointCrossover(self, x1, x2):
        x1, x2 = list(x1), list(x2)
        nVar=len(x1)
        c = random.randint(1,nVar-1)
        y1=x1[0:c]
        y1=y1+x2[c:nVar]
        y2=x2[0:c]
        y2=y2+x1[c:nVar]
        return tuple(deepcopy(y1)), tuple(deepcopy(y2))

    def DoublePointCrossover(self, x1, x2):
        x1, x2 = list(x1), list(x2)
        nVar = len(x1)
        cc = random.sample(range(1,nVar), 2)
        c1 = min(cc)
        c2 = max(cc)
        y1 = x1[0:c1]+x2[c1:c2]+x1[c2:nVar]
        y2 = x2[0:c1]+x1[c1:c2]+x2[c2:nVar]
        return tuple(deepcopy(y1)), tuple(deepcopy(y2))

    def UniformCrossover(self, x1, x2):
        parents = [x1,x2]
        ind1, ind2 = [parents[random.randint(0, 1)][i] for i in range(len(x1))], [parents[random.randint(0, 1)][i] for i in range(len(x1))]
        return tuple(deepcopy(ind1)), tuple(deepcopy(ind2))

    def blend(self, ind1, ind2, fitness_dict, z=0.4, alpha=0.5, beta=0.1):
        # rank all individuals in fitness dict
        ranked_list = self.select([i for i in fitness_dict], fitness_dict, len(fitness_dict.items()), choice="best")
        # determine the better individual among the two for implementing the alpha-beta crossover
        if ranked_list.index(ind1) < ranked_list.index(ind2): better, worse = deepcopy(list(ind1)), deepcopy(list(ind2))
        else: better, worse = deepcopy(list(ind2)), deepcopy(list(ind1))

        for i in range(self.chromosome_length):
            # for type 'choice' apply roulette wheel selection by biasing the wheel towards the two current values
            if self.chromosome_type[i] == 'choice':
                if len(self.bit_limits[i]) == 2: better[i], worse[i] = worse[i], better[i]
                else:
                    scores = []
                    for inx, ch in enumerate(self.bit_limits[i]):
                        if ch == better[i] or ch == worse[i]: scores.append(3)
                        else: scores.append(1)
                    sc_dict = {self.bit_limits[i][j]: scores[j] for j in range(len(scores))}
                    better[i], worse[i] = self.select(self.bit_limits[i], sc_dict, 2)
            else:
                while True:
                    d = abs(better[i] - worse[i])
                    if (better[i] <= worse[i]):
                        min = better[i] - d * alpha
                        max = worse[i] + d * beta
                    else:
                        min = worse[i] - d * beta
                        max = better[i] + d * alpha
                    tm1, tm2 = min + random.random() * (max - min), min + random.random() * (max - min)
                    # check new values with user-defined bounds
                    if self.bit_limits[i][0] <= tm1 <= self.bit_limits[i][1] and self.bit_limits[i][0] <= tm2 <= self.bit_limits[i][1]: break
                better[i], worse[i] = tm1, tm2
                if self.chromosome_type[i] == 'int':
                    better[i], worse[i] = int(better[i]), int(worse[i])
        return tuple(deepcopy(better)), tuple(deepcopy(worse))

    def fused(self, ind1, ind2, fitness_dict):
        ind1, ind2 = list(ind1), list(ind2)
        x_ind1, y_ind1 = ind1[:self.fused_cutoff], ind1[self.fused_cutoff:]
        x_ind2, y_ind2 = ind2[:self.fused_cutoff], ind2[self.fused_cutoff:]
        x_ind1, x_ind2 = self.blend(x_ind1, x_ind2, fitness_dict)
        y_ind1, y_ind2 = self.UniformCrossover(y_ind1, y_ind2)
        return tuple(deepcopy(list(x_ind1) + list(y_ind1))), tuple(deepcopy(list(x_ind2) + list(y_ind2)))

    def select(self, population, fit_dict, num, choice="Roulette"):
        if num >= len(population): return population
        o_fits = [fit_dict[i] for i in population]

        df_fits = pd.DataFrame(o_fits)
        # scale all values in range 1-2
        df2 = [((df_fits[i] - df_fits[i].min()) / (df_fits[i].max() - df_fits[i].min())) + 1 for i in range(df_fits.shape[1])]
        # inverse min columns
        df2 = pd.DataFrame([df2[i]**self.fit_val[i] for i in range(len(df2))]).T
        # rescale all values in range 1-2
        df2 = pd.DataFrame([((df2[i] - df2[i].min()) / (df2[i].max() - df2[i].min())) + 1 for i in range(df2.shape[1])])

        fitnesses = list(df2.sum())

        if choice == "Roulette":
            total_fitness = float(sum(fitnesses))
            rel_fitness = [f/total_fitness for f in fitnesses]
            # Generate probability intervals for each individual
            probs = [sum(rel_fitness[:i+1]) for i in range(len(rel_fitness))]
            # Draw new population
            new_population = []
            for _ in range(num):
                r = random.random()
                for i, individual in enumerate(population):
                    if r <= probs[i]:
                        new_population.append(deepcopy(individual))
                        break
            return new_population
        else:
            fits_sort = sorted(fitnesses, reverse=True)
            best = [deepcopy(population[fitnesses.index(fits_sort[i])]) for i in range(min(num, len(population)))]
            return best

    def custom_mutate(self, indi, fitness_dict):
        # rank all individuals in fitness dict
        ranked_list = self.select([i for i in fitness_dict], fitness_dict, len(fitness_dict.items()), choice="best")
        # calculate parameter to adjust Gaussian distribution according to individual's rank for uniform type
        parent_fit_param = (ranked_list.index(indi) + 1)*2/len(ranked_list)
        indi = list(indi)
        # if there is no uniform type hyperparamter in the space variable, run the 'if' condition below to select from a pre-defined superlist of mutations and crossovers.
        if self.global_cm_list is not None:
            try:
                if self.chromosome_length == 1:
                    indi = random.choice(list(set(self.global_cm_list)-set([i[0] for i in fitness_dict.keys()])))
                    if isinstance(indi, int): return indi,
                    else: return tuple(indi)
                else:
                    indi = random.choice(list(set(self.global_cm_list)-set(list(fitness_dict.keys()))))
                    return tuple(indi)
            except:
                return None
        for i in range(self.chromosome_length):
            if self.chromosome_type[i] == 'uniform':
                if random.random() < self.mutation_prob:
                    while True:
                        # modify the gaussian mean and standard deviation according to individual's rank
                        add = random.gauss(self.mutation_params[i][0], parent_fit_param*self.mutation_params[i][1]) + indi[i]
                        # check validity of new value in the user-defined range
                        if self.bit_limits[i][0] <= add <= self.bit_limits[i][1]: break
                    indi[i] = add
            elif self.chromosome_type[i] == 'int':
                if random.random() < self.mutation_prob:
                    indi[i] = random.randint(self.bit_limits[i][0],
                                            self.bit_limits[i][1])
            elif self.chromosome_type[i] == 'choice':
                if random.random() < self.mutation_prob:
                    indi[i] = random.choice(list(set(self.bit_limits[i]) - set([indi[i]])))
        if tuple(indi) in fitness_dict.keys(): indi = self.custom_mutate(tuple(indi), fitness_dict)
        return tuple(indi)

    def search(self, n_generations=20, early_stopping=10, init_ratio = 0.35, crossover_ratio = 0.35):
        """
        Parameters
        ----------
        n_generations : integer, optional (default = 20)
                An integer for the number of generations to evolve the population for.

        early_stopping : int, optional (default=10)
                Integer specifying the maximum number of generations for which the algorithm can select the same best individual, after which
                the search terminates.

        init_ratio : float, optional (default = 0.4)
            Fraction of initial population to select for next generation. Required only for algorithm 3.

        crossover_ratio : float, optional (default = 0.3)
            Fraction of crossover population to select for next generation. Required only for algorithm 3.

        Attributes
        ----------
        population : list,
            list of individuals from the final generation

        fitness_dict : dict,
            dictionary of all individuals evaluated by the algorithm

        Returns
        -------
        best_ind_df :  pandas dataframe
            A pandas dataframe of best individuals of each generation

        best_ind :  dict,
            The best individual after the last generation.

        """
        def fit_eval(invalid_ind, fitness_dict):
            if invalid_ind:
                invalid_ind = [i for i in invalid_ind if i not in fitness_dict.keys()]
                fitnesses = list(map(self.evaluate, invalid_ind))
                for ind, fit in zip(invalid_ind, fitnesses):
                    fitness_dict[tuple(ind)] = fit
            return fitness_dict

        if init_ratio >=1 or crossover_ratio >=1 or (init_ratio+crossover_ratio)>=1: raise Exception("Sum of parameters init_ratio and crossover_ratio should be in the range (0,1)")
        if self.population is not None:
            pop = self.population
            fitness_dict = self.fitness_dict
        else:
            pop = self.pop_generator(n=self.pop_size)       # list of tuples
            fitness_dict = {}

        # Evaluate the initial population
        fitness_dict = fit_eval(pop, fitness_dict)

        best_indi_per_gen, best_indi_fitness_values, timer, total_pop, convergence, flag = [], [], [], [], 0, False
        for c_gen in range(n_generations):
            if convergence >= early_stopping:
                print("The search converged with convergence criteria = ", early_stopping)
                break
            else:
                st_time = time.time()
                cross_pop, mutant_pop, co_pop, psum = [], [], [], len(list(fitness_dict.items()))

                # Generate crossover population
                co_pop = self.select(pop, fitness_dict, int(math.ceil(self.crossover_size)))
                co_pop = list(itertools.combinations(list(set(co_pop)), 2))
                combi = list(itertools.combinations(list(set(pop + total_pop)), 2))
                co_pop += combi
                for child1, child2 in co_pop:
                    if (len(list(fitness_dict.items())) - psum) >= int(math.ceil(self.crossover_size)): break
                    if self.crossover_type == "SinglePoint":
                        c1, c2 = self.SinglePointCrossover(child1, child2)
                    elif self.crossover_type == "DoublePoint":
                        c1, c2 = self.DoublePointCrossover(child1, child2)
                    elif self.crossover_type == "Blend":
                        c1, c2 = self.blend(child1, child2, fitness_dict)
                    elif self.crossover_type == "Fused":
                        c1, c2 = self.fused(child1, child2, fitness_dict)
                    elif self.crossover_type == "Uniform":
                        c1, c2 = self.UniformCrossover(child1, child2)
                    if c1 in fitness_dict.keys() or c2 in fitness_dict.keys() or c1==c2: continue
                    fitness_dict = fit_eval([c1, c2], fitness_dict)
                    cross_pop.extend([c1, c2])

                # Generate mutation population
                if self.algo == 4:
                    mu_pop = self.select(cross_pop, fitness_dict, int(math.ceil(self.mutation_size)))
                else:
                    mu_pop = self.select(pop, fitness_dict, int(math.ceil(self.mutation_size)))

                for mutant in mu_pop:
                    a = self.custom_mutate(mutant, fitness_dict)
                    if a is not None:
                        mutant_pop.append(a)
                        fitness_dict = fit_eval([a], fitness_dict)
                    else:
                        print("All combinations exhausted. Stopping genetic algorithm iterations.")
                        flag = True
                        break

                # Select the next generation individuals
                total_pop = pop + cross_pop + mutant_pop
                if self.algo == 2 and c_gen != n_generations - 1:
                    pop = self.select(total_pop, fitness_dict, self.pop_size)
                elif self.algo == 3:
                    p1 = self.select(pop, fitness_dict, int(init_ratio*self.pop_size), choice="best")
                    p2 = self.select(cross_pop, fitness_dict, int(crossover_ratio*self.pop_size), choice="best")
                    p3 = self.select(mutant_pop, fitness_dict, self.crossover_size+self.mutation_size-len(p1)-len(p2), choice="best")
                    pop = p1 + p2 + p3
                else: pop = self.select(total_pop, fitness_dict, self.pop_size, choice="best")

                # Storing the best individuals after each generation
                best_individual = pop[0]
                if len(best_indi_per_gen)>0:
                    if best_individual==best_indi_per_gen[-1]: convergence += 1
                    else: convergence = 0
                best_indi_per_gen.append(best_individual)
                best_indi_fitness_values.append(fitness_dict[best_individual])
                tot_time = (time.time() - st_time)/(60*60)
                timer.append(tot_time)
                b1 = pd.Series(best_indi_per_gen, name='Best_individual')
                b2 = pd.Series(best_indi_fitness_values, name='Fitness_values')
                b3 = pd.Series(timer, name='Time (hours)')
                best_ind_df = pd.concat([b1, b2, b3], axis=1)
                if flag: break

        self.population = pop    # stores best individuals of last generation
        self.fitness_dict = fitness_dict
        best_ind_dict = {}
        for name, val in zip(self.var_names, best_individual):
            best_ind_dict[name] = val
        return best_ind_df, best_ind_dict

max_issues_repo_name: hachmannlab/ChemML

max_forks_repo_name: hachmannlab/ChemML


text:
import argparse
import os
from tqdm import tqdm
from solver import Solver
from data_loader_2 import get_loader
from torch.backends import cudnn
import random
import torch
import cv2
from network import U_Net,R2U_Net,AttU_Net,R2AttU_Net
from torch import optim
import torch.nn.functional as F
import matplotlib.pyplot as plt
import numpy as np
def compute_integrated_gradient(batch_x, model,px,py):
    mean_grad = 0
    n = 100
    batch_blank=torch.zeros_like(batch_x)
    for i in tqdm(range(1, n + 1)):
        x = batch_blank + i / n * (batch_x - batch_blank)
        x.requires_grad = True
        y = model(x)[0,0,px,py]
        (grad,) = torch.autograd.grad(y, x)
        mean_grad += grad / n

    integrated_gradients = (batch_x - batch_blank) * mean_grad
    return integrated_gradients

def plot_images(images,output,gdmap,output_path):
    fig, axs = plt.subplots(1,3)
    fig.set_figheight(10)
    fig.set_figwidth(16)
    axs[0].imshow(images)
    axs[2].imshow(output,cmap = plt.cm.gray)
    axs[1].imshow(gdmap)
    fig.tight_layout()
    plt.savefig(output_path)

def plot_grad(data):
    abs_max=max(data.max(),-data.min())
    img=np.ones((256,256,3))*255.0
    img[:,:,0]+=data*(data<0)/abs_max*255
    img[:,:,1]+=data*(data<0)/abs_max*255
    img[:,:,1]-=data*(data>0)/abs_max*255
    img[:,:,2]-=data*(data>0)/abs_max*255
    return img.astype(np.uint8)

def options():
    parser = argparse.ArgumentParser()

    # model hyper-parameters
    parser.add_argument('--image_size', type=int, default=1024)
    parser.add_argument('--t', type=int, default=3, help='t for Recurrent step of R2U_Net or R2AttU_Net')

    # training hyper-parameters
    parser.add_argument('--img_ch', type=int, default=3)
    parser.add_argument('--output_ch', type=int, default=1)
    parser.add_argument('--num_epochs', type=int, default=50)
    parser.add_argument('--num_epochs_decay', type=int, default=70)
    parser.add_argument('--batch_size', type=int, default=1)
    parser.add_argument('--num_workers', type=int, default=8)
    parser.add_argument('--lr', type=float, default=0.0002)
    parser.add_argument('--beta1', type=float, default=0.5)  # momentum1 in Adam
    parser.add_argument('--beta2', type=float, default=0.999)  # momentum2 in Adam
    parser.add_argument('--augmentation_prob', type=float, default=0.4)

    parser.add_argument('--log_step', type=int, default=2)
    parser.add_argument('--val_step', type=int, default=2)

    # misc
    parser.add_argument('--model_type', type=str, default='U_Net', help='U_Net/R2U_Net/AttU_Net/R2AttU_Net')
    parser.add_argument('--model_name', type=str, default='U_Net-149.pkl', help='U_Net-149.pkl')
    parser.add_argument('--prediction_path', type=str, default='/raid/zsj/srh/IS/prediction/')
    parser.add_argument('--model_path', type=str, default='/home/srh/Image_Segmentation/models')
    parser.add_argument('--test_path', type=str, default='/home/srh/PyTorch-U-Net/ISBI2012/imgs/test/')
    parser.add_argument('--result_path', type=str, default='/raid/srh/IS/result/')

    config = parser.parse_args()

    return config

class Solver(object):
    def __init__(self, config, test_loader):

        # Data loader
        self.test_loader = test_loader

        # Models
        self.unet = None
        self.optimizer = None
        self.img_ch = config.img_ch
        self.output_ch = config.output_ch
        self.criterion = torch.nn.BCELoss()
        self.augmentation_prob = config.augmentation_prob

        # Hyper-parameters
        self.lr = config.lr
        self.beta1 = config.beta1
        self.beta2 = config.beta2

        # Training settings
        self.num_epochs = config.num_epochs
        self.num_epochs_decay = config.num_epochs_decay
        self.batch_size = config.batch_size

        # Step size
        self.log_step = config.log_step
        self.val_step = config.val_step

        # Path
        self.model_path = config.model_path
        self.model_name = config.model_name
        self.result_path = config.result_path
        self.prediction_path = config.prediction_path
        self.test_path = config.test_path

        self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
        self.model_type = config.model_type
        self.t = config.t
        self.build_model()

    def build_model(self):
        """Build generator and discriminator."""
        if self.model_type =='U_Net':
            self.unet = U_Net(img_ch=3,output_ch=1)
        elif self.model_type =='R2U_Net':
            self.unet = R2U_Net(img_ch=3,output_ch=1,t=self.t)
        elif self.model_type =='AttU_Net':
            self.unet = AttU_Net(img_ch=3,output_ch=1)
        elif self.model_type == 'R2AttU_Net':
            self.unet = R2AttU_Net(img_ch=3,output_ch=1,t=self.t)

        self.optimizer = optim.Adam(list(self.unet.parameters()),
                                      self.lr, [self.beta1, self.beta2])
        self.unet.to(self.device)
        self.unet = torch.nn.DataParallel(self.unet)
        self.unet = self.unet.cuda()

        # self.print_network(self.unet, self.model_type)

    def print_network(self, model, name):
        """Print out the network information."""
        num_params = 0
        for p in model.parameters():
            num_params += p.numel()
        print(model)
        print(name)
        print("The number of parameters: {}".format(num_params))

    def to_data(self, x):
        """Convert variable to tensor."""
        if torch.cuda.is_available():
            x = x.cuda()
        return x.data

    def update_lr(self, g_lr, d_lr):
        for param_group in self.optimizer.param_groups:
            param_group['lr'] = d_lr

    def reset_grad(self):
        """Zero the gradient buffers."""
        self.unet.zero_grad()

    def tensor2img(self,x):
        img = (x[:,0,:,:]>x[:,1,:,:]).float()
        img = img*255
        return img

    def test(self):
        unet_path = os.path.join(self.model_path, self.model_name)
        save_path = self.prediction_path + self.model_type + '/' + self.model_name.split('.')[0] +'/'
        if not os.path.exists(save_path):
            os.makedirs(save_path)

        self.build_model()
        print("build_success")
        self.unet.load_state_dict(torch.load(unet_path))

        self.unet.train(False)
        self.unet.eval()

        for i, (images, GT, filename) in enumerate(self.test_loader):
            images = images.to(self.device)
            GT = GT.to(self.device)
            SR = F.sigmoid(self.unet(images))
            out = SR.cpu().detach().numpy()
            print(filename[0])

            # draw by pixel
            for j in range(0,256,8):
                for k in range(0,256,8):
                    ig=compute_integrated_gradient(images,self.unet,j,k)
                    ig=plot_grad(ig[0,0].cpu().detach().numpy())
                    pic=np.zeros((256,256,3)).astype(np.uint8)
                    pic[:,:,0]=out[0,0,:,:]*255
                    pic[:,:,1]=out[0,0,:,:]*255
                    pic[:,:,2]=out[0,0,:,:]*255
                    for u in range(-2,2):
                        for v in range(-2,2):
                            if j+u>=0 and j+u<256 and k+v>=0 and k+v<256:
                                pic[j+u,k+v,0]=255
                                pic[j+u,k+v,1]=0
                                pic[j+u,k+v,2]=0
                    input_img = cv2.imread(self.test_path+filename[0])
                    # plot_images(images[0].cpu().permute(1,2,0),pic,ig,'/raid/zsj/srh/IS/IGtest/'+str(j)+'_'+str(k)+'.png')
                    plot_images(input_img,pic,ig,'/raid/srh/IS/IGtest_2012/'+str(j)+'_'+str(k)+'.png')

            # # draw by TP
            # TP = np.zeros((256,256)).astype(np.uint8)+255
            # new_out, new_gt = np.zeros((256,256)).astype(np.uint8),np.zeros((256,256)).astype(np.uint8)
            # new_out[out[0,0,:,:]>0.5] = 1
            # new_gt[GT[0,0,:,:].cpu().detach().numpy()>0.5] = 1
            # all_nodes = []
            # for m in range(new_gt.shape[0]):
            #     for n in range(new_gt.shape[1]):
            #         if new_gt[m][n] == 0 and new_out[m][n] == 0:
            #             TP[m][n] == 0
            #             all_nodes.append((m,n))
            # print('TP nodes num is: ', len(all_nodes))

            # all_ig = np.zeros((256,256)).astype(np.uint8)

            # # for j in range(len(all_nodes)):
            # for j in range(3):
            #     print(all_nodes[j][0],all_nodes[j][1])
            #     ig_i=compute_integrated_gradient(images,self.unet,all_nodes[j][0],all_nodes[j][1])
            #     all_ig += ig_i[0,0].cpu().detach().numpy().astype(np.uint8)

            # np.save('all_ig.npy', all_ig)
            # print(all_ig.max())
            # ig=plot_grad(all_ig)

            # pic=np.zeros((256,256,3)).astype(np.uint8)
            # pic[:,:,0]=TP
            # pic[:,:,1]=TP
            # pic[:,:,2]=TP
            # # for u in range(-2,2):
            # #     for v in range(-2,2):
            # #         if j+u>=0 and j+u<256 and k+v>=0 and k+v<256:
            # #             pic[j+u,k+v,0]=255
            # #             pic[j+u,k+v,1]=0
            # #             pic[j+u,k+v,2]=0

            # input_img = cv2.imread(self.test_path+filename[0])
            # # save_ = '/raid/zsj/srh/IS/'+filename[0].split('/')[0]+'_'+filename[0].split('/')[1]+'.png'
            # save_ = '/raid/srh/IS/'+filename[0].split('/')[0]+'_'+filename[0].split('/')[1]+'.png'
            # # plot_images(images[0].cpu().permute(1,2,0),pic,ig,'/raid/zsj/srh/IS/IGtest/'+str(j)+'_'+str(k)+'.png')
            # plot_images(input_img,pic,ig,save_)

            break

if __name__ == '__main__':
    config = options()
    cudnn.benchmark = True

    lr = random.random() * 0.0005 + 0.0000005
    augmentation_prob = random.random() * 0.7
    epoch = random.choice([100, 150, 200, 250])
    decay_ratio = random.random() * 0.8
    decay_epoch = int(epoch * decay_ratio)

    config.augmentation_prob = augmentation_prob
    config.num_epochs = epoch
    config.lr = lr
    config.num_epochs_decay = decay_epoch

    #==========  load test data  =====#
    test_loader = get_loader(image_path=config.test_path,
                             image_size=config.image_size,
                             batch_size=config.batch_size,
                             num_workers=config.num_workers,
                             mode='test',
                             augmentation_prob=0.)

    solver = Solver(config, test_loader)
    print("start_test")
    solver.test()

ext: py

max_stars_repo_name: EmmaSRH/U-RISC-Data-Code

max_issues_repo_name: EmmaSRH/U-RISC-Data-Code

max_forks_repo_name: EmmaSRH/U-RISC-Data-Code

reason: import numpy


text:
#!/usr/bin/env python

"""
m2g.register
~~~~~~~~~~~~~
Contains m2g's registration classes, organized as full registration workflows.
Used for the majority of the registration described here: https://neurodata.io/talks/ndmg.pdf#page=20
"""

# standard library imports
import os
from argparse import ArgumentParser
import subprocess
from pathlib import Path

# package imports
import nibabel as nib
import numpy as np
from nilearn.image import load_img
from nilearn.image import math_img
from dipy.tracking.streamline import deform_streamlines
from dipy.io.streamline import load_trk
from dipy.tracking import utils

# m2g imports
from m2g.utils import gen_utils
from m2g.utils import reg_utils
from m2g.scripts import m2g_bids
from m2g.stats.qa_skullstrip import gen_overlay_pngs
from m2g.stats.qa_reg import reg_mri_pngs
from m2g.stats.qa_fast import qa_fast_png

@gen_utils.timer
def direct_streamline_norm(streams, fa_path, outdir: Path):
    """Applys the Symmetric Diffeomorphic Registration (SyN) Algorithm onto the streamlines to the atlas space defined by .../atlases/reference_brains/FSL_HCP1065_FA_2mm.nii.gz
    Parameters
    ----------
    streams : str
        Path to streamlines.trk file to be transformed
    fa_path : str
        Path to subject's FA tensor image
    outdir : Path
        Output directory location.

    Returns
    -------
    ArraySequence
        Transformed streamlines
    str
        Path to tractogram streamline file: streamlines_dsn.trk
    """

    atlas_dir = m2g_bids.get_atlas_dir()
    template_path = atlas_dir + "/atlases/reference_brains/FSL_HCP1065_FA_2mm.nii.gz"

    # Run SyN and normalize streamlines
    fa_img = nib.load(fa_path)
    vox_size = fa_img.get_header().get_zooms()[0]
    template_img = nib.load(template_path)

    # SyN FA->Template
    mapping, affine_map = reg_utils.wm_syn(template_path, fa_path, str(outdir / "tmp"))
    streamlines = load_trk(streams, reference="same")

    # Warp streamlines
    adjusted_affine = affine_map.affine.copy()
    adjusted_affine[1][3] = -adjusted_affine[1][3] / vox_size ** 2
    mni_streamlines = deform_streamlines(
        streamlines.streamlines,
        deform_field=mapping.get_forward_field()[-1:],
        stream_to_current_grid=template_img.affine,
        current_grid_to_world=adjusted_affine,
        stream_to_ref_grid=template_img.affine,
        ref_grid_to_world=np.eye(4),
    )

    # Save streamlines
    hdr = fa_img.header
    trk_affine = np.eye(4)
    trk_hdr = nib.streamlines.trk.TrkFile.create_empty_header()
    trk_hdr["hdr_size"] = 1000
    trk_hdr["dimensions"] = hdr["dim"][1:4].astype("float32")
    trk_hdr["voxel_sizes"] = hdr["pixdim"][1:4]
    trk_hdr["voxel_to_rasmm"] = trk_affine
    trk_hdr["voxel_order"] = "RAS"
    trk_hdr["pad2"] = "RAS"
    trk_hdr["image_orientation_patient"] = np.array(
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    ).astype("float32")
    trk_hdr["endianness"] = "<"
    trk_hdr["_offset_data"] = 1000
    trk_hdr["nb_streamlines"] = len(mni_streamlines)
    tractogram = nib.streamlines.Tractogram(mni_streamlines, affine_to_rasmm=trk_affine)
    trkfile = nib.streamlines.trk.TrkFile(tractogram, header=trk_hdr)
    streams_mni = streams.split(".trk")[0] + "_dsn.trk"
    nib.streamlines.save(trkfile, streams_mni)
    return mni_streamlines, streams_mni

class DmriReg:
    """Class containing relevant paths and class methods for analysing tractography
    Parameters
    ----------
    outdir : Path
        Output directory location.
    nodif_B0 : str
        path to mean b0 image
    nodif_B0_mask : str
        path to mean b0 mask (nodif_B0....nii.gz)
    t1w_in : str
        path to t1w file
    vox_size : str
        voxel resolution ('2mm' or '1mm')
    skull : str, optional
        skullstrip parameter pre-set. Default is "none"
    simple : bool, optional
        Whether you want to attempt non-linear registration when transforming between mni, t1w, and dwi space. Default is False
    Raises
    ------
    ValueError
        FSL atlas for ventricle reference not found
    """

    def __init__(
        self,
        outdir: Path,
        nodif_B0,
        nodif_B0_mask,
        t1w_in,
        vox_size,
        skull=None,
        simple=False,
    ):

        atlas_dir = m2g_bids.get_atlas_dir()
        FSLDIR = os.environ["FSLDIR"]

        if vox_size == "2mm":
            vox_dims = "2x2x2"
        elif vox_size == "1mm":
            vox_dims = "1x1x1"

        # TODO : clean up all these attributes
        self.outdir = outdir
        self.reg_a = str(outdir / "tmp/reg_a")
        self.reg_m = str(outdir / "tmp/reg_m")
        self.reg_anat = str(outdir / "anat/registered")
        self.prep_anat = str(outdir / "anat/preproc")
        self.qa = str(outdir / "qa")
        self.qa_reg = str(Path(self.qa) / "reg")
        self.simple = simple
        self.nodif_B0 = nodif_B0
        self.nodif_B0_mask = nodif_B0_mask
        self.t1w = t1w_in
        self.vox_size = vox_size
        self.t1w_name = "t1w"
        self.dwi_name = "dwi"
        self.skull = skull
        self.t12mni_xfm_init = f"{self.reg_m}" + "/xfm_t1w2mni_init.mat"
        self.mni2t1_xfm_init = f"{self.reg_m}" + "/xfm_mni2t1w_init.mat"
        self.t12mni_xfm = f"{self.reg_m}" + "/xfm_t1w2mni.mat"
        self.mni2t1_xfm = f"{self.reg_m}" + "/xfm_mni2t1.mat"
        self.mni2t1w_warp = f"{self.reg_a}/mni2t1w_warp.nii.gz"
        self.warp_t1w2mni = f"{self.reg_a}/warp_t12mni.nii.gz"
        self.t1w2dwi = f"{self.reg_anat}/{self.t1w_name}_in_dwi.nii.gz"
        self.t1_aligned_mni = f"{self.prep_anat}/{self.t1w_name}_aligned_mni.nii.gz"
        self.t1w_brain = f"{self.prep_anat}/{self.t1w_name}_brain.nii.gz"
        self.dwi2t1w_xfm = f"{self.reg_m}" + "/dwi2t1w_xfm.mat"
        self.t1w2dwi_xfm = f"{self.reg_m}" + "/t1w2dwi_xfm.mat"
        self.t1w2dwi_bbr_xfm = f"{self.reg_m}" + "/t1w2dwi_bbr_xfm.mat"
        self.dwi2t1w_bbr_xfm = f"{self.reg_m}" + "/dwi2t1w_bbr_xfm.mat"
        self.t1wtissue2dwi_xfm = f"{self.reg_m}" + "/t1wtissue2dwi_xfm.mat"
        self.xfm_atlas2t1w_init = (
            f"{self.reg_m}" + f"/{self.t1w_name}_xfm_atlas2t1w_init.mat"
        )
        self.xfm_atlas2t1w = f"{self.reg_m}/{self.t1w_name}_xfm_atlas2t1w.mat"
        self.temp2dwi_xfm = f"{self.reg_m}/{self.dwi_name}_xfm_temp2dwi.mat"

        self.input_mni = f"{FSLDIR}/data/standard/MNI152_T1_{vox_size}_brain.nii.gz"
        self.input_mni_mask = (
            f"{FSLDIR}/data/standard/MNI152_T1_{vox_size}_brain_mask.nii.gz"
        )
        self.temp2dwi_xfm = f"{self.reg_m}/{self.dwi_name}_xfm_temp2dwi.mat"
        self.map_path = f"{self.prep_anat}/{self.t1w_name}_seg"
        self.wm_mask = f"{self.prep_anat}/{self.t1w_name}_wm.nii.gz"
        self.wm_mask_thr = f"{self.prep_anat}/{self.t1w_name}_wm_thr.nii.gz"
        self.wm_edge = f"{self.reg_a}/{self.t1w_name}_wm_edge.nii.gz"
        self.csf_mask = f"{self.prep_anat}/{self.t1w_name}_csf.nii.gz"
        self.gm_mask = f"{self.prep_anat}/{self.t1w_name}_gm.nii.gz"
        self.xfm_roi2mni_init = f"{self.reg_m}/roi_2_mni.mat"
        self.lvent_out_file = f"{self.reg_a}/LVentricle.nii.gz"
        self.rvent_out_file = f"{self.reg_a}/RVentricle.nii.gz"
        self.csf_mask_dwi = f"{self.reg_anat}/{self.t1w_name}_csf_mask_dwi.nii.gz"
        self.gm_in_dwi = f"{self.reg_anat}/{self.t1w_name}_gm_in_dwi.nii.gz"
        self.wm_in_dwi = f"{self.reg_anat}/{self.t1w_name}_wm_in_dwi.nii.gz"
        self.csf_mask_dwi_bin = f"{self.reg_a}/{self.t1w_name}_csf_mask_dwi_bin.nii.gz"
        self.gm_in_dwi_bin = f"{self.reg_a}/{self.t1w_name}_gm_in_dwi_bin.nii.gz"
        self.wm_in_dwi_bin = f"{self.reg_a}/{self.t1w_name}_wm_in_dwi_bin.nii.gz"
        self.vent_mask_dwi = f"{self.reg_a}/{self.t1w_name}_vent_mask_dwi.nii.gz"
        self.vent_csf_in_dwi = f"{self.reg_a}/{self.t1w_name}_vent_csf_in_dwi.nii.gz"
        self.vent_mask_mni = f"{self.reg_a}/vent_mask_mni.nii.gz"
        self.vent_mask_t1w = f"{self.reg_a}/vent_mask_t1w.nii.gz"
        self.wm_gm_int_in_dwi = (
            f"{self.reg_anat}/{self.t1w_name}_wm_gm_int_in_dwi.nii.gz"
        )
        self.wm_gm_int_in_dwi_bin = (
            f"{self.reg_anat}/{self.t1w_name}_wm_gm_int_in_dwi_bin.nii.gz"
        )
        self.input_mni_sched = f"{FSLDIR}/etc/flirtsch/T1_2_MNI152_2mm.cnf"
        self.mni_atlas = f"{atlas_dir}/atlases/label/Human/HarvardOxfordsub-maxprob-thr25_space-MNI152NLin6_label_all_res-{vox_dims}.nii.gz"
        self.mni_vent_loc = f"{atlas_dir}/atlases/mask/HarvardOxford-thr25_space-MNI152NLin6_variant-lateral-ventricles_res-{vox_dims}_descr-brainmask.nii.gz"
        self.corpuscallosum = (
            f"{atlas_dir}/atlases/mask/CorpusCallosum_res_{vox_size}.nii.gz"
        )
        self.corpuscallosum_mask_t1w = (
            f"{self.reg_anat}/{self.t1w_name}_corpuscallosum.nii.gz"
        )
        self.corpuscallosum_dwi = (
            f"{self.reg_anat}/{self.t1w_name}_corpuscallosum_dwi.nii.gz"
        )

    @gen_utils.timer
    def gen_tissue(self):
        """Extracts the brain from the raw t1w image (as indicated by self.t1w), uses it to create WM, GM, and CSF masks,
        reslices all 4 files to the target voxel resolution and extracts the white matter edge. Each mask is saved to
        location indicated by self.map_path
        """
        # BET needed for this, as afni 3dautomask only works on 4d volumes
        print("Extracting brain from raw T1w image...")
        reg_utils.t1w_skullstrip(self.t1w, self.t1w_brain, self.skull)

        #  QA part of skull strip
        skullstrip_qa = Path(self.qa) / "skull_strip"
        if not os.path.exists(skullstrip_qa):
            skullstrip_qa.mkdir(parents=True, exist_ok=True)
        print('QA_skullstrip_path  ', skullstrip_qa)
        gen_overlay_pngs(self.t1w_brain, self.t1w, skullstrip_qa)

        # Segment the t1w brain into probability maps
        self.maps = reg_utils.segment_t1w(self.t1w_brain, self.map_path)
        self.wm_mask = self.maps["wm_prob"]
        self.gm_mask = self.maps["gm_prob"]
        self.csf_mask = self.maps["csf_prob"]

        # Generates quality analysis pictures of white matter, gray matter and cerebrospinal fluid
        qa_fast_png(self.csf_mask, self.gm_mask, self.wm_mask, str(Path(self.qa_reg)/"qa_fast.png"))

        self.t1w_brain = gen_utils.match_target_vox_res(
            self.t1w_brain, self.vox_size, self.outdir, sens="anat"
        )
        self.wm_mask = gen_utils.match_target_vox_res(
            self.wm_mask, self.vox_size, self.outdir, sens="anat"
        )
        self.gm_mask = gen_utils.match_target_vox_res(
            self.gm_mask, self.vox_size, self.outdir, sens="anat"
        )
        self.csf_mask = gen_utils.match_target_vox_res(
            self.csf_mask, self.vox_size, self.outdir, sens="anat"
        )

        # Threshold WM to binary in dwi space
        self.t_img = load_img(self.wm_mask)
        self.mask = math_img("img > 0.2", img=self.t_img)
        self.mask.to_filename(self.wm_mask_thr)

        # Extract wm edge
        # TODO : this should be a function in reg_utils so that we can print it to the log
        cmd = f"fslmaths {self.wm_mask_thr} -edge -bin -mas {self.wm_mask_thr} {self.wm_edge}"
        print("Extracting white matter edge...")
        gen_utils.run(cmd)

    @gen_utils.timer
    def t1w2dwi_align(self):
        """Alignment from t1w to mni, making t1w_mni, and t1w_mni to dwi. A function to perform self alignment. Uses a local optimisation cost function to get the
        two images close, and then uses bbr to obtain a good alignment of brain boundaries. Assumes input dwi is already preprocessed and brain extracted.
        """

        # Create linear transform/ initializer T1w-->MNI
        reg_utils.align(
            self.t1w_brain,
            self.input_mni,
            xfm=self.t12mni_xfm_init,
            bins=None,
            interp="spline",
            out=None,
            dof=12,
            cost="mutualinfo",
            searchrad=True,
        )

        # Attempt non-linear registration of T1 to MNI template
        if self.simple is False:
            try:
                print("Running non-linear registration: T1w-->MNI ...")
                # Use FNIRT to nonlinearly align T1 to MNI template
                reg_utils.align_nonlinear(
                    self.t1w_brain,
                    self.input_mni,
                    xfm=self.t12mni_xfm_init,
                    out=self.t1_aligned_mni,
                    warp=self.warp_t1w2mni,
                    ref_mask=self.input_mni_mask,
                    config=self.input_mni_sched,
                )

                # Get warp from MNI -> T1
                reg_utils.inverse_warp(
                    self.t1w_brain, self.mni2t1w_warp, self.warp_t1w2mni
                )

                # Get mat from MNI -> T1
                cmd = f"convert_xfm -omat {self.mni2t1_xfm_init} -inverse {self.t12mni_xfm_init}"
                print(cmd)
                gen_utils.run(cmd)

            except RuntimeError("Error: FNIRT failed!"):
                pass
        else:
            # Falling back to linear registration
            reg_utils.align(
                self.t1w_brain,
                self.input_mni,
                xfm=self.t12mni_xfm,
                init=self.t12mni_xfm_init,
                bins=None,
                dof=12,
                cost="mutualinfo",
                searchrad=True,
                interp="spline",
                out=self.t1_aligned_mni,
                sch=None,
            )
        reg_mri_pngs(self.t1_aligned_mni, self.input_mni, self.qa_reg)

        # Align T1w-->DWI
        reg_utils.align(
            self.nodif_B0,
            self.t1w_brain,
            xfm=self.t1w2dwi_xfm,
            bins=None,
            interp="spline",
            dof=6,
            cost="mutualinfo",
            out=None,
            searchrad=True,
            sch=None,
        )
        cmd = f"convert_xfm -omat {self.dwi2t1w_xfm} -inverse {self.t1w2dwi_xfm}"
        print(cmd)
        gen_utils.run(cmd)

        if self.simple is False:
            # Flirt bbr
            try:
                print("Running FLIRT BBR registration: T1w-->DWI ...")
                reg_utils.align(
                    self.nodif_B0,
                    self.t1w_brain,
                    wmseg=self.wm_edge,
                    xfm=self.dwi2t1w_bbr_xfm,
                    init=self.dwi2t1w_xfm,
                    bins=256,
                    dof=7,
                    searchrad=True,
                    interp="spline",
                    out=None,
                    cost="bbr",
                    finesearch=5,
                    sch="${FSLDIR}/etc/flirtsch/bbr.sch",
                )
                cmd = f"convert_xfm -omat {self.t1w2dwi_bbr_xfm} -inverse {self.dwi2t1w_bbr_xfm}"
                gen_utils.run(cmd)

                # Apply the alignment
                reg_utils.align(
                    self.t1w_brain,
                    self.nodif_B0,
                    init=self.t1w2dwi_bbr_xfm,
                    xfm=self.t1wtissue2dwi_xfm,
                    bins=None,
                    interp="spline",
                    dof=7,
                    cost="mutualinfo",
                    out=self.t1w2dwi,
                    searchrad=True,
                    sch=None,
                )
            except RuntimeError("Error: FLIRT BBR failed!"):
                pass
        else:
            # Apply the alignment
            reg_utils.align(
                self.t1w_brain,
                self.nodif_B0,
                init=self.t1w2dwi_xfm,
                xfm=self.t1wtissue2dwi_xfm,
                bins=None,
                interp="spline",
                dof=6,
                cost="mutualinfo",
                out=self.t1w2dwi,
                searchrad=True,
                sch=None,
            )
        reg_mri_pngs(self.t1w2dwi, self.nodif_B0, self.qa_reg)

    def atlas2t1w2dwi_align(self, atlas, dsn=True):
        """alignment from atlas to t1w to dwi. A function to perform atlas alignmet. Tries nonlinear registration first, and if that fails, does a liner
        registration instead.
        Note: for this to work, must first have called t1w2dwi_align.
        Parameters
        ----------
        atlas : str
            path to atlas file you want to use
        dsn : bool, optional
            is your space for tractography native-dsn, by default True
        Returns
        -------
        str
            path to aligned atlas file
        """

        self.atlas = atlas
        self.atlas_name = gen_utils.get_filename(self.atlas)
        self.aligned_atlas_t1mni = (
            f"{self.reg_a}/{self.atlas_name}_aligned_atlas_t1w_mni.nii.gz"
        )
        self.aligned_atlas_skull = (
            f"{self.reg_a}/{self.atlas_name}_aligned_atlas_skull.nii.gz"
        )
        self.dwi_aligned_atlas = (
            f"{self.reg_anat}/{self.atlas_name}_aligned_atlas.nii.gz"
        )

        reg_utils.align(
            self.atlas,
            self.t1_aligned_mni,
            init=None,
            xfm=None,
            out=self.aligned_atlas_t1mni,
            dof=12,
            searchrad=True,
            interp="nearestneighbour",
            cost="mutualinfo",
        )

        if (self.simple is False) and (dsn is False):
            try:
                # Apply warp resulting from the inverse of T1w-->MNI created earlier
                reg_utils.apply_warp(
                    self.t1w_brain,
                    self.aligned_atlas_t1mni,
                    self.aligned_atlas_skull,
                    warp=self.mni2t1w_warp,
                    interp="nn",
                    sup=True,
                )

                # Apply transform to dwi space
                reg_utils.align(
                    self.aligned_atlas_skull,
                    self.nodif_B0,
                    init=self.t1wtissue2dwi_xfm,
                    xfm=None,
                    out=self.dwi_aligned_atlas,
                    dof=6,
                    searchrad=True,
                    interp="nearestneighbour",
                    cost="mutualinfo",
                )
            except:
                print(
                    "Warning: Atlas is not in correct dimensions, or input is low quality,\nusing linear template registration."
                )
                # Create transform to align atlas to T1w using flirt
                reg_utils.align(
                    self.atlas,
                    self.t1w_brain,
                    xfm=self.xfm_atlas2t1w_init,
                    init=None,
                    bins=None,
                    dof=6,
                    cost="mutualinfo",
                    searchrad=True,
                    interp="spline",
                    out=None,
                    sch=None,
                )
                reg_utils.align(
                    self.atlas,
                    self.t1_aligned_mni,
                    xfm=self.xfm_atlas2t1w,
                    out=None,
                    dof=6,
                    searchrad=True,
                    bins=None,
                    interp="spline",
                    cost="mutualinfo",
                    init=self.xfm_atlas2t1w_init,
                )

                # Combine our linear transform from t1w to template with our transform from dwi to t1w space to get a transform from atlas ->(-> t1w ->)-> dwi
                reg_utils.combine_xfms(
                    self.xfm_atlas2t1w, self.t1wtissue2dwi_xfm, self.temp2dwi_xfm
                )

                # Apply linear transformation from template to dwi space
                reg_utils.applyxfm(
                    self.nodif_B0, self.atlas, self.temp2dwi_xfm, self.dwi_aligned_atlas
                )
        elif dsn is False:
            # Create transform to align atlas to T1w using flirt
            reg_utils.align(
                self.atlas,
                self.t1w_brain,
                xfm=self.xfm_atlas2t1w_init,
                init=None,
                bins=None,
                dof=6,
                cost="mutualinfo",
                searchrad=None,
                interp="spline",
                out=None,
                sch=None,
            )
            reg_utils.align(
                self.atlas,
                self.t1w_brain,
                xfm=self.xfm_atlas2t1w,
                out=None,
                dof=6,
                searchrad=True,
                bins=None,
                interp="spline",
                cost="mutualinfo",
                init=self.xfm_atlas2t1w_init,
            )

            # Combine our linear transform from t1w to template with our transform from dwi to t1w space to get a transform from atlas ->(-> t1w ->)-> dwi
            reg_utils.combine_xfms(
                self.xfm_atlas2t1w, self.t1wtissue2dwi_xfm, self.temp2dwi_xfm
            )

            # Apply linear transformation from template to dwi space
            reg_utils.applyxfm(
                self.nodif_B0, self.atlas, self.temp2dwi_xfm, self.dwi_aligned_atlas
            )
        else:
            pass

        # Set intensities to int
        if dsn is False:
            self.atlas_img = nib.load(self.dwi_aligned_atlas)
        else:
            self.atlas_img = nib.load(self.aligned_atlas_t1mni)
        self.atlas_data = np.around(self.atlas_img.get_data()).astype("int16")
        node_num = len(np.unique(self.atlas_data))
        self.atlas_data[self.atlas_data > node_num] = 0

        t_img = load_img(self.wm_gm_int_in_dwi)
        mask = math_img("img > 0", img=t_img)
        mask.to_filename(self.wm_gm_int_in_dwi_bin)

        if dsn is False:
            nib.save(
                nib.Nifti1Image(
                    self.atlas_data.astype(np.int32),
                    affine=self.atlas_img.affine,
                    header=self.atlas_img.header,
                ),
                self.dwi_aligned_atlas,
            )
            reg_mri_pngs(self.dwi_aligned_atlas, self.nodif_B0, self.qa_reg)
            return self.dwi_aligned_atlas
        else:
            nib.save(
                nib.Nifti1Image(
                    self.atlas_data.astype(np.int32),
                    affine=self.atlas_img.affine,
                    header=self.atlas_img.header,
                ),
                self.aligned_atlas_t1mni,
            )
            reg_mri_pngs(self.aligned_atlas_t1mni, self.t1_aligned_mni, self.qa_reg)
            return self.aligned_atlas_t1mni

    @gen_utils.timer
    def tissue2dwi_align(self):
        """alignment of ventricle and CC ROI's from MNI space --> dwi and CC and CSF from T1w space --> dwi
        A function to generate and perform dwi space alignment of avoidance/waypoint masks for tractography.
        First creates ventricle and CC ROI. Then creates transforms from stock MNI template to dwi space.
        NOTE: for this to work, must first have called both t1w2dwi_align and atlas2t1w2dwi_align.
        Raises
        ------
        ValueError
            Raised if FSL atlas for ventricle reference not found
        """

        # Create MNI-space ventricle mask
        print("Creating MNI-space ventricle ROI...")
        if not os.path.isfile(self.mni_atlas):
            raise ValueError("FSL atlas for ventricle reference not found!")
        cmd = f"fslmaths {self.mni_vent_loc} -thr 0.1 -bin {self.mni_vent_loc}"
        gen_utils.run(cmd)

        cmd = f"fslmaths {self.corpuscallosum} -bin {self.corpuscallosum}"
        gen_utils.run(cmd)

        cmd = f"fslmaths {self.corpuscallosum} -sub {self.mni_vent_loc} -bin {self.corpuscallosum}"
        gen_utils.run(cmd)

        # Create a transform from the atlas onto T1w. This will be used to transform the ventricles to dwi space.
        reg_utils.align(
            self.mni_atlas,
            self.input_mni,
            xfm=self.xfm_roi2mni_init,
            init=None,
            bins=None,
            dof=6,
            cost="mutualinfo",
            searchrad=True,
            interp="spline",
            out=None,
        )

        # Create transform to align roi to mni and T1w using flirt
        reg_utils.applyxfm(
            self.input_mni, self.mni_vent_loc, self.xfm_roi2mni_init, self.vent_mask_mni
        )

        if self.simple is False:
            # Apply warp resulting from the inverse MNI->T1w created earlier
            reg_utils.apply_warp(
                self.t1w_brain,
                self.vent_mask_mni,
                self.vent_mask_t1w,
                warp=self.mni2t1w_warp,
                interp="nn",
                sup=True,
            )

            # Apply warp resulting from the inverse MNI->T1w created earlier
            reg_utils.apply_warp(
                self.t1w_brain,
                self.corpuscallosum,
                self.corpuscallosum_mask_t1w,
                warp=self.mni2t1w_warp,
                interp="nn",
                sup=True,
            )

        # Applyxfm tissue maps to dwi space
        reg_utils.applyxfm(
            self.nodif_B0,
            self.vent_mask_t1w,
            self.t1wtissue2dwi_xfm,
            self.vent_mask_dwi,
        )
        reg_mri_pngs(self.vent_mask_dwi, self.nodif_B0, self.qa_reg)
        reg_utils.applyxfm(
            self.nodif_B0,
            self.corpuscallosum_mask_t1w,
            self.t1wtissue2dwi_xfm,
            self.corpuscallosum_dwi,
        )
        reg_mri_pngs(self.corpuscallosum_dwi, self.nodif_B0, self.qa_reg)
        reg_utils.applyxfm(
            self.nodif_B0, self.csf_mask, self.t1wtissue2dwi_xfm, self.csf_mask_dwi
        )
        reg_mri_pngs(self.csf_mask_dwi, self.nodif_B0, self.qa_reg)
        reg_utils.applyxfm(
            self.nodif_B0, self.gm_mask, self.t1wtissue2dwi_xfm, self.gm_in_dwi
        )
        reg_mri_pngs(self.gm_in_dwi, self.nodif_B0, self.qa_reg)
        reg_utils.applyxfm(
            self.nodif_B0, self.wm_mask, self.t1wtissue2dwi_xfm, self.wm_in_dwi
        )
        reg_mri_pngs(self.wm_in_dwi, self.nodif_B0, self.qa_reg)

        # Threshold WM to binary in dwi space
        thr_img = nib.load(self.wm_in_dwi)
        thr_img.get_data()[thr_img.get_data() < 0.15] = 0
        nib.save(thr_img, self.wm_in_dwi_bin)

        # Threshold GM to binary in dwi space
        thr_img = nib.load(self.gm_in_dwi)
        thr_img.get_data()[thr_img.get_data() < 0.15] = 0
        nib.save(thr_img, self.gm_in_dwi_bin)

        # Threshold CSF to binary in dwi space
        thr_img = nib.load(self.csf_mask_dwi)
        thr_img.get_data()[thr_img.get_data() < 0.99] = 0
        nib.save(thr_img, self.csf_mask_dwi)

        # Threshold WM to binary in dwi space
        self.t_img = load_img(self.wm_in_dwi_bin)
        self.mask = math_img("img > 0", img=self.t_img)
        self.mask.to_filename(self.wm_in_dwi_bin)

        # Threshold GM to binary in dwi space
        self.t_img = load_img(self.gm_in_dwi_bin)
        self.mask = math_img("img > 0", img=self.t_img)
        self.mask.to_filename(self.gm_in_dwi_bin)

        # Threshold CSF to binary in dwi space
        self.t_img = load_img(self.csf_mask_dwi)
        self.mask = math_img("img > 0", img=self.t_img)
        self.mask.to_filename(self.csf_mask_dwi_bin)

        # Create ventricular CSF mask
        print("Creating ventricular CSF mask...")
        cmd = f"fslmaths {self.vent_mask_dwi} -kernel sphere 10 -ero -bin {self.vent_mask_dwi}"
        gen_utils.run(cmd)
        print("Creating Corpus Callosum mask...")
        cmd = f"fslmaths {self.corpuscallosum_dwi} -mas {self.wm_in_dwi_bin} -bin {self.corpuscallosum_dwi}"
        gen_utils.run(cmd)
        cmd = f"fslmaths {self.csf_mask_dwi} -add {self.vent_mask_dwi} -bin {self.vent_csf_in_dwi}"
        gen_utils.run(cmd)

        # Create gm-wm interface image
        cmd = (
            f"fslmaths {self.gm_in_dwi_bin} -mul {self.wm_in_dwi_bin} -add {self.corpuscallosum_dwi} "
            f"-sub {self.vent_csf_in_dwi} -mas {self.nodif_B0_mask} -bin {self.wm_gm_int_in_dwi}"
        )
        gen_utils.run(cmd)


text:
"""
@Time   :   2021-01-21 10:57:33
@File   :   csc_trainer.py
@Author :   Abtion
@Email  :   abtion{at}outlook.com
"""
import torch
import operator
import numpy as np
from transformers import BertTokenizer
from bbcm.utils import flatten
from bbcm.utils.evaluations import compute_corrector_prf, compute_sentence_level_prf
from .bases import BaseTrainingEngine

class CscTrainingModel(BaseTrainingEngine):
    """
        用于CSC的BaseModel, 定义了训练及预测步骤
        """

    def __init__(self, cfg, *args, **kwargs):
        super().__init__(cfg, *args, **kwargs)
        # loss weight for cor & det
        self.w = cfg.MODEL.HYPER_PARAMS[0]
        self.tokenizer = BertTokenizer.from_pretrained(cfg.MODEL.BERT_CKPT)
        # threshold for prediction judgment
        self.judge_line = 0.5

    @staticmethod
    def pt2list(_tensor):
        return _tensor.cpu().numpy().tolist()

    def training_step(self, batch, batch_idx):
        ori_text, cor_text, det_labels = batch
        outputs = self.forward(ori_text, cor_text, det_labels)
        loss = self.w * outputs[1] + (1 - self.w) * outputs[0]
        return loss

    def validation_step(self, batch, batch_idx):
        ori_text, cor_text, det_labels = batch

        # 检错loss，纠错loss，检错输出，纠错输出
        outputs = self.forward(ori_text, cor_text, det_labels)
        loss = self.w * outputs[1] + (1 - self.w) * outputs[0]
        det_y_hat = (outputs[2] > self.judge_line).long()
        cor_y_hat = torch.argmax((outputs[3]), dim=-1)
        encoded_x = self.tokenizer(cor_text, padding=True, return_tensors='pt')
        encoded_x.to(self._device)
        cor_y = encoded_x['input_ids']
        cor_y_hat *= encoded_x['attention_mask']

        results = []
        det_acc_labels = []
        cor_acc_labels = []
        for src, tgt, predict, det_predict, det_label in zip(ori_text, cor_y, cor_y_hat, det_y_hat, det_labels):
            _src = self.tokenizer(src, add_special_tokens=False)['input_ids']

            # whether a token from correction is the same as the truth token.
            c_tgt = self.pt2list(tgt[1:len(_src) + 1])
            c_predict = self.pt2list(predict[1:len(_src) + 1])
            cor_acc_arr = [int(_t == _p) for _t, _p in zip(c_tgt, c_predict)]
            cor_acc_labels.append(cor_acc_arr)

            # whether a token from detection is the same as the truth token.
            d_tgt = self.pt2list(det_label[1:len(_src) + 1])
            d_predict = self.pt2list(det_predict[1:len(_src) + 1])
            det_acc_arr = [int(_t == _p) for _t, _p in zip(d_tgt, d_predict)]
            det_acc_labels.append(det_acc_arr)

            # a tuple of lists
            results.append((_src, c_tgt, c_predict, d_predict))

        return loss.cpu().item(), det_acc_labels, cor_acc_labels, results

    def validation_epoch_end(self, outputs):
        det_acc_labels = []
        cor_acc_labels = []
        results = []

        # loss, det_acc_labels, cor_acc_labels, results
        for out in outputs:
            det_acc_labels += flatten(out[1])
            cor_acc_labels += flatten(out[2])
            results += out[3]
        loss = np.mean([out[0] for out in outputs])
        self.log('val_loss', loss)
        self._logger.info(f'loss: {loss}')
        det_acc = np.array(det_acc_labels).mean()
        # print(det_acc)
        self._logger.info(f'Detection:\n'
                          f'acc: {det_acc:.4f}')
        cor_acc = np.array(cor_acc_labels).mean()
        # print(cor_acc)
        self._logger.info(f'Correction:\n'
                          f'acc: {cor_acc:.4f}')
        det_f1, cor_f1, details = compute_corrector_prf(
            results, self._logger, on_detected=True)
        sent_acc, sent_p, sent_r, sent_f1 = compute_sentence_level_prf(
            results, self._logger)
        details.update({
            'sent_acc': sent_acc,
            'sent_f1': sent_f1,
            'sent_p': sent_p,
            'sent_r': sent_r,
        })
        self.log_dict(details, logger=True)
        return det_acc, cor_acc, det_f1, cor_f1

    def test_step(self, batch, batch_idx):
        return self.validation_step(batch, batch_idx)

    def test_epoch_end(self, outputs):
        return self.validation_epoch_end(outputs)

    def evaluate_from_loader(self, loader):
        outputs = []
        if isinstance(loader, str):
            from bbcm.data.loaders.collator import DataCollatorForCsc
            _collate_fn = DataCollatorForCsc(tokenizer=self.tokenizer)

            from bbcm.data.build import make_loaders
            trn, dev, tst = make_loaders(cfg=self.cfg,
                                         _collate_fn=_collate_fn)
            loader = {
                'train': trn,
                'dev': dev,
                'valid': dev,
                'test': tst
            }.get(loader, tst)
        for b_idx, batch in enumerate(loader):
            ori_texts, cor_texts, det_labels = batch
            batch = (ori_texts, cor_texts, det_labels.to(self.cfg.MODEL.DEVICE))
            outputs.append(self.validation_step(batch, b_idx))
        return self.validation_epoch_end(outputs)

    def predict(self, texts, detail=False):
        inputs = self.tokenizer(texts, padding=True, return_tensors='pt')
        inputs.to(self.cfg.MODEL.DEVICE)
        with torch.no_grad():
            # 检测输出，纠错输出
            outputs = self.forward(texts)
            y_hat = torch.argmax(outputs[1], dim=-1)
            expand_text_lens = torch.sum(inputs['attention_mask'], dim=-1) - 1
        rst = []
        for t_len, _y_hat in zip(expand_text_lens, y_hat):
            rst.append(self.tokenizer.decode(_y_hat[1:t_len]).replace(' ', ''))
        if detail:
            return outputs, rst
        return rst


text:
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2020
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
Clifford operator class.
"""
import re
import numpy as np

from qiskit.exceptions import QiskitError
from qiskit.circuit import QuantumCircuit, Instruction
from qiskit.circuit.library.standard_gates import IGate, XGate, YGate, ZGate, HGate, SGate
from qiskit.quantum_info.operators.base_operator import BaseOperator
from qiskit.quantum_info.operators.operator import Operator
from qiskit.quantum_info.operators.scalar_op import ScalarOp
from qiskit.quantum_info.synthesis.clifford_decompose import decompose_clifford
from qiskit.quantum_info.operators.mixins import generate_apidocs, AdjointMixin
from .stabilizer_table import StabilizerTable
from .clifford_circuits import _append_circuit

class Clifford(BaseOperator, AdjointMixin):
    """An N-qubit unitary operator from the Clifford group.

    **Representation**

    An *N*-qubit Clifford operator is stored as a length *2N*
    :class:`~qiskit.quantum_info.StabilizerTable` using the convention
    from reference [1].

    * Rows 0 to *N-1* are the *destabilizer* group generators
    * Rows *N* to *2N-1* are the *stabilizer* group generators.

    The internal :class:`~qiskit.quantum_info.StabilizerTable` for the Clifford
    can be accessed using the :attr:`table` attribute. The destabilizer or
    stabilizer rows can each be accessed as a length-N Stabilizer table using
    :attr:`destabilizer` and :attr:`stabilizer` attributes.

    A more easily human readable representation of the Clifford operator can
    be obtained by calling the :meth:`to_dict` method. This representation is
    also used if a Clifford object is printed as in the following example

    .. jupyter-execute::

        from qiskit import QuantumCircuit
        from qiskit.quantum_info import Clifford

        # Bell state generation circuit
        qc = QuantumCircuit(2)
        qc.h(0)
        qc.cx(0, 1)
        cliff = Clifford(qc)

        # Print the Clifford
        print(cliff)

        # Print the Clifford destabilizer rows
        print(cliff.destabilizer)

        # Print the Clifford stabilizer rows
        print(cliff.stabilizer)

    **Circuit Conversion**

    Clifford operators can be initialized from circuits containing *only* the
    following Clifford gates: :class:`~qiskit.circuit.library.IGate`,
    :class:`~qiskit.circuit.library.XGate`, :class:`~qiskit.circuit.library.YGate`,
    :class:`~qiskit.circuit.library.ZGate`, :class:`~qiskit.circuit.library.HGate`,
    :class:`~qiskit.circuit.library.SGate`, :class:`~qiskit.circuit.library.SdgGate`,
    :class:`~qiskit.circuit.library.CXGate`, :class:`~qiskit.circuit.library.CZGate`,
    :class:`~qiskit.circuit.library.SwapGate`.
    They can be converted back into a :class:`~qiskit.circuit.QuantumCircuit`,
    or :class:`~qiskit.circuit.Gate` object using the :meth:`~Clifford.to_circuit`
    or :meth:`~Clifford.to_instruction` methods respectively. Note that this
    decomposition is not necessarily optimal in terms of number of gates.

    .. note::

        A minimally generating set of gates for Clifford circuits is
        the :class:`~qiskit.circuit.library.HGate` and
        :class:`~qiskit.circuit.library.SGate` gate and *either* the
        :class:`~qiskit.circuit.library.CXGate` or
        :class:`~qiskit.circuit.library.CZGate` two-qubit gate.

    Clifford operators can also be converted to
    :class:`~qiskit.quantum_info.Operator` objects using the
    :meth:`to_operator` method. This is done via decomposing to a circuit, and then
    simulating the circuit as a unitary operator.

    References:
        1. S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
           Phys. Rev. A 70, 052328 (2004).
           `arXiv:quant-ph/0406196 <https://arxiv.org/abs/quant-ph/0406196>`_
    """

    def __array__(self, dtype=None):
        if dtype:
            return np.asarray(self.to_matrix(), dtype=dtype)
        return self.to_matrix()

    def __init__(self, data, validate=True):
        """Initialize an operator object."""

        # Initialize from another Clifford by sharing the underlying
        # StabilizerTable
        if isinstance(data, Clifford):
            self._table = data._table

        # Initialize from ScalarOp as N-qubit identity discarding any global phase
        elif isinstance(data, ScalarOp):
            if not data.num_qubits or not data.is_unitary():
                raise QiskitError("Can only initialize from N-qubit identity ScalarOp.")
            self._table = StabilizerTable(np.eye(2 * data.num_qubits, dtype=bool))

        # Initialize from a QuantumCircuit or Instruction object
        elif isinstance(data, (QuantumCircuit, Instruction)):
            self._table = Clifford.from_circuit(data)._table

        # Initialize StabilizerTable directly from the data
        else:
            self._table = StabilizerTable(data)

            # Validate table is a symplectic matrix
            if validate and not Clifford._is_symplectic(self._table.array):
                raise QiskitError(
                    "Invalid Clifford. Input StabilizerTable is not a valid symplectic matrix."
                )

        # Initialize BaseOperator
        super().__init__(num_qubits=self._table.num_qubits)

    def __repr__(self):
        return f"Clifford({repr(self.table)})"

    def __str__(self):
        return "Clifford: Stabilizer = {}, Destabilizer = {}".format(
            str(self.stabilizer.to_labels()), str(self.destabilizer.to_labels())
        )

    def __eq__(self, other):
        """Check if two Clifford tables are equal"""
        return super().__eq__(other) and self._table == other._table

    # ---------------------------------------------------------------------
    # Attributes
    # ---------------------------------------------------------------------
    def __getitem__(self, key):
        """Return a stabilizer Pauli row"""
        return self._table.__getitem__(key)

    def __setitem__(self, key, value):
        """Set a stabilizer Pauli row"""
        self._table.__setitem__(key, value)

    @property
    def table(self):
        """Return StabilizerTable"""
        return self._table

    @table.setter
    def table(self, value):
        """Set the stabilizer table"""
        # Note this setter cannot change the size of the Clifford
        # It can only replace the contents of the StabilizerTable with
        # another StabilizerTable of the same size.
        if not isinstance(value, StabilizerTable):
            value = StabilizerTable(value)
        self._table._array[:, :] = value._table._array
        self._table._phase[:] = value._table._phase

    @property
    def stabilizer(self):
        """Return the stabilizer block of the StabilizerTable."""
        return StabilizerTable(self._table[self.num_qubits : 2 * self.num_qubits])

    @stabilizer.setter
    def stabilizer(self, value):
        """Set the value of stabilizer block of the StabilizerTable"""
        inds = slice(self.num_qubits, 2 * self.num_qubits)
        self._table.__setitem__(inds, value)

    @property
    def destabilizer(self):
        """Return the destabilizer block of the StabilizerTable."""
        return StabilizerTable(self._table[0 : self.num_qubits])

    @destabilizer.setter
    def destabilizer(self, value):
        """Set the value of destabilizer block of the StabilizerTable"""
        inds = slice(0, self.num_qubits)
        self._table.__setitem__(inds, value)

    # ---------------------------------------------------------------------
    # Utility Operator methods
    # ---------------------------------------------------------------------

    def is_unitary(self):
        """Return True if the Clifford table is valid."""
        # A valid Clifford is always unitary, so this function is really
        # checking that the underlying Stabilizer table array is a valid
        # Clifford array.
        return Clifford._is_symplectic(self.table.array)

    # ---------------------------------------------------------------------
    # BaseOperator Abstract Methods
    # ---------------------------------------------------------------------

    def conjugate(self):
        return Clifford._conjugate_transpose(self, "C")

    def adjoint(self):
        return Clifford._conjugate_transpose(self, "A")

    def transpose(self):
        return Clifford._conjugate_transpose(self, "T")

    def tensor(self, other):
        if not isinstance(other, Clifford):
            other = Clifford(other)
        return self._tensor(self, other)

    def expand(self, other):
        if not isinstance(other, Clifford):
            other = Clifford(other)
        return self._tensor(other, self)

    @classmethod
    def _tensor(cls, a, b):
        # Pad stabilizers and destabilizers
        destab = b.destabilizer.expand(a.num_qubits * "I") + a.destabilizer.tensor(
            b.num_qubits * "I"
        )
        stab = b.stabilizer.expand(a.num_qubits * "I") + a.stabilizer.tensor(b.num_qubits * "I")

        # Add the padded table
        return Clifford(destab + stab, validate=False)

    def compose(self, other, qargs=None, front=False):
        if qargs is None:
            qargs = getattr(other, "qargs", None)
        # If other is a QuantumCircuit we can more efficiently compose
        # using the _append_circuit method to update each gate recursively
        # to the current Clifford, rather than converting to a Clifford first
        # and then doing the composition of tables.
        if not front and isinstance(other, (QuantumCircuit, Instruction)):
            ret = self.copy()
            _append_circuit(ret, other, qargs=qargs)
            return ret

        if not isinstance(other, Clifford):
            other = Clifford(other)

        # Validate compose dimensions
        self._op_shape.compose(other._op_shape, qargs, front)

        # Pad other with identities if composing on subsystem
        other = self._pad_with_identity(other, qargs)

        if front:
            table1 = self.table
            table2 = other.table
        else:
            table1 = other.table
            table2 = self.table

        num_qubits = self.num_qubits

        array1 = table1.array.astype(int)
        phase1 = table1.phase.astype(int)

        array2 = table2.array.astype(int)
        phase2 = table2.phase.astype(int)

        # Update Pauli table
        pauli = StabilizerTable(array2.dot(array1) % 2)

        # Add phases
        phase = np.mod(array2.dot(phase1) + phase2, 2)

        # Correcting for phase due to Pauli multiplication
        ifacts = np.zeros(2 * num_qubits, dtype=int)

        for k in range(2 * num_qubits):

            row2 = array2[k]
            x2 = table2.X[k]
            z2 = table2.Z[k]

            # Adding a factor of i for each Y in the image of an operator under the
            # first operation, since Y=iXZ

            ifacts[k] += np.sum(x2 & z2)

            # Adding factors of i due to qubit-wise Pauli multiplication

            for j in range(num_qubits):
                x = 0
                z = 0
                for i in range(2 * num_qubits):
                    if row2[i]:
                        x1 = array1[i, j]
                        z1 = array1[i, j + num_qubits]
                        if (x | z) & (x1 | z1):
                            val = np.mod(np.abs(3 * z1 - x1) - np.abs(3 * z - x) - 1, 3)
                            if val == 0:
                                ifacts[k] += 1
                            elif val == 1:
                                ifacts[k] -= 1
                        x = np.mod(x + x1, 2)
                        z = np.mod(z + z1, 2)

        p = np.mod(ifacts, 4) // 2

        phase = np.mod(phase + p, 2)

        return Clifford(StabilizerTable(pauli, phase), validate=False)

    # ---------------------------------------------------------------------
    # Representation conversions
    # ---------------------------------------------------------------------

    def to_dict(self):
        """Return dictionary representation of Clifford object."""
        return {
            "stabilizer": self.stabilizer.to_labels(),
            "destabilizer": self.destabilizer.to_labels(),
        }

    @staticmethod
    def from_dict(obj):
        """Load a Clifford from a dictionary"""
        destabilizer = StabilizerTable.from_labels(obj.get("destabilizer"))
        stabilizer = StabilizerTable.from_labels(obj.get("stabilizer"))
        return Clifford(destabilizer + stabilizer)

    def to_matrix(self):
        """Convert operator to Numpy matrix."""
        return self.to_operator().data

    def to_operator(self):
        """Convert to an Operator object."""
        return Operator(self.to_instruction())

    def to_circuit(self):
        """Return a QuantumCircuit implementing the Clifford.

        For N <= 3 qubits this is based on optimal CX cost decomposition
        from reference [1]. For N > 3 qubits this is done using the general
        non-optimal compilation routine from reference [2].

        Return:
            QuantumCircuit: a circuit implementation of the Clifford.

        References:
            1. S. Bravyi, D. Maslov, *Hadamard-free circuits expose the
               structure of the Clifford group*,
               `arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_

            2. S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
               Phys. Rev. A 70, 052328 (2004).
               `arXiv:quant-ph/0406196 <https://arxiv.org/abs/quant-ph/0406196>`_
        """
        return decompose_clifford(self)

    def to_instruction(self):
        """Return a Gate instruction implementing the Clifford."""
        return self.to_circuit().to_gate()

    @staticmethod
    def from_circuit(circuit):
        """Initialize from a QuantumCircuit or Instruction.

        Args:
            circuit (QuantumCircuit or ~qiskit.circuit.Instruction):
                instruction to initialize.

        Returns:
            Clifford: the Clifford object for the instruction.

        Raises:
            QiskitError: if the input instruction is non-Clifford or contains
                         classical register instruction.
        """
        if not isinstance(circuit, (QuantumCircuit, Instruction)):
            raise QiskitError("Input must be a QuantumCircuit or Instruction")

        # Convert circuit to an instruction
        if isinstance(circuit, QuantumCircuit):
            circuit = circuit.to_instruction()

        # Initialize an identity Clifford
        clifford = Clifford(np.eye(2 * circuit.num_qubits), validate=False)
        _append_circuit(clifford, circuit)
        return clifford

    @staticmethod
    def from_label(label):
        """Return a tensor product of single-qubit Clifford gates.

        Args:
            label (string): single-qubit operator string.

        Returns:
            Clifford: The N-qubit Clifford operator.

        Raises:
            QiskitError: if the label contains invalid characters.

        Additional Information:
            The labels correspond to the single-qubit Cliffords are

            * - Label
              - Stabilizer
              - Destabilizer
            * - ``"I"``
              - +Z
              - +X
            * - ``"X"``
              - -Z
              - +X
            * - ``"Y"``
              - -Z
              - -X
            * - ``"Z"``
              - +Z
              - -X
            * - ``"H"``
              - +X
              - +Z
            * - ``"S"``
              - +Z
              - +Y
        """
        # Check label is valid
        label_gates = {
            "I": IGate(),
            "X": XGate(),
            "Y": YGate(),
            "Z": ZGate(),
            "H": HGate(),
            "S": SGate(),
        }
        if re.match(r"^[IXYZHS\-+]+$", label) is None:
            raise QiskitError("Label contains invalid characters.")
        # Initialize an identity matrix and apply each gate
        num_qubits = len(label)
        op = Clifford(np.eye(2 * num_qubits, dtype=bool))
        for qubit, char in enumerate(reversed(label)):
            _append_circuit(op, label_gates[char], qargs=[qubit])
        return op

    # ---------------------------------------------------------------------
    # Internal helper functions
    # ---------------------------------------------------------------------

    @staticmethod
    def _is_symplectic(mat):
        """Return True if input is symplectic matrix."""
        # Condition is
        # table.T * [[0, 1], [1, 0]] * table = [[0, 1], [1, 0]]
        # where we are block matrix multiplying using symplectic product

        dim = len(mat) // 2
        if mat.shape != (2 * dim, 2 * dim):
            return False

        one = np.eye(dim, dtype=int)
        zero = np.zeros((dim, dim), dtype=int)
        seye = np.block([[zero, one], [one, zero]])
        arr = mat.astype(int)
        return np.array_equal(np.mod(arr.T.dot(seye).dot(arr), 2), seye)

    @staticmethod
    def _conjugate_transpose(clifford, method):
        """Return the adjoint, conjugate, or transpose of the Clifford.

        Args:
            clifford (Clifford): a clifford object.
            method (str): what function to apply 'A', 'C', or 'T'.

        Returns:
            Clifford: the modified clifford.
        """
        ret = clifford.copy()
        if method in ["A", "T"]:
            # Apply inverse
            # Update table
            tmp = ret.destabilizer.X.copy()
            ret.destabilizer.X = ret.stabilizer.Z.T
            ret.destabilizer.Z = ret.destabilizer.Z.T
            ret.stabilizer.X = ret.stabilizer.X.T
            ret.stabilizer.Z = tmp.T
            # Update phase
            ret.table.phase ^= clifford.dot(ret).table.phase
        if method in ["C", "T"]:
            # Apply conjugate
            ret.table.phase ^= np.mod(np.sum(ret.table.X & ret.table.Z, axis=1), 2).astype(bool)
        return ret

    def _pad_with_identity(self, clifford, qargs):
        """Pad Clifford with identities on other subsystems."""
        if qargs is None:
            return clifford

        padded = Clifford(StabilizerTable(np.eye(2 * self.num_qubits, dtype=bool)), validate=False)

        inds = list(qargs) + [self.num_qubits + i for i in qargs]

        # Pad Pauli array
        pauli = clifford.table.array
        for i, pos in enumerate(qargs):
            padded.table.array[inds, pos] = pauli[:, i]
            padded.table.array[inds, self.num_qubits + pos] = pauli[:, clifford.num_qubits + i]

        # Pad phase
        padded.table.phase[inds] = clifford.table.phase

        return padded

# Update docstrings for API docs
generate_apidocs(Clifford)