From a44f6c75d3353abf0d36dfba009eab9486f5c4e0 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Thu, 3 Jul 2025 11:01:22 -0400 Subject: [PATCH 01/19] V0.01 --- scripts/quac_matrix_adj_build.py | 12 +++++++++ ...timal_path_find_max_intensity_diffusion.py | 27 +++++++++++++++++++ 2 files changed, 39 insertions(+) diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 0a2da80..0c88ea3 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -4,6 +4,7 @@ import nibabel as nib import numpy as np import matplotlib.pyplot as plt +from pathlib import Path from quactography.adj_matrix.reconst import ( build_adjacency_matrix, @@ -15,6 +16,7 @@ ) from quactography.image.utils import slice_along_axis from quactography.adj_matrix.io import save_graph +from quac_optimal_path_find_max_intensity_diffusion import rap_funct """ @@ -120,6 +122,16 @@ def main(): # print("node indices", node_indices) # save output save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) + path = Path(args.out_graph) + + rap_funct( + path, + starting_node=0, + ending_node=len(node_indices) - 1, + output_file="rap_output", + plt_cost_landscape=False, + save_only=True + ) print("Graph saved") diff --git a/scripts/quac_optimal_path_find_max_intensity_diffusion.py b/scripts/quac_optimal_path_find_max_intensity_diffusion.py index 3c6444f..a644079 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -74,6 +74,33 @@ def _build_arg_parser(): return p +def rap_funct(graph, starting_node, ending_node, output_file, save_only, output_directory="data/output_graphs", alphas=[1.2], + reps=[1], number_processors=2, optimizer="Differential", plt_cost_landscape=False): + """ + Function to return the path from starting node to ending node in the graph. + """ + weighted_graph, _, _ = load_graph(graph) + graph = Graph(weighted_graph, starting_node, ending_node) + + # 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......................") + for i in range(len(reps)): + multiprocess_qaoa_solver_edge( + hamiltonians, + reps[i], + number_processors, + output_file, + output_directory, + optimizer, + plt_cost_landscape, + save_only, + ) def main(): parser = _build_arg_parser() From d9d1ab14e42ee576ea95ebe2fe68324a6df9ade5 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Thu, 3 Jul 2025 13:58:34 -0400 Subject: [PATCH 02/19] V0.01 of receving masks and output streamline 2D --- quactography/adj_matrix/reconst.py | 2 +- .../hamiltonian/hamiltonian_qubit_edge.py | 34 ++-- quactography/solver/qaoa_solver_qu_edge.py | 39 ++-- scripts/quac_matrix_adj_build.py | 12 +- ...timal_path_find_max_intensity_diffusion.py | 166 +++++++++--------- 5 files changed, 134 insertions(+), 119 deletions(-) diff --git a/quactography/adj_matrix/reconst.py b/quactography/adj_matrix/reconst.py index ad6c4a4..dc6492a 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -98,7 +98,7 @@ def build_weighted_graph(adj_matrix, node_indices, sh, axis_name): """ # Get directions depending if we are in axial, coronal or sagittal : sphere = _get_sphere_for_axis(axis_name) - sf = sh_to_sf(sh, sphere, sh_order=12, basis_type="tournier07") + sf = sh_to_sf(sh, sphere, sh_order=8, basis_type="tournier07") sf[sf < 0.0] = 0.0 sf /= np.max(sf, axis=(0, 1), keepdims=True) # sf = sf / np.max(sf, axis=(-1), keepdims=True) diff --git a/quactography/hamiltonian/hamiltonian_qubit_edge.py b/quactography/hamiltonian/hamiltonian_qubit_edge.py index 9d27356..4597a47 100644 --- a/quactography/hamiltonian/hamiltonian_qubit_edge.py +++ b/quactography/hamiltonian/hamiltonian_qubit_edge.py @@ -54,7 +54,7 @@ def mandatory_cost(self): for i in range(self.graph.number_of_edges): str1 = ( "I" * (self.graph.number_of_edges - i - 1) + "Z" + "I" * i, - -self.graph.weights[0][i] / 2, + -self.graph.weights[i] / 2, ) pauli_weight_first_term.append(str1) @@ -272,22 +272,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( @@ -351,14 +349,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/qaoa_solver_qu_edge.py b/quactography/solver/qaoa_solver_qu_edge.py index ba6ed51..ca1b3b4 100644 --- a/quactography/solver/qaoa_solver_qu_edge.py +++ b/quactography/solver/qaoa_solver_qu_edge.py @@ -215,6 +215,7 @@ def multiprocess_qaoa_solver_edge( optimizer, cost_landscape, save_only, + nb_edges ): """ Solve the optimization problem using the QAOA algorithm @@ -244,19 +245,31 @@ def multiprocess_qaoa_solver_edge( ------- None """ - pool = multiprocessing.Pool(nbr_processes) - pool.map( - find_longest_path, - zip( - hamiltonians, - itertools.repeat(reps), - itertools.repeat(output_file), - itertools.repeat(optimizer), - itertools.repeat(cost_landscape), - itertools.repeat(output_folder), - itertools.repeat(save_only), - ), - ) + # pool = multiprocessing.Pool(nbr_processes) + # pool.map( + # find_longest_path, + # zip( + # hamiltonians, + # itertools.repeat(reps), + # itertools.repeat(output_file), + # itertools.repeat(optimizer), + # itertools.repeat(cost_landscape), + # itertools.repeat(output_folder), + # itertools.repeat(save_only), + # ), + # ) + exact_path = ((hamiltonians[0].exact_path[0][::-1]).zfill(nb_edges)) + print(f"Exact path: {exact_path}... length: {list(exact_path)}") + all_coords = [] + value_to_find = '1' + + for i, x in enumerate(exact_path): + if x == value_to_find: + all_coords.append(np.unravel_index(i, (hamiltonians[0].graph.number_of_edges, hamiltonians[0].graph.number_of_edges))) + + + print(f"All occurrences of {value_to_find} are at indices: {all_coords}") + print( "------------------MULTIPROCESS SOLVER FINISHED-------------------------" diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 0c88ea3..0480d35 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -16,7 +16,7 @@ ) from quactography.image.utils import slice_along_axis from quactography.adj_matrix.io import save_graph -from quac_optimal_path_find_max_intensity_diffusion import rap_funct +from scripts.quac_optimal_path_find_max_intensity_diffusion import rap_funct """ @@ -65,6 +65,7 @@ def _build_arg_parser(): def main(): + parser = _build_arg_parser() args = parser.parse_args() nodes_mask_im = nib.load(args.in_nodes_mask) @@ -76,7 +77,7 @@ def main(): keep_node_indices = None if args.keep_mask: - keep_mask_im = nib.load(args.keep_mask) + keep_mask_im = nib.load(keep_mask) keep_mask = slice_along_axis( keep_mask_im.get_fdata().astype(bool), args.axis_name, args.slice_index ) @@ -122,17 +123,18 @@ def main(): # print("node indices", node_indices) # save output save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) + print("Graph saved") path = Path(args.out_graph) rap_funct( path, - starting_node=0, - ending_node=len(node_indices) - 1, + starting_node=node_indices[0], + ending_node=node_indices[3], output_file="rap_output", plt_cost_landscape=False, save_only=True ) - print("Graph saved") + if __name__ == "__main__": diff --git a/scripts/quac_optimal_path_find_max_intensity_diffusion.py b/scripts/quac_optimal_path_find_max_intensity_diffusion.py index a644079..06ecf2d 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -14,65 +14,65 @@ """ -def _build_arg_parser(): - p = argparse.ArgumentParser( - description=__doc__, formatter_class=argparse.RawTextHelpFormatter - ) - p.add_argument( - "in_graph", - help="Adjacency matrix which graph we want path that maximizes weights in graph, (npz file)", - type=str, - ) - p.add_argument("starting_node", - help="Starting node of the graph", type=int) - p.add_argument("ending_node", - help="Ending node of the graph", type=int) - p.add_argument("output_file", - help="Output file name (npz file)", type=str) - p.add_argument("output_directory", - help="Directory where the files will be outputed", type=str, - default="data/output_graphs/" - ) - p.add_argument( - "--alphas", - nargs="+", - type=float, - help="List of alphas", - default=[1.2] - ) - p.add_argument( - "--reps", - nargs="+", - type=int, - help="List of repetitions to run for the QAOA algorithm", - default=[1], - ) - p.add_argument( - "-npr", - "--number_processors", - help="Number of cpu to use for multiprocessing", - default=1, - type=int, - ) - p.add_argument( - "--optimizer", - help="Optimizer to use for the QAOA algorithm", - default="Differential", - type=str, - ) - p.add_argument( - "--plt_cost_landscape", - help="True or False, Plot 3D and 2D of the cost landscape" - "(for gamma and beta compact set over all possible angles-0.1 incrementation)", - action="store_false", - ) - p.add_argument( - "--save_only", - help="Save only the figure without displaying it", - action="store_true", - ) - - return p +# def _build_arg_parser(): +# p = argparse.ArgumentParser( +# description=__doc__, formatter_class=argparse.RawTextHelpFormatter +# ) +# p.add_argument( +# "in_graph", +# help="Adjacency matrix which graph we want path that maximizes weights in graph, (npz file)", +# type=str, +# ) +# p.add_argument("starting_node", +# help="Starting node of the graph", type=int) +# p.add_argument("ending_node", +# help="Ending node of the graph", type=int) +# p.add_argument("output_file", +# help="Output file name (npz file)", type=str) +# p.add_argument("output_directory", +# help="Directory where the files will be outputed", type=str, +# default="data/output_graphs/" +# ) +# p.add_argument( +# "--alphas", +# nargs="+", +# type=float, +# help="List of alphas", +# default=[1.2] +# ) +# p.add_argument( +# "--reps", +# nargs="+", +# type=int, +# help="List of repetitions to run for the QAOA algorithm", +# default=[1], +# ) +# p.add_argument( +# "-npr", +# "--number_processors", +# help="Number of cpu to use for multiprocessing", +# default=1, +# type=int, +# ) +# p.add_argument( +# "--optimizer", +# help="Optimizer to use for the QAOA algorithm", +# default="Differential", +# type=str, +# ) +# p.add_argument( +# "--plt_cost_landscape", +# help="True or False, Plot 3D and 2D of the cost landscape" +# "(for gamma and beta compact set over all possible angles-0.1 incrementation)", +# action="store_false", +# ) +# p.add_argument( +# "--save_only", +# help="Save only the figure without displaying it", +# action="store_true", +# ) + +# return p def rap_funct(graph, starting_node, ending_node, output_file, save_only, output_directory="data/output_graphs", alphas=[1.2], reps=[1], number_processors=2, optimizer="Differential", plt_cost_landscape=False): @@ -86,6 +86,7 @@ def rap_funct(graph, starting_node, ending_node, output_file, save_only, output_ # then optimize with QAOA/scipy: hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in alphas] + print(hamiltonians[0]) # print(hamiltonians[0].total_hamiltonian.simplify()) @@ -100,35 +101,36 @@ def rap_funct(graph, starting_node, ending_node, output_file, save_only, output_ optimizer, plt_cost_landscape, save_only, + graph.number_of_edges ) -def main(): - parser = _build_arg_parser() - args = parser.parse_args() +# def main(): +# parser = _build_arg_parser() +# args = parser.parse_args() - weighted_graph, _, _ = load_graph(args.in_graph) - graph = Graph(weighted_graph, args.starting_node, args.ending_node) +# weighted_graph, _, _ = load_graph(args.in_graph) +# graph = Graph(weighted_graph, args.starting_node, args.ending_node) - # Construct Hamiltonian when qubits are set as edges, - # then optimize with QAOA/scipy: +# # Construct Hamiltonian when qubits are set as edges, +# # then optimize with QAOA/scipy: - hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in args.alphas] +# hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in args.alphas] - # print(hamiltonians[0].total_hamiltonian.simplify()) +# # print(hamiltonians[0].total_hamiltonian.simplify()) - print("\n Calculating qubits as edges......................") - for i in range(len(args.reps)): - multiprocess_qaoa_solver_edge( - hamiltonians, - args.reps[i], - args.number_processors, - args.output_file, - args.output_directory, - args.optimizer, - args.plt_cost_landscape, - args.save_only, - ) +# print("\n Calculating qubits as edges......................") +# for i in range(len(args.reps)): +# multiprocess_qaoa_solver_edge( +# hamiltonians, +# args.reps[i], +# args.number_processors, +# args.output_file, +# args.output_directory, +# args.optimizer, +# args.plt_cost_landscape, +# args.save_only, +# ) -if __name__ == "__main__": - main() +# if __name__ == "__main__": +# main() From 2a6c8001a1f3583d06dadc6bab50d64d7413e30a Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Thu, 3 Jul 2025 15:07:42 -0400 Subject: [PATCH 03/19] command script to function --- scripts/quac_matrix_adj_build.py | 115 ++++++++++++++++--------------- 1 file changed, 59 insertions(+), 56 deletions(-) diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 0480d35..e18eec0 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -24,74 +24,75 @@ """ -def _build_arg_parser(): - p = argparse.ArgumentParser( - description=__doc__, formatter_class=argparse.RawTextHelpFormatter - ) - p.add_argument("in_nodes_mask", - help="Input nodes mask image (.nii.gz file)") - p.add_argument("in_sh", - help="Input SH image. (nii.gz file)") - p.add_argument("out_graph", - help="Output graph file name (npz file)") - - p.add_argument( - "--keep_mask", - help="Nodes that must not be filtered out." - ) - p.add_argument( - "--threshold", - default=0.0, - type=float, - help="Cut all weights below a given threshold. [%(default)s]", - ) - p.add_argument( - "--slice_index", - type=int, - help="If None, midslice is taken." - ) - p.add_argument( - "--axis_name", - default="axial", - choices=["sagittal", "coronal", "axial"], - help="Axis along which a slice is taken.", - ) - p.add_argument( - "--save_only", - action="store_true", - help="Does not plot the matrix, only saves a copy of npz adjacency matrix file" - ) - return p - - -def main(): +# def _build_arg_parser(): +# p = argparse.ArgumentParser( +# description=__doc__, formatter_class=argparse.RawTextHelpFormatter +# ) +# p.add_argument("in_nodes_mask", +# help="Input nodes mask image (.nii.gz file)") +# p.add_argument("in_sh", +# help="Input SH image. (nii.gz file)") +# p.add_argument("out_graph", +# help="Output graph file name (npz file)") + +# p.add_argument( +# "--keep_mask", +# help="Nodes that must not be filtered out." +# ) +# p.add_argument( +# "--threshold", +# default=0.0, +# type=float, +# help="Cut all weights below a given threshold. [%(default)s]", +# ) +# p.add_argument( +# "--slice_index", +# type=int, +# help="If None, midslice is taken." +# ) +# p.add_argument( +# "--axis_name", +# default="axial", +# choices=["sagittal", "coronal", "axial"], +# help="Axis along which a slice is taken.", +# ) +# p.add_argument( +# "--save_only", +# action="store_true", +# help="Does not plot the matrix, only saves a copy of npz adjacency matrix file" +# ) +# return p + + +def build_mat(in_nodes_mask, in_sh, out_graph, + keep_mask=None, threshold=0.2, slice_index=None, + axis_name="axial", save_only=False): + - parser = _build_arg_parser() - args = parser.parse_args() - nodes_mask_im = nib.load(args.in_nodes_mask) - sh_im = nib.load(args.in_sh) + nodes_mask_im = nib.load(in_nodes_mask) + sh_im = nib.load(in_sh) nodes_mask = slice_along_axis( - nodes_mask_im.get_fdata().astype(bool), args.axis_name, args.slice_index + nodes_mask_im.get_fdata().astype(bool), axis_name, slice_index ) keep_node_indices = None - if args.keep_mask: + if keep_mask: keep_mask_im = nib.load(keep_mask) keep_mask = slice_along_axis( - keep_mask_im.get_fdata().astype(bool), args.axis_name, args.slice_index + keep_mask_im.get_fdata().astype(bool), axis_name, slice_index ) keep_node_indices = np.flatnonzero(keep_mask) # !!Careful, we remove a dimension, but the SH amplitudes still exist in 3D - sh = slice_along_axis(sh_im.get_fdata(), args.axis_name, args.slice_index) + sh = slice_along_axis(sh_im.get_fdata(), axis_name, slice_index) # adjacency graph adj_matrix, node_indices = build_adjacency_matrix(nodes_mask) # assign edge weights weighted_graph, node_indices = build_weighted_graph( - adj_matrix, node_indices, sh, args.axis_name + adj_matrix, node_indices, sh, axis_name ) # # Could be added in the code if needed: @@ -104,7 +105,7 @@ def main(): # # _______________________________________________________________________ # filter graph edges by weight - weighted_graph[weighted_graph < args.threshold] = 0.0 + weighted_graph[weighted_graph < threshold] = 0.0 # remove intermediate nodes that connect only two nodes weighted_graph = remove_intermediate_connections( @@ -116,15 +117,15 @@ def main(): weighted_graph, node_indices, keep_node_indices ) - if not args.save_only: + if not save_only: plt.imshow(np.log(weighted_graph + 1)) plt.show() # print("node indices", node_indices) # save output - save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) + save_graph(weighted_graph, node_indices, nodes_mask.shape, out_graph) print("Graph saved") - path = Path(args.out_graph) + path = Path(out_graph) rap_funct( path, @@ -137,5 +138,7 @@ def main(): -if __name__ == "__main__": - main() +# if __name__ == "__main__": +# parser = _build_arg_parser() +# args = parser.parse_args() +# main(args.in_nodes_mask, args.in_sh, args.out_graph) From e595e4128ca7b1c7e5e64e407279e45ad85b1a34 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Thu, 3 Jul 2025 16:05:50 -0400 Subject: [PATCH 04/19] copy of 3D pull request --- quactography/adj_matrix/filter.py | 26 +++++- quactography/adj_matrix/reconst.py | 98 ++++++++++------------ quactography/solver/qaoa_solver_qu_edge.py | 8 +- scripts/quac_matrix_adj_build.py | 37 +++++--- 4 files changed, 101 insertions(+), 68 deletions(-) diff --git a/quactography/adj_matrix/filter.py b/quactography/adj_matrix/filter.py index 2cbbf53..a8264f1 100644 --- a/quactography/adj_matrix/filter.py +++ b/quactography/adj_matrix/filter.py @@ -34,6 +34,30 @@ def remove_orphan_nodes(graph, node_indices, keep_indices=None): return out_graph, out_indices +def extract_slice_at_index(graph, node_indices, shape, slice_index, axis_name): + out_graph = [] + out_it = [] + for it, graph_row in enumerate(graph): + if _is_on_target_slice(node_indices[it], shape, slice_index, axis_name): + out_graph.append(graph_row) + out_it.append(it) + out_graph = np.take(np.asarray(out_graph), out_it, axis=1) + out_indices = node_indices[np.asarray(out_it)] + + return out_graph, out_indices + + +def _is_on_target_slice(node_idx, shape, slice_index, axis_name): + x, y, z = np.unravel_index(node_idx, shape) + if axis_name == 'sagittal': + return x == slice_index + if axis_name == 'coronal': + return y == slice_index + if axis_name == 'axial': + return z == slice_index + raise ValueError((f'Unknown axis name {axis_name}.')) + + def remove_intermediate_connections(graphR, node_indices=None, keep_indices=None): """ Remove nodes that do not add a change in direction between the two nodes @@ -136,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 dc6492a..602f16c 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -30,8 +30,6 @@ def build_adjacency_matrix(nodes_mask): adj_matrix : np.ndarray Adjacency matrix of the nodes. """ - assert len(nodes_mask.shape) == 2 - # 1st. Assign labels to non-zero voxels (nodes) node_indices = np.flatnonzero(nodes_mask) @@ -50,32 +48,23 @@ def build_adjacency_matrix(nodes_mask): label = node_indices[i] # Coordinates of voxel - x, y = np.unravel_index(label, nodes_mask.shape) - - # Adds possibility of an edge to the actual node closest neighbour in 8 directions - for offset in [-1, 1]: - adj_matrix = _add_edge_perhaps( - x + offset, y, i, nodes_mask, labels_volume, adj_matrix - ) - adj_matrix = _add_edge_perhaps( - x, y + offset, i, nodes_mask, labels_volume, adj_matrix - ) - for y_offset in [-1, 1]: - adj_matrix = _add_edge_perhaps( - x + offset, y + y_offset, i, nodes_mask, labels_volume, adj_matrix - ) - # # Verify number of non-zero potential connexions (max) for closest neighbours : - # num_non_zero = 0 - # adj_matrix_1 = adj_matrix.flatten() - # for i in adj_matrix_1: - # if i != 0: - # num_non_zero += 1 - # print(num_non_zero) + x, y, z = np.unravel_index(label, nodes_mask.shape) + + # 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 -def build_weighted_graph(adj_matrix, node_indices, sh, axis_name): +def build_weighted_graph(adj_matrix, node_indices, sh, sh_order=12): """ Build the weighted graph for a given set of nodes in a 2D image. @@ -86,9 +75,9 @@ def build_weighted_graph(adj_matrix, node_indices, sh, axis_name): node_indices : np.ndarray Indices of the nodes in the graph. sh : np.ndarray - Spherical harmonics coefficients, set to 12th order here, needs to be changed if needed. - axis_name : str - Name of the axis, can be "sagittal", "coronal" or "axial". + Spherical harmonics coefficients. + sh_order: int, optional + Spherical harmonics maximum order. Returns ------- weighted_graph : np.ndarray @@ -97,8 +86,8 @@ def build_weighted_graph(adj_matrix, node_indices, sh, axis_name): Indices of the nodes in the graph. """ # Get directions depending if we are in axial, coronal or sagittal : - sphere = _get_sphere_for_axis(axis_name) - sf = sh_to_sf(sh, sphere, sh_order=8, basis_type="tournier07") + sphere = _get_sphere() + sf = sh_to_sf(sh, sphere, sh_order=sh_order, basis_type="tournier07") sf[sf < 0.0] = 0.0 sf /= np.max(sf, axis=(0, 1), keepdims=True) # sf = sf / np.max(sf, axis=(-1), keepdims=True) @@ -106,28 +95,28 @@ def build_weighted_graph(adj_matrix, node_indices, sh, axis_name): # print(sh.shape) weighted_graph = np.zeros_like(adj_matrix) - x, y = np.unravel_index(node_indices, sh.shape[:2]) + x, y, z = np.unravel_index(node_indices, sh.shape[:3]) # node traversal for it, node_row in enumerate(adj_matrix): nb_connections = np.count_nonzero(node_row) if nb_connections > 0: - start_x, start_y = x[it], y[it] + start_x, start_y, start_z = x[it], y[it], z[it] # which nodes are connected to every starting node: - connected_xs, connected_ys = x[node_row > 0], y[node_row > 0] + connected_xs, connected_ys, connected_zs = x[node_row > 0], y[node_row > 0], z[node_row > 0] w_list = [] for conn_idx in range(nb_connections): - conn_x, conn_y = connected_xs[conn_idx], connected_ys[conn_idx] + conn_x, conn_y, conn_z = connected_xs[conn_idx], connected_ys[conn_idx], connected_zs[conn_idx] - direction = np.array([[conn_x, conn_y]], dtype=float) - np.array( - [[start_x, start_y]], dtype=float + direction = np.array([[conn_x, conn_y, conn_z]], dtype=float) - np.array( + [[start_x, start_y, start_z]], dtype=float ) # The directions : - dir_id = np.argmax(np.dot(direction, DIRECTIONS_2D.T)) + dir_id = np.argmax(np.dot(direction, sphere.vertices.T)) - w = sf[start_x, start_y, dir_id] + sf[conn_x, conn_y, dir_id] + w = sf[start_x, start_y, start_z, dir_id] + sf[conn_x, conn_y, conn_z, dir_id] w_list.append(w) weights = np.zeros((len(node_row),)) @@ -137,33 +126,34 @@ def build_weighted_graph(adj_matrix, node_indices, sh, axis_name): return weighted_graph, node_indices -def _get_sphere_for_axis(axis_name): - directions = np.zeros((len(DIRECTIONS_2D), 3)) - if axis_name == "sagittal": - directions[:, 1] = DIRECTIONS_2D[:, 0] - directions[:, 2] = DIRECTIONS_2D[:, 1] - elif axis_name == "coronal": - directions[:, 0] = DIRECTIONS_2D[:, 0] - directions[:, 2] = DIRECTIONS_2D[:, 1] - elif axis_name == "axial": - directions[:, 0] = DIRECTIONS_2D[:, 0] - directions[:, 1] = DIRECTIONS_2D[:, 1] - # print(directions) +def _get_sphere(): + directions = [] + for i in [-1, 0, 1]: + for j in [-1, 0, 1]: + for k in [-1, 0, 1]: + if i == 0 and j == 0 and k == 0: + continue + vec = np.array([i, j, k], dtype=float) + vec /= np.linalg.norm(vec) + directions.append(vec) + directions = np.asarray(directions) return Sphere(xyz=directions) def _add_edge_perhaps( - pos_x, pos_y, current_node, nodes_mask, labels_volume, adj_matrix + pos_x, pos_y, pos_z, current_node, nodes_mask, labels_volume, adj_matrix ): - if _is_valid_pos(pos_x, pos_y, nodes_mask): - neighbor_label = labels_volume[pos_x, pos_y] + if _is_valid_pos(pos_x, pos_y, pos_z, nodes_mask): + neighbor_label = labels_volume[pos_x, pos_y, pos_z] adj_matrix[current_node, neighbor_label] = 1 return adj_matrix -def _is_valid_pos(pos_x, pos_y, nodes_mask): +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 if pos_y < 0 or pos_y >= nodes_mask.shape[1]: return False - return nodes_mask[pos_x, pos_y] + if pos_z < 0 or pos_z >= nodes_mask.shape[2]: + return False + return nodes_mask[pos_x, pos_y, pos_z] \ No newline at end of file diff --git a/quactography/solver/qaoa_solver_qu_edge.py b/quactography/solver/qaoa_solver_qu_edge.py index ca1b3b4..aae4471 100644 --- a/quactography/solver/qaoa_solver_qu_edge.py +++ b/quactography/solver/qaoa_solver_qu_edge.py @@ -267,7 +267,13 @@ def multiprocess_qaoa_solver_edge( if x == value_to_find: all_coords.append(np.unravel_index(i, (hamiltonians[0].graph.number_of_edges, hamiltonians[0].graph.number_of_edges))) - + if os.path.exists(output_folder) == False: + os.makedirs(output_folder) + + np.savez( + "streamline_coords", + all_coords=all_coords, + ) print(f"All occurrences of {value_to_find} are at indices: {all_coords}") diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index e18eec0..10de0aa 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -13,6 +13,7 @@ from quactography.adj_matrix.filter import ( remove_orphan_nodes, remove_intermediate_connections, + extract_slice_at_index ) from quactography.image.utils import slice_along_axis from quactography.adj_matrix.io import save_graph @@ -56,6 +57,18 @@ # choices=["sagittal", "coronal", "axial"], # help="Axis along which a slice is taken.", # ) + # 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, midslice is taken." + # help="If None, a 3D graph is built." + # ) # p.add_argument( # "--save_only", # action="store_true", @@ -66,33 +79,28 @@ def build_mat(in_nodes_mask, in_sh, out_graph, keep_mask=None, threshold=0.2, slice_index=None, - axis_name="axial", save_only=False): + axis_name="axial", sh_order=8, save_only=False): nodes_mask_im = nib.load(in_nodes_mask) sh_im = nib.load(in_sh) - nodes_mask = slice_along_axis( - nodes_mask_im.get_fdata().astype(bool), axis_name, slice_index - ) + nodes_mask = nodes_mask_im.get_fdata().astype(bool) + keep_node_indices = None if keep_mask: - keep_mask_im = nib.load(keep_mask) - keep_mask = slice_along_axis( - keep_mask_im.get_fdata().astype(bool), axis_name, slice_index - ) + keep_mask = nib.load(keep_mask).get_fdata().astype(bool) keep_node_indices = np.flatnonzero(keep_mask) - # !!Careful, we remove a dimension, but the SH amplitudes still exist in 3D - sh = slice_along_axis(sh_im.get_fdata(), axis_name, slice_index) + sh = sh_im.get_fdata() # adjacency graph adj_matrix, node_indices = build_adjacency_matrix(nodes_mask) # assign edge weights weighted_graph, node_indices = build_weighted_graph( - adj_matrix, node_indices, sh, axis_name + adj_matrix, node_indices, sh, sh_order ) # # Could be added in the code if needed: @@ -116,6 +124,10 @@ def build_mat(in_nodes_mask, in_sh, out_graph, 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 + ) if not save_only: plt.imshow(np.log(weighted_graph + 1)) @@ -135,7 +147,8 @@ def build_mat(in_nodes_mask, in_sh, out_graph, plt_cost_landscape=False, save_only=True ) - + + # if __name__ == "__main__": From 488baeba49cb484ef7963eb602839e786593806a Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 8 Jul 2025 09:52:30 -0400 Subject: [PATCH 05/19] RAP WIP v0.02 --- data/hurdle_mask.nii.gz | Bin 0 -> 1665 bytes quactography/graph/undirected_graph.py | 6 +- .../hamiltonian/hamiltonian_qubit_edge.py | 2 +- scripts/quac_matrix_adj_build.py | 166 ++++++++++++------ 4 files changed, 115 insertions(+), 59 deletions(-) create mode 100644 data/hurdle_mask.nii.gz diff --git a/data/hurdle_mask.nii.gz b/data/hurdle_mask.nii.gz new file mode 100644 index 0000000000000000000000000000000000000000..6d9dd2b3395bf86b3f4708280a1e835f603d007d GIT binary patch literal 1665 zcmd6n{ZrBh9LLSmop!gaG)uR#WvQ2*=hD;?Q`_xC=cW}YK7_VRy+CyOO%YMPd3(64 z6i-bGQ?Xf-h4@{f7+TOy6GMfW0RcXAGhoS6`4totowmE&^%wN{>GR|3{dsupz3ODQ zQ4VW*`FO1|blXQ~n5iPdu$16jD19!*<3xEpe9^ zh(8Hn#hH?0WbTStubp+t9Jn3GCu^59;t|y}ALI{d2XuQ&tSV81GJHrd%^Oh}XV%<( z5x&uv?6%lDIBoml!)WgBe_f~v3EaMk^BQa4gIG~?kEmwXpgtr%4TGuM#C=d7t@7mC z=(h~xU4C#{wzvh|?tA%Nq({U`v^)$g^#JcN;`8DppyLH%v71yzSlyHB^EP*BV5LF~ z|G^aM`m4!@2?VcpzX0;U=Li5-M(DHfd$@u||}23-fR^VY?b?=q5~ z1MJ*RmWa8=6i-?eIz?q~V zJrG$qMLinkm%)x|Bl18V9p0)@So(&v=EXvM30@Jrx|VgVDGCa`3B}rPU`8ajT7AKa z44u$_M;e|hfEWM`jmpFqa|iCG;L&1^>D zJ&daC$FM*WB0?|tIVTyZ=u{}&d41(CbVgY|cy5xc0`xR)R@umB@ZT$wk>8IBwaxK3 zwt$w3@!5-)xp5Yx_4@Jdd%WQ%n%N6k8Ob}7TAtvXwE{~|Cx8_V&3T#1j4l58xxY+t zz(mE79Cc&|d#=ta%rKXAO4dM+Vs;wJy>nbPwpvL?x9W+VdQzuJkGdZTmlGwlB1bk(c|!B=l~@c zK$ImWj*4n&U^$ESm!i;I%G*s^IN~>%kp(9BysX;B=r~lLtDF-CLHVpUmF@~a#{5iU zc)vaO8C0vaz9Za7eO#*@L=L?U9$aLGk=t@F_^u~@CTqMc8Bnle_}f<{c3-RQJtKxzd{ zV`HK#`8|?Xc$I|Ki=CO^^IS@OsEnbm5c6Fvjw<8%edyw*^GU%^27c^8tg2^%w9(s zcP%o@Z#U7k=~8xtVV1HR@5(v04Ckoxp{|5k%KOJ5sn%m+`9>S^)^)aMs=CQW)~7)^HpB7iwjFPM^F!hb zc~n6g@A-vi&l*$lyHeqQmd6_hM{mvq3#>tb3QI@65Kw!pZCIaT4Mi(sCnd3ajmZVi{FT<^PtHbc4wf#zTq+j39Qd@^I5w zstCdugDl$x#WP?HBSDRDD*coq`*-yB#Sk5wbZ(5* y-%-uIggrz5Gtzseyz96!)s}YJ^Fz|;ras?;-8APr<6o>?v#YFQrkyoyvHKS%fAfO? literal 0 HcmV?d00001 diff --git a/quactography/graph/undirected_graph.py b/quactography/graph/undirected_graph.py index 9ed3a19..ef0621f 100644 --- a/quactography/graph/undirected_graph.py +++ b/quactography/graph/undirected_graph.py @@ -81,9 +81,9 @@ def connexions_edges(self, mat_adj, num_nodes): # print(f"End :{ending_nodes}") weights = [] - for _ in range(number_of_edges): - for _ in starting_nodes: - for _ in ending_nodes: + for _ in range(10): + for _ in range(10): + for _ in range(10): weight_qubit = mat_adj[starting_nodes, ending_nodes] weights.append(weight_qubit) # weights[0][0] diff --git a/quactography/hamiltonian/hamiltonian_qubit_edge.py b/quactography/hamiltonian/hamiltonian_qubit_edge.py index 4597a47..24ab28b 100644 --- a/quactography/hamiltonian/hamiltonian_qubit_edge.py +++ b/quactography/hamiltonian/hamiltonian_qubit_edge.py @@ -54,7 +54,7 @@ def mandatory_cost(self): for i in range(self.graph.number_of_edges): str1 = ( "I" * (self.graph.number_of_edges - i - 1) + "Z" + "I" * i, - -self.graph.weights[i] / 2, + -self.graph.weights[0][i] / 2, ) pauli_weight_first_term.append(str1) diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 10de0aa..63a3887 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -15,7 +15,6 @@ remove_intermediate_connections, extract_slice_at_index ) -from quactography.image.utils import slice_along_axis from quactography.adj_matrix.io import save_graph from scripts.quac_optimal_path_find_max_intensity_diffusion import rap_funct @@ -25,56 +24,50 @@ """ -# def _build_arg_parser(): -# p = argparse.ArgumentParser( -# description=__doc__, formatter_class=argparse.RawTextHelpFormatter -# ) -# p.add_argument("in_nodes_mask", -# help="Input nodes mask image (.nii.gz file)") -# p.add_argument("in_sh", -# help="Input SH image. (nii.gz file)") -# p.add_argument("out_graph", -# help="Output graph file name (npz file)") - -# p.add_argument( -# "--keep_mask", -# help="Nodes that must not be filtered out." -# ) -# p.add_argument( -# "--threshold", -# default=0.0, -# type=float, -# help="Cut all weights below a given threshold. [%(default)s]", -# ) -# p.add_argument( -# "--slice_index", -# type=int, -# help="If None, midslice is taken." -# ) -# p.add_argument( -# "--axis_name", -# default="axial", -# choices=["sagittal", "coronal", "axial"], -# help="Axis along which a slice is taken.", -# ) - # 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, midslice is taken." - # help="If None, a 3D graph is built." - # ) -# p.add_argument( -# "--save_only", -# action="store_true", -# help="Does not plot the matrix, only saves a copy of npz adjacency matrix file" -# ) -# return p +def _build_arg_parser(): + p = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawTextHelpFormatter + ) + p.add_argument("in_nodes_mask", + help="Input nodes mask image (.nii.gz file)") + p.add_argument("in_sh", + help="Input SH image. (nii.gz file)") + p.add_argument("out_graph", + help="Output graph file name (npz file)") + + p.add_argument( + "--keep_mask", + help="Nodes that must not be filtered out." + ) + p.add_argument( + "--threshold", + default=0.0, + type=float, + help="Cut all weights below a given threshold. [%(default)s]", + ) + p.add_argument( + "--axis_name", + default="axial", + choices=["sagittal", "coronal", "axial"], + help="Axis along which a slice is taken.", + ) + 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", + help="Does not plot the matrix, only saves a copy of npz adjacency matrix file" + ) + return p def build_mat(in_nodes_mask, in_sh, out_graph, @@ -148,10 +141,73 @@ def build_mat(in_nodes_mask, in_sh, out_graph, save_only=True ) +def main(): + parser = _build_arg_parser() + args = parser.parse_args() + nodes_mask_im = nib.load(args.in_nodes_mask) + sh_im = nib.load(args.in_sh) + + nodes_mask = nodes_mask_im.get_fdata().astype(bool) + + keep_node_indices = None + if args.keep_mask: + keep_mask = nib.load(args.keep_mask).get_fdata().astype(bool) + keep_node_indices = np.flatnonzero(keep_mask) + + sh = sh_im.get_fdata() + + # adjacency graph + adj_matrix, node_indices = build_adjacency_matrix(nodes_mask) + + # assign edge weights + weighted_graph, node_indices = build_weighted_graph( + adj_matrix, node_indices, sh, args.sh_order + ) + + # # Could be added in the code if needed: + # # Select sub-graph and filter:___________________________________________ + # select_1 = 30 + # select_2 = 200 + + # weighted_graph = weighted_graph[select_1:select_2, select_1:select_2] + # # print(weighted_graph) + # # _______________________________________________________________________ + + # filter graph edges by weight + weighted_graph[weighted_graph < args.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 args.slice_index is not None: + weighted_graph, node_indices = extract_slice_at_index( + weighted_graph, node_indices, nodes_mask.shape, args.slice_index, args.axis_name + ) + + if not args.save_only: + plt.imshow(np.log(weighted_graph + 1)) + plt.show() + + # print("node indices", node_indices) + # save output + save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) + rap_funct( + args.out_graph, + starting_node=node_indices[0], + ending_node=node_indices[3], + output_file="rap_output", + plt_cost_landscape=False, + save_only=True + ) + print("Graph saved") -# if __name__ == "__main__": -# parser = _build_arg_parser() -# args = parser.parse_args() -# main(args.in_nodes_mask, args.in_sh, args.out_graph) +if __name__ == "__main__": + main() From 3d946017ca9dc7f7f07bb8f0fd27648e7b4032a5 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 15 Jul 2025 10:45:30 -0400 Subject: [PATCH 06/19] Quacto package installation --- setup.cfg | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/setup.cfg b/setup.cfg index d9da413..729e6d4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,10 +3,12 @@ name = quactography version = 0.0.1 [options] -packages = - quactography +packages = quactography.solver + quactography.adj_matrix + quactography.graph + quactography.hamiltonian scripts - data + install_requires = dipy==1.10.* matplotlib==3.6.* @@ -31,6 +33,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 - - + From c7e3567bfdef59ccff7ce387f85c6f932d36af07 Mon Sep 17 00:00:00 2001 From: desk3505 Date: Thu, 24 Jul 2025 11:52:00 -0400 Subject: [PATCH 07/19] Changed function name, added docs, allowed for command launch of scripts for tests if need be --- quactography/solver/qaoa_solver_qu_edge.py | 192 +++++++++++++-- scripts/quac_matrix_adj_build.py | 115 +++++---- ...timal_path_find_max_intensity_diffusion.py | 226 ++++++++++-------- 3 files changed, 360 insertions(+), 173 deletions(-) diff --git a/quactography/solver/qaoa_solver_qu_edge.py b/quactography/solver/qaoa_solver_qu_edge.py index aae4471..63e21b6 100644 --- a/quactography/solver/qaoa_solver_qu_edge.py +++ b/quactography/solver/qaoa_solver_qu_edge.py @@ -190,7 +190,27 @@ 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): + """ + 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 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 + """ path = Path(in_folder) glob_path = path.glob(f'*_alpha_{alpha}_reps_{rep}.npz') @@ -214,8 +234,7 @@ def multiprocess_qaoa_solver_edge( output_folder, optimizer, cost_landscape, - save_only, - nb_edges + save_only ): """ Solve the optimization problem using the QAOA algorithm @@ -245,17 +264,67 @@ def multiprocess_qaoa_solver_edge( ------- None """ + pool = multiprocessing.Pool(nbr_processes) + pool.map( + find_longest_path, + zip( + hamiltonians, + itertools.repeat(reps), + itertools.repeat(output_file), + itertools.repeat(optimizer), + itertools.repeat(cost_landscape), + itertools.repeat(output_folder), + itertools.repeat(save_only), + ), + ) + print( + "------------------MULTIPROCESS SOLVER FINISHED-------------------------" + ) + +#Only returns solution by classic solution as it is +def multiprocess_qaoa_solver_edge_rap( + hamiltonians, + reps, + nbr_processes, + nb_edges, + optimizer="Differential", +): + """ + 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 + ------- + line : list + List of coordinates for the streamline. + """ # pool = multiprocessing.Pool(nbr_processes) # pool.map( - # find_longest_path, + # find_longest_path_rap, # zip( # hamiltonians, # itertools.repeat(reps), - # itertools.repeat(output_file), # itertools.repeat(optimizer), - # itertools.repeat(cost_landscape), - # itertools.repeat(output_folder), - # itertools.repeat(save_only), # ), # ) exact_path = ((hamiltonians[0].exact_path[0][::-1]).zfill(nb_edges)) @@ -266,17 +335,100 @@ def multiprocess_qaoa_solver_edge( for i, x in enumerate(exact_path): if x == value_to_find: all_coords.append(np.unravel_index(i, (hamiltonians[0].graph.number_of_edges, hamiltonians[0].graph.number_of_edges))) - - if os.path.exists(output_folder) == False: - os.makedirs(output_folder) - - np.savez( - "streamline_coords", - all_coords=all_coords, - ) - print(f"All occurrences of {value_to_find} are at indices: {all_coords}") + return all_coords - 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 + + # Save the minimum cost and the corresponding parameters + min_cost = cost_func(resx, estimator, ansatz, h.total_hamiltonian) # type: ignore + + # Scatter optimal point on cost Landscape ---------------------------- + if args[4]: + if reps == 1: + fig, ax1, ax2 = plt_cost_func(estimator, ansatz, h) + ax1.scatter( # type: ignore + resx[0], resx[1], min_cost, color="red", + marker="o", s=100, label="Optimal Point" # type: ignore + ) + ax2.scatter( # type: ignore + resx[0], resx[1], s=100, color="red", + marker="o", label="Optimal Point" # type: ignore + ) + plt.savefig("Opt_point_visu.png") + print("Optimal point saved in Opt_point_visu.png") + if not args[5]: + plt.show() + else: + pass + # ----------------------------------------------------- + + 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 \ No newline at end of file diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 63a3887..722802b 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -70,20 +70,18 @@ def _build_arg_parser(): return p -def build_mat(in_nodes_mask, in_sh, out_graph, - keep_mask=None, threshold=0.2, slice_index=None, - axis_name="axial", sh_order=8, save_only=False): - - - nodes_mask_im = nib.load(in_nodes_mask) - sh_im = nib.load(in_sh) +#main used to test from command instructions +def main(): + parser = _build_arg_parser() + args = parser.parse_args() + nodes_mask_im = nib.load(args.in_nodes_mask) + sh_im = nib.load(args.in_sh) nodes_mask = nodes_mask_im.get_fdata().astype(bool) - keep_node_indices = None - if keep_mask: - keep_mask = nib.load(keep_mask).get_fdata().astype(bool) + if args.keep_mask: + keep_mask = nib.load(args.keep_mask).get_fdata().astype(bool) keep_node_indices = np.flatnonzero(keep_mask) sh = sh_im.get_fdata() @@ -93,7 +91,7 @@ def build_mat(in_nodes_mask, in_sh, out_graph, # assign edge weights weighted_graph, node_indices = build_weighted_graph( - adj_matrix, node_indices, sh, sh_order + adj_matrix, node_indices, sh, args.sh_order ) # # Could be added in the code if needed: @@ -106,7 +104,7 @@ def build_mat(in_nodes_mask, in_sh, out_graph, # # _______________________________________________________________________ # filter graph edges by weight - weighted_graph[weighted_graph < threshold] = 0.0 + weighted_graph[weighted_graph < args.threshold] = 0.0 # remove intermediate nodes that connect only two nodes weighted_graph = remove_intermediate_connections( @@ -117,41 +115,70 @@ def build_mat(in_nodes_mask, in_sh, out_graph, weighted_graph, node_indices = remove_orphan_nodes( weighted_graph, node_indices, keep_node_indices ) - if slice_index is not None: + + if args.slice_index is not None: weighted_graph, node_indices = extract_slice_at_index( - weighted_graph, node_indices, nodes_mask.shape, slice_index, axis_name + weighted_graph, node_indices, nodes_mask.shape, args.slice_index, args.axis_name ) - if not save_only: + if not args.save_only: plt.imshow(np.log(weighted_graph + 1)) plt.show() # print("node indices", node_indices) # save output - save_graph(weighted_graph, node_indices, nodes_mask.shape, out_graph) - print("Graph saved") - path = Path(out_graph) + save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) + # rap_funct( + # args.out_graph, + # starting_node=node_indices[0], + # ending_node=node_indices[3], + # output_file="rap_output", + # plt_cost_landscape=False, + # save_only=True + # ) - rap_funct( - path, - starting_node=node_indices[0], - ending_node=node_indices[3], - output_file="rap_output", - plt_cost_landscape=False, - save_only=True - ) -def main(): - parser = _build_arg_parser() - args = parser.parse_args() - nodes_mask_im = nib.load(args.in_nodes_mask) - sh_im = nib.load(args.in_sh) + +def quack_rap(in_nodes_mask, in_sh, out_graph, + keep_mask=None, threshold=0.2, slice_index=None, + axis_name="axial", sh_order=8): + """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). + out_graph : str + Output graph file name (npz file). + 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. + + 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 = nodes_mask_im.get_fdata().astype(bool) + keep_node_indices = None - if args.keep_mask: - keep_mask = nib.load(args.keep_mask).get_fdata().astype(bool) + if keep_mask: + keep_mask = nib.load(keep_mask).get_fdata().astype(bool) keep_node_indices = np.flatnonzero(keep_mask) sh = sh_im.get_fdata() @@ -161,7 +188,7 @@ def main(): # assign edge weights weighted_graph, node_indices = build_weighted_graph( - adj_matrix, node_indices, sh, args.sh_order + adj_matrix, node_indices, sh, sh_order ) # # Could be added in the code if needed: @@ -174,7 +201,7 @@ def main(): # # _______________________________________________________________________ # filter graph edges by weight - weighted_graph[weighted_graph < args.threshold] = 0.0 + weighted_graph[weighted_graph < threshold] = 0.0 # remove intermediate nodes that connect only two nodes weighted_graph = remove_intermediate_connections( @@ -185,28 +212,18 @@ def main(): weighted_graph, node_indices = remove_orphan_nodes( weighted_graph, node_indices, keep_node_indices ) - - if args.slice_index is not None: + if slice_index is not None: weighted_graph, node_indices = extract_slice_at_index( - weighted_graph, node_indices, nodes_mask.shape, args.slice_index, args.axis_name + weighted_graph, node_indices, nodes_mask.shape, slice_index, axis_name ) - if not args.save_only: - plt.imshow(np.log(weighted_graph + 1)) - plt.show() - - # print("node indices", node_indices) - # save output - save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) + #function to process the graph before quantum path finding rap_funct( - args.out_graph, + out_graph, starting_node=node_indices[0], ending_node=node_indices[3], - output_file="rap_output", plt_cost_landscape=False, - save_only=True ) - print("Graph saved") if __name__ == "__main__": diff --git a/scripts/quac_optimal_path_find_max_intensity_diffusion.py b/scripts/quac_optimal_path_find_max_intensity_diffusion.py index 06ecf2d..f2d2080 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -5,7 +5,7 @@ from quactography.graph.undirected_graph import Graph from quactography.adj_matrix.io import load_graph from quactography.hamiltonian.hamiltonian_qubit_edge import Hamiltonian_qubit_edge -from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge +from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge, multiprocess_qaoa_solver_edge_rap """ @@ -13,124 +13,142 @@ parameters found if only one reps, and returns the optimization results. """ - -# def _build_arg_parser(): -# p = argparse.ArgumentParser( -# description=__doc__, formatter_class=argparse.RawTextHelpFormatter -# ) -# p.add_argument( -# "in_graph", -# help="Adjacency matrix which graph we want path that maximizes weights in graph, (npz file)", -# type=str, -# ) -# p.add_argument("starting_node", -# help="Starting node of the graph", type=int) -# p.add_argument("ending_node", -# help="Ending node of the graph", type=int) -# p.add_argument("output_file", -# help="Output file name (npz file)", type=str) -# p.add_argument("output_directory", -# help="Directory where the files will be outputed", type=str, -# default="data/output_graphs/" -# ) -# p.add_argument( -# "--alphas", -# nargs="+", -# type=float, -# help="List of alphas", -# default=[1.2] -# ) -# p.add_argument( -# "--reps", -# nargs="+", -# type=int, -# help="List of repetitions to run for the QAOA algorithm", -# default=[1], -# ) -# p.add_argument( -# "-npr", -# "--number_processors", -# help="Number of cpu to use for multiprocessing", -# default=1, -# type=int, -# ) -# p.add_argument( -# "--optimizer", -# help="Optimizer to use for the QAOA algorithm", -# default="Differential", -# type=str, -# ) -# p.add_argument( -# "--plt_cost_landscape", -# help="True or False, Plot 3D and 2D of the cost landscape" -# "(for gamma and beta compact set over all possible angles-0.1 incrementation)", -# action="store_false", -# ) -# p.add_argument( -# "--save_only", -# help="Save only the figure without displaying it", -# action="store_true", -# ) - -# return p - -def rap_funct(graph, starting_node, ending_node, output_file, save_only, output_directory="data/output_graphs", alphas=[1.2], - reps=[1], number_processors=2, optimizer="Differential", plt_cost_landscape=False): +def rap_funct(weighted_graph, starting_node, ending_node, alphas=[1.2], + reps=1, number_processors=2, optimizer="Differential"): """ - Function to return the path from starting node to ending node in the graph. + 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. """ - weighted_graph, _, _ = load_graph(graph) graph = Graph(weighted_graph, starting_node, ending_node) # 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]) # print(hamiltonians[0].total_hamiltonian.simplify()) print("\n Calculating qubits as edges......................") - for i in range(len(reps)): - multiprocess_qaoa_solver_edge( - hamiltonians, - reps[i], - number_processors, - output_file, - output_directory, - optimizer, - plt_cost_landscape, - save_only, - graph.number_of_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, + optimizer, + graph.number_of_edges + ) + return line + +def _build_arg_parser(): + p = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawTextHelpFormatter + ) + p.add_argument( + "in_graph", + help="Adjacency matrix which graph we want path that maximizes weights in graph, (npz file)", + type=str, + ) + p.add_argument("starting_node", + help="Starting node of the graph", type=int) + p.add_argument("ending_node", + help="Ending node of the graph", type=int) + p.add_argument("output_file", + help="Output file name (npz file)", type=str) + p.add_argument("output_directory", + help="Directory where the files will be outputed", type=str, + default="data/output_graphs/" + ) + p.add_argument( + "--alphas", + nargs="+", + type=float, + help="List of alphas", + default=[1.2] + ) + p.add_argument( + "--reps", + nargs="+", + type=int, + help="List of repetitions to run for the QAOA algorithm", + default=[1], + ) + p.add_argument( + "-npr", + "--number_processors", + help="Number of cpu to use for multiprocessing", + default=1, + type=int, + ) + p.add_argument( + "--optimizer", + help="Optimizer to use for the QAOA algorithm", + default="Differential", + type=str, + ) + p.add_argument( + "--plt_cost_landscape", + help="True or False, Plot 3D and 2D of the cost landscape" + "(for gamma and beta compact set over all possible angles-0.1 incrementation)", + action="store_false", + ) + p.add_argument( + "--save_only", + help="Save only the figure without displaying it", + action="store_true", + ) + + return p + +def main(): + parser = _build_arg_parser() + args = parser.parse_args() + + weighted_graph, _, _ = load_graph(args.in_graph) + graph = Graph(weighted_graph, args.starting_node, args.ending_node) -# def main(): -# parser = _build_arg_parser() -# args = parser.parse_args() - -# weighted_graph, _, _ = load_graph(args.in_graph) -# graph = Graph(weighted_graph, args.starting_node, args.ending_node) - -# # Construct Hamiltonian when qubits are set as edges, -# # then optimize with QAOA/scipy: + # Construct Hamiltonian when qubits are set as edges, + # then optimize with QAOA/scipy: -# hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in args.alphas] + hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in args.alphas] -# # print(hamiltonians[0].total_hamiltonian.simplify()) + # print(hamiltonians[0].total_hamiltonian.simplify()) -# print("\n Calculating qubits as edges......................") -# for i in range(len(args.reps)): -# multiprocess_qaoa_solver_edge( -# hamiltonians, -# args.reps[i], -# args.number_processors, -# args.output_file, -# args.output_directory, -# args.optimizer, -# args.plt_cost_landscape, -# args.save_only, -# ) + print("\n Calculating qubits as edges......................") + for i in range(len(args.reps)): + multiprocess_qaoa_solver_edge( + hamiltonians, + args.reps[i], + args.number_processors, + args.output_file, + args.output_directory, + args.optimizer, + args.plt_cost_landscape, + args.save_only, + ) -# if __name__ == "__main__": -# main() +if __name__ == "__main__": + main() From 3ad646c98983da12a9e08f123b934cd80db208d4 Mon Sep 17 00:00:00 2001 From: desk3505 Date: Thu, 24 Jul 2025 12:19:55 -0400 Subject: [PATCH 08/19] correction to setup, and return line to main RAp function --- scripts/quac_matrix_adj_build.py | 4 ++-- setup.cfg | 7 +++---- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 722802b..c51bba0 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -218,13 +218,13 @@ def quack_rap(in_nodes_mask, in_sh, out_graph, ) #function to process the graph before quantum path finding - rap_funct( + line = rap_funct( out_graph, starting_node=node_indices[0], ending_node=node_indices[3], plt_cost_landscape=False, ) - + return line if __name__ == "__main__": main() diff --git a/setup.cfg b/setup.cfg index 729e6d4..5494086 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,11 +3,10 @@ name = quactography version = 0.0.1 [options] -packages = quactography.solver - quactography.adj_matrix - quactography.graph - quactography.hamiltonian +packages = + quactography scripts + data install_requires = dipy==1.10.* From 67c0c26295e1f9289738c2c4ed3a1f355fb47399 Mon Sep 17 00:00:00 2001 From: desk3505 Date: Wed, 30 Jul 2025 10:34:15 -0400 Subject: [PATCH 09/19] Add end node connected to final nodes of ROI (untested) --- quactography/adj_matrix/reconst.py | 31 ++++++++++++++++ quactography/graph/utils.py | 59 ++++++++++++++++++++++++++++++ scripts/quac_matrix_adj_build.py | 30 ++++++++++----- 3 files changed, 111 insertions(+), 9 deletions(-) create mode 100644 quactography/graph/utils.py diff --git a/quactography/adj_matrix/reconst.py b/quactography/adj_matrix/reconst.py index 602f16c..f01e506 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -3,6 +3,7 @@ from dipy.core.sphere import Sphere + DIRECTIONS_2D = np.array( [ [-1.0, 0.0], @@ -30,6 +31,10 @@ 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) @@ -148,6 +153,32 @@ def _add_edge_perhaps( adj_matrix[current_node, neighbor_label] = 1 return adj_matrix +def _add_end_point_edge(adj_matrix, end, node_indes): + """ + 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. + """ + labels = np.unravel_index(node_indes, adj_matrix.shape) + new_shape = (adj_matrix.shape[0] + 1, adj_matrix.shape[1] + 1) + adj_matrix = np.lib.pad(adj_matrix, new_shape, 'constant', constant_values=(0)) + + start = labels.index(end) + 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]: diff --git a/quactography/graph/utils.py b/quactography/graph/utils.py new file mode 100644 index 0000000..ece7a16 --- /dev/null +++ b/quactography/graph/utils.py @@ -0,0 +1,59 @@ +from scipy.ndimage import binary_erosion +import numpy as np + + +def get_output_nodes(mask, entry_node, propagation_direction, angle_rad): + """ + Get output nodes from a hotspot mask. The input mask is expected + to be hole-free. Output nodes are edge nodes found within angle_rad + radians from the propagation direction. + Parameters + ---------- + mask: ndarray + Mask defining the hotspot mask. + entry_node: ndarray (1, 3) + Node we enter the hotspot from. + propagation_direction: ndarray (1, 3) + Direction along which propagation is done. + angle_rad: float + Aperture angle in radians. + Returns + ------- + indices: ndarray (N, 3) + Array of output node indices. + """ + if not mask[entry_node[0], entry_node[1], entry_node[2]]: + raise ValueError('Entry node should be included in mask.') + edges = np.logical_and(mask, ~binary_erosion(mask)) + edge_indices = np.argwhere(edges) + direction_to_edge = edge_indices - np.reshape(entry_node, (1, 3)) + direction_norms = np.linalg.norm(direction_to_edge, axis=-1) + + # remove edge position corresponding to entry_point + direction_to_edge = direction_to_edge[direction_norms > 0] + edge_indices = edge_indices[direction_norms > 0] + direction_norms = direction_norms[direction_norms > 0] + + direction_to_edge = direction_to_edge / direction_norms[:, None] + + propagation_direction_norm = np.linalg.norm(propagation_direction) + if propagation_direction_norm != 1.0: + propagation_direction = propagation_direction / propagation_direction_norm + + directions_cos_angle = direction_to_edge.dot(propagation_direction.reshape((3, 1))) + propagation_cos_angle = np.cos(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) \ No newline at end of file diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index c51bba0..9e23e9b 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -8,13 +8,15 @@ from quactography.adj_matrix.reconst import ( build_adjacency_matrix, - build_weighted_graph + 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.utils import get_output_nodes from quactography.adj_matrix.io import save_graph from scripts.quac_optimal_path_find_max_intensity_diffusion import rap_funct @@ -139,9 +141,9 @@ def main(): -def quack_rap(in_nodes_mask, in_sh, out_graph, +def quack_rap(in_nodes_mask, in_sh, keep_mask=None, threshold=0.2, slice_index=None, - axis_name="axial", sh_order=8): + 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 ---------- @@ -149,8 +151,6 @@ def quack_rap(in_nodes_mask, in_sh, out_graph, Input nodes mask image (.nii.gz file). in_sh : str Input SH image (.nii.gz file). - out_graph : str - Output graph file name (npz file). keep_mask : str, optional Nodes that must not be filtered out. If None, all nodes are filtered. threshold : float, optional @@ -183,9 +183,17 @@ def quack_rap(in_nodes_mask, in_sh, out_graph, sh = sh_im.get_fdata() + # Get end points of the streamline + end_points = get_output_nodes( + nodes_mask, + entry_node=np.array(node_indices[0]), + propagation_direction=prev_direction, + angle_rad=theta + ) + # adjacency graph - adj_matrix, node_indices = build_adjacency_matrix(nodes_mask) - + adj_matrix, node_indices = build_adjacency_matrix(nodes_mask,end_points) + # assign edge weights weighted_graph, node_indices = build_weighted_graph( adj_matrix, node_indices, sh, sh_order @@ -216,14 +224,18 @@ def quack_rap(in_nodes_mask, in_sh, out_graph, weighted_graph, node_indices = extract_slice_at_index( weighted_graph, node_indices, nodes_mask.shape, slice_index, axis_name ) + + # Add end point edges to the adjacency matrix + weighted_graph = _add_end_point_edge(weighted_graph, end_points, node_indices) #function to process the graph before quantum path finding line = rap_funct( - out_graph, + weighted_graph, starting_node=node_indices[0], - ending_node=node_indices[3], + ending_node= weighted_graph.shape[0] -1, plt_cost_landscape=False, ) + line.pop() return line if __name__ == "__main__": From 8ffb69d4e83c86de2796d82a4ede0614a67a0f73 Mon Sep 17 00:00:00 2001 From: desk3505 Date: Fri, 1 Aug 2025 14:29:24 -0400 Subject: [PATCH 10/19] added Size check, and args from scilpy rap call. Finished addading fictional end point for shortest path --- quactography/adj_matrix/reconst.py | 36 +++++++---- scripts/quac_matrix_adj_build.py | 59 ++++++++----------- ...timal_path_find_max_intensity_diffusion.py | 6 +- 3 files changed, 55 insertions(+), 46 deletions(-) diff --git a/quactography/adj_matrix/reconst.py b/quactography/adj_matrix/reconst.py index f01e506..0a1a5f6 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -1,6 +1,7 @@ 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 @@ -66,7 +67,7 @@ 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 + return adj_matrix, node_indices, labels_volume def build_weighted_graph(adj_matrix, node_indices, sh, sh_order=12): @@ -153,7 +154,8 @@ def _add_edge_perhaps( adj_matrix[current_node, neighbor_label] = 1 return adj_matrix -def _add_end_point_edge(adj_matrix, end, node_indes): + +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. @@ -171,13 +173,14 @@ def _add_end_point_edge(adj_matrix, end, node_indes): np.ndarray Updated adjacency matrix with end point edges added. """ - labels = np.unravel_index(node_indes, adj_matrix.shape) - new_shape = (adj_matrix.shape[0] + 1, adj_matrix.shape[1] + 1) - adj_matrix = np.lib.pad(adj_matrix, new_shape, 'constant', constant_values=(0)) - - start = labels.index(end) - adj_matrix[start, -1] = 1 - adj_matrix[-1, start] = 1 # Assuming undirected graph + #labels = np.unravel_index(node_indes, adj_matrix.shape) + #new_shape = (adj_matrix.shape[0] + 1, adj_matrix.shape[1] + 1) + 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): @@ -187,4 +190,17 @@ def _is_valid_pos(pos_x, pos_y, pos_z, nodes_mask): return False if pos_z < 0 or pos_z >= nodes_mask.shape[2]: return False - return nodes_mask[pos_x, pos_y, pos_z] \ No newline at end of file + 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/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 9e23e9b..c6ddeaa 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -9,7 +9,7 @@ from quactography.adj_matrix.reconst import ( build_adjacency_matrix, build_weighted_graph, - _add_end_point_edge + add_end_point_edge ) from quactography.adj_matrix.filter import ( remove_orphan_nodes, @@ -130,18 +130,10 @@ def main(): # print("node indices", node_indices) # save output save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) - # rap_funct( - # args.out_graph, - # starting_node=node_indices[0], - # ending_node=node_indices[3], - # output_file="rap_output", - # plt_cost_landscape=False, - # save_only=True - # ) -def quack_rap(in_nodes_mask, in_sh, +def quack_rap(in_nodes_mask, in_sh, 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). @@ -151,6 +143,8 @@ def quack_rap(in_nodes_mask, in_sh, 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 @@ -162,6 +156,11 @@ def quack_rap(in_nodes_mask, in_sh, 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 ------- @@ -182,32 +181,15 @@ def quack_rap(in_nodes_mask, in_sh, keep_node_indices = np.flatnonzero(keep_mask) sh = sh_im.get_fdata() - - # Get end points of the streamline - end_points = get_output_nodes( - nodes_mask, - entry_node=np.array(node_indices[0]), - propagation_direction=prev_direction, - angle_rad=theta - ) # adjacency graph - adj_matrix, node_indices = build_adjacency_matrix(nodes_mask,end_points) + 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 ) - # # Could be added in the code if needed: - # # Select sub-graph and filter:___________________________________________ - # select_1 = 30 - # select_2 = 200 - - # weighted_graph = weighted_graph[select_1:select_2, select_1:select_2] - # # print(weighted_graph) - # # _______________________________________________________________________ - # filter graph edges by weight weighted_graph[weighted_graph < threshold] = 0.0 @@ -225,18 +207,29 @@ def quack_rap(in_nodes_mask, in_sh, weighted_graph, node_indices, nodes_mask.shape, slice_index, axis_name ) + # Get end points of the streamline + end_points = get_output_nodes( + nodes_mask, + entry_node=np.array(node_indices[0]), + propagation_direction=prev_direction, + angle_rad=theta + ) # Add end point edges to the adjacency matrix - weighted_graph = _add_end_point_edge(weighted_graph, end_points, node_indices) + weighted_graph = add_end_point_edge(weighted_graph, end_points, labels=labes) + if len(np.flatnonzero(weighted_graph))> 17: + print("RAPGraph: max number of points exceeded") + is_line_valid = False + return line, prev_direction, is_line_valid #function to process the graph before quantum path finding line = rap_funct( weighted_graph, - starting_node=node_indices[0], - ending_node= weighted_graph.shape[0] -1, - plt_cost_landscape=False, + starting_node = start_point, + alphas=[alpha], + reps=reps, ) line.pop() - return line + return line, prev_direction, True 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 f2d2080..b187449 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -13,8 +13,8 @@ parameters found if only one reps, and returns the optimization results. """ -def rap_funct(weighted_graph, starting_node, ending_node, alphas=[1.2], - reps=1, number_processors=2, optimizer="Differential"): +def rap_funct(weighted_graph, starting_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. @@ -41,7 +41,7 @@ def rap_funct(weighted_graph, starting_node, ending_node, alphas=[1.2], line : list List of coordinates for the streamline. """ - graph = Graph(weighted_graph, starting_node, ending_node) + graph = Graph(weighted_graph, starting_node) # Construct Hamiltonian when qubits are set as edges, # then optimize with QAOA/scipy: From 29ac1c033f3f887b50188ea7a4851cf893cb4826 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Thu, 7 Aug 2025 10:49:30 -0400 Subject: [PATCH 11/19] functionnal --- scripts/quac_matrix_adj_build.py | 18 ++++++++++-------- ...ptimal_path_find_max_intensity_diffusion.py | 8 ++++---- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index c6ddeaa..033208e 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -208,25 +208,27 @@ def quack_rap(in_nodes_mask, in_sh, start_point, reps, alpha, ) # Get end points of the streamline + end_points = get_output_nodes( nodes_mask, - entry_node=np.array(node_indices[0]), + entry_node=start_point, propagation_direction=prev_direction, angle_rad=theta ) # Add end point edges to the adjacency matrix weighted_graph = add_end_point_edge(weighted_graph, end_points, labels=labes) - if len(np.flatnonzero(weighted_graph))> 17: - print("RAPGraph: max number of points exceeded") - is_line_valid = False - return line, prev_direction, is_line_valid + end = np.flatnonzero(weighted_graph) + if len(end)> 17: + raise Exception("RAPGraph: max number of points exceeded") + #function to process the graph before quantum path finding line = rap_funct( weighted_graph, - starting_node = start_point, - alphas=[alpha], - reps=reps, + starting_node = labes[start_point[0], start_point[1], start_point[2]], + end_node = end[-1], + alphas = [alpha], + reps = reps, ) line.pop() return line, prev_direction, True diff --git a/scripts/quac_optimal_path_find_max_intensity_diffusion.py b/scripts/quac_optimal_path_find_max_intensity_diffusion.py index b187449..37a6caa 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -13,7 +13,7 @@ parameters found if only one reps, and returns the optimization results. """ -def rap_funct(weighted_graph, starting_node, alphas, +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 @@ -41,7 +41,7 @@ def rap_funct(weighted_graph, starting_node, alphas, line : list List of coordinates for the streamline. """ - graph = Graph(weighted_graph, starting_node) + graph = Graph(weighted_graph, starting_node, ending_node) # Construct Hamiltonian when qubits are set as edges, # then optimize with QAOA/scipy: @@ -57,8 +57,8 @@ def rap_funct(weighted_graph, starting_node, alphas, hamiltonians, reps, number_processors, - optimizer, - graph.number_of_edges + graph.number_of_edges, + optimizer ) return line From c1cd161e4cc149010082167126bd81568c9fdb27 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 12 Aug 2025 11:33:40 -0400 Subject: [PATCH 12/19] removed func from script --- quactography/solver/track_from_scilpy.py | 115 +++++++++++++++++++++++ scripts/quac_matrix_adj_build.py | 106 --------------------- 2 files changed, 115 insertions(+), 106 deletions(-) create mode 100644 quactography/solver/track_from_scilpy.py diff --git a/quactography/solver/track_from_scilpy.py b/quactography/solver/track_from_scilpy.py new file mode 100644 index 0000000..559284f --- /dev/null +++ b/quactography/solver/track_from_scilpy.py @@ -0,0 +1,115 @@ +import nibabel as nib +import numpy as np + +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.utils import get_output_nodes +from scripts.quac_optimal_path_find_max_intensity_diffusion import rap_funct + +def quack_rap(in_nodes_mask, in_sh, 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 = nodes_mask_im.get_fdata().astype(bool) + + + keep_node_indices = None + if keep_mask: + keep_mask = nib.load(keep_mask).get_fdata().astype(bool) + keep_node_indices = np.flatnonzero(keep_mask) + + sh = sh_im.get_fdata() + + # 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_points = get_output_nodes( + nodes_mask, + entry_node=start_point, + propagation_direction=prev_direction, + angle_rad=theta + ) + # Add end point edges to the adjacency matrix + weighted_graph = add_end_point_edge(weighted_graph, end_points, labels=labes) + end = np.flatnonzero(weighted_graph) + if len(end)> 17: + raise Exception("RAPGraph: max number of points exceeded") + + + #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]], + end_node = end[-1], + alphas = [alpha], + reps = reps, + ) + line.pop() + return line, prev_direction, True \ No newline at end of file diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 033208e..89833b8 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -4,21 +4,17 @@ import nibabel as nib import numpy as np import matplotlib.pyplot as plt -from pathlib import Path 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.utils import get_output_nodes from quactography.adj_matrix.io import save_graph -from scripts.quac_optimal_path_find_max_intensity_diffusion import rap_funct """ @@ -131,107 +127,5 @@ def main(): # save output save_graph(weighted_graph, node_indices, nodes_mask.shape, args.out_graph) - - -def quack_rap(in_nodes_mask, in_sh, 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 = nodes_mask_im.get_fdata().astype(bool) - - - keep_node_indices = None - if keep_mask: - keep_mask = nib.load(keep_mask).get_fdata().astype(bool) - keep_node_indices = np.flatnonzero(keep_mask) - - sh = sh_im.get_fdata() - - # 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_points = get_output_nodes( - nodes_mask, - entry_node=start_point, - propagation_direction=prev_direction, - angle_rad=theta - ) - # Add end point edges to the adjacency matrix - weighted_graph = add_end_point_edge(weighted_graph, end_points, labels=labes) - end = np.flatnonzero(weighted_graph) - if len(end)> 17: - raise Exception("RAPGraph: max number of points exceeded") - - - #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]], - end_node = end[-1], - alphas = [alpha], - reps = reps, - ) - line.pop() - return line, prev_direction, True - if __name__ == "__main__": main() From cb011628c2280fd8f3ec1c96471b02d495874ed9 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 12 Aug 2025 11:47:47 -0400 Subject: [PATCH 13/19] in order to clena the scripts --- quactography/solver/track_from_scilpy.py | 55 ++++++++++++++++++- ...timal_path_find_max_intensity_diffusion.py | 51 +---------------- 2 files changed, 54 insertions(+), 52 deletions(-) diff --git a/quactography/solver/track_from_scilpy.py b/quactography/solver/track_from_scilpy.py index 559284f..cdd0d07 100644 --- a/quactography/solver/track_from_scilpy.py +++ b/quactography/solver/track_from_scilpy.py @@ -11,8 +11,10 @@ remove_intermediate_connections, extract_slice_at_index ) +from quactography.graph.undirected_graph import Graph from quactography.graph.utils import get_output_nodes -from scripts.quac_optimal_path_find_max_intensity_diffusion import rap_funct +from quactography.hamiltonian.hamiltonian_qubit_edge import Hamiltonian_qubit_edge +from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge_rap def quack_rap(in_nodes_mask, in_sh, start_point, reps, alpha, keep_mask=None, threshold=0.2, slice_index=None, @@ -112,4 +114,53 @@ def quack_rap(in_nodes_mask, in_sh, start_point, reps, alpha, reps = reps, ) line.pop() - return line, prev_direction, True \ No newline at end of file + return line, 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. + """ + graph = Graph(weighted_graph, starting_node, ending_node) + + # 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 + ) + return line \ No newline at end of file diff --git a/scripts/quac_optimal_path_find_max_intensity_diffusion.py b/scripts/quac_optimal_path_find_max_intensity_diffusion.py index 37a6caa..bcf02a7 100644 --- a/scripts/quac_optimal_path_find_max_intensity_diffusion.py +++ b/scripts/quac_optimal_path_find_max_intensity_diffusion.py @@ -5,7 +5,7 @@ from quactography.graph.undirected_graph import Graph from quactography.adj_matrix.io import load_graph from quactography.hamiltonian.hamiltonian_qubit_edge import Hamiltonian_qubit_edge -from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge, multiprocess_qaoa_solver_edge_rap +from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge """ @@ -13,55 +13,6 @@ parameters found if only one reps, and returns the optimization results. """ -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. - """ - graph = Graph(weighted_graph, starting_node, ending_node) - - # 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 - ) - return line - def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawTextHelpFormatter From 7bfcf3caffb06f855989d1aaa1be7ccb0cef461d Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 12 Aug 2025 11:54:21 -0400 Subject: [PATCH 14/19] more significant name --- quactography/solver/{track_from_scilpy.py => rap_tracking.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename quactography/solver/{track_from_scilpy.py => rap_tracking.py} (100%) diff --git a/quactography/solver/track_from_scilpy.py b/quactography/solver/rap_tracking.py similarity index 100% rename from quactography/solver/track_from_scilpy.py rename to quactography/solver/rap_tracking.py From 9f1a9b1b0a6131653ea8a9ed876dbe25ae653ec0 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 12 Aug 2025 13:47:05 -0400 Subject: [PATCH 15/19] classic method, for testing --- quactography/adj_matrix/reconst.py | 2 +- quactography/solver/Dijkstra | 75 +++++++++++++++++++++++++++++ quactography/solver/rap_tracking.py | 57 ++++++++++++---------- 3 files changed, 108 insertions(+), 26 deletions(-) create mode 100644 quactography/solver/Dijkstra diff --git a/quactography/adj_matrix/reconst.py b/quactography/adj_matrix/reconst.py index 0a1a5f6..5683920 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -178,7 +178,7 @@ def add_end_point_edge(adj_matrix, end, labels): 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]] + start = labels[i[0], i[1], i[2]] -1 adj_matrix[start, -1] = 1 adj_matrix[-1, start] = 1 # Assuming undirected graph return adj_matrix diff --git a/quactography/solver/Dijkstra b/quactography/solver/Dijkstra new file mode 100644 index 0000000..36879aa --- /dev/null +++ b/quactography/solver/Dijkstra @@ -0,0 +1,75 @@ +import heapq +import networkx as nx + + +def dijkstra_stepwise(Graph, start, target, diagonal_mode="nondiagonal"): + """ + 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 + + else: + 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/rap_tracking.py b/quactography/solver/rap_tracking.py index cdd0d07..541f49a 100644 --- a/quactography/solver/rap_tracking.py +++ b/quactography/solver/rap_tracking.py @@ -15,8 +15,9 @@ 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(in_nodes_mask, in_sh, start_point, reps, alpha, +def quack_rap(in_nodes_mask_img, in_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). @@ -52,10 +53,10 @@ def quack_rap(in_nodes_mask, in_sh, start_point, reps, alpha, """ - nodes_mask_im = nib.load(in_nodes_mask) - sh_im = nib.load(in_sh) + #nodes_mask_im = nib.load(in_nodes_mask) + #sh_im = nib.load(in_sh) - nodes_mask = nodes_mask_im.get_fdata().astype(bool) + nodes_mask = in_nodes_mask_img.get_fdata().astype(bool) keep_node_indices = None @@ -63,7 +64,7 @@ def quack_rap(in_nodes_mask, in_sh, start_point, reps, alpha, keep_mask = nib.load(keep_mask).get_fdata().astype(bool) keep_node_indices = np.flatnonzero(keep_mask) - sh = sh_im.get_fdata() + sh = in_sh_img.get_fdata() # adjacency graph adj_matrix, node_indices, labes = build_adjacency_matrix(nodes_mask) @@ -101,22 +102,19 @@ def quack_rap(in_nodes_mask, in_sh, start_point, reps, alpha, # Add end point edges to the adjacency matrix weighted_graph = add_end_point_edge(weighted_graph, end_points, labels=labes) end = np.flatnonzero(weighted_graph) - if len(end)> 17: - raise Exception("RAPGraph: max number of points exceeded") - #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]], - end_node = end[-1], + ending_node = end[-1], alphas = [alpha], reps = reps, ) line.pop() return line, prev_direction, True -def rap_funct(weighted_graph, starting_node,ending_node, alphas, +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 @@ -144,23 +142,32 @@ def rap_funct(weighted_graph, starting_node,ending_node, alphas, line : list List of coordinates for the streamline. """ - graph = Graph(weighted_graph, starting_node, ending_node) + # 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: + # # Construct Hamiltonian when qubits are set as edges, + # # then optimize with QAOA/scipy: - hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in alphas] + # hamiltonians = [Hamiltonian_qubit_edge(graph, alpha) for alpha in alphas] - # print(hamiltonians[0].total_hamiltonian.simplify()) + # # 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 - ) + # 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 + # ) + + line = dijkstra_stepwise( + weighted_graph, + starting_node, + ending_node + ) return line \ No newline at end of file From 5d1af10fcaf99524a37668a99131bc8bc95fc615 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 12 Aug 2025 13:50:26 -0400 Subject: [PATCH 16/19] missing .py --- quactography/solver/{Dijkstra => Dijkstra.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename quactography/solver/{Dijkstra => Dijkstra.py} (100%) diff --git a/quactography/solver/Dijkstra b/quactography/solver/Dijkstra.py similarity index 100% rename from quactography/solver/Dijkstra rename to quactography/solver/Dijkstra.py From 2ec5cb82351064f913a2d085a3328ee00ba8d796 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Tue, 12 Aug 2025 14:21:57 -0400 Subject: [PATCH 17/19] changing the Dijkstra ouput into coordinates for the stream line --- quactography/solver/rap_tracking.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/quactography/solver/rap_tracking.py b/quactography/solver/rap_tracking.py index 541f49a..2e12ac1 100644 --- a/quactography/solver/rap_tracking.py +++ b/quactography/solver/rap_tracking.py @@ -112,7 +112,12 @@ def quack_rap(in_nodes_mask_img, in_sh_img, start_point, reps, alpha, reps = reps, ) line.pop() - return line, prev_direction, True + sline = np.unravel_index(line, labes.shape) + llist = [] + for i in range(len(sline[0])): + llist.append([sline[0][i], sline[1][i], sline[2][i]]) + sline = llist + return sline, prev_direction, True def rap_funct(weighted_graph, starting_node, ending_node, alphas, reps, number_processors=2, optimizer="Differential"): From b2923f052081e04ccddc320c933391174f9bbcf3 Mon Sep 17 00:00:00 2001 From: Desk3505 Date: Mon, 18 Aug 2025 10:13:07 -0400 Subject: [PATCH 18/19] preparing to for upstream merge --- data/hurdle_mask.nii.gz | Bin 1665 -> 0 bytes quactography/graph/undirected_graph.py | 6 +- quactography/solver/qaoa_solver_qu_edge.py | 152 --------------------- scripts/quac_matrix_adj_build.py | 2 +- setup.cfg | 3 +- 5 files changed, 5 insertions(+), 158 deletions(-) delete mode 100644 data/hurdle_mask.nii.gz diff --git a/data/hurdle_mask.nii.gz b/data/hurdle_mask.nii.gz deleted file mode 100644 index 6d9dd2b3395bf86b3f4708280a1e835f603d007d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1665 zcmd6n{ZrBh9LLSmop!gaG)uR#WvQ2*=hD;?Q`_xC=cW}YK7_VRy+CyOO%YMPd3(64 z6i-bGQ?Xf-h4@{f7+TOy6GMfW0RcXAGhoS6`4totowmE&^%wN{>GR|3{dsupz3ODQ zQ4VW*`FO1|blXQ~n5iPdu$16jD19!*<3xEpe9^ zh(8Hn#hH?0WbTStubp+t9Jn3GCu^59;t|y}ALI{d2XuQ&tSV81GJHrd%^Oh}XV%<( z5x&uv?6%lDIBoml!)WgBe_f~v3EaMk^BQa4gIG~?kEmwXpgtr%4TGuM#C=d7t@7mC z=(h~xU4C#{wzvh|?tA%Nq({U`v^)$g^#JcN;`8DppyLH%v71yzSlyHB^EP*BV5LF~ z|G^aM`m4!@2?VcpzX0;U=Li5-M(DHfd$@u||}23-fR^VY?b?=q5~ z1MJ*RmWa8=6i-?eIz?q~V zJrG$qMLinkm%)x|Bl18V9p0)@So(&v=EXvM30@Jrx|VgVDGCa`3B}rPU`8ajT7AKa z44u$_M;e|hfEWM`jmpFqa|iCG;L&1^>D zJ&daC$FM*WB0?|tIVTyZ=u{}&d41(CbVgY|cy5xc0`xR)R@umB@ZT$wk>8IBwaxK3 zwt$w3@!5-)xp5Yx_4@Jdd%WQ%n%N6k8Ob}7TAtvXwE{~|Cx8_V&3T#1j4l58xxY+t zz(mE79Cc&|d#=ta%rKXAO4dM+Vs;wJy>nbPwpvL?x9W+VdQzuJkGdZTmlGwlB1bk(c|!B=l~@c zK$ImWj*4n&U^$ESm!i;I%G*s^IN~>%kp(9BysX;B=r~lLtDF-CLHVpUmF@~a#{5iU zc)vaO8C0vaz9Za7eO#*@L=L?U9$aLGk=t@F_^u~@CTqMc8Bnle_}f<{c3-RQJtKxzd{ zV`HK#`8|?Xc$I|Ki=CO^^IS@OsEnbm5c6Fvjw<8%edyw*^GU%^27c^8tg2^%w9(s zcP%o@Z#U7k=~8xtVV1HR@5(v04Ckoxp{|5k%KOJ5sn%m+`9>S^)^)aMs=CQW)~7)^HpB7iwjFPM^F!hb zc~n6g@A-vi&l*$lyHeqQmd6_hM{mvq3#>tb3QI@65Kw!pZCIaT4Mi(sCnd3ajmZVi{FT<^PtHbc4wf#zTq+j39Qd@^I5w zstCdugDl$x#WP?HBSDRDD*coq`*-yB#Sk5wbZ(5* y-%-uIggrz5Gtzseyz96!)s}YJ^Fz|;ras?;-8APr<6o>?v#YFQrkyoyvHKS%fAfO? diff --git a/quactography/graph/undirected_graph.py b/quactography/graph/undirected_graph.py index ef0621f..9ed3a19 100644 --- a/quactography/graph/undirected_graph.py +++ b/quactography/graph/undirected_graph.py @@ -81,9 +81,9 @@ def connexions_edges(self, mat_adj, num_nodes): # print(f"End :{ending_nodes}") weights = [] - for _ in range(10): - for _ in range(10): - for _ in range(10): + for _ in range(number_of_edges): + for _ in starting_nodes: + for _ in ending_nodes: weight_qubit = mat_adj[starting_nodes, ending_nodes] weights.append(weight_qubit) # weights[0][0] diff --git a/quactography/solver/qaoa_solver_qu_edge.py b/quactography/solver/qaoa_solver_qu_edge.py index 63e21b6..472e24f 100644 --- a/quactography/solver/qaoa_solver_qu_edge.py +++ b/quactography/solver/qaoa_solver_qu_edge.py @@ -280,155 +280,3 @@ def multiprocess_qaoa_solver_edge( print( "------------------MULTIPROCESS SOLVER FINISHED-------------------------" ) - -#Only returns solution by classic solution as it is -def multiprocess_qaoa_solver_edge_rap( - hamiltonians, - reps, - nbr_processes, - nb_edges, - optimizer="Differential", -): - """ - 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 - ------- - line : list - List of coordinates for the streamline. - """ - # pool = multiprocessing.Pool(nbr_processes) - # pool.map( - # find_longest_path_rap, - # zip( - # hamiltonians, - # itertools.repeat(reps), - # itertools.repeat(optimizer), - # ), - # ) - exact_path = ((hamiltonians[0].exact_path[0][::-1]).zfill(nb_edges)) - print(f"Exact path: {exact_path}... length: {list(exact_path)}") - all_coords = [] - value_to_find = '1' - - for i, x in enumerate(exact_path): - if x == value_to_find: - all_coords.append(np.unravel_index(i, (hamiltonians[0].graph.number_of_edges, hamiltonians[0].graph.number_of_edges))) - - return all_coords - -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 - - # Save the minimum cost and the corresponding parameters - min_cost = cost_func(resx, estimator, ansatz, h.total_hamiltonian) # type: ignore - - # Scatter optimal point on cost Landscape ---------------------------- - if args[4]: - if reps == 1: - fig, ax1, ax2 = plt_cost_func(estimator, ansatz, h) - ax1.scatter( # type: ignore - resx[0], resx[1], min_cost, color="red", - marker="o", s=100, label="Optimal Point" # type: ignore - ) - ax2.scatter( # type: ignore - resx[0], resx[1], s=100, color="red", - marker="o", label="Optimal Point" # type: ignore - ) - plt.savefig("Opt_point_visu.png") - print("Optimal point saved in Opt_point_visu.png") - if not args[5]: - plt.show() - else: - pass - # ----------------------------------------------------- - - 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 \ No newline at end of file diff --git a/scripts/quac_matrix_adj_build.py b/scripts/quac_matrix_adj_build.py index 89833b8..7f0135e 100644 --- a/scripts/quac_matrix_adj_build.py +++ b/scripts/quac_matrix_adj_build.py @@ -7,7 +7,7 @@ from quactography.adj_matrix.reconst import ( build_adjacency_matrix, - build_weighted_graph, + build_weighted_graph ) from quactography.adj_matrix.filter import ( remove_orphan_nodes, diff --git a/setup.cfg b/setup.cfg index 5494086..c44e367 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,11 +3,10 @@ name = quactography version = 0.0.1 [options] -packages = +packages = quactography scripts data - install_requires = dipy==1.10.* matplotlib==3.6.* From fa09f4c749184e416f68367cbffc295b0b820b42 Mon Sep 17 00:00:00 2001 From: desk3505 Date: Thu, 28 Aug 2025 11:35:36 -0400 Subject: [PATCH 19/19] cleaned, updated libs --- pyproject.toml | 20 ++++ quactography/adj_matrix/reconst.py | 6 +- quactography/solver/Dijkstra.py | 5 +- quactography/solver/qaoa_solver_qu_edge.py | 113 +++++++++++++++++++++ quactography/solver/rap_tracking.py | 63 +++++++----- setup.cfg | 20 ++-- 6 files changed, 187 insertions(+), 40 deletions(-) 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/reconst.py b/quactography/adj_matrix/reconst.py index bdc1ddf..af43ebe 100644 --- a/quactography/adj_matrix/reconst.py +++ b/quactography/adj_matrix/reconst.py @@ -167,12 +167,10 @@ def add_end_point_edge(adj_matrix, end, labels): np.ndarray Updated adjacency matrix with end point edges added. """ - #labels = np.unravel_index(node_indes, adj_matrix.shape) - #new_shape = (adj_matrix.shape[0] + 1, adj_matrix.shape[1] + 1) - adj_matrix = np.lib.pad(adj_matrix, (0, 1), 'constant', constant_values=(0)) + 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]] -1 + start = labels[i[0], i[1], i[2]] adj_matrix[start, -1] = 1 adj_matrix[-1, start] = 1 # Assuming undirected graph return adj_matrix diff --git a/quactography/solver/Dijkstra.py b/quactography/solver/Dijkstra.py index 36879aa..d800fad 100644 --- a/quactography/solver/Dijkstra.py +++ b/quactography/solver/Dijkstra.py @@ -2,7 +2,7 @@ import networkx as nx -def dijkstra_stepwise(Graph, start, target, diagonal_mode="nondiagonal"): +def dijkstra_stepwise(Graph, start, target): """ Perform Dijkstra's algorithm in a stepwise manner on a given graph. Parameters @@ -61,8 +61,7 @@ def dijkstra_stepwise(Graph, start, target, diagonal_mode="nondiagonal"): if current_node == target: break - else: - neighbors = list(G.neighbors(current_node)) + neighbors = list(G.neighbors(current_node)) for neighbor in neighbors: if neighbor not in evaluated_nodes: diff --git a/quactography/solver/qaoa_solver_qu_edge.py b/quactography/solver/qaoa_solver_qu_edge.py index 6cb7eee..164248c 100644 --- a/quactography/solver/qaoa_solver_qu_edge.py +++ b/quactography/solver/qaoa_solver_qu_edge.py @@ -280,3 +280,116 @@ def multiprocess_qaoa_solver_edge( 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 index 2e12ac1..302e18b 100644 --- a/quactography/solver/rap_tracking.py +++ b/quactography/solver/rap_tracking.py @@ -1,5 +1,9 @@ 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, @@ -17,7 +21,8 @@ from quactography.solver.qaoa_solver_qu_edge import multiprocess_qaoa_solver_edge_rap from quactography.solver.Dijkstra import dijkstra_stepwise -def quack_rap(in_nodes_mask_img, in_sh_img, start_point, reps, alpha, + +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). @@ -56,19 +61,17 @@ def quack_rap(in_nodes_mask_img, in_sh_img, start_point, reps, alpha, #nodes_mask_im = nib.load(in_nodes_mask) #sh_im = nib.load(in_sh) - nodes_mask = in_nodes_mask_img.get_fdata().astype(bool) - + 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) - - sh = in_sh_img.get_fdata() - + # 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 @@ -77,47 +80,52 @@ def quack_rap(in_nodes_mask_img, in_sh_img, start_point, reps, alpha, # 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 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 + # ) - # 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 end point edges to the adjacency matrix - weighted_graph = add_end_point_edge(weighted_graph, end_points, labels=labes) - end = np.flatnonzero(weighted_graph) + # 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[-1], + ending_node = end, alphas = [alpha], reps = reps, ) + + # reconstruct the streamline coordinates from the node indices line.pop() - sline = np.unravel_index(line, labes.shape) llist = [] - for i in range(len(sline[0])): - llist.append([sline[0][i], sline[1][i], sline[2][i]]) - sline = llist - return sline, prev_direction, True + 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"): @@ -147,6 +155,8 @@ def rap_funct(weighted_graph, starting_node, ending_node, alphas, 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") @@ -170,8 +180,9 @@ def rap_funct(weighted_graph, starting_node, ending_node, alphas, # optimizer # ) + #---------------CLASSIC DIJKSTRA SOLVER----------------- line = dijkstra_stepwise( - weighted_graph, + -weighted_graph, starting_node, ending_node ) diff --git a/setup.cfg b/setup.cfg index c44e367..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 + 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