Repository: Technion-Kishony-lab/quibbler (language: python)

Issue:
Feature/apply along axis



Discussion / hints:

Fix patch:
diff --git a/examples/compare_images/compare_images.py b/examples/compare_images/compare_images.py
index c90c5f1ed..40c813339 100644
--- a/examples/compare_images/compare_images.py
+++ b/examples/compare_images/compare_images.py
@@ -1,39 +1,60 @@
 # THIS DEMO DOES NOT WORK 100%
 import weakref
-
 import numpy as np
+from functools import partial
+
+# import objgraph
 from matplotlib import pyplot as plt, widgets
+from matplotlib.widgets import AxesWidget
 from mpl_toolkits.axes_grid1 import ImageGrid
 
-from pyquibbler import iquib, override_all, q, quibbler_user_function
+from pyquibbler import iquib, override_all, q
 
 override_all()
 
+axes_widgets = []
+prev = AxesWidget.__init__
+
+
+def wrap_and_save(obj, *args, **kwargs):
+    axes_widgets.append(weakref.ref(obj))
+    return prev(obj, *args, **kwargs)
+
+
+AxesWidget.__init__ = wrap_and_save
+
 
-@quibbler_user_function(lazy=False, receive_quibs=True)
+@partial(np.vectorize, signature='(extents),()->()', pass_quibs=True)
 def create_roi(roi, axes):
     print("Creating ROI")
-    widgets.RectangleSelector(axes, extents=roi)
+    widgets.RectangleSelector(axes, extents=roi, allow_resize=False)
 
 
-@quibbler_user_function(lazy=False)
+@partial(np.vectorize, signature='(w,h,c),(extents)->(w2,h2,c)')
 def cut_image(image, roi):
     print("Cutting image")
-    im = image[roi[2]:roi[3], roi[0]:roi[1]]
-    return im
+    return image[roi[2]:roi[3], roi[0]:roi[1]]
 
 
-@quibbler_user_function(lazy=False)
+@partial(np.vectorize, signature='(w,h,c),(w,h,c)->()')
 def image_distance(img1, img2):
     print("Comparing images")
-    res = np.linalg.norm(np.average(img1, axis=(0, 1)) - np.average(img2, axis=(0, 1))) / 255
-    return res
+    return np.linalg.norm(np.average(img1, axis=(0, 1)) - np.average(img2, axis=(0, 1))) / 255
+
+
+@np.vectorize
+def show_adjacency(axes, x, y, adjacent):
+    axes.scatter(x, y,
+                 s=adjacent * 100 + 1,
+                 marker='x',
+                 color='r',
+                 linewidths=2)
 
 
 file_name = iquib('./pipes.jpg')
 image = plt.imread(file_name)
 
-images_count = iquib(6)
+images_count = iquib(3)
 images_count.set_assignment_template(0, 10, 1)
 
 roi_default = iquib([[10, 110, 10, 110]])
@@ -43,9 +64,8 @@ def image_distance(img1, img2):
 rois.set_assignment_template(0, 1000, 1)
 rois.allow_overriding = True
 
-rois_itered = list(rois.iter_first(images_count.get_value()))
 similiarity_threshold = iquib(.1)
-cut_images_lst = [cut_image(image, roi) for roi in rois_itered]
+cut_images = cut_image(image, rois)
 
 
 def create_figure_1():
@@ -53,8 +73,7 @@ def create_figure_1():
     plt.axes([.1, .2, .8, .7])
     plt.imshow(image, aspect="auto")
 
-    for roi in rois_itered:
-        create_roi(roi, plt.gca())
+    create_roi(rois, plt.gca())
 
     axfreq = plt.axes([0.25, 0.1, 0.65, 0.03])
     widgets.Slider(
@@ -71,53 +90,44 @@ def create_figure_1():
         ax=axfreq,
         label=q("Image count ".format, images_count),
         valmin=1,
-        valmax=6,
-        valstep=.1,
+        valmax=9,
+        valstep=1,
         valinit=images_count
     )
 
 
 def create_figure_2():
-
     # Plot images
     fig = plt.figure(2)
-    grid = ImageGrid(fig, 111,  # similar to subplot(111)
-                     nrows_ncols=(3, 3),  # creates 2x2 grid of axes
-                     axes_pad=0.1,  # pad between axes in inch.
-                     )
+    grid_axes = iquib(np.array(list(ImageGrid(fig, 111,
+                                              nrows_ncols=(3, 3),
+                                              axes_pad=0.1,
+                                              ))))
 
-    for ax, im in zip(grid, cut_images_lst):
-        # Iterating over the grid returns the Axes.
-        ax.imshow(im)
+    np.vectorize(lambda ax, im: ax.imshow(im), signature='(),(w,h,c)->()')(grid_axes[:q(int, images_count)], cut_images)
 
 
 def create_figure_3():
     # Compare sub images
-    image_distances = np.array([[image_distance(img1, img2) for img2 in cut_images_lst] for img1 in cut_images_lst])
+    image_distances = image_distance(np.expand_dims(cut_images, 1), cut_images)
     adjacents = image_distances < similiarity_threshold
 
     # Plot distance matrix
     plt.figure(3)
-    plt.axis([0.5, images_count+0.5, 0.5, images_count+0.5 ])
+    plt.axis([-0.5, images_count - 0.5, -0.5, images_count - 0.5])
     plt.title('pairwise distance between images')
     plt.xlabel('Image number')
     plt.ylabel('Image number')
 
-
-    #
-    for i in range(1, images_count.get_value() + 1):
-        adjacents_for_image = adjacents[i - 1]
-        ss = [(adjacents_for_image[i] * 100 + 1) for i in range(images_count.get_value())]
-        plt.scatter(list(range(1, images_count.get_value() + 1)), np.repeat([i], images_count),
-                    s=ss,
-                    marker='x',
-                    color='r',
-                    linewidths=2)
+    show_adjacency(plt.gca(), np.expand_dims(np.arange(images_count), 1), np.arange(images_count), adjacents)
 
 
 create_figure_1()
 create_figure_2()
 create_figure_3()
 
-
-plt.show(block=True)
+# def save():
+#     objgraph.show_backrefs(axes_widgets[0](), filename='aa.png', max_depth=10)
+# TODO: fix slider
+# TODO: fix scatter
+plt.show()
diff --git a/examples/default_overriding_two_levels.py b/examples/default_overriding_two_levels.py
index 065ef4d23..9f8f19cc5 100644
--- a/examples/default_overriding_two_levels.py
+++ b/examples/default_overriding_two_levels.py
@@ -20,7 +20,7 @@
 plt.bar(x, plate_factor, color=(0.7, 0.7, 0.7))
 plt.plot([-0.5, n_plates - 0.5], default_factor[[0, 0]], linewidth=5, picker=True)
 
-for i in range(len(x)):
+for i in range(x.get_shape()[0]):
     color = q(tuple, [(0, 1, 0), (1, 0, 0)])[plate_factor.get_override_mask()[i]]
     plt.plot(x[i], plate_factor[i],
              marker='s', markerfacecolor=color, markersize=9,
@@ -40,7 +40,7 @@
 plt.bar(xx, well_factor, color=(0.7, 0.7, 0.7), width=0.1)
 plt.axis([-0.5, n_plates - 0.5, 0, 100])
 
-for i in range(len(xx)):
+for i in range(xx.get_shape()[0]):
     color = q(tuple, [(0, 1, 0), (1, 0, 0)])[well_factor.get_override_mask()[i]]
     plt.plot(xx[i], well_factor[i],
              marker='s', markerfacecolor=color, markersize=9,
diff --git a/pyquibbler/decorators/quibbler_user_function.py b/pyquibbler/decorators/quibbler_user_function.py
index 57c1f5985..0416d7023 100644
--- a/pyquibbler/decorators/quibbler_user_function.py
+++ b/pyquibbler/decorators/quibbler_user_function.py
@@ -1,9 +1,7 @@
 import functools
 
-from pyquibbler.quib.graphics import GraphicsFunctionQuib
 
-
-def quibbler_user_function(lazy=True, receive_quibs=False):
+def quibbler_user_function(lazy=True, pass_quibs=False):
     """
     Decorate your function with this in order for quibbler to automatically unpack quibs sent as arguments to this
     function, while reruninng this function every time any argument quib changes.
@@ -12,13 +10,16 @@ def quibbler_user_function(lazy=True, receive_quibs=False):
 
     :param lazy: Should this function be run immediately (lazy=False), or only when it's needed down the line?
     (lazy=False) - if you do any graphics in this function you should probably set lazy=False
-    :param receive_quibs: Should this function receive quibs or their values?
+    :param pass_quibs: Should this function receive quibs or their values?
     """
+    from pyquibbler.quib.graphics import GraphicsFunctionQuib
+
     def _decorator(func):
         @functools.wraps(func)
         def _wrapper(*args, **kwargs):
             return GraphicsFunctionQuib.create(func=func, func_args=args, func_kwargs=kwargs, lazy=lazy,
-                                               receive_quibs=receive_quibs)
+                                               pass_quibs=pass_quibs)
 
         return _wrapper
+
     return _decorator
diff --git a/pyquibbler/env.py b/pyquibbler/env.py
index 208874330..53a242906 100644
--- a/pyquibbler/env.py
+++ b/pyquibbler/env.py
@@ -4,6 +4,8 @@
 
 LAZY = Flag(True)
 
+GRAPHICS_LAZY = Flag(False)
+
 # Assignment restrictions are incomplete at the moment -
 # they consider changes to the same quib contradictory, while
 # only changes to the same paths should be.
diff --git a/pyquibbler/quib/__init__.py b/pyquibbler/quib/__init__.py
index 23b711198..52aa5193c 100644
--- a/pyquibbler/quib/__init__.py
+++ b/pyquibbler/quib/__init__.py
@@ -4,8 +4,9 @@
 from .input_quib import iquib
 from .operator_overriding import override_quib_operators
 from .quib import Quib, OverridingNotAllowedException, QuibIsNotNdArrayException
-from .assignment import Assignment
+from .assignment import Assignment, PathComponent
 from .override_choice import AssignmentNotPossibleException, get_overrides_for_assignment, \
     get_overrides_for_assignment_group
+from .proxy_quib import ProxyQuib
 
 override_quib_operators()
diff --git a/pyquibbler/quib/function_quibs/cache/shallow/nd_cache/nd_unstructured_array_cache.py b/pyquibbler/quib/function_quibs/cache/shallow/nd_cache/nd_unstructured_array_cache.py
index aaedc0726..a11a8c5df 100644
--- a/pyquibbler/quib/function_quibs/cache/shallow/nd_cache/nd_unstructured_array_cache.py
+++ b/pyquibbler/quib/function_quibs/cache/shallow/nd_cache/nd_unstructured_array_cache.py
@@ -29,8 +29,7 @@ def _get_all_uncached_paths(self) -> List[List[PathComponent]]:
     def _is_completely_invalid(self):
         return np.all(self._invalid_mask)
 
-    def _get_uncached_paths_at_path_component(self,
-                                              path_component: PathComponent) -> List[List[PathComponent]]:
+    def _get_uncached_paths_at_path_component(self, path_component: PathComponent) -> List[List[PathComponent]]:
         boolean_mask_of_indices = create_empty_array_with_values_at_indices(
             self._value.shape,
             indices=path_component.component,
diff --git a/pyquibbler/quib/function_quibs/default_function_quib.py b/pyquibbler/quib/function_quibs/default_function_quib.py
index 183f3a8a7..352a7d59b 100644
--- a/pyquibbler/quib/function_quibs/default_function_quib.py
+++ b/pyquibbler/quib/function_quibs/default_function_quib.py
@@ -1,14 +1,14 @@
 from __future__ import annotations
+import numpy as np
 from sys import getsizeof
 from time import perf_counter
 from typing import Callable, Any, Mapping, Tuple, Optional, List, TYPE_CHECKING
 
-import numpy as np
-
 from pyquibbler.quib.function_quibs.cache import create_cache
 from .cache.cache import CacheStatus
 from .cache.holistic_cache import PathCannotHaveComponentsException
 from .cache.shallow.indexable_cache import transform_cache_to_nd_if_necessary_given_path
+
 from .function_quib import FunctionQuib, CacheBehavior
 from ..assignment import AssignmentTemplate
 from ..assignment.utils import get_sub_data_from_object_in_path
@@ -53,8 +53,9 @@ def _ensure_cache_matches_result(self, new_result: Any):
             self._cache = create_cache(new_result)
         return self._cache
 
-    def _invalidate_self(self, path: List['PathComponent']):
+    def _invalidate_self(self, path: List[PathComponent]):
         if len(path) == 0:
+            self._on_type_change()
             self._reset_cache()
 
         if self._cache is not None:
@@ -76,14 +77,13 @@ def _should_cache(self, result: Any, elapsed_seconds: float):
         return elapsed_seconds > self.MIN_SECONDS_FOR_CACHE \
             and getsizeof(result) / elapsed_seconds < self.MAX_BYTES_PER_SECOND
 
-    def _get_uncached_paths_matching_path(self, path: Optional[List['PathComponent']]):
+    def _get_uncached_paths_matching_path(self, path: List[Optional[List['PathComponent']]]):
         """
         Get a list of paths that are uncached within the given path- these paths must be a subset of the given path
         (or the path itself)
         """
 
         if self._cache is not None:
-
             if path is None:
                 # We need to be valid at no paths, so by definitions we also have no uncached paths that match no paths
                 return []
@@ -124,7 +124,7 @@ def _truncate_path_to_match_shallow_caches(self, path: Optional[List['PathCompon
         return new_path
 
     def _run_func_on_uncached_paths(self, truncated_path: List[PathComponent],
-                                    uncached_paths: List[List[PathComponent]]):
+                                    uncached_paths: List[Optional[List[PathComponent]]]):
         """
         Run the function a list of uncached paths, given an original truncated path, storing it in our cache
         """
diff --git a/pyquibbler/quib/function_quibs/elementwise_function_quib.py b/pyquibbler/quib/function_quibs/elementwise_function_quib.py
index d415b795f..191f7ee14 100644
--- a/pyquibbler/quib/function_quibs/elementwise_function_quib.py
+++ b/pyquibbler/quib/function_quibs/elementwise_function_quib.py
@@ -172,7 +172,7 @@ def _forward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> A
          [True, False, False]]
         """
         bool_mask = self._get_source_shaped_bool_mask(quib, indices)
-        return np.broadcast_to(bool_mask, self.get_shape().get_value())
+        return np.broadcast_to(bool_mask, self.get_shape())
 
     def raise_if_multiple_args_have_common_ancestor(self):
         """
@@ -200,9 +200,9 @@ def _get_indices_to_change(self, argument_quib: Quib, working_indices) -> Any:
         and broadcast it's index grid to the shape of the result, so we can see the corresponding quib indices for the
         result indices
         """
-        index_grid = np.indices(argument_quib.get_shape().get_value())
+        index_grid = np.indices(argument_quib.get_shape())
         broadcasted_grid = np.broadcast_to(index_grid,
-                                           (index_grid.shape[0], *self.get_shape().get_value()))
+                                           (index_grid.shape[0], *self.get_shape()))
         return tuple([
             dimension[working_indices]
             for dimension in broadcasted_grid
diff --git a/pyquibbler/quib/function_quibs/function_quib.py b/pyquibbler/quib/function_quibs/function_quib.py
index f6929dc37..f4c2bc3a5 100644
--- a/pyquibbler/quib/function_quibs/function_quib.py
+++ b/pyquibbler/quib/function_quibs/function_quib.py
@@ -1,7 +1,5 @@
 from __future__ import annotations
-
 import functools
-
 import numpy as np
 import types
 from enum import Enum
@@ -59,7 +57,7 @@ def parents(self) -> Set[Quib]:
         return set(iter_quibs_in_args(self.args, self.kwargs))
 
     @classmethod
-    def create(cls, func, func_args=(), func_kwargs=None, cache_behavior=None, **kwargs):
+    def create(cls, func, func_args=(), func_kwargs=None, cache_behavior=None, **init_kwargs):
         """
         Public constructor for FunctionQuib.
         """
@@ -78,7 +76,7 @@ def create(cls, func, func_args=(), func_kwargs=None, cache_behavior=None, **kwa
                        for k, v in func_kwargs.items()}
         func_args = deep_copy_without_quibs_or_artists(func_args)
         self = cls(func=func, args=func_args, kwargs=func_kwargs,
-                   cache_behavior=cache_behavior, **kwargs)
+                   cache_behavior=cache_behavior, **init_kwargs)
         for arg in iter_quibs_in_args(func_args, func_kwargs):
             arg.add_child(self)
         if not LAZY:
@@ -99,12 +97,12 @@ def _add_ufunc_members_to_wrapper(cls, func: Callable, quib_supporting_func_wrap
                 setattr(quib_supporting_func_wrapper, member, val)
 
     @classmethod
-    def _wrapper_call(cls, func, args, kwargs):
+    def _wrapper_call(cls, func, args, kwargs, **create_kwargs):
         """
         The actual code that runs when a quib-supporting function wrapper is called.
         """
         if is_there_a_quib_in_args(args, kwargs):
-            return cls.create(func=func, func_args=args, func_kwargs=kwargs)
+            return cls.create(func, args, kwargs, **create_kwargs)
 
         return func(*args, **kwargs)
 
@@ -206,8 +204,7 @@ def _prepare_args_for_call(self, valid_path: Optional[List[PathComponent]]):
         quibs_to_paths = {} if valid_path is None else self._get_source_paths_of_quibs_given_path(valid_path)
         replace_func = functools.partial(self._replace_sub_argument_with_value, quibs_to_paths)
         new_args = [recursively_run_func_on_object(replace_func, arg) for arg in self.args]
-        new_kwargs = {key: recursively_run_func_on_object(replace_func, val)
-                      for key, val in self.kwargs.items()}
+        new_kwargs = {key: recursively_run_func_on_object(replace_func, val) for key, val in self.kwargs.items()}
         return new_args, new_kwargs
 
     def _call_func(self, valid_path: Optional[List[PathComponent]]):
@@ -221,7 +218,7 @@ def _call_func(self, valid_path: Optional[List[PathComponent]]):
         valid at valid_path
         """
         new_args, new_kwargs = self._prepare_args_for_call(valid_path)
-        return self._func(*new_args, **new_kwargs)
+        return self.func(*new_args, **new_kwargs)
 
     def _forward_translate_invalidation_path(self, invalidator_quib: Quib, path: List[PathComponent]) -> \
             List[List[PathComponent]]:
diff --git a/pyquibbler/quib/function_quibs/indices_translator_function_quib.py b/pyquibbler/quib/function_quibs/indices_translator_function_quib.py
index 2eefc5f92..20845986c 100644
--- a/pyquibbler/quib/function_quibs/indices_translator_function_quib.py
+++ b/pyquibbler/quib/function_quibs/indices_translator_function_quib.py
@@ -2,7 +2,7 @@
 from abc import abstractmethod
 from dataclasses import dataclass
 from functools import lru_cache
-from typing import Set, Optional, Dict, Callable, List, Any, Union
+from typing import Set, Optional, Dict, Callable, List, Any, Union, Tuple
 
 from pyquibbler.quib.assignment import Assignment
 from pyquibbler.quib.function_quibs.utils import create_empty_array_with_values_at_indices, ArgsValues
@@ -11,14 +11,14 @@
 from pyquibbler.quib.assignment import PathComponent, QuibWithAssignment
 from pyquibbler.quib.utils import iter_objects_of_type_in_object_shallowly
 
-Args = List[Any]
+Args = Tuple[Any, ...]
 Kwargs = Dict[str, Any]
 
 
 @dataclass
 class SupportedFunction:
     data_source_indices: Union[Set[Union[int, str, slice]],
-                               Callable[[Args, Kwargs, Kwargs], List[Any]]]
+                               Callable[[ArgsValues], List[Any]]]
 
     def get_data_source_args(self, args_values: ArgsValues) -> List[Any]:
         if callable(self.data_source_indices):
@@ -38,7 +38,7 @@ def _get_source_shaped_bool_mask(self, invalidator_quib: Quib, indices: Any) ->
             value=True,
             empty_value=False,
             indices=indices,
-            shape=invalidator_quib.get_shape().get_value()
+            shape=invalidator_quib.get_shape()
         )
 
     def _get_bool_mask_representing_indices_in_result(self, indices) -> Union[np.ndarray, bool]:
@@ -46,7 +46,7 @@ def _get_bool_mask_representing_indices_in_result(self, indices) -> Union[np.nda
         Get a boolean mask representing where the indices that were changed are in the result- this will be in
         same shape as the result
         """
-        return create_empty_array_with_values_at_indices(self.get_shape().get_value(),
+        return create_empty_array_with_values_at_indices(self.get_shape(),
                                                          indices=indices, value=True,
                                                          empty_value=False)
 
@@ -56,7 +56,7 @@ def _get_representative_result(self, working_component, value):
         not a `PathComponent`) set to a value
         """
         return create_empty_array_with_values_at_indices(
-            self.get_shape().get_value(),
+            self.get_shape(),
             indices=working_component,
             value=value,
         )
@@ -80,9 +80,6 @@ def _get_data_source_quibs(self) -> Set:
     def _forward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
         pass
 
-    def _get_source_paths_of_quibs_given_path(self, filtered_path_in_result: List[PathComponent]):
-        return {}
-
     def _get_quibs_to_relevant_result_values(self, assignment: Assignment):
         return {}
 
diff --git a/pyquibbler/quib/function_quibs/transpositional/getitem_function_quib.py b/pyquibbler/quib/function_quibs/transpositional/getitem_function_quib.py
index feeaf549b..5cff8e5ae 100644
--- a/pyquibbler/quib/function_quibs/transpositional/getitem_function_quib.py
+++ b/pyquibbler/quib/function_quibs/transpositional/getitem_function_quib.py
@@ -59,6 +59,7 @@ def _forward_translate_invalidation_path(self, quib: Quib,
         # we simply want to check if our getitem's item is equal to the invalidations item (ie it's path).
         # If so, invalidate. This is true for field arrays as well (We do need to
         # add support for indexing multiple fields).
+        assert not isinstance(working_component.component, np.ndarray)
         if self.args[1] == working_component.component:
             return [rest_of_path]
 
diff --git a/pyquibbler/quib/function_quibs/transpositional/transpositional_function_quib.py b/pyquibbler/quib/function_quibs/transpositional/transpositional_function_quib.py
index 07290c2a8..4b8603baf 100644
--- a/pyquibbler/quib/function_quibs/transpositional/transpositional_function_quib.py
+++ b/pyquibbler/quib/function_quibs/transpositional/transpositional_function_quib.py
@@ -25,7 +25,10 @@ class TranspositionalFunctionQuib(DefaultFunctionQuib, IndicesTranslatorFunction
         getitem: SupportedFunction({0}),
         np.reshape: SupportedFunction({0}),
         np.ravel: SupportedFunction({0}),
-        np.squeeze: SupportedFunction({0})
+        np.squeeze: SupportedFunction({0}),
+        np.array: SupportedFunction({0}),
+        np.swapaxes: SupportedFunction({0}),
+        np.expand_dims: SupportedFunction({0})
     }
 
     def _replace_quibs_in_arguments_which_can_potentially_change(self, replace_func: Callable[['Quib'], Any]):
@@ -59,7 +62,7 @@ def _get_quib_ids_mask(self) -> np.ndarray:
 
         def replace_quib_with_id(obj):
             if isinstance(obj, Quib) and obj in self._get_data_source_quibs():
-                return np.full(obj.get_shape().get_value(), quibs_to_ids[obj])
+                return np.full(obj.get_shape(), quibs_to_ids[obj])
             return obj
 
         new_arguments = recursively_run_func_on_object(
@@ -74,7 +77,7 @@ def _get_quibs_to_index_grids(self) -> Dict[Quib, np.ndarray]:
         Get a mapping between quibs and their indices grid
         """
         return {
-            quib: np.indices(quib.get_shape().get_value())
+            quib: np.indices(quib.get_shape())
             for quib in self._get_data_source_quibs()
         }
 
@@ -125,12 +128,12 @@ def get_quibs_to_indices_in_quibs(self, filtered_indices_in_result: Any) -> Dict
         """
         relevant_indices_mask = create_empty_array_with_values_at_indices(
             indices=filtered_indices_in_result,
-            shape=self.get_shape().get_value(),
+            shape=self.get_shape(),
             value=True,
             empty_value=False
         )
         quibs = self._get_data_source_quibs()
-        max_shape_length = max([len(quib.get_shape().get_value())
+        max_shape_length = max([len(quib.get_shape())
                                 for quib in quibs]) if len(quibs) > 0 else 0
         # Default is to set all - if we have a shape we'll change this
         quibs_to_indices_in_quibs = {
@@ -194,7 +197,7 @@ def _replace_arg_with_corresponding_mask_or_arg(q):
                 if q is quib:
                     return self._get_source_shaped_bool_mask(q, indices)
                 else:
-                    return np.full(q.get_shape().get_value(), False)
+                    return np.full(q.get_shape(), False)
             return q
 
         new_arguments = recursively_run_func_on_object(
diff --git a/pyquibbler/quib/graphics/__init__.py b/pyquibbler/quib/graphics/__init__.py
index 360cbe8cd..1cdab594a 100644
--- a/pyquibbler/quib/graphics/__init__.py
+++ b/pyquibbler/quib/graphics/__init__.py
@@ -1,7 +1,8 @@
 from .event_handling import CanvasEventHandler
 from .graphics_function_quib import GraphicsFunctionQuib
-from .redraw import redraw_axes
-from .global_collecting import ArtistsCollector, is_within_artists_collector
+from .redraw import redraw_axeses
+from .global_collecting import ArtistsCollector, AxesWidgetsCollector
 from .plot_graphics_function_quib import PlotGraphicsFunctionQuib
-from .axiswise_function_quib import ReductionAxisWiseGraphicsFunctionQuib, AlongAxisGraphicsFunctionQuib, \
-    VectorizeGraphicsFunctionQuib
+from .axiswise_function_quibs import ReductionAxisWiseGraphicsFunctionQuib
+from .axiswise_function_quibs import ApplyAlongAxisGraphicsFunctionQuib
+from .axiswise_function_quibs import VectorizeGraphicsFunctionQuib, QVectorize
diff --git a/pyquibbler/quib/graphics/axiswise_function_quib.py b/pyquibbler/quib/graphics/axiswise_function_quib.py
deleted file mode 100644
index fa0588817..000000000
--- a/pyquibbler/quib/graphics/axiswise_function_quib.py
+++ /dev/null
@@ -1,304 +0,0 @@
-import numpy as np
-from dataclasses import dataclass
-from itertools import chain
-from functools import cached_property
-from abc import abstractmethod
-from typing import Any, Dict, Optional, List
-
-from pyquibbler.quib.quib import Quib
-from pyquibbler.quib.utils import copy_and_replace_quibs_with_vals
-
-from .graphics_function_quib import GraphicsFunctionQuib
-from ..assignment import PathComponent
-from ..function_quibs.indices_translator_function_quib import IndicesTranslatorFunctionQuib, SupportedFunction
-from ..function_quibs.utils import ArgsValues
-from ..utils import iter_args_and_names_in_function_call
-
-
-def get_vectorize_call_data_args(args_values: ArgsValues) -> List[Any]:
-    """
-    Given a call to a vectorized function, return the arguments which act as data sources.
-    We are using args_values.args and args_values.kwargs instead of the full args dict on purpose,
-    to match vectorize function behavior.
-    """
-    vectorize, *args = args_values.args
-    return [val for key, val in chain(enumerate(args), args_values.kwargs.items()) if key not in vectorize.excluded]
-
-
-class VectorizeGraphicsFunctionQuib(GraphicsFunctionQuib, IndicesTranslatorFunctionQuib):
-    """
-    np.vectorize, without the signature parameter, turns normal python functions into elementwise functions.
-    The signature parameter allows specifying the amount of core dimensions each argument has (and the result)
-    and thereby turns the function into an axiswise function. When a vectorized function is called,
-    the loop dimensions of all non-excluded arguments are broadcast together.
-    The dimension of the result will be the broadcast loop dimension, plus the result core dimension specified
-    in the signature, () by default.
-    """
-    SUPPORTED_FUNCTIONS = {
-        np.vectorize.__call__: SupportedFunction(get_vectorize_call_data_args)
-    }
-
-    @cached_property
-    def _arg_values_in_wrapped_function(self):
-        return list(dict(iter_args_and_names_in_function_call(self._vectorize.pyfunc,
-                                                              self.args[1:], self.kwargs, False)).values())
-
-    @cached_property
-    def _vectorize(self) -> np.vectorize:
-        """
-        Get the vectorize object we were called with
-        """
-        return self.args[0]
-
-    def _get_arg_core_ndim(self, arg: Any) -> int:
-        """
-        Return the core ndim of the argument in the given index.
-        """
-        arg_index = next(i for i, val in enumerate(self._arg_values_in_wrapped_function) if val is arg)
-        if self._vectorize._in_and_out_core_dims is not None:
-            return len(self._vectorize._in_and_out_core_dims[0][arg_index])
-        return 0
-
-    def _get_result_core_ndim(self):
-        """
-        Assuming the result is not a tuple, return the core ndim in the result.
-        """
-        if self._vectorize._in_and_out_core_dims is not None:
-            out_core_dims = self._vectorize._in_and_out_core_dims[1]
-            assert len(out_core_dims) == 1
-            return len(out_core_dims[0])
-        return 0
-
-    @property
-    def _loop_shape(self):
-        """
-        Get the shape of the vectorize loop.
-        This could also be done be broadcasting together the loop dimensions of all arguments.
-        """
-        if self._get_tuple_output_len():
-            result_shape = self.get_value()[0].shape
-        else:
-            result_shape = self.get_shape().get_value()
-        in_and_out_core_dims = self._vectorize._in_and_out_core_dims
-        if in_and_out_core_dims is None:
-            out_core_dims = 0
-        else:
-            out_core_dims = len(in_and_out_core_dims[1][0])
-        return result_shape[:-out_core_dims] if out_core_dims else result_shape
-
-    def _get_sample_elementwise_func_result(self):
-        """
-        Assuming no argument has core dimensions (and therefore the vectorized function must be elementwise),
-        run the vectorized function on one set of arguments and return the result.
-        """
-        args = [np.asarray(copy_and_replace_quibs_with_vals(arg)).flat[0] for arg in self.args[1:]]
-        kwargs = {key: np.asarray(copy_and_replace_quibs_with_vals(val)).flat[0] for key, val in
-                  self.kwargs.items()}
-        return self._vectorize.pyfunc(*args, **kwargs)
-
-    def _get_tuple_output_len(self) -> Optional[int]:
-        """
-        If this vectorize function returns a tuple, return its length. Otherwise, return None.
-        """
-        in_and_out_core_dims = self._vectorize._in_and_out_core_dims
-        if in_and_out_core_dims is not None:
-            out_core_dims_len = len(in_and_out_core_dims[1])
-            return None if out_core_dims_len == 1 else out_core_dims_len
-
-        # If no signature is given, args core dimensions are all () so we can just flatten them.
-        result = self._get_sample_elementwise_func_result()
-        if isinstance(result, tuple):
-            return len(result)
-        return None
-
-    def _forward_translate_indices_to_bool_mask(self, invalidator_quib: Quib, indices: Any) -> Any:
-        source_bool_mask = self._get_source_shaped_bool_mask(invalidator_quib, indices)
-        core_dims = self._get_arg_core_ndim(invalidator_quib)
-        if core_dims > 0:
-            source_bool_mask = np.any(source_bool_mask, axis=tuple(range(source_bool_mask.ndim)[-core_dims:]))
-        return np.broadcast_to(source_bool_mask, self._loop_shape)
-
-    def _forward_translate_invalidation_path(self, invalidator_quib: Quib,
-                                             path: List[PathComponent]) -> List[Optional[List[PathComponent]]]:
-        working_component, *rest_of_path = path
-        bool_mask_in_output_array = self._forward_translate_indices_to_bool_mask(invalidator_quib,
-                                                                                 working_component.component)
-        if np.any(bool_mask_in_output_array):
-            starting_path = [PathComponent(self.get_type(), bool_mask_in_output_array), *rest_of_path]
-        else:
-            starting_path = None
-
-        tuple_len = self._get_tuple_output_len()
-        if starting_path is None:
-            return []
-        elif tuple_len is None:
-            return [starting_path]
-        else:
-            return [[PathComponent(tuple, i), *starting_path] for i in range(tuple_len)]
-
-    def _backward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
-        quib_loop_shape = quib.get_shape().get_value()[:-self._get_arg_core_ndim(quib) or None]
-        result_bool_mask = self._get_bool_mask_representing_indices_in_result(indices)
-        result_core_axes = tuple(range(-1, -1 - self._get_result_core_ndim(), -1))
-        quib_result_loop_ndim_before_broadcast = result_bool_mask.ndim - len(quib_loop_shape) - len(result_core_axes)
-        broadcast_loop_dimensions_to_remove = tuple(range(0, quib_result_loop_ndim_before_broadcast))
-        reduced_bool_mask = np.any(result_bool_mask, axis=result_core_axes + broadcast_loop_dimensions_to_remove)
-        broadcast_loop_dimensions_to_reduce = tuple(i for i, (result_len, quib_len) in
-                                                    enumerate(zip(reduced_bool_mask.shape, quib_loop_shape))
-                                                    if result_len != quib_len)
-        reduced_bool_mask = np.any(reduced_bool_mask, axis=broadcast_loop_dimensions_to_reduce, keepdims=True)
-        return np.broadcast_to(reduced_bool_mask, quib_loop_shape)
-
-    def _get_source_path_in_quib(self, quib: Quib, filtered_path_in_result: List[PathComponent]):
-        if len(filtered_path_in_result) == 0 or filtered_path_in_result[0].indexed_cls is tuple:
-            return []
-        working_component, *rest_of_path = filtered_path_in_result
-        indices_in_data_source = self._backward_translate_indices_to_bool_mask(quib, working_component.component)
-        return [PathComponent(self.get_type(), indices_in_data_source)]
-
-    def _get_source_paths_of_quibs_given_path(self, filtered_path_in_result: List[PathComponent]):
-        return {quib: self._get_source_path_in_quib(quib, filtered_path_in_result)
-                for quib in self._get_data_source_quibs()}
-
-
-@dataclass
-class Arg:
-    name: str
-
-    def get_value(self, arg_dict: Dict[str, Any]) -> Any:
-        return arg_dict[self.name]
-
-
-@dataclass
-class ArgWithDefault(Arg):
-    default: Any
-
-    def get_value(self, arg_dict: Dict[str, Any]) -> Any:
-        return arg_dict.get(self.name, self.default)
-
-
-class AxisWiseGraphicsFunctionQuib(GraphicsFunctionQuib, IndicesTranslatorFunctionQuib):
-    """
-    Axiswise functions are functions that loop over an input array,
-    and perform a calculation on each sub-item.
-    The looping can be done on any combination of the dimensions of the input array.
-    The dimensions used for looping are called "loop dimensions", and the dimensions of the sub-items
-    passed to the underlying calculation are called "core dimensions".
-    In numpy terminology, when a function is applied "along an axis" or "over axes",
-    it means that those axes will be the core dimensions in the calculation.
-    For example, in the function np.sum(data, axis=x), x could be:
-     - None: the core dimensions are ()
-     - n (an integer): the core dimensions are (n,)
-     - t (a tuple of integers): the core dimensions are t
-
-    In axiswise function, we could also separate the result ndarrays into loop dimensions and core dimensions.
-    In reduction functions like min, max and sum, the loop dimensions will be the outer dimensions in the resulting
-     array.
-    In functions applied along an axis, like cumsum and apply_along_axis, the loop dimensions remains the same, but the
-    core dimension is potentially expanded to multiple dimensions, depending on the calculation result dimensions.
-    """
-    # Some numpy functions that allow not passing a keyword argument don't actually have a default value for it.
-    # The default value is an instance of _NoValueType, and it is passed to an underlying function that
-    # actually calculates the default value. In the future we could try and get the default value per
-    # data type and function, but right now we assume we can know the default in advance.
-    TRANSLATION_RELATED_ARGS: List[Arg]
-
-    @abstractmethod
-    def _forward_translate_bool_mask(self, args_dict, boolean_mask, invalidator_quib: Quib):
-        """
-        Do the actual forward index translation.
-        """
-
-    def _get_translation_related_arg_dict(self):
-        arg_dict = {key: val for key, val in self._get_args_values(include_defaults=True).arg_values_by_name.items()
-                    if not isinstance(val, np._globals._NoValueType)}
-        return {arg.name: arg.get_value(arg_dict) for arg in self.TRANSLATION_RELATED_ARGS}
-
-    def _forward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
-        source_bool_mask = self._get_source_shaped_bool_mask(quib, indices)
-        args_dict = self._get_translation_related_arg_dict()
-        return self._forward_translate_bool_mask(args_dict, source_bool_mask, quib)
-
-    @abstractmethod
-    def _backward_translate_bool_mask(self, args_dict, bool_mask, quib: Quib):
-        pass
-
-    def _backward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
-        result_bool_mask = self._get_bool_mask_representing_indices_in_result(indices)
-        args_dict = self._get_translation_related_arg_dict()
-        return self._backward_translate_bool_mask(args_dict, result_bool_mask, quib)
-
-    def _get_source_path_in_quib(self, quib: Quib, filtered_path_in_result: List[PathComponent]):
-        if len(filtered_path_in_result) == 0:
-            return []
-        working_component, *rest_of_path = filtered_path_in_result
-        indices_in_data_source = self._backward_translate_indices_to_bool_mask(quib, working_component.component)
-        return [PathComponent(self.get_type(), indices_in_data_source)]
-
-    def _get_source_paths_of_quibs_given_path(self, filtered_path_in_result: List[PathComponent]):
-        return {quib: self._get_source_path_in_quib(quib, filtered_path_in_result)
-                for quib in self._get_data_source_quibs()}
-
-
-class ReductionAxisWiseGraphicsFunctionQuib(AxisWiseGraphicsFunctionQuib):
-    SUPPORTED_FUNCTIONS = {
-        np.sum: SupportedFunction({0}),
-        np.min: SupportedFunction({0}),
-        np.amin: SupportedFunction({0}),
-        np.max: SupportedFunction({0}),
-        np.amax: SupportedFunction({0}),
-        np.mean: SupportedFunction({0}),
-        np.average: SupportedFunction({0}),
-        np.linalg.norm: SupportedFunction({0})
-    }
-    TRANSLATION_RELATED_ARGS = [Arg('axis'), ArgWithDefault('keepdims', False), ArgWithDefault('where', True)]
-
-    def _forward_translate_bool_mask(self, args_dict, boolean_mask, invalidator_quib: Quib):
-        """
-        Calculate forward index translation for reduction functions by reducing the boolean arrays
-        with the same reduction params.
-        """
-        return np.logical_or.reduce(boolean_mask, **args_dict)
-
-    def _backward_translate_bool_mask(self, args_dict, bool_mask, quib: Quib):
-        keepdims = args_dict['keepdims']
-        if not keepdims:
-            input_core_dims = args_dict['axis']
-            if input_core_dims is not None:
-                bool_mask = np.expand_dims(bool_mask, input_core_dims)
-        bool_mask = np.broadcast_to(bool_mask, quib.get_shape().get_value())
-        return np.logical_and(bool_mask, args_dict['where'])
-
-
-class AlongAxisGraphicsFunctionQuib(AxisWiseGraphicsFunctionQuib):
-    SUPPORTED_FUNCTIONS = {
-        np.apply_along_axis: SupportedFunction({2}),
-    }
-    TRANSLATION_RELATED_ARGS = [Arg('axis')]
-
-    def _get_expanded_dims(self, axis, result_shape, source_shape):
-        func_result_ndim = len(result_shape) - len(source_shape) + 1
-        assert func_result_ndim >= 0, func_result_ndim
-        return tuple(range(axis, axis + func_result_ndim) if axis >= 0 else
-                     range(axis, axis - func_result_ndim, -1))
-
-    def _forward_translate_bool_mask(self, args_dict, boolean_mask, invalidator_quib: Quib):
-        """
-        Calculate forward index translation for apply_along_axis by applying np.any on the boolean mask.
-        After that we expand and broadcast the reduced mask to match the actual result shape, which is dependent
-        on the applied function return type.
-        """
-        axis = args_dict.pop('axis')
-        result_shape = self.get_shape().get_value()
-        dims_to_expand = self._get_expanded_dims(axis, result_shape, invalidator_quib.get_shape().get_value())
-        applied = np.apply_along_axis(np.any, axis, boolean_mask, **args_dict)
-        expanded = np.expand_dims(applied, dims_to_expand)
-        broadcast = np.broadcast_to(expanded, result_shape)
-        return broadcast
-
-    def _backward_translate_bool_mask(self, args_dict, bool_mask, quib: Quib):
-        axis = args_dict.pop('axis')
-        source_shape = quib.get_shape().get_value()
-        expanded_dims = self._get_expanded_dims(axis, bool_mask.shape, source_shape)
-        mask = np.expand_dims(np.any(bool_mask, axis=expanded_dims), axis)
-        return np.broadcast_to(mask, source_shape)
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/__init__.py b/pyquibbler/quib/graphics/axiswise_function_quibs/__init__.py
new file mode 100644
index 000000000..6ecdeab42
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/__init__.py
@@ -0,0 +1,3 @@
+from .vectorize import VectorizeGraphicsFunctionQuib, QVectorize
+from .apply_along_axis_function_quib import ApplyAlongAxisGraphicsFunctionQuib
+from .axiswise_function_quib import ReductionAxisWiseGraphicsFunctionQuib
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/apply_along_axis_function_quib.py b/pyquibbler/quib/graphics/axiswise_function_quibs/apply_along_axis_function_quib.py
new file mode 100644
index 000000000..22a239c3f
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/apply_along_axis_function_quib.py
@@ -0,0 +1,194 @@
+from typing import Optional, List, Any, Callable, Tuple, Mapping
+
+import numpy as np
+from numpy import ndindex, s_
+
+from pyquibbler.quib.function_quibs.utils import ArgsValues
+from pyquibbler.quib.graphics.utils import remove_created_graphics
+from pyquibbler.quib.proxy_quib import ProxyQuib
+from pyquibbler.quib.assignment import PathComponent
+from pyquibbler.quib.quib import Quib
+from pyquibbler.quib.function_quibs.indices_translator_function_quib import SupportedFunction
+from pyquibbler.quib.graphics.axiswise_function_quibs.axiswise_function_quib import AxisWiseGraphicsFunctionQuib, Arg
+from pyquibbler.quib.quib import cache_method_until_full_invalidation
+
+
+class ApplyAlongAxisGraphicsFunctionQuib(AxisWiseGraphicsFunctionQuib):
+    SUPPORTED_FUNCTIONS = {
+        np.apply_along_axis: SupportedFunction({2}),
+    }
+    TRANSLATION_RELATED_ARGS = [Arg('axis')]
+
+    @classmethod
+    def create(cls, func, func_args=(), func_kwargs=None, cache_behavior=None, lazy=None, **init_kwargs):
+        func_kwargs = func_kwargs or {}
+        pass_quibs = func_kwargs.pop('pass_quibs', False)
+        return super(ApplyAlongAxisGraphicsFunctionQuib, cls).create(func=func, func_args=func_args,
+                                                                     cache_behavior=cache_behavior,
+                                                                     lazy=lazy,
+                                                                     func_kwargs=func_kwargs,
+                                                                     pass_quibs=pass_quibs,
+                                                                     **init_kwargs)
+
+    def _get_expanded_dims(self, axis, result_shape, source_shape):
+        func_result_ndim = len(result_shape) - len(source_shape) + 1
+        assert func_result_ndim >= 0, func_result_ndim
+        return tuple(range(axis, axis + func_result_ndim) if axis >= 0 else
+                     range(axis, axis - func_result_ndim, -1))
+
+    def _get_loop_shape(self) -> Tuple[int, ...]:
+        return tuple([s for i, s in enumerate(self.arr.get_shape()) if i != self.core_axis])
+
+    def _forward_translate_bool_mask(self, args_dict, boolean_mask, invalidator_quib: Quib):
+        """
+        Calculate forward index translation for apply_along_axis by applying np.any on the boolean mask.
+        After that we expand and broadcast the reduced mask to match the actual result shape, which is dependent
+        on the applied function return type.
+        """
+        axis = args_dict.pop('axis')
+        result_shape = self.get_shape()
+        dims_to_expand = self._get_expanded_dims(axis, result_shape, invalidator_quib.get_shape())
+        applied = np.apply_along_axis(np.any, axis, boolean_mask, **args_dict)
+        expanded = np.expand_dims(applied, dims_to_expand)
+        broadcast = np.broadcast_to(expanded, result_shape)
+        return broadcast
+
+    def _run_func1d(self, arr: np.ndarray, *args, **kwargs) -> Any:
+        """
+        Run the one dimensional function on args and kwargs and potentially get value of the result if it's a quib
+        """
+        res = self.func1d(arr, *args, **kwargs)
+        if isinstance(res, Quib):
+            return res.get_value()
+        return res
+
+    @cache_method_until_full_invalidation
+    def _get_sample_result(self) -> Any:
+        """
+        Get a result representing a possible return value of the user defined func1d.
+        In practice, we simply run on index 0 of every dimension and take the oned slice.
+        """
+        input_array_shape = self.arr.get_shape()
+        item = tuple([slice(None) if i == self.core_axis else 0 for i in range(len(input_array_shape))])
+
+        if self._pass_quibs:
+            input_array = ProxyQuib(self.arr)
+        else:
+            input_array = self.arr.get_value_valid_at_path([PathComponent(component=item, indexed_cls=np.ndarray)])
+
+        oned_slice = input_array[item]
+        args, kwargs = self._prepare_args_for_call(None)
+        args_values = ArgsValues.from_function_call(func=self.func, args=args, kwargs=kwargs, include_defaults=False)
+
+        with remove_created_graphics():
+            return self._run_func1d(
+                oned_slice,
+                *args_values.arg_values_by_name.get('args', []),
+                **args_values.arg_values_by_name.get('kwargs', {})
+            )
+
+    @cache_method_until_full_invalidation
+    def _get_invalid_value_at_correct_shape_and_dtype(self) -> np.ndarray:
+        """
+        Get a value that represents the real result in both shape and dtype.
+        Because the returned value can potentially be an ndarray, this will potentially execute the func1d once in
+        order to get the results shape (but will remove any artists created by running it).
+        """
+        input_array_shape = self.arr.get_shape()
+        sample_result_arr = np.asarray(self._get_sample_result())
+        dims_to_expand = list(range(0, self.core_axis))
+        dims_to_expand += list(range(-1, -(len(input_array_shape) - self.core_axis), -1))
+        expanded = np.expand_dims(sample_result_arr, dims_to_expand)
+        new_shape = [dim
+                     for i, d in enumerate(input_array_shape)
+                     for dim in ([d] if i != self.core_axis else sample_result_arr.shape)]
+        return np.array(np.broadcast_to(expanded, new_shape))
+
+    @property
+    def core_axis(self) -> int:
+        axis = self._get_args_values()['axis']
+        if isinstance(axis, Quib):
+            # since this is a parameter quib, we always need it completely valid in order to run anything
+            axis = axis.get_value()
+        return axis if axis >= 0 else len(self.arr.get_shape()) + axis
+
+    @property
+    def arr(self) -> Quib:
+        from pyquibbler import q
+        arr_ = self._get_args_values()['arr']
+        # ensure we're dealing with an ndarray- because we're not always running apply_along_axis which takes care of
+        # this for us (for example, when getting a sample result) we do this on any access to the array to ensure no
+        # issues if we were passed a list
+        return q(np.array, arr_)
+
+    @property
+    def func1d(self) -> Callable:
+        return self._get_args_values()['func1d']
+
+    def _get_oned_slice_for_running_func1d(self, indices: Tuple):
+        """
+        Get the proper slice as a parameter for the func- this depends on whether the user specified that he wants a
+        quib representing the result (and if so, we create a proxy quib so as not to invalidate inner quibs he creates)
+        or the values themselves of the slice
+        """
+        return ProxyQuib(self.arr[indices]) if self._pass_quibs else self.arr[indices].get_value()
+
+    def _get_result_at_indices(self,
+                               requested_indices_bool_mask: np.ndarray,
+                               indices_before_axis: Tuple,
+                               indices_after_axis: Tuple,
+                               func1d_args: Tuple[Any, ...],
+                               func1d_kwargs: Mapping[str, Any]):
+        """
+        Get a result at the indices given as arguments- this does not necessarily mean that func1d will be run; if the
+        bool mask is not True within the given indices, then the result of this iteration was not requested (at
+        valid_path), and so we will simply return the sample result as a placeholder
+        """
+        indices = indices_before_axis + s_[(...,)] + indices_after_axis
+        if np.any(requested_indices_bool_mask[indices]):
+            with self._call_func_context(self._graphics_collection_ndarr[indices_before_axis + indices_after_axis]):
+                res = self._run_func1d(self._get_oned_slice_for_running_func1d(indices),
+                                       *func1d_args,
+                                       **func1d_kwargs)
+        else:
+            res = self._get_sample_result()
+        return res
+
+    def _apply_along_axis(self, valid_path):
+        """
+        Run "apply_along_axis"- in reality, we need to map several different ndarrays, and so running apply_along_axis
+        itself would be problematic (as we need the indices themselves of the 1d slice). Given this, we simply run on
+        two loops of ndindex- the outer loop representing indices *before* the loop dimension, the inner loop
+        representing indices *after* the loop dimension. We then select everything in between the two index tuples,
+        which is a 1d slice.
+        """
+        indices = True if len(valid_path) == 0 else valid_path[0].component
+        ni, nk = self.arr.get_shape()[:self.core_axis], self.arr.get_shape()[self.core_axis + 1:]
+        out = self.get_value_valid_at_path(None)
+        args_values = ArgsValues.from_function_call(func=self.func, args=self.args, kwargs=self.kwargs,
+                                                    include_defaults=False)
+        args_by_name = args_values.arg_values_by_name
+        bool_mask = self._get_bool_mask_representing_indices_in_result(indices)
+        for ii in ndindex(ni):
+            for kk in ndindex(nk):
+                out[ii + s_[(...,)] + kk] = self._get_result_at_indices(bool_mask,
+                                                                        indices_before_axis=ii,
+                                                                        indices_after_axis=kk,
+                                                                        func1d_args=args_by_name.get('args', []),
+                                                                        func1d_kwargs=args_by_name.get('kwargs', {}))
+
+        return out
+
+    def _call_func(self, valid_path: Optional[List[PathComponent]]) -> Any:
+        if valid_path is None:
+            return self._get_invalid_value_at_correct_shape_and_dtype()
+
+        self._initialize_graphics_ndarr()
+        return self._apply_along_axis(valid_path)
+
+    def _backward_translate_bool_mask(self, args_dict, bool_mask, quib: Quib):
+        axis = args_dict.pop('axis')
+        source_shape = quib.get_shape()
+        expanded_dims = self._get_expanded_dims(axis, bool_mask.shape, source_shape)
+        mask = np.expand_dims(np.any(bool_mask, axis=expanded_dims), axis)
+        return np.broadcast_to(mask, source_shape)
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/axiswise_function_quib.py b/pyquibbler/quib/graphics/axiswise_function_quibs/axiswise_function_quib.py
new file mode 100644
index 000000000..fc295a3f5
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/axiswise_function_quib.py
@@ -0,0 +1,121 @@
+from __future__ import annotations
+import numpy as np
+from dataclasses import dataclass
+from abc import abstractmethod
+from typing import Any, Dict, List
+
+from pyquibbler.quib.quib import Quib
+
+from pyquibbler.quib.graphics.graphics_function_quib import GraphicsFunctionQuib
+from pyquibbler.quib.assignment import PathComponent
+from pyquibbler.quib.function_quibs.indices_translator_function_quib import IndicesTranslatorFunctionQuib, \
+    SupportedFunction
+
+
+@dataclass
+class Arg:
+    name: str
+
+    def get_value(self, arg_dict: Dict[str, Any]) -> Any:
+        return arg_dict[self.name]
+
+
+@dataclass
+class ArgWithDefault(Arg):
+    default: Any
+
+    def get_value(self, arg_dict: Dict[str, Any]) -> Any:
+        return arg_dict.get(self.name, self.default)
+
+
+class AxisWiseGraphicsFunctionQuib(GraphicsFunctionQuib, IndicesTranslatorFunctionQuib):
+    """
+    Axiswise functions are functions that loop over an input array,
+    and perform a calculation on each sub-item.
+    The looping can be done on any combination of the dimensions of the input array.
+    The dimensions used for looping are called "loop dimensions", and the dimensions of the sub-items
+    passed to the underlying calculation are called "core dimensions".
+    In numpy terminology, when a function is applied "along an axis" or "over axes",
+    it means that those axes will be the core dimensions in the calculation.
+    For example, in the function np.sum(data, axis=x), x could be:
+     - None: the core dimensions are ()
+     - n (an integer): the core dimensions are (n,)
+     - t (a tuple of integers): the core dimensions are t
+
+    In axiswise function, we could also separate the result ndarrays into loop dimensions and core dimensions.
+    In reduction functions like min, max and sum, the loop dimensions will be the outer dimensions in the resulting
+     array.
+    In functions applied along an axis, like cumsum and apply_along_axis, the loop dimensions remains the same, but the
+    core dimension is potentially expanded to multiple dimensions, depending on the calculation result dimensions.
+    """
+    # Some numpy functions that allow not passing a keyword argument don't actually have a default value for it.
+    # The default value is an instance of _NoValueType, and it is passed to an underlying function that
+    # actually calculates the default value. In the future we could try and get the default value per
+    # data type and function, but right now we assume we can know the default in advance.
+    TRANSLATION_RELATED_ARGS: List[Arg]
+
+    @abstractmethod
+    def _forward_translate_bool_mask(self, args_dict, boolean_mask, invalidator_quib: Quib):
+        """
+        Do the actual forward index translation.
+        """
+
+    def _get_translation_related_arg_dict(self):
+        arg_dict = {key: val for key, val in self._get_args_values(include_defaults=True).arg_values_by_name.items()
+                    if not isinstance(val, np._globals._NoValueType)}
+        return {arg.name: arg.get_value(arg_dict) for arg in self.TRANSLATION_RELATED_ARGS}
+
+    def _forward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
+        source_bool_mask = self._get_source_shaped_bool_mask(quib, indices)
+        args_dict = self._get_translation_related_arg_dict()
+        return self._forward_translate_bool_mask(args_dict, source_bool_mask, quib)
+
+    @abstractmethod
+    def _backward_translate_bool_mask(self, args_dict, bool_mask, quib: Quib):
+        pass
+
+    def _backward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
+        result_bool_mask = self._get_bool_mask_representing_indices_in_result(indices)
+        args_dict = self._get_translation_related_arg_dict()
+        return self._backward_translate_bool_mask(args_dict, result_bool_mask, quib)
+
+    def _get_source_path_in_quib(self, quib: Quib, filtered_path_in_result: List[PathComponent]):
+        if len(filtered_path_in_result) == 0:
+            return []
+        working_component, *rest_of_path = filtered_path_in_result
+        indices_in_data_source = self._backward_translate_indices_to_bool_mask(quib, working_component.component)
+        return [PathComponent(self.get_type(), indices_in_data_source)]
+
+    def _get_source_paths_of_quibs_given_path(self, filtered_path_in_result: List[PathComponent]):
+        return {quib: self._get_source_path_in_quib(quib, filtered_path_in_result)
+                for quib in self._get_data_source_quibs()}
+
+
+class ReductionAxisWiseGraphicsFunctionQuib(AxisWiseGraphicsFunctionQuib):
+    SUPPORTED_FUNCTIONS = {
+        np.sum: SupportedFunction({0}),
+        np.min: SupportedFunction({0}),
+        np.amin: SupportedFunction({0}),
+        np.max: SupportedFunction({0}),
+        np.amax: SupportedFunction({0}),
+        np.mean: SupportedFunction({0}),
+        np.average: SupportedFunction({0}),
+        np.linalg.norm: SupportedFunction({0})
+    }
+    TRANSLATION_RELATED_ARGS = [Arg('axis'), ArgWithDefault('keepdims', False), ArgWithDefault('where', True)]
+
+    def _forward_translate_bool_mask(self, args_dict, boolean_mask, invalidator_quib: Quib):
+        """
+        Calculate forward index translation for reduction functions by reducing the boolean arrays
+        with the same reduction params.
+        """
+        return np.logical_or.reduce(boolean_mask, **args_dict)
+
+    def _backward_translate_bool_mask(self, args_dict, bool_mask, quib: Quib):
+        keepdims = args_dict['keepdims']
+        if not keepdims:
+            input_core_dims = args_dict['axis']
+            if input_core_dims is not None:
+                bool_mask = np.expand_dims(bool_mask, input_core_dims)
+        bool_mask = np.broadcast_to(bool_mask, quib.get_shape())
+        return np.logical_and(bool_mask, args_dict['where'])
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/__init__.py b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/__init__.py
new file mode 100644
index 000000000..148eaa3fe
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/__init__.py
@@ -0,0 +1,1 @@
+from .vectorize_graphics_function_quib import VectorizeGraphicsFunctionQuib, QVectorize
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/utils.py b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/utils.py
new file mode 100644
index 000000000..876c71847
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/utils.py
@@ -0,0 +1,94 @@
+from __future__ import annotations
+import numpy as np
+from typing import Tuple, Iterable, Optional, Dict, Union, Any, Callable, TYPE_CHECKING
+from string import ascii_letters
+from itertools import islice, chain
+
+from pyquibbler.quib.function_quibs.indices_translator_function_quib import Args, Kwargs
+
+if TYPE_CHECKING:
+    from .vectorize_metadata import ArgsMetadata, ArgId
+
+Shape = Tuple[int, ...]
+
+
+def get_core_axes(core_ndim: int) -> Tuple[int, ...]:
+    """
+    Return the numbers of the core axes given the amount of core dimensions.
+    """
+    return tuple(range(-1, -1 - core_ndim, -1))
+
+
+def construct_core_dims_str(ndim: int, names: Iterable[str]) -> str:
+    """
+    Given an dimension amount and an iterable of usable dimension names,
+    return a core dimension specification as used in a numpy ufunc signature.
+    """
+    return f'({",".join(islice(names, ndim))})'
+
+
+def iter_arg_ids_and_values(args: Args, kwargs: Kwargs) -> Iterable[Tuple[ArgId, Any]]:
+    return chain(enumerate(args), kwargs.items())
+
+
+def construct_signature(args_core_ndims: Iterable[int], results_core_ndims: Iterable[int]) -> str:
+    """
+    Create a numpy ufunc signature given args core ndims and results core ndims.
+    """
+    names = iter(ascii_letters)
+    construct_core_dims_strs = lambda ndims: ",".join(construct_core_dims_str(ndim, names) for ndim in ndims)
+    return f'{construct_core_dims_strs(args_core_ndims)}->{construct_core_dims_strs(results_core_ndims)}'
+
+
+def alter_signature(args_metadata: ArgsMetadata, results_core_ndims: Iterable[int],
+                    arg_ids_to_new_core_ndims: Optional[Dict[Union[str, int], int]]) -> str:
+    """
+    Create a new ufuc signature according the current args metadata, result core ndims,
+    and a mapping between arg ids to their new core dimensions.
+    """
+    args_core_ndims = [arg_ids_to_new_core_ndims.pop(arg_id, arg_meta.core_ndim)
+                       for arg_id, arg_meta in args_metadata.items()]
+    assert not arg_ids_to_new_core_ndims, f'Invalid arg ids: {set(arg_ids_to_new_core_ndims.keys())}'
+    return construct_signature(args_core_ndims, results_core_ndims)
+
+
+def copy_vectorize(vectorize, func=None, otypes=None, excluded=None, signature=None) -> np.vectorize:
+    """
+    Copy a vectorize object while allowing to replace some attributes.
+    """
+    if func is None:
+        func = vectorize.pyfunc
+    if otypes is None:
+        otypes = vectorize.otypes
+    if excluded is None:
+        excluded = vectorize.excluded
+    if signature is None:
+        signature = vectorize.signature
+    return np.vectorize(func, otypes=otypes, doc=vectorize.__doc__, excluded=excluded, cache=vectorize.cache,
+                        signature=signature, pass_quibs=vectorize.pass_quibs)
+
+
+class Indices:
+    """
+    A small wrapper of indices, used to keep in an array without having core dimensions.
+    """
+    __slots__ = ['indices']
+
+    def __init__(self, indices):
+        self.indices = tuple(indices)
+
+
+def get_indices_array(shape: Shape) -> np.ndarray:
+    """
+    Given a shape, return an ndarray in which each cell holds an Indices object with
+    indices pointing to that cell.
+    """
+    return np.apply_along_axis(Indices, -1, np.moveaxis(np.indices(shape), 0, -1))
+
+
+def convert_args_and_kwargs(converter: Callable, args: Args, kwargs: Kwargs):
+    """
+    Apply the given converter on all given arg and kwarg values.
+    """
+    return (tuple(converter(i, val) for i, val in enumerate(args)),
+            {name: converter(name, val) for name, val in kwargs.items()})
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/vectorize_graphics_function_quib.py b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/vectorize_graphics_function_quib.py
new file mode 100644
index 000000000..eee69751b
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/vectorize_graphics_function_quib.py
@@ -0,0 +1,258 @@
+from __future__ import annotations
+import numpy as np
+from functools import cached_property
+from typing import Any, Optional, List, Tuple, Set
+
+from pyquibbler.quib.quib import Quib, cache_method_until_full_invalidation
+from pyquibbler.quib.proxy_quib import ProxyQuib
+from pyquibbler.quib.graphics.graphics_function_quib import GraphicsFunctionQuib
+from pyquibbler.quib.assignment import PathComponent
+from pyquibbler.quib.function_quibs.indices_translator_function_quib import IndicesTranslatorFunctionQuib, \
+    SupportedFunction
+from pyquibbler.quib.function_quibs.utils import ArgsValues
+from pyquibbler.quib.utils import copy_and_replace_quibs_with_vals
+
+from .utils import copy_vectorize, get_core_axes, get_indices_array, iter_arg_ids_and_values, alter_signature, \
+    convert_args_and_kwargs
+from .vectorize_metadata import VectorizeMetadata, ArgId, VectorizeCall
+from ...utils import remove_created_graphics
+
+
+def get_vectorize_call_data_args(args_values: ArgsValues) -> List[Any]:
+    """
+    Given a call to a vectorized function, return the arguments which act as data sources.
+    We are using args_values.args and args_values.kwargs instead of the full args dict on purpose,
+    to match vectorize function behavior.
+    """
+    vectorize, *args = args_values.args
+    return [val for key, val in iter_arg_ids_and_values(args, args_values.kwargs) if key not in vectorize.excluded]
+
+
+class VectorizeGraphicsFunctionQuib(GraphicsFunctionQuib, IndicesTranslatorFunctionQuib):
+    """
+    A Quib for wrapping np.vectorize.__call__.
+    Support forward and backwards translation, and partial calculation.
+    """
+    SUPPORTED_FUNCTIONS = {
+        np.vectorize.__call__: SupportedFunction(get_vectorize_call_data_args)
+    }
+
+    @classmethod
+    def _wrapper_call(cls, func, args, kwargs, **create_kwargs):
+        """
+        When create a quib function wrapper for np.vectorize.__call__, pass the QVectorize pass_quibs
+        attribute as a kwargs to create().
+        """
+        return super()._wrapper_call(func, args, kwargs, **create_kwargs, pass_quibs=args[0].pass_quibs)
+
+    @cached_property
+    def _vectorize(self) -> np.vectorize:
+        """
+        Get the vectorize object we were called with.
+        """
+        return self.args[0]
+
+    def _get_arg_ids_for_quib(self, quib: Quib) -> Set[ArgId]:
+        """
+        Given a parent quib, get all arg_ids it is given with.
+        For example, in the call `f(q, x=q)`, `_get_arg_ids_for_quib(q)` will return `{0, 'x'}`
+        """
+        return {arg_id for arg_id, arg in iter_arg_ids_and_values(self.args[1:], self.kwargs) if quib is arg}
+
+    def _wrap_vectorize_call_to_pass_quibs(self, call: VectorizeCall, args_metadata,
+                                           results_core_ndims) -> VectorizeCall:
+        """
+        Given a vectorize call and some metadata, return a new vectorize call that passes quibs to the vectorize pyfunc.
+        This is done by wrapping the original pyfunc. The quibs in args and kwargs are replaced with
+        indices arrays, and when the wrapper receives an index it uses it to get a GetItemQuib from the passed quib.
+        """
+        # We convert quibs to numpy arrays so we can slice them with tuples even if they are originally lists
+        quib_args = {arg_id: ProxyQuib(np.array(val)) for arg_id, val in iter_arg_ids_and_values(call.args, call.kwargs)
+                     if isinstance(val, Quib) and arg_id not in call.vectorize.excluded}
+
+        def convert_quibs_to_indices(arg_id, arg_val):
+            if arg_id in quib_args:
+                return get_indices_array(args_metadata[arg_id].loop_shape)
+            elif isinstance(arg_val, Quib):
+                return ProxyQuib(arg_val)
+            return arg_val
+
+        args, kwargs = convert_args_and_kwargs(convert_quibs_to_indices, call.args, call.kwargs)
+        # Indices arrays have 0 core dimensions, so if the signature is None, it just stays None.
+        # Otherwise, we need to construct a new signature in which the core dimensions of the quib
+        # args are zero.
+        # If indices core dimensions were one, then we couldn't keep using an empty signature when calling
+        # vectorize - and that means we would have needed to call the original function to get the tuple
+        # length - which is bad because we need to call it with quibs.
+        # To solve that we make sure our core dimensions are always 0 by using the Indices class.
+        signature = None if call.vectorize.signature is None else \
+            alter_signature(args_metadata, results_core_ndims, {arg_id: 0 for arg_id in quib_args})
+
+        def convert_indices_to_quibs(arg_id, arg_val):
+            quib = quib_args.get(arg_id)
+            if quib is None:
+                return arg_val
+            return quib[arg_val.indices]
+
+        def wrapper(*args, **kwargs):
+            args, kwargs = convert_args_and_kwargs(convert_indices_to_quibs, args, kwargs)
+            result = call.vectorize.pyfunc(*args, **kwargs)
+            return copy_and_replace_quibs_with_vals(result)
+
+        return VectorizeCall(copy_vectorize(call.vectorize, func=wrapper, signature=signature), args, kwargs)
+
+    def _get_vectorize_call(self, args_metadata, results_core_ndims, valid_path) -> VectorizeCall:
+        """
+        Get a VectorizeCall object to actually call vectorize with.
+        If asked to pass quibs, return a modified call that passes quibs.
+        """
+        if self._pass_quibs:
+            (vectorize, *args), kwargs = self.args, self.kwargs
+            call = VectorizeCall(vectorize, args, kwargs)
+            call = self._wrap_vectorize_call_to_pass_quibs(call, args_metadata, results_core_ndims)
+        else:
+            (vectorize, *args), kwargs = self._prepare_args_for_call(valid_path)
+            call = VectorizeCall(vectorize, args, kwargs)
+        return call
+
+    def _get_sample_result(self, args_metadata, results_core_ndims):
+        """
+        Run the vectorized function on sample arguments to get s single sample result.
+        Remove any graphics that were created during the call.
+        """
+        call = self._get_vectorize_call(args_metadata, results_core_ndims, None)
+        with remove_created_graphics():
+            return call.get_sample_result(args_metadata)
+
+    @property
+    @cache_method_until_full_invalidation
+    def _vectorize_metadata(self) -> VectorizeMetadata:
+        """
+        Get and cache metadata for the vectorize call.
+        """
+        (vectorize, *args), kwargs = self._prepare_args_for_call(None)
+        return VectorizeCall(vectorize, args, kwargs).get_metadata(self._get_sample_result)
+
+    def _forward_translate_indices_to_bool_mask(self, invalidator_quib: Quib, indices: Any) -> Any:
+        """
+        Translate source bool mask to result bool mask by calling any on the source core dimensions,
+        and broadcasting the result to the loop shape.
+        """
+        vectorize_metadata = self._vectorize_metadata
+        source_bool_mask = self._get_source_shaped_bool_mask(invalidator_quib, indices)
+        core_ndim = max(vectorize_metadata.args_metadata[arg_id].core_ndim
+                        for arg_id in self._get_arg_ids_for_quib(invalidator_quib))
+        source_bool_mask = np.any(source_bool_mask, axis=get_core_axes(core_ndim))
+        return np.broadcast_to(source_bool_mask, vectorize_metadata.result_loop_shape)
+
+    def _forward_translate_invalidation_path(self, invalidator_quib: Quib,
+                                             path: List[PathComponent]) -> List[Optional[List[PathComponent]]]:
+        working_component, *rest_of_path = path
+        bool_mask_in_output_array = self._forward_translate_indices_to_bool_mask(invalidator_quib,
+                                                                                 working_component.component)
+        if not np.any(bool_mask_in_output_array):
+            return []
+        starting_path = [PathComponent(self.get_type(), bool_mask_in_output_array), *rest_of_path]
+
+        vectorize_metadata = self._vectorize_metadata
+        if vectorize_metadata.is_result_a_tuple:
+            return [[PathComponent(tuple, i), *starting_path] for i in range(vectorize_metadata.tuple_length)]
+        else:
+            return [starting_path]
+
+    def _backward_translate_indices_to_bool_mask(self, quib: Quib, indices: Any) -> Any:
+        """
+        Translate indices in result backwards to indices in a data source quib, by calling any on result
+        core dimensions and de-broadcasting the loop dimensions into the argument loop dimension.
+        """
+        vectorize_metadata = self._vectorize_metadata
+        quib_arg_id = self._get_arg_ids_for_quib(quib).pop()
+        quib_loop_shape = vectorize_metadata.args_metadata[quib_arg_id].loop_shape
+        result_bool_mask = self._get_bool_mask_representing_indices_in_result(indices)
+        result_core_axes = vectorize_metadata.result_core_axes
+        new_broadcast_ndim = result_bool_mask.ndim - len(quib_loop_shape) - len(result_core_axes)
+        assert new_broadcast_ndim >= 0
+        new_broadcast_axes = tuple(range(0, new_broadcast_ndim))
+        reduced_bool_mask = np.any(result_bool_mask, axis=result_core_axes + new_broadcast_axes)
+        broadcast_loop_dimensions_to_reduce = tuple(i for i, (result_len, quib_len) in
+                                                    enumerate(zip(reduced_bool_mask.shape, quib_loop_shape))
+                                                    if result_len != quib_len)
+        reduced_bool_mask = np.any(reduced_bool_mask, axis=broadcast_loop_dimensions_to_reduce, keepdims=True)
+        return np.broadcast_to(reduced_bool_mask, quib_loop_shape)
+
+    def _get_source_path_in_quib(self, quib: Quib, filtered_path_in_result: List[PathComponent]):
+        """
+        Given a path in the result, return the path in the given quib on which the result path depends.
+        """
+        if len(filtered_path_in_result) == 0 or filtered_path_in_result[0].indexed_cls is tuple:
+            return []
+        working_component, *rest_of_path = filtered_path_in_result
+        indices_in_data_source = self._backward_translate_indices_to_bool_mask(quib, working_component.component)
+        return [PathComponent(self.get_type(), indices_in_data_source)]
+
+    def _get_source_paths_of_quibs_given_path(self, filtered_path_in_result: List[PathComponent]):
+        return {quib: self._get_source_path_in_quib(quib, filtered_path_in_result)
+                for quib in self._get_data_source_quibs()}
+
+    def _get_loop_shape(self) -> Tuple[int, ...]:
+        if self._vectorize_metadata.is_result_a_tuple:
+            return ()
+        return self._vectorize_metadata.result_loop_shape
+
+    def _wrap_vectorize_call_to_calc_only_needed(self, call: VectorizeCall, valid_path, otypes):
+        """
+        1. Create a bool mask with the same shape as the broadcast loop dimensions
+        2. Create a wrapper for the pyfunc that receives the same arguments, and an additional boolean that
+           instructs it to run the user function or return an empty result.
+        3. Vectorize the wrapper and run it with the boolean mask in addition to the original arguments
+        """
+        valid_indices = True if len(valid_path) == 0 else valid_path[0].component
+        bool_mask = np.any(self._get_bool_mask_representing_indices_in_result(valid_indices),
+                           axis=self._vectorize_metadata.result_core_axes)
+        # The logic in this wrapper should be as minimal as possible (including attribute access, etc.)
+        # as it is run for every index in the loop.
+        pyfunc = call.vectorize.pyfunc
+        empty_result = np.zeros(self._vectorize_metadata.result_core_shape, dtype=self._vectorize_metadata.result_dtype)
+
+        def wrapper(graphics_collection, should_run, *args, **kwargs):
+            if should_run:
+                with self._call_func_context(graphics_collection):
+                    return pyfunc(*args, **kwargs)
+            return empty_result
+
+        args_to_add = (self._graphics_collection_ndarr, bool_mask)
+        wrapper_excluded = {i + len(args_to_add) if isinstance(i, int) else i for i in self._vectorize.excluded}
+        wrapper_signature = call.vectorize.signature if call.vectorize.signature is None \
+            else '(),' * len(args_to_add) + call.vectorize.signature
+        vectorize = copy_vectorize(call.vectorize, func=wrapper, excluded=wrapper_excluded, signature=wrapper_signature,
+                                   otypes=otypes)
+        return VectorizeCall(vectorize, (*args_to_add, *call.args), call.kwargs)
+
+    def _call_func(self, valid_path: Optional[List[PathComponent]]):
+        """
+        Call vectoriz, but call the internal function only when required by the given valid_path.
+        """
+        self._initialize_graphics_ndarr()
+        vectorize_metadata = self._vectorize_metadata
+        if vectorize_metadata.is_result_a_tuple:
+            return super()._call_func(valid_path)
+        if valid_path is None:
+            return np.zeros(vectorize_metadata.result_shape, dtype=vectorize_metadata.result_dtype)
+        call = self._get_vectorize_call(vectorize_metadata.args_metadata,
+                                        vectorize_metadata.result_or_results_core_ndims, valid_path)
+        call = self._wrap_vectorize_call_to_calc_only_needed(call, valid_path, vectorize_metadata.otypes)
+        return call()
+
+
+class QVectorize(np.vectorize):
+    """
+    A small wrapper to the np.vectorize class, adding options to __init__ and wrapping __call__
+    with a quib function wrapper.
+    """
+
+    def __init__(self, *args, pass_quibs=False, signature=None, cache=False, **kwargs):
+        # We don't need the underlying vectorize object to cache, we are doing that ourselves.
+        super().__init__(*args, signature=signature, cache=False, **kwargs)
+        self.pass_quibs = pass_quibs
+
+    __call__ = VectorizeGraphicsFunctionQuib.create_wrapper(np.vectorize.__call__)
diff --git a/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/vectorize_metadata.py b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/vectorize_metadata.py
new file mode 100644
index 000000000..24bf730de
--- /dev/null
+++ b/pyquibbler/quib/graphics/axiswise_function_quibs/vectorize/vectorize_metadata.py
@@ -0,0 +1,244 @@
+from __future__ import annotations
+
+from functools import partial
+
+import numpy as np
+from dataclasses import dataclass
+from typing import Any, Dict, Optional, List, Tuple, Union, Callable
+
+from pyquibbler.quib.function_quibs.indices_translator_function_quib import Args, Kwargs
+
+from .utils import Shape, get_core_axes, iter_arg_ids_and_values, convert_args_and_kwargs
+
+ArgId = Union[int, str]
+ArgsMetadata: Dict[ArgId, VectorizeArgMetadata]
+
+
+@dataclass
+class VectorizeArgMetadata:
+    """
+    Holds metadata regarding core and loop dimensions of an argument to vectorize.
+    """
+    shape: Shape
+    core_shape: Shape
+    loop_shape: Shape
+
+    @classmethod
+    def from_arg_and_core_ndim(cls, arg: Any, core_ndim: int) -> VectorizeArgMetadata:
+        """
+        Given an argument and its core dimension amount, create and return a VectorizeArgMetadata object.
+        """
+        shape = np.shape(arg)
+        assert len(shape) >= core_ndim
+        if core_ndim == 0:
+            return cls(shape, (), shape)
+        return cls(shape, shape[-core_ndim:], shape[:-core_ndim])
+
+    @property
+    def core_ndim(self) -> int:
+        return len(self.core_shape)
+
+    @property
+    def loop_ndim(self) -> int:
+        return len(self.loop_shape)
+
+
+@dataclass
+class VectorizeMetadata:
+    """
+    Metadata regarding a vectorize call.
+    Depending on the arguments to vectorize, we might know or not know some of the metadata
+    in advance.
+    Sometimes, to get all necessary details, we must call the vectorize pyfunc to get a sample result.
+    This class does so lazily, by providing properties for all relevant metadata details that get a sample
+    result and update the metadata if needed before returning.
+    """
+    _get_sample_result: Callable
+
+    args_metadata: ArgsMetadata
+    result_loop_shape: Shape
+
+    _is_result_a_tuple: Optional[bool]
+    _results_core_ndims: Optional[List[int]]
+    _results_dtypes: Optional[List[np.dtype]]
+    _results_core_shapes: Optional[List[Shape]] = None
+
+    def _run_sample_and_update_metadata(self):
+        """
+        We need some piece of information we can get only by calling the pyfunc -
+        this function calls the pyfunc using the given _get_sample_result and updates all
+        relevant metadata, while asserting new details don't contradict previously known ones.
+        """
+        sample_result = self._get_sample_result(self.args_metadata, self._results_core_ndims)
+
+        # update _is_result_a_tuple
+        is_tuple = isinstance(sample_result, tuple)
+        if self._is_result_a_tuple is None:
+            self._is_result_a_tuple = is_tuple
+        else:
+            assert self._is_result_a_tuple == is_tuple
+
+        # update _results_core_shapes
+        assert self._results_core_shapes is None
+        self._results_core_shapes = list(map(np.shape, sample_result)) if is_tuple else [np.shape(sample_result)]
+
+        # update _results_core_ndims
+        results_core_ndims = list(map(len, self._results_core_shapes))
+        if self._results_core_ndims is None:
+            self._results_core_ndims = results_core_ndims
+        else:
+            assert self._results_core_ndims == results_core_ndims
+
+        # update _results_dtypes
+        results_dtypes = [np.asarray(result).dtype for result in (sample_result if is_tuple else (sample_result,))]
+        if self._results_dtypes is None:
+            self._results_dtypes = results_dtypes
+        else:
+            assert self._results_dtypes == results_dtypes
+
+    @property
+    def is_result_a_tuple(self) -> bool:
+        if self._is_result_a_tuple is None:
+            self._run_sample_and_update_metadata()
+        return self._is_result_a_tuple
+
+    @property
+    def results_core_ndims(self) -> List[int]:
+        if self._results_core_ndims is None:
+            self._run_sample_and_update_metadata()
+        assert self._is_result_a_tuple is True
+        return self._results_core_ndims
+
+    @property
+    def result_core_ndim(self) -> int:
+        if self._results_core_ndims is None:
+            self._run_sample_and_update_metadata()
+        assert self._is_result_a_tuple is False
+        return self._results_core_ndims[0]
+
+    @property
+    def result_or_results_core_ndims(self) -> List[int]:
+        if self._results_core_ndims is None:
+            self._run_sample_and_update_metadata()
+        return self._results_core_ndims
+
+    @property
+    def result_core_axes(self) -> Tuple[int, ...]:
+        return get_core_axes(self.result_core_ndim)
+
+    @property
+    def result_core_shape(self) -> Shape:
+        if self.result_core_ndim == 0:
+            return ()
+        if self._results_core_shapes is None:
+            self._run_sample_and_update_metadata()
+        assert not self.is_result_a_tuple
+        return self._results_core_shapes[0]
+
+    @property
+    def result_shape(self) -> Shape:
+        return self.result_loop_shape + self.result_core_shape
+
+    @property
+    def result_dtype(self) -> np.dtype:
+        if self._results_dtypes is None:
+            self._run_sample_and_update_metadata()
+        assert not self.is_result_a_tuple
+        return self._results_dtypes[0]
+
+    @property
+    def results_dtypes(self) -> List[np.dtype]:
+        if self._results_dtypes is None:
+            self._run_sample_and_update_metadata()
+        assert self.is_result_a_tuple
+        return self._results_dtypes
+
+    @property
+    def otypes(self) -> str:
+        if self._results_dtypes is None:
+            self._run_sample_and_update_metadata()
+        return ''.join(dtype.char for dtype in self._results_dtypes)
+
+    @property
+    def tuple_length(self):
+        """
+        Assuming vectorize returns a tuple, return its length.
+        """
+        return len(self.results_dtypes)
+
+
+@dataclass
+class VectorizeCall:
+    vectorize: np.vectorize
+    args: Args
+    kwargs: Kwargs
+
+    @staticmethod
+    def get_sample_arg_core(args_metadata: ArgsMetadata, arg_id: Union[str, int], arg_value: Any) -> Any:
+        """
+        Get a sample core value from an array to call a non-vectorized function with.
+        """
+        meta = args_metadata.get(arg_id)
+        if meta is None:
+            return arg_value
+        # We should only use the loop shape and not the core shape, as the core shape changes with pass_quibs=True
+        return np.asarray(arg_value)[(0,) * meta.loop_ndim]
+
+    def get_sample_result(self, args_metadata: ArgsMetadata) -> Any:
+        """
+        Get one sample result from the inner function of a vectorize
+        """
+        args, kwargs = convert_args_and_kwargs(partial(self.get_sample_arg_core, args_metadata), self.args, self.kwargs)
+        return self.vectorize.pyfunc(*args, **kwargs)
+
+    def __call__(self):
+        # If we pass quibs to the wrapper, we will create a new graphics quib, so we use the original vectorize
+        return np.vectorize.__overridden__.__call__(self.vectorize, *self.args, **self.kwargs)
+
+    def _get_args_and_results_core_ndims(self):
+        """
+        Get the args and results core dimensions in a vectorize call.
+        """
+        num_not_excluded = len((set(range(len(self.args))) | set(self.kwargs)) - self.vectorize.excluded)
+        if self.vectorize._in_and_out_core_dims is None:
+            args_core_ndims = [0] * num_not_excluded
+            results_core_ndims = None
+            is_tuple = None
+        else:
+            in_core_dims, out_core_dims = self.vectorize._in_and_out_core_dims
+            args_core_ndims = list(map(len, in_core_dims))
+            assert len(args_core_ndims) == num_not_excluded
+            is_tuple = len(out_core_dims) > 1
+            results_core_ndims = list(map(len, out_core_dims))
+        return args_core_ndims, results_core_ndims, is_tuple
+
+    def _get_args_metadata(self, args_core_ndims) -> ArgsMetadata:
+        """
+        Get the args metadata in a vectorize call.
+        """
+        arg_index = 0
+        args_metadata = {}
+        for i, arg in iter_arg_ids_and_values(self.args, self.kwargs):
+            if i not in self.vectorize.excluded:
+                args_metadata[i] = VectorizeArgMetadata.from_arg_and_core_ndim(arg, args_core_ndims[arg_index])
+                arg_index += 1
+        assert arg_index == len(args_core_ndims)
+        return args_metadata
+
+    def get_metadata(self, sample_result_callback=None):
+        """
+        Create a metadata object from a vectorize call.
+        If sample_result_callback is given, use it to get a sample result instead of calling vectorize.pyfunc
+        directly.
+        """
+        args_core_ndims, results_core_ndims, is_tuple = self._get_args_and_results_core_ndims()
+        args_metadata = self._get_args_metadata(args_core_ndims)
+        # Calculate the result shape like done in np.function_base._parse_input_dimensions
+        dummy_arrays = [np.lib.stride_tricks.as_strided(0, arg_metadata.loop_shape)
+                        for arg_metadata in args_metadata.values()]
+        result_loop_shape = np.lib.stride_tricks._broadcast_shape(*dummy_arrays)
+        results_dtypes = None if self.vectorize.otypes is None else list(map(np.dtype, self.vectorize.otypes))
+        if sample_result_callback is None:
+            sample_result_callback = lambda args_metadata, results_core_ndims: self.get_sample_result(args_metadata)
+        return VectorizeMetadata(sample_result_callback, args_metadata, result_loop_shape, is_tuple, results_core_ndims,
+                                 results_dtypes)
diff --git a/pyquibbler/quib/graphics/event_handling/plot_inverse_assigner.py b/pyquibbler/quib/graphics/event_handling/plot_inverse_assigner.py
index a30bea0ca..bb69c661b 100644
--- a/pyquibbler/quib/graphics/event_handling/plot_inverse_assigner.py
+++ b/pyquibbler/quib/graphics/event_handling/plot_inverse_assigner.py
@@ -63,7 +63,7 @@ def get_quibs_with_assignments_for_axes(args: List[Any],
         if isinstance(quib, Quib):
             # Support indexing of lists when more than one marker is dragged
             for data_index in data_indices:
-                shape = quib.get_shape().get_value()
+                shape = quib.get_shape()
                 if len(shape) == 0:
                     path = []
                 elif len(shape) == 1:
diff --git a/pyquibbler/quib/graphics/global_collecting.py b/pyquibbler/quib/graphics/global_collecting.py
index 192c2d2e9..22a3d86e6 100644
--- a/pyquibbler/quib/graphics/global_collecting.py
+++ b/pyquibbler/quib/graphics/global_collecting.py
@@ -2,10 +2,10 @@
 import contextlib
 import functools
 import threading
-from typing import List, Callable
+from abc import ABC
+from typing import List, Callable, Type
 from matplotlib.artist import Artist
-
-from pyquibbler.exceptions import PyQuibblerException
+from matplotlib.widgets import AxesWidget
 
 ORIGINAL_ARTIST_INIT = Artist.__init__
 
@@ -14,57 +14,51 @@
 OVERRIDDEN_GRAPHICS_FUNCTIONS_RUNNING = 0
 
 
-def is_within_artists_collector():
-    """
-    Returns whether there's an active artists collector at this moment
-    """
-    return GLOBAL_ARTIST_COLLECTORS > 0
-
-
-class AlreadyCollectingArtistsException(PyQuibblerException):
-    pass
-
-
-class ArtistsCollector:
+class GraphicsCollector(ABC):
     """
     A context manager to begin global collecting of artists.
     After exiting the context manager use `get_artists_collected` to get a list of artists created during this time
     """
 
-    def __init__(self, artists_collected: List[Artist] = None, thread_id: int = None,
-                 raise_if_within_collector: bool = False, previous_init: Callable = None):
-        self._artists_collected = artists_collected or []
+    cls: Type = NotImplemented
+
+    def __init__(self, objects_collected: List = None, thread_id: int = None,
+                 previous_init: Callable = None):
+        self._objects_collected = objects_collected or []
         self._thread_id = thread_id or threading.get_ident()
-        self._raise_if_within_collector = raise_if_within_collector
         self._previous_init = previous_init
 
     @property
-    def artists_collected(self):
-        return self._artists_collected
+    def objects_collected(self):
+        return self._objects_collected
 
-    def wrap_artist_creation(self, func):
+    def wrap_object_creation(self, func):
         @functools.wraps(func)
-        def wrapped_init(artist, *args, **kwargs):
-            res = func(artist, *args, **kwargs)
+        def wrapped_init(obj, *args, **kwargs):
+            res = func(obj, *args, **kwargs)
             if threading.get_ident() == self._thread_id and OVERRIDDEN_GRAPHICS_FUNCTIONS_RUNNING:
-                self._artists_collected.append(artist)
+                self._objects_collected.append(obj)
             return res
 
         return wrapped_init
 
     def __enter__(self):
-        global GLOBAL_ARTIST_COLLECTORS
-        if self._raise_if_within_collector and GLOBAL_ARTIST_COLLECTORS:
-            raise AlreadyCollectingArtistsException()
-        GLOBAL_ARTIST_COLLECTORS += 1
-        self._previous_init = Artist.__init__
-        Artist.__init__ = self.wrap_artist_creation(Artist.__init__)
+        self._previous_init = self.cls.__init__
+        self.cls.__init__ = self.wrap_object_creation(self.cls.__init__)
         return self
 
     def __exit__(self, *_, **__):
-        global GLOBAL_ARTIST_COLLECTORS
-        GLOBAL_ARTIST_COLLECTORS -= 1
-        Artist.__init__ = self._previous_init
+        self.cls.__init__ = self._previous_init
+
+
+class ArtistsCollector(GraphicsCollector):
+
+    cls = Artist
+
+
+class AxesWidgetsCollector(GraphicsCollector):
+
+    cls = AxesWidget
 
 
 @contextlib.contextmanager
diff --git a/pyquibbler/quib/graphics/graphics_collection.py b/pyquibbler/quib/graphics/graphics_collection.py
new file mode 100644
index 000000000..26a747189
--- /dev/null
+++ b/pyquibbler/quib/graphics/graphics_collection.py
@@ -0,0 +1,8 @@
+from dataclasses import dataclass, field
+from typing import Set
+
+
+@dataclass
+class GraphicsCollection:
+    widgets: Set = field(default_factory=set)
+    artists: Set = field(default_factory=set)
diff --git a/pyquibbler/quib/graphics/graphics_function_quib.py b/pyquibbler/quib/graphics/graphics_function_quib.py
index c21a83f3a..094b1eca1 100644
--- a/pyquibbler/quib/graphics/graphics_function_quib.py
+++ b/pyquibbler/quib/graphics/graphics_function_quib.py
@@ -1,33 +1,27 @@
-from typing import List, Callable, Tuple, Any, Mapping, Dict, Optional, Iterable
+import numpy as np
+from contextlib import contextmanager
+from typing import List, Callable, Tuple, Any, Mapping, Dict, Optional, Iterable, Set
 from matplotlib.artist import Artist
 from matplotlib.axes import Axes
-from matplotlib.collections import Collection
-from matplotlib.image import AxesImage
-from matplotlib.lines import Line2D
-from matplotlib.patches import Patch
-from matplotlib.spines import Spine
-from matplotlib.table import Table
-from matplotlib.text import Text
-
-from . import global_collecting
-from .event_handling import CanvasEventHandler
+from matplotlib.widgets import AxesWidget
+
+from .global_collecting import ArtistsCollector, AxesWidgetsCollector
+from .graphics_collection import GraphicsCollection
 from .quib_guard import QuibGuard
+from .utils import save_func_and_args_on_artists, get_axeses_to_array_names_to_starting_indices_and_artists, \
+    remove_artist, get_axeses_to_array_names_to_artists, get_artist_array, ArrayNameToArtists, track_artist
 from ..assignment import AssignmentTemplate, PathComponent
 from ..function_quibs import DefaultFunctionQuib, CacheBehavior
-from ..utils import iter_object_type_in_args, recursively_run_func_on_object
+from ..utils import recursively_run_func_on_object, iter_object_type_in_args, iter_quibs_in_args
+from ...env import GRAPHICS_LAZY
 
 
-def save_func_and_args_on_artists(artists: List[Artist], func: Callable, args: Iterable[Any]):
-    """
-    Set the drawing func and on the artists- this will be used later on when tracking artists, in order to know how to
-    inverse and handle events.
-    If there's already a creating func, we assume the lower func that already created the artist is the actual
-    drawing func (such as a user func that called plt.plot)
-    """
-    for artist in artists:
-        if getattr(artist, '_quibbler_drawing_func', None) is None:
-            artist._quibbler_drawing_func = func
-            artist._quibbler_args = args
+def proxify_args(args, kwargs):
+    from pyquibbler.quib import Quib, ProxyQuib
+    replace_quibs_with_proxy_quibs = lambda arg: ProxyQuib(arg) if isinstance(arg, Quib) else arg
+    args = recursively_run_func_on_object(replace_quibs_with_proxy_quibs, args)
+    kwargs = {k: recursively_run_func_on_object(replace_quibs_with_proxy_quibs, v) for k, v in kwargs.items()}
+    return args, kwargs
 
 
 class GraphicsFunctionQuib(DefaultFunctionQuib):
@@ -40,16 +34,6 @@ class GraphicsFunctionQuib(DefaultFunctionQuib):
     # Avoid unnecessary repaints
     _DEFAULT_CACHE_BEHAVIOR = CacheBehavior.ON
 
-    TYPES_TO_ARTIST_ARRAY_NAMES = {
-        Line2D: "lines",
-        Collection: "collections",
-        Patch: "patches",
-        Text: "texts",
-        Spine: "spines",
-        Table: "tables",
-        AxesImage: "images"
-    }
-
     # Some attributes, like 'color', are chosen by matplotlib automatically if not specified.
     # Therefore, if these attributes were not specified, we need to copy them
     # from old artists to new artists.
@@ -59,6 +43,7 @@ class GraphicsFunctionQuib(DefaultFunctionQuib):
         '_color': {'color'},
         '_facecolor': {'facecolor'},
     }
+
     # Note that this current implementation does not account for implicit color specification,
     # such as with plot(x,y,c), where c = iquib('r').
     # also it does not account for artiststs created within user-defined functions.
@@ -68,57 +53,36 @@ def __init__(self,
                  args: Tuple[Any, ...],
                  kwargs: Mapping[str, Any],
                  cache_behavior: Optional[CacheBehavior],
-                 artists: List[Artist],
                  assignment_template: Optional[AssignmentTemplate] = None,
                  had_artists_on_last_run: bool = False,
-                 receive_quibs: bool = False
-                 ):
+                 pass_quibs: bool = False):
         super().__init__(func, args, kwargs, cache_behavior, assignment_template=assignment_template)
-        self._artists = artists
         self._had_artists_on_last_run = had_artists_on_last_run
-        self._receive_quibs = receive_quibs
+        self._pass_quibs = pass_quibs
+        self._graphics_collection_ndarr = None
 
     @classmethod
-    def create(cls, func, func_args=(), func_kwargs=None, cache_behavior=None, lazy=False, **kwargs):
-        self = super().create(func, func_args, func_kwargs, cache_behavior, artists=[], **kwargs)
+    def create(cls, func, func_args=(), func_kwargs=None, cache_behavior=None, lazy=None, **init_kwargs):
+        self = super().create(func, func_args, func_kwargs, cache_behavior, **init_kwargs)
+        if lazy is None:
+            lazy = GRAPHICS_LAZY
         if not lazy:
             self.get_value()
         return self
 
-    def persist_self_on_artists(self):
+    def persist_self_on_artists(self, artists: Set[Artist]):
         """
         Persist self on all artists we're connected to, making sure we won't be garbage collected until they are
         off the screen
         We need to also go over args as there may be a situation in which the function did not create new artists, but
         did perform an action on an existing one, such as Axes.set_xlim
         """
-        artists_to_persist_on = self._artists if len(self._artists) > 0 else iter_object_type_in_args(Artist, self.args,
-                                                                                                      self.kwargs)
+        artists_to_persist_on = artists if artists else iter_object_type_in_args(Artist, self.args, self.kwargs)
         for artist in artists_to_persist_on:
             quibs = getattr(artist, '_quibbler_graphics_function_quibs', set())
             quibs.add(self)
             artist._quibbler_graphics_function_quibs = quibs
 
-    def _get_artist_array_name(self, artist: Artist):
-        return next((
-            array_name
-            for type_, array_name in self.TYPES_TO_ARTIST_ARRAY_NAMES.items()
-            if isinstance(artist, type_)
-        ), "artists")
-
-    def _get_artist_array(self, artist: Artist):
-        return getattr(artist.axes, self._get_artist_array_name(artist))
-
-    def _get_axeses_to_array_names_to_artists(self) -> Dict[Axes, Dict[str, List[Artist]]]:
-        """
-        Creates a mapping of axes -> artist_array_name (e.g. `lines`) -> artists
-        """
-        axeses_to_array_names_to_artists = {}
-        for artist in self._artists:
-            array_names_to_artists = axeses_to_array_names_to_artists.setdefault(artist.axes, {})
-            array_names_to_artists.setdefault(self._get_artist_array_name(artist), []).append(artist)
-        return axeses_to_array_names_to_artists
-
     def _update_position_and_attributes_of_created_artists(
             self,
             axes: Axes,
@@ -141,8 +105,8 @@ def _update_position_and_attributes_of_created_artists(
             for previous_artist, new_artist in zip(previous_artists, new_artists):
                 for attribute in self.ATTRIBUTES_TO_COPY_FROM_ARTIST_TO_ARTIST_UNLESS_SPECIFED.keys():
                     if hasattr(previous_artist, attribute) \
-                        and not (self.kwargs.keys() &
-                                 self.ATTRIBUTES_TO_COPY_FROM_ARTIST_TO_ARTIST_UNLESS_SPECIFED[attribute]):
+                            and not (self.kwargs.keys() &
+                                     self.ATTRIBUTES_TO_COPY_FROM_ARTIST_TO_ARTIST_UNLESS_SPECIFED[attribute]):
                         setattr(new_artist, attribute, getattr(previous_artist, attribute))
 
     def _update_new_artists_from_previous_artists(self,
@@ -150,7 +114,7 @@ def _update_new_artists_from_previous_artists(self,
                                                       Axes, Dict[str, Tuple[int, List[Artist]]]
                                                   ],
                                                   current_axeses_to_array_names_to_artists: Dict[
-                                                      Axes, Dict[str, List[Artist]]
+                                                      Axes, ArrayNameToArtists
                                                   ]):
         """
         Updates the positions and attributes of the new artists from old ones
@@ -169,120 +133,104 @@ def _update_new_artists_from_previous_artists(self,
                 # If the array name isn't in previous_array_names_to_indices_and_artists,
                 # we don't need to update positions, etc
 
-    def save_drawing_func_on_artists(self):
-        """
-        Set the drawing func on the artists- this will be used later on when tracking artists, in order to know how to
-        inverse and handle events.
-        If there's already a creating func, we assume the lower func that already created the artist is the actual
-        drawing func (such as a user func that called plt.plot)
-        """
-        for artist in self._artists:
-            if getattr(artist, '_quibbler_drawing_func', None) is None:
-                artist._quibbler_drawing_func = self._func
-
-    def _run_function_on_quibs(self):
-        proxy_quibs = set()
-
-        def _replace_quibs_with_proxy_quibs(arg):
-            from pyquibbler.quib import Quib
-            from pyquibbler.quib.proxy_quib import ProxyQuib
-            if isinstance(arg, Quib):
-                proxy_quib = ProxyQuib(quib=arg)
-                proxy_quibs.add(proxy_quib)
-                return proxy_quib
-            return arg
-
-        args = recursively_run_func_on_object(_replace_quibs_with_proxy_quibs, self.args)
-        kwargs = {k: recursively_run_func_on_object(_replace_quibs_with_proxy_quibs, v) for k, v in self.kwargs.items()}
-        with QuibGuard(proxy_quibs):
-            return self.func(*args, **kwargs)
+    def _remove_artists(self, artist_set: Set[Artist]):
+        for artist in artist_set:
+            remove_artist(artist)
+        artist_set.clear()
+
+    def _disable_widgets(self, widget_set: Set[AxesWidget], always_disable=False):
+        from pyquibbler.quib.graphics.widgets import QRectangleSelector
+
+        for widget in widget_set:
+            if isinstance(widget, QRectangleSelector):
+                if not always_disable and widget.event_is_relevant_to_current_selector():
+                    widget.set_should_deactivate_after_release()
+                    continue
+                else:
+                    widget.set_active(False)
+                    widget.set_visible(False)
+            else:
+                widget.set_active(False)
 
-    def _create_new_artists(self,
-                            valid_path,
-                            previous_axeses_to_array_names_to_indices_and_artists:
-                            Dict[Axes, Dict[str, Tuple[int, List[Artist]]]] = None):
-        """
-        Create the new artists, then update them (if appropriate) with correct attributes (such as color) and place them
-        in the same place they were in their axes
-        """
+        widget_set.clear()
+
+    @contextmanager
+    def _handle_new_graphics_collection(self, graphics_collection: GraphicsCollection):
+        self._remove_artists_that_were_removed_from_axes(graphics_collection.artists)
+        # Get the *current* artists together with their starting indices (per axes per artists array) so we can
+        # place the new artists we create in their correct locations
         previous_axeses_to_array_names_to_indices_and_artists = \
-            previous_axeses_to_array_names_to_indices_and_artists or {}
+            get_axeses_to_array_names_to_starting_indices_and_artists(graphics_collection.artists)
+        self._remove_artists(graphics_collection.artists)
+        self._disable_widgets(graphics_collection.widgets)
 
-        with global_collecting.ArtistsCollector() as collector:
-            if self._receive_quibs:
-                func_res = self._run_function_on_quibs()
-            else:
-                with QuibGuard(self.parents):
-                    func_res = super()._call_func(valid_path)
+        with ArtistsCollector() as artists_collector, AxesWidgetsCollector() as widgets_collector:
+            yield
 
-        self._artists = collector.artists_collected
-        self._had_artists_on_last_run = len(self._artists) > 0
+        graphics_collection.artists.update(artists_collector.objects_collected)
+        graphics_collection.widgets.update(widgets_collector.objects_collected)
 
-        save_func_and_args_on_artists(self._artists, func=self.func, args=self.args)
-        self.persist_self_on_artists()
-        self.track_artists()
+        self._had_artists_on_last_run = len(graphics_collection.artists) > 0
 
-        current_axeses_to_array_names_to_artists = self._get_axeses_to_array_names_to_artists()
+        for artist in graphics_collection.artists:
+            save_func_and_args_on_artists(artist, func=self.func, args=self.args)
+            track_artist(artist)
+        self.persist_self_on_artists(graphics_collection.artists)
+
+        current_axeses_to_array_names_to_artists = get_axeses_to_array_names_to_artists(graphics_collection.artists)
         self._update_new_artists_from_previous_artists(previous_axeses_to_array_names_to_indices_and_artists,
                                                        current_axeses_to_array_names_to_artists)
-        return func_res
 
-    def _remove_current_artists(self):
+    def _remove_artists_that_were_removed_from_axes(self, artist_set: Set[Artist]):
         """
-        Remove all the current artists (does NOT redraw)
+        Remove any artists that we created that were removed by another means other than us (for example, cla())
         """
-        for artist in self._artists:
-            self._get_artist_array(artist).remove(artist)
+        artist_set.difference_update([artist for artist in artist_set if artist not in get_artist_array(artist)])
 
-    def _get_axeses_to_array_names_to_starting_indices_and_artists(self) -> Dict[
-        Axes, Dict[str, Tuple[int, List[Artist]]]
-    ]:
-        """
-        Creates a mapping of axes -> artists_array_name -> (starting_index, artists)
-        """
-        axeses_to_array_names_to_indices_and_artists = {}
-        for axes, array_names_to_artists in self._get_axeses_to_array_names_to_artists().items():
-            array_names_to_indices_and_artists = {}
-            axeses_to_array_names_to_indices_and_artists[axes] = array_names_to_indices_and_artists
+    def _iter_artist_sets(self) -> Iterable[Set[Artist]]:
+        return [] if self._graphics_collection_ndarr is None else map(lambda g: g.artists,
+                                                                      self._graphics_collection_ndarr.flat)
 
-            for array_name, artists in array_names_to_artists.items():
-                exemplifying_artist = artists[0]
-                array = getattr(exemplifying_artist.axes, array_name)
-                array_names_to_indices_and_artists[array_name] = (array.index(exemplifying_artist), artists)
+    def _iter_widget_sets(self) -> Iterable[Set[AxesWidget]]:
+        return [] if self._graphics_collection_ndarr is None else map(lambda g: g.widgets,
+                                                                      self._graphics_collection_ndarr.flat)
 
-        return axeses_to_array_names_to_indices_and_artists
+    def _iter_artists(self) -> Iterable[Artist]:
+        return (artist for artists in self._iter_artist_sets() for artist in artists)
 
-    def get_axeses(self) -> List[Axes]:
-        return list(self._get_axeses_to_array_names_to_artists().keys())
+    def get_axeses(self) -> Set[Axes]:
+        return {artist.axes for artist in self._iter_artists()}
 
-    def track_artists(self):
-        for artist in self._artists:
-            CanvasEventHandler.get_or_create_initialized_event_handler(artist.figure.canvas)
+    def _get_loop_shape(self) -> Tuple[int, ...]:
+        return ()
 
-    def _remove_artists_from_self_that_were_removed_from_axes(self):
-        """
-        We want to remove any artists that we created from self._artists that were removed by another means
-        other than us (for example, cla())
-        """
-        for artist in self._artists[:]:
-            if artist not in self._get_artist_array(artist):
-                self._artists.remove(artist)
+    def _initialize_graphics_ndarr(self):
+        loop_shape = self._get_loop_shape()
+        if self._graphics_collection_ndarr is not None and self._graphics_collection_ndarr.shape != loop_shape:
+            for artist_set in self._iter_artist_sets():
+                self._remove_artists(artist_set)
+            for widget_set in self._iter_widget_sets():
+                self._disable_widgets(widget_set, always_disable=True)
+            self._graphics_collection_ndarr = None
+        if self._graphics_collection_ndarr is None:
+            self._graphics_collection_ndarr = np.vectorize(lambda _: GraphicsCollection())(np.empty(loop_shape))
+
+    @contextmanager
+    def _call_func_context(self, graphics_collection):
+        with self._handle_new_graphics_collection(graphics_collection):
+            with QuibGuard(set(iter_quibs_in_args(self.args, self.kwargs))):
+                yield
 
     def _call_func(self, valid_path: Optional[List[PathComponent]]) -> Any:
         """
-        The main entrypoint- this reruns the function that created the artists in the first place,
-        and replaces the current artists with the new ones
-        :param valid_path:
-        """
-        self._remove_artists_from_self_that_were_removed_from_axes()
-        if self._had_artists_on_last_run and len(self._artists) == 0 and len(self.children) == 0:
-            # If we both do not have artists and do not have children, we do not want to run, as in this situation
-            # our artists have been cleared (by someone other than us)
-            # and nothing is actually pointing to us- we could have just as well been garbage collected
-            return
-
-        # Get the *current* artists together with their starting indices (per axes per artists array) so we can
-        # place the new artists we create in their correct locations
-        axeses_to_array_names_to_indices_and_artists = self._get_axeses_to_array_names_to_starting_indices_and_artists()
-        self._remove_current_artists()
-        return self._create_new_artists(valid_path, axeses_to_array_names_to_indices_and_artists)
+        Create the new artists, then update them (if appropriate) with correct attributes (such as color) and place them
+        in the same place they were in their axes.
+        Return the function's result.
+        """
+        self._initialize_graphics_ndarr()
+        # This implementation does not support partial calculation
+        assert self._graphics_collection_ndarr.ndim == 0
+        args, kwargs = proxify_args(self.args, self.kwargs) if self._pass_quibs \
+            else self._prepare_args_for_call(valid_path)
+        with self._call_func_context(self._graphics_collection_ndarr[()]):
+            return self.func(*args, **kwargs)
diff --git a/pyquibbler/quib/graphics/plot_graphics_function_quib.py b/pyquibbler/quib/graphics/plot_graphics_function_quib.py
index f783df576..c6c78c778 100644
--- a/pyquibbler/quib/graphics/plot_graphics_function_quib.py
+++ b/pyquibbler/quib/graphics/plot_graphics_function_quib.py
@@ -6,8 +6,8 @@ class PlotGraphicsFunctionQuib(GraphicsFunctionQuib):
     Made to add the _index_in_plot attribute to plot artists.
     """
 
-    def _create_new_artists(self, *args, **kwargs):
-        func_res = super()._create_new_artists(*args, **kwargs)
+    def _call_func(self, *args, **kwargs):
+        func_res = super()._call_func(*args, **kwargs)
         for i, artist in enumerate(func_res):
             artist._index_in_plot = i
         return func_res
diff --git a/pyquibbler/quib/graphics/quib_guard.py b/pyquibbler/quib/graphics/quib_guard.py
index b3298472a..38138ea81 100644
--- a/pyquibbler/quib/graphics/quib_guard.py
+++ b/pyquibbler/quib/graphics/quib_guard.py
@@ -1,10 +1,13 @@
 from __future__ import annotations
-
 import contextlib
-from typing import Set
+from dataclasses import dataclass
+from typing import Set, TYPE_CHECKING
 
 from pyquibbler.exceptions import PyQuibblerException
 
+if TYPE_CHECKING:
+    from pyquibbler.quib import Quib
+
 QUIB_GUARDS = []
 
 
@@ -12,14 +15,16 @@ def is_within_quib_guard():
     return len(QUIB_GUARDS) > 0
 
 
-def get_current_quib_guard():
+def get_current_quib_guard() -> QuibGuard:
     return QUIB_GUARDS[-1]
 
 
+@dataclass
 class CannotAccessQuibInScopeException(PyQuibblerException):
+    quib: Quib
 
     def __str__(self):
-        return "You cannot access quibs from the global scope within a user function"
+        return f"Illegal access to {self.quib}. Note that access to quibs from the global scope is not allowed"
 
 
 class AnotherQuibGuardIsAlreadyActiveException(PyQuibblerException):
@@ -32,6 +37,7 @@ class QuibGuard:
     A quib guard allows you to specify places in which only certain quibs (and their recursive children)
     are allowed to be accessed
     """
+
     def __init__(self, quibs_allowed: Set, allowed_quib_get_values_count: int = 0):
         self._quibs_allowed = quibs_allowed
         self._allowed_quib_get_values_count = allowed_quib_get_values_count
@@ -44,10 +50,13 @@ def __enter__(self):
     def get_value_context_manager(self, quib):
         all_allowed_quibs = set().union(*[{quib, *quib._get_children_recursively()} for quib in self._quibs_allowed])
         if self._allowed_quib_get_values_count == 0 and quib not in all_allowed_quibs:
-            raise CannotAccessQuibInScopeException()
+            raise CannotAccessQuibInScopeException(quib)
         self._allowed_quib_get_values_count += 1
         yield
         self._allowed_quib_get_values_count -= 1
 
+    def add_allowed_quib(self, quib: Quib):
+        self._quibs_allowed.add(quib)
+
     def __exit__(self, exc_type, exc_val, exc_tb):
         QUIB_GUARDS.pop()
diff --git a/pyquibbler/quib/graphics/redraw.py b/pyquibbler/quib/graphics/redraw.py
index 5c453028f..f21c6c699 100644
--- a/pyquibbler/quib/graphics/redraw.py
+++ b/pyquibbler/quib/graphics/redraw.py
@@ -1,4 +1,6 @@
 import contextlib
+from typing import Set
+
 from matplotlib.axes import Axes
 
 from pyquibbler.logger import logger
@@ -18,17 +20,18 @@ def aggregate_redraw_mode():
     IN_AGGREGATE_REDRAW_MODE = True
     yield
     IN_AGGREGATE_REDRAW_MODE = False
-    for axes in AXESES_TO_REDRAW:
-        redraw_axes(axes)
+    redraw_axeses(AXESES_TO_REDRAW)
     AXESES_TO_REDRAW.clear()
 
 
-def redraw_axes(axes: Axes, force: bool = False):
+def redraw_axeses(axeses: Set[Axes], force: bool = False):
     """
     Actual redrawing of axes- this should be WITHOUT rendering anything except for the new artists
     """
     if IN_AGGREGATE_REDRAW_MODE and not force:
-        AXESES_TO_REDRAW.add(axes)
+        AXESES_TO_REDRAW.update(axeses)
     else:
+        canvases = {axes.figure.canvas for axes in axeses}
         with timer(name="redraw", callback=lambda t: logger.info(f"redraw canvas {t}")):
-            axes.figure.canvas.draw()
+            for canvas in canvases:
+                canvas.draw()
diff --git a/pyquibbler/quib/graphics/replacing_graphics_function_quib.py b/pyquibbler/quib/graphics/replacing_graphics_function_quib.py
index 55d9bf60b..8a807b190 100644
--- a/pyquibbler/quib/graphics/replacing_graphics_function_quib.py
+++ b/pyquibbler/quib/graphics/replacing_graphics_function_quib.py
@@ -1,5 +1,5 @@
 from itertools import chain
-from typing import Optional
+from typing import Optional, Set
 from matplotlib.artist import Artist
 
 from pyquibbler.quib import GraphicsFunctionQuib
@@ -15,8 +15,8 @@ def _remove_previous_quib_from_parents(self, previous_quib: 'ReplacingGraphicsFu
         for parent in previous_quib.parents:
             parent.remove_child(previous_quib)
 
-    def persist_self_on_artists(self):
-        for artist in chain(self._artists, iter_object_type_in_args(Artist, self.args, self.kwargs)):
+    def persist_self_on_artists(self, artists: Set[Artist]):
+        for artist in chain(artists, iter_object_type_in_args(Artist, self.args, self.kwargs)):
             name = f'_quibbler_{self.func.__name__}'
             current_quib: Optional['ReplacingGraphicsFunctionQuib'] = getattr(artist, name, None)
             if current_quib is not self and current_quib is not None:
diff --git a/pyquibbler/quib/graphics/utils.py b/pyquibbler/quib/graphics/utils.py
new file mode 100644
index 000000000..babf210e8
--- /dev/null
+++ b/pyquibbler/quib/graphics/utils.py
@@ -0,0 +1,103 @@
+from typing import List, Callable, Tuple, Any, Dict, Iterable, Set
+from contextlib import contextmanager
+
+from matplotlib.axes import Axes
+from matplotlib.artist import Artist
+from matplotlib.collections import Collection
+from matplotlib.image import AxesImage
+from matplotlib.lines import Line2D
+from matplotlib.patches import Patch
+from matplotlib.spines import Spine
+from matplotlib.table import Table
+from matplotlib.text import Text
+
+from .event_handling import CanvasEventHandler
+from .global_collecting import ArtistsCollector, AxesWidgetsCollector
+
+ArrayNameToArtists = Dict[str, List[Artist]]
+
+TYPES_TO_ARTIST_ARRAY_NAMES = {
+    Line2D: "lines",
+    Collection: "collections",
+    Patch: "patches",
+    Text: "texts",
+    Spine: "spines",
+    Table: "tables",
+    AxesImage: "images"
+}
+
+
+def get_artist_array_name(artist: Artist):
+    return next((
+        array_name
+        for type_, array_name in TYPES_TO_ARTIST_ARRAY_NAMES.items()
+        if isinstance(artist, type_)
+    ), "artists")
+
+
+def get_artist_array(artist: Artist):
+    return getattr(artist.axes, get_artist_array_name(artist))
+
+
+def get_axeses_to_array_names_to_artists(artists: Iterable[Artist]) -> Dict[Axes, ArrayNameToArtists]:
+    """
+    Creates a mapping of axes -> artist_array_name (e.g. `lines`) -> artists
+    """
+    axeses_to_array_names_to_artists = {}
+    for artist in artists:
+        array_names_to_artists = axeses_to_array_names_to_artists.setdefault(artist.axes, {})
+        array_names_to_artists.setdefault(get_artist_array_name(artist), []).append(artist)
+    return axeses_to_array_names_to_artists
+
+
+def remove_artist(artist: Artist):
+    """
+    Remove an artist from its artist array (does NOT redraw)
+    """
+    get_artist_array(artist).remove(artist)
+
+
+def get_axeses_to_array_names_to_starting_indices_and_artists(artists: Set[Artist]) \
+        -> Dict[Axes, Dict[str, Tuple[int, List[Artist]]]]:
+    """
+    Creates a mapping of axes -> artists_array_name -> (starting_index, artists)
+    """
+    axeses_to_array_names_to_indices_and_artists = {}
+    for axes, array_names_to_artists in get_axeses_to_array_names_to_artists(artists).items():
+        array_names_to_indices_and_artists = {}
+        axeses_to_array_names_to_indices_and_artists[axes] = array_names_to_indices_and_artists
+
+        for array_name, array_artists in array_names_to_artists.items():
+            exemplifying_artist = array_artists[0]
+            array = getattr(exemplifying_artist.axes, array_name)
+            array_names_to_indices_and_artists[array_name] = (array.index(exemplifying_artist), array_artists)
+
+    return axeses_to_array_names_to_indices_and_artists
+
+
+def save_func_and_args_on_artists(artist: Artist, func: Callable, args: Iterable[Any]):
+    """
+    Set the drawing func and on the artist- this will be used later on when tracking an artist, in order to know how to
+    inverse and handle events.
+    If there's already a creating func, we assume the lower func that already created the artist is the actual
+    drawing func (such as a user func that called plt.plot)
+    """
+    if getattr(artist, '_quibbler_drawing_func', None) is None:
+        artist._quibbler_drawing_func = func
+        artist._quibbler_args = args
+
+
+def track_artist(artist: Artist):
+    CanvasEventHandler.get_or_create_initialized_event_handler(artist.figure.canvas)
+
+
+@contextmanager
+def remove_created_graphics():
+    with ArtistsCollector() as collector, AxesWidgetsCollector() as widgets_collector:
+        yield
+    for artist in collector.objects_collected:
+        remove_artist(artist)
+
+    for widget in widgets_collector.objects_collected:
+        widget.set_active(False)
+        widget.set_visible(False)
diff --git a/pyquibbler/quib/graphics/widgets/rectangle_selector.py b/pyquibbler/quib/graphics/widgets/rectangle_selector.py
index d8b314e6a..718f135d4 100644
--- a/pyquibbler/quib/graphics/widgets/rectangle_selector.py
+++ b/pyquibbler/quib/graphics/widgets/rectangle_selector.py
@@ -30,32 +30,45 @@ def unlock(self):
 class QRectangleSelector(RectangleSelector):
     CURRENT_SELECTOR = Locked(Mutable(None))
 
-    def __init__(self, ax, onselect=None, extents=None, interactive=None, **kwargs):
+    def __init__(self, ax, onselect=None, extents=None, allow_resize=True, interactive=None, **kwargs):
         self.changed_callback = None
         if interactive is None:
             interactive = True
         if onselect is None:
             onselect = lambda *args, **kwargs: None
         super().__init__(ax, onselect, interactive=interactive, **kwargs)
+        self.allow_resize = allow_resize
         if extents is not None:
             self.extents = extents
+        self._should_deactivate_after_release = False
+
+    def set_should_deactivate_after_release(self):
+        self._should_deactivate_after_release = True
+
+    def is_current_event_a_move_event(self):
+        return 'move' in self.state or self.active_handle == 'C'
 
     def event_is_relevant_to_current_selector(self) -> bool:
-        return (self.active_handle and self.active_handle != 'C') or (
-            (('move' in self.state or self.active_handle == 'C')))
+        return (self.active_handle and self.active_handle != 'C') or self.is_current_event_a_move_event()
 
     def _onmove(self, event):
         if self.event_is_relevant_to_current_selector():
-            with self.CURRENT_SELECTOR.unlock() as current_selector:
-                if current_selector.val is None or current_selector.val is self:
-                    current_selector.val = self
-                    return super()._onmove(event)
+            if self.allow_resize or self.is_current_event_a_move_event():
+                with self.CURRENT_SELECTOR.unlock() as current_selector:
+                    if current_selector.val is None or current_selector.val is self:
+                        current_selector.val = self
+                        return super()._onmove(event)
 
     def _release(self, event):
         with self.CURRENT_SELECTOR.unlock() as current_selector:
             if self.event_is_relevant_to_current_selector() and current_selector.val is self:
                 release_result = super()._release(event)
                 current_selector.val = None
+
+                if self._should_deactivate_after_release:
+                    self.set_active(False)
+                    self.set_visible(False)
+
                 return release_result
 
     @property
@@ -76,10 +89,31 @@ class RectangleSelectorGraphicsFunctionQuib(WidgetGraphicsFunctionQuib):
     """
     WIDGET_CLS = RectangleSelector
 
+    def _widget_is_attempting_to_resize_when_not_allowed(self, extents):
+        """
+        There is a bug in the matplotlib widget that causes it to sometimes give incorrect extents,
+        and attempt to resize even when the user did not request to resize- we here check if the widget attempted to
+        resize when it should not have been able to
+        """
+        init_val = self._get_args_values().get('extents')
+        allow_resize = self._get_args_values().get('allow_resize')
+        if isinstance(init_val, Quib):
+            previous_value = init_val.get_value()
+        else:
+            previous_value = self.get_value().extents
+
+        return not allow_resize and (
+                previous_value[1] - previous_value[0] != extents[1] - extents[0] or
+                previous_value[3] - previous_value[2] != extents[3] - extents[2]
+        )
+
     def _on_changed(self, extents):
         init_val = self._get_args_values().get('extents')
+
         with timer("selector_change", lambda x: logger.info(f"selector change {x}")):
             if isinstance(init_val, Quib):
+                if self._widget_is_attempting_to_resize_when_not_allowed(extents):
+                    return
                 init_val[:] = extents
             else:
                 # We only need to invalidate children if we didn't assign
diff --git a/pyquibbler/quib/graphics/widgets/widget_graphics_function_quib.py b/pyquibbler/quib/graphics/widgets/widget_graphics_function_quib.py
index 89acd8580..ac2335875 100644
--- a/pyquibbler/quib/graphics/widgets/widget_graphics_function_quib.py
+++ b/pyquibbler/quib/graphics/widgets/widget_graphics_function_quib.py
@@ -5,6 +5,7 @@
 
 
 class WidgetGraphicsFunctionQuib(GraphicsFunctionQuib):
+
     WIDGET_CLS: Type[Widget]
 
     def get_axeses(self):
diff --git a/pyquibbler/quib/input_quib.py b/pyquibbler/quib/input_quib.py
index a6e5e7f20..599c0f8f9 100644
--- a/pyquibbler/quib/input_quib.py
+++ b/pyquibbler/quib/input_quib.py
@@ -31,6 +31,11 @@ def __init__(self, value: Any, assignment_template: Optional[AssignmentTemplate]
             if is_there_a_quib_in_object(value, force_recursive=True):
                 raise CannotNestQuibInIQuibException(self)
 
+        from .graphics.quib_guard import is_within_quib_guard, get_current_quib_guard
+        if is_within_quib_guard():
+            quib_guard = get_current_quib_guard()
+            quib_guard.add_allowed_quib(self)
+
     def _get_inner_value_valid_at_path(self, path: List[PathComponent]) -> Any:
         """
         No need to do any calculation, this is an input quib.
diff --git a/pyquibbler/quib/override_choice/choice_cache.py b/pyquibbler/quib/override_choice/choice_cache.py
deleted file mode 100644
index 22d0d0cf9..000000000
--- a/pyquibbler/quib/override_choice/choice_cache.py
+++ /dev/null
@@ -1,44 +0,0 @@
-from __future__ import annotations
-from typing import TYPE_CHECKING, Optional, Dict, Any
-
-from .override_dialog import OverrideChoice
-
-if TYPE_CHECKING:
-    from pyquibbler.quib import Quib
-    from .override_choice import OverrideOptionsTree
-
-
-class ChoiceCache:
-    """
-    Implements caching for user override choices.
-    A choice can be reused in the context of a new assignment,
-    if the inversed quib (the quib that was originally assigned to) is the same, and the options are the same.
-    The override options for a inversed quib might change because when different indices are assigned to, the
-    inversion tree might change (for example if one of the quib's parents is concat).
-    """
-
-    def __init__(self):
-        self._map: Dict[Any, OverrideChoice] = {}
-
-    def _get_key(self, inversed_quib: Quib, options_tree: OverrideOptionsTree):
-        """
-        Returns a key for cache lookup.
-        """
-        return (inversed_quib,
-                tuple(option.override.quib for option in options_tree.options),
-                options_tree.can_diverge)
-
-    def store(self, choice: OverrideChoice, inversed_quib: Quib, options_tree: OverrideOptionsTree):
-        """
-        Store the given user choice according to the inversed quib and the options.
-        """
-        self._map[self._get_key(inversed_quib, options_tree)] = choice
-
-    def load(self, inversed_quib: Quib, options_tree: OverrideOptionsTree) -> Optional[OverrideChoice]:
-        """
-        Try to find a cached choice for the given inversed quib and override options.
-        """
-        return self._map.get(self._get_key(inversed_quib, options_tree))
-
-    def clear(self):
-        self._map.clear()
diff --git a/pyquibbler/quib/override_choice/choice_context.py b/pyquibbler/quib/override_choice/choice_context.py
new file mode 100644
index 000000000..39ae194c2
--- /dev/null
+++ b/pyquibbler/quib/override_choice/choice_context.py
@@ -0,0 +1,21 @@
+from __future__ import annotations
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Tuple
+
+if TYPE_CHECKING:
+    from pyquibbler.quib import Quib
+
+
+@dataclass(frozen=True)
+class ChoiceContext:
+    """
+    The context in which a choice was taken.
+    Implements __hash__ and __eq__ to allow caching of user override choices.
+    A choice can be reused in the context of a new assignment,
+    if the inversed quib (the quib that was originally assigned to) is the same, and the options are the same.
+    The override options for a inversed quib might change because when different indices are assigned to, the
+    inversion tree might change (for example if one of the quib's parents is concat).
+    """
+    inversed_quib: Quib
+    quibs_available_for_override: Tuple[Quib, ...]
+    can_diverge: bool
diff --git a/pyquibbler/quib/override_choice/override_choice.py b/pyquibbler/quib/override_choice/override_choice.py
index e646a7bdf..ac68fbeca 100644
--- a/pyquibbler/quib/override_choice/override_choice.py
+++ b/pyquibbler/quib/override_choice/override_choice.py
@@ -1,14 +1,14 @@
 from __future__ import annotations
 from dataclasses import dataclass
-from typing import List, TYPE_CHECKING, Optional
+from typing import List, TYPE_CHECKING, Optional, Dict, ClassVar
 
 from pyquibbler.exceptions import PyQuibblerException
 from pyquibbler.quib.assignment import QuibWithAssignment, CannotReverseException, \
     Assignment
 
+from .choice_context import ChoiceContext
 from .override_dialog import OverrideChoiceType, OverrideChoice, choose_override_dialog
 from .types import OverrideRemoval, OverrideGroup, OverrideWithOverrideRemovals
-from .choice_cache import ChoiceCache
 from ...env import ASSIGNMENT_RESTRICTIONS
 
 if TYPE_CHECKING:
@@ -47,7 +47,7 @@ class OverrideOptionsTree:
     - override_removals: all the override removals that will have to be applied if we
                          diverge and use children overrides
     """
-    _CHOICE_CACHE = ChoiceCache()
+    _CHOICE_CACHE: ClassVar[Dict[ChoiceContext, OverrideChoice]] = {}
 
     inversed_quib: FunctionQuib
     options: List[OverrideWithOverrideRemovals]
@@ -76,17 +76,23 @@ def __bool__(self):
         """
         return len(self.options) > 0 or self.can_diverge
 
+    @property
+    def choice_context(self):
+        return ChoiceContext(self.inversed_quib,
+                             tuple(option.override.quib for option in self.options),
+                             self.can_diverge)
+
     def _try_load_choice(self) -> Optional[OverrideChoice]:
         """
         If a choice fitting the current options has been cached, return it. Otherwise return None.
         """
-        return self._CHOICE_CACHE.load(self.inversed_quib, self)
+        return self._CHOICE_CACHE.get(self.choice_context)
 
     def _store_choice(self, choice: OverrideChoice):
         """
         Store a user override choice in the cache for future use.
         """
-        self._CHOICE_CACHE.store(choice, self.inversed_quib, self)
+        self._CHOICE_CACHE[self.choice_context] = choice
 
     def get_override_choice(self) -> OverrideChoice:
         """
diff --git a/pyquibbler/quib/proxy_quib.py b/pyquibbler/quib/proxy_quib.py
index 7d8e0a53a..65898f781 100644
--- a/pyquibbler/quib/proxy_quib.py
+++ b/pyquibbler/quib/proxy_quib.py
@@ -1,8 +1,11 @@
-from typing import Set, Optional, List, Any
+from __future__ import annotations
+from typing import Set, Optional, List, Any, TYPE_CHECKING
 
-from pyquibbler import Assignment
-from pyquibbler.quib import Quib
-from pyquibbler.quib.assignment import PathComponent, QuibWithAssignment
+from pyquibbler.quib.assignment import QuibWithAssignment
+from pyquibbler.quib.quib import Quib
+
+if TYPE_CHECKING:
+    from pyquibbler.quib import Assignment, PathComponent
 
 
 class ProxyQuib(Quib):
@@ -16,6 +19,7 @@ class ProxyQuib(Quib):
     def __init__(self, quib):
         super().__init__()
         self._quib: Quib = quib
+        quib.add_child(self)
 
     def _get_inner_value_valid_at_path(self, path: Optional[List[PathComponent]]) -> Any:
         return self._quib.get_value_valid_at_path(path=path)
diff --git a/pyquibbler/quib/quib.py b/pyquibbler/quib/quib.py
index f93023a69..44b40543e 100644
--- a/pyquibbler/quib/quib.py
+++ b/pyquibbler/quib/quib.py
@@ -1,13 +1,12 @@
 from __future__ import annotations
-
 import contextlib
-
 import numpy as np
+from functools import wraps
 from functools import cached_property
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
 from operator import getitem
-from typing import Set, Any, TYPE_CHECKING, Optional, Tuple, Type, List
+from typing import Set, Any, TYPE_CHECKING, Optional, Tuple, Type, List, Callable, Dict
 from weakref import ref as weakref
 
 from pyquibbler.exceptions import PyQuibblerException
@@ -20,7 +19,6 @@
 from .utils import quib_method, Unpacker, recursively_run_func_on_object
 from .assignment import PathComponent
 
-
 from ..env import LEN_RAISE_EXCEPTION
 from ..logger import logger
 
@@ -46,6 +44,31 @@ def __str__(self):
         return f'Cannot override {self.quib} with {self.override} as it does not allow overriding.'
 
 
+@dataclass(frozen=True)
+class FunctionCall:
+    func: Callable
+    args: Tuple[Any, ...]
+    kwargs: Tuple[Tuple[str, Any], ...]
+
+    @classmethod
+    def create(cls, func: Callable, args: Tuple[Any, ...], kwargs: Dict[str, Any]):
+        return cls(func, args, tuple(kwargs.items()))
+
+
+def cache_method_until_full_invalidation(func: Callable) -> Callable:
+    @wraps(func)
+    def wrapper(self: Quib, *args, **kwargs):
+        call = FunctionCall.create(func, (self, *args), kwargs)
+        try:
+            return self.method_cache[call]
+        except KeyError:
+            result = func(self, *args, **kwargs)
+            self.method_cache[call] = result
+            return result
+
+    return wrapper
+
+
 class Quib(ABC):
     """
     An abstract class to describe the common methods and attributes of all quib types.
@@ -62,6 +85,7 @@ def __init__(self,
         if allow_overriding is None:
             allow_overriding = self._DEFAULT_ALLOW_OVERRIDING
         self.allow_overriding = allow_overriding
+        self.method_cache = {}
 
     @property
     def children(self) -> Set[Quib]:
@@ -97,21 +121,18 @@ def __redraw(self) -> None:
         """
         Redraw all artists that directly or indirectly depend on this quib.
         """
-        from pyquibbler.quib.graphics import redraw_axes
+        from pyquibbler.quib.graphics import redraw_axeses
         quibs = self._get_graphics_function_quibs_recursively()
         quibs_that_are_invalid = [quib for quib in quibs if quib.cache_status != CacheStatus.ALL_VALID]
         logger.info(f"redrawing {len(quibs_that_are_invalid)} quibs")
         for graphics_function_quib in quibs_that_are_invalid:
             graphics_function_quib.get_value()
 
-        axeses = set()
-        for graphics_function_quib in quibs_that_are_invalid:
-            for axes in graphics_function_quib.get_axeses():
-                axeses.add(axes)
+        axeses = {axes for graphics_function_quib in quibs_that_are_invalid
+                  for axes in graphics_function_quib.get_axeses()}
 
         logger.info(f"redrawing {len(axeses)} axeses")
-        for axes in axeses:
-            redraw_axes(axes)
+        redraw_axeses(axeses)
 
     def invalidate_and_redraw_at_path(self, path: Optional[List[PathComponent]] = None) -> None:
         """
@@ -143,6 +164,9 @@ def _get_paths_for_children_invalidation(self, invalidator_quib: Quib,
         """
         return [[]]
 
+    def _on_type_change(self):
+        self.method_cache.clear()
+
     def _invalidate_self(self, path: List[PathComponent]):
         """
         This method is called whenever a quib itself is invalidated; subclasses will override this with their
@@ -203,7 +227,7 @@ def _invalidate_quib_with_children_at_path(self, invalidator_quib, path: List[Pa
         for new_path in new_paths:
             if new_path is not None:
                 self._invalidate_self(new_path)
-                if not self._is_completely_overridden_at_first_component(new_path):
+                if len(path) == 0 or not self._is_completely_overridden_at_first_component(new_path):
                     self._invalidate_children_at_path(new_path)
 
     def add_child(self, quib: Quib) -> None:
@@ -236,6 +260,8 @@ def override(self, assignment: Assignment, allow_overriding_from_now_on=True):
         if not self.allow_overriding:
             raise OverridingNotAllowedException(self, assignment)
         self._overrider.add_assignment(assignment)
+        if len(assignment.path) == 0:
+            self._on_type_change()
 
         self.invalidate_and_redraw_at_path(assignment.path)
 
@@ -244,6 +270,8 @@ def remove_override(self, path: List[PathComponent], invalidate_and_redraw: bool
         Remove overriding in a specific path in the quib.
         """
         self._overrider.remove_assignment(path)
+        if len(path) == 0:
+            self._on_type_change()
         if invalidate_and_redraw:
             self.invalidate_and_redraw_at_path(path=path)
 
@@ -331,7 +359,8 @@ def get_value_valid_at_path(self, path: Optional[List[PathComponent]]) -> Any:
             context = contextlib.nullcontext()
 
         with context:
-            return self._overrider.override(self._get_inner_value_valid_at_path(path), self._assignment_template)
+            inner_value = self._get_inner_value_valid_at_path(path)
+            return self._overrider.override(inner_value, self._assignment_template)
 
     def get_value(self) -> Any:
         """
@@ -348,13 +377,16 @@ def get_override_list(self) -> Overrider:
         """
         return self._overrider
 
+    # We cache the type, so quibs without cache will still remember their types.
+    @cache_method_until_full_invalidation
     def get_type(self) -> Type:
         """
-        Get the type of wrapped value
+        Get the type of wrapped value.
         """
         return type(self.get_value_valid_at_path(None))
 
-    @quib_method
+    # We cache the shape, so quibs without cache will still remember their shape.
+    @cache_method_until_full_invalidation
     def get_shape(self) -> Tuple[int, ...]:
         """
         Assuming this quib represents a numpy ndarray, returns a quib of its shape.
@@ -368,26 +400,18 @@ def get_shape(self) -> Tuple[int, ...]:
             raise
 
     @quib_method
-    def _get_override_mask(self, shape: Tuple[int, ...]) -> np.ndarray:
+    def get_override_mask(self):
         """
-        Return an override mask based on assignments in self._overrider, in the given shape.
-        This is an internal method so when called with a quib instead of the shape, the resulting quib
-        will be dependent on both self and the shape quib.
+        Assuming this quib represents a numpy ndarray, return a quib representing its override mask.
+        The override mask is a boolean array of the same shape, in which every value is
+        set to True if the matching value in the array is overridden, and False otherwise.
         """
         if issubclass(self.get_type(), np.ndarray):
-            mask = np.zeros(shape, dtype=np.bool)
+            mask = np.zeros(self.get_shape(), dtype=np.bool)
         else:
             mask = recursively_run_func_on_object(func=lambda x: False, obj=self.get_value())
         return self._overrider.fill_override_mask(mask)
 
-    def get_override_mask(self) -> Quib:
-        """
-        Assuming this quib represents a numpy ndarray, return a quib representing its override mask.
-        The override mask is a boolean array of the same shape, in which every value is
-        set to True if the matching value in the array is overridden, and False otherwise.
-        """
-        return self._get_override_mask(self.get_shape())
-
     def iter_first(self, amount: Optional[int] = None):
         """
         Returns an iterator to the first `amount` elements of the quib.
diff --git a/pyquibbler/third_party_overriding/overriding.py b/pyquibbler/third_party_overriding/overriding.py
index a11528c4c..46ece83e7 100644
--- a/pyquibbler/third_party_overriding/overriding.py
+++ b/pyquibbler/third_party_overriding/overriding.py
@@ -14,7 +14,8 @@
 from pyquibbler.quib.function_quibs.elementwise_function_quib import ElementWiseFunctionQuib
 from pyquibbler.quib.function_quibs.transpositional.transpositional_function_quib import TranspositionalFunctionQuib
 from pyquibbler.quib.graphics import global_collecting, ReductionAxisWiseGraphicsFunctionQuib, \
-    AlongAxisGraphicsFunctionQuib, VectorizeGraphicsFunctionQuib
+    ApplyAlongAxisGraphicsFunctionQuib
+from pyquibbler.quib.graphics import QVectorize
 from pyquibbler.quib.graphics.plot_graphics_function_quib import PlotGraphicsFunctionQuib
 from pyquibbler.quib.graphics.widgets import SliderGraphicsFunctionQuib, CheckButtonsGraphicsFunctionQuib, \
     RadioButtonsGraphicsFunctionQuib, RectangleSelectorGraphicsFunctionQuib, QRectangleSelector, QRadioButtons
@@ -27,35 +28,48 @@
 __eq__
 '''
 
-WIDGET_OVERRIDES = {
-    'RectangleSelector': QRectangleSelector,
-    'RadioButtons': QRadioButtons
-}
+
+def wrap_overridden_graphics_function(func: Callable) -> Callable:
+    @wraps(func)
+    def _wrapper(*args, **kwargs):
+        # We need overridden funcs to be run in `overridden_graphics_function` context manager
+        # so artists will be collected
+        with global_collecting.overridden_graphics_function():
+            return func(*args, **kwargs)
+
+    return _wrapper
+
+
+NON_QUIB_OVERRIDES = [
+    (widgets, {
+        'RectangleSelector': wrap_overridden_graphics_function(QRectangleSelector),
+        'RadioButtons': wrap_overridden_graphics_function(QRadioButtons)
+    }),
+    (np, {
+        'vectorize': QVectorize
+    })
+]
 
 NUMPY_OVERRIDES = [
     (np, [
         (DefaultFunctionQuib, {"abs", "arange", "polyfit",
-                               "linspace", "polyval", "array", "genfromtxt", 'prod', 'corrcoef', 'mean'}),
+                               "linspace", "polyval", "genfromtxt", 'prod', 'corrcoef', 'mean'}),
         (GraphicsFunctionQuib, {'apply_over_axes'}),
         (TranspositionalFunctionQuib, {'reshape', 'rot90', 'ravel', 'concatenate', 'repeat', 'full', 'concatenate',
-                                       'squeeze'}),
+                                       'squeeze', 'array', 'swapaxes', 'expand_dims'}),
         (ElementWiseFunctionQuib, {'add', 'square', "sin", "cos", "tan", "sinh", "cosh", "tanh", "real", "imag",
                                    "arcsin", "arccos", "arctan", "arcsinh", "arccosh", "arctanh",
                                    "exp", "exp2", "expm1", "log", "log2", "log1p", "log10",
                                    "sqrt", "int", "float", "ceil", "floor", "round", 'around', 'hypot'}),
         (ReductionAxisWiseGraphicsFunctionQuib, {"max", "amax", "min", "amin", "sum", "average"}),
-        (AlongAxisGraphicsFunctionQuib, {'apply_along_axis'}),
+        (ApplyAlongAxisGraphicsFunctionQuib, {'apply_along_axis'}),
     ]),
     (np.random, [
         (ImpureFunctionQuib, {'rand', 'randn', 'randint'})
     ]),
-    (np.vectorize, [
-        (VectorizeGraphicsFunctionQuib, {'__call__'})
-    ]),
     (np.linalg, [
         (ReductionAxisWiseGraphicsFunctionQuib, {'norm'})
-    ]
-    )
+    ])
 ]
 
 MPL_OVERRIDES = [
@@ -79,17 +93,6 @@
 ]
 
 
-def wrap_overridden_graphics_function(func: Callable) -> Callable:
-    @wraps(func)
-    def _wrapper(*args, **kwargs):
-        # We need overridden funcs to be run in `overridden_graphics_function` context manager
-        # so artists will be collected
-        with global_collecting.overridden_graphics_function():
-            return func(*args, **kwargs)
-
-    return _wrapper
-
-
 def override_func(obj: Any, name: str, quib_type: Type[FunctionQuib],
                   function_wrapper: Optional[Callable[[Callable], Callable]] = None):
     """
@@ -100,7 +103,6 @@ def override_func(obj: Any, name: str, quib_type: Type[FunctionQuib],
     if function_wrapper is not None:
         func = function_wrapper(func)
     func = quib_type.create_wrapper(func)
-
     setattr(obj, name, func)
 
 
@@ -112,9 +114,12 @@ def apply_quib_creating_overrides(override_list: List[Tuple[Any, List[Tuple[Type
                 override_func(obj, func_name, quib_type, function_wrapper)
 
 
-def override_widgets():
-    for name, new_cls in WIDGET_OVERRIDES.items():
-        setattr(widgets, name, wrap_overridden_graphics_function(new_cls))
+def apply_non_quib_overrides():
+    for module, overrides in NON_QUIB_OVERRIDES:
+        for name, replacement in overrides.items():
+            overridden = getattr(module, name)
+            setattr(module, name, replacement)
+            replacement.__overridden__ = overridden
 
 
 @ensure_only_run_once_globally
@@ -122,6 +127,6 @@ def override_all():
     """
     Overrides all modules (such as numpy and matplotlib) to support quibs
     """
-    override_widgets()
+    apply_non_quib_overrides()
     apply_quib_creating_overrides(NUMPY_OVERRIDES)
     apply_quib_creating_overrides(MPL_OVERRIDES, wrap_overridden_graphics_function)


Test patch:
diff --git a/tests/functional/decorators/test_quibbler_user_function.py b/tests/functional/decorators/test_quibbler_user_function.py
index 4a3462b70..df700bc5b 100644
--- a/tests/functional/decorators/test_quibbler_user_function.py
+++ b/tests/functional/decorators/test_quibbler_user_function.py
@@ -27,7 +27,7 @@ def test_quibbler_user_function_non_lazy():
 
 def test_quibbler_user_function_with_quibs():
     mock_func = mock.Mock()
-    user_function = quibbler_user_function(receive_quibs=True, lazy=False)(mock_func)
+    user_function = quibbler_user_function(pass_quibs=True, lazy=False)(mock_func)
     quib = iquib(6)
 
     res = user_function(quib)
diff --git a/tests/functional/quib/assignment/inverse_functions/test_transpositional_inversal.py b/tests/functional/quib/assignment/inverse_functions/test_transpositional_inversal.py
index 4a235f473..434312614 100644
--- a/tests/functional/quib/assignment/inverse_functions/test_transpositional_inversal.py
+++ b/tests/functional/quib/assignment/inverse_functions/test_transpositional_inversal.py
@@ -288,3 +288,12 @@ def test_inverse_with_resulting_int_and_changing_value_shape():
     b[:, :] = 0
 
     assert np.array_equal(a.get_value(), np.full((2, 3), fill_value=0))
+
+
+def test_inverse_np_array():
+    a = iquib([[1, 2, 3, 4]])
+    b = np.array(a)
+
+    b[(0, 0)] = 20
+
+    assert a[0][0].get_value() == 20
diff --git a/tests/functional/quib/function_quibs/conftest.py b/tests/functional/quib/function_quibs/conftest.py
index 2de8c0cd6..94d7f9260 100644
--- a/tests/functional/quib/function_quibs/conftest.py
+++ b/tests/functional/quib/function_quibs/conftest.py
@@ -12,7 +12,7 @@ def _create(shape=None, get_value_result=None, children=None):
         get_value_result = get_value_result or [[1, 2, 3]]
         mock_quib = mock.Mock(spec=Quib)
         mock_quib.get_value_valid_at_path.return_value = get_value_result
-        mock_quib.get_shape.return_value.get_value.return_value = shape
+        mock_quib.get_shape.return_value = shape
         mock_quib._get_children_recursively.return_value = children or set()
         return mock_quib
     return _create
diff --git a/tests/functional/quib/function_quibs/test_elementwise_function_quib.py b/tests/functional/quib/function_quibs/test_elementwise_function_quib.py
index 9e40997d8..6d767e00b 100644
--- a/tests/functional/quib/function_quibs/test_elementwise_function_quib.py
+++ b/tests/functional/quib/function_quibs/test_elementwise_function_quib.py
@@ -43,7 +43,7 @@ def test_elementwise_function_quib_invalidation_with_broadcast_numpy_array():
 def test_elementwise_function_quib_does_not_request_unneeded_indices_on_get_value():
     fake_quib = mock.Mock(spec=Quib)
     fake_quib.get_value_valid_at_path.return_value = np.array([1, 2])
-    fake_quib.get_shape.return_value.get_value.return_value = (2,)
+    fake_quib.get_shape.return_value = (2,)
     b = ElementWiseFunctionQuib.create(
         func=np.add,
         func_args=(fake_quib, 1)
diff --git a/tests/functional/quib/function_quibs/test_transpositional_quib.py b/tests/functional/quib/function_quibs/test_transpositional_quib.py
index 3d1031642..fa2b8a9f5 100644
--- a/tests/functional/quib/function_quibs/test_transpositional_quib.py
+++ b/tests/functional/quib/function_quibs/test_transpositional_quib.py
@@ -256,6 +256,13 @@ def filter_out_none_calls(mock_calls):
     ]
 
 
+def create_mock_quib(shape, get_value_result):
+    mock_quib = mock.Mock(spec=Quib)
+    mock_quib.get_value_valid_at_path.return_value = get_value_result
+    mock_quib.get_shape.return_value = shape
+    return mock_quib
+
+
 @pytest.mark.regression
 def test_transpositional_concatenate_does_diverge(create_mock_quib):
     first = create_mock_quib((1, 3), ([[1, 2, 3]]))
diff --git a/tests/functional/quib/graphics/conftest.py b/tests/functional/quib/graphics/conftest.py
new file mode 100644
index 000000000..d1aa71b64
--- /dev/null
+++ b/tests/functional/quib/graphics/conftest.py
@@ -0,0 +1,50 @@
+import pytest
+from unittest import mock
+
+from pyquibbler.quib.graphics import graphics_function_quib, utils
+from pyquibbler.quib.graphics.axiswise_function_quibs import apply_along_axis_function_quib
+from pyquibbler.quib.graphics.axiswise_function_quibs.vectorize import vectorize_graphics_function_quib
+
+
+class MockArtistsCollector:
+    def __init__(self, mock_axes):
+        self.mock_axes = mock_axes
+        self.all_mock_artists_created = []
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, exc_type, exc_val, exc_tb):
+        pass
+
+    def _create_mock_artist(self):
+        artist = mock.Mock()
+        artist.axes = self.mock_axes
+        self.mock_axes.artists.append(artist)
+        return artist
+
+    @property
+    def objects_collected(self):
+        artist = self._create_mock_artist()
+        self.all_mock_artists_created.append(artist)
+        return [artist]
+
+    @property
+    def mock_artists_in_axes(self):
+        return [artist for artist in self.all_mock_artists_created if artist in artist.axes.artists]
+
+
+@pytest.fixture()
+def mock_axes():
+    axes = mock.Mock()
+    axes.figure.canvas.supports_blit = False
+    axes.artists = []
+    return axes
+
+
+@pytest.fixture()
+def mock_artists_collector(mock_axes, monkeypatch):
+    collector = MockArtistsCollector(mock_axes)
+    monkeypatch.setattr(graphics_function_quib, "ArtistsCollector", lambda: collector)
+    monkeypatch.setattr(utils, "ArtistsCollector", lambda: collector)
+    return collector
diff --git a/tests/functional/quib/graphics/test_apply_along_axis_function_quib.py b/tests/functional/quib/graphics/test_apply_along_axis_function_quib.py
new file mode 100644
index 000000000..9dc7dfb15
--- /dev/null
+++ b/tests/functional/quib/graphics/test_apply_along_axis_function_quib.py
@@ -0,0 +1,275 @@
+import dataclasses
+import itertools
+from typing import Callable, Any
+
+import numpy as np
+import pytest
+from pytest import mark
+
+from pyquibbler import iquib
+from pyquibbler.env import GRAPHICS_LAZY
+from pyquibbler.quib import PathComponent, Quib
+from pyquibbler.quib.assignment.utils import get_sub_data_from_object_in_path
+from pyquibbler.quib.graphics import ApplyAlongAxisGraphicsFunctionQuib
+from tests.functional.quib.graphics.test_axiswise_graphics_function_quib import parametrize_indices_to_invalidate, \
+    parametrize_data
+from tests.functional.quib.utils import check_invalidation, check_get_value_valid_at_path, get_func_mock
+
+
+@parametrize_indices_to_invalidate
+@parametrize_data
+@mark.parametrize('axis', [0, 1, 2, -1, -2])
+@mark.parametrize('func_out_dims', [0, 1, 2])
+def test_apply_along_axis_invalidation(indices_to_invalidate, axis, func_out_dims, data):
+    func1d = lambda slice: np.sum(slice).reshape((1,) * func_out_dims)
+    check_invalidation(lambda quib: np.apply_along_axis(func1d, axis, quib), data, indices_to_invalidate)
+
+
+@parametrize_data
+@mark.parametrize('axis', [0, 1, 2, -1, -2])
+@mark.parametrize('func_out_dims', [0, 1, 2])
+@mark.parametrize('indices_to_get_value_at', [0, (0, 0), (-1, ...)])
+def test_apply_along_axis_get_value_valid_at_path(indices_to_get_value_at, axis, func_out_dims, data):
+    func1d = lambda slice: np.sum(slice).reshape((1,) * func_out_dims)
+    path_to_get_value_at = [PathComponent(np.ndarray, indices_to_get_value_at)]
+
+    check_get_value_valid_at_path(lambda quib: np.apply_along_axis(func1d, axis, quib), data, path_to_get_value_at)
+
+
+def create_lazy_apply_along_axis_quib(func, arr, axis, args=None, kwargs=None, pass_quibs=False):
+    with GRAPHICS_LAZY.temporary_set(True):
+        return np.apply_along_axis(func, axis, iquib(arr) if not isinstance(arr, Quib) else arr,
+                                   *(args or []), **(kwargs or {}), pass_quibs=pass_quibs)
+
+
+@mark.parametrize('shape, axis, func1d_res', [
+        (tuple(dimensions), axis, res)
+        for shape_size in range(0, 3)
+        for axis in range(-shape_size, shape_size)
+        for dimensions in itertools.product(*[range(1, 4) for _ in range(shape_size)])
+        for res in [1, np.arange(9).reshape((3, 3)), None]
+])
+@mark.parametrize('pass_quibs', [True, False])
+def test_apply_along_axis_get_shape(shape, axis, func1d_res, pass_quibs):
+    func = get_func_mock(lambda x: func1d_res)
+    arr = np.arange(np.prod(shape)).reshape(shape)
+    quib = create_lazy_apply_along_axis_quib(func=func, arr=arr, axis=axis, pass_quibs=pass_quibs)
+    expected_input_arr = arr[tuple([slice(None) if i == quib.core_axis else 0 for i in range(len(arr.shape))])]
+
+    res = quib.get_shape()
+
+    assert func.call_count == 1
+    call = func.mock_calls[0]
+    oned_slice = call.args[0]
+    if pass_quibs:
+        assert isinstance(oned_slice, Quib)
+        oned_slice = oned_slice.get_value()
+    else:
+        assert not isinstance(oned_slice, Quib)
+    assert np.array_equal(oned_slice, expected_input_arr)
+    assert res == np.apply_along_axis(func, axis=axis, arr=arr).shape
+
+
+def get_sample_indexing_paths(input_shape, apply_shape, axis):
+    res = np.apply_along_axis(func1d=lambda _: np.zeros(apply_shape), arr=np.zeros(input_shape), axis=axis)
+
+    paths = [[()], []]
+
+    if len(res.shape) > 0:
+        paths.append([np.array([0 for _ in res.shape])])
+
+    if len(res.shape) > 1:
+        paths.append([(0, np.arange(res.shape[1]))])
+        paths.append([(0, res.shape[1] - 1)])
+    return paths
+
+
+def assert_all_apply_calls_with_slices_were_relevant(func, axis, input_arr, path, applied_slices):
+    whole_result = np.apply_along_axis(func, axis=axis, arr=input_arr)
+    current_result = get_sub_data_from_object_in_path(whole_result, path)
+
+    for slice_ in applied_slices:
+        for num in np.ravel(slice_):
+            new_arr = np.array(input_arr)
+            new_arr[new_arr == num] = 999
+            new_result = np.apply_along_axis(func, axis=axis, arr=new_arr)
+            new_result_at_path = get_sub_data_from_object_in_path(new_result, path)
+            assert not np.array_equal(new_result_at_path, current_result )
+
+
+@mark.parametrize('input_shape, apply_result_shape, axis, components', [
+        (tuple(input_dimensions), tuple(apply_dimensions), axis, components)
+        for input_shape_size in range(0, 3)
+        for axis in range(-input_shape_size, input_shape_size)
+        for apply_result_shape_size in range(0, 3)
+        for input_dimensions in itertools.product(*[range(1, 3) for _ in range(input_shape_size)])
+        for apply_dimensions in [[], *itertools.product(*[range(1, 3) for _ in range(apply_result_shape_size)])]
+        for components in get_sample_indexing_paths(tuple(input_dimensions), tuple(apply_dimensions), axis)
+])
+@mark.parametrize('pass_quibs', [True, False])
+def test_apply_along_axis_get_value(input_shape, apply_result_shape, axis, components, pass_quibs):
+    assert len(components) <= 1, "Sanity: No support for testing multiple components (also irrelevant)"
+    slices = []
+    running_in_quib = False
+
+    def apply(oned_slice):
+        if running_in_quib:
+            if pass_quibs:
+                assert isinstance(oned_slice, Quib)
+                oned_slice = oned_slice.get_value()
+            else:
+                assert not isinstance(oned_slice, Quib)
+
+            slices.append(oned_slice)
+        return np.full(apply_result_shape, np.sum(oned_slice))
+
+    func = get_func_mock(apply)
+    arr = np.arange(np.prod(input_shape)).reshape(input_shape)
+    quib = create_lazy_apply_along_axis_quib(arr=arr, func=func, axis=axis, pass_quibs=pass_quibs)
+    quib.get_shape()  # We need to call get_shape to cache it as get_shape is a zero cost operation in overall scheme
+    # and is allowed to be called without consequence by the quib
+    path = [PathComponent(component=component, indexed_cls=np.ndarray) for component in components]
+    running_in_quib = True
+
+    res = quib.get_value_valid_at_path(path)
+
+    running_in_quib = False
+    whole_apply_axis_result = np.apply_along_axis(func, axis=axis, arr=arr)
+    expected_result = get_sub_data_from_object_in_path(whole_apply_axis_result, path)
+    res_at_components = get_sub_data_from_object_in_path(res, path)
+    assert np.array_equal(res_at_components, expected_result)
+    assert_all_apply_calls_with_slices_were_relevant(
+        func=func,
+        axis=axis,
+        input_arr=arr,
+        applied_slices=slices,
+        path=path
+    )
+
+
+def test_apply_along_axis_returning_quib():
+    res_arr = np.array([100])
+    func = lambda x: iquib(res_arr)
+    arr = np.array([[1]])
+    axis = 0
+    quib = create_lazy_apply_along_axis_quib(func=func, arr=arr, axis=axis)
+    expected_res = np.apply_along_axis(lambda _: np.array([100]), arr=arr, axis=axis)
+
+    res = quib.get_value()
+
+    assert np.array_equal(res, expected_res)
+
+
+@pytest.fixture
+def mock_func_for_args_kwargs():
+    return get_func_mock(lambda x, *_, **__: 1)
+
+
+@pytest.fixture
+def args():
+    return [11, 1, 1]
+
+
+@pytest.fixture
+def kwargs():
+    return {
+        'funtimes': 1
+    }
+
+
+@pytest.fixture
+def quib_with_args_and_kwargs(mock_func_for_args_kwargs, args, kwargs):
+    return create_lazy_apply_along_axis_quib(func=mock_func_for_args_kwargs,
+                                             arr=np.array([[1]]),
+                                             axis=0,
+                                             args=args,
+                                             kwargs=kwargs)
+
+
+def test_apply_along_axis_get_value_with_args_and_kwargs(quib_with_args_and_kwargs, mock_func_for_args_kwargs,
+                                                         args, kwargs):
+    quib_with_args_and_kwargs.get_value()
+
+    mock_call = mock_func_for_args_kwargs.mock_calls[-1]
+    assert mock_call.kwargs == kwargs
+    assert mock_call.args[1:] == tuple(args)
+
+
+def test_apply_along_axis_get_shape_with_args_and_kwargs(quib_with_args_and_kwargs, mock_func_for_args_kwargs,
+                                                         args, kwargs):
+    quib_with_args_and_kwargs.get_shape()
+
+    assert len(mock_func_for_args_kwargs.mock_calls) == 1
+    mock_call = mock_func_for_args_kwargs.mock_calls[0]
+    assert mock_call.kwargs == kwargs
+    assert mock_call.args[1:] == tuple(args)
+
+
+def test_apply_along_axis_get_shape_with_looping_axis_quib():
+    axis_quib = iquib(0)
+    quib = create_lazy_apply_along_axis_quib(
+        arr=np.array([[1]]),
+        axis=axis_quib,
+        func=lambda x: 1
+    )
+
+    shape = quib.get_shape()
+
+    assert shape == (1,)
+
+
+def test_apply_along_axis_get_shape_with_list():
+    arr_quib = iquib([[1, 2], [3, 4]])
+    quib = create_lazy_apply_along_axis_quib(
+        arr=arr_quib,
+        axis=0,
+        func=lambda x: 1,
+    )
+
+    assert quib.get_shape() == (2,)
+
+
+def test_apply_along_axis_get_shape_does_not_create_artists(mock_artists_collector):
+    quib = create_lazy_apply_along_axis_quib(
+        arr=np.array([[1, 2], [3, 4]]),
+        axis=0,
+        func=lambda x: 1
+    )
+
+    quib.get_shape()
+
+    assert len(mock_artists_collector.mock_artists_in_axes) == 0
+
+
+def test_apply_along_axis_removes_and_recreates_artists(mock_artists_collector):
+    parent = iquib([[1, 2], [3, 4]])
+    quib = create_lazy_apply_along_axis_quib(
+        arr=parent,
+        axis=0,
+        func=lambda x: 1
+    )
+
+    quib.get_value()
+    mock_artists = mock_artists_collector.mock_artists_in_axes
+    parent.invalidate_and_redraw_at_path([PathComponent(component=0, indexed_cls=np.ndarray)])
+    quib.get_value()
+
+    assert len(mock_artists_collector.mock_artists_in_axes) == 2
+    assert set(mock_artists_collector.mock_artists_in_axes) != set(mock_artists)
+
+
+def test_apply_along_axis_does_not_remove_artists_that_are_not_his(mock_artists_collector):
+    quib = create_lazy_apply_along_axis_quib(
+        arr=np.array([[1, 2], [3, 4]]),
+        axis=0,
+        func=lambda x: 1
+    )
+    quib.get_value_valid_at_path([PathComponent(component=0, indexed_cls=np.ndarray)])
+    artists = mock_artists_collector.mock_artists_in_axes
+    assert len(artists) == 1, "Sanity check"
+
+    quib.get_value_valid_at_path([PathComponent(component=1, indexed_cls=np.ndarray)])
+
+    assert len(mock_artists_collector.mock_artists_in_axes) == 2
+    assert len(set(artists) & set(mock_artists_collector.mock_artists_in_axes)) == 1
+
diff --git a/tests/functional/quib/graphics/test_axiswise_graphics_function_quib.py b/tests/functional/quib/graphics/test_axiswise_graphics_function_quib.py
index 06cfa4e16..cf7cc16c6 100644
--- a/tests/functional/quib/graphics/test_axiswise_graphics_function_quib.py
+++ b/tests/functional/quib/graphics/test_axiswise_graphics_function_quib.py
@@ -1,10 +1,15 @@
+import itertools
+from unittest import mock
+
 import numpy as np
+import pytest
 from pytest import mark
 
-from pyquibbler import iquib, CacheBehavior
+from pyquibbler import CacheBehavior, iquib
+from pyquibbler.env import GRAPHICS_LAZY
 from pyquibbler.quib.assignment import PathComponent
 
-from ..utils import check_invalidation, check_get_value_valid_at_path, MockQuib, PathBuilder
+from ..utils import check_invalidation, check_get_value_valid_at_path, MockQuib, get_func_mock
 
 # A 3d array in which every dimension has a different size
 parametrize_data = mark.parametrize('data', [np.arange(24).reshape((2, 3, 4))])
@@ -70,111 +75,3 @@ def test_reduction_axiswise_get_value_valid_at_path(axis, data, keepdims, where,
         kwargs['where'] = where
     path_to_get_value_at = [PathComponent(np.ndarray, indices_to_get_value_at)]
     check_get_value_valid_at_path(lambda quib: np.sum(quib, **kwargs), data, path_to_get_value_at)
-
-
-@parametrize_indices_to_invalidate
-@parametrize_data
-@mark.parametrize('axis', [0, 1, 2, -1, -2])
-@mark.parametrize('func_out_dims', [0, 1, 2])
-def test_apply_along_axis_invalidation(indices_to_invalidate, axis, func_out_dims, data):
-    func1d = lambda slice: np.sum(slice).reshape((1,) * func_out_dims)
-    check_invalidation(lambda quib: np.apply_along_axis(func1d, axis, quib), data, indices_to_invalidate)
-
-
-@parametrize_data
-@mark.parametrize('axis', [0, 1, 2, -1, -2])
-@mark.parametrize('func_out_dims', [0, 1, 2])
-@mark.parametrize('indices_to_get_value_at', [0, (0, 0), (-1, ...)])
-def test_apply_along_axis_get_value_valid_at_path(indices_to_get_value_at, axis, func_out_dims, data):
-    func1d = lambda slice: np.sum(slice).reshape((1,) * func_out_dims)
-    path_to_get_value_at = [PathComponent(np.ndarray, indices_to_get_value_at)]
-    check_get_value_valid_at_path(lambda quib: np.apply_along_axis(func1d, axis, quib), data, path_to_get_value_at)
-
-
-@parametrize_indices_to_invalidate
-@parametrize_data
-@mark.parametrize('excluded', [{0}, set(), None])
-@mark.parametrize('func', [lambda x: np.sum(x), lambda x: (np.sum(x), np.sum(x))])
-def test_vectorize_invalidation(indices_to_invalidate, data, excluded, func):
-    kwargs = {}
-    if excluded is not None:
-        kwargs['excluded'] = excluded
-    check_invalidation(np.vectorize(func, **kwargs), data, indices_to_invalidate)
-
-
-@parametrize_data
-@mark.parametrize('func', [lambda x: np.sum(x)])
-@mark.parametrize('indices_to_get_value_at', [0, (0, 0), (-1, ...)])
-def test_vectorize_get_value_valid_at_path(data, func, indices_to_get_value_at):
-    path_to_get_value_at = [PathComponent(np.ndarray, indices_to_get_value_at)]
-    check_get_value_valid_at_path(np.vectorize(func), data, path_to_get_value_at)
-
-
-def test_vectorize_get_value_valid_at_path_with_excluded_quib():
-    excluded = MockQuib([1, 2, 3])
-    func = np.vectorize(lambda a, b: np.array([1, 2, 3]), excluded={1}, signature='(n)->(m)')
-    fquib = func([0, 1], excluded)
-    fquib.set_cache_behavior(CacheBehavior.OFF)
-    path = PathBuilder(fquib)[0].path
-
-    with excluded.collect_valid_paths() as valid_paths:
-        fquib.get_value_valid_at_path(path)
-
-    assert valid_paths == [[]]
-
-
-@mark.parametrize(['quib_data', 'non_quib_data'], [
-    (np.zeros((2, 3)), np.zeros((4, 2, 3))),
-    (np.zeros((3, 3)), np.zeros((1, 3))),
-    (np.zeros((1, 3, 3)), np.zeros((4, 1, 3))),
-    (np.zeros((4, 2, 3)), np.zeros((2, 3))),
-])
-def test_vectorize_get_value_valid_at_path_when_args_have_different_loop_dimensions(quib_data, non_quib_data):
-    func = lambda quib: np.vectorize(lambda x, y: x + y)(quib, quib_data)
-    check_get_value_valid_at_path(func, non_quib_data, [PathComponent(np.ndarray, 0)])
-
-
-@mark.parametrize('indices_to_get_value_at', [0, (1, 1), (1, ..., 2)])
-def test_vectorize_get_value_at_path_with_core_dims(indices_to_get_value_at):
-    quib_data = np.zeros((2, 2, 3, 4))
-    non_quib_data = np.zeros((2, 3, 4, 5))
-    func = lambda a, b: np.array([np.sum(a) + np.sum(b)] * 6)
-    vec = np.vectorize(func, signature='(a,b),(c)->(d)')
-    check_get_value_valid_at_path(lambda quib: vec(non_quib_data, quib), quib_data,
-                                  [PathComponent(np.ndarray, indices_to_get_value_at)])
-
-
-def test_vectorize_invalidation_with_non_numpy_return_value():
-    vec = np.vectorize(lambda a: int(np.sum(a)), signature='(a)->()')
-    check_invalidation(lambda quib: vec(quib), [1, 2, 3], 0)
-
-
-@parametrize_indices_to_invalidate
-@parametrize_data
-def test_vectorize_invalidation_with_multiple_params(data, indices_to_invalidate):
-    quib_arg = iquib(np.arange(12).reshape((3, 4)))
-    non_quib_arg = np.arange(8).reshape((2, 1, 4))
-    quib_kwarg = iquib(3)
-    non_quib_kwarg = np.arange(3).reshape((3, 1))
-
-    vec = np.vectorize(lambda a, b, c, d, e=0: a + b + c + d + e)
-    check_invalidation(lambda quib: vec(quib, quib_arg, non_quib_arg, d=quib_kwarg, e=non_quib_kwarg),
-                       data, indices_to_invalidate)
-
-
-@mark.parametrize('excluded', [{1}, {'b'}, set(), None])
-def test_vectorize_invalidation_with_quib_as_kwarg(excluded):
-    kwargs = {}
-    if excluded is not None:
-        kwargs['excluded'] = excluded
-    vec = np.vectorize(lambda a, b: a + b)
-    check_invalidation(lambda quib: vec(5, b=quib), [1, 2, 3], 0)
-
-
-@parametrize_indices_to_invalidate
-@parametrize_data
-def test_vectorize_invalidation_with_different_core_dims(data, indices_to_invalidate):
-    data2 = np.arange(100, 122).reshape(1, 22)
-    func = lambda a, b: (np.array([np.sum(a) + np.sum(b)] * 4), np.sum(a) + np.sum(b))
-    vec = np.vectorize(func, signature='(a,b),(c)->(d),()')
-    check_invalidation(lambda quib: vec(quib, data2), data, indices_to_invalidate)
diff --git a/tests/functional/quib/graphics/test_global_collecting.py b/tests/functional/quib/graphics/test_global_collecting.py
index 1db702855..7cde9929f 100644
--- a/tests/functional/quib/graphics/test_global_collecting.py
+++ b/tests/functional/quib/graphics/test_global_collecting.py
@@ -6,7 +6,6 @@
 
 from pyquibbler import iquib
 from pyquibbler.quib.graphics import global_collecting
-from pyquibbler.quib.graphics.global_collecting import AlreadyCollectingArtistsException
 
 
 @pytest.mark.parametrize("data", [
@@ -17,42 +16,16 @@ def test_global_graphics_collecting_mode_happy_flow(axes, data):
     with global_collecting.ArtistsCollector() as collector:
         plt.plot(data)
 
-    assert len(collector.artists_collected) == 1
-    assert isinstance(collector.artists_collected[0], Line2D)
-
-
-def test_global_graphics_collecting_raises_exception_when_within_collector():
-    with pytest.raises(AlreadyCollectingArtistsException):
-        with global_collecting.ArtistsCollector(), global_collecting.ArtistsCollector(raise_if_within_collector=True):
-            pass
-
-
-def test_global_graphics_collecting_when_within_collector():
-    with global_collecting.ArtistsCollector():
-        assert global_collecting.is_within_artists_collector()
-
-
-def test_global_graphics_collecting_when_not_within_collector():
-    with global_collecting.ArtistsCollector():
-        pass
-
-    assert global_collecting.is_within_artists_collector() is False
+    assert len(collector.objects_collected) == 1
+    assert isinstance(collector.objects_collected[0], Line2D)
 
 
 def test_global_graphics_collecting_within_collecting(axes):
     with global_collecting.ArtistsCollector() as collector, \
             global_collecting.ArtistsCollector() as nested_collector:
         plt.plot([1, 2, 3])
-    assert len(collector.artists_collected) == 1
-    assert len(nested_collector.artists_collected) == 1
-
-
-@pytest.mark.regression
-def test_global_graphics_collecting_within_collecting_is_within_collector(axes):
-    with global_collecting.ArtistsCollector():
-        with global_collecting.ArtistsCollector():
-            pass
-        assert global_collecting.is_within_artists_collector()
+    assert len(collector.objects_collected) == 1
+    assert len(nested_collector.objects_collected) == 1
 
 
 @pytest.mark.regression
@@ -62,4 +35,4 @@ def test_overridden_graphics_function_within_overridden_graphics_function(axes):
             pass
         artist = Artist()
 
-    assert collector.artists_collected == [artist]
+    assert collector.objects_collected == [artist]
diff --git a/tests/functional/quib/graphics/test_graphics_function_quib.py b/tests/functional/quib/graphics/test_graphics_function_quib.py
index 83d2f8bfe..d3ef0cb92 100644
--- a/tests/functional/quib/graphics/test_graphics_function_quib.py
+++ b/tests/functional/quib/graphics/test_graphics_function_quib.py
@@ -3,39 +3,9 @@
 from unittest import mock
 
 from pyquibbler import iquib
+from pyquibbler.env import GRAPHICS_LAZY
 from pyquibbler.quib.assignment.assignment import PathComponent
-from pyquibbler.quib.graphics import GraphicsFunctionQuib, global_collecting
-
-
-class MockArtistsCollector:
-    def __init__(self, mock_axes):
-        self.mock_axes = mock_axes
-        self.all_mock_artists_created = []
-
-    def __enter__(self):
-        return self
-
-    def __exit__(self, exc_type, exc_val, exc_tb):
-        pass
-
-    def _create_mock_artist(self):
-        artist = mock.Mock()
-        artist.axes = self.mock_axes
-        self.mock_axes.artists.append(artist)
-        return artist
-
-    @property
-    def artists_collected(self):
-        artist = self._create_mock_artist()
-        self.all_mock_artists_created.append(artist)
-        return [artist]
-
-
-@pytest.fixture()
-def mock_axes():
-    axes = mock.Mock()
-    axes.artists = []
-    return axes
+from pyquibbler.quib.graphics import GraphicsFunctionQuib
 
 
 def get_graphics_quib(func):
@@ -43,8 +13,7 @@ def get_graphics_quib(func):
         args=tuple(),
         kwargs={},
         cache_behavior=None,
-        func=func,
-        artists=[]
+        func=func
     )
 
 
@@ -57,16 +26,13 @@ def test_graphics_function_quib_get_value_returns_value():
     assert res == mock_func.return_value
 
 
-def test_graphics_function_quib_rerun_removes_artists_created(monkeypatch, mock_axes):
-    mock_artists_collector = MockArtistsCollector(mock_axes)
-    monkeypatch.setattr(global_collecting, "ArtistsCollector", lambda: mock_artists_collector)
-
+def test_graphics_function_quib_rerun_removes_artists_created(mock_axes, mock_artists_collector):
     father_quib = iquib(1)
     quib = get_graphics_quib(mock.Mock())
     father_quib.add_child(quib)
 
     quib.get_value()
-    father_quib.invalidate_and_redraw_at_path(path=[PathComponent(component=..., indexed_cls=np.ndarray)])
+    father_quib.invalidate_and_redraw_at_path()
 
     assert len(mock_artists_collector.all_mock_artists_created) == 2
     assert len(mock_axes.artists) == 1
@@ -104,3 +70,11 @@ def test_graphics_function_quib_does_not_copy_color(axes):
     assert artist_color_upon_creation == [1, 0, 0]
     assert artist_color_after_color_change == [1, 1, 0]
 
+
+def test_graphics_function_quib_does_not_run_when_lazy_flag_set_to_true():
+    func = mock.Mock()
+
+    with GRAPHICS_LAZY.temporary_set(True):
+        GraphicsFunctionQuib.create(func=func)
+
+    assert func.call_count == 0
diff --git a/tests/functional/quib/graphics/test_graphics_function_quib_graphics.py b/tests/functional/quib/graphics/test_graphics_function_quib_graphics.py
index ee5483506..2b3d68ad7 100644
--- a/tests/functional/quib/graphics/test_graphics_function_quib_graphics.py
+++ b/tests/functional/quib/graphics/test_graphics_function_quib_graphics.py
@@ -4,7 +4,7 @@
 from matplotlib.testing.decorators import image_comparison
 
 from pyquibbler import iquib
-from pyquibbler.quib.graphics import PlotGraphicsFunctionQuib, redraw_axes
+from pyquibbler.quib.graphics import PlotGraphicsFunctionQuib, redraw_axeses
 
 # This import is critical!! Do not delete!!
 # This imports a fixture (it has autouse=True) that cleans up matplotlib stuffs our tests do
@@ -13,8 +13,7 @@
 
 def run_graphics_function_quib_and_redraw(q: PlotGraphicsFunctionQuib):
     q.get_value()
-    for axes in q.get_axeses():
-        redraw_axes(axes)
+    redraw_axeses(q.get_axeses())
 
 
 quibbler_image_comparison = functools.partial(image_comparison, remove_text=True, extensions=['png'],
diff --git a/tests/functional/quib/graphics/test_redraw.py b/tests/functional/quib/graphics/test_redraw.py
index 7277d03ca..6967b9311 100644
--- a/tests/functional/quib/graphics/test_redraw.py
+++ b/tests/functional/quib/graphics/test_redraw.py
@@ -1,30 +1,31 @@
-from unittest import mock
+from unittest.mock import Mock
 
-import pytest
-
-from pyquibbler.quib.graphics import redraw_axes
+from pyquibbler.quib.graphics import redraw_axeses
 from pyquibbler.quib.graphics.redraw import aggregate_redraw_mode
 
 
-@pytest.fixture()
-def mock_axes():
-    axes = mock.Mock()
-    axes.figure.canvas.supports_blit = False
-    return axes
-
-
 def test_redraw_axes_happy_flow(mock_axes):
-
-    redraw_axes(mock_axes)
+    redraw_axeses({mock_axes})
 
     mock_axes.figure.canvas.draw.assert_called_once()
 
 
 def test_redraw_in_aggregate_mode(mock_axes):
-
     with aggregate_redraw_mode():
-        redraw_axes(mock_axes)
-        redraw_axes(mock_axes)
-        redraw_axes(mock_axes)
+        redraw_axeses({mock_axes})
+        redraw_axeses({mock_axes})
+        redraw_axeses({mock_axes})
 
     mock_axes.figure.canvas.draw.assert_called_once()
+
+
+def test_redraw_axeses_does_not_redraw_same_canvas_twice():
+    mock_axes1 = Mock()
+    mock_axes2 = Mock()
+    mock_canvas = Mock()
+    mock_axes1.figure.canvas = mock_canvas
+    mock_axes2.figure.canvas = mock_canvas
+
+    redraw_axeses({mock_axes1, mock_axes2})
+
+    mock_canvas.draw.assert_called_once()
diff --git a/tests/functional/quib/graphics/test_replacing_graphics_function_quib.py b/tests/functional/quib/graphics/test_replacing_graphics_function_quib.py
index f763433ea..f917e89c9 100644
--- a/tests/functional/quib/graphics/test_replacing_graphics_function_quib.py
+++ b/tests/functional/quib/graphics/test_replacing_graphics_function_quib.py
@@ -1,7 +1,5 @@
-from unittest import mock
-
 import pytest
-from matplotlib.artist import Artist
+from unittest import mock
 from matplotlib.axes import Axes
 
 from pyquibbler import iquib
diff --git a/tests/functional/quib/graphics/vectorize/__init__.py b/tests/functional/quib/graphics/vectorize/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/functional/quib/graphics/vectorize/test_utils.py b/tests/functional/quib/graphics/vectorize/test_utils.py
new file mode 100644
index 000000000..f8fdac6af
--- /dev/null
+++ b/tests/functional/quib/graphics/vectorize/test_utils.py
@@ -0,0 +1,34 @@
+import numpy as np
+from pytest import mark
+
+from pyquibbler.quib.graphics.axiswise_function_quibs.vectorize.utils import copy_vectorize, alter_signature
+from pyquibbler.quib.graphics.axiswise_function_quibs.vectorize.vectorize_metadata import VectorizeCall
+
+
+@mark.parametrize(['args', 'kwargs'], [((np.zeros((2, 2, 2)),), dict(y=np.zeros((2, 1, 2))))])
+@mark.parametrize(['func', 'signature', 'expected_core_ndims', 'expected_result_ndims', 'arg_ids_to_new_core_ndims'], [
+    (None, '(x),(y,z)->()', [1, 2], 0, {}),
+    (None, '(),()->(x),()', [0, 0], [1, 0], {}),
+    (lambda x, y: (x, y), None, [0, 0], [0, 0], {}),
+    (lambda x, y: x + y, None, [0, 0], 0, {}),
+    (None, '(x),(y,z)->()', [0, 2], 0, {0: 0}),
+    (None, '(x),(y,z)->()', [1, 1], 0, {'y': 1}),
+    (None, '(x),(y,z)->()', [3, 3], 0, {0: 3, 'y': 3}),
+    (lambda x, y: (x, y), None, [1, 0], [0, 0], {0: 1}),
+    (lambda x, y: (x, y), None, [0, 2], [0, 0], {'y': 2}),
+    (lambda x, y: (x, y), None, [3, 3], [0, 0], {0: 3, 'y': 3}),
+])
+def test_alter_signature(func, signature, args, kwargs, expected_core_ndims, expected_result_ndims,
+                         arg_ids_to_new_core_ndims):
+    vectorize = np.vectorize(func, signature=signature)
+    metadata = VectorizeCall(vectorize, args, kwargs).get_metadata()
+    altered_signature = alter_signature(metadata.args_metadata, metadata.result_or_results_core_ndims,
+                                        arg_ids_to_new_core_ndims)
+    vectorize_with_new_signature = copy_vectorize(vectorize, signature=altered_signature)
+    metadata_with_new_signature = VectorizeCall(vectorize_with_new_signature, args, kwargs).get_metadata()
+
+    assert metadata_with_new_signature.is_result_a_tuple == metadata.is_result_a_tuple
+    assert [meta.core_ndim for meta in metadata_with_new_signature.args_metadata.values()] == expected_core_ndims
+    assert (metadata_with_new_signature.results_core_ndims == expected_result_ndims) \
+        if metadata.is_result_a_tuple else \
+        (metadata_with_new_signature.result_core_ndim == expected_result_ndims)
diff --git a/tests/functional/quib/graphics/vectorize/test_vectorize_graphics_function_quib.py b/tests/functional/quib/graphics/vectorize/test_vectorize_graphics_function_quib.py
new file mode 100644
index 000000000..ed7d3e025
--- /dev/null
+++ b/tests/functional/quib/graphics/vectorize/test_vectorize_graphics_function_quib.py
@@ -0,0 +1,175 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from functools import partial
+from pytest import mark, fixture
+
+from pyquibbler import iquib, CacheBehavior
+from pyquibbler.env import GRAPHICS_LAZY
+from pyquibbler.quib import ProxyQuib
+from pyquibbler.quib.assignment import PathComponent
+
+from ...utils import check_invalidation, check_get_value_valid_at_path, MockQuib, PathBuilder, get_func_mock
+
+parametrize_data = mark.parametrize('data', [np.arange(24).reshape((2, 3, 4))])
+parametrize_indices_to_invalidate = mark.parametrize('indices_to_invalidate',
+                                                     [-1, 0, (0, 0), (0, 1, 2), (0, ...), [True, False]])
+
+
+@fixture
+def temp_axes():
+    ax = plt.gca()
+    yield ax
+    ax.clear()
+
+
+@parametrize_indices_to_invalidate
+@parametrize_data
+@mark.parametrize('excluded', [{0}, set(), None])
+@mark.parametrize('func', [lambda x: np.sum(x), lambda x: (np.sum(x), np.sum(x))])
+def test_vectorize_invalidation(indices_to_invalidate, data, excluded, func):
+    kwargs = {}
+    if excluded is not None:
+        kwargs['excluded'] = excluded
+    check_invalidation(np.vectorize(func, **kwargs), data, indices_to_invalidate)
+
+
+@parametrize_data
+@mark.parametrize('func', [lambda x: np.sum(x)])
+@mark.parametrize('indices_to_get_value_at', [0, (0, 0), (-1, ...)])
+def test_vectorize_get_value_valid_at_path(data, func, indices_to_get_value_at):
+    path_to_get_value_at = [PathComponent(np.ndarray, indices_to_get_value_at)]
+    check_get_value_valid_at_path(np.vectorize(func), data, path_to_get_value_at)
+
+
+@mark.parametrize('pass_quibs', [True, False])
+def test_vectorize_get_value_valid_at_path_with_excluded_quib(pass_quibs):
+    excluded = MockQuib(np.array([1, 2, 3]))
+
+    @partial(np.vectorize, excluded={1}, signature='(n)->(m)', pass_quibs=pass_quibs)
+    def func(_a, b):
+        if pass_quibs:
+            assert isinstance(b, ProxyQuib) and b._quib is excluded
+        return b
+
+    fquib = func([0, 1], excluded)
+    fquib.set_cache_behavior(CacheBehavior.OFF)
+    path = PathBuilder(fquib)[0].path
+
+    with excluded.collect_valid_paths() as valid_paths:
+        fquib.get_value_valid_at_path(path)
+
+    assert valid_paths == [[]]
+
+
+@mark.parametrize(['quib_data', 'non_quib_data'], [
+    (np.zeros((2, 3)), np.zeros((4, 2, 3))),
+    (np.zeros((3, 3)), np.zeros((1, 3))),
+    (np.zeros((1, 3, 3)), np.zeros((4, 1, 3))),
+    (np.zeros((4, 2, 3)), np.zeros((2, 3))),
+])
+def test_vectorize_get_value_valid_at_path_when_args_have_different_loop_dimensions(quib_data, non_quib_data):
+    func = lambda quib: np.vectorize(lambda x, y: x + y)(quib, quib_data)
+    check_get_value_valid_at_path(func, non_quib_data, [PathComponent(np.ndarray, 0)])
+
+
+@mark.parametrize('indices_to_get_value_at', [0, (1, 1), (1, ..., 2)])
+def test_vectorize_get_value_at_path_with_core_dims(indices_to_get_value_at):
+    quib_data = np.zeros((2, 2, 3, 4))
+    non_quib_data = np.zeros((2, 3, 4, 5))
+    func = lambda a, b: np.array([np.sum(a) + np.sum(b)] * 6)
+    vec = np.vectorize(func, signature='(a,b),(c)->(d)')
+    check_get_value_valid_at_path(lambda quib: vec(non_quib_data, quib), quib_data,
+                                  [PathComponent(np.ndarray, indices_to_get_value_at)])
+
+
+def test_vectorize_invalidation_with_non_numpy_return_value():
+    vec = np.vectorize(lambda a: int(np.sum(a)), signature='(a)->()')
+    check_invalidation(lambda quib: vec(quib), [1, 2, 3], 0)
+
+
+@parametrize_indices_to_invalidate
+@parametrize_data
+def test_vectorize_invalidation_with_multiple_params(data, indices_to_invalidate):
+    quib_arg = iquib(np.arange(12).reshape((3, 4)))
+    non_quib_arg = np.arange(8).reshape((2, 1, 4))
+    quib_kwarg = iquib(3)
+    non_quib_kwarg = np.arange(3).reshape((3, 1))
+
+    vec = np.vectorize(lambda a, b, c, d, e=0: a + b + c + d + e)
+    check_invalidation(lambda quib: vec(quib, quib_arg, non_quib_arg, d=quib_kwarg, e=non_quib_kwarg),
+                       data, indices_to_invalidate)
+
+
+@mark.parametrize('excluded', [{1}, {'b'}, set(), None])
+def test_vectorize_invalidation_with_quib_as_kwarg(excluded):
+    kwargs = {}
+    if excluded is not None:
+        kwargs['excluded'] = excluded
+    vec = np.vectorize(lambda a, b: a + b)
+    check_invalidation(lambda quib: vec(5, b=quib), [1, 2, 3], 0)
+
+
+@parametrize_indices_to_invalidate
+@parametrize_data
+def test_vectorize_invalidation_with_different_core_dims(data, indices_to_invalidate):
+    data2 = np.arange(100, 122).reshape(1, 22)
+    func = lambda a, b: (np.array([np.sum(a) + np.sum(b)] * 4), np.sum(a) + np.sum(b))
+    vec = np.vectorize(func, signature='(a,b),(c)->(d),()')
+    check_invalidation(lambda quib: vec(quib, data2), data, indices_to_invalidate)
+
+
+def test_vectorize_partial_calculation():
+    def f(x):
+        return x
+
+    func_mock = get_func_mock(f)
+    with GRAPHICS_LAZY.temporary_set(True):
+        quib = np.vectorize(func_mock)(iquib(np.arange(3)))
+    # Should call func_mock twice
+    quib.get_value_valid_at_path(PathBuilder(quib)[0].path)
+    assert func_mock.call_count == 2, func_mock.mock_calls
+    # Should call func_mock one time
+    quib.get_value_valid_at_path(PathBuilder(quib)[1].path)
+    assert func_mock.call_count == 3, func_mock.mock_calls[2:]
+    # Should not call func_mock
+    quib.get_value_valid_at_path(PathBuilder(quib)[0].path)
+    assert func_mock.call_count == 3, func_mock.mock_calls[3:]
+
+
+def test_vectorize_get_value_valid_at_path_none():
+    quib = np.vectorize(lambda x: x)(iquib([1, 2, 3]))
+    value = quib.get_value_valid_at_path(None)
+    assert len(value) == 3
+
+
+def test_vectorize_with_pass_quibs():
+    @partial(np.vectorize, pass_quibs=True)
+    def vectorized(x):
+        return iquib(x.get_value() + 1)
+
+    result = vectorized(iquib(np.arange(2)))
+    assert np.array_equal(result.get_value(), [1, 2])
+
+
+def test_vectorize_with_pass_quibs_and_core_dims():
+    @partial(np.vectorize, pass_quibs=True, signature='(a)->(x)')
+    def vectorized(x):
+        return iquib(x.get_value() + 1)[:-1]
+
+    result = vectorized(iquib(np.zeros((2, 3))))
+    assert np.array_equal(result.get_value(), np.ones((2, 2)))
+
+
+@mark.parametrize('pass_quibs', [True, False])
+def test_vectorize_does_not_redraw_valid_artists(temp_axes, pass_quibs):
+    parent = iquib([[1, 2], [3, 4]])
+    vectorized_plot = np.vectorize(plt.plot, signature='(x)->()', otypes=[np.object], pass_quibs=pass_quibs)
+    vectorized_plot(parent)
+    assert len(temp_axes.lines) == 2
+    ids = list(map(id, temp_axes.lines))
+
+    parent[0] = [5, 6]
+
+    assert len(temp_axes.lines) == 2
+    assert id(temp_axes.lines[0]) != ids[0]
+    assert id(temp_axes.lines[1]) == ids[1]
diff --git a/tests/functional/quib/test_quib.py b/tests/functional/quib/test_quib.py
index e20ce03c6..55096248d 100644
--- a/tests/functional/quib/test_quib.py
+++ b/tests/functional/quib/test_quib.py
@@ -1,6 +1,6 @@
 import operator
 from math import trunc, floor, ceil
-from typing import Set
+from typing import Set, Optional, List
 from unittest import mock
 from unittest.mock import Mock
 import numpy as np
@@ -58,7 +58,7 @@ def test_quib_invalidate_and_redraw_calls_graphics_function_quib_children(exampl
     mock_func = mock.Mock()
     mock_func.return_value = []
 
-    quib = GraphicsFunctionQuib(func=mock_func, args=tuple(), kwargs={}, artists=[], cache_behavior=None)
+    quib = GraphicsFunctionQuib(func=mock_func, args=tuple(), kwargs={}, cache_behavior=None)
     example_quib.add_child(quib)
     quib.get_value()  # we want to set cache to valid
 
@@ -235,7 +235,7 @@ def test_quib_get_shape():
     arr = np.array([1, 2, 3])
     quib = ExampleQuib(arr)
 
-    assert quib.get_shape().get_value() == arr.shape
+    assert quib.get_shape() == arr.shape
 
 
 @mark.parametrize(['data', 'overrides', 'expected_mask'], [
@@ -360,7 +360,7 @@ def __init__(self, parent):
         def parents(self) -> Set[Quib]:
             return {self._parent}
 
-        def _invalidate_self(self, path):
+        def _invalidate_self(self, path: List[PathComponent]):
             for parent in self.parents:
                 parent.remove_child(self)
 
diff --git a/tests/functional/quib/utils/__init__.py b/tests/functional/quib/utils/__init__.py
index 03b3048b9..399897698 100644
--- a/tests/functional/quib/utils/__init__.py
+++ b/tests/functional/quib/utils/__init__.py
@@ -1,3 +1,3 @@
 from .check_get_value_at_path import check_get_value_valid_at_path, MockQuib
 from .check_invalidation import check_invalidation
-from .utils import PathBuilder, slicer, get_mock_with_repr
+from .utils import PathBuilder, slicer, get_mock_with_repr, get_func_mock
diff --git a/tests/functional/quib/utils/check_invalidation.py b/tests/functional/quib/utils/check_invalidation.py
index d260eb17e..db4921fe4 100644
--- a/tests/functional/quib/utils/check_invalidation.py
+++ b/tests/functional/quib/utils/check_invalidation.py
@@ -11,7 +11,7 @@
 def breakdown_quib(quib: Quib) -> Set[Quib]:
     quib_type = quib.get_type()
     if issubclass(quib_type, np.ndarray):
-        return {quib[idx] for idx in np.ndindex(quib.get_shape().get_value())}
+        return {quib[idx] for idx in np.ndindex(quib.get_shape())}
     if issubclass(quib_type, (np.generic, int)):
         return set()
     if issubclass(quib_type, (tuple, list)):
@@ -35,5 +35,5 @@ def check_invalidation(func, data, indices_to_invalidate):
     input_quib[indices_to_invalidate] = 999
 
     invalidated_children = {child for child in children if child.cache_status == CacheStatus.ALL_INVALID}
-    changed_children = {child for child in children if child.get_value() != original_values[child]}
+    changed_children = {child for child in children if not np.array_equal(child.get_value(), original_values[child])}
     assert invalidated_children == changed_children
diff --git a/tests/functional/quib/utils/utils.py b/tests/functional/quib/utils/utils.py
index 50675a0c6..17f957e42 100644
--- a/tests/functional/quib/utils/utils.py
+++ b/tests/functional/quib/utils/utils.py
@@ -1,3 +1,4 @@
+from unittest import mock
 from unittest.mock import Mock
 from dataclasses import dataclass, field
 from typing import List
@@ -22,3 +23,7 @@ class PathBuilder:
 
     def __getitem__(self, item):
         return PathBuilder(self.quib[item], [*self.path, PathComponent(self.quib.get_type(), item)])
+
+
+def get_func_mock(func):
+    return mock.create_autospec(func, side_effect=func)