Merge pull request #13 from jschnab/dev

Add more data fitting models and a tutorial

README.md

@@ -56,7 +56,7 @@ Available subcommands:
 * `default`: rollback to the original configuration (the one provided when you install this package)
 
 When you enter the configuration mode, instructions will be displayed before you are prompted for input. The name of the parameter is displayed followed by its current value between ( parentheses ).
-You have three possible actions (validate any of them by pressing <Enter>):
+You have three possible actions (validate any of them by pressing \<Enter\>):
 * keep the current value: leave the input field empty
 * modify the current value: type your input (if a list is expected, use space to separate values)
 * enter a void value (no value): enter 'none' (without quotes)
@@ -78,17 +78,17 @@ This subcommand plots absorption spectra.
 eda plot spectrum [arguments ...]
 ```
 File names are positional arguments, they should be passed before optional argument. Optional arguments include:
-* `-l` or `--label` specify labels on data for the plot legend
-* `--figure-size` specify width and height in inches
+* `-l` or `--label` labels for the plot legend
+* `--figure-size` width and height in inches
 * `--xcolumn` name of the column containing x-axis values
 * `--ycolumn` name of the column containing y-axis values
-* `--xlabel` label on the plot's x-axis
-* `--ylabel` label on the plot's y-axis
-* `--xlimit` specify left and right values for x-axis limits
-* `--ylimit` specify bottom and top values for y-axis limits
+* `--xlabel` label on the x-axis
+* `--ylabel` label on the y-axis
+* `--xlimit` left and right values for x-axis limits
+* `--ylimit` bottom and top values for y-axis limits
 * `--skip-header` number of rows to skip at the beginning of the file
 * `--legend-location` run `eda plot spectrum -h` for more information
-* `--title` specify the title of the plot
+* `--title` title of the plot
 
 For example:
 ```
@@ -108,25 +108,83 @@ eda plot kinetics [arguments ...]
 ```
 
 File names are positional arguments and should be passed before optional arguments. Optional arguments include:
-* `-l` or `--label` specify labels on data for the plot legend
-* `-m` or `--model` is a flag indicating an exponential model should be fitted on the data
-* `--figure-size` specify width and height in inches
+* `-l` or `--label` labels for the plot legend
+* `-f` or `--fit` fit the data with a mathematical model
+* `-m` or `--model` specify the mathematical model used to fit the data. Choices include:
+	- `exp` (default) fit both first-order and second-order exponential models and selects the best
+	- `exp1` fit a first-order exponential
+	- `exp2` fit a second-order exponential
+	- `linear` fit a linear model
+* `--init-params` provide initial parameters for the curve-fitting algorithm
+* `--skip-header` number of rows to skip at the beginning of the file
 * `--xcolumn` name of the column containing x-axis values
 * `--ycolumn` name of the column containing y-axis values
-* `--xlabel` label on the plot's x-axis
-* `--ylabel` label on the plot's y-axis
-* `--xlimit` specify left and right values for x-axis limits
-* `--ylimit` specify bottom and top values for y-axis limits
-* `--skip-header` number of rows to skip at the beginning of the file
+* `--xlabel` label on the x-axis
+* `--ylabel` label on the y-axis
+* `--xlimit` left and right values for x-axis limits
+* `--ylimit` bottom and top values for y-axis limits
+* `--figure-size` width and height in inches
 * `--legend-location` run `eda plot spectrum -h` for more information
-* `--title` specify the title of the plot
+* `--title` title of the plot
 
 For example:
 ```
-eda plot kinetics file1.csv file2.csv -l experiment1 experiment2 -m
+eda plot kinetics file1.csv file2.csv -l experiment1 experiment2 -f
 ```
 
 For more information:
 ```
 eda plot kinetics -h
 ```
+
+## Tutorials
+
+You should do these two steps prior to the tutorial:
+* install the `easy-data-analysis` package
+* download the CSV data files from `eda/docs/`
+
+### Plot a kinetics curve
+
+
+Plotting a kinetics curve is done by providing the file names(s) to `eda plot kinetics`:
+```
+eda plot kinetics first_expo.csv
+```
+
+You should see this output:
+
+<img src="eda/docs/images/first_order.png" width="600" alt="first order exponential plot">
+
+
+You do not need to provide optional arguments because the format of `first_order.csv` matches the default configuration. Files with various formats can be dealt with by changing or reviewing the configuration by running `eda configure kinetics` or by providing arguments to `eda plot kinetics`.
+
+Let's plot and fit `second_expo.csv`. This file contains no rows to be skipped (the first line contains column names).
+```
+eda plot kinetics second_expo.csv -f -m exp2 --xcolumn time --ycolumn absorbance --xlabel time --skip-header 0
+```
+
+You should see this output:
+
+<img src="eda/docs/images/second_order.png" width="600" alt="second order exponential plot">
+
+When fitting a curve you will be provided with fitting results:
+```
+second_expo.csv
+y = a1 * exp(k1 * x) + a2 * exp(k2 * x)
+------------------------------
+Parameter      Value   Std Err
+------------------------------
+a1           -0.1169    0.0069
+a2           +0.7081    0.0051
+k1           -1.0201    0.1330
+k2           +0.0063    0.0010
+R-square     0.91967
+t1 (sec)       40.77
+t2 (sec)     6600.02
+```
+
+The first line is the file name, the second line is the equation of the data model. The following lines show the value and standard error of the parameters of the equation. The R-square value indicates the goodness of fit and varies from 0 (poor fit) to 1 (perfect fit). The parameters `t1` and `t2` are the doubling times (or halving times for exponential decay) of the first and second components of the equation. They are calculated as:
+
+<a href="https://www.codecogs.com/eqnedit.php?latex=t_x=\frac{\ln2}{|k_x|}" target="_blank"><img src="https://latex.codecogs.com/gif.latex?t_x=\frac{\ln2}{|k_x|}" title="t_x=\frac{\ln2}{|k_x|}" /></a>
+
+`t1` and `t2` are shown in seconds and assume that you provide data in *minute* by default. If the time unit of your data is *second*, you can specify it by running `eda configure kinetics` and modify the parameter `time_unit` or by using the parameter `--time-unit` when calling `eda plot kinetics`.

eda/cli.py

@@ -81,10 +81,7 @@ def plot(self):
 
     def spectrum(self):
         parser = ArgumentParser(
-            description=(
-                "Plot an absorbance spectrum. "
-                "For more information: eda plot spectrum -h"
-            ),
+            description=("Plot an absorbance spectrum"),
             usage=(
                 "eda plot spectrum file [file ...] [-h] [-l] [--figure-size] "
                 "[--xcolumn] [--ycolumn] [--xlabel] [--ylabel] [--xlimit] "
@@ -190,26 +187,48 @@ def spectrum(self):
     def kinetics(self):
         parser = ArgumentParser(
             description=(
-                "Plot an absorbance kinetics curve. "
-                "For more information: eda plot kinetics -h"
-            ),
+                "Plot an absorbance kinetics curve"),
             usage=(
-                "eda plot spectrum file [file ...] [-h] [-m] [-l] "
+                "eda plot spectrum file [file ...] [-h] [-f] [-m] [-l] "
                 "[--figure-size] [--xcolumn] [--ycolumn] [--xlabel] [--ylabel]"
                 " [--xlimit] [--ylimit] [--skip-header] [--legend-location] "
-                "[--title]"
+                "[--title] [--time-unit] [--init-params]"
             ),
         )
         parser.add_argument(
             "file",
             nargs="+",
             help="CSV files storing data to plot",
         )
+        parser.add_argument(
+            "-f",
+            "--fit",
+            action="store_true",
+            help=(
+                "Fit the data using a model, specify model with the '-f' or "
+                "'--fit' argument)"
+            ),
+        )
         parser.add_argument(
             "-m",
             "--model",
-            action="store_true",
-            help="Model the data using exponential decay",
+            choices=["linear", "exp1", "exp2", "exp"],
+            help=(
+                "Specify the model to use when fitting data, choose 'linear' "
+                "(linear model), 'exp1' (first-order exponential), 'exp2' "
+                "(second-order exponential), or 'exp' (select the best fit "
+                "between 'exp1' and 'exp2')"
+            ),
+        )
+        parser.add_argument(
+            "--init-params",
+            nargs="+",
+            type=float,
+            help=(
+                "The initial parameters to use when fitting data (make sure "
+                "the number of parameters is appropriate for the selected "
+                "model)"
+            ),
         )
         parser.add_argument(
             "-l",
@@ -284,12 +303,19 @@ def kinetics(self):
             "--title",
             nargs="?",
             type=str,
-            help="Specify the title of the plot",
+            help="Title of the plot",
+        )
+        parser.add_argument(
+            "--time-unit",
+            choices=["minute", "second"],
+            help="Time unit of the kinetics experiment",
         )
         args = parser.parse_args(sys.argv[3:])
         kinetics.run(
             input_files=args.file,
+            fit=args.fit,
             model=args.model,
+            init_params=args.init_params,
             labels=args.label,
             fig_size=args.fig_size,
             x_col=args.xcolumn,
@@ -301,6 +327,7 @@ def kinetics(self):
             skip_header=args.skip_header,
             legend_loc=args.legend_loc,
             title=args.title,
+            time_unit=args.time_unit,
         )
 
 

eda/config_default.yaml

@@ -12,7 +12,12 @@ plot:
         legend_location: "upper left"
         title:
     kinetics:
-        model: True 
+        fit: True
+        model: "exp"
+        init_params:
+            exp1: [1, -1, 1]
+            exp2: [1, 1, 0, 0]
+            linear: [1, 1]
         use_cols: [0, 1]
         figure_size: [7, 5]
         xcolumn: "Time (min)"
@@ -24,4 +29,4 @@ plot:
         skip_header: 1
         legend_location: "lower right"
         title:
-
+        time_unit: "minute"

eda/configure.py

@@ -17,7 +17,10 @@
     "skip_header": lambda x: int(x),
     "legend_location": lambda x: x,
     "title": lambda x: x,
-    "model": lambda x: True if x.lower() == 'true' else False,
+    "fit": lambda x: True if x.lower() == 'true' else False,
+    "model": lambda x: x,
+    "init_params": lambda x: list(map(float, x.split())),
+    "time_unit": lambda x: x,
 }
 
 
@@ -32,7 +35,7 @@ def __init__(self):
 
     def print_info(self):
         print()
-        print(f"Configuration instructions")
+        print("Configuration instructions")
         print("--------------------------")
         print("The current value is displayed between ( parentheses ).")
         print("The absence of value is indicated by 'None'.")
@@ -52,6 +55,24 @@ def record(self, subcommand):
         inputs = {}
         for key, value in config["plot"][subcommand].items():
             while True:
+                # special case: init_params is a nested dictionary
+                if key == "init_params":
+                    print("init_params:")
+                    init_params = {}
+                    for model, params in value.items():
+                        while True:
+                            user_input = input(f"    {model} ( {params} ):")
+                            if not user_input or user_input.lower() == "none":
+                                init_params[model] = params
+                                break
+                            try:
+                                t = transform_input["init_params"](user_input)
+                            except ValueError:
+                                print("Invalid value")
+                                continue
+                            init_params[model] = t
+                    inputs[key] = init_params
+                    break
                 user_input = input(f"{key} ( {value} ): ")
                 if not user_input:
                     inputs[key] = value