Repository: ekwan/PyPlate (language: python)

Issue:
Experiment design



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index f89675c..1c07649 100644
--- a/.gitignore
+++ b/.gitignore
@@ -160,4 +160,5 @@ cython_debug/
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
 .idea/
 
-docs/
\ No newline at end of file
+docs/
+.DS_Store
diff --git a/API.md b/API.md
index cbabe7f..d961bda 100644
--- a/API.md
+++ b/API.md
@@ -14,12 +14,20 @@ The *PyPlate* Python API defines a set of objects and operations for implementin
 
 Four simple HTE classes will be exposed to the user: `Substance`, `Container`, `Plate`, and `Recipe`.  *All classes are immutable.*  (An immutable object is one whose fields cannot be changed once it has been constructed.)
 
+Note: all quantity, volume, and max_volume parameters are given as strings. For example, '10 mL', '5 g', '1 mol', or '11 U'.
+
+The following are set based on preferences read `pyplate.yaml`:
+
+  - Units in which moles and volumes are stored internally. `moles_storage` and `volume_storage`
+  - Density of solids in g/mL. `solid_density`
+  - Units for '%w/v' concentrations ('g/mL'). `default_weight_volume_units:`
+
 ---
 
 ### Substance
 
 #### Definition: 
-- An abstract chemical or biological entity (e.g., reagent, enzyme, solvent, etc.).  Immutable.  Solids and enzymes are assumed to require zero volume.
+- An abstract chemical or biological entity (e.g., reagent, enzyme, solvent, etc.).  Immutable.  Solids and enzymes are assumed to have a density of 1.0.
 
 #### Constructors/Factory Methods
 
@@ -53,21 +61,36 @@ Four simple HTE classes will be exposed to the user: `Substance`, `Container`, `
 #### Attributes:
 
 - name (str): Name of this Container
-- contents (dict): map from `Substances` to amounts. Amounts are stored as mL for liquids, moles for solids, and U for enzymes
-- max_volume (float): in microlitres.
+- contents (dict): map from `Substances` to amounts. Amounts are stored in liters for liquids, moles for solids, and activity units for enzymes
+- max_volume (float): in storage format (determined by volume_storage from `pyplate.yaml`).
+
+#### Methods
 
+- has_liquid():
+  - Returns true if any substance in the container is a liquid.
+- remove(what):
+  - Removes substances from this Container. Defaults to removing all liquids.
+- dilute(solute, concentration, solvent, new_name)
+  - Creates a new diluted solution with respect to `solute`
+  - Concentration can be any of "0.1 M", "0.1 m", "0.1 g/mL", "0.01 umol/10 uL", "5 %v/v", "5 %w/v", "5 %w/w"
+  - Name of new container is optionally set to `new_name`
+- fill_to(solvent, quantity)
+  - Returns new container filled with `solvent` up to `quantity`.
+- get_concentration(solute, units)
+  - Returns the current concentration of `solute` in `units`.
 
 #### Static Methods:
 
-- add(source, destination, volume):
-  - Move the given quantity of the *source* substance to the *destination* container. A new copy of *destination* will be returned.
-- transfer(source, destination, how_much):
-  - Move *volume* from *source* to *destination* container, returning copies of the objects with amounts adjusted accordingly.
-  - Note that all `Substances` in the source will be transferred in proportion to their volumetric ratios.
+- `transfer(source, destination, quantity)`:
+  - Move *quantity* from *source* to *destination* container, returning copies of the objects with amounts adjusted accordingly.
+  - Note that all `Substances` in the source will be transferred in proportion to their appropriate ratios.
   - *source* can be a container, a plate, or a slice of a plate.
-- create_stock_solution(what, concentration, solvent, volume)
-  - Create a new container with the desired volume and containing the desired concentration of `what`.
-  - If `what` is a liquid, volumes will be calculated appropriately.
+- `create_solution(solute, solvent, name, concentration?, quantity?, total_quantity?)`
+  - Create a new container with the desired solution based on given arguments.
+  - Two of concentration, quantity, total_quantity must be specified
+  - Concentration can be any of "0.1 M", "0.1 m", "0.1 g/mL", "0.01 umol/10 uL", "5 %v/v", "5 %w/v", "5 %w/w"
+  - If `solute` is a liquid, volumes will be calculated appropriately.
+  - name is optional. If none is given, an appropriate name will be applied.
 
 ---
 
@@ -79,13 +102,14 @@ Four simple HTE classes will be exposed to the user: `Substance`, `Container`, `
 
 #### Constructors/Factory Methods:
 
-- Plate(name, max_capacities, n_rows, n_columns, row_labels, column_labels)
-  - name (str): name of this Container
-  - max_capacities (float): assumed to be the same volume for all wells, in uL
-  - n_rows (int): how many rows there are (default, 8)
-  - n_cols (int): how many columns there are (default, 12)
-  - row_labels (list of str): by default, A, B, C, ..., AA, AB, ...
-  - column_labels (list of str): by default, "1", "2", ...
+- Plate(name, max_volume_per_well, make, rows, cols)
+  - name (str): name of this plate
+  - max_capacities (str): assumed to be the same volume for all wells
+  - make (str): name of this kind of plate
+  - rows (int or list of str): Either how many rows there are or labels for the rows
+  - cols (int or list of str): Either how many columns there are or labels for the columns
+  - Row names default to "A", "B", ..., "AA", "AB", etc.
+  - Column names default to "1", "2", etc.
 
 
 #### Methods:
@@ -93,21 +117,32 @@ Four simple HTE classes will be exposed to the user: `Substance`, `Container`, `
 - Plate[slice]
   - Returns a slice of the plate.
 
-- volumes()
-  - Returns a `numpy` array of used volumes
+- remove(what):
+  - Removes substances from all wells in this plate. Defaults to removing all liquids.
 - substances()
   - Returns a set of all substances used
-- moles(substance)
+- volumes(substance, unit)
+  - Returns a `numpy` array of used volumes
+  - If substance is given, volumes will be restricted to volumes of substance
+  - If unit is given, volumes will be given in unit, otherwise in `default_volume_unit` defined in `pyplate.yaml`
+- volumes_dataframe(substance, unit, cmap)
+  - Returns a shaded dataframe of volumes in each well
+  - Unit defaults to `default_volume_unit` defined in `pyplate.yaml`
+  - cmap defaults to `default_colormap` defined in `pyplate.yaml`
+- moles(substance, unit)
   - Returns a `numpy` array of moles of given substance
-- concentrations(substance)
-
-** Volumes, substances, moles, and concentrations can all be called on a slice of the plate
+  - If unit is given, moles will be return in unit, otherwise in `default_moles_unit` defined in `pyplate.yaml`
+- moles_dataframe(substance, unit, cmap)
+  - Returns a shaded dataframe of moles in each well
+  - Unit defaults to `default_moles_unit` defined in `pyplate.yaml`
+  - cmap defaults to `default_colormap` defined in `pyplate.yaml`
+  - 
+** Remove, substances, volumes, and moles can all be called on a slice of the plate
 
 
 #### Static Methods:
 
-- add(source, destination, volume): Move the given quantity of the *source* substance to the *destination* container. A new copy of *destination* will be returned.
-- transfer(source, destination, how_much): Move *volume* from *source* to *destination* plate or slice, returning copies of the objects with amounts adjusted accordingly.
+- transfer(source, destination, volume): Move *volume* from *source* to *destination* plate or slice, returning copies of the objects with amounts adjusted accordingly.
   - Note that all `Substances` in the source will be transferred in proportion to their volumetric ratios.
   - *source* can be a container, a plate, or a slice of a plate.
 
@@ -131,19 +166,28 @@ uses (list): a list of *Containers* that will be used in this `Recipe`.  An exce
 
 - uses(*containers)
   - declare `*containers` (iterable of `Containers`) as being used in the recipe.
-- add(source, destination, volume):
-  - Adds a step to the recipe which will move the given quantity of the *source* substance to the *destination*.
-- transfer(source, destination, how_much):
-  - Adds a step to the recipe which will move *volume* from *source* to *destination*.
-  - Note that all `Substances` in the source will be transferred in proportion to their volumetric ratios.
+- transfer(source, destination, quantity):
+  - Adds a step to the recipe which will move *quantity* from *source* to *destination*.
+  - Note that all `Substances` in the source will be transferred in proportion to their respective ratios.
 - create_container(name, max_volume, initial_contents)
   - Keep track of steps to create container in recipe
   - Adds a step that creates a container as above and adds it to the used list.
   - Returns new container so that it can be used later in the same recipe.
-- create_stock_solution(what, concentration, solvent, volume)
-  - Adds a step to the recipe with will create a new container with the desired volume and containing the desired concentration of `what`.
-  - If `what` is a liquid, volumes will be calculated appropriately.
+- create_solution(solute, solvent, name, concentration?, quantity?, total_quantity?)
+  - Adds a step to the recipe with will create a new container with the desired solution based on given arguments.
+  - Two of concentration, quantity, total_quantity must be specified
+  - Concentration can be any of "0.1 M", "0.1 m", "0.1 g/mL", "0.01 umol/10 uL", "5 %v/v", "5 %w/v", "5 %w/w"
+  - If `solute` is a liquid, volumes will be calculated appropriately.
+  - name is optional. If none is given, an appropriate name will be applied.
   - Returns new container so that it can be used later in the same recipe.
+- dilute(destination, solute, concentration, solvent, new_name)
+  - Adds a step to create a new container diluted to a certain `concentration` of `solute` from `destination`
+  - Concentration can be any of "0.1 M", "0.1 m", "0.1 g/mL", "0.01 umol/10 uL", "5 %v/v", "5 %w/v", "5 %w/w"
+  - Name of new container is optionally set to `new_name`
+- fill_to(destination, solvent, quantity)
+  - Adds a step to fill `destination` container with `solvent` up to `quantity`.
+- remove(destination, what)
+  - Adds a step to removes substances from destination. Defaults to removing all liquids.
 - bake()
   - Checks the validity of each step and ensures all Containers are used.
   - Returns all new Containers and Plates in the order they were defined in `uses()`.
diff --git a/Experimental-Design-Proposal.md b/Experimental-Design-Proposal.md
new file mode 100644
index 0000000..555f306
--- /dev/null
+++ b/Experimental-Design-Proposal.md
@@ -0,0 +1,105 @@
+﻿
+# Proposal for PyPlate Experimental Design
+
+## Factor
+
+Represents an experimental condition that is under the experimenter's control.
+
+### Instance Variables:
+
+- `name`: A string representing the name of the factor, must be unique within an Experimental Space.
+
+- `possible_values`: A list of strings, numerics, or Substances representing the possible values of the factor.
+
+- `verifier`: A function that takes a Well object and a possible value. It returns a boolean value if the contents of the well match the desired value provided. In addition, it verifies the provided value is in `possible_values`.
+
+## Experiment
+
+Represents a single experiment.
+
+### Instance Variables:
+
+- `factors`: A dictionary with factor names as keys and corresponding values.
+
+- `experiment_id`: Some arbitrary identifier for the experiment, may be auto generated or specified by the user.
+
+- `replicate_idx`: An identifier to distinguish between repeated experiments with the same factors, by default is `experiment_id-{numeric idx starting from 1}`. 
+
+- `well`: A pointer to the `Well` an experiment corresponds to
+
+### Methods
+
+-`map_well(Well)`: Maps a given `Well` to this `Experiment`.
+
+
+## ExperimentalSpace
+
+Represents the experimental space, including fixed and variable factors.  
+
+### Instance Variables:
+
+- `factors`: A list of Factor objects representing the fixed and variable factors in the experimental space.
+
+- `blocks`: A dictionary with names as keys and values as lists containing experiments representing the generated experiments.
+-`experiments`: A list of all experiments registered with the space
+- `experiment_id_generator`: A function that returns an identifier for a given Experiment in the space. The function will take in an `Experiment` object, and return a string that will serve as the ID for a given experiment.
+
+### Methods:
+
+- `generate_experiments(factors, n_replicates, block_on)`: Fully enumerates the tree of variable factors, generates experiments with well/plate indices that respect blocking rules, and stores them in the instance variable. `n_replicates` is number of replicates to be generated. If experiments with this combination of factors exist, `n_replicates` is the number of *additional* replicates to generated. `block_on` is a list of sets of factor names to block on (analogous to `GROUP BY` in SQL). Automatically adds all experiments to the space.
+
+- `is_valid(Experiment)`: Checks if a given experiment is valid. If `generate_experiments` has run, simply checks if the experiment is in the list. Otherwise, manually checks if values are in a Factor's possible values.
+
+- `add_experiment(Experiment)`: Manually adds an experiment to the space. Checks for validity of the experiment.
+
+  
+
+## Example Usage
+
+```python
+
+# Create Factors for all experimental conditions
+
+f1 = Factor(name="Baking Time", possible_values=[10, 20, 30], verifier=verify_time)
+f2 = Factor(name="Baking Temperature", possible_values=[300, 350, 400], verifier=verify_temperature)
+f3 = Factor(name="Flavor", possible_values=["Chocolate Chip", "Oatmeal Raisin", "Peanut Butter"], verifier=verify_flavor)
+
+# Create Experimental Space
+space = ExperimentalSpace(factors=[f1, f2, f3])
+space.generate_experiments(
+    factors={
+        "Baking Time": [10, 30],
+        "Baking Temperature": [300, 350],
+        "Flavor": "all"
+    },
+    n_replicates=2,
+    blocking_factors=["Baking Temperature"]
+)
+
+# returns:
+[
+    # block of 300 baking temperature
+    [
+        Experiment(factors={"Baking Time": 10, "Baking Temperature": 300, "Flavor": "Chocolate Chip"}),
+        Experiment(factors={"Baking Time": 10, "Baking Temperature": 300, "Flavor": "Oatmeal Raisin"}),
+        Experiment(factors={"Baking Time": 10, "Baking Temperature": 300, "Flavor": "Peanut Butter"}),
+        Experiment(factors={"Baking Time": 30, "Baking Temperature": 300, "Flavor": "Chocolate Chip"}),
+        Experiment(factors={"Baking Time": 30, "Baking Temperature": 300, "Flavor": "Oatmeal Raisin"}),
+        Experiment(factors={"Baking Time": 30, "Baking Temperature": 300, "Flavor": "Peanut Butter"})
+    ],
+    
+    # block of 350 baking temperature
+    [
+        Experiment(factors={"Baking Time": 10, "Baking Temperature": 350, "Flavor": "Chocolate Chip"}),
+        Experiment(factors={"Baking Time": 10, "Baking Temperature": 350, "Flavor": "Oatmeal Raisin"}),
+        Experiment(factors={"Baking Time": 10, "Baking Temperature": 350, "Flavor": "Peanut Butter"}),
+        Experiment(factors={"Baking Time": 30, "Baking Temperature": 350, "Flavor": "Chocolate Chip"}),
+        Experiment(factors={"Baking Time": 30, "Baking Temperature": 350, "Flavor": "Oatmeal Raisin"}),
+        Experiment(factors={"Baking Time": 30, "Baking Temperature": 350, "Flavor": "Peanut Butter"})
+        
+    ]
+]
+
+# Add a single experiment to the space manually
+space.add_experiment(Experiment(factors={"Baking Time": 20, "Baking Temperature": 350, "Flavor": "Chocolate Chip"}))
+```
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..f490053
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,15 @@
+PYTHONPATH := $(shell pwd)
+
+docs: pyplate/pyplate.py
+	pdoc pyplate/pyplate.py -d google -o docs/ --no-include-undocumented
+
+test:
+	pytest tests
+
+coverage:
+	pytest --cov=pyplate --cov-report html tests
+
+all: docs coverage
+
+clean:
+	rm -rf docs htmlcov
\ No newline at end of file
diff --git a/README.md b/README.md
index d6a3d73..ce8ac31 100644
--- a/README.md
+++ b/README.md
@@ -1,74 +1,209 @@
-# PyPlate: Substance and Mixture Conversion Library
+# PyPlate
 
-PyPlate is an innovative Python library designed for scientists, researchers, and enthusiasts working in chemistry and biochemistry. It provides a comprehensive set of tools to define substances, create mixtures, and handle complex conversions between different units, all within a well-structured plate format for experiments.
+An open-source Python tool for high-throughput chemistry.
 
-## Key Features
+### Introduction
 
-- **Substance Definition**: Easily define substances with properties such as state (solid, liquid, enzyme) and molecular weight. Each substance is treated as an immutable entity, ensuring consistency throughout your experiments.
-- **Mixture Creation**: Combine different substances into mixtures and specify quantities in various units. PyPlate handles the complex backend calculations for these combinations.
-- **Unit Conversion**: Seamlessly convert between different units like grams, liters, moles, molar concentrations, and arbitrary units (AUs). This feature is invaluable for maintaining accuracy in experimental protocols.
-- **Plate Experiment Design**: Organize your substances and mixtures in a plate format, perfect for high-throughput screening or analysis. The library supports a spatially ordered collection of containers in a rectangular arrangement, commonly used in 96-well plates.
-- **Recipe Transformation**: Not just limited to static designs, PyPlate allows you to create dynamic recipes for transforming one set of containers into another, enabling advanced experimental workflows.
+PyPlate provides tools for the design and implementation of high-throughput chemistry experiments (in particular, reaction screening).  It allows the user to define a space of experimental parameters to be explored, select points in that space for experimentation, and design liquid/solid handling steps to implement those experiments in 96 well plates.
 
-## Installation
+### Installation
 
-Before you start using PyPlate, you need to install it. Here's how you can install it:
+PyPlate requires Python 3.10 or later.
 
-\```bash
-pip install pyplate
-\```
+`pip install pyplate-hte`
 
+To view plate visualizations, you will need an interactive Python shell like Jupyter.
 
+### Philosophy
 
-## Getting Started
+All experiments are divided into a *design* phase and an *implementation* phase.
 
-Dive into PyPlate by defining substances, creating mixtures, and more! Here's how:
+**Design Phase**: TBD
 
-### 1. Defining Substances
+**Implementation Phase**: PyPlate mimics the physical process of dispensing solids or liquids into plates.  `Substance`s are placed into `Container`s and dispensed into `Plate`s.  The instructions for creating a set of plates are grouped into `Recipe` objects.
+
+- add something about how to check that the implementation and design are consistent
+- All objects in PyPlate are considered immutable.
+
+## Quick Start
 
 ```python
-from pyplate import Substance
 
-# Define a solid substance
-sodium_chloride = Substance.solid(name='Sodium Chloride', mol_weight=58.44)
+from pyplate import Substance, Container, Plate, Recipe
+
+triethylamine = Substance.liquid(name="triethylamine", mol_weight=101.19, density=0.726)
+DMSO = Substance.liquid(name="DMSO", mol_weight=78.13, density=1.1004)
+
+triethylamine_50mM = Container.create_solution(solute=triethylamine, solvent=DMSO, concentration='50 mM', total_quantity='10 mL')
+plate = Plate(name='plate', max_volume_per_well='50 uL')
+
+recipe = Recipe().uses(triethylamine_50mM, plate)
+recipe.transfer(source=triethylamine_50mM, destination=plate[2:7, 2:11], quantity='10 uL')
+triethylamine_50mM, plate = recipe.bake()
 
-# Define a liquid substance
-water = Substance.liquid(name='Water', mol_weight=18.01528, density=1.0)
+recipe.steps[-1].visualize(what='destination', mode='final', unit='uL')[0]
 
-# Define an enzyme substance
-enzyme = Substance.enzyme(name='Enzyme X')
 ```
+![img.png](images/simple_visualization.png)
+## User Guide
 
-### 2. Creating Mixtures
+### Substances
+
+A `Substance` is a solid, liquid, or enzyme:
 
 ```python
-from pyplate import Mixture
+# solid
+sodium_chloride = Substance.solid(name="Sodium Chloride", mol_weight=58.44)
+
+# liquid
+water = Substance.liquid(name="Water", mol_weight=18.01528, density=1.0)
+
+# enzyme
+enzyme_X = Substance.enzyme(name="Enzyme X")
+```
 
-# Create a new mixture named 'Solution 1'
-solution_1 = Mixture('Solution 1')
+Enzymes are like solids, but are specified in *units of activity* instead of moles.  The default density for substances is 1.0.
 
-# Add substances to the mixture
-solution_1 += (sodium_chloride, '10 g')  # Adding 10 grams of Sodium Chloride
-solution_1 += (water, '1 L')  # Adding 1 liter of Water
+### Containers
+
+A `Container` holds defined amounts of a set of `Substance`s and has a maximum volume.
+
+Placing 10 g of NaCl in a 1 L container:
+
+```python
+sodium_chloride_stock = Container(name = "sodium_chloride_stock",
+                                  max_volume = "1 L",
+                                  initial_contents =
+                                      [(sodium_chloride, "10 g")])
 ```
 
-### 3. Converting Units
+Making 10 mL of a 0.01 M stock solution of NaCl in water:
 
 ```python
-# Convert grams to moles for solid substance
-moles_of_sodium_chloride = sodium_chloride.convert_to_unit_value('10 g')
+sodium_chloride_10mM = Container.create_solution(
+                           solute = sodium_chloride,
+                           solvent = water,
+                           concentration = "0.01 M",
+                           total_quantity = "10 mL")
+```
 
-# Convert liters to milliliters for liquid substance
-milliliters_of_water = water.convert_to_unit_value('1 L')
+Diluting the above solute to 0.005 M:
+
+```python
+sodium_chloride_10mM, sodium_chloride_5mM = Container.create_solution_from(source=sodium_chloride_10mM, solute=sodium_chloride, concentration="5 mM", solvent=water, quantity="10 mL")
 ```
 
-## Configuration
+### Manipulating Containers
 
-PyPlate's functionality can be enhanced with a configuration file. For example, you can specify how volumes and moles are stored using a "pyplate.yaml" file.
+`Container`s can be diluted to a desired concentration or filled to a desired volume.
 
-```yaml
-volume_storage: 'uL'  # store volumes in microliters
-moles_storage: 'umol'  # store moles in micromoles
+```python
+salt_water = salt_water.dilute(solute=sodium_chloride_10mM, concentration='5 mM', solvent=water)
+salt_water = salt_water.fill_to(solvent=water, quantity='50 mL')
 ```
 
-Adjust `config.internal_precision` for computation rounding to suit the precision needs of your experiments.
+### Creating a Plate
+
+`Plate`s are rectangular arrays of `Container`s.
+
+To create a generic 96 well plate:
+
+`plate = Plate(name="test plate", max_volume_per_well="50 uL")`
+
+Custom plates can be created with different numbers of rows and columns as well as different labels:
+
+`plate = Plate(name="custom plate", max_volume_per_well="50 uL", rows=16, columns=24)`
+`plate = Plate(name="custom plate", max_volume_per_well="50 uL", rows=['i', 'ii', 'iii'], columns=['a', 'b', c'])`
+
+### Locations on a Plate and slices
+
+PyPlate follows the `pandas` convention of having both integer- and label-based indices for referencing wells in `Plate`s.  When row or column specifiers are provided as integers, they are assumed to be integer indices (1, 2, 3, ...).  When specifiers are provided as strings, they are assumed to be label indices ("A", "B", "C", ...).
+
+By default, rows in plates are given alphabetical labels "A", "B", "C", ... and columns in plates are given numerical labels 1, 2, 3.  However, rows and columns are always given integer indices 1, 2, 3, ....  For example, "B:3" and (2,3) both refer to well B3.
+
+Here are some ways to refer to a specific well:
+
+- **String Method**: "A:1"
+- **Tuple Method**: ('A', 1)
+
+You can refer to multiple wells as a list:
+
+`plate[[('A', 1), ('B', 2), ('C', 3), 'D:4']]`
+
+
+Slicing syntax is supported:
+
+
+- In addition, you can provide python slices of wells with 1-based indexes:
+  - `plate[:3], plate[:, :3], plate['C':], plate[1, '3':]`
+
+### Recipes
+
+A `Recipe` is a set of instructions for transforming one set of containers into another.
+
+```python
+plate = Plate(name='plate', max_volume_per_well='50 uL')
+recipe = Recipe()
+```
+
+`Container`s and `Plate`s must be declared in the `Recipe` before use:
+
+```python
+recipe.uses(plate)
+```
+
+It can be convenient to create solutions and declare them for use in the same step:
+
+```python
+triethylamine_50mM = recipe.create_solution(solute=triethylamine, solvent=DMSO, concentration='0.05 M', total_quantity='10.0 mL')
+```
+
+(Note that solutions made in this way are not actually created until `recipe.bake()` is called.)  Performing transfer steps:
+
+```python
+recipe.transfer(source=triethylamine_50mM, destination=plate[:3], quantity='10 uL')
+```
+
+When `recipe.bake()` is called, the resulting `Container` and `Plate` objects are returned (leaving the input objects, if any, unchanged) in order of declaration or creation:
+
+```
+plate, triethylamine_50mM = recipe.bake()
+```
+
+
+Each operation called on a recipe is stored as a step. You can retrieve the steps for a recipe using `recipe.steps`.
+Each step has instructions as to what happened during the step (`step.instruction`) as well as visualizations.
+
+
+### Visualizations
+
+Each step in a recipe that involves a plate will yield visualizations. We can visualize information about the source or
+destination plate (as applicable) or both. Visualization can be determined for two different modes, either the 'final' state of the plate or a 'delta' of what changed during the step. Finally, the unit in which data is returned must be provided.
+
+```
+...
+last_step = recipe.steps[-1]
+for df in last_step.visualize(what='destination', mode='delta', unit='umol'):
+    display.display_html(df)
+```
+![Example Visualizations](images/example_visualization.png)
+
+### Units
+
+The basic units of pyplate are moles, grams, liters, and activity units. ('mol', 'g', 'L', 'U')
+
+Any time units are required, metric prefixes may be specified. ('mg', 'umol', 'dL', ...)
+
+All quantities are specified as strings with a value and a unit. ('1 mmol', '10 g', '10 uL' ...)
+
+### Building documentation
+
+Documentation can be built by executing `make docs`. The resulting documentation will be in the `docs` folder.
+
+### Running tests
+
+Tests can be run by executing `pytest`.
+## License
+
+Licensed under the Apache License, Version 2.0 (the "License")
+You may obtain a copy of the License at https://www.apache.org/licenses/LICENSE-2.0
\ No newline at end of file
diff --git a/examples/Example.py b/examples/Example.py
index edc7ef4..01114f6 100644
--- a/examples/Example.py
+++ b/examples/Example.py
@@ -6,30 +6,32 @@
 print("reagents:")
 sodium_sulfate = Substance.solid("sodium sulfate", mol_weight=142.04)
 triethylamine = Substance.liquid("triethylamine", mol_weight=101.19, density=0.726)
+water_tap = Substance.liquid("tap water", 18.0153, 1)
+water_DI = Substance.liquid("DI water", 18.0153, 1)
+DMSO = Substance.liquid("DMSO", 78.13, 1.1004)
 print(sodium_sulfate)
 print(triethylamine)
+print(water_tap)
+print(water_DI)
+print(DMSO)
 print()
 
 # create solvents
 print("solvents:")
-water_DI = Substance.liquid("DI water", 18.0153, 1)
-water_DI_stock = Container("DI water", max_volume=10, initial_contents=[(water_DI, "10 mL")])
-water_tap = Substance.liquid("tap water", 18.0153, 1)
-water_tap_stock = Container("tap water", max_volume=20, initial_contents=[(water_tap, "20 mL")])
-DMSO = Substance.liquid("DMSO", 78.13, 1.1004)
-DMSO_stock = Container("DMSO", max_volume=15, initial_contents=[(DMSO, "15 mL")])
-print(water_DI)
-print(water_tap)
-print(DMSO)
+water_DI_container = Container("DI water", max_volume='10 mL', initial_contents=[(water_DI, "10 mL")])
+water_tap_container = Container("tap water", max_volume='20 mL', initial_contents=[(water_tap, "20 mL")])
+DMSO_container = Container("DMSO", max_volume='15 mL', initial_contents=[(DMSO, "15 mL")])
+print(water_DI_container)
+print(water_tap_container)
+print(DMSO_container)
 print()
 
 # create stocks
 # concentrations in M
-# volumes in M
 print("stock solutions:")
-sodium_sulfate_halfM = Container.create_stock_solution(sodium_sulfate, 0.5, solvent=water_DI, volume=10.0)
-triethylamine_10mM = Container.create_stock_solution(triethylamine, 0.01, DMSO, volume=10.0)
-triethylamine_50mM = Container.create_stock_solution(triethylamine, 0.05, DMSO, volume=10.0)
+sodium_sulfate_halfM = Container.create_stock_solution(sodium_sulfate, 0.5, solvent=water_DI, volume='10.0 mL')
+triethylamine_10mM = Container.create_stock_solution(triethylamine, 0.01, DMSO, volume='10.0 mL')
+triethylamine_50mM = Container.create_stock_solution(triethylamine, 0.05, DMSO, volume='10.0 mL')
 print(sodium_sulfate_halfM)
 print(triethylamine_10mM)
 print(triethylamine_50mM)
@@ -37,7 +39,7 @@
 
 # create plate
 print("plate:")
-plate = Plate("test plate", max_volume_per_well=500.0)
+plate = Plate("test plate", max_volume_per_well='500.0 uL')
 print(plate)
 
 # add stuff to the plate
@@ -47,15 +49,17 @@
     recipe.transfer(sodium_sulfate_halfM, plate[k], f"{v} uL")
 for k, v in {"D:10": 1.0, (5, 10): 2.0, (5, "11"): 3.0}.items():
     recipe.transfer(triethylamine_10mM, plate[k], f"{v} uL")
-# recipe.transfer(triethylamine_10mM, plate["D:10"], "20_000 uL")
+# recipe.transfer(triethylamine_10mM, plate["D:10"], "20 mL")
 
 plate, sodium_sulfate_halfM, triethylamine_10mM = recipe.bake()
 
-print(plate.volumes())
+print('first recipe:')
+print('volumes in uL:')
+print(plate.volumes(unit='uL'))
 print(sodium_sulfate_halfM)
-# exit()
+print()
 
-recipe2 = Recipe().uses(plate, triethylamine_10mM, water_DI_stock, DMSO_stock)
+recipe2 = Recipe().uses(plate, triethylamine_10mM, water_DI_container, DMSO_container)
 for i in range(1, plate.n_columns+1):
     recipe2.transfer(triethylamine_10mM, plate[1, i], f"{30 * i} uL")
 recipe2.transfer(triethylamine_10mM, plate[6], "2 uL")
@@ -63,10 +67,14 @@
 recipe2.transfer(triethylamine_10mM, plate[:, "10"], "8 uL")
 recipe2.transfer(triethylamine_10mM, plate[:, 11:12], "9 uL")
 
-recipe2.transfer(water_DI_stock, plate[1:8], "20 uL")
-recipe2.transfer(DMSO_stock, plate[:, 1:12], "1 uL")
+recipe2.transfer(water_DI_container, plate[1:8], "20 uL")
+recipe2.transfer(DMSO_container, plate[:, 1:12], "1 uL")
 
-plate, triethylamine_10mM, water_DI_stock, DMSO_stock = recipe2.bake()
+plate, triethylamine_10mM, water_DI_container, DMSO_container = recipe2.bake()
+print('second recipe:')
+print('volumes in uL:')
+print(plate.volumes(unit='uL'))
+print(DMSO_container)
 print()
 
 # print total volumes
@@ -81,16 +89,35 @@
         print(f"{substance.name} : {plate.volumes(substance).sum():.1f} uL")
 print()
 
-exit()
 
 # print moles
 print("micromoles:\n")
-for i in range(len(plate.moles)):
-    print(plate.reagents[i])
-    print(plate.moles[i, :, :])
+for substance in plate.substances():
+    print(f"{substance.name}:")
+    print(plate.moles(substance, 'umol'))
     print()
 print()
 
+
+plate = plate.remove()
+print(plate.substances())
+print(plate.volumes(unit='mL'))
+
+triethylamine_10mM = Container.create_stock_solution(triethylamine, 0.01, DMSO, volume='10.0 mL')
+print(triethylamine_10mM)
+print("Diluting to 0.005 M")
+result = triethylamine_10mM.dilute(triethylamine, '0.005 M', DMSO)
+print(result)
+print("New concentration:", result.contents[triethylamine]/result.volume)
+
+sodium_sulfate_halfM = Container.create_stock_solution(sodium_sulfate, 0.5, solvent=water_DI, volume='10.0 mL')
+print(sodium_sulfate_halfM)
+result = sodium_sulfate_halfM.dilute(sodium_sulfate, '0.25 M', water_DI)
+print(result)
+print(result.contents[sodium_sulfate]/result.volume)
+
+exit()
+
 # dump plate to excel
 filename = "plate.xlsx"
 plate.to_excel(filename)
diff --git a/examples/Test_Example.py b/examples/Test_Example.py
index 3b59673..efc41ff 100644
--- a/examples/Test_Example.py
+++ b/examples/Test_Example.py
@@ -4,14 +4,14 @@
 salt = Substance.solid('NaCl', 58.4428)
 
 recipe = Recipe()
-recipe.create_container('halfM salt water', 1000, ((water, '100 mL'), (salt, '50 mmol')))
+recipe.create_container('halfM salt water', '1 L', ((water, '100 mL'), (salt, '50 mmol')))
 salt_water_halfM, = recipe.bake()
 
 triethylamine = Substance.liquid("triethylamine", mol_weight=101.19, density=0.726)
-triethylamine_halfM = Container.create_stock_solution(triethylamine, 0.5, water, 10)
-print(triethylamine_halfM)
+triethylamine_5mM = Container.create_stock_solution(triethylamine, 0.005, water, '10 mL')
+print(triethylamine_5mM)
 
-salt_water_oneM = Container('oneM salt water', 1000)
+salt_water_oneM = Container('oneM salt water', '1 L')
 recipe2 = Recipe().uses(water, salt, salt_water_oneM)
 recipe2.transfer(water, salt_water_oneM, '100 mL')
 recipe2.transfer(salt, salt_water_oneM, '100 mmol')
@@ -21,7 +21,7 @@
 recipe3.transfer(salt_water_halfM, salt_water_oneM, '5 mL')
 salt_water_halfM, salt_water_oneM = recipe3.bake()
 
-plate = Plate('plate', 50_000.0)
+plate = Plate('plate', '50 mL')
 recipe4 = Recipe().uses(plate, salt_water_oneM, salt_water_halfM)
 recipe4.transfer(salt_water_oneM, plate[:4], '0.5 mL')
 recipe4.transfer(salt_water_halfM, plate, '0.5 mL')
@@ -37,3 +37,11 @@
 print(plate.substances())
 print(salt_water_halfM)
 print(salt_water_oneM)
+
+recipe6 = Recipe()
+recipe6.create_stock_solution(salt, 0.5, water, '10 mL')
+salt_water2, = recipe6.bake()
+
+salt_water3 = salt_water2.remove()
+print(salt_water2)
+print(salt_water3)
diff --git a/images/example_visualization.png b/images/example_visualization.png
new file mode 100644
index 0000000..0db5855
Binary files /dev/null and b/images/example_visualization.png differ
diff --git a/images/simple_visualization.png b/images/simple_visualization.png
new file mode 100644
index 0000000..787b3e8
Binary files /dev/null and b/images/simple_visualization.png differ
diff --git a/pyplate.yaml b/pyplate.yaml
index 36b24d3..10fa910 100644
--- a/pyplate.yaml
+++ b/pyplate.yaml
@@ -6,4 +6,12 @@ external_precision: 3
 # 1e-6 means we will store volumes as microliters
 volume_storage: uL
 # 1e-6 means we will store moles as micromoles
-moles_storage: umol
\ No newline at end of file
+moles_storage: umol
+# density in g/mL
+solid_density: 1
+# units for %w/v
+default_weight_volume_units: g/mL
+default_moles_unit: mol
+default_volume_unit: uL
+default_colormap: Purples
+default_diverging_colormap: PuOr
\ No newline at end of file
diff --git a/pyplate/__init__.py b/pyplate/__init__.py
index 63d54e5..5cee6ed 100644
--- a/pyplate/__init__.py
+++ b/pyplate/__init__.py
@@ -2,6 +2,7 @@
 import os
 import yaml
 
+# import pyplate
 
 class Config:
     def __init__(self):
@@ -25,29 +26,16 @@ def __init__(self):
 
         self.internal_precision = yaml_config['internal_precision']
         self.external_precision = yaml_config['external_precision']
-        self.volume_storage = Config.convert_prefix_to_multiplier(yaml_config['volume_storage'][:-1])
-        self.moles_storage = Config.convert_prefix_to_multiplier(yaml_config['moles_storage'][:-3])
-
-    @staticmethod
-    def convert_prefix_to_multiplier(prefix: str) -> float:
-        """
-
-        Converts an SI prefix into a multiplier.
-        Example: "m" -> 1e-3, "u" -> 1e-6
-
-        Arguments:
-            prefix:
-
-        Returns:
-             Multiplier (float)
-
-        """
-        if not isinstance(prefix, str):
-            raise TypeError("SI prefix must be a string.")
-        prefixes = {'n': 1e-9, 'u': 1e-6, 'µ': 1e-6, 'm': 1e-3, 'c': 1e-2, 'd': 1e-1, '': 1, 'k': 1e3, 'M': 1e6}
-        if prefix in prefixes:
-            return prefixes[prefix]
-        raise ValueError(f"Invalid prefix: {prefix}")
-
-
-config = Config()
+        self.moles_prefix = yaml_config['moles_storage']
+        assert self.moles_prefix[1:] == 'mol'
+        self.volume_prefix = yaml_config['volume_storage']
+        assert self.volume_prefix[1:] == 'L'
+        self.solid_density = float(yaml_config['solid_density'])
+        self.default_weight_volume_units = yaml_config['default_weight_volume_units']
+        self.default_moles_unit = yaml_config['default_moles_unit']
+        self.default_volume_unit = yaml_config['default_volume_unit']
+        self.default_colormap = yaml_config['default_colormap']
+        self.default_diverging_colormap = yaml_config['default_diverging_colormap']
+
+
+from .pyplate import Substance, Container, Plate, Recipe
\ No newline at end of file
diff --git a/pyplate/experiment_design.py b/pyplate/experiment_design.py
new file mode 100644
index 0000000..e7984f5
--- /dev/null
+++ b/pyplate/experiment_design.py
@@ -0,0 +1,318 @@
+from typing import Union, Optional, List
+
+from pyplate import Substance, Container
+import itertools
+
+
+class Factor:
+    """"
+    A factor is a variable that is under the control of the experimenter.
+    It has a name (must be unique to a given experimental space), a list
+    of possible values, and a verifier function.
+    verifier takes a well object and a desired value and returns True if
+    the well is consistent with the description.
+
+    :param name: The name of the factor
+    :type name: str
+    :param possible_values: A list of possible values for the factor
+    :type possible_values: list[str | Substance]
+    :param verifier: A function that takes a well and a value and returns True if the well is consistent with the value
+    :type verifier: callable
+    """
+
+    name: str
+    possible_values: list[str | Substance]
+    verifier: callable
+
+    def __init__(self, name: str, possible_values: list[str | Substance], verifier: callable):
+        self.name = name
+        self.possible_values = possible_values
+        self.verifier = verifier
+
+    def __repr__(self):
+        return f"Factor({self.name}, {self.possible_values}, {self.verifier})"
+
+    def __str__(self):
+        return f"Factor: {self.name} with possible values {self.possible_values}"
+
+    def __eq__(self, other):
+        return (self.name == other.name
+                and self.possible_values == other.possible_values
+                and self.verifier == other.verifier)
+
+    def __hash__(self):
+        return hash((self.name, tuple(self.possible_values), self.verifier))
+
+
+class Experiment:
+    """
+    An experiment represents a single experiment within an experimental space. It has
+    a dictionary of factors and their values. It's experiment_id is a unique identifier
+    for the experiment within the experimental space. The replicate_idx distinguishes
+    between multiple replicates of the same experiment. The well object is a reference
+    to the well in which the experiment was performed.
+
+    :param factors: A dictionary of factors and their values
+    :type factors: dict
+    :param experiment_id: A unique identifier for the experiment within the experimental space
+    :type experiment_id: int
+    :param replicate_idx: A unique identifier for the replicate of the experiment
+    :type replicate_idx: int
+    :param well: A reference to the well in which the experiment was performed
+    :type well: Container
+    """
+
+    factors: dict[str, [str | Substance]]
+    experiment_id: int
+    replicate_idx: int
+    well: Optional[Container]
+
+    def __init__(self, factors: dict[str, [str | Substance]], experiment_id: int, replicate_idx: int,
+                 well: Optional[Container] = None):
+        self.factors = factors
+        self.experiment_id = experiment_id
+        self.replicate_idx = replicate_idx
+        self.well = well
+
+    def map_well(self, well: Container) -> None:
+        """
+        Map the experiment to a well. This is useful for keeping track of
+        which well an experiment was performed in.
+
+        :param well: The well to map to
+        :type well: Container
+        :return: None
+        """
+        self.well = well
+
+    def __repr__(self):
+        return f"Experiment({self.factors}, {self.experiment_id}, {self.replicate_idx}, {self.well})"
+
+    def __str__(self):
+        return (f"Experiment: {self.factors} "
+                f"with experiment_id {self.experiment_id} "
+                f"and replicate_idx {self.replicate_idx}")
+
+    def __getitem__(self, key):
+        return self.factors[key]
+
+    def __setitem__(self, key, value):
+        self.factors[key] = value
+
+    def __contains__(self, key):
+        return key in self.factors
+
+    def __iter__(self):
+        return iter(self.factors)
+
+    def __len__(self):
+        return len(self.factors)
+
+    def __eq__(self, other):
+        return (self.factors == other.factors
+                and self.experiment_id == other.experiment_id
+                and self.replicate_idx == other.replicate_idx
+                and self.well == other.well)
+
+    def __hash__(self):
+        return hash((frozenset(self.factors.items()), self.experiment_id, self.replicate_idx, self.well))
+
+
+class ExperimentalSpace:
+    """
+    An experimental space is a collection of experiments that share the same factors.
+    It a list of factors, a list of blocks, and a list of all experiments. The blocks
+    allow for grouping of experiments with similar fixed factors.
+
+    :param factors: A set of factors that define the experimental space
+    :type factors: set[Factor]
+    :param experiment_id_generator: A function that generates unique experiment ids
+    :type experiment_id_generator: callable
+    """
+
+    factors: set[Factor]
+    blocks: list[Experiment]
+    experiments: list[Experiment]
+    blocks: Optional[list[Experiment]]
+    experiments: Optional[list[Experiment]]
+
+    def __init__(self, factors: set[Factor], experiment_id_generator: callable):
+        self.factors = factors
+        self.experiment_id_generator = experiment_id_generator
+        self.blocks = None
+        self.experiments = None
+
+    def register_factor(self, factor: Factor) -> None:
+        """
+        Register a factor with the experimental space. Ensures that factor
+        names are unique.
+
+        :param factor: The factor to register
+        :type factor: Factor
+        :return: None
+        """
+        if factor in self.factors:
+            raise ValueError(f"Factor {factor.name} already exists in experimental space")
+        self.factors.add(factor)
+
+    def add_experiment(self, experiment: Experiment) -> None:
+        """
+        Add an experiment to the experimental space. This is useful for manually
+        adding experiments without `generate_experiments`.
+
+        :param experiment: The experiment to add
+        :type experiment: Experiment
+        :return: None
+        """
+        if not self.experiments:
+            self.experiments = []
+        if experiment.factors.keys() != {x.name for x in self.factors}:
+            raise ValueError(f"Experiment factors {experiment.factors.keys()}"
+                             f" do not match experimental space factors {self.factors}")
+        for factor_name in experiment.factors:
+            factor_object = self.get_registered_factor(factor_name)
+            if experiment[factor_name] not in factor_object.possible_values:
+                raise ValueError(f"Experiment factor {factor_name} value {experiment[factor_name]}"
+                                 f" not in possible values {factor_object.possible_values}")
+        self.experiments.append(experiment)
+
+    def filter_experiments(self, filter_function: callable) -> None:
+        """
+        Filter the experiments in the experimental space. This is useful for
+        removing experiments that do not meet certain criteria.
+
+        :param filter_function: A function that takes an experiment and returns True if it should be kept
+        :type filter_function: callable
+        :return: None
+        """
+        self.experiments = [x for x in self.experiments if filter_function(x)]
+
+    def get_registered_factor(self, factor_name) -> Factor:
+        """
+        Get a registered factor by name.
+
+        :param factor_name: The name of the factor to get
+        :type factor_name: str
+        :return: The factor object
+        :rtype: Factor
+        """
+        for factor in self.factors:
+            if factor.name == factor_name:
+                return factor
+        raise ValueError(f"Factor {factor_name} not found in experimental space")
+
+    def generate_experiments(self, factors: dict[str, [str | Substance]],
+                             n_replicates: int, blocking_factors: list[str]) -> list[list[Experiment]]:
+        """
+        Generate experiments for the experimental space. This is useful for
+        generating a full factorial design of a subset of the experimental space.
+
+        :param factors: A list of factors to generate experiments for (must be a subset of the experimental space factors)
+        :type factors: dict[str, [str | Substance]]
+        :param n_replicates: The number of replicates to generate for each experiment
+        :type n_replicates: int
+        :param blocking_factors: A list of factors to use for blocking
+        :type blocking_factors: list[str]
+        :return: A list of blocks, where each block is a list of experiments
+        """
+        filtered_factors = {}
+        for factor in self.factors:
+            if factors[factor.name] == "all":
+                filtered_factors[factor.name] = factor.possible_values
+            else:
+                filtered_factors[factor.name] = [value for value in factor.possible_values if value in factors[factor.name]]
+        
+        # Initialize blocks
+        blocks = {}
+
+        """
+        The following lines return a list of tuples shown below, assuming you block on Factor 1 and Factor 2:
+        [
+            [
+                [("Factor 1", "value 1"), ("Factor 2", "value 1"), ("Factor 3", "value 1")],
+                [("Factor 1", "value 1"), ("Factor 2", "value 1"), ("Factor 3", "value 2")],
+            ],
+            [
+                [("Factor 1", "value 1"), ("Factor 2", "value 2"), ("Factor 3", "value 1")],
+                [("Factor 1", "value 1"), ("Factor 2", "value 2"), ("Factor 3", "value 2")],
+            ],
+            ...
+        ]
+        
+        Easy to use this to create dicts and thus experiments
+        """
+
+        # If needed form blocking combinations
+        blocking_combinations = list(itertools.product(*[filtered_factors[name] for name in blocking_factors]))
+
+        # Iterate over each combination of blocking factors
+        for block_combination in blocking_combinations:
+            block = []
+            
+            # Filter out the current blocking factors and create combinations of other factors
+            non_blocking_factors = {name: values for name, values in filtered_factors.items() if name not in blocking_factors}
+            other_combinations = list(itertools.product(*non_blocking_factors.values()))
+
+            # Generate experiments for each combination
+            for comb in other_combinations:
+                factors_dict = dict(zip(non_blocking_factors.keys(), comb))
+                
+                # Set the values of the blocking factors
+                for bf, value in zip(blocking_factors, block_combination):
+                    factors_dict[bf] = value
+
+                # Create replicates for each unique combination
+                for rep in range(n_replicates):
+                    experiment = Experiment(factors=factors_dict, replicate_idx=rep+1, experiment_id=self.experiment_id_generator())
+                    block.append(experiment)
+
+            # Use the blocking factor combination as the key for the blocks dictionary
+            blocks[block_combination] = block
+
+        return blocks
+    
+        #Another implementation
+        # experiment_factor_values = {}
+        # for factor in self.factors:
+        #     if factors[factor.name] == "all":
+        #         experiment_factor_values[factor.name] = factor.possible_values
+        #     else:
+        #         experiment_factor_values[factor.name] = [value for value in factor.possible_values if
+        #                                                  value in factors[factor.name]]
+
+        # non_blocking_factor_values = dict(filter(lambda x: x[0] not in blocking_factors, experiment_factor_values.items()))
+        # blocking_factor_values = dict(filter(lambda x: x[0] in blocking_factors, experiment_factor_values.items()))
+
+        # # Initialize blocks
+        # blocks = []
+
+        """
+        The following lines return a list of tuples shown below, assuming you block on Factor 1 and Factor 2:
+        [
+            [
+                [("Factor 1", "value 1"), ("Factor 2", "value 1"), ("Factor 3", "value 1")],
+                [("Factor 1", "value 1"), ("Factor 2", "value 1"), ("Factor 3", "value 2")],
+            ],
+            [
+                [("Factor 1", "value 1"), ("Factor 2", "value 2"), ("Factor 3", "value 1")],
+                [("Factor 1", "value 1"), ("Factor 2", "value 2"), ("Factor 3", "value 2")],
+            ],
+            ...
+        ]
+        
+        Easy to use this to create dicts and thus experiments
+        """
+        # blocking_factor_combos = [list(zip(blocking_factor_values.keys(), values)) for values in itertools.product(*blocking_factor_values.values())]
+        # non_blocking_factor_combos = [list(zip(non_blocking_factor_values.keys(), values)) for values in itertools.product(*non_blocking_factor_values.values())]
+
+        # for blocking_combo in blocking_factor_combos:
+        #     block = []
+        #     for non_blocking_combo in non_blocking_factor_combos:
+        #         for rep in range(n_replicates):
+        #             experiment = Experiment(dict(blocking_combo + non_blocking_combo), self.experiment_id_generator(), rep + 1)
+        #             block.append(experiment)
+        #     blocks.append(block)
+
+        # return blocks
+    
+
diff --git a/pyplate/pyplate.py b/pyplate/pyplate.py
index 2386014..d31455e 100644
--- a/pyplate/pyplate.py
+++ b/pyplate/pyplate.py
@@ -30,12 +30,17 @@
 
 # Allow typing reference while still building classes
 from __future__ import annotations
-import re
+
+from functools import cache
 from typing import Tuple, Dict, Iterable
 from copy import deepcopy, copy
 import numpy
+import pandas
+
 from pyplate.slicer import Slicer
-from pyplate import config
+from . import Config
+
+config = Config()
 
 
 class Unit:
@@ -61,7 +66,8 @@ def convert_prefix_to_multiplier(prefix: str) -> float:
         """
         if not isinstance(prefix, str):
             raise TypeError("SI prefix must be a string.")
-        prefixes = {'n': 1e-9, 'u': 1e-6, 'µ': 1e-6, 'm': 1e-3, 'c': 1e-2, 'd': 1e-1, '': 1, 'k': 1e3, 'M': 1e6}
+        prefixes = {'n': 1e-9, 'u': 1e-6, 'µ': 1e-6, 'm': 1e-3, 'c': 1e-2, 'd': 1e-1, '': 1, 'da': 1e1, 'k': 1e3,
+                    'M': 1e6}
         if prefix in prefixes:
             return prefixes[prefix]
         raise ValueError(f"Invalid prefix: {prefix}")
@@ -81,13 +87,17 @@ def parse_quantity(quantity: str) -> Tuple[float, str]:
 
         """
         if not isinstance(quantity, str):
-            raise TypeError("Amount must be a string.")
-        # (floating point number in 1.0e1 format) possibly some white space (alpha string)
-        match = re.fullmatch(r"([_\d]+(?:\.\d+)?(?:e-?\d+)?)\s*([a-zA-Z]+)", quantity)
-        if not match:
-            raise ValueError("Invalid quantity. Quantity should be in the format '10 mL'.")
-        value, unit = match.groups()
-        value = float(value)
+            raise TypeError("Quantity must be a string.")
+
+        if quantity.count(' ') != 1:
+            raise ValueError("Value and unit must be separated by a single space.")
+
+        value, unit = quantity.split(' ')
+        try:
+            value = float(value)
+        except ValueError as exc:
+            raise ValueError("Value is not a valid float.") from exc
+
         if unit == 'U':
             return value, unit
         for base_unit in ['mol', 'g', 'L', 'M']:
@@ -98,84 +108,155 @@ def parse_quantity(quantity: str) -> Tuple[float, str]:
         raise ValueError("Invalid unit {base_unit}.")
 
     @staticmethod
-    def convert_to_unit_value(substance: Substance, quantity: str, volume: float = 0.0) -> float:
+    def parse_concentration(concentration):
         """
+        Parses concentration string to (value, numerator, denominator).
+        Args:
+            concentration: concentration, '1 M', '1 umol/uL', '0.1 umol/10 uL'
 
-        Converts amount to standard unit.
+        Returns: Tuple of value, numerator, denominator. (0.01, 'mol', 'L')
 
+        """
+        if '/' not in concentration:
+            if concentration[-1] == 'm':
+                concentration = concentration[:-1] + 'mol/kg'
+            elif concentration[-1] == 'M':
+                concentration = concentration[:-1] + 'mol/L'
+            else:
+                raise ValueError("Only m and M are allowed as concentration units.")
+        replacements = {'%v/v': 'L/L', '%w/w': 'g/g', '%w/v': config.default_weight_volume_units}
+        if concentration[-4:] in replacements:
+            concentration = concentration[:-4] + replacements[concentration[-4:]]
+            numerator, denominator = map(str.split, concentration.split('/'))
+            numerator[0] = float(numerator[0]) / 100  # percent
+        else:
+            numerator, denominator = map(str.split, concentration.split('/'))
+        if len(numerator) < 2 or len(denominator) < 1:
+            raise ValueError("Concentration must be of the form '1 umol/mL'.")
+        try:
+            numerator[0] = float(numerator[0])
+            if len(denominator) > 1:
+                numerator[0] /= float(denominator.pop(0))
+        except ValueError as exc:
+            raise ValueError("Value is not a float.") from exc
+        units = ('mol', 'L', 'g')
+        for unit in units:
+            if numerator[1].endswith(unit):
+                numerator[0] *= Unit.convert_prefix_to_multiplier(numerator[1][:-len(unit)])
+                numerator[1] = unit
+            if denominator[0].endswith(unit):
+                numerator[0] /= Unit.convert_prefix_to_multiplier(denominator[0][:-len(unit)])
+                denominator[0] = unit
+        if numerator[1] not in ('U', 'mol', 'L', 'g') or denominator[0] not in ('mol', 'L', 'g'):
+            raise ValueError("Concentration must be of the form '1 umol/mL'.")
+        return round(numerator[0], config.internal_precision), numerator[1], denominator[0]
 
-        Arguments:
-            substance: Substance in question.
-            quantity: Quantity of substance.
-            volume: Volume of containing Mixture in mL
+    @staticmethod
+    def convert_from(substance: Substance, quantity: float, from_unit: str, to_unit: str) -> float:
+        """
+                    Convert quantity of substance between units.
 
-        Returns: Value in standard unit, converted to storage format.
+                    Arguments:
+                        substance: Substance in question.
+                        quantity: Quantity of substance.
+                        from_unit: Unit to convert quantity from ('mL').
+                        to_unit: Unit to convert quantity to ('mol').
 
-        ---
+                    Returns: Converted value.
 
-        Standard units:
+                """
 
-        - SOLID: moles
+        if not isinstance(substance, Substance):
+            raise TypeError(f"Invalid type for substance, {type(substance)}")
+        if not isinstance(quantity, (int, float)):
+            raise TypeError("Quantity must be a float.")
+        if not isinstance(from_unit, str) or not isinstance(to_unit, str):
+            raise TypeError("Unit must be a str.")
+
+        if from_unit == 'U' and not substance.is_enzyme():
+            raise ValueError("Only enzymes can be measured in activity units. 'U'")
 
-        - LIQUID: liter
+        for suffix in ['U', 'L', 'g', 'mol']:
+            if from_unit.endswith(suffix):
+                prefix = from_unit[:-len(suffix)]
+                quantity *= Unit.convert_prefix_to_multiplier(prefix)
+                from_unit = suffix
+                break
 
-        - ENZYME: activity unit
+        result = None
+        if substance.is_enzyme() and not to_unit.endswith('U'):
+            return 0
+        if to_unit.endswith('U'):
+            prefix = to_unit[:-1]
+            if not substance.is_enzyme():
+                return 0
+            if not from_unit == 'U':
+                raise ValueError("Enzymes can only be measured in activity units. 'U'")
+            result = quantity
+        elif to_unit.endswith('L'):
+            prefix = to_unit[:-1]
+            if from_unit == 'L':
+                result = quantity
+            elif from_unit == 'mol':
+                # mol * g/mol / (g/mL)
+                result_in_mL = quantity * substance.mol_weight / substance.density
+                result = result_in_mL / 1000
+            elif from_unit == 'g':
+                # g / (g/mL)
+                result_in_mL = quantity / substance.density
+                result = result_in_mL / 1000
+        elif to_unit.endswith('mol'):
+            prefix = to_unit[:-3]
+            if from_unit == 'L':
+                value_in_mL = quantity * 1000
+                # mL * g/mL / (g/mol)
+                result = value_in_mL * substance.density / substance.mol_weight
+            elif from_unit == 'mol':
+                result = quantity
+            elif from_unit == 'g':
+                # g / (g/mol)
+                result = quantity / substance.mol_weight
+        elif to_unit.endswith('g'):
+            prefix = to_unit[:-1]
+            if from_unit == 'L':
+                value_in_mL = quantity * 1000
+                # mL * g/mL
+                result = value_in_mL * substance.density
+            elif from_unit == 'mol':
+                # mol * g/mol
+                result = quantity * substance.mol_weight
+            elif from_unit == 'g':
+                result = quantity
+        else:
+            raise ValueError("Only L, U, g, and mol are valid units.")
 
+        assert result is not None, f"{substance} {quantity} {from_unit} {to_unit}"
+
+        return result / Unit.convert_prefix_to_multiplier(prefix)
+
+    @staticmethod
+    def convert(substance: Substance, quantity: str, unit: str) -> float:
         """
+            Convert quantity of substance to unit.
 
-        #          +--------+--------+--------+--------+--------+----------+
-        #          |              Valid Input                   | Standard |
-        #          |    g   |    L   |   mol  |    M   |    U   |   Unit   |
-        # +--------+--------+--------+--------+--------+--------+----------+
-        # |SOLID   |   Yes  |   No   |   Yes  |   Yes  |   No   |    mol   |
-        # +--------+--------+--------+--------+--------+--------+----------+
-        # |LIQUID  |   Yes  |   Yes  |   Yes  |   No   |   No   |    mL    |
-        # +--------+--------+--------+--------+--------+--------+----------+
-        # |ENZYME  |   No   |   No   |   No   |   No   |   Yes  |    U     |
-        # +--------+--------+--------+--------+--------+--------+----------+
+            Arguments:
+                substance: Substance in question.
+                quantity: Quantity of substance ('10 mL').
+                unit: Unit to convert quantity to ('mol').
 
-        # TODO: Explain why moles and mL were chosen as units.
-        # TODO: Enzyme could be in solution.
+            Returns: Converted value.
+
+        """
 
         if not isinstance(substance, Substance):
             raise TypeError(f"Invalid type for substance, {type(substance)}")
         if not isinstance(quantity, str):
             raise TypeError("Quantity must be a str.")
-        if not isinstance(volume, float):
-            raise TypeError("Volume, if provided, must be a float.")
-
-        value, unit = Unit.parse_quantity(quantity)
-        if substance.is_solid():  # Convert to moles
-            if unit == 'g':  # mass
-                result = value / substance.mol_weight
-            elif unit == 'mol':  # moles
-                result = value
-            elif unit == 'M':  # molar
-                # A molar concentration with zero volume would be undefined.
-                if volume <= 0.0:
-                    raise ValueError('Must have a liquid in which to dissolve the solid ' +
-                                     'in order to create a molar concentration')
-                # value = molar concentration in mol/L, volume = volume in mL
-                result = value * volume / 1000
-            else:
-                raise ValueError("We only measure solids in grams and moles.")
-            return Unit.convert_to_storage(result, 'mol')
-        if substance.is_liquid():
-            if unit == 'g':  # mass
-                # g -> mL
-                result = value / substance.density
-            elif unit == 'L':  # volume
-                result = value * 1000
-            elif unit == 'mol':  # moles
-                # mol -> mL
-                result = value / substance.concentration
-            else:
-                raise ValueError
-            return Unit.convert_to_storage(result, 'mL')
-        if substance.is_enzyme():
-            if unit == 'U':
-                return value
-            raise ValueError
+        if not isinstance(unit, str):
+            raise TypeError("Unit must be a str.")
+
+        value, quantity_unit = Unit.parse_quantity(quantity)
+        return Unit.convert_from(substance, value, quantity_unit, unit)
 
     @staticmethod
     def convert_to_storage(value: float, unit: str) -> float:
@@ -198,10 +279,10 @@ def convert_to_storage(value: float, unit: str) -> float:
 
         if unit[-1] == 'L':
             prefix_value = Unit.convert_prefix_to_multiplier(unit[:-1])
-            result = value * prefix_value / config.volume_storage
+            result = value * prefix_value / Unit.convert_prefix_to_multiplier(config.volume_prefix[0])
         else:  # moles
             prefix_value = Unit.convert_prefix_to_multiplier(unit[:-3])
-            result = value * prefix_value / config.moles_storage
+            result = value * prefix_value / Unit.convert_prefix_to_multiplier(config.moles_prefix[0])
         return round(result, config.internal_precision)
 
     @staticmethod
@@ -218,24 +299,161 @@ def convert_from_storage(value: float, unit: str) -> float:
         Returns: Converted value.
 
         """
-        if not isinstance(value, float):
+        if not isinstance(value, (int, float)):
             raise TypeError("Value must be a float.")
         if not isinstance(unit, str):
             raise TypeError("Unit must be a str.")
 
         if unit[-1] == 'L':
             prefix_value = Unit.convert_prefix_to_multiplier(unit[:-1])
-            result = value * config.volume_storage / prefix_value
+            result = value * Unit.convert_prefix_to_multiplier(config.volume_prefix[0]) / prefix_value
         else:  # moles
             prefix_value = Unit.convert_prefix_to_multiplier(unit[:-3])
-            result = value * config.moles_storage / prefix_value
+            result = value * Unit.convert_prefix_to_multiplier(config.moles_prefix[0]) / prefix_value
         return round(result, config.internal_precision)
 
+    @staticmethod
+    def convert_from_storage_to_standard_format(what, quantity: float):
+        """
+        Converts a quantity of a substance or container to a standard format.
+        Example: (water, 1e6) -> (18.015, 'mL'), (NaCl, 1e6) -> (58.443, 'g')
+        Args:
+            what: Substance or Container
+            quantity: Quantity in storage format.
+
+        Returns: Tuple of quantity and unit.
+
+        """
+        if isinstance(what, Substance):
+            if what.is_enzyme():
+                unit = 'U'
+            elif what.is_solid():
+                unit = 'g'
+                # convert moles to grams
+                # molecular weight is in g/mol
+                quantity *= Unit.convert_prefix_to_multiplier(config.moles_prefix[0]) * what.mol_weight
+            elif what.is_liquid():
+                unit = 'L'
+                # convert moles to liters
+                # molecular weight is in g/mol
+                # density is in g/mL
+                quantity *= (Unit.convert_prefix_to_multiplier(config.moles_prefix[0])
+                             * what.mol_weight / what.density / 1e3)
+        elif isinstance(what, Container):
+            # Assume the container contains a liquid
+            unit = 'L'
+            quantity *= Unit.convert_prefix_to_multiplier(config.volume_prefix[0])
+        else:
+            raise TypeError("Invalid type for what.")
+
+        multiplier = 1
+        while quantity < 1:
+            quantity *= 1e3
+            multiplier /= 1e3
+        quantity = round(quantity, config.external_precision)
+        unit = {1: '', 1e-3: 'm', 1e-6: 'u'}[multiplier] + unit
+        return quantity, unit
+
+    @staticmethod
+    def get_human_readable_unit(value: float, unit: str):
+        """
+        Returns a more human-readable value and unit.
+
+        Args:
+            value: Value to work with.
+            unit:  Unit to determine type and default unit if value is zero.
+
+        Returns: Tuple of new value and unit
+
+        """
+        if value == 0:
+            return value, unit
+        value = abs(value)
+        if unit[-1] == 'L':
+            unit = 'L'
+        elif unit[-3:] == 'mol':
+            unit = 'mol'
+        elif unit[-1] == 'g':
+            unit = 'g'
+        elif unit[-1] == 'U':
+            unit = 'U'
+        multiplier = 1.0
+        while value < 1:
+            value *= 1e3
+            multiplier /= 1e3
+
+        multiplier = max(multiplier, 1e-6)
+
+        return value, {1: '', 1e-3: 'm', 1e-6: 'u'}[multiplier] + unit
+
+    @staticmethod
+    def calculate_concentration_ratio(solute: Substance, concentration: str, solvent: Substance):
+        """
+        Helper function for dealing with concentrations.
+
+        Returns: ratio of moles or Activity Units per mole storage unit ('umol', etc.).
+
+        """
+        c, numerator, denominator = Unit.parse_concentration(concentration)
+        if numerator not in ('g', 'L', 'mol', 'U'):
+            raise ValueError("Invalid unit in numerator.")
+        if denominator not in ('g', 'L', 'mol'):
+            raise ValueError("Invalid unit in denominator.")
+
+        ratio = None  # ration of solute to solvent in moles
+        if numerator == 'g':
+            if denominator == 'g':
+                ratio = c * solvent.mol_weight / (1 - c) / solute.mol_weight
+            elif denominator == 'mol':
+                ratio = c / (solute.mol_weight - c)
+            elif denominator == 'L':
+                c /= 1000  # g/mL
+                ratio = c * solvent.mol_weight / (solute.mol_weight * solvent.density * (1 - c / solute.density))
+        elif numerator == 'L':
+            if denominator == 'g':
+                c *= 1000  # mL/g
+                ratio = c * solvent.mol_weight / (solute.mol_weight * (1 / solute.density - c))
+            elif denominator == 'mol':
+                c *= 1000  # mL/mol
+                ratio = c / (solute.mol_weight / solute.density - c)
+            elif denominator == 'L':
+                ratio = c * solvent.mol_weight / solvent.density / (solute.mol_weight / solute.density) / (1 - c)
+        elif numerator == 'mol':
+            if denominator == 'g':
+                ratio = c * solvent.mol_weight / (1 - c * solute.mol_weight)
+            elif denominator == 'mol':
+                ratio = c / (1 - c)
+            elif denominator == 'L':
+                c /= 1000  # mol/mL
+                ratio = c * solvent.mol_weight / solvent.density / (1 - c * solute.mol_weight / solute.density)
+        elif numerator == 'U':
+            if denominator == 'g':
+                ratio = c * solvent.mol_weight
+            elif denominator == 'L':
+                # density is g/mL
+                ratio = c * solvent.mol_weight / solvent.density / 1000.0
+            else:
+                ratio = c
+            # ratio can be multiplied by a stored value of moles to get number of U
+            ratio *= Unit.convert_from_storage(1, 'mol')
+        return ratio, numerator, denominator
+
+    @staticmethod
+    def calculate_concentration_ratio_moles(solute: Substance, quantity: str, solvent: Substance):
+        q, q_unit = Unit.parse_quantity(quantity)
+        if q_unit not in ('g', 'L', 'mol'):
+            raise ValueError("Invalid unit in quantity.")
+        if q_unit == 'g':
+            q /= solute.mol_weight
+        elif q_unit == 'L':
+            q *= solute.density
+        return Unit.calculate_concentration_ratio(solute, f"{q} mol", solvent)
+
 
 class Substance:
     """
     An abstract chemical or biological entity (e.g., reagent, enzyme, solvent, etc.). Immutable.
-    Solids and enzymes are assumed to require zero volume.
+    Enzymes are assumed to require zero volume.
 
     Attributes:
         name: Name of substance.
@@ -244,6 +462,7 @@ class Substance:
         concentration: Calculated concentration if `Substance` is a liquid (mol/mL).
         molecule: `cctk.Molecule` if provided.
     """
+    # An attribute `classes` could be added to support classes of substances.
     SOLID = 1
     LIQUID = 2
     ENZYME = 3
@@ -270,15 +489,12 @@ def __init__(self, name: str, mol_type: int, molecule=None):
 
         self.name = name
         self._type = mol_type
-        self.mol_weight = self.density = self.concentration = None
+        self.mol_weight = self.concentration = None
+        self.density = float('inf')
         self.molecule = molecule
 
     def __repr__(self):
-        if self.is_solid():
-            return f"SOLID ({self.name}: {self.mol_weight})"
-        if self.is_liquid():
-            return f"LIQUID ({self.name}: {self.mol_weight}, {self.density})"
-        return f"ENZYME ({self.name})"
+        return f"{self.name} ({'SOLID' if self.is_solid() else 'LIQUID' if self.is_liquid() else 'ENZYME'})"
 
     def __eq__(self, other):
         if not isinstance(other, Substance):
@@ -312,6 +528,7 @@ def solid(name: str, mol_weight: float, molecule=None) -> Substance:
 
         substance = Substance(name, Substance.SOLID, molecule)
         substance.mol_weight = mol_weight
+        substance.density = config.solid_density
         return substance
 
     @staticmethod
@@ -393,7 +610,7 @@ class Container:
         max_volume: Maximum volume Container can hold in storage format.
     """
 
-    def __init__(self, name: str, max_volume: float = float('inf'),
+    def __init__(self, name: str, max_volume: str = 'inf L',
                  initial_contents: Iterable[Tuple[Substance, str]] = None):
         """
         Create a Container.
@@ -403,20 +620,20 @@ def __init__(self, name: str, max_volume: float = float('inf'),
             max_volume: Maximum volume that can be stored in the container in mL
             initial_contents: (optional) Iterable of tuples of the form (Substance, quantity)
         """
-        # TODO: make max_volume a str
         if not isinstance(name, str):
             raise TypeError("Name must be a str.")
         if len(name) == 0:
             raise ValueError("Name must not be empty.")
 
-        if not isinstance(max_volume, (int, float)):
-            raise TypeError("Maximum volume must be a float.")
-        if not max_volume > 0:
+        if not isinstance(max_volume, str):
+            raise TypeError("Maximum volume must be a str, ('10 mL').")
+        max_volume, _ = Unit.parse_quantity(max_volume)
+        if max_volume <= 0:
             raise ValueError("Maximum volume must be positive.")
         self.name = name
         self.contents: Dict[Substance, float] = {}
         self.volume = 0.0
-        self.max_volume = Unit.convert_to_storage(max_volume, 'mL')
+        self.max_volume = Unit.convert_to_storage(max_volume, 'L')
         if initial_contents:
             if not isinstance(initial_contents, Iterable):
                 raise TypeError("Initial contents must be iterable.")
@@ -427,6 +644,22 @@ def __init__(self, name: str, max_volume: float = float('inf'),
                 if not isinstance(substance, Substance) or not isinstance(quantity, str):
                     raise TypeError("Element in initial_contents must be a (Substance, str) tuple.")
                 self._self_add(substance, quantity)
+            contents = []
+            for substance, quantity in self.contents.items():
+                quantity, unit = Unit.convert_from_storage_to_standard_format(substance, quantity)
+                contents.append(f"{quantity} {unit} of {substance.name}")
+            self.instructions = f"Add {', '.join(contents)} "
+            if self.max_volume != float('inf'):
+                max_volume, unit = Unit.convert_from_storage_to_standard_format(self, self.max_volume)
+                self.instructions += f" to a {max_volume} {unit} container."
+            else:
+                self.instructions += " to a container."
+        else:
+            if self.max_volume != float('inf'):
+                max_volume, unit = Unit.convert_from_storage_to_standard_format(self, self.max_volume)
+                self.instructions = f"Create a {max_volume} {unit} container."
+            else:
+                self.instructions = "Create a container."
 
     def __eq__(self, other):
         if not isinstance(other, Container):
@@ -452,62 +685,93 @@ def _self_add(self, source: Substance, quantity: str) -> None:
             raise TypeError("Source must be a Substance.")
         if not isinstance(quantity, str):
             raise TypeError("Quantity must be a str.")
-        if quantity.endswith('M') and not source.is_solid():
-            # TODO: molarity from liquids?
-            raise ValueError("Molarity solutions can only be made from solids.")
-        volume = Unit.convert_from_storage(self.volume, 'mL')
-        amount_to_transfer = round(Unit.convert_to_unit_value(source, quantity, volume), config.internal_precision)
-        # add source Substance to self.contents
-        self.contents[source] = round(self.contents.get(source, 0) + amount_to_transfer, config.internal_precision)
-        if source.is_liquid():
-            self.volume = round(self.volume + amount_to_transfer, config.internal_precision)
-        if self.volume > self.max_volume:
+
+        if source.is_enzyme():
+            volume_to_add = 0
+            amount_to_add = Unit.convert(source, quantity, 'U')
+        else:
+            volume_to_add = Unit.convert(source, quantity, config.volume_prefix)
+            amount_to_add = Unit.convert(source, quantity, config.moles_prefix)
+        if self.volume + volume_to_add > self.max_volume:
             raise ValueError("Exceeded maximum volume")
+        self.volume = round(self.volume + volume_to_add, config.internal_precision)
+        self.contents[source] = round(self.contents.get(source, 0) + amount_to_add, config.internal_precision)
 
-    def _transfer(self, source_container: Container, volume: str):
+    def _transfer(self, source_container: Container, quantity: str):
         """
-        Move volume ('10 mL') from container to self.
+        Move quantity ('10 mL', '5 mg') from container to self.
 
         Arguments:
             source_container: `Container` to transfer from.
-            volume: How much to transfer.
+            quantity: How much to transfer.
 
         Returns: New source and destination container.
         """
 
         if not isinstance(source_container, Container):
             raise TypeError("Invalid source type.")
-        volume_to_transfer, unit = Unit.parse_quantity(volume)
-        volume_to_transfer = Unit.convert_to_storage(volume_to_transfer, 'L')
-        volume_to_transfer = round(volume_to_transfer, config.internal_precision)
-        if unit != 'L':
-            raise ValueError("We can only transfer volumes from other containers.")
-        if volume_to_transfer > source_container.volume:
-            raise ValueError("Not enough mixture left in source container." +
-                             f" {Unit.convert_from_storage(volume_to_transfer, 'mL')} mL "
-                             f" needed of {source_container.name}" +
-                             f" out of {Unit.convert_from_storage(source_container.volume, 'mL')} mL.")
-        # source_container, to = source_container.copy(), self.copy()
+        quantity_to_transfer, unit = Unit.parse_quantity(quantity)
+
+        if unit == 'L':
+            volume_to_transfer = Unit.convert_to_storage(quantity_to_transfer, 'L')
+            volume_to_transfer = round(volume_to_transfer, config.internal_precision)
+
+            if volume_to_transfer > source_container.volume:
+                raise ValueError(f"Not enough mixture left in source container ({source_container.name}). " +
+                                 f"Only {Unit.convert_from_storage(source_container.volume, 'mL')} mL available, " +
+                                 f"{Unit.convert_from_storage(volume_to_transfer, 'mL')} mL needed.")
+            ratio = volume_to_transfer / source_container.volume
+
+        elif unit == 'g':
+            mass_to_transfer = round(quantity_to_transfer, config.internal_precision)
+            total_mass = sum(Unit.convert(substance, f"{amount} {config.moles_prefix}", "g") for
+                             substance, amount in source_container.contents.items())
+            ratio = mass_to_transfer / total_mass
+        elif unit == 'mol':
+            moles_to_transfer = Unit.convert_to_storage(quantity_to_transfer, 'mol')
+            total_moles = sum(amount for substance, amount in source_container.contents.items()
+                              if not substance.is_enzyme())
+            ratio = moles_to_transfer / total_moles
+        else:
+            raise ValueError("Invalid quantity unit.")
+
         source_container, to = deepcopy(source_container), deepcopy(self)
-        ratio = volume_to_transfer / source_container.volume
         for substance, amount in source_container.contents.items():
-            to.contents[substance] = round(to.contents.get(substance, 0) + amount * ratio,
+            to_transfer = amount * ratio
+            to.contents[substance] = round(to.contents.get(substance, 0) + to_transfer,
                                            config.internal_precision)
-            source_container.contents[substance] = round(source_container.contents[substance] - amount * ratio,
+            source_container.contents[substance] = round(source_container.contents[substance] - to_transfer,
                                                          config.internal_precision)
-        to.volume = round(to.volume + volume_to_transfer, config.internal_precision)
+            # if quantity to remove is the same as the current amount plus a very small delta,
+            # we will get a negative 0 answer.
+            if source_container.contents[substance] == -0.0:
+                source_container.contents[substance] = 0.0
+        if source_container.has_liquid():
+            transfer = Unit.convert_from_storage(ratio * source_container.volume, 'L')
+            transfer, unit = Unit.get_human_readable_unit(transfer, 'L')
+        else:
+            # total mass in source container times ratio
+            mass = sum(Unit.convert(substance, f"{amount} {config.moles_prefix if not substance.is_enzyme() else 'U'}",
+                                    "mg") for substance, amount in source_container.contents.items())
+            transfer, unit = Unit.get_human_readable_unit(mass * ratio, 'mg')
+
+        to.instructions += f"\nTransfer {transfer} {unit} of {source_container.name} to {to.name}"
+        to.volume = round(sum(Unit.convert(substance, f"{amount} {config.moles_prefix}", config.volume_prefix) for
+                              substance, amount in to.contents.items()), config.internal_precision)
         if to.volume > to.max_volume:
             raise ValueError(f"Exceeded maximum volume in {to.name}.")
-        source_container.volume = round(source_container.volume - volume_to_transfer, config.internal_precision)
+        source_container.volume = sum(Unit.convert(substance, f"{amount} {config.moles_prefix}", config.volume_prefix)
+                                      for substance, amount in source_container.contents.items())
+        source_container.volume = round(source_container.volume, config.internal_precision)
         return source_container, to
 
-    def _transfer_slice(self, source_slice: Plate or PlateSlicer, volume: str):
+    def _transfer_slice(self, source_slice: Plate or PlateSlicer, quantity: str):
         """
-        Move volume ('10 mL') from slice to self.
+        Move quantity ('10 mL', '5 mg') from each well in a slice to self.
 
         Arguments:
             source_slice: Slice or Plate to transfer from.
-            volume: How much to transfer.
+            quantity: How much to transfer.
 
         Returns:
             A new plate and a new container, both modified.
@@ -515,7 +779,7 @@ def _transfer_slice(self, source_slice: Plate or PlateSlicer, volume: str):
 
         def helper_func(elem):
             """ Moves volume from elem to to_array[0]"""
-            elem, to_array[0] = Container.transfer(elem, to_array[0], volume)
+            elem, to_array[0] = Container.transfer(elem, to_array[0], quantity)
             return elem
 
         if isinstance(source_slice, Plate):
@@ -525,54 +789,59 @@ def helper_func(elem):
         to = deepcopy(self)
         source_slice = copy(source_slice)
         source_slice.plate = deepcopy(source_slice.plate)
-        volume_to_transfer, unit = Unit.parse_quantity(volume)
-        volume_to_transfer = Unit.convert_to_storage(volume_to_transfer, 'L')
-        volume_to_transfer = round(volume_to_transfer, config.internal_precision)
-        if unit != 'L':
-            raise ValueError("We can only transfer volumes from other containers.")
-        if source_slice.get().size * volume_to_transfer > (to.max_volume - to.volume):
-            raise ValueError(f"Exceeded maximum volume in {to.name}.")
 
         to_array = [to]
-        result = numpy.vectorize(helper_func, cache=True)(source_slice.get())
-        source_slice.set(result)
-        return source_slice.plate, to_array[0]
+        source_slice.apply(helper_func)
+        to = to_array[0]
+        return source_slice.plate, to
 
     def __repr__(self):
+        contents = []
+        for substance, value in sorted(self.contents.items(), key=lambda elem: (elem[0]._type, -elem[1])):
+            if substance.is_enzyme():
+                contents.append(f"{substance}: {value} U")
+            else:
+                value, unit = Unit.get_human_readable_unit(Unit.convert_from_storage(value, 'mol'), 'mmol')
+                contents.append(
+                    f"{substance}: {round(value, config.external_precision)} {unit}")
+
         max_volume = ('/' + str(Unit.convert_from_storage(self.max_volume, 'mL'))) \
             if self.max_volume != float('inf') else ''
         return f"Container ({self.name}) ({Unit.convert_from_storage(self.volume, 'mL')}" + \
-            f"{max_volume} mL of ({self.contents})"
+            f"{max_volume} mL of ({contents})"
 
-    @staticmethod
-    def add(source: Substance, destination: Container, quantity: str) -> Container:
+    @cache
+    def has_liquid(self):
+        """
+        Returns: True if any substance in the container is a liquid.
+        """
+        return any(substance.is_liquid() for substance in self.contents)
+
+    def _add(self, source: Substance, quantity: str) -> Container:
         """
-        Add the given quantity ('10 mol') of the source substance to the destination container.
+        Add the given quantity ('10 mol') of the source substance to the container.
 
         Arguments:
             source: Substance to add to `destination`.
-            destination: Container to add to.
             quantity: How much `Substance` to add.
 
         Returns:
-            A new copy of `destination` container.
+            A new container with added substance.
         """
-        if not isinstance(destination, Container):
-            raise TypeError("You can only use Container.add to add to a Container")
-        destination = deepcopy(destination)
+        destination = deepcopy(self)
         destination._self_add(source, quantity)
         return destination
 
     @staticmethod
-    def transfer(source: Container or Plate or PlateSlicer, destination: Container, volume):
+    def transfer(source: Container | Plate | PlateSlicer, destination: Container, quantity: str):
         """
-        Move volume ('10 mL') from source to destination container,
+        Move quantity ('10 mL', '5 mg') from source to destination container,
         returning copies of the objects with amounts adjusted accordingly.
 
         Arguments:
             source: Container, plate, or slice to transfer from.
             destination: Container to transfer to:
-            volume: How much to transfer.
+            quantity: How much to transfer.
 
         Returns:
             A tuple of (T, Container) where T is the type of the source.
@@ -580,49 +849,377 @@ def transfer(source: Container or Plate or PlateSlicer, destination: Container,
         if not isinstance(destination, Container):
             raise TypeError("You can only use Container.transfer into a Container")
         if isinstance(source, Container):
-            return destination._transfer(source, volume)
+            return destination._transfer(source, quantity)
         if isinstance(source, (Plate, PlateSlicer)):
-            return destination._transfer_slice(source, volume)
+            return destination._transfer_slice(source, quantity)
         raise TypeError("Invalid source type.")
 
+    def get_concentration(self, solute, units='M'):
+        """
+        Get the concentration of solute in the current solution.
+
+        Args:
+            solute: Substance interested in.
+            units: Units to return concentration in, defaults to Molar.
+
+        Returns: Concentration
+
+        """
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(units, str):
+            raise TypeError("Units must be a str.")
+
+        mult, *units = Unit.parse_concentration('1 ' + units)
+
+        numerator = Unit.convert(solute, f"{self.contents.get(solute, 0)} {config.moles_prefix}", units[0])
+        if units[1].endswith('L'):
+            denominator = Unit.convert_from_storage(self.volume, units[1])
+        else:
+            denominator = sum(Unit.convert(substance, f"{amount} {config.moles_prefix}", units[1]) for substance, amount in self.contents.items())
+
+        return round(numerator/denominator/mult, config.external_precision)
 
     @staticmethod
-    def create_stock_solution(solute: Substance, concentration: float, solvent: Substance, volume: float):
+    def create_solution(solute, solvent, name=None, **kwargs):
         """
-        Create a stock solution.
+        Create a solution.
 
-        Note: solids are assumed to have zero volume.
+        Two out of concentration, quantity, and total_quantity must be specified.
 
         Arguments:
             solute: What to dissolve.
-            concentration: Desired concentration in mol/L.
             solvent: What to dissolve with.
-            volume: Desired total volume in mL.
+            name: Optional name for new container.
+            concentration: Desired concentration. ('1 M', '0.1 umol/10 uL', etc.)
+            quantity: Desired quantity of solute. ('3 mL', '10 g')
+            total_quantity: Desired total quantity. ('3 mL', '10 g')
+
 
         Returns:
             New container with desired solution.
         """
+
         if not isinstance(solute, Substance):
             raise TypeError("Solute must be a Substance.")
-        if not isinstance(concentration, (int, float)):
-            raise TypeError("Concentration must be a float.")
         if not isinstance(solvent, Substance):
             raise TypeError("Solvent must be a Substance.")
+        if name and not isinstance(name, str):
+            raise TypeError("Name must be a str.")
+
+        if 'concentration' in kwargs and not isinstance(kwargs['concentration'], str):
+            raise TypeError("Concentration must be a str.")
+        if 'quantity' in kwargs and not isinstance(kwargs['quantity'], str):
+            raise TypeError("Quantity must be a str.")
+        if 'total_quantity' in kwargs and not isinstance(kwargs['total_quantity'], str):
+            raise TypeError("Total quantity must be a str.")
+        if ('concentration' in kwargs) + ('total_quantity' in kwargs) + ('quantity' in kwargs) != 2:
+            raise ValueError("Must specify two values out of concentration, quantity, and total quantity.")
+
+        if not name:
+            name = f"solution of {solute.name} in {solvent.name}"
+
+        if 'concentration' in kwargs and 'total_quantity' in kwargs:
+            concentration = kwargs['concentration']
+            total_quantity = kwargs['total_quantity']
+            quantity, quantity_unit = Unit.parse_quantity(total_quantity)
+            if quantity <= 0:
+                raise ValueError("Quantity must be positive.")
+
+            ratio, numerator, denominator = Unit.calculate_concentration_ratio(solute, concentration, solvent)
+
+            if ratio <= 0:
+                raise ValueError("Solution is impossible to create.")
+
+            if numerator == 'U':
+                if not solute.is_enzyme():
+                    raise TypeError("Solute must be an enzyme.")
+                solvent_quantity = Unit.convert(solvent, f"{quantity} {quantity_unit}", config.moles_prefix)
+                units = ratio * solvent_quantity
+                return Container(name,
+                                 initial_contents=((solute, f"{units} U"), (solvent, f"{quantity} {quantity_unit}")))
+
+            if quantity_unit == 'g':
+                ratio *= solute.mol_weight / solvent.mol_weight
+            elif quantity_unit == 'mol':
+                pass
+            elif quantity_unit == 'L':
+                ratio *= (solute.mol_weight / solute.density) / (solvent.mol_weight / solvent.density)
+            else:
+                raise ValueError("Invalid quantity unit.")
+
+            # x is quantity of solute in moles, y is quantity of solvent in moles
+
+            y = quantity / (1 + ratio)
+            x = quantity - y
+
+            assert x >= 0 and y >= 0
+            solution = Container(name,
+                                 initial_contents=((solute, f"{x} {quantity_unit}"), (solvent, f"{y} {quantity_unit}")))
+            return solution
+        if 'quantity' in kwargs and 'total_quantity' in kwargs:
+            result = Container(name, initial_contents=[(solute, kwargs['quantity'])])
+            result = result.fill_to(solvent, kwargs['total_quantity'])
+        else:  # 'quantity' and 'concentration'
+            concentration = Unit.calculate_concentration_ratio(solute, kwargs['concentration'], solvent)
+            quantity = Unit.convert(solute, kwargs['quantity'], concentration[1])
+            result = Container(name, initial_contents=[(solute, kwargs['quantity'])])
+            result._self_add(solvent, f"{quantity / concentration[0]} {concentration[1]}")
+        contents = []
+        for substance, value in result.contents.items():
+            value, unit = Unit.convert_from_storage_to_standard_format(substance, value)
+            contents.append(f"{value} {unit} of {substance.name}")
+        result.instructions = "Add " + ", ".join(contents) + " to a container."
+        return result
+
+    @staticmethod
+    def create_solution_from(source: Container, solute: Substance, concentration: str, solvent: Substance,
+                             quantity: str, name=None):
+        """
+        Create a diluted solution from an existing solution.
+
+
+        Arguments:
+            source: Solution to dilute.
+            solute: What to dissolve.
+            concentration: Desired concentration. ('1 M', '0.1 umol/10 uL', etc.)
+            solvent: What to dissolve with.
+            quantity: Desired total quantity. ('3 mL', '10 g')
+            name: Optional name for new container.
+
+        Returns:
+            Residual from the source container and a new container with desired solution.
+        """
+
+        if not isinstance(source, Container):
+            raise TypeError("Source must be a Container.")
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(concentration, str):
+            raise TypeError("Concentration must be a str.")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a Substance.")
+        if not isinstance(quantity, str):
+            raise TypeError("Quantity must be a str.")
+        if name and not isinstance(name, str):
+            raise TypeError("Name must be a str.")
+
+        quantity_value, quantity_unit = Unit.parse_quantity(quantity)
+        if quantity_value <= 0:
+            raise ValueError("Quantity must be positive.")
+
+        if solute not in source.contents:
+            raise ValueError(f"Container does not contain {solute.name}.")
+
+        if not name:
+            name = f"solution of {solute.name} in {solvent.name}"
+
+        new_ratio, numerator, denominator = Unit.calculate_concentration_ratio(solute, concentration, solvent)
+        current_ratio = source.contents[solute] / sum(source.contents[substance] for
+                                                      substance in source.contents if not substance.is_enzyme())
+        if new_ratio <= 0:
+            raise ValueError("Solution is impossible to create.")
+
+        if abs(new_ratio - current_ratio) <= 1e-6:
+            return deepcopy(source)
+
+        if new_ratio > current_ratio:
+            raise ValueError("Desired concentration is higher than current concentration.")
+
+        potential_solution = Container.create_solution(solute, solvent, concentration=concentration,
+                                                       total_quantity=quantity)
+        ratio = potential_solution.contents[solute] / source.contents[solute]
+        solution = Container(name, max_volume=f"{source.max_volume} {config.volume_prefix}")
+
+        residual, solution = Container.transfer(source, solution, f"{source.volume * ratio} {config.volume_prefix}")
+        solution = solution.fill_to(solvent, quantity)
+
+        contents = [(*Unit.convert_from_storage_to_standard_format(solution, solution.volume), source.name),
+                    (*Unit.convert_from_storage_to_standard_format(solvent, quantity), solvent.name)]
+        max_volume, unit = Unit.convert_from_storage_to_standard_format(solution, solution.max_volume)
+        solution.instructions = ("Add "
+                                 + ", ".join(f"{value} {unit} of {substance}" for value, unit, substance in contents)
+                                 + f" to a {max_volume} {unit} container.")
+
+        return residual, solution.fill_to(solvent, quantity)
+
+    def remove(self, what=Substance.LIQUID):
+        """
+        Removes substances from `Container`
+
+        Arguments:
+            what: What to remove. Can be a type of substance or a specific substance. Defaults to LIQUID.
+
+        Returns: New Container with requested substances removed.
+
+        """
+        new_container = deepcopy(self)
+        new_container.contents = {substance: value for substance, value in self.contents.items()
+                                  if what not in (substance._type, substance)}
+        new_container.volume = sum(Unit.convert_from(substance, value, 'U' if substance.is_enzyme() else
+        config.moles_prefix, config.volume_prefix) for substance, value in
+                                   new_container.contents.items())
+        new_container.instructions = self.instructions
+        classes = {Substance.SOLID: 'solid', Substance.LIQUID: 'liquid', Substance.ENZYME: 'enzyme'}
+        if what in classes:
+            new_container.instructions += f"Remove all {classes[what]}s."
+        else:
+            new_container.instructions += f"Remove all {what.name}s."
+        return new_container
+
+    def dilute(self, solute: Substance, concentration: str, solvent: Substance, name=None):
+        """
+        Dilutes `solute` in solution to `concentration`.
+
+        Args:
+            solute: Substance which is subject to dilution.
+            concentration: Desired concentration.
+            solvent: What to dilute with.
+            name: Optional name for new container.
+
+        Returns: A new container containing a solution with the desired concentration of `solute`.
+
+        """
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(concentration, str):
+            raise TypeError("Concentration must be a str.")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a substance.")
+        if name and not isinstance(name, str):
+            raise TypeError("New name must be a str.")
+        if solute not in self.contents:
+            raise ValueError(f"Container does not contain {solute.name}.")
+
+        new_ratio, numerator, denominator = Unit.calculate_concentration_ratio(solute, concentration, solvent)
+
+        if numerator == 'U':
+            if not solute.is_enzyme():
+                raise TypeError("Solute must be an enzyme.")
+
+        current_ratio = self.contents[solute] / sum(self.contents[substance] for
+                                                    substance in self.contents if not substance.is_enzyme())
+
+        if new_ratio <= 0:
+            raise ValueError("Solution is impossible to create.")
+
+        if abs(new_ratio - current_ratio) <= 1e-6:
+            return deepcopy(self)
+
+        if new_ratio > current_ratio:
+            raise ValueError("Desired concentration is higher than current concentration.")
+
+        current_umoles = Unit.convert_from_storage(self.contents.get(solvent, 0), 'umol')
+        required_umoles = Unit.convert_from_storage(self.contents[solute], 'umol') / new_ratio - current_umoles
+        new_volume = self.volume + Unit.convert(solvent, f"{required_umoles} umol", config.volume_prefix)
+
+        if new_volume > self.max_volume:
+            raise ValueError("Dilute solution will not fit in container.")
+
+        if name:
+            # Note: this copies the container twice
+            destination = deepcopy(self)
+            destination.name = name
+        else:
+            destination = self
+        needed_umoles = f"{required_umoles} umol"
+        result = destination._add(solvent, needed_umoles)
+        needed_volume, unit = Unit.get_human_readable_unit(Unit.convert(solvent, needed_umoles, 'umol'), 'L')
+        result.instructions += f"Dilute with {needed_volume} {unit} of {solvent.name}."
+        return result
+
+    def dilute_mols(self, solute: Substance, quantity_in_moles: float, solvent: Substance, volume: float, name=None):
+        """
+        Dilutes `solute` to achieve the desired quantity in moles and volume.
+
+        Args:
+            solute: Substance to be diluted.
+            quantity_in_moles: Desired quantity of the solute in moles.
+            solvent: Substance to dilute with.
+            volume: Desired volume of the solution in liters.
+            name: Optional name for the new container.
+
+        Returns: A new container containing a solution with the desired quantity of `solute` and volume.
+        """
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(quantity_in_moles, (int, float)):
+            raise TypeError("Quantity in moles must be a number.")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a substance.")
         if not isinstance(volume, (int, float)):
-            raise TypeError("Volume must be a float.")
+            raise TypeError("Volume must be a number.")
+        if name and not isinstance(name, str):
+            raise TypeError("New name must be a str.")
+        if solute not in self.contents:
+            raise ValueError(f"Container does not contain {solute.name}.")
 
-        container = Container(f"{solute.name} {concentration:.2}M", max_volume=volume)
-        moles_to_add = volume * concentration
-        if solute.is_enzyme():
-            raise TypeError("You can't add enzymes by molarity.")
-        if solute.is_solid():
-            container = container.add(solvent, container, f"{volume} mL")
-            container = container.add(solute, container, f"{round(moles_to_add, config.internal_precision)} mol")
-        else:  # Liquid
-            volume_to_add = round(moles_to_add * solute.mol_weight / (solute.density * 1000), config.external_precision)
-            container = container.add(solvent, container, f"{volume - volume_to_add} mL")
-            container = container.add(solute, container, f"{volume_to_add} mL")
-        return container
+        new_ratio, numerator, denominator = Unit.calculate_concentration_ratio_moles(solute, quantity_in_moles)
+        if numerator == 'U':
+            if not solute.is_enzyme():
+                raise TypeError("Solute must be an enzyme.")
+
+        current_ratio = self.contents[solute] / sum(
+            self.contents[substance] for substance in self.contents if not substance.is_enzyme())
+
+        if new_ratio <= 0:
+            raise ValueError("Solution is impossible to create.")
+
+        if abs(new_ratio - current_ratio) <= 1e-6:
+            return deepcopy(self)
+
+        if new_ratio > current_ratio:
+            raise ValueError("Desired concentration is higher than current concentration.")
+
+        current_umoles = Unit.convert_from_storage(self.contents[solvent], 'umol')
+        required_umoles = Unit.convert_from_storage(self.contents[solute], 'umol') / new_ratio - current_umoles
+        new_volume = self.volume + Unit.convert(solvent, f"{required_umoles} umol", config.volume_prefix)
+
+        if new_volume > self.max_volume:
+            raise ValueError("Dilute solution will not fit in the container.")
+
+        if name:
+            # Note: this copies the container twice
+            destination = deepcopy(self)
+            destination.name = name
+        else:
+            destination = self
+        result = destination._add(solvent, f"{required_umoles} umol")
+        needed_volume, unit = Unit.get_human_readable_unit(Unit.convert(solvent, f"{required_umoles} umol", 'umol'), 'L')
+        result.instructions += f"Dilute with {needed_volume} {unit} of {solvent.name}."
+        return result
+
+    def fill_to(self, solvent: Substance, quantity: str):
+        """
+        Fills container with `solvent` up to `quantity`.
+
+        Args:
+            solvent: Substance to use to fill.
+            quantity: Desired final quantity in container.
+
+        Returns: New Container with desired final `quantity`
+
+        """
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a Substance.")
+        if not isinstance(quantity, str):
+            raise TypeError("Quantity must be a str.")
+
+        quantity, quantity_unit = Unit.parse_quantity(quantity)
+        if quantity <= 0:
+            raise ValueError("Quantity must be positive.")
+        if quantity_unit not in ('L', 'g', 'mol'):
+            raise ValueError("We can only fill to mass or volume.")
+
+        current_quantity = sum(Unit.convert(substance, f"{value} {config.moles_prefix}", quantity_unit)
+                               for substance, value in self.contents.items() if not substance.is_enzyme())
+
+        required_quantity = quantity - current_quantity
+        result = self._add(solvent, f"{required_quantity} {quantity_unit}")
+        required_volume = Unit.convert(solvent, f"{required_quantity} {quantity_unit}", 'L')
+        required_volume, unit = Unit.get_human_readable_unit(required_volume, 'L')
+        result.instructions += f"Fill with {required_volume} {unit} of {solvent.name}."
+        return result
 
 
 class Plate:
@@ -631,13 +1228,13 @@ class Plate:
     The spatial arrangement must be rectangular. Immutable.
     """
 
-    def __init__(self, name: str, max_volume_per_well: float, make: str = "generic", rows=8, columns=12):
+    def __init__(self, name: str, max_volume_per_well: str, make: str = "generic", rows=8, columns=12):
         """
             Creates a generic plate.
 
             Attributes:
                 name: name of plate
-                max_volume_per_well: maximum volume of each well in uL
+                max_volume_per_well: maximum volume of each well. (50 uL)
                 make: name of this kind of plate
                 rows (int or list): number of rows or list of names of rows
                 columns (int or list): number of columns or list of names of columns
@@ -651,17 +1248,21 @@ def __init__(self, name: str, max_volume_per_well: float, make: str = "generic",
             raise ValueError("invalid plate make")
         self.make = make
 
+        if not isinstance(max_volume_per_well, str):
+            raise TypeError("Maximum volume must be a str, ('10 mL').")
+        max_volume_per_well, _ = Unit.parse_quantity(max_volume_per_well)
+
         if isinstance(rows, int):
             if rows < 1:
                 raise ValueError("illegal number of rows")
             self.n_rows = rows
             self.row_names = []
-            for n in range(1, rows + 1):
+            for row_num in range(1, rows + 1):
                 result = []
-                while n > 0:
-                    n -= 1
-                    result.append(chr(ord('A') + n % 26))
-                    n //= 26
+                while row_num > 0:
+                    row_num -= 1
+                    result.append(chr(ord('A') + row_num % 26))
+                    row_num //= 26
                 self.row_names.append(''.join(reversed(result)))
         elif isinstance(rows, list):
             if len(rows) == 0:
@@ -680,13 +1281,9 @@ def __init__(self, name: str, max_volume_per_well: float, make: str = "generic",
         else:
             raise ValueError("rows must be int or list")
 
-        try:
-            max_volume_per_well = float(max_volume_per_well)
-            if max_volume_per_well <= 0:
-                raise ValueError("max volume per well must be greater than zero")
-            self.max_volume_per_well = Unit.convert_to_storage(max_volume_per_well, 'uL')
-        except (ValueError, OverflowError) as exc:
-            raise ValueError(f"invalid max volume per well {max_volume_per_well}") from exc
+        if max_volume_per_well <= 0:
+            raise ValueError("max volume per well must be greater than zero")
+        self.max_volume_per_well = Unit.convert_to_storage(max_volume_per_well, 'L')
 
         if isinstance(columns, int):
             if columns < 1:
@@ -710,9 +1307,9 @@ def __init__(self, name: str, max_volume_per_well: float, make: str = "generic",
         else:
             raise ValueError("columns must be int or list")
 
-        self.wells = numpy.array([[Container(f"well {row + 1},{col + 1}",
-                                             max_volume=Unit.convert_from_storage(self.max_volume_per_well, 'mL'))
-                                   for col in range(self.n_columns)] for row in range(self.n_rows)])
+        self.wells = numpy.array([[Container(f"well {row},{col}",
+                                             max_volume=f"{max_volume_per_well} L")
+                                   for col in self.column_names] for row in self.row_names])
 
     def __getitem__(self, item) -> PlateSlicer:
         return PlateSlicer(self, item)
@@ -720,106 +1317,198 @@ def __getitem__(self, item) -> PlateSlicer:
     def __repr__(self):
         return f"Plate: {self.name}"
 
-    def volumes(self, substance: Substance = None) -> numpy.ndarray:
+    def volumes(self, substance: Substance = None, unit: str = None) -> numpy.ndarray:
         """
 
         Arguments:
+            unit: unit to return volumes in.
             substance: (optional) Substance to display volumes of.
 
         Returns:
-            numpy.ndarray of volumes for each well in uL
+            numpy.ndarray of volumes for each well in desired unit.
 
         """
-        if substance is None:
-            return numpy.vectorize(lambda elem: Unit.convert_from_storage(elem.volume, 'uL'))(self.wells)
-
-        if not isinstance(substance, Substance):
-            raise TypeError(f"Substance is not a valid type, {type(substance)}.")
-
-        def helper(elem):
-            """ Returns volume of elem. """
-            if substance.is_liquid() and substance in elem.contents:
-                return round(Unit.convert_from_storage(elem.contents[substance], 'uL'),
-                             config.external_precision)
-            return 0
-
-        return numpy.vectorize(helper)(self.wells)
+        if unit is None:
+            unit = config.default_volume_unit
+        return self[:].volumes(substance=substance, unit=unit)
 
     def substances(self):
         """
 
-        Returns: A set of substances present in the plate.
+        Returns: A set of substances present in the slice.
 
         """
-        substances_arr = numpy.vectorize(lambda elem: elem.contents.keys())(self.wells)
-        return set.union(*map(set, substances_arr.flatten()))
+        return self[:].substances()
 
-    def moles(self, substance: Substance) -> numpy.ndarray:
+    def moles(self, substance: Substance, unit: str = None) -> numpy.ndarray:
         """
+
         Arguments:
+            unit: unit to return moles in. ('mol', 'mmol', 'umol', etc.)
             substance: Substance to display moles of.
 
         Returns: moles of substance in each well.
         """
+        if unit is None:
+            unit = config.default_moles_unit
+        return self[:].moles(substance=substance, unit=unit)
 
-        if not isinstance(substance, Substance):
-            raise TypeError(f"Substance is not a valid type, {type(substance)}.")
+    def moles_dataframe(self, substance: Substance, unit: str = None, cmap: str = None):
+        """
 
-        def helper(elem):
-            """ Returns moles of substance in elem. """
-            if substance not in elem.contents:
-                return 0
-            if substance.is_liquid():
-                return round(Unit.convert_from_storage(elem.contents[substance], 'mL') * substance.density /
-                             substance.mol_weight, config.external_precision)
-            if substance.is_solid():
-                return round(Unit.convert_from_storage(elem.contents[substance], 'mol'), config.external_precision)
+        Arguments:
+            substance: Substance to display moles of.
+            unit: unit to return moles in. ('mol', 'mmol', 'umol', etc.)
+            cmap: Colormap to shade dataframe with.
 
-        return numpy.vectorize(helper, cache=True)(self.wells)
+        Returns: Shaded dataframe of number of moles of substance in each well.
 
-    def volume(self):
         """
-        Returns: total volume stored in plate in uL.
+        if unit is None:
+            unit = config.default_moles_unit
+        if cmap is None:
+            cmap = config.default_colormap
+        moles = self.moles(substance, unit)
+        dataframe = pandas.DataFrame(moles, columns=self.column_names, index=self.row_names)
+        return dataframe.style.format('{:.3f}').background_gradient(cmap, vmin=0, vmax=moles.max())
+
+    def volumes_dataframe(self, substance: Substance = None, unit: str = None, cmap: str = None):
         """
-        return self.volumes().sum()
 
-    @staticmethod
-    def add(source: Substance, destination: Plate | PlateSlicer, quantity: str):
+        Arguments:
+            unit: unit to return volumes in.
+            substance: (optional) Substance to display volumes of.
+            cmap: Colormap to shade dataframe with.
+
+        Returns: Shaded dataframe of volumes in each well.
+
         """
-        Add the given quantity ('10 mol') of the source substance to the destination.
+        if unit is None:
+            unit = config.default_volume_unit
+        if cmap is None:
+            cmap = config.default_colormap
+        volumes = self.volumes(substance, unit)
+        dataframe = pandas.DataFrame(volumes, columns=self.column_names, index=self.row_names)
+        return dataframe.style.format('{:.3f}').background_gradient(cmap, vmin=0, vmax=volumes.max())
 
+    def volume(self, unit: str = 'uL'):
+        """
         Arguments:
-            source: `Substance` to add.
-            destination: Plate or slice of a plate to add to.
-            quantity: How much to add.
-        Returns:
-            A new copy of `destination` plate.
+            unit: unit to return volumes in.
+
+        Returns: total volume stored in slice in uL.
         """
-        if not isinstance(destination, (Plate, PlateSlicer)):
-            raise TypeError("You can only use Plate.add to add to a Plate")
-        if isinstance(destination, Plate):
-            destination = destination[:]
-        # noinspection PyProtectedMember
-        return PlateSlicer._add(source, destination, quantity)
+        return self.volumes(unit=unit).sum()
 
     @staticmethod
-    def transfer(source, destination: Plate | PlateSlicer, volume):
+    def transfer(source: Container | Plate | PlateSlicer, destination: Plate | PlateSlicer, quantity: str):
         """
-        Move volume ('10 mL') from source to destination,
+        Move quantity ('10 mL', '5 mg') from source to destination,
         returning copies of the objects with amounts adjusted accordingly.
 
         Arguments:
             source: What to transfer.
             destination: Plate or slice of a plate to transfer to.
-            volume: How much to transfer.
+            quantity: How much to transfer.
 
         Returns:
             A tuple of (T, Plate) where T is the type of the source.
         """
         if not isinstance(destination, (Plate, PlateSlicer)):
             raise TypeError("You can only use Plate.transfer into a Plate")
+        if isinstance(destination, Plate):
+            destination = destination[:]
         # noinspection PyProtectedMember
-        return PlateSlicer._transfer(source, destination, volume)
+        return PlateSlicer._transfer(source, destination, quantity)
+
+    def remove(self, what=Substance.LIQUID):
+        """
+        Removes substances from `Plate`
+
+        Arguments:
+            what: What to remove. Can be a type of substance or a specific substance. Defaults to LIQUID.
+
+        Returns: New Plate with requested substances removed.
+
+        """
+        return self[:].remove(what)
+
+    def fill_to(self, solvent, quantity):
+        """
+        Fills all wells in plate with `solvent` up to `quantity`.
+
+        Args:
+            solvent: Substance to use to fill.
+            quantity: Desired final quantity in each well.
+
+        Returns: New Plate with desired final `quantity` in each well.
+
+        """
+        return self[:].fill_to(solvent, quantity)
+
+
+class RecipeStep:
+    def __init__(self, operator, frm, to, *operands):
+        self.operator = operator
+        self.frm = [frm]
+        self.to = [to]
+        self.operands = operands
+
+    def visualize(self, what, mode, unit, substance='all', cmap=None):
+        """
+
+        Provide visualization of what happened during the step.
+
+        Args:
+            what: 'source', 'destination', or 'both'
+            mode: 'delta', or 'final'
+            unit: Unit we are interested in. ('mmol', 'uL', 'mg')
+            substance: Substance we are interested in. ('all', 'water', 'ATP')
+            cmap: Colormap to use. Defaults to default_colormap from config.
+
+        Returns: A dataframe with the requested information or a list of dataframes if what is 'both'.
+        """
+
+        def helper(elem):
+            """ Returns volume of elem. """
+            if substance == 'all':
+                total = 0
+                for subs, amount in elem.contents.items():
+                    if subs.is_enzyme():
+                        total += Unit.convert(subs, f"{amount} U", unit)
+                    else:
+                        total += Unit.convert(subs, f"{amount} {config.moles_prefix}", unit)
+                return total
+            if substance in elem.contents:
+                quantity = f"{elem.contents[substance]} {config.moles_prefix if not substance.is_enzyme() else 'U'}"
+                return Unit.convert(substance, quantity, unit)
+            return 0
+
+        if what == 'both':
+            return [self.visualize('source', mode, unit), self.visualize('destination', mode, unit)]
+        if what == 'source':
+            what = self.frm
+        elif what == 'destination':
+            what = self.to
+        else:
+            raise ValueError("What must be source, destination, or both.")
+
+        assert mode == 'delta' or mode == 'final'
+
+        if not isinstance(what[0], Plate):
+            return None
+        data = numpy.vectorize(helper, cache=True, otypes='d')(what[1].wells)
+        if mode == 'delta':
+            data -= numpy.vectorize(helper, cache=True, otypes='d')(what[0].wells)
+            if cmap is None:
+                cmap = config.default_diverging_colormap
+        if cmap is None:
+            cmap = config.default_colormap
+        dataframe = pandas.DataFrame(data, columns=what[0].column_names, index=what[0].row_names)
+        extreme = max(abs(numpy.min(data)), abs(numpy.max(data)))
+        return dataframe.style.format('{:.3f}').background_gradient(cmap, vmin=-extreme,
+            vmax=extreme).set_caption(
+            (substance.name if isinstance(substance, Substance) else substance) + f"  ({unit})")
 
 
 class Recipe:
@@ -841,7 +1530,7 @@ class Recipe:
     def __init__(self):
         self.indexes = {}
         self.results = []
-        self.steps = []
+        self.steps: list[RecipeStep] = []
         self.locked = False
         self.used = set()
 
@@ -853,37 +1542,16 @@ def uses(self, *args):
             raise RuntimeError("This recipe is locked.")
         for arg in args:
             if arg not in self.indexes:
-                self.indexes[arg] = len(self.results)
-                if isinstance(arg, Substance):
-                    self.results.append(arg)
-                else:
+                # TODO: Do we need to keep track of substances?
+                if isinstance(arg, (Container, Plate)):
+                    self.indexes[arg] = len(self.results)
                     self.results.append(deepcopy(arg))
         return self
 
-    def add(self, source, destination, quantity):
+    def transfer(self, source: Container, destination: Container | Plate, quantity: str):
         """
-        Adds a step to the recipe which will move the given quantity ('10 mol')
-        of the source substance to the destination.
-
-        """
-        if self.locked:
-            raise RuntimeError("This recipe is locked.")
-        if not isinstance(source, Substance):
-            raise TypeError("Invalid source type.")
-        if not isinstance(destination, (Container, Plate, PlateSlicer)):
-            raise TypeError("Invalid destination type.")
-        if (destination.plate if isinstance(destination, PlateSlicer) else destination) not in self.indexes:
-            name = destination.plate.name if isinstance(destination, PlateSlicer) else destination.name
-            raise ValueError(f"Destination {name} has not been previously declared for use.")
-        if isinstance(destination, Plate):
-            destination = destination[:]
-        self.steps.append(('add', source, destination, quantity))
-        return self
-
-    def transfer(self, source: Container, destination: Container | Plate, volume: str):
-        """
-        Adds a step to the recipe which will move volume from source to destination.
-        Note that all Substances in the source will be transferred in proportion to their volumetric ratios.
+        Adds a step to the recipe which will move quantity from source to destination.
+        Note that all Substances in the source will be transferred in proportion to their respective ratios.
 
         """
         if self.locked:
@@ -897,21 +1565,22 @@ def transfer(self, source: Container, destination: Container | Plate, volume: st
         if (destination.plate if isinstance(destination, PlateSlicer) else destination) not in self.indexes:
             name = destination.plate.name if isinstance(destination, PlateSlicer) else destination.name
             raise ValueError(f"Destination {name} has not been previously declared for use.")
+        if not isinstance(quantity, str):
+            raise TypeError("Volume must be a str. ('5 mL')")
         if isinstance(source, Plate):
             source = source[:]
         if isinstance(destination, Plate):
             destination = destination[:]
-        self.steps.append(('transfer', source, destination, volume))
-        return self
+        self.steps.append(RecipeStep('transfer', source, destination, quantity))
 
-    def create_container(self, name, max_volume, initial_contents=None):
+    def create_container(self, name: str, max_volume: str = 'inf L', initial_contents=None):
 
         """
         Adds a step to the recipe which creates a container.
 
         Arguments:
             name: Name of container
-            max_volume: Maximum volume that can be stored in the container in mL
+            max_volume: Maximum volume that can be stored in the container. ('10 mL')
             initial_contents: (optional) Iterable of tuples of the form (Substance, quantity)
 
         Returns:
@@ -919,37 +1588,202 @@ def create_container(self, name, max_volume, initial_contents=None):
         """
         if self.locked:
             raise RuntimeError("This recipe is locked.")
-        new_container = Container(name, max_volume)
-        self.uses(new_container)
+        if not isinstance(name, str):
+            raise TypeError("Name must be a str.")
+        if not isinstance(max_volume, str):
+            raise TypeError("Maximum volume must be a str.")
+
         if initial_contents:
+            if not isinstance(initial_contents, Iterable):
+                raise TypeError("Initial contents must be iterable.")
+            if not all(isinstance(elem, tuple) and len(elem) == 2 for elem in initial_contents):
+                raise TypeError("Elements of initial_contents must be of the form (Substance, quantity.)")
             for substance, quantity in initial_contents:
                 if not isinstance(substance, Substance):
-                    raise ValueError("Containers can only be created from substances.")
-                self.steps.append(('add', substance, new_container, quantity))
+                    raise TypeError("Containers can only be created from substances.")
+                if not isinstance(quantity, str):
+                    raise TypeError("Quantity must be a str. ('10 mL')")
+        new_container = Container(name, max_volume)
+        self.uses(new_container)
+        self.steps.append(RecipeStep('create_container', None, new_container, max_volume, initial_contents))
+
         return new_container
 
-    def create_stock_solution(self, what: Substance, concentration: float, solvent: Substance, volume: float):
+    def create_solution(self, solute, solvent, name=None, **kwargs):
         """
+        Adds a step to the recipe which creates a solution.
 
-        Adds a step to the recipe which creates a stock solution.
+        Two out of concentration, quantity, and total_quantity must be specified.
 
         Arguments:
-            what: What to dilute.
-            concentration: Desired concentration in mol/L
-            solvent: What to dilute with.
-            volume: Desired total volume in mL.
+            solute: What to dissolve.
+            solvent: What to dissolve with.
+            name: Optional name for new container.
+            concentration: Desired concentration. ('1 M', '0.1 umol/10 uL', etc.)
+            quantity: Desired quantity of solute. ('3 mL', '10 g')
+            total_quantity: Desired total quantity. ('3 mL', '10 g')
+
 
         Returns:
             A new Container so that it may be used in later recipe steps.
-
         """
-        if self.locked:
-            raise RuntimeError("This recipe is locked.")
-        new_container = Container("{what.name} {concentration:.2}M", max_volume=volume)
+
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a Substance.")
+        if name and not isinstance(name, str):
+            raise TypeError("Name must be a str.")
+
+        if 'concentration' in kwargs and not isinstance(kwargs['concentration'], str):
+            raise TypeError("Concentration must be a str.")
+        if 'quantity' in kwargs and not isinstance(kwargs['quantity'], str):
+            raise TypeError("Quantity must be a str.")
+        if 'total_quantity' in kwargs and not isinstance(kwargs['total_quantity'], str):
+            raise TypeError("Total quantity must be a str.")
+        if ('concentration' in kwargs) + ('total_quantity' in kwargs) + ('quantity' in kwargs) != 2:
+            raise ValueError("Must specify two values out of concentration, quantity, and total quantity.")
+
+        if not name:
+            name = f"solution of {solute.name} in {solvent.name}"
+
+        new_container = Container(name)
         self.uses(new_container)
-        self.steps.append(('stock', new_container, what, concentration, solvent, volume))
+        self.steps.append(RecipeStep('solution', None, new_container, solute, solvent, kwargs))
+
         return new_container
 
+    def create_solution_from(self, source: Container, solute: Substance, concentration: str, solvent: Substance,
+                             quantity: str, name=None):
+        """
+        Adds a step to create a diluted solution from an existing solution.
+
+
+        Arguments:
+            source: Solution to dilute.
+            solute: What to dissolve.
+            concentration: Desired concentration. ('1 M', '0.1 umol/10 uL', etc.)
+            solvent: What to dissolve with.
+            quantity: Desired total quantity. ('3 mL', '10 g')
+            name: Optional name for new container.
+
+        Returns:
+            A new Container so that it may be used in later recipe steps.
+        """
+
+        if not isinstance(source, Container):
+            raise TypeError("Source must be a Container.")
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(concentration, str):
+            raise TypeError("Concentration must be a str.")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a Substance.")
+        if not isinstance(quantity, str):
+            raise TypeError("Quantity must be a str.")
+        if name and not isinstance(name, str):
+            raise TypeError("Name must be a str.")
+
+        quantity_value, quantity_unit = Unit.parse_quantity(quantity)
+        if quantity_value <= 0:
+            raise ValueError("Quantity must be positive.")
+
+        if not name:
+            name = f"solution of {solute.name} in {solvent.name}"
+
+        new_ratio, numerator, denominator = Unit.calculate_concentration_ratio(solute, concentration, solvent)
+        if new_ratio <= 0:
+            raise ValueError("Solution is impossible to create.")
+
+        new_container = Container(name, max_volume=f"{source.max_volume} {config.volume_prefix}")
+        self.uses(new_container, source)
+        self.steps.append(RecipeStep('solution_from', source, new_container, solute, concentration, solvent, quantity))
+
+        return new_container
+
+    def remove(self, destination: Container | Plate | PlateSlicer, what=Substance.LIQUID):
+        """
+        Adds a step to removes substances from destination.
+
+        Arguments:
+            destination: What to remove from.
+            what: What to remove. Can be a type of substance or a specific substance. Defaults to LIQUID.
+        """
+
+        if isinstance(destination, PlateSlicer):
+            if destination.plate not in self.indexes:
+                raise ValueError(f"Destination {destination.plate.name} has not been previously declared for use.")
+        elif isinstance(destination, (Container, Plate)):
+            if destination not in self.indexes:
+                raise ValueError(f"Destination {destination.name} has not been previously declared for use.")
+        else:
+            raise TypeError(f"Invalid destination type: {type(destination)}")
+
+        self.steps.append(RecipeStep('remove', None, destination, what))
+
+    def dilute(self, destination: Container, solute: Substance,
+               concentration: str, solvent: Substance, new_name=None):
+        """
+        Adds a step to dilute `solute` in `destination` to `concentration`.
+
+        Args:
+            destination: Container to dilute.
+            solute: Substance which is subject to dilution.
+            concentration: Desired concentration in mol/L.
+            solvent: What to dilute with.
+            new_name: Optional name for new container.
+        """
+
+        if not isinstance(solute, Substance):
+            raise TypeError("Solute must be a Substance.")
+        if not isinstance(concentration, str):
+            raise TypeError("Concentration must be a float.")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a substance.")
+        if new_name and not isinstance(new_name, str):
+            raise TypeError("New name must be a str.")
+        if not isinstance(destination, Container):
+            raise TypeError("Destination must be a container.")
+        if destination not in self.indexes:
+            raise ValueError(f"Destination {destination.name} has not been previously declared for use.")
+        if solute not in destination.contents:
+            raise ValueError(f"Container does not contain {solute.name}.")
+
+        ratio, *_ = Unit.calculate_concentration_ratio(solute, concentration, solvent)
+        if ratio <= 0:
+            raise ValueError("Concentration is impossible to create.")
+
+        if solute.is_enzyme():
+            # TODO: Support this.
+            raise ValueError("Not currently supported.")
+
+        self.steps.append(RecipeStep('dilute', None, destination, solute, concentration, solvent, new_name))
+
+    def fill_to(self, destination: Container, solvent: Substance, quantity: str):
+        """
+        Adds a step to fill `destination` container/plate/slice with `solvent` up to `quantity`.
+
+        Args:
+            destination: Container/Plate/Slice to fill.
+            solvent: Substance to use to fill.
+            quantity: Desired final quantity in container.
+
+        """
+        if isinstance(destination, PlateSlicer):
+            if destination.plate not in self.indexes:
+                raise ValueError(f"Destination {destination.plate.name} has not been previously declared for use.")
+        elif isinstance(destination, (Container, Plate)):
+            if destination not in self.indexes:
+                raise ValueError(f"Destination {destination.name} has not been previously declared for use.")
+        else:
+            raise TypeError(f"Invalid destination type: {type(destination)}")
+        if not isinstance(solvent, Substance):
+            raise TypeError("Solvent must be a substance.")
+        if not isinstance(quantity, str):
+            raise TypeError("Quantity must be a str.")
+
+        self.steps.append(RecipeStep('fill_to', None, destination, solvent, quantity))
+
     def bake(self):
         """
         Completes steps stored in recipe.
@@ -960,31 +1794,31 @@ def bake(self):
             Copies of all used objects in the order they were declared.
 
         """
-        for operation, *rest in self.steps:
-            if operation == 'add':
-                frm, to, quantity = rest
-                to_index = self.indexes[to] if not isinstance(to, PlateSlicer) else self.indexes[to.plate]
-
-                # containers and such can change while building the recipe
-                if isinstance(to, PlateSlicer):
-                    new_to = deepcopy(to)
-                    new_to.plate = self.results[to_index]
-                    to = new_to
-                else:
-                    to = self.results[to_index]
-
-                if isinstance(to, Container):
-                    to = Container.add(frm, to, quantity)
-                elif isinstance(to, PlateSlicer):
-                    to = Plate.add(frm, to, quantity)
-
-                self.results[to_index] = to
+        if self.locked:
+            raise RuntimeError("Recipe has already been baked.")
+
+        # for operation, *rest in self.steps:
+        for step in self.steps:
+            operator = step.operator
+            if operator == 'create_container':
+                dest = step.to[0]
+                step.frm.append(None)
+                max_volume, initial_contents = step.operands
+                to_index = self.indexes[dest]
+                step.to[0] = self.results[to_index]
                 self.used.add(to_index)
-            elif operation == 'transfer':
-                frm, to, quantity = rest
+                self.results[to_index] = Container(dest.name, max_volume, initial_contents)
+                step.to.append(self.results[to_index])
+            elif operator == 'transfer':
+                frm = step.frm[0]
+                dest = step.to[0]
+                quantity, = step.operands
                 # used items can change in a recipe
                 frm_index = self.indexes[frm] if not isinstance(frm, PlateSlicer) else self.indexes[frm.plate]
-                to_index = self.indexes[to] if not isinstance(to, PlateSlicer) else self.indexes[to.plate]
+                to_index = self.indexes[dest] if not isinstance(dest, PlateSlicer) else self.indexes[dest.plate]
+
+                step.frm[0] = self.results[frm_index]
+                step.to[0] = self.results[to_index]
 
                 self.used.add(frm_index)
                 self.used.add(to_index)
@@ -997,30 +1831,93 @@ def bake(self):
                 else:
                     frm = self.results[frm_index]
 
-                if isinstance(to, PlateSlicer):
-                    new_to = deepcopy(to)
+                if isinstance(dest, PlateSlicer):
+                    new_to = deepcopy(dest)
                     new_to.plate = self.results[to_index]
-                    to = new_to
+                    dest = new_to
                 else:
-                    to = self.results[to_index]
-
-                if isinstance(frm, Substance):  # Adding a substance is handled differently
-                    if isinstance(to, Container):
-                        to = Container.add(frm, to, quantity)
-                    elif isinstance(to, PlateSlicer):
-                        to = Plate.add(frm, to, quantity)
-                elif isinstance(to, Container):
-                    frm, to = Container.transfer(frm, to, quantity)
-                elif isinstance(to, PlateSlicer):
-                    frm, to = Plate.transfer(frm, to, quantity)
+                    dest = self.results[to_index]
+
+                if isinstance(dest, Container):
+                    frm, dest = Container.transfer(frm, dest, quantity)
+                elif isinstance(dest, PlateSlicer):
+                    frm, dest = Plate.transfer(frm, dest, quantity)
 
                 self.results[frm_index] = frm
-                self.results[to_index] = to
-            elif operation == 'stock':
-                to, what, concentration, solvent, volume = rest
-                to_index = self.indexes[to]
+                self.results[to_index] = dest
+                step.frm.append(frm)
+                step.to.append(dest)
+            elif operator == 'solution':
+                dest = step.to[0]
+                step.frm.append(None)
+                solute, solvent, kwargs = step.operands
+                # solute, concentration, solvent, quantity = step.operands
+                to_index = self.indexes[dest]
+                step.to[0] = self.results[to_index]
+                self.used.add(to_index)
+                self.results[to_index] = Container.create_solution(solute, solvent, dest.name, **kwargs)
+                step.to.append(self.results[to_index])
+            elif operator == 'solution_from':
+                source = step.frm[0]
+                dest = step.to[0]
+                solute, concentration, solvent, quantity = step.operands
+                frm_index = self.indexes[source]
+                to_index = self.indexes[dest]
+                step.frm[0] = self.results[frm_index]
+                step.to[0] = self.results[to_index]
+
+                self.used.add(frm_index)
+                self.used.add(to_index)
+                source = self.results[frm_index]
+                self.results[frm_index], self.results[to_index] = Container.create_solution_from(source, solute,
+                                                                                                 concentration, solvent,
+                                                                                                 quantity, dest.name)
+                step.frm.append(self.results[frm_index])
+                step.to.append(self.results[to_index])
+
+            elif operator == 'remove':
+                dest = step.to[0]
+                step.frm.append(None)
+                what, = step.operands
+                to_index = self.indexes[dest] if not isinstance(dest, PlateSlicer) else self.indexes[dest.plate]
+                step.to[0] = self.results[to_index]
+                self.used.add(to_index)
+
+                if isinstance(dest, PlateSlicer):
+                    new_to = deepcopy(dest)
+                    new_to.plate = self.results[to_index]
+                    dest = new_to
+                else:
+                    dest = self.results[to_index]
+
+                self.results[to_index] = dest.remove(what)
+                step.to.append(self.results[to_index])
+            elif operator == 'dilute':
+                dest = step.to[0]
+                solute, concentration, solvent, new_name = step.operands
+                step.fxrm.append(None)
+                to_index = self.indexes[dest]
+                step.to[0] = self.results[to_index]
                 self.used.add(to_index)
-                self.results[to_index] = Container.create_stock_solution(what, concentration, solvent, volume)
+                self.results[to_index] = dest.dilute(solute, concentration, solvent, new_name)
+                step.to.append(self.results[to_index])
+            elif operator == 'fill_to':
+                dest = step.to[0]
+                solvent, quantity = step.operands
+                step.frm.append(None)
+                to_index = self.indexes[dest] if not isinstance(dest, PlateSlicer) else self.indexes[dest.plate]
+                step.to[0] = self.results[to_index]
+                self.used.add(to_index)
+
+                if isinstance(dest, PlateSlicer):
+                    new_to = deepcopy(dest)
+                    new_to.plate = self.results[to_index]
+                    dest = new_to
+                else:
+                    dest = self.results[to_index]
+
+                self.results[to_index] = dest.fill_to(solvent, quantity)
+                step.to.append(self.results[to_index])
 
         if len(self.used) != len(self.indexes):
             raise ValueError("Something declared as used wasn't used.")
@@ -1048,8 +1945,7 @@ def array(self, array):
     def _add(frm, to, quantity):
         to = copy(to)
         to.plate = deepcopy(to.plate)
-        result = numpy.vectorize(lambda elem: Container.add(frm, elem, quantity), cache=True)(to.get())
-        to.set(result)
+        to.apply(lambda elem: elem._add(frm, quantity))
         return to.plate
 
     @staticmethod
@@ -1064,10 +1960,7 @@ def helper_func(elem):
                 return elem
 
             frm_array = [frm]
-            result = numpy.vectorize(helper_func, cache=True)(to.get())
-            if isinstance(to.get(), Container):  # A zero-dim array was returned.
-                result = result.item()
-            to.set(result)
+            to.apply(helper_func)
             return frm_array[0], to.plate
         if not isinstance(frm, (Plate, PlateSlicer)):
             raise TypeError("Invalid source type.")
@@ -1076,44 +1969,57 @@ def helper_func(elem):
         frm = copy(frm)
 
         if to.plate != frm.plate:
+            different = True
             to.plate = deepcopy(to.plate)
             frm.plate = deepcopy(frm.plate)
         else:
+            different = False
             to.plate = frm.plate = deepcopy(to.plate)
 
         if frm.size == 1:
             # Source from the single element in frm
-            if frm.shape != ():
-                raise RuntimeError("Shape of source should have been ()")
+            if frm.shape != (1, 1):
+                raise RuntimeError("Shape of source should have been (1, 1)")
 
             def helper_func(elem):
                 """ @private """
-                frm_array[0], elem = Container.transfer(frm_array[0], elem, quantity)
+                frm_array[0][0], elem = Container.transfer(frm_array[0][0], elem, quantity)
+                if different:
+                    instructions = elem.instructions.splitlines()
+                    instructions[-1] = instructions[-1].replace(frm_array[0][0].name,
+                                                                frm.plate.name + " " + frm_array[0][0].name, 1)
+                    elem.instructions = "\n".join(instructions)
+
                 return elem
 
-            frm_array = [frm.get()]
-            result = numpy.vectorize(helper_func, cache=True)(to.get())
-            to.set(result)
-            frm.set(frm_array[0])
+            frm_array = frm.get()
+            to.apply(helper_func)
 
         elif to.size == 1:
             #  Replace the single element in self
-            if to.shape != ():
-                raise RuntimeError("Shape of source should have been ()")
+            if to.shape != (1, 1):
+                raise RuntimeError("Shape of source should have been (1, 1)")
 
             def helper_func(elem):
                 """ @private """
-                elem, to_array[0] = to_array[0].transfer(elem, quantity)
+                elem, to_array[0][0] = to_array[0][0].transfer(elem, quantity)
+                instructions = to_array[0][0].instructions.splitlines()
+                instructions[-1] = instructions[-1].replace(elem.name, frm.plate.name + " " + elem.name, 1)
+                elem.instructions = "\n".join(instructions)
                 return elem
 
-            to_array = [to.get()]
-            frm.set(numpy.vectorize(helper_func, cache=True)(frm.get()))
-            to.set(to_array[0])
+            to_array = to.get()
+            frm.apply(helper_func)
 
         elif frm.size == to.size and frm.shape == to.shape:
             def helper(elem1, elem2):
                 """ @private """
-                return Container.transfer(elem1, elem2, quantity)
+                elem1, elem2 = Container.transfer(elem1, elem2, quantity)
+                if different:
+                    instructions = elem2.instructions.splitlines()
+                    instructions[-1] = instructions[-1].replace(elem1.name, frm.plate.name + " " + elem1.name, 1)
+                    elem2.instructions = "\n".join(instructions)
+                return elem1, elem2
 
             func = numpy.frompyfunc(helper, 2, 2)
             frm_result, to_result = func(frm.get(), to.get())
@@ -1123,3 +2029,90 @@ def helper(elem1, elem2):
             raise ValueError("Source and destination slices must be the same size and shape.")
 
         return frm.plate, to.plate
+
+    def volumes(self, substance: Substance = None, unit: str = 'uL') -> numpy.ndarray:
+        """
+
+        Arguments:
+            unit:  unit to return volumes in.
+            substance: (optional) Substance to display volumes of.
+
+        Returns:
+            numpy.ndarray of volumes for each well in uL
+
+        """
+        if substance is None:
+            return numpy.vectorize(lambda elem: Unit.convert_from_storage(elem.volume, unit), cache=True,
+                                   otypes='d')(self.get())
+
+        if not isinstance(substance, Substance):
+            raise TypeError(f"Substance is not a valid type, {type(substance)}.")
+
+        def helper(elem):
+            """ Returns volume of elem. """
+            if substance in elem.contents:
+                quantity = f"{elem.contents[substance]} {config.moles_prefix}"
+                return Unit.convert(substance, quantity, unit)
+            return 0
+
+        return numpy.vectorize(helper, cache=True, otypes='d')(self.get())
+
+    def substances(self):
+        """
+
+        Returns: A set of substances present in the plate.
+
+        """
+        substances_arr = numpy.vectorize(lambda elem: elem.contents.keys(), cache=True)(self.get())
+        return set.union(*map(set, substances_arr.flatten()))
+
+    def moles(self, substance: Substance, unit: str = 'mol') -> numpy.ndarray:
+        """
+        Arguments:
+            unit: unit to return moles in. ('mol', 'mmol', 'umol', etc.)
+            substance: Substance to display moles of.
+
+        Returns: moles of substance in each well.
+        """
+
+        if not isinstance(substance, Substance):
+            raise TypeError(f"Substance is not a valid type, {type(substance)}.")
+
+        def helper(elem):
+            """ Returns moles of substance in elem. """
+            if substance not in elem.contents:
+                return 0
+            quantity = f"{elem.contents[substance]} {config.moles_prefix}"
+            return Unit.convert(substance, quantity, unit)
+
+        return numpy.vectorize(helper, cache=True, otypes='d')(self.get())
+
+    def remove(self, what=Substance.LIQUID):
+        """
+        Removes substances from slice
+
+        Arguments:
+            what: What to remove. Can be a type of substance or a specific substance. Defaults to LIQUID.
+
+        Returns: New Plate with requested substances removed.
+
+        """
+        self.plate = deepcopy(self.plate)
+        self.apply(lambda elem: elem.remove(what))
+        return self.plate
+
+    def fill_to(self, solvent, quantity):
+        """
+        Fills all wells in slice with `solvent` up to `quantity`.
+
+        Args:
+            solvent: Substance to use to fill.
+            quantity: Desired final quantity in each well.
+
+        Returns: New Plate with desired final `quantity` in each well.
+
+        """
+        self.plate = deepcopy(self.plate)
+        self.apply(lambda elem: elem.fill_to(solvent, quantity))
+
+        return self.plate
diff --git a/pyplate/sample_notebook.ipynb b/pyplate/sample_notebook.ipynb
new file mode 100644
index 0000000..bdf0449
--- /dev/null
+++ b/pyplate/sample_notebook.ipynb
@@ -0,0 +1,60 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "vscode": {
+     "languageId": "plaintext"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "        filtered_factors = {}\n",
+    "        for factor in self.factors:\n",
+    "            if factors[factor.name] == \"all\":\n",
+    "                filtered_factors[factor.name] = factor.possible_values\n",
+    "            else:\n",
+    "                filtered_factors[factor.name] = [value for value in factor.possible_values if value in factors[factor.name]]\n",
+    "        \n",
+    "        # Initialize blocks\n",
+    "        blocks = []\n",
+    "\n",
+    "        # If needed form blocking combinations\n",
+    "        #blocking_combinations = list(itertools.product(*[filtered_factors[name] for name in blocking_factors]))\n",
+    "\n",
+    "        # Iterate over each combination of blocking factors\n",
+    "        for block_comb in blocking_factors:\n",
+    "            block = []\n",
+    "            \n",
+    "            # Create a dictionary for the current blocking combination\n",
+    "            block_dict = dict(zip(blocking_factors, block_comb))\n",
+    "\n",
+    "            # Filter out blocking factors and create combinations of other factors\n",
+    "            other_factors = {name: values for name, values in filtered_factors.items() if name not in blocking_factors}\n",
+    "            other_combinations = list(itertools.product(*other_factors.values()))\n",
+    "\n",
+    "            # Generate experiments for each combination\n",
+    "            for comb in other_combinations:\n",
+    "                factors_dict = dict(zip(other_factors.keys(), comb))\n",
+    "                factors_dict.update(block_dict)  # Add blocking factor values\n",
+    "\n",
+    "                # Create replicates for each unique combination\n",
+    "                for rep in range(n_replicates):\n",
+    "                    experiment = Experiment(factors=factors_dict, replicate_idx=rep+1)\n",
+    "                    block.append(experiment)\n",
+    "\n",
+    "            blocks.append(block)\n",
+    "\n",
+    "        return blocks"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/pyplate/slicer.py b/pyplate/slicer.py
index cf5c4ac..51e9191 100644
--- a/pyplate/slicer.py
+++ b/pyplate/slicer.py
@@ -1,59 +1,57 @@
+from copy import copy
 from functools import cached_property
 from typing import Union, Tuple, List
+
+import numpy
 import numpy as np
 
 __docformat__ = "google"
 
-# Valid slices:
-#
-# Label = Union(str, int)
-#
-# Slice                                             ->      item
-#       shape                                               size
-#
-# array['A:1'] or array['A', '1'] or array[1, 1]    ->      str or Tuple[Label, Label]
-#       ()                                                  1
-#
-# array[i:j]                                        ->      Slice
-#       (n, n_cols)                                         n * n_cols
-#
-# array[i:j, 1]                                     ->      Tuple[Slice, Label]
-#       (n,)                                                n
-#
-# array[1, i:j]                                     ->      Tuple[Label, Slice]
-#       (n,)                                                n
-#
-# array[i:j, k:l]                                   ->      Tuple[Slice, Slice]
-#       (n, o)                                              n * o
-#
-# array[['A:1', ('A', '1'), (1,1)]]                 ->      List[Union[Label, Tuple[Label, Label]]
-#       (3,)                                                3
-
-
 Label = Union[str, int]
 """@private"""
 Single = Union[Label, Tuple[Label, Label]]
 """@private"""
 Singles = Union[List[Single], Single]
 """@private"""
-Slice = Union[Singles, slice, Tuple[slice, Label], Tuple[Label, slice], Tuple[slice, slice]]
+Slice = Union[Singles, slice, Tuple[slice, Label], Tuple[Label, slice],
+Tuple[slice, slice], Tuple[List[Single], slice], Tuple[slice, List[Single]]]
 """@private"""
 
 
 class Slicer:
     """
-    Helper class to provide slicing for composed numpy ndarrays.
-
-    Extend the class to make it useful for a particular application
-
-    All integer indexes are one-based. i.e. array[1,1] is the upper left element.
-    Does not currently support negative indexing.
+    This is a helper class designed to facilitate slicing operations on composed numpy ndarrays.
+    To use this class for a specific application, it needs to be extended.
+    The class supports one-based integer indexing, meaning array[1,1] refers to the upper left element.
+    Negative indexing is not currently supported.
+    The item to be sliced can be specified in several ways:
+    - A string in the format 'A:1' (for a single element)
+    - A tuple in the format (1, 1) (for a single element)
+    - A list of strings in the format 'A:1', or tuples in the format (1,1) (for multiple elements)
+    - A single label or integer (for a row)
+    - A label or integer and a slice (for a row and multiple columns)
+    - A slice and a label or integer (for multiple rows and a column)
+    - A slice
+    - Two slices (for a row and column)
+
+    Examples of valid slicing operations include:
+    - Single string: ['A:1'], ['1:1']
+    - Single tuple: [(1, 1)], [('A', '1')]
+    - List of strings or tuples: [['A:1', 'B:2', (1,1), (2,2), ('A','1'), ('B','2')]]
+    - Single label or integer: [1], ['A']
+    - Label or integer and a slice: [1, 1:], ['A', 1:], [1, '1':]
+    - Slice and a label or integer: [1:, 1], ['A':, 1], [1:, '1'], ['A':, '1']
+    - A slice: [1:], ['A':]
+    - Two slices: [1:, 1:], ['A':, 1:], [1:, '1':], ['A':, '1':]
+
+    The resulting slices will be either a single view or a list of views.
     """
 
-    def __init__(self, array_obj: np.ndarray, row_labels, col_labels, item):
+    def __init__(self, array_obj: np.ndarray, row_labels: tuple | list, col_labels: tuple | list,
+                 item: Slice | list):
         """
 
-        Args:\
+        Args:
             array_obj (numpy.ndarry): Array to slice.
             row_labels (list): Row labels.
             col_labels (list): Column labels.
@@ -77,15 +75,56 @@ def __init__(self, array_obj: np.ndarray, row_labels, col_labels, item):
             raise ValueError("Length of labels must match shape of array.")
 
         if isinstance(item, str):
-            self.slices = self.parse_single(item)
+            if ':' in item:
+                row, col = self.parse_single(item)
+                self.slices = (slice(row, row + 1), slice(col, col + 1))
+            else:
+                row = self.resolve_labels(item, self.row_labels)
+                self.slices = (slice(row, row + 1), slice(None))
         elif isinstance(item, list):
-            self.slices = list(map(self.parse_single, item))
+            self.slices = []
+            for elem in item:
+                if isinstance(elem, tuple):
+                    row, col = self.parse_tuple(elem)
+                    self.slices.append((slice(row, row + 1), slice(col, col + 1)))
+                elif isinstance(elem, str) and ':' in elem:
+                    row, col = self.parse_single(elem)
+                    self.slices.append((slice(row, row + 1), slice(col, col + 1)))
+                else:
+                    raise TypeError("Invalid slice.")
         elif isinstance(item, int):
             if not 1 <= item <= self.n_rows:
                 raise ValueError("Row index out of range.")
-            self.slices = item - 1
+            self.slices = (slice(item - 1, item), slice(None))
+        elif isinstance(item, slice):
+            self.slices = (self.parse_slice(item, self.row_labels), slice(None))
+        elif isinstance(item, tuple):
+            if len(item) == 1 and isinstance(item[0], slice):
+                self.slices = (Slicer.parse_slice(item[0], self.row_labels), slice(None))
+            if len(item) == 2:
+                if isinstance(item[0], (int, str)) and isinstance(item[1], (int, str)):
+                    row = self.resolve_labels(item[0], self.row_labels)
+                    col = self.resolve_labels(item[1], self.col_labels)
+                    self.slices = (slice(row, row + 1), slice(col, col + 1))
+                elif isinstance(item[0], slice) and isinstance(item[1], slice):
+                    self.slices = (Slicer.parse_slice(item[0], self.row_labels),
+                                   Slicer.parse_slice(item[1], self.col_labels))
+                elif isinstance(item[0], slice) and isinstance(item[1], (int, str)):
+                    col = self.resolve_labels(item[1], self.col_labels)
+                    if not 0 <= col < self.n_cols:
+                        raise ValueError("Column index out of range.")
+                    self.slices = (self.resolve_labels(item[0], self.row_labels), slice(col, col + 1))
+                elif isinstance(item[0], (int, str)) and isinstance(item[1], slice):
+                    row = self.resolve_labels(item[0], self.row_labels)
+                    if not 0 <= row < self.n_rows:
+                        raise ValueError("Row index out of range.")
+                    self.slices = (slice(row, row + 1), self.resolve_labels(item[1], self.col_labels))
+                else:
+                    raise TypeError("Invalid slice.")
+            else:
+                raise TypeError("Invalid slice.")
         else:
-            self.slices = self.parse_tuple(item)
+            raise TypeError("Invalid slice.")
 
     def copy(self):
         return Slicer(self.array, self.row_labels, self.col_labels, self.item)
@@ -119,71 +158,102 @@ def parse_tuple(self, item):
         if isinstance(item, tuple) and len(item) == 1:
             item = item[0]
         if isinstance(item, tuple) and len(item) == 2:
-            return (self.resolve_labels(item[0], self.row_labels),
-                    self.resolve_labels(item[1], self.col_labels))
+            if isinstance(item[0], (str, int)) and isinstance(item[1], (str, int)):
+                return self.resolve_labels(item[0], self.row_labels), self.resolve_labels(item[1], self.col_labels)
+            else:
+                raise TypeError("Invalid slice.")
         elif isinstance(item, tuple):
             raise ValueError("We only have two axes.")
         else:
-            return self.resolve_labels(item, self.row_labels)
+            raise TypeError("Invalid slice.")
 
     @staticmethod
-    def resolve_labels(item, labels: List[str]) -> Union[int, slice]:
+    def parse_slice(item, labels: List[str]) -> slice:
         """
-        Convert argument passed into __getitem__ to only contain integer indices.
+        Convert a slice to a tuple of integers.
 
         Args:
-            item: Slice(s) passed to __getitem__.
             labels: Row or column labels.
+            item: Slice to parse, i.e. slice(1, None, None) or slice('A', None, None).
 
-        Raises:
-            ValueError: Label not found, or index out of range.
-            TypeError:  Invalid type for start, stop, or step.
-                        Invalid type for item.
         """
-        if isinstance(item, int):
-            if not 1 <= item <= len(labels):
-                raise ValueError("Index out of range")
-            return item - 1  # Use one-based indexing
-        if isinstance(item, str):
-            try:
-                return labels.index(item)
-            except ValueError:
-                raise ValueError(f"Label not found: {item}")
         if isinstance(item, slice):
             start, stop, step = item.start, item.stop, item.step
             if not (step is None or isinstance(step, int)):
                 raise TypeError("Step must be None or an integer")
-            if isinstance(start, str):
-                start = labels.index(start)
-            elif isinstance(start, int):
-                if not 1 <= start <= len(labels):
-                    raise ValueError("Index out of range")
-                start -= 1  # Use one-based indexing
-            elif start is not None:
-                raise TypeError("Start value must be a label or int.")
-
-            if isinstance(stop, str):
-                stop = labels.index(stop)
-            elif isinstance(stop, int):
-                if not 1 <= stop <= len(labels):
-                    raise ValueError("Index out of range")
-                stop -= 1  # Use one-based indexing
-            elif stop is not None:
-                raise TypeError("Stop value must be a label or int.")
-
+            if start is not None:
+                if isinstance(start, str):
+                    start = Slicer.resolve_labels(start, labels)
+                elif isinstance(start, int):
+                    if not 1 <= start <= len(labels):
+                        raise ValueError("Index out of range")
+                    start -= 1
+                else:
+                    raise TypeError("Invalid type for start.")
             if stop is not None:
-                stop += 1  # We want values up to and including stop.
+                if isinstance(stop, str):
+                    if stop not in labels:
+                        raise ValueError(f"Label not found: {stop}")
+                    stop = labels.index(stop) + 1
+                elif isinstance(stop, int):
+                    if not 1 <= stop <= len(labels):
+                        raise ValueError("Index out of range")
+                    # stop is not decremented because slice() is exclusive
+                else:
+                    raise TypeError("Invalid type for stop.")
             return slice(start, stop, step)
-        raise TypeError("Invalid type for a slice.")
+        raise TypeError("Invalid slice.")
+
+    @staticmethod
+    def resolve_labels(item: str | slice, labels: List[str]) -> int | slice:
+        """
+        Convert argument passed into __getitem__ to only contain integer indices.
+
+        Args:
+            item: String or slice passed to __getitem__.
+            labels: Row or column labels.
+
+        Returns: Converted item.
+
+        """
+        if isinstance(item, int):
+            if not 1 <= item <= len(labels):
+                raise ValueError("Index out of range")
+            return item - 1
+        elif isinstance(item, str):
+            if item in labels:
+                return labels.index(item)
+            else:
+                raise ValueError(f"Label not found: {item}")
+        elif isinstance(item, slice):
+            return Slicer.parse_slice(item, labels)
 
     def get(self):
         """
         Get data pointed to by slices.
         """
         if isinstance(self.slices, list):
-            return np.array(list(map(self.array.__getitem__, self.slices)))
+            return np.array(list(map(self.array.__getitem__, self.slices))).flatten()
         return self.array.__getitem__(self.slices)
 
+    def apply(self, func):
+        """
+        Apply function to data pointed to by slices.
+
+        Args:
+            func: Function to apply.
+
+        Returns:
+            Result of function.
+
+        """
+        if isinstance(self.slices, list):
+            for elem in self.slices:
+                result = numpy.vectorize(func, cache=True)(self.array.__getitem__(elem))
+                self.array.__setitem__(elem, result)
+        else:
+            self.array.__setitem__(self.slices, numpy.vectorize(func, cache=True)(self.array.__getitem__(self.slices)))
+
     def set(self, values):
         """
         Replace data pointed to by slice(s).
@@ -194,15 +264,19 @@ def set(self, values):
             values: Data to store.
 
         Raises:
-            ValueError: Shape or size don't match.
+            ValueError: Shape or size doesn't match.
         """
-        if np.shape(values) != self.shape or np.size(values) != self.size:
-            raise ValueError("Shape or size of values doesn't match.")
-        if isinstance(self.slices, list):
+        if isinstance(values, list):
+            if len(values) != len(self.slices):
+                raise ValueError("Shape or size of values doesn't match.")
             for index, value in zip(self.slices, values):
-                self.array.__setitem__(index, value)
-        else:
+                self.array.__setitem__(index, [[value]])
+        elif isinstance(values, np.ndarray):
+            if np.shape(values) != self.shape or np.size(values) != self.size:
+                raise ValueError("Shape or size of values doesn't match.")
             self.array.__setitem__(self.slices, values)
+        else:
+            self.array.__setitem__(self.slices, [[values]])
 
     @cached_property
     def shape(self):
@@ -220,3 +294,64 @@ def size(self):
 
     def __repr__(self):
         return f"Slice: [{self.slices}]"
+
+    @staticmethod
+    def _process_sub_slice(curr_slice: slice, sub_slice: slice, labels):
+        start, stop, step = curr_slice.start, curr_slice.stop, curr_slice.step
+        if start is None:
+            start = 0
+        if stop is None:
+            stop = len(labels)
+        if step is None:
+            step = 1
+        if sub_slice.start is None:
+            sub_slice = slice(0, sub_slice.stop, sub_slice.step)
+
+        if sub_slice.start is not None and sub_slice.start < 0:
+            sub_slice = slice(stop + sub_slice.start, sub_slice.stop, sub_slice.step)
+        if sub_slice.stop is not None and sub_slice.stop < 0:
+            sub_slice = slice(sub_slice.start, stop + sub_slice.stop, sub_slice.step)
+        if sub_slice.stop is None:
+            length = float('inf')
+        else:
+            length = sub_slice.stop - sub_slice.start
+
+        start = start + sub_slice.start * step
+        stop = min(start + length * step, stop)
+        if sub_slice.step is not None:
+            step *= sub_slice.step
+
+        if start == 0:
+            start = None
+        if stop >= len(labels):
+            stop = None
+        if step == 1:
+            step = None
+        return slice(start, stop, step)
+
+    def __getitem__(self, item):
+        # negative indexing not supported
+        if isinstance(self.slices, list):
+            new_slicer = copy(self)
+            new_slicer.slices = new_slicer.slices.__getitem__(item)
+            return new_slicer
+
+        if isinstance(item, (int, slice)):
+            item = (item, slice(None))
+        if isinstance(item, tuple) and len(item) == 2:
+            if isinstance(item[0], int):
+                row = item[0]
+                item = (slice(row, row + 1), item[1])
+            if isinstance(item[1], int):
+                col = item[1]
+                item = (item[0], slice(col, col + 1))
+            if isinstance(item[0], slice) and isinstance(item[1], slice):
+                new_slicer = copy(self)
+                new_slicer.items = item
+                new_slicer.slices = (Slicer._process_sub_slice(self.slices[0], item[0], self.row_labels),
+                                     Slicer._process_sub_slice(self.slices[1], item[1], self.col_labels))
+                return new_slicer
+            else:
+                raise TypeError("Invalid slice.")
+        else:
+            raise TypeError("Invalid slice.")
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..47d0864
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,29 @@
+[build-system]
+requires      = ["setuptools>=63.0.0", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyplate-hte"
+version = "0.2"
+description = "A Python tool for designing chemistry experiments in plate format"
+readme = "README.md"
+license = {file = "LICENSE.txt"}
+authors = [
+    {name = "Eugene Kwan and James Marvin"},
+    {email = "ekwan16@gmail.com"}
+]
+classifiers = [
+    "License :: OSI Approved :: GNU General Public License v3 (GPLv3)"
+]
+dependencies = [
+    "numpy",
+    "pandas",
+    "pyyaml"
+]
+
+
+[tool.setuptools]
+packages = ["pyplate"]
+
+[project.urls]
+Repository = "https://github.com/ekwan/PyPlate"
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index 7321845..fcf65ef 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,1 +1,6 @@
-numpy>=1.23.5
\ No newline at end of file
+numpy>=1.23.5
+pytest>=7.1.2
+pytest-cov>=4.1.0
+pyyaml>=6.0.1
+pandas==1.5.3
+pdoc==14.0.0
diff --git a/solution_formulas.md b/solution_formulas.md
new file mode 100644
index 0000000..16f9085
--- /dev/null
+++ b/solution_formulas.md
@@ -0,0 +1,123 @@
+x = solute in moles, y = solvent in moles
+
+$MW_x$ = molecular weight of solute, $d_x$ = density of solute
+
+$MW_y$ = molecular weight of solvent, $d_y$ = density of solvent
+
+---
+c = $\frac{g}{g}$ = $\frac{x\cdot MW_x}{x\cdot MW_x + y\cdot MW_y}$
+
+$c\cdot x\cdot MW_x + c\cdot y\cdot MW_y = x\cdot MW_x$
+
+$c\cdot y\cdot MW_y = x\cdot MW_x (1-c)$
+
+$\frac{x}{y} = \frac{c\cdot MW_y}{(1-c)MW_x}$
+---
+c = $\frac{g}{mol}$ = $\frac{x\cdot MW_x}{x+y}$
+
+$c\cdot x + c\cdot y = x\cdot MW_x$
+
+$c\cdot y = x(MW_x-c)$
+
+$\frac{x}{y} = \frac{x}{MW_x -c}$
+---
+c = $\frac{g}{L}\cdot \frac{1 L}{1000 mL} = \frac{g}{mL}$
+
+c = $\frac{x\cdot MW_x}{x\cdot MW_x/d_x + y\cdot MW_y/d_y}$
+
+$c\cdot x\cdot MW_x/d_y + c\cdot y\cdot MW_y/d_y = x\cdot MW_x$
+
+$c\cdot y\cdot MW_y/d_y = x\cdot MW_x(1-c/d_y)$
+
+$\frac{x}{y} = \frac{c\cdot MW_y}{MW_x\cdot d_y(1-c/d_x)}$
+---
+
+c = $\frac{L}{g}\cdot \frac{1000 mL}{1 L} = \frac{mL}{g}$
+
+c = $\frac{x\cdot MW_x/d_x}{x\cdot MW_x + y\cdot MW_y}$
+
+$c\cdot x\cdot MW_x + c\cdot y\cdot MW_y = x\cdot MW_x/d_x$
+
+$c\cdot y\cdot MW_y = x\cdot MW_x(1/d_x-c)$
+
+$\frac{x}{y} = \frac{c\cdot MW_y}{MW_x(1/d_x-c)}$
+
+---
+
+c = $\frac{L}{mol}\cdot \frac{1000 mL}{1 L} = \frac{mL}{mol}$
+
+c = $\frac{x\cdot MW_x/d_x}{x + y}$
+
+$c\cdot x + c\cdot y = x\cdot MW_x/d_x$
+
+$c\cdot y = x(MW_x/d_x - c)$
+
+$\frac{x}{y} = \frac{c}{MW_x/d_x-c}$
+
+---
+
+c = $\frac{L}{L} = \frac{x\cdot MW_x/d_x}{x\cdot MW_x/d_x + y\cdot MW_y/d_y}$
+
+$c\cdot x\cdot MW_x/d_x + c\cdot y\cdot MW_y/d_y = x\cdot MW_x/d_x$
+
+$c\cdot y\cdot MW_y/d_y = (1-c)x\cdot MW_x/d_x$
+
+$\frac{x}{y} = \frac{c\cdot MW_y/d_y}{(1-c)MW_x/d_x}$
+
+---
+
+c = $\frac{mol}{g} = \frac{x}{x\cdot MW_x + y\cdot MW_y}$
+
+$x = c\cdot x\cdot MW_x + c\cdot y\cdot MW_y$
+
+$x(1-c\cdot MW_x) = c\cdot y\cdot MW_y$
+
+$\frac{x}{y} = \frac{c\cdot MW_y}{1-c\cdot MW_x}$
+
+---
+
+c = $\frac{mol}{mol} = \frac{x}{x+y}$
+
+$x = c\cdot x + c\cdot y$
+
+$c\cdot y = x(1-c)$
+
+$\frac{x}{y} = \frac{c}{1-c}$
+
+---
+
+c = $\frac{mol}{L}\cdot \frac{1 L}{1000 mL} = \frac{mol}{mL} = \frac{x}{x\cdot MW_x/d_x + y\cdot MW_y/d_y}$
+
+$c\cdot x\cdot MW_x/d_x + c\cdot y\cdot MW_y/d_y = x$
+
+$c\cdot y\cdot MW_y/d_y = x(1 - c\cdot MW_x/d_x)$
+
+$\frac{x}{y} = \frac{x\cdot MW_y/d_y}{1 - c\cdot MW_x/d_x}$
+
+---
+---
+
+#### For quantities of grams, we convert $\frac{x}{y}$ to $\frac{g}{g}$
+
+$\frac{x}{y}\cdot \frac{MW_x}{MW_y}$
+
+#### For quantities of liters, we convert $\frac{x}{y}$ to $\frac{mL}{mL}$. Since it is a ratio, it doesn't matter if we use mL or L.
+
+$\frac{x}{y}\cdot \frac{MW_x/d_x}{MW_y/d_y}$
+
+#### For quantities of moles, we do nothing since the ratio is already in $\frac{mol}{mol}$
+
+---
+
+$ratio = quantity_{solute}/quantity_{solvent}$
+
+$quantity_{solute} = ratio\cdot quantity_{solvent}$
+
+$quantity = quantity_{solute} + quantity_{solvent} = ratio\cdot quantity_{solvent} + quantity_{solvent}$
+
+$= quantity_{solvent}(1 + ratio)$
+
+$quantity_{solvent} = \frac{quantity}{1+ratio}$
+
+$quantity_{solute} = quantity - quantity_{solvent}$
+


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..8aab1a1
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,46 @@
+import pytest
+from pyplate.pyplate import Substance, Container, Plate
+
+
+@pytest.fixture
+def salt() -> Substance:
+    return Substance.solid('NaCl', 58.4428)
+
+
+@pytest.fixture
+def water() -> Substance:
+    return Substance.liquid('H2O', mol_weight=18.0153, density=1)
+
+
+@pytest.fixture
+def water_stock(water) -> Container:
+    return Container('water', initial_contents=((water, '10 mL'),))
+
+
+@pytest.fixture
+def salt_water(water, salt) -> Container:
+    return Container('salt water', initial_contents=((water, '100 mL'), (salt, '50 mmol')))
+
+
+@pytest.fixture
+def dmso() -> Substance:
+    return Substance.liquid('DMSO', 78.13, 1.1004)
+
+
+@pytest.fixture
+def sodium_sulfate() -> Substance:
+    return Substance.solid('Sodium sulfate', 142.04)
+
+
+@pytest.fixture
+def triethylamine() -> Substance:
+    return Substance.liquid("triethylamine", mol_weight=101.19, density=0.726)
+
+@pytest.fixture
+def lipase() -> Substance:
+    return Substance.enzyme('lipase')
+
+
+@pytest.fixture
+def empty_plate() -> Plate:
+    return Plate('plate', '200 uL')  # 100 uL
diff --git a/tests/test_Container.py b/tests/test_Container.py
index 7e24424..b4503d1 100644
--- a/tests/test_Container.py
+++ b/tests/test_Container.py
@@ -1,101 +1,91 @@
 import pytest
-from pyplate.pyplate import Substance, Container, Unit
+from pyplate import Container
+from pyplate.pyplate import config, Unit
 
 
-@pytest.fixture
-def water() -> Substance:
-    return Substance.liquid('H2O', mol_weight=18.0153, density=1)
-
-
-@pytest.fixture
-def salt() -> Substance:
-    return Substance.solid('NaCl', 58.4428)
-
-
-@pytest.fixture
-def water_stock(water) -> Container:
-    return Container('water', initial_contents=((water, '10 mL'),))
-
-
-@pytest.fixture
-def salt_water(water, salt) -> Container:
-    return Container('salt water', initial_contents=((water, '100 mL'), (salt, '50 mol')))
+def test_make_Container(water, salt):
+    """
 
+    Tests that `Container` constructor checks the type of all arguments.
 
-def test_make_Container(water, salt):
+    """
     # Argument types checked
     with pytest.raises(TypeError, match="Name must be a str"):
         Container(1)
     with pytest.raises(ValueError, match="Name must not be empty"):
         Container('')
-    with pytest.raises(TypeError, match="Maximum volume must be a float"):
-        Container('container', 'max_volume')
+    with pytest.raises(TypeError, match='Maximum volume must be a str'):
+        Container('container', 1)
+    with pytest.raises(ValueError, match="Value is not a valid float"):
+        Container('container', 'max_volume L')
     with pytest.raises(ValueError, match="Maximum volume must be positive"):
-        Container('container', -1)
+        Container('container', '-1 L')
     with pytest.raises(ValueError, match="Maximum volume must be positive"):
-        Container('container', 0)
+        Container('container', '0 L')
     with pytest.raises(TypeError, match="Initial contents must be iterable"):
-        Container('container', 1, 1)
+        Container('container', '1 L', 1)
     with pytest.raises(TypeError, match="Element in initial_contents must be"):
-        Container('container', 1, [1])
+        Container('container', '1 L', [1])
     with pytest.raises(TypeError, match="Element in initial_contents must be"):
-        Container('container', 1, [water, salt])
+        Container('container', '1 L', [water, salt])
 
 
-def test_Container_add(water):
-    container = Container('container')
-    with pytest.raises(TypeError, match='Source must be a Substance'):
-        Container.add(1, container, '10 mL')
-    with pytest.raises(TypeError, match='You can only use Container.add to add to a Container'):
-        Container.add(water, 1, '10 mL')
-    with pytest.raises(TypeError, match='Quantity must be a str'):
-        Container.add(water, container, 1)
-
-    result = Container.add(water, container, '10 mL')
-    # 10 mL of water should have been added.
-    assert isinstance(result, Container)
-    assert result.volume == Unit.convert_to_storage(10, 'mL')
-    assert water in result.contents and result.contents[water] == Unit.convert_to_storage(10, 'mL')
-    # Original should be unchanged
-    assert container.volume == 0
-    assert water not in container.contents
+def test_Container_transfer(water, salt, water_stock, salt_water):
+    """
 
+    Tests if transfers between `Container`s works properly.
 
-def test_Container_transfer(water, salt, water_stock, salt_water):
+    """
+    # Argument types checked
     with pytest.raises(TypeError, match='into a Container'):
         Container.transfer(1, 1, '10 mL')
     with pytest.raises(TypeError, match='Invalid source type'):
         Container.transfer(1, water_stock, '10 mL')
-    with pytest.raises(TypeError, match='Amount must be a str'):
+    with pytest.raises(TypeError, match='Quantity must be a str'):
         Container.transfer(salt_water, water_stock, 1)
 
     initial_hashes = hash(water_stock), hash(salt_water)
-    # water_stock is 10 mL, salt_water is 100 mL and 50 mol
-    container1, container2 = Container.transfer(salt_water, water_stock, '10 mL')
+    # water_stock is 10 mL, salt_water is 100 mL and 50 mmol
+    salt_water_volume = Unit.convert_from_storage(salt_water.volume, 'mL')
+    container1, container2 = Container.transfer(salt_water, water_stock, f"{salt_water_volume*0.1} mL")
     # 10 mL of water and 5 mol of salt should have been transferred
-    assert container1.volume == Unit.convert_to_storage(90, 'mL')
-    assert container1.contents[water] == Unit.convert_to_storage(90, 'mL')
-    assert container1.contents[salt] == Unit.convert_to_storage(45, 'mol')
-    assert container2.volume == Unit.convert_to_storage(20, 'mL')
-    assert salt in container2.contents and container2.contents[salt] == Unit.convert_to_storage(5, 'mol')
-    assert container2.contents[water] == Unit.convert_to_storage(20, 'mL')
+    assert container1.volume == Unit.convert(water, '90 mL', config.volume_prefix)\
+           + Unit.convert(salt, '45 mmol', config.volume_prefix)
+    assert container1.contents[water] == Unit.convert(water, '90 mL', config.moles_prefix)
+    assert container1.contents[salt] == Unit.convert(salt, '45 mmol', config.moles_prefix)
+    assert container2.volume == Unit.convert(water, '20 mL', config.volume_prefix)\
+           + Unit.convert(salt, '5 mmol', config.volume_prefix)
+    assert salt in container2.contents and container2.contents[salt] == \
+           Unit.convert(salt, '5 mmol', config.moles_prefix)
+    assert container2.contents[water] == pytest.approx(Unit.convert(water, '20 mL', config.moles_prefix))
 
     # Original containers should be unchanged.
     assert initial_hashes == (hash(water_stock), hash(salt_water))
 
+    salt_stock = Container('salt stock', initial_contents=[(salt, '10 g')])
+    container1, container2 = Container.transfer(salt_stock, salt_water, '1 g')
+    assert container2.contents[salt] == \
+           pytest.approx(salt_water.contents[salt] + Unit.convert(salt, '1 g', config.moles_prefix))
+
 
 def test_create_stock_solution(water, salt, salt_water):
+    """
+
+    Tests if transfers between `Container`s works properly.
+
+    """
+    # Argument types checked
     with pytest.raises(TypeError, match='Solute must be a Substance'):
-        Container.create_stock_solution('salt', 0.5, water, 100)
-    with pytest.raises(TypeError, match='Concentration must be a float'):
-        Container.create_stock_solution(salt, '0.5', water, 100)
+        Container.create_solution('salt', water, concentration='0.5 M', total_quantity='100 mL')
+    with pytest.raises(TypeError, match='Concentration must be a str'):
+        Container.create_solution(salt, water, concentration=.5, total_quantity='100 mL')
     with pytest.raises(TypeError, match='Solvent must be a Substance'):
-        Container.create_stock_solution(salt, 0.5, 'water', 100)
-    with pytest.raises(TypeError, match='Volume must be a float'):
-        Container.create_stock_solution(salt, 0.5, water, '100')
+        Container.create_solution(salt, 'water', concentration='0.5 M', total_quantity='100 mL')
+    with pytest.raises(TypeError, match='Total quantity must be a str'):
+        Container.create_solution(salt, water, concentration='0.5 M', total_quantity=100.0)
 
-    stock = Container.create_stock_solution(salt, 0.5, water, 100)
+    stock = Container.create_solution(salt, water, concentration='0.5 M', total_quantity='100 mL')
+    # stock should have 100 mL of water and 50 mmol of salt
     assert water in stock.contents and salt in stock.contents
-    assert stock.contents[water] == salt_water.contents[water]
-    assert stock.contents[salt] == salt_water.contents[salt]
-
+    assert stock.volume == Unit.convert_to_storage(100, 'mL')
+    assert stock.contents[salt] == pytest.approx(salt_water.contents[salt])
diff --git a/tests/test_Experiment.py b/tests/test_Experiment.py
new file mode 100644
index 0000000..06916cb
--- /dev/null
+++ b/tests/test_Experiment.py
@@ -0,0 +1,62 @@
+import pytest
+from pyplate.experiment_design import Factor, Experiment, ExperimentalSpace
+from pyplate.pyplate import Container
+
+
+@pytest.fixture
+def example_experiment():
+    factors = {"Factor1": "Value1", "Factor2": "Value2"}
+    return Experiment(factors, 1, 0)
+
+
+# Tests for the Experiment class
+def test_experiment_creation(example_experiment):
+    assert example_experiment.experiment_id == 1
+    assert example_experiment.replicate_idx == 0
+    assert example_experiment.well is None
+
+
+def test_experiment_mapping_well(example_experiment):
+    well = Container("96-well")
+    example_experiment.map_well(well)
+    assert example_experiment.well == well
+
+
+def test_experiment_repr(example_experiment):
+    expected_repr = "Experiment({'Factor1': 'Value1', 'Factor2': 'Value2'}, 1, 0, None)"
+    assert repr(example_experiment) == expected_repr
+
+
+def test_experiment_str(example_experiment):
+    expected_str = "Experiment: {'Factor1': 'Value1', 'Factor2': 'Value2'} with experiment_id 1 and replicate_idx 0"
+    assert str(example_experiment) == expected_str
+
+
+def test_experiment_getitem(example_experiment):
+    assert example_experiment["Factor1"] == "Value1"
+
+
+def test_experiment_setitem(example_experiment):
+    example_experiment["Factor1"] = "ModifiedValue"
+    assert example_experiment["Factor1"] == "ModifiedValue"
+
+
+def test_experiment_contains(example_experiment):
+    assert "Factor1" in example_experiment
+    assert "NonexistentFactor" not in example_experiment
+
+
+def test_experiment_iter(example_experiment):
+    factors_iterator = iter(example_experiment)
+
+    # Check that the iterator returns the factors in the expected order
+    assert next(factors_iterator) == "Factor1"
+    assert next(factors_iterator) == "Factor2"
+
+    # Check that trying to get the next element raises StopIteration
+    with pytest.raises(StopIteration):
+        next(factors_iterator)
+
+
+def test_experiment_len(example_experiment):
+    assert len(example_experiment) == 2
diff --git a/tests/test_ExperimentalSpace.py b/tests/test_ExperimentalSpace.py
new file mode 100644
index 0000000..16ea98a
--- /dev/null
+++ b/tests/test_ExperimentalSpace.py
@@ -0,0 +1,106 @@
+import pytest
+from pyplate.experiment_design import Factor, Experiment, ExperimentalSpace
+import re
+import collections
+
+
+@pytest.fixture
+def example_factor():
+    return Factor("ExampleFactor", ["Value1", "Value2"], lambda x, y: True)
+
+
+@pytest.fixture
+def another_factor():
+    return Factor("AnotherFactor", ["ValueA", "ValueB"], lambda x, y: True)
+
+
+@pytest.fixture
+def example_experiment(example_factor):
+    factors = {"ExampleFactor": "Value1"}
+    return Experiment(factors, 1, 0)
+
+
+@pytest.fixture
+def example_experimental_space(example_factor):
+    return ExperimentalSpace(set(), lambda: 1)
+
+
+@pytest.fixture
+def cookie_experimental_space():
+    f1 = Factor(name="Baking Time", possible_values=[10, 20, 30], verifier=lambda: True)
+    f2 = Factor(name="Baking Temperature", possible_values=[300, 350, 400], verifier=lambda: True)
+    f3 = Factor(name="Flavor", possible_values=["Chocolate Chip", "Oatmeal Raisin", "Peanut Butter"],
+                verifier=lambda: True)
+
+    return ExperimentalSpace(factors={f1, f2, f3}, experiment_id_generator=lambda: 1)
+
+
+# Tests for the ExperimentalSpace class
+def test_experimental_space_register_factor_verification(example_experimental_space, example_factor):
+    example_experimental_space.register_factor(example_factor)
+
+    # Verify that registering a factor with the same name raises ValueError
+    with pytest.raises(Exception, match=f"Factor {example_factor.name} already exists in experimental space"):
+        example_experimental_space.register_factor(example_factor)
+
+
+def test_experimental_space_add_experiment_verification(example_experimental_space, example_experiment, another_factor):
+    example_experimental_space.register_factor(another_factor)
+
+    # Verify that adding an experiment with factors not matching experimental space factors raises ValueError
+    invalid_experiment = Experiment({"InvalidFactor": "Value1"}, 2, 0)
+    with pytest.raises(ValueError, match=re.escape(f"Experiment factors {invalid_experiment.factors.keys()} "
+                                                   f"do not match experimental space factors {example_experimental_space.factors}")):
+        example_experimental_space.add_experiment(invalid_experiment)
+
+    # Verify that adding an experiment with a factor value not in possible values raises ValueError
+    invalid_experiment = Experiment({"AnotherFactor": "InvalidValue"}, 3, 0)
+    with pytest.raises(ValueError, match=re.escape("Experiment factor AnotherFactor value InvalidValue "
+                                                   f"not in possible values ['ValueA', 'ValueB']")):
+        example_experimental_space.add_experiment(invalid_experiment)
+
+
+def test_experimental_space_generate_experiments(cookie_experimental_space):
+    expected_exps = \
+        {(300, 'Chocolate Chip'): [
+            Experiment({'Baking Time': 10, 'Baking Temperature': 300, 'Flavor': 'Chocolate Chip'}, 1, 1, None),
+            Experiment({'Baking Time': 10, 'Baking Temperature': 300, 'Flavor': 'Chocolate Chip'}, 1, 2, None),
+            Experiment({'Baking Time': 30, 'Baking Temperature': 300, 'Flavor': 'Chocolate Chip'}, 1, 1, None),
+            Experiment({'Baking Time': 30, 'Baking Temperature': 300, 'Flavor': 'Chocolate Chip'}, 1, 2, None)],
+         (300, 'Oatmeal Raisin'): [
+             Experiment({'Baking Time': 10, 'Baking Temperature': 300, 'Flavor': 'Oatmeal Raisin'}, 1, 1, None),
+             Experiment({'Baking Time': 10, 'Baking Temperature': 300, 'Flavor': 'Oatmeal Raisin'}, 1, 2, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 300, 'Flavor': 'Oatmeal Raisin'}, 1, 1, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 300, 'Flavor': 'Oatmeal Raisin'}, 1, 2, None)],
+         (300, 'Peanut Butter'): [
+             Experiment({'Baking Time': 10, 'Baking Temperature': 300, 'Flavor': 'Peanut Butter'}, 1, 1, None),
+             Experiment({'Baking Time': 10, 'Baking Temperature': 300, 'Flavor': 'Peanut Butter'}, 1, 2, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 300, 'Flavor': 'Peanut Butter'}, 1, 1, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 300, 'Flavor': 'Peanut Butter'}, 1, 2, None)],
+         (350, 'Chocolate Chip'): [
+             Experiment({'Baking Time': 10, 'Baking Temperature': 350, 'Flavor': 'Chocolate Chip'}, 1, 1, None),
+             Experiment({'Baking Time': 10, 'Baking Temperature': 350, 'Flavor': 'Chocolate Chip'}, 1, 2, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 350, 'Flavor': 'Chocolate Chip'}, 1, 1, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 350, 'Flavor': 'Chocolate Chip'}, 1, 2, None)],
+         (350, 'Oatmeal Raisin'): [
+             Experiment({'Baking Time': 10, 'Baking Temperature': 350, 'Flavor': 'Oatmeal Raisin'}, 1, 1, None),
+             Experiment({'Baking Time': 10, 'Baking Temperature': 350, 'Flavor': 'Oatmeal Raisin'}, 1, 2, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 350, 'Flavor': 'Oatmeal Raisin'}, 1, 1, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 350, 'Flavor': 'Oatmeal Raisin'}, 1, 2, None)],
+         (350, 'Peanut Butter'): [
+             Experiment({'Baking Time': 10, 'Baking Temperature': 350, 'Flavor': 'Peanut Butter'}, 1, 1, None),
+             Experiment({'Baking Time': 10, 'Baking Temperature': 350, 'Flavor': 'Peanut Butter'}, 1, 2, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 350, 'Flavor': 'Peanut Butter'}, 1, 1, None),
+             Experiment({'Baking Time': 30, 'Baking Temperature': 350, 'Flavor': 'Peanut Butter'}, 1, 2, None)]}
+
+    generated_exps = cookie_experimental_space.generate_experiments(
+        factors={
+            "Baking Time": [10, 30],
+            "Baking Temperature": [300, 350],
+            "Flavor": "all"
+        },
+        n_replicates=2,
+        blocking_factors=["Baking Temperature", "Flavor"]
+    )
+
+    assert generated_exps == expected_exps
diff --git a/tests/test_Factor.py b/tests/test_Factor.py
new file mode 100644
index 0000000..b5df713
--- /dev/null
+++ b/tests/test_Factor.py
@@ -0,0 +1,32 @@
+import pytest
+from pyplate.experiment_design import Factor
+
+
+def dummy_verifier():
+    return True
+
+
+@pytest.fixture
+def example_factor():
+    return Factor("ExampleFactor", ["Value1", "Value2"], dummy_verifier)
+
+
+def test_factor_creation(example_factor):
+    assert example_factor.name == "ExampleFactor"
+    assert example_factor.possible_values == ["Value1", "Value2"]
+
+
+def test_factor_representation(example_factor):
+    assert str(example_factor) == "Factor: ExampleFactor with possible values ['Value1', 'Value2']"
+
+
+def test_factor_equality(example_factor):
+    # Create a copy of the example_factor
+    copied_factor = Factor("ExampleFactor", ["Value1", "Value2"], dummy_verifier)
+
+    # Check that the __eq__ method correctly identifies the factors as equal
+    assert example_factor == copied_factor
+
+    # Change a property and check that the factors are not equal
+    copied_factor.possible_values = ["Value1", "Value2", "Value3"]
+    assert example_factor != copied_factor
diff --git a/tests/test_Plate.py b/tests/test_Plate.py
new file mode 100644
index 0000000..758d7aa
--- /dev/null
+++ b/tests/test_Plate.py
@@ -0,0 +1,113 @@
+import pytest
+import numpy
+from pyplate.pyplate import Plate, Unit, config, Container
+
+
+def test_make_Plate():
+    """
+
+    Test that all argument types are checked in Plate constructor.
+
+    """
+    with pytest.raises(ValueError, match='invalid plate name'):
+        Plate(1, '10 uL')
+    with pytest.raises(ValueError, match='invalid plate name'):
+        Plate('', '10 uL')
+    with pytest.raises(TypeError, match='Maximum volume must be a str'):
+        Plate('plate', 10)
+    with pytest.raises(ValueError, match='invalid plate make'):
+        Plate('plate', '10 uL', make=1)
+    with pytest.raises(ValueError, match='invalid plate make'):
+        Plate('plate', '10 uL', make='')
+    with pytest.raises(ValueError, match='rows must be int or list'):
+        Plate('plate', '10 uL', rows='8')
+    with pytest.raises(ValueError, match='illegal number of rows'):
+        Plate('plate', '10 uL', rows=0)
+    with pytest.raises(ValueError, match='columns must be int or list'):
+        Plate('plate', '10 uL', columns='8')
+    with pytest.raises(ValueError, match='illegal number of columns'):
+        Plate('plate', '10 uL', columns=0)
+
+    with pytest.raises(ValueError, match='must have at least one row'):
+        Plate('plate', '10 uL', rows=[])
+    with pytest.raises(ValueError, match='row names must be strings'):
+        Plate('plate', '10 uL', rows=[1])
+    with pytest.raises(ValueError, match='duplicate row names found'):
+        Plate('plate', '10 uL', rows=['a', 'a'])
+
+    with pytest.raises(ValueError, match='must have at least one column'):
+        Plate('plate', '10 uL', columns=[])
+    with pytest.raises(ValueError, match='column names must be strings'):
+        Plate('plate', '10 uL', columns=[1])
+    with pytest.raises(ValueError, match='duplicate column names found'):
+        Plate('plate', '10 uL', columns=['a', 'a'])
+
+
+def test_volume_and_volumes(salt, water, dmso, empty_plate):
+    """
+
+    Test volume() and volumes() for a plate.
+
+    """
+    epsilon = 1e-3
+    with pytest.raises(TypeError, match="Substance is not a valid type"):
+        empty_plate.volumes('1')
+
+    zeros = numpy.zeros(empty_plate.wells.shape)
+    uL = numpy.ones(empty_plate.wells.shape)
+
+    salt_volume = round(Unit.convert(salt, '5 umol', 'uL'), 20)
+    assert empty_plate.volume() == 0
+    assert numpy.array_equal(empty_plate.volumes(), zeros)
+    assert numpy.array_equal(empty_plate.volumes(water), zeros)
+
+    salt_container = Container('salt', initial_contents=((salt, '1 mol'),))
+    salt_container, new_plate = Plate.transfer(salt_container, empty_plate, '5 umol')
+    assert new_plate.volume() == pytest.approx((salt_volume * uL).sum(), epsilon)
+    assert numpy.allclose(new_plate.volumes(), salt_volume * uL)
+    assert numpy.array_equal(new_plate.volumes(water), zeros)
+    assert numpy.allclose(new_plate.volumes(salt), salt_volume * uL, atol=epsilon)
+
+    water_container = Container('water', initial_contents=((water, '1 L'),))
+    water_container, new_plate = Plate.transfer(water_container, new_plate, '50 uL')
+    assert new_plate.volume() == pytest.approx((salt_volume * uL).sum() + (50 * uL).sum(), epsilon)
+    assert numpy.allclose(new_plate.volumes(), (salt_volume + 50) * uL)
+    assert numpy.allclose(new_plate.volumes(water), 50 * uL)
+    assert numpy.allclose(new_plate.volumes(salt), salt_volume * uL, atol=1e-3)
+
+    dmso_container = Container('dmso', initial_contents=((dmso, '1 L'), ))
+    dmso_container, new_plate = Plate.transfer(dmso_container, new_plate, '25 uL')
+    assert new_plate.volume() == pytest.approx((salt_volume * uL).sum() + (75 * uL).sum(), epsilon)
+    assert numpy.allclose(new_plate.volumes(), (salt_volume + 75) * uL)
+    assert numpy.allclose(new_plate.volumes(water), 50 * uL)
+    assert numpy.allclose(new_plate.volumes(dmso), 25 * uL)
+    assert numpy.allclose(new_plate.volumes(salt), salt_volume * uL, atol=epsilon)
+    assert numpy.allclose(new_plate.volumes(water, unit='mL'), 0.05 * uL)
+
+
+def test_moles(salt, water, empty_plate):
+    """
+
+    Test moles() for a plate.
+
+    """
+    with pytest.raises(TypeError, match="Substance is not a valid type"):
+        empty_plate.moles('1')
+
+    zeros = numpy.zeros(empty_plate.wells.shape)
+    ones = numpy.ones(empty_plate.wells.shape)
+    assert numpy.array_equal(empty_plate.moles(salt), zeros)
+
+    salt_container = Container('salt', initial_contents=((salt, '1 mol'),))
+    salt_container, new_plate = Plate.transfer(salt_container, empty_plate, '5 umol')
+    assert numpy.array_equal(new_plate.moles(salt, 'umol'), 5 * ones)
+    assert numpy.array_equal(new_plate.moles(water), zeros)
+
+    water_container = Container('water', initial_contents=((water, '1 L'),))
+    water_container, new_plate = Plate.transfer(water_container, new_plate, '1 mmol')
+    assert numpy.array_equal(new_plate.moles(salt, 'umol'), 5 * ones)
+    assert numpy.array_equal(new_plate.moles(water), 0.001 * ones)
+
+    water_container = Container('water', initial_contents=((water, '1 L'),))
+    water_container, new_plate = Plate.transfer(water_container, new_plate, '5 uL')
+    new_plate.moles(water)
diff --git a/tests/test_Recipe.py b/tests/test_Recipe.py
new file mode 100644
index 0000000..7ab90e2
--- /dev/null
+++ b/tests/test_Recipe.py
@@ -0,0 +1,137 @@
+import pytest
+from pyplate import Plate, Container, Recipe
+from pyplate.pyplate import Unit, config
+from tests.conftest import salt_water
+
+
+def test_locked():
+    """
+
+    Ensures once a recipe is baked, it cannot be baked again.
+
+    """
+    recipe = Recipe()
+    recipe.bake()
+    with pytest.raises(RuntimeError, match="Recipe has already been baked"):
+        recipe.bake()
+
+
+def test_uses(salt_water):
+    """
+
+    Test uses() for a Recipe.
+    Ensures the object is added and used.
+    uses() cannot be called once the recipe is baked.
+
+    """
+    recipe = Recipe()
+    recipe.uses(salt_water)
+    index = recipe.indexes.get(salt_water, None)
+    assert index is not None
+    assert recipe.results[index] == salt_water
+
+    with pytest.raises(ValueError, match="Something declared as used wasn't used"):
+        _ = recipe.bake()
+
+    recipe = Recipe()
+    recipe.bake()
+    with pytest.raises(RuntimeError, match="This recipe is locked"):
+        recipe.uses(salt_water)
+
+
+def test_transfer(salt_water):
+    """
+
+    Test transfers within a Recipe.
+
+    """
+    container = Container('container')
+    recipe = Recipe()
+    recipe.uses(salt_water, container)
+    # Argument types checked
+    with pytest.raises(TypeError, match="Invalid source type"):
+        recipe.transfer('1', container, '10 mL')
+    with pytest.raises(TypeError, match="Invalid destination type"):
+        recipe.transfer(salt_water, '1', '10 mL')
+    with pytest.raises(TypeError, match="Volume must be a str"):
+        recipe.transfer(salt_water, container, 1)
+    recipe.transfer(salt_water, container, '10 mL')
+    # Container should still be empty until we bake.
+    assert container.volume == 0
+    new_salt_water, container = recipe.bake()
+    assert container.volume == Unit.convert_to_storage(10, 'mL')
+    assert salt_water.volume - new_salt_water.volume == Unit.convert_to_storage(10, 'mL')
+
+    # transfer() cannot ba called once the recipe is baked.
+    with pytest.raises(RuntimeError, match="This recipe is locked"):
+        recipe.transfer(salt_water, container, '10 mL')
+    recipe = Recipe()
+    plate = Plate('plate', max_volume_per_well='1 mL')
+    recipe.uses(salt_water, plate)
+    recipe.transfer(salt_water, plate, '5 uL')
+    recipe.transfer(salt_water, plate[1, 1], '1 uL')
+    new_salt_water, new_plate = recipe.bake()
+    # 5 uL transferred to 96 wells and 1 uL to one well.
+    assert salt_water.volume - new_salt_water.volume == pytest.approx(Unit.convert_to_storage(5 * 96 + 1, 'uL'))
+    volumes = new_plate.volumes(unit='uL')
+    # plate[1, 1] should have 6 uL, all others 5 uL
+    assert volumes[0, 0] == 6
+    assert volumes[0, 1] == 5
+
+    container = Container('container')
+    recipe = Recipe()
+    recipe.uses(container, new_plate)
+    recipe.transfer(new_plate[1, 1], container, '2 uL')
+    recipe.transfer(new_plate[1:2, 1:2], container, '1 uL')
+    container, new_plate2 = recipe.bake()
+    # 2 uL from plate[1, 1] and 1 uL from four wells.
+    assert container.volume == Unit.convert_to_storage(2 + 4, 'uL')
+    assert new_plate.volume('uL') - new_plate2.volume('uL') == 2 + 4
+
+
+def test_create_container(water, salt):
+    """
+
+    Test create_container within a Recipe.
+
+    """
+    recipe = Recipe()
+    # Argument types checked
+    with pytest.raises(TypeError, match='Name must be a str'):
+        recipe.create_container(1)
+    with pytest.raises(TypeError, match='Maximum volume must be a str'):
+        recipe.create_container('container', max_volume=1.0)
+    with pytest.raises(TypeError, match='Initial contents must be iterable'):
+        recipe.create_container('container', '10 mL', 1)
+    with pytest.raises(TypeError, match='Initial contents must be iterable'):
+        recipe.create_container('container', '10 mL', 1)
+    for contents in [1, (1,), (salt,), (salt, 1, 1)]:
+        with pytest.raises(TypeError, match='Elements of initial_contents must be of the form'):
+            recipe.create_container('container', '10 mL', [contents])
+    with pytest.raises(TypeError, match='Containers can only be created from substances'):
+        recipe.create_container('container', '10 mL', [(1, 1)])
+    with pytest.raises(TypeError, match='Quantity must be a str'):
+        recipe.create_container('container', '10 mL', [(salt, 1)])
+
+    recipe = Recipe()
+    recipe.bake()
+    # cannot create_container once the recipe is baked.
+    with pytest.raises(RuntimeError, match="This recipe is locked"):
+        recipe.create_container('container')
+
+    recipe = Recipe()
+    # Name set correctly
+    container = recipe.create_container('container')
+    assert container.name == 'container'
+
+    # Maximum volume set correctly
+    container2 = recipe.create_container('container2', max_volume='10 mL')
+    assert Unit.convert_to_storage(10, 'mL') == container2.max_volume
+
+    # Initial contents added correctly only at bake
+    container3 = recipe.create_container('container3', initial_contents=[(water, '10 mL'), (salt, '5 mmol')])
+    recipe.used = set(recipe.indexes.values())
+    assert container3.volume == 0
+    container, container2, container3, = recipe.bake()
+    assert Unit.convert_to_storage(10, 'mL') + Unit.convert(salt, '5 mmol', config.volume_prefix) == container3.volume
+    assert container3.contents.get(salt, None) == Unit.convert_to_storage(5, 'mmol')
diff --git a/tests/test_Slicer.py b/tests/test_Slicer.py
index 8af59ad..c10b16e 100644
--- a/tests/test_Slicer.py
+++ b/tests/test_Slicer.py
@@ -1,6 +1,6 @@
 import pytest
-from pyplate.slicer import Slicer
 import numpy as np
+from pyplate.slicer import Slicer
 
 
 @pytest.fixture
@@ -8,18 +8,30 @@ def array():
     return np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
 
 
+@pytest.fixture
+def row_labels():
+    return ['A', 'B', 'C']
+
+
+@pytest.fixture
+def col_labels():
+    return ['1', '2', '3', '4']
+
+
 def test_resolve_labels_zero():
-    # Index myst be in range 1..n
+    # Index must be in range 1..n
     with pytest.raises(ValueError):
         Slicer.resolve_labels(0, [])
 
 
 def test_resolve_labels_high():
+    # Index must be in range 1..n
     with pytest.raises(ValueError):
         Slicer.resolve_labels(100, ['1'])
 
 
 def test_resolve_labels_high_slices():
+    # Slice indexes must be in range 1..n
     with pytest.raises(ValueError):
         Slicer.resolve_labels(slice(100, None, None), ['1'])
     with pytest.raises(ValueError):
@@ -38,30 +50,145 @@ def test_resolve_labels_invalid_step():
         Slicer.resolve_labels(slice(None, None, 'label'), [])
 
 
-def test_resolve_labels_typerror():
-    # item must be an int, str, or slice
-    with pytest.raises(TypeError):
-        Slicer.resolve_labels(1.0, [])
-
-
 def test_Slicer_arraytype():
+    # array_obj must be a numpy.ndarray
     with pytest.raises(TypeError, match="must be a numpy.ndarray"):
         Slicer(None, None, None, None)
 
 
 def test_Slicer_row_labels(array):
+    # row_labels must be a list or tuple
     with pytest.raises(TypeError, match="row_labels"):
         Slicer(array, None, None, None)
+    # row_labels must not be empty
     with pytest.raises(TypeError, match="row_labels"):
         Slicer(array, (), None, None)
-    with pytest.raises(TypeError, match="row_labels"):
-        Slicer(array, [1], None, None)
 
 
 def test_Slicer_col_labels(array):
-    with pytest.raises(TypeError, match="row_labels"):
-        Slicer(array, None, None, None)
-    with pytest.raises(TypeError, match="row_labels"):
-        Slicer(array, (), None, None)
-    with pytest.raises(TypeError, match="row_labels"):
-        Slicer(array, [1], None, None)
+    # col_labels must be a list or tuple
+    with pytest.raises(TypeError, match="col_labels"):
+        Slicer(array, [], None, None)
+    # col_labels must not be empty
+    with pytest.raises(TypeError, match="col_labels"):
+        Slicer(array, ['1'], (), None)
+
+
+def test_single_string(array, row_labels, col_labels):
+    # Test a single string
+    s = Slicer(array, row_labels, col_labels, 'A:1')
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(0, 1), slice(0, 1))
+    assert np.array_equal(s.get(), np.array([[array[0, 0]]]))
+
+
+def test_single_tuple(array, row_labels, col_labels):
+    # Test a single tuple
+    s = Slicer(array, row_labels, col_labels, ('A', '1'))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(0, 1), slice(0, 1))
+    assert np.array_equal(s.get(), np.array([[array[0, 0]]]))
+
+
+def test_list_of_strings(array, row_labels, col_labels):
+    # Test a list of strings
+    s = Slicer(array, row_labels, col_labels, ['A:1', 'B:2'])
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == [(slice(0, 1), slice(0, 1)), (slice(1, 2), slice(1, 2))]
+    assert np.array_equal(s.get(), np.array([array[0, 0], array[1, 1]]))
+
+
+def test_single_label(array, row_labels, col_labels):
+    # Test a single label
+    s = Slicer(array, row_labels, col_labels, 'A')
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(0, 1), slice(None))
+    assert np.array_equal(s.get(), array[:1])
+
+
+def test_single_int(array, row_labels, col_labels):
+    # Test a single integer
+    s = Slicer(array, row_labels, col_labels, 1)
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(0, 1), slice(None))
+    assert np.array_equal(s.get(), array[:1])
+
+
+def test_single_label_and_slice(array, row_labels, col_labels):
+    # Test a single label and a slice
+    s = Slicer(array, row_labels, col_labels, ('A', slice(None)))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(0, 1), slice(None))
+    assert np.array_equal(s.get(), array[:1, :])
+
+
+def test_single_int_and_slice(array, row_labels, col_labels):
+    # Test a single label and a slice
+    s = Slicer(array, row_labels, col_labels, (1, slice(None)))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(0, 1), slice(None))
+    assert np.array_equal(s.get(), array[:1, :])
+
+
+def test_slice_and_single_label(array, row_labels, col_labels):
+    # Test a single label and a slice
+    s = Slicer(array, row_labels, col_labels, (slice(None), '1'))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(None), slice(0, 1))
+    assert np.array_equal(s.get(), array[:, :1])
+
+
+def test_slice_and_single_int(array, row_labels, col_labels):
+    # Test a single label and a slice
+    s = Slicer(array, row_labels, col_labels, (slice(None), 1))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(None), slice(0, 1))
+    assert np.array_equal(s.get(), array[:, :1])
+
+
+def test_test_two_slices(array, row_labels, col_labels):
+    # Test two slices
+    s = Slicer(array, row_labels, col_labels, (slice(None), slice(None)))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == np.s_[:, :]
+    assert np.array_equal(s.get(), array[:, :])
+
+
+def test_stop_label(array, row_labels, col_labels):
+    # Test a stop label
+    s = Slicer(array, row_labels, col_labels, (slice(None, 'B'), slice(None)))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(None, 2), slice(None))
+    assert np.array_equal(s.get(), array[:2, ])
+
+
+def test_stop_int(array, row_labels, col_labels):
+    # Test a stop label
+    s = Slicer(array, row_labels, col_labels, (slice(None, 2), slice(None)))
+    assert s.array is array
+    assert s.row_labels is row_labels
+    assert s.col_labels is col_labels
+    assert s.slices == (slice(None, 2), slice(None))
+    assert np.array_equal(s.get(), array[:2, ])
diff --git a/tests/test_Substance.py b/tests/test_Substance.py
index 2fe23a0..02ce67a 100644
--- a/tests/test_Substance.py
+++ b/tests/test_Substance.py
@@ -2,32 +2,12 @@
 from pyplate.pyplate import Substance
 
 
-@pytest.fixture
-def water() -> Substance:
-    return Substance.liquid('H2O', mol_weight=18.0153, density=1)
-
-
-@pytest.fixture
-def salt() -> Substance:
-    return Substance.solid('NaCl', 58.4428)
-
-
-@pytest.fixture
-def dmso() -> Substance:
-    return Substance.liquid('DMSO', 78.13, 1.1004)
-
-
-@pytest.fixture
-def sodium_sulfate() -> Substance:
-    return Substance.solid('Sodium sulfate', 142.04)
-
-
-@pytest.fixture
-def lipase() -> Substance:
-    return Substance.enzyme('lipase')
+def test_make_solid():
+    """
 
+    Tests creating a solid `Substance`.
 
-def test_make_solid():
+    """
     # Argument types checked
     with pytest.raises(TypeError, match="Name must be a str"):
         Substance.solid(1, 1)
@@ -43,11 +23,21 @@ def test_make_solid():
 
 
 def test_solid(salt):
+    """
+
+    Tests that members of a solid `Substance` are correct.
+
+    """
     assert salt.name == 'NaCl'
     assert salt.mol_weight == 58.4428
 
 
 def test_make_liquid():
+    """
+
+    Tests creating a liquid `Substance`.
+
+    """
     # Argument types checked
     with pytest.raises(TypeError, match="Name must be a str"):
         Substance.liquid(1, 1, 1)
@@ -69,12 +59,22 @@ def test_make_liquid():
 
 
 def test_liquid(water):
+    """
+
+    Tests that members of a liquid `Substance` are correct.
+
+    """
     assert water.name == 'H2O'
     assert water.mol_weight == 18.0153
     assert water.density == 1
 
 
 def test_make_enzyme():
+    """
+
+    Tests creating an enzyme `Substance`.
+
+    """
     # Argument types checked
     with pytest.raises(TypeError, match="Name must be a str"):
         Substance.enzyme(1)
@@ -83,22 +83,42 @@ def test_make_enzyme():
 
 
 def test_enzyme(lipase):
+    """
+
+    Tests that members of an enzyme `Substance` are correct.
+
+    """
     assert lipase.name == 'lipase'
 
 
 def test_is_solid(salt, water, lipase):
+    """
+
+    Tests that is_solid() returns the correct values.
+
+    """
     assert salt.is_solid() is True
     assert water.is_solid() is False
     assert lipase.is_solid() is False
 
 
 def test_is_liquid(salt, water, lipase):
+    """
+
+    Tests that is_liquid() returns the correct values.
+
+    """
     assert salt.is_liquid() is False
     assert water.is_liquid() is True
     assert lipase.is_liquid() is False
 
 
 def test_is_enzyme(salt, water, lipase):
+    """
+
+    Tests that is_enzyme() returns the correct values.
+
+    """
     assert salt.is_enzyme() is False
     assert water.is_enzyme() is False
     assert lipase.is_enzyme() is True
diff --git a/tests/test_convert.py b/tests/test_convert.py
new file mode 100644
index 0000000..f898c79
--- /dev/null
+++ b/tests/test_convert.py
@@ -0,0 +1,74 @@
+import pytest
+from pyplate.pyplate import Unit
+
+
+def test_convert(salt, water, lipase, dmso):
+    """
+
+    Test internal conversion utilities.
+
+    """
+    # Argument types checked
+    with pytest.raises(TypeError, match='Invalid type for substance'):
+        Unit.convert(None, '', '')
+    with pytest.raises(TypeError, match='Quantity must be a str'):
+        Unit.convert(salt, None, '')
+    with pytest.raises(TypeError, match='Unit must be a str'):
+        Unit.convert(salt, '10 mL', None)
+
+    with pytest.raises(ValueError, match='Only enzymes can be measured in activity units'):
+        Unit.convert(water, '1 U', 'mL')
+    with pytest.raises(ValueError, match='Enzymes can only be measured in activity units'):
+        Unit.convert(lipase, '1 g', 'U')
+
+    # Only enzymes have activity units
+    assert Unit.convert(water, '1 mL', 'U') == 0
+    assert Unit.convert(lipase, '1 U', 'U') == 1
+
+    # prefixes applied correctly
+    assert Unit.convert(water, '1 mL', 'L') == 0.001
+    assert Unit.convert(water, '1 L', 'mL') == 1000
+
+    # from grams
+    assert Unit.convert(water, '1 g', 'g') == 1
+    assert Unit.convert(water, '1 g', 'mL') == 1 / water.density
+    assert Unit.convert(water, '1 g', 'mol') == 1 / water.mol_weight
+
+    # from moles
+    assert Unit.convert(water, '1 mol', 'g') == water.mol_weight
+    assert Unit.convert(water, '1 mol', 'mol') == 1
+    # mol * g/mol / (g/mL)
+    assert Unit.convert(water, '1 mol', 'mL') == water.mol_weight / water.density
+
+    # from L
+    assert Unit.convert(water, '1 L', 'g') == 1000 * water.density
+    assert Unit.convert(water, '1 L', 'mol') == 1000 * water.density / water.mol_weight
+    assert Unit.convert(water, '1 L', 'L') == 1
+
+    # Repeat for a different liquid
+    # from grams
+    assert Unit.convert(dmso, '1 g', 'g') == 1
+    assert Unit.convert(dmso, '1 g', 'mL') == 1 / dmso.density
+    assert Unit.convert(dmso, '1 g', 'mol') == 1 / dmso.mol_weight
+
+    # from moles
+    assert Unit.convert(dmso, '1 mol', 'g') == dmso.mol_weight
+    assert Unit.convert(dmso, '1 mol', 'mol') == 1
+    # mol * g/mol / (g/mL)
+    assert Unit.convert(dmso, '1 mol', 'mL') == dmso.mol_weight / dmso.density
+
+    # from L
+    assert Unit.convert(dmso, '1 L', 'g') == 1000 * dmso.density
+    assert Unit.convert(dmso, '1 L', 'mol') == 1000 * dmso.density / dmso.mol_weight
+    assert Unit.convert(dmso, '1 L', 'L') == 1
+
+    # Repeat for a solid
+    # from grams
+    assert Unit.convert(salt, '1 g', 'g') == 1
+    assert Unit.convert(salt, '1 g', 'mL') == 1 / salt.density
+    assert Unit.convert(salt, '1 g', 'mol') == 1 / salt.mol_weight
+
+    # from moles
+    assert Unit.convert(salt, '1 mol', 'g') == salt.mol_weight
+    assert Unit.convert(salt, '1 mol', 'mol') == 1
+    assert Unit.convert(salt, '1 mol', 'mL') == pytest.approx(salt.mol_weight / salt.density)
diff --git a/tests/test_convert_to_unit_value.py b/tests/test_convert_to_unit_value.py
deleted file mode 100644
index 8f34a50..0000000
--- a/tests/test_convert_to_unit_value.py
+++ /dev/null
@@ -1,102 +0,0 @@
-from pyplate.pyplate import Substance, Unit
-import pytest
-
-EPSILON = 1e-6
-
-
-@pytest.fixture
-def water() -> Substance:
-    return Substance.liquid('H2O', mol_weight=18.0153, density=1)
-
-
-@pytest.fixture
-def salt() -> Substance:
-    return Substance.solid('NaCl', 58.4428)
-
-
-@pytest.fixture
-def dmso() -> Substance:
-    return Substance.liquid('DMSO', 78.13, 1.1004)
-
-
-@pytest.fixture
-def sodium_sulfate() -> Substance:
-    return Substance.solid('Sodium sulfate', 142.04)
-
-
-def test_convert_to_unit_value_solid_g(salt):
-    moles = Unit.convert_to_unit_value(salt, '10 mg')
-    assert abs(moles - Unit.convert_to_storage(0.01 / 58.4428, 'mol')) < EPSILON
-
-    moles = Unit.convert_to_unit_value(salt, '10 g')
-    assert abs(moles - Unit.convert_to_storage(10 / 58.4428, 'mol')) < EPSILON
-
-
-def test_convert_to_unit_value_solid_L(salt):
-    with pytest.raises(ValueError):
-        Unit.convert_to_unit_value(salt, '1 L')
-
-    with pytest.raises(ValueError):
-        Unit.convert_to_unit_value(salt, '1 mL')
-
-
-def test_convert_to_unit_value_solid_mol(salt):
-    moles = Unit.convert_to_unit_value(salt, '1 mol')
-    assert moles == Unit.convert_to_storage(1, 'mol')
-
-
-def test_convert_to_unit_value_solid_U(salt):
-    with pytest.raises(ValueError):
-        Unit.convert_to_unit_value(salt, '1 U')
-
-
-def test_convert_to_unit_value_solid_garbage(salt):
-    with pytest.raises(ValueError):
-        Unit.convert_to_unit_value(salt, '1 garbage')
-
-
-###########################
-def test_convert_to_unit_value_liquid_g(water, dmso):
-    volume = Unit.convert_to_unit_value(water, '1 g')
-    assert volume == Unit.convert_to_storage(1, 'mL')
-
-    volume = Unit.convert_to_unit_value(water, '1 mg')
-    assert volume == Unit.convert_to_storage(1e-3, 'mL')
-
-    volume = Unit.convert_to_unit_value(dmso, '1 g')
-    assert volume == Unit.convert_to_storage(1/1.1004, 'mL')
-
-
-def test_convert_to_unit_value_liquid_L(water):
-    volume = Unit.convert_to_unit_value(water, '1 L')
-    assert volume == Unit.convert_to_storage(1000, 'mL')
-
-    volume = Unit.convert_to_unit_value(water, '1 mL')
-    assert volume == Unit.convert_to_storage(1, 'mL')
-
-
-def test_convert_to_unit_value_liquid_mol(water, dmso):
-    volume = Unit.convert_to_unit_value(water, '1 mol')
-    assert abs(volume - Unit.convert_to_storage(18.0153, 'mL')) < EPSILON
-    volume = Unit.convert_to_unit_value(dmso, '1 mol')
-    assert abs(volume - Unit.convert_to_storage(78.13/1.1004, 'mL')) < EPSILON
-
-
-def test_convert_to_unit_value_liquid_U(water):
-    with pytest.raises(ValueError):
-        Unit.convert_to_unit_value(water, '1 U')
-
-
-def test_convert_to_unit_value_liquid_garbage(water):
-    with pytest.raises(ValueError):
-        Unit.convert_to_unit_value(water, '1 garbage')
-
-
-def test_convert_to_unit_value_enzyme():
-    lactase = Substance.enzyme('Lactase')
-    for amount in ['1 mg', '1 mL', '1 mol', '1 garbage']:
-        with pytest.raises(ValueError):
-            Unit.convert_to_unit_value(lactase, amount)
-
-    assert Unit.convert_to_unit_value(lactase, '1 U') == 1
-    assert Unit.convert_to_unit_value(lactase, '10 U') == 10
diff --git a/tests/test_create_solution.py b/tests/test_create_solution.py
new file mode 100644
index 0000000..ed038ef
--- /dev/null
+++ b/tests/test_create_solution.py
@@ -0,0 +1,61 @@
+import math
+from itertools import product
+from pyplate.pyplate import Unit, Container, config
+
+epsilon = 1e-6
+
+
+def test_create_solution(salt, water, triethylamine, dmso, sodium_sulfate, lipase):
+    """
+
+    Create a solution using each a quantity of each solvent and solute in each unit.
+    Try "0.001 numerator/denominator" and "0.01 numerator/10 denominator"
+    Ensure the correct amount of solvent, solute, and total solution is applied.
+    Also try to create a "1.1 U/denominator" enzyme solution for each solvent.
+
+    """
+    solvents = [water, dmso]
+    solutes = [salt, triethylamine, sodium_sulfate]
+    units = ['g', 'mol', 'mL']
+    for numerator, denominator, quantity_unit in product(units, repeat=3):
+        for solute in solutes:
+            for solvent in solvents:
+                if numerator == 'mL' and solute.is_solid() and config.solid_density == float('inf'):
+                    continue
+                con = Container.create_solution(solute, solvent, concentration=f"0.001 {numerator}/{denominator}",
+                                                total_quantity=f"10 {quantity_unit}")
+                assert all(value > 0 for value in con.contents.values())
+                total = sum(Unit.convert(substance, f"{value} {config.moles_prefix}", quantity_unit) for
+                            substance, value in con.contents.items())
+                assert abs(total - 10) < epsilon, f"Making 10 {quantity_unit} of a 0.001 {numerator}/{denominator}" \
+                                                  f" solution of {solute} and {solvent} failed."
+                conc = Unit.convert(solute, f"{con.contents[solute]} {config.moles_prefix}", numerator) / \
+                    sum(Unit.convert(substance, f"{value} {config.moles_prefix}", denominator)
+                        for substance, value in con.contents.items())
+                assert abs(conc - 0.001) < epsilon
+
+                con = Container.create_solution(solute, solvent, concentration=f"0.01 {numerator}/10 {denominator}",
+                                                total_quantity=f"10 {quantity_unit}")
+                total = sum(Unit.convert(substance, f"{value} {config.moles_prefix}", quantity_unit) for
+                            substance, value in con.contents.items())
+                assert abs(total - 10) < epsilon
+                conc = Unit.convert(solute, f"{con.contents[solute]} {config.moles_prefix}", numerator) / \
+                    sum(Unit.convert(substance, f"{value} {config.moles_prefix}", denominator)
+                        for substance, value in con.contents.items())
+                assert abs(conc - 0.01/10) < epsilon
+
+    # Solute is an enzyme, concentration has U in the numerator
+    solute = lipase
+    for denominator, quantity_unit in product(units, repeat=2):
+        for solvent in solvents:
+            con = Container.create_solution(solute, solvent, concentration=f"1.1 U/{denominator}",
+                                            total_quantity=f"10 {quantity_unit}")
+            assert all(value > 0 for value in con.contents.values())
+            total = sum(Unit.convert(substance, f"{value} {config.moles_prefix}", quantity_unit) for
+                        substance, value in con.contents.items())
+            assert abs(total - 10) < epsilon, f"Making 10 {quantity_unit} of a 1.1 U/{denominator}" \
+                                              f" solution of {solute} and {solvent} failed."
+            conc = con.contents[solute] / \
+                sum(Unit.convert(substance, f"{value} {config.moles_prefix}", denominator)
+                    for substance, value in con.contents.items())
+            assert abs(conc - 1.1) < epsilon
diff --git a/tests/test_dilute.py b/tests/test_dilute.py
new file mode 100644
index 0000000..dd1d194
--- /dev/null
+++ b/tests/test_dilute.py
@@ -0,0 +1,26 @@
+from itertools import product
+import pytest
+from pyplate.pyplate import Container, Unit, config
+
+
+def test_dilute(salt, sodium_sulfate, triethylamine, water, dmso):
+    """
+    Test diluting solutions.
+    Create a solution of 10 quantity_unit of 0.001 numerator/denominator for each solute
+    and solvent and then dilute to 0.0005 numerator/denominator.
+    Ensure the new concentration is correct.
+
+    """
+    solvents = [water, dmso]
+    solutes = [salt, triethylamine, sodium_sulfate]
+    units = ['g', 'mol', 'mL']
+    for numerator, denominator, quantity_unit in product(units, repeat=3):
+        for solute in solutes:
+            for solvent in solvents:
+                con = Container.create_solution(solute, solvent, concentration=f"0.001 {numerator}/{denominator}",
+                                                total_quantity=f"10 {quantity_unit}")
+                con2 = con.dilute(solute, f'0.0005 {numerator}/{denominator}', solvent)
+                new_concentration = Unit.convert(solute, f"{con2.contents[solute]} {config.moles_prefix}", numerator) / \
+                                    sum(Unit.convert(substance, f"{value} {config.moles_prefix}", denominator)
+                                        for substance, value in con2.contents.items())
+                assert new_concentration == pytest.approx(0.0005)
diff --git a/tests/test_transfers.py b/tests/test_transfers.py
index f574c95..cea61d4 100644
--- a/tests/test_transfers.py
+++ b/tests/test_transfers.py
@@ -1,118 +1,69 @@
-from pyplate.pyplate import Plate, Substance, Container, Unit
 import pytest
 import numpy
-
-
-@pytest.fixture
-def water() -> Substance:
-    return Substance.liquid('H2O', mol_weight=18.0153, density=1)
-
-
-@pytest.fixture
-def salt() -> Substance:
-    return Substance.solid('NaCl', 58.4428)
-
-
-@pytest.fixture
-def dmso() -> Substance:
-    return Substance.liquid('DMSO', 78.13, 1.1004)
-
-
-@pytest.fixture
-def sodium_sulfate() -> Substance:
-    return Substance.solid('Sodium sulfate', 142.04)
+from pyplate import Plate, Container
+from pyplate.pyplate import Unit, config
 
 
 @pytest.fixture
 def plate1() -> Plate:
-    return Plate('plate1', 50_000)
+    return Plate('plate1', '50 mL')
 
 
 @pytest.fixture
 def plate2() -> Plate:
-    return Plate('plate2', 50_000)
+    return Plate('plate2', '50 mL')
 
 
 @pytest.fixture
 def solution1(water, salt) -> Container:
-    return Container('sol1', max_volume=1_000, initial_contents=[(water, '100 mL'), (salt, '50 mol')])
+    return Container('sol1', initial_contents=[(water, '100 mL'), (salt, '50 mmol')])
 
 
 @pytest.fixture
 def solution2(dmso, sodium_sulfate) -> Container:
-    return Container('sol2', max_volume=1_000, initial_contents=[(dmso, '100 mL'), (sodium_sulfate, '50 mol')])
+    return Container('sol2', initial_contents=[(dmso, '100 mL'), (sodium_sulfate, '50 mmol')])
 
 
-def test_add_to_container(salt, water):
-    """
-    Tests adding a substance to a container.
-    """
-    container = Container('container', max_volume=10)
-    container = Container.add(salt, container, '10 mol')
-
-    # container should contain 10 moles of salt
-    assert salt in container.contents and Unit.convert_from_storage(container.contents[salt], 'mol') == 10
-
-    container1 = Container('container', max_volume=20)
-    container1 = Container.add(water, container1, '10 mL')
-    # container1 should contain 10 mL of water and its volume should be 10 mL
-    assert water in container1.contents and Unit.convert_from_storage(container1.contents[water], 'mL') == 10
-    assert Unit.convert_from_storage(container1.volume, 'mL') == 10
-
-    container2 = Container.add(salt, container1, '5 mol')
-    # container2's volume shouldn't be changed by adding salt
-    assert Unit.convert_from_storage(container2.volume, 'mL') == 10
-
-    container3 = Container.add(water, container2, '5 mL')
-    # water should stack in the contents and the volume should be updated
-    assert Unit.convert_from_storage(container3.contents[water], 'mL') == 15 and \
-           Unit.convert_from_storage(container3.volume, 'mL') == 15
-
-    # the original container should be unchanged
-    assert container.volume == 0 and \
-           Unit.convert_from_storage(container.contents[salt], 'mol') == 10
-
-
-def test_transfer_between_containers(solution1, solution2, water, salt):
+def test_transfer_between_containers(solution1, solution2, water, salt, dmso, sodium_sulfate):
     """
     Tests transferring from one container to another.
     """
-    solution3, solution4 = Container.transfer(solution1, solution2, '10 mL')  # solution2.transfer(solution1, '10 mL')
+    # solution1 has 100mL of water and 50 nmol of solt
+    # solution2 has 100mL of dmso and 50 nmol of sodium sulfate
+    solution1_volume = 100 + Unit.convert(salt, '50 mmol', 'mL')
+    solution2_volume = 100 + Unit.convert(sodium_sulfate, '50 mmol', 'mL')
+    solution3, solution4 = Container.transfer(solution1, solution2, f"{solution1_volume*0.1} mL")
     # original solutions should be unchanged
-    assert Unit.convert_from_storage(solution2.volume, 'mL') == 100 and \
-           Unit.convert_from_storage(solution1.volume, 'mL') == 100
+    assert solution1.volume == Unit.convert_to_storage(solution1_volume, 'mL')
+    assert solution2.volume == Unit.convert_to_storage(solution2_volume, 'mL')
     # 10 mL of water and 5 moles of salt should have been transferred
-    assert Unit.convert_from_storage(solution3.volume, 'mL') == 90 and \
-           Unit.convert_from_storage(solution4.volume, 'mL') == 110
-    assert Unit.convert_from_storage(solution3.contents[water], 'mL') == 90 and \
-           Unit.convert_from_storage(solution3.contents[salt], 'mol') == 45
-    assert Unit.convert_from_storage(solution4.contents[water], 'mL') == 10 and \
-           Unit.convert_from_storage(solution4.contents[salt], 'mol') == 5
+    assert solution3.volume == Unit.convert_to_storage(solution1_volume*0.9, 'mL')
+    assert solution4.volume == Unit.convert_to_storage(solution2_volume + solution1_volume*0.1, 'mL')
+    assert solution3.contents[water] == Unit.convert(water, '90 mL', config.moles_prefix)
+    assert solution3.contents[salt] == Unit.convert(salt, '45 mmol', config.moles_prefix)
+    assert solution4.contents[water] == pytest.approx(Unit.convert(water, '10 mL', config.moles_prefix))
+    assert solution4.contents[salt] == Unit.convert(salt, '5 mmol', config.moles_prefix)
 
 
-def test_add_to_slice(plate1, salt):
-    """
-    Tests adding a substance to each well in a slice.
-    """
-    plate3 = Plate.add(salt, plate1[:], '10 mol')
-    # 10 moles of salt should be in each well
-    assert numpy.array_equal(plate3.moles(salt), numpy.full(plate3.wells.shape, 10))
-    # Original plate should be unchanged
-    assert numpy.array_equal(plate1.moles(salt), numpy.zeros(plate1.wells.shape))
-
-
-def test_transfer_to_slice(plate1, solution1):
+def test_transfer_to_slice(plate1, solution1, salt):
     """
     Tests transferring from a container to each well in a slice.
     """
-    solution3, plate3 = Plate.transfer(solution1, plate1[:], '1 mL')  # plate1[:].transfer(solution1, '1 mL')
-    # 1 mL should have been transferred to each well in the plate
-    assert plate3.volume() == plate3[:].size * 1000         # volume() is in uL
-    assert numpy.all(numpy.vectorize(lambda elem: Unit.convert_from_storage(elem.volume, 'mL') == 1)(plate3.wells))
+    solution1_volume = 100 + Unit.convert(salt, '50 mmol', 'mL')
+    to_transfer = solution1_volume / 100
+    solution3, plate3 = Plate.transfer(solution1, plate1[:], f"{to_transfer} mL")
+    # 1 mL of water should have been transferred to each well in the plate
+    assert plate3.volume() == pytest.approx(plate3[:].size * 1000 * to_transfer)  # volume() is in uL
+    assert numpy.allclose(plate3.volumes(unit='mL'), numpy.ones(plate3.wells.shape) * to_transfer)
+    # assert numpy.all(numpy.vectorize(lambda elem: abs(elem.volume - to_transfer) < epsilon)(plate3.wells))
     # Original solution and plate should be unchanged
-    assert Unit.convert_from_storage(solution1.volume, 'mL') == 100
+    assert solution1.volume == Unit.convert_to_storage(solution1_volume, 'mL')
     assert plate1.volume() == 0
 
+    solution4, plate4 = Plate.transfer(solution1, plate1[1, 1], f"{to_transfer} mL")
+    assert plate4.volume() == pytest.approx(1000 * to_transfer)
+    assert solution1.volume - solution4.volume == pytest.approx(Unit.convert_to_storage(to_transfer, 'mL'))
+
 
 def test_transfer_between_slices(plate1, plate2, solution1, solution2):
     """
@@ -154,8 +105,7 @@ def test_transfer_from_slice(plate1, solution1):
     """
     # solution3, plate3 = plate1[:].transfer(solution1, '1 mL')
     solution3, plate3 = Plate.transfer(solution1, plate1[:], '1 mL')
-    destination_solution = Container('destination', 100)
-    # plate4, destination_solution = destination_solution.transfer_slice(plate3, '0.5 mL')
+    destination_solution = Container('destination', '100 mL')
     plate4, destination_solution = Container.transfer(plate3, destination_solution, '0.5 mL')
     # Original plate should have 1 mL in each well
     assert numpy.array_equal(plate3.volumes(), numpy.ones(plate3[:].shape) * 1000)
@@ -163,3 +113,6 @@ def test_transfer_from_slice(plate1, solution1):
     assert destination_solution.volume == Unit.convert_to_storage(0.5 * plate3[:].size, 'mL')
     # Source wells should all have 0.5 mL
     assert numpy.array_equal(plate4.volumes(), numpy.ones(plate4[:].shape) * 0.5 * 1000)
+
+    plate4, destination_solution = Container.transfer(plate3[1, :], destination_solution, '0.5 mL')
+    assert plate3.volume() - plate4.volume() == 500 * plate3.n_columns