Specifying user parameters for GC SAFT methods

[Garren-H]

Hi all

I just wanted to know if it would be possible to specify name-value pairs for GC SAFT parameters similar to that when using non-GC methods, something like Clapeyron.SAFTgammaMie(compounds, userlocations=(;sigma=updated_pars, ...)), I did not see any examples of these in the docs. I need to define something along these lines to perform parameter estimations on GC methods, utilizing some ML method. I do have an implementation where I define string variables, representitive of csv files, and pass these as the userlocations argument. This does produce output for my objective function but seems to be incompatable with ForwardDiff.jl

Any help will be much appreciated!

[longemen3000]

Hello,

On the parameter inputs, the answer is that we don't have a clean interface to passing GC user parameters. you see, when we create a GC model we do the following things internally:

1. create a GroupParam from the input components field

components = [("propanol",["CH2" => 2, "CH3" => 1, "OH" => 1]), ("water",["H2O" => 1])]
groups = GroupParam(components)

2. the GroupParam creates a list of unique, sorted groups

unique_groups = groups.flattenedgroups #["CH2", "CH3", "OH", "H2O"]

3. we use this list of unique groups to search the database. (including any user parameters). that is, instead of searching the database for components, we search the database for the groups.

What does this mean? you can technically pass GC parameters directly if you know the order of your unique groups beforehand.
If you found this a little confusing, is because it is. I'm really open to suggestions on a possible easier API for this, but if you are fitting parameters within a fixed group set per optimization step. Then using that list of ordered groups instead of components should help you.

On the incompatibility of your function with ForwardDiff, this is expected. The use of parametric number types is restricted (as of now) to PCSAFTparam and SAFTVRMieparam (and all models that use it). With that said, and taking into the account the recent advances on AD through equilibria, (I did not said anything on the last discussion, But what are you trying to do is intrinsically related to #343), we can start adding the support for parametric parameters for more EoS, on a case by case basis. May i ask, what EoS are you using? gcsPCSAFT?

[Garren-H]

Thanks for the fast response, it is truely much appreciated!

What does this mean? you can technically pass GC parameters directly if you know the order of your unique groups beforehand.

This makes sense. So upon knowing the ordering of the groups, we pass either vectors as inputs for SingleParam types and matrices for PairParam types? What about the association sites/params, how would this be specified?

I'm really open to suggestions on a possible easier API for this

A user friendly API would be utilizing either Dicts or Dataframes as inputs. Having one of these would make it easier to specify a full database within the code and updating the relevant fields once per iteration

1. Dicts: For SingleParam types we can specify the inputs as Dict((species)=>val) and for PairParam types we can use Dict((species1, species2)=>val). For association parameters it would probably be best to use the notation currently employed Dict((species1{site1}, species1{site2})=>val) (I saw this is the paper on Clapeyron when viewing model.vrmodel.sites). Sample code would then look something like:

# SingleParam types. Also applicabe when using combining rules for parameters
segment_dict = Dict(("CH2")=>segment_val1,
                    ("CH3")=>segment_val2,
                    ...)

# PairParam types
sigma_dict = Dict(("CH2", "CH2")=>val1,
                  ("CH2", "CH3")=>val2,
                  ("CH3", "CH2")=>val3, # only specify if assymetric params are used, otherwise, assume symmetric, and same as above
                  ...)

# Association parameters. Not sure what the most efficient approach would be, but I would think we can combine the n_sites, either as a vector of integers, or as a named tuple
n_sites = Dict(("OH")=>[1, 2, 0], ...) # Vector of int corresponds to [n_H, n_e1, n_e2]
n_sites = Dict(("OH")=>("H"=>1, "e1"=>2), ...) # specifying the named tuple, we only need the sites which do occur per model

epsilon_assoc_dict = Dict(("OH{H}", "OH{e1}")=>val)
bondvol_dict = Dict(("OH{H}", "OH{e1}")=>val)

# Specifying Dicts as input params
model = SAFTgammaMie(components, userlocations=(;Mw=Mw_dict,
                                                 segment = segment_dict,
                                                 sigma_dict = sigma_dict, 
                                                 ...))

2. DataFrames: We would need to specify two dataframes, one for the like parameters and one for the unlike parameters, keeping the column names consistent with what would be expected in CSV files. Using this there would probably need to be a redefinition of the SAFTgammaMie function to accept these explicit inputs, specifying something like SAFTgammaMie(compounds; userlocations=[df_like, df_unlike], userlocations_is_df=true) (probably not the most efficient and only applicable if DataFrames are typesafe?).

These are just rough ideas which I think would work, but did not think of the implications of this. These approaches would make it easier for a user to specify databases inline and let Clapeyron do the search. If we have a model consisting of OH, CH2 and CH3 groups, we can specify a database consisting of OH, CH2, CH3, NH2 and Clapeyron should search only for the groups which are present, OH, CH2, CH3 (ignoring NH2)

but if you are fitting parameters within a fixed group set per optimization step. Then using that list of ordered groups instead of components should help you.

This does help, but if we want to specify more then one type of property for different pure compounds/mixtures it becomes more difficult. For instance, if we pure compound data for water, ethanol and 1-propanol and we have mixture data for water + ethanol, water + ethanol and water + 1-propanol, we would have to redefine the model each time, by considering the groups within the current model and getting the ordering of the groups correct. If we extend our database to have for instance 100 different compounds, hence having any combination of these mixture data it would result in performance issues. Hence specifying a database once in-line per itteration would be prefereable. We can then just use that database throughout, without having to redefine and match groups to the correct order for each dataset considered.

On the incompatibility of your function with ForwardDiff, this is expected. The use of parametric number types is restricted (as of now) to PCSAFTparam and SAFTVRMieparam

I should've been clear on what I did here. My sample code was as follows for SAFTgammaMie:

unique_groups = # specify unique groups here to use, across all compounds to be considered

# Base csv-like strings to parse parameters
base_csv_like = """Clapeyron Database File
SAFTgammaMie Like Parameters[csvtype = like, grouptype = SAFTgammaMie]
species, vst, S, lambda_r, lambda_a, sigma, epsilon, n_H, n_e1, ne2
""";

base_csv_unlike = """Clapeyron Database File
SAFTgammaMie Unlike Parameters [csvtype = unlike, grouptype = SAFTgammaMie]
species1, species2, lambda_r, sigma, epsilon
""";

function transform_ml_to_model_pars(theta)
   # Perform some computation with theta to get parameter values
   return shape_factor, lambda_r, sigma, epsilon # shape_factor is a vector, all other params are matrices corresponding to the order of unique_groups
end;

function get_csv_files(theta) # uses the base csv-like strings, and generates string databases for the full set of groups considered once
    shape_factor, lambda_r, sigma, epsilon = transform_ml_to_model_pars(theta)
    
    # generate string of `species, vst, S, lambda_r, lambda_a, sigma, epsilon, n_H, n_e1, ne2` for each group and join with \n
    append_like = join(["$(unique_groups[i]), 1, $(shape_factor[i]), $(lambda_r[i,i]), 6, $(sigma[i,i]), $(epsilon[i,i]), 0, 0, 0" for i in 1:length(unique_groups)], "\n")
    # gerate string of `species1, species2, lambda_r, sigma, epsilon` for each combination of unlike groups and join with \n
    append_unlike = join([join(["$(unique_groups[i]), $(unique_groups[j]), $(lambda_r[i,j]), $(sigma[i,j]), $(epsilon[i,j])" for j in i+1:length(unique_groups)], "\n") for i in 1:length(unique_groups)], "\n")

    return [base_csv_like * append_like, base_csv_unlike * append_unlike]
end;

function obj_crit(theta) # train on pure compound critical properties
    csv_files = get_csv_files(theta)
    
    obj = 0.0
    for i in 1:length(comps_of_interest)
        model = Clapeyron.SAFTgammaMie([comps_of_interest[i]], userlocations=csv_files)
        crit_exp = (Tc[i], Pc[i], Vc[i])
        crit_calc = Clapeyron.crit_pure(model)
        obj += sum(((crit_exp .- crit_calc) ./ (0.01 .* crit_exp)).^2)
    end
    return obj
end;

# Add other functions obj functions below, for instance for mixture properties and add obj together to obtain final form

The conversion of the parameters to string types in the get_csv_files function above results in incomptabilities with ForwardDiff, resulting in the error when running ForwardDiff.gradient(obj_crit, theta)

MethodError: no method matching *(::Float64, ::String)
The function `*` exists, but no method is defined for this combination of argument types.

Which arises due to the string concatenation using *. Just running obj_crit(theta) does produce real as output

I did not said anything on the last discussion, But what are you trying to do is intrinsically related to #343

This is closely related to the issue you outline yes, but before the model = build_pcsaft(x) in #343 (comment) we would need to have a function which transforms our database of ML parameters to that interpretable by Clapeyron. I.e. if we have some generic function f(theta) which takes ML parameters theta as input and for instance outputs sigma. The function f(theta) will be different depending on what ML technique we use, so building a generic function to do this will probably not be feasable?

May i ask, what EoS are you using? gcsPCSAFT?

The end goal is to parameterize SAFTgammaMie and sSAFTgammaMie but starting with a simpler model such as gcsPCSAFT would be prefferable as a proof of concept, hence the currently implemented code should be fine? The code structure I would like to implement is as above with my example, or more generically

# Transform ML parameters to SAFT parameters
function transform_ml_pars_to_SAFT_pars(theta)
  # perform calculations
  return SAFT_pars # as vectors or matrices
end

# Update parameter databse only once
function update_param_databases(theta)
   # only update parameter databases once
   SAFT_pars =  transform_ml_pars_to_SAFT_pars(theta) # transform ML parameters to SAFT parameters once
   updated_database = # use SAFT_pars and transform it to a valid interpretable databse for Clapeyron
   return updated_database
end

# Define all objective functions
function obj_crit(theta; updated_database)
   # should loop through all compounds for which data is known for and initialize the model each time
   # return some objective function with critical properties
end

# Pure compound saturation properties
function obj_sat_pure(theta; updated_database)
    ...
end

# Add some obj for mixture properties; probably VLE/VLLE/LLE and liquid excess enthalpies

# final regularized objective functions
function final_reg_obj(theta)
   updated_database = update_param_databases(theta) # update only once
   obj = # add all previous obj, using updated_database as the additional input
   reg_term = # some function of theta for regularization terms
   return obj + reg_term
end

# Perform some optimization using final_reg_obj and FowardDiff

[Garren-H]

Hi all

I did a deep dive into the source code and think I have a way to do this, without changing much to the current code. Ultimately we only want to change the parameters (pure, mixture or assoc), whilst keeping the sites, molecular weights, the decomposition of the compounds into their functional groups and vst. So everything which are not parameters to be tuned, should remain constant in the models, without any redefinitions.

To this end, we can create a new definition for the SAFTgammaMie function (and other GC models?) which allows for the groups and params as input, and hence ignores the current assignment of these functions within the current definition. This would just be copying the existing code into the new definition without the assignment of groups and params:

function SAFTgammaMie(groups::GroupParam, params::Dict{String,ClapeyronParam};
    idealmodel = BasicIdeal,
    ideal_userlocations = String[],
    reference_state = nothing,
    verbose = false,
    epsilon_mixing = :default,
    assoc_options = AssocOptions())

    sites = params["sites"]
    components = groups.components
    ... # remaining code in the existing definition
end

If this is incorporated, for efficiency we can redefine the current defintion to:

function SAFTgammaMie(components;
    idealmodel = BasicIdeal,
    userlocations = String[],
    group_userlocations = String[],
    ideal_userlocations = String[],
    reference_state = nothing,
    verbose = false,
    epsilon_mixing = :default,
    assoc_options = AssocOptions())

    groups = GroupParam(components, ["SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv"]; group_userlocations = group_userlocations,verbose = verbose)
    params = getparams(groups, ["SAFT/SAFTgammaMie","properties/molarmass_groups.csv"]; userlocations = userlocations, verbose = verbose)

    return SAFTgammaMie(groups, params;
                        idealmodel = idealmodel,
                        ideal_userlocations = ideal_userlocations ,
                        reference_state = reference_state,
                        verbose = verbose,
                        epsilon_mixing = epsilon_mixing,
                        assoc_options = assoc_options)
end

We can then generate our own groups (from the csv files) using the GroupParam function and obtain the order of the groups using groups.flattenedgroups field. We can then use the remaining skeleton csv files (potentially only containing headers, or a lot of missing params) to generate the params Dict using the getparams function. If there are missing values for single params in the skeleton csv files, we can just ignore these using the ParamOptions(;ignore_missing_singleparams=...) function.

Some skeleton code is as follows:

using Clapeyron

groulist = ["comp1", "comp2", "comp3", ...]; # compi should be replaced by the compounds considered
groups = GroupParam(grouplist, [csv_to_SAFTgammaMie_groups]); # the default csv_to_SAFTgammaMie_groups = "SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv". The comps should be defined here in terms of their groups

options = ParamOptions(;ignore_missing_singleparams=["sigma"]); # in this case assume sigma is missing. 
params = getparams(groups, [path_to_skeleton_like_unlike_and_assoc_csvs_as_well_as_molar_masses_of_groups],options); # 'vst' should always be known in the like file.

The params is now a Dict with single, pair and assoc pairs. Single and pair params have the fields values and ismissingvalues which are either vectors (single) or matrices(pair). values is type Float, containing the parameters values and ismissingvalues is type Bool, true where missing values are encountered. We can hence overwrite values for single and pair params by using matrices/vectors and updating the ismissingvalues field. If the value for epsilon at position [1,2] is for instance missing, we can do:

params["epsilon"].values[1,2] = params["epsilon"].values[2,1] = some_val; # assuming symmetric values
params["epsilon"].ismissingvalues[1,2] = params["epsilon"].ismissingvalues[2,1] = false; 

Assoc params are a bit more complicated as these are constructed using the compressed4DMatrix. The relevant field to consider in the AssocParam struct are values.values, values.outer_indices and values.inner_indices. values.values is a vector of floats of known pairs of sites. values.outer_indices contains the indices corresponding to the compounds (as ordered in unique_groups) and values.inner_indices contains the indices corresponding to the sites in the different groups; arranged according to the site names in params["sites"].n_sites. That is if params["sites"].n_sites[1] = [e2, e1, H] then the first index of values.inner_indices correspond to e2, the second to e1 and so on. This a bit inconsistent since if we only have H and e1 sites than the order is [H, e1] (but this should not effect the implementation, so it is fine). If we for instance have missing values for group1 site H and group2 site e1 we can add the values to bond_vol and epsilon_assoc. Some sample code:

idx_group1_outer = findfirst(unique_groups .== group1); # find the group index = outer_index
idx_group2_outer = findfirst(unique_groups .== group2); 

idx_group1_inner = findfirst(params["sites"].sites[idx_group1_outer] .== "e1"); # find the site index in the group = inner_index
idx_group2_inner = findfirst(params["sites"].sites[idx_group2_outer] .== "H");

outer_index_app = [(idx_group1_outer, idx_group2_outer)]; # outer index to append
inner_index_app = [(idx_group1_inner, idx_group2_inner)]; # inner index to append

append!(params["epsilon_assoc"].values.values, val_epsilon_assoc);
append!(params["epsilon_assoc"].values.outer_indices, outer_indices_app);
append!(params["epsilon_assoc"].values.inner_indices, inner_indices_app);

append!(params["bondvol"].values.values, val_bondvol); 
append!(params["bondvol"].values.outer_indices, outer_indices_app);
append!(params["bondvol"].values.inner_indices, inner_indices_app);

The final step is then to call the SAFTgammaMie function

model = SAFTgammaMie(groups, params);

This implementation should circumvent the issues associated with autodiff for parameter estimation (should not be limited to SAFT since this is generalizable to other GC approaches). This approach is probably sub-optimal but should be fine for now until a better solution comes up?

Something I am uncertain of however is if the symmetry between sites are maintained in the GC implemetations? That is if we have for instance group1 with [H, e1] sites and group2 with [H,e1] sites, if we assign some value to (group1, "H") >=< (group2, "e1") does Clapeyron automatically assume that (group1, "e1") >=< (group2, "H")? This is a typical assumption when using non-GC methods, so just need some clarity. Any comment here will be much appreciated!

Many thanks
Garren

[longemen3000]

Hello,

Lots to unpack in your comment.

First of all, i've been working on supporting AD on SAFTgammaMie and the last release is almost ready for that. there is already another user of this feature, and we are debugging the final kinks in #366. I mentioned parametric because when calculating a gradient of a function with respect to SAFTgammaMie parameters, you are actually asking the model to hold ForwardDiff.Dual numbers. that was not possible until the last release. So it is a good step towards what you want to do.

On the last release, i separated some fields in SAFTgammaMie that were used in EoS Evaluation (mixed_segment::MixedSegmentGCParam). So there is additional work to do in terms of params, but that shouldn't be user facing at a high level interface.

Now, on the strategy presented here. Seems like a great idea. GroupParam has all the necessary orderings to reorder any parameter in a GC setting. (on a component basis setting, the components field assumes that role.

We already have functions to build skeletons of parameters. For example, Calling SingleParam("paramname",["a","b"]) will create an empty SingleParam{Float64}. same thing with PairParam(name,components) (there is a method missing for AssocParam, but it we can already create Compressed4DMatrix from sites, so that function is easy to add). What we don't have right now, is a way to build skeletons of collections of parameters, and that could be a lot more useful.

On AssocParam, i want to mention that, if the index does not exist, trying to set a new index will fail, but getting an index will result in zero. What we do when we are in need of new AssocParam is creating a "dense" Compressed4DMatrix (via Clapeyron.assoc_similar(sites, T), filling the matrix with the necessary values, and then cleaning the zero values (via SparseArrays.dropzeros!).

On SAFTgammaMie, additionally, we do a transformation of the sites from a group base to a component base. this allows us to use the same association solver for all SAFTs. (that is the main function of the inner vrmodel field.

For most EoS. MyParam(components) should build an skeleton of MyParam. the most general method should something along the lines of:

#pseudocode, we can do this efficiently with a generated function for almost all types, and manually do this for types that don't conform to the API.
function (::Type{P})(compdata::Union{GroupParam,Vector{String},sites::Union{nothing,SiteParam}) where P <: Clapeyron.EoSParam
res = []
comps = compdata isa GroupParam ? : compdata.flattenedgroups : compdata
for name in fieldnames(P)
  f = fieldtype(P,name)
    if f isa AssocParam
      push!(res,f(String(name), sites)
    else
      push!(res,f(String(name), comps )
    end
  end
  return P(res...)
end

function (::Type{P})(compdata::Union{GroupParam,Vector{String}) where P <: Clapeyron.EoSParam
  return P(compdata,nothing)
end

Now, we have a skeleton, and we can fill it with data. we do have some utility functions to do this:

• SingleParam support get/set index with strings (a["CH2"] = 2)
• PairParam also support get/set index with strings (a["CH2","CH3"] = 4, note that the default is symmetric setindex, so for asymmetric models, you should do a["CH2","CH3",false] = 5)
• AssocParam also supports get/set index with strings. the syntax is a little different (param[("H2O","e1"),("H2O","H"),false])

Finally, to answer some questions.

If we assign some value to (group1, "H") >=< (group2, "e1") does Clapeyron automatically assume that (group1, "e1") >=< (group2, "H")

No, we don't assume this in GC nor component-based models, you need to explicitly add both interactions.

Happy to hear what do you think of this.

[Garren-H]

Thanks for this! I did not know that a new version was available. Since I built my installation from the master branch, it didn't show an update was available in the package manager. I have since updated

We already have functions to build skeletons of parameters. For example, Calling SingleParam("paramname",["a","b"]) will create an empty SingleParam{Float64}. same thing with PairParam(name,components) (there is a method missing for AssocParam, but it we can already create Compressed4DMatrix from sites, so that function is easy to add). What we don't have right now, is a way to build skeletons of collections of parameters, and that could be a lot more useful.

I did not know this utility existed. We still need the constant parameters however, such as the molecular weights, vst and the number of sites which should come from csv files. I have a (like) csv file with headers: "species, Mw, vst, n_H, n_e1, n_e2, source". We can use this data to manually construct the SiteParam and use this in the generation of AssocParam. I have written the following code, which performs the required function:

path_of_GC = ["path_to_excel_with_like_values_containing_headers_above.csv"];
groups = GroupParam(my_comps, path_of_GC); # generate the group param
params = getparams(groups, path_of_GC); # generates SingleParam with the headers above

unique_groups = groups.flattenedgroups; # unique groups

mat_num_sites = hcat(params["n_H"].values, params["n_e1"].values, params["n_e2"].values); # extract number of sites in a matrix form
unique_sites = String["H", "e1", "e2"]; # unique sites (always in this order)

sites = [unique_sites[row .> 0] for row in eachrow(mat_num_sites)]; # extract unique sites per species
n_sites = [row[row .> 0] for row in eachrow(mat_num_sites)]; # number of each site per species

# generate SiteParam
site_param = SiteParam(unique_groups,
                       sites,
                       n_sites,
                       String[],
                       nothing);

# construct base compressed_4D_matrix
assoc_compressed_4D = Clapeyron.assoc_similar(site_param, Float64);

# Assign non-zero values for the outer and inner indices considered. Will change to automate based on the outer and inner indices
assoc_compressed_4D[1,1][1,2] = 1;
SparseArrays.dropzeros!(assoc_compressed_4D); # drop all non-zero values, such that we can form a matrix

# Function to transform model pars to SAFT pars
function transform_ML_pars_to_SAFT(x)
    S = # vector
    sigma = # Vector
    epsilon = # Matrix
    lambda_a = # Matrix; could be set to 6.0 and excluded
    lambda_r = # Matrix
    epsilon_assoc = # Flattened vector; assuming correct odering of site pairs
    bondvol = # Flattened vector; assuming correct ordering of sites pairs
    return S, sigma, epsilon, lambda_a, lambda_r, epsilon_assoc, bondvol
end;

# function to generate Dict{String, Clapeyron.Param}
function transform_SAFT_to_params(x)
    S, sigma, epsilon, lambda_a, lambda_r, epsilon_assoc, bondvol = transform_ML_pars_to_SAFT(x)
    
    # can probably automate this, but works for now
    S_param = SingleParam("S", unique_groups, S);
    sigma_param = SingleParam("sigma", unique_groups, sigma);
    epsilon_param = PairParam("epsilon", unique_groups, epsilon);
    lambda_a_param = PairParam("lambda_a", unique_groups, lambda_a);
    lambda_r_param = PairParam("lambda_r", unique_groups, lambda_r);

    # Use deepcopy to adjust compressed4D for each AssocParam. Assuming correct ordering of site pairs
    epsilon_assoc_param = AssocParam("epsilon_assoc", unique_groups, deepcopy(assoc_compressed_4D), sites)
    epsilon_assoc_param.values.values .= epsilon_assoc
    bondvol_param = AssocParam("bondvol", unique_groups, deepcopy(assoc_compressed_4D), sites)
    bondvol_param.values.values. = bondvol
    
    Clap_pars = Dict([("Mw",params["Mw"]), # Mw from original csv sheet
                    ("vst",params["vst"]), # from original csv sheet
                    ("sites",site_param), # generated site_param
                    ("S",S_param),
                    ("sigma",sigma_param),
                    ("epsilon",epsilon_param),
                    ("lambda_a",lambda_a_param),
                    ("lambda_r",lambda_r_param),
                    ("epsilon_assoc",epsilon_assoc_param),
                    ("bondvol",bondvol_param)])

    return Clap_pars
end;

We can then test this using SAFTgammaMie(groups, Clap_pars(x)) using the new definition, which works on my side

• SingleParam support get/set index with strings (a["CH2"] = 2)
• PairParam also support get/set index with strings (a["CH2","CH3"] = 4, note that the default is symmetric setindex, so for asymmetric models, you should do a["CH2","CH3",false] = 5)
• AssocParam also supports get/set index with strings. the syntax is a little different (param[("H2O","e1"),("H2O","H"),false])

No, we don't assume this in GC nor component-based models, you need to explicitly add both interactions.

Noted, thanks!

First of all, i've been working on supporting AD on SAFTgammaMie and the last release is almost ready for that. there is already another user of this feature, and we are debugging the final kinks in #366.

I shall have a look at this thread. The only thing which doesn't work now is AD (with v0.6.12 currently on master)

Also, this methods breaks when specifying Real as the eltype. I have traced this back to several functions. I shall raise an issue for this

[longemen3000]

Hello,

I'm interested in the errors, you experience with AD, please open an issue for those (with a minimum reproducible example).

on the Real eltype, that is known in julia as an abstract type, it could affect performance and cause unexpected issues.

Finally, it seems like we are missing some convenience methods for SiteParam, I'm gonna check what can i do.

[Garren-H]

I have opened this issue #369