Bug in class ReferenceIO

CADET Troubleshooting
1

Hi CADET-Devs,

I encountered an issue with line 94 in CADETProcess\reference.py:

solution = np.array(solution, ndmin=2, dtype=np.float64).reshape(len(time), -1) # current code
solution = np.array(solution, ndmin=2, dtype=np.float64).transpose() # my suggestion

While building a minimal reproducable code for my application based on the example process I could not store references for multiple species in the same reference. When I change the code line, it works.

My minimal example
(plots show both created peaks for both components instead of one each):

from examples.batch_elution.process import process
import numpy as np

from CADETProcess.simulator import Cadet
simulator = Cadet()
simulation_result = simulator.simulate(process)
from CADETProcess.reference import ReferenceIO

# gauss peaks as reference data
def gaussian(x, peak_time, peak_height, std_dev):
    return peak_height * np.exp(-((x - peak_time) ** 2) / (2 * std_dev ** 2))

time = np.linspace(0, process.cycle_time, 1000)

data1 = gaussian(time, 6*60, 15, 10)
data2 = gaussian(time, 7.5*60, 5, 30)

reference = ReferenceIO('ref',
                    time, (data1, data2))
reference.component_system = process.component_system

# print references
for ic, comp in enumerate(process.component_system):
    reference.plot(components = str(comp))

edit: code formatting

2

Hi @ortler,

I believe the issue comes from the “wrong” dimension of the solution vector you provide to the ReferenceIO object which expects shape (n_time, n_comp).

In your case,

data1 = gaussian(time, 6*60, 15, 10)
data2 = gaussian(time, 7.5*60, 5, 30)

generates two arrays with shape (1000,). When they are just put in a tuple like you did, this will not crash, however the reshaping will behave unexpectedly.

Instead, you should do the following:

data = np.stack((data1, data1)).T

which has the correct shape (1000, 2).

Admittedly, this is not well documented. Here we will add the following to the docstring of the ReferenceIO class:

solution : array-like
    The reference solution values with shape = (n_time, n_comp).
1 Like
3

Hey Johannes,

thank you for the hint with the right dimensions. It works now to store multiple reference solutions in one reference. :slight_smile:

A subsequent issue appeared while comparing to the simulation results. It works to add all the solutions stored in the reference at once to the comparator, but when I select one component as mentioned here: CADET-Process comparator.add_difference_metric() - CADET Troubleshooting - CADET (cadet-web.de), the comparison fails.

Both reference solutions are still stored in the comparator, I’m not sure what the intended behavior is and where it fails.

from CADETProcess.simulator import Cadet
from CADETProcess.reference import ReferenceIO
from examples.batch_elution.process import process
import numpy as np

simulator = Cadet()
simulation_result = simulator.simulate(process)

# gauss peaks as reference data
def gaussian(x, peak_time, peak_height, std_dev):
    return peak_height * np.exp(-((x - peak_time) ** 2) / (2 * std_dev ** 2))

time = np.linspace(0, process.cycle_time, 1000)
data1 = gaussian(time, 6*60, 15, 10)
data2 = gaussian(time, 7.5*60, 5, 30)
data = np.stack((data1, data2)).T

reference = ReferenceIO('ref', time, data, component_system = process.component_system)

# # print references
# print(process.component_system)
# for ic, comp in enumerate(process.component_system):
#     reference.plot(components = str(comp)) # this works now :-)

from CADETProcess.comparison import Comparator
metric = 'NRMSE'

comparator = Comparator('comp_test')
comparator.add_reference(reference)
comparator.add_difference_metric(
                            metric,
                            reference, # reference for all components
                            'outlet.inlet',
                            components =['B'] # components to be considered
                            )

# print comparators
# comparator.plot_comparison(simulation_result, plot_individual= True) # this fails
# comparator.evaluate(simulation_result) # this fails

# print(comparator.references['ref'].component_system)
# print(comparator.references['ref'].solution)
4

We’ve got a fix. It needs testing but it should be ready in a day or so.

2 Likes
5

The fix is in this pull-request. If you @ortler could confirm for us that this solves the issue we can merge it.

Thank you for bringing this to our attention.

1 Like
6

Hey Ronald,

thank you for the quick fix, it is working now for me! :smiley:

2 Likes