Add docs for microwave cycle fitting (#18)

* Fix on docs what is already in README

* Minor revision

* Add a script, still needs docs comments

* Put some better text in the example for documentation

* Touch up docs elsewhere

Author: Ickaser

docs/example/2023-03-02_RF_Mannitol_process.csv

@@ -1,10 +1,9 @@
-Run started: 6/4/2024 10:38:37 AM,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-Product Name: Man5 controlled anneal,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-Product Number: 1,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-Operator: ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-Recipe File Name: man5.rcp,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+Run started: 6/4/2024 10:38:37 AM
+
+Product Name: Man5 controlled anneal
+Product Number: 1
+Recipe File Name: man5.rcp
+
 CycleTime,Cycle,Phase,Step,VacSetPT,VacPirani,VacCPM,ShelfSetPT,ShelfInlet,Condenser,TPAvg,TP1,TP2,TP3,TP4,TP5,TP6,TP7,TP8,TShelfSurface,TProdConsider,HeatFlux,KvShelf,MassFlowVial,MassFlowHR,RProd,PercentDried,GroupAvg,Qshelf,KvSys_SS,KvFinal_SS,HFSystem,Percent_FRZ,LyoSIM,LyoSIM_SetPT,Percent_RM
 10:38:38,9,1,0,500,738690,0,-25,2.6,9.6,3.4,5.2,3.2,999.9,1.9,999.9,999.9,999.9,999.9,0,5.2,-4004027,0,0,0,0,0,5.2,100,0,0,-4004027,NaN,1.3,-25,0
 10:39:39,9,1,1,0,463110,0,13.1,1.3,9.9,3.2,5.1,2.9,999.9,1.6,999.9,999.9,999.9,999.9,-0.6,5.1,0,0,0,0,0,0,5.1,100,0,0,0,NaN,-2.4,3,0

docs/example/2023-03-02_RF_Mannitol_temperature.csv

@@ -25,9 +25,11 @@ using Accessors
 
 # # Read Process Data
 
+# This needs to point to the right file, which for documentation is kinda wonky;
+# adjust this in to refer to an appropriate file. 
 ## Data start at 8th row of CSV file.
-## This needs to point to the right file, which for documentation is kinda wonky
-procdata_raw = CSV.read(joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv"), Table, header=8)
+file_loc = joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv")
+procdata_raw = CSV.read(file_loc, Table, header=7)
 t = uconvert.(u"hr", procdata_raw.CycleTime .- procdata_raw.CycleTime[1])
 ## At midnight, timestamps revert to zero, so catch that case
 for i in eachindex(t)[begin+1:end]

docs/example/2024-06-04-10_MFD_AH.csv

@@ -3,9 +3,10 @@
 # `using LyoPronto`, you have effectively also done `using Unitful` and a few others.
 using LyoPronto
 
-# These are other packages that I use in the test suite,
+# These packages are used in the test suite,
 # but you can use others in their place.
-# TypedTables provides a lightweight table structure, not as broadly flexible as a DataFrame but great for our needs
+
+# TypedTables provides a lightweight table structure, not as broadly flexible as a DataFrame but great for our needs.
 using TypedTables, CSV
 # TransformVariables provides tools for mapping optimization parameters to sensible ranges.
 using TransformVariables
@@ -25,7 +26,8 @@ using LaTeXStrings
 
 ## Data start at 8th row of CSV file.
 ## This needs to point to the right file, which for documentation is kinda wonky
-procdata_raw = CSV.read(joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv"), Table, header=8)
+file_loc = joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv")
+procdata_raw = CSV.read(file_loc, Table, header=7)
 ## MicroFD, used for this experiment, has a column indicating primary drying
 t = uconvert.(u"hr", procdata_raw.CycleTime .- procdata_raw.CycleTime[1])
 ## At midnight, timestamps revert to zero, so catch that case
@@ -35,6 +37,12 @@ for i in eachindex(t)[begin+1:end]
     end
 end
 
+cp(file_loc, "./2024-06-04-10_MFD_AH.csv"); #md #hide
+cp(file_loc[1:end-4] * ".yaml", "./2024-06-04-10_MFD_AH.yaml"); #md #hide
+# If you want to follow along exactly, you can download the CSV file 
+# [here](2024-06-04-10_MFD_AH.csv), and some metadata about the cycle 
+# [here](2024-06-04-10_MFD_AH.yaml).
+
 ## Rename the columns we will use, and add units
 procdata = map(procdata_raw) do row
     ## In the anonymous `do` function, `row` is a row of the table.

docs/example/all_recipes.jl

@@ -0,0 +1,216 @@
+# # Imports
+using LyoPronto 
+
+# NonlinearSolve and TransformVariables are used for parameter fitting
+using NonlinearSolve
+using TransformVariables
+
+# CSV and TypedTables load and store experimental data
+using TypedTables, CSV
+
+# Plots is a frontend to several plotting packages, defaulting to GR
+using Plots
+
+## All we need from StatsPlots is the @df macro
+using StatsPlots: @df
+
+## Latexify helps with LaTeX formatting in plots
+using Latexify
+using LaTeXStrings
+## Set a default for Latexify labels
+set_default(labelformat=:square) # Latexify, not Plots
+
+
+
+# # Load process data 
+# The process data here are from an actual microwave-assisted lyophilization cycle.
+# Since thermocouples cannot be used in a microwave field, the product temperature
+# measurements (from fiber optic probes) are in a separate file, and we must take care to 
+# synchronize the time.
+
+# As in the other examples, the file locations here are wonky in order to execute this 
+# documentation remotely; to follow along, adjust the file paths to match your local setup.
+
+filesdir = joinpath(@__DIR__, "..", "..", "example")
+
+dat1 = CSV.read(joinpath(filesdir, "2023-03-02_RF_Mannitol_temperature.csv"), Table)
+dat2 = CSV.read(joinpath(filesdir, "2023-03-02_RF_Mannitol_process.csv"), Table,
+    comment="#", stripwhitespace=true, dateformat="mm/dd/yyyy H:M:S")
+
+cp(joinpath(filesdir, "2023-03-02_RF_Mannitol_temperature.csv"), "./2023-03-02_RF_Mannitol_temperature.csv"); #md #hide
+cp(joinpath(filesdir, "2023-03-02_RF_Mannitol_process.csv"), "./2023-03-02_RF_Mannitol_process.csv"); #md #hide
+# To follow along, you can use the same data files [here](2023-03-02_RF_Mannitol_process.csv)
+# and [here](2023-03-02_RF_Mannitol_temperature.csv).
+
+time1 = dat1.var"Elapsed [s]"
+
+lyo_full_pre = map(dat2) do row
+    nt = (tstamp = row.Timestamp,
+          pch_sp = row.var"SPLYO.VACUUM_SP.F_CV"*u"mTorr",
+          pch_pir = row.var"SPLYO.CHAMBER_PIRANI.F_CV"*u"mTorr",
+          pch_cm = row.var"SPLYO.CHAMBER_CM.F_CV"*u"mTorr",
+          Tsh_sp = row.var"SPLYO.SHELF_SP.F_CV"*u"°C",
+          Tsh_i = row.var"SPLYO.SHELF_INLET.F_CV"*u"°C",
+          Tsh_o = row.var"SPLYO.SHELF_OUTLET.F_CV"*u"°C",
+    )
+    return nt
+end
+time2 = uconvert.(u"hr", lyo_full_pre.tstamp .- lyo_full_pre.tstamp[1])
+lyo_full = Table(lyo_full_pre; t=time2)
+
+# Here we are concerned with the primary drying step, so we need to identify that.
+istart_lyo_pd = findfirst(lyo_full.pch_sp .== 100u"mTorr")
+lyo_pd = Table(lyo_full[istart_lyo_pd:end])
+## Set start of primary drying as zero time
+lyo_pd.t .-= lyo_pd.t[1]
+nothing #md #hide
+
+# It's a good idea at this point to plot temperatures and pressures, 
+# to make sure everything looks right.
+plT = plot(lyo_pd.t, lyo_pd.Tsh_sp)
+plp = plot(lyo_pd.t, lyo_pd.pch_pir, ylim=(0, 300))
+plot!(plp, lyo_pd.t, lyo_pd.pch_cm)
+plot(plT, plp, layout=(2,1), link=:x)
+
+# We will use this to set the shelf temperature and pressure that are used later in fitting.
+Tsh = RampedVariable(uconvert.(u"K", [-40.0u"°C", 10.0u"°C"]), 0.5u"K/minute")
+pch = RampedVariable(100u"mTorr")
+
+# We will also need to know when primary drying *actually* ends, 
+# so we will use the second-derivative test in Pirani-CM convergence. 
+t_end = identify_pd_end(lyo_pd.t, lyo_pd.pch_pir, Val(:der2))
+@df lyo_pd plot(:t, :pch_pir)
+@df plot!(:t, pch_pir_sm)
+tendplot!(t_end)
+
+
+# Next, we need to read the temperature data from the fiber optic probes.
+thm_full = map(dat1) do row
+    nt = (tstamp = row.var"Elapsed [s]"*u"s",
+          T1 = row.var"T1 [C]"*u"°C",
+          T2 = row.var"T2 [C]"*u"°C",
+          T3 = row.var"T3 [C]"*u"°C",
+          T4 = row.var"T4 [C]"*u"°C",
+    )
+    return nt
+end
+@df thm_full exptfplot(:tstamp, :T1, :T2, :T3, :T4, nmarks=40, xunit=u"hr")
+
+# Clearly we can ignore T2, which is not reading anything physical.
+# Also, note that one of our three remaining thermal probes is on the outside wall of a vial,
+# while the other two are in frozen product. Based on the behavior we see in plots, it is 
+# clear that the series denoted T3 is the odd one out.
+
+# Next, we need to identify the portion of the temperature data that corresponds to primary drying.
+# After some trial and error, I decided the following achieves that:
+istart_thm_pd = argmin(thm_full.T1)
+iend_thm_pd = argmax(thm_full.T4) + 600
+
+thm_pd_sub = Table(thm_full[istart_thm_pd:iend_thm_pd])
+thm_pd = Table(thm_pd_sub; t = uconvert.(u"hr", thm_pd_sub.tstamp .- thm_pd_sub.tstamp[1]))
+
+## Plot our temperatures to make sure everything looks right
+plot(u"hr", u"°C", xlabel="Time", ylabel="Temperature", size=(400,300), legend=false)
+@df thm_pd exptfplot!(:t, :T1, :T4, nmarks=40)
+@df thm_pd exptvwplot!(:t, :T3, nmarks=40)
+plot!(Tsh, c=:black, label="shelf", tmax=13.5u"hr")
+
+# Based on that plot, we identify the part that we want to fit:
+# everything up until the temperature minimum near 10 hours.
+iend_T4 = argmin(thm_pd.T4[1000:end]) + 1000
+iend_T1 = argmin(thm_pd.T1[1000:end]) + 1000
+## Factor of 6 is because the temperature data are every 10 seconds, 
+## compared to every minute for process data
+fitdat = @df thm_pd[begin:6:end] PrimaryDryFit(:t, (:T4[begin:iend_T4÷6], 
+    :T1[begin:iend_T1÷6]), :T3, t_end);
+plot(fitdat)
+
+# # Set up other model parameters
+## Vial parameters
+vialsize = "6R"
+rad_i, rad_o = get_vial_radii(vialsize)
+A_p = π*rad_i^2  # cross-sectional area inside the vial
+A_v = π*rad_o^2 # vial bottom area
+m_v = get_vial_mass(vialsize)
+## Formulation and fill
+c_solid = 0.05u"g/mL" # g solute / mL solution
+ρ_solution = 1u"g/mL" # g/mL total solution density
+R0 = 1.4u"cm^2*hr*Torr/g"
+A1 = 16.0u"cm*hr*Torr/g"
+A2 = 0.0u"1/cm"
+Rp = RpFormFit(R0, A1, A2)
+Vfill = 5u"mL"
+## Heat transfer
+KC = 2.75e-4u"cal/s/K/cm^2"
+KP = 8.93e-4u"cal/s/K/cm^2/Torr"
+KD = 0.46u"1/Torr"
+K_shf_f = RpFormFit(KC, KP, KD)
+## Geometry
+h_f0 = Vfill/A_p
+m_f0 = Vfill * ρ_solution
+## RF fit parameters (dummy values, will be replaced in fitting)
+Bf = 2.0e7u"Ω/m^2"
+Bvw = 0.9e7u"Ω/m^2"
+Kvwf = 10.0u"W/K/m^2"
+## Controllable inputs
+f_RF = 8u"GHz"
+## Total of 10W nominal power, multiplied by 0.54 to account for system losses
+P_per_vial = RampedVariable(10u"W"/17 * 0.54) # actual power / vial
+
+# We combine these into a `ParamObjRF` which will bundle up all the parameters we need.
+# This object needs to be constructed with all these parameters in this order,
+# so it has a constructor taking this tuple-of-tuples form that helps logical grouping.
+params_base = ParamObjRF((
+    (Rp, h_f0, c_solid, ρ_solution),
+    (K_shf_f, A_v, A_p),
+    (pch, Tsh, P_per_vial),
+    (m_f0, LyoPronto.cp_ice, m_v, LyoPronto.cp_gl),
+    (f_RF, LyoPronto.eppf, LyoPronto.epp_gl),
+    (Kvwf, Bf, Bvw),
+))
+# Estimates of some physical properties we will often need are included in LyoPronto,
+# including the heat capacity of ice and glass. A literature correlation for the dielectric
+# loss of ice is provided as `eppf`, and for glass we provide a more-uncertain 
+# single value of `epp_gl`.
+
+# # Parameter fitting
+# First, we set up the fit problem, and double check that our guess values are close-ish.
+## Transform variables to map from a 3-component vector to bounded physical values
+trans_KBB = KBB_transform_bounded(Kvwf, Bf, Bvw)
+## Create nonlinear least-squares function
+nls_M1 = NonlinearFunction{true}(nls_pd!, resid_prototype=zeros(num_errs(fitdat)))
+## Guess values for fit parameters, in log space
+## In practice, these often need tinkering with
+p0 = [3.0, 3.0, 0.3]
+## Check the physical solution for these guess values
+tsol = gen_sol_pd(p0, trans_KBB, params_base)
+modrftplot(tsol)
+plot!(fitdat, nmarks=40)
+
+# Now, we actually run the fit, which is a nonlinear least squares problem
+
+opt1 = solve(NonlinearLeastSquaresProblem(nls_M1, p0, (trans_KBB, params_base, fitdat)), LevenbergMarquardt())
+## Get the fitted solution
+prof_RF = gen_sol_pd(opt1.u, trans_KBB, params_base)
+## Get the fitted parameters in parameter space
+transform(trans_KBB, opt1.u)
+
+# Now, we can plot the fit results against the experimental data.
+
+## Set up the plot
+plT = plot(u"hr", u"°C", xlabel="Time", ylabel="Temperature")
+## Experimental data
+@df thm_pd exptfplot!(:t, :T4, :T1, nmarks=40) 
+@df thm_pd exptvwplot!(:t, :T3, nmarks=40, label=L"$T_\mathrm{vw1}$, exp.")
+## Model results
+modrftplot!(prof_RF, markeralpha=0, trimend=1)
+## Shelf temperature
+plot!(Tsh, c=:black, label=L"T_\mathrm{sh}")
+## End of primary drying
+tendplot!(fitdat.t_end, label="", ls=:dash)
+## Mark the end of drying with a nice label
+annotate!(9, -32, Plots.text("end of drying,\nRF off", 12, "Computer Modern"))
+plot!([fitdat.t_end-1.5u"hr", fitdat.t_end-0.2u"hr"], [-30, -30], arrow=:arrow, c=:black, linewidth=1, label="")
+## Set other plot attributes
+plot!(legend=:topleft, ylim=(-40, 50), xlim=(0,13))
+plot!(size=(600,400), left_margin=15Plots.px)

docs/example/fitting_mannitol.jl

@@ -21,7 +21,7 @@ set_default(labelformat=:square) # Sets a Latexify default
 
 ## Data start at 8th row of CSV file.
 ## This needs to point to the right file, which for documentation is kinda wonky
-procdata_raw = CSV.read(joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv"), Table, header=8)
+procdata_raw = CSV.read(joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv"), Table, header=7)
 t = uconvert.(u"hr", procdata_raw.CycleTime .- procdata_raw.CycleTime[1])
 ## At midnight, timestamps revert to zero, so catch that case
 for i in eachindex(t)[begin+1:end]

docs/example/fitting_rf_mannitol.jl

@@ -3,10 +3,10 @@ using LyoPronto
 using Documenter
 using Literate
 
-@info "Using Literate to generate example"
-Literate.markdown((@__DIR__)*"/example/fitting_mannitol.jl", (@__DIR__)*"/src/generated", documenter=true)
-Literate.markdown((@__DIR__)*"/example/all_recipes.jl", (@__DIR__)*"/src/generated", documenter=true)
-Literate.markdown((@__DIR__)*"/example/utilities.jl", (@__DIR__)*"/src/generated", documenter=true)
+@info "Using Literate to generate examples"
+for file in ["fitting_mannitol.jl", "fitting_rf_mannitol.jl", "all_recipes.jl", "utilities.jl"]
+    Literate.markdown((@__DIR__)*"/example/$file", (@__DIR__)*"/src/generated", documenter=true)
+end
 
 @info "Building Documentation"
 makedocs(;
@@ -17,6 +17,7 @@ makedocs(;
     pages = [
         "Home" => "index.md",
         "Example, conventional lyo" => "generated/fitting_mannitol.md",
+        "Example, microwave-assisted lyo" => "generated/fitting_rf_mannitol.md",
         "Other tools" => "generated/utilities.md",
         "Plot recipes" => "generated/all_recipes.md",
         "Reference" => "alldocstrings.md",

docs/example/utilities.jl

@@ -19,27 +19,44 @@ Some key advantages this has over the original version of LyoPRONTO are:
 ## Installation
 As a Julia package, this code can be easily installed with the Julia package manager. 
 
-From the Julia REPL's Pkg mode (open a REPL and type `]` so that the prompt turns blue), add this package as a Git repo:
+You can add LyoPronto from the Julia General registry (so just like most other packages), using the Julia REPL's Pkg mode (open the REPL and type `]` so the prompt turns blue):
 ```
-add https://github.com/LyoHUB/LyoPronto.jl.git
+add LyoPronto
 ```
 `dev` can be substituted for `add` if you want to make changes to this package yourself, as explained in the [Julia Pkg manual](https://pkgdocs.julialang.org/v1/managing-packages/).
 
 ## Dependencies and Reexports
-This package leverages the strengths of the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) ecosystem to solve equations quickly and efficiently, although it only directly depends on `OrdinaryDiffEqRosenbrock` and `DiffEqCallbacks`, which are both reexported.
 
-Also provided are plot recipes for [Plots.jl](https://docs.juliaplots.org/stable), although this package only depends on `RecipesBase`.
+The following are reexported by this package (so that you don't need to import them after importing LyoPronto):
+- [OrdinaryDiffEqRosenbrock](https://docs.sciml.ai/DiffEqDocs/stable/); used for solving the DAEs and ODEs inherent here
+- [DiffEqCallbacks](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/), used for ending DAE and ODE solves when drying ends
+- [Unitful](https://juliaphysics.github.io/Unitful.jl/stable/); specifically, the `u""` macro, `ustrip`, `uconvert`, and `NoUnits`, which is all the API surface needed for regular usage of LyoPronto.
+
+Other noteworthy dependencies:
+- [TransformVariables.jl](https://tpapp.github.io/TransformVariables.jl/stable/), which is used to map vector spaces onto realistic parameter values for the inevitable parameter fitting step
+- LyoPronto provides [plot recipes](https://docs.juliaplots.org/stable/recipes/) for [Plots.jl](https://docs.juliaplots.org/stable), although it does not depend on Plots in full.  
+
+Plots.jl is used for plotting in this documentation, with the following defaults:
+```@example plot_defaults
+using Plots # hide
+default(:fontfamily, "Computer Modern")
+default(:framestyle, :box)
+default(:lw, 2)
+default(:markersize, 4)
+default(:markerstrokewidth, 0.5)
+default(:unitformat, :square)
+resetfontsizes(); scalefontsizes(1.2)
+```
 
-Heavy use is made of [Unitful.jl](https://painterqubits.github.io/Unitful.jl/stable/), which is reexported.
 
 
 ## Authors
 
-Written by Isaac S. Wheeler, a PhD student at Purdue University.
+Written by Isaac S. Wheeler, a PhD student at Purdue University, advised by Vivek Narsimhan and Alina Alexeenko.
 This work was supported in part by funding for NIIMBL project PC4.1-307 .
 
-## License
+## Licensing
 
-None yet. My intentions are to use the MIT license once this has been published in a scientific journal.
+This package is released with the MIT license.