optimeo.doe

  1# Copyright (c) 2025 Colin BOUSIGE
  2# Contact: colin.bousige@cnrs.fr
  3#
  4# This program is free software: you can redistribute it and/or modify
  5# it under the terms of the MIT License as published by
  6# the Free Software Foundation, either version 3 of the License, or
  7# any later version. 
  8
  9"""
 10This module provides a class for creating and visualizing a design of experiments (DoE).
 11It supports various types of designs including Full Factorial, Sobol sequence, Fractional Factorial, Definitive Screening Design, Space Filling Latin Hypercube, Randomized Latin Hypercube, Optimal, Plackett-Burman, and Box-Behnken.
 12The class allows the user to specify the parameters, their types, and values, and generates the design accordingly.
 13It also provides a method to plot the design using scatter plots.
 14
 15You can see an example notebook [here](../examples/doe.html).
 16"""
 17
 18
 19import warnings
 20warnings.simplefilter(action='ignore', category=FutureWarning)
 21warnings.simplefilter(action='ignore', category=DeprecationWarning)
 22warnings.simplefilter(action='ignore', category=UserWarning)
 23warnings.simplefilter(action='ignore', category=RuntimeError)
 24import numpy as np
 25import pandas as pd
 26from typing import Any, Dict, List, Optional
 27from dexpy.optimal import build_optimal
 28from dexpy.model import ModelOrder
 29from dexpy.design import coded_to_actual
 30from doepy import build
 31from sklearn.preprocessing import LabelEncoder
 32from pyDOE3 import *
 33import definitive_screening_design as dsd
 34import plotly.express as px
 35from itertools import combinations
 36import plotly.graph_objects as go
 37
 38class DesignOfExperiments:
 39    """
 40    Class to create a design of experiments (DoE) for a given model.
 41    This class allows the user to specify the type of design, the parameters,
 42    and various options for the design generation.
 43    The design can be visualized using scatter plots.
 44    
 45    
 46    Parameters
 47    ----------
 48    
 49    type : str
 50        The type of design to create. Must be one of:
 51        `'Full Factorial'`, `'Sobol sequence'`, `'Fractional Factorial'`,
 52        `'Definitive Screening'`, `'Space Filling Latin Hypercube'`,
 53        `'Randomized Latin Hypercube'`, `'Optimal'`, `'Plackett-Burman'`,
 54        `'Box-Behnken'` or `'Central Composite'`.
 55    parameters : List[Dict[str, Dict[str, Any]]]
 56        List of parameters for the design, each with a dictionary of properties.
 57        Each dictionary should contain 'name', 'type', and 'values'.
 58        'values' should be a list of possible values for the parameter.
 59        'type' should be either "int", "integer", "float", "<other>". 
 60        Any <other> will be considered as "categorical".
 61        'values' should be a list of possible values for the parameter.
 62    Nexp : int, optional
 63        Number of experiments in the design, when applicable. Default is 4.
 64    order : int, optional
 65        Order of the model (for 'Optimal' design). Default is 2.
 66    randomize : bool, optional
 67        Whether to randomize the run order. Default is True.
 68    reduction : int, optional
 69        Reduction factor for 'Fractional Factorial' designs. Default is 2.
 70    feature_constraints : Optional[List[Dict[str, Any]]], optional
 71        Feature constraints of the experiment for Sobol sequence. Default is None.
 72        If a single dictionary is provided, it will be converted to a list.
 73        If a string is provided, it will be converted to a list with one element.
 74        If a list is provided, it will be used as is.
 75        If None, no constraints will be applied.
 76    
 77    Attributes
 78    ----------
 79    
 80    type : str
 81        The type of design.
 82    parameters : List[Dict[str, Dict[str, Any]]]
 83        The parameters for the design.
 84    Nexp : int
 85        Number of experiments in the design.
 86    order : int
 87        Order of the model.
 88    randomize : bool
 89        Whether to randomize the run order.
 90    reduction : int
 91        Reduction factor for `'Fractional Factorial'` designs.
 92    design : pd.DataFrame
 93        The design DataFrame.
 94    lows : Dict[str, float]
 95        Lower bounds for the parameters.
 96    highs : Dict[str, float]
 97        Upper bounds for the parameters.
 98    
 99    Methods
100    -------
101    - **create_design()**:
102        Create the design of experiments based on the specified type and parameters.
103    - **plot()**:
104        Plot the design of experiments using plotly.
105    
106    
107    Example
108    -------
109    
110    ```python
111    from doe import DesignOfExperiments
112    parameters = [
113        {'name': 'Temperature', 'type': 'integer', 'values': [20, 30, 40]},
114        {'name': 'Pressure', 'type': 'float', 'values': [1, 2, 3]},
115        {'name': 'Catalyst', 'type': 'categorical', 'values': ['A', 'B', 'C']}
116    ]
117    doe = DesignOfExperiments(
118        type='Full Factorial',
119        parameters=parameters
120    )
121    design = doe.design
122    print(design)
123    figs = doe.plot()
124    for fig in figs:
125        fig.show()
126    ```
127    
128
129    """
130
131    def __init__(self, 
132                 type: str,
133                 parameters: List[Dict[str, Dict[str, Any]]],
134                 Nexp: int = 4, 
135                 order: int = 2, 
136                 randomize: bool = True, 
137                 reduction: int = 2,
138                 feature_constraints: Optional[List[Dict[str, Any]]] = None,
139                 center=(2,2),
140                 alpha='o',
141                 face='ccc'):
142        self.type = type
143        self.parameters = parameters
144        self.Nexp = Nexp
145        self.order = order
146        self.randomize = randomize
147        self.reduction = reduction
148        self.center = center
149        self.alpha = alpha
150        self.face = face
151        self.design = None
152        self.lows = {}
153        self.feature_constraints = feature_constraints
154        self.highs = {}
155        self.create_design()
156
157    def __repr__(self):
158        return self.__str__()
159
160    def __str__(self):
161        """
162        Return a string representation of the DesignOfExperiments instance.
163        """
164        printpar = "\n".join([str(par) for par in self.parameters])
165        return f"""
166- Design of Experiments type: {self.type}
167- Parameters:
168{printpar}
169- Lows: {self._lows}
170- Highs: {self._highs}
171- If applicable:
172    - Randomize: {self.randomize}
173    - Number of Experiments: {self.Nexp}
174    - Order: {self.order}
175    - Reduction: {self.reduction}
176- Design:
177{self.design}
178"""
179
180    @property
181    def type(self) -> str:
182        """The type of design to create. Must be one of: `'Full Factorial'`, `'Sobol sequence'`, `'Fractional Factorial'`, `'Definitive Screening'`, `'Space Filling Latin Hypercube'`, `'Randomized Latin Hypercube'`, `'Optimal'`, `'Plackett-Burman'`, `'Box-Behnken'` or `'Central Composite'`."""
183        return self._type
184
185    @type.setter
186    def type(self, value: str):
187        """Set the type of design."""
188        self._type = value
189
190    @property
191    def parameters(self) -> List[Dict[str, Dict[str, Any]]]:
192        """List of parameters for the design, each with a dictionary of properties.
193            Each dictionary should contain the keys `"name"`, `"type"`, and `"values"`.
194            `"values"` should be a list of possible values for the parameter.
195            `"type"` should be either `"int"`, `"integer"`, `"float"`, `"<other>"`. 
196            Any `"<other>"` will be considered as `"categorical"`.
197            `values` should be a list of possible values for the parameter."""
198        return self._parameters
199
200    @parameters.setter
201    def parameters(self, value: List[Dict[str, Dict[str, Any]]]):
202        """Set the parameters for the design."""
203        self._parameters = value
204
205    @property
206    def Nexp(self) -> int:
207        """Number of experiments in the design, when applicable. Default is `4`."""
208        return self._Nexp
209
210    @Nexp.setter
211    def Nexp(self, value: int):
212        """Set the number of experiments."""
213        self._Nexp = value
214
215    @property
216    def order(self) -> int:
217        """Order of the model (for `'Optimal'` design). Default is `2`."""
218        return self._order
219    
220    @property
221    def center(self) -> tuple:
222        """Center for the Central Composite Design. Must be a tuple of two values."""
223        return self._center
224    
225    @center.setter
226    def center(self, value: tuple):
227        """Set the center of the design."""
228        if not isinstance(value, tuple):
229            raise ValueError("Center must be a tuple of two values.")
230        if len(value) != 2:
231            raise ValueError("Center must be a tuple of two values.")
232        if not all(isinstance(i, (int, float)) for i in value):
233            raise ValueError("Center must be a tuple of two numeric values.")
234        self._center = value
235    
236    @property
237    def alpha(self) -> str:
238        """Alpha for the Central Composite Design. Default is `'o'` (orthogonal).
239        Can be either `'o'` or `'r'` (rotatable)."""
240        return self._alpha
241    
242    @alpha.setter
243    def alpha(self, value: str):
244        """Set the alpha of the design."""
245        if value not in ['o', 'r']:
246            raise ValueError("Alpha must be either 'o' (orthogonal) or 'r' (rotatable).")
247        self._alpha = value
248    
249    @property
250    def face(self) -> str:
251        """The relation between the start points and the corner (factorial) points for the Central Composite Design.
252        
253        There are three possible options for this input:
254        
255        1. 'circumscribed' or 'ccc' (Default)
256        2. 'inscribed' or 'cci'
257        3. 'faced' or 'ccf'"""
258        return self._face
259
260    @face.setter
261    def face(self, value: str):
262        """Set the face of the design."""
263        if value not in ['ccc', 'cci', 'ccf']:
264            raise ValueError("Face must be either 'ccc' (circumscribed), 'cci' (inscribed), or 'ccf' (faced).")
265        self._face = value
266    
267    @property
268    def lows(self) -> Dict[str, float]:
269        """Get the lower bounds for the parameters."""
270        return self._lows
271    
272    @lows.setter
273    def lows(self, value: Dict[str, float]):
274        """Set the lower bounds for the parameters."""
275        self._lows = value
276    
277    @property
278    def highs(self) -> Dict[str, float]:
279        """Get the upper bounds for the parameters."""
280        return self._highs
281    
282    @highs.setter
283    def highs(self, value: Dict[str, float]):
284        """Set the upper bounds for the parameters."""
285        self._highs = value
286
287    @order.setter
288    def order(self, value: int):
289        """Set the order of the model."""
290        self._order = value
291
292    @property
293    def randomize(self) -> bool:
294        """Whether to randomize the run order. Default is `True`."""
295        return self._randomize
296
297    @randomize.setter
298    def randomize(self, value: bool):
299        """Set the randomize flag."""
300        self._randomize = value
301
302    @property
303    def reduction(self) -> int:
304        """Reduction factor for `'Fractional Factorial'` designs. Default is `2`."""
305        return self._reduction
306
307    @reduction.setter
308    def reduction(self, value: int):
309        """Set the reduction factor."""
310        self._reduction = value
311
312    @property
313    def design(self) -> pd.DataFrame:
314        """Get the design DataFrame."""
315        return self._design
316
317    @design.setter
318    def design(self, value: pd.DataFrame):
319        """Set the design DataFrame."""
320        self._design = value
321    
322    @property
323    def feature_constraints(self):
324        """
325        Get the feature constraints of the experiment for Sobol sequence.
326        """
327        return self._feature_constraints
328
329    @feature_constraints.setter
330    def feature_constraints(self, value):
331        """
332        Set the feature constraints of the experiment with validation for Sobol sequence.
333        """
334        if isinstance(value, dict):
335            self._feature_constraints = [value]
336        elif isinstance(value, list):
337            self._feature_constraints = value if len(value) > 0 else None
338        elif isinstance(value, str):
339            self._feature_constraints = [value]
340        else:
341            self._feature_constraints = None
342
343# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
344# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
345# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
346
347    def create_design(self):
348        """
349        Create the design of experiments based on the specified type and parameters.
350        """
351        for par in self.parameters:
352            if par['type'].lower() == "categorical":
353                if self.type != 'Sobol sequence':
354                    le = LabelEncoder()
355                    label = le.fit_transform(par['values'])
356                    par['values'] = label
357                    par['encoder'] = le
358                    self.lows[par['name']] = np.min(par['values'])
359                    self.highs[par['name']] = np.max(par['values'])
360                else:
361                    par['encoder'] = None
362            else:
363                self.lows[par['name']] = np.min(par['values'])
364                self.highs[par['name']] = np.max(par['values'])
365
366        pars = {par['name']: par['values'] for par in self.parameters}
367
368        if self.type == 'Full Factorial':
369            self.design = build.full_fact(pars)
370        elif self.type == 'Sobol sequence':
371            from ax.modelbridge.generation_strategy import GenerationStep, GenerationStrategy
372            from ax.modelbridge.registry import Models
373            from ax.service.ax_client import AxClient, ObjectiveProperties
374
375            ax_client = AxClient()
376            params = []
377            for par in self.parameters:
378                if par['type'].lower() == "float":
379                    params.append({'name': par['name'],
380                                   'type': 'range',
381                                   'value_type': 'float',
382                                   'bounds': [float(np.min(par['values'])), float(np.max(par['values']))]})
383                elif par['type'].lower() in ["integer", 'int']:
384                    params.append({'name': par['name'],
385                                   'type': 'range',
386                                   'value_type': 'int',
387                                   'bounds': [int(np.min(par['values'])), int(np.max(par['values']))]})
388                else:
389                    params.append({'name': par['name'],
390                                   'type': 'choice',
391                                   'values': par['values']})
392            ax_client.create_experiment(
393                name="DOE",
394                parameters=params,
395                objectives={"response": ObjectiveProperties(minimize=False)},
396                parameter_constraints=self._feature_constraints
397            )
398            gs = GenerationStrategy(
399                steps=[GenerationStep(
400                    model=Models.SOBOL,
401                    num_trials=-1,
402                    should_deduplicate=True,
403                    model_kwargs={"seed": 165478},
404                    model_gen_kwargs={},
405                )]
406            )
407            generator_run = gs.gen(
408                experiment=ax_client.experiment,
409                data=None,
410                n=self.Nexp
411            )
412            if self.Nexp == 1:
413                ax_client.experiment.new_trial(generator_run)
414            else:
415                ax_client.experiment.new_batch_trial(generator_run)
416            trials = ax_client.get_trials_data_frame()
417            self.design = trials[trials['trial_status'] == 'CANDIDATE']
418            self.design = self._design.drop(columns=['trial_index',
419                                                      'trial_status',
420                                                      'arm_name',
421                                                      'generation_method',
422                                                      'generation_node'])
423        elif self.type == 'Fractional Factorial':
424            for par in range(len(self.parameters)):
425                if self.parameters[par]['type'] == "Numerical":
426                    self.parameters[par]['type'] = "Categorical"
427                    le = LabelEncoder()
428                    label = le.fit_transform(self.parameters[par]['values'])
429                    self.parameters[par]['values'] = label
430                    self.parameters[par]['encoder'] = le
431            design = gsd([len(par['values']) for par in self.parameters], self.reduction)
432            self.design = pd.DataFrame(design, columns=[par['name'] for par in self.parameters])
433        elif self.type == 'Definitive Screening':
434            params = {par['name']: [np.min(par['values']), np.max(par['values'])] for par in self.parameters}
435            self.design = dsd.generate(factors_dict=params)
436        elif self.type == 'Space Filling Latin Hypercube':
437            self.design = build.space_filling_lhs(pars, num_samples=self.Nexp)
438        elif self.type == 'Randomized Latin Hypercube':
439            self.design = build.lhs(pars, num_samples=self.Nexp)
440        elif self.type == 'Optimal':
441            reaction_design = build_optimal(
442                len(self.parameters),
443                order=ModelOrder(self.order),
444                run_count=self.Nexp)
445            reaction_design.columns = [par['name'] for par in self.parameters]
446            self.design = coded_to_actual(reaction_design, self._lows, self._highs)
447        elif self.type == 'Plackett-Burman':
448            self.design = build.plackett_burman(pars)
449        elif self.type == 'Box-Behnken':
450            if len(self.parameters) < 3 or any([len(par['values']) < 3 for par in self.parameters]):
451                self.design = pd.DataFrame({})
452                raise Warning("Box-Behnken design is not possible with less than 3 parameters and with less than 3 levels for any parameter.")
453            else:
454                self.design = build.box_behnken(d=pars)
455        elif self.type == 'Central Composite':
456            self.design = build.central_composite(pars, 
457                                                  center=self.center,
458                                                  alpha=self.alpha,
459                                                  face=self.face)
460        else:
461            raise Warning("Unknown design type. Must be one of: 'Full Factorial', 'Sobol sequence', 'Fractional Factorial', 'Definitive Screening', 'Space Filling Latin Hypercube', 'Randomized Latin Hypercube', 'Optimal', 'Plackett-Burman', 'Box-Behnken' or 'Central Composite'.")
462
463        for par in self.parameters:
464            if par['type'] == "Categorical" and self.type != 'Sobol sequence':
465                vals = self._design[par['name']].to_numpy()
466                self.design[par['name']] = par['encoder'].inverse_transform([int(v) for v in vals])
467
468        # randomize the run order
469        self.design['run_order'] = np.arange(len(self._design)) + 1
470        if self.randomize:
471            ord = self._design['run_order'].to_numpy()
472            self.design['run_order'] = np.random.permutation(ord)
473        cols = self._design.columns.tolist()
474        cols = cols[-1:] + cols[:-1]
475        self.design = self._design[cols]
476        # apply the column types
477        for col in self._design.columns:
478            for par in self.parameters:
479                if col == par['name']:
480                    if par['type'].lower() == "float":
481                        self.design[col] = self._design[col].astype(float)
482                    elif par['type'].lower() in ["int", "integer"]:
483                        self.design[col] = self._design[col].astype(int)
484                    else:
485                        self.design[col] = self._design[col].astype(str)
486        return self._design
487
488    def plot(self):
489        """
490        Plot the design of experiments.
491
492        Returns
493        -------
494        List of plotly.graph_objs._figure.Figure
495            A list of Plotly figures representing the design of experiments.
496        """
497        fig = []
498        count = 0
499        if len(self.design) > 0:
500            if len(self.parameters) <= 2:
501                # Create 2D scatter plots
502                for i, faci in enumerate(self.parameters):
503                    for j, facj in enumerate(self.parameters):
504                        if j > i:
505                            fig.append(px.scatter(
506                                self.design,
507                                x=facj['name'],
508                                y=faci['name'],
509                                title=f"""{faci['name']} vs {facj['name']}""",
510                                labels={facj['name']: facj['name'], faci['name']: faci['name']}
511                            ))
512                            fig[count].update_traces(marker=dict(size=10))
513                            fig[count].update_layout(
514                                margin=dict(l=10, r=10, t=50, b=50),
515                                xaxis=dict(
516                                    showgrid=True,
517                                    gridcolor="lightgray",
518                                    zeroline=False,
519                                    showline=True,
520                                    linewidth=1,
521                                    linecolor="black",
522                                    mirror=True
523                                ),
524                                yaxis=dict(
525                                    showgrid=True,
526                                    gridcolor="lightgray",
527                                    zeroline=False,
528                                    showline=True,
529                                    linewidth=1,
530                                    linecolor="black",
531                                    mirror=True
532                                ),
533                            )
534                            count += 1
535            else:
536                # Create 3D scatter plots
537                for k, (faci, facj, fack) in enumerate(combinations(self.parameters, 3)):
538                    fig.append(go.Figure(data=[go.Scatter3d(
539                        x=self.design[facj['name']],
540                        y=self.design[faci['name']],
541                        z=self.design[fack['name']],
542                        mode='markers',
543                        marker=dict(size=10, color='royalblue', opacity=0.7),
544                    )]))
545                    fig[count].update_layout(
546                        template='ggplot2',
547                        height=500,
548                        width=500,
549                        scene=dict(
550                            xaxis_title=facj['name'],
551                            yaxis_title=faci['name'],
552                            zaxis_title=fack['name'],
553                        ),
554                        title=f"{faci['name']} vs {facj['name']}<br>vs {fack['name']}",
555                        margin=dict(l=10, r=10, t=50, b=50)
556                    )
557                    count += 1
558        return fig