diff --git a/pyproject.toml b/pyproject.toml index 241aca3..f532d52 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,23 @@ [build-system] requires = ["setuptools"] build-backend = "setuptools.build_meta" + +[project] +name = "quactography" +version = "0.1" +description = "Quactography : quantum tracking" +authors = [{ name = "SCIL Team" }] +readme = "README.md" +requires-python = ">=3.10, <3.13" +license-files = ["LICENSE"] +dynamic = ["dependencies", "scripts"] +classifiers = [ + "Development Status :: 3 - Alpha", + "Environment :: Console", + "Intended Audience :: Science/Research", + "Operating System :: OS Independent", + "Programming Language :: Python", + "Topic :: Scientific/Engineering" +] + + diff --git a/quactography/adj_matrix/filter.py b/quactography/adj_matrix/filter.py index f55af77..a8264f1 100644 --- a/quactography/adj_matrix/filter.py +++ b/quactography/adj_matrix/filter.py @@ -160,4 +160,4 @@ def remove_zero_columns_rows(mat: np.ndarray): zero_rows = np.all(mat == 0, axis=1) non_zero_cols = np.where(~zero_cols)[0] non_zero_rows = np.where(~zero_rows)[0] - return mat[np.ix_(non_zero_rows, non_zero_cols)] + return mat[np.ix_(non_zero_rows, non_zero_cols)] \ No newline at end of file diff --git a/quactography/adj_matrix/reconst.py b/quactography/adj_matrix/reconst.py index fef970b..af43ebe 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -1,6 +1,8 @@ import numpy as np from dipy.reconst.shm import sh_to_sf from dipy.core.sphere import Sphere +from quactography.graph.utils import get_output_nodes + def build_adjacency_matrix(nodes_mask): @@ -16,6 +18,7 @@ def build_adjacency_matrix(nodes_mask): adj_matrix : np.ndarray Adjacency matrix of the nodes. """ + # 1st. Assign labels to non-zero voxels (nodes) node_indices = np.flatnonzero(nodes_mask) @@ -47,7 +50,18 @@ def build_adjacency_matrix(nodes_mask): x + x_offset, y + y_offset, z + z_offset, i, nodes_mask, labels_volume, adj_matrix ) - return adj_matrix, node_indices + # Adds possibility of an edge to the actual node closest neighbour in 26 directions + for x_offset in [-1, 0, 1]: + for y_offset in [-1, 0, 1]: + for z_offset in [-1, 0, 1]: + if x_offset == 0 and y_offset == 0 and z_offset == 0: + continue + else: + adj_matrix = _add_edge_perhaps( + x + x_offset, y + y_offset, z + z_offset, i, nodes_mask, labels_volume, adj_matrix + ) + + return adj_matrix, node_indices, labels_volume def build_weighted_graph(adj_matrix, node_indices, sh, sh_order=12): @@ -135,6 +149,33 @@ def _add_edge_perhaps( return adj_matrix +def add_end_point_edge(adj_matrix, end, labels): + """ + Add a node and edges to the end points of the ROI in the adjacency matrix. + + Parameters + ---------- + adj_matrix : np.ndarray + Adjacency matrix of the graph. + end_points : list of tuples + List of end points as (x, y, z) coordinates. + labels_volume : np.ndarray + List of every column in the image: + + Returns + ------- + np.ndarray + Updated adjacency matrix with end point edges added. + """ + + adj_matrix = np.lib.pad(adj_matrix, (0, 1), 'constant', constant_values=(0)) + for i in end: + start = labels[i[0], i[1], i[2]] + adj_matrix[start, -1] = 1 + adj_matrix[-1, start] = 1 # Assuming undirected graph + return adj_matrix + + def _is_valid_pos(pos_x, pos_y, pos_z, nodes_mask): if pos_x < 0 or pos_x >= nodes_mask.shape[0]: return False @@ -143,3 +184,16 @@ def _is_valid_pos(pos_x, pos_y, pos_z, nodes_mask): if pos_z < 0 or pos_z >= nodes_mask.shape[2]: return False return nodes_mask[pos_x, pos_y, pos_z] + +if __name__ == "__main__": + # little demo code here and proto test + mask = np.zeros((5, 5, 5), dtype=bool) + mask[1:3, 1:3, 1:3] = True + mat, nodes, labels = build_adjacency_matrix(mask) + entry_node = np.array([1, 1, 1]) + propagation_direction = np.array([0, 1, 1]) + angle_rad = np.pi / 8 + + indices = get_output_nodes(mask, entry_node, propagation_direction, angle_rad) + new = add_end_point_edge(mat,indices,labels) + print(new) diff --git a/quactography/graph/utils.py b/quactography/graph/utils.py index 5498fec..0f25795 100644 --- a/quactography/graph/utils.py +++ b/quactography/graph/utils.py @@ -48,3 +48,15 @@ def get_output_nodes(mask, entry_node, propagation_direction, angle_rad): is_valid_edge = directions_cos_angle > propagation_cos_angle valid_indices = edge_indices[is_valid_edge.reshape((-1,))] return valid_indices + + +if __name__ == '__main__': + # little demo code here + mask = np.zeros((25, 25, 25), dtype=bool) + mask[5:10, 5:10, 5:10] = True + entry_node = np.array([5, 5, 5]) + propagation_direction = np.array([0, 1, 1]) + angle_rad = np.pi / 8 + + indices = get_output_nodes(mask, entry_node, propagation_direction, angle_rad) + print(indices) diff --git a/quactography/hamiltonian/hamiltonian_qubit_edge.py b/quactography/hamiltonian/hamiltonian_qubit_edge.py index 4848e3e..abff716 100644 --- a/quactography/hamiltonian/hamiltonian_qubit_edge.py +++ b/quactography/hamiltonian/hamiltonian_qubit_edge.py @@ -328,22 +328,20 @@ def get_exact_sol(self): # # # # TEST:-------------------------------------------------------------------------- +# # from quactography.hamiltonian.hamiltonian_qubit_node import Hamiltonian_qubit_node +# import numpy as np -# # # mat = np.array([[0, 1, 1, 0], [1, 0, 0, 5], [1, 0, 0, 6], [0, 5, 6, 0]]) +# from quactography.adj_matrix.io import save_graph +# mat = np.array([[0, 1, 1, 0], [1, 0, 0, 5], [1, 0, 0, 6], [0, 5, 6, 0]]) # # # # # This is the given format you should use to save the graph, for mat: -# # # # save_graph(mat, np.array([0, 1, 2, 3]), np.array([4, 4]), "rand_graph.npz") +# save_graph(mat, np.array([0, 1, 2, 3]), np.array([4, 4]), "rand_graph.npz") # import sys -# sys.path.append(r"C:\Users\harsh\quactography") - # from quactography.graph.undirected_graph import Graph # from quactography.adj_matrix.io import load_graph -# # from quactography.hamiltonian.hamiltonian_qubit_node import Hamiltonian_qubit_node -# import numpy as np -# from quactography.adj_matrix.io import save_graph # my_graph_class = Graph( # np.array( @@ -407,14 +405,14 @@ def get_exact_sol(self): # print_hamiltonian_circuit(h.hint_c, ["10101"]) -# # # print("total2") -# # # print_hamiltonian_circuit(h.total_hamiltonian, ["11111"]) -# # # print("mandatory2") -# # # print_hamiltonian_circuit(h.mandatory_c, ["11111"]) -# # # print("start2") -# # # print_hamiltonian_circuit(h.starting_node_c, ["11111"]) -# # # print("finish2") -# # # print_hamiltonian_circuit(h.ending_node_c, ["11111"]) -# # # print("int2") -# # # print_hamiltonian_circuit(h.hint_c, ["11111"]) +# print("total2") +# print_hamiltonian_circuit(h.total_hamiltonian, ["11111"]) +# print("mandatory2") +# print_hamiltonian_circuit(h.mandatory_c, ["11111"]) +# print("start2") +# print_hamiltonian_circuit(h.starting_node_c, ["11111"]) +# print("finish2") +# print_hamiltonian_circuit(h.ending_node_c, ["11111"]) +# print("int2") +# print_hamiltonian_circuit(h.hint_c, ["11111"]) # # # ------------------------------------------------------------------------------------------------------------------------------------------ diff --git a/quactography/solver/Dijkstra.py b/quactography/solver/Dijkstra.py new file mode 100644 index 0000000..d800fad --- /dev/null +++ b/quactography/solver/Dijkstra.py @@ -0,0 +1,74 @@ +import heapq +import networkx as nx + + +def dijkstra_stepwise(Graph, start, target): + """ + Perform Dijkstra's algorithm in a stepwise manner on a given graph. + Parameters + ---------- + G : networkx.Graph + The input graph where each node is typically a tuple (e.g., (x, y)). + Edges may have a 'weight' attribute indicating traversal cost + (default is 1). + start : hashable + The starting node for the pathfinding algorithm. + target : hashable + The target node to reach. + diagonal_mode : str, optional + Neighbor retrieval mode. Must be either: + - "nondiagonal": only 4-directional neighbors (up, down, left, right) + - "diagonal": include diagonal neighbors (8 directions total) + Default is "nondiagonal". + Returns + ------- + evaluated_nodes : list + The list of nodes in the order they were evaluated by the algorithm. + path_to_current : list of list + A list containing the reconstructed path from the start node to each + node evaluated so far, in evaluation order. + current_distance : float + The final distance to the target node if found, or to the last + evaluated node otherwise. Returns None if no path is found. + """ + + G = nx.from_numpy_array(Graph) + distances = {node: float('inf') for node in G.nodes()} + distances[start] = 0 + previous_nodes = {node: None for node in G.nodes()} + evaluated_nodes = [] + path_to_current = [] + priority_queue = [(0, start)] + heapq.heapify(priority_queue) + + while priority_queue: + current_distance, current_node = heapq.heappop(priority_queue) + + if current_node not in evaluated_nodes: + evaluated_nodes.append(current_node) + + temp_path = [] + node = current_node + while node is not None: + temp_path.append(node) + if node in previous_nodes: + node = previous_nodes[node] + else: + break + temp_path.reverse() + path_to_current.append(temp_path) + + if current_node == target: + break + + neighbors = list(G.neighbors(current_node)) + + for neighbor in neighbors: + if neighbor not in evaluated_nodes: + edge_weight = G[current_node][neighbor].get("weight", 1) + new_distance = current_distance + edge_weight + if new_distance < distances[neighbor]: + distances[neighbor] = new_distance + previous_nodes[neighbor] = current_node + heapq.heappush(priority_queue, (new_distance, neighbor)) + return path_to_current[-1] diff --git a/quactography/solver/qaoa_solver_qu_edge.py b/quactography/solver/qaoa_solver_qu_edge.py index 89c1058..164248c 100644 --- a/quactography/solver/qaoa_solver_qu_edge.py +++ b/quactography/solver/qaoa_solver_qu_edge.py @@ -190,18 +190,24 @@ def find_longest_path(args): opt_params=resx # type: ignore ) # type: ignore -def find_max_cost(in_folder, alpha, rep): +def find_min_cost(in_folder, alpha, rep): """ - Remove all optimization results files in the input folder - that do not have the minimum cost for a given alpha and repetition. + Find the minimum cost in the optimization results for given alpha and rep, + and remove all files with a cost lower than the minimum cost. + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + QUACK_RAP does not save any files, nor have iteration with different alphas or repetitions, + it only returns a list of coordinates for the streamline. + Parameters ---------- - in_folder: str - The input folder containing the optimization results files in .npz format. - alpha: float - The alpha value to filter the optimization results files. - rep: int - The repetition value to filter the optimization results files. + in_folder : str + The folder containing the optimization results files. + alpha : float + The alpha value used in the optimization. + rep : int + The number of repetitions for the QAOA algorithm. + Returns ------- None @@ -228,7 +234,7 @@ def multiprocess_qaoa_solver_edge( output_folder, optimizer, cost_landscape, - save_only, + save_only ): """ Solve the optimization problem using the QAOA algorithm @@ -271,7 +277,119 @@ def multiprocess_qaoa_solver_edge( itertools.repeat(save_only), ), ) + print( + "------------------MULTIPROCESS SOLVER FINISHED-------------------------" + ) + + + +def multiprocess_qaoa_solver_edge_rap( + hamiltonians, + reps, + nbr_processes, + optimizer +): + """ + Solve the optimization problem using the QAOA algorithm + with multiprocessing on the alpha values. + Parameters + ---------- + hamiltonians : list + List of Hamiltonian objects from quactography library, Hamiltonian_qubit_edge. + batch_count : int + Number of time the command will be ran + reps : int + Number of repetitions for the QAOA algorithm, + determines the number of sets of gamma and beta angles. + nbr_processes : int + Number of cpu to use for multiprocessing. default=1 + output_file : str + The output file name for the optimization results in .npz format. + optimizer : str + Optimizer to use for the QAOA algorithm. default="Differential" + cost_landscape : bool + Plot the cost landscape with the optimal point if reps=1. default=False + save_only : bool + If True, the figure is saved without displaying it. default=False + + Returns + ------- + None + """ + pool = multiprocessing.Pool(nbr_processes) + pool.map( + find_longest_path, + zip( + hamiltonians, + itertools.repeat(reps), + itertools.repeat(optimizer), + ), + ) print( "------------------MULTIPROCESS SOLVER FINISHED-------------------------" ) + +def find_longest_path_rap(args): + """ + Find the longest path in a graph using the QAOA algorithm, + with a plot of the cost landscape if reps=1 and + the optimal point if cost_landscape=True. + Parameters + ---------- + args : tuple + Tuple containing the Hamiltonian object from quactography library, + Hamiltonian_qubit_edge, the number of repetitions for the QAOA algorithm, + the output file name for the optimization results in .npz format, the optimizer + to use for the QAOA algorithm, a boolean to plot the cost landscape with + the optimal point if reps=1, and a boolean to save the figure + without displaying it. + Returns + ------- + None + """ + h = args[0] + reps = args[1] + optimizer = args[2] + mixerf = mixer(h) + # Create QAOA circuit. + ansatz = QAOAAnsatz(h.total_hamiltonian, reps, mixer_operator=mixerf, name="QAOA", flatten=True) + # Plot the circuit layout: + # ansatz.decompose(reps=3).draw() + # ----------------------------------------------------------------RUN LOCALLY: ----- + # Run on local estimator and sampler: + # Save output file name diffrerent for each alpha and loop: + estimator = Estimator(options={"shots": 1000000, "seed": 43}) + sampler = Sampler(options={"shots": 1000000, "seed": 43}) + # ----------------------------------------------------------------------------------- + if optimizer == "Differential": + # Reference: https://www.youtube.com/watch?v=o-OPrQmS1pU + # Define fixed arguments + cost_func_with_args = partial( + cost_func, + estimator=estimator, + ansatz=ansatz, + hamiltonian=h.total_hamiltonian, + ) + # Call differential evolution with the modified cost function + bounds = [[0, 2 * np.pi], [0, np.pi]] * reps + res = differential_evolution(cost_func_with_args, bounds, disp=False) + resx = res.x + + # ----------------------------------------------------- + circ = ansatz.copy() + circ.measure_all() + + dist = sampler.run(circ, resx).result().quasi_dists[0] # type: ignore + bin_str = list(map(int, max(dist.binary_probabilities(), + key=dist.binary_probabilities().get))) # type: ignore + + bin_str_reversed = bin_str[::-1] + bin_str_reversed = np.array(bin_str_reversed) # type: ignore + + # Concatenate the binary path to a string: + str_path_reversed = ["".join(map(str, bin_str_reversed))] # type: ignore + str_path_reversed = str_path_reversed[0] # type: ignore + # Save parameters alpha and min_cost with path in csv file: + opt_path = str_path_reversed + return opt_path diff --git a/quactography/solver/rap_tracking.py b/quactography/solver/rap_tracking.py new file mode 100644 index 0000000..302e18b --- /dev/null +++ b/quactography/solver/rap_tracking.py @@ -0,0 +1,189 @@ +import nibabel as nib +import numpy as np +from scipy.sparse.csgraph import dijkstra +import matplotlib.pyplot as plt +import networkx as nx + + +from quactography.adj_matrix.reconst import ( + build_adjacency_matrix, + build_weighted_graph, + add_end_point_edge +) +from quactography.adj_matrix.filter import ( + remove_orphan_nodes, + remove_intermediate_connections, + extract_slice_at_index +) +from quactography.graph.undirected_graph import Graph +from quactography.graph.utils import get_output_nodes +from quactography.hamiltonian.hamiltonian_qubit_edge import Hamiltonian_qubit_edge +from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge_rap +from quactography.solver.Dijkstra import dijkstra_stepwise + + +def quack_rap(rap_mask_img, rap_sh_img, start_point, reps, alpha, + keep_mask=None, threshold=0.2, slice_index=None, + axis_name="axial", sh_order=8, prev_direction=[0,0,0], theta=45): + """Build adjacency matrix from diffusion data (white matter mask and fodf peaks). + Parameters + ---------- + in_nodes_mask : str + Input nodes mask image (.nii.gz file). + in_sh : str + Input SH image (.nii.gz file). + start_point : int + Starting node index in the graph. + keep_mask : str, optional + Nodes that must not be filtered out. If None, all nodes are filtered. + threshold : float, optional + Cut all weights below a given threshold. Default is 0.2. + slice_index : int, optional + If None, a 3D graph is built. If an integer, a slice is extracted + along the specified axis. + axis_name : str, optional + Axis along which a slice is taken. Default is "axial". + sh_order : int, optional + Maximum SH order. Default is 8. + prev_direction : list, optional + Previous direction of the streamline, used to determine the propagation direction. + Default is [0, 0, 0]. + theta : float, optional + Aperture angle in degrees for the propagation direction. Default is 45. + + Returns + ------- + line : list + List of coordinates for the streamline. + """ + + + #nodes_mask_im = nib.load(in_nodes_mask) + #sh_im = nib.load(in_sh) + + nodes_mask = rap_mask_img.get_fdata().astype(bool) + sh = rap_sh_img.get_fdata() + + keep_node_indices = None + if keep_mask: + keep_mask = nib.load(keep_mask).get_fdata().astype(bool) + keep_node_indices = np.flatnonzero(keep_mask) + + # adjacency graph + adj_matrix, node_indices, labes = build_adjacency_matrix(nodes_mask) + + # assign edge weights + weighted_graph, node_indices = build_weighted_graph( + adj_matrix, node_indices, sh, sh_order + ) + + # filter graph edges by weight + weighted_graph[weighted_graph < threshold] = 0.0 + + + # # remove intermediate nodes that connect only two nodes + # weighted_graph = remove_intermediate_connections( + # weighted_graph, node_indices, keep_node_indices + # ) + + # # remove nodes without edges + # weighted_graph, node_indices = remove_orphan_nodes( + # weighted_graph, node_indices, keep_node_indices + # ) + + if slice_index is not None: + weighted_graph, node_indices = extract_slice_at_index( + weighted_graph, node_indices, nodes_mask.shape, slice_index, axis_name + ) + + # Get end points of the streamline + end_pointX = np.array([[30,17,2]]) + end_points = get_output_nodes( + nodes_mask, + entry_node=start_point, + propagation_direction=prev_direction, + angle_rad=theta + ) + # Add a node that connects to all end points with neutral weight in order to know where to get out the region + weighted_graph = add_end_point_edge(weighted_graph, end_pointX, labels=labes) + end = weighted_graph.shape[0]-1 + + #function to process the graph before quantum path finding + line = rap_funct( + weighted_graph, + starting_node = labes[start_point[0], start_point[1], start_point[2]], + ending_node = end, + alphas = [alpha], + reps = reps, + ) + + # reconstruct the streamline coordinates from the node indices + line.pop() + llist = [] + for i in line: + sline = np.unravel_index(node_indices[i],labes.shape) + llist.append(np.array(sline)) + llist.append(np.add(llist[-1], prev_direction)) + + return llist, prev_direction, True + +def rap_funct(weighted_graph, starting_node, ending_node, alphas, + reps, number_processors=2, optimizer="Differential"): + """ + Process he Graph in order to create the Hamiltonian matrix before optimization + with QAOA algorithm. The Hamiltonian is constructed with qubits as edges. + + Parameters + ---------- + graph : str + Path to the input graph file (npz file). + starting_node : int + Starting node of the graph. + ending_node : int + Ending node of the graph. + alphas : list of float, optional + List of alpha values for the Hamiltonian. Default is [1.2]. + reps : int, optional + Number of repetitions for the QAOA algorithm, determines the number + of sets of gamma and beta angles. Default is 1. + number_processors : int, optional + Number of CPU to use for multiprocessing. Default is 2. + optimizer : str, optional + Optimizer to use for the QAOA algorithm. Default is "Differential". + Returns + ------- + line : list + List of coordinates for the streamline. + """ + + #---------------QUTANTUM GRAPH SOLVER (QAOA)----------------- + # graph = Graph(weighted_graph, starting_node, ending_node) + # if graph.number_of_edges > 17: + # raise Exception("RAPGraph: max number of points exceeded") + + + # # Construct Hamiltonian when qubits are set as edges, + # # then optimize with QAOA/scipy: + + # hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in alphas] + + # # print(hamiltonians[0].total_hamiltonian.simplify()) + + # print("\n Calculating qubits as edges......................") + # # Run the multiprocess QAOA solver for edge qubits + # # and return the optimal path as a list of coordinates. + # line = multiprocess_qaoa_solver_edge_rap( + # hamiltonians, + # reps, + # number_processors, + # graph.number_of_edges, + # optimizer + # ) + + #---------------CLASSIC DIJKSTRA SOLVER----------------- + line = dijkstra_stepwise( + -weighted_graph, + starting_node, + ending_node + ) + return line \ No newline at end of file diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index dd6bede..a575460 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -60,6 +60,17 @@ def _build_arg_parser(): choices=["sagittal", "coronal", "axial"], help="Axis along which a slice is taken. Ignored when slice_index is None.", ) + p.add_argument( + '--sh_order', + type=int, + default=12, + help='Maximum SH order. [%(default)s]' + ) + p.add_argument( + "--slice_index", + type=int, + help="If None, a 3D graph is built.", + ) p.add_argument( "--save_only", action="store_true", @@ -68,6 +79,7 @@ def _build_arg_parser(): return p +#main used to test from command instructions def main(): parser = _build_arg_parser() args = parser.parse_args() @@ -124,8 +136,6 @@ def main(): # save output save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) - print("Graph saved") - if __name__ == "__main__": main() diff --git a/scripts/quac_optimal_path_find_max_intensity_diffusion.py b/scripts/quac_optimal_path_find_max_intensity_diffusion.py index 3c6444f..bcf02a7 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -13,7 +13,6 @@ parameters found if only one reps, and returns the optimization results. """ - def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawTextHelpFormatter @@ -74,7 +73,6 @@ def _build_arg_parser(): return p - def main(): parser = _build_arg_parser() args = parser.parse_args() diff --git a/setup.cfg b/setup.cfg index d9da413..78aeab3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -4,16 +4,22 @@ version = 0.0.1 [options] packages = - quactography + quactography.solver + quactography.hamiltonian + quactography.graph + quactography.image + quactography.visu + quactography.adj_matrix scripts data + install_requires = - dipy==1.10.* - matplotlib==3.6.* + dipy==1.11.* + matplotlib==3.10.* networkx==3.4.* - nibabel==5.2.* - numpy==1.25.* - scipy==1.11.* + nibabel==5.3.* + numpy==1.26.* + scipy==1.15.* qiskit==1.1.1 pandas == 2.2.3 seaborn == 0.13.2 @@ -31,6 +37,5 @@ console_scripts = quac_randmatrix_adj_draw.py = scripts.quac_randmatrix_adj_draw:main quac_heatmap_rep_alpha.py = scripts.quac_heatmap_rep_alpha:main quac_params_comparison.py = scripts.quac_params_comparison:main - - +