YOUR_MATERIAL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages\n",
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import sys\n",
    "from multiprocessing import Process, Queue\n",
    "import pandas as pd\n",
    "import optuna\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "\n",
    "sys.path.append('/opt/conda/GSASII/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configurations\n",
    "\n",
    "### Change here ###\n",
    "STUDY_NAME = 'YOUR_MATERIAL'\n",
    "RANDOM_SEED = 1024\n",
    "\n",
    "DATA_DIR = '/bbo_rietveld/data/' + STUDY_NAME\n",
    "# all output files include GSAS project file (*.gpx) will be saved in WORK_DIR\n",
    "WORK_DIR = '/bbo_rietveld/work/' + STUDY_NAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make directories\n",
    "! rm -f $WORK_DIR/$STUDY_NAME*\n",
    "! mkdir -p $WORK_DIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ProjectBBO:\n",
    "    def __init__(self, trial_number):\n",
    "        import GSASIIscriptable as G2sc\n",
    "        import shutil\n",
    "        \n",
    "        # Create a project with a default project name\n",
    "        ### Change here ###\n",
    "        shutil.copyfile(DATA_DIR+'/'+'YOUR_PROJECT_FILE.gpx', \n",
    "                        WORK_DIR+'/'+'BBO_seed{0}_trial_{1}.gpx'.format(RANDOM_SEED, trial_number))\n",
    "        \n",
    "        self.gpx = G2sc.G2Project(gpxfile=os.path.join(WORK_DIR, 'BBO_seed{0}_trial_{1}.gpx'.format(RANDOM_SEED, trial_number)))\n",
    "\n",
    "        # Add two histograms to the project\n",
    "        self.hist1 = self.gpx.histograms()[0]\n",
    "        self.phase0 = self.gpx.phases()[0]\n",
    "        self.hist1.data['Instrument Parameters'][0]['I(L2)/I(L1)'] = [0.5, 0.5, 0]\n",
    "\n",
    "        # Set to use iso\n",
    "        for val in self.phase0.data['Atoms']:\n",
    "            val[9] = 'I'\n",
    "\n",
    "    def refine_and_calc_Rwp(self, param_dict):\n",
    "        self.gpx.do_refinements([param_dict])\n",
    "        for hist in self.gpx.histograms():\n",
    "            _, Rwp = hist.name, hist.get_wR()\n",
    "        return Rwp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(trial):\n",
    "    \"\"\"\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    trial : optuna.trial object\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Rwp : float\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # Here, you should define search space and perform the refinement.\n",
    "    # Please see other notebook.\n",
    "    # Copy and paste from other notebook and some modifications would be enough.\n",
    "    \n",
    "    return Rwp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Optuna study\n",
    "study = optuna.create_study(study_name=STUDY_NAME + '_seed%s' % (RANDOM_SEED),\n",
    "                            sampler=optuna.samplers.TPESampler(n_startup_trials=20, seed=RANDOM_SEED))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run 200 refinements to find the best configuration. It may take abount an hour to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimize\n",
    "study.optimize(objective, n_trials=200, n_jobs=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Results\n",
    "df = study.trials_dataframe()\n",
    "df.columns = [' '.join(col).replace('params', '').strip() for col in df.columns.values]\n",
    "df.rename(columns={'value':'Rwp', 'number':'trial'}, inplace=True)\n",
    "df.drop(columns=['state', 'system_attrs _number'], inplace=True)\n",
    "df.sort_values('Rwp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Best configuration\n",
    "study.best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Best Rwp\n",
    "study.best_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rwp plot\n",
    "def rwp_plot():\n",
    "    minvalues = [df.iloc[0]['Rwp']]\n",
    "    for i in range(1, df.shape[0]):\n",
    "        minvalues.append(min(minvalues[-1], df.iloc[i]['Rwp']))\n",
    "    minvalues = pd.DataFrame(minvalues)\n",
    "    \n",
    "    minvalues.plot(legend=None)\n",
    "#     plt.ylim([6, 16])\n",
    "    plt.grid(color='#cccccc')\n",
    "    plt.ylabel('$R_{wp}$')\n",
    "    plt.xlabel('Number of trials')\n",
    "    plt.show()\n",
    "    \n",
    "rwp_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rietveld plot\n",
    "def rietveld_plot():\n",
    "    import GSASIIscriptable as G2sc\n",
    "\n",
    "    gpx = G2sc.G2Project(\n",
    "        '%s/%s_seed%s_trial_%s.gpx' % (WORK_DIR, STUDY_NAME, RANDOM_SEED, study.best_trial.number))\n",
    "\n",
    "    hist1 = gpx.histograms()[0]\n",
    "    phase0 = gpx.phases()[0]\n",
    "\n",
    "    hist = hist1\n",
    "    i = 5\n",
    "    two_theta = hist.getdata(\"X\")[::i]\n",
    "    Yobs = hist.getdata(\"Yobs\")[::i]\n",
    "    Ycalc = hist.getdata(\"Ycalc\")[::i]\n",
    "    bg = hist.getdata(\"Background\")[::i]\n",
    "    residual = hist.getdata(\"Residual\")[::i]\n",
    "\n",
    "    fig = plt.figure()\n",
    "    gs = GridSpec(5, 1, figure=fig)\n",
    "    ax1 = fig.add_subplot(gs[:4, :])\n",
    "    ax2 = fig.add_subplot(gs[4, :])\n",
    "    fig.subplots_adjust(hspace=0)\n",
    "    ax1.grid(color='#cccccc')\n",
    "\n",
    "    ax1.scatter(two_theta, Yobs, marker='P', lw=0.0001, c='Black', label='XRD (Obs)')\n",
    "    ax1.plot(two_theta, Ycalc, label='XRD (Calc)')\n",
    "    ax1.plot(two_theta, bg, color='red', label='Background (Calc)')\n",
    "    ax1.set_ylabel('Intensity')\n",
    "    ax1.legend()\n",
    "    ax2.plot(two_theta, residual, color='blue')\n",
    "    plt.setp(ax1.get_xticklabels(), visible=False);\n",
    "    # ax2.set_ylim(-6600, 6600)\n",
    "    plt.xlabel(r'$2\\theta$ (deg.)')\n",
    "    ax2.set_ylabel('Residual')\n",
    "    # change 2theta range according to your data\n",
    "    ax1.set_xlim(15, 150)\n",
    "    ax2.set_xlim(15, 150)\n",
    "    plt.show()\n",
    "    \n",
    "rietveld_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}