Repository: dugucrypter/georunes (language: python)

Issue:
Dev



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index f7aac7e..249296b 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,7 @@ This code is written in **Python 3.6**, and depends on **matplotlib 3.3**, **pan
 
 ```python
 import matplotlib.pyplot as plt
-from georunes.binary.versus import DiagramVs
+from georunes.plot.binary.versus import DiagramVs
 
 # Initialize diagram class
 diag_nb_ta = DiagramVs(datasource="path/to/data.xls", sheet="sheet1",
diff --git a/examples/WAC_granitoids.csv b/examples/diagrams/WAC_granitoids.csv
similarity index 100%
rename from examples/WAC_granitoids.csv
rename to examples/diagrams/WAC_granitoids.csv
diff --git a/examples/WAC_granitoids_comp.xls b/examples/diagrams/WAC_granitoids_comp.xls
similarity index 99%
rename from examples/WAC_granitoids_comp.xls
rename to examples/diagrams/WAC_granitoids_comp.xls
index 915e9ba..6a430d3 100644
Binary files a/examples/WAC_granitoids_comp.xls and b/examples/diagrams/WAC_granitoids_comp.xls differ
diff --git a/examples/binary_diagrams.py b/examples/diagrams/binary_diagrams.py
similarity index 83%
rename from examples/binary_diagrams.py
rename to examples/diagrams/binary_diagrams.py
index 8751c13..bdb32e0 100644
--- a/examples/binary_diagrams.py
+++ b/examples/diagrams/binary_diagrams.py
@@ -1,10 +1,12 @@
 import os
+
 import matplotlib.pyplot as plt
-from georunes.binary.ratiosversus import DiagramRatiosVs
-from georunes.binary.sialkali import DiagramSiAlkali
-from georunes.binary.versus import DiagramVs
+
+from georunes.plot.binary.pearce import DiagramPearceRYN
+from georunes.plot.binary.ratiosversus import DiagramRatiosVs
+from georunes.plot.binary.sialkali import DiagramSiAlkali
+from georunes.plot.binary.versus import DiagramVs
 from georunes.tools.language import format_chemical_formula
-from georunes.binary.pearce import DiagramPearceRYN
 
 directory = os.path.dirname(os.path.realpath("WAC_granitoids_comp.xls"))
 source = os.path.join(directory, "WAC_granitoids_comp.xls")
@@ -43,8 +45,8 @@
                                sheet=sheet,
                                group_name='Lithology',
                                legend_ncol=4,
-                               padding={"bottom":0.17}
-)
+                               padding={"bottom": 0.17}
+                               )
 
 cox.plot()
 ntzh.plot()
diff --git a/examples/preprocess_files.py b/examples/diagrams/preprocess_files.py
similarity index 97%
rename from examples/preprocess_files.py
rename to examples/diagrams/preprocess_files.py
index a0a3adc..fe1ee69 100644
--- a/examples/preprocess_files.py
+++ b/examples/diagrams/preprocess_files.py
@@ -1,6 +1,6 @@
 import os
 import matplotlib.pyplot as plt
-from georunes.binary.versus import DiagramVs
+from georunes.plot.binary.versus import DiagramVs
 from georunes.tools.preprocessing import populate_file, populate_file_with_config
 from georunes.tools.language import format_chemical_formula as _fml
 
diff --git a/examples/spider_diagrams.py b/examples/diagrams/spider_diagrams.py
similarity index 97%
rename from examples/spider_diagrams.py
rename to examples/diagrams/spider_diagrams.py
index b981979..6baee65 100644
--- a/examples/spider_diagrams.py
+++ b/examples/diagrams/spider_diagrams.py
@@ -1,6 +1,6 @@
 import os
 import matplotlib.pyplot as plt
-from georunes.base.spider import DiagramSpider
+from georunes.plot.spider import DiagramSpider
 
 # 'WAC_granitoids_pop2.csv' is created from WAC_granitoids.csv in the script preprocess_files.py
 directory = os.path.dirname(os.path.realpath("WAC_granitoids_pop2.csv"))
diff --git a/examples/ternary_diagrams.py b/examples/diagrams/ternary_diagrams.py
similarity index 85%
rename from examples/ternary_diagrams.py
rename to examples/diagrams/ternary_diagrams.py
index 3a92eae..105ba46 100644
--- a/examples/ternary_diagrams.py
+++ b/examples/diagrams/ternary_diagrams.py
@@ -1,8 +1,8 @@
 import os
 import matplotlib.pyplot as plt
-from georunes.ternary.anabor import DiagramAnAbOr
-from georunes.ternary.jensen import DiagramJensen
-from georunes.ternary.qap import DiagramQAP
+from georunes.plot.ternary.anabor import DiagramAnAbOr
+from georunes.plot.ternary.jensen import DiagramJensen
+from georunes.plot.ternary.qap import DiagramQAP
 
 # 'WAC_granitoids_pop2.csv' is created from WAC_granitoids.csv in the script preprocess_files.py
 directory = os.path.dirname(os.path.realpath("WAC_granitoids_pop2.csv"))
diff --git a/examples/modal mineralogy/minerals.csv b/examples/modal mineralogy/minerals.csv
new file mode 100644
index 0000000..56ba54f
--- /dev/null
+++ b/examples/modal mineralogy/minerals.csv
@@ -0,0 +1,9 @@
+Mineral,SiO2,TiO2,Al2O3,Fe2O3,MnO,MgO,CaO,Na2O,K2O,P2O5,LOI
+Albite,68.74,,19.44,,,,,11.82,,,
+Anorthite,43.19,,36.65,,,,20.16,,,,
+Quartz,100,,,,,,,,,,
+Hematite,,,,100,,,,,,,
+Biotite,36.42,2.39,15.4,20.76,0.33,10.85,,,9.44,,4.41
+Muscovite,46.32,0.82,31.87,5.66,,1.01,,0.37,10.48,,3.47
+Orthoclase,64.76,,18.32,,,,,,16.92,,
+Apatite-F,,,,,,,55.6,,,42.22,2.18
\ No newline at end of file
diff --git a/georunes/base/__init__.py b/georunes/modmin/__init__.py
similarity index 100%
rename from georunes/base/__init__.py
rename to georunes/modmin/__init__.py
diff --git a/georunes/binary/__init__.py b/georunes/modmin/norm/__init__.py
similarity index 100%
rename from georunes/binary/__init__.py
rename to georunes/modmin/norm/__init__.py
diff --git a/georunes/modmin/norm/cipw.py b/georunes/modmin/norm/cipw.py
new file mode 100644
index 0000000..358bf38
--- /dev/null
+++ b/georunes/modmin/norm/cipw.py
@@ -0,0 +1,824 @@
+import warnings
+import numpy as np
+from pandas import Series
+from georunes.tools.chemistry import ratio_el_to_ox, molar_mass, el_molar_mass
+from georunes.modmin.optim.base import BaseOptimizer
+
+list_major_ox = ['SiO2', 'Na2O', 'K2O', 'CaO', 'MnO', 'FeO', 'Fe2O3', 'MgO', 'TiO2', 'Al2O3', 'P2O5']
+list_minor_el = ['F', 'Cl', 'S', 'Ni', 'Co', 'Ba', 'Sr', 'Rb', 'Cs', 'Li', 'Zr', 'Cr', 'V']
+list_ox_from_minor_el = ['SO3', 'NiO', 'CoO', 'BaO', 'SrO', 'Rb2O', 'Cs2O', 'Li2O', 'ZrO2', 'Cr2O3', 'V2O3']
+list_all_ox_el = [*list_major_ox, 'CO2', 'F', 'Cl', 'S', *list_ox_from_minor_el]
+final_list_min = ['Q', 'Or', 'Ab', 'An', 'C', 'Wo', 'Cs', 'Tn', 'Pf', 'Cc', 'Mt', 'Hm', 'Il', 'Ru',
+                  'Ne', 'Th', 'Nc', 'Ac', 'Ns', 'Lc', 'Kp', 'Ks', 'Z', 'Hl', 'Fr', 'Pr', 'Cm']
+
+oxide_corresp = {
+    'Ni': 'NiO',
+    'Co': 'CoO',
+    'Ba': 'BaO',
+    'Sr': 'SrO',
+    'Rb': 'Rb2O',
+    'Cs': 'Cs2O',
+    'Li': 'Li2O',
+    'Zr': 'ZrO2',
+    'Cr': 'Cr2O3',
+    'V': 'V2O3'
+}
+
+_mol_w = dict()
+for _oxel in list_all_ox_el:
+    if _oxel in molar_mass:
+        _mol_w[_oxel] = molar_mass[_oxel]
+    elif _oxel in el_molar_mass:
+        _mol_w[_oxel] = el_molar_mass[_oxel]
+for _oxel in ['O', 'Cl', 'S']:
+    _mol_w[_oxel] = el_molar_mass[_oxel]
+
+# Algorithm updated from Verma et al., 2003
+# Verma, S.P., Torres-Alvarado, I.S. and Velasco-Tapia, F., 2003. A revised CIPW norm. Swiss Bulletin of Mineralogy
+# and Petrology, 83(2), pp.197-216.
+
+class CIPWNorm(BaseOptimizer):
+    def __init__(self, **kwargs):
+        BaseOptimizer.__init__(self, **kwargs)
+        self.notif = ">>>>>> CIPW norm"
+
+    def compute(self, raw_data, skip_cols, normalize_entry=False, minor_included=False, to_round=4,
+                co2_cancrinite=False, co2_calcite=False):
+
+        # 1 / Preparing data
+        if self.verbose > 1: print("Step 1 - Prepare data")
+
+        co2_cancrinite = float(co2_cancrinite)
+        co2_calcite = float(co2_calcite)
+        data = raw_data.iloc[:, :skip_cols].copy()
+        raw_data_keys = raw_data.iloc[:, skip_cols:].keys()
+
+        for oxel in raw_data_keys:
+            if oxel in [*list_all_ox_el, *list_minor_el]:
+                data[oxel] = raw_data[oxel]
+
+        if 'Total' in raw_data_keys or 'Sum' in raw_data_keys:
+            data['Total'] = raw_data['Total']
+        else:
+            warnings.warn("Total column not found.")
+
+        data_keys = data.keys()
+        filtered_columns = [oxel for oxel in raw_data_keys if oxel not in data_keys]
+        if self.verbose: print('Ignored columns :', filtered_columns)
+
+        for ox in list_major_ox:
+            if ox not in data_keys:
+                data[ox] = 0
+
+        n_phase = data.iloc[:, :skip_cols].copy()
+        next_step = Series(np.array([""] * len(data.index)))
+        si_deff = data.iloc[:, :skip_cols].copy()
+        prop = data.iloc[:, :skip_cols].copy()
+        free = data.iloc[:, :skip_cols].copy()
+        partitions = data.iloc[:, :skip_cols].copy()
+        suppl = data.iloc[:, :skip_cols].copy()
+
+        # Initializing columns for SiO2 deficiency calculation
+        si_deff['Y'] = 0
+        si_deff['D'] = 0
+
+        # 2 / CO2 handling options
+        if self.verbose > 1: print("Step 2 - CO2 handling options")
+
+        if minor_included:
+            if co2_cancrinite + co2_calcite != 1 and not (co2_cancrinite == co2_calcite == 0):
+                raise Exception('Check the proportions of co2_cancrinite and co2_calcite.')
+            if self.verbose > 1:
+                print('CO2 proportion in calcite : ', co2_calcite, ' --- CO2 proportion in cancrinite:', co2_cancrinite)
+
+        # 3 / Conversion of trace elements units
+        if self.verbose > 1: print("Step 3 - Conversion of trace elements units")
+
+        if minor_included:
+            data_keys = data.keys()
+            for oxel in data_keys:
+                if oxel in list_minor_el and oxel not in ['F', 'Cl', 'S']:
+                    oxide = oxide_corresp[oxel]
+                    data[oxide] = data[oxel] / (ratio_el_to_ox[oxide] * 10000)
+                    data = data.drop(columns=oxel)
+                elif oxel in ['F', 'Cl', 'S']:
+                    data[oxel] = data[oxel] / 10000
+
+            # If some minor elements does not exist in data, add the corresponding oxide with zero values for
+            # treatment. Also add other required elements that are missing.
+
+            data_keys = data.keys()
+            for ox in list_ox_from_minor_el:
+                if ox not in data_keys:
+                    data[ox] = 0
+            if 'F' not in data_keys:
+                data['F'] = 0
+
+        # 4-5 / Adjust components to 100%
+        if self.verbose > 1: print("Step 4, 5 - Adjust components to 100%")
+
+        data_keys = data.keys()
+        total_calculated = data.iloc[:, skip_cols:].sum(axis=1) - data['Total']
+        if self.verbose > 1:
+            print(">>> Total calculated")
+            print(total_calculated.to_string())
+
+        if normalize_entry:
+            oxel_to_sum = list_major_ox if minor_included else list_major_ox
+            for oxel in oxel_to_sum:
+                if oxel in data_keys:
+                    data[oxel] = 100 * data[oxel] / total_calculated
+            total_normalized = data.iloc[:, skip_cols:].sum(axis=1) - data['Total']
+
+            if self.verbose > 1:
+                print(">>> Total normalized")
+                print(total_normalized.to_string())
+        else:
+            total_normalized = total_calculated
+
+        # Calculation of some petrogenetical parameters
+        suppl['pp_FeOt/MgO'] = ((2 * _mol_w['FeO'] / _mol_w['Fe2O3']) * data['Fe2O3'] + data['FeO']) / data['MgO']
+        suppl['pp_ratio_K2O_Na2O'] = data['K2O'] / data['Na2O']
+        suppl['pp_SI'] = 100 * data['MgO'] / (data['MgO'] + data['FeO'] + data['Fe2O3'] + data['Na2O'] + data['K2O'])
+        suppl['pp_Mg#'] = 100 * (data['MgO'] / _mol_w['MgO']) / (
+                data['MgO'] / _mol_w['MgO'] + data['FeO'] / _mol_w['FeO'])
+
+        for i in data.index:
+            if data.loc[i, 'SiO2'] > 5 and 1 < suppl.loc[i, 'pp_ratio_K2O_Na2O'] < 2.5:
+                suppl.loc[i, 'pp_AR'] = (data.loc[i, 'Al2O3'] + data.loc[i, 'CaO'] + 2 * data.loc[i, 'Na2O']) / (
+                        data.loc[i, 'Al2O3'] + data.loc[i, 'CaO'] - 2 * data.loc[i, 'Na2O'])
+            else:
+                suppl.loc[i, 'pp_AR'] = (data.loc[i, 'Al2O3'] + data.loc[i, 'CaO'] + data.loc[i, 'Na2O'] + data.loc[
+                    i, 'K2O']) / (data.loc[i, 'Al2O3'] + data.loc[i, 'CaO'] - data.loc[i, 'Na2O'] - data.loc[i, 'K2O'])
+
+        # 6 / Mole computations
+        if self.verbose > 1: print("Step 6 - Mole computations")
+
+        data_n = data.copy().drop(['Total'], axis=1).fillna(0)
+
+        if self.verbose > 1:
+            print(">>> Initial compositions (wt %)")
+            print(data_n.to_string())
+
+        for oxel in data_n.keys():
+            if oxel in list_all_ox_el:
+                data_n[oxel] = data_n[oxel] / _mol_w[oxel]
+
+        if self.verbose:
+            print(">>> Initial molar concentrations (mol)")
+            print(data_n.to_string())
+
+        # 7-8 / Minor oxides combination and oxide molecular weight computations
+        if self.verbose > 1: print(
+            "Step 7, 8 - Minor oxides combination and oxide molecular weight computations")
+
+        if minor_included:
+            FeO_corr = data_n['FeO'] + data_n['MnO'] + data_n['NiO'] + data_n['CoO']
+            CaO_corr = data_n['CaO'] + data_n['BaO'] + data_n['SrO']
+            K2O_corr = data_n['K2O'] + data_n['Rb2O'] + data_n['Cs2O']
+            Na2O_corr = data_n['Na2O'] + data_n['Li2O']
+            Cr2O3_corr = data_n['Cr2O3'] + data_n['V2O3']
+            xFeO = data_n['FeO'] / FeO_corr
+            xMnO = data_n['MnO'] / FeO_corr
+            xNiO = data_n['NiO'] / FeO_corr
+            xCoO = data_n['CoO'] / FeO_corr
+            xCaO = data_n['CaO'] / CaO_corr
+            xBaO = data_n['BaO'] / CaO_corr
+            xSrO = data_n['SrO'] / CaO_corr
+            xK2O = data_n['K2O'] / K2O_corr
+            xRb2O = data_n['Rb2O'] / K2O_corr
+            xCs2O = data_n['Cs2O'] / K2O_corr
+            xNa2O = data_n['Na2O'] / Na2O_corr
+            xLi2O = data_n['Li2O'] / Na2O_corr
+            xCr2O3 = data_n['Cr2O3'] / Cr2O3_corr
+            xV2O3 = data_n['V2O3'] / Cr2O3_corr
+
+            data_n['FeO'] = data_n['FeO'] + data_n['MnO'] + data_n['NiO'] + data_n['CoO']
+            data_n['CaO'] = data_n['CaO'] + data_n['BaO'] + data_n['SrO']
+            data_n['K2O'] = data_n['K2O'] + data_n['Rb2O'] + data_n['Cs2O']
+            data_n['Na2O'] = data_n['Na2O'] + data_n['Li2O']
+            data_n['Cr2O3'] = data_n['Cr2O3'] + data_n['V2O3']
+
+            # Other oxides to zero
+            for ox in ['MnO', 'NiO', 'CoO', 'BaO', 'SrO', 'Rb2O', 'Cs2O', 'Li2O', 'V2O3']:
+                data_n[ox] = 0
+
+        else:
+            FeO_corr = data_n['MnO'] + data_n['FeO']
+            xFeO = data_n['FeO'] / FeO_corr
+            xMnO = data_n['MnO'] / FeO_corr
+            data_n['FeO'] = data_n['FeO'] + data_n['MnO']
+            data_n['MnO'] = 0
+
+        if self.verbose:
+            print(">>> Corrected molar concentrations (mol)")
+            print(data_n.to_string())
+
+        # Molecular weight correction
+        corr_mol_w = _mol_w.copy()  # Each value will be a number of a list of number (for FeO for example)
+        if minor_included:
+            corr_mol_w['FeO'] = xFeO * _mol_w['FeO'] + xMnO * _mol_w['MnO'] + \
+                                xNiO * _mol_w['NiO'] + xCoO * _mol_w['CoO']
+            corr_mol_w['CaO'] = xCaO * _mol_w['CaO'] + xBaO * _mol_w['BaO'] + xSrO * _mol_w['SrO']
+            corr_mol_w['K2O'] = xK2O * _mol_w['K2O'] + xRb2O * _mol_w['Cs2O'] + xCs2O * _mol_w['Cs2O']
+            corr_mol_w['Na2O'] = xNa2O * _mol_w['Na2O'] + xLi2O * _mol_w['Li2O']
+            corr_mol_w['Cr2O3'] = xCr2O3 * _mol_w['Cr2O3'] + xV2O3 * _mol_w['V2O3']
+
+            corr_mol_w['FeO'] = corr_mol_w['FeO'].fillna(0)
+            corr_mol_w['CaO'] = corr_mol_w['CaO'].fillna(0)
+            corr_mol_w['K2O'] = corr_mol_w['K2O'].fillna(0)
+            corr_mol_w['Na2O'] = corr_mol_w['Na2O'].fillna(0)
+            corr_mol_w['Cr2O3'] = corr_mol_w['Cr2O3'].fillna(0)
+
+        else:
+            corr_mol_w['FeO'] = xFeO * _mol_w['FeO'] + xMnO * _mol_w['MnO']
+            corr_mol_w['FeO'] = corr_mol_w['FeO'].fillna(0)
+
+        # 9-10 / Correction of normative mineral molecular weights
+        if self.verbose > 1: print("Step 9, 10 - Correction of normative mineral molecular weights")
+
+        mol_w_min = {
+            'Ab': corr_mol_w['Na2O'] + corr_mol_w['Al2O3'] + 6 * corr_mol_w['SiO2'],
+            'Ne': corr_mol_w['Na2O'] + corr_mol_w['Al2O3'] + 2 * corr_mol_w['SiO2'],
+            'Th': corr_mol_w['Na2O'] + corr_mol_w['SO3'],
+            'Nc': corr_mol_w['Na2O'] + corr_mol_w['CO2'],
+            'An': corr_mol_w['CaO'] + corr_mol_w['Al2O3'] + 2 * corr_mol_w['SiO2'],
+            'Di-mg': corr_mol_w['CaO'] + corr_mol_w['MgO'] + 2 * corr_mol_w['SiO2'],
+            'Di-fe': corr_mol_w['CaO'] + corr_mol_w['FeO'] + 2 * corr_mol_w['SiO2'],
+            'Ac': corr_mol_w['Na2O'] + corr_mol_w['Fe2O3'] + 4 * corr_mol_w['SiO2'],
+            'Ns': corr_mol_w['Na2O'] + corr_mol_w['SiO2'],
+            'Or': corr_mol_w['K2O'] + corr_mol_w['Al2O3'] + 6 * corr_mol_w['SiO2'],
+            'Lc': corr_mol_w['K2O'] + corr_mol_w['Al2O3'] + 4 * corr_mol_w['SiO2'],
+            'Kp': corr_mol_w['K2O'] + corr_mol_w['Al2O3'] + 2 * corr_mol_w['SiO2'],
+            'Ks': corr_mol_w['K2O'] + corr_mol_w['SiO2'],
+            'Wo': corr_mol_w['CaO'] + corr_mol_w['SiO2'],
+            'Cs': 2 * corr_mol_w['CaO'] + corr_mol_w['SiO2'],
+            'Mt': corr_mol_w['FeO'] + corr_mol_w['Fe2O3'],
+            'Il': corr_mol_w['FeO'] + corr_mol_w['TiO2'],
+            'Tn': corr_mol_w['CaO'] + corr_mol_w['TiO2'] + corr_mol_w['SiO2'],
+            'Pf': corr_mol_w['CaO'] + corr_mol_w['TiO2'],
+            'Cc': corr_mol_w['CaO'] + corr_mol_w['CO2'],
+            'Hy-fe': corr_mol_w['FeO'] + corr_mol_w['SiO2'],
+            'Ol-fe': 2 * corr_mol_w['FeO'] + corr_mol_w['SiO2'],
+            'Q': corr_mol_w['SiO2'],
+            'C': corr_mol_w['Al2O3'],
+            'Z': corr_mol_w['SiO2'] + corr_mol_w['ZrO2'],
+            'Hy-mg': corr_mol_w['MgO'] + corr_mol_w['SiO2'],
+            'Ol-mg': 2 * corr_mol_w['MgO'] + corr_mol_w['SiO2'],
+            'Hm': corr_mol_w['Fe2O3'],
+            'Ru': corr_mol_w['TiO2'],
+            'Ap-F': 3 * corr_mol_w['CaO'] + corr_mol_w['P2O5'] + 1 / 3 * (
+                    (corr_mol_w['CaO'] - _mol_w['O']) + 2 * corr_mol_w['F']),
+            'Ap-O': 3 * corr_mol_w['CaO'] + corr_mol_w['P2O5'] + corr_mol_w['CaO'] * 1 / 3,
+            'Cm': corr_mol_w['FeO'] + corr_mol_w['Cr2O3'],
+            'Hl': (corr_mol_w['Na2O'] - _mol_w['O']) / 2 + _mol_w['Cl'],
+            'Fr': (corr_mol_w['CaO'] - _mol_w['O']) + 2 * _mol_w['F'],
+            'Pr': (corr_mol_w['FeO'] - _mol_w['O']) + 2 * _mol_w['S']
+        }
+
+        # 11 / Normative zircon
+        if self.verbose > 1: print("Step 11 - Normative zircon")
+
+        if minor_included:
+            for i in data_n.index:
+                if data_n.loc[i, 'SiO2'] > data_n.loc[i, 'ZrO2']:
+                    n_phase.loc[i, 'Z'] = data_n.loc[i, 'ZrO2']
+                    si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + n_phase.loc[i, 'Z']
+                    data_n.loc[i, 'ZrO2'] = 0
+                else:
+                    print('WARNING : No further SiO2 after zircon attribution for composition ' + str(
+                        i) + '. Check data.')
+
+        # 12 / Normative apatite
+        if self.verbose > 1: print("Step 12 - Normative apatite")
+
+        temp_ap = np.zeros(shape=[len(data_n.index), 1])
+        for i in data_n.index:
+            if data_n.loc[i, 'CaO'] >= (3 + 1 / 3) * data_n.loc[i, 'P2O5']:
+                temp_ap[i] = data_n.loc[i, 'P2O5']
+                data_n.loc[i, 'CaO'] = data_n.loc[i, 'CaO'] - (3 + 1 / 3) * data_n.loc[i, 'P2O5']
+            else:
+                temp_ap[i] = data_n.loc[i, 'CaO'] / (3 + 1 / 3)
+                data_n.loc[i, 'CaO'] = 0
+                free.loc[i, 'P2O5'] = data_n.loc[i, 'P2O5'] - temp_ap[i]
+            free.loc[i, 'P2O5'] = data_n.loc[i, 'P2O5'] - temp_ap[i]
+            data_n.loc[i, 'P2O5'] = 0
+
+        n_phase['Ap-F'] = 0
+        n_phase['Ap-O'] = 0
+        if minor_included:
+            for i in data_n.index:
+                if data_n.loc[i, 'F'] >= 2 / 3 * temp_ap[i]:
+                    n_phase.loc[i, 'Ap-F'] = temp_ap[i]
+                    free.loc[i, 'O_12b'] = 1 / 3 * n_phase.loc[i, 'Ap-F']
+                    data_n.loc[i, 'F'] = data_n.loc[i, 'F'] - 2 / 3 * n_phase.loc[i, 'Ap-F']
+                else:
+                    n_phase.loc[i, 'Ap-F'] = 1.5 * data_n.loc[i, 'F']
+                    n_phase.loc[i, 'Ap-O'] = temp_ap[i] - 1.5 * data_n.loc[i, 'F']
+                    free.loc[i, 'O_12c'] = data_n.loc[i, 'F'] / 2
+                    data_n.loc[i, 'F'] = 0
+            n_phase['Ap'] = n_phase['Ap-F'] + n_phase['Ap-O']
+        else:
+            n_phase['Ap-O'] = temp_ap
+        n_phase['Ap'] = n_phase['Ap-O'] + n_phase['Ap-F']
+
+        # 13 / Normative fluorite
+        if self.verbose > 1: print("Step 13 - Normative fluorite")
+
+        if minor_included:
+            for i in data_n.index:
+                if data_n.loc[i, 'CaO'] >= data_n.loc[i, 'F'] / 2:
+                    n_phase.loc[i, 'Fr'] = data_n.loc[i, 'F'] / 2
+                    data_n.loc[i, 'CaO'] = data_n.loc[i, 'CaO'] - data_n.loc[i, 'F'] / 2
+                    data_n.loc[i, 'F'] = 0
+                else:
+                    n_phase.loc[i, 'Fr'] = data_n.loc[i, 'CaO']
+                    data_n.loc[i, 'CaO'] = 0
+                    free.loc[i, 'F'] = data_n.loc[i, 'F'] - 2 * n_phase.loc[i, 'Fr']  # Unused F
+                    data_n.loc[i, 'F'] = 0
+                free.loc[i, 'O_13'] = n_phase.loc[i, 'Fr']
+
+        # 14 / Normative halite
+        if self.verbose > 1: print("Step 14 - Normative halite")
+
+        if minor_included and 'Cl' in data_n.keys():
+            for i in data_n.index:
+                if data_n.loc[i, 'Na2O'] >= 2 * data_n.loc[i, 'Cl']:
+                    n_phase.loc[i, 'Hl'] = data_n.loc[i, 'Cl']
+                    data_n.loc[i, 'Na2O'] = data_n.loc[i, 'Na2O'] - n_phase.loc[i, 'Hl'] / 2
+                    free.loc[i, 'O_14'] = n_phase.loc[i, 'Hl'] / 2
+                    data_n.loc[i, 'Cl'] = 0
+                else:
+                    n_phase.loc[i, 'Hl'] = data_n.loc[i, 'Na2O'] / 2
+                    free.loc[i, 'Cl'] = data_n.loc[i, 'Cl'] - n_phase.loc[i, 'Hl']
+                    data_n.loc[i, 'Na2O'] = 0
+
+        # 15 / Normative thenardite
+        if self.verbose > 1: print("Step 15 - Normative thenardite")
+
+        if minor_included and 'SO3' in data_n.keys():
+            for i in data_n.index:
+                if data_n.loc[i, 'Na2O'] >= 2 * data_n.loc[i, 'SO3']:
+                    n_phase.loc[i, 'Th'] = data_n.loc[i, 'SO3']
+                    data_n.loc[i, 'Na2O'] = data_n.loc[i, 'Na2O'] - n_phase.loc[i, 'Th']
+                    data_n.loc[i, 'SO3'] = 0
+                else:
+                    n_phase.loc[i, 'Th'] = data_n.loc[i, 'Na2O']
+                    free.loc[i, 'SO3'] = data_n.loc[i, 'SO3'] - n_phase.loc[i, 'Th']
+                    data_n.loc[i, 'Na2O'] = 0
+
+        # 16 / Normative pyrite
+        if self.verbose > 1: print("Step 16 - Normative pyrite")
+
+        if minor_included and 'S' in data_n.keys():
+            for i in data_n.index:
+                if data_n.loc[i, 'FeO'] >= 2 * data_n.loc[i, 'S']:
+                    n_phase.loc[i, 'Pr'] = data_n.loc[i, 'S'] / 2
+                    data_n.loc[i, 'FeO'] = data_n.loc[i, 'FeO'] - n_phase.loc[i, 'Pr']
+                    data_n.loc[i, 'S'] = 0
+                else:
+                    n_phase.loc[i, 'Pr'] = data_n.loc[i, 'FeO']
+                    free.loc[i, 'S'] = data_n.loc[i, 'S'] - 2 * n_phase.loc[i, 'Pr']
+                    data_n.loc[i, 'FeO'] = 0
+                free.loc[i, 'O_16'] = n_phase.loc[i, 'Pr']
+
+        # 17 / Normative sodium carbonate or calcite
+        if self.verbose > 1: print("Step 17 - Normative sodium carbonate or calcite")
+
+        if minor_included:
+            if 'CO2' in data_keys:
+                free['CO2'] = 0
+                if co2_cancrinite == co2_calcite == 0:
+                    for i in data_n.index:
+                        free.loc[i, 'CO2'] = data_n.loc[i, 'CO2']
+                        data_n.loc[i, 'CO2'] = 0
+                else:
+                    # Cancrinite
+                    for i in data_n.index:
+                        if data_n.loc[i, 'CO2'] > 0:
+                            if data_n.loc[i, 'Na2O'] >= data_n.loc[i, 'CO2'] * co2_cancrinite:
+                                n_phase.loc[i, 'Nc'] = data_n.loc[i, 'CO2'] * co2_cancrinite
+                                data_n.loc[i, 'Na2O'] = data_n.loc[i, 'Na2O'] - n_phase.loc[i, 'Nc']
+                            else:
+                                n_phase.loc[i, 'Nc'] = data_n.loc[i, 'Na2O']
+                                data_n.loc[i, 'Na2O'] = 0
+                                free.loc[i, 'CO2'] = free.loc[i, 'CO2'] + data_n.loc[i, 'CO2'] * co2_cancrinite \
+                                                     - n_phase.loc[i, 'Nc']  # Free CO2
+                            data_n.loc[i, 'CO2'] = data_n.loc[i, 'CO2'] * (1 - co2_cancrinite)
+
+                    # Calcite, with all the rest of CO2
+                    for i in data_n.index:
+                        if data_n.loc[i, 'CO2'] > 0:
+                            if data_n.loc[i, 'CaO'] > data_n.loc[i, 'CO2']:
+                                n_phase.loc[i, 'Cc'] = data_n.loc[i, 'CO2']
+                                data_n.loc[i, 'CaO'] = data_n.loc[i, 'CaO'] - n_phase.loc[i, 'Cc']
+                            else:
+                                n_phase.loc[i, 'Cc'] = data_n.loc[i, 'CaO']
+                                data_n.loc[i, 'CaO'] = 0
+                                free.loc[i, 'CO2'] = free.loc[i, 'CO2'] + data_n.loc[i, 'CO2'] - n_phase.loc[
+                                    i, 'Cc']  # CO2 to free
+                            data_n.loc[i, 'CO2'] = 0
+
+        # 18 / Normative chromite
+        if self.verbose > 1: print("Step 18 - Normative chromite")
+
+        if minor_included:
+            for i in data_n.index:
+                if 0 < data_n.loc[i, 'Cr2O3'] <= data_n.loc[i, 'FeO']:
+                    n_phase.loc[i, 'Cm'] = data_n.loc[i, 'Cr2O3']
+                    data_n.loc[i, 'FeO'] = data_n.loc[i, 'FeO'] - n_phase.loc[i, 'Cm']
+                    data_n.loc[i, 'Cr2O3'] = 0
+                elif data_n.loc[i, 'FeO'] < data_n.loc[i, 'Cr2O3']:
+                    n_phase.loc[i, 'Cm'] = data_n.loc[i, 'FeO']
+                    free.loc[i, 'Cr2O3'] = data_n.loc[i, 'Cr2O3'] - n_phase.loc[i, 'Cm']
+                    data_n.loc[i, 'FeO'] = 0
+                    data_n.loc[i, 'Cr2O3'] = 0
+
+        # 19 / Normative ilmenite
+        if self.verbose > 1: print("Step 19 - Normative ilmenite")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'FeO'] >= data_n.loc[i, 'TiO2']:
+                n_phase.loc[i, 'Il'] = data_n.loc[i, 'TiO2']
+                data_n.loc[i, 'FeO'] = data_n.loc[i, 'FeO'] - n_phase.loc[i, 'Il']
+                data_n.loc[i, 'TiO2'] = 0
+            else:
+                n_phase.loc[i, 'Il'] = data_n.loc[i, 'FeO']
+                data_n.loc[i, 'TiO2'] = data_n.loc[i, 'TiO2'] - n_phase.loc[i, 'Il']
+                data_n.loc[i, 'FeO'] = 0
+
+        # 20 / Normative orthoclase or potassium metasilicate
+        if self.verbose > 1: print("Step 20 - Normative orthoclase or potassium metasilicate")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'Al2O3'] >= data_n.loc[i, 'K2O']:
+                n_phase.loc[i, 'Orp'] = data_n.loc[i, 'K2O']
+                data_n.loc[i, 'Al2O3'] = data_n.loc[i, 'Al2O3'] - n_phase.loc[i, 'Orp']
+                data_n.loc[i, 'K2O'] = 0
+                si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 6 * n_phase.loc[i, 'Orp']
+            else:
+                n_phase.loc[i, 'Orp'] = data_n.loc[i, 'Al2O3']
+                n_phase.loc[i, 'Ks'] = data_n.loc[i, 'K2O'] - n_phase.loc[i, 'Orp']  # Rest of K2O to Ks
+                data_n.loc[i, 'Al2O3'] = 0
+                data_n.loc[i, 'K2O'] = 0
+                si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 6 * n_phase.loc[i, 'Orp'] + n_phase.loc[i, 'Ks']
+
+        # 21 / Normative albite
+        if self.verbose > 1: print("Step 21 - Normative albite")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'Al2O3'] >= data_n.loc[i, 'Na2O']:
+                n_phase.loc[i, 'Abp'] = data_n.loc[i, 'Na2O']
+                data_n.loc[i, 'Al2O3'] = data_n.loc[i, 'Al2O3'] - n_phase.loc[i, 'Abp']
+                data_n.loc[i, 'Na2O'] = 0
+            else:
+                n_phase.loc[i, 'Abp'] = data_n.loc[i, 'Al2O3']
+                data_n.loc[i, 'Na2O'] = data_n.loc[i, 'Na2O'] - n_phase.loc[i, 'Abp']  # Rest of Na2O to Ab
+                data_n.loc[i, 'Al2O3'] = 0
+            si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 6 * n_phase.loc[i, 'Abp']
+
+        # 22 / Normative acmite or sodium metasilicate
+        if self.verbose > 1: print("Step 22 - Normative acmite or sodium metasilicate")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'Na2O'] >= data_n.loc[i, 'Fe2O3']:
+                n_phase.loc[i, 'Ac'] = data_n.loc[i, 'Fe2O3']
+                n_phase.loc[i, 'Ns'] = data_n.loc[i, 'Na2O'] - n_phase.loc[i, 'Ac']  # Rest of Na2O to Ns
+                data_n.loc[i, 'Na2O'] = 0
+                data_n.loc[i, 'Fe2O3'] = 0
+                si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 4 * n_phase.loc[i, 'Ac'] + n_phase.loc[i, 'Ns']
+            else:
+                n_phase.loc[i, 'Ac'] = data_n.loc[i, 'Na2O']
+                data_n.loc[i, 'Fe2O3'] = data_n.loc[i, 'Fe2O3'] - n_phase.loc[i, 'Ac']
+                data_n.loc[i, 'Na2O'] = 0
+                si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 4 * n_phase.loc[i, 'Ac']
+
+        # 23 / Normative anorthite or corundum
+        if self.verbose > 1: print("Step 23 - Normative anorthite or corundum")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'Al2O3'] >= data_n.loc[i, 'CaO']:
+                n_phase.loc[i, 'An'] = data_n.loc[i, 'CaO']
+                n_phase.loc[i, 'C'] = data_n.loc[i, 'Al2O3'] - data_n.loc[i, 'CaO']  # Rest of Al2O3 in C
+                data_n.loc[i, 'CaO'] = 0
+                data_n.loc[i, 'Al2O3'] = 0
+            else:
+                n_phase.loc[i, 'An'] = data_n.loc[i, 'Al2O3']
+                data_n.loc[i, 'CaO'] = data_n.loc[i, 'CaO'] - n_phase.loc[i, 'An']
+                data_n.loc[i, 'Al2O3'] = 0
+            si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 2 * n_phase.loc[i, 'An']
+
+        # 24 / Normative sphene / rutile
+        if self.verbose > 1: print("Step 24 - Normative sphene / rutile")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'CaO'] >= data_n.loc[i, 'TiO2']:
+                n_phase.loc[i, 'Tnp'] = data_n.loc[i, 'TiO2']
+                data_n.loc[i, 'CaO'] = data_n.loc[i, 'CaO'] - n_phase.loc[i, 'Tnp']
+                data_n.loc[i, 'TiO2'] = 0
+            else:
+                n_phase.loc[i, 'Tnp'] = data_n.loc[i, 'CaO']
+                n_phase.loc[i, 'Ru'] = data_n.loc[i, 'TiO2'] - n_phase.loc[i, 'Tnp']  # Rest of TiO2 in Ru
+                data_n.loc[i, 'CaO'] = 0
+                data_n.loc[i, 'TiO2'] = 0
+            si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + n_phase.loc[i, 'Tnp']
+
+        # 25 / Normative magnetite or hematite
+        if self.verbose > 1: print("Step 25 - Normative magnetite or hematite")
+
+        for i in data_n.index:
+            if data_n.loc[i, 'Fe2O3'] >= data_n.loc[i, 'FeO']:
+                n_phase.loc[i, 'Mt'] = data_n.loc[i, 'FeO']
+                n_phase.loc[i, 'Hm'] = data_n.loc[i, 'Fe2O3'] - n_phase.loc[i, 'Mt']  # Rest of Fe2O3 in Mt
+                data_n.loc[i, 'FeO'] = 0
+                data_n.loc[i, 'Fe2O3'] = 0
+            else:
+                n_phase.loc[i, 'Mt'] = data_n.loc[i, 'Fe2O3']
+                data_n.loc[i, 'FeO'] = data_n.loc[i, 'FeO'] - n_phase.loc[i, 'Mt']
+                data_n.loc[i, 'Fe2O3'] = 0
+
+        # 26 / Subdivision of Mg and Fe in some minerals
+        if self.verbose > 1: print("Step 26 - Repartition of Mg and Fe in minerals")
+
+        data_n['FeMgO'] = data_n['MgO'] + data_n['FeO']
+        prop['xMg'] = data_n['MgO'] / (data_n['MgO'] + data_n['FeO'])
+        prop['xFe'] = data_n['FeO'] / (data_n['MgO'] + data_n['FeO'])
+        data_n['MgO'] = 0
+        data_n['FeO'] = 0
+
+        # 27 / Provisional normative diopside, wollastonite or hypersthene
+        if self.verbose > 1: print("Step 27 - normative diopside, wollastonite or hypersthene")
+
+        n_phase['Wop'] = 0
+        for i in data_n.index:
+            if data_n.loc[i, 'CaO'] >= data_n.loc[i, 'FeMgO']:
+                n_phase.loc[i, 'Dip'] = data_n.loc[i, 'FeMgO']
+                n_phase.loc[i, 'Wop'] = data_n.loc[i, 'CaO'] - n_phase.loc[i, 'Dip']
+                data_n.loc[i, 'CaO'] = 0
+                data_n.loc[i, 'FeMgO'] = 0
+                si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 2 * n_phase.loc[i, 'Dip'] + n_phase.loc[i, 'Wop']
+                n_phase.loc[i, 'Hyp'] = 0
+            else:
+                n_phase.loc[i, 'Dip'] = data_n.loc[i, 'CaO']
+                n_phase.loc[i, 'Hyp'] = data_n.loc[i, 'FeMgO'] - n_phase.loc[i, 'Dip']
+                data_n.loc[i, 'CaO'] = 0
+                data_n.loc[i, 'FeMgO'] = 0
+                si_deff.loc[i, 'Y'] = si_deff.loc[i, 'Y'] + 2 * n_phase.loc[i, 'Dip'] + n_phase.loc[i, 'Hyp']
+
+        # 28 / Normative quartz and Si deficiency
+        if self.verbose > 1: print("Step 28 - Normative quartz and Si deficiency")
+        n_phase['Q'] = 0
+        for i in data_n.index:
+            if data_n.loc[i, 'SiO2'] >= si_deff.loc[i, 'Y']:
+                n_phase.loc[i, 'Q'] = data_n.loc[i, 'SiO2'] - si_deff.loc[i, 'Y']
+                data_n.loc[i, 'SiO2'] = 0
+                next_step[i] = '36b'  # Instead of 36a which is obligatory in this algorithm
+                if self.verbose: print("Si saturated for the composition", i)
+            else:
+                si_deff.loc[i, 'D'] = si_deff.loc[i, 'Y'] - data_n.loc[i, 'SiO2']
+                data_n.loc[i, 'SiO2'] = 0
+                next_step[i] = '29'
+
+        # 29 / Normative olivine or hypersthene
+        if self.verbose > 1: print("Step 29 - Normative olivine or hypersthene")
+        for i in data_n.index:
+            if next_step[i] == '29':
+                if si_deff.loc[i, 'D'] < n_phase.loc[i, 'Hyp'] / 2:
+                    n_phase.loc[i, 'Ol'] = si_deff.loc[i, 'D']
+                    n_phase.loc[i, 'Hy'] = n_phase.loc[i, 'Hyp'] - 2 * si_deff.loc[i, 'D']
+                    next_step[i] = '36b'
+                else:
+                    n_phase.loc[i, 'Ol'] = n_phase.loc[i, 'Hyp'] / 2
+                    n_phase.loc[i, 'Hy'] = 0
+                    si_deff.loc[i, 'D1'] = si_deff.loc[i, 'D'] - n_phase.loc[i, 'Hyp'] / 2
+                    next_step[i] = '30'
+            else:
+                n_phase.loc[i, 'Hy'] = n_phase.loc[i, 'Hyp']
+                n_phase.loc[i, 'Ol'] = 0
+
+        # 30 / Normative sphene or perovskite
+        if self.verbose > 1: print("Step 30 - Normative sphene or perovskite")
+
+        for i in data_n.index:
+            if next_step[i] == '30':
+                if si_deff.loc[i, 'D1'] < n_phase.loc[i, 'Tnp']:
+                    n_phase.loc[i, 'Tn'] = n_phase.loc[i, 'Tnp'] - si_deff.loc[i, 'D1']
+                    n_phase.loc[i, 'Pf'] = si_deff.loc[i, 'D1']
+                    next_step[i] = '36c'
+                else:
+                    n_phase.loc[i, 'Pf'] = n_phase.loc[i, 'Tnp']
+                    n_phase.loc[i, 'Tn'] = 0
+                    si_deff.loc[i, 'D2'] = si_deff.loc[i, 'D1'] - n_phase.loc[i, 'Tnp']
+                    next_step[i] = '31'
+
+        # 31 / Normative nepheline or albite
+        if self.verbose > 1: print("Step 31 - Normative nepheline or albite")
+
+        for i in data_n.index:
+            if next_step[i] == '31':
+                if si_deff.loc[i, 'D2'] < 4 * n_phase.loc[i, 'Abp']:
+                    n_phase.loc[i, 'Ab'] = n_phase.loc[i, 'Abp'] - si_deff.loc[i, 'D2'] / 4
+                    n_phase.loc[i, 'Ne'] = si_deff.loc[i, 'D2'] / 4
+                    next_step[i] = '36d'
+                else:
+                    n_phase.loc[i, 'Ne'] = n_phase.loc[i, 'Abp']
+                    n_phase.loc[i, 'Ab'] = 0
+                    si_deff.loc[i, 'D3'] = si_deff.loc[i, 'D2'] - 4 * n_phase.loc[i, 'Abp']
+                    next_step[i] = '32'
+
+        # 32 / Normative leucite or orthoclase
+        if self.verbose > 1: print("Step 32 - Normative leucite or orthoclase")
+
+        n_phase['Lcp'] = 0
+        for i in data_n.index:
+            if next_step[i] == '32':
+                if si_deff.loc[i, 'D3'] < 2 * n_phase.loc[i, 'Orp']:
+                    n_phase.loc[i, 'Or'] = n_phase.loc[i, 'Orp'] - si_deff.loc[i, 'D3'] / 2
+                    n_phase.loc[i, 'Lc'] = si_deff.loc[i, 'D3'] / 2
+                    next_step[i] = '36e'
+                else:
+                    n_phase.loc[i, 'Lcp'] = n_phase.loc[i, 'Orp']
+                    n_phase.loc[i, 'Or'] = 0
+                    si_deff.loc[i, 'D4'] = si_deff.loc[i, 'D3'] - 2 * n_phase.loc[i, 'Orp']
+                    next_step[i] = '33'
+
+        # 33 / Normative dicalcium silicate or wollastonite
+        if self.verbose > 1: print("Step 33 - Normative dicalcium silicate or wollastonite")
+
+        for i in data_n.index:
+            if next_step[i] == '33':
+                if si_deff.loc[i, 'D4'] < 2 * n_phase.loc[i, 'Wop'] / 2:
+                    n_phase.loc[i, 'Wo'] = n_phase.loc[i, 'Wop'] - 2 * si_deff.loc[i, 'D4']
+                    n_phase.loc[i, 'Cs'] = si_deff.loc[i, 'D4']
+                    next_step[i] = '36e'  # Instead of 36f
+                else:
+                    n_phase.loc[i, 'Cs'] = n_phase.loc[i, 'Wop'] / 2
+                    n_phase.loc[i, 'Wo'] = 0
+                    si_deff.loc[i, 'D5'] = si_deff.loc[i, 'D4'] - n_phase.loc[i, 'Wop'] / 2
+                    next_step[i] = '34'
+
+        # 34 / Normative dicalcium silicate or olivine
+        if self.verbose > 1: print("Step 34 - Normative diopside or olivine adjustment")
+
+        for i in data_n.index:
+            if next_step[i] == '34':
+                if si_deff.loc[i, 'D5'] < n_phase.loc[i, 'Dip']:
+                    n_phase.loc[i, 'Cs'] = n_phase.loc[i, 'Cs'] + si_deff.loc[i, 'D5'] / 2
+                    n_phase.loc[i, 'Ol'] = n_phase.loc[i, 'Ol'] + si_deff.loc[i, 'D5'] / 2
+                    n_phase.loc[i, 'Di'] = n_phase.loc[i, 'Dip'] - si_deff.loc[i, 'D5']
+                    next_step[i] = '36g'
+                else:
+                    n_phase.loc[i, 'Cs'] = n_phase.loc[i, 'Cs'] + n_phase.loc[i, 'Dip'] / 2
+                    n_phase.loc[i, 'Ol'] = n_phase.loc[i, 'Ol'] + n_phase.loc[i, 'Dip'] / 2
+                    n_phase.loc[i, 'Di'] = n_phase.loc[i, 'Di-fe'] = n_phase.loc[i, 'Di-mg'] = 0
+                    si_deff.loc[i, 'D6'] = si_deff.loc[i, 'D5'] - n_phase.loc[i, 'Dip']
+                    next_step[i] = '35'
+            else:
+                n_phase['Di'] = n_phase['Dip']
+
+        # 35 / Normative kaliophilite or leucite
+        if self.verbose > 1: print("Step 35 - Normative kaliophilite or leucite")
+
+        for i in data_n.index:
+            if next_step[i] == '35':
+                if n_phase.loc[i, 'Lcp'] >= si_deff.loc[i, 'D6'] / 2:
+                    n_phase.loc[i, 'Kp'] = si_deff.loc[i, 'D6'] / 2
+                    n_phase.loc[i, 'Lc'] = n_phase.loc[i, 'Lcp'] - si_deff.loc[i, 'D6'] / 2
+                    next_step[i] = '36g'
+                else:
+                    n_phase.loc[i, 'Kp'] = n_phase.loc[i, 'Lcp']
+                    n_phase.loc[i, 'Lc'] = 0
+                    suppl.loc[i, 'defSiO2'] = (si_deff.loc[i, 'D6'] - 2 * n_phase.loc[i, 'Kp']) * corr_mol_w['SiO2']
+                    next_step[i] = '36g'
+
+        # 36 / Definite mineral proportions
+
+        for i in data_n.index:
+            if self.verbose > 1: print("Step 36 - Solution", i)
+
+            if next_step[i] == '36b':
+                if self.verbose > 1: print("     36b")
+                n_phase.loc[i, 'Tn'] = n_phase.loc[i, 'Tnp']
+                next_step[i] = '36c'
+
+            if next_step[i] == '36c':
+                if self.verbose > 1: print("     36c")
+                n_phase.loc[i, 'Ab'] = n_phase.loc[i, 'Abp']
+                next_step[i] = '36d'
+
+            if next_step[i] == '36d':
+                if self.verbose > 1: print("     36d")
+                n_phase.loc[i, 'Or'] = n_phase.loc[i, 'Orp']
+                n_phase.loc[i, 'Lc'] = n_phase.loc[i, 'Lcp']
+                next_step[i] = '36e'
+
+            if next_step[i] == '36e':
+                if self.verbose > 1: print("     36e")
+                if 'Wop' in n_phase.keys():
+                    n_phase.loc[i, 'Wo'] = n_phase.loc[i, 'Wop']
+
+        # Steps of 36a and 36f, obligatory for this algorithm
+        if self.verbose > 1: print("     Obilgatory 36a, 36f and Fe-Mg distribution in olivine")
+        n_phase['Hy-fe'] = n_phase['Hy'] * prop['xFe']
+        n_phase['Hy-mg'] = n_phase['Hy'] * prop['xMg']
+        n_phase['Di-fe'] = n_phase['Di'] * prop['xFe']
+        n_phase['Di-mg'] = n_phase['Di'] * prop['xMg']
+        n_phase['Ol-fe'] = n_phase['Ol'] * prop['xFe']
+        n_phase['Ol-mg'] = n_phase['Ol'] * prop['xMg']
+
+        # 37 / Conversion of normative minerals in %, normative sum
+        if self.verbose > 1: print("Step 37 - Conversion of normative minerals in %, normative sum")
+
+        ignored_phases = []
+        for mineral in n_phase.keys():
+            if mineral in final_list_min:
+                partitions[mineral] = n_phase[mineral] * mol_w_min[mineral]
+            elif mineral == 'Ap':
+                suppl['Ap-F'] = n_phase['Ap-F'] * mol_w_min['Ap-F']
+                suppl['Ap-O'] = n_phase['Ap-O'] * mol_w_min['Ap-O']
+                partitions['Ap'] = suppl['Ap-F'] + suppl['Ap-O']
+            elif mineral == 'Ol':
+                suppl['Ol-mg'] = n_phase['Ol-mg'] * mol_w_min['Ol-mg']
+                suppl['Ol-fe'] = n_phase['Ol-fe'] * mol_w_min['Ol-fe']
+                partitions['Ol'] = suppl['Ol-mg'] + suppl['Ol-fe']
+            elif mineral == 'Hy':
+                suppl['Hy-mg'] = n_phase['Hy-mg'] * mol_w_min['Hy-mg']
+                suppl['Hy-fe'] = n_phase['Hy-fe'] * mol_w_min['Hy-fe']
+                partitions['Hy'] = suppl['Hy-mg'] + suppl['Hy-fe']
+            elif mineral == 'Di':
+                suppl['Di-mg'] = n_phase['Di-mg'] * mol_w_min['Di-mg']
+                suppl['Di-fe'] = n_phase['Di-fe'] * mol_w_min['Di-fe']
+                partitions['Di'] = suppl['Di-mg'] + suppl['Di-fe']
+            else:
+                ignored_phases.append(mineral)
+
+        if self.verbose > 1: print("Temporary phases ignored in final composition :", *ignored_phases[skip_cols:])
+
+        # Add free CO2
+        free = free.fillna(0)
+        free_phases = free.keys()
+        if 'CO2' in free_phases:
+            partitions['CO2'] = free['CO2'] * corr_mol_w['CO2']
+
+        # Add free O
+        if minor_included:
+            # Free O
+            free['O_wt%'] = 0
+            if 'O_12b' in free_phases:
+                free['O_wt%'] = free['O_wt%'] + (1 + (0.1 * (mol_w_min['Ap-F'] / 328.8691887) - 1)) * \
+                                corr_mol_w['O'] * free['O_12b']
+            if 'O_12c' in free_phases:
+                free['O_wt%'] = free['O_wt%'] + (
+                        1 + 0.1 * (n_phase['Ap-F'] / n_phase['Ap']) * (mol_w_min['Ap-F'] / 328.8691887) - 1) * \
+                                corr_mol_w['O'] * free['O_12c']
+            if 'O_13' in free_phases:
+                free['O_wt%'] = free['O_wt%'] + (1 + (corr_mol_w['CaO'] / 56.0774 - 1)) * corr_mol_w['O'] * free['O_13']
+            if 'O_14' in free_phases:
+                free['O_wt%'] = free['O_wt%'] + (1 + 0.5 * (corr_mol_w['Na2O'] / 61.97894 - 1)) * corr_mol_w['O'] * \
+                                free['O_14']
+            if 'O_16' in free_phases:
+                free['O_wt%'] = free['O_wt%'] + (1 + (corr_mol_w['FeO'] / 71.8444 - 1)) * corr_mol_w['O'] * free['O_16']
+            partitions['O'] = free['O_wt%']
+
+        # Add free oxides
+        partitions['free_Ox'] = 0
+        free_keys = free.keys()
+        for oxel in ['P2O5', 'F', 'Cl', 'SO3', 'Cr2O3']:
+            if oxel in free_keys:
+                partitions['free_oxides'] = partitions['free_Ox'] + free[oxel] * _mol_w[oxel]
+
+        # Totals
+        partitions['Sum_norm'] = partitions.iloc[:, skip_cols:].sum(axis=1)
+        partitions = partitions.round(to_round)
+
+        # 38 / Correctness of normative sum
+        if self.verbose > 1: print("Step 38 - Correctness of normative sum")
+
+        suppl['diff_sum'] = total_normalized - partitions['Sum_norm']
+
+        # 39 / Other petrogenetical parameters
+        if self.verbose > 1: print("Step 39 - Other petrogenetical parameters")
+
+        suppl['pp_salic'] = n_phase['Q'] + n_phase['Or'] + n_phase['Ab'] + n_phase['An']
+        suppl['pp_femic'] = suppl['Di-mg'] + suppl['Di-fe'] + suppl['Hy-mg'] + suppl['Hy-fe'] + suppl['Ol-mg'] + \
+                            n_phase['Ol-fe'] + n_phase['Mt'] + n_phase['Il'] + n_phase['Hm']
+        suppl['pp_CI'] = n_phase['An'] + 2.1570577 * suppl['Di-mg'] + suppl['Ol-mg'] + 0.7007616 * suppl['Hy-fe']
+        suppl['pp_DI'] = n_phase['Q'] + n_phase['Or'] + n_phase['Ab'] + n_phase['An']
+
+        partitions = partitions.fillna(0)
+        partitions = partitions.loc[:, (partitions != 0).any(axis=0)]
+        free = free.fillna(0)
+        free = free.round(to_round)
+        free = free.loc[:, (free != 0).any(axis=0)]
+        suppl = suppl.round(to_round)
+        suppl = suppl.fillna(0)
+
+        if self.verbose:
+            print(">>> Final compositions (wt %)")
+            print(partitions.to_string())
+            print(">>> Free phases (mol)")
+            print(free.to_string())
+            print(">>> Supplementary data")
+            print(suppl.to_string())
+
+        return partitions, (free, suppl)
diff --git a/georunes/ternary/__init__.py b/georunes/modmin/optim/__init__.py
similarity index 100%
rename from georunes/ternary/__init__.py
rename to georunes/modmin/optim/__init__.py
diff --git a/georunes/modmin/optim/base.py b/georunes/modmin/optim/base.py
new file mode 100644
index 0000000..d14f608
--- /dev/null
+++ b/georunes/modmin/optim/base.py
@@ -0,0 +1,218 @@
+import numpy as np
+from numpy import linalg
+from random import random
+
+
+def is_in_bounds(partition, max_minerals_prop, min_minerals_prop):
+    if any(partition[i] > max_minerals_prop[i] or partition[i] < min_minerals_prop[i] for i in range(len(partition))):
+        return False
+    return True
+
+
+def random_part_in_hypercube(current_part, semiedge, max_minerals_prop=None, min_minerals_prop=None, total=1,
+                             max_iter=100, verbose=0):
+    nb_minerals = len(current_part)
+    if not min_minerals_prop:
+        min_minerals_prop = [0] * nb_minerals
+    if not max_minerals_prop:
+        max_minerals_prop = [1] * nb_minerals
+    if len(max_minerals_prop) != nb_minerals or len(min_minerals_prop) != nb_minerals:
+        raise Exception("Wrong dimensions for max or min mineral compositions")
+
+    it = 0
+    prop = random()  # A random factor to allow smaller steps during the research of new values
+    while it < max_iter:
+        it += 1
+        dec = np.random.uniform(-semiedge * prop, semiedge * prop, nb_minerals - 1)
+        dec = np.append(dec, total - sum(current_part) - sum(dec))  # To make sure that sum(new_partition) = total
+        new_partition = [current_part[i] + dec[i] for i in range(nb_minerals)]
+        if is_in_bounds(new_partition, max_minerals_prop, min_minerals_prop):
+            # The new composition is found
+            return new_partition
+
+        if verbose > 1 and it == max_iter:
+            print("WARNING : Maximum iteration reached in the research of a random composition.")
+
+    return current_part
+
+
+def random_part_with_bounds(nb_minerals, max_minerals_prop=None, min_minerals_prop=None, total=1, verbose=0):
+    if not min_minerals_prop:
+        min_minerals_prop = [0] * nb_minerals
+    if not max_minerals_prop:
+        max_minerals_prop = [1] * nb_minerals
+
+    if len(max_minerals_prop) != nb_minerals or len(min_minerals_prop) != nb_minerals:
+        raise Exception("Wrong dimensions for max or min mineral compositions")
+
+    # Loop
+    found = False
+    k = 0
+    while not found:
+        k += 1
+        new_partition = []
+        for i in range(nb_minerals):
+            new_partition.append(min_minerals_prop[i] + (max_minerals_prop[i] - min_minerals_prop[i]) * random())
+        new_total = sum(new_partition)
+        found = True
+        if new_total > total:
+            new_partition = [i * total / new_total for i in new_partition]
+            # Check if minima conditions are respected
+            for j in range(nb_minerals):
+                if new_partition[j] < min_minerals_prop[j]:
+                    found = False
+                    break
+        else:
+            new_partition = [i * total / new_total for i in new_partition]
+            for j in range(nb_minerals):
+                if new_partition[j] > max_minerals_prop[j]:
+                    found = False
+                    break
+
+        if k > 1000:
+            raise Exception("The algorithm struggles to find a coherent random partition. Please check the entry data.")
+    if verbose > 2:
+        print("Partition found after", k, "iterations.")
+
+    return new_partition
+
+
+def random_part(nb_minerals, verbose=0):
+    return random_part_with_bounds(nb_minerals, verbose=verbose)
+
+
+class BaseOptimizer:
+    def __init__(self, verbose=None):
+        if isinstance(verbose, (int, bool)):
+            self.verbose = int(verbose)
+        else:
+            self.verbose = 0
+
+    def deviation(self, part_a, part_b):
+        if self.dist_func in ("euclidian", "L2-norm"):
+            dist = np.linalg.norm(part_a - part_b)
+        elif self.dist_func == "max":
+            dist = np.linalg.norm(part_a - part_b, ord=np.inf)
+        elif self.dist_func == "min":
+            dist = np.linalg.norm(part_a - part_b, ord=-np.inf)
+        elif self.dist_func == "MAE":  # Mean absolute error
+            dist = np.absolute(part_a - part_b).mean()
+        elif self.dist_func == "MSE":  # Mean square error
+            dist = np.square(part_a - part_b).mean()
+        elif self.dist_func == "RMSE":  # Root mean square error
+            dist = np.sqrt(np.square(part_a - part_b).mean())
+        elif self.dist_func == "SMAPE":  # Symmetric mean absolute percentage error
+            sum_abs_a_b = np.abs(part_a) + np.abs(part_b)
+            terms = np.divide(2 * np.abs(part_a - part_b), sum_abs_a_b, out=np.zeros_like(part_a),
+                              where=sum_abs_a_b != 0)  # Put division to 0 when A = B = 0
+            dist = 100 * terms.mean()
+        else:
+            raise Exception("Unknown parameter for distance function.")
+        return dist
+
+
+class Optimizer(BaseOptimizer):
+    def __init__(self, **kwargs):
+        BaseOptimizer.__init__(self, **kwargs)
+        # self.dist_func defined in children classes
+
+    def prepare_data(self, raw_data, skip_cols, raw_minerals_data, ignore_oxides):
+        raw_data = raw_data.fillna(0)
+        if 'Sum' in raw_data.keys():
+            raw_data = raw_data.rename(columns={'Sum': 'Total'})
+        if 'Total' in raw_data.keys():
+            self.init_total = raw_data["Total"]
+        else:
+            raise Exception("Total column not found.")
+
+        if ignore_oxides:
+            el_ox_to_drop = [*ignore_oxides, 'Total']
+        else:
+            el_ox_to_drop = ['Total', ]
+
+        for ox in el_ox_to_drop:
+            if ox in raw_minerals_data.keys():
+                raw_minerals_data = raw_minerals_data.drop(columns=ox)
+            if ox in raw_data.keys():
+                raw_data = raw_data.drop(columns=ox)
+
+        raw_minerals_data = raw_minerals_data.fillna(0)
+        raw_minerals_data = raw_minerals_data.set_index(raw_minerals_data.keys()[0])
+        self.data = raw_data.iloc[:, skip_cols:].copy()
+        self.list_bulk_ox = self.data.keys().tolist()
+        self.nb_oxides = len(self.list_bulk_ox)
+        raw_minerals_data = raw_minerals_data[[*list(self.list_bulk_ox)]]  # Order oxides as in source
+        self.minerals_data = raw_minerals_data.transpose()
+        self.list_minerals = self.minerals_data.keys().tolist()
+        self.nb_minerals = len(self.list_minerals)
+
+    def compute(self, raw_data, skip_cols, raw_minerals_data, to_round=4, ignore_oxides=None, ratios=None, **kwargs):
+        if self.verbose:
+            print(self.notif)
+
+        if ratios:
+            notif_cfg = ""
+            rmd_t = raw_minerals_data.transpose()
+            rmd_t.columns = rmd_t.iloc[0]
+            df_t = rmd_t[1:].copy()
+            for sol, cfg in ratios.items():
+                if len(cfg) != 2:
+                    print("WARNING : Incorrect configuration for the solution :", sol)
+                elif len(cfg[0]) != len(cfg[1]):
+                    print("WARNING : Missing data in the solution :", sol)
+                else:
+                    mins, parts = cfg[0], cfg[1]
+                    if notif_cfg:
+                        notif_cfg += " / " + str([mins, parts, ])
+                    else:
+                        notif_cfg += str([mins, parts, ])
+
+                    # Calculation of mineral solutions and replacement in minerals_data
+                    tot = sum(parts)
+                    for miner, pa in zip(mins, parts):
+                        df_t[miner] = df_t[miner] * pa / tot
+                    df_t[sol] = df_t[mins].sum(axis=1)
+                    df_t = df_t.drop(columns=mins)
+
+            if self.verbose and len(notif_cfg) > 0:
+                print("Considered mineral solutions :", notif_cfg)
+
+                # Reset dataframe as raw table
+                sub_data = df_t.transpose()
+                sub_data.reset_index(drop=False, inplace=True)
+
+                # Compute partitions
+                partitions, suppl = self._compute(raw_data, skip_cols, sub_data, to_round=to_round,
+                                                  ignore_oxides=ignore_oxides, **kwargs)
+
+                # Component reconstitution
+                for sol, cfg in ratios.items():
+                    mins, parts = cfg[0], cfg[1]
+                    tot = sum(parts)
+                    suppl[sol] = partitions[sol]
+                    for miner, pa in zip(mins, parts):
+                        partitions[miner] = partitions[sol] * pa / tot
+                    partitions = partitions.drop(columns=sol, axis=1)
+                # Move Total to the last column
+                cols = partitions.columns.tolist()
+                cols.remove("Total")
+                cols.append("Total")
+                partitions = partitions[cols]
+            else:
+                partitions, suppl = self._compute(raw_data, skip_cols, raw_minerals_data, to_round=to_round,
+                                                  ignore_oxides=ignore_oxides, **kwargs)
+
+        else:
+            partitions, suppl = self._compute(raw_data, skip_cols, raw_minerals_data, to_round=to_round,
+                                              ignore_oxides=ignore_oxides, **kwargs)
+
+        if self.verbose:
+            self.show_results(partitions, suppl)
+
+        return partitions, suppl
+
+    def show_results(self, partitions, suppl):
+        print(">>> Final compositions (wt %)")
+        print(partitions.to_string())
+        print(">>> Supplementary data")
+        print(suppl.to_string())
diff --git a/georunes/modmin/optim/bvls.py b/georunes/modmin/optim/bvls.py
new file mode 100644
index 0000000..0a200d4
--- /dev/null
+++ b/georunes/modmin/optim/bvls.py
@@ -0,0 +1,88 @@
+import numpy as np
+from pandas import DataFrame
+from scipy.optimize import lsq_linear
+from georunes.modmin.optim.base import Optimizer
+
+
+class BVLS(Optimizer):
+    def __init__(self, **kwargs):
+        Optimizer.__init__(self, **kwargs)
+        self.dist_func = "euclidian"
+        self.notif = ">>>>>> Bounded Variable Least Squares Method"
+
+    def _compute(self, raw_data, skip_cols, raw_minerals_data, to_round=4, max_minerals=None, min_minerals=None,
+                 ignore_oxides=None):
+        if self.verbose > 1:
+            print("Round digits :", to_round, " --  Oxides to ignore : " + str(ignore_oxides) if ignore_oxides else "")
+        self.prepare_data(raw_data, skip_cols, raw_minerals_data, ignore_oxides)
+        partitions = DataFrame(columns=self.list_minerals)
+        deviation_name = "deviation_" + self.dist_func
+        suppl = DataFrame(columns=[deviation_name])
+
+        for i in self.data.index:
+            if self.verbose: print(">>> Composition", i)
+            list_minerals_i = self.list_minerals.copy()
+            minerals_data_i = self.minerals_data.copy()
+            bulk_chem = self.data.iloc[i].to_numpy()
+
+            # Get maximum possible proportion for each mineral
+            max_minerals_prop = []
+            unnecessary_minerals = []
+            for mineral in self.list_minerals:
+                # Maximum proportion for each oxide equals Ox wt% in bulk chemistry / Ox wt% in mineral composition
+                max_oxides_prop = self.data.iloc[i] / minerals_data_i[mineral].transpose()
+                mineral_prop = max_oxides_prop.min()
+                if np.isinf(mineral_prop) or mineral_prop == 0:
+                    minerals_data_i = minerals_data_i.drop(columns=mineral)
+                    list_minerals_i.remove(mineral)
+                    unnecessary_minerals.append(mineral)
+                else:
+                    max_minerals_prop.append(mineral_prop)
+            nb_minerals_i = len(list_minerals_i)
+            min_minerals_prop = [0] * nb_minerals_i
+
+            # Update minerals proportions with defined settings if provided
+            if max_minerals:
+                for key, value in max_minerals.items():
+                    idx = list_minerals_i.index(key)
+                    if idx:
+                        max_minerals_prop[idx] = min(max_minerals_prop[idx], max_minerals[key])
+            if min_minerals:
+                for key, value in min_minerals.items():
+                    idx = list_minerals_i.index(key)
+                    if idx:
+                        min_minerals_prop[idx] = max(min_minerals_prop[idx], min_minerals[key])
+            if sum(max_minerals_prop) < 1.05:  # 5% tolerance
+                raise Exception("The authorized maximum proportions do not allow to fill the composition to 100. "
+                                "Check the entry data.")
+
+            if self.verbose and unnecessary_minerals:
+                print("Unnecessary minerals :", *unnecessary_minerals)
+
+            # Direct calculation of the result
+            result = lsq_linear(minerals_data_i, bulk_chem, (min_minerals_prop, max_minerals_prop), method='bvls',
+                                verbose=min(self.verbose, 2))
+            dict_res = {list_minerals_i[i]: result.x[i] for i in range(len(list_minerals_i))}
+
+            idx_new_row = len(partitions)
+            for miner, prop in dict_res.items():
+                partitions.loc[idx_new_row, miner] = prop
+            partitions = partitions.fillna(0)
+            partitions.iloc[idx_new_row] = 100 * partitions.iloc[idx_new_row]
+            partitions = partitions.round(to_round)
+
+            deviation = self.deviation(bulk_chem, np.dot(self.minerals_data, partitions.iloc[idx_new_row] / 100))
+            suppl.loc[idx_new_row, deviation_name] = deviation
+
+            if self.verbose:
+                print("Solution", idx_new_row)
+                print(partitions.loc[idx_new_row].to_dict())
+                print("Corresponding composition")
+                chem = np.dot(self.minerals_data, partitions.loc[idx_new_row] / 100).round(to_round)
+                print([str(self.list_bulk_ox[i]) + " : " + str(chem[i]) for i in range(self.nb_oxides)])
+                print("Deviation :", round(deviation, to_round), "%" if self.dist_func == "SMAPE" else "", "\n------")
+
+        partitions["Total"] = partitions.sum(axis=1).round(to_round)
+        suppl["Diff_total"] = partitions["Total"] - self.init_total
+
+        return partitions, suppl
diff --git a/georunes/modmin/optim/gd.py b/georunes/modmin/optim/gd.py
new file mode 100644
index 0000000..db26b4c
--- /dev/null
+++ b/georunes/modmin/optim/gd.py
@@ -0,0 +1,163 @@
+import numpy as np
+from pandas import DataFrame
+from georunes.modmin.optim.base import Optimizer, is_in_bounds, random_part_with_bounds
+
+
+class GradientDescent(Optimizer):
+
+    def __init__(self, dist_func="euclidian", **kwargs):
+        Optimizer.__init__(self, **kwargs)
+        self.dist_func = dist_func
+        self.notif = ">>>>>> Gradient Descent method / deviation function : " + dist_func
+
+    def grad(self, x, y, A, total=1):
+        B = A.copy()
+        B.loc[len(B.index)] = np.array([1] * len(x))  # Row added to minerals_data (A) to keep a sum equals to total
+        Bt = B.transpose()
+        y = np.append(y, total)
+        out = (np.dot(np.dot(Bt, B), x) - np.dot(Bt, y)) / (len(B.index))
+        return out
+
+    def _compute(self, raw_data, skip_cols, raw_minerals_data, to_round=4,
+                 max_minerals=None, min_minerals=None, ignore_oxides=None, max_iter=int(1e6),
+                 learn_rate=None, tolerance=1e-08, starting_partition=None, force_totals=False):
+
+        if not learn_rate:
+            learn_rate = 1 * pow(10, -to_round)
+
+        if self.verbose > 1:
+            print("Round digits :", to_round, " --  Maximum iterations :", max_iter, " --  Learn rate :", learn_rate,
+                  " --  Tolerance :", tolerance,
+                  " --  Oxides to ignore : " + str(ignore_oxides) if ignore_oxides else "")
+            if starting_partition: print("Starting partition :", starting_partition)
+
+        self.prepare_data(raw_data, skip_cols, raw_minerals_data, ignore_oxides)
+        partitions = DataFrame(columns=self.list_minerals)
+        deviation_name = "deviation_" + self.dist_func
+        suppl = DataFrame(columns=[deviation_name])
+
+        if force_totals:
+            target_totals = self.init_total
+        else:
+            target_totals = [100] * len(self.data.index)
+
+        # Get maximum proportion for each mineral
+        for i in self.data.index:
+            if self.verbose: print(">>> Composition", i)
+            list_minerals_i = self.list_minerals.copy()
+            minerals_data_i = self.minerals_data.copy()
+            bulk_chem = self.data.iloc[i].to_numpy()
+
+            # Get maximum possible proportion for each mineral
+            max_minerals_prop = []
+            unnecessary_minerals = []
+            for mineral in self.list_minerals:
+                # Maximum proportion for each oxide equals Ox wt% in bulk chemistry / ox wt% in mineral composition
+                max_oxides_prop = self.data.iloc[i] / minerals_data_i[mineral].transpose()
+                mineral_prop = max_oxides_prop.min()
+                if np.isinf(mineral_prop) or mineral_prop == 0:
+                    minerals_data_i = minerals_data_i.drop(columns=mineral)
+                    list_minerals_i.remove(mineral)
+                    unnecessary_minerals.append(mineral)
+                else:
+                    max_minerals_prop.append(mineral_prop)
+            nb_minerals_i = len(list_minerals_i)
+            min_minerals_prop = [0] * nb_minerals_i
+
+            # Update minerals proportions with defined settings if provided
+            if max_minerals:
+                for key, value in max_minerals.items():
+                    idx = list_minerals_i.index(key)
+                    if idx:
+                        max_minerals_prop[idx] = min(max_minerals_prop[idx], max_minerals[key])
+            if min_minerals:
+                for key, value in min_minerals.items():
+                    idx = list_minerals_i.index(key)
+                    if idx:
+                        min_minerals_prop[idx] = max(min_minerals_prop[idx], min_minerals[key])
+            if sum(max_minerals_prop) < 1.05:  # 5% tolerance
+                raise Exception("The authorized maximum proportions do not allow to fill the composition to 100. "
+                                "Check the entry data.")
+
+            if self.verbose and unnecessary_minerals:
+                print("Unnecessary minerals :", *unnecessary_minerals)
+
+            # Gradient descent
+            if self.verbose: print("Start calculations. Expected number of iterations :", max_iter)
+
+            # Starting value
+            if isinstance(starting_partition, dict):
+                if sum(starting_partition.values()) != 100:
+                    print("WARNING : The starting partition does not complete to 100. The calculations will be "
+                          "started with a random composition.")
+                    candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                        verbose=self.verbose)
+                else:
+                    # Check if all minerals in the config are present in minerals_data_i
+                    if all(el in list_minerals_i for el in starting_partition.keys()):
+                        candidate = []
+                        for key in list_minerals_i:
+                            val = starting_partition[key] / 100. if key in starting_partition.keys() else 0
+                            candidate.append(val)
+                        if not is_in_bounds(candidate, max_minerals_prop, min_minerals_prop):
+                            print("WARNING : The starting partition does not respect the provided or the calculated "
+                                  "bounds data : Minerals =", list_minerals_i, ", Max =", max_minerals_prop, ", Min =",
+                                  min_minerals_prop, "The calculations will be started with a random composition.")
+                            candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                                verbose=self.verbose)
+                    else:
+                        print("WARNING : Some minerals in starting partition are not present in mineral chemistry "
+                              "data. The calculations will be started with a random composition.")
+                        candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                            verbose=self.verbose)
+            else:
+                candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                    verbose=self.verbose)
+
+            if self.verbose: print("Starting partition", dict(zip(list_minerals_i, candidate)))
+
+            # Loop
+            for k in range(max_iter):
+                diff = -learn_rate * np.array(
+                    self.grad(candidate, bulk_chem, minerals_data_i, total=target_totals[i] / 100))
+                new_candidate = candidate + diff
+                new_candidate = np.where(new_candidate < min_minerals_prop, min_minerals_prop, new_candidate)
+                new_candidate = np.where(new_candidate > max_minerals_prop, max_minerals_prop, new_candidate)
+                clean_diff = new_candidate - candidate
+
+                if np.all(np.abs(clean_diff) <= tolerance):
+                    if self.verbose: print("Tolerance condition reached after iteration", k)
+                    break
+
+                candidate = new_candidate
+                if self.verbose > 2:
+                    print("New solution at iteration", k)
+                    print(dict(zip(list_minerals_i, candidate)))
+                    corresp_chem = np.dot(minerals_data_i, candidate).round(decimals=to_round)
+                    print("New deviation :", self.deviation(bulk_chem, corresp_chem),
+                          "%" if self.dist_func == "SMAPE" else "")
+
+            dict_res = {list_minerals_i[i]: round(candidate[i], to_round) for i in range(len(list_minerals_i))}
+
+            idx_new_row = len(partitions)
+            for miner, prop in dict_res.items():
+                partitions.loc[idx_new_row, miner] = prop
+            partitions = partitions.fillna(0)
+            partitions.iloc[idx_new_row] = 100 * partitions.iloc[idx_new_row]
+            partitions = partitions.round(to_round)
+
+            deviation = self.deviation(bulk_chem, np.dot(self.minerals_data, partitions.iloc[idx_new_row] / 100.))
+            suppl.loc[idx_new_row, deviation_name] = deviation
+
+            if self.verbose:
+                print("Solution", idx_new_row)
+                print(partitions.loc[idx_new_row].to_dict())
+                print("Corresponding composition")
+                chem = np.dot(self.minerals_data, partitions.loc[idx_new_row] / 100.).round(to_round)
+                print([str(self.list_bulk_ox[i]) + " : " + str(chem[i]) for i in range(self.nb_oxides)])
+                print("Deviation :", round(deviation, to_round), "%" if self.dist_func == "SMAPE" else "", "\n------")
+
+        partitions["Total"] = partitions.sum(axis=1).round(to_round)
+        suppl["Diff_total"] = partitions["Total"] - self.init_total
+
+        return partitions, suppl
diff --git a/georunes/modmin/optim/nnls.py b/georunes/modmin/optim/nnls.py
new file mode 100644
index 0000000..e9500ef
--- /dev/null
+++ b/georunes/modmin/optim/nnls.py
@@ -0,0 +1,66 @@
+import numpy as np
+from pandas import DataFrame
+from scipy.optimize import nnls
+from georunes.modmin.optim.base import Optimizer
+
+
+class NNLS(Optimizer):
+    def __init__(self, verbose=False, **kwargs):
+        Optimizer.__init__(self, verbose=verbose, **kwargs)
+        self.dist_func = "euclidian"
+        self.notif = ">>>>>> Non-Negative Least Squares method"
+
+    def _compute(self, raw_data, skip_cols, raw_minerals_data, to_round=4, ignore_oxides=None,
+                 max_iter=None):
+        if self.verbose > 1:
+            print("Round digits :", to_round, " --  Maximum iterations :", max_iter,
+                  " --  Oxides to ignore : " + str(ignore_oxides) if ignore_oxides else "")
+        self.prepare_data(raw_data, skip_cols, raw_minerals_data, ignore_oxides)
+
+        partitions = DataFrame(columns=self.list_minerals)
+        deviation_name = "deviation_" + self.dist_func
+        suppl = DataFrame(columns=[deviation_name])
+        for i in self.data.index:
+            if self.verbose: print(">>> Composition", i)
+            list_minerals_i = self.list_minerals.copy()
+            minerals_data_i = self.minerals_data.copy()
+            bulk_chem = self.data.iloc[i].to_numpy()
+
+            # Ignore minerals missing some oxides
+            unnecessary_minerals = []
+            for ox in self.list_bulk_ox:
+                if self.data.loc[i, ox] == 0:
+                    for mineral in minerals_data_i:
+                        if minerals_data_i.loc[ox, mineral] > 0:
+                            minerals_data_i = minerals_data_i.drop(columns=mineral)
+                            list_minerals_i.remove(mineral)
+                            unnecessary_minerals.append(mineral)
+            if self.verbose and unnecessary_minerals:
+                print("Unnecessary minerals :", *unnecessary_minerals)
+
+            # Direct calculation of the result
+            result = nnls(minerals_data_i, bulk_chem, max_iter)
+            dict_res = {list_minerals_i[i]: result[0][i] for i in range(len(list_minerals_i))}
+
+            idx_new_row = len(partitions)
+            for miner, prop in dict_res.items():
+                partitions.loc[idx_new_row, miner] = prop
+            partitions = partitions.fillna(0)
+            partitions.iloc[idx_new_row] = 100 * partitions.iloc[idx_new_row]
+            partitions = partitions.round(to_round)
+
+            deviation = self.deviation(bulk_chem, np.dot(self.minerals_data, partitions.iloc[idx_new_row] / 100))
+            suppl.loc[idx_new_row, deviation_name] = deviation
+
+            if self.verbose:
+                print("Solution", idx_new_row)
+                print(partitions.loc[idx_new_row].to_dict())
+                print("Corresponding composition")
+                chem = np.dot(self.minerals_data, partitions.loc[idx_new_row] / 100).round(to_round)
+                print([str(self.list_bulk_ox[i]) + " : " + str(chem[i]) for i in range(self.nb_oxides)])
+                print("Deviation :", round(deviation, to_round), "%" if self.dist_func == "SMAPE" else "", "\n------")
+
+        partitions["Total"] = partitions.sum(axis=1).round(to_round)
+        suppl["Diff_total"] = partitions["Total"] - self.init_total
+
+        return partitions, suppl
diff --git a/georunes/modmin/optim/randsearch.py b/georunes/modmin/optim/randsearch.py
new file mode 100644
index 0000000..b4881f9
--- /dev/null
+++ b/georunes/modmin/optim/randsearch.py
@@ -0,0 +1,172 @@
+import numpy as np
+from pandas import DataFrame
+from georunes.modmin.optim.base import Optimizer, is_in_bounds, random_part_with_bounds, random_part_in_hypercube
+
+
+class RandomSearch(Optimizer):
+
+    def __init__(self, dist_func="euclidian", **kwargs):
+        Optimizer.__init__(self, **kwargs)
+        self.dist_func = dist_func
+        self.notif = ">>>>>> Random research Method / deviation function : " + dist_func
+
+    def _compute(self, raw_data, skip_cols, raw_minerals_data, to_round=4,
+                 max_minerals=None, min_minerals=None, ignore_oxides=None, max_iter=1000000,
+                 limit_deviation=1e-03, search_semiedge=1, scale_semiedge=None,
+                 starting_partition=None, force_totals=False):
+
+        if self.verbose > 1:
+            print("Round digits :", to_round, " --  Maximum iterations :", max_iter, " --  Limit deviation:",
+                  limit_deviation,
+                  " --  Semiegde of hypercube research :", search_semiedge,
+                  " --  Oxides to ignore : " + str(ignore_oxides) if ignore_oxides else "")
+            if starting_partition: print("Starting partition :", starting_partition)
+
+        self.prepare_data(raw_data, skip_cols, raw_minerals_data, ignore_oxides)
+        partitions = DataFrame(columns=self.list_minerals)
+        deviation_name = "deviation_" + self.dist_func
+        suppl = DataFrame(columns=[deviation_name])
+
+        if force_totals:
+            target_totals = self.init_total
+        else:
+            target_totals = [100] * len(self.data.index)
+
+        if not scale_semiedge or scale_semiedge < 0 or scale_semiedge > 1:
+            if self.verbose > 0:
+                print("WARNING : The scale must be superior to 0 and inferior or equals to 1. Value set to 1.")
+            scale_semiedge = 1
+
+        # Get maximum proportion for each mineral
+        for i in self.data.index:
+            if self.verbose: print(">>> Composition", i)
+            list_minerals_i = self.list_minerals.copy()
+            minerals_data_i = self.minerals_data.copy()
+            bulk_chem = self.data.iloc[i].to_numpy()
+            search_semiedge_i = search_semiedge
+
+            # Get maximum possible proportion for each mineral
+            max_minerals_prop = []
+            unnecessary_minerals = []
+            for mineral in self.list_minerals:
+                # Maximum proportion for each oxide equals Ox wt% in bulk chemistry / ox wt% in mineral composition
+                max_oxides_prop = self.data.iloc[i] / minerals_data_i[mineral].transpose()
+                mineral_prop = max_oxides_prop.min()
+                if np.isinf(mineral_prop) or mineral_prop == 0:
+                    minerals_data_i = minerals_data_i.drop(columns=mineral)
+                    list_minerals_i.remove(mineral)
+                    unnecessary_minerals.append(mineral)
+                else:
+                    max_minerals_prop.append(mineral_prop)
+            nb_minerals_i = len(list_minerals_i)
+            min_minerals_prop = [0] * nb_minerals_i
+
+            # Update minerals proportions with defined settings if provided
+            if max_minerals:
+                for key, value in max_minerals.items():
+                    idx = list_minerals_i.index(key)
+                    if idx:
+                        max_minerals_prop[idx] = min(max_minerals_prop[idx], max_minerals[key])
+            if min_minerals:
+                for key, value in min_minerals.items():
+                    idx = list_minerals_i.index(key)
+                    if idx:
+                        min_minerals_prop[idx] = max(min_minerals_prop[idx], min_minerals[key])
+            if sum(max_minerals_prop) < 1.05:  # 5% tolerance
+                raise Exception("The authorized maximum proportions do not allow to fill the composition to 100. "
+                                "Check the entry data.")
+
+            if self.verbose and unnecessary_minerals:
+                print("Unnecessary minerals :", *unnecessary_minerals)
+
+            # Random calculation
+            if self.verbose: print("Start calculations. Expected number of iterations :", max_iter)
+            min_deviation = 1e10
+            # Starting value
+            if isinstance(starting_partition, dict):
+                if sum(starting_partition.values()) != 100:
+                    print("WARNING : The starting partition does not complete to 100. The calculations will be "
+                          "started with a random composition.")
+                    candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                        verbose=self.verbose)
+                else:
+                    # Check if all minerals in the config are present in minerals_data_i
+                    if all(el in list_minerals_i for el in starting_partition.keys()):
+                        candidate = []
+                        for key in list_minerals_i:
+                            val = starting_partition[key] / 100 if key in starting_partition.keys() else 0
+                            candidate.append(val)
+                        if not is_in_bounds(candidate, max_minerals_prop, min_minerals_prop):
+                            print("WARNING : The starting partition does not respect the provided and calculated "
+                                  "bounds data : Minerals ,", list_minerals_i, "Max ,", max_minerals_prop, "Min ,",
+                                  min_minerals_prop, "The calculations will be started with a random composition.")
+                            candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                                verbose=self.verbose)
+                    else:
+                        print("WARNING : Some minerals in starting partition are not present in mineral chemistry "
+                              "data. The calculations will be started with a random composition.")
+                        candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                            verbose=self.verbose)
+            else:
+                candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                    verbose=self.verbose)
+
+            if self.verbose: print("Starting partition", dict(zip(list_minerals_i, candidate)))
+
+            # Loop
+            for k in range(max_iter):
+                if search_semiedge_i > 1:
+                    raise Exception("The variable search_semiedge_i must be inferior or equal to 1.")
+                elif search_semiedge_i < 1 and candidate:
+                    new_candidate = random_part_in_hypercube(candidate, search_semiedge_i, max_minerals_prop,
+                                                             min_minerals_prop, total=target_totals[i] / 100,
+                                                             verbose=self.verbose)
+                    if set(new_candidate) == set(candidate) and search_semiedge_i < 1:
+                        search_semiedge_i = search_semiedge_i * scale_semiedge
+                        if self.verbose > 0:
+                            print("In iteration", k, ", search_semiedge_i changed to", search_semiedge_i)
+
+                else:
+                    new_candidate = random_part_with_bounds(nb_minerals_i, max_minerals_prop, min_minerals_prop,
+                                                            total=target_totals[i] / 100, verbose=self.verbose)
+
+                corresp_chem = np.dot(minerals_data_i, new_candidate).round(decimals=to_round)
+                dist = self.deviation(bulk_chem, corresp_chem)
+                if dist < min_deviation:
+                    min_deviation = dist
+                    candidate = new_candidate
+                    if min_deviation < limit_deviation:
+                        if self.verbose: print("Bottom deviation condition reached after iteration", k)
+                        break
+                    if self.verbose > 1:
+                        print("Better solution at iteration", k)
+                        print(dict(zip(list_minerals_i, candidate)))
+                        print("New deviation :", min_deviation, "%" if self.dist_func == "SMAPE" else "")
+
+                if self.verbose and k == max_iter - 1:
+                    print("Max iterations reached")
+
+            dict_res = {list_minerals_i[i]: round(candidate[i], to_round) for i in range(len(list_minerals_i))}
+
+            idx_new_row = len(partitions)
+            for miner, prop in dict_res.items():
+                partitions.loc[idx_new_row, miner] = prop
+            partitions = partitions.fillna(0)
+            partitions.iloc[idx_new_row] = 100 * partitions.iloc[idx_new_row]
+            partitions = partitions.round(to_round)
+
+            deviation = self.deviation(bulk_chem, np.dot(self.minerals_data, partitions.iloc[idx_new_row] / 100))
+            suppl.loc[idx_new_row, deviation_name] = deviation
+
+            if self.verbose:
+                print("Solution", idx_new_row)
+                print(partitions.loc[idx_new_row].to_dict())
+                print("Corresponding composition")
+                chem = np.dot(self.minerals_data, partitions.loc[idx_new_row] / 100).round(to_round)
+                print([str(self.list_bulk_ox[i]) + " : " + str(chem[i]) for i in range(self.nb_oxides)])
+                print("Deviation :", round(deviation, to_round), "%" if self.dist_func == "SMAPE" else "", "\n------")
+
+        partitions["Total"] = partitions.sum(axis=1).round(to_round)
+        suppl["Diff_total"] = partitions["Total"] - self.init_total
+
+        return partitions, suppl
diff --git a/georunes/plot/__init__.py b/georunes/plot/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/georunes/base/base.py b/georunes/plot/base.py
similarity index 100%
rename from georunes/base/base.py
rename to georunes/plot/base.py
diff --git a/georunes/plot/binary/__init__.py b/georunes/plot/binary/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/georunes/binary/debonlefortba.py b/georunes/plot/binary/debonlefortba.py
similarity index 93%
rename from georunes/binary/debonlefortba.py
rename to georunes/plot/binary/debonlefortba.py
index 0f885e1..c93cab8 100644
--- a/georunes/binary/debonlefortba.py
+++ b/georunes/plot/binary/debonlefortba.py
@@ -1,11 +1,14 @@
 import numpy as np
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.chemistry import val_ox_to_mc
 from georunes.tools.language import get_translator
 from georunes.tools.plotting import get_spline
 
 
+# Debon, F. and Le Fort, P., 1983. A chemical–mineralogical classification of common plutonic rocks and associations.
+# Earth and Environmental Science Transactions of the Royal Society of Edinburgh, 73(3), pp.135-149.
+
 class DiagramBA(DiagramBase, ArrowDrawer, LegendDrawer):
     def __init__(self, datasource, title="B-A multicationic classification (Debon and Le Fort 1983)",
                  annotation=None, decor_set="Debon", legend_loc="upper right", **kwargs):
diff --git a/georunes/binary/debonlefortpq.py b/georunes/plot/binary/debonlefortpq.py
similarity index 94%
rename from georunes/binary/debonlefortpq.py
rename to georunes/plot/binary/debonlefortpq.py
index a5eebb8..7c6e47c 100644
--- a/georunes/binary/debonlefortpq.py
+++ b/georunes/plot/binary/debonlefortpq.py
@@ -1,9 +1,12 @@
 import numpy as np
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import ArrowDrawer, LegendDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import ArrowDrawer, LegendDrawer
 from georunes.tools.chemistry import val_ox_to_mc
 
 
+# Debon, F. and Le Fort, P., 1983. A chemical–mineralogical classification of common plutonic rocks and associations.
+# Earth and Environmental Science Transactions of the Royal Society of Edinburgh, 73(3), pp.135-149.
+
 class DiagramPQ(DiagramBase, ArrowDrawer, LegendDrawer):
     def __init__(self, datasource, title="P-Q multicationic classification (Debon and Le Fort 1983)",
                  padding={"bottom": 0.20}, annotation=None, **kwargs):
diff --git a/georunes/binary/delaroche.py b/georunes/plot/binary/delaroche.py
similarity index 97%
rename from georunes/binary/delaroche.py
rename to georunes/plot/binary/delaroche.py
index 82eaa39..bdec48b 100644
--- a/georunes/binary/delaroche.py
+++ b/georunes/plot/binary/delaroche.py
@@ -1,10 +1,14 @@
 import numpy as np
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.chemistry import val_ox_to_mc
 from georunes.tools.language import get_translator
 
 
+# De la Roche, H.D., Leterrier, J.T., Grandclaude, P. and Marchal, M., 1980. A classification of volcanic and
+# plutonic rocks using R1R2-diagram and major-element analyses—its relationships with current nomenclature. Chemical
+# geology, 29(1-4), pp.183-210.
+
 class DiagramR1R2(DiagramBase, ArrowDrawer, LegendDrawer):
     def __init__(self, datasource, title="R1-R2 multicationic classification (De La Roche et al. 1980)",
                  padding={"bottom": 0.20}, annotation=None, **kwargs):
diff --git a/georunes/plot/binary/frost.py b/georunes/plot/binary/frost.py
new file mode 100644
index 0000000..c8b6a17
--- /dev/null
+++ b/georunes/plot/binary/frost.py
@@ -0,0 +1,175 @@
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import ArrowDrawer, LegendDrawer
+from georunes.tools.language import format_chemical_formula as _fml, get_translator
+
+
+# Frost, B.R., Barnes, C.G., Collins, W.J., Arculus, R.J., Ellis, D.J. and Frost, C.D., 2001. A geochemical
+# classification for granitic rocks. Journal of petrology, 42(11), pp.2033-2048.
+
+class DiagramFrostSiFeNb(DiagramBase, ArrowDrawer, LegendDrawer):
+    def __init__(self, datasource, title=_fml("SiO2 vs Fe# Frost diagram"), **kwargs):
+        DiagramBase.__init__(self, datasource=datasource, title=title, **kwargs)
+        if 'FeO' not in self.data.keys():
+            raise Exception("FeO data missing. Check the table.")
+        self.xlabel = _fml("SiO2 (wt %)")
+        self.ylabel = "FeO/(FeO + MgO)"
+
+    def set_decoration(self):
+        if not self.no_title:
+            self.ax.set_title(self.title, size=self.title_fs)
+        _ = get_translator(self.lang_cfg)
+        self.ax.set_xlim(50, 80)
+        self.ax.set_ylim(0.4, 1)
+
+        self.ax.set_xlabel(self.xlabel, fontsize=self.fontsize)
+        self.ax.set_ylabel(self.ylabel, fontsize=self.fontsize)
+
+        x1 = (50, 80)
+        y1 = (0.676, 0.814)
+
+        self.ax.plot(x1, y1, self.decor_line_col, linewidth=0.6)
+
+        self.ax.text(52, 0.8, _("ferroan"), size="small")
+        self.ax.text(52, 0.6, _("magnesian"), size="small")
+
+    def plot(self):
+        DiagramBase.plot_config(self)
+        self.set_decoration()
+
+        # Categorize by group
+        groups = self.data.groupby(self.group_name)
+        for name, group in groups:
+
+            if self.is_group_allowed(name):
+
+                fenb = group["FeO"] / (group["FeO"] + group["MgO"])
+                si = (group["SiO2"])
+
+                zorder = 4
+                if self.drawing_order:
+                    zorder = list(group[self.drawing_order])[0]
+
+                self.ax.scatter(si, fenb, edgecolors=group["color"],
+                                marker=list(group["marker"])[0], label=name, facecolors=group["color"],
+                                s=self.markersize,
+                                alpha=0.7, zorder=zorder)
+
+        self.plot_arrows()
+        self.plot_legend()
+
+
+class DiagramFrostSiFeTotNb(DiagramBase, ArrowDrawer, LegendDrawer):
+    def __init__(self, datasource, title=_fml("SiO2 vs Fe* Frost diagram"), calc_total_Fe=False, **kwargs):
+        DiagramBase.__init__(self, datasource=datasource, title=title, **kwargs)
+        if 'FeOt' not in self.data.keys():
+            if calc_total_Fe:
+                if 'FeO' not in self.data.keys():
+                    raise Exception("FeO data missing. Cannot calculate total FeO.")
+                elif 'Fe2O3' not in self.data.keys():
+                    raise Exception("Fe2O3 data missing. Cannot calculate total FeO.")
+                else:
+                    self.data['FeOt'] = self.data['FeO'] + self.data["Fe2O3"] * 0.8998
+            else:
+                raise Exception("FeOt data missing. Check the table.")
+        self.xlabel = _fml("SiO2 (wt %)")
+        self.ylabel = "FeO$^{tot}$/(FeO$^{tot}$ + MgO)"
+
+    def set_decoration(self):
+        if not self.no_title:
+            self.ax.set_title(self.title, size=self.title_fs)
+        _ = get_translator(self.lang_cfg)
+        self.ax.set_xlim(50, 80)
+        self.ax.set_ylim(0.4, 1)
+
+        self.ax.set_xlabel(self.xlabel, fontsize=self.fontsize)
+        self.ax.set_ylabel(self.ylabel, fontsize=self.fontsize)
+
+        x1 = (50, 80)
+        y1 = (0.716, 0.854)
+
+        self.ax.plot(x1, y1, self.decor_line_col, linewidth=0.6)
+
+        self.ax.text(52, 0.8, _("ferroan"), size="small")
+        self.ax.text(52, 0.6, _("magnesian"), size="small")
+
+    def plot(self):
+        DiagramBase.plot_config(self)
+        self.set_decoration()
+
+        # Categorize by group
+        groups = self.data.groupby(self.group_name)
+        for name, group in groups:
+
+            if self.is_group_allowed(name):
+
+                fenb = group["FeOt"] / (group["FeOt"] + group["MgO"])
+                si = (group["SiO2"])
+
+                zorder = 4
+                if self.drawing_order:
+                    zorder = list(group[self.drawing_order])[0]
+
+                self.ax.scatter(si, fenb, edgecolors=group["color"],
+                                marker=list(group["marker"])[0], label=name, facecolors=group["color"],
+                                s=self.markersize,
+                                alpha=0.7, zorder=zorder)
+
+        self.plot_arrows()
+        self.plot_legend()
+
+
+class DiagramFrostSiMALI(DiagramBase, ArrowDrawer, LegendDrawer):
+    def __init__(self, datasource, title=_fml("SiO2 vs modified alkali-lime index Frost diagram"), **kwargs):
+        DiagramBase.__init__(self, datasource=datasource, title=title, **kwargs)
+        self.xlabel = _fml("SiO2 (wt %)")
+        self.ylabel = _fml("Na2O + K2O - CaO (wt %)")
+
+    def set_decoration(self):
+        if not self.no_title:
+            self.ax.set_title(self.title, size=self.title_fs)
+        _ = get_translator(self.lang_cfg)
+        self.ax.set_xlim(50, 80)
+        self.ax.set_ylim(-8, 12)
+
+        self.ax.set_xlabel(self.xlabel, fontsize=self.fontsize)
+        self.ax.set_ylabel(self.ylabel, fontsize=self.fontsize)
+
+        lx = [i for i in range(50, 75, 5)]
+
+        ly1 = [-41.86 + 1.112 * x - 0.00572 * x * x for x in lx]
+        ly2 = [-44.72 + 1.094 * x - 0.00527 * x * x for x in lx]
+        ly3 = [-45.36 + 1.0043 * x - 0.00427 * x * x for x in lx]
+
+        self.ax.plot(lx, ly1, self.decor_line_col, linewidth=0.6)
+        self.ax.plot(lx, ly2, self.decor_line_col, linewidth=0.6)
+        self.ax.plot(lx, ly3, self.decor_line_col, linewidth=0.6)
+
+        self.ax.text(51, 2.5, _("alkalic"), size="small", )
+        self.ax.text(51, -1.3, _("alkali-calcic"), size="small", rotation=30)
+        self.ax.text(51, -3.5, _("calc-alkalic"), size="small", rotation=30)
+        self.ax.text(51, -6.5, _("calcic"), size="small", )
+
+    def plot(self):
+        DiagramBase.plot_config(self)
+        self.set_decoration()
+
+        # Categorize by group
+        groups = self.data.groupby(self.group_name)
+        for name, group in groups:
+
+            if self.is_group_allowed(name):
+
+                mali = group["Na2O"] + group["K2O"] - group["CaO"]
+                si = (group["SiO2"])
+
+                zorder = 4
+                if self.drawing_order:
+                    zorder = list(group[self.drawing_order])[0]
+
+                self.ax.scatter(si, mali, edgecolors=group["color"],
+                                marker=list(group["marker"])[0], label=name, facecolors=group["color"],
+                                s=self.markersize,
+                                alpha=0.7, zorder=zorder)
+
+        self.plot_arrows()
+        self.plot_legend()
diff --git a/georunes/binary/pearce.py b/georunes/plot/binary/pearce.py
similarity index 97%
rename from georunes/binary/pearce.py
rename to georunes/plot/binary/pearce.py
index 7b49f73..b4a8c32 100644
--- a/georunes/binary/pearce.py
+++ b/georunes/plot/binary/pearce.py
@@ -1,10 +1,12 @@
 import numpy as np
-
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.language import get_translator
 
 
+# Pearce, J.A., Harris, N.B.W. and Tindle, A.G., 1984. Trace Element Discrimination Diagrams for the Tectonic
+# Interpretation of Granitic Rocks. Journal of Petrology, 25(4), pp.956-983.
+
 class PearceBase(DiagramBase, ArrowDrawer, LegendDrawer):
     def __init__(self, datasource,
                  x_formatter=None, y_formatter=None,
diff --git a/georunes/binary/peccetaylor.py b/georunes/plot/binary/peccetaylor.py
similarity index 90%
rename from georunes/binary/peccetaylor.py
rename to georunes/plot/binary/peccetaylor.py
index 333e5d0..f2471b9 100644
--- a/georunes/binary/peccetaylor.py
+++ b/georunes/plot/binary/peccetaylor.py
@@ -1,9 +1,12 @@
 import numpy as np
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.language import format_chemical_formula as _fml, get_translator
 
 
+# Peccerillo, A. and Taylor, S.R., 1976. Geochemistry of Eocene calc-alkaline volcanic rocks from the Kastamonu area,
+# northern Turkey. Contributions to mineralogy and petrology, 58(1), pp.63-81.
+
 class DiagramPecceTaylor(DiagramBase, ArrowDrawer, LegendDrawer):
     def __init__(self, datasource, title=_fml("K2O vs SiO2") + "classification diagram (Peccerillo and Taylor 1976)",
                  padding={"bottom": 0.20}, annotation=None, **kwargs):
diff --git a/georunes/binary/ratiosversus.py b/georunes/plot/binary/ratiosversus.py
similarity index 97%
rename from georunes/binary/ratiosversus.py
rename to georunes/plot/binary/ratiosversus.py
index d507144..65a69ff 100644
--- a/georunes/binary/ratiosversus.py
+++ b/georunes/plot/binary/ratiosversus.py
@@ -1,6 +1,6 @@
 from matplotlib.colors import to_rgba
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.plotting import normalize_marker_size
 
 
diff --git a/georunes/binary/shand.py b/georunes/plot/binary/shand.py
similarity index 92%
rename from georunes/binary/shand.py
rename to georunes/plot/binary/shand.py
index 035b9e4..d94e914 100644
--- a/georunes/binary/shand.py
+++ b/georunes/plot/binary/shand.py
@@ -1,9 +1,12 @@
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.chemistry import molar_ratio
 from georunes.tools.language import format_chemical_formula as _fml, get_translator
 
 
+# Shand, S.J., 1943. Eruptive rocks: their genesis, composition, and classification, with a chapter on meteorites. J.
+# Wiley & sons, Incorporated.
+
 class DiagramShand(DiagramBase, ArrowDrawer, LegendDrawer):
     def __init__(self, datasource, title="A/NK vs A/CNK diagram of Shand (1943)",
                  annotation=None,
diff --git a/georunes/binary/sialkali.py b/georunes/plot/binary/sialkali.py
similarity index 94%
rename from georunes/binary/sialkali.py
rename to georunes/plot/binary/sialkali.py
index 6201b85..3ded3da 100644
--- a/georunes/binary/sialkali.py
+++ b/georunes/plot/binary/sialkali.py
@@ -1,9 +1,14 @@
 import numpy as np
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.language import format_chemical_formula as _fml, get_translator
 
 
+# Cox, K.G., Bell, J.D. and Pankhurst, R.J., 1979. Petrographic aspects of plutonic rocks. In The Interpretation of
+# Igneous Rocks (pp. 283-307). Springer, Dordrecht.
+# Le Maitre, R.W., Streckeisen, A., Zanettin, B., Le Bas, M.J., Bonin, B. and Bateman, P., 2005. Igneous Rocks: A
+# Classification and Glossary of Terms. Igneous Rocks: A Classification and Glossary of Terms, p.252.
+
 class DiagramSiAlkali(DiagramBase, ArrowDrawer, LegendDrawer):
 
     def __init__(self, datasource, decor_set="Cox", title="Total alkali versus silica diagram",
diff --git a/georunes/binary/versus.py b/georunes/plot/binary/versus.py
similarity index 97%
rename from georunes/binary/versus.py
rename to georunes/plot/binary/versus.py
index ad2e4b6..c6708e2 100644
--- a/georunes/binary/versus.py
+++ b/georunes/plot/binary/versus.py
@@ -1,6 +1,6 @@
 from matplotlib.colors import to_rgba
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer, ArrowDrawer
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer, ArrowDrawer
 from georunes.tools.chemistry import val_el_to_mol
 from georunes.tools.plotting import normalize_marker_size
 
diff --git a/georunes/base/helpers.py b/georunes/plot/helpers.py
similarity index 100%
rename from georunes/base/helpers.py
rename to georunes/plot/helpers.py
diff --git a/georunes/base/spider.py b/georunes/plot/spider.py
similarity index 97%
rename from georunes/base/spider.py
rename to georunes/plot/spider.py
index 23c58b8..acab176 100644
--- a/georunes/base/spider.py
+++ b/georunes/plot/spider.py
@@ -1,8 +1,10 @@
 import warnings
+
 import matplotlib.patches as patches
 import matplotlib.pyplot as plt
 from matplotlib import lines
-from georunes.base.base import DiagramBase
+
+from georunes.plot.base import DiagramBase
 from georunes.tools.chemistry import val_ox_to_el_ppm, name_el_to_def_ox
 from georunes.tools.data import row_min, row_max
 from georunes.tools.preprocessing import check_geochem_res
@@ -23,7 +25,7 @@ def __init__(self, datasource,
                  enclosed_in_bg=("",),
                  markersize=8,
                  xlabel="", ylabel="",
-                 drawing_order = "zorder",
+                 drawing_order="zorder",
                  **kwargs
                  ):
 
@@ -49,7 +51,7 @@ def __init__(self, datasource,
             self.norm = res.compos["CI"]
         self.fillmode = fillmode
         self.enclosed_in_bg = enclosed_in_bg
-        self.drawing_order = drawing_order
+        self.drawing_order = drawing_order if drawing_order in self.data.columns else None
         self.xlabel = xlabel
         self.ylabel = ylabel
         self.ylim = ylim
@@ -183,7 +185,8 @@ def plot(self):
                                                  )
 
                 if custom_fillmode in ("enclosed", "enclosed-lines", "mixed"):
-                    self.ax.fill_between(self.listing, min_vals, max_vals, facecolor=list(group["color"])[0], zorder=zorder,
+                    self.ax.fill_between(self.listing, min_vals, max_vals, facecolor=list(group["color"])[0],
+                                         zorder=zorder,
                                          alpha=0.25)
 
                     # For legend
diff --git a/georunes/plot/ternary/__init__.py b/georunes/plot/ternary/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/georunes/ternary/afm.py b/georunes/plot/ternary/afm.py
similarity index 92%
rename from georunes/ternary/afm.py
rename to georunes/plot/ternary/afm.py
index ff6ea7a..e9c1792 100644
--- a/georunes/ternary/afm.py
+++ b/georunes/plot/ternary/afm.py
@@ -1,9 +1,12 @@
-from georunes.base.helpers import ArrowDrawerTernary
-from georunes.base.ternary import DiagramTernaryBase
+from georunes.plot.helpers import ArrowDrawerTernary
+from georunes.plot.ternary.ternbase import DiagramTernaryBase
 from georunes.tools.language import get_translator
 from georunes.tools.plotting import normalize_marker_size
 
 
+# Irvine, T.N. and Baragar, W.R.A., 1971. A Guide to the Chemical Classification of the Common Volcanic Rocks.
+# Canadian Journal of Earth Sciences, 8(5), pp.523-548.
+
 class DiagramAFM(DiagramTernaryBase, ArrowDrawerTernary):
     def __init__(self, datasource, title="AFM diagram (Irvine and Baragar 1971)", **kwargs):
         DiagramTernaryBase.__init__(self, datasource=datasource,
diff --git a/georunes/ternary/anabor.py b/georunes/plot/ternary/anabor.py
similarity index 93%
rename from georunes/ternary/anabor.py
rename to georunes/plot/ternary/anabor.py
index 9559779..3e0b4a8 100644
--- a/georunes/ternary/anabor.py
+++ b/georunes/plot/ternary/anabor.py
@@ -1,9 +1,14 @@
-from georunes.base.helpers import ArrowDrawerTernary
-from georunes.base.ternary import DiagramTernaryBase
+from georunes.plot.helpers import ArrowDrawerTernary
+from georunes.plot.ternary.ternbase import DiagramTernaryBase
 from georunes.tools.language import get_translator
 from georunes.tools.plotting import normalize_marker_size
 
 
+# Barker, F., 1979. Trondhjemite: definition, environment and hypotheses of origin. In Developments in petrology (
+# Vol. 6, pp. 1-12). Elsevier.
+# O'Connor, J.T., 1965. A classification of quartz-rich igneous rocks based on feldspar ratios.
+# US Geol. Surv. Prof. Paper, 82.
+
 class DiagramAnAbOr(DiagramTernaryBase, ArrowDrawerTernary):
     def __init__(self, datasource, title="Ternary feldspar plot (anorthite, albite and orthoclase)",
                  decor_set="Barker79", padding={"bottom": 0.124}, **kwargs):
@@ -22,8 +27,8 @@ def __init__(self, datasource, title="Ternary feldspar plot (anorthite, albite a
 
     def set_decoration(self, ):
         _ = get_translator(self.lang_cfg)
-        # Order for points : Or, An, Ab
 
+        # Order for points : Or, An, Ab
         if self.decor_set == "blank":
             pass
 
diff --git a/georunes/ternary/jensen.py b/georunes/plot/ternary/jensen.py
similarity index 95%
rename from georunes/ternary/jensen.py
rename to georunes/plot/ternary/jensen.py
index 9507b5f..3415967 100644
--- a/georunes/ternary/jensen.py
+++ b/georunes/plot/ternary/jensen.py
@@ -1,11 +1,14 @@
-from georunes.base.helpers import ArrowDrawerTernary
-from georunes.base.ternary import DiagramTernaryBase
+from georunes.plot.helpers import ArrowDrawerTernary
+from georunes.plot.ternary.ternbase import DiagramTernaryBase
 from georunes.tools.chemistry import val_ox_to_el, val_ox_to_mc
 from georunes.tools.colors import darken_color
 from georunes.tools.language import get_translator
 from georunes.tools.plotting import normalize_marker_size
 
 
+# Jensen, L.S., 1976. A new plot for classifying subalkalic volcanic rocks. Ontario Division of Mines,
+# Miscellaneous Paper, 66, pp.1-22.
+
 class DiagramJensen(DiagramTernaryBase, ArrowDrawerTernary):
     def __init__(self, datasource, title="Ternary plot of Jensen (1976)", **kwargs):
         DiagramTernaryBase.__init__(self, datasource=datasource, title=title,
diff --git a/georunes/plot/ternary/nesbitt_young.py b/georunes/plot/ternary/nesbitt_young.py
new file mode 100644
index 0000000..3d31af7
--- /dev/null
+++ b/georunes/plot/ternary/nesbitt_young.py
@@ -0,0 +1,140 @@
+import math
+import numpy as np
+from georunes.plot.helpers import ArrowDrawerTernary
+from georunes.plot.ternary.ternbase import DiagramTernaryBase
+from georunes.tools.chemistry import molar_ratio, molar_ratio_specified
+from georunes.tools.language import format_chemical_formula as _fml
+from georunes.tools.language import get_translator
+from georunes.tools.plotting import normalize_marker_size
+
+sin60 = math.sqrt(3) / 2
+tan30 = 1 / math.sqrt(3)
+tan60 = math.sqrt(3)
+
+
+# Nesbitt, H.W. and Young, G.M., 1984.
+# Prediction of some weathering trends of plutonic and volcanic rocks based on thermodynamic and
+# kinetic considerations. Geochimica et cosmochimica acta, 48(7).
+
+class DiagramNesbittYoung(DiagramTernaryBase, ArrowDrawerTernary):
+    def __init__(self, datasource, title="Nesbitt and Young diagram",
+                 markersize=20, ignore_apatite_correction=False, **kwargs):
+        DiagramTernaryBase.__init__(self, datasource=datasource, title=title,
+                                    markersize=markersize,
+                                    top_label=_fml('A (Al2O3)'), left_label=_fml('CN (CaO+Na2O)'),
+                                    right_label=_fml('K (K2O)'),
+                                    no_ticks_label=True,
+                                    **kwargs)
+
+        self.ignore_apatite_correction = ignore_apatite_correction
+
+    def set_decoration(self, ):
+        _ = get_translator(self.lang_cfg)
+        # Order for points : right var, top var, left var
+
+        _ax = self.tax.get_axes()
+
+        _ax.plot((0, 50 + 35 * sin60 * tan30), (sin60 * 65, sin60 * 65), linewidth=0.8, color=self.decor_line_col,
+                 linestyle="--")
+        _ax.plot((0, 50 + 15 * sin60 * tan30), (sin60 * 85, sin60 * 85), linewidth=0.8, color=self.decor_line_col,
+                 linestyle="--")
+        _ax.plot((0, 50 + 50 * sin60 * tan30), (sin60 * 50, sin60 * 50), linewidth=0.8, color=self.decor_line_col,
+                 linestyle="--")
+
+        self.ax.arrow(35, 50 * sin60, (30 * sin60) / tan60, 30 * sin60, linewidth=1, color="#222222",
+                      head_width=3, head_length=3)
+
+        _ax.annotate(_("Weak\nweathering"), (3, 45,), fontsize="small", color=self.decor_text_col, fontweight='bold')
+        _ax.annotate(_("Intermediate\nweathering"), (3, 60,), fontsize="small", color=self.decor_text_col,
+                     fontweight='bold')
+        _ax.annotate(_("Strong\nweathering"), (3, 75,), fontsize="small", color=self.decor_text_col, fontweight='bold')
+        self.tax.annotate(_("Hornblende"), (3, 16,), color=self.decor_text_col, fontsize="small")
+        self.tax.annotate(_("Pyroxene"), (3, 6,), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("Plagioclase"), (3, sin60 * 46,), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("Smectite"), (25, sin60 * 80,), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("Kaolinite, gibbsite, chlorite"), (55, 85,), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("Illite"), (64, int(80 * sin60),), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("Muscovite"), (68, int(72 * sin60),), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("Biotite"), (74, int(58 * sin60),), color=self.decor_text_col, fontsize="small")
+        _ax.annotate(_("K-feldspar"), (79, int(50 * sin60),), color=self.decor_text_col, fontsize="small")
+
+        minerals = [
+            (0, 8, 92),  # Clinopyroxene
+            (0, 18, 82),  # Hornblende
+            (0, 50, 50),  # Plagioclase
+            (0, 80, 20),  # Smectite
+            (0, 100, 0),  # Kaolininte, gibbsite, Chlorite
+            (20, 80, 0),  # Illite
+            (28, 72, 0),  # Muscovite
+            (42, 58, 0),  # Biotite
+            (50, 50, 0)  # K-feldspar
+        ]
+
+        self.tax.scatter(minerals, marker="o", facecolors="#333333", zorder=4, s=20)
+        self.ax.arrow(0, 0, 0, 50 * tan60, head_width=3, head_length=3, fc='k', ec='k')
+
+        # Draw vertical axis ticks
+        cia_val = np.arange(0, 101, 10)
+        for tick_value in cia_val:
+            y_tick = tick_value * 50 * tan60 / 100
+            self.ax.plot((0, 0.5), (y_tick, y_tick), linewidth=1, color=self.decor_line_col, )
+            _ax.annotate(tick_value, (-3, y_tick - 2,), color=self.decor_text_col, fontsize="small", va="bottom",
+                         ha="right")
+
+        _ax.annotate("CIA", (2, 22,), color=self.decor_text_col, rotation=90, )
+
+    def plot(self):
+        DiagramTernaryBase.plot_config(self, no_gridline=True, no_ticks=self.no_ticks,
+                                       no_ticks_label=self.no_ticks_label)
+        self.set_decoration()
+
+        # Categorize by group
+        groups = self.data.groupby(self.group_name)
+        for name, group in groups:
+
+            if self.is_group_allowed(name):
+                top_var = molar_ratio(group['Al2O3'])
+                if self.ignore_apatite_correction:
+                    left_var = molar_ratio_specified(group['CaO'], "CaO") + molar_ratio(group['Na2O'])
+                else:
+                    left_var = molar_ratio_specified(group['CaO'] - 55.07 * group["P2O5"] / 41.85, "CaO") + molar_ratio(
+                        group['Na2O'])
+
+                right_var = molar_ratio(group['K2O'])
+
+                var_sum = top_var + left_var + right_var
+                norm_z = 100 * left_var / var_sum
+                norm_x = 100 * right_var / var_sum
+                norm_y = 100 * top_var / var_sum
+
+                if self.marker != '':
+                    mrk = self.marker
+                else:
+                    mrk = list(group["marker"])[0]
+
+                if self.marker_size_scaled():
+                    size = normalize_marker_size(group[self.markersize['var_scale']], self.markersize['val_max'],
+                                                 self.markersize['val_min'], self.markersize['size_max'],
+                                                 self.markersize['size_min'])
+                else:
+                    size = self.markersize
+
+                # Convert points to tuples
+                points = [(norm_x.get(i), norm_y.get(i), norm_z.get(i)) for i, sample in top_var.iteritems()]
+
+                zorder = 4
+                if self.drawing_order:
+                    zorder = list(group[self.drawing_order])[0]
+
+                self.tax.scatter(points, edgecolors=group["color"],
+                                 marker=mrk, label=name, facecolors=group["color"], s=size,
+                                 alpha=self.alpha_color, zorder=zorder)
+
+                if self.annotation:
+                    for i, sample in group[self.top_var].iteritems():
+                        self.tax.annotate(group[self.annotation].get(i), (norm_x.get(i), norm_y.get(i)),
+                                          fontsize='xx-small')
+
+        self.plot_arrows()
+        self.fig.subplots_adjust(right=0.92)
+        self.plot_legend()
diff --git a/georunes/ternary/qap.py b/georunes/plot/ternary/qap.py
similarity index 92%
rename from georunes/ternary/qap.py
rename to georunes/plot/ternary/qap.py
index 069b119..961d3db 100644
--- a/georunes/ternary/qap.py
+++ b/georunes/plot/ternary/qap.py
@@ -1,9 +1,15 @@
-from georunes.base.helpers import ArrowDrawerTernary
-from georunes.base.ternary import DiagramTernaryBase
+from georunes.plot.helpers import ArrowDrawerTernary
+from georunes.plot.ternary.ternbase import DiagramTernaryBase
 from georunes.tools.language import get_translator
 from georunes.tools.plotting import normalize_marker_size
 
 
+# Streckeisen, A., 1974. Classification and nomenclature of plutonic rocks recommendations of the IUGS subcommission
+# on the systematics of igneous rocks. Geologische Rundschau, 63(2), pp.773-786.
+# Streckeisen, A., 1978. IUGS Subcommission on the Systematics of Igneous Rocks. Classification and Nomenclature of
+# Volcanic Rocks, Lamprophyres, Carbonatites and Melilite Rocks. Recommendations and Suggestions. Neues Jahrbuch fur
+# Mineralogie. Stuttgart. Abhandlungen, 143, pp.1-14.
+
 class DiagramQAP(DiagramTernaryBase, ArrowDrawerTernary):
     def __init__(self, datasource, title="QAP ternary diagram", decor_set="Plut", decor_text_col="#777777", **kwargs):
         if decor_set not in ("Plut", "Volc"):
diff --git a/georunes/base/ternary.py b/georunes/plot/ternary/ternbase.py
similarity index 92%
rename from georunes/base/ternary.py
rename to georunes/plot/ternary/ternbase.py
index 7c0bb68..f5065e5 100644
--- a/georunes/base/ternary.py
+++ b/georunes/plot/ternary/ternbase.py
@@ -1,8 +1,10 @@
 import math
+
 import matplotlib.pyplot as plt
 import ternary
-from georunes.base.base import DiagramBase
-from georunes.base.helpers import LegendDrawer
+
+from georunes.plot.base import DiagramBase
+from georunes.plot.helpers import LegendDrawer
 from georunes.tools.plotting import normalize_marker_size
 
 
@@ -64,7 +66,10 @@ def plot_config(self, no_gridline=False, no_ticks=False, no_ticks_label=False):
         lw = 0 if no_ticks else 1
 
         if not no_ticks:
-            self.tax.ticks(multiple=10, axis='lbr', linewidth=lw, fontsize=10, linestyle='-', offset=0.02)
+            if no_ticks_label:
+                self.tax.ticks(multiple=10, axis='lbr', linewidth=lw, linestyle='-', tick_formats="")
+            else:
+                self.tax.ticks(multiple=10, axis='lbr', linewidth=lw, linestyle='-', fontsize=10, offset=0.02)
 
         if self.vertical_ticks:
             self.ax.set_xlim()
diff --git a/georunes/tools/chemistry.py b/georunes/tools/chemistry.py
index 05c8186..fdb655a 100644
--- a/georunes/tools/chemistry.py
+++ b/georunes/tools/chemistry.py
@@ -50,16 +50,22 @@
     'MnO': 0.774,
     'MgO': 0.603,
     'CaO': 0.715,
-    'Li2O': 0.465,
+    'Li2O': 0.464570,
     'Na2O': 0.742,
     'K2O': 0.830,
-    'Rb2O': 0.914,
-    'Cs2O': 0.943,
-    'BaO': 0.896,
-    'Cr2O3': 0.684,
+    'Rb2O': 0.914412,
+    'Cs2O': 0.943226,
+    'BaO': 0.895651,
+    'Cr2O3': 0.684202,
     'P2O5': 0.436,
     'SnO': 0.881,
-    'Ta2O5': 0.819
+    'Ta2O5': 0.819,
+    'NiO': 0.785797,
+    'CoO': 0.786483,
+    'SrO': 0.845595,
+    'ZrO2': 0.740318,
+    'V2O3': 0.684324,
+    'SO3': 0.400504,
 }
 
 nb_cation = {
@@ -86,29 +92,40 @@
 }
 
 molar_mass = {
-    'SiO2': 60.0843,
-    'TiO2': 79.8988,
-    'Al2O3': 101.9613,
-    'Fe2O3': 159.692,
-    'Fe2O3t': 159.692,
-    'FeO': 71.8464,
-    'FeOt': 71.8464,
-    'MnO': 70.9374,
-    'MgO': 40.3114,
-    'CaO': 56.0794,
-    'Li2O': 29.882,
-    'Na2O': 61.9788,
-    'K2O': 94.2034,
-    'Rb2O': 186.9356,
-    'Cs2O': 281.8109,
-    'BaO': 153.3391,
-    'Cr2O3': 151.989,
-    'P2O5': 141.943,
+    'SiO2': 60.0843,  # a
+    'TiO2': 79.8658,  # a
+    'Al2O3': 101.961276,  # a
+    'Fe2O3': 159.6882,  # a
+    'Fe2O3t': 159.6882,  # a
+    'FeO': 71.8444,  # a
+    'FeOt': 71.8444,  # a
+    'MnO': 70.937449,  # a
+    'MgO': 40.3044,  # a
+    'CaO': 56.0774,  # a
+    'Li2O': 29.8814,  # a
+    'Na2O': 61.97894,  # a
+    'K2O': 94.1960,  # a
+    'SrO': 103.6194,  # a
+    'BaO': 153.3264,  # a
+    'Rb2O': 186.935,  # a
+    'Cs2O': 281.8103,  # a
+    'ZrO2': 123.228,  # a
+    'P2O5': 141.944522,  # a
+    'NiO': 74.6928,  # a
+    'CoO': 74.9326,  # a
+    'Cr2O3': 151.9904,  # a
+    'V2O3': 149.8812,  # a
+    'CO2': 44.0095,  # a
+    'SO3': 80.0642,  # a
     'SnO': 134.689,
     'Ta2O5': 441.891
 }
 
 el_molar_mass = {
+    'F': 18.9984,  # a
+    'O': 15.9994,  # a
+    'Cl': 35.4527,  # a
+    'S': 32.066,  # a
     'Nb': 92.90638,
     'Ta': 180.9479,
     'Nb_Ta': 92.90638 / 180.9479,
@@ -148,3 +165,6 @@ def molar_ratio(group):
 
 def molar_ratio_specified(group, group_name):
     return group / molar_mass[group_name]
+
+# a. Verma, S.P., Torres-Alvarado, I.S. and Velasco-Tapia, F., 2003.
+# A revised CIPW norm. Swiss Bulletin of Mineralogy and Petrology, 83(2), pp.197-216.
diff --git a/georunes/tools/filemanager.py b/georunes/tools/filemanager.py
index e5303d3..aee2cea 100644
--- a/georunes/tools/filemanager.py
+++ b/georunes/tools/filemanager.py
@@ -39,7 +39,7 @@ def read_file(self, datasource, sheet_name=None, sep=",", delimiter=None, index=
             if fext in (".xls", ".xlsx"):
                 data = pd.read_excel(datasource, sheet_name=sheet_name, index_col=index)
             elif fext in (".csv", ".txt"):
-                data = pd.read_csv(datasource, sep, delimiter, index_col=index)
+                data = pd.read_csv(datasource, sep=sep, delimiter=delimiter, index_col=index)
             else:
                 msg = "Extension file " + fext + " not recognized."
                 raise Exception(msg)
diff --git a/georunes/tools/formatters.py b/georunes/tools/formatters.py
index 52f06a7..d529fd3 100644
--- a/georunes/tools/formatters.py
+++ b/georunes/tools/formatters.py
@@ -9,7 +9,6 @@ def __init__(self, max_num=1000000, min_num=0.001,
         self.min_num = min_num
         LogFormatter.__init__(self, **kwargs)
 
-
     def _num_to_string(self, x, vmin, vmax):
         if x > self.max_num:
             s = '%1.0e' % x
diff --git a/georunes/tools/plotting.py b/georunes/tools/plotting.py
index e4973c6..c25e1f6 100644
--- a/georunes/tools/plotting.py
+++ b/georunes/tools/plotting.py
@@ -1,4 +1,5 @@
 import math
+
 import numpy as np
 from scipy import interpolate
 
@@ -16,13 +17,13 @@ def get_spline(x, y, prec=50):
     # resample it at smaller distance intervals
     interp_d = np.linspace(dist_along[0], dist_along[-1], prec)
     interp_x, interp_y = interpolate.splev(interp_d, spline)
-    return (interp_x, interp_y)
+    return interp_x, interp_y
 
 
 def tern_coords_to_bin_coords(rvar, tvar, lvar=None, scale=100):
     xvar = tvar / 2 + rvar
     yvar = tvar * sin60
-    return (xvar, yvar)
+    return xvar, yvar
 
 
 def normalize_marker_size(serie, val_max, val_min, size_max, size_min):
@@ -41,4 +42,4 @@ def normalize_marker_size(serie, val_max, val_min, size_max, size_min):
 
 def is_in_canvas(x, y, xlim, ylim):
     # Check if element is in canvas
-    return ((x > xlim[0]) and (x < xlim[1]) and (y > ylim[0]) and (y < ylim[1]))
+    return (x > xlim[0]) and (x < xlim[1]) and (y > ylim[0]) and (y < ylim[1])
diff --git a/georunes/tools/preprocessing.py b/georunes/tools/preprocessing.py
index 3b30b99..fbf48a1 100644
--- a/georunes/tools/preprocessing.py
+++ b/georunes/tools/preprocessing.py
@@ -1,9 +1,11 @@
 import itertools
 import warnings
 from pathlib import Path
+
 import matplotlib.cm as cm
 import pandas as pd
 from matplotlib.colors import to_hex
+
 from georunes.tools.data import unique
 
 _vtc = ('color', 'marker', 'order', 'zorder', 'label',)
@@ -59,7 +61,7 @@ def populate_file(file, sheet=0, group_name='group', cmap='viridis',
     if fext in (".xls", ".xlsx"):
         data = pd.read_excel(file, sheet_name=sheet)
     elif fext in (".csv", ".txt"):
-        data = pd.read_csv(file, sep, delimiter)
+        data = pd.read_csv(file, sep=sep, delimiter=delimiter)
     if group_name not in data.columns:
         raise ValueError("The specified 'group_name' parameter is not in the excel sheet.")
 
@@ -131,14 +133,15 @@ def populate_file_with_config(file, config, sheet=0, group_name='group', output_
     card = len(config['group'])
     for key in attribs:
         if key not in ('group', 'label') and len(config[key]) != card:
-            msg = "The length of the parameter '" + key + "' is different from the number of group defined in the config."
+            msg = "The length of the parameter '" + key + \
+                  "' is different from the number of group defined in the config."
             raise ValueError(msg)
 
     # Load data
     if fext in (".xls", ".xlsx"):
         data = pd.read_excel(file, sheet_name=sheet)
     elif fext in (".csv", ".txt"):
-        data = pd.read_csv(file, sep, delimiter)
+        data = pd.read_csv(file, sep=sep, delimiter=delimiter)
     if group_name not in data.columns:
         raise ValueError("The specified 'group_name' parameter is not in the excel sheet.")
 


Test patch:
diff --git a/examples/modal mineralogy/cipw_test.csv b/examples/modal mineralogy/cipw_test.csv
new file mode 100644
index 0000000..d4aaaf9
--- /dev/null
+++ b/examples/modal mineralogy/cipw_test.csv
@@ -0,0 +1,8 @@
+Sample,SiO2,TiO2,Al2O3,Fe2O3,FeO,MnO,MgO,CaO,Na2O,K2O,P2O5,CO2,SO3,LOI,Cr,F,Cl,Total
+Granite-1,70.06,0.47,15.1,2.7,0,0.03,0.66,2.08,3.66,4.2,0.12,,,,,8200.00,,99.90
+Basalte-1,50.33,1.3,13.67,14.66,0,0.22,5.85,8.1,2.58,0.31,0.08,,,,,,,97.10
+Basalte-2,50.87,0.97,12.99,9.2,0,0.19,6.45,10.5,2.25,1.65,0.3,,0.03,,,,,95.40
+Peridotite,40.33,0.18,3.99,13.86,0,0.18,30.57,2.36,0.25,0.06,0.01,,,7.82,3550,,,99.97
+Syenite-1,57.99,0.86,15.96,3.14,2.11,0.1,1.17,4.39,5.7,6.88,0.13,0.2,,,,,,98.63
+Syenite-2,58.2,0.49,18.54,4.51,0,0.07,0.59,2.13,3.06,11.43,0.11,,,0.34,,,,99.47
+Syenite-3,56.27,0.22,21.31,2.79,0,0.02,0.14,0.39,0.32,17.43,0.03,,,0.69,,,,99.61
\ No newline at end of file
diff --git a/examples/modal mineralogy/cipw_test.py b/examples/modal mineralogy/cipw_test.py
new file mode 100644
index 0000000..32a1a30
--- /dev/null
+++ b/examples/modal mineralogy/cipw_test.py
@@ -0,0 +1,11 @@
+import os
+from georunes.modmin.norm.cipw import CIPWNorm
+from georunes.tools.filemanager import FileManager
+
+filemanager = FileManager.get_instance()
+directory = os.path.dirname(os.path.realpath("cipw_test.csv"))
+source = os.path.join(directory, "cipw_test.csv")
+data = filemanager.read_file(source)
+
+cipw = CIPWNorm(verbose=1)
+p, s = cipw.compute(data, skip_cols=1, minor_included=True, co2_calcite=0.8, co2_cancrinite=0.2)
diff --git a/examples/modal mineralogy/gradient_descent_test.py b/examples/modal mineralogy/gradient_descent_test.py
new file mode 100644
index 0000000..2ac8981
--- /dev/null
+++ b/examples/modal mineralogy/gradient_descent_test.py
@@ -0,0 +1,27 @@
+import os
+from georunes.tools.filemanager import FileManager
+from georunes.modmin.optim.gd import GradientDescent
+
+filemanager = FileManager.get_instance()
+directory = os.path.dirname(os.path.realpath("modalmin_test.csv"))
+source_comp = os.path.join(directory, "modalmin_test.csv")
+source_minerals = os.path.join(directory, "minerals.csv")
+data = filemanager.read_file(source_comp)
+minerals = filemanager.read_file(source_minerals)
+
+gd = GradientDescent(verbose=2, dist_func="SMAPE")  # SMAPE = Symmetric mean absolute percentage error
+
+# A starting composition used for the gradient descent method
+starting_partition = {
+    "Albite": 25,
+    "Anorthite": 10,
+    "Quartz": 28.999,
+    "Hematite": 2,
+    "Biotite": 6,
+    "Muscovite": 8,
+    "Orthoclase": 20,
+    'Apatite-F': 0.001,
+}
+
+p, s = gd.compute(data, skip_cols=1, raw_minerals_data=minerals,
+                  max_iter=100000, learn_rate=0.0001, starting_partition=starting_partition, force_totals=True)
diff --git a/examples/modal mineralogy/ls_tests.py b/examples/modal mineralogy/ls_tests.py
new file mode 100644
index 0000000..3a5c536
--- /dev/null
+++ b/examples/modal mineralogy/ls_tests.py
@@ -0,0 +1,18 @@
+import os
+from georunes.tools.filemanager import FileManager
+from georunes.modmin.optim.bvls import BVLS
+from georunes.modmin.optim.nnls import NNLS
+
+filemanager = FileManager.get_instance()
+directory = os.path.dirname(os.path.realpath("modalmin_test.csv"))
+source_comp = os.path.join(directory, "modalmin_test.csv")
+source_minerals = os.path.join(directory, "minerals.csv")
+data = filemanager.read_file(source_comp)
+minerals = filemanager.read_file(source_minerals)
+
+bvls = BVLS(verbose=1)  # Bounded-variable least squares
+nnls = NNLS(verbose=1)  # Non negative least squares
+
+p1, s1 = bvls.compute(data, skip_cols=1, raw_minerals_data=minerals)
+print("\n\n")
+p2, s2 = nnls.compute(data, skip_cols=1, raw_minerals_data=minerals)
diff --git a/examples/modal mineralogy/modalmin_test.csv b/examples/modal mineralogy/modalmin_test.csv
new file mode 100644
index 0000000..970fa9a
--- /dev/null
+++ b/examples/modal mineralogy/modalmin_test.csv
@@ -0,0 +1,3 @@
+Sample,SiO2,TiO2,Al2O3,Fe2O3,MnO,MgO,CaO,Na2O,K2O,P2O5,Total
+Granite-1,70.06,0.47,15.1,2.7,0.03,0.66,2.08,3.66,4.2,0.12,99.08
+Syenite-1,58.2,0.49,18.54,4.51,0.07,0.59,2.13,3.06,11.43,0.11,99.13
\ No newline at end of file
diff --git a/examples/modal mineralogy/random_search_test.py b/examples/modal mineralogy/random_search_test.py
new file mode 100644
index 0000000..e94e5ab
--- /dev/null
+++ b/examples/modal mineralogy/random_search_test.py
@@ -0,0 +1,19 @@
+import os
+from georunes.tools.filemanager import FileManager
+from georunes.modmin.optim.randsearch import RandomSearch
+
+filemanager = FileManager.get_instance()
+directory = os.path.dirname(os.path.realpath("modalmin_test.csv"))
+source_comp = os.path.join(directory, "modalmin_test.csv")
+source_minerals = os.path.join(directory, "minerals.csv")
+data = filemanager.read_file(source_comp)
+minerals = filemanager.read_file(source_minerals)
+
+rs = RandomSearch(verbose=1, dist_func='MAE')  # MAE = Mean absolute error
+ratios = {"Plagioclase": [["Albite", "Anorthite"],
+                          [76, 24]], }  # Define plagioclase as weigh proportions of 76% albite and 24% anorthite
+
+# Optimize solution by checking random compositions in a defined hypercube. When no new random solution can be found
+# (because of bound limits), the hypercube is reduced by multiplying the semi-edge with the parameter 'scale_semiedge'
+p, s = rs.compute(data, skip_cols=1, raw_minerals_data=minerals, ratios=ratios,
+                  max_iter=100000, search_semiedge=0.1, scale_semiedge=0.5, force_totals=True)